CONFIDENTIAL EXTENDED ABSTRACT.
DO NOT DISTRIBUTE ANYWHERE.

The 23rd Meeting on Image Recognition and Understanding

Self-calibrating Near-light Photometric Stereo underAnisotropic Light Emission
H ENG G UO1,a)B OXIN S HI2M ICHAEL WAECHTER3YASUYUKI M ATSUSHITA1

Abstract
In conventional calibrated photometric stereo settings, a lightsource is assumed to be directional or to have an isotropic radiant intensity distribution (RID). These assumptions are hard toachieve in practice, and the deviation brings inaccurate shape recovery. In this work, we propose a self-calibrating near light photometric stereo method to jointly estimate the object shape and theanisotropic point light RID. Speciﬁcally, our method takes image
observations, point light positions and sparse depth of the target object as input, and iterates the procedures of the anisotropic
RID estimation and shape reconstruction. Experimental resultson both synthetic and real data demonstrate the effectiveness of
our method.

1.

Introduction

Photometric stereo recovers surface orientations and albedos
from a set of images taken from a ﬁxed viewpoint under different
light conditions. Early photometric stereo approaches [2, 10] assume distant and uniform light source, which is hard to perfectly
achieve in reality. Near-light photometric stereo [5, 8] setup relieves such assumptions by taking the effects of light falloff andspatially varying light directions into account.
To achieve accurate shape recovery, precise calibration of theinvolved point light sources is required, which includes their 3D
positions and the radiant intensity distributions (RIDs), i.e., radiant intensity of a light source under different emitted directions. Existing near light photometric stereo methods either assume isotropic point light emission [4, 12], or use a calibration
object to estimate the parametric RID models before applying
photometric stereo algorithms [7, 8, 11]. However, existing calibration approaches have some disadvantages: 1) the radianceemitted by a point light typically changes over time due to rampup time and camera auto-exposure [5], leading to different photometric observations between the RID calibration step and thefollowing image acquisition step in photometric stereo. 2) current RID representation models used in RID calibration [7, 11]
and near-light photometric stereo [5, 6, 8] are limited to radiallysymmetric RIDs, which is not ﬂexible enough to handle diversity
of RIDs in the real world. To address these problems, we propose
a self-calibrating photometric stereo framework in which object
1
2
3
a)

Osaka University
Peking University
Facebookheng.guo@ist.osaka-u.ac.jp

Object

𝒙𝒙

𝝎𝝎𝑠𝑠
RID

𝑅𝑅1

𝒔𝒔 Point light

Image
plane
𝑅𝑅2
CameraFigure 1: Illustration of a near-light scene under an anisotropic
radiant intensity distribution (RID).
shape reconstruction and RID calibration are conducted simultaneously without the use of reference objects. Our key observationis that the observed RID as shown in Fig. 1 is located in a local
region when illuminating an object. We assume this local RID
is smooth and can be interpolated by spherical harmonic basis
which is widely used for low-frequency illumination approximation. Based on the RID representation, we aim to estimate theshape of a Lambertian object with given calibrated point light positions and an initial estimate of sparse depth map for the target
object.
This paper makes the following contributions.
• We propose a near-light photometric stereo method with RID
auto-calibration without using a reference object.
• We propose an easy-to-ﬁt RID representation in a local region for real-world anisotropic light emission.
• We show our self-calibrated near-light photometric stereo
solution achieves higher accuracy on shape recovery thanstate-of-the-arts.

2.

Proposed method

As shown in Fig. 2, our method takes image observations,point light positions, and sparse depth as inputs, and output object
shape, albedo and the local RID for each point light source. Ouralgorithm includes two modules: radiant intensity distribution estimation and dense surface reconstruction, i.e., we ﬁrst calculate
the light RID, then apply it to the surface reconstruction process
to optimize dense surface normal and depth alternatively. In thefollowing sections, we will introduce the two modules in detail.
1

The 23rd Meeting on Image Recognition and Understanding

2.1 Radiant intensity distribution estimation
Under lambertian reﬂectance, a measured irradiance illumiatedby the i-th near point light source si can be written asmi (p) = b⊤ (p)li (p), i ∈ {1, · · · , f }Near-light Image
Formation Model

RID Representation
Model

(1)
Image observations

where b(p) = ρ(p)n(p) ∈ R is an albedo-scaled surface normal vector at image coordinate p. Considering light falloff andanisotropic radiant intensity distribution, the lighting vector at position p induced by the point light si can be deﬁned as
3li (p) = ei (p)s − x(p),∥si − x(p)∥32

si − x(p)),∥si − x(p)∥

(3)

h
∑
si − x(p)cik yk ( i),∥s − x(p)∥

(4)

k=1

