Fetal whole-heart 4D imaging using motion-corrected multi-planar real-time MRI

Purpose: To develop a MRI acquisition and reconstruction framework for volumetric cine visualisation of the fetal heart and great vessels in the presence of maternal and fetal motion. Methods: Four-dimensional depiction was achieved using a highly-accelerated multi-planar real-time balanced steady state free precession acquisition combined with retrospective image-domain techniques for motion correction, cardiac synchronisation and outlier rejection. The framework was evaluated and optimised using a numerical phantom, and evaluated in a study of 20 mid- to late-gestational age human fetal subjects. Reconstructed cine volumes were evaluated by experienced cardiologists and compared with matched ultrasound. A preliminary assessment of flow-sensitive reconstruction using the velocity information encoded in the phase of dynamic images is included. Results: Reconstructed cine volumes could be visualised in any 2D plane without the need for highly-specific scan plane prescription prior to acquisition or for maternal breath hold to minimise motion. Reconstruction was fully automated aside from user-specified masks of the fetal heart and chest. The framework proved robust when applied to fetal data and simulations confirmed that spatial and temporal features could be reliably recovered. Expert evaluation suggested the reconstructed volumes can be used for comprehensive assessment of the fetal heart, either as an adjunct to ultrasound or in combination with other MRI techniques. Conclusion: The proposed methods show promise as a framework for motion-compensated 4D assessment of the fetal heart and great vessels.


| INTRODUCTION
Fetal cardiac MRI has been limited by the challenges associated with imaging a small, rapidly beating heart that is subject to various regular and spontaneous movements within the maternal torso. 1,2 However, as motion-tolerant fetal cardiac MR methods are introduced 3,4 there is an increasing capacity for MRI to visualize the fetal heart and great vessels. The success of retrospective methods for high-resolution 3D depiction of the fetal brain 5,6 suggests the potential of this approach for 4D reconstruction of the fetal heart in the presence of motion, allowing for robust cardiac evaluation.
Recent published studies have taken advantage of the flexibility of continuous golden angle radial sampling and the resolution that can be achieved with segmented cine imaging using compressed sensing for 2D imaging of the fetal heart, [7][8][9] with impressive results. In-plane motion-correction in k-space is possible with this approach, 3 however through-plane fetal and maternal movements cannot currently be corrected and in-plane motion across the entire field of view must be dealt with prior to or during cine MR image reconstruction.
This work builds on the motion-tolerant image-domain cine reconstruction framework previously described for 2D cardiac fetal imaging 4 and methods for 3D image volume reconstruction, 5 with the aim of developing an acquisition and reconstruction approach to generate a 4D cine representation of the fetal cardiovascular anatomy in utero from 2D multiplanar real-time MR images, without the need for maternal breath-hold or significant manual processing. Whole heart coverage presents challenges beyond 2D imaging, as motion must be corrected or compensated and the cardiac cycle synchronized for all acquired data is needed to allow for volumetric reconstruction.
In this work, 4D capability was achieved using a multiplanar real-time acquisition combined with retrospective image-domain techniques for motion correction, cardiac synchronization and outlier rejection. Preliminary results showed the proposed framework was capable of reconstructing 4D data from real-time MR images acquired without maternal breath-hold and without requiring precise scan plane prescription during acquisition. 10 In this work, the framework was validated using numerically simulated cardiac MRI data and evaluated in a study of 20 human fetal subjects, including comparison with matched ultrasound. A preliminary assessment of velocity-sensitive reconstruction was also performed.

| METHODS
The proposed framework for 4D whole-heart reconstruction is shown in Figure 1 and consists of motion correction using 'static' temporal mean images to achieve rough spatial alignment, followed by cardiac synchronization and further motion correction using 'dynamic' real-time images, and concluding with a final 4D reconstruction including outlier rejection. The symbols used to denote the acquired real-time F I G U R E 1 Framework for 4D cine volumetric reconstruction, consisting of (A) acquisition and reconstruction of multi-planar realtime 2D MRI; (B) an initial motion correction stage to achieve rough spatial alignment of the fetal heart using temporal mean (i.e., static) images for stack-stack registration followed by slice-volume registration interleaved with static volume (3D) reconstruction; (C) cardiac synchronization, including heart rate estimation and slice-slice cardiac cycle alignment; and (D) further motion-correction using realtime (i.e., dynamic) image frames interleaved with 4D reconstruction; and (E) a final volumetric reconstruction, including outlier rejection. User-specified ROIs, and identification of a target stack for stackstack registration are the only manual preparations required for reconstruction MR images and reconstructed 4D cine MRI throughout this manuscript follow the convention commonly used for sliceto-volume reconstruction (SVR) methods, 5,6 which differs from those used previously to describe the purely 2D cine framework. 4

