## 2 Methodology
### Algorithm Based on B-Spline  
An alternative model is the free form deformation (FFD) based on B-spline, defined as a mapping from \(\mathbb{R}^{3}\) to \(\mathbb{R}^{3}\) of a tensor product of Bernstein polynomials with three variables [24]. The basic idea of FFDs is to deform an object by manipulating an underlying mesh of control points as seen in Figure 2. The resulting deformation controls the object's shape and produces a smooth and continuous \(C^{2}\) transformation [25; 26]. To define a FFD based on B-spline, we denote the domain of the volume image as \(\Omega=\{(x,y,z)\mid 0\leq x<X,\;0\leq y<Y,\;0\leq z<Z\}\). Let \(\Phi\) denote the mesh \(n_{x}\times n_{y}\times n_{z}\) of control points \(\phi_{i,j,k}\) with uniform spacing \(\delta\). Then the FFD can be written as a 3-D tensor product of the 1-D family of cubic B-splines  
\[T(x,y,z)=\sum_{l=0}^{3}\sum_{m=0}^{3}\sum_{n=0}^{3}B_{l}(u)\cdot B_{m}(v) \cdot B_{n}(w)\cdot\phi_{i+l,j+m,k+n} \tag{4}\]  
where \(i=\left\lfloor\frac{x}{n_{x}}\right\rfloor-1,j=\left\lfloor\frac{y}{n_{y}} \right\rfloor-1,z=\left\lfloor\frac{z}{n_{z}}\right\rfloor-1\), \(u=\frac{x}{n_{x}}-\left\lfloor\frac{x}{n_{x}}\right\rfloor\), \(v=\frac{y}{n_{y}}-\left\lfloor\frac{y}{n_{y}}\right\rfloor\), \(w=\frac{z}{n_{z}}-\left\lfloor\frac{z}{n_{z}}\right\rfloor\) and where \(B_{l}\), \(B_{m}\), and \(B_{n}\) represent the basis functions for the B-spline for parameters \(u\), \(v\) and \(w\)[3], as shown in one of the cases below  
Figure 2: The images (**a**–**d**) are the CT scans with their corresponding meshes of control points (intersection between lines) with the corresponding deformation. \[B_{0}(u) =\frac{(1-u)^{3}}{6}\] \[B_{1}(u) =\frac{3\cdot u^{3}-6\cdot u^{2}+4}{6}\] \[B_{2}(u) =\frac{-3\cdot u^{3}+3\cdot u^{2}+3\cdot u+1}{6}\] \[B_{3}(u) =\frac{u^{3}}{6}\]  
[PERSON] [27] makes an interesting FFD proposal where a set of control points in the fixed image and moving image are established. Their respective deformation describes the displacements of these points and constructs a DVF \(\mathbf{d}\). The goal is to find a measure of the total variation for the regularization term in the objective function [27]. In effect, given the images \(I_{F}\) and \(I_{M}\), the displacement vectors of the control points are described as \(\mathbf{d}=\{\mathbf{d}_{1},\dots,\mathbf{d}_{N}\}\) for a N-dimensional voxel and their gradient is described as \(\
abla\mathbf{d}=\{\
abla\mathbf{d}_{1},\dots,\
abla\mathbf{d}_{N}\}\). The initial model is to find \(\mathbf{d}^{*}\) such that  
\[\mathbf{d}^{*}=\arg\min_{\mathbf{d}}\ E_{D}(\mathbf{d};I_{F},I_{M})+\lambda E _{R}(\mathbf{d}) \tag{5}\]  
where \(E_{D}\) corresponds to the similarity measure and \(E_{R}\) to the regularization term. The central idea is to propose a term that estimates the transformation \(\mathbf{d}^{*}\). The most used regularization term is based on the \(\ell_{2}\) norm and is the smooth regularization model described as [3]:  
\[E_{R}^{\text{sm}}(\mathbf{d})=v\cdot\sum_{n\leq N}\|\
abla\mathbf{d}_{n}\|_{2} ^{2} \tag{6}\]  
where \(v\) is the volume of the voxel. This norm allows any continuous optimization method to minimize (5), given that it is smooth and allows for analytical differentiation with respect to the transformation parameters. On the other hand, we propose how to solve the discontinuity problem for elastic registrations. Starting off with the anisotropic total variation regularization, based on the norm \(\ell_{1}\) given by  
\[E_{R}^{\text{sTV}}(\mathbf{d})=v\cdot\sum_{n\leq N}\|\
abla\mathbf{d}_{n}\|_{1} \tag{7}\]  
we build a total variation regularization term based on (6) and (7) for the \(l\) control points:  
\[T(x,y,z)=\sum_{l=0}^{3}\sum_{m=0}^{3}\sum_{n=0}^{3}B_{l}(u)\cdot B_{m}(v) \cdot B_{n}(w)\cdot\phi_{i+l,j+m,k+n} \tag{4}\]0  
for the DVF given by \(\mathcal{D}(\mathbf{d})=[\
abla_{1}\mathbf{d}_{1}\ \
abla_{2}\mathbf{d}_{1}\ \dots\ \
abla_{1}\mathbf{d}_{2}\ \dots\ \
abla_{N}\mathbf{d}_{N}]\).