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

code for storeys division cleaned

Browse Source
This commit is contained in:
Pilar 2021-03-31 14:17:53 -04:00
parent 809dc1b3b6
commit 1498731275
5 changed files with 124 additions and 154 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
city_model_structure/attributes/polygon.py
View File

@ -6,8 +6,7 @@ Copyright © 2020 Project Author Pilar Monsalvete Álvarez de Uribarri pilar.mon
import sys
import numpy as np
from helpers.geometry_helper import GeometryHelper as gh
import math
class Polygon:
@ -51,7 +50,7 @@ class Polygon:
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 += gh.angle_between_vectors(vec_1, vec_2)
alpha += self._angle_between_vectors(vec_1, vec_2)
if alpha == 0:
sys.stderr.write('Warning: the area of a line or point cannot be calculated 2. Area = 0\n')
return 0
@ -80,7 +79,7 @@ class Polygon:
for point in self.points:
horizontal_points.append([point[0], point[1], 0])
else:
alpha = gh.angle_between_vectors(normal_vector, z_vector)
alpha = self._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)
@ -117,7 +116,7 @@ class Polygon:
cross_product = np.cross(vector_1, vector_2)
if np.linalg.norm(cross_product) != 0:
cross_product = cross_product / np.linalg.norm(cross_product)
alpha = gh.angle_between_vectors(vector_1, vector_2)
alpha = self._angle_between_vectors(vector_1, vector_2)
else:
# todo modify here
cross_product = [0, 0, 0]
@ -144,7 +143,7 @@ class Polygon:
cross_product_next = np.cross(vector_1, vector_2)
if np.linalg.norm(cross_product_next) != 0:
cross_product_next = cross_product_next / np.linalg.norm(cross_product_next)
alpha = gh.angle_between_vectors(vector_1, vector_2)
alpha = Polygon._angle_between_vectors(vector_1, vector_2)
else:
cross_product_next = [0, 0, 0]
alpha = 0
@ -161,6 +160,10 @@ class Polygon:
triangulates a polygon following the ear clipping methodology
:return: list[triangles]
"""
# todo: review triangulate_polygon in
# https://github.com/mikedh/trimesh/blob/dad11126742e140ef46ba12f8cb8643c83356467/trimesh/creation.py#L415,
# it had a problem with a class called 'triangle', but, if solved,
# it could be a very good substitute of this method
points_list = self.points_list
normal = self.normal
# are points concave or convex?
@ -266,7 +269,8 @@ class Polygon:
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)
rows = points.size // 3
points = points.reshape(rows, 3)
triangle = Polygon(points)
return triangle
@ -340,7 +344,8 @@ class Polygon:
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)
rows = new_points.size // 3
new_points = new_points.reshape(rows, 3)
new_triangle = Polygon(new_points)
area_points += new_triangle.area
if abs(area_points - area_ear) < 1e-6:
@ -401,7 +406,8 @@ class Polygon:
accepted_error = 0.1
points = np.append(previous_point, point)
points = np.append(points, next_point)
points = gh.to_points_matrix(points)
rows = points.size // 3
points = points.reshape(rows, 3)
triangle = Polygon(points)
error_sum = 0
for i in range(0, len(normal)):
@ -409,3 +415,22 @@ class Polygon:
if np.abs(error_sum) < accepted_error:
is_concave = True
return is_concave
@staticmethod
def _angle_between_vectors(vec_1, vec_2):
"""
angle between vectors in radians
:param vec_1: vector
:param vec_2: vector
:return: float
"""
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
cosine = np.dot(vec_1, vec_2) / np.linalg.norm(vec_1) / np.linalg.norm(vec_2)
if cosine > 1 and cosine-1 < 1e-5:
cosine = 1
elif cosine < -1 and cosine+1 > -1e-5:
cosine = -1
alpha = math.acos(cosine)
return alpha

17
city_model_structure/attributes/polyhedron.py
View File

@ -5,10 +5,9 @@ Copyright © 2020 Project Author Guille Gutierrez guillermo.gutierrezmorote@conc
Contributors Pilar Monsalvete pilar_monsalvete@yahoo.es
"""
import numpy as np
import math
from trimesh import Trimesh
from helpers.geometry_helper import GeometryHelper
from helpers.configuration_helper import ConfigurationHelper
from city_model_structure.attributes.polygon import Polygon
class Polyhedron:
@ -31,13 +30,17 @@ class Polyhedron:
self._min_z = None
self._min_y = None
self._min_x = None
self._geometry = GeometryHelper()
def _position_of(self, point, face):
vertices = self.vertices
for i in range(len(vertices)):
# ensure not duplicated vertex
if i not in face and GeometryHelper.distance_between_points(vertices[i], point) == 0:
power = 0
vertex2 = vertices[i]
for dimension in range(0, 3):
power += math.pow(vertex2[dimension] - point[dimension], 2)
distance = math.sqrt(power)
if i not in face and distance == 0:
return i
return -1
@ -54,7 +57,11 @@ class Polyhedron:
found = False
for vertex_2 in self._vertices:
found = False
if GeometryHelper.distance_between_points(vertex_1, vertex_2) == 0:
power = 0
for dimension in range(0, 3):
power += math.pow(vertex_2[dimension] - vertex_1[dimension], 2)
distance = math.sqrt(power)
if distance == 0:
found = True
break
if not found:

