# Orthogonal Illuminations in Two Light-Source Photometric Stereo

Ryszard Kozera, Alexander Prokopenya

► To cite this version:

Ryszard Kozera, Alexander Prokopenya. Orthogonal Illuminations in Two Light-Source Photometric Stereo. 15th IFIP International Conference on Computer Information Systems and Industrial Management (CISIM), Sep 2016, Vilnius, Lithuania. pp.402-415, 10.1007/978-3-319-45378-1_36 hal-01637521

HAL Id: hal-01637521

https://inria.hal.science/hal-01637521

Submitted on 17 Nov 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Orthogonal Illuminations in Two Light-Source Photometric Stereo

Ryszard Kozera and Alexander Prokopenya

Warsaw University of Life Sciences - SGGW
Faculty of Applied Informatics and Mathematics
Nowoursynowska str. 159, 02-776 Warsaw, Poland
ryszard.kozera@gmail.com alexander_prokopenya@sggw.pl

**Abstract.** In this paper we investigate the case of ambiguous shape reconstruction from two light-source photometric stereo based on illuminating the unknown Lambertian surface. So-far this problem is merely well-understood for two linearly independent light-source directions with one illumination assumed as overhead. As already established, a necessary and sufficient condition to disambiguate the entire shape reconstruction process is controlled by the satisfaction of the corresponding second-order linear PDE with constant coefficients in two independent variables. This work extends the latter to an arbitrary pair of light-source directions transforming the above constraint into a special nonlinear PDE. In addition, a similar ambiguity analysis is also performed for a special configuration of two light-source directions assumed this time as orthogonal and contained in the vertical plane. Finally, this work is supplemented by illustrative examples exploiting symbolic computation used within a framework of continuous reflectance map model (i.e. an image irradiance equation) and applied to a genuine Lambertian surfaces.

**Keywords:** Photometric stereo, image irradiance equation, ambiguity in shape reconstruction, computer vision.

## 1 Problem formulation

The task of surface reconstruction from its image (or images) poses one of the major challenges in computer vision. There are two common approaches adopted to tackle this problem. The first one (termed a *stereo-like method*) involves projective geometry tools applied to multiple camera images [1–3]. Solving the so-called *matching problem* constitutes here the gist of the entire reconstruction process based on incorporating the triangulation-like technique [1, 2, 4]. The second approach, coined *shape from shading* [3–5], relies on a *priori* knowledge of the physical properties of light reflectance inherent to the specific material coating the unknown surface *S*. In this model the surface *S* is assumed to be formed as a graph *S* = *graph(u)* of the unknown *C*<sup>*k*</sup> (*k* = 0, 1) function *u*: Ω → ℝ defined over an image domain Ω ⊂ ℝ<sup>2</sup>. Recall that the normal *n(s)* = (*n*<sub>1(s)</sub>, *n*<sub>2(s)</sub>, *n*<sub>3(s)</sub>) to *S* at the point *s* = (*x*, *y*, *u(x,y)*)) ∈ *S* (with
$n_3(s) < 0$) reads as $n(s) = (u_x(x, y), u_y(x, y), -1)$. Assuming a distant light-source (yielding a parallel light-beam along the direction $p = (p_1, p_2, p_3)$) the so-called image irradiance equation

$$R(u_x(x, y), u_y(x, y)) = E(x, y) \qquad (1)$$

encapsulates the relation between reflected light from surface S and absorbed light in the image $\Omega$ - see [4]. Here $E(x, y)$ denotes the intensity of the absorbed light at the image point $(x, y) \in \Omega$. On the other hand, R called reflectance map refers to the intensity of the reflected light from $s \in S$ with normal $n(s)$. In the special case of a Lambertian surface $S_L$, $R_L$ is proportional to $\cos(\alpha)$, where $\alpha$ represents the angle between the normal $n(s)$ and a given light-source direction $p = (p_1, p_2, p_3)$ - see [4]. Consequently, for $S_L$ the corresponding image irradiance equation (1) (over $\Omega = \{(x, y) \in \mathbb{R}^2 : 0 \le E(x, y) \le 1\}$) reads as:

$$\frac{\langle p|n(s) \rangle}{\|p\|\|n(s)\|} = \frac{p_1 u_x(x, y) + p_2 u_y(x, y) - p_3}{\sqrt{p_1^2 + p_2^2 + p_3^2} \sqrt{u_x^2(x, y) + u_y^2(x, y) + 1}} = E(x, y). \quad (2)$$

By Cauchy-Schwartz inequality $|E(x, y)| \le 1$ holds. The extra condition required $E(x, y) \ge 0$ excludes shadowed $(x, y)$-points. The equation (2) yields the first order non-linear PDE (in two variables) modelling a single image shape from shading for $S_L$. Commonly, in solving (2) one usually searches for $u \in C^k$ ($k = 1, 2$), modulo a vertical shift $v = u + c$ (here $c$ is an arbitrary constant). A single image irradiance equation (2) renders generically an ill-posed problem (see [6–10]). Though various mathematical extra assumptions can disambiguate (2) or limit its number of solutions [4, 11–15], they often turn unrealistic or difficult to be met from real camera images.

A feasible remedy is to employ a photometric stereo technique for which the unknown surface $S_L$ (or S) is consecutively illuminated by distant multiple light-sources positioned along linearly independent directions [4, 16]. In this set-up, contrary to the classical stereo method relying on images taken from multiple cameras, only a single camera is deployed.

As demonstrated in [4, 16, 17], three-light source photometric stereo suffices to guarantee a unique surface $S_L = graph(u)$ (up to a vertical shift) with $u \in C^1$ satisfying the corresponding system of three image irradiance equations (2) (formulated for $p, q$ and $r$) over an image $\Omega = \Omega_1 \cap \Omega_2 \cap \Omega_3$. The respective right-hand sides of (2) (i.e. the intensities $E_i \ge 0$ over $\Omega_i$, for $i = 1, 2, 3$) are obtained by consecutive illuminations of $S_L$ along three linearly independent directions $p, q, r \in \mathbb{R}^3$. The entire reconstruction process is decomposed here into two following steps:

