In practice, we often encounter situations where two point clouds with different properties overlap and need to be combined, or where a triangulated object, such as a road, must be integrated into a terrain. These tasks require careful approaches for best possible results.
This article discusses two practical solutions that you can adapt to your project. You can find the complete C++ source code of this example in the download package under the path 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 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.
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;
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, meaning that the triangulation of the fine point cloud can also include triangles that do not represent its 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
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.
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.
Second Demo: Merging a Triangulated Object with a Point Cloud
Scenario Overview:
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();
pObjZoneStep 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’.
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);
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 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);
