3D model-based super-resolution motion-corrected cardiac T1 mapping

Objective. To provide 3D high-resolution cardiac T1 maps using model-based super-resolution reconstruction (SRR). Approach. Due to signal-to-noise ratio limitations and the motion of the heart during imaging, often 2D T1 maps with only low through-plane resolution (i.e. slice thickness of 6–8 mm) can be obtained. Here, a model-based SRR approach is presented, which combines multiple stacks of 2D acquisitions with 6–8 mm slice thickness and generates 3D high-resolution T1 maps with a slice thickness of 1.5–2 mm. Every stack was acquired in a different breath hold (BH) and any misalignment between BH was corrected retrospectively. The novelty of the proposed approach is the BH correction and the application of model-based SRR on cardiac T1 Mapping. The proposed approach was evaluated in numerical simulations and phantom experiments and demonstrated in four healthy subjects. Main results. Alignment of BH states was essential for SRR even in healthy volunteers. In simulations, respiratory motion could be estimated with an RMS error of 0.18 ± 0.28 mm. SRR improved the visualization of small structures. High accuracy and precision (average standard deviation of 69.62 ms) of the T1 values was ensured by SRR while the detectability of small structures increased by 40%. Significance. The proposed SRR approach provided T1 maps with high in-plane and high through-plane resolution (1.3 × 1.3 × 1.5–2 mm3). The approach led to improvements in the visualization of small structures and precise T1 values.


