### Install Dependencies with Pip

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/README.md

Installs all project dependencies listed in the requirements.txt file. Ensure you have pip installed and are in the project's root directory.

```bash
pip install -r requirements.txt
```

--------------------------------

### Download and Process 3D BAG Tile

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/README.md

A step-by-step tutorial demonstrating how to download a tile from 3D BAG, run cityStats to compute metrics, and optionally integrate a val3dity report for validation. This provides a practical example of using the tool.

```bash
wget --header='Accept-Encoding: gzip' https://data.3dbag.nl/cityjson/v210908_fd2cee53/3dbag_v210908_fd2cee53_5910.json
```

```bash
python cityStats.py 3dbag_v210908_fd2cee53_5910.json -o 5910.csv
```

```bash
python cityStats.py 3dbag_v210908_fd2cee53_5910.json -v report.json -o 5910.csv
```

--------------------------------

### Mesh Preprocessing and Initial Analysis Example

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/pyvista_tests.ipynb

This snippet demonstrates a basic preprocessing step for two meshes, trimesh1 and trimesh2, by centering them around the origin. It then shows an example of clustering these meshes and visually plotting the results using pyvista. This is a typical workflow for preparing meshes for further analysis.

```python
import numpy as np
import pyvista as pv

# Assuming trimesh1, trimesh2, and cluster_meshes are defined elsewhere

t = np.mean(trimesh1.points, axis=0)
trichimesh1.points -= t
trichimesh2.points -= t

labels = cluster_meshes([trimesh1, trimesh2])
labels

# p = pv.Plotter()

# p.add_mesh(trimesh1, scalars=labels[0])
# p.add_mesh(trimesh2, scalars=labels[1])

# p.show()

intersect_surfaces([trimesh1, trimesh2])[0].plot()
```

--------------------------------

### Load and Sample Clustering Data using Pandas

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/Figures.ipynb

This snippet shows how to load clustering data from a CSV file using pandas. It includes commented-out code for sampling data and visualizing meshes based on cluster assignments, which would require additional setup and data loading functions.

```python
clustering = pd.read_csv("clustering_200k_30n_11f_average.csv")

# show_n = 5

# p = pv.Plotter(shape=(1, show_n))

# sample = clustering[clustering["cluster"] != 5].sample(n = show_n)

# selected_ids = sample["id"]

# i = 0
# for objid in selected_ids:
#     p.subplot(0, i)

#     mesh, geom, verts = load_building(objid, tile_csv=df)
#     p.add_title(str(sample.iloc[i]["cluster"]))
#     p.add_mesh(mesh)
    
#     i += 1

# p.show()
```

--------------------------------

### Filter and Save Example Buildings

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/Analysis.ipynb

Selects a predefined list of building IDs from the loaded DataFrame and saves them to a new CSV file named 'selected.csv'. This is used for in-depth analysis of specific example buildings.

```python
selected_ids = [
    "NL.IMBAG.Pand.0599100000702379-0",
    "NL.IMBAG.Pand.0518100001635181-0",
    "NL.IMBAG.Pand.0599100000701103-0",
#     "NL.IMBAG.Pand.0518100000225439-0",
    "NL.IMBAG.Pand.0518100000273015-0",
    "NL.IMBAG.Pand.0363100012075730-0",
    "NL.IMBAG.Pand.0363100012185598-0",
    "NL.IMBAG.Pand.0344100000031226-0",
    "NL.IMBAG.Pand.0344100000077683-0",
    "NL.IMBAG.Pand.0344100000099499-0",
    "NL.IMBAG.Pand.0599100000080428-0",
    "NL.IMBAG.Pand.0518100000230634-0",
#     "NL.IMBAG.Pand.0518100000206625-0",
    "NL.IMBAG.Pand.0518100000226316-0",
    "NL.IMBAG.Pand.0518100000282020-0",
    "NL.IMBAG.Pand.0599100000432858-0",
    "NL.IMBAG.Pand.0629100000020777-0",
    "NL.IMBAG.Pand.0363100012236081-0",
    "NL.IMBAG.Pand.0518100000222277-0"
]

df[df["id"].isin(selected_ids)].to_csv("selected.csv")
```

--------------------------------

### Batch Process CityJSON Files with val3dity Validation

Source: https://context7.com/tudelft3d/3d-building-metrics/llms.txt

This bash script demonstrates batch processing of multiple CityJSON files. It first uses `val3dity` to validate the geometry of each file and generate a validation report. Subsequently, it uses `cityStats.py` to compute metrics, incorporating the validation report. The script also shows examples of multi-threaded processing for large datasets and combining results from multiple CSV files into a single output.

```bash
# Batch process all CityJSON files in directory with validation
for i in *.json; do
    # Run val3dity validation
    val3dity $i --report "${i%.json}_v3.json"

    # Compute metrics with validation report
    python cityStats.py $i -o "${i%.json}.csv" -v "${i%.json}_v3.json"
done

# Multi-threaded batch processing for large datasets
for i in tile_*.json; do
    val3dity $i --report "${i%.json}_v3.json" &
done
wait

for i in tile_*.json; do
    python cityStats.py $i -o "${i%.json}.csv" \
        -v "${i%.json}_v3.json" \
        -j 8 \
        --density-2d 0.5 \
        --density-3d 0.25 &
done
wait

# Combine results into single file
head -n 1 tile_0001.csv > all_results.csv
tail -n +2 -q tile_*.csv >> all_results.csv
```

