diff --git a/city_model_structure/attributes/polygon.py b/city_model_structure/attributes/polygon.py index 08584db2..a745b191 100644 --- a/city_model_structure/attributes/polygon.py +++ b/city_model_structure/attributes/polygon.py @@ -15,18 +15,27 @@ class Polygon: Polygon class """ - def __init__(self, vertices): - self._vertices = vertices + def __init__(self, points): self._area = None - self._points = None + self._points = points + self._points_list = None self._normal = None @property def points(self) -> np.ndarray: - if self._points is None: - self._points = self._vertices return self._points + @property + def points_list(self) -> np.ndarray: + """ + Solid surface point coordinates list [x, y, z, x, y, z,...] + :return: np.ndarray + """ + if self._points_list is None: + s = self.points + self._points_list = np.reshape(s, len(s) * 3) + return self._points_list + @property def area(self): """ diff --git a/city_model_structure/attributes/polyhedron.py b/city_model_structure/attributes/polyhedron.py index 655361b0..d9d6a2bb 100644 --- a/city_model_structure/attributes/polyhedron.py +++ b/city_model_structure/attributes/polyhedron.py @@ -9,18 +9,17 @@ import numpy as np from trimesh import Trimesh from helpers.geometry_helper import GeometryHelper from helpers.configuration_helper import ConfigurationHelper -from city_model_structure.attributes.surface import Surface +from city_model_structure.attributes.polygon import Polygon +from helpers.geometry_helper import GeometryHelper as gh class Polyhedron: - # todo: modify class to call polygons better than surfaces """ Polyhedron class """ - def __init__(self, surfaces): - self._surfaces = list(surfaces) - self._polygons = [s.polygon for s in surfaces] + def __init__(self, polygons): + self._polygons = polygons self._polyhedron = None self._triangulated_polyhedron = None self._volume = None @@ -52,7 +51,7 @@ class Polyhedron: """ if self._vertices is None: vertices, self._vertices = [], [] - _ = [vertices.extend(s.points) for s in self._surfaces] + _ = [vertices.extend(s.points) for s in self._polygons] for vertex_1 in vertices: found = False for vertex_2 in self._vertices: @@ -65,18 +64,14 @@ class Polyhedron: self._vertices = np.asarray(self._vertices) return self._vertices - @staticmethod - def _point(coordinates): - return coordinates[0], coordinates[1], coordinates[2] - - def _triangulate(self, surface): + def _triangulate(self, polygon) -> [Polygon]: """ triangulates a polygon following the ear clipping methodology - :param surface: surface + :param polygon: polygon :return: list[triangles] """ - points_list = surface.points_list - normal = surface.normal + points_list = polygon.points_list + normal = polygon.normal # are points concave or convex? total_points_list, concave_points, convex_points = self._starting_lists(points_list, normal) @@ -87,7 +82,7 @@ class Polyhedron: 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])) + rest_points.append(list(polygon.points[p])) if self._is_ear(ear, rest_points): ears.append(ear) point_to_remove = concave_points[i] @@ -115,8 +110,8 @@ class Polyhedron: continue break if len(total_points_list) <= 3 and len(convex_points) > 0: - sys.stderr.write(f'Not able to triangulate surface type {surface.type}\n') - return [surface] + sys.stderr.write(f'Not able to triangulate polygon\n') + return [polygon] last_ear = self._triangle(points_list, total_points_list, concave_points[1]) ears.append(last_ear) return ears @@ -215,16 +210,18 @@ class Polyhedron: """ 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) + points = np.append(previous_point, point) + points = np.append(points, next_point) + points = gh.to_points_matrix(points) + triangle = Polygon(points) error_sum = 0 for i in range(0, len(normal)): - error_sum += triangle.perimeter_surface.normal[i] - normal[i] + error_sum += triangle.normal[i] - normal[i] if np.abs(error_sum) < accepted_error: is_concave = True return is_concave - def _triangle(self, points_list, total_points_list, point_position): + def _triangle(self, points_list, total_points_list, point_position) -> Polygon: """ creates a triangular polygon out of three points :param points_list: points_list @@ -234,10 +231,11 @@ class Polyhedron: """ index = point_position * 3 previous_point_index, next_point_index = self._enveloping_points_indices(point_position, total_points_list) - 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) + points = points_list[previous_point_index:previous_point_index + 3] + points = np.append(points, points_list[index:index + 3]) + points = np.append(points, points_list[next_point_index:next_point_index + 3]) + points = gh.to_points_matrix(points) + triangle = Polygon(points) return triangle @staticmethod @@ -288,11 +286,11 @@ class Polyhedron: def _is_ear(ear, points) -> bool: """ finds whether a triangle is an ear of the polygon - :param ear: surface + :param ear: polygon :param points: [point] :return: boolean """ - area_ear = ear.perimeter_surface.area + area_ear = ear.area for point in points: area_points = 0 point_is_not_vertex = True @@ -303,11 +301,16 @@ class Polyhedron: 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[:]]) + new_points = ear.points[i][:] + new_points = np.append(new_points, ear.points[i + 1][:]) + new_points = np.append(new_points, 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.perimeter_surface.area + new_points = ear.points[i][:] + new_points = np.append(new_points, point[:]) + new_points = np.append(new_points, ear.points[0][:]) + new_points = gh.to_points_matrix(new_points) + new_triangle = Polygon(new_points) + area_points += new_triangle.area if abs(area_points - area_ear) < 1e-6: # point_inside_ear = True return False @@ -324,17 +327,17 @@ class Polyhedron: if self._faces is None: self._faces = [] - for surface in self._surfaces: + for polygon in self._polygons: face = [] - points = surface.points + points = polygon.points if len(points) != 3: - sub_surfaces = self._triangulate(surface) + sub_polygons = self._triangulate(polygon) # todo: I modified this! To be checked @Guille - if len(sub_surfaces) >= 1: - for sub_surface in sub_surfaces: + if len(sub_polygons) >= 1: + for sub_polygon in sub_polygons: face = [] - points = sub_surface.points + points = sub_polygon.points for point in points: face.append(self._position_of(point, face)) self._faces.append(face) @@ -371,8 +374,8 @@ class Polyhedron: """ if self._max_z is None: self._max_z = ConfigurationHelper().min_coordinate - for surface in self._surfaces: - for point in surface.points: + for polygon in self._polygons: + for point in polygon.points: if self._max_z < point[2]: self._max_z = point[2] return self._max_z @@ -385,8 +388,8 @@ class Polyhedron: """ if self._max_y is None: self._max_y = ConfigurationHelper().min_coordinate - for surface in self._surfaces: - for point in surface.points: + for polygon in self._polygons: + for point in polygon.points: if self._max_y < point[1]: self._max_y = point[1] return self._max_y @@ -399,8 +402,8 @@ class Polyhedron: """ if self._max_x is None: self._max_x = ConfigurationHelper().min_coordinate - for surface in self._surfaces: - for point in surface.points: + for polygon in self._polygons: + for point in polygon.points: self._max_x = max(self._max_x, point[0]) return self._max_x @@ -412,8 +415,8 @@ class Polyhedron: """ if self._min_z is None: self._min_z = self.max_z - for surface in self._surfaces: - for point in surface.points: + for polygon in self._polygons: + for point in polygon.points: if self._min_z > point[2]: self._min_z = point[2] return self._min_z @@ -426,8 +429,8 @@ class Polyhedron: """ if self._min_y is None: self._min_y = self.max_y - for surface in self._surfaces: - for point in surface.points: + for polygon in self._polygons: + for point in polygon.points: if self._min_y > point[1]: self._min_y = point[1] return self._min_y @@ -440,8 +443,8 @@ class Polyhedron: """ if self._min_x is None: self._min_x = self.max_x - for surface in self._surfaces: - for point in surface.points: + for polygon in self._polygons: + for point in polygon.points: if self._min_x > point[0]: self._min_x = point[0] return self._min_x diff --git a/city_model_structure/attributes/surface.py b/city_model_structure/attributes/surface.py index eab0d10f..b1be28d8 100644 --- a/city_model_structure/attributes/surface.py +++ b/city_model_structure/attributes/surface.py @@ -27,6 +27,7 @@ class Surface: self._holes_points = None self._holes_points_list = None self._perimeter_points = None + self._perimeter_points_list = None self._azimuth = None self._inclination = None self._area_above_ground = None @@ -37,9 +38,9 @@ class Surface: self._min_z = None self._shared_surfaces = [] self._global_irradiance = dict() - self._perimeter_surface = None - self._hole_surfaces = None - self._solid_surface = None + self._perimeter_polygon = None + self._hole_polygons = None + self._solid_polygons = None def parent(self, parent, surface_id): """ @@ -103,6 +104,37 @@ class Surface: self._holes_points.append(hole_points) return self._holes_points + @property + def perimeter_points(self) -> np.ndarray: + """ + Matrix of points of the perimeter in the same order as in coordinates [[x, y, z],[x, y, z],...] + :return: np.ndarray + """ + if self._perimeter_points is None: + if self.holes_points is None: + self._perimeter_points = self.points + else: + _perimeter_coordinates = self._coordinates + for hole_points in self.holes_points: + _hole = np.append(hole_points, hole_points[0]) + _closed_hole = ' '.join(str(e) for e in [*_hole[:]]) + # add a mark 'M' to ensure that the recombination of points does not provoke errors in finding holes + _perimeter_coordinates = _perimeter_coordinates.replace(_closed_hole, 'M') + _perimeter_coordinates = _perimeter_coordinates.replace('M', '') + self._perimeter_points = np.fromstring(_perimeter_coordinates, dtype=float, sep=' ') + self._perimeter_points = gh.to_points_matrix(self._perimeter_points) + # remove duplicated points + pv = np.array([self._perimeter_points[0]]) + for point in self._perimeter_points: + duplicated_point = False + for p in pv: + if gh().almost_equal(0.0, p, point): + duplicated_point = True + if not duplicated_point: + pv = np.append(pv, [point], axis=0) + self._perimeter_points = pv + return self._perimeter_points + @property def points_list(self) -> np.ndarray: """ @@ -128,6 +160,17 @@ class Surface: np.add(self._holes_points_list, hole_points_list) return self._holes_points_list + @property + def perimeter_points_list(self) -> np.ndarray: + """ + Solid surface point coordinates list [x, y, z, x, y, z,...] + :return: np.ndarray + """ + if self._perimeter_points_list is None: + s = self.perimeter_points + self._perimeter_points_list = np.reshape(s, len(s) * 3) + return self._perimeter_points_list + def _min_coord(self, axis): if axis == 'x': axis = 0 @@ -180,7 +223,7 @@ class Surface: :return: float """ if self._area_above_ground is None: - self._area_above_ground = self.perimeter_surface.area - self.area_below_ground + self._area_above_ground = self.perimeter_polygon.area - self.area_below_ground return self._area_above_ground # todo: to be implemented @@ -199,7 +242,7 @@ class Surface: :return: float """ if self._azimuth is None: - normal = self.perimeter_surface.normal + normal = self.perimeter_polygon.normal self._azimuth = np.arctan2(normal[1], normal[0]) return self._azimuth @@ -210,7 +253,7 @@ class Surface: :return: float """ if self._inclination is None: - self._inclination = np.arccos(self.perimeter_surface.normal[2]) + self._inclination = np.arccos(self.perimeter_polygon.normal[2]) return self._inclination @property @@ -236,7 +279,7 @@ class Surface: :param intersection_area: :return: """ - percent = intersection_area / self.perimeter_surface.area + percent = intersection_area / self.perimeter_polygon.area self._shared_surfaces.append((percent, surface)) # todo reimplement @@ -267,27 +310,27 @@ class Surface: self._global_irradiance = value @property - def perimeter_surface(self) -> Polygon: + def perimeter_polygon(self) -> Polygon: """ total surface defined by the perimeter, merging solid and holes :return: Polygon """ - if self._perimeter_surface is None: - self._perimeter_surface = Polygon(self.perimeter_points) - return self._perimeter_surface + if self._perimeter_polygon is None: + self._perimeter_polygon = Polygon(self.perimeter_points) + return self._perimeter_polygon @property - def solid_surface(self) -> Polygon: + def solid_polygon(self) -> Polygon: """ solid surface :return: Polygon """ - if self._solid_surface is None: - self._solid_surface = Polygon(self.points) - return self._solid_surface + if self._solid_polygons is None: + self._solid_polygons = Polygon(self.points) + return self._solid_polygons @property - def hole_surfaces(self) -> [Polygon]: + def hole_polygons(self) -> [Polygon]: """ hole surfaces, a list of hole polygons found in the surface :return: None, [] or [Polygon] @@ -295,42 +338,11 @@ class Surface: [] -> no holes in the surface [Polygon] -> one or more holes in the surface """ - if self._hole_surfaces is None: + if self._hole_polygons is None: if self.holes_points is None: - self._hole_surfaces = None + self._hole_polygons = None else: - self._hole_surfaces = [] + self._hole_polygons = [] for hole_points in self.holes_points: - self._hole_surfaces.append(Polygon(hole_points)) - return self._hole_surfaces - - @property - def perimeter_points(self) -> np.ndarray: - """ - Matrix of points of the perimeter in the same order as in coordinates [[x, y, z],[x, y, z],...] - :return: np.ndarray - """ - if self._perimeter_points is None: - if self.holes_points is None: - self._perimeter_points = self.points - else: - _perimeter_coordinates = self._coordinates - for hole_points in self.holes_points: - _hole = np.append(hole_points, hole_points[0]) - _closed_hole = ' '.join(str(e) for e in [*_hole[:]]) - # add a mark 'M' to ensure that the recombination of points does not provoke errors in finding holes - _perimeter_coordinates = _perimeter_coordinates.replace(_closed_hole, 'M') - _perimeter_coordinates = _perimeter_coordinates.replace('M', '') - self._perimeter_points = np.fromstring(_perimeter_coordinates, dtype=float, sep=' ') - self._perimeter_points = gh.to_points_matrix(self._perimeter_points) - # remove duplicated points - pv = np.array([self._perimeter_points[0]]) - for point in self._perimeter_points: - duplicated_point = False - for p in pv: - if gh().almost_equal(0.0, p, point): - duplicated_point = True - if not duplicated_point: - pv = np.append(pv, [point], axis=0) - self._perimeter_points = pv - return self._perimeter_points + self._hole_polygons.append(Polygon(hole_points)) + return self._hole_polygons diff --git a/city_model_structure/attributes/thermal_boundary.py b/city_model_structure/attributes/thermal_boundary.py index c4eb90e8..1316d24e 100644 --- a/city_model_structure/attributes/thermal_boundary.py +++ b/city_model_structure/attributes/thermal_boundary.py @@ -19,7 +19,7 @@ class ThermalBoundary: self._surface = surface self._delimits = delimits # ToDo: up to at least LOD2 will be just one thermal opening per Thermal boundary, review for LOD3 and LOD4 - self._thermal_openings = [ThermalOpening] + self._thermal_openings = [ThermalOpening()] self._layers = None self._outside_solar_absorptance = None self._outside_thermal_absorptance = None @@ -83,7 +83,7 @@ class ThermalBoundary: Thermal boundary area in square meters :return: float """ - return self._surface.area + return self._surface.solid_polygon.area @property def area_above_ground(self): diff --git a/city_model_structure/attributes/thermal_zone.py b/city_model_structure/attributes/thermal_zone.py index f686f291..a1345312 100644 --- a/city_model_structure/attributes/thermal_zone.py +++ b/city_model_structure/attributes/thermal_zone.py @@ -57,7 +57,7 @@ class ThermalZone: self._floor_area = 0 for s in self._surfaces: if s.type == 'Ground': - self._floor_area += s.area + self._floor_area += s.perimeter_polygon.area return self._floor_area @property diff --git a/city_model_structure/city_object.py b/city_model_structure/city_object.py index 2810b954..d9f6f488 100644 --- a/city_model_structure/city_object.py +++ b/city_model_structure/city_object.py @@ -19,7 +19,8 @@ class CityObject: self._name = name self._lod = lod self._surfaces = surfaces - self._polyhedron = None + self._detailed_polyhedron = None + self._simplified_polyhedron = None self._geometry = GeometryHelper() @property @@ -37,17 +38,36 @@ class CityObject: City object volume in cubic meters :return: float """ - return self.polyhedron.volume + return self.detailed_polyhedron.volume @property - def polyhedron(self): + def detailed_polyhedron(self): """ - City object polyhedron + City object polyhedron including details such as holes :return: Polyhedron """ - if self._polyhedron is None: - self._polyhedron = Polyhedron(self.surfaces) - return self._polyhedron + if self._detailed_polyhedron is None: + polygons = [] + for surface in self.surfaces: + polygons.append(surface.solid_polygon) + if surface.hole_polygons is not None: + for hole_polygon in surface.hole_polygons: + polygons.append(hole_polygon) + self._detailed_polyhedron = Polyhedron(polygons) + return self._detailed_polyhedron + + @property + def simplified_polyhedron(self): + """ + City object polyhedron, just the simple lod2 representation + :return: Polyhedron + """ + if self._simplified_polyhedron is None: + polygons = [] + for surface in self.surfaces: + polygons.append(surface.perimeter_polygon) + self._simplified_polyhedron = Polyhedron(polygons) + return self._simplified_polyhedron @property def surfaces(self) -> List[Surface]: @@ -74,32 +94,7 @@ class CityObject: City object location :return: [x,y,z] """ - if self._polyhedron is None: - self._polyhedron = Polyhedron(self.surfaces) - return self._polyhedron.center - - def stl_export(self, path): - """ - Export the city object to stl file (city_object_name.stl) to the given path - :param path: str - :return: None - """ - # todo: this is a method just for debugging, it will be soon removed - if self._polyhedron is None: - self._polyhedron = Polyhedron(self.surfaces) - full_path = (Path(path) / (self._name + '.stl')).resolve() - self._polyhedron.stl_export(full_path) - - def obj_export(self, path): - if self._polyhedron is None: - self._polyhedron = Polyhedron(self.surfaces) - full_path = (Path(path) / (self._name + '.obj')).resolve() - self._polyhedron.obj_export(full_path) - - def show(self): - if self._polyhedron is None: - self._polyhedron = Polyhedron(self.surfaces) - self._polyhedron.show() + return self.simplified_polyhedron.center @property def max_height(self): @@ -107,6 +102,4 @@ class CityObject: City object maximal height in meters :return: float """ - if self._polyhedron is None: - self._polyhedron = Polyhedron(self.surfaces) - return self._polyhedron.max_z + return self.simplified_polyhedron.max_z diff --git a/imports/geometry_feeders/city_gml.py b/imports/geometry_feeders/city_gml.py index 2e7561a8..9cff2ac7 100644 --- a/imports/geometry_feeders/city_gml.py +++ b/imports/geometry_feeders/city_gml.py @@ -120,7 +120,7 @@ class CityGml: terrains = [] for curve in curves: curve_points = np.fromstring(curve, dtype=float, sep=' ') - curve_points = self._geometry.to_points_matrix(curve_points, True) + curve_points = self._geometry.to_points_matrix(curve_points) terrains.append(curve_points) return terrains @@ -136,7 +136,7 @@ class CityGml: @staticmethod def _lod1_multi_surface(o): - surfaces = [Surface(CityGml._remove_last_point((s['Polygon']['exterior']['LinearRing']['posList'])) + surfaces = [Surface(CityGml._remove_last_point(s['Polygon']['exterior']['LinearRing']['posList'])) for s in o['Building']['lod1MultiSurface']['MultiSurface']['surfaceMember']] return surfaces @@ -178,4 +178,4 @@ class CityGml: def _remove_last_point(points): array = points.split(' ') res = " " - return res.join(array[0:len(array) -3]) + return res.join(array[0:len(array) - 3]) diff --git a/tests/test_geometry_factory.py b/tests/test_geometry_factory.py index 935e6a5e..97e2901e 100644 --- a/tests/test_geometry_factory.py +++ b/tests/test_geometry_factory.py @@ -85,7 +85,6 @@ class TestGeometryFactory(TestCase): 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) def test_citygml_buildings(self): @@ -123,7 +122,6 @@ class TestGeometryFactory(TestCase): :return: None """ city = self._get_citygml() - print(city) for building in city.buildings: for surface in building.surfaces: self.assertIsNotNone(surface.name, 'surface name is none') @@ -139,8 +137,6 @@ class TestGeometryFactory(TestCase): self.assertIsNotNone(surface.min_x, 'surface min_x is none') self.assertIsNotNone(surface.min_y, 'surface min_y is none') self.assertIsNotNone(surface.min_z, 'surface min_z is none') - print('points', surface.points) - print('coordinates', surface._coordinates) def test_citygml_thermal_zone(self): """ @@ -257,10 +253,10 @@ class TestGeometryFactory(TestCase): print(surface.holes_points) print('perimeter:', surface.perimeter_points) for i in range(0, len(holes_coordinates)): - print(surface.hole_surfaces[i].area) - print('perimeter:', surface.perimeter_surface.area) - print('solid:', surface.solid_surface.area) + print(surface.hole_polygons[i].area) + print('perimeter:', surface.perimeter_polygon.area) + print('solid:', surface.solid_polygon.area) for i in range(0, len(holes_coordinates)): - print(surface.hole_surfaces[i].normal) - print('perimeter:', surface.perimeter_surface.normal) - print('solid:', surface.solid_surface.normal) + print(surface.hole_polygons[i].normal) + print('perimeter:', surface.perimeter_polygon.normal) + print('solid:', surface.solid_polygon.normal) diff --git a/tests_data/kelowna_test_case.pickle b/tests_data/kelowna_test_case.pickle index 4897d569..af0a71b1 100644 Binary files a/tests_data/kelowna_test_case.pickle and b/tests_data/kelowna_test_case.pickle differ