On the Use of Historical Flights for the Urban Growth Analysis of Cities Through Time: The Case Study of Avila (Spain)

: Historical aerial images are a unique and relatively unexplored means of deriving spatio-temporal information for scenes and landscapes. Such historical imagery can be combined with photointerpretation and image-based 3D modelling techniques, providing the fourth dimension of time to 3D geometrical representations. This allows urban planners, historians, and other specialists to identify, describe, and analyse changes in scenes and landscapes. Urban growth has an important impact on the sustainable development of cities. An important step for the analysis of urban growth is the identiﬁcation of di ﬀ erent urban sectors. To this end, this paper proposes a methodology for the 4D urban growth analysis of cities through time using a free and open source software developed by the authors. This approach uses the latest advances in photogrammetry, including the so-called incremental Structure from Motion, to evaluate the urbanistic changes of a city by means of confronting two-point clouds from di ﬀ erent eras. The objectives of this paper are twofold: (i) ﬁrst, the processing of historical aerial images using modern photogrammetric techniques; (ii) second, deriving spatio-temporal information for urban cities, o ﬀ ering a method for researchers to identify changes over time. In order to validate this method, the urban growth of the city of Avila between 1956 and 2017 was assessed taking the historical American ﬂight of 1956 and the digital aerial ﬂight of 2017. The results were statistically assessed according to georeferencing quality, conﬁrming that the approach developed can be used to support urban growth analysis through time and providing relevant data in 2D and 3D.


Introduction
Cities and towns are temporary expressions of the social, economic, and cultural history of the population that lives and works within them. These spaces have been formed throughout history and are represented by the construction of new urban spaces; the creation of new infrastructures (e.g., roads or bridges); or even the expansion and remodelling of existing spaces, mainly historical districts [1]. These new urban spaces play an important role for the future development and sustainability of these cities and towns [2,3]. Trying to find a balance for the urban development of cities, urban growth Sustainability 2020, 12, 4673 3 of 17 The paper has been structured as follows: after this introduction, Section 2 describes the materials and methods used and Section 3 outlines and discussed the results obtained. A final section is devoted to highlighting the main conclusions and future trends.

Study Area
The city of Avila is located in the south of the community of Castilla y Leon (Spain) with a population of 58,149 people in 2017, according to the National Statistic Institute.
The urban framework of the city (Figure 1) was established considering five different sectors ( Figure 1): (i) the intramural sector which encloses the historical center and is bordered by the medieval wall of Avila; (ii) the extramural sector surrounding the intramural sector; (iii) the industrial sector which spreads out in the east of the city; (iv) the old residential sector  beside the extramural sector in the east and south parts of the city; (v) the modern residential sector , which is spread out along the peripherical part of the city, especially towards the south and east. On the other hand, towards the north and considering the barrier of the railway, the predominant buildings present are semi-detached houses, with areas of low density around the Catholic University of Avila, as well as the area that appears around the river banks in the south of the city.
New buildings of one or two floors, such as urbanizations of single-family houses, appeared in the new residential areas located to the east of the city together with the main hospital and several shopping centres ( Figure 11). Analysis of urban growth in the industrial area: The industrial area of Avila emerged in the east, close to the city's main routes of communication routes-the railway station and the highway to Madrid. In its surroundings, new residential areas were created and hospitals and commercial areas also appeared. This industrial area is highlighted when comparing both flights through buildings of large surface areas and single ground floors (Figure 11), being the result of the application of the PGOU carried out in 2005 and 2007, where a total of 150 Ha were classified as industrial areas.
After analysing the urban growth in the different areas of the city using the proposed method, we can conclude that the city of Avila has undergone a gradual process, expanding the built space from the 1950s to the present, going through different phases related to the geomorphological characteristics of the territory and the populational movements between the centre and the periphery. In this way, from the intramural zone and its surroundings outside the walls, traditional In general, the city of Avila expanded peripherally between 1956 and 2017, except to the west, where the river acts as a natural barrier to urban expansion.

