Smooth and Inextensible Surface Reconstruction

Although the strategem of maximizing the sum of depths $\sum_{i=1}^{n_c} \mu_i$ described in the previous section gives reasonable results, it is merely a heuristic, not based on any valid principle related to surface properties. We therefore consider next a new formulation based on the principle of surface inextensibility.

Let the surface be modelled as a function $\mathcal {W} : \mathbb{R}^2 \rightarrow \mathbb{R}^3$, mapping the planar template to $3$-dimensional space. The inextensibility constraint is equivalent to saying that the map $\mathcal {W}$ must be everywhere a local isometry. This condition may be expressed in terms of its Jacobian. Let $\mathsf{J}(\mathbf{q}) \in \mathbb{R}^{3 \times 2}$ be the Jacobian matrix $\partial \mathcal {W} / \partial \mathbf{q}$ evaluated at the point $\mathbf{q}$. The map  $\mathcal {W}$ is an isometry at $\mathbf{q}$ if the columns of $\mathsf{J}(\mathbf{q})$ are orthonormal. This local isometry can be enforced for the whole surface with the following least-squares constraint:

$$\iint \left\Vert \mathsf{J}(\mathbf{q})^\mathsf{T}\mathsf{J}(\mathbf{q}) - \mathsf{I}_2 \right\Vert^2 \mathrm d\mathbf{q} = 0.$$ (7.7)

In practice, we consider a discretization of the quantity in equation (7.7), namely

$$\mathcal {E}_i(\mathcal {W}) = \sum_{j=1}^{n_j} \left\Vert \maths... ...(\mathbf{g}_j)^\mathsf{T}\mathsf{J}(\mathbf{g}_j) - \mathsf{I}_2 \right\Vert^2,$$ (7.8)

where $\{\mathbf{g}_j\}_{j=1}^{n_j}$ is a set of 2D points in the template image space taken on a fine and regular grid (for instance, a grid of size $30 \times 30$). This term $\mathcal {E}_i(\mathcal {W})$ measures the departure from inextensibility of the surface $\mathcal {W}$.

Our minimization problem is then to minimize this quantity, over all possible surfaces, subject to the projection constraints, namely that point $\mathcal {W}(\mathbf{q}_i)$ projects to (or near to) the image point $\mathbf{q}_i'$, for all $i$.

Parametric Surface Model

The problem just described involves a minimization over all possible surfaces. Instead of considering this as a variational problem over all possible surfaces, we consider a parametrized family of surfaces. For this purpose, we chose Free-Form Deformations (FFD) (162) based on uniform cubic B-splines (61). All the details on this parametric model have been given in section 2.3.2. Here, we just precise the notation used in this chapter. Let  $\mathcal {W}_{\mbox{\boldmath ${\ell}$}} : \mathbb{R}^2 \rightarrow \mathbb{R}^3$ be the parametric FFD, parametrized by a family of 3D points ${\ell}$$_{jk}; j \in \{1,\ldots,n_u\}, k\in\{1, \ldots ,n_v\}$, which act as `attractors' for the surface.

For a point $\mathbf{q} = (u, v)$ in the template, the surface point is explicitly given as

$$\mathcal {W}_{\mbox{\boldmath${\ell}$}}(\mathbf{q}) = \sum_{j=1}^{n_u} \sum_{k=1}^{n_v} \mbox{\boldmath${\ell}$}_{jk} N_j(u) N_k(v).$$ (7.9)

The functions $N_j$ are the B-spline basis functions (61) which are polynomials of degree $3$. If point $\mathbf{q}_i = (u_i, v_i)$ is fixed and known then the surface point $\mathcal {W}_{\mbox{\boldmath ${\ell}$}}(\mathbf{q}_i)$ is expressed as a linear combination of the points ${\ell}$$_{jk}$, and hence can be written in the form $\mathcal {W}_{\mbox{\boldmath ${\ell}$}}(\mathbf{q}_i) = \mathsf{W}_i \mbox{\boldmath ${\ell}$}$, where $\mathsf{W}_i$ is a $3\times n_u n_v$ matrix depending only on the point $\mathbf{q}_i$, and ${\ell}$ is the vector obtained by concatenating all the points ${\ell}$$_{jk}$. Thus, the 3D point is a linear expression in terms of the parameter vector ${\ell}$. Since the polynomials $N_j$ and $N_k$ depend only on a local set of the attractor points ${\ell}$$_{jk}$, the matrix $\mathsf{W}_i$ is sparse, which is important for computational efficiency.

Surface Reconstruction as a Least-Squares Problem

By replacing $\mathbf{Q}_i$ by $\mathsf{W}_i$   ${\ell}$ in equation (7.6) we may arrive at a constraint:

$$\left\Vert \left( \left[\begin{array}{c} \mathbf{p}_1^\mathsf{... ...on_{\mathcal {I}} \mathbf{p}_3^\mathsf{T}\mathsf{W}_i \mbox{\boldmath${\ell}$}.$$ (7.10)

We may then formulate the optimization problem as minimizing the inextensibility cost $\mathcal {E}_i(\mathcal {W}_{\mbox{\boldmath ${\ell}$}})$ given in equation (7.8) over all choices of parameters ${\ell}$, subject to constraints equation (7.10). The constraints are SOCP constraints, but the cost function equation (7.8) is of higher degree in the parameters. To avoid the difficulties of constrained non-linear optimization, we choose a different course, by including the reprojection error into the cost function, leading to an unconstrained problem.

To simplify the formulation of the reprojection error, we introduce the depths $\mu_i$ as subsidiary variables, for reasons that become evident below. This is not strictly necessary, but reduces the degree of the reprojection-error term. The minimization problem now takes the form:

$$\min_{\mathbold{\mu}, \mbox{\boldmath${\ell}$}} \mathcal {E}_d(... ...}_i(\mbox{\boldmath${\ell}$}) + \beta \mathcal {E}_s(\mbox{\boldmath${\ell}$}),$$ (7.11)

where  $\mathcal{E}_d$, $\mathcal {E}_i$, $\mathcal{E}_s$ are the data (reprojection error), inextensibilty, and smoothing terms respectively. The data term ensures the consistency of the point correspondences with the reconstructed surface. $\mathcal {E}_i$ forces the inextensibility of the surface. $\mathcal {E}_b$ promotes smooth surface in order to cope with, for instance, lack of data. The relative influence of these three terms are controlled with the weights  $\alpha\in\mathbb{R}_+$ and  $\beta\in\mathbb{R}_+$. Note that the choice of $\alpha$ and $\beta$ is generally not critical.

The inextensibility term has been described previously. We now describe the two other terms in equation (7.11).

Data term.

Replacing $\mathbf{Q}_i$ by $\mathsf{W}_i$   ${\ell}$ in equation (7.5) gives an expression for the reprojection error associated with some point. However, the resulting expression is non-linear with respect to the parameters  ${\ell}$. We thus prefer a linear data term expressed in terms of `3D errors', which is the reason why we introduced the depths  $\mathbold{\mu}$ of the data points in the optimization problem. The data term is then defined by:

