Bimodal Stereo: Joint Shape and Pose Estimation from Color-Depth Image Pair

TL;DR

This paper introduces a novel iterative framework for jointly estimating shape and pose from uncalibrated color and depth image pairs, improving accuracy in multi-modal data registration.

Contribution

It proposes a Shape-from-Shading based method that refines shape and pose simultaneously, addressing noise and holes in depth maps.

Findings

Effective shape and pose recovery demonstrated

Improved registration accuracy in experiments

Robustness to noisy depth data

Abstract

Mutual calibration between color and depth cameras is a challenging topic in multi-modal data registration. In this paper, we are confronted with a "Bimodal Stereo" problem, which aims to solve camera pose from a pair of an uncalibrated color image and a depth map from different views automatically. To address this problem, an iterative Shape-from-Shading (SfS) based framework is proposed to estimate shape and pose simultaneously. In the pipeline, the estimated shape is refined by the shape prior from the given depth map under the estimated pose. Meanwhile, the estimated pose is improved by the registration of estimated shape and shape from given depth map. We also introduce a shading based refinement in the pipeline to address noisy depth map with holes. Extensive experiments showed that through our method, both the depth map, the recovered shape as well as its pose can be desirably refined and recovered.

Tables

Table 1. Table 1: Estimated pose paramters of Figure  4

	s	$\alpha$	$\beta$	$\gamma$
Esimated	1.0000	0.03216∘	18.7013∘	0.01812∘
Groundtruth	1.0000	0∘	20∘	0∘

Table 2. Table 2: Error of estimated rotation under different rotate angles and different wideness of overlap area.

$R^{\ast }$ / $P_w$	0.125	0.25	0.375	0.5	0.75	1
20	0.022806	0.015999	0.006442	0.006005	0.007384	0.008664
40	0.022903	0.011479	0.006860	0.005768	0.005786	0.006653
60	0.022379	0.018440	0.012250	0.007660	0.004729	0.004509
80	0.010623	0.008897	0.008637	0.005853	0.002500	0.001922
90	0.007870	1.631535	0.006505	0.004082	0.000622	0.000649

Table 3. Table 3: Estimated pose paramters of Figure  6

	s	$\alpha$	$\beta$	$\gamma$
Esimated Pose	0.9972	-0.4813∘	17.3492∘	9.3220∘

Equations

$$\underset{{\mathbf{n}}^{(k)},{\mathbf{z}}^{(k)},{\mathcal{T}}^{(k)}(\cdot )}{\min }\ell ({\mathbf{n}}^{(k)},{\mathbf{z}}^{(k)},{\mathcal{T}}^{(k)}(\cdot ))+f({\mathbf{z}}^{(k)},{\mathcal{T}}^{(k)}(\cdot ))+g({\mathbf{n}}^{(k)},{\mathcal{T}}^{(k)}(\cdot )),$$

$$h(\mathbf{L})={\lambda }_L(\mathbf{L}-{\mu }_L{)}^T{\mathbf{\Sigma }}_L^{-1}(\mathbf{L}-{\mu }_L),$$

$$\log S_i={{\hat{\mathbf{n}}}_i}^T{\mathbf{M}}_{\mathbf{j}}{\hat{\mathbf{n}}}_i$$

$${\left[\begin{array}{cccc}-\frac{p_i}{\sqrt{1+p_i^2+q_i^2}},& -\frac{q_i}{\sqrt{1+p_i^2+q_i^2}},& \frac{1}{\sqrt{1+p_i^2+q_i^2}},& 1\end{array}\right]}^T.$$

$${\mathbf{M}}_j=\left[\begin{array}{cccc}{c}_1{L}_{j,9}& c_1{L}_{j,5}& c_1{L}_{j,8}& c_2{L}_{j,4}\\ c_1{L}_{j,5}& -c_1{L}_{j,9}& c_1{L}_{j,6}& c_2{L}_{j,2}\\ c_1{L}_{j,8}& c_1{L}_{j,6}& c_3{L}_{j,7}& c_2{L}_{j,3}\\ c_2{L}_{j,4}& c_2{L}_{j,2}& c_2{L}_{j,3}& c_4{L}_{j,1}-c_5{L}_{j,7}\end{array}\right]$$

$$c_1=0.429043,c_2=0.511664,c_3=0.743125,$$