a) a gradient computation (an algebraic step):

$$u_x(x, y) = f_1(x, y, E_1, E_2, E_3), \quad u_y(x, y) = f_2(x, y, E_1, E_2, E_3). \qquad (3)$$

In case of three light-sources the resulting vector field $\nabla u = (u_x, u_y)$ satisfying three equations of type (2) is uniquely determined in (3) in terms of image intensities $E_i$ ($i = 1, 2, 3$) and light-source directions $p, q$ and $r$ - see [4, 17].
b) a gradient integration (an analytic step):

$$u(x, y) = u(x_0, y_0) + \int_{\gamma} u_x(x, y)dx + u_y(x, y)dy, \quad (4)$$

where $(x_0, y_0)$, $(x, y) \in \Omega$ and $\gamma: [t_0, t_1] \to \Omega$ is an arbitrary piecewise-$C^1$ curve such that $\gamma(t_0) = (x_0, y_0)$ and $\gamma(t_1) = (x, y)$. Note that $(x_0, y_0)$ is fixed here and $u(x_0, y_0)$ represents any constant $c$. The choice of $\gamma$ joining $(x_0, y_0)$ with varying $(x, y)$ is arbitrary (at least in continuous model of (1)), provided $\Omega$ is simply connected and the vector field (3) fulfills the so-called integrability condition. The latter reads for $u \in C^2$ as

$$u_{xy}(x, y) = u_{yx}(x, y), \quad (5)$$

and for $u \in C^1$ as

$$\int_{\gamma_c} u_x(x, y)dx + u_y(x, y)dy = 0, \quad (6)$$

holding for any piecewise-$C^1$ loop $\gamma_c \subset \Omega$ (with $\gamma_c(t_0) = \gamma_c(t_1)$). Assuming $E_i$ ($i = 1, 2, 3$) are formed by a genuine Lambertian surface $S_L = graph(u)$, a unique vector field $\nabla u$ obtained from (3) is automatically integrable and hence the formula (4) determines an unambiguous $u \in C^1$ (up to a constant).

The case of two-light source photometric stereo requires more intricate analysis (see [17–19]). The respective system of two image irradiance equations:

$$\begin{aligned} \frac{p_1 u_x(x, y) + p_2 u_y(x, y) - p_3}{\sqrt{p_1^2 + p_2^2 + p_3^2} \sqrt{u_x^2(x, y) + u_y(x, y)^2 + 1}} &= E_1(x, y), \\ \frac{q_1 u_x(x, y) + q_2 u_y(x, y) - q_3}{\sqrt{q_1^2 + q_2^2 + q_3^2} \sqrt{u_x^2(x, y) + u_y(x, y)^2 + 1}} &= E_2(x, y) \end{aligned} \quad (7)$$

is solved by the following vector field (see [17, 18]):

$$\begin{aligned} u_x &= \frac{\|p\|(q_1\langle p|q \rangle - p_1\|q\|^2)E_1 + \|q\|(p_1\langle p|q \rangle - q_1\|p\|^2)E_2 + (p_3q_2 - p_2q_3)\varepsilon\sqrt{\Lambda}}{\|p\|(p_3\|q\|^2 - q_3\langle p|q \rangle)E_1 + \|q\|(q_3\|p\|^2 - p_3\langle p|q \rangle)E_2 + (p_1q_2 - p_2q_1)\varepsilon\sqrt{\Lambda}}, \\ u_y &= \frac{\|p\|(q_2\langle p|q \rangle - p_2\|q\|^2)E_1 + \|q\|(p_2\langle p|q \rangle - q_2\|p\|^2)E_2 + (p_1q_3 - p_3q_1)\varepsilon\sqrt{\Lambda}}{\|p\|(p_3\|q\|^2 - q_3\langle p|q \rangle)E_1 + \|q\|(q_3\|p\|^2 - p_3\langle p|q \rangle)E_2 + (p_1q_2 - p_2q_1)\varepsilon\sqrt{\Lambda}}, \end{aligned} \quad (8)$$

where

$$\Lambda = \|p\|^2 \|q\|^2 (1 - E_1^2(x, y) - E_2^2(x, y)) - \langle p|q \rangle (\langle p|q \rangle - 2\|p\|\|q\| E_1(x, y) E_2(x, y)) \quad (9)$$

with the function $\varepsilon(x, y)$ taking values $\pm 1$ so that $f(x, y) = \varepsilon(x, y)\sqrt{\Lambda(x, y)}$ is continuous (for $u \in C^1$) or smooth (for $u \in C^2$).
As discussed in [17, 18] an image domain $\Omega$ is often decomposed into $\Omega = \Omega^{(1)} \cup \Omega^{(2)} \cup \Gamma$ with $\Lambda > 0$ over disjoint sub-domains $\Omega^{(j)} \subset \Omega$ (here $j = 1, 2$) and with $\Lambda = 0$ satisfied along some smooth curve $\Gamma$. Consequently, over each component $\Omega^{(j)}$ either $\varepsilon(x, y) \equiv 1$ or $\varepsilon(x, y) \equiv -1$ hold, yielding the respective two vector fields $(u_x^\pm, u_y^\pm)$ in (8) satisfying (7). Assuming $E_1$ and $E_2$ are generated by a genuine Lambertian surface $S_L$, at least one of these two pairs of vector fields (over each $\Omega^{(j)}$) is integrable.

The ambiguous case of more than one integrable vector field from (8) is well-understood merely for the special configuration of light-source directions, namely when $p = (0, 0, -1)$ and $q = (q_1, q_2, q_3)$ (here $q_1^2 + q_2^2 > 0$) - see [17, 18]. In fact, generically there is only one integrable vector field over each $\Omega^{(j)}$ (and thus over $\Omega$). The latter is governed by the fulfillment of the tight condition (10). Indeed, given two simply connected sub-images $\Omega^{(j)}$, a necessary and sufficient condition enforcing both vector fields $(u_x^+, u_y^+)$ and $(u_x^-, u_y^-)$ to be integrable (and thus yielding the existence of $u^+, u^- \in C^2(\Omega^{(i)})$ solving (7)) reads as:

