In this chapter,^{1} we present a new method for generating synthetic MRI data. The quantitative validation of reconstruction algorithms requires reliable data. Rasterized simulations are popular but they are tainted by an aliasing component that impacts the assessment of the performance of reconstruction. We introduce analytical simulation tools that are suited to parallel magnetic resonance imaging and allow one to build realistic phantoms. The proposed phantoms are composed of ellipses and regions with piecewisepolynomial boundaries, including spline contours, Bézier contours, and polygons. In addition, they take the channel sensitivity into account, for which we investigate two possible models. Our analytical formulations provide welldefined data in both the spatial and kspace domains. Our main contribution is the closedform determination of the Fourier transforms that are involved. Experiments validate the proposed implementation. In a typical parallel MRI reconstruction experiment, we quantify the bias in the overly optimistic results obtained with rasterized simulations—the inversecrime situation. We provide a package that implements the different simulations and contains tools to guide the design of realistic phantoms.
An active area of research in magnetic resonance imaging (MRI) is the development of reconstruction algorithms. In particular, the inverseproblem approach is getting popular [58], where one relies on an accurate model of the measurement process and possibly on additional information about the object being imaged.
In general, the development of any reconstruction approach requires that it be evaluated and compared to other existing methods. There are several reasons to rely on simulations in a first step
However, for the results to be meaningful, simulations must be accomplished carefully. For instance, the inversecrime situation, where exactly the same discrete model is used for simulation and reconstruction, leads to artificially good results. In the context of MRI, many developers of algorithms base their simulations on rasterized images. One should just be aware that such testing does not account for the full continuousdomain reality, because it neglects the aliasing that is inherent to spatial discretization. More realistic simulations are required to remove this bias and to ensure that the methods will perform adequately in practice.
A method to obtain resolutionindependent simulations is to formulate the simulation analytically in the continuous domain. This approach goes back to Shepp and Logan [59], who introduced an ellipsebased phantom (SL) for Xray tomography. For MRI, several analytical phantoms have been proposed. The first works, based on the SL phantom, are by Smith et al. [60], followed by Van de Walle et al. [61]. More recently, Koay et al. [62] worked out the MR contribution of an ellipsoid for the 3D extension of the SL phantom. Gach et al. [29] adapted these elliptical phantoms specifically for MRI, introducing realistic physical parameters as well as T_{1} and T_{2} relaxation times. The family of analytical phantoms is extended by two recent works by Greengard and Stucchio [63] that use Gaussian functions, and Ngo et al. [64] that introduce 3D polyhedra.
The attractiveness of currently known analytical phantoms is limited for two reasons. First, the vast majority of currently available phantoms (except [64]) use ellipses as basic elements. While such simple shapes have the advantage of mathematical tractability, they do not lend themselves well to the generation of images with realistic anatomical features. Secondly, to the best of our knowledge, no analytical phantom has been proposed that would take into account MRI receivingcoil sensitivities in the context of the simulation of parallel MRI experiments [3].
In this work, we extend the class of available analytical phantoms by introducing regions parameterized by spline contours which are general enough to reproduce polygons and Bézier contours. Our shapes are well suited to the description of realistic anatomical regions [65]. To accurately simulate image formation in parallel MRI, we also make use of analytical models for the coil sensitivity maps. Specifically, we investigate the use of two classes of basis functions—polynomials [3] and complex sinusoids—which both have the ability to generate maps that are physically realistic. These parametric forms are used to derive closedform solutions for the MRI coil data. We have implemented and tested both models.
This chapter is organized as follows: in Section 4.2, we present the different models considered for the parallel MRI measurement process, the analytical phantom, and the coil sensitivities. We motivate and compare the polynomial and the proposed sinusoidal models. In Section 4.3, we propose the main theoretical elements that make the analytical MRI simulation possible, deferring the more technical considerations until Appendices A.1, A.2, and A.3. Finally, we present in Section 4.4 the experiments that validate our implementation of the theoretical tools and an application that quantifies the bias of rasterized simulations on linear and nonlinear reconstructions, in a typical parallel MRI setup.
In this section, we present the MRI measurement model and building blocks that are used to define our phantom.
We use the wellestablished linear model for parallel MRI that relates the object ρ to the kspace signal m_{Sn} observed by each receiving coil, via the Fourier integral
m_{Sn}(k) =  ∫  S_{n}(r)ρ(r)e^{−j 2π k·r}d r, (1) 
where S_{n} accounts for the sensitivity map of the nth receiving channel. We refer to Chapter 2 for more details, in particular, the relation with the BiotSavart law.
We mathematically define the phantom ρ as a simple function, involving R regions R_{i} of constant intensity ρ_{i}
ρ(r) = 
 ρ_{i}χ_{ Ri}(r). (2) 