$$c_4=0.886227,c_5=0.247708$$

$$\underset{{\mathbf{n}}_i^{\ast }}{\min }\ell ({\mathbf{n}}_i^{\ast },z_i^{\ast },\mathcal{T}(\cdot ))$$

$$s.t.{{\mathbf{n}}_i^{\ast }}^T{\mathbf{n}}_i^{\ast }=1$$

$$\underset{{\mathbf{n}}_i^{\ast }}{\min }\ell ({\mathbf{n}}_i^{\ast },z_i^{\ast },\mathcal{T}(\cdot ))$$

$$+{\lambda }_3({{\mathbf{n}}_i^{\ast }}^T{\mathbf{n}}_i^{\ast }-1)$$

$${\ell }_{\text{1}}({\mathbf{n}}_i^{\ast })={(\log S_{ij}-{{\hat{\mathbf{n}}}_i}^T{\mathbf{M}}_j{\hat{\mathbf{n}}}_i)}^2,$$

$${\ell }_2({\mathbf{n}}_i^{\ast },{\mathbf{n}}_{R_h},\mathcal{T}(\cdot ))=\Vert {\mathbf{n}}_i^{\ast }-\mathcal{T}({\mathbf{n}}_{R_h})\Vert ,$$

$${\mathbf{n}}_i^{\ast }={\left[\begin{array}{ccc}-\frac{p_i}{\sqrt{1+p_i^2+q_i^2}},& -\frac{q_i}{\sqrt{1+p_i^2+q_i^2}},& \frac{1}{\sqrt{1+p_i^2+q_i^2}}\end{array}\right]}^T,$$

$$\{\begin{array}{c}\frac{\partial z_i^{\ast }}{\partial x}=p_i\\ \frac{\partial z_i^{\ast }}{\partial y}=q_i\end{array}$$

$$f(z^{\ast },s,\mathbf{R},\mathbf{t})=\underset{s,\mathbf{R},\mathbf{t}}{\mathrm{minimize}}{\Vert P_c^{\ast }-s\mathbf{R}{P}_c-\mathbf{t}\Vert }^2$$

$$T_{\theta }(\mathbf{y})=s\mathbf{Ry}+\mathbf{t}=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]\mathbf{y}+\left[\begin{array}{c}{t}_1\\ t_2\\ t_3\end{array}\right],$$

$$\theta =$$

$$\begin{array}{ccc}{a}_{23}& a_{31}& \begin{array}{ccc}\begin{array}{ccc}{a}_{32}& a_{33}& t_1\end{array}& t_2& t_3\end{array}\end{array}{]}^T$$

$$\mathbf{b}=\left[\begin{array}{c}{x}^{(1)}\\ ⋮\\ ⋮\\ x^{(n)}\end{array}\right]=\left[\begin{array}{c}{T}_{\theta }(y^{(1)})\\ ⋮\\ ⋮\\ T_{\theta }(y^{(n)})\end{array}\right]=\mathbf{A}\theta$$

$$\mathbf{R}=s{\mathbf{R}}_x(\alpha ){\mathbf{R}}_y(\beta ){\mathbf{R}}_z(\gamma )$$

$$g(\alpha ,\beta ,\gamma )=\underset{\alpha ,\beta ,\gamma }{\min }{\Vert {\mathbf{R}}_x(\alpha ){\mathbf{R}}_y(\beta ){\mathbf{R}}_z(\gamma )-\mathbf{R}\Vert }^2$$

$$\min g({\mathbf{n}}_{R_h}^{\ast })$$

$$+{\lambda }_{g3}({{\mathbf{n}}_{R_h}^{\ast }}^T{\mathbf{n}}_{R_h}^{\ast }-1).$$

$$g_{\text{1}}({\mathbf{n}}_{R_h}^{\ast })=(\log S_{i}j-({\hat{\mathbf{n}}}_{R_h}^{\ast }{)}^T{\mathbf{M}}_j({\hat{\mathbf{n}}}_{R_h}^{\ast }){)}^2,$$

$$g_2({\mathbf{n}}_{R_h}^{\ast })=\Vert {\mathbf{n}}_{R_h}^{\ast }-{\mathcal{T}}^{-1}({\mathbf{n}}_{\mathbf{i}}^{\ast })\Vert ,$$

