Efficient combined SSIM- and landmark-driven image registration in a variational framework

<!-- metadata commented in wiki content
==Efficient combined SSIM- and landmark-driven image registration in a variational framework==

'''Jorge Larrey-Ruiz<sup>a</sup>, Juan Morales-Sánchez<sup>a</sup>, Laura Larrey-Ruiz<sup>b</sup>'''

{|
|-
|<sup>a</sup> Tecnologías de la Información y las Comunicaciones, Universidad Politécnica de Cartagena, 30202 Cartagena, Spain
|}

{|
|-
|<sup>b</sup> Servicio de Patología Digestiva, Consorcio Hospital General Universitario de Valencia, 46014 Valencia, Spain
|}
-->

==Abstract==

Image registration methods are used to establish geometrical correspondences between different datasets. Various characteristics of image data can be exploited to drive image registration algorithms. Thus, the currently available schemes can be roughly divided into two classes: intensity-based and feature-based registration schemes. In this paper, we present a mathematical framework, based on the calculus of variations, for combining these two classes in order to benefit from the advantages of both strategies. The goal is to obtain a registration algorithm which achieves a good matching of datasets near landmark locations but also away from them (by matching the corresponding intensities). The proposed approach includes the novel formulation of a disparity term which simultaneously takes into account the structural similarity index (a similarity measure which considers spatial dependencies in the images) and the location of outstanding points. Since the iteration which results of the variational formulation is translated into the frequency domain, the implementation of the proposed algorithm offers a good speed-performance trade-off when compared to other state-of-the-art image registration implementations. Experimental results show the advantages, in the medical setting, of the combined SSIM- and landmark-based approach over well-established registration techniques which use either landmark or intensity information alone. In particular, the registration of triple-phase 3D computed tomographies of the liver under injection of a contrast agent has been chosen for such a comparison. The datasets are acquired at different times depending on the arrival time of the contrast agent in arteries, portal and hepatic veins, so they have to be registered in order to show the liver structures acquired at each phase in a common framework. These multi-phase studies provide tumor enhancement on the arterial and portal venous phases that support differential diagnosis of lesions in the liver.

'''keywords'''

Image registration, variational methods, Fourier domain, structural similarity index, landmark-driven registration, contrast-enhanced liver CT

==1 Introduction==

