Guide Maître : Traitement de Nuages de Points avec Python

1. Préparer l'environnement

Un traitement 3D dépend de bibliothèques compilées en C++, qui se contredisent facilement d'une version à l'autre. Un environnement isolé n'est pas une précaution de confort : c'est ce qui rend un résultat reproductible six mois plus tard. mamba résout le graphe de dépendances nettement plus vite que conda, pour un comportement identique.

# Un environnement dédié au calcul 3D
mamba create -n spatial-ai python=3.11 -y
mamba activate spatial-ai

# laspy[lazrs] : le décodeur LAZ écrit en Rust, et le support COPC
pip install open3d "laspy[lazrs]" numpy matplotlib

L'extension lazrs mérite une ligne à elle seule. Sans ce décodeur, laspy retombe sur une implémentation Python de la décompression .LAZ, et la lecture devient le poste le plus coûteux de toute la chaîne : devant les algorithmes eux-mêmes.

L'environnement est prêt. Reste la question qui décide de tout le reste : comment lire des fichiers que la mémoire ne peut pas contenir ?

2. Lire sans tout charger

Voici le problème dans sa forme concrète. La dalle LHD_FXX_0877_6382 du programme LiDAR HD de l'IGN couvre un kilomètre carré : elle pèse 229 Mo compressée et contient 39 025 611 points. Décompressée en float64, elle occuperait près d'un gigaoctet de mémoire vive : pour une seule dalle, sur les quelque 500 000 que compte la couverture nationale.

La réponse tient dans le format. Un fichier .copc.laz (Cloud Optimized Point Cloud) est un LAZ ordinaire dont les points sont rangés dans un octree écrit à l'intérieur du fichier. On peut donc demander une emprise, ou un niveau de détail, et ne rapatrier que les octets correspondants : y compris à travers le réseau, sans jamais télécharger le fichier.

import numpy as np
import laspy
from laspy import CopcReader

URL = ("https://data.geopf.fr/telechargement/download/LiDARHD-NUALID/"
       "NUALHD_1-0__LAZ_LAMB93_OO_2025-03-14/"
       "LHD_FXX_0877_6382_PTS_LAMB93_IGN69.copc.laz")

with CopcReader.open(URL) as lecteur:
    # L'en-tête seul suffit à connaître l'emprise : quelques kilo-octets.
    entete = lecteur.header
    print(entete.point_count, entete.mins, entete.maxs)

    # Une fenêtre de 140 m de côté autour d'une ferme, en Lambert-93.
    cx, cy, rayon = 877132.7, 6381971.0, 70.0
    emprise = laspy.copc.Bounds(mins=np.array([cx - rayon, cy - rayon]),
                                maxs=np.array([cx + rayon, cy + rayon]))
    points = lecteur.query(emprise)

# laspy applique déjà scale et offset : x, y, z sont en mètres, en float64.
xyz = np.vstack((points.x, points.y, points.z)).transpose()
classes = np.asarray(points.classification)
print(xyz.shape, np.unique(classes))

Ce passage rend 474 033 points (un centième de la dalle) en une dizaine de secondes, sans qu'aucun fichier ne touche le disque. C'est la différence entre un pipeline qui tient sur un portable et un pipeline qui exige une station de travail.

Les points sont en mémoire. Ils sont encore beaucoup trop nombreux pour ce qui va suivre.

3. Voxelisation : réduire sans déformer

Le sous-échantillonnage par voxels ne supprime pas des points au hasard. Il découpe l'espace en cubes de côté fixe et remplace le contenu de chaque cube occupé par son barycentre. Deux conséquences : la densité devient uniforme (un capteur en produit une très inégale, dense sous le vol, clairsemée en bord de bande) et le nombre de points cesse de dépendre du nombre de mesures pour ne plus dépendre que du volume occupé.

import open3d as o3d

# Le passage de NumPy à Open3D. C'est ici que le nuage devient un objet
# manipulable : recentré pour rester en précision confortable.
centre = xyz.mean(axis=0)
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(xyz - centre)

reduit = pcd.voxel_down_sample(voxel_size=0.5)
print(f"{len(pcd.points)} -> {len(reduit.points)} points")
# 474033 -> 127881 points

La taille du voxel est le seul réglage, et c'est un arbitrage explicite entre le détail conservé et le coût de tout ce qui suit. Les valeurs ci-dessous ont été mesurées sur la fenêtre de 474 033 points de la section précédente.

À 10 cm, la réduction n'est que de 6 % : le LiDAR HD délivre environ dix points au mètre carré, donc la plupart des voxels de 10 cm ne contiennent déjà qu'un seul point. Le réglage ne devient utile qu'à partir du moment où le voxel est plus grand que l'espacement moyen des points. C'est la règle à retenir, et elle explique pourquoi une même valeur ne se transpose pas d'un jeu de données à l'autre.

