Method of reconstructing magnetic resonance image data

Abstract

A method of reconstructing magnetic resonance (MR) image data from k-space data. The method includes obtaining k-space data of an image region of a subject; and reconstructing, using a sparse image coding procedure, the MR image data from the k-space data by performing an iterative optimization method. The optimization method includes a data consistency iteration step and a denoising iteration step applied to MR image data generated by the data consistency iteration step. The denoising iteration step incorporates a sparsifying operation to provide a sparse representation of the MR image data for the imaged region as an input to the data consistency iteration step.

Claims

1. A method of reconstructing a magnetic resonance (MR) image from k-space data, the method comprising: providing a computer with k-space data of an image of a subject; and in said computer, performing a reconstructing algorithm, using a sparse image coding procedure in order to reconstruct, the MR image from the k-space data with iterative optimization comprising: a data consistency iteration step comprising a determination of a difference between the reconstructed MR image transformed into k-space and the provided k-space data, wherein the determination of the difference comprises determining an l.sub.2-norm of the result of the difference between the provided k-space data and the reconstructed MR image transformed into k-space; and a denoising iteration step applied to an MR image generated by the data consistency iteration step, wherein the denoising iteration step incorporates a sparsifying operation to provide a sparse representation of the MR image for the imaged region as an input to the data consistency iteration, wherein the denoising iteration step is of the form $\hat{x}=\arg \underset{x}{\min }h\left(x\right),h\left(x\right),$ h(x) is a function defined as {sparse consistency component in image space+λ(sparse coefficients associated with a synthesis dictionary for the sparse representation of the MR image data)}, and λ is a sparsity-control parameter a value of which is adjustable to adjust a strength of sparsity for the sparse representation of the MR image.

2. A method as claimed in claim 1, wherein the data consistency iteration step comprises minimising an Augmented Lagrangian with respect to the reconstructed MR image.

3. A method as claimed in claim 1, wherein the denoising iteration step comprises minimising an Augmented Lagrangian with respect to sparse coefficients associated with a synthesis dictionary for the sparse representation of the MR image.

4. A method as claimed in claim 1, wherein the denoising iteration step comprises a determination of the difference between the reconstructed MR image and the sparse representation of the MR image.

5. A method as claimed in claim 4, wherein the determination of the difference comprises determining a l.sub.2-norm of a result of a subtraction of the sparse representation of the MR image from the reconstructed MR image.

6. A method as claimed in claim 1, wherein the data consistency iteration step of an optimization problem is of the form $\hat{x}=\arg \underset{x}{\min }f\left(x\right),$ wherein ƒ(x) is a function defined as {data consistency component in k-space+μ(sparse consistency component in image space)}, wherein μ is a regularization parameter.

7. A method as claimed in claim 1, wherein the sparsifying operation of the denoising iteration step is applied to sparse coefficients associated with a synthesis dictionary for the sparse representation of the MR image, whereby sparse coefficients having a value below a threshold value are set to zero.

8. A method as claimed in claim 7, wherein the threshold value is equal to the value of the sparsity-control parameter.

9. A non-transitory computer readable medium comprising instructions recorded thereon which, when executed by a computer, cause the computer to perform the method as claimed in claim 1.

10. A magnetic resonance (MR) apparatus, comprising: a gradient system configured to apply a magnetic field gradient; an excitation system configured to apply an excitation pulse to a subject and to receive signals from the subject; and a computing system configured: to receive the signals from the excitation system, to execute program code to control the gradient system and the excitation system to acquire k-space data of an image of a subject, and a plurality of navigator signals for the image of the subject, to reconstruct, using a sparse image coding procedure, an MR image from the acquired k-space data by performing an iterative optimization method comprising: a data consistency iteration step comprising a determination of a difference between the reconstructed MR image transformed into k-space and the provided k-space data, wherein the determination of the difference comprises determining an l.sub.2-norm of the result of the difference between the provided k-space data and the reconstructed MR image transformed into k-space; and a denoising iteration step applied to MR image generated by the data consistency iteration step, wherein the denoising iteration step incorporates a sparsifying operation to provide a sparse representation of the MR image for the image of the subject as an input to the data consistency iteration step.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) Examples of the present disclosure will now be described with reference to the accompanying drawings, in which:

(2) FIG. 1 shows a process diagram for an example method according to the first aspect of the present disclosure; and

(3) FIG. 2 shows an example MR apparatus according to the second aspect of the present disclosure;

(4) FIG. 3 shows a schematic representation of a synthesis dictionary in a sparse representation of image data;

(5) FIG. 4 shows an implementation of an embodiment of the disclosure upon 3-dimensional (3D) MRI data; and

