tri_triangulation.cpp

三角网的生成算法

	if( nPoints <= i || nPoints <= j ) return false;

	if( !bCheckCone_TRI( Pnt2List, i, j ) ) return false;
	return bCheckDiagonalie_TRI( Pnt2List, i, j );
}


//----------------------------------------------------------------------
// devide polygon into triangles
// 
// Returns:
// true or false
//
// Notes:
// polygon must be counter-clockwise
//
int iTriangulatePolygon_TRI( 
		vector<Cp_dbl2> *Pnt2List,	// [i/o] polygon
		vector<Cp_dtri> *Tri2List )	// [i/o] list of triangles
{
	int i, j, k, nTri, nPoints;
	Cp_dtri tri;

	nPoints = (int)Pnt2List->size( );
	if( 3 >= nPoints ) goto PIX_EXIT;

	for( i = 0 ; i < nPoints ; i++ ) {
		j = ( i+1 ) % nPoints;
		k = ( i+2 ) % nPoints;
		if( bCheckDiagonal_TRI( Pnt2List, i, k ) ) {
			tri.Set( (int)Tri2List->size( ), Pnt2List->at( i ).GetID( ), 
				Pnt2List->at( j ).GetID( ), Pnt2List->at( k ).GetID( ) );
			Tri2List->push_back( tri );
			Pnt2List->erase( Pnt2List->begin( ) + j );
			nTri = iTriangulatePolygon_TRI( Pnt2List, Tri2List );
			break;
		}
	}

	// --- Done ---
PIX_EXIT:
	return nPoints;
}


//----------------------------------------------------------------------
// devide polygon into triangles
// 
// Returns:
// true or false
//
// Notes:
// polygon must be counter-clockwise
//

void StartPolygonTri_TRI(
		vector<Cp_dbl2> *Pnt2List,	// [i]
		vector<Cp_dtri> *Tri2List )
{
	int nTri;
	Cp_dbl2 ep;
	vector<Cp_dbl2> PntList;
	vector<Cp_dbl2>::iterator Iter;

	// --- copy the polygon ---
	PntList.clear( );
	for ( Iter = Pnt2List->begin( ); Iter != Pnt2List->end( ); Iter++ ) {
		ep = *Iter;
		PntList.push_back( ep );
	}

	// --- do the triangulation ---
	Tri2List->clear( );
	nTri = iTriangulatePolygon_TRI( &PntList, Tri2List );

	// --- Done ---
	
	PntList.clear( );
	return;
}


/**********************************************************************/
/**********************************************************************/
/*** triangulated Irregular Network ***********************************/
/**********************************************************************/
/**********************************************************************/

//======================================================================
// 
//
bool bInitNeighbourList( 
		int				npoints,
		CDT_Matrix<int>	*nbr )
{
	int i, j, n;
	n = 0;
	for ( i = 0; i < npoints; i++ ) {
		for( j = 0; j < npoints; j++ ) {
			nbr->Set( i, j, n );
		}
	}

	return true;
}