| Multi-planar real-time MRI
Multi-planar data were acquired in stacks of parallel slices, as in previous SVR reconstruction methods. 5 Stacks were acquired in multiple orientations, ideally with the initial three stacks planned mutually orthogonal to each other to ensure full coverage of cardiac anatomy. For each slice a series of real-time images (or frames) was obtained, thus ensuring dense sampling in both space and time. Fetal and maternal movement were expected to cause shifting of the fetal heart relative to the imaging field of view. These shifts can be large on the pixel scale and in any direction.
The acquired data form a set of N k real-time MR image frames, Y = Y k k = 1, …, N k , where each frame has an associated acquisition time t k and consists of elements y jk at 2D spatial coordinates indexed by j. Subsets of Y for a single slice are defined as In the previous 2D+time framework 4 real-time MR images were reconstructed using k-t SENSE 11 with spatiallyadaptive regularization to preserve the full temporal resolution of the accelerated acquisition. These complex-valued real-time images were then combined as cine image series using robust statistics based on complex-valued errors to suppress voxels corrupted by artefact. This was possible because the phase of the complex-valued real-time images was fairly consistent in a single slice. However, in this 4D framework, the phase of the real-time images is inconsistent due to changes in slice orientation and position. As a result, these complex-valued images cannot be directly combined as coherent complex-valued volumetric data and magnitudevalued images were used. Consequently, magnitude-valued images from a k-t SENSE reconstruction were used to simplify the reconstruction methods and focus on the key challenge of locating slices in 3D space and time to reconstruct 4D cine representations of the fetal cardiovascular system.
Motion correction, cardiac synchronization and volumetric reconstruction were all performed on a region of interest (ROI) centred on the fetal heart. In this work static ROIs covering the heart, great vessels and arterial arches were manually prescribed for each slice l. A fetal chest volume of interest and target stack used for motion correction were also user-specified. These were the only manual preparations required. Fetal heart ROIs were prescribed in a few seconds each, while the fetal chest volume of interest was specified in a few minutes.

| Volumetric reconstruction
Reconstruction of a 4D cine from 2D MR images requires information about the spatial and temporal position of the heart in each acquired 2D image frame. While the slice position and acquisition times are known, the relative spatial displacement due to fetal-maternal motion and the cardiac phase are not. However, if the spatial and temporal positions can be accurately estimated, the spatial location of the image frames can be transformed using rigid transformation matrices A = A k k = 1, …, N k , assuming rigid-body displacement, and acquisition times can be mapped to cardiac phases, Image domain volumetric reconstruction is based on forward modelling of the image acquisition process, with volumetric reconstruction formulated as the solution to an inverse problem. In this work the aim was to recover a 4D representation the beating fetal heart. The MR image acquisition model previously used for static anatomy 5 was modified to include a temporal component, assuming the periodic motion of the cardiac anatomy can be characterized by a cardiac phase so that all acquired real-time images can be combined as a single cardiac cycle. This MR image acquisition model describes the relationship between MR image frames, Y k , and a high-resolution 4D cine, X = X h h = 1, …, N h , with three spatial dimensions and a fourth periodic temporal dimension, where X h has elements x ih for spatial index i and temporal index h corresponding to cardiac phases = h h = 1, …, N h , and matrices W hk are the product of a spatial weight, M k , and a temporal weight, d hk , so that w hk ij = d hk m k ij . Each row of M k is made of coefficients {m k ij } i = 1, …, N i that relate the spatial locations of 4D cine X and MR image frame Y k , taking in to account spatial blurring and down-sampling, as well as the movement defined by spatial transformation A k . A point spread function (PSF), determined by the MR acquisition, forms the basis for W hk . The spatial in-plane and through-plane PSF were approximated as Gaussian, 5,12 while the temporal PSF for the band-limited k-t SENSE reconstruction was a sinc. 4 An initial estimate of X was obtained by PSF-weighted interpolation of the scattered data. 5,13 The error between acquired image frames and the image frames predicted by the acquisition model (Equation 1), given an estimate of X, was calculated as where y * jk are the elements of acquired image frames intensitymatched to account for variation in signal scaling and differential bias fields. 5 Robust statistics were employed on both a voxel-and image frame-wise basis to reduce the impact of data corrupted by artefact or being misplaced in space or time. The probabilities of a voxel, p voxel jk , or image frame, p frame k , being an inlier rather than outlier were combined as p jk = p voxel jk p frame k .
PSF-weighted interpolation leads to blurring in the reconstructed volumetric data due to the thick slices of the acquired images. However, super-resolution (SR) methods 5,6 can be used to recover high-resolution spatial data. This was accomplished using gradient descent to minimize the error in Equation (2) by solving which includes a spatial edge-preserving regularization term, R(X), to stabilize the reconstruction and a regularization controlling parameter, λ. Intensity matching, bias correction, robust statistics and X were updated at each iteration. 5

