Fiabilise le pipeline de volume et corrige plusieurs failles de securite

Pipeline: corrige le biais de sous-echantillonnage sur median_nn (facteur
sqrt(total/echantillon), gonflait le pas spatial jusqu'a x20 sur les gros
nuages), applique un decalage global X/Y pour la precision float32 en
coordonnees projetees (Lambert-93), fixe la portee de remplissage des trous
sur la grille du niveau 0 (au lieu de croitre a chaque niveau), utilise un
percentile robuste pour zref au lieu du minimum brut, et filtre plus
precisement le fichier de grille parasite genere par -VOLUME.

Securite: valide l'id de run (uuid) sur toutes les routes avant de construire
un chemin filesystem (traversee de repertoire post-auth via zip/csv/bin/images),
assainit le nom de fichier uploade avant multer, retire enableCors() (surface
inutile), passe le conteneur en utilisateur non-root.

Documente l'ensemble de la methode et ses limites dans CONCEPT.md.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
This commit is contained in:
Nicolas Fryder 2026-07-10 18:13:27 +02:00
parent 051209b022
commit 9d6ebf5177
12 changed files with 557 additions and 27 deletions

407
CONCEPT.md Normal file
View file

@ -0,0 +1,407 @@
# 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) )
```
`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 `n1`**,
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 15 |
|---|---|---|
| `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 ≈ 100200 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 9699 % 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 n1
├─ 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.

View file

@ -52,10 +52,24 @@ COPY public ./public
COPY docker/entrypoint.sh /entrypoint.sh
RUN chmod +x /entrypoint.sh && mkdir -p /data
# Utilisateur non-root : l'image debian:trixie-slim + nodejs (nodesource) n'a pas d'utilisateur
# "node" preexistant (contrairement a l'image officielle node:*). /app et /data sont chown AVANT
# la directive VOLUME pour que l'appartenance soit reprise par le volume nomme Coolify/Compose a
# sa premiere creation (Docker copie le contenu existant du repertoire dans un volume vide neuf).
# /tmp/.X11-unix doit exister et etre world-writable (mode 1777, meme convention que /tmp) AVANT
# de passer en non-root : Xvfb refuse de creer ce repertoire lui-meme si euid != 0
# ("_XSERVTransmkdir: ERROR: euid != 0"), il ne peut alors plus creer son socket X11 - verifie
# empiriquement, le healthcheck HTTP passe quand meme (Xvfb echoue silencieusement en arriere-plan)
# mais tout appel CloudCompare plante ensuite faute de display.
RUN useradd --create-home --uid 1000 --shell /usr/sbin/nologin appuser \
&& mkdir -p /tmp/.X11-unix && chmod 1777 /tmp/.X11-unix \
&& chown -R appuser:appuser /app /data
VOLUME ["/data"]
EXPOSE 3000
HEALTHCHECK --interval=30s --timeout=5s --start-period=30s --retries=3 \
CMD curl -fsS "http://localhost:${PORT:-3000}/health" || exit 1
USER appuser
ENTRYPOINT ["/entrypoint.sh"]

View file

@ -7,11 +7,21 @@ ne fournit pas le plugin LAS/LAZ - seulement le plugin "Core I/O").
Un .copc.laz est simplement lu comme un LAZ standard (les points sont
identiques, seule l'indexation octree COPC, inutile ici, est ignoree).
Decalage global X/Y : CloudCompare stocke les coordonnees en float32 en interne. En coordonnees
projetees (ex. Lambert-93 / EPSG:2154, Y ~ 6.8e6 m), l'ULP de float32 a cette magnitude vaut ~0.5m,
ce qui quantifie silencieusement la position des points sur une grille de 0.5m - incompatible avec
un pas de calcul de l'ordre du decimetre (voir CONCEPT.md §7.2). On soustrait donc un decalage fixe
(coin bas de la bbox, arrondi au metre, lu directement dans l'entete LAS - pas de second passage sur
les points) de X et Y avant ecriture. Le volume 2.5D etant invariant par translation, ce decalage
n'affecte aucun resultat ; il est renvoye en sortie pour tracabilite. Z n'est pas decale : en
alti-metrique (~10^2-10^3 m), la precision float32 y est deja largement suffisante.
Usage:
python3 las_to_xyz.py --input nuage.laz --output nuage.xyz
"""
import argparse
import json
import math
import sys
import laspy
@ -25,10 +35,16 @@ def main() -> int:
args = parser.parse_args()
try:
with laspy.open(args.input) as reader:
offset_x = math.floor(reader.header.mins[0])
offset_y = math.floor(reader.header.mins[1])
total = 0
with laspy.open(args.input) as reader, open(args.output, "w", encoding="ascii") as out:
with open(args.output, "w", encoding="ascii") as out:
for chunk in reader.chunk_iterator(2_000_000):
xyz = np.column_stack([np.asarray(chunk.x), np.asarray(chunk.y), np.asarray(chunk.z)])
x = np.asarray(chunk.x) - offset_x
y = np.asarray(chunk.y) - offset_y
xyz = np.column_stack([x, y, np.asarray(chunk.z)])
np.savetxt(out, xyz, fmt="%.6f")
total += xyz.shape[0]
except Exception as exc:
@ -39,7 +55,7 @@ def main() -> int:
print(json.dumps({"error": "le fichier LAS/LAZ ne contient aucun point"}))
return 1
print(json.dumps({"ok": True, "count": total}))
print(json.dumps({"ok": True, "count": total, "offsetX": offset_x, "offsetY": offset_y}))
return 0

View file

@ -44,6 +44,12 @@ def main() -> int:
"zmin": float(np.min(pts[:, 2])),
"zmax": float(np.max(pts[:, 2])),
}
# Percentiles bas de Z : contrairement a zmin (sensible a un unique point aberrant, ex. bruit
# multi-echo), un percentile est un estimateur robuste du plan de reference et n'est PAS biaise
# par le sous-echantillonnage (la distribution marginale de Z ne depend pas de la densite de
# points, au contraire de la distance au plus proche voisin). Voir CONCEPT.md §7.7.
z_p01 = float(np.percentile(pts[:, 2], 0.1))
z_p1 = float(np.percentile(pts[:, 2], 1.0))
tree = cKDTree(pts)
# k=2 : le plus proche voisin de chaque point est le point lui-meme (distance 0),
@ -62,6 +68,8 @@ def main() -> int:
"median_nn": float(np.median(nn_dist)),
"std_nn": float(np.std(nn_dist)),
"bbox": bbox,
"z_p01": z_p01,
"z_p1": z_p1,
"approximate": True,
}
print(json.dumps(result))

View file

@ -87,7 +87,7 @@ export class CcRunnerService {
* Recherche recursive : CloudCompare ecrit les rapports/rasters relativement au dossier du
* nuage charge via -O (ex: "levels/"), pas forcement au cwd du process.
*/
async findLatestFile(rootDir: string, opts: { prefix?: string; suffix?: string; after: number }): Promise<string | null> {
async findLatestFile(rootDir: string, opts: { prefix?: string; suffix?: string; mustContain?: string; after: number }): Promise<string | null> {
const candidates: Array<{ full: string; mtime: number }> = [];
const walk = async (dir: string): Promise<void> => {
@ -100,6 +100,7 @@ export class CcRunnerService {
}
if (opts.prefix && !entry.name.startsWith(opts.prefix)) continue;
if (opts.suffix && !entry.name.endsWith(opts.suffix)) continue;
if (opts.mustContain && !entry.name.includes(opts.mustContain)) continue;
const stat = await fs.stat(full);
const mtime = stat.mtimeMs;
if (mtime < opts.after - 2000) continue; // marge de securite horloge

View file

@ -21,8 +21,10 @@ export const config = {
// dynamiquement a chaque niveau, pas figee sur la resolution du niveau 0 - voir CLAUDE.md).
gridStepMultiplier: 2,
// Remplissage des trous (occlusions de scan) avant -VOLUME : distance max d'interpolation
// Delaunay = gridStep du niveau x ce multiplicateur. Empeche de combler les parties concaves
// du contour de l'amas (qui doivent rester "vides") tout en bouchant les petits trous de scan.
// Delaunay = gridStep DU NIVEAU 0 (initialStep x gridStepMultiplier) x ce multiplicateur, FIXE
// sur les 6 niveaux (pas recalcule a partir du gridStep de chaque niveau - sinon la portee de
// remplissage double a chaque niveau et finit par combler des concavites reelles du contour de
// l'amas, biais systematique croissant avec la decimation, voir CONCEPT.md §7.3).
// Voir CLAUDE.md : -VOLUME n'a PAS d'option de remplissage native, on passe par -RASTERIZE.
maxEdgeLengthMultiplier: 3,
// Defauts des parametres CloudCompare exposes dans l'UI (mirroir de la boite de dialogue

View file

@ -6,8 +6,10 @@ import { basicAuthMiddleware } from './auth/basic-auth.middleware';
import { config } from './config';
async function bootstrap() {
// Pas de CORS : app privee mono-utilisateur, servie et consommee depuis la meme origine (UI
// statique + API sur le meme port). enableCors() ferait aussi court-circuiter le preflight
// OPTIONS avant le middleware d'auth ci-dessous - surface inutile a garder ouverte.
const app = await NestFactory.create(AppModule, { bodyParser: true });
app.enableCors();
app.use(basicAuthMiddleware());
await app.listen(config.port, '0.0.0.0');
new Logger('Bootstrap').log(`Application demarree sur http://0.0.0.0:${config.port} (data dir: ${config.dataDir})`);

