High resolution historical topography: getting more from archival aerial photographs

High resolution elevation data is a fundamental information for multiple applications in geomorphology spanning from visual analyses to modeling. Nowadays, gathering of high-quality elevation data relies on multiple sensors and technologies which can be mounted on terrestrial, aerial and satellite platforms. In the last years, the Structure from Motion (SfM) algorithms have made possible the acquisition of high and very-high resolution elevation data from optical images acquired with high overlapping rates at virtually no cost. Such a feature made it possible to exploit remote sensing archival imagery to build historical topographic information with unprecedented detail. However, despite the increasing number of applications of SfM algorithms in the scientific literature, still little has been done in terms optimization and quality evaluation of the results. We have applied the SfM algorithm developed in the photogrammetric open source software MicMac to six black and white aerial photographs taken in 1954 at 1:33.000 in a mountainous and steep area in Central Italy. The aim of the experiment consists in a quantitative evaluation of the digital surface models obtained for different scanning settings.


Introduction
Photogrammetry allows the reconstruction of the shape and metric of an investigated object, starting from measurements of points on photographs. The recent development of photogrammetric algorithms such as Structure from Motion (SfM) and multi-view-stereo (MVS) algorithms opens an unprecedented possibility to derive high-resolution Digital Surface Models (DSMs) from archival images despite the common lack of the calibration certificates [1][2][3][4].
Although the number of papers dealing with SfM algorithms in the scientific literature is increasing, more work still needs to be done to evaluate the quality of the resulting DSMs and to optimize the digitalization process [1,5]. In this work, we test how photographs digitalization settings (geometric and radiometric resolution) affect the accuracy of the output DSMs obtained applying the same photogrammetric processing chain. We processed six archival aerial photographs taken in 1954 at 1:30.000 scale in a mountainous and steep area in Central Italy, which topography represents a challenge for producing photogrammetric DSMs.
Both proprietary and open source photogrammetric software applications exist, which are used for work and research [6][7][8]. Unlike proprietary software, however, open source gives the possibility to examine [9,10], modify and redistribute the algorithms, and is particularly suitable for research reproducibility and replicability [11] and closer to the general concept of open science [12]. For this work, we have adopted MicMac, an open source photogrammetric software that implements SfM and and MVS algorithms [13].

Materials
The area of interest is covered by six black and white 23x23cm aerial photographs (printings of original negative films), acquired in September 1954 at approximately 1:33.000 scale with a calibrated focal length of 153.78mm, using a FAIRCHILD camera. The photographs belong to two strips, which have an overlap of 60% and a sidelap of 40%.
For the computation of the DSMs, we used an Asus G750JX laptop, with a Intel c Core TM i7 4700HQ 8-cores processor, 24GB of DDR3 RAM memory and NVIDIA R GeForce R GTX 770M GPU graphic chip-set with 3GB DDR5 dedicated memory.
The digitalization of the aerial photographs was carried out by the Epson expression 10000 XL scanner. Its maximum geometric resolution is 2.400 dpi at a maximum of 16 bit in gray-scale and 32 bit in color mode.
The survey to collect ground control points (GCPs) and control points in kinematic mode was carried out using the GNSS receiver Leica Zeno 20, mounting the external antenna Leica AS10, and it was used in RTK mode, which allowed for a centimetric accuracy.
The software used at different steps of this work were: (i) SilverFast 6.6, (ii) MicMac [14].

Procedure to obtain the photogrammetric DSM
The first step of the procedure consisted in scanning the photographs. Each photograph was acquired at 400, 800, 1600 and 2400 dpi in gray-scale both at 8 and 16 bit using SilverFast software, for a total of 8 datasets. Images were then cropped using a rectangle mask, with the sides tangent to the fiducial marks.
In the second step, the images exif proprieties were edited (MicMac command SetExif ) to provide information on focal length F mm and equivalent focal length F 35 for the internal calibration. The focal value was available from the distributor website (photograph example). The 35 mm equivalent focal length (F 35 ) was computed as follows: where W 35 is the film width (in millimeters) for a 35 mm camera, and W f ilm is the film width of the acquisition camera (in millimeters). In this case, F 35 was equal to 23.86 mm.
The third step consists in the tie points definition (homologous pixels automatically detected across different images) through an image matching process performed by the MicMac tool Tapioca.
In the next step masks can be created to prevent MicMac from unnecessary computation. The graphical interface tool used is SaisieMasqQT. Then, tie-points are filtered out by HomolFilterMasq.
The bundle block adjustment (BBA) is performed in the fifth step by the MicMac command Tapas, an interface to the more complex tool Apero [15].
In the sixth step, GCPs have to be identified and selected on each image by a graphic interface (tool SaisieAppuisInitQT). The GCP ground absolute coordinates and the user-defined uncertainties are stored in a .txt file, then converted in a .xml file by the tool GCPConvert.
In the next step, a rigid transformation from the model coordinate system Arbitrary to the absolute coordinate system GroundInt is performed by the tool GCPBascule. Subsequently, a second BBA that fits tie points and GCPs is performed by the command Campari which calls the tool Apero.
Lastly, once the sparse point cloud has been generated, MicMac can proceed to the dense matching by the tool Malt. Finally, the dense point cloud is exported in .ply format by Nuage2Ply.

