arXiv:1802.08793v1 [cs.CV] 24 Feb 2018

Multispectral Image Intrinsic Decomposition via Low Rank Constraint
Qian Huang
Nanjing University

Weixin Zhu
Nanjing University

Yang Zhao
Nanjing University

cianhwang@gmail.com

njuzhwx@163.com

zhaoyang@smail.nju.edu.cn

Linsen Chen
Nanjing University

Yao Wang
New York University

Tao Yue
Nanjing University

Xun Cao
Nanjing University

njucls@163.com

yw523@nyu.edu

yuetao@nju.edu.cn

caoxun@nju.edu.cn

Abstract

for images in general. This paper addressed the problem
of the recovery of reflectance and shading from a single
multspectral image, namely, the Intrinsic Image Decomposition problem of a whole multispectral image captured under general spectral illumination, hereafter referred to as the
IID problem. This problem is worth exploring since geometry and color information are useful under certain circumstances, but one of them always interferes the detection of
the other one. Unfortunately, growing dimension of data
makes this problem harder to cope with.

Multispectral images contain many clues of surface
characteristics of the objects, thus can be widely used in
many computer vision tasks, e.g., recolorization and segmentation. However, due to the complex illumination and
the geometry structure of natural scenes, the spectra curves
of a same surface can look very different. In this paper,
a Low Rank Multispectral Image Intrinsic Decomposition
model (LRIID) is presented to decompose the shading and
reflectance from a single multispectral image. We extend the
Retinex model, which is proposed for RGB image intrinsic
decomposition, for multispectral domain. Based on this, a
low rank constraint is proposed to reduce the ill-posedness
of the problem and make the algorithm solvable. A dataset
of 12 images is given with the ground truth of shadings and
reflectance, so that the objective evaluations can be conducted. The experiments demonstrate the effectiveness of
proposed method.

The low rank constraint that we proposed is based on
the low rank prior, or low rank nature of both shading and
reflectance images. According to the inherent nature of
the multispectral image, we derive shading basis on the
knowledge of illumination condition and derive reflectance
subspace bases by means of principle component analysis
(PCA) of Munsell color board. Assuming that the basis of
Retinex theory would continue to take effect on multispectral domain, we proposed an low-rank-based model so that
deriving reflectance and shading from a multispectral image can be modelled as an convex optimization problem.
In a significant departure from the conventional approaches
which operate in the logarithmic domain, we directly operate on the image domain to avoid adding additional noise
or breaking low rank nature. The flowchart of our proposed
algorithm for LRIID is shown in Fig. 1.

1. Introduction
The observed spectrum of a single pixel is determined
by illumination, reflectance and shading. Shading image
contains illumination condition and geometry information,
while reflectance image contains the color information and
material reflectance property, which are invariable to light
condition and shadow effect. The decomposition problem
has been a long standing problem of assorted areas such as
both computer graphics and computer vision applications.
For instance, shape-from-shading algorithms could benefit
from an image with only shading effects, while image segmentation would be easier in a world without cast shadows.
Obviously, intrinsic image decomposition is an ill-posed
problem, since there are more unknowns than observations.
In order to solve this problem, many work [29, 30, 12] focus on sparse representation spatially, but this does not hold

Suffering from the lack of ground truth data of shading
and reflectance, we provide a ground-truth dataset for multispectral intrinsic images, which provides us a possible way
to judge the quality of the decomposition results. Quantitative and qualitative experiments on our dataset have demonstrated that the performance of our work are better than
prior work in multispectral domain. Our work can bring
merits to multiple applications, such as recolorization, relighting, scene reconstruction and image segmentation.
Our major contribution can be summarized as follows:
(1) we extend the Retinex model to multispectral image in1

trinsic decomposition problem, and propose a low rank constraint to handle the ill-posedness of the problem; (2) we
provide a ground-truth dataset for multispectral intrinsic images, which can facilitate future evaluation and comparison
of multispectral image intrinsic decomposition algorithms;
(3) the proposed method achieves promising results on various of images both quantitatively and qualitatively.

