new triangulate function in polyhedron.py working. v1.0

city_model_structure/attributes/polyhedron.py

@@ -66,61 +66,170 @@ class Polyhedron:
     return coordinates[0], coordinates[1], coordinates[2]
 
   def _triangulate(self, surface):
-    print(surface.type)
-    triangles = []
-    triangles_count = len(surface.points) - 2
     points_list = surface.points_list
     normal = surface.normal
     concave_points = []
     convex_points = []
-    # case 1: first point as center ear
+    # lists of concave and convex points
-    points = ' '.join(str(e) for e in [*points_list[len(points_list) - 3:], *points_list[0:6]])
+    # case 1: first point
-    triangle = Surface(points, remove_last=False)
+    point = points_list[0:3]
-    if self._point_is_concave(normal, triangle):
+    previous_point = points_list[len(points_list) - 3:]
-      concave_points.append(points_list[0: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(points_list[0:3])
+      convex_points.append(index)
     # case 2: all points except first and last
     for i in range(0, int((len(points_list)-6)/3)):
-      points = ' '.join(str(e) for e in [*points_list[i*3:(i+3)*3]])
+      point = points_list[(i+1)*3:(i+2)*3]
-      triangle = Surface(points, remove_last=False)
+      previous_point = points_list[i*3:(i+1)*3]
-      if self._point_is_concave(normal, triangle):
+      next_point = points_list[(i+2)*3:(i+3)*3]
-        concave_points.append(points_list[(i+1)*3:(i+2)*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(points_list[(i+1)*3:(i+2)*3])
+        convex_points.append(index)
-    # case 3: last point as center ear
+    # case 3: last point
-    points = ' '.join(str(e) for e in [*points_list[len(points_list) - 6:], *points_list[0:3]])
+    point = points_list[len(points_list) - 3:]
-    triangle = Surface(points, remove_last=False)
+    previous_point = points_list[len(points_list) - 6:len(points_list) - 3]
-    if self._point_is_concave(normal, triangle):
+    next_point = points_list[0:3]
-      concave_points.append(points_list[len(points_list) - 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(points_list[len(points_list) - 3:])
+      convex_points.append(index)
-#    print('point list', points_list)
+    # list of ears
-    print('concave', concave_points)
+    ears = []
-    print('convex', convex_points)
+    # todo remove counter
+    counter = 0
+    while len(concave_points) > 3 or len(convex_points) != 0 and counter < 40:
+      counter += 1
+      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]
+          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]
+          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
+          # todo CHECK PREVIOUS AND NEXT POINTS IF OUT OF RANGE
+          for j in range(0, len(convex_points)):
+            if convex_points[j] == previous_point_in_list:
+              point = points_list[previous_point_in_list*3:(previous_point_in_list + 1)*3]
+              pointer = total_points_list.index(previous_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(previous_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_points[j]:
+                  concave_points.insert(0, convex_points[j])
+                elif concave_points[len(concave_points)-1] < convex_points[j]:
+                  concave_points.append(convex_points[j])
+                else:
+                  for p in range(0, len(concave_points)-1):
+                    if concave_points[p] < convex_points[j] < concave_points[p+1]:
+                      concave_points.insert(p, convex_points[j])
+                convex_points.remove(convex_points[j])
+                break
+              continue
+            if convex_points[j] == next_point_in_list:
+              point = points_list[next_point_in_list*3:(next_point_in_list + 1)*3]
+              pointer = total_points_list.index(next_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(next_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):
+                concave_points.insert(0, convex_points[j])
+                convex_points.remove(convex_points[j])
+                break
+              continue
+          break
+    last_ear = self._triangle(points_list, total_points_list, concave_points[1])
+    ears.append(last_ear)
+    return ears
 
-    # todo: recursive function, not good solution
+  @staticmethod
-#    for point in concave_points:
+  def _point_is_concave(normal, point, previous_point, next_point) -> bool:
-#      is_ear_point = self._is_ear_point(point)
-#      if is_ear_point:
-#        ear = self._extract_ear()
-#        triangles.append(ear)
-#        self._remove_point()
-#        continue
-    return triangles
-
-  def _point_is_concave(self, normal, triangle) -> 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
-    print('normal', normal)
-    print('normal triangle', triangle.normal)
     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
 
+  @staticmethod
+  def _triangle(points_list, total_points_list, point_position):
+    index = point_position * 3
+    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
+    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 _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]]:
 

city_model_structure/attributes/surface.py

@@ -279,29 +279,51 @@ class Surface:
     if self._normal is None:
       points = self.points
       accepted_error = 0.01
+      # todo: IF THE FIRST ONE IS 0, START WITH THE NEXT
       cross_product = np.cross(points[len(points)-1] - points[len(points)-2],
                                points[0] - points[len(points)-2])
       cross_product_next = np.cross(points[0] - points[len(points)-2],
                                     points[1] - points[len(points)-2])
+      if np.linalg.norm(cross_product) != 0:
+        cross_product = cross_product / np.linalg.norm(cross_product)
+        alpha = GeometryHelper.angle_between_vectors(points[len(points)-1] - points[len(points)-2],
+                                                     points[0] - points[len(points)-2])
+      else:
+        cross_product = [0, 0, 0]
+        alpha = 0
+      if len(points) == 3:
+        return cross_product
+      if np.linalg.norm(cross_product_next) != 0:
+        cross_product_next = cross_product_next / np.linalg.norm(cross_product_next)
+        beta = GeometryHelper.angle_between_vectors(points[0] - points[len(points)-2],
+                                                    points[1] - points[len(points)-2])
+      else:
+        cross_product_next = [0, 0, 0]
+        beta = 0
       error_sum = 0
       for j in range(0, 3):
         error_sum += cross_product[j] - cross_product_next[j]
       if np.abs(error_sum) < accepted_error:
-        clockwise = 1
+        alpha += beta
-        counter_clockwise = 0
       else:
-        clockwise = 0
+        alpha -= beta
-        counter_clockwise = 1
+      for i in range(0, len(points)-4):
-      for i in range(0, len(points)-2):
         cross_product_next = np.cross(points[i+1] - points[len(points)-2], points[i+2] - points[len(points)-2])
+        if np.linalg.norm(cross_product_next) != 0:
+          cross_product_next = cross_product_next / np.linalg.norm(cross_product_next)
+          beta = GeometryHelper.angle_between_vectors(points[i+1] - points[len(points)-2],
+                                                      points[i+2] - points[len(points)-2])
+        else:
+          cross_product_next = [0, 0, 0]
+          beta = 0
         error_sum = 0
         for j in range(0, 3):
           error_sum += cross_product[j] - cross_product_next[j]
         if np.abs(error_sum) < accepted_error:
-          clockwise += 1
+          alpha += beta
         else:
-          counter_clockwise += 1
+          alpha -= beta
-      if clockwise < counter_clockwise:
+      if alpha < 0:
         cross_product = np.cross(points[0] - points[len(points) - 2],
                                  points[len(points) - 1] - points[len(points) - 2])
       else:

tests/test_geometry_factory.py

@@ -10,8 +10,7 @@ from unittest import TestCase
 from city_model_structure.city import City
 from factories.geometry_factory import GeometryFactory
 from helpers.geometry_helper import GeometryHelper
-import math
+from datetime import datetime
-import numpy as np
 
 
 class TestGeometryFactory(TestCase):
@@ -25,7 +24,7 @@ class TestGeometryFactory(TestCase):
     """
     self._city_gml = None
     self._example_path = (Path(__file__).parent.parent / 'tests_data').resolve()
-    self._pickle_file = (self._example_path / 'city.pickle').resolve()
+    self._pickle_file = (self._example_path / 'city_new.pickle').resolve()
     self._kelowna_pickle_file = (self._example_path / 'kelowna_test_case.pickle').resolve()
     self._output_path = (Path(__file__).parent / 'surface_outputs').resolve()
 
@@ -79,9 +78,10 @@ class TestGeometryFactory(TestCase):
     :return:
     """
     city = City.load(self._kelowna_pickle_file)
-    helper = GeometryHelper(delta=0.0, area_delta=0.5)
+    helper = GeometryHelper(delta=0.0, area_delta=0.0)
     errors = 0
     for building in city.buildings:
+      print(building.name)
       building._polyhedron._geometry = helper
       if str(building.volume) == 'inf':
         building.stl_export(self._output_path)
@@ -93,13 +93,30 @@ class TestGeometryFactory(TestCase):
 
     :return:
     """
-    building = 'bld100087.gml'
+    now = datetime.now()
-    file_path = (self._example_path / building).resolve()
+    print("start =", now)
-    city = GeometryFactory('citygml', file_path).city
+#    file_path = (self._example_path / 'gml_17_12_2020.gml').resolve()
+#    city = GeometryFactory('citygml', file_path).city
+#    city.save(self._pickle_file)
+#    now = datetime.now()
+#    print("end pickle =", now)
+
+    city = City.load(self._pickle_file)
+    helper = GeometryHelper(delta=0.0, area_delta=0.0)
+
+    counter = 0
     for building in city.buildings:
-      print('volume', building.volume)
+      if building.name != 'BLD121958':
-      #print('volume om', building.volume_om)
+#      if building.name == 'BLD101288':
-      building.stl_export(self._output_path)
+        building._polyhedron._geometry = helper
+        print('building name', building.name)
+        print('volume', building.name, building.volume)
+        if str(building.volume) == 'inf':
+          counter += 1
+        building.stl_export(self._output_path)
+    print('total number of buildings with volume inf', counter)
+    now = datetime.now()
+    print("end =", now)
 
   def test_citygml_buildings(self):
     """
@@ -238,10 +255,3 @@ class TestGeometryFactory(TestCase):
             self.assertIsNone(thermal_opening.thickness, 'thermal_opening thickness_m was initialized')
             self.assertRaises(Exception, lambda: thermal_opening.u_value, 'thermal_opening u_value was initialized')
 
-  def tests_with_building_bld100087(self):
-    building = 'bld100087.gml'
-    file_path = (self._example_path / building).resolve()
-    city = GeometryFactory('citygml', file_path).city
-    for building in city.buildings:
-      print('volume', building.volume)
-