$$err=\Vert \frac{{\text{R}}^{\ast }}{\Vert {\text{ R}}^{\ast }\Vert }-\frac{{\text{R}}_g}{\Vert {\text{ R}}_g\Vert }\Vert$$

Full text

††journal: Machine Vision and Applications

1 Introduction

With the development of depth sensors (time-of-flight (ToF) camera, laser range scanner, structured light scanner, etc.), researchers have been dedicating themselves on utilizing both color and depth information for 3D scene reconstruction and related topics. The calibration of color and depth cameraHanTOC13 , especially extrinsic calibration (i.e. pose estimation between cameras), thus becomes a prerequisite in building such applications. Existing manual and semi-automatic methodsHerreraPAMI12 ; ZhangICME11 are solving this problem in the settings where multiple color images and checkerboard are required.

In this paper, we are confronted with a similar but different problem, which aims to solve camera pose from a pair of an uncalibrated color image and a depth map from different views automatically (Illustrated in Figure 1). Such problem has varieties of difficulties that existing method may be failed. Firstly, 3D coordinates of feature points in image coordinate system could not be easily obtained from an uncalibrated color image. Moreover, point-wise correspondence may be difficult to generate between this pair of multi-modal images.

To solve such “Bimodal Stereo” problem, we take the assumption of Spherical Harmonics (SH) illumination RamamoorthiSIGGRAPH01 based image formation process and propose an iterative Shape-from-Shading (SfS) based framework to estimate shape and pose simultaneously, shown in Figure 2. In each iteration, the input data is processed through a pipeline of three procedures, namely Shape from shading, Point clouds registration and Shading based shape refinement. This method is based on the notation that given sufficient shape prior the SfS reconstruction could be largely improved (see details in Section 7.1.1). In further, more accurate pose can be estimated through registration between reasonable reconstruction and given depth map , which in return benefits the SfS process. Therefore, given the estimated lighting parameters, an uncalibrated color image and depth map, the proposed EM-like algorithm is able to obtain plausible shape and pose.

The contribution of our work is listed as follows:

Firstly, an interesting problem of Bimodal Stereo is introduced. This problem is non-trivial since it aims to automatically solve the pose between a pair of a single uncalibrated image and a depth map.

Besides, Spherical Harmonics Shape from shading method is utilized to solve the Kinect extrinsic calibration fully automatically and targetlessly. To our knowledge, it is the first attempt to register multi-modal images by image formation model, other than feature-based methods HeinrichMIA11 or mutual information (MI) based methodsPluimMICCAI00 . With Spherical Harmonics modelling the illumination, our method becomes a well-posed one, which can be effectively solved by non-linear least square methods.

Moreover, the cost function of SH Shape from shading is regularized by imposing a simple shape prior term, which bridges the shading information from color image and the shape information from depth map by the estimated pose between two cameras.

This paper will proceed as follows: In Section 2, related works are reviewed and compared with the proposed method. In Section 3 our method will be formulated as a minimization problem over Shape from Shading, point clouds registration and shading based shape refinement. In Sections 4, 5, and 6 the cost function and solution of each process are described respectively. In Section 7 we will not only show the results on synthesized data and real data, but also conduct some quantitative experiments to validate the assumption and the robustness of our method. In the last section, this paper will be concluded while the limitations of proposed method are discussed.

2 Related Works

Although the physical process of image formation is still under discussion yet, the majority of researchers in computer vision area would agree that images convey shape information implicitly. In 1975 Horn proposed the well-known ill-posed inverse problem called ”Shape from Shading (SfS)” ,which aims to recover shape from a single grayscale image. The problem is based upon the assumption that the shading information of an image is formed by a known light source and surface normals when the albedo of the surface is constant.

From then on, many efforts have been made to solve this problem under simplified assumptions on illumination. The most classical setting is to recover the shape of a Lambertian surface, i.e. a diffusely reflecting surface, under the lighting of a distant point source. From early 80’s to now, various method proposed, including directly numerical solutionsIkeuchiAI81 ; PradosICCV03 ,linear approximationTsaiIVC94 and minimizationCrouzilTPAMI03 ,achieved success in solving the problem under this oversimplified assumption. Since Basri et al.BasriTPAMI03 have proved images of a convex Lambertian object under a wide variety of lighting conditions can be approximated accurately by a 9D SH coefficients vector, researcher starts to pay attention on SH illumination modelRamamoorthiSIGGRAPH01 which achieved great success in high quality rendering. Instead of solving the surface gradients or depth in existing method, our method seek for the unit surface normal which satisfies the RGB-shading model, thus the problem become well-posed if the three coefficients vectors for R, G and B channel respectively are whitenedBarronPAMI15 . Moreover, we posed simple shape prior term in our SfS method as a bridge between the two multi-modal images instead of other various regularization termsZhangTPAMI99 , which makes our cost function parsimony while still strong enough.