--------------------------------

### Process Multiple Files and Validate with Val3dity

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/README.md

Iterates through all .json files in the current directory, generates a val3dity report for each, and then runs cityStats to compute metrics and integrate validation results. This is useful for batch processing and quality assurance.

```bash
for i in *.json; do val3dity $i --report "${i%.json}_v3.json"; python cityStats.py $i -o "${i%.json}.csv" -v "${i%.json}_v3.json"; done
```

--------------------------------

### Visualize 3D Building Model with CityPlot

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/README.md

Launches a visualization tool to display a specific 3D building model. This is helpful for inspecting the geometry and understanding the spatial characteristics of the data.

```python
python cityPlot.py [file_path]
```

--------------------------------

### Measure Face Clustering Time

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/pyvista_tests.ipynb

This code measures the execution time of the `cluster_faces` function, providing a performance metric for the face clustering operation. It measures the time taken to cluster faces of the `main_mesh` using a threshold. The `time` module is used for timing, and the difference between start and end times is printed.

```python
import time

start = time.time()

labels, n_clusters = cluster_faces(face_planes(main_mesh))

end = time.time()
print(end - start)
```

--------------------------------

### Measure Mesh Clustering Time

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/pyvista_tests.ipynb

This code measures the execution time of the `cluster_meshes` function, which is assumed to perform some form of mesh clustering. It uses the `time` module to record the start and end times and prints the duration in seconds, providing a basic performance metric for the clustering operation. The `meshes` variable is assumed to be defined elsewhere and contains the meshes to cluster.

```python
import time

start = time.time()

cluster_meshes(meshes)

end = time.time()
print(end - start)
```

--------------------------------

### Calculate Distance Matrix for Clustering

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/pyvista_tests.ipynb

This code computes a distance matrix using scipy's `distance_matrix` to prepare data for clustering. It calculates two distance matrices, one for the original data and another for the negated data, and then takes the minimum of the two. This setup is presumably for handling orientations or to deal with potential ambiguities. The computed matrix is used in a clustering algorithm.

```python
import scipy
import numpy as np

data = np.array(face_planes(main_mesh))

dm1 = scipy.spatial.distance_matrix(data, data)
dm2 = scipy.spatial.distance_matrix(data, -data)

np.minimum(dm1, dm2)
```

--------------------------------

### Visualize CityJSON Buildings - Python

Source: https://context7.com/tudelft3d/3d-building-metrics/llms.txt

Displays 3D building models from CityJSON files with interactive visualization using PyVista. The script can visualize an entire city model, prompt for LoD selection if multiple exist, save visualizations to file formats like VTK, and allows programmatic visualization with custom styling, including coloring buildings by type.

```python
# Visualize entire city model
python cityPlot.py buildings.json

# If multiple LoDs exist, you'll be prompted to select one
# The script will display all buildings in an interactive 3D viewer

# Save visualization to file format
python cityPlot.py buildings.json --save
# Saves as cm.vtm (VTK multiblock format)

# Programmatic visualization with custom styling
import cityjson
import pyvista as pv
import json
import numpy as np

with open('buildings.json', 'r') as f:
    cm = json.load(f)

if "transform" in cm:
    s = cm["transform"]["scale"]
    t = cm["transform"]["translate"]
    verts = [[v[0]*s[0]+t[0], v[1]*s[1]+t[1], v[2]*s[2]+t[2]]
             for v in cm["vertices"]]
else:
    verts = cm["vertices"]

vertices = np.array(verts)

# Create plotter with custom settings
p = pv.Plotter()
p.background_color = "white"

# Add buildings with different colors by type
for obj_id in cm["CityObjects"]:
    obj = cm["CityObjects"][obj_id]
    if len(obj["geometry"]) > 0:
        mesh = cityjson.to_triangulated_polydata(obj["geometry"][0], vertices)

        # Color by building type
        if obj["type"] == "Building":
            p.add_mesh(mesh, color="lightblue", opacity=0.8)
        elif obj["type"] == "BuildingPart":
            p.add_mesh(mesh, color="orange", opacity=0.8)

p.show()
```

--------------------------------

### Compare Volumes of Fixed OBB and Original Mesh

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/pyvista_tests.ipynb

Compares the volumes of the original mesh and a fixed oriented bounding box mesh. This comparison helps in evaluating the accuracy of the bounding box representation.

```python
m = MeshFix(obb.clean().triangulate())
m.repair()
fixed_obb = m.mesh

print(f"Volume: {clean.volume}")
print(f"OBB: {obb.volume}")

# p = pv.Plotter()

# p.add_mesh(obb.clean(), show_edges=True, opacity=0.3)
# p.add_mesh(trimesh)

# p.show()
```

--------------------------------

