modified triangulation method

This commit is contained in:
Pilar 2023-02-10 05:42:57 -05:00
parent b61722db2e
commit a258d33dc9
3 changed files with 79 additions and 414 deletions

View File

@ -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:

View File

@ -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

View File

@ -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:
"""