105
city_model_structure/building.py
View File

@ -6,8 +6,10 @@ contributors: Pilar Monsalvete pilar_monsalvete@yahoo.es
"""
import sys
from typing import List
import numpy as np
import math
from city_model_structure.attributes.surface import Surface
from city_model_structure.attributes.thermal_boundary import ThermalBoundary
@ -16,12 +18,8 @@ from city_model_structure.attributes.usage_zone import UsageZone
from city_model_structure.city_object import CityObject
from helpers.geometry_helper import GeometryHelper as gh
from helpers.configuration_helper import ConfigurationHelper
import math
from pathlib import Path
from trimesh import intersections
from trimesh import Trimesh
from city_model_structure.attributes.polygon import Polygon
from city_model_structure.attributes.polyhedron import Polyhedron
class Building(CityObject):
"""
@ -52,6 +50,7 @@ class Building(CityObject):
self._min_y = ConfigurationHelper().max_coordinate
self._min_z = ConfigurationHelper().max_coordinate
self._centroid = None
self._eave_height = None
self._grounds = []
self._roofs = []
@ -369,65 +368,50 @@ class Building(CityObject):
return self._building_lower_corner
@property
def storeys(self):
storeys = []
height = self.average_storey_height
if self.storeys_above_ground is not None:
number_of_storeys = self.storeys_above_ground
else:
number_of_storeys = math.floor(float(self.max_height) / height) + 1
print('number_of_storeys', number_of_storeys)
last_storey_height = float(self.max_height) - height*(number_of_storeys-1)
print('last_storey_height', last_storey_height)
if last_storey_height < 0.9*height:
number_of_storeys -= 1
print('number storeys', number_of_storeys)
def eave_height(self):
"""
building eave height in meters
:return: float
"""
if self._eave_height is None:
self._eave_height = 0
for wall in self.walls:
self._eave_height = max(self._eave_height, wall.max_z)
return self._eave_height
@property
def storeys(self) -> [Trimesh]:
"""
subsections of building trimesh by storage in case of no interiors defined
:return: [Trimesh]
"""
trimesh = self.simplified_polyhedron.trimesh
rest_trimesh = trimesh
if self.average_storey_height is None:
if self.storeys_above_ground is None or self.storeys_above_ground <= 0:
sys.stderr.write('Warning: not enough information to divide building into storeys, '
'either number of storeys or average storey height must be provided.\n')
return [trimesh]
else:
number_of_storeys = int(self.storeys_above_ground)
height = self.eave_height / number_of_storeys
else:
height = self.average_storey_height
if self.storeys_above_ground is not None:
number_of_storeys = int(self.storeys_above_ground)
else:
number_of_storeys = math.floor(float(self.eave_height) / height) + 1
last_storey_height = float(self.eave_height) - height*(number_of_storeys-1)
if last_storey_height < 0.3*height:
number_of_storeys -= 1
storeys = []
for n in range(0, number_of_storeys - 1):
print(n)
point_plane = [self.building_lower_corner[0], self.building_lower_corner[1],
self.building_lower_corner[2] + height*(n+1)]
print('point plane', point_plane)
print('rest trimesh', rest_trimesh.volume)
vertices, faces = intersections.slice_faces_plane(rest_trimesh.vertices, rest_trimesh.faces, [0, 0, -1],
point_plane)
lines = list(intersections.mesh_plane(rest_trimesh, [0, 0, -1], point_plane))
line_points = gh.segment_list_to_point_cloud(lines)
polyhedron = Polyhedron(Polygon(line_points).triangulate())
tri = Trimesh(polyhedron.vertices, polyhedron.faces)
new_faces = []
for face in tri.faces:
new_face = []
for point in face:
point += len(vertices)
new_face.append(point)
new_faces.append(new_face)
vertices = np.append(vertices, tri.vertices, axis=0)
faces = np.append(faces, new_faces, axis=0)
storey = Trimesh(vertices, faces)
file_name = 'storey_' + str(n) + '.obj'
path_name = (Path(__file__).parent.parent / 'tests' / 'tests_outputs' / file_name).resolve()
with open(path_name, 'w') as file:
file.write(storey.export(file_type='obj'))
vertices, faces = intersections.slice_faces_plane(rest_trimesh.vertices, rest_trimesh.faces, [0, 0, 1],
point_plane)
new_faces = []
for face in tri.faces:
new_face = []
for point in face:
point += len(vertices)
new_face.append(point)
new_faces.append(new_face)
vertices = np.append(vertices, tri.vertices, axis=0)
faces = np.append(faces, new_faces, axis=0)
rest_trimesh = Trimesh(vertices, faces)
file_name = 'rest_trimesh_' + str(n) + '.obj'
path_name = (Path(__file__).parent.parent / 'tests' / 'tests_outputs' / file_name).resolve()
with open(path_name, 'w') as file:
file.write(rest_trimesh.export(file_type='obj'))
self.building_lower_corner[2] + height * (n + 1)]
normal = [0, 0, -1]
storey, trimesh = gh.divide_mesh_by_plane(trimesh, normal, point_plane)
storeys.append(storey)
storeys.append(rest_trimesh)
storeys.append(trimesh)
return storeys
@property
@ -442,7 +426,6 @@ class Building(CityObject):
if 355 > grads > 5:
self._roof_type = 'pitch'
break
print (self._roof_type)
return self._roof_type
@property

