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 piecewise-polynomial 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 well-defined data in both the spatial and k-space domains. Our main contribution is the closed-form 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 inverse-crime 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 inverse-problem 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 inverse-crime 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 continuous-domain 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 resolution-independent simulations is to formulate the simulation analytically in the continuous domain. This approach goes back to Shepp and Logan [59], who introduced an ellipse-based phantom (SL) for X-ray 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 3-D extension of the SL phantom. Gach et al. [29] adapted these elliptical phantoms specifically for MRI, introducing realistic physical parameters as well as T1 and T2 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 3-D 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 receiving-coil 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 closed-form 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 well-established linear model for parallel MRI that relates the object ρ to the k-space signal mSn observed by each receiving coil, via the Fourier integral
| mSn(k) = | ∫ | Sn(r)ρ(r)e−j 2π k·rd r, (1) |
where Sn accounts for the sensitivity map of the n-th receiving channel. We refer to Chapter 2 for more details, in particular, the relation with the Biot-Savart law.
We mathematically define the phantom ρ as a simple function, involving R regions Ri 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=1R Ri) 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 quadratic-spline curves. We show in Figure 4.1 three such phantoms that we use in our experiments.
Figure 4.1: Phantoms parameterized by elliptical and Bézier-defined regions. From left to right: rectangle phantom used for validation, the Shepp-Logan phantom, and a proposed brain phantom. The PDF and SVG versions of the SL and brain phantoms are available online.
For computations, we need to parameterize the complex sensitivity maps. It is commonly admitted that they are smooth and slowly-varying spatially. It is therefore possible to generate physically-realistic sensitivity maps using a reasonably small number of low-pass basis functions. Here, we discuss two models that are well-suited 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 multi-index notation rα defined as zα=∏ziα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) = |
|
| sd,α rα, ∀r∈ROI. (4) |
As the degree D increases, the model will reproduce sharper transitions. The number of polynomial coefficients in 2-D is Np=(D+1)(D+2)/2.
The corresponding MR response is given by
| mS(k) = |
| ρi |
|
| sd,αf Riα(2π k). (5) |
Alternatively, the coil sensitivity is defined by the linear combination of complex exponentials
| S(r) = |
| sv ejr·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 low-frequency 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 2-D is given by Ns = L2. The corresponding MR response is given by
| mS(k) = |
| ρi |
| svf Ri0(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 Biot-Savart’s law (9), we simulated the complex sensitivity maps of a 24-channel 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 2-D parameterized by its center rc, the angle θ formed between its semimajor axis A and the abscissa, and its semiminor axis B. The linear transformation
| r ↦ u= D−1RT | ⎛ ⎝ | r − rc | ⎞ ⎠ | , (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
| Gn(x) = Jn | ⎛ ⎝ | ⎪⎪ ⎪⎪ | x | ⎪⎪ ⎪⎪ | ⎞ ⎠ | / | ⎪⎪ ⎪⎪ | x | ⎪⎪ ⎪⎪ | n, (9) |
where Jn denotes the n-th order Bessel function of the first kind [66].
Using the sinusoidal sensitivity model, the integral f E0 can be worked out [61] as
| f E0(ω) = 2π|D|e−jω·rcG1 | ⎛ ⎝ | DRTω | ⎞ ⎠ | , (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 ω·rd r = 2π|D|j|α|e−jω·rc | ⎛ ⎜ ⎜ ⎝ |
| ⎞ ⎟ ⎟ ⎠ | ⎛ ⎝ | DRTω | ⎞ ⎠ | . (11) |
The interesting point is that the partial derivatives ∂|α|G1/∂xα can be decomposed recursively as a sum of Gn thanks to the property
| ∇Gn(x) = −xGn+1(x). (12) |
The coefficients of the polynomial in terms of the new coordinates (8) are required to satisfy
| S(r) = |
|
| sd,α rα= |
|
| td,α uα. (13) |
They can be computed by inverting the matrix that relates the Np coefficients to the sensitivity values at N≥ Np 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.
Table 4.1: Closed-form MR contribution of elliptical and quadratic Bézier regions for the proposed sensitivity models.
Model Contribution of the region to MR measurements Ellipse Sinusoidal 2π|D|e−2πjk·rc∑v svejv·rcG1( DRT(2πk−v)) Polynomial 2π|D|e−2πjk·rc∑d=0D ∑|α| = d j|α|td,α ∂|α|G1/∂xα( 2πDRTk) Bézier Sinusoidal ∑v sv f B0 (2πk−v) Polynomial ∑d=0D ∑|α| = d sd,α f Bα (2πk)
In this section, we first provide relations for the computation of the d-dimensional 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 d-dimensional 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 quadratic-spline curves involve a family of 1-D 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αmg Bm (ω), (16) |
| f Bα (0) = g Bα(0). (17) |
The consequence of Theorem 1 is that the d-dimensional 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 2-D spline contours.
The region B is defined by its boundary, the contour ∂ B. In 2-D , a convenient way to parameterize the contour is by the use of a B-spline generating function ϕ such that
| ∀ r∈∂ B,∃ t∈ℝ, r(t) = |
| cpϕ(t−p). (18) |
The considered contour is closed. Consequently, the vector-valued function r must be periodic. In addition, the number N of coefficients cp that characterize the curve must be finite. The simplest way to satisfy these constraints is to impose that the sequence of coefficients cp be N-periodic. This enforces the N-periodicity of r.
If we note ϕp the N-periodized version of ϕ, the contour is parameterized either globally as
| ∀ t∈ | ⎡ ⎣ | 0,N | ⎡ ⎣ | , r(t) = |
| cqϕp(t−q) (19) |
or piecewise, with 0≤ t=n+λ<M, n∈ 0… N−1 and λ∈[0,1[, as
| r(λ +n)= |
| cn−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) |
Figure 4.5: Example of a quadratic-spline-defined region with N=6 control points. The boundary is counter-clockwise and the normal vector is outward-pointing. The contour is represented piecewise by quadratic Bézier curves.
In the sequel, we focus on contours represented by linear and quadratic B-splines. 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 |
| dn,i h(i)(ω·βn,ω·γn) (23) |
| g Bα (0) = |
|
| d′n,i h(i)(0,0), (24) |
|
| ω· |
| , (25) |
|
| ek· |
| . (26) |
The values h(m)(a,b) follow a three-term 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 = rn+1−rn 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 64-bit 8-core 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 spline-defined 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 vector-graphics formats SVG 1.1 and PDF (via the PGF/Tikz LATEX package) is supported. The package is distributed in order to provide sensitivity fitting, phantom-design 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 three-term 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 Ns=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 k-space 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 | ⎪⎪ ⎪⎪ | 22 + λ 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 spline-defined contours. We test the accuracy of our proposed simulation method and of the rasterized approach against the closed-form solution. To do so, we consider the MR response associated with a homogeneous receiving coil sensitivity and a 256× 256 Cartesian k-space sampling. The simulation errors are reported in Tables 4.3 and 4.4.
Table 4.3: Errors of our analytical simulations with the Rectangle phantom.
NRMSE max. error max. error in k-space inverse DFT 1.5e-15 2.8e-16 7.0e-15
Table 4.4: Rasterized simulation errors with the Rectangle phantom.
Resolution NRMSE max. error max. error in k-space inverse DFT 256 5.58e-02 5.5e-03 5.5e-01 352 2.51e-02 2.0e-03 1.9e-01 400 2.01e-02 1.5e-03 1.6e-01 512 1.25e-02 1.1e-03 1.0e-01 704 7.45e-03 5.5e-04 6.2e-02 800 6.04e-03 5.3e-04 5.4e-02 1024 3.85e-03 3.6e-04 3.9e-02 1408 1.59e-03 1.2e-04 1.5e-02 1600 1.27e-03 1.0e-04 1.2e-02 2048 8.32e-04 5.8e-05 6.7e-03
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 closed-form 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 Biot-Savart’s law and is approximated on the support of each phantom with the sinusoidal model. The k-space is on a 128× 128 Cartesian grid. Errors are reported in Tables 4.5 and 4.6.
Table 4.5: Rasterized simulation errors versus our analytical simulations with the Brain phantom.
Resolution NRMSE max. error max. error in k-space inverse DFT 128 1.45e-01 1.1e-02 2.3e-01 176 9.26e-02 5.7e-03 1.4e-01 256 5.45e-02 3.3e-03 1.1e-01 352 3.48e-02 2.6e-03 7.9e-02 512 2.13e-02 1.4e-03 5.1e-02 704 1.02e-02 6.6e-04 2.0e-02 1024 6.70e-03 4.1e-04 2.1e-02 1408 4.06e-03 2.4e-04 1.4e-02 2048 2.03e-03 1.5e-04 4.8e-03 2816 1.49e-03 9.5e-05 5.6e-03
Table 4.6: Rasterized simulation errors versus our analytical simulations with the SL phantom.
Resolution NRMSE max. error max. error in k-space inverse DFT 128 2.76e-01 2.9e-02 4.7e-01 176 1.79e-01 1.6e-02 3.0e-01 256 9.74e-02 8.8e-03 1.6e-01 352 5.38e-02 4.9e-03 1.1e-01 512 2.85e-02 2.6e-03 6.5e-02 704 2.01e-02 1.7e-03 3.9e-02 1024 1.28e-02 1.0e-03 3.3e-02 1408 6.23e-03 6.1e-04 1.3e-02 2048 3.34e-03 3.0e-04 7.0e-03 2816 2.03e-03 1.7e-04 5.0e-03
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 mS(k) = |M|f^(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 (non-necessarily uniform) discrete Fourier transform (DFT)
| mM(k) = |M|F | ⎛ ⎝ | e−2πjMk | ⎞ ⎠ | , (28) |
with ||Mk ||∞≤ 1/2 and
| F | ⎛ ⎝ | e2jπν | ⎞ ⎠ | = |
| f(p)e−2jπp·ν = |
|
| ⎛ ⎝ | ν+q | ⎞ ⎠ | . (29) |
The right-hand 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 anti-aliasing 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 lowpass-filtered 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 k-space 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 inverse-crime 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 head-coils that are uniformly distributed around the phantom. The corresponding sensitivity maps are computed thanks to Biot-Savart’s law. Spiral and EPI k-space 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 ground-truth 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 inverse-crime situation (the 256× 256 rasterized simulation) in the different scenarios.
Table 4.7: Reconstruction bias of rasterized simulations on Quadratic and TV regularization reconstructions of the brain phantom for the spiral SENSE experiment. Optimized SER and corresponding bias are shown in dB.
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
Table 4.8: Reconstruction bias of rasterized simulations on Quadratic and TV regularization reconstructions of the brain phantom for the EPI SENSE experiment. Optimized SER and corresponding bias are shown in dB.
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
Figure 4.7: Reconstructed brain phantom and error maps for the spiral SENSE experiments.
Figure 4.8: Reconstructed brain phantom and error maps for the EPI SENSE experiments.
The reconstructions in the spiral experiment are penalized compared to the EPI ones, in the sense that the high-frequency corners of the k-space 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 inverse-crime 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 piecewise-constant 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.