Ground Control Points
To further constrain the model, GCPs were collected by a GNSS survey. When collecting GCPs, care was taken to avoid cluster effects, and covering as homogeneously as possible the study area, also in elevation [16]. Particular attention was paid to the portion of the area covered by the black outlined rectangle in Figure 1. Due to the poorly built up landscape around the only municipality in the area (Visso) in 1954, and the difficulty in accurately finding the corresponding location of GCPs in today's landscape, most of the GCPs was selected at crossroads or on structures (Figure 1). For each of the 7 GCPs collected, longitude, latitude and elevation above ETM 2008 geoid were acquired using the available GNSS equipment in RTK mode. Elevation was converted to ITA 2005 geoid, more accurate for Italy.
To collect CPs for the DSMs validation, a GNSS survey was conducted in the field using the true kinematic method. The antenna was mounted on a car, and the frequency of acquisition set to 1 point every 10 meters. The over 1,000 points in the area were collected along two roads approximately running in the N-S and E-W directions (Figure 1). The points with uncertainties greater than 10 cm were filtered out, obtaining a robust set of elevation data used as control points (CP K ) to validate the output DSMs. Table 1 summarizes the results of the experiment presented. For reference, the first four rows describe the characteristics of the eight datasets.

Results and discussion
When considering the number of tie points generated, no appreciable difference is observed across color depths. The highest number of tie points were obtained for 800 dpi (both 8 and 16 bit) images. Furthermore, a decreasing trend in BBA error is observed across scanning resolutions. However, the improvement between the 1600 and 2400 dpi datasets is negligible, both at 8 and 16 bit.  Table 1 also reports the running time of Tapioca and Malt, the most demanding (parallelized) tools of the entire pipeline. The time increase trend is steeper for Malt compared to Tapioca. Therefore, processing images at high geometric resolutions implies more time expense due to the dense matching process. On the other hand, it is observed that increasing the radiometric resolution does not lead to running time increase, neither for Tapioca, nor for Malt.
The distance between the seven GCPs and the location of the same GCPs computed by the model [13] can be considered as an indicator of its internal coherence. In detail, Table 1 shows the two indicators MAE and RMSE. Inspection of Table 1 reveals that (i) g8 datasets have residuals comparable to or slightly smaller than the corresponding g16 datasets, (ii) the 800g8 dataset shows the lowest residuals, (iii) the minimum values of RMSE and MAE were obtained for 800 dpi (both for 16 bit and 8 bit scanned images). Unexpectedly, we notice that increasing color depth does not necessarily correspond to an increased internal coherence. We hypothesize that this can be due to the fact that the original images are grey-scale printed photographs. Since previous research indicate that 16 bit acquisition is advised as a best practice [17], based on our results, further tests should be carried out to further look into this topic.
The total 8 dense point clouds were projected into DSMs, which resolution was set close to the GSD of the starting images (Table 1).
To evaluate the accuracy of the eight DSMs, we computed the difference between the elevation of each CP K and the elevation of each DSM at the planimetric location of each CP K . Outliers were singled out and filtered according to a three-sigma (3 σ) error criterion. Statistics of residuals distributions is shown in the last section of Table 1.
Inspection of model validation data (Table 1) shows that 800g8 is the best dataset in terms of mean value proximity to zero, which can be taken as indicator of small systematic error. Anyway, as already mentioned above, 1600g8 has the smallest standard deviation, which can be taken as indicator of good relative accuracy. Moreover, 1600g8 has the smallest root mean squared error (RMSE), which takes into account both random and systematic error [18]. Finally, 1600g8 has also the smallest mean absolute error (MAE), which is an even more important quality indicator than RMSE, according to some authors [19]. To check for normality of the two best distribution candidates 800g8 and 1600g8, we made a one-sample Kolmogorov-Smirnov test. The 800g8 dataset seems better normally distributed, showing a maximum distance value of 0.20 compared to the 0.30 of the 1600g8 dataset.
Overall, considering the model internal coherence, the mean point density of the point clouds, and the results of the model validation (Table 1), it emerges that the optimal performance was obtained at 1600 dpi for images scanned at 8 bit, whereas it was obtained at 800 dpi for images scanned at 16 bit. When compared to the computation time, it is also clear that a low increase in the internal model of 2400 dpi (8 and 16 bit) corresponds to a computation time way too large compared to the 1600 dpi datasets. Such an evidence underlines that acquisition at 2400 dpi appears inconvenient, even for high-performance machines. Moreover, Table 1 shows that the number of tie points automatically detected by Tapioca decreases with the increasing scanning resolution higher than 800 dpi.

Conclusions
This work quantified the influence of acquisition of archive aerial images on the accuracy of DSMs produced applying photogrammetric algorithms, and hence can help to define the best acquisition mode to get the most from archival aerial photographs.
Limitedly to the historical aerial photographs used in this study (available for the entire Italian territory), results revealed that the optimal scanning geometrical resolution is between 800 and 1600 dpi, and that the radiometric resolution does neither affect computation time nor provide any quality improvement of the final DSMs.
The photographs used in this experiment belong to an extensive aerial photographic survey conducted between 1954 and 1956 over all the Italian territory. Based on this study, such large set of historical aerial photographs can be exploited to produce DSMs at a resolution between 0.5 and 1 m with errors close to 2 m in elevation. Such an invaluable dataset could provide accurate high-resolution quantitative information on topography, which is of great value for multiple studies in Earth Science.