ARCHIVE AND WARTIME AERIAL PHOTOGRAPHS AND PROCEDURES OF THEIR TREATMENT

The article presents with the use of archive aerial photographs. The first task was to search and identify drainage detail from archive aerial photographs. The second task is to create procedures for processing aerial reconnaissance images (from WWII) to identify sites with potential pyrotechnic load. Both of these tasks are connecting by the effort to determine the internal orientation parameters of the cameras for using and easier calculation of exterior parameters by image correlation. Complete automation process searching of fiducial mark (FM) identification was implemented. The coordinates of all FM are calculated automatically from archive aerial photographs. In addition, the edges of the photographs are automatically found and a program was created to minimize of the cropping of the archive aerial photographs. The next part of the paper describes the procedures of averaging the values of the relative position of FM and transforming archive aerial photographs to a uniform dimension from a set of images taken with the same camera. The second part of the paper describes the process of creating a historical ortophoto with the standard calculation of bundle adjustment performed by an external process in the background of the OrthoEngine module using the Celery library installed as a python service. Finding of external image orientation parameters through bundle adjustment calculation are parameters, in the first, defined in the local system and then transformed into the national geodetic system of the Czech Republic. This entire section is available and free to use for on the internet. The third part of the article describes the practical procedure of the interpretation of archive and wartime photographs with aim of identification of the drainage detail and the procedures leading to the interpretation, identification, location and calculation of the position of unexploded air ammunition. * Corresponding author


Objective
Between April 2019 and December 2021, the Department of Geodesy and Mine Surveying with the Faculty of Mining and Geology at VSB-Technical University of Ostrava https://www.hgf.vsb.cz/544/en and the company Primis spol. s r.o. from Brno http://www.primis.cz/index.php/en are jointly investigating a project supported by the Ministry of the Interior of the Czech Republic, called "Finding unexploded aerial ammunition of World War II (PID-VI3VS/778)". Since the declassification of the contents of the National archive of aerial photos in 1993 with the Military Geographic and Hydrometeorological Institute (formerly Military Topographic Institute in Dobruška) http://www.mapy.army.cz/, the use and processing of archival photos in the Czech Republic have been of interest with many independent investigators, geodetic companies, legal offices and research institutions. Between 2003 and 2005 two of the co-authors participated in the creation of a historical orthophotomap of the Czech Republic made from photos taken for the purposes of topographic maps production (scale 1:25 000). The historical orthophotomap was made from approximately 21 000 aerial photos taken between 1952 and 1957. To date, the orthophotomap has been the most complex piece of work made from archival aerial photos (AAP) of the whole Czech Republic (73 000 kmsq). The historical orthophotomap is available to public at https://geoportal.gov.cz/web/guest/map (layers -Orthophoto map (50´s)). Currently, it is primarily the company PRIMIS in the Czech Republic who uses exact procedures of AAP processing using digital aerial photogrammetry. Within their applied research activities, PRIMIS are developing optimised procedures for the processing of AAP. The company also process photos from World War II and they closely cooperate with the major processor of war photos in Europe, the company Luftbilddatenbank Dr. Carls GmbH https://www.luftbilddatenbank.de/main/index.php?webcode=ho me&language=english . The companies mainly cooperate on the research in the interpretation of unexploded ammunition based on orthophotomaps made from war photos. Within ISPRS literature, recently a number of works have been published on the topic (Jae Sung Kim et al, 2010;Patrik Meixner et al, 2016;Chen et al, 2016 ). P. Meixner et al, 2016 have been written in cooperation with the above stated German company. Nowadays, Luftbilddatenbank Dr. Carls has approximately 106 000 photos of the Czech territory back from 1939 to 1946.

Solution procedure
Through gradual steps, the authors made a semi-automatic system of AAP in order to provide a practical self-service system to make high-quality and position precise (RMSE xy =2.4m) orthophotomaps to identify unexploded ammunition. The task the company is dealing with is the identification of the drainage details in the AAP.

Concept of solution procedure
The solution comprises several interlinked programmes. In the first step, searching for the drainage detail, the users select suitable photos for processing based on the access to the public AAP portal of the Czech State Administration of Land Surveying and Cadastre (CUZK) https://geoportal.cuzk.cz/(S(m4hjwcxvau1lkhzitbwwcuwq))/Def ault.aspx?lng=EN&head_tab=sekce-00 gp&mode=TextMeta& text=uvod_uvod&menu=01&news=yes&UvodniStrana=yes. The photos are selected by overlaying the area of land drain circumference and the circumference of photos overlaying the selected land drain. A simple condition applies for the choice of the photos from the photo database: the date of photo must be one year + following the recording of completion of drainage construction (both the dates are included in the databases). In case of searching for unexploded ammunition, based on the agreement with Luftbilddatenbank Dr. Carls GmbH, timerelevant photos of the locality of interest are selected. Usually, it is necessary to select three time periods of the given locality (1 before an air attack and two after an air attack). The following procedures for both tasks may be identical.