$$q_1q_2(u_{yy}(x, y) - u_{xx}(x, y)) + (q_1^2 - q_2^2)u_{xy}(x, y) = 0, \quad (10)$$

which is to be satisfied by either $u^+$ or $u^-$. In addition (see [17, 18]), if one of $u^+$ (or of $u^-$) satisfies (10) then so does the complementary one i.e. $u^-$ (or $u^+$). In a rare situation of (10) holding, both pairs $(u^+, u^-)$, determined by (4) over $\Omega^{(j)}$ ($j = 1, 2$), can bifurcate (i.e. can be glued together) along the curve $\Gamma$ to yield either zero or two or four global solutions $u \in C^2$ over the whole image $\Omega$ (see [17, 18]). As it turns out, there are extra geometrical relations between both graphs of $u^+$ and $u^-$. Indeed, let $K_u(x, y)$ denote the Gaussian curvature of the surface $S_L = graph(u)$ taken at the point $(x, y, u(x, y)) \in S_L$ and determined by the formula (see [20]):

$$K_u(x, y) = \frac{u_{xx}(x, y)u_{yy}(x, y) - u_{xy}^2(x, y)}{(1 + u_x^2(x, y) + u_y^2(x, y))^2}. \quad (11)$$

Interestingly, it is proved in [17, 18] that if both (7) and (10) hold then:

$$K_{u^+}(x, y) = -K_{u^-}(x, y). \quad (12)$$

Consequently, if $s_1 = (x, y, u^\pm(x, y)) \in S_L^\pm$ yields a hyperbolic point (having negative Gaussian curvature) then $s_2 = (x, y, u^\mp(x, y)) \in S_L^\mp$ renders an elliptic point (having positive Gaussian curvature). Thus for $p = (0, 0, -1)$ and $q = (q_1, q_2, q_3)$ the convexity/concavity ambiguity is automatically excluded. Noticeably, such ambiguity eventuates for single image shape from shading with $p = (0, 0, -1)$ (see e.g. [4, 8]). Finally, if $s_1$ is a parabolic point (with vanishing Gaussian curvature) then so is $s_2$.

This paper extends the above results with the following (see Section 2):

1. The necessary and sufficient condition (10) for testing the ambiguity in (7) is extended to arbitrary pairs of linearly independent light-source directions $p = (p_1, p_2, p_3)$ and $q = (q_1, q_2, q_3)$. Ultimately, such general case leads to the non-linear second-order PDE (see (18)), which inherently constitutes a challenging problem.
2. To alleviate the above difficulty, special configurations of two light-source directions (different from $p = (0, 0, -1)$ and $q = (q_1, q_2, q_3)$) are here admitted. Namely, the family of orthogonal unit vectors (parallel to OXZ-plane):

$$p_{\alpha} = (-\sqrt{1 - \alpha^2}, 0, -\alpha), \quad q_{\alpha} = (\alpha, 0, -\sqrt{1 - \alpha^2}) \quad (13)$$

(with $0 < \alpha < 1$) is considered. The corresponding analysis for solving (7) (with $p = p_{\alpha}$ and $q = q_{\alpha}$) addresses the uniqueness issue and establishes, in rare ambiguity cases, intrinsic geometrical inter-relations between multiple solutions in two light-source photometric stereo (see also Th. 1).

The theoretical Section 2 is also supplemented (see Section 3) by illustrative examples supported by *Mathematica* numerical and symbolic computation applied to the continuous Lambertian model with pixels represented as ideal points $(x, y) \in \Omega$ and image intensities $(E_1, E_2)$ simulated here from $u$ as left-hand sides in (7). Finally, the closing Section 4 summarizes the main thrust of this paper together with indicating its extensions and hints their possible solutions.

## 2 Orthogonal Illumination Directions in Vertical Plane

In the first part of this section (i.e in (i)) we consider an arbitrary configuration of linearly independent light-source directions $p$ and $q$ with $\|p\| = \|q\| = 1$. Under such general assumption, condition (10) (testing whether the second vector field from (8) over each $\Omega^{(j)}$ is also integrable) is subsequently extended to the non-linear PDE expressed by (18). Noticeably, the latter constitutes a difficult task for further theoretical analysis. In order to deal with the latter somehow, a special case of orthogonal vectors $p$ and $q$ contained in the OXZ-plane is here admitted (see (ii)). The analysis to follow complements already established results in [17], covering a different special choice of $p = (0, 0, -1)$ and $q = (q_1, q_2, q_3)$ with $q_1^2 + q_2^2 > 0$.

(i) Assume now that $E_1$ and $E_2$ introduced in (7) are generated by a genuine Lambertian surface $S_L = graph(u)$, where $u \in C^2$. Evidently, still by (7) both image intensities are also expressible in terms of $u_x$ and $u_y$. Using symbolic computation in *Mathematica* [21–23] (or alternatively see complicated proof in [17]) one arrives at:

$$\Lambda(x, y) = \frac{(u_x(p_2q_3 - p_3q_2) - u_y(p_1q_3 - p_3q_1) + p_2q_1 - p_1q_2)^2}{u_x^2 + u_y^2 + 1} = \frac{(\langle p \times q | n \rangle)^2}{u_x^2 + u_y^2 + 1}. \quad (14)$$

Combining the latter with (8) (over each $\Omega^{(j)}$ determined by $\Lambda > 0$) leads to:

$$u_x = \begin{cases} u_x, & \text{if } \operatorname{sgn}(\varepsilon)\operatorname{sgn}(\theta) > 0; \\ \frac{(a^2 - b^2 - c^2)u_x + 2acu_y + 2ab}{2abu_x + 2bcu_y + b^2 - a^2 - c^2}, & \text{if } \operatorname{sgn}(\varepsilon)\operatorname{sgn}(\theta) < 0 \end{cases} \quad (15)$$

and

$$u_y = \begin{cases} u_y, & \text{if } \operatorname{sgn}(\varepsilon)\operatorname{sgn}(\theta) > 0; \\ \frac{2acu_x + (c^2 - a^2 - b^2)u_y + 2bc}{2abu_x + 2bcu_y + b^2 - a^2 - c^2}, & \text{if } \operatorname{sgn}(\varepsilon)\operatorname{sgn}(\theta) < 0, \end{cases} \quad (16)$$
where

