# Registration of close-range, multi-lens multispectral imagery by retrieving the scene 3D structure

Sylvain Jay, Frédéric Baret, Samuel Thomas, Marie Weiss

► To cite this version:

Sylvain Jay, Frédéric Baret, Samuel Thomas, Marie Weiss. Registration of close-range, multi-lens multispectral imagery by retrieving the scene 3D structure. ISPRS Journal of Photogrammetry and Remote Sensing, 2025, 227, pp.125-144. 10.1016/j.isprsjprs.2025.06.001. hal-05126456

HAL Id: hal-05126456

https://hal.inrae.fr/hal-05126456v1

Submitted on 23 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Registration of close-range, multi-lens multispectral imagery by retrieving the scene 3D structure

Sylvain Jay<sup>a, *</sup>, Frédéric Baret<sup>a</sup>, Samuel Thomas<sup>b</sup>, Marie Weiss<sup>a</sup>

## ARTICLE INFO

**Keywords:**

Image registration

Multispectral

Close range

Stereovision

3D point cloud

Vegetation

Crop phenotyping

## ABSTRACT

Multispectral, multi-lens cameras, which acquire spectral images from different individual cameras equipped with different optical filters, are among the most widely used multispectral cameras available on the market. However, their use for close-range sensing is limited by the lack of registration algorithms capable of handling the strong parallax effects observed on scenes with non-negligible relief. In this paper, we propose a method based on stereo camera calibration and disparity estimation to register a close-range multispectral image while retrieving the corresponding 3D point cloud. The method takes advantage of the rigidity of these cameras and the synchronized capture of multispectral bands, both of which are thus compulsory. The algorithm is three-fold. First, the optimal combination of band pair alignments is found. Then, the semi-global matching stereovision algorithm combined with a robust matching cost function are used to align these band pairs and to compute the point cloud. Finally, a pixel filling step that exploits the spectral covariances of the different classes of materials in the image is implemented to limit the number of missing pixels, e.g., due to occlusions.

The method was tested on Airphen multispectral images of four plant crops (wheat, sunflower, cover crops and maize) acquired at a distance to the ground ranging from 1.5 to 3 m, thus encompassing a large variability in 3D structure and parallax effects. The results demonstrate that the proposed method achieves better registration performance than six state-of-the-art existing methods, while maintaining a reasonable processing time. Further, the point cloud provides accurate information on the 3D structure of the imaged scene, as shown by the centimetric plant height estimation accuracy. As the point cloud is aligned with the registered multispectral bands, the method provides a 4D (spectral and spatial) description of the scene with a single image, i.e., a multispectral point cloud. This opens up interesting prospects for several applications in close-range sensing including, but not restricted to, vegetation characterization.

## 1. Introduction

Multispectral imagery is widely used in remote and proximal sensing, as it offers a good compromise in terms of spatial resolution, spectral information, ease of use and cost as compared to RGB (red–green–blue) or hyperspectral imagery (Berger et al., 2022). For example, such spectral information can be used to monitor phytoplankton based on satellite kilometer-scale observations (O'Reilly et al., 1998), to characterize crop biophysical properties based on unmanned aerial vehicle (UAV) centimeter-scale observations (Jay et al., 2019), to assess degradation of historic monuments (Adamopoulos and Rinaudo, 2021a; Lerma et al., 2012) or to discriminate sedimentary rocks (Pozo

et al., 2015) based on close-range images. In addition, 3D information like plant height (Li et al., 2021; Lu et al., 2022; Madec et al., 2017) can also be retrieved using Structure from Motion techniques applied on multiple overlapping UAV multispectral images. When acquired at one to few meters from the target, most multispectral sensors provide millimeter-scale images that may be used to characterize biophysical properties of organs in an individual plant for phenotyping or precision agriculture.

The different multispectral sensors available on the market are usually based either on single-lens cameras or multi-lens cameras (Jhan et al., 2018). Single-lens cameras include RGB cameras with modified optical filters to get the near-infrared (NIR) signal (Hunt et al., 2010;

<footnote><sup>a</sup> INRAE UMR 1114 EMMAH, Domaine Saint-Paul, Site Agroparc, F-84914 Avignon, France</footnote>

<footnote><sup>b</sup> Arvalis Institut du végétal, F-84914 Avignon, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: sylvain.jay@inrae.fr (S. Jay).</footnote>
Rabatel et al., 2014), snapshot multispectral cameras based on a custom Bayer-like mosaic filter and a classical CMOS sensor,<fnref n="1" /> or cameras based on Fabry-Pérot interferometry (Honkavaara et al., 2013). Multi-lens cameras are also very common, and include the Airphen,<fnref n="2" /> Parrot Sequoia<fnref n="3" /> or MicaSense RedEdge-P<fnref n="4" /> and Altum-PT<fnref n="5" /> cameras. In this case, each spectral band is acquired from a slightly different point of view corresponding to a different optical center. Therefore, images in different spectral bands must be registered before further multispectral image analysis. As compared to images of the same nature, it is necessary to account for the spectral signature of the target, which might be quite contrasted between one band and another, e.g., vegetation in the visible and NIR domains.

In the computer vision and remote sensing communities, many methods have been proposed to register multispectral images under the assumption of planar scene, i.e., when the distance between the camera and the scene is assumed to be much greater than the depth of the imaged scene (Cao et al., 2020; Chen et al., 2018; Fernández et al., 2021; Kerkech et al., 2020; Meng et al., 2021; Rabatel and Labbé, 2016; Shen et al., 2019). When three bands or more have to be registered, the first step is to find a combination of N<sub>b</sub>-1 band pair alignments (where N<sub>b</sub> is the number of bands). For example, N<sub>b</sub>-1 bands can be registered with the remaining reference band, or more complex combinations based on intermediate reference band(s) can be used to facilitate the registration process (Dias Junior et al., 2020; Vayssade et al., 2020; Yasir et al., 2018). This combination of N<sub>b</sub>-1 band pairs can be set empirically (Banerjee et al., 2018; Cao et al., 2020; Chen et al., 2018; Dandrifosse et al., 2021; Fernández et al., 2021; Laveglia et al., 2025; Meng et al., 2021; Rana et al., 2024; Shen et al., 2019) or optimized by maximizing the total number of matched feature points over the band pairs (Dias Junior et al., 2020; Vayssade et al., 2020; Yasir et al., 2018). Then, for every band pair to be registered, these methods estimate a two-dimensional (2D) affine or homographic transformation that relate pixel coordinates in the fixed image (the reference, also called left image in the literature) and pixel coordinates in the moving image (the image to be registered with the reference, also called right image in the literature). This transformation is retrieved using either feature- (Banerjee et al., 2018; Dias Junior et al., 2020; Kerkech et al., 2020; Laveglia et al., 2025; Lowe, 2004; Rana et al., 2024; Vayssade et al., 2020; Yasir et al., 2018; Ye and Shan, 2014) or area-based (Cao et al., 2020; Chen et al., 2018; Douarre et al., 2019; Sharma et al., 2023; Shen et al., 2019; Wang et al., 2010) methods (Zitova and Flusser, 2003) that are robust to radiometric variations across bands.

While the above methods usually perform well for UAV and satellite images, the assumption of planar scene that they make is violated when the depth of the imaged scene and the distance between the camera and the scene are of the same order of magnitude, e.g., when imaging plants at one to few meters distance. In this case, 2D affine or homographic transformations are not sufficient to model parallax effects, which vary within the image as a function of the pixel depth (Debaque et al., 2022; Krotosky and Trivedi, 2007; Liu et al., 2018; Szeliski, 2010). More complex representations of the motion field should then be used to deal with close-range multispectral images, such as parametric spline-based models or more general optical flow methods (Szeliski, 2010). Widely used for medical image registration, spline-based models represent small pixel displacements using a linear combination of 2D spline basis functions regularly distributed within the image. Additionally, they are able to cope with radiometric variations when using mutual information as similarity measure (Rueckert et al., 1999). These local deformation

models are often used to refine an initial alignment obtained with a global affine or homographic transformation. For example, Dandrifosse et al. (2021) compared different strategies to register close-range multimodal images of wheat crops, and found that a combination of a depth-dependent homography and a cubic B-spline model based on mutual information provided the best results. Unlike 2D planar transformations and spline-based models, optical flow methods are non-parametric and estimate an independent pixel displacement for each pixel, e.g., between two adjacent frames in a video sequence (Szeliski, 2010). For example, Barjaktarovic et al. (2024) showed the interest of dense optical flow to refine the coarse registration results obtained with a phase correlation technique. While classic optical flow methods establish dense, pixel-to-pixel correspondences based on raw pixel values, the SIFT flow algorithm (Liu et al., 2011) efficiently matches Scale-Invariant Feature Transform (SIFT) descriptors (Lowe, 2004) computed at each pixel. Inspired from the SIFT flow method, Kim et al. (2017) proposed to replace the SIFT descriptor with their dense adaptive self-correlation (DASC) descriptor in the SIFT flow optimization framework to better deal with multispectral and multimodal images. They showed that the combination of DASC features and SIFT flow optimization could handle complex radiometric and geometric variations between images. However, one disadvantage of spline-based and optical flow models is that several regularization terms must be added to the main loss function to better constrain the optimization problem, like for example a smoothness term on the motion field (Liu et al., 2011; Rueckert et al., 1999). This regularization process may induce registration problems near discontinuities (Szeliski, 2010) and require additional post-processing such as pixel removal using image erosion (Dandrifosse et al., 2021).

Whenever possible, another way to constrain the registration problem is to exploit the rigidity of multi-lens cameras and the synchronized capture of multispectral bands. For any band pair combination, optical flow model retrieval can then be transformed into a simpler stereo matching problem involving the retrieval of the cameras' extrinsic parameters, which are the relative positions and orientations of the two cameras corresponding to the two spectral bands, and a per-pixel estimation of disparity, which is the distance between the projections of the same object point on the two band images (Szeliski, 2010; Tosi et al., 2022). Camera extrinsics can be first estimated through a calibration experiment using multiple images of a chessboard, and used to rectify new image pairs, i.e., to warp the images in order to make corresponding epipolar lines coincident. This improves the matching process since the 2D search is transformed in a 1D one (Szeliski, 2010), and allows the conversion of disparity values (in pixels) to distance values (e.g., in metric unit). This is particularly useful to characterize the 3D structure of the imaged scene in a single multispectral image capture.

Disparity relating two spectral bands can be estimated using either traditional stereo matching methods or deep learning based methods (Guo et al., 2019). On the one hand, traditional stereo pipelines are composed of four steps, i.e., matching cost computation, cost aggregation, disparity optimization and disparity refinement (Scharstein and Szeliski, 2002). However, most studies dealing with cross-spectral stereo matching only focused on designing an as robust as possible matching cost that can mitigate differences in intensity values between the two bands. Examples of such robust matching costs include the mutual information (MI) (Heo et al., 2013; Hirschmüller and Scharstein, 2007; Hirschmuller and Scharstein, 2009; Krotosky and Trivedi, 2007), the L1 norm between dense Histograms of Oriented Gradient features (Pinggera et al., 2012), the census transform combined with Hamming distance (Hirschmuller and Scharstein, 2009) and the normalized total gradient (NTG) (Chen et al., 2018). On the other hand, more recent deep learning based stereo methods exploit one or several dedicated deep neural network(s) trained with real or synthetic datasets. Zbontar and LeCun (2016) first proposed to compute the matching cost between two small RGB image patches using a Siamese network, consisting of two shared-weight convolutional neural networks (CNN) joined at the head.

<footnote n="1">¹ https://www.silios.com/toucan-camera.</footnote>

<footnote n="2">² https://www.hiphen-plant.com.</footnote>

<footnote n="3">³ https://www.parrot.com/fr/shop/accessoires-pieces-detachees/autres-drones/sequoia.</footnote>

<footnote n="4">⁴ https://ageagle.com/drone-sensors/rededge-p/.</footnote>

<footnote n="5">⁵ https://ageagle.com/drone-sensors/altum-pt-camera/.</footnote>
Zhao et al. (2017) adapted this approach to cross-spectral stereo matching by training the same Siamese network on heterogeneous pairs of red, green and blue channels extracted from the training images used by Zbontar and LeCun (2016). Beaupre and Bilodeau (2021) proposed to extract features using a Siamese network without shared weights to better account for radiometric differences between spectral bands. Alternatively, unsupervised methods have been developed due to the lack of annotated cross-spectral datasets with ground-truth disparity (Duan and Jung, 2022). They generally combine a spectral translation network that makes the two spectral images more similar (e.g., it generates a NIR image based on a RGB image), with a disparity prediction network that ensures disparity smoothness as well as left-right disparity image consistency and left-right spectral image consistency after warping using corresponding disparity images (Chang et al., 2022; Duan and Jung, 2022; Liang et al., 2019; Zhi et al., 2018). Deep learning based stereo matching methods now outperform their traditional counterparts (Cournet et al., 2020), as illustrated by the method intercomparison results regularly updated on the Middlebury stereovision page<fnref n="6" /> (Scharstein and Szeliski, 2002). However, they require huge datasets for training (with or without annotations) that are difficult to obtain for multispectral imagery, e.g., due to the differences in spectral configuration across camera brands which implies the constitution of a training dataset per camera configuration (Tosi et al., 2022). This may be overcome by simulating multispectral bands based on publicly available RGB stereo datasets (Genser et al., 2020), but simulated imagery might not easily transfer to the true cross-spectral case (Pinggera et al., 2012).

