Shape Reconstruction Based on Similarity inRadiance Changes under Varying Illumination
Imari Sato
National Institute of Informatics

Takahiro Okabe Qiong Yu Yoichi SatoThe University of Tokyo

imarik@nii.ac.jp{takahiro,qiong,ysato}@iis.u-tokyo.ac.jp

Abstract
This paper presents a technique for determining an object’s shape based on the similarity of radiance changes observed at points on its surface under varying illumination. To
examine the similarity, we use an observation vector that represents a sequence of pixel intensities of a point on the surface
under different lighting conditions. Assuming convex objects
under distant illumination and orthographic projection, weshow that the similarity between two observation vectors is
closely related to the similarity between the surface normalsof the corresponding points. This enables us to estimate theobject’s surface normals solely from the similarity of radiance changes under unknown distant lighting by using dimensionality reduction. Unlike most previous shape reconstruction methods, our technique neither assumes particular reﬂection models nor requires reference materials. This makes
our method applicable to a wide variety of objects made ofdifferent materials.

1. Introduction
The appearance of an object is determined by several factors including illumination, viewing geometry, and the shapeand surface material of the object. Changing any one of thesefactors should change the object’s appearance. These factors
are, in general, nonlinearly related to the object’s appearance,and consequently an inverse problem of estimating them becomes difﬁcult unless there is some a priori knowledge of thescene. Most of the previous techniques for solving this inverse problem thus try to estimate some of these factors fromimages of a scene by assuming that the others are known.
Previous studies of shape reconstruction have demonstrated that the shape of an object can be recovered from multiple images of the object taken under various lighting conditions assuming that there is some knowledge of the scene,e.g., illumination and/or surface materials [19, 20, 7, 12].1
1 Most of the shape-from-shading approaches estimate an object’s shapefrom a single image assuming distant illumination and uniform Lambertian
reﬂectance [21].

978-1-4244-1631-8/07/$25.00 ©2007 IEEE

Radiance changes

Varying illumination
Figure 1. Determining an object’s shape based on similarity of radiance changes observed at its surface points.

The classical photometric stereo approach presented in [20, 7]
recovers the shape of a Lambertian object from multiple images of the object taken under known light sources. Later,
this approach was extended to several directions including photometric stereo methods for non-Lambertian surfaces[8, 3, 12, 4] (see [19] for a good survey) and uncalibrated
photometric stereo methods [5, 1, 2].
However, real-world materials sometimes have complex
appearances that prevent us from extracting their shapes using analytic reﬂectance models. To cope with this problem,
Hertzmann and Seitz proposed the use of a calibration object[6], which had also been suggested in early work on photometric stereo [7]. Instead of computing a reﬂectance map
based on an analytic reﬂection model, images of a calibration object of a known shape made of the same material arecaptured under the same lighting conditions as a target object
and used as an empirical reﬂectance map. Later, Hertzmann
and Seitz extended their method so that instead of using asingle calibration object, it uses a small number of them toestimate the shape of objects with non-uniform surface materials. However, the issue of what kinds of reference materials are suitable for representing the complex appearancesof a wide variety of real-world materials still remains to beaddressed.

Recently, Koppal and Narasimhan presented a clusteringbased technique to ﬁnd iso-surface normal clusters of scene
points without requiring knowledge about materials and lighting [9]. Their approach shows how effective it is to analyze
the temporal variation in the appearance of a scene for obtaining its meaningful geometric structure; e.g., each face ofa planar object is found regardless of its texture. However,since the method ﬁnds iso-surface normal clusters as the geometric structure of the scene, the 3D shape of the scene cannot be obtained unless another shape reconstruction techniquesuch as uncalibrated photometric stereo is used, as suggested
in their work.
In this work, we present a method for recovering an object’s shape directly from input images of the object based on
the similarity of the radiance changes observed at points on
its surface under varying illumination (Figure 1). To examine
the similarity, we consider an observation vector that represents a sequence of pixel intensities of a surface point under
different lighting conditions. Assuming convex objects under distant illumination and orthographic projection, we showthat the similarity between two observation vectors is closely
related to the similarity between surface normals of the corresponding points.
This gives rise to the interesting theory that within thespace containing a set of observation vectors, all of the observation vectors should lie on a manifold that represents thesurface normals of the corresponding surface points. Accordingly, we can ﬁnd the surface normals of an object by using
a dimensionality reduction technique that best preserves theintrinsic structure in the input space of observation vectors
as illustrated in Figure 2. Our method is capable of estimating the dense distribution of an object’s surface normals by
analyzing the similarities of radiance changes between scene
points.
As far as we know, this is the ﬁrst attempt to discover an
object’s shape solely from a set of images of the object basedon manifold learning techniques. Unlike most previous shapereconstruction methods, our technique neither uses particular reﬂection models, nor requires empirical reﬂection maps.
This makes our method applicable to a wide variety of objects
with different types of material.
The rest of this paper is organized as follows. Section
2 considers the similarity of radiance changes observed atpoints on an object’s surface under varying illumination, andSection 3 outlines our idea of estimating an object’s surface
normals from the similarity of radiance changes. Section 4
discusses the issue of how the similarity of two observation
vectors is related to the similarity of their surface normals.
Section 5 describes our algorithm for recovering the object’s
shape by using a dimensionality reduction technique. Section
6 shows experimental results for synthetic and real data andSection 7 presents concluding remarks.

N-d space

3-d space
z

nq

||

Iw

Iq

||

nw

x
y
Preserve intrinsic structure

2-d manifold

2-d hemisphere

Figure 2. 3-d embedding of observation vectors.

Ir

Iq

Iw

Is

Euclidean
Geodesic
Figure 3. Geodesic distance measured along a manifold.

2. Similarity of Radiance Changes
Let us ﬁrst analyze the similarity of radiance changes observed at points on an object’s surface under n different lightsources distributed around the object to illuminate it from various directions. Suppose Ipk is the brightness of each pixel or
each corresponding surface point p seen under the kth illumination with unit radiance.
A vector containing the set of Ipk (k = 1, 2, · · · , n) in its
elements is normalized( 1)T
Ip , · · · , Ipn
Ip = √∑
(1)
n
k )2
(I
k=1 p
and referred to as an observation vector. From input images
with m surface points (p = 1, · · · , m), we obtain m observation vectors. In the following, a n-dimensional space containing a set of m observation vectors is referred to as an observation space.
Assuming distant light, orthographic projection, isotropic
reﬂection, and no local effects such as cast shadows and interreﬂections, we consider the brightness of p under the kth
illumination described in the form:
Ipk = ρp g(θk′ , θv′ , φ′k − φ′v ),(2)

where (θk′ , φ′k ) and (θv′ , φ′v ) are the incident and reﬂection directions deﬁned in the local coordinates with respect to thepoint p. Function g represents how much of the incident lightfrom the direction (θk′ , φ′k ) is reﬂected on the point toward(θv′ , φ′v ). ρp is a reﬂection coefﬁcient that varies over the surface.

Let np be the surface normal of p and lk and v denote theincident direction of the kth illumination and the direction to
the viewer represented by the global coordinate, respectively.
Note that lk and v are constant over all surface points of theobject under the assumptions of distant illumination and orthographic projection.
The angles θk′ , θv′ , and φ′k − φ′v can be computed from np ,v, and lk ,2 and (2) is thus rewritten as
Ipk = ρp g ′ (np , lk , v),(3)

where the function g ′ represents how much of the incident
light lk is reﬂected on the point of its normal np toward v.
From the normalization process in (1), reﬂectance coefﬁcient ρp comes to have no inﬂuence on the observation vector Ip . From this, we see that differences in the observation
vectors among all surface points are caused by differencesin their surface normals. Namely, Iq = Iw is equivalent to
nq = nw for surface points q and w, which was used as theorientation-consistency cue in [6].
As we will see later in Section 4, the similarity between
two observation vectors is closely related to the similarity
between the surface normals of the corresponding surface
points. This enables us to estimate the distribution of an object’s surface normals by analyzing the similarities of its observation vectors.
Note that the form in (2) can represent a wide variety of
BRDFs including commonly used reﬂection models such as
Lambertian, Oren-Nayar [13], Phong, Torrance-Sparrow[16],and isotropic Gauss[17] models. On the other hand, some
BRDFs do not follow this form exactly: e.g., a dichromatic
reﬂection model whose relative strength between its diffuse
and specular components varies over the object surface. However, as shown later in our experimental results, our method
is still able to recover object shapes reasonably well even for
BRDFs that do not follow (2) exactly.

3. Shape Recovery by Embedding of Observation Vectors
The set of m observation vectors seen under n different
light sources can be thought of as m points in n-dimensional
space. As we see in the previous section, differences in theobservation vectors among all surface points are caused bydifferences in their surface normals in the case where their nelements are described in the form (2).
It can be further said from this that within the observation space, all of the m points lie on an intrinsically twodimensional manifold that represents the surface normals of
the corresponding m surface points because surface normals
have two degrees of freedom.3
2 Since cos θ ′ = nT l , cos θ ′ = nT v, cos(φ′ − φ′ ) = [lT v −
p k
v
p
v
k
k
k
T
−1 (nT l ) sin cos−1 (nT v)].
(np lk )(nTp v)]/[sin cos
p kp
3 Note that this manifold has two dimensions because surface normals(nx , ny , nz ) are described in elevation and azimuthal angles (θ, φ).

The key idea of our work is that since surface normals aredeﬁned on a unit sphere in a three-dimensional space, if we
ﬁnd three-dimensional embedding of the observation vectors
using the dimensionality reduction that best preserves its intrinsic structure, the resulting three-dimensional space should
reveal the distribution of the object’s surface normals as illustrated in Figure 2.
It is important to note that what we count on here is thetwo-dimensional manifold hidden in the observation space
and the similarity of two observation vectors needs to be determined by their shortest distances measured along a manifold. Such a distance is called geodesic distances.
In our method, a nonlinear embedding method, called isometric feature mapping (Isomap) is used for dimensionality
reduction [15].4 Isomap assumes that, for neighboring data
points, distances deﬁned in an input space such as a Euclidean
distance provides good approximations of the geodesic distances. Accordingly geodesic distances for faraway points
are computed by adding up a sequence of distances between
neighboring points.
For instance, in Figure 3, the geodesic distance between
two observation vectors Iq and Is denoted to as dG (q, s) isapproximated by adding up a sequence of Euclidean distances
between neighboring points as
dG (q, s) ≅ |Iq − Iw | + |Iw − Ir | + |Ir − Is |.
Tenenbaum et al. demonstrated that Isomap could ﬁnd
low-dimensional representations of natural images, e.g., images with different poses and lighting directions [15]. Previous studies have also shown that scene conditions underwhich input images are captured could be recovered throughdimensionality reduction: Pless and Simon estimated viewing directions or object poses in [14], and Winnemoller et al.
recovered the light source directions in [18]. However, notheoretical insight into why dimensionality reduction is effective in extracting such information from the input image
space has been provided.
In the next section, we will show that the Euclidean distance between two neighboring observation vectors is indeed
a good approximation of the difference between their surface
normals. This means that Isomap is suitable for revealing the
intrinsic structure of the observation space.
It should be also noted that previously proposed methodsof analyzing images using dimensionality reduction aimed atrecovering scene conditions under which images of unknown
subjects are captured, e.g., viewing direction or light source
directions. On the other hand, our goal is to extract the information about a photographed subject from images captured
under unknown conditions.
4 The set of observation vectors generally has intrinsically nonlinear structures, and classical techniques for dimensional reduction, such as principal
component analysis (PCA) and multidimensional scaling (MDS) are not effective in revealing such structures.

4. Theoretical Insight into Similarity of Observation Vectors
This section discusses the important issue of how the similarity of two observation vectors is related to the similarity of
their surface normals. To get some insight into good similarity measures for ﬁnding a manifold representing surface normals, we derive a distance measure for two commonly used
reﬂectance models.
We show that the obtained measures are approximately
equal to or closely related to the Euclidean distance |Ii − Ij |when the distance between two observation vectors of two
surface points i and j is small (Ii ≅ Ij ). Note that, eventhough the two similarity measures are derived from parametric reﬂection models, our method is by no means restricted toparticular reﬂection models.
Recall that Isomap assumes that, for neighboring data
points, a distance deﬁned in an input space provides good
approximations of the geodesic distances. If two observation
vectors are detected as neighbors with a small distance, thisshould also indicate that the distances between their surface
normals are small.

4.1. Discussion of Body Reﬂection
First, we consider radiance changes of body reﬂection as
an example. As shown in Appendix A, assuming that a convex object with a Lambertian surface is illuminated by a setof directional light sources distributed around the object withuniform density, the inner product of observation vectors Ii
and Ij is given by
ITi Ij ≅ nTi nj ,(4)when the corresponding surface normals satisfy ni ≅ nj .
Because the observation vectors are normalized so thatthey have unit length, it is clear that
ITi Ij = 1 − (|Ii − Ij |2 /2).

(5)

In addition, when Ii ≅ Ij , the angle between ni and nj issmall. Then, the following approximationnTi nj = cos[arccos(nTi nj )] ≅ 1 −2

arccos (nTi nj )2

(6)holds. From (4), (5), and (6), the Euclidean distance between
observation vectors Ii and Ij is given by|Ii − Ij | ≅ arccos(nTi nj ),(7)when ni ≅ nj . This means that |Ii − Ij | is approximately
equal to the ideal distance that is the angle between ni and nj
that represents their relationship.

4.2. Discussion of Surface Reﬂection
Second, we investigate radiance changes of surface reﬂection as another example. As shown in Appendix B, the innerproduct of observation vectors Ii and Ij is given by[]arccos2 (nTi nj )TIi Ij ≅ exp −
2σ 2

(8)

for the simpliﬁed Torrance-Sparrow model, when ni ≅ nj
and surface roughness σ ≪ π.
Because ni ≅ nj is equivalent to Ii ≅ Ij , the following
approximation()|Ii − Ij |2|Ii − Ij |2
T
ln(Ii Ij ) = ln 1 −
≅−
(9)
2
2
holds. Hence, the Euclidean distance between observation
vectors Ii and Ij is given by|Ii − Ij | ≅arccos(nTi nj )
.
σ

(10)

In other words, the distance |Ii − Ij | is proportional to the angle between ni and nj . This means that the geodesic distanceis approximately equal to the ideal distance up to some scale
factor. Note that this scale factor can be removed by normalizing geodesic distances between all pairs of observation
vectors so that 0 ≤ dG (i, j) ≤ π.

5. Algorithm for Shape Recovery
Given m sets of n-dimensional observation vectors, thesurface normals and heights of the corresponding m surface
points are estimated as described in the following steps:
1. Provide a geodesic distance matrix that contains geodesic
distances between all pairs of observation vectors (Section
5.1).
2. Find the three-dimensional embedding of the observation
vectors based on the geodesic distance matrix (Section 5.2).
3. Recover surface normals and object heights. The occluding boundary is used as a reference to transform the output
from Isomap to the true distribution of the object’s surface
normals (Section 5.3).

5.1. Geodesic Distance Matrix
The geodesic distances between all pairs of observation
vectors are approximated by adding up a sequence of Euclidean distances between neighboring points as follows:
1. Compute Euclidean distances d(i, j) between all pairsof observation vectors Ii and Ij asd(i, j) = |Ii − Ij |.

(11)

2. Initialize a graph G where points i and j are connected if
i is one of the k nearest neighbors of j based on d(i, j).
Geodesic distances dG (i, j) for all i, j pairs of m observation vectors are initialized such that dG (i, j) = d(i, j)
if i and j are connected as neighbors and dG (i, j) = ∞otherwise.

3. Compute the shortest path distances between all pairsof observation vectors in G by using Floyd’s algorithm
and update the geodesic distance dG (i, j). This provides
a matrix of the geodesic distances between all pairs ofobservation vectors Dg = {dG (i, j)}.

5.2. Dimensionality Reduction Based on Geodesic
Distance
After the matrix of geodesic distances Dg is obtained,classical MDS is applied to this matrix to ﬁnd a threedimensional embedding of the observation vectors such thatthe estimated intrinsic geometry is best preserved through dimensionality reduction.
Let ep be a three-dimensional vector corresponding to observation vector Ip after dimensionality reduction. Then ep
are determined by minimizing the cost function
E = |τ (Dg ) − τ (D3 )|,(12)

where D3 is the matrix of Euclidean distances in the threedimensional√space d3 (i, j) = |ei − ej | and |A| is the Frobe∑
2
nius norm
ij Aij of matrix A. The τ operator converts
distances to inner products to support efﬁcient optimization.5
By setting the coordinates of ep to the top three eigenvectorsof the matrix τ (Dg ), we are able to ﬁnd the global minimumof the cost function (12) [11].

5.3. Solving Ambiguity of Object’s Surface Normals
There are distance-preserving transformations such astranslations, reﬂections, and rotations in a three-dimensional
space, and the estimated three-dimensional vector ep does notnecessarily correspond to the true distribution of the object’s
surface normals. In other words, what Isomap computes is
the relative relationships of these surface normals. In addition, ep does not necessarily has a unit length. For this reason,
we need a further step to solve this ambiguity and obtain thetrue distribution of the object’s surface normals. We project
ep onto a unit sphere to ensure |eb | = 1 and ﬁnd a valid transformation that projects ep onto the unit sphere based on the
known surface normals along the occluding boundary .
Under the orthographic projection viewed from the direction (0, 0, 1), the true surface normals of the occluding
boundary can be obtained as their gradient directions in the
2D input image coordinate systems. Let cb be the obtainedtrue surface normals of those boundary points and eb be theircorresponding embedded vectors. Then a transformation ς()that projects the embedded vectors eb onto
∑ a unit sphere isfound by minimizing the cost function all b |cb − ς(eb )|.
Here the center and orientation (θ = 0 and φ = 0 directions)of the unit sphere deﬁnes this transformation, and ς() are thus
estimated with respect to those parameters.
5 τ (D) = −HSH/2, where S is the matrix of squared distances S
ij =2 and H is the centering matrix H = δ − 1/m. See Tenenbaum [15]
Dij
ij
ij
and Mardia [11] for details.

After projecting all embedded vectors ep onto the unit
sphere based on the estimated ς(), the recovered surface normals are stored in their original pixel locations to provide anormal ﬁeld, as shown in Figure 6(b). Then the height ﬁeldof the object is estimated from the recovered normal ﬁeld byusing the shapelet technique [10].

6. Experimental Results
6.1. Synthetic Data
We evaluated the accuracy of our technique for shape reconstruction using synthetic data. The images of objects withdifferent surface materials were synthesized under 450 distant light sources that were randomly distributed around theobjects with (1) diffuse (textured Lambertian), (2) specular(Torrance-Sparrow), (3) diffuse (uniform Lambertian) + specular (Torrance-Sparrow), and (4) diffuse (textured Lambertian) + specular (Torrance-Sparrow) surfaces.
The residual variance after the dimensionality reduction isshown in Figures 4. We can clearly see that residual variances decrease as the dimensionality increases and the curve
levels off after three dimensions. This indicates that threedimensional embedding of the observation vectors were successfully found. It should be noted that in the process of dimensionality reduction, there is only one parameter k representing the number of neighbors, which was adjusted to minimize the cost function (12).
The estimated surface normals and shapes of the objects
and some of the input images used for this experiments areshown in Figure 5. In this ﬁgure, the estimated surface normals are stored directly in the RGB values of the surface normal map images where absolute values of x, y, and z coordinates are shown in red, green, and blue, respectively: x andy axes correspond to the horizontal and vertical axis of theimage plane, and z axis is deﬁned to be perpendicular to theimage plane.
Compared with the ground truth provided in Figure 5 (a),we clearly see that our technique was able to estimate thesurface normals and shapes of the object quite well for all
four cases in spite of the fact that they were recovered solely
from input images without requiring light source directions
or reﬂection models.
Table 1 provides a quantitative evaluation of the techniquethat summarizes the root-mean-square errors in the estimated
Case
(1)

Diffuse Texture Specular
⋄
⋄
⋄

(2)
(3)
(4)

⋄
⋄

⋄

⋄
⋄

Distance
d
dd
d
ds
d
d

Normal
5.7
5.7
8.5
6.9
5.2
6.0

Height
0.066
0.066
0.093
0.052
0.059
0.072

Table 1. Root-mean-square errors in the estimated surface normals
and object’s heights.

(a)(a) synthetic data(b) real data

Figure 4. Residual variance after the dimensionality reduction.

(b)(a) ground truth





(c)



(b) some input images(a)(c) case 1(d) case 2(b)(e) case 3(f) case 4(a)Figure 5. Synthetic data: estimated surface normals and object
heights.

surface normal directions and heights for all cases. Here,the error in surface normal directions are shown in degrees,
and the heights are normalized between 0 and 1. For cases(1) and (2), we also tested the ideal distances analyticallyderived in the previous section. The distance between two
neighboring observation vectors d(i, j) = |Ii − Ij | deﬁned
in (11) was replaced with the ideal distances derived fromT(4) and (8):√ dd (i, j) = arccos(Ii Ij ) for Lambertian, andTds (i, j) = − ln(Ii Ij ) for the simpliﬁed Torrance-Sparrowmodel. The column Distance in this table indicates which
distance was used for the estimation.
In Table 1, for all cases, our method could estimate surface
normals and object’s heights with reasonably high accuracy.
This shows that our algorithm successfully found the intrin-

(b)

(c)

(c)

Figure 6. Results for real images: (a) input images, (b) surface normals, and (c) object’s heights.

sic structure of the observation space to be used for shapereconstruction. Especially, in the cases of (1) and (2), about
the same accuracy was achieved for the surface normals using
the distance d and the ideal distances dd and ds to measure thesimilarity between two neighboring observation vectors. This
indicates that the distance d well approximates the geodesic
distance between these vectors.

On the other hand, strictly speaking, case (4) does not satisfy the necessary conditions of our approach since its radiance changes cannot be described in the form (2). Nevertheless, we could achieve about the same accuracy as in case (3).
Body reﬂection is often observed at most points on an object’s
surface under varying illumination, while specular reﬂection
are less likely to be observed. We consider that the diffuse
reﬂection components contributed more to the shape reconstruction than did the specular reﬂection components in thiscase.

Appendix A
Assuming the Lambertian model, let us consider radiance changes observed
at a point i on an object’s surface under a set of directional light sources{l1 , l2 , · · ·, lF }. Here, the light source vector lk is a product of its brightness and a unit vector representing the direction
q of the light source, so theunnormalized observation vector I′i (Ii = I′i /T
T
I′i = (max(ρi nT
i l1 , 0), · · ·, max(ρi ni lF , 0)) ,I′i I′jT

7. Conclusion
We presented a novel technique for recovering an object’s
shape directly from input images of the object that is basedon the similarity of the radiance changes observed at points
on its surface under varying illumination. The main contribution of our work is that we have shown that the similarity of the radiance changes is closely related to the similarity
between surface normals of the corresponding points and introduced a new theory that within the observation space, allof the observation vectors of these points lie on a manifoldthat represent their surface normals. This enables us to ﬁnd
surface normals of an object by using a dimensionality reduction technique that best preserves the intrinsic structurein the observation space as captured in geodesic distances ona manifold. A future research direction of this work would
be to analyze the manifold in more detail by using objects ofknown shapes with various kinds of surface materials and toinvestigate how many images are generally sufﬁcient in orderto obtain an object’s shape correctly.

(13)

where ρi and ni are the albedo and normal, respectively. Note that we represent attached shadows using operator max(∗, 0) but do not take account ofcast shadows.
The inner product of vectors I′i and I′j is derived as

6.2. Real Objects
We also tested our technique using real images of several
objects with different surface materials. Figure 6(a) shows theobjects used in this experiment: cat (painted bisque), lemon(plastic), and sheep (ceramic) ornaments. Figure 6 (b) and (c)show the estimated surface normals and the objects’ heights.
For each object, 130 ∼ 150 images were captured by moving
a point light source around the object, which was roughly 1
m away from the object of 5 ∼ 10 cm in diameter.
In this ﬁgure, we see that our method was able to estimate
their shapes regardless of the complexity of their surface materials and shapes. In the case of the cat ornament, this object
has both convex and concave shapes and thus some cast shadows were observed in the images. Strictly speaking, this does
not follow our assumption of a convex object. However, it
was found through the experiment that cast shadows did notcause that much trouble in estimating an object’s surface normals as long as its unshadowed appearances were observed
enough in other input images under different lighting conditions. It seems that our method is less sensitive to this kindof obstruction in cases where observation vectors contain correct values for most of their elements.

|I′i |2 ) is given by

=

F
X
k=1

=

T
max(ρi nT
i lk , 0) max(ρj nj lk , 0)

0

1

X

B
ρi ρj nT
i @

T
{k|nT
i lk >0,nj lk >0}

where

(14)

C
lk lT
k A nj ,(15)

Pmeans the summation over light source vecT
{k|nT
i lk >0,nj lk >0}Ttors satisfying both nT
i lk > 0 and nj lk > 0. That is, the summation istaken when neither i nor j are in attached shadows.
When ni ≅ nj , the summation is approximated as
0
1
X
B
C
T
I′i I′j ≅ ρi ρj nT
lk lT
(16)
i @
k A nj{k|(ni +nj )T lk >0}

0

=

2
lx
B
ρi ρj nTl
i @ x ly
lx lz

lx ly
ly2
ly lz

1
lx lz
C
ly lz A nj ,
lz2

(17)

because both i and j are in attached shadows for almost the same lightP 2
2 =sources. Here, we denote lk = (lkx , lky , lkz )T and deﬁne lx
k lkx .
Assuming that directional light sources with brightness l distributearound the object with uniform density, the above summations can be converted to integrals such as
Z π/2 Z 2π2π 22 = l2(sin θ cos φ)2 sin θdθdφ =
l .
(18)
lx
3
0
0
By computing the other elements of the matrix in (17) similarly, we obtain
0
1
1 0 0
2π 2
′T ′
T @
0 1 0 A nj
Ii Ij ≅l ρi ρj ni
(19)
3
0 0 1
The norm of I′i is given by
Z π/2 Z 2π2π 2 2|I′i |2 = ρ2i l2(cos θ)2 sin θdθdφ =l ρi .
30
0(20)Therefore, the inner products of the normalized observation vectors are
related to the inner products of the surface normals at the corresponding
points as
I′i T I′j
q
IT
I
=
= nT
(21)
j
i
i nj ,|I′i |2 |I′j |2qbecause |I′i |2 |I′j |2 = 2πl2 ρi ρj .
3

Appendix B
We derive the relationship between the similarity of observation vectors
and the similarity of surface normals for the simpliﬁed Torrance-Sparrowmodel. Let us consider a polar coordinate system (Figure 7) whose northpole coincides with the viewing direction v, and denote surface normals ni ,

Therefore, the inner products of the normalized observation vectors are
related to the angles between surface normals as













„
«
I′i T I′j
β2
q
I
=
IT
≅
exp
−
.
i j
2σ 2
|I′i |2 |I′j |2






(29)









References

Figure 7. Polar coordinate system.

nj , and a light source vector l by (θi , φi ), (θj , φj ), and (θ, φ). Then, half
vector h, which bisects v and l, is given by (θ/2, φ). We denote the angles
between h and ni , h and nj , and ni and nj by αi , αj , and β.
For the sake of simplicity, we do not take the Fresnel and geometric attenuation terms into consideration. We also obtain radiance Ii at a point iilluminated by a single directional light source l with brightness l as
„
«
α2
ρi l
Ii = T exp − 2i .
(22)
σ
ni v
Here, ρi is the reﬂectance of surface reﬂection, and σ represents surface
roughness. The angle αi satisﬁesθθsin θi cos(φ − φi ) + cos cos θi .
(23)
2
2
Assuming that the directional light sources with brightness l are distributed around the object with uniform density, the inner product of observation vectors I′i and I′j is represented by an integral with respect to l
as
Z Z
ρi l ρj l
T
I′i I′j =
T
T
nT
{l|nT
i v nj v
i l>0,nj l>0,|l|=1}!
22
αi + αj
dl
(24)
exp −
σ2
cos αi = hT ni = sinin a similar manner to the Lambertian model.
As is often the case with real materials, we assume that the specular reﬂection has a sharp peak and then σ ≪ π. Thus, only combinations of ni ,nj , and l that satisﬁes α2i + α2j ≪ π 2 can contribute to the integral in (24).
Therefore, we assume ni ≅ nj , and consider l such that h ≅ ni ≅ nj .
Accordingly, taking Taylor series expansion of (23) and ignoring higher
order terms with respect to αi , (θ/2 − θi ), and (φ − φi ), we obtain
«2
„
θ
θ
α2i ≅
− θi
+ sin sin θi (φ − φi )2 .
(25)
2
2
Similarly, we can derive
β 2 ≅ (θi − θj )2 + sin θi sin θj (φi − φj )2 .

(26)

Substituting these approximations into (24), we can derive
„
«Z Z
β2
ρi l ρj l
T
exp−
I′i I′j ≅T
2σ 2nT{l|(ni +nj )T l>0,|l|=1}
i v nj v»
11(θ − θi − θj )2
exp − 2σ2„« )#φi + φj 2
θ+ sin (sin θi + sin θj ) φ −dl. (27)
2
2
The norm of the observation vector is calculated as
Z Z(ρi l)2|I′i |2 ≅T
2(ni v){l|nT
i l>0,|l|=1}