2. Related Work
Intrinsic Image Decomposition. The problem of Intrinsic Image Decomposition (IID) was first introduced by Barrow and Tenenbaum [4]. The reflectance describes the
illumination-invariant albedo of the surface, while the shading contains surface geometric information and illumination
condition.
Assorted methods take advantage of additional information, including images sequences [33, 22] and videos [24]
to avoid shadow effect in poor lighting condition. With
the improvement of sensing devices like kinect, depth
cue [2, 11, 24] or surface normal [27] have been applied
to strengthen their assumption. More recently, Bousseau et
al. [8] proposed a user-assisted method to further improve
the result of separation.
Lots of methods with single input image are also proposed for the separation task. Bell et al. [6] developed a
dense conditional random field (CRF) based intrinsic image algorithm for images in the wild. Barron et al. [3]
introduced shape, illumination and reflectance from shading (SIRFS) model which performs well on images of segmented objects. Sai et al. [7] proposed L1 Image Transform
model for scene-level intrinsic decomposition. Entropy
method [13] raised by Finlayson et al. offered us a new
viewpoint to understand this problem. With the abundance
and availability of datasets and the development of computational equipment, training-based models [5, 32, 31, 37]
have been built to derive reflectance and shading from images.
An especially well-know and wide-employed model
called Retinex [23] made an assumption that the large chromatically change is generally caused by changes in reflectance rather than shading. With Retinex theory, we are
able to pinpoint where the reflectance changes in local area.
Horn, and Funt and Drew [16] analyzed local derivatives
for distinguishing between image variations that are due to
shading or reflectance. But it neglects the connection between pixels sharing the same neighborhood. On the basis
of Retinex thoery, we follow work of Chen et al. [12] to
handle intrinsic image decomposition task in multispectral
domain.
Sparse Representation. Researchers have also extended
trichromatic color constancy model to multispectral images

in order to separate reflectance and shading in higher spectral dimensions. A lot of trials have been made to explore this area. For example, Ikari et al. [18] showed us
the possibility of separating dozens of multispectral signals.
Huynh et al. [17] assumed that the scene could be segmented into several homogeneous surface patches, and were
able to estimate the illumination and reflectance spectra under the dichromatic reflectance model. In remote sensing
area, Kang et al. [20] fit multispectral data into trichromatic
model to extract features.
To overcome the drawbacks of ambiguity in local analysis, lots of research have been done to reduce the ambiguity
of both reflectance and shading. Shen et al. [29] proposed a
global optimization algorithm which combines Retinex theory and non-texture constraint to obtain global consistency
of image structures. Shen [30] further applied sparse representation of reflectance as global constraint of their observation. Material cues [26] has also been introduced. As
to multispectral images, Chen et al. [12] used super-pixel
to cut down the number of unknown parameters in this underdetermined problem. Unlike the approaches above, we
assumed that both shading and reflectance live in low dimensional subspace. the low rank nature of shading space
is widely acknowledged and exploited in prior work [16],
and the low-dimensional subspace model of reflectance is
introduced by [25, 28, 36]. With the help of training data
of reflectance bases and illumination spectra, we can solve
this problem effectively.
Dataset. As far as establishing ground-truth for intrinsic images, Tappen et al. [32] created small sets of
both computer-generated and real intrinsic images. The
computer-generated images consisted of shaded ellipsoids
with piecewise-constant reflectance. The real images were
created using green marker on crumpled paper [31]. Grosse
et al. [15] provided an all-rounded dataset which is widely
used in following analysis. Bell et al. [6] also introduced Intrinsic Images in the Wild, a large scale, public dataset for
evaluating intrinsic image decompositions of indoor scenes.
In multispectral field, [34, 10] provided a set of multispectral images of various objects without ground truth, and
Chen et al. [12] built a dataset but the number, spectral resolution and diversity of images are limited. There is no other
attempts to establish ground truth for multispectral intrinsic
images.

3. Our Model
We assume the object surface as Lambertian and hence
has diffuse reflection. In most prior work on intrinsic image decomposition, the captured luminance spectrum at every point lp is modelled as the product of Lambertian reflectance spectrum rp and shading spectrum sp , where sp

LR Shading

Original
Image

Refined
Shading

Reconstruction

Spectra
Distance

Iteration
Initial Shading

Weight
Map

Initial
Reflectance
Refined
Reflectance
Alternating
Projection

LR Reflectance
Figure 1. The flowchart of our proposed LRIID algorithm.

is used to characterize the combined effect of object geometry, illumination, occlusion and shadowing. Mathematically, this model can be expressed as

2) When the expected reflectance difference between two
pixels is small, the recovered reflectance difference between the two pixels should be small.

lp = sp . ∗ rp

First we look at whether we can separate reflectance and
shading from the measured luminance signal independently.
Take reflectance for example. By recognizing two adjacent
pixels which have the same shading, the ratio relationship
can be written as lp . ∗ rq = lq . ∗ rp , or Lp rq = Lq rp where
Lp is a diagonal matrix consisting of spectral elements in
lp , we formulate the energy functions in terms of the reflectance vectors directly:

(1)