$$a = p_3q_2 - p_2q_3, \quad b = p_1q_2 - p_2q_1, \quad \text{and} \quad c = p_1q_3 - p_3q_1, \qquad (17)$$

and $\theta = u_x(p_2q_3 - p_3q_2) - u_y(p_1q_3 - p_3q_1) + p_2q_1 - p_1q_2$. Here the function $\varepsilon = \varepsilon(x, y)$ is everywhere constant taking values $\pm 1$. The first pair $(u_x, u_y)$ from (15) and (16) satisfies integrability condition (5) over each $\Omega^{(j)}$ as $u \in C^2$. On the other hand, the integrability of the second vector field from (15) and (16) eventuates, if and only if:

$$\left( \frac{(a^2 - b^2 - c^2)u_x + 2acu_y + 2ab}{2abu_x + 2bcu_y + b^2 - a^2 - c^2} \right)_y = \left( \frac{2acu_x + (c^2 - a^2 - b^2)u_y + 2bc}{2abu_x + 2bcu_y + b^2 - a^2 - c^2} \right)_x .$$

Upon resorting to the symbolic computation in *Mathematica*, the last equation is transformable into the following non-linear PDE:

$$c(a - bu_x(x, y)) u_{yy}(x, y) + (a^2 - c^2 + bcu_y(x, y) - abu_x(x, y)) u_{xy}(x, y) \\ +a(bu_y(x, y) - c) u_{xx}(x, y) = 0, \qquad (18)$$

which is to be satisfied by $u$. Evidently, the latter does not hold generically and therefore the resulting integrability condition (18) disambiguates almost always two-source photometric stereo modelled by (7) - see also Ex. 1. Still, in a pursue of solving a rare ambiguity in (7), one ought to deal with (18) which inevitably leads to a non-trivial task. Thus in the next step (ii) of this section, a tighter constraint imposed on illumination directions $p$ and $q$ is considered.

(ii) Suppose now that two light-source directions $p_\alpha$ and $q_\alpha$ are introduced according to (13). The resulting two image irradiance equations coincide with:

$$\frac{\alpha - \sqrt{1 - \alpha^2} u_x(x, y)}{\sqrt{u_x^2(x, y) + u_y(x, y)^2 + 1}} = E_1(x, y), \quad \frac{\alpha u_x(x, y) + \sqrt{1 - \alpha^2}}{\sqrt{u_x^2(x, y) + u_y(x, y)^2 + 1}} = E_2(x, y). \qquad (19)$$

Since $p_\alpha \perp q_\alpha$ (i.e. are orthogonal), the function $\Lambda$ defined in (9) simplifies into:

$$\Lambda(x, y) = 1 - E_1^2(x, y) - E_2^2(x, y) = \frac{u_y^2(x, y)}{1 + u_x^2(x, y) + u_y^2(x, y)}. \qquad (20)$$

Thus $\Omega_{\Lambda>0} = \{(x, y) \in \Omega : u_y(x, y) \neq 0\}$ and $\Omega_{\Lambda=0} = \{(x, y) \in \Omega : u_y(x, y) = 0\}$. Similarly to the special case of $p = (0, 0, -1)$ and $q = (q_1, q_2, q_3)$ discussed in [17, 18], often the triples $u, p_\alpha$ and $q_\alpha$ (see Ex. 1) yield $\Omega_{\Lambda>0} = \Omega^{(1)} \cup \Omega^{(2)}$, with $\Omega^{(j)}$ ($j = 1, 2$) standing for two disjoint sub-domains of $\Omega$ and $\Omega_{\Lambda=0} = \Gamma$ representing a smooth curve in $\Omega$. Here a bifurcation curve (along which solutions over $\Omega^{(1)}$ and $\Omega^{(2)}$ are glued) coincides with the curve (overlapping with $\Omega_{\Lambda=0}$):

$$\Gamma = \{(x, y) \in \Omega : u_y(x, y) = 0\}. \qquad (21)$$

Furthermore in (17), since $p_2 = q_2 = 0$ then $a = b = 0$ and since $q_1 = -p_3 = \alpha$ and $p_1 = q_3 = -\sqrt{1 - \alpha^2}$ then $c = 1$. Consequently, both formulae (15) and (16)
are reducible into:

$$ u_x = \begin{cases} u_x, & \text{if } \operatorname{sgn}(\varepsilon)\operatorname{sgn}(\theta) > 0; \\ u_x, & \text{if } \operatorname{sgn}(\varepsilon)\operatorname{sgn}(\theta) < 0; \end{cases} \qquad u_y = \begin{cases} u_y, & \text{if } \operatorname{sgn}(\varepsilon)\operatorname{sgn}(\theta) > 0; \\ -u_y, & \text{if } \operatorname{sgn}(\varepsilon)\operatorname{sgn}(\theta) < 0. \end{cases} \qquad (22) $$

Thus for the non-generic case of the second vector field ($u_x, -u_y$) also integrable (over each $\Omega^{(i)}$), the function $u$ should satisfy the following linear PDE:

$$ u_{xy}(x, y) = 0. \qquad (23) $$

Note that the last equation can be independently reached by substituting $a = b = 0$ and $c = 1$ into (18). The generic case of (23) not fulfilled is illustrated also in Ex. 1. Upon double integration (first over $x$ and then over $y$) of (23), the following representation for $u$ (if two vector fields from (22) are to be integrable over $\Omega^{(i)}$) holds:

$$ u(x, y) = \phi_1(x) + \psi_1(y), \qquad (24) $$

