system_assignation/city_model_structure/attributes/polygon.py

"""
Polygon module
SPDX - License - Identifier: LGPL - 3.0 - or -later
Copyright © 2020 Project Author Pilar Monsalvete Álvarez de Uribarri pilar.monsalvete@concordia.ca
"""

import sys
import numpy as np

from helpers.geometry_helper import GeometryHelper as gh


class Polygon:
  """
  Polygon class
  """

  def __init__(self, vertices):
    self._vertices = vertices
    self._remove_last = True
    self._area = None
    self._points = None
    self._normal = None

  @property
  def points(self) -> np.ndarray:
    """
    Surface point matrix
    :return: np.ndarray
    """
    if self._points is None:
      self._points = np.fromstring(self._vertices, dtype=float, sep=' ')
      self._points = gh.to_points_matrix(self._points, self._remove_last)
    return self._points

  @property
  def area(self):
    """
    Surface area in square meters
    :return: float
    """
    # New method to calculate area
    if self._area is None:
      if len(self.points) < 3:
        sys.stderr.write('Warning: the area of a line or point cannot be calculated 1. Area = 0\n')
        return 0
      alpha = 0
      vec_1 = self.points[1] - self.points[0]
      for i in range(2, len(self.points)):
        vec_2 = self.points[i] - self.points[0]
        alpha += gh.angle_between_vectors(vec_1, vec_2)
      if alpha == 0:
        sys.stderr.write('Warning: the area of a line or point cannot be calculated 2. Area = 0\n')
        return 0
      horizontal_points = self.rotate_surface_to_horizontal
      area = 0
      for i in range(0, len(horizontal_points)-1):
        point = horizontal_points[i]
        next_point = horizontal_points[i+1]
        area += (next_point[1] + point[1]) / 2 * (next_point[0] - point[0])
      next_point = horizontal_points[0]
      point = horizontal_points[len(horizontal_points)-1]
      area += (next_point[1] + point[1]) / 2 * (next_point[0] - point[0])
      self._area = abs(area)
    return self._area

  @property
  def rotate_surface_to_horizontal(self):
    z_vector = [0, 0, 1]
    normal_vector = self.normal
    horizontal_points = []
    x = normal_vector[0]
    y = normal_vector[1]

    if x == 0 and y == 0:
      # Already horizontal
      for point in self.points:
        horizontal_points.append([point[0], point[1], 0])
    else:
      alpha = gh.angle_between_vectors(normal_vector, z_vector)
      rotation_line = np.cross(normal_vector, z_vector)
      third_axis = np.cross(normal_vector, rotation_line)
      w_1 = rotation_line / np.linalg.norm(rotation_line)
      w_2 = normal_vector
      w_3 = third_axis / np.linalg.norm(third_axis)
      rotation_matrix = np.array([[1, 0, 0],
                                 [0, np.cos(alpha), -np.sin(alpha)],
                                 [0, np.sin(alpha), np.cos(alpha)]])
      base_matrix = np.array([w_1, w_2, w_3])
      rotation_base_matrix = np.matmul(base_matrix.transpose(), rotation_matrix.transpose())
      rotation_base_matrix = np.matmul(rotation_base_matrix, base_matrix)

      if rotation_base_matrix is None:
        sys.stderr.write('Warning: rotation base matrix returned None\n')
      else:
        for point in self.points:
          new_point = np.matmul(rotation_base_matrix, point)
          horizontal_points.append(new_point)
    return horizontal_points

  @property
  def normal(self) -> np.ndarray:
    """
    Surface normal vector
    :return: np.ndarray
    """
    if self._normal is None:
      points = self.points
      # todo: IF THE FIRST ONE IS 0, START WITH THE NEXT
      point_origin = points[len(points)-2]
      vector_1 = points[len(points)-1] - point_origin
      vector_2 = points[0] - point_origin
      vector_3 = points[1] - point_origin
      cross_product = np.cross(vector_1, vector_2)
      if np.linalg.norm(cross_product) != 0:
        cross_product = cross_product / np.linalg.norm(cross_product)
        alpha = gh.angle_between_vectors(vector_1, vector_2)
      else:
        # todo modify here
        cross_product = [0, 0, 0]
        alpha = 0
      if len(points) == 3:
        return cross_product
      alpha += self._angle(vector_2, vector_3, cross_product)
      for i in range(0, len(points)-4):
        vector_1 = points[i+1] - point_origin
        vector_2 = points[i+2] - point_origin
        alpha += self._angle(vector_1, vector_2, cross_product)
      vector_1 = points[len(points) - 1] - point_origin
      vector_2 = points[0] - point_origin
      if alpha < 0:
        cross_product = np.cross(vector_2, vector_1)
      else:
        cross_product = np.cross(vector_1, vector_2)
      self._normal = cross_product / np.linalg.norm(cross_product)
    return self._normal

  @staticmethod
  def _angle(vector_1, vector_2, cross_product):
    accepted_normal_difference = 0.01
    cross_product_next = np.cross(vector_1, vector_2)
    if np.linalg.norm(cross_product_next) != 0:
      cross_product_next = cross_product_next / np.linalg.norm(cross_product_next)
      alpha = gh.angle_between_vectors(vector_1, vector_2)
    else:
      cross_product_next = [0, 0, 0]
      alpha = 0
    delta_normals = 0
    for j in range(0, 3):
      delta_normals += cross_product[j] - cross_product_next[j]
    if np.abs(delta_normals) < accepted_normal_difference:
      return alpha
    else:
      return -alpha