Morphological Band Registration of Multispectral Cameras for Water Quality Analysis with Unmanned Aerial Vehicle

Multispectral imagery contains abundant spectral information on terrestrial and oceanic targets, and retrieval of the geophysical variables of the targets is possible when the radiometric integrity of the data is secured. Multispectral cameras typically require the registration of individual band images because their lens locations for individual bands are often displaced from each other, thereby generating images of different viewing angles. Although this type of displacement can be corrected through a geometric transformation of the image coordinates, a mismatch or misregistration between the bands still remains, owing to the image acquisition timing that differs by bands. Even a short time difference is critical for the image quality of fast-moving targets, such as water surfaces, and this type of deformation cannot be compensated for with a geometric transformation between the bands. This study proposes a novel morphological band registration technique, based on the quantile matching method, for which the correspondence between the pixels of different bands is not sought by their geometric relationship, but by the radiometric distribution constructed in the vicinity of the pixel. In this study, a Micasense Rededge-M camera was operated on an unmanned aerial vehicle and multispectral images of coastal areas were acquired at various altitudes to examine the performance of the proposed method for different spatial scales. To assess the impact of the correction on a geophysical variable, the performance of the proposed method was evaluated for the chlorophyll-a concentration estimation. The results showed that the proposed method successfully removed the noisy spatial pattern caused by misregistration while maintaining the original spatial resolution for both homogeneous scenes and an episodic scene with a red tide outbreak.


Introduction
In multispectral images, precise registration of multispectral bands is critical for subsequent quantitative data analysis, which relies on the "spectrum" of the signal (e.g., radiance or reflectance) from the target. If the bands are not perfectly aligned with each other, due to the reasons such as lens distortion, displacement in lens location, and inaccurate geometry transformation between the locations, spectral radiometric values in a fixed pixel location may originate from different targets. Algorithms that depend on the band ratios, or the band difference, are particularly sensitive to the quality of the band registration and may produce significant errors for an inhomogeneous target area if misregistration exists.
There are various sources for misregistration: a difference in the lens locations for each band, a difference in the image acquisition timing accompanied by fast-moving targets, etc. Commercial multispectral cameras typically include individual lenses for the multi-bands and have different exposure times to maximize the effective radiometric range for various targets. The differences in the lens locations for the multi-bands can be modeled via a projective transformation if the image is assumed to be free of nonlinear image distortions, such as radial distortion; the differences in the viewing geometry can be corrected to a reference band using band-to-band projective transformations that specifies eight parameters for rotation (1), translation (2), isotropic scaling (1), anisotropic scaling (1), skew (1), and perspective shortening (2) (figures in the parenthesis denotes the number of free variables for the quantity) [1,2]. The methods based on Fourier transform does not require the time-consuming process of finding matching points between two images, and effectively register multiple bands solely based on its spatial frequency pattern [3].
However, such transformation approaches that rely on geometric characteristics of the scene cannot effectively address the cases having non-rigid body target deformation, where the forms of targets may vary between the bands. This issue is prominent when analyzing the color of water where the targets (i.e., water surface) move or deform quickly during a short time interval (<1 s) between the acquisition of different bands. As shown in Figure 1, multispectral band images for ocean surface in the normal coastal area in Korea reveal that the differences in water reflectance, and its spatial pattern, clearly do not correspond to rigid-body transformation; thus, the differences cannot be resolved by a projective transformation. Note that the band images shown in Figure 1 have already undergone band-to-band registration through a projective transformation.
Remote Sens. 2020, x, x FOR PEER REVIEW 2 of 21 quality of the band registration and may produce significant errors for an inhomogeneous target area if misregistration exists. There are various sources for misregistration: a difference in the lens locations for each band, a difference in the image acquisition timing accompanied by fast-moving targets, etc. Commercial multispectral cameras typically include individual lenses for the multi-bands and have different exposure times to maximize the effective radiometric range for various targets. The differences in the lens locations for the multi-bands can be modeled via a projective transformation if the image is assumed to be free of nonlinear image distortions, such as radial distortion; the differences in the viewing geometry can be corrected to a reference band using band-to-band projective transformations that specifies eight parameters for rotation (1), translation (2), isotropic scaling (1), anisotropic scaling (1), skew (1), and perspective shortening (2) (figures in the parenthesis denotes the number of free variables for the quantity) [1,2]. The methods based on Fourier transform does not require the time-consuming process of finding matching points between two images, and effectively register multiple bands solely based on its spatial frequency pattern [3].
However, such transformation approaches that rely on geometric characteristics of the scene cannot effectively address the cases having non-rigid body target deformation, where the forms of targets may vary between the bands. This issue is prominent when analyzing the color of water where the targets (i.e., water surface) move or deform quickly during a short time interval (< 1 s) between the acquisition of different bands. As shown in Figure 1, multispectral band images for ocean surface in the normal coastal area in Korea reveal that the differences in water reflectance, and its spatial pattern, clearly do not correspond to rigid-body transformation; thus, the differences cannot be resolved by a projective transformation. Note that the band images shown in Figure 1 have already undergone band-to-band registration through a projective transformation. This type of image registration, which involves a non-rigid body transformation, has been investigated for morphological image registration [4][5][6][7]. It is applicable when the images of two different targets, or of one target object that experienced a non-rigid body deformation, are expected to have a similar internal structure and most of the image contents have common features. The basic mathematical tools for morphological image transformation include calculus of variations, optimization with regularization or constraints, and the derivation of invariant features. The methods are typically used in medical image processing, such as computed tomography and This type of image registration, which involves a non-rigid body transformation, has been investigated for morphological image registration [4][5][6][7]. It is applicable when the images of two different targets, or of one target object that experienced a non-rigid body deformation, are expected to have a similar internal structure and most of the image contents have common features. The basic mathematical tools for morphological image transformation include calculus of variations, optimization with regularization or constraints, and the derivation of invariant features. The methods are typically used in medical image processing, such as computed tomography and magnetic resonance imaging, for diagnostic purposes; however, applications to remote sensing are rare, particularly for the band registration, because typical observation targets in optical remote sensing are stationary objects, such as land surfaces and man-made structures. However, because water surfaces move fast, owing to ocean currents and wave motion, multispectral images from ships and low-altitude platforms, such as unmanned aerial vehicles (UAVs), clearly exhibit locational errors.
In this study, we develop a novel morphological band registration technique, designed for high-resolution water quality analysis, which preserves the true spectrum of fast-moving targets. The proposed method exploits the quantile plots between the bands to accurately determine the radiometric correspondence of the pixels in different bands. For multispectral images, a Micasense Rededge-M camera was operated onboard a UAV, DJI Inspire-2, and ocean surface images were acquired from coastal areas of Korea, of various altitudes and biological conditions. The specification and radiometric properties of the multispectral camera and the image data are described in Section 2 (Materials). The radiometric/geometric preprocessing, derivation of the "remote sensing reflectance" for water quality analysis are presented in Section 3 (Methodology and Analysis), along with the analysis on the adverse impact of the residual misregistration on the water quality variable estimation. In Section 4 (Algorithm Development and Assessment), the proposed morphological registration scheme is described in detail and the development of the entire correction procedure is presented. The algorithm results are demonstrated for multiple test images, taken at various altitudes. In Section 5 (Discussion and Conclusion), the correction results and remaining tasks are discussed.