for some twice continuously differentiable functions $\phi_1$ and $\psi_1$ in a single variable. It is not difficult to show that $\phi_1(x) = u(x, 0) - c$, where $\psi_1(0) = c$. Similarly, $\psi(y) = u(0, y) - \phi_1(0) = u(0, y) - u(0, 0) + c$. Naturally, an analogous argument applies to the second solution $v \in C^2$ to (19) resulting in $v$ satisfying (23) and hence $v(x, y) = \phi_2(x) + \psi_2(y)$, where $\phi_2$ and $\psi_2$ are defined similarly to the introduction of $\phi_1$ and $\psi_1$. Furthermore, by (22), the second function $v \in C^2$ (over each $\Omega^{(j)}$) fulfills $(v_x, v_y) = (u_x, -u_y)$. Combining the latter with $u_x(x, y) = \phi'_1(x)$, $u_y(x, y) = \psi'_1(y)$, $v_x(x, y) = \phi'_2(x)$ and $v_y(x, y) = \psi'_2(y)$ yields "a conjugate-like" relation between $u$ and $v$:

$$ u(x, y) = \phi_1(x) + \psi_1(y), \qquad v(x, y) = \phi_1(x) - \psi_1(y) + c_1, \qquad (25) $$

with $c_1$ being a constant. Formula (25) determines specific analytic representations of the solutions to (19) over each $\Omega^{(j)}$ in a rare situation of the ambiguous two light-source photometric stereo. A straightforward verification shows that the geometrical constraint from (12) is also preserved for $p_\alpha$ and $q_\alpha$. Indeed, combining (11) with (25) leads to $K_u(x, y) = (\phi_1''(x)\psi_1''(y))(1 + (\phi_1'(x))^2 + (\psi_1'(y))^2)^{-2} = -K_v(x, y)$. Note that the condition (21) coupled with (25) implies that $\Gamma$ (or more general $\Omega_{\Lambda=0}$) represents a line $L = \{(x, y*) \in \Omega : \psi_1'(y*) = 0\}$ (or a collection of lines) parallel to the $OX$-axis. In addition, as $\nabla u = (u_x, u_y)$ and $\nabla v = (u_x, -u_y)$ visibly any critical point of $u$ is also a critical point of $v$. Moreover, if such point represents a local minimum (maximum, saddle) for $u$ then it is also a local maximum (minimum, saddle) for $v$. Note also that by (21) any critical point of $u$ (and thus of $v$) belongs to the set $\Omega_{\Lambda=0}$ and thus to the potential bifurcation curve $\Gamma$. The non-generic ambiguity case discussed above for $p_\alpha$ and $q_\alpha$ is illustrated in Ex. 2.

Evidently, as (18) or (23) are generically not fulfilled, there exists only one solution $u \in C^2$ to (7) (or (19)) over each $\Omega^{(j)}$. Upon gluing $u$ together, only one global solution $u \in C^2$ prevails over entire image $\Omega$. On the other hand, the rare
scenario of the existence of two solutions $u, v \in C^2$ over each $\Omega^{(j)}$ (i.e. satisfying (23)) leads to the possible bifurcations along $\Gamma$ rendering 0, 2 or 4 global solutions of class $C^2$ over entire image $\Omega$. Recall that we assume here $\Omega = \Omega^{(1)} \cup \Omega^{(2)} \cup \Gamma$. The other decomposition topologies of $\Omega = \Omega_{\Lambda>0} \cup \Omega_{\Lambda=0}$ are discussed in [17]. The detailed analysis justifying necessary and sufficient conditions to guarantee successful $C^k$ bifurcations ($k = 0, 1, 2$) (in case of $\Omega = \Omega^{(1)} \cup \Omega^{(2)} \cup \Gamma$) holding along $\Gamma$ exceeds the scope of this paper and therefore is here omitted.

Taking into account the above argument, the main theoretical contribution of this paper can be summarized into the following:

**Theorem 1.** Assume that $u \in C^2$ (which graph $S_L = graph(u)$ represents an illuminated genuine Lambertian surface) satisfies (19) (or (7)) with $p_\alpha$ and $q_\alpha$ determined by (13) (or with arbitrary linearly independent $p$ and $q$). Suppose, moreover that $\Lambda$ in (9) satisfies $\Lambda > 0$ over simply-connected $\Omega^{(j)}$ ($j = 1, 2$). In order that there exists only one more solution $v \in C^2$ to (19) over $\Omega^{(j)}$ it is necessary and sufficient for $u$ to satisfy (23) (or (18)). In addition, the ambiguous case for $p_\alpha$ and $q_\alpha$ yields $u(x, y) = \phi(x) + \psi(y)$ and $v(x, y) = \phi(x) - \psi(y)$, with $\phi$ and $\psi$ determined as in (24). Finally, the Gaussian curvatures of graphs of $u$ and $v$ at respective points $(x, y, u(x, y))$ and $(x, y, v(x, y))$ (see (11)) satisfy $K_u(x, y) = -K_v(x, y)$.

As already pointed out, a local ambiguity to (19) or (7) (over $\Omega_{\Lambda>0}$) can even be more proliferated to a global one (over entire image $\Omega$) due to possible bifurcations of $u$ and $v$ along $\Omega_{\Lambda=0}$. The matter gets more complicated if $\Omega_{\Lambda=0}$ forms an open subset of $\Omega$ (see [17, 18]). The latter occurs once $\langle p \times q | n \rangle = 0$ as implied by (14). Due to these intricacies, the respective discussion on bifurcation issue in two light-source photometric stereo is here left out.

## 3 Experiments

This section includes two examples illustrating the main results established in Section 2 (see also Th. 1). The experiments presented here are carried out with the aid of Mathematica symbolic computation. The corresponding pictures of images (with the respective intensities $E_1$ and $E_1$) are simulated synthetically upon admitting arbitrary or specific illumination directions and assuming a genuine $u$ as temporarily initially given.

**Example 1.** (a) Consider a Lambertian hemi-sphere $S_L^1 = graph(u_1)$ with $u_1 \in C^2(\hat{\Omega})$ defined as $u_1(x, y) = \sqrt{R^2 - x^2 - y^2}$, where $\hat{\Omega} = \{(x, y) \in \mathbb{R}^2 : x^2 + y^2 \le R^2\}$. For two linearly independent normalized light-source directions $p$ and $q$ the respective image irradiance equations read as:

$$ \frac{p_1 u_x + p_2 u_y - p_3}{\sqrt{u_x^2 + u_y^2 + 1}} = E_{11}(x, y) = \frac{-p_1 x - p_2 y - p_3 \sqrt{R^2 - x^2 - y^2}}{R}, \\ \frac{q_1 u_x + q_2 u_y - q_3}{\sqrt{u_x^2 + u_y^2 + 1}} = E_{21}(x, y) = \frac{-q_1 x - q_2 y - q_3 \sqrt{R^2 - x^2 - y^2}}{R}, \quad (26) $$
over $\Omega = \Omega_1 \cap \Omega_2 \cap \hat{\Omega}$, where $\Omega_i = \{(x, y) \in \mathbb{R}^2 : E_{i1}^1(x, y) \ge 0\}$ (for $i = 1, 2$). The negative values of $E_{i1}$ represent shadowed subareas of $\hat{\Omega}$. The bifurcation curve $\Gamma$ (see (14)) reduces into a planar quadratic determined by:

$$ (p_3q_2 - p_2q_3)x - (p_3q_1 - p_1q_3)y + (p_2q_1 - p_1q_2)\sqrt{R^2 - x^2 - y^2} = 0. \quad (27) $$

Furthermore, the condition (18) ascertaining the existence of exactly one solution $u_1 \in C^2(\Omega)$ (modulo its vertical shift) to (26) enforces $u_1$ to satisfy:

$$ (cx - ay) \left( ax + cy - b\sqrt{R^2 - x^2 - y^2} \right) = 0, \quad (28) $$

over $\Omega^{(j)}$, with $a, b$ and $c$ defined as in (17). Clearly, the equation (28) is not satisfied by $u_1$ for all $(x, y) \in \Omega^{(j)}$. This yields uniqueness in solving (26) (i.e. $u = u_1$) within $u \in C^2(\Omega)$. In particular, for $p_\alpha$ and $q_\alpha$ (see (13)), (17) combined with $a = b = 0$ and $c = 1$ reduce (27) into $y = 0$ rendering $\Gamma = \{(x, y) \in \hat{\Omega} : y = 0\}$, $\Omega^{(1)} = \{(x, y) \in \Omega : y < 0\}$ and $\Omega^{(2)} = \{(x, y) \in \Omega : y > 0\}$. The ambiguity condition (28) is transformed into $xy = 0$ merely fulfilled along both X- and Y-axes. Hence, again for arbitrary $p_\alpha$ and $q_\alpha$, uniqueness of $u_1$ prevails. This is expected since $p_\alpha$ and $q_\alpha$ represents a special case of general positions of $p$ and $q$ analyzed above.

(b) Let a Lambertian hill-like surface $S_L^2 = \text{graph}(u_2)$ with $u_2 \in C^2(\hat{\Omega})$ be defined according to $u_2(x, y) = (2(1 + x^2 + y^2))^{-1}$, over e.g. $\hat{\Omega} = \{(x, y) \in \mathbb{R}^2 : |x| \le 1 \text{ and } |y| \le 1\}$. The respective two image irradiance equations read as:

$$ \begin{aligned} \frac{p_1 u_x + p_2 u_y - p_3}{\sqrt{u_x^2 + u_y^2 + 1}} &= E_{12}(x, y) = \frac{-p_1 x - p_2 y - p_3 (x^2 + y^2 + 1)^2}{\sqrt{x^2 + y^2 + (x^2 + y^2 + 1)^4}}, \\ \frac{q_1 u_x + q_2 u_y - q_3}{\sqrt{u_x^2 + u_y^2 + 1}} &= E_{22}(x, y) = \frac{-q_1 x - q_2 y - q_3 (x^2 + y^2 + 1)^2}{\sqrt{x^2 + y^2 + (x^2 + y^2 + 1)^4}}, \end{aligned} \quad (29) $$

over unshadowed $\Omega = \Omega_1 \cap \Omega_2 \cap \hat{\Omega}$, where $\Omega_i = \{(x, y) \in \mathbb{R}^2 : E_{i2}^1(x, y) \ge 0\}$ (for $i = 1, 2$). The bifurcation curve $\Gamma$ from (14) is defined by the following equation:

$$ (p_2q_1 - p_1q_2) + \frac{(p_3q_2 - p_2q_3)x}{2(1 + x^2 + y^2)^{3/2}} + \frac{(p_3q_1 - p_1q_3)y}{2(1 + x^2 + y^2)^{3/2}} = 0. \quad (30) $$

On the other hand, the integrability condition (18) stipulates $u_2$ to satisfy:

$$ (cx - ay) \left( b + 6(ax + cy)\sqrt{1 + x^2 + y^2} \right) = 0, \quad (31) $$

over $\Omega^{(j)}$, with $a, b$ and $c$ introduced as in (17). Again, (31) is not fulfilled by $u_2$ for all $(x, y) \in \Omega^{(j)}$. Hence there exists a unique solution of class $C^2$ to (29) (i.e. $u = u_2$) over $\Omega$. In the special case of $p_\alpha$ and $q_\alpha$ (see (13)), (17) coupled with $a = b = 0$ and $c = 1$ transform (30) into $y = 0$. Hence again $\Gamma = \{(x, y) \in \hat{\Omega} : y = 0\}$, $\Omega^{(1)} = \{(x, y) \in \Omega : y < 0\}$ and $\Omega^{(2)} = \{(x, y) \in \Omega : y > 0\}$. Furthermore the ambiguity constraint (31) is reduced into $xy = 0$ which holds again only along both X- and Y-axes. Thus for $p_\alpha$ and $q_\alpha$ the uniqueness of $u_2$ to (29) eventuates (as also follows from the above general $(p, q)$-position analysis).
Example 2. Consider now a Lambertian paraboloid $S_L^3 = graph(u_3)$ with $u_3 \in C^2(\hat{\Omega})$ defined as $u_3(x, y) = (x^2+y^2)/2$ over $\hat{\Omega} = \{(x, y) \in \mathbb{R}^2 : |x| \le 1 \text{ and } |y| \le 1\}$. For a general position of $p$ and $q$, the integrability condition (18) reads here as $b(ay - cx) = 0$. Provided $b \ne 0$, such constraint never holds over any open subset of $\Omega$ and thus uniqueness of $u_3$ follows. Noticeably the bifurcation curve $\Gamma$ (a line) is determined here by $p_2q_1 - p_1q_2 + (p_2q_3 - p_3q_2)x + (p_3q_1 - p_1q_3)y = 0$.

For $p_\alpha$ and $q_\alpha$ from (13) the resulting two image irradiance equations are:

$$ \frac{\alpha - \sqrt{1 - \alpha^2}u_x(x, y)}{\sqrt{u_x^2(x, y) + u_y(x, y)^2 + 1}} = E_{13}(x, y) = \frac{\alpha - \sqrt{1 - \alpha^2}x}{\sqrt{x^2 + y^2 + 1}}, \\ \frac{\alpha u_x(x, y) + \sqrt{1 - \alpha^2}}{\sqrt{u_x^2(x, y) + u_y(x, y)^2 + 1}} = E_{23}(x, y) = \frac{\alpha x + \sqrt{1 - \alpha^2}}{\sqrt{x^2 + y^2 + 1}}, \quad (32) $$

over $\Omega_\alpha = \hat{\Omega}$ (for $\alpha = 1/\sqrt{2}$) and over $\Omega_\alpha = \{(x, y) \in \hat{\Omega}_\alpha : -1 \le x \le 0.204\}$ (for $\alpha = 1/5$). Fig. 1 shows images of $S_L^3$ with $\alpha = 1/\sqrt{2}$ over $\Omega_{1/\sqrt{2}} = \hat{\Omega} = [-1, 1] \times [-1, 1]$. The bifurcation curve $\Gamma_\alpha$ in (21) coincides here with the X-axis (i.e. here $y = 0$) trimmed either to $-1 \le x \le 1$ or to $-1 \le x \le 0.204$, for $\alpha = 1/\sqrt{2}$ or $\alpha = 1/5$, respectively. The remaining decomposition components of $\Omega_\alpha$ (along which $\Lambda > 0$) read as $\Omega_\alpha^{(1)} = \{(x, y) \in \Omega_\alpha : y > 0\}$ and $\Omega_\alpha^{(2)} = \{(x, y) \in \Omega_\alpha : y < 0\}$ - see Fig. 2. Furthermore, since by (17) the constant $b = 0$, the ambiguity condition $b(ay - cx) = 0$ is now satisfied. Hence, upon combining (25) with $\phi(x) = x^2/2$ and $\psi(y) = y^2/2$, the only one another $C^2$ solution to (32) equals to $u_4(x, y) = (x^2 - y^2)/2$ (modulo a vertical shift), over each $\Omega_\alpha^{(j)}$ (here $j = 1, 2$). Once $u_k$ are glued with itself ($k = 3, 4$) along X-axis, two $C^2$ class global solutions to (32) over entire image $\Omega_\alpha$ are defined (i.e. $u_3$ and $u_4$) - see Fig. 3 for $\alpha = 1/\sqrt{2}$ (with gluing curve). On the other hand the local solutions $u_3$ (or $u_4$) over $\Omega_\alpha^{(1)}$ cross-bifurcate along $\Gamma_{1/\sqrt{2}}$ with $u_4$ (or with $u_3$) over $\Omega_\alpha^{(2)}$ to yield next two only $C^1$ class solutions $u_{34}$ and $u_{43}$ to (32) over $\Omega_\alpha$. Indeed, a $C^2$ differentiability is excluded as $\lim_{(x,y) \to (x,0)} u_{3yy}(x, y) = 2 \ne -2 = \lim_{(x,y) \to (x,0)} u_{4yy}(x, y)$. The remaining two $C^1$ global solutions to (32) over $\Omega_{1/\sqrt{2}}$ (i.e. $u_{34}$ and $u_{43}$) are plotted in Fig. 4. The case for $\alpha = 1/5$ differs merely by different $\Omega_{1/5}$. Note that by (11) the respective Gaussian curvatures yield $K_{u_3}(x, y) = 4(1 + x^2 + y^2)^{-2} = -K_{u_4}(x, y)$ which is consistent with Th. 1. A unique critical point $(0, 0)$ of $u_3, u_4, u_{34}$ and $u_{43}$ belongs, as expected to $\Gamma$.

## 4 Conclusion

This paper extends the claims of [17, 18], where a special configuration of two light-source directions $p = (0, 0, -1)$ and $q = (q_1, q_2, q_3)$ (for $q_1^2 + q_2^2 > 0$) in photometric stereo is studied. In particular, a respective integrability condition expressed by the non-linear PDE (18) is derived here for a general configuration of linearly independent $p$ and $q$. Not unexpectedly, such PDE forms a difficult theoretical problem and consequently a new family of orthogonal illumination directions $p_\alpha$ and $q_\alpha$ is introduced in (13). The corresponding system of two image irradiance equation (19) is subsequently analyzed and the resulting ambiguity
<figure><img src="image_2.png" /><figcaption>Fig. 1. Two images of $S_L^3 = graph(u_3)$ illuminated along $p_{1/\sqrt{2}}$ and $q_{1/\sqrt{2}}$ directions.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Decomposition of $\Omega_\alpha = \Omega_\alpha^{(1)} \cup \Omega_\alpha^{(2)} \cup \Gamma_\alpha$ for (32) with a) $\alpha = 1/\sqrt{2}$, b) $\alpha = 1/5$.</figcaption></figure>

