city_retrofit/city_model_structure/attributes/polyhedron.py

"""
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
    print(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()