View file

@ -95,7 +95,7 @@ export class PipelineService {
// ---------- 1) Statistiques ----------
onProgress('Estimation du pas spatial (echantillonnage + KD-tree)...');
const t0 = Date.now();
const primaryStats = await this.computeStatsFor(primaryNorm, workDir, 'sample_primary.xyz');
const primaryStats = await this.computeStatsFor(primaryNorm.path, workDir, 'sample_primary.xyz');
let initialStep: number;
let zref: number | null = null;
@ -104,14 +104,17 @@ export class PipelineService {
if (!isCompare) {
initialStep = params.initialStepOverride ?? primaryStats.stats.median_nn;
zref = params.zrefOverride ?? primaryStats.stats.bbox.zmin;
// z_p1 (1er percentile) plutot que bbox.zmin : un minimum absolu est ecrase par un seul point
// aberrant (bruit multi-echo) et gonfle le volume de Δz x surface. Le percentile n'est en plus
// pas biaise par le sous-echantillonnage (contrairement a median_nn) - voir CONCEPT.md §7.7.
zref = params.zrefOverride ?? primaryStats.stats.z_p1;
zrefSource = params.zrefOverride !== undefined ? 'override' : 'auto';
onProgress(
`Pas initial=${initialStep.toFixed(5)}m (${params.initialStepOverride ? 'override' : 'auto, mediane NN'}), ` +
`Zref=${zref.toFixed(4)} (${zrefSource}). [${Date.now() - t0}ms]`,
`Zref=${zref.toFixed(4)} (${zrefSource}, percentile 1%). [${Date.now() - t0}ms]`,
);
} else {
const secondaryStats = await this.computeStatsFor(secondaryNorm!, workDir, 'sample_secondary.xyz');
const secondaryStats = await this.computeStatsFor(secondaryNorm!.path, workDir, 'sample_secondary.xyz');
const topSpacing = primaryStats.stats.median_nn;
const bottomSpacing = secondaryStats.stats.median_nn;
const matchedSpacing = Math.max(topSpacing, bottomSpacing);
@ -141,6 +144,16 @@ export class PipelineService {
const rolesForThisRun: CloudRole[] = isCompare ? ['top', 'bottom'] : ['unique'];
// Portee de remplissage des trous FIXE sur tous les niveaux (derivee de la grille du niveau 0),
// et non recalculee a partir du gridStep de CHAQUE niveau. Avec un multiplicateur applique au
// gridStep de chaque niveau, la portee double a chaque niveau (x32 au niveau 5) et finit par
// combler des concavites reelles du contour de l'amas au lieu des seuls trous de scan -
// biais systematique croissant avec la decimation, voir CONCEPT.md §7.3. Le niveau 0 a le meme
// gridStep (= initialStep * gridStepMultiplier) en mode const_height (spatialStep=null) et en
// mode cloud_compare (spatialStep = initialStep * factor^0 = initialStep).
const gridStep0 = initialStep * gridStepMultiplier;
const maxEdgeLengthDefault = gridStep0 * config.maxEdgeLengthMultiplier;
for (let i = 0; i <= config.decimationSteps; i++) {
if (stopped) {
levels.push(this.skippedLevel(i, stopReason));
@ -151,7 +164,7 @@ export class PipelineService {
// Niveau 0 (cloud_compare) = pas apparie = initialStep (factor^0) ; niveaux suivants x2.
const spatialStep = !isCompare && i === 0 ? null : initialStep * Math.pow(config.decimationFactor, i);
const gridStep = (spatialStep ?? initialStep) * gridStepMultiplier;
const maxEdgeLength = params.maxEdgeLengthOverride ?? gridStep * config.maxEdgeLengthMultiplier;
const maxEdgeLength = params.maxEdgeLengthOverride ?? maxEdgeLengthDefault;
onProgress(
`Niveau ${i} : ${spatialStep === null ? 'nuage complet' : `decimation (pas=${spatialStep.toFixed(5)}m)`}, ` +
@ -161,7 +174,7 @@ export class PipelineService {
try {
const cloudResults: CloudPrepResult[] = [];
for (const role of rolesForThisRun) {
const sourcePath = i === 0 ? (role === 'bottom' ? secondaryNorm! : primaryNorm) : prevRaw[role]!;
const sourcePath = i === 0 ? (role === 'bottom' ? secondaryNorm!.path : primaryNorm.path) : prevRaw[role]!;
const isFirstConversion = !isCompare && i === 0;
const res = await this.prepareCloudAtLevel({
level: i,
@ -253,8 +266,11 @@ export class PipelineService {
await fs.rename(reportFile, path.join(rawDir, `L${i}_volume_report.txt`));
// AUTO_SAVE ON sauvegarde aussi une grille "*_HEIGHT_DIFFERENCE_*.bin" qu'on ne veut pas
// garder (les nuages combles sont maintenant en ASCII, donc tout .bin frais ici est ce fichier)
const strayGrid = await this.cc.findLatestFile(workDir, { suffix: '.bin', after: beforeVolume });
// garder. Filtrage par prefixe (pas juste ".bin") : sur un run avec peu de points, les
// L{n}_*.bin bruts fraichement ecrits (chaines depuis prepareCloudAtLevel) tombent aussi
// dans la fenetre "after" de findLatestFile - un filtre ".bin" nu risquerait de supprimer
// le nuage brut du niveau au lieu de la grille parasite, cassant le chainage de decimation.
const strayGrid = await this.cc.findLatestFile(workDir, { prefix: '', suffix: '.bin', after: beforeVolume, mustContain: '_HEIGHT_DIFFERENCE_' });
if (strayGrid) {
await fs.unlink(strayGrid).catch(() => undefined);
}
@ -381,6 +397,10 @@ export class PipelineService {
medianNn: primaryStats.stats.median_nn,
bboxApprox: primaryStats.stats.bbox,
},
coordinateOffset: {
primary: primaryNorm.offsetX !== null ? { x: primaryNorm.offsetX, y: primaryNorm.offsetY! } : null,
secondary: secondaryNorm?.offsetX !== null && secondaryNorm?.offsetX !== undefined ? { x: secondaryNorm.offsetX, y: secondaryNorm.offsetY! } : null,
},
},
levels,
mergedBinFile,
@ -399,14 +419,23 @@ export class PipelineService {
return report;
}
/** Convertit LAS/LAZ/COPC en ASCII (le paquet apt CloudCompare n'a pas le plugin LAS) ; laisse .bin tel quel. */
private async normalizeInput(inputPath: string, workDir: string, tag: string, onProgress: (msg: string) => void): Promise<string> {
/**
* Convertit LAS/LAZ/COPC en ASCII (le paquet apt CloudCompare n'a pas le plugin LAS) ; laisse .bin
* tel quel. Applique un decalage global X/Y pour preserver la precision float32 de CloudCompare
* en coordonnees projetees (voir CONCEPT.md §7.2).
*/
private async normalizeInput(
inputPath: string,
workDir: string,
tag: string,
onProgress: (msg: string) => void,
): Promise<{ path: string; offsetX: number | null; offsetY: number | null }> {
const lower = inputPath.toLowerCase();
if (!lower.endsWith('.las') && !lower.endsWith('.laz')) return inputPath;
if (!lower.endsWith('.las') && !lower.endsWith('.laz')) return { path: inputPath, offsetX: null, offsetY: null };
onProgress(`Conversion du nuage (${tag}) LAS/LAZ/COPC en ASCII...`);
const xyzPath = path.join(workDir, `input_${tag}.xyz`);
await this.py.convertLasToXyz(inputPath, xyzPath);
return xyzPath;
const { offsetX, offsetY } = await this.py.convertLasToXyz(inputPath, xyzPath);
return { path: xyzPath, offsetX, offsetY };
}
private async computeStatsFor(normalizedPath: string, workDir: string, sampleFileName: string): Promise<{ totalPointCount: number; stats: CloudStats }> {
@ -420,7 +449,22 @@ export class PipelineService {
}
const samplePath = path.join(workDir, sampleFileName);
await this.assertExists(samplePath, 'echantillonnage du nuage', res.stdout);
const stats = await this.py.computeStats(samplePath);
const rawStats = await this.py.computeStats(samplePath);
// Correction du biais de sous-echantillonnage (voir CONCEPT.md §7.1) : la distance mediane au
// plus proche voisin depend de la densite de points (median_nn ~ 1/sqrt(densite)). Un
// sous-echantillon de taille fixe (statsSampleCap) a une densite plus faible que le nuage
// complet des que ce dernier depasse statsSampleCap points, ce qui SURESTIME systematiquement
// median_nn/mean_nn d'un facteur sqrt(totalPointCount/sampleCount). Sans cette correction, un
// nuage LiDAR HD de 20M points (echantillon 50k, ratio 400x) voit son pas spatial estime
// gonfle d'un facteur ~20 - toute la grille de calcul et le remplissage de trous en heritent.
// No-op si le nuage complet est deja <= statsSampleCap (echantillon = nuage entier).
const correctionFactor = rawStats.count > 0 ? Math.sqrt(totalPointCount / rawStats.count) : 1;
const stats: CloudStats = {
...rawStats,
mean_nn: rawStats.mean_nn * correctionFactor,
median_nn: rawStats.median_nn * correctionFactor,
std_nn: rawStats.std_nn * correctionFactor,
};
return { totalPointCount, stats };
}

View file

@ -45,6 +45,7 @@ describe('reportToCsv', () => {
steps: 1,
densityMatch: null,
sampleStats: { sampleCount: 100, meanNn: 0.02, medianNn: 0.02, bboxApprox: { xmin: 0, xmax: 1, ymin: 0, ymax: 1, zmin: 100, zmax: 101 } },
coordinateOffset: { primary: null, secondary: null },
},
levels: [
{

View file

@ -85,6 +85,12 @@ export interface RunReport {
medianNn: number;
bboxApprox: { xmin: number; xmax: number; ymin: number; ymax: number; zmin: number; zmax: number };
};
/**
* Decalage (m, arrondi) soustrait de X/Y a la conversion LAS->ASCII pour preserver la precision
* float32 de CloudCompare en coordonnees projetees (voir CONCEPT.md §7.2). null si le fichier
* d'entree n'etait pas LAS/LAZ (ex. .bin deja normalise) ou en l'absence de second nuage.
*/
coordinateOffset: { primary: { x: number; y: number } | null; secondary: { x: number; y: number } | null };
};
levels: LevelResult[];
mergedBinFile: string | null;

View file

@ -11,6 +11,9 @@ export interface CloudStats {
median_nn: number;
std_nn: number;
bbox: { xmin: number; xmax: number; ymin: number; ymax: number; zmin: number; zmax: number };
/** 0.1e et 1er percentile de Z : estimateurs robustes du plan de reference, voir CONCEPT.md §7.7 */
z_p01: number;
z_p1: number;
approximate: boolean;
}
@ -83,7 +86,14 @@ export class PyHelperService {
* On convertit donc les entrees LAS/LAZ/COPC en ASCII XYZ via laspy avant de les transmettre
* a CloudCompare, qui lit nativement le format ASCII sans plugin.
*/
async convertLasToXyz(inputPath: string, outputXyzPath: string): Promise<{ count: number }> {
/**
* offsetX/offsetY : decalage (en metres, arrondi au metre) soustrait de X/Y avant ecriture, voir
* CONCEPT.md §7.2. Necessaire en coordonnees projetees (ex. Lambert-93, X/Y ~ 10^6-10^7 m) ou la
* precision float32 interne de CloudCompare degraderait sinon la position des points a l'echelle
* du decimetre. Le volume 2.5D etant invariant par translation, ce decalage n'affecte pas le
* resultat - seul zref (mode const_height) reste en coordonnees Z non decalees.
*/
async convertLasToXyz(inputPath: string, outputXyzPath: string): Promise<{ count: number; offsetX: number; offsetY: number }> {
return this.runJson('las_to_xyz.py', ['--input', inputPath, '--output', outputXyzPath]);
}
}

View file

@ -24,6 +24,16 @@ import { RunsService } from './runs.service';
const ALLOWED_EXT = ['.las', '.laz', '.copc.laz', '.bin'];
// Format uuidv4 genere par RunsService.createRun - toute route parametree par ":id" doit valider
// contre ce format AVANT de construire un chemin filesystem, sinon un id du type "..%2F..%2Fetc"
// (decode par Express dans les params de route) permet une traversee de repertoire post-auth
// (ex: telechargement/zip d'un dossier arbitraire du conteneur).
const RUN_ID_RE = /^[0-9a-f]{8}-[0-9a-f]{4}-[0-9a-f]{4}-[0-9a-f]{4}-[0-9a-f]{12}$/i;
function assertValidRunId(id: string): void {
if (!RUN_ID_RE.test(id)) throw new BadRequestException('Identifiant de run invalide');
}
function isAllowedFilename(name: string): boolean {
const lower = name.toLowerCase();
return ALLOWED_EXT.some((ext) => lower.endsWith(ext));
@ -36,7 +46,10 @@ const uploadStorage = diskStorage({
cb(null, dir);
},
filename: (_req, file, cb) => {
cb(null, `${Date.now()}_${Math.random().toString(36).slice(2)}__${file.originalname}`);
// path.basename retire tout composant de repertoire (../, chemin absolu) d'un originalname
// fourni par le client, avant que multer ne le concatene au dossier de destination.
const safeName = path.basename(file.originalname);
cb(null, `${Date.now()}_${Math.random().toString(36).slice(2)}__${safeName}`);
},
});
@ -122,16 +135,19 @@ export class RunsController {
@Get(':id')
detail(@Param('id') id: string) {
assertValidRunId(id);
return this.runs.getDetail(id);
}
@Delete(':id')
async remove(@Param('id') id: string) {
assertValidRunId(id);
return this.runs.delete(id);
}
@Get(':id/download/bin')
async downloadBin(@Param('id') id: string, @Res() res: Response) {
assertValidRunId(id);
const row = this.runs.getRowOrThrow(id);
const report = row.report_json ? JSON.parse(row.report_json) : null;
if (!report?.mergedBinFile) throw new NotFoundException('Fichier .bin fusionne non disponible pour ce run');
@ -141,6 +157,7 @@ export class RunsController {
@Get(':id/download/csv')
async downloadCsv(@Param('id') id: string, @Res() res: Response) {
assertValidRunId(id);
const filePath = path.join(this.runs.workDirFor(id), 'report.csv');
if (!fs.existsSync(filePath)) throw new NotFoundException('report.csv non disponible pour ce run');
res.download(filePath, `report_${id}.csv`);
@ -148,6 +165,7 @@ export class RunsController {
@Get(':id/download/zip')
async downloadZip(@Param('id') id: string, @Res() res: Response) {
assertValidRunId(id);
const dir = this.runs.workDirFor(id);
if (!fs.existsSync(dir)) throw new NotFoundException('Resultats non disponibles pour ce run');
@ -165,6 +183,7 @@ export class RunsController {
@Get(':id/images/:filename')
async image(@Param('id') id: string, @Param('filename') filename: string, @Res() res: Response) {
assertValidRunId(id);
const safeName = path.basename(filename);
if (safeName !== filename) throw new BadRequestException('Nom de fichier invalide');
const filePath = path.join(this.runs.workDirFor(id), 'images', safeName);