Currently, while some of the researchers tried to add some constraints making SfS problem well-posedPradosCVPR05 , others YuCVPR13 ; WuTOG14 turned to embed SfS techniques to refine the noisy and holey depth map obtained by Kinect. In Yu’s methodYuCVPR13 the SH coefficients vectors for each channel are not whitened hence cannot ensure the property well-posed and generally slower with more terms than ours. Wu’s method WuTOG14 is the first real-time shading based shape refinement method. However, since the method directly refines the depth information, it is hard to impose our regularization of shape prior. Our method, instead, can be viewed as a symmetric version of proposed SfS method. Therefore it is also a well-posed problem which solves surface normal at each pixel.

As stated before, the proposed Bimodal stereo problem can be viewed as an automatically extrinsic calibration between color camera and depth camera of a Kinect. Though factory calibrated, the calibration of Kinect is still important to applications in high-accuracy 3d reconstruction and adding new cameras. Generally, existing methods HerreraPAMI12 ; ZhangICME11 are derived from Zhang’s RGB camera calibration methodZZhangTPAMI00 . Both methods requires calibration card to calibrate the color image first before extrinsic calibration by feature matching on both images. Our method, however, does not require RGB camera calibration. Moreover, the pose between cameras is solved automatically, which meets the requirement of Bimodal Stereo problem.

3 Problem Formulation

3.1 Method Overview

In this paper an iterative computational framework is proposed to solve problem Eq. (1), as shown in Algorithm 1. and illustrated in Figure 2. The input of the system is depth map $z_{R}$ and RGB log-shading $\{S_{ij}\}_{i=1,...,M,j=1,2,3}$, where $M$ is the total amount of pixels and $j=1,2,3$ denotes the R, G or B log-values of the i-th pixel. The goal is to estimate normal $\mathbf{n}^{(k)}$ and depth $\mathbf{z}^{(k)}$ from the shading information and the refined depth $z_{R}^{(k)}$ as well as the motion $\mathcal{T}^{(k)}$ between RGB camera and depth camera. In the loop, we first estimate $\mathbf{n}$ in a shape-from-shading way based on Spherical Harmonics Illumination model. Then the depth $\mathbf{z}^{(k)}$ can be obtained by solving the partial differential equations of estimated surface gradients $p,q$ from $\mathbf{n}^{(k)}$. Afterwards, the transformation $\mathcal{T}^{(k)}(\mathbf{P})=s\mathbf{RP}+\mathbf{t}$ is estimated by registering the point clouds $P_{c}^{(k)}=\{(x,y,z^{(k)})\}$ from $z^{(k)}$and $P_{c}$ from $z_{R}^{(k)}$. Also, a mapping describes the point wise correspondence between the two point clouds is established during the registration. The last procedure of one iteration is to refine the overlapping area of given depth $z_{R}$ and its normal $n_{R}$ based on the rotation $\mathbf{R}$ and shading $\{S\}$ in this area. Then the improved normal $\mathbf{n}_{R}^{(k)}$ is utilized as the initial input in the next iteration. The loop will end when the norm of the difference two matrix $\mathbf{R}^{(k)}$ and $\mathbf{R}^{(k-1)}$ is less than a threshold set by user.

3.2 Formulation

We call our problem formulation for iteratively estimating shape information from the unknown correspondence of a RGB image and a Kinect depth map from different views as ”Bimodal Stereo”. Bimodal Stereo can be thought of as a shape completion and refinement by multi modal data registration. It is an integrated iterative process including three stages Shape from Shading(SfS), Point Cloud Registration and Shape refining. The ”Bimodal Stereo” problem formulation is :

[TABLE]

