From research to production: radiometric block adjustment to automate and improve tonal adjustment of orthomosaics created from airborne images

Production of orthophotos and orthomosaics from airborne imaging has been standard procedure of National mapping agencies for decades. Various aspects affect to the colour, or tones of the images: time of the day and year, atmosphere, illumination conditions, view and illumination angle, BRDF-effects (Bi-directional reflectance distribution function), object, sensor, and whole imaging system. Quantitative and automated solution for creating evenly coloured image mosaics is to use method based on radiometric block adjustment, that has similarities to more well-known geometric block adjustment. We have applied radiometric block adjustment method called RadBA, originally developed for drone image blocks collected with hyperspectral sensor, to image block collected with large-format photogrammetric camera. Goals of our work are: 1) to speed up orthophoto deliveries to Finish Food Authority used for EU farming subsidies monitoring, 2) improve the quality of images delivered to Finnish Forest Centre used for forest inventory, and 3) to automate, fasten and improve the current colour balancing process of orthophotos at National Land Survey of Finland. RadBA-based tonal processing with 4-parameter BRDF-correction was able to clearly improve the tonal quality of the image mosaic collected during three different dates. We still need to automate various steps of the RadBA-workflow and improve final steps in converting 16-bits per band mosaics to visually good looking 8-bits per band mosaics. Our long-term goal is to create tonally high quality, seamless orthomosaic of whole Finland.


Introduction 1.1 Orthophoto production in Finland
Production of orthophotos and orthomosaics from airborne imaging has been standard procedure of National mapping agencies for decades.In Finland, the imaging campaign covers whole Finland at 40 cm ground sampling distance (GSD) in three years cycle.The airborne imaging work is governed by KALLIO coalition, with core members National Land Survey of Finland (NLS), Finnish Food Authority (FFA), Finnish Forest Centre (FFC), Metsähallitus (Forest Administration of Finland), and Finnish Environment Institute (Syke).NLS is responsible for the practical organization of the imaging program.Additionally, Ministry of Agriculture and Forestry, Ministry of the Environment and Ministry of Defence are represented in the KALLIO coalition.
Finland is divided to 147 imaging areas, production blocks, with average size of 54 km x 54 km i.e 2900 km 2 (Figure 1).Each block is imaged every three years, except the most northern blocks at Lapland where the imaging cycle is 12 years (shown in brown colour in Figure 1).Each year, NLS performs the imaging of half of the areas, about 20 areas yearly, and rest of the areas are imaged by two imaging contractors based on tendering process.Images are collected with large-format photogrammetric cameras storing information in four bands, red (R), green (G), blue (B) and near-infrared (NIR).
During an imaging season, images are collected from early spring since the snow has melted to late summer until deciduous trees start their autumn foliage.Imaging flights are done only during clear sky conditions, and the requirement for solar elevation angle is 25° above horizon.When the weather is good, optimally one image block can be collected in one flight approximately in one hour.When the weather is more challenging, one image block may require two or more flights during several days.From the acquired images, NLS and contractors produces 8-bits per band RGB and false-colour CIR (NIR, R, G) orthophoto mosaics with 50 cm GSD.Orthophotos and stereo models are used in topographic map updating at NLS, and orthophotos are also published as open data via website and interface services (NLS 2024a(NLS , 2024b)).In addition to NLS, collected images are also delivered to FFA and FFC.FFC uses imagery for collecting forest inventory data for private forest owners.FFC uses both automatic and interactive interpretation methods that are based on raw 16 bit per band images and on radiometrically adjusted CIR forest orthos.FFA operates as the EU's paying agency for farming subsidies and uses CIR orthophotos for monitoring the use of farming subsidies.EU sets strict deadlines for performing the subsidies monitoring each year.
The current orthophoto production workflow at NLS is as follows: image block aerial triangulation is performed with Inpho MatchAT software based on GNSS-IMU data collected during flight.Orthophoto colours are tuned using UltraMap software.This colour balancing part is done by expert operator and is an iterative, partly subjective process.Goals of the colour balancing are to produce visually realistic, good looking RGB orthomosaic, good looking CIR orthomosaic, and to contain as much as possible of the colour information of the image.Colour balancing is performed using 16 bit per band images, but final orthomosaic is produced in 8 bits per band.The amount of colour information of the end product is monitored based on image histogram statistics with in-house software called HISTOQC.It calculates histogram width and number of pixels at 0 and 255 values, i.e. dark pixels and saturation.These statistics where originally developed for quality control of scanned film images (Markelin and Honkavaara, 2004).Finally, the orthomosaics are visually compared to NLS reference mosaics, i.e. previous mosaics that have been judged to be visually excellent.Imaging contractors may have different brands and versions of photogrammetric cameras compared to NLS, and contractors perform orthophoto production with software tools they have.
Even though the quality criteria's and HISTOQC control for imaging contractors are same as for NLS own production, variations in the tonal quality of orthomosaics is unavoidable.
Based on the imaging conditions within one year and between years, due to three different imaging operators with different sensors, and due to the manual colour balancing process, visual appeal of neighbouring image blocks may vary significantly.There can also be changes in hue inside one image block due to limits of the manual colour balancing procedure.Example of tonal differences of four RGB image mosaics is shown in Figure 2.

