Solving Uncalibrated Photometric Stereo Using Fewer Images by Jointly Optimizing Low-rank Matrix Completion and Integrability

TL;DR

This paper presents an integrated optimization method for uncalibrated photometric stereo that effectively reconstructs 3D shapes from as few as 4-6 images by jointly enforcing low-rank and integrability constraints.

Contribution

It introduces a novel unified optimization framework using ADMM that incorporates rank, integrability, and missing data constraints for better 3D reconstruction.

Findings

Improved 3D reconstruction quality with fewer images (4-6)

Effective handling of missing data in photometric stereo

Significant performance gains over prior methods in experiments

Abstract

We introduce a new, integrated approach to uncalibrated photometric stereo. We perform 3D reconstruction of Lambertian objects using multiple images produced by unknown, directional light sources. We show how to formulate a single optimization that includes rank and integrability constraints, allowing also for missing data. We then solve this optimization using the Alternate Direction Method of Multipliers (ADMM). We conduct extensive experimental evaluation on real and synthetic data sets. Our integrated approach is particularly valuable when performing photometric stereo using as few as 4-6 images, since the integrability constraint is capable of improving estimation of the linear subspace of possible solutions. We show good improvements over prior work in these cases.

Equations

$$\underset{\hat{M}}{\min }\Vert \hat{M}-M{\Vert }_F^{2}s.t.\mathrm{rank}(\hat{M})=3.$$

$$\underset{\hat{M}}{\min }\Vert \hat{M}-M{\Vert }_F^2$$

$$s.t.$$

$$M_{ij}=\max (0,{\rho }_j{l}_i^T{\hat{n}}_j)$$

$$M=LS.$$

$$S=N\Lambda ,$$

$$X=\left[\begin{array}{cc}{X}^I& X^N\\ X^L& X^M\end{array}\right]=\left[\begin{array}{cc}I& N\\ L& M{\Lambda }^{-1}\end{array}\right],$$

$$X=\left[\begin{array}{c}{A}^{-1}\\ L{A}^{-1}\end{array}\right]\left[\begin{array}{cc}A& AN\end{array}\right].$$

$$X^N={\left[\begin{array}{ccc}{D}_x\mathbf{z},& D_y\mathbf{z},& -1\end{array}\right]}^T,$$

$$-1\le {\lambda }_j\le 0$$

$$X^I=I_{3\times 3}.$$

$$f_{data}(X,\Lambda )=\frac{1}{2}||W\odot (M-X^M\Lambda )|{|}_F^2,$$

$$\underset{X,\Lambda ,\mathbf{z}}{\min }$$

$$s.t.$$

$$f_{tnn}(X)=||X|{|}_{\ast }-\sum _{k=1}^3{\sigma }_k(X),$$

$$\underset{X,\Lambda ,\mathbf{z}}{\min }$$

$$s.t.$$

$$f_{maj}^{(k)}(X)=\Vert X{\Vert }_{\ast }-\mathrm{trace}(U_3^{T}X{V}_3).$$

$$\underset{X,\Lambda ,\mathbf{z}}{\min }$$

$$s.t.$$

$$\underset{\Gamma }{\max }\underset{X,\Lambda ,\mathbf{z},Y}{\min }$$

$$c\left(||Y|{|}_{\ast }-\mathrm{trace}(U_3^{T}Y{V}_3)\right)+\frac{\tau }{2}||Y-X+\Gamma |{|}_F^2$$

$$s.t.X^I=I_{3\times 3},$$

$$X^{L(k+1)}$$

$$(X^{N(k+1)},{\mathbf{z}}^{(k+1)})$$

$$s.t.X^N$$

$$(X^{M(k+1)},{\Lambda }^{(k+1)})=$$

$$+\frac{\tau }{2}\Vert Y^{M(k)}-X^M+{\Gamma }^{M(k)}{\Vert }_F^2$$

$$s.t.-1\le {\lambda }_j\le 0\forall j.$$

$$X_{ij}^{M(k+1)}=Y_{ij}^{M(k)}+{\Gamma }_{ij}^{M(k)}.$$

$$(X_j^{M(k+1)},{\lambda }_j^{(k+1)})=$$