where $\mathbf{n}^{(k)}$ and $\mathbf{z}^{(k)}$ is the estimated normal and depth map of the surface from the RGB shading information. $\mathcal{T}(\mathbf{P})=s\mathbf{RP}+\mathbf{t}$ is the estimated 3D transformation for a point cloud $\mathbf{P}$. In transformation $\mathcal{T}$, $s$ is the scale factor, $\mathbf{R}$ is the $3\times 3$ rotation matrix between point clouds $P_{c}^{(k)}=\{(x,y,z^{(k)})\}$ and $P_{c}$ from depth map and $\mathbf{t}$ is the $3\times 1$ translation vector. These parameters are estimated by minimizing the sum of cost functions $\ell(\mathbf{n}^{(k)},\mathbf{z}^{(k)},\mathcal{T}(\cdot))$,$f(\mathbf{z}^{(k)},\mathcal{T}^{(k)}(\cdot))$ and $g(\mathbf{n}^{(k)},\mathcal{T}^{(k)}(\cdot))$ for SfS, registration and shape refining respectively.

Therefore the detailed description of our method is divided into 3 parts. SH based Shape from Shading is explained in the next chapter while registration and shape refinement are discussed in Chapter 4 and 5 respectively. In Chapter 6 we’ll show the validation and quantitative analysis of our method.

4 Spherical Harmonics Shape from Shading

In this section, we present our Shape from Shading method. Based on the assumption that the RGB image is formed by Spherical Harmonics Lighting model and surface normals, the SfS process can be formulated as an inverse problem of this image formation model. By posing the shape prior in the SfS formulation, the correlation between estimated shape and estimated pose is set up. Such problem is well-posed and can be solved by nonlinear least square method. Once obtained the solution normal, the 3D point cloud of estimated shape can be computed through the least squares surface reconstruction.

4.1 Spherical Harmonics Lighting Model

Traditional shape from shading methods assume the illumination model as a sufficiently remote point source. Such simple model though effective in shape reconstruction from synthetic images, could not perform well with images acquired in the natural world. The Spherical Harmonics (SH) Illumination however provide a more accurate while not complex description of the natural lighting. Instead of modelling the illuminant in whole, the SH model assumes that at an arbitrary point in the 3D space, the light direction could be anisotropy. Such complex lighting condition could be parametrized by a 27 (color, 9 dimensions per RGB channel) dimensional vector, $L$.

After fitting a multivariate Gaussian to the SH illuminations in the training set, the estimated SH coefficients vector $L$ is obtained by minimizing the following cost functionBarronPAMI15

[TABLE]

where $\mu_{L}$ and $\Sigma_{L}$ are learned parameters of the Gaussian, while $\lambda_{L}$ is the multiplier learned on dataset. To ensure the coefficients vector of each channel to be independent with each other, we preprocess these vectors by whitening transformation described in BarronPAMI15 .

The log-intensity at pixel $i$, $\log S_{i}$ of a Lambertian surface under SH illumination could be described in the following equation:

[TABLE]

where $\mathbf{\hat{n}}_{i}$ is a 4-dimensional vector:

[TABLE]

In this Equation, $p_{i}=\frac{\partial z_{i}}{\partial x}$ and $q_{i}=\frac{\partial z_{i}}{\partial y}$ is the partial derivative of the depth map $z$ at pixel $i$ along $x$ and $y$ axis, respectively. It is worth to point that $\mathbf{\hat{n}_{i}}=\left[\begin{matrix}\mathbf{n}_{i}^{T},&1\end{matrix}\right]^{T}$ where $\mathbf{n}_{i}^{T}$ is the surface normal at pixel $i$.

$\mathbf{M}_{j}$ is a $4\times 4$ matrix for the j-th RGB channel, described by

[TABLE]

where

[TABLE]

and $\mathbf{L}_{j}$ denote the SH coefficients of j-th RGB channel.

4.2 Shape from Shading with shape prior

Recovering surface normal $\mathbf{n}_{i}^{*}$ from the corresponding brightness $S_{ij},j=1,2,3$ could be viewed as a constrained optimization problem for $\mathbf{n}_{i}^{*}$ on the following equation:

[TABLE]