$$\mathcal {E}_d(\mathbold{\mu}, \mbox{\boldmath${\ell}$} ) = \sum_{i=1}^{n_c} \left\Vert \mathcal {W}_{\mbox{\boldmath${\ell}$}}(\mathbf{q}_i) - \mu_i \mathsf{P}^{-1}\bar{\mathbf{q}}_i' \right\Vert^2,$$ (7.12)

which measures the distance between the point $\mathcal {W}_{\mbox{\boldmath ${\ell}$}}$ on the surface and the point at depth $\mu_i$ along the ray defined by $\mathbf{q}_i'$.

Smoothing term.

In some cases, the point correspondences and the hypothesis of an inextensible surface are not sufficient. For instance, imagine that there is no point correspondence in a corner of the surface. In this case, there is nothing that indicates how the surface should behave. The corners of the surface can bend freely as long as they do not extend or shrink (like the corners of a piece of paper). To overcome this difficulty, we can add a third term (the smoothing term) in our cost function that favours non-bending surfaces. Note that usually, such terms are used to compensate for the undesirable effects of under-fitting and over-fitting. Doing so is usually a problem because it requires one to determine a correct value for the weight associated to the smoothing term (value $\beta$ in equation (7.11)). This is a sensible and critical way of balancing the effective complexity of the surface against the complexity of the data. Here, we do not have to care too much. Indeed, the complexity of the surface is limited by the fact that it is inextensible. Any small value (but big enough to be not negligible, for instance  $\beta = 10^{-4}$) is thus suitable for the weight of the smoothing term. We define our smoothing term using the bending energy:

$$\mathcal {E}_s(\mathbold{\mu}, \mbox{\boldmath${\ell}$} ) = \sum_{i=1}^3 \iint \left\Vert \frac{\partial^2 \mathcal {... ...q})}{\partial \mathbf{q}^2} \right\Vert^2_{\mathcal {F}} \mathrm d\mathbf{q}.$$ (7.13)

where $\mathcal {W}_{\mbox{\boldmath ${\ell}$}}^i(\mathbf{q})$ is the $i$-th coordinate of the point, and $\Vert\cdot\Vert _{\mathcal{F}}$ is the Frobenius norm of the Hessian matrix. With FFD, there exists a simple and linear closed-form expression for the bending energy:

$$\mathcal {E}_s( \mbox{\boldmath${\ell}$} ) = \Vert \mathsf{B}^{1/2} \mbox{\boldmath${\ell}$} \Vert^2 = \mbox{\boldmath${\ell}$} ^\mathsf{T}\mathsf{B} \mbox{\boldmath${\ell}$}$$ (7.14)

where $\mathsf{B}\in\mathbb{R}^{3p \times 3p}$ is a symmetric, positive, and semi-definite matrix which can be easily computed from the second derivatives of the B-spline basis functions.

Initial solution.

The problem of equation (7.11) is a non-linear least-squares minimization problem typically solved using an iterative scheme such as Levenberg-Marquardt. Such an algorithm requires a correct initial solution. We used an FFD surface fitted to the 3D points reconstructed with one of the point-wise methods presented in section 7.3. Subsequently, since we use a surface model which is linear with respect to its parameters, the initial parameters ${\ell}$ can be found by solving the least-squares problem:

$$\min_{\mbox{\boldmath${\ell}$}} \sum_{i=1}^{n_c} \left\Vert \ma... ... \left\Vert \mathsf{W}_i\mbox{\boldmath${\ell}$}- \mathbf{Q}_i \right\Vert^2 .$$ (7.15)

An alternative is to modify the problem (7.6), expressing $\mathbf{Q}_i$ in terms of the required parameters ${\ell}$, according to $\mathbf{Q}_i = \mathsf{W}_i$   ${\ell}$. Then one may solve for ${\ell}$ directly using SOCP. If necessary, the linear smoothing term of equation (7.13) can be included in equation (7.15).