denoise : un nouvel algorithme/programme pour éliminer le bruit des MNT

sam 26-11-2011 Martin Laloux

SIG OpenSource


	Niveau	Intermédiaire
	Logiciels utilisés	denoise GDAL GRASS GIS Quantum GIS ImageJ
	Plateforme	Windows \| Mac \| Linux \| FreeBSD

Enlever ou du moins, réduire le bruit aléatoire d'un MNT, est crucial pour l'interprétation géomorphologique. Ce bruit introduit en effet de nombreuses erreurs dans les traitements et interprétations.

Les méthodes, les filtres pour le faire, sont nombreux et plus ou moins complexes (filtres médians, passe-haut, coupe-bas, par moyenne mobile, gaussien, par transformée d'ondelette ou géostatistiques). Ils ont fait l'objet de nombreux travaux et thèses et il existe un grand nombre d'entre eux publiés. Le but est toujours de tenter de préserver les structures essentielles du relief (crêtes, etc. ) en modérant le lissage.

Je ne vais pas rentrer ici dans le détail, mais me contenter de présenter un nouvel algorithme / programme beaucoup plus simple pour le faire : denoise

Il a été présenté en 2007 par Sun, X., Rosin, P.L., Martin, R.R. et Langbein, F.C. dans "Fast and Effective Feature-Preserving Mesh Denoising", ralph.cs.cf.ac.uk/papers/Geometry/MeshDenoising.pdf) pour supprimer le bruit aléatoire sur des maillages (mesh) de surface en trois dimensions tels ceux utilisés en CAD et autres tout en préservant l'information réelle du modèle comme les arêtes et les coins, à l'inverse du lissage.

Ils ont démontré l'efficacité de leur algorithme sur des modèles 3D divers où un bruit aléatoire avait été introduit, mais ces modèles sont relativement simples.

----->

avant après

image tirée de www.cs.cf.ac.uk/meshfiltering/index_files/Page342.htm

Ce n'est pas le cas avec les surfaces topographiques naturelles qui sont beaucoup plus complexes. Stevenson, J.A., Sun, X. et Mitchell, N.C., "Despeckling SRTM and other topographic data with a denoising algorithm" dans la revue Geomorphology (Volume 114, Issue 3, 15, 2010, p. 238-252, www.sciencedirect.com/science/article/pii/S0169555X09002967 ) ont publié une étude récente sur l'application de cet algorithme aux DEMs SRTM et aux données LIDAR.

Mitchell propose aussi sur sa page personnelle personalpages.manchester.ac.uk/staff/neil.mitchell/mdenoise/ des explications détaillées et le module r.denoise (grass.osgeo.org/wiki/GRASS_AddOns#r.denoise) a été placé sur le site des AddOns de GRASS GIS.

Principes de l'algorithme

Toutes les grilles de point sont préalablement transformées en grilles triangulaires puis l'algorithme fonctionne de manière itérative en modifiant à chaque étape les positions des sommets des facettes triangulaires qui forment le maillage (noeuds) jusqu'à arriver à une solution. L'ensemble est contrôlé par une valeur seuil.

a) à chaque itération, la meilleure orientation du vecteur normal à chaque facette (flèche rouge sur la figure suivante) est calculée à partir de la moyenne pondérée des orientations des vecteurs normaux des facettes adjacentes.
b) les noeuds sont ensuite éventuellement déplacés en fonction de cette meilleure orientation.

Le nombre d'itérations est défini par le paramètre n

La préservation des caractéristiques est obtenue en pondérant les normales des facettes adjacentes à chaque itération. Ceci se fait avec le paramètre t, qui correspond au cosinus de l'angle maximum autorisé entre les normales de deux facettes adjacentes.

si les orientations de deux facettes adjacentes sont supérieures au seuil autorisé (t) elles sont écartées du traitement (poids 0, points blancs dans la figure suivante)
dans les autres cas, elles sont plus ou moins pondérées suivant leurs similitudes et les noeuds sont éventuellement déplacés (points rouges sur la figure suivante).

figure adaptée de la Fig.1 de Stevenson et al. (2010)

Ce paramètre t affecte la préservation des structures :

si le paramètre t est élevé (angle de différence faible), le bruit sera préservé.
à l'inverse, un paramètre bas (angle de différence élevé) produira un lissage général, ce qui n'est pas le but de l'algorithme. Ainsi, à titre d'exemple, avec un t = 0.71, seules les ruptures de pentes supérieures à 45° seront préservées.

Tout le problème est alors de trouver les bons paramètres pour un cas donné :

c'est un des buts du travail de Stevenson et al. (2010). Après la justification de l'algorithme et la comparaison avec d'autres méthodes de « débruitage », ils cherchent à trouver, théoriquement et pratiquement, les « meilleurs » paramètres pour traiter les DEM SRTM. (le problème des données SRTM est qu'elles sont préalablement filtrées avant d'être soumises et que les bruits sont spatialement corrélés (voir "Accuracy and resolution of shuttle radar topography mission data", Smith, B., Sandwell, D., 2003, Geophysical Research Letters 30 (9), 1467).
De même, A. Muriy, sur la liste de diffusion de GRASS GIS (osgeo-org.1803224.n2.nabble.com/ASTER-GDEM-Version-2-gt-filtering-and-smoothing-td6924936.html), propose déjà les « meilleurs » paramètres n et t pour les nouveaux DEM Aster (V2) (www.portailsig.org/content/sortie-du-gdem-aster-v2-modele-numerique-d-elevation)

