2.5D Terrain Triangulation

Fade2.5D is a 2D Delaunay triangulation for 3D points. Or more precisely, a 2D Delaunay triangulation with elevations at the vertices. It triangulates 3D points from terrains like lidar data, surface laser scans and other surface metrology techniques. The example source code described below is contained in the download, see examples_25D/terrain.cpp. This example

  • reads 3D points (x y z)-lines from an ASCII file
  • computes a 2.5D triangulation
  • visualizes the triangulated surface
  • computes ISO contours (polylines at a certain surface elevation level)
Delaunay Terrain with Fade2.5D

Delaunay Terrain with Fade2.5D


Fade2.5D and Fade2D are closely related but technically the two are independent libraries. The class and function names of Fade2D exist identically in Fade2.5D but the namespace is different, namely GEOM_FADE25D instead of GEOM_FADE2D.

Reading the terrain points

When using Visual Studio please run all examples in a command line window, not in a VS debug window. This example starts with a simple ASCII file reader for xyz data: Three numbers per line, whitespace-separated.

// Read an ASCII file containing xyz coordinates
void readTerrain(std::vector<Point2>& vInPoints,const char* filename)
{
    std::ifstream f(filename);
    if(f.is_open())
    {
        std::cout<<"reading "<<filename<<std::endl;
    }
    else
    {
        std::cout<<filename<<" is unreadable"<<std::endl; exit(1); } // For min/max statistics double minx(DBL_MAX),miny(DBL_MAX),minz(DBL_MAX); double maxx(DBL_MIN),maxy(DBL_MIN),maxz(DBL_MIN); double x,y,z; while(f>>x>>y>>z)
    {
        vInPoints.push_back(Point2(x,y,z));
        if(x<minx) minx=x;
        if(y<miny) miny=y;
        if(z<minz) minz=z; if(x>maxx) maxx=x;
        if(y>maxy) maxy=y;
        if(z>maxz) maxz=z;
    }

    std::cout<<vInPoints.size()<<" points"<<std::endl;
    std::cout<<"RANGE X:  "<<minx<<" to "<<maxx<<endl;
    std::cout<<"RANGE Y:  "<<miny<<" to "<<maxy<<endl;
    std::cout<<"RANGE Z:  "<<minz<<" to "<<maxz<<endl;
}

This common point file format does actually not require a handcrafted file reader like above. Have a look at the module File I/O where various read/write functions for ASCII and binary files are documented.

Triangulate the Terrain Points

The code below creates a new 2.5D triangulation. Then is reads points from an ASCII file. Then it inserts the points into the triangulation.

	Fade_2D* pDt=new Fade_2D();

	// 2) Read terrain points from an ASCII file (xyz-coordinates)
	std::vector<Point2> vInputPoints;
	readTerrain(vInputPoints,"../examples_25D/gaeta_small.xyz");

	// 3) Insert the points	
	pDt->insert(vInputPoints);

Delaunay benchmark Fade2D

Fade is very fast, one million points (two million triangles) take half a second on a recent desktop CPU


Visualization

Let’s show the surface triangulation. We can choose from various visualization methods:

  • Use the Postscript writer for a simple 2D triangulation
  • Write the 2.5D scene as Wavefront *.obj file which is supported by many viewers.
  • Interactive 3D scenes for web presentation can also be created
	
	// Write a wavefront *.obj file 
	pDt->writeObj("gaeta.obj");
	
	// Write a visualization for (technically up-to-date) web browsers
	pDt->writeWebScene("gaeta_webScene");

	// Write a file for the Geomview viewer
	pDt->showGeomview("gaeta_geomview.list");

	// 5) Write a postscript file that a 2D triangulation
	pDt->show("gaeta_2d_postscript.ps");

Here is for example the web scene output. You can rotate, zoom, shift and switch to different view modes.
Interactive 3D Viewer

ISO contour computation

	// 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);

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 above source code creates a global zone (consisting of all triangles), then it fetches the triangles from the zone and calls the isoContours() function defined below.

Function isoContours()

The isoContours(..) function below 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;
		}
	}
}

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.

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

Spread the word. Share this post!

2 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

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