$$+\frac{\tau }{2}\Vert Y_j^{M(k)}-X_j^M+{\Gamma }_j^{M(k)}{\Vert }_2^2$$

$$s.t.-1\le {\lambda }_j\le 0.$$

$$\widetilde{X_j}$$

$$=\frac{{\lambda }_j\widetilde{M_j}+\tau {\tilde{A}}_j^{M(k)}}{{\lambda }_j^2+\tau }.$$

$${\lambda }_j$$

Full text

1 Introduction

Uncalibrated photometric stereo (UPS) is the problem of recovering the 3D shape of an object and associated lighting conditions, given images taken with varying, unknown illumination. In this work we replace the existing pipeline for solving UPS with an integrated approach. This paper, like much prior work  [15, 4, 31, 29, 2, 11, 22], focuses on Lambertian objects illuminated by a single distant point light source in each image. Existing methods, pioneered by [15], formulate UPS as the problem of finding a low-rank factorization of the measurements. Specifically, given $m$ images each with $p$ pixels, let $M$ denote the $m\times p$ matrix containing the pixel intensities. These methods optimize

[TABLE]

This problem can be solved by SVD, from which we produce a family of solutions, each consisting of a set of light sources, albedos, and surface normals. These solutions are related by a $3\times 3$ ambiguity matrix. The surface normals provided by SVD are in general inconsistent with the partial derivatives of the surface (i.e. they are not integrable). Consequently, existing methods apply an additional sequence of steps aimed at reducing the ambiguity and fitting a surface to the recovered normals.

In this paper we propose instead to optimize:

[TABLE]

Eq. (1) optimizes over rank 3 matrices, which can represent sets of images produced by any set of surface normals. In contrast, in (2) we optimize over only those rank 3 matrices that correspond to integrable surfaces.

Intuitively, a single optimization over all constraints will have a better global optimum than a sequence of optimizations in which constraints are used one at a time to increasingly narrow the solution (see illustration in Figure 1). A similar intuition has motivated the use of bundle adjustment [14] as the dominant approach to large scale structure-from-motion. Specifically in UPS the measurement matrix may contain many errors due to shadows and specular effects. Therefore, while in theory UPS can be solved with as few as three images, SVD can properly handle these modeling errors only when many images are supplied. Indeed, current methods [2, 11] typically use 10 or more images. With fewer images SVD results tend to provide noisy solutions. Our method incorporates integrability into this estimation, providing valuable additional constraints that can reduce this noise. Our experiments indicate that our method can produce reasonable reconstructions with as few as 4 images and good reconstruction with 6 images, significantly improving over state-of-the-art methods with these few images.

For our approach we optimize a cost function based on (2) over the surface, lighting, normals, and (restored) error-free observations. The cost ensures that normals and lighting are consistent with the measurements, which must have low rank. We use constraints that ensure integrability. This is somewhat tricky because rank constraints apply to the measurements while integrability constraints apply to the normals. We show that by constructing a rank 3 matrix that contains normals, measurements, and lighting, we can impose the rank and integrability constraints together. Specifically, we use a truncated nuclear norm approach [16] to enforce the rank constraint, while integrability is represented by linear equalities. This leads to a single non-convex problem that we solve using a series of Alternate Direction Method of Multipliers (ADMM) operations [5, 13].

Our formulation allows us to easily account for missing data in the measurement matrix. This commonly occurs when pixels are dark due to shadows, or saturated due to specularities. In some of the prior approaches, this can be solved with a preprocessing step, which may lead to a pipeline with yet another optimization [29]. We handle missing data using matrix completion based on the rank constraint. We initialize our optimization using prior approaches, in much the same way that bundle adjustment is initialized using simpler, but non-optimal algorithms [14].

2 Background and Previous Work

In this section we introduce in detail the problem of uncalibrated photometric stereo for Lambertian objects and review past work. We assume that we view an object in multiple images from a fixed viewpoint. In each image the object is illuminated by a single, distant point light source. We represent lighting in image $i$ with $l_{i}\in R^{3}$, in which the direction of $l_{i}$ represents the direction to the lighting, and $\|l_{i}\|$ represents its magnitude. We represent the object using a set of surface normals $\hat{n}_{j}\in R^{3}$, and albedos $\rho_{j}\in R$ for each pixel. We then obtain images with the equation:

[TABLE]

where $M_{ij}$ represents the $j$-th pixel of the $i$-th image. We define the surface normal $\hat{n}_{j}=\frac{n_{j}}{\|n_{j}\|}$, $n_{j}=(-z_{x},-z_{y},1)^{T}$, where $z_{x}$ and $z_{y}$ represent partial derivatives of the surface $z(x,y)$ at pixel $j$. Negative values of $\rho_{j}l_{i}^{T}\hat{n}_{j}$ are set to 0; these appear as attached shadows.

We now describe the creation of all images using matrix operations. We define $S$ to be a $3\times p$ matrix in which column $j$ contains $\rho_{j}\hat{n}_{j}$. Given $m$ images, we can stack the light into the matrix $L$ of dimension $m\times 3$, where each row denotes one light per image. We concatenate all the images to form an observation matrix $M$ of dimension $m\times p$, where $p$ is the number of pixels. Now, in the absence of shadows, we can write the equation of UPS as:

[TABLE]

Classical work on photometric stereo (e.g. [28], see a recent review in [1]) has assumed that known lighting is obtained by careful calibration. With $L$ known, (4) can be solved as a linear least squares problem. A more general and challenging case is unconstrained photometric stereo, in which the $L$ is unknown. A common approach, which we use as a baseline algorithm, follows the steps in Algorithm 1.

Belheumer et al. [4] showed that in UPS the integrable set of surface normals can only be recovered up to a Generalized Bas-Relief transformation (GBR). A number of recent papers have concentrated on methods of solving the GBR ambiguity. Researchers have used priors on the albedo distribution [2], reflectance extrema [11], grouping based on image appearance and color [24], inter-reflections [9], isotropy and symmetries [27], and specularity [10] as constraints while solving for the GBR. All of these methods have first used the above mentioned baseline described in Algorithm 1 to obtain a solution up to the GBR.

Recent works have explored a variety of other research directions in photometric stereo. Mecca et al. [18] proposed an integrated, PDE based approach to calibrated photometric stereo that uses a mere two images under perspective projection. It is not clear how to extend this to uncalibrated photometric stereo. Basri et al. [3] extended the baseline to handle multiple light sources in each image using a spherical harmonics formulation. Chandraker et al. [8] proposed a method to handle attached and cast shadows in the case of multiple light sources per image. In [25] the authors determine the visibility subspace for a set of images to remove the cast and attached shadows for performing UPS. Various works have addressed non-Lambertian materials (e.g., Georghiades et al. [12] and Okabe et al. [20]).

In the context of Lambertian UPS, Georghiades et al. [12] proposed to remove shadows and specularities and recover the missing pixel values using matrix completion algorithms, e.g., using the damped Wiberg [21] or Cabral’s algorithm [6]. Wu et al. [29] proposed a Robust PCA formulation as preprocessing for calibrated photometric stereo. Their approach seeks a low-rank (not necessarily rank 3) approximation to $M$ while removing outlier pixels (corresponding to shadows and specularities). Oh et al. [19] applied Robust PCA in the context of calibrated photometric Stereo, replacing the Nuclear Norm with a Truncated Nuclear Norm (TNN) regularizer [16]. In [11], Favaro et al. have used Robust PCA as preprocessing to the baseline algorithm for UPS.

3 Our Approach

In this section we introduce our integrated formulation that enforces integrability of surface normals in solving the uncalibrated photometric stereo problem. We recall from (4) that the measurement matrix $M$ can be factored into $M=LS$. To access the derivatives of $z(x,y)$ we write $S$ as a product

[TABLE]

where $N$ is a $3\times p$ matrix whose $j$’th column is $n_{j}=(-z_{x},-z_{y},1)^{T}$ and $\Lambda=\mathrm{diag}(\lambda_{i},\lambda_{2},...,\lambda_{p})$ with $\lambda_{j}=-\rho_{j}/\|n_{j}\|$. We next define the matrix:

[TABLE]