Partial tasks for recreating AAP
The further steps are significantly automated work procedures of AAP processing and their conversion to be subjected to bundle adjustment and matching calculations to produce the final orthophoto to be able to interpret the contained situation. AAP are scanned using special photogrammetric scanners with track-in internal accuracy of the identical spot on an image of approx. 2 micrometres, i.e. one seventh of a common scanning element size (14 micrometres per TD1). Images are placed in the scanner per discrete shots as they have been preserved in the archives of aerial photos in the Czech Republic (or war photos in reels). This implies that despite careful placement of images in the scanner, the beginning of scanning of AAP is always general towards the fiducial marks (FM). The programme developed by authors to unify and convert the images automatically searches for the fiducial marks on the discrete AAP and modifies the images to prepare them for matching and further digital processing. The programme is conceived as noinstaller one. To run the programme, only a Microsoft Visual C++ library is required, which makes part of the Microsoft package. The programme supports three types of fiducial marks of AAP. Other fiducial marks can be added into the library any time. The authors adhered to all common photogrammetric standards when converting the AAP. The steps are as follows: 1. Constructing the AAP fiducial mark library contents 2. Distribution of FM in AAP 3. Automated identification of FM in the AAP images 4. Determination of FM coordinates in image pixels 5. Calculation of FM values and their averaging towards all images in a set 6.Identification of analogue image margins and calculation of image size cropping 7.Recalculation of image size and conversion of the original AAP 8.Storing the final converted image data for further processing or third-party software's processing.

Constructing the AAP fiducial mark library contents:
A digital AAP is obtained via scanning the real analogue photo using a photogrammetric scanner. AAP is usually scanned margin to margin. The geometrical parameters of AAP are ensured by FM. Figure 1 below shows three types of fiducial marks that were selected for the FM library with respect to their high frequency.

Automated identification of FM in the AAP images:
The automatic identification of FM image in AAP is based on the feature surroundings evaluation methods in an image and its specification within FM surroundings. First, according to FMs and in the AAP image contents, key features and their descriptors are calculated using SURF (Speeded-Up Robust Features) method. Next, correct pairs of features are searched for using RANSAC (RANdomSAmpleConsen) algorithm. The last step is to calculate homography and the position of the detected mark is determined. Figure 4 shows the final result of the calculation and the identification of the surroundings based on identifying the mark in the FM library.

Determination of FM coordinates in image pixels:
Based on the identification of FM and the geometrization of its internal part (the inner circle of the mark in Figure 4), we determine the focal point of the mark and calculate the inner part pixel in image pixel coordinatessee Figure 5. All FMs of each AAP are stored in a temporary file for further calculations.

Calculation of FM values and their averaging towards all images in a set:
Having calculated all FMs of all AAP from a set of photos, the calculated values are averaged, and a systematic repositioning is determined, which eliminates the difference from the discrete FM median. A "basic image size matrix" is determined in the temporary files.
2.1.6 Identification of analogue image margins and calculation of image size cropping: AAP images that are supposed to enter the matching calculations to produce an orthophoto of the territory must be "unmasked", i.e. the parts with the fiducial marks or frames are cropped. Automatic steps are sequenced to find the frames, and crop these non-image parts so that they showed only the image.

Recalculation of image size and conversion of the original AAP:
Having averaged the size of AAP croppings to get rid of frame data, FM locations in the images and own frames, we recalculated the image sizes towards FMs, or their average positions calculated as per step 2.1.4. This procedure ensures that the recalculation of the scanned AAP and conversion of AAP into formats suitable for processing in the next step is correct with respect to the best available techniques. This way, the image set still maintains an identical focal distance and image metrics (with scanning element of 14 micrometres) is maintained. The principal feature of the image as well as the symmetry point are identical in all the aerial photos of the given set, and the dimensions of all images in pixels will be identical for all images that enter all further calculations.
2.1.8 Storing the final converted image data for further processing or third-party software's processing: The converted AAP are stored automatically in the final converted image directory.

Automatic calculation of bundle adjustment
The aim of creating the module for triangulation and inlaying of images was to provide users with a simple web interface, where they upload automatically transformed photos into detected fiducial marks along with manually measured reference points. The chain of applications analyses the images and automatically creates the final orthophoto, without any user's interference. This way, the user does not need any software, only a web browser and Internet connection.
Processing phases: 1. Detection and calculation of key features 2. Matching the key featuresdetermination of relative orientation 3. Calculation of incremental bundle adjustment 4. Transformation of the image bundle into S-JTSK coordinate system 5. Orthogonalisation and inlaying into the final orthophoto Currently, apart from commercial applications, the issue of automatic processing of disorganised sets of images has been approached via a number of open source projects. For the first three phases, OpenMVG library (Moulon, 2016) was used, which also integrates an open source Ceres library http://ceressolver.org/ to calculate bundle adjustment.