Historical American flight
The so-called American flight dates back as the oldest flight in the national territory of Spain. This flight was executed between 1929 and 1930 with a scale of 1:10,000. Although its initial spatial coverage was limited, more photogrammetric flights were carried out by the United States and the United Kingdom air forces (USAF and RAF) after the end of the Second World War, with the aim of mapping out European territories. The first flights, known as "Series A", were made between 1945 and 1946 with a scale of 1:43,000. "Series B" flights were carried out between 1956 and 1957 with a scale of 1:33,000, and "Series C" between 1967 and 1968 with a scale of 1:18,000 [31].
Spain was a key territory in the eyes of the United States due to its geostrategic position, establishing specific agreements between both countries on the availability of photogrammetric flights.
Regarding the case study evaluated, three historical black and white images with a scale of 1:33,000 were analysed. These American flight images were executed by the USAF in collaboration with the SGE and the IGN between March 1956 and September 1957 ( Figure 2). A Beech-craft RC-45 turboprop aircraft was used for these flights. The photogrammetric camera employed was a Fairchild T-11 with a format size of 22.8 cm X 22.8 cm. This camera was equipped with a Metrogon II type lens with a focal length of 152 mm, a fixed focus, a wide angle (73º field of view), a fixed aperture f:6.3, a central shutter type Rapidyne (1/10 to 1/500), and a yellow filter with built-in anti-vignetting. Regarding the metadata existing on the historical images, the focal length (153.55 mm), lens serial number (XF 599), camera serial number (52049), and flight height (about 5500 m) are known. The approximate scale used was 1:33,000, with a stereoscopic forward overlap of 60% and a side overlap of 30%. The photographic film of this flight is kept safe at the Cartographic and Photogrammetric Center of the Air Force (CECAF) and was digitized in September 2011 with a spatial resolution of 21 microns.
These types of flights were planned in order to obtain a good accuracy in the restitution of points on the ground (i.e., a good base/height ratio to guarantee a good intersection of the perspective rays). Therefore, the uncertainty associated with the internal parameters of the cameras used at that time, especially the physical parameters (e.g., radial lens distortion), together with the large baselines (i.e., separation between images), would not contribute to the application of new photogrammetric methods based on dense matching, where dealing with small baselines between images and large overlaps (>80%) is crucial.
In this paper, the American flight "Series B" (1956)(1957) was used as a basis to quantify urban changes through time, and it is possible to extrapolate this methodology to other urban centres.

PNOA Flight
The National Plan of Aerial Orthophotography (PNOA) is a cooperative project co-financed by the general administration of the Spanish state and the autonomous communities that were born following the trend set by the Infrastructure for Spatial Information in Europe (INSPIRE) in 2004. These projects aim to obtain photogrammetric products, orthoimages, along the entire national territory, with a ground sample distance (GSD) of 25 or 50 cm, allowing the obtaining of cartographic scales of 1:5,000 and 1:10,000, respectively, and high-resolution digital elevation models (DEM) [32].
For the area of study, a total of 21 colour aerial images were used from the PNOA 2017 flight with a GSD of 25 cm (pixel size) ( Figure 3). A Beech-craft RC-45 turboprop aircraft was used for these flights. The photogrammetric camera employed was a Fairchild T-11 with a format size of 22.8 cm X 22.8 cm. This camera was equipped with a Metrogon II type lens with a focal length of 152 mm, a fixed focus, a wide angle (73º field of view), a fixed aperture f:6.3, a central shutter type Rapidyne (1/10 to 1/500), and a yellow filter with built-in anti-vignetting. Regarding the metadata existing on the historical images, the focal length (153.55 mm), lens serial number (XF 599), camera serial number (52049), and flight height (about 5500 m) are known. The approximate scale used was 1:33,000, with a stereoscopic forward overlap of 60% and a side overlap of 30%. The photographic film of this flight is kept safe at the Cartographic and Photogrammetric Center of the Air Force (CECAF) and was digitized in September 2011 with a spatial resolution of 21 microns.
These types of flights were planned in order to obtain a good accuracy in the restitution of points on the ground (i.e., a good base/height ratio to guarantee a good intersection of the perspective rays). Therefore, the uncertainty associated with the internal parameters of the cameras used at that time, especially the physical parameters (e.g., radial lens distortion), together with the large baselines (i.e., separation between images), would not contribute to the application of new photogrammetric methods based on dense matching, where dealing with small baselines between images and large overlaps (>80%) is crucial.
In this paper, the American flight "Series B" (1956)(1957) was used as a basis to quantify urban changes through time, and it is possible to extrapolate this methodology to other urban centres.

PNOA Flight
The National Plan of Aerial Orthophotography (PNOA) is a cooperative project co-financed by the general administration of the Spanish state and the autonomous communities that were born following the trend set by the Infrastructure for Spatial Information in Europe (INSPIRE) in 2004. These projects aim to obtain photogrammetric products, orthoimages, along the entire national territory, with a ground sample distance (GSD) of 25 or 50 cm, allowing the obtaining of cartographic scales of 1:5,000 and 1:10,000, respectively, and high-resolution digital elevation models (DEM) [32].
For the area of study, a total of 21 colour aerial images were used from the PNOA 2017 flight with a GSD of 25 cm (pixel size) ( Figure 3). The large format digital camera used was the VEXCEL-Ultracam-XP with a focal length of 100.5 mm; a gyro-stabilized platform, GSM3000; a global navigation satellite system (GNSS) supported by an inertial measurement unit (IMU); an AeroControl Iid with a 256 Hz frequency and a GPS positional unit with a 2Hz frequency. The configuration of the imagery panchromatic sensor was a frame with a resolution of 17310 x 11310 pixels and a physical pixel size of 6 microns. The following table (Table 1) describes the characteristics of both flights.

