Flat-field calibration method for hyperspectral frame cameras

This paper presents a method for characterising spatial responsivity of hyperspectral cameras. Knowing the responsivity of the camera as a function of pixel coordinates allows applying a flat-field correction on image data. The method is based on scanning the field of view of the camera with a broadband radiance source, based on an integrating sphere, and combining the captured frames to synthesise a uniform radiance source filling the whole field of view of the camera at the focus distance. The method was compared with a traditional approach where the aperture of an integrating sphere is imaged from a close distance, filling the entire field of view of the camera. The measurement setup was tested with a hyperspectral camera, based on a tunable Fabry–Pérot interferometer. Without the flat-field correction, the average standard deviation of the pixel responsivities across all the spectral channels of the camera was 3.78%. After the correction, the average standard deviation was reduced to 0.40% and 3.25% for the aperture-scanning method and the close-distance method, respectively. The expanded uncertainty (k  =  2) for the flat-field correction obtained using the scanning method was 0.68%–0.78%, depending on the spectral channel of the camera.


Introduction
In the case of an ideal camera, uniform radiance input would produce a spatially uniform image. However, due to imperfections and physical limitations in imaging hardware, camera systems have non-uniform optical throughput and spatially variable sensor array responsivity [1][2][3][4][5][6][7]. Even with fixed camera parameters, these effects will vary with individual units of the same camera model due to imperfections in the comp onents and the final assembly.
In recent years, miniaturised multispectral and hyperspectral cameras have entered the commercial markets to serve rapidly evolving drone and machine vision applications [8][9][10][11][12][13][14][15][16]. Various filter technologies and multi-sensor solutions employed in these types of cameras inflict additional spatial non-uniformity to the responsivity of the camera, thus further increasing measurement uncertainty.
All components causing non-uniform spatial responsivity of the imaging system can be collectively taken into account by using a flat-field correction [17][18][19][20][21][22]. The flat-field correction matrix for a camera can be obtained by imaging a spatially uniform radiance L e (unit W sr −1 m −2 ) source. Multi-and hyperspectral cameras have been previously calibrated using various approaches, including integrating spheres irradiated with a tungsten-halogen lamp [9-11, 23, 24], using a white plane surface in the laboratory [12,25,26] or in the field [11,14].
This paper presents a flat-field characterisation method based on creating a synthetic uniform radiance source by Metrologia Flat-field calibration method for hyperspectral frame cameras scanning the aperture of an integrating sphere with the camera under test, mounted on a two-axis rotary stage. The method shares the same principle of producing uniform irradiance as described in [27][28][29]. The flat-field characterisation method is tested with a hyperspectral camera, based on a Fabry-Pérot etalon filter, and compared with a commonly used method in which the camera is placed at, or close to, the aperture of an integrating sphere.

Flat-field correction
The flat-field correction can be applied to account for the nonuniform spatial responsivity of the imaging hardware. Figure 1 shows an example photograph illustrating conspicuous radial intensity fall-off towards the edges. This type of imaging artefact is typically caused by a combination of different sources of vignetting. Vignetting is caused by angular dependency of the radiation intensity passing through the camera objective, angular sensitivity of the digital detector array, and parts of the radiation from the peripheries of the field of view being blocked by optical or mechanical elements [7]. Moreover, the spatial responsivity of the camera is affected by the non-uniform responsivity of individual sensor pixels and the transmittance of any additional optical elements, such as filters and beam splitting mirrors.
All these imperfections in captured images can be collectively mitigated with the flat-field correction using the equation where G is the raw image data, D G is the dark signal frame, F is the flat-field correction matrix, and G is the corrected image matrix. Indices i and j are the row and column coordinates of the image matrices, respectively, and index n is the spectral channel of the image. Optimally, flat-field correction matrix F and dark signal frame D G should be acquired using the same exposure time as that employed for capturing image G. In the case of a camera with zoom optics or non-fixed focus, matrices G and F should also be captured with the same zoom and focus settings. Figure 2 shows the measurement setup for obtaining flatfield correction matrix F. The setup consists of an integrating sphere, acting as the broadband radiance source, a two-axis rotary stage for rotating the camera, and a 1.5 m distance rail for altering the measurement distance. The camera was enclosed in a black cabinet with two apertures to reduce stray light in the field of view of the camera.

