# Guide Maître : Traitement de Nuages de Points avec Python > Lire une dalle LiDAR HD de l'IGN sans la télécharger, voxeliser, estimer les normales, segmenter au RANSAC et au DBSCAN, mesurer : la chaîne complete en Python, code execute et resultats mesures. - Adresse : https://fr.learngeodata.eu/blog/tutoriel-nuage-de-points-python - Rubrique : Articles techniques - Langue : fr - Auteur : Dr. Florent Poux (https://fr.learngeodata.eu/formateur/florent-poux) - Éditeur : 3D Geodata Academy (organisme de formation, France) - Index du site pour les assistants : https://fr.learngeodata.eu/llms.txt - Publié le 2026-01-10, mis à jour le 2026-09-10 --- Une dalle LiDAR HD de l'IGN, coloriée par altitude. C'est le point de départ de la chaîne décrite ici, et l'article la déroule sur cette scène exacte. ## 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 474 033 points effectivement lus, coloriés par leur classe ASPRS : sol (285 178), végétation basse, moyenne et haute (155 694 au total), bâtiment (28 856). Cette classification est livrée par l'IGN dans le fichier ; elle n'a pas été calculée ici. Dalle `LHD_FXX_0877_6382`, Lambert-93, altitudes IGN69. 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. | Côté du voxel | Points restants | Réduction | Mémoire (XYZ float64) | Durée | | --- | --- | --- | --- | --- | | brut | 474 033 | - | 11,4 Mo | - | | 10 cm | 447 103 | 6 % | 10,7 Mo | 204 ms | | 25 cm | 311 424 | 34 % | 7,5 Mo | 187 ms | | 50 cm | 127 881 | 73 % | 3,1 Mo | 60 ms | | 1 m | 38 592 | 92 % | 0,9 Mo | 29 ms | À 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 Trente mètres autour du bâtiment, même point de vue et même échelle. À gauche, les lignes de balayage du capteur : la densité varie du simple au triple selon l'angle de vue. À droite, la trame régulière produite par la voxelisation. C'est cette régularité, et non le gain de mémoire, qui compte pour les algorithmes de la suite. 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 ``` Les normales d'un scan de bureau (196 133 points, calculées en 0,11 s), représentées selon la convention usuelle : les trois composantes du vecteur deviennent les trois canaux rouge, vert et bleu. Toutes les surfaces de même orientation prennent donc la même couleur : le sol en vert, les deux familles de murs en bleu et en rouge. Une teinte franche signale une surface bien définie ; un bruit de couleur signale un endroit où la géométrie locale est ambiguë. 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. Un seul appel à `segment_plane` sépare le sol du reste : 11 371 points retenus sur 46 074, soit 24,7 %, en 0,02 s. Le plan trouvé a pour normale (0,00 ; 1,00 ; 0,03) : l'axe vertical de ce scan est *Y*, et le sol est horizontal à trois centièmes près. Remarquez que le pied de la chaise reste en violet : le seuil de 2 cm ne l'a pas absorbé. 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. Les 24 grappes trouvées par `DBSCAN` en 0,02 s, une couleur par grappe, le bruit (271 points) en gris clair. La chaise ressort d'un bloc, en orange : accoudoirs, assise et dossier forment un ensemble continu, donc une seule grappe. Aucun de ces objets n'a été annoté ; ils émergent de la seule géométrie. 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") ``` Les cinq plus grosses grappes et leur boîte orientée. La chaise (4 547 points) mesure **0,94 × 0,68 × 0,51 m** : des dimensions de chaise de bureau, obtenues sans qu'aucun modèle n'ait été entraîné. Remarquez que les boîtes ne sont pas alignées sur les axes de l'image : c'est exactement ce que l'`OBB` apporte sur l'`AABB`. | Boîte | Calcul | Coût | Dimensions rendues | Usage | | --- | --- | --- | --- | --- | | AABB | min / max par axe | O(n) | dépendantes du repère | test de collision, indexation, tuilage | | OBB | analyse en composantes principales | O(n) + décomposition 3×3 | propres à l'objet | Scan-to-BIM, métré, classement par gabarit | 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. | Étape | Entrée | Complexité | Durée mesurée | Ce qui limite | | --- | --- | --- | --- | --- | | Lecture COPC | 474 033 pts | O(k) sur la fenêtre | ≈ 10 s | le réseau, pas le CPU | | Voxelisation | 474 033 pts | O(n) | 60 ms | bande passante mémoire | | Normales (KD-Tree) | 196 133 pts | O(n log n) | 110 ms | CPU, parallélisable | | RANSAC (1 plan) | 46 074 pts | O(i · n) | 20 ms | le nombre d'itérations | | DBSCAN | 10 768 pts | O(n log n) avec index | 20 ms | la densité, pas le nombre | 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](https://fr.learngeodata.eu/blog/pointnet-tutoriel-python) et sur la [segmentation](https://fr.learngeodata.eu/blog/segmentation-nuage-de-points-python). ### Aller plus loin avec le Dr. Poux Cet article couvre la géométrie. La formation « Architecte IA de Systèmes 3D » couvre la chaîne entière (acquisition, apprentissage, déploiement) sur vos propres données, en 120 heures. [Ou recevez d'abord le guide complet, gratuitement →](https://fr.learngeodata.eu/guide-architecte-ia-3d)