where $j=1,2,3$ denotes the R,G,B channel and the corresponding point of the $i$-th pixel is $h$-th point on the given depth map. $\mathcal{T}(\cdot)$ is the estimated pose. The constraint is the norm regularization which enforces normals to be of unit length. We can pose the constraint ${\mathbf{n}_{i}^{*}}^{T}\mathbf{n}_{i}^{*}=1$ into $\ell(\mathbf{n}_{i}^{*},z_{i}^{*},\mathcal{T}(\cdot))$ so that the solution of Eq. (6) can be converted to the following unconstrained optimization problem:

[TABLE]

In Eq. (6), ${{\ell}_{1}}(\mathbf{n}_{i}^{*})$ is the Brightness Constraint measuring the total squared brightness error over the grayscale image $S$:

[TABLE]

The second term ${{\ell}_{2}}(\mathbf{n}_{i}^{*},\mathbf{n}_{R_{h}},\mathcal{T}(\cdot)))$ serves as Shape Prior, which enforces the estimated normals $\mathbf{n}_{i}^{*}$ of $z_{i}^{*}$ to be close to the rotated ground truth normals $\mathbf{n}_{R_{h}}$ from depth map $z_{R_{h}}$. This term is modelled as

[TABLE]

where $\mathbf{R}$ is the estimated rotation matrix between RGB image and depth image.

Our chose to solve normal $\mathbf{n^{*}}$ instead of $p$ and $q$ is of advantages for two reasons. Firstly, $p$ and $q$ is unbounded, i.e. , the absolute value of $p$ or $q$ can be much larger than 1 while the elements in the normalized vector $\mathbf{n^{*}}$ vary among $\left[-1,1\right]$. Secondly, in traditional SfS method, the relationship between $p$ and $q$ is not fully utilized. However, in the proposed method, $p$ and $q$ are solved together in the sense of $\mathbf{n}_{i}^{*}$, namely,

[TABLE]

which is more suitable for the SH based shape from shading.

4.3 Optimization

Since for $i$-th pixel a system of 3 equations can be set up for the 3 unknown variables $\mathbf{n_{i}^{*}}=\begin{bmatrix}n_{i1}^{*}&n_{i2}^{*}&n_{i3}^{*}\end{bmatrix}^{T}$, trust-region-reflective algorithm is chosen to solve this well-posed nonlinear least square system. Compared with gradient descent methods, it does not need to compute actual partial derivatives of the function and runs faster generally. The optimization toolbox is integrated and provided in Matlab. The maximum number of iteration is set to be 1000, while in most cases the method will converge in less than 300 times of iteration.

4.4 Postprocessing

Once the estimated unit normal $\mathbf{n_{i}^{*}}$ is obtained, we utilized the method on least squares surface reconstruction from gradients proposed by Harker et al HarkerCVPR08 ; HarkerCVPR11 to solve the following system for $z_{i}^{*}$:

[TABLE]

After obtaining the depth map $z_{i}^{*}(x,y)$, a point cloud $\mathbf{P}_{c}^{*}=\{(x,y,z^{*}))\}$ can be formed for further registration.

5 Point Cloud Registration

In this section, we will discuss how to estimate the pose $\mathcal{T(\cdot)}$ from $P_{c}$ to $P_{c}^{*}$, where $P_{c}$ is the point cloud from Kinect depth map $z_{R}$. This problem to estimate the pose can be viewed as estimating the scale factor $\mathbf{s}$, rotation matrix $\mathbf{R}$ and translation vector $\mathbf{t}$ by minimizing the correspondent point wise distance between estimated point clouds $P_{c}^{*}$ and $P_{c}$:

[TABLE]

In one iteration, the estimated point cloud is fixed hence the problem is simplified as 3D point clouds registration. Iterative Closest Point (ICP) method is chosen to obtain an initial transform and set up the primary point-wise correspondences while RANSAC is chosen to reject the weak correspondences and get the refined pose $\mathcal{T(\cdot)}$.

5.1 Pose Estimation

The following form of RST transformation model is utilized in the RANSAC process:

[TABLE]

so that parameter vector is

[TABLE]

Given $n$ point correspondences ${{\left\{{{x}_{i}},{{y}_{i}}\right\}}_{i=1,..,n}}$ these transformation relationship can be grouped as following linear system

[TABLE]

Then the solution of Eq. (14), $\mathbf{\hat{\theta}}$ can be estimated through linear least square method.

