diff --git a/city_model_structure/attributes/polyhedron.py b/city_model_structure/attributes/polyhedron.py index 3b6145a1..0897545c 100644 --- a/city_model_structure/attributes/polyhedron.py +++ b/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 - points = ' '.join(str(e) for e in [*points_list[len(points_list) - 3:], *points_list[0:6]]) - triangle = Surface(points, remove_last=False) - if self._point_is_concave(normal, triangle): - concave_points.append(points_list[0:3]) + # 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(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]]) - triangle = Surface(points, remove_last=False) - if self._point_is_concave(normal, triangle): - concave_points.append(points_list[(i+1)*3:(i+2)*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(points_list[(i+1)*3:(i+2)*3]) - # case 3: last point as center ear - points = ' '.join(str(e) for e in [*points_list[len(points_list) - 6:], *points_list[0:3]]) - triangle = Surface(points, remove_last=False) - if self._point_is_concave(normal, triangle): - concave_points.append(points_list[len(points_list) - 3:]) + 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(points_list[len(points_list) - 3:]) -# print('point list', points_list) - print('concave', concave_points) - print('convex', convex_points) + convex_points.append(index) + # list of ears + ears = [] + # 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 -# for point in concave_points: -# 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: + @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 - 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]]: diff --git a/city_model_structure/attributes/surface.py b/city_model_structure/attributes/surface.py index 29123816..396d7928 100644 --- a/city_model_structure/attributes/surface.py +++ b/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 - counter_clockwise = 0 + alpha += beta else: - clockwise = 0 - counter_clockwise = 1 - for i in range(0, len(points)-2): + alpha -= beta + for i in range(0, len(points)-4): 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 - if clockwise < counter_clockwise: + alpha -= beta + if alpha < 0: cross_product = np.cross(points[0] - points[len(points) - 2], points[len(points) - 1] - points[len(points) - 2]) else: diff --git a/tests/test_geometry_factory.py b/tests/test_geometry_factory.py index 764ca5ac..939b4f63 100644 --- a/tests/test_geometry_factory.py +++ b/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 -import numpy as np +from datetime import datetime 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' - file_path = (self._example_path / building).resolve() - city = GeometryFactory('citygml', file_path).city + now = datetime.now() + print("start =", now) +# 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) - #print('volume om', building.volume_om) - building.stl_export(self._output_path) + if building.name != 'BLD121958': +# if building.name == 'BLD101288': + 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) -