Implicit Neural Networks with Fourier-Feature Inputs for Free-breathing
  Cardiac MRI Reconstruction
Cardiac magnetic resonance imaging (MRI) requires reconstructing a real-time
video of a beating heart from continuous highly under-sampled measurements.
This task is challenging since the object to be reconstructed (the heart) is
continuously changing during signal acquisition. In this paper, we propose a
reconstruction approach based on representing the beating heart with an
implicit neural network and fitting the network so that the representation of
the heart is consistent with the measurements. The network in the form of a
multi-layer perceptron with Fourier-feature inputs acts as an effective signal
prior and enables adjusting the regularization strength in both the spatial and
temporal dimensions of the signal. We study the proposed approach for 2D
free-breathing cardiac real-time MRI in different operating regimes, i.e., for
different image resolutions, slice thicknesses, and acquisition lengths. Our
method achieves reconstruction quality on par with or slightly better than
state-of-the-art untrained convolutional neural networks and superior image
quality compared to a recent method that fits an implicit representation
directly to Fourier-domain measurements. However, this comes at a relatively
high computational cost. Our approach does not require any additional patient
data or biosensors including electrocardiography, making it potentially
applicable in a wide range of clinical scenarios.
Johannes F. Kunz, Stefan Ruschke, Reinhard Heckel
# Implicit Neural Networks with Fourier-Feature Inputs for Free-breathing Cardiac MRI Reconstruction

###### Abstract

In this paper, we propose an approach for cardiac magnetic resonance imaging (MRI), which aims to reconstruct a real-time video of a beating heart from continuous highly under-sampled measurements. This task is challenging since the object to be reconstructed (the heart) is continuously changing during signal acquisition. To address this challenge, we represent the beating heart with an implicit neural network and fit the network so that the representation of the heart is consistent with the measurements. The network in the form of a multi-layer perceptron with Fourier-feature inputs acts as an effective signal prior and enables adjusting the regularization strength in both the spatial and temporal dimensions of the signal. We examine the proposed approach for 2D free-breathing cardiac real-time MRI in different operating regimes, i.e., for different image resolutions, slice thicknesses, and acquisition lengths. Our method achieves reconstruction quality on par with or slightly better than state-of-the-art untrained convolutional neural networks and superior image quality compared to a recent method that fits an implicit representation directly to Fourier-domain measurements. However, this comes at a higher computational cost. Our approach does not require any additional patient data or biosensors including electrocardiography, making it potentially applicable in a wide range of clinical scenarios.

## 1 Introduction

Real-time magnetic resonance imaging (MRI) is a dynamic imaging technique for assessing the anatomic structure and function of moving organs in the human body. The reconstruction of cardiac function is challenging due to the relatively fast organ movements compared to the achievable encoding speed. In this paper, we consider the reconstruction of free-breathing cardiac real-time MRI data, where the goal is to reconstruct a video of the beating heart from highly under-sampled data continuously acquired over multiple cardiac and respiratory cycles. The video depicts the cardiac function in real-time meaning that arrhythmia and unprecedented motion can be imaged without assuming periodicity of the motion. This task is especially difficult when only very few samples are collected at any given time instance during the cardiac and respiratory cycles.

Various reconstruction methods for cardiac MRI are based on gated or triggered data [1, 2]. Such methods bin the measurements according to different times in the cardiac and respiratory cycle. The bins are then reconstructed independently using methods for static MRI, or they are reconstructed jointly using signal priors for regularizing among bins [1, 2, 3]. However, binning can lead to artifacts since the heart may not be exactly in the same position at a given point in the cardiac and respiratory cycle, and periodicity of the motion needs to be assumed.
Another class of methods does not rely on gating or triggering, but reconstructs images from continuously acquired un-gated data by regularizing in the spatial and temporal dimensions with hand-crafted signal priors, for example with total-variation norm penalties [14] and through imposing low-rank priors [15, 16, 17, 18, 19]. Thereby, the heart can be imaged in real-time without binning-related inaccuracies.

With the advent of deep learning, signal priors can be learned from data. Supervised methods have been proposed for dynamic MRI that are trained on pairs of measurement data and corresponding ground-truth images [15, 16, 17, 18, 19]. However, ground-truth fully sampled images are generally not available for real-time MRI without assuming periodicity of the motion.

A recent elegant approach to avoid the use of training data but benefit from the prior inherent in neural networks are untrained networks that regularize by fitting a neural network to the measurement data. Untrained neural networks are typically based on convolutional neural networks (CNNs) [19, 16, 17, 18], since CNNs have a provable bias towards fitting smooth signals [13, 14].

In this paper, we propose an untrained approach for cardiac real-time MRI that is based on representing the beating heart with a multi-layer perceptron (MLP) with Fourier-feature inputs and fitting the network so that the representation is consistent with the measurements. Using a Fourier-feature MLP is critical for our approach since it enables us to impose a strong spatial and temporal prior on the heart to be reconstructed. Our contributions are as follows:

* We propose an approach for cardiac real-time MRI based on parameterizing the beating heart with an MLP with Fourier-feature inputs. The method allows controlling the spatial and temporal regularization strength through the setup of the Fourier-feature inputs.
* We compare our approach to the time-dependent deep image prior (t-DIP) [19], a state-of-the-art untrained method based on CNNs. We study the reconstruction quality and computational cost on various datasets in different operating regimes, including different image resolutions, slice thicknesses, and acquisition lengths, and find our method performs on par in terms of image quality, sometimes even marginally better, but at a higher computational cost.
* Compared to two very recently proposed approaches [19, 14] also based on implicit neural representations, we achieve better image quality at a higher computational cost. Specifically, Huang et al. [19] enable cardiac MRI by fitting a Fourier-feature MLP directly to the measurements in the Fourier domain. Fitting the MLP in the Fourier domain has the advantage of being faster compared to fitting it in the image domain. However, this comes at a significant loss in image quality, as we show later, since Fourier-feature MLPs are not a good model for representing an image in the Fourier domain. Specifically, an image is relatively smooth, whereas the spectrum of an image is not smooth, and a Fourier-feature MLP has a bias towards fitting smooth signals. Feng et al. [14]'s method is based on fitting an MLP with hash encoding to the data. Using a hash encoding also has computational benefits over our approach, but a hash encoding is not an effective image prior, thus the method requires additional explicit regularization.
Related Work

