# On the Regularising Levenberg-Marquardt Method for Blinn-Phong Photometric
Stereo

Georg Radow and Michael Breuß
Chair of Applied Mathematics
Brandenburg University of Technology Cottbus-Senftenberg
<EMAIL_ADDRESS>
###### Abstract

Photometric stereo refers to the process to compute the 3D shape of an object
using information on illumination and reflectance from several input images
from the same point of view. The most often used reflectance model is the
Lambertian reflectance, however this does not include specular highlights in
input images. In this paper we consider the arising non-linear optimisation
problem when employing Blinn-Phong reflectance for modeling specular effects.
To this end we focus on the regularising Levenberg-Marquardt scheme. We show
how to derive an explicit bound that gives information on the convergence
reliability of the method depending on given data, and we show how to gain
experimental evidence of numerical correctness of the iteration by making use
of the Scherzer condition. The theoretical investigations that are at the
heart of this paper are supplemented by some tests with real-world imagery.

##  1 Introduction

The photometric stereo (PS) problem is a fundamental task in computer vision
[5]. The aim of PS is to infer the 3D shape of an object from a set of
multiple images. Thereby the images depict an object from the same
perspective, but the illumination direction changes throughout the images. An
important information besides the illumination is the light reflectance of the
object. The classic PS model [14, 13] is formulated in terms of Lambertian
light reflectance. A Lambertian surface is characterised by diffuse
reflectance and the independence of perceived shading from the viewing angle.
The Lambertian set-up is certainly convenient for modeling, as it represents
the most simple mathematical model for reflectance, and thus resulting formula
and inverse problems are relatively simple. However, it is quite well known
that in PS specular highlights [6] as well as non-Lambertian diffuse effects
[7] may have an important impact on 3D reconstruction.

Let us also comment on some other basic characteristics of PS. Depending on
the knowledge on the lighting, one discerns between calibrated and
uncalibrated PS. In this work we consider only the calibrated case, where
lighting directions and intensities are known. Furthermore, the final goal of
PS is to obtain a depth map, such that for each relevant image pixel three-
dimensional information of the depicted object is obtained. While some
approaches tackle this problem directly in terms of depth values [8], the more
common strategy is to divide depth computation into two sub-problems. In doing
so at first a map of normal vectors is computed, from which the (relative)
depth is obtained in a second step. See for instance [11] for a survey on
surface normal integration. In this paper we only consider the first of the
latter tasks, that is to find the normal vectors. Another aspect is sometimes
the projection performed by the camera during image acquisition, often leading
to orthographic or perspective models, respectively. In this work we address
effectively both settings.

Our contribution. In this paper, we consider some theoretical aspects of
practical value in the optimisation of PS when using Blinn-Phong reflectance.
Here we extend in several ways upon previous work; let us especially refer to
[6], where the Blinn-Phong model is employed in a similar way as here.
Thereby, we consider to include the potentially most important specularity
parameter, the so-called shininess, as an unknown in the optimisation, which
is in contrast to [6] and many other works in the field. The approximate
solution of the non-linear optimisation problem arising pixel-wise is
performed by the regularising Levenberg-Marquardt method, see especially [2].
As this is an iterative method, it is important to assess the influence of
initialisation on the convergence and to give a rigorous bound as a stopping
criterion. Furthermore as the problem is non-linear, one can observe in
practical examples, that it may be difficult to minimise the underlying
residual. To address this issue we investigate the use of a coarse-to-fine
(CTF) scheme as well as an initialisation obtained through classical PS. We
show how to explore Scherzer’s criterion [3], which appeared in [4] for the
first time. This criterion is considered for theoretical purposes within the
construction of the method, in order to assess the convergence property in our
PS problem experimentally.

##  2 Classical Photometric Stereo