In this study, we propose a calibration-based method that embeds a traditional pairwise, cross-spectral stereo matching method within a three-step framework, to register a multispectral image with three or more bands. This method is automatic, can be applied to any multi-lens, multispectral camera directly (i.e., without constituting a training dataset beforehand), and is especially dedicated to register multispectral images acquired over a scene with non-negligible relief. The key points of our method are (1) the traditional semi-global matching (SGM) stereo algorithm (Hirschmüller, 2008) using census transform as matching cost for pairwise band registration and pixelwise depth estimation, (2) the Normalized Total Gradient to quantify the pairwise registration accuracy, (3) an optimal combination of N<sub>b</sub>-1 band pairs according to which the N<sub>b</sub> bands are registered, and (4) a pixel filling step to limit the number of invalid pixels, e.g., due to occlusions. In Section 2, we present our method as well as the field data and the six state-of-the-art methods that are used for comparative performance assessment. Since plant phenotyping is one of the main applications of close-range multispectral imagery, these methods were tested on images of four contrasted crops (wheat, sunflower, cover crops and maize) acquired at 1.5 m to 3 m distance from the ground level. The results obtained with the different methods in terms of registration accuracy, depth estimation accuracy and processing time are compared in Section 3, and discussed in Section 4. We finally draw some conclusions and provide some perspectives in Section 5.

## 2. Materials and methods

### 2.1. Proposed algorithm

The proposed algorithm<fnref n="7" /> consists of three steps. In step 1, the combination of band pair alignments that optimizes the overall registration accuracy is found (section 2.1.1). In step 2, bands are aligned thanks to the optimal combination determined in step 1, and the point cloud is computed (section 2.1.2). Step 3 allows the filling of missing band pixels using regression models (section 2.1.3).

#### 2.1.1. Step 1: Finding the optimal combination of band pair alignments

The algorithm starts by aligning each possible band pair using corresponding disparity (Fig. 1). Note that bands are here subsampled by a factor of two for faster processing. Each band pair alignment is performed similarly (Fig. 2). The fixed (also called left) and moving (also called right) images are first rectified using parameters obtained from prior stereo calibration so that the projections of the same object point are located on the same row on the fixed and moving images (section 2.2). By definition, disparity represents the remaining horizontal distance between the projections of that same object point on both images. Pixel (i, j) in the moving image, noted I<sub>m</sub>(i, j), can then be aligned with its counterpart in the fixed image by shifting it horizontally by D(i, j) columns, where D(i, j) is the disparity value in pixel (i, j):

$$I_{m,aligned}(i,j) = I_m(i,j - D(i,j)) \qquad (1)$$

where I<sub>m,aligned</sub>(i, j) is the value of aligned moving image at row i and column j. In this study, disparity image is computed using the SGM stereovision algorithm (Hirschmüller, 2008), which is one of the top-ranked traditional stereo algorithms (Cournet et al., 2020). As explained earlier, the census transform combined with Hamming distance (Zabih and Woodfill, 1994) is used as matching cost within the SGM algorithm to decrease the influence of radiometric differences between bands (Hirschmüller and Scharstein, 2009). After computing disparity image D and aligned moving image I<sub>m,aligned</sub> (Eq. 1), both are finally warped using inverse rectification transforms obtained from prior calibration (section 2.2).

The registration accuracy of each possible band pair is then evaluated by computing the NTG value with the fixed image and the aligned moving image, because of its robustness against radiometric differences (Chen et al., 2018). The NTG value is computed from the L1-norms of image gradients and relies on the assumption that the gradient of the difference between two aligned images is sparser than that between two misaligned ones (Chen et al., 2018). The NTG tends towards zero when the two images are perfectly aligned. A NTG matrix is therefore calculated, where the NTG value at row i and column j is obtained by using the i<sup>th</sup> band as fixed image and the j<sup>th</sup> band as moving image (Fig. 1).

Finally, we select the optimal combination of N<sub>b</sub> – 1 band pair alignments among all candidates that respect the following constraints (Fig. 1):

- Each combination must include N<sub>b</sub> – 1 pairs of a fixed band and a moving band;

- Each of the N<sub>b</sub> bands must be present in at least one of these pairs;

- Each combination must have a master band with which every other band will be aligned; in the pairs where it is present, the master band is always the fixed image;

- Each combination must be of depth one (i.e., N<sub>b</sub> – 1 moving bands are aligned with the remaining, fixed master band) or of depth two (i.e., first, 1 to N<sub>b</sub> – 2 moving band(s) are aligned with the remaining non-master band(s), here considered as fixed images, and second, these N<sub>b</sub> – 1 non-master bands, now considered as moving images, are aligned with the fixed master band);

- Each band pair alignment is characterized by the NTG value computed with the fixed band and the aligned moving band.

Examples of such combinations of N<sub>b</sub> – 1 band pair alignments are shown in Fig. 1. The depth of these combinations of band pair alignments is lower than or equal to two, which facilitates the registration of strongly different bands thanks to the use of an intermediate band (Dias Junior et al., 2020; Vayssade et al., 2020; Yasir et al., 2018; Zhao et al., 2017), while limiting the numbers of occluded pixels (as detected by the SGM algorithm) that add up when aligning several bands on top of each other.

Among all candidate band pair combinations (e.g., a total of 1176 combinations if N<sub>b</sub> = 6 bands), we select the one that has both low

<footnote n="6">⁶ https://vision.middlebury.edu/stereo/.</footnote>

<footnote n="7">⁷ Source code is available here: https://forgemia.inra.fr/sylvain.jay/close_range_multispectral Registration.git</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Flowchart describing how to find the optimal combination of band pair alignments based on Normalized Total Gradient (NTG) (step 1 of proposed algorithm). A multispectral image of maize plants acquired with a six-band Airphen camera is used in this example (section 2.3). Stereo calibration parameters here include rectification transforms and extrinsic camera parameters obtained from prior calibration. Each arrow for each represented band pair goes from the moving image to the fixed image, and is associated with the corresponding NTG value computed with the fixed band and aligned moving band. For a given combination of five band pairs, $\overline{NTG}_c$ is the average NTG value computed over these five band pairs.</figcaption></figure>

average NTG value (i.e., high registration accuracy) and low average distance between band optical centers (leading to fewer occluded areas), the two averages being computed over the $N_b$ – 1 band pairs (Fig. 1). More specifically, we first identify the best band pair combinations with average NTG values falling within the $[NTG_{min}; NTG_{min} + \Delta NTG]$ interval, where $NTG_{min}$ is the minimum average NTG value obtained over all possible combinations and $\Delta NTG$ is empirically set to 5 % of the range of possible average NTG values ($\Delta NTG \approx 0.01$ in practice). Among these few best band pair combinations, we finally select the one with minimum average distance between band optical centers, these distances being obtained from prior stereo calibration (section 2.2).

### 2.1.2. Step 2: Aligning bands and computing point cloud

In step 2, each of the optimal $N_b$ – 1 band pairs determined in step 1 (section 2.1.1) is aligned using corresponding disparity image (Fig. 3). Unlike step 1, full resolution spectral bands are used here. If the optimal combination of band pair alignments is of depth two, the 1 to $N_b$ – 2 depth-two band(s) is(are) only aligned with its(their) corresponding depth-one, intermediate band(s) and associated disparity image(s). Then, the disparity image(s) already used to register this(these) depth-one band(s) with the master band are used again to align the pre-aligned depth-two band(s) and associated disparity image(s) with the master band. For example, in Fig. 3, the 730–850 disparity image is first
<figure><img src="image_5.png" /><figcaption>Fig. 2. Registration process applied for each band pair. Rectification transforms are obtained from prior stereo calibration (section 2.2). In this example, the fixed image is the 730 nm band, while the moving image to be aligned is the 850 nm band.</figcaption></figure>

used to align the depth-two 850 nm band with the depth-one 730 nm band. The pre-aligned 850 nm band and the 730–850 disparity image are then aligned with the master band using the 710–730 disparity image.

This process allows us to register not only the N<sub>b</sub> spectral bands, but also the N<sub>b</sub> – 1 disparity images. These disparity images are then combined into one disparity image of better quality before conversion to 3D point cloud (Fig. 3). However, as illustrated by the differences between 710–730 and 730–850 disparity images in Fig. 3, disparity not only depends on depth, but also on focal length and baseline, i.e., the distance between the two optical centers (Szeliski, 2010):

$$D_n(i, j) = f_n^* \frac{b_n}{Z(i, j)} \qquad (2)$$

where D<sub>n</sub>(i, j) is the disparity value (in pixels) at pixel (i, j) for the n<sup>th</sup> band pair, Z(i, j) is the depth value (in unit of length) at pixel (i, j), and f<sub>n</sub> and b<sub>n</sub> are, respectively, the focal length (in pixels) and baseline (in unit of length) for the n<sup>th</sup> band pair. Before being combined into one by computing their median, the disparity images were corrected from differences in baseline and focal length. One of the N<sub>b</sub> – 1 disparity images was arbitrarily chosen as reference, with baseline b<sub>ref</sub> and focal length f<sub>ref</sub>. Eq. 2 allows us to relate two disparity images of the same 3D scene, represented by Z(i, j). Here it is used to relate the disparity image D<sub>n</sub> obtained with baseline b<sub>n</sub> and focal length f<sub>n</sub>, and a corrected disparity image D<sub>n,corr</sub> that would be obtained with baseline b<sub>ref</sub> and focal length f<sub>ref</sub>:

$$f_n^* \frac{b_n}{D_n(i, j)} = f_{ref}^* \frac{b_{ref}}{D_{n,corr}(i, j)}. \qquad (3)$$

Eq. 3 shows that a multiplicative factor can be applied to each disparity image D<sub>n</sub> to correct for differences in baseline and focal length across the N<sub>b</sub> – 1 band pairs:

$$D_{n,corr}(i, j) = D_n(i, j)^* \frac{b_{ref}^* f_{ref}}{b_n^* f_n}. \qquad (4)$$

The median disparity image is finally computed over all corrected disparity images, and converted to 3D point cloud using stereo parameters derived from prior calibration (section 2.2).

### 2.1.3. Step 3: Filling missing band values in the aligned multispectral image
After step 2, there can be a non-negligible number of missing pixels

in each band of the registered multispectral image, especially for depth-two bands that are aligned with master band using two successive disparity images (section 2.1.2). These missing pixels correspond to objects that are not seen from all band optical centers, i.e., objects outside the camera field-of-view for at least one band, and objects located behind other objects in the scene (occlusions). For a given pixel with at least one missing band value, it may be possible to predict this missing value with sufficient accuracy based on one or several combination(s) of other available band values, thanks to the covariance between spectral bands.

Step 3 starts with the calibration of clustering and regression models using randomly selected half of the image pixels with all bands available (CAL set in Fig. 4), after subsampling the images by a factor of three for faster processing (Fig. 4). As normal distributions can accurately represent pixel spectral variability for each class of material (Jay et al., 2012; Manolakis et al., 2003; Manolakis and Shaw, 2002), a Gaussian mixture model (GMM) is used to cluster the pixels into N<sub>c</sub> classes (the N<sub>c</sub> value is optimized as described in section 2.4), in which the spectral variability is reduced and the covariance between bands is larger. Then, for each class k, multiple linear regression (MLR) models relating every band to be filled (output) B<sup>k</sup>(λ<sub>i</sub>), and every possible subset of remaining bands (input) S, are calibrated (Fig. 4). Calibrating MLR models on groups of similar pixels clustered with GMM is expected to be more effective since, for example, it should be easier to predict the red band value from the other visible band values if we consider only green vegetation pixels or only soil pixels.

These models are validated using the other half of pixels with all bands available (VAL set in Fig. 4). A subset of missing band(s) S<sub>m</sub> is simulated, and the class label is predicted for each VAL pixel all using remaining available band(s) S<sub>a</sub> and a truncated version of the GMM, in which missing band(s) has(have) been removed from the N<sub>c</sub> mean vectors and N<sub>c</sub> spectral covariance matrices. Then, for each band B<sub>(λ<sub>i</sub>)</sub> ∈ S<sub>m</sub> and each class k, an estimate, B<sub>(λ<sub>i</sub>)</sub>, is computed using every possible MLR models calibrated previously for this class and taking as input a subset S of S<sub>a</sub>. The validation relative root mean square error (RMSE<sub>VAL</sub>) is computed for each subset S, and the model with lowest RMSE<sub>VAL</sub> is selected. This validation procedure is repeated for every possible subset of missing band(s).

Finally, for every pixel with at least one missing band value (TEST set in Fig. 4), the class label is predicted using all bands available, and each missing band value is estimated using the best model for this class that has been identified on the VAL set, provided that the corresponding
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Flowchart describing how to register all the spectral bands together based on the optimal band pair combination found at step 1, and stereo calibration parameters that include rectification transforms, intrinsic and extrinsic camera parameters. The multispectral image used to illustrate this process is the same as in Fig. 1. The bands presented here were chosen according to the optimal band pair combination presented in Fig. 1 (master band = 710 nm, depth-one band = 730 nm and depth-two band = 850 nm).</figcaption></figure>

