IMPROVING LOCAL ADAPTIVE FILTERING METHOD EMPLOYED IN RADIOMETRIC CORRECTION OF ANALOGUE AIRBORNE CAMPAIGNS

An orthophotomosaic is as a single image that can be layered on a map. It is produced from a set of aerial images impaired by radiometric inhomogeneity mostly due to atmospheric phenomena, like hotspot, haze or high altitude clouds shadows as well as the camera itself, like lens vignetting. These create some unsightly radiometric inhomogeneity in the mosaic that could be corrected by using a local adaptive filter, also named Wallis filter. Yet this solution leads to a significant loss of contrast at small scales. This current work introduces two elementary studies. In a first time, in order to quantify the loss of contrast due to the use of Wallis filter, a simple multi-scale score is proposed based on mathematical morphology operations. In a second time, an optimal window size for the filter is identified by considering some systematic radiometric behaviours in the images forming the mosaic through Principal Component Analysis (PCA). These two elementary studies are preliminary steps leading to a method of radiometric correction combining Wallis filtering and PCA.


CONTEXT
The data considered in the present work are a set of 143 orthoimages that have been computed, from scanned analogue airborne shots, with the open-source software MicMac (Rupnik et al., 2017). The resulting mosaic is displayed on Figure 1 and presents some radiometric inhomogeneities that seem to follow a similar pattern confirmed by considering the mean of the images as shown on Figure 2. This pattern is due to a combination of several factors: some of them are linked to external parameters, like the hotspot (Teng et al., 1997), a sum of retro-specular phenomena caused by the way aerosol and ground are interacting with sunbeams, and some to internal parameters, like the lens vignetting (Dehos et al., 2012).
Knowing these external and internal parameters for each image of the airborne campaign allows to perform a physical rigorous radiometric correction (Chandelier and Martinoty, 2009). Yet, these parameters are not always available. It is unfortunately the case of the images constituting the data-set of the current study * Corresponding author. where the lens vignetting or the Sun position are unknown. In this case, it is not possible to rely on physical models in order to perform radiometric correction. Another method (Wallis, 1976) is then applied consisting of giving locally the same mean and the same standard-deviation, which is the definition of Wallis filtering: (1) where µw(x, y) and σw(x, y) are respectively the mean and the standard-deviation computed on the pixel values of the image Img contained in a w×w pixels square centred on (x, y). µ0 and σ0 are respectively the desired mean and standard-deviation.
In practice, https://remonterletemps.ign.fr (IGN, 2016), a governmental French website releasing, among other geographical data, orthophotomosaic of the whole country from different dates in order to show evolution of the lands, displays both recent (digital) orthophotomosaics, where the radiometry is corrected using physical method, and old (analogue) ones corrected by Wallis filtering. A screenshot of the website shown on Figure 3 illustrates the loss of contrast caused by Wallis filtering at small scales. The next part of the work introduces a multi-scale score differentiating these two behaviours.

MULTI-SCALE ANALYSIS
The experiment is quite basic: the same location is considered on both orthophotomosaics, at different scales (7 in the following example). The 14 level of greys orthoimages (2 dates at 7 scales) are thresholded using the median value of each image in order to obtain binary images with half white and half black pixels. The score for an image, i.e. a correction at a given scale, is define as the proportion of pixel, in the binary image, that are invariant by opening and closing in the sense of mathematical morphology. In the presented case, the structuring element (s.e.) is a square window of 3, 5 and 7 pixels displayed respectively in the layer red, green, and blue of the invariant map shown as part of Figure 4. Pixels that are invariant with a s.e. of size 2n + 3 are also invariant with a s.e. of size 2n + 1: the yellow pixels of the maps are corresponding to pixels invariant with 5 × 5 s.e. and consequently 3 × 3 s.e. too but not 7 × 7 s.e. that are represented in white.
The curves or Figure 4 seems to suggest that the multi-scale score defined here is decreasing for small scales in the case of a Wallis filtering radiometric correction whereas it stays almost constant for physical corrections. To confirm this observation, this multi-scale score is performed on the 143 orthoimages of our data-set, actually converted in level of grays with the formula max(R, G, B), and their mean multi-scale scores are displayed for each Wallis filtering size in the Figure 5.
The influence of the Wallis filter size is double: it has an influence on small scales, where the large binary area of same value are vanishing but also on big scales where the proportion of invariant pixels is decreasing while the filtering size gets smaller. Eventually, one should insist on the fact that the choice of µ0 and σ0 in Wallis filtering has no influence on the resulting binary image (the threshold is the median of the image). Then the score only depends on the filter size.
A small filter size results into a degradation of the image. This leads to the following question: which is the largest Wallis filtering size that removes hotspot and is the resulting correction acceptable at small scales? Wallis filtering. In the middle, their respective multi-scale scores: in red and brighter, proportion of pixels that are invariant by opening and closing with a 3×3 pixels structuring element, in yellow and brighter, the same with 5×5 and in white, with 7×7.