Detection and calculation of key features:
The key features ( Figure 6) unambiguously characterize the image area so that the area could be found and compared with the identical area on another image. To detect and compare key features in an image, the OpenMVG library uses a SIFT (Scale Invariant Feature Transform) detector in the default setting (Lowe, D.G., 2004). Unlike a simple correlation of two areas in photos, this detector is partially invariant towards changes in the view geometry, i.e. rotation (approx. 15 degrees), changes in the scale and noise. If key features are detected in each image, including the descriptors, it is possible to proceed to their matching and identification of the corresponding pairsthe correspondences arising from the projection of a point in 3D space into each image, which shall have very similar descriptors. The level of agreement of two key features is clearly defined on the grounds of Euclidian distance of their SIFT descriptors. In case the collection of images is not prearranged and relationships among the different images are unknown, all images must be compared one-to-one. The corresponding sets of key features obtained via matching are usually burdened by an error and false correspondences that are caused by changes in the camera position, in the lighting, digital image noise, etc. These false correspondences may be eliminated using a geometrical criterionepipolar conditionsee Figure 7.

Figure 7. Simple representation of epipolar coordinates
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition) Point X in the 3D space forms an epipolar plane along with the projection centres C and C'. Thanks to the line of the epipolar plane and the projection plane, epipolar lines (epipolars) are formed, which pass through the points x and x' being the projections of point X into the projection plane. These lines crosscut the epipoles e and e', where the epipole is the projection of the projection centre of one camera into the second camera's projection plane. The algebraic formulation of the epipolar condition is the formula below: Where F is the fundamental matrix of 3 x 3 size and the ranks 2 define the relative relation between two cameras independent of the scene structure. To calculate the fundamental matrix, it is not always necessary to know the parameters of the discrete cameras' internal orientation. In the OpenMVG library (Moulon, P., 2016) there is a RANSAC algorithm (Fischler, M 1981) to select the key features correspondences complying with the epipolar condition. The implemented algorithm allows for iterative identification of the best solution matching the given model, i.e. the equation x' F x = 0, and exclude wrongly detected correspondences. The fundamental matrix is calculated using a 7/8-point algorithm (Hartley, R, 2003).

Calculation of incremental bundle adjustment:
Based on the calculated relative orientations for the discrete images, it is possible to determine the approximate spatial data corresponding to the image coordinates of the detected correspondences that are used as an estimate of the input parameters entering the complex bundle adjustment. The aim of bundle adjustment is to find optimal parameters of internal and external orientation, including the radial distortion coefficients of the lens and such spatial data, for which the distance between the point projections in the space into the image and their detected image coordinates is minimised. Adjustment is attributed to points in the 3-D space and thus the parameters of external and internal orientation. The partial derivation values may be determined analytically via a function derivation according to the variables as well as numerically. Due to significant complexity of the functional relationships to calculate the image coordinates, where the analytical development is difficult, the Ceres library http://ceressolver.org/ uses a numerical solution. At the start of the computation, the most suitable pair of images is chosen, e.g. based on the number of detected key features. The projection centre of the first image in the pair defines the start of local coordinate system, and the external orientation rotation matrix of the first image is selected as the unit matrix. Each image from the pair contains the key features also detected in other images. Thanks to these correspondences, other images are "connected" into the local coordinate system. Bundle adjustment is carried out after each iterationsee Figure 8 for the result.