relative RMSE<sub>VAL</sub> is lower than an arbitrary threshold of 10 % (otherwise, the band value is not filled).

## 2.2. Algorithm implementation

Intrinsic camera parameters (e.g., focal length and distortion coefficients) for each band, extrinsic camera parameters (translation vector and rotation matrix) and rectification transforms for each band pair were estimated using the OpenCV library (Bradski, 2000) and several dozens of multispectral images of a chessboard acquired under various angles and at distances that are similar to the ones encountered in the field, i.e., between 1 m and 4 m.

The stereo matching framework Pandora<sup>8</sup> (Cournet et al., 2020), version 1.5.0, was used to compute every disparity map. As explained by Cournet et al. (2020), Pandora is based on the work of Scharstein and Szeliski (2002) and its modularity makes it possible to define our own processing pipeline, from the choice of the matching cost to the detection of occlusions. As explained earlier, census transform combined with

Hamming distance was selected as matching cost, and SGM was selected as stereovision method. Every other Pandora parameter was set to its default value. For disparity computation in step 2 of the proposed algorithm (Fig. 3), additional steps were implemented to refine the obtained disparity image, i.e., median filtering with a window size of 5 pixels, quadratic disparity interpolation and cross-checking to detect and to mask occlusions.

A preliminary step was implemented to automatically find the range of disparity values to be explored for every pair, which is a required Pandora input. It consists in (1) selecting the pair of bands that have the most similar pixel value distributions, (2) computing the corresponding disparity image using images subsampled by a factor of three and a disparity range corresponding to a 0.5 m–5 m depth range that is much larger than the expected one, (3) extracting minimum and maximum values in this disparity image and extending this range for three pixels on both sides, and (4) converting this disparity range to a depth range using Eq. 2. For any disparity computation in steps 1 and 2 of the registration algorithm, this depth range is then transformed to a specific disparity range based on appropriate baseline and focal length values (Eq. 2).

<footnote n="8">8 https://github.com/CNES/Pandora.git.</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Flowchart describing how to fill the band values missing in the aligned multispectral image after step 2. The multispectral image used to illustrate this process is the same as in Figs. 1 and 3 ($N_b = 6$, master band = 710 nm). The CAL set contains half of the pixels with all bands available, the VAL set contains the other half, and the TEST set contains pixels with at least one missing band value. GMM and MLR refer to Gaussian mixture model and multiple linear regression, respectively.</figcaption></figure>

## 2.3. Field data used for algorithm optimization and assessment

The proposed algorithm was assessed using 40 images acquired with two multi-lens Airphen multispectral cameras with slightly different spectral configurations (Table 1). These two Airphen cameras were composed of six individual cameras equipped with filters of 10 nm spectral resolution and spaced from one to four centimeters (see (Vayssade et al., 2020) for a schematic representation). For the first camera, the six filters were centered on 450, 570, 675, 710, 730 and 850 nm, while for the second camera, they were centered on 450, 530, 570, 675, 730 and 850 nm. Both Airphen cameras were equipped with an 8 mm focal length lens and acquired 16-bit TIFF images with 1280 × 960 pixels. The integration time of each of the six individual bands was adjusted automatically to minimize saturation and maximize the dynamics.

Six datasets were collected in three sites in France (Table 1). Datasets #1 and #2 were acquired in Avignon (43°54′58″N, 4°52′44″E) in october 2020, and encompassed 20 images of a maize canopy at flowering stage with some weeds, with both green and senescent leaves. These images were acquired under natural illumination at nadir by fixing the first Airphen camera on a tripod. The maize and weed plants were imaged at 1.50 m and 1.90 m distances from the ground level, resulting in distances between the camera and the plants ranging, respectively, from 0.50 to 1.50 m, and from 0.90 to 1.90 m (Table 1). To quantify the registration and depth estimation accuracies, 40 small crosses were drawn on maize leaves using a black marker before acquiring the images, and used as control points. The height of each cross from the ground level was then measured using a ruler. This height varied between 5 cm and 72 cm, with a mean of 30 cm and a standard deviation of 17 cm. Eleven and nine images were, respectively, acquired for 1.50 m
<table><caption>**Table 1**

Overview of the six datasets used to assess the performance of the proposed registration algorithm.</caption><thead><tr><th>Dataset ID</th><th>#1</th><th>#2</th><th>#3</th><th>#4</th><th>#5</th><th>#6</th></tr></thead><tbody><tr><td>Species</td><td>Maize</td><td>Maize</td><td>Wheat</td><td>Wheat</td><td>Sunflower</td><td>Cover crops</td></tr><tr><td>Vector</td><td>Tripod</td><td>Tripod</td><td>Tripod</td><td>Phenomobile</td><td>Phenomobile</td><td>Phenomobile</td></tr><tr><td>Illumination</td><td>Passive</td><td>Passive</td><td>Passive</td><td>Active</td><td>Active</td><td>Active</td></tr><tr><td>Minimum and maximum distances between camera and scene (in m)</td><td>0.5 – 1.5</td><td>0.9 – 1.9</td><td>1.4 – 2.2</td><td>1.3 – 2.4</td><td>1.2 – 3.0</td><td>1.3 – 2.0</td></tr><tr><td>Band set (in nm)</td><td>450, 570, 675, 710, 730, 850</td><td>450, 570, 675, 710, 730, 850</td><td>450, 570, 675, 675, 730, 850</td><td>450, 570, 675, 710, 730, 850</td><td>450, 570, 675, 710, 730, 850</td><td>450, 570, 675, 710, 730, 850</td></tr><tr><td>Number of images</td><td>11</td><td>9</td><td>5</td><td>5</td><td>5</td><td>5</td></tr><tr><td>Used to...</td><td>... Compare best band pair combinations</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td></tr><tr><td></td><td>... Optimize number of classes for pixel filling</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td></tr><tr><td></td><td>... Assess registration accuracy based on global MI values</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td></tr><tr><td></td><td>... Assess registration accuracy based on local ROIs</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td></tr><tr><td></td><td>... Assess registration accuracy based on control points</td><td>x</td><td>x</td><td>—</td><td>—</td><td>—</td></tr><tr><td></td><td>... Assess height estimation</td><td>x</td><td>x</td><td>—</td><td>—</td><td>—</td></tr><tr><td></td><td>... Assess processing runtime</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td></tr></tbody></table>

and 1.90 m distances, such that (1) nearly every cross could be seen in at least one image, and (2) some crosses could be seen in several images at different locations.

