SYSTEM AND METHOD FOR CONVOLUTION OPERATIONS FOR DATA ESTIMATION FROM COVARIANCE IN MAGNETIC RESONANCE IMAGING
20170336488 · 2017-11-23
Inventors
Cpc classification
G01R33/5611
PHYSICS
G01R33/5608
PHYSICS
International classification
G01R33/561
PHYSICS
G01R33/56
PHYSICS
Abstract
Described here are systems and methods for reconstructing images of a subject using a magnetic resonance imaging (“MRI”) system. As part of the reconstruction, synthesized data are estimated at arbitrarily specified k-space locations from measured data at known k-space locations. In general, the synthesized data is estimated using a convolution operation that is based on measured or estimated covariances in the acquired data. The systems and methods described here can thus be referred to as Convolution Operations for Data Estimation from Covariance (“CODEC”).
Claims
1. A method for producing an image of a subject using a magnetic resonance imaging (MRI) system, the steps of the method comprising: (a) directing the MRI system to acquire data from the subject using multiple receive coils; (b) producing covariance maps that define covariances of the acquired data as a function of differences of k-space locations between pairs of the multiple receive coils; (c) estimating intermediate data using the acquired data and the produced covariance maps; (d) selecting a desired k-space sampling pattern; (e) producing synthesized data on the selected k-space sampling pattern by convolving the intermediate data with the covariance maps; and (f) reconstructing an image of the subject from at least the synthesized data.
2. The method as recited in claim 1, wherein step (f) includes producing combined data by combining the synthesized data with the acquired data and reconstructing the image of the subject from the combined data.
3. The method as recited in claim 1, wherein step (b) includes producing the covariance maps from the acquired data.
4. The method as recited in claim 3, wherein step (b) includes reconstructing low-resolution images from the acquired data and estimating the covariance maps from the low-resolution images.
5. The method as recited in claim 4, wherein estimating the covariance maps includes Fourier transforming a product calculated by multiplying a low-resolution image associated with a first receive coil with a conjugate of a low-resolution image associated with a second receive coil.
6. The method as recited in claim 1, wherein step (b) includes producing the covariance maps from coil sensitivity maps associated with the multiple receive coils.
7. The method as recited in claim 1, wherein step (b) includes directing the MRI system to acquire additional data, and the covariance maps are produced from the additional data.
8. The method as recited in claim 1, wherein step (c) includes estimating the intermediate data according to, x=R.sub.xxz, wherein x is the acquired data, R.sub.xx is a matrix defined by the covariance maps, and z is the intermediate data.
9. The method as recited in claim 8, wherein the intermediate data are estimated using a conjugate gradient method.
10. The method as recited in claim 1, wherein the acquired data are acquired using a non-Cartesian acquisition.
11. The method as recited in claim 10, wherein the non-Cartesian acquisition includes sampling k-space using a variable density spiral trajectory.
12. The method as recited in claim 11, wherein the k-space sampling pattern selected in step (d) includes spiral trajectories that are interleaved with the variable density spiral trajectory used to acquire data in step (a).
13. The method as recited in claim 1, wherein the k-space sampling pattern selected in step (d) includes a Cartesian grid pattern.
14. The method as recited in claim 1, wherein the synthesized data are estimated in step (e) for different coils than the acquired data.
15. A magnetic resonance imaging (MRI) system comprising: a magnet system configured to generate a polarizing magnetic field about at least a portion of a subject arranged in the MRI system; a plurality of gradient coils configured to apply a gradient field to the polarizing magnetic field; a radio frequency (RF) system configured to apply an excitation field to the subject and acquire MR image data from a ROI; a computer system programmed to: control the plurality of gradient coils and RF system to acquire data from the subject using multiple receive coils; produce covariance maps that define covariances of the acquired data as a function of differences of k-space locations between pairs of the multiple receive coils; estimate intermediate data using the acquired data and the produced covariance maps; select a desired k-space sampling pattern; produce synthesized data on the selected k-space sampling pattern by convolving the intermediate data with the covariance maps; and reconstruct an image of the subject from at least the synthesized data.
16. The system as recited in claim 15, wherein the computer system is further programmed to combine the synthesized data with the acquired data to form combined data and reconstruct the image of the subject from the combined data.
17. The system as recited in claim 15, wherein the computer system is further programmed to produce the covariance maps from the acquired data, reconstruct low-resolution images from the acquired data, estimate the covariance maps from the low-resolution images.
18. The system as recited in claim 17, wherein to estimate the covariance maps, the computer system is further programmed to perform a Fourier transform on a product calculated by multiplying a low-resolution image associated with a first receive coil with a conjugate of a low-resolution image associated with a second receive coil.
19. The system as recited in claim 15, wherein the computer system is further programmed to produce the covariance maps from coil sensitivity maps associated with the multiple receive coils.
20. The system as recited in claim 15, wherein the computer is further programmed to estimate the intermediate data according to, x=R.sub.xxz, wherein x is the acquired data, R.sub.xx is a matrix defined by the covariance maps, and z is the intermediate data.
21. The system as recited in claim 15, wherein the computer is further programmed to use a non-Cartesian acquisition to acquire the acquired data.
22. The system as recited in claim 15, wherein the computer is further programmed to control the plurality of gradient coils and RF system to acquire data from the subject using multiple receive coils using a parallel imaging process.
Description
BRIEF DESCRIPTION OF THE DRAWINGS
[0015]
[0016]
DETAILED DESCRIPTION
[0017] Described here are systems and methods for reconstructing images of a subject using a magnetic resonance imaging (“MRI”) system. As part of the reconstruction, data, y.sub.m, can be estimated or otherwise synthesized at arbitrarily specified k-space locations from measured data, x.sub.n, at known k-space locations. In general, the synthesized data are estimated using a convolution operation that is based on measured or estimated covariances in the acquired data. The systems and methods described here can thus be referred to as Convolution Operations for Data Estimation from Covariance (“CODEC”).
[0018] Both the measured and synthesized data may correspond to one or more receive coils, either real or virtual. This method can be applied in a straightforward manner to data collected from, or synthesized onto, any type of trajectory. One general advantage of this class of methods is in reducing scan time by collecting a fraction of the desired data and estimating the rest; although, this class of methods could also be used for enforcing data consistency, or other purposes.
[0019] The systems and methods described here are advantageous for data synthesis for non-Cartesian imaging, which is currently achieved by two major approaches: GRAPPA (and its variants) and iterative SENSE. The general advantage of CODEC over GRAPPA is that CODEC does not require a fixed geometry between sampled and undersampled points and does not require training. The general advantage of CODEC over iterative SENSE is speed, as the iterations in the CG approach for iterative SENSE require gridding, FFT, IFFT, and degridding, while for CODEC require only a single convolution. It is contemplated that for 3D non-Cartesian imaging in particular, CODEC will be substantially faster and require far less memory than iterative SENSE.
[0020] The estimate of synthesized data can be framed as follows. Having sampled points x.sub.1, . . . , x.sub.n in k-space, it is desired to estimate a different point, y.sub.1, with a linear combination, s, of the sampled points, x.sub.1, . . . , x.sub.n. First, covariance can be defined as,
r.sub.mn=E[x.sub.mx.sub.n*] (1);
[0021] variance can be defined as,
σ.sub.n.sup.2=r.sup.m (2);
[0022] and the correlation coefficient can be defined as,
[0023] A conditional estimation based on correlation can be defined as follows:
[0024] Letting,
[0025] where c and x are vectors of the coefficients of c.sub.n and x.sub.n, respectively, the correlation of s with y.sub.1, which it is desired to maximize, can be calculated as,
[0026] where r.sub.yn is a vector of the covariances of y.sub.1 and x.sub.n, and R.sub.xx is the covariance matrix of the measured data, x. The correlation in Eqn. (6) is maximized using,
[0027] where d is introduced as a shortcut for the radicand. Eqns. (4)-(7) provide an estimate of y.sub.1 as,
[0028] Eqn. (8) can be expanded to estimate a vector of points, y, using measured points, x, as follows:
ŷ=R.sub.yxR.sub.xx.sup.−1x (9);
[0029] where R.sub.yx is a covariance matrix of the covariances between each sampled point, x.sub.n, and each point to be estimated, y.sub.m, whose elements {m,n} are the covariance values, E[y.sub.mx.sub.n*]; and R.sub.xx is a covariance matrix of the sampled data whose elements {m,n} are the covariance values, E[x.sub.mx.sub.n*].
[0030] For data collected in a particular scan, these covariances are generally calculated for every pair of coils c.sub.a (for a given y.sub.m) and c.sub.b (for a given x.sub.n) as a function of the distance in k-space from point m to point n. The covariance may include information from the coil sensitivities, the object being scanned, or any other source of covariance.
[0031] As one example, the covariances can be approximated from fully sampled images, f.sub.a and f.sub.b, from coils c.sub.a and c.sub.b, respectively. By way of example, the images, f.sub.a and f.sub.b, can be low resolution images reconstructed from a central portion of k-space. To estimate the covariances, the image f.sub.a can be multiplied by the conjugate image f.sub.b* and the result Fourier transformed,
[0032] This forms a matrix of covariance maps for each coil pair, which also includes object-induced covariance. Convolution in k-space by these maps and then resampling at locations corresponding to x or y is equivalent to multiplying a vector of data by R.sub.xx or R.sub.xy, respectively. These covariance maps may be windowed before convolution.
[0033] Because of the inverse matrix in Eqn. (9), it is generally impractical to estimate the synthesized data, ŷ, directly, even with exact knowledge of R.sub.xx and R.sub.xy. To overcome this limitation, an intermediate step can be performed to relax the computational burden of estimating the synthesized data. An intermediate parameter, z, can be defined via the following:
x=R.sub.xxz (11).
[0034] Using this parameter, the synthesized data can be estimated using a two-step process. First, Eqn. (11) can be solved using a linear algebra technique, such as the conjugate gradient method. In this example, to force the operating matrix to be positive semi-definite and to ensure convergence, Eqn. (11) can be rewritten as,
R.sub.xx.sup.Hx=R.sub.xx.sup.HR.sub.xxz (12);
[0035] where multiplication by R.sub.xx.sup.H is achieved by flipping the convolution kernel used for R.sub.xx for each individual coil pair, but not taking its conjugate. That is, r.sub.ab[−dk] can be used. The values of the vector, z, correspond to the same corresponding k-space locations as those for the sampled data vector, x. This is used for the dk index of the convolution kernels, r.sub.ab.
[0036] With an estimate of z, Eqn. (9) now becomes,
ŷ=R.sub.yxz (13);
[0037] which is a simple convolution process from z (located at the sampled data locations) to the desired locations for y. The values for y can thus be defined on arbitrary trajectories, including possibly Cartesian k-space locations, which would obviate the need for regridding the data after this process.
[0038] Solving Eqn. (12) can be the most time-consuming part of the above-described methods. Assuming their are N coils for both y and x, their are N.sup.2 convolutions that must be performed and N.sup.2 convolution kernels that must be stored in memory. If their are also M data sets (e.g. for diffusion imaging or temporal imaging) over which correlations will be estimated, then the number of stored kernels and required convolutions becomes (M*N).sup.2, which can become computationally burdensome, particularly for 3D imaging. Note that the multiplication by the matrix R.sub.xx is implemented by convolution with a kernel rab, given in Eqn. (10). That is, for a given iteration between data sets (e.g. coils) a and b:
R.sub.xxzz
r.sub.ab (14).
[0039] From Eqn. (10), r.sub.ab can be substituted with the Fourier transforms of f.sub.a(x) and f*.sub.b(x). That is,
R.sub.xxzz
[F*.sub.b(−k)
F.sub.a(k)] (15).
[0040] The vector, z, can then be broken into the constituent set b, and Eqn. (15) can be rearranged as:
[0041] The convolution in square brackets of Eqn. (16) is performed for all M*N sets, and can be stored (i.e. convolved onto) a sufficiently dense Cartesian grid. The second M*N convolutions (with F.sub.a) can then be performed off of this single Cartesian grid, for a complete process that requires only 2*M*N convolutions and requires storage of M*N kernels (instead of the original [M*N].sup.2 kernels and convolutions required).
[0042] Referring now to
[0043] An example of a variable density spiral trajectory includes radial data undersampling that begins at a k-space radius equal to ten percent of the maximum radius; increases to some value, K, at thirty percent of the maximum radius; and then remains constant at K for the remainder of k-space sampling.
[0044] After the data are acquired, low-resolution images are reconstructed from the acquired data, as indicated at step 104. For instance, low-resolution images are reconstructed from the central portion of k-space, which in the variable density spiral acquisition corresponds to a fully sampled region of k-space.
[0045] Covariance maps are then estimated from the low-resolution images, as indicated at step 106. In one example, the covariance maps are estimated based on pairwise multiplication of the low-resolution images corresponding to pairs of receive coils, as described above. For instance, the low-resolution image, f.sub.a, associated with a first coil, c.sub.a, is multiplied with the conjugate of the low-resolution image, f.sub.b*, associated with a second coil, c.sub.b. The resulting product, f.sub.af.sub.b*, is then Fourier transformed to obtain the estimate of the covariance maps, r.sub.ab[dk].
[0046] In some embodiments, the covariance maps are cropped to a desired convolution radius and then windowed, such as by using a Hanning window. This windowing of r.sub.ab appears to advantageously stabilize the solution to Eqn. (10). To understand this stabilization, multiplying by R.sub.xx.sup.−1 can be viewed as equivalent to deconvolving by r.sub.ab in k-space, which is in turn the same as dividing by its Fourier transform, R.sub.ab, in image space. To keep R.sub.ab from have zero-crossings in the image, r.sub.ab can be windowed as described above.
[0047] A desired sampling pattern for the synthesized data, y, is provided, selected, or otherwise computed, as indicated at step 108. As mentioned above, the sampling pattern for the synthesized data can be any arbitrary sampling pattern, whether Cartesian or non-Cartesian. As one example, from the region where undersampling starts in a variable density spiral acquisition (e.g., ten percent maximum radius) to the edge, coordinates for y can be generated as spiral trajectories matching those used for data measurement, but rotated to be equally spaced between the measured arms of x. The number of y arms between each pair of x arms is equal to the ceiling of K−1, where the ceiling function is used should K be a non-integer value, so that the radial spacing of the data including both y and x meets or exceeds the Nyquist sampling constraint of 1/FOV.
[0048] In some embodiments, the synthesized data can be defined on a Cartesian grid, obviating the need for subsequent gridding prior to image formation. The synthesized data can also be defined on a single, potentially homogenous, receive coil. This latter embodiment would obviate the need for subsequent coil combination and reducing the gridding operation to a single data set prior to image formation.
[0049] The synthesized data are estimated next, as indicated at step 110. As described above, the synthesized data can be estimated in a two-step process. First, the covariances computed between the measured data are utilized in Eqn. (11) or (12) to estimate an intermediate vector, z. Eqn. (13) can then be used to convolve this intermediate vector, z, to the provided sampling pattern for the synthesized data, y, which thus estimates the synthesized data. As mentioned above, in some embodiments, the measured and synthesized data can be on the same or different coils, which in turn can be real or virtual.
[0050] Based at least on the synthesized data, y, an image of the subject is then reconstructed, as indicated at step 112. In some embodiments, the image can be reconstructed from the synthesized data, y, alone. In some other embodiments, however, the image can be reconstructed from a combination of the measured data, x, and the synthesized data, y.
[0051] In some embodiments, the systems and methods described here can be used in any dimensionality, including space, time, and frequency, where covariances can be estimated.
[0052] The systems and methods described here can be incorporated with other encoding that affects covariance. For example, the phase imposed from motion during the application of diffusion weighting gradients can be incorporated into the algorithm described here.
[0053] Referring particularly now to
[0054] The pulse sequence server 210 functions in response to instructions downloaded from the operator workstation 202 to operate a gradient system 218 and a radiofrequency (“RF”) system 220. Gradient waveforms necessary to perform the prescribed scan are produced and applied to the gradient system 218, which excites gradient coils in an assembly 222 to produce the magnetic field gradients G.sub.x, G.sub.y, and G.sub.z used for position encoding magnetic resonance signals. The gradient coil assembly 222 forms part of a magnet assembly 224 that includes a polarizing magnet 226 and a whole-body RF coil 228.
[0055] RF waveforms are applied by the RF system 220 to the RF coil 228, or a separate local coil (not shown in
[0056] The RF system 220 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 228 to which it is connected, and a detector that detects and digitizes the I 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 I and Q components:
M=√{square root over (I.sup.2+Q.sup.2)} (14);
[0057] and the phase of the received magnetic resonance signal may also be determined according to the following relationship:
[0058] The pulse sequence server 210 also optionally receives patient data from a physiological acquisition controller 230. By way of example, the physiological acquisition controller 230 may receive signals from a number of different sensors connected to the patient, such as electrocardiograph (“ECG”) signals from electrodes, or respiratory signals from a respiratory bellows or other respiratory monitoring device. Such signals are typically used by the pulse sequence server 210 to synchronize, or “gate,” the performance of the scan with the subject's heart beat or respiration.
[0059] The pulse sequence server 210 also connects to a scan room interface circuit 232 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 232 that a patient positioning system 234 receives commands to move the patient to desired positions during the scan.
[0060] The digitized magnetic resonance signal samples produced by the RF system 220 are received by the data acquisition server 212. The data acquisition server 212 operates in response to instructions downloaded from the operator workstation 202 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 212 does little more than pass the acquired magnetic resonance data to the data processor server 214. However, in scans that require information derived from acquired magnetic resonance data to control the further performance of the scan, the data acquisition server 212 is programmed to produce such information and convey it to the pulse sequence server 210. For example, during prescans, magnetic resonance data is acquired and used to calibrate the pulse sequence performed by the pulse sequence server 210. As another example, navigator signals may be acquired and used to adjust the operating parameters of the RF system 220 or the gradient system 218, or to control the view order in which k-space is sampled. In still another example, the data acquisition server 212 may also be employed to process magnetic resonance signals used to detect the arrival of a contrast agent in a magnetic resonance angiography (“MRA”) scan. By way of example, the data acquisition server 212 acquires magnetic resonance data and processes it in real-time to produce information that is used to control the scan.
[0061] The data processing server 214 receives magnetic resonance data from the data acquisition server 212 and processes it in accordance with instructions downloaded from the operator workstation 202. 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.
[0062] Images reconstructed by the data processing server 214 are conveyed back to the operator workstation 202 where they are stored. Real-time images are stored in a data base memory cache (not shown in
[0063] The MRI system 200 may also include one or more networked workstations 242. By way of example, a networked workstation 242 may include a display 244; one or more input devices 246, such as a keyboard and mouse; and a processor 248. The networked workstation 242 may be located within the same facility as the operator workstation 202, or in a different facility, such as a different healthcare institution or clinic.
[0064] The networked workstation 242, whether within the same facility or in a different facility as the operator workstation 202, may gain remote access to the data processing server 214 or data store server 216 via the communication system 240. Accordingly, multiple networked workstations 242 may have access to the data processing server 214 and the data store server 216. In this manner, magnetic resonance data, reconstructed images, or other data may be exchanged between the data processing server 214 or the data store server 216 and the networked workstations 242, such that the data or images may be remotely processed by a networked workstation 242. This data may be 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.
[0065] 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.