### Load CityJSON File and Process Vertices - Python

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/pyvista_tests.ipynb

Sets up configurations for numerical output formatting and defines a dictionary of CityJSON file paths. It then loads a specified CityJSON file, processes its vertices by applying scale and translate transformations if present, and stores the vertices in a NumPy array.

```python
float_formatter = "{:.3f}".format
np.set_printoptions(formatter={'float_kind':float_formatter})

bag_tile = 7173

models = {
    "DenHaag": "",
    "Helsinki": rpath("~/Dropbox/CityJSON/Helsinki/CityGML_BUILDINGS_LOD2_NOTEXTURES_672496x2.json"),
    "Vienna": rpath("~/Dropbox/CityJSON/Vienna/Vienna_102081.json"),
    "Montreal": rpath("~/Dropbox/CityJSON/Montreal/VM05_2009.json"),
    "random": rpath("~/Downloads/random10_1.json"),
    "Delfshaven": rpath("~/Dropbox/CityJSON/rotterdam/Version Old/3-20-DELFSHAVEN_uncompressed.json"),
    "NYC": rpath("~/Dropbox/CityJSON/NewYork/NYCsubset.json"),
    "bag_tile": rpath(f"~/3DBAG_09/LoD2.2/{bag_tile}.json")
}

filename = models[bag_tile]

with open(filename) as file:
    cm = json.load(file)

if "transform" in cm:
    s = cm["transform"]["scale"]
    t = cm["transform"]["translate"]
    verts = [[v[0] * s[0] + t[0], v[1] * s[1] + t[1], v[2] * s[2] + t[2]]
            for v in cm["vertices"]]
else:
    verts = cm["vertices"]

# mesh points
vertices = np.array(verts)
```

--------------------------------

### Compute Building Statistics from CityJSON File (Python)

Source: https://context7.com/tudelft3d/3d-building-metrics/llms.txt

Processes a CityJSON file to compute geometric metrics for buildings, outputting results to a CSV file. Supports optional validation reports and multi-threaded execution. Key parameters include input file, output file, number of threads, and flags for validation or specific metric exclusions.

```bash
python cityStats.py building_data.json -o metrics.csv -v val3dity_report.json
python cityStats.py building_data.json -o metrics.csv -j 8 --density-2d 1.0 --density-3d 0.5
python cityStats.py building_data.json -p -f "building_id_123"
python cityStats.py building_data.json -o metrics.csv --without-indices
```

--------------------------------

### Compute Full 3D Oriented Bounding Box (OBB)

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/pyvista_tests.ipynb

Attempts to compute a fully 3D oriented bounding box for a set of points using the `pyobb` library. Note that this library might not be consistently reliable.

```python
from pyobb.obb import OBB

obb_full_3d = OBB.build_from_points(dataset.clean().points)
```

--------------------------------

### Visualize 3D Mesh with Extruded OBB

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/pyvista_tests.ipynb

Displays the original 3D mesh and its corresponding extruded 3D oriented bounding box using PyVista. This allows for visual assessment of how well the bounding box fits the object.

```python
p = pv.Plotter()

p.add_mesh(obb, opacity=0.3)
p.add_mesh(trimesh)

p.show()
```

--------------------------------

### Create Surface Grid for Mesh (Python)

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/pyvista_tests.ipynb

Generates a 2D grid of points that conforms to the surface of a 3D mesh. It projects mesh faces onto a 2D plane defined by the face's normal, creates a grid on this plane, and then transforms the grid points back into 3D space. Requires helper functions for geometry and axes calculation.

```python
from helpers.geometry import surface_normal
from shapely.geometry import Polygon

def to_3d(points, normal, origin):
    """Translate local 2D coordinates to 3D"""
    
    x_axis, y_axis = axes_of_normal(normal)
  
    return (np.repeat([origin], len(points), axis=0)
        + np.matmul(points, [x_axis, y_axis]))    

def axes_of_normal(normal):
    """Returns an x-axis and y-axis on a plane of the given normal"""
    if normal[2] > 0.001 or normal[2] < -0.001:
        x_axis = [1, 0, -normal[0]/normal[2]];
    elif normal[1] > 0.001 or normal[1] < -0.001:
        x_axis = [1, -normal[0]/normal[1], 0];
    else:
        x_axis = [-normal[1] / normal[0], 1, 0];
    
    x_axis = x_axis / np.linalg.norm(x_axis)
    y_axis = np.cross(normal, x_axis)

    return x_axis, y_axis

def project_2d(points, normal):
    origin = points[0]

    x_axis, y_axis = axes_of_normal(normal)
     
    return [[np.dot(p - origin, x_axis), np.dot(p - origin, y_axis)] for p in points]

def create_surface_grid(mesh, density=1):
    """Create a 2-dimensional grid along the surface of a 3D mesh"""
    
    result = []
    
    sized = mesh.compute_cell_sizes()
    
    for i in range(mesh.n_cells):
        if not mesh.cell_type(i) in [5, 6, 7, 9, 10]:
            continue
        
        pts = mesh.cell_points(i)
        
        normal = surface_normal(pts)
        
        pts_2d = project_2d(pts, normal)
        poly_2d = Polygon(pts_2d)
        
        grid = create_grid_2d(poly_2d, density)
        grid = MultiPoint(grid).intersection(poly_2d)
        
        if grid.is_empty:
            continue
        elif grid.geom_type == "Point":
            grid = np.array(grid.coords)
        else:
            grid = np.array([list(p.coords[0]) for p in grid.geoms])
        
        # TODO: Randomise the origin
        result.extend(list(to_3d(grid, normal, pts[0])))
            
    return result

s_grid = pv.PolyData(create_surface_grid(trimesh, 1))

p = pv.Plotter()

p.add_mesh(dataset, opacity=0.9)
# p.add_mesh(clean.extract_cells(82))
p.add_mesh(s_grid)

p.show()
```

