section{Point-based Elastic Shape Registration}label{sec:Point-based-Elastic}Our objective is to find a function that gives the pointcorrespondences between the two given domains (source and target).Let us define the $2D$ shape registration as follows:A map $mathbf{C}^s(p):0,1in R

ightarrow R^2$ defines a planarsource curve with a parameter $p$. The target is defined by$mathbf{C}^t(p):0,1in R

ightarrow R^2$. The cartesiancoordinates of the point vector can be defined by$mathbf{C}(p)=x~y^T$ where $0leq xleq X$ and $0leq yleq Y$.subsection{Finding Correspondences}Assume that $mathbf{C}^t(p)$ is the corresponding point of$mathbf{C}^s(p)$ (the criteria for finding correspondences willcome in the following sections). The output will be a $C^0$ function$f:R^2

ightarrow R^2$ with$f(mathbf{C}^s(p))=mathbf{C}^t(p),~forall pin0,1$. Differentinterpolation functions have been proposed to handle this problemcite{Farin_2004}. We choose the free form deformation $FFD$ model,based on B-splines cite{Parry_FFD_1986,Hawkes_FFD_1999}, which is apowerful tool for modeling deformable objects and has beenpreviously applied to the tracking and motion analysis problems. Thebasic idea is to deform the shape by manipulating a mesh of controlpoints. The resulting deformation controls the shape of the objectand produces a smooth and continuous transformation.Consider a lattice of control points$mathbf{P}=mathbf{P}_{m,n};~m=1,dots,M.~n=1,dots,N$, each pointon the source shape will have the following form of deformation:egin{equation}label{e:Deformation}mathbf{L}(p)=sum_{k=0}^3sum_{l=0}^3B_k(u)B_l(v)deltamathbf{P}_{i+k,j+l}end{equation}where $deltamathbf{P}=deltamathbf{P}_{m,n}$ is the control pointdeformation, $i=lfloor(x.(M-1)/X)

floor+1$,$j=lfloor(y.(N-1)/Y)

floor+1$, $u=x.M/X-lfloor(x.M/X)

floor$,$v=y.N/Y-lfloor(y.N/Y)