| Motion correction
To achieve the spatial coherence required for volumetric reconstruction, rigid body image registration was used to estimate transformations A that align real-time image frames Y with 4D cine X. Registration was performed in three stages to facilitate convergence, as depicted in Figure 1. The temporal mean of all real-time images in a slice, Y l , was used as a static reference free from cardiac pulsation. Use of static images for the initial stages of motion correction allowed for gross spatial alignment, facilitating cardiac synchronization between slices and providing an initial estimate of A prior to motion correction of individual image frames.
Stacks were first aligned by stack-stack registration with a target stack as in other volumetric reconstruction work, 5,12,13 using static temporal mean images Y.
Using the transformations estimated by stack-stack registration, an initial static volume was reconstructed from Y and then each slice Y l was registered to the volume. Interleaved volume reconstruction and slice-volume registration was repeated over multiple iterations to establish slice-wise alignment.
Frame-wise spatial alignment was performed following cardiac synchronization. An initial 4D cine was reconstructed using estimated cardiac phases, θ, and slice-wise transformations, A. Subsequently, frame-volume registration was performed between each real-time image frame, Y k and the 4D cine frame, X h , with matched cardiac phase, i.e., h = k , followed by volumetric reconstruction using the estimated frame-wise transformations. Interleaved 4D reconstruction and frame-volume registration was repeated over multiple iterations.
Displacement was used to assess the change in position of the voxels in Y subject to estimated transformations A. Global displacement, disp(A), was calculated as where dist y jk , A k (y jk ) is the spatial distance between the original position of voxel y jk and the position of y jk transformed by A k . Similarly, deviation from the average slice transformation was used to assess the dispersion of estimated transformations. Global deviation, dev(A), was calculated as where the average slice transformation, A l , is the p frameweighted Frechét mean. 14 Slice-wise deviation, dev(A l ), was calculated from Equation (5) by restricting the summations to image frames k acquired at slice location l.

| Cardiac synchronization
The cardiac phase of each image frame must be known for the data to be combined as a single cardiac cycle. The heart rate was first estimated independently for each slice and then the cardiac cycle was synchronized between slices, as shown in Figure 2.
A constant heart rate was estimated for each slice from the temporal frequency spectra, i.e., the temporal Fourier transform of the real-time image series, by identifying the maxima in the spatial mean of the temporal frequency spectrum in the fetal heart ROI 4 ( Figure 2A). These maxima were not highly sensitive to the ROI definition, however displacement of the heart and variation of the fetal heart rate resulted in diminished peaks ( Figure 2B). Unreliable heart rate estimates were identified as those with peak height or width more than three normal standard deviations from the median, estimated from the median absolute deviation. 15 If a peak could not be identified in the mean temporal frequency spectrum, the heart rate for that slice was instead calculated as the linear interpolation of reliable heart rate estimates for the slices acquired before and after that slice.
The cardiac cycle was synchronized across slices by estimating a temporal offset for each slice that aligns the pulsation of the heart in overlapping slices. To determine the overlap between slices in 3D space, real-time MR images were cast as cines in volume space. Specifically, a 4D cine, X l , was reconstructed for each slice, l, from the real-time images in that slice, Y l . Volume weights, V l , representing the weighted contribution of the image frames acquired at slice used to determine the spatial intersection of slices. The overlap between slice l and slice l ′ ( Figure 2D) was calculated as Application of a temporal offset, offset l , to X l was achieved by applying a cyclic Fourier time shift, offset, i.e., addition of a linear phase in the temporal frequency domain. Aligning the cardiac cycles of all the slices was equivalent to determining { offset l } l = 1, …, N l that maximized the overlap-weighted similarity between all offset l X l . A group-wise optimization was computationally intensive in practice, so a bootstrap approach was adopted instead, where one temporal offset was estimated at each iteration, i.e., where similarity measure ρ is Pearson's correlation and {l �� } is a set of target slices. The slice l ′ with the greatest overlap with all other slices was assigned offset l � = 0, and used as an initial target slice. A different slice with the greatest overlap with l′ was then identified and offset l was estimated by Equation (6) and cardiac phases were adjusted as k = ( k + offset l ) mod 2 . This slice was then added to the set of fixed target slices {l �� } and the process was repeated using the next slice with the greatest overlap with {l �� }. This continued until all offset l were estimated.

| Fetal study
Multi-planar real-time MR imaging was acquired in a cohort of 11 consecutive singleton pregnancies as part of a larger research project when examination time permitted acquisition of a minimum of three stacks with a combined total of at least 30 slices to ensure good coverage of cardiovascular anatomy. The cases in this initial cohort were used for method development and evaluation. Subsequently, data were acquired in a second cohort of 9 fetuses, 8 singleton and 1 twin pregnancy (ID13), at the discretion of the attending cardiologist without a physicist present at the examination to facilitate robust stack prescription. These cases were used to assess the utility of the methods. Fetal subjects ranged from 23 to 33 weeks gestational age. The study was conducted with the approval of the local research ethics committee and all participants gave written informed consent prior to enrolment. Details of all fetal cases are listed in Table 1, comprising 3 healthy cases and 17 with cardiac anomalies.