--------------------------------

### Batch Plotting Clusters to Directory in Python

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/Figures.ipynb

This code iterates through sorted cluster labels and plots each one to a specified output directory. It includes logic to create the output directory if it doesn't exist and customizes the plotting for each cluster, such as setting cell size and disabling edge colors.

```python
import os

# Assuming 'rpath', 'plot_cluster', 'np', and 'pd' are defined elsewhere
# Example placeholder definitions for context:
def rpath(path_str):
    return path_str.replace("~", "/home/user")

# Mocking plot_cluster for standalone execution example
def plot_cluster(cluster_label, cell_size=None, filename=None, show_class=True, show_ids=True, edge_color=None):
    print(f"Plotting cluster {cluster_label} to {filename} with cell_size={cell_size}, edge_color={edge_color}")

# Placeholder data
clustering = pd.DataFrame({'cluster': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 'id': range(10)})

out_folder = rpath("~/figures/clusters/sample_200k_30n_average")

if not os.path.exists(out_folder):
    os.makedirs(out_folder)

labels = np.sort(pd.unique(clustering["cluster"]))

for i, c in enumerate(labels[6:]):
    print(f"Plotting {c}...")
    plot_cluster(c, cell_size=[683, 683], filename=f"{out_folder}/{c}.pdf", show_class=False, show_ids=False, edge_color=None)

print("Done!")

```

--------------------------------

### Compute Grid-Based Shape Indices (Python)

Source: https://context7.com/tudelft3d/3d-building-metrics/llms.txt

Calculates various shape indices (proximity, spin, depth, dispersion, girth, roughness) for 2D polygons and 3D meshes using interior grid points. Requires PyVista and shape_index libraries. Handles both 2D and 3D cases with adjustable grid density.

```python
import pyvista as pv
import shape_index as si
from shapely.geometry import Polygon

# 2D grid-based indices with custom density
footprint = Polygon([(0,0), (10,0), (12,5), (10,10), (0,10)])
density_2d = 0.5  # grid spacing

# Proximity index: how close points are to centroid
proximity_2d = si.proximity_2d(footprint, density=density_2d)
# Returns: normalized value (2/3 * √(area/π) / mean_distance_to_centroid)

# Spin index: rotational distribution measure
spin_2d = si.spin_2d(footprint, density=density_2d)
# Returns: normalized value (0.5 * (area/π) / mean_squared_distance)

# Depth index: how far interior points are from boundary
depth_2d = si.depth_2d(footprint, density=density_2d)
# Returns: normalized depth (3 * mean_distance_to_boundary / √(area/π))

# Dispersion index: how uniformly shape expands from center
dispersion_2d = si.dispersion_2d(footprint, density=0.2)
# Returns: 0-1 (1 = uniform expansion, like a circle)

# 3D grid-based indices
mesh = pv.Cube()
density_3d = 0.25

# Create 3D interior grid
grid = si.create_grid_3d(mesh, density=density_3d, check_surface=False)
# Returns: array of 3D points inside the volume

# Proximity 3D: distance to center of mass
proximity_3d = si.proximity_3d(mesh, grid=grid, density=density_3d)
# Returns: normalized value

# Girth index: largest inscribed sphere ratio
girth_3d = si.girth_3d(mesh, grid=grid, density=density_3d)
# Returns: radius_inscribed / radius_equal_volume_sphere

# Roughness index: surface irregularity measure
roughness_3d = si.roughness_index_3d(mesh, grid, density=0.5)
# Returns: normalized roughness based on distance variance
```

--------------------------------

### Detect and Visualize Open Edges in a Mesh

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/pyvista_tests.ipynb

Identifies and visualizes open edges in a 3D mesh using Trimesh and PyVista to check for watertightness. Open edges are displayed as red lines.

```python
edges = trimesh.extract_feature_edges(boundary_edges=True,
                           feature_edges=False,
                           manifold_edges=False)

p = pv.Plotter()

p.add_mesh(trimesh, opacity=1.0, show_edges=True)
if trimesh.n_open_edges:
    p.add_mesh(edges, color='red', line_width=10)

p.add_title(f"{obj} {'is watertight' if trimesh.n_open_edges == 0 else f'has {trimesh.n_open_edges} open edges'}", 8)    

p.show()
```

--------------------------------

### Run CityStats for Geometric Metrics

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/README.md