ﬀ–
»
1
θ
1
(θ−2θi )2 +2 sin sin θi (φ−φi )2
dl. (28)
exp − 2σ22[1] P. N. Belhumeur, D. J. Kriegman, and A. L. Yuille. The bas-relief
ambiguity. Int’l J. Computer Vision, 35(1):33–44, 1999. 1[2] M. K. Chandraker, F. Kahl, and D. J. Kriegman. Reﬂections on thegeneralized bas-relief ambiguity. In Proc. IEEE Conf. Computer Vision
and Pattern Recognition, pages I–788–795, 2005. 1[3] E. N. Coleman and R. Jain. Obtaining 3-dimensional hape of textured
and specular surfaces using four-source photometry. Computer Vision,
Graphics, and Image Processing, 18(4):309–328, 1982. 1[4] A. S. Georghiades. Incorporating the torrance and sparrow model ofreﬂectance in uncalibrated photometric stereo. In Proc. Int’l. Conf.
Computer Vision, pages 816–825, 2003. 1[5] H. Hayakawa. Photometric stereo under a light-source with arbitrary
motion. J. Optical Soc. Am. A, 11(11):3079–3089, 1994. 1[6] A. Hertzmann and S. M. Seitz. Example-based photometric stereo:
Shape reconstruction with general, varying brdfs. IEEE Trans. Pattern
Analalysis and Machine Intelligence, 27(8):1254–1264, 2005. 1, 3
[7] B. Horn. Robot vision. MIT Press, Cambridge, MA, 1986. 1[8] K. Ikeuchi. Determining surface orientations of specular surfaces byusing the photometric stereo method. IEEE Trans. Pattern Analysis
and Machince Intelligence, 3(6):661–669, 1981. 1[9] S. Koppal and S. G. Narasimhan. Clustering appearance for scene analysis. In Proc. IEEE Conf. Computer Vision and Pattern Recognition,
pages II–1323–1330, 2006. 2[10] P. Kovesi. Shapelets correlated with surface normals produce surfaces.
In Proc. Int’l. Conf. Computer Vision, pages 994–1001, 2005. 5[11] K. V. Mardia, J. T. Kent, and J. M. Bibby. Multivariate Analysis. Academic Press, London, 1979. 5[12] S. K. Nayar, K. Ikeuchi, and T. Kanade. Determining shape and reﬂectance of hybrid surfaces by photometric sampling. IEEE Trans.
Robotics and Automation, 6(4):418–443, 1990. 1[13] M. Oren and S. Nayar. Generalization of the lambertian model andimplications for machine vision. Int’l J. Computer Vision, 14(3):227–
251, 1995. 3[14] R. Pless and I. Simon. Embedding images in non-ﬂat spaces. In
Proc. Conf. Imaging Science, Systems, and Technology, pages 182–
188, 2002. 3[15] J. B. Tenenbaum, V. Silva, and J. C. Langford. A global geometricframework for nonlinear dimensionality reduction. Science, 290:2319–
2323, 2000. 3, 5[16] K. E. Torrance and E. M. Sparrow. Theory for off-specular reﬂection
from roughened surface. J. Optical Society of America, 57:1105–1114,
1967. 3[17] G. J. Ward. Measuring and modeling anisotropic reﬂection. Proc.
SIGGRAPH ’92, pages 265–272, 1992. 3[18] H. Winnemoller, A. Mohan, J. Tumblin, and B. Gooch. Light waving:
estimating light positions from photographs alone. In Proc. Eurographics 2005, pages 433–488, 2005. 3[19] L. Wolff, S. Shafer, and G. E. Healey. Physics-based vision: Principals
and Practice, Shape Recovery. Jones and Bartlett., 1992. 1[20] R. J. Woodham. Analysing images of curved surfaces. Artif. Intell.,17(1-3):117–140, 1981. 1[21] R. Zhang, P.-S. Tsai, J. E. Cryer, and M. Shah. Shape from shading: A survey. IEEE Trans. Pattern Analysis and Machine Intelligence,21(8):690–706, 1999. 1