FOSS4G-fr tutorial "La photogrammétrie pour tous : MicMac pour la reconstruction 3D de scènes géoréférencées à partir de photographies numériques."

Résumé court :

Un outil de traitement photogrammétrique pour la reconstruction de nuages de points 3D d'une scène observée de divers points de vues, applicable aux modèles numériques d'élévation ou bâtiments.

Résumé long :

La photogrammétrie connait un nouvel essor avec la prolifération des outils de prises de vues géoréférencées (eg. téléphone mobile). L'analyse d'une multitude de photographies visant un meme objet, bâtiment ou scène permet de construire un nuage de points représentant la scène considérée en 3D. Ce nuage de point s'insère dans son contexte géographique dans un outil de SIG tel que QGis pour analyse. La présentation se portera sur MicMac, un outil opensource de traitement photogrammétrique permettant la génération de nuages de points géoréferencés et colorisés. Deux exemples pratiques sont proposés : un cas de reconstruction de scène en intérieur, et un cas appliqué à un bâtiment en extérieur permettant d'exploiter le géoréférencement des photographies pour passer d'un référentiel arbitraire au référentiel géographique. L'application aux photographies aériennes prises par drone sera illustré sur des clichés acquis par les présentateurs.

Programme

  1. Pésentation de MicMac : généralités sur la photogrammétrie
  2. Conditions de prises de vues, sur les objectifs, sur le cadrage et sur le réglage de l'appareil photo
  3. Exemple de traitement sur une séquence de prises de vues en intérieur (fisheye) : prise en main de MicMac (interface)
    1. Identification des points homologues
    2. Identification des propriétés optiques de l'objectif et de la caméra
    3. Identification de la position de prise de vue des photographies
    4. Validation sur un nuage de points grossier
    5. Créatio du nuage de points dense
  4. Exemple de traitement sur des prises de vues géoréférencées en extérieur.
    1. En plus des cas précédents, passage du référentiel arbitraire au référentiel géographique (WGS84 -> WGS84UTM31N ou Lambert93 -> centrage de la photographie près de la coordonnée (0,0))
    2. insertion du nuage de points dans QGis
    3. colorisation du nuage de points
    4. modèle numérique d'élévation.

Au travail ...

Le jeu de données qui sera traité est disponible ICI. Ces données ont été acquises avec un téléphone Sony Xperia Z5 exécutant l'application OSMTracker (pour avoir la trace du chemin parcouru lors des prises de vues) et sa fonction Appareil Photo. Les photographies prises par ce téléphone sont automatiquement géoréférencées (tag GPS dans la sortie de exiftool) :

Compilation des outils :

MicMac à partir des sources :

git clone https://github.com/micmacIGN/micmac
cd micmac
mkdir build
cd build
cmake ../ -DBUILD_POISSON=1 -DBUILD_RNX2RTKP=1
cd UILD_RNX2RTKP make -j4 && make install

Autres

Vérifier le bon fonctionnement de exiftool et qgis, ainsi que la disponibilité de ImageMagick et GNU/Octave. Micmac nécessite par ailleurs exiv2.
quelques Goctets seront nécessaires sur le support de stockage non-volatil pour conserver les fichiers temporaires.

Version simplifiée

Nous fournissons une machine virtuelle (2 GB !) assemblée et testée contre l'archive (v5.0.16) de VirtualBox (sous Debian/GNU Linux Jessie) complétée des ExtraPacks (en administrateur, VBoxManage extpack install Oracle_VM_VirtualBox_Extension_Pack-5.0.16-105871.vbox-extpack).
La machine virtuelle permet de se connecter en root (pas de mot de passe) en mode console, ou login=passwd=foss4g en utilisateur. Pour lancer l'interface graphique, startx lance icewm (bouton de droite de la souris pour lancer un terminal). La visualisation des images est possible par mirage. Le répertoire /home/foss4g contiendra les fichiers temporaires, tandis que les images à traiter sont dans /home/data.

Tutorial :

Référentiel arbitraire :