Micasense Rededge-M and Data
The Rededge-M camera has five spectral bands, the center wavelengths of which are located at 475, 550, 668, 717, and 840 nm ( Figure 2). The red (668 nm) and near-infrared (NIR) (717 nm) bands are designed to capture the red edge feature, which is salient in vegetation, and to quantify the photosynthetic pigments via indices, such as the normalized differenced vegetation index [8] and soil-adjusted vegetation index [9]. The blue band (475 nm) is useful when quantifying pigment absorption in water when it is referenced with respect to the green band (550 nm) [10]. In water color analysis, the last NIR band (840 nm) is particularly useful for quantifying atmospheric scattering between the target and the sensor because clear water theoretically has zero water-leaving radiance in the 840 nm band [11]. For turbid waters, the radiance at 840 nm is often large and can thus be used to detect the existence of suspended sediments in water [12]; however, as a result, the estimation of atmospheric effects becomes more complicated [13,14]. The radiometric sensitivity of the Rededge-M was tested for water color analysis in Kim et al. (2019) [15] in a brief experiment in which the radiometric data from a Rededge-M camera was compared with that from a hyperspectral radiometer, TriOS RAMSES. The study showed that the Rededge-M was able to retrieve a comparable spectral shape for a water body (which typically has a low radiance level compared to terrestrial targets) when calculated using remote sensing reflectance (R rs ).
In this study, four Rededge-M image sets from multiple field campaigns were used for the development of a morphological registration algorithm. Table 1 shows the dates, locations, and altitudes of the camera images that were used for the analysis, and the study site is presented in Figure 3a. All images were captured by a drone, DJI Inspire-2. The camera body and the downwelling irradiance sensor were installed on the drone using a simple bracket and a damper (Figure 3b). The Rededge-M camera was installed with a fixed viewing zenith angle of 40 • and the camera attitude was controlled to head north to constrain the relative azimuth angle within 90-135 • with respect to the sun direction, to minimize the surface reflectance [16,17]. The Zenmuse-X5s camera was installed in front of the Rededge-M to capture a wider-angle overview of the target scene with higher spatial resolution. RGB images of water reflectance are presented in Figure 4 for all four scenes. In this study, four Rededge-M image sets from multiple field campaigns were used for the development of a morphological registration algorithm. Table 1 shows the dates, locations, and altitudes of the camera images that were used for the analysis, and the study site is presented in Figure 3a. All images were captured by a drone, DJI Inspire-2. The camera body and the downwelling irradiance sensor were installed on the drone using a simple bracket and a damper ( Figure 3b). The Rededge-M camera was installed with a fixed viewing zenith angle of 40° and the camera attitude was controlled to head north to constrain the relative azimuth angle within 90°-135° with respect to the sun direction, to minimize the surface reflectance [16,17]. The Zenmuse-X5s camera was installed in front of the Rededge-M to capture a wider-angle overview of the target scene with higher spatial resolution. RGB images of water reflectance are presented in Figure 4 for all four scenes.    In this study, four Rededge-M image sets from multiple field campaigns were used for the development of a morphological registration algorithm. Table 1 shows the dates, locations, and altitudes of the camera images that were used for the analysis, and the study site is presented in Figure 3a. All images were captured by a drone, DJI Inspire-2. The camera body and the downwelling irradiance sensor were installed on the drone using a simple bracket and a damper ( Figure 3b). The Rededge-M camera was installed with a fixed viewing zenith angle of 40° and the camera attitude was controlled to head north to constrain the relative azimuth angle within 90°-135° with respect to the sun direction, to minimize the surface reflectance [16,17]. The Zenmuse-X5s camera was installed in front of the Rededge-M to capture a wider-angle overview of the target scene with higher spatial resolution. RGB images of water reflectance are presented in Figure 4 for all four scenes.