Relative radiometric normalization
Various aspects affect to the colour, or tones of the images collected from manned and unmanned aircrafts and form satellites: time of the day and year, atmosphere, illumination conditions, view and illumination angle, BRDF-effects (Bidirectional reflectance distribution function), object, sensor and whole imaging system (Chandelier andMartinoty, 2009, Li et al. 2019).Due to these reasons, creating an evenly coloured mosaic from multiple, overlapping images is a challenge (Gehrke and Beshah, 2016).
In satellite imagery research, method called Relative radiometric normalization, RRN, is a common topic (Syariz et al. 2019).
RRN can be defined as the procedure of adjusting various distortions from grey levels of a multispectral subject image based on a multispectral reference image.(Moghimi et al. 2022).
RRN is often based on pseudo-invariant features (PIFs), that are expected to have stable reflectance or colour on images taken under different conditions.
However, methods based on PIFs are not so well suitable for aerial images due to great difference in adjacent image view angles, and potential specular reflections on images collected from low altitudes (Pan et al. 2010).Major idea presented to airborne images, has been to use radiometric tie-points in overlapping areas between images, and to solve radiometric model parameters through least-squares minimization (Collings et al. 2011, López et al. 2011).This method has been named as radiometric aerial triangulation or radiometric block adjustment, RBB, as it has similarities to more well-known geometric block adjustment.(Chandelier and Martinoty 2009, Honkavaara et al. 2012, Pros et al. 2013, Gehrke and Beshah, 2016).Radiometric block adjustment -based methods have also been applied to images collected with drones (Honkavaara et al. 2012, Miyoshi et al. 2018, Pastucha et al. 2022).