//======================================================================
// prepare PSEUDO_POINT_COUNT pseudo points
//
bool bSetPseudoPoints_TRI( 	
		int				opoints,	// [i] number of original points
		vector<Cp_dbl2> *PntList )	// [i/o] list of input points (including pseudo points )
{
	int		i;
	double	xsize, ysize;
	Cp_dbl2	ep, p_min, p_max;

	/*------------------------------------------------------------------
	#### pseudo points.
	#### prepare PSEUDO_POINT_COUNT pseudo points which are added around 
	#### the input to points to suppport the quick configulation of TIN. 
	#### points are located at the following positions.
	####  X--------------X--------------X--------------X ... p_max.y
	####  |3              4              5            6|      /\
	####  |                                            |      ||
	####  |        +----------------*---------+        | ... p_max.y
	####  |        |                        * |        |
	####  |        |     *                    |        |
	####  X 2      |                          |       7X
	####  |        |           *              *        |
	####  |        |                          |        |
	####  |        |       *           *      |        |
	####  |        |                          |        |
	####  |        |                *         |        |
	####  |        |      *                   |        |
	####  X 1      |                          |       8X
	####  |        |                *       * |        |
	####  |        |  *                       |        |
	####  |        +-----*--------------------+        | ... p_min.y
	####  |      Minimum rectangle which includes      |      ||
	####  |0       .     all input points.    .       9|      \/
	####  X--------------X--------------X--------------X ... p_min.y
	####  .        .     11             10    .        .
	####  . <===== .                          . =====> .
	####  .        .                          .        .
	####  p_min.x                          p_max.x    
	####           p_min.x                          p_max.x 
	------------------------------------------------------------------*/
	PixMiniMax( opoints, PntList, &p_min, &p_max );
	xsize = p_max.x - p_min.x;
	ysize = p_max.y - p_min.y;
	
	p_min.x = p_min.x - xsize / 2;
	p_min.y = p_min.y - ysize / 2;
	p_max.x = p_max.x + xsize / 2;
	p_max.y = p_max.y + ysize / 2;

	i = opoints;
	ep.Set( i, p_min.x, p_min.y );
	PntList->push_back( ep );

	i++;
	ep.Set( i, p_min.x, p_min.y + ( p_max.y - p_min.y ) / 3 );
	PntList->push_back( ep );

	i++;
	ep.Set( i, p_min.x, p_min.y + ( p_max.y - p_min.y ) * 2 / 3 );
	PntList->push_back( ep );

	i++;
	ep.Set( i, p_min.x, p_max.y );
	PntList->push_back( ep );

	i++;
	ep.Set( i, p_min.x + ( p_max.x - p_min.x ) / 3, p_max.y );
	PntList->push_back( ep );

	i++;
	ep.Set( i, p_min.x + ( p_max.x - p_min.x ) * 2 / 3, p_max.y );
	PntList->push_back( ep );

	i++;
	ep.Set( i, p_max.x, p_max.y );
	PntList->push_back( ep );

	i++;
	ep.Set( i, p_max.x, p_min.y + ( p_max.y - p_min.y ) * 2 / 3 );
	PntList->push_back( ep );

	i++;
	ep.Set( i, p_max.x, p_min.y + ( p_max.y - p_min.y ) / 3 );
	PntList->push_back( ep );

	i++;
	ep.Set( i, p_max.x, p_min.y );
	PntList->push_back( ep );

	i++;
	ep.Set( i, p_min.x + ( p_max.x - p_min.x ) * 2 / 3, p_min.y );
	PntList->push_back( ep );

	i++;
	ep.Set( i, p_min.x + ( p_max.x - p_min.x ) / 3, p_min.y );
	PntList->push_back( ep );

	return true;
}


