223 lines
7.1 KiB
Python
223 lines
7.1 KiB
Python
"""
|
|
Geometry helper
|
|
SPDX - License - Identifier: LGPL - 3.0 - or -later
|
|
Copyright © 2020 Project Author Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
|
|
Contributors Pilar Monsalvete Álvarez de Uribarri pilar.monsalvete@concordia.ca
|
|
"""
|
|
import math
|
|
import numpy as np
|
|
import requests
|
|
from trimesh import Trimesh
|
|
from trimesh import intersections
|
|
from helpers.configuration_helper import ConfigurationHelper
|
|
from city_model_structure.attributes.polygon import Polygon
|
|
from city_model_structure.attributes.polyhedron import Polyhedron
|
|
|
|
|
|
class GeometryHelper:
|
|
"""
|
|
Geometry helper class
|
|
"""
|
|
|
|
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, a1, a2):
|
|
"""
|
|
Compare two areas and decides if they are almost equal (absolute error under delta)
|
|
:param a1
|
|
:param a2
|
|
:return: Boolean
|
|
"""
|
|
if a1 == 0 or a2 == 0:
|
|
return False
|
|
delta = math.fabs(a1 - a2)
|
|
return delta <= self._area_delta
|
|
|
|
def almost_equal(self, delta_max, v1, v2):
|
|
"""
|
|
Compare two points and decides if they are almost equal (distance under delta_max)
|
|
:param delta_max: maximum distance to be considered same point
|
|
:param v1: [x,y,z]
|
|
:param v2: [x,y,z]
|
|
:return: Boolean
|
|
"""
|
|
delta = self.distance_between_points(v1, v2)
|
|
return delta <= delta_max
|
|
|
|
def is_almost_same_surface(self, s1, s2):
|
|
"""
|
|
Compare two surfaces and decides if they are almost equal (quadratic error under delta)
|
|
:param s1: Surface
|
|
:param s2: Surface
|
|
:return: Boolean
|
|
"""
|
|
|
|
# delta is grads an need to be converted into radians
|
|
delta = np.rad2deg(self._delta)
|
|
difference = (s1.inclination - s2.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 = s2.polygon.get_parametric()
|
|
n2 = s2.normal
|
|
for point in s1.points:
|
|
distance = abs(
|
|
(point[0] * parametric[0]) + (point[1] * parametric[1]) + (point[2] * parametric[2]) + parametric[3])
|
|
normal_module = math.sqrt(pow(n2[0], 2) + pow(n2[1], 2) + pow(n2[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 s1.intersect(s2) is None:
|
|
return False
|
|
else:
|
|
return True
|
|
|
|
@staticmethod
|
|
def to_points_matrix(points):
|
|
"""
|
|
Transform a point vector into a point matrix
|
|
:param points: [x, y, z, x, y, z ...]
|
|
:return: [[x,y,z],[x,y,z]...]
|
|
"""
|
|
rows = points.size // 3
|
|
points = points.reshape(rows, 3)
|
|
return points
|
|
|
|
@staticmethod
|
|
def segment_list_to_trimesh(lines) -> 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
|
|
elif 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 gml_surface_to_libs(surface):
|
|
if surface == 'WallSurface':
|
|
return 'Wall'
|
|
elif surface == 'GroundSurface':
|
|
return 'Ground'
|
|
else:
|
|
return 'Roof'
|
|
|
|
@staticmethod
|
|
def get_location(latitude, longitude):
|
|
url = 'https://nominatim.openstreetmap.org/reverse?lat={latitude}&lon={longitude}&format=json'
|
|
resp = requests.get(url.format(latitude=latitude, longitude=longitude))
|
|
if resp.status_code != 200:
|
|
# This means something went wrong.
|
|
raise Exception('GET /tasks/ {}'.format(resp.status_code))
|
|
else:
|
|
response = resp.json()
|
|
# todo: this is wrong, remove in the future
|
|
city = 'new_york_city'
|
|
country = 'us'
|
|
if 'city' in response['address']:
|
|
city = response['address']['city']
|
|
if 'country_code' in response['address']:
|
|
country = response['address']['country_code']
|
|
return [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
|