CityBEM-CityLayers-SaeedRayegan/hub/exports/formats/stl.py

"""
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:
      with open(stl_file_path.parent/'buildingSurfaceId.txt', 'w') as textfile:
        for building in self._city.buildings:
          stl.write(f"solid building{building.name}\n")
          # write the building id in the buildingSurfaceId.txt file
          textfile.write(building.name + " ")
          for surface in building.surfaces:
            surfaceID=surface.id
            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:
              #write the id of the surface (the triangles that belong to one surface should have similar id). Note that surface id is number like 2,3,4,
              textfile.write(surfaceID + " ")
              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")
          textfile.write("\n") # End the line after all surface IDs are written for a given building
          stl.write(f"endsolid building{building.name}\n")