where $X$ is $(3+m)\times(3+p)$. The matrices $X$, $\Lambda,$ and the depth values ($z(x,y)$) form the unknowns in our optimization. Note that, because $LN=M\Lambda^{-1}$, the following holds for any $3\times 3$ non-degenerate matrix $A$

[TABLE]

This shows that $X$ is rank 3. The matrix $A$ represents a linear ambiguity. However, forcing the normals in $N$ to be integrable will reduce this ambiguity to the GBR.

To force integrability we denote by $\mathbf{z}=(z_{1},...,z_{p})^{T}$ the vector of unknown depth values and require

[TABLE]

where $D_{x},D_{y}$ denote respectively the $x$- and $y$-derivative operators and $\mathbf{1}$ denotes the vector of all 1’s.

Additional constraints are obtained by noticing that, because $0\leq\rho_{j}\leq 1$ and $\|n_{j}\|\geq 1$,

[TABLE]

and

[TABLE]

We are now ready to define our optimization function. Let $W$ be a binary, $m\times p$ matrix so that $W_{ij}=0$ if $M_{ij}$ is missing and $W_{ij}=1$ otherwise, and let

[TABLE]

where $\odot$ denotes element-wise multiplication. Then (2) can be written as

[TABLE]

Handling the rank-3 constraint: Enforcing the non-convex constraint $\mathrm{rank}(X)=3$ can be challenging. In the context of matrix completion a recent paper [16] proposed using the Truncated Nuclear Norm (TNN) regularization term:

[TABLE]

where $||X||_{*}$ denotes the nuclear norm of $X$ and $\sigma_{k}(X)$ is the $k$-th largest singular value of $X$. Clearly, $f_{tnn}(X)=0$ if and only if $\mathrm{rank}(X)\leq 3$. We use $f_{tnn}$ as a regularizer and solve

[TABLE]

where $c$ is a preset scalar.

4 Optimization using ADMM

In this section we introduce a method for solving (3). This is a challenging problem because both $f_{data}$ and $f_{tnn}$ are non-convex. Specifically, $f_{data}$ (11) is bilinear in $X$ and $\Lambda$, while $f_{tnn}$ (13) is a difference between two convex functions. Our solution is based on a nested iteration in which the outer loop uses majorization to decrease $f_{tnn}$ whereas the inner loop uses the scaled ADMM with alternation to decrease $f_{data}$.

Outer loop: Following [16] at each iteration of the outer loop we replace $f_{tnn}(X)$ with a majorizer. Specifically, at iteration $k$ let $X^{(k)}=U\Sigma V^{T}$ be the singular value decomposition of $X^{(k)}$, and let $U_{3}$ (and $V_{3}$) be the matrices containing the left (right) singular vectors corresponding to the three largest singular values of $X^{(k)}$. $U_{3}$ and $V_{3}$ are determined in the outer loop and are held constant throughout the inner loop. We then define

[TABLE]

It was shown in [16] that $f_{maj}^{(k)}(X)\geq f_{tnn}(X)$ for all $X$ and that $f_{maj}^{(k)}(X^{(k)})=f_{tnn}(X^{(k)})$, and so decreasing $f_{maj}$ leads to decreasing $f_{tnn}$.

Inner loop: In the inner loop we seek to minimize

[TABLE]

We use scaled ADMM, a variant of the augmented Lagrangian method that splits the objective function and aims to solve the different subproblems separately. We maintain a second copy of $X$, which we denote by $Y$ and form the augmented Lagrangian of (16) as follows

[TABLE]

where $||Y-X+\Gamma||^{2}_{F}$, denotes the Lagrangian penalty; $\tau$ is a constant, and $\Gamma$ is a matrix of Lagrange multipliers the same size as $X$ that is updated by the ADMM steps [5, 13]. We next describe the ADMM steps (applied iteratively).

Step 1**: Solving for $(X,\Lambda,\mathbf{z})$**.

In each iteration, $k$, we solve the following sub-problems:

Optimize w.r.t. $X^{I}$: $X^{I\,(k+1)}=I_{3\times 3}$. 2. 2.

Optimize w.r.t. $X^{L}$:

[TABLE] 3. 3.