where lp , rp and sp are all vectors with dimensions equal to
the number of spectral bands of the captured image, .∗ denotes element-wise multiplication. The problem is to derive
sp and rp from observed multispectral luminance vector lp :
In this project, we will focus on recovering the reflectance
spectrum using this model. Once rp is determined, the shading image can be derived by point-wise division.
Different from the conventional approaches which operate in the logarithmic domain, we directly formulate the
problem in the image domain, and this can overcome numerical problems caused by the logarithmic transformation
of the image values, where noise in pixels with low intensity values can lead to large variations. Besides, although
there have been substantial evidence of the low rank nature
of the reflectance space, it is not clear whether the logarithmically transformed reflectance space is still low rank. This
makes it hard to incorporate the low rank prior in formulations based on log-transformed images.

3.1. Estimate Reflectance or Shading Independently
The Retinex model makes following two important observations:
1) When there is significant reflectance change between
two adjacent pixels p and q, the shading is typically constant. This leads to the relation lp ./lq = rp ./rq , where ./
denotes element-wise division;

Esc =

P

kwp,q (Lp rq − Lq rp )kd

p,q∈Nsc

Erc = kvp,q (rp − rq )k

(2)
d

where Nsc and Nrc denote neighborhood pair sets and
wp,q and vp,q denote weights. wp,q would be large but vp,q
will be small when the expected reflectance difference between two adjacent pixels p and q are large, and vice verse.
To make the formulation general, we use d to indicate the
error norm, with d = 2 for L2 norm, and d = 1 for L1
norm. L1 norm is more difficult to solve, but the solution
can be more robust to outliers.
If we directly solve for rp , the above energy function can
be written as the sum of K terms, one for each spectral component and each term can be separately minimized. With a
little exercise, it can be shown that the minimal is achieved
exactly when rp = lp ; This is due to the inherent ambiguity
of the problem, when no other constraints are imposed on
rp . we reduce the ambiguity by exploiting the fact that the
reflectance spectra of typical object surfaces live in a low

dimensional subspace of RK , so that any reflectance vector can be written as a linear combination of Jr basis, with
Jr < K.
Let Br represent the K × J basis matrix for representing
the reflectance vector, rp can be written as rp = Br ∗ rep .
The energy function in Eq.(2) now becomes:
P
Esc =
kwp,q (Lp Br req − Lq Br rep )kd
p,q∈Nsc

(3)
d

Erc = kvp,q (Br rep − Br req )k

The combined energy can be represented in a matrix
form as:
e d + λ1 kVB Rk
e d
E = kWL,Br Rk
r

(4)

p,q∈Nsc

e consists of rep for all pixels in a vector. The
where R
matrix WL,Br depends on the neighborhood Nsc considered, the weight wp,q , the reflectance basis Br used, and
importantly the luminance data lp ; whereas the matrix VBr
depends on the neighborhood Nrc considered, the weight
vp,q and the reflectance basis Br used. Therefore, with this
non-logarithmic formulation, we encode the constraint due
to the measured luminance data in the matrix WL,Br .
The ambiguity with the scaling factor is inherent in all intrinsic image decomposition problems since only the product of reflectance and shading is known. To circumvent the
ambiguity about the scaling factor, we also explored another
solution, where we express the generic constraint on the coe = C, and augment the original energy
efficient sum as M R
function to enforce this constraint:
e d + λ2 kMr R
e − Ckd
e d + λ1 kVB Rk
Erefl = kWL,Br Rk
r
(5)
Similarly, the low rank nature of the shading space is
also widely acknowledged and exploited in prior work [14].
Shading is inherently low rank, because there are usually
only a few lighting sources with different illumination spectra acting in each captured scene, and the shading effect due
to geometry and shadowing only modifies the spectra by a
location-dependent scalar. If there is a single illumination
source and its spectrum is known or is able to be identified
by method of [36], we will use this spectrum(after normalization) as the only shading basis vector(Js = 1 and Bs
equals to this normalized spectrum). Likewise the problem
can be formulate as minimize
e d + λ1 kVL,Bs Sk
e d + λ2 kMs Se − Ckd
Eshad = kWBs Sk

shading and reflectance are known, represented by basis
matrices Bs and Br , respectively, so that sp = Bs sep and
rp = Br rep . We will use Se to denote the long vector cone
sisting of shading coefficient vectors sep at all pixels, and R
the long vector consisting of reflectance coefficient vectors
e
rep . We propose to solve sep and rep , or equivalently Se and R,
by minimizing a weighted average of the following energy
terms.
When shading is expected to be similar in pixels p and
q, we have sp ≈ sq and lp . ∗ rq ≈ lq . ∗ rp , or Lp rq =
Lq rp , where Lp is a diagonal matrix consisting of spectral
elements in lp . We formulate the energy functions directly:
X
kwp,q (Lp rq − Lq rp )kd + kwp,q (sp − sq )kd
Esc =

(6)