PRINCIPAL COMPONENT ANALYSIS
To bring an answer to the previous question, two concepts are introduced: Principal Component Analysis (PCA) and Moran's I. PCA focuses on the systematic behaviours of images (typic- The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France ally, the hotspot) and Moran's I measures the spatial autocorrelation of a phenomenon, or in other words, the behaviour of a variable compared to its neighbours. The PCA, also known as Karhunen-Loève transform (Turk and Pentland, 1991), and referred to KLT in the rest of this document, is illustrated by the Figure 6: a set of images, considered as vectors of same dimension, are decomposed into eigen images (associated to eigen values) that form a basis. The first vectors of the KLT basis are related to the main behaviours of the image set. In order to apply this transform to our data-set, each image is cropped to have the same size. More specifically, they have been sub-sampled from a factor 15 then cropped to form 600 × 600 pixels images. Yet an issue appears after performing the KLT: the eigen images ( Figure 7) are focusing on the image border masks more than on the image content. In order to avoid this, an extrapolation of image values (nearest neighbourhoods) is performed on the border masks. The resulting eigen vectors ( Figure 8) are no more disturbed by the images border masks. Each image can be reconstructed as a linear combination of the eigen images forming the KLT basis. By construction, the first KLT axes contain most information of images and each eigen image could be interpreted as a pattern shared by several images. By considering the images coefficients for the four first KLT axes as displayed on Figure 9, some spatial correlation seems to occur. The first KLT axes may be interpreted as a systematic and spatially autocorrelated phenomenon (e.g. the hotspot). The following assertion could be made: removing the hotspot with Wallis filtering would have an impact on the KLT reconstruction coefficients by reducing their spatial autocorrelation. Here, the coefficients u n i are designing the projection of the images Imgi on the KLT axes n. By definition, the sum of the KLT axes (i.e. the eigen images) respectively multiplied by the coefficients (u n i )n=1...N returns the image Imgi.
In order to quantify the spatial autocorrelation of the coeffi- where N is the number of images, 143 in this work, ui the coefficient related to image i,ū the mean value of the coefficients and pij the proportion if pixels belonging to i in j (by convention pii = 0). Moran's I could be interpreted the following way: if I is close to +1, the spatial correlation is perfect, if I is close to −1, the spatial dispersion is perfect and if I is close to 0, the spatial model is random. Figure 10. Evolution of Moran's I for different Wallis filter size. Figure 10 shows the impact of Wallis filter size on spatial correlation of the coefficients. A Wallis filtering with a windows The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France size of 100% of the image is a global contrast and luminosity applied to the whole image. In this case, only the three first KLT axes are showing spatial correlation and may be related to hotspot and vignetting. This spatial correlation seems to vanish for Wallis filter size smaller than 9% of the image size. Yet, as seen previously on Figure5, this filter size leads to some degradation of the image, more specifically on small scales. To confirm this observation, the resulting orthophotomosaic is displayed on Figure 11: the Wallis filtering of size w = 9% is done on each colour channel of the image, where µ0 is the mean of the channel and σ0 its standard-deviation. Figure 11. Orthophotomosaic with basic Wallis filtering correction using optimal window size (9%).
Hopefully, the results seen in the previous paragraphs may lead to some improvement of Wallis filtering developed in the following section.