Methodology
The photogrammetric process for the dense reconstruction was performed with the open source software GRAPHOS [33], applying four consecutive steps ( Figure 4): (i) the preprocessing of the images, (ii) the extraction and matching of key points, (iii) the orientation and self-calibration of the images, and (iv) dense reconstruction. The large format digital camera used was the VEXCEL-Ultracam-XP with a focal length of 100.5 mm; a gyro-stabilized platform, GSM3000; a global navigation satellite system (GNSS) supported by an inertial measurement unit (IMU); an AeroControl Iid with a 256 Hz frequency and a GPS positional unit with a 2Hz frequency. The configuration of the imagery panchromatic sensor was a frame with a resolution of 17310 x 11310 pixels and a physical pixel size of 6 microns. The following table (Table 1) describes the characteristics of both flights.

Methodology
The photogrammetric process for the dense reconstruction was performed with the open source software GRAPHOS [33], applying four consecutive steps ( Figure 4): (i) the preprocessing of the images, (ii) the extraction and matching of key points, (iii) the orientation and self-calibration of the images, and (iv) dense reconstruction.

Preprocessing of Images
Image preprocessing is an important step, since it can provide a means for better feature extraction and matching. This is specifically relevant in those cases where the texture quality is unfavourable due to low radiometric resolution, rotations, and large temporal differences, especially with the case of the American flight. In order to homogenize the results, all the input images were preprocessed based on the Recursively Separated and Weighted Histogram Equalization (RSWHE) algorithm [34], since RSWHE preserves the image brightness more accurately and produces images with better contrast enhancement. While this step is considered optional in several photogrammetric processing pipelines, it is highly recommended when applied for feature extraction. This algorithm has particular success in obtaining better results in 3D point clouds.

Feature Extraction and Matching
The extraction and matching of the key points between images are essential requirements for dense reconstructions. Feature extraction consists of the identification of several meaningful features within each image. This is additionally important in ensuring that the points are re-detectable in different images. Detectors are operators which search for 2D locations within images (key points), whereas descriptors analyse the surroundings of the detected key points and produce a corresponding 2D vector of information. This information can then be used to quickly classify the extracted points or matching processes [35]. Several algorithms based on detectors and descriptors have been proposed over the last decades in order to reliably detect features among images with geometric and radiometric transformations. Among the different detectors/descriptors available in GRAPHOS [33], the algorithm that best adapted to the type of historical images used was Scale Invariant Fine Transform (SIFT) [36]. Once the key points were extracted, the homologous points were matched between the images. For this purpose, the descriptor also provided by SIFT was used. The features codified by the descriptor were matched, first according to the Euclidean distance (L2-Norm) and secondly filtered by an optimized random sampling algorithm (ORSA) [37]. This algorithm is a variant of the random sample consensus (RANSAC) [38], with an adaptive criterion to filter the mismatches by the employment of epipolar geometry constraints.

Preprocessing of Images
Image preprocessing is an important step, since it can provide a means for better feature extraction and matching. This is specifically relevant in those cases where the texture quality is unfavourable due to low radiometric resolution, rotations, and large temporal differences, especially with the case of the American flight. In order to homogenize the results, all the input images were preprocessed based on the Recursively Separated and Weighted Histogram Equalization (RSWHE) algorithm [34], since RSWHE preserves the image brightness more accurately and produces images with better contrast enhancement. While this step is considered optional in several photogrammetric processing pipelines, it is highly recommended when applied for feature extraction. This algorithm has particular success in obtaining better results in 3D point clouds.