Let us reiterate the classic PS approach of Woodham [14, 13]. Given is a set
of $m\geq 3$ images
$\left(\mathcal{I}_{1},\dots,\mathcal{I}_{m}\right)^{\top}\eqqcolon\mathcal{I}$,
so that $\mathcal{I}:\Omega\rightarrow\mathbb{R}^{m}$, along with the
corresponding lighting directions $L_{k}\in\mathbb{R}^{3}$ with $\left\lVert
L_{k}\right\rVert=1$ for $k=1,\dots,m$, with associated intensities $l_{k}\geq
0$. Throughout the paper $\lVert\>\cdot\>\rVert$ denotes the Euclidean norm or
the induced spectral norm. The object to be reconstructed is depicted usually
as a non-rectangular domain $\Omega\in\mathbb{R}^{2}$, which is embedded in
the image domain.

The surface normal vectors $\mathcal{N}:\Omega\rightarrow\mathbb{R}^{3}$ with
$\lVert\mathcal{N}(x,y)\rVert=1$ for all $(x,y)^{\top}\in\mathbb{R}^{2}$ and
the albedo $\rho^{d}:\Omega\rightarrow\mathbb{R}$ are fitted through a least
squares approach, by minimising

$\iint_{\Omega}\left\lVert\mathcal{R}^{\text{L}}(x,y)-\mathcal{I}(x,y)\right\rVert^{2}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!x\mathop{\mathrm{\mathstrut d}}\nolimits\\!y,$ (1)

with reflectance function
$\mathcal{R}^{\text{L}}\coloneqq\left(\mathcal{R}^{\text{L}}_{1},\dots,\mathcal{R}^{\text{L}}_{m}\right)^{\top}$,
consisting of components

$\mathcal{R}^{\text{L}}_{k}\coloneqq\rho^{d}l_{k}L^{\top}_{k}\mathcal{N},\quad
k=1,\dots,m.$ (2)

In practice this boils down to finding a local solution $N\in\mathbb{R}^{3}$
at every sample location $(x,y)^{\top}$ for the problem

$\min_{N}\left\lVert LN-I\right\rVert^{2},\quad
L\coloneqq\begin{pmatrix}l_{1}L_{1}^{\top}\\\ \vdots\\\
l_{m}L_{m}^{\top}\end{pmatrix},\quad I\coloneqq\mathcal{I}(x,y).$ (3)

This, in turn, leads to the computation of the normal vectors and, as a
byproduct, the albedo according to

$N=\left(L^{\top}L\right)^{-1}L^{\top}I,\qquad\rho^{d}(x,y)=\left\lVert
N\right\rVert,\qquad\mathcal{N}(x,y)=N/\left\lVert N\right\rVert\,.$ (4)

##  3 Blinn-Phong Photometric Stereo

In the general least squares approach Eq. 1, we can modify the reflectance
function to account for non-Lambertian effects. To this end we investigate the
Blinn-Phong (BP) model [9, 1], which has the form
$\mathcal{R}^{\text{BP}}\coloneqq\left(\mathcal{R}^{\text{BP}}_{1},\dots,\mathcal{R}^{\text{BP}}_{m}\right)^{\top}$
with components

$\mathcal{R}^{\text{BP}}_{k}\coloneqq\rho^{d}l_{k}L^{\top}_{k}\mathcal{N}+\rho^{s}h_{k}\max\left\\{0,\mathcal{H}^{\top}_{k}\mathcal{N}\right\\}^{\alpha},$
(5)

$k=1,\dots,m$. We observe by (5) that in the BP model, diffuse reflection as
in (2) is supplemented by a specular reflection term. Here
$\rho^{s}:\Omega\rightarrow\mathbb{R}$ denotes the specular albedo. Another
material parameter is the specular sharpness or shininess
$\alpha:\Omega\rightarrow\mathbb{R}$. The halfway vectors
$\mathcal{H}_{k}:\Omega\rightarrow\mathbb{R}^{3}$ depend on the viewing
directions $\mathcal{V}:\Omega\rightarrow\mathbb{R}^{3}$ and are computed for
$k=1,\dots,m$ as

$\mathcal{H}_{k}(x,y)\coloneqq H_{k}/\left\lVert H_{k}\right\rVert,\qquad
H_{k}\coloneqq L_{k}+\mathcal{V}(x,y).$ (6)

