1.

Introduction

Recently, high-fidelity color image reproduction has become increasingly important because of its many potential applications in textiles, medicine, digital archives, etc. To view objects under various illumination conditions, multispectral imaging has been extensively studied to estimate spectral reflectance of object surfaces.^{1, 2, 3, 4, 5, 6, 7} Multispectral images are usually acquired by trichromatic or monochrome cameras, accompanied by a set of color filters. The imaging process is often modeled by a linear system when the number of imaging channels is large. In recovery of spectral reflectance, Wiener estimation is deduced under the condition that the data conform to normal distribution.⁸ It works fairly well when the linearity condition is satisfied, and becomes a standard technique. However, it is likely that measured data are not in accordance to the normal distribution,^{2, 3} and the optoelectronic conversion function (OECF) of the camera is nonlinear.⁹ Nonlinearity usually degrades the estimation accuracy of linear methods. Adaptive methods, which use local statistics instead of global ones, provide feasible solutions for improving accuracy.^{2, 4} However, as adaptive methods usually need to recalculate the transform between responses and reflectance,⁴ they are computationally expensive and hence unsuitable for time-critical applications. An alternative way is to introduce nonlinear variables such as high-order polynomials. Hong, Luo, and Rhodes¹⁰ applied ordinary polynomial regression to predict colorimetric stimulus values from three-channel camera responses. The same technique has also been adopted in multispectral imaging.⁵ However, the extension of polynomial responses causes overfitting and collinearity problems when the number of imaging channels is large. Heikkinen ⁶ introduced regularized polynomial modeling methods and a more general regularization framework for robust reflectance estimation.

We propose a global method for spectral estimation based on polynomial extension of camera responses and partial least-squares (PLS) regression.^{11, 12} The PLS is implemented in an iterative manner; its dimension (or number of PLS components) is determined by the spectral error distribution. The accuracy of the PLS method is compared with Wiener estimation and also polynomial regressions solved by ordinary least squares (OLS) and regularized least squares (RLS).

2.

Previous Methods

Suppose that the continuous visible spectrum is uniformly sampled at $N$ (usually $N=31$ ) discrete wavelengths, and the number of imaging channels is $C=16$ . Let $\mathcal{F}(\cdot )$ be the OECF of camera, $\mathbf{r}∊{\mathcal{R}}^{N\times 1}$ be the reflectance of imaged object surface, and $\mathbf{M}∊{\mathcal{R}}^{C\times N}$ be the spectral responsivity of the imaging system, then camera response $\mathbf{u}∊{\mathcal{R}}^{C\times 1}$ is computed as

By ignoring $\mathcal{F}(\cdot )$ , the linear method tries to find a transform matrix $\mathbf{W}∊{\mathcal{R}}^{N\times C}$ such that the estimate of reflectance can be computed as

$\mathbf{W}$ can also be solved under the least-squares criterion. Let $K$ be the number of training samples. We can construct reflectance matrix $\mathbf{R}=[\mathbf{r}\left(1\right),\mathbf{r}\left(2\right),\dots ,\mathbf{r}\left(K\right)]∊{\mathcal{R}}^{N\times K}$ and response matrix $\mathbf{U}=[\mathbf{u}\left(1\right),\mathbf{u}\left(2\right),\dots ,\mathbf{u}\left(K\right)]∊{\mathcal{R}}^{C\times K}$ , and then calculate the transform matrix as^{1, 10}

Nonlinearity can degrade the performance of the linear reflectance estimation methods. It is natural to define a two-order polynomial response vector $\stackrel{\mathrm{̃}}{\mathbf{u}}∊{\mathcal{R}}^{J\times 1}$ to deal with the nonlinearity due to OECF and non-Gaussian data distribution:

Alternatively, $\mathbf{W}$ can also be computed by introducing a regularization (or penalization) term as

3.

Partial Least-Squares-Based Reflectance Estimation

As mentioned, the overfitting problem can occur in polynomial regression solved by OLS when the number of parameters in the mathematical model is greater than the number of dimensions of data variation. The overfitting problem is probably related to the increasing colinearity between the extended polynomial responses. In this regard, we propose to deal with it by dimensionality reduction using the PLS technique.¹¹

The polynomial response matrix $\stackrel{\mathrm{̃}}{\mathbf{U}}$ can be decomposed into a score matrix $\mathbf{T}∊{\mathcal{R}}^{K\times L}$ and a loading matrix $\mathbf{P}∊{\mathcal{R}}^{J\times L}$ , with $L$ being the number of PLS components, as

The goal of PLS is to extract the common structure between ${\stackrel{\mathrm{̃}}{\mathbf{U}}}^{\mathsf{T}}$ and ${\mathbf{R}}^{\mathsf{T}}$ by searching a projection such that the covariance between the score matrices $\mathbf{T}$ and $\mathbf{D}$ is maximized. In matrix form, this relationship is written as