Introduction
Cardiovascular MR is a well-established technique for the diagnosis of cardiac diseases. Over the last years, T1 mapping has been translated into clinical application as an important quantitative approach for cardiac tissue differentiation (Guo et al 2022). It has been demonstrated that T1 mapping can be used to diagnose a wide range of different cardiac pathologies, including entities with preserved ejection fraction (Haaf et al 2016, Schelbert and Messroghli 2016, Al-Wakeel-Marquard et al 2021. However, one of the major challenges in cardiac T1 mapping is that the achievable image resolution is often restricted due to respiratory and cardiac motion, signal-to-noise ratio (SNR) and limited acquisition time. In clinical practice, T1 maps are acquired with a 2D acquisition scheme resulting in one slice per breath hold (BH) that has a high in-plane but a poor through-plane resolution of 6-8 mm (Becker et al 2019). 3D T1 Mapping has been proposed but with a long acquisition time (Qi et al 2019(Qi et al , 2020. With the work proposed in Becker et al (2020), 80% of the cardiac cycle can be used for T1 mapping, allowing the acquisition of six slices with a slice thickness of 6-8 mm per BH. Nonetheless, image resolution is compromised by partial volume effects. This can impair the accurate detection of subtle fibrosis in the myocardium in different entities and limit the capability to differentiate myocardial injury within the thin myocardial wall of young patients.
Super-resolution reconstruction (SRR) has been proposed to improve the tradeoff between spatial resolution, acquisition time and SNR (Greenspan et al 2002, Shuzhou et al 2006 de Senneville et al 2020, Ebner et al 2020, Sui et al 2019, 2021. The resolution is thereby increased by acquiring various lowresolution (LR) images with complementary information about the object. This is ensured by shifting the image positions of the LR stacks along the slice direction or by changing their slice orientations. Subsequently, the LR stacks are combined into a high-resolution (HR) dataset by solving an inverse problem. The reconstructed HR image thus benefits from the high SNR of the LR images while providing HR diagnostic information. For quantitative MRI, the parametric model of the mapping can be combined with the SRR model which has been demonstrated on the brain (Van Steenkiste et al 2017, Bano et al 2020). Such a model-based SRR enables the direct estimation of HR T1 maps from LR T1-weighted images (dynamics).
The principle of SRR is based on knowledge about the geometric relationship between different LR datasets. Motion leads to misalignment and strongly impairs the achievable image quality of SRR ( The application of SRR to quantitative cardiac MRI data is especially challenging due to cardiac and respiratory motion. Next to that, the acquisition during BH imposes severe limitations on the acquisition time and thus limits the number of slices per LR stack. So far, no model-based SRR T1 mapping has been applied on cardiac data, which requires advanced acquisition and motion correction (moco) schemes.
In this study, we present a model-based SRR for cardiac T1 mapping, providing precise HR T1 maps with improved visualization of small structures compared to the direct LR acquisitions. It combined multiple stacks of 2D acquisitions with 6 to 8 mm slice thickness and generated 3D HR T1 maps with a target slice thickness of 1.5-2 mm in six to ten BH. Cardiac and residual respiratory motion was corrected. The approach was evaluated in native T1 mapping in numerical simulations and phantom experiments and feasibility was demonstrated in four healthy volunteers.

Methods
The proposed workflow to achieve motion-corrected model-based SRR T1 maps is depicted in figure 1: multiple stacks of 2D slices were acquired continuously with one stack per BH. In a first step, non-rigid cardiac motion was estimated and used in a model-based T1 reconstruction (Becker et al 2019) resulting in the dynamics w g and parameter maps m g (6 slices à 6-8 mm per stack) which were all in the same cardiac motion state. In a second step, the stacks were registered to each other to estimate and compensate for different BH positions. After the motion alignment, the maps were then used to calculate the first estimate of the HR map m 0 G as initialization of the SRR. Finally, a HR T1 map m final G was generated by SRR.

Data acquisition
Data was acquired using a Golden-angle radial sampling scheme on 3 Tesla MR scanner (Verio, Siemens Healthineers, Erlangen, Germany) with a commercial 32-channel cardiac coil. After a slice-selective radiofrequency inversion pulse, data was continuously acquired in multiple stacks with six slices each resulting in an acquisition time of 16.8 s for a single stack (2.8 s for each slice) with the following parameters: flip angle α = 5°, resolution 1.3 × 1.3 × 6.0-8.0 mm 3 , field of view (FOV) 320 × 320 × 84-105 mm 3 , TE/TR: 2.19/4.9 ms, orientation short-axis-view, subject specific slice gap of 4-9 mm to cover the desired FOV while avoiding slice interference from the radio-frequency inversion and excitation radio-frequency pulses. Six to ten stacks (one stack per BH) were acquired in total with an offset of 1.5-2 mm between stacks along the slice direction. Due to the short acquisition time, a slice-selective inversion pulse was used in combination with an interleaved multislice ordering. The ECG was recorded for retrospective cardiac moco.

Model-based T1 reconstruction
Dynamic cardiac motion-resolved images were reconstructed with a temporal resolution of 44.1 ms. Spatial and temporal total variation regularization (regularization parameters l along time and space were 0.5) was applied to suppress undersampling artefacts (Block et al 2007). To accelerate the motion estimation, a subject specific rectangular region of interest covering both ventricles was selected. In an iterative fashion, the non-rigid cardiac motion was estimated using the MIRTK Toolkit (Rueckert et al 1999).
The estimated cardiac motion information was used in an iterative model-based T1 reconstruction (Becker et al 2019(Becker et al , 2020. Data from the entire cardiac cycle was used, except for the 30% of the systole with the greatest through-plane motion. A Look-Locker model q was used in an iterative reconstruction scheme to estimate m g with the quantitative parameter m p, , T1 [ ] a = and w g with a temporal resolution of 83.3 ms, p denoting the equilibrium magnetization and a the flip angle. In the following, only the T1 parameter is mainly considered because it is clinically the most relevant.

BH registration
Each stack was acquired in a different BH. To correct for potential misalignments of BH positions, the stacks were registered to each other using a cross-correlation approach (Padfield 2012). A two-stage process was developed for this purpose (figure 2). In the first step of the motion estimation, the rigid motion in the in-plane direction of the slices s m g was determined. For that, the T1 maps of the LR slices of stacks s were registered to each other: each slice of each stack was registered to the slice which was closest (i.e. smallest distance along the slice direction) to it. The stacks were acquired in an overlapping fashion, therefore the closest slice was part of another stack and hence, s T1 g was registered to s 1 T1 ( ) gusing a phase-cross-correlation registration. That yielded information about the in-plane motion of every slice of every stack. The median of the motion detected in its six slices was finally assigned to the entire stack of LR slices . In the second step, the LR stacks were registered with respect to shifts along the slice encoding direction. For that, s T1 g was interpolated along the slice encoding direction using bicubic spline interpolation, which also filled the gaps between the LR slices. The interpolated T1 maps of the LR stacks were then combined and an average stack avg T1 g was calculated. In an iterative process, each stack was then registered to . avg T1 g In the next iteration a new avg T1 g was calculated taking the estimated motion into account. Only translational shifts were considered.
Two iterations were used in total.

Model-based SRR
For SRR, several LR stacks acquired with an offset to each other were combined to a HR volume. The model q calculates dynamics from given parameter maps. The SRR used here is model-based and thus q was incorporated into SRR to at the end obtain a HR T1 map m G from . w g As initialization m 0 G of the SRR, m g were calculated from w g using a voxel-wise three-parameter T1 fit and combined: Figure 1. Comparison of the proposed motion-corrected model-based SRR workflow and the common approach. In this schematic comparison, data was acquired over eight breath holds (BH). In the common approach, one slice could be reconstructed per BH. In the proposed approach, one stack per BH with six 2D slices each was acquired. Cardiac motion was estimated and included in a model-based T1 reconstruction of the k-space data k yielding the dynamics w g and the parameter maps m g of the LR stacks. Then, the different stacks were registered to each other. The motion corrected m g were used to calculate the first estimate of the HR map m 0 G and SRR yielded the final 3D HR parameter map .
m final G In this example, the proposed approach led to six times more slices with the same number of BH, with a slice thickness reduced by a factor of four compared to the common approach. g ) and a total variation based regularization term was minimized, which could be described by the following minimization problem: where k describes the regularization parameter and G corresponds to the forward finite differences operator. As a s h l , describes the relationship between a HR and LR slice, by solving problem (3) an estimate of the HR slices could be recovered. Since solving problem (3) directly is challenging due to the non-smoothness of the L1-norm as well as the non-linear function q, a variable splitting (Wang et al 2008, Bano et al 2020 approach was used. This allowed solving the resulting sub-problems with suitable algorithms. By introducing auxiliary variables x q : t t m ( ) = G for all t and u: m =G the problem was reformulated as a joint minimization problem. These equalities were relaxed by including two quadratic penalty terms, weighted by l and , m yielding:  The solution of problem (4) was approached by alternating the minimization of (4) with respect to one of the variables and keeping the other two fixed. For fixed u , , Subproblem (4a) was minimized with respect to x, assuming m G and u are fixed. Solving (4a) involved solving a linear system for which a conjugate gradient approach was used. For was solved using the iterative algorithm proposed in Chambolle (2004).
For fixed u and x, updating m G in problem (4) corresponded to solving Due to the non-linearity of function q, the Limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm (Liu and Nocedal 1989) was used for solving problem (4c). To solve problem (4), the subproblems were alternated eight times and the solution of problem (3) was referred to as .

Simulation experiments
Simulated data was generated using the XCAT phantom (Segars et al 2010). A dataset X orig with the voxel size 1.3 × 1.3 × 0.5 mm was generated. From this, eight stacks of LR dynamics were simulated with the same parameters used for the phantom and the in vivo experiments (slice thickness of 6 mm, a gap between the LR slices of 6 mm and an offset between the stacks of 1.5 mm). As reference X ref , a dataset with a slice thickness of 1.5 mm was generated from X orig . In X orig , two cubical fibrotic structures were simulated in the septum with a width of 6 mm each along the slice encoding direction. They were separated by a 6 mm gap of healthy myocardium. Cardiac motion was simulated using the default settings of the XCAT phantom. Data acquisition was simulated with multiple receiver coils using the same acquisition parameters as for the phantom and the in vivo experiments. Zero-mean noise was added. This allowed the application of the entire pipeline including model-based T1 reconstruction and cardiac moco on the simulated data.
In the simulations, a T1 time of 1300 ms was assigned to the myocardium, 400 ms to fat, 800 ms to the liver, 900 ms to muscle and 1800 ms to the simulated fibrosis. Blood was simulated with an apparent T1 time of 350 ms, as it could not be estimated due to the in-flow effect caused by the slice-selective inversion pulse (Keith et al 2017).
Two different types of simulations were carried out. First, a simulation assuming perfect BH (i.e. no misalignment between different BH) was performed to evaluate the possible improvement that could be achieved with SRR. In a second step, misalignment between the BH was included into the simulation. Different BH positions were simulated by applying translation shifts. 20 configurations with different breath-hold positions of the stacks were simulated. The simulated motion was in the range of (3.5, 1.9, 8.2) mm in the (anterior-posterior, right-left, food-head) direction, based on half of the motion range between end expiration and end inspiration measured in (Scott et al 2009). For reasons of computational time, no heart motion was included in this simulation.
To assess the outcome of the SRR applied on simulated data, the detectability d between the simulated fibrosis and surrounding healthy myocardium was measured, using the following formula: where the mean T1 value structure m was measured in a region-of-interest (ROI) within the fibrosis, the mean T1 value nextStructure m was calculated in a ROI next to the fibrosis and the standard deviation (SD) background s was calculated from a ROI in the surrounding healthy myocardium. The ROI was calculated from the position of the simulated fibrotic structures in X . orig To evaluate the breathing moco, the RMS error Î  between the originally simulated motion and the estimated motion was calculated in mm. The LR data in the proposed approach was always acquired in the same orientation according to the shortaxis-view of the heart. To assess through-plane resolution, images in this publication are often presented orthogonally from the side, resulting in a four-chamber view.

Phantom experiments
To evaluate the proposed approach in phantom measurements, imaging was performed with the abovedescribed scan parameter in a 'T1MES-phantom' with nine tubes with different T1 times developed for cardiac imaging (Captur et al 2016). In order to cover the whole phantom with every stack, 12 slices per stack were acquired. Furthermore, a scan m orth g orthogonal to the LR data was acquired, where the slice encoding direction of m g became an in-plane direction of . m orth g As a reference, an inversion recovery spin-echo ref SE was acquired, also in orthogonal direction to the LR data with 7 TIs between 25 and 4800 ms (TE/TR: 12/8000 ms, FOV: 143 × 160 mm 2 , spatial resolution: 0.8 × 0.8 × 5 mm 3 ).
To assess the outcome of the SRR applied on the T1MES phantom, a ROI was drawn in every tube in ,  G is compared with an orthogonal acquisition . orth T1 g A line plot through three tubes (brown line) along the slice encoding direction (SE) shows an improved differentiation (pink arrows) between tubes and background after SRR as shown by the reference in green.

In vivo experiments
To evaluate the proposed approach in in vivo measurements, data was obtained from four healthy subjects (4 males, aged 34.0 ± 11.7 y). All subjects gave written informed consent before participation, in accordance with the institution's ethical committee. For the in vivo data, an orthogonal scan m orth g was acquired. For reference, a 3(3)3(3)5 modified Look-Locker inversion recovery (MOLLI) scan MOLLI G was acquired with the following scan parameter: FOV: 360 × 306 mm 2 , TE/TR: 1.12/2.7 ms, flip angle: 35°, and spatial resolution: 2.1 × 1.4 × 6 mm 3 once in four chamber (4CH) and once in two chamber (2CH) orientation. The T1 values of the SRR result and MOLLI reference were compared using a ROI placed in the septum.
To assess the outcome of the SRR applied on in vivo data, m orth g was qualitatively compared to . m final G The precision of the T1 values was evaluated quantitatively by comparing the bull's eye plots (Weissman et al 2002) before and after the SRR, using four selected slices (apex, apical, mid-cavity and basal) and calculating the SD over four healthy volunteers. No fibrotic tissue was present in the healthy volunteers and therefore the detectability of the right ventricle was calculated to assess the effect of SRR on small structures. Edge sharpness was calculated by manually drawing a line along the edge of interest, gathering the intensities perpendicular to the line, and generating an average edge profile. The first-order derivative of the edge profile was calculated and edge sharpness of 100% referred to the case when the maximum derivative of the average edge profile was equal to the maximum intensity difference in the average edge profile, similar to Etienne et al (2002). Figure 3 shows the results of the numerical simulations assuming perfect BH positions. In the LR stacks the two different fibrotic structures could not be distinguished along the slice encoding direction. The apex was inaccurately depicted in .   Figure 5 shows the in-plane view and the orthogonal reformation of , G are shown without moco, when the estimated motion was applied and when the reference motion was applied during moco. SRR without moco shows motion artefacts, which could be removed after applying the calculated moco. The visual result after applying the estimated motion shifts is similar to applying the reference motion shifts during moco. The RMSE between estimated and reference motion was (0.03, 0.04, 0.61) mm. Without BH alignment, motion artefacts could be seen in the form of a discontinuous septum in the orthogonal view and an ambiguous delineation of the myocardium in the in-plane view, which is highlighted by the pink arrows in the figure. The motion artefacts were less visible in the initialization of the SRR compared to its output. Using the proposed BH registration and subsequent correction, the motion artefacts after SRR could be reduced. Figure 9 shows the in-plane view and the orthogonal reformation of ,

In vivo experiments
G and compares it to an orthogonal acquisition orth T1 g and .

MOLLI
G Due to the slice-selective inversion pulse, blood appeared with a low T1 value. The visualization of the apex as well as the differentiation between the right ventricle and blood improved after SRR. Due to scan time limitations, orth T1 g could not be acquired for one volunteer. The mean T1 value across all volunteers in a ROI in the septum in final T1 G was 1211.49 ± 75.17 ms and in MOLLI G 1276.11 ± 38.77 ms. One volunteer had to be excluded from the calculation because no MOLLI scan was available.
In figure 10, four selected slices (apex, apical, mid-cavity and basal) of , G are compared in-plane.
The apex was more clearly visible in the apex slice after SRR. The visualization of the right ventricle also improved after SRR. In addition, the combination of multiple LR slices in the SRR also reduced artefacts (black arrow in figure 10) and improved e.g. the quantification of the inferolateral segment of the basal slice. d in the right ventricle increased from 2.4 ± 1.35 in , initialization than in the LR images, which could be attributed to the mixing of the partial volume effects in the individual LR images when combining them for the initialization. SRR was able to restore the original edge sharpness of the LR images. T1 maps of the three other volunteers can be seen in Supporting Information figure S1. Figure 11 shows the bulls-eye plots of , T1 g 0 T1 G and , final T1 G averaged over four healthy volunteers. The SD before and after the SRR remained comparable, indicating that SRR did not affect the precision of the T1 values. The T1 values in the segments varied in final T1 G by an average of 63.72 ms across the four healthy volunteers. The T1 intensities of the apex segment were underestimated before the SRR and showed a high SD. This was compensated by SRR.

Discussion
In this study, a novel motion-corrected model-based SRR approach was presented, providing 3D HR T1 maps in six to ten 17 s BH. The proposed model-based SRR scheme improved the visibility of small structures while the accuracy and precision of the T1 values after SRR remained comparably high. An alignment of different BH states showed great improvement of the SRR result. G are reformatted orthogonally and compared to a direct orthogonal acquisition orth T1 g and to a MOLLI reference scan . MOLLI G All the results shown were obtained with the proposed moco approach. The visualization of the apex and the right ventricle improved after SRR (pink arrows). Due to the slice-selective inversion pulse, blood appeared with a low T1 value. SE indicates slice encoding direction.
Small structures, as e.g. the simulated fibrosis, the differentiation between phantom tubes and background or the right ventricle could be better visualized using SRR. Furthermore, anatomic information which was impaired in some LR stacks due to partial volume effects, as e.g. the apex, was successfully recovered by the proposed SRR approach.
An improvement in the imaging of small features by SRR can be concluded from the improved visualization of small structures in all volunteers of the in vivo experiments, such as in the right ventricle.
The accurate mapping of the right ventricular myocardium poses a great challenge due to its small thickness but would help to improve the diagnosis of e.g. right ventricular myocarditis or arrhythmogenic right ventricular cardiomyopathy. Its assessment could be improved by SRR, moving towards whole heart T1 mapping in the future.
The results were compared to a clinical reference scan and the T1 values after SRR were in good agreement with both the reference values resulting from the modified Look-Locker inversion recovery reference scan and those presented in literature (von Knobelsdorff-Brenkenhoff et al 2013). The small underestimation of the myocardial T1 values after the SRR compared to reference values was probably due to the use of a slice-selective inversion pulse. A similar underestimation of the T1 values was reported in Huang et al (2020), which was attributed to magnetization transfer effects. A direct comparison of the SRR result to an in vivo reference scan was however difficult since this was acquired in another BH, hence showing a different motion state. Thus, the accuracy of the T1 values could only be determined in phantom measurement but not in the volunteer scans.
The proposed approach was not compared to previously published 3D T1 mapping frameworks as for example MR multitasking (Christodoulou et al 2018)  The precision of the T1 values was not calculated with a retest but as the SD over several healthy volunteers. It was assumed, that the T1 values of the myocardium were similar in all healthy volunteers.
One limitation of this approach is that T1 values of voxels representing blood could not be estimated and appeared shortened, due to the in-flow effect caused by the slice-selective inversion pulse. The slice-selective inversion pulse only inverted the blood spins that were in the corresponding slice at the time of the inversion. In the course of the cardiac cycle, these were replaced by inflowing, non-inverted blood spins, which then lead to an Figure 10. Four selected slices (apex, apical, mid-cavity and basal) before and after SRR. The visualization of the apex and the right ventricle improved in the SRR result final T1 G compared to its initialization 0 T1 G and a single LR stack T1 g (pink arrow in apex and midcavity slice). SRR reduced artefacts (black arrow in basal slice). All the results shown were obtained with the proposed moco approach. apparent shortening of the T1 value of the blood (Keith et al 2017). To still be able to calculate the extracellular volume, a further fast acquisition of a single LR slice with a global inversion pulse would provide the necessary information about the blood pool T1 values.
Every stack was acquired in a separate BH. Due to variations in BH, an alignment of different BH states was necessary before SRR. In agreement with Van Reeth et al (2012), motion estimation was a key step in the SRR process and significantly affected the quality of the SRR result. Imperfectness in moco could lead to artefacts after SRR. G is compared to a LR slice T1 g and the SRR initialization . 0 T1 G In clinical practice, 17 s BH are sometimes not feasible. To adapt the BH duration, the acquisition time per stack would need to be reduced and compensated for by acquiring more stacks in total. Due to the higher number of stacks, the proposed moco approach would have an even greater influence on the SRR result.
Compared to brain T1 mapping, cardiac imaging is restricted with respect to the number of LR slices per stack, due to limited BH time. To still cover a specific field-of-view in the slice encoding direction, gaps needed to be introduced between the LR slices. To compensate for these gaps, more stacks of LR slices needed to be acquired. According to Rahman and Wesarg (2010a), the more stacks used, the greater the degrading influence of inaccuracies in the motion registration on the SRR.
The stacks were planned such that they overlapped with each other by 1.5-2 mm. As each stack was obtained in a different BH position, the original distribution of stacks was impaired, even if the respiratory motion was correctly detected and estimated. Overall, the detectability of the simulated fibrosis was high in most of the motion corrected SRR maps. Nevertheless, depending on the distribution of stacks relative to each other after moco, the depiction of the fibrosis could still be impaired. This effect could in future be reduced by orientating the LR stacks differently (e.g. rotated to each other) along the slice encoding direction, which has also previously shown to improve SRR reconstructions compared to shifting the LR stacks (Shilling et al 2009).
The results of the moco could be improved by its integration into the optimization scheme of the SRR, as described in Dzyubachyk et al (2015). Furthermore, integrating rotation and deformation into the moco would probably further improve the alignment of BH states. Next to that, registering the slices within one stack separately to the HR volume would also account for inter stack motion due to poor breath holding. In addition, instead of retrospective BH correction, the position of the slices could be tracked prospectively and the acquisition adjusted accordingly, for example using the Pilot tone (Ludwig et al 2021). Next to that, SRR is not limited to T1 mapping, but could be extended to other quantitative parameters such as T2 (Giri et al 2009), for example using MR multitasking (Christodoulou et al 2018) or MR fingerprinting (Ma et al 2013). To improve the overall result of the SRR in future approaches, the SRR optimization scheme could be integrated in a modelbased reconstruction framework as performed in Bano et al (2020). By that, the SRR would incorporate the acquired raw data in the entire reconstruction optimization scheme instead of using it only in the model-based T1 reconstruction as in the presented approach.
This work was only evaluated in healthy volunteers, nevertheless, from the improved visualization of pathologies in the simulated data, it can be concluded that SRR might lead to an improved image quality in patients as well.

Conclusion
In this study, a novel motion-corrected cardiac model-based SRR approach was presented, providing 3D HR T1 maps in six to ten 17 s BH. The proposed approach was successfully applied in four healthy volunteers leading to improved visualization of small structures and precise T1 values. In future studies, an integration of the BH alignment and the T1 reconstruction into the optimization scheme could further improve the results.