Merging Triangulations and Point Clouds

The focus here is on merging: Two point clouds with different properties may overlap and need to be combined. Alternatively, you might need to merge triangulated objects with point clouds, such as a road with a terrain point cloud. This article presents practical solutions for merging triangulations and point clouds in 2.5D. The C++ source code explained in both this article and the next one, which covers merging two triangulations, is available in the download package under examples_real/mergingTriangulations.

First Demo: Merging Two Point Clouds

Scenario Overview:

Let’s assume that we have two point clouds: a coarse point cloud representing a larger area of terrain and a dense point cloud that covers a smaller, but more detailed region. Simply merging the two point clouds by copying can lead to issues: First, because the coarse cloud represents the overlapping area with less accuracy than then dense cloud. Second, in case of imperfect alignment of the two point clouds, simply merging them by copying would result in noise and steep artificial slopes.

Merging Point Clouds: Two overlapping point clouds of a terrain: One is coarse and represents a large terrain, the other one represents only a part of this terrain but is finer. — Two overlapping point clouds: The coarse point cloud (green) and the dense point cloud (red) used as inputs for this demo

A practical solution is to remove points from the coarse point cloud (green) that fall within the area covered by the dense point cloud (red) before combining them. However, if points from both clouds are too close to each other in the boundary region, this could result in steep artificial slopes. To address this issue, it is advisable to slightly enlarge the area designated for removal. Next, we will discuss the implementation of the first demo function, which demonstrates this approach.

Step 1: Loading Data from Files

We load the coarse and fine point clouds from .PLY files and visualize them as .VTK files using the VtkWriter. Files of this type can be inspected with ParaView and other 3D viewers.

// * Step 1 *   Load two input point clouds from files
cout<<"\nStep 1: Loading the input files"<<endl;
vector<Point2> vCoarseCloud,vFineCloud;
bool bCoarseOK= readPointsPLY("coarseCloud.ply",true,vCoarseCloud);
bool bFineOK= readPointsPLY("fineCloud.ply",true,vFineCloud);
if(!bCoarseOK || !bFineOK)
{
	cout<<"Input file(s) not found"<<std::endl;
	return;
}
cout << "- Coarse point cloud: " << vCoarseCloud.size() << endl;
cout << "- Fine point cloud: " << vFineCloud.size() << endl;
VtkWriter vtkw("d1_inputClouds.vtk"); // Visualize the two clouds and save them as vtk file for Paraview
vtkw.addPoints(vCoarseCloud, VTK_GREEN);
vtkw.addPoints(vFineCloud, VTK_RED);
vtkw.writeFile();

Step 2: Creating two Separate Triangulations

We create two separate triangulations for the coarse and fine point clouds, visualizing them as .VTK files. Additionally, we determine the median edge lengths of each triangulation, which we will need for the subsequent steps.

// * Step 2 *   Triangulate the two point clouds and visualize the resulting triangulations
cout<<"\nStep 2: Triangulating the point clouds"<<endl;
// For the coarse point cloud
Fade_2D coarseDt;
coarseDt.insert(vCoarseCloud);
coarseDt.showVtk("d1_coarseTriangulation.vtk",VTK_GREEN);
double coarseMedianEdgeLen2D(coarseDt.getMedianEdgeLength2D());
// For the fine point cloud
Fade_2D fineDt;
fineDt.insert(vFineCloud);
fineDt.showVtk("d1_fineTriangulation.vtk",VTK_RED);
double fineMedianEdgeLen2D(fineDt.getMedianEdgeLength2D());
cout<<"- Coarse: "<<coarseDt.numberOfTriangles()<<" triangles, median len: "<<coarseMedianEdgeLen2D<<endl;
cout<<"- Fine:   "<<fineDt.numberOfTriangles()<<" triangles, median len: "<<fineMedianEdgeLen2D<<endl;

Triangulation of a coarse point cloud — Triangulation of the coarse point cloud

Triangulation of a fine point cloud — Triangulation of the fine point cloud

Step 3: Computing the Shape of the Input Point Cloud

