Skip to content

Commit 8cca66a

Browse files
DamienGilliardDamienGilliard
committed
UPDATE-FIX: segmentation methods more robust, and IsPointOnFace method now method of DFMesh
1 parent e8f7f0c commit 8cca66a

File tree

4 files changed

+147
-98
lines changed

4 files changed

+147
-98
lines changed

‎src/diffCheck/geometry/DFMesh.cc‎

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -95,6 +95,46 @@ namespace diffCheck::geometry
9595
return this->NormalsFace[0];
9696
}
9797

98+
bool DFMesh::IsPointOnFace(Eigen::Vector3d point, double associationThreshold)
99+
{
100+
/*
101+
To check if the point is in the face, we take into account all the triangles forming the face.
102+
We calculate the area of each triangle, then check if the sum of the areas of the tree triangles
103+
formed by two of the points of the referencr triangle and our point is equal to the reference triangle area
104+
(within a user-defined margin). If it is the case, the triangle is in the face.
105+
*/
106+
std::vector<Eigen::Vector3i> faceTriangles = this->Faces;
107+
for (Eigen::Vector3i triangle : faceTriangles)
108+
{
109+
Eigen::Vector3d v0 = this->Vertices[triangle[0]];
110+
Eigen::Vector3d v1 = this->Vertices[triangle[1]];
111+
Eigen::Vector3d v2 = this->Vertices[triangle[2]];
112+
Eigen::Vector3d n = (v1 - v0).cross(v2 - v0);
113+
double normOfNormal = n.norm();
114+
n.normalize();
115+
116+
Eigen::Vector3d projectedPoint = point - n * (n.dot(point - v0)) ;
117+
118+
double referenceTriangleArea = normOfNormal*0.5;
119+
Eigen::Vector3d n1 = (v1 - v0).cross(projectedPoint - v0);
120+
double area1 = n1.norm()*0.5;
121+
Eigen::Vector3d n2 = (v2 - v1).cross(projectedPoint - v1);
122+
double area2 = n2.norm()*0.5;
123+
Eigen::Vector3d n3 = (v0 - v2).cross(projectedPoint - v2);
124+
double area3 = n3.norm()*0.5;
125+
double res = (area1 + area2 + area3 - referenceTriangleArea) / referenceTriangleArea;
126+
127+
// arbitrary value to avoid false positives (points that, when projected on the triangle, are in it, but that are actually located too far from the mesh to actually belong to it)
128+
double maxProjectionDistance = std::min({(v1 - v0).norm(), (v2 - v1).norm(), (v0 - v2).norm()});
129+
130+
if (std::abs(res) < associationThreshold && (projectedPoint - point).norm() < maxProjectionDistance)
131+
{
132+
return true;
133+
}
134+
}
135+
return false;
136+
}
137+
98138
void DFMesh::LoadFromPLY(const std::string &path)
99139
{
100140
std::shared_ptr<diffCheck::geometry::DFMesh> tempMesh_ptr = diffCheck::io::ReadPLYMeshFromFile(path);

‎src/diffCheck/geometry/DFMesh.hh‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -78,6 +78,14 @@ namespace diffCheck::geometry
7878
*/
7979
Eigen::Vector3d GetFirstNormal();
8080

81+
/**
82+
* @brief Check if a point is on a face of the mesh
83+
*
84+
* @param point the point to check
85+
* @param associationThreshold the threshold to consider the point associable to the mesh. It is the ratio between the surface of the closest mesh triangle of the mesh face, and the sum of the areas of the three triangles described by the point projected on the mesh face and two of the mesh triangle vertices. The lower the number, the more strict the association will be and some poinnts on the mesh face might be wrongfully excluded. In theory, in a perfect case, a value of 0 could be used, but in practice, values of 0.05-0.2 are more realistic, depending on the application.
86+
*/
87+
bool IsPointOnFace(Eigen::Vector3d point, double associationThreshold = 0.1);
88+
8189
public: ///< I/O loader
8290
/**
8391
* @brief Read a mesh from a file

‎src/diffCheck/segmentation/DFSegmentation.cc‎

Lines changed: 99 additions & 86 deletions
Original file line numberDiff line numberDiff line change
@@ -106,8 +106,8 @@ namespace diffCheck::segmentation
106106
// iterate through the mesh faces given as function argument
107107
if (referenceMesh.size() == 0)
108108
{
109-
DIFFCHECK_WARN("No mesh faces to associate with the clusters. Returning the first clusters untouched.");
110-
return clusters[0];
109+
DIFFCHECK_WARN("No mesh faces to associate with the clusters. Returning an empty point cloud.");
110+
return unifiedPointCloud;
111111
}
112112
for (std::shared_ptr<diffCheck::geometry::DFMesh> face : referenceMesh)
113113
{
@@ -123,8 +123,8 @@ namespace diffCheck::segmentation
123123
double faceDistance = std::numeric_limits<double>::max();
124124
if (clusters.size() == 0)
125125
{
126-
DIFFCHECK_WARN("No clusters to associate with the mesh faces. Returning the first mesh face untouched.");
127-
return clusters[0];
126+
DIFFCHECK_WARN("No clusters to associate with the mesh faces. Returning an empty point cloud.");
127+
return unifiedPointCloud;
128128
}
129129
for (auto segment : clusters)
130130
{
@@ -163,7 +163,7 @@ namespace diffCheck::segmentation
163163
for (Eigen::Vector3d point : correspondingSegment->Points)
164164
{
165165
bool pointInFace = false;
166-
if (DFSegmentation::IsPointOnFace(face, point, associationThreshold))
166+
if (face->IsPointOnFace(point, associationThreshold))
167167
{
168168
unifiedPointCloud->Points.push_back(point);
169169
unifiedPointCloud->Normals.push_back(
@@ -190,16 +190,6 @@ namespace diffCheck::segmentation
190190
point),
191191
correspondingSegment->Points.end());
192192
}
193-
// removing the corresponding segment if it is empty after the point transfer
194-
// if (correspondingSegment->GetNumPoints() == 0)
195-
// {
196-
// clusters.erase(
197-
// std::remove(
198-
// clusters.begin(),
199-
// clusters.end(),
200-
// correspondingSegment),
201-
// clusters.end());
202-
// }
203193
}
204194
return unifiedPointCloud;
205195
}
@@ -218,8 +208,9 @@ namespace diffCheck::segmentation
218208
}
219209
for (std::shared_ptr<geometry::DFPointCloud> cluster : unassociatedClusters)
220210
{
211+
std::shared_ptr<geometry::DFMesh> correspondingMeshFace;
221212
Eigen::Vector3d clusterCenter;
222-
Eigen::Vector3d clusterNormal;
213+
Eigen::Vector3d clusterNormal = Eigen::Vector3d::Zero();
223214

224215
if (cluster->GetNumPoints() == 0)
225216
{
@@ -231,33 +222,40 @@ namespace diffCheck::segmentation
231222
clusterCenter += point;
232223
}
233224
clusterCenter /= cluster->GetNumPoints();
225+
if (cluster->GetNumNormals() == 0)
226+
{
227+
DIFFCHECK_WARN("Empty normals in the cluster. Skipping the cluster.");
228+
continue;
229+
}
234230
for (Eigen::Vector3d normal : cluster->Normals)
235231
{
236232
clusterNormal += normal;
237233
}
238234
clusterNormal.normalize();
239235

240-
int meshIndex = std::numeric_limits<int>::max();
241-
int faceIndex = std::numeric_limits<int>::max() ;
236+
int meshIndex = 0;
237+
int faceIndex = 0 ;
238+
int goodMeshIndex = 0;
239+
int goodFaceIndex = 0;
242240
double distance = std::numeric_limits<double>::max();
243241

244242
if (meshes.size() == 0)
245243
{
246244
DIFFCHECK_WARN("No meshes to associate with the clusters. Skipping the cluster.");
247245
continue;
248246
}
249-
for (std::vector<std::shared_ptr<geometry::DFMesh>> piece : meshes)
247+
for (std::vector<std::shared_ptr<geometry::DFMesh>> mesh : meshes)
250248
{
251-
if (piece.size() == 0)
249+
if (mesh.size() == 0)
252250
{
253251
DIFFCHECK_WARN("Empty piece in the meshes vector. Skipping the mesh face vector.");
254252
continue;
255253
}
256-
for (std::shared_ptr<geometry::DFMesh> meshFace : piece)
254+
for (std::shared_ptr<geometry::DFMesh> meshFace : mesh)
257255
{
258-
Eigen::Vector3d faceCenter;
259-
Eigen::Vector3d faceNormal;
260-
256+
Eigen::Vector3d faceCenter = Eigen::Vector3d::Zero();
257+
Eigen::Vector3d faceNormal = Eigen::Vector3d::Zero();
258+
261259
std::shared_ptr<open3d::geometry::TriangleMesh> o3dFace = meshFace->Cvt2O3DTriangleMesh();
262260

263261
faceNormal = meshFace->GetFirstNormal();
@@ -273,88 +271,103 @@ namespace diffCheck::segmentation
273271
The following two lines are not super optimized but more readable than the optimized version
274272
*/
275273

276-
double clusterNormalToJunctionLineAngle = std::abs(std::acos(clusterNormal.dot((clusterCenter - faceCenter).normalized())));
274+
double dotProduct = clusterNormal.dot((clusterCenter - faceCenter).normalized());
275+
dotProduct = std::max(-1.0, std::min(1.0, dotProduct));
276+
double clusterNormalToJunctionLineAngle = std::acos(dotProduct);
277277

278-
double currentDistance = (clusterCenter - faceCenter).norm() * std::abs(std::cos(clusterNormalToJunctionLineAngle)) / std::abs(clusterNormal.dot(faceNormal));
278+
double currentDistance = (clusterCenter - faceCenter).norm() * std::abs(std::cos(clusterNormalToJunctionLineAngle))
279+
/ std::min(std::abs(clusterNormal.dot(faceNormal)), 0.05) ;
279280
if (std::abs(sin(acos(faceNormal.dot(clusterNormal)))) < angleThreshold && currentDistance < distance)
280281
{
282+
goodMeshIndex = meshIndex;
283+
goodFaceIndex = faceIndex;
281284
distance = currentDistance;
282-
meshIndex = std::distance(meshes.begin(), std::find(meshes.begin(), meshes.end(), piece));
283-
faceIndex = std::distance(piece.begin(), std::find(piece.begin(), piece.end(), meshFace));
285+
correspondingMeshFace = meshFace;
284286
}
287+
faceIndex++;
285288
}
289+
meshIndex++;
286290
}
287-
if (meshIndex >= meshes.size() || faceIndex >= meshes[meshIndex].size())
291+
292+
if (correspondingMeshFace == nullptr)
288293
{
289-
// this one generates a lot of warnings
290294
DIFFCHECK_WARN("No mesh face found for the cluster. Skipping the cluster.");
291295
continue;
292296
}
293-
std::shared_ptr<geometry::DFPointCloud> completed_segment = existingPointCloudSegments[meshIndex];
297+
if (goodMeshIndex >= existingPointCloudSegments.size())
298+
{
299+
DIFFCHECK_WARN("No segment found for the face. Skipping the face.");
300+
continue;
301+
}
302+
std::shared_ptr<geometry::DFPointCloud> completed_segment = existingPointCloudSegments[goodMeshIndex];
303+
294304
for (Eigen::Vector3d point : cluster->Points)
295305
{
296-
std::vector<Eigen::Vector3i> faceTriangles = meshes[meshIndex][faceIndex]->Faces;
297-
if (IsPointOnFace(meshes[meshIndex][faceIndex], point, associationThreshold))
306+
if (correspondingMeshFace->IsPointOnFace(point, associationThreshold))
298307
{
299308
completed_segment->Points.push_back(point);
300309
completed_segment->Normals.push_back(cluster->Normals[std::distance(cluster->Points.begin(), std::find(cluster->Points.begin(), cluster->Points.end(), point))]);
301310
}
302311
}
303-
// std::vector<int> indicesToRemove;
304-
305-
// for (int i = 0; i < cluster->Points.size(); ++i)
306-
// {
307-
// if (std::find(completed_segment->Points.begin(), completed_segment->Points.end(), cluster->Points[i]) != completed_segment->Points.end())
308-
// {
309-
// indicesToRemove.push_back(i);
310-
// }
311-
// }
312-
// for (auto it = indicesToRemove.rbegin(); it != indicesToRemove.rend(); ++it)
313-
// {
314-
// std::swap(cluster->Points[*it], cluster->Points.back());
315-
// cluster->Points.pop_back();
316-
// std::swap(cluster->Normals[*it], cluster->Normals.back());
317-
// cluster->Normals.pop_back();
318-
// }
319-
}
320-
};
312+
std::vector<int> indicesToRemove;
321313

