From 4c7eea5524177fc391cb672ce2c12bd182738558 Mon Sep 17 00:00:00 2001 From: Guille Date: Mon, 21 Dec 2020 09:42:54 -0500 Subject: [PATCH] Correct surface triangulation algorithms --- city_model_structure/attributes/polyhedron.py | 76 +++++++++++++++---- city_model_structure/attributes/surface.py | 34 +++------ factories/geometry_feeders/city_gml.py | 40 +++++++--- 3 files changed, 101 insertions(+), 49 deletions(-) diff --git a/city_model_structure/attributes/polyhedron.py b/city_model_structure/attributes/polyhedron.py index a8ed3d42..7996b106 100644 --- a/city_model_structure/attributes/polyhedron.py +++ b/city_model_structure/attributes/polyhedron.py @@ -69,6 +69,46 @@ class Polyhedron: for point in surface.points: print(point[0] - self.min_x, point[1] - self.min_y, point[2] - self.min_z) + @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:]]) + if len(triangle_right) < 3: + print(f'error right!!!! {triangle_right}') + if len(triangle_left) < 3: + print(f'error right!!!! {triangle_left}') + 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 _print(self, surface): + for point in surface.points: + print(point[0]-surface.min_x, point[1]-surface.min_y, point[2]-surface.min_z) + def _triangulate(self, surface): triangles = [] triangles_count = len(surface.points) - 2 @@ -76,33 +116,43 @@ class Polyhedron: point_index = 0 area = surface.area while len(triangles) < triangles_count: - # select a triangle starting at point index - triangle_points = ' '.join(str(e) for e in [*points_list[point_index:point_index + 9]]) - # remove the middle vertex from previous triangle - rest_points = ' '.join(str(e) for e in [*points_list[0:point_index+3], *points_list[point_index+6:]]) - triangular_surface = Surface(triangle_points, remove_last=False) - rest_surface = Surface(rest_points, remove_last=False) - t_area = triangular_surface.area - if t_area == 0 and self._geometry.almost_same_area(area, (t_area + rest_surface.area)): - area = rest_surface.area - triangles.append(triangular_surface) - points_list = rest_surface.points_list - if len(rest_surface.points) == 3: - triangles.append(rest_surface) + # 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): + area = left_direction[1].area point_index = 0 + triangles.append(left_direction[0]) + points_list = left_direction[1].points_list + 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 else: point_index = point_index + 3 + if point_index > len(points_list): + raise Exception('not ear found!') + if len(points_list) == 9: + # the rest point's are already a triangle + triangles.append(Surface(points_list, remove_last=False)) return triangles @property def faces(self) -> [[int]]: + """ Polyhedron faces :return: [[int]] """ + if self._faces is None: self._faces = [] + for surface in self._surfaces: + face = [] points = surface.points if len(points) != 3: diff --git a/city_model_structure/attributes/surface.py b/city_model_structure/attributes/surface.py index e3bd1bfc..759091be 100644 --- a/city_model_structure/attributes/surface.py +++ b/city_model_structure/attributes/surface.py @@ -198,19 +198,6 @@ class Surface: self._ground_polygon = None return self._ground_polygon - @property - def area_geoms_class(self): - """ - Surface area in square meters - :return: float - """ - if self._area is None: - try: - self._area = self.polygon.get_area() - except AttributeError: - self._area = 0 - return self._area - @property def area(self): """ @@ -230,15 +217,14 @@ class Surface: if alpha == 0: print('Warning: the area of a line or point cannot be calculated. Area = 0') return 0 - points_2d = self.rotate_surface_to_horizontal(self) - polygon_2d = pn.Polygon(np.array(points_2d)) + horizontal_points = self.rotate_surface_to_horizontal(self) area = 0 - for i in range(0, len(polygon_2d.points)-1): - point = polygon_2d.points[i] - next_point = polygon_2d.points[i+1] + 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 = polygon_2d.points[0] - point = polygon_2d.points[len(polygon_2d.points)-1] + 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 @@ -247,14 +233,14 @@ class Surface: def rotate_surface_to_horizontal(surface): z_vector = [0, 0, 1] normal_vector = surface.normal - points_2d = [] + horizontal_points = [] 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]) + 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) @@ -274,8 +260,8 @@ class Surface: else: for point in surface.points: new_point = np.matmul(rotation_base_matrix, point) - points_2d.append(new_point) - return points_2d + horizontal_points.append(new_point) + return horizontal_points @staticmethod def angle_between_vectors(vec_1, vec_2): diff --git a/factories/geometry_feeders/city_gml.py b/factories/geometry_feeders/city_gml.py index 92028b41..3b3673ea 100644 --- a/factories/geometry_feeders/city_gml.py +++ b/factories/geometry_feeders/city_gml.py @@ -68,6 +68,7 @@ class CityGml: :return: City """ if self._city is None: + # todo: refactor this method to clearly choose the gml type self._city = City(self._lower_corner, self._upper_corner, self._srs_name) i = 0 for o in self._gml['CityModel']['cityObjectMember']: @@ -79,10 +80,11 @@ class CityGml: surfaces = CityGml._lod1_solid(o) elif 'lod1MultiSurface' in o['Building']: lod += 1 - surfaces = CityGml._lod1_multisurface(o) + surfaces = CityGml._lod1_multi_surface(o) elif 'lod2MultiSurface' in o['Building']: + # todo: check if this is a real case or a miss-formed citygml lod = 2 - surfaces = surfaces + CityGml._lod2_multisurface(o) + surfaces = surfaces + CityGml._lod2_solid_multi_surface(o) else: for bound in o['Building']['boundedBy']: surface_type = next(iter(bound)) @@ -133,27 +135,41 @@ class CityGml: return surfaces @staticmethod - def _lod1_multisurface(o): + def _lod1_multi_surface(o): surfaces = [Surface(s['Polygon']['exterior']['LinearRing']['posList']) for s in o['Building']['lod1MultiSurface']['MultiSurface']['surfaceMember']] return surfaces @staticmethod - def _lod2_multisurface(o): + def _lod2_solid_multi_surface(o): surfaces = [Surface(s['Polygon']['exterior']['LinearRing']['posList']) for s in o['Building']['lod2MultiSurface']['MultiSurface']['surfaceMember']] return surfaces + @staticmethod + def _lod2_composite_surface(s): + surfaces = [Surface(sm['Polygon']['exterior']['LinearRing']['posList']) + for sm in s['CompositeSurface']['surfaceMember']] + return surfaces + + @staticmethod + def _lod2_multi_surface(s, surface_type): + # todo: this need to be changed into surface bounded? + try: + surfaces = [Surface(s['Polygon']['exterior']['LinearRing']['posList']['#text'], + surface_type=GeometryHelper.gml_surface_to_libs(surface_type))] + except TypeError: + surfaces = [Surface(s['Polygon']['exterior']['LinearRing']['posList'], + surface_type=GeometryHelper.gml_surface_to_libs(surface_type))] + return surfaces + @staticmethod def _lod2(bound): surfaces = [] for surface_type in iter(bound): - try: - surfaces = [Surface(s['Polygon']['exterior']['LinearRing']['posList']['#text'], - surface_type=GeometryHelper.gml_surface_to_libs(surface_type)) - for s in bound[surface_type]['lod2MultiSurface']['MultiSurface']['surfaceMember']] - except TypeError: - surfaces = [Surface(s['Polygon']['exterior']['LinearRing']['posList'], - surface_type=GeometryHelper.gml_surface_to_libs(surface_type)) - for s in bound[surface_type]['lod2MultiSurface']['MultiSurface']['surfaceMember']] + for s in bound[surface_type]['lod2MultiSurface']['MultiSurface']['surfaceMember']: + if 'CompositeSurface' in s: + surfaces = surfaces + CityGml._lod2_composite_surface(s) + else: + surfaces = surfaces + CityGml._lod2_multi_surface(s, surface_type) return surfaces