Planar Registration is done in a 2D area. The displacement is first rotation, then translation. There are different types of displacement. First, there is the Rigid Linear Transform. It is a linear transformation that preserves the length of a point’s vector. The orientation of the latter may change, but the magnitude does not.

Then, there is the Affine Transformation which is an extension of the linear transform. The rotation in this case is Matrix which is multiplied to the moving points. Afterwards, the translation is added as a vector. This results in following formula $//<![CDATA[ \begin{array}{l}\vec{b} = R\vec{a}+\vec{d}\end{array} //]]>$ where $//<![CDATA[ \begin{array}{l}\vec{b}\end{array} //]]>$ is the moving image, $//<![CDATA[ \begin{array}{l}R\end{array} //]]>$ is the rotation matrix and $//<![CDATA[ \begin{array}{l}\vec{d}\end{array} //]]>$ is the translation. If the rotation matrix is the identity matrix the process is called pure translation because there is no rotation. When we now try to invert with given $//<![CDATA[ \begin{array}{l}\vec{b}\end{array} //]]>$ and solve for $//<![CDATA[ \begin{array}{l}\vec{a}\end{array} //]]>$, we get $//<![CDATA[ \begin{array}{l}\vec{a} = R^T \vec{b} - \left[ R^T \vec{d} \right]\end{array} //]]>$. So we can say that inverse displacements are also displacements.

There is also the Least Squares Translation. It is the vector between the means ofand. With the mean the zero-mean version of the data can be calculated:

Figure 2 shows some example data which is rotated by the rotation matrix $//<![CDATA[ \begin{array}{l}R\end{array} //]]>$ and translated by the vector $//<![CDATA[ \begin{array}{l}\vec{d}\end{array} //]]>$. With the least squares approach the data is converted to the zero-mean version, which can be seen in figure 3. Now the only thing that is left is the rotation. Therefore, you have to find the best $//<![CDATA[ \begin{array}{l}\Theta\end{array} //]]>$ such that $//<![CDATA[ \begin{array}{l}\vec{q_j} \approx \left[ \begin{matrix} cos(\Theta) & -sin(\Theta)\\ sin(\Theta) & cos(\Theta) \end{matrix} \right] \vec{p}\end{array} //]]>$. But this is more difficult than it looks because with optimization there is no guarantee of a global minimum. The estimation of the nearest orthogonal matrix is very hard. So we try to convert planar vectors to complex numbers. This was already discovered by Argand in 1806. He proposed to convert the zero-mean vectors to row vectors with complex entries:

So the rotation is just the multiplication by a complex number which results in $//<![CDATA[ \begin{array}{l}\vec{q} \approx r \vec{p}\end{array} //]]>$. The rotation can then be calculated by $//<![CDATA[ \begin{array}{l}r \approx \vec{q} \left[ \vec{p}^T \right]\end{array} //]]>$ [1].

Intensity-based Registration

The process is started with the moving image volume $//<![CDATA[ \begin{array}{l}I_m\end{array} //]]>$ on top left and the initial transformation $//<![CDATA[ \begin{array}{l}T\end{array} //]]>$. The moving image is transformed with the initial transformation: $//<![CDATA[ \begin{array}{l}T(I_m)=I'_m\end{array} //]]>$. Then the new moving image $//<![CDATA[ \begin{array}{l}I'_m\end{array} //]]>$ is compared to the fixed image volume $//<![CDATA[ \begin{array}{l}I_f\end{array} //]]>$. For this comparison different similarity measures are available:

• Sum of Squared Differences (SSD): $//<![CDATA[ \begin{array}{l}SSD = {1 \over N}\sum_{i \in \Omega} (I_f(i)-I'_m(i))^2\end{array} //]]>$

• Sum of Absolute Differences (SAD): $//<![CDATA[ \begin{array}{l}SAD={1\over N} \sum_{i \in \Omega} \left| I_f(i) - I'_m(i) \right|\end{array} //]]>$

• Normalized Cross Correlation (NCC): $//<![CDATA[ \begin{array}{l}NCC= {1\over N} {{\sum_{i\in\Omega}(I_f(i) - \bar{I_f})(I'_m(i)-\bar{I'_m})}\over{\sigma_{I_f}\sigma_{I'_m}}}\end{array} //]]>$

• Mutual Information (MI): $//<![CDATA[ \begin{array}{l}MI(X,Y)=\sum_i \sum_j p_{xy}(i,j) log {{p_{xy}(i,j)}\over{p_x(i)p_y(j)}}\end{array} //]]>$

SSD and SAD can only be used for intra-modal registration. This is the case because they totally depend on the intensities of the image. The NCC is not independent of the intensities but is already more detached. The last step of the iterative process is optimization of the transformation. In this step the new $//<![CDATA[ \begin{array}{l}T\end{array} //]]>$ for the transformation is calculated: $//<![CDATA[ \begin{array}{l}T^\star=\underset{T}{argmin}\:\sigma(I_f, T(I'_m))\end{array} //]]>$. For this problem there are many approaches like simplex method, gradient descent, newton’s method or least-squares method [1].

Demons Registration

For demons registration you need a source image $//<![CDATA[ \begin{array}{l}I_s\end{array} //]]>$ and a target image $//<![CDATA[ \begin{array}{l}I_t\end{array} //]]>$. The goal is to find a deformation field $//<![CDATA[ \begin{array}{l}u\end{array} //]]>$. The problem can be formulated as variational problem $//<![CDATA[ \begin{array}{l}\underset{u}{min}\:E(u)\end{array} //]]>$ where the energy $//<![CDATA[ \begin{array}{l}E\end{array} //]]>$ can be decomposed into a data term and a regularizer. The data term is based on the intesity of the source and target image. Using the diffusion equation the energy $//<![CDATA[ \begin{array}{l}E\end{array} //]]>$ can be minimized by $//<![CDATA[ \begin{array}{l}\partial u(x,t)=-(I_w(x)-I_t(x))\nabla I_w(x)+\lambda\Delta u(x)\end{array} //]]>$. By updating the equation with a Gaussian filter ($//<![CDATA[ \begin{array}{l}u(x,t+\tau)=G_{\sigma_e} * (u(x,t)+\tau(I_t(x)-I_w(x))\nabla I_w(x))\end{array} //]]>$ where $//<![CDATA[ \begin{array}{l}\sigma_e = 2 \sqrt{\lambda\tau}\end{array} //]]>$) the model assumes an elastic tissue deformation resulting a the so called elastic demons algorithm. If the Gaussian filter is applied before the update instead of after the update ($//<![CDATA[ \begin{array}{l}u(x,t+\tau)=u(x,t)+\tau G_{\sigma_f} * (I_t(x)-I_w(x))\nabla I_w(x)\end{array} //]]>$) the algorithm is called fluid demons. The combination of both versions is called elastic-fluid demons and uses this update equation $//<![CDATA[ \begin{array}{l}u(x,t+\tau)=G_{\sigma_e} * (u(x,t) + \tau G_{\sigma_f} * F(u, I_w,I_t))\end{array} //]]>$ where $//<![CDATA[ \begin{array}{l}F(u,I_w,I_t)=(I_t(x)-I_w(x))\nabla I_w(x)\end{array} //]]>$ [2].

For diffeomorphic demons you cannot use the additive updates. These have to change to compositive because diffeomorphisms are invertible functions that map one smooth surface to another such that both the function and its inverse are smooth. The update equation changes to $//<![CDATA[ \begin{array}{l}u(x,t+\tau)=G_{\sigma_e} * (u(x,t) \circ (Id + \tau G_{\sigma_f} * F(u, I_w;I_t)))\end{array} //]]>$. Diffeomorphisms are elements of a smooth but curved space (hypersurface). Simple addition does not work for diffeomorphisms as it may lead to non-invertible solutions (see figure). That is why computed forces have to be mapped with exponential map resulting in $//<![CDATA[ \begin{array}{l}u(x,t+\tau)=G_{\sigma_e} * (u(x,t) \circ exp(\tau G_{\sigma_f} * F(u, I_w;I_t)))\end{array} //]]>$ [2].