Feature Extraction and Matching
The extraction and matching of the key points between images are essential requirements for dense reconstructions. Feature extraction consists of the identification of several meaningful features within each image. This is additionally important in ensuring that the points are re-detectable in different images. Detectors are operators which search for 2D locations within images (key points), whereas descriptors analyse the surroundings of the detected key points and produce a corresponding 2D vector of information. This information can then be used to quickly classify the extracted points or matching processes [35]. Several algorithms based on detectors and descriptors have been proposed over the last decades in order to reliably detect features among images with geometric and radiometric transformations. Among the different detectors/descriptors available in GRAPHOS [33], the algorithm that best adapted to the type of historical images used was Scale Invariant Fine Transform (SIFT) [36]. Once the key points were extracted, the homologous points were matched between the images. For this purpose, the descriptor also provided by SIFT was used. The features codified by the descriptor were matched, first according to the Euclidean distance (L2-Norm) and secondly filtered by an optimized random sampling algorithm (ORSA) [37]. This algorithm is a variant of the random sample consensus (RANSAC) [38], with an adaptive criterion to filter the mismatches by the employment of epipolar geometry constraints.
In order to optimize this matching process, low-resolution matching was carried out in the first instance. Subsequently and considering the previously obtained results through matching, a full resolution matching was performed.

Camera Orientation and Self-Calibration
Knowing the 2D coordinates of the matching points in two or more images, it was possible to compute the external parameters of the camera (spatial and angular position). For this purpose, a hierarchical procedure defined by [39] was applied: (i) a relative orientation between the images through the fundamental matrix [40] and (ii) an absolute orientation (georeferencing) of the images through an iterative minimization process (bundle adjustment) that used the condition of collinearity (Equation (1)): Collinearity equations establish a functional and stochastic model where the internal camera parameters (f, x´0, y´0, ∆x´, ∆y´), the external camera parameters (X 0C , Y 0C , Z 0C , r ij ), and the object coordinates (X,Y,Z) are the unknown variables to be solved, being (x',y') the known coordinates of the matched points. Please note that the bundle adjustment is non-linear; thus, an iterative process is required. For each iteration, the estimated results (the solution and its co-factor matrix) for all the unknown registration parameters are updated [41].
Regarding the absolute orientation, a series of ground control points (GCPs) was considered to compute the georeferencing of the images, as well as ground checkpoints (GChPs) that assessed the validity of the adjustment. More specifically, GCPs and GChPs corresponding to urban and rustic areas were used as the reference elements, trying to guarantee their static permanence over time. In rustic areas, elements such as rocks and hermitages were considered as references, while in urban areas existing historical monuments in the city of Avila were taken as reference. The different GCPs and GChPs were surveyed with a dual frequency GPS system.

Dense Point Cloud Generation
The photogrammetric workflow ended with dense matching, where 3D points of the surface were computed for each image pixel. In our case, dense point cloud generation was carried out using a Semi-Global Matching (SGM) approach [42]. Most of these strategies are supported by the use of the camera pose together with the epipolar geometry, minimizing the energy function (Equation (2)). The dense matching approach used was the MicMac algorithm [43], ideal for aerial vertical images, compared to the Patch-based Multi-view Stereo (PMVS) algorithm [44], which is more suitable for oblique images.
This strategy consists of minimizing the energy function E (Equation (2)) along the eight basic directions that a pixel can take (every 45º). This energy function is supported by a matching cost function, M (pixel matching cost), that accounts for the degree of similarity of the pixels between two images, x and x', together with the incorporation of two constraints, P 1 and P 2 , which account for the possible presence of gross errors in the SGM process. Likewise, a third constraint was added to the SGM process, consisting of epipolar geometry from photogrammetry [45] that allows the reduction of the search space of each pixel and thus reduces the computational cost.
where E(D) corresponds to the energy function to be minimized based on the disparity (parallax) between the matching points, M (pixel matching cost) assesses the degree of similarity between pixels x and its homologous x´through its disparity or parallax D p , while the terms P 1 and P 2 correspond to two constraints that avoid the presence of gross errors in the dense matching process.

Accuracy Assessment
Once both dense point clouds were obtained and georeferenced in the same coordinate system, an analytical comparison of both point clouds was performed to observe the urban growth of the city between 1956 and 2017. In addition, a detailed analysis of the Z coordinates allowed for the analysis of the emergence of new buildings and industries, as well as the classification of buildings according to the number of floors in the city of Ávila between 1956 and 2017.
Since the point cloud obtained from the American Flight was less dense than the point cloud generated with the PNOA flight, the former was converted into a mesh to perform a cloud-to-mesh comparison. The approach used for mesh triangulation was the Poisson surface reconstruction [46], in which the octree level was established according to the spatial resolution of the point cloud [47]. The result of the comparison was the triple scalar field (X,Y,Z) corresponding to the three components of the deviation vector, which represented the absolute distance between each compared point and its closest point in the reference mesh.
For the assessment of accuracy within this data, discrepancies associated with every 3D photogrammetric point can be assessed by their central tendency and dispersion. Nevertheless, errors or discrepancies derived from photogrammetric data do not usually follow a Gaussian distribution, which is attributable to the presence of gross errors [48]. This lack of normality can be checked by the use of visual methods, such as quantile-quantile plots [49], or numerically by the use of normality tests [50]. Once the nature of the error distributions has been established, accuracy assessments are carried out on the basis of robust estimators such as the median, m; the normalized median absolute deviation, NMAD (Equation (3)); the square root of the biweight midvariance, BWMV (Equation (4)); and the inter-percentile ranges, (IPR). These can be described by [51]: where the median absolute deviation (MAD) (Equation (7))-that is, the median (m) of the absolute deviations from the data's median (m x )-is defined as: Please note that due to asymmetric distribution, a plus-minus range could not be provided, and therefore an absolute inter-percentile range at the 90% confidence interval was provided.