The area defined by a point cloud can be non-convex and may contain holes. In contrast, the Delaunay triangulation of any point cloud is always convex, and as a consequence the the triangulation of our fine point cloud includes triangles that do not represent real data, as shown in the image of the triangulation (red) above, where triangles with large edges do not correspond to data from the fine point cloud.
Therefore, we extract those triangles from the fine triangulation, fineDt, whose edges are sufficiently small because they truly represent the data in the point cloud. We use three times the median edge length of the fine triangulation as a threshold, classifying triangles with edges shorter than this threshold as valid data triangles, as displayed in the image below.

// * Step 3 *   Extract the 'data' triangles from fineDt: A Delaunay triangulation is always convex
// while the area that the input data represents is likely non-convex and may contain holes. We want
// to extract only the triangles that truly represent data.
cout<<"\nStep 3: Extracting the 'data' triangles from the fine triangulation"<<endl;
vector<Triangle2*> vAllFineT; // All triangles from fineDt
vector<Triangle2*> vDataFineT; // Extracted 'data' triangles from fineDt.
fineDt.getTrianglePointers(vAllFineT);
double lenThreshold(3.0 * fineMedianEdgeLen2D); // Edge length threshold for usable triangles
double sqThreshold(lenThreshold*lenThreshold); // Squared edge length threshold
for(Triangle2* pT:vAllFineT)
{
	if (pT->getSquaredEdgeLength2D(0) < sqThreshold &&
		pT->getSquaredEdgeLength2D(1) < sqThreshold &&
		pT->getSquaredEdgeLength2D(2) < sqThreshold)
	{
		vDataFineT.emplace_back(pT); // All three edges are small enough, store this triangle!
	}
}
cout<<"- Have extracted "<<vDataFineT.size()<<" usable triangles from the fine zone"<<endl;
if(vDataFineT.empty()) // Check!
{
	cout<<"vDataFineT.empty(), exit"<<endl;
	exit(1);
}

Step 4: Creating a Data Zone and its Offset Boundary

Filtered data triangles of the fine triangulation, surrounded by a red offset polygon. — The extracted data triangles of the fine triangulation, with an offset polygon outlined in red.

We create a zone from the valid triangles extracted in the previous step. Then, we compute an offset polygon around this zone, using the median edge length of the coarse triangulation as the offset distance. The purpose of this offset zone is to establish clearance, ensuring that points in the fine zone and points in the coarse zone are not too close together. This helps to avoid steep slopes that could arise from poor point quality or imperfect alignment between the two point clouds.

// * Step 4 *   Create a zone from the fine data triangles and extract its offset boundary.
// The offset is necessary to ensure that points from the two clouds are not too close
// together to avoid steep slopes.
cout<<"\nStep 4: Computing an offset boundary of the fine zone"<<endl;
Zone2* pFineZone(fineDt.createZone(vDataFineT));
vector<Segment2> vFineOffBoundary;
pFineZone->getOffsetBoundary(coarseMedianEdgeLen2D,vFineOffBoundary,1.0,15); // merge-angle=1.0, angle-step=15.0
cout<<"- Offset boundary: "<<vFineOffBoundary.size()<<" segments"<<endl;
VtkWriter vtkw2("d1_fineZone.vtk"); // Visualize the fine triangle area and the offset polygon
pFineZone->showVtk(&vtkw2,VTK_WHITE);
vtkw2.addSegments(vFineOffBoundary,VTK_RED);
vtkw2.writeFile();

Step 5: Using the Cookie Cutter

We use the Cookie-Cutter method to create a Zone2* pRemoveZone in the coarse triangulation that corresponds to the offset polygon from the fine triangulation. Then, we subtract this area from the coarse triangulation to obtain the remaining triangles in Zone2* pCoarseKeptZone. We then visualize the result as a VTK file.

Coarse triangulation with a section removed, which corresponds to the offset polygon of the fine triangulation. — The remaining coarse triangulation, `pCoarseKeptZone`, after the removal of the area corresponding to the offset polygon of the fine triangulation.

