WO2015152957A1 - Inverse imaging with magnetic resonance imaging using blipped gradient encoding - Google Patents

Abstract

Systems and methods for magnetic resonance imaging ("MRI") using an inverse imaging ("Inl") data acquisition that incorporates a blipped partition-encoding gradient scheme are provided. The Inl image reconstruction is modified to incorporate additional phase shifting of the acquired data caused by the blipped gradient encoding scheme. The blipped gradient encoding scheme reduces spatial blurring in the acquired data, thereby improving the spatial resolution achievable as compared with other Inl techniques.

Description

INVERSE IMAGING WITH MAGNETIC RESONANCE IMAGING USING

BLIPPED GRADIENT ENCODING

CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] The present application is based on, claims priority to, and incorporates herein by reference in its entirety US provisional Application Serial No. 61 /972,822 filed March 31, 2014 and entitled "SYSTEMS AND METHODS FOR IN VERSE IMAGING WITH MAGN ETIC RESONANCE IMAGING USING BLIPPED GRADIENT ENCODING."

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH

[0002] This invention was made with government support under EB012107 and

MH093765 awarded by the National Institutes of Health. The government has certain rights in the invention.

BACKGROUND OF THE INVENTION

[0003] The field of the invention is systems and methods for magnetic resonance imaging ("MRI"). More particularly, the invention relates to systems and methods for inverse imaging ("Inl") using a blipped gradient encoding scheme for improved spatial resolution.

[0004] The acquisition time for whole-brain MR] data is limited by the time required to traverse k-space. Data acquisition can be completed by either a set of two- dimensional k-space trajectories across slices, or a single three-dimensional k-space trajectory covering the whole imaging volume.

[000S] Conventional functional MRI ("fMRI") using blood-oxygen-level- dependent ("BOLD") contrast is usually performed using an echo-p!anar imaging ("EPS") pulse sequence. Considering the state-of-the-art gradient slew rate and maximal gradient strength, EPI can complete a two-dimensional k-space traversal in approximately 40 ms per slice, or 2 seconds for whole brain coverage.

[0006] The desire for higher temporal resolution in fMRI has been motivated by its potential to detect and to suppress physiological fluctuations in order to increase the sensitivity of detecting brain activity. Likewise, high-speed echo-volumar image ("EVI") integrated with parallel imaging could resolve the physiological signal fluctuation to increase the sensitivity in event-related fMRI. In addition, higher temporal resolution of fMRI has been suggested to be beneficial for improving the power of detecting neuronal timing information and connectivity among brain areas. [0007] One approach to accelerating fMRI data acquisition is to optimize the k- space trajectory and the corresponding reconstruction method. This can be achieved by using partial k-space sampling, compressed sensing, or exploiting a priori information, such as by using key-hole imaging or and singu lar-value-decomposition MRI.

[0008] As the technology of radio-frequency (RF) receiver coil array advances, parallel MRI methods, which simultaneously acquire MRI data from multiple channels of an RF coil array, have become a method of reducing the scanning time. Parallel MRI methods, such as the k-space SMASH and GRAPPA, or the image domain SENSE method, reduce data acquisition times by reducing the k-space traversal at a cost of reduced signal-to-noise ratio ("SNR"). In fMRI, paraiiei MR] has been successfully combined with gradient-echo EPI to achieve accelerated acquisitions. It has also been demonstrated that incorporating static a priori information can further improve the sensitivity of fMRI.

[0009] The inverse imaging ("Inl") method is a further generalized parallel MRI method for three-dimensional volumetric acquisition that includes leaving out all partition-encoding steps. Consequently, the volumetric brain is projected along the partition-encoding direction onto a single plane. Inl reconstructs a three-dimensional image from a set of two-dimensional projection images from different channels of an RF coil array using the coil sensitivity information. Mathematically, the image reconstruction is performed by solving a set of underdetermined linear systems. Combined with an echo shifting technique, the sampling rate of whole-brain Inl can become as high as 40 Hz.

[0010] While Inl allows for a very high temporal resolution, the attainable spatial resolution depends on the available spatial information in the RF coil array. Correlated coil patial information will cause spatial blurring in the Inl reconstruction. Previous attempts at improving the spatial resolution of Inl have included using sophisticated reconstruction algorithms, such as reconstructing the images in k-space, or using spatial filtering. Another strategy is to modify the data acquisition by collecting data from multiple projections instead of one single projection. These approaches are not without their own limitations, however.

[001.1] It is therefore desirable to provide systems and methods for MRI that are capable of achieving high temporal resolution and high spatial resolution for applications, such as fMRI, where it is preferable to rapidly image a large volume-of- interest, such as the human brain.

SUMMARY OF THE INVENTION

