Implement pca patch/plane detec

src/diffCheck/geometry/DFPointCloud.cc

@@ -216,74 +216,21 @@ namespace diffCheck::geometry
             this->Normals.push_back(normal);
     }
 
-    std::vector<Eigen::Vector3d> DFPointCloud::GetPrincipalAxes(int nComponents)
+    Eigen::Vector3d DFPointCloud::FitPlaneRANSAC(
+            double distanceThreshold,
+            int ransacN,
+            int numIterations)
     {
-        std::vector<Eigen::Vector3d> principalAxes;
-
-        if (! this->HasNormals())
-        {
-            DIFFCHECK_WARN("The point cloud has no normals. Normals will be estimated with knn = 20.");
-            this->EstimateNormals(true, 20);
-        }
-
-        // Convert normals to Eigen matrix
-        Eigen::Matrix<double, 3, Eigen::Dynamic> normalMatrix(3, this->Normals.size());
-        for (size_t i = 0; i < this->Normals.size(); ++i)
-        {
-            normalMatrix.col(i) = this->Normals[i].cast<double>();
-        }
-
-        cilantro::KMeans<double, 3> kmeans(normalMatrix); 
-        kmeans.cluster(nComponents);
-
-        const cilantro::VectorSet3d& centroids = kmeans.getClusterCentroids();
-        const std::vector<size_t>& assignments = kmeans.getPointToClusterIndexMap();
-        std::vector<int> clusterSizes(nComponents, 0);
-        for (size_t i = 0; i < assignments.size(); ++i) 
-        {
-            clusterSizes[assignments[i]]++;
-        }
-        // Sort clusters by size
-        std::vector<std::pair<int, Eigen::Vector3d>> sortedClustersBySize(nComponents);
-        for (size_t i = 0; i < nComponents; ++i) 
+        if (this->Points.size() < ransacN)
         {
-            sortedClustersBySize[i] = {clusterSizes[i], centroids.col(i)};
+            DIFFCHECK_ERROR("Not enough points to fit a plane with RANSAC.");
+            return Eigen::Vector3d::Zero();
         }
-        std::sort(sortedClustersBySize.begin(), sortedClustersBySize.end(), [](const auto& a, const auto& b) 
-        {
-            return a.first > b.first;
-        });
-
-        for(size_t i = 0; i < nComponents; ++i) 
-        {
-            if(principalAxes.size() == 0)
-            {
-                principalAxes.push_back(sortedClustersBySize[i].second);
-            }
-            else
-            {
-                bool isAlreadyPresent = false;
-                for (const auto& axis : principalAxes)
-                {
-                    double dotProduct = std::abs(axis.dot(sortedClustersBySize[i].second));
-                    if (std::abs(dotProduct) > 0.7) // Threshold to consider as similar direction
-                    {
-                        isAlreadyPresent = true;
-                        break;
-                    }
-                }
-                if (!isAlreadyPresent)
-                {
-                    principalAxes.push_back(sortedClustersBySize[i].second);
-                }
-            }
-        }
-        if (principalAxes.size() < 2) // Fallback to OBB if k-means fails to provide enough distinct axes
-        {
-            open3d::geometry::OrientedBoundingBox obb = this->Cvt2O3DPointCloud()->GetOrientedBoundingBox();
-            principalAxes = {obb.R_.col(0), obb.R_.col(1), obb.R_.col(2)};
-        }
-        return principalAxes;
+
+        auto O3DPointCloud = this->Cvt2O3DPointCloud();
+        std::tuple< Eigen::Vector4d, std::vector<size_t>> planeModel = O3DPointCloud->SegmentPlane(distanceThreshold, ransacN, numIterations);
+        Eigen::Vector3d planeParameters = std::get<0>(planeModel).head<3>();
+        return planeParameters;
     }
 
     void DFPointCloud::Crop(const Eigen::Vector3d &minBound, const Eigen::Vector3d &maxBound)

src/diffCheck/geometry/DFPointCloud.hh

