Breakline Insertion, ISO Contours and the Cookie Cutter

Breakline insertion

In the field of land surveying constraint edges are maybe better known as “breaklines”. Similarly to the 2D case you can insert them into the triangle mesh, forcing segment subdivision or not. In addition you can insert a segment with its original height or you can adjust its height to the triangulation. This results in the four possibilities:

1. A breakline is inserted at its original elevation
2. It is inserted at its original elevation and subdivided
3. A breakline is projected to the surface and then it is inserted
4. Or it is projected to the surface, then it is subdivided and inserted (this is most likely what you want)

Breakline insertion at breakline elevation

The first two methods pose the simplest case. They simply insert breaklines at their own elevation. The code below:

1. creates a triangulation as described in the previous examples. Then it creates a vector of line segments that we will insert.
2. calls Fade_2D::createConstraint()to insert the segments as breaklines. Optionally makeDelaunay() can be called to subdivide the segments.

“Subdivision with the makeDelaunay() method is a technique that subdivides the line segment to be inserted until all sub-edges satisfy the empty cirlce property. This property is applied to the projection of the geometry to the xy-plane i.e., the height of the points is ignored. Thus makeDelaunay() is probably more useful in combination with breakline-insertion at surface elevation further below.”

“Be aware that any (x,y) coordinate pair can only have one height value z. Therefore, if it happens that an endpoint of your breakline exactly hits an already existing point, the height of the already existing point has higher priority.”

// 2: Insert breaklines (with subdivision, no drape)
void terrain_breaklines_withSubdiv_noDraping()
{
	// * 1 *   Create a surface and some segments
	Fade_2D dt;
	createTriangulation(dt);
    vector<Segment2> vSegments;
    createCircularPolygon(vSegments);

	// * 2 *   Insert the segments ('breaklines')
    ConstraintGraph2* pCG=dt.createConstraint(vSegments,CIS_CONSTRAINED_DELAUNAY);
	double minLen(0.1);
	pCG->makeDelaunay(minLen); // Subdivision

	// * 3 *   Visualize
	Visualizer3 vis("c2_terrain_withSubdiv_noDraping.list"); // For the Geomview viewer (Linux)
	dt.showGeomview(&vis,Visualizer3::CWHEAT);
	vis.writeSegments(vSegments,"0 0 1 0.5",true);
	dt.show("c2_terrain_withSubdiv_noDraping.ps"); // Postscript for gv, gsview or online ps-viewers
}

Breakline insertion at surface elevation

We want to insert breaklines (blue) at surface elevation (red). For this purpose we need one preparatory function call: We drape the segments (blue) onto the surface (red).

// 3: Insert breaklines (with subdivision, with drape)
void terrain_breaklines_withSubdiv_withDraping()
{
    ...
    // * 2 *   Insert the segments ('breaklines')
	double zTolerance(.1);
	vector<Segment2> vDrapedSegments;
	dt.drape(vSegments,vDrapedSegments,zTolerance);
	ConstraintGraph2* pCG=dt.createConstraint(vDrapedSegments,CIS_CONSTRAINED_DELAUNAY);
	double minLen(0.1);
	pCG->makeDelaunay(minLen); // Subdivision
    ...
}

The method Fade_2D::drape() takes the original segments and drapes them (subdivided if necessary) onto the TIN (Triangulated Irregular Network). This could lead to a large number of subdivisions. The zTolerance parameter is used to adjust the number of subdivisions to the user’s requirements. If a negative zTolerance is provided then the original segments are splitted at all intersections with triangulation edges. Otherwise they are only subdivided when required to maintain the maximum height difference zTolerance between the segment and the surface. The draped segments are then inserted into the TIN using Fade_2D::createConstraint() as before, see the source code below.

The Cookie Cutter

An often required function is the cutting of triangle meshes. You can do this with the above described breaklines, but there are special cases to consider and for convenience a separate function is available, with which you can easily cut out parts of a mesh like with a cookie cutter. Have a look at the below code snippet:

1. Create a triangulation. Create a circual polygon consisting of a few Segment2 objects
2. Create a Zone2* pGlobalZone consisting of all triangles.
3. Create a Zone2* pCookie using createZone_cookieCutter() with the before created polygon as parameter.
4. Through the boolean operation zoneDifference() create also a Zone2* pNotCookie consisting of the complementary set of triangles.

We have already discussed Boolean operations (4) on polygons i.e., union, intersection, difference and symmetric difference in this 2D example. The same applies analogously to 2.5D meshes.