WALLIS FILTERING IMPROVEMENT
The local mean µw(x, y) and local standard-deviation σw(x, y) defined in equation 1 could actually be considered as "images" of the same size than the image Img(x, y). For example, µw(x, y) can be considered as Img(x, y) convolved with a square mean blur of w×w and σw(x, y) as the square root of the image (Img(x, y) − µw(x, y)) 2 convolved by the same w×w square blur. The computation time of µw(x, y) and σw(x, y) can be significantly shortened by following a method developed by (Phan et al., 2012) and using the concept of integral images of power m defined as: The local mean can then be computed from the first order integral image with this formula: µ2n+1(x, y) = ( Img 1 (x + n, y + n) − Img 1 (x + n, y − n) − Img 1 (x − n, y + n) + Img 1 (x − n, y − n))/(2n + 1) 2 .
For the local standard-deviation, the formula introduced by (Phan et al., 2012) is: .
The computed local means and the local standard-deviations are then forming two sets of N = 143 images each. A KLT is then performed on each set (µ)1...N and (σ)1...N . Theμ k (respectively theσ k ) are defined as the images reconstructed with the k first KLT axes from the set (µ)1...N (respectively the (σ)1...N ). In the previous section, it has been shown (Figure 10) that the radiometric distortion caused by the hotspot and the lenses vignetting is mostly linked to the three first KLT axes. The following assertion will be made:μ 3 andσ 3 are only related to systematic phenomena, like hotspot and vignetting, and, therefore, independent from the ground (forest, fields, urban area, etc.). The new correction introduced in the present work, and inspired by Wallis filtering equation 1, is given by the formula: (6) where µ0 and σ0 are respectively defined as the image global mean and global standard-deviation (for a given colour channel). Figure 12. Orthophotomosaic with improved Wallis filtering correction using optimal window size.
The resulting orthophotomosaic is displayed on Figure 12. Compared to Figure 11, the new orthophotomosaic seems to conserve contrast at small scale. This visual observation is confirmed quantitatively with the multi-scale score of the corrected images (max(R, G, B)) shown on Figure 13. Compared to the images corrected with basic Wallis filtering ( Figure 5: multiscale score related to a window size of 9%) the images corrected with an improved Wallis filtering show better multi-scale behaviour.
Yet, a last issue remains, the images do not show the same colour. These inhomogeneities could be caused by several factors, Figure 13. Multi-scale score of the 143 images corrected with the improved Wallis filtering.
like atmospheric veil or, more likely, some very little variations in the complex chemical process of analogue film developing. In order to remove these colour inhomogeneities, the last section of this work focuses on a simple colour correction method.  In order to make the initial image look like its equivalent image extracted from the mosaic, a histogram transfer is applied:

FINAL EQUALISATION
where CHImg is the cumulative distribution function of the histogram of Img or, in other words, CHImg(v) is proportional to the number of pixel of Img having a value smaller or equal to v. CH −1 Img is defined as the inverse function of CH Img .
Using the histogram transfer method instead of a simpler one, as for example a global Wallis "filter", is justified here because of the nature of the analogue images: their response is not linear like a CCD sensor (Litvinov and Schechner, 2005) so a more complex approach may be required.
After performing this colour equalisation for all the images, the same step is repeated in an iterative way with the Img ′ and their corresponding Img ′ . The final result is reached in a few iterations only (Figure 15) and gives a more aesthetic result, at small scales, than the previous approaches (Figures 11 and 12).

COMPARISON WITH SMALL-SCALE IMAGERY
In this last section, a little experiment is proposed in order to test the limits of the presented method at small scales: a large image crop of our correction result shown at Figure 15 is compared to a small-scale orthoimage of the same area provided by the French website https://planetobserver.com/ (Masselot, 2000). Both images have the same size, our orthophotomosaic has been sub-sampled to fit the small-scale image. In stead of performing a multi-scale score on both images, the score will be computed only for one scale.
The maps displayed on Figure 16 are defined the same way than in Figure 4 and 4: the level of greys images are defined by max(R, G, B). In our image the proportions of invariant pixels for structuring elements of size 3 × 3 (red), 5 × 5 (yellow) and 7 × 7 (white) are respectively of 80%, 62% and 49%. At the same scale, the small-scale image provided by PlanetObserver returns a score triplet of 95%, 85% and 75%. The proportion of invariant pixels is quite lower in the case of the improved Wallis filtering correction which still betrays some loss of contrast at small scale.
Even if our approach shows some real improvement, specifically for small scale visualisation, in comparison to the Wallis filtering that is currently used to correct the radiometric of scanned analogue airborne images it does not reach the level of readability of "satellite" images at very small scales.

PERSPECTIVES
The current work mostly focuses on radiometric correction of orthoimages with the aim of displaying aesthetic orthophotomosaics. Yet some further experiments are scheduled: as applying the same process directly on raw images in their initial camera geometry (interior orientation) or trying to determine the centre of the hotspot in order to estimate the position of the Sun in the sky.
One undiscussed issue in this work is the fact that, on one hand, the varying sizes of the orthoimages and, moreover, their (large) border mask needs to be managed with a time consuming extrapolation and might even lead to some unknown artefacts. One the other hand, raw images, i.e. in their initial camera geometry, are sharing the same size and the same narrow border mask. In addition, the campaign are usually following the same scheme of acquisition, by band, which may be considered as a new systematic phenomenon. Transposing our approach to this kind of images may lead to some interesting results.
The other issue of our approach is that, even if the resulting corrected images seem better than the ones provided by a simple Wallis filtering, the multi-scale behaviour of our image is still far from the quality provided by physically based corrections. Unfortunately, the metadata needed to perform a physical correction are not available for old campaigns. An ambitious perspective would be to modify our approach in order to guess some of the metadata, like the Sun position, the atmospheric veil, or the camera lenses vignetting and, eventually, perform a physically based correction like it is done on more recent campaign.
To sum up, this work could be considered as the first steps of future studies on the radiometry of analogue images provided by historical campaign and archived in some National Mapping Agencies, or any equivalent organisations, since the late nineteenth century.