Computes and saves geometric metrics for a given 3D building file. Supports single-threaded execution by default, with options for multithreading and val3dity report integration. Output is saved to a CSV file.

```python
python cityStats.py [file_path] -o [output.csv] [-v val3dity_report.json]
```

```python
python cityStats.py [file_path] -j [number]
```

```python
python cityStats.py [file_path] -p -f [unique_id]
```

--------------------------------

### Visualize 2D PCA Variance Ratio with Plotnine

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/Analysis.ipynb

Utilizes the plotnine library to create a horizontal bar chart showing the variance explained by each principal component from the 2D PCA. This visualization helps in interpreting the contribution of each component.

```python
variance_ratio_2d = pd.DataFrame(np.transpose(np.array([indices_2d, pca_2d.explained_variance_ratio_])), columns=["metric", "value"]) #.set_index("metric")
variance_ratio_2d["value"] = variance_ratio_2d["value"].astype("float")

from plotnine import *
%matplotlib inline

(ggplot(variance_ratio_2d)         # defining what data to use
 + aes(x="reorder(metric, value)", y="value")    # defining what variable to use
 + geom_bar(stat="identity")  # defining the type of plot to use
 + coord_flip()
)
```

--------------------------------

### Compare Buildings Across Different Levels of Detail - Python

Source: https://context7.com/tudelft3d/3d-building-metrics/llms.txt

Compares the same building represented at different Levels of Detail (LoDs) by computing volumetric differences, intersections, and symmetric differences. This script facilitates the assessment of geometric changes or data quality between different LoD representations. The output CSV includes source volume, destination volume, intersection volume, symmetric difference volume, and destination minus source volume.

```python
# Compare buildings between two CityJSON files at different LoDs
python cityCompare.py source_lod12.json destination_lod22.json \
    --lod-source "1.2" \
    --lod-destination "2.2" \
    -o comparison_results.csv \
    -j 4 \
    --repair

# Compare with visualization for specific building
python cityCompare.py source.json destination.json \
    --lod-source "1.2" \
    --lod-destination "2.2" \
    -f "building_id_456" \
    --plot

# Export geometry differences for further analysis
python cityCompare.py source.json destination.json \
    --lod-source "1.2" \
    --lod-destination "2.2" \
    -e \
    -o results.csv

# Output CSV contains:
# - source_volume: volume of source LoD geometry
# - destination_volume: volume of destination LoD geometry
# - intersection_volume: overlapping volume
# - symmetric_difference_volume: total non-overlapping volume
# - destination_minus_source: volume added in destination
```

--------------------------------

### Visualize Voxel Centers and Mesh with PyVista

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/pyvista_tests.ipynb

Renders a voxelized mesh, its cell centers, the original mesh, and the average of voxel centers using PyVista. This visualization helps in understanding the spatial distribution of voxels and their relation to the original geometry.

```python
p = pv.Plotter()

p.add_mesh(voxel, opacity=0.2, show_edges=True, color='yellow')
p.add_mesh(voxel.cell_centers(), color='black')
p.add_mesh(clean, color='grey')
p.add_mesh(pv.PolyData(np.mean(voxel.cell_centers().points, axis=0)), color='white')

p.show()
```

--------------------------------

### Compute Height Statistics using cityStats

Source: https://context7.com/tudelft3d/3d-building-metrics/llms.txt

This Python script demonstrates how to compute various statistical measures for building heights using the `compute_stats` function from the `cityStats` module. It calculates mean, median, mode, range, standard deviation, percentiles, and percentage of range for a given list of height values. Dependencies include `cityStats` and `numpy`.

```python
from cityStats import compute_stats
import numpy as np

roof_heights = [12.5, 12.8, 13.0, 13.2, 13.5, 13.5, 13.5, 13.8, 14.0, 14.2]

stats = compute_stats(roof_heights, percentile=90, percentage=75)

print(f"Mean height: {stats['Mean']:.2f} m")
print(f"Median height: {stats['Median']:.2f} m")
print(f"Max height: {stats['Max']:.2f} m")
print(f"Min height: {stats['Min']:.2f} m")
print(f"Height range: {stats['Range']:.2f} m")
print(f"Standard deviation: {stats['Std']:.2f} m")
print(f"90th percentile: {stats['Percentile']:.2f} m")
print(f"75% of range: {stats['Percentage']:.2f} m")

if stats['ModeStatus'] == 'Y':
    print(f"Mode (most common): {stats['Mode']:.2f} m")
else:
    print("No mode found (all values unique)")
```

--------------------------------

### Visualize 3D PCA Variance Ratio with Plotnine

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/Analysis.ipynb

Uses the plotnine library to create a bar chart visualizing the variance explained by each principal component from the 3D PCA. The chart is horizontally oriented for better readability of component labels.

```python
variance_ratio_3d = pd.DataFrame(np.transpose(np.array([["PC%d" % i for i in range(pca_3d.n_components)], pca_3d.explained_variance_ratio_])), columns=["PC", "value"]) #.set_index("metric")
variance_ratio_3d["value"] = variance_ratio_3d["value"].astype("float")

from plotnine import *
%matplotlib inline

(ggplot(variance_ratio_3d)         # defining what data to use
 + aes(x="reorder(PC, value)", y="value")    # defining what variable to use
 + geom_bar(stat="identity")  # defining the type of plot to use
 + coord_flip()
#  + ggtitle("")
)
```