| Imaging
Stacks of multi-planar real-time images were acquired on a 1.5 T Ingenia MR system (Philips, Netherlands) using a 2D balanced steady-state free precession (bSSFP) sequence. Highlyprecise scan plane prescription was not required, as specific views could be later obtained from the reconstructed 4D data.
Cardiac synchronization in healthy 30 +3 week gestational age fetus. Heart rate estimation: (A) Temporal mean image, Y l , for a slice with little motion (dev (A l ) = 1.5 mm) with ROI over the fetal heart (dotted yellow line), rotated to show the fetus in radiographic orientation, and plot of temporal frequency spectrum (spatial mean in the fetal heart ROI of the temporal Fourier transform of the real-time image series) showing the clear maxima in the range of expected heart rates (grey band, 115-180 bpm). Peak height (light blue vertical line) and width (dark blue horizontal line) are shown. (B) Fetal movement during the acquisition of a different slice from the same stack (dev(A l ) = 9.4 mm) leads to blurring in the temporal mean image, as well as a lower and broader peak in the temporal frequency spectrum. Unreliable heart rates, such as (B), were interpolated from temporally adjacent slices in the same stack. Sliceslice synchronization: (C) Single frame of the cine, X l , corresponding to slice shown in (A) and (D) volume weights, V l , (blue region), as in-plane and orthogonal through-plane views showing intersection with a slice from a different stack (red region). Slice-slice cardiac synchronization was performed by temporally aligning the cardiac cycles between X l constructed from Y l while using volume weights, Vl, to determine the overlap between slices  However, to achieve a good distribution of scan plane orientations, the first three stacks were prescribed roughly transverse, sagittal and coronal to the fetal trunk for all cases in cohort 1, with additional stacks prescribed in scanner transverse, sagittal, or coronal planes. In-plane rotation was used to reduce the size of the field of view in the phase encode direction while avoiding alignment of pulsatile maternal anatomy (e.g., maternal descending aorta) with the fetal heart in the phase encode direction to avoid complications in the k-t SENSE reconstruction.
Sequence parameters were selected to balance competing goals of good signal and high spatio-temporal resolution with full coverage of the fetomaternal anatomy in the field of view, minimal scan time and safety constraints. Scanner operation was constrained to ensure fetal and maternal safety with respect to specific absorption rate (whole body SAR <2.0 W/kg 16 ), peripheral nerve stimulation (low PNS mode), and sound pressure level (SPL <85 dB(A), accounting for >30 dB attenuation in utero 17 ). While necessary, safety constraints placed a limit on scanner performance and, consequently, achievable resolution and image quality.
Multi-planar bSSFP data were collected with regular Cartesian k-t undersampling 18 using the following default parameters: TR/TE 3.8/1.9 ms, flip angle 60 • , FOV 400 × 304 mm, voxel size 2.0 × 2.0 × 6.0 mm, 8 × acceleration, 72 ms temporal resolution, 96 images per slice, slice overlap 2-3 mm. A steady state was established using an α/2-TR/2 preparation pulse 19 and dummy excitations. Coil calibration data were acquired in a scan immediately preceding the multi-planar acquisition, and low spatial-resolution training data was acquired immediately following the under-sampled data. Acquisition of a typical stack took 155 s. In some cases the default FOV was either decreased or increased in the phase encode direction to accommodate maternal anatomy, with a proportional change in the temporal resolution (cohort 1 median decrease 7.6 ms, median increase 3.8 ms).
Matched ultrasound (US) was acquired on the same day or three days prior to the MRI examination for 7 cases for fetuses in cohort 1 that were part of a multimodality research protocol, as indicated in Table 1. US was performed using an EPIQ V7G system with C9-2 curved array, V6-2 curved volume and X6-1 matrix array transducers (Philips, Netherlands). The echocardiography protocol included comprehensive 2D M-mode imaging as well as B-mode sweeps covering the fetal trunk. B-mode images were combined using the Spatio-Temporal Image Correlation (STIC) method, which uses image-based estimates of the fetal heart rate to produce 4D data. 20

| Reconstruction
All k-t SENSE reconstructions were performed in MATLAB (Mathworks, USA), with additional functionality from ReconFrame 3.0.535 (GyroTools, Switzerland) to place the data in a common 3D space. Slice-slice synchronization (Equation 6) was achieved by constrained non-linear multivariate minimization using the interior point approach implemented in MATLAB's Optimization Toolbox.
Volume reconstruction was performed using the Image Registration Toolkit (BioMedIA, UK), building on the framework implemented by Kuklisova et al 5 with additional functionality for dynamic data. Reconstruction parameters were based on those previously described, 5 and validated for this work using simulation experiments. Volumetric reconstruction was performed at a spatial resolution of 1.25 mm, or 5/8 of the actual acquired in-plane resolution, similar to the ratio used for reconstruction of fetal brain MRI, 5 and a temporal resolution corresponding to N h = 25 cardiac phases, to provide frame-volume registration targets across the entire cardiac cycle. However, the nominal resolution of the reconstructed data was 2 × 2 × 2 mm and 72 ms. A user-specified volume of interest covering the fetal chest was used for slice-volume registration, making use of the anatomy surrounding the heart to facilitate the registration process, while user-specified 2D fetal heart ROIs were combined as a volume and used for frame-volume registration. Reconstruction code can be found at github.com/mriphysics/fetal_cmr_4d (SHA1:b6c94f073c).

