# 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;
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

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
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.

# 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.