Introduction

The light spectrum is a continuous intensity distribution across wavelengths. This spectral information is commonly used to help determine and/or discriminate the physical properties of object surfaces, for example in remote sensing [25,8,14,23] and medical imaging [28,29]. Also, in various practical applications, the devices (sensors or displays), light sources and object surfaces are characterized by spectral measurements [1,9,16,27].

Despite the advantages of spectral capture, almost all images that we record contain just 3 measurements -the 3 weighted spectral averages over the Red, Green and Blue spectral regions. Perforce, much spectral information is therefore lost in this RGB image formation process. Indeed, it is a classical result in color science, that there are many spectra -called metamers [12] -which integrate to the same RGB, and of course, given only one RGB measurement we cannot tell which physical spectrum induced it. Still, by adopting learning approaches we can estimate the spectrum that likely corresponds to a given RGB. Estimating spectra from RGBs is called spectral reconstruction (SR). In Fig. 1 we illustrate RGB image formation and the SR process. In the top-left panel we see a single radiance spectrum measured at one location in the hyperspectral image (bottom left). This spectrum is sampled by 3 sensors, resulting in the 3-value RGB response (top-right). Repeating this process for all image locations, the corresponding RGB image in the bottom right is derived from the hyperspectral image. Then, the spectral reconstruction algorithms attempt to recover the hyperspectral image back from the RGB image (or an approximation thereof).

Historically, this SR problem is effectively solved by least-squares regression [15,11,17,18,2], where the map from RGBs (or the non-linear RGB features) to spectra is modelled as a simple linear transform. More recently, deep learning approaches [4,22,6] have been developed that provide even better SR performance. Effectively, this performance increment is achieved by regressing an RGB in the context of its neighborhood to its corresponding spectra. Clearly, this patch-based idea has merit. For example, if the algorithm can identify a patch in the scene as a 'skin region' then spectral recovery is plausibly easier to solve, since skin spectra have characteristic spectral shapes [7,19].

Despite clear rationale behind the deep-learning approach, Aeschbacher et al. [2] show that the regression-based A+ algorithm provides very competitive performance. Moreover, Lin and Finlayson [17] shows that several regression methods actually generalize better than the leading deep-learning models when the scene exposure changes.

The main concern of this paper is in the regularization step of regressionbased SR algorithms. The classical (multivariate) regression problem from statis-tics is written as

MA ≈ B ,(1)

where A is an m×N matrix as the table of measured data (m is the dimension of the measured data and N is the number of data samples, with N m), and B is the corresponding target data matrix, of dimension k × N (k is the dimension of the target data). The aim is to find the k × m linear mapping M that makes the approximation as good as possible. Now, let us suppose small fluctuations in the target data, denoted as a matrix E of very small numbers (all entries are close to 0). The following regression is almost identical as Equation ( 1):

M A ≈ B + E . (2)

And yet, often we find that the best solution for M and M are very different from one another. The reason for this is that some dimensions of the measured data (the rows of the data matrix A) could be highly correlated, such that there can be very different M's that fit B equally well.

Regularization theory [24] is a way of dealing with this kind of non-robustness. Given the example of Equation ( 1) and ( 2), we may ask that among all plausible (near optimized) solutions, which one is more likely to generalize to the 'unseen data' better. Typically, the principle of regularization follows the idea that the best fitting function (i.e. M) should be the simplest possible solution that can still fits the data well.

In Fig. 2 we show a 1-D example. In the least-squares sense, the wiggly red curve is found to best fit the training data (black data points). Yet intuitively, this is not the correct fit, as the data points appear to follow a much simpler distribution. In contrast, the regularized fit (blue curve) seems to model the data better.

Fig. 2. Example of a least-squares fit (red curve) and a regularized least-squares fit (blue curve). The black data points are the training data. In the least-squares sense, the overall distance between the data points and the red wiggly fit is less than the blue smooth fit, but the latter looks 'more correct'.

Returning to our spectral reconstruction problem, taking linear regression as an example, the matrix A corresponds to a set of image RGBs (m = 3), and B refers to the spectra we are trying to recover (k is the number of spectral channels we have measured). The insight that we explore in this paper is that there is no reason why the fits for different spectral channels should be regularized altogether. Rather, we seek to consider the per-channel regression problem, where each spectral channel is recovered in turn, and correspondingly, each row of M is solved in turn. This simple modification allows us to carry out a per-channel regularization that ensures individual optimizations for all spectral channels in the spectral reconstruction problem.

