triangulate.cpp

由TIN生成DEM

/*
TRIANGULATE

Triangulate takes a file containing a random distribution of x y z
points and creates a facet gridded representation of the surface
at a user specified resolution. The grid is generated over the complex
hull of the spot heights.

This is an implementation of the original program that run on the Macintosh
and supports a number of export formats. This version acts like a generator
for the Radiance program.

Call as: 
   triangulate mat name xyzfile xcells ycells [zscale]
where
   xyzfile  the file of spot heights arranged as x y z values
   xcells   the resolution of the gridded mesh in the x direction
   ycells   the resolution of the gridded mesh in the y direction
   zscale   height scaling
The Radiance description of the surface is written to standard output.

*/
#include "stdafx.h"

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "triangulate.h"


//unsigned int  trimax =20000;

PLANE   *demGrid;
int nPlane;

void GenDEM(char filename[255])
{
   unsigned int i,j,k,nv,found[4],corner,duplicate;
   double x,y,z,xmin,xmax,ymin,ymax,zmin,zmax,dx,dy;
   unsigned int ntri; 
   //double xp,yp,xp1,yp1;
   XYZ ptemp,p[5],p1,p2,p3;
	
   XYZ *pxyz;
   TRIANGLE *v;   

   FILE *fptr;  
   fptr = fopen(filename,"r");
   nv = 0;

   char str[255];
    
   while(fscanf(fptr,"%s %lf %lf %lf\n",str,&x,&y,&z) == 4)
      nv++;

   fclose(fptr);
	
   /* Must be 3 or more vertices */
   if (nv < 3) 
   {
      //printf("TRIANGULATE - Insufficient data points\n");
      return;
   }
	
   /* Allocate any necessary memory */
   if ((pxyz = (XYZ *)malloc((nv+4)*(long)sizeof(XYZ))) == NULL) 
   {
     //printf("TRIANGULATE - Unable to allocate enough memory\n");
     return;
   }
   /*trimax = 2 * (nv + 10);*/
   if ((v = (TRIANGLE *)malloc(trimax*(long)sizeof(TRIANGLE))) == NULL) 
   {
      //printf("TRIANGULATE - Unable to allocate enough memory\n");
      free(pxyz);
      return;
   }
	
   /* Read coordinates while computing the bounds */
   fptr = fopen(filename,"r");
   nv = 0;
   xmin =   INFINITY;	xmax = - INFINITY;
   ymin =   INFINITY;	ymax = - INFINITY;
   zmin =   INFINITY;	zmax = - INFINITY;
  while(fscanf(fptr,"%s %lf %lf %lf\n",str,&x,&y,&z) == 4)
   {
      pxyz[nv+1].x = x;
      pxyz[nv+1].y = y;
      pxyz[nv+1].z = z;
      nv++;
      xmin = MIN(xmin,x);	xmax = MAX(xmax,x);
      ymin = MIN(ymin,y);	ymax = MAX(ymax,y);
      zmin = MIN(zmin,z);	zmax = MAX(zmax,z);
   }
   fclose(fptr);

	
   /* 以X坐标先排序一下 */
   for (i=1;i<nv;i++) 
   {
      for (j=i+1;j<=nv;j++) 
	  {
         if (pxyz[j].x < pxyz[i].x) 
		 {
       	    ptemp   = pxyz[i];
            pxyz[i] = pxyz[j];
            pxyz[j] = ptemp;
         }
      }
   }

   /* 检查是否存在重复点，如果有则删除 */
   duplicate = 0;
   for (i=1;i<nv;i++) 
   {
      for (j=i+1;j<=nv;j++) 
	  {
         if (pxyz[j].x == pxyz[i].x) 
		 {
       	    if (pxyz[j].y == pxyz[i].y) 
			{
                duplicate++;
				pxyz[i] = pxyz[nv];
				nv--;
			 }
		  }
      }
   }


   XYZ  m_Scale,m_Center;  

   m_Center.x=(xmax+xmin)/2;
   m_Center.y=(ymax+ymin)/2;
  //m_Center.z=(zmax+zmin)/2;

   m_Scale.x=RANG/(xmax-xmin);
   m_Scale.y=RANG/(ymax-ymin);
   m_Scale.z=HIGH; //


   xmin =   INFINITY;	xmax = - INFINITY;
   ymin =   INFINITY;	ymax = - INFINITY;
   zmin =   INFINITY;	zmax = - INFINITY;

   for (i=1;i<nv;i++) 
   {
	  pxyz[i].x=m_Scale.x*(pxyz[i].x-m_Center.x);
	  pxyz[i].y=m_Scale.y*(pxyz[i].y-m_Center.y);
	  pxyz[i].z=m_Scale.z*pxyz[i].z;

      xmin = MIN(xmin,pxyz[i].x);	xmax = MAX(xmax,pxyz[i].x);
      ymin = MIN(ymin,pxyz[i].y);	ymax = MAX(ymax,pxyz[i].y);

   }



   triangulate(nv,pxyz,v,&ntri);  
   
   dx = (xmax - xmin) / xcells;
   dy = (ymax - ymin) / ycells;

   
   nPlane=0;
   for (j=0;j<=ycells;j++) 
   {			
      p[0].y = ymin + j * dy;
      p[1].y = p[0].y;
      p[2].y = p[0].y + dy;
      p[3].y = p[2].y;
	  	
      for (i=0;i<=xcells;i++) 
	  {
          p[0].x = xmin + i * dx;
          p[1].x = p[0].x + dx;
    	  p[2].x = p[1].x;
	      p[3].x = p[0].x;
	       		
          for (corner=0;corner<4;corner++) 
		  {
				found[corner] = FALSE;

				for (k=1;k<=ntri;k++) 
				{
					p1 = pxyz[v[k].p1];
					p2 = pxyz[v[k].p2];
					p3 = pxyz[v[k].p3];

					if (intriang(p[corner].x,p[corner].y,p1,p2,p3)) 
					{
						p[corner].z = planepoint(p[corner].x,p[corner].y,p1,p2,p3);
						found[corner] = TRUE;
						break;
					}
				}
	   }
	       		
	   if (found[0] && found[1] && found[2] && found[3]) //面的四个点都找到位置
	       AddFace(p);
       }
   }

   free(v);
   free(pxyz);
}