(6) FIG. 5 shows examples for the application of the embodiment (“3D-PROST”) of FIG. 4 on a high-resolution phantom. Fully sampled acquisition (FS-left column) is compared to 5-fold and 9-fold accelerated acquisitions (middle and right columns respectively), for both zero-filled (ZF-top row) and 3D-PROST (middle row) reconstructions. The variable density under-sampled Cartesian images with ZF reconstruction show significant blurring and contrast loss, while the 3D-PROST images exhibit sharp edges with faithful preservation of small details (as shown on the cross-section profiles). Reconstructed resolution for ZF and the proposed 3D-PROST technique are shown on the top-left corner (MTF-profiles taken in the region-of-interest box). Differences with statistical significance are identified by *P<0.05 (versus FS).

DETAILED DESCRIPTION

(7) The present disclosure relates to compressed sensing procedures which comprise solving optimization problems comprising a data consistency component and a transform sparsity component.

(8) By way of example, a compressed sensing procedure suitable for use with the present disclosure will now be described which comprises the following steps (1) to (3).

(9) In step (1) of this example sparse image coding procedure, the original acquired k-space data is transformed into MR image data (e.g. a MR image), e.g. through use of a standard iterative SENSE reconstruction, which consist of solving equation (3) without priors (b=Dα=0). SENSE (SENSitivity Encoding) is a well-known imaging technique, readily available to the skilled person, which incorporates knowledge of the sensitivities of the coil elements. The MR image data may suffer from sub-sampling artefacts. However, if an incoherent sampling scheme is used to acquire the k-space, these artefacts may appear to be ‘smeared’ over the image. That is, almost homogeneous noise-like artefacts may appear in the MR image data. The reconstructed MR image data may serve as a starting point for the optimization problem. The optimization problem seeks to find a better solution (a better set of MR image data) that is also consistent with the original acquired k-space data.

(10) In step (2) a data consistency iteration step is performed according to an Augmented Lagrangian using as input the image data obtained from a denoising iteration step (4), below as image data input.

(11) In step (3) a denoising iteration step is performed according to an Augmented Lagrangian using as input the image data obtained from step (2), above. A sparse representation of the image data is applied in step (3) by removing low-value coefficients of the coefficient vector associated with a synthesis dictionary representation of the image data. The purpose of the sparse representation is to attempt to locally separate the desired signal from the noise artefacts. A higher sparsity means that the useful image content is concentrated in a few pixels, while most of the pixels have only a very low signal.

(12) In step (3), the coefficient vector is thresholded to remove noise.

(13) One example thresholding procedure which may be used is referred to as ‘hard thresholding’ in which all coefficients of the coefficient vector with a value less than the threshold are set to zero. As a result of the thresholding procedure, many coefficients of the coefficient vector have a value of 0, and as such there are fewer non-zero coefficients of the coefficient vector. Consequently, the sparsity of the image represented in terms of the coefficients of the coefficient vector is increased. This is schematically illustrated in FIG. 3 in which most of the coefficients of the coefficient vector (α) are zero-valued, even though this represents the image data vector (x) which has no zero valued coefficients/elements. In FIG. 3, zero-values are represented by non-shaded (empty) elements in the schematic vectors (x, and α) and schematic synthesis dictionary (matrix D), and non-zero values are represented by shaded (occupied) elements.

(14) The following examples may use the sparse image coding procedure of the form as described generally above. It will further be appreciated that the present disclosure is not directed towards a particular sparse image coding procedure, but rather the beneficial incorporation sparsity in images of the subject obtained from the plurality of navigator signals into the data consistency and denoising/sparsity iteration steps of the optimization problem of the sparse image coding procedure.

(15) Referring to FIG. 1, there is shown a method of reconstructing MR image data from k-space data.

(16) Step 101 comprises obtaining k-space data of an image region of a subject.

(17) In the example of FIG. 1, the k-space data is three-dimensional (3D) whole-heart coronary MR angiography data (CMRA) acquired using a Cartesian trajectory. The Cartesian trajectory samples the k.sub.y−k.sub.z plane using a Cartesian trajectory with spiral profile ordering. The MR image sequence for generating the 3D CMRA data in this example implementation is a 3D balanced steady-state free precession (bSSFP) sequence, which may be ECG-triggered.

(18) The present disclosure is not limited to the particular k-space acquisition procedure of the example of FIG. 1. In particular, the present disclosure is not limited to 3D whole heart CMRA. Imaging other regions of the subject, such as the liver is within the scope of the present disclosure. Further, obtaining 2D image regions rather than 3D image regions is within the scope of the present disclosure. Further, the present disclosure is not limited to Cartesian sampling, other forms of sampling such as radial sampling and radial sampling with golden angle separation are within the scope of the present disclosure. Further, the present disclosure is not limited to bSSFP sequences. Other MR image sequences such as gradient echo sequences are within the scope of the present disclosure.