In this section, we briefly discuss classical approaches to cardiac real-time MRI, implicit neural networks and CNN-based neural networks for medical image reconstruction.

Approaches to cardiac MRI with sparse and low-rank regularization techniquesA variety of classical approaches to cardiac MRI are based on regularization in the spatial and temporal dimensions with hand-crafted signal priors. For example, Feng et al. [14] impose a total-variation norm penalty over time and Feng et al. [14] impose a total-variation penalty over cardiac phases.

Another line of work incorporates temporal relations by assuming low-rank structures within the temporal series of images, through low-rank tensor decompositions [15], through decompositions into image patches of low-rank [16], and through multi-scale low-rank decompositions [21].

Besides these linear signal models, methods have been developed that assume low-rankness in a kernel space [17, 18] or assume that the images lie on a smooth low-dimensional manifold [19, 1].

Implicit neural networksOur work relies on representing time-varying signals with an implicit neural network, i.e., a function parameterized by a neural network that maps a spatial coordinate vector and a time variable to a pixel value. Implicit neural networks are successfully used for representing images and volumes in a variety of applications. The architecture of the network is important, for images and volumes transforming the coordinates with a Fourier-feature map before passing them through a multi-layer perceptron works well, and most current architectures including NeRF [19], SIREN [18], and MLPs with Fourier-feature inputs [16] use this architecture.

Implicit networks for medical imagingImplicit neural networks are excellent image models and can therefore be used as an image prior for image reconstruction tasks, as demonstrated by Tancik et al. [16] on a toy dataset. Subsequently, Shen et al. [20] have improved the reconstruction quality for static CT and MRI imaging by pre-training the network on a previously reconstructed image of the same subject. Shen et al. [20] also reconstructed a temporal series of 3D images, and used the first image for pre-training the network. Contrary to our approach, the method does not model temporal variations.

Implicit networks have also been used for super-resolution tasks in medical imaging, for example for obtaining high-resolution 3D MRI [22] from low-resolution 2D images, and for building a scale-agnostic model for MRI [23].

Implicit networks for dynamic MRI reconstructionAs mentioned before, most related to our work are two recent methods using implicit networks for dynamic MRI reconstruction. Huang et al. [11] introduced the method NIK which uses a SIREN network [18] with additional explicit regularization to represent the signal in k-space. It outputs a k-space value given the k-space coordinate, the time, and the measurement coil index. Once the network has been trained on the measured frequencies, it can be evaluated at the missing k-space frequencies to yield a reconstructed k-space. We compare our approach to a variant of NIK called k-space Fourier-feature MLP (KFMLP) that we introduce in the methodology section.
Feng et al. [14] learn a neural representation of the dynamic object in the image domain. The implicit network is implemented by an MLP with a hash-encoding for the input coordinates. The MLP is small in size compared to our methods and can be evaluated very efficiently. However, the method relies on additional explicit regularization of the images' temporal total variation and an additional loss term to enforce low-rankness of the solution. Thus, the method relies on explicit regularization and not on an implicit bias inherent to the network architecture as we do.

Untrained methods based on CNNsOur approach relies on regularization with a neural network in the form of an MLP with Fourier-feature inputs. We find that the network has a regularizing effect that is very similar to state-of-the-art untrained CNNs.

Untrained neural networks enable image reconstruction by fitting a neural network to the measurement data, as observed by Ulyanov et al. in the seminal DIP paper [13]. Very simple CNNs enable excellent signal reconstruction performance without any training since CNNs can be excellent image priors [12], and provably reconstruct smooth signals [15, 16]. Untrained networks can outperform classical sparsity based reconstruction for accelerated MRI [17].

The method time-dependent deep image prior (t-DIP) [21] extends the deep image prior [13] to dynamic MRI. Like our approach, the t-DIP fits a single model to the entire series of images, and time is used as an input coordinate. In contrast to using a CNN to regularize in the spatial dimension, we model both the temporal and the spatial dimensions through an implicit network.

The method Gen-SToRM [18] trains a CNN to generate a temporal series of images from a temporal series of low-dimensional network inputs. The network inputs are trained jointly with the CNN. The temporal derivative of the network inputs and the CNNs Jacobian are regularized to favour reconstructions that are on a smooth low-dimensional manifold in image space and that vary smoothly in time. A conceptually similar method has been applied for video reconstruction [1] aside from the context of MRI reconstruction.

## 3 Problem Formulation

Our goal is to reconstruct a video of a moving object from sequentially acquired under-sampled linear measurements. We consider an object \(\mathbf{x}(\sigma)\in\mathbb{R}^{w\times h}\) in the form of an image with width \(w\) and height \(h\) that is parameterized by a motion state \(\sigma\). In our setup, the object \(\mathbf{x}(\sigma)\) is the image of a heart at a certain state in the cardiac and respiratory cycle.

We consider the standard MRI measurement model, where we collect noisy linear measurements

\[\mathbf{y}_{\tau,c}=\mathbf{M}_{\tau}\mathbf{FS}_{c}\mathbf{x}(\sigma_{\tau})+\mathbf{n}_{\tau,c} \tag{1}\]