The term region refers to a connected and bounded set of ℝ^{d}. The symbol χ_{ R} denotes the characteristic function of a region R. Such a phantom has a limited spatial support (∪_{i=1}^{R} R_{i}) that we call a region of interest (ROI).
This model allows us to render realistic phantoms of two kinds
We investigate the first approach in this chapter. The contours that are considered are ellipses, polygons, and quadraticspline curves. We show in Figure 4.1 three such phantoms that we use in our experiments.
For computations, we need to parameterize the complex sensitivity maps. It is commonly admitted that they are smooth and slowlyvarying spatially. It is therefore possible to generate physicallyrealistic sensitivity maps using a reasonably small number of lowpass basis functions. Here, we discuss two models that are wellsuited for this task. They both relate linearly the parameters to the complex sensitivity values. Moreover, their corresponding MRI models involve the Fourier integrals of monomials over the regions of the phantom.
Here, we adopted the multiindex notation r^{α} defined as z^{α}=∏z_{i}^{αi} ∈ ℝ.
This model, first proposed in [3] to represent the local behavior of the sensitivity, assumes that the coil sensitivity S is represented by a polynomial of degree D inside the ROI as
S(r) = 

 s_{d,α} r^{α}, ∀r∈ROI. (4) 
As the degree D increases, the model will reproduce sharper transitions. The number of polynomial coefficients in 2D is N_{p}=(D+1)(D+2)/2.
The corresponding MR response is given by
m_{S}(k) = 
 ρ_{i} 

 s_{d,α}f_{ Ri}^{α}(2π k). (5) 
Alternatively, the coil sensitivity is defined by the linear combination of complex exponentials
S(r) = 
 s_{v} e^{jr·v}, ∀r∈ROI. (6) 
We propose to constrain the problem to the angular frequencies v on a Cartesian grid with spacings that correspond to twice the considered field of view (FOV). The lowfrequency properties are ensured by only considering the L× L angular frequencies around the origin (see Figure 4.2).
Similarly to the effect of the polynomial degree D, an increase in the parameter L allows one to reproduce sharper transitions. The number of coefficients in 2D is given by N_{s} = L^{2}. The corresponding MR response is given by
m_{S}(k) = 
 ρ_{i} 
 s_{v}f_{ Ri}^{0}(2π k−v). (7) 
In order to evaluate and compare the ability of the two models to describe realistic sensitivity maps, we considered a 256× 256 rasterization of the SL phantom and the 27 648 pixels of its ROI. Using BiotSavart’s law (9), we simulated the complex sensitivity maps of a 24channel circular head coil array (FOV: 28 cm, distance to center: 17 cm, radius: 5 cm) distributed around the phantom. In Figure 4.3, the average fitting properties of the two models are presented as a function of the number of parameters.
We observe that the fitting accuracy of both models rapidly increases with the number of parameters, with a sensible advantage for the sinusoidal model. The downside is an increased condition number for the fitting operations. With respect to that criterion, the sinusoidal model behaves also better. The maximal spatial errors are comparable for both models.
In this section, we present the theoretical tools that are necessary to derive the analytical expression of the MRI measurements. Proofs and additional calculation details are provided in Appendices A.1, A.2, and A.3.
The models presented in the previous section allow us to decompose the analytical MRI measurements into Fourier integrals of the sensitivity over the regions that compose the phantom. Depending on the type of region or sensitivity model, we propose tailored methods to decompose the analytical response as a sum of special functions that can be computed accurately and rapidly. In Figure 4.4, we present the road map of these decompositions that are defined and worked out in the sequel.
Let us consider an elliptical region E in 2D parameterized by its center r_{c}, the angle θ formed between its semimajor axis A and the abscissa, and its semiminor axis B. The linear transformation
r ↦ u= D^{−1}R^{T}  ⎛ ⎝  r − r_{c}  ⎞ ⎠  , (8) 
with D =diag(A,B) and R the rotation matrix of angle θ, maps E into a unit disk, that is to say, E = {u  u ≤ 1}. The Fourier transform of the unit disk involves the functions
G_{n}(x) = J_{n}  ⎛ ⎝  ⎪⎪ ⎪⎪  x  ⎪⎪ ⎪⎪  ⎞ ⎠  /  ⎪⎪ ⎪⎪  x  ⎪⎪ ⎪⎪  ^{n}, (9) 
where J_{n} denotes the nth order Bessel function of the first kind [66].
Using the sinusoidal sensitivity model, the integral f_{ E}^{0} can be worked out [61] as
f_{ E}^{0}(ω) = 2πDe^{−jω·rc}G_{1}  ⎛ ⎝  DR^{T}ω  ⎞ ⎠  , (10) 
where D represents the absolute value of the determinant of matrix D.
When considering the polynomial sensitivity model, we suggest to first consider the change of variables (8), rather than computing f_{ E}^{α} directly. We write that
∫ 
 u^{α}e^{−j ω·r}d r = 2πDj^{α}e^{−jω·rc}  ⎛ ⎜ ⎜ ⎝ 
 ⎞ ⎟ ⎟ ⎠  ⎛ ⎝  DR^{T}ω  ⎞ ⎠  . (11) 
