Remove Border Triangles

How To Remove Badly Shaped Border Triangles

A Delaunay triangulation is convex but that does not always hold for your data. As a consequence your triangulation may contain triangles close to the convex hull that do not really belong to the desired non-convex result. This post shows how to remove them.

Minimal Example in C++

A small example shall demonstrate the problem and the solution.

The Good Case

Let’s start with a good case: 6 points that have a certain elevation (z) are triangulated, a Zone2 object is constructed from the triangulation and as expected it holds exactly 4 triangles.

// * 1 *   Create 6 surface points
 vector<Point2> vSurfacePoints;
 vSurfacePoints.push_back(Point2(0,0,10));
 vSurfacePoints.push_back(Point2(0,50,0));
 vSurfacePoints.push_back(Point2(0,100,10));
 vSurfacePoints.push_back(Point2(20,0,10));
 vSurfacePoints.push_back(Point2(20,50,0));
 vSurfacePoints.push_back(Point2(20,100,10));

 // * 2 *   Triangulate and create a zone
 cout<<"Constructing pZone0"<<endl; 
 Fade_2D dt0; 
 dt0.insert(vSurfacePoints); 
 Zone2* pZone0(dt0.createZone(NULL,ZL_GLOBAL,false)); 

 // * 3 * Show  
 pZone0->showGeomview("pZone0.list","1 0 0 0.5");
 cout<<"Zone0, numT: "<<pZone0->getNumberOfTriangles()<<endl;

4 triangles, as expected


The Bad Case

Let’s move one of the above points slightly inwards, triangulate and create another Zone2 object again. Now the zone holds 5 triangles.

	vSurfacePoints[1]=Point2(1e-5,50,0);

	// * 5 *   Triangulate and create a zone again
	Fade_2D dt1;
	dt1.insert(vSurfacePoints);
	Zone2* pZone1(dt1.createZone(NULL,ZL_GLOBAL,false));

	// * 6 *   Show
	pZone1->showGeomview("pZone1.list","0 1 0 0.5");
	cout<<"Zone1, numT: "<<pZone1->getNumberOfTriangles()<<endl;

The additional triangle has almost no 2D area (area in the XY plane) but its 2.5D area may be significant because such triangles can deviate by almost 90 degrees from the XY plane.

An additional triangle appears


The Solution

The solution takes a rule according to which triangles shall be removed. But clearly there can’t be one general rule. Instead you can define your own rule by means of the UserPredicateT class and apply this rule to the triangulation. It’s easy. Let’s derive a class from UserPredicateT and implement a custom rule there. The rule is: “Any border triangle that deviates by more than 88 degrees from the XY plane shall be removed”.

class PeelPredicate:public UserPredicateT
{
public:
	bool operator()(const Triangle2* pT)
	{
		Vector2 nv(pT->getNormalVector());
		Vector2 up(0,0,1);

		double angle(-1);
		double cosPhi(nv*up);
		if(cosPhi>1.0) angle=0; // Robustness in case of numeric inaccuracy
			else if(cosPhi<-1.0) angle=180.0; // Robustness in case of numeric inaccuracy
				else angle=acos(cosPhi)*180/3.14159;

		if(angle>88.0) return true;
			else return false;
	}
};

...

	// * 7 *   Peel unwanted border triangles off
	PeelPredicate pp;
	Zone2* pZone1peeled=peelOffIf(pZone1,&pp,false);

	pZone1peeled->showGeomview("pZone1peeled.list","0 0 1 0.5");
	cout<<"Zone1peeled, numT: "<<pZone1peeled->getNumberOfTriangles()<<endl;

The additional triangle has been removed

Details

The peelOffIf() function removes only border triangles. But former interior triangles may become border triangles in a later step and thus the rule is applied iteratively. Above we have used the deviation of a triangle to the XY plane. Such a value does not exist for (planar) 2D triangulations but other criteria or even a combination of them can be used: The minimum interior angle, the 2D area of a triangle or its maximum edge length.

Spread the word. Share this post!

Leave A Reply

Your email address will not be published. Required fields are marked *

*

By continuing to use the site, you agree to the use of cookies. more information

The cookie settings on this website are set to "allow cookies" to give you the best browsing experience possible. If you continue to use this website without changing your cookie settings or you click "Accept" below then you are consenting to this.

Close