mm3d Tapioca Line "2016-04-02_12-1[2-6]{1}....jpg" 1500 3
mm3d Tapas RadialStd "2016-04-02_12-12-11.jpg|2016-04-02_12-12-26.jpg|2016-04-02_12-12-44.jpg|2016-04-02_12-12-58.jpg|2016-04-02_12-13-11.jpg" Out=Cal
mm3d Tapas AutoCal "2016-04-02_12-1[2-6]{1}....jpg" InCal=Cal Out=All
mm3d AperiCloud "2016-04-02_12-1[2-6]{1}....jpg" All Out=PosCams_All.ply
mm3d SaisieMasq 2016-04-02_12-15-08.jpg   # bouton de gauche pour le masque, SHIFT-bouton gauche pour fermer, CTRL-bouton droite pour exit
mm3d Malt GeomImage "2016-04-02_12-1[2-6]{1}....jpg" All Master="2016-04-02_12-15-08.jpg" DirMEC=Arbitraire UseTA=1 ZoomF=4 ZoomI=32 Purge=true AffineLast=false # DefCor=0.01
mm3d Nuage2Ply Arbitraire/NuageImProf_STD-MALT_Etape_6.xml Attr=2016-04-02_12-15-08.jpg RatioAttrCarte=4

Positions des prises de vues :

En shell bash, recherche des coordonnées de chaque photographie à partir des entêtes EXIF : 
for i in *.jpg; do echo  $i ; done > noms
for i in *.jpg; do echo -n $i ; exiftool $i | grep atitu | grep -v R; done > latitude
for i in *.jpg; do echo -n $i ; exiftool $i | grep ngitu | grep -v R; done > longitude
for i in *.jpg; do echo -n $i ; exiftool $i | grep Altit | grep -v R; done > altitude
cat latitude  | cut -c 58-100 | cut -d\  -f1,3,4 | sed "s/'//g" | sed 's/"//g' > latitude.txt
cat longitude | cut -c 58-100 | cut -d\  -f1,3,4 | sed "s/'//g" | sed 's/"//g' > longitude.txt
cat altitude | cut -c 58- | cut -d\  -f1 > altitude.txt
Alternativement, le fichier GPX issu de OSMTracker contient les mêmes informations (balise wpt) : 
grep "/jpg/g' | tr '\n' ' ' | sed 's/jpg/jpg\n/g' | sed 's/">//g' | sed 's///g' | sed 's/<\/ele>//g' | expand | sed 's/  */ /g' | sed 's/^ //g'

Degrés décimaux :

Sous GNU/Octave : 
load latitude.txt
l=latitude(:,1)+latitude(:,2)/60+latitude(:,3)/3600;
save -text latitude.dec l
load longitude.txt
l=longitude(:,1)+longitude(:,2)/60+longitude(:,3)/3600;
save -text longitude.dec l

Coordonnées projetées :

Retirer les premières lignes des fichiers issus de GNU/Octave (*.dec), puis 
paste  latitude.dec longitude.dec altitude.txt noms > for_qgis.txt
On pourra se simplifier la vie en éditant le fichier for_qgis.txt et en ajoutant une ligne X Y Z N comme labels pour que QGis reconnaisse quelle colonne contient quelle information.
Dans QGis, charger le fichier ASCII (icone en forme de virgule), sélectionner WGS84 comme CRS, puis sauver en UTM31N (ou Lambert93) pour travailler dans un référentiel projeté. Une fois sorti de QGis, dans le shell nous ne conservons que les colonnes utiles et tronquons la valeur moyenne des colonnes pour évitre les erreurs de précision sur les calculs 
cat from_qgis.csv  | cut -d, -f1-2,5,6 | sed 's/,/ /g' | sed 's/^418//g' | sed 's/ 530/ /g' > positions_31N.final
Finalement, nous modifions l'entête pour le rendre compatible avec Micmac : 
#F=X Y Z N
Nous avons fini de créer le fichier de coordonnées GPS des prises de vues, qu'on suppose avoir nommé positions_31N.final. Micmac ne travaille pas sur un fichier ASCII brut mais sur un fichier au format XML, qui est généré par 
mm3d OriConvert OriTxtInFile positions_31N.final Nav-Brut-RTL MTD1=1 NameCple=FileImagesNeighbour.xml CalcV=1

Renseigner les propriétés du capteur optique :