Making use of focal length $f$, the viewing directions $\mathcal{V}^{\perp}$
and $\mathcal{V}^{\angle}$ in the orthographic and perspective setting
respectively are

$\mathcal{V}^{\perp}=(0,0,1)^{\top},\qquad\mathcal{V}^{\angle}(x,y)=(x,y,f)^{\top}.$
(7)

We reinterpret $l_{k}$ as diffuse intensity of the light source and denote
$h_{k}\geq 0$ as specular intensity. To ensure that image intensities are only
increased due to diffuse and specular terms, it is reasonable to enforce
$\rho_{d},\rho_{s}\geq 0$. Furthermore $\rho_{d},\rho_{s}\leq 1$ ensures that
at most as much image intensity is added as light intensity is supplied by
each light source. Finally, it is reasonable to enforce $\alpha>1$ to actually
produce specular highlights through the specular term.

The BP model was originally proposed for computer graphics. It is not based on
physical laws, but it enables to create plausible images with a still simple
model compared to other possible approaches. Despite its simplicity, for use
in inverse problems in computer vision, the non-linearities in Eq. 5 may pose
considerable hurdles.

Let us now discuss the modeling of the components in Eq. 5 along with a few
adaptations we employ. First we turn our attention to the normal vectors
$\mathcal{N}$. One may model them through derivatives of the depth or its
logarithm. In this approach we may parametrise them at a specific location
through depth derivatives $p,q$ as

$\mathcal{N}(x,y)=\frac{N(p,q)}{\left\lVert N(p,q)\right\rVert}.$ (8)

However the step of obtaining a normal vector of length 1 in Eq. 8 adds
another layer of non-linearity to the model. In numerical experiments we found
this approach to be not very reliable. Therefore we opt for an approach in
analogy to classical PS. In Eq. 5 we replace $\rho^{d}\mathcal{N}=N$
introducing the auxiliary variable $r=\rho^{s}/(\rho^{d})^{\alpha}$.By
furthermore replacing $\alpha=1+\exp(a)$ we ensure that
$\mathcal{R}^{\text{BP}}$ has continuous first derivatives. Eq. 5 then takes
the form

$\mathcal{R}^{\text{BP}}_{k}(N,r,a)=l_{k}L^{\top}_{k}N+rh_{k}\max\left\\{0,\mathcal{H}^{\top}_{k}N\right\\}^{1+\exp(a)},$
(9)

with $r,a\in\mathbb{R}$ and $N\in\mathbb{R}^{3}$.

##  4 On the Optimisation Strategy

With BP reflectance, we have to solve a non-linear least squares problem, to
which end we utilise the regularising Levenberg-Marquardt (RLM) scheme [2, 3].
Writing the underlying task in standard notation, with this algorithm one may
aim to find a solution $\vec{x}$ of the problem

$F(\vec{x})=\vec{y},\qquad F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m},$ (10)

with a known differentiable function $F$. Let us note that the description and
discussion of the RLM algorithm in [3] is in a more general setting. For
simplicity we only give an overview of the algorithm based on finite
dimensional spaces, as is fitting for the problem at hand.

It is furthermore assumed that the original data $\vec{y}$ is not known, but
with some $\delta>0$ an estimate is required on how good the given data
$\vec{y}^{\delta}$ approximates the original data, according to

$\left\lVert\vec{y}^{\delta}-\vec{y}\right\rVert\leq\delta.$ (11)

Then with some starting point $\vec{x}_{0}$ the iterative rule takes the form

$\vec{x}_{k+1}=\vec{x}_{k}+\left(F^{\prime}(\vec{x}_{k})^{\top}F^{\prime}(\vec{x}_{k})+\alpha_{k}I_{n}\right)^{-1}F^{\prime}(\vec{x}_{k})^{\top}\left(\vec{y}^{\delta}-F(\vec{x})\right)$
(12)

with Jacobian matrix $F^{\prime}$, $n\times n$-dimensional identity matrix
$I_{n}$ and a regularisation weight $\alpha_{k}>0$ such that with a
preassigned $\rho\in(0,1)$ the new iterate $\vec{x}_{k+1}$ fulfils

