[Dynamo += Python] Subdivision de polygones en rectangles

Subdivision de polygones en rectangles : force brute ou algorithme d'optimisation ?

#DynamoBIM #Python #shapely #scipy #polygon #AutodeskExpertElite #AutodeskCommunity 

if this article is not in your language, use the Google Translate widget ⬈ (bottom of page for Mobile version ⬇)

Nous allons disséquer deux approches algorithmiques avancées pour découper automatiquement un polygone complexe en sous-rectangles (la décomposition soustractive gloutonne). 

L'un utilise une grille géométrique discrète, tandis que l'autre fait appel à l'intelligence computationnelle (algorithmes évolutionnaires via Scipy).

• Code 1 : L'approche par grille discrète (Grid-Based Search)

Ce premier script utilise les sommets du polygone pour tisser une grille virtuelle, puis utilise cette dernière pour générer les rectangles.

import clr
clr.AddReference('ProtoGeometry')
import Autodesk.DesignScript.Geometry as DS

import numpy as np
from itertools import combinations
from shapely import affinity
from shapely.geometry import Polygon, box
from shapely.ops import unary_union
from shapely.prepared import prep


def convert_poly_toDSPoly(geo):
    pts = [DS.Point.ByCoordinates(x, y, zvalue) for x, y in geo.exterior.coords]
    try:
        return DS.Polygon.ByPoints(DS.Point.PruneDuplicates(pts))
    except:
        return None


def get_largest_internal_rect(poly):
    rings = [poly.exterior] + list(poly.interiors)
    coords = np.vstack([np.asarray(r.coords)[:, :2] for r in rings])
    xs, ys = map(np.unique, coords.T)
    pp = prep(poly)

    rects = (
        box(x1, y1, x2, y2)
        for x1, x2 in combinations(xs, 2)
        for y1, y2 in combinations(ys, 2)
    )
    return max((r for r in rects if pp.contains(r)), key=lambda r: r.area, default=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

ds_polygon = IN[0]
v = DS.Vector.XAxis()
cs = DS.BoundingBox.ByMinimumVolume([ds_polygon]).ContextCoordinateSystem
angle = v.AngleWithVector(cs.XAxis) * (1 if v.Cross(cs.XAxis).Z > 0 else -1)
zvalue = ds_polygon.Points[0].Z

poly = Polygon([(p.X, p.Y) for p in ds_polygon.Points]).buffer(0)
poly = affinity.rotate(poly, -angle, origin=(0, 0))
# grid-based search algo
final_decomposition = subtractive_decomposition(poly)
final_decomposition = [affinity.rotate(x, angle, origin=(0, 0)) for x in final_decomposition]

OUT = [convert_poly_toDSPoly(x) for x in final_decomposition if x.geom_type == 'Polygon' and len(x.exterior.coords) <= 5]

Zoom sur les fonctions :

• convert_poly_toDSPoly(geo) : Transforme un polygone calculé en Python (shapely) en un vrai polygone exploitable graphiquement dans l'interface de Dynamo (Autodesk.DesignScript.Geometry).

• get_polys(geo) : Une fonction de sécurité pour s'assurer que si une opération de découpe sépare un morceau en deux, on récupère bien une liste propre de polygones individuels, en ignorant les résidus géométriques (lignes, points).

• get_largest_internal_rect(poly) : Elle extrait toutes les coordonnées X et Y uniques de tous les sommets du polygone pour former un maillage. Elle teste ensuite de manière exhaustive toutes les combinaisons de rectangles possibles formées par ce maillage (combinations), vérifie s'ils restent bien à l'intérieur du bâtiment ( prep(poly).contains ) et renvoie le plus grand en termes de surface.

• subtractive_decomposition(poly, max_rooms) : la boucle trie les morceaux restants par taille, appelle get_largest_internal_rect pour extraire le plus grand rectangle du plus grand morceau, et met à jour la forme restante via target.difference(rect) .
  
  

Résultat : L'algorithme privilégie la taille maximale de chaque rectangle.

• Code 2 : L'approche Differential Evolution (via scipy)

Le second script remplace la recherche par grille par une méthode d'optimisation à l'aide de la bibliothèque SciPy.

import clr
clr.AddReference('ProtoGeometry')
import Autodesk.DesignScript.Geometry as DS

import numpy as np
from itertools import combinations
from scipy.optimize import differential_evolution
from shapely import affinity
from shapely.geometry import Polygon, box
from shapely.ops import unary_union
from shapely.prepared import prep


def convert_poly_toDSPoly(geo):
    pts = [DS.Point.ByCoordinates(x, y, zvalue) for x, y in geo.exterior.coords]
    try:
        return DS.Polygon.ByPoints(DS.Point.PruneDuplicates(pts))
    except:
        return None


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


ds_polygon = IN[0]
v = DS.Vector.XAxis()
cs = DS.BoundingBox.ByMinimumVolume([ds_polygon]).ContextCoordinateSystem
angle = v.AngleWithVector(cs.XAxis) * (1 if v.Cross(cs.XAxis).Z > 0 else -1)
zvalue = ds_polygon.Points[0].Z

poly = Polygon([(p.X, p.Y) for p in ds_polygon.Points]).buffer(0)
poly = affinity.rotate(poly, -angle, origin=(0, 0))
# subtractive_decomposition via differential_evolution algo
final_decomposition = subtractive_decomposition(poly)
final_decomposition = [affinity.rotate(x, angle, origin=(0, 0)) for x in final_decomposition]

OUT = [convert_poly_toDSPoly(x) for x in final_decomposition if x.geom_type == 'Polygon' and len(x.exterior.coords) <= 5]

La principale modification vient de la fonction get_largest_internal_rect :

Au lieu de tester une grille fixe, cette fonction lâche une "population" de rectangles aléatoires (définis par leur centre cx, cy et leurs demi-dimensions hx, hy ).

• On utilise la méthode SciPy differential_evolution couplée à une fonction de perte ( obj(v) ) qui applique une pénalité si les rectangles dépassent le contour du bâtiment.

Résultat : L'algorithme privilégie la couverture totale maximale (remplir tout l'intérieur) plutôt que la surface de chaque rectangle individuel.

vidéo

« Le repos rend l'esprit plus libre et plus sain pour réfléchir »

George Sand

...