Low cost DEM acquisition from UAV and signal processing using Micmac

Summary of the information gathered on the micmac forum from knowledgeable users answering my questions about UAV acquired image processing for DEM generation.
  1. Acquire data. In our case, a DJI Phantom3 Pro was used, with each image geotagged by the C/A L1-only GPS onboard receiver. No GPS-postprocessing was performed. Flight elevation was 90 m above takeoff altitude, with an horizontal speed set to maximum (about 10 m/s) and one picture manually acquired every 1 s (more or less). Since the flight plan was not implemented at the time (Sept. 2015), all flights were performed manually with an attempt to overlap successive passes from visual features seen in the real time transmission to the pilot. The inset map was also used to try and match successive flights over nearby regions. Each flight duration was about 20 min, so the whole dataset of about 1000 useful images requires a storage space of 4-5 GB.
  2. All images related to a given flight are collected in a single directory on a GNU/Linux computer. We first extract GPS coordinates (WGS84 spherical coordinates) from each picture and converted, using QGis, to the local projected coordinate system WGS84/UTM33N. I enjoy manual processing of headers using shell script that look something like (assuming the copied files have an extension .JPG in capital letters): 
    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 56-100 | cut -d\  -f1,3,4 | sed "s/'//g" | sed 's/"//g' > latitude.txt
cat longitude | cut -c 56-100 | cut -d\  -f1,3,4 | sed "s/'//g" | sed 's/"//g' > longitude.txt
cat altitude | cut -c 56- | cut -d\  -f1 > altitude.txt
paste noms latitude.txt longitude.txt altitude.txt
    I sometimes have to tune the 
    56-100
    argument of cut depending on filename length and exiftool organization. The value (56) must be the index of the number value of interest (latitude, longitude or altitude).
    Because WGS84 coordinates are in degree.second format and I want decimal degrees, I convert using 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
    and finally the resulting dataset needed for QGis is generated, after removing the GNU/Octave header from the .dec files, by combining all these columns in a single file 
    paste  latitude.dec longitude.dec altitude.txt noms > for_qgis.txt
    In QGis, I load the files using the comma-shaped icon, with latitude and longitude assigned respectively to the Y and X columns. Right click on the item in the item-list, Save As and select CRS WGS84-UTM33N (EPSG:32633). Save as positions_31N.csv and clean the resulting file from unwanted columns and comma separator with 
    cat positions_33N.csv | cut -d, -f1-2,5,6 | sed 's/,/ /g' | sed 's/^43//g' | sed 's/ 875/ /g' > positions_33N.final
    Fields 1-2,5,6 are longitude, latitude, altitude, filename, so I manually (vi) add the header 
    #F=X Y Z N
    In my particular case, WGS84-UTM33N latitudes are around 8750000 and longitudes are around 430000 so I remove the leading 875 and 43 to match the region I work in and help Micmac's computation by preventing excessively large coordinates from being used. THIS MUST BE ADAPTED TO EACH REGION. I my case, some remaining coordinates might be left with 876xxxx so I also replace (sed) all 
    \ 876
    with 
     1
    .
  3. Finally the MicMac part: 
    mm3d OriConvert OriTxtInFile positions_33N.final jmfgps MTD1=1 NameCple=FileImagesNeighbour.xml CalcV=1 ImC=DJI_0115.JPG NbImC=25
    generates a FileImagesNeighbour.xml image pairs, and a PATC list which must be adapted to shell scripting (remove spaces around = and surround the list with double quotes).
    I have no reason of trusting the image matching algorithm, so I want to assess how appropriate the use of the picture GPS coordinates for identifying image pairs is. For that I wrote a couple of scripts for converting the image name pairs to coordinate pairs and encapsulate the coordinate set into a WKT file for QGis. This is done first by 
    cat FileImagesNeighbour.xml | sed 's///g' | sed 's/<\/Cple>//g' | grep -v S | grep -v x | sed 's/     //g' > imgpairs
    I generate a script to be executed later by executing an awk script (named script.awk) containing 
    { print "sed 's/"$4"/43"$1" "875$2"/g' | \\" }
    where p*final includes the coordinates we created earlier (it was called positions_33N.final): 
    cat p*final | awk -f script.awk > positions_31N.sh
    In the script positions_31N.sh, I replace the first line with cat imgpairs | \ and I remove in the last line the trailing | \. I finally execute this script to generate the wkt file to be loaded in QGis (using the coma-shaped icon): 
    sh ./positions_31N.sh | sed 's/^/LINESTRING(/g' | sed 's/$/)/g' | sed 's/\ /,/2' > aida.wkt
    The result is for example the following image (for a couple of datasets):
    PATC="DJI_0515.JPG|DJI_0514.JPG|DJI_0516.JPG|DJI_0172.JPG|DJI_0173.JPG|DJI_0517.JPG|DJI_0513.JPG|DJI_0171.JPG|DJI_0174.JPG|DJI_0170.JPG|DJI_0175.JPG|DJI_0518.JPG|DJI_0512.JPG|DJI_0169.JPG|DJI_0176.JPG|DJI_0428.JPG|DJI_0429.JPG|DJI_0168.JPG|DJI_0177.JPG|DJI_0511.JPG|DJI_0427.JPG|DJI_0430.JPG|DJI_0496.JPG|DJI_0495.JPG|DJI_0167.JPG"
