From 43c32bfceb2ed2e9cc5d61f20de40e7bf13c166b Mon Sep 17 00:00:00 2001
From: pilar <pilar_monsalvete@yahoo.es>
Date: Tue, 8 Dec 2020 16:48:07 -0500
Subject: [PATCH] actualizing dynamic workflow

---
 city_model_structure/attributes/surface.py | 83 +++++++++++++---------
 1 file changed, 51 insertions(+), 32 deletions(-)

diff --git a/city_model_structure/attributes/surface.py b/city_model_structure/attributes/surface.py
index 461bff60..f41dfd61 100644
--- a/city_model_structure/attributes/surface.py
+++ b/city_model_structure/attributes/surface.py
@@ -219,40 +219,17 @@ class Surface:
     """
     # New method to calculate area
     if self._area is None:
-      if len(self.points) < 3:
+      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 += self.angle_between_vectors(vec_1, vec_2)
+      if len(self.points) < 3 or alpha == 0:
         area = 0
+        print('Warning: the area of a line or point cannot be calculated. Area = 0')
       else:
-        z_vector = [0, 0, 1]
-        normal_vector = self.normal
-        points_2d = []
-        x = normal_vector[0]
-        y = normal_vector[1]
-
-        if x != 0 or y != 0:
-          alpha = math.acos(np.dot(normal_vector, z_vector) / np.linalg.norm(normal_vector) / np.linalg.norm(z_vector))
-          turning_line = np.cross(normal_vector, z_vector)
-          third_axis = np.cross(normal_vector, turning_line)
-          w_1 = turning_line / np.linalg.norm(turning_line)
-          w_2 = normal_vector
-          w_3 = third_axis / np.linalg.norm(third_axis)
-          # turning_base_matrix
-          turning_matrix = np.array([[1, 0, 0],
-                                    [0, math.cos(alpha), -math.sin(alpha)],
-                                    [0, math.sin(alpha), math.cos(alpha)]])
-          base_matrix = np.array([w_1, w_2, w_3])
-          turning_base_matrix = np.matmul(base_matrix.transpose(), turning_matrix.transpose())
-          turning_base_matrix = np.matmul(turning_base_matrix, base_matrix)
-
-          if turning_base_matrix is None:
-            print('Error processing turning base matrix')
-          else:
-            for point in self.points:
-              new_point = np.matmul(turning_base_matrix, point)
-              points_2d.append(new_point)
-        else:
-          for point in self.points:
-            points_2d.append([point[0], point[1], 0])
-
+        # todo: Guilleeeee!!!!
+        points_2d = self.rotate_surface_to_horizontal(self)
         polygon_2d = pn.Polygon(np.array(points_2d))
         area = 0
         for i in range(0, len(polygon_2d.points)-1):
@@ -265,6 +242,48 @@ class Surface:
       self._area = abs(area)
     return self._area
 
+  @staticmethod
+  def rotate_surface_to_horizontal(surface):
+    z_vector = [0, 0, 1]
+    normal_vector = surface.normal
+    points_2d = []
+    x = normal_vector[0]
+    y = normal_vector[1]
+
+    if x == 0 and y == 0:
+      # Already horizontal
+      for point in surface.points:
+        points_2d.append([point[0], point[1], 0])
+    else:
+      alpha = surface.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, math.cos(alpha), -math.sin(alpha)],
+                                 [0, math.sin(alpha), math.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:
+        print('Error processing rotation base matrix')
+      else:
+        for point in surface.points:
+          new_point = np.matmul(rotation_base_matrix, point)
+          points_2d.append(new_point)
+    return points_2d
+
+  @staticmethod
+  def angle_between_vectors(vec_1, vec_2):
+    alpha = math.acos(np.dot(vec_1, vec_2) / np.linalg.norm(vec_1) / np.linalg.norm(vec_2))
+    if np.linalg.norm(vec_1) == 0 or np.linalg.norm(vec_2) == 0:
+      print('ERROR')
+      # todo: rise error
+    return alpha
+
   def _is_almost_same_terrain(self, terrain_points, ground_points):
     equal = 0
     for terrain_point in terrain_points: