""" Geometry helper SPDX - License - Identifier: LGPL - 3.0 - or -later Copyright © 2022 Concordia CERC group Project Coder Guille Gutierrez guillermo.gutierrezmorote@concordia.ca Code contributors: Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca """ import math import numpy as np import requests from trimesh import Trimesh from trimesh import intersections from hub.city_model_structure.attributes.polygon import Polygon from hub.city_model_structure.attributes.polyhedron import Polyhedron from hub.helpers.location import Location from hub.helpers.configuration_helper import ConfigurationHelper class GeometryHelper: """ Geometry helper class """ srs_transformations = { 'urn:adv:crs:ETRS89_UTM32*DE_DHHN92_NH': 'epsg:25832' } def __init__(self, delta=0, area_delta=0): self._delta = delta self._area_delta = area_delta @staticmethod def adjacent_locations(location1, location2): """ Determine when two attributes may be adjacent or not based in the dis :param location1: :param location2: :return: Boolean """ 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 """ line_points = [lines[0][0], lines[0][1]] lines.remove(lines[0]) while len(lines) > 1: i = 0 for line in lines: i += 1 if GeometryHelper.distance_between_points(line[0], line_points[len(line_points) - 1]) < 1e-8: line_points.append(line[1]) lines.pop(i - 1) break if GeometryHelper.distance_between_points(line[1], line_points[len(line_points) - 1]) < 1e-8: line_points.append(line[0]) lines.pop(i - 1) break polyhedron = Polyhedron(Polygon(line_points).triangulate()) trimesh = Trimesh(polyhedron.vertices, polyhedron.faces) return trimesh @staticmethod def _merge_meshes(mesh1, mesh2): v_1 = mesh1.vertices f_1 = mesh1.faces v_2 = mesh2.vertices f_2 = mesh2.faces length = len(v_1) v_merge = np.concatenate((v_1, v_2)) f_merge = np.asarray(f_1) for item in f_2: point1 = item.item(0) + length point2 = item.item(1) + length point3 = item.item(2) + length surface = np.asarray([point1, point2, point3]) 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(trimesh, normal_plane, point_plane): """ Divide a mesh by a plane :param trimesh: Trimesh :param normal_plane: [x, y, z] :param point_plane: [x, y, z] :return: [Trimesh] """ # 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] 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 get_location(latitude, longitude) -> Location: """ Get Location from latitude and longitude """ url = 'https://nominatim.openstreetmap.org/reverse?lat={latitude}&lon={longitude}&format=json' response = requests.get(url.format(latitude=latitude, longitude=longitude)) if response.status_code != 200: # This means something went wrong. raise Exception('GET /tasks/ {}'.format(response.status_code)) response = response.json() city = 'Unknown' country = 'ca' if 'city' in response['address']: city = response['address']['city'] if 'country_code' in response['address']: country = response['address']['country_code'] return Location(country, city) @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