where cik and yk (·) denote the coefﬁcient and the basis of the k-th
spherical harmonics for the i-th point light source. As shown in
Fig. 1, the light direction from point light to the object fall into alimited region R1 . For photometric stereo, we mainly care about
the radiant intensity distribution in this local region and it is reasonable to assume the RID in this patch is smooth enough to beinterpolated by limited order of spherical harmonic (SH) basis.
Experiments shown in Fig. 3 support our observations. Here
we use the top-9 coefﬁcients of SH to ﬁt real world RIDs undera local region deﬁned in the range of 50 degrees deviated fromthe light principal direction. Compared with cosine-power based
method used in [8], our SH-based RID representation model
achieve higher ﬁtting accuracy in anisotropic RID cases.
If we know surface point position x(p) and the light position*1(𝑥𝑥2 , 𝑦𝑦2 , 𝑧𝑧2 )Sparsenormal map(𝑥𝑥𝑝𝑝 , 𝑦𝑦𝑝𝑝 , 𝑧𝑧𝑝𝑝 )Light RID

Near-light Image
Formation Model

Light position

Normal(2)where ω i is the principal direction of point light, µi is the attenuation coefﬁcient which can be estimated from light source
datasheet. This RID representation model is limited to radiallysymmetric RIDs and also requires light source manufacturing information and additional calibration for the point light principle
direction. To handle these problems, we propose a new parametric RID representation model which is able to handle various
RIDs and has no need of extra calibration.
Since RIDF can be deﬁned as a function in the spherical coordinate system, it is convenient to represent it with the linear
combination of spherical harmonic (SH) basis, as shown below,
ei (p) =RID Estimation(𝑥𝑥1 , 𝑦𝑦1 , 𝑧𝑧1 )Albedoiwhere x(p) ∈ R3 is the surface point position corresponding to
pixel position p, and ei (p) ∈ R represents the irradiance on thispoint emitted from i-th point light. We describe the light irradiance by the radiant intensity distribution fucntion (RIDF) ei (·).
In general, each point light source has its independent RIDF
ei (·). An isotropic point light source emits light equally in all
directions. Therefore, the RIDF ei (·) is a constant function. In
real-world, most RIDs are anisotropic with regard to the light directions. Existing methods [6, 8] use an exponential function withcosine as basis to represent the RID, as shown in Eq. (3)ei (p) = cosµi (ω i ·Object shapehttps://www.nichia.co.jp/en/product/led.html

Depth

Poisson Surface
Reconstruction
Sparse depth

Dense Surface
Reconstruction

Local RID

Input

MethodOutputFigure 2: The pipeline of our method.
SK6812RGBW

NCSB119B-V1

0°

-50°

0°

50°

-50°

1.0 0.8 0.6 0.4 0.2 0.0

1.0 0.8 0.6 0.4 0.2 0.0

NSPWF50DS

NESB146A0°-50°

1.0 0.8 0.6 0.4 0.2 0.0

50°

0°

50°

QY18 [8]

-50°

Ours

50°

1.0 0.8 0.6 0.4 0.2 0.0
GTIsotropic

Figure 3: The representation power comparison of different approaches on four real world RIDs extracted from NACHIA Corp
datasheet *1 . The product name of the point lights are shown on
the top of each RID plot. Notice that the point light “NESB146A”
is not radially symmetric.
si , we formulate the image formation constraint about albedoscaled normal vector b(p) and inverse RID f i (p) following
Eq. (1) and Eq. (2).
Ef {b(p), f i (p)} = {mi (p)∥si −x(p)∥32 f i (p)−b⊤ (p)(si −x(p))}2(5)where the inverse RID f i (p) = 1/ei (p) is also a function deﬁned in the spherical coordinate, which can be interpolated by
SH basis. Similarly, we deﬁne the inverse RID constraint about
the inverse RID f i (p) and the corresponding SH coefﬁcient of
dik from Eq. (4)Er {f i (p), dik } = {f i (p) −

h
∑
k=1

dik yk (si − x(p))}2 (6)∥si − x(p)∥Since we have the input depth z at sparse positions p̂, we canrecover the corresponding 3D sparse surface point position basedon perspective camera model as shown in Eq. (7).

2The 23rd Meeting on Image Recognition and Understandingx(p̂) = z(p̂)K−1 p̂h ,(7)where K and p̂h are camera intrinsic matrix and homogeneous
image coordinates.
With point light positions si , sparse surface points x(p̂) andimage observations mi (p), we solve sparse albedo scaled surface normal, inverse RID and its correspnding SH coefﬁcients foreach point light by addressing both image formation constraint
Eq. (5) and RID constraint Eq. (6).
{b(p̂), f i (p̂), dik } =f∑∑(Ef {b(p̂), f i (p̂)} + λEr {f i (p̂), dik })(8)p̂∈Ω i=1

λ is the weight to balance the two constraints. The above energy function can be formulated as a homogeneous linear equation, which can be easily solved by singular value decomposition(SVD) up to a scale ambiguity. The surface normal and albedo
can be then calculated from b by applying per-pixel normalization.
2.2 Dense surface reconstruction
From last section we obtain sparse depth z(p̂) and n(p̂). Following [3], the complete surface normal ﬁeld can be calculatedon the distinct patches containing such sparse oriented points. We
obtained dense surface normal from [3] and refer it as initial surface normal n0 (p). Next, we reﬁne the surface normal and depth
iteratively guided by the image formation constraint and Poisson
equation.
Let’s denote depth in logarithmic scale as zl (p) = ln z(p). As
stated in [9], the relationship between depth zl (p) and surface
normal n(p) under the perspective camera model is,∇zl (p) = [ψ(p), γ(p)]⊤ ,






n1 (p)
ψ(p) = −
,
u(p)n
(p)
+
v(p)n
1
2 (p) + fl n3 (p)
n2 (p)

γ(p) = −,u(p)n1 (p) + v(p)n2 (p) + fl n3 (p)z ξ+1 (p) with updated surface normal map. We conduct this alternative optimization process until convergence or the iteration
times reaching the predeﬁned maximum number.

3.

Experiment

We evaluate our method on both synthetic and real data in thesetting of anisotropic near-light illumination.
3.1 Result on synthetic data
We render three synthetic objects named “Bunny”, “Fattony”
and “Jolie” with the size of 14 cm, and illuminate the scene with
49 virtual point lights ﬁxed on a planer board placed 30 cm awayfrom the object. The same real-world RID are assigned to all the
point lights on the board. We evaluate the depth surface normal
and RID estimated by our method and the semi-calibrated photometric stereo algorithm proposed in [8], in which the cosinebased model is used for RID representation. To get a fair comparison, we interpolate the input sparse depth to a complete depth
map as the initial depth in [8]. To calculate the accuracy of depth,
normal and RID, we apply the mean absolute error, mean angular error and mean squared error between the estimation and theground truth, respectively. From Table 1 to Table 3 we can see,our self-calibrating photometric stereo method outperforms [8]on both shape and RID estimation.
Table 1: Comparison on depth estimation evaluated in mm
MethodRID

Bunny

FatTony

Julia

Ours

NCSB119B-V1
NSPWF50DS
SK6812RGBW

0.031
0.016
0.016

0.034
0.017
0.017

0.025
0.015
0.016

QY18 [8]NCSB119B-V1
NSPWF50DS
SK6812RGBW

0.767
1.381
1.309

1.171
0.530
1.301

2.104
1.154
1.376

Table 2: Comparison on normal estimation evaluated in degree.
(9)

where fl is the camera focal length, u(p) and v(p) represent
the horizontal and vertical image coordinate at pixel position p,n(p) = [n1 (p), n2 (p), n3 (p)]⊤ .
Applying the divergence operator, we can transform the above
equation into a Poisson equation. By solving the following linear system, we can recover the complete depth map based on theinitial surface normal map n0 (p) and sparse depth z(p̂).
)
( )
(
ηI
ηzl (p̂).
(10)zl (p) = ∂ψ∂γ∆∂u + ∂v
where η deﬁnes the weight for the sparse depth constraint, experimentally we set it to 1. After obtaining log-scaled zl (p), the realcomplete depth map can be obtained by z(p) = ezl (p) .
We refer the complete depth map computed from Eq. (10) as
ξz (p). Fixing the RIDF, we then update the albedo scaled surface
normal b(ξ+1) (p) by addressing the image formation constraint
shown in Eq. (8). After that, we use Eq. (10) to calculate the depth

MethodRID

Bunny

FatTony

Julia

Ours

NCSB119B-V1
NSPWF50DS
SK6812RGBW

0.022
0.238
0.097

0.039
0.013
0.191

0.096
0.283
0.319

QY18 [8]NCSB119B-V1
NSPWF50DS
SK6812RGBW

1.495
2.249
2.505

3.509
2.017
3.031

5.194
3.008
3.454

Table 3: Comparison on RID estimation.
MethodRID

Bunny

FatTony

Julia

Ours

NCSB119B-V1
NSPWF50DS
SK6812RGBW

1.023e-6
2.018e-51.325e-5

1.227e-6
1.285e-7
3.675e-57.963e-6
1.407e-53.467e-5

QY18 [8]NCSB119B-V1
NSPWF50DS
SK6812RGBW

4.717e-51.559e-4
1.565e-43.383e-4
6.125e-52.590e-4

3.135e-48.416e-51.603e-4

3.2 Result on real data
Figure 4 shows our capture setup. It includes two Canon EOS
5D Mark IV cameras with 85mm lens and 256 point lights with
their positions calibrated. The baseline between two cameras is
165mm. We apply semi-global stereo matching algorithm [1] to

3The 23rd Meeting on Image Recognition and Understanding

Barrel

Ours

LED
array

QY18 [8]

Pot

Ours

StereocameraQY18 [8]Input images

Figure 4: Our capture setup.

[5]

[6]

[7]

[8]

Acknowledgement

This work is supported by JSPS KAKENHI Grant Number
JP19H01123.

References[1] Hirschmuller, H.: Stereo processing by semiglobal matching and mutual information, IEEE Transactions on patternanalysis and machine intelligence, Vol. 30, No. 2, pp. 328–
341 (2007). 3[2] Ikehata, S., Wipf, D., Matsushita, Y. and Aizawa, K.:
Robust photometric stereo using sparse regression, 2012
IEEE Conference on Computer Vision and Pattern Recognition, IEEE, pp. 318–325 (2012). 1[3] Kazhdan, M., Bolitho, M. and Hoppe, H.: Poisson sur*2[4]

Conclusion

In this paper, we propose a self-calibrating near-light photometric stereo method to handle shape recovery under anisotropic
point lights. Our key observation is that in the photometric stereo
task the interest RID is located in a local region. Based on thisobservation, we propose a RID representation model based on
spherical harmonic basis which is ﬂexible to various type of RIDs
and without extra calibration. Given the sparse depth of the target object as input, our recovered RID and complete object shapeachieve higher accuracy compared to the existing methods.

5.

Front view

AlbedoNormalFigure 5: Shape and albedo comparison between our method and Quéau et al. [8]on real data

calculate the disparity between the stereo image pairs and further obtain the initial sparse depth of the target. Our system isequipped with point light “SK6812RGBW” from OPSCO Optoelectronics Corp*2 with its theoretical RID shown in Fig. 3. The
distance between the object and light sources is roughly 0.8m.
As shown in Fig. 5, we capture two objects “Barrel” and “Pot”
illuminated by 40 of 256 point lights, and then compare with [8]on shape and albedo estimation. From the sideview in Fig. 5 we
can see that our method with SH-based RID representation model
outputs more accurate surface shape in “Barrel” case and achieves
comparable results on “Pot” object.

4.

Side view

[9]

[10]

[11]

[12]

face reconstruction, Proceedings of the Eurographics symposium on Geometry processing, Vol. 7 (2006). 3
Liu, C., Narasimhan, S. G. and Dubrawski, A. W.: Nearlight photometric stereo using circularly placed point lightsources, 2018 IEEE International Conference on Computational Photography (ICCP), IEEE, pp. 1–10 (2018). 1
Logothetis, F., Mecca, R. and Cipolla, R.: Semi-calibrated
near ﬁeld photometric stereo, Proceedings of the Conference on Computer Vision and Pattern Recognition (CVPR),
pp. 941–950 (2017). 1
Logothetis, F., Mecca, R., Quéau, Y. and Cipolla, R.: Nearﬁeld photometric stereo in ambient light, Proceedings ofthe British Machine Vision Conference (BMVC) (2016). 1,
2
Park, J., Sinha, S. N., Matsushita, Y., Tai, Y.-W. and
So Kweon, I.: Calibrating a non-isotropic near point lightsource using a plane, Proceedings of the Conference on
Computer Vision and Pattern Recognition (CVPR), pp.
2259–2266 (2014). 1Quéau, Y., Durix, B., Wu, T., Cremers, D., Lauze, F.
and Durou, J.-D.: Led-based photometric stereo: modeling, calibration and numerical solution, Journal of Mathematical Imaging and Vision, Vol. 60, No. 3, pp. 313–340(2018). 1, 2, 3, 4Quéau, Y., Durou, J.-D. and Aujol, J.-F.: Normal integration: a survey, Journal of Mathematical Imaging and Vision, Vol. 60, No. 4, pp. 576–593 (2018). 3
Shi, B., Matsushita, Y., Wei, Y., Xu, C. and Tan, P.: Selfcalibrating photometric stereo, 2010 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, IEEE, pp. 1118–1125 (2010). 1
Visentini-Scarzanella, M. and Kawasaki, H.: Simultaneous
camera, light position and radiant intensity distribution calibration, Image and Video Technology, Springer, pp. 557–
571 (2015). 1
Xie, L., Xu, Y., Zhang, X., Bao, W., Tong, C. and Shi,
B.: A self-calibrated photo-geometric depth camera, The
Visual Computer, Vol. 35, No. 1, pp. 99–108 (2019). 1

https://www.opscoled.com/en/about/index.html

4