Avant · 26 390 points

Après, voxels de 50 cm · 5 659 points

Le nuage est réduit et régulier. Pour aller plus loin, il faut donner à chaque point une notion de la surface à laquelle il appartient.

4. Les normales : donner un sens aux surfaces

Une normale est un vecteur unitaire perpendiculaire à la surface locale. Elle ne se mesure pas : elle se calcule, à partir des voisins de chaque point. Open3D construit un KD-Tree, prend les voisins dans un rayon donné, en calcule la matrice de covariance, et retient le vecteur propre associé à la plus petite valeur propre : la direction dans laquelle les voisins s'étalent le moins, c'est-à-dire la perpendiculaire au plan qu'ils dessinent.

# On change de jeu de données ici, et ce n'est pas un détail : les normales
# et la segmentation par plans se lisent sur des objets, pas sur une forêt
# vue d'avion. Ce scan de bureau est livré avec Open3D.
piece = o3d.io.read_point_cloud(o3d.data.DemoCropPointCloud().point_cloud_path)
print(len(piece.points))  # 196133

# Recherche hybride : un rayon de 10 cm, mais au plus 30 voisins.
# Le rayon donne une échelle physique, le plafond borne le coût.
recherche = o3d.geometry.KDTreeSearchParamHybrid(radius=0.1, max_nn=30)
piece.estimate_normals(search_param=recherche)

# Le calcul donne une direction, pas un sens : la normale peut pointer
# vers l'intérieur comme vers l'extérieur. On tranche explicitement.
piece.orient_normals_to_align_with_direction(orientation_reference=[0, 0, 1])

normales = np.asarray(piece.normals)  # (n, 3), norme 1

Chaque point sait maintenant à quelle surface il appartient. Il reste à découper la scène.

5. RANSAC : trouver le sol

RANSAC (Random Sample Consensus) procède par paris répétés : il tire trois points au hasard, en déduit un plan, compte combien de points du nuage tombent à moins d'un seuil de ce plan, et recommence. Après quelques milliers de tirages, le plan qui a rassemblé le plus large consensus est retenu. La méthode ne suppose rien sur le bruit et supporte qu'une majorité de points n'appartiennent pas au modèle : ce qui en fait l'outil de choix pour extraire un sol ou une façade.

# RANSAC tire au hasard : sans graine, deux exécutions ne donnent pas
# exactement le même compte. Toutes les valeurs de cet article sont
# celles de la graine 42.
o3d.utility.random.seed(42)

scene = piece.voxel_down_sample(voxel_size=0.015)  # 46074 points

modele, inliers = scene.segment_plane(
    distance_threshold=0.02,   # tolérance : 2 cm autour du plan
    ransac_n=3,                 # trois points définissent un plan
    num_iterations=2000,       # nombre de tirages
)
a, b, c, d = modele            # équation ax + by + cz + d = 0

sol = scene.select_by_index(inliers)
reste = scene.select_by_index(inliers, invert=True)
print(f"sol : {len(sol.points)} points, normale ({a:.2f}, {b:.2f}, {c:.2f})")

Le distance_threshold est le paramètre sensible : il doit correspondre à l'épaisseur réelle de la surface, c'est-à-dire au bruit du capteur. Trop serré, le sol se fragmente ; trop large, RANSAC « mange » le bas des murs et les premiers centimètres des objets posés dessus.

Le sol est isolé. Comment séparer maintenant les objets restants les uns des autres ?

6. DBSCAN : séparer par la densité

Contrairement au K-Means, DBSCAN n'a pas besoin qu'on lui annonce le nombre de groupes : il les découvre. Son principe est le contact : deux points assez proches appartiennent au même groupe, et un groupe ne se déclare qu'à partir d'un nombre minimal de points, ce qui évacue le bruit isolé au passage.

# Retirer le sol ne suffit pas en intérieur : les murs touchent tout ce
# qui est posé contre eux. On enlève les plans dominants tant qu'ils pèsent.
objets = reste
while len(objets.points) > 3000:
    _, inl = objets.segment_plane(0.02, 3, 2000)
    if len(inl) < 1500:
        break
    objets = objets.select_by_index(inl, invert=True)

etiquettes = np.array(objets.cluster_dbscan(eps=0.05, min_points=12))
# -1 signale le bruit ; les grappes sont numérotées à partir de 0.
print(f"{etiquettes.max() + 1} grappes, {(etiquettes < 0).sum()} points de bruit")

Cette boucle mérite qu'on s'y arrête, parce qu'elle corrige l'erreur la plus fréquente sur ce sujet. En appliquant DBSCAN directement après le retrait du seul sol, on obtient une unique grappe de 27 183 points : les murs relient physiquement la chaise, le meuble et les montants entre eux, et un algorithme qui regroupe par contact les considère à juste titre comme un seul objet. Après retrait des sept plans dominants, il reste 10 768 points et 24 grappes : les objets réels de la scène.