--------------------------------

### Compute Surface Areas by Semantic Type (Python)

Source: https://context7.com/tudelft3d/3d-building-metrics/llms.txt

Calculates and aggregates areas for different semantic surface types (GroundSurface, WallSurface, RoofSurface) from a PyVista mesh derived from CityJSON. Requires cityjson, geometry, and numpy libraries. Also provides counts of points and surfaces for each type and extracts specific surface points.

```python
import cityjson
import geometry
import numpy as np

# Assume cm, vertices loaded from CityJSON
building = cm["CityObjects"]["building_id_123"]
geom = building["geometry"][0]

# Create mesh with semantic surface labels
mesh = cityjson.to_polydata(geom, vertices)

# Compute areas by surface type
area, point_count, surface_count = geometry.area_by_surface(mesh)

print(f"Ground surface area: {area['GroundSurface']:.2f} m²")
print(f"Wall surface area: {area['WallSurface']:.2f} m²")
print(f"Roof surface area: {area['RoofSurface']:.2f} m²")

print(f"Ground points: {point_count['GroundSurface']}")
print(f"Wall points: {point_count['WallSurface']}")
print(f"Roof points: {surface_count['RoofSurface']}")

print(f"Ground surfaces: {surface_count['GroundSurface']}")
print(f"Wall surfaces: {surface_count['WallSurface']}")
print(f"Roof surfaces: {surface_count['RoofSurface']}")

# Get points of specific surface type
roof_points = geometry.get_points_of_type(mesh, "RoofSurface")
if len(roof_points) > 0:
    max_height = np.max(roof_points[:, 2])
    min_height = np.min(roof_points[:, 2])
    print(f"Roof height range: {min_height:.2f} - {max_height:.2f} m")
```

--------------------------------

### Calculate 3D Shape Indices (Python)

Source: https://context7.com/tudelft3d/3d-building-metrics/llms.txt

Computes volumetric shape metrics for 3D building meshes using PyVista and custom functions. Includes hemisphericality, cubeness, 3D fractality, 3D convexity, and circumference index, quantifying the volume's shape characteristics relative to its surface area and convex hull.

```python
import pyvista as pv
import shape_index as si
import numpy as np
import scipy.spatial as ss

# Create or load a 3D building mesh
points = np.array([[0,0,0], [1,0,0], [1,1,0], [0,1,0],
                   [0,0,1], [1,0,1], [1,1,1], [0,1,1]])
faces = np.array([4,0,1,2,3, 4,4,5,6,7, 4,0,1,5,4,
                  4,2,3,7,6, 4,0,3,7,4, 4,1,2,6,5])
mesh = pv.PolyData(points, faces)

# Hemisphericality: 3D equivalent of circularity
hemisphericality = si.hemisphericality(mesh)
# Returns: ~0.9 (3√2 * √π * volume / area^(3/2))

# Cubeness: measures how cube-like the volume is
cubeness = si.cubeness(mesh)
# Returns: ~1.0 (6 * volume^(2/3) / area)

# Fractality 3D: surface complexity measure
fractality_3d = si.fractality_3d(mesh)
# Returns: complexity value (1 - log(volume) / (3/2 * log(area)))

# Convexity 3D: volume ratio to convex hull
points_array = mesh.clean().points
ch_volume = ss.ConvexHull(points_array).volume
convexity_3d = mesh.volume / ch_volume
# Returns: 1.0 (perfect convex volume)

# Circumference index: normalized surface area
circumference = si.circumference_index_3d(mesh)
# Returns: ~1.0 (4π * (3V/4π)^(2/3) / area)
```

--------------------------------

### Voxelizing and Plotting a Mesh

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/Figures.ipynb

This snippet demonstrates how to voxelize a given mesh and then visualize the resulting voxel grid. It utilizes PyVista's `voxelize` function to create a voxel representation of the mesh and then plots this voxelized mesh with customizable appearance, including line width, color, and edge visibility. The output can be saved to a file.

```python
import pyvista as pv
from os.path import expanduser as rpath

# Assuming 'mesh_dest' is a pre-defined PyVista mesh object
# voxel = pv.voxelize(mesh_dest, density=0.5, check_surface=False)

p = pv.Plotter(window_size=[2048, 2048])

p.background_color = "white"

# Add the voxelized mesh to the plotter
# p.add_mesh(voxel, line_width=2.0, color='lightgrey', ambient=0.5, diffuse=0.8, show_edges=True)

p.show()

# Save the visualization to a PDF file
# p.save_graphic(rpath("~/figures/example-grid.pdf"))

```

--------------------------------

### Generate 2D Grid for Polygon

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/pyvista_tests.ipynb

Creates a grid of points within a given 2D polygon. It takes a shapely Polygon object and a density value as input, returning a list of (x, y) coordinates.