Acquisition Time Difference in Rededge-M Band Images
The Rededge-M employs a proprietary Auto Gain Control (AGC) algorithm that works to minimize the number of overexposed pixels, but the number of overexposed pixels will never be zero because the AGC also wants there to be a maximal number of properly exposed pixels. The AGC optimizes the gain (ISO) and exposure of each capture for each of the five imagers such that the resulting picture is properly exposed. The Rededge-M has five imagers that trigger the top of each frame together, and one "capture" is created during each triggering, which is represented by five different frames, one for each wavelength band. Each imager has a different filter as to only capture data from the wavelength band of interest. Because of the AGC, the ISO Speed and exposure time (which can be inspected in the metadata of each frame) may vary by band on an individual capture, and different captures taken during the same flight will also vary. The small differences between the exposure times among frames of a capture usually don't make a difference. However, motion blur may occur when there is a large difference (i.e., 1 ms vs. 5 ms) in exposure time between multiple frames in a single capture. If the camera is mounted on an aircraft and is in motion, the frames with longer exposures will have motion blur, and will have a slightly different geometric offset compared to the shorter exposure time. The geometric offset between frames with different exposure times on a single capture will be a function of the angular rate of the camera and the respective exposure times of the frames. This will manifest as motion blur, and will result in a difference in average pointing angle of -, where is the pointing angle from a frame with a longer exposure time and is the pointing angle from a frame with a shorter exposure time. With all this information in mind, while each frame will have had the top of the frame triggered at the

Acquisition Time Difference in Rededge-M Band Images
The Rededge-M employs a proprietary Auto Gain Control (AGC) algorithm that works to minimize the number of overexposed pixels, but the number of overexposed pixels will never be zero because the AGC also wants there to be a maximal number of properly exposed pixels. The AGC optimizes the gain (ISO) and exposure of each capture for each of the five imagers such that the resulting picture is properly exposed. The Rededge-M has five imagers that trigger the top of each frame together, and one "capture" is created during each triggering, which is represented by five different frames, one for each wavelength band. Each imager has a different filter as to only capture data from the wavelength band of interest. Because of the AGC, the ISO Speed and exposure time (which can be inspected in the metadata of each frame) may vary by band on an individual capture, and different captures taken during the same flight will also vary. The small differences between the exposure times among frames of a capture usually don't make a difference. However, motion blur may occur when there is a large difference (i.e., 1 ms vs. 5 ms) in exposure time between multiple frames in a single capture. If the camera is mounted on an aircraft and is in motion, the frames with longer exposures will have motion blur, and will have a slightly different geometric offset compared to the shorter exposure time. The geometric offset between frames with different exposure times on a single capture will be a function of the angular rate of the camera and the respective exposure times of the frames. This will manifest as motion blur, and will result in a difference in average pointing angle of θ 2 -θ 1 , where θ 2 is the pointing angle from a frame with a longer exposure time and θ 1 is the pointing angle from a frame with a shorter exposure time. With all this information in mind, while each frame will have had the top of the frame triggered at the same instant, the exposure time metadata may be different, resulting in slightly different exposure times among a single capture. For example, the exposure time of the 5 band images of Scene-A are 1/741, 1/585, 1/780, 1/367, and 1/356 s, respectively, for 475, 550, 668, 717, and 840 nm, causing the targets to be captured in different status.

Radiometric and Geometric Calibration
Basic radiometric preprocessing was performed using the processing modules provided on the Micasense Github page [18]. A Vignette correction was first performed for each band image and these Vignette-corrected band images were input to the radiometric calibration process to produce the radiance data. The radial and tangential distortion were corrected according to following formula, where u and v are image coordinates, r = √ u 2 + v 2 , k 1 , k 2 , k 3 are coefficients for radial distortion, and p 1 , p 2 are for tangential distortion. As shown in Figure 1, the Rededge-M acquires radiance at five wavelengths, through five individual lenses, inevitably leading to misalignment in the band images. The misalignment can be corrected using a projective transformation, constructed by matching numerous matching points between two band images. The projective transformation between the bands is where x and x are 3 × 1 homogeneous vectors of image coordinates in two band images and M is a 3 × 3 non-singular matrix for the projective transformation.