108
helpers/geometry_helper.py
View File

@ -7,11 +7,12 @@ Contributors Pilar Monsalvete Álvarez de Uribarri pilar.monsalvete@concordia.ca
import sys
import math
import numpy as np
import open3d as o3d
import requests
from trimesh import Trimesh
from trimesh import intersections
from helpers.configuration_helper import ConfigurationHelper
from city_model_structure.attributes.polygon import Polygon
from city_model_structure.attributes.polyhedron import Polyhedron
class GeometryHelper:
@ -23,10 +24,6 @@ class GeometryHelper:
self._delta = delta
self._area_delta = area_delta
@property
def config(self):
print(f'delta {self._delta} area {self._area_delta}')
@staticmethod
def adjacent_locations(location1, location2):
"""
@ -110,7 +107,7 @@ class GeometryHelper:
return points
@staticmethod
def segment_list_to_point_cloud(lines):
def segment_list_to_trimesh(lines) -> Trimesh:
line_points = [lines[0][0], lines[0][1]]
lines.remove(lines[0])
while len(lines) > 1:
@ -125,22 +122,9 @@ class GeometryHelper:
line_points.append(line[0])
lines.pop(i - 1)
break
return line_points
@staticmethod
def _point_cloud_to_mesh(point_list, normal_list):
# todo @Guille: I think this method should be removed (and create_mesh??)
# Return a mesh composed only by triangles
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(point_list)
pcd.normals = o3d.utility.Vector3dVector(normal_list)
distances = pcd.compute_nearest_neighbor_distance()
avg_dist = np.mean(distances)
radius = 3 * avg_dist
bpa_mesh = o3d.geometry.TriangleMesh().create_from_point_cloud_ball_pivoting(
pcd, o3d.utility.DoubleVector([radius, radius * 2]))
mesh_result = Trimesh(vertices=np.asarray(bpa_mesh.vertices), faces=np.asarray(bpa_mesh.triangles))
return mesh_result
polyhedron = Polyhedron(Polygon(line_points).triangulate())
trimesh = Trimesh(polyhedron.vertices, polyhedron.faces)
return trimesh
@staticmethod
def _merge_meshes(mesh1, mesh2):
@ -160,14 +144,15 @@ class GeometryHelper:
f_merge = np.concatenate((f_merge, [surface]))
mesh_merge = Trimesh(vertices=v_merge, faces=f_merge)
mesh_merge.fix_normals()
return mesh_merge
@staticmethod
def divide_mesh_by_plane(mesh, normal_plane, point_plane):
def divide_mesh_by_plane(trimesh, normal_plane, point_plane):
"""
Divide a mesh by a plane
:param mesh: Trimesh
:param trimesh: Trimesh
:param normal_plane: [x, y, z]
:param point_plane: [x, y, z]
:return: [Trimesh]
@ -175,28 +160,26 @@ class GeometryHelper:
# The first mesh returns the positive side of the plane and the second the negative side.
# If the plane does not divide the mesh (i.e. it does not touch it or it is coplanar with one or more faces),
# then it returns only the original mesh.
# todo: review split method in https://github.com/mikedh/trimesh/issues/235,
# once triangulate_polygon in Polygon class is solved
normal_plane_opp = [None] * len(normal_plane)
for i in range(0, len(normal_plane)):
normal_plane_opp[i] = - normal_plane[i]
normal = [normal_plane, normal_plane_opp]
normal_opp = [normal_plane_opp, normal_plane]
mesh_final = []
for i in range(0, 2):
mesh_1 = intersections.slice_mesh_plane(mesh, normal[i], point_plane)
mesh_1_segments = intersections.mesh_plane(mesh, normal[i], point_plane)
mesh.difference(mesh_1, engine='blender')
if len(mesh_1_segments) <= 0 or len(mesh_1.faces) == len(mesh.faces):
mesh_final.append(mesh)
break
else:
points = GeometryHelper.segment_list_to_point_cloud(mesh_1_segments)
points_normals = [[None] * 3] * len(points)
for j in range(0, len(points_normals)):
points_normals[j] = normal_opp[i]
mesh_2 = GeometryHelper._point_cloud_to_mesh(points, points_normals)
mesh_final.append(GeometryHelper._merge_meshes(mesh_1, mesh_2))
return mesh_final
section_1 = intersections.slice_mesh_plane(trimesh, normal_plane, point_plane)
if section_1 is None:
return [trimesh]
lines = list(intersections.mesh_plane(trimesh, normal_plane, point_plane))
cap = GeometryHelper.segment_list_to_trimesh(lines)
trimesh_1 = GeometryHelper._merge_meshes(section_1, cap)
section_2 = intersections.slice_mesh_plane(trimesh, normal_plane_opp, point_plane)
if section_2 is None:
return [trimesh_1]
trimesh_2 = GeometryHelper._merge_meshes(section_2, cap)
return [trimesh_1, trimesh_2]
@staticmethod
def gml_surface_to_libs(surface):
@ -225,45 +208,16 @@ class GeometryHelper:
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))
@staticmethod
def distance_between_points(vertex1, vertex2):
"""
distance between points in an n-D Euclidean space
:param vertex1: point or vertex
:param vertex2: point or vertex
:return: float
"""
power = 0
for dimension in range(0, len(vertex1)):
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
cosine = np.dot(vec_1, vec_2) / np.linalg.norm(vec_1) / np.linalg.norm(vec_2)
if cosine > 1 and cosine-1 < 1e-5:
cosine = 1
elif cosine < -1 and cosine+1 > -1e-5:
cosine = -1
alpha = math.acos(cosine)
return alpha

5
tests/test_building.py
View File

@ -30,8 +30,9 @@ class MyTestCase(TestCase):
def test_storeys_division(self):
file = 'kelowna.gml'
city = self._get_citygml(file)
i = 0
for building in city.buildings:
if building.name == 'BLD126221':
if i < 5:
building.average_storey_height = 1.5
print(building.name)
print(building.volume)
@ -39,4 +40,4 @@ class MyTestCase(TestCase):
print(building.max_height)
print(building.centroid)
print(len(building.storeys))
i += 1