| Evaluation
Whole-heart 4D cine data were assessed to establish the quality of the reconstruction method. Three observers (DL, KP, MvP) with 4, 7, and 1 years experience reading cardiac MRI, respectively, were asked to independently navigate the reconstructed 4D cine MR data using the Medical Image Interaction Toolkit Workbench (DKFZ, Germany) and score them in 11 categories, based on the segmental approach to defining cardiac anatomy and pathology. 2,21,22 Data were presented to the readers without context, though the readers may have been present at the time of scanning and may have recalled unique cases. A five-point scale was used for scoring, as follows: 4: high-image quality and distinct appearance of cardiac structures; 3: adequate image quality to determine most details; 2: sufficient image quality to determine some details; 1: poor image quality with significant lack of detail; and 0: inadequate image quality to visualize any cardiac structure.
Ultrasound images acquired during the broader research project were used for comparison with 4D cine data in cases where matched data were available. Left and right ventricular length and width were measured at end-diastole and end-systole in 4-chamber view in the US and MR data by two readers (DL, KP) using the Medical Image Interaction Toolkit Workbench (DKFZ, Germany). Smaller features, such as vessel diameter or septal thicknesses, were not measured as they were depicted by only a few voxels in the real-time MR images. Readers performed measurements on the same 2D M-mode US clip, but were free to select the frames for enddiastolic and end-systolic measurements based on ventricular dilation and contraction. For the 4D data, the plane of the 4-chamber view and frames were determined independently by each reader. In the MRI data, length was measured using the valve annulus as the basal point and width was measured half-way along the length of the ventricle. The MR data was presented to the readers without context and in random order. Measurements were compared using Bland-Altman analysis, paired two-tailed Student's t-test and linear regression to evaluate inter-observer variability and agreement between modalities.

| Simulation
A numerically simulated phantom was used to optimize and evaluate the proposed methods under controlled conditions, as demonstrated in Figure 3. The MRXCAT phantom 23 was used to generate a 4D cine scaled to the size of the fetal heart, χ, with high spatio-temporal resolution (0.44 mm isotropic, 25 cardiac phases) and tissue properties adapted to the in utero environment, i.e., the signal of air in the lungs, trachea and the space surrounding the body in the postnatal simulation were replaced with the signal characteristics of amniotic fluid.
Real-time MR images with 0.5 mm in-plane resolution and 6 mm through-plane resolution were generated from χ using Equation (1) and further down-sampled to 2 mm in-plane resolution by truncating k-space while adding noise. Cardiac phases and displacements were based on in utero measurements of the fetal heart made during the 1062 | VAN AMEROM Et Al. fetal study. Transformation matrices were scaled to change the amount of displacement by a scalar factor 24 to assess varying degrees of motion. Volumetric cine data were reconstructed from simulated real-time MR images and the results were compared with the ground truth data to evaluate the 4D reconstruction framework.

| Velocity-sensitive volume reconstruction
Volume reconstruction from complex-valued real-time images is not a straightforward extension of the proposed framework, but has the potential to generate 4D cine velocity data, i.e., 4D velocity maps, from the same images. Phase contrast bSSFP methods have been proposed, 25 using the inherent velocity-sensitivity of the bSSFP sequence to encode velocity in the phase images. 26 For complex-valued image frames, Ỹ , the phase of an acquired voxel, ∠ỹ jk , is the sum of a velocity-encoded phase, jk , and background phase. The smoothly varying background phase can be estimated by fitting a low-order polynomial to the phase in static tissues and subtracted from ∠ỹ jk so that only jk remains. 27 The jk can then be related to the underlying velocity as where γ is the gyromagnetic ratio, and υ are the velocities and M 1 are the first moments of the gradient waveforms in the readout, phase-encode and slice directions. 25 To reconstruct volumetric velocity maps, the velocityencoded phase would need to be related to three-dimensional space, and the image acquisition model (Equation 1) would have to be adjusted to include velocity terms, precipitating a number of changes to the framework. Instead, as a proof of principle, the framework was modified to reconstruct complex-valued X . Background phase was removed from Ỹ by subtracting a third-order 3D polynomial fit to the phase in static amniotic fluid and tissue in each acquired stack. A phase sign correction was applied so that the signs of all jk were aligned with the target stack, to account for differences in the direction of velocity encoding. All steps in the framework were performed as described previously, with the modification that real-and imaginary-valued 4D cines were also generated separately from the real and imaginary components of Ỹ .