The interesting point is that the partial derivatives ∂^{α}G_{1}/∂x^{α} can be decomposed recursively as a sum of G_{n} thanks to the property
∇G_{n}(x) = −xG_{n+1}(x). (12) 
The coefficients of the polynomial in terms of the new coordinates (8) are required to satisfy
S(r) = 

 s_{d,α} r^{α}= 

 t_{d,α} u^{α}. (13) 
They can be computed by inverting the matrix that relates the N_{p} coefficients to the sensitivity values at N≥ N_{p} randomly chosen points in terms of the new coordinates.
The MR contribution of such an elliptical contour is presented in the upper part of Table 4.1.
Model Contribution of the region to MR measurements Ellipse Sinusoidal 2πDe^{−2πjk·rc}∑_{v} s_{v}e^{jv·rc}G_{1}( DR^{T}(2πk−v)) Polynomial 2πDe^{−2πjk·rc}∑_{d=0}^{D} ∑_{α = d} j^{α}t_{d,α} ∂^{α}G_{1}/∂x^{α}( 2πDR^{T}k) Bézier Sinusoidal ∑_{v} s_{v} f_{ B}^{0} (2πk−v) Polynomial ∑_{d=0}^{D} ∑_{α = d} s_{d,α} f_{ B}^{α} (2πk)
In this section, we first provide relations for the computation of the ddimensional Fourier transform of a monomial delimited by a connected subset B of ℝ^{d}. With methods that are similar to the ones employed in [67], we show how to decompose the ddimensional Fourier integral into a sum of integrals over the contour ∂ B. These summed integrals are of reduced dimensionality. In a second step, we show how quadraticspline curves involve a family of 1D integrals.
We show that the surface integral f_{ B}^{α} in (3) can be decomposed into a sum of contour integrals.

f_{ B}^{α} (ω) = j 
 ⎛ ⎜ ⎜ ⎜ ⎝ 
 ⎞ ⎟ ⎟ ⎟ ⎠ 
 α−m!C_{α}^{m}g_{ B}^{m} (ω), (16) 
f_{ B}^{α} (0) = g_{ B}^{α}(0). (17) 
The consequence of Theorem 1 is that the ddimensional integral f_{ B}^{α} can be decomposed into a sum of (d−1)dimensional integrals. The proof is provided in Appendix A.1.
Note that the case ω=0, which corresponds to the calculation of the moments of the region, has been worked out first by Jacob et al. [68] for parametric 2D spline contours.
The region B is defined by its boundary, the contour ∂ B. In 2D , a convenient way to parameterize the contour is by the use of a Bspline generating function ϕ such that
∀ r∈∂ B,∃ t∈ℝ, r(t) = 
 c_{p}ϕ(t−p). (18) 
The considered contour is closed. Consequently, the vectorvalued function r must be periodic. In addition, the number N of coefficients c_{p} that characterize the curve must be finite. The simplest way to satisfy these constraints is to impose that the sequence of coefficients c_{p} be Nperiodic. This enforces the Nperiodicity of r.
If we note ϕ_{p} the Nperiodized version of ϕ, the contour is parameterized either globally as
∀ t∈  ⎡ ⎣  0,N  ⎡ ⎣  , r(t) = 
 c_{q}ϕ_{p}(t−q) (19) 
or piecewise, with 0≤ t=n+λ<M, n∈ 0… N−1 and λ∈[0,1[, as
r(λ +n)= 
 c_{n−q}ϕ_{p}(λ+q). (20) 
We introduce the notation z^{⊥} for the vector perpendicular to z with same norm and pointing outwards the region B at the considered point (see Figure 4.5). We write r′(t) = ∂ r/∂ t(t). The piecewise representation of the contour (20) can be exploited to decompose the contour integral of interest, for instance (14) or (15), which leads to the form
∫ 
 F(r)·nd σ= 
 ∫ 

 ·r′^{⊥}(q+λ)d λ. (21) 
In the sequel, we focus on contours represented by linear and quadratic Bsplines. The former describe polygons while the latter give a piecewise description of quadratic Bézier curves. Three equivalent piecewise representations can be useful and are given in Table 4.2 with their relationships.
g_{ B}^{α} (ω) = 

 e^{−j ω·rn} 
 d_{n,i} h^{(i)}(ω·β_{n},ω·γ_{n}) (23) 
g_{ B}^{α} (0) = 

 d′_{n,i} h^{(i)}(0,0), (24) 

 ω· 
 , (25) 

 e_{k}· 
 . (26) 
The values h^{(m)}(a,b) follow a threeterm recurrence relation [69]. More details on their numerical computation are given in Appendix A.2.
Note that the piecewise parameterization of the contour of a polygon corresponds to the particular case of a quadratic parameterization with β_{n} = r_{n+1}−r_{n} and γ_{n} = 0. Such simpler polygonal models with homogeneous sensitivities have been considered in prior work [63, Prop. 3.2] using a similar formulation.
Our implementation uses Matlab 7.12 (Mathworks, Natick). The experiments run on a 64bit 8core computer, clock rate 2.8 GHz, 8 GiB RAM (DDR2 at 800 MHz), Mac OS X 10.6.7.
We implemented the analytical computations as described by the scheme in Figure 4.4, with double float precision. For efficient computations of the error function of a complex variable, we coded the critical parts of erfz
in C++/MEX, with POSIX multithreading, following Marcel Leutenegger’s recommendations. The code implementing Theorem 1 utilizes Matt Fig’s npermutek
. The rasterization of splinedefined regions, which is performed without approximation, partly relies on Bruno Luong’s MEX implementation of insidepoly
.
Our package also includes graphical tools to design the analytical phantoms. For purposes of adequate visualization, export to the popular vectorgraphics formats SVG 1.1 and PDF (via the PGF/Tikz L^{A}T_{E}X package) is supported. The package is distributed in order to provide sensitivity fitting, phantomdesign interface, analytical simulation tools, and to allow replication of the experiments of this section.
Unlike the sinusoidal model which is very robust to numerical errors, our current implementation of the threeterm recurrence relation (see Appendix A.2) leads to instabilities when using the polynomial model. The theoretical relation h^{(m)}(a,b)≤ 1/(m+1) is sometimes violated for orders m≥ 2 and large values of the first argument. This prevented us to present valid simulations of piecewise quadratic contours combined with a polynomial sensitivity. Given the comparison of the two models in Section 4.2.3, we considered the sinusoidal model with parameter L=7, that is N_{s}=49 in Figure 4.3, which lead to accurate representations of the physical sensitivities and numerically tractable inversions.
As an alternative to our analytical method, we consider the traditional simulation procedure that consists in i) sampling the phantom with a grid of a given size and ii) resampling the DFT of this discrete image according to the desired kspace trajectory. We call this procedure a rasterized simulation. It is expected to be consistent with our analytical method only when considering an infinitely dense sampling.
For reconstructions, we consider an optimization problem of the form
x^{⋆}= arg 
 ⎪⎪ ⎪⎪  m−Ex  ⎪⎪ ⎪⎪  _{2}^{2} + λ P(x), (27) 