Measurement setup
The integrating sphere is 30 cm in diameter with an 8.5 cm opening at which the radiance surface is formed. The sphere surface has 97% reflectance barium sulphate coating. The radiant flux was introduced into the sphere using four halogen lamps installed on the same hemisphere as the output aperture, so that none of the flux exited the sphere without first reflecting off the sphere surface. The lamps were connected in series, and supplied with a constant current source. The correlated colour temperature of the radiance source was 2132 K, extending over the spectral range of the camera under test.
Unless either the field of view or the distance of the camera to the sphere is small, the aperture area with emission is not large enough to fill the sensor array of the camera at once. Therefore, the camera was rotated so that all the areas of the array got exposed to the radiance from the source, effectively scanning the whole field of view of the camera.

Camera under test
A commercial hyperspectral camera based on a tunable Fabry-Pérot interferometer (FPI) operating in the visible and near-infrared wavelength range of 502 nm-907 nm was used in this study. The camera is based on the technology described in [30][31][32]. The device uses two single-colour complementary metal oxide semiconductor (CMOS) image sensors and a beam splitting component [33]. One sensor is optimised for  the wavelength range from 500 nm to 636 nm, and the other one from 650 nm onwards. Spectral bands can be selected with 1 nm resolution, and the number of active bands can be up to 380. Camera configuration resulting in 46 spectral channels was used for the measurements.
A Fabry-Pérot interferometer acts as a tunable filter, where the air gap of the interferometer determines the centre wavelength of the transmission band. The final spectral response of the camera is dependent on the signal passing through the FPI and the spectral characteristics of the detector. During data collection, the size of the air gap is varied to acquire a sequence of spectral images in the desired spectral range.
The camera was equipped with fixed optics with a focal length of 9.0 mm and an f number of 2.8, producing 37 • field of view. The pixel resolution of both detector arrays was 1010 × 1010 pixels with a pixel size of 5.5 µm. The minimum distance for the camera at which the imaged objects are in focus is 1 m.

Measurement sequence
The measurement scan was performed at a distance of 1.0 m from the opening of the integrating sphere, with the angular step size of 1.2 • around both axes producing large overlaps between the areas with direct signal from the radiance source. Figure 3 shows the number of overlapping frames averaged for every pixel of the camera under test. The x i and y i (subscript i for image) are the columns and rows in the pixel coordinate system, respectively. The range of overlapping frames per pixel was 6-13, with an average of 11 frames per pixel. The figure also illustrates the size of single sphere-aperture images on the camera sensor, when imaged from the 1 m distance. At that distance, one image of the sphere aperture covered approximately 5% of the camera sensor area.
In addition to the scanning measurement, stationary images were captured close to the aperture of the sphere to enable comparison of the two methods used for the flat-field correction. The close-up measurements were performed from the distance of 4 cm with the image of the sphere opening covering the entire field of view of the camera. The total distance to the back of the sphere was approximately 34 cm, which is less than the minimum focus distance of the camera.
Thirty dark frames were captured before and after each measurement sequence with the same exposure time as the actual measurements. Before capturing any images, the radiation source and the camera were allowed to stabilise for a minimum of two hours.

Calculation of the flat-field matrix
To obtain the flat-field correction matrix F, the data obtained using the scanning procedure were merged. The merging algorithm consists of five steps: removing the dark signal, thresholding out the pixel areas not irradiated directly by the source, removing the edges of the irradiated area in the images, combining the individual frames from the scan into one image, and finally filtering the image to remove high-frequency components.
First, the dark signal frames are averaged and subtracted from each frame of the scan. Next, the pixel values not irradiated by the source are thresholded out by replacing all lowintensity pixel values with a non-numerical value. By using non-numerical values for the areas without direct signal from the radiance source, these areas are prevented from contributing to the average pixel intensities of the combined frames. This is required because the number of frames with signal is not constant across the pixels.
After thresholding, the very edges of the aperture, whose radiance may considerably deviate from that of the rest of the aperture area, are removed by convolving each frame with the following kernel In computing, such an operation results in non-numerical values for all those pixels which are adjacent to a non-numerical-value pixel, while leaving the rest of the pixel values unchanged. The number of padding zeroes can be increased in accordance with the pixel resolution of the image and the thickness of the non-uniform edge area. In this study, the kernel size for the edge removal was 9 × 9. After removing the edges, all the frames are averaged into a single image matrix for each spectral channel n, ignoring all non-numerical values. Figure 4 shows the principle of combining the frames. The smaller the step when turning the camera during the scan, the more the irradiated areas overlap, thus resulting in increased averaging of the signal, which reduces the impact of the spatial non-uniformity of the radiance source [28]. Finally, flat-field correction matrix F is obtained by convolving each spectral channel with a Gaussian filter kernel. This spatial low-pass filtering reduces the impact of random noise in the measurement data and the spatial non-uniformity of the radiance source output. The filtering is a compromise between reducing the impact of these imperfections and obtaining the explicit responsivities of the individual pixels.

