# B-spline Shape from Motion & Shading: An Automatic Free-form Surface
Modeling for Face Reconstruction

Weilong Peng
School of Computer Science, Tianjin University
<EMAIL_ADDRESS>Zhiyong Feng
School of Computer Science, Tianjin University
<EMAIL_ADDRESS>Chao Xu
School of Software, Tianjin University
Tianjin, 300072, China
<EMAIL_ADDRESS>
###### Abstract

Recently, many methods have been proposed for face reconstruction from
multiple images, most of which involve fundamental principles of Shape from
Shading and Structure from motion. However, a majority of the methods just
generate discrete surface model of face. In this paper, B-spline Shape from
Motion and Shading (BsSfMS) is proposed to reconstruct continuous B-spline
surface for multi-view face images, according to an assumption that shading
and motion information in the images contain 1st- and 0th-order derivative of
B-spline face respectively. Face surface is expressed as a B-spline surface
that can be reconstructed by optimizing B-spline control points. Therefore,
normals and 3D feature points computed from shading and motion of images
respectively are used as the 1st- and 0th- order derivative information, to be
jointly applied in optimizing the B-spline face. Additionally, an IMLS
(iterative multi-least-square) algorithm is proposed to handle the difficult
control point optimization. Furthermore, synthetic samples and LFW dataset are
introduced and conducted to verify the proposed approach, and the experimental
results demonstrate the effectiveness with different poses, illuminations,
expressions etc., even with wild images.

##  1 Introduction

3D face has been widely used in areas of face recognition and verification,
expression recognition and facial animations. Face reconstruction is usually
based on a scene that the subject is fixed and scanned with special equipment
under a lab condition, video sequences or multi-view photographs. Recently,
more challenging works are reconstruction from wild images with unknown camera
calibration.

As a whole, there are two kinds of reconstruction technologies: high-cost
equipment relied technologies under lab condition and low-cost 2D image based
techniques out of lab. The former ones have implemented high quality face
reconstruction such as laser [23], stereo [2, 3, 24], structural lighting
[27], and light stages [9]. These technologies always rely on high resolution
camera, synchronized data or calibration. The latter kinds of reconstructions
are based on 2D informations, such as single-image method [8], video sequences
based method [19], multi-view photographs [16, 11], and even wild images based
methods [15, 18]. In contrast, the latter are low-cost and convenient to be
distributed. It is greatly significant but challengingly to research on face
reconstruction from images.

Generally, 2D methods involves several kinds of fundamental methods including
Shape from Shading (SFS) [28], Structure from Motion(SFM) [21], and 3D
Morphable Model (3DMM) [4, 5]. SFS assumes that the subject has a Lambertian
surface and relatively distant illumination, to compute surface normals. SFM
is based on assumption that it exists a coordinate transformation between
image coordinate system and camera coordinate system, to compute the positions
of surface points. And the idea of 3DMM is that human faces are within a
linear subspace, and that any novel face shape can be represented by linear
combination of shape eigenvectors deduced by PCA. SFS and SFM give the
geometrical and physical descriptions of face shape and imaging, and 3DMM
concentrates on the statistical explanation of 3D meshes or skeletons. We
believe face reconstruction is rather a geometrical optimization problem than
a statistical problem, as 3DMM is just an assist of geometrical method when
building detailed shape,eg.[25].

Motivated by the geometric and physical thought in SFS and SFM, we propose to
use B-spline surface to represent face shape, because it is a continual
surface model, while 3D mesh or skeleton model are just approximation of
continuous face by using discontinuous fragments. Deep in relation between
B-spline face and face images, we find that motion and shading in images are
from 0th- and 1st- order derivatives, respectively, of B-spline face surface,
by which a complete geometrical optimization for B-spline surface
reconstruction can be introduced: 1) the dense shading information can be
applied to optimize local B-spline shape; and 2) the sparse motion can be used
to optimize the global. During optimization, 2nd-order informations are also
used to guarantee a smooth warping and local information keeping. The
processes 1) and 2) iterate until convergence. Significantly, the final
B-spline face is high-order differential, and can approximate the true shape
in arbitrary precision.

Shape in Shading and Structure in Motion. SFS has been widely used for
reconstruction, e.g. single-view reconstruction [14], multiple frontal images
based reconstruction [22], and wild image based reconstruction [15, 18]. As
single-view is ill posed [17], a reference is always needed [14]. For wild
images, photometric stereo is applied to obtain accurate normals locally [15,
18]. SFM uses multiple frame or images to recover sparse 3D structure of
feature points of an object [21]. Spatial-transformation approach [20] only
estimates the depth of facial points. Bundle adjustment [1] fits the large
scale rigid object reconstruction, but it cannot generate the dense model of
non-rigid face. Incremental SFM approach [10] is proposed to build a 3D
generic face model for non-rigid face. The paper [18] optimizes the local
information with normals from shading, based on a 3D feature points-driven
global warping. Therefore, shading and motion are very important and different
geometric information of face, will enhance the reconstruction when they are
combined. In our method, they are integrated into B-spline face model as 1st-
and 0th- order derivatives of spline surface.

Face Surface Modeling. The surface modeling depends on the data input (point
cloud, noise, outlier, etc), output (point cloud, mesh, skeleton), and types
of shape(man-made shape, organic shape). Point cloud, skeleton and mesh grid
are the widely used man-made shape type for face reconstruction. Spatial
transformation method [20] estimates positions of sparse facial feature
points. Bundle adjustment [1] builds the dense point cloud for large scale
rigid object with a great number of images. Jingu et al [11] reconstruct face
dense mesh based on skeleton and 3DMM. Kemelmacher et al [15] apply
integration of normals to get discrete surface points. Roth et al [15]
reconstruct face based on Laplace mesh editing. The method of surface-
smoothness priors is also needed for reconstructing smooth mesh on point
cloud, e.g. radial basis function [7] and Poisson surface reconstruction [13].
Due to the fact that point cloud and 3D mess are discontinuous geometric
shape, they cannot approximate the true shape of face in arbitrary precision.
In contrast, our technology is based on a continuous free-form surface that
approximates a organic shape, and reconstruct B-spline face from images
directly, instead of intermediate point data. What’s more, because of B-spline
surface is a case of NURBS (Non-Uniform Rational B-Spline), the B-spline
reconstruction also can be imported to 3D modeling software like Rhino3D for
further editing, analysis, and transformation conveniently by adjusting the
B-spline control points.