[0012] The present invention overcomes the aforementioned drawbacks by providing systems and methods for improving the spatial resolution achievable with inverse imaging ("Inl") techniques using a magnetic resonance imaging ("M I") system, [0013] it is an aspect of the invention to provide a method for producing a plurality of images of a subject with an MRI system. The method includes directing the MRI system to perform a pulse sequence that includes applying a radio frequency (RF) excitation field to a volume-of-interest; establishing at least one magnetic field gradient along a frequency-encoding direction following the application of the RF excitation field in order to form echo signals; sequentially producing a plurality of magnetic field gradient blips along a partition-encoding direction while the at least one magnetic field gradient is established; and acquiring, with an array of RF receiver coils, data indicative of the formed echo signals. A plurality of projection images are reconstructed from the acquired data. Each projection image includes a superposition of partitions in the volume-of-interest resulting from not establishing a partition-encoding gradient before the frequency-encoding gradient. The plurality of images of the subject are reconstructed by applying an inverse operator to the plurality of projection images. The inverse operator is based at least in part on a spatial encoding provided by the plurality of magnetic field gradient blips.

[0014] The foregoing and other aspects and advantages of the invention will appear from the following description. In the description, reference is made to the accompanying drawings that form a part hereof, and in which there is shown by way of illustration a preferred embodiment of the invention. Such embodiment does not necessarily represent the full scope of the invention, however, and reference is made therefore to the claims and herein for interpreting the scope of the invention.

BRIEF DESCRIPTION OF THE DRAWINGS

[0015] FIG. 1A is an example of a pulse sequence that includes applying a series of magnetic field gradient blips along a partition-encoding direction concurrently with phase-encoding gradient blips and a frequency-encoding gradient; [0016] FIG. IB is another example of a pulse sequence that includes applying a series of magnetic field gradient blips along a partition-encoding direction concurrently with phase-encoding gradient blips and a frequency-encoding gradient;

[0017] FIGS. 2A and 2B illustrate an example of the principles behind the blipped inverse imaging technique described here;

[0018] FIG. 3 is a flowchart setting forth the steps of an example of a method for performing blipped inverse imaging; and

[0019] FIG. 4 is a block diagram of an example of a magnetic resonance imaging

("MR I") system.

DETAILED DESCRIPTION OF THE INVENTION

[0020] Described here are systems and methods for magnetic resonance imaging

("MRI") using an inverse imaging ("Inl") data acquisition and image reconstruction technique that is modified to incorporate additional phase shifting of the acquired data implemented using a blipped gradient encoding scheme. The systems and methods of the present invention employ a radio frequency ("RF") coil array with separate coil elements positioned at different locations relative to a subject positioned in the field-of- vievv ("FQV"). Each coil element receives magnetic resonance signals that are separately amplified, digitized, and processed according to the methods for reconstructing an image described below.

[0021] Using the systems and methods described here, the spatial resolution of

Inl can be improved by modifying the spatial encoding used in the pulse sequence. Because the systems and methods of the present invention utilize a blipped partition- encoding gradient scheme in an Inl framework, they may be referred to as blipped-Inl, or simply "talnl." This blipped gradient encoding scheme is applied to an Inl acquisition in order to reduce spatial blurring, thereby increasing the attainable spatial resolution.

[0022] The succeeding description is provided with reference to the following orthogonal spatial encoding directions common to MRI: a slice-encoding direction or partition-encoding direction; a phase-encoding direction; and a frequency-encoding direction. By way of example, as referred to herein, the partition-encoding direction corresponds to the z-direction in the image domain, which is aligned along the longitudinal axis of the bore of an exemplary MRI system, and the /i, -direction in k- space. Accordingly, as referred to herein, the phase-encoding direction corresponds to the y-direction in the image domain, and the k_v -direction in k-space; and the frequency-encoding direction corresponds to the x-direction in the image domain, and the k_x -direction in k-space, it will be appreciated by those skilled in the art that the choice of these directions is arbitrary and any suitable permutation of these directions, or any set of orthogonal oblique directions, is possible and within the scope of the present invention.

[0023] An example of a pulse sequence that can be implemented to direct an MR! system to acquire image data in accordance with some embodiments of the present invention is illustrated in FIG. 1A. This example pulse sequence is similar to a single- slice two-dimensional echo planar imaging ("EPi") acquisition, except for the additional encoding gradient blips played out along the G_ -axis and the slab-selective RF pulse. As described below, the G, blips are played out synchronously with phase-encoding gradient blips in order to provide extra spatial encoding along the z-axis,

[0024] The pulse sequence includes a spatially selective radio frequency ("RF") excitation pulse 100 that is played out in the presence of a slab-selective gradient 102 in order to produce transverse magnetization in an imaging volume, which may contain a plurality of different imaging slices, or partitions. The slab-selective gradient 102 includes a repbasmg lobe 104 that acts to rephase unwanted phase dispersions introduced by the slab-selective gradient 102 such that signal losses resultant from these phase dispersions are .mitigated.

