The 2014 liver ultrasound tracking benchmark

Abstract The Challenge on Liver Ultrasound Tracking (CLUST) was held in conjunction with the MICCAI 2014 conference to enable direct comparison of tracking methods for this application. This paper reports the outcome of this challenge, including setup, methods, results and experiences. The database included 54 2D and 3D sequences of the liver of healthy volunteers and tumor patients under free breathing. Participants had to provide the tracking results of 90% of the data (test set) for pre-defined point-landmarks (healthy volunteers) or for tumor segmentations (patient data). In this paper we compare the best six methods which participated in the challenge. Quantitative evaluation was performed by the organizers with respect to manual annotations. Results of all methods showed a mean tracking error ranging between 1.4 mm and 2.1 mm for 2D points, and between 2.6 mm and 4.6 mm for 3D points. Fusing all automatic results by considering the median tracking results, improved the mean error to 1.2 mm (2D) and 2.5 mm (3D). For all methods, the performance is still not comparable to human inter-rater variability, with a mean tracking error of 0.5–0.6 mm (2D) and 1.2–1.8 mm (3D). The segmentation task was fulfilled only by one participant, resulting in a Dice coefficient ranging from 76.7% to 92.3%. The CLUST database continues to be available and the online leader-board will be updated as an ongoing challenge.

six methods which participated in the challenge. Quantitative evaluation was performed by the organizers with respect to manual annotations. Results of all methods showed a mean tracking error ranging between 1.4 mm and 2.1 mm for 2D points, and between 2.6 mm and 4.6 mm for 3D points. Fusing all automatic results by considering the median tracking results, improved the mean error to 1.2 mm (2D) and 2.5 mm (3D). For all methods, the performance is still not comparable to human inter-rater variability, with a mean tracking error of 0.5-0.6 mm (2D) and 1.2-1.8 mm (3D). The segmentation task was fulfilled only by one participant, resulting in a Dice coefficient ranging from 76.7% to 92.3%. The CLUST database continues to be available and the online leader-board will be updated as an ongoing challenge.
Keywords: respiratory motion, ultrasound, tracking, image registration, image guidance, motion estimation, challenge (Some figures may appear in colour only in the online journal)

Introduction
Ultrasound (US) imaging is a widely used medical imaging technique. As US has high temporal resolution and is non-ionizing, it is an appealing choice for applications which require tracking and tissue motion analysis, such as motion compensation in image-guided intervention and therapy. Conformal and minimally invasive tumor treatments, such as high intensity focused ultrasound and intensity-modulated radiation and proton therapy, deliver highly localized dose into the target tissues. Yet, the motion of the organs in the treatment region is a critical limitation. Specifically, we want to address the issue of respiratory motion in the liver (Keall et al 2006, Shirato et al 2007. Note that liver tumors are not necessarily visible in US images. Instead, other visible structures (e.g. vessels) are tracked and these tracking results are used as input surrogate data to 4D motion models to predict the tumor position , McClelland et al 2013.
Despite the rapid development of image-guided therapy, intervention systems and medical imaging tools, the translation into clinical practice of automated motion estimation is limited. One of the main reasons for algorithms not being integrated in clinical practice is the lack of adequate validation. Open datasets for designing and testing tracking algorithms are missing, and private datasets differ in size, image dimension and sequence length. The variations in tracking objective (full organ, anatomical landmarks, tumor) and validation strategies are additional impediments to strategy comparisons. For image-guided therapies, tracking methods should have high accuracy, robustness over the duration of the therapy and real-time capability.
Several methods have been proposed for tracking human liver structures on US sequences. Yet quantitative evaluation of tracking the human liver under free breathing was reported only by Banerjee et al (2015), Harris et al (2010), Lediju Bell et al (2012) and Vijayan et al (2014) for 3D sequences and by Cifor et al (2012), (2013), De Luca et al (2012 and De Luca et al (2013) for 2D sequences. Intensity-based and hybrid approaches achieved good accuracy (∼1.4 mm mean tracking error  and ∼90% mean overlap ratio (Cifor et al 2013)). The non-rigid registration method of Vijayan et al (2014) estimated liver motion with an error of 1 mm (75% percentile of a root-mean squared metric over all datasets), which was lower than the inter-observer variability of 1.4 mm. More recently, a fast 3D affine block-matching algorithm with an outlier rejection strategy achieved a mean tracking error of 1.8 mm (Banerjee et al 2015). However, these methods were only tested off-line on short sequences (<1 min). For longer sequences (5-10 min long), a learning-based block-matching algorithm (De Luca et al 2013) achieved a mean tracking accuracy of ± 1.0 0.6 mm. In this paper we present the outcome of the open challenge for the validation of liver motion tracking algorithms. The reported methods are selected from the ones presented at the open challenge on liver US tracking (CLUST, http://clust14.ethz.ch), held in conjunction with the 2014 international conference on medical image computing and computer assisted intervention (MICCAI 2014). The aim of CLUST was to present the current state-of-the-art in automated tracking of anatomical landmarks in the liver (vessel centers (2D), vessel bifurcations (3D) and tumor contours (2D)) and enable comparison between different methods. For the test set, the annotations of the first images were provided, which needed to be tracked over time. This paper reports the results for the full test set, while the challenge proceedings exclude 20% of the test data, which were distributed shortly before the MICCAI conference. Furthermore this publication reports results after imposing restrictions against adjusting parameters per sequence by visual inspection (which is not realistic). Thus, method parameters were either automatically determined or generally fixed for at most each US scanner and task. In addition, we investigate the inter-observer variability of the annotations, analyze various aspects influencing the tracking performance, and explore the tracking performance when fusing results.
The paper is organized as follows. In section 2 we describe the challenge data and tracking objectives. In section 3 the methods proposed by the 6 participant groups are presented. The fusion of all methods, by considering the median of their results, is also considered. The evaluation criteria are described in section 4 and the tracking results are reported and compared in section 5. Discussion and conclusions of the challenge outcome are provided in sections 6 and 7.

Ultrasound data
The collected database included a total of 54 US sequences of the liver of patients and volunteers under free breathing. The data were provided by 6 groups, namely the Computer Vision Laboratory, ETH Zurich, Switzerland (ETH) (De Luca et al 2013, Preiswerk et al 2014; mediri GmbH, Heidelberg, Germany (MED); Institute of Biomedical Engineering, University of Oxford, UK (OX) (Cifor et al 2013); Biomedical Imaging Group, Departments of Radiology and Medical Informatics, Erasmus MC, Rotterdam, The Netherlands (EMC) (Banerjee et al 2014); Joint Department of Physics, Institute of Cancer Research & Royal Marsden NHS Foundation Trust, London and Sutton, UK (ICR) (Lediju et al 2010, Lediju Bell et al 2012; and SINTEF Medical Technology, Image Guided Therapy, Trondheim, Norway (SMT) (Vijayan et al 2013).
The sequences were acquired with 6 US scanners, 7 types of transducer and different acquisition settings. An overview of the data is given in table 1. The length of the sequences ranges from 4 s to 10 min, with a temporal resolution in the range of 6-25 Hz. The dataset is divided into three subsets, according to the image dimension and annotation type. The first subset is composed of 28 2D sequences from healthy volunteers with point-landmark annotations. The second subset contains 10 2D sequences from 5 patients with segmentation annotations. The third subset consists of 16 3D sequences with point-landmark annotations from healthy volunteers. The data were anonymized and divided into a training set (10% of the sequences) and a test set (90%). For the training set annotations were released, to allow for some tuning of the tracking 4V-D

2.5
Note: The test set is listed in black font. The training sequences, for which all available annotations were provided, are highlighted in red. algorithm. For the test set, the annotations of the first images were provided. These needed to be . For all sequences, annotations of the first frame (P j (0) or S j (0)) were provided.

Methods
The challenge raised interest from 55 individuals worldwide. All of them successfully downloaded the data. Only 8 from the downloaders submitted their results to CLUST14. After the withdrawal of one of the groups, a total of 7 papers were accepted to the MICCAI workshop. In this paper we included the 6 contributions with the highest mean accuracy, namely (in alphabetical order of the abbreviations) from Konica Minolta Inc., Osaka, Japan  In the following we briefly describe each algorithm. An overview of the methods' main features is presented in table 2. All 6 presented methods were tested on 2D landmark tracking. Of these, 3 were extended to 3D tracking, while only one submission covered all the challenge tasks. A more detailed description of each method can be found in the workshop proceedings 15 .

KM: template matching
Kondo, Konica Minolta Inc. (KM), developed a multiple-template matching method, based assuming that all regions of interest (ROIs) from the same image move along almost the same direction, and the motion has high periodicity due to breathing. Similar to De Luca et al (2013) and Matthews et al (2004), templates are selected from a plurality of recent frames to exploit the highly periodicity. The method consists of five steps: Step1: selection of global and long-term templates. From the first image I(0) two templates (i.e. unchanging subimages) are selected: a global template ( ) T x G , which is used to determine the motion on the entire frame and is defined by the largest rectangle which can be inscribed in the US image; and a long-term template , which is a squared ROI around the annotated point P j (0), whose size is based on minimizing the variance of the pixel values inside the ROI. The set of pixel coordinates included in and ( ) T x g are denoted as B L,j and B G , respectively. Step2: global motion estimation. For each image I(t), t > 0, the global displacement ( ) t d G , of point P j (0), is estimated by maximizing the normalized cross-correlation (NCC) over ( ) T x G using exhaustive search: Step3: long-term motion estimation. Similarly the long-term motion component is computed as where S L is the set of tested pixel displacements. It depends on the global motion estimation (Step2) as follows: are stored in the short-term buffer. Step4: short-term motion estimation. Firstly, the cycle length c is estimated from the past tracking results. Then, two short-term templates for the t-th frame are selected: the ROI with the maximum ; and the ROI with the minimum . Principal component analysis (PCA) is applied to the 2D trajectory of the tracking positions, and the motion estimation is performed only in the first principal direction. The resulting is determined by providing the maximum NCC for these templates ( where α = 0.95 is a weight that prioritizes the long-term motion estimation.

Run-time.
This method was implemented in C++ using OpenCV and OpenMP on a computer with Intel Core i7 3.3 GHz GPU, 6 cores and 64 GB memory. The average processing time was approximately 84 ms per frame, which can be improved by optimizing the motion estimation routine, currently based on an OpenCV function.

MEVIS: variational real-time registration
König et al Fraunhofer MEVIS Project Group Image Registration, Lübeck, Germany, (MEVIS) proposed a novel scheme for point tracking in long US sequences, based on an efficient, state-of-the-art variational non-linear registration method (Modersitzki 2009, König and. It is extended by a moving window scheme with additional fallback strategy to minimize tracking failures. (2015) 5571 fields (NGF) (Haber and Modersitzki 2006), and S is the second order curvature regularizer (Fischer and Modersitzki 2003). See ,  and Modersitzki (2009) for details on the numerical optimization.

Registration algorithm. Transformation T between images A and B was determined by a variational approach with the objective function
Tracking scheme. The non-linear registration algorithm is embedded in a tracking framework. By calculating registrations of moving windows on each frame, as the center of a new window, this process is then repeated for all frames.
As a safeguard, the original intensity at position P j (0) is compared to the current intensity Parametrization. Parameters for optimization, deformation resolution ( × 17 17) and multilevel scheme (2 levels) were kept constant. The window size was 50 mm in each dimension. Weight β is automatically determined before each registration as The parameters α and η (edge parameter of NGF) were manually chosen per scanner (ETH datasets: 2). The threshold was θ = 0.5 for ETH and MED2, and θ = 0.75 for MED1.

Run-time.
The algorithm achieved close to real-time performance in all cases (23-90 ms per frame), exceeding acquisition rate in 48% of all cases, computed on a three year old Intel i7-2600 PC with 3.40 GHz. Thus real-time performance is easily within reach when using recent hardware.

MEVIS + FOKUS: high performance optical flow
The method of Lübke, Fraunhofer MEVIS, Bremen, and Grozea, Fraunhofer FOKUS, Berlin, Germany, (MEVIS + FOKUS) is based on a two-frame motion estimation by an optical flow approach using polynomial expansion (Farnebäck 2003) and builds on the OpenCV function calcOpticalFlowFarneback. It yields a motion vector field for the entire frame which allows to track multiple points in parallel without any computational overhead. The method consists of approximating images locally with quadratic polynomials and then obtaining the dense displacement field by inferring the local displacement analytically. This is based on the coefficients of the fitted polynomial surfaces, with smoothing coming from the assumption of a global parametrized displacement model. In contrast to the original method, the OpenCV implementation did not provide a certainty for each pixel. This reduced to some extent the accuracy when the tracked point approached repeatedly the border of the acquisition region.
Real-time capability is achieved by downscaling the images before tracking and filtering, and subsequent upscaling of the result to the original resolution. A fixed downscaling factor to 30% has been chosen to match the computation time with the frame-rate. Fixed-parameter bilateral filtering and histogram equalization are used to obtain more stable results. Outliers, which are defined as motion vector component changes greater than 12 pixels, are discarded and the previous positions are used instead. Linear fitting (sliding window) is used to smooth the trajectory for each dimension to eliminate high-frequency motion.
This method was applied to 2D and 3D datasets. In the latter case it is in fact a 2.5D method, as we are evaluating the motion in two orthogonal slices intersecting at the manual annotation. This yields two separate tracking results with redundant information on the intersection. The method is sensitive to out-of-plane-motion as any 2D tracking method, since no adjustment based on the combined information has yet been done.
Another approach was proposed by this group and applied on 2D datasets, see the CLUST14 proceedings for details (Lübke and Grozea 2014). It was based on patch-matching, maximizing NCC, random sampling, explicit masking to improve accuracy at the border, and GPU implementation. This approach was not included in this paper, as it was on average slower (84%) and had a similar mean accuracy (3% better) as the group's included method.
Run-time. The method was implemented using the Mevislab software (MeVis Medical Solutions, Bremen, Germany), Python, OpenCV and Numpy. The run-time was measured on a Intel Core i7-4770k with 32 GB RAM. The average run-time for all 2D points is 40 ms. As the 3D tracking is performed on two orthogonal slices, the average processing time per volume of 61 ms is slightly higher than in the 2D datasets but still below the 3D images frame-rate (real-time), except for the ICR data subset with temporal resolution of 24 Hz.  Particle Filter. A particle filter based algorithm (Isard and Blake 1996, Arulampalam et al 2002, Zhang et al 2010 is used to track P j (0) in a 2D or 3D sequence. A small region around target P j (0) is considered to translate through the sequence, yielding a d = 2,3 dimensional state space for 2D and 3D translations respectively. The state of the target is represented by a probability density function in state space, approximated by a fixed size set of N particles. Each particle n holds the state ( ) = { ( )} t s t s n n d of a position hypothesis and the associated weight w n (t), normalized by ∑ ( ) = w t 1 n n , for = … n N 1, , . With each tracking step s n and w n are updated by estimations and observations, respectively.
Target Description. The target region is described by an under-sampled grid, characterized by radius R 0 around P j (0) and grid constant R 1 (see figure 3), resulting in N P grid points r i , = … i N 1, , P . Each r j has associated information for belonging to the bright p i brt and dark where b max and b min are the maximal and minimal intensity among all ∈ ( ) b I r 0, i i . Estimation. States are updated by re-sampling (Isard and Blake 1996) and application of stochastic drift and diffusion (equation (5)), taking into account the direct predecessor state (Markov property): When coinciding with bright image areas, bright grid points increase the weight while dark grid points decrease it. The final weighting function is with Heaviside function Θ(⋅) to set negative weights to zero, and with squaring of weights as this improved performance. The tracking result is computed to be the centroid of the target model, transformed into the image by the weighted mean sample. The tracking step was repeated M times to allow the distribution's mean to settle, leading to reduced lag.
Execution. Per target, R 0 and R 1 were adapted to fit the given target scenario and resolution of I (0) Bremen, Germany) for high level evaluation routines using Python scripts. The code was executed single threaded on a Windows 7 machine with an Intel Core i7-2600 CPU @ 3.4 GHz and 32 GB RAM.

PhR-sparse demons
Somphone et al Medisys Lab, Philips Research, Suresnes, France, (PhR) based their method on the Sparse Demons framework, where a dense displacement field is found by minimizing an energy E defined only on a small number P of points of interest where R and T are the reference and template images respectively, Ω is the image domain and δ is the Dirac function. The dissimilarity between the images is measured by D ( ∫ Minimizing E w.r.t. v is done by gradient descent. Calculus of variations results in the following evolution equation: The subsequent algorithm has similarities with the Demons algorithm (Thirion 1998, Mansi et al 2011. Its computational complexity is however lower since image forces are only computed at points x i .

2D point tracking. D( ) =
x x /2 2 is chosen as dissimilarity measure. The reference points of interest are chosen in the neighborhood of the landmarks, based on having a high amplitude of the image gradient. The tracking strategy consists of two phases. In the initial (t − 1)-to-t strategy for ∈ [ ] t 0, 100 , R = I(t − 1) and T = I(t), and a mean reference patch (of size × 30 30 ixels 2 ) is built around P j (0) from the patches centered on the tracked landmark positions. To prevent error accumulation, the 0-to-t patch tracking is used from t > 100.

3D point tracking.
Reference points x i are chosen in the neighborhood of the landmark in a × × 80 80 80 mm 3 region, regularly spaced by 10 mm. The baseline tracking follows the (t − 1)-to-t scheme. To prevent drifting, a 0-to-t tracking is enabled if the difference between the histograms of I(t) and I(0) is lower than 20% of the histogram difference between I(1) and I(0), or landmarks positions are closer than 1.8 mm to P j (0).

2D segmentation tracking.
The entropy of the difference between the reference and the transformed template is used as dissimilarity measure, since the OX sequences display large intensity changes between frames. The points of interest x i are selected in the neighborhood of the segmentation boundary, again based on the image gradient amplitude. As the OX sequences are short, only the (t − 1)-to-t scheme is used.

Run-time
The average processing time of 25 ms for 2D and 100 ms for 3D was obtained on a multithreaded PC platform. Yet, the method was not specifically optimized for run-time. Each pixel contribution to a histogram bin u is weighted based on the radially symmetric Epanechnikov kernel (Epanechnikov 1969), which assigns smaller weights to pixel locations farther away from the center.
In each frame I(t) the goal is to find ( ) p x that best matches q. The discrete Bhattacharyya coefficient (Comaniciu et al 2003) is used as similarity measure between q and ( ) p x . In each frame t the procedure to find the location x that maximizes ρ( ) x is started at location x 0 , which initially is set to the previous solution − x t 1 . After linearization through a Taylor series expansion around x 0 , ρ( ) x can be maximized using the mean shift procedure (Fukunaga and Hostetler 1975). When using the Epanechnikov kernel, the mean shift iteration step to move the kernel center position from x 0 to its new position x 1 is The weight w i for each pixel depends on a comparison of the histogram bins of q u and ( ) p x u that intensity ( ) I t x , i falls into: where δ is the Kronecker delta function and ( ( )) b I t x , i a binning function that maps intensity ( ) I t x , i to a histogram bin number. After each mean shift iteration, convergence is checked (number of iterations >20 or mean shift vector length <0.05 pixels) and if met, x t is set to x 1 , otherwise x 0 is set to x 1 and the mean shift procedure is repeated.
This method incorporates modifications proposed by Ning et al (2012) to make it adaptive to scale and orientation. First the target's scale is estimated based on the sum of the weights of all pixels in the search region and adjusted by the current Bhattacharyya coefficient. The estimated area is then used to adjust an ellipsoidal target descriptor to match the width, height and orientation of the ellipsoidal target descriptor, which is manually initialized in I(0). Between two adjacent frames the kernel size is enlarged by parameter Δd, which is set to 1 and 3 for all ETH and MED sequences respectively.
Assuming motion periodicity, two failure recovery strategies are employed. First, if the Bhattacharyya coefficient ρ drops below 0.8, the found target position is discarded and tracking is repeated using the target search area defined in I(0). Second, if the target search area in the current frame is three times larger than the initial target's size, the search area and its position is reset to the one in the first frame.
Run-time. The algorithm was implemented in MATLAB 2013b, and the experiments were conducted on a machine with an Intel i5-3320M processor at 2.6 GHz clock speed and 8 GB RAM. Tracking speed using this hardware set up was approximately 33 ms.

Tracking by decision fusion
Fusing the results from different methods or annotations has shown improvements for various applications (Sinha et al 2008) including classification (Kittler et al 1998), segmentation (Rohlfing et al 2003, Heckemann et al 2006, and tracking. Therefore we included an investigation of such a fusion approach. The tracking results of all six previously described methods were combined by computing for each frame t the median position of the tracked points P j (t) from the automatic methods. Using the median helps to reduce the influence of outlier results. Furthermore we determined the performance when fusing only the two methods which had the lowest mean tracking errors on the training data.

Evaluation
We compared the performance of the methods described in section 3 on the test data, consisting of 66 point-landmarks (vessel centers) in 26 2D sequences and 28 point-landmarks (vessel bifurcations) in 13 3D sequences, which the observers were confident to be able to reliably annotate. We evaluated the performance of the segmentation method in section 3.5 (PhR) on 11 manually segmented tumors in 9 2D sequences. In the following we describe the evaluation scheme used to validate and quantify the tracking accuracy.

Point-landmark tracking error
We randomly selected at least 10% of the images from each sequence and manually annotated the corresponding position of the initial point P j (0) in each of these selected images ( ) I t , denoted as ¯( ) P t j . The number of annotated frames per sequence is listed in table 1. For the annotated frame ( ) I t and landmark j, we calculated the tracking error (TE) as the Euclidean distance between the estimated landmark position ( ) P t j and its annotation (ground truth) ¯( ) P t j : V De Luca et al Phys. Med. Biol. 60 (2015) 5571 We summarized the results by the mean (MTE), standard deviation (STD) and 95th percentile of the single distribution including all ( ) t TEĵ belonging to a particular subgroup. These subgroups were the individual landmarks j (MTE j ), the image groups (MTE ETH , MTE MED1 and MTE MED2 for 2D sequences, and MTE EMC , MTE ICR and MTE SMT for 3D sequences) and the landmark dimensionality (MTE 2D and MTE 3D ).
We also included the motion magnitude of the landmarks, defined as: We estimated the inter-observer variability of the results. For this two additional experts annotated a randomly selected subset of the images marked by the first observer amounting to a total of 3% of all images. We then defined as ground truth the mean position over the three annotations and computed the tracking error as described above.
Per tracking task, the median results from all employed methods were pair-wise tested for statistically significantly differences at the probability level p = 0.001 using the sign test. The sign test was used as it is a non-parametric test, which neither assumes a normal distribution nor a symmetric distribution.

Segmentation accuracy
A clinical expert segmented the visible boundaries of the tumors corresponding to the reference segmentation in each frame t of the 2D OX sequences. The results of the segmentation tracking were quantified by computing the Dice coefficient: where S j and S j denote the manually segmented and the estimated tumor region, respectively, in each frame t. ∈ [ ] Dice 0, 100 % measures the overlap ratio of the two regions and is at best 100%. We summarized the results by the mean, standard deviation and 5th percentile of the Dice coefficient distribution for each segmented tumor in the entire sequence. For comparison, we calculated the initial Dice coefficient, before tracking, as follows:

2D landmarks
The results for the 2D point-landmark tracking are summarized in table 3. The MTE 2D ranges from 1.4 mm to 2.1 mm for the automatic methods, with best results achieved by MEVIS.
Fusing the results of all tracking methods improved accuracy by 15-41% in comparison to the individual results, achieving a MTE 2D of 1.2 mm. Three of the proposed methods had for some sequences a higher mean error (MTE j ) than the landmark motion M j .
To assess the robustness of all methods, we quantified the percentage of failures, i.e. the percentage of annotated frames per landmark for which > TE 3 mm or > TE 5 mm. These results, summarized in figure 6, show that the percentage of failures ranges between 6.3% and 15.8% (1.6% and 7.2%) for > TE 3 mm ( > TE 5 mm), with the fusion method being the best one.
We illustrate the difference in tracking performance between all methods for the point which had the highest variance of MTE (24.85 mm) across all methods, namely P 1 from MED-07 (sequence details are listed in table 1). From figure 4 it can be observed that high tracking errors occur for some of the proposed methods when the landmark position drifted after a deep inhalation (at frame t c ) of the subject under investigation. Comparing the associated ROIs, it seems that methods KM and PhR maintain a similar distance to the diaphragm, rather than stay with the vessel.
Lower errors were obtained with respect to the mean annotation of the 3 observers with MTE 2D reduced by 0.06 to 0.23 mm for all methods, while the mean motion M was similar, see table 5. Note that now the order changed between MEVIS + MED and TUM and between PhR and MEVIS + FOKUS. The MTE 2D of the 3 observers was more than 50% lower than any individual method. An overview of these results is given in the box-plots in figure 5.
There was low correlation between the motion magnitude of the landmarks and the tracking errors. In details, the sample Pearson correlation coefficients (ρ) between landmark motion and tracking errors for the individual methods ranged from 0.11 to 0.37, while for the observers there was no correlation (ρ ∈ [ ] 0.06, 0.11 ). We also studied if tracking errors were influenced by the change in image quality due to the range of center frequencies used during the US acquisitions (see table 1). Only low correlation (ρ < 0.43) was found between MTE ∈ j D 2 of the individual methods and center frequency of landmark j.
Tracking landmarks close to the acquisition border can be difficult. We analyzed the correlation between the landmark distance to the acquisition border and the tracking error for the landmark which comes closest to the acquisition border (P 2 from MED-08), and found a high correlation (ρ > 0.70) for all methods except TUM (ρ = 0.12).
The median results with respect to the mean annotation of the 3 observers from all methods were significantly different to each other, except for . The median tracking errors of the observers were similar, apart from Obs1's being significantly higher.

3D landmarks
Results for all methods on the 3D sequences are shown in table 4. On average, the most accurate results were achieved by PhR, with MTE 3D of 2.55 mm. The fusion method slightly improved the average results by 3%. The percentages of failures are shown in figure 8 for all methods. These ranged between 18.7% and 51.3% (6.8% and 29.6%) for > TE 3 mm (5 mm). The tracking errors with respect to the mean manual annotation of 3 observers are increased by 0.05-0.28 mm, while the mean motion is reduced by 0.36 mm. The error obtained by fusing  all results is slightly reduced, see table 5. An overview of these results is given in the box-plots in figure 7. Figure 9 shows the result with the highest MTE j for MEVIS + MED and PhR. Tracking methods and observers disagreed mainly along the vessel (x-coordinate), showing the challenge of tracking landmarks within elongated structures without a stable bifurcation.
All MTE 3D were larger than the estimated equivalent 3D error from MTE 2D 16 (e.g. +70% Obs1, 84% Fusion, 55% MEVIS + MED), while the estimated 3D motion was on average lower (−29%). Yet this clear relationship was not sustained for error measures calculated in pixels (see table 6). In particular Fusion resulted in very similar results.
The median error (see table 5) of MEVIS + FOKUS was statistically different (p < 0.001) with respect to the one of Fusion and PhR. No other significant difference between the methods existed. The median tracking errors of the 3 observers were statistically significantly lower than the ones from the 4 tracking methods, but not significantly different amongst each other.
We investigated if the landmark error is correlated with the motion and found little evidence for such a correlation, as the sample Pearson correlation coefficients ρ was in the range of 0.25-0.45 for all methods and observers, apart from MEVIS + FOKUS (ρ = 0.62).
For 3D, a low to moderate correlation (ρ ∈ [ ] 0.26, 0.51 ) was found between tracking errors (MTE ∈ j D 3 ) and center frequency of the US acquisition protocol (see table 1) for each landmark j.
We analyzed the correlation between the mean tracking error and the distance to the acquisition border for a 3D vessel bifurcation landmark, which was moving close to this border and had good local image contrast throughout the sequence. We found that the distance was moderately correlated with the TE of landmark P 1 from SMT-02 for all methods (ρ ∈ [ ] 0.47, 0.58 ).

2D segmentations
The results of the tumor segmentation task, performed by PhR are summarized in table 7. The segmentation accuracy, expressed in mean Dice coefficient, ranges between 76.7% and 92.3%. The mean accuracy is 8% higher than the initial overlap ratio, but lower in 3 sequences. Figure 10 shows the results of the segmentation on the sequence with the highest variance of the Dice coefficient, i.e. S 1 from OX-6. The overlap gradually decreases, likely due to lower image contrast compared to the other sequences. , which assumes equal error components.

Discussion
When considering the individual methods (described in sections 3.1-3.6), determining the local translations of the vessels by a Bayesian approach (MEVIS + MED) was computationally the most efficient algorithm and was runner-up (by 6%) in both landmark tasks. The Sparse Demons method on patches (PhR) performed best for the 3D landmarks, but had a relative poor performance for the 2D landmarks. For the case illustrated in figure 5, the method tracks the diaphragm as this is the main bright feature in the reference patch, resulting in high errors when the distance to the vessel changes. The variational registration approach within relative large regions (MEVIS) provided on average the most accurate results for the 2D landmarks, but was unfortunately not tested for the 3D landmarks. Registration to the initial frame and if needed to the previous frame, allowing smooth deformations, as well as including NGF in the image similarity likely contributed to this good performance.
Fusion of the tracking results from all the independent methods by taking their median value provided on average the highest landmark tracking accuracy, with improvements by at least 15% (3%) for 2D (3D) landmarks. This reduced mean motion by 82% (60%) resulting in a mean accuracy of 1.2 mm (2.5 mm). Yet it requires the run-time of the slowest method and n times the computing power, where n is the number of fused methods (in this paper n = 6 for 2D and n = 3 for 3D). To save resources, fusing only the 2 methods, which performed on average best on the training set, led to worse results for 2D (MTE = 1.6 2D mm for median of MEVIS + MED and MEVIS + FOKUS) as the training performance of MEVIS + FOKUS was not representative. More training data will hopefully enable selecting instead one of the better methods (e.g. MEVIS), which would have led to a slightly better fusion result (1.4 mm). Alternatively fusing more methods might circumvent this problem (e.g. incorporating MEVIS as third best from training data would result in an MTE of 1.3 mm). For 3D, similar results were achieved, with the median of MEVIS + MED and PhR lead to MTE = 2.5 3D mm. The tracking performance was generally not dependent on the motion magnitude. Comparing the errors in millimeters, tracking the 3D landmarks appears to be a harder task than tracking the 2D landmarks. This can mostly be explained by the lower image resolution, since 2D and 3D results from Fusion are very similar when considering the pixel resolution and assuming equal errors in each dimension (2D: 2.7 pixels, estimated 3D: 3.3 pixels, 3D: 3.4 pixels). Hence it seems very likely that advances in 3D image resolution will improve tracking Note: The results are in millimeters and ranked according to alphabetical order of the methods. results. Other factors include some worse image quality (e.g. EMC), a more cumbersome annotation task and greater difficulties of visually inspecting the 3D images for parameter optimization.
Methods that adjusted their parameters per sequence for the MICCAI 2014 challenge were TUM, MEVIS, MEVIS + MED and MEVIS + FOKUS (only 2D case). The mean performance of these 'tuned' methods for all 2D test data was worse (TUM 1.84 mm, MEVIS 1.51 mm), similar (MEVIS + MED 1.52 mm), or better (MEVIS + FOKUS 1.91 mm) than with fixed parameters (table 3). For the 3D case, resubmission was needed only from MEVIS + MED, for which the original submission resulted in worse accuracy (MTE = 2.80 3D mm) than fixed parameters (table 4).
We compared the methods to previously published results (Cifor et al 2013, De Luca et al 2013, Vijayan et al 2014 by re-computing the tracking error for the same subset (excluding any training data), the same error measure and annotations. For 2D landmark tracking, two of the individual methods performed on average better than previously published results, namely MEVIS + MED achieved an MTE (95th TE) of 0.6 (1.4) mm and MEVIS of 0.7 (1.6) mm, versus 0.8 (1.7) mm from De Luca et al (2013). This comparison was based on ETH-01 to ETH-09 without ETH-05. For 3D landmark tracking, comparison was done on the SMT subset. The method proposed in Vijayan et al (2014) achieved an MTE (95th TE) of 3.6 (14.9) mm. PhR and MEVIS + MED obtained on average better results, with an MTE (95th TE) of 2.4 (6.8) mm and 2.6 (7.3) mm respectively. The single method attempting the tracking of tumor regions performed similar to the previous baseline method (Cifor et al 2013). The comparison was possible only for 3 sequences, where PhR was on average 5% (OX-01) and 1% (OX-04) worse, but 4% better (OX-02) than the method from Cifor et al (2013).  (see table 4). TE is evaluated with respect to one observer.
On average, the tracking performance was worse than the observer annotation accuracy, which hints at potentials for further method improvements. Error reduction beyond the observer accuracy will require improvements in image resolution and quality. The clinically acceptable tracking error depends on the application. The American Association of Physicists in Medicine recommends for external-beam radiation therapy that respiratory management techniques are warranted when the target motion is greater than 5 mm (Keall et al 2006). This implies that techniques for predicting the target motion should have a maximum error of at most 5 mm. Tracking will need to be more accurate as errors are also introduced by temporal prediction to compensate system latencies and spatial prediction if the target cannot be directly tracked. Overall this was clearly not yet achieved by any automatic tracking methods (see figures 4, 6, 7 and 8) Figure 9. Illustration of tracking performance for landmark P 1 from sequence EMC-03. Tracking errors MTE ∈ − 1 EMC 03 with respect to the mean of 3 observers were 7.77 mm (MEVIS + FOKUS), 5.63 mm (MEVIS + MED) and 9.93 mm (PhR). The mean motion for the landmark was 11.71 mm. Inter-observer errors were 4.06 mm (Obs1), 11.59 mm (Obs2) and 9.76 mm (Obs3). The tracking results and annotations are shown at time * t for the same ROI( * t ) with planes cut at the corresponding ( *) P t 1 from each method. and also observers had some difficulties. Only for one method (MEVIS + MED) and the ICR sequences the tracking errors stayed always within the 5 mm limit (maximum TE: 4.7 mm), likely due to good image quality, smaller field of view, no disappearing structures, and high temporal and spatial resolution.  The high 95th percentile errors and number of tracking failures indicate that the methods are not very robust. Visual inspection of failure cases hint at problems due to accumulation of errors, disappearing structures at the image acquisition border, out-of-plane motion, substantial changes in image appearance or insufficient motion capture ranges. These might be the reasons of the high variance in the individual tracking performances across the sequences.
Tracking run-times per frame range from 1.8 ms to 84 ms for 2D and from 51 ms to about 100 ms for 3D. Processing times generally varied depending on the number of landmarks and the image appearance. Real-time performance was achieved by 4 methods for 2D and 3 for 3D landmarks (MEVIS + MED, MEVIS + FOKUS, PhR and TUM), assuming average acquisition times of 50 ms (20 Hz) and 125 ms (8 Hz). The remaining methods are not far off, indicating that further code optimization and hardware improvements are likely to provide real-time performance for the current dataset. Furthermore, off-line applications, such as dose accumulation during radiation therapy, do not require real-time performance. Hence, the fusion of several tracking results can be a promising approach, despite having the cumulative computational burden of all fused methods (which is however easily parallelizable), as it achieved the highest accuracy.
For the purpose of evaluating different tracking approaches, specific anatomical landmarks in the first frame were provided in CLUST. These landmarks were manually selected after inspecting the US images, to ensure that they would not disappear during freebreathing. In clinics, real-time applications should require as little manual interventions as possible. Therefore, tracking algorithms might benefit from incorporating the automatic detection of stable features to track throughout the treatment and update tracking regions when necessary. Feature detection might be based on analyzing the tracking results from an initial set of images covering a few breathing cycles. Testing this functionality is outside the CLUST challenge, as the aim is to compare tracking performances for the same landmarks.

Conclusion
This paper describes the results of the MICCAI 2014 challenge on liver US tracking (CLUST14), which enabled for the first time the quantitative, direct comparison of tracking methods for this application.
The challenge data included a large number of realistic sequences, which varied in length, spatial and temporal resolution, acquisition settings and US scanner.
Quantitative evaluation of all results showed a mean tracking error from 1.4 mm to 2.1 mm for 2D points, and from 2.6 mm to 4.6 mm for 3D points. Considering the median tracking results of all methods improved the mean error to 1.2 mm (2D) and 2.5 mm (3D). The segmentation task, fulfilled only by one participant, resulted in a mean Dice coefficient of 84.1%. All best approaches are comparable or better than the state-of-the-art.
Applicability for therapy guidance still requires general improvement of the 3D landmark tracking accuracy as well as reduction of tracking failures ( > TE 5 mm). Advances in image resolution and quality will support this task. Furthermore the diaphragm is a prominent feature in liver US images which should also be tracked.
For some sequences, the variability of the observers was particularly high, due to the difficulty in manually annotating 3D volumes. Having more observers might lead to a more reliable reference for computing the tracking results.
Most of the methods achieved, or are close to, real-time, while running on standard machines.
In conclusion, this CLUST14 challenge provided a good basis for a first comparison of US tracking methods for the liver. Its accompanying workshop facilitated lively discussions of the involved researchers. The research community can benefit from this benchmark and the CLUST challenge remains open for future participants to evaluate their method and be included in the online leader board. Future work includes the increase of the number of sequences and tracking of other available structures in the liver, e.g. the diaphragm.