2.5D Terrain Triangulation (TIN)

Fade2.5D creates 2.5D triangulations from point clouds. These are like 2D Delaunay triangulations but with elevated vertices. 2.5D TINs are useful for lidar data and surface laser scans. The below example source code is contained in the download, see examples_25D/terrain.cpp. This example

  • creates points on a random surface (OR you can edit the example so that it reads your points from an (x y z)-ASCII file)
  • optionally prunes the point cloud
  • computes a 2.5D triangulation and visualizes it
  • inserts breaklines (constraint segments) into the TIN
  • computes ISO contours (polylines at a certain surface elevation level)
Delaunay Triangulation (TIN) with Fade2.5D

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);
Have a look at the module File I/O where various read/write functions for ASCII and binary files are documented.

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 could pruning to represent a TIN more efficiently

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 in a TIN, original breakline elevation

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 in a TIN, adapted to surface elevation

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

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);
Technically Fade2D and Fade2.5D are two independent libraries. But they share a common code base and the class and function names of Fade2D exist also in Fade2.5D. Thus they live in different namespaces, namely GEOM_FADE2D and GEOM_FADE25D.

Spread the word. Share this post!

4 Comments

  1. Reply Robert

    Hello,

    I am a R&D Engineer at Penn State interested in extracting iso-contours from triangulated meshes.
    We use CGAL in our research and would like to know more about terrain modeling and finding iso-contours from your experience.

    Would you be able to discuss the process for terrain/bathymetric analysis or to help guide our approach?

    We take in DEM data in an .xyz file, attain a 3D delaunay triangulation, “convert” the triangulation into a mesh, and use CGAL’s mesh-slicing functions. However this seems dreadfully slow. We are working on the scale of number of triangles (~500K) seen in your example, and would like to extract perhaps up to 100 contour levels at differing heights.

    Would you be able to lend your experience/insight?

    Thank you,
    Robert

    • Reply Bernhard Kornberger

      Hi,

      I am not sure whether I have correctly understood your approach: Do you really put 3D points into a 3D Delaunay triangulation? The output is a convex tetrahedral complex and even if you use a fake point at infinity to extract an approximate surface representation, that will be very slow for two reasons: First, a tetrahedral mesh has more elements than a triangular mesh. Second, the math is more complex, especially due to the high number of almost degenerate tetrahedra that DEM data probably causes.

      Fade2.5D should compute the ISO lines within a second. If your project is non-commercial I can give you a free research license, just send me an email with the details (research group, project leader, a short description of link to the project website)

      Best regards
      Bernhard

  2. Reply Jinming Zhang

    Hello,
    I am a phd student of Wuhan University in point cloud processing.
    Is there a function that can get all Vertices indice linked with a chosen point?

    • Reply Bernhard Kornberger

      Hello,

      you are right, such a function would probably be convenient. I can implement that for the next version. Meantime, let me provide this little piece of source:

      void incidentIndices()
      {
      	// * 1 *   Create random points
      	vector<Point2> vInputPoints;
      	generateRandomPoints(25,0,10,vInputPoints,1);
      
      	// * 2 *   Set indices
      	for(size_t i=0;i<vInputPoints.size();++i)
      	{
      		vInputPoints[i].setCustomIndex(i);
      	}
      
      	// * 3 *   Insert the points, retrieve the vertex handles
      	Fade_2D dt;
      	vector<Point2*> vHandles;
      	dt.insert(vInputPoints,vHandles);
      
      	// * 4 *   For a certain point handle, retrieve the incident triangles
      	Point2* pVtx(vHandles[0]);
      	vector<Triangle2*> vIncidentT;
      	dt.getIncidentTriangles(pVtx,vIncidentT);
      
      	// * 5 *   Iterate over the triangles to fetch the indicent points
      	set<Point2*> vIncidentVtx;
      	for(vector<Triangle2*>::iterator it(vIncidentT.begin());it!=vIncidentT.end();++it)
      	{
      		Triangle2* pT(*it);
      		for(int i=0;i<3;++i)
      		{
      			Point2* pCorner(pT->getCorner(i));
      			if(pCorner!=pVtx) vIncidentVtx.insert(pCorner);
      		}
      	}
      
      	// * 6 *   Print the indices
      	cout<<"Vertices incident to "<<pVtx->getCustomIndex()<<": ";
      	for(set<Point2*>::iterator it(vIncidentVtx.begin());it!=vIncidentVtx.end();++it)
      	{
      		Point2* pIncidentVtx(*it);
      		cout<<pIncidentVtx->getCustomIndex()<<" ";
      	}
      
      	// * 7 *   Draw it
      	Visualizer2 v("triangles.ps");
      	dt.show(&v);
      	v.addObject(Label(*pVtx,std::to_string(pVtx->getCustomIndex()),true,25),Color(CRED));
      	for(set<Point2*>::iterator it(vIncidentVtx.begin());it!=vIncidentVtx.end();++it)
      	{
      		Point2* pIncidentVtx(*it);
      		v.addObject(Label(*pIncidentVtx,std::to_string(pIncidentVtx->getCustomIndex()),true,25),Color(CGREEN));
      	}
      	v.writeFile();
      }
      

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