CORRECTION OF SYSTEMATIC RADIOMETRIC INHOMOGENEITY IN SCANNED AERIAL CAMPAIGNS USING PRINCIPAL COMPONENT ANALYSIS

Orthophotomosaic is defined as a single image that can be layered on a map. The term “mosaic” implies that it is produced from a set of images, usually aerial images. Even if these images are taken during cloudless period, they are impaired by radiometric inhomogeneity mostly due to atmospheric phenomena, like hotspot, haze or high altitude clouds shadows as well as imaging device systematisms, like lens vignetting. These create some unsightly radiometric inhomogeneity in the orthophotomosaic that could be corrected by using a Wallis filter. Yet this solution leads to a significant loss of contrast at small scales. This work introduces an alternative to Wallis filter by considering some systematic radiometric behaviours in the images through a principal component analysis process. Figure 1a: Orthophotomosaic with overlaps


Context
Working on former aerial campaigns has become a major issue in remote sensing and, more generally, in data sciences like for examples the European project Time Machine (Time Machine, 2019) or more local project as Swisstopo historical scanned maps (SwissTopo, 2018) and French National Mapping Agency scanned aerial campaigns available on the website IGN Remonter le Temps (IGN, 2016).
These data are interesting as they are a unique mean to access to past information about land cover. have generally be scanned by national mapping agencies and are thus available as digital data compatible with automatic processing. Thus, they open new research opportunity, for example in the studies of time varying phenomena. Yet, their exploitation leads to several geometric and radiometric processing issues. This work will specifically focus on the radiometric issues in scanned analogue airborne images.

Dataset
Datasets are only composed of images. Unlike recent aerial campaigns, no information about their exposure time nor the camera * Corresponding author. calibration is provided with images, yet it is here assumed that they were shot with the same camera and following the same acquisition pattern (Abdullah et al., 2013) which consists on several stripes, with overlaps to ensure stereoscopic analysis. Geometric processing of these images have been done using MicMac photogrammetric software (ENSG, 2016) (Rupnik et al., 2017): images are oriented and orthorectified, making it possible to generate raw orthophotomosaic renderings as shown on Figure 1.
There are many ways of mosaicking per image orthophotos, some using more advanced blending methods than others. Yet, in order to focus this work on image radiometric corrections, it was decided to consider only two very simple rendering solutions: using the overlaps (Figure 1a) where each pixel is the mean of the different contributing images or choosing the contribution of only one image per pixel, usually the one which the centre is the closest to the considered pixel, which leads to a Voronoi cells partition (Figure 1b).
Both rendering methods have their pros and cons: the radiometry of Figure 1a seems more homogeneous yet small shift between the different images could be observed at full resolution whereas the inhomogeneity of Figure 1b is more pronounced. In any case, both rendering are presenting radiometric heterogeneity requiring a correction step.
ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume V-2-2020, 2020 XXIV ISPRS Congress (2020 edition)  Figure 2, one can notice that the contrast perception is not the same at each scale: small scales suffer from a loss of contrast whereas large scales (full resolution zoom) do not present any visible contrast alteration.
Some more recent works using Contourlet transform (Li et al., 2016) are also returning unsatisfying results when considered at small scales.
1.3.2 Multi-scales approaches: In order to perform correction at different scales, one could consider multi-resolution blending (Ogden et al., 1985) mostly used in panoramic image rendering (Brown and Lowe, 2007) yet requiring a preliminary vignetting (image radial darkening) correction. This requirement focuses on the main issue that is here dealt with: the vignetting parameters are unknown in the present case.

Physical approaches:
Vignetting is mostly due to two factors: the constant lens vignetting and the hotspot which could be determined by the Sun position. Knowing these parameters and using a physical model (Chandelier and Martinoty, 2009) could usually provide some relevant image corrections. Unfortunately, the Sun position, given by the date and the time of the aerial shot should be determined for each image, in addition of other parameters (like the aerosol composition of the atmosphere or the film response curve) that cannot be obtained on old scanned airborne missions. Consequently, physically based radiometric correction cannot be performed in our case.

PROPOSED APPROACH
The proposed approach is purely image based and could be come down to a Wallis filtering method regularized by all the images of a mission, the parameters of the filter being obtained by principal component analysis (PCA), also known as Karhunen-Loève transform (KLT).

Wallis filtering
The principle is illustrated by the Figure 3. There are only three parameters: the size of the window, the desired mean (e.g. the mean value 128, in the case of an 8bit image) and the desired standard deviation (e.g. a value taken between 32 or 64). The mean and the standard deviation would respectively influence the luminosity and the contrast of the resulting image.
The influence of the window size is shown in Figure 4. A small window will remove the vignetting effect and all the low frequencies of the image whereas a large window will keep the image vignetting uncorrected. The size of the window is chosen empirically to obtain an acceptable compromise between vignetting correction and image legibility. In other words, this method is equivalent to calculate a gain and an offset for each pixel according to it neighbourhood.
Yet, even after having chosen manually an optimal neighbourhood window, the final orthophotomosaic rendering is not satisfying at small scales as shown previously on Figure 2 where large land structures are not clearly visible.