Once got the optimal $\mathbf{\theta}$ and $\mathbf{R}$ when $CS$ between two point sets is the largest, the global rotation matrix $\mathbf{R}$ can be decomposed into the rotations along the three axis respectively:

[TABLE]

where $\alpha$,$\beta$ and $\gamma$ are rotated angles, namely yaw, pitch, and roll along x-axis, y-axis and z-axis respectively, $\mathbf{R}(\cdot)$ is the corresponding rotation matrix.

Such problem can be easily formulated as a nonlinear least square estimation problem, which aims to obtain $\hat{\alpha}$,$\hat{\beta}$ and $\hat{\gamma}$ by minimizing

[TABLE]

Levenberg-Marquardt method is leveraged to optimize the Eq. (16). Compared to the trust-region-reflective algorithms, LM method is generally faster and hence will converge in several iterations in most cases.

5.2 Point-wise Correspondences Set Up

Theoretically, the estimated pose $\mathcal{T(\cdot)}$ can determine a pair of quadrilateral overlapping areas on depth maps $\mathbf{Z}^{*}$ and $\mathbf{Z}_{R}$ respectively. However, in practice depth maps are represented by point arrays. Therefore we set up the one-to-one mapping between points of overlap areas on the different depth maps in following three steps:

• 1

Map $P_{c}$ onto the space of $P_{c}^{*}$ by apply the transformation with parameters $\mathcal{T}(P_{c})$ of estimated pose.

• 2

For each point in $\{\mathcal{T}(P_{c})\}$, find the 1st nearest neighbour in the point set $\{P_{c}^{*}\}$ and hence to form an initial correspondence.

• 3

Accept the point wise correspondence based on the distance $d(P_{c}^{*},\mathcal{T}(P_{c}))$ between $P_{c}^{*}$ and $\mathcal{T}(P_{c})$.

A naive modelling of acceptance criterion can be hard-counting, namely accepting the points pair if the correspondent distance $d(P_{c}^{*},\mathcal{T}(P_{c}))$ is less than a threshold. In our method the threshold is set to be 1.

6 Shading-based Depth Refining

Given depth map may contain holes or noises. Since normal $\mathbf{n}_{R_{h}}$ and depth $z_{R_{h}}$ of $h$-th pixel in the overlap area of given depth map are subjected to the shading model Eq. (3) at $i$-th pixel on the RGB image, if the pose $\mathcal{T}(\cdot)$ is determined, we can refine the normal $\mathbf{n}_{R_{h}}$ in the a symmetric way of proposed SfS process, i.e. ,

[TABLE]

The brightness constraint at $k$-th normal on the given depth map can be formed as

[TABLE]

where $\mathbf{\hat{n}}_{R_{h}}=\left[\begin{matrix}\mathbf{R}^{-1}\mathbf{n}_{R_{h}}^{*},&1\end{matrix}\right]^{T}$.

The shape prior term in the normal refinement is nothing but the symmetric version of Eq. (9):

[TABLE]

where $\mathcal{T}^{-1}(\cdot)$ is the inverse transformation of the estimated pose $\mathcal{T}(\cdot)$, which maps $P_{c}^{*}$ onto $P_{c}$. The refined depth at $h$-th pixel on the depth map $z_{R_{h}}^{*}$ can also be solved in the same manner of shape from surface gradients method.

7 Experiments

In this section, we will validate the proposed method on both synthesized data and real RGB-D data. Both qualitative and quantitative analysis are conducted not only to show the results of our method, but also to explore the reason why our method may work and how the method perform under different pose and different wideness of overlap area.

7.1 Synthetic Study

Two quantitative experiments are designed and conducted to explore how and to what extent the shape prior information can be utilized in the shape from shading process as well as to validate the performance of our methods on the $32\times 32$ synthesized pair of RGB image and rotated depth map of Mozart’s face. The synthesis data can be obtained in the following steps, shown in Figure 3:

a. Given certain SH coefficients and a frontal depth map of Mozart’s face, the log-RGB image can be obtained through Eq. (3).

b. Given certain rotation angle $\alpha,\beta,\gamma$ along x-,y- and z-axis respectively, the rotated depth map can be generated.

c. The wideness of the baseline is determined by the grid we use to sample the rotated depth map.