P=".*JPG"
    PATC might or might not be relevant, depending on the convergence of the camera model. Sometimes the camera lens model (Tapas) will not converge (residue = NaN) and I have to manually select a subset of the pictures better suited for lens characteristics characterization. I iterate on 20-image sequences until the RadiaStd model converges, usually with some success after the 3rd attempt at most ?! 
    mm3d Tapioca File "FileImagesNeighbour.xml" -1
mm3d Tapas RadialStd "$PATC" Out=Cal
mm3d Tapas AutoCal "$P" InCal=Cal Out=Init
CenterBascule "$P" Init jmfgps tmp 
mm3d AperiCloud "$P" tmp
    Here we check on the coarse point cloud that all camera positions are correct and that the coordinate system is coherent. If elevations are much more than 90 m or camera spacing much more than 10 m, then something was wrong in the coordinate extraction. If one camera is not in the UAV path, I throw away the associated picture. 
    mm3d Malt Ortho "$P" tmp "DirMEC=Resultats" UseTA=1 ZoomF=2 ZoomI=32 Purge=true AffineLast=false
mm3d Tawny Ortho-Resultats/ RadiomEgal=0
Nuage2Ply Resultats/NuageImProf_STD-MALT_Etape_7.xml Attr="Ortho-Resultats/Ortho-Eg-Test-Redr.tif"
    results in the orthorectified images in Ortho-Resultats and the DEMs in Resultats.
    4. From now on we have the fun part: inserting the orthorectified image and DEM in QGis.
    1. for the orthorectified image, we are only interested in 
      Ortho-Resultats/*edr*
      : the TIF format is surely nice but huge, so I first convert all useful images to PNG 
      for i in Orth*/*Redr*.tif; do nom=`echo $i|sed 's/tif/png/g'`;convert $i $nom;done
      and combine the resulting image in a single huge PNG orthorectified image since I do not know how to trivially compute the corner coordinates of the sub-frames (should be easy to compute ... homework). 
      export MAGICK_TMPDIR=/home/jmfriedt/t
montage -tile 1x2 -geometry +0+0 Ortho-Eg-Test-Redr_Tile_0_0.png Ortho-Eg-Test-Redr_Tile_0_1.png  out.png
      The TMPDIR is due to the fact that my 512 MB tmpfs /tmp directory is not large enough for the imagemagick processing. The tile order seems rationale from the imagemagick web page, but I still need a few trial and errors every time I combine the images to get them in the right position. This is a lengthy computation: in my case, the 210 JPG images acquired from the UAV generate a 14300x40500 pixel large, 251 MB PNG image.
      The associated world file out.pngw 
      0.045
0
0
-0.045
8839.62
10127.52
      is appended with the offset coordinate of the corner we initially removed to help micmac's computation, resulting in 
      0.045
