hub/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
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:
  """
  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()

  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):
    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):
    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):
    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:
    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.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):
    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):
    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):
    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:
    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 = ' '.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.area
      if abs(area_points - area_ear) < 1e-6:
        # point_inside_ear = True
        return False
    return True

  @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)
          # 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_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()