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
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:
  """
  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._max_y = None
    self._max_x = None
    self._min_z = None
    self._min_y = None
    self._min_x = None
    self._geometry = GeometryHelper(delta=5.0, area_delta=0.01)

  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 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, True):
            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 _get_regions(point_index, points_list):
    if point_index == 0:
      # first point in the polygon so the triangle is the points n-1, n, 0
      triangle_left = ' '.join(str(e) for e in [*points_list[len(points_list) - 6:], *points_list[0:3]])
      # remove point n
      rest_points_left = ' '.join(str(e) for e in [*points_list[:len(points_list) - 3]])
    elif point_index == 3:
      # second point in the polygon so the triangle is the points n, 0, 1
      triangle_left = ' '.join(str(e) for e in [*points_list[len(points_list) - 3:], *points_list[0:6]])
      # remove point 0
      rest_points_left = ' '.join(str(e) for e in [*points_list[3:]])
    else:
      # normal point index-2¸index-1, index
      triangle_left = ' '.join(str(e) for e in [*points_list[point_index - 6:point_index + 3]])
      # remove middle point (index - 1)
      rest_points_left = ' '.join(str(e) for e in [*points_list[0:point_index - 3], *points_list[point_index:]])
    if point_index < len(points_list) - 6:
      # normal point index, index+1, index+2
      triangle_right = ' '.join(str(e) for e in [*points_list[point_index:point_index + 9]])
      rest_points_right = ' '.join(str(e) for e in [*points_list[0:point_index + 3], *points_list[point_index + 6:]])
    elif point_index == (len(points_list) - 6):
      # last two points in the polygon so the triangle is the points n-1, n, 0
      triangle_right = ' '.join(str(e) for e in [*points_list[point_index:], *points_list[0:3]])
      rest_points_right = ' '.join(str(e) for e in [*points_list[:len(points_list - 3)]])
    else:
      # last point in the polygon so the triangle is n, 0, 1
      triangle_right = ' '.join(str(e) for e in [*points_list[len(points_list) - 3:], *points_list[0:6]])
      rest_points_right = ' '.join(str(e) for e in [*points_list[3:]])
    return (Surface(triangle_left, remove_last=False), Surface(rest_points_left, remove_last=False)), \
           (Surface(triangle_right, remove_last=False), Surface(rest_points_right, remove_last=False))

  def _triangulate(self, surface):
    triangles = []
    complementary_surface = surface
    triangles_count = len(surface.points) - 2
    points_list = surface.points_list
    point_index = 0
    area = surface.area
    while len(triangles) < triangles_count:
      # get triangles and regions in both direction to find ears
      left_direction, right_direction = Polyhedron._get_regions(point_index, points_list)
      # todo: use enum to describe triangle or rest instead 0 1
      right_area = right_direction[0].area + right_direction[1].area
      left_area = left_direction[0].area + left_direction[1].area
      if self._geometry.almost_same_area(area, left_area) and self._geometry.almost_same_area(area, right_area):
        # Both seems to be an ear, choose the more precise
        if np.abs(left_area-area) < np.abs(right_area-area):
          area = left_direction[1].area
          point_index = 0
          triangles.append(left_direction[0])
          points_list = left_direction[1].points_list
          complementary_surface = left_direction[1]
        else:
          area = right_direction[1].area
          point_index = 0
          triangles.append(right_direction[0])
          points_list = right_direction[1].points_list
          complementary_surface = right_direction[1]
      elif self._geometry.almost_same_area(area, left_area):
        area = left_direction[1].area
        point_index = 0
        triangles.append(left_direction[0])
        points_list = left_direction[1].points_list
        complementary_surface = left_direction[1]
      elif self._geometry.almost_same_area(area, right_area):
        area = right_direction[1].area
        point_index = 0
        triangles.append(right_direction[0])
        points_list = right_direction[1].points_list
        complementary_surface = right_direction[1]
      else:
        point_index = point_index + 3
        if point_index >= len(points_list):
          # print('wroooong', complementary_surface.type)
          # print(complementary_surface.points)
          return triangles
      if len(points_list) == 9:
        # the rest point's are already a triangle
        triangles.append(complementary_surface)

    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:
          sub_surfaces = self._triangulate(surface)
          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_mesh(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_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
    """
    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_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()