Optimize w.r.t. $X^{N}$ and $\mathbf{z}$:

[TABLE]

The problem is solved by setting the third row of $X^{N\,(k+1)}$ to $-\mathbf{1}$ and by substituting $D_{x}\mathbf{z}$ and $D_{y}\mathbf{z}$ for the first two rows of $X^{N}$ in the objective, obtaining linear least squares equations in $\mathbf{z}$ that can be solved directly. 4. 4.

Optimize w.r.t. $X^{M}$ and $\Lambda$:

[TABLE]

We will separate this into the known and unknown pixels based on $W$. For an unknown pixel $j$ in frame $i$ ($W_{ij}=0$) the first term vanishes and the minimization only determines the respective entry of $X^{M}$ so that:

[TABLE]

For the known pixels, since $\Lambda$ is diagonal we can write these equations separately for each column $j$ (corresponding to the $j$-th pixel):

[TABLE]

The problem (21) is non-convex. We will solve it with alternate optimization. $X^{M}$ and $\Lambda$ are updated by the following steps until convergence.

$\mathbf{X^{M}}$ : Let $\tilde{M_{j}}=W_{j}\odot M_{j}$ , $\tilde{X_{j}}=W_{j}\odot X^{M}_{j}$ and $\tilde{A}^{M\,(k)}_{j}=W_{j}\odot(Y_{j}^{M\,(k)}+\Gamma_{j}^{M\,(k)})$. Then,

[TABLE]

$\mathbf{\Lambda}$:

[TABLE]

Step 2**: Solving for $Y$.** Solving for $Y$ requires a solution to

[TABLE]

Below we show that this problem can be solved in closed form by applying the shrinkage operator, obtaining

[TABLE]

where the shrinkage operator $D_{t}(.)$ is defined as follows. For a scalar $s$ we define $D_{t}(s)=\mathrm{sign}(s)\times\max(|s|-t,0)$. For a diagonal matrix $S=\mathrm{diag}(s_{1},s_{2},...)$ with non-negative entries we define $D_{t}(S)=\mathrm{diag}(D_{t}(s_{1}),D_{t}(s_{2}),...)$. Finally, for a general matrix $\Upsilon$, let $\Upsilon=\tilde{U}S\tilde{V}^{T}$ be its singular value decomposition, then $D_{t}(\Upsilon)=\tilde{U}D_{t}(S)\tilde{V}^{T}$.

To derive (25), we rewrite (4) as:

[TABLE]

where $T=\mathrm{trace}(V_{3}U_{3}^{T}(X^{(k+1)}-\Gamma^{(k)}))+\dfrac{c}{2\tau}\|U_{3}V_{3}^{T}\|_{F}^{2}$ is independent of $Y$. Equation (26) is of the general form $\underset{Y}{\min}\|Y\|_{*}+\dfrac{1}{2t}\|Y-C\|^{2}_{F}$, for which the solution is $D_{t}(C)$, as is shown in [7], implying (25).

Step 3**: Update of $\Gamma$.** The matrix $\Gamma$ contains Lagrange multipliers that are used in the saddle-point formulation (17) to enforce the equality constraint $X=Y$. The following update is a gradient ascent step that acts to maximize the augmented Lagrangian (17) for $\Gamma.$ For details, see [5, 13].

[TABLE]

The entire optimization process is listed in Algorithm 2. We will make the code available.

5 Experimental Results

In this section we evaluate and compare the performance of our algorithm with two versions of the baseline algorithm, in both real world and synthetic examples. We compare the following methods:

Baseline: Algorithm 1 described in Section 2. This method is used in  [2, 11, 24, 9, 27, 10].

RPCA: Images are preprocessed using Robust PCA [29], parameters are chosen as suggested by [11]. Then we apply the baseline algorithm to the obtained matrix. This method is used in [11]. RPCA solves a sparse low rank optimization to detect shadows and other non-Lambertian effects. The method uses $L_{1}$ regularization to identify outlier pixels, even when they do not result in intensities near 0 or 1.

Our(NC): Our proposed formulation as described in Section 4 using $W=1$, i.e., no completion. This allows comparison to Baseline, which also does not perform matrix completion.