(19) Steps 102 and 103 comprise using a sparse image coding procedure to reconstruct the MR image data from the k-space data, the sparse image coding procedure comprises using a sparse image coding procedure to reconstruct the MR image data from the k-space data by iteratively solving an optimization problem comprising a data consistency iteration step 102, and a denoising iteration step 103 applied to MR image data generated by the data consistency iteration step

(20) The first time step 102 is implemented there is no prior information (from step 103), however, each subsequent application of step 102 incorporates the output of step 103 as input 104 in an iterative reconstruction The denoising iteration step 103 incorporates a sparsifying operation to provide a sparse representation of the MR image data for the imaged region as an input (104) to the data consistency iteration step of the optimization problem.

(21) The present disclosure is not limited to the correcting procedure of FIG. 1. Other ways for correcting for motion are within the scope of the present disclosure. Further, the present disclosure does not require that the motion is corrected for in all embodiments.

(22) Two-dimensional image based navigators (iNAVs) were acquired at each cardiac cycle. These low spatial resolution iNAV images are acquired in every heartbeat before the MRI data acquisition. This technique provides a set of images with high-temporal resolution which are then used to estimate respiratory motion of the heart. This approach has the advantage to spatially isolate the moving heart from the surrounding static tissues and also avoid the use of gating window, therefore increasing the total scan efficiency to or close to 100%. Higher scan efficiency also means predictable and shorter scan times. The k-space data is motion-corrected by applying phase shift correction, with the shift values being estimated in the image domain from the iNAVs directly using rigid motion correction.

(23) Referring to FIG. 2, there is shown an example MR apparatus 200 according the second aspect of the disclosure. The MR apparatus 200 comprises a computing system 201, a gradient system 203, and an excitation system 205.

(24) The gradient system 203 is configured to apply a magnetic field gradient. The gradient system 203 may be configured to apply magnetic field gradients along three spatial axes.

(25) The excitation system 205 may comprise a transmitter (not shown) and a receiver (not shown). The excitation system 205 can be an RF system with one or more RF coils (not shown). The excitation system 205 is configured to apply an excitation pulse to the subject and to receive signals from the subject.

(26) The MR apparatus 200 includes a magnet (not shown) for establishing a stationary magnetic field. The magnet can include a permanent magnet, a superconducting magnet or other type of magnet.

(27) The computing system 201 is in communication with the excitation system 205, and the gradient system 203 for controlling these components. The computing system 201 is configured to receive the signals from the excitation system 205.

(28) The computing system 201 is further configured to execute program code to control the gradient system 203 and the excitation system 205 to acquire k-space data of an image region of a subject, and a plurality of navigator signals for the image region of the subject.

(29) The computing system 201 is further configured to use a sparse coding procedure to reconstruct the MR image data from the k-space data. The sparse coding procedure comprises solving an optimization problem comprising a data consistency iteration component and a denoising/sparsity component.

(30) The computing system 201 is configured to execute program code for controlling the MR apparatus 200 to perform the method of the first aspect. The computing system 201 could be an integrated component of the MR apparatus 200. The computing system 201 could be a unit of a desktop computer, a workstation, a server, or a laptop computer.

(31) According to aspects of the disclosure, there is also provided a computer-readable medium having instructions recorded thereon which, when executed by a processing unit, cause the processing unit to perform the method of the first or second aspect.

(32) In other words, computing system is in communication with the excitation system, and the gradient system for controlling these components, and is configured to receive the signals from the excitation system. The computing system is further configured to execute program code to control the gradient system and the excitation system to acquire k-space data of an image region of a subject. The computing system is configured to use a sparse image coding procedure to reconstruct the MR image data from the acquired k-space data by iteratively solving an optimization problem comprising a data consistency iteration step and a denoising iteration step applied to MR image data generated by the data consistency iteration step. The denoising iteration step incorporates a sparsifying operation to provide a sparse representation of the MR image data for said imaged region as an input to the data consistency iteration step of the optimization problem. The computing system may be further configured to obtain a plurality of navigator signals for the image region of the subject.

(33) Data Acquisition:

(34) A prototype under-sampled Cartesian variable density 3D spiral-like Cartesian sampling (VD-CASPR) was implemented on a 1.5T scanner (Siemens, Aera). A 2D image-navigator precedes each VD-CASPR acquisition to enable beat-to-beat 2D translational respiratory motion estimation/correction without data rejection. Five healthy volunteers (3 females, 30±4 years) underwent free-breathing CMRA. Data were acquired in middiastole with an ECG-gated 3D bSSFP sequence, 18-channel body and 32-channel spine coils, fat-saturation and T2-preparation pulses (FOV=320×320×120 mm3, TR/TE=3.35/1.47 ms, FA=90°, bandwidth=890 Hz/pixel, T2-preparation duration=40 ms) with 1.2 mm.sup.3 isotropic resolution and under-sampling factors of 5 and 9.