3.2. Simultaneous Recovery of Reflectance and
Shading
Based on the formulation that solves reflectance and
shading respectively, we proposed an optimization algorithm that simultaneously solves both shading and reflectance. We assume that the low rank subspace of the

e d + kWB Sk
e d
= kWL,Br Rk
s
(7)
When reflectance is expected to be similar in pixels p
and q, we have rp ≈ rq and lp . ∗ sq ≈ lq . ∗ sp , leading to a
regularization energy
X
Erc =
kvp,q (Lp sq − Lq sp )kd + kvp,q (rp − rq )kd
p,q∈Nrc

e d + kVB Rk
e d
= kVL,Bs Sk
r
(8)
The inherent data constraint lp = sp .∗rp leads to another
energy function:
X
d
e
Edata =
ksp . ∗ rp − lp kd = kQS,B
e s ,Br R − Lk
p
(9)
d
e
= kQR,B
e r ,Bs S − Lk
where QS,B
e s ,Br is a block diagonal matrix that dee and the basis matrices Bs and
pends on the solution for R
Br (likewise QR,B
).
e r ,Bs
e that minimizes a
The problem is to find Se and R
weighted average of the three energy functions:
E = Esc + λ1 Erc + 2λdata Edata


e d + kWBs Sk
e d + λ1 kVL,Bs Sk
e d + kVBr Rk
e d
= kWL,Br Rk
d
d
e
e
+ λdata kQS,B
e s ,Br R − Lk + λdata kQR,B
e r ,Bs S − Lk

(10)

Direct solution of the above problem solving sp and rp
simultaneously is hard because of the bilinear nature of the
e
data term. We applied iterative solution, where we solve R
and Se iteratively using alternating projection. As the dimension of the shading subspace is likely to be smaller than the
dimension of the reflectance subspace, we solve the shading
Se first. Also there are typically more subregions in an image
with similar reflectance, where it is easier to use the constant reflectance constraint to resolve the ambiguity about
shading.

The only problem is that we need to give a initial estimae in order to effectuate the data constraint in
tion of Se and R
Eq.(9). With the generic constraint in Part 3.1, Eq.(5) and
Eq.(6), we give a previous estimation of Se first, and then
e and find reflectance basis matrix, finally R
e
we can update R
e
and S are estimated respectively using the data constraint in
Eq.(9).
More specifically, the whole recovery algorithm can be
classified as Algorithm 1:

4.2. Low Rank Constraint

Algorithm 1: LRIID algorithm
1

2

3

4
5

6

7
8

Step 1: Assign constant-shading weights wp,q and
constant-reflectance weights vp,q .
Step 2: Solve an initial low rank estimate of the
shading image Se using a generic constraint.
Step 3: Solve an initial low rank reflectance estimate
e using a generic constraint and the data constraint
R
defined by the previous shading estimate.
repeat
Step 4: Solve Se using the data constraint defined
by the previous reflectance estimate.
e using the data constraint defined
Step 5: Solve R
by the previous shading estimate.
e converge;
until until the solution for Se and R
Step 6: Reconstruct S and R to get the refined shading
and reflectance.

4. Details

An important step to formulate a low rank constraint is
to derive the low rank basis. With the help of multispectral
imaging systems as PMIS [9] and CASSI [1], we can successfully get grouth-truth illumination spectra. Also, there
are plenty of work referring to how to extract illumination
from the image. For example, [36] can be applied in multispectral domain and performs well in implementation. In
order not to complicate our method, We simplify the process of getting normalized illumination spectra Bs by using
the grouth-truth illumination data.
When it comes to the low rank basis of reflectance, the
authors of [25, 28] have found Jr to be around 8 so as to
reach the best trade-off between expression power and noise
resistance in the process of fitting reflectance spectra.we set
Jr to be 8 and use the matrix introduced by [25] and perform Principle Component Analysis (PCA) to derive Br
from it.

4.3. Initial Estimation

4.1. Weight Choice
As we know, images suffering from poor light condition
may contain shadow area, which would in turn bring in unnecessary edges that confuse the algorithm. Various methods are used to determined weights wp,q and vp,q , including
pixel gradient [14, 21, 23], hue [35], correlation between
vectors [19] and learning [32]. we proposed a illuminationrobust and compute-friendly distance – normalized cosine
distance, to signify the differences between spectra of pixels in one neighborhood. This distance can be formulated
as
lp0 lq
(11)
|lp | · |lq |
dp,q tends to approximate to 0 when pixel p and q have
same spectra, and depart from 0 when spectra of p and q
are different. In order to derive weight wp,q , we need to
further magnify the difference between homogeneous and
heterogeneous pixels and make it more robust to the noise.
In our implementation
dp,q∈Nsc = 1 −