$\left\lVert\vec{y}^{\delta}-F(\vec{x}_{k})-F^{\prime}(\vec{x}_{k})\left(\vec{x}_{k+1}-\vec{x}_{k}\right)\right\rVert=\rho\left\lVert\vec{y}^{\delta}-F(\vec{x}_{k})\right\rVert.$
(13)

The _stopping criterion_ of the RLM scheme depends explicitly on the noise
level $\delta$ in the given data. To stop at an iterate $\vec{x}_{k}$, it has
to fulfil

$\left\lVert\vec{y}^{\delta}-F(\vec{x}_{k})\right\rVert\leq\tau\delta,$ (14)

with a preassigned $\tau>2$, fulfilling $\rho\tau>1$. For numerical
experiments we set $\rho=0.5,~{}\tau=2.5$, following [3].

The discussion of the RLM scheme in [3] relies on the strong Scherzer
condition [4]. For the Jacobian matrices at two points
$\vec{x}_{1},\vec{x}_{2}\in\mathbb{R}^{n}$ there exists a matrix
$R=R(\vec{x}_{1},\vec{x}_{2})$ such that
$F^{\prime}(\vec{x}_{1})=RF^{\prime}(\vec{x}_{2})$ and

$\left\lVert R-I_{m}\right\rVert\leq
C^{R}\left\lVert\vec{x}_{1}-\vec{x}_{2}\right\rVert$ (15)

with some $C^{R}>0$, which is constant for all
$\vec{x}_{1},\vec{x}_{2}\in\mathbb{R}^{n}$. This condition imposes a certain
regularity of the Jacobian matrix $F^{\prime}$. In this context we are
interested in a local approximation of $C^{R}$. For two consecutive iterations
$\vec{x}_{k},\vec{x}_{k+1}$ we estimate $R$ as a solution of
$F^{\prime}(\vec{x}_{k})=R(\vec{x}_{k},\vec{x}_{k+1})F^{\prime}(\vec{x}_{k+1})$
with minimal norm. Then we can locally approximate the constant in Eq. 15 as

$C_{k}^{R,\text{loc}}=\frac{\left\lVert
R(\vec{x}_{k},\vec{x}_{k+1})-I_{m}\right\rVert}{\left\lVert\vec{x}_{k}-\vec{x}_{k+1}\right\rVert}.$
(16)

Since $F$ in Eq. 10 is nonlinear, we employ a CTF framework. In doing so the
data is scaled to a coarser scale, _i.e._ to a lower resolution. The obtained
result is then used as initialisation on the next finer scale, until we arrive
at the original resolution.

Let us focus on the assumption Eq. 11. The noise level $\delta$ governs the
stopping criterion of the RLM scheme. If Eq. 11 is not fulfilled then the
iterates may actually diverge.

At this point we make the assumption that our data $\mathcal{I}(x,y)$ is a
realisation of the BP model corrupted by additive white Gaussian noise, _i.e._
it can be modelled as

$\mathcal{I}(x,y)=\mathcal{R}(x,y)+\varepsilon(x,y),\quad\text{for
}(x,y)^{\top}\in\Omega.$ (17)

Here $\varepsilon(x,y)$ is a realisation of a multivariate normal
distribution, such that the $m$ components are independent and identically
distributed (i.i.d.) with mean zero and standard deviation $\sigma>0$, the
corresponding density function is

$f(X)=\frac{1}{\sqrt{2\pi}^{m}\sigma^{m}}\exp\left(-\frac{1}{2\sigma^{2}}\sum_{i=1}^{m}X_{i}^{2}\right),$
(18)

_cf._ [10]. The probability that Eq. 11 holds can be computed with the
following result. The proof, which is technical but straightforward, is
included for the readers convenience. The following result is also related to
the Chi distribution.

######  Proposition 1.

Let $m\in\mathbb{N}$, $\delta>0$ and let $\varepsilon$ be a realisation of an
$m$-dimensional multivariate normal distribution with mean zero, standard
deviation $\sigma>0$ and density Eq. 18. The probability of $P\coloneqq
P(\lVert\varepsilon\rVert\leq\delta|\sigma,m)$ can be computed as follows:

  1. (i)