Water Color Analysis
To conduct water color analysis for the estimation of in-water constituents (e.g., chlorophyll-a concentration), remote sensing reflectance must first be derived from the radiance measurements. Remote sensing reflectance (R rs ) can be calculated as where L wT is the total radiance from water, L sky is the downward radiance from the sky, E d is the downward irradiance, and ρ is the Fresnel reflectance factor [19,20]. Setting aside ρ, to derive R rs in the field, three radiometric measurements are required for each scene: L wT , L sky , and E d . The first two radiance variables-L wT and L sky -are acquired by capturing the water surface and sky using the Rededge-M, following the measurement protocol suggested in the ocean color analysis [16]. For both observations, the recommended azimuth angle is 135 • , with respect to the sun azimuth, and the recommended zenith angles are 45 • and −45 • for the water and sky, respectively. The measurement protocol is intended to minimize the variation in ρ, which varies from 0.02 to 0.07, depending on wind speed and sun-sensor-target geometry [17], where the factor is confined to an approximate range of 0.02 to 0.025 when the aforementioned measurement protocol is observed. However, high-altitude drone images with a wide viewing angle lead to a wide range of viewing zenith and azimuth angles, in which case, ρ varies significantly outside the 0.02 to 0.025 range, requiring a Remote Sens. 2020, 12, 2024 7 of 19 pixel-based adaptive ρ estimation. A simple method to determine ρ adaptively for different locations is by exploiting the fact that the total water radiance at 840 nm is solely attributed to the surface-reflected radiance (not to the water-leaving radiance from the water body). This assumption holds when the water is clear; thus, the water-leaving radiance at 840 nm is nearly zero, where L w and L s f c denote the water-leaving radiance and the surface-reflectance radiance, respectively. If L w (840) = 0, then L wT (840) = L s f c , leading to ρ = L wT (840)/L sky (840). There are two options to determine the downwelling irradiance: (1) using the DLS and (2) the reference panel. The irradiance from the panel reflectance can be calculated as where L re f is the radiance from the reference panel and r panel is the reflectance of the panel. If the two instruments are located and perfectly calibrated, the results of the two calculations should theoretically match. In this experiment, the DLS was attached to the UAV and the reference panel measurements were made on the ship, causing a difference in the altitude. In this experiment, the irradiance difference caused by the atmospheric conditions (water vapor, aerosol, etc.) was approximately 15-20%. Because our focus is on the water surface, we opted to use the irradiance measured at the ship level, via the reference panel.
After R rs is obtained, it can be used to derive bio-geochemical variables, such as chlorophyll-a(Chla). To examine the effect of misregistration between the bands and the performance of the proposed morphological registration method, retrieval results are computed for chlorophyll concentrations, which is one of the most central biological quantities in the water quality analysis. For Chla concentrations in a non-turbid ocean condition, the OC2 algorithm, which utilizes one blue and one green band, was used [10,21].
where R rs (λ) is the remote sensing reflectance for the wavelength λ and a i 's are the algorithm coefficients. The OC2 coefficients for the Landsat-8 operational land imager were used (482 nm for blue and 561 nm for green) for the Rededge-M, whose blue band is centered at 475 nm and the green at 550 nm [21]. For the tested scenes with a red tide outbreak, the red-to-blue ratio (RBR) algorithm [22], developed for the geostationary ocean color imager (490 nm for blue and 680 nm for green), was used to retrieve the chlorophyll contents in the bloom. It is important to note that the algorithm coefficients for OC2 and RBR were not specifically tuned to the Rededge-M in this study because the focus of the study is not on the precise retrieval of Chla concentrations but the analysis of the impact of misregistration (particularly spatial pattern). The band centers in the original OC2 and RBR algorithms do not significantly differ from those of the Rededge-M, from which we can reasonably assume that the spatial anomaly pattern would be similar, even after the fine calibration of the algorithms to Rededge-M. Figure 5a,b shows the RGB images of R rs for Scene-A, for the fixed Fresnel factor (ρ =0.025) and adaptive Fresnel factor cases, respectively. R rs with a fixed ρ demonstrates that slant viewing angles cause a high surface reflectance in the upper right corner of the image. It can be observed that wave facets of different surface normals also led to varying viewing geometry, resulting in a variation of the R rs estimation, which exhibits residual sky reflectance on the surface. On the contrary, the adaptive approach demonstrates that the variation caused by the viewing geometry is significantly reduced with less across-image R rs variation and smaller residual surface reflectances by the wave facets. Figure 5c,d are the OC2 Chla estimates, derived from the data, with a fixed and adaptive , respectively. The large residual surface reflectance caused by the slant viewing angles in the upper right corner led to significant underestimates of Chla and inflation in the bottom left corner. The anomalies are less significant in the adaptive case; however, they have not been completely removed, even with the adaptive scheme. This implies that the surface-reflectance mechanism is more complex than what is described by the adaptive scheme model (e.g., the existence of a nonlinear band-by-band behavior). The phenomenon to focus on here is the large and noise-like Chla variation in a small-scale window. For a more detailed analysis, the subset areas marked by the red rectangles in Figure 5c,d were displayed in Figure 6. The Chla subset images show that the pixel-to-pixel variation is significantly large for both the fixed and adaptive approaches and such a high-frequency pattern is not caused by the real Chla spatial variation in the field. The adaptive approach exhibited a similar degree of variation to the fixed approach, which reveals that the noise-like pattern is not from the variation of wave facets. The images of Band 4 and Band 5 support this interpretation because the spatial pattern of the reflectance in the two bands are consistent with each other and it reflects the wave facet distribution (note that the two NIR bands have values of almost zero; therefore, the residual surface reflectance mostly contributes to the reflectance of the bands). The spatial variation of Chla and NIR appears to have no clear spatial correlation, implying that the high Chla variation in the images is not caused by the residual surface reflectance (equivalently, the viewing geometry) but from a factor related to the image quality or pixel registration. Figure 5c,d are the OC2 Chla estimates, derived from the R rs data, with a fixed ρ and adaptive ρ, respectively. The large residual surface reflectance caused by the slant viewing angles in the upper right corner led to significant underestimates of Chla and inflation in the bottom left corner. The anomalies are less significant in the adaptive ρ case; however, they have not been completely removed, even with the adaptive scheme. This implies that the surface-reflectance mechanism is more complex than what is described by the adaptive scheme model (e.g., the existence of a nonlinear band-by-band behavior). The phenomenon to focus on here is the large and noise-like Chla variation in a small-scale window. For a more detailed analysis, the subset areas marked by the red rectangles in Figure 5c,d were displayed in Figure 6. The Chla subset images show that the pixel-to-pixel variation is significantly large for both the fixed and adaptive approaches and such a high-frequency pattern is not caused by the real Chla spatial variation in the field. The adaptive approach exhibited a similar degree of variation to the fixed approach, which reveals that the noise-like pattern is not from the variation of wave facets. The images of Band 4 and Band 5 support this interpretation because the spatial pattern of the reflectance in the two bands are consistent with each other and it reflects the wave facet distribution (note that the two NIR bands have R rs values of almost zero; therefore, the residual surface reflectance mostly contributes to the reflectance of the bands). The spatial variation of Chla and NIR R rs appears to have no clear spatial correlation, implying that the high Chla variation in the images is not caused by the residual surface reflectance (equivalently, the viewing geometry) but from a factor related to the image quality or pixel registration. A regression analysis was performed on a further subset area (40 × 40 pixels) of the scene. Figure 7 shows the scatter plots of the pixel-to-pixel water reflectance of four bands (Bands 1 and 3-5), with respect to the reference band, Band 2. In all bands, a general linear relationship was identified; however, it showed a low correlation ( < 0.75). To assess the spatial pattern of the misregistration, the Band 2 image was regressed to Band 1, using the slope and offset estimated in the regression analysis. The difference between the original Band 1 reflectance and the regressed Band 1 image (Figure 7b) clearly shows a pixel-wise mismatch and the spatial patterns differ from that of the wave facets. It can be observed that the Chla estimation from the mismatch data ( Figure  7b) shows a similar spatial variation to that of the difference image. A regression analysis was performed on a further subset area (40 × 40 pixels) of the scene. Figure 7 shows the scatter plots of the pixel-to-pixel water reflectance of four bands (Bands 1 and 3-5), with respect to the reference band, Band 2. In all bands, a general linear relationship was identified; however, it showed a low correlation (R 2 < 0.75). To assess the spatial pattern of the misregistration, the Band 2 image was regressed to Band 1, using the slope and offset estimated in the regression analysis. The difference between the original Band 1 reflectance and the regressed Band 1 image (Figure 7b) clearly shows a pixel-wise mismatch and the spatial patterns differ from that of the wave facets. It can be observed that the Chla estimation from the mismatch data ( Figure 7b) shows a similar spatial variation to that of the difference image. Remote Sens. 2020, x, x FOR PEER REVIEW 10 of 20 (a) (b) Figure 7. (a) Regression results between the water reflectance for multispectral bands, with respect to Band 2, and (b) water reflectance for a 40 × 40 pixel subset area of Band 1, Band 1 estimates regressed from Band 2, the difference, and Chla estimates from the two bands.