(35) Reconstruction:

(36) The disclosure, in a preferred embodiments, was implemented and is herein denoted the ‘3D-Patch-based low-rank ReconstructiOn with Self-similariTy learning (3D-PROST)’ method. This integrates anatomical structure information from 3D patch neighbourhoods through sparse representation, exploiting the redundancy of 3D patches in the acquired data itself. The optimization problem iterates between a data consistency step, which reconstructs a high-resolution isotropic MR volume x, and a low-complexity 3D patch-based denoising step, which provides a reconstructed volume as prior for the next step (FIG. 4).

(37) The two sub-problems are solved iteratively in an effective Augmented Lagrangian (AL):

(38) $ℒ\left(x,\alpha \right)=\arg \underset{x,\alpha }{\min }\left\{\frac{1}{2}{.Math.\mathrm{Ex}-k.Math.}_2^2+\lambda {.Math.\alpha .Math.}_0+\frac{\mu }{2}{.Math.x-D\alpha -b.Math.}_2^2\right\}$

(39) Here, λ is a sparsity-control parameter the value of which is adjustable to adjust the strength of sparsity for the sparse representation of the MR image data, μ is a regularization parameter, α is a (sparse) coefficient vector associated with a dictionary matrix (D), and b is the Augmented Lagrangian multiplier. Here, E is the encoding operator (including coils, Fourier operator and sampling), x denotes the undersampled data, D is also the patch-grouping operator (synthesis dictionary). The ∥ ∥.sub.0 norm counts the number of non-zeros element in α and λ>0 controls the strength of sparsity.

(40) The two problems are solved iteratively:

(41) Iteration Stage 1:

(42) MR reconstruction (FIG. 4: Stage 1) is performed using a conjugate gradient (CG) descent and uses the 3D denoised volume (x=Dα) obtained from Stage 2 as prior knowledge (x=b=0 initially, then updated by Stage 2). The data consistency iteration step is of the form of an Augmented Lagrangian:

(43) ${ℒ}_1\left(x\right)=\arg \underset{x}{\min }\left\{\frac{1}{2}{.Math.\mathrm{Ex}-k.Math.}_2^2+\frac{\mu }{2}{.Math.x-D\alpha -b.Math.}_2^2\right\}$
Iteration Stage 2:

(44) A 3D block-matching 2,5 algorithm is used to exploit redundancies in the volume x. The most similar 3D blocks to a reference block are extracted, vectorised and concatenated into a 2D matrix. Sparsity of this low-rank matrix is enforced using complex Singular Value Decomposition (SVD) and by hard thresholding the singular values below a specific threshold A. The de-noised 3D blocks are then placed back to their original positions by averaging (FIG. 4: Stage 2). The denoising iteration step takes the of the form of an Augmented Lagrangian:

(45) 0${ℒ}_2\left(\alpha \right)=\arg \underset{\alpha }{\min }\left\{\frac{1}{2}{.Math.x-D\alpha -b.Math.}_2^2++\lambda {.Math.\alpha .Math.}_0\right\}$
Implementation and Analysis:

(46) The following parameters were empirically selected to provide the best reconstruction quality: patch_size=5×5×5 voxels, window_search=14×14×14 voxels, AL_iterations (stages 1 and 2 above)=4, CG_iterations=30, λ=4, μ=1. The proposed approach was compared to conventional zero-filled (ZF), iterative-SENSE, and Wavelet-based Compressed Sensing (CS) reconstructions.

(47) A prototype free-breathing 3D whole-heart, electrocardiogram (ECG)-triggered, balanced steady-state free precession (bSSFP) sequence with variable density Cartesian under-sampling was implemented. The k.sub.y−k.sub.z phase-encoding plane was sampled following approximate spiral interleaves on the Cartesian grid with variable density along each spiral arm. One spiral arm is acquired per cardiac cycle and is then rotated from one cardiac cycle to the next one. The k.sub.y−k.sub.z plane is then segmented in two sets of concentric rings, the first defining the fully-sampled k-space centre and the second representing the accelerated spiral branches.

(48) The phase-encoding lines within each ring were sorted according to a defined increment angle from 0 to 360°. The sampling of the accelerated branches is accelerated exponentially from the k-space centre to its periphery. The size of the fully-sampled k-space centre was optimized on several datasets and was set to 20% the size of k.sub.y and k.sub.z encoding directions. This under-sampled trajectory ensures a pseudo-random pattern through the cardiac cycle, resulting in incoherent aliasing which spreads irregularly in a noise-like fashion.