If $m$ is even, then

$P=1-\exp\left(-\frac{\delta^{2}}{2\sigma^{2}}\right)\sum_{i=0}^{\frac{m}{2}-1}\left(\frac{\delta^{2}}{2\sigma^{2}}\right)^{i}\frac{1}{i!}.$
(19)

  2. (ii)

If $m$ is odd, then

$P=\sqrt{\frac{2}{\pi}}\Bigg{(}\frac{1}{\sigma}\int_{0}^{\delta}\exp\left(-\frac{r^{2}}{2\sigma^{2}}\right)\mathop{\mathrm{\mathstrut
d}}\nolimits\\!r\\\
-\exp\left(-\frac{\delta^{2}}{2\sigma^{2}}\right)\sum_{i=1}^{\frac{m-1}{2}}\Bigg{(}\left(\frac{\delta}{\sigma}\right)^{m-2i}\prod_{j=1}^{\frac{m+1}{2}-i}\left(\frac{1}{2j-1}\right)\Bigg{)}\Bigg{)}.$
(20)

###### Proof.

For any continuous probability density $f$ we have

$P=P(\lVert\varepsilon\rVert\leq\delta|\sigma,m)=\int_{\lVert
X\rVert\leq\delta}f(X)\mathop{\mathrm{\mathstrut d}}\nolimits\\!X.$ (21)

Since the density function in Eq. 18 is radially symmetric, this simplifies to

$P=\int_{0}^{\delta}O_{m}(r)f(r,0,\dots,0)\mathop{\mathrm{\mathstrut
d}}\nolimits\\!r,$ (22)

where

$O_{m}(r)=2r^{m-1}\frac{\pi^{\frac{m}{2}}}{\Gamma\left(\frac{m}{2}\right)}$
(23)

denotes the surface area of a sphere with radius $r$ around the origin in
$\mathbb{R}^{m}$. $\Gamma$ denotes the gamma function. Inserting Eq. 18, we
write

$P=\frac{2^{1-\frac{m}{2}}}{\sigma^{m}\Gamma\left(\frac{m}{2}\right)}\int_{0}^{\delta}r^{m-1}\exp\left(-\frac{r^{2}}{2\sigma^{2}}\right)\mathop{\mathrm{\mathstrut
d}}\nolimits\\!r.$ (24)

Since $\int r\exp(r^{2}/(2a))\mathop{\mathrm{\mathstrut
d}}\nolimits\\!r=a\exp(r^{2}/(2a))+c$, for $m>2$ the integral in Eq. 24 can be
simplified by partial integration, _i.e._

$\int_{0}^{\delta}r^{m-2}\cdot
r\exp\left(-\frac{r^{2}}{2\sigma^{2}}\right)\mathop{\mathrm{\mathstrut
d}}\nolimits\\!r\\\
=-\sigma^{2}\left[r^{m-2}\exp\left(-\frac{r^{2}}{2\sigma^{2}}\right)\right]_{r=0}^{\delta}+\sigma^{2}(m-2)\int_{0}^{\delta}r^{m-4}\cdot
r\exp\left(-\frac{r^{2}}{2\sigma^{2}}\right)\mathop{\mathrm{\mathstrut
d}}\nolimits\\!r.$ (25)

We now consider the two cases of $m$ being even or odd.

Let $m\in\mathbb{N}$ be even. Then repeated partial integration of the
integral Eq. 24 leads to