In our work, initial local normals and motion estimation are obtained from
images by referring to the work [15], and they are updated during B-spline
optimization. The B-spline face is optimized and updated locally by using
normals and globally by 3D feature points respectively, and the process
iterates until convergence.

In summary, there are two primary contributions:

1\. The face is described as a B-spline face whose 1st- and 0th- order
derivatives can be discovered from shading and motion informations
respectively in the images; and iterative multi-least-square (IMLS)
optimization algorithm is also given.

2\. We introduce an automatic free-form surface modeling for face
reconstruction from multiple images with variation such as different poses,
illuminations, and expressions etc., even from images in the wild.

##  2 Methodology

The proposed method is operated on a set of face images with different poses,
illuminations, and expression, of an individual. In this section, face is
modeled as a uniform B-spline surface, whose 1st- and 0th- order derivatives
have intimate relations with normals in shading and 3D feature points from
motion respectively in the images in Section 2.1. Therefore, local and global
B-spline optimization by shading and motion are formulated in Section 2.2 and
2.3 severally. Finally, an IMLS optimization algorithm is given in Section
2.4.

###  2.1 Modeling B-spline Face

Assuming human face is a uniform B-spline surface, surface curvature and
normal at each point are described by basis function and control points, and
face can be reconstructed conveniently by computing the B-spline control
points.

We model face as a uniform B-spline surface F of degree $4\times 4$, with
$\emph{{B}}=\\{\emph{{b}}_{mn}\\}_{M\times N}$ as the control points, and
$U=\\{u_{m}\\}_{m=1}^{M+4},V=\\{v_{n}\\}_{n=1}^{N+4}$ as knots spliting $uv$
parameter plane in $u$ and $v$ direction. Then

$\textbf{{F}}(u,v)=\sum\limits_{m=1}^{M}{\sum\limits_{n=1}^{N}{R_{m,n}(u,v)\emph{{b}}_{mn}}},$
(1)

with $R_{m,n}(u,v)=N_{m,4}(u)\cdot N_{n,4}(v)$ and