Aim of this work
In this paper, we apply radiometric block adjustment method using RadBA software (Honkavaara et al. 2012(Honkavaara et al. , 2014(Honkavaara et al. , 2018) ) to photogrammetric airborne images collected during national mapping campaigns in Finland.We describe modifications to original drone-based RadBA-process created for hyperspectral sensor, differences due to image types, and challenges faced when the goal is to apply the process to operational orthomosaic production workflow at national mapping agency.
We have three main motivations / goals for this work: 1. To speed up CIR orthophoto deliveries to Finish Food Authority used for EU farming subsidies monitoring 2. Improve the quality of images delivered to Finnish Forest Centre used for forest inventory 3. To automate, fasten and improve the current colour balancing process of orthophotos at NLS

Materials
We used photogrammetric image block collected during three days in August-September 2023 from Lappeenranta, Finland, as test set (60.957N, 28.151E, Figures 1 and 3, Table 1).Additionally, radiometric block adjustment process requires grid of XYZ-points used as radiometric tie points, and digital surface model DSM for final orthomosaic calculation.These were taken from the existing national 2 m grid digital elevation model maintained by NLS.Radiometric tie point file was created with 1 km intervals meaning that there were approximately 6 x 10 points per area that one image covers on the ground.GSD of the DSM for RadBA-based orthomosaic calculation was 15 meters, as it this RadBA-product is only used for visual inspection.

Radiometric block adjustment with RadBA
RadBA -software was originally developed at FGI (Finnish Geospatial Research Institute, part of NLS since 2015) for the quantitative radiometric processing of hyperspectral frame images and creation of image mosaics for images collected from drones (Honkavaara et al. 2012(Honkavaara et al. , 2014(Honkavaara et al. , 2018)).Specifically, various RadBA-processing steps were designed to handle data collected with tuneable Fabry-Pérot Interferometer (FPI) based hyperspectral camera (Saari et al. 2011).In addition to Finnish landscapes, RadBA has also been used in a Brazilian environment constituted of tropical forest and sugarcane plantation areas (Miyoshi et al. 2018).
RadBA adjusts images relatively compared to reference image, and adjustment is done individually and independently for all bands of the sensor.Radiometric tie points from image overlaps and tie point statistics are used for the creation of radiometric model parameters, that are solved through non-linear leastsquares optimization.
, where bBRDFm, m = 1, …, 4, are the BRDF parameters, θi is the incident illumination zenith angle; θr is the reflection view zenith angle; φ = φrφi is the relative azimuth angle, where φr and φi are the azimuth angles of reflected and incident light, respectively.In this study, we calculated individual solar elevation and azimuth angles for each image based on their location and acquisition time.These angles were given as input for RadBA to be used for solving 4-parameter BRDF parameters for the image block.In this study, initial values for the BRDF parameters for R, G and B bands were 0, 0, 0, 0.1, and 0, 0, 0, 0.2 for the NIR band.Standard deviations for the parameters were 0.25, 0.25, 0.25, 0.05 for RGB and 0.25, 0.25, 0.25, 0.1 for NIR band respectively.
RadBA is one exe-file that takes images and input parameters as text files as inputs and produces orthomosaics and several text files as output.We have coded set of python scripts for creating, modifying, and analysing the RadBA input and output text files.These python scripts were modified and simplified be suitable for photogrammetric images.Simplification was possible, as some of the RadBA-processing steps were specific to FPIhyperspectral camera not needed for standard multispectral photogrammetric sensor such as Vexcel UltraCam Eagle.RadBA was run individually for each R, G, B, N bands and PAN-band.

Final steps
As a result, RadBA gives image-and band-wise relative parameter rel_aij values and block-wise BRDF-parameters.It also produces 32 bits per band orthomosaic that can used for further analysis and or visual quality inspection of the process.To achieve goals 1 and 2 mentioned in section 1.3, there are couple of steps needed to produce final images.For the goal 2, original level2 16-bit image pixel values must be adjusted with relative parameter rel_aij and with BRDF-parameters, if BRDFoption was used.
For the goal 1, the image block tones must be adjusted to be visually "good looking", with weight on the CIR orthomosaic.Current procedure is to look at 16-bit histograms of the RadBAproduced CIR orthomosaic (Figure 4), to see how much the histograms can be stretched to better cover the whole 16-bit value range.Same histogram stretch per band is done for each image.Additionally, gamma correction per band is done to improve contrast of the images.Then, 16 bits per pixel images are converted to 8 bits per pixel 4-band images, and GDALpansharpen -algorithm was used to create 4-band hi-resolution 23010 x 14790-pixel size images.Finally, based on these images, orthophoto mosaics were created with standard NLS process with Inpho-software.

Radiometric block adjustment
The RadBA processing with four iterations took time about 10 hours, all bands were calculated in parallel.The processing computer had Windows Server 2016 Standard operating system, 384 GB RAM, and Intel® Xeon® Gold 6148 CPU @ 2.40GHz with 40 cores.
Figure 4 shows the RadBA-processed CIR orthomosaic with 4parameter BRDF correction.This figure is created directly from the 32-bit ENVI-format float mosaics without final steps described in section 2.3, with just min max -clip per band for visualization.Compared to original mosaic (Figure 3), tonal differences due to three different flights have disappeared, and adjacent flight lines that were clearly visible for 15.8.images in the original mosaic are now seamless.
RadBA calculates band-wise average radiometric tie point pixel value standard deviation before and after radiometric block adjustment.From the before and after values, so called homogenization factor is calculated that describes the overall improvement of the mosaic quality.Radiometric tie point pixel standard deviations before and after, and homogenization factors in percent for the test block are shown in Figure 5. RadBAprocess clearly improved the homogeneity of image mosaic for each band; homogenization factors ranged from 0.4 of the NIR band to 0.7 of the blue band.

Final steps
Applying RadBA results, image-and band-wise rel_aij parameters to original level2 16-bit images is a straightforward process of performing just one multiplication per band.Applying 4-parameter BRDF correction (Equation 1) requires more computational effort and is not yet implemented in the operational workflow.After each image and its bands are tuned with rel_aij, contrast stretch and gamma, they were converted from 16 to 8 bits per band as follows: potential values lower than 0 or greater than 2 16 were clipped, and range [0, 65535] was truncated to [0,255].Finally, images were pansharpened with GDAL-pansharpen -algorithm.
After RadBA-results are applied to images, it is important that all the rest of the image operations must be the same for each band, so that the final adjustments do not break the RadBA results.We evaluated different methods to produce visually "good looking" images.First option was to adjust tones of the reference image manually to be visually as good looking as possible, then measure pixel values of bright objects from both the original and adjusted reference image and calculate ratio of adjusted/original to produce band-wise conversion factors.This process is basically one-parameter empirical line method with only scale factor used.The drawback of this method was that it only shifts the locations of the image histograms but does not widen the histogram.
Because of this, after 16 to 8 bits conversion images appeared to be low contrast.After this, band wise histogram stretch and gamma workflow was implemented, as gamma adjusts image tones non-linearly.

Future work
Even though we achieved our goal in performing radiometric block adjustment for airborne photogrammetric images with RadBA-software, there are still several steps to be improved before the workflow can be fully implemented into the NLS orthophoto production.
Speed-up for delivering CIR orthophotos to FFA comes from two reasons: faster orthophoto process with RadBA and possibility of skipping the HISTOQC quality acceptance procedure of our current production workflow.It would be better to have only one, final orthophoto product that has passed the quality acceptance procedure.We still need to test how RadBA-processed images perform in the HISTOQC quality control, and if HISTOQCbased quality control could be improved or replaced with something else.
Current RadBA-process with supporting python codes includes still many manual steps and requires some expert knowledge of the RadBA-software.In the future, we will evaluate all processing steps and seek ways to automate them as much as possible, so that the process could be implemented in the operational NLS orthophoto production workflow.
Workflow presented in this paper performs radiometric block adjustment for single NLS production area (block).To reduce the tonal differences between production areas (Figure 2) it would be desirable to tie radiometric block adjustment of new block to the already adjusted neighbouring block or blocks.This would be possible by adding flight line with images and their rel_aij parameters from the adjusted block to the one under process and set the std_rel_aij for these images to very low value.Now, when performing RBB, these images would basically stay unadjusted (as they already have rel_aij values form earlier process) and tie overlapping images of the new block to them.We will test this procedure of tying neighbouring blocks tonally together in the future works.
In the final steps after RadBA-processing, there are various steps that can be further improved.Main issue is to make the tonal adjustment to visually "good looking" mosaics and to make matching to NLS reference image more automate and robust.It would also be desirable to have just one or a few reference images for the entire Finland, rather than using new ones in each image block.Finally, the GDAL pansharpening algorithm is just a basic one, and there are many other algorithms that could produce higher quality pansharpening.

Conclusions
We have shown that RadBA-based radiometric block adjustment method, originally developed for hyperspectral drone images, can be used for image blocks collected with large-format photogrammetric camera.RadBA-based tonal processing with 4parameter BRDF-correction was able to clearly improve the tonal quality of the image mosaic collected during three different dates.
Results of this work are used to speed up CIR orthophoto deliveries to Finish Food Authority used for EU farming subsidies monitoring, and in improving the quality of images delivered to Finnish Forest Centre used for forest inventory.Our goal is to implement RadBA-based tonal processing to operational NLS orthophoto production workflow.To achieve that, we still need to automate, fasten, and improve various steps in the RadBA-workflow, including optimizing the creation of 8 bits per band orthomosaics from 16 bits per band data.Our longterm goal is to create tonally high quality, seamless orthomosaic of whole Finland.

Figure 1 .
Figure 1.KALLIO image blocks covering whole Finland.Planned three-year cycle 2023-2025 in grey, and 12year cycle 2020-2031 for Lapland in brown.Lappeenranta block used in this work is marked in red.

Figure 2 .
Figure 2. Example of tonal differences between four RGB image mosaics.Image taken from the website https://kartta.paikkatietoikkuna.fi/?lang=en .All mosaics have passed the HISTOQC test.
Central parameter and result of the RRB-process is image-and band-wise relative parameter rel_aij, where i = 1, …, N, N is the number of images and j = 1, …,4 the band, that describes how much pixel values are adjusted relative to reference image/band.As input, standard deviation, std_ rel_aij, for rel_aij is given; typical values for std_rel_aij range from 0.02 for image blocks collected under clear, even illumination conditions to 0.1 for image block collected under varying illumination conditions and/or during multiple days.Another input parameter for RadBA is the standard deviation of the image pixel values in the radiometric tie points, std_dn; typical values for std_dn range from 0.05 for uniform agricultural areas to 0.1 or 0.2 for images containing forested and nonhomogeneous land cover.In this study, values 0.1 for std_rel_aij and 0.05 for std_dn were used, as images were collected during the three days, and typical Finnish forested landscape can be considered homogenous from the 7 km flying height.Size of one radiometric tie point was set to 100 x 100 pixels.RadBA solves values for rel_aij iteratively until desired range is achieved, or maximum numbers of iterations is reached.Honkavaara and Khoramshahi (2018) presented implementation of simple empirical four parameter BRDF-model to RadBA.The 3-parameter model is feasible for data collected during a short period of time, but an extended 4-parameter version is used for cases with a varying sun zenith angle (Equation 4 in Honkavaara and Khoramshahi 2018):

Figure 5 .
Figure 5. Band-wise Radiometric tie point standard deviations before and after RadBA-processing, and homogenization factors in percent for the test block.Central results of the RadBA-process are the band-and imagewise rel_aij parameter values that describe how much pixel values of each image and band is adjusted in relation to the reference image.Rel_aij values for the test block are shown in Figure 6.Rel_aij values ranged from 0.65 to 1.35.From the Figure 6 in can be seen, that the rel_aij values are typically below 1 for the 10.8. and 7.9.flights and close to or greater than 1 for the 15.8.flight.Especially the difference between 15.8.and 7.9.flight rel_aij values is clearly visible as a drop from 1.1 to 0.9.As the test block was collected during three days i.e. with different solar angles, 4-parameter BRDF correction model was used to improve the tonal quality of the image mosaics.BRDF model parameter values and their standard deviations for all bands are shown in Figure7.BRDF-parameter values and their standard deviations were similar for RGB bands and were in the range of 0 and 0.1.Lowest value, -0.022 was for parameter 2 on NIR band, and largest value 0.19 for parameter 4 on NIR band.Values for panchromatic band were similar to RGB bands, and values for NIR band were clearly different from other bands.

Figure 6 .
Figure 6.Relative parameter rel_aij values for all images and bands.Image number 502 (1 st image on the x-axis) was chosen as reference image, so all rel_aij values for 502 are 1.Vertical lines separate images between different dates.

Figure 7
Figure 7. 4-parameter BRDF model parameter values and their standard deviations for all bands.

Table 1 .
Dates, times, average solar angles, and number of images for airborne campaigns for the used image block.Figure 3. Image block used in this study.Original CIR mosaic without any tonal adjustments, mosaicked with most nadir -method.Dots show the image acquisition locations.Date 10.8. in red, 15.8. in green and 7.9. in blue.Background map provided by OpenStreetMap.