$\begin{split}&\int_{0}^{\delta}r^{m-1}\exp\left(-\frac{r^{2}}{2\sigma^{2}}\right)\mathop{\mathrm{\mathstrut
d}}\nolimits\\!r\\\
&=-\sum_{i=1}^{\frac{m}{2}-1}\sigma^{2i}\prod_{j=1}^{i-1}(m-2j)\left[r^{m-2i}\exp\left(-\frac{r^{2}}{2\sigma^{2}}\right)\right]_{r=0}^{\delta}\\\
&\quad+\sigma^{m-2}\prod_{j=1}^{\frac{m}{2}-1}(m-2j)\int_{0}^{\delta}r\exp\left(-\frac{r^{2}}{2\sigma^{2}}\right)\mathop{\mathrm{\mathstrut
d}}\nolimits\\!r\\\
&=-\sum_{i=1}^{\frac{m}{2}}\sigma^{2i}\prod_{j=1}^{i-1}(m-2j)\left[r^{m-2i}\exp\left(-\frac{r^{2}}{2\sigma^{2}}\right)\right]_{r=0}^{\delta}\\\
&=-\sum_{i=1}^{\frac{m}{2}}\sigma^{2i}\frac{2^{i-1}\left(\frac{m}{2}-1\right)!}{\left(\frac{m}{2}-i\right)!}\left[r^{m-2i}\exp\left(-\frac{r^{2}}{2\sigma^{2}}\right)\right]_{r=0}^{\delta}\\\
&=\sigma^{m}2^{\frac{m}{2}-1}\left(\frac{m}{2}-1\right)!-\sum_{i=1}^{\frac{m}{2}}\sigma^{2i}\frac{2^{i-1}\left(\frac{m}{2}-1\right)!}{\left(\frac{m}{2}-i\right)!}\delta^{m-2i}\exp\left(-\frac{\delta^{2}}{2\sigma^{2}}\right).\end{split}$
(26)

This formula can easily be verified for $m=2$, as in this case the initial
integral simplifies to the form $\int
r\exp(r^{2}/(2a))\mathop{\mathrm{\mathstrut d}}\nolimits\\!r$. Inserting Eq.
26 and $\Gamma(m/2)=(m/2-1)!$ into Eq. 24, we obtain after an index shift Eq.
19.

Now let $m\in\mathbb{N}$ be odd. Again we use repeated partial integration on
the integral in Eq. 24, until we arrive at

$\int_{0}^{\delta}r^{m-1}\exp\left(-\frac{r^{2}}{2\sigma^{2}}\right)\mathop{\mathrm{\mathstrut
d}}\nolimits\\!r=\sigma^{m-1}\prod_{j=1}^{\frac{m-1}{2}}(2j-1)\int_{0}^{\delta}\exp\left(-\frac{r^{2}}{2\sigma^{2}}\right)\mathop{\mathrm{\mathstrut
d}}\nolimits\\!r\\\
-\sum_{i=1}^{\frac{m-1}{2}}\sigma^{2i}\frac{\prod_{j=1}^{\frac{m-1}{2}}(2j-1)}{\prod_{j=1}^{\frac{m+1}{2}-i}(2j-1)}\delta^{m-2i}\exp\left(-\frac{\delta^{2}}{2\sigma^{2}}\right).$
(27)

Plugging this together with

$\Gamma\left(\frac{m}{2}\right)=\Gamma\left(\frac{m-1}{2}+\frac{1}{2}\right)=\frac{\left(m-1\right)!\,\sqrt{\pi}}{\left(\frac{m-1}{2}\right)!\,2^{m-1}}\\\
=\frac{\prod_{j=1}^{\frac{m-1}{2}}\left((2j)(2j-1)\right)\sqrt{\pi}}{2^{\frac{m-1}{2}}\prod_{j=1}^{\frac{m-1}{2}}\left(2j\right)}=\frac{\sqrt{\pi}}{2^{\frac{m-1}{2}}}\prod_{j=1}^{\frac{m-1}{2}}(2j-1)$
(28)

into Eq. 24 leads to Eq. 20. ∎

##  5 Experiments

|  |  |
---|---|---|---
Figure 1:  _(left-to-right:)_ one of the input images of the _sphere_ rendered
using the BP model; colour coded vector field of ground truth normal vectors;
classical PS with preprocessing [15], AAE $1.02$; developed BP framework with
CTF, AAE $0.37$