Either for the conventional global regularization or for our per-channel approach, care must be taken not to overly tune the terms in M to the data at hand. This led us to develop a modified cross validation procedure. Our method separates the data on hand into three subsets, respectively for training, regularization and testing, which is a novel adjustment from the standard methodology [20] and is another contribution of this paper.

Background

Hyperspectral Images and RGB Simulation

In a hyperspectral image, spectra are measured discretely at some sampled wavelengths. Suppose the visible spectrum runs from 400 to 700 nanometers and the spectral sampling is every 10 nanometers, we get a 31-dimensional discrete representation of spectra, denoted as r ∈ R 31 .

Correspondingly, the spectral sensitivities of the R, G and B camera sensors can also be represented in discrete vector form (i.e. as 31-dimensional vectors), respectively denoted as s R , s G and s B . Then, as per the illustration in Fig. 1 we can write image formation as [26]:

x =    R G B    =    s T R s T G s T B    r ,(3)

where x = (R, G, B) T is the 3-value RGB camera response.

In the SR problem (the bottom of Fig. 1) we seek to recover hyperspectral images from the RGB images. Denote an SR algorithm as

Ψ : R 3 → R 31 , Ψ (x) ≈ r . (4)

Regression-Based Spectral Reconstruction

The general regression-based formulation of the SR problem is written as

Ψ (x) = Mϕ(x) ,(5)

where ϕ(•) is a feature transformation which maps each RGB to a corresponding p-term feature vector, and in turn it is mapped by the regression matrix M. Each of the various regression-based models [15,11,17,18,2] adopts a bespoke definition of ϕ(x). For details of the considered models, including Linear, Root-Polynomial and Adjusted Anchored Neighborhood Regression [15,17,2], see Appendix A.

Least-Squares Optimization The most common least-squares optimization seeks to minimize the sum of squared errors between the ground-truth training spectral data and the reconstruction from their corresponding RGBs: Ψ (x).

Given the formulation of Ψ (x) in Equation ( 5), the least-squares optimization of M is formulated as:

M = arg min M N i=1 ||r i − Mϕ(x i )|| 2 2 , (6)

where N is the number of data points in the training set and i indexes an individual spectrum. Collating all spectral training data in a data matrix R = (r 1 , r 2 , ..., r N ) and the corresponding feature vector matrix Φ = (ϕ(x 1 ), ϕ(x 2 ), ..., ϕ(x N )), Equation ( 6) can be written as:

M = arg min M ||R − MΦ|| 2 F .(7)

Here || • || 2 F is the squared Frobenius norm, which is exactly the sum-of-squares of all entries of the enclosed matrix.

Tikhonov Regularization

In regression-based SR, the most common method to regularize a model is Tikhonov Regularization [15,24], which hypothesizes that a more natural fit is obtained when the 'magnitude' (or the 'matrix norm') of M is bounded to some extent. Based on this assumption, the least-squares optimization in Equation ( 7) is extended to incorporate a regularization term:

M = arg min M ||R − MΦ|| 2 F + γ||M|| 2 F .(8)

Here, the ||M|| 2 F term (the regularization term, or penalty term) is controlled by a user-defined regularization parameter γ ≥ 0, which is usually determined empirically [13,21].

Equation ( 8) is solved in closed form [15]:

M = RΦ T (ΦΦ T + γI) −1 , (9)

where I is the p × p identity matrix (recall that p is the dimension of the feature vectors ϕ(x)).

Proposed Method

In spectral reconstruction, we wish to recover spectral measurements in the range from 400 to 700 nanometers (the visible spectrum). Suppose we know the intensity of light entering the camera at 400 nanometers. Given this knowledge if we wished to predict the value of the spectrum at 410 nanometers, it makes sense to assume a similar value as the one at 400 nanometers. Indeed, the fact that the intensity values at close-by wavelengths are similar is why we can represent the continuous visible spectrum at discrete wavelengths. Conversely, one could not use the knowledge of light at 400 nanometers to predict the spectral value at, say, 700 nanometers. And yet, in the literature when we regularize the regression-based SR models, we are -in some sense -assuming that all wavelengths are related. Our new per-channel reformulation of Tikhonov regularization for spectral reconstruction effectively allows the recovery of the spectral values at distant wavelengths to be considered more independently from one another.

