Try mybimapi/triangle

mybimapi/triangle_heatmap

IFC Geometry Density Heatmap

tag v2

person mybimapi

This function generates a high-resolution (2048x2048), top-down 2D heatmap image that visualizes the geometric complexity and density of a given 3D IFC building model.

The script extracts all 3D mesh triangles from the model, projects them onto a flat 2D plane (top view), and rasterizes them into a grid. It calculates a "density weight" for every triangle based on its size—smaller, highly detailed triangles contribute more weight than large, simple ones. The result is mapped to an "inferno" color scale (ranging from dark purple for low detail to bright yellow for high detail) with a transparent background.

This tool is exceptionally useful for quickly identifying geometric bottlenecks, heavily tessellated objects (like complex furniture or detailed curved elements), and overall model complexity without needing to manually inspect the 3D file.

Parameters

model (ifcopenshell.file) The loaded IFC file object you wish to analyze.

Returns

PIL.Image.Image Returns a Python Imaging Library (PIL) Image object in RGBA format.
Size: 2048 x 2048 pixels.
Visuals: Uses the matplotlib inferno colormap. Empty space where no geometry exists is completely transparent, allowing the image to be easily overlaid onto floor plans or web canvases.

How It Works Under the Hood

Extraction: Iterates through the IFC model using multiprocessing (4 threads) to extract all geometric shapes as triangulated meshes.
2D Projection: Discards the Z-axis (height) and normalizes the X and Y coordinates to fit perfectly within the 2048x2048 pixel boundaries, preserving the model's original aspect ratio.
Density Rasterization: Uses a high-performance, JIT-compiled (Numba) algorithm to draw the triangles onto a grid. Small, clustered triangles accumulate higher density values.
Normalization & Rendering: Applies a logarithmic scale to the density map to balance out extreme hotspots, then converts the data into an RGBA heatmap image.

code Source code

import ifcopenshell
import ifcopenshell.geom
import numpy as np
from matplotlib import colormaps
from numba import njit
from PIL import Image

SIZE = 2048

def run(model: ifcopenshell.file) -> Image:
    triangles = read_triangles(model)
    triangles = fit_triangles(triangles)

    density = np.zeros((SIZE, SIZE), dtype=np.float32)
    rasterize_triangles(triangles, density)

    norm = normalize_density(density)
    rgba = heatmap_rgba(norm, density)

    return Image.fromarray(rgba, mode="RGBA")


def shape_triangles(shape):
    geometry = shape.geometry

    vertices = np.asarray(geometry.verts, dtype=np.float32).reshape(-1, 3)
    faces = np.asarray(geometry.faces, dtype=np.int32).reshape(-1, 3)

    return vertices[faces][:, :, :2]


def read_triangles(model):
    settings = ifcopenshell.geom.settings()
    settings.set("use-world-coords", True)

    chunks = []

    for shape in ifcopenshell.geom.iterate(
        settings,
        model,
        num_threads=4,
        geometry_library="hybrid-cgal-simple-opencascade",
    ):
        triangles = shape_triangles(shape)

        if len(triangles):
            chunks.append(triangles)

    if not chunks:
        return np.empty((0, 3, 2), dtype=np.float32)

    return np.concatenate(chunks, axis=0).astype(np.float32)


def fit_triangles(triangles):
    if len(triangles) == 0:
        return triangles

    xy = triangles.reshape(-1, 2)

    lo = xy.min(axis=0)
    hi = xy.max(axis=0)

    extent = hi - lo
    span = float(max(extent[0], extent[1]))

    if span == 0.0:
        return triangles * 0.0

    scale = (SIZE - 1) / span

    fitted = triangles.copy()

    fitted[:, :, 0] = (fitted[:, :, 0] - lo[0]) * scale
    fitted[:, :, 1] = (hi[1] - fitted[:, :, 1]) * scale

    used = extent * scale
    offset = (SIZE - 1 - used) * 0.5

    fitted[:, :, 0] += offset[0]
    fitted[:, :, 1] += offset[1]

    return fitted.astype(np.float32)


@njit(cache=True)
def edge(ax, ay, bx, by, px, py):
    return (px - ax) * (by - ay) - (py - ay) * (bx - ax)


@njit(cache=True)
def length(x0, y0, x1, y1):
    dx = x1 - x0
    dy = y1 - y0
    return np.sqrt(dx * dx + dy * dy)


@njit(cache=True)
def rasterize_triangles(triangles, density):
    h, w = density.shape

    grow_pixels = 1.0
    area_epsilon = 1.0

    for i in range(triangles.shape[0]):
        x0 = triangles[i, 0, 0]
        y0 = triangles[i, 0, 1]
        x1 = triangles[i, 1, 0]
        y1 = triangles[i, 1, 1]
        x2 = triangles[i, 2, 0]
        y2 = triangles[i, 2, 1]

        signed_area2 = edge(x0, y0, x1, y1, x2, y2)
        projected_area = abs(signed_area2) * 0.5
        weight = 1.0 / np.sqrt(projected_area + area_epsilon)

        l0 = length(x1, y1, x2, y2)
        l1 = length(x2, y2, x0, y0)
        l2 = length(x0, y0, x1, y1)

        t0 = grow_pixels * l0
        t1 = grow_pixels * l1
        t2 = grow_pixels * l2

        xmin = max(0, int(np.floor(min(x0, x1, x2) - grow_pixels)))
        xmax = min(w - 1, int(np.ceil(max(x0, x1, x2) + grow_pixels)))
        ymin = max(0, int(np.floor(min(y0, y1, y2) - grow_pixels)))
        ymax = min(h - 1, int(np.ceil(max(y0, y1, y2) + grow_pixels)))

        if signed_area2 >= 0.0:
            for y in range(ymin, ymax + 1):
                py = y + 0.5

                for x in range(xmin, xmax + 1):
                    px = x + 0.5

                    e0 = edge(x1, y1, x2, y2, px, py)
                    e1 = edge(x2, y2, x0, y0, px, py)
                    e2 = edge(x0, y0, x1, y1, px, py)

                    if e0 >= -t0 and e1 >= -t1 and e2 >= -t2:
                        density[y, x] += weight
        else:
            for y in range(ymin, ymax + 1):
                py = y + 0.5

                for x in range(xmin, xmax + 1):
                    px = x + 0.5

                    e0 = edge(x1, y1, x2, y2, px, py)
                    e1 = edge(x2, y2, x0, y0, px, py)
                    e2 = edge(x0, y0, x1, y1, px, py)

                    if e0 <= t0 and e1 <= t1 and e2 <= t2:
                        density[y, x] += weight


def normalize_density(density):
    norm = np.log1p(density)
    vmax = norm.max()

    if vmax == 0.0:
        return norm

    return norm / vmax


def heatmap_rgba(norm, density):
    cmap = colormaps["inferno"]

    rgba = cmap(norm, bytes=True)
    rgba[:, :, 3] = np.where(density > 0.0, 255, 0).astype(np.uint8)

    return rgba

Output