Proposed Approach for the Morphological Registration
Because the surface of the water is not a rigid body, no geometric transformation can find its appropriate pixel-to-pixel correspondence. For a solution, the overall reflectance distribution must be conserved in a sufficiently small area (referred to as "window" hereafter), even if we do not know which pixel in a window corresponds to a pixel in the other band. The images of all five bands contain five instances of the scene at slightly different timings. Consequently, it can be assumed that the distribution of physical quantity, such as the radiance, does not significantly vary in the short period. By setting the image acquisition time of one band as a reference time frame, the radiometric values of the other bands can be matched to the reference band, according to the radiometric distribution, not the pixel location. Figure 8 shows the quantile-to-quantile plot (QQ plot) of the four bands, with respect to Band 2, exhibiting that the radiometric relationship can be established with a nearly perfect correlation when the reflectance of the two bands are compared based on the quantile in reflectance, not on the pixel location. A comparison of the QQ plot with the previous pixel-to-pixel scattered plots (Figure 7) reveals how the pixel-to-pixel misregistration, based on the location, degraded the correlation between the bands. In all four cases of the QQ plots, is nearly one, producing band-dependent slopes and offsets for the linear relationship. Using this new linear model, the Band 2 image was regressed to Band 1 and compared with the previous results from the regression analysis. Figure 9 shows that Band 1, regressed from the QQ plot, exhibits reflectance levels that are more similar to the original Band 1 image than the location-based regression case. The application of the linear models, derived by the QQ plots, to all four bands, serves as the morphological registration between bands and the subsequent Chla estimation produces a significant improvement in the image quality ( Figure 10). Figure 10 shows that the noise pattern for Chla in the original data was significantly reduced and the corrected Chla image contains only the spatial variation caused by wave facets, without being affected by the pixel-to-pixel misregistration. The mean and median values are comparable between the two results; however, the standard deviation and the coefficient of variation reduce from 1.03 to 0.28 and from 29% to 8%, respectively.

