Multispectral Photometric Stereo for Spatially-Varying Spectral Reflectances

Abstract

Multispectral photometric stereo (MPS) aims at recovering the surface normal of a scene measured under multiple light sources with different wavelengths. While it opens up a capability of a single-shot measurement of surface normal, the problem has been known ill-posed. To make the problem well-posed, existing MPS methods rely on restrictive assumptions, such as shape prior, surfaces having a monochromatic with uniform albedo. This paper alleviates these restrictive assumptions in existing methods. We show that the problem becomes well-posed for surfaces with uniform chromaticity but spatially-varying albedos based on our new formulation. Specifically, if at least three (or two) scene points share the same chromaticity, the proposed method uniquely recovers their surface normals with the illumination of no less than four (or five) spectral lights in a closed-form. In addition, we show that a more general setting of spatially-varying both chromaticities and albedos can become well-posed if the light spectra and camera spectral sensitivity are calibrated. For this general setting, we derive a unique and closed-form solution for MPS using the linear bases extracted from a spectral reflectance database. Experiments on both synthetic and real captured data with spatially-varying reflectance demonstrate the effectiveness of our method and show the potential applicability for multispectral heritage preservation.

Appendix: Limitation of 3-Band Multispectral Photometric Stereo Methods

As described in Sec. 2, existing methods Ozawa et al. (2018); Chakrabarti and Sunkavalli (2016) provide a unique solution for the monochromatic surface with a uniform albedo (SRT II) using 3-channel RGB images. These methods are limited to 3 channels, and it is not straightforward to extend them to take more channels as input. However, the capability of taking more channels is favorable because it allows us to use robust estimation techniques, effectively neglecting outliers such as shadows and specular reflections. Here we show the reason why the existing methods are limited to 3-channel input.

Without loss of generality, we fix the common albedo ${\tilde{\rho }}$ for all the scene points on SRT II surface to be 1 and define a diagonal matrix $\mathbf {Q}$ as $\mathbf {Q} = \mathrm{diag}(\mathbf {q})$. Following Eq. (4), the image observation for a pixel can be represented as

$$\begin{aligned} \mathbf {m} = \mathbf {QL}\mathbf {n}. \end{aligned}$$

(22)

Defining Moore-Penrose inverse matrix $\mathbf {K} \in {\mathbb {R}}^{3 \times f}$ as $\mathbf {K} = (\mathbf {QL})^{\dagger }$, the surface normal is then calculated by

$$\begin{aligned} \mathbf {n} = \mathbf {K}\mathbf {m}. \end{aligned}$$

(23)

The existing methods Ozawa et al. (2018); Chakrabarti and Sunkavalli (2016) use a unit norm constraint about a surface normal as

$$\begin{aligned} \mathbf {n}^\top \mathbf {n} = \mathbf {m}^\top \mathbf {K}^{\top } \mathbf {K}\mathbf {m} = 1. \end{aligned}$$

(24)

As shown in Eq. (25), by defining $\mathbf {E} = \mathbf {K}^{\top } \mathbf {K} \in {\mathbb {R}}^{f \times f}$, each one of the p scene points provides an equation about $\mathbf {E}$ as