Our(MC): Our proposed formulation as in Section 4 with $w_{ij}\in\{0,1\}$, allowing for matrix completion. In both versions of our algorithm we use $c=1$ and use RPCA to initialize optimization. We identify missing pixels as those with normalized intensity outside the range of $(0.02,0.98)$.

All the tested methods solve for the surface only up to a GBR ambiguity. To compare the results with ground truth, we find the GBR that optimizes the fit to ground truth, and measure the residual error.

In the presence of a large number of images with noise and non-Lambertian effects, we expect the sequential pipeline of Baseline and RPCA, involving SVD, to produce accurate solutions, because the problem solved by SVD is heavily overconstrained. In the presence of fewer images, our integrated method will be able to produce a more accurate decomposition by using both rank and integrability constraints to find the right linear subspace. Thus we expect our integrated approach to improve over the Baseline and RPCA as we reduce the number of images. In the following sub-section we will show results with synthetic and real world data that supports our claim.

5.1 Experiments on Synthetic Data

We use five real objects (“cat”, “owl”, “rock”, “horse”, “buddha”) to produce synthetic images, their shape is obtained by applying calibrated photometric stereo to a publicly available dataset [17]. We use the normals and albedos from these objects to generate images. Each image is generated by a randomly selected light source which lies at 30 degrees of the viewing direction on average. All images are of size 512 $\times$ 340 with objects occupying 29-72K pixels. A segmentation mask is also supplied. To show the variation of performance with the number of images $N_{I}$, we use sets of 4, 6, 8, 10, 15, 20, 25 and 30 images respectively. We add Gaussian noise with standard deviation ranging from 1% to 7% (in steps of 2%) of the maximum intensity. For each choice of noise, we run 5 different trials with random noise and lighting to generate the synthetic images. Thus we have 5 objects, 4 levels of noise and 5 random simulations, making a total of 100 experiments for each of the 8 different sets of 4, 6, 8, 10, 15, 20, 25 and 30 images. As a measure of performance, we calculate the error in the reconstructed depth map. Let the ground truth surface be $Z_{T}$ and the reconstructed surface be $Z_{rec}$. We measure error in depth as $Z_{err}=100\times\frac{||Z_{T}-Z_{rec}||}{||Z_{T}||}$. To compare two algorithms (say, algorithm A vs. algorithm B), we define the following two terms :

Relative Improvement (in $\%$) : Denote $e^{a}_{k}$ and $e^{b}_{k}$ as the depth error for each trial $k$ by using algorithm A and B respectively. The Relative Improvement of algorithm B over A is the average of $\frac{(e^{b}_{k}-e^{a}_{k})}{e^{b}_{k}}$ over all trials $K$ for each choice of $N_{I}$ expressed in percentage.

Percent of Improved Trials : This denotes the number of trials in which algorithm B improves over A. In terms of notation introduced previously, this is $\frac{1}{K}\sum_{k=1}^{K}\mathbb{I}(e^{\mathrm{a}}_{k}<e^{\mathrm{b}}_{k})$, where $\mathbb{I}(.)$ is in indicator variable and $K$ is the total number of trials for each choice of $N_{I}$. The measure is expressed in percentage.

In Figure 2(a) we compare performance of Our(MC) with Baseline and RPCA, on synthetic data in the presence of Gaussian noise. We initialize our methods with RPCA. We observe that as the number of images decreases, our method improves compared to Baseline and RPCA. With simple Gaussian noise RPCA doesn’t produce additional advantages as there are no outliers.

In Figure 2(b) we compare the performance of our methods on synthetic data with Gaussian noise and with specularities generated by the Phong reflectance model [26, 23]. Mathematically each image $M_{i}$ can be represented as :

[TABLE]

where $V$ is the viewing direction and $R$ denotes the directions of perfect reflection for incoming light $L_{i}$ for each pixel $j$. Larger $\alpha$ produces sharper specularities, while larger $k_{s}$ causes more light to be reflected as specularity. We use $k_{s}=0.2$ and $\alpha=10$. We observe that the advantage of Our(MC) degrades as the number of images increases, as expected. This experiment shows that even though our method is designed specifically for Lambertian objects it can tolerate a certain amount of model irregularities such as specularity. With 4 images our method beats RPCA in 85% of the all trials with a relative improvement of 22.12%.