@@ -92,12 +92,17 @@ namespace diffCheck::geometry
         void RemoveStatisticalOutliers(int nbNeighbors, double stdRatio);
 
         /**
-         * @brief Get the nCompoments principal axes of the normals of the point cloud
-         * It is used to compute the pose of "boxy" point clouds. It relies on KMeans clustering to find the main axes of the point cloud.
-         *  @param nComponents the number of components to compute (default 6, each of 3 main axes in both directions)
-         * @return std::vector<Eigen::Vector3d> the principal axes of the point cloud ordered by number of normals
+         * @brief Fit a plane to the point cloud using RANSAC
+         * 
+         * @param distanceThreshold the distance threshold to consider a point as an inlier
+         * @param ransacN the number of points to sample for each RANSAC iteration
+         * @param numIterations the number of RANSAC iterations
+         * @return The Normal vector of the fitted plane as an Eigen::Vector3d
          */
-        std::vector<Eigen::Vector3d> GetPrincipalAxes(int nComponents = 6);
+        Eigen::Vector3d FitPlaneRANSAC(
+            double distanceThreshold = 0.01,
+            int ransacN = 3,
+            int numIterations = 100);
 
         /**
          *  @brief Crop the point cloud to a bounding box defined by the min and max bounds

src/diffCheck/segmentation/DFSegmentation.cc

@@ -330,7 +330,7 @@ namespace diffCheck::segmentation
     void DFSegmentation::CleanUnassociatedClusters(
         bool isCylinder,
         std::vector<std::shared_ptr<geometry::DFPointCloud>> &unassociatedClusters,
-        std::vector<std::shared_ptr<geometry::DFPointCloud>> &existingPointCloudSegments,
+        std::vector<std::vector<std::shared_ptr<geometry::DFPointCloud>>> &existingPointCloudSegments,
         std::vector<std::vector<std::shared_ptr<geometry::DFMesh>>> meshes,
         double angleThreshold,
         double associationThreshold)
@@ -459,12 +459,12 @@ namespace diffCheck::segmentation
                     DIFFCHECK_WARN("No mesh face found for the cluster. Skipping the cluster.");
                     continue;
                 }
-                if (goodMeshIndex >= existingPointCloudSegments.size())
+                if (goodMeshIndex >= existingPointCloudSegments.size() || goodFaceIndex >= existingPointCloudSegments[goodMeshIndex].size())
                 {
                     DIFFCHECK_WARN("No segment found for the face. Skipping the face.");
                     continue;
                 }
-                std::shared_ptr<geometry::DFPointCloud> completed_segment = existingPointCloudSegments[goodMeshIndex];
+                std::shared_ptr<geometry::DFPointCloud> completed_segment = existingPointCloudSegments[goodMeshIndex][goodFaceIndex];
 
                 for (Eigen::Vector3d point : cluster->Points)
                 {

src/diffCheck/segmentation/DFSegmentation.hh

@@ -44,15 +44,15 @@ namespace diffCheck::segmentation
         /** @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 existingPointCloudSegments the already associated segments per mesh face.
          * @param meshes the mesh faces for all the model. This is used to associate the clusters to the mesh faces.
          * * @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. 
          * @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<std::shared_ptr<geometry::DFPointCloud>> &unassociatedClusters,
-            std::vector<std::shared_ptr<geometry::DFPointCloud>> &existingPointCloudSegments,
+            std::vector<std::vector<std::shared_ptr<geometry::DFPointCloud>>> &existingPointCloudSegments,
             std::vector<std::vector<std::shared_ptr<geometry::DFMesh>>> meshes,
             double angleThreshold = 0.1,
             double associationThreshold = 0.1);

src/diffCheckBindings.cc

@@ -61,8 +61,11 @@ PYBIND11_MODULE(diffcheck_bindings, m) {
         .def("remove_statistical_outliers", &diffCheck::geometry::DFPointCloud::RemoveStatisticalOutliers, 
             py::arg("nb_neighbors"), py::arg("std_ratio"))
 
-        .def("get_principal_axes", &diffCheck::geometry::DFPointCloud::GetPrincipalAxes,
-            py::arg("n_components") = 6)
+        .def("fit_plane_ransac", &diffCheck::geometry::DFPointCloud::FitPlaneRANSAC,
+            py::arg("distance_threshold") = 0.01,
+            py::arg("ransac_n") = 3,
+            py::arg("num_iterations") = 100)
+
         .def("crop",
             (void (diffCheck::geometry::DFPointCloud::*)(const Eigen::Vector3d&, const Eigen::Vector3d&))
             &diffCheck::geometry::DFPointCloud::Crop,

src/gh/components/DF_CAD_segmentator/code.py

@@ -5,6 +5,7 @@
 import Rhino
 from ghpythonlib.componentbase import executingcomponent as component
 from Grasshopper.Kernel import GH_RuntimeMessageLevel as RML
+import ghpythonlib.treehelpers as th
 
 
 from diffCheck.diffcheck_bindings import dfb_segmentation
@@ -19,7 +20,7 @@ def RunScript(self,
         i_clouds: System.Collections.Generic.IList[Rhino.Geometry.PointCloud],
         i_assembly,
         i_angle_threshold: float = 0.1,
-        i_association_threshold: float = 0.1) -> Rhino.Geometry.PointCloud:
+        i_association_threshold: float = 0.1):
 
         if i_clouds is None or i_assembly is None:
             self.AddRuntimeMessage(RML.Warning, "Please provide a cloud and an assembly to segment.")
@@ -29,20 +30,18 @@ def RunScript(self,
         if i_association_threshold is None:
             i_association_threshold = 0.1
 
-        o_clusters = []
+        o_face_clusters = []
         df_clusters = []
         # we make a deepcopy of the input clouds
         df_clouds = [df_cvt_bindings.cvt_rhcloud_2_dfcloud(cloud.Duplicate()) for cloud in i_clouds]
 
         df_beams = i_assembly.beams
-        df_beams_meshes = []
-        rh_beams_meshes = []
 
         for df_b in df_beams:
+            o_face_clusters.append([])
+
             rh_b_mesh_faces = [df_b_f.to_mesh() for df_b_f in df_b.side_faces]
             df_b_mesh_faces = [df_cvt_bindings.cvt_rhmesh_2_dfmesh(rh_b_mesh_face) for rh_b_mesh_face in rh_b_mesh_faces]
-            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(
@@ -53,27 +52,30 @@ def RunScript(self,
                 association_threshold=i_association_threshold
             )
 
-            df_asssociated_cluster = dfb_geometry.DFPointCloud()
-            for df_associated_face in df_asssociated_cluster_faces:
-                df_asssociated_cluster.add_points(df_associated_face)
-
             dfb_segmentation.DFSegmentation.clean_unassociated_clusters(
                 is_roundwood=df_b.is_roundwood,
                 unassociated_clusters=df_clouds,
-                associated_clusters=[df_asssociated_cluster],
+                associated_clusters=[df_asssociated_cluster_faces],
                 reference_mesh=[df_b_mesh_faces],
                 angle_threshold=i_angle_threshold,
                 association_threshold=i_association_threshold
             )
 
+            o_face_clusters[-1] = [df_cvt_bindings.cvt_dfcloud_2_rhcloud(cluster) for cluster in df_asssociated_cluster_faces]
+
+            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)
 
-        o_clusters = [df_cvt_bindings.cvt_dfcloud_2_rhcloud(cluster) for cluster in df_clusters]
+        o_beam_clouds = [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
+        for i, o_beam_cloud in enumerate(o_beam_clouds):
+            if not o_beam_cloud.IsValid:
+                o_beam_clouds[i] = None
                 ghenv.Component.AddRuntimeMessage(RML.Warning, "Some beams could not be segmented and were replaced by 'None'")  # noqa: F821
 
+        o_face_clouds = th.list_to_tree(o_face_clusters)
 
-        return o_clusters
+        return [o_beam_clouds, o_face_clouds]

src/gh/components/DF_CAD_segmentator/metadata.json

@@ -64,9 +64,17 @@
         ],
         "outputParameters": [
             {
-                "name": "o_clusters",
-                "nickname": "o_clusters",
-                "description": "The clouds associated to each beam.",
+                "name": "o_beam_clouds",
+                "nickname": "o_beam_clouds",
+                "description": "The merged clouds associated to each beam.",
+                "optional": false,
+                "sourceCount": 0,
+                "graft": false
+            },
+            {
+                "name": "o_face_clouds",
+                "nickname": "o_face_clouds",
+                "description": "The datatree of clouds associated to each face.",
                 "optional": false,
                 "sourceCount": 0,
                 "graft": false

src/gh/components/DF_pose_estimation/code.py

@@ -1,24 +1,28 @@
+"""This compoment calculates the pose of a data tree of point clouds."""
 #! python3
 
 from diffCheck import df_cvt_bindings
 from diffCheck import df_poses
+from diffCheck.diffcheck_bindings import dfb_geometry
 
 import Rhino
 from Grasshopper.Kernel import GH_RuntimeMessageLevel as RML
+import Grasshopper
+import ghpythonlib.treehelpers as th
 
 from ghpythonlib.componentbase import executingcomponent as component
-import System
 
 
 class DFPoseEstimation(component):
     def RunScript(self,
-            i_clouds: System.Collections.Generic.List[Rhino.Geometry.PointCloud],
+            i_face_clouds: Grasshopper.DataTree[Rhino.Geometry.PointCloud],
             i_assembly,
-            i_save: bool,
-            i_reset: bool):
+            i_reset: bool,
+            i_save: bool):
 
+        clusters_per_beam = th.tree_to_list(i_face_clouds)
         # ensure assembly has enough beams
-        if len(i_assembly.beams) < len(i_clouds):
+        if len(i_assembly.beams) < len(clusters_per_beam):
             ghenv.Component.AddRuntimeMessage(RML.Warning, "Assembly has fewer beams than input clouds")  # noqa: F821
             return None, None
 
@@ -29,32 +33,34 @@ def RunScript(self,
             return None, None
 
         all_poses_this_time = []
-        for i, cloud in enumerate(i_clouds):
+        for i, face_clouds in enumerate(clusters_per_beam):
             try:
-                df_cloud = df_cvt_bindings.cvt_rhcloud_2_dfcloud(cloud)
-                if df_cloud is None:
-                    return None, None
-                if not df_cloud.has_normals():
-                    ghenv.Component.AddRuntimeMessage(RML.Error, f"Point cloud {i} has no normals. Please compute the normals.")  # noqa: F821
+                df_cloud = dfb_geometry.DFPointCloud()
 
-                df_points = df_cloud.get_axis_aligned_bounding_box()
-                df_point = (df_points[0] + df_points[1]) / 2
-                rh_point = Rhino.Geometry.Point3d(df_point[0], df_point[1], df_point[2])
+                rh_face_normals = []
+                for face_cloud in face_clouds:
+                    df_face_cloud = df_cvt_bindings.cvt_rhcloud_2_dfcloud(face_cloud)
+                    df_cloud.add_points(df_face_cloud)
+                    plane_normal = df_face_cloud.fit_plane_ransac()
+                    if all(plane_normal) == 0:
+                        ghenv.Component.AddRuntimeMessage(RML.Warning, f"There was a missing face in the cloud of beam {i}: the face was skipped in the pose estimation of that beam")  # noqa: F821
+                        continue
+                    rh_face_normals.append(Rhino.Geometry.Vector3d(plane_normal[0], plane_normal[1], plane_normal[2]))
 
-                axes = df_cloud.get_principal_axes(3)
-                vectors = []
-                for axe in axes:
-                    vectors.append(Rhino.Geometry.Vector3d(axe[0], axe[1], axe[2]))
+                df_bb_points = df_cloud.get_axis_aligned_bounding_box()
+                df_bb_centroid = (df_bb_points[0] + df_bb_points[1]) / 2
+                rh_bb_centroid = Rhino.Geometry.Point3d(df_bb_centroid[0], df_bb_centroid[1], df_bb_centroid[2])
 
-                new_xDirection, new_yDirection = df_poses.select_vectors(vectors, i_assembly.beams[i].plane.XAxis, i_assembly.beams[i].plane.YAxis)
+                new_xDirection, new_yDirection = df_poses.select_vectors(rh_face_normals, i_assembly.beams[i].plane.XAxis, i_assembly.beams[i].plane.YAxis)
 
                 pose = df_poses.DFPose(
-                    origin = [rh_point.X, rh_point.Y, rh_point.Z],
+                    origin = [rh_bb_centroid.X, rh_bb_centroid.Y, rh_bb_centroid.Z],
                     xDirection = [new_xDirection.X, new_xDirection.Y, new_xDirection.Z],
                     yDirection = [new_yDirection.X, new_yDirection.Y, new_yDirection.Z])
                 all_poses_this_time.append(pose)
-                plane = Rhino.Geometry.Plane(origin = rh_point, xDirection=new_xDirection, yDirection=new_yDirection)
+                plane = Rhino.Geometry.Plane(origin = rh_bb_centroid, xDirection=new_xDirection, yDirection=new_yDirection)
                 planes.append(plane)
+
             except Exception as e:
                 # Any unexpected error on this cloud, skip it and keep going
                 ghenv.Component.AddRuntimeMessage(RML.Error, f"Cloud {i}: processing failed ({e}); skipping.")  # noqa: F821

src/gh/components/DF_pose_estimation/metadata.json

@@ -14,9 +14,9 @@
         "iconDisplay": 2,
         "inputParameters": [
             {
-                "name": "i_clouds",
-                "nickname": "i_clouds",
-                "description": "clouds whose pose is to be calculated",
+                "name": "i_face_clouds",
+                "nickname": "i_face_clouds",
+                "description": "datatree of beam clouds whose pose is to be calculated",
                 "optional": false,
                 "allowTreeAccess": true,
                 "showTypeHints": true,

tests/unit_tests/DFPointCloudTest.cc

@@ -221,10 +221,12 @@ TEST_F(DFPointCloudTestFixture, Transform) {
 // Others
 //-------------------------------------------------------------------------
 
-TEST_F(DFPointCloudTestFixture, KMeansClusteringOfNormals) {
-    std::string path = diffCheck::io::GetTwoConnectedPlanesPlyPath();
-    diffCheck::geometry::DFPointCloud dfPointCloud2Planes;
-    dfPointCloud2Planes.LoadFromPLY(path);
-    std::vector<Eigen::Vector3d> axes = dfPointCloud2Planes.GetPrincipalAxes(2);
-    EXPECT_TRUE((axes[0] - Eigen::Vector3d(0, 0, 1)).norm() < 1e-2 || (axes[1] - Eigen::Vector3d(0, 0, 1)).norm() < 1e-2);
-}
+TEST_F(DFPointCloudTestFixture, FitPlaneRANSAC) {
+    std::shared_ptr<diffCheck::geometry::DFPointCloud> dfPointCloudPlane = std::make_shared<diffCheck::geometry::DFPointCloud>();
+    dfPointCloudPlane->LoadFromPLY(diffCheck::io::GetPlanePCWithOneOutliers());
+    Eigen::Vector3d planeNormal = dfPointCloudPlane->FitPlaneRANSAC(0.01, 3, 100);
+    // plane model should be close to (0, 0, 1, d)
+    EXPECT_NEAR(planeNormal[0], 0.0, 1e-2);
+    EXPECT_NEAR(planeNormal[1], 0.0, 1e-2);
+    EXPECT_NEAR(std::abs(planeNormal[2]), 1.0, 1e-2);
+}