(49) The acquisition of each spiral arm was preceded by a 2D iNAV to enable beat-to-beat 2D translational (SI: superior-inferior, RL: right-left) respiratory motion estimation/compensation and 100% scan efficiency. iNAVs were obtained by spatially encoding 14 startup echoes of the bSSFP CMRA sequence. 2D translational motion may be estimated using a template-matching algorithm, with the template selected around the heart during acquisition planning. Motion compensation may be performed by modulating the k-space data with a linear phase shift to a reference position at end-expiration. Motion estimation/compensation may be performed before 3D-PROST reconstruction and may be implemented inline in the scanner software. The acquisition framework is depicted in the sequence diagram of FIG. 4a.

(50) In FIG. 4—‘Data Acquisition’:

(51) Motion correction is performed with an under-sampled 3D variable density spiral-like Cartesian trajectory, (VD-CASPR) preceded by 2D image navigators (iNAV), T2 preparation and fat saturation pulses. iNAVs were used to estimate and correct for the beat-to-beat 2D translational respiratory-induced motion of the heart (Moco).

(52) In FIG. 4—‘Reconstruction’:

(53) The 3D-PROST reconstruction involves two stages of an Augmented Lagrangian optimization scheme as described above. In stage 1, image reconstruction is performed with patch prior and data consistency enforcement. In stage 2, image denoising is performed using 3D block-matching, which groups similar 3D patches in the image, followed by a low-rank approximation of each group using 2D-SVD shrinkage. The de-noised volume from stage 2 is used in the reconstruction process in stage 1 as prior knowledge to guide the reconstruction problem and further reduce the noise.

(54) The general formulation of sparse representation in terms of a redundant dictionary considers a 3D image I as the approximation I≈Dα, where D is a fixed dictionary and α a sparse vector satisfying the sparsity-inducing condition:
∥α∥.sub.0≤T,
where the l.sub.0-norm counts the number of nonzero elements in α and T is a predefined threshold. Strictly speaking, the image I can be represented with a minimum number of sparse coefficients α in the redundant dictionary D. In the present embodiment, a novel reconstruction algorithm is used which iteratively exploits the structure redundancy in the acquired data I to construct a specific dictionary D for each group of similar 3D patches. The 3D-PROST scheme for isotropic CMRA reconstruction using 3D patch redundancy is formulated as the following unconstrained optimization on the sparse coefficients α:

(55) $\arg \underset{\alpha }{\min }\left\{\frac{1}{2}{.Math.{\mathrm{AFS}}_{c}D\alpha -k.Math.}_2^2+\lambda {.Math.\alpha .Math.}_0\right\}$

(56) Where S.sub.c are the known coil sensitivities for channel c, F is the Fourier transform, A is the sampling operator, k is the acquired multi-channel k-space data, and λ is the regularization parameter.

(57) The data acquired by each coil (the c.sup.th coil) can be expressed in the form:
m.sub.c=AFS.sub.cI

(58) Here, I=Dα, is the image series to be reconstructed and m.sub.c is the acquired data for the c.sup.th coil element of MRI apparatus, and the coil sensitivities S.sub.c for the coil/channel c may be in the form of a sensitivity map in the space (e.g. Cartesian: x-y) of the image. To implement multi-coil imaging (parallel image reconstruction) the image data from the individual coil elements/channels are concatenated such that:

(59) $m=\left[\begin{array}{c}{m}_1\\ \mathrm{.Math.}\\ m_j\end{array}\right]\mathrm{and}E=A.Math.F.Math.\left[\begin{array}{c}{S}_1\\ \mathrm{.Math.}\\ S_j\end{array}\right],$
where j is the number of coils, then one may employ the multi-coil imaging (parallel image reconstruction) of data acquired by all j coils:
m=EI

(60) One may transform the above minimisation equation into its equivalent constrained minimization problem:

(61) $\arg \underset{I,\alpha }{\min }\left\{\frac{1}{2}{.Math.\mathrm{EI}-k.Math.}_2^2+\lambda {.Math.\alpha .Math.}_0\right\}\mathrm{subject}\mathrm{to}I=D\alpha$

(62) Where I denotes the image to recover. This optimization problem is difficult to solve directly due to the constrained form of the functional. The present embodiment solves this minimization problem using its Augmented Lagrangian formulation, which implies minimizing the augmented Lagrangian L, defined below, with regards to I and α:

(63) ${ℒ}_{3\mathrm{DPROST}}\left(I,\alpha \right)=\arg \underset{x\alpha }{\min }\left\{\frac{1}{2}{.Math.\mathrm{EI}-k.Math.}_2^2+\lambda {.Math.\alpha .Math.}_0+\frac{\mu }{2}{.Math.I-D\alpha -b.Math.}_2^2\right\}$

