\section{Introduction}
Image guided radiation therapy has been a key component to improve the accuracy of radiation therapy \cite{xing2006overview}. Daily imaging of the current anatomical state informs whether to modify the treatment plan or not. With integration of MRI with the treatment machines is today possible to do real-time magnetic resonance imaging (MRI) during a treatment session \cite{paganelli2018mri}. Since MRI is ideally suited for imaging soft tissues, this allows identification of the shape and position of tumors and organs at risk and adjusting of the treatment accordingly. 

In this paper we will describe a new method for medical image registration. Our method is suitable for 2D sequential images and it is inspired by recent progress within the domain of computer vision. We first estimate a sparse displacement field using discriminative correlation filters (DCF)~\cite{danelljan2018learning} to track a fixed number of points in the image sequence. Then, we use a trained sparse-to-dense interpolation scheme in the form of a neural network that uses normalized convolutional layers~\cite{eldesokey2019confidence}. Normalized convolutions~\cite{knutsson1993normalized} use confidence and are well-suited for irregularly sampled data. Advantages of our method include i) it is modality agnostic due to the fact that the model is trained online and customized for the current image type, ii) it allows flexibility in trading off computation versus accuracy due to the freedom in selecting tracking points and iii) it is designed with real-time requirements in mind since each tracker is fast and the trackers are independent of each other. Possible disadvantages might be if a tracker loses its target and the filter is learn online to track a divergent point or if the tracking points are too sparse and movements are missed.

Even though conventional methods have shown good result they are computationally heavy as each registration problem is typically formulated as a large optimization problem that needs to be solved iteratively. Image registration methods can be categorized as physics-based methods, e.g. like Demons~\cite{thirion1998image}, interpolation-based methods or knowledge-based methods \cite{Sotiras2013}. Our method is categorized as an interpolation-based image registration method \cite{Sotiras2013}. The two most reputable interpolation-based techniques are B-splines~\cite{rueckert1999nonrigid} and thin-plate splines (TPS)~\cite{bookstein1989principal}. Revaud et al.\cite{revaud2015epicflow} proposed an interpolation-based method that uses an edge-aware geodesic distance to weight the sample points  when interpolating from a sparse to dense representation.

Recently, the interest for learning-based methods has increased and the results have been impressive~\cite{balakrishnan2019voxelmorph,dalca2019unsupervised}. By training a neural network to predict the solution, the method no longer needs to solve an optimization problem at every time step which drastically reduces the execution time. However, those methods require training data, which directly determines their applicability in the sense that the trained model is specific to the modality of the training data.
\section{Method}
In the following sections we describe the different steps of the method. We start with the initial location of the trackers, then continue with specific tracker techniques and finally the sparse-to-dense interpolation scheme and how this network is trained. The workflow we suggest is shown in Fig.~\ref{fig:workflow}.
\begin{figure}[ht]
    \centering
    \fbox{
    \begin{tikzpicture}
    \node[text width=3cm] at (-2,0.5) {Images:};
    \node[text width=3cm] (track) at (-2,-1.5) {Tracking:};
    \node[text width=3cm] at (-2,-3) {Sparse-to-dense \\ interpolation:};
    \node[text width=3cm] at (-2,-5.0) {Displacement field};
    
    \node[inner sep=0pt,label={$t=0$}] (img0) at (0,0.5) {\includegraphics[scale=0.15]{Images/workflow/img_0.eps}};
    \node[inner sep=0pt,label=below:{Initialization}] (dots0) at (0,-1.5) {\includegraphics[scale=0.15]{Images/workflow/dots_0.eps}};
    \draw[->,thick] (img0.south) -- (dots0.north);
   
    \foreach 0 [count=\i] in {0,1}
    {
        \node[inner sep=0pt,label={$t=\i$}] (img\i) at (\i * 2,0.5) {\includegraphics[scale=0.15]{Images/workflow/img_\i.eps}};
        \node[inner sep=0pt] (dots\i) at (\i * 2,-1.5){\includegraphics[scale=0.15]{Images/workflow/dots_\i.eps}};
        \node[inner sep=0pt] (sparse\i) at (\i * 2,-3.0) {\includegraphics[scale=0.15]{Images/workflow/sparse_flow_\i.eps}};
        \node[inner sep=0pt] (dense\i) at (\i * 2,-5) {\includegraphics[scale=0.15]{Images/workflow/dense_flow_\i.eps}};
        \draw[->,thick] (img\i.south) -- (dots\i.north);
        \draw[->,thick] (dots0.east) -- (dots\i.west);
        \draw[->,thick] (dots\i.south) -- (sparse\i.north);
        \draw[->,thick] (sparse\i.south) -- (dense\i.north);
    }
    
    \node[text width=0.5cm](end) at (\i * 2 + 1.5,0.5) {$\hdots$};
    
    \node[draw,inner sep=9pt, rounded corners=0.5cm,fit=(track)(dots0)(sparse\i), label=right:\rotatebox{-90}{\textbf{Suggested method}}](box){};
   
    \end{tikzpicture}
    }
    \caption{From the first image ($t=0$), a fixed number of trackers are initialized at specific points. When a new image ($t=1,2 \dots)$ arrives the tracking positions are updated and a sparse displacement field is estimated. A sparse-to-dense interpolation scheme is used to get the desired dense representation.}
    \label{fig:workflow}