Proposed Approach for the Morphological Registration
Because the surface of the water is not a rigid body, no geometric transformation can find its appropriate pixel-to-pixel correspondence. For a solution, the overall reflectance distribution must be conserved in a sufficiently small area (referred to as "window" hereafter), even if we do not know which pixel in a window corresponds to a pixel in the other band. The images of all five bands contain five instances of the scene at slightly different timings. Consequently, it can be assumed that the distribution of physical quantity, such as the radiance, does not significantly vary in the short period. By setting the image acquisition time of one band as a reference time frame, the radiometric values of the other bands can be matched to the reference band, according to the radiometric distribution, not the pixel location. Figure 8 shows the quantile-to-quantile plot (QQ plot) of the four bands, with respect to Band 2, exhibiting that the radiometric relationship can be established with a nearly perfect correlation when the reflectance of the two bands are compared based on the quantile in reflectance, not on the pixel location. A comparison of the QQ plot with the previous pixel-to-pixel scattered plots (Figure 7) reveals how the pixel-to-pixel misregistration, based on the location, degraded the correlation between the bands. In all four cases of the QQ plots, R 2 is nearly one, producing band-dependent slopes and offsets for the linear relationship. Using this new linear model, the Band 2 image was regressed to Band 1 and compared with the previous results from the regression analysis. Figure 9 shows that Band 1, regressed from the QQ plot, exhibits reflectance levels that are more similar to the original Band 1 image than the location-based regression case. The application of the linear models, derived by the QQ plots, to all four bands, serves as the morphological registration between bands and the subsequent Chla estimation produces a significant improvement in the image quality ( Figure 10). Figure 10 shows that the noise pattern for Chla in the original data was significantly reduced and the corrected Chla image contains only the spatial variation caused by wave facets, without being affected by the pixel-to-pixel misregistration. The mean and median values are comparable between the two results; however, the standard deviation and the coefficient of variation reduce from 1.03 to 0.28 and from 29% to 8%, respectively. Remote Sens. 2020, x, x FOR PEER REVIEW 11 of 20    Because the determination of the window size can be critical for this approach, the sensitivity of the window size has been investigated. Six different window sizes-15, 25, 51, 101, 251, and 501 pixels-were tested for the QQ plots between Bands 1 and 2 ( Figure 11). The window areas for the various sizes were displayed in the corrected Band 1 reflectance image. The plots showed that a high correlation between Bands 1 and 2 was maintained throughout all window sizes; however, the derived linear relationships were different for different window sizes. As the window size increased, the slopes increased and the y-intercepts decreased. This reveals that the reflectance distribution may change depending on the areas used for the QQ calculation and the local characteristics (in a small window) may be lost when the QQ is derived for a large area. To achieve the goal of the proposed morphological registration, the window size must be kept as small as the local distribution because fetching the reflectance value from a distant location may not guarantee that the two values are from the continuum of targets.
Remote Sens. 2020, x, x FOR PEER REVIEW 12 of 20 pixels-were tested for the QQ plots between Bands 1 and 2 ( Figure 11). The window areas for the various sizes were displayed in the corrected Band 1 reflectance image. The plots showed that a high correlation between Bands 1 and 2 was maintained throughout all window sizes; however, the derived linear relationships were different for different window sizes. As the window size increased, the slopes increased and the y-intercepts decreased. This reveals that the reflectance distribution may change depending on the areas used for the QQ calculation and the local characteristics (in a small window) may be lost when the QQ is derived for a large area. To achieve the goal of the proposed morphological registration, the window size must be kept as small as the local distribution because fetching the reflectance value from a distant location may not guarantee that the two values are from the continuum of targets.