```python
from shapely.geometry import Point, MultiPoint, Polygon
import numpy as np

def create_grid_2d(shape, density):
    """Return the grid for a given polygon"""
    
    x_min, y_min, x_max, y_max = shape.bounds
    x = np.arange(x_min, x_max, density)
    y = np.arange(y_min, y_max, density)
    x, y = np.meshgrid(x, y)
    
    x = np.hstack(x)
    y = np.hstack(y)
    
    return [(x[i], y[i]) for i in range(len(x))]
```

--------------------------------

### Calculate 2D Shape Indices (Python)

Source: https://context7.com/tudelft3d/3d-building-metrics/llms.txt

Computes various 2D shape complexity and form metrics for building footprints represented as Shapely Polygons. Metrics include circularity, convexity, fractality, squareness, rectangularity, and elongation, characterizing the shape's properties based on its area and perimeter.

```python
import shape_index as si
from shapely.geometry import Polygon

# Create a building footprint polygon
footprint = Polygon([(0, 0), (10, 0), (10, 8), (0, 8)])

# Circularity: measures how close the shape is to a circle (0-1, 1 = perfect circle)
cirularity = si.circularity(footprint)
# Returns: 0.785 (4π * area / perimeter²)

# Convexity: ratio of area to convex hull area
convexity_2d = footprint.area / footprint.convex_hull.area
# Returns: 1.0 (perfect convex shape)

# Fractality: measures boundary complexity
fractality = si.fractality_2d(footprint)
# Returns: ~0.0 (low complexity for simple rectangle)

# Squareness: measures how square-like the shape is
squareness = si.squareness(footprint)
# Returns: 0.894 (4√area / perimeter)

# Rectangularity: ratio of area to minimum rotated rectangle
rectangularity = footprint.area / footprint.minimum_rotated_rectangle.area
# Returns: 1.0 (perfect rectangle)

# Elongation: measure of length-to-width ratio
S, L = si.get_box_dimensions(footprint.minimum_rotated_rectangle)
elongation = si.elongation(S, L)
# Returns: 0.2 (1 - min_dim/max_dim)
```

--------------------------------

### Parse and Convert CityJSON Geometries (Python)

Source: https://context7.com/tudelft3d/3d-building-metrics/llms.txt

Loads CityJSON files and converts geometries into PyVista meshes or Shapely polygons. Handles coordinate transformations and extracts semantic information. Requires json, numpy, cityjson, and pyvista libraries. Outputs include triangulated meshes, 2D footprints, and raw point data.

```python
import json
import numpy as np
import cityjson
import pyvista as pv

# Load CityJSON file
with open('buildings.json', 'r') as f:
    cm = json.load(f)

# Handle coordinate transformations if present
if "transform" in cm:
    s = cm["transform"]["scale"]
    t = cm["transform"]["translate"]
    verts = [[v[0]*s[0]+t[0], v[1]*s[1]+t[1], v[2]*s[2]+t[2]]
             for v in cm["vertices"]]
else:
    verts = cm["vertices"]

vertices = np.array(verts)

# Get a specific building
building = cm["CityObjects"]["building_id_001"]
geom = building["geometry"][0]

# Convert to PyVista mesh (preserves semantic surfaces)
mesh = cityjson.to_polydata(geom, vertices)
print(f"Mesh has {mesh.n_points} points, {mesh.n_cells} cells")
if "semantics" in mesh.cell_data:
    print(f"Semantics: {mesh.cell_data['semantics']}")

# Convert to triangulated mesh for volume calculations
tri_mesh = cityjson.to_triangulated_polydata(geom, vertices)
print(f"Volume: {tri_mesh.volume:.2f} cubic meters")

# Extract 2D footprint as Shapely polygon (ground surfaces only)
footprint = cityjson.to_shapely(geom, vertices, ground_only=True)
print(f"Footprint area: {footprint.area:.2f} square meters")

# Get all building points
points = cityjson.get_points(geom, verts)
print(f"Total points in geometry: {len(points)}")

# Get bounding boxxmin, xmax, ymin, ymax, zmin, zmax = cityjson.get_bbox(geom, verts)
print(f"BBox: X[{xmin:.1f}, {xmax:.1f}] Y[{ymin:.1f}, {ymax:.1f}] Z[{zmin:.1f}, {zmax:.1f}]")
```

--------------------------------

### Calculate Plane Parameters for Mesh Faces

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/pyvista_tests.ipynb

These functions calculate the parameters (a, b, c, d) of the plane equation for each face of a given mesh. It utilizes numpy for calculations and assumes mesh face normals and cell points are available. The output is a list of plane parameter arrays.

```python
import numpy as np

def plane_params(normal, origin, rounding=2, absolute=True):
    """Returns the params (a, b, c, d) of the plane equation"""
    a, b, c = np.round_(normal, 3)
    x0, y0, z0 = origin
    
    d = -(a * x0 + b * y0 + c * z0)
    
    if rounding >= 0:
        d = round(d, rounding)
    
    return np.array([a, b, c, d])

def face_planes(mesh):
    return [plane_params(mesh.face_normals[i], mesh.cell_points(i)[0]) for i in range(mesh.n_cells)]
```