$\left\\{{\begin{array}[]{*{20}{c}}{{N_{i,1}}(u)=\left\\{{\begin{array}[]{*{20}{l}}{\begin{array}[]{*{20}{c}}1&{{u_{i}}\leq
u<{u_{i+1}},}\end{array}}\\\
{\begin{array}[]{*{20}{c}}0&{otherwise,}\end{array}}\end{array}}\right.}\\\
{{N_{i,4}}(u)=\frac{{(u-{u_{i}})\cdot{N_{i,3}}(u)}}{{{u_{i+3}}-{u_{i}}}}+\frac{{({u_{i+4}}-u)\cdot{N_{i+1,3}}(u)}}{{{u_{i+4}}-{u_{i+1}}}}}.\end{array}}\right.$

Particularly, $\textbf{{F}}(u,v)$ is $C^{2}$-continuous, and z-component of a
B-spline face surface and its high order partial derivatives w.r.t. $u$ and
$v$ are shown by the maps in Fig.1. The continuous derivability means that F
can approximate the true shape in arbitrary $uv$ precision with deterministic
$k$-ordered partial derivative
$\frac{{{\partial^{k}}\textbf{{F}}}}{{\partial{u^{k}}}}$ and
$\frac{{{\partial^{k}}\textbf{{F}}}}{{\partial{v^{k}}}}$, $k=1,2$, and
$\frac{{{\partial^{2}}\textbf{{F}}}}{\partial u\partial v}$.

Figure 1: $C^{2}$ continuous characteristics of a $4\times 4$ order B-spline
face: maps of component z, $\frac{{\partial z}}{{\partial u}}$,
$\frac{{\partial z}}{{\partial v}}$, $\frac{{\partial^{2}z}}{{\partial
u^{2}}}$, $\frac{{\partial^{2}z}}{{\partial v^{2}}}$, and
$\frac{{\partial^{2}z}}{{\partial u}{\partial v}}$ of the face are shown.

####  2.1.1 Surface normals from 1st-order derivative

The 1st-order partial derivatives of F w.r.t $u$ and $v$ are

$\begin{array}[]{l}\textbf{{F}}_{u}^{\prime}(u,v)=\sum\limits_{m=1}^{M}{\sum\limits_{n=1}^{N}{N_{m,4}^{\prime}(u)\cdot
N_{n,4}(v)\textbf{{b}}_{mn}}}\\\
=\sum\limits_{m=1}^{M}{\sum\limits_{n=1}^{N}{(\frac{4}{{u_{m+4}-u_{i}}}N_{m,3}^{\prime}(u)-\frac{4}{{u_{m+5}-u_{m+1}}}N_{m+1,3}^{\prime}(u))\cdot
N_{n,4}(v)\textbf{b}_{mn}}}\\\ \end{array}$

and

$\begin{array}[]{l}\textbf{{F}}_{v}^{\prime}(u,v)=\sum\limits_{m=1}^{M}{\sum\limits_{n=1}^{N}{N_{m,4}(u)\cdot
N_{n,4}^{\prime}(v)\textbf{{b}}_{mn}}}\\\
=\sum\limits_{m=1}^{M}{\sum\limits_{n=1}^{N}{N_{m,4}(u)\cdot(\frac{4}{{u_{n+4}-u_{n}}}N_{n,3}^{\prime}(v)-\frac{4}{{u_{n+5}-u_{n+1}}}N_{n+1,3}^{\prime}(v))\textbf{b}_{mn}}}\\\
\end{array}$

respectively. And a more brief form of 1st-order derivatives at $(u,v)$ can be
written linearly as

$\left\\{{\begin{array}[]{*{20}{c}}{\textbf{{F}}{{}^{\prime}_{u}}{|_{u,v}}={\textbf{{T}}_{1}}{|_{u,v}}\cdot\textbf{{b}}}\\\
{\textbf{{F}}{{}^{\prime}_{v}}{|_{u,v}}={\textbf{{T}}_{2}}{|_{u,v}}\cdot\textbf{{b}}}\end{array}},\right.$
(2)

where vector $\textbf{\emph{b}}_{3MN\times 1}$ represents B-spline control
points, and $\textbf{\emph{T}}_{1}|_{u,v}$ and $\textbf{\emph{T}}_{2}|_{u,v}$
are 3-by-$3MN$ matrixes stacking the 1st-order coefficient of control points,
being referred as 1st-order coefficient matrixes.

Therefore, the surface normal vector at $(u,v)$ can be computed by the cross
product

$\emph{{H}}|_{u,v}=\frac{{\textbf{\emph{F}}_{u}^{\prime}|_{u,v}\times\textbf{\emph{F}}_{v}^{\prime}|_{u,v}}}{{\left\|{\textbf{\emph{F}}_{u}^{\prime}|_{u,v}\times\textbf{\emph{F}}_{v}^{\prime}|_{u,v}}\right\|}}=s|_{u,v}\cdot\textbf{\emph{F}}_{u}^{\prime}|_{u,v}\times\textbf{\emph{F}}_{v}^{\prime}|_{u,v}$
(3)

####  2.1.2 3D structure by 0th-order derivative

Particularly, 3D structure defined by $f$ feature points on the surface, as
shown in Fig.2, is given intuitively as

$\begin{array}[]{*{20}c}F(\emph{{B}},u(i),v(i))=\\\
\begin{array}[]{c}\left[{\begin{array}[]{*{20}c}{N_{m-3,3}(u(i))}\\\
{N_{m-2,3}(u(i))}\\\ {N_{m-1,3}(u(i))}\\\ {N_{m,3}(u(i))}\\\
\end{array}}\right]^{T}\cdot\emph{{B}}_{m,n}^{(4,4)}\cdot\left[{\begin{array}[]{*{20}c}{N_{n-3,3}(v(i))}\\\
{N_{n-2,3}(v(i))}\\\ {N_{n-1,3}(v(i))}\\\ {N_{n,3}(v(i))}\\\
\end{array}}\right],\\\ \end{array}\\\ \end{array}$ (4)

$i=1,2,...,f,$ and

$\emph{{B}}_{m,n}^{(4,4)}=\left[{\begin{array}[]{*{20}c}{\textbf{\emph{b}}_{m-3,n-3}}&{\textbf{\emph{b}}_{m-3,n-2}}&{\textbf{\emph{b}}_{m-3,n-1}}&{\textbf{\emph{b}}_{m-3,n}}\\\
{\textbf{\emph{b}}_{m-2,n-3}}&{\textbf{\emph{b}}_{m-2,n-2}}&{\textbf{\emph{b}}_{m-2,n-1}}&{\textbf{\emph{b}}_{m-2,n}}\\\
{\textbf{\emph{b}}_{m-1,n-3}}&{\textbf{\emph{b}}_{m-1,n-2}}&{\textbf{\emph{b}}_{m-1,n-1}}&{\textbf{\emph{b}}_{m-1,n}}\\\
{\textbf{\emph{b}}_{m,n-3}}&{\textbf{\emph{b}}_{m,n-2}}&{\textbf{\emph{b}}_{m,n-1}}&{\textbf{\emph{b}}_{m,n}}\\\
\end{array}}\right].$

Figure 2:  Face structure defined by 40 feature points: the left side shows
the point positions in a face image; the right side shows the structure
topology with eye center points of $O_{1}$ (-25, 0, 0) and $O_{2}$ (25, 0, 0)
in 3$D$ space, which looks like a frontal 2D face structure from the direction
of normal (0, 0, 1).

Equ.4 can be rewritten linearly as

$\textbf{\emph{F}}|_{u(i),v(i)}=\textbf{\emph{T}}|_{u(i),v(i)}\cdot\textbf{\emph{b}},i=1,2,...,f,$
(5)

where $\textbf{\emph{T}}|_{u(i),v(i)}$ is a 3-by-$3MN$ 0th-order coefficient
matrix. $\textbf{\emph{T}}|_{u(i),v(i)}$ is also a sparse matrix as every
feature point is determined by only 16 control points according to
$\emph{{B}}_{m,n}^{(4,4)}$.

Surface normals are the geometric details of a surface, and structure is a
sparse but global shape of a face. They are 1st- and 0th- order information
that can be represented by control parameter b, seen in Equ.2, Equ.3 and
Equ.5. Conversely, B-spline parameter b can be rebuilt from normals and
structure, to build the face surface of an individual.

###  2.2 Local B-spline Shape From Shading

The shading in images covers normals and albedo of a surface, which can be
used to optimize B-spline control points of a B-spline face. Firstly,
photometric stereo method is used to estimate the normal information by
assuming that face images are captured from a B-spline shape by a camera model
of scaled orthographic projection, relatively distant illumination and
Lambertian reflectance. Secondly, a normals driven control points optimization
is proposed to optimize the local shape of B-spline face.

####  2.2.1 Photometric stereo

Photometric stereo is also based on assumption that a set of images of a same
person are captured under different illumination conditions, and that the
shape of each point can be solved by the observed variation in shading of the
images.

Being faced with the variations in the wild photos, it is difficult to define
the ”shape” of a person’s face, for which the normal information is also
undetermined. Our solution is to solve a face shape that is locally similar to
as many photos as possible from the normalized wild images as done in the
paper [15]. Firstly, $n$ images are normalized and clipped to frontal-posed
with $p$ pixels, and input into $M_{n\times p}$. Secondly, we estimate the
initial shape $\tilde{S}$ and lighting $\tilde{L}$ by factorizing $M=LS$ via
SVD[26]. $\tilde{L}=U\sqrt{D}$ and $\tilde{S}=\sqrt{D}V^{T}$, where
$M=UD{V^{T}}$.

To approach the true normal information, we estimate the shape $S$ and
ambiguity $A$ by following the approach in paper [15]. Firstly, we update
$\tilde{S}$ firstly by minimizing
${\min_{\tilde{S}}}\left\|{M-\tilde{L}\tilde{S}}\right\|+{S^{T}}GS$ to get the
initial normal $S$, where $G=diag(-1,1,1,1)$, with initial $\tilde{L}$.
Secondly, we solve ambiguity $A$ by minimizing
${\min_{A}}{\left\|{S^{t}-A\tilde{S}}\right\|^{2}}$, where $S^{t}$ is the
normal shape of a template and it is consistent with $S$ at level of pixel
position. At last, the surface normals are obtained by $S=A\tilde{S}$.

Particularly, the consistent normal shape $S^{t}$ are from the normals of
template
$\frac{{\textbf{\emph{F}}_{u}^{\prime}\times\textbf{\emph{F}}_{v}^{\prime}}}{{\left\|{\textbf{\emph{F}}_{u}^{\prime}\times\textbf{\emph{F}}_{v}^{\prime}}\right\|}}=\Lambda\cdot\textbf{\emph{F}}_{u}^{\prime}\otimes\textbf{\emph{F}}_{v}^{\prime}$,
where symbol $\otimes$ is compound operator of cross-product $\times$ in
Equ.3.

The definition of $\otimes$ operation is
$\textbf{{w}}\otimes\textbf{{v}}=\left[{\begin{array}[]{*{20}{c}}{{\textbf{{w}}_{1}}\times{\textbf{{v}}_{1}}}\\\
{{\textbf{{w}}_{2}}\times{\textbf{{v}}_{2}}}\\\ \vdots\\\
{{\textbf{{w}}_{p}}\times{\textbf{{v}}_{p}}}\end{array}}\right]$, where w and
v are $3p\times 1$ vectors containing $p$ normals.

Therefore, the vecorization of $S^{t}$ is

$\textbf{{h}}_{t}=\Lambda|_{\textbf{{b}}_{t}}\cdot(\textbf{\emph{F}}_{u}^{\prime}\otimes\textbf{\emph{F}}_{v}^{\prime})^{(sn)}=\Lambda|_{\textbf{{b}}_{t}}\cdot((\textbf{\emph{T}}_{1}^{(sn)}\textbf{\emph{b}}_{t})\otimes(\textbf{\emph{T}}_{2}^{(sn)}\textbf{\emph{b}}_{t}))$

, where vector $\textbf{{h}}_{t}$ is $3p\times 1$, and $(*)^{(sn)}$ is a
selection operator that selects $3p$ rows of 1st-order coefficients
corresponding to $p$ normals consistent with $S$; and $\Lambda$, whose value
depends on $\textbf{{b}}_{t}$ is a $3p\times 3p$ diagonal matrix that stores
$3p$ reciprocals of lengths of the $p$ selected normals; and
$\textbf{\emph{b}}_{t}$ represent the control points of B-spline template
face.

Similarly, the obtained $S$ can also be vectorized as h that is $3p\times 1$.
Then h can be used to optimize control points b for B-spline objective face
conversely, which will be introduced in the following subsection.

####  2.2.2 Local B-spline optimization

According to the relation of normals vector h and B-spline control points b,
B-spline face can be solved by

$\mathop{\min}\limits_{\textbf{{b}}}\textbf{O}_{1}=\mathop{\min}\limits_{\textbf{{b}}}\left\|{\textbf{{h}}-\Lambda|_{\textbf{{b}}}\cdot((\textbf{{T}}_{1}^{(sn)}\textbf{{b}})\otimes((\textbf{{T}}_{2}^{(sn)}\textbf{{b}}))}\right\|,$

which is also referred to as 1st-order optimization.

However, there exists two issues: 1) the low-dimension h may not guarantee an
unique solution of high-dimension b; and 2) the system is not simply linear,
which is difficult to be solved.

Therefore, a frontal constraint with $\textbf{{b}}_{t}$ is applied to make a
unique solution; And a strategy of approximating to linearization is also
proposed to make a linear solution.

Frontal Constraint. As the normals are estimated actually from frontal-face
images, a frontal constraint will enhance an effective solution of 1st-order
optimization surely by

${\textbf{O}_{2}}=\left\|{\textbf{{T}}^{(sxy)}\textbf{{b}}-\textbf{{T}}^{(sxy)}{\textbf{{b}}_{t}}}\right\|<\epsilon,$

where matrix T is the stacked coefficients of control points, and
$(*)^{(sxy)}$ is a selection operator that sets the coefficients corresponding
to z- components of control points to zeros, retaining the coefficients
corresponding to x- and y- components.

This constraint is also referred as 0th-order regularization for local
B-spline optimization as following,

$\mathop{\min}\limits_{\textbf{{b}}}{\textbf{{O}}_{1}},\ \ \ with\ \ \
{\textbf{{O}}_{2}}<\epsilon.$

Particularly, the item $\textbf{{O}}_{2}$ is a linear form, but the first item
$\textbf{{O}}_{1}$ is not a simple linear form. Therefore, an approximating to
linearization is proposed.

Approximating to Linearization. According to the characteristics of the cross-
product “$\otimes$”, the first item $\textbf{{O}}_{1}$ can be rewritten as a
linear-like formulation:

$\left\|{\textbf{{h}}-{\boldsymbol{L}}{|_{\textbf{{b}}}}\cdot\textbf{{b}}}\right\|\
\ or\ \
\left\|{\textbf{{h}}-{\boldsymbol{R}}{|_{\textbf{{b}}}}\cdot\textbf{{b}}}\right\|,$

where

${\boldsymbol{L}}{|_{\textbf{{b}}}}=\Lambda{|_{\textbf{{b}}}}\cdot{\left[{\textbf{{T}}_{1}^{(sn)}\textbf{{b}}}\right]_{\otimes}}\cdot\textbf{{T}}_{2}^{(sn)}$

and

$\textbf{{R}}{|_{\textbf{{b}}}}=-\Lambda{|_{\textbf{{b}}}}\cdot{\left[{\textbf{{T}}_{2}^{(sn)}\textbf{{b}}}\right]_{\otimes}}\cdot\textbf{{T}}_{1}^{(sn)}.$

Particularly, the operation $[*]_{\otimes}$ makes vector
$\textbf{{w}}_{3p\times 1}$ a $3p\times 3p$ matrix
$\left[\textbf{{w}}\right]_{\otimes}=diag(\left[{{\textbf{{w}}_{1}}}\right]_{\times},\left[{{\textbf{{w}}_{2}}}\right]_{\times},...,\left[{{\textbf{{w}}_{p}}}\right]_{\times})$,
where
$\textbf{{w}}=\left[\textbf{{w}}_{1};\textbf{{w}}_{2};...;\textbf{{w}}_{p}\right]$,
and

$\left[{{\textbf{{w}}_{i}}}\right]_{\times}=\left[{\begin{array}[]{*{20}{c}}0&{-{\textbf{{w}}_{i}^{z}}}&{{\textbf{{w}}_{i}^{y}}}\\\
{{\textbf{{w}}_{i}^{z}}}&0&{-{\textbf{{w}}_{i}^{x}}}\\\
{-{\textbf{{w}}_{i}^{y}}}&{{\textbf{{w}}_{i}^{x}}}&0\end{array}}\right],i=1,2,...p.$

If b is a known parameter, e.g as $\textbf{{b}}_{0}$, for
${\boldsymbol{L}}|_{\textbf{{b}}}$ and ${\boldsymbol{R}}|_{\textbf{{b}}}$, the
minimization of both
$\left\|{\textbf{{h}}-{\boldsymbol{L}}{|_{\textbf{{b}}}}_{0}\cdot\textbf{{b}}}\right\|\
\ and\ \
\left\|{\textbf{{h}}-{\boldsymbol{R}}{|_{\textbf{{b}}}}_{0}\cdot\textbf{{b}}}\right\|$
will be a linear system.

In fact, formulation
$\left\|{\textbf{{h}}-{\boldsymbol{L}}{|_{\textbf{{b}}}}_{0}\cdot\textbf{{b}}}\right\|$
is used to optimize the control points in parameter space of $v$ by fixing
$u$, while
$\left\|{\textbf{{h}}-{\boldsymbol{R}}{|_{\textbf{{b}}}}_{0}\cdot\textbf{{b}}}\right\|$
is used to optimize the control points in parameter space of $u$ by fixing
$v$.

Figure 3:  Iterative adjustment on two partial derivatives: Process (1) to (2)
adjusts $\frac{{\partial\textbf{{F}}}}{{\partial u}}$ by fixing
$\frac{{\partial\textbf{{F}}}}{{\partial v}}$, and process (3) to (4) adjusts
$\frac{{\partial\textbf{{F}}}}{{\partial v}}$ by fixing
$\frac{{\partial\textbf{{F}}}}{{\partial u}}$, … until that
$\frac{{\partial\textbf{{F}}}}{{\partial
u}}\times\frac{{\partial\textbf{{F}}}}{{\partial v}}$ is infinitely close to
objective n; Process A implements a practically and iteratively linear handle
for B-spline surface adjustment in B.

A practical skill is to optimize the control points on the two parameter
spaces alternately and iteratively. The two iteration items for
$\textbf{O}_{1}$ can be rewritten as

$\left\|{\textbf{{h}}-{\boldsymbol{L}}{|_{{\textbf{{b}}_{0}}}}\cdot\textbf{{b}}}\right\|+\lambda\left\|{\textbf{{T}}_{1}^{(sn)}\cdot\textbf{{b}}-\textbf{{T}}_{1}^{(sn)}\cdot{\textbf{{b}}_{t}}}\right\|$

and

$\left\|{\textbf{{h}}-{\boldsymbol{R}}{|_{{\textbf{{b}}_{0}}}}\cdot\textbf{{b}}}\right\|+\lambda\left\|{\textbf{{T}}_{2}^{(sn)}\cdot\textbf{{b}}-\textbf{{T}}_{2}^{(sn)}\cdot{\textbf{{b}}_{t}}}\right\|.$

As shown in Fig.3, two partial derivatives
$\frac{{\partial\textbf{{F}}}}{{\partial v}}$ and
$\frac{{\partial\textbf{{F}}}}{{\partial u}}$ at $(u,v)$ are updated until
$\frac{{\partial\textbf{{F}}}}{{\partial
v}}\times\frac{{\partial\textbf{{F}}}}{{\partial u}}$ converges to n. Finally,
the whole surface approaches to the optimal.

By integrating with $\textbf{O}_{2}$, the final formulation for local
optimization is the two iteration items as following:

$\footnotesize\mathop{\min}\limits_{\textbf{{b}}}\left\|{\left[{\begin{array}[]{*{20}{c}}\textbf{{h}}\\\
{{\textbf{{T}}^{(sxy)}}{\textbf{{b}}_{t}}}\end{array}}\right]-\left[{\begin{array}[]{*{20}{c}}{{\boldsymbol{L}}{|_{{\textbf{{b}}_{0}}}}}\\\
{{\textbf{{T}}^{(sxy)}}}\end{array}}\right]\textbf{{b}}}\right\|+\lambda\left\|{\textbf{{T}}_{1}^{(sn)}\textbf{{b}}-\textbf{{T}}_{1}^{(sn)}{\textbf{{b}}_{t}}}\right\|$
(6)

and

$\footnotesize\mathop{\min}\limits_{\textbf{{b}}}\left\|{\left[{\begin{array}[]{*{20}{c}}\textbf{{h}}\\\
{{\textbf{{T}}^{(sxy)}}{\textbf{{b}}_{t}}}\end{array}}\right]-\left[{\begin{array}[]{*{20}{c}}{{\boldsymbol{R}}{|_{{b_{0}}}}}\\\
{{\textbf{{T}}^{(sxy)}}}\end{array}}\right]\textbf{{b}}}\right\|+\lambda\left\|{\textbf{{T}}_{2}^{(sn)}\textbf{{b}}-\textbf{{T}}_{2}^{(sn)}{\textbf{{b}}_{t}}}\right\|$
(7)

The $\textbf{{b}}_{0}$ is inited by value of $\textbf{{b}}_{t}$. Then we can
solve b and update $\textbf{{b}}_{0}$ orderly by Minimization (6) and (7)
iteratively until convergence. Please refer to line 4\- 9 in Algorithm 1 for
details.

###  2.3 Global B-spline Shape From Motion

The motion informations across images cover the geometry projections of
feature points of a 3D object on 2D projection planes, which also can be used
to optimize B-spline control points of a B-spline face. What’s important, the
motion informations feature their sparsity and globality. Firstly, we use SVD
method to re-estimate the sparse motion information by assuming that face
images are captured from a B-spline shape by a camera model of a scaled
orthographic projection. Secondly, the re-estimated motion information is used
to optimize the shape of B-spline face globally. The two steps also are
conducted iteratively.

####  2.3.1 Re-estimating the Motion

According to relation between facial points $q$ on image and $Q$ on B-spline
face: $q=sRQ+t$, pose parameters $s$,$R$ and $t$ for each image are recovered
based on the linear transformation and SVD method [15].

Note a vector $\textbf{{f}}_{2nf\times 1}$ stacking $f$ feature points $(q-t)$
of $n$ images as the the re-estimated 2D positions, and a projection matrix
$\textbf{{P}}_{2nf\times 3f}$ stacking $n$ views of parameters $sR$.

####  2.3.2 Global B-spline Warping

According to the projection relation between 2D feature points and and 3D
feature points, warping is based on an optimized B-spline face
$\textbf{{b}}_{0}$ that contains local information of objective face:

$\mathop{\min}\limits_{\textbf{{b}}}\left\|{\textbf{{f}}-\textbf{{P}}\cdot{\textbf{{T}}^{(sf)}}\textbf{{b}}}\right\|+\zeta\cdot
l(\textbf{{b}},\textbf{{b}}_{0}),$

where $(*)^{(sf)}$ is a selection operator that selects $3f$ rows of
coefficients corresponding to $f$ feature points, and $\textbf{{T}}^{(sf)}$ is
a $3f\times 3MN$ matrix that stacks $\textbf{{T}}|_{u(i),v(i)},i=1,2,...,f$.
Regularity item $l(\textbf{{b}},\textbf{{b}}_{0})$ is the distance of local
information between faces b and $\textbf{{b}}_{0}$, and $\zeta$ is a weight
factor. The first item is also referred to 0st-order item.

To eliminate affect of geometric rotation brought by 0st-order warping which
changes the local normals greatly, and to make sure a smoothness of shape, we
regularize the 0st-order item by using 2nd-order derivative information.

2nd-order Regularization. According to Equ.2 that represents 1st-order partial
derivative as
$\textbf{{F}}{{}^{\prime}_{u}}={\textbf{{T}}_{1}}\cdot\textbf{{b}}$, the 2nd-
order partial derivative can be introduced as

$\left\\{{\begin{array}[]{*{20}{c}}{\textbf{{F}}^{\prime\prime}_{uu}={\textbf{{T}}_{11}}\cdot\textbf{{b}}}\\\
{\textbf{{F}}^{\prime\prime}_{vv}={\textbf{{T}}_{22}}\cdot\textbf{{b}}}\\\
{\textbf{{F}}^{\prime\prime}_{uv}={\textbf{{T}}_{12}}\cdot\textbf{{b}}}\end{array}}\right.,$

in a similar way. The regularization item representing the 2nd-order
approximation between template and objective face is given as:

$l(\textbf{{b}},\textbf{{b}}_{0})=\left\|{\left[{\begin{array}[]{*{20}{c}}{{\textbf{{T}}_{11}}}\\\
{{\textbf{{T}}_{22}}}\\\
{{\textbf{{T}}_{12}}}\end{array}}\right]\textbf{{b}}-\left[{\begin{array}[]{*{20}{c}}{{\textbf{{T}}_{11}}}\\\
{{\textbf{{T}}_{22}}}\\\
{{\textbf{{T}}_{12}}}\end{array}}\right]{\textbf{{b}}_{0}}}\right\|,$

which is also referred to as 2nd-order regularization.

With the 0th- and 2nd- order information, a globally warping of B-spline face
can be formulated as

$\mathop{\min}\limits_{\textbf{{b}}}\left\|{\textbf{{f}}-\textbf{{P}}\cdot{\textbf{{T}}^{(sf)}}\textbf{{b}}}\right\|+\zeta\left\|{\left[{\begin{array}[]{*{20}{c}}{{\textbf{{T}}_{11}}}\\\
{{\textbf{{T}}_{22}}}\\\
{{\textbf{{T}}_{12}}}\end{array}}\right]\textbf{{b}}-\left[{\begin{array}[]{*{20}{c}}{{\textbf{{T}}_{11}}}\\\
{{\textbf{{T}}_{22}}}\\\
{{\textbf{{T}}_{12}}}\end{array}}\right]{\textbf{{b}}_{0}}}\right\|$ (8)

Motion re-estimation and B-spline warping are also carried on iteratively,
during which 2D position f, projection P, objective b and template
$\textbf{{b}}_{t}$ are updated in turn, until that template $\textbf{{b}}_{t}$
is converged. Please refer to line 10-14 in Algorithm 1 for details.

###  2.4 Algorithm

Algorithm 1 IMLS Algorithm for B-spline face optimization from Shading &
Motion

0: A set of face images, B-spline template face $\textbf{{b}}_{t}$, and $uv$
parameters $(u_{i},v_{i})$, $i=1,2,...,f$ for its feature points.

1: Detect the feature points of images

2: while b is not converged do

3: Let $\textbf{{b}}_{0}=\textbf{{b}}_{t}$. // Inited with template// Compute
local B-spline from shading

4: while $\textbf{{b}}_{0}$ is not converged do

5: Conduct photometric stereo to obtain normal vector h using template
$\textbf{{b}}_{0}$. // Optimize shape in parameter space of $v$

6: With known $\textbf{{b}}_{0}$, solve Minimization (6) by
$\textbf{{b}}^{(1)}={({{\boldsymbol{G}}^{T}}{\boldsymbol{G}}+\lambda{\left[{\textbf{{T}}_{1}^{(sn)}}\right]^{T}}\textbf{{T}}_{1}^{(sn)})^{-1}}\cdot({{\boldsymbol{G}}^{T}}\textbf{{h}}+\lambda{\left[{\textbf{{T}}_{1}^{(sn)}}\right]^{T}}{\textbf{{T}}_{1}^{(sn)}}{\textbf{{b}}_{0}}),$
where
${\boldsymbol{G}}=\left[{\begin{array}[]{*{20}{c}}{{\boldsymbol{L}}{|_{{\textbf{{b}}_{0}}}}}\\\
{{\textbf{{T}}^{(sxy)}}}\end{array}}\right]$. // Optimize shape in parameter
space of $u$

7: Let $\textbf{{b}}_{0}=\textbf{{b}}^{(1)}$, and solve Minimization (7) by
$\textbf{{b}}^{(2)}={({{\boldsymbol{H}}^{T}}{\boldsymbol{H}}+\lambda{\left[{\textbf{{T}}_{2}^{(sn)}}\right]^{T}}T_{2}^{(sn)})^{-1}}\cdot({{\boldsymbol{H}}^{T}}\textbf{{h}}+\lambda{\left[{\textbf{{T}}_{2}^{(sn)}}\right]^{T}}{\textbf{{T}}_{2}^{(sn)}}{\textbf{{b}}_{0}}),$
where
${\boldsymbol{H}}=\left[{\begin{array}[]{*{20}{c}}{{\boldsymbol{R}}{|_{{b_{0}}}}}\\\
{{\textbf{{T}}^{(sxy)}}}\end{array}}\right]$.

8: Let $\textbf{{b}}_{0}=\textbf{{b}}^{(2)}$. // Update local shape.

9: end while // Warp B-spline globally via motion

10: while $\textbf{{b}}_{t}$ is not converged do

11: Re-estimate motion information f and P using the B-spline template
$\textbf{{b}}_{t}$.

12: With known $\textbf{{b}}_{0}$, solve Minimization (8) by:
$\textbf{{b}}={({(\textbf{{P}}{\textbf{{T}}^{(sf)}})^{T}}\textbf{{P}}{\textbf{{T}}^{(sf)}}+\zeta{\boldsymbol{K}^{T}}\boldsymbol{K})^{-1}}\cdot({(\textbf{{P}}{\textbf{{T}}^{(sf)}})^{T}}f+\zeta{\boldsymbol{K}^{T}}\boldsymbol{K}{\textbf{{b}}_{0}}),$
where
${\boldsymbol{K}}=\left[{\begin{array}[]{*{20}{c}}{{\textbf{{T}}_{11}}}\\\
{{\textbf{{T}}_{22}}}\\\ {{\textbf{{T}}_{12}}}\end{array}}\right]$.

13: Let $\textbf{{b}}_{t}=\textbf{{b}}$. // Update template

14: end while

15: end while

15: Solution of B-spline objective face b.

An IMLS (iterative multi-least-square) algorithm is proposed to optimize the
B-spline face, seen in Algorithm 1. Processes of computing local B-spline
shape from shading and warping B-spline globally via motion are carried on
iteratively to get the final surface, differing from the work in [18] where
the uses of shading and motion are in two separate processes. Particularly,
Minimization (6) and Minimization (7) during local optimization, and
Minimization (8) during global optimization are linear system with closed-from
solutions of least square, seen in line 6,7,12 in Algorithm 1.

Figure 4:  Visual reconstructions and comparisons for sessions S1, S2, S3, S4
and S5: for each session of reconstructions, the left shows the results by
reimplemented method[15], and the right shows our results. Figure 5: Sample
data simulated by the spaces-times faces[27] : images and 3D model with
various poses and illuminations are available; samples S1,S2,S3, and S4 are
used for evaluation.  Figure 6: Visual reconstructions and comparisons for
Arnold Schwarzenegger, Colin Powell, Donald Rumsfeld, George W.Bush, Hugo
Chavez and Gloria Macapagal Arroyo in LFW [12]: 1) Images in column 1 and 8
are mean textures and true images respectively; 2) Columns 2, 4, and 6 show
the results by the reimplemented method[15], and columns 3, 5, and 7 show our
results; 3) columns 6 and 7 are rendered with mean texture.