(64) Here, b represents the Lagrange multiplier associated with the constraint “I=Dα”, and μ≥0 is the penalty parameter. This equation is solved using a variable splitting approach by alternating the minimization with respect to the image I (Iteration Stage 1) and the sparse coefficients α (Iteration Stage 2), followed by an update of the augmented multiplier b, and repeating these three steps until a convergence criterion is satisfied. The rationale behind this splitting approach is that each sub-problem of this alternating minimization is much simpler to solve than the original unconstrained problem.

(65) ITERATION STAGE 1: MR Reconstruction Update

(66) The first sub-problem with regards to the variable I is a conventional MR reconstruction that incorporate the denoised volume Dα (obtained at the end of stage 2), as prior information:

(67) ${ℒ}_1\left(I\right)=\arg \underset{x}{\min }\left\{\frac{1}{2}{.Math.\mathrm{EI}-k.Math.}_2^2+\frac{\mu }{2}{.Math.I-D\alpha -b.Math.}_2^2\right\}$

(68) Differentiating with respect to I, we find that the residual gradient step is:
m=E.sup.HEI−E.sup.HK+μ(I−ω−b)

(69) Where the operator E.sup.H denotes the Hermitian transpose of E and ω=Dα represents the truncated SVD reconstruction obtained at iteration stage 2. For the initialization, ω and b were set to 0, which reduce the reconstruction problem to a standard iterative SENSE with Tikhonov regularization. The gradient descent optimization method may be used to iteratively update the reconstructed volume I:
(t+1)←I.sup.(t)−βm.sup.(t)

(70) Here, the relaxation parameter β can be updated iteratively or be set to a specific value (e.g. β=0.1) to ensure convergence.

(71) ITERATION STAGE 2: 3D Patch-Based Denoising Update

(72) The second sub-problem minimizes with respect to the sparse coefficients α and is given by:

(73) ${ℒ}_2\left(\alpha \right)=\arg \underset{\alpha }{\min }\left\{\frac{1}{2}{.Math.I-D\alpha -b.Math.}_2^2++\lambda {.Math.\alpha .Math.}_0\right\}$

(74) Considering that sparse image coding is a local model representation and that neighbouring patches in CMRA images are highly redundant, this optimization can be performed on an image patch basis. A 3D patch I.sub.k in the previously reconstructed (at Iteration Stage 1) volume I of size N voxels is defined as a small 3D block of size n.sup.3 voxels around the voxel at index k.

(75) Define the operator R.sub.k that extracts the patch I.sub.k from the image I: I.sub.k=(I). Inversely, one can recover the image I from its set of patches in a straightforward manner:
I=(ΣR(I.sub.k))./(ΣR.sub.k.sup.HR.sub.k)=(ΣR.sub.k.sup.H(D.sub.kα.sub.k))./(ΣR.sub.k.sup.HR.sub.k), for k=1 to N.

(76) Where the operator R.sub.k.sup.HR.sub.k is a matrix of same size as I.sub.k with all elements being 1 (i.e. averaging matrix), (α,) represents the concatenation of all (α.sub.k,k) and the operator “./” denotes the element-wise division. The above equation can thus be rewritten as the following patch-based minimization given by:

(77) $\arg \underset{\alpha }{\min }\left\{\frac{1}{2}{.Math.I_k-D_k{\alpha }_k-b_k.Math.}_2^2++\lambda {.Math.{\alpha }_k.Math.}_0\right\}\mathrm{for}k=1,\mathrm{.Math.},N$

(78) Efficient optimization is achieved through the choice of the dictionary D.sub.k which induces the highest sparsity in its associated group of similar patches, and how accurately the sparse coefficients α.sub.k can be recovered from this dictionary. In 3D the present patch-based representation, the self-similarity is considered as a 4D set of similar 3D patches [I.sub.1, . . . , I.sub.k] selected in a large 3D window d, and its associated sparse coefficients [α.sub.1, . . . , α.sub.L] obtained from a dictionary D.sub.k. In order to account for 3D patches, the present embodiment reduces the complexity of the problem by concatenating each similar vectorized 3D patch into a 2D matrix. This 2D matrix, containing a high degree of similarity, exhibits a low-rank structure which can be sparsely approximated using singular value decomposition (SVD). Using the unitary property of the SVD, the minimization problem of the above equation is equivalent to minimizing with regards to the sparse coefficients the following equation,

(79) ${ℒ}_2\left(\alpha \right)=\arg \underset{{\alpha }_k}{\min }\left\{\frac{1}{2}{.Math.{\tilde{\alpha }}_k-{\alpha }_k.Math.}_2^2++\lambda {.Math.{\alpha }_k.Math.}_0\right\}$

