Delaunay Terrain with Fade2.5D

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

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

Fade is very fast, one million points (two million triangles) take 0.16 seconds 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.

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