UAV based DEM generation

A set of 173 pictures was collected above the Besancon obervatory using a DJI Phantom 3 Pro UAV. All pictures are georeferenced using the on-board GPS receiver, with elevation above ground (at take off) given in the altitude field.
First in bash script: extract all picture georeferencing information using exiftool by selecting the appropriate GPS fields. Notice that selecting the right field (47-100) requires tuning the initial value if subdirectories are considered.
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 47-100 | cut -d\  -f1,3,4 | sed "s/'//g" | sed 's/"//g' > latitude.txt
cat longitude | cut -c 47-100 | cut -d\  -f1,3,4 | sed "s/'//g" | sed 's/"//g' > longitude.txt
cat altitude | cut -c 47- | cut -d\  -f1 > altitude.txt
Latitude and longitudes are in degree.min.second format, so we convert to fractional degrees using GNU/Octave:
load latitude.txt 
save -text latitude.dec l
load longitude.txt 
save -text longitude.dec l
Each .dec file is prefixed with the GNU/Octave header which we delete using our favorite text editor (i.e. vi). Having done so, the columns are pasted in the single file in the unix shell:
paste  latitude.dec longitude.dec altitude.txt noms > pour_qgis.txt
This file is fed to QGIS for conversion from WGS84 spherical coordinates to projected WGS84/UTM31N coordinates (open the file using the colon icon of QGIS, select input CRS as WGS84, right click on the layer name and Save As, select WGS84/UTM31N as output CRS and save CSV file -- will be positions_31N.csv). We then remove the position offset to prevent rounding errors:
cat positions_31N.csv | cut -d, -f1-2,5,6 | sed 's/,/ /g' | sed 's/^726//g' | sed 's/ 5237/ /g' >
The last step is to add the header expected by OriConvert when reading a text file including the coordinates and filename of each image: add as the first line of the line #F=X Y Z N which states that the comments start with #, that we give longitude, latitude and altitude followed by the file name.
The micmac processing sequence is then
Pr="DJI_041[2-9]{1}.JPG"   # subset for lens calibration

mm3d OriConvert OriTxtInFile jmfgps

mm3d Tapioca MulScale "$P" 300 1500
mm3d Tapas RadialBasic "$Pr" Out=Cal
mm3d Tapas AutoCal "$P" InCal=Cal Out=Init
CenterBascule "$P" Init jmfgps tmp 
mm3d AperiCloud "$P" tmp

mm3d Malt Ortho "$P" tmp "DirMEC=Resultats" ZoomF=1 ZoomI=32 Purge=true AffineLast=false DefCor=0.001
mm3d Tawny Ortho-Resultats/
Nuage2Ply Resultats/NuageImProf_STD-MALT_Etape_8.xml Attr="Ortho-Resultats/Ortho-Eg-Test-Redr.tif"

Left: reference image provided by Clement overlaid with the flight path and the position at which each picture was acquired. Right: coarse point cloud.
After a 24-hour computation, the dense cloud result also generates an orthophoto (full resolution: 115 MB)and DEM (full resolution:: 11.3 MB).

The consistency of the orthophoto is assessed by overlapping with a reference image recovered from IGN: the error on the observatory tower position is 6.8 m.

Left: elevation map from the point cloud. Middle: orthophoto over the reference image. Right: transparent orthophoto over reference image.

In the previous example we used the MulScale option of Tapioca, which might take a significant amount of time by matching all pictures with each other when the number of pictures rises significantly. When georeferenced pictures are used, OriConvert can generate a list of matching pictures to ease the job of Tapioca which is then fed the input file with the File method: using this strategy, we restrict to 2550 analysis instead of 14878 (C^173_2) combinations. However, we do not know how the sequence
mm3d OriConvert OriTxtInFile jmfgps MTD1=1 NameCple=FileImagesNeighbour.xml CalcV=1 ImC=DJI_0554.JPG NbImC=9
mm3d Tapioca File "FileImagesNeighbour.xml" -1
selects matching images. We have analyzed the output file FileImagesNeighbour.xml using some awk and sed scripts to generate a WKT file fed to QGIS (following the documentation at this web page) for drawing the lines between coordinates of associated pictures. The result is given (left: IGN image ; right: our image) below


The GPR records (traces recorded every 10 cm along the fast axis, and every 40 cm along the slow axis) are processed with the following GNU/Octave script, despite the RD3 library for Matlab being broken with the latest (4.0) release of GNU/Octave. The script files have been modified to be functional with GNU/Octave: run go_archaeology to plot the A-scan (returned voltage v.s time), B-scan (voltage is coded as a color, fast axis is vertical and slow axis is horizontal) and some cross sections at various depths.
Note that trace 12 is erroneous and should be deleted from the directory including all files being processed.
The result of the GNU/Octave script is given below (from left to right: A-scan, B-scan and C-scan, with depth indication assuming a 100 m/us velocity). The buried path leaving the library and going towards the road is nicely visible on the top left part of each C-scan in the shallowest cross sections.

The last picture is a pointcloud generated from the GPR data and displayed in CloudCompare.

Data fusion

After proper georeferencing, both point clouds can be integrated in a common framework for display with CloudCompare: