Fix/plane detection and joint detection

src/diffCheck/segmentation/DFSegmentation.cc

@@ -258,22 +258,21 @@ namespace diffCheck::segmentation
                     Eigen::Vector3d segmentCenter;
                     Eigen::Vector3d segmentNormal;
 
-                    for (auto point : segment->Points){segmentCenter += point;}
-                    if (segment->GetNumPoints() > 0)
-                    {
-                        segmentCenter /= segment->GetNumPoints();
-                    }
-                    else
+                    if (segment->GetNumPoints() == 0)
                     {
                         DIFFCHECK_WARN("Empty segment. Skipping the segment.");
                         continue;
                     }
+                    segmentCenter = segment->GetAxixAlignedBoundingBox()[0] + (segment->GetAxixAlignedBoundingBox()[1] - segment->GetAxixAlignedBoundingBox()[0])/2.0;
+
                     for (auto normal : segment->Normals){segmentNormal += normal;}
                     segmentNormal.normalize();
                     double currentDistance = (faceCenter - segmentCenter).norm();
                     double currentDitanceOrthogonalToFace = std::abs((faceCenter - segmentCenter).dot(faceNormal));
                     double currentAngle = std::abs(sin(acos(faceNormal.dot(faceCenter - segmentCenter))));
-                    if (std::abs(sin(acos(faceNormal.dot(segmentNormal)))) < angleThreshold && currentDitanceOrthogonalToFace < maximumFaceSegmentDistance  && currentDitanceOrthogonalToFace < faceDistance)
+                    if (std::abs(sin(acos(faceNormal.dot(segmentNormal)))) < angleThreshold 
+                        && currentDitanceOrthogonalToFace < maximumFaceSegmentDistance  
+                        && currentDitanceOrthogonalToFace < faceDistance)
                     {
                         correspondingSegment = segment;
                         faceDistance = currentDitanceOrthogonalToFace;
@@ -435,7 +434,9 @@ namespace diffCheck::segmentation
                             double currentDistance = (center - clusterCenter).norm() ;
                             double adaptedDistance = currentDistance * std::abs(dotProduct);
 
-                            if (std::abs(dotProduct) < angleThreshold && adaptedDistance < distance && currentDistance < (max - min).norm()*associationThreshold)
+                            if (std::abs(dotProduct) < angleThreshold 
+                                && adaptedDistance < distance 
+                                && currentDistance < (max - min).norm()*associationThreshold)
                             {
                                 goodMeshIndex = meshIndex;
                                 goodFaceIndex = faceIndex;

src/gh/components/DF_build_assembly/code.py

@@ -14,7 +14,8 @@ class DFBuildAssembly(component):
     def RunScript(self,
             i_assembly_name,
             i_breps : System.Collections.Generic.IList[Rhino.Geometry.Brep],
-            i_is_roundwood : bool):
+            i_is_roundwood : bool,
+            i_allow_curved_joint_faces : bool):
         beams: typing.List[DFBeam] = []
 
         if i_assembly_name is None or i_breps is None:
@@ -23,8 +24,12 @@ def RunScript(self,
         if i_is_roundwood is None:
             i_is_roundwood = False
 
+        if i_allow_curved_joint_faces is None:
+            i_allow_curved_joint_faces = False
+
         for brep in i_breps:
-            beam = DFBeam.from_brep_face(brep, i_is_roundwood)
+            brep.Faces.ShrinkFaces()
+            beam = DFBeam.from_brep_face(brep, i_is_roundwood, i_allow_curved_joint_faces)
             beams.append(beam)
 
         o_assembly = DFAssembly(beams, i_assembly_name)

src/gh/components/DF_build_assembly/metadata.json

@@ -48,6 +48,18 @@
                 "wireDisplay": "default",
                 "sourceCount": 0,
                 "typeHintID": "bool"
+            },
+            {
+                "name": "i_allow_curved_joint_faces",
+                "nickname": "i_allow_curved_joint_faces",
+                "description": "Whether to allow curved faces as joint faces.",
+                "optional": true,
+                "allowTreeAccess": false,
+                "showTypeHints": true,
+                "scriptParamAccess": "item",
+                "wireDisplay": "default",
+                "sourceCount": 0,
+                "typeHintID": "bool"
             }
         ],
         "outputParameters": [

src/gh/components/DF_pose_estimation/code.py

@@ -52,7 +52,7 @@ def RunScript(self,
                 rh_tentative_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(rh_face_normals, i_assembly.beams[i].plane.XAxis, i_assembly.beams[i].plane.YAxis)
-                rh_tentative_plane = Rhino.Geometry.Plane(rh_tentative_bb_centroid, new_yDirection, new_xDirection)
+                rh_tentative_plane = Rhino.Geometry.Plane(rh_tentative_bb_centroid, new_xDirection, new_yDirection)
 
                 rh_beam_cloud = Rhino.Geometry.PointCloud()
                 for face_cloud in face_clouds:

src/gh/diffCheck/diffCheck/df_geometries.py

@@ -178,8 +178,17 @@ def from_brep_face(cls,
             loop_vertices = loop_curve.Points
             loop = []
             for l_v in loop_vertices:
-                vertex = DFVertex(l_v.X, l_v.Y, l_v.Z)
-                loop.append(vertex)
+                rg_pt = rg.Point3d(l_v.X, l_v.Y, l_v.Z)
+                res = loop_curve.ClosestPoint(rg_pt)
+                if res:
+                    t = res[1]
+                else:
+                    t = 0 # this is a fallback, but it should not happen since the point is on the curve
+                point_on_curve = loop_curve.PointAt(t)
+                distance = rg.Point3d.DistanceTo(rg_pt, point_on_curve)
+                if distance < 10 * Rhino.RhinoDoc.ActiveDoc.ModelAbsoluteTolerance:
+                    vertex = DFVertex(l_v.X, l_v.Y, l_v.Z)
+                    loop.append(vertex)
             all_loops.append(loop)
 
         df_face = cls(all_loops, joint_id)
@@ -521,13 +530,18 @@ def compute_plane(self) -> rg.Plane:
 
         :return plane: The plane of the beam
         """
-        beam_direction = self.axis.Direction
+        bounding_geometry = diffCheck.df_util.compute_oriented_bounding_box(self.to_brep())
+        center = Rhino.Geometry.AreaMassProperties.Compute(bounding_geometry).Centroid
+        edge_lengths = [edge.GetLength() for edge in bounding_geometry.Edges]
+        longest_edge = bounding_geometry.Edges[edge_lengths.index(max(edge_lengths))]
+        z_axis = rg.Vector3d(longest_edge.PointAt(1) - longest_edge.PointAt(0))
+
         df_faces = [face for face in self.faces]
         sorted_df_faces = sorted(df_faces, key=lambda face: Rhino.Geometry.AreaMassProperties.Compute(face._rh_brepface).Area if face._rh_brepface else 0, reverse=True)
         largest_side_face_normal = sorted_df_faces[0].normal
         rh_largest_side_face_normal = rg.Vector3d(largest_side_face_normal[0], largest_side_face_normal[1], largest_side_face_normal[2])
 
-        return rg.Plane(self.center, rg.Vector3d.CrossProduct(beam_direction, rh_largest_side_face_normal), rh_largest_side_face_normal)
+        return rg.Plane(center, rg.Vector3d.CrossProduct(z_axis, rh_largest_side_face_normal), rh_largest_side_face_normal)
 
     def compute_joint_distances_to_midpoint(self) -> typing.List[float]:
         """
@@ -590,13 +604,13 @@ def compute_joint_angles(self) -> typing.List[float]:
         return jointface_angles
 
     @classmethod
-    def from_brep_face(cls, brep, is_roundwood=False):
+    def from_brep_face(cls, brep, is_roundwood=False, allow_curved_joint_faces=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, is_roundwood).run()
+        data_faces = diffCheck.df_joint_detector.JointDetector(brep, is_roundwood).run(allow_curved_joint_faces)
         for data in data_faces:
             face = DFFace.from_brep_face(data[0], data[1])
             faces.append(face)

src/gh/diffCheck/diffCheck/df_joint_detector.py

@@ -7,8 +7,6 @@
 import diffCheck.df_util
 import diffCheck.df_transformations
 
-import numpy as np
-
 
 @dataclass
 class JointDetector:
@@ -78,7 +76,7 @@ def _find_largest_cylinder(self):
 
         return largest_cylinder
 
-    def _find_joint_faces(self, bounding_geometry):
+    def _find_joint_faces(self, bounding_geometry, allow_curved_joint_faces=False):
         """
         Finds the brep faces that are joint faces.
 
@@ -100,7 +98,11 @@ def _find_joint_faces(self, bounding_geometry):
                 face_centroid = rg.AreaMassProperties.Compute(face).Centroid
                 coord = face.ClosestPoint(face_centroid)
                 projected_centroid = face.PointAt(coord[1], coord[2])
-                faces[idx] = (face,
+                if allow_curved_joint_faces:
+                    faces[idx] = (face,
+                                bounding_geometry.IsPointInside(projected_centroid, sc.doc.ModelAbsoluteTolerance, True))
+                else:
+                    faces[idx] = (face,
                             bounding_geometry.IsPointInside(projected_centroid, sc.doc.ModelAbsoluteTolerance, True)
                             * face.IsPlanar(1 * sc.doc.ModelAbsoluteTolerance))
 
@@ -132,32 +134,28 @@ def _compute_adjacency_of_faces(self, faces):
 
         return adjacency_of_faces
 
-    def run(self):
+    def run(self, allow_curved_joint_faces=False):
             """
                 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())
+                bounding_geometry = diffCheck.df_util.compute_oriented_bounding_box(self.brep)
 
             # 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
+            rh_Bounding_geometry_center = Rhino.Geometry.AreaMassProperties.Compute(bounding_geometry).Centroid
             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
+            scale_factor = 0.15
             xform = rg.Transform.Scale(
                 rh_Bounding_geometry_plane,
                 1 - scale_factor,
@@ -166,7 +164,7 @@ def run(self):
             )
             bounding_geometry.Transform(xform)
 
-            faces = self._find_joint_faces(bounding_geometry)
+            faces = self._find_joint_faces(bounding_geometry, allow_curved_joint_faces)
             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()]

src/gh/diffCheck/diffCheck/df_poses.py

@@ -134,23 +134,21 @@ def select_vectors(vectors, previous_xDirection, previous_yDirection):
         new_xDirection = sorted_vectors_by_alignment[0]
     else:
         new_xDirection = vectors[0]
-
-    condidates_for_yDirection = []
-    for v in vectors:
-        if compute_dot_product(v, new_xDirection) ** 2 < 0.5:
-            condidates_for_yDirection.append(v)
-
-    if not condidates_for_yDirection:
-        return new_xDirection, None
+    new_xDirection.Unitize()
 
     if previous_xDirection is not None and previous_yDirection is not None:
-        sorted_vectors_by_perpendicularity = sorted(condidates_for_yDirection, key=lambda v: abs(compute_dot_product(v, previous_yDirection)), reverse=True)
-        new_xDirection = sorted_vectors_by_alignment[0]
+        sorted_vectors_by_perpendicularity = sorted(vectors, key=lambda v: abs(compute_dot_product(v, previous_xDirection)))
         new_yDirection = sorted_vectors_by_perpendicularity[0] - compute_dot_product(sorted_vectors_by_perpendicularity[0], new_xDirection) * new_xDirection
+        if compute_dot_product(new_xDirection, previous_xDirection) < 0:
+            new_xDirection = -sorted_vectors_by_alignment[0]
+        if compute_dot_product(new_yDirection, previous_yDirection) < 0:
+            new_yDirection = -new_yDirection
         new_yDirection.Unitize()
     else:
-
-        sorted_vectors = sorted(vectors[1:], key=lambda v: compute_dot_product(v, new_xDirection)**2)
+        sorted_vectors = sorted(vectors[1:], key=lambda v: abs(compute_dot_product(v, new_xDirection)))
         new_yDirection = sorted_vectors[0] - compute_dot_product(sorted_vectors[0], new_xDirection) * new_xDirection
+        if previous_yDirection is not None and compute_dot_product(new_yDirection, previous_yDirection) < 0:
+            new_yDirection = -new_yDirection
         new_yDirection.Unitize()
+
     return new_xDirection, new_yDirection

src/gh/diffCheck/diffCheck/df_util.py

@@ -2,6 +2,10 @@
 import Rhino.Geometry as rg
 import scriptcontext as sc
 
+import diffCheck.diffcheck_bindings
+import diffCheck.df_cvt_bindings
+import numpy as np
+
 import typing
 
 
@@ -180,3 +184,24 @@ def merge_shared_indexes(original_dict):
         if not intersection_found:
             new_dict[key] = (face, indexes)
     return new_dict
+
+def compute_oriented_bounding_box(brep):
+    """
+    Computes the oriented bounding box of a brep.
+    We use the point cloud of the vertices of the brep's 4 largest faces to compute the bounding box.
+
+    :param brep: the brep to compute the bounding box of
+    :return: the oriented bounding box of the brep
+    """
+    df_cloud = diffCheck.diffcheck_bindings.dfb_geometry.DFPointCloud()
+    sorted_faces = sorted(brep.Faces, key=lambda f : Rhino.Geometry.AreaMassProperties.Compute(f.ToBrep()).Area, reverse = True)
+    if len(sorted_faces) > 4:
+        largest_faces = sorted_faces[:4]
+    else:
+        largest_faces = sorted_faces
+    bb_vertices = []
+    for face in largest_faces:
+        for v in face.ToBrep().Vertices:
+            bb_vertices.append(Rhino.Geometry.Point3d(v.Location.X, v.Location.Y, v.Location.Z))
+    df_cloud.points = [np.array([vertex.X, vertex.Y, vertex.Z]).reshape(3, 1) for vertex in bb_vertices]
+    return diffCheck.df_cvt_bindings.cvt_dfOBB_2_rhbrep(df_cloud.get_tight_bounding_box())