En pratique : téléchargement

denoise

www.cs.cf.ac.uk/meshfiltering/index_files/Page342.htm propose le téléchargement du programme denoise sous forme d'exe pour Windows et de code source à compiler pour Linux et Mac OS X (très facile). Le programme est disponible gratuitement.

r.denoise, module pour GRASS GIS

r.denoise est disponible sur trac.osgeo.org/grass/browser/grass-addons/raster/r.denoise

En pratique : exécution du programme

Le programme s'exécute à partir de la ligne de commande. Le « débruitage » est contrôlé par les deux paramètres, n, le nombre d'itérations et t, le seuil. L'ensemble est bien expliqué par Mitchel sur sa page personnelle personalpages.manchester.ac.uk/staff/neil.mitchell/mdenoise/ (notamment l'utilisation de GDAL pour la transformation de fichiers en grilles ESRI (.asc) et la création des fichiers de projection (.prj) nécessaires au traitement).

fichiers asc

Les fichiers .asc doivent être projetés dans un système de projection cartésien. Ils doivent donc nécessairement être accompagnés de fichiers de projection .prj

MDenoise(.exe) -i monMNT.asc -n 5 -t 0.99 -o monnouvMNT.asc

La commande donnera pour résultat les fichiers monnouvMNT.asc et monnouvMNT.prj

fichiers xyz

Les points doivent être encodés sous la forme de fichiers ASCII se terminant par .xyz. Ils doivent contenir les valeurs x, y et z séparées par des espaces dans les trois premières colonnes. Le reste est ignoré (éventuelles colonnes supplémentaires)

MDenoise(.exe) -i mespoints.xyz -n 5 -t 0.99 -o mesnouvpoints.xyz

Durant le processus, les points sont susceptibles d'être déplacés dans les trois dimensions. Si on ne veut travailler que sur l'altitude, il suffit d'ajouter le drapeau '-z' à la commande.

En pratique : application au nouveau GDEM Aster V2 ("ASTER GDEM is a product of METI and NASA")

Les essais préalablement cités utilisent des zones à reliefs fort contrastés (montagnes avec de fortes pentes et plaines) mais qu'en est-il des autres régions. Le traitement a été effectué sur la région de Pondrôme dans le sud de la Belgique avec les GDEM Aster V2 (www.portailsig.org/content/sortie-du-gdem-aster-v2-modele-numerique-d-elevation). Les paramètres donnés par A. Muriy sont t = 0.8 et n = 10-20. Nous présentons ici des résultats exemplatifs sans justification théorique ou autre et sans essai d'interprétation, simplement pour illustrer des résultats de l'application de l'algorithme.

DEM traité avec GRASS GIS à partir des données ASTER-GDEM version 2 (ici avec une altitude x 2 pour améliorer la visibilité du relief) :

Carte des pentes (en degrés) réalisée avec GRASS GIS pour illustrer les contrastes du relief:

Variation du nombre d'itérations avec une valeur t constante

Résultats du traitement avec t = 0,8 et n = 5 - 10 - 15 - 20 - 25 -30 (GIF animé)

On remarque que le résultat du traitement se stabilise rapidement après un certain nombre d'itérations et ne bouge plus par après.

Variation de la valeur t avec un même nombre d'itérations

Résultats du traitement avec n = 20 et t = 0,9 - 0,8 - 0,5 - 0,1 (GIF animé)

ici aussi, le résultat se stabilise rapidement.

Examen éventuel des modifications

Les modifications apportées peuvent être observées en calculant les différences entre le DEM original et un des DEMs modifiés (ici, celui calculé avec les paramètres n = 10 et t = 0.95 (GRASS GIS : r.mapcalc : DEM original - DEM résultant, Quantum GIS avec la calculatrice raster)

Dem original Dem résultant n = 10 et t = 0.95

Résultat :

mais la manière la plus parlante pour visualiser les différences est, à mon avis, d'utiliser ImageJ qui permet de visualiser (et de traiter) les fichiers .asc :

Dem original (avec Look Up table 16 colours) Dem résultant n = 10 et t = 0.95

Détail :

Dem original Dem résultant n = 10 et t = 0.95

Conclusions

Les conclusions de l'article de Stevenson et al. (2010) montrent que cet algorithme constitue une très bonne alternative aux méthodes traditionnelles. D'après eux, il donne même de meilleurs résultats du fait de ses caractéristiques, pondérations et itérations convergentes (stabilisation rapide du résultat).

Je vous laisse juge, mais il est surtout librement téléchargeable et très facile d'utilisation sur toutes les plateformes, sans nécessité de logiciel spécialisé (Stevenson et al., 2010, ont utilisé GRASS GIS pour tous leurs traitements complémentaires, mais ce peut être fait avec Qgis et autres logiciels SIG). Pour moi, l'élément perturbant est qu'il faut « jouer » avec les paramètres pour obtenir les résultats adéquats, sans cadre théorique. La rapide convergence des résultats est le seul critère.

Tous les traitements ont été effectués sur Mac OS X avec denoise, GRASS GIS 6.4.1 et r.denoise, ImageMagick pour les GIFs animés et ImageJ version 1.45s avec le plugin ESRIASCIIGrid (terraincartography.com/imagej/imagej.html) pour traiter les fichiers .asc