481 lines
16 KiB
Python
481 lines
16 KiB
Python
"""
|
|
Polyhedron module
|
|
SPDX - License - Identifier: LGPL - 3.0 - or -later
|
|
Copyright © 2020 Project Author Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
|
|
Contributors Pilar Monsalvete pilar_monsalvete@yahoo.es
|
|
"""
|
|
import sys
|
|
import numpy as np
|
|
from trimesh import Trimesh
|
|
from helpers.geometry_helper import GeometryHelper
|
|
from helpers.configuration_helper import ConfigurationHelper
|
|
from city_model_structure.attributes.polygon import Polygon
|
|
from helpers.geometry_helper import GeometryHelper as gh
|
|
|
|
|
|
class Polyhedron:
|
|
"""
|
|
Polyhedron class
|
|
"""
|
|
|
|
def __init__(self, polygons):
|
|
self._polygons = polygons
|
|
self._polyhedron = None
|
|
self._triangulated_polyhedron = None
|
|
self._volume = None
|
|
self._faces = None
|
|
self._vertices = None
|
|
self._mesh = None
|
|
self._centroid = None
|
|
self._max_z = None
|
|
self._max_y = None
|
|
self._max_x = None
|
|
self._min_z = None
|
|
self._min_y = None
|
|
self._min_x = None
|
|
self._geometry = GeometryHelper()
|
|
|
|
def _position_of(self, point, face):
|
|
vertices = self.vertices
|
|
for i in range(len(vertices)):
|
|
# ensure not duplicated vertex
|
|
if i not in face and GeometryHelper.distance_between_points(vertices[i], point) == 0:
|
|
return i
|
|
return -1
|
|
|
|
@property
|
|
def vertices(self) -> np.ndarray:
|
|
"""
|
|
Polyhedron vertices
|
|
:return: np.ndarray(int)
|
|
"""
|
|
if self._vertices is None:
|
|
vertices, self._vertices = [], []
|
|
_ = [vertices.extend(s.points) for s in self._polygons]
|
|
for vertex_1 in vertices:
|
|
found = False
|
|
for vertex_2 in self._vertices:
|
|
found = False
|
|
if GeometryHelper.distance_between_points(vertex_1, vertex_2) == 0:
|
|
found = True
|
|
break
|
|
if not found:
|
|
self._vertices.append(vertex_1)
|
|
self._vertices = np.asarray(self._vertices)
|
|
return self._vertices
|
|
|
|
def _triangulate(self, polygon) -> [Polygon]:
|
|
"""
|
|
triangulates a polygon following the ear clipping methodology
|
|
:param polygon: polygon
|
|
:return: list[triangles]
|
|
"""
|
|
points_list = polygon.points_list
|
|
normal = polygon.normal
|
|
# are points concave or convex?
|
|
total_points_list, concave_points, convex_points = self._starting_lists(points_list, normal)
|
|
|
|
# list of ears
|
|
ears = []
|
|
while len(concave_points) > 3 or len(convex_points) != 0:
|
|
for i in range(0, len(concave_points)):
|
|
ear = self._triangle(points_list, total_points_list, concave_points[i])
|
|
rest_points = []
|
|
for p in total_points_list:
|
|
rest_points.append(list(polygon.points[p]))
|
|
if self._is_ear(ear, rest_points):
|
|
ears.append(ear)
|
|
point_to_remove = concave_points[i]
|
|
previous_point_in_list, next_point_in_list = self._enveloping_points(point_to_remove, total_points_list)
|
|
total_points_list.remove(point_to_remove)
|
|
concave_points.remove(point_to_remove)
|
|
# Was any of the adjacent points convex? -> check if changed status to concave
|
|
# for j in range(0, len(convex_points)):
|
|
for convex_point in convex_points:
|
|
if convex_point == previous_point_in_list:
|
|
concave_points, convex_points, end_loop = self._if_concave_change_status(normal, points_list,
|
|
convex_point, total_points_list,
|
|
concave_points, convex_points,
|
|
previous_point_in_list)
|
|
if end_loop:
|
|
break
|
|
continue
|
|
if convex_point == next_point_in_list:
|
|
concave_points, convex_points, end_loop = self._if_concave_change_status(normal, points_list,
|
|
convex_point, total_points_list,
|
|
concave_points, convex_points,
|
|
next_point_in_list)
|
|
if end_loop:
|
|
break
|
|
continue
|
|
break
|
|
if len(total_points_list) <= 3 and len(convex_points) > 0:
|
|
sys.stderr.write(f'Not able to triangulate polygon\n')
|
|
return [polygon]
|
|
last_ear = self._triangle(points_list, total_points_list, concave_points[1])
|
|
ears.append(last_ear)
|
|
return ears
|
|
|
|
def _if_concave_change_status(self, normal, points_list, convex_point, total_points_list,
|
|
concave_points, convex_points, point_in_list):
|
|
"""
|
|
checks whether an convex specific point change its status to concave after removing one ear in the polygon
|
|
returning the new convex and concave points lists together with a flag advising that the list of total points
|
|
already 3 and, therefore, the triangulation must be finished.
|
|
:param normal: normal
|
|
:param points_list: points_list
|
|
:param convex_point: int
|
|
:param total_points_list: [point]
|
|
:param concave_points: [point]
|
|
:param convex_points: [point]
|
|
:param point_in_list: int
|
|
:return: list[points], list[points], boolean
|
|
"""
|
|
end_loop = False
|
|
point = points_list[point_in_list * 3:(point_in_list + 1) * 3]
|
|
pointer = total_points_list.index(point_in_list) - 1
|
|
if pointer < 0:
|
|
pointer = len(total_points_list) - 1
|
|
previous_point = points_list[total_points_list[pointer] * 3:total_points_list[pointer] * 3 + 3]
|
|
pointer = total_points_list.index(point_in_list) + 1
|
|
if pointer >= len(total_points_list):
|
|
pointer = 0
|
|
next_point = points_list[total_points_list[pointer] * 3:total_points_list[pointer] * 3 + 3]
|
|
if self._point_is_concave(normal, point, previous_point, next_point):
|
|
if concave_points[0] > convex_point:
|
|
concave_points.insert(0, convex_point)
|
|
elif concave_points[len(concave_points) - 1] < convex_point:
|
|
concave_points.append(convex_point)
|
|
else:
|
|
for p in range(0, len(concave_points) - 1):
|
|
if concave_points[p] < convex_point < concave_points[p + 1]:
|
|
concave_points.insert(p + 1, convex_point)
|
|
convex_points.remove(convex_point)
|
|
end_loop = True
|
|
return concave_points, convex_points, end_loop
|
|
|
|
def _starting_lists(self, points_list, normal):
|
|
"""
|
|
creates the list of vertices (points) that define the polygon (total_points_list), together with other two lists
|
|
separating points between convex and concave
|
|
:param points_list: points_list
|
|
:param normal: normal
|
|
:return: list[point], list[point], list[point]
|
|
"""
|
|
concave_points = []
|
|
convex_points = []
|
|
# lists of concave and convex points
|
|
# case 1: first point
|
|
point = points_list[0:3]
|
|
previous_point = points_list[len(points_list) - 3:]
|
|
next_point = points_list[3:6]
|
|
index = 0
|
|
total_points_list = [index]
|
|
if self._point_is_concave(normal, point, previous_point, next_point):
|
|
concave_points.append(index)
|
|
else:
|
|
convex_points.append(index)
|
|
# case 2: all points except first and last
|
|
for i in range(0, int((len(points_list)-6)/3)):
|
|
point = points_list[(i+1)*3:(i+2)*3]
|
|
previous_point = points_list[i*3:(i+1)*3]
|
|
next_point = points_list[(i+2)*3:(i+3)*3]
|
|
index = i+1
|
|
total_points_list.append(index)
|
|
if self._point_is_concave(normal, point, previous_point, next_point):
|
|
concave_points.append(index)
|
|
else:
|
|
convex_points.append(index)
|
|
# case 3: last point
|
|
point = points_list[len(points_list) - 3:]
|
|
previous_point = points_list[len(points_list) - 6:len(points_list) - 3]
|
|
next_point = points_list[0:3]
|
|
index = int(len(points_list)/3) - 1
|
|
total_points_list.append(index)
|
|
if self._point_is_concave(normal, point, previous_point, next_point):
|
|
concave_points.append(index)
|
|
else:
|
|
convex_points.append(index)
|
|
return total_points_list, concave_points, convex_points
|
|
|
|
@staticmethod
|
|
def _point_is_concave(normal, point, previous_point, next_point) -> bool:
|
|
"""
|
|
returns whether a point is concave
|
|
:param normal: normal
|
|
:param point: point
|
|
:param previous_point: point
|
|
:param next_point: point
|
|
:return: boolean
|
|
"""
|
|
is_concave = False
|
|
accepted_error = 0.1
|
|
points = np.append(previous_point, point)
|
|
points = np.append(points, next_point)
|
|
points = gh.to_points_matrix(points)
|
|
triangle = Polygon(points)
|
|
error_sum = 0
|
|
for i in range(0, len(normal)):
|
|
error_sum += triangle.normal[i] - normal[i]
|
|
if np.abs(error_sum) < accepted_error:
|
|
is_concave = True
|
|
return is_concave
|
|
|
|
def _triangle(self, points_list, total_points_list, point_position) -> Polygon:
|
|
"""
|
|
creates a triangular polygon out of three points
|
|
:param points_list: points_list
|
|
:param total_points_list: [point]
|
|
:param point_position: int
|
|
:return: polygon
|
|
"""
|
|
index = point_position * 3
|
|
previous_point_index, next_point_index = self._enveloping_points_indices(point_position, total_points_list)
|
|
points = points_list[previous_point_index:previous_point_index + 3]
|
|
points = np.append(points, points_list[index:index + 3])
|
|
points = np.append(points, points_list[next_point_index:next_point_index + 3])
|
|
points = gh.to_points_matrix(points)
|
|
triangle = Polygon(points)
|
|
return triangle
|
|
|
|
@staticmethod
|
|
def _enveloping_points_indices(point_position, total_points_list):
|
|
"""
|
|
due to the fact that the lists are not circular, a method to find the previous and next points
|
|
of an specific one is needed
|
|
:param point_position: int
|
|
:param total_points_list: [point]
|
|
:return: int, int
|
|
"""
|
|
previous_point_index = None
|
|
next_point_index = None
|
|
if point_position == total_points_list[0]:
|
|
previous_point_index = total_points_list[len(total_points_list) - 1] * 3
|
|
next_point_index = total_points_list[1] * 3
|
|
if point_position == total_points_list[len(total_points_list) - 1]:
|
|
previous_point_index = total_points_list[len(total_points_list) - 2] * 3
|
|
next_point_index = total_points_list[0] * 3
|
|
for i in range(1, len(total_points_list)-1):
|
|
if point_position == total_points_list[i]:
|
|
previous_point_index = total_points_list[i - 1] * 3
|
|
next_point_index = total_points_list[i + 1] * 3
|
|
return previous_point_index, next_point_index
|
|
|
|
@staticmethod
|
|
def _enveloping_points(point_to_remove, total_points_list):
|
|
"""
|
|
due to the fact that the lists are not circular, a method to find the previous and next points
|
|
of an specific one is needed
|
|
:param point_to_remove: point
|
|
:param total_points_list: [point]
|
|
:return: point, point
|
|
"""
|
|
index = total_points_list.index(point_to_remove)
|
|
if index == 0:
|
|
previous_point_in_list = total_points_list[len(total_points_list) - 1]
|
|
next_point_in_list = total_points_list[1]
|
|
elif index == len(total_points_list) - 1:
|
|
previous_point_in_list = total_points_list[len(total_points_list) - 2]
|
|
next_point_in_list = total_points_list[0]
|
|
else:
|
|
previous_point_in_list = total_points_list[index - 1]
|
|
next_point_in_list = total_points_list[index + 1]
|
|
return previous_point_in_list, next_point_in_list
|
|
|
|
@staticmethod
|
|
def _is_ear(ear, points) -> bool:
|
|
"""
|
|
finds whether a triangle is an ear of the polygon
|
|
:param ear: polygon
|
|
:param points: [point]
|
|
:return: boolean
|
|
"""
|
|
area_ear = ear.area
|
|
for point in points:
|
|
area_points = 0
|
|
point_is_not_vertex = True
|
|
for i in range(0, 3):
|
|
if abs(np.linalg.norm(point) - np.linalg.norm(ear.points[i])) < 0.0001:
|
|
point_is_not_vertex = False
|
|
break
|
|
if point_is_not_vertex:
|
|
for i in range(0, 3):
|
|
if i != 2:
|
|
new_points = ear.points[i][:]
|
|
new_points = np.append(new_points, ear.points[i + 1][:])
|
|
new_points = np.append(new_points, point[:])
|
|
else:
|
|
new_points = ear.points[i][:]
|
|
new_points = np.append(new_points, point[:])
|
|
new_points = np.append(new_points, ear.points[0][:])
|
|
new_points = gh.to_points_matrix(new_points)
|
|
new_triangle = Polygon(new_points)
|
|
area_points += new_triangle.area
|
|
if abs(area_points - area_ear) < 1e-6:
|
|
# point_inside_ear = True
|
|
return False
|
|
return True
|
|
|
|
@property
|
|
def faces(self) -> [[int]]:
|
|
|
|
"""
|
|
Polyhedron triangular faces
|
|
:return: [[int]]
|
|
"""
|
|
|
|
if self._faces is None:
|
|
self._faces = []
|
|
|
|
for polygon in self._polygons:
|
|
|
|
face = []
|
|
points = polygon.points
|
|
if len(points) != 3:
|
|
sub_polygons = self._triangulate(polygon)
|
|
# todo: I modified this! To be checked @Guille
|
|
if len(sub_polygons) >= 1:
|
|
for sub_polygon in sub_polygons:
|
|
face = []
|
|
points = sub_polygon.points
|
|
for point in points:
|
|
face.append(self._position_of(point, face))
|
|
self._faces.append(face)
|
|
else:
|
|
for point in points:
|
|
face.append(self._position_of(point, face))
|
|
self._faces.append(face)
|
|
return self._faces
|
|
|
|
@property
|
|
def polyhedron_trimesh(self):
|
|
if self._mesh is None:
|
|
self._mesh = Trimesh(vertices=self.vertices, faces=self.faces)
|
|
return self._mesh
|
|
|
|
@property
|
|
def volume(self):
|
|
"""
|
|
Polyhedron volume in cubic meters
|
|
:return: float
|
|
"""
|
|
if self._volume is None:
|
|
if not self.polyhedron_trimesh.is_volume:
|
|
self._volume = np.inf
|
|
else:
|
|
self._volume = self.polyhedron_trimesh.volume
|
|
return self._volume
|
|
|
|
@property
|
|
def max_z(self):
|
|
"""
|
|
Polyhedron maximal z value
|
|
:return: float
|
|
"""
|
|
if self._max_z is None:
|
|
self._max_z = ConfigurationHelper().min_coordinate
|
|
for polygon in self._polygons:
|
|
for point in polygon.points:
|
|
if self._max_z < point[2]:
|
|
self._max_z = point[2]
|
|
return self._max_z
|
|
|
|
@property
|
|
def max_y(self):
|
|
"""
|
|
Polyhedron maximal y value
|
|
:return: float
|
|
"""
|
|
if self._max_y is None:
|
|
self._max_y = ConfigurationHelper().min_coordinate
|
|
for polygon in self._polygons:
|
|
for point in polygon.points:
|
|
if self._max_y < point[1]:
|
|
self._max_y = point[1]
|
|
return self._max_y
|
|
|
|
@property
|
|
def max_x(self):
|
|
"""
|
|
Polyhedron maximal x value
|
|
:return: float
|
|
"""
|
|
if self._max_x is None:
|
|
self._max_x = ConfigurationHelper().min_coordinate
|
|
for polygon in self._polygons:
|
|
for point in polygon.points:
|
|
self._max_x = max(self._max_x, point[0])
|
|
return self._max_x
|
|
|
|
@property
|
|
def min_z(self):
|
|
"""
|
|
Polyhedron minimal z value
|
|
:return: float
|
|
"""
|
|
if self._min_z is None:
|
|
self._min_z = self.max_z
|
|
for polygon in self._polygons:
|
|
for point in polygon.points:
|
|
if self._min_z > point[2]:
|
|
self._min_z = point[2]
|
|
return self._min_z
|
|
|
|
@property
|
|
def min_y(self):
|
|
"""
|
|
Polyhedron minimal y value
|
|
:return: float
|
|
"""
|
|
if self._min_y is None:
|
|
self._min_y = self.max_y
|
|
for polygon in self._polygons:
|
|
for point in polygon.points:
|
|
if self._min_y > point[1]:
|
|
self._min_y = point[1]
|
|
return self._min_y
|
|
|
|
@property
|
|
def min_x(self):
|
|
"""
|
|
Polyhedron minimal x value
|
|
:return: float
|
|
"""
|
|
if self._min_x is None:
|
|
self._min_x = self.max_x
|
|
for polygon in self._polygons:
|
|
for point in polygon.points:
|
|
if self._min_x > point[0]:
|
|
self._min_x = point[0]
|
|
return self._min_x
|
|
|
|
@property
|
|
def center(self):
|
|
"""
|
|
Polyhedron centroid
|
|
:return: [x,y,z]
|
|
"""
|
|
x = (self.max_x + self.min_x) / 2
|
|
y = (self.max_y + self.min_y) / 2
|
|
z = (self.max_z + self.min_z) / 2
|
|
return [x, y, z]
|
|
|
|
def stl_export(self, full_path):
|
|
"""
|
|
Export the polyhedron to stl given file
|
|
:param full_path: str
|
|
:return: None
|
|
"""
|
|
self.polyhedron_trimesh.export(full_path, 'stl_ascii')
|
|
|
|
def obj_export(self, full_path):
|
|
"""
|
|
Export the polyhedron to obj given file
|
|
:param full_path: str
|
|
:return: None
|
|
"""
|
|
self.polyhedron_trimesh.export(full_path, 'obj')
|
|
|
|
def show(self):
|
|
self.polyhedron_trimesh.show()
|