versus uniqueness question is addressed in Th. 1. As proved in this paper, similar ambiguity results are obtainable for $p_\alpha$ and $q_\alpha$ defined in (13) (see Th. 1 and Section 1) as compared to those already established in [17, 18] for $p = (0, 0, -1)$ and $q = (q_1, q_2, q_3)$. The experiments reported in Section 3 illustrate the main results from Section 2. A possible extension of this paper includes a relevant ambiguity analysis performed for the general positions of two linearly independent light-source directions. Inevitably, any pending argument would rely on characterizing the multiple solutions to the non-linear PDE defined in (18). Another worth investigation venue is to complete a global uniqueness-ambiguity analysis for two image photometric stereo over entire image $\Omega$ (see also [17, 18]). The latter should first cover the case of $p_\alpha$ and $q_\alpha$ defined in (13) and then should refer to an arbitrary selection of two linearly independent illumination directions. In particular, an extra attention should be paid here in derivation of analytical and numerical methods designed to localize possible bifurcation curve(s) $\Gamma$. Some relevant clues concerning this task can be found in [24–28, 30]. Finally, a similar analysis based on real Lambertian images, where pixelization and noise occur, forms another vital extension topic of this paper. Usually, handling noisy and digitized image data requires a more robust integration techniques dealing computationally with $u$ recovery phase determined by theoretical formula (4). Such analytic step is often substituted by the pertinent optimization task designed to compute numerically the closest discrete integrable vector field - for more see e.g. [4, 5, 29].
<figure><img src="image_4.png" /><figcaption>Fig. 3. Two $C^2$ class global solutions $u_3$ and $u_4$ to (32), over $\Omega_{1/\sqrt{2}} = [-1, 1] \times [-1, 1]$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Two $C^1$ class global solutions $u_{34}$ and $u_{43}$ to (32), over $\Omega_{1/\sqrt{2}} = [-1, 1] \times [-1, 1]$.</figcaption></figure>

