Categories

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

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
{
// * 1 *   Create a surface and a polygon
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

// * 4 *   Compute also the complementary piece
}
```

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