2.5D Terrain triangulation

2.5D Surface triangulation

Fade 2.5D triangulates surface elevation data, surface metrology data, laser scans and more. Fade2.5D has Fade2D as code base and is a 2.5D extension to Fade2D but technically the two are independent libraries. New Fade2.5D users should work through the 2D examples first because they also apply to Fade2.5D. In many cases the only difference is the z-coordinate. Additional features of Fade2.5D will be covered in the present article.

When you program with Fade2.5D the class names are identically to the ones of Fade2D but the namespace is different, namely GEOM_FADE25D instead of GEOM_FADE2D.

Triangulating terrain data

Please run this example in a terminal and not in a debug window of Visual Studio, otherwise it won’t find the path to the input file and exit with code 1.

File reader

The function below is a simple ASCII file reader for xyz point data: Three numbers per line, separated by whitespace.

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

The present example creates a new Fade triangulation, reads xyz data from an ASCII file 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);

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