Image registration is the process of finding the spatial correspondence between two different datasets (images or volumes), which typically represent different views of the same object or similar ones; in other words, register is equivalent to ''align'' in this context. Particularly, in medical imaging there are several applications that require a registration step (e.g. image fusion, atlas matching, pathological diagnosis). Classical reviews of image registration can be found, for example, in the references <span id='citeF-6'></span>[[#cite-6|[6]]], <span id='citeF-28'></span>[[#cite-28|[28]]], <span id='citeF-47'></span>[[#cite-47|[47]]] or <span id='citeF-29'></span>[[#cite-29|[29]]]. More recent comprehensive overviews about image registration methods can be found in <span id='citeF-37'></span>[[#cite-37|[37]]] and <span id='citeF-36'></span>[[#cite-36|[36]]], where the state-of-the-art advances in this field are described in a systematic fashion, and the techniques applied to medical datasets are emphasized as well. When classified according to the image characteristics used to estimate the transformation which geometrically relates the involved datasets, registration algorithms are usually divided into three groups: geometric feature-based, which rely on the segmentation, done generally before the registration process itself, of part or all of the images, therefore obtaining objects that are registered by minimizing some geometrical distance between them; intensity-based, which maximize a intensity similarity measure computed between points lying at the same spatial position; and iconic feature-based, which can be understood as intermediate between the two previous categories, since they use explicitly some type of geometrical distance in addition to the intensity similarity measure, see <span id='citeF-8'></span>[[#cite-8|[8]]]. The performance of intensity-based image registration methods is highly influenced by the choice of the similarity measure. This measure can be defined directly on image intensities as, e.g., sum of squared differences (SSD), correlation coefficient (CC), correlation ratio (CR), or mutual information (MI). A complete analysis of all these metrics can be found in  <span id='citeF-27'></span>[[#cite-27|[27]]]. These measures often rely on the assumption of independence and stationarity of the intensities from pixel to pixel without considering their spatial dependencies. Alternative metrics have been proposed in order to consider intensity nonstationarities and complex spatially-varying intensity distortions. For instance, in <span id='citeF-31'></span>[[#cite-31|[31]]] a similarity measure based on the residual complexity is described. Another example is the structural similarity index (SSIM), introduced by <span id='citeF-43'></span>[[#cite-43|[43]]], which achieves a good trade-off between speed and performance and is able to cope with complex relationships between image intensity values. Due to the mathematical properties of SSIM <span id='citeF-7'></span>[[#cite-7|[7]]] and its potentials in both theoretical development and practical applications, it is being incorporated into optimization frameworks in order to improve perceived image/video quality in a number of image processing problems, as e.g. <span id='citeF-17'></span>[[#cite-17|[17]]] or <span id='citeF-34'></span>[[#cite-34|[34]]]. This similarity measure has been previously used in image registration in <span id='citeF-4'></span>[[#cite-4|[4]]], <span id='citeF-1'></span>[[#cite-1|[1]]], <span id='citeF-41'></span>[[#cite-41|[41]]], <span id='citeF-42'></span>[[#cite-42|[42]]], <span id='citeF-45'></span>[[#cite-45|[45]]], or <span id='citeF-19'></span>[[#cite-19|[19]]], but its use is still not widespread |for example, the evaluation addressed in <span id='citeF-33'></span>[[#cite-33|[33]]] does not include SSIM and the review in <span id='citeF-36'></span>[[#cite-36|[36]]] does not deal with similarity measures based on the SSIM|.

Intensity-based approaches are in general fully automatic and usually yield good registration results. However, they may perform poorly for specific, important locations such as anatomical landmarks. On the other hand, geometric feature-based registration techniques are designed to accurately match user specified landmarks, see e.g. <span id='citeF-35'></span>[[#cite-35|[35]]]. Recent examples of landmark-based registration methods in medical imaging scenarios can be found in <span id='citeF-44'></span>[[#cite-44|[44]]] or <span id='citeF-21'></span>[[#cite-21|[21]]]. A common drawback of most landmark-driven registration methods is the fact that the intensities of the datasets are completely neglected. Consequently, the registration outcome away from the landmarks may be very poor. In the literature, the registration results are typically obtained by interpolating the transformation with a ''thin-plate spline'' (TPS) model <span id='citeF-5'></span>[[#cite-5|[5]]] or by combining iterative closest point (ICP) registration and parametric relaxation <span id='citeF-38'></span>[[#cite-38|[38]]], among other techniques. Although these approaches produce a smooth transformation from one dataset to another, they do not define a consistent correspondence between the two datasets except at the landmark points. Some previous attempts to combine the sparse correspondences with a intensity-based registration method can be found in <span id='citeF-18'></span>[[#cite-18|[18]]] or <span id='citeF-39'></span>[[#cite-39|[39]]]. However, to the authors knowledge, there has been no attempt to formulate a registration method which combines in the distance term both a SSIM-based term and a landmark-based term. In recent works, such as <span id='citeF-32'></span>[[#cite-32|[32]]] or <span id='citeF-26'></span>[[#cite-26|[26]]], the authors of this paper dealt with intensity-based registration approaches within the variational framework first presented in <span id='citeF-25'></span>[[#cite-25|[25]]], whose formulation in the frequency domain allowed for implementations of high efficiency <span id='citeF-40'></span>[[#cite-40|[40]]]. Please refer to <span id='citeF-2'></span>[[#cite-2|[2]]] or <span id='citeF-46'></span>[[#cite-46|[46]]] for further details regarding variational image registration. The aforementioned approaches involved similarity terms derived from SSD- or CR-based measures, so none of them included a SSIM-based nor a landmark-based distance term. Regardless, the method presented in <span id='citeF-26'></span>[[#cite-26|[26]]] achieved better results than publicly available state-of-the art approaches such as ''Elastix'' <span id='citeF-20'></span>[[#cite-20|[20]]] and ANTs (<math display="inline">http://stnava.github.io/ANTs/</math>).

The paper is organized as follows: Section [[#2 Mathematical framework|2]] addresses the formulation, within the variational framework, of a novel registration method which combines in the distance term both a SSIM-based term and a landmark-based term. Section [[#3 Results|3]] shows the results of applying the proposed registration algorithm to a medical imaging scenario; in particular, the image registration of triple-phase 3D computed tomography datasets of the liver under injection of a contrast agent has been chosen; the suitability of our approach in such a scenario is also proved in this section through the comparison of the results it provides and the outcome obtained with the CR-based version of the algorithm, which was recently reported to outperform other well-known approaches in the medical setting <span id='citeF-26'></span>[[#cite-26|[26]]]. The conclusions of the paper as well as the future lines of research are discussed in Section [[#4 Conclusions and discussion|4]]. Finally, an appendix closes the paper with the detailed derivation of the external forces field related to the proposed SSIM-based distance term (i.e., the first addend of the corresponding Euler-Lagrange equation).

==2 Mathematical framework==

Let <math display="inline">R</math> and <math display="inline">T</math> be two datasets, a reference and a template respectively, which represent the same object (or similar ones) by using the same or different imaging modalities. These <math display="inline">d</math>-dimensional datasets are defined as <math display="inline">R,T:\mathbf{\Psi }\rightarrow \mathcal{I}</math>, with <math display="inline">\mathbf{\Psi }\subset \mathbb{R}^d</math> being the domain where they are supported, and <math display="inline">\mathcal{I}\subset \mathbb{R}</math> representing the intensity level in each spatial coordinate <math display="inline">\mathbf{x}=(x_1,\ldots ,x_d)</math>. We search for a displacement field <math display="inline">\mathbf{u}:\mathbf{\Psi }\rightarrow \mathbb{R}^d</math> that makes the transformed template <math display="inline">T_\mathbf{u}:\mathbf{\Psi }\rightarrow \mathcal{I}</math> similar to the reference dataset in the geometrical sense, i.e., <math display="inline">\mathbf{x}_{T_\mathbf{u}}=\mathbf{x}_T+\mathbf{u}(\mathbf{x}_T)\approx \mathbf{x}_R</math>, where <math display="inline">\mathbf{u}(\mathbf{x})=\left(u_1(\mathbf{x}),\ldots ,u_d(\mathbf{x})\right)^\top </math>. This problem can be approached in terms of the variational calculus by defining a joint energy functional to be minimized:

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{J}[\mathbf{u}]=\mathcal{D}[R,T;\mathbf{u}]+\alpha \,\mathcal{S}[\mathbf{u}]\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

where <math display="inline">\mathcal{D}</math> is an energy term which quantifies the disparity between the datasets, <math display="inline">\mathcal{S}</math> is a penalty term which measures the roughness of the solution <math display="inline">\mathbf{u}(\mathbf{x})</math>, and <math display="inline">\alpha{>0}</math> is a scalar parameter, usually referred to as the regularization parameter, which is used to control and weight the influence of the penalty term versus the distance term. It should be noted that the present formulation considers the different functions to be continuous in the domain <math display="inline">\mathbf{\Psi }</math> where they are supported; the natural spatial (and temporal) discretization is tackled later in this section.

In the literature, several choices for the dissimilarity measure <math display="inline">\mathcal{D}</math> can be found. Depending on the particular datasets to be compared, statistical measures derived from the MI or the CR are usually the most appropriate candidates. As an alternative, in this work we propose the novel approach of incorporating both the SSIM <span id='citeF-43'></span>[[#cite-43|[43]]] and the Euclidean distance between identifiable corresponding points (i.e., landmarks) into the variational framework which is being presented. The resulting disparity term is the following:

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{D}[R,T;\mathbf{u}] =\frac{-(2\,\mu _R\,\mu _T+C_1)(2\,\sigma _{RT}+C_2)}{(\mu _R^2+\mu _T^2+C_1)(\sigma _{RR}+\sigma _{TT}+C_2)}+\frac{\beta }{2}\int _{\mathbf{\Psi }}\sum _{j=1}^N\left\|\mathbf{x}_{R_j}-\mathbf{x}_{T_j}-\mathbf{u}(\mathbf{x}_{T_j})\right\|^2d\mathbf{x}\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

where <math display="inline">\mu _R</math>, <math display="inline">\mu _T</math>, <math display="inline">\sigma _{RR}</math> and <math display="inline">\sigma _{TT}</math> denote the mathematical expectations and variances of the intensity levels of <math display="inline">R</math> and <math display="inline">T_\mathbf{u}</math>, respectively; <math display="inline">\sigma _{RT}</math> represents the covariance; <math display="inline">C_1</math> and <math display="inline">C_2</math> are small constants which aim to characterize the saturation effects of the visual system at low luminance and contrast regions and which assure numerical stability when the expectations or variances are very close to zero; <math display="inline">\beta{>0}</math> is a scalar used to weight the influence of the landmark-based term versus the SSIM-based term; <math display="inline">N</math> is the number of identified landmarks; and <math display="inline">\mathbf{x}_{R_j}</math> and <math display="inline">\mathbf{x}_{T_j}</math> denote the spatial position of the <math display="inline">j</math>-th landmark in <math display="inline">R</math> and <math display="inline">T</math>, respectively.

The regularization term <math display="inline">\mathcal{S}</math> is used to add some prior knowledge on the displacement field, thus preferentially obtaining more likely solutions and giving the smoothness characteristics to <math display="inline">\mathbf{u}(\mathbf{x})</math>. In this work, the regularizer is defined by

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{S}[\mathbf{u}]=\frac{1}{2}\sum _{l=1}^d\int _\mathbf{\Psi }\left\|\nabla ^\lambda u_l\right\|^2\,d\mathbf{x}\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

where <math display="inline">\nabla </math> denotes the <math display="inline">d</math>-dimensional gradient operator, and <math display="inline">\lambda \in \{{}1,2\}{}</math> is the regularization order: if <math display="inline">\lambda=1</math>, the resulting term is the diffusion smoother <span id='citeF-14'></span>[[#cite-14|[14]]]; if <math display="inline">\lambda=2</math>, the regularizer is the curvature smoother <span id='citeF-15'></span>[[#cite-15|[15]]].

According to the variational calculus, a necessary condition for a minimizer <math display="inline">\mathbf{u}</math> of the joint energy functional ([[#eq-1|1]]) is that the first variation of <math display="inline">\mathcal{J}[\mathbf{u}]</math> in any direction (also known as the ''Gateaux'' derivative) vanishes for all possible perturbations <math display="inline">\mathbf{z}\in \mathbb{R}^d</math>:

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathrm{d}\mathcal{J}[\mathbf{u};\mathbf{z}]=\frac{\partial }{\partial \epsilon }\Bigl.\mathcal{J}[\mathbf{u}+\epsilon \mathbf{z}]\Bigr|_{\epsilon=0}=0\,,\quad \forall \mathbf{z}\in \mathbb{R}^d\,. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

The computation of the previous ''Gateaux'' derivative leads to the following Euler-Lagrange (E-L) equation:

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{1}{\mathcal{V}_\mathbf{\Psi }}\left(\hat{G}_\rho \ast \frac{\partial L_\mathbf{u}^{\mbox{SSIM}}}{\partial i_t}\right)\nabla T_{\mathbf{u}}+\beta \sum _{j=1}^N\left(\mathbf{x}_{T_{\mathbf{u}j}}-\mathbf{x}_{R_j}\right)+\alpha \,(-\Delta )^\lambda \mathbf{u}=\mathbf{0}\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

where <math display="inline">\mathcal{V}_\mathbf{\Psi }</math> represents the hypervolume of <math display="inline">\mathbf{\Psi }</math>; <math display="inline">\hat{G}_\rho :\mathcal{I}^2\rightarrow \mathbb{R}</math> is a Gaussian kernel (strictly positive and differentiable); <math display="inline">\ast </math> denotes the <math display="inline">d</math>-dimensional convolution operator; <math display="inline">\Delta </math> is the Laplacian operator; finally, we find the term

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{\partial L_\mathbf{u}^{\mbox{SSIM}}}{\partial i_t}=\frac{\gamma _0+C_1\,\gamma _1+C_2\,\gamma _2+C_1C_2\,\gamma _3}{(\mu _R^2+\mu _T^2+C_1)^2(\sigma _{RR}+\sigma _{TT}+C_2)^2}\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

where the variables <math display="inline">\gamma _0</math>, <math display="inline">\gamma _1</math>, <math display="inline">\gamma _2</math> and <math display="inline">\gamma _3</math> are defined as

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\gamma _0=4\,\mu _R\Big(\mu _R^2(\sigma _{RR}\,\sigma _{RT}+\sigma _{TT}\,\sigma _{RT}-\mu _R\,\mu _T\,\sigma _{RR}-\mu _R\,\mu _T\,\sigma _{TT}+2\,\mu _T^2\,\sigma _{RT}+\mu _T\,\sigma _{RR}\,i_r+\mu _T\,\sigma _{TT}\,i_r-2\,\mu _T\,\sigma _{RT}\,i_t)</math>
|-
| style="text-align: center;" | <math> -\mu _T^2(\sigma _{RR}\,\sigma _{RT}+\sigma _{TT}\,\sigma _{RT}+\mu _R\,\mu _T\,\sigma _{RR}+\mu _R\,\mu _T\,\sigma _{TT}-2\,\mu _T^2\,\sigma _{RT}-\mu _T\,\sigma _{RR}\,i_r-\mu _T\,\sigma _{TT}\,i_r+2\,\mu _T\,\sigma _{RT}\,i_t)\Big)\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\gamma _1=2\Big(\mu _R(-\mu _R^2\,\sigma _{RR}-\mu _R^2\,\sigma _{TT}+2\,\sigma _{RR}\,\sigma _{RT}+2\,\sigma _{TT}\,\sigma _{RT}+\mu _R\,\sigma _{RR}\,i_r+\mu _R\,\sigma _{TT}\,i_r-\sigma _{RR}\,C_1-\sigma _{TT}\,C_1)</math>
|-
| style="text-align: center;" | <math> -\mu _T(2\,\mu _R^2\,\sigma _{RR}+2\,\mu _R^2\,\sigma _{TT}+\mu _R\,\mu _T\,\sigma _{RR}+\mu _R\,\mu _T\,\sigma _{TT}-2\,\mu _R\,\sigma _{RR}\,i_r-2\,\mu _R\,\sigma _{TT}\,i_r-\mu _T\,\sigma _{RR}\,i_r-\mu _T\,\sigma _{TT}\,i_r)</math>
|-
| style="text-align: center;" | <math> +2\,\mu _T\,\sigma _{RT}(\mu _R^2+\mu _T^2-\sigma _{RR}-\sigma _{TT}+2\,\mu _R\,\mu _T-2\,\mu _R\,i_t-\mu _T\,i_t+C_1)+C_1\,i_r(\sigma _{RR}+\sigma _{TT})-2\,\sigma _{RT}\,i_t(\mu _R^2+C_1)\Big)\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\gamma _2=2\,\mu _R\Big(\mu _R^2(\sigma _{RR}+\sigma _{TT}+2\,\sigma _{RT}-2\,\mu _R\,\mu _T+2\,\mu _T^2+2\,\mu _T\,i_r-2\,\mu _T\,i_t+C_2)</math>
|-
| style="text-align: center;" | <math> -\mu _T^2(\sigma _{RR}+\sigma _{TT}+2\,\sigma _{RT}+2\,\mu _R\,\mu _T-2\,\mu _T^2-2\,\mu _T\,i_r+2\,\mu _T\,i_t+C_2)\Big)\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\gamma _3=2\Big(\mu _R(\sigma _{RR}+\sigma _{TT}-2\,\sigma _{RT}-\mu _R^2+\mu _T^2+\mu _R\,i_r-\mu _R\,i_t+2\,\mu _T\,i_r+C_2)</math>
|-
| style="text-align: center;" | <math> -\mu _T(\sigma _{RR}+\sigma _{TT}+2\,\sigma _{RT}+\mu _R^2-\mu _T^2+2\,\mu _R\,i_t-\mu _T\,i_r+\mu _T\,i_t-C_1+C_2)-C_1\,i_t\Big)\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

with <math display="inline">i_r\,,i_t\in \mathcal{I}\subset \mathbb{R}</math> being the intensities of <math display="inline">R</math> and <math display="inline">T_\mathbf{u}</math>, respectively.

In [[#5 Variational formulation of SSIM|5]], the first addend in equation ([[#eq-5|5]]) is deduced. The proof of the internal forces field of the Euler-Lagrange equation (i.e., the last addend in eq.([[#eq-5|5]])) can be found in <span id='citeF-24'></span>[[#cite-24|[24]]]. Finally, the addend which corresponds to the landmark-based term can be easily obtained by applying the definition of the ''Gateaux'' derivative to the last addend in eq.([[#eq-2|2]]):

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathrm{d}\mathcal{D}^{\mbox{LM}}[R,T;\mathbf{u};\mathbf{z}]=\frac{\partial }{\partial \epsilon }\Bigl.\frac{\beta }{2}\int _{\mathbf{\Psi }}\sum _{j=1}^N\left\|\mathbf{x}_{R_j}-\mathbf{x}_{T_j}-\mathbf{u}(\mathbf{x}_{T_j})-\epsilon \mathbf{z}(\mathbf{x}_{T_j})\right\|^2d\mathbf{x}\Bigr|_{\epsilon=0}</math>
|-
| style="text-align: center;" | <math> =\Bigl.\int _{\mathbf{\Psi }}\Bigl\langle\beta \sum _{j=1}^N\mathbf{x}_{R_j}-\mathbf{x}_{T_j}-\mathbf{u}(\mathbf{x}_{T_j})-\epsilon \mathbf{z}(\mathbf{x}_{T_j}),-\mathbf{z}(\mathbf{x}_{T_j})\Bigr\rangle_{\mathbb{R}^d}d\mathbf{x}\Bigr|_{\epsilon=0} =\int _{\mathbf{\Psi }}\Bigl\langle\beta \sum _{j=1}^N\left(\mathbf{x}_{T_{\mathbf{u}j}}-\mathbf{x}_{R_j}\right),\mathbf{z}(\mathbf{x}_{T_j})\Bigr\rangle_{\mathbb{R}^d}d\mathbf{x}\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

where <math display="inline">\langle \cdot ,\cdot \rangle _{\mathbb{R}^d}</math> denotes the inner (or ''dot'') product in <math display="inline">\mathbb{R}^d</math>.

In order to solve the Euler-Lagrange equation ([[#eq-5|5]]), an iterative time-marching scheme can be used. For doing so, an artificial time <math display="inline">t</math> has to be added to the E-L equation and then the steady-state solution has to be computed, i.e., <math display="inline">\frac{\partial }{\partial t}\mathbf{u}+\mathbf{f}+\alpha \,(-\Delta )^\lambda \mathbf{u}=\mathbf{0}</math> (where <math display="inline">\mathbf{f}</math> gathers the first two addends in eq.([[#eq-5|5]])). In the steady-state, the displacement field <math display="inline">\mathbf{u}</math> converges and hence <math display="inline">\frac{\partial }{\partial t}\mathbf{u}=\mathbf{0}</math>. The time <math display="inline">t</math> is discretized, <math display="inline">t=\xi \,\tau </math> (with <math display="inline">\xi \in \mathbb{N}</math> being the iteration index and where <math display="inline">\tau{>0}</math> is the time-step), and finally the temporal derivative is replaced with its discrete approximation (first backward difference). Due to the fact that digital datasets are usually encoded by uniformly distributed spatial elements in each dimension (e.g. pixels if <math display="inline">d=2</math> or voxels if <math display="inline">d=3</math>), the discretization of the spatial variable becomes a natural approach too. Therefore in the following the notation <math display="inline">\mathbf{u}^{(\xi )}[\mathbf{n}]=\mathbf=\mathbf{n},\xi\,\tau</math> is used, where <math display="inline">\mathbf{n}=(n_1,\ldots ,n_d)</math> is the index of the discrete spatial position, with <math display="inline">n_i=0,\ldots \,N_i-1</math>, and <math display="inline">N_i</math> being the number of discrete spatial elements in the <math display="inline">i</math>-th dimension. Using this notation, the resulting semi-implicit iteration for the <math display="inline">l</math>-th component of the displacement field is the following:

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u_l^{(\xi )}[\mathbf{n}]=u_l^{(\xi{-1)}}[\mathbf{n}]-\tau \,f_l^{(\xi{-1)}}[\mathbf{n}]-\tau \,\alpha \,q[\mathbf{n}]\ast u_l^{(\xi )}[\mathbf{n}]\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

where the discrete operator <math display="inline">q[\mathbf{n}]</math> stands for a kernel which performs the discrete approximation of the spatial derivatives of <math display="inline">(-\Delta )^{\lambda }</math>.

As an alternative to the spatial approach, solving the iteration ([[#eq-12|12]]) in the frequency domain provides a stable implementation for the computation of a numerical solution for the displacement field <math display="inline">\mathbf=\mathbf{n}</math>, and in a more efficient way than other existing approaches if the <math display="inline">d</math>-dimensional fast Fourier transform (FFT) is taken into account <span id='citeF-40'></span>[[#cite-40|[40]]]. Translating the semi-implicit iteration ([[#eq-12|12]]) into the frequency domain results in

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>U_l^{(\xi )}(\boldsymbol{\omega })=U_l^{(\xi{-1)}}(\boldsymbol{\omega })-\tau \,F_l^{(\xi{-1)}}(\boldsymbol{\omega })-\tau \,\alpha \,Q(\boldsymbol{\omega })\,U_l^{(\xi )}(\boldsymbol{\omega })\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

where <math display="inline">U_l^{(\xi )}</math>, <math display="inline">U_l^{(\xi{-1)}}</math> and <math display="inline">F_l^{(\xi{-1)}}</math> are the <math display="inline">d</math>-dimensional Fourier transforms of <math display="inline">u_l^{(\xi )}</math>, <math display="inline">u_l^{(\xi{-1)}}</math> and <math display="inline">f_l^{(\xi{-1)}}</math>, respectively, and <math display="inline">\boldsymbol{\omega }=(\omega _1,\ldots ,\omega _d)</math> is the <math display="inline">d</math>-dimensional counterpart in the frequency domain associated to the discrete spatial variable <math display="inline">\mathbf{n}</math>. The operator <math display="inline">Q(\boldsymbol{\omega })</math> performs the spatial derivatives in the frequency domain and allows for their calculation by means of the product of spectra (i.e., the translation of the spatial convolution into the frequency domain); as deduced in <span id='citeF-40'></span>[[#cite-40|[40]]], its analytical expression is the following:

<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>Q(\boldsymbol{\omega })=\Bigl(2\sum _{m=1}^d(1-\cos \omega _m)\Bigr)^{\lambda }\,. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

It should be noted that, mathematically, any value <math display="inline">\lambda \in [1,2]</math> makes sense in the frequency domain (but not in the spatial domain). This allows for the design of ''hybrid'' diffusion-curvature regularizers with adaptable properties, please refer to <span id='citeF-25'></span>[[#cite-25|[25]]] for further details.

As explained in previous paragraphs, it is taken into account that the datasets are discrete and then the spatial variable <math display="inline">\mathbf{x}</math> gives rise to the discrete spatial index <math display="inline">\mathbf{n}</math>. Similarly, instead of handling continuous spectra, the frequency domain is also discretized and only the <math display="inline">N_m</math> uniformly sampled frequencies <math display="inline">\omega _m=\bigl(0,\frac{2\pi }{N_m},\ldots ,\frac{2\pi }{N_m}(N_m-1)\bigr)</math> are evaluated in each dimension. This way the iteration ([[#eq-13|13]]) can be implemented efficiently using the fast algorithms for the computation of the DFT which are provided by many programming languages (e.g., MATLAB's <math display="inline">fftn</math> built-in function which is based on the FFTW library by <span id='citeF-16'></span>[[#cite-16|[16]]]).

Finally, the iteration of the proposed registration algorithm is the following:

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u_l^{(\xi )}[\mathbf{n}]=\mathrm{IFFT}\left\{H[\mathbf{k}]\,\,\mathrm{FFT}\{{}u_l^{(\xi{-1)}}[\mathbf{n}]-\tau \,f_l^{(\xi{-1)}}[\mathbf{n}]\}{}\right\}\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

where <math display="inline">H[\mathbf{k}]=(1+\tau \,\alpha \,Q[\mathbf{k}])^{-1}</math>, <math display="inline">Q[\mathbf{k}]=\bigl(2\sum _{m=1}^d(1-\cos (\frac{2\pi }{N_m}k_m))\bigr)^{\lambda }</math>, <math display="inline">\mathbf{k}=(k_1,\ldots ,k_d)</math>, and <math display="inline">k_m=0,\ldots ,N_m-1</math>. Due to the fact that <math display="inline">Q[\mathbf{k}]</math> results in a circulant block matrix, all the spectra products and pseudo-inversions become Hadamard (i.e., pointwise) products and divisions, respectively <span id='citeF-10'></span>[[#cite-10|[10]]].

==3 Results==

In order to assess the performance of the proposed methodology, it is tested on 15 medical experiments involving triple-phase 3D computed tomography (CT) volumes of the liver under contrast agent injection. Although the acquisition of the different phases is continuous in most cases, there is no exact correspondence between them, so they must be registered in order to show the results in a common 3D framework. First, a study without contrast agent is acquired. Once the contrast agent is injected, it reaches the arteries (arterial phase), then the portal veins and next the hepatic veins. Hepatic and portal veins are then visible in the portal venous phase. Image registration allows to locate exactly the vessels in 3D at each phase of contrast-enhanced CT data in order to measure distances and volumes and provide objective parameters of the pathology which facilitate comparisons between patients, the tracking of tumors <span id='citeF-9'></span>[[#cite-9|[9]]], to make calculations on the volume of the liver to be preserved prior to a liver resection <span id='citeF-12'></span>[[#cite-12|[12]]], or to generate vessel models for planning surgical procedures <span id='citeF-22'></span>[[#cite-22|[22]]]. The experiments carried out involve data from 5 patients. The acquisition device is a GE LightSpeed VTC (General Electric Medical Systems). The size of the volumes vary from <math display="inline">512\times{512}\times{34}</math> to <math display="inline">512\times{512}\times{51}</math> voxels, with a resolution of <math display="inline">0.9102\times{0.9102}\times{5}</math> millimeters.

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Review_341874659107-MI.png|100px|Mean similarity between images in terms of mutual information for the three considered scenarios (arterial-portal, arterial-non-contrast and portal-non-contrast) before and after the registration.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' Mean similarity between images in terms of mutual information for the three considered scenarios (arterial-portal, arterial-non-contrast and portal-non-contrast) before and after the registration.
|}

<div id='img-2a'></div>
<div id='img-2b'></div>
<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Review_341874659107-error_box.png|100px|Box plot of the registration error]]
|[[Image:Review_341874659107-error_box_zoom.png|100px|Close-up view of (a)]]
|- style="text-align: center; font-size: 75%;"
| (a) Box plot of the registration error
| (b) Close-up view of (a)
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 2:''' Spatial distances (in millimeters) between corresponding landmarks before and after the registration process.
|}

<div id='img-3a'></div>
<div id='img-3b'></div>
<div id='img-3c'></div>
<div id='img-3d'></div>
<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Review_341874659107-gimenez_art_vac_R.png|100px|Arterial phase (reference dataset, R)]]
|[[Image:Review_341874659107-gimenez_art_vac_T.png|100px|Non-contrast phase (template dataset, T)]]
|[[Image:Review_341874659107-gimenez_art_vac_Tu_CR.png|100px|Registered template, T_u (CR-based method)]]
|[[Image:Review_341874659107-gimenez_art_vac_Tu_SSIM_LM.png|100px|Registered template, T_u (proposed method)]]
|- style="text-align: center; font-size: 75%;"
| (a) Arterial phase (reference dataset, <math display="inline">R</math>)
| (b) Non-contrast phase (template dataset, <math display="inline">T</math>)
| (c) Registered template, <math display="inline">T_\mathbf{u}</math> (CR-based method)
| (d) Registered template, <math display="inline">T_\mathbf{u}</math> (proposed method)
|- style="text-align: center; font-size: 75%;"
| colspan="4" | '''Figure 3:''' Visual outcomes of experiment 1 (slice 9): registration of arterial and non-contrast phases of patient 1.
|}

<div id='img-4b'></div>
<div id='img-4c'></div>
<div id='img-4d'></div>
<div id='img-4'></div>
<div id='img-4a'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Review_341874659107-lopez_art_por_T.png|100px|Portal venous phase (template dataset, T)]]
|[[Image:Review_341874659107-lopez_art_por_Tu_CR.png|100px|Registered template, T_u (CR-based method)]]
|[[Image:Review_341874659107-lopez_art_por_Tu_SSIM_LM.png|100px|Registered template, T_u (proposed method)Visual outcomes of experiment 5 (slice 12): registration of arterial and portal venous phases of patient 2.]]
|[[Image:Review_341874659107-lopez_art_vac_R.png|100px|Arterial phase (reference dataset, R)]]
|- style="text-align: center; font-size: 75%;"
| (b) Portal venous phase (template dataset, <math display="inline">T</math>)
| (c) Registered template, <math display="inline">T_\mathbf{u}</math> (CR-based method)
| (d) Registered template, <math display="inline">T_\mathbf{u}</math> (proposed method)
|- style="text-align: center; font-size: 75%;"
| colspan="4" | '''Figure a:''' Arterial phase (reference dataset, <math display="inline">R</math>)
|}

For the assessment of the proposed method, we carry out a comparison with the purely intensity-based variational method recently presented by the authors of this work in <span id='citeF-26'></span>[[#cite-26|[26]]]. This CR-based approach was reported to outperform publicly available state-of-the art methods such as ''Elastix'' and ANTs in the medical setting. As can be seen in Fig.[[#img-1|1]], the proposed framework shows excellent results for the three considered registration scenarios (arterial-portal, arterial-non-contrast and portal-non-contrast), reaching average values of 1.47, 1.44 and 1.52 bits in terms of mutual information, corresponding to the arterial-portal, arterial-non-contrast and portal-non-contrast cases, respectively; this represents a mean improvement of 28.9%, 48.45% and 51.16% in relative terms of mutual information, thus outperforming the CR-based registration algorithm, which achieves a mean improvement of 26.48%, 44.22% and 43.25%, respectively. Additionally, due to the analogous behavior (i.e., comparable final values of mutual information) of the proposed method in the three scenarios, all available experiments can be grouped into one ensemble in order to assess a more comprehensive validation of the actual registration error. A ground truth was established by an expert in the form of identifiable anatomical locations (landmarks) for all experiments. The registration errors were then obtained by computing the spatial distance between the corresponding landmarks in the reference and registered template datasets. Figures Fig.[[#img-2|2]](a) and Fig.[[#img-2|2]](b) show through box plots the registration error (in millimeters) achieved by the methods under comparison, gathering the results from the three considered registration scenarios. These box plots collect the final spatial distances between corresponding landmarks, along with the median distance error and its statistical significance (notch showing the 95% confidence interval of the true median). According to Fig.[[#img-2|2]](b), the proposed method significantly improves on the registration error of the CR-based approach, since it reduces the initial median error from 9.50 mm to a residual median distance between landmarks of 1.41 mm, decreasing at the same time the outliers occurrence.

In addition to the previous measurements, the visual outcomes of two of the experiments are shown in figures Fig.[[#img-3|3]] and Fig.[[#img-4|4]], whose purpose is to highlight the most illustrative differences (from a medical point of view) between the results provided by the compared methods. In Fig.[[#img-3|3]], we observe a normal size of the liver, with discretely irregular contours and homogeneous signal intensity. In hepatic segment II, there is a lesion of 40 mm of maximum axis, encapsulated and with well-defined contours and heterogeneous enhancement in arterial phase (after administration of intravenous contrast), suggestive of hepatocellular carcinoma (HCC). In this slice of the CT scan, we can also observe the aorta that shines in the arterial phase, the lower area of the stomach and the upper area of the spleen. In Fig.[[#img-4|4]], the liver has a normal size with discretely irregular contours in relation to changes due to chronic liver disease. In hepatic segment IV, a 36 mm diameter focal lesion is identified, which has arterial phase enhancement with a small area of necrosis of 13 mm; it corresponds to a HCC previously chemoembolized with partial necrosis. In this slice of CT, we can also observe the aorta, the gastric chamber and the spleen. When comparing the two methods under study, it can be seen how in Fig.[[#img-3|3]] the resulting registered datasets are very similar. However, looking closely, it can be noticed that in the right part of the image (left side of the patient) the shape and width of the structures corresponding to the stomach and the spleen in Fig.[[#img-3|3]](d) match better those in the reference dataset. Likewise, the part of the rib at the upper right of the image is more similar to the same region in the reference dataset by using the proposed method. Regarding the experiment shown in Fig.[[#img-4|4]], it can be easily appreciated how the geometrical matching (with respect to the reference dataset, Fig.[[#img-4|4]](a)) of the structures in the right side of the image (specially the gastric chamber) is visually more satisfactory in Fig.[[#img-4|4]](d). Moreover, the area of tumor necrosis which results from the proposed method is also slightly better aligned.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 Mean execution time per iteration (in seconds) of the considered variational methods and ''Elastix'', involving datasets of different sizes. These times were obtained on a PC with Intel Core i5-2500K, 3.3 GHz, 16 GB RAM, running Windows 8.1 (64 bits). All methods have been implemented in C++.
|- style="border-top: 2px solid;"
|       Size                      
| CR     
| SSIM+LM 
| ''Elastix'' 
|- style="border-top: 2px solid;"
| <math display="inline">512\times{512}\times{16}</math>    
| 2.40   
| 2.52    
| 5.69 
|-
| <math>512\times{512}\times{32}</math>
| 5.03   
| 5.27    
| 11.20 
|- style="border-bottom: 2px solid;"
| <math>512\times{512}\times{64}</math>
| 10.42  
| 10.93   
| 22.01 

|}

Throughout this work, the SSIM-based term uses the following parameters setting: <math display="inline">C_1=0.01</math> and <math display="inline">C_2=0.03</math>. These are the empirically obtained values suggested in the reference paper <span id='citeF-43'></span>[[#cite-43|[43]]]. Regardless, we find that in the current scenario, the performance of the proposed registration algorithm is fairly insensitive to variations of these values. Regarding the remaining parameters, which are exclusively related to our variational approach, the values of <math display="inline">\lambda=2</math> (i.e., curvature smoother, since the deformations to be corrected are notable in some cases), <math display="inline">\alpha=200</math>, <math display="inline">\beta=5</math> (only applicable to the combined SSIM- and landmark-based method) and <math display="inline">\tau=1</math> were used for all experiments |we recommend a common set which achieves a good balance between generalization and performance|. This parameter setting was obtained following the guidelines first introduced in <span id='citeF-23'></span>[[#cite-23|[23]]]. As for the number of iterations, a value of <math display="inline">\xi _{\max }=80</math> grants convergence in all cases; the cost function stabilizes after 50&#8211;70 iterations for both the CR-based and the proposed algorithm. From the results gathered in Table [[#table-1|1]], it can be concluded that the computational overhead introduced by the novel approach does not increase significantly the execution time of the CR-based algorithm. Moreover, both variational approaches outperform well-established registration methods such as ''Elastix'' in terms of efficiency. ''Elastix'' was the fastest method which provided better results among the open source methods of image registration which entered in the second phase of Empire10 Challenge <span id='citeF-30'></span>[[#cite-30|[30]]]. The overall complexity of each iteration of the compared algorithms is <math display="inline">\mathcal{O}(N)</math> |where <math display="inline">N</math> is the number of voxels of the datasets|, since doubling <math display="inline">N</math> means doubling the computational time.

==4 Conclusions and discussion==

In this work, a theoretical framework has been presented for approaching the image registration problem. In particular, we have formalized the variational formulation of image registration, making it valid even for a general <math display="inline">d</math>-dimensional case. The resulting semi-implicit iteration, which is solved in the frequency domain, allows for an efficient implementation if fast algorithms for the computation of the DFT are taken into account. The main novelty of this paper lies in the inclusion within such a framework of a disparity energy which combines the information provided by both the structural similarity index (SSIM) and the location of geometrical identifiable points (landmarks). With this purpose, the corresponding external forces fields of the resulting Euler-Lagrange equation have been deduced, and their analytical expressions are explicitly provided.

The suitability of the proposed methodology in the medical setting has been validated by means of illustrative experiments involving liver CT data under different contrast agent injection. When compared with an intensity-based method which was recently reported to outperform other popular state-of-the-art methods, it has been shown that the novel approach achieves higher values for both the similarity measure considered (mutual information) and the actual registration error (distance in space between corresponding landmarks). Moreover, the results provided by the method proposed in this paper are subjectively considered more satisfying after being visually inspected by an expert. However, the main advantage of the proposed registration method is its efficiency. As expected, the FFT-based, C++ implementation of our approach outperforms fast methods such as ''Elastix'' in terms of computational cost, while keeping the execution times in the range of its main contestant (the previously proposed CR-based algorithm). In a clinical environment, where computational times should be kept as low as possible, this feature can be of paramount importance. Regardless, it should be noted that comparing the performance of different image registration algorithms is not a trivial task, since each approach has its own set of user-defined parameters, which may heavily influence the final outcome. For instance, if they are tested on datasets with different characteristics (e.g., modality or anatomical region), new optimal parameters need to be determined. In this sense, open challenges are an interesting way to categorize image registration software packages in a common context.

It should be noted that the proposed method considers the input directional field as a periodic signal. Unless this is taken into account, the vectors located near the boundaries will be erroneous. As stated by other authors |see e.g. <span id='citeF-13'></span>[[#cite-13|[13]]]|, this is hardly a problem if the data is contained within a uniform background (e.g., when dealing with medical images or volumes, as in the current scenario), or else it can be overcome by using a folded algorithm which extends the dataset symmetrically, thus resulting in additional reflections of the original input, see e.g. <span id='citeF-3'></span>[[#cite-3|[3]]].

The most notable limitation of the present method is that the SSIM-based term would not be an appropriate choice for the combined disparity energy unless matching voxels of the datasets have similar intensities (i.e., in pure monomodal or pseudo-monomodal scenarios). In fact, this is an inherent drawback of the structural similarity index, not a problem which arises with our methodology. However, even if the aim was to register multimodal datasets, one could still use the proposed variational approach, but considering an alternative distance term instead (e.g., a combined CR- and landmark-based disparity energy).

In ongoing research, we would like to address the application of our methodology to other scenarios, such as the fusion of datasets from different modalities: for example, anatomical (CT and MRI), metabolic (spectroscopic MRI and PET) and functional (f-MRI) volumes in the medical setting. Alternative definitions of the joint functional, including additional energy terms |e.g., a measure of the differences between forward and reverse transformations, so that the resulting displacement field is ''consistent''|, are currently being explored within the presented theoretical framework.

figuretableAppendix &nbsp;equation

==5 Variational formulation of SSIM==

In this work, we propose the opposite of the structural similarity index as the intensity-based part of the disparity term of the joint energy functional (first addend in equation ([[#eq-2|2]])). Thus, a maximization problem is transformed into the minimization of the following cost function:

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{D}^{\mbox{SSIM}}[R,T;\mathbf{u}]:=-\mathrm=R,T_\mathbf{u}=\frac{-(2\,\mu _R\,\mu _T+C_1)(2\,\sigma _{RT}+C_2)}{(\mu _R^2+\mu _T^2+C_1)(\sigma _{RR}+\sigma _{TT}+C_2)}\,. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

The mathematical expectations, variances and covariance <math display="inline">\mu _R</math>, <math display="inline">\mu _T</math>, <math display="inline">\sigma _{RR}</math>, <math display="inline">\sigma _{TT}</math> and <math display="inline">\sigma _{RT}</math> are defined as follows:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mu _R=\mathrm{E}\{{\lbracechar}R\}{\rbracechar}=\int _\mathcal{I}i_r\,p_R(i_r)di_r\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mu _T=\mathrm{E}\{{\lbracechar}T_\mathbf{u}\}{\rbracechar}=\int _\mathcal{I}i_t\,p_{T_\mathbf{u}}(i_t)di_t\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma _{RR}=\mathrm{Var}\{{\lbracechar}R\}{\rbracechar}=\int _\mathcal{I}i^2_r\,p_R(i_r)di_r-\mathrm{E}^2\{{\lbracechar}R\}{\rbracechar}\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma _{TT}=\mathrm{Var}\{{\lbracechar}T_\mathbf{u}\}{\rbracechar}=\int _\mathcal{I}i^2_t\,p_{T_\mathbf{u}}(i_t)di_t-\mathrm{E}^2\{{\lbracechar}T_\mathbf{u}\}{\rbracechar}\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma _{RT}=\mathrm{Cov}\{{\lbracechar}R,T_\mathbf{u}\}{\rbracechar}=\int _{\mathcal{I}^2}i_r\,i_t\,p_{R,T_\mathbf{u}}(i_r,i_t)di_r\,di_t-\mathrm{E}\{{\lbracechar}R\}{\rbracechar}\,\mathrm{E}\{{\lbracechar}T_\mathbf{u}\}{\rbracechar}\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

where <math display="inline">p_R</math> and <math display="inline">p_{T_\mathbf{u}}</math> denote the marginal intensity distributions estimated from <math display="inline">R</math> and <math display="inline">T_\mathbf{u}</math>, respectively, and <math display="inline">p_{R,T_\mathbf{u}}</math> stands for an estimate of the joint intensity distribution. The intensities of the datasets, <math display="inline">i_r\,,i_t\in \mathcal{I}\subset \mathbb{R}</math>, are considered as random variables whose probability densities are respectively given by <math display="inline">p_R</math> and <math display="inline">p_{T_\mathbf{u}}</math>.

In order to obtain the external forces field related to the SSIM-based energy term (i.e., the first addend in the Euler-Lagrange equation ([[#eq-5|5]])), the computation of the corresponding ''Gateaux'' derivative <math display="inline">\mathrm{d}\mathcal{D}^{\mbox{SSIM}}</math> is required. Therefore, we carry out an explicit computation from the definition of the ''Gateaux'' derivative (see eq.([[#eq-4|4]])):

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathrm{d}\mathcal{D}^{\mbox{SSIM}}[R,T;\mathbf{u};\mathbf{z}]=\frac{\partial }{\partial \epsilon }\left.\mathcal{D}^{\mbox{SSIM}}[R,T;\mathbf{u}+\epsilon \mathbf{z}]\right|_{\epsilon=0}</math>
|-
| style="text-align: center;" | <math> =-\frac{(\mu _R^2+\mu _T^2+C_1)(\sigma _{RR}+\sigma _{TT}+C_2)(2\,\mu _R\,\mu _T'(2\,\sigma _{RT}+C_2)+2\,\sigma _{RT}'(2\,\mu _R\,\mu _T+C_1))}{(\mu _R^2+\mu _T^2+C_1)^2(\sigma _{RR}+\sigma _{TT}+C_2)^2}</math>
|-
| style="text-align: center;" | <math> +\left.\frac{(2\,\mu _R\,\mu _T+C_1)(2\,\sigma _{RT}+C_2)(2\,\mu _T\,\mu _T'(\sigma _{RR}+\sigma _{TT}+C_2)+\sigma _{TT}'(\mu _R^2+\mu _T^2+C_1))}{(\mu _R^2+\mu _T^2+C_1)^2(\sigma _{RR}+\sigma _{TT}+C_2)^2}\right|_{\epsilon=0}</math>
|-
| style="text-align: center;" | <math> =\int _{\mathcal{I}^2}\left.L^{\mbox{SSIM}}_{\mathbf{u}+\epsilon \mathbf{z}}(i_r,i_t)\frac{\partial }{\partial \epsilon }p_{R,T_{\mathbf{u}+\epsilon \mathbf{z}}}(i_r,i_t)di_r\,di_t\right|_{\epsilon=0}\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

with <math display="inline">L^{\mbox{SSIM}}_{\mathbf{u}+\epsilon \mathbf{z}}</math> being an intensity operator |whose definition, not necessary at this point, is not explicitly provided for the sake of simplicity|, and where the notation <math display="inline">(\cdot )'=\frac{\partial }{\partial \epsilon }(\cdot )</math> is used.

Assuming a displacement <math display="inline">\mathbf{u}+\epsilon \mathbf{z}</math>, the joint intensity distribution estimated from <math display="inline">R</math> and <math display="inline">T_{\mathbf{u}+\epsilon \mathbf{z}}</math> is provided by a non-parametric Parzen-Rozenblatt density model <span id='citeF-11'></span>[[#cite-11|[11]]]:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p_{R,T_{\mathbf{u}+\epsilon \mathbf{z}}}(i_r,i_t)=\frac{1}{\mathcal{V}_\mathbf{\Psi }}\int _\mathbf{\Psi }\hat{G}_\rho \left(R(\mathbf{x})-i_r,T(\mathbf{x}+\mathbf{u}(\mathbf{x})+\epsilon \mathbf{z}(\mathbf{x})-i_t)\right)d\mathbf{x}\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

where <math display="inline">\mathcal{V}_\mathbf{\Psi }</math> denotes the hypervolume of <math display="inline">\mathbf{\Psi }</math> and <math display="inline">\hat{G}_\rho :\mathcal{I}^2\rightarrow \mathbb{R}</math> is a Gaussian kernel (strictly positive and differentiable):

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hat{G}_\rho (i_r,i_t)=\frac{1}{2\pi \rho ^2}\,e^{-\frac{i_r^2+i_t^2}{2\rho ^2}}\,. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

At this point, derivatives of <math display="inline">p_{R,T_{\mathbf{u}+\epsilon \mathbf{z}}}</math> can be easily computed. In particular,

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{\partial }{\partial \epsilon }p_{R,T_{\mathbf{u}+\epsilon \mathbf{z}}}(i_r,i_t)= \frac{1}{\mathcal{V}_\mathbf{\Psi }}\int _\mathbf{\Psi }\frac{\partial }{\partial  i_t}\hat{G}_\rho \left(R(\mathbf{x})-i_r,T\left(\mathbf{x}+\mathbf{u}(\mathbf{x})+\epsilon \mathbf{z}(\mathbf{x})\right)-i_t\right)\nabla T\left(\mathbf{x}+\mathbf{u}(\mathbf{x})+\epsilon \mathbf{z}(\mathbf{x})\right)\mathbf{z}(\mathbf{x})d\mathbf{x}\,. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

We now replace ([[#eq-10|10]]) in ([[#eq-7|7]]), and then let <math display="inline">\epsilon=0</math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathrm{d}\mathcal{D}^{\mbox{SSIM}}[R,T;\mathbf{u};\mathbf{z}]=\frac{1}{\mathcal{V}_\mathbf{\Psi }}\int _\mathbf{\Psi }\int _{\mathcal{I}^2} L_\mathbf{u}^{\mbox{SSIM}}(i_r,i_t)\frac{\partial }{\partial  i_t}\hat{G}_\rho \left(R(\mathbf{x})-i_r,T\left(\mathbf{x}+\mathbf{u}(\mathbf{x})\right)-i_t\right)\nabla T\left(\mathbf{x}+\mathbf{u}(\mathbf{x})\right)\mathbf{z}(\mathbf{x})di_r\,di_t\,d\mathbf{x}</math>
|-
| style="text-align: center;" | <math> =\int _\mathbf{\Psi }\Bigl\langle\frac{1}{\mathcal{V}_\mathbf{\Psi }}\left(\hat{G}_\rho \ast \frac{\partial L_\mathbf{u}^{\mbox{ SSIM}}}{\partial i_t}\Bigl(R(\mathbf{x}),T\left(\mathbf{x}+\mathbf{u}(\mathbf{x})\right)\Bigr)\right)\nabla  T\left(\mathbf{x}+\mathbf{u}(\mathbf{x})\right),\mathbf{z}(\mathbf{x})\Bigr\rangle_{\mathbb{R}^d}d\mathbf{x}</math>
|-
| style="text-align: center;" | <math> =\int _\mathbf{\Psi }\bigl\langle\mathbf{f}^{\mbox{ SSIM}}(\mathbf{x};\mathbf{u}),\mathbf{z}(\mathbf{x})\bigr\rangle_{\mathbb{R}^d}d\mathbf{x}\,. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

In the previous equation, a bidimensional convolution with respect to the intensity variable <math display="inline">\mathbf{i}=(i_r,i_t)</math> appears naturally. This convolution conmutes with the derivative <math display="inline">\frac{\partial }{\partial i_t}</math> (since both are linear operators). Moreover, by identifying the resulting expression with an inner product in <math display="inline">\mathbb{R}^d</math>, the external forces field of the Euler-Lagrange equation ([[#eq-5|5]]) which corresponds to the SSIM-based energy term can be finally obtained:

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{f}^{\mbox{SSIM}}=\frac{1}{\mathcal{V}_\mathbf{\Psi }}\left(\hat{G}_\rho \ast \frac{\partial L_\mathbf{u}^{\mbox{SSIM}}}{\partial i_t}\right)\nabla T_{\mathbf{u}}\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

with

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{\partial L_\mathbf{u}^{\mbox{SSIM}}}{\partial i_t}=\frac{\gamma _0+C_1\,\gamma _1+C_2\,\gamma _2+C_1C_2\,\gamma _3}{(\mu _R^2+\mu _T^2+C_1)^2(\sigma _{RR}+\sigma _{TT}+C_2)^2}\,, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

and where the definitions of <math display="inline">\gamma _0</math>, <math display="inline">\gamma _1</math>, <math display="inline">\gamma _2</math> and <math display="inline">\gamma _3</math> can be found in eqs.([[#eq-7|7]])-([[#eq-10|10]]).

===BIBLIOGRAPHY===

<code>#1</code> [ ]#2 #1 doi: arXiv: URL:  pmid: [http://dx.doi.org/#1 #1] [pmid:#1 #1] #2  <div id="cite-1"></div>
'''[[#citeF-1|[1]]]'''  authorAmintoosi, M., authorFathy, M.,   authorMozayani, N., year2009. titlePrecise image registration with structural similarity   error measurement applied to superresolution. journalEURASIP Journal on Advances in Signal Processing   volume2009, pages305479. <div id="cite-2"></div>
'''[[#citeF-2|[2]]]'''  authorAmit, Y., year1994. titleA nonlinear variational problem for image matching. journalSIAM Journal of Scientific Computing   volume15, pages207&#8211;224. <div id="cite-3"></div>
'''[[#citeF-3|[3]]]'''  authorBai, J., authorFeng, X.C., year2007. titleFractional-order anisotropic diffusion for image   denoising. journalIEEE Transactions on Image Processing   volume16, pages2492&#8211;2502. <div id="cite-4"></div>
'''[[#citeF-4|[4]]]'''  authorBen&nbsp;Sassi, O., authorDelleji, T.,   authorTaleb-Ahmed, A., authorFeki, I.,   authorBen&nbsp;Hamida, A., year2008. titleMr image monomodal registration using structure   similarity index, in: booktitleImage Processing Theory, Tools and   Applications, 2008. IPTA 2008. First Workshops on, pp.   pages1&#8211;5. <div id="cite-5"></div>
'''[[#citeF-5|[5]]]'''  authorBookstein, F.L., year1989. titlePrincipal warps: thin-plate splines and the   decomposition of deformations. journalIEEE Transactions on Pattern Analysis and Machine   Intelligence volume11, pages567&#8211;585. <div id="cite-6"></div>
'''[[#citeF-6|[6]]]'''  authorBrown, L.G., year1992. titleA survey of image registration techniques. journalACM Computing Surveys volume24,   pages325&#8211;376. <div id="cite-7"></div>
'''[[#citeF-7|[7]]]'''  authorBrunet, D., authorVrscay, E.,   authorWang, Z., year2012. titleOn the mathematical properties of the structural   similarity index. journalImage Processing, IEEE Transactions on   volume21, pages1488&#8211;1499. <div id="cite-8"></div>
'''[[#citeF-8|[8]]]'''  authorCachier, P., authorBardinet, E.,   authorDormont, D., authorPennec, X.,   authorAyache, N., year2003. titleIconic feature based nonrigid registration: the   PASHA algorithm. journalComputer Vision and Image Understanding   volume89, pages272&#8211;298. <div id="cite-9"></div>
'''[[#citeF-9|[9]]]'''  authorCharnoz, A., authorAgnus, V.,   authorMalandain, G., authorForest, C.,   authorTajine, M., authorSoler, L.,   year2005. titleLiver registration for the follow-up of hepatic   tumors, in: editorDuncan, J., editorGerig, G.   (Eds.), booktitleMedical Image Computing and Computer-Assisted   Intervention MICCAI. publisherSpringer Berlin / Heidelberg.   volume volume3750 of ''seriesLecture Notes in   Computer Science'', pp. pages155&#8211;162. <div id="cite-10"></div>
'''[[#citeF-10|[10]]]'''  authorDavis, P.J., year1979. titleCirculant Matrices. publisherWiley-Interscience, NY. <div id="cite-11"></div>
'''[[#citeF-11|[11]]]'''  authorDuda, R., authorHart, P., year1973. titlePattern Classification and Scene Analysis. publisherWiley. <div id="cite-12"></div>
'''[[#citeF-12|[12]]]'''  authorElhawary, H., authorOguro, S.,   authorTuncali, K., authorMorrison, P.,   authorShyn, P., authorTatli, S.,   authorSilverman, S., authorHata, N.,   year2009. titleIntra-operative multimodal non-rigid registration of   the liver for navigated tumor ablation, in: editorYang, G.Z.,   editorHawkes, D., editorRueckert, D.,   editorNoble, A., editorTaylor, C. (Eds.),   booktitleMedical Image Computing and Computer-Assisted   Intervention MICCAI. publisherSpringer Berlin / Heidelberg.   volume volume5761 of ''seriesLecture Notes in   Computer Science'', pp. pages837&#8211;844. <div id="cite-13"></div>
'''[[#citeF-13|[13]]]'''  authorFischer, B., authorModersitzki, J.,   year1999. titleFast inversion of matrices arising in image   processing. journalNumerical Algorithms volume22,   pages1&#8211;11. <div id="cite-14"></div>
'''[[#citeF-14|[14]]]'''  authorFischer, B., authorModersitzki, J.,   year2002. titleFast diffusion registration. journalM.Z. Nashed, O. Scherzer (eds), Contemporary   Mathematics 313, Inverse Problems, Image Analysis, and Medical Imaging,   AMS , pages117&#8211;129. <div id="cite-15"></div>
'''[[#citeF-15|[15]]]'''  authorFischer, B., authorModersitzki, J.,   year2004. titleA unified approach to fast image registration and a   new curvature based registration technique. journalLinear Algebra and its Applications   volume308, pages107&#8211;124. <div id="cite-16"></div>
'''[[#citeF-16|[16]]]'''  authorFrigo, M., authorJohnson, S.G.,   year2005. titleThe design and implementation of FFTW3. journalProceedings IEEE volume93,   pages216&#8211;231. <div id="cite-17"></div>
'''[[#citeF-17|[17]]]'''  authorHuang, Y.H., authorOu, T.S., authorSu,   P.Y., authorChen, H., year2010. titlePerceptual rate-distortion optimization using   structural similarity index as quality metric. journalCircuits and Systems for Video Technology, IEEE   Transactions on volume20, pages1614&#8211;1624. <div id="cite-18"></div>
'''[[#citeF-18|[18]]]'''  authorJohnson, H.J., authorChristensen, G.E.,   year2002. titleConsistent landmark and intensity-based image   registration. journalIEEE Transactions on Medical Imaging   volume21, pages450&#8211;461. <div id="cite-19"></div>
'''[[#citeF-19|[19]]]'''  authorKalinic, H., authorLoncaric, S.,   authorBijnens, B., year2011. titleA novel image similarity measure for image   registration, in: booktitleImage and Signal Processing and   Analysis (ISPA), 2011 7th International Symposium on, pp.   pages195&#8211;199. <div id="cite-20"></div>
'''[[#citeF-20|[20]]]'''  authorKlein, S., authorStaring, M.,   authorMurphy, K., authorViergever, M.A.,   authorPluim, J.P.W., year2010. titleelastix: A toolbox for intensity-based medical image   registration. journalIEEE Transactions on Medical Imaging   volume29, pages196&#8211;205. <div id="cite-21"></div>
'''[[#citeF-21|[21]]]'''  authorKolesov, I., authorLee, J.,   authorSharp, G., authorVela, P.,   authorTannenbaum, A., year2016. titleA stochastic approach to diffeomorphic point set   registration with landmark constraints. journalIEEE Transactions on Pattern Analysis and Machine   Intelligence volume38, pages238&#8211;251. <div id="cite-22"></div>
'''[[#citeF-22|[22]]]'''  authorLange, T., authorWenckebach, T.H.,   authorLamecker, H., authorSeebass, M.,   authorHünerbein, M., authorEulenstein, S.,   authorGebauer, B., authorSchlag, P.M.,   year2005. titleRegistration of different phases of contrast-enhanced   ct/mri data for computer-assisted liver surgery planning: Evaluation of   state-of-the-art methods. journalThe International Journal of Medical Robotics and   Computer Assisted Surgery volume1, pages6&#8211;20. <div id="cite-23"></div>
'''[[#citeF-23|[23]]]'''  authorLarrey-Ruiz, J., authorMorales-Sánchez, J.,   year2006. titleOptimal parameters selection for non-parametric image   registration methods. journalLecture Notes in Computer Science   volume4179, pages564&#8211;575. <div id="cite-24"></div>
'''[[#citeF-24|[24]]]'''  authorLarrey-Ruiz, J., authorMorales-Sánchez, J.,   authorVerdú-Monedero, R., year2007. titleGeneralized regularization term for non-parametric   multimodal image registration. journalSignal Processing volume87,   pages2837&#8211;2842. <div id="cite-25"></div>
'''[[#citeF-25|[25]]]'''  authorLarrey-Ruiz, J., authorVerdú-Monedero, R.,   authorMorales-Sánchez, J., year2008. titleA fourier domain framework for variational image   registration. journalJournal of Mathematical Imaging and Vision   volume32, pages57&#8211;72. <div id="cite-26"></div>
'''[[#citeF-26|[26]]]'''  authorLegaz-Aparicio, A.G., authorVerdú-Monedero, R.,   authorLarrey-Ruiz, J., authorMorales-Sánchez, J.,   authorLópez-Mir, F., authorNaranjo, V.,   authorBernabéu, A., year2017. titleEfficient variational approach to multimodal   registration of anatomical and functional intra-patient tumorous brain data. journalInternational Journal of Neural Systems   volume(27), pages1&#8211;10. <div id="cite-27"></div>
'''[[#citeF-27|[27]]]'''  authorLoizou, C.P., authorTheofanous, C.,   authorPantziaris, M., authorKasparis, T.,   year2014. titleDespeckle filtering software toolbox for ultrasound   imaging of the common carotid artery. journalComputer Methods and Programs in Biomedicine   volume114, pages109 &#8211; 124. <div id="cite-28"></div>
'''[[#citeF-28|[28]]]'''  authorMaintz, J., authorViergever, M.,   year1998. titleA survey of medical image registration. journalMedical Image Analysis volume2,   pages1&#8211;36. <div id="cite-29"></div>
'''[[#citeF-29|[29]]]'''  authorModersitzki, J., year2004. titleNumerical Methods for Image Registration. publisherOxford University Press, NY. <div id="cite-30"></div>
'''[[#citeF-30|[30]]]'''  authorMurphy, K., authorvan Ginneken, B.,   authorReinhardt, J., authorKabus, S.,   authorDing, K., authorDeng, X., authorCao,   K., authorDu, K., authorChristensen, G.,   authorGarcia, V., authorVercauteren, T.,   authorAyache, N., authorCommowick, O.,   authorMalandain, G., authorGlocker, B.,   authorParagios, N., authorNavab, N.,   authorGorbunova, V., authorSporring, J.,   authorde&nbsp;Bruijne, M., authorHan, X.,   authorHeinrich, M., authorSchnabel, J.,   authorJenkinson, M., authorLorenz, C.,   authorModat, M., authorMcClelland, J.,   authorOurselin, S., authorMuenzing, S.,   authorViergever, M., authorDe&nbsp;Nigris, D.,   authorCollins, D., authorArbel, T.,   authorPeroni, M., authorLi, R.,   authorSharp, G., authorSchmidt-Richberg, A.,   authorEhrhardt, J., authorWerner, R.,   authorSmeets, D., authorLoeckx, D.,   authorSong, G., authorTustison, N.,   authorAvants, B., authorGee, J.,   authorStaring, M., authorKlein, S.,   authorStoel, B., authorUrschler, M.,   authorWerlberger, M., authorVandemeulebroucke, J.,   authorRit, S., authorSarrut, D.,   authorPluim, J., year2011. titleEvaluation of registration methods on thoracic ct:   The empire10 challenge. journalMedical Imaging, IEEE Transactions on   volume30, pages1901&#8211;1920. <div id="cite-31"></div>
'''[[#citeF-31|[31]]]'''  authorMyronenko, A., authorSong, X.,   year2010. titleIntensity-based image registration by minimizing   residual complexity. journalMedical Imaging, IEEE Transactions on   volume29, pages1882&#8211;1891. <div id="cite-32"></div>
'''[[#citeF-32|[32]]]'''  authorNaranjo, V., authorMorales, S.,   authorLegaz-Aparicio, A.G., authorLarrey-Ruiz, J.,   authorBernabéu, A., authorFuentes-Hurtado, F.,   year2015. titleA software for surgical and radiotherapy planning   through multimodal brain image registration and fusion. journalInternational Journal of Computer Assisted   Radiology and Surgery volume10, pages15&#8211;17. <div id="cite-33"></div>
'''[[#citeF-33|[33]]]'''  authorOu, Y., authorAkbari, H.,   authorBilello, M., authorDa, X.,   authorDavatzikos, C., year2014. titleComparative evaluation of registration algorithms in   different brain databases with varying difficulty: Results and insights. journalMedical Imaging, IEEE Transactions on   volume33, pages2039&#8211;2065. <div id="cite-34"></div>
'''[[#citeF-34|[34]]]'''  authorRehman, A., authorWang, Z.,   authorBrunet, D., authorVrscay, E.,   year2011. titleSsim-inspired image denoising using sparse   representations, in: booktitleAcoustics, Speech and Signal   Processing (ICASSP), 2011 IEEE International Conference on, pp.   pages1121&#8211;1124. <div id="cite-35"></div>
'''[[#citeF-35|[35]]]'''  authorRohr, K., year2001. titleLandmark-based image analysis: using geometric and   intensity models. journalComputational Imaging and Vision Series, Kluwer   Academic Publishers, Dordrecht volume21. <div id="cite-36"></div>
'''[[#citeF-36|[36]]]'''  authorSotiras, A., authorDavatzikos, C.,   authorParagios, N., year2013. titleDeformable medical image registration: A survey. journalMedical Imaging, IEEE Transactions on   volume32, pages1153&#8211;1190. <div id="cite-37"></div>
'''[[#citeF-37|[37]]]'''  authorSotiras, A., authorOu, Y.,   authorParagios, N., authorDavatzikos, C.,   year2015. titleGraph-based deformable image registration, in:   editorParagios, N., editorDuncan, J.,   editorAyache, N. (Eds.), booktitleHandbook of   Biomedical Imaging; Methodologies and Clinical Research.   publisherSpringer. Lecture Notes in Computer Science (Book 779). <div id="cite-38"></div>
'''[[#citeF-38|[38]]]'''  authorTosun, D., authorRettmann, M.E.,   authorPrince, J.L., year2004. titleMapping techniques for aligning sulci across multiple   brains. journalMedical Image Analysis volume8,   pages295&#8211;309. <div id="cite-39"></div>
'''[[#citeF-39|[39]]]'''  authorUrschler, M., authorZach, C.,   authorDitt, H., authorBischof, H.,   year2006. titleAutomatic point landmark matching for regularizing   nonlinear intensity registration: application to thoracic CT images. journalLecture Notes in Computer Science   volume4191, pages710&#8211;717. <div id="cite-40"></div>
'''[[#citeF-40|[40]]]'''  authorVerdú-Monedero, R., authorLarrey-Ruiz, J.,   authorMorales-Sánchez, J., year2008. titleFrequency implementation of the Euler-Lagrange   equations for variational image registration. journalIEEE Signal Processing Letters   volume15, pages321&#8211;324. <div id="cite-41"></div>
'''[[#citeF-41|[41]]]'''  authorVishnukumar, S., authorWilscy, M.,   year2009. titleA comparative study on adaptive local image   registration methods, in: booktitle2009 International Conference   on Signal Processing Systems, pp. pages140&#8211;145. <div id="cite-42"></div>
'''[[#citeF-42|[42]]]'''  authorWang, J., authorZhang, J.Q.,   year2009. titleAn iterative refinement dsa image registration   algorithm using structural image quality measure, in:   booktitleIntelligent Information Hiding and Multimedia Signal   Processing, 2009. IIH-MSP '09. Fifth International Conference on, pp.   pages973&#8211;976. <div id="cite-43"></div>
'''[[#citeF-43|[43]]]'''  authorWang, Z., authorBovik, A.C.,   authorSheikh, H.R., authorSimoncelli, E.P.,   year2004. titleImage quality assessment: from error visibility to   structural similarity. journalIEEE Transactions on Image Processing   volume13, pages600&#8211;612. <div id="cite-44"></div>
'''[[#citeF-44|[44]]]'''  authorWen, C., authorWang, D., authorShi,   L., authorChu, W.C.W., authorCheng, J.C.Y.,   authorLui, L.M., year2015. titleLandmark constrained registration of high-genus   surfaces applied to vestibular system morphometry. journalComputerized Medical Imaging and Graphics   volume44, pages1&#8211;12. <div id="cite-45"></div>
'''[[#citeF-45|[45]]]'''  authorWirth, M., authorBobier, B.,   year2009. titleRegistration of the prokudin-gorskii colour   photographs using a multiresolution ssim algorithm, in:   editorKamel, M., editorCampilho, A. (Eds.),   booktitleImage Analysis and Recognition.   publisherSpringer Berlin Heidelberg. volume   volume5627 of ''seriesLecture Notes in Computer   Science'', pp. pages917&#8211;926. <div id="cite-46"></div>
'''[[#citeF-46|[46]]]'''  authorZhang, Z., authorJiang, Y.,   authorTsui, H., year2006. titleConsistent multi-modal non-rigid registration based   on a variational approach. journalPattern Recognition Letters volume27,   pages715&#8211;725. <div id="cite-47"></div>
'''[[#citeF-47|[47]]]'''  authorZitová, B., authorFlusser, J.,   year2003. titleImage registration methods: a survey. journalImage and Vision Computing volume21,   pages997&#8211;1000.