where x represents an image, x^{⋆} is the reconstructed one, m is the concatenated scanner data vector, and P is a regularization function. The MRI encoding matrix E is formed as defined in (7) using the usual implicit choice of Dirac’s delta for the generating function.
We used two types of regularizations in our experiments
Please refer to Chapter 3 for more details on these regularization schemes.
As first validation, we consider the simple phantom composed of a rectangular region that is represented in Figure 4.1. Under a proper change of variables, it yields a square and its Fourier transform is given by a product of sinc functions. This phantom is composed of a polygon and consequently falls in the category of the splinedefined contours. We test the accuracy of our proposed simulation method and of the rasterized approach against the closedform solution. To do so, we consider the MR response associated with a homogeneous receiving coil sensitivity and a 256× 256 Cartesian kspace sampling. The simulation errors are reported in Tables 4.3 and 4.4.
Resolution NRMSE max. error max. error in kspace inverse DFT 256 5.58e02 5.5e03 5.5e01 352 2.51e02 2.0e03 1.9e01 400 2.01e02 1.5e03 1.6e01 512 1.25e02 1.1e03 1.0e01 704 7.45e03 5.5e04 6.2e02 800 6.04e03 5.3e04 5.4e02 1024 3.85e03 3.6e04 3.9e02 1408 1.59e03 1.2e04 1.5e02 1600 1.27e03 1.0e04 1.2e02 2048 8.32e04 5.8e05 6.7e03
As expected, the error of rasterized simulations decreases when the sampling density increases. Meanwhile, the accuracy of our analytical implementation is as good as the machine double float
precision would allow. Thus, we conclude that we can indistinctly use the closedform ground truth or our proposed analytical model in the conditions of Section 4.4.2.
We now use our analytical phantom as a gold standard to evaluate the accuracy the measurements obtained from rasterized simulations. We consider the SL and brain phantoms. The single sensitivity map is computed using BiotSavart’s law and is approximated on the support of each phantom with the sinusoidal model. The kspace is on a 128× 128 Cartesian grid. Errors are reported in Tables 4.5 and 4.6.
Resolution NRMSE max. error max. error in kspace inverse DFT 128 1.45e01 1.1e02 2.3e01 176 9.26e02 5.7e03 1.4e01 256 5.45e02 3.3e03 1.1e01 352 3.48e02 2.6e03 7.9e02 512 2.13e02 1.4e03 5.1e02 704 1.02e02 6.6e04 2.0e02 1024 6.70e03 4.1e04 2.1e02 1408 4.06e03 2.4e04 1.4e02 2048 2.03e03 1.5e04 4.8e03 2816 1.49e03 9.5e05 5.6e03
Resolution NRMSE max. error max. error in kspace inverse DFT 128 2.76e01 2.9e02 4.7e01 176 1.79e01 1.6e02 3.0e01 256 9.74e02 8.8e03 1.6e01 352 5.38e02 4.9e03 1.1e01 512 2.85e02 2.6e03 6.5e02 704 2.01e02 1.7e03 3.9e02 1024 1.28e02 1.0e03 3.3e02 1408 6.23e03 6.1e04 1.3e02 2048 3.34e03 3.0e04 7.0e03 2816 2.03e03 1.7e04 5.0e03
We observe that the errors decrease with the same trend as in the rectangle case, which strongly suggests that our gold standard is accurate. Meanwhile, for a given sampling density, the errors occurring with the SL phantom are consistently larger than the ones corresponding to the brain phantom. This is explained by the fact that the SL phantom presents edge transitions of larger intensity.
Let us consider the function f(u)=Sρ(Mu) which depends on the spatial sampling step matrix M. According to (1), the analytical MR data are given by m_{S}(k) = Mf^{^}(2πMk).
When the benefits of an analytical model are forsaken, the MRI data are generated from a rasterized version of the phantom and the sensitivity, using the (nonnecessarily uniform) discrete Fourier transform (DFT)
m_{M}(k) = MF  ⎛ ⎝  e^{−2πjMk}  ⎞ ⎠  , (28) 
with Mk _{∞}≤ 1/2 and
F  ⎛ ⎝  e^{2jπν}  ⎞ ⎠  = 
 f(p)e^{−2jπp·ν} = 

 ⎛ ⎝  ν+q  ⎞ ⎠  . (29) 
The righthand side of (29) can be worked out using Poisson’s summation formula. The terms with q≠ 0 represent the aliasing that occurs with rasterized simulations. Due to the intrinsically discontinuous nature of the phantom ρ, the Fourier transform f^{^} decreases slowly, leading to significant aliasing artifacts. However, as the sampling density increases (Tr(M)→ 0), the impact of aliasing is reduced, as we saw in Section 4.4.2.
Let us define an ideal antialiasing filter h in the Fourier domain as
 (ν) = 
 (30) 
For normalized frequencies ν such that ν _{∞}≤ 1/2, the analytical simulation (unaliased) is characterized as the DFT of the samples of the lowpassfiltered continuous signal
 ⎛ ⎝  ν  ⎞ ⎠  = 
 ⎛ ⎝  h*f  ⎞ ⎠  (p)e^{−2jπp·ν}, (31) 
