diff --git a/hub/city_model_structure/attributes/polygon.py b/hub/city_model_structure/attributes/polygon.py index 94d47725..789a2f2b 100644 --- a/hub/city_model_structure/attributes/polygon.py +++ b/hub/city_model_structure/attributes/polygon.py @@ -12,6 +12,9 @@ from typing import List import numpy as np from trimesh import Trimesh import trimesh.intersections +import trimesh.creation +import trimesh.geometry +from shapely.geometry.polygon import Polygon as shapley_polygon from hub.city_model_structure.attributes.plane import Plane from hub.city_model_structure.attributes.point import Point @@ -22,6 +25,7 @@ class Polygon: """ Polygon class """ + # todo: review with @Guille: Points, Coordinates, Vertices, Faces def __init__(self, coordinates): self._area = None @@ -66,20 +70,6 @@ class Polygon: """ return self._coordinates - @staticmethod - def _module(vector): - x2 = vector[0] ** 2 - y2 = vector[1] ** 2 - z2 = vector[2] ** 2 - return math.sqrt(x2+y2+z2) - - @staticmethod - def _scalar_product(vector_0, vector_1): - x = vector_0[0] * vector_1[0] - y = vector_0[1] * vector_1[1] - z = vector_0[2] * vector_1[2] - return x+y+z - def contains_point(self, point): """ Determines if the given point is contained by the current polygon @@ -98,9 +88,9 @@ class Polygon: vector_1[0] = vector_1[0] - point.coordinates[0] vector_1[1] = vector_1[1] - point.coordinates[1] vector_1[2] = vector_1[2] - point.coordinates[2] - module = Polygon._module(vector_0) * Polygon._module(vector_1) + module = np.linalg.norm(vector_0) * np.linalg.norm(vector_1) - scalar_product = Polygon._scalar_product(vector_0, vector_1) + scalar_product = np.dot(vector_0, vector_1) angle = np.pi/2 if module != 0: angle = abs(np.arcsin(scalar_product / module)) @@ -150,69 +140,17 @@ class Polygon: Get surface area in square meters :return: float """ - # New method to calculate area if self._area is None: - if len(self.points) < 3: - sys.stderr.write('Warning: the area of a line or point cannot be calculated 1. Area = 0\n') - return 0 - alpha = 0 - vec_1 = self.points[1].coordinates - self.points[0].coordinates - for i in range(2, len(self.points)): - vec_2 = self.points[i].coordinates - self.points[0].coordinates - 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 - horizontal_points = self._points_rotated_to_horizontal - area = 0 - for i in range(0, len(horizontal_points) - 1): - point = horizontal_points[i] - next_point = horizontal_points[i + 1] - area += (next_point[1] + point[1]) / 2 * (next_point[0] - point[0]) - next_point = horizontal_points[0] - point = horizontal_points[len(horizontal_points) - 1] - area += (next_point[1] + point[1]) / 2 * (next_point[0] - point[0]) - self._area = abs(area) + self._area = 0 + for triangle in self.triangles: + ab = np.zeros(3) + ac = np.zeros(3) + for i in range(0, 3): + ab[i] = triangle.coordinates[1][i] - triangle.coordinates[0][i] + ac[i] = triangle.coordinates[2][i] - triangle.coordinates[0][i] + self._area += np.linalg.norm(np.cross(ab, ac)) / 2 return self._area - @property - def _points_rotated_to_horizontal(self): - """ - polygon points rotated to horizontal - :return: [float] - """ - z_vector = [0, 0, 1] - normal_vector = self.normal - horizontal_points = [] - x = normal_vector[0] - y = normal_vector[1] - - if x == 0 and y == 0: - # Already horizontal - for point in self.points: - horizontal_points.append([point.coordinates[0], point.coordinates[1], 0]) - else: - 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) - w_2 = normal_vector - w_3 = third_axis / np.linalg.norm(third_axis) - rotation_matrix = np.array([[1, 0, 0], - [0, np.cos(alpha), -np.sin(alpha)], - [0, np.sin(alpha), np.cos(alpha)]]) - base_matrix = np.array([w_1, w_2, w_3]) - rotation_base_matrix = np.matmul(base_matrix.transpose(), rotation_matrix.transpose()) - rotation_base_matrix = np.matmul(rotation_base_matrix, base_matrix) - - if rotation_base_matrix is None: - sys.stderr.write('Warning: rotation base matrix returned None\n') - else: - for point in self.points: - new_point = np.matmul(rotation_base_matrix, point.coordinates) - horizontal_points.append(new_point) - return horizontal_points - @property def normal(self) -> np.ndarray: """ @@ -275,284 +213,67 @@ class Polygon: return alpha return -alpha - def triangulate(self) -> List[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 - # this method is very dirty and has an infinite loop solved with a counter!! + @staticmethod + def triangle_mesh(vertices, normal): + min_x = 1e16 + min_y = 1e16 + min_z = 1e16 + for vertex in vertices: + if vertex[0] < min_x: + min_x = vertex[0] + if vertex[1] < min_y: + min_y = vertex[1] + if vertex[2] < min_z: + min_z = vertex[2] + + new_vertices = [] + for vertex in vertices: + vertex = [vertex[0]-min_x, vertex[1]-min_y, vertex[2]-min_z] + new_vertices.append(vertex) + + transformation_matrix = trimesh.geometry.plane_transform(origin=new_vertices[0], normal=normal) + + coordinates = [] + for vertex in vertices: + transformed_vertex = [vertex[0]-min_x, vertex[1]-min_y, vertex[2]-min_z, 1] + transformed_vertex = np.dot(transformation_matrix, transformed_vertex) + coordinate = [transformed_vertex[0], transformed_vertex[1]] + coordinates.append(coordinate) + + polygon = shapley_polygon(coordinates) + + vertices_2d, faces = trimesh.creation.triangulate_polygon(polygon, engine='triangle') + mesh = Trimesh(vertices=vertices, faces=faces) + + # check orientation + normal_sum = 0 + for i in range(0, 3): + normal_sum += normal[i] + mesh.face_normals[0][i] + + if abs(normal_sum) <= 1E-10: + new_faces = [] + for face in faces: + new_face = [] + for i in range(0, len(face)): + new_face.append(face[len(face)-i-1]) + new_faces.append(new_face) + mesh = Trimesh(vertices=vertices, faces=new_faces) + + return mesh + + @property + def triangles(self) -> List[Polygon]: if self._triangles is None: - points_list = self.points_list - normal = self.normal - if np.linalg.norm(normal) == 0: - sys.stderr.write('Not able to triangulate polygon\n') - return [self] - # are points concave or convex? - total_points_list, concave_points, convex_points = self._starting_lists(points_list, normal) - - # list of ears - ears = [] - j = 0 - while (len(concave_points) > 3 or len(convex_points) != 0) and j < 100: - j += 1 - for i in range(0, len(concave_points)): - ear = self._triangle(points_list, total_points_list, concave_points[i]) - rest_points = [] - for points in total_points_list: - rest_points.append(list(self.coordinates[points])) - if self._is_ear(ear, rest_points): - ears.append(ear) - point_to_remove = concave_points[i] - previous_point_in_list, next_point_in_list = self._enveloping_points(point_to_remove, - total_points_list) - total_points_list.remove(point_to_remove) - concave_points.remove(point_to_remove) - # Was any of the adjacent points convex? -> check if changed status to concave - for convex_point in convex_points: - if convex_point == previous_point_in_list: - concave_points, convex_points, end_loop = self._if_concave_change_status(normal, - points_list, - convex_point, - total_points_list, - concave_points, - convex_points, - previous_point_in_list) - if end_loop: - break - continue - if convex_point == next_point_in_list: - concave_points, convex_points, end_loop = self._if_concave_change_status(normal, - points_list, - convex_point, - total_points_list, - concave_points, - convex_points, - next_point_in_list) - if end_loop: - break - continue - break - if len(total_points_list) <= 3 and len(convex_points) > 0: - sys.stderr.write('Not able to triangulate polygon\n') - return [self] - if j >= 100: - sys.stderr.write('Not able to triangulate polygon\n') - return [self] - last_ear = self._triangle(points_list, total_points_list, concave_points[1]) - ears.append(last_ear) - self._triangles = ears + self._triangles = [] + _mesh = self.triangle_mesh(self.coordinates, self.normal) + for face in _mesh.faces: + points = [] + for vertex in face: + points.append(self.coordinates[vertex]) + polygon = Polygon(points) + self._triangles.append(polygon) return self._triangles - @staticmethod - def _starting_lists(points_list, normal) -> [List[float], List[float], List[float]]: - """ - creates the list of vertices (points) that define the polygon (total_points_list), together with other two lists - separating points between convex and concave - :param points_list: points_list - :param normal: normal - :return: list[point], list[point], list[point] - """ - concave_points = [] - convex_points = [] - # lists of concave and convex points - # case 1: first point - point = points_list[0:3] - previous_point = points_list[len(points_list) - 3:] - next_point = points_list[3:6] - index = 0 - total_points_list = [index] - if Polygon._point_is_concave(normal, point, previous_point, next_point): - concave_points.append(index) - else: - convex_points.append(index) - # case 2: all points except first and last - for i in range(0, int((len(points_list) - 6) / 3)): - point = points_list[(i + 1) * 3:(i + 2) * 3] - previous_point = points_list[i * 3:(i + 1) * 3] - next_point = points_list[(i + 2) * 3:(i + 3) * 3] - index = i + 1 - total_points_list.append(index) - if Polygon._point_is_concave(normal, point, previous_point, next_point): - concave_points.append(index) - else: - convex_points.append(index) - # case 3: last point - point = points_list[len(points_list) - 3:] - previous_point = points_list[len(points_list) - 6:len(points_list) - 3] - next_point = points_list[0:3] - index = int(len(points_list) / 3) - 1 - total_points_list.append(index) - if Polygon._point_is_concave(normal, point, previous_point, next_point): - concave_points.append(index) - else: - convex_points.append(index) - return total_points_list, concave_points, convex_points - - @staticmethod - def _triangle(points_list, total_points_list, point_position) -> Polygon: - """ - creates a triangular polygon out of three points - :param points_list: points_list - :param total_points_list: [point] - :param point_position: int - :return: polygon - """ - index = point_position * 3 - previous_point_index, next_point_index = Polygon._enveloping_points_indices(point_position, total_points_list) - 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]) - rows = points.size // 3 - points = points.reshape(rows, 3) - triangle = Polygon(points) - return triangle - - @staticmethod - def _enveloping_points_indices(point_position, total_points_list): - """ - due to the fact that the lists are not circular, a method to find the previous and next points - of an specific one is needed - :param point_position: int - :param total_points_list: [point] - :return: int, int - """ - previous_point_index = None - next_point_index = None - if point_position == total_points_list[0]: - previous_point_index = total_points_list[len(total_points_list) - 1] * 3 - next_point_index = total_points_list[1] * 3 - if point_position == total_points_list[len(total_points_list) - 1]: - previous_point_index = total_points_list[len(total_points_list) - 2] * 3 - next_point_index = total_points_list[0] * 3 - for i in range(1, len(total_points_list) - 1): - if point_position == total_points_list[i]: - previous_point_index = total_points_list[i - 1] * 3 - next_point_index = total_points_list[i + 1] * 3 - return previous_point_index, next_point_index - - @staticmethod - def _enveloping_points(point_to_remove, total_points_list): - """ - due to the fact that the lists are not circular, a method to find the previous and next points - of an specific one is needed - :param point_to_remove: point - :param total_points_list: [point] - :return: point, point - """ - index = total_points_list.index(point_to_remove) - if index == 0: - previous_point_in_list = total_points_list[len(total_points_list) - 1] - next_point_in_list = total_points_list[1] - elif index == len(total_points_list) - 1: - previous_point_in_list = total_points_list[len(total_points_list) - 2] - next_point_in_list = total_points_list[0] - else: - previous_point_in_list = total_points_list[index - 1] - next_point_in_list = total_points_list[index + 1] - return previous_point_in_list, next_point_in_list - - @staticmethod - def _is_ear(ear, points) -> bool: - """ - finds whether a triangle is an ear of the polygon - :param ear: polygon - :param points: [point] - :return: boolean - """ - area_ear = ear.area - for point in points: - area_points = 0 - point_is_not_vertex = True - for i in range(0, 3): - if abs(np.linalg.norm(point) - np.linalg.norm(ear.coordinates[i])) < 0.0001: - point_is_not_vertex = False - break - if point_is_not_vertex: - for i in range(0, 3): - if i != 2: - new_points = ear.coordinates[i][:] - new_points = np.append(new_points, ear.coordinates[i + 1][:]) - new_points = np.append(new_points, point[:]) - else: - new_points = ear.coordinates[i][:] - new_points = np.append(new_points, point[:]) - new_points = np.append(new_points, ear.coordinates[0][:]) - 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: - # point_inside_ear = True - return False - return True - - @staticmethod - def _if_concave_change_status(normal, points_list, convex_point, total_points_list, - concave_points, convex_points, point_in_list) -> [List[float], List[float], bool]: - """ - checks whether an convex specific point change its status to concave after removing one ear in the polygon - returning the new convex and concave points lists together with a flag advising that the list of total points - already 3 and, therefore, the triangulation must be finished. - :param normal: normal - :param points_list: points_list - :param convex_point: int - :param total_points_list: [point] - :param concave_points: [point] - :param convex_points: [point] - :param point_in_list: int - :return: list[points], list[points], boolean - """ - end_loop = False - point = points_list[point_in_list * 3:(point_in_list + 1) * 3] - pointer = total_points_list.index(point_in_list) - 1 - if pointer < 0: - pointer = len(total_points_list) - 1 - previous_point = points_list[total_points_list[pointer] * 3:total_points_list[pointer] * 3 + 3] - pointer = total_points_list.index(point_in_list) + 1 - if pointer >= len(total_points_list): - pointer = 0 - next_point = points_list[total_points_list[pointer] * 3:total_points_list[pointer] * 3 + 3] - if Polygon._point_is_concave(normal, point, previous_point, next_point): - if concave_points[0] > convex_point: - concave_points.insert(0, convex_point) - elif concave_points[len(concave_points) - 1] < convex_point: - concave_points.append(convex_point) - else: - for point_index in range(0, len(concave_points) - 1): - if concave_points[point_index] < convex_point < concave_points[point_index + 1]: - concave_points.insert(point_index + 1, convex_point) - convex_points.remove(convex_point) - end_loop = True - return concave_points, convex_points, end_loop - - @staticmethod - def _point_is_concave(normal, point, previous_point, next_point) -> bool: - """ - returns whether a point is concave - :param normal: normal - :param point: point - :param previous_point: point - :param next_point: point - :return: boolean - """ - is_concave = False - accepted_error = 0.1 - points = np.append(previous_point, point) - points = np.append(points, next_point) - rows = points.size // 3 - points = points.reshape(rows, 3) - triangle = Polygon(points) - error_sum = 0 - for i in range(0, len(normal)): - error_sum += triangle.normal[i] - normal[i] - if np.abs(error_sum) < accepted_error: - is_concave = True - return is_concave - @staticmethod def _angle_between_vectors(vec_1, vec_2): """ @@ -652,12 +373,12 @@ class Polygon: @property def vertices(self) -> np.ndarray: """ - Get polyhedron vertices + Get polygon vertices :return: np.ndarray(int) """ if self._vertices is None: vertices, self._vertices = [], [] - _ = [vertices.extend(s.coordinates) for s in self.triangulate()] + _ = [vertices.extend(s.coordinates) for s in self.triangles] for vertex_1 in vertices: found = False for vertex_2 in self._vertices: @@ -677,17 +398,17 @@ class Polygon: @property def faces(self) -> List[List[int]]: """ - Get polyhedron triangular faces + Get polygon triangular faces :return: [face] """ if self._faces is None: self._faces = [] - for polygon in self.triangulate(): + for polygon in self.triangles: face = [] points = polygon.coordinates if len(points) != 3: - sub_polygons = polygon.triangulate() + sub_polygons = polygon.triangles # todo: I modified this! To be checked @Guille if len(sub_polygons) >= 1: for sub_polygon in sub_polygons: diff --git a/hub/helpers/geometry_helper.py b/hub/helpers/geometry_helper.py index 1c67bf4a..24c202e3 100644 --- a/hub/helpers/geometry_helper.py +++ b/hub/helpers/geometry_helper.py @@ -39,59 +39,12 @@ class GeometryHelper: max_distance = ConfigurationHelper().max_location_distance_for_shared_walls return GeometryHelper.distance_between_points(location1, location2) < max_distance - def almost_same_area(self, area_1, area_2): - """ - Compare two areas and decides if they are almost equal (absolute error under delta) - :param area_1 - :param area_2 - :return: Boolean - """ - if area_1 == 0 or area_2 == 0: - return False - delta = math.fabs(area_1 - area_2) - return delta <= self._area_delta - - def is_almost_same_surface(self, surface_1, surface_2): - """ - Compare two surfaces and decides if they are almost equal (quadratic error under delta) - :param surface_1: Surface - :param surface_2: Surface - :return: Boolean - """ - - # delta is grads an need to be converted into radians - delta = np.rad2deg(self._delta) - difference = (surface_1.inclination - surface_2.inclination) % math.pi - if abs(difference) > delta: - return False - # s1 and s2 are at least almost parallel surfaces - # calculate distance point to plane using all the vertex - # select surface1 value for the point (X,Y,Z) where two of the values are 0 - minimum_distance = self._delta + 1 - parametric = surface_2.polygon.get_parametric() - normal_2 = surface_2.normal - for point in surface_1.points: - distance = abs( - (point[0] * parametric[0]) + (point[1] * parametric[1]) + (point[2] * parametric[2]) + parametric[3]) - normal_module = math.sqrt(pow(normal_2[0], 2) + pow(normal_2[1], 2) + pow(normal_2[2], 2)) - - if normal_module == 0: - continue - distance = distance / normal_module - if distance < minimum_distance: - minimum_distance = distance - if minimum_distance <= self._delta: - break - - if minimum_distance > self._delta or surface_1.intersect(surface_2) is None: - return False - return True - @staticmethod def segment_list_to_trimesh(lines) -> Trimesh: """ Transform a list of segments into a Trimesh """ + # todo: trimesh has a method for this line_points = [lines[0][0], lines[0][1]] lines.remove(lines[0]) while len(lines) > 1: @@ -106,7 +59,7 @@ class GeometryHelper: line_points.append(line[0]) lines.pop(i - 1) break - polyhedron = Polyhedron(Polygon(line_points).triangulate()) + polyhedron = Polyhedron(Polygon(line_points).triangles) trimesh = Trimesh(polyhedron.vertices, polyhedron.faces) return trimesh diff --git a/hub/imports/geometry_factory.py b/hub/imports/geometry_factory.py index b8ba7e68..2c870085 100644 --- a/hub/imports/geometry_factory.py +++ b/hub/imports/geometry_factory.py @@ -9,7 +9,6 @@ import geopandas from hub.city_model_structure.city import City from hub.imports.geometry.citygml import CityGml from hub.imports.geometry.obj import Obj -from hub.imports.geometry.osm_subway import OsmSubway from hub.imports.geometry.rhino import Rhino from hub.imports.geometry.gpandas import GPandas from hub.imports.geometry.geojson import Geojson @@ -83,14 +82,6 @@ class GeometryFactory: self._function_field, self._function_to_hub).city - @property - def _osm_subway(self) -> City: - """ - Enrich the city by using OpenStreetMap information as data source - :return: City - """ - return OsmSubway(self._path).city - @property def _rhino(self) -> City: """