floor$, and the spline basis functions($B$) are defined in cite{Hawkes_FFD_1999}. So the cubic B-splineis used as an approximation function for our interpolation problem.We propose the following energy to measure the difference betweenthe deformed contour and its target points based on the Euclideandistance:egin{equation}label{e:Energy}E(deltamathbf{P})=int_0^1|mathbf{C}^s(p)+mathbf{L}(p)-mathbf{C}^t(p)|^2dpend{equation}Also, we need to avoid any undesired distortion in the shape andpreserve the regularity of the registration grid flow. Anotherweighted term (by $lambdain R^+$) is added for smoothnessconstrain as follows cite{Terzoplos_1996}:egin{equation}label{e:EnergySmoothed}egin{array}{l}E(deltamathbf{P})=int_0^1|mathbf{C}^s(p)+mathbf{L}(p)-mathbf{C}^t(p)|^2dp+lambdaint_0^1(|mathbf{L}_p|^2+|mathbf{L}_{pp}|^2)dpend{array}end{equation}where $mathbf{L}_p$ and $mathbf{L}_{pp}$ are the first and secondderivatives respectively of the deformation vector with respect tothe parameterizations $p$. The above objective function is requiredto be minimal to calculate the control points locations and hencethe deformation field at each point in the domain. Gradient descentand calculus of variation are used to optimize the given function asfollows:egin{equation}label{e:EnergySmoothedOptimize}egin{array}{l}\frac{partial E}{partialdeltamathbf{P}}=2int_0^1(mathbf{L}(p)-mathbf{C}(p))^Tfrac{partialmathbf{L}}{partialdeltamathbf{P}}dp+2lambdaint_0^1(mathbf{L}_p^Tfrac{partialmathbf{L}_p}{partialdeltamathbf{P}}+mathbf{L}_{pp}^Tfrac{partialmathbf{L}_{pp}}{partialdeltamathbf{P}})dpend{array}end{equation}where $mathbf{C}(p)=mathbf{C}^t(p)-mathbf{C}^s(p)$. A detailedillustration for the local deformation derivatives are found incite{Paragios_PAMI_2006}.By setting the above equation to zero, we get a linear system of thecontrol points coordinates:egin{equation}label{e:Linear System}egin{array}{l}\int_0^1(mathbf{C}(p))^Tfrac{partialmathbf{L}}{partialdeltamathbf{P}}dp=int_0^1mathbf{L}^Tfrac{partialmathbf{L}}{partialdeltamathbf{P}}dp+lambdaint_0^1(mathbf{L}_p^Tfrac{partialmathbf{L}_p}{partialdeltamathbf{P}}+mathbf{L}_{pp}^Tfrac{partialmathbf{L}_{pp}}{partialdeltamathbf{P}})dpend{array}end{equation}The left hand side is free from the coordinates of the controlpoints while the other side is linearly a function of theseunknowns. The following linear equation holds:egin{equation}label{e:Linear Eq}PsiTheta=Lambdaend{equation}where$Theta=deltamathbf{P}_{1,1}^x,dots,deltamathbf{P}_{M,N}^x,deltamathbf{P}_{1,1}^y,dots,deltamathbf{P}_{M,N}^y^T$and $x,~y$ are the coordinates of the embedding space. The othermatrices elements are calculated as follows:egin{equation}label{e:Psi}egin{array}{l}\Psi_{r,c}=int_0^1(mathbf{L}^{r,c})^Tfrac{partialmathbf{L}}{partial heta_r}dp+lambdaint_0^1((mathbf{L}_p^{r,c})^Tfrac{partialmathbf{L}_p}{partial heta_r}+(mathbf{L}_{pp}^{r,c})^Tfrac{partialmathbf{L}_{pp}}{partial heta_r})dpend{array}end{equation}egin{equation}label{e:Lambda}egin{array}{l}\Lambda_r=int_0^1(mathbf{C}(p))^Tfrac{partialmathbf{L}}{partial heta_r}dpend{array}end{equation}where $ heta_rinTheta$ and $mathbf{L}^{r,c}$ stands for thecoefficient of the control point vector calculated from the cubicspline interpolation. Row and column are represented by $r~and~c$respectively.The resulting matrix equation size depends on the number of controlpoints and the space dimensions. Size of the data points of thecurve does not have any impact on the matrices sizes whichguarantees its computational efficiency (this is a huge advantageover the work proposed in cite{Farin_2004} where number of linearequations is equal to number of data points). Also note that thesmoothing constrains do not add extra load to the solution of thematrix equation.subsection{A Coarse to Fine Strategy (Incremental Free Form Deformation)}label{sec:Coarse_Fine}The above objective function is required to be minimized by movingthe control points to get the correct correspondence over shapeboundaries. A very small error can be achieved when using a highresolution control lattice because the number of degrees of freedomincreases. However this is not enough. Such sudden movement willresult in unnecessary cross overs of the domain grid lines and theregistration process will be meaningless. This will result inchanging and corrupting the object topology. A better way is to movethe grid step by step towards the target.To illustrate the problem, we formed a contour which representssides of a square (this will stand for the source). Another contouris obtained by rotating the first one by an angle of $pi/2$ (thiswill stand for the target). Correspondences between the source andtarget shapes are established by the rotation effect. A controllattice is constructed and computed using the above formula for thepoint-based registration. When using a relatively high resolution of$12 imes 12$ lattice, fold-overs occur as shown in Fig.

ef{fig:Good_vs_Bad} (top left image). To avoid this, a coarse tofine strategy is used (equivalent to the incremental free formdeformations used in cite{Paragios_PAMI_2006}). We start with aresolution of $4 imes 4$ and solve for the deformation. Iterativelywe increase the resolution as $5 imes 5$, $6 imes 6$, and so onand so forth. In each step, the control points positions arecomputed and the contour moves to the new position until asatisfactory error distance is obtained. The result is smooth andthe correspondence is achieved accurately. The final griddeformation is shown in the right image. This process handlesextremely well the error and gives impressive infinitesimal energyfunction and smooth grid deformations at the same time.