0
0
-0.045
438839.62
8760127.52
    2. We validate the orthoimage is properly located over a Formosat background reference image:

      Left to right: raw raster orthophoto loaded with undefined areas in black -- undefined areas removed by defining the (0, 0, 0) color as transparent -- zoom on the glacier front and outgoing stream.
      To remove the black background by converting it to a transparent layer, right click on the layer name in QGis, Properties - Transparency and "+" icon to add RGB=0 0 0 to the transparency list.
  4. Finally, the DEM is processed quite the same way. The corner coordinate in the tfw image associated with Resultats/Z_Num7_DeZoom2_STD-MALT.tif is also appended with the position offset as was done above, and the raster file is loaded in QGis. However the 16-bit TIF file does not encode an altitude: a conversion must be made using the information provided in the (altitude offset) and (DEM altitude resolution) of the Resultats/Z_Num7_DeZoom2_STD-MALT.xml. The actual DEM elevation is TIF_value*ResolutionAlti+OrigineAlti as explained at this web page. For this computation, I use in QGis the Raster - Raster Calculator with the Raster band the DEM I loaded as a geotif file, Click on Current Layer Extent, and generate the Raster Calculator Expression such as 
    "Z_Num7_DeZoom2_STD-MALT@1"*0.09-9.36
    (since in my example the ResolutionAlti is 9 cm and OrigineAlti is -9.36 m).

    Left to right: raw formosat background image (taken in 2009), greyscale DEM overlay, colored DEM overlay

    Left to right: cropped colored DEM, and raw background. Notice the hill on the top-right part of the DEM.
Having acquired images of the same area over multiple days (2-day separation between dataset acquisition), I wonder how reproducible the measurement is. The 2015/09/24 DEM is considered as reference, and I am moving the 2015/09/26 reference to visually match as best some if the features, namely a canyon carved during the summer by a stream flowing in the newly de-iced area by glacier retreat. The manual translation (Rasmover plugin from QGis) is about 6 meters to match some GCP manually placed on some curves of the stream on the 09/24 DEM. Having moved the 2015/09/26 DEM, the result is saved in a geotif file (somehow processing the moved DEM always results in NaN if the dataset is not saved first), I add 60 m to the 2015/09/26 DEM to bring is close to 0. This offset might be due to the fact that DEMs are referenced to take off altitude, and both DEMs starting point were not at the same altitude. The following cross section (Profile plugin) illustrates how reproducible the elevation measurement is. The standard deviation over the whole length on the elevation difference is 80 cm, which drops to 50 cm if we only consider the first 200 m of the profile and exclude the last 50 m where there is an obvious position mismatch.

Left to right: elevations along profile (screenshot) ; elevation along profile and elevation differences ; aerial picture of the canyon being investigated.
Campari, as explained at http://www.wiki.gnewarchaeology.it/w/MicMac, is used "for compensation of heterogeneous measures, that is: tie points and ground control points". My understanding: camapari will use the camera position extracted from photogrammetry and the (noisy/biased) GPS coordinates to correct the DEM by combining, in a weighted average, both information (had I had GCPs, they could have been thrown in as well).
The difference in DEM with and without Campari, used through the command 
Campari "$P" tmp tmp-Campari EmGPS=[jmfgps,2] FocFree=1 PPFree=1 
mm3d AperiCloud "$P" tmp-Campari
mm3d Malt Ortho "$P" tmp-Campari "DirMEC=Resultats-Campari" UseTA=1 ZoomF=2 ZoomI=32 Purge=true AffineLast=false
mm3d Tawny Ortho-Resultats-Campari/ RadiomEgal=0
Nuage2Ply Resultats-Campari/NuageImProf_STD-MALT_Etape_7.xml Attr="Ortho-Resultats-Campari/Ortho-Eg-Test-Redr.tif"
(ie executing Campari to generate a new Orientation directory called tmp-Campari, and from there running all the last processing sequences on this new orientation) yields the following chart, with differences in the +/- 30 cm elevation range.