Per-Channel Regularization

Let us split the regression matrix M by row: M = (m 1 , m 2 , ..., m 31 ) T , such that the general form of regression-based SR formulated in Equation (5) becomes

Ψ (x) = Mϕ(x) =       m T 1 m T 2 . . . m T 31       ϕ(x) =       r1 r2 . . . r31       ,(10)

where (r 1 , r2 , ..., r31 ) T = r is the reconstructed spectrum. For an arbitrary k th spectral channel, the estimated intensity value rk is given by rk = m T k ϕ(x). Note that as we represent the regression model by channel, we do not alter the original model. This says that the regression-based spectral reconstruction has always been in such a way that the reconstruction for each spectral channel depends exclusively on the corresponding row of M.

Given this fact, we might expect that each row of M would be optimized independently. However, this was not the case for the conventional regularized least-squares solution in Equation (9). Indeed, we see in Equation ( 8) the regularization term is only controlled by one single regularization parameter γ and the fits for all spectral channels (all rows of M) are regularized by the same γ. Regardless of how we optimize this γ parameter, this setting clearly contradicts to the inherent independence between the rows of M.

Let us split the spectral reconstruction problem into 31 independent problems, where the function Ψ k : R 3 → R reconstructs the k th -channel values of the reconstructed spectra by the k th row of M:

Ψ k (x) = m T k ϕ(x) .(11)

Then, we are to determine m T k , as the least-squares fit for the k th channel data. Recall the training spectral data matrix R whose columns are individual training spectra, now we split R by spectral channel instead:

R = (r 1 , r 2 , ..., r N ) =       ρ T 1 ρ T 2 . . . ρ T 31       ,(12)

where ρ T k includes the k th channel values of all training spectral data. m T k is then optimized following the regularized least-squares optimization:

m T k = arg min m T k ||ρ T k − m T k Φ|| 2 2 + γ k ||m k || 2 2 ,(13)

and solved in closed form:

m T k = ρ T k Φ T (ΦΦ T + γ k I) −1 . (14)

Here γ k represents the channel-wise regularization parameter that only controls the regularization for the k th channel. Clearly, our per-channel regularization scheme (Equation ( 14)) solves the regression matrix M row-by-row, such that each row is ready to be regularized independently. The remaining question is then how we are going to optimize these regularization parameters.

Modified Cross Validation

