三角形等值线算法，生成三角形网格，后画等值线

#include<stdio.h> #include<graphics.h> #include<math.h> #include<malloc.h> typedef struct tagVERTEX { double x; double y; double z; int v0,v1,b0,b1; /* v0,v1 标志高程点左右的三角型 b0,b1为点在那个边上 */ }VERTEX; typedef struct tagTRIANGLE { int vv0; int vv1; int vv2; }TRIANGLE; #define MAX_VERTEX 500 #define MAX_TRIANGLE 500 void GetData(char fileName[],VERTEX Vertex[],int *pointNumber) { int i=1,tmp; FILE *fp; double x,y,z; if((fp=fopen(fileName,"r"))==NULL) { printf("Open file Error!\n"); exit(0); } while(!feof(fp)) { fscanf(fp,"%d,,%lf,%lf,%lf",&tmp,&x,&y,&z); Vertex[i].x=x;Vertex[i].y=y;Vertex[i].z=z; i++; } fclose(fp); *pointNumber=i-2; } int FvsTriangleInCircle(double xp, double yp, double x1, double y1, double x2, double y2, double x3, double y3, double *xc,double *yc, double *r) { int Ret; double eps; double m1; double m2; double mx1; double mx2; double my1; double my2; double dx; double dy; double rsqr; double drsqr; eps = 0.000001; Ret = 0; if ( fabs(y1 - y2) < eps && fabs(y2 - y3) < eps ) { printf("INCIRCUM - F - Points are coincident !\n"); return Ret; } if ( fabs(y2 - y1) < eps ) { m2 = -(x3 - x2) / (y3 - y2); mx2 = (x2 + x3) / 2; my2 = (y2 + y3) / 2; *xc = (x2 + x1) / 2; *yc = m2 * (*xc - mx2) + my2; } else if ( fabs(y3 - y2) < eps ) { m1 = -(x2 - x1) / (y2 - y1); mx1 = (x1 + x2) / 2; my1 = (y1 + y2) / 2; *xc = (x3 + x2) / 2; *yc = m1 * ((*xc) - mx1) + my1; } else { m1 = -(x2 - x1) / (y2 - y1); m2 = -(x3 - x2) / (y3 - y2); mx1 = (x1 + x2) / 2; mx2 = (x2 + x3) / 2; my1 = (y1 + y2) / 2; my2 = (y2 + y3) / 2; *xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2); *yc = m1 * ((*xc) - mx1) + my1; } dx = x2 - (*xc); dy = y2 - (*yc); rsqr = dx * dx + dy * dy; *r = sqrt(rsqr); dx = xp - (*xc); dy = yp - (*yc); drsqr = dx * dx + dy * dy; if ( drsqr <= rsqr ) return 1; else return 0; } int FvsTriangleWhichSide(double xp, double yp, double x1, double y1, double x2, double y2) { double equation; equation = ((yp - y1) * (x2 - x1)) - ((y2 - y1) * (xp - x1)); if ( equation > 0 ) return -1; else if ( equation == 0 ) return 0; else return 1; } int FvsTrianglate(int nvert,VERTEX Vertex[],TRIANGLE Triangle[]) /* 构成Delaunay */ { int Complete[MAX_VERTEX]; long Edges[3][MAX_VERTEX]; long Nedge; long xmin; long xmax; long ymin; long ymax; long xmid; long ymid; double dx; double dy; double dmax; int i; int j; int k; int ntri; double xc; double yc; double r; int inc; for(i=1; i <= 200; i++) Complete[i] = 0; xmin = Vertex[1].x; ymin = Vertex[1].y; xmax = xmin; ymax = ymin; for ( i = 2; i <= nvert; i++) /* 求最大坐标与最小坐标 */ { if (Vertex[i].x < xmin) xmin = Vertex[i].x; if (Vertex[i].x > xmax) xmax = Vertex[i].x; if (Vertex[i].y < ymin) ymin = Vertex[i].y; if (Vertex[i].y > ymax) ymax = Vertex[i].y; } dx = xmax - xmin; dy = ymax - ymin; if ( dx > dy ) dmax = dx; else dmax = dy; xmid = (xmax + xmin) / 2; ymid = (ymax + ymin) / 2; Vertex[nvert + 1].x = xmid - 2 * (long)dmax; Vertex[nvert + 1].y = ymid - (long)dmax; Vertex[nvert + 2].x = xmid; Vertex[nvert + 2].y = ymid + 2 * (long)dmax; Vertex[nvert + 3].x = xmid + 2 * (long)dmax; Vertex[nvert + 3].y = ymid - (long)dmax; Triangle[1].vv0 = nvert + 1; Triangle[1].vv1 = nvert + 2; Triangle[1].vv2 = nvert + 3; Complete[1] = 0; ntri = 1; for ( i = 1; i <= nvert; i++) { Nedge = 0; j = 0; do { j = j + 1; if ( Complete[j] != 1 ) { inc = FvsTriangleInCircle(Vertex[i].x, Vertex[i].y, Vertex[Triangle[j].vv0].x,Vertex[Triangle[j].vv0].y, Vertex[Triangle[j].vv1].x, Vertex[Triangle[j].vv1].y, Vertex[Triangle[j].vv2].x, Vertex[Triangle[j].vv2].y,&xc, &yc, &r); if ( inc==1 ) { Edges[1][Nedge + 1] = Triangle[j].vv0; Edges[2][Nedge + 1] = Triangle[j].vv1; Edges[1][Nedge + 2] = Triangle[j].vv1; Edges[2][Nedge + 2] = Triangle[j].vv2; Edges[1][Nedge + 3] = Triangle[j].vv2; Edges[2][Nedge + 3] = Triangle[j].vv0; Nedge = Nedge + 3; Triangle[j].vv0 = Triangle[ntri].vv0; Triangle[j].vv1 = Triangle[ntri].vv1; Triangle[j].vv2 = Triangle[ntri].vv2; Complete[j] = Complete[ntri]; j = j - 1; ntri = ntri - 1; } } }while(j < ntri); for ( j = 1; j <= Nedge - 1; j++) { if ( !(Edges[1][j] == 0) && !(Edges[2][j] == 0) ) { for ( k = j + 1; k <= Nedge; k++) { if ( !(Edges[1][k] == 0 ) && !(Edges[2][k] == 0) ) { if ( Edges[1][j] == Edges[2][k] ) { if ( Edges[2][j] == Edges[1][k] ) { Edges[1][j] = 0; Edges[2][j] = 0; Edges[1][k] = 0; Edges[2][k] = 0; } } } } } } for ( j = 1; j <= Nedge; j++) { if ( !(Edges[1][j] == 0) && !(Edges[2][j] == 0 ) ) { ntri = ntri + 1; Triangle[ntri].vv0 = Edges[1][j]; Triangle[ntri].vv1 = Edges[2][j]; Triangle[ntri].vv2 = i; Complete[ntri] =0; } } } /* end of For */ i = 0; do { i = i + 1; if ( Triangle[i].vv0 > nvert || Triangle[i].vv1 > nvert || Triangle[i].vv2 > nvert ) { Triangle[i].vv0 = Triangle[ntri].vv0; Triangle[i].vv1 = Triangle[ntri].vv1; Triangle[i].vv2 = Triangle[ntri].vv2; i = i - 1; ntri = ntri - 1; } }while( i < ntri ); return ntri; } void GetSignPoint(VERTEX p1, VERTEX p2,int *nPoint,int n,int b0,int b1,VERTEX SignPoint[]) /*内插*/ { double eps = 0.000001,dz,dx,dy; int zp1,t=1,i; VERTEX max,min; if(p1.z-p2.z>eps) /* 比较两点的高程值 */ { max=p1; min=p2; } else { max=p2; min=p1; } dz=max.z-min.z; zp1=min.z; if(min.z-zp1<eps) /*处理端点为整数的点*/ { SignPoint[*nPoint]=min; SignPoint[*nPoint].v0=n; SignPoint[*nPoint].v1=0; SignPoint[*nPoint].b0=b0; SignPoint[*nPoint].b1=b1; (*nPoint)++; } while(zp1+t-max.z<eps) /* 内插高程整数据点 */ { SignPoint[*nPoint].x=min.x+(max.x-min.x)*(zp1+t-min.z)/dz; SignPoint[*nPoint].y=min.y+(max.y-min.y)*(zp1+t-min.z)/dz; SignPoint[*nPoint].z=zp1+t; SignPoint[*nPoint].v0=n; SignPoint[*nPoint].v1=0; SignPoint[*nPoint].b0=b0; SignPoint[*nPoint].b1=b1; (*nPoint)++; t++; } } void signTriangle(int b0,int b1,VERTEX D[],int n,int v) /* 标记点所归的三角形 */ { int i; for(i=0;i<n;i++) { if(((D[i].b0==b0)&&(D[i].b1==b1))||((D[i].b1==b0)&&(D[i].b0==b1))) D[i].v1=v; } } int IsMarker(int a[][2],int *n,int vv0,int vv1,int v,VERTEX D[],int np) /* 判断该边 是否已经 插值 等高线点*/ { int i; if(*n==0) { a[*n][0]=vv0; a[*n][1]=vv1;