The other four datasets (#3 – #6) contained 20 images of wheat, sunflower and two cover crops (faba bean and vetch) mainly acquired at flowering stage. Half of the wheat images (dataset #3) were acquired under natural illumination in Gréoux-les-bains (43°45′01″N, 5°51′17″E) in June 2019, using the second Airphen camera fixed on a tripod at 2.20 m from the ground (Table 1). The remaining images of wheat (dataset #4), sunflower (dataset #5) and cover crops (dataset #6) were acquired in Toulouse (43°31′42″N, 1°30′19″E) between May 2020 and September 2021, using the first Airphen camera embedded on the measurement

head of the Phenomobile unmanned ground vehicle (Table 1) (De Solan et al., 2015). The camera was synchronized with four Phoxene FR60 Xenon flashes with a tunable energy level ranging from 5 to 100 J, in order to make the image acquisition independent from natural illumination conditions. Thanks to a Sick LMS 400 LiDAR, the height of the measurement head was adjusted automatically to approximately 1.50 m above the average top of the canopy, resulting in distances between the camera and the ground ranging from 2 m for cover crops, to 3 m for sunflower crops (Table 1). For more details about the Phenomobile, please see De Solan et al. (2015) and Jay et al. (2020).

These six datasets allowed us to include a large variability in vegetation architecture, leaf color and image acquisition. They were used to

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Examples of Airphen images in datasets (a) #1, (b) #2, (c) #3, (d) #4, (e) #5 and (f) #6, illustrated with their 850 nm bands. Several plant species with variable canopy structures were considered: (a-b) maize and weeds, (c-d) wheat, (e) sunflower, and (f) faba bean. Images on the first row (a-c) were acquired under natural, passive illumination with a tripod, while images on the second row (d-f) were acquired under active illumination (flashes) with the Phenomobile. In (a-b), the same maize canopy was imaged at 1.50 m (a) and 1.90 m (b) from the ground level. The black small crosses used to assess registration and depth estimation accuracies can be observed on maize leaves.</figcaption></figure>
compare the best selected band pair combinations, to optimize the number of classes used in the GMM clustering algorithm for pixel filling, and to assess the processing runtime and the registration accuracy based on global MI measure and local, manually delineated regions of interest (ROIs) (Table 1, section 2.5.2). In addition, datasets #1 and #2 were used to assess the registration and depth estimation accuracies based on the small crosses of known height (Table 1). Examples of images are shown in Fig. 5.

## 2.4. Field data-based optimization of the proposed algorithm

The six datasets presented in section 2.3 were exploited to optimize our registration method. First, we investigated the variability in the optimal band pair combination selected in step 1 (section 2.1.1) to determine if the same combination could be used for all images, which could save some processing time. For every image to be registered, we retained the few best band pair combinations with average NTG values falling within the $[NTG_{min}; NTG_{min} + \Delta NTG]$ interval, as explained in section 2.1.1 (i.e., the best combinations identified before filtering based on average distance between band optical centers). For each dataset, the occurrence of a particular combination in this set of best combinations was then estimated over all images of this dataset, by dividing the number of images for which this combination was one of these best combinations, by the total number of images.

Second, the optimal number of classes $N_c$ used for pixel clustering and pixel filling at step 3 (section 2.1.3) was investigated. For each dataset, the average percentage of filled missing pixels over all images in this dataset was computed for $N_c$ ranging from one to twenty classes. The $N_c$ value was chosen by making a compromise between high percentage of filled missing pixels and low processing runtime, and fixed for comparison with other methods.

## 2.5. Assessment of the proposed algorithm

### 2.5.1. Methods used for comparative performance assessment

The proposed algorithm was compared to six state-of-the-art published registration methods (Barjaktarovic et al., 2024; Dandrifosse et al., 2021; Kim et al., 2017; Sharma et al., 2023; Vayssade et al., 2020; Yasir et al., 2018) using the criteria presented in section 2.5.2. To avoid potential differences in performance due to geometric correction, every band image was undistorted before registration, using intrinsic camera parameters derived from prior stereo calibration (section 2.2).

The first three methods (Sharma et al., 2023; Vayssade et al., 2020; Yasir et al., 2018) rely on the use of affine or homographic global transforms. Sharma et al. (2023) first use depth-dependent affine transformations to roughly align band images. These transformations are derived from prior calibration and their uses require approximate knowledge on the distance between the camera and the imaged scene. Then, the results are refined by optimizing other affine transformations based on Mattes mutual information to minimize the influence of spectral differences between bands. The methods of Vayssade et al. (2020) and Yasir et al. (2018) share similar steps, i.e., (1) feature extraction, (2) identification of band pairs used for registration, and (3) band registration. Like Sharma et al. (2023), Vayssade et al. (2020) add an initial coarse-registration step to the above framework, based on depth-dependent affine transformations. For step (1), both methods test several types of feature points but only the best ones are used in our comparative study, i.e., FAST (FAST Algorithm for Corner Detection) features (Trajković and Hedley, 1998) extracted from normalized gradient band images for Vayssade et al. (2020), and SIFT features extracted from original band images for Yasir et al. (2018). For step (2), both methods identify the optimal combination of band pair alignments based on the number of matched feature points between band images, this combination being of depth one, i.e., $N_b - 1$ bands registered with the remaining master band for Vayssade et al. (2020), or of any depth for Yasir et al. (2018). Step (3) is then the same for both methods, and

consists in finding and applying the $N_b - 1$ homographic transforms between previously identified band pairs to register band images.

The other three methods used for comparison (Barjaktarovic et al., 2024; Dandrifosse et al., 2021; Kim et al., 2017) do not assume that the scene is planar. Of all the methods tested by Dandrifosse et al. (2021), the best one combines a pre-registration step based on homographies, and a refinement step based on 2D spline-based models. Pre-registration is performed using calibration-derived, distance-dependent homographies (Berenstein et al., 2015), similarly to Vayssade et al. (2020) and Sharma et al. (2023). Registration results are then refined by optimizing a local deformation model based on cubic B-splines (Rueckert et al., 1999), using normalized mutual information as similarity measure and 16-pixel spacing between B-splines control points. Barjaktarovic et al. (2024) also combine a coarse registration step and a fine registration step. The former estimates a global translation using phase correlation, while the latter exploits dense optical flow to refine the results. Unlike the above five methods, the algorithm developed by Kim et al. (2017) does not rely on parametric global models. Instead, it first computes the proposed DASC feature for every image pixel, using guided filtering for more efficient processing. DASC features are then densely matched using non-parametric SIFT flow optimization (Liu et al., 2011).

Besides the above six published registration methods, two simpler versions of our method were assessed to see the impact of per-image optimization of band pair combination (step 1, section 2.1.1). These methods use the same band pair combination for every multispectral image. The first version uses the best combination of depth one identified as described in section 2.4, while the second version uses the best combination of depth two.

### 2.5.2. Evaluation criteria

The above registration methods were compared in terms of registration accuracy and processing time (Table 1). In addition, our proposed method was evaluated in terms of crop height estimation accuracy.

Registration performance: three criteria were used to quantify the registration accuracy at different spatial scales. First, for every registered image, the MI global image similarity measure for every band pair, $M_{ij}$, was computed as

$$M_{ij} = \sum_{x=1}^{N_{bins}} \sum_{y=1}^{N_{bins}} h_{ij}(x, y)^* \log \left[ \frac{h_{ij}(x, y)}{\sum_{x=1}^{N_{bins}} h_{ij}(x, y)^* \sum_{y=1}^{N_{bins}} h_{ij}(x, y)} \right] \quad (4)$$

where $h_{ij}(x, y)$ is the normalized joint histogram derived from the $i^{th}$ and $j^{th}$ band images and computed over $N_{bins} = 20$ bins. The average $M_{ij}$ value over all the band pairs, $\overline{M}$, was then computed to have a single value per image. Finally, the average and the standard deviation of $\overline{M}$ over all images in each dataset were computed for each dataset and each registration method. Note that the robust MI measure was selected instead of the NTG measure, to make the evaluation of performance independent from the registration algorithm itself (Dandrifosse et al., 2021; Rohlfing, 2012).

As greater global image similarity does not necessarily mean better registration result (Rohlfing, 2012), two more local criteria were used as a complement. The first one was the average Intersection Over Union (IOU, also called Jaccard Index) computed over three, sufficiently local ROIs for every multispectral image (leading to 120 ROIs in total). For this purpose, the same three ROIs were manually delineated on the six non-registered Airphen spectral bands. In general, these three ROIs, respectively, corresponded to one object on the ground (e.g., a stone), one leaf and one reproductive organ (i.e., a wheat head or a flower). Their sizes strongly differed between objects and datasets, with higher sizes obtained for large maize and sunflower leaves and flowers (5000–10000 pixels in general), and lower sizes obtained for smaller wheat and cover crops leaves and reproductive organs (500–1000 pixels in general) (see Table A1 in supplementary material). For each band,
these three ROIs were used to generate a mask image, on which the same transformations computed for image registration were applied. This allowed us to register the six mask images, corresponding to the six bands. The average IOU over the three ROIs for every band pair, $\overline{IOU}_{i,j}$, was computed as

$$ \overline{IOU}_{i,j} = \frac{1}{3} \sum_{k=1}^{3} \frac{|ROI_{k,i} \cap ROI_{k,j}|}{|ROI_{k,i} \cup ROI_{k,j}|} \qquad (5) $$

where $ROI_{k,i}$ and $ROI_{k,j}$, respectively, refer to the positions of the k th ROI in the i th and j th mask images. Similarly to MI, the average $\overline{IOU}_{i,j}$ value over all the band pairs, $\overline{IOU}$, was computed, and the average and the standard deviation of $\overline{IOU}$ over all images in each dataset were computed for each dataset and each registration method.

The second local criterion used to quantify registration performance was the point-based registration error computed using black crosses drawn on maize leaves in datasets #1 and #2 (section 2.3). For every registered maize image, the position of the center of the k th visible cross, $(x_{i,k}, y_{i,k})$, in the i th band was manually recorded and used as control point. The registration error, $\bar{e}_{i,j}$, was then computed for bands i and j, using all images in datasets #1 and #2. It was defined as the average Euclidean distance between the positions of control points in registered bands i and j (Cao et al., 2020; Chen et al., 2018; Dandrifosse et al., 2021):

$$ \bar{e}_{i,j} = \frac{1}{N_{crosses}} \sum_{k=1}^{N_{crosses}} \sqrt{(x_{i,k} - x_{j,k})^2 + (y_{i,k} - y_{j,k})^2} \qquad (6) $$

where $N_{crosses}$ is the total number of crosses that are visible in the 20 registered maize images. The average $\bar{e}_{i,j}$ value over all the band pairs, $\bar{e}$, was finally computed to provide a single error value for each registration method. In addition, the same errors were computed by filtering control points based on their distances to the camera (defined as the differences between known camera altitude and known cross heights), to observe potential effects of non-negligible relief on registration performance.

**Processing time:** every multispectral image was registered using a 2.20 GHz Intel Core i9-13950HX CPU with 24 physical cores. The average and the standard deviation of processing time was computed over all images for each dataset and each method. Since the codes of the above six state-of-the-art methods (section 2.5.1) were not parallelized, we compared their processing times with those obtained with both non-parallel and parallel (based on Ray Python library) implementations of our algorithm.

**Crop height estimation performance:** Since depth is inversely proportional to disparity (Eq. 2), which is itself used to register the moving band image (Eq. 1), depth estimation accuracy can be considered as a proxy of registration accuracy. The depth image corresponding to the z-component of the 3D point cloud was converted to a height image by (1) manually extracting several sets of soil pixels distributed within the whole scene, (2) fitting a plane to these soil pixels, and (3) subtracting the depth values to this plane. Crop height estimation was first assessed qualitatively, by visual inspection of the height images for every crop. Then, it was quantified using the black crosses of known height that had been drawn on maize leaves (section 2.3): the median height value was computed over a 20-pixel square window around each imaged cross, and compared to the ground-truth height measurement. In total, 304 points were available over all images of both datasets #1 (134 points) and #2 (170 points).

# 3. Results

## 3.1. Identification of best combinations of band pair alignments

The results obtained for the first camera band set (450, 570, 675, 710, 730 and 850 nm) and the corresponding five datasets indicate that

a unique, depth-two band pair combination is nearly always present among the optimal ones: it is selected for 100 % of maize, wheat and sunflower images and for 80 % of cover crop images, for an average occurrence of 96 % (Table 2). According to this combination, the 450 and 675 nm bands should be first registered with the 570 nm band, while the 850 nm band should be registered with the 730 nm band. The 570 and 730 nm bands should be then registered with the 710 nm master band. This best combination shares the latter three band pairs, i.e., (710, 570), (710, 730) and (730, 850), with the second and third best combinations, which are also of depth two (Table 2). The best depth-one combination also uses the 710 nm band as master band, with which every other band is registered (not shown).

Similarly, a unique, depth-two band pair combination is always selected among the optimal ones for the five wheat images acquired with the second camera band set (450, 530, 570, 675, 730 and 850 nm, Table 2). According to this combination, the 675 and 850 nm bands should be first registered with the 450 and 730 nm bands, respectively. The 450, 570 and 730 nm bands should then be registered with the 530 nm master band (which is also the master band selected for the best depth-one band pair combination, result not shown). This best depth-two combination also shares three common band pairs with the second and third best combinations, i.e., (530, 450), (530, 730) and (730, 850) (Table 2).

## 3.2. Optimization of the number of classes for pixel filling

For every dataset, the percentage of missing band values that can be filled in aligned multispectral images increase with the number of classes used in the GMM clustering algorithm (Fig. 6). A large increase is generally observed between one and four classes. After four classes, the evolution of the percentage is variable across bands and datasets: no evolution is observed in some cases (e.g., for maize images and all bands except the blue band), while significant increases are still visible in other cases (e.g., the red band in wheat images). In general, a larger proportion of missing pixels can be filled for maize and sunflower images compared to wheat and cover crop images. For example, it is greater than 50 % for all bands except the blue band in the case of maize images, while it remains lower than 40 % for all bands in the case of wheat images.

Based on these results, we set the number of classes to $N_c = 10$ for the rest of the study. This choice was made to ensure a large number of filled missing pixels in most cases (especially for the blue and red bands), while avoiding an unnecessary increase in processing time, the latter increasing linearly with $N_c$. The example presented in Fig. A2 (supplementary material) illustrates the interest of pixel filling, not only to retrieve objects located behind other objects, but also objects outside the camera field-of-view for one or several band(s).

## 3.3. Registration performance

The registration results expressed in terms of $\overline{Ml}$ global image similarity indicate that the full version of the proposed algorithm (i.e., using an optimal band pair combination for every image) perform similarly or better than other state-of-the-art methods for the six datasets, except Dataset #5 for which the method of Dandrifosse et al. (2021) performs slightly better (Table 3). An average $\overline{Ml}$ of 0.70 is obtained with our algorithm when merging all datasets, followed by Dandrifosse et al. (2021) and Kim et al. (2017) with 0.56, Barjaktarovic et al. (2024) with 0.54, Vayssade et al. (2020) and Yasir et al. (2018) with 0.49 and 0.50, respectively, and Sharma et al. (2023) with 0.40. Using the same combination of band pair alignments for every multispectral image shows slightly decreased or similar performance, with average $\overline{Ml}$ of 0.68 and 0.70 obtained, respectively, when using the best depth-one and depth-two combinations identified in Section 3.1.

Similar observations can be made on $\overline{IOU}$ values calculated over local
<caption>**Table 2**</caption>

Most frequent optimal combinations of band pair alignments selected over the six datasets for the two camera band sets. Each combination in column 2 is represented as in Fig. 1, with each arrow going from the moving image to the fixed image (master band is at the top). From columns 3 to 8, each box corresponding to a dataset D (column) and a band pair combination C (row) indicates the percentage of images of D for which C appears among the best combinations selected (section 2.4). For each camera band set, the band pairs shared by all best combinations are in orange, and the combinations are sorted from the most frequent to the least frequent when considering all datasets (last column), keeping only the best three combinations.

<table><thead><tr><th>Band set (in nm)</th><th>Best combinations of band pair alignments</th><th>Dataset #1<br>Maize<br>Tripod<br>0.5-1.5 m</th><th>Dataset #2<br>Maize<br>Tripod<br>0.9-1.9 m</th><th>Dataset #3<br>Wheat<br>Tripod<br>1.4-2.2 m</th><th>Dataset #4<br>Wheat<br>Phenom.<br>1.3-2.4 m</th><th>Dataset #5<br>Sunflower<br>Phenom.<br>1.2-3.0 m</th><th>Dataset #6<br>Cover Crops<br>Phenom.<br>1.3-2.0 m</th><th>Average</th></tr></thead><tbody><tr><td rowspan="3">450, 570, 675, 710, 730, 850</td><td>Image 1</td><td>100</td><td>100</td><td>–</td><td>100</td><td>100</td><td>80</td><td>96</td></tr><tr><td>Image 2</td><td>82</td><td>33</td><td>–</td><td>100</td><td>0</td><td>0</td><td>43</td></tr><tr><td>Image 3</td><td>0</td><td>0</td><td>–</td><td>0</td><td>100</td><td>0</td><td>20</td></tr><tr><td rowspan="3">450, 530, 570, 675, 730, 850</td><td>Image 1</td><td>–</td><td>–</td><td>100</td><td>–</td><td>–</td><td>–</td><td>100</td></tr><tr><td>Image 2</td><td>–</td><td>–</td><td>80</td><td>–</td><td>–</td><td>–</td><td>80</td></tr><tr><td>Image 3</td><td>–</td><td>–</td><td>80</td><td>–</td><td>–</td><td>–</td><td>80</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Percentages of filled missing pixels as functions of the number of classes used in the GMM, for each of the six datasets and for each band to be filled. Markers and error bars, respectively, show the averages and the standard deviations of percentage values obtained over all images of the dataset. Note that 710 nm and 530 nm bands are not shown, since these two bands are always selected as master bands for camera band sets 1 and 2, respectively (section 3.1), so they do not have missing pixels.</figcaption></figure>
<caption>**Table 3**</caption>

Averages and standard deviations of $\overline{\text{MI}}$ obtained for every dataset and every tested method. $\overline{\text{MI}}$ values obtained by merging all datasets are shown in the ‘All’ column. For each column, the best average $\overline{\text{MI}}$ value(s) is(are) in bold.

<table><thead><tr><th>Methods</th><th>Dataset #1<br>Maize<br>Tripod<br>0.5–1.5 m</th><th>Dataset #2<br>Maize<br>Tripod<br>0.9–1.9 m</th><th>Dataset #3<br>Wheat<br>Tripod<br>1.4–2.2 m</th><th>Dataset #4<br>Wheat<br>Phenom.<br>1.3–2.4 m</th><th>Dataset #5<br>Sunflower<br>Phenom.<br>1.2–3.0 m</th><th>Dataset #6<br>Cover Crops<br>Phenom.<br>1.3–2.0 m</th><th>All</th></tr></thead><tbody><tr><td>Sharma et al (2023)</td><td>0.40 ± 0.05</td><td>0.35 ± 0.05</td><td>0.40 ± 0.02</td><td>0.29 ± 0.05</td><td>0.59 ± 0.07</td><td>0.38 ± 0.05</td><td>0.40 ± 0.10</td></tr><tr><td>Vayssade et al (2020)</td><td>0.57 ± 0.10</td><td>0.61 ± 0.10</td><td>0.41 ± 0.04</td><td>0.29 ± 0.04</td><td>0.50 ± 0.04</td><td>0.39 ± 0.08</td><td>0.49 ± 0.14</td></tr><tr><td>Yasir et al (2018)</td><td>0.56 ± 0.10</td><td>0.62 ± 0.10</td><td>0.42 ± 0.04</td><td>0.28 ± 0.04</td><td>0.55 ± 0.11</td><td>0.42 ± 0.09</td><td>0.50 ± 0.14</td></tr><tr><td>Dandrifosse et al (2021)</td><td>0.52 ± 0.07</td><td>0.56 ± 0.08</td><td>0.54 ± 0.05</td><td>0.47 ± 0.04</td><td><strong>0.78 ± 0.06</strong></td><td>0.54 ± 0.10</td><td>0.56 ± 0.11</td></tr><tr><td>Barjaktarovic et al (2024)</td><td>0.52 ± 0.13</td><td>0.56 ± 0.10</td><td>0.53 ± 0.04</td><td>0.44 ± 0.04</td><td>0.70 ± 0.07</td><td>0.50 ± 0.03</td><td>0.54 ± 0.11</td></tr><tr><td>Kim et al (2017)</td><td>0.60 ± 0.07</td><td>0.61 ± 0.07</td><td>0.44 ± 0.06</td><td>0.40 ± 0.03</td><td>0.71 ± 0.06</td><td>0.48 ± 0.12</td><td>0.56 ± 0.12</td></tr><tr><td>This study (same depth-one band pair combination for every image)</td><td>0.75 ± 0.07</td><td>0.74 ± 0.07</td><td>0.56 ± 0.05</td><td>0.53 ± 0.05</td><td><strong>0.78 ± 0.07</strong></td><td><strong>0.59 ± 0.09</strong></td><td>0.68 ± 0.12</td></tr><tr><td>This study (same depth-two band pair combination for every image)</td><td>0.78 ± 0.07</td><td><strong>0.77 ± 0.08</strong></td><td><strong>0.60 ± 0.06</strong></td><td><strong>0.55 ± 0.05</strong></td><td>0.76 ± 0.07</td><td>0.58 ± 0.09</td><td><strong>0.70 ± 0.12</strong></td></tr><tr><td>This study (optimal band pair combination for every image)</td><td><strong>0.80 ± 0.08</strong></td><td>0.76 ± 0.08</td><td>0.59 ± 0.06</td><td><strong>0.55 ± 0.05</strong></td><td>0.77 ± 0.07</td><td><strong>0.59 ± 0.09</strong></td><td><strong>0.70 ± 0.12</strong></td></tr></tbody></table>

ROIs (Table 4). The full version of the proposed algorithm leads to mean $\overline{\text{IOU}}$ values ranging from 0.78 to 0.87 for the six datasets, with an average $\overline{\text{IOU}}$ of 0.82 when merging all datasets. It generally performs better than the method of Barjaktarovic et al. (2024) that leads to an average $\overline{\text{IOU}}$ of 0.67, followed by the methods of Kim et al. (2017) with an average $\overline{\text{IOU}}$ of 0.65, Dandrifosse et al. (2021) and Sharma et al. (2023) with average $\overline{\text{IOU}}$ of around 0.57, and Vayssade et al. (2020) and Yasir et al. (2018) with average $\overline{\text{IOU}}$ of around 0.49. It can be noted that the proposed algorithm also provides significantly less variable results, with a standard deviation three to four times lower than the other algorithms. Finally, Table 4 shows that using the same depth-two, combination of band pair alignments for every multispectral image provides similar results as compared to using an optimal combination per image (average $\overline{\text{IOU}}$ of 0.82 when merging all datasets for both methods). On the other hand, using the same depth-one combination slightly decreases the overall performance (average $\overline{\text{IOU}}$ of 0.80).

The different algorithms tested generally rank similarly when expressing registration performance in terms of point-based registration error computed using the black crosses drawn on maize leaves in datasets #1 and #2 (Fig. 7). Without filtering control points based on their distances to the camera, the full version of the proposed method obtains the lowest error ($\bar{\epsilon} = 0.5$ pixel), closely followed by the other two simplified versions ($\bar{\epsilon} = 0.6$ pixel and $\bar{\epsilon} = 0.7$ pixel when using the same depth-two and depth-one combinations of band pair alignments, respectively) and the method of Kim et al. (2017) ($\bar{\epsilon} = 1.0$ pixel). The algorithms proposed by Dandrifosse et al. (2021), Barjaktarovic et al. (2024), Sharma et al. (2023), Yasir et al. (2018) and Vayssade et al. (2020) perform significantly worse, with $\bar{\epsilon} \ge 5.4$ pixels and strongly

variable performance over the range of distances investigated: the error increases from around 2.5 pixels for distances close to 1.60 m, to more than 15 pixels for distances close to 0.80 m (Fig. 7). On the other hand, the performances of other methods are more stable, although the method of Kim et al. (2017) shows an increase in the error for the lowest distance only.

The registration results quantified in Tables 3 and 4 and Fig. 7 can be qualitatively assessed in Fig. 8, using some examples of registered true color composite images that were extracted from five out of the six images presented in Fig. 5. The superior performance of the proposed algorithm can be observed, especially for pixels corresponding to the top of canopy, e.g., for maize and sunflower flowers. Note the presence of white, missing pixels in the images that were registered with our algorithm. These remaining pixels could not be filled with sufficient accuracy for at least one of the red, green and blue bands used for image composition, and correspond to occluded areas as identified by the stereovision algorithm.

### 3.4. Processing time performance

The method proposed by Vayssade et al. (2020), Barjaktarovic et al. (2024), Yasir et al. (2018) and Sharma et al. (2023) were significantly faster than the other ones, with average processing times of 2, 4, 5, and 6 s, respectively, per multispectral image, when merging all datasets (Table 5). The three non-parallelized versions of the proposed algorithm took around 101 s per image, the two methods without optimization of band pair combination being slightly faster than the full version. They were faster than the methods of Dandrifosse et al. (2021) and Kim et al.

<caption>**Table 4**</caption>

Averages and standard deviations of $\overline{\text{IOU}}$ obtained for every dataset and every tested method. $\overline{\text{IOU}}$ values obtained by merging all datasets are shown in the ‘All’ column. For each column, the best average $\overline{\text{IOU}}$ value(s) is(are) in bold.

<table><thead><tr><th>Methods</th><th>Dataset #1<br>Maize<br>Tripod<br>0.5–1.5 m</th><th>Dataset #2<br>Maize<br>Tripod<br>0.9–1.9 m</th><th>Dataset #3<br>Wheat<br>Tripod<br>1.4–2.2 m</th><th>Dataset #4<br>Wheat<br>Phenom.<br>1.3–2.4 m</th><th>Dataset #5<br>Sunflower<br>Phenom.<br>1.2–3.0 m</th><th>Dataset #6<br>Cover Crops<br>Phenom.<br>1.3–2.0 m</th><th>All</th></tr></thead><tbody><tr><td>Sharma et al (2023)</td><td>0.56 ± 0.12</td><td>0.58 ± 0.12</td><td>0.39 ± 0.05</td><td>0.57 ± 0.11</td><td>0.55 ± 0.21</td><td>0.58 ± 0.11</td><td>0.55 ± 0.14</td></tr><tr><td>Vayssade et al (2020)</td><td>0.57 ± 0.10</td><td>0.58 ± 0.20</td><td>0.31 ± 0.07</td><td>0.37 ± 0.17</td><td>0.50 ± 0.13</td><td>0.38 ± 0.19</td><td>0.48 ± 0.18</td></tr><tr><td>Yasir et al (2018)</td><td>0.56 ± 0.10</td><td>0.59 ± 0.19</td><td>0.35 ± 0.06</td><td>0.43 ± 0.14</td><td>0.52 ± 0.17</td><td>0.43 ± 0.16</td><td>0.50 ± 0.17</td></tr><tr><td>Dandrifosse et al (2021)</td><td>0.65 ± 0.10</td><td>0.62 ± 0.20</td><td>0.43 ± 0.06</td><td>0.54 ± 0.21</td><td>0.68 ± 0.17</td><td>0.42 ± 0.23</td><td>0.58 ± 0.19</td></tr><tr><td>Barjaktarovic et al (2024)</td><td>0.56 ± 0.16</td><td>0.75 ± 0.10</td><td>0.46 ± 0.04</td><td>0.78 ± 0.06</td><td>0.67 ± 0.24</td><td><strong>0.83 ± 0.06</strong></td><td>0.67 ± 0.18</td></tr><tr><td>Kim et al (2017)</td><td>0.66 ± 0.11</td><td>0.69 ± 0.09</td><td>0.37 ± 0.09</td><td>0.76 ± 0.07</td><td>0.63 ± 0.18</td><td>0.75 ± 0.15</td><td>0.65 ± 0.16</td></tr><tr><td>This study (same depth-one band pair combination for every image)</td><td>0.78 ± 0.05</td><td>0.83 ± 0.04</td><td>0.77 ± 0.03</td><td>0.78 ± 0.05</td><td><strong>0.88 ± 0.02</strong></td><td>0.80 ± 0.06</td><td>0.80 ± 0.06</td></tr><tr><td>This study (same depth-two band pair combination for every image)</td><td><strong>0.81 ± 0.04</strong></td><td><strong>0.84 ± 0.04</strong></td><td><strong>0.79 ± 0.03</strong></td><td><strong>0.80 ± 0.04</strong></td><td><strong>0.88 ± 0.02</strong></td><td>0.79 ± 0.05</td><td><strong>0.82 ± 0.05</strong></td></tr><tr><td>This study (optimal band pair combination for every image)</td><td><strong>0.81 ± 0.05</strong></td><td><strong>0.84 ± 0.04</strong></td><td>0.78 ± 0.03</td><td><strong>0.80 ± 0.04</strong></td><td>0.87 ± 0.02</td><td>0.80 ± 0.05</td><td><strong>0.82 ± 0.05</strong></td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>Fig. 7. Registration error computed from black crosses drawn on maize leaves in datasets #1 and #2, as a function of the distance between the camera and the black crosses. All the tested methods are presented in (a), and a zoom on the most accurate methods is shown in (b). A distance sampling interval of 20 cm is used to compute the error. The total error value, $\bar{\epsilon}$, is shown in the legend.</figcaption></figure>

(2017), that had average processing times of 185 and 648 s, respectively. The processing time required by the methods of Dandrifosse et al. (2021), Kim et al. (2017) and Sharma et al. (2023) was relatively stable across images, with relative standard deviations lower than 17 % (Table 5). It was more variable for the proposed algorithms, with relative standard deviations of around 27%; the processing time was higher for maize and sunflower images, and lower for wheat images. The methods of Vayssade et al. (2020) and, to a lower extent, Yasir et al. (2018) show the most variable processing times relatively to the averages, with relative standard deviations of 100 % and 60 %.

Parallel implementation of the proposed method allowed us to reduce the processing time by about 60 % overall (Table 5). On average, the parallelized versions of the proposed algorithm required 38 s and 40 s per image when using the same depth-one and depth-two band pair combinations, respectively, while the full version required 42 s per image.

### 3.5. Height estimation

The height images derived from image registration provide detailed information on the 3D structure of the imaged scene (Fig. 9). For example, large and planophile individual sunflower leaves or smaller wheat heads can be well identified. Validation using black crosses of known height on maize datasets #1 and #2 demonstrate the centimetric accuracy of such height images (Fig. 10). The RMSE is 1.9 cm when computed over the two datasets, and even reaches 1.2 cm when focusing on dataset #1 acquired from shorter distance.

## 4. Discussion

### 4.1. The proposed method offers a good compromise between accuracy and efficiency for close-range multispectral image registration

The results obtained with the three global ($\bar{Ml}$) and local ($\bar{IOU}$ and $\bar{\epsilon}$) registration measures (Tables 3 and 4, Fig. 7) demonstrate that our stereovision-based method provides the highest registration accuracy,

followed by the optical flow methods of Kim et al. (2017) and Barjaktarovic et al. (2024), the spline-based method of Dandrifosse et al. (2021), and the three methods of Sharma et al. (2023), Vayssade et al. (2020), and Yasir et al. (2018) based on affine or homographic transformations. The poorer performances obtained with the latter three methods seemingly come from the assumption of planar scene that is made when modeling changes in point of view using affine or homographic transformations. This assumption is here violated since the heights of the plants and the distances between the camera and the plants are of the same order of magnitude. This is particularly visible in Fig. 7, in which these methods show a strong increase in registration error when the distance between the camera and the plants gets closer to the depth of the imaged scene.

To correct for residual registration errors that remain after applying global transforms, Barjaktarovic et al. (2024) and Dandrifosse et al. (2021), respectively, used a dense optical flow algorithm and a grid of 2D spline basis functions, to warp locally the moving images so that they match the fixed image. Tables 3 and 4, and Fig. 7 show that this generally improves the registration accuracy, compared to affine- and homography-based methods. However, a strong decrease in performance is still observed when the distance between the camera and the plants decreases (Fig. 7), indicating that such residuals may be too large to be properly corrected with the 2D spline basis functions and dense optical flow algorithm used.

On the other hand, our stereovision-based method and the method of Kim et al. (2017) provide accurate registration results for all the distances investigated (Fig. 7). Overall, the best results obtained with our method based on the three $\bar{Ml}$, $\bar{IOU}$ and $\bar{\epsilon}$ measures may be due to the fact that we exploit the rigidity of the multi-lens camera and the synchronized capture of multispectral bands to transform the 2D matching problem into a better-constrained 1D one, thanks to prior stereo calibration. The added value of exploiting such constraints is particularly visible on the IOU results obtained on 120 ROIs, which are significantly better and less variable over all tested datasets and crops, even for smallest and closest objects. However, it should be noted that, unlike the method of Kim et al. (2017), our method cannot be applied without a
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Examples of true color composite images obtained after image registration with every tested method for maize (column 1), wheat (columns 2 and 3), sunflower (column 4) and faba bean cover crops (column 5). For greater visiblity, these images are 300 x 300 pixel sub-images extracted in zones with high depth variability. The images obtained in the row « This study » are obtained with the full version of the proposed algorithm, i.e., using an optimal combination of band pair alignments for every image. In the last row, white pixels are the remaining, unfilled missing pixels for at least one of the red, green and blue bands used for image composition. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<caption>**Table 5**</caption>

Averages and standard deviations of processing time (in seconds) obtained for every dataset and every tested method. Processing times obtained by merging all datasets are shown in the ‘All’ column. For each column, the best average processing time is in bold.

<table><thead><tr><th>Methods</th><th></th><th>Dataset #1<br>Maize<br>Tripod<br>0.5–1.5 m</th><th>Dataset #2<br>Maize<br>Tripod<br>0.9–1.9 m</th><th>Dataset #3<br>Wheat<br>Tripod<br>1.4–2.2 m</th><th>Dataset #4<br>Wheat<br>Phenom.<br>1.3–2.4 m</th><th>Dataset #5<br>Sunflower<br>Phenom.<br>1.2–3.0 m</th><th>Dataset #6<br>Cover<br>Crops<br>Phenom.<br>1.3–2.0 m</th><th>All</th></tr></thead><tbody><tr><td rowspan="5">Non-parallelized</td><td>Sharma et al (2023)</td><td>5 ± 0</td><td>6 ± 0</td><td>6 ± 1</td><td>6 ± 0</td><td>7 ± 0</td><td>6 ± 1</td><td>6 ± 1</td></tr><tr><td>Vayssade et al (2020)</td><td><b>1 ± 0</b></td><td><b>1 ± 0</b></td><td>6 ± 1</td><td>4 ± 0</td><td><b>1 ± 0</b></td><td><b>3 ± 1</b></td><td><b>2 ± 2</b></td></tr><tr><td>Yasir et al (2018)</td><td>2 ± 0</td><td>5 ± 1</td><td>7 ± 0</td><td>6 ± 0</td><td>3 ± 0</td><td>9 ± 1</td><td>5 ± 3</td></tr><tr><td>Dandrifosse et al (2021)</td><td>198 ± 14</td><td>180 ± 1</td><td>180 ± 2</td><td>179 ± 0</td><td>184 ± 6</td><td>179 ± 0</td><td>185 ± 11</td></tr><tr><td>Barjaktarovic et al (2024)</td><td>3 ± 0</td><td>4 ± 1</td><td><b>4 ± 1</b></td><td><b>3 ± 0</b></td><td>6 ± 1</td><td>4 ± 1</td><td>4 ± 1</td></tr><tr><td rowspan="5">Parallelized</td><td>Kim et al (2017)</td><td>663 ± 58</td><td>652 ± 22</td><td>636 ± 2</td><td>639 ± 4</td><td>637 ± 4</td><td>637 ± 2</td><td>648 ± 34</td></tr><tr><td>This study (same depth-one band pair combination for every image)</td><td>109 ± 25</td><td>91 ± 12</td><td>62 ± 1</td><td>76 ± 2</td><td>104 ± 14</td><td>88 ± 13</td><td>92 ± 22</td></tr><tr><td>This study (same depth-two band pair combination for every image)</td><td>118 ± 28</td><td>98 ± 13</td><td>66 ± 2</td><td>82 ± 3</td><td>110 ± 15</td><td>94 ± 14</td><td>98 ± 25</td></tr><tr><td>This study (optimal band pair combination for every image)</td><td>134 ± 32</td><td>122 ± 45</td><td>69 ± 2</td><td>90 ± 3</td><td>119 ± 12</td><td>102 ± 13</td><td>112 ± 35</td></tr><tr><td>This study (same depth-one band pair combination for every image)</td><td>46 ± 11</td><td>38 ± 6</td><td>26 ± 1</td><td>31 ± 1</td><td>41 ± 9</td><td>35 ± 5</td><td>38 ± 10</td></tr><tr><td rowspan="2">This study (same depth-two band pair combination for every image)</td><td>50 ± 11</td><td>38 ± 8</td><td>28 ± 1</td><td>34 ± 2</td><td>40 ± 7</td><td>40 ± 6</td><td>40 ± 11</td></tr><tr><td>This study (optimal band pair combination for every image)</td><td>53 ± 11</td><td>40 ± 8</td><td>27 ± 1</td><td>36 ± 2</td><td>47 ± 6</td><td>40 ± 6</td><td>42 ± 11</td></tr></tbody></table>

good synchronization of image capture, a camera calibration preliminary step and a rigid stereo rig that ensures that the estimated extrinsic parameters do not change over time. However, these constraints can be generally met with most multi-lens cameras.

The proposed method estimates a disparity image per band pair so it is, of course, less efficient than the methods of Sharma et al. (2023), Vayssade et al (2020) and Yasir et al (2018) based on affine or homographic transformations (Table 5), for which only a few parameters per band pair have to be retrieved based on extracted sparse feature points. The non-parallelized version of the proposed method is, however, more efficient than the methods of Kim et al (2017) and Dandrifosse et al (2021) (Table 5), probably because it solves a simpler, computationally less expensive 1D matching problem, compared to the 2D problem considered by these other two methods. Unlike these latter two methods, the processing time required by our method is image-dependent: it is proportional to the range of disparity values investigated, which is itself related not only to the average distance between the camera and the imaged scene but also to the depth of the imaged scene (Eq. 2). This explains why the processing time is higher for tall crops seen at short distance, i.e., maize and sunflower, compared to smaller crops seen at a greater distance, i.e., wheat and cover crops (Table 5). Note that, like other methods, our method can be easily parallelized for further gain in processing time, this gain being maximum if the number of physical cores is equal or greater than the number of band pairs to be aligned.

## 4.2. The proposed method properly handles occluded areas using pixel masking and filling

When the average distance between the camera and the imaged scene decreases and/or the depth of the imaged scene increases, the performance of close-range multispectral image registration may be severely hindered by the presence of objects that are not seen from all band optical centers, i.e., objects outside the camera field-of-view for at least one band, and objects located behind other objects in the scene (occlusions). Unlike the other methods tested in this paper, the proposed method properly handles these two kinds of objects using pixel masking and pixel filling.

First, these pixels are detected and masked thanks to the consistency check that is performed by the stereovision algorithm and that is based on a comparison between the left-right and right-left disparity images

(Hirschmüller, 2008). Note that a few non-occluded pixels may sometimes be wrongly detected as occluded, when disparity cannot be accurately estimated (e.g., see the leaves on the 850 nm band in Fig. A2).

Then, unlike stereovision algorithms that fill holes in disparity maps to retrieve depth in occlusions (Hirschmüller, 2008; Zbontar and LeCun, 2016), the proposed method fills holes in the aligned multispectral bands to retrieve band values in occlusions. It exploits the high spectral covariance observed within a class of material, which allows a missing band value to be retrieved based on other band value(s) and a simple MLR model. Pixel filling is an important step of our method because the number of missing pixels can be relatively high, especially for bands at depth two that are aligned with master band using two successive disparity maps (see how missing pixels adds up when registering the 850 nm band with the 710 nm master band through the 730 nm one in Fig. A2). Clustering in groups of spectrally similar pixels (e.g., soil versus green vegetation versus yellow vegetation, shaded versus sunlit, etc) before applying MLR models is effective to improve pixel filling, as shown by the increasing performance with the number of classes (Fig. 6). That is due to the resulting increasing correlation between bands that facilitates the prediction of one unknown band value based on other available band values. In general, good filling performances are already obtained for the sparsest canopies tested in this paper (sunflower and maize) using a GMM with a few classes, roughly corresponding to green and yellow vegetation, soil and/or reproductive organs. Adding more classes allows us to improve regression model performance and pixel filling for denser canopies with more complex illumination variations within the scene, e.g., wheat canopies.

## 4.3. The proposed method can be simplified by using a unique, image-independent combination of band pair alignments for a given type of scene

When imaging the same type of scene (vegetation and soil background in our study), our findings suggest that our method can be simplified by using a unique, image-independent combination of band pair alignments, without significant decrease in registration accuracy. In this paper, this optimal combination has been identified over 35 images of various crops for the first camera band set, and over 5 images of wheat crops for the second camera band set. Given the already large diversity in canopy structure observed over the 35 images (more diversity is still probably required for the second camera band set), the optimal depth-
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** True color composite images (column 1) and estimated height images (column 2) obtained with the proposed registration algorithm for maize (row 1), wheat (rows 2 and 3), sunflower (row 4) and faba bean cover crops (row 5). The white pixels are the remaining, unfilled missing pixels.</figcaption></figure>

one and depth-two band pair combinations identified in section 3.1 for the first camera band set can be used directly to register new images of crops acquired with this camera, thus skipping the first step of the proposed method (section 2.1.1). This can decrease the processing time by up to 18 %, with a slightly more important gain when using a unique depth-one combination instead of a unique depth-two combination (Table 5).

The best depth-one and depth-two combinations identified in section 3.1 are consistent with the literature focusing on vegetation scenes (Dias Junior et al., 2020; Fernández et al., 2021; Meng et al., 2021; Vayssade et al., 2020; Yasir et al., 2018), and reflect similarities between spectral

bands for soil and vegetation. For example, green or red-edge bands are always selected as master bands due to the intermediate response of soil and vegetation at these wavelengths compared to blue, red and NIR wavelengths (corresponding to either strong or low chlorophyll absorption). Optimal combinations also prevent two bands with strongly different spectral responses to be selected as a band pair, e.g., red and NIR bands in the case of vegetation. Note that these optimal band combinations are expected to change when considering other types of scenes with different optical properties, such as historic monuments (Adamopoulos and Rinaudo, 2021b), wall paintings (Grifoni et al., 2022), sedimentary rocks (Pozo et al., 2015) or aquatic habitats (Cortesi
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Height estimation results obtained on the 40 black crosses drawn on maize leaves (datasets #1 and #2). The total RMSE computed over both datasets is 1.9 cm.</figcaption></figure>

et al., 2021).

The choice of a depth-one or a depth-two combination depends on the scene being imaged, and results from a compromise between registration accuracy, processing time and number of occlusions. Compared to using a depth-one combination, using a depth-two combination leads to slightly higher registration accuracy thanks to the use of intermediate bands (Dias Junior et al., 2020; Vayssade et al., 2020; Yasir et al., 2018), however, at the price of higher processing time (Table 5) and larger number of occlusions before pixel filling (e.g., compare the number of occluded pixels in 730 nm and 850 nm bands in Fig. A2).

### 4.4. The proposed method provides a multispectral point cloud of the scene in a single image capture

The proposed method outputs $N_b$ multispectral bands and a 3D point cloud in metric unit, which are all aligned with each other. Even if the $N_b$ viewing directions (corresponding to the $N_b$ band optical centers) are not so different, the 3D point cloud is of good quality, as demonstrated by the detailed height images in Fig. 9 and the centimetric height estimation accuracy in Fig. 10 (note that manual height measurements are also prone to errors (Madec et al., 2017)). This good quality is not only due to the use of robust criteria against radiometric differences (census transform combined with Hamming distance for the stereo matching cost, NTG measure to select the best band pairs), but also to the fusion of the $N_b$ – 1 disparity images corresponding to different points of view (section 2.1.2).

The proposed method thus provides a multispectral point cloud of the imaged scene, i.e., three spatial dimensions and one spectral dimension. This 4D description is here easily obtained through a single multispectral image capture, e.g., compared to other multimodal solutions based on multi/hyperspectral imagery and LiDAR (Lerma et al., 2012; Rincón et al., 2022; Tusa et al., 2019) or RGB-D data (Sun et al., 2019; Vuletić et al., 2023, 2022; Xie et al., 2023), that are more expensive and complex to implement. This opens up interesting prospects for close-range sensing, in particular for vegetation characterization since both biochemical (leaf pigments) and structural (height, volume, leaf area index...) vegetation traits can be estimated with a single multispectral image. In addition, knowledge of camera footprint

in metric unit may help normalize some of these traits with respect to soil area (e.g., plant count to plant density). Finally, knowing the 3D position of each pixel makes it possible to focus the spectral analysis on objects located within a specific depth range. For example, all vegetation pixels (i.e., both green and senescent) can be easily discriminated from background pixels, a task that might be challenging when exploiting only spectral information due to the similar signatures of soil and senescent vegetation pixels (Fig. 11.a). The height image can also be useful to estimate the vertical distribution of leaf biochemical traits, or simply to compute the spectra of plant organs such as leaves and wheat ears (Fig. 11.b-d), which is critical to improve the spectral sensitivity to their biochemistry (Jay et al., 2019, 2017; Moorthy et al., 2008; Zarco-Tejada et al., 2004).

## 5. Conclusions and perspectives

In this study, we propose a method to register a multispectral image acquired with a multi-lens camera over a non-planar scene, i.e., whose depth is comparable to the distance between the camera and the scene. The method exploits the rigidity of these cameras and the synchronized capture of multispectral bands, both of which are therefore compulsory constraints. For each band pair to be aligned, it transforms the 2D matching problem into a simpler stereo matching problem that is composed of a camera extrinsic parameter estimation and a per-pixel disparity estimation. Images of four contrasted crops acquired with two camera band sets were used to compare the performance of the proposed method with those of six state-of-the-art methods. The results demonstrate that our method leads to the highest registration accuracy, while keeping reasonable processing time. It also provides an accurate 3D point cloud aligned with the multispectral bands, thus constituting a multispectral point cloud of the imaged scene, which is classically obtained using more acquisition complex systems.

Immediate prospects include using this algorithm to register close-range multispectral images for biochemistry estimation in vegetation studies, possibly by exploiting the 3D point cloud as explained in section 4.4. Since the proposed method has been tested only on vegetation images in this study, it will be interesting to assess its performance on other types of images, e.g., to assess deterioration of historic monuments and structures (Adamopoulos and Rinaudo, 2021b), to digitize wall paintings (Grifoni et al., 2022), to discriminate sedimentary rocks (Pozo et al., 2015) or to detect plastic pollution in aquatic habitats (Cortesi et al., 2021). Also, the code implementation should be further optimized to improve the processing time, which might become a limiting factor when applying the method to multispectral images with higher resolution compared to the 1.2 MP Airphen images used in this study. Alternatively, the method could be used to generate a set of pseudo ground-truth disparity images, which could be then used to train a deep learning model that could be much faster to implement once it has been trained (Tosi et al., 2022).

For those interested, the codes and the datasets used in this study are publicly available at https://forgemia.inra.fr/sylvain.jay/close-range_multispectral-registration.git.

### CRediT authorship contribution statement

**Sylvain Jay:** Writing – review & editing, Writing – original draft, Visualization, Validation, Methodology, Investigation, Formal analysis, Conceptualization. **Frédéric Baret:** Writing – review & editing, Supervision, Funding acquisition, Conceptualization. **Samuel Thomas:** Methodology. **Marie Weiss:** Writing – review & editing, Supervision, Methodology, Funding acquisition.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence
<figure><img src="image_14.png" /><figcaption>Fig. 11. Extraction of background, ear and leaf spectra based on estimated height image for a wheat scene imaged under natural illumination (row 2 in Fig. 9): (a) Pixels with height lower than 10 cm are assigned to the background class; (b) Pixels with height higher than 84 % of canopy height (d'Argaignon et al., 2024) are assigned to the ear class; (c) Other pixels are assigned to the leaf class; (d) Spectra can then be averaged over only the pixels of interest.</figcaption></figure>

the work reported in this paper.

## Acknowledgments

This work was supported by the French National Research Agency in the framework of the "Investissements d'avenir" program PHENOME (ANR-11-INBS-012). We would like to thank Sebastien Dandrifosse, Seungryong Kim, Neelesh Sharma and Jehan-Antoine Vayssade for sharing their codes, as well as Remy Marandel and Mathieu Roy for their help in the experiments.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.isprsjprs.2025.06.001.

## References

Adamopoulos, E., Rinaudo, F., 2021a. Combining multiband imaging, photogrammetric techniques, and FOSS GIS for affordable degradation mapping of stone monuments. Buildings 11, 304. https://doi.org/10.3390/buildings11070304.

Adamopoulos, E., Rinaudo, F., 2021b. Close-range sensing and data fusion for built heritage inspection and monitoring—A review. Remote Sens. 13, 3936. https://doi.org/10.3390/rs13193936.

Banerjee, B.P., Raval, S.A., Cullen, P.J., 2018. Alignment of UAV-hyperspectral bands using keypoint descriptors in a spectrally complex environment. Remote Sens. Lett. 9, 524–533. https://doi.org/10.1080/2150704X.2018.1446564.

Barjaktarovic, M., Santoni, M., Bruzzone, L., 2024. Design and verification of a low-cost multispectral Camera for Precision agriculture application. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 17, 6945–6957. https://doi.org/10.1109/JSTARS.2024.3377104.

Beaupre, D.-A., Bilodeau, G.-A., 2021. Domain siamese CNNs for sparse multispectral disparity estimation. In: 2020 25th International Conference on Pattern Recognition (ICPR). Presented at the 2020 25th International Conference on Pattern Recognition (ICPR), pp. 3667–3674, 10.1109/ICPR48806.2021.9412723.

Berenstein, R., Hocevar, M., Godesa, T., Edan, Y., Ben-Shahar, O., 2015. Distance-dependent multimodal image registration for agriculture tasks. Sensors 15, 20845–20862. https://doi.org/10.3390/s150820845.

Berger, K., Machwitz, M., Kycko, M., Kefauver, S.C., Van Wittenberghe, S., Gerhards, M., Verrelst, J., Atzberger, C., van der Tol, C., Damm, A., Rascher, U., Herrmann, I., Paz, V.S., Fahrner, S., Pieruschka, R., Prikaziuk, E., Buchaillot, Ma.L., Halabuk, A., Celesti, M., Koren, G., Gormus, E.T., Rossini, M., Foerster, M., Siegmann, B., Abdelbaki, A., Tagliabue, G., Hank, T., Darvishzadeh, R., Asen, H., Garcia, M., Poças, I., Bandopadhyay, S., Sulis, M., Tomelleri, E., Rozenstein, O., Filchev, L., Stancile, G., Schlerf, M., 2022. Multi-sensor spectral synergies for crop stress detection and monitoring in the optical domain: a review. Remote Sens. Environ. 280, 113198. https://doi.org/10.1016/j.rse.2022.113198.

Bradski, G., 2000. The openCV library. Dr Dobbs. J. Softw. Tools Prof. Program. 25, 120–123.

Cao, S.Y., Shen, H.L., Chen, S.J., Li, C., 2020. Boosting structure consistency for multispectral and multimodal image registration. IEEE Trans. Image Process. 29, 5147–5162. https://doi.org/10.1109/TIP.2020.2980972.

Chang, Y.-J., Lee, B.-G., Jeon, M., 2022. Warping-based spectral translation network for unsupervised cross-spectral stereo matching. Inf. Sci. 588, 214–230. https://doi.org/10.1016/j.ins.2021.12.075.

Chen, S.-J., Shen, H.-L., Li, C., Xin, J.H., 2018. Normalized Total gradient: a new measure for multispectral image registration. IEEE Trans. Image Process. 27, 1297–1310. https://doi.org/10.1109/TIP.2017.2776753.

Cortesi, I., Masiero, A., De Giglio, M., Tucci, G., Dubbini, M., 2021. Random forest-based river plastic detection with a handheld multispectral camera. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 9–14. https://doi.org/10.5194/isprsjournals-XLIII-B1-2021-9-2021.

Cournet, M., Sarrazin, E., Dumas, L., Michel, J., Guinet, J., Youssefi, D., Defonte, V., Fardet, Q., 2020. Ground truth generation and disparity estimation for optical satellite imagery. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. - ISPRS Arch. 43, 127–134. https://doi.org/10.5194/isprsjournals-XLIII-B2-2020-127-2020.

d'Argaignon, M.-P., Lopez Lozano, R., Jay, S., Ausset, A., Berthon, B., Chapuis, R., De Solan, B., Grau, A., Larue, F., Le Roy, R., Marandel, R., Mercier, V., Roy, M., Tison, G., Venault, F., Martre, P., 2024. Estimation of the average wheat ear size and ear density at the microplot level from RGB images: characterizing the variability of yield components from high-throughput plant phenotyping. 8th International Plant Phenotyping Symposium. Presented at the 8th International Plant Phenotyping Symposium, Lincoln, Nebraska, USA.

Dandrifosse, S., Carlier, A., Dumont, B., Mercatoris, B., 2021. Registration and fusion of close-range multimodal wheat images in field conditions. Remote Sens. 13, 1380. https://doi.org/10.3390/rs13071380.

De Solan, B., Baret, F., Thomas, S., Weiss, M., Beauchene, K., Comar, A., Fournier, A., 2015. Development and use of a fully automated PHENOMOBILE for field phenotyping. Presented at the EPPN Plant Phenotyping Symposium.

Debaque, B., Perreault, H., Mercier, J.-P., Drouin, M.-A., David, R., Chatelais, B., Duclos-Hindié, N., Roy, S., 2022. Thermal and visible image registration using deep homography. In: International Conference on Information Fusion (FUSION). Presented at the 2022 25th International Conference on Information Fusion (FUSION), pp. 1–8, 10.23919/FUSION49751.2022.9841256.

Dias Junior, J.D., Backes, A.R., Escarpinati, M.C., Pinto Silva, L.H.F., Costa, B.C.S., Avelar, M.H.F., 2020. UAV-m multispectral sensed data band co-registration framework. In: Int. Conf. Syst. Signals Image Process. 2020-July, pp. 223–228, 10.1109/IWSSIP48289.2020.9145079.

Douarre, C., Crispim-Junior, C., Gelibert, A., Rousseau, D., Tougne, L., 2019. A strategy for multimodal canopy images registration: Case study. 7th International Workshop on Image Analysis Methods in the Plant Sciences. Lyon, France.

Duan, Z., Jung, C., 2022. Joint disparity estimation and pseudo NIR generation from cross spectral image pairs. IEEE Access 10, 7153–7163. https://doi.org/10.1109/ACCESS.2022.3140435.

Fernández, C.I., Haddadi, A., Leblon, B., Wang, J., Wang, K., 2021. Comparison between three registration methods in the case of non-georeferenced close range of multispectral images. Remote Sens. 13, 396. https://doi.org/10.3390/rs13030396.

Genser, N., Spruck, A., Seiler, J., Kaup, A., 2020. Deep learning based cross-spectral disparity estimation for stereo imaging. In: Proc. IEEE International Conference on Image Processing (ICIP).

Grifoni, E., Gargano, M., Melada, J., Interlenghi, M., Castiglioni, I., Gosetti, R., di Sturmeck, S., Ludwig, N., 2022. Documenting cultural heritage in very hostile friction contexts: the synoptic visualization of giottesque frescoes by multispectral and 3D close-range imaging. J. Phys. Conf. Ser. 2204, 012060. https://doi.org/10.1088/1742-6596/2204/1/012060.
Guo, X., Yang, K., Yang, W., Wang, X., Li, H., 2019. Group-wise correlation stereo network. In: 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR). IEEE, Long Beach, CA, USA, pp. 3268-3277, 10.1109/ CVPR.2019.00339.

Heo, Y.S., Lee, K.M., Lee, S.U., 2013. Joint depth map and color consistency estimation for stereo images with different illuminations and cameras. IEEE Trans. Pattern Anal. Mach. Intell. 35, 1094-1106. https://doi.org/10.1109/TPAMI2012.167.

Hirschmüller, H., 2008. Stereo processing by semiglobal matching and mutual information. IEEE Trans. Pattern Anal. Mach. Intell. 30, 328-341. https://doi.org/ 10.1109/TPAMI.2007.1166.

Hirschmüller, H., Scharstein, D., 2009. Evaluation of stereo matching costs on images with radiometric differences. IEEE Trans. Pattern Anal. Mach. Intell. 31, 1582-1599. https://doi.org/10.1109/TPAMI.2008.221.

Hirschmüller, H., Scharstein, D., 2007. Evaluation of cost functions for stereo matching. Proc. IEEE Comput. Soc. Conf. Comput. Vis. Pattern Recognit, 10.1109/ CVPR.2007.383248.

Honkavaara, E., Saari, H., Kaivosoja, J., Pölönen, I., Hakala, T., Litkey, P., Mäkynen, J., Pesonen, L., 2013. Processing and assessment of spectrometric, stereoscopic imagery collected using a lightweight UAV spectral camera for precision agriculture. Remote Sens. 5, 5006-5039. https://doi.org/10.3390/rs15105006.

Hunt, E.R., Hively, W.D., Fujikawa, S.J., Linden, D.S., Daughtry, C.S.T., McCarty, G.W., 2010. Acquisition of NIR-green-blue digital photographs from unmanned aircraft for crop monitoring. Remote Sens. 2, 290-305.

Jay, S., Baret, F., Dutartre, D., Malatesta, G., Héno, S., Comar, A., Weiss, M., Maupas, F., 2019. Exploiting the centimeter resolution of UAV multispectral imagery to improve remote-sensing estimates of canopy structure and biochemistry in sugar beet crops. Remote Sens. Environ. 231. https://doi.org/10.1016/j.rse.2018.09.011.

Jay, S., Comar, A., Benicio, R., Beauvois, J., Dutartre, D., Daubige, G., Li, W., Labrosse, J., Thomas, S., Henry, N., Weiss, M., Baret, F., 2020. Scoring cercospora leaf spot on sugar beet: comparison of UGV and UAV phenotyping systems. Plant Phen. 2020, 1-18. https://doi.org/10.34133/2020/9452123.

Jay, S., Gorretta, N., Morel, J., Maupas, F., Dutendoula, R., Rabatel, G., Dutartre, D., Comar, A., Baret, F., 2017. Estimating leaf chlorophyll content in sugar beet canopies using millimeter- to centimeter-scale reflectance imagery. Remote Sens. Environ. 198, 173-186. https://doi.org/10.1016/j.rse.2017.06.008.

Jay, S., Guillaume, M., Blanc-Talon, J., 2012. Underwater target detection with hyperspectral data: solutions for both known and unknown water quality. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 5, 1213-1221. https://doi.org/10.1109/ JSTARS.2012.2185488.

Jhan, J.P., Rau, J.Y., Haala, N., 2018. Robust and adaptive band-to-band image transform of UAS miniature multi-lens multispectral camera. ISPRS J. Photogramm. Remote Sens. 137, 47-60. https://doi.org/10.1016/j.isprsjprs.2017.12.009.

Kerkech, M., Hafiane, A., Canals, R., 2020. Vine disease detection in UAV multispectral images using optimized image registration and deep learning segmentation approach. Comput. Electron. Agric. 174, 105446. https://doi.org/10.1016/j. compag.2020.105446.

Kim, S., Min, D., Ham, B., Do, M.N., Sohn, K., 2017. DASC: robust dense descriptor for multi-modal and multi-spectral correspondence estimation. IEEE Trans. Pattern Anal. Mach. Intell. 39, 1712-1729. https://doi.org/10.1109/TPAMI.2016.2615619.

Krotosky, S.J., Trivedi, M.M., 2007. Mutual information based registration of multimodal stereo videos for person tracking. Comput. Vis. Image Underst. 106, 270-287. https://doi.org/10.1016/j.cviu.2006.10.008.

Laveglia, S., Altieri, G., Genovese, F., Matera, A., Scarano, L., Di Renzo, G.C., 2025. Development of a short-range multispectral camera calibration method for geometric image correction and health assessment of baby crops in greenhouses. Appl. Sci. 15, 2893. https://doi.org/10.3390/app15062893.

Lerma, J.L., Cabrelles, M., Akasheh, T.S., Haddad, N.A., 2012. Documentation of weathered architectural heritage with visible, near Infrared, thermal and laser scanning data. Int. J. Herit. Digit. Era 1, 251-275. https://doi.org/10.1260/2047- 4970.1.2.251.

Li, W., Comar, A., Weiss, M., Jay, S., Colombeau, G., Lopez-Lozano, R., Madec, S., Baret, F., 2021. A double swath configuration for improving throughput and Accuracy of trait estimate from UAV images. Plant Phenomics 2021, 1-11. https:// doi.org/10.34133/2021/9892647.

Liang, M., Guo, X., Li, H., Wang, X., Song, Y., 2019. Unsupervised cross-spectral stereo matching by learning to synthesize. Proc. AAAI Conf. Artif. Intell. 33, 8706-8713. https://doi.org/10.1609/aaii.v33i01.33018706.

Liu, C., Yuen, J., Torralba, A., 2011. SIFT flow: dense correspondence across scenes and its applications. IEEE Trans. Pattern Anal. Mach. Intell. 33, 978-994. https://doi. org/10.1109/TPAMI.2010.147.

Liu, H., Lee, S.H., Chah, J.S., 2018. Registration of multispectral 3D points for plant inspection. Precis. Agric. 19, 513-536. https://doi.org/10.1007/s11119-017-9536-3.

Lowe, D.G., 2004. Distinctive image features from scale-invariant keypoints. Int. J. Comput. Vis. 60, 91-110. https://doi.org/10.3390/rs14010078.

Lu, W., Okayama, T., Komatsuzaki, M., 2022. Rice height monitoring between different estimation models using UAV photogrammetry and multispectral technology. Remote Sens. 14, 78. https://doi.org/10.3390/rs14010078.

Madec, S., Baret, F., de Solan, B., Thomas, S., Dutartre, D., Jezequel, S., Hemmerlé, M., Colombeau, G., Comar, A., 2017. High-throughput phenotyping of Plant height: Comparing unmanned aerial vehicles and ground LiDAR estimates. Front. Plant Sci. 8.

Manolakis, D., Marden, D., Shaw, G.A., 2003. Hyperspectral image processing for automatic target detection applications. Linc. Lab. J. 14, 79-116.

Manolakis, D., Shaw, G., 2002. Detection algorithms for hyperspectral imaging applications. IEEE Signal Process Mag. 19, 29-43.

Meng, L., Zhou, J., Liu, S., Ding, L., Zhang, J., Wang, S., Lei, T., 2021. Investigation and evaluation of algorithms for unmanned aerial vehicle multispectral image registration. Int. J. Appl. Earth Obs. Geoinform. 102, 102403. https://doi.org/ 10.1016/j.jag.2021.102403.

Moorthy, I., Miller, J.R., Noland, T.L., 2008. Estimating chlorophyll concentration in conifer needles with hyperspectral data: an assessment at the needle and canopy level. Remote Sens. Environ. 112, 2824-2838. https://doi.org/10.1016/j. rse.2008.01.013.

O'Reilly, J.E., Maritorena, S., Mitchell, B.G., Siegel, D.A., Carder, K.L., Garver, S.A., Kahru, M., McClain, C., 1998. Ocean color chlorophyll algorithms for SeaWiFS. J. Geophys. Res. Oceans 103, 24937-24953. https://doi.org/10.1029/98JC02160.

Pinggera, P., Breckon, T., Bischof, H., 2012. On cross-spectral stereo matching using dense gradient features. BMVC 2012 - electron. Proc. Br. Mach. vis. Conf. 2012, 1-12. https://doi.org/10.5244/C.26.26.

Pozo, S.D., Lindenbergh, R., Rodríguez-González, P., Blom, J.K., González-Aguilera, D., 2015. Discrimination between Sedimentary rocks from close-range visible and very- near-infrared images. PLoS One 10, e0132471. https://doi.org/10.1371/journal. pone.0132471.

Rabatel, G., Gorretta, N., Labbe, S., 2014. Getting simultaneous red and near-infrared band data from a single digital camera for plant monitoring applications: theoretical and practical study. Biosyst. Eng. 117, 2-14.

Rabatel, G., Labbé, S., 2016. Registration of visible and near infrared unmanned aerial vehicle images based on Fourier-Mellin transform. Precis. Agric. 17, 564-587. https://doi.org/10.1007/s11119-016-9437-x.

Rana, S., Gerbino, S., Crimaldi, M., Cirillo, V., Carillo, P., Sarghini, F., Maggio, A., 2024. Comprehensive evaluation of multispectral image registration strategies in heterogeneous agriculture environment. J. Imaging 10, 61. https://doi.org/10.3390/ jimaging10030061.

Rincón, M.G., Mendez, D., Colorado, J.D., 2022. Four-dimensional Plant phenotyping model integrating low-density LiDAR data and multispectral images. Remote Sens. 14, 356. https://doi.org/10.3390/rs14020356.

Rohlfing, T., 2012. Image Similarity and tissue overlaps as surrogates for image registration Accuracy: widely used but unreliable. IEEE Trans. Med. Imaging 31, 153-163. https://doi.org/10.1109/TMI.2011.2163944.

Rueckert, D., Sonoda, L.I., Hayes, C., Hill, D.L.G., Leach, M.O., Hawkes, D.J., 1999. Nonrigid registration using free-form deformations: application to breast MR images. IEEE Trans. Med. Imaging 18, 712-721. https://doi.org/10.1109/42.796284.

Scharstein, D., Szeliski, R., 2002. A taxonomy and evaluation of dense two-frame stereo correspondence algorithms. Int. J. Comput. Vis. 47, 7-42.

Sharma, N., Banerjee, B.P., Hayden, M., Kant, S., 2023. An open-source package for thermal and multispectral image analysis for plants in glasshouse. Plants 12, 317. https://doi.org/10.3390/plants12020317.

Shen, H.L., Zou, Z., Zhu, Y., Li, S., 2019. Block-based multispectral image registration with application to spectral color measurement. Opt. Commun. 451, 46-54. https:// doi.org/10.1016/j.optcom.2019.06.041.

Sun, G., Wang, X., Sun, Y., Ding, Y., Lu, W., 2019. Measurement method based on multispectral three-dimensional imaging for the chlorophyll contents of greenhouse tomato plants. Sensors 19, 3345. https://doi.org/10.3390/s19153345.

Szeliski, R., 2010. Computer Vision : Algorithms And Applications. Springer Science & Business Media.

Tosi, F., Ramirez, P.Z., Poggi, M., Salti, S., Mattoccia, S., Di Stefano, L., 2022. RGB- multispectral matching: dataset, Learning methodology, evaluation. In: Presented at the Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 15958-15968.

Trajković, M., Hedley, M., 1998. Fast corner detection. Image Vis. Comput. 16, 75-87. https://doi.org/10.1016/S0262-8856(97)00056-5.

Tusa, E., Laybros, A., Monnet, J.-M., Dalla Mura, M., Barré, J.-B., Vincent, G., Dalponte, M., Férét, J.-B., Chanussot, J., 2019. Chapter 2.11 - Fusion of hyperspectral imaging and LiDAR for forest monitoring. In: Amigo, J.M. (Ed.), Data Handling in Science and Technology, Hyperspectral Imaging. Elsevier, pp. 281-303. DOI: 10.1016/B978-0-444-63977-6.00013-4.

Vayssade, J.A., Jones, G., Paoli, J.N., Gee, C., 2020. Two-step multi-spectral registration via key-point detector and gradient similarity: application to agronomic scenes for proxy-sensing. VISIGRAPP 2020 - proc. 15th int. jt. conf. Comput. vis. Imaging Comput. Graph. Theory Appl. 4, 103-110. https://doi.org/10.5220/ 00091693010301100.

Vuletić, J., Car, M., Orsag, M., 2023. Close-range multispectral imaging with multispectral-depth (MS-D) system. Biosyst. Eng. 231, 178-194. https://doi.org/ 10.1016/j.biosystemseng.2023.06.002.

Vuletić, J., Polić, M., Orsag, M., 2022. Introducing multispectral-depth (MS-D): sensor fusion for close range multispectral imaging. In: 2022 IEEE 18th International Conference on Automation Science and Engineering (CASE). Presented at the 2022 IEEE 18th International Conference on Automation Science and Engineering (CASE), pp. 1218-1223, 10.1109/CASE49997.2022.9926611.

Wang, X., Yang, W., Wheaton, A., Cooley, N., Moran, B., 2010. Efficient registration of optical and IR images for automatic plant water stress assessment. Comput. Electron. Agric. 74, 230-237. https://doi.org/10.1016/j.compag.2010.08.004.

Xie, P., Du, R., Ma, Z., Cen, H., 2023. Generating 3D multispectral point clouds of plants with fusion of snapshot spectral and RGB-D images. Plant Phenomics 5, 0040. https://doi.org/10.34133/plantphenomics.0040.

Yasir, R., Eramian, M., Stavness, I., Shirliffe, S., Duddu, H., 2018. Data-driven multispectral image registration. proc. - 2018 15th conf. Comput. Robot vis. CRV 2018, 230-237. https://doi.org/10.1109/CRV.2018.00040.

Ye, Y., Shan, J., 2014. A local descriptor based registration method for multispectral remote sensing images with non-linear intensity differences. ISPRS J. Photogramm. Remote Sens. 90, 83-95. https://doi.org/10.1016/j.isprsjprs.2014.01.009.
Zabih, R., Woodfill, J., 1994. Non-parametric local transforms for computing visual correspondence. In: European Conference on Computer Vision. Springer, pp. 151–158.

Zarco-Tejada, P.J., Miller, J.R., Morales, A., Berjón, A., Agüera, J., 2004. Hyperspectral indices and model simulation for chlorophyll estimation in open-canopy tree crops. *Remote Sens. Environ.* **90**, 463–476. https://doi.org/10.1016/j.rse.2004.01.017.

Zbontar, J., LeCun, Y., 2016. Stereo matching by training a convolutional neural network to compare image patches. *J. Mach. Learn. Res.* **17**, 2287–2318.

Zhao, Y., Yue, T., Chen, L., Wang, H., Ma, Z., Brady, D.J., Cao, X., 2017. Heterogeneous camera array for multispectral light field imaging. *Opt. Express* **25**, 14008. https://doi.org/10.1364/oe.25.014008.

Zhi, T., Pires, B.R., Hebert, M., Narasimhan, S.G., 2018. Deep material-aware cross-spectral stereo matching. In: *Proc. IEEE Comput. Soc. Conf. Comput. Vis. Pattern Recognit.*, pp. 1916–1925 doi: 10.1109/CVPR.2018.00205.

Zitova, B., Flusser, J., 2003. Image registration methods: a survey. *Image Vis. Comput.* **21**, 977–1000.