at the receiver coils \(c=1,\dots,C\). Here, \(\mathbf{F}\) is the 2-D Fourier matrix, \(\mathbf{S}_{c}\) is a diagonal matrix containing the sensitivities of receiver coil \(c\), and \(\mathbf{M}_{\tau}\) is a binary sampling mask encoding the frequencies collected during time \(\tau\). In this work, we consider Cartesian sampling patterns. The measurements are distorted by additive noise \(\mathbf{n}_{\tau,c}\).

Note that the object \(\mathbf{x}(\sigma)\) is changing continuously. The measurement model (1) assumes that the object is essentially constant during a very short time frame which we index by \(\tau\). During this short time frame, only very few measurements (frequencies as selected by the mask \(\mathbf{M}_{\tau}\)) are collected.
The measurements from different coils can be stacked into a single measurement vector

\[\mathbf{y}_{\tau}=\mathbf{A}_{\tau}\mathbf{x}(\sigma_{\tau})+\mathbf{n}_{\tau},\]

where \(\mathbf{y}_{\tau}=[\mathbf{y}_{\tau,1}^{T},\ldots,\mathbf{y}_{\tau,C}^{T}]^{T}\) and \(\mathbf{n}_{\tau}=[\mathbf{n}_{\tau,1}^{T},\ldots,\mathbf{n}_{\tau,C}^{T}]^{T}\). The matrix \(\mathbf{A}_{\tau}\) is the resulting forward map for frame \(\tau\).

Our goal is to reconstruct the image series (or video) \(\mathbf{x}(\sigma_{\tau}),\,\tau=1,\ldots,T\) from the noisy measurements \(\mathbf{y}_{\tau}\). During each time frame indexed by \(\tau\), only very few measurements are collected. Therefore, reconstructing each frame individually results in poor image quality, even if we take prior information about the image into account. Successful reconstruction relies on using prior information about the images in spatial directions as well as in the temporal direction.

Note that in our problem formulation, the index \(\tau\) refers to a short time frame that is used for binning the measurement data such that motion can be neglected. This time frame may not necessarily correspond to the time frame that is used in the video that is ultimately rendered. Thus, the generation of the final video might require re-sampling at a different frame rate. An advantage of our method based on implicit representation is that we can use different frame rates for reconstruction and for the output showed to the end-user.

## 4 Methods

Our FMLP method is based on parameterizing the object \(\mathbf{x}(\sigma)\) with a Fourier-feature MLP and fitting the neural network to the measurement data. Once fitted, we query the neural network to generate a video as an output. The Fourier-feature MLP we use is a good prior for a (slowly) time-varying smooth object, and we use no other regularization for reconstruction.

We parameterize the object \(\mathbf{x}(\sigma)\) by a neural network \(f_{\mathbf{\theta}}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{2}\) mapping a spatio-temporal coordinate vector \(\mathbf{c}=[x,y,t]^{T}\) to a complex image intensity value, represented by a real and an imaginary part. We focus on 2D images, but the method can easily be extended to 3D by adding another spatial coordinate so that \(\mathbf{c}=[x,y,z,t]\). The coordinates \(x\) and \(y\) denote physical locations within the field of view measured in meters, and the time-coordinate \(t\) is time measured in seconds. The origins of the spatial and temporal coordinate axes can be set arbitrarily as the particular neural network that we use is shift-invariant. We define \(t=0\,\mathrm{s}\) as the beginning of the measurement process. For computing a pixel-based image \(\mathbf{f}_{\mathbf{\theta}}(t_{\tau})\) at time \(t_{\tau}\), the network \(f_{\mathbf{\theta}}\) is evaluated on a regular grid of spatial coordinates within the field of view with the temporal coordinate fixed to \(t_{\tau}\)

Figure 1: The FMLPs’ network consists of separate spatial and temporal Fourier-feature embeddings that are concatenated, followed by an MLP that outputs the complex image intensity at the specified coordinate.
We fit a network by minimizing the reconstruction loss defined as

\[\mathcal{L}(\mathbf{\theta})=\frac{1}{T}\sum_{\tau=1}^{T}\left\|\mathbf{A}_{\tau}\mathbf{f}_{ \mathbf{\theta}}(t_{\tau})-\mathbf{y}_{\tau}\right\|_{2}^{2}.\]

The time \(t_{\tau}\) is chosen to be at the center of the acquisition time window of frame \(\tau\) and the spatial coordinate grid used for computing \(\mathbf{f}_{\mathbf{\theta}}(t_{\tau})\) is matched to the Cartesian sampling grid in k-space. After fitting, reconstructed images can be sampled from the implicit representation as \(\mathbf{\hat{x}}(\sigma_{\tau})=\mathbf{f}_{\mathbf{\theta}}(t_{\tau})\).

The network architecture is critical for performance. Our architecture is based on using Fourier-features as inputs to a ReLU-MLP and is depicted in Figure 1. The Fourier-feature inputs are used in a variety of implicit neural network architectures [11, 12, 13]. The Fourier-feature layer embeds the spatial coordinates and the temporal coordinate separately as \(\mathbf{\gamma}(\mathbf{c})=[\mathbf{\gamma}_{\text{spatial}}([x,y])\), \(\mathbf{\gamma}_{\text{temporal}}(t)]^{T}\) with

\[\mathbf{\gamma}_{\text{spatial}}([x,y]) =\left[\sin\left(\mathbf{B}\begin{bmatrix}s_{\text{x}}x\\ s_{\text{y}}y\end{bmatrix}\right),\,\cos\left(\mathbf{B}\begin{bmatrix}s_{\text{ x}}x\\ s_{\text{y}}y\end{bmatrix}\right)\right]^{T},\] \[\mathbf{\gamma}_{\text{temporal}}(t) =\left[\sin(\mathbf{b}s_{\text{t}}t),\,\cos(\mathbf{b}s_{\text{t}}t) \right]^{T}.\]