Algorithm Development for the Entire Image
The QQ plot approach is iterated over the image dimension to process the entire image. However, the calculation of quantiles, which is essentially an order statistics, for all pixels requires exhaustive computation with the complexity of ( • ). A Rededge-M image consists of 1280 ×

Algorithm Development for the Entire Image
The QQ plot approach is iterated over the image dimension to process the entire image. However, the calculation of quantiles, which is essentially an order statistics, for all pixels requires exhaustive computation with the complexity of O(n·log n). A Rededge-M image consists of 1280 × 960 pixels, which totals to~1.2 × 10 6 pixels. Thus, we employ an alternate fast approach, where the QQ calculation is performed for subsampled pixels (e.g. every n-th pixel), and the resultant linear model coefficients (i.e., slope and y-intercept) are propagated to the vicinity of the subsampled pixels with distanced weights assigned by a 2-dimensional Gaussian filter. Figure 12a,b displays the slope and the y-intercept images that were calculated at every pixel, which were then compared with the images obtained with the subsampling scheme, involving the Gaussian filers (Figure 12c,d) (the window size of the Gaussian filter (w gauss ) was 25, and the step size (d step ) was 12). While the spatial patterns do not significantly deviate from each other, the computation time scales down from 12 min to 1 min per band, saving more than 90% of the computation time (computation done with Intel®Core TM i-5-8265U CPU@1.60GHz, and 8GB RAM). The overall flow chart of the algorithm is presented in Figure 13.
Remote Sens. 2020, x, x FOR PEER REVIEW 13 of 20 960 pixels, which totals to ~1.2 × 10 6 pixels. Thus, we employ an alternate fast approach, where the QQ calculation is performed for subsampled pixels (e.g. every n-th pixel), and the resultant linear model coefficients (i.e., slope and y-intercept) are propagated to the vicinity of the subsampled pixels with distanced weights assigned by a 2-dimensional Gaussian filter. Figure 12a,b displays the slope and the y-intercept images that were calculated at every pixel, which were then compared with the images obtained with the subsampling scheme, involving the Gaussian filers (Figure 12c,d) (the window size of the Gaussian filter ( ) was 25, and the step size ( ) was 12). While the spatial patterns do not significantly deviate from each other, the computation time scales down from 12 min to 1 min per band, saving more than 90% of the computation time (computation done with Intel® Core TM i-5-8265U CPU@1.60GHz, and 8GB RAM). The overall flow chart of the algorithm is presented in Figure 13.

Results
The proposed algorithm was applied to the four data sets listed in Table 1. Figure 14 shows the comparison between the Chla estimates, before and after the application of the algorithm for Scene-A, taken at an altitude of 85.4 m. The high-frequency Chla variation before the correction was significantly and consistently reduced after the correction, over all areas of the image, enhancing the sharpness of the wave features on the surface. Note that the extremely high Chla values under the ship is caused by the ship shadows that made the blue-to-green ratio significantly altered

