""" Polyhedron module SPDX - License - Identifier: LGPL - 3.0 - or -later Copyright © 2020 Project Author Guille Gutierrez guillermo.gutierrezmorote@concordia.ca """ import numpy as np import matplotlib.pyplot as plt from trimesh import Trimesh from helpers.geometry_helper import GeometryHelper from city_model_structure.attributes.surface import Surface class Polyhedron: """ Polyhedron class """ def __init__(self, surfaces): self._surfaces = list(surfaces) self._polygons = [s.polygon for s in surfaces] self._polyhedron = None self._volume = None self._faces = None self._vertices = None self._mesh = None self._centroid = None self._max_z = None self._geometry = GeometryHelper() def _position_of(self, point): vertices = self.vertices for i in range(len(vertices)): if self._geometry.almost_equal(vertices[i], point): 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 self._geometry.almost_equal(vertex_1, vertex_2): 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] @staticmethod def _ground_weird_shape(points, original, triangle=False): x = [] y = [] xo = [] yo = [] for o in original: xo.append(o[0]) yo.append(o[1]) for point in points: x.append(point[0]) y.append(point[1]) x = [val - min(xo) for val in x] y = [val - min(yo) for val in y] if triangle: print(len(points)) for point in range(0, len(x)): print('point', x[point], y[point]) def _triangulate(self, surface): Polyhedron._ground_weird_shape(surface.points, surface.points) triangles = [] triangles_count = len(surface.points) - 2 points_list = surface.points_list point_index = 0 area = surface.area print(f'area is {area}') while len(triangles) < triangles_count: print(f'try vertex {point_index}') # select a triangle starting at point index triangle_points = ' '.join(str(e) for e in [*points_list[point_index:point_index + 9]]) # remove the middle vertex from previous triangle rest_points = ' '.join(str(e) for e in [*points_list[0:point_index+3], *points_list[point_index+6:]]) triangular_surface = Surface(triangle_points, remove_last=False) rest_surface = Surface(rest_points, remove_last=False) if self._geometry.almost_same_area(area, (triangular_surface.area + rest_surface.area)): area = rest_surface.area print(f'ok! new area is {area}') print('Triangle-----------------------------------------') Polyhedron._ground_weird_shape(triangular_surface.points, surface.points) print('rest---------------------------------------------') Polyhedron._ground_weird_shape(rest_surface.points, surface.points) triangles.append(triangular_surface) points_list = rest_surface.points_list if len(rest_surface.points) == 3: triangles.append(rest_surface) point_index = 0 else: print(f'nok! area {(triangular_surface.area + rest_surface.area)} is not {area}') print('Triangle-----------------------------------------') Polyhedron._ground_weird_shape(triangular_surface.points, surface.points) print('rest---------------------------------------------') Polyhedron._ground_weird_shape(rest_surface.points, surface.points) point_index = point_index + 3 return triangles @property def faces(self) -> [[int]]: """ Polyhedron faces :return: [[int]] """ if self._faces is None: self._faces = [] for surface in self._surfaces: face = [] points = surface.points if len(points) != 3: # three vertices print(f'Number of vertex {len(points)}') sub_surfaces = self._triangulate(surface) print(f'Transformed {len(points)} vertex into {len(sub_surfaces)} triangles') print(f'Total area: {surface.area}') sum_area = 0.0 for sub_surface in sub_surfaces: sum_area = sum_area + sub_surface.area points = sub_surface.points for point in points: face.append(self._position_of(point)) self._faces.append(face) print(f'Accumulated area: {sum_area}') else: for point in points: face.append(self._position_of(point)) self._faces.append(face) return self._faces @property def _cloud_mesh(self): if self._mesh is None: self._mesh = GeometryHelper.create_mesh(self._surfaces) return self._mesh @property def _polyhedron_mesh(self): if self._mesh is None: try: print("create mesh") self._mesh = Trimesh(vertices=self.vertices, faces=self.faces) except SyntaxError: self._mesh = self._cloud_mesh return self._mesh @property def volume(self): """ Polyhedron volume in cubic meters :return: float """ if self._volume is None: if not self._polyhedron_mesh.is_volume: self._volume = np.inf else: self._volume = self._polyhedron_mesh.volume return self._volume @property def max_z(self): """ Polyhedron maximal z value :return: float """ bounds = self._polyhedron_mesh.bounds z_max = max(bounds[:, 2]) return z_max @property def centroid(self): """ Polyhedron centroid :return: [x,y,z] """ return self._polyhedron_mesh.centroid def stl_export(self, full_path): """ Export the polyhedron to stl given file :param full_path: str :return: None """ self._polyhedron_mesh.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_mesh.export(full_path, 'obj') def show(self): self._polyhedron_mesh.show()