Figure 4 shows the experiment conducting on the synthesis RGB image under given lighting conditions and depth map of Mozart’s face from another view. The depth map is obtained by rotating the ground truth depth for 20 degrees along y-axis. The loop terminated after 16 iterations. During the process the estimated depth $z^{*}$ is obviously improved, especially the right upper part of the shape. The estimated pose paramters are listed in Table 1.

Besides, we conducted two quantitative experiments to further study the capability of the proposed method. One is to explore the reason why shape prior can improve the reconstruction result of the proposed SfS method, the other is to discuss the performance under different rotations and wideness of baseline.

7.1.1 Experiments on Shape Prior

In the first experiment, we will validate the underlying assumption of our method, i.e. by posing the point-wise shape prior term in the original brightness equation Eq. (8), the SfS reconstruction can be improved.

The experiment is conducted by controlling the percentage $P_{er}$ of the shape prior provided. $P_{er}=0$ when no prior is offered while $P_{er}=1$ means every point is given the prior information.

The result of this experiment can be shown in Figure 5. According to the results, recovered shape is improved and approaching the ground truth when $P_{er}$ is increased. In this sense we can conclude that the SfS reconstruction result of certain image area can be improved by given shape prior of this area.

7.1.2 Experiments on Different Rotation and Baseline

The next experiment is to discover the method performance in ”Bimodal” settings, i.e. the depth map and RGB images are synthesized from different views. We’d like to simulate the situation where the depth camera and RGB camera shares partial overlap views and study the relationship between estimation accuracy and parameters, i.e. rotation and baseline, of given depth map. The term ’baseline’ from the Stereopsis to describe the distance of two optical centers, which describes the wideness of overlap area.

In the experiment, the rotation angle along y-axis $\beta$ is chosen from $\{20,40,60,80,90\}$ and the wideness of overlap area $P_{w}$, i.e. the percentage of the total area, is chosen from $\{0.125,0.25,0.375,0.5,0.75,1\}$. The normalized difference of estimated rotation $R^{*}$ and ground truth $R$ is computed in the following equation:

[TABLE]

The numerical results are shown in Table 2 while in Figure 7 the illustration of wideness of the overlap area is shown in (a). In (b) and (c) we show the error variation at 40*∘* when wideness changes as well as the visualization (right side in (c)) and registration (left side in (c)) of estimated shape and given point cloud. The error is decreasing when the overlap area becomes larger. However in (d) of Figure 7, the failure case appears when $\beta$ reaches 90 and $P_{w}$ is 0.25, which is visualized in (e).

As shown, since the intrinsics of the depth camera are taken as known parameters and hence can obtain the point cloud easily, the method is robust in most rotations and wideness of overlap area.

7.2 Real Data Study

The proposed Bimodal stereo method is validated on RGB-D Object Dataset proposed by Lai* et al* LaiICRA11 . This dataset contains over 300 instances in 51 categories, where the dataset of each instance is composed of cropped RGBD images of the object from different views. The experiments are conducted on the data of foodbox, mug and lemon, etc.

The reconstruction and refinement on the instance foodbox is displayed in Figure 6 while the estimated pose parameters are shown in Table 3.

In this case, we chose foodbox-1-1-5 as RGB input and foodbox-1-1-18 as depth input, the pose between which could be observed obviously. Besides, the depth map contains small holes on the surface of foodbox, which shouldn’t have been there.

The improvement of the given normal and depth are shown in the results part of Figure 6. Moreover, after 19 times of the iteration, a plausible depth map of foodbox-1-1-5 is obtained. Such results can validate the usefulness of the proposed method.

8 Discussion and conclusion

Limitations. Although the capability of our method has been proved robust in the previous experiments, it should be pointed out that like other Shape from shading method or shading based shape refinement solved by non-linear least square method, the most limited part is the time and space costs of the optimization . For example, a single iteration on $32\times 32$ images may take 10 minutes. Thus for an experiment on images of such size can take several hours.

Conclusion. In this paper, a new problem called ”Bimodal Stereo” is established, which literately stereopsis from two different views of different modality, aims to recover the global shape and estimating rigid transformation between RGB camera and depth camera simultaneously. To solve this problem, we proposed an iterative computational framework which regularizes SfS process with shape prior from depth map while estimates the pose parameters between recovered shape and rotated depth map. Quantitative experiments showed that through our method, both SfS results and depth map can be refined as well as the pose parameters are robustly estimated.