/*****************************************************************************
	Triangulation subroutine
	
	Takes as input NV vertices in array pxyz
	Returned is a list of ntri triangular faces in the array v
	These triangles are arranged in clockwise order.
	*/   

/* 生成 delaunay 三角形 */
/*
		input : nv 点的数目
		pxyz 点数据
		output :
		v    三角形数据
		ntri 三角形数目
*/

void triangulate(int nv,XYZ *pxyz,TRIANGLE *v,unsigned int *ntri)
{
	int *complete;
	EDGE *edges;
	int nedge;
	
	int inside;
	int i,j,k;
	double xp,yp,x1,y1,x2,y2,x3,y3,xc,yc,r;
	double xmin,xmax,ymin,ymax,xmid,ymid;
	double dx,dy,dmax;
	
	/* Allocate memory for the edge list */
	if ((edges = (EDGE *)malloc(EMAX*(long)sizeof(EDGE))) == NULL) {
	   //printf("TRIANGULATE - unable to allocate enough memory\n");
	   return;
	}
	
	/* Allocate memory for the completeness list, flag for each triangle */
	if ((complete = (int *)malloc(trimax*(long)sizeof(int))) == NULL) {
	   //printf("TRIANGULATE - unable to allocate enough memory\n");
	   free(edges);
	   exit(0);
	}
	
	/*
	   Find the maximum and minimum vertex bounds.
	   This is to allow calculation of the bounding triangle
	*/
	xmin = pxyz[1].x;
	ymin = pxyz[1].y;
	xmax = xmin;
	ymax = ymin;
	for (i=2;i<=nv;i++) {
	   if (pxyz[i].x < xmin) xmin = pxyz[i].x;
	   if (pxyz[i].x > xmax) xmax = pxyz[i].x;
	   if (pxyz[i].y < ymin) ymin = pxyz[i].y;
	   if (pxyz[i].y > ymax) ymax = pxyz[i].y;
	}
	dx = xmax - xmin;
	dy = ymax - ymin;
	dmax = (dx > dy) ? dx : dy;
	xmid = (xmax + xmin) / 2.0;
	ymid = (ymax + ymin) / 2.0;
	
	/*
	   Set up the supertriangle
	   This is a triangle which encompasses all the sample points.
	   The supertriangle coordinates are added to the end of the 
	   vertex list. The supertriangle is the first triangle in
	   the triangle list.
	*/
	pxyz[nv + 1].x = xmid - 2.0 * dmax;
	pxyz[nv + 1].y = ymid - dmax;
	pxyz[nv + 1].z = 0.0;
	pxyz[nv + 2].x = xmid;
	pxyz[nv + 2].y = ymid + 2.0 * dmax;
	pxyz[nv + 2].z = 0.0;
	pxyz[nv + 3].x = xmid + 2.0 * dmax;
	pxyz[nv + 3].y = ymid - dmax;
	pxyz[nv + 3].z = 0.0;

	v[1].p1 = nv + 1;
	v[1].p2 = nv + 2;
	v[1].p3 = nv + 3;
	complete[1] = FALSE;
	*ntri = 1;

	/*
	   Include each point one at a time into the existing mesh
	*/
	for (i=1;i<=nv;i++) 
	{
		
	   xp = pxyz[i].x;
	   yp = pxyz[i].y;
	   nedge = 0;

	   /*
	      Set up the edge buffer. 
       	      If the point (xp,yp) lies inside the circumcircle then the
	      three edges of that triangle are added to the edge buffer
	      and that triangle is removed.
	   */
	   j = 0;
	   do {
	      j++;
	      if (!complete[j]) {
	         x1 = pxyz[v[j].p1].x;
	         y1 = pxyz[v[j].p1].y;
	         x2 = pxyz[v[j].p2].x;
	         y2 = pxyz[v[j].p2].y;
	         x3 = pxyz[v[j].p3].x;
	         y3 = pxyz[v[j].p3].y;
	         inside = circumcircle(xp,yp,x1,y1,x2,y2,x3,y3,&xc,&yc,&r);
	         if (xc + r < xp)
	       	    complete[j] = TRUE;
	            if (inside) {
	               /* Check that we haven't exceeded the edge list size */
		       if (nedge+3 > EMAX) {
	     	          //printf("TRIANGULATE - Internal edge list exceeded\n");
              	       	  exit(0);
       	       	       }
		       edges[nedge+1].p1 = v[j].p1;
		       edges[nedge+1].p2 = v[j].p2;
		       edges[nedge+2].p1 = v[j].p2;
		       edges[nedge+2].p2 = v[j].p3;
		       edges[nedge+3].p1 = v[j].p3;
		       edges[nedge+3].p2 = v[j].p1;
		       nedge += 3;
		       v[j] = v[*ntri];
		       complete[j] = complete[*ntri];
		       j--;
		       (*ntri)--;
		    }
 		}
	    } while (j < *ntri);

	    /*
	        Tag multiple edges
	        Note: if all triangles are specified anticlockwise then all
	              interior edges are opposite pointing in direction.
	    */
	    for (j=1;j<nedge;j++) 
		{
	        for (k=j+1;k<=nedge;k++) 
			{
	       	   if ((edges[j].p1 == edges[k].p2) && (edges[j].p2 == edges[k].p1)) 
			   {
					edges[j].p1 = 0;	
					edges[j].p2 = 0;
					edges[k].p1 = 0;	
					edges[k].p2 = 0;
			   }

				/* Shouldn't need the following, see note above */
				if ((edges[j].p1 == edges[k].p1) && (edges[j].p2 == edges[k].p2)) 
				{
						edges[j].p1 = 0;	
						edges[j].p2 = 0;
						edges[k].p1 = 0;	
						edges[k].p2 = 0;
				}
			}  
	    }

	    /*
	       Form new triangles for the current point
	       Skipping over any tagged edges.
	       All edges are arranged in clockwise order.
	    */
	    for (j=1;j<=nedge;j++) 
		{
	       if (edges[j].p1 != 0 && edges[j].p2 != 0) 
		   {
	          if (*ntri > trimax) 
			  {
				//printf("TRIANGULATE - %d triangles exceeds maximum\n",*ntri);
	             exit(0);
			  }

			(*ntri)++;
			v[*ntri].p1 = edges[j].p1;
			v[*ntri].p2 = edges[j].p2;
			v[*ntri].p3 = i;
			complete[*ntri] = FALSE;
	       }
	   }
       }

       /*
	   Remove triangles with supertriangle vertices
	   These are triangles which have a vertex number greater than nv
       */
       i = 0;
       do 
	   {
			i++;
			if ((v[i].p1 > nv) || (v[i].p2 > nv) || (v[i].p3 > nv)) 
			{
				v[i] = v[*ntri];
				i--;
				(*ntri)--;
			}
	   } while (i < *ntri);
	
    free(edges);
	free(complete);
}