##  3 Experiment

In this section experiments is resented to verify our method of MsSfMS. We
first describe the pipeline to prepare a collection of face images of a same
person for B-spline face reconstruction. Then quantitative and qualitative
comparisons are demonstrated with baseline method on ground truth data [27]
with various expressions, illuminations and poses. Finally, challenging
reconstructions are also conducted and compared based on real unconstrained
data taken from the Labeled Faces in Wild (LFW) [12] database.

###  3.1 Pipeline of experimental data

Ground truth data with expression The ground truth data is from the space-
times faces [27] which contains 3D face models with different expressions. And
different poses and illuminations can also be simulated by the spaces-times
faces, seen in Fig.5. Images with various poses and illuminations are
collected, and feature points manually labeled. The reconstruction is
evaluated by the error to the ground truth model.

LFW data Images of each person are collected and input into a facial point
detector [6] that has a similar high performance to human, to find the 40
facial points shown in Fig.1. The initial B-spline template face is computed
from a neutral model of space-time faces.

###  3.2 Standard images

We conduct 5 session of reconstructions: the first four is reconstruction
expression S1, S2, S3 and S4 by using their corresponding images, and the
firth reconstruction S5 is based on images with different expressions. Each
session contains 40 images with various illumination and different poses.
Reconstruction result is compared with the reimplemented method[15].