Transformation of image bundles into the Czech coordinate system called S-JTSK:
The last phase in the determination of external orientation parameters is the transformation of the geodetic coordinate system via 3D similarity transformation using reference points. Because the image coordinates are always measured in two and more images, it is possible to determine 3D coordinates in the relative system of coordinates. These coordinates are used along with the coordinates of the reference control points in the geodetic system to calculate the parameters of 3D similarity transformation. Using a given transformation key, the 3D point coordinates are transformed, including the calculated parameters of external orientation into the geodetic system. In the experiments we obtained a higher number of consistent results (image correspondences in the final orthophoto) during bundle adjustment in the local coordinate system and the subsequent 3D transformation than from the final bundle adjustment with the reference points of S-JTSK system. We assume this was caused by certain internal stiffness of the image bundle and discrepancies between the terrain model used for reference point level interpolation (carried out from the current DTM) and the reality captured in the archival aerial images decades ago. The bundle adjustment described in the previous Subsection (2.2.3) is solved in the local system. For the transformations into the S-JTSK geodetic system it is vital to find Ground Control Points (GCP). GCP for photogrammetric calculations and creation of an orthophoto may be selected manually or set up through a required text file using the application http://www.vugtk.cz/euradin/gcp/ .The manual setup of a GCP file is based on the selection of GCP from the (CUZK -https://geoportal.cuzk.cz/) data, either from the cadastral map or from a ready orthophoto of Ground sample distance (GSD) in the size at least identical to the size of the scanned AAP's GSD or using the trigonometric point (TP) coordinates of the CR. The last option is theoretically the most accurate but our capacities to identify TP in historical images is very limited, except for the triangulation towers and point signalisation of IV. and V. order, as well as through the required density of the reference points to produce an orthophoto, the usual identifiable number of TP is not sufficient. To select reference points, it is possible to use points in the archival images with points in the cadastre, e.g. in Figure 9 (left) and its detail in Figure 9 (right) from the situation in AAP. Figure 10a is an identical situation from the current imaging and Figure 10b shows a cut-out from a cadastral map. The coordinates are obtained reading them from a digital cadastral map. However, on the territory of the Czech Republic it is currently more advisable to find identical points in AAP and the current orthophoto as, in general, the quality of the orthophotomap in the urban areas is higher than RMSE xy of the cadastral map. Therefore, we recommend users to read coordinates of GCP from produced orthophotos, particularly in such points, in which we may assume that they have not changed their position since the AAP was taken. In case of doubt it is recommended to measure the reference points geodetically in the terrain from the points interpreted from AAP. The automated procedure is described below. As mentioned above, bundle adjustment is first carried out in the local system. For the transformation into the national geodetic system (S-JTSK), the GCP in the national system are needed.
Having read the GCP coordinates, the application was cloned from repository https://github.com/posm/posm-gcpi, where the sources for the underlayers and the output coordinate system were modified. The images are uploaded into the application only locally in the browser, they are not sent anywhere. The measurement outcome is a text file containing the measured image coordinates in the image coordinate system, and their 2D equivalents in the S-JTSK coordinate system. Such a file is sent by the user to authors' server, where a digital model of the Czech Republic is stored, and which is used to interpolate the Z coordinate of the selected GCP. Each GCP must be measured in at least two images to be able to calculate its 3D coordinate and the triangulation in the local coordinate system, in which the initial bundle adjustment is carried out. Figure 11 shows the application to measure the reference points; on the left there are the aerial photos, and on the right, there are the measured 2D reference points in the geodetic system of the tested locality. Figure 11 Web application to measure coordinate reference points

Orthogonalisation and inlaying into the final orthophoto:
The final step of AAP processing, after the conversion of images, bundle adjustment and transformation into S-JTSK, is orthogonalization and inlaying of the images into the final orthophoto. For the orthogonalization, we use the digital terrain model of the Czech Republic of a regular step of 20m. Another alternative for the future is the calculation of DMT via the correlation of AAP used for the production of orthophotos. One of the authors' aims was also to create a module for automatic inlaying. As photographic aerial surveying is usually carried out with sufficient coverage, it is clear that each place on the orthophotomap is captured on two or more images. Thus, it is necessary to define a function to clearly and optimally choose such parts of the images that are least burdened by lens' optical errors. In the inlaying module there is an implemented function, where an image is chosen for each pixel of the final orthophotomap, where the distance of the 3D coordinate projection from the optical axis in the focal plane y is minimal, using the projection of its 3D coordinate in the geodetic system into the plane of the image. For the purposes of visualisation, the discrete images may be replaced by grey levels and represent only parts used in the final orthophotomapsee Figure 12. The choice operation for the most suitable image for each pixel is automatic during the programme run. This way, there is no need for orthorectification and storing of all full images. To maintain the image quality, the R, G, B values of the different colour channels are determined by bilinear transformation from the nearest pixel's surroundings in the source image. For the OrthoEngine module there is also a web interface, where the input is the digital terrain model, the parameters of external and internal orientation in the native format of OpenMVG library, and 2D coordinates of reference points. Figure 12 Mask used for inlaying images into the final orthophotomap

CONCLUSIONS
The aim to produce a historical orthophoto is to allow for the interpretations and display of conditions at the time of image exposure. This way, it is possible to interpret unexploded ammunition, or the structures of the drainage system, which are difficult to identify in the contemporary images. The final interpretations and vectorisations of the objects in demand based on the produced historical orthophoto are very valuable due to the data on the interpreted object position. The orthophotos have been produced using the described exact methods, and verifiable and controlled procedures. According to copyright rules, the final data interpretation and the orthophoto itself are a cartographic representation of original research results obtained by authors. However, the identification of objects on historical orthophotos does not end by simple production of the orthophoto, but the accuracy of the whole system of the historical orthophoto production is verified by digging out the drainage system piping or the unexploded bombs back from World War II. Figure 13 shows an example of eight processed images from the testing locality having applied the procedure described herein. Figure 14 shows the final processing of wartime images.