133 lines
4.6 KiB
Python
133 lines
4.6 KiB
Python
import math
|
|
import numpy as np
|
|
from trimesh import Trimesh
|
|
from trimesh import intersections
|
|
import open3d as o3d
|
|
|
|
|
|
class Geometry:
|
|
def __init__(self, delta=0.5):
|
|
self._delta = delta
|
|
|
|
def almost_equal(self, v1, v2):
|
|
delta = math.sqrt(pow((v1[0] - v2[0]), 2) + pow((v1[1] - v2[1]), 2) + pow((v1[2] - v2[2]), 2))
|
|
return delta <= self._delta
|
|
|
|
def is_almost_same_surface(self, s1, s2):
|
|
# delta is grads an need to be converted into radians
|
|
delta = np.rad2deg(self._delta)
|
|
difference = (s1.inclination - s2.inclination) % math.pi
|
|
if abs(difference) > delta:
|
|
return False
|
|
# s1 and s2 are at least almost parallel surfaces
|
|
# calculate distance point to plane using all the vertex
|
|
# select surface1 value for the point (X,Y,Z) where two of the values are 0
|
|
minimum_distance = self._delta + 1
|
|
parametric = s2.polygon.get_parametric()
|
|
n2 = s2.normal
|
|
for point in s1.points:
|
|
distance = abs(
|
|
(point[0] * parametric[0]) + (point[1] * parametric[1]) + (point[2] * parametric[2]) + parametric[3])
|
|
normal_module = math.sqrt(pow(n2[0], 2) + pow(n2[1], 2) + pow(n2[2], 2))
|
|
|
|
if normal_module == 0:
|
|
continue
|
|
distance = distance / normal_module
|
|
if distance < minimum_distance:
|
|
minimum_distance = distance
|
|
if minimum_distance <= self._delta:
|
|
break
|
|
if minimum_distance > self._delta or s1.intersect(s2) is None:
|
|
return False
|
|
else:
|
|
return True
|
|
|
|
@staticmethod
|
|
def to_points_matrix(points, remove_last=False):
|
|
rows = points.size // 3
|
|
points = points.reshape(rows, 3)
|
|
if remove_last:
|
|
points = np.delete(points, rows - 1, 0)
|
|
return points
|
|
|
|
@staticmethod
|
|
def _segment_list_to_point_cloud(segment_list):
|
|
point_list = np.asarray(segment_list[0])
|
|
for segment in segment_list:
|
|
for new_point in segment:
|
|
found = False
|
|
for point in point_list:
|
|
same_point = np.allclose(new_point, point)
|
|
if same_point:
|
|
found = True
|
|
break
|
|
if not found:
|
|
point_list = np.concatenate((point_list, [new_point]))
|
|
return point_list
|
|
|
|
@staticmethod
|
|
def _point_cloud_to_mesh(point_list, normal_list):
|
|
# Return a mesh composed only by triangles
|
|
pcd = o3d.geometry.PointCloud()
|
|
pcd.points = o3d.utility.Vector3dVector(point_list)
|
|
pcd.normals = o3d.utility.Vector3dVector(normal_list)
|
|
distances = pcd.compute_nearest_neighbor_distance()
|
|
avg_dist = np.mean(distances)
|
|
radius = 3 * avg_dist
|
|
bpa_mesh = o3d.geometry.TriangleMesh().create_from_point_cloud_ball_pivoting(
|
|
pcd, o3d.utility.DoubleVector([radius, radius * 2]))
|
|
mesh_result = Trimesh(vertices=np.asarray(bpa_mesh.vertices), faces=np.asarray(bpa_mesh.triangles))
|
|
return mesh_result
|
|
|
|
@staticmethod
|
|
def _merge_meshes(mesh1, mesh2):
|
|
v_1 = mesh1.vertices
|
|
f_1 = mesh1.faces
|
|
v_2 = mesh2.vertices
|
|
f_2 = mesh2.faces
|
|
length = len(v_1)
|
|
v_merge = np.concatenate((v_1, v_2))
|
|
f_merge = np.asarray(f_1)
|
|
|
|
for item in f_2:
|
|
point1 = item.item(0) + length
|
|
point2 = item.item(1) + length
|
|
point3 = item.item(2) + length
|
|
surface = np.asarray([point1, point2, point3])
|
|
f_merge = np.concatenate((f_merge, [surface]))
|
|
|
|
mesh_merge = Trimesh(vertices=v_merge, faces=f_merge)
|
|
|
|
return mesh_merge
|
|
|
|
@staticmethod
|
|
def divide_mesh_by_plane(mesh, normal_plane, point_plane):
|
|
# The first mesh returns the positive side of the plane and the second the negative side.
|
|
# If the plane does not divide the mesh (i.e. it does not touch it or it is coplanar with one or more faces),
|
|
# then it returns only the original mesh.
|
|
normal_plane_opp = [None] * len(normal_plane)
|
|
for i in range(0, len(normal_plane)):
|
|
normal_plane_opp[i] = - normal_plane[i]
|
|
|
|
normal = [normal_plane, normal_plane_opp]
|
|
normal_opp = [normal_plane_opp, normal_plane]
|
|
mesh_final = []
|
|
for i in range(0, 2):
|
|
mesh_1 = intersections.slice_mesh_plane(mesh, normal[i], point_plane)
|
|
mesh_1_segments = intersections.mesh_plane(mesh, normal[i], point_plane)
|
|
boo = mesh.difference(mesh_1, engine='blender')
|
|
print(boo)
|
|
quit()
|
|
if len(mesh_1_segments) <= 0 or len(mesh_1.faces) == len(mesh.faces):
|
|
mesh_final.append(mesh)
|
|
break
|
|
else:
|
|
points = Geometry._segment_list_to_point_cloud(mesh_1_segments)
|
|
points_normals = [[None] * 3] * len(points)
|
|
for j in range(0, len(points_normals)):
|
|
points_normals[j] = normal_opp[i]
|
|
mesh_2 = Geometry._point_cloud_to_mesh(points, points_normals)
|
|
mesh_final.append(Geometry._merge_meshes(mesh_1, mesh_2))
|
|
|
|
return mesh_final
|