// * Step 5 *   Use the 'cookie cutter' method to create the area to be removed
// from the coarse triangulation, which is referred to as pCoarseRemoveZone. Then
// subtract this zone from the global coarse zone to obtain the remaining triangles,
// i.e., the area outside the offset boundary that should be preserved.
cout<<"\nStep 5: Determine the coarse area outside the fine offset boundary"<<endl;
Zone2* pCoarseRemoveZone(coarseDt.createZone_cookieCutter(vFineOffBoundary,true));
Zone2* pCoarseGlobalZone(coarseDt.createZone(NULL,ZL_GLOBAL));
Zone2* pCoarseKeptZone(zoneDifference(pCoarseGlobalZone,pCoarseRemoveZone));
if(pCoarseKeptZone==NULL)
{
	cout<<"No remaining triangles in the coarse zone, stop."<<endl;
	return;
}
cout<<"Keeping "<<pCoarseKeptZone->getNumberOfTriangles()<<" triangles"<<endl;
pCoarseKeptZone->showVtk("d1_coarseKeptZone.vtk",VTK_BROWN);

Step 6: Combining the Point Clouds

We extract the remaining coarse points from the Zone2* pCoarseKeptZone and insert them into the fine triangulation. Then we visualize the final result as a VTK file.

// * Step 6 *   Extract the points from pCoarseKeepZone and insert them into the
// fine triangulation (fineDt). Finally, visualize the resulting merged triangulation.
vector<Point2*> vCoarseZonePoints;
pCoarseKeptZone->getVertices(vCoarseZonePoints);
vector<Point2> vCoarseZonePoints2;
for(Point2* pVtx:vCoarseZonePoints) vCoarseZonePoints2.emplace_back(*pVtx);
fineDt.insert(vCoarseZonePoints2);
fineDt.showVtk("d1_mergedResult.vtk",VTK_GREEN);

The result is a high-quality triangulation of the merged point clouds that maintains the detailed accuracy of the dense point cloud while providing a smooth transition to the lower-sampled area.

Triangulation of merged point clouds showing a smooth transition between the areas, achieved by creating clearance with the offset polygon. — Merged point clouds: Despite imperfect alignment, the transition between the two areas is smooth due to the clearance established by the offset polygon.

Second Demo: Merging a Triangulated Object with a Point Cloud

Second Scenario:

In this second scenario, we merge a triangulated object, such as a road or another artificial structure, with a point cloud representing a terrain. This example differs from the first demo, where we worked with two point clouds, as we now have a point cloud and an existing triangulation. Consequently, we must ensure that the edges of the triangulated object do not get flipped during the insertion of points from the coarse point cloud, which requires us to protect its edges beforehand.

Step 1: Loading a Triangulation and a Point Cloud

We load a point cloud from a .PLY file and perform triangulation. Additionally, a predefined triangulation of the letter ‘F’ is loaded from a .PLY file. Both triangulations are then visualized and saved as .VTK files suitable for use in 3D viewers like ParaView and others.

vector<Point2> vCloud;
bool bCloudOK= readPointsPLY("coarseCloud.ply",true,vCloud);
Fade_2D objDt;
Zone2* pObjZone=objDt.importTrianglesFromPly("letterF.ply");
if(!bCloudOK || pObjZone==NULL)
{
	cout<<"Input file(s) not found"<<std::endl;
	return;
}
Fade_2D cloudDt;
cloudDt.insert(vCloud);
cout << "- Point cloud triangulation: " << cloudDt.numberOfTriangles() << " triangles"<<endl;
cout << "- Loaded object: " << objDt.numberOfTriangles() << " triangles"<<endl;
VtkWriter vtkw("d2_inputs.vtk"); // Visualize the two clouds and save them as vtk file for Paraview
cloudDt.showVtk(&vtkw,VTK_GREEN);
pObjZone->showVtk(&vtkw,VTK_RED);
vtkw.writeFile();

Overlapping triangulations of a point cloud (green) and a letter 'F' shape (red). — Two overlapping inputs: The green triangulation represents the input point cloud, while the red triangulation represents the letter “F” shape stored in `pObjZone`