322-
bool DFSegmentation::IsPointOnFace(
323-
std::shared_ptr<diffCheck::geometry::DFMesh> face,
324-
Eigen::Vector3d point,
325-
double associationThreshold)
326-
{
327-
/*
328-
To check if the point is in the face, we take into account all the triangles forming the face.
329-
We calculate the area of each triangle, then check if the sum of the areas of the tree triangles
330-
formed by two of the points of the referencr triangle and our point is equal to the reference triangle area
331-
(within a user-defined margin). If it is the case, the triangle is in the face.
332-
*/
333-
std::vector<Eigen::Vector3i> faceTriangles = face->Faces;
334-
for (Eigen::Vector3i triangle : faceTriangles)
335-
{
336-
Eigen::Vector3d v0 = face->Vertices[triangle[0]];
337-
Eigen::Vector3d v1 = face->Vertices[triangle[1]];
338-
Eigen::Vector3d v2 = face->Vertices[triangle[2]];
339-
Eigen::Vector3d n = (v1 - v0).cross(v2 - v0);
340-
double normOfNormal = n.norm();
341-
n.normalize();
342-
343-
Eigen::Vector3d projectedPoint = point - n * (n.dot(point - v0)) ;
344-
345-
double referenceTriangleArea = normOfNormal*0.5;
346-
Eigen::Vector3d n1 = (v1 - v0).cross(projectedPoint - v0);
347-
double area1 = n1.norm()*0.5;
348-
Eigen::Vector3d n2 = (v2 - v1).cross(projectedPoint - v1);
349-
double area2 = n2.norm()*0.5;
350-
Eigen::Vector3d n3 = (v0 - v2).cross(projectedPoint - v2);
351-
double area3 = n3.norm()*0.5;
352-
double res = (area1 + area2 + area3 - referenceTriangleArea) / referenceTriangleArea;
353-
if (std::abs(res) < associationThreshold)
314+
for (int i = 0; i < cluster->Points.size(); ++i)
354315
{
355-
return true;
316+
if (std::find(completed_segment->Points.begin(), completed_segment->Points.end(), cluster->Points[i]) != completed_segment->Points.end())
317+
{
318+
indicesToRemove.push_back(i);
319+
}
320+
}
321+
for (auto it = indicesToRemove.rbegin(); it != indicesToRemove.rend(); ++it)
322+
{
323+
std::swap(cluster->Points[*it], cluster->Points.back());
324+
cluster->Points.pop_back();
325+
std::swap(cluster->Normals[*it], cluster->Normals.back());
326+
cluster->Normals.pop_back();
356327
}
357328
}
358-
return false;
359-
}
329+
};
330+
331+
// bool DFSegmentation::IsPointOnFace(
332+
// std::shared_ptr<diffCheck::geometry::DFMesh> face,
333+
// Eigen::Vector3d point,
334+
// double associationThreshold)
335+
// {
336+
// /*
337+
// To check if the point is in the face, we take into account all the triangles forming the face.
338+
// We calculate the area of each triangle, then check if the sum of the areas of the tree triangles
339+
// formed by two of the points of the referencr triangle and our point is equal to the reference triangle area
340+
// (within a user-defined margin). If it is the case, the triangle is in the face.
341+
// */
342+
// std::vector<Eigen::Vector3i> faceTriangles = face->Faces;
343+
// for (Eigen::Vector3i triangle : faceTriangles)
344+
// {
345+
// Eigen::Vector3d v0 = face->Vertices[triangle[0]];
346+
// Eigen::Vector3d v1 = face->Vertices[triangle[1]];
347+
// Eigen::Vector3d v2 = face->Vertices[triangle[2]];
348+
// Eigen::Vector3d n = (v1 - v0).cross(v2 - v0);
349+
// double normOfNormal = n.norm();
350+
// n.normalize();
351+
352+
// Eigen::Vector3d projectedPoint = point - n * (n.dot(point - v0)) ;
353+
354+
// double referenceTriangleArea = normOfNormal*0.5;
355+
// Eigen::Vector3d n1 = (v1 - v0).cross(projectedPoint - v0);
356+
// double area1 = n1.norm()*0.5;
357+
// Eigen::Vector3d n2 = (v2 - v1).cross(projectedPoint - v1);
358+
// double area2 = n2.norm()*0.5;
359+
// Eigen::Vector3d n3 = (v0 - v2).cross(projectedPoint - v2);
360+
// double area3 = n3.norm()*0.5;
361+
// double res = (area1 + area2 + area3 - referenceTriangleArea) / referenceTriangleArea;
362+
363+
// // arbitrary value to avoid false positives (points that, when projected on the triangle, are in it, but that are actually located too far from the mesh to actually belong to it)
364+
// double maxProjectionDistance = std::min({(v1 - v0).norm(), (v2 - v1).norm(), (v0 - v2).norm()});
365+
366+
// if (std::abs(res) < associationThreshold && (projectedPoint - point).norm() < maxProjectionDistance)
367+
// {
368+
// return true;
369+
// }
370+
// }
371+
// return false;
372+
// }
360373
} // namespace diffCheck::segmentation

