m_osman/city_retrofit
1
0
Fork 0
forked from s_ranjbar/city_retrofit
Code Pull Requests Activity

Partial implementation for non-triangular surfaces

Browse Source
This commit is contained in:
Guille 2020-11-27 11:31:25 -05:00
parent 9f4dc2680d
commit 5bfac5b01b
8 changed files with 151 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

92
city_model_structure/attributes/polyhedron.py
View File

@ -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
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:
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

1
city_model_structure/city.py
View File

@ -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

29
factories/geometry_feeders/city_gml.py
View File

@ -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

1
factories/occupancy_feeders/helpers/occupancy_helper.py
View File

@ -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',

1
factories/physics_feeders/helpers/us_pluto_to_function.py
View File

@ -240,3 +240,4 @@ class UsPlutoToFunction:
:return: str
"""
return UsPlutoToFunction.building_function[building_pluto_function]

42
helpers/geometry_helper.py
View File

@ -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))

2
tests/test_geometry_factory.py
View File

@ -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

2
tests/test_idf.py
View File

@ -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)