wp,q =

α and β are parameters of sigmoid function. To set α and
β, We sample 20 value of α within [1000, 10000] and 50
values of β within [10−5 , 10−2 ] and choose which perform
best.
Fig. 2 shows different separation results when β changes.
Local mean squared error (LMSE) is valued in 3 results. If
β is too small, shading tends to be more blurred; while β is
too big, reflectance would be blurred.

1
1+eα(dp,q −β)

vp,q = 1 − wp,q

(12)

In our implementation, only horizonally adjacent and
vertically adjacent pixels are considered in the neighborhood Nsc and Nrc . Suppose there are N pixels in the image.
e are 4N K × N Jr ,
The size of matrices WL,Br , VBr and R
4N K × N Jr and N Jr × 1 respectively in Eq.(3). Similarly, The size of matrices WBs , VL,Bs and Se in Eq.(6) are
4N K × N Js , 4N K × N Js and N Js × 1 respectively.
To avoid ambiguity, we further require that the shading
image has small deviation from the input image. In Eq.(5),
let Mr be a identity matrix, and C be a long vector concatenating all the pixels of the original image.
We use L2 form for all terms. so that the above problem for solving Eq.(5)or Eq.(6) is a quadratic programming
problem, and can be solved efficiently using conjugated gradient method. In Algorithm 1, The solution to the unconstrained optimization problem in Step 2 satisfies the following linear equation


T
T
e = Qs S
e
λ2 MsT C = WB
WBs + λ1 VL,B
V
+ λ2 MsT Ms S
s
s L,Bs
(13)

In Step 3, a data constraint defined by shading estimate
need to be added to Eq.(5), so that the linear equation can

(a) LMSE = 0.050, α = 1000, β = 0.01

(b) LMSE = 0.026, α = 5000, β = 0.0032

(c) LMSE = 0.031, α = 8000, β = 7.9e − 4

Figure 2. Results of threshold α and β. (b) achieves good result while shading and reflectance overlap clearly in (a) and (c).

5. Experimental Results

be written as
T
λdata QT
e L + λ2 M r C =
S


T
e = Qr R
e
Q
R
WL,Br + λ1 VBTr VBr + λ2 MrT Mr + λdata QT
WL,B
e
e
S
r
S
(14)

Because the matrix Qr and Qs is self-adjoint and sparse,
we can solve this equation iteratively, which typically converges very fast. Here, λ1 and λ2 are positive weights for
combining three different objective functions. In our implementation, we set λ1 = 2, λ2 = 0.01 and λdata = 1
empirically.

4.4. Iteration Performance
We use alternating projection to get refined shading and
reflectance. Just like what we stated in Step 4 and Step 5
in Algorithm 1, we we update shading first and then the
reflectance. In each round of iteration, gradient descent
method is applied. Eq. 10 converges whenever it reaches
the maximum iteration times 1000 or ∇E < 0.01.
Fig. 3 demonstrates the iteration performance of our algorithm. Both cost function and LMSE decrease via iteration. Before the iteration, some shadow which should be the
component of the shading image remains on the reflectance,
while the brightness of the reflectance image tends to be
uniform after iteration.

After
Iteration

Before
Iteration

Figure 3. Iteration Performance. The algorithm converges within
50 times. Two images demonstrate the reflectance image which
generated before and after iteration.

In this section, we provide extensive experimental validation of the proposed method. For the better visualization,
we show the result in pseudo-rgb and linearly normalized
the image to the range [0, 1]. We first show the performance
of our algorithm. This is followed by our method on on-line
dataset and visual results. Finally, we test on our dataset
with ground truth and compare our method with [12].

5.1. Experiments on proposed dataset
We provide a benchmark dataset with ground truth for
the performance evaluation of multispectral image intrinsic
decomposition problem. Following [15], we use local mean
squared error (LMSE) from the ground truth to measure
shading and reflectance image. We also compare with [12],
which proposed intrinsic image decomposition algorithm in
multispectral domain.
A benchmark dataset with ground-truth illumination,
shading, reflectance and specularity was presented in [12]
for performance evaluation of multispectral image intrinsic decomposition. Inspired by their ideas, we build up the
newest multispectral intrinsic ground-truth dataset including 12 scenes under the same environmental conditions.
We apply the updated mobile multispectral imaging camera to acquire the multispectral scenes, which could provide
higher resolutions in spectral data ranging from 450nm700nm with 118 spectral channels. Compared with the
dataset provided by [12], ours is a similar but flowering in
the diversity of projects with more detailed and bumpy scenarios which enables the dataset show further and potential
applications in other vision researches (e.g. IRSS, segmentation and recognition).
Here, we evaluate our algorithm via our proposed dataset
and use LMSE from the ground truth to validate our algorithm quantitatively. Compared with ground-truth, decomposition results that we achieved are desirable in terms of
both the LMSE score (0.018 in average for the entire data
set and the visual quality of the decomposed reflectance and
shading results).
We display 4 examples from our dataset. For all the input diffuse images equipped with multispectral data with-