Here, the \(\sin\) and \(\cos\) functions are applied element-wise. The coordinates are scaled by hyper-parameters \(s_{\text{x}}\), \(s_{\text{y}}\), and \(s_{\text{t}}\) with units \(\frac{1}{\text{m}}\) and \(\frac{1}{\text{s}}\), respectively, which control the variance of the angular frequencies in the spatial and temporal dimensions. The elements in the matrix \(\mathbf{B}\) and the vector \(\mathbf{b}\) are drawn independently from a zero-mean Gaussian distribution with unit variance \(\mathcal{N}\left(0,\,1\right)\) and are fixed and not optimized over. We choose \(\mathbf{B}\in\mathbb{R}^{256\times 2}\) and \(\mathbf{b}\in\mathbb{R}^{64}\) for a total feature size of \(512+128\).

We observed that embedding the spatial and temporal coordinates separately leads to sharper reconstructions compared to the intuitive setup of having one Fourier-feature embedding for all three coordinates together, see Section A.4 in the supplementary material.

The MLP takes the Fourier-feature vector \(\mathbf{\gamma}(\mathbf{c})\) as input. After each linear hidden layer, the ReLU activation is applied and the output feature vector is normalized to zero mean and unit variance. The output layer has two neurons with a linear activation function and without normalization, which output the real and imaginary part of the complex image intensity value, respectively. The weights of the linear layers are initialized by drawing from \(\mathcal{N}\left(0,\,\sigma_{\text{linear}}^{2}\right)\), where the variance \(\sigma_{\text{linear}}^{2}\) can be tuned. The final output is multiplied by a constant \(s_{\text{out}}\) that is tuned as a hyperparameter depending on the scaling of the k-space data and the sensitivity maps. In our implementation, we set \(\sigma_{\text{linear}}=0.01\) and the MLP consists of 7 hidden layers with 512 neurons each if not explicitly stated otherwise.

### Baseline methods

In our experimental results, we consider two method described in more detail below for comparison. The first method is a variant of Huang et al. [14]'s NIK method that is based on fitting a Fourier-feature MLP in the Fourier domain. We consider a variant of the method and not the original one, since code for the original is not available, our sampling setup is different, and we wanted to get best performance when using a Fourier-feature MLP in the Fourier domain, which
required us to deviate from NIK. The second method we compare to is the time-dependent DIP [20], a state-of-the-art untrained method based on CNNs.

#### 4.1.1 Fourier-feature MLP as implicit k-space prior (KFMLP)

The KFMLP is conceptually similar to the NIK [14] that has been proposed recently. It parameterizes the time-varying k-space \(\mathbf{FSx}(\sigma)\) with a neural network. The network \(f_{\mathbf{\theta}}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{2C}\) maps frequency-time coordinates \([k_{\mathrm{x}},k_{\mathrm{y}},t]\) to complex k-space values that are represented by a real and an imaginary part for each of the \(C\) receiver coils. The frequency-coordinates \(k_{\mathrm{x}}\) and \(k_{\mathrm{y}}\) are normalized to \([-\pi,\pi)\) and the time-coordinate \(t\), as for the FMLP, measures time in seconds. The network is fitted to the measured k-space data at frequency-time coordinates along the sampling trajectory. For computational purposes, we segment the trajectory into frames that are treated as samples of the dataset. Specifically, we evaluate the network \(f_{\mathbf{\theta}}\) at the k-space coordinates and sampling times that correspond to the measurements in the vectors \(\mathbf{y}_{\tau}\). The resulting reconstructed k-space along the trajectory of frame \(\tau\) is denoted as \(\mathbf{f}_{\mathbf{\theta},\tau}^{\mathrm{traj}}\). We minimize a similar loss as proposed for the NIK [14],

\[\mathcal{L}(\mathbf{\theta})=\frac{1}{T}\sum_{\tau=1}^{T}\mathcal{L}_{\mathrm{ HDR}}(\mathbf{f}_{\mathbf{\theta},\tau}^{\mathrm{traj}},\mathbf{y}_{\tau})+R(\mathbf{f}_{\mathbf{ \theta},\tau}^{\mathrm{traj}}).\]

The loss is composed by a high dynamic range reconstruction loss and an explicit k-space regularization term. The high dynamic range loss is defined as

\[\mathcal{L}_{\mathrm{HDR}}(\hat{\mathbf{y}}_{\tau},\mathbf{y}_{\tau})=\left\|\frac{ \hat{\mathbf{y}}_{\tau}-\mathbf{y}_{\tau}}{|\operatorname{sg}(\hat{\mathbf{y}}_{\tau})|+ \varepsilon}\right\|_{2}^{2}.\]

Here, the operator \(\operatorname{sg}(\cdot)\) stops the propagation of the gradient during back-propagation, \(|\cdot|\) computes the absolute value accounting for the positive and negative k-space values, and \(\varepsilon>0\) is a hyperparameter that adjusts the compression of the dynamic range. The division is computed element-wise. Note that our definition of the high dynamic range loss differs from the NIK that does not take the absolute value in the denominator.

The k-space regularization term \(R(\mathbf{f}_{\mathbf{\theta},\tau}^{\mathrm{traj}})\) is adopted from the NIK and is defined as

\[R(\mathbf{f}_{\mathbf{\theta},\tau}^{\mathrm{traj}})=\lambda_{\mathrm{denoiser}} \mathcal{L}_{\mathrm{HDR}}(\mathbf{f}_{\mathbf{\theta},\tau}^{\mathrm{traj}},\mathbf{K}_{ \tau}\mathbf{f}_{\mathbf{\theta},\tau}^{\mathrm{traj}}),\]