Perforce, regularization parameters (conventionally the single γ and our perchannel γ k 's) are determined empirically. In the literature, a grid-search approach is adopted, where different parameters are tried to regularize the model. These 'intermediate models' are then used to recover spectra from a set of unseen RGB images, and the model that minimizes the evaluation criteria is selected. For example, see [11,17,2].

As usual we would like to train, regularize and test a model using the images from the same database, we must partition the database into several subsets for these different usages. All (to our knowledge) deep-learning models simply separate the image database into 3 subsets randomly (respectively for training, validation and testing, in the parlance of deep learning), see [22,6]. However, this setting can potentially create so-called 'unfair' separations, such that if the database is separated differently, the results may vary.

A better practice is using a cross-validation process. In this paper we develop our own cross validation scheme, which is modified from the conventional Kfold cross validation [20]. This is because the conventional K-fold only seeks to separate a dataset into training and testing data, and here we need an additional partition as the regularization data.

In Fig. 3 we show the comparison between the conventional 4-fold cross validation (left) and our method (right). For both methods the same experiment is conducted 4 times. In each trial, the conventional method selects 3 out of 4 portions of data for training (marked in blue) and the remaining part is for testing (marked in orange). In our method, however, only 2 out of 4 portions of data are for training, which allows 1 portion of data (marked in green) used for regularization, that is to determine the γ and γ k parameters. Subject to these terms we solve for the best regression model for the training (blue) data. The performance statistics are calculated based on the recovery errors on the testing (orange) data and averaged over the 4 trials. Notice for our cross validation method there actually exists more possible permutations than the presented 4-trial setting. To be exact, there should be 12 different permutations. We remark that according to our empirical study, experimenting with more trials (than the presented setting) does not make significant difference in the performance statistics.

Experiments

Considered Models

In this paper we consider 3 regression-based models:

-Linear Regression (LR) [15] -Root-Polynomial Regression (RPR) [17] -Adjusted Anchored Neighborhood Regression (A+) [2] . For all the above models we adopt both the original regularization (as reported in their citations and as per Equation ( 9)) and our per-channel regularization (as per Equation ( 14)).

Database

We use the ICVL hyperspectral image database [3] (Fig. 4), which provides 201 hyperspectral images of spatial dimension 1392 × 1300 and with 31 spectral channels. The spectral channels represent narrow-band intensity measurements, respectively at every 10 nanometers (nm) between 400 and 700 nm.

The corresponding RGB images are simulated following the linear RGB simulation setting (Equation ( 3)), and the CIE 1964 color matching functions [10] are used as the spectral sensitivities.

Evaluation Criteria

The selected evaluation metric is Relative Absolute Error (RAE) [4,5], which is defined per channel as:

RAE(r k , rk ) = r k − rk r k ,(15)

where r k and rk are respectively the k th -channel values of the ground-truth and reconstructed spectra. Effectively, this metric measures the percentage absolute error. RAE is the most common performance measure used in recent research, and the rationale of using this metric can be found in [4].

Results and Discussion

In Table 1, we present the per-channel error statistics of LR (left table), RPR (middle table) and A+ (right table) under the original settings -where a single penalty term is used in the regularization [15,17,2] -and our per-channel regularization method. The spectral channels are represented by the wavelengths λ (nm). We also calculate the percentage 'improvements' as: which is presented in the rightmost column of each table. In the bottom of each table, the Mean RAE results (averaging over all spectral channels) are shown. First, we see that for all considered models, our method improves the RAE in multiple channels by over 10% (marked in bold and with underlines), with maximal improvements around 16-17%.

Improve (%) = 100 × RAE original − RAE ours RAE original ,(16)

Secondly, in terms of Mean RAE performance, our method improves the RPR model the most, by an 8.6% increment, compared to 3.2% for A+ and 3.1% for LR. Significantly, the A+ model is the leading sparse coding model, which is shown able to perform equally-well as some deep-learning solutions [2]. By improving the A+ model, we effectively bring forward the shallow-learned baseline. Moreover, our method reduces the gap between RPR and A+. Relative to A+, RPR model is much simpler (with significantly less model parameters) [17], which ensures more effective model re-training and shorter runtime.

Lastly, for the A+ model, it seems curious that the per-channel performances in the first three channels (400, 410 and 420 nm) degrade by minute differences. Indeed, this means the regularization parameters we chose for these channels are not actually optimized for the test-set data. We remark that this is most likely originated from the unequal separation of data subsets in cross validation, such that the best regularization parameter for the regularization-set data does not correspond to the best for the test-set data. We are investigating how to remedy this issue.

For one example image in the ICVL database [3], we visualize the spectral reconstruction errors as the Mean RAE error maps in Fig. 5. It is clear that for all models our method improves the Mean RAE in various parts of the image. For example, the tree stem for LR and RPR, and the grassy ground for A+.

Conclusion

In the spectral reconstruction (SR) problem, hyperspectral images are reconstructed from RGB images. Many approaches are based on least-squares regression, where the fitting function is modelled by a simple linear transformation, and a Tikhonov regularization process is applied to improve model generalizability. Conventionally, the fits for all spectral channels are jointly regularized. We demonstrate that the fit for each spectral channel can be formulated independently, such that the fit for each channel is regularized (therefore optimized) independently. We also provide a novel modification of K-fold cross validation so that the models can be fairly trained, regularized and tested with a single image database. Compared to the original models, our per-channel regularization method improves the accuracy of recovery for individual spectral channels by up to 17% increments, and by 3-9% in mean improvements over all spectral channels.

The spectral reconstruction then seeks to linearly map these root-polynomial vectors to spectra:

Ψ (x) = Mϕ α (x) .(18)

In this paper we set α = 6, which is the 6 th -order RPR.

A.3 Adjusted Anchored Neighborhood Regression (A+ Sparse

Coding)

The leading sparse coding method 'A+' [2] assumes linear maps from RGBs to spectra (effectively, operates LR in every neighborhood). Denote Ψ i (x) as the spectral reconstruction mapping for the data in the i th neighborhood. On input of an RGB x, the reconstruction is written as:

neighborhood(x) = i ⇒ Ψ i (x) = M i x .(19)

See [2] for more details about the model.