2.5D Terrain Triangulation

Fade 2.5D triangulates xyz-points from terrains like lidar data, surface laser scans and other surface metrology techniques. The library is very fast, one million points (two million triangles) take half a second on a recent desktop CPU.

Fade2.5D and Fade2D are closely related but technically the two are independent libraries. When you program with Fade2.5D the class and function names are identically to the ones of Fade2D but the namespace is different, namely GEOM_FADE25D instead of GEOM_FADE2D.

Reading the terrain points

The example source code below is contained in the download. When using Visual Studio please run it on the command line, not in a VS debug window. The example starts with a simple ASCII file reader for xyz data: Three numbers per line, whitespace-separated.

void readTerrain(std::vector<Point2>& vInPoints,const char* filename)
{
	std::ifstream f(filename);
	while(!f.eof())
	{
		double x,y,z;
		f>>x>>y>>z;
		vInPoints.push_back(Point2(x,y,z));
	}
}

Triangulation of Surface Points

The example code below creates a new Fade triangulation and 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);

Visualization

Show what we have. We can still use the Postscript writer class Visualizer2, but it shows only 2D. To show a 3D scene Fade2.5D writes Wavefront *.obj files which are supported by many viewers. Moreover it can write a web scene and it writes files for the Geomview viewer.

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

Drag mouse to rotate, use SHIFT to
zoom and CTRL to translate

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