where \(\mathbf{K}_{\tau}\) is a diagonal matrix with entries \(e^{-d/2\sigma^{2}}\) that weight the k-space values according to their respective coordinate distances \(d=\sqrt{k_{\mathrm{x}}^{2}+k_{\mathrm{y}}^{2}}\) to the k-space center. The hyper-parameters \(\lambda_{\mathrm{denoiser}}>0\) and \(\sigma>0\) adjust the regularization strength. By ablation studies, we find that \(\varepsilon=10^{4}\), \(\sigma=10\) and \(\lambda_{\mathrm{denoiser}}=0.1\) yield good results for our datasets. Note that the high dynamic range reconstruction loss and the k-space regularization term are not essential, but they yield a small improvement compared to using the standard \(l_{2}\) reconstruction loss without any explicit k-space regularization, see Section B.2 and Section B.3 in the supplementary materials.

After training, a reconstructed image of at time \(t_{\tau}\) is obtained by evaluating the network at all frequencies on the Cartesian grid and combining the reconstructed k-spaces via coil combination. The reconstructed k-space on the Cartesian grid, with time coordinate fixed to \(t_{\tau}\), is denoted
by \(\mathbf{f}_{\mathbf{\theta}}^{\text{grid}}(t_{\tau})\) that stacks the reconstructed k-spaces \(\mathbf{f}_{\mathbf{\theta},c}^{\text{grid}}(t_{\tau})\) of the receiver coils \(c=1,\ldots,C\). A coil-combined image is computed as

\[\mathbf{\hat{x}}(\sigma_{\tau})=\sum_{c=1}^{C}\bar{\mathbf{S}}_{c}^{*}\mathbf{F}^{-1}\mathbf{f}_ {\mathbf{\theta},c}^{\text{grid}}(t_{\tau}),\]

where \(\bar{\mathbf{S}}_{c}^{*}\) denotes the complex conjugate of the normalized sensitivity map of coil \(c\) given as

\[\bar{\mathbf{S}}_{c}=\left(\sum_{j=1}^{C}\mathbf{S}_{j}^{*}\mathbf{S}_{j}\right)^{-1}\mathbf{S }_{c}.\]

The KFMLP uses the same network architecture as the FMLP with Fourier-features that are given as input to an MLP with ReLU activation function and normalization. The k-space coordinates are embedded in a separate Fourier-feature vector from the time coordinate. Again, the coordinate scales \(s_{\text{x}}\), \(s_{\text{y}}\), and \(s_{\text{t}}\) are tuned as hyper-parameters.

The KFMLP architecture differs slightly from the NIK architecture that uses a sinusoidal activation function, jointly embeds the time and k-space coordinates, and treats the coils \(c=1,\ldots,C\) as an additional input dimension.

#### 4.1.2 Time-dependent DIP

We implement the time-dependent DIP [20] which is a state-of-the-art method for untrained reconstruction. It uses a CNN with time-varying network inputs as implicit prior. The time-varying network inputs are generated by embedding a temporal input coordinate on a helix trajectory and then mapping the embedding to the CNN's input features using a small fully connected network.

The helix is a trajectory in \(\mathbb{R}^{3}\) with unit radius and slope \(z_{\text{slack}}\) that we tune as a hyperparameter. The angular frequency of the helix introduces prior information about the heart rate and periodicity of motion. For matching the angular frequency of the helix with the heart rate, an estimate of the number of cardiac cycles is required, which we obtain from an electrocardiogram (ECG).

Due to the validation lines that are extracted at random positions, the frames are not equally spaced in time. We account for the deviation by modifying the sampled positions on the helix trajectory. Instead of sampling equally spaced points based on the frame index \(k\) as proposed by the authors, we use the measurement time \(t_{\tau}\) to ensure that the extraction of validation lines does not diminish the performance of the t-DIP. The sampled positions on the helix are mapped to a more expressive latent space by a small MLP called MapNet. The output of the MapNet is reshaped as input for the subsequent CNN.

We adopt the CNN architecture but modify it to fit our image resolution. The proposed architecture of the CNN up-samples the feature resolution in powers of two using nearest-neighbor interpolation. We adjust the input feature resolution and the number of up-sampling layers to obtain an output resolution that is close to our target resolution. The output image is then center-cropped to the target resolution. Between interpolation layers, two convolutional layers are applied. We tune the channel depth of the convolutional layers as a hyperparameter.
Experiments and Results