La taille du capteur optique doit être renseignée dans include/XML_User/DicoCamera.xml, avec un identifiant identique à celui fourni par exiftool. Pour des téléphones Samsung S3 ou Sony Experia, nous utilisons <MMCameraDataBase> <CameraEntry> <Name> E5823 </Name> <SzCaptMm> 6.16 4.62 </SzCaptMm> <ShortName> E5823 </ShortName> </CameraEntry> <CameraEntry> <Name> GT-I9300 </Name> <SzCaptMm> 4.8 3.6 </SzCaptMm> <ShortName> GT-I9300 </ShortName> </CameraEntry> </MMCameraDataBase> dont la taille du capteur CMOS est plus ou moins exacte selon les informations glanées sur internet, et dont micmac semble parfaitement s'accomoder.

Trouver les points homologues sur les photographies :

La seule subtilité ici est que rien ne ressemble plus à une façade de cathédrale qu'une autre façade de cathédrale. Micmac n'arrive donc pas à positionner correctement sur le nuage dense de points des points qu'il croît appartenir au même mur mais qui sont sur des faces différentes. Pour palier à ce problème, nous pourrions limiter la recherche de points homologues entre images adjacentes, au lieu de rechercher sur l'ensemble de tous les couples possibles. Alternativement, nous exploitons les informations de position GPS issues du traitement précédent pour indiquer à micmac quels sont les couples les plus probables : 
mm3d Tapioca File FileImagesNeighbour.xml 1500
Noter la résolution de 1500, soit environ 1/3 de la résolution selon le plus grand axe de chaque image, pour ne pas s'accrocher sur des détails insignifiants. La validation de l'estimation de la position des prises de vues est un point clé des traitements qui vont suivre : le respect médiocre des règles de prises de vues imposées par micmac lors de l'acquisition des images impose de rechercher les paramètres de convergence des algorithmes d'identification des propriétés optiques de l'objectif et conditions de prises de vues.

Gauche : cas Mulscale dans lequel des caméras sont placées à des positions aberrantes. Droite : résultat correct dans le cas de File.

Créer le nuage de points grossier :

Nous trouvons les propriétés de la caméra supposée fixe au cours des prises de vues 
mm3d Tapas RadialStd ".*jpg" Out=Cal
Si l'étape précédente ne se fait que sur un sous ensemble des points, nous propageons le calcul de position de caméra pour toutes les prises de vues par 
mm3d Tapas AutoCal ".*jpg" InCal=Cal Out=All-Rel
qui n'est pas nécessaire ici (et qui au contraire induit une estimation erronée de la position des prises de vues dans notre cas).
Le passage du référentiel arbitraire de Mimcac au référentiel absolu (tenant compte des positions GPS des prises de vues) s'obtient par 
CenterBascule ".*jpg" Cal Nav-Brut-RTL tmp # CalcV=1
avant de finalement générer le nuage de points grossier contenant les positions de caméras 
mm3d AperiCloud ".*jpg" tmp Out=PosCams_GPS.ply

Créer le nuage de points dense :

Le nuage de point dense représentant la cathédrale est formée de plusieurs sous-nuages denses, chacun calculé autour d'une image maîtresse (Master=). De cette façon, toutes les façades sont couvertes, et pour la première, par deux photographies prises selon des directions différentes pour éviter trop de zones non-couvertes. 
mm3d Malt GeomImage "2016-04-02_12-09-53.jpg|2016-04-02_12-10-04.jpg|2016-04-02_12-10-16.jpg|2016-04-02_12-10-27.jpg|2016-04-02_12-10-49.jpg|2016-04-02_12-10-59.jpg|2016-04-02_12-11-13.jpg|2016-04-02_12-11-44.jpg|2016-04-02_12-11-56.jpg|2016-04-02_12-12-11.jpg|2016-04-02_12-12-26.jpg|2016-04-02_12-12-44.jpg|2016-04-02_12-12-58.jpg" tmp Master="2016-04-02_12-10-04.jpg" DirMEC=Result1GPS UseTA=1 ZoomF=2 ZoomI=32 Purge=true AffineLast=false # DefCor=0.01
mm3d Nuage2Ply Result1GPS/NuageImProf_STD-MALT_Etape_7.xml Attr=./2016-04-02_12-10-04.jpg RatioAttrCarte=2

