\documentclass{article}

\newcommand{\ba}{\left( \begin{array}}
\newcommand{\ea}{\end{array} \right)}
\newcommand{\bq}{\begin{eqnarray*}}
\newcommand{\eq}{\end{eqnarray*}}
\newcommand{\bqn}{\begin{equation}}
\newcommand{\eqn}{\end{equation}}

\newcommand{\ot}{\displaystyle \circ}

\newcommand {\diff}[1] {\frac{\partial}{\partial #1}}
\newcommand{\jacob}[3]{\frac{\partial^2 #3}{\partial #1 \partial #2}}
\newcommand{\der}[2]{\frac{\partial #2}{\partial #1}}
\newcommand{\p}[1]{\partial_{#1}}

%\renewcommand{\baselinestretch}{1.5}

\begin{document}
\newtheorem{lem}{Lemma}
\newtheorem{thm}{Theorem}
\newtheorem{defn}{Definition}

\title{\bf B-Splines}

\author{Moo K. Chung \\
\it{Department of Mathematics and Statistics}\\
\it{McGill University}}
\date{April 17, 2000}
\maketitle
%\large

\section{B-Spline Interpolation of Curves and Surfaces}
A special B-Spline called {\it Bezier curves} was developed by
Bezier and De Casteljau (Gu\'{e}ziec, 1996) in the 1970's for
CAD/CAM operations for car manufacturing industry. It became the
foundation of the entire Adobe Postscript drawing model. We can
construct a parametric curve interpolating given points
$V_0,\dots, V_m$ using Bezier curves. The simplest Bezier curve is
a line segment $V_{0,1}(\lambda)$ joining two points $V_0$ and
$V_1$ and parameterized by $\lambda \in (0,1)$
 such that
 $$V_{0,1}(\lambda) = (1-\lambda)V_0 + \lambda V_1.$$
For three points $V_0,V_1$ and $V_2$, we have a quadratic Bezier
curve $$ V_{0,2}(\lambda) = (1-\lambda) V_{0,1} + \lambda V_{1,2}=
(1-\lambda)^2V_0 + 2(1-\lambda)\lambda V_1 + \lambda^2 V_2.$$ In
general, for $m+1$ points $V_0,\dots,V_m$, $m$-th order Bezier
curve is
$$V_{0,m}(\lambda) = \sum_{j=0}^m
{m \choose j} (1-\lambda)^{m-j} \lambda^jV_j$$ Note that the basis functions
$\{ {m \choose j} (1-\lambda)^{m-j} \lambda^j,j=0,\dots,m \}$ are
{\it Bernstein polynomials} of degree $m$. Cubic Bezier curves are the
most often used Bezier curve because quadratic are necessarily
planner whereas cubics are not and cubics satisfy the minimum
curvature property(Gu\'{e}ziec, 1996) or minimum strain energy(McLeod, 1998)
which make them a more suitable
tool for data approximation. To avoid increasing the degree of
the Bezier curve, we need to piece together Bezier curves. If
continuity conditions are satisfied for each Bezier curve
segments, the result is a B-spline curve. In general, a B-spline
curve of degree $K-1$ with $m+1$ vertices $V_0,\dots V_m$ is
defined as $$S(u) = \sum_{j=0}^m B_{j,K}(u)V_j,$$ where the
B-spline functions $\{B_{j,k}(u),j=0,\dots, m \}$ are defined
recursively with $m+K+1$ knot values $u_j$ such that
\bq B_{j,1}(u) &=& 1, u_j \leq u < u_{j+1} \\
               & & 0, \mbox{ otherwise}
\eq
$$B_{j,k} = \frac{u-u_j}{u_{j+k-1} - u_j}B_{j,k-1}(u) +
\frac{u_{j+k}-u}{u_{j+k} - u_{j+1}}B_{j+1,k-1}(u).$$
For knot values $\{0,0,0,1,1,1\}$, we have the quadratic Bezier curve and
for knot values $\{0,0,0,0,1,1,1,1\}$, we have the cubic Bezier curve. \\

Then we can use these B-spline curves to generate surfaces.
Consider a rectangular mesh of vertices $\{V_{ij}\}$. A B-spline
patch $X$ of degree $K-1$ can be defined by the tensor product:
$$X(u,v) = \sum_{i,j}B_{i,K}(u)B_{j,K}(v)V_{ij}.$$
The advantage for using spline to represent the surface is that it is
easy to evaluate the curvature of a surface or other geometric
characteristics of the surface because polynomial functions can
be differentiated easily.

\section{Polynomial Regression of Surfaces}
Instead of using B-splines to form a parametric surface, there
is another method called {\it polynomial regression} (Ripple, 1981).
This is a smoothing technique to fit the given vertices $V_0,\dots,V_m$
by least squares method to a polynomial function of the form
$$f(x,y) = \sum_{i+j \leq p} a_{ij}x^iy^j$$

Then for $V_k=(x_k,y_k,z_k) \in \Re^3$, $P=\frac{(p+1)(p+2)}{2}$
unknown coefficients $a_{ij}$'s are chosen to minimize the residual
$$\sum_{k=0}^m [z_k - f(x_k,y_k)]^2.$$ The drawback of the polynomial
regression is that there is tendency to wave the outer most vertices
to find vertices in the center. So if we are warping the whole cortical
surface using a single polynomial function $f(x,y)$, the fit may break
down near the extreme points in $D$. However, from the purely statistical
point of view, we do not need to warp the cortex by a single
function $f$ over all domain $D$. What we really need is the parametric
representation in the neighborhood of each voxel. Later we will show that
the quadratic polynomial regression at each voxel with its
neighboring voxels is sufficient for analyzing the variability in cortex using
the Riemannian metric tensor approach. \\

Let $$f_1(X)=1,f_2(X)=x,f_3(X)=y,f_4(X) = x^2,$$
$$f_5(X)=xy,f_6(X)=y^2,\dots,f_P(X)=y^P.$$ and the coefficients
$a_{ij}$'s as $\beta_0,\beta_1,\beta_2,\dots, \beta_P$. Then the
problem becomes a standard {\it multiple regression analysis}:
$$z=X\beta + \epsilon,$$ where the matrix $X=(X_{ij})
=(f_j(x_i,y_i))$, the
parameter $\beta = (\beta_1,\dots,\beta_P)$ and $z=(z_1,\dots,z_k)$.
Then it could be given to any statistical package to estimate the
parameter $\beta$ automatically. Orthogonal polynomials are
often used for polynomial regression, especially with regularly spaced
data but modern multiple regression techniques rely on orthogonalizing
the functions $f_1,\dots,fp$ and make the orthogonal polynomial
technique unnecessary. Other families of surfaces has been suggested
but are rarely used. One of them is to use finite Fourier series: \bq
f(x,y) &=& \sum_{0 \leq i,j \leq r}
[a_{ij}\sin(i\omega_1x)\sin(j\omega_2y)
+ b_{ij}\sin(i\omega_1x)\cos(j\omega_2y) \\
&& + c_{ij}\cos(i\omega_ix)\sin(j\omega_2y)
+ d_{ij}\cos(i\omega_1x)\cos(j\omega_2y)]
\eq
Such a surface requires a large number of coefficients,
i.e. $4r^2 +4r+1$, so even with $r=3$ there are 49 coefficients
plus the fundamental frequencies $\omega_1$ and $\omega_2$ to be
estimated. Hence, surface warping based on Fourier series on the large
data set is computationally intensive although the result might be better than
the polynomial regression. With small values of $r$, the surface shows ripples at
the fundamental frequencies in the $x$ and $y$ directions and so is
really suitable for periodic surface. In certain sense, the cortical
surface is somewhat periodic in a sense that there is a continuous ups
and downs of the cortical folds. So the above Fourier series expansion can be
used to represent the cortex.
Instead of using the flat domain $D \in \Re^2$, we can also
find a smooth mapping from the cortical surface to a unit sphere
using the spherical harmonics,
which allows points on the sphere to be represented by two angles
$\theta$ and $\psi$ (Bakircioglu {it et al.}, Thompson {\it et al.}, 1996)


\end{document}
\begin{thebibliography} \\
\bibitem{Bakircioglu} MM Bakircioglu, S Joshi, MI Miller, {\sl Landmark Matching on Brain Surfaces via Large Deformation Diffeomorphisms on the Sphere}
\bibitem{Dale} AM Dale, B Fischl, {\sl Cortical Surface-Based Analysis I: Segmentation and Surface Reconstruction.}
\bibitem{Fischl} B Fischl, MI Sereno, {\sl Cortical Surface-Based Analysis II: Inflation, Flattening, and a Surface-Based Coordinate System.}
\bibitem{Davatzikos} C Davatzikos RN Bryan, {\sl Using a Deformable Surface Model to Obtain a Shape Representation of the Cortex}, IEEE Trans. Medical Imaging, Vol. 15, No. 6, 1996
\bibitem{Davatzikos2} C Davatzikos, {\sl Nonlinear Registration of Brain Images Using Deformable Models}, IEEE, 1996
\bibitem{Davatzikos3} C Davatzikos, M Vaillant, SM Resnick, JL Prince, S Letovsk, RN Bryan,{\sl A Computerized Approach for Morphological Analysis o the Corpus Callosum}
\bibitem{Gueziec} A Gu\'{e}ziec, {\sl Spline Curves and Surfaces for Data Modeling}, Advances in Morphometrics, Plenum Press, New Work, 1996
\bibitem{Joshi} SC Joshi, J Wang, MI Miller, {\sl On the differential geometry of the cortical surface}
\bibitem{Khaneja} N Khaneja, MI Miller, U. Grenander {\sl Dynamic Programming Generation of Curves on Brain Surfaces}
\bibitem{Ripley} BD Ripley, {\sl Spatial Statistics}, John Wiley $\&$ Sons, Inc, 1981
\bibitem{Kreyszig} E Kreyszig, {\sl Differential Geometry}. University of Toronto Press, Toronto, 1964
\bibitem{MacDonald} JD MacDonald, {\sl A Method for Identifying Geometrically Simple Surfaces from Three Dimensional Images}, Ph.D. Thesis, McGill University, 1997
\bibitem{McLeod} RJY McLeod, ML Baart, {\sl Geometry and Interpolation of Curves and Surfaces}, Cambridge University Press, 1998
\bibitem{Thompson} P Thompson, AW Toga, {\sl A Surface-Based Technique for Warping Three-Dimensional Images of the Brain}, 1996

\end{thebibliography}