""" 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.surface import Surface class Polyhedron: # todo: modify class to call polygons better than surfaces """ Polyhedron class """ def __init__(self, surfaces): self._surfaces = list(surfaces) self._polygons = [s.polygon for s in surfaces] 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._surfaces] 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 @staticmethod def _point(coordinates): return coordinates[0], coordinates[1], coordinates[2] def _triangulate(self, surface): """ triangulates a polygon following the ear clipping methodology :param surface: surface :return: list[triangles] """ points_list = surface.points_list normal = surface.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(surface.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 surface type {surface.type}\n') return [surface] 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 = ' '.join(str(e) for e in [*previous_point[:], *point[:], *next_point[:]]) triangle = Surface(points, remove_last=False) error_sum = 0 for i in range(0, len(normal)): error_sum += triangle.perimeter_surface.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): """ 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 = ' '.join(str(e) for e in [*points_list[previous_point_index:previous_point_index + 3], *points_list[index:index + 3], *points_list[next_point_index:next_point_index + 3]]) triangle = Surface(points, remove_last=False) 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: surface :param points: [point] :return: boolean """ area_ear = ear.perimeter_surface.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 = ' '.join(str(e) for e in [*ear.points[i][:], *ear.points[i + 1][:], *point[:]]) else: new_points = ' '.join(str(e) for e in [*ear.points[i][:], *point[:], *ear.points[0][:]]) new_triangle = Surface(new_points, remove_last=False) area_points += new_triangle.perimeter_surface.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 surface in self._surfaces: face = [] points = surface.points if len(points) != 3: sub_surfaces = self._triangulate(surface) # todo: I modified this! To be checked @Guille if len(sub_surfaces) >= 1: for sub_surface in sub_surfaces: face = [] points = sub_surface.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 surface in self._surfaces: for point in surface.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 surface in self._surfaces: for point in surface.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 surface in self._surfaces: for point in surface.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 surface in self._surfaces: for point in surface.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 surface in self._surfaces: for point in surface.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 surface in self._surfaces: for point in surface.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()