out the effect from illumination, the corresponding visualized RGB images for reflectance and shading are listed
in Fig. 4 together with the LMSE results. It is clear that
our method could produce pleasing visually decompositions
based on multispectral images. In [12] generic constraint
forces the sum value of a group of pixels to be a certain
constant C, while in our method we let each pixel in shading image approximates to the pixel value of original image
in the same position. Our generic constraint matrix is more
sparse than that in [12], so that our algorithm will converge
more quickly. What we have to emphasize here is that,
we did all the image processing and metric computations
of LMSE on down-sampled 30 out of 118 spectral channels
(because [12] failed to process more bands on our 32G bytes
memory) but all multispectral images are finally visualized
above using corresponding synthesized RGB data.
Input Images

Reflectance

Table 1. Performance statistics for dataset image

box
cup
car1
bottle1
bottle2
bottle3
bus
car2
dinosaur
minion
plane
train
Avg.

Time(s)
SIID
LRIID
9.221
1.351
7.610
0.675
14.771 0.876
11.275 0.666
8.342
2.709
7.993
2.337
13.667 2.669
7.643
1.807
13.465 2.306
17.794 2.517
10.763 2.125
9.108
0.832
10.971 1.739

LMSE
SIID LRIID
0.032 0.023
0.016 0.012
0.025 0.012
0.062 0.031
0.005 0.008
0.009 0.009
0.030 0.031
0.030 0.024
0.021 0.023
0.020 0.018
0.024 0.015
0.017 0.015
0.024 0.018

Shading

BOX

CUP

PLASTER

tively more efficient than that of SIID, and the running time
approximates to the time in RGB cases in [29]. Moreover,
We note that our method is memory-friendly and is able to
process larger images with more spectral bands.
As to the reflectance, since it is not visualized from the
spectral perspective, we try to compare the spectral curves
of the reflectance from the ground truth and our algorithm.
we choose patches in some scenes of our ground truth. In
Fig. 5, it is obvious that our reflectance matches well with
the ground truth which means we could gain accurate spectral reflectance with better performance in computation.
0.55

CAR

0.5

Pixel Value

0.45

Figure 4. Results on sample images from the benchmark data set.
we only show reflectance and shading color images synthesized
from spectral data of 4 examples.

0.4
0.35
0.3
0.25
0.2
0.15
440

460

480

500

520

540

560

580

600

620

Wavelength(nm)

(a) GT

(b) Ours

0.24
0.22
0.2

Pixel Value

We display the computation time for examples in Fig. 4.
Here SIID denotes spectral intrinsic image decomposition
from the latest method based on Retinex in [12]. LRIID
denotes the low rank multispectral image intrinsic decomposition from our algorithm. We compare the computation
results obtained with and without low rank constraint. It
can be seen that the low rank constraint helps to improve
the computation time and decomposition results at the same
time.
In table 1, we demonstrate the performance statics for
dataset images. The last row shows the average metrics for
the ground truth dataset. On average, the SIID method requires almost 11 seconds to process a multispectral image
and generates results with 0.024 LMSE. While our method
LRIID takes 1.739 seconds and produce results with 0.018
LMSE. Although we reduce the average LMSE slightly
from 0.024 to 0.018, but the computation time is qualita-

(c) Spectra Curve

0.18
0.16
0.14
0.12
0.1
0.08
450

500

550

600

650

700

Wavelength(nm)

(d) GT

(e) Ours

(f) Spectra Curve

Figure 5. (a) and (d) are ground truth of our dataset, (b) and (e)
are reflectance images of our results, and the (c) and (f) are spectra
curve of marked area(solid red: curve of ours; dotted black: curve
of ground-truth).

5.2. Experiments on Nayar multispectral image
database[34]
A great variety of problems assume low-dimensional
subspace structure and have been solved by adding low-

(a)

(b)

(c)

(d)

(e)

Figure 6. Comparison of decomposition without and with the low-rank constraint. (a) A single image for five different scenarios in Nayar
dataset [34]. (b)-(c) and (d)-(e) show the reflectance and shading components computed by our LRIID solution without and with low-rank
constraint, respectively. Note that both the original multispectral images (a) and multispectral reflectance (d) are integrated into 3-channel
images by using the response curves of RGB sensors for visualization.