The PLS algorithm is carried out in an iterative manner.¹² To obtain an orthogonal score matrix $\mathbf{T}$ , a weight matrix $\mathbf{G}∊{\mathcal{R}}^{J\times L}$ is introduced in the iterative procedure. Let $j$ be the iteration index, and $\mathbf{g}$ , $\mathbf{t}$ , $\mathbf{q}$ , and $\mathbf{d}$ be the $j$ ’th column vectors of matrices $\mathbf{G}$ , $\mathbf{T}$ , $\mathbf{Q}$ , and $\mathbf{D}$ , respectively. Before starting the iteration, let $\mathbf{E}={\mathbf{U}}^{\mathsf{T}}$ and $\mathbf{F}={\mathbf{R}}^{\mathsf{T}}$ . Matrices $\mathbf{E}$ and $\mathbf{F}$ are then column centered and normalized so that each variable has zero mean and unit variance. Let $j=1$ and $\mathbf{d}$ be any column of $\mathbf{F}$ , then compute Eqs. 11, 12, 13, 14 iteratively:

The residual matrices $\mathbf{E}$ and $\mathbf{F}$ needed for the next iteration are calculated as

Note that Eqs. 15, 16 remove the variance associated with the obtained score and loading vectors before the next iteration. If $j<L$ , let $j=j+1$ and continue the iteration starting from Eq. 11; otherwise, stop the iteration and compute the PLS transform¹²

As $L$ controls the number of iterations, its value is influential to PLS. If $L=\min (J,K)$ , $\mathbf{E}$ and $\mathbf{F}$ become zeros and PLS reduces to OLS; otherwise, if $L<\min (J,K)$ , the colinearity of matrix $\stackrel{\mathrm{̃}}{\mathbf{U}}$ is reduced.¹² In this work, $L$ is determined according to the spectral accuracy of the reflectance estimation, as is discussed in the following section.

4.

Experiment

In the multispectral imaging system, we used a monochrome digital camera (model Cool-SNAP HQ2, Roper Scientific Incorporated, Ottobrunn, Germany) with $14\text{-bit}$ digitization and 16 narrowband filters ( $10\mathrm{nm}$ half-width, product of Andover Company, Salem, New Hampshire) that uniformly cover the visible spectrum ranging from $400\text{to}700\mathrm{nm}$ . The response of the camera deviates from linearity by less than 1%.¹⁴ As it is difficult to accurately acquire OECF at such a nonlinear level by imaging a number of gray samples, we consider it appropriate to treat it by extending the polynomial camera responses, as discussed in previous sections. We used 414 textile Pantone patches as the color targets, with half for training and half for testing. The surfaces of these patches contain weak textures and some gloss reflection, and hence are not ideal diffusers. The reflectance data of these patches were measured by the spectrophotometer GretagMacBeth (Grand Rapids, Michigan) 7000A; the 16-channel multispectral images were acquired by the imaging system. To reduce imaging noise, three sequential images were captured and averaged for the same scene, and the responses of each color sample were averaged in spatial areas with approximately $60\times 60\text{pixels}$ . The transform between camera responses and reflectance was calculated from the training set and evaluated on the test set. The estimation accuracy was examined by spectral root-mean-square (rms)¹ error and CIEDE2000 color difference¹⁵ error under CIE standard illuminants. We tried different randomization strategies for selecting training and test samples, and found that the estimation accuracies were quite close. Hence we only present the experimental results of the case where the odd numbered samples were used for training and the rest for testing.

The performance of the PLS method is compared with the OLS and RLS methods. For the RLS method, we set the regularization parameter $\eta =0.001$ , which approximately produces the minimum spectral rms error. For the PLS method, it is found that the spectral error is not very sensitive to the number of PLS component $L$ . The suitable value of $L$ is in the range from 40 to 60, and we adopted $L=40$ for computational efficiency. Table 1 gives the spectral and colorimetric errors of the PLS method, compared with the Wiener estimation, OLS, and RLS methods. It is clear that in terms of both colorimetric and spectral error metrics, the proposed method outperforms the Wiener estimation and OLS methods while being close to the RLS method. This is expected, as Wiener estimation cannot account for nonlinearity, and the OLS method has the inherent overfitting problem. The RLS method exhibits improved estimation accuracy through the introduction of the regularization term. By removing the common variance in the response and reflectance matrices, PLS also yields good accuracy.

Table 1

Spectral rms errors and CIEDE2000 (ΔE00) errors of Wiener estimation, OLS, RLS, and PLS methods.

Due to the iterative nature of PLS, the computation time of the training procedure (calculation of ${\mathbf{W}}_{\mathrm{PLS}}$ ) is approximately proportional to $L$ . As the transform matrix ${\mathbf{W}}_{\mathrm{PLS}}$ is calculated before reflectance estimation, this computation time does not affect algorithm efficiency. The computation time of the test procedure is determined by the size of the matrix ${\mathbf{W}}_{\mathrm{PLS}}$ and is irrelevant to $L$ . The PLS method was implemented under the MATLAB^® environment and run on a PC with an Intel Core 2 CPU at $1.86\mathrm{GHz}$ and with $3\text{-}\mathrm{GB}$ memory. For a multispectral image with $1392\times 1040\text{pixels}$ , the PLS method costs about $9\mathrm{s}$ and can be much faster if programmed using C language. This indicates that the computational efficiency of the PLS method is acceptable to many practical applications.

5.

Conclusions

We propose a method for estimating reflectance from multichannel camera responses based on high-order polynomials and partial least squares. The proposed method is capable of dealing with nonlinearity in the imaging process. The appropriate number of PLS components is determined based on spectral rms error distribution. In terms of spectral and colorimetric error metrics, the proposed technique is superior to Wiener estimation and polynomial regression solved by ordinary least squares, and is close to polynomial regression solved by regularized least squares.