/*****************************************************************************
	Return TRUE if a point (xp,yp) is inside the circumcircle made up
	of the points (x1,y1), (x2,y2), (x3,y3)
	The circumcircle centre is returned in (xc,yc) and the radius r
	NOTE: A point on the edge is inside the circumcircle
*/
int circumcircle(double xp,double yp,
				 double x1,double y1,
				 double x2,double y2,
				 double x3,double y3,
				 double *xc,double *yc,
				 double *r)
{
   double m1,m2,mx1,mx2,my1,my2;
   double dx,dy,rsqr,drsqr;

   if (ABS(y1-y2) < EPSILON && ABS(y2-y3) < EPSILON) 
   {
       //printf("TRIANGULATE - Coincident points\n");
       return(FALSE);
   }

   if (ABS(y2-y1) < EPSILON) 
   {
       	m2 = - (x3-x2) / (y3-y2);
       	mx2 = (x2 + x3) / 2.0;
       	my2 = (y2 + y3) / 2.0;
       	*xc = (x2 + x1) / 2.0;
       	*yc = m2 * (*xc - mx2) + my2;
    } 
   else if (ABS(y3-y2) < EPSILON) 
   {
       	m1 = - (x2-x1) / (y2-y1);
       	mx1 = (x1 + x2) / 2.0;
       	my1 = (y1 + y2) / 2.0;
       	*xc = (x3 + x2) / 2.0;
       	*yc = m1 * (*xc - mx1) + my1;
    } 
   else 
   {
       	m1 = - (x2-x1) / (y2-y1);
       	m2 = - (x3-x2) / (y3-y2);
       	mx1 = (x1 + x2) / 2.0;
       	mx2 = (x2 + x3) / 2.0;
       	my1 = (y1 + y2) / 2.0;
       	my2 = (y2 + y3) / 2.0;
       	*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;
	
    return((drsqr <= rsqr) ? TRUE : FALSE);
}

/*****************************************************************************
   Returns TRUE if the point (xp,yp) lies inside the projection triangle
   on the x-y plane with the vertices p1,p2,p3

   A point is in the centre if it is on the same side of all the edges
   or if it lies on one of the edges.
*/
int intriang(double xp,double yp,XYZ p1,XYZ p2,XYZ p3)
{
   int side1,side2,side3;

   side1 = whichside(xp,yp,p1.x,p1.y,p2.x,p2.y);
   side2 = whichside(xp,yp,p2.x,p2.y,p3.x,p3.y);
   side3 = whichside(xp,yp,p3.x,p3.y,p1.x,p1.y);
	
   if (side1 == 0 && side2 == 0)
      return(TRUE);
   if (side1 == 0 && side3 == 0)
      return(TRUE);
   if (side2 == 0 && side3 == 0)
      return(TRUE);

   if (side1 == 0 && (side2 == side3))
      return(TRUE);
   if (side2 == 0 && (side1 == side3))
      return(TRUE);
   if (side3 == 0 && (side1 == side2))
      return(TRUE);
	
   if ((side1 == side2) && (side2 == side3))
      return(TRUE);
		
   return(FALSE);
}

/*****************************************************************************
   Given a plane through the point p1,p2,p3 and a point on the x-y
   plane determine the intersection of the perpendicular.
*/

/*
	ax+by+cz+d=0 平面方程
    所以在位置(xp,yp)出的z为:
	z=- (a * xp + b * yp + d) / c
*/
double planepoint(double xp,double yp,XYZ p1,XYZ p2,XYZ p3)
{
   double a,b,c,d;

   a = p1.y * (p2.z - p3.z) + p2.y * (p3.z - p1.z) + p3.y * (p1.z - p2.z);
   b = p1.z * (p2.x - p3.x) + p2.z * (p3.x - p1.x) + p3.z * (p1.x - p2.x);
   c = p1.x * (p2.y - p3.y) + p2.x * (p3.y - p1.y) + p3.x * (p1.y - p2.y);
   d = - p1.x * (p2.y * p3.z - p3.y * p2.z) 
       - p2.x * (p3.y * p1.z - p1.y * p3.z) 
       - p3.x * (p1.y * p2.z - p2.y * p1.z);
	
   if (ABS(c) > EPSILON)
      return(- (a * xp + b * yp + d) / c);
   else
      return(0.0);
}

/*****************************************************************************
   Determines which side of a line the point (x,y) lies.
   The line goes from (x1,y1) to (x2,y2)

   Return codes are 	-1 for points to the left
		      	         0 for points on the line
		       	        +1 for points to the right
*/
int whichside(double x,double y,double x1,double y1,double x2,double y2)
{
    double dist;

    dist = (y - y1) * (x2 - x1) - (y2 - y1) * (x - x1);
    if (dist > 0)
       return(-1);
    else if (dist < 0)
       return(1);
    else
       return(0);
}

void AddFace(XYZ *p)
{

  demGrid[nPlane].P0.x=p[0].x;
  demGrid[nPlane].P0.y=p[0].y;
  demGrid[nPlane].P0.z=p[0].z;

  demGrid[nPlane].P1.x=p[1].x;
  demGrid[nPlane].P1.y=p[1].y;
  demGrid[nPlane].P1.z=p[1].z;

  demGrid[nPlane].P2.x=p[2].x;
  demGrid[nPlane].P2.y=p[2].y;
  demGrid[nPlane].P2.z=p[2].z;

  demGrid[nPlane].P3.x=p[3].x;
  demGrid[nPlane].P3.y=p[3].y;
  demGrid[nPlane].P3.z=p[3].z;

  nPlane++;
}