rank constraint, so is our method. To demonstrate the benefitd of this constraint, we compare results with and without
it in Fig. 6 using Nayar Multispectral Image Database [34].
The original images are shown in Fig. 6(a), the decomposed
reflectance and shading without and with the low-rank constraint are shown in Fig. 6(b)-(e) respectively. From the
comparison, we can see that the low-rank constraint helps
to maintain global structures and improves decompostion
results.

6. Conclusion
We have addressed the problem of the recovery of reflectance and shading from a single multspectral image captured under general spectral illumination. We have applied
a low rank constraint to settle the multispectral image intrinsic decomposition problem, which significantly reduced
the ambiguity. gradient descent has been used to give the
initial estimation of reflectance and shading, and alternating projection method has been applied to solve the bilinear
problem. Experiments on our dataset have demonstrated
that the performance of our work are better than prior work
in multispectral domain.

Our work has left out depth information. In fact, Retinex
theory fails to take effect when both shading and reflectance
change extensively in local area. Shading depends on the
object surface geometry, which can be derived from depth
information. We hope that we can make more accurate hypothesis about shading variation using the depth and surface
normal information in the future.

References
[1] G. R. Arce, D. J. Brady, L. Carin, H. Arguello, and D. S.
Kittle. Compressive coded aperture spectral imaging: An
introduction. IEEE Signal Processing Magazine, 31(1):105–
115, 2014. 5
[2] J. T. Barron and J. Malik. Intrinsic scene properties from a
single rgb-d image. In Proceedings of the IEEE Conference
on Computer Vision and Pattern Recognition, pages 17–24,
2013. 2
[3] J. T. Barron and J. Malik. Shape, illumination, and reflectance from shading. IEEE Transactions on Pattern Analysis and Machine Intelligence, 37(8):1670–1687, 2015. 2
[4] H. Barrow and J. Tenenbaum. Computer vision systems.
Computer vision systems, 2, 1978. 2
[5] M. Bell and E. Freeman. Learning local evidence for shading and reflectance. In IEEE International Conference on
Computer Vision, volume 1, pages 670–677. IEEE, 2001. 2
[6] S. Bell, K. Bala, and N. Snavely. Intrinsic images in the wild.
ACM Transactions on Graphics, 33(4):159, 2014. 2
[7] S. Bi, X. Han, and Y. Yu. An l1 image transform for edgepreserving smoothing and scene-level intrinsic decomposition. ACM Transactions on Graphics, 34(4):78, 2015. 2
[8] A. Bousseau, S. Paris, and F. Durand. User-assisted intrinsic
images. In ACM Transactions on Graphics, volume 28, page
130. ACM, 2009. 2
[9] X. Cao, H. Du, X. Tong, Q. Dai, and S. Lin. A prismmask system for multispectral video acquisition. IEEE
Transactions on Pattern Analysis and Machine Intelligence,
33(12):2423–2435, 2011. 5
[10] A. Chakrabarti and T. Zickler. Statistics of Real-World Hyperspectral Images. In Proceedings of the IEEE Conference
on Computer Vision and Pattern Recognition, pages 193–
200, 2011. 2
[11] Q. Chen and V. Koltun. A simple model for intrinsic image
decomposition with depth cues. In Proceedings of the IEEE
International Conference on Computer Vision, pages 241–
248, 2013. 2
[12] X. Chen, W. Zhu, Y. Zhao, Y. Yu, Y. Zhou, T. Yue, S. Du,
and X. Cao. Intrinsic decomposition from a single spectral
image. Applied Optics, 56(20):5676–5684, 2017. 1, 2, 6, 7
[13] G. D. Finlayson, M. S. Drew, and C. Lu. Intrinsic images by
entropy minimization. In European Conference on Computer
Vision, pages 582–595. Springer, 2004. 2
[14] B. V. Funt, M. S. Drew, and M. Brockington. Recovering shading from color images. In European Conference on
Computer Vision, pages 124–132. Springer, 1992. 5
[15] R. Grosse, M. K. Johnson, E. H. Adelson, and W. T. Freeman. Ground truth dataset and baseline evaluations for in-

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]
[24]

[25]

[26]

[27]

[28]

[29]

[30]

