diff --git a/CMakeLists.txt b/CMakeLists.txt index 79f25866..a68c648f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -10,6 +10,9 @@ include(options) # disabling warnings if (MSVC) add_compile_options(/wd4003) + # Disable _SECURE_SCL to avoid issues with stdext::checked_array_iterator + # in newer MSVC versions with Open3D's bundled fmt library + add_compile_definitions(_SECURE_SCL=0 _HAS_STDEXT_CHECKED_ARRAY_ITERATOR=0) endif() # check that the -DCMAKE_BUILD_TYPE is set diff --git a/cmake/copy_dlls_script.cmake b/cmake/copy_dlls_script.cmake new file mode 100644 index 00000000..15123be0 --- /dev/null +++ b/cmake/copy_dlls_script.cmake @@ -0,0 +1,18 @@ +# This script is executed at build time to copy DLLs +# It uses variables passed from the main CMake configuration + +# Remove old DLLs from the destination directory +file(GLOB old_dlls "${DIR_TO_CLEAN}/*.dll") +if(old_dlls) + file(REMOVE ${old_dlls}) +endif() + +# Get all DLLs from the source directory (evaluated at build time) +file(GLOB dll_files "${SRC_DIR}/*.dll") + +# Copy each DLL to the destination directory +foreach(dll_file ${dll_files}) + get_filename_component(dll_name ${dll_file} NAME) + message(STATUS "Copying ${dll_name} to ${DST_DIR}") + file(COPY ${dll_file} DESTINATION ${DST_DIR}) +endforeach() diff --git a/cmake/external_tools.cmake b/cmake/external_tools.cmake index 3b0c5714..f9ff2cfe 100644 --- a/cmake/external_tools.cmake +++ b/cmake/external_tools.cmake @@ -232,19 +232,19 @@ endfunction() # ------------------------------------------------------------------------------ function (copy_dlls directory_to_copy_dlls post_build_target) - message (STATUS "Erasing old DLLs and copy new ones to ${directory_to_copy_dlls}") - file(GLOB files ${directory_to_copy_dlls}/*.dll) - foreach(file ${files}) - message(STATUS "Removing ${file}") - file(REMOVE ${file}) - endforeach() - file(GLOB files ${CMAKE_BINARY_DIR}/bin/${CMAKE_BUILD_TYPE}/*.dll) - foreach(file ${files}) - message(STATUS "Copying ${file} to ${directory_to_copy_dlls}") - add_custom_command(TARGET ${post_build_target} POST_BUILD - COMMAND ${CMAKE_COMMAND} -E copy - ${file} - ${directory_to_copy_dlls} - ) - endforeach() + message (STATUS "Configuring DLL copy to ${directory_to_copy_dlls} at build time") + + # Get the path to the script relative to the project source dir + set(COPY_DLLS_SCRIPT ${PROJECT_SOURCE_DIR}/cmake/copy_dlls_script.cmake) + + # Add a post-build command that will copy DLLs at build time + # This ensures the file list is evaluated at build time, not configure time + add_custom_command(TARGET ${post_build_target} POST_BUILD + COMMAND ${CMAKE_COMMAND} + -DDIR_TO_CLEAN="${directory_to_copy_dlls}" + -DSRC_DIR="${CMAKE_BINARY_DIR}/bin/${CMAKE_BUILD_TYPE}" + -DDST_DIR="${directory_to_copy_dlls}" + -P "${COPY_DLLS_SCRIPT}" + COMMENT "Copying DLLs to ${directory_to_copy_dlls}" + ) endfunction() \ No newline at end of file diff --git a/src/diffCheck/geometry/DFMesh.cc b/src/diffCheck/geometry/DFMesh.cc index 875fc254..249b46f8 100644 --- a/src/diffCheck/geometry/DFMesh.cc +++ b/src/diffCheck/geometry/DFMesh.cc @@ -112,10 +112,12 @@ 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); - n.normalize(); + double n2 = n.squaredNorm(); - // Project the point onto the plane of the triangle - Eigen::Vector3d projectedPoint = point - n * (n.dot(point - v0)); + // Handle degenerate triangle + if (n2 < 1e-20){continue;}// skip this triangle + + Eigen::Vector3d projectedPoint = point - n * (n.dot(point - v0) / n2); // Compute vectors Eigen::Vector3d v0v1 = v1 - v0; @@ -130,16 +132,21 @@ namespace diffCheck::geometry double dot12 = v0v1.dot(v0p); // create u,v isoparametric mapping to the triangle where (u,v) = (1,0) if projectedPoint = v2, (u,v) = (0,1) if projectedPoint = v1 and (u,v) = (0,0) if projectedPoint = v0 - double invDenom = 1.0 / (dot00 * dot11 - dot01 * dot01); + double denom = dot00 * dot11 - dot01 * dot01; + if (std::abs(denom) < 1e-20) + continue; + + double invDenom = 1.0 / denom; double u = (dot11 * dot02 - dot01 * dot12) * invDenom; double v = (dot00 * dot12 - dot01 * dot02) * invDenom; // Check if point is in triangle - if ((u >= -associationThreshold / 100) && (v >= -associationThreshold / 100) && (u + v <= 1 + associationThreshold / 100)) + double epsilon = 1e-6; + + if ((u >= -epsilon) && (v >= -epsilon) && (u + v <= 1 + epsilon)) { // Check if the point is close enough to the face - double maxProjectionDistance = associationThreshold * std::min({(v1 - v0).norm(), (v2 - v1).norm(), (v0 - v2).norm()}) ; - if ((projectedPoint - point).norm() < maxProjectionDistance) + if ((projectedPoint - point).squaredNorm() < associationThreshold * associationThreshold) { return true; } diff --git a/src/diffCheck/geometry/DFPointCloud.cc b/src/diffCheck/geometry/DFPointCloud.cc index e7021f70..ea4662d5 100644 --- a/src/diffCheck/geometry/DFPointCloud.cc +++ b/src/diffCheck/geometry/DFPointCloud.cc @@ -156,7 +156,7 @@ namespace diffCheck::geometry } for (auto &normal : O3DPointCloud->normals_) { - if(normal.z() < -0.8) + if(normal.z() < -0.1) { normal = -normal; } @@ -173,7 +173,7 @@ namespace diffCheck::geometry this->Normals.clear(); for (int i = 0; i < cilantroPointCloud->normals.cols(); i++) { - if(cilantroPointCloud->normals.col(i).z() < -0.8) + if(cilantroPointCloud->normals.col(i).z() < -0.1) { cilantroPointCloud->normals.col(i) = -cilantroPointCloud->normals.col(i); } diff --git a/src/diffCheck/geometry/DFPointCloud.hh b/src/diffCheck/geometry/DFPointCloud.hh index fb1367fa..b317c5a4 100644 --- a/src/diffCheck/geometry/DFPointCloud.hh +++ b/src/diffCheck/geometry/DFPointCloud.hh @@ -145,13 +145,13 @@ namespace diffCheck::geometry * * @param targetSize the target size of the cloud */ - void DownsampleBySize(int targetSize); + /** * @brief Get the tight bounding box of the point cloud * - * @return std::vector A vector of two Eigen::Vector3d, with the first one being the minimum - * point and the second one the maximum point of the bounding box. + * @return std::vector A vector of eight Eigen::Vector3d, representing the corners of the bounding box. + * The order of the corners is as follows: * /// ------- x * /// /| * /// / | diff --git a/src/diffCheck/segmentation/DFSegmentation.cc b/src/diffCheck/segmentation/DFSegmentation.cc index d79cba5b..144ce311 100644 --- a/src/diffCheck/segmentation/DFSegmentation.cc +++ b/src/diffCheck/segmentation/DFSegmentation.cc @@ -255,25 +255,26 @@ namespace diffCheck::segmentation } for (auto segment : clusters) { - Eigen::Vector3d segmentCenter; - Eigen::Vector3d segmentNormal; + Eigen::Vector3d segmentNormal = Eigen::Vector3d::Zero(); - 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; } + Eigen::Vector3d segmentCenter = segment->GetAxixAlignedBoundingBox()[0] + (segment->GetAxixAlignedBoundingBox()[1] - segment->GetAxixAlignedBoundingBox()[0])/2.0; + for (auto normal : segment->Normals){segmentNormal += normal;} + if (segmentNormal.norm() == 0) + { + DIFFCHECK_WARN("Segment normal is zero. Skipping the segment."); + continue; + } 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; @@ -288,64 +289,47 @@ namespace diffCheck::segmentation } bool hasColors = correspondingSegment->GetNumColors() > 0; - for (Eigen::Vector3d point : correspondingSegment->Points) + std::vector indicesToRemove; + for (size_t i = 0; i < correspondingSegment->Points.size(); i++) { + const Eigen::Vector3d& point = correspondingSegment->Points[i]; + if (discriminatePoints) { - 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))] - ); + facePoints->Normals.push_back(correspondingSegment->Normals[i]); if (hasColors) { - facePoints->Colors.push_back( - correspondingSegment->Colors[std::distance( - correspondingSegment->Points.begin(), - std::find(correspondingSegment->Points.begin(), - correspondingSegment->Points.end(), - point))] - ); + facePoints->Colors.push_back(correspondingSegment->Colors[i]); } + indicesToRemove.push_back(i); } } else { 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))] - ); + facePoints->Normals.push_back(correspondingSegment->Normals[i]); if (hasColors) { - facePoints->Colors.push_back( - correspondingSegment->Colors[std::distance( - correspondingSegment->Points.begin(), - std::find(correspondingSegment->Points.begin(), - correspondingSegment->Points.end(), - point))] - ); + facePoints->Colors.push_back(correspondingSegment->Colors[i]); } + indicesToRemove.push_back(i); } } - for(Eigen::Vector3d point : facePoints->Points) + for (auto it = indicesToRemove.rbegin(); it != indicesToRemove.rend(); ++it) { - correspondingSegment->Points.erase( - std::remove( - correspondingSegment->Points.begin(), - correspondingSegment->Points.end(), - point), - correspondingSegment->Points.end()); + int i = *it; + + correspondingSegment->Points.erase(correspondingSegment->Points.begin() + i); + correspondingSegment->Normals.erase(correspondingSegment->Normals.begin() + i); + + if (hasColors) + { + correspondingSegment->Colors.erase(correspondingSegment->Colors.begin() + i); + } } faceSegments.push_back(facePoints); } @@ -373,7 +357,7 @@ namespace diffCheck::segmentation for (std::shared_ptr cluster : unassociatedClusters) { std::shared_ptr correspondingMeshFace; - Eigen::Vector3d clusterCenter; + Eigen::Vector3d clusterCenter = Eigen::Vector3d::Zero(); Eigen::Vector3d clusterNormal = Eigen::Vector3d::Zero(); if (cluster->GetNumPoints() == 0) @@ -391,7 +375,7 @@ namespace diffCheck::segmentation DIFFCHECK_WARN("No meshes to associate with the clusters. Skipping the cluster."); continue; } - for (Eigen::Vector3d point : cluster->Points) + for (const Eigen::Vector3d& point : cluster->Points) { clusterCenter += point; } @@ -435,7 +419,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; @@ -465,11 +451,13 @@ namespace diffCheck::segmentation 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 < maximumFaceSegmentDistance && currentDistance * (std::abs(faceNormal.dot((faceCenter - clusterCenter) / (faceCenter - clusterCenter).norm()))) < distance) + double anglePenalty = 100*std::abs(clusterNormal.dot(faceNormal)); + double currentDistance = (clusterCenter - faceCenter).norm() * (.1 + std::abs(dotProduct)) / std::max(anglePenalty, 1.0); + double normalAlignment = std::abs(faceNormal.dot(clusterNormal)); + if (std::abs(std::sqrt(1.0 - normalAlignment * normalAlignment)) < angleThreshold + && currentDistance < maximumFaceSegmentDistance + && currentDistance < distance) { goodMeshIndex = meshIndex; goodFaceIndex = faceIndex; @@ -494,13 +482,16 @@ namespace diffCheck::segmentation } std::shared_ptr completed_segment = existingPointCloudSegments[goodMeshIndex][goodFaceIndex]; - for (Eigen::Vector3d point : cluster->Points) + std::vector indicesToRemove; + for (size_t i = 0; i < cluster->Points.size(); i++) { + const Eigen::Vector3d& point = cluster->Points[i]; 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))]); + completed_segment->Normals.push_back(cluster->Normals[i]); + completed_segment->Colors.push_back(cluster->Colors[i]); + indicesToRemove.push_back(i); } else { @@ -509,27 +500,20 @@ namespace diffCheck::segmentation 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))]); + completed_segment->Normals.push_back(cluster->Normals[i]); + completed_segment->Colors.push_back(cluster->Colors[i]); + indicesToRemove.push_back(i); } } else { 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))]); + completed_segment->Normals.push_back(cluster->Normals[i]); + completed_segment->Colors.push_back(cluster->Colors[i]); + indicesToRemove.push_back(i); } } } - 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()) - { - indicesToRemove.push_back(i); - } - } for (auto it = indicesToRemove.rbegin(); it != indicesToRemove.rend(); ++it) { std::swap(cluster->Points[*it], cluster->Points.back()); diff --git a/src/gh/components/DF_CAD_segmentator/code.py b/src/gh/components/DF_CAD_segmentator/code.py index 13371761..533c3957 100644 --- a/src/gh/components/DF_CAD_segmentator/code.py +++ b/src/gh/components/DF_CAD_segmentator/code.py @@ -43,6 +43,7 @@ def RunScript(self, i_maximum_face_segment_distance = 0.1 o_face_clusters = [] + o_poses_from_icp = [] transforms = [] df_clusters = [] # we make a deepcopy of the input clouds @@ -88,6 +89,11 @@ def RunScript(self, df_asssociated_cluster_faces_per_beam = [] for i, df_b in enumerate(df_beams): + beam_detected_pose = Rhino.Geometry.Plane(df_b.plane) + if beam_detected_pose.Transform(transforms[i]): + o_poses_from_icp.append(beam_detected_pose) + else: + o_poses_from_icp.append(None) rh_b_mesh_faces = [df_b_f.to_mesh() for df_b_f in df_b.side_faces] rh_test_mesh = Rhino.Geometry.Mesh() for j in range(len(rh_b_mesh_faces)): @@ -144,4 +150,4 @@ def RunScript(self, o_face_clouds = th.list_to_tree(o_face_clusters) - return [o_beam_clouds, o_face_clouds] + return [o_beam_clouds, o_face_clouds, o_poses_from_icp] diff --git a/src/gh/components/DF_CAD_segmentator/metadata.json b/src/gh/components/DF_CAD_segmentator/metadata.json index 2146f6a1..795b2c26 100644 --- a/src/gh/components/DF_CAD_segmentator/metadata.json +++ b/src/gh/components/DF_CAD_segmentator/metadata.json @@ -126,6 +126,14 @@ "optional": false, "sourceCount": 0, "graft": false + }, + { + "name": "o_poses_from_icp", + "nickname": "o_poses_from_icp", + "description": "The list of poses resulting from the ICP registration. If i_make_registration is False, this output will contain the original assembly beam planes (no ICP applied).", + "optional": false, + "sourceCount": 0, + "graft": false } ] } diff --git a/src/gh/components/DF_build_assembly/code.py b/src/gh/components/DF_build_assembly/code.py index 0be4fa94..9b554b4f 100644 --- a/src/gh/components/DF_build_assembly/code.py +++ b/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) diff --git a/src/gh/components/DF_build_assembly/metadata.json b/src/gh/components/DF_build_assembly/metadata.json index 18c7a446..662edd9d 100644 --- a/src/gh/components/DF_build_assembly/metadata.json +++ b/src/gh/components/DF_build_assembly/metadata.json @@ -48,13 +48,25 @@ "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": [ { "name": "o_assembly", "nickname": "o_assembly", - "description": "The create DFAssembly object representing the timber elements.", + "description": "The created DFAssembly object representing the timber elements.", "optional": false, "sourceCount": 0, "graft": false diff --git a/src/gh/components/DF_csv_exporter/code.py b/src/gh/components/DF_csv_exporter/code.py index b27f5671..3f49774c 100644 --- a/src/gh/components/DF_csv_exporter/code.py +++ b/src/gh/components/DF_csv_exporter/code.py @@ -8,7 +8,7 @@ from ghpythonlib.componentbase import executingcomponent as component import Grasshopper as gh -from diffCheck.df_error_estimation import DFInvalidData, DFVizResults +from diffCheck.df_error_estimation import DFInvalidData, DFVizResults, DFPoseResults def add_bool_toggle(self, @@ -170,21 +170,50 @@ def RunScript(self, if i_dump: os.makedirs(i_export_dir, exist_ok=True) - self.prefix = i_result.analysis_type - - if i_export_seperate_files: - for idx in range(len(i_result.source)): - element_id = self._get_id(idx, i_result) - csv_analysis_path = os.path.join(i_export_dir, f"{i_file_name}_{self.prefix}_{element_id}.csv") - rows = [self._prepare_row(idx, i_result)] - self._write_csv(csv_analysis_path, rows) + if isinstance(i_result, DFVizResults): + self.prefix = i_result.analysis_type + + if i_export_seperate_files: + for idx in range(len(i_result.source)): + element_id = self._get_id(idx, i_result) + csv_analysis_path = os.path.join(i_export_dir, f"{i_file_name}_{self.prefix}_{element_id}.csv") + rows = [self._prepare_row(idx, i_result)] + self._write_csv(csv_analysis_path, rows) + if i_export_distances: + csv_distances_path = os.path.join(i_export_dir, f"{i_file_name}_{self.prefix}_{element_id}_distances.csv") + self._write_csv(csv_distances_path, rows, is_writing_only_distances=True) + else: + csv_analysis_path = os.path.join(i_export_dir, f"{i_file_name}.csv") + merged_rows = [self._prepare_row(idx, i_result) for idx in range(len(i_result.source))] + self._write_csv(csv_analysis_path, merged_rows) if i_export_distances: - csv_distances_path = os.path.join(i_export_dir, f"{i_file_name}_{self.prefix}_{element_id}_distances.csv") - self._write_csv(csv_distances_path, rows, is_writing_only_distances=True) - else: - csv_analysis_path = os.path.join(i_export_dir, f"{i_file_name}.csv") - merged_rows = [self._prepare_row(idx, i_result) for idx in range(len(i_result.source))] - self._write_csv(csv_analysis_path, merged_rows) - if i_export_distances: - csv_distances_path = os.path.join(i_export_dir, f"{i_file_name}_distances.csv") - self._write_csv(csv_distances_path, merged_rows, is_writing_only_distances=True) + csv_distances_path = os.path.join(i_export_dir, f"{i_file_name}_distances.csv") + self._write_csv(csv_distances_path, merged_rows, is_writing_only_distances=True) + + elif isinstance(i_result, DFPoseResults): + data_dicts = [] + csv_analysis_path = os.path.join(i_export_dir, f"{i_file_name}_pose_data.csv") + elem_names, elem_last_dist_err, elem_last_rot_err, assembly_dist_err_hist, assembly_rot_err_hist = i_result.compute_history_pose_errors() + for i in range(len(elem_names)): + dist_err_list = [] + rot_err_list = [] + for data in assembly_dist_err_hist[i]: + if data is not None: + dist_err_list.append(float(data)) + else: + dist_err_list.append("nan") + for data in assembly_rot_err_hist[i]: + if data: + rot_err_list.append(float(data)) + else: + rot_err_list.append("nan") + + data_dict = { + "element_name": elem_names[i], + "last_distance_error": elem_last_dist_err[i], + "last_rotation_error": elem_last_rot_err[i], + "assembly_distance_error_history": dist_err_list, + "assembly_rotation_error_history": rot_err_list + } + data_dicts.append(data_dict) + self._write_csv(csv_analysis_path, data_dicts) diff --git a/src/gh/components/DF_pose_comparison/code.py b/src/gh/components/DF_pose_comparison/code.py index 267afd9d..53b62974 100644 --- a/src/gh/components/DF_pose_comparison/code.py +++ b/src/gh/components/DF_pose_comparison/code.py @@ -7,6 +7,8 @@ import ghpythonlib.treehelpers as th import diffCheck.df_geometries +import diffCheck.df_poses +import diffCheck.df_error_estimation import numpy def compute_comparison(measured_pose, cad_pose): @@ -57,6 +59,7 @@ def RunScript(self, i_assembly: diffCheck.df_geometries.DFAssembly, i_measured_p o_distances[beam_id].append(dist) o_angles[beam_id].append(angle) o_transforms_cad_to_measured[beam_id].append(transform_cad_to_measured) + else: i_measured_planes.Flatten() measured_plane_list = th.tree_to_list(i_measured_planes) @@ -67,7 +70,13 @@ def RunScript(self, i_assembly: diffCheck.df_geometries.DFAssembly, i_measured_p o_angles.append(angle) o_transforms_cad_to_measured.append(transform_cad_to_measured) + + df_poses_assembly = diffCheck.df_poses.DFPosesAssembly() + df_poses_assembly.from_gh_tree(i_measured_planes) + o_result = diffCheck.df_error_estimation.DFPoseResults(i_assembly) + o_result.add_history(df_poses_assembly.poses_per_element_dictionary) + if bc == 1: - return o_distances, o_angles, o_transforms_cad_to_measured + return o_distances, o_angles, o_transforms_cad_to_measured, o_result else: - return th.list_to_tree(o_distances), th.list_to_tree(o_angles), th.list_to_tree(o_transforms_cad_to_measured) + return th.list_to_tree(o_distances), th.list_to_tree(o_angles), th.list_to_tree(o_transforms_cad_to_measured), o_result diff --git a/src/gh/components/DF_pose_comparison/metadata.json b/src/gh/components/DF_pose_comparison/metadata.json index c9ad3a8e..e3130c95 100644 --- a/src/gh/components/DF_pose_comparison/metadata.json +++ b/src/gh/components/DF_pose_comparison/metadata.json @@ -61,6 +61,14 @@ "optional": false, "sourceCount": 0, "graft": false + }, + { + "name": "o_result", + "nickname": "o_result", + "description": "A DFPoseResults object containing detailed error information for each pose comparison.", + "optional": false, + "sourceCount": 0, + "graft": false } ] } diff --git a/src/gh/components/DF_pose_estimation/code.py b/src/gh/components/DF_pose_estimation/code.py index 60af0e56..b650869f 100644 --- a/src/gh/components/DF_pose_estimation/code.py +++ b/src/gh/components/DF_pose_estimation/code.py @@ -17,6 +17,7 @@ class DFPoseEstimation(component): def RunScript(self, i_face_clouds: Grasshopper.DataTree[Rhino.Geometry.PointCloud], i_assembly, + i_poses_from_icp: list, i_reset: bool, i_save: bool): @@ -35,40 +36,47 @@ def RunScript(self, all_poses_this_time = [] for i, face_clouds in enumerate(clusters_per_beam): try: - df_cloud = dfb_geometry.DFPointCloud() - - 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])) - - df_bb_points = df_cloud.get_tight_bounding_box() - df_bb_centroid = sum(df_bb_points)/len(df_bb_points) - 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_beam_cloud = Rhino.Geometry.PointCloud() - for face_cloud in face_clouds: - rh_beam_cloud.Merge(face_cloud) - - rh_bbox = rh_beam_cloud.GetBoundingBox(rh_tentative_plane) - rh_bbox.Transform(Rhino.Geometry.Transform.PlaneToPlane(Rhino.Geometry.Plane.WorldXY, rh_tentative_plane)) - rh_bb_centroid = rh_bbox.Center - - pose = df_poses.DFPose( - 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_bb_centroid, xDirection=new_xDirection, yDirection=new_yDirection) - planes.append(plane) + if i_poses_from_icp and len(i_poses_from_icp) > i and i_poses_from_icp[i] is not None: + # if there is a pose from ICP for this beam, use it directly without processing the cloud + planes.append(i_poses_from_icp[i]) + all_poses_this_time.append(df_poses.DFPose.from_rh_plane(i_poses_from_icp[i])) + continue + + else: + df_cloud = dfb_geometry.DFPointCloud() + + 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])) + + df_bb_points = df_cloud.get_tight_bounding_box() + df_bb_centroid = sum(df_bb_points)/len(df_bb_points) + 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_xDirection, new_yDirection) + + rh_beam_cloud = Rhino.Geometry.PointCloud() + for face_cloud in face_clouds: + rh_beam_cloud.Merge(face_cloud) + + rh_bbox = rh_beam_cloud.GetBoundingBox(rh_tentative_plane) + rh_bbox.Transform(Rhino.Geometry.Transform.PlaneToPlane(Rhino.Geometry.Plane.WorldXY, rh_tentative_plane)) + rh_bb_centroid = rh_bbox.Center + + pose = df_poses.DFPose( + 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_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 diff --git a/src/gh/components/DF_pose_estimation/metadata.json b/src/gh/components/DF_pose_estimation/metadata.json index f7c780ae..f09bc53f 100644 --- a/src/gh/components/DF_pose_estimation/metadata.json +++ b/src/gh/components/DF_pose_estimation/metadata.json @@ -37,6 +37,18 @@ "sourceCount": 0, "typeHintID": "ghdoc" }, + { + "name": "i_poses_from_icp", + "nickname": "i_poses_from_icp", + "description": "The optional poses detected in the CAD segmentation component though ICP. If provided, the pose will not be re-calculated using the faces normals.", + "optional": true, + "allowTreeAccess": true, + "showTypeHints": true, + "scriptParamAccess": "list", + "wireDisplay": "default", + "sourceCount": 0, + "typeHintID": "ghdoc" + }, { "name": "i_reset", "nickname": "i_reset", diff --git a/src/gh/diffCheck/diffCheck/df_error_estimation.py b/src/gh/diffCheck/diffCheck/df_error_estimation.py index b9ab03c0..e61d79b4 100644 --- a/src/gh/diffCheck/diffCheck/df_error_estimation.py +++ b/src/gh/diffCheck/diffCheck/df_error_estimation.py @@ -19,6 +19,7 @@ from diffCheck import diffcheck_bindings # type: ignore from diffCheck import df_cvt_bindings from diffCheck.df_geometries import DFAssembly +from diffCheck.df_poses import DFPosesBeam, DFPose class NumpyEncoder(json.JSONEncoder): @@ -268,6 +269,64 @@ def analysis_type(self): self._analysis_type = self._compute_dfresult_type() return self._analysis_type +class DFPoseResults(): + """ + This class compiles the results of the pose estimation into one object + """ + def __init__(self, assembly: DFAssembly): + self.assembly = assembly + self.pose_history : dict[str, DFPosesBeam] = dict() + self.last_poses : dict[str, DFPose] = dict() + + def add_history(self, pose_history : dict[str, DFPosesBeam]): + """ + The pose history is a dictionnary where the keys are the element names ("element_0", "element_1", etc), + and the values are DFPosesBeam objects containing a dictionnary of poses for each element. + """ + self.pose_history = pose_history + for element in pose_history: + poses_dict = pose_history[element].poses_dictionary + self.last_poses[element] = poses_dict[next(reversed(poses_dict))] if poses_dict else None + + def add_last_poses(self, last_poses : dict[str, DFPose]): + """ + Adds a dictionnary of the last poses for each element. The keys are the element names_0", "element_1", etc), + """ + self.last_poses = last_poses + + def compute_history_pose_errors(self): + """ + This function computes the error of the pose estimation for each element at each step, compared to the poses defined in the assembly. + """ + element_names = [] + element_last_dist_errors = [] + element_last_rot_errors = [] + assembly_dist_error_history = [] + assembly_rot_error_history = [] + + for element_name, df_poses_beam in self.pose_history.items(): + element_index = int(element_name.split("_")[1]) + df_beam = self.assembly.beams[element_index] + df_beam_pose_plane = df_beam.plane + dist_error_history = [] + rot_error_history = [] + for pose_name, pose in df_poses_beam.poses_dictionary.items(): + if pose is None: + dist_error_history.append(None) + rot_error_history.append(None) + continue + else: + dist, angle, transform_error = pose.compare_to_rh_plane(df_beam_pose_plane) + dist_error_history.append(dist) + rot_error_history.append(angle) + assembly_dist_error_history.append(dist_error_history) + assembly_rot_error_history.append(rot_error_history) + element_names.append(element_name) + element_last_dist_errors.append(dist_error_history[-1]) + element_last_rot_errors.append(rot_error_history[-1]) + return element_names, element_last_dist_errors, element_last_rot_errors, assembly_dist_error_history, assembly_rot_error_history + + # FIXME: ths is currently broken, we need to fix it def df_cloud_2_df_cloud_comparison( assembly: DFAssembly, diff --git a/src/gh/diffCheck/diffCheck/df_geometries.py b/src/gh/diffCheck/diffCheck/df_geometries.py index 99d4c531..4931a5e8 100644 --- a/src/gh/diffCheck/diffCheck/df_geometries.py +++ b/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) @@ -232,7 +241,7 @@ def to_mesh(self): for mesh_part in mesh_parts: mesh.Append(mesh_part) mesh.Faces.ConvertQuadsToTriangles() - # mesh.Compact() + mesh.Compact() return mesh @@ -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) diff --git a/src/gh/diffCheck/diffCheck/df_joint_detector.py b/src/gh/diffCheck/diffCheck/df_joint_detector.py index b8fa4678..93a1d454 100644 --- a/src/gh/diffCheck/diffCheck/df_joint_detector.py +++ b/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,7 +134,7 @@ 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. @@ -140,24 +142,20 @@ def run(self): :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()] diff --git a/src/gh/diffCheck/diffCheck/df_poses.py b/src/gh/diffCheck/diffCheck/df_poses.py index 98f836c4..8455c84a 100644 --- a/src/gh/diffCheck/diffCheck/df_poses.py +++ b/src/gh/diffCheck/diffCheck/df_poses.py @@ -3,6 +3,7 @@ import Rhino import json +import numpy from dataclasses import dataclass, field # use a key and not all the sticky @@ -21,6 +22,19 @@ class DFPose: xDirection: list yDirection: list + @staticmethod + def from_rh_plane(rh_plane): + """ + Create a DFPose object from a Rhino Plane object. + + :param rh_plane: the Rhino Plane to convert + :return: a DFPose object representing the same pose as the input Rhino Plane + """ + return DFPose( + origin = [rh_plane.Origin.X, rh_plane.Origin.Y, rh_plane.Origin.Z], + xDirection = [rh_plane.XAxis.X, rh_plane.XAxis.Y, rh_plane.XAxis.Z], + yDirection = [rh_plane.YAxis.X, rh_plane.YAxis.Y, rh_plane.YAxis.Z]) + def to_rh_plane(self): """ Convert the pose to a Rhino Plane object. @@ -30,6 +44,31 @@ def to_rh_plane(self): yDirection = Rhino.Geometry.Vector3d(self.yDirection[0], self.yDirection[1], self.yDirection[2]) return Rhino.Geometry.Plane(origin, xDirection, yDirection) + def compare_to_rh_plane(self, rh_plane): + """ + Compare this pose to another pose and return the differences in origin, xDirection and yDirection. + + :param rh_plane: the Rhino Plane to compare to + :return: a tuple containing the distance between the origins, the angle between the xDirections and the rhino transform to go from the compared pose to this pose. + """ + other_origin = rh_plane.Origin + measured_origin = self.to_rh_plane().Origin + distance = other_origin.DistanceTo(measured_origin) + + # Compare the orientations using the formula: $$ \theta = \arccos\left(\frac{\text{trace}(R_{\text{pred}}^T R_{\text{meas}}) - 1}{2}\right) $$ + transform_o_to_other = Rhino.Geometry.Transform.PlaneToPlane(Rhino.Geometry.Plane.WorldXY, rh_plane) + transform_o_to_current = Rhino.Geometry.Transform.PlaneToPlane(Rhino.Geometry.Plane.WorldXY, self.to_rh_plane()) + np_transform_o_to_other = numpy.array(transform_o_to_other.ToDoubleArray(rowDominant=True)).reshape((4, 4)) + np_transform_o_to_measured = numpy.array(transform_o_to_current.ToDoubleArray(rowDominant=True)).reshape((4, 4)) + + R_other = np_transform_o_to_other[:3, :3] + R_measured = np_transform_o_to_measured[:3, :3] + R_rel = numpy.dot(R_other.T, R_measured) + theta = numpy.arccos(numpy.clip((numpy.trace(R_rel) - 1) / 2, -1.0, 1.0)) + + transform_other_to_current_plane = Rhino.Geometry.Transform.PlaneToPlane(rh_plane, self.to_rh_plane()) + return distance, theta, transform_other_to_current_plane + @dataclass class DFPosesBeam: """ @@ -117,6 +156,32 @@ def to_gh_tree(self): list_of_poses.append(list_of_pose_of_element) return th.list_to_tree(list_of_poses) + def from_gh_tree(self, gh_tree): + """ + Load the assembly poses from a Grasshopper tree structure. + + :param gh_tree: the Grasshopper tree containing the poses in the form of Rhino Planes + """ + self.reset() + bc = gh_tree.BranchCount + if bc > 1: + list_of_poses = th.tree_to_list(gh_tree) + else: + gh_tree.Flatten() + list_of_poses = [th.tree_to_list(gh_tree)] + n_poses = len(list_of_poses[0]) if list_of_poses else 0 + for i in range(n_poses): + new_poses = [] + for poses_of_element in list_of_poses: + if poses_of_element[i] is None: + new_poses.append(None) + continue + new_poses.append(DFPose( + origin = [poses_of_element[i].Origin.X, poses_of_element[i].Origin.Y, poses_of_element[i].Origin.Z], + xDirection = [poses_of_element[i].XAxis.X, poses_of_element[i].XAxis.Y, poses_of_element[i].XAxis.Z], + yDirection = [poses_of_element[i].YAxis.X, poses_of_element[i].YAxis.Y, poses_of_element[i].YAxis.Z])) + self.add_step(new_poses) + def compute_dot_product(v1, v2): """ @@ -134,23 +199,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 diff --git a/src/gh/diffCheck/diffCheck/df_util.py b/src/gh/diffCheck/diffCheck/df_util.py index ba35f5e4..ce079f27 100644 --- a/src/gh/diffCheck/diffCheck/df_util.py +++ b/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())