[Dynamo += Python] la bibliothèque SciPy

Saviez-vous que la librairie Python SciPy est inclus dans Dynamo ?

#DynamoBIM #Python #SciPy #RevitAPI     #AutodeskExpertElite #AutodeskCommunity 

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

Construite sur les fondations de la bibliothèque NumPy (qui est aussi incluse dans Dynamo), qui fournit un support pour les tableaux et les matrices multidimensionnelles, SciPy étend ces capacités avec des modules dédiés à des tâches scientifiques courantes.

Pour en bénéficier dans Dymamo, il suffit d'utiliser les nœuds Python CPython3 ou encore mieux, PythonNet3 (Dynamo 3+)

Module	Description	Exemples d'utilisation
scipy.linalg	Algèbre linéaire	Résolution de systèmes d'équations linéaires, calcul de valeurs et vecteurs propres, décompositions de matrices (SVD, LU, etc.).
scipy.optimize	Optimisation	Recherche de minima et de maxima de fonctions, ajustement de courbes (curve fitting), résolution de problèmes de programmation linéaire.
scipy.stats	Statistiques	Fonctions de densité de probabilité, tests statistiques (test t, ANOVA), statistiques descriptives, génération de nombres aléatoires.
scipy.integrate	Intégration numérique	Calcul d'intégrales définies, résolution d'équations différentielles ordinaires (EDO).
scipy.interpolate	Interpolation	Estimation de nouvelles valeurs à partir de points de données existants.
scipy.fft	Transformée de Fourier	Analyse de signaux dans le domaine fréquentiel, filtrage.
scipy.signal	Traitement du signal	Filtrage de signaux, analyse spectrale, conception de filtres.
scipy.sparse	Matrices creuses	Stockage et manipulation efficaces de matrices contenant une grande proportion de zéros.
scipy.spatial	Algorithmes spatiaux	Calcul de distances, triangulation de Delaunay, enveloppes convexes.
scipy.ndimage	Traitement d'images	Filtrage, morphologie, segmentation d'images n-dimensionnelles.

Nous ne pouvons pas passer en revus toutes les fonctionnalités de SciPy, mais voici quelques exemples de possibilités avec Dynamo sans rajout de bibliothèques supplémentaires.

 

• Courbe de Bézier
  

Dans cet exemple, nous créons/modifions une courbe de Bézier grâce à des points de contrôles

code

import sys
import clr
clr.AddReference('ProtoGeometry')
from Autodesk.DesignScript.Geometry import *

import numpy as np
from scipy.special import comb

def bernstein(n, k, t):
    """compute Bernstein polynomial at t using scipy comb."""
    return comb(n, k) * (t ** k) * ((1 - t) ** (n - k))

def bezier_curve(points, num_points=100):
    """compute a Bézier curve from control points using Bernstein polynomials."""
    n = len(points) - 1  # Degree of the Bézier curve
    t_values = np.linspace(0, 1, num_points)
    
    curve = np.zeros((num_points, 2))
    for i in range(n + 1):
        curve += np.outer(bernstein(n, i, t_values), points[i])
    
    return curve

ds_points = IN[0]
control_points = [[p.X, p.Y] for p in ds_points]
curve = bezier_curve(control_points)
out_points =[Point.ByCoordinates(x, y) for x, y in curve]
OUT = NurbsCurve.ByPoints(out_points)

• k-dimensional tree (KDTree)

L'objectif ici est de trouver (sur un plan 2D) le point le plus près d'une ligne grâce à l'algorithme k-dimensional tree. 

Note : scipy.spatial.KDTree ne fonctionne qu'en 2D

code

import sys
import clr
import System
clr.AddReference('ProtoGeometry')
from Autodesk.DesignScript.Geometry import *

import sysconfig
# standard library path
sys.path.append(sysconfig.get_path("platstdlib"))
# site-package library path
sys.path.append(sysconfig.get_path("platlib"))
import numpy as np
from scipy import spatial


def get_points_onCurve(curve, step = 0.2):
    out = []
    i = step
    while i <= curve.Length :
        try:
            po = curve.PointAtSegmentLength(i)
            i += step
            out.append([po.X, po.Y, po.Z])
        except:
            break
    return out
    
def get_groupIdx_by_flatIdx(idx):
    curr = 0
    for i, grp in enumerate(lstgroup_pts):
        for j , pt in enumerate(grp):
            if curr == idx:
                return i
            else:
                curr += 1

lstgroup_pts = IN[0]
lstLines = IN[1]

a_arr = np.array([[p.X, p.Y, p.Z] for group in lstgroup_pts for p in group])
ktree = spatial.KDTree(a_arr)

out = []
for l in lstLines:
    pts_arr = get_points_onCurve(l)
    distances, lstindex = ktree.query(pts_arr, k=1)
    counts = np.bincount(lstindex)
    idx = np.argmax(counts)
    # search the group index 
    gpIdx = get_groupIdx_by_flatIdx(idx)
    out.append(gpIdx)
OUT = out

• Trouver la courbe médiane entre 2 courbes
  

Voir article précédent

