106 lines
4.3 KiB
Python
106 lines
4.3 KiB
Python
|
"""
|
||
|
|
export a city into STL format. (Each building is a solid, suitable for RC models such as CityBEM)
|
||
|
|
SPDX - License - Identifier: LGPL - 3.0 - or -later
|
||
|
|
Copyright © 2024 Concordia CERC group
|
||
|
|
Project Coder Saeed Rayegan sr283100@gmail.com
|
||
|
|
"""
|
||
|
|
from pathlib import Path
|
||
|
|
import numpy as np
|
||
|
|
from scipy.spatial import Delaunay
|
||
|
|
|
||
|
|
class Stl:
|
||
|
|
"""
|
||
|
|
Export to stl format
|
||
|
|
"""
|
||
|
|
def __init__(self, city, path):
|
||
|
|
self._city = city
|
||
|
|
self._path = path
|
||
|
|
self._export()
|
||
|
|
|
||
|
|
def _triangulate_stl(self, points_2d, height):
|
||
|
|
#This function requires a set of 2D points for triangulation
|
||
|
|
# Assuming vertices is a NumPy array
|
||
|
|
tri = Delaunay(points_2d)
|
||
|
|
triangles2D = points_2d[tri.simplices]
|
||
|
|
triangles3D = []
|
||
|
|
|
||
|
|
# Iterate through each triangle in triangles2D
|
||
|
|
for triangle in triangles2D:
|
||
|
|
# Extract the existing x and y coordinates
|
||
|
|
x1, y1 = triangle[0]
|
||
|
|
x2, y2 = triangle[1]
|
||
|
|
x3, y3 = triangle[2]
|
||
|
|
|
||
|
|
# Create a 3D point with the specified height
|
||
|
|
point3D=[[x1, height, y1],[x2, height, y2],[x3, height, y3]]
|
||
|
|
|
||
|
|
# Append the 3D points to the triangle list
|
||
|
|
triangles3D.append(point3D)
|
||
|
|
|
||
|
|
return triangles3D
|
||
|
|
|
||
|
|
def _ground(self, coordinate):
|
||
|
|
x = coordinate[0] - self._city.lower_corner[0]
|
||
|
|
y = coordinate[1] - self._city.lower_corner[1]
|
||
|
|
z = coordinate[2] - self._city.lower_corner[2]
|
||
|
|
return x, y, z
|
||
|
|
|
||
|
|
def _to_vertex_stl(self, coordinate):
|
||
|
|
x, y, z = self._ground(coordinate)
|
||
|
|
return [x, z, -y] # Return as a list # to match opengl expectations (check it later)
|
||
|
|
|
||
|
|
def _to_normal_vertex_stl(self, coordinates):
|
||
|
|
ground_vertex = []
|
||
|
|
for coordinate in coordinates:
|
||
|
|
x, y, z = self._ground(coordinate)
|
||
|
|
ground_vertex.append(np.array([x, y, z]))
|
||
|
|
# recalculate the normal to get grounded values
|
||
|
|
edge_1 = ground_vertex[1] - ground_vertex[0]
|
||
|
|
edge_2 = ground_vertex[2] - ground_vertex[0]
|
||
|
|
normal = np.cross(edge_1, edge_2)
|
||
|
|
normal = normal / np.linalg.norm(normal)
|
||
|
|
# Convert normal to list for easier handling in the write operation
|
||
|
|
return normal.tolist()
|
||
|
|
|
||
|
|
|
||
|
|
def _export(self):
|
||
|
|
if self._city.name is None:
|
||
|
|
self._city.name = 'unknown_city'
|
||
|
|
stl_name = f'{self._city.name}.stl'
|
||
|
|
stl_file_path = (Path(self._path).resolve() / stl_name).resolve()
|
||
|
|
with open(stl_file_path, 'w', encoding='utf-8') as stl:
|
||
|
|
for building in self._city.buildings:
|
||
|
|
stl.write(f"solid building{building.name}\n")
|
||
|
|
for surface in building.surfaces:
|
||
|
|
vertices = []
|
||
|
|
normal = self._to_normal_vertex_stl(surface.perimeter_polygon.coordinates) #the normal vector should be calculated for every surface
|
||
|
|
for coordinate in surface.perimeter_polygon.coordinates:
|
||
|
|
vertex = self._to_vertex_stl(coordinate)
|
||
|
|
if vertex not in vertices:
|
||
|
|
vertices.append(vertex)
|
||
|
|
vertices = np.array(vertices)
|
||
|
|
#After collecting the unique vertices of a surface, there is a need to identify if it is located on the roof, floor, or side walls
|
||
|
|
roofStatus=1 #multiplication of the height of all vertices in a surface
|
||
|
|
heightSum=0 #summation of the height of all vertices in a surface
|
||
|
|
for vertex in vertices:
|
||
|
|
roofStatus *= vertex[1]
|
||
|
|
heightSum += vertex[1]
|
||
|
|
if roofStatus>0:
|
||
|
|
#this surface is the roof (first and third elements of vertices should be passed to the triangulation function)
|
||
|
|
triangles=self._triangulate_stl(vertices[:, [0, 2]], vertices[0][1])
|
||
|
|
elif roofStatus==0 and heightSum==0:
|
||
|
|
# this surface is the floor
|
||
|
|
triangles=self._triangulate_stl(vertices[:, [0, 2]], vertices[0][1])
|
||
|
|
elif roofStatus==0 and heightSum>0:
|
||
|
|
# this surface is a vertical wall (no need for triangulation as it can be done manually)
|
||
|
|
triangles = [[vertices[0],vertices[1],vertices[2]], [vertices[2], vertices[3], vertices[0]]]
|
||
|
|
|
||
|
|
# write the facets (triangles) in the stl file
|
||
|
|
for triangle in triangles:
|
||
|
|
stl.write(f"facet normal {normal[0]} {normal[2]} {normal[1]}\n") #following the idea that y axis is the height
|
||
|
|
stl.write(" outer loop\n")
|
||
|
|
for vertex in triangle:
|
||
|
|
stl.write(f" vertex {vertex[0]} {vertex[1]} {vertex[2]}\n")
|
||
|
|
stl.write(" endloop\n")
|
||
|
|
stl.write("endfacet\n")
|
||
|
|
stl.write(f"endsolid building{building.name}\n")
|