| RESULTS
Simulation experiments were used to evaluate the proposed method and select appropriate reconstruction parameters (Supporting Information Figures S1 and S2). The volume reconstructed from simulated MR images using the established parameters ( Figure 3D) showed cardiac anatomy and pulsation matching the numerically simulated 4D cine, providing evidence of geometric accuracy. However, there was some blurring in all reconstructed 4D cines compared to the high-resolution numerical phantom ( Figure 3B), due to the low spatio-temporal resolution sampling of the real-time images, as well as signal inhomogeneity within the blood pool in the heart.
Fetal results are summarized in Table 2, including reader scores as an indication of image quality. 4D cines were successfully reconstructed in all 11 fetal cases in the first cohort.
Reconstructed 4D cines could be viewed in any arbitrary plane and at any phase in the cardiac cycle to examine cardiovascular morphology and connectivity. An example is shown in Figure 4 where the volume is viewed in variety of double-oblique cardiac planes. The relative size of the chambers of the heart can be seen in horizontal long axis (4chamber view) and mid-short axis views, with dilation and contraction of the left (LV) and right (RV) ventricles and left (LA) and right (RA) atria as the heart beats. Outflow tract and other off-axis views reveal arterial and venous connections of the aorta (Ao), pulmonary artery (PA) and the superior (SVC) and inferior (IVC) vena cava, as well as the positions of the Ao and ductal arch (DA). The Supporting Information includes a video of the reconstructed 4D cine (Supporting Information Video S1) and intermediate results for cardiac synchronization, motion-correction and outlier rejection (Supporting Information Figures S3-S6) for this case. Figure 5 shows a fetus with a cardiac tumor (ID11). The ability to visualize the relationship of the mass and the myocardium in both space and time allowed for qualitative assessment of the impact of the tumour on cardiac function. A particularly large outlier class was estimated in this case (44% of frames rejected) compared to all other cases in cohort 1 (median 10% rejected, range 5-28%); however this particular heart was grossly different in shape and appearance compared to the other fetal cases. It may be that the current implementation of outlier rejection is better suited to fine, high-contrast anatomy, such as the normal fetal brain and heart, than to anatomy containing large regions with homogeneous signal.
The scoring evaluation is summarized in Figure 6 for fetal cases in cohort 1. Scores were generally high in all categories suggesting a potential to assess a range of cardiac anatomy in the 4D cines. Mean scores were between 3 and 4 for all fetal cases (median 3.4) in cohort 1 ( Table 2) indicating that 4D cine MRI data were of adequate quality to determine most anatomical details. The lowest mean scores across all readers were for pulmonary venous connections, atrioventricular connections, and arch anatomy, suggesting that small anatomical features were the most challenging for the readers. The reconstructed 4D cines of the second cohort of fetal subjects were visually assessed by one reader (DL) and determined to be of similar high quality as the first cohort in the majority of cases, while insufficient data resulted in reduced quality in 2 cases for which 1920 (ID16) and 2590 (ID18) image frames were acquired covering the fetal heart, compared to a minimum of 3072 frames for all cases in cohort 1.
All steps of the proposed method were performed as described, with the exception of 2 cases in cohort 1 that required some motion-corrupted data to be manually excluded prior to reconstruction (ID04,ID07). Prior to excluding motioncorrupted slices these 2 cases had dev(A) = 3.8 mm (ID04) and 2.5 mm (ID07), larger than than all other cases in cohort 1 (median 1.3 mm, range 1.0-2.3 mm). Manual intervention resulted in dev(A) = 1.8 mm (ID04) and 2.3 mm (ID07), and yielded visually improved 4D cine reconstructions in both cases, though both received low mean reader scores and had large outlier classes, as expected. Results from motion-corrupted data are shown in Supporting Information Figure S7.
Reconstructed 4D cine MRI was compared with 2D and STIC ultrasound in the cases where matched data was acquired. An example is shown in Figure 7 where the three imaging methods are compared in matched views. Though STIC quality has been reported to be reduced in fetuses in the third trimester, 28 the STIC volume of this 32 +1 week gestational age fetus was the best quality of those acquired. In the other 6 fetal cases with matched ultrasound data, 2D echocardiography images were of comparable quality to those shown. However STIC data in 4 cases were clearly of poor quality, particularly in the through-plane direction, presumably due to fetal motion, while MRI 4D reconstructions were of high quality in all 7 cases with matched US data.
Both 2D US and STIC showed clear definition of vessel and chamber boundaries, including valves, with high spatial and temporal resolution. Reconstructed cine MRI had good contrast between blood and surrounding tissue. However, visualization of valves and other small and rapidly moving anatomy was limited due to the spatio-temporal resolution of the acquisition. While the real-time aspect of 2D US made it robust to fetal motion, interpretation of anatomy was limited to the views acquired at the time of examination. By contrast, both STIC and 4D cine MRI allowed for offline analysis in arbitrary planes with full coverage of the heart and great vessels, facilitating understanding of the spatial relationships between cardiovascular structures.
Analyses of cardiac dimensions measured on US and MR are shown in Supporting Information Figure S8. There was good agreement between 4D MR and 2D US measurements for both readers (p > 0.10), with relatively small difference (bias < 5% of the mean distance) but large 95% limits of agreement (LOA ≈ 5 mm, or 35% of the mean distance), comparable to those observed in other studies of fetal MR and US. 8 STIC data was not measured due to the poor quality in 4 of 7 cases. There was only moderate agreement between readers for both 2D US (p < 0.01, bias = −1.1 mm, LOA = 3.3 mm) and 4D MR (p ≈ 0.05, bias = −0.3 mm, LOA = 2.0 mm), which may also contribute to the large limits of agreement observed in the comparison of 4D MRI and 2D US. The variability of measurements on 2D US made by the two readers is in line with large inter-observer errors in linear measurements in repeatably studies of fetal echocardiography, 29 but may also be due to differences between diastolic/systolic frames selected. The variability of measurements made on 4D MRI may relate to differences in the plane used for the 4-chamber view by the two readers. Note: ID, fetal case number; Reader Score, mean of scores assigned across all categories by all reviewers (Cohort 1) or reviewer 1 only (Cohort 2); Heart Rate, mean and standard deviation of estimated fetal heart rates; Transformations, global displacement, disp(A), and deviation from average slice transformation, dev(A); Outliers, percent of voxel, p voxel , and image frame, p frame , probabilities below 0.5; No. Frames, number of image frames contributing to volume reconstruction, prior to outlier rejection. Images at the edges of the volume, i.e., those contributing fewer than median(N j )∕10 voxels, were excluded from these summary statistics as they had minimal impact on the reconstructed volume but were more likely to be misaligned in space or time; Cohort 1, consecutive singletons scanned for method development and evaluation; Cohort 2, fetuses scanned for clinical assessment; Median, Cohort 1 median values.
A complex-valued 4D cine is shown in Figure 8 with the velocity-encoded phase most sensitive to flow in the fetal superior-inferior direction and the magnitude-valued 4D cine shown as an anatomical reference. High flow can be seen in the great vessels in systole, with reduced flow in diastole.