where (h*f) represents the spatial continuous convolution of h and f.
When using a full Cartesian kspace sampling, the classical approach to reconstruction is to perform an inverse DFT. In this case, the samples of the signal f will be perfectly recovered out of the rasterized simulation (29) which is not desired because it conceals the existence of the Gibbs phenomenon due to the antialiasing filter (see, for instance, [71]). By contrast, the data provided by our analytical model lead to a fairer reconstruction where the Gibbs phenomenon appears. This effect is illustrated in Figure 4.6.
Counterintuitively, the reconstructions out of rasterized simulations lead to aliasing effects that have a positive impact on visual quality. This situation, which occurs when the same model is used for both simulation and reconstruction, is sometimes referred to as “inverse crime”. It arises because of the artificially imposed consistency between the computational forward models used for simulation and reconstruction. In such an inversecrime situation, the continuous nature of the underlying physical model is not taken into account.
We consider a plausible pMRI setup. It involves an array of 8 receiver headcoils that are uniformly distributed around the phantom. The corresponding sensitivity maps are computed thanks to BiotSavart’s law. Spiral and EPI kspace trajectories are considered, both supporting a 256× 256 reconstruction matrix with reduction factor R=4. The channel data are generated using our analytical method as well as 256× 256 and 512× 512 rasterized simulations (see Section 4.4.3). The same realization of complex Gaussian noise is added to the simulated data with different intensities, according to three scenarios: very low noise (40 dB SNR), normal data (30 dB SNR), and very noisy data (20 dB SNR). Reconstructions are performed using quadratic (Tikhonov linear solution) and TV regularizations. The reconstruction algorithms exploit the same forward model, in the form of the same encoding matrix E. The experiments only differ in terms of the input data. The regularization parameter is tuned to optimize the SER with respect to the groundtruth phantom (256× 256 rasterization of the phantom). We report our results in Table 4.7 for the spiral trajectory and in Table 4.8 for the EPI experiments. Reconstructed images are shown in Figures 4.7 and 4.8, together with their error maps, in order to illustrate the impact of the inversecrime situation (the 256× 256 rasterized simulation) in the different scenarios.
Channel data SNR 40dB 30dB 20dB Sampling density 256 512 256 512 256 512 Linear SER 24.61 19.92 20.31 17.99 14.09 13.45 Bias 5.07 0.37 2.56 0.24 0.75 0.11 TV SER 33.75 20.80 27.60 20.26 19.61 17.72 Bias 13.45 0.49 7.75 0.42 2.43 0.54
Channel data SNR 40dB 30dB 20dB Sampling density 256 512 256 512 256 512 Linear SER 36.25 20.77 26.30 19.79 16.73 15.31 Bias 16.02 0.54 6.95 0.44 1.61 0.19 TV SER 42.25 20.98 32.75 20.73 23.92 19.29 Bias 21.85 0.58 12.57 0.55 5.02 0.39
The reconstructions in the spiral experiment are penalized compared to the EPI ones, in the sense that the highfrequency corners of the kspace are not sampled which leads to slightly inferior resolution. This explains that, all other parameters remaining constant, the EPI reconstructions outperform the spiral ones qualitatively and quantitatively.
We observe that the reconstructions from rasterized simulations consistently outperform the ones obtained from analytical measurements. While large differences can occur between the inversecrime scenario (the 256× 256 rasterized simulations) and the analytical simulation data, the 512× 512 simulations yield much closer performance, with at most a 0.6 dB SER difference. This is explained by the reduced aliasing artifacts when doubling the sampling density (see Section 4.4.3). As expected for this type of piecewiseconstant phantom, the TV reconstructions consistently outperform the linear ones. Whatever the simulation method is, TV brings a significant improvement in the very noisy scenario. However, for the other scenarios (SNR 30dB and 40dB), the improvement over linear reconstruction is modest when using the analytic measurements, whereas it is artificially spectacular using the 256× 256 rasterized simulations. We believe that our quality assessment, obtained analytically, offers fairer predictions of the practical worth of a reconstruction method than its overly optimistic rasterized version.
We proposed a method to develop realistic analytical phantoms for parallel MRI. Our analytical phantom approach offers strong advantages for the quantitative validation of MRI and pMRI reconstruction software: it is flexible enough to represent general imaging targets, it provides highly accurate representation of the physical continuous model and avoids overly optimistic reconstructions. This kind of framework is also applicable to the assessment of advanced MRI reconstruction methods such as autocalibrating parallel imaging, B0 correction [2], motion correction [4,5], or higher order field imaging [7].
Implementations of the phantom are made available to the community.