To compare the approaches numerically, we compute the shortest point-to-point
distance from ground truth to reconstructed surface. Point clouds are sampled
from B-spline face and aligned according to absolute orientation problem. Mean
euclidean distance and the root mean square (RMS) of the distances, after
normalized by the eye-to-eye distance, are reported in Table.1. Visual
comparison is shown in Fig.4. Our results demonstrate promise performance. An
important fact is that method[15] cannot make a credible depth information and
global shape, but our method solves the problem by local optimization and
global warping.

Table 1: Distances of the reconstruction to the ground truth Method | Index | S1 | S2 | S3 | S4 | S5
---|---|---|---|---|---|---
Method | Mean(%) | 8.53 | 8.67 | 8.92 | 10.49 | 9.29
[15] | RMS(%) | 6.92 | 7.34 | 4.68 | 6.21 | 7.58
Ours | Mean(%) | 7.09 | 7.86 | 7.11 | 8.92 | 8.49
| RMS(%) | 4.57 | 5.61 | 4.18 | 5.61 | 5.60

###  3.3 Real unconstrained data

Our method is also tested based on real unconstrained data. Unlike the
experiment in paper[15] using hundreds of images, we conduct reconstruction
with limited number of images, as so many images are not always available for
some task like criminal investigation and small simple size problem. 35 images
are collected from LFW for each of six persons (Arnold Schwarzenegger, Colin
Powell, Donald Rumsfeld, George W.Bush, Hugo Chavez and Gloria Macapagal
Arroyo), covering different pose and expressions. Reconstructions and
comparisons are shown in Fig.6. Our method steadily produce good looking model
from different viewpoints. Especially for case of Gloria Macapagal Arroyo,
while the method[15] is not working, our method generates satisfactory result.
It is convinced that the reconstructions express the individual face structure
well and fit well with pose-variety photos.