Experimental Results
Three main steps were followed to study the urban growth in the city of Avila between 1956 and 2017 ( Figure 3): (i) the generation of dense point clouds from both flights, (ii) the comparison of the georeferenced 3D models, (iii) an urban growth analysis.

Suitability Assessment
Twenty-one aerial images coming from the 2017 PNOA flight and three aerial images coming from the American flight were used to obtain the 3D reconstruction of the city of Avila between 1956 and 2017 by means of the approach defined in Section 3. The full photogrammetric workflow was carried out using the free and open source photogrammetric software, GRAPHOS [33], supported by a scientific project of the International Society of Photogrammetry and Remote Sensing (ISPRS).
First, a SIFT algorithm was applied in order to obtain the key points (Table 2). Then, a robust matching strategy to identify homologous key points was carried out based on L2-Norm, ORSA, and epipolar geometry. These key points were later used to perform the relative orientation and self-calibration of both the photogrammetric flights (Table 2). During the absolute orientation, several GCPs and GChPs were used with the aim of georeferencing the final product and assessing the accuracy, respectively ( Figure 5). Regarding the PNOA flight, 16 GCPs and 25 GChPs were used. These GCPs and GChPs were surveyed with the Leica 1200GPS system using the real time kinematic (RTK) approach and taking the European Terrestrial Reference System 1989 (ETRS89) and the Universal Transverse Mercator (UTM) 30N projection as a reference system and national map projection, respectively. A total of 10 GCPs and 26 GChPs were used for the historical American flight. Three main steps were followed to study the urban growth in the city of Avila between 1956 and 2017 ( Figure 3): (i) the generation of dense point clouds from both flights, (ii) the comparison of the georeferenced 3D models, (iii) an urban growth analysis.

Suitability Assessment
Twenty-one aerial images coming from the 2017 PNOA flight and three aerial images coming from the American flight were used to obtain the 3D reconstruction of the city of Avila between 1956 and 2017 by means of the approach defined in Section 3. The full photogrammetric workflow was carried out using the free and open source photogrammetric software, GRAPHOS [33], supported by a scientific project of the International Society of Photogrammetry and Remote Sensing (ISPRS).
First, a SIFT algorithm was applied in order to obtain the key points (Table 2). Then, a robust matching strategy to identify homologous key points was carried out based on L2-Norm, ORSA, and epipolar geometry. These key points were later used to perform the relative orientation and selfcalibration of both the photogrammetric flights (Table 2). During the absolute orientation, several GCPs and GChPs were used with the aim of georeferencing the final product and assessing the accuracy, respectively ( Figure 5). Regarding the PNOA flight, 16 GCPs and 25 GChPs were used. These GCPs and GChPs were surveyed with the Leica 1200GPS system using the real time kinematic (RTK) approach and taking the European Terrestrial Reference System 1989 (ETRS89) and the Universal Transverse Mercator (UTM) 30N projection as a reference system and national map projection, respectively. A total of 10 GCPs and 26 GChPs were used for the historical American flight. According to an ideal equilateral triangular distribution for a circular neighbourhood [35], the average point density achieved for the PNOA flight was 89.7 cm (confidence interval of 85.0 -95.3 cm), whereas for the historical American flight the point density was 250.5 cm (confidence interval of 232.5 -273.6 cm). According to an ideal equilateral triangular distribution for a circular neighbourhood [35], the average point density achieved for the PNOA flight was 89.7 cm (confidence interval of 85.0-95.3 cm), whereas for the historical American flight the point density was 250.5 cm (confidence interval of 232.5-273.6 cm).
As a result, both of the point cloud densities were registered and georeferenced under the same reference system, ETRS89 ( Figure 6).
With the aim of establishing a range of confidence during the comparison stage of both the photogrammetric point clouds, a statistical analysis was carried out (Tables 3 and 4). The Gaussian parameters (i.e., the mean and standard deviation) were estimated with a 68% confidence interval for the different axes, as well as the 3D error. Regarding the robust parameters, the central tendency of the error was estimated by the median and the error dispersion as the square root of the biweight midvariance (Equation (4)).
Sustainability 2020, 12, x FOR PEER REVIEW 10 of 18 As a result, both of the point cloud densities were registered and georeferenced under the same reference system, ETRS89 (Figure 6). With the aim of establishing a range of confidence during the comparison stage of both the photogrammetric point clouds, a statistical analysis was carried out (Tables 3 and 4). The Gaussian parameters (i.e., the mean and standard deviation) were estimated with a 68% confidence interval for the different axes, as well as the 3D error. Regarding the robust parameters, the central tendency of the error was estimated by the median and the error dispersion as the square root of the biweight midvariance (Equation 4).