We assess the reconstruction quality and computational demands of FMLP on real data that we collected as well as on phantom data in the supplementary material. We compare to the KFMLP method and the time-dependent DIP, and find that FMLP performs as well as the t-DIP and significantly better than KFMLP in terms of image quality, but at a higher computational cost than both methods. Code to reproduce the results is available at [https://github.com/MLI-lab/cinemri](https://github.com/MLI-lab/cinemri).

### Datasets

We conduct experiments on datasets that were acquired on a 3 T Edition X scanner (Philips Healthcare, The Netherlands). The scans are taken of a healthy 30 years old male volunteer. This study was approved by the local ethics committee and written informed consent was obtained.

Measurement data was acquired at two different resolutions \(2.27\times 2.26\,\mathrm{mm}^{2}\) (\(264\times 186\) acquisition matrix size) and \(1.25\times 1.26\,\mathrm{mm}^{2}\) (\(480\times 334\)), with \(10\,\mathrm{mm}\) and \(5\,\mathrm{mm}\) slice thickness, respectively. We refer to those as 'low-resolution high-SNR' and 'high-resolution' datasets. A third dataset was acquired at low resolution with a reduced (isotropic) slice thickness of \(2.25\,\mathrm{mm}\), decreasing the signal strength and the signal-to-noise ratio (SNR). This dataset is referred to as the 'low-resolution low-SNR' dataset. For each of those configurations, a breath-hold triggered scan and a free-breathing scan were taken. The acquisition length was set retrospectively for the free-breathing scan and we conduct experiments with \(4.0\,\mathrm{s}\), \(8.0\,\mathrm{s}\), and \(15.9\,\mathrm{s}\) of acquired data.

For the breath-hold scan, the data acquisition is triggered by an ECG such that measurements of similar cardiac phases are binned. The dynamics are then reconstructed with a sparsity-based reconstruction method for static MRI which the MRI scanner performs by default. The breath-hold reconstructions serve as a visual reference for the free-breathing reconstructions.

During the free-breathing scan, measurements are acquired continuously (ungated). We use a partial-Fourier Cartesian sampling pattern, where the k-space is fully sampled along the \(k_{\mathrm{x}}\)-dimension (frequency-encoding) and randomly under-sampled along the \(k_{\mathrm{y}}\)-dimension (phase-encoding). Thus, measurements are taken sequentially along \(k_{\mathrm{y}}\)-lines in the k-space. The measured \(k_{\mathrm{y}}\)-lines are binned into frames retrospectively. Further details on the measurement parameters are listed in Section E in the supplementary materials.

### Performance metrics and quality criteria

It is common practice to measure the image quality of reconstruction methods with full-reference image quality metrics, such as the structural similarity index measure (SSIM) [23] or the visual information fidelity (VIF) [15]. For our synthetically generated phantom dataset in the supplement (see Section C), we measure performance with those metrics.

For our original data (and original data in the context of real-time MRI in general), ground-truth images are not available. We therefore measure performance in terms of an estimate of the MSE as well as through visual comparisons.

Specifically, we visually compare the reconstructed images to the ECG-triggered breath-hold reconstructions obtained from the scanner by binning-based reconstruction methods using sparsity regularization. We focus on anatomic details with the region of interest, the heart, such as papillary muscles in the left ventricle that are rapidly changing shape during systoles.
We estimate a normalized version of the mean squared error based on a hold-out set of randomly subsampled frequencies as follows. We randomly extract \(5\,\%\) of the measured \(k_{\text{y}}\)-lines for validation. The remaining \(k_{\text{y}}\)-lines are binned into frames \(\tau=1,\dots,T\) for training, such that the k-space of every frame contains \(N_{\text{lines}}=6\) many \(k_{\text{y}}\)-lines. Each validation line \(\mathbf{y}_{v},v=1,\dots,V\) is assigned to the frame \(\tau_{v}\) that is closest in sample time. Let \(\hat{\mathbf{y}}_{v}\) be the reconstructed k-space lines predicted by a reconstruction algorithm (FMLP, KFMLP, and t-DIP in our setup). The lines are obtained by first applying the sensitivity maps and Fourier transform to the reconstructed images \(\hat{\mathbf{x}}(\sigma_{\tau_{v}})\) and then extracting the \(k_{\text{y}}\)-line corresponding to \(\mathbf{y}_{v}\). We estimate the signal-to-error ratio (SER) as a metric for quantifying data consistency

\[\text{SER}=10\log_{10}\frac{\sum\limits_{v=1}^{V}\left\lVert\mathbf{y}_{v}\right \rVert_{2}^{2}}{\sum\limits_{v=1}^{V}\left\lVert\hat{\mathbf{y}}_{v}-\mathbf{y}_{v} \right\rVert_{2}^{2}}.\]

### Implementation details

We implemented the FMLP, KFMLP, and t-DIP using PyTorch. The models are optimized using Adam [1] with learning rate \(2\times 10^{-4}\) for the FMLP and the KFMLP, and learning rate \(1\times 10^{-4}\) for the t-DIP. The models are trained until no new SER high-score has been reached for 200 epochs. We evaluate the reconstruction quality at the epoch of maximum SER. The experiments are conducted on a server equipped with a RTX 6000 GPU with \(24\,\mathrm{GB}\) VRAM and an Intel(R) Core(TM) i9-9940X CPU with \(128\,\mathrm{GB}\) available RAM.

### Quantitative reconstruction quality analysis

In this section, we quantitatively compare the reconstruction quality of the FMLP, the KFMLP, and the t-DIP in terms of the estimated normalized MSE (specifically, in terms of the SER defined in the previous section).

Different operating regimesWe start by comparing the image reconstruction quality measured in terms of the SER in different operating regimes that feature different image resolutions and slice thicknesses. Specifically, we fit all three methods to the low-resolution high-SNR dataset, the low-resolution low-SNR dataset, and the high-resolution dataset. The hyper-parameters of the models are tuned on each dataset individually and the methods are trained on \(T=225\) frames.

The results in Table 1 show that the FMLP and the t-DIP achieve a similar SER with marginal differences in all evaluated operating regimes. While the t-DIP achieves a slightly higher SER than the FMLP on the low-resolution high-SNR dataset, the FMLP performs slightly better on the low-resolution low-SNR and the high-resolution datasets. The KFMLP, learning a representation of the k-space, achieves much lower SER scores than the FMLP or the t-DIP. The difference is substantial on all evaluated datasets.

Performance as a function of acquisition lengthNext, we study the reconstruction quality for different measurement acquisition length. We measure the image quality in terms of the SER within the first \(4\,\mathrm{s}\) of reconstructed images, while training the methods on increasing acquisition
lengths, i.e., on 4 s, 8 s, and 16 s of acquired data, respectively. The experiment is conducted on the low-resolution high-SNR dataset.

Figure 2 shows that the SER of the FMLP and the t-DIP improves when training the methods on more data, indicating that both methods not only take advantage of temporal correlations between adjacent frames but also over the entire acquisition length. The KFMLP, by contrast, does not improve substantially. Remarkably, the FMLP and the t-DIP improve by a similar margin even though the t-DIP uses prior information about the heart rate that is encoded by the angular frequency of the helix trajectory. Information about the cardiac phase might be advantageous for correlating frames over multiple cardiac cycles. The FMLP does not incorporate such prior but nevertheless improves at a similar rate as a function of the acquisition length.

### Visual quality assessment

We compare the reconstructed images of the FMLP, the KFMLP, and the t-DIP to the ECG-triggered breath-hold reconstructions that serve as a visual reference. The methods are compared on the low-resolution high-SNR, the low-resolution low-SNR, and the high-resolution dataset. For each

\begin{table}
\begin{tabular}{c|c|c|c} \hline \hline dataset & FMLP & KFMLP & t-DIP \\ \hline low-resolution, high-SNR & 17.16 & 12.53 & **17.19** \\ low-resolution, low-SNR & **10.62** & 8.64 & 10.59 \\ high-resolution & **9.00** & 6.51 & 8.97 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Maximum SER of the models on datasets in different operating regimes. FMLP and t-DIP perform similarly, while KFMLP performs significantly worse.

Figure 2: The validation SER within the first 4 s of acquisition time improves with increasing the amount of measurement data beyond 4 s. The reconstruction quality can be improved by training the model on a longer acquisition time. The models were trained on the low-resolution high-SNR dataset with the same configurations as in Figure 3 and \(z_{\text{slack}}=0.1,0.2\), and 0.4 for \(T=225\) (4 s), 450 (8 s), and 900 (16 s), respectively.
Figure 3: On the low-resolution high-SNR dataset, the image quality of the FMLP and the t-DIP are similar, whereby the FMLP recovers anatomic details such as the papillary muscles in the left ventricle (red arrow) more accurately. The reconstructions by the KFMLP are distorted by aliasing-like artifacts and fine-structured noise such that anatomic details are well not recognizable. The models were trained on \(T=225\) frames. The hyperparameter configuration are \(s_{\rm x}=s_{\rm y}=30\,\rm m^{-1}\), \(s_{\rm t}=1\,\rm s^{-1}\), 7 hidden layers, 512 neurons per hidden layer for the FMLP, \(s_{\rm x}=s_{\rm y}=15\), \(s_{\rm t}=1\,\rm s^{-1}\), 7 hidden layers, 512 neurons per hidden layer for the KFMLP, and 256 channels and \(z_{\rm slack}=0.1\) for the t-DIP.
dataset, we show reconstructions of four frames that have been selected to cover a cardiac cycle.

The reconstructions in the Figures 3, 4, and 5, show that FMLP performs on par with the t-DIP in terms of perceived image quality. Although the overall differences in quality are marginal, the FMLP recovers anatomic details slightly more accurately on the low-resolution high-SNR dataset. The differences are most apparent in small moving anatomic structures of the heart, such as the papillary muscles, see Figure 3. On the low-resolution low-SNR and the high-resolution datasets, the FMLP and the t-DIP exhibit similar artifacts, see Figures 4, and 5.

We find that the FMLP achieves a better image quality compared to the KFMLP on all datasets. On the low-resolution high-SNR dataset, the KFMLP suffers aliasing-like artifacts that blur anatomic details, such as the papillary muscles, see Figure 3. On the low-resolution low-SNR dataset, the aliasing-like artifacts are superimposed by the noise that further degrades the image quality compared to the FMLP (Figure 4). However, the noise of the KFMLP becomes most prominent on the high-resolution dataset, where it obscures any meaningful details in the cardiac

Figure 4: For an isotropic slice thickness (low-resolution low-SNR dataset) the reconstruction quality of all methods degrades and new artifacts are introduced. The FMLP and the t-DIP achieve a similar image quality and outperform the KFMLP. The models were trained on \(T=225\) frames. The hyperparameter configuration are: \(s_{\mathrm{x}}=s_{\mathrm{y}}=30\,\mathrm{m}^{-1}\), \(s_{\mathrm{t}}=2\,\mathrm{s}^{-1}\), 5 hidden layers, and 256 neurons per hidden layer for the FMLP, \(s_{\mathrm{x}}=s_{\mathrm{y}}=15\), \(s_{\mathrm{t}}=1\,\mathrm{s}^{-1}\), 7 hidden layers, and 512 neurons per hidden layer for the KFMLP, and 256 channels and \(z_{\mathrm{slack}}=0.5\) for the t-DIP.
region (Figure 5).

### Computational cost

The exact computational costs depend on the hyperparameter configuration, but for most of our tested hyperparameter configurations, training the FMLP requires more epochs and takes longer to train per epoch than the t-DIP and the KFMLP. The KFMLP and the t-DIP reach the maximum SER in a similar amount of training time, whereas the t-DIP converges in fewer epochs.

The computational performance characteristics of the configurations with the highest SER are listed in Table 2. We report the model sizes, memory requirements during training, training time per epoch, and training time until the maximum SER is reached on the low-resolution high-SNR dataset and on the high-resolution dataset, each with \(T=225\) frames. The total training time is measured by the wall-clock time and includes data-loading, computation of the SER for validation, and logging the training progress. The time per epoch is measured by the wall-clock time over 100 epochs of

Figure 5: On the high-resolution dataset, the FMLP and the t-DIP achieve a similar reconstruction quality and both methods suffer from similar artifacts. The quality of the KFMLP is decreased substantially by noise-like artifacts. The models were trained on \(T=225\) frames. The hyperparameter configuration are: \(s_{\rm x}=s_{\rm y}=30\,\rm m^{-1}\), \(s_{\rm t}=1.0\,\rm s^{-1}\), 7 hidden layers, 512 neurons per hidden layer for the FMLP, \(s_{\rm x}=s_{\rm y}=10\), \(s_{\rm t}=1.0\,\rm s^{-1}\), 7 hidden layers, 512 neurons per hidden layer for the KFMLP, and 256 channels and \(z_{\rm slack}=0.5\) for the t-DIP.
training without additional overhead. The hyperparameter configurations with maximum SER are chosen that are specified in Figure 3 and 5. Note that the t-DIP's CNN uses five up-sampling layers on the low-resolution dataset, whereas six up-sampling layers are required on the high-resolution dataset.

We note that the computational demands required for training scale differently with image resolution for the three methods. An advantage of the FMLP is that the network size can stay constant as the image resolution is increased. However, the network has to be evaluated for all pixels and, thus, the computational demands scale linearly with image resolution. The CNN of the t-DIP, by contrast, needs to be adapted to an increasing image resolution by adding up-sampling layers and convolutional layers. Thus, the number of trainable parameters increases with image resolution. However, doubling the image resolution merely requires adding one more stage of up-sampling and convolutional layers. Thus, the cost of evaluating the CNN scales efficiently with image resolution. The network size of the KFMLP can also be set independently of the image resolution. Moreover, the computational demands do not scale with the image resolution but with the number of measured coordinates. The network is evaluated on the full Cartesian grid only once after training.

## 6 Discussion and Conclusion

We proposed an untrained reconstruction method based on implicit networks, called FMLP, for cardiac real-time MRI. The method uses an ReLU-MLP with Fourier-feature inputs to encode spatio-temporal coordinate vectors for representing an image of the beating heart. We evaluated the method on experimental datasets for 2D free-breathing cardiac real-time MRI covering different operating regimes, i.e., different image resolutions, slice thicknesses, and acquisition lengths.

We find that FMLP can achieve reconstruction quality on par or slightly better with the best CNN-based untrained methods (t-DIP), as measured with quantitative metrics (SER) and as determined through visual comparisons. FMLP tends to improve slightly faster as a function of the acquisition time (see Figure 2), and slightly surpassed t-DIP for longer acquisition times. The FMLP only relies on the implicit neural network for regularization.

We note that the similarity in performance to the t-DIP is somewhat expected since MLPs with Fourier-feature inputs and CNNs both have a bias towards fitting smooth images and low-frequency

\begin{table}
\begin{tabular}{c|c|c|c|c} \hline \hline res. & property & FMLP & KFMLP & t-DIP \\ \hline \multirow{4}{*}{\begin{tabular}{c} FMLP \\ \end{tabular} } & \# parameters & \(1.91\times 10^{6}\) & \(1.93\times 10^{6}\) & \(5.94\times 10^{6}\) \\  & memory on GPU & \(3.55\,\mathrm{GB}\) & \(1.14\,\mathrm{GB}\) & \(1.94\,\mathrm{GB}\) \\  & time per epoch & \(19.92\,\mathrm{s}\) & \(4.06\,\mathrm{s}\) & \(8.56\,\mathrm{s}\) \\  & number of epochs & \(744\) & \(309\) & \(252\) \\  & total training time & \(510\,\mathrm{min}\) & \(35\,\mathrm{min}\) & \(42\,\mathrm{min}\) \\ \hline \multirow{4}{*}{
\begin{tabular}{c} FMLP \\ \end{tabular} } & \# parameters & \(1.91\times 10^{6}\) & \(1.93\times 10^{6}\) & \(7.12\times 10^{6}\) \\  & memory on GPU & \(9.21\,\mathrm{GB}\) & \(1.22\,\mathrm{GB}\) & \(4.03\,\mathrm{GB}\) \\ \cline{1-1}  & time per epoch & \(63.82\,\mathrm{s}\) & \(4.11\,\mathrm{s}\) & \(26.65\,\mathrm{s}\) \\ \cline{1-1}  & number of epochs & \(1872\) & \(601\) & \(176\) \\ \cline{1-1}  & total training time & \(3237\,\mathrm{min}\) & \(110\,\mathrm{min}\) & \(91\,\mathrm{min}\) \\ \hline \hline \end{tabular}
\end{table}
Table 2: Computational performance characteristics.
contents of the signal first [14, 21]. This implicit bias suppresses high-frequency aliasing artifacts, thereby interpolating the missing frequencies in the k-space and reconstructing the image.