Step 2: Computing an Offset Polygon

As in the previous example, we aim to maintain clearance between the triangulated object and the point cloud. To achieve this, we calculate the median edge length of the point cloud’s triangulation and use this value as the distance to compute an offset polygon around the triangulation of the letter ‘F’.

Triangulation of a letter 'F' shape with a red offset polygon. — Triangulation of the letter ‘F’ shape, loaded from a file, along with an offset polygon in red.

vector<Segment2> vOffBoundary;
pObjZone->getOffsetBoundary(cloudDt.getMedianEdgeLength2D(),vOffBoundary,1.0,15.0); // merge-angle=1.0, angle-step=15.0
cout<<"- Offset boundary: "<<vOffBoundary.size()<<" segments"<<endl;
VtkWriter vtkw2("d2_objZone_with_offset.vtk"); // Visualize the fine triangle area and the offset polygon
pObjZone->showVtk(&vtkw2,VTK_WHITE);
vtkw2.addSegments(vOffBoundary,VTK_RED);
vtkw2.writeFile();

Step 3: Computing the Remaining Zone

We use the Cookie Cutter method to create a Zone2* pRemoveZone in the triangulation of the point cloud that corresponds to the offset polygon of the letter ‘F’. Then, we compute the difference of the full triangulation and pRemoveZone to obtain the remaining area, pKeptZone. The points in this zone will be preserved.

Zone2* pRemoveZone(cloudDt.createZone_cookieCutter(vOffBoundary,true));
Zone2* pGlobalZone(cloudDt.createZone(NULL,ZL_GLOBAL));
Zone2* pKeptZone(zoneDifference(pGlobalZone,pRemoveZone));
if(pKeptZone==NULL)
{
	cout<<"No remaining triangles in the cloud zone, stop."<<endl;
	return;
}
cout<<"- Keeping "<<pKeptZone->getNumberOfTriangles()<<" cloud triangles"<<endl;
pKeptZone->showVtk("d2_keptZone.vtk",VTK_BROWN);

Triangulation of the loaded point cloud minus the area corresponding to the offset polygon. — The point cloud triangulation minus the area corresponding to the offset polygon.

Step 4: Protecting the Edges of the Object

We cannot simply insert the remaining points from the point cloud into the triangulation of the letter ‘F’, as this could cause the edges of the ‘F’ to flip. Therefore, we first need to protect its boundary edges. The code snippet below demonstrates how the boundary edges of the ‘F’ are extracted, converted to Segment2 objects, and subsequently inserted as constraint edges.

std::vector<Edge2> vExactBoundaryEdges;
pObjZone->getBoundaryEdges(vExactBoundaryEdges);
vector<Segment2> vExactBoundarySegments;
for(Edge2& e:vExactBoundaryEdges)
{
	vExactBoundarySegments.emplace_back(Segment2(*e.getSrc(),*e.getTrg()));
}
cout<<"- Protecting "<<vExactBoundaryEdges.size()<<" object edges"<<endl;
objDt.createConstraint(vExactBoundarySegments,CIS_CONSTRAINED_DELAUNAY);

Step 5: Integrating the Point Cloud into the Existing Triangulation

We extract the points from the remaining point cloud triangles stored in pKeptZone. Then, we insert these points into the existing triangulation of the letter ‘F’ and visualize the merged result as a .VTK file.

vector<Point2*> vRmgCloudVertices;
pKeptZone->getVertices(vRmgCloudVertices);
vector<Point2> vPoints;
for(Point2* pVtx:vRmgCloudVertices) vPoints.emplace_back(*pVtx);
cout<<"- Inserting the "<<vPoints.size()<<" remaining point cloud points"<<endl;
objDt.insert(vPoints);
objDt.showVtk("d2_mergedResult.vtk",VTK_GREEN);

Merging a triangulation and a point cloud: Resulting triangulation of the letter 'F' shape and the point cloud outside the offset zone. — Result: The triangulation of the letter ‘F’ and the point cloud outside the offset zone.

Need to merge two triangulated shapes? Be sure to check out the next article about merging two triangulations.