--------------------------------

### Compare Voxelized and Actual Mesh Volumes

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/pyvista_tests.ipynb

Calculates and prints the volume of a voxelized mesh and the original mesh. This comparison serves as an indicator of the mesh's validity and the accuracy of the voxelization process.

```python
print(f"Voxel: {voxel.volume}")
print(f"Actual: {clean.volume}")
```

--------------------------------

### Load and Prepare Building Data

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/Figures.ipynb

Loads building data from CityJSON files, processes transformations, and converts CityJSON geometry to a triangulated mesh using PyVista. It handles potential scale and translation transformations within the CityJSON structure. Dependencies include pyvista, pymesh, cityjson, and pandas.

```python
import json
import pyvista as pv
import pymesh
from helpers.mesh import *
import cityjson
import pandas as pd
import os

def rpath(path):
    return os.path.expanduser(path)

df = pd.read_csv(rpath("~/3DBAG_09/all/lod2.2.csv"))

def load_citymodel(file):
    cm = json.load(file)

    if "transform" in cm:
        s = cm["transform"]["scale"]
        t = cm["transform"]["translate"]
        verts = [[v[0] * s[0] + t[0], v[1] * s[1] + t[1], v[2] * s[2] + t[2]]
                for v in cm["vertices"]]
    else:
        verts = cm["vertices"]

    vertices = np.array(verts)

    return cm, vertices

def get_geometry(co, lod):
    """Returns the geometry of the given LoD. If lod is None then it returns the first one."""

    if len(co["geometry"]) == 0:
        return None

    if lod is None:
        return co["geometry"][0]

    for geom in co["geometry"]:
        if str(geom["lod"]) == str(lod):
            return geom

def is_valid(mesh):
    return mesh.volume > 0 and mesh.n_open_edges == 0

def load_building(objid, tile_id=None, tile_csv=None, lod="2.2"):
    
    if tile_id is None:
        tile_id = tile_csv.set_index("id").loc[objid]["tile_id"]

    filename = rpath(f"~/3DBAG_09/{tile_id}.json")

    with open(filename, "rb") as file:
        cm, verts = load_citymodel(file)

    building = cm["CityObjects"][objid]

    geom = get_geometry(building, lod)
    mesh = cityjson.to_triangulated_polydata(geom, verts)
    
    return mesh, geom, verts

```

--------------------------------

### Import Libraries and Define Path Function - Python

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/pyvista_tests.ipynb

Imports necessary libraries for data manipulation, CityJSON processing, and 3D visualization. It also defines a helper function `rpath` to resolve user-specific paths, which is useful for accessing local data files.

```python
import json
import numpy as np
import pyvista as pv
from pymeshfix import MeshFix
import cityjson
import shapely
import math
from tqdm import tqdm
import os

def rpath(path):
    return os.path.expanduser(path)
```

--------------------------------

### Extrude 2D OBB to 3D Bounding Box using PyVista

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/pyvista_tests.ipynb

Extrudes a 2D oriented bounding box into a 3D bounding box using PyVista. The extrusion is performed along the Z-axis based on the calculated height of the object. Points are temporarily moved to the origin to mitigate potential precision issues in VTK.

```python
ground_z = np.min(dataset.clean().points[:, 2])
height = np.max(dataset.clean().points[:, 2]) - ground_z
box = np.array([[p[0], p[1], ground_z] for p in list(obb_2d.corner_points)])

obb = pv.PolyData(box).delaunay_2d()
pts = obb.points

t = np.mean(pts, axis=0)

# We need to move the points to the origin before extruding due to VTK's precision issues
obb.points = obb.points - t
obb = obb.extrude([0.0, 0.0, height])
obb.points = obb.points + t
```

--------------------------------

### Creating and Plotting Surface Grid on Mesh

Source: https://github.com/tudelft3d/3d-building-metrics/blob/main/Figures.ipynb

This code generates a surface grid for a given mesh and overlays it for visualization. It uses a hypothetical `create_surface_grid` function (presumably from `shape_index` module) to compute the grid points. The resulting grid points are then plotted as spheres on top of the mesh and its edges using PyVista. The visualization can be saved to a file.

```python
import pyvista as pv
from os.path import expanduser as rpath
# Assuming 'mesh_dest' is a pre-defined PyVista mesh object
# from shape_index import create_surface_grid

# sgrid = pv.PolyData(create_surface_grid(mesh_dest, 0.6))

p = pv.Plotter(window_size=[2048, 2048])

# Add the original mesh with specific color and ambient/diffuse properties
# p.add_mesh(mesh_dest, color="lightgrey", ambient=1, diffuse=0)
# Add mesh edges with specified color and line width
# p.add_mesh(mesh_dest.extract_feature_edges(), color="black", line_width=10)
# Add the surface grid points rendered as spheres
# p.add_mesh(sgrid, render_points_as_spheres=True, point_size=10, color="black")

p.show()

# Save the visualization to a PDF file
# p.save_graphic(rpath("~/figures/example-surface-grid.pdf"))

```