diff --git a/city_model_structure/attributes/plane.py b/city_model_structure/attributes/plane.py index d6efe1f3..71e0f6d9 100644 --- a/city_model_structure/attributes/plane.py +++ b/city_model_structure/attributes/plane.py @@ -15,12 +15,13 @@ class Plane: Plane class """ - def __init__(self, origin=None, normal=None): + def __init__(self, origin, normal): # todo: other options to define the plane: # by two lines # by three points self._origin = origin self._normal = normal + self._equation = None self._opposite_normal = None @property @@ -29,8 +30,6 @@ class Plane: Get plane origin point :return: Point """ - if self._origin is None: - raise NotImplementedError return self._origin @property @@ -39,10 +38,36 @@ class Plane: Get plane normal [x, y, z] :return: np.ndarray """ - if self._normal is None: - raise NotImplementedError return self._normal + @property + def equation(self) -> (float, float, float, float): + """ + Get the plane equation components Ax + By + Cz + D = 0 + :return: (A, B, C, D) + """ + if self._equation is None: + + a = self.normal[0] + b = self.normal[1] + c = self.normal[2] + d = ((-1 * self.origin.coordinates[0]) * self.normal[0]) + d += ((-1 * self.origin.coordinates[1]) * self.normal[1]) + d += ((-1 * self.origin.coordinates[2]) * self.normal[2]) + self._equation = (a, b, c, d) + return self._equation + + def distance(self, point): + """ + Distance between the given point and the plane + :return: float + """ + p = point + e = self.equation + denominator = np.abs((p[0] * e[0]) + (p[1] * e[1]) + (p[2] * e[2]) + e[3]) + numerator = np.sqrt((e[0]**2) + (e[1]**2) + (e[2]**2)) + return denominator/numerator + @property def opposite_normal(self): """ diff --git a/city_model_structure/attributes/polygon.py b/city_model_structure/attributes/polygon.py index d10dec16..67ff827f 100644 --- a/city_model_structure/attributes/polygon.py +++ b/city_model_structure/attributes/polygon.py @@ -11,423 +11,491 @@ from typing import List import numpy as np from trimesh import Trimesh import trimesh.intersections -from city_model_structure.attributes.point import Point +from city_model_structure.attributes.plane import Plane +from city_model_structure.attributes.point import Point +import helpers.constants as cte class Polygon: - """ + """ Polygon class """ - def __init__(self, coordinates): + def __init__(self, coordinates): + self._area = None + self._points = None + self._points_list = None + self._normal = None + self._inverse = None + self._edges = None + self._coordinates = coordinates + self._triangles = None + self._vertices = None + self._faces = None + self._plane = None - self._area = None - self._points = None - self._points_list = None - self._normal = None - self._inverse = None - self._edges = None - self._coordinates = coordinates - self._triangles = None - self._vertices = None - self._faces = None - - @property - def points(self) -> List[Point]: - """ + @property + def points(self) -> List[Point]: + """ Get the points belonging to the polygon [[x, y, z],...] :return: [Point] """ - if self._points is None: - self._points = [] - for coordinate in self.coordinates: - self._points.append(Point(coordinate)) - return self._points + if self._points is None: + self._points = [] + for coordinate in self.coordinates: + self._points.append(Point(coordinate)) + return self._points - @property - def coordinates(self) -> List[np.ndarray]: - """ + @property + def plane(self) -> Plane: + """ + Get the polygon plane + :return: Plane + """ + if self._plane is None: + self._plane = Plane(normal=self.normal, origin=self.points[0]) + return self._plane + + @property + def coordinates(self) -> List[np.ndarray]: + """ Get the points in the shape of its coordinates belonging to the polygon [[x, y, z],...] :return: [np.ndarray] """ - return self._coordinates + return self._coordinates - @property - def points_list(self) -> np.ndarray: - """ + @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 + :return: boolean + """ + # fixme: This method doesn't seems to work. + n = len(self.vertices) + angle_sum = 0 + for i in range(0, n): + vector_0 = self.vertices[i] + vector_1 = self.vertices[(i+1) % n] + # set to origin + vector_0[0] = vector_0[0] - point.coordinates[0] + vector_0[1] = vector_0[1] - point.coordinates[1] + vector_0[2] = vector_0[2] - point.coordinates[2] + 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) + + scalar_product = Polygon._scalar_product(vector_0, vector_1) + angle = np.pi/2 + if module != 0: + angle = abs(np.arcsin(scalar_product / module)) + angle_sum += angle + print(angle_sum) + return abs(angle_sum - math.pi*2) < cte.EPSILON + + def contains_polygon(self, polygon): + """ + Determines if the given polygon is contained by the current polygon + :return: boolean + """ + print('contains') + for point in polygon.points: + print(point.coordinates, self.contains_point(point)) + + if not self.contains_point(point): + return False + print('Belong!') + return True + + @property + def points_list(self) -> np.ndarray: + """ Get the solid surface point coordinates list [x, y, z, x, y, z,...] :return: np.ndarray """ - if self._points_list is None: - s = self.coordinates - self._points_list = np.reshape(s, len(s) * 3) - return self._points_list + if self._points_list is None: + s = self.coordinates + self._points_list = np.reshape(s, len(s) * 3) + return self._points_list - @property - def edges(self) -> List[List[Point]]: - """ + @property + def edges(self) -> List[List[Point]]: + """ Get polygon edges list :return: [[Point]] """ - if self._edges is None: - self._edges = [] - for i in range(0, len(self.points) - 1): - point_1 = self.points[i] - point_2 = self.points[i + 1] - self._edges.append([point_1, point_2]) - self._edges.append([self.points[len(self.points) - 1], self.points[0]]) - return self._edges + if self._edges is None: + self._edges = [] + for i in range(0, len(self.points) - 1): + point_1 = self.points[i] + point_2 = self.points[i + 1] + self._edges.append([point_1, point_2]) + self._edges.append([self.points[len(self.points) - 1], self.points[0]]) + return self._edges - @property - def area(self): - """ + @property + def area(self): + """ 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) - return self._area + # 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) + return self._area - @property - def _points_rotated_to_horizontal(self): - """ + @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] + 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 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 + 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: - """ + @property + def normal(self) -> np.ndarray: + """ Get surface normal vector :return: np.ndarray """ - if self._normal is None: - points = self.coordinates - # todo: IF THE FIRST ONE IS 0, START WITH THE NEXT - point_origin = points[len(points) - 2] - vector_1 = points[len(points) - 1] - point_origin - vector_2 = points[0] - point_origin - vector_3 = points[1] - point_origin - 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 = self._angle_between_vectors(vector_1, vector_2) - else: - cross_product = [0, 0, 0] - alpha = 0 - if len(points) == 3: - return cross_product - if np.linalg.norm(cross_product) == 0: - return cross_product - alpha += self._angle(vector_2, vector_3, cross_product) - for i in range(0, len(points) - 4): - vector_1 = points[i + 1] - point_origin - vector_2 = points[i + 2] - point_origin - alpha += self._angle(vector_1, vector_2, cross_product) - vector_1 = points[len(points) - 1] - point_origin - vector_2 = points[0] - point_origin - if alpha < 0: - cross_product = np.cross(vector_2, vector_1) - else: - cross_product = np.cross(vector_1, vector_2) - self._normal = cross_product / np.linalg.norm(cross_product) - return self._normal + if self._normal is None: + points = self.coordinates + # todo: IF THE FIRST ONE IS 0, START WITH THE NEXT + point_origin = points[len(points) - 2] + vector_1 = points[len(points) - 1] - point_origin + vector_2 = points[0] - point_origin + vector_3 = points[1] - point_origin + 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 = self._angle_between_vectors(vector_1, vector_2) + else: + cross_product = [0, 0, 0] + alpha = 0 + if len(points) == 3: + return cross_product + if np.linalg.norm(cross_product) == 0: + return cross_product + alpha += self._angle(vector_2, vector_3, cross_product) + for i in range(0, len(points) - 4): + vector_1 = points[i + 1] - point_origin + vector_2 = points[i + 2] - point_origin + alpha += self._angle(vector_1, vector_2, cross_product) + vector_1 = points[len(points) - 1] - point_origin + vector_2 = points[0] - point_origin + if alpha < 0: + cross_product = np.cross(vector_2, vector_1) + else: + cross_product = np.cross(vector_1, vector_2) + self._normal = cross_product / np.linalg.norm(cross_product) + return self._normal - @staticmethod - def _angle(vector_1, vector_2, cross_product): - """ + @staticmethod + def _angle(vector_1, vector_2, cross_product): + """ alpha angle in radians :param vector_1: [float] :param vector_2: [float] :param cross_product: [float] :return: float """ - accepted_normal_difference = 0.01 - 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 = Polygon._angle_between_vectors(vector_1, vector_2) - else: - cross_product_next = [0, 0, 0] - alpha = 0 - delta_normals = 0 - for j in range(0, 3): - delta_normals += cross_product[j] - cross_product_next[j] - if np.abs(delta_normals) < accepted_normal_difference: - return alpha - return -alpha + accepted_normal_difference = 0.01 + 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 = Polygon._angle_between_vectors(vector_1, vector_2) + else: + cross_product_next = [0, 0, 0] + alpha = 0 + delta_normals = 0 + for j in range(0, 3): + delta_normals += cross_product[j] - cross_product_next[j] + if np.abs(delta_normals) < accepted_normal_difference: + return alpha + return -alpha - def triangulate(self) -> List[Polygon]: - """ + 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!! - 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) + # 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!! + 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 - return self._triangles + # 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 + return self._triangles - @staticmethod - def _starting_lists(points_list, normal) -> [List[float], List[float], List[float]]: - """ + @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 + 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: - """ + @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 + 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): - """ + @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 + 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): - """ + @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 + 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: - """ + @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 + 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]: - """ + @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. @@ -440,32 +508,32 @@ class Polygon: :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 + 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: - """ + @staticmethod + def _point_is_concave(normal, point, previous_point, next_point) -> bool: + """ returns whether a point is concave :param normal: normal :param point: point @@ -473,182 +541,182 @@ class Polygon: :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 + 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): - """ + @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 + 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 - @property - def inverse(self): - """ + @property + def inverse(self): + """ Get the polygon coordinates in reversed order :return: [np.ndarray] """ - if self._inverse is None: - self._inverse = self.coordinates[::-1] - return self._inverse + if self._inverse is None: + self._inverse = self.coordinates[::-1] + return self._inverse - def divide(self, plane): - """ + def divide(self, plane): + """ Divides the polygon in two by a plane :param plane: plane that intersects with self to divide it in two parts (Plane) :return: Polygon, Polygon, [Point] """ - tri_polygons = Trimesh(vertices=self.vertices, faces=self.faces) - intersection = trimesh.intersections.mesh_plane(tri_polygons, plane.normal, plane.origin.coordinates) - polys_1 = trimesh.intersections.slice_mesh_plane(tri_polygons, plane.opposite_normal, plane.origin.coordinates) - polys_2 = trimesh.intersections.slice_mesh_plane(tri_polygons, plane.normal, plane.origin.coordinates) - triangles_1 = [] - for triangle in polys_1.triangles: - triangles_1.append(Polygon(triangle)) - polygon_1 = self._reshape(triangles_1) - triangles_2 = [] - for triangle in polys_2.triangles: - triangles_2.append(Polygon(triangle)) - polygon_2 = self._reshape(triangles_2) - return polygon_1, polygon_2, intersection + tri_polygons = Trimesh(vertices=self.vertices, faces=self.faces) + intersection = trimesh.intersections.mesh_plane(tri_polygons, plane.normal, plane.origin.coordinates) + polys_1 = trimesh.intersections.slice_mesh_plane(tri_polygons, plane.opposite_normal, plane.origin.coordinates) + polys_2 = trimesh.intersections.slice_mesh_plane(tri_polygons, plane.normal, plane.origin.coordinates) + triangles_1 = [] + for triangle in polys_1.triangles: + triangles_1.append(Polygon(triangle)) + polygon_1 = self._reshape(triangles_1) + triangles_2 = [] + for triangle in polys_2.triangles: + triangles_2.append(Polygon(triangle)) + polygon_2 = self._reshape(triangles_2) + return polygon_1, polygon_2, intersection - def _reshape(self, triangles) -> Polygon: - edges_list = [] - for i in range(0, len(triangles)): - for edge in triangles[i].edges: - if not self._edge_in_edges_list(edge, edges_list): - edges_list.append(edge) - else: - edges_list = self._remove_from_list(edge, edges_list) - points = self._order_points(edges_list) - return Polygon(points) + def _reshape(self, triangles) -> Polygon: + edges_list = [] + for i in range(0, len(triangles)): + for edge in triangles[i].edges: + if not self._edge_in_edges_list(edge, edges_list): + edges_list.append(edge) + else: + edges_list = self._remove_from_list(edge, edges_list) + points = self._order_points(edges_list) + return Polygon(points) - @staticmethod - def _edge_in_edges_list(edge, edges_list): - for edge_element in edges_list: - if (edge_element[0].distance_to_point(edge[0]) == 0 and edge_element[1].distance_to_point(edge[1]) == 0) or \ - (edge_element[1].distance_to_point(edge[0]) == 0 and edge_element[0].distance_to_point( - edge[1]) == 0): - return True - return False + @staticmethod + def _edge_in_edges_list(edge, edges_list): + for edge_element in edges_list: + if (edge_element[0].distance_to_point(edge[0]) == 0 and edge_element[1].distance_to_point(edge[1]) == 0) or \ + (edge_element[1].distance_to_point(edge[0]) == 0 and edge_element[0].distance_to_point( + edge[1]) == 0): + return True + return False - @staticmethod - def _order_points(edges_list): - # todo: not sure that this method works for any case -> RECHECK - points = edges_list[0] - for _ in range(0, len(points)): - for i in range(1, len(edges_list)): - point_1 = edges_list[i][0] - point_2 = points[len(points) - 1] - if point_1.distance_to_point(point_2) == 0: - points.append(edges_list[i][1]) - points.remove(points[len(points) - 1]) - array_points = [] - for point in points: - array_points.append(point.coordinates) - return np.array(array_points) + @staticmethod + def _order_points(edges_list): + # todo: not sure that this method works for any case -> RECHECK + points = edges_list[0] + for _ in range(0, len(points)): + for i in range(1, len(edges_list)): + point_1 = edges_list[i][0] + point_2 = points[len(points) - 1] + if point_1.distance_to_point(point_2) == 0: + points.append(edges_list[i][1]) + points.remove(points[len(points) - 1]) + array_points = [] + for point in points: + array_points.append(point.coordinates) + return np.array(array_points) - @staticmethod - def _remove_from_list(edge, edges_list): - new_list = [] - for edge_element in edges_list: - if not ((edge_element[0].distance_to_point(edge[0]) == 0 and edge_element[1].distance_to_point( - edge[1]) == 0) or - (edge_element[1].distance_to_point(edge[0]) == 0 and edge_element[0].distance_to_point( - edge[1]) == 0)): - new_list.append(edge_element) - return new_list + @staticmethod + def _remove_from_list(edge, edges_list): + new_list = [] + for edge_element in edges_list: + if not ((edge_element[0].distance_to_point(edge[0]) == 0 and edge_element[1].distance_to_point( + edge[1]) == 0) or + (edge_element[1].distance_to_point(edge[0]) == 0 and edge_element[0].distance_to_point( + edge[1]) == 0)): + new_list.append(edge_element) + return new_list - @property - def vertices(self) -> np.ndarray: - """ + @property + def vertices(self) -> np.ndarray: + """ Get polyhedron vertices :return: np.ndarray(int) """ - if self._vertices is None: - vertices, self._vertices = [], [] - _ = [vertices.extend(s.coordinates) for s in self.triangulate()] - for vertex_1 in vertices: - found = False - for vertex_2 in self._vertices: - found = False - 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: - self._vertices.append(vertex_1) - self._vertices = np.asarray(self._vertices) - return self._vertices + if self._vertices is None: + vertices, self._vertices = [], [] + _ = [vertices.extend(s.coordinates) for s in self.triangulate()] + for vertex_1 in vertices: + found = False + for vertex_2 in self._vertices: + found = False + 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: + self._vertices.append(vertex_1) + self._vertices = np.asarray(self._vertices) + return self._vertices - @property - def faces(self) -> List[List[int]]: - """ + @property + def faces(self) -> List[List[int]]: + """ Get polyhedron triangular faces :return: [face] """ - if self._faces is None: - self._faces = [] + if self._faces is None: + self._faces = [] - for polygon in self.triangulate(): - face = [] - points = polygon.coordinates - if len(points) != 3: - sub_polygons = polygon.triangulate() - # todo: I modified this! To be checked @Guille - if len(sub_polygons) >= 1: - for sub_polygon in sub_polygons: - face = [] - points = sub_polygon.coordinates - for point in points: - face.append(self._position_of(point, face)) - self._faces.append(face) - else: - for point in points: - face.append(self._position_of(point, face)) - self._faces.append(face) - return self._faces + for polygon in self.triangulate(): + face = [] + points = polygon.coordinates + if len(points) != 3: + sub_polygons = polygon.triangulate() + # todo: I modified this! To be checked @Guille + if len(sub_polygons) >= 1: + for sub_polygon in sub_polygons: + face = [] + points = sub_polygon.coordinates + for point in points: + face.append(self._position_of(point, face)) + self._faces.append(face) + else: + for point in points: + face.append(self._position_of(point, face)) + self._faces.append(face) + return self._faces - def _position_of(self, point, face): - """ + def _position_of(self, point, face): + """ position of a specific point in the list of points that define a face :return: int """ - vertices = self.vertices - for i in range(len(vertices)): - # ensure not duplicated vertex - 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 + vertices = self.vertices + for i in range(len(vertices)): + # ensure not duplicated vertex + 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 diff --git a/city_model_structure/attributes/polyhedron.py b/city_model_structure/attributes/polyhedron.py index 220dcdcc..3f292f79 100644 --- a/city_model_structure/attributes/polyhedron.py +++ b/city_model_structure/attributes/polyhedron.py @@ -14,240 +14,240 @@ from helpers.configuration_helper import ConfigurationHelper class Polyhedron: - """ + """ Polyhedron class """ - def __init__(self, polygons): - self._polygons = polygons - self._polyhedron = None - self._triangulated_polyhedron = None - self._volume = None - self._faces = None - self._vertices = None - self._trimesh = None - self._centroid = None - self._max_z = None - self._max_y = None - self._max_x = None - self._min_z = None - self._min_y = None - self._min_x = None + def __init__(self, polygons): + self._polygons = polygons + self._polyhedron = None + self._triangulated_polyhedron = None + self._volume = None + self._faces = None + self._vertices = None + self._trimesh = None + self._centroid = None + self._max_z = None + self._max_y = None + self._max_x = None + self._min_z = None + self._min_y = None + self._min_x = None - def _position_of(self, point, face): - """ + def _position_of(self, point, face): + """ position of a specific point in the list of points that define a face :return: int """ - vertices = self.vertices - for i in range(len(vertices)): - # ensure not duplicated vertex - 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 + vertices = self.vertices + for i in range(len(vertices)): + # ensure not duplicated vertex + 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 - @property - def vertices(self) -> np.ndarray: - """ + @property + def vertices(self) -> np.ndarray: + """ Get polyhedron vertices :return: np.ndarray(int) """ - if self._vertices is None: - vertices, self._vertices = [], [] - _ = [vertices.extend(s.coordinates) for s in self._polygons] - for vertex_1 in vertices: - found = False - for vertex_2 in self._vertices: - found = False - 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: - self._vertices.append(vertex_1) - self._vertices = np.asarray(self._vertices) - return self._vertices + if self._vertices is None: + vertices, self._vertices = [], [] + _ = [vertices.extend(s.coordinates) for s in self._polygons] + for vertex_1 in vertices: + found = False + for vertex_2 in self._vertices: + found = False + 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: + self._vertices.append(vertex_1) + self._vertices = np.asarray(self._vertices) + return self._vertices - @property - def faces(self) -> List[List[int]]: - """ + @property + def faces(self) -> List[List[int]]: + """ Get polyhedron triangular faces :return: [face] """ - if self._faces is None: - self._faces = [] + if self._faces is None: + self._faces = [] - for polygon in self._polygons: + for polygon in self._polygons: - face = [] - points = polygon.coordinates - if len(points) != 3: - sub_polygons = polygon.triangulate() - # todo: I modified this! To be checked @Guille - if len(sub_polygons) >= 1: - for sub_polygon in sub_polygons: - face = [] - points = sub_polygon.coordinates - for point in points: - face.append(self._position_of(point, face)) - self._faces.append(face) - else: - for point in points: - face.append(self._position_of(point, face)) - self._faces.append(face) - return self._faces + face = [] + points = polygon.coordinates + if len(points) != 3: + sub_polygons = polygon.triangulate() + # todo: I modified this! To be checked @Guille + if len(sub_polygons) >= 1: + for sub_polygon in sub_polygons: + face = [] + points = sub_polygon.coordinates + for point in points: + face.append(self._position_of(point, face)) + self._faces.append(face) + else: + for point in points: + face.append(self._position_of(point, face)) + self._faces.append(face) + return self._faces - @property - def trimesh(self) -> Union[Trimesh, None]: - """ + @property + def trimesh(self) -> Union[Trimesh, None]: + """ Get polyhedron trimesh :return: Trimesh """ - if self._trimesh is None: - for face in self.faces: - if len(face) != 3: - sys.stderr.write('Not able to generate trimesh\n') - return None - self._trimesh = Trimesh(vertices=self.vertices, faces=self.faces) - return self._trimesh + if self._trimesh is None: + for face in self.faces: + if len(face) != 3: + sys.stderr.write('Not able to generate trimesh\n') + return None + self._trimesh = Trimesh(vertices=self.vertices, faces=self.faces) + return self._trimesh - @property - def volume(self): - """ + @property + def volume(self): + """ Get polyhedron volume in cubic meters :return: float """ - if self._volume is None: - if self.trimesh is None: - self._volume = np.inf - elif not self.trimesh.is_volume: - self._volume = np.inf - else: - self._volume = self.trimesh.volume - return self._volume + if self._volume is None: + if self.trimesh is None: + self._volume = np.inf + elif not self.trimesh.is_volume: + self._volume = np.inf + else: + self._volume = self.trimesh.volume + return self._volume - @property - def max_z(self): - """ + @property + def max_z(self): + """ Get polyhedron maximal z value in meters :return: float """ - if self._max_z is None: - self._max_z = ConfigurationHelper().min_coordinate - for polygon in self._polygons: - for point in polygon.coordinates: - self._max_z = max(self._max_z, point[2]) - return self._max_z + if self._max_z is None: + self._max_z = ConfigurationHelper().min_coordinate + for polygon in self._polygons: + for point in polygon.coordinates: + self._max_z = max(self._max_z, point[2]) + return self._max_z - @property - def max_y(self): - """ + @property + def max_y(self): + """ Get polyhedron maximal y value in meters :return: float """ - if self._max_y is None: - self._max_y = ConfigurationHelper().min_coordinate - for polygon in self._polygons: - for point in polygon.coordinates: - if self._max_y < point[1]: - self._max_y = point[1] - return self._max_y + if self._max_y is None: + self._max_y = ConfigurationHelper().min_coordinate + for polygon in self._polygons: + for point in polygon.coordinates: + if self._max_y < point[1]: + self._max_y = point[1] + return self._max_y - @property - def max_x(self): - """ + @property + def max_x(self): + """ Get polyhedron maximal x value in meters :return: float """ - if self._max_x is None: - self._max_x = ConfigurationHelper().min_coordinate - for polygon in self._polygons: - for point in polygon.coordinates: - self._max_x = max(self._max_x, point[0]) - return self._max_x + if self._max_x is None: + self._max_x = ConfigurationHelper().min_coordinate + for polygon in self._polygons: + for point in polygon.coordinates: + self._max_x = max(self._max_x, point[0]) + return self._max_x - @property - def min_z(self): - """ + @property + def min_z(self): + """ Get polyhedron minimal z value in meters :return: float """ - if self._min_z is None: - self._min_z = self.max_z - for polygon in self._polygons: - for point in polygon.coordinates: - if self._min_z > point[2]: - self._min_z = point[2] - return self._min_z + if self._min_z is None: + self._min_z = self.max_z + for polygon in self._polygons: + for point in polygon.coordinates: + if self._min_z > point[2]: + self._min_z = point[2] + return self._min_z - @property - def min_y(self): - """ + @property + def min_y(self): + """ Get polyhedron minimal y value in meters :return: float """ - if self._min_y is None: - self._min_y = self.max_y - for polygon in self._polygons: - for point in polygon.coordinates: - if self._min_y > point[1]: - self._min_y = point[1] - return self._min_y + if self._min_y is None: + self._min_y = self.max_y + for polygon in self._polygons: + for point in polygon.coordinates: + if self._min_y > point[1]: + self._min_y = point[1] + return self._min_y - @property - def min_x(self): - """ + @property + def min_x(self): + """ Get polyhedron minimal x value in meters :return: float """ - if self._min_x is None: - self._min_x = self.max_x - for polygon in self._polygons: - for point in polygon.coordinates: - if self._min_x > point[0]: - self._min_x = point[0] - return self._min_x + if self._min_x is None: + self._min_x = self.max_x + for polygon in self._polygons: + for point in polygon.coordinates: + if self._min_x > point[0]: + self._min_x = point[0] + return self._min_x - @property - def centroid(self) -> Union[None, List[float]]: - """ + @property + def centroid(self) -> Union[None, List[float]]: + """ Get polyhedron centroid :return: [x,y,z] """ - if self._centroid is None: - trimesh = self.trimesh - if trimesh is None: - return None - self._centroid = self.trimesh.centroid - return self._centroid + if self._centroid is None: + trimesh = self.trimesh + if trimesh is None: + return None + self._centroid = self.trimesh.centroid + return self._centroid - def stl_export(self, full_path): - """ + def stl_export(self, full_path): + """ Export the polyhedron to stl given file :param full_path: str :return: None """ - self.trimesh.export(full_path, 'stl_ascii') + self.trimesh.export(full_path, 'stl_ascii') - def obj_export(self, full_path): - """ + def obj_export(self, full_path): + """ Export the polyhedron to obj given file :param full_path: str :return: None """ - self.trimesh.export(full_path, 'obj') + self.trimesh.export(full_path, 'obj') - def show(self): - """ + def show(self): + """ Auxiliary function to render the polyhedron :return: None """ - self.trimesh.show() + self.trimesh.show() diff --git a/city_model_structure/building.py b/city_model_structure/building.py index 423d5ab8..951eab05 100644 --- a/city_model_structure/building.py +++ b/city_model_structure/building.py @@ -14,14 +14,14 @@ from city_model_structure.building_demand.usage_zone import UsageZone from city_model_structure.building_demand.storey import Storey from city_model_structure.city_object import CityObject from city_model_structure.building_demand.household import Household +from city_model_structure.attributes.polyhedron import Polyhedron class Building(CityObject): """ Building(CityObject) class """ - def __init__(self, name, lod, surfaces, year_of_construction, function, - city_lower_corner, terrains=None): + def __init__(self, name, lod, surfaces, year_of_construction, function, city_lower_corner, terrains=None): super().__init__(name, lod, surfaces, city_lower_corner) self._households = None self._basement_heated = None @@ -34,6 +34,7 @@ class Building(CityObject): self._floor_area = None self._roof_type = None self._storeys = None + self._geometrical_zones = None self._thermal_zones = [] self._thermal_boundaries = None self._usage_zones = [] @@ -60,6 +61,15 @@ class Building(CityObject): else: self._internal_walls.append(surface) + @property + def geometrical_zones(self) -> List[Polyhedron]: + if self._geometrical_zones is None: + polygons = [] + for surface in self.surfaces: + polygons.append(surface.perimeter_polygon) + self._geometrical_zones = [Polyhedron(polygons)] + return self._geometrical_zones + @property def grounds(self) -> List[Surface]: """ @@ -256,10 +266,6 @@ class Building(CityObject): if value is not None: self._storeys_above_ground = int(value) - @staticmethod - def _tuple_to_point(xy_tuple): - return [xy_tuple[0], xy_tuple[1], 0.0] - @property def heating(self) -> dict: """ diff --git a/city_model_structure/building_demand/surface.py b/city_model_structure/building_demand/surface.py index cd8a6502..1e2178c8 100644 --- a/city_model_structure/building_demand/surface.py +++ b/city_model_structure/building_demand/surface.py @@ -8,7 +8,7 @@ contributors Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca from __future__ import annotations import uuid import numpy as np -from typing import Union +from typing import List, Union from city_model_structure.attributes.polygon import Polygon from city_model_structure.attributes.plane import Plane from city_model_structure.attributes.point import Point @@ -233,8 +233,16 @@ class Surface: """ return self._solid_polygon + @solid_polygon.setter + def solid_polygon(self, value): + """ + Set the solid surface + :return: Polygon + """ + self._solid_polygon = value + @property - def holes_polygons(self) -> Union[Polygon, None]: + def holes_polygons(self) -> Union[List[Polygon], None]: """ Get hole surfaces, a list of hole polygons found in the surface :return: None, [] or [Polygon] @@ -244,6 +252,14 @@ class Surface: """ return self._holes_polygons + @holes_polygons.setter + def holes_polygons(self, value): + """ + Set the hole surfaces + :param value: [Polygon] + """ + self._holes_polygons = value + @property def pv_system_installed(self) -> PvSystem: """ diff --git a/exports/exports_factory.py b/exports/exports_factory.py index 854ff584..6ea6b18b 100644 --- a/exports/exports_factory.py +++ b/exports/exports_factory.py @@ -30,6 +30,10 @@ class ExportsFactory: """ raise NotImplementedError + @property + def _collada(self): + raise NotImplementedError + @property def _energy_ade(self): """ diff --git a/helpers/constants.py b/helpers/constants.py index 00416177..5841af8d 100644 --- a/helpers/constants.py +++ b/helpers/constants.py @@ -87,3 +87,7 @@ HVAC_AVAILABILITY = 'HVAC Avail' INFILTRATION = 'Infiltration' COOLING_SET_POINT = 'ClgSetPt' HEATING_SET_POINT = 'HtgSetPt' + +# Geometry + +EPSILON = 0.0000001 diff --git a/imports/geometry/rhino.py b/imports/geometry/rhino.py index 66f200d8..f9f41338 100644 --- a/imports/geometry/rhino.py +++ b/imports/geometry/rhino.py @@ -3,27 +3,40 @@ Rhino module parses rhino files and import the geometry into the city model stru SPDX - License - Identifier: LGPL - 3.0 - or -later Copyright © 2020 Project Author Guille Gutierrez guillermo.gutierrezmorote@concordia.ca """ +from numpy import inf from rhino3dm import * from rhino3dm._rhino3dm import MeshType +from city_model_structure.attributes.point import Point import numpy as np -import helpers.configuration_helper +from helpers.configuration_helper import ConfigurationHelper from city_model_structure.attributes.polygon import Polygon from city_model_structure.building import Building from city_model_structure.city import City from city_model_structure.building_demand.surface import Surface as LibsSurface +from helpers.constants import EPSILON from imports.geometry.helpers.geometry_helper import GeometryHelper -from helpers.configuration_helper import ConfigurationHelper class Rhino: def __init__(self, path): self._model = File3dm.Read(str(path)) - max_float = 1.7976931348623157e+308 - min_float = -1.7976931348623157e+308 + max_float = float(ConfigurationHelper().max_coordinate) + min_float = float(ConfigurationHelper().min_coordinate) self._min_x = self._min_y = self._min_z = max_float self._max_x = self._max_y = self._max_z = min_float - print('init') + + @staticmethod + def _in_perimeter(wall, corner): + res = wall.contains_point(Point(corner)) + print(f'belong: {res} wall:({wall.coordinates}) corner: ({corner})') + return res + + @staticmethod + def _add_hole(solid_polygon, hole): + first = solid_polygon.points[0] + points = first + hole.points + solid_polygon.points + return Polygon(points) @staticmethod def _solid_points(coordinates) -> np.ndarray: @@ -47,8 +60,10 @@ class Rhino: @property def city(self) -> City: + rhino_objects = [] buildings = [] - print('city') + windows = [] + holes = [] for obj in self._model.Objects: name = obj.Attributes.Id surfaces = [] @@ -56,20 +71,45 @@ class Rhino: if face is None: break _mesh = face.GetMesh(MeshType.Default) + polygon_points = None for i in range(0, len(_mesh.Faces)): faces = _mesh.Faces[i] _points = '' - for index in faces: self._lower_corner(_mesh.Vertices[index].X, _mesh.Vertices[index].Y, _mesh.Vertices[index].Z) _points = _points + f'{_mesh.Vertices[index].X} {_mesh.Vertices[index].Y} {_mesh.Vertices[index].Z} ' polygon_points = Rhino._solid_points(_points.strip()) - print(_points) surfaces.append(LibsSurface(Polygon(polygon_points), Polygon(polygon_points))) - buildings.append(Building(name, 3, surfaces, 'unknown', 'unknown', (self._min_x, self._min_y, self._min_z), [])) + rhino_objects.append(Building(name, 3, surfaces, 'unknown', 'unknown', (self._min_x, self._min_y, self._min_z), [])) lower_corner = (self._min_x, self._min_y, self._min_z) upper_corner = (self._max_x, self._max_y, self._max_z) - city = City(lower_corner, upper_corner, 'Montreal') + city = City(lower_corner, upper_corner, 'EPSG:26918') + for rhino_object in rhino_objects: + if rhino_object.volume is inf: + # is not a building but a window! + for surface in rhino_object.surfaces: + # add to windows the "hole" with the normal inverted + print('window') + windows.append(Polygon(surface.perimeter_polygon.inverse)) + else: + buildings.append(rhino_object) + + # todo: this method will be pretty inefficient + for hole in windows: + corner = hole.coordinates[0] + for building in buildings: + for surface in building.surfaces: + plane = surface.perimeter_polygon.plane + if plane.distance(corner) <= EPSILON: + # todo: this is a hack for dompark project it should not be done this way + # check if the window is in the right high. + if surface.upper_corner[2] >= corner[2] >= surface.lower_corner[2]: + if surface.holes_polygons is None: + surface.holes_polygons = [] + surface.holes_polygons.append(hole) + + + for building in buildings: city.add_city_object(building) return city diff --git a/imports/geometry_factory.py b/imports/geometry_factory.py index da35c334..d24199bf 100644 --- a/imports/geometry_factory.py +++ b/imports/geometry_factory.py @@ -8,7 +8,7 @@ from city_model_structure.city import City from imports.geometry.citygml import CityGml from imports.geometry.obj import Obj from imports.geometry.osm_subway import OsmSubway -#from imports.geometry.rhino import Rhino +from imports.geometry.rhino import Rhino class GeometryFactory: @@ -43,13 +43,13 @@ class GeometryFactory: """ return OsmSubway(self._path).city -# @property -# def _rhino(self) -> City: -# """ -# Enrich the city by using OpenStreetMap information as data source -# :return: City -# """ -# return Rhino(self._path).city + @property + def _rhino(self) -> City: + """ + Enrich the city by using OpenStreetMap information as data source + :return: City + """ + return Rhino(self._path).city @property def city(self) -> City: @@ -59,10 +59,10 @@ class GeometryFactory: """ return getattr(self, self._file_type, lambda: None) -# @property -# def city_debug(self) -> City: -# """ -# Enrich the city given to the class using the class given handler -# :return: City -# """ -# return Rhino(self._path).city + @property + def city_debug(self) -> City: + """ + Enrich the city given to the class using the class given handler + :return: City + """ + return Rhino(self._path).city