Performance assessment
To assess the performance of the obtained flat-field calibration, and to compare the methods, the flat-field matrices obtained using the scanning and the close-up methods were used to correct a reference image according to equation (1). The reference image was created by scanning the aperture of the integrating sphere from 1.4 m distance. For comparison, the reference image was also corrected using the factory calibration of the camera. To quantify the results, the channelwise standard deviations, relative to the mean pixel intensity of the channel, were calculated for the reference image pixel intensity values before and after the flat-field corrections for the three different methods. Another reference image for the correction was created by imaging a 0.5 cm polytetrafluoroethylene (PTFE) diffuser plate backlit by a matrix of halogen lamps. The diffuser plate was placed close to the camera. The channel-wise standard deviations were also calculated for this second reference image before and after the corrections. Figure 5 shows a comparison of the normalised spatial responsivity, or the flat-field correction matrix, of an example channel, obtained using the scanning method at 1 m distance and by imaging the aperture of the integrating sphere from a close distance of 4 cm. Essentially, the close-up method produces an out-of-focus image of the back wall of the sphere. The spatial uniformity of the radiance of this area is influenced by irradiation of the sphere wall and the spatial nonuniformity of the surface. In this measurement setup, uneven irradiation is particularly severe because the four lamps irradiate the imaged area directly, which means that the angular intensity distributions of the lamps influence the image [34]. These factors have a smaller impact when the output aperture of the sphere is imaged from a larger distance and the overlapping images are averaged. Figure 6 shows the variations of the pixel intensities per spectral channel when the combined frames captured from 1.4 m distance were flat-field corrected using the correction matrices obtained from 1.0 m scan and from the close-up measurements. The results have been calculated as a standard deviation across the whole image area for each spectral channel. The channel wavelengths are the nominal central wavelengths according to the camera setting. The mean standard deviation was 3.78% for the raw uncorrected image data, 0.40% for the data corrected using the scanning method, and 3.25% for the data corrected using the close-up method. For the corrected pixel intensities obtained using the factory calibration of the camera, the mean standard deviation was 1.72%.