int iListNeighbour_TRI(
		int 			p,				// [i] point currently being looked at
		int 			opoints, 
		vector<Cp_dbl2> *PntList,
		int 			*iBelngCount,
		CDT_Matrix<int>	*iNbrList,
		CDT_Matrix<int>	*iBlgList, 
		vector<Cp_dtri> *TriVector )
{
	bool bRet;
	int i, j, k, nneib, ptr0, ptr1, ptr2, ptr3, ninside;
	int bt, npoints, nright, ptr2_in_trilist;
	short sStatus;
	double dMaxAngle, dRadiusStep, dRadius, dDistance, angle, angle2, dx, dy;
	Cp_dtri tri;

	int *iPntSorted = NULL;
	int *iRight = NULL;
	int *iIncluded = NULL;

	Cp_Circle2	Circl, Circl2;
	Cp_dbl2		ep, t_cross, t_midd, t_plus;
	Cp_lin2		TmpLine, line_ptr01;

	/*------------------------------------------------------------------
	#### iIncluded[i]	TRUE if i-th points ( include the pseudo points )  
	#### 				is already in the neighbour list around the 
	#### 				input point "p".
	#### iRight[i]		TRUE if i-th points ( include the pseudo points ) 
	#### 				is on the "right" side of the line connecting 
	#### 				"p" ( neighbour list[p][0] ) and first connected  
	#### 				point in the neighbour list. FALSE if not on the 
	#### 				right side or unknown.
	------------------------------------------------------------------*/
	
	npoints = (int)PntList->size( ); 

	iIncluded = new int[npoints];
	if( NULL == iIncluded ) {
		nneib = -1;
		goto PIX_EXIT;
	}

	iRight = new int[npoints];
	if( NULL == iRight ) {
		nneib = -1;
		goto PIX_EXIT;
	}

	for ( i = 0; i < npoints; i++ ) {
		iIncluded[i] = false;
		iRight[i] = false;
	}

	// --- create a sorted list according to the distance from p ---
	
	iPntSorted = new int[npoints];
	bRet = bSortByDistance( p, PntList, iPntSorted );
	if( !bRet ) goto PIX_EXIT;

	/*------------------------------------------------------------------
	#### set "p" itself as the first point in the neighbour list around 
	#### "p". it is also set to "ptr0".
	------------------------------------------------------------------*/
	ptr0 = p;
	nneib = 0;
	iNbrList->Set( p, nneib, ptr0 );
	iIncluded[ptr0] = true;
	nneib++;

	/*------------------------------------------------------------------
	#### decide the second point to be connected to the first. the 
	#### point is decided to be the second if no other point is 
	#### included in the circle, of which diamter is defined by the 
	#### first and the current point. 
	#### 
	#### any point which first meets this condition will be the second 
	#### point. the second point obtained is named "ptr2" and it is set 
	#### in the "neighbour list". 
	#### 
	#### ptr0: first point.
	#### ptr1: first neighbour.
	#### get_circle_end: get and return the circle in the two dimensional space, of which 
	#### two end of the diamter are on two given points.
	------------------------------------------------------------------*/
	
	ptr1 = iPntSorted[1];
	if( ptr1 >= npoints ) ptr1 = 0;
	Circl.GetCircleEnd( &PntList->at( ptr0 ), &PntList->at( ptr1 ) );

	while( true ) {
		ninside = 0;
		for ( i = 0; i < npoints; i++ ) {
			k = iPntSorted[i];
			if( k == ptr0 ) continue;
			if( k == ptr1 ) continue;
			if( Circl.bPointInCircle( &PntList->at( k ) ) ) {
				ptr1 = k;
				Circl.GetCircleEnd( &PntList->at( ptr0 ), &PntList->at( ptr1 ) );
				ninside++;
				break;
			}
		}

		if( ninside == 0 ) break;
	}

	iNbrList->Set( p, nneib, ptr1 );
	iIncluded[ptr1] = true;
	nneib++;
	
	/*------------------------------------------------------------------
	#### t_plus mugen enten
	------------------------------------------------------------------*/
	t_plus.x = PntList->at( ptr0 ).x + 10000000.0;
	t_plus.y = PntList->at( ptr0 ).y;

	while( true ) {
		/*--------------------------------------------------------------
		#### look for the triangle which include "ptr0" as vertex in 
		#### the triangle list around "ptr1". if found, and if the 
		#### third vertex of the triangle is not in the neighbour list 
		#### of ptr0, and it also is on the right side of the vector 
		#### connecting ptr0 and ptr1 ( it is equivalent to that the 
		#### angle defined by ptr0 -> ptr1 -> ( the third vertex ) is 
		#### smaller than PI ), the third vertex is decided to 
		#### be the next point in the neighbour list. 
		--------------------------------------------------------------*/
		ptr2 = -1;
		ptr2_in_trilist = false;
		for ( i = 0; i < iBelngCount[ptr1]; i++ ) {
			ptr3 = -1;
		
			bt = iBlgList->Get( ptr1, i );			
			tri = TriVector->at( bt );
	
			if( tri.GetPnt(0) == ptr0 ) {
				if( tri.GetPnt(1) == ptr1 ) {
					ptr3 = tri.GetPnt(2);
				} else if( tri.GetPnt(2) == ptr1 ) {
					ptr3 = tri.GetPnt(1);
				}
			} else if( tri.GetPnt(1) == ptr0 ) {
				if( tri.GetPnt(0) == ptr1 ) {
					ptr3 = tri.GetPnt(2);
				} else if( tri.GetPnt(2) == ptr1 ) {
					ptr3 = tri.GetPnt(0);
				}
			} else if( tri.GetPnt(2) == ptr0 ) {
				if( tri.GetPnt(0) == ptr1 ) {
					ptr3 = tri.GetPnt(1); 
				} else if( tri.GetPnt(1) == ptr1 ) {
					ptr3 = tri.GetPnt(0);
				}
			}

			if( ptr3 < 0 ) continue;
			if( bCmpNbrList( p, 1, ptr3, iNbrList ) && !bCmpNbrList( p, 2, ptr1, iNbrList ) ) {
				ptr2_in_trilist = true;
				ptr2 = ptr3;
				break;
			}

			if( iIncluded[ptr3] ) continue;
			angle = get_3angle( &PntList->at( ptr0 ), &PntList->at( ptr1 ), &PntList->at( ptr3 ) );
			if( angle < PIE ) {
				ptr2_in_trilist = true;
				ptr2 = ptr3;
				break;
			}
		}

		/*--------------------------------------------------------------
		#### searching for virgins. if no connecting point is found 
		#### using the above ( triangle ) step, try to find it from the 
		#### gemetrical conditions.
		--------------------------------------------------------------*/
		if( ptr2 < 0 ) {
			/*----------------------------------------------------------
			#### check for all points, if it is on the right side of 
			#### the vector connecting "ptr0"->"ptr1".
			----------------------------------------------------------*/
			nright = 0;
			line_ptr01.PlaneLine( &PntList->at( ptr0 ), &PntList->at( ptr1 ) );
			for ( j = 0; j < npoints; j++ ) {
				iRight[j] = false;
				if( j == ptr0 ) continue;
				if( j == ptr1 ) continue;
				if( RIGHT == line_ptr01.iWhichSideOfLine( &PntList->at( j ) ) ) {
					nright++;
					iRight[j] = true;
				}
			}

			/*----------------------------------------------------------
			#### if nopoint is found on the right side of "ptr0" -> 
			#### "ptr1", stop the searching.
			----------------------------------------------------------*/
			if( nright <= 0 ) break;

			/*----------------------------------------------------------
			### define a straight line to pass the center of ptr0 and 
			### ptr1 at the right angle to ptr0->ptr1 
			----------------------------------------------------------*/
			angle = get_3angle( &t_plus, &PntList->at( ptr0 ), &PntList->at( ptr1 ) );
			t_midd.Center( &PntList->at( ptr0 ), &PntList->at( ptr1 ) );
			angle2 = GoodRadian( angle - PIE05 );
			TmpLine.OnePointAngleLine( &t_midd, angle2 );
			dRadius = dbl2_distance( &PntList->at( ptr0 ), &PntList->at( ptr1 ) ) / 2;
			//dDegStep = 0.0;
			dRadiusStep = dRadius * DISTANCE_PITCH;
			dDistance = 0.0;
			
			/*----------------------------------------------------------
			----------------------------------------------------------*/
			sStatus = true;
			while( sStatus ) {
				/*------------------------------------------------------
				#### define a circle which pass "ptr0" and "ptr1" which 
				#### is "little bit" widened to the right of "ptr0" -> 
				#### "ptr1" fromn the previous circle.
				------------------------------------------------------*/
				dDistance = dDistance + dRadiusStep;
				t_cross = TmpLine.PointFromLineStart( dDistance );
				dx = PntList->at( ptr0 ).x - t_cross.x;
				dy = PntList->at( ptr0 ).y - t_cross.y;
				Circl2.Set( t_cross, sqrt( dx * dx + dy * dy ) );

				/*------------------------------------------------------
				#### look for the point which is on the right side of 
				#### the vector and included in the just defined circle.
				#### among these points, select the one which minimize 
				#### the angle "ptr0"->"new point"->"ptr1".
				------------------------------------------------------*/
				ninside = 0;
				ptr2 = -1;
				for ( i = 0; i < npoints; i++ ) {
					if( i == ptr0 ) continue;
					if( i == ptr1 ) continue;
					if( !iRight[i] ) continue;
					if( !Circl2.bPointInCircle( &PntList->at( i ) ) ) continue;
					angle = get_3angle( &PntList->at( ptr1 ), &PntList->at( i ), &PntList->at( ptr0 ) );
					if( 0 == ninside ) {
						dMaxAngle = angle;
						ptr2 = i;
						ninside++;
					} else if( angle >= dMaxAngle ) {
						ptr2 = i;
						dMaxAngle = angle;
						ninside++;
					}
				}

				/*------------------------------------------------------
				#### if any point which can be connected is found, break. 
				#### if not, restart the step with widened circle.
				------------------------------------------------------*/
				if( ptr2 >= 0 ) sStatus = false;
			}
		}

		/*--------------------------------------------------------------
		#### welcome back!!!
		#### if the next connected point is the second point in the 
		#### neighbour list, stop the search and exit from the loop.
		--------------------------------------------------------------*/