\end{figure}

\subsection{Selection of tracking points}
One factor that may have a strong influence on the performance of our method is how the set of tracking points is selected. In computer vision, there has been a range of work on automatically generating candidates such as edges and corners \cite{schmid2000evaluation,shi1994good,stephens1988combined}. On the other hand, it is also possible to use task-based prior knowledge, such as a segmentation mask, when selecting candidates. If the aim is to estimate the deformation of a tumor or organs at risk, a task-based candidate might be preferable considering the fact that the interpolation accuracy tends to deteriorate with the distance~\cite{Sotiras2013}. Therefore, we suggest to sample candidates inversely proportional to the distance from a segmentation mask.
\subsection{Tracking}
In an image registration problem, the motion vector from one image to the another is referred to as the displacement field. In conventional image registration methods, the problem is most often formulated as a regression problem where the goal is to find a displacement field that minimizes some metric. A common choice is to use static points in space and estimate the displacement field by the corresponding points in the next image. These points are often referred to as control points. We, on the other hand, define the displacement field from a sparse set of tracking points where the spatial location of these points is constantly updated. In the next section, we describe how to estimate a dense displacement field from the location of the tracking points.

We suggest discriminative correlation filters, (DCF) \cite{danelljan2018learning} as trackers which have lately proven extremely useful for visual tracking, achieving state-of-the-art results on several benchmarks \cite{bhat2019learning,bolme2010visual,danelljan2016beyond}. The idea is to train a correlation filter $f$ that distinguishes the tracking target from the background. A continuous score variable $y$ is associated with each image patch $x$ where a low score, $y \approx 0$, represents background while high score, $y \approx 1$, represents the target. The typical choice of the score is a sampled Gaussian function with its mean centered in the target position. The filter is initially trained on a set of image patches, $\{x_i\}_{i=1, \ldots,N}$, wherein the target and the desired output $y_i$ are defined. In practice it is convenient to generate data from the first frame by applying different augmentation techniques, such as rotating, blurring and flipping, etc. In the simple case where the image patch has a single channel the filter is determined by solving a regularized least-squares problem of the form~\cite{henriques2012exploiting}
 \begin{align}
    \label{eq_l2_mosse}
    f^\star = \arg \min_{f} \sum_{i} ||x_i \star f - y_i||^2 + \lambda||f||^2,
\end{align}
but more general versions, including both multi-channel and multi-resolution representations, exist~\cite{danelljan2018learning}.

Once trained, the filter can be used to predict scores by convolving it with subsequent images. For speed, the convolution can be restricted to a smaller search region. For each tracking point, the displacement vector is extracted as the target location corresponding to the highest score. One advantage with DCF is that the filter can easily be updated with new training data (online) without retraining the entire model. 

In this paper, we use the ECO tracker \cite{danelljan2017eco} which uses deep convolutional features together with a dimensionality-reduction step that makes it faster and more robust against overfitting.

\subsection{Sparse-to-dense interpolation}
We suggest a sparse-to-dense interpolation scheme based on normalized convolutions~\cite{knutsson1993normalized}. Traditional convolutions are well-suited for regularly sampled data, but less so if the data is irregularly sampled. Normalized convolution uses a confidence map derived from the data where a value of $0$ indicates no data is given at that point and values close to $1$ indicate a high certainty of the value. A dense signal can then be estimated from a sparse representation based on a non-negative applicability function and some basis functions. Estimating the applicability functions for several normalized convolutional filters using constant basis functions has shown good results~\cite{eldesokey2019confidence,hua2018normalized} when interpolating dense depth scenes from LiDAR data. 

We design our interpolation scheme as a neural network inspired by U-net~\cite{ronneberger2015u} with normalized convolutional layers. A detailed illustration of the network is shown in Fig.~\ref{fig:Nconv}. We propagate the sparse displacement field and the confidence map through the network by using non-negative constraints on the weights and as a loss function we use a combination of Huber loss, to minimize the estimation error, and a term that maximizes the output confidence~\cite{eldesokey2019confidence}. 

We generate synthetic displacement fields to train our model. We assume that the displacement field is diffeomorphic and use a geodesic shooting method~\cite{miller2006geodesic} to  generate synthetic displacement fields independently of any imaging data. We artificially create a sparse representation by uniformly sampling coordinates in the generated displacement field. We use a binary confidence map based on the sampled coordinates. To avoid overfitting, we generate unique displacement fields for every iteration. The model is trained offline and only executed when interpolating the sparse representation from the trackers to the dense displacement field. 
\subfile{nconv}
\section{Experiments}
For evaluation we used a publicly available cardiac MRI dataset\footnote{\url{http://www.cse.yorku.ca/~mridataset/}}. The dataset contains measurements from 32 subjects (one was removed due to poor image quality). The data of each subject contains an image sequence of 20 images where the ventricle starts in an expanded position, after a few time steps the ventricle has contract an then expands again to its starting position. Each time step has $8$ to $15$ slices of $256 \times 256$ pixels and the slices include manual segmentation of the epicardium and endocardium. We used the manual segmentation masks to segment the myocardium (region in between the epicardium and endocardium) and the ventricle (region inside the endocardium). Since our method is 2D based we select a specific slice in the dataset. We compared our method with several others by applying the estimated displacement fields on the segmentation and calculating the Dice score~\cite{dice1945measures}, 
\begin{align}
    \text{Dice} = \dfrac{2 |S_{t} \cap \hat{S}_{t_0 \to t} |}{|S_t| + |\hat{S}_{t_0 \to t}|}
\end{align}
where $S_t$ is the manual segmentation at time $t$ and $\hat{S}_{t_0 \to t}$ is the warped segmentation using the image at time step $t_0$ and the estimated displacement field between the image at time step $t_0$ and the image at time step $t$. If this paper is accepted the implementation will be publicly available here\footnote{\ifthenelse
{\boolean{anonymous}}
{***}{\url{https://github.com/}}}.
\subsection{Setup}
We initialized 150 different trackers by sampling points inversely proportional to the distance of the union of the two segmentation masks. For each tracker we used an image patch of size $10 \times 10$ pixels and a search area of $4.5$ times the patch size. The trackers were trained on the first image and the filter coefficients are updated online after every new image. For the sparse-to-dense interpolation scheme we used a binary confidence map with $1$ for tracking coordinates and $0$ elsewhere.
\subsection{Result}
For a specific image plane ($z = 4$) the Dice scores between the warped and ground truth segmentation masks were calculated. For comparison we used the following methods i) the same tracker (ECO) but replacing the spare-to-dense interpolation scheme with thin-plate splines (TPS)~\cite{bookstein1989principal}, ii) B-spline image registration~\cite{klein2007evaluation} with mean square error using a grid size of $20 \times 20$ pixels and iii) the Demons algorithm~\cite{thirion1998image}. Table~\ref{tab_result} shows the mean and standard deviation of the result for all image frames and subjects.
\begin{table*}
\centering
\caption{The result for all images and subjects. The Dice mean ($\mu$) and standard deviation ($\sigma$) is shown for each method and segmentation regions (myocardium and ventricle).}
\label{tab_result}
\ra{1.3}
\begin{tabular}{@{}lrrrrcrrrcrrr@{}}\toprule
& \phantom{abc} & \multicolumn{2}{c}{Myocardium (Dice)} & \phantom{abc} && \multicolumn{2}{c}{Ventricle (Dice)} &
\phantom{abc} \\
\cmidrule{3-4} \cmidrule{6-8}
& & $\mu$ & $\sigma$ &&& $\mu$ & $\sigma$\\ \midrule
\textbf{ECO + NCONV} & & $0.810$ & $0.066$ &&& $0.870$ & $0.086$\\

ECO + TPS\cite{bookstein1989principal} & & $0.791$ & $0.062$ &&& $0.866$ & $0.081$\\

B-Spline\cite{klein2007evaluation} & & $0.730$ & $0.148$ &&& $0.837$ & $0.125$\\

Demons\cite{thirion1998image} & & $0.760$ & $0.111$ &&& $0.853$ & $0.111$\\
\bottomrule
\end{tabular}
\end{table*}

In Fig.~\ref{fig:average} we show the mean Dice score at each time step, averaged across all subjects. In Fig.~\ref{fig:compare} we show the warped segmentation for each method between two time steps where the ventricle first is in its expanded position and then in its contracted position. This is illustrated by time steps $t=0$ and $t=8$ for subject $15$.
\begin{figure}
    \centering
    \includegraphics[width=.75\textwidth]{Images/average.eps}
    \caption{The figure shows the mean Dice score at each time step, averaged across all subjects for both the myocardium (left) and ventricle (right) segmented areas.}
    \label{fig:average}
\end{figure}
\begin{figure}
\centering
\begin{subfigure}{.20\textwidth}
  \centering
  \includegraphics[width=1.0\linewidth]{Images/mov.eps}  
  \caption{}
  \label{fig:sub-mov}
\end{subfigure}
\begin{subfigure}{.20\textwidth}
  \centering
  \includegraphics[width=1.0\linewidth]{Images/ref.eps}  
  \caption{}
  \label{fig:sub-fixed}
\end{subfigure}
\\
\begin{subfigure}{.19\textwidth}
  \centering
  \includegraphics[width=1.0\linewidth]{Images/manual.eps}  
  \caption{}
  \label{fig:sub-manual}
\end{subfigure}
\begin{subfigure}{.19\textwidth}
  \centering
  \includegraphics[width=1.0\linewidth]{Images/eco_nconv.eps}  
  \caption{}
  \label{fig:sub-eco_nconv}
\end{subfigure}
\begin{subfigure}{.19\textwidth}
  \centering
  \includegraphics[width=1.0\linewidth]{Images/eco_tps.eps}  
  \caption{}
  \label{fig:sub-eco_tps}
\end{subfigure}
\begin{subfigure}{.19\textwidth}
  \centering
  \includegraphics[width=1.0\linewidth]{Images/bspline.eps}  
  \caption{}
  \label{fig:sub-bspline}
\end{subfigure}
\begin{subfigure}{.19\textwidth}
  \centering
  \includegraphics[width=1.0\linewidth]{Images/demon.eps}  
  \caption{}
  \label{fig:sub-demons}
\end{subfigure}
\caption{Illustration of the segmented areas. The ventricle (inside the green contour) and myocardium (between green and blue contour) segmentation areas are warped using the estimated displacement field for the different methods. Fig.~\subref{fig:sub-mov}~-~\subref{fig:sub-fixed} show the image at two time steps  ($t=0$ and $t=8$) for subject $15$. Fig.~\subref{fig:sub-manual}~-~\subref{fig:sub-demons} show the zoomed (yellow rectangle) segmented areas for the manual followed by the warped using the estimated displacement fields from time step $t=0$ to $t=8$.}
\label{fig:compare}
\end{figure}
\section{Discussion and conclusion}
Our method shows an overall better performance in terms of Dice score when compared to conventional methods. As shown in Fig.~\ref{fig:average} and Fig.~\ref{fig:compare} our method handles larger movements better than conventional methods and it appears more clearly for smaller deformable region such as the myocardium. The method allows a trade-off between accuracy and computational capacity because of the flexibility in choosing the number and position of the tracking points. Although our current implementation is not fast enough for real-time estimations, we remain hopeful. A single tracker can be executed in real-time \cite{danelljan2017eco} but in the current implementation all trackers run sequentially. Since the trackers are independent of each other and thus the potential for significant speed improvements through parallellization exists. We also expect further gains from relating the tracking score to the confidence, and from tweaking the data augmentation for medical applications. 
\section*{Acknowledgement}
This research was partially supported by the 
\ifthenelse
{\boolean{anonymous}}
{***}
{\emph{Wallenberg AI, Autonomous Systems and Software Program (WASP)} funded by Knut and Alice Wallenberg Foundation}. The authors would like to thank \ifthenelse
{\boolean{anonymous}}
{***}{Dr. David Tilly} for helpful comments.