mm3d Malt GeomImage "2016-04-02_12-12-44.jpg|2016-04-02_12-12-58.jpg|2016-04-02_12-13-11.jpg|2016-04-02_12-13-32.jpg|2016-04-02_12-13-44.jpg|2016-04-02_12-14-06.jpg|2016-04-02_12-14-26.jpg|2016-04-02_12-14-35.jpg|2016-04-02_12-14-54.jpg|2016-04-02_12-15-06.jpg|2016-04-02_12-15-08.jpg|2016-04-02_12-15-27.jpg|2016-04-02_12-15-40.jpg" tmp Master="2016-04-02_12-13-32.jpg" DirMEC=Resultats2GPS UseTA=1 ZoomF=2 ZoomI=32 Purge=true AffineLast=false
mm3d Nuage2Ply Resultats2GPS/NuageImProf_STD-MALT_Etape_7.xml Attr=./2016-04-02_12-13-32.jpg RatioAttrCarte=2

mm3d Malt GeomImage "2016-04-02_12-14-35.jpg|2016-04-02_12-14-54.jpg|2016-04-02_12-15-06.jpg|2016-04-02_12-15-08.jpg|2016-04-02_12-15-27.jpg|2016-04-02_12-15-40.jpg|2016-04-02_12-15-59.jpg|2016-04-02_12-16-04.jpg|2016-04-02_12-16-20.jpg|2016-04-02_12-16-36.jpg|2016-04-02_12-16-40.jpg" tmp Master="2016-04-02_12-16-20.jpg" DirMEC=Resultats3GPS UseTA=1 ZoomF=2 ZoomI=32 Purge=true AffineLast=false
mm3d Nuage2Ply Resultats3GPS/NuageImProf_STD-MALT_Etape_7.xml Attr=./2016-04-02_12-16-20.jpg RatioAttrCarte=2

mm3d Malt GeomImage "2016-04-02_12-09-53.jpg|2016-04-02_12-10-04.jpg|2016-04-02_12-10-16.jpg|2016-04-02_12-10-27.jpg|2016-04-02_12-10-49.jpg|2016-04-02_12-10-59.jpg|2016-04-02_12-11-13.jpg|2016-04-02_12-11-44.jpg|2016-04-02_12-11-56.jpg|2016-04-02_12-12-11.jpg|2016-04-02_12-12-26.jpg|2016-04-02_12-12-44.jpg|2016-04-02_12-12-58.jpg" tmp Master="2016-04-02_12-10-59.jpg" DirMEC=Result4GPS UseTA=1 ZoomF=2 ZoomI=32 Purge=true AffineLast=false # DefCor=0.01
mm3d Nuage2Ply Result4GPS/NuageImProf_STD-MALT_Etape_7.xml Attr=./2016-04-02_12-10-59.jpg RatioAttrCarte=2
M. Pierrot-Deseilligny nous informe d'une nouvelle fonctionnalité de MicMac, à savoir la capacité à automatiquement calculer les masques et assembler ces divers nuages de points, et ce en une seule commande : 
mm3d C3DC QuickMac ".*jpg" Init
avec QuickMac la recherche rapide de solution, qui se décline aussi en MicMac et BigMac (pour la plus haute résolution).

Inclusion dans QGis :

Les nuages de points sont insérés dans CloudCompare, fusionnés en un seul nuage de points global, et sauvé en format ASCII avec 6 colonnes/point, les coordonnées (X, Y, Z) et la couleur (R, G, B).
Ce nuage de points est translaté dans le référentiel absolu pour compenser la soustraction de valeur moyenne que nous avons retiré auparavant : 
cat nuage.txt | sed 's/^/418/g' | sed 's/ /530/1' > nuage_pos.txt
Ce nuage de points est inséré dans QGis (icone en forme de virgule), et superposé avec le fond de carte OpenStreeMap par Vector -> OpenStreetMap -> Dowload Data .... Le résultat est ci-dessous

Gauche : nuage de points inséré dans QGis sur fond OpenStreeMap. Droite : coloration des points (RGB).

Résultats :

Ressources :

  1. M. Pierrot-Deseilligny, MicMac , une solution libre de photogrammétrie, RMLL 2015
  2. J.-M Friedt, E. Bernard, F. Tolle, Utilisation de Micmac pour la génération de modèle numérique d'élévation par traitement d'images acquises par microdrone GNU/Linux Magazine 191 (March 2016), pp.48-57
  3. J.-M Friedt, Reconstruction de structures tridimensionnelles par photographies : le logiciel MicMac OpenSilicium 12 (Sept-Oct-Nov 2014)