Chaque objet est isolé. Quelles sont ses dimensions ?

7. Boîtes englobantes : passer à la mesure

Deux boîtes existent, et le choix entre elles n'est pas cosmétique. L'AABB (Axis-Aligned Bounding Box) s'obtient par un simple minimum et maximum sur chaque axe : elle est immédiate, mais ses dimensions dépendent de l'orientation de l'objet par rapport au repère. L'OBB (Oriented Bounding Box) aligne d'abord la boîte sur les axes principaux de l'objet, par analyse en composantes principales, et rend donc une longueur, une largeur et une hauteur qui ont un sens métier.

points_objets = np.asarray(objets.points)

for i in range(etiquettes.max() + 1):
    grappe = o3d.geometry.PointCloud()
    grappe.points = o3d.utility.Vector3dVector(points_objets[etiquettes == i])
    if len(grappe.points) < 100:
        continue
    boite = grappe.get_oriented_bounding_box()
    L, l, h = boite.extent
    print(f"grappe {i} : {len(grappe.points)} pts, {L:.2f} x {l:.2f} x {h:.2f} m")

Ces dimensions sont exploitables telles quelles : un objet de 0,90 × 0,45 × 1,80 m debout contre un mur est un radiateur ou une armoire, pas un véhicule. C'est le premier étage d'une classification, et il ne coûte rien.

8. Du nuage au tenseur : le passage vers PyTorch

Les réseaux qui opèrent sur des points (PointNet++, Point Transformer, PTv3) attendent des tenseurs, pas des objets Open3D. La conversion est simple, mais elle est le lieu de deux erreurs coûteuses : la copie inutile et la précision perdue.

import torch

appareil = "cuda" if torch.cuda.is_available() else "cpu"

def vers_tenseur(pcd, appareil):
    # np.asarray ne copie pas : il ouvre une vue sur la mémoire d'Open3D.
    points = np.asarray(pcd.points)          # float64, (n, 3)

    # torch.from_numpy ne copie pas non plus : le tenseur partage ce tampon.
    # En revanche .float() copie : float64 vers float32 change la taille.
    # C'est voulu : les réseaux travaillent en float32, et le nuage a déjà
    # été recentré, donc sept chiffres significatifs suffisent.
    return torch.from_numpy(points).float().to(appareil)

tenseur = vers_tenseur(reduit, appareil)
print(tenseur.shape, tenseur.dtype, tenseur.device)

Le commentaire du milieu corrige un raccourci que l'on lit partout. torch.from_numpy ne copie effectivement rien, mais l'appel à .float() qui suit, lui, copie l'intégralité du tableau. Le gain est ailleurs : il est dans le fait de n'écrire qu'une seule copie, au lieu des deux ou trois qu'enchaîne un torch.tensor(np.array(...), dtype=...) écrit sans y penser.

La chaîne est complète. Que retenir de son coût ?

9. Ce que coûte chaque étape

Voici le bilan du passage complet dont sont issues toutes les illustrations de cet article. Les durées proviennent d'un seul cœur CPU, sur les jeux indiqués ; elles servent à comparer les ordres de grandeur entre étapes, pas à établir un record.

Une conclusion s'impose, et elle est contre-intuitive : sur une chaîne géométrique correctement écrite, l'entrée-sortie domine tout le reste de deux ordres de grandeur. Les algorithmes ne sont pas le goulot d'étranglement ; la lecture l'est. C'est là qu'il faut investir, et c'est précisément ce que résout un format indexé comme le COPC.

Conclusion

Nous sommes partis d'un fichier de 229 Mo posé sur un serveur et nous avons obtenu, sans le télécharger, des objets identifiés et mesurés au centimètre. Aucun modèle n'a été entraîné : tout repose sur de la géométrie, et c'est ce socle qui rend ensuite l'apprentissage profond utile plutôt que nécessaire.

• Lecture : un format indexé change la nature du problème. Cherchez le .copc.laz avant de chercher une plus grosse machine.
• Précision : restez en float64 en coordonnées absolues, recentrez avant de passer en float32.
• Voxelisation : elle vaut pour la régularité qu'elle apporte, davantage que pour la mémoire qu'elle économise.
• Segmentation : RANSAC pour les surfaces planes, DBSCAN pour ce qui reste, et retirez tous les plans dominants avant de densifier, sinon tout se touche.
• Mesure : une boîte orientée par grappe, jamais une pour l'ensemble.

La suite naturelle est la sémantique : nommer ce que la géométrie a séparé. C'est le sujet des articles sur PointNet et sur la segmentation.