‎src/diffCheck/segmentation/DFSegmentation.hh‎

Lines changed: 0 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -52,17 +52,5 @@ namespace diffCheck::segmentation
5252
std::vector<std::vector<std::shared_ptr<geometry::DFMesh>>> meshes,
5353
double angleThreshold = 0.1,
5454
double associationThreshold = 0.1);
55-
56-
// FIXME: this one should go into DFMesh
57-
private: ///< helper methods
58-
/** @brief private method to check if a point is on a face of a triangle mesh triangle, within a certain association threshold. This takes into account the fact that, in 3D, a point can be "above" a triangle of a triangle mesh but still considered as being on the mesh face.
59-
* @param face the triangle mesh face to check the point against
60-
* @param point the point to check
61-
* @param associationThreshold the threshold to consider the point associable to the mesh. It is the ratio between the surface of the closest mesh triangle and the sum of the areas of the three triangles that form the rest of the pyramid described by the mesh triangle and the point we want to associate or not. The lower the number, the more strict the association will be and some poinnts on the mesh face might be wrongfully excluded.
62-
*/
63-
static bool DFSegmentation::IsPointOnFace(
64-
std::shared_ptr<diffCheck::geometry::DFMesh> face,
65-
Eigen::Vector3d point,
66-
double associationThreshold);
6755
};
6856
} // namespace diffCheck::segmentation

0 commit comments

Comments
 (0)