diff --git a/doc/gh_DFMergeAssemblies.rst b/doc/gh_DFMergeAssemblies.rst new file mode 100644 index 00000000..bf650361 --- /dev/null +++ b/doc/gh_DFMergeAssemblies.rst @@ -0,0 +1,8 @@ +.. image:: ../src/gh/components/DF_merge_assemblies/icon.png + :align: left + :width: 40px + +``DFMergeAssemblies`` component +=============================== + +.. ghcomponent_to_rst:: ../src/gh/components/DF_merge_assemblies \ No newline at end of file diff --git a/doc/gh_components.rst b/doc/gh_components.rst index 0dad3459..0b6781d1 100644 --- a/doc/gh_components.rst +++ b/doc/gh_components.rst @@ -81,6 +81,10 @@ DF has a Grasshopper_ plugin with a set of components that allows the user to in - .. image:: ../src/gh/components/DF_remove_statistical_outliers/icon.png - `gh_DFRemoveStatisticalOutliers `_ + * - .. image:: ../src/gh/components/DF_merge_assemblies/icon.png + - `gh_DFMergeAssemblies `_ + - + - .. toctree:: @@ -114,4 +118,5 @@ DF has a Grasshopper_ plugin with a set of components that allows the user to in gh_DFRemoveBeam gh_DFColorizeCloud gh_DFBrepToCloud - gh_DFRemoveStatisticalOutliers \ No newline at end of file + gh_DFRemoveStatisticalOutliers + gh_DFMergeAssemblies \ No newline at end of file diff --git a/src/diffCheck/geometry/DFMesh.cc b/src/diffCheck/geometry/DFMesh.cc index 71b8b16f..0669b56d 100644 --- a/src/diffCheck/geometry/DFMesh.cc +++ b/src/diffCheck/geometry/DFMesh.cc @@ -112,31 +112,67 @@ namespace diffCheck::geometry Eigen::Vector3d v1 = this->Vertices[triangle[1]]; Eigen::Vector3d v2 = this->Vertices[triangle[2]]; Eigen::Vector3d n = (v1 - v0).cross(v2 - v0); - double normOfNormal = n.norm(); n.normalize(); - Eigen::Vector3d projectedPoint = point - n * (n.dot(point - v0)) ; + // Project the point onto the plane of the triangle + Eigen::Vector3d projectedPoint = point - n * (n.dot(point - v0)); - double referenceTriangleArea = normOfNormal*0.5; - Eigen::Vector3d n1 = (v1 - v0).cross(projectedPoint - v0); - double area1 = n1.norm()*0.5; - Eigen::Vector3d n2 = (v2 - v1).cross(projectedPoint - v1); - double area2 = n2.norm()*0.5; - Eigen::Vector3d n3 = (v0 - v2).cross(projectedPoint - v2); - double area3 = n3.norm()*0.5; - double res = (area1 + area2 + area3 - referenceTriangleArea) / referenceTriangleArea; + // Compute vectors + Eigen::Vector3d v0v1 = v1 - v0; + Eigen::Vector3d v0v2 = v2 - v0; + Eigen::Vector3d v0p = projectedPoint - v0; - // 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) - double maxProjectionDistance = std::min({(v1 - v0).norm(), (v2 - v1).norm(), (v0 - v2).norm()}) / 2; + // Compute dot products + double dot00 = v0v2.dot(v0v2); + double dot01 = v0v2.dot(v0v1); + double dot02 = v0v2.dot(v0p); + double dot11 = v0v1.dot(v0v1); + double dot12 = v0v1.dot(v0p); - if (std::abs(res) < associationThreshold && (projectedPoint - point).norm() < maxProjectionDistance) + // Compute barycentric coordinates + double invDenom = 1.0 / (dot00 * dot11 - dot01 * dot01); + double u = (dot11 * dot02 - dot01 * dot12) * invDenom; + double v = (dot00 * dot12 - dot01 * dot02) * invDenom; + + // Check if point is in triangle + if ((u >= -associationThreshold) && (v >= -associationThreshold) && (u + v <= 1 + associationThreshold)) { - return true; + // Check if the point is close enough to the face + double maxProjectionDistance = std::min({(v1 - v0).norm(), (v2 - v1).norm(), (v0 - v2).norm()}) ; + if ((projectedPoint - point).norm() < maxProjectionDistance) + { + return true; + } } } return false; } + std::tuple DFMesh::ComputeOBBCenterAndAxis() + { + Eigen::Vector3d center = Eigen::Vector3d::Zero(); + Eigen::Vector3d axis = Eigen::Vector3d::Zero(); + + std::vector tightBoundingBox = this->GetTightBoundingBox(); + + Eigen::Vector3d deltaFirstDir = tightBoundingBox[1] - tightBoundingBox[0]; + Eigen::Vector3d deltaSecondDir = tightBoundingBox[2] - tightBoundingBox[0]; + Eigen::Vector3d deltaThirdDir = tightBoundingBox[3] - tightBoundingBox[0]; + + for (Eigen::Vector3d direction : {deltaFirstDir, deltaSecondDir, deltaThirdDir}) + { + if (direction.norm() > axis.norm()) + { + axis = direction; + } + } + axis.normalize(); + + center = tightBoundingBox[0] + deltaFirstDir/2 + deltaSecondDir/2 + deltaThirdDir/2; + + return std::make_tuple(center, axis); + } + void DFMesh::LoadFromPLY(const std::string &path) { std::shared_ptr tempMesh_ptr = diffCheck::io::ReadPLYMeshFromFile(path); diff --git a/src/diffCheck/geometry/DFMesh.hh b/src/diffCheck/geometry/DFMesh.hh index 2e865f34..19712411 100644 --- a/src/diffCheck/geometry/DFMesh.hh +++ b/src/diffCheck/geometry/DFMesh.hh @@ -86,6 +86,13 @@ namespace diffCheck::geometry */ bool IsPointOnFace(Eigen::Vector3d point, double associationThreshold = 0.1); + /** + * @brief Get the center and main axis of oriented boundung box of the mesh. It was developped for the cylinder case, but can be used for other shapes. + * + * @return std::tuple the first element is the center of the obb of the mesh, the second element is the main axis of the obb of the mesh + */ + std::tuple ComputeOBBCenterAndAxis(); + public: ///< I/O loader /** * @brief Read a mesh from a file diff --git a/src/diffCheck/segmentation/DFSegmentation.cc b/src/diffCheck/segmentation/DFSegmentation.cc index 37bb549a..3d09b8a6 100644 --- a/src/diffCheck/segmentation/DFSegmentation.cc +++ b/src/diffCheck/segmentation/DFSegmentation.cc @@ -96,12 +96,13 @@ namespace diffCheck::segmentation } std::vector> DFSegmentation::AssociateClustersToMeshes( + bool isCylinder, std::vector> referenceMesh, std::vector> &clusters, double angleThreshold, double associationThreshold) { - std::vector> jointFaces = std::vector>(); + std::vector> faceSegments = std::vector>(); // iterate through the mesh faces given as function argument if (referenceMesh.size() == 0) @@ -109,101 +110,223 @@ namespace diffCheck::segmentation DIFFCHECK_WARN("No mesh faces to associate with the clusters. Returning an empty point cloud."); return std::vector>(); } - for (std::shared_ptr face : referenceMesh) + + //differentiate between cylinder and other shapes + if (isCylinder) { - std::shared_ptr correspondingSegment; - std::shared_ptr facePoints = std::make_shared(); + for (std::shared_ptr face : referenceMesh) + { + std::tuple centerAndAxis = face->ComputeOBBCenterAndAxis(); + Eigen::Vector3d cylinderCenter = std::get<0>(centerAndAxis); + Eigen::Vector3d cylinderAxis = std::get<1>(centerAndAxis); - // Getting the center of the mesh face - Eigen::Vector3d faceCenter = face->Cvt2O3DTriangleMesh()->GetCenter(); + double faceDistance = std::numeric_limits::max(); + std::shared_ptr correspondingSegment; + std::shared_ptr facePoints = std::make_shared(); - // Getting the normal of the mesh face - Eigen::Vector3d faceNormal = face->GetFirstNormal(); - faceNormal.normalize(); - - double faceDistance = std::numeric_limits::max(); - if (clusters.size() == 0) - { - DIFFCHECK_WARN("No clusters to associate with the mesh faces. Returning an empty point cloud."); - return std::vector>(); - } - for (auto segment : clusters) - { - Eigen::Vector3d segmentCenter; - Eigen::Vector3d segmentNormal; + Eigen::Vector3d min = face->Vertices[0]; + Eigen::Vector3d max = face->Vertices[0]; + for (auto vertex : face->Vertices) + { + if(vertex.x() < min.x()){min.x() = vertex.x();} + if(vertex.y() < min.y()){min.y() = vertex.y();} + if(vertex.z() < min.z()){min.z() = vertex.z();} + if(vertex.x() > max.x()){max.x() = vertex.x();} + if(vertex.y() > max.y()){max.y() = vertex.y();} + if(vertex.z() > max.z()){max.z() = vertex.z();} + + } - for (auto point : segment->Points){segmentCenter += point;} - if (segment->GetNumPoints() > 0) + if (clusters.size() == 0) { - segmentCenter /= segment->GetNumPoints(); + DIFFCHECK_WARN("No clusters to associate with the mesh faces. Returning an empty point cloud."); + return std::vector>(); + } + for (auto segment : clusters) + { + Eigen::Vector3d segmentCenter; + Eigen::Vector3d segmentNormal; + + for (auto point : segment->Points) + { + segmentCenter += point; + } + if (segment->GetNumPoints() > 0) + { + segmentCenter /= segment->GetNumPoints(); + } + else + { + DIFFCHECK_WARN("Empty segment. Skipping the segment."); + continue; + } + + for (auto normal : segment->Normals) + { + segmentNormal += normal; + } + segmentNormal.normalize(); + + // we consider the distance to the cylinder axis, not the cylinder center + Eigen::Vector3d projectedSegmentCenter = (segmentCenter - cylinderCenter).dot(cylinderAxis) * cylinderAxis + cylinderCenter; + double currentDistance = (cylinderCenter - projectedSegmentCenter).norm(); + double absoluteDistance = (segmentCenter - cylinderCenter).norm(); + if (std::abs(cylinderAxis.dot(segmentNormal)) < angleThreshold && currentDistance < faceDistance && absoluteDistance < (max - min).norm()*associationThreshold) + { + correspondingSegment = segment; + faceDistance = currentDistance; + } } - else + if (correspondingSegment == nullptr) { - DIFFCHECK_WARN("Empty segment. Skipping the segment."); + DIFFCHECK_WARN("No segment found for the face. Skipping the face."); continue; } + bool hasColors = correspondingSegment->GetNumColors() > 0; + for (Eigen::Vector3d point : correspondingSegment->Points) + { + bool pointInFace = false; + if (face->IsPointOnFace(point, associationThreshold)) + { + facePoints->Points.push_back(point); + facePoints->Normals.push_back( + correspondingSegment->Normals[std::distance( + correspondingSegment->Points.begin(), + std::find(correspondingSegment->Points.begin(), + correspondingSegment->Points.end(), + point))] + ); + if (hasColors) + { + facePoints->Colors.push_back( + correspondingSegment->Colors[std::distance( + correspondingSegment->Points.begin(), + std::find(correspondingSegment->Points.begin(), + correspondingSegment->Points.end(), + point))] + ); + } + } + } + for(Eigen::Vector3d point : facePoints->Points) + { + correspondingSegment->Points.erase( + std::remove( + correspondingSegment->Points.begin(), + correspondingSegment->Points.end(), + point), + correspondingSegment->Points.end()); + } + if (correspondingSegment->GetNumPoints() == 0) + { + DIFFCHECK_WARN("No point was left in the segment. Deleting the segment."); + clusters.erase( + std::remove( + clusters.begin(), + clusters.end(), + correspondingSegment), + clusters.end()); + } + faceSegments.push_back(facePoints); + } + + } + else + { + for (std::shared_ptr face : referenceMesh) + { + std::shared_ptr correspondingSegment; + std::shared_ptr facePoints = std::make_shared(); - for (auto normal : segment->Normals){segmentNormal += normal;} - segmentNormal.normalize(); + // Getting the center of the mesh + Eigen::Vector3d faceCenter = face->Cvt2O3DTriangleMesh()->GetCenter(); - double currentDistance = (faceCenter - segmentCenter).norm(); - // if the distance is smaller than the previous one, update the distance and the corresponding segment - if (std::abs(sin(acos(faceNormal.dot(segmentNormal)))) < angleThreshold && currentDistance < faceDistance) + // Getting the normal of the mesh face + Eigen::Vector3d faceNormal = face->GetFirstNormal(); + faceNormal.normalize(); + + double faceDistance = std::numeric_limits::max(); + if (clusters.size() == 0) { - correspondingSegment = segment; - faceDistance = currentDistance; + DIFFCHECK_WARN("No clusters to associate with the mesh faces. Returning an empty point cloud."); + return std::vector>(); } - - } + for (auto segment : clusters) + { + Eigen::Vector3d segmentCenter; + Eigen::Vector3d segmentNormal; - if (correspondingSegment == nullptr) - { - DIFFCHECK_WARN("No segment found for the face. Skipping the face."); - continue; - } - bool hasColors = correspondingSegment->GetNumColors() > 0; + for (auto point : segment->Points){segmentCenter += point;} + if (segment->GetNumPoints() > 0) + { + segmentCenter /= segment->GetNumPoints(); + } + else + { + DIFFCHECK_WARN("Empty segment. Skipping the segment."); + continue; + } + for (auto normal : segment->Normals){segmentNormal += normal;} + segmentNormal.normalize(); + double currentDistance = (faceCenter - segmentCenter).norm(); + // if the distance is smaller than the previous one, update the distance and the corresponding segment + if (std::abs(sin(acos(faceNormal.dot(segmentNormal)))) < angleThreshold && currentDistance < faceDistance) + { + correspondingSegment = segment; + faceDistance = currentDistance; + } + } - for (Eigen::Vector3d point : correspondingSegment->Points) - { - bool pointInFace = false; - if (face->IsPointOnFace(point, associationThreshold)) + if (correspondingSegment == nullptr) { - facePoints->Points.push_back(point); - facePoints->Normals.push_back( - correspondingSegment->Normals[std::distance( - correspondingSegment->Points.begin(), - std::find(correspondingSegment->Points.begin(), - correspondingSegment->Points.end(), - point))] - ); - if (hasColors) + DIFFCHECK_WARN("No segment found for the face. Skipping the face."); + continue; + } + bool hasColors = correspondingSegment->GetNumColors() > 0; + + for (Eigen::Vector3d point : correspondingSegment->Points) + { + bool pointInFace = false; + if (face->IsPointOnFace(point, associationThreshold)) { - facePoints->Colors.push_back( - correspondingSegment->Colors[std::distance( + facePoints->Points.push_back(point); + facePoints->Normals.push_back( + correspondingSegment->Normals[std::distance( correspondingSegment->Points.begin(), std::find(correspondingSegment->Points.begin(), correspondingSegment->Points.end(), point))] ); + if (hasColors) + { + facePoints->Colors.push_back( + correspondingSegment->Colors[std::distance( + correspondingSegment->Points.begin(), + std::find(correspondingSegment->Points.begin(), + correspondingSegment->Points.end(), + point))] + ); + } } } + + for(Eigen::Vector3d point : facePoints->Points) + { + correspondingSegment->Points.erase( + std::remove( + correspondingSegment->Points.begin(), + correspondingSegment->Points.end(), + point), + correspondingSegment->Points.end()); + } + faceSegments.push_back(facePoints); } - - for(Eigen::Vector3d point : facePoints->Points) - { - correspondingSegment->Points.erase( - std::remove( - correspondingSegment->Points.begin(), - correspondingSegment->Points.end(), - point), - correspondingSegment->Points.end()); - } - jointFaces.push_back(facePoints); } - return jointFaces; + return faceSegments; } void DFSegmentation::CleanUnassociatedClusters( + bool isCylinder, std::vector> &unassociatedClusters, std::vector> &existingPointCloudSegments, std::vector>> meshes, @@ -215,128 +338,166 @@ namespace diffCheck::segmentation DIFFCHECK_WARN("No unassociated clusters. Nothing is done"); return; } - for (std::shared_ptr cluster : unassociatedClusters) + else { - std::shared_ptr correspondingMeshFace; - Eigen::Vector3d clusterCenter; - Eigen::Vector3d clusterNormal = Eigen::Vector3d::Zero(); - - if (cluster->GetNumPoints() == 0) - { - DIFFCHECK_WARN("Empty cluster. Skipping the cluster."); - continue; - } - for (Eigen::Vector3d point : cluster->Points) - { - clusterCenter += point; - } - clusterCenter /= cluster->GetNumPoints(); - if (cluster->GetNumNormals() == 0) - { - DIFFCHECK_WARN("Empty normals in the cluster. Skipping the cluster."); - continue; - } - for (Eigen::Vector3d normal : cluster->Normals) + for (std::shared_ptr cluster : unassociatedClusters) { - clusterNormal += normal; - } - clusterNormal.normalize(); + std::shared_ptr correspondingMeshFace; + Eigen::Vector3d clusterCenter; + Eigen::Vector3d clusterNormal = Eigen::Vector3d::Zero(); - int meshIndex = 0; - int faceIndex = 0 ; - int goodMeshIndex = 0; - int goodFaceIndex = 0; - double distance = std::numeric_limits::max(); - - if (meshes.size() == 0) - { - DIFFCHECK_WARN("No meshes to associate with the clusters. Skipping the cluster."); - continue; - } - for (std::vector> mesh : meshes) - { - if (mesh.size() == 0) + if (cluster->GetNumPoints() == 0) { - DIFFCHECK_WARN("Empty piece in the meshes vector. Skipping the mesh face vector."); + DIFFCHECK_WARN("Empty cluster. Skipping the cluster."); continue; } - for (std::shared_ptr meshFace : mesh) + if (cluster->GetNumNormals() == 0) { - Eigen::Vector3d faceCenter = Eigen::Vector3d::Zero(); - Eigen::Vector3d faceNormal = Eigen::Vector3d::Zero(); - - std::shared_ptr o3dFace = meshFace->Cvt2O3DTriangleMesh(); - - faceNormal = meshFace->GetFirstNormal(); - faceNormal.normalize(); - faceCenter = o3dFace->GetCenter(); - /* - To make sure we select the right meshFace, we add another metric: - Indeed, from experimentation, sometimes the wrong mesh face is selected, because it is parallel to the correct one - (so the normal don't play a role) and the center of the face is closer to the cluster center than the correct face. - To prevent this, we take into the account the angle between the line linking the center of the meshFace considered - and the center of the point cloud cluster and the normal of the cluster. This value should be close to pi/2 - - The following two lines are not super optimized but more readable than the optimized version - */ - - double dotProduct = clusterNormal.dot((clusterCenter - faceCenter).normalized()); - dotProduct = std::max(-1.0, std::min(1.0, dotProduct)); - double clusterNormalToJunctionLineAngle = std::acos(dotProduct); - - double currentDistance = (clusterCenter - faceCenter).norm() * std::abs(std::cos(clusterNormalToJunctionLineAngle)) - / std::min(std::abs(clusterNormal.dot(faceNormal)), 0.05) ; - if (std::abs(sin(acos(faceNormal.dot(clusterNormal)))) < angleThreshold && currentDistance < distance) + DIFFCHECK_WARN("Empty normals in the cluster. Skipping the cluster."); + continue; + } + if (meshes.size() == 0) + { + DIFFCHECK_WARN("No meshes to associate with the clusters. Skipping the cluster."); + continue; + } + for (Eigen::Vector3d point : cluster->Points) + { + clusterCenter += point; + } + clusterCenter /= cluster->GetNumPoints(); + for (Eigen::Vector3d normal : cluster->Normals) + { + clusterNormal += normal; + } + clusterNormal.normalize(); + + int meshIndex = 0; + int faceIndex = 0 ; + int goodMeshIndex = 0; + int goodFaceIndex = 0; + double distance = std::numeric_limits::max(); + + + for (std::vector> mesh : meshes) + { + if (mesh.size() == 0) + { + DIFFCHECK_WARN("Empty piece in the meshes vector. Skipping the mesh face vector."); + continue; + } + for (std::shared_ptr meshFace : mesh) { - goodMeshIndex = meshIndex; - goodFaceIndex = faceIndex; - distance = currentDistance; - correspondingMeshFace = meshFace; + Eigen::Vector3d faceCenter = Eigen::Vector3d::Zero(); + Eigen::Vector3d faceNormal = Eigen::Vector3d::Zero(); + if (isCylinder) + { + std::vector minmax = meshFace->GetTightBoundingBox(); + Eigen::Vector3d min = minmax[0]; + Eigen::Vector3d max = minmax[1]; + + std::tuple centerAndAxis = mesh[0]->ComputeOBBCenterAndAxis(); + Eigen::Vector3d center = std::get<0>(centerAndAxis); + Eigen::Vector3d axis = std::get<1>(centerAndAxis); + double dotProduct = clusterNormal.dot(axis); + dotProduct = std::max(-1.0, std::min(1.0, dotProduct)); + + double currentDistance = (center - clusterCenter).norm() ; + double adaptedDistance = currentDistance * std::abs(dotProduct); + + if (std::abs(dotProduct) < angleThreshold && adaptedDistance < distance && currentDistance < (max - min).norm()*associationThreshold) + { + goodMeshIndex = meshIndex; + goodFaceIndex = faceIndex; + distance = adaptedDistance; + correspondingMeshFace = meshFace; + } + + } + else + { + + + std::shared_ptr o3dFace = meshFace->Cvt2O3DTriangleMesh(); + + faceNormal = meshFace->GetFirstNormal(); + faceNormal.normalize(); + faceCenter = o3dFace->GetCenter(); + /* + To make sure we select the right meshFace, we add another metric: + Indeed, from experimentation, sometimes the wrong mesh face is selected, because it is parallel to the correct one + (so the normal don't play a role) and the center of the face is closer to the cluster center than the correct face. + To prevent this, we take into the account the angle between the line linking the center of the meshFace considered + and the center of the point cloud cluster and the normal of the cluster. This value should be close to pi/2 + + The following two lines are not super optimized but more readable than the optimized version + */ + + double dotProduct = clusterNormal.dot((clusterCenter - faceCenter).normalized()); + dotProduct = std::max(-1.0, std::min(1.0, dotProduct)); + double clusterNormalToJunctionLineAngle = std::acos(dotProduct); + + double currentDistance = (clusterCenter - faceCenter).norm() * std::abs(std::cos(clusterNormalToJunctionLineAngle)) + / std::min(std::abs(clusterNormal.dot(faceNormal)), 0.05) ; + if (std::abs(sin(acos(faceNormal.dot(clusterNormal)))) < angleThreshold && currentDistance < distance) + { + goodMeshIndex = meshIndex; + goodFaceIndex = faceIndex; + distance = currentDistance; + correspondingMeshFace = meshFace; + } + } + faceIndex++; } - faceIndex++; + meshIndex++; } - meshIndex++; - } - if (correspondingMeshFace == nullptr) - { - DIFFCHECK_WARN("No mesh face found for the cluster. Skipping the cluster."); - continue; - } - if (goodMeshIndex >= existingPointCloudSegments.size()) - { - DIFFCHECK_WARN("No segment found for the face. Skipping the face."); - continue; - } - std::shared_ptr completed_segment = existingPointCloudSegments[goodMeshIndex]; + if (correspondingMeshFace == nullptr) + { + DIFFCHECK_WARN("No mesh face found for the cluster. Skipping the cluster."); + continue; + } + if (goodMeshIndex >= existingPointCloudSegments.size()) + { + DIFFCHECK_WARN("No segment found for the face. Skipping the face."); + continue; + } + std::shared_ptr completed_segment = existingPointCloudSegments[goodMeshIndex]; - for (Eigen::Vector3d point : cluster->Points) - { - if (correspondingMeshFace->IsPointOnFace(point, associationThreshold)) + for (Eigen::Vector3d point : cluster->Points) { - completed_segment->Points.push_back(point); - completed_segment->Normals.push_back(cluster->Normals[std::distance(cluster->Points.begin(), std::find(cluster->Points.begin(), cluster->Points.end(), point))]); - if (cluster->GetNumColors() > 0) + if(isCylinder) { + completed_segment->Points.push_back(point); + completed_segment->Normals.push_back(cluster->Normals[std::distance(cluster->Points.begin(), std::find(cluster->Points.begin(), cluster->Points.end(), point))]); completed_segment->Colors.push_back(cluster->Colors[std::distance(cluster->Points.begin(), std::find(cluster->Points.begin(), cluster->Points.end(), point))]); } + else + if (correspondingMeshFace->IsPointOnFace(point, associationThreshold)) + { + completed_segment->Points.push_back(point); + completed_segment->Normals.push_back(cluster->Normals[std::distance(cluster->Points.begin(), std::find(cluster->Points.begin(), cluster->Points.end(), point))]); + completed_segment->Colors.push_back(cluster->Colors[std::distance(cluster->Points.begin(), std::find(cluster->Points.begin(), cluster->Points.end(), point))]); + } } - } - std::vector indicesToRemove; + std::vector indicesToRemove; - for (int i = 0; i < cluster->Points.size(); ++i) - { - if (std::find(completed_segment->Points.begin(), completed_segment->Points.end(), cluster->Points[i]) != completed_segment->Points.end()) + for (int i = 0; i < cluster->Points.size(); ++i) { - indicesToRemove.push_back(i); + if (std::find(completed_segment->Points.begin(), completed_segment->Points.end(), cluster->Points[i]) != completed_segment->Points.end()) + { + indicesToRemove.push_back(i); + } + } + for (auto it = indicesToRemove.rbegin(); it != indicesToRemove.rend(); ++it) + { + std::swap(cluster->Points[*it], cluster->Points.back()); + cluster->Points.pop_back(); + std::swap(cluster->Normals[*it], cluster->Normals.back()); + cluster->Normals.pop_back(); + std::swap(cluster->Colors[*it], cluster->Colors.back()); + cluster->Colors.pop_back(); } - } - for (auto it = indicesToRemove.rbegin(); it != indicesToRemove.rend(); ++it) - { - std::swap(cluster->Points[*it], cluster->Points.back()); - cluster->Points.pop_back(); - std::swap(cluster->Normals[*it], cluster->Normals.back()); - cluster->Normals.pop_back(); } } }; diff --git a/src/diffCheck/segmentation/DFSegmentation.hh b/src/diffCheck/segmentation/DFSegmentation.hh index d9a6d913..f88b3e71 100644 --- a/src/diffCheck/segmentation/DFSegmentation.hh +++ b/src/diffCheck/segmentation/DFSegmentation.hh @@ -27,6 +27,7 @@ namespace diffCheck::segmentation public: ///< segmentation refinement methods /** @brief Associates point cloud segments to mesh faces and merges them. It uses the center of mass of the segments and the mesh faces to find correspondances. For each mesh face it then iteratively associate the points of the segment that are actually on the mesh face. + * @param isCylinder a boolean to indicate if the model is a cylinder. If true, the method will use the GetCenterAndAxis method of the mesh to find the center and axis of the mesh. based on that, we only want points that have normals more or less perpendicular to the cylinder axis. * @param referenceMesh the vector of mesh faces to associate with the segments. It is a representation of a beam and its faces. * @param clusters the vector of clusters from cilantro to associate with the mesh faces of the reference mesh * @param angleThreshold the threshold to consider the a cluster as potential candidate for association. the value passed is the minimum sine of the angles. A value of 0 requires perfect alignment (angle = 0), while a value of 0.1 allows an angle of 5.7 degrees. @@ -34,12 +35,14 @@ namespace diffCheck::segmentation * @return std::shared_ptr The unified segments */ static std::vector> DFSegmentation::AssociateClustersToMeshes( + bool isCylinder, std::vector> referenceMesh, std::vector> &clusters, double angleThreshold = 0.1, double associationThreshold = 0.1); /** @brief Iterated through clusters and finds the corresponding mesh face. It then associates the points of the cluster that are on the mesh face to the segment already associated with the mesh face. + * @param isCylinder a boolean to indicate if the model is a cylinder. If true, the method will use the GetCenterAndAxis method of the mesh to find the center and axis of the mesh. based on that, we only want points that have normals more or less perpendicular to the cylinder axis. * @param unassociatedClusters the clusters from the normal-based segmentatinon that haven't been associated yet. * @param existingPointCloudSegments the already associated segments * @param meshes the mesh faces for all the model. This is used to associate the clusters to the mesh faces. @@ -47,6 +50,7 @@ namespace diffCheck::segmentation * @param associationThreshold the threshold to consider the points of a segment and a mesh face as associable. 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. */ static void DFSegmentation::CleanUnassociatedClusters( + bool isCylinder, std::vector> &unassociatedClusters, std::vector> &existingPointCloudSegments, std::vector>> meshes, diff --git a/src/diffCheckBindings.cc b/src/diffCheckBindings.cc index fff0697c..606e31d7 100644 --- a/src/diffCheckBindings.cc +++ b/src/diffCheckBindings.cc @@ -200,12 +200,14 @@ PYBIND11_MODULE(diffcheck_bindings, m) { py::arg("color_clusters") = false) .def_static("associate_clusters", &diffCheck::segmentation::DFSegmentation::AssociateClustersToMeshes, + py::arg("is_roundwood"), py::arg("reference_mesh"), py::arg("unassociated_clusters"), py::arg("angle_threshold") = 0.1, py::arg("association_threshold") = 0.1) .def_static("clean_unassociated_clusters", &diffCheck::segmentation::DFSegmentation::CleanUnassociatedClusters, + py::arg("is_roundwood"), py::arg("unassociated_clusters"), py::arg("associated_clusters"), py::arg("reference_mesh"), diff --git a/src/gh/components/DF_CAD_segmentator/code.py b/src/gh/components/DF_CAD_segmentator/code.py index f0cc9678..b2f67286 100644 --- a/src/gh/components/DF_CAD_segmentator/code.py +++ b/src/gh/components/DF_CAD_segmentator/code.py @@ -44,7 +44,9 @@ def RunScript(self, df_beams_meshes.append(df_b_mesh_faces) rh_beams_meshes.append(rh_b_mesh_faces) + # different association depending on the type of beam df_asssociated_cluster_faces = dfb_segmentation.DFSegmentation.associate_clusters( + is_cylinder=df_b.is_cylinder, reference_mesh=df_b_mesh_faces, unassociated_clusters=df_clouds, angle_threshold=i_angle_threshold, @@ -54,16 +56,24 @@ def RunScript(self, df_asssociated_cluster = dfb_geometry.DFPointCloud() for df_associated_face in df_asssociated_cluster_faces: df_asssociated_cluster.add_points(df_associated_face) - df_clusters.append(df_asssociated_cluster) - dfb_segmentation.DFSegmentation.clean_unassociated_clusters( + dfb_segmentation.DFSegmentation.clean_unassociated_clusters( + is_cylinder=df_b.is_cylinder, unassociated_clusters=df_clouds, - associated_clusters=df_clusters, - reference_mesh=df_beams_meshes, + associated_clusters=[df_asssociated_cluster], + reference_mesh=[df_b_mesh_faces], angle_threshold=i_angle_threshold, association_threshold=i_association_threshold ) + df_clusters.append(df_asssociated_cluster) + o_clusters = [df_cvt_bindings.cvt_dfcloud_2_rhcloud(cluster) for cluster in df_clusters] + for o_cluster in o_clusters: + if not o_cluster.IsValid: + o_cluster = None + ghenv.Component.AddRuntimeMessage(RML.Warning, "Some beams could not be segmented and were replaced by 'None'") # noqa: F821 + + return o_clusters diff --git a/src/gh/components/DF_build_assembly/code.py b/src/gh/components/DF_build_assembly/code.py index 7a9af9fc..0761d778 100644 --- a/src/gh/components/DF_build_assembly/code.py +++ b/src/gh/components/DF_build_assembly/code.py @@ -13,10 +13,14 @@ class DFBuildAssembly(component): def RunScript(self, i_assembly_name, - i_breps : System.Collections.Generic.IList[Rhino.Geometry.Brep]): + i_breps : System.Collections.Generic.IList[Rhino.Geometry.Brep], + i_is_roundwood : bool): beams: typing.List[DFBeam] = [] + if i_is_roundwood is None: + i_is_roundwood = False + for brep in i_breps: - beam = DFBeam.from_brep_face(brep) + beam = DFBeam.from_brep_face(brep, i_is_roundwood) beams.append(beam) o_assembly = DFAssembly(beams, i_assembly_name) diff --git a/src/gh/components/DF_build_assembly/metadata.json b/src/gh/components/DF_build_assembly/metadata.json index a0355271..18c7a446 100644 --- a/src/gh/components/DF_build_assembly/metadata.json +++ b/src/gh/components/DF_build_assembly/metadata.json @@ -36,6 +36,18 @@ "wireDisplay": "default", "sourceCount": 0, "typeHintID": "brep" + }, + { + "name": "i_is_roundwood", + "nickname": "i_is_roundwood", + "description": "Whether the beams are of roundwood type or not.", + "optional": true, + "allowTreeAccess": false, + "showTypeHints": true, + "scriptParamAccess": "item", + "wireDisplay": "default", + "sourceCount": 0, + "typeHintID": "bool" } ], "outputParameters": [ diff --git a/src/gh/components/DF_cloud_mesh_distance/code.py b/src/gh/components/DF_cloud_mesh_distance/code.py index c466a8ff..5bdd235d 100644 --- a/src/gh/components/DF_cloud_mesh_distance/code.py +++ b/src/gh/components/DF_cloud_mesh_distance/code.py @@ -35,10 +35,15 @@ def RunScript(self, return None, None, None, None, None, None # conversion - df_cloud_source_list = [df_cvt_bindings.cvt_rhcloud_2_dfcloud(i_cl_s) for i_cl_s in i_cloud_source] + siffed_df_cloud_source_list = [] + siffed_rh_mesh_target_list = [] + for i in range(len(i_cloud_source)): + if i_cloud_source[i] is not None: + siffed_df_cloud_source_list.append(df_cvt_bindings.cvt_rhcloud_2_dfcloud(i_cloud_source[i])) + siffed_rh_mesh_target_list.append(rh_mesh_target_list[i]) # calculate distances - o_result = df_error_estimation.df_cloud_2_rh_mesh_comparison(df_cloud_source_list, rh_mesh_target_list, i_signed_flag, i_swap) + o_result = df_error_estimation.df_cloud_2_rh_mesh_comparison(siffed_df_cloud_source_list, siffed_rh_mesh_target_list, i_signed_flag, i_swap) return o_result.distances, o_result.distances_rmse, o_result.distances_max_deviation, o_result.distances_min_deviation, o_result.distances_sd_deviation, o_result diff --git a/src/gh/components/DF_joint_segmentator/code.py b/src/gh/components/DF_joint_segmentator/code.py index 5446b278..0ee32aad 100644 --- a/src/gh/components/DF_joint_segmentator/code.py +++ b/src/gh/components/DF_joint_segmentator/code.py @@ -66,7 +66,7 @@ def RunScript(self, # find the corresponding clusters and merge them df_joint_cloud = diffcheck_bindings.dfb_geometry.DFPointCloud() - df_joint_face_segments = diffcheck_bindings.dfb_segmentation.DFSegmentation.associate_clusters(df_joint, df_cloud_clusters, i_angle_threshold, i_distance_threshold) + df_joint_face_segments = diffcheck_bindings.dfb_segmentation.DFSegmentation.associate_clusters(i_assembly.contains_cylinders, df_joint, df_cloud_clusters, i_angle_threshold, i_distance_threshold) for df_joint_face_segment in df_joint_face_segments: df_joint_cloud.add_points(df_joint_face_segment) diff --git a/src/gh/components/DF_merge_assemblies/code.py b/src/gh/components/DF_merge_assemblies/code.py new file mode 100644 index 00000000..49b4a1bd --- /dev/null +++ b/src/gh/components/DF_merge_assemblies/code.py @@ -0,0 +1,23 @@ +#! python3 + +import diffCheck +from diffCheck.df_geometries import DFBeam, DFAssembly + +from ghpythonlib.componentbase import executingcomponent as component + +import System + +class DFMergeAssemblies(component): + def RunScript(self, + i_new_name: str, + i_assemblies: System.Collections.Generic.IList[diffCheck.df_geometries.DFAssembly] + ) -> diffCheck.df_geometries.DFAssembly: + + beams = System.Collections.Generic.List[DFBeam]() + for assembly in i_assemblies: + for beam in assembly.beams: + beams.Add(beam) + + o_assembly = DFAssembly(beams, i_new_name) + + return o_assembly diff --git a/src/gh/components/DF_merge_assemblies/icon.png b/src/gh/components/DF_merge_assemblies/icon.png new file mode 100644 index 00000000..f42959c9 Binary files /dev/null and b/src/gh/components/DF_merge_assemblies/icon.png differ diff --git a/src/gh/components/DF_merge_assemblies/metadata.json b/src/gh/components/DF_merge_assemblies/metadata.json new file mode 100644 index 00000000..8f6f7be7 --- /dev/null +++ b/src/gh/components/DF_merge_assemblies/metadata.json @@ -0,0 +1,52 @@ +{ + "name": "DFMergeAssemblies", + "nickname": "DFMergeAssemblies", + "category": "diffCheck", + "subcategory": "Structure", + "description": "This component parse a series of DFAssemblies into a unique DFAssembly object.", + "exposure": 4, + "instanceGuid": "f6f0f6ca-0944-4311-828f-7308321b289f", + "ghpython": { + "hideOutput": true, + "hideInput": true, + "isAdvancedMode": true, + "marshalOutGuids": true, + "iconDisplay": 2, + "inputParameters": [ + { + "name": "i_new_name", + "nickname": "i_new_name", + "description": "The new name of the assembly to export.", + "optional": false, + "allowTreeAccess": true, + "showTypeHints": true, + "scriptParamAccess": "item", + "wireDisplay": "default", + "sourceCount": 0, + "typeHintID": "str" + }, + { + "name": "i_assemblies", + "nickname": "i_assemblies", + "description": "The assemblies to merge.", + "optional": true, + "allowTreeAccess": true, + "showTypeHints": true, + "scriptParamAccess": "list", + "wireDisplay": "default", + "sourceCount": 0, + "typeHintID": "ghdoc" + } + ], + "outputParameters": [ + { + "name": "o_assembly", + "nickname": "o_assembly", + "description": "The create DFAssembly object representing the timber elements.", + "optional": false, + "sourceCount": 0, + "graft": false + } + ] + } +} \ No newline at end of file diff --git a/src/gh/diffCheck/diffCheck/df_geometries.py b/src/gh/diffCheck/diffCheck/df_geometries.py index 75bfac29..f3fead03 100644 --- a/src/gh/diffCheck/diffCheck/df_geometries.py +++ b/src/gh/diffCheck/diffCheck/df_geometries.py @@ -93,7 +93,7 @@ def __post_init__(self): # if df_face is created from a rhino brep face, we store the rhino brep face self._rh_brepface = None - self.is_cylinder = False + self.is_roundwood = False def __repr__(self): return f"Face id: {(self.id)}, IsJoint: {self.is_joint} Loops: {len(self.all_loops)}" @@ -135,11 +135,6 @@ def from_brep_face(cls, all_loops = [] df_face: DFFace = cls([], joint_id) - if brep_face.IsCylinder(): - cls.is_cylinder = True - df_face._rh_brepface = brep_face - return df_face - for idx, loop in enumerate(brep_face.Loops): loop_curve = loop.To3dCurve() loop_curve = loop_curve.ToNurbsCurve() @@ -164,7 +159,7 @@ def to_brep_face(self): if self._rh_brepface is not None: return self._rh_brepface - if self.is_cylinder: + if self.is_roundwood: ghenv.Component.AddRuntimeMessage( # noqa: F821 RML.Warning, "The DFFace was a cylinder created from scratch \n \ , it cannot convert to brep.") @@ -285,6 +280,8 @@ class DFBeam: def __post_init__(self): self.name = self.name or "Unnamed Beam" self.faces = self.faces or [] + self.is_roundwood = False + self._joint_faces = [] self._side_faces = [] self._vertices = [] @@ -297,22 +294,25 @@ def __post_init__(self): self._index_assembly = None self._center = None + self.__id = uuid.uuid4().int def deepcopy(self): return DFBeam(self.name, [face.deepcopy() for face in self.faces]) + @classmethod - def from_brep_face(cls, brep): + def from_brep_face(cls, brep, is_roundwood=False): """ Create a DFBeam from a RhinoBrep object. It also removes duplicates and creates a list of unique faces. """ faces : typing.List[DFFace] = [] - data_faces = diffCheck.df_joint_detector.JointDetector(brep).run() + data_faces = diffCheck.df_joint_detector.JointDetector(brep, is_roundwood).run() for data in data_faces: face = DFFace.from_brep_face(data[0], data[1]) faces.append(face) beam = cls("Beam", faces) + beam.is_roundwood = is_roundwood return beam def to_brep(self): @@ -422,6 +422,13 @@ def __post_init__(self): self._all_vertices: typing.List[DFVertex] = [] self._all_joints: typing.List[DFJoint] = [] + for beam in self.beams: + if beam.is_roundwood: + self.contains_cylinders = True + break + else: + self.contains_cylinders = False + self._mass_center = None def __repr__(self): diff --git a/src/gh/diffCheck/diffCheck/df_joint_detector.py b/src/gh/diffCheck/diffCheck/df_joint_detector.py index ebc604d3..bef18712 100644 --- a/src/gh/diffCheck/diffCheck/df_joint_detector.py +++ b/src/gh/diffCheck/diffCheck/df_joint_detector.py @@ -16,6 +16,7 @@ class JointDetector: This class is responsible for detecting joints in a brep """ brep: Rhino.Geometry.Brep + is_roundwood: bool def __post_init__(self): self._faces = [] @@ -47,44 +48,132 @@ def _assign_ids(self, joint_face_ids): return extended_ids - def run(self): + def _find_largest_cylinder(self): + """ + Finds and returns the largest cylinder in the brep + + :return: the largest cylinder if beam detected as a cylinder, None otherwise """ - Run the joint detector. We use a dictionary to store the faces of the cuts based wethear they are cuts or holes. - - for cuts: If it is a cut we return the face, and the id of the joint the faces belongs to. - - for sides: If it is a face from the sides, we return the face and None. - :return: a list of faces from joins and faces + # extract all cylinders from the brep + candidate_open_cylinders = [] + for face in self.brep.Faces: + if not face.IsPlanar(): + candidate_open_cylinder = face.ToBrep() + candidate_open_cylinders.append(candidate_open_cylinder) + + # find largest cylinder + largest_cylinder = None + largest_srf = 0 + for candidate_cylinder in candidate_open_cylinders: + area = candidate_cylinder.GetArea() + if area > largest_srf and candidate_cylinder.CapPlanarHoles(sc.doc.ModelAbsoluteTolerance): + largest_srf = area + largest_cylinder = candidate_cylinder.CapPlanarHoles(sc.doc.ModelAbsoluteTolerance) + + # Check if the cylinder exists + if largest_cylinder is None: + print("No cylinder found") + return None + + return largest_cylinder + + def _find_joint_faces(self, bounding_geometry): + """ + Finds the brep faces that are joint faces. + + the structure of the disctionnary is as follows: + { + face_id: (face, is_inside) + ... + } + face_id is int + face is Rhino.Geometry.BrepFace + is_inside is bool + """ + faces = {} + if self.is_roundwood: + for idx, face in enumerate(self.brep.Faces): + faces[idx] = (face, face.IsPlanar(1000 * sc.doc.ModelAbsoluteTolerance)) + else: + for idx, face in enumerate(self.brep.Faces): + face_centroid = rg.AreaMassProperties.Compute(face).Centroid + coord = face.ClosestPoint(face_centroid) + projected_centroid = face.PointAt(coord[1], coord[2]) + faces[idx] = (face, + bounding_geometry.IsPointInside(projected_centroid, sc.doc.ModelAbsoluteTolerance, True) + * face.IsPlanar(1000 * sc.doc.ModelAbsoluteTolerance)) + + return faces + + def _compute_adjacency_of_faces(self, faces): """ - # brep vertices to cloud - df_cloud = diffCheck.diffcheck_bindings.dfb_geometry.DFPointCloud() - df_cloud.points = [np.array([vertex.Location.X, vertex.Location.Y, vertex.Location.Z]).reshape(3, 1) for vertex in self.brep.Vertices] - rh_OBB = diffCheck.df_cvt_bindings.cvt_dfOBB_2_rhbrep(df_cloud.get_tight_bounding_box()) - - # scale the box in the longest edge direction by 1.5 from center on both directions - rh_OBB_center = rh_OBB.GetBoundingBox(True).Center - edges = rh_OBB.Edges - edge_lengths = [edge.GetLength() for edge in edges] - longest_edge = edges[edge_lengths.index(max(edge_lengths))] - - rh_OBB_zaxis = rg.Vector3d(longest_edge.PointAt(1) - longest_edge.PointAt(0)) - rh_OBB_plane = rg.Plane(rh_OBB_center, rh_OBB_zaxis) - scale_factor = 0.09 - xform = rg.Transform.Scale( - rh_OBB_plane, - 1 - scale_factor, - 1 - scale_factor, - 1 + scale_factor - ) - rh_OBB.Transform(xform) - - # check if face's centers are inside the OBB - faces = {idx: (face, rh_OBB.IsPointInside(rg.AreaMassProperties.Compute(face).Centroid, sc.doc.ModelAbsoluteTolerance, True)) for idx, face in enumerate(self.brep.Faces)} - - # get the proximity faces of the joint faces - joint_face_ids = [[key] + [adj_face for adj_face in value[0].AdjacentFaces() if faces[adj_face][1] and adj_face != key] for key, value in faces.items() if value[1]] - - face_ids = self._assign_ids(joint_face_ids) - - self._faces = [(face, face_ids[idx]) for idx, face in enumerate(self.brep.Faces)] - - return self._faces + Computes the adjacency list of each face + + the structure of the dictionnary is as follows: + { + face_id: (face, [adj_face_id_1, adj_face_id_2, ...]) + ... + } + face_id is int + face is Rhino.Geometry.BrepFace + adj_face_id_1, adj_face_id_2, ... are int + """ + adjacency_of_faces = {} + for idx, face in faces.items(): + if not face[1]: + continue + adjacency_of_faces[idx] = (face[0], + [adj_face + for adj_face in face[0].AdjacentFaces() + if faces[adj_face][1] + and faces[adj_face][0].IsPlanar(1000 * sc.doc.ModelAbsoluteTolerance) + and adj_face != idx]) + + return adjacency_of_faces + + def run(self): + """ + Run the joint detector. We use a dictionary to store the faces of the cuts based wethear they are cuts or holes. + - for cuts: If it is a cut we return the face, and the id of the joint the faces belongs to. + - for sides: If it is a face from the sides, we return the face and None. + + :return: a list of faces from joins and faces + """ + + # brep vertices to cloud + df_cloud = diffCheck.diffcheck_bindings.dfb_geometry.DFPointCloud() + df_cloud.points = [np.array([vertex.Location.X, vertex.Location.Y, vertex.Location.Z]).reshape(3, 1) for vertex in self.brep.Vertices] + if self.is_roundwood: + bounding_geometry = self._find_largest_cylinder() + else: + bounding_geometry = diffCheck.df_cvt_bindings.cvt_dfOBB_2_rhbrep(df_cloud.get_tight_bounding_box()) + + # scale the bounding geometry in the longest edge direction by 1.5 from center on both directions + rh_Bounding_geometry_center = bounding_geometry.GetBoundingBox(True).Center + edges = bounding_geometry.Edges + edge_lengths = [edge.GetLength() for edge in edges] + longest_edge = edges[edge_lengths.index(max(edge_lengths))] + + rh_Bounding_geometry_zaxis = rg.Vector3d(longest_edge.PointAt(1) - longest_edge.PointAt(0)) + rh_Bounding_geometry_plane = rg.Plane(rh_Bounding_geometry_center, rh_Bounding_geometry_zaxis) + scale_factor = 0.1 + xform = rg.Transform.Scale( + rh_Bounding_geometry_plane, + 1 - scale_factor, + 1 - scale_factor, + 1 + scale_factor + ) + bounding_geometry.Transform(xform) + + faces = self._find_joint_faces(bounding_geometry) + adjacency_of_faces = self._compute_adjacency_of_faces(faces) + adjacency_of_faces = diffCheck.df_util.merge_shared_indexes(adjacency_of_faces) + joint_face_ids = [[key] + value[1] for key, value in adjacency_of_faces.items()] + + # get the proximity faces of the joint faces + face_ids = self._assign_ids(joint_face_ids) + + self._faces = [(face, face_ids[idx]) for idx, face in enumerate(self.brep.Faces)] + + return self._faces diff --git a/src/gh/diffCheck/diffCheck/df_util.py b/src/gh/diffCheck/diffCheck/df_util.py index 2a1b2717..b928984c 100644 --- a/src/gh/diffCheck/diffCheck/df_util.py +++ b/src/gh/diffCheck/diffCheck/df_util.py @@ -159,3 +159,44 @@ def get_doc_2_meters_unitf(): elif RhinoDoc.ModelUnitSystem == Rhino.UnitSystem.Yards: unit_scale = 0.9144 return unit_scale + +def merge_shared_indexes(original_dict): + """ + Merge the shared indexes of a dictionary + + Assume we have a dictionary with lists of indexes as values. + We want to merge the lists that share some indexes, in order to have a dictionary with, for each key, indexes that are not present under other keys. + + :param original_dict: the dictionary to merge + :return: the merged dictionary + """ + merged_dict = {} + index_to_key = {} + + for key, (face, indexes) in original_dict.items(): + merged_indexes = set(indexes) + keys_to_merge = set() + + for index in indexes: + if index in index_to_key: + keys_to_merge.add(index_to_key[index]) + + for merge_key in keys_to_merge: + merged_indexes.update(merged_dict[merge_key][1]) + # del merged_dict[merge_key] + + for index in merged_indexes: + index_to_key[index] = key + + merged_dict[key] = (face, list(merged_indexes)) + + keys_with_duplicates = {} + + for key in merged_dict.keys(): + for other_key, (face, indexes) in merged_dict.items(): + if key in indexes: + if key not in keys_with_duplicates: + keys_with_duplicates[key] = [] + keys_with_duplicates[key].append(other_key) + + return merged_dict