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

The source code described in this post is contained in the download, see examples_25D/terrain.cpp. 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.

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

You do not need to create handcrafted reader/writer functions like above. Have a look at the module File I/O where various file-reader and -writer functions for points and segments are documented.

Triangulate the Terrain Points

Let’s call the above described readTerrain() function to read the surface points. Then create a new Fade2.5D triangulation and insert the points.

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