##  4 Conclusions

We presented an B-spline free-form surface modeling for face reconstruction,by
leveraging the spline cues of motion and shading in a set of face images. The
method works well for data with different poses and expressions, even for wild
cases. The key contributions are that the proposed $C^{2}$ B-spline shape
model describes the motion and shading geometrically in 0th- and 1st- order
derivatives respectively, and that an optimization algorithm is also supplied
for reconstruction. In the future, the model also can be extented for other
research areas, such as expression tracking and facial animation by adjusting
the control points locally, and recognition and verification based on local
B-spline shape information.

## References

  * [1] S. Agarwal, Y. Furukawa, N. Snavely, I. Simon, B. Curless, S. M. Seitz, and R. Szeliski.  Building rome in a day.  Communications of the ACM, 54(10):105–112, 2011.
  * [2] T. Beeler, B. Bickel, P. Beardsley, B. Sumner, and M. Gross.  High-quality single-shot capture of facial geometry.  ACM Transactions on Graphics (TOG), 29(4):40, 2010.
  * [3] T. Beeler, F. Hahn, D. Bradley, B. Bickel, P. Beardsley, C. Gotsman, R. W. Sumner, and M. Gross.  High-quality passive facial performance capture using anchor frames.  In ACM Transactions on Graphics (TOG), volume 30, page 75. ACM, 2011\.
  * [4] V. Blanz, A. Mehl, T. Vetter, and H. P. Seidel.  A statistical method for robust 3d surface reconstruction from sparse data.  In 3D Data Processing, Visualization and Transmission, 2004. 3DPVT 2004. Proceedings. 2nd International Symposium on, pages 293–300, 2004\.
  * [5] V. Blanz and T. Vetter.  A morphable model for the synthesis of 3d faces.  In Proceedings of the 26th annual conference on Computer graphics and interactive techniques, pages 187–194, 1999.
  * [6] X. P. Burgos-Artizzu, P. Perona, and P. Dollár.  Robust face landmark estimation under occlusion.  In Computer Vision (ICCV), 2013 IEEE International Conference on, pages 1513–1520. IEEE, 2013.
  * [7] J. C. Carr, R. K. Beatson, J. B. Cherrie, T. J. Mitchell, W. R. Fright, B. C. Mccallum, and T. R. Evans.  Reconstruction and representation of 3d objects with radial basis functions.  Computer Graphics Pages–76 Acm Siggraph, pages 67–76, 2001.
  * [8] J. Fortuna and A. M. Martinez.  Rigid structure from motion from a blind source separation perspective.  International journal of computer vision, 88(3):404–424, 2010.
  * [9] A. Ghosh, G. Fyffe, B. Tunwattanapong, J. Busch, X. Yu, and P. Debevec.  Multiview face capture using polarized spherical gradient illumination.  ACM Transactions on Graphics (TOG), 30(6):129, 2011.
  * [10] J. Gonzalez-Mora, F. De la Torre, N. Guil, and E. L. Zapata.  Learning a generic 3d face model from 2d image databases using incremental structure-from-motion.  Image and Vision Computing, 28(7):1117–1129, 2010.
  * [11] J. Heo and M. Savvides.  In between 3D Active Appearance Models and 3D Morphable Models.  In Computer Vision and Pattern Recognition, 2009.
  * [12] G. B. Huang, M. Ramesh, T. Berg, and E. Learned-Miller.  Labeled faces in the wild: A database for studying face recognition in unconstrained environments.  Technical Report 07-49, University of Massachusetts, Amherst, October 2007\.
  * [13] M. Kazhdan, M. Bolitho, and H. Hoppe.  Poisson surface reconstruction.  In Proceedings Symposium on Geometry Processing (SGP ’06, 32(3):10–10, 2006.
  * [14] I. Kemelmacher, Shlizerman and R. Basri.  3d face reconstruction from a single image using a single reference face shape.  Pattern Analysis and Machine Intelligence, IEEE Transactions on, 33(2):394–405, 2011.
  * [15] I. Kemelmacher, Shlizerman and S. M. Seitz.  Face reconstruction in the wild.  In Computer Vision (ICCV), 2011 IEEE International Conference on, pages 1746–1753. IEEE, 2011.
  * [16] H.-S. Koo and K.-M. Lam.  Recovering the 3d shape and poses of face images based on the similarity transform.  Pattern Recognition Letters, 29(6):712–723, 2008.
  * [17] E. Prados and O. Faugeras.  Shape from shading: a well-posed problem?  In Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on, volume 2, pages 870–877. IEEE, 2005.
  * [18] J. Roth, Y. Tong, and X. Liu.  Unconstrained 3d face reconstruction.  Trans. Graph, 33(4):43, 2014.
  * [19] A. K. Roy-Chowdhury and R. Chellappa.  Statistical bias in 3-d reconstruction from a monocular video.  Image Processing, IEEE Transactions on, 14(8):1057–1062, 2005.
  * [20] Z.-L. Sun, K.-M. Lam, and Q.-W. Gao.  Depth estimation of face images using the nonlinear least-squares model.  Image Processing, IEEE Transactions on, 22(1):17–30, 2013.
  * [21] C. Tomasi and T. Kanade.  Shape and motion from image streams under orthography: a factorization method.  International Journal of Computer Vision, 9(2):137–154, 1992.
  * [22] H. Wang, H. Wei, and Y. Wang.  Face representation under different illumination conditions.  In Multimedia and Expo, 2003. ICME ’03. Proceedings. 2003 International Conference on, pages 285–288, 2003.
  * [23] T. Weise, B. Leibe, and L. Van Gool.  Fast 3d scanning with automatic motion compensation.  In Computer Vision and Pattern Recognition, 2007. CVPR’07. IEEE Conference on, pages 1–8. IEEE, 2007.
  * [24] C. Wu, C. Stoll, L. Valgaerts, and C. Theobalt.  On-set performance capture of multiple actors with a stereo camera.  ACM Transactions on Graphics (TOG), 32(6):161, 2013.
  * [25] C. Yang, J. Chen, N. Su, and G. Su.  Improving 3d face details based on normal map of hetero-source images.  In 2014 IEEE Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), pages 9–14, 2014.
  * [26] A. L. Yuille, D. Snow, R. Epstein, and P. N. Belhumeur.  Determining generative models of objects under varying illumination: Shape and albedo from multiple images using svd and integrability.  International Journal of Computer Vision, 35(3):203–222, 1999.
  * [27] L. Zhang, N. Snavely, B. Curless, and S. M. Seitz.  Spacetime faces: High-resolution capture for~ modeling and animation.  In Data-Driven 3D Facial Animation, pages 248–276. Springer, 2007\.
  * [28] R. Zhang, P.-S. Tsai, J. E. Cryer, and M. Shah.  Shape-from-shading: a survey.  Pattern Analysis and Machine Intelligence, IEEE Transactions on, 21(8):690–706, 1999.