Since we focus on the computed vector fields of surface normals, it appears
adequate to employ colour coding of surface normals for visual assessment,
_cf._ Figure 1. For quantitative evaluation we consider here the standard AAE,
where the averaging is performed over the object domain. Let us note that we
use the result obtained through classical PS as an initialisation for the BP
model. Throughout the experiments we computed $\delta$ according to
Proposition 1, such that Eq. 11 is fulfilled with a probability of $95\%$. We
observed that the choice of this confidence level is not critical for the
outcome of our experiments.

### Synthetic Test Example.

As a synthetic experiment for our investigations we consider the _sphere_
example, see Figure 1. Let us note that we consider an orthographic setting
for all the _sphere_ experiments. As we observe in Figure 1, in this
experiment the developed computational model and set-up enables to obtain a
nearly perfect result. For optimisation we employed in total $5$ input images,
of which we show here just one example. For comparison, we give here the
corresponding result obtained by Lambertian PS applied at analogous input
images where we filtered the specular highlights by the subspace technique
proposed in [15], which is supposed to make the input nearly Lambertian. As is
confirmed here visually as well as quantitatively, it appears favorable (at
least in this example) to explore an explicit modeling like with the proposed
BP framework.

Let us note that in fact this test example may not be too easy, as can be
observed by the results obtained by preprocessing and Lambertian PS. The
reason is that the specular highlights in the input are not perfectly
distributed over the sphere and may result in distortions if not being
accounted for sufficiently accurate in the model.

### Evaluation of Scherzer’s Condition.

As discussed in Sec. 4, between two iterates of the RLM scheme we observe the
local approximation $C_{k}^{R,\text{loc}}$ of the constant in Eq. 15 according
to Eq. 16. As the Scherzer condition is an important assumption for the
results in [3], we opt to add a break condition, where the algorithm stop if
the estimate grows too large. In practice the algorithm is halted if we
observe an iterate with $C_{k}^{R,\text{loc}}\geq 2000$. As can be seen in
Figs. 2 and 3 this is usually the case at locations where specular highlights
may occur, as the angle between halfway vectors and surface normals becomes
small. One may interprete this result in the way, that the energy that is
minimised features at highlights many small variations that makes it difficult
to obtain a reliable local minimum.

We evaluated the restarting of the RLM scheme with a larger parameter $\rho$
in Eq. 13, if it stopped before an iterate fulfils Eq. 14. This may lead to a
smaller trust region and to a more stable behaviour of the algorithm. However
we did in general not observe a significant increase in quality. The results
displayed here were thus computed without restarting the RLM scheme, giving an
account of the unstabilised version of the method.

---
Figure 2:  Algorithmic behaviour in the _sphere_ experiment _(left)_ and an
example from the _DiLiGent_ data set [12] _(right)_. White depicts the
locations where the RLM scheme stopped due to the $C_{k}^{R,\text{loc}}\geq
2000$ criterion.

### Real World Test Example.

In order to assess the properties and usefulness of the developed numerical BP
framework, we exploit here a selected variety of examples taken from the
_DiLiGent_ data set [12] which gives an account of photographed real-world
objects with different reflectance properties. Here we do not employ a CTF
scheme, as we rely on the initialisation obtained with classical PS. Let us
note that the underlying model is now (in practice, weakly) perspective.

As can be visually assessed by means of Figure 3, the proposed model along
with its adaptations performs very reasonably but in some details not perfect,
depending on the actual example. For clarifying thereby the zones of influence
of the specular terms we depict masks showing the object parts where the BP
model gives an effective contribution. When taking into account the properties
of the considered examples, it appears especially that the broad specularities
as appearing in the input (teddy bear, goblet) may result in a certain
inaccuracy. In turn, when highlights appear but are not too strong (cat, tea
pot), results are quite convincing, given that the underlying reflectance in
these cases is supposed to be non-linear in the diffuse reflectance as the
underlying material is rough. In the tested real world setting from _DiLiGent_
the results are overall of similar quality to the preprocessed Lambertian
method. Therefore we conjecture that our numerical BP framework appears to be
especially suited for dealing with objects with not too strong highlights,
being at the same time able to tackle a certain range of diffuse reflectance
of rough materials.