Results
The proposed algorithm was applied to the four data sets listed in Table 1. Figure 14 shows the comparison between the Chla estimates, before and after the application of the algorithm for Scene-A, taken at an altitude of 85.4 m. The high-frequency Chla variation before the correction was significantly and consistently reduced after the correction, over all areas of the image, enhancing the sharpness of the wave features on the surface. Note that the extremely high Chla values under the ship is caused by the ship shadows that made the blue-to-green ratio significantly altered compared to the sun shed areas, thus the anomalously high Chla values are not artifacts of the proposed morphological registration method.
occupied with fairly homogeneous radiance values. However, even if the surface features have a significantly smaller scale than the low-altitude images, the noisy pattern is still strong in all areas of the image as well as in the boxed area, which was effectively removed by the morphological band registration. In the boxed area, while the mean and the median values stay similar (mean: 2.248 to 2.229, media: 2.225 to 2.224), the scatter statics greatly improved (s.d.: 0.243 to 0.0066, c.v.: 11% to 3%) (Figure 16b). The proposed algorithm was also tested on a scene with an event in a part of the image (Figure 17). Scene-D was acquired for a strong red tide outbreak of Cochlodinium polykrikoides (the species was confirmed from microscopic analysis during the field campaign) and the Chla contents were estimated with the RBR algorithm. The morphological registration reduced the number of pixels that had extremely high Chla estimates (> 50 mg/m 3 ), which are considered to be generated misregistration pixels. While the overall extent and concentration level in the morphological registration is comparable to that of the 32 × 32 filter results, a spatial pattern of a higher frequency is still observable in the improved result.  The evaluation of images captured at various altitudes is important because the effect of misregistration may vary with the spatial frequency of surface reflectance features. Data sets from 8.1 m (Scene-B) and 196 m (Scene-C) were tested and the results are displayed in Figures 15 and 16, respectively. The figures display four types of Chla estimates, each of which is derived from (1) R rs data before the correction, (2) R rs data after low pass filtering with a 2 × 2 average window, (3) R rs data after low pass filtering with a 32 × 32 average window, and (4) R rs data after correction using the proposed morphological registration. Figure 15a shows the results for Scene-B (altitude 8.1 m), where the Chla estimates before the correction exhibits a very noisy spatial pattern, which remains even after the 2 × 2 mean filter. The noisy pattern was not minimized until the size of the mean filter increased to 32 × 32, as shown in the figure. The Chla estimates, with the morphological registration, exhibited a noise-free retrieval while maintaining the sharpness of the image. Note that the 32 × 32 mean filter removed the noise at the expense of losing spatial details, or sharpness. For a more detailed evaluation, a boxed area, marked in the figure, is displayed in Figure 15b. The figure demonstrates that the misregistered pixels caused many spikes with Chla estimates exceeding 3.0 mg/m 3 in the original resolution. In the 2 × 2 mean filter case, the Chla values of the corrected estimates are in the level of approximately 2.0 mg/m 3 , with the lowest value at approximately 1.0 mg/m 3 . It is neither realistic that Chla varies with a factor of three in such a small area (< 1 m 2 ) nor true that the viewing or reflecting geometry changes with such a high frequency. The scatter statics such as standard deviation and coefficient of variation decreased significantly after the morphological registration compared to the 2 × 2 mean filter (s.d.: 601 to 0.2, c.v.: 14224% to 15%), while the median value stays similar (1.745 to 1.768), implying that the outliers caused by misregistration significantly deteriorated the image quality. It is shown that the mean values were also significantly affected by the noise (4.228 to 1.745).   The results for the scene of much higher altitude (Scene-C (196 m)) are presented in Figure 16. It shows that the scene is in general more homogeneous than the previous low-altitude image (Scene-B (8.1 m)), and it exhibits a gradual radiance change across the scene with the boxed area occupied with fairly homogeneous radiance values. However, even if the surface features have a significantly smaller scale than the low-altitude images, the noisy pattern is still strong in all areas of the image as well as in the boxed area, which was effectively removed by the morphological band registration. In the boxed area, while the mean and the median values stay similar (mean: 2.248 to 2.229, media: 2.225 to 2.224), the scatter statics greatly improved (s.d.: 0.243 to 0.0066, c.v.: 11% to 3%) (Figure 16b). The proposed algorithm was also tested on a scene with an event in a part of the image (Figure 17). Scene-D was acquired for a strong red tide outbreak of Cochlodinium polykrikoides (the species was confirmed from microscopic analysis during the field campaign) and the Chla contents were estimated with the RBR algorithm. The morphological registration reduced the number of pixels that had extremely high Chla estimates (> 50 mg/m 3 ), which are considered to be generated misregistration pixels. While the overall extent and concentration level in the morphological registration is comparable to that of the 32 × 32 filter results, a spatial pattern of a higher frequency is still observable in the improved result.

Discussion and Conclusions
This study analyzed the impact of residual misregistration between the bands of a multispectral camera, which is particularly problematic for fast-moving targets such as the water surface. The analysis suggested that the residual misregistration is difficult to correct using geometric coordinate transformation (e.g., projective transformation). It produces abnormal band

Discussion and Conclusions
This study analyzed the impact of residual misregistration between the bands of a multispectral camera, which is particularly problematic for fast-moving targets such as the water surface. The analysis suggested that the residual misregistration is difficult to correct using geometric coordinate transformation (e.g., projective transformation). It produces abnormal band ratios and band differences, resulting in significant anomalies in the water quality variables, such as the Chla concentration. The proposed registration algorithm succeeded in effectively removing noisy spatial patterns caused by the misregistration while maintaining the original spatial resolution of the image, unlike the smoothing approach which significantly degrades the sharpness of the images. Contrary to the intuition that high-altitude images will be less affected by pixel-level misregistration (because a scene appears more homogeneous when observed from far), the test results for various altitudes showed that the residual misregistration exist at all tested altitudes (8-390 m). This suggests that there exist various frequencies of the surface reflectance feature on the water surface. The proposed algorithm was robust to local events that occurred in a partial section of the image, which may have distinct spectral characteristics of the remaining image area (usually normal water surface), as shown in the red tide image.
The proposed method is expected to improve estimation of other water quality variables such as colored dissolved organic matter (CDOM) and total suspended sediments (TSM) as many of the CDOM and TSM retrieval algorithms rely on the band ratio of reflectance.
Future work includes analysis on the effects on other water quality variables, and also the further correction of residual sky reflectance that is often caused by the spatially varying normals of the wave facets. The residual sky reflectance still existed even after the morphological band registration. The research requires a comprehensive understanding of the reflecting mechanism and more detailed modeling of the surface reflectance, associated with the analysis on the downward sky radiance incident on the water surface. Funding: Funding: This work was performed with the support of Jeollannam-do through the research project titled "Development of a monitoring system for damage reduction in aquacultures using marine environment observation drones" (B0080802000131) operated by the Jeonnam Technopark (JNTP), and was also supported by the National Research Foundation (NRF) project "Developing an atmospheric correction module for retrieving water quality from hyperspectral images" (NRF-2019R1F1A1062585).