(80) Where {tilde over (α)}.sub.k are the sparse coefficients associated with I.sub.k−b.sub.k=Ĩ.sub.k and λ>0 controls the strength of sparsity. The lower the parameter λ, the more accurate the reconstructed solution, at the price of reducing the sparseness. The solution of the above equation has a closed form, and the optimum is obtained using hard-thresholding:
α.sub.k*=H.sub.√{square root over (2λ)}({tilde over (α)}.sub.k)

(81) Here, H.sub.0(.) is the element-wise hard-thresholding operator, defined for a scalar v as: H.sub.0(v)=0.1.sub.|v|>0. In other words, any singular values {tilde over (α)}.sub.k below √{square root over (2λ)} is set to 0. Note that only the singular values are modified, but the singular vectors are unperturbed. The value α.sub.k* is indirectly an estimate of the noise Ĩ.sub.k=I.sub.k−b.sub.k, and thus the hard-thresholding step effectively shrinks the noise.

(82) This step is repeated for each voxel in the volume I and the final de-noised 3D volume ω=Dα is obtained by aggregating the multiple estimates D.sub.kα.sub.k* at each voxel location k=1, . . . . The Lagrange multiplier is then updated at iteration t+1 as:
(t+1)←(t)+r(ω.sup.(t+1)−I.sup.(t+1))

(83) Where the penalty parameter r is set to an appropriate value (e.g. set to 0.1). Optimizations for Stage 1 and Stage 2 are processed iteratively to improve the accuracy of the reconstructed volume. A flowchart illustrating all steps is shown in FIG. 4b.

(84) Image Reconstruction Implementation

(85) Important parameters of interest are the size of patch n, size of neighbourhood window d, number of selected patches L, regularization parameters λ and μ, as well as the number of outer iterations. The size of patch n controls the degree of structural information within each patch. On one hand, a large value of n would capture the most geometric information and leads to a higher level of denoising, while a small value would act as a local filter and would potentially reduce the denoising performance of the algorithm.

(86) An example is to set the size of patches to be 5×5×5 voxels and the selection window d to 14, which gives a good trade-off between reconstruction quality and computation time. The number of selected similar patches L may be set to 40 to avoid high computation cost and excessive memory requirements associated with large values of L. The performance of 3D-PROST may be used, for example, with the values λ=0.1 and μ=0.3.

(87) 3D-PROST reconstruction was compared to iterative SENSE (itSENSE) and a CS reconstruction with wavelet regularization, as implemented in the BART toolbox. The regularization parameter was carefully tuned and set to λCS=0.01 in all studies and the algorithm was stopped after 30 iterations (respectively five for itSENSE), since preliminary testing revealed that these numbers of iterations led to the best reconstructions.

(88) Phantom Study

(89) Acquisition:

(90) A phantom acquisition was performed to evaluate the impact of highly under-sampled acquisitions on the reconstructed resolution. The phantom was composed of parallel bar patterns with in-plane resolution of 1 mm. The relevant imaging parameters were as follows: 3D bSSFP sequence, field-of-view (FOV)=306×306×80 mm3, TE=1.63 ms, TR=3.7 ms, flip angle (FA)=90°, 0.9 mm.sup.3 isotropic resolution, sampling bandwidth=875 Hz/pixel, T2-preparation duration=40 ms.

(91) Acquisitions were performed with under-sampling factors of 5 and 9. Additionally, a fully-sampled acquisition was performed for comparison purposes.

(92) Acquired motion-corrected data was reconstructed with itSENSE, CS and the 3D-PROST of the present embodiment. Under-sampled reconstructions were compared against the fully-sampled reference images. Reconstructed images were reformatted along the right (RCA) and left anterior descendent (LAD) coronary arteries and visible vessel length and vessel sharpness (first 4 cm and full length) were measured. A vessel sharpness of 100% marked an abrupt change in signal intensity whereas a sharpness of 0% indicated the absence of an edge.

(93) Data Analysis:

(94) The modulation transfer function (MTF) was used to quantify the spatial resolution for the under-sampled and fully sampled reconstructions. The MTF was directly obtained using a slanted edge method. A profile (edge spread function) was taken along a sharp edge of the phantom, then differentiated to obtain the associated line spread function, and finally fast Fourier transformed to get the final MTF, which was then normalized. The reconstructed resolution was given by the value of the MTF at 10%. The MTF was obtained from 20 intensity profiles taken on the reference fully-sampled image, the zero-filled (ZF) reconstruction and the 3D-PROST reconstructions. The average MTFs were compared for both accelerations. A paired two-tailed Student t-test was used for statistical analysis, with P<0.05 considered statistically significant.

(95) Results:

(96) Phantom Study:

(97) The total scan times (min:sec) were 11:30 (fully-sampled), 2:13 (5-fold acceleration) and 1:15 (9-fold acceleration) with isotropic resolution of 0.9 mm.sup.3.

(98) The reconstructed images using ZF and 3D-PROST are shown in FIG. 5 in comparison with the reference fully-sampled image. The measured MTFs are shown on the top of each image. Significant aliasing artifacts can be observed on the ZF images, and were particularly pronounced for high undersampling where the depiction of the bar pattern becomes indistinguishable. Conversely, high image quality can be appreciated using the 3D-PROST reconstruction of the present embodiment, with clear and sharp depiction of the phantom structures with high spatial resolution, even for high acceleration (×9), providing excellent visual agreement with the reference fully-sampled image. This observation is confirmed and illustrated in FIG. 5 (bottom) where the intensity profiles taken perpendicularly across the pattern are plotted for each reconstruction. The fine and sharp structures of the bars are well preserved with 3D-PROST reconstruction while ZF reconstruction exhibits blurring and loss of contrast. The spatial resolution, as quantified by the MTF, was preserved with the proposed approach, showing no statistical differences with the reference fully-sampled image (5-fold: 0.90±0.01, P=0.31 and 9-fold: 0.90±0.02, P=0.17).

(99) 2D-PROST:

(100) The disclosure may be implemented to reconstruct 2-dimensional (2D) MRI data, herein referred to as the 2D-PROST reconstruction embodiment. This may be performed with the same parameters as 3D-PROST except for the size of the patches (for example: 2D: 5×5 voxels vs. 3D: 5×5×5 voxels) and the selection window d (for example: 2D: 5×5 voxels vs. 3D: 5×5×5 voxels).

DISCUSSION

(101) The proposed 3D-PROST reconstruction embodiment integrates self-similarity information, by grouping 3D patches with similar structures. Low-rank properties and sparsity of the group are enforced to reduce the noise of the reconstructed volume while the MR reconstruction step was used to recover an isotropic 3D volume and enforce data fidelity. An Augmented Lagrangian formulation efficiently decomposes the main cost function into two sub-problems that have straightforward solutions. The improved performance can be explained by the fact that CMRA images contain a rich amount of correlated 3D structures and therefore high sparsity degree can be achieved by merging this information. The increased sparsity promotes superior denoising and structures recovery in comparison to other established approaches. Similar image quality was achieved for an isotropic 1.2 mm.sup.3 resolution.

(102) Under-sampled reconstruction according to the disclosure, may be used to obtain isotropic sub-millimetre 3D coronary images under free-breathing in ˜5 min predictable scan time, which is crucial for its integration in routine clinical examinations.

(103) At least some of the example embodiments described herein may be constructed, partially or wholly, using dedicated special-purpose hardware. Terms such as ‘component’, ‘module’ or ‘unit’ used herein may include, but are not limited to, a hardware device, such as circuitry in the form of discrete or integrated components, a Field Programmable Gate Array (FPGA) or Application Specific Integrated Circuit (ASIC), which performs certain tasks or provides the associated functionality. In some embodiments, the described elements may be configured to reside on a tangible, persistent, addressable storage medium and may be configured to execute on one or more processors. These functional elements may in some embodiments include, by way of example, components, such as software components, object-oriented software components, class components and task components, processes, functions, attributes, procedures, subroutines, segments of program code, drivers, firmware, microcode, circuitry, data, databases, data structures, tables, arrays, and variables. Although the example embodiments have been described with reference to the components, modules and units discussed herein, such functional elements may be combined into fewer elements or separated into additional elements.

(104) The described and illustrated embodiments are to be considered as illustrative and not restrictive in character, it being understood that only the preferred embodiments have been shown and described and that all changes and modifications that come within the scope of the embodiments as defined in the claims are desired to be protected. It should be understood that while the use of words such as “preferable”, “preferably”, “preferred” or “more preferred” in the description suggest that a feature so described may be desirable, it may nevertheless not be necessary and embodiments lacking such a feature may be contemplated as within the scope of the disclosure. In relation to the claims, it is intended that when words such as “a,” “an,” “at least one,” or “at least one portion” are used to preface a feature there is no intention to limit the claim to only one such feature unless specifically stated to the contrary in the claim. When the language “at least a portion” and/or “a portion” is used the item can include a portion and/or the entire item unless specifically stated to the contrary.

(105) Although a few preferred embodiments have been shown and described, it will be appreciated by those skilled in the art that various changes and modifications might be made without departing from the scope of the disclosure.

(106) Each feature disclosed in this specification may be replaced by alternative features serving the same, equivalent or similar purpose, unless expressly stated otherwise. Thus, unless expressly stated otherwise, each feature disclosed is one example only of a generic series of equivalent or similar features.