trinsic image algorithms. In IEEE International Conference
on Computer Vision, pages 2335–2342. IEEE, 2009. 2, 6
J. Ho, B. V. Funt, and M. S. Drew. Separating a color signal
into illumination and surface reflectance components: Theory and applications. IEEE Transactions on Pattern Analysis
and Machine Intelligence, 12(10):966–977, 1990. 2
C. P. Huynh and A. Robles-Kelly. A solution of the dichromatic model for multispectral photometric invariance. International Journal of Computer Vision, 90(1):1–27, 2010. 2
A. Ikari, R. Kawakami, R. T. Tan, and K. Ikeuchi. Separating
illumination and surface spectral from multiple color signals.
Digitally Archiving Cultural Objects, pages 297–321, 2008.
2
X. Jiang, A. J. Schofield, and J. L. Wyatt. Correlationbased intrinsic image extraction from a single image. In
European Conference on Computer Vision, pages 58–71.
Springer, 2010. 5
X. Kang, S. Li, L. Fang, and J. A. Benediktsson. Intrinsic
image decomposition for feature extraction of hyperspectral
images. IEEE Transactions on Geoscience and Remote Sensing, 53(4):2241–2253, 2015. 2
R. Kimmel, M. Elad, D. Shaked, R. Keshet, and I. Sobel. A
variational framework for retinex. International Journal of
Computer Vision, 52(1):7–23, 2003. 5
P.-Y. Laffont, A. Bousseau, and G. Drettakis. Rich intrinsic
image decomposition of outdoor scenes from multiple views.
IEEE Transactions on Visualization and Computer Graphics,
19(2):210–224, 2013. 2
E. H. Land and J. J. McCann. Lightness and retinex theory.
Josa, 61(1):1–11, 1971. 2, 5
K. J. Lee, Q. Zhao, X. Tong, M. Gong, S. Izadi, S. U. Lee,
P. Tan, and S. Lin. Estimation of intrinsic image sequences
from image+ depth video. In European Conference on Computer Vision, pages 327–340. Springer, 2012. 2
L. T. Maloney. Evaluation of linear models of surface spectral reflectance with small numbers of parameters. JOSA A,
3(10):1673–1683, 1986. 2, 5
A. Nadian-Ghomsheh, Y. Hassanian, and K. Navi. Intrinsic image decomposition via structure-preserving image smoothing and material recognition.
PloS one,
11(12):e0166772, 2016. 2
R. A. Newcombe, S. Izadi, O. Hilliges, D. Molyneaux,
D. Kim, A. J. Davison, P. Kohi, J. Shotton, S. Hodges, and
A. Fitzgibbon. Kinectfusion: Real-time dense surface mapping and tracking. In IEEE International Symposium on
Mixed and Augmented Reality, pages 127–136. IEEE, 2011.
2
J. P. Parkkinen, J. Hallikainen, and T. Jaaskelainen. Characteristic spectra of munsell colors. JOSA A, 6(2):318–322,
1989. 2, 5
L. Shen, P. Tan, and S. Lin. Intrinsic image decomposition with non-local texture cues. In IEEE Conference on
Computer Vision and Pattern Recognition, pages 1–7. IEEE,
2008. 1, 2, 7
L. Shen, C. Yeo, and B.-S. Hua. Intrinsic image decomposition using a sparse representation of reflectance. IEEE
Transactions on Pattern Analysis and Machine Intelligence,
35(12):2904–2915, 2013. 1, 2

[31] M. F. Tappen, E. H. Adelson, and W. T. Freeman. Estimating intrinsic component images using non-linear regression.
In IEEE Computer Society Conference on Computer Vision
and Pattern Recognition, volume 2, pages 1992–1999. IEEE,
2006. 2
[32] M. F. Tappen, W. T. Freeman, and A. E. H. Recovering intrinsic images from a single image. In IEEE Transactions
on Pattern Analysis and Machine Intelligence, volume 27,
pages 1459–1472, 2005. 2, 5
[33] Y. Weiss. Deriving intrinsic images from image sequences.
In IEEE International Conference on Computer Vision, volume 2, pages 68–75. IEEE, 2001. 2
[34] F. Yasuma, T. Mitsunaga, D. Iso, and S. K. Nayar. Generalized assorted pixel camera: postcapture control of resolution,
dynamic range, and spectrum. IEEE Transactions on Image
Processing, 19(9):2241–2253, 2010. 2, 7, 8
[35] Q. Zhao, P. Tan, Q. Dai, L. Shen, E. Wu, and S. Lin. A
closed-form solution to retinex with nonlocal texture constraints. IEEE Transactions on Pattern Analysis and Machine
Intelligence, 34(7):1437–1444, 2012. 5
[36] Y. Zheng, I. Sato, and Y. Sato. Illumination and reflectance
spectra separation of a hyperspectral image meets low-rank
matrix factorization. In Proceedings of the IEEE Conference
on Computer Vision and Pattern Recognition, pages 1779–
1787, 2015. 2, 4, 5
[37] T. Zhou, P. Krahenbuhl, and A. A. Efros. Learning datadriven reflectance priors for intrinsic image decomposition.
In Proceedings of the IEEE International Conference on
Computer Vision, pages 3469–3477, 2015. 2

