differential_subdivision_script.py

# /// script
# dependencies = [
#     "shapely==2.2.0",
#     "numpy==2.5.3",
#     "scipy"
# ]
# ///

from math import atan2

import numpy as np
from scipy.optimize import differential_evolution
from System import Guid
from System.Collections.Generic import List
from shapely import affinity
from shapely.geometry import Polygon, box
from shapely.ops import unary_union
from Autodesk.Revit import DB, UI
from Autodesk.Revit.UI.Selection import ISelectionFilter, ObjectType
from System.Diagnostics import Trace


class RoomSpaceSelectionFilter(ISelectionFilter):
    __namespace__ = str(Guid.NewGuid())  # pythonnet3: unique namespace for interface impl

    def AllowElement(self, element):
        return isinstance(element, (DB.Architecture.Room, DB.Mechanical.Space))

    def AllowReference(self, reference, position):
        return False


def _curve_to_xy_points(curve):
    return [(p.X, p.Y) for p in curve.Tessellate()]


def _ring_from_boundary_loop(loop):
    pts = []
    for seg in loop:
        seg_pts = _curve_to_xy_points(seg.GetCurve())
        if not seg_pts:
            continue
        if pts:
            pts.extend(seg_pts[1:])
        else:
            pts.extend(seg_pts)

    if len(pts) < 3:
        return None
    if pts[0] != pts[-1]:
        pts.append(pts[0])
    return pts


def _spatial_element_to_polygon(spatial_element):
    options = DB.SpatialElementBoundaryOptions()
    loops = spatial_element.GetBoundarySegments(options)
    if loops is None:
        return None

    ring_data = []
    for loop in loops:
        ring = _ring_from_boundary_loop(loop)
        if ring is None:
            continue

        poly = Polygon(ring).buffer(0)
        if poly.is_empty or poly.area < 1e-8:
            continue
        ring_data.append((ring, poly.area))

    if not ring_data:
        return None

    ring_data.sort(key=lambda x: x[1], reverse=True)
    shell = ring_data[0][0]
    holes = [ring for ring, _ in ring_data[1:]]
    return Polygon(shell, holes).buffer(0)


def _get_alignment_angle(poly):
    mrr = poly.minimum_rotated_rectangle
    if mrr.is_empty:
        return 0.0

    coords = list(mrr.exterior.coords)
    if len(coords) < 3:
        return 0.0

    edges = []
    for i in range(4):
        x1, y1 = coords[i]
        x2, y2 = coords[i + 1]
        dx = x2 - x1
        dy = y2 - y1
        edges.append((dx * dx + dy * dy, dx, dy))

    _, dx, dy = max(edges, key=lambda e: e[0])
    return np.degrees(atan2(dy, dx))


def _polygon_to_lines(poly, z):
    lines = List[DB.Line]()

    coords = list(poly.exterior.coords)
    for i in range(len(coords) - 1):
        x1, y1 = coords[i]
        x2, y2 = coords[i + 1]
        p1 = DB.XYZ(x1, y1, z)
        p2 = DB.XYZ(x2, y2, z)
        lines.Add(DB.Line.CreateBound(p1, p2))

    return lines




def get_largest_internal_rect(poly, maxiter=120, popsize=15, seed=1):
    minx, miny, maxx, maxy = poly.bounds
    sx, sy = maxx - minx, maxy - miny
    pen = max(poly.area, sx * sy, 1.0) * 1e6

    def make_rect(v):
        cx, cy, hx, hy = v
        return box(cx - hx, cy - hy, cx + hx, cy + hy)

    def obj(v):
        r = make_rect(v)
        return -r.area + pen * r.difference(poly).area

    res = differential_evolution(
        obj,
        bounds=[(minx, maxx), (miny, maxy), (0, sx / 2.0), (0, sy / 2.0)],
        maxiter=maxiter,
        popsize=popsize,
        tol=1e-7,
        polish=True,
        seed=seed,
        workers=1
    )

    cx, cy, hx, hy = res.x
    for _ in range(80):
        r = box(cx - hx, cy - hy, cx + hx, cy + hy)
        if r.area > 0.01 and poly.covers(r):
            return r
        hx *= 0.995
        hy *= 0.995

    return None


def get_polys(geo):
    if geo.is_empty:
        return []
    if geo.geom_type == 'Polygon':
        return [geo]
    return [g for g in getattr(geo, 'geoms', []) if g.geom_type == 'Polygon']


def subtractive_decomposition(poly, max_rooms=15):
    remaining, rectangles = poly, []

    for _ in range(max_rooms):
        if remaining.is_empty or remaining.area < .01:
            break

        pieces = sorted(get_polys(remaining), key=lambda p: p.area, reverse=True)
        if not pieces:
            break

        target, others = pieces[0], pieces[1:]
        rect = get_largest_internal_rect(target)
        ok = rect is not None and rect.area >= .01

        rectangles.append(rect if ok else target)
        rest = [target.difference(rect)] if ok else []
        remaining = unary_union(rest + others) if rest or others else Polygon()

    return rectangles


def main():
    uiapp: UI.UIApplication = __revit__  # type: ignore  # noqa: F821
    uidoc = uiapp.ActiveUIDocument
    doc = uidoc.Document

    try:
        ref = uidoc.Selection.PickObject(
            ObjectType.Element,
            RoomSpaceSelectionFilter(),
            "Select a Room or Space",
        )
    except Exception as exc:
        print("Selection cancelled or failed: {}".format(exc))
        return

    spatial = doc.GetElement(ref)
    poly = _spatial_element_to_polygon(spatial)
    if poly is None or poly.is_empty or poly.area < 0.01:
        print("Selected element has no valid boundary polygon.")
        return

    angle = _get_alignment_angle(poly)
    aligned = affinity.rotate(poly, -angle, origin=(0, 0))

    final_decomposition = subtractive_decomposition(aligned)
    final_decomposition = [affinity.rotate(x, angle, origin=(0, 0)) for x in final_decomposition]
    final_decomposition = [x for x in final_decomposition if x.geom_type == "Polygon" and len(x.exterior.coords) <= 5]

    all_lines = List[DB.Line]()
    zvalue = spatial.Location.Point.Z if hasattr(spatial, "Location") and spatial.Location else 0.0

    for region in final_decomposition:
        region_lines = _polygon_to_lines(region, zvalue)
        Trace.Write(region_lines)
        for line in region_lines:
            all_lines.Add(line)

    print("Subdivided into {} regions and exported {} lines via Trace.Write.".format(len(final_decomposition), all_lines.Count))


if __name__ == "__main__":
    main()