[0025] Following excitation of the nuclear spins in the prescribed imaging volume, k-space data is acquired by sampling a series of gradient-recalled echo signals in the presence of an alternating readout gradient 106. The alternating readout gradient is preceded by the application of a pre-winding gradient 108 that acts to move the first sampling point along the frequency-encoding, or readout, direction by a distance &k_x in k-space. Spatial encoding of the echo signals along a phase-encoding direction is performed by a series of phase encoding gradient "blips" 110, which are each played out in between the successive signals readouts such that each echo signal is separately phase encoded. The phase encoding gradient blips 110 are preceded by the application of a pre-winding gradient 112 that acts to move the first sampling point along the phase-encoding direction a distance Ak„ in k-space. Together, the pre- winding gradients 108 and 112 serve to begin the sampling of k-space at a defined k- space location [ k_x,k ] .

[0026] During the application of each phase encoding gradient blip 110, a magnetic field gradient blip is also played out along the partition-encoding direction. These partition-encoding gradient blips 114 act to impart phase shifts on the acquired data, as will described in more detail below. Each successive slice-encoding gradient blip 114 is equal in magnitude and opposite in polarity than the preceding blip. In this manner, the partition-encoding gradient blips 114 do not produce phase accumulations in the phase-encoding direction in k-space because each successive blip 114 serves to unwind the phase accrued by the previous blip 114.

[0027] A variation of the pulse sequence in FIG. 1A is illustrated in FIG. IB. The pulse sequence of FIG. 1A achieves an in-plane shift of FQV72, whereas the pulse sequence of FIG. IB achieves an ΐη-p!ane shift of FOV/3. The gradient moment of the G, blips in the bin! pulse sequence can be expressed as, γ· O _z

[0028] where β denotes a real-number scale factor, γ denotes the gyromagnetic ratio, FOV_z denotes the length along partition encoding direction, and Ak_z is the minimum spacing in k-space along the &_z -direction.

[0029] For the simultaneous muitislice acquisition, the G_ blip encoding creates an inter-slice image shift along the phase-encoding direction between simultaneously excited slices, as described in U.S. Patent No. 8,405,395, which is herein incorporated by reference in its entirety. The effect of this blip-encoding on a volumetric bin! acquisition, however, is different from the multi-slice acquisition.

[0030] FIGS. 2A and 2B show the principle of blnl encoding for the case of an in- plane FOV/2 shift and of 1. The dashed lines in the left-most image in FIG. 2A delineate image partitions (obtained after all partition encoding steps). Because the G_z blips always have zero gradient moments in odd phase encoding lines, the G_z blips only introduce phase shifts in even phase encoding lines. Moreover, different partitions receive different phase offsets introduced by the G_z blips. As an example, the partition S2 is at the scanner iso-center and thus has no phase offset

[0031] On the contrary, Gz blips introduce -7τ and +ff/4 at partitions 51 and

S3, respectively, which are at the position 3/8 and 5/8 of the FOV in the partition- encoding direction. These phase modulations on different phase encoding lines of the different partitions are marked at the right margins of the k-space diagrams for partitions SI, S2, and S3 in FIG. 2A. In image space, these k-space phase modulations cause spatial shifts of FOV/2 along the phase encoding direction for all partitions, but with different weighting for different partitions, as shown in FIG 2B.

[0032] The top row in FIG. 2B shows the images in representative partitions without the phase modulation introduced by G . blips; the middle row shows the phase offsets introduced by G. blips; and the bottom row shows the modulated slice images.

Using such G„ blips, strong /2 ghosting are observed in the partitions toward the edge of the excitation volume, while central slices show relatively weak N /2 ghosting. In accelerated bin! acquisition, all the partition-encoding steps are removed and consequently all the partitions are integrated, as shown in the right-most images in FIG. 2B. it is noted that a similar analysis can be applied to other types of shifts, such as a FOV/3 shift.

[0033] The general Inl acquisition and reconstruction method, upon which the present invention is developed, is described in U.S. Patent No. 7,394,251, which is herein incorporated by reference in its entirety. In the previous Inl method, each vector, which includes image pixels of the same frequency/phase encoding but different partition encoding, is integrated into one projection voxel. These projection voxels are independent to each other and can be reconstructed independently. In blnl, however, the projection voxels are mutually dependent because of the different spatial shifts introduced along the phase encoding direction in different partitions. Because of this mutual dependence, a single voxel in 3D is related to multiple phase encoding and partition encoding positions in the accelerated projection image. As a result, the spatial reconstruction implemented in bin! accounts for these correlations.

[0034] For a projection image at time instant, t , all phase encoding image voxels are cascaded with their frequency encoding index X_; i < Χ_Ί≤ N_x ) from all channels in a coil array into one f N_y X n_c 1 -by-1 vector, y (t) , where N_y denotes the number of image pixels along the phase encoding direction and n„ denotes the number of channels in an RF coil array. Each image plane corresponding to the same frequency encoding indices, but different phase encoding and partition encoding indices, and can be independently reconstructed. To simplify the notation, the explicit x-coordinate dependency (χ_Ί) is thus left out in the following description.