https://voltadynabim.blogspot.com/2024/05/dynamo-python-trouver-la-courbe-mediane.html

• identifier des clusters de points

L'objectif ici est de regrouper des points par proximité (clusters).

Nous créons un modèle de "clustering"  en utilisant la méthode cKDTree.query_ball_point (méthode similaire à DBSCAN de Scikit-learn).

code

import sys
import clr
clr.AddReference('ProtoGeometry')
from Autodesk.DesignScript.Geometry import *

clr.AddReference('GeometryColor')
from Modifiers import GeometryColor
clr.AddReference('DSCoreNodes') 
import DSCore
from DSCore import Color as DSColor

import numpy as np
import random
from scipy.spatial import cKDTree

random.seed(40)
inputs_points = IN[0]
array_pt = np.array([[p.X, p.Y] for p in inputs_points])
groups = []
    

    
def dbscan(x_pts, eps, min_samples):
    """
    dbscan “from scratch” with scipy and numpy

    returns
    -------
    labels : ndarray[int]
        each integer represent a label
        -2 = unvisited, -1 = noise, >=0 = cluster id
    """
    x_pts = np.asarray(x_pts)
    n = x_pts.shape[0]

    # precompute neighbors for each point
    tree = cKDTree(x_pts)
    neighbors = tree.query_ball_point(x_pts, eps)

    # initialization
    labels = np.full(n, -2, dtype=int)   # -2 = unvisited
    cluster_id = 0

    # scan through all points
    for i in range(n):
        if labels[i] != -2:
            continue                    # already visited
        if len(neighbors[i]) < min_samples:
            labels[i] = -1              # mark as noise
            continue

        # i is a core point → start a new cluster
        labels[i] = cluster_id
        stack = neighbors[i].copy()    # points to expand in this cluster

        # expand cluster
        while stack:
            j = stack.pop()
            if labels[j] == -1:
                labels[j] = cluster_id  # reclassify noise to border
            if labels[j] != -2:
                continue                # already processed
            labels[j] = cluster_id     # assign to cluster
            if len(neighbors[j]) >= min_samples:
                # j is a core point → add its neighbors
                stack.extend(neighbors[j])

        cluster_id += 1

    return labels
    
###### MAIN #####
eps = 3  
min_samples = 15

labels = dbscan(array_pt, eps, min_samples)
#print(labels)

cluster_indices = {c: np.where(labels==c)[0].tolist()
            for c in np.unique(labels) if c>=0}
noise_indices = np.where(labels==-1)[0].tolist()

# color results
for lab, idxs in cluster_indices.items():
    subpts = np.take(inputs_points, idxs, axis=0)
    col = DSColor.ByARGB(
        255,
        random.randint(0, 255),
        random.randint(0, 255),
        random.randint(0, 255)
    )
    for p in subpts:
        groups.append( GeometryColor.ByGeometryColor(p, col) )

# (optionnel) on affiche aussi les points bruit en gris
grey = DSColor.ByARGB(128, 180, 180, 180)
for i in noise_indices:
    groups.append( GeometryColor.ByGeometryColor(inputs_points[i], grey) )

OUT = groups
    

.....

• alpha-shape

