4D MRI: Robust sorting of free breathing MRI slices for use in interventional settings

Purpose: We aim to develop a robust 4D MRI method for large FOVs enabling the extraction of irregular respiratory motion that is readily usable with all MRI machines and thus applicable to support a wide range of interventional settings. Method: We propose a 4D MRI reconstruction method to capture an arbitrary number of breathing states. It uses template updates in navigator slices and search regions for fast and robust vessel cross-section tracking. It captures FOVs of 255 mm x 320 mm x 228 mm at a spatial resolution of 1.82 mm x 1.82 mm x 4mm and temporal resolution of 200ms. To validate the method, a total of 38 4D MRIs of 13 healthy subjects were reconstructed. A quantitative evaluation of the reconstruction rate and speed of both the new and baseline method was performed. Additionally, a study with ten radiologists was conducted to assess the subjective reconstruction quality of both methods. Results: Our results indicate improved mean reconstruction rates compared to the baseline method (79.4\% vs. 45.5\%) and improved mean reconstruction times (24s vs. 73s) per subject. Interventional radiologists perceive the reconstruction quality of our method as higher compared to the baseline (262.5 points vs. 217.5 points, p=0.02). Conclusions: Template updates are an effective and efficient way to increase 4D MRI reconstruction rates and to achieve better reconstruction quality. Search regions reduce reconstruction time. These improvements increase the applicability of 4D MRI as base for seamless support of interventional image guidance in percutaneous interventions.


Introduction
During the last decade, 4D MRI has gained considerable interest in research, because it promises access to information on the respiratory motion of the thorax and abdomen free of radiation. Respiratory motion information is vital for many medical applications in diagnostic [1], treatment planning [2] and execution [3]. Our application scenario are MRI guided percutaneous interventions on the liver like radio frequency-, microwaveand cryoablation, biopsies or brachytherapy where the challenge of a moving target exist. 4D MRI methods have been proposed but none satisfy all needs for our interventional application. These needs are first, physiological correctness of the 4D sequence, and second, robustness against out-of-plane motion. In this study, we propose a new 4D MRI reconstruction method. It utilizes retrospective sorting of dynamic 2D TRUFI MRI slices and is capable of imaging the whole liver during free breathing and capturing organ deformations caused by respiration. It reconstructs physiologically meaningful sequence of respiratory states by utilizing a dedicated navigator frame and copes with out-of-plain motion.

Related Work
To our knowledge, there exist two approaches to acquiring 4D MRI, each with its unique advantages and disadvantages. The first is to acquire 3D MRI sequences in real-time as done by Kim et al. [4] and Bled et al. [5]. The advantages of this approach are that it does not rely on gating and thus supports imaging events that do not occur repeatedly, i.e., events that are not periodic. The disadvantages of this approach are its low temporal and spatial resolution [6,7] and its relative small FOV, rendering it impossible to capture the respiratory motion of large organs like the liver.
The second approach is to reconstruct volumes for different organ states or breathing phases in retrospection by binning previously acquired data. Two main types of this approach exist. In the first type, the k-space data is sparsely sampled and binned before reconstructing a volume for a given organ state [8][9][10]. The strength of this type lies in capturing periodic organ state changes with a large FOV within a few minutes depending on the length of the motion cycle. Its weaknesses are its assumption of strictly periodic organ motion. It thus can only reconstruct an average motion cycle of the target organ, which is not ensured to be physiologically meaningful. Furthermore, this type introduces image artifacts [11,12] that could hinder motion estimation from the reconstructed 4D MRI.
The second type of the second approach reconstructs fast dynamic 2D sequences at all slice positions to cover the organ of interest. Then retrospective gating is applied to the resulting 2D images, binning them by different organ states, i.e., breathing states, and sorting them in their respective volumes. Its advantages are its applicability for non-periodic or quasi-periodic changes in the organ state and its high temporal and spatial resolution. Hence it is well-suited to capture motion variation, e.g., deep or shallow, abdominal or thoracic breaths within one session. It can work with a navigator or respiratory signal to ensure physiological correctness of reconstructed motion. A further advantage of the binning strategy is its availability because it is readily usable with all MRI machines and all 2D sequences. Its disadvantages are that it is more time-intensive than the k-space binning and that much of the acquired data is redundant. The latter, however, can advantageously be used to increase the SNR of the reconstructed 4D images.
For both types, the surrogate can be intrinsic, relying on image or k-space information, or extrinsic, relying on externally recorded signals, e.g., from using a breathing belt or form tracking markers that are placed on the abdomen of the subject.
Siebenthal et al. [13,14] utilize navigator slices as surrogate and vessel cross-section tracking as matching criterion. Cai et al. [15] use the body area. Lee et al. [16] use sagittal diaphragm profiles and reconstruct one breathing cycle. Tong et al. [17] propose a graph-based sorting where the weights are based on image information and semi-automatic assigned respiratory phase although, they are only able to reconstruct one best breathing cycle and not a variety of breathing cycles. Romaguera et al. [18] propose a graph-based approach using pseudo-navigators. A drawback of the graph-based navigator-less approach is that physiological correctness cannot be ensured even if temporal coherence is ensured,

Materials and methods
We decided to follow the retrospective sorting approach because, as set out in the related work section, it is the only one suited for capturing physiologically meaningful, non-periodic organ motion with high temporal and spatial resolution and large field of views. Its only disadvantage is the long acquisition time, which can be overcome as shown in this work. Specifically, we build upon the proposed method of von Siebenthal et al. [13,14]. First, we describe the general concept behind the baseline method. Second, we describe our improvements that overcome the named drawbacks.

MRI Acquisition
Our MR data was acquired on a MAGNETOM Skyra MRI scanner (Siemens Medical Solutions, Erlangen, Germany). All images were acquired with a TRUFI sequence (TR = 40.0 ms, TE = 1.5 ms, flip angle = 30 degree, in-plane resolution 1.82mm x 1.82mm, out of plane resolution 4 mm, FOV: 255 mm x 320 mm, matrix size: 140 x 176). Using this setup, we achieve acquisition times of 200 ms per slice. The acquisition setup was chosen to mimic an interventional setup as closely as possible. This specifically means high acquisition speed and just good enough contrast to detect the respiratory motion. No body coil was used, which makes this 4D MRI method compatible with a wide range of external surrogates, including those that need a free line of sight to the abdomen of the subject. This includes but is not limited to surrogates that are based on a scan of the abdomen surface, or the tracking of markers on the abdomen. This is important to make the gathered motion information available for a wide range of interventional scenarios where different surrogates may be used to track breathing. In total, 19 data sets of 13 healthy subjects were acquired. One subject was imaged three times, four subjects were imaged twice, and eight subjects were imaged once. If a subject was imaged multiple times, then each data set acquisition was performed on different days to include variations that might occur in between imaging sessions. Each data set consists of two reference sequences and several interleaved sequences, one for each slice position, capturing the volume of the liver.
A reference sequence is a dynamic 2D MRI sequence of so-called navigator frames. A schematic depiction can be found in Fig 1. The navigator frames picture an image plane, in which the respiratory motion is clearly visible. In our case, we used a slice in sagittal orientation that intersects the target organ -the liver -and shows vessel cross-sections, because their spatial distribution describes the breathing state well. This sequence is the reference for the 4D reconstruction. The reference contains natural succession of different breathing patterns, like shallow or deep, thoracic or abdominal breathing and is thus physiologically highly meaningful. One reference sequence was acquired at the beginning and one at the end of each session. A reference sequence comprises 513 images (time points) covering a time of 102 seconds (about 20 breathing cycles). The former are sorted into the 4D MRI sequences based on information extracted from the latter. Data slices and navigator slices were imaged alternatingly, facilitating the interleaved character of the sequence. The navigator slices are positioned exactly as in the reference sequence, rendering temporal reconstruction possible. The data slice sweeps over the target organ in 4 mm gaps during acquisition (see Fig 3), rendering spatial reconstruction possible. The total number of interleaved sequences ranges between 38 and 57 sequences, depending on the size of the subjects target organ. Thus, the total acquisition time ranged between 40 min and 80 min, excluding planning time. Planning took roughly 15 min per subject. The total acquisition time can be reduced by optimizing the acquisition scheme and giving in-time breathing instructions to the subject to use the acquisition time more efficiently. half when using matching criteria that do not depend on a navigator slice. Schematic depiction of an interleaved sequence. An interleaved sequence consists of navigator frames and data frames that were imaged alternatingly. It shows a different breathing curve than the navigator sequence but contains similar breathing patterns.

4D MRI Reconstruction
Our method uses the reference sequence as grounds for temporal reconstruction of a 4D MRI sequence showing a physiologically meaningful course of breathing states. The general scheme of the reconstruction process is depicted in Fig 4. For each time point in the reference sequence, i.e., for each frame, a volume is reconstructed. First, the breathing state of the frame is determined. Second, in each interleaved sequence, all data frames are found that match the breathing state, using a matching criterion, see Fig 5. Third, the found frames are averaged (binned) to one slice to improve the SNR (signal-to-noise ratio). Fourth, the averaged slice is inserted (sorted) into the volume at its designated position, which is unique and known for each interleaved sequence. Doing this for all reference frames results in a continuous 4D MRI sequence. The reconstructed FOV's range from 255 mm x 320 mm x 152 mm to 228 mm (140 x 176 x 38 to 57 voxels) depending on the size of the target organ. In the next section the matching criterion is described in detail.

Matching Criterion
A matching criterion is used to find all data slices showing the reference breathing state within an interleaved sequence. The respiratory state of a frame is determined by its enclosing navigator frames. Hence, the matching criterion acts on pairs of navigator frames that encase another frame (navigator or data frame); see brackets in Fig 5. It is based on the displacement of tracked vessels within the navigator frames. Assume a navigator frame n ti at time point t i in the reference sequence that shows a reference breathing state BS r . We want to find a data frame d tj with the same breathing state as n ti . To this end, the enclosing navigator frames of both d tj and n ti are used. The enclosing navigator frames of d tj are n tj−1 and n tj+1 and the enclosing frames of n ti are n ti−1 and n ti+1 . The vessel displacements from n tj−1 to n ti−1 and from n tj+1 to n ti+1 are calculated. When the summed displacements of all vessels for two pairs of navigator frames is under a certain threshold, then the two enclosed frames are assumed to be a match, i.e., to show the same breathing state. The threshold is the only parameter of the method. It determines the maximally allowed displacements between two frames to  Scheme of finding data slices that match specific breathing state. On the left hand the reference sequence is depicted. The red bracket represents the third breathing state. It is found in the interleaved sequence, depicted on the right, by comparing the enclosing navigator slices. be counted as a match.
The vessel tracking is realized via template matching using OpenCV [19] and its similarity measure TM CCOEFF NORMED. The templates are manually defined for each tracked vessel cross-section in the reference sequence. To this end, a user identifies well trackable vessels in one slice of the reference sequence prior to the 4D reconstruction, which takes only a few seconds. In our case well trackable means that the vessel will be visible in most navigator frames and has a high contrast to the surrounding tissue.

Template updates and search region
One of the challenges for the template matching is the out-of-plane motion of the vessel cross-sections in the navigator frames. In these cases, the searched-for regions are changing their appearance throughout breathing; hence, they are difficult to find using fixed templates.
To increase robustness against the out-of-plane motion, we propose to apply template updates within the reference sequence. In Fig 6, one can see how the appearance of a vessel cross-section can change during a breathing cycle. The method starts with the templates that were defined manually on reference frame n t0 . Then, for each following navigator frame n ti that was captured at time point t i , the templates get automatically updated, as follows: The positions of all tracked vessels in n ti are found with subpixel precision using the templates from time point t i−1 . Then a new set of templates is cut from n ti based on the position of the matched templates. The template position is updated with floating-point precision. The updates ensure that changes in the appearance of the tracked vessel are represented in the updated templates. The subpixel precision in the updates is needed to avoid drift during the update. Fig 6. Out-of-plane motion and template updates. The figure shows a series of navigator slices. The green rectangle denotes a typical ROI that was manually determine as template. In the red rectangles one can see how the vessel cross-section changes its appearance during the breathing cycle. For viewing purposes only, the images gradation curve was altered globally to enhance contras.
Another concern of the reconstruction approach is speed. In its original form, the method matches each template against each navigator frame, resulting in a substantial computational burden. We propose to speed up the vessel tracking by exploiting spatial coherence between temporally adjacent navigator frames. The underlying assumption is that the next searched-for match is in a small spatial neighborhood around the previously found match, which is the case due to fast and continuous acquisition. Therefore, we only use a small neighborhood around the last matched template position as search area.
More over we automatically detect breathing states that cannot be reconstructed entirely and use that knowledge to inform where (temporally and spatially) the 4D sequence is incomplete. This is important for the later application, because of the visual feedback that can be provided to the physician in real-time when the motion information is insufficient to fuse the planning data to the interventional data.

Evaluation
We compare our method with the baseline method of Siebenthal et al. by means of reconstruction rate and image quality. We define the reconstruction rate as the percentage of the number of slices in the volume that could be reconstructed by the method. Note that this does not account for false positives or false negatives because the ground truth for this is not available to us. We also investigate how the acquisition order of reference sequence and interleaved sequence influences the methods ability to find matching data frames. We evaluate the point of false positives indirectly using a qualitative assessment of both approaches. The image quality is assessed in a double-blind study with interventional radiologists.

Reconstruction rate
We compare the reconstruction rate of both methods for different parameterizations. This is possible because the baseline method uses the same parameters in its matching criterion. We tested the parameters shown in table 1. We tested the threshold, for the values 0.5, 1, and 2. Evaluating different thresholds from a quantitative point-of-view allows us to judge which method will be more suitable for different applications that In the axial and sigittal orientation one can see, that our method is capable to reconstruct smooth and continuous volumes. differ in the kind of trade-off between precision and coverage that is preferable in the application. With lower (stricter) thresholds, the coverage goes down and the precision increases. With higher thresholds, the coverage increases and the precision decreases. We tested two similarity measures from openCV, namely TM CCORR NORMED and TM CCOEFF NORMED and we tested the influence of the chosen reference sequence, ref. 1

Reconstruction quality
We conducted a double-blind study with ten interventional radiologists to compare the reconstruction quality of both methods. Participants were recruited from a General Radiology clinic. Participants' professional experience ranged from 4 months to 20 years (median: 16 months). Each radiologist was shown a set of 48 slice pairs, where a pair contained a reconstructed slice from each method. The radiologists had to decide which of the images in a pair shows the anatomy of the target organ more faithfully, i.e., with fewer image artifacts. Participants did not see the two slices from each pair simultaneously but could switch back and forth between them as often as they wanted before picking one. Participants were asked to select the slice they considered better. A neutral option was provided. The slices were sampled from the reconstructed 4D MRI of both methods , respectively. For the evaluation, the parameter set chosen to be 1 px threshold and TM CCOEFF NORMED as a similarity measure for both methods. Only those volumes were considered for comparison for which both methods had a reconstruction rate of at least 80%. From these, 48 volume pairs were chosen randomly. Each radiologist was shown a different set.
Furthermore, in both volumes, we automatically masked slices out, where either of the methods did not find a matching data frame. We made both volumes identical in the amount and distribution of black slices. This was done because it is likely that a reduced reconstruction rate for a volume would be detrimental to its perceived reconstruction quality. Each slice pair was sampled at a random orientation and position that was chosen within a range, such that the sampled slice would show the target organ. Slices were sampled either in sagittal, coronal or axial orientation. Due to a software error, the number of slices for different planes was slightly imbalanced: Overall, 100 slices were shown for the sagittal and axial orientation each, and 280 slices were shown for the coronal orientation. For each of the 480 image pairs shown to participants, we recorded which method, if either, was preferred. For data analysis, the two methods were appointed one 'point' each for every time they had been preferred. For each neutral vote, both methods were appointed a half 'point'. This led to a dichotomous variable that allows for a direct comparison of the two methods' scores. A one-sided binomial test was conducted (H 0 : p our method ≤ 0.5, H 1 : p our method > 0.5). Table 2, shows the mean reconstruction rates for all parameter combinations.  The four-factorial ANOVA showed significant main effects for all four factors and one significant interaction effect for the reconstruction method and the threshold used (Table 3). This interaction effect describes that, while our method performs better than the baseline method at all threshold levels, it achieves greater improvements at higher thresholds (see also Fig 8 and Fig 9. Our method has a consistently higher reconstruction rate than the baseline (about twice as high) for all parameter sets. On the tested data, it was also more robust against the chosen similarity measure used for the template matching and also more robust against whether the reference sequence was acquired in the beginning or at the end of the session. Though, these interaction effects could not be shown to be significant in the ANOVA.   A correlation between acquisition order of the slice positions relative to the reference sequence and the ability of the methods to reconstruct these slice position can be seen in Fig 10. With increasing temporal distance between the acquisition of an interleaved sequence and the reference sequence, both methods find fewer similar slices for the corresponding slice position. Reference sequence one (red graphs) is acquired before the interleaved sequences. Here both methods find more slices for the earlier slice positions. Reference sequence two (blue graphs) is acquired after all interleaved sequences. Here both methods find more slices for the later slice positions. Blue graphs represent matches for the second reference sequence. Graphs with dots represent our method, graphs with crosses represent the baseline method.

Results
In the double-blind study, overall, participants selected our method in 156 trials, the baseline method in 111 trials, and had no preference in 213 trials (see Fig 11). Following our analysis method, this yielded 262.5 'points' for our method and 217.5 'points' for the baseline method (p=0.02).
The study shows that radiologists perceive the reconstruction quality of our method as significantly better than the baseline method, although the effect seems to be small.

Discussion and conclusion
The particular acquisition scheme shows difficulties with changes in breathing patterns that arise over a more extended period, like the typical flattening of the resting breath. Slice positions to the left are imaged only at the end of acquisition time, whereas slices to the right are only imaged at the beginning. As a consequence, if the reference sequence was captured in the beginning, it can show breathing states that do not occur later, when slice positions to the left are imaged. Deep breaths often can not be fully reconstructed since image data of the left slice positions was not acquired for deep breathing states. Generally speaking, the scheme has difficulties with breathing states that are less frequent. This problem can be solved in changing the acquisition scheme. Instead of first acquiring all slices in one position before moving on to the next slice position, it is beneficial to move the slice position after each acquisition while keeping the navigator position fixed. This rotating acquisition scheme could also be combined with intermediate reference sequences. This would directly counter the problem with flattening breath over time. Furthermore, with the new scheme, it is feasible to give a few commands, so the subject can take a few more deep breaths in the beginning before starting to relax more. The rotating acquisition scheme was used by Siebenthal et al. on a 1.5T Philips Intera whole-body MRI system [14]. However, Siemens MRI machines do not allow this kind of scheme. A solution to the problem that is independent of the scanner used, is to use external respiratory signals instead of navigator frames. Preiswerk et al. [20] correlated 1D ultrasound with 2D and multiplanar MRI. This allows for a continuous rotating acquisition of the data slices on any MRI machine. This also would reduce the acquisition time by half. Celicanin et al. [21] propose a MRI sequence that allows for the simultaneous acquisition of navigator and data frames increasing the temporal coherence of navigator and data frame. This would integrate well with the rotational acquisition scheme.
Regarding the reconstruction rate, because of the lack of ground truth, it is not possible to account for false negatives and false positives in the evaluation. Taking this fact into account, the reconstruction rate of both methods is likely to be higher than measured in this study. This is because in our test data, the number of reconstructable slice positions in each volume is lower than the number of slices in a volume, due to the aforementioned acquisition scheme.
An open issue arises when vessel cross-sections in the navigator frame are not continually visible. This frequently happens depending on blood flow. To solve this, one could detect outliers in the template matching step and omit those for the calculation of the summed displacement.
No body coil is used in our method to ensure a line of sight for external marker tracking. However, it would be ideal to use a flat and flexible interventional receiver coil to have better image contrast while still maintaining a line of sight for the marker tracking. This will become more feasible in the future when those coils become more widely used in interventional settings.
In summary, our results clearly show that template updates are an effective and efficient means to increase reconstruction rates and image quality of the reconstruction result for template based 4D MRI reconstruction methods. We reported that employing search regions greatly reduce reconstruction time. The results suggest that our method is preferable compared to the baseline. This is regardless of the applications favorable trade-off between precision and coverage because in all cases reconstruction rates are higher than the baseline.