| DISCUSSION
Whole-heart 4D reconstruction of the fetal heart from multiplanar real-time MRI was successful using the proposed framework in all fetal cases with at least 3000 image frames covering the heart. These reconstructed 4D cines could be visualized in any 2D plane without the need for highly-specific scan plane prescription prior to acquisition or maternal breath hold to minimize motion. The three-dimensional spatial aspect of these 4D cines allowed for interpretation of the connections of the cardiac chambers and vessels, while the temporal aspect improved the depiction of pulsatile anatomy.
Reconstructions of simulated MR images confirmed that spatial and temporal features could be reliably recovered, within the limits of the resolution of the acquisition. There was some blurring in the 4D cine data reconstructed from simulated MR images compared to the high-resolution ground truth, which can be attributed to the low spatio-temporal resolution sampling of the real-time data. There was also a ringing-like signal intensity effect in the 4D cine data that was apparent in the simulation results ( Figure 3D), but had a more subtle impact on the fetal results. This effect may arise from the approximated in-plane spatial PSF, however testing is not straightforward due to the complexity of implementing and computing a sinc PSF.
Overall high scores were assigned by three readers in an evaluation based on the segmental approach to defining cardiac anatomy and pathology. The lowest scores were assigned in categories focused on small anatomical features, such as pulmonary venous connections and arch anatomy, where spatio-temporal resolution limited the depiction of fine details, as was also the case for small septal defects and all atrioventricular and ventriculo-arterial valves. Though a category such as ventricular morphology received high scores in this evaluation, as there was enough other information to give the readers high confidence, particular abnormalities, such as valvular defects, were more visible in ultrasound images. However, some cardiac anatomy was obscured by artefact in the ultrasound images due to acoustic shadowing and, in the case of STIC data, blurring and misalignment. There was good agreement between two-point measurements of ventricular lengths and widths on 2D M-mode US and reconstructed 4D cine MRI, with little bias but large limits of agreement, and moderate inter-observer variability for both modalities. However, the accuracy of MR measurements may be limited by spatiotemporal resolution.
The results of the evaluation suggest that there is potential utility of the method to generate 4D cines that can be used for a comprehensive assessment of the fetal heart, either as an adjunct to ultrasound or in combination with other MRI techniques. For example, while aortic and ductal arch branching patterns may be difficult to determine in the 4D reconstruction, a static volume reconstruction of T 2 -weighted spin echo MRI 30 may more clearly depict the F I G U R E 5 Reconstructed 4D cine of the heart of a 33 +3 week gestational age fetus (ID11) with large fibroma, measuring 24 × 24 × 37 mm, involving the lateral wall of the right ventricle. A volume surface rendering is shown in (A) superior and (B) anterior views, for reference of relative shape and position of the cardiac blood pool (blue surface), fibroma (lined mesh) and pericardial effusion (dotted mesh). The reconstructed 4D cine is shown in (C) 4-chamber, (D) short axis, (E) right ventricular outflow tract and (F) right 3-chamber views. Extensive associated pericardial effusion (PE) can be seen as bright signal surrounding the heart. The fibroma (asterisk) is encapsulated within myocardial tissue and confined to the margin of the right ventricular (RV) cavity, causing external compression of the right side of the heart, including the RV, right atrium (RA) and right ventricular outflow tract (RVOT), while the left ventricle (LV), aorta (Ao), and both the superior (SVC) and inferior (IVC) vena cava appear normal. The blood pool volume rendering in (A) and (B) may be smaller than the true blood pool. All boxes bounding the re-sliced views in (C-F) measure 80 × 80 mm and the fetal heart are shown in radiological orientation, i.e., image axes towards the top and right of the page correspond to left, anterior and superior anatomical directions. Views are shown using spatial B-spline interpolation to avoid voxel distortions details of the arch anatomy. However, 4D reconstruction captures cardiac movement across the entire cardiac cycle, which may improve visualization of anatomy affected by cardiac motion. In the future, the use of higher acquisition under-sampling factors could be explored, possibly in combination with the use of flexible receive coil arrays, 31 to determine how well the reconstruction copes with increased x-f aliasing. These improvements will allow for higher resolution and strengthen the role of MRI in the assessment of the fetal heart in utero. The proposed framework is fully automated aside from anatomical ROIs and identification of a target for initial stack-stack registration. While these user-interactions are not particularly laborious, they could be replaced with automated techniques for segmentation 32,33 and target stack identification using a measure of slice alignment. 34 Manual intervention was required to exclude data with significant motion in 2 cases in cohort 1. In the original 4D cines reconstructed using these motion-corrupted data, the anatomy was blurred and inconsistent as errors in stack-and slice-wise motion-correction lead to slice-slice cardiac synchronization errors, and the reconstruction did not converge to a good depiction of the underlying fetal heart. A motion metric 34 could be used in a preprocessing stage to identify slices with significant motion prior to any motion correction or cardiac synchronization without need for intervention. The focus of this work was on volumetric depiction of the whole fetal heart, and, consequently, the acquisition and reconstruction of k-t SENSE MRI was not fundamentally changed from the 2D implementation 4 and many of the potential improvements previously suggested still apply.
The use of a constant heart rate for each slice does not accommodate beat-to-beat variation resulting in small timing errors. 4 Cardiac monitoring using doppler ultrasound gating 35 or image-based self-gating 33,36,37 could reduce these timing errors, either in combination with or instead of the current approach.
Simulation experiments showed slice-slice cardiac synchronization error of less than 5% of the cardiac cycle for the range of displacements measured in the majority of fetal cases ( Figure S2D), suggesting the proposed method provided reliable results. Cardiac synchronization was performed prior to volumetric reconstruction in the proposed framework and synchronization errors may be improved if the two steps were interleaved in a similar manner to framevolume registration.

