From 1fccac8f4c407cfed93ebc5c98ebf7c861d96a40 Mon Sep 17 00:00:00 2001 From: Pilar Date: Wed, 13 Jan 2021 16:41:45 -0500 Subject: [PATCH] first attempt to new triangulate function in polyhedron.py --- city_model_structure/attributes/polyhedron.py | 121 +++--- city_model_structure/attributes/surface.py | 115 +++--- helpers/geometry_helper.py | 12 +- tests/test_geometry_factory.py | 4 +- tests_data/bld100087.gml | 373 ++++++++++++++++++ 5 files changed, 501 insertions(+), 124 deletions(-) create mode 100644 tests_data/bld100087.gml diff --git a/city_model_structure/attributes/polyhedron.py b/city_model_structure/attributes/polyhedron.py index 1c501f39..3b6145a1 100644 --- a/city_model_structure/attributes/polyhedron.py +++ b/city_model_structure/attributes/polyhedron.py @@ -65,87 +65,62 @@ class Polyhedron: def _point(coordinates): return coordinates[0], coordinates[1], coordinates[2] - @staticmethod - def _get_regions(point_index, points_list): - if point_index == 0: - # first point in the polygon so the triangle is the points n-1, n, 0 - triangle_left = ' '.join(str(e) for e in [*points_list[len(points_list) - 6:], *points_list[0:3]]) - # remove point n - rest_points_left = ' '.join(str(e) for e in [*points_list[:len(points_list) - 3]]) - elif point_index == 3: - # second point in the polygon so the triangle is the points n, 0, 1 - triangle_left = ' '.join(str(e) for e in [*points_list[len(points_list) - 3:], *points_list[0:6]]) - # remove point 0 - rest_points_left = ' '.join(str(e) for e in [*points_list[3:]]) - else: - # normal point index-2¸index-1, index - triangle_left = ' '.join(str(e) for e in [*points_list[point_index - 6:point_index + 3]]) - # remove middle point (index - 1) - rest_points_left = ' '.join(str(e) for e in [*points_list[0:point_index - 3], *points_list[point_index:]]) - if point_index < len(points_list) - 6: - # normal point index, index+1, index+2 - triangle_right = ' '.join(str(e) for e in [*points_list[point_index:point_index + 9]]) - rest_points_right = ' '.join(str(e) for e in [*points_list[0:point_index + 3], *points_list[point_index + 6:]]) - elif point_index == (len(points_list) - 6): - # last two points in the polygon so the triangle is the points n-1, n, 0 - triangle_right = ' '.join(str(e) for e in [*points_list[point_index:], *points_list[0:3]]) - rest_points_right = ' '.join(str(e) for e in [*points_list[:len(points_list - 3)]]) - else: - # last point in the polygon so the triangle is n, 0, 1 - triangle_right = ' '.join(str(e) for e in [*points_list[len(points_list) - 3:], *points_list[0:6]]) - rest_points_right = ' '.join(str(e) for e in [*points_list[3:]]) - return (Surface(triangle_left, remove_last=False), Surface(rest_points_left, remove_last=False)), \ - (Surface(triangle_right, remove_last=False), Surface(rest_points_right, remove_last=False)) - def _triangulate(self, surface): + print(surface.type) triangles = [] - complementary_surface = surface triangles_count = len(surface.points) - 2 points_list = surface.points_list - point_index = 0 - area = surface.area - while len(triangles) < triangles_count: - # get triangles and regions in both direction to find ears - left_direction, right_direction = Polyhedron._get_regions(point_index, points_list) - # todo: use enum to describe triangle or rest instead 0 1 - right_area = right_direction[0].area + right_direction[1].area - left_area = left_direction[0].area + left_direction[1].area - if self._geometry.almost_same_area(area, left_area) and self._geometry.almost_same_area(area, right_area): - # Both seems to be an ear, choose the more precise - if np.abs(left_area-area) < np.abs(right_area-area): - area = left_direction[1].area - point_index = 0 - triangles.append(left_direction[0]) - points_list = left_direction[1].points_list - complementary_surface = left_direction[1] - else: - area = right_direction[1].area - point_index = 0 - triangles.append(right_direction[0]) - points_list = right_direction[1].points_list - complementary_surface = right_direction[1] - elif self._geometry.almost_same_area(area, left_area): - area = left_direction[1].area - point_index = 0 - triangles.append(left_direction[0]) - points_list = left_direction[1].points_list - complementary_surface = left_direction[1] - elif self._geometry.almost_same_area(area, right_area): - area = right_direction[1].area - point_index = 0 - triangles.append(right_direction[0]) - points_list = right_direction[1].points_list - complementary_surface = right_direction[1] + 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]) + else: + convex_points.append(points_list[0:3]) + # 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]) else: - point_index = point_index + 3 - if point_index >= len(points_list): - return triangles - if len(points_list) == 9: - # the rest point's are already a triangle - triangles.append(complementary_surface) + 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:]) + else: + convex_points.append(points_list[len(points_list) - 3:]) +# print('point list', points_list) + print('concave', concave_points) + print('convex', convex_points) + # 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: + is_concave = False + accepted_error = 0.1 + 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 + @property def faces(self) -> [[int]]: diff --git a/city_model_structure/attributes/surface.py b/city_model_structure/attributes/surface.py index f38a571a..29123816 100644 --- a/city_model_structure/attributes/surface.py +++ b/city_model_structure/attributes/surface.py @@ -211,17 +211,17 @@ class Surface: # 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. Area = 0\n') + 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 += self.angle_between_vectors(vec_1, vec_2) + alpha += GeometryHelper.angle_between_vectors(vec_1, vec_2) if alpha == 0: - sys.stderr.write('Warning: the area of a line or point cannot be calculated. Area = 0\n') + 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(self) + horizontal_points = self.rotate_surface_to_horizontal area = 0 for i in range(0, len(horizontal_points)-1): point = horizontal_points[i] @@ -233,48 +233,6 @@ 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 - horizontal_points = [] - x = normal_vector[0] - y = normal_vector[1] - - if x == 0 and y == 0: - # Already horizontal - for point in surface.points: - horizontal_points.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: - sys.stderr.write('Warning: rotation base matrix returned None\n') - else: - for point in surface.points: - new_point = np.matmul(rotation_base_matrix, point) - horizontal_points.append(new_point) - return horizontal_points - - @staticmethod - def angle_between_vectors(vec_1, vec_2): - if np.linalg.norm(vec_1) == 0 or np.linalg.norm(vec_2) == 0: - sys.stderr.write("Warning: impossible to calculate angle between planes' normal. Return 0\n") - return 0 - alpha = math.acos(np.dot(vec_1, vec_2) / np.linalg.norm(vec_1) / np.linalg.norm(vec_2)) - return alpha - def _is_almost_same_terrain(self, terrain_points, ground_points): equal = 0 for terrain_point in terrain_points: @@ -320,7 +278,35 @@ class Surface: """ if self._normal is None: points = self.points - cross_product = np.cross(points[1] - points[0], points[2] - points[0]) + accepted_error = 0.01 + 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]) + 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 + else: + clockwise = 0 + counter_clockwise = 1 + 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]) + 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 + else: + counter_clockwise += 1 + if clockwise < counter_clockwise: + cross_product = np.cross(points[0] - points[len(points) - 2], + points[len(points) - 1] - points[len(points) - 2]) + else: + cross_product = np.cross(points[len(points) - 1] - points[len(points) - 2], + points[0] - points[len(points) - 2]) self._normal = cross_product / np.linalg.norm(cross_product) return self._normal @@ -495,3 +481,38 @@ class Surface: self._is_planar = False break return self._is_planar + + @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 = GeometryHelper.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: + 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 + diff --git a/helpers/geometry_helper.py b/helpers/geometry_helper.py index 589f1c31..c1772ea9 100644 --- a/helpers/geometry_helper.py +++ b/helpers/geometry_helper.py @@ -4,6 +4,7 @@ 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 math import numpy as np import open3d as o3d @@ -18,7 +19,7 @@ class GeometryHelper: Geometry helper class """ - def __init__(self, delta=0.5, area_delta=0.5): + def __init__(self, delta=0, area_delta=0): self._delta = delta self._area_delta = area_delta @@ -252,3 +253,12 @@ class GeometryHelper: power += math.pow(vertex2[dimension]-vertex1[dimension], 2) distance = math.sqrt(power) return distance + + @staticmethod + def angle_between_vectors(vec_1, vec_2): + if np.linalg.norm(vec_1) == 0 or np.linalg.norm(vec_2) == 0: + sys.stderr.write("Warning: impossible to calculate angle between planes' normal. Return 0\n") + return 0 + alpha = math.acos(np.dot(vec_1, vec_2) / np.linalg.norm(vec_1) / np.linalg.norm(vec_2)) + return alpha + diff --git a/tests/test_geometry_factory.py b/tests/test_geometry_factory.py index 081aa3d4..764ca5ac 100644 --- a/tests/test_geometry_factory.py +++ b/tests/test_geometry_factory.py @@ -238,12 +238,10 @@ 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 test_is_planar(self): + 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: - for surface in building.surfaces: - print(surface.type, surface.area) print('volume', building.volume) diff --git a/tests_data/bld100087.gml b/tests_data/bld100087.gml new file mode 100644 index 00000000..ab6abdf4 --- /dev/null +++ b/tests_data/bld100087.gml @@ -0,0 +1,373 @@ + + + + +326011.03601000085 5526048.416990001 -1.6000000000058208 +329466.6600299999 5529018.72205 9.80000000000291 + + + + + +148 + + +m2 + +residential +2019 +4.4 +1 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +327900.0092300009 5527617.35561 3.73359000000346 327903.1820100006 5527620.4350000005 2.8999999999941792 327896.92998999916 5527620.528000001 2.8999999999941792 327900.0092300009 5527617.35561 3.73359000000346 + + + + + + + + + + + + + + + + +327903.6182000004 5527610.50165 4.076969999994617 327904.02800000086 5527606.081 2.8999999999941792 327907.9790000003 5527606.880999999 2.8999999999941792 327903.6182000004 5527610.50165 4.076969999994617 + + + + + + + + + + + + + + + + +327906.6085899994 5527617.37174 3.30688000000373 327908.15699999966 5527618.875 2.8999999999941792 327905.10500000045 5527618.92 2.8999999999941792 327906.6085899994 5527617.37174 3.30688000000373 + + + + + + + + + + + + + + + + +327903.1820100006 5527620.4350000005 0 327903.1820100006 5527620.4350000005 2.8999999999941792 327903.13698999956 5527617.434 2.8999999999941792 327903.13698999956 5527617.434 0 327903.1820100006 5527620.4350000005 0 + + + + + + + + + + + + + + + + +327896.92998999916 5527620.528000001 0 327896.92998999916 5527620.528000001 2.8999999999941792 327903.1820100006 5527620.4350000005 2.8999999999941792 327903.1820100006 5527620.4350000005 0 327896.92998999916 5527620.528000001 0 + + + + + + + + + + + + + + + + +327896.716 5527606.1899999995 0 327896.716 5527606.1899999995 2.8999999999941792 327896.92998999916 5527620.528000001 2.8999999999941792 327896.92998999916 5527620.528000001 0 327896.716 5527606.1899999995 0 + + + + + + + + + + + + + + + + +327904.02800000086 5527606.081 0 327904.02800000086 5527606.081 2.8999999999941792 327896.716 5527606.1899999995 2.8999999999941792 327896.716 5527606.1899999995 0 327904.02800000086 5527606.081 0 + + + + + + + + + + + + + + + + +327907.9790000003 5527606.880999999 0 327907.9790000003 5527606.880999999 2.8999999999941792 327904.02800000086 5527606.081 2.8999999999941792 327904.02800000086 5527606.081 0 327907.9790000003 5527606.880999999 0 + + + + + + + + + + + + + + + + +327908.15699999966 5527618.875 0 327908.15699999966 5527618.875 2.8999999999941792 327907.9790000003 5527606.880999999 2.8999999999941792 327907.9790000003 5527606.880999999 0 327908.15699999966 5527618.875 0 + + + + + + + + + + + + + + + + +327905.10500000045 5527618.92 0 327905.10500000045 5527618.92 2.8999999999941792 327908.15699999966 5527618.875 2.8999999999941792 327908.15699999966 5527618.875 0 327905.10500000045 5527618.92 0 + + + + + + + + + + + + + + + + +327905.08299 5527617.404999999 0 327905.08299 5527617.404999999 2.8999999999941792 327905.10500000045 5527618.92 2.8999999999941792 327905.10500000045 5527618.92 0 327905.08299 5527617.404999999 0 + + + + + + + + + + + + + + + + +327903.13698999956 5527617.434 0 327903.13698999956 5527617.434 2.8999999999941792 327905.08299 5527617.404999999 2.8999999999941792 327905.08299 5527617.404999999 0 327903.13698999956 5527617.434 0 + + + + + + + + + + + + + + + + +327903.13698999956 5527617.434 0 327905.08299 5527617.404999999 0 327905.10500000045 5527618.92 0 327908.15699999966 5527618.875 0 327907.9790000003 5527606.880999999 0 327904.02800000086 5527606.081 0 327896.716 5527606.1899999995 0 327896.92998999916 5527620.528000001 0 327903.1820100006 5527620.4350000005 0 327903.13698999956 5527617.434 0 + + + + + + + + + + + + + + + + +327903.1820100006 5527620.4350000005 2.8999999999941792 327900.0092300009 5527617.35561 3.73359000000346 327899.9643300008 5527614.354800001 3.733559999993304 327903.13698999956 5527617.434 2.8999999999941792 327903.1820100006 5527620.4350000005 2.8999999999941792 + + + + + + + + + + + + + + + + +327896.716 5527606.1899999995 2.8999999999941792 327902.4248399995 5527611.731040001 4.399999999994179 327902.4261399992 5527611.81845 4.399999999994179 327899.9643300008 5527614.354800001 3.733559999993304 327900.0092300009 5527617.35561 3.73359000000346 327896.92998999916 5527620.528000001 2.8999999999941792 327896.716 5527606.1899999995 2.8999999999941792 + + + + + + + + + + + + + + + + +327904.02800000086 5527606.081 2.8999999999941792 327903.6182000004 5527610.50165 4.076969999994617 327902.4248399995 5527611.731040001 4.399999999994179 327896.716 5527606.1899999995 2.8999999999941792 327904.02800000086 5527606.081 2.8999999999941792 + + + + + + + + + + + + + + + + +327908.15699999966 5527618.875 2.8999999999941792 327906.6085899994 5527617.37174 3.30688000000373 327906.58634999953 5527615.85674 3.3068099999945844 327902.4261399992 5527611.81845 4.399999999994179 327902.4248399995 5527611.731040001 4.399999999994179 327903.6182000004 5527610.50165 4.076969999994617 327907.9790000003 5527606.880999999 2.8999999999941792 327908.15699999966 5527618.875 2.8999999999941792 + + + + + + + + + + + + + + + + +327905.08299 5527617.404999999 2.8999999999941792 327906.58634999953 5527615.85674 3.3068099999945844 327906.6085899994 5527617.37174 3.30688000000373 327905.10500000045 5527618.92 2.8999999999941792 327905.08299 5527617.404999999 2.8999999999941792 + + + + + + + + + + + + + + + + +327903.13698999956 5527617.434 2.8999999999941792 327899.9643300008 5527614.354800001 3.733559999993304 327902.4261399992 5527611.81845 4.399999999994179 327906.58634999953 5527615.85674 3.3068099999945844 327905.08299 5527617.404999999 2.8999999999941792 327903.13698999956 5527617.434 2.8999999999941792 + + + + + + + + + + +