From 8f8f49111826b7f8805340eef31f00a158c0b21c Mon Sep 17 00:00:00 2001 From: PedroCorcaque Date: Wed, 1 May 2024 16:19:21 -0300 Subject: [PATCH 1/5] Added AsTriangulation as a module of the GeometricPrimitiveFittingTools --- dataset_generator.py | 23 ++++ lib/surfaces_projector.py | 104 +++++++++++++++++ lib/triangulation.py | 229 ++++++++++++++++++++++++++++++++++++++ requirements.txt | 8 ++ 4 files changed, 364 insertions(+) create mode 100644 lib/surfaces_projector.py create mode 100644 lib/triangulation.py create mode 100644 requirements.txt diff --git a/dataset_generator.py b/dataset_generator.py index 5eb30b3..6aad1d7 100644 --- a/dataset_generator.py +++ b/dataset_generator.py @@ -13,6 +13,7 @@ from lib.utils import loadFeatures, computeLabelsFromFace2Primitive, savePCD, downsampleByPointIndices, rayCastingPointCloudGeneration, funif, computeFeaturesPointIndices from lib.writers import DatasetWriterFactory +from lib import triangulation from asGeometryOCCWrapper.surfaces import SurfaceFactory @@ -58,8 +59,22 @@ parser.add_argument('-m_np', '--min_number_points', type=float, default = 0.0001, help='filter geometries by number of points.') parser.add_argument('-ls', '--leaf_size', type=float, default = 0.0, help='') + parser.add_argument('--generate_mesh', action='store_true', default=False, help='Generate mesh from point cloud') + parser.add_argument('--mesh_foldername', type=str, default='mesh', help='') + parser.add_argument('--triangulation_foldername', type=str, default='triangulation', help='') + parser.add_argument('--triangulation_features_foldername', type=str, default='triangulation_features', help='') + parser.add_argument('-gf', "--generate_triangulation_features", action="store_true", help="[Optional] Used in mesh generation.") + parser.add_argument('-pp', "--project_points", action="store_true", help="[Optional] Used in mesh generation.") + parser.add_argument('-r', "--reegenerate", action="store_true", help="[Optional] Used in mesh generation.") + parser.add_argument('-nf', "--no_filter", action="store_true", help="[Optional] Used in mesh generation.") + parser.add_argument('--n_neighbors_mesh', type=int, default=N_NEIGHBORS, help="[Optional] Used in mesh generation.") + parser.add_argument('--leaf_size_mesh', type=float, default=LEAF_SIZE, help="[Optional] Used in mesh generation.") + parser.add_argument('--max_workers_mesh', type=int, default=MAX_WORKERS, help="[Optional] Used in mesh generation.") + args = vars(parser.parse_args()) + generate_mesh = args.generate_mesh + folder_name = args['folder'] formats = [s.lower() for s in args['formats'].split(',')] curve_types = [s.lower() for s in args['curve_types'].split(',')] @@ -175,6 +190,14 @@ labels, features_point_indices = computeLabelsFromFace2Primitive(labels_mesh.copy(), features_data['surfaces']) print('Done.\n') + if generate_mesh: + print("Generating mesh from point cloud...") + if (generate_mesh_from_pointcloud(args, pc_filename, features_data)): + print("\tMesh generated on: {}".format(mesh_folder_name)) + else: + print("\tThe mesh could not be generated.") + print("Done") + elif mesh_filename is not None: print('Opening Mesh...') mesh = o3d.io.read_triangle_mesh(mesh_filename, print_progress=True) diff --git a/lib/surfaces_projector.py b/lib/surfaces_projector.py new file mode 100644 index 0000000..a92082c --- /dev/null +++ b/lib/surfaces_projector.py @@ -0,0 +1,104 @@ +from OCC.Core.Geom import Geom_CylindricalSurface, Geom_ConicalSurface, Geom_Plane, Geom_SphericalSurface, Geom_BSplineSurface, Geom_SurfaceOfLinearExtrusion, Geom_SurfaceOfRevolution, Geom_ToroidalSurface +from OCC.Core.GeomAPI import GeomAPI_ProjectPointOnSurf +from OCC.Core.gp import gp_Ax3, gp_Pnt, gp_Dir +from OCC.Core.TColgp import TColgp_Array2OfPnt +from OCC.Core.TColStd import TColStd_Array2OfReal, TColStd_Array1OfReal, TColStd_Array1OfInteger + +import numpy as np + +class SurfacesProjector: + PROJECTOR = GeomAPI_ProjectPointOnSurf() + + def buildCylinder(features): + radius = features['radius'] + location = features['location'] + z_axis = features['z_axis'] + axis = gp_Ax3(gp_Pnt(location[0], location[1], location[2]), gp_Dir(z_axis[0], z_axis[1], z_axis[2])) + surf = Geom_CylindricalSurface(axis, radius) + return surf + + def buildBSpline(features): + poles = np.asarray(features['poles']) + poles_tcol = TColgp_Array2OfPnt(1, poles.shape[0], 1, poles.shape[1]) + for i in range(poles.shape[0]): + for j in range(poles.shape[1]): + print(i, j) + e = poles[i, j, :] + #poles_tcol.SetValue(i, j, gp_Pnt(e[0], e[1], e[2]) ) + + weights = features['weights'] + u_knots = features['u_knots'] + v_knots = features['v_knots'] + u_degree = features['u_degree'] + v_degree = features['v_degree'] + surf = Geom_BSplineSurface(poles_tcol, TColStd_Array2OfReal(weights), TColStd_Array1OfReal(u_knots), TColStd_Array1OfReal(v_knots), TColStd_Array1OfInteger([1 for i in range(len(u_knots))]), TColStd_Array1OfInteger([1 for i in range(len(v_knots))]), u_degree, v_degree) + return surf + + def buildCone(features): + radius = features['radius'] + location = features['location'] + z_axis = features['z_axis'] + angle = features['angle'] + axis = gp_Ax3(gp_Pnt(location[0], location[1], location[2]), gp_Dir(z_axis[0], z_axis[1], z_axis[2])) + surf = Geom_ConicalSurface(axis, angle, radius) + return surf + + def buildPlane(features): + location = features['location'] + normal = features['z_axis'] + axis = gp_Ax3(gp_Pnt(location[0], location[1], location[2]), gp_Dir(normal[0], normal[1], normal[2])) + surf = Geom_Plane(axis) + return surf + + def buildSphere(features): + radius = features['radius'] + location = features['location'] + z_axis = features['z_axis'] + axis = gp_Ax3(gp_Pnt(location[0], location[1], location[2]), gp_Dir(z_axis[0], z_axis[1], z_axis[2])) + surf = Geom_SphericalSurface(axis, radius) + return surf + + def buildSurfaceRevolution(features): + pass + + def buildTorus(features): + min_radius = features['min_radius'] + max_radius = features['max_radius'] + location = features['location'] + z_axis = features['z_axis'] + axis = gp_Ax3(gp_Pnt(location[0], location[1], location[2]), gp_Dir(z_axis[0], z_axis[1], z_axis[2])) + surf = Geom_ToroidalSurface(axis, max_radius, min_radius) + return surf + + def buildSurfaceExtrusion(features): + pass + + BUILD_FUNCTION_MAP = { + 'cylinder': buildCylinder, + 'cone': buildCone, + 'plane': buildPlane, + 'sphere': buildSphere, + 'torus': buildTorus + } + + @staticmethod + def projectPointsOnSurfaceFeatures(points, features): + tp = features['type'].lower() + if tp in SurfacesProjector.BUILD_FUNCTION_MAP.keys(): + if len(points) == 0: + return [], [] + surface = SurfacesProjector.BUILD_FUNCTION_MAP[tp](features) + uvs = np.zeros((points.shape[0], 2), dtype=np.float64) + points_projected = np.zeros(points.shape, dtype=np.float64) + SurfacesProjector.PROJECTOR.Init(gp_Pnt(points[0, 0], points[0, 1], points[0, 2]), surface) + uvs[0] = np.array(SurfacesProjector.PROJECTOR.LowerDistanceParameters()) + pnt = SurfacesProjector.PROJECTOR.NearestPoint() + points_projected[0] = np.array([pnt.X(), pnt.Y(), pnt.Z()]) + for idx, pt in enumerate(points[1:, :]): + SurfacesProjector.PROJECTOR.Perform(gp_Pnt(pt[0], pt[1], pt[2])) + uvs[idx+1] = np.array(SurfacesProjector.PROJECTOR.LowerDistanceParameters()) + pnt = SurfacesProjector.PROJECTOR.NearestPoint() + points_projected[idx+1] = np.array([pnt.X(), pnt.Y(), pnt.Z()]) + return uvs, points_projected + else: + return None \ No newline at end of file diff --git a/lib/triangulation.py b/lib/triangulation.py new file mode 100644 index 0000000..f2c36ab --- /dev/null +++ b/lib/triangulation.py @@ -0,0 +1,229 @@ +import os +import numpy as np +import open3d as o3d +from scipy.spatial import Delaunay +from functools import partial +from tqdm.contrib.concurrent import thread_map +from pypcd import pypcd + +from .surfaces_projector import SurfacesProjector + +# Constants +N_NEIGHBORS = 10 +FILTER = True # Apply voxel grid and remove statistical outlier? +LEAF_SIZE = 0.02 # Voxel size +PROJECT_POINTS = False +MAX_WORKERS = 8 + +def compute_labels_from_face_2_primitive(labels: list, features_data: dict): + max_face = np.max(labels) + for feature in features_data: + max_face = max(0 if len(feature["face_indices"]) == 0 else max(feat["face_indices"]), max_face) + face_2_primitive = np.zeros(shape=(max_face+1,), dtype=np.int32) - 1 + face_primitive_count = np.zeros(shape=(max_face+1,), dtype=np.int32) + for feature_idx, feature in enumerate(features_data): + for face_id in feature["face_indices"]: + face_2_primitive[face_id] = feature_idx + face_primitive_count[face_id] += 1 + assert len(np.unique(face_primitive_count)) <= 2 + + features_point_indices = [[] for i in range(0, len(features_data)+1)] + for i in range(0, len(labels)): + index = face_2_primitive[labels[i]] + features_point_indices[index].append(i) + labels[i] = index + features_point_indices.pop(-1) + + for i in range(0, len(features_point_indices)): + features_point_indices[i] = np.array(features_point_indices[i], dtype=np.int64) + + return labels, features_point_indices + +def compute_local_densities(pcd: o3d.geometry.PointCloud, k: int = N_NEIGHBORS): + if len(pcd.points) <= 1: + return 0, 0, 0 + + pcd_tree = o3d.geometry.KDTreeFlann(pcd) + points = np.asarray(pcd.points) + distances = [] + for i in range(len(points)): + A = points[i] + _, idx, _ = pcd_tree.search_knn_vector_3d(A, k) + points_query = points[np.asarray(idx)] + B = points_query[-1] + distances.append(np.linalg.norm(B-A, ord=2)) + dist_arr = np.array(distances) + return np.min(dist_arr, axis=0), np.mean(dist_arr, axis=0), np.max(dist_arr, axis=0) + + +def filter_pcd(pcd: o3d.geometry.PointCloud) -> o3d.geometry.PointCloud: + if FILTER: + pcd = pcd.voxel_down_sample(voxel_size=LEAF_SIZE) + pcd, _ = pcd.remove_statistical_outlier(nb_neighbors=N_NEIGHBORS, std_ratio=2.0) + return pcd + +# def delaunay_uv_triangulation(pcd: o3d.geometry.PointCloud, surface) -> o3d.geometry.TriangleMesh: +# pcd = filter_pcd(pcd) +# if len(pcd.points) <= 2: +# return o3d.geometry.TriangleMesh() + +# _, _, max_r = compute_local_densities(pcd) + +# points = np.asarray(pcd.points) + +# response = SurfacesProjector.projectPointsOnSurfaceFeatures(points, surface) + +# if response is None: +# return o3d.geometry.TriangleMesh() + +# uvs, _ = response + +# tri = Delaunay(uvs) +# triangles = tri.simplices + +# triangles_filter = [] +# index_perm = [(0,1), (0,2), (1,2)] +# for triangle in triangles: +# keep = True +# for i, j in index_perm: +# A = points[triangle[i]] +# B = points[triangle[j]] +# dist = np.linalg.norm(B-A, ord=2) +# if dist > max_r: +# keep = False +# break +# triangles_filter.append(keep) +# triangles_filter = np.array(triangles_filter) + +# triangles = triangles[triangles_filter] + +# mesh = o3d.geometry.TriangleMesh() +# mesh.vertices = o3d.utility.Vector3dVector(points) +# mesh.triangles = o3d.utility.Vector3iVector(triangles) + +# plt.scatter(uvs[:,0], uvs[:,1]) +# plt.show() + +# return mesh + +def bpa_triangulation(pcd: o3d.geometry.PointCloud, surface = None) -> o3d.geometry.TriangleMesh: + pcd = filter_pcd(pcd) + if len(pcd.points) <= 2: + return o3d.geometry.TriangleMesh() + + _, _, max_r = compute_local_densities(pcd) + radii = [max_r/2, max_r] + mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_ball_pivoting( + pcd, o3d.utility.DoubleVector(radii) + ) + return mesh + +def poisson_triangulation(pcd: o3d.geometry.PointCloud, surface = None) -> o3d.geometry.TriangleMesh: + pcd = filter_pcd(pcd) + if len(pcd.points) <= 2: + return o3d.geometry.TriangleMesh() + + mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson( + pcd, depth=9 + ) + return mesh + +def triangulation(pcd: o3d.geometry.PointCloud, data: dict) -> o3d.geometry.TriangleMesh: + pcd_surf = pcd.select_by_index(data["point_indices"]) + + mesh = o3d.geometry.TriangleMesh() + + if len(pcd_surf.points) > 0: + result = SurfacesProjector.projectPointsOnSurfaceFeatures(np.asarray(pcd_surf.points), data) + if result is not None: + uvs, points = result + data["point_parameters"] = uvs + if PROJECT_POINTS: + pcd_surf.points = o3d.utility.Vector3dVector(points) + mesh = bpa_triangulation(pcd_surf, surface=data) + + return mesh + +def triangulation_by_surface(pcd: o3d.geometry.PointCloud, surfaces_data: dict) -> o3d.geometry.TriangleMesh: + func = partial(triangulation, pcd) + result_ = thread_map(func, surfaces_data, max_workers=MAX_WORKERS, chunk_size=1) + + final_mesh = o3d.geometry.TriangleMesh() + for result_idx, result in enumerate(result_): + surfaces_data[result_idx]["vert_indices"] = list(range(len(final_mesh.vertices), + len(final_mesh.vertices)+len(result.vertices))) + del surfaces_data[result_idx]["vert_parameters"] + res = SurfacesProjector.projectPointsOnSurfaceFeatures(np.asarray(result.vertices), + surfaces_data[result_idx]) + if res is not None: + uvs, _ = res + surfaces_data[result_idx]["vert_parameters"] = uvs if type(uvs) is list else uvs.tolist() + surfaces_data[result_idx]["face_indices"] = list(range(len(final_mesh.triangles), + len(final_mesh.trangles)+len(result.triangles))) + surfaces_data[result_idx]["point_indices"] = surfaces_data[result_idx]["point_indices"] \ + if type(surfaces_data[result_idx]["point_indices"]) is list \ + else surfaces_data[i]['point_indices'].tolist() + if "point_parameters" in surfaces_data[result_idx]: + surfaces_data[result_idx]["point_parameters"] = surfaces_data[result_idx]["point_parameters"] \ + if type(surfaces_data[result_idx]["point_parameters"]) is list \ + else surfaces_data[result_idx]["point_parameters"].tolist() + + final_mesh += result + final_mesh.compute_triangle_normals() + return final_mesh + +def generate_mesh_from_pointcloud(args, pc: pypcd.PointCloud, pc_filename: str, geometry_data: dict) -> bool: + filename = pc_filename + foldername = args['mesh_foldername'] + triangulation_foldername = args['triangulation_foldername'] + triangulation_features_foldername = args['triangulation_features_foldername'] + GENERATE_TRIANGULATION_FEATURES = args['generate_triangulation_features'] + PROJECT_POINTS = args['project_points'] + N_NEIGHBORS = args['n_neighbors_mesh'] + LEAF_SIZE = args['leaf_size_mesh'] + MAX_WORKERS = args['max_workers_mesh'] + FILTER = not args['no_filter'] + reegenerate = args['reegenerate'] + + print() + print("Generating mesh for the file: {}".format(filename)) + + pointcloud = o3d.geometry.PointCloud() + + if 'x' in pc.dtype.names and 'y' in pc.dtype.names and 'z' in pc.dtype.names: + pointcloud.points = o3d.utility.Vector3dVector(np.vstack((pc['x'], pc['y'], pc['z'])).T) + else: + print("There is no point in the input pointcloud") + return False + + if 'normal_x' in pc.dtype.names and 'normal_y' in pc.dtype.names and 'normal_z' in pc.dtype.names: + pointcloud.normals = o3d.utility.Vector3dVector(np.vstack((pc['normal_x'], pc['normal_y'], pc['normal_z'])).T) + else: + print('There is no normal in the input pointcloud') + + has_labels = False + if 'label' in pc.dtype.names: + labels = pc['label'] + has_labels = True + else: + print("there is no label in the input pointcloud") + + surfaces_data = geometry_data['surfaces'] + + if has_labels: + labels, fpi = compute_labels_from_face_2_primitive(labels, surfaces_data) + for i in range(len(fpi)): + surfaces_data[i]['point_indices'] = fpi[i] + + if not pointcloud.has_normals(): + pointcloud.estimate_normals() + + mesh = triangulation_by_surface(pointcloud, surfaces_data) + igl.write_triangle_mesh(os.path.join(folder_name, triangulation_foldername, filename + '.obj'), np.asarray(mesh.vertices), np.asarray(mesh.triangles)) + + if GENERATE_TRIANGULATION_FEATURES: + features = {'surfaces': surfaces_data} + with open(os.path.join(folder_name, triangulation_features_foldername, filename + '.json'), 'w') as f: + json.dump(features, f, indent=4) + + return True diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..b84e8bb --- /dev/null +++ b/requirements.txt @@ -0,0 +1,8 @@ +numba +numpy +pypcd @ git+https://github.com/klintan/pypcd.git +open3d +scikit-learn +matplotlib +tqdm +h5py \ No newline at end of file From 331a699f66310986cb82acafa5d498e350a8db10 Mon Sep 17 00:00:00 2001 From: PedroCorcaque Date: Wed, 1 May 2024 16:29:23 -0300 Subject: [PATCH 2/5] Change the args object to dict --- dataset_generator.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/dataset_generator.py b/dataset_generator.py index 6aad1d7..5778e05 100644 --- a/dataset_generator.py +++ b/dataset_generator.py @@ -67,13 +67,13 @@ parser.add_argument('-pp', "--project_points", action="store_true", help="[Optional] Used in mesh generation.") parser.add_argument('-r', "--reegenerate", action="store_true", help="[Optional] Used in mesh generation.") parser.add_argument('-nf', "--no_filter", action="store_true", help="[Optional] Used in mesh generation.") - parser.add_argument('--n_neighbors_mesh', type=int, default=N_NEIGHBORS, help="[Optional] Used in mesh generation.") - parser.add_argument('--leaf_size_mesh', type=float, default=LEAF_SIZE, help="[Optional] Used in mesh generation.") - parser.add_argument('--max_workers_mesh', type=int, default=MAX_WORKERS, help="[Optional] Used in mesh generation.") + parser.add_argument('--n_neighbors_mesh', type=int, default=-1, help="[Optional] Used in mesh generation.") + parser.add_argument('--leaf_size_mesh', type=float, default=-1, help="[Optional] Used in mesh generation.") + parser.add_argument('--max_workers_mesh', type=int, default=-1, help="[Optional] Used in mesh generation.") args = vars(parser.parse_args()) - generate_mesh = args.generate_mesh + generate_mesh = args["generate_mesh"] folder_name = args['folder'] formats = [s.lower() for s in args['formats'].split(',')] From 55a90c9b43c3b6dbc9ac79f14f4c935162274e62 Mon Sep 17 00:00:00 2001 From: PedroLimaCorcaque Date: Mon, 13 May 2024 14:21:11 -0300 Subject: [PATCH 3/5] Change the way to call the mesh_generator --- lib/triangulation.py | 229 ------------------------------- mesh_generator.py | 311 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 311 insertions(+), 229 deletions(-) delete mode 100644 lib/triangulation.py create mode 100644 mesh_generator.py diff --git a/lib/triangulation.py b/lib/triangulation.py deleted file mode 100644 index f2c36ab..0000000 --- a/lib/triangulation.py +++ /dev/null @@ -1,229 +0,0 @@ -import os -import numpy as np -import open3d as o3d -from scipy.spatial import Delaunay -from functools import partial -from tqdm.contrib.concurrent import thread_map -from pypcd import pypcd - -from .surfaces_projector import SurfacesProjector - -# Constants -N_NEIGHBORS = 10 -FILTER = True # Apply voxel grid and remove statistical outlier? -LEAF_SIZE = 0.02 # Voxel size -PROJECT_POINTS = False -MAX_WORKERS = 8 - -def compute_labels_from_face_2_primitive(labels: list, features_data: dict): - max_face = np.max(labels) - for feature in features_data: - max_face = max(0 if len(feature["face_indices"]) == 0 else max(feat["face_indices"]), max_face) - face_2_primitive = np.zeros(shape=(max_face+1,), dtype=np.int32) - 1 - face_primitive_count = np.zeros(shape=(max_face+1,), dtype=np.int32) - for feature_idx, feature in enumerate(features_data): - for face_id in feature["face_indices"]: - face_2_primitive[face_id] = feature_idx - face_primitive_count[face_id] += 1 - assert len(np.unique(face_primitive_count)) <= 2 - - features_point_indices = [[] for i in range(0, len(features_data)+1)] - for i in range(0, len(labels)): - index = face_2_primitive[labels[i]] - features_point_indices[index].append(i) - labels[i] = index - features_point_indices.pop(-1) - - for i in range(0, len(features_point_indices)): - features_point_indices[i] = np.array(features_point_indices[i], dtype=np.int64) - - return labels, features_point_indices - -def compute_local_densities(pcd: o3d.geometry.PointCloud, k: int = N_NEIGHBORS): - if len(pcd.points) <= 1: - return 0, 0, 0 - - pcd_tree = o3d.geometry.KDTreeFlann(pcd) - points = np.asarray(pcd.points) - distances = [] - for i in range(len(points)): - A = points[i] - _, idx, _ = pcd_tree.search_knn_vector_3d(A, k) - points_query = points[np.asarray(idx)] - B = points_query[-1] - distances.append(np.linalg.norm(B-A, ord=2)) - dist_arr = np.array(distances) - return np.min(dist_arr, axis=0), np.mean(dist_arr, axis=0), np.max(dist_arr, axis=0) - - -def filter_pcd(pcd: o3d.geometry.PointCloud) -> o3d.geometry.PointCloud: - if FILTER: - pcd = pcd.voxel_down_sample(voxel_size=LEAF_SIZE) - pcd, _ = pcd.remove_statistical_outlier(nb_neighbors=N_NEIGHBORS, std_ratio=2.0) - return pcd - -# def delaunay_uv_triangulation(pcd: o3d.geometry.PointCloud, surface) -> o3d.geometry.TriangleMesh: -# pcd = filter_pcd(pcd) -# if len(pcd.points) <= 2: -# return o3d.geometry.TriangleMesh() - -# _, _, max_r = compute_local_densities(pcd) - -# points = np.asarray(pcd.points) - -# response = SurfacesProjector.projectPointsOnSurfaceFeatures(points, surface) - -# if response is None: -# return o3d.geometry.TriangleMesh() - -# uvs, _ = response - -# tri = Delaunay(uvs) -# triangles = tri.simplices - -# triangles_filter = [] -# index_perm = [(0,1), (0,2), (1,2)] -# for triangle in triangles: -# keep = True -# for i, j in index_perm: -# A = points[triangle[i]] -# B = points[triangle[j]] -# dist = np.linalg.norm(B-A, ord=2) -# if dist > max_r: -# keep = False -# break -# triangles_filter.append(keep) -# triangles_filter = np.array(triangles_filter) - -# triangles = triangles[triangles_filter] - -# mesh = o3d.geometry.TriangleMesh() -# mesh.vertices = o3d.utility.Vector3dVector(points) -# mesh.triangles = o3d.utility.Vector3iVector(triangles) - -# plt.scatter(uvs[:,0], uvs[:,1]) -# plt.show() - -# return mesh - -def bpa_triangulation(pcd: o3d.geometry.PointCloud, surface = None) -> o3d.geometry.TriangleMesh: - pcd = filter_pcd(pcd) - if len(pcd.points) <= 2: - return o3d.geometry.TriangleMesh() - - _, _, max_r = compute_local_densities(pcd) - radii = [max_r/2, max_r] - mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_ball_pivoting( - pcd, o3d.utility.DoubleVector(radii) - ) - return mesh - -def poisson_triangulation(pcd: o3d.geometry.PointCloud, surface = None) -> o3d.geometry.TriangleMesh: - pcd = filter_pcd(pcd) - if len(pcd.points) <= 2: - return o3d.geometry.TriangleMesh() - - mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson( - pcd, depth=9 - ) - return mesh - -def triangulation(pcd: o3d.geometry.PointCloud, data: dict) -> o3d.geometry.TriangleMesh: - pcd_surf = pcd.select_by_index(data["point_indices"]) - - mesh = o3d.geometry.TriangleMesh() - - if len(pcd_surf.points) > 0: - result = SurfacesProjector.projectPointsOnSurfaceFeatures(np.asarray(pcd_surf.points), data) - if result is not None: - uvs, points = result - data["point_parameters"] = uvs - if PROJECT_POINTS: - pcd_surf.points = o3d.utility.Vector3dVector(points) - mesh = bpa_triangulation(pcd_surf, surface=data) - - return mesh - -def triangulation_by_surface(pcd: o3d.geometry.PointCloud, surfaces_data: dict) -> o3d.geometry.TriangleMesh: - func = partial(triangulation, pcd) - result_ = thread_map(func, surfaces_data, max_workers=MAX_WORKERS, chunk_size=1) - - final_mesh = o3d.geometry.TriangleMesh() - for result_idx, result in enumerate(result_): - surfaces_data[result_idx]["vert_indices"] = list(range(len(final_mesh.vertices), - len(final_mesh.vertices)+len(result.vertices))) - del surfaces_data[result_idx]["vert_parameters"] - res = SurfacesProjector.projectPointsOnSurfaceFeatures(np.asarray(result.vertices), - surfaces_data[result_idx]) - if res is not None: - uvs, _ = res - surfaces_data[result_idx]["vert_parameters"] = uvs if type(uvs) is list else uvs.tolist() - surfaces_data[result_idx]["face_indices"] = list(range(len(final_mesh.triangles), - len(final_mesh.trangles)+len(result.triangles))) - surfaces_data[result_idx]["point_indices"] = surfaces_data[result_idx]["point_indices"] \ - if type(surfaces_data[result_idx]["point_indices"]) is list \ - else surfaces_data[i]['point_indices'].tolist() - if "point_parameters" in surfaces_data[result_idx]: - surfaces_data[result_idx]["point_parameters"] = surfaces_data[result_idx]["point_parameters"] \ - if type(surfaces_data[result_idx]["point_parameters"]) is list \ - else surfaces_data[result_idx]["point_parameters"].tolist() - - final_mesh += result - final_mesh.compute_triangle_normals() - return final_mesh - -def generate_mesh_from_pointcloud(args, pc: pypcd.PointCloud, pc_filename: str, geometry_data: dict) -> bool: - filename = pc_filename - foldername = args['mesh_foldername'] - triangulation_foldername = args['triangulation_foldername'] - triangulation_features_foldername = args['triangulation_features_foldername'] - GENERATE_TRIANGULATION_FEATURES = args['generate_triangulation_features'] - PROJECT_POINTS = args['project_points'] - N_NEIGHBORS = args['n_neighbors_mesh'] - LEAF_SIZE = args['leaf_size_mesh'] - MAX_WORKERS = args['max_workers_mesh'] - FILTER = not args['no_filter'] - reegenerate = args['reegenerate'] - - print() - print("Generating mesh for the file: {}".format(filename)) - - pointcloud = o3d.geometry.PointCloud() - - if 'x' in pc.dtype.names and 'y' in pc.dtype.names and 'z' in pc.dtype.names: - pointcloud.points = o3d.utility.Vector3dVector(np.vstack((pc['x'], pc['y'], pc['z'])).T) - else: - print("There is no point in the input pointcloud") - return False - - if 'normal_x' in pc.dtype.names and 'normal_y' in pc.dtype.names and 'normal_z' in pc.dtype.names: - pointcloud.normals = o3d.utility.Vector3dVector(np.vstack((pc['normal_x'], pc['normal_y'], pc['normal_z'])).T) - else: - print('There is no normal in the input pointcloud') - - has_labels = False - if 'label' in pc.dtype.names: - labels = pc['label'] - has_labels = True - else: - print("there is no label in the input pointcloud") - - surfaces_data = geometry_data['surfaces'] - - if has_labels: - labels, fpi = compute_labels_from_face_2_primitive(labels, surfaces_data) - for i in range(len(fpi)): - surfaces_data[i]['point_indices'] = fpi[i] - - if not pointcloud.has_normals(): - pointcloud.estimate_normals() - - mesh = triangulation_by_surface(pointcloud, surfaces_data) - igl.write_triangle_mesh(os.path.join(folder_name, triangulation_foldername, filename + '.obj'), np.asarray(mesh.vertices), np.asarray(mesh.triangles)) - - if GENERATE_TRIANGULATION_FEATURES: - features = {'surfaces': surfaces_data} - with open(os.path.join(folder_name, triangulation_features_foldername, filename + '.json'), 'w') as f: - json.dump(features, f, indent=4) - - return True diff --git a/mesh_generator.py b/mesh_generator.py new file mode 100644 index 0000000..cae9b51 --- /dev/null +++ b/mesh_generator.py @@ -0,0 +1,311 @@ +import argparse +from os.path import join +from pypcd import pypcd +import numpy as np +from scipy.spatial import Delaunay + +import matplotlib.pyplot as plt + +from scipy.spatial import Delaunay + +import open3d as o3d +import json + +from tqdm.contrib.concurrent import thread_map + +from functools import partial + +from lib.surfaces_projector import SurfacesProjector + +import os +from pathlib import Path + +import igl + +import json + +EPS = np.finfo(float).eps +LEAF_SIZE = 0.02 +N_NEIGHBORS = 10 +MAX_WORKERS = 8 +GENERATE_TRIANGULATION_FEATURES = False +PROJECT_POINTS = False +FILTER = True + +def computeLabelsFromFace2Primitive(labels, features_data): + max_face = np.max(labels) + for feat in features_data: + max_face = max(0 if len(feat['face_indices']) == 0 else max(feat['face_indices']), max_face) + face_2_primitive = np.zeros(shape=(max_face+1,), dtype=np.int32) - 1 + face_primitive_count = np.zeros(shape=(max_face+1,), dtype=np.int32) + for i, feat in enumerate(features_data): + for face in feat['face_indices']: + face_2_primitive[face] = i + face_primitive_count[face] += 1 + assert len(np.unique(face_primitive_count)) <= 2 + features_point_indices = [[] for i in range(0, len(features_data) + 1)] + for i in range(0, len(labels)): + index = face_2_primitive[labels[i]] + features_point_indices[index].append(i) + labels[i] = index + features_point_indices.pop(-1) + + for i in range(0, len(features_point_indices)): + features_point_indices[i] = np.array(features_point_indices[i], dtype=np.int64) + + return labels, features_point_indices + +def compute_local_densities(pcd, k=N_NEIGHBORS): + if len(pcd.points) > 1: + pcd_tree = o3d.geometry.KDTreeFlann(pcd) + points = np.asarray(pcd.points) + distances = [] + for i in range(len(points)): + A = points[i] + _, idx, _ = pcd_tree.search_knn_vector_3d(A, k) + points_query = points[np.asarray(idx)] + B = points_query[-1] + dist = np.linalg.norm(B-A, ord=2) + distances.append(dist) + dist_arr = np.array(distances) + return np.min(dist_arr, axis=0), np.mean(dist_arr, axis=0), np.max(dist_arr, axis=0) + else: + return 0, 0, 0 + +def filter_pcd(pcd): + if FILTER: + pcd = pcd.voxel_down_sample(voxel_size=LEAF_SIZE) + pcd, _ = pcd.remove_statistical_outlier(nb_neighbors=N_NEIGHBORS, std_ratio=2.0) + return pcd + +def delaunay_uv_triangulation(pcd, surface): + pcd = filter_pcd(pcd) + if len(pcd.points) > 2: + _, _, max_r = compute_local_densities(pcd) + + points = np.asarray(pcd.points) + + response = SurfacesProjector.projectPointsOnSurfaceFeatures(points, surface) + + if response is None: + return o3d.geometry.TriangleMesh() + + uvs, _ = response + + tri = Delaunay(uvs) + triangles = tri.simplices + + triangles_filter = [] + index_perm = [(0,1), (0,2), (1,2)] + for t in triangles: + keep = True + for i, j in index_perm: + A = points[t[i]] + B = points[t[j]] + dist = np.linalg.norm(B - A, ord=2) + if dist > max_r: + keep = False + break + triangles_filter.append(keep) + triangles_filter = np.array(triangles_filter) + + triangles = triangles[triangles_filter] + + mesh = o3d.geometry.TriangleMesh() + mesh.vertices = o3d.utility.Vector3dVector(points) + mesh.triangles = o3d.utility.Vector3iVector(triangles) + + plt.scatter(uvs[:, 0], uvs[:, 1]) + plt.show() + + return mesh + else: + return o3d.geometry.TriangleMesh() + +def bpa_triangulation(pcd, surface=None): + pcd = filter_pcd(pcd) + if len(pcd.points) > 2: + _, _, max_r = compute_local_densities(pcd) + radii = [max_r/2, max_r] + mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_ball_pivoting( + pcd, o3d.utility.DoubleVector(radii)) + return mesh + else: + return o3d.geometry.TriangleMesh() + +def poisson_triangulation(pcd, surface=None): + pcd = filter_pcd(pcd) + if len(pcd.points) > 2: + mesh, _ = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson( + pcd, depth=9) + return mesh + else: + return o3d.geometry.TriangleMesh() + +def triangulation(pcd, data): + pcd_surf = pcd.select_by_index(data['point_indices']) + + mesh = o3d.geometry.TriangleMesh() + + if len(pcd_surf.points) > 0: + result = SurfacesProjector.projectPointsOnSurfaceFeatures(np.asarray(pcd_surf.points), data) + if result is not None: + uvs, points = result + data['point_parameters'] = uvs + if PROJECT_POINTS: + pcd_surf.points = o3d.utility.Vector3dVector(points) + + mesh = bpa_triangulation(pcd_surf, surface=data) + + return mesh + +def triangulation_by_surface(pcd, surfaces_data): + func = partial(triangulation, pcd) + result = thread_map(func, surfaces_data, max_workers=MAX_WORKERS, chunksize=1) + + final_mesh = o3d.geometry.TriangleMesh() + for i, r in enumerate(result): + surfaces_data[i]['vert_indices'] = list(range(len(final_mesh.vertices), len(final_mesh.vertices) + len(r.vertices))) + del surfaces_data[i]['vert_parameters'] #missing + res = SurfacesProjector.projectPointsOnSurfaceFeatures(np.asarray(r.vertices), surfaces_data[i]) + if res is not None: + uvs, _ = res + surfaces_data[i]['vert_parameters'] = uvs if type(uvs) is list else uvs.tolist() + surfaces_data[i]['face_indices'] = list(range(len(final_mesh.triangles), len(final_mesh.triangles) + len(r.triangles))) + surfaces_data[i]['point_indices'] = surfaces_data[i]['point_indices'] if type(surfaces_data[i]['point_indices']) is list else surfaces_data[i]['point_indices'].tolist() + if 'point_parameters' in surfaces_data[i]: + surfaces_data[i]['point_parameters'] = surfaces_data[i]['point_parameters'] if type(surfaces_data[i]['point_parameters']) is list else surfaces_data[i]['point_parameters'].tolist() + + final_mesh += r + + final_mesh.compute_triangle_normals() + return final_mesh + +def list_files(input_dir: str, return_str=False) -> list: + files = [] + path = Path(input_dir) + for file_path in path.glob('*'): + files.append(file_path if not return_str else str(file_path)) + return sorted(files) + +def list_filenames(input_dir): + files = list_files(input_dir, return_str=True) + filenames = [] + for f in files: + start_index = f.rfind('/') + start_index = start_index if start_index != -1 else 0 + final_index = f.rfind('.') + filenames.append(f[(start_index+1):final_index]) + return filenames + + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description='') + parser.add_argument('folder', type=str, help='') + parser.add_argument('-gf','--generate_triangulation_features', action='store_true', help='') + parser.add_argument('-pp','--project_points', action='store_true', help='') + parser.add_argument('-r','--reegenerate', action='store_true', help='') + parser.add_argument('-nf','--no_filter', action='store_true', help='') + parser.add_argument('--filename', type=str, default='', help='') + parser.add_argument('--pc_foldername', type=str, default='pc', help='') + parser.add_argument('--triangulation_foldername', type=str, default='triangulation', help='') + parser.add_argument('--triangulation_features_foldername', type=str, default='triangulation_features', help='') + parser.add_argument('--n_neighbors', type=int, default=N_NEIGHBORS, help='') + parser.add_argument('--leaf_size', type=float, default=LEAF_SIZE, help='') + parser.add_argument('--max_workers', type=int, default=MAX_WORKERS, help='') + + + args = vars(parser.parse_args()) + + folder_name = args['folder'] + filename = args['filename'] + pc_foldername = args['pc_foldername'] + triangulation_foldername = args['triangulation_foldername'] + triangulation_features_foldername = args['triangulation_features_foldername'] + GENERATE_TRIANGULATION_FEATURES = args['generate_triangulation_features'] + PROJECT_POINTS = args['project_points'] + N_NEIGHBORS = args['n_neighbors'] + LEAF_SIZE = args['leaf_size'] + MAX_WORKERS = args['max_workers'] + FILTER = not args['no_filter'] + reegenerate = args['reegenerate'] + + if not os.path.exists(folder_name) or not os.path.isdir(folder_name): + print('[Generator Error] Input path not found') + exit() + + filenames = [] + if filename != '': + filenames.append(filename) + else: + filenames = list_filenames(os.path.join(folder_name, 'features')) + + if not reegenerate: + filenames_triang = list_filenames(os.path.join(folder_name, triangulation_foldername)) + filenames_triang_features = list_filenames(os.path.join(folder_name, triangulation_features_foldername)) + + + diff_1 = set(filenames) - set(filenames_triang) + diff_2 = set(filenames) - set(filenames_triang_features) + + if GENERATE_TRIANGULATION_FEATURES: + diff_1 = diff_1.union(diff_2) + + filenames = list(diff_1) + + os.makedirs(os.path.join(folder_name, triangulation_foldername), exist_ok=True) + + if GENERATE_TRIANGULATION_FEATURES: + os.makedirs(os.path.join(folder_name, triangulation_features_foldername), exist_ok=True) + + for filename in filenames: + print() + print('Processing file {}:'.format(filename)) + + pc_filename = join(folder_name, 'pc', filename + '.pcd') + + geometry_filename = join(folder_name, 'features', filename + '.json') + + pc = pypcd.PointCloud.from_path(pc_filename).pc_data + + pointcloud = o3d.geometry.PointCloud() + + if 'x' in pc.dtype.names and 'y' in pc.dtype.names and 'z' in pc.dtype.names: + pointcloud.points = o3d.utility.Vector3dVector(np.vstack((pc['x'], pc['y'], pc['z'])).T) + else: + print('there is no point in the input pointcloud') + exit() + + if 'normal_x' in pc.dtype.names and 'normal_y' in pc.dtype.names and 'normal_z' in pc.dtype.names: + pointcloud.normals = o3d.utility.Vector3dVector(np.vstack((pc['normal_x'], pc['normal_y'], pc['normal_z'])).T) + else: + print('there is no normal in the input pointcloud') + + has_labels = False + + if 'label' in pc.dtype.names: + labels = pc['label'] + has_labels = True + else: + print('there is no label in the input pointcloud') + + + with open(geometry_filename, 'r') as f: + geometry_data = json.load(f) + surfaces_data = geometry_data['surfaces'] + + if has_labels: + labels, fpi = computeLabelsFromFace2Primitive(labels, surfaces_data) + for i in range(len(fpi)): + surfaces_data[i]['point_indices'] = fpi[i] + + if not pointcloud.has_normals(): + pointcloud.estimate_normals() + + mesh = triangulation_by_surface(pointcloud, surfaces_data) + igl.write_triangle_mesh(os.path.join(folder_name, triangulation_foldername, filename + '.obj'), np.asarray(mesh.vertices), np.asarray(mesh.triangles)) + + if GENERATE_TRIANGULATION_FEATURES: + features = {'surfaces': surfaces_data} + with open(os.path.join(folder_name, triangulation_features_foldername, filename + '.json'), 'w') as f: + json.dump(features, f, indent=4) \ No newline at end of file From 79c910bf76b3631885c5609e0cbda56ab7c8a6bb Mon Sep 17 00:00:00 2001 From: PedroCorcaque Date: Tue, 11 Jun 2024 19:37:30 -0300 Subject: [PATCH 4/5] rename the executable file to 'triangulator' --- mesh_generator.py => triangulator.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename mesh_generator.py => triangulator.py (100%) diff --git a/mesh_generator.py b/triangulator.py similarity index 100% rename from mesh_generator.py rename to triangulator.py From 0430ffa83c9fdb298eb7c6e9d8e08be1dfd420a5 Mon Sep 17 00:00:00 2001 From: PedroCorcaque Date: Mon, 24 Jun 2024 11:14:14 -0300 Subject: [PATCH 5/5] Removed non-used method; Change the surface projector method to use from AsGeometryOCCWrapper --- dataset_generator.py | 23 --------- lib/surfaces_projector.py | 104 -------------------------------------- triangulator.py | 66 +++++------------------- 3 files changed, 12 insertions(+), 181 deletions(-) delete mode 100644 lib/surfaces_projector.py diff --git a/dataset_generator.py b/dataset_generator.py index 5778e05..5eb30b3 100644 --- a/dataset_generator.py +++ b/dataset_generator.py @@ -13,7 +13,6 @@ from lib.utils import loadFeatures, computeLabelsFromFace2Primitive, savePCD, downsampleByPointIndices, rayCastingPointCloudGeneration, funif, computeFeaturesPointIndices from lib.writers import DatasetWriterFactory -from lib import triangulation from asGeometryOCCWrapper.surfaces import SurfaceFactory @@ -59,22 +58,8 @@ parser.add_argument('-m_np', '--min_number_points', type=float, default = 0.0001, help='filter geometries by number of points.') parser.add_argument('-ls', '--leaf_size', type=float, default = 0.0, help='') - parser.add_argument('--generate_mesh', action='store_true', default=False, help='Generate mesh from point cloud') - parser.add_argument('--mesh_foldername', type=str, default='mesh', help='') - parser.add_argument('--triangulation_foldername', type=str, default='triangulation', help='') - parser.add_argument('--triangulation_features_foldername', type=str, default='triangulation_features', help='') - parser.add_argument('-gf', "--generate_triangulation_features", action="store_true", help="[Optional] Used in mesh generation.") - parser.add_argument('-pp', "--project_points", action="store_true", help="[Optional] Used in mesh generation.") - parser.add_argument('-r', "--reegenerate", action="store_true", help="[Optional] Used in mesh generation.") - parser.add_argument('-nf', "--no_filter", action="store_true", help="[Optional] Used in mesh generation.") - parser.add_argument('--n_neighbors_mesh', type=int, default=-1, help="[Optional] Used in mesh generation.") - parser.add_argument('--leaf_size_mesh', type=float, default=-1, help="[Optional] Used in mesh generation.") - parser.add_argument('--max_workers_mesh', type=int, default=-1, help="[Optional] Used in mesh generation.") - args = vars(parser.parse_args()) - generate_mesh = args["generate_mesh"] - folder_name = args['folder'] formats = [s.lower() for s in args['formats'].split(',')] curve_types = [s.lower() for s in args['curve_types'].split(',')] @@ -190,14 +175,6 @@ labels, features_point_indices = computeLabelsFromFace2Primitive(labels_mesh.copy(), features_data['surfaces']) print('Done.\n') - if generate_mesh: - print("Generating mesh from point cloud...") - if (generate_mesh_from_pointcloud(args, pc_filename, features_data)): - print("\tMesh generated on: {}".format(mesh_folder_name)) - else: - print("\tThe mesh could not be generated.") - print("Done") - elif mesh_filename is not None: print('Opening Mesh...') mesh = o3d.io.read_triangle_mesh(mesh_filename, print_progress=True) diff --git a/lib/surfaces_projector.py b/lib/surfaces_projector.py deleted file mode 100644 index a92082c..0000000 --- a/lib/surfaces_projector.py +++ /dev/null @@ -1,104 +0,0 @@ -from OCC.Core.Geom import Geom_CylindricalSurface, Geom_ConicalSurface, Geom_Plane, Geom_SphericalSurface, Geom_BSplineSurface, Geom_SurfaceOfLinearExtrusion, Geom_SurfaceOfRevolution, Geom_ToroidalSurface -from OCC.Core.GeomAPI import GeomAPI_ProjectPointOnSurf -from OCC.Core.gp import gp_Ax3, gp_Pnt, gp_Dir -from OCC.Core.TColgp import TColgp_Array2OfPnt -from OCC.Core.TColStd import TColStd_Array2OfReal, TColStd_Array1OfReal, TColStd_Array1OfInteger - -import numpy as np - -class SurfacesProjector: - PROJECTOR = GeomAPI_ProjectPointOnSurf() - - def buildCylinder(features): - radius = features['radius'] - location = features['location'] - z_axis = features['z_axis'] - axis = gp_Ax3(gp_Pnt(location[0], location[1], location[2]), gp_Dir(z_axis[0], z_axis[1], z_axis[2])) - surf = Geom_CylindricalSurface(axis, radius) - return surf - - def buildBSpline(features): - poles = np.asarray(features['poles']) - poles_tcol = TColgp_Array2OfPnt(1, poles.shape[0], 1, poles.shape[1]) - for i in range(poles.shape[0]): - for j in range(poles.shape[1]): - print(i, j) - e = poles[i, j, :] - #poles_tcol.SetValue(i, j, gp_Pnt(e[0], e[1], e[2]) ) - - weights = features['weights'] - u_knots = features['u_knots'] - v_knots = features['v_knots'] - u_degree = features['u_degree'] - v_degree = features['v_degree'] - surf = Geom_BSplineSurface(poles_tcol, TColStd_Array2OfReal(weights), TColStd_Array1OfReal(u_knots), TColStd_Array1OfReal(v_knots), TColStd_Array1OfInteger([1 for i in range(len(u_knots))]), TColStd_Array1OfInteger([1 for i in range(len(v_knots))]), u_degree, v_degree) - return surf - - def buildCone(features): - radius = features['radius'] - location = features['location'] - z_axis = features['z_axis'] - angle = features['angle'] - axis = gp_Ax3(gp_Pnt(location[0], location[1], location[2]), gp_Dir(z_axis[0], z_axis[1], z_axis[2])) - surf = Geom_ConicalSurface(axis, angle, radius) - return surf - - def buildPlane(features): - location = features['location'] - normal = features['z_axis'] - axis = gp_Ax3(gp_Pnt(location[0], location[1], location[2]), gp_Dir(normal[0], normal[1], normal[2])) - surf = Geom_Plane(axis) - return surf - - def buildSphere(features): - radius = features['radius'] - location = features['location'] - z_axis = features['z_axis'] - axis = gp_Ax3(gp_Pnt(location[0], location[1], location[2]), gp_Dir(z_axis[0], z_axis[1], z_axis[2])) - surf = Geom_SphericalSurface(axis, radius) - return surf - - def buildSurfaceRevolution(features): - pass - - def buildTorus(features): - min_radius = features['min_radius'] - max_radius = features['max_radius'] - location = features['location'] - z_axis = features['z_axis'] - axis = gp_Ax3(gp_Pnt(location[0], location[1], location[2]), gp_Dir(z_axis[0], z_axis[1], z_axis[2])) - surf = Geom_ToroidalSurface(axis, max_radius, min_radius) - return surf - - def buildSurfaceExtrusion(features): - pass - - BUILD_FUNCTION_MAP = { - 'cylinder': buildCylinder, - 'cone': buildCone, - 'plane': buildPlane, - 'sphere': buildSphere, - 'torus': buildTorus - } - - @staticmethod - def projectPointsOnSurfaceFeatures(points, features): - tp = features['type'].lower() - if tp in SurfacesProjector.BUILD_FUNCTION_MAP.keys(): - if len(points) == 0: - return [], [] - surface = SurfacesProjector.BUILD_FUNCTION_MAP[tp](features) - uvs = np.zeros((points.shape[0], 2), dtype=np.float64) - points_projected = np.zeros(points.shape, dtype=np.float64) - SurfacesProjector.PROJECTOR.Init(gp_Pnt(points[0, 0], points[0, 1], points[0, 2]), surface) - uvs[0] = np.array(SurfacesProjector.PROJECTOR.LowerDistanceParameters()) - pnt = SurfacesProjector.PROJECTOR.NearestPoint() - points_projected[0] = np.array([pnt.X(), pnt.Y(), pnt.Z()]) - for idx, pt in enumerate(points[1:, :]): - SurfacesProjector.PROJECTOR.Perform(gp_Pnt(pt[0], pt[1], pt[2])) - uvs[idx+1] = np.array(SurfacesProjector.PROJECTOR.LowerDistanceParameters()) - pnt = SurfacesProjector.PROJECTOR.NearestPoint() - points_projected[idx+1] = np.array([pnt.X(), pnt.Y(), pnt.Z()]) - return uvs, points_projected - else: - return None \ No newline at end of file diff --git a/triangulator.py b/triangulator.py index cae9b51..2f41bc8 100644 --- a/triangulator.py +++ b/triangulator.py @@ -15,7 +15,7 @@ from functools import partial -from lib.surfaces_projector import SurfacesProjector +from asGeometryOCCWrapper.surfaces import SurfaceFactory import os from pathlib import Path @@ -78,50 +78,6 @@ def filter_pcd(pcd): pcd, _ = pcd.remove_statistical_outlier(nb_neighbors=N_NEIGHBORS, std_ratio=2.0) return pcd -def delaunay_uv_triangulation(pcd, surface): - pcd = filter_pcd(pcd) - if len(pcd.points) > 2: - _, _, max_r = compute_local_densities(pcd) - - points = np.asarray(pcd.points) - - response = SurfacesProjector.projectPointsOnSurfaceFeatures(points, surface) - - if response is None: - return o3d.geometry.TriangleMesh() - - uvs, _ = response - - tri = Delaunay(uvs) - triangles = tri.simplices - - triangles_filter = [] - index_perm = [(0,1), (0,2), (1,2)] - for t in triangles: - keep = True - for i, j in index_perm: - A = points[t[i]] - B = points[t[j]] - dist = np.linalg.norm(B - A, ord=2) - if dist > max_r: - keep = False - break - triangles_filter.append(keep) - triangles_filter = np.array(triangles_filter) - - triangles = triangles[triangles_filter] - - mesh = o3d.geometry.TriangleMesh() - mesh.vertices = o3d.utility.Vector3dVector(points) - mesh.triangles = o3d.utility.Vector3iVector(triangles) - - plt.scatter(uvs[:, 0], uvs[:, 1]) - plt.show() - - return mesh - else: - return o3d.geometry.TriangleMesh() - def bpa_triangulation(pcd, surface=None): pcd = filter_pcd(pcd) if len(pcd.points) > 2: @@ -148,12 +104,12 @@ def triangulation(pcd, data): mesh = o3d.geometry.TriangleMesh() if len(pcd_surf.points) > 0: - result = SurfacesProjector.projectPointsOnSurfaceFeatures(np.asarray(pcd_surf.points), data) - if result is not None: - uvs, points = result - data['point_parameters'] = uvs + surface = SurfaceFactory.fromDict(data) + if surface is not None: + proj_points, _, proj_params = surface.projectPointsOnGeometry(np.asarray(pcd_surf.points)) + data['point_parameters'] = proj_params if PROJECT_POINTS: - pcd_surf.points = o3d.utility.Vector3dVector(points) + pcd_surf.points = o3d.utility.Vector3dVector(proj_points) mesh = bpa_triangulation(pcd_surf, surface=data) @@ -167,10 +123,12 @@ def triangulation_by_surface(pcd, surfaces_data): for i, r in enumerate(result): surfaces_data[i]['vert_indices'] = list(range(len(final_mesh.vertices), len(final_mesh.vertices) + len(r.vertices))) del surfaces_data[i]['vert_parameters'] #missing - res = SurfacesProjector.projectPointsOnSurfaceFeatures(np.asarray(r.vertices), surfaces_data[i]) - if res is not None: - uvs, _ = res - surfaces_data[i]['vert_parameters'] = uvs if type(uvs) is list else uvs.tolist() + + surface = SurfaceFactory.fromDict(surfaces_data[i]) + if surface is not None: + _, _, proj_params = surface.projectPointsOnGeometry(np.asarray(r.vertices)) + surfaces_data[i]['vert_parameters'] = proj_params + surfaces_data[i]['face_indices'] = list(range(len(final_mesh.triangles), len(final_mesh.triangles) + len(r.triangles))) surfaces_data[i]['point_indices'] = surfaces_data[i]['point_indices'] if type(surfaces_data[i]['point_indices']) is list else surfaces_data[i]['point_indices'].tolist() if 'point_parameters' in surfaces_data[i]: