Multiresolution radial MRI to reduce IDLE time in pre-beam imaging on an MR-Linac (MR-RIDDLE)

Online adaptive MR-guided radiation therapy improves treatment quality at the expense of considerable longer treatment time. The treatment lengthening partially originates from the preparatory (pre-beam) MR imaging required to encode all the information needed for contour propagation, contour adaptation and replanning. MRI requires several minutes of scan time before the encoded information is converted to usable images, which results in long idle times before the first clinical tasks are performed. In this study we propose a novel imaging sequence, called MR-RIDDLE, that reduces the idle time and therefore speeds-up the workflow in online MR-guided radiation therapy. MR-RIDDLE enables multiresolution image reconstruction to commence during data acquisition where low resolution images are available within one minute, after which the data collection continuous for subsequent high-resolution image updates. We demonstrate that the low resolution images can be used to accurately propagate contours from the pre-treatment scan. For abdominothoracic tumours MR-RIDDLE inherently captures a motion-blurred representation of the mid-position, which we were able to deblur using a combination of an internal motion surrogate and auto-adaptive soft-gating filters. Our results demonstrate that MR-RIDDLE provides a robust, flexible and time-efficient strategy for pre-beam imaging, even for cases with large respiratory movements or baseline shifts within the acquisition. We anticipate that this novel concept of parallelising the MR imaging and the clinical tasks has the potential to considerably speed-up and streamline the online MR-guided radiation therapy workflow.

before starting radiation therapy. Further, for abdominothoracic tumors the anatomy needs to be captured in a representative motion state (e.g. mid-position) (Wolthaus et al 2006(Wolthaus et al , 2008. Typical (diagnostic) sequences are unable to meet these requirements simultaneously. The ability to encode all the required information in a single MR scan that is also able to reconstruct the desired information at the correct moment in time would be invaluable for MR-guided radiation therapy.
In this work we describe the design and implementation of such a single MR scan that is optimised in conjunction with the clinical workflow to facilitate fast online treatment planning. We named the combination of the MR scan and the corresponding processing MR-RIDDLE: multiresolution radial MRI to reduce IDLE time in pre-beam imaging on an MR-Linac. The scan is based on a continuous 3D gradient echo golden angle stack-ofstars trajectory that enables reconstruction with flexible spatial and temporal resolution (Winkelmann et al 2007, Block et al 2014. The golden angle stack-of-stars trajectory enables image reconstruction to commence during data acquisition, where the reconstruction resolution is flexible and optimised with respect to the amount of data acquired at each instance in time. As a consequence, low resolution images are available after just a few seconds of scanning, which we hypothesize can be used to propagate the pre-treatment contours and therefore accelerate the clinical workflow (Glitzner et al 2015, Pandey et al 2017. The proposed method extends to anatomies affected by respiratory motion by combining a free-breathing acquisition with a motion-weighted image reconstruction. The stack-of-stars trajectory inherently portrays the time-averaged (blurred) position (mid-position) (Stemkens et al 2016), while simultaneously providing a selfnavigated motion surrogate (Spincemaille et al 2011). The combination of these two enables a motion-weighted image reconstruction that is able to reconstruct a deblurred mid-position volume with minimal processing time (Johnson et al 2012. The method will be especially attractive for patients that are unable to suspend respiration for a sufficient amount of time. We envision that our proposed method will provide a robust and relatively simple approach for online mid-position-based MR-guided radiation therapy treatment planning. Figure 2 provides an overview of the MR-RIDDLE pre-beam imaging method, which consist of radial k-space data acquisition, multiresolution image reconstruction and motion-weighted image reconstruction. K-space data are continuously acquired during free-breathing using a golden angle stack-of-stars trajectory. The image reconstruction commences at multiple time-points during the data acquisition and aliasing-free images are obtained by optimizing the reconstruction resolution with respect to the amount of acquired k-space data. Upon reaching the maximum spatial resolution, respiratory motion surrogates are estimated directly from the k-space data to perform motion-weighted image reconstructions that reduce motion-induced blurring. Note that the proposed method is a form of an expanding window image reconstruction and not a (conventional) sliding window reconstruction.

MRI data acquisition
MRI data were acquired using a golden angle stack-of-stars trajectory where all phase encodes along k z were scanned sequentially. The trajectory has four major advantages for the radiation therapy workflow. First, the radial trajectory repeatedly samples the k-space centre and therefore is inherently robust against motion- Figure 1. Overview of the online MR-guided workflow for a patient with a spine tumor. First the daily MR scan is required (≈ 4 min). The MR scan is registered to the pre-treatment CT or MR scan (≈1 min). The obtained deformation is used to propagate the contours to the daily MR. The contours are manually adapted to the current anatomy (≈10 min, dependent on anatomy). The updated contours are used to generate a new treatment plan (≈5 min). When the treatment plan is almost finished a second MR scan is acquired for position verification (≈2 min). If no large shifts occurred between the first and the second MR scan the treatment can commence.
induced artefacts (Glover andPauly 1992, Block et al 2014). Second, each readout contributes equivalently in k-space content to the complete image and therefore inherently portrays the time-averaged (blurred) position of the anatomy (Stemkens et al 2016). Third, the repeated sampling of the k-space centre enables extraction of a motion-state surrogate signal that allows a motion-weighted image reconstruction to effectively deblur the motion-blurred anatomy (Johnson et al 2012). Fourth, the golden angular increment ensures near uniform k-space filling at any instance in time, which enables image reconstruction from arbitrary imaging windows (Winkelmann et al 2007).

MRI data processing
Radially acquired MR data are prone to geometrical distortions and image artefacts due to gradient systems imperfections (Block and Uecker 2011, Campbell-Washburn et al 2014, Moussavi et al 2014. To mitigate the artefacts, the zeroth and first spatial order gradient impulse response functions of the gradient system were characterised and used for k-space trajectory correction and phase corrected image reconstructions (Vannesjo et al 2013, Bruijnen et al 2018. The centre k-space sample along the partition direction (kz) was used as a surrogate for the global motion (Spincemaille et al 2011). The surrogate was calculated by concatenating the projections for all the receivers coils into one 2D matrix and applying principal component analysis to detect the main variation in the frequency range (0.0-0.5 Hz) (Pang et al 2014). The motion surrogate was determined for each radial angle in the partition and thus the temporal resolution of the surrogate depends on the number of partitions (kz) times the repetition time.

MR-RIDDLE reconstruction for non-moving organs
The proposed method reconstructs aliasing-free images at multiple time-points during the acquisition by optimizing the reconstruction resolution with respect to the amount of acquired data. The reconstruction resolution is controlled by filtering the in-plane k-space with a block filter with isotropic width (W) that corresponds to reconstruction resolution (N). The relationship between the width of the block filter and the number of acquired projections (N proj ) for golden angle radial stack-of-stars is defined by equation (1), assuming no undersampling in the partition direction and isotropic in-plane image size. The final image resolution will equal the acquisition resolution (N acq ).
(1) Figure 2. Overview of the MR-RIDDLE protocol: K-space data are acquired during free breathing using a golden angle stackof-stars trajectory. Top row: The acquired k-space data is block-filtered (red) such that the residual k-space (blue) complies with the Nyquist criterion and aliasing-free (multiresolution) images can be reconstructed at any moment in time. Upon reaching the acquisition spatial resolution, the k-space data is motion-weighted using a soft-gating filter to reduce motion-induced blurring (middle row). The soft-gating filter is parametrised with the coefficients c 1 and c 2 . The light-blue in the soft-gating filters indicates the downweighting of the data, which is also shown in the middle row.
After in-plane resolution filtering, the k-space data of the individual coils were multiplied by a fixed Ram-Lak filter for density compensation followed by the non-uniform fast Fourier transform (NUFFT) (Fessler et al 2003). The coil images were combined using a Roemer coil combination (Roemer et al 1990) where the coil sensitivity maps were estimated from the low resolution images using ESPIRiT .

MR-RIDDLE reconstruction for moving organs
For abdominothoracic anatomies the multiresolution approach extends to the motion dimension in the form of a motion-weighted image reconstruction (soft-gating) (Johnson et al 2012). Imaging data that are acquired far away from the target motion position are downweighted in the reconstruction using soft-weights. These soft-weights (sw) are optimised with respect to the acquired k-space data (y) for a certain set of parameters (c = [c 1 , c 2 , c 3 ]) for the exponential weighting function shown in equation (2)  .
The soft-weights are derived from both the set of parameters c and the normalized distance d(y) between the current motion state and the target motion state. Here d(y) is defined as a normalized proxy (d(y) ∈ [0, 1]) for the physical distance to the target motion state. The threshold c 1 defines the region that is considered motionfree and parameter c 2 controls the downweighting in the regions where |d(y)| > c 1 . In addition, the penalty term c 3 k r (k 2 r = k 2 x + k 2 y ) differentiates between low and high frequencies within a readout to equalise their contribution to the total artefact power (Weiger et al 1997). The parameters c 1 and c 3 were fixed to c 1 = 0.5 and c 3 = 0.1 based on empirical observations, while c 2 was adapted dynamically over the acquisition. Soft-gating is in essence a sub-sampling operator and therefore c 2 can be optimised to fulfil the Nyquist theorem at any point in time, which for the stack-of-stars case can be defined as in equation (3). Here, N samp is defined as the number of samples acquired per k-space readout. The optimization effectively imposes the maximum amount of soft-gating without sub-Nyquist sampling. The search for the optimal c 2 can be approximately solved very fast using a bruteforce search in the parameter space ∈ [0 : 0.01 : 20]. The k-space data were then multiplied with the soft-weights followed by the image reconstruction as explained in the previous paragraph.

Digital phantom (XCAT) simulations on non-moving organs: proof of concept
To establish whether MR-RIDDLE reconstructs aliasing-free images from arbitrary imaging windows, we performed in-silico experiments using the 4D extended cardiac-torso (XCAT) digital phantom (Segars et al 2010). In short, the simulations were single receive channel acquisitions on discrete temporal volumetric XCAT data and did not include a magnetization model. The k-space data were obtained by applying the adjoint NUFFT and adding complex Gaussian noise ℵ(0, 0.5% |k(0, 0, 0)|). A detailed description of the simulation setup can be found in supplementary information I (stacks.iop.org/PMB/64/055011/mmedia).
Multiresolution images were reconstructed at 10 s intervals, using all the available data, and were qualitatively analysed for residual aliasing artefacts. The multiresolution images were effectively reconstructed with a different aperture and therefore show different point-spread functions.

Digital phantom (XCAT) simulations on moving anatomy: mid-position validation
To assess the capability of MR-RIDDLE to cope with respiratory motion and tumour drift in the thorax, we performed in silico experiments similar to the non-moving experiment. The XCAT phantom was parametrised such that the maximum respiratory amplitude of the diaphragm in the feet-head direction was at the higher end of physiological reported motion (3 cm) (Stam et al 2013, Heerkens et al 2017. In addition, we superimposed a gradual baseline shift of up to 2 cm. We inserted a spherical tumour (1 cm diameter) in the left lung that moved along with the prescribed motion. A detailed description of the simulation setup can be found in supplementary information I.
Images were reconstructed every 10 s using all the available imaging data, where the spatial resolution and extent of soft-gating increased in time. These reconstructions were used to assess the ability of MR-RIDDLE to effectively deblur the motion-induced image artefacts while portraying the tumour mid-position. The midposition was defined as the average position over the whole examined period. The mid-position of the tumour was estimated by registering the reconstructed volumes to the XCAT exhale volume. The registration was done using a masked, cross-correlation, based rigid registration method (Guizar-Sicairos et al 2008) and the tumour position in feet-head was extracted. The measured tumour position from the reconstructed (multiresolution) volumes was then compared against the analytical tumour positions provided by XCAT (supplementary information II). Note that the proposed method is a form of an expanding window image reconstruction and not a (conventional) sliding window reconstruction.

MRL imaging: contour propagation to low resolution images
To demonstrate one example of the potential workflow gain of MR-RIDDLE we implemented the sequence on an Elekta Unity 1.5 T MR-Linac (Elekta AB, Stockholm, Sweden). The key idea is that a low resolution image already suffices to accurately propagate the organ contours from the pre-treatment scan and therefore speed-up the clinical workflow. We acquired multiple MR-RIDDLE scans in six volunteers in the upper abdomen with different contrasts. The sequence parameters for the scans are specified in supplementary information I. MR-RIDDLE images were reconstructed at 10s intervals, using all the available data, without any image acceleration such as parallel imaging. The final reconstructed image was then selected as the reference and used to delineate the liver, the right kidney and the aorta. Note that the reference is from the same acquisition, i.e. not from a separate scan. All the reconstructed images were then deformable registered to the reference image using an in-house validated optical flow algorithm (α = 0.3) (Zachiu et al 2015a, 2015b). The resulting deformation vector fields were used to propagate the contours to the MR-RIDDLE images. The propagated contours were quantitatively compared on each time-point versus the reference using the Dice coefficient and the mean deformation vector (|DVF xyz |) (equation (4)). Note that in the case that no significant drift occurred within the 180 s scanning the Dice should be close to 1 and the |DVF xyz | close to 0 for all time-points.

Digital phantom (XCAT) simulations on non-moving anatomy: proof of concept
The relationship between the acquisition time and the in-plane (kx, ky) spatial resolution of the reconstructed images is shown in figure 3(A). Images with 5 mm resolution were reconstructed within 25 s of scanning time, while 180 s were required to achieve acquisition resolution (0.8 mm). The effect of the multiresolution filter (i.e. block filter width) is shown in figure 3(B), which shows reduced Gibbs ringing and a narrower point spread function with increasing scan time. The multiresolution images showed no residual aliasing artefacts at any reconstructed time-points (figure 3(C)). The signal-to-noise (SNR) increased with acquisition time linearly until reaching the acquisition resolution. The SNR increases linearly because both, the widening of the block filter and the prolonged scanning, include more samples for image reconstruction.

Digital phantom (XCAT) simulations on moving organs: mid-position validation
Multiresolution images were reconstructed and showed no aliasing artefacts (figure 4(A)). From left to right the spatial resolution increases and the motion-induced blurring reduces by gradually increasing the softgating. The reconstructions at time-points <10 s showed differences in tumour position of more than 5 mm compared to the true mid-position in the feet-head direction ( figure 4(B)). This difference in tumour position is presumably caused by the inability of the reconstruction and registration algorithms to function at very coarse image resolution (12 mm). Note that the proposed method was able to follow the systematic tumour drift inferred by the motion waveform shown in the top of figure 4(B). The tumour position was within 1 mm of the true mid-position after 20 s of acquisition time with mean difference of 0.6 ± 0.3 mm in the feet-head direction.

MRL imaging: contour propagation to low resolution images
Multiresolution images were reconstructed on a five second interval and typical image quality for spoiled gradient echo and balanced gradient echo are shown in figure 5. The final reconstruction was used to delineate the liver, the right kidney and the aorta. The reconstructions on different time-points show the same improvements in image quality as observed in the simulations. The image acquisition time is shown in the left bottom corner of the axial slices while the additional reconstructions times were approximately 10 s for the first volume and 20 s for 3 min of acquisition. The multiresolution images were deformable registered to the final reconstructed image using optical flow and the obtained DVFs were used to analyse the propagated contours for the aorta, kidney and the liver. The mean |DVF xyz | were around 0.5 mm after 30 s of scanning time and gradually approached zero over all six vol-unteers. Note that the mean |DVF xyz | should converge to zero, because the final image is the reference. The DVFs were used to warp to propagate the contours and the corresponding Dice coefficients were larger than 0.97 after 30 s of scanning and were approximately 1 after 150 s.

Discussion
Here, we are the first to design and to implement a method that encodes all the information required for treatment planning in a single pre-beam scan to reduce the idle time in the clinical workflow. In particular, we showed in digital simulations that MR-RIDDLE produces usable (low-resolution) images after 30 s of imaging, after which higher-resolution image-updates are available on demand. Then we showed with a 4D digital phantom that MR-RIDDLE inherently captures the mid-position of the anatomy, with a precision better than 1 mm, even for cases with large respiratory movements or baseline shifts. Finally, we demonstrated on a 1.5 T MR-Linac, in six volunteers and multiple contrasts, that MR-RIDDLE reliably produces high quality mid-position volumes for initialization of contouring in well under a minute without introducing additional registration uncertainties.
The major advantage of MR-RIDDLE is the ability to provide usable low-resolution images for contouring within one minute. These low-resolution images reduce the time required to start recontouring with 2:34 min for the protocol used in Raaymakers et al (2017), while simultaneously allowing even higher resolution acquisitions (supplementary information III). In addition, MR-RIDDLE enables images to be reconstructed from arbitrary temporal windows, which could be used for position verification scans that are commonly performed just prior to irradiation. These position verification scans would be instantly available on demand with our proposed method, which further reduces the idle time in the workflow (1:40 min for the example in figure 1). Moreover, MR-RIDDLE does not rely on prospective motion compensation techniques (e.g. gating). Instead, our free-breathing approach captures all the information and as a consequences makes the method applicable even to patients that are not able to suspend their breath. In addition, the soft-gated image reconstruction is able to reconstruct a deblurred mid-position volume with minimal processing time (<5 s). The proposed method considerably simplifies the currently used mid-position method, which requires the reconstruction of a 4D volume followed by multiple 3D image registrations (Wolthaus et al 2008). We envision that our proposed method will provide a robust and relative simple approach for online mid-position based MR-guided radiation therapy treatment planning.
Aside from pre-beam imaging, MR-RIDDLE could potentially be used as a real-time method to track the mid-position of slower moving anatomy, such as the prostate. In these type of applications we could increase the temporal resolution (decrease spatial resolution) together with sliding window reconstructions to provide robust and consistent large field-of-view 3D volumes for tracking. The sliding windows used for these reconstructions are completely flexible and can be tailored in terms of temporal and spatial resolution for targeted applications. For these reconstructions we could use a 95 DICE score cut-off (figure 6) and get accurate deformation vector fields with a temporal window of only 22 s without acceleration methods.
Besides online MR-guided radiation therapy, MR-RIDDLE could be used in diagnostic scans for patients that are unable to remain sufficiently still. For these patients the longest scanning period without bulk motion could be detected in retrospect and the motion-free images can be recovered using a lower-resolution reconstruction, which would reduce the number of failed scans. The current implementation of MR-RIDDLE supports spoiled gradient echo or balanced gradient echo contrast, but not pure T2 contrast. However, emerging approaches to generate T2 contrast with stack-of-stars sequences are promising (Benkert et al 2017).
MR-RIDDLE has several limitations that require discussion. First, the method enables high resolution prebeam acquisitions at a cost of a minor increase in repetition time compared to conventional low-resolution acquisitions. The increase in repetition time is a function of machine hardware parameters and is discussed in supplementary information IV. Second, the quality of the soft-gated image reconstruction to deblur the anatomy in the mid-position depends largely on the motion surrogate. Systematic shifts of this motion surrogate could introduce a position bias, however the inclusion of multiple surrogates could tackle this issue. Note that the MR- RIDDLE reconstructions will always have a small delay in detecting these systematic shifts due to the intrinsic blurring caused by the temporal window. Third, MR-RIDDLE is not easily translatable to spin-echo type of sequences that are commonly used in radiation therapy.
So far, MR-RIDDLE was only evaluated in volunteers in an offline setting. In future work we aim to embed the image reconstruction in the online workflow to perform treatment simulations. Using the treatment simulations we can begin to determine optimal reconstruction timings, investigate the convergence to sharp reconstructions in abdominal/thoracic regions and to examine the overall practicality of the method in the clinical workflow. Moreover, MR-RIDDLE could be extended to reconstruct respiratory resolved 4D MRIs from the exactly the same dataset (Stemkens et al 2015. These respiratory-resolved 4D MRIs do not capture systematic baseline shifts and therefore could be supplemented with a low spatial/high temporal resolution sliding window (as discussed in the real-time method) for complete motion estimation. Both these extensions are excellent examples of the potential of the flexibility of the golden angle sampling scheme that allows the user to encode information of multiple temporal and spatial resolutions in a single scan.
The current implementation of MR-RIDDLE did not use acceleration techniques such parallel imaging due to the long processing times required for 3D non-Cartesian sensitivity encoding. The longer processing time would counteract the saved acquisition time and therefore not speed-up the workflow. An alternative method to accelerate the imaging would be to use partial Fourier, Cartesian sensitivity encoding in the partition direction or use variable density or asymmetric field-of-view stack-of-stars sampling schemes. Using these methods we believe to gain a factor of 2-3 acceleration on top of the imaging times reported in this work.

Conclusion
MR-RIDDLE is capable of producing a mid-position volume for initialization of contouring or contour propagation in well under a minute, after which the image acquisition continues to collect data for highresolution image updates. Moreover, continuing the data acquisition during the recontouring enables the final image reconstruction to be very high resolution, which is infeasible with conventional sequences. MR-RIDDLE creates a more streamlined and time-efficient workflow where continuous data acquisition and human interaction are performed in parallel. Figure 6. Analysis of the propagated contours for the multiresolution MR-RIDDLE reconstructions. (A) The methodology used to assess the contours. All the reconstructed images were registered to the final reconstruction using optical flow (α = 0.3). The obtained DVFs were analysed using the total displacement vector in the aorta, kidney and liver. Additionally, the DVFs were used to deform the organ segmentations to calculate the corresponding Dice coefficients. (B) The mean displacement vector and mean Dice coefficient across six volunteers are shown as a function of the acquisition time.