Ici, nous utilisons la classe Delaunay pour générer une alpha-shape (une enveloppe polygonale 2D d'un ensemble de points).
Voir mon article précédent

import clr
import sys
clr.AddReference('ProtoGeometry')
from Autodesk.DesignScript.Geometry import *
import numpy as np
from scipy.spatial import Delaunay

# --- helper to compute all triangle circumradii ---
def circumradii(pts, simplices):
    pa, pb, pc = [pts[s] for s in simplices.T]
    a = np.linalg.norm(pb - pc, axis=1)
    b = np.linalg.norm(pa - pc, axis=1)
    c = np.linalg.norm(pa - pb, axis=1)
    # vectorized triangle area
    area = 0.5 * np.abs(
        (pb[:,0] - pa[:,0]) * (pc[:,1] - pa[:,1]) -
        (pb[:,1] - pa[:,1]) * (pc[:,0] - pa[:,0])
    )
    # avoid div-by-zero → inf circumradius
    return np.where(area > 0, a * b * c / (4.0 * area), np.inf)

def alpha_shape(pts, alpha_norm, percentile=99):
    # triangulate & get radii
    tri = Delaunay(pts)
    simp = tri.simplices
    R = circumradii(pts, simp)

    # pick 99th percentile as the “max radius” and clamp alpha_norm to [0,1]
    maxR = np.nanpercentile(R[np.isfinite(R)], percentile)
    alpha = np.clip(float(alpha_norm), 0.0, 1.0) * maxR

    # collect “small” triangles’ edges and count occurrences
    edges = {}
    for simplex, r in zip(simp, R):
        if r < alpha:
            for i,j in ((0,1),(1,2),(2,0)):
                e = tuple(sorted((simplex[i], simplex[j])))
                edges[e] = edges.get(e, 0) + 1

    # boundary edges appear only once
    boundary = [e for e,c in edges.items() if c == 1]

    # build adjacency and walk the loop
    adj = {}
    for i,j in boundary:
        adj.setdefault(i, []).append(j)
        adj.setdefault(j, []).append(i)

    start = boundary[0][0]
    loop = [start]
    prev = None
    curr = start
    j = 0
    while j < len(pts): # avoid infinite loop
        j += 1
        nbrs = adj[curr]
        nxt = nbrs[0] if nbrs[0] != prev else nbrs[1]
        if nxt == start:
            break
        loop.append(nxt)
        prev, curr = curr, nxt
    return loop

# read 
dynamo_pts = IN[0]
alpha_norm = IN[1]
pts = np.array([[p.X, p.Y] for p in dynamo_pts])

loop_idxs = alpha_shape(pts, alpha_norm)
contour_pts = np.take(dynamo_pts, loop_idxs, axis=0)
# output Dynamo points in the right order
OUT = contour_pts # [dynamo_pts[i] for i in loop_idxs]

...

• Résolution d'équations

Dans ce dernier exemple, nous utilisons la classe scipy.linalg.solve pour calculer l'intersection projetée de deux lignes consécutive.

Chaque ligne est convertie en sa forme générale d'équation cartésienne, qui est ax+by=c. 

Pour une ligne passant par les points (x1, y1) et (x2, y2), les coefficients sont :

a=y1 − y2

b=x2 − x 1

c=x2*y1 − y2*x1

Le point d'intersection (x, y) est le point unique qui satisfait les deux équations simultanément. On a donc un système de deux équations à deux inconnues (x et y).

Ce système peut être écrit sous forme matricielle comme Ax=b :

Il s'agit ici d'un problème simple (qui d'ailleurs peut être résolu sans scipy voir code), mais il est possible de résoudre des équations beaucoup plus complexes.

# Load the Python Standard and DesignScript Libraries
import sys
import clr
import numpy as np
from scipy.linalg import solve
clr.AddReference('ProtoGeometry')
from Autodesk.DesignScript.Geometry import *

def compute_intersection_scipy(lineA, lineB):
    """Compute intersection point of 2 lines using SciPy's linear solver."""
    x1, y1 = lineA.StartPoint.X, lineA.StartPoint.Y
    x2, y2 = lineA.EndPoint.X, lineA.EndPoint.Y
    x3, y3 = lineB.StartPoint.X, lineB.StartPoint.Y
    x4, y4 = lineB.EndPoint.X, lineB.EndPoint.Y
    # linear equation of lines
    # Line A: (y1 - y2)x + (x2 - x1)y = x2*y1 - y2*x1
    # Line B: (y3 - y4)x + (x4 - x3)y = x4*y3 - y4*x3
    # create equilatent of matrix system : Ax = b
    A_matrix  = np.array([
        [y1 - y2, x2 - x1],
        [y3 - y4, x4 - x3]
    ])

    b_vector  = np.array([
        x2*y1 - y2*x1,
        x4*y3 - y4*x3
    ])
    #
    # check if lines are parallals / coincident
    determinant = np.linalg.det(A_matrix)
    if abs(determinant) < 0.01:
        return None, None
    else:
        # Solve the linear system Ax = b
        px, py  = solve(A_matrix, b_vector ) # intersection_point
        return px, py



def compute_intersection(lineA, lineB):
    """alternative computation intersection point of 2 lines without scipy"""
    x1, y1 = lineA.StartPoint.X, lineA.StartPoint.Y
    x2, y2 = lineA.EndPoint.X, lineA.EndPoint.Y
    x3, y3 = lineB.StartPoint.X, lineB.StartPoint.Y
    x4, y4 = lineB.EndPoint.X, lineB.EndPoint.Y

    denom = (x1 - x2)*(y3 - y4) - (y1 - y2)*(x3 - x4)
    print(denom)
    if abs(denom) < 0.001:
        return None, None  # Lines are parallel

    px = ((x1*y2 - y1*x2)*(x3 - x4) - (x1 - x2)*(x3*y4 - y3*x4)) / denom
    py = ((x1*y2 - y1*x2)*(y3 - y4) - (y1 - y2)*(x3*y4 - y3*x4)) / denom
    return float(px), float(py)

def ordered_curves(lines):
    ordered_curves = [lines.pop(0)]
    n = 1
    while lines and n < 200:
        n += 1
        lines.sort(key = lambda x : ordered_curves[-1].DistanceTo(x))
        ordered_curves.append(lines.pop(0))
    return ordered_curves

lines = IN[0]
perimeter_points = []
lines_ordered = ordered_curves(lines)
lines_ordered.append(lines_ordered[0])
for idx, lineB in enumerate(lines_ordered):
    if idx > 0:
        lineA = lines_ordered[idx - 1]
        x, y = compute_intersection_scipy(lineA, lineB)
        if x is None:
            perimeter_points.extend([lineA.StartPoint, lineA.EndPoint, lineB.StartPoint, lineB.EndPoint])
        else:
            perimeter_points.append(Point.ByCoordinates(x, y, lineB.StartPoint.Z))

polygon = Polygon.ByPoints(perimeter_points)

OUT = polygon, perimeter_points

...

 

« Un jour j'irai vivre en Théorie, car en Théorie tout se passe bien. »

Pierre Desproges