F I G U R E 6
Reader score by evaluation category for cardiac 4D cines of 11 fetuses in cohort 1. Scores are shown for reader 1 (left-pointing triangles), reader 2 (squares), and reader 3 (right-pointing triangles), with 4, 7, and 1 years experience reading cardiac MRI, respectively. Median reader scores are shown as large grey markers while colored markers denote scores for individual fetal cases Outlier rejection was performed as implemented by Kuklisova et al, 5 with robust statistics shown to effectively reduce the influence of inconsistent data. As magnitudevalued images were used, sensitivity to image artefacts that manifest in the phase of complex-valued images was reduced compared to the complex-valued outlier rejection used in the 2D framework. 4 However, this did not appear to have a dramatic effect as real-time images were reconstructed using spatially uniform k-t SENSE regularization, thereby suppressing image artefact to some extent. Complex-valued outlier rejection should be possible in a velocity-sensitive reconstruction framework. F I G U R E 7 Comparison of MRI and ultrasound (US) in a 32 +1 week gestational age fetus (ID01) with moderate dilation of the ascending aortic root and rightward tortuosity of the proximal ascending aorta. The fetal heart is shown in matched views in 2D echocardiogram (left column), 4D spatio-temporal image correlation (STIC) ultrasound (centre column) and reconstructed 4D cine MRI (right column). Four-chamber views at

| CONCLUSION
In summary, four-dimensional representation of the fetal heart and great vessels was achieved using a highlyaccelerated multi-planar real-time MRI acquisition combined with retrospective motion correction, cardiac synchronization, outlier rejection and volumetric cine reconstruction in the image domain. The motion-tolerant framework did not require maternal breath-hold or precise scan planning during acquisition, and reconstruction was fully automated aside from user-specified fetal heart and chest ROIs. The framework proved to be robust when used on fetal data and successfully generated good quality 4D cines in all fetal cases where sufficient data was collected covering the entire fetal heart. Readers had an overall high confidence in a comprehensive evaluation of the fetal heart using reconstructed 4D cine MRI and there was good agreement, but large limits of agreement, between two-point measurements of ventricular dimensions on 2D M-mode ultrasound and reconstructed 4D cine MRI. Reconstructions from simulated MR images confirmed that correct spatial and temporal features could be reliably recovered, though there was some blurring due to the spatiotemporal resolution of the acquired real-time images. The use of image domain motion correction methods conferred significant advantages in terms of outlier rejection and the ability to accommodate quite large fetal motion, but had the downside of presenting major challenges for increasing the spatial and temporal resolution of the real-time MR acquisition.
The proposed methods show promise as a framework for motion-corrected reconstruction and 4D assessment of the fetal heart and great vessels. Promising results from a preliminary assessment of velocity-sensitive volume reconstruction suggest potential for simultaneous reconstruction of a fourdimensional cine and fully-encoded velocity information from multi-planar real-time bSSFP MRI. F I G U R E 8 Pseudo 4D velocity mapping in a healthy 30 +3 week gestational age fetus (ID02). Complex-valued cine volume magnitude, |X|, and velocity-encoded phase, ∠X, are shown in (A)-(C) sagittal aortic arch plane, (D)-(F) 3-vessel view, and (G)-(I) coronal plane perpendicular to 3-vessel view. The phase of the reconstructed 4D cine was most sensitive to velocities in the fetal superior-inferior direction, perpendicular to the target stack acquired transverse to the fetal trunk. Flow in the inferior direction (blue), is seen in the descending aorta (DAo) and superior vena cava (SVC), while flow in the superior direction (red) is seen in the ascending aorta (Ao) and pulmonary artery (PA). Peak flows can be seen in systole, sys , with reduced flow during diastole, dia . The right atrium (RA) is labelled in (G) for reference. The fetal heart is oriented in radiological convention with double-oblique planes shown using spatial B-spline interpolation to avoid voxel distortion