diff --git a/city_model_structure/attributes/polygon.py b/city_model_structure/attributes/polygon.py index 779aa956..566a401b 100644 --- a/city_model_structure/attributes/polygon.py +++ b/city_model_structure/attributes/polygon.py @@ -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 diff --git a/city_model_structure/attributes/polyhedron.py b/city_model_structure/attributes/polyhedron.py index dc0f48ca..2d033a74 100644 --- a/city_model_structure/attributes/polyhedron.py +++ b/city_model_structure/attributes/polyhedron.py @@ -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: diff --git a/city_model_structure/building.py b/city_model_structure/building.py index 0b7909b8..6b241125 100644 --- a/city_model_structure/building.py +++ b/city_model_structure/building.py @@ -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 diff --git a/helpers/geometry_helper.py b/helpers/geometry_helper.py index 10cf3e11..9bf4973c 100644 --- a/helpers/geometry_helper.py +++ b/helpers/geometry_helper.py @@ -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 - diff --git a/tests/test_building.py b/tests/test_building.py index 14940467..d83b06da 100644 --- a/tests/test_building.py +++ b/tests/test_building.py @@ -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