diff --git a/city_model_structure/attributes/surface.py b/city_model_structure/attributes/surface.py index 461bff60..f41dfd61 100644 --- a/city_model_structure/attributes/surface.py +++ b/city_model_structure/attributes/surface.py @@ -219,40 +219,17 @@ class Surface: """ # New method to calculate area if self._area is None: - if len(self.points) < 3: + alpha = 0 + 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 += self.angle_between_vectors(vec_1, vec_2) + if len(self.points) < 3 or alpha == 0: area = 0 + print('Warning: the area of a line or point cannot be calculated. Area = 0') else: - z_vector = [0, 0, 1] - normal_vector = self.normal - points_2d = [] - x = normal_vector[0] - y = normal_vector[1] - - if x != 0 or y != 0: - alpha = math.acos(np.dot(normal_vector, z_vector) / np.linalg.norm(normal_vector) / np.linalg.norm(z_vector)) - turning_line = np.cross(normal_vector, z_vector) - third_axis = np.cross(normal_vector, turning_line) - w_1 = turning_line / np.linalg.norm(turning_line) - w_2 = normal_vector - w_3 = third_axis / np.linalg.norm(third_axis) - # turning_base_matrix - turning_matrix = np.array([[1, 0, 0], - [0, math.cos(alpha), -math.sin(alpha)], - [0, math.sin(alpha), math.cos(alpha)]]) - base_matrix = np.array([w_1, w_2, w_3]) - turning_base_matrix = np.matmul(base_matrix.transpose(), turning_matrix.transpose()) - turning_base_matrix = np.matmul(turning_base_matrix, base_matrix) - - if turning_base_matrix is None: - print('Error processing turning base matrix') - else: - for point in self.points: - new_point = np.matmul(turning_base_matrix, point) - points_2d.append(new_point) - else: - for point in self.points: - points_2d.append([point[0], point[1], 0]) - + # todo: Guilleeeee!!!! + points_2d = self.rotate_surface_to_horizontal(self) polygon_2d = pn.Polygon(np.array(points_2d)) area = 0 for i in range(0, len(polygon_2d.points)-1): @@ -265,6 +242,48 @@ class Surface: self._area = abs(area) return self._area + @staticmethod + def rotate_surface_to_horizontal(surface): + z_vector = [0, 0, 1] + normal_vector = surface.normal + points_2d = [] + x = normal_vector[0] + y = normal_vector[1] + + if x == 0 and y == 0: + # Already horizontal + for point in surface.points: + points_2d.append([point[0], point[1], 0]) + else: + alpha = surface.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, math.cos(alpha), -math.sin(alpha)], + [0, math.sin(alpha), math.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: + print('Error processing rotation base matrix') + else: + for point in surface.points: + new_point = np.matmul(rotation_base_matrix, point) + points_2d.append(new_point) + return points_2d + + @staticmethod + def angle_between_vectors(vec_1, vec_2): + alpha = math.acos(np.dot(vec_1, vec_2) / np.linalg.norm(vec_1) / np.linalg.norm(vec_2)) + if np.linalg.norm(vec_1) == 0 or np.linalg.norm(vec_2) == 0: + print('ERROR') + # todo: rise error + return alpha + def _is_almost_same_terrain(self, terrain_points, ground_points): equal = 0 for terrain_point in terrain_points: