script-volume-auto/CONCEPT.md

Estimation de volume 2.5D par nuage de points et sensibilité à la décimation spatiale

Document de référence méthodologique, destiné à être adapté en section « Matériel et méthodes » d'une publication. Il décrit le protocole implémenté par ce dépôt, les grandeurs mesurées, et — section 7 — les limites et biais connus qui doivent être déclarés ou corrigés avant publication.

---

1. Question de recherche

L'estimation du volume d'un amas (stock de matériaux, tas, dépôt) à partir d'un nuage de points aéroporté est une opération courante mais dont la robustesse à la densité de points est rarement quantifiée. La question posée ici est double :

1. Le volume estimé à partir d'un nuage LiDAR HD (IGN) est-il comparable à une mesure de référence ?
2. Dans quelle mesure varie-t-il lorsque la densité de points est dégradée de façon contrôlée ?

Le second point conditionne le premier : si le volume est stable sur deux ordres de grandeur de densité, l'estimation ne dépend pas de la campagne d'acquisition ; s'il dérive, tout chiffre publié doit être accompagné de la densité à laquelle il a été obtenu.

---

2. Grandeur mesurée

Le volume est calculé selon la méthode 2.5D par différence de surfaces rasterisées, telle qu'implémentée par CloudCompare 2.13.2 (commande -VOLUME).

Soit une direction de projection verticale d (par défaut Z), et une grille régulière de pas g dans le plan orthogonal à d. Pour chaque cellule c de la grille, on définit :

• z_ceil(c) : hauteur de la surface supérieure (le « plafond ») ;
• z_ground(c) : hauteur de la surface inférieure (le « sol ») ;

et le volume est la somme de Riemann

V = g² · Σ_{c ∈ C_match} ( z_ceil(c) − z_ground(c) )

où C_match est l'ensemble des cellules pour lesquelles les deux hauteurs sont définies. La méthode décompose V en un volume ajouté (cellules où z_ceil > z_ground) et un volume retiré (cellules où z_ceil < z_ground) ; la grandeur V rapportée est le volume net, c'est-à-dire leur différence signée. Cette distinction importe dès que les deux surfaces s'intersectent.

Il s'agit donc d'une quadrature par la règle du point milieu sur la fonction de hauteur h(x,y) = z_ceil − z_ground. L'erreur de quadrature est en O(g²·‖∇²h‖) : elle est nulle pour une surface plane, et croît avec la courbure. Cette méthode suppose que la surface est un graphe de (x,y) : elle ne peut représenter ni surplomb, ni cavité fermée. Cette hypothèse est admissible pour un amas de matériaux, et constitue une limite explicite du protocole.

2.1. Deux modes de définition du sol

Mode « plan de référence constant » (const_height). Un seul nuage est fourni ; il tient lieu de surface supérieure. La surface inférieure est le plan z = z_ref, avec z_ref par défaut égal à l'altitude minimale du nuage.

Mode « comparaison de deux nuages » (cloud_compare). Deux nuages sont fournis : la surface supérieure de l'amas, et une surface de base (socle, plan d'ajustement, terrain naturel avant dépôt). La surface inférieure est la seconde.

Résultat expérimental déterminant. Le mode const_height n'est valide que si le nuage ne contient que l'objet mesuré. Dès que le fichier inclut du terrain environnant, z_ref (altitude minimale globale) n'a plus aucun rapport avec la base réelle de l'amas, et le volume obtenu est arbitraire. Sur un amas industriel de ~150 m d'emprise, le mode const_height a donné ≈ 50 m³, contre ≈ 14 244 m³ en mode cloud_compare avec une base ajustée sous l'amas — soit un facteur 285. Une inspection de la carte de hauteurs a confirmé que les deux nuages ont bien la même emprise au sol et que le second chiffre est le bon. Toute publication utilisant le mode const_height doit justifier que le nuage a été préalablement segmenté.

---

3. Protocole de dégradation contrôlée

3.1. Estimation de la résolution native

La résolution spatiale native s₀ du nuage est estimée par la médiane des distances au plus proche voisin (structure KD-tree), calculée sur un échantillon de points du nuage. Pour un semis approximativement poissonnien de densité surfacique ρ (points·m⁻²), cette médiane vaut

median_NN = sqrt( ln 2 / (π ρ) ) ≈ 0.4697 / sqrt(ρ)

Elle constitue donc un estimateur direct de l'espacement inter-points.

Voir §7.1 : l'implémentation actuelle calcule cette médiane sur un sous-échantillon aléatoire de taille fixe, ce qui introduit un biais multiplicatif connu qu'il faut corriger.

3.2. Chaîne de décimation

On construit six niveaux de densité, indexés n = 0…5. Le pas spatial cible du niveau n est

s(n) = s₀ · 2ⁿ

Le niveau n est obtenu en appliquant un sous-échantillonnage spatial (contrainte de distance minimale s(n) entre points conservés ; -SS SPATIAL de CloudCompare) au nuage du niveau n−1, et non au nuage original. Ce chaînage garantit que chaque niveau est un sous-ensemble strict du précédent, donc que la comparaison isole bien l'effet d'une perte d'information monotone.

Point à ne pas mal interpréter. Le pas spatial double à chaque niveau. Comme les points échantillonnent une surface (dimension 2), le nombre de points est divisé par ≈ 4 par niveau. Au niveau 5 le pas est ×32 mais le nombre de points est de l'ordre de 1/1000 du niveau 0. Le nombre de points effectivement obtenu n'est jamais garanti a priori : -SS SPATIAL impose une distance minimale, pas un taux de décimation.

Le niveau 0 diffère selon le mode :

Mode	Niveau 0	Niveaux 1–5
const_height	nuage intégral, non décimé	s(n) = s₀·2ⁿ
cloud_compare	nuages appariés en densité au pas s₀	s(n) = s₀·2ⁿ, appliqué aux deux nuages

Appariement de densité (mode cloud_compare). Comparer une surface supérieure à 10 pts·m⁻² avec une base à 0.5 pt·m⁻² fait porter au résultat la résolution de la plus grossière, sans que cela soit visible. On calcule donc median_NN séparément pour chaque nuage, on retient s₀ = max(median_NN_top, median_NN_bottom), et on ramène le nuage le plus dense à ce pas par un sous-échantillonnage spatial unique. L'opération est un no-op sur le nuage le plus clairsemé (un filtre de distance minimale ne peut pas densifier). Le niveau 0 du mode cloud_compare est donc lui-même déjà décimé — les écarts relatifs sont exprimés par rapport à ce niveau.

La décimation des deux nuages est ensuite synchronisée (même s(n)). Si l'un des deux passe sous un seuil de 25 points, la progression s'arrête pour les deux : les deux surfaces sont nécessaires à une comparaison valide.

3.3. Grille de calcul adaptative

Le pas de la grille de quadrature du niveau n est

g(n) = s(n) · k,    k = 2 par défaut

Ce choix est essentiel et contre-intuitif. Une première version du protocole figeait g à la résolution du niveau 0 « pour garder les niveaux comparables ». C'est exactement l'inverse de ce qu'il faut faire : dès que la décimation devient plus grossière que cette grille fixe, chaque cellule cesse de contenir un point, et l'indicateur de couverture (§4) s'effondre mécaniquement (99 % → 58 % → 13 % → 3 % → 0.9 % → 0.2 % sur les cinq niveaux). Cet effondrement mesure un désalignement grille/densité, pas une perte d'information : il ne dit rien sur la robustesse du volume.

En liant g(n) à s(n), chaque niveau est intégré sur une grille adaptée à sa propre densité. La contrepartie, qui doit être déclarée, est que le volume V(n) varie alors sous l'effet de deux causes confondues : la perte de points, et le changement de pas de quadrature (erreur en O(g²)). Le protocole ne les sépare pas. Une expérience de contrôle à g fixé permet de borner la seconde.

3.4. Traitement des occlusions (remplissage de cellules vides)

Un nuage aéroporté comporte des cellules vides (occlusions, zones d'ombre). Laissées vides, elles sont exclues de C_match et le volume est sous-estimé. La commande -VOLUME de CloudCompare n'expose aucune option de remplissage : la vérification du code source (qCC/ccCommandRaster.cpp, CommandVolume25D::process) montre que la stratégie LEAVE_EMPTY y est codée en dur pour les deux surfaces, alors même que les constantes des autres stratégies existent dans le fichier et sont consommées par CommandRasterize::process.

Le protocole contourne cette limite en rasterisant et comblant chaque nuage séparément avant le calcul de volume (-RASTERIZE -GRID_STEP g -EMPTY_FILL … -OUTPUT_CLOUD), puis en soumettant à -VOLUME la grille comblée plutôt que le nuage original. Le nuage brut, non comblé, reste utilisé pour le chaînage de la décimation, afin que l'interpolation d'un niveau ne se propage pas au suivant.

Six stratégies sont exposées, reproduisant la boîte de dialogue « Compute Volume » de l'interface graphique : LEAVE_EMPTY, MIN_H, MAX_H, CUSTOM_H, INTERP (interpolation de Delaunay bornée par une longueur d'arête maximale L), KRIGING. La stratégie retenue par défaut est INTERP.

