Delaunay Triangulation (TIN) with Fade2.5D

Create/Read a Surface Point Cloud

Create random surface points. Alternatively points can be read from an ASCII file.

void createTerrain(Fade_2D& dt)
{
  std::vector<Point2> vInputPoints;

  // * 1 *   Create random terrain points
  int seed(1); // seed>0 is for repeatable randomness
  generateRandomSurfacePoints(
    150, // numPointsX
    100, // numPointsY
    2, // numCenters
    0,0,0,100,100,50, // Bounds xmin,ymin,zmin,xmax,ymax,zmax
    vInputPoints, // Output vector
    seed  // Seed
  );

  // * 1 *   Alternative: Read terrain points from an ASCII file
  //readXYZ("yourFileName.xyz",vInputPoints);

Optional Point Cloud Pruning

As an optional step the point cloud is pruned in order to represent the TIN more efficiently. This is achieved by a call to EfficientModel::extract(double maxError, …). Thereby maxError specifies the allowed height error between the original and the pruned surface.

Point Cloud Pruning, Left: Original data, 15 000 points. Right: Pruned version, 655 points

  // * 2 *   Optional pruning
  cout<<"\nOriginal number of points: "<<vInputPoints.size()<<endl;
  EfficientModel em(vInputPoints);
  vInputPoints.clear();
  double maxError(.1);
  em.extract(maxError,vInputPoints);
  cout<<"Efficient number of points: "<<vInputPoints.size()<<endl<<endl;

TIN Visualization

The input points are inserted into a triangulation (Fade_2D object). Then, for development purposes it is often useful to quickly visualize the current state and thus Fade provides several visualization options: The common *.obj file format, visualization for web browsers, *.list files for the Geomview viewer (Linux) and 2D visualization using Postscript.

  // * 3 *   Insert the points
  dt.insert(vInputPoints);

  // * 4 *   Several options to visualize
  // Write a wavefront *.obj file
  dt.writeObj("terrain.obj");
  // Write a visualization for web browsers
  dt.writeWebScene("terrain_webScene");
  // Write a file for the Geomview viewer
  dt.showGeomview("terrain.list");
  // Write a postscript file (2D)
  dt.show("terrain.ps")

Breakline insertion

Constraint edges (‘Breaklines’) can be inserted into a TIN. Like in the 2D case (see Constraint Edges – Example3) segments can be inserted with or without subdivision. The purpose of subdivision is to achieve better shaped triangles (with larger minimal interior angles). Segments can be inserted at their original height (z-coordinate) or they can be adapted to the elevation of the surface. This gives four possible combinations, namely brealine insertion at…

• breakline elevation
• breakline elevation with subdivision
• surface elevation
• surface elevation with subdivision

Breakline insertion at breakline elevation

The simplest case is breakline insertion at breakline elevation. It requires just one call to Fade_2D::createConstraint(). ConstraintGraph2::makeDelaunay(double minLen) can additionally be called when subdivision is desired. Thereby the minLen parameter is used to define that segments smaller than minLen shall not be further subdivided. This helps to avoid an insane number of tiny triangles in narrow settings.

Breaklines at breakline elevation. Left without subdivision, right with subdivision

// A vector of line segments to be inserted
std::vector<Segment2> vSegments; 
...

// Insert vSegments without subdivision except when an existing vertex
// or another constraint segment is crossed.
ConstraintGraph2* pCG(NULL);
pCG=dt.createConstraint(vSegments,CIS_CONSTRAINED_DELAUNAY);

// If desired, subdivide the breaklines to achieve better triangles.
// Thereby minLen specifies that segments smaller than minLen shall
// not be further splitted.
double minLen(0.1);
pCG->makeDelaunay(minLen); // Create better shaped triangles

Breakline insertion at surface elevation

Inserting breaklines at surface elevation requires one preparatory function call: The segments (blue) must be draped onto the surface (red).

The method Fade_2D::drape(..) takes the original segments and drapes them (subdivided if necessary) onto the TIN. It takes a parameter double zTolerance. When zTolerance is negative then the original segments are splitted at all intersections with triangulation edges. Otherwise they are only subdivided when required to maintain the maximum height difference zTolerance between the segment and the surface. When the draped segments are ready they are inserted into the TIN using Fade_2D::createConstraint() as before, see the source code below.

Breaklines at surface elevation. Blue: Original segments. Red: Segments projected to the surface. Right: Subdivided to achieve better shaped triangles

std::vector<Segment2> vSegments; // Input segments
std::vector<Segment2> vDrapedSegments; // Draped segments
double zTolerance(0.1); // Allowed geometric error
dt.drape(vSegments,vDrapedSegments,zTolerance);

// Insert the draped segments
ConstraintGraph2* pCG(NULL);
pCG=dt.createConstraint(vDrapedSegments,CIS_CONSTRAINED_DELAUNAY);

// Subdivide if desired
double minLen(0.1);
pCG->makeDelaunay(minLen); // Create better shaped triangles

ISO contour computation

Iso contours are lines consisting of all points of a terrain having a certain height. Simply imagine that we cut with a horizontal plane through the height field. The below source code creates a global zone (consisting of all triangles), then it fetches the triangles from the zone and calls the isoContours() function.

Terrain triangulation 2.5D, ISO contours of 550 000 triangles at 30 different levels computed in 0.38 seconds.

	// 1) Create a global zone 
	Zone2* pZone(pDt->createZone(NULL,ZL_GLOBAL));	
	assert(pZone!=NULL);
	
	// 2) Get its triangles
	std::vector<Triangle2*> vTriangles;
	pZone->getTriangles(vTriangles);
	std::cout<<"Number of triangles: "<<vTriangles.size()<<std::endl;
	
	// 3) Compute ISO contours
	isoContours(vTriangles);

Function isoContours()

The isoContours(..) function creates an instance of the IsoContours class using the triangles of the global zone as argument. Then it determines three heights and calls getContours(height,...)

void isoContours(std::vector<Triangle2*>& vTriangles)
{
	// Create an instance of the IsoContours class and 3 arbitrary heights
	IsoContours isoManager(vTriangles);
	double minZ(isoManager.getMinHeight());
	double maxZ(isoManager.getMaxHeight());
	vector<double> vHeights;
	vHeights.push_back(minZ+(maxZ-minZ)*1/4);
	vHeights.push_back(minZ+(maxZ-minZ)*2/4);
	vHeights.push_back(minZ+(maxZ-minZ)*3/4);

	// Open a result file
	std::ofstream f("contours.txt");

	// For all (3) height values
	for(size_t i=0;i<vHeights.size();++i)
	{
		double height(vHeights[i]);
		
		// The result will be stored in a vector of polylines
		std::vector<std::vector<Segment2> > vvIsoContours;
	
		// IsoContours::getContours() intersects the triangulation with
		// a horizontal plane at height z. At certain heights degenerate
		// intersections exist and then IsoContours::getContours() 
		// returns false and a different height value must be chosen 
		// (small random perturbation is sufficient).
		while(!isoManager.getContours(height,vvIsoContours,true))
		{
			double rnd=(1e-4*rand()/(RAND_MAX+1.0));
			std::cout<<"Degenerate intersection at height "<<height<<std::endl;
			height+=rnd;
		}
	
		f<<"\n\n\n* NumContours at height "<<height<<": "<< \
                   vvIsoContours.size()<<"\n\n";
		
		// Write the result
		for(size_t i=0;i<vvIsoContours.size();++i)
		{
			std::vector<Segment2>& vContour(vvIsoContours[i]);
			f<<"Contour no. "<<i<<" at z="<<height<<" consists of "<<vContour.size()<<" segments"<<std::endl;
			Point2 sourcePoint(vContour[0].getSrc());
			Point2 targetPoint(vContour.back().getTrg());
			if(sourcePoint==targetPoint)
			{
				f<<"the contour is a closed polygon"<<std::endl;
			}
			else
			{
				f<<"the contour is an open polyline"<<std::endl;
			}
	
			for(std::vector<Segment2>::iterator it(vContour.begin());it!=vContour.end();++it)
			{
				Segment2& seg(*it);
				f<<seg<<std::endl;
			}
			f<<std::endl;
		}
	}
}

Height queries

Fade2.5D supports fast height queries with arbitrary (x,y) coordinate pairs. The getHeight(...) method locates the triangle that contains (x,y) and returns the height value z of the point. If the point (x,y) is outside the convex hull of the existing data points then getHeight(...) returns false.

double x(2390000);
double y(4569100);
double height;
bool hasHeight=pDt->getHeight(x,y,height);