|  |  |
---|---|---|---
|  |  |
|  |  |
Figure 3:  _(left-to-right:)_ Examples from _Diligent_ data sets. _(top-to-
bottom:)_ Visualisation of ground truth normals; normal fields based on BP
(where we note the effect of the not satisfied Scherzer condition at some
highlights at the goblet); mask based on half directions. White depicts
locations where the maximum of the cosines between halfway vectors and the
normal vector obtained with classical PS is $\geq 0.99$.

##  6 Conclusion

We discussed the BP reflectance in the context of PS. The augmentation of
classical PS with this reflectance model is straightforward, but solving the
arising optimisation problem is less so. This task can be tackled with the RLM
scheme, which leads to satisfactory results.

The findings for the implementation of the RLM scheme may be translated to
other problems, since the assumption that the data follows a normal
distribution is very common. The application of the BP model to more complex
data sets poses considerable hurdles, which may be adressed in future work.

## References

  * [1] Blinn, J.F.: Models of light reflection for computer synthesized pictures. In: Proceedings of the 4th annual conference on Computer graphics and interactive techniques - SIGGRAPH '77. ACM Press (1977)
  * [2] Hanke, M.: A regularizing Levenberg - Marquardt scheme, with applications to inverse groundwater filtration problems. Inverse Problems 13(1), 79–95 (1997)
  * [3] Hanke, M.: The regularizing Levenberg-Marquardt scheme is of optimal order. Journal of Integral Equations and Applications 22(2), 259–283 (2010)
  * [4] Hanke, M., Neubauer, A., Scherzer, O.: A convergence analysis of the Landweber iteration for nonlinear ill-posed problems. Numerische Mathematik 72(1), 21–37 (1995). https://doi.org/10.1007/s002110050158
  * [5] Horn, B.K.P.: Robot Vision. MIT Electrical Engineering and Computer Science, MIT Press (1986)
  * [6] Khanian, M., Boroujerdi, A.S., Breuß, M.: Photometric stereo for strong specular highlights. Computational Visual Media 4(1), 83–102 (2018)
  * [7] McGunnigle, G., Dong, J., Wang, X.: Photometric stereo applied to diffuse surfaces that violate lambert’s law. Journal of the Optical Society of America A 29(4), 627 (mar 2012). https://doi.org/10.1364/josaa.29.000627
  * [8] Mecca, R., Tankus, A., Wetzler, A., Bruckstein, A.M.: A direct differential approach to photometric stereo with perspective viewing. SIAM Journal on Imaging Sciences 7(2), 579–612 (2014). https://doi.org/10.1137/120902458
  * [9] Phong, B.T.: Illumination for computer generated pictures. Communications of the ACM 18(6), 311–317 (1975)
  * [10] Prince, S.J.: Computer Vision: Models, Learning, and Inference. Cambridge University Press (2012)
  * [11] Quéau, Y., Durou, J.D., Aujol, J.F.: Normal integration: A survey. Journal of Mathematical Imaging and Vision 60(4), 576–593 (2017)
  * [12] Shi, B., Mo, Z., Wu, Z., Duan, D., Yeung, S.K., Tan, P.: A benchmark dataset and evaluation for non-Lambertian and uncalibrated photometric stereo. IEEE Transactions on Pattern Analysis and Machine Intelligence pp. 1–14 (2018)
  * [13] Woodham, R.J.: Photometric method for determining surface orientation from multiple images. Optical Engineering 19(1), 134–144 (1980)
  * [14] Woodham, R.J.: Photometric stereo: A reflectance map technique for determining surface orientation from image intensity. In: Nevatiam, R. (ed.) Image Understanding Systems and Industrial Applications. Proceedings of the Society of Photo-Optical Instrumentation Engineers, vol. 155, pp. 136–143. SPIE (1978)
  * [15] Wu, L., Ganesh, A., Shi, B., Matsushita, Y., Wang, Y., Ma, Y.: Robust photometric stereo via low-rank matrix completion and recovery. In: Asian Conference on Computer Vision (ACCV), Lecture Notes in Computer Science, vol. 6494, pp. 703–717. Springer Berlin Heidelberg (2010)