Rôle de L (longueur d'arête maximale). L borne la portée du remplissage : un trou de diamètre supérieur à L n'est pas comblé. C'est ce qui empêche l'interpolation de « recoudre » les concavités réelles du contour de l'amas, qui doivent rester vides. Vérification sur un nuage synthétique comportant une occlusion circulaire de rayon 0.6 m :

Configuration	Cellules appariées	Trou comblé
Sans remplissage	96.4 %	non (volume sous-estimé)
INTERP, L = 0.30 m (> diamètre du trou)	98.9 %	oui
INTERP, L = 0.08 m (< diamètre du trou)	96.4 %	non (identique au cas sans remplissage)

Le paramètre se comporte donc bien comme une borne de portée. Sa valeur par défaut est L(n) = 3·g(n).

Voir §7.3 : L(n) croissant avec n introduit un biais systématique qui doit être neutralisé pour l'étude de sensibilité (L fixé en absolu).

3.5. Paramètres exposés

Paramètre	Rôle	Défaut
s₀ (initialStep)	résolution native	médiane NN
k (gridStepMultiplier)	g(n) = s(n)·k	2
L (maxEdgeLength)	portée du remplissage INTERP	3·g(n)
z_ref	plan de référence (mode const_height)	altitude min du nuage
stratégie de remplissage	LEAVE_EMPTY…KRIGING	INTERP
projection par cellule	AVG / MIN / MAX / MED	AVG
direction verticale	X / Y / Z	Z
facteur, nombre de niveaux	fixés par le protocole	2, 5

Note d'implémentation : -PROJ (hauteur de cellule) n'existe que pour -RASTERIZE — -VOLUME a PROJ_AVERAGE_VALUE codé en dur. Comme le nuage soumis à -VOLUME est déjà une grille comblée à un point par cellule, le moyennage interne de -VOLUME est une passe-plat : c'est bien le -PROJ du remplissage qui détermine le résultat. -VERT_DIR en revanche est accepté par les deux commandes et leur est passé de façon cohérente.

---

4. Indicateurs rapportés

Pour chaque niveau n :

• V(n) : volume net (m³), et sa décomposition en volumes ajouté / retiré ;
• S(n) : surface projetée des cellules appariées (m²) ;
• N(n) : nombre de points de chaque nuage, et N(n)/N(0) ;
• g(n), s(n), L(n) ;
• matching cells % : proportion de cellules de la grille où les deux surfaces sont simultanément définies, avec sa décomposition en cellules « ground » et « ceil » manquantes ;
• l'écart relatif 100·(V(n) − V(0))/V(0), indicateur principal de robustesse.

Une carte de comparaison (différence de hauteur z_ceil − z_ground, dégradé bleu→rouge, cellules non appariées en gris) est produite par niveau. Elle est reconstruite en Python (binned_statistic_2d + matplotlib) à partir d'exports ASCII des grilles comblées, plutôt que par l'export GeoTIFF de CloudCompare (voir §6).

Tous les nuages exploitables des six niveaux sont assemblés dans un unique fichier .bin, comme entités distinctes (le format .bin de CloudCompare supporte nativement une hiérarchie de nuages ; il ne faut surtout pas fusionner les géométries), pour inspection visuelle a posteriori.

Voir §7.4 : matching cells % est, dans l'implémentation actuelle, largement neutralisé par l'étape de remplissage. Il ne doit pas être présenté comme l'indicateur de robustesse principal.

---

5. Résultat de référence

Sur l'amas industriel de ~150 m d'emprise (Amas 1 Nuage.las vs Plan d'ajustement NUAGE.las), en mode cloud_compare : V ≈ 14 244 m³, stable à ± 4 % jusqu'au niveau 4 (facteur 16 sur le pas spatial, nuage réduit à 129 points).

Ce résultat, s'il se confirme sur un échantillon d'amas plus large et avec les corrections de §7, est le message principal : le volume 2.5D d'un amas est une grandeur intégrale, donc peu sensible à la densité de points tant que la géométrie d'ensemble reste échantillonnée. Les erreurs locales de hauteur, positives et négatives, se compensent dans la somme. La densité conditionne la précision du contour et des détails de surface, pas celle du volume.

---

6. Contraintes d'implémentation (reproductibilité)

Ces points ne relèvent pas de la méthode mais conditionnent la reproductibilité de l'exécution. Ils ont été vérifiés empiriquement sur CloudCompare 2.13.2 (build officiel Windows et paquet apt Debian trixie) et contredisent parfois la documentation.

1. -VOLUME ne produit son rapport texte que si -AUTO_SAVE ON. Avec -AUTO_SAVE OFF, la commande s'exécute, journalise [2.5D VOLUME CALCULATION] finished, et n'écrit rien — sans erreur. Il faut réactiver l'auto-sauvegarde pour ce seul appel, et supprimer ensuite la grille *_HEIGHT_DIFFERENCE_*.bin produite au passage.
2. CloudCompare écrit ses sorties automatiques relativement au dossier du fichier chargé par -O, pas au répertoire de travail du processus. La recherche du rapport doit être récursive.
3. Les grands nombres sont formatés avec une virgule comme séparateur de milliers dans le rapport (Volume: 14,244.464657). Un parsing naïf tronque silencieusement à la virgule : 14 244 devient 14. Ce bug a réellement produit un volume erroné avant d'être identifié.
4. Le paquet Debian de CloudCompare n'embarque pas le plugin LAS/LAZ. Les fichiers LAS/LAZ/COPC sont convertis en ASCII XYZ via laspy en amont ; l'import ASCII fonctionne sans plugin.
5. L'export GeoTIFF (-RASTERIZE -OUTPUT_RASTER_Z) provoque une assertion et un core dump sur ce paquet (GDAL non compilé). Les cartes sont donc rendues en Python.
6. CloudCompare reste une application Qt même en -SILENT et exige un serveur X (Xvfb).
7. Le CLI est une machine à états séquentielle sur un seul processus : chaque étape logique est un appel indépendant avec son propre rechargement.
8. Chaque appel -VOLUME du protocole est journalisé, et le rapport texte brut de CloudCompare est conservé (raw/L{n}_volume_report.txt), ce qui permet de rejouer le parsing indépendamment.

---

7. Limites, biais connus et corrections requises

Cette section doit être lue avant toute publication de chiffres issus de ce pipeline. Les points 7.1 et 7.2 sont des défauts d'implémentation à corriger ; 7.3 à 7.6 sont des limites de protocole à déclarer ou à traiter.

7.1. Biais de sous-échantillonnage sur l'estimation de s₀ (à corriger)

median_NN est actuellement calculé sur un sous-échantillon aléatoire de taille fixe (50 000 points) tiré du nuage complet. La distance au plus proche voisin n'est pas invariante par sous-échantillonnage : prélever une fraction f divise la densité par f et multiplie median_NN par 1/√f.

Sur une dalle LiDAR HD de 20 M de points (~10 pts·m⁻², s₀ réel ≈ 0.15 m), f = 2.5·10⁻³ et l'estimation obtenue est ≈ 3 m, soit un facteur 20. Tout le protocole en hérite : g(0), L(0) et l'ensemble des s(n). Le biais est nul pour les nuages de moins de 50 000 points et croît en √N au-delà — il est donc invisible sur les petits jeux de test.

Correction. Soit corriger analytiquement, s₀ = median_NN_échantillon · √(N_échantillon / N_total) (le nombre total de points est déjà lu par le pipeline), soit — préférable car exact plutôt qu'asymptotique — calculer median_NN sur une emprise spatiale restreinte (un pavé de la bbox contenant ≤ 50 000 points), qui préserve la densité locale.

7.2. Précision numérique en coordonnées projetées (à corriger)

CloudCompare stocke les coordonnées en simple précision (float32). En Lambert-93 (EPSG:2154), les ordonnées valent ≈ 6.8·10⁶ m, où l'ULP de float32 est de 0.5 m. Sans décalage global (-GLOBAL_SHIFT, absent du pipeline), les coordonnées X/Y sont quantifiées sur un réseau de 0.5 m, ce qui vide de son sens toute grille de pas inférieur au mètre. CloudCompare n'émet qu'un avertissement de journal (Entity has very big coordinates: original accuracy may be lost!).

Le biais est passé inaperçu sur les jeux de test, dont les coordonnées semblent exprimées dans un repère local. Il est rédhibitoire sur du LiDAR HD IGN brut. Noter que seules X et Y sont affectées : Z ≈ 100–200 m reste représentable au micromètre.

Correction. Le volume 2.5D étant invariant par translation, soustraire un décalage fixe (coin de la bbox, arrondi) lors de la conversion LAS→ASCII, et consigner ce décalage dans le rapport.

7.3. Le remplissage devient plus agressif à chaque niveau (à traiter)

Avec le défaut L(n) = 3·g(n) = 6·s(n), la portée du remplissage double à chaque niveau. Au niveau 5, l'interpolation comble des lacunes 32 fois plus larges qu'au niveau 0, y compris des concavités réelles du contour — précisément ce que L était censé empêcher. Le volume dérive donc vers le haut avec n pour une raison étrangère à la décimation.

Pour l'étude de sensibilité, L doit être fixé en valeur absolue, identique sur les six niveaux (le paramètre existe déjà en surcharge). Toute courbe V(n) produite avec L(n) proportionnel à g(n) mélange deux effets et ne peut être publiée telle quelle.

7.4. matching cells % n'est plus un indicateur de couverture

Depuis que les nuages sont rasterisés et comblés au pas g(n) avant d'être soumis à -VOLUME, qui les rasterise à nouveau au même pas g(n), chaque cellule reçoit environ un point par construction. En mode const_height, le sol est un plan constant, défini partout. matching cells % ne mesure alors plus la couverture des données, mais essentiellement le désalignement d'une demi-cellule entre la grille du -RASTERIZE (ancrée sur la bbox du nuage brut) et celle du -VOLUME (ancrée sur la bbox de la grille comblée). Les valeurs de 96–99 % observées sont cohérentes avec cette interprétation.

L'indicateur reste utile avec la stratégie LEAVE_EMPTY, où il retrouve son sens originel. Il ne doit pas être présenté comme le critère de robustesse quand le remplissage est actif. L'indicateur de robustesse à publier est l'écart relatif V(n)/V(0) − 1, accompagné de sa dispersion (§7.6).

7.5. Co-registration des grilles en mode cloud_compare

Les deux nuages sont rasterisés indépendamment, donc chacun sur une grille ancrée sur sa propre bbox. Les centres de cellules des deux surfaces sont ainsi décalés d'une fraction de cellule non contrôlée. -VOLUME les reprojette ensuite sur une grille commune : deux points verticalement superposés peuvent tomber dans des cellules distinctes. L'erreur induite sur z_ceil − z_ground est d'ordre pente · g(n) et croît donc avec le niveau.

CloudCompare n'expose pas l'origine de la grille de rasterisation. La correction propre consiste à effectuer la rasterisation des deux nuages en Python sur une grille commune explicitement ancrée. À défaut, l'amplitude de l'effet doit être bornée expérimentalement (par exemple en translatant l'un des nuages d'une fraction de g et en observant la dispersion de V).

7.6. Absence de quantification d'incertitude (à traiter)

Le protocole produit une seule valeur de V par niveau, sans barre d'erreur. Le sous-échantillonnage spatial étant quasi déterministe (il dépend de l'ordre de parcours de l'octree), une simple répétition ne produirait aucune variance.

Deux sources de dispersion doivent être échantillonnées pour publier un intervalle :

• la réalisation de la décimation : répéter chaque niveau avec un sous-échantillonnage aléatoire à nombre de points équivalent, sur M tirages ;
• la phase de la grille : translater le nuage d'un décalage uniforme dans [0, g)² avant rasterisation, sur M tirages, ce qui échantillonne directement l'erreur de quadrature.

V(n) doit alors être rapporté comme moyenne ± écart-type sur ces M réalisations.

Par ailleurs, le pipeline n'est pas reproductible en l'état : le sous-échantillonnage aléatoire servant à estimer s₀ et z_ref n'est pas graîné, ces deux quantités varient donc d'une exécution à l'autre sur un même fichier. Une graine explicite doit être fixée et consignée.

7.7. Sensibilité de z_ref aux points aberrants

En mode const_height, z_ref est l'altitude minimale du nuage. Un unique point aberrant bas (bruit, écho multiple) abaisse z_ref de Δz et gonfle le volume de Δz · S. Sur une emprise de 10 000 m², un aberrant à 10 cm sous le sol ajoute 1 000 m³. Un percentile robuste (p₀.₁ ou p₁) doit être substitué au minimum, et la valeur retenue doit figurer dans le rapport. (Défaut supplémentaire de l'implémentation actuelle : z_ref est tiré de la bbox du sous-échantillon, ce qui le rend en outre biaisé vers le haut et aléatoire.)

7.8. Hypothèses structurelles du modèle 2.5D

• La surface est un graphe de (x,y) : ni surplomb, ni cavité fermée n'est représentable.
• L'emprise diagonale d'un amas dans sa bounding box (observée : ~150 × 113 m) n'invalide pas le calcul mais gonfle le nombre de cellules vides, donc l'influence de la stratégie de remplissage.
• La grandeur rapportée est le volume net. Dès que les deux surfaces s'intersectent, il faut présenter aussi les volumes ajouté et retiré ; le net seul est trompeur.
• Le format E57 est hors périmètre.

---

8. Synthèse du pipeline

LAS / LAZ / COPC.LAZ / BIN
  └─ conversion ASCII XYZ (laspy)                    [+ décalage global — §7.2]
  └─ estimation de s₀ = médiane NN (KD-tree)         [+ correction de biais — §7.1]
  └─ [cloud_compare] appariement de densité : s₀ = max(NN_top, NN_bottom)
  │
  └─ pour n = 0…5 :
       ├─ décimation spatiale au pas s(n) = s₀·2ⁿ, chaînée depuis le niveau n−1
       ├─ rasterisation + remplissage des occlusions au pas g(n) = 2·s(n)   [L fixe — §7.3]
       ├─ calcul de volume 2.5D sur les grilles comblées (CloudCompare -VOLUME)
       ├─ carte de comparaison z_ceil − z_ground (matplotlib)
       └─ arrêt si un nuage passe sous 25 points
  │
  └─ assemblage des 6 niveaux dans un .bin unique (entités distinctes)
  └─ rapport CSV / JSON + graphiques de synthèse

Résultat attendu : la courbe V(n) et son écart relatif à V(0), qui quantifient la robustesse de l'estimation de volume à la dégradation de densité — sous réserve des corrections de la section 7.