Flat-field correction matrix
For the reference image which was flat-field corrected using the scanning method, the standard deviation of pixel intensities decreases as the spectral channel wavelength increases. This is likely caused by the higher signal-to-noise ratio at the longer wavelengths, due to the halogen-lamp-based radiance source. The factory calibration of the camera improved the uniformity of almost every channel, but the channels with high initial non-uniformity retained residuals of the gradients in their spatial responsivity.
For the second reference image, captured through the diffuser plate, the mean standard deviations were 6.16%, 4.02%, 5.97%, and 4.29% for the raw uncorrected image data, for the data corrected using the scanning method, for the data corrected using the close-up method, and for the factory calibration, respectively. Although the flat-field matrix obtained from the 1 m scan reduced the variations of the pixel intensities the most, the general increase in the standard deviations, when compared with the first reference image, suggests that spatially and angularly non-uniform properties of the diffuser, and the possible non-uniformities of the irradiation on the diffuser plate, add to the intensity variation of the image data.
The flat-field matrix obtained using the close-up method reduced the responsivity variation for the channels with large non-uniformities, i.e. the channels for 502 nm and 628 nm-646 nm. For relatively uniform channels, the close-up method can have an adverse effect, as can be seen in figure 6. The nearly constant standard deviation shows that the non-uniformity of the spatial responsivity of the camera is outweighed by the non-uniformity of the imaged surface. Figure 7 shows the uneven irradiation of the back wall of the sphere. For the figure, an image captured with the close-up method has been flat-field corrected, using the matrix from the scanning method according to equation (1), and rotated 90 • counterclockwise, as the camera was sideways in the close-up measurements. The increased pixel intensity levels in three of the corners are caused by the reflections of the direct irradiance coming from the sphere lamps. The small, slightly darker area in figure 7 at ( x i , y i ) = (450, 550) is caused by an internal reflection within the camera during the scanning measurement of the flat-field matrix. The reflection is view-angle dependent, and was thresholded out elsewhere except for this one area where it overlapped with the direct radiance signal, creating an area with 1%-2% error in the obtained flat-field matrix. In the case that the camera under test produces this type of artefact, it can be corrected, for example, by shifting the camera laterally and rescanning that area. Figure 8 shows the responsivity of seven example spectral channels across a diagonal of the obtained flat-field matrix F. The three spectral channels for 628-646 nm that coincide with the change of the sensor exhibit strong gradiental behaviour across the vertical axis of the image area. The responsivity of the 628 nm channel (longest wavelength of the first sensor) decreases towards the bottom of the detector array, while the 640 nm and 646 nm channels (shortest wavelengths of the second sensor) have increasing responsivity towards the bottom of the sensor. This may be caused by the fact that the cut-off wavelength of the dichroic mirror, used to split the radiation, is dependent on the incident angle of the radiation.
The 502 nm channel has a vertical responsivity gradient decreasing towards the bottom of the sensor. The 502 nm channel of the camera has a filter leak at approximately 640 nm-650 nm. Due to the considerably lower radiance level of the source at 502 nm compared with 640 nm-650 nm, the combination of the leak and the dichroic mirror is likely the reason for this gradiental responsivity.
Besides the filter and mirror specific non-uniformities, in all the channels, including the ones with strong responsivity gradient, there are features which can be attributed to vignetting. For instance, the intensity fall-off towards the edges, and the abrupt intensity drop in the corner (x i , y i ) = (1, 1).

Uncertainty analysis
The expanded uncertainty (k = 2) for the flat-field correction matrix F obtained using the scanning method is 0.68%-0.78%, depending on the spectral channel of the camera. The uncertainty budget is presented in table 1. The largest source of uncertainty for the flat-field correction is the spatial nonuniformity of the radiance source. The uncertainty components due to the temporal stability of the source and the sensor noise of the camera are of type A. The components due to spatial uniformity and temperature sensitivity are of type B in individual frames. However, due to the averaging nature of the data processing algorithm, the former uncertainties reduce statistically, and are thus of type A in the obtained flat-field matrix F.
The uncertainties were estimated using Monte Carlo simulation by deviating pixel intensity values of individual frames measured during the 1 m scan, before processing them according to the steps described in section 2.5. The effects of the deviations were evaluated by employing the original flat-field matrix F as the image to be corrected, or the numerator in equation (1), and dividing it pixel-wise by the matrix obtained when the randomly deviated frames were combined. The standard deviations of each quotient matrix G n were used to quantify the combined measurement uncertainty. The impact of individual components was estimated by excluding the other ones from the simulation.
The standard deviation of the spatial uniformity of the radiance source was evaluated to be 1.2%. The uniformity was determined from the pixel intensity values of a flat-field corrected hyperspectral image of the integrating sphere aperture. Figure 9 shows one channel of the flat-field corrected image of the sphere aperture. The image was obtained by averaging thirty stationary images captured at a 1 m distance from the source, and removing the edge area of the aperture. The non-uniformity is mostly due to the increased intensity in the lower left area of the aperture, which is caused by the first order reflection of a light source. The range of the intensity values across the aperture was 11%. The uncertainty component for the budget was determined using the Monte Carlo simulation by adding a linear-intensity gradient with 11% range to the aperture area of the individual scan frames. The direction of the gradient was randomly selected for every execution of the Monte Carlo simulation.
The uncertainty component due to the spatial non-uniformity of the radiance source depends on the density of the captured frames during the scan. The smaller the angular step of the rotary stage, the higher is the number of captured images and the number of achieved overlaps per one single pixel. For instance, when averaging 20 × 20 evenly distributed image circles shown in figure 9, which cover the whole 1010 × 1010 pixel image area, the standard deviation of the pixel values of the resulting image is 0.50%. If the number of individual frames is increased to 40 × 40, the standard deviation decreases to 0.24%.
The temporal stability of the source was studied by monitoring the output of the source for four hours with a temperature-stabilised luminance meter after the four lamps had already been on for two hours. The standard deviation of the recorded signal was 0.003%, which also includes the impact of the luminance meter stability. This value was used as the standard deviation for the random offset level of the radiance source intensity.
The impact of sensor noise of the camera on the measurement uncertainty was studied by capturing thirty images of  0.68-0.78 Figure 9. Flat-field corrected image of the integrating sphere aperture captured using 772 nm channel. The edge area of the aperture has been removed from the image.
the sphere aperture with the camera stabilised and stationary. The frequency of capturing the images was the same as during the scans in order to warm up the camera similarly. The magnitude of the sensor noise was determined from the changes of the raw pixels values for the constant radiance signal. The standard deviation of the sensor noise ranged from 1.2% to 3.2% depending on the camera channel. Gaussian white noise with the corresponding standard deviation of each channel was added to the respective scan-frame channels in the Monte Carlo simulation. The captured images together with the temperature values recorded by the camera were used to determine the temperature dependency of the responsivity of the camera. The temper ature sensitivity of the camera was evaluated by turning on the camera and capturing images of the radiance source in rapid succession to heat up the camera from the room temperature to 38 • C. During the flat-field scans, the temperature of the camera stayed within 33.5 °C-34.3 °C. In that range, the temperature sensitivity of the camera varied from 0.0% per 1 • C to 0.3% per 1 • C, depending on the spectral channel. A drift of 0.8 • C in temperature during the flat-field scan causes less than a 0.25% gradient of systematic error (peak-to-peak) in the obtained flat-field correction matrix. Half of the latter value was used as the standard uncertainty related to temperature sensitivity in table 1.