While an MLP with Fourier-features is a good image model, it is not a good model for representing an image in the Fourier domain. We found that KFMLP, which relies on representing the signal in the Fourier domain with a Fourier-feature MLP, suffered from aliasing and noise-like artifacts that distorted the images considerably. Although the KFMLP learned to represent the measurements along the sampling trajectory well, it poorly generalized to frequencies that had not been measured. We note that the partial-Fourier random Cartesian sampling pattern used for the experiments may pose an especially difficult challenge to the reconstruction methods as a small subset of k-space lines is measured repetitively. Thus, a large fraction of the k-space is not measured at any time and needs to be interpolated with an efficient image model. Huang et al. [12] tested the NIK (the KFMLP is a variant of the NIK) on retrospective data with a golden-angle radial trajectory. On such data, the method did not suffer from strong artifacts compared to the KFMLP on our datasets.

We hasten to add that the imaging quality comes at a high computational cost. The FMLP is more expensive to train than the t-DIP and the KFMLP. However, we think that future research (such as pre-training the FMLP) can substantially improve the running time of the FMLP, and different feature encoding altogether might also improve the running time.

A major advantage of the FMLP over the t-DIP and other approaches is it's flexibility: To generalized generalize to 3D data, we only need to add an additional spatial input.

## Acknowledgements

The authors would like to thank you Kilian Weiss (Philips GmbH Market DACH, Hamburg, Germany) for helpful discussion. The work is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) 11 - 456465471, 464123524.