Karhunen-Loève transform
This method, based on principal component analysis, is illustrated by Figure 5. A set of vector of the same length, in the present case images represented on a canonical basis are processed to generate a new basis, called KLT basis, where most of the information is contained in the first vectors while the other vectors would give more specific information. In another way, the first vectors will show the systematic behaviours of the vector set. In the case of images, they correspond to common global behaviours shared by all the images.
In this case the input vectors would be of two kind: the local mean of each aerial images and their local standard deviation (for a given Wallis filtering window chosen empirically). For each image a new local mean and a new local standard deviation are ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume V-2-2020, 2020 XXIV ISPRS Congress (2020 edition) determined using only the first vectors of the KLT basis. Using all the vectors from the KLT basis to reconstruct the local means and standard deviation would obviously lead to the initial Wallis filtering shown in Figure 2.
Principal component analysis acts here, to some extent, as a low pass filter: the first axes of the KLT basis are displaying the common information somehow corresponding to low frequencies or small scales behaviour of the dataset.

Image preparation
To obtain a KLT basis, the vectors, derived from the aerial orthorectified images, should have the same length, which is not the case with the "raw" orthorectified images from different size.
In order to give all the images the same size without losing information, instead of cropping them, it was decided to include them in a bigger square frame completed by nearest neighbourhood extrapolation. Each new resized image is given a new mask and a new orientation file taking into account the transformation. Some examples are given in Figure 6. The extrapolation is justified by the fact that using directly the images with their mask would follow some irrelevant KLT basis illustrated in Appendix A.
As the phenomena that have to be corrected are related to very low frequencies, in other word very small scales, working on low resolution images is justified. In this example, considered images size is 800×800pixels.

RESULTS
The data set is composed of 143 orthorectified images (1.5 km footprint each) shot with an analogue camera over the Larzac area (France) in 1996. Computing the KLT basis of this set leads to 143 "eigen" images (KLT axis) associated to respective eigenvalues. Only the ones with stronger eigenvalues will be considered.

Choosing the relevant number of KLT axis
The choice of the "threshold" will be empirical (yet, justified by some observations). Let's consider the correlation of each image with the vector axis of the KLT basis of which the main axes are shown in Figure 8.
The results are shown in Figures 7. All the images are strongly correlated with the first axis: this is the main behaviour of our images set. The second axis is more interesting: two distinct behaviours are observed on Figure 7b between the North and the South of the mission. After considering the eigen image relative to the second axis, this could be interpreted as a change of the ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume V-2-2020, 2020 XXIV ISPRS Congress (2020 edition) hotspot shape, probably due to a different position of the Sun.
Further investigation showed that the images were not taken the same day; it was a two days mission. The correlation with the third or fourth axis (Figures 7c and 7d) is more difficult to interpret: the variation seems somehow more or less correlated along one strip. The correlations on other axis do no show any structured behaviour. So only the three (or four) first axis seem relevant to model the local means and standard deviations of each image, the other orders axis being not correlated with systematic phenomena that we try to remove.

Orthophotomosaic rendering
In order to compare and discuss the results in a more visual and efficient way by artificially increasing the dynamic of the images, the Figure 9 proposes a colour scale representation of the resulting orthophotomosaics rather than a greyscale one.
Compared to uncorrected overlap rendered orthophotomosaic (Figure 9a), obtained results are more homogeneous even if some little discontinuities seem to be present in their Voronoi cells rendering shown on Figure 9c (that could be compared with the uncorrected version given by Figure 1b). These discontinuities are almost invisible when considering overlaps rendering (Figure 9d).
A compromise is obtained by merging the two reconstructions: in practice, applying another Wallis filter on the corrected images where the desired mean and standard deviation are no more fixed but given by the overlaps rendered orthophotomosaic: the details coming from the Voronoi cells rendering and the low frequencies from the overlaps rendering.
Eventually, the proposed method avoids the loss of contrast following the Wallis filtering approach (Figure 9b). Yet, some initial choices could be discussed and developed.

Discussion and further developments
This method is justified by the fact that the image set is composed of a large amount of images, yet on small sets of images, this method would probably show some limits.
Also the fact orthorectified images are used instead of considering directly the camera geometry images might be justified by the assumption that oriented images would give better results, yet restarting the experiment directly on the "original" scanned images before orthorectification could be an interesting alternative approach to this method. The image preparation illustrated by Figure 6 would also be much faster.

PERSPECTIVES
Compared to existing radiometric equalization methods, without any context information, the proposed simple approach gives better visual results. Yet, this method is only limited to panchromatic images. Applying it directly to colour images (separately on red, green and blue bands) leads to some drab and unsightly results shown in Appendix B. The colour approach should be considered in a different way, typically by first avoiding using fixed global mean and standard deviation as it is done in this work.
Another perspective would be to perform the method directly on "raw" images (e.g. before orthorectification) and see its influence on the photogrammetric process leading to final product as digital terrain model and orthorectified images. Eventually, the last step would consist on testing our method on larger datasets, like for example a region or a whole country in order to discuss it possible limits.

APPENDIX B
From top to bottom, colour orthophotomosaic rendered: without correction, with Wallis filtering on the RGB layer, with our method applied directly on the RGB layers: ... This contribution has been peer-reviewed. The double-blind peer-review was conducted on the basis of the full paper.