$$\begin{aligned} \left\{ \begin{array}{lr} \mathbf {m}_0^\top \mathbf {E} \mathbf {m}_0 = 1, \\ \mathbf {m}_1^\top \mathbf {E} \mathbf {m}_1 = 1, \\ \quad \quad \vdots \\ \mathbf {m}_p^\top \mathbf {E} \mathbf {m}_p = 1. \\ \end{array} \right. \end{aligned}$$

(25)

Defining $\mathbf {m} \otimes \mathbf {m} = \mathrm{vec}(\mathbf {m}\mathbf {m}^\top )$, we rewrite Eq. (25) in a matrix form,

$$\begin{aligned} \underbrace{ \left[ \begin{array}{l} \mathbf {m}_0 \otimes \mathbf {m}_0 \\ \mathbf {m}_1 \otimes \mathbf {m}_1 \\ \vdots \\ \mathbf {m}_p \otimes \mathbf {m}_p \\ \end{array} \right] }_{\mathbf {G}} \underbrace{ \left[ \begin{array}{l} \mathrm{vec}(\mathbf {E}) \end{array} \right] }_{\mathbf {y}} = \mathbf {1}, \end{aligned}$$

(26)

where $\otimes$ represents the Kronecker product, and $\mathbf {H}$ forms a $p \times f^2$ matrix. Since $\mathbf {E} \in {\mathbb {R}}^{f \times f}$ is symmetric, $\mathbf {y}$ only has at most $\frac{f(f+1)}{2}$ distinct elements. We extract the elements of $\mathbf {y}$ that correspond to the upper triangle elements from $\mathbf {E}$ as $\mathbf {z} \in {\mathbb {R}}^{\frac{f(f+1)}{2}}$ and the corresponding columns from $\mathbf {H}$ as $\hat{\mathbf {H}} \in {\mathbb {R}}^{p \times \frac{f(f+1)}{2}}$. Then we rewrite Eq. (26) as

$$\begin{aligned} \hat{\mathbf {H}}\mathbf {z} = \mathbf {1}. \end{aligned}$$

(27)

The necessary condition to obtain a unique approximate solution for $\mathbf {z}$ is $\hat{\mathbf {H}}$ to have full-rank, i.e., assuming $p \ge f(f+1)$,

$$\begin{aligned} \mathrm{rank}(\hat{\mathbf {H}}) = \frac{f(f+1)}{2}. \end{aligned}$$

(28)

On the other hand, since the image observations for all the scene points under Lambertian reflectance has the rank of 3, we can represent any irradiance measurements with three independent basis $\{\mathbf {e}_1, \mathbf {e}_2, \mathbf {e}_3 \in {\mathbb {R}}^f\}$, i.e.,

$$\begin{aligned} \mathbf {m} = c_1\mathbf {e}_1 + c_2\mathbf {e}_2 + c_3\mathbf {e}_3. \end{aligned}$$

(29)

With this expression, we can represent $\mathbf {m} \otimes \mathbf {m}$ as

$$\begin{aligned} \begin{aligned} \mathbf {m} \otimes \mathbf {m}&= (c_1\mathbf {e}_1 + c_2\mathbf {e}_2 + c_3\mathbf {e}_3) \otimes (c_1\mathbf {e}_1 + c_2\mathbf {e}_2 + c_3\mathbf {e}_3)\\&= c_1^2(\mathbf {e}_1 \otimes \mathbf {e}_1) + 2c_1c_2(\mathbf {e}_1 \otimes \mathbf {e}_2) + c_2^2(\mathbf {e}_2 \otimes \mathbf {e}_2) \\&\quad + 2c_1c_3(\mathbf {e}_1 \otimes \mathbf {e}_3) + 2c_2c_3(\mathbf {e}_2 \otimes \mathbf {e}_3) + c_3^2(\mathbf {e}_3 \otimes \mathbf {e}_3). \end{aligned} \end{aligned}$$

(30)

It indicates that $\mathbf {m} \otimes \mathbf {m}$ can be represented by at most 6 independent f-dimensional basis vectors $\mathbf {e}_i \otimes \mathbf {e}_j$. Since $\mathbf {H}$ in Eq. (26) is a stack of $\mathbf {m} \otimes \mathbf {m}$, the rank of $\mathbf {H}$ should satisfy

$$\begin{aligned} \mathrm{rank}(\mathbf {H}) \le 6. \end{aligned}$$

(31)

Together with the necessary condition in Eq. (28) for solving Eq. (27), it leads to the following inequality,

$$\begin{aligned} \frac{f(f+1)}{2} = \mathrm{rank}(\hat{\mathbf {H}}) \le \mathrm{rank}(\mathbf {H}) \le 6, \end{aligned}$$

(32)

which indicates that the number of spectral channels f of the input multispectral image should be no more than 3. Therefore, these existing method Ozawa et al. (2018); Chakrabarti and Sunkavalli (2016) cannot be adapted to multispectral images with more than three bands. On the other hand, our method is free from this restriction.