From 5bfac5b01b1cfdb9944262aaadaebed27a28a2f7 Mon Sep 17 00:00:00 2001 From: Guille Date: Fri, 27 Nov 2020 11:31:25 -0500 Subject: [PATCH] Partial implementation for non-triangular surfaces --- city_model_structure/attributes/polyhedron.py | 100 ++++++++++++++++-- city_model_structure/city.py | 1 + factories/geometry_feeders/city_gml.py | 29 +++-- .../helpers/occupancy_helper.py | 1 + .../helpers/us_pluto_to_function.py | 1 + helpers/geometry_helper.py | 42 +++++++- tests/test_geometry_factory.py | 2 +- tests/test_idf.py | 2 +- 8 files changed, 151 insertions(+), 27 deletions(-) diff --git a/city_model_structure/attributes/polyhedron.py b/city_model_structure/attributes/polyhedron.py index 0c23875f..463aec50 100644 --- a/city_model_structure/attributes/polyhedron.py +++ b/city_model_structure/attributes/polyhedron.py @@ -4,14 +4,17 @@ SPDX - License - Identifier: LGPL - 3.0 - or -later Copyright © 2020 Project Author Guille Gutierrez guillermo.gutierrezmorote@concordia.ca """ import numpy as np +import matplotlib.pyplot as plt from trimesh import Trimesh from helpers.geometry_helper import GeometryHelper +from city_model_structure.attributes.surface import Surface class Polyhedron: """ Polyhedron class """ + def __init__(self, surfaces): self._surfaces = list(surfaces) self._polygons = [s.polygon for s in surfaces] @@ -52,26 +55,110 @@ class Polyhedron: self._vertices = np.asarray(self._vertices) return self._vertices + @staticmethod + def _point(coordinates): + return coordinates[0], coordinates[1], coordinates[2] + + @staticmethod + def _ground_weird_shape(points, original, triangle=False): + x = [] + y = [] + xo = [] + yo = [] + for o in original: + xo.append(o[0]) + yo.append(o[1]) + for point in points: + x.append(point[0]) + y.append(point[1]) + x = [val - min(xo) for val in x] + y = [val - min(yo) for val in y] + if triangle: + print(len(points)) + for point in range(0, len(x)): + print('point', x[point], y[point]) + + def _triangulate(self, surface): + Polyhedron._ground_weird_shape(surface.points, surface.points) + triangles = [] + triangles_count = len(surface.points) - 2 + points_list = surface.points_list + point_index = 0 + area = surface.area + print(f'area is {area}') + while len(triangles) < triangles_count: + print(f'try vertex {point_index}') + # 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) + if self._geometry.almost_same_area(area, (triangular_surface.area + rest_surface.area)): + area = rest_surface.area + print(f'ok! new area is {area}') + print('Triangle-----------------------------------------') + Polyhedron._ground_weird_shape(triangular_surface.points, surface.points) + print('rest---------------------------------------------') + Polyhedron._ground_weird_shape(rest_surface.points, surface.points) + triangles.append(triangular_surface) + points_list = rest_surface.points_list + if len(rest_surface.points) == 3: + triangles.append(rest_surface) + point_index = 0 + else: + print(f'nok! area {(triangular_surface.area + rest_surface.area)} is not {area}') + print('Triangle-----------------------------------------') + Polyhedron._ground_weird_shape(triangular_surface.points, surface.points) + print('rest---------------------------------------------') + Polyhedron._ground_weird_shape(rest_surface.points, surface.points) + point_index = point_index + 3 + return triangles + @property - def faces(self) -> np.ndarray: + def faces(self) -> [[int]]: """ Polyhedron faces - :return: np.ndarray([int]) + :return: [[int]] """ if self._faces is None: self._faces = [] for surface in self._surfaces: face = [] points = surface.points - for point in points: - face.append(self._position_of(point)) - self._faces.append(face) + if len(points) != 3: # three vertices + print(f'Number of vertex {len(points)}') + sub_surfaces = self._triangulate(surface) + print(f'Transformed {len(points)} vertex into {len(sub_surfaces)} triangles') + print(f'Total area: {surface.area}') + sum_area = 0.0 + for sub_surface in sub_surfaces: + sum_area = sum_area + sub_surface.area + points = sub_surface.points + for point in points: + face.append(self._position_of(point)) + self._faces.append(face) + print(f'Accumulated area: {sum_area}') + else: + for point in points: + face.append(self._position_of(point)) + self._faces.append(face) return self._faces + @property + def _cloud_mesh(self): + if self._mesh is None: + self._mesh = GeometryHelper.create_mesh(self._surfaces) + return self._mesh + @property def _polyhedron_mesh(self): if self._mesh is None: - self._mesh = Trimesh(vertices=self.vertices, faces=self.faces) + try: + print("create mesh") + self._mesh = Trimesh(vertices=self.vertices, faces=self.faces) + except SyntaxError: + self._mesh = self._cloud_mesh return self._mesh @property @@ -82,7 +169,6 @@ class Polyhedron: """ if self._volume is None: if not self._polyhedron_mesh.is_volume: - print('The geometry is not a closed volume') self._volume = np.inf else: self._volume = self._polyhedron_mesh.volume diff --git a/city_model_structure/city.py b/city_model_structure/city.py index 654428e5..3441f82a 100644 --- a/city_model_structure/city.py +++ b/city_model_structure/city.py @@ -142,6 +142,7 @@ class City: for surface in building.surfaces: for surface2 in new_city_object.surfaces: surface.shared(surface2) + self._buildings.append(new_city_object) @property diff --git a/factories/geometry_feeders/city_gml.py b/factories/geometry_feeders/city_gml.py index f9ac05b2..14bf99f2 100644 --- a/factories/geometry_feeders/city_gml.py +++ b/factories/geometry_feeders/city_gml.py @@ -39,7 +39,7 @@ class CityGml: 'http://www.opengis.net/citygml/relief/2.0 http://schemas.opengis.net/citygml/relief/2.0/relief.xsd" ' 'xmlns="http://www.opengis.net/citygml/2.0': None, 'http://www.opengis.net/citygml/2.0': None - }, force_list=('cityObjectMember', 'curveMember', 'boundedBy', 'surfaceMember')) + }, force_list=('cityObjectMember', 'curveMember', 'boundedBy', 'surfaceMember', 'CompositeSurface')) self._city_objects = None self._geometry = GeometryHelper() for bound in self._gml['CityModel']['boundedBy']: @@ -76,7 +76,10 @@ class CityGml: surfaces = [] if 'lod1Solid' in o['Building']: lod += 1 - surfaces = self._lod1(o) + surfaces = CityGml._lod1_solid(o) + elif 'lod1MultiSurface' in o['Building']: + lod += 1 + surfaces = CityGml._lod1_multisurface(o) else: for bound in o['Building']['boundedBy']: surface_type = next(iter(bound)) @@ -92,7 +95,6 @@ class CityGml: terrains = [] if lod_terrain_str in o['Building']: terrains = self._terrains(o, lod_terrain_str) - year_of_construction = None function = None if 'yearOfConstruction' in o['Building']: @@ -117,25 +119,20 @@ class CityGml: terrains.append(curve_points) return terrains - def _lod1(self, o): + @staticmethod + def _lod1_solid(o): try: surfaces = [Surface(s['Polygon']['exterior']['LinearRing']['posList']['#text']) for s in o['Building']['lod1Solid']['Solid']['exterior']['CompositeSurface']['surfaceMember']] except TypeError: surfaces = [Surface(s['Polygon']['exterior']['LinearRing']['posList']) for s in o['Building']['lod1Solid']['Solid']['exterior']['CompositeSurface']['surfaceMember']] - terrains = [] - for s in surfaces: - ground = True - ground_points = [] - for row in s.points: - # all surface vertex are in the lower z - ground_point = [row[0], row[1], self._lower_corner[2]] - ground_points = np.concatenate((ground_points, ground_point), axis=None) - ground = ground and self._geometry.almost_equal(row, ground_point) - if ground: - ground_points = self._geometry.to_points_matrix(ground_points) - terrains.append(ground_points) + return surfaces + + @staticmethod + def _lod1_multisurface(o): + surfaces = [Surface(s['Polygon']['exterior']['LinearRing']['posList']) + for s in o['Building']['lod1MultiSurface']['MultiSurface']['surfaceMember']] return surfaces @staticmethod diff --git a/factories/occupancy_feeders/helpers/occupancy_helper.py b/factories/occupancy_feeders/helpers/occupancy_helper.py index baab5423..12c199fd 100644 --- a/factories/occupancy_feeders/helpers/occupancy_helper.py +++ b/factories/occupancy_feeders/helpers/occupancy_helper.py @@ -8,6 +8,7 @@ Copyright © 2020 Project Author Guille Gutierrez guillermo.gutierrezmorote@conc class OccupancyHelper: occupancy_function = { 'C1': 'C-12 Residential', + 'C5': 'C-12 Residential', 'D3': 'C-12 Residential', 'D6': 'C-12 Residential', 'I1': 'C-2 Health', diff --git a/factories/physics_feeders/helpers/us_pluto_to_function.py b/factories/physics_feeders/helpers/us_pluto_to_function.py index 87352841..4f628b08 100644 --- a/factories/physics_feeders/helpers/us_pluto_to_function.py +++ b/factories/physics_feeders/helpers/us_pluto_to_function.py @@ -240,3 +240,4 @@ class UsPlutoToFunction: :return: str """ return UsPlutoToFunction.building_function[building_pluto_function] + diff --git a/helpers/geometry_helper.py b/helpers/geometry_helper.py index a0df287a..62e0ad78 100644 --- a/helpers/geometry_helper.py +++ b/helpers/geometry_helper.py @@ -4,7 +4,6 @@ SPDX - License - Identifier: LGPL - 3.0 - or -later Copyright © 2020 Project Author Guille Gutierrez guillermo.gutierrezmorote@concordia.ca """ import math - import numpy as np import open3d as o3d import requests @@ -28,11 +27,22 @@ class GeometryHelper: :param location2: :return: Boolean """ + print("check adjacency") max_distance = ConfigurationHelper().max_location_distance_for_shared_walls x = location1[0] - location2[0] y = location1[1] - location2[1] return math.sqrt(math.pow(x, 2) + math.pow(y, 2)) < max_distance + def almost_same_area(self, a1, a2): + """ + Compare two areas and decides if they are almost equal (absolute error under delta) + :param a1 + :param a2 + :return: Boolean + """ + delta = math.fabs(a1 - a2) + return delta <= self._delta *3 + def almost_equal(self, v1, v2): """ Compare two points and decides if they are almost equal (quadratic error under delta) @@ -199,4 +209,32 @@ class GeometryHelper: raise Exception('GET /tasks/ {}'.format(resp.status_code)) else: response = resp.json() - return [response['address']['country_code'], response['address']['city']] + # todo: this is wrong, remove in the future + city = 'new_york_city' + country = 'us' + if 'city' in response['address']: + city = response['address']['city'] + if 'country_code' in response['address']: + country = response['address']['country_code'] + return [country, city] + + @staticmethod + def create_mesh(surfaces): + point_list = [] + # todo: ensure use almost equal in the points to avoid duplicated points + normal_list = [] + for surface in surfaces: + for point in surface.points: + point_list.append(point) + normal_list.append(surface.normal) + pcd = o3d.geometry.PointCloud() + points = np.asarray(point_list) + normals = np.asarray(normal_list) + pcd.points = o3d.utility.Vector3dVector(points) + pcd.normals = o3d.utility.Vector3dVector(normals) + alpha = 0.5 + tetra_mesh, pt_map = o3d.geometry.TetraMesh.create_from_point_cloud(pcd) + mesh = o3d.geometry.TriangleMesh().create_from_point_cloud_alpha_shape(pcd, alpha, tetra_mesh, pt_map) + mesh.compute_vertex_normals() + o3d.visualization.draw_geometries([mesh], mesh_show_back_face=True) + return Trimesh(vertices=np.asarray(mesh.vertices), faces=np.asarray(mesh.triangles)) diff --git a/tests/test_geometry_factory.py b/tests/test_geometry_factory.py index d293b5e7..2b5f2e27 100644 --- a/tests/test_geometry_factory.py +++ b/tests/test_geometry_factory.py @@ -24,7 +24,7 @@ class TestGeometryFactory(TestCase): def _get_citygml(self): if self._city_gml is None: - file_path = (self._example_path / '20buildings.gml').resolve() + file_path = (self._example_path / 'EngHT_Flat_model_lod1.gml').resolve() self._city_gml = GeometryFactory('citygml', file_path).city self.assertIsNotNone(self._city_gml, 'city is none') return self._city_gml diff --git a/tests/test_idf.py b/tests/test_idf.py index 6919048b..6cf40006 100644 --- a/tests/test_idf.py +++ b/tests/test_idf.py @@ -30,7 +30,7 @@ class TestIdf(TestCase): def _get_city(self): if self._city_gml is None: - file_path = (self._example_path / 'buildings.gml').resolve() + file_path = (self._example_path / '20buildings.gml').resolve() self._city_gml = GeometryFactory('citygml', file_path).city PhysicsFactory('us_new_york', self._city_gml) UsageFactory('us_new_york', self._city_gml)