Conclusion
This paper presents a method for characterising the spatial responsivity of hyperspectral cameras to allow applying flatfield corrections to measured image data. The method is based on scanning the field of view of the camera with the opening of an integrating sphere in order to synthesise a uniform radiance source. The combined scan data are used to obtain the flat-field correction matrix. The measurement setup built for the study consisted of a radiance source based on a 30 cm integrating sphere and a two-axis rotary stage for turning the camera.
The method was tested by characterising a tunable Fabry-Pérot interferometer hyperspectral camera, and flat-field correcting a uniform reference image, which was obtained by scanning the radiance source from a different distance. The obtained correction was compared with the results produced by an alternative method where the flat-field correction matrix is measured by imaging the aperture of the integrating sphere at a close distance, exposing the entire camera sensor to the radiance source at once.
The uniformity of the flat-field corrected reference image was quantified by calculating the average standard deviations of the pixel intensities from the mean intensity across all the spectral channels of the camera. For the raw, uncorrected data the mean standard deviation was 3.78%. Using the correction matrix obtained with the scanning method, the mean standard deviation was reduced to 0.40%. For comparison, the flat-field correction matrix measured from a close distance from the aperture of the sphere reduced the mean standard deviation to 3.25%. The expanded uncertainty (k = 2) for the flat-field correction obtained using the scanning method was 0.68%-0.78%, depending on the spectral channel of the camera.
In addition to hyperspectral frame cameras, which are increasingly used in various applications, such as environmental measurements or industrial process control, the method is also suitable for other types of frame cameras. Compared with the approach where the aperture of an integrating sphere is imaged at a close distance, the scanning method is particularly advantageous for characterising the spatial responsivity of wide-angle-lens cameras. For instance a fisheye-lens camera placed at the aperture of an integrating sphere would inevitably see the structural elements of the integrating sphere within its field of view, which would affect the obtained flatfield matrix.

Acknowledgment
The work leading to this study is partly funded by the Academy of Finland project Quantitative remote sensing by 3D hyperspectral UAVs-From theory to practice. The authors thank E Vaskonen for the work on the measurement setup, and O Pekkala, N Viljanen, and T Hakala for the insights on the par ticulars of the device under test. A Kokka acknowledges the funded position of Aalto ELEC doctoral school and the support of KAUTE Foundation.