In Figure 3 we compare Our(MC) with Baseline and RPCA with variation of noise for different subsets of images (4,6,10 and 15). We can conclude that our method is robust to noise and its advantages do not degrade with an increase in noise.

5.2 Experiments on Real World Data

To test our approach on real data, we used the two publicly available data sets [17] and [30] consisting of 5 and 7 objects respectively. The datasets provide calibrated lighting, which we use to perform calibrated photometric stereo. The obtained depth map, albedo, and surface normals are considered as ground-truth for photometric stereo with unknown lighting similar to [2]. To show the variation of performance with the number of images, we select subset of 4, 6, 8 and 10 images for each object. We perform 10 random selections of subset of images for each of the 12 objects. Thus we have 120 experiments for every subset of images.

In Figure 4 we compare the performance of our methods, Our(MC) and Our(NC), with Baseline and RPCA with variation in the number of images. We see that for fewer images our methods outperform Baseline and RPCA by a significant amount and are comparable to RPCA for more images. For 4 images Our(MC) outperforms Baseline in 84.9% cases with a relative improvement of 30.6% and outperforms RPCA in 81.4% cases with a relative improvement of 12%. However for 10 images we beat Baseline in 75% cases with a relative improvement of 10.7% and beat RPCA in only 47.3% cases with a relative improvement of -7.2%.

Figure 5 shows the average reconstruction error obtained by Our(MC), RPCA and Baseline on 12 real-world objects over 10 random simulations. We observe that Our(MC) outperforms RPCA on 11 out of 12 objects for 4 images and 10 out of 12 objects for 6 images (and is comparable in 1). With 10 images the average reconstruction error using Our(MC) over all objects and all trials is 4.6%. This increases to 8.1% with four images, and is only 5.4% with six images. This shows that we have reasonable reconstruction with 4 images and good reconstruction with as few as 6 images.

In Figure 6 we compare the error in surface reconstruction between Baseline, RPCA, and Our(MC) on some of our real world examples. Figure 7 shows two views of surfaces reconstructed using the Our(MC) algorithm using 4 images, showing reasonable surface reconstruction. These results suggest that our joint approach to enforcing rank and integrability constraints can significantly improve the performance of photometric stereo in the presence of a few images.

In general, we see that incorporating matrix completion into our formulation results in a slight improvement, with Our(MC) somewhat outperforming Our(NC). This indicates that the improvement of our method compared to RPCA or Baseline is mostly due to the joint optimization formulation and not due to matrix completion. We further note that RPCA seems to significantly improve over Baseline. RPCA is able to identify outliers and use that extra information for better recovery. This also suggests that the robust error function used by RPCA is important. However our integrated approach, which does not have a robust cost function like RPCA, still outperforms RPCA for 4 and 6 images and is almost equal for 8 or 10 images. This shows that an integrated approach is very useful for a small number of images and provides similar gain compared to RPCA for more images. It would be an interesting topic of future work to amend the cost function of Our(MC) to include RPCA’s robust handling of error, to see if this further improves its performance.

For an image of size $512\times 340$ with an object occupying an area of 30K pixels, our algorithm takes 20 minutes on a 2.7 GHz Intel Core i5 machine.

6 Conclusion and Future Work

In this paper we have introduced a new low-rank constrained optimization method for solving uncalibrated photometric stereo using fewer images. The key to this approach is to combine rank and integrability constraints in a single optimization problem. This relies on a novel formulation that exposes both depth and surface normals to the optimization, linking them with an integrability constraint. We then show how to perform this optimization using a truncated nuclear norm and ADMM. Our joint formulation produces better solutions, compared to other methods that use SVD, for fewer images. We have shown promising results compared to baseline approaches using both real and synthetic examples. We also observe that our method can handle certain degrees of model irregularities as it has outperformed RPCA in synthetic examples with specularities generated using the Phong model.

In the future, it will be interesting to apply the idea of Robust PCA to our formulation. We would also like to extend this work to handle more general lighting configurations, e.g., using spherical harmonic approximations to lighting.