Gaussian Estimation (m) Robust Estimation (m)
Number  Table 4. Georeferencing error estimation of the American flight according to the Gaussian and robust parameters. Units in meters.  Table 3 presents the comparison of the Gaussian parameters (mean and standard deviation) and the robust parameters (median and biweight midvariance) for the georeferenced PNOA flight. Please note the overestimation of the 3D error for the Gaussian estimation, both in the GCPs and GChPs. Moreover, the robust parameters provide a clearer view for the error distribution. For instance, along the x-axis the GCPs and GChPs dispersions were overestimated (±15 cm vs. ±10 cm and ±21 cm vs. ±14 cm for the GCPs and GChPs, respectively). In general terms, the final georeferencing error was compatible with the technical specifications of the PNOA flight (Table 1).

Gaussian Estimation (m) Robust Estimation (m)
Regarding the georeferencing of the American flight (Table 4), the same overestimation can be detected. This is seen through the largest errors in 3D measurements using GChPs, with a discrepancy of almost 30 cm (3.15 m vs. 2.87 m). Due to the particular characteristics of the data obtained from the American flight (scale, pixel quality) and the uncertainty of the GChPs selection, the final error rose to values of almost 2 m for each axis with an unavoidable bias. To contextualize these values, please note that the spatial resolution of the scanned images produced, at least, an error source of 0.70 m (approximated by the scale size). This error was amplified by the uncertainty related to the manual identification of the GCPs and the GChPs, which could vary empirically up to 2 to 3 times. These error values were like those obtained by similar studies [52]. Similarly, in [53] 1947 historical images were scanned with a resolution of 21 microns to generate an orthomosaic with a mean error of 1.97 m and root mean square error (RMSE) of 1.00 m, with errors ranging between 0.90 m to 5.88 m. Table 5 presents MAD and NMAD (Equation (3)) values, as well as the absolute inter-percentile range at the 90% confidence interval (useful for an asymmetric error distribution). It is worth noting that the difference between the percentiles 5% and 95% (90% of confidence level) is ca. 4.5 m. This error implies an acceptable value to be used for a mapping analysis at a scale of 1:22,500 or smaller.