// 4: Cookie Cutter
void cookieCutter()
{
	// * 1 *   Create a surface and a polygon
	Fade_2D dt;
	createTriangulation(dt);
    vector<Segment2> vPolygon;
    createCircularPolygon(vPolygon);

	// * 2 *   Show the triangulation and the polygon
	Visualizer3 vis("c4_surface.list");
	Zone2* pGlobalZone(dt.createZone(NULL,ZL_GLOBAL));
	pGlobalZone->showGeomview(&vis,Visualizer3::CWHEAT);
	vis.writeSegments(vPolygon,Visualizer3::CRED);

	// * 3 *   Cut out a piece with the cookieCutter
	Zone2* pCookie(dt.createZone_cookieCutter(vPolygon,true));
    pCookie->showGeomview("c4_cookie.list",Visualizer3::CGREEN);

	// * 4 *   Compute also the complementary piece
	Zone2* pNotCookie(zoneDifference(pGlobalZone,pCookie));
	pNotCookie->showGeomview("c4_notCookie.list",Visualizer3::CWHEAT);
}

ISO contour computation

Iso contours are lines consisting of the intersection of a terrain with a horizontal plane at a certain height i.e., at a certain z-value. Fade has an exceptionally fast algorithm to calculate iso lines. We use it like this:

1. Create a triangulation and fetch the triangles in vTriangles
2. Create IsoContours icoc(vTriangles)
3. Call isoc.getContours(), the parameters are the query-height as well as a vector of vector<Segment2> which is used to return the polylines.

// 5: Compute ISO lines (intersection with a horizontal plane)
void isoContours()
{
	// * 1 *   Create a surface and fetch the triangles
	Fade_2D dt;
	createTriangulation(dt);
	std::vector<Triangle2*> vTriangles;
	dt.getTrianglePointers(vTriangles);

    // * 2 *   Create IsoContours isoc
    IsoContours isoc(vTriangles);
    double minZ(isoc.getMinHeight());
    double maxZ(isoc.getMaxHeight());
    double height(minZ+0.45*(maxZ-minZ));

	// * 3 *   Compute the contours
	std::vector<std::vector<Segment2> > vvContours;
	isoc.getContours(height,vvContours,true,true);
	cout<<vvContours.size()<<"polylines at height "<<height<<endl;

    // Write the result
	for(size_t i=0;i<vvContours.size();++i)
	{
		std::vector<Segment2>& vContour(vvContours[i]);
		cout<<"\nContour no. "<<i<<" consists of "<<vContour.size()<<" segments"<<std::endl;
		vis.writeSegments(vContour,"1 0 0 0.5",true);
		for(size_t i=0;i<vContour.size();++i)
		{
			cout<<vContour[i]<<endl;
		}
    }
}

Height queries

Fade2.5D supports fast 2.5D height queries with arbitrary (x,y) coordinate pairs. The method bool getHeight(x,y,z,...) locates the triangle that contains (x,y) and finds the height value z of the point. If (x,y) is outside the triangulation then z is not set and the method returns false.

// 6: Find the height at arbitrary (x,y) coordinates
void heightQueries()
{
	// * 1 *   Create a surface and fetch the triangles
	Fade_2D dt;
	createTriangulation(dt);

	// * 2 *   Retrieve the height values at arbitrary (x,y) coordinates
	for(int i=0;i<10;++i)
	{
		double x((100.0*rand())/RAND_MAX);
		double y((100.0*rand())/RAND_MAX);
		double z;
		bool bHaveHeight(dt.getHeight(x,y,z));
		if(bHaveHeight)
		{
			cout<<"x="<<x<<", y="<<y<<", z="<<z<<endl;
		}
		else
		{
			cout<<"x="<<x<<", y="<<y<<" is out of bounds"<<endl;
		}
	}
}

* Height Queries

x=15.6679, y=40.0944, z=27.2008

x=12.979, y=10.8809, z=15.9557

x=99.8925, y=21.8257, z=17.4128

x=51.2932, y=83.9112, z=64.6294

x=61.264, y=29.6032, z=23.9884

x=63.7552, y=52.4287, z=35.941

x=49.3583, y=97.2775 is out of bounds

x=29.2517, y=77.1358, z=52.2147

x=52.6745, y=76.9914, z=55.0812

x=40.0229, y=89.1529 is out of bounds

“getHeight() has two optional parameters: Triangle2* pApproxT=NULL and double tolerance=0. pApproxT is used to speed up the point location in case that you know a nearby triangle. The tolerance parameter is used to calculate a height even if the query-point is slightly out of bounds, as it can happen due to rounding errors.”