[0035] The vector, y(t) , is linearly related to the image vector to be reconstructed,

i i ) + n ( 3; where s{ tj is a j_^N, X N₇ j -by-1 image vector to be reconstructed, and is a j N_v X N_C j-tay- 1 vector denoting the contamination noise. The matrix, A, is a

N.. xN by- N_v X N. j forward matrix that maps the signals from the phase- partition encoding plane (here, the y-z plane) to one vector of observed signal in the accelerated image across n_c -channel in the coil array.

[0037] The forward matrix, A , includes three components: spatial integration along the partition encoding direction due to leaving out partition encoding steps in the accelerated acquisition (projection matrix), spatial encoding due to G₇ blips, and RF coil sensitivity. Let A_f. denote a part of the forward matrix related to channel C with 1 < i≤ n_c. The matrix, A_c can be written explicitly as a product of three matrices,

A . zzz A^^'A^A ⁰'⁷ (3);

[0038] where A^pmj denotes a spatial projection matrix, A°_c'^ip denotes a G_ -biip encoding matrix, and A™" denotes an RF coii sensitivity matrix. Each of these matrices will now be discussed in turn.

The coil-sensitivity matrix coil, A^('°" , is diagonal with dimensions of

N . X N. -bv- N.. x N. , Diagonal elements of A^c _r ^od can be empirically measured by a reference scan. As an example, the reference scan can be performed using a conventional Ini acquisition without the blipped G_z gradient encoding scheme described above.

[0040] The projection matrix, Α '°³ , describes the projection operation due to omitting all partition encoding steps. This matrix can be formulated as,

where f ->'j is an identity " with dimensions N > .. -by- N . The dimension of

A is .Y..-by- [A _r x A ] .

[0042] The aliasing matrix, A°"^p describes the blip encoding introduced by G. - blips. This matrix can either be derived empirically from inl reference scans, in order to account for possible imperfections in the G_ gradient encoding, or can be derived theoretically from a given G blip encoding scheme. The dimensions of A^" is also

The first step in deriving A^bJ^lp is to calculate the phase modulation in the k-space for each partition. As an example of an empirical approach for deriving the phase modulations imparted by the G_z blip encoding scheme, the phase modulation can be estimated from the full three-dimensional k-space data of bin! and Inl reference scans (with ail partition encoding steps) by subtracting the phase of the even k-space lines between the Inl and bin I data. In some instances, To achieve a good estimate, only the central two even k-space lines may be used because they are likely to have higher SNR.

[0044] As an example of a theoretical approach for deriving the phase modulations imparted by the G_z blip encoding scheme, the phase modulation (at even phase encoding lines) at partition,∑.., can be calculated as,

2 ί

0 = (5)_;

FOV_z

[0045] where d_z denotes the distance away from the isocenter in the partition encoding direction and z. = 1,... ,N_z , which denotes the number of partitions.

[0046] Proceeding with the derivation of the blip-encoding matrix, let k^" and k ^~ni denote N„-by-l data vectors of Inl and blnl in hybrid space at partition encoding index, z_{, respectively. Here, the hybrid space is defined as the three-dimensional space with y-dimension in frequency domain and x- and z-dimensions in image domain. The k^&/ and k^M vectors are related to each other as,

[0047] where Θ_Ζ o k^*"⁷ denotes the element-by-element Hadamard product between the vectors θ_ζ and , and Θ_ζ denotes a vector of complex phase offsets.

For the case of a FOV/2 shift, the odd elements of Θ are equal to one and the even elements of Θ_ζ are equal to ex ( i _z ) . For the case of a FOV/3 shift, the (3/? + l)^r" ,

( 3/7 + 2 )^"" , and ( 3/> ÷ ) elements, with p— 0, 1, 2, . . ., are equal to one, ex ( y#_ , and exp( ]2Θ_Ζ ) , respectively.

[0048] After applying a discrete Fourier transform, Eqn. (6) becomes,

-::: 1ΤΤ{θ } ® p*^'' (7);

[0049] where p^'^J and ^'¹ denote image vectors of blnl and Inl, respectively;

FFT^Q_z j denotes the discrete Fourier transform; and <¾ denotes a circular convolution operator. Mathematically, the circular convolution operation can also be rewritten as a matrix-vector product,

p':^/ - / <| Θ . } <8> p - Aj* · p (8);

[0050] where Α'^'^ψ ₇ denotes the blip encoding matrix at partition z. with dimension of N_y-by- N_v. By simply replacing ¹"¹ with an identity matrix, ϊ_γ , in Eqn. (8), the blip encoding matrix can be derived analytically, as follows:

A* J · = AjJ = /7 ^'{θ . } ® I (9).

[0051] Combining all the Α'^' with z, = i,..., N_T , the blip encoding matrix for channel c is,

Taking A^P"'^J , A^Mp , A ^o" , and Eqn. (3) together, the dimension of A,, is

JV„ X JV, j . Finally, the forward matrix, A, is a vertical concatenation of A,. with dimensions of N , X w, -bv- N X N ¹

[0053] As an example, the random variability and inter-channel correlations of multi- channel projection images can be characterized by a noise covariance matrix, C. Using the noise covariance, the forward matrix, A, can be whitened to equalize the sensitivity of all channels. For each frequency encoding index, %. , the image series can be formatted as an image matrix with dimension of T denotes

the number of time points.

[0054] By calculating the sample covariance of the resting state image matrix, the noise covariance, C , is a ! N^ X W_cj -by- ^N_y X n_c j matrix. Singular value decomposition ("SVD^") decomposes the noise covariance matrix,

C = UAU" ^' = (c )(c ⁱ ) (12);

[0055] where U is a complex unitary matrix and Λ is a non-negative real-valued rectangular diagonal matrix. The superscript, H , indicates the complex conjugate and transpose. By using,

C^!/2 - UA^V2 ( 13);

[0056] The spatially whitened signal equation becomes,

y^w (i) = A^ws(i ) + iT (i) (14);

[0057] where,

y^w (/) - C^'-^¾/2y i» (15);

A^,v = tr^,/2A (16); and

where ^■ · ·) represents an ensemble average operation identity matrix of size N„xn, -by- N xn

[0060] Solving for the images, § i ) , in Eqn. (14) is an ill-posed inverse problem when the number of unknowns is larger than the observations, in order to obtain a unique solution, it is preferable to provide additional constraints. One common choice is the minimum-norm estimate ("MNE"}:

% ( / ) - ar mi n { | ^' ( / ) ···· Λ " s ( / ) _, - ?: ||s ( / )||^ j (18);

[0061] where jj- · -| is the square of the ² -norm and ?< ^' is a regularization parameter. Although MNE is employed to reconstruct the spatial information, it will be appreciated that other suitable reconstruction methods can also be employed. For instance, more sophisticated reconstruction methods, such as linear constrained minimum variance ("LCMV") (Lin et a!,, 2008) and k-space-reconstruction Inl (Lin et al., 2010) may be employed to bring further spatial resolution improvement.

[0062] The solution to Eqn. (18) can be obtained by a time-invariant linear inverse operator, W, as follows:

A- j A- R i / ^' + Λ , . _γ ) ^' (19); j where R is a source covariance matrix. If no other spatial prior information is to be incorporated into R, then,

R = I,v yx_;V ^z (20);

[0064] corresponding to the assumption that a spatially uniform image contributes to the measurement. The regularization parameter, A² , can be calculated

[006S] where C ^" denotes the inverse of the pre-defined SNR of the whitened data and tr( - - -) represents the trace of a matrix. Eqn. (19) can be applied to reconstruct a y-z image plane over all time points according to,

WV (/) - R ( /V ) ( A^WR(A^W ) ^/ + J ¹ y^w (/ } (22);

Repeating this inversion procedure over y(f) of the different frequency encoding indices yields the volumetric reconstruction of dynamic images. The time series of the reconstructed brain can be reshaped to a matrix, S, with dimensions of N

-by- N x ^r _x N_z j , where N_t denotes the number of time instants.

[0067] Referring now to FIG. 3, a flowchart setting forth the steps of an example of a method for reconstructing images for use in a functional MRi ("fMRl") applications. Images are acquired as a series of image frames using a volumetric biipped-Inl scan; thus, the method begins with the acquisition of data, as indicated at step 302. As an example, this data acquisition can be accomplished using the pulse sequence of FIG. 1A or FIG IB.

[0068] As indicated at step 306, a time series of complex-valued projection images is reconstructed from the acquired data using a standard two-dimensional Fourier transform. These images are stored for later use in the reconstruction. Optional preprocessing can be applied to the acquired data, as indicated at step 306. In one example, preprocessing may include suppressing physiological noise. For instance, an algorithm, such as the DRIFTER algorithm described by Sarkka, et al., 2012, can be used to estimate and separate the cardiac (0.75 Hz - 2.5Hz) and respiratory (0.083 Hz - 0.5 Hz) fluctuations from the fMRl signals.

[0069] After the optional preprocessing, volumetric spatial reconstruction is performed, as generally indicated at 308. This reconstruction includes the application of an inverse operator to the projection images reconstructed in step 304, and as generally described above with respect to Eqn. (22). To this end, an inverse operator is computed, as indicated at step 310, and described above. For instance, computing the inverse operator may include forming a spatial projection matrix, a blip encoding matrix, and an RF coil sensitivity matrix, as indicated at steps 31.2, 314, and 316, respectively. Computing the coil sensitivity matrix may include calculating coil sensitivity maps or provided previously calculated coil sensitivity maps. As an example, the spatial sensitivity patterns of a coil can be estimated from low resolution images (.magnitude and phase) acquired with minimal tissue contrast. These matrices can be used to compute the inverse operator as described above. As an example, the inverse operator, W , can be computed using Eqns. (11), (16), and (19)-(21), or suitable combinations thereof.

[0070] After the images have been reconstructed., they can be stored or displayed to a user, as indicated at step 318. In some embodiments, additional processing of the images may be desired, such as in fMRl applications where it may be desirable to produce neuronal activation maps from the reconstructed images.

[0071] As an example of additional processing consistent with fMRi applications, for each voxel in the reconstructed 3D volume, the hemodynamic response ("HDR") can be separately estimated by using the standard general linear mode! ("GLM") framework (Friston et al., 1995a, Friston et a!., 1.995b, Friston et al., 1995c). A finite impulse response ("FIR") basis functions can be used to model the HDR to avoid bias in estimating the shape of HDRs. By way of example, the FIR bases can be temporally synchronized to the onsets of the stimuli, spanning over a 30 second period, including a 6 second pre-stimulus baseline and a 24 second post-stimulus interval. Because Inl and bin! data are sampled at 10 Hz, the FIR bases in this example include 300 temporally shifted discrete delta functions. The GLM can also model confounds of temporal constants, a linear trend, and low-frequency drifts. The coefficient for each FIR basis is calculated by generalized least squares.

[0072] More specifically for this example, let s denote the r^th column vector of

S and let denote the design matrix of the GLM. The estimated HDR basis coefficients are, h_r = (X^7'x)^"] X⁷ S_r (23); [0073] where l> denotes estimated HDR coefficients.

[0074] To obtain statistical inferences, the noise level of the reconstructed volumetric images can be estimated from the pre-stimulus baseline. Using these noise estimates, dynamic statistical parametric maps ("dSPM") can be derived as the ratio of the MNE values over the baseline noise estimates at each voxel. The dSPM follows a t- like distribution under the null hypothesis of no hemodynamic response (Dale et al, 2000). Since the number of samples used to calculate the noise variance is high, the t- distribution approaches the norma! distribution and the t-statistics approximates the z- scores.

[0075] Described here is a refined version of the inverse imaging ("Inl") technique to achieve fast sampling rate (e.g., 10 Hz with whole brain coverage) in fMRI experiments with improved spatial resolution, inl is fast because of the removal of partition encoding, but has relatively low spatial resolution toward the center of the brain due to spatial information in channels of a coil array that is insufticien to separate between partitions. Additional G_z gradient encoding blips are incorporated into the Inl acquisition to improve the forward matrix in the image reconstruction, which results in a higher spatial resolution and localization accuracy.

[0076] Blipped Ini employs whole brain volume excitation and acquisition, and is therefore a type of single-shot 3D imaging. In the 3D k-space view, blnl has a zig-zag k- space trajectory in contrast to the standard inl method, which acquired data on a single k-plane trajectory at —0.

[0077] in some embodiments, it may be sufficient to assume that the phase encodings of the G, gradient blips are constant across the voxel dimension along the image partition. However, additional linear phase encoding components may arise due to the finite thickness, Δζ , of the image partition. This finite thickness may cause through-partition intra-voxel dephasing, similar to the through-slice dephasing in the biipped-CAIPi acquisition. This intra-voxel dephasing results in signal attenuation, which can be calculated as,

[0078] Stronger blip gradient moments, that is, those with a larger β , will cause more signal attenuation. For the case of β— 4 and N, = 64 , the signal attenuation is around 0.64%. Therefore, the through-partition signal attenuation is very much negligible when the gradient blip strength is not stronger than 4Ak_ .

[0079] it is contemplated that blnl can be used to reduce spatial blurring and localization error of the previously described inl method. The performance of b!nl is expected to improve further with larger gradient blip sizes, such as 4Ak_z . The spatiotemporal pattern of BOLD activation in blnl has been shown to be consistent with the canonical BOLD response. Given its demonstrated benefits, it is contemplated that blnl can be a useful tool for investigating human brain function in cortical and subcortical areas at high spatiotemporal resolution.

[0080] Referring particularly now to FIG. 4, an example of a magnetic resonance imaging ("MRI"j system 400 is illustrated. The MRi system 400 includes an operator workstation 402, which will typically include a display 404; one or more input devices 406, such as a keyboard and mouse; and a processor 408. The processor 408 may include a commercially available programmable machine running a commercially available operating system. The operator workstation 402 provides the operator interface that enables scan prescriptions to be entered into the MR! system 400. In general, the operator workstation 402 may be coupled to four servers: a pulse sequence server 410; a data acquisition server 412; a data processing server 414; and a data store server 416. The operator workstation 402 and each server 410, 412, 414, and 416 are connected to communicate with each other. For example, the servers 41.0, 412, 414, and 416 may be connected via a communication system 440, which may include any suitable network connection, whether wired, wireless, or a combination of both. As an example, the communication system 440 may include both proprietary or dedicated networks, as well as open networks, such as the internet

[0081] The pulse sequence server 410 functions in response to instructions downloaded from the operator workstation 402 to operate a gradient system 418 and a radiofrequency ("RF") system 420. Gradient waveforms necessary to perform the prescribed scan are produced and applied to the gradient system 418, which excites gradient coils in an assembly 422 to produce the magnetic field gradients G_y, G , and

G_z used for position encoding magnetic resonance signals. The gradient coil assembly

422 forms part of a magnet assembly 424 that includes a polarizing magnet 426 and a whole-body RF coil 428.

[0082] RF waveforms are applied by the RF system 420 to the RF coil 428, or a separate local coil (not shown in FIG. 4), in order to perform the prescribed magnetic resonance pulse sequence. Responsive magnetic resonance signals detected by the RF coil 428, or a separate local coil (not shown in F3G. 4), are received by the RF system 420, where they are amplified, demodulated, filtered, and digitized under direction of commands produced by the pulse sequence server 410. The RF system 420 includes an RF transmitter for producing a wide variety of RF pulses used in MRI pulse sequences. The RF transmitter is responsive to the scan prescription and direction from the pulse sequence server 410 to produce RF pulses of the desired frequency, phase, and pulse amplitude waveform. The generated RF pulses may be applied to the whole-body RF coil 428 or to one or more local coils or coil arrays (not shown in FIG. 4).

[0083] The RF system 420 also includes one or more RF receiver channels. Each

RF receiver channel includes an RF preamplifier that amplifies the magnetic resonance signal received by the coil 428 to which it is connected, and a detector that detects and digitizes the / and Q quadrature components of the received magnetic resonance signal The magnitude of the received magnetic resonance signal may, therefore, be determined at any sampled point by the square root of the sum of the squares of the / and Q components: = ^/² + β² (25);

[0084] and the phase of the received magnetic resonance signal may also be determined according to the following relationship:

[0085] The pulse sequence server 410 also optionally receives patient data from a physiological acquisition controller 430. By way of example, the physiological acquisition controller 430 may receive signals from a number of different sensors connected to the patient, such as electrocardiograph ("EGG") signals from electrodes, or respirator}' signals from a respiratory bellows or other respiratory monitoring device. Such signals are typically used by the pulse sequence server 410 to synchronize, or "gate," the performance of the scan with the subject's heart beat or respiration.

[0086] The pulse sequence server 410 also connects to a scan room interface circuit 432 that receives signals from various sensors associated with the condition of the patient and the magnet system. It is also through the scan room interface circuit 432 that a patient positioning system 434 receives commands to move the patient to desired positions during the scan.

[0087] The digitized magnetic resonance signal samples produced by the RF system 420 are received b the data acquisition server 412. The data acquisition server 412 operates in response to instructions downloaded from the operator workstation 402 to receive the real-time magnetic resonance data and provide buffer storage, such that no data is lost by data overrun, in some scans, the data acquisition server 412 does little more than pass the acquired magnetic resonance data to the data processor server 414. However, in scans that require information derived from acquired magnetic resonance data to control the further performance of the scan, the data acquisition server 412 is programmed to produce such information and convey it to the pulse sequence server 410. For example, during prescans, magnetic resonance data is acquired and used to calibrate the pulse sequence performed by the pulse sequence server 410. As another example, navigator signals may be acquired and used to adjust the operating parameters of the RF system 420 or the gradient system 418, or to control the view order in which k-space is sampled. In still another example, the data acquisition server 412 may also be employed to process magnetic resonance signals used to detect the arrival of a contrast agent in a magnetic resonance angiography ("MR A"] scan. By way of example, the data acquisition server 412 acquires magnetic resonance data and processes it in real-time to produce information that is used to control the scan.

[0088] The data processing server 414 receives magnetic resonance data from the data acquisition server 412 and processes it in accordance with instructions downloaded from the operator workstation 402. Such processing may, for example, include one or more of the following: reconstructing two-dimensional or three- dimensional images by performing a Fourier transformation of raw k-space data; performing other image reconstruction algorithms, such as iterative or backprojection reconstruction algorithms; applying filters to raw k-space data or to reconstructed images; generating functional magnetic resonance images; calculating motion or flow images; and so on.

[0089] images reconstructed by the data processing server 414 are conveyed back to the operator workstation 402 where they are stored. Real-time images are stored in a data base memory cache (not shown in FIG. 4), from which they may be output to operator display 412 or a display 436 that is located near the magnet assembly 424 for use by attending physicians. Batch mode images or selected real time images are stored in a host database on disc storage 438. When such images have been reconstructed and transferred to storage, the data processing server 414 notifies the data store server 416 on the operator workstation 402. The operator workstation 402 may be used by an operator to archive the images, produce films, or send the images via a network to other facilities.

[0090] The MR1 system 400 may also include one or more networked workstations 442. By way of example, a networked workstation 442 may include a display 444; one or more input devices 446, such as a keyboard and mouse; and a processor 448. The networked workstation 442 may be located within the same facility as the operator workstation 402, or in a different facility, such as a different healthcare institution or clinic.

[0091] The networked workstation 442, whether within the same facility or in a different facility as the operator workstation 402, may gain remote access to the data processing server 414 or data store server 416 via the communication system 440. Accordingly, multiple networked workstations 442 may have access to the data processing server 414 and the data store server 416. in this manner, magnetic resonance data, reconstructed images, or other data may exchanged between the data processing server 414 or the data store server 416 and the networked workstations 442, such that the data or images may be remotely processed by a networked workstation 442. This data may he exchanged in any suitable format, such as in accordance with the transmission control protocol ("TCP"), the internet protocol (^""IP"), or other known or suitable protocols.

[0092] The present invention has been described in terms of one or more preferred embodiments, and it should be appreciated that many equivalents, alternatives, variations, and modifications, aside from those expressly stated, are possible and within the scope of the invention.

Claims

1. A method for producing a plurality of images of a subject with a magnetic resonance imaging (MR!) system, the steps of the method comprising:

a) directing the R3 system to perform a pulse sequence that includes:

applying a radio frequency (RF) excitation field to a vo!ume-of-inferest; establishing, with the MRI system, at least one magnetic field gradient along a frequency-encoding direction following the application of the

RF excitation field in order to form echo signals;

sequentially producing a plurality of magnetic field gradient blips along a partition-encoding direction while the at least one magnetic field gradient is established;

acquiring with an array of RF receiver coils, data indicative of the formed echo signals

b) reconstructing a plurality of projection images from the acquired data, each projection image comprising a superposition of partitions in the volume-of- interest resulting from not establishing a partition-encoding gradient before the frequency-encoding gradient;

c) reconstructing the plurality of images of the subject by applying an inverse operator to the plurality of projection images, the inverse operator being based at least in part on a spatial encoding provided by the plurality of magnetic field gradient blips.

2. The method as recited in claim 1, wherein step c) includes calculating the inverse operator based on a blip-encoding matrix that describes the spatial encoding provided by the plurality of magnetic field gradient blips; a coil sensitivity matrix that describes coil sensitivities of the RF receiver coils; and a projection matrix that describes for each projection image, a spatial integration of the image slices in the volume-of-interest.

3. The method as recited in claim 2, wherein the blip-encoding matrix is formed by calculating phase modulations in the acquired data associated with each partition.

4. The method as recited in claim 3, wherein the phase modulations associated with a given partition are calculated based on a distance of the partition from an isocenter of the MR I system along the partition-encoding direction, and a size of a field-of-view along the partition-encoding direction.

5. The method as recited in claim 3, wherein the phase modulations associated with a given partition are estimated based on the acquired data and on reference data acquired without sequentially producing the plurality of magnetic field gradient blips along the partition-encoding direction.

6. The method as recited in claim 5, wherein estimating the phase modulations includes subtracting a phase of even k-space lines in the acquired data from a phase of even k-space lines in the reference data.

7. The method as recited in claim 3, wherein the blip-encoding matrix is formed as a circular convolution between an identity matrix and a Fourier transform of the phase modulations.

8. The method as recited in claim 2, wherein the coil sensitivity matrix is formed based on calculating a coil sensitivity map for each of the plurality of RF receiver coils.

9. The method as recited in claim 2, wherein the projection matrix is formed as a concatenation of identity matrices.

10. The method as recited in claim 2, wherein calculating the inverse operator includes:

forming a forward matrix using the blip-encoding .matrix, the coil sensitivity matrix, and the projection matrix;

calculating a noise covariance matrix among the plurality of RF receiver coils; calculating a spatially whitened forward operator using the forward .matrix and the covariance matrix; and

producing the inverse operator from the spatially whitened forward operator.

11. The method as recited in claim 2, wherein calculating the inverse operator includes:

forming a forward matrix using the blip-encoding matrix, the coil sensitivity matrix, and the projection matrix;

calculating a source covariance matrix using known information about the subject being imaged; and

wherein the inverse operator is produced using the source covariance matrix and the forward matrix.

12. The method as recited in claim 11, wherein calculating the inverse operator includes spatially weightening the forward matrix by calculating a noise covariance matrix among the plurality of RF receiver coils and applying the covariance matrix to the forward matrix.

13, The method as recited in claim 11, wherein calculating the inverse operator includes calculating a regularization parameter,

14, The method as recited in claim 13, wherein the regularization parameter is based at least in part on the forward matrix and the source covariance .matrix.