## References

1. Hartley, R., Zisserman R.: Multiple View Geometry in Computer Vision. Cambridge University Press, Cambridge, UK, (2003).

2. Faugeras, O.: Three-Dimensional Computer Vision - A Geometric View Point. MIT Press, Cambridge, Massachusetts, London, England, (2001).

3. Trucco, E.: Introductory Techniques for 3-D Computer Vision. Prentice Hall, Englewood Cliffs, New York, (1998).

4. Horn, B.K.P.: Robot Vision. MIT Press, Cambridge, Massachusetts, London, England, (2001).

5. Horn, B.K.P., Brooks M.J.: Shape from Shading. MIT Press, Cambridge, Massachusetts, London, England, (1989).

6. Brooks, M.J., Chojnacki, W., Kozera, R.: Shading without shape. Quarterly of Applied Mathematics. 50(1), 27-38 (1992).

7. Brooks, M.J., Chojnacki, W., Kozera, R.: Circularly symmetrical eikonal equations and non-uniqueness in computer vision. Journal of Mathematical Analysis and Applications. 165 (1), 192-215 (1992).

8. Brooks, M.J., Chojnacki, W., Kozera, R.: Impossible and ambiguous shading patterns. International Journal of Computer Vision. 7(1), 119-126 (1992).

9. Bruss, A.R.: The eikonal equation: some results applicable to computer vision. Journal of Mathematical Physics. 5(23), 890-896 (1982).