Urban Growth Analysis
Beginning with the two photogrammetric point clouds (1965-2017) previously obtained, an analysis of the urban growth of the city of Avila, based on height discrepancies, was carried out (Figure 7). To this end, the following approach was applied: (i) extraction of the area of overlap between both point clouds, (ii) comparison between the PNOA point cloud and the mesh generated from the American flight point cloud, and (iii) classification by height. From the classification results obtained, we can see how the urban structure of the city of Avila has considerably changed from 1956 to 2017. In 1956, Avila had a population of around 23,600 inhabitants and approximately 2000 residential buildings, according to the data provided by the provincial registry. In 2017, the population of Avila reached around 59,000 inhabitants, with a total of 8000 residential buildings. In the year 2000, during the "real-estate bubble", almost 2300 residential buildings were built. From the classification results displayed in Figure 7, it can be highlighted that this urbanistic revolution was concentrated mainly towards the north and east areas of the city, in connection with the railway station and the main roads that connect Avila with the capital city of Madrid (Figure 7).
Next, a detailed analysis of urban growth is performed in the five different sectors outlined in Figure 4, trying to assess what this method can contribute in terms of urban development of the city between 1965 and 2017. From the classification results obtained, we can see how the urban structure of the city of Avila has considerably changed from 1956 to 2017. In 1956, Avila had a population of around 23,600 inhabitants and approximately 2000 residential buildings, according to the data provided by the provincial registry. In 2017, the population of Avila reached around 59,000 inhabitants, with a total of 8000 residential buildings. In the year 2000, during the "real-estate bubble", almost 2300 residential buildings were built. From the classification results displayed in Figure 7, it can be highlighted that this urbanistic revolution was concentrated mainly towards the north and east areas of the city, in connection with the railway station and the main roads that connect Avila with the capital city of Madrid (Figure 7).
Next, a detailed analysis of urban growth is performed in the five different sectors outlined in Figure 4, trying to assess what this method can contribute in terms of urban development of the city between 1965 and 2017.
Analysis of urban growth in the historical centre (intramural and extramural sector): In the historical centre, new buildings appeared occupying areas previously used for gardens, crops, and private corrals. Likewise, a modification of the existing ground floor buildings by higher buildings took place due to activities related to replacement, restoration, renovation, or rehabilitation. Nevertheless, these new buildings were limited in height by the general plan of urban planning (PGOU) of the city of Avila. Due to this, within the walled enclosure only buildings up to three floors appeared. Considering these intramural and extramural sectors, the most important modifications occurred to the west of the walled enclosure, an area where properties with unbuilt lots, gardens, and private corrals existed in the 1950s. Likewise, modifications of one or two heights can be observed in some buildings of the historical centre, reaching the maximum height allowed by the PGOU (Figure 8). Analysis of urban growth in the old residential area : The need to expand the city due to the increase in population in the 1960s was resolved by the appearance of new neighbourhoods near the historical centre, in which the constructive typology consisted of multi-family buildings in blocks of three and four floors. Construction at this time was characterized by the low quality of the buildings and the lack of infrastructure. Furthermore, with the natural barrier of the river in the west, the city continued growing up to the railway, which constituted a new barrier for the expansion of the city. As a result, new residential neighbourhoods emerged in the east and then in the south during the 1950s and 1980s. It was in these areas where new buildings with three and four floors appeared ( Figure 9). In addition to this, the comparison highlights the construction of two important public buildings relevant to the city, the Catholic University of Avila and the Palace of Congress, both in in the northwest of the city (red rectangle marked in Figure 9). Analysis of urban growth in the modern residential area after 1980: Subsequently and around the main surroundings of the city centre, new industrial polygons and residential areas (single-family homes) were built to the east. Towards the south of the city, medium-density residential elements Analysis of urban growth in the old residential area : The need to expand the city due to the increase in population in the 1960s was resolved by the appearance of new neighbourhoods near the historical centre, in which the constructive typology consisted of multi-family buildings in blocks of three and four floors. Construction at this time was characterized by the low quality of the buildings and the lack of infrastructure. Furthermore, with the natural barrier of the river in the west, the city continued growing up to the railway, which constituted a new barrier for the expansion of the city. As a result, new residential neighbourhoods emerged in the east and then in the south during the 1950s and 1980s. It was in these areas where new buildings with three and four floors appeared ( Figure 9). In addition to this, the comparison highlights the construction of two important public buildings relevant to the city, the Catholic University of Avila and the Palace of Congress, both in in the northwest of the city (red rectangle marked in Figure 9).
Analysis of urban growth in the modern residential area after 1980: Subsequently and around the main surroundings of the city centre, new industrial polygons and residential areas (single-family homes) were built to the east. Towards the south of the city, medium-density residential elements (multi-family buildings) were built together with the main city facilities (e.g., municipal stadium, hypermarkets, provincial hospitals, and universities, among others).
buildings and the lack of infrastructure. Furthermore, with the natural barrier of the river in the west, the city continued growing up to the railway, which constituted a new barrier for the expansion of the city. As a result, new residential neighbourhoods emerged in the east and then in the south during the 1950s and 1980s. It was in these areas where new buildings with three and four floors appeared ( Figure 9). In addition to this, the comparison highlights the construction of two important public buildings relevant to the city, the Catholic University of Avila and the Palace of Congress, both in in the northwest of the city (red rectangle marked in Figure 9). Analysis of urban growth in the modern residential area after 1980: Subsequently and around the main surroundings of the city centre, new industrial polygons and residential areas (single-family homes) were built to the east. Towards the south of the city, medium-density residential elements (multi-family buildings) were built together with the main city facilities (e.g., municipal stadium, hypermarkets, provincial hospitals, and universities, among others).
These areas congregated more diversity in terms of building types and the height of buildings. The tallest buildings with seven and eight floors appeared in the southern area of these modern residential neighbourhoods after 1980 ( Figure 10). It was also in these neighbourhoods where most of the buildings with five and six floors were congregated. These areas congregated more diversity in terms of building types and the height of buildings. The tallest buildings with seven and eight floors appeared in the southern area of these modern residential neighbourhoods after 1980 ( Figure 10). It was also in these neighbourhoods where most of the buildings with five and six floors were congregated. On the other hand, towards the north and considering the barrier of the railway, the predominant buildings present are semi-detached houses, with areas of low density around the Catholic University of Avila, as well as the area that appears around the river banks in the south of the city. New buildings of one or two floors, such as urbanizations of single-family houses, appeared in the new residential areas located to the east of the city together with the main hospital and several shopping centres ( Figure 11).  On the other hand, towards the north and considering the barrier of the railway, the predominant buildings present are semi-detached houses, with areas of low density around the Catholic University of Avila, as well as the area that appears around the river banks in the south of the city. New buildings of one or two floors, such as urbanizations of single-family houses, appeared in the new residential areas located to the east of the city together with the main hospital and several shopping centres ( Figure 11).
Analysis of urban growth in the industrial area: The industrial area of Avila emerged in the east, close to the city's main routes of communication routes-the railway station and the highway to Madrid. In its surroundings, new residential areas were created and hospitals and commercial areas also appeared. This industrial area is highlighted when comparing both flights through buildings of large surface areas and single ground floors (Figure 11), being the result of the application of the PGOU carried out in 2005 and 2007, where a total of 150 Ha were classified as industrial areas.
After analysing the urban growth in the different areas of the city using the proposed method, we can conclude that the city of Avila has undergone a gradual process, expanding the built space from the 1950s to the present, going through different phases related to the geomorphological characteristics of the territory and the populational movements between the centre and the periphery. In this way, from the intramural zone and its surroundings outside the walls, traditional neighbourhoods were created between the 1950s and 1980s, always to the east of the city. The expansion of the city continued from the 1980s, looking for open spaces, new building typologies, and properties of better quality and more affordable prices than in the downtown areas.
On the other hand, towards the north and considering the barrier of the railway, the predominant buildings present are semi-detached houses, with areas of low density around the Catholic University of Avila, as well as the area that appears around the river banks in the south of the city. New buildings of one or two floors, such as urbanizations of single-family houses, appeared in the new residential areas located to the east of the city together with the main hospital and several shopping centres (Figure 11). After analysing the urban growth in the different areas of the city using the proposed method, we can conclude that the city of Avila has undergone a gradual process, expanding the built space from the 1950s to the present, going through different phases related to the geomorphological characteristics of the territory and the populational movements between the centre and the periphery. In this way, from the intramural zone and its surroundings outside the walls, traditional

Conclusions
Historical images can be considered an important legacy of enormous value for multiple possible analyses. Nevertheless, the processing of these images presents a number of issues involving the need for located control points, the finding of which is often a tedious task and difficult to automate. This article presents a method supported by the latest advances in photogrammetry and computer vision as well as advances in the automatic analysis of 3D point clouds, with the aim of exploiting the historical images of cities for urban growth analysis. The present case study shows how the city of Ávila (Spain) has undergone a significant urban development in recent decades.
Although there is no doubt about the impact and value that historical data can provide in many applications, there are three important bottlenecks that can be found in the exploitation of historical flights: (i) the loss or absence of metadata related to camera calibration; (ii) the availability of control points, mandatory for georeferencing the images in the same coordinate system and performing comparison through time; and (iii) the lack of a robust automatic method that allows the generation and comparison of models without user intervention.
In our case, the implemented method allows us to deal with the automatic orientation and self-calibration of the historical images without having metadata related to the calibration parameters or the orientation of the images, using to this end an open source solution [33]. In addition to this, the results obtained can be statistically assessed in accordance with georeferencing errors, which is not easy considering the problems inherent to historical flights. Last but not least, the method developed can be extrapolated to other historical flights and thus be used to derive 3D metric information (elevation) from 2D historical images.
In any case, there is still a lot of work to be done in relation to the conservation, publication, and valorisation of historical aerial images. Although there are methods and guidelines in the geo-community and even methods that automate the geoprocessing and 3D exploitation of this data, such as those presented in this article, these types of studies still require human intervention since georeferencing remains necessary. This can be considered the main weakness of this type of analysis. Regarding future perspectives, it could be interesting to incorporate automatic learning methods (e.g., Machine Learning or Deep Learning) capable of automatically extracting semantic entities that allow progress in the georeferencing of images and in the detection of changes. Last but not least, the method presented in this paper raises many questions about governmental policies, from analogue image distribution policies and their related products to the land use and cultural changes that occur. For a more complete study of urban growth, each of these factors should be incorporated in future research, hopefully aiding in the improvement of disciplines such as urban planning, history, and geography.