10. Kozera, R.: On complete integrals and uniqueness in shape from shading. Applied Mathematics and Computation. 73(1), 1-37 (1995).
11. Deift, P., Sylvester, J.: Some remarks on shape-from-shading in computer vision. Journal of Mathematical Analysis and Applications. 1(84), 235-248 (1991).

12. Oliensis, J. Uniqueness in shape from shading. International Journal of Computer Vision. 6(2), 75-104 (1991).

13. Chojnacki, W., Brooks, M.J.: A direct computation of shape from shading. Proceedings of 12-th International Conference on Artificial Intelligence and Pattern Recognition. Haifa, Israel, 114-119 (1994).

14. Kimmel, R., Bruckstein, A.: Tracking level sets by level sets: a method of solving the shape from shading problem. Computer Vision, Graphics and Image Understanding. 62, 47-58 (1995).

15. Kozera, R.: Uniqueness in shape from shading revisited. International Journal of Mathematical Imaging and Vision. 7, 123-138 (1997).

16. Woodham, R.J.: Photometric stereo: a reflectance map technique for determining surface orientation from multiple images. Optical Engineering. 19(1), 139-144 (1980).

17. Kozera, R.: Existence and uniqueness in photometric stereo. Applied Mathematics and Computation. 44(1), 1-104 (1991).

18. Kozera, R.: On shape recovery from two shading patterns. International Journal of Pattern Recognition and Artificial Intelligence. 6(4), 673-698 (1992).

19. Onn, R., Bruckstein, A.M.: Integrability disambiguates surface recovery in two-image photometric stereo. International Journal of Computer Vision. 5(1), 105-113 (1990).

20. Do Carmo, M.P.: Differential Geometry of Curves and Surfaces. Prentice-Hall, Englewood Cliffs, New York, (1986).

21. Wolfram Mathematica 9, Documentation Center, URL link: reference.wolfram.com/mathematica/guide/Mathematica.html.

22. Budzko, D.A., Prokopenya, A.N.: Symbolic-numerical methods for searching equilibrium states in a restricted four-body problem. Programming and Computer Software 39(2), 74-80 (2013).

23. Budzko, D.A., Prokopenya, A.N.: Symbolic-numeric analysis of the equilibrium solutions in the restricted four-body problem. Programming and Computer Software 36(2), 68-74 (2010).

24. Kozera, R.: Curve modeling via interpolation based on multidimensional reduced data. Studia Informatica 25(4B-61), 1-140 (2004).

25. de Boor, C.: A Practical Guide to Spline. Springer-Verlag, New York Heidelberg, Berlin, (1985).

26. Kozera, R., Noakes, L., Klette, R.: External versus internal parameterization for length of curves with nonuniform samplings. Lecture Notes in Computer Science 2615, Springer, Geometry, Morphology and Computational Imaging, 403-418 (2003).

27. Kvasov, B.I.: Methods of Shape-Preserving Spline Approximation. World Scientific, Singapore, (2000).

28. Kozera, R., Noakes, L.: Piecewise-quadratics and exponential parameterization for reduced data. Applied Mathematics and Computation 221, 620-638 (2013).

29. Noakes, L., Kozera, R.: Non-linearities and noise reduction in 3-source photometric stereo. Journal of Mathematical Imaging and Vision, 18(2), 119-127 (2003).

30. Atamanyuk, I.P., Kondratenko, V.Y., Kozlov, O.V., Kondratenko, Y.P.: The algorithm for optimal polynomial extrapolation of random processes. Lecture Notes in Business Information Processing 115, Springer, Berlin-Heidelberg, 78-87 (2012).