WO2009037711A2

WO2009037711A2 - Mri image reconstruction method involving nonlinear encoding gradients or phase cycling

Info

Publication number: WO2009037711A2
Application number: PCT/IL2008/001261
Authority: WO
Inventors: Yuval Zur
Original assignee: Topspin Medical (Israel) Ltd.
Priority date: 2007-09-20
Filing date: 2008-09-21
Publication date: 2009-03-26
Also published as: WO2009037710A2; WO2009037709A3; WO2009037709A2; WO2009037711A3; WO2009037710A3

Abstract

A method of reconstructing an image from a set of MRI data points in space acquired from an imaging slice, using space encoding with one or more gradient fields that are each nonlinear as functions of position in a space encoding direction, the method comprising: a) multiplying a vector of the data points by a reconstruction matrix which differs from a fourier transform matrix at least due to taking into account the nonlinearity of the one or more gradient fields; and b) obtaining a vector of voxel densities for the image.

Description

IMAGE RECONSTRUCTION METHODS FOR MRI PROBES

RELATED APPICATIONS

The present application claims benefit under 35 USC 119(e) from US provisional patent applications 60/960,212, and 60/960,213, both filed on September 20, 2007. This application is related to two PCT applications filed on even date, the first titled "MRI Probe," having attorney docket number 44842, and the second titled "MRI Magnet and Coil Configurations," having attorney docket number 44841. The disclosures of all of these applications are incorporated herein by reference. This application is also related to the following applications: U.S. Patent

Application No. 10/968,853 filed October 18, 2004, entitled "Magnet and Coil Configurations for MRI probes" and published as US2006/0084861; U.S. Patent Application No. 10/597,325 filed July 20, 2006 and entitled "MRI Probe for Prostate Imaging"; PCT Patent Application No. PCT/IL2005/000074 filed January 20, 2005, entitled "MRI Probe for Prostate Imaging" and published as WO 2005/067392; and U.S. Provisional Patent Application No. 60/537,030 filed January 20, 2004 and entitled "MRI Probe for Prostate Imaging". The disclosures of the above applications are each fully incoφorated herein by reference.

FIELD AND BACKGROUND OF THE INVENTION

The present invention, in some embodiments thereof, relates to self-contained MRI probes and methods of image reconstruction for them, and, more particularly, but not exclusively, to methods of image reconstruction for probes which have substantially nonuniform field gradients, and/or have limited signal to noise ratio (SNR), which is often the case for MRI probes small enough to use inside the body..

Conventional medical MRI systems have a field of view with very uniform static magnetic field in the bore of a large magnet, and use gradient coils that produce highly linear gradient fields in three dimensions, for gradient encoding of images. Conventional MRI systems sometimes use RF receiver probes that can be used inside blood vessels, the rectum (for prostate imaging) and other body cavities, in order to improve signal to noise ratio (SNR) when imaging regions adjacent to the probe, as described, for example, by US patent 5,699,801 to Atalar, and in US patent 5,476,095 to Schnall et al.

In order to reduce the cost of MRI systems, while allowing imaging of small regions inside the body with high resolution and SNR, fully self-contained MRI probes have been proposed. In these "inside out" MRI devices, inspired by NMR probes used for well- logging, the field of view is outside the probe. US 5,572,132 to Pulyer et al describes a self- contained MRI probe, for use in blood vessels or other body cavities, with permanent magnets configured to produce a static magnetic field with a saddle point, outside the probe. The imaging region of the probe is a small region around the saddle point, where the magnetic field is locally very uniform. Gradient coils produce field gradients in the imaging region, in the radial, azimuthal, and longitudinal directions, for gradient encoding of images.

US patent 5,304,930 to Crowley et al describes a larger self-contained "inside-out" MRI device, designed to be used outside the body, to image a part of the body. In this device, the static magnetic field is not nearly uniform, but has a substantial gradient in one direction, for example the x-direction, which is used for gradient encoding of images. This field gradient is required to be very linear, and the surfaces of constant field are required to be very flat, in the y-z plane.

US patent 6,489,767 to Prado et al.describes an "NMR-MOUSE" type of self- cotained MRI probe, too large to fit inside the body, which uses gradient coils for phase encoding.

K. P. Pruessmann et al, "SENSE: Sensitivity Encoding for Fast MRI," Magnetic Resonance in Medicine 42:952-962 (1999), describes a method of using multiple RF receiving antennas in conventional MRI in a very uniform magnetic field and very linear gradient fields. The different receiving antennas have different spatial patterns of sensitivity, which can be used to decrease acquisition time when used in conjunction with gradient phase encoding. Pruessmann et al describe how to find a reconstruction matrix which is optimized for a given set of receiving antennas, in that it tranforms MRI data in phase space into voxel densities, for voxels which are as close as possible to a set of ideal voxel functions, but can be expressed as a linear combination of phase encoding functions, including the spatial sensitivity of the receiving coils.

Y. Zur (the inventor on the present application) and S. Stokar, "A Phase-Cycling

Technique for Canceling Spurious Echoes in NMR Imaging," J. of Magnetic Resonance 71, 212-228 (1987) describes a technique, in conventional MRI, for reducing the effects of image artifacts caused by imperfect RF pulses which produce spurious echoes. The technique uses phase cycling to separate echoes of the correct and the wrong parity.

The contents of all of the above documents are incorporated by reference as if fully set forth herein.

SUMMARY OF THE INVENTION

An aspect of some embodiments of the invention concerns methods of MRI image reconstruction which provide improved accuracy when the gradient fields are mildly or very nonlinear functions of position, and/or provide improved SNR when SNR is marginal.

These are typical characteristics of self-contained MRI probes, especially those small enough to use inside the body, so the techniques may be especially useful for such probes. There is thus provided, in accordance with an exemplary embodiment of the invention, a method of reconstructing an image from a set of MRI data points in phase space acquired from an imaging slice, using phase encoding with one or more gradient fields that are each nonlinear as functions of position in a phase encoding direction, the method comprising: a) multiplying a vector of the data points by a reconstruction matrix which differs from a fourier transform matrix at least due to taking into account the nonlinearity of the one or more gradient fields; and b) obtaining a vector of voxel densities for the image.

Optionally, the reconstruction matrix is the product of an encoding matrix and the inverse of a correlation matrix of phase encoding functions for the phase encoding with the one or more nonlinear gradient fields, each element of the phase encoding matrix being an integral, over at least one phase encoding direction, of a product of one of the phase encoding functions and one of a set of ideal voxel functions, each of which is compact around a different grid point in at least one of the one or more phase encoding directions. Optionally, the ideal voxel functions are delta functions, boxcar functions, half-sin functions, or cubic spline functions of position in at least one of the one or more phase encoding directions.

Optionally, the MRI data points consist of K points with phase encoding differing only in characteristics of a single gradient field pulse which varies in one phase encoding direction, the voxel densities consist of densities of N voxel functions of the one phase encoding direction, the reconstruction matrix is a K by N matrix, and the vector of the voxel densities is the product of the vector of the K data points, and the reconstruction matrix for the one phase encoding direction.

In an embodiment of the invention, the imaging slice has two phase encoding directions, and the MRI data points are phase encoded by two gradient fields, one for each of the directions. Optionally, the MRI data points comprise an array of Kx by Ky points in two dimensional phase space, the voxel densities comprise an array of Nx by Ny densities for an array of two-dimensional voxel functions each of which is a product of one of a set of Nx voxel functions of x, a first of the two phase encoding directions, and one of a set of Ny voxel functions of y, a second of the two phase encoding directions, and transforming the data points comprises: a) multiplying a matrix of the Kx by Ky data points by a Kx by Nx reconstruction matrix for the first of the two phase encoding directions; and b) multiplying the product by a Ky by Ny reconstruction matrix for the second of the two phase encoding directions, to obtain an Nx by Ny matrix of the voxel densities. Optionally, the one-dimensional voxel functions of x have a minimum of a measure of difference from a set of Nx ideal compact voxel functions localized around an array of grid points, for any set of Nx one-dimensional functions that can be expressed as linear combinations of Kx phase encoding functions of x for Kx evenly spaced amplitudes of the gradient field for x, the reconstruction matrix for x being a linear transformation from the phase encoding functions to the voxel functions of x, and wherein the one-dimensional voxel functions of y have a minimum of a measure of difference from a set of Ny ideal compact voxel functions localized around an array of grid points, for any set of Ny one- dimensional functions that can expressed as linear combinations of Ky phase encoding functions of y for Ky evenly spaced amplitudes of the gradient field for y, the reconstruction matrix for y being a linear transformation from the phase encoding functions to the voxel functions of y.

Alternatively, the MRI data points comprise K data points, each characterized by a different pair (kx,ky), the voxel functions comprise N voxel functions, each peaked at a different pair of values (x,y) of the two dimensions x and y, and transforming the data points comprises multiplying a vector of the K data points by a K by N reconstruction matrix for the two dimensions, to obtain a vector of the N voxel densities in the two dimensions.

There is further provided, in accordance with an exemplary embodiment of the invention, a method of calculating an image reconstruction matrix for MRI data acquired from an imaging slice using phase encoding with one or more gradient fields that are each a nonlinear function of position in a phase encoding direction, the method comprising: a) choosing a set of ideal voxel functions; b) calculating a set of phase encoding functions, one for each value of each of the one or more gradient fields at which the MRI data are acquired; and c) finding a set of voxel functions, each of which has a minimum of a measure of difference from one of the ideal voxel functions, subject to the constrain that the voxel function can be expressed as a linear combination of the phase encoding functions; wherein the reconstruction matrix is the linear transformation from the phase encoding functions to the voxel functions.

Optionally, finding the set of voxel functions comprises: a) calculating an encoding matrix, each element of which is an integral over the one or more phase encoding directions of one of the phase encoding functions and one of the ideal voxel functions; b) calculating a correlation matrix of the phase encoding functions; c) multiplying the encoding matrix by the inverse of the correlation matrix, thereby obtaining the reconstruction matrix which is a linear transformation from the phase encoding functions to the voxel functions.

There is further provided, in accordance with an exemplary embodiment of the in invention, a method of reconstructing an image in an imaging slice from a set of MRI data points in phase space acquired using phase encoding by one or more gradient fields, at least one of which has a field gradient that is non-uniform in magnitude or direction or both, in the imaging slice, the method comprising: a) transforming the set of data points to a set of voxel densities, each voxel being identified by a different set of values of the one of more gradient fields; b) for each voxel, determining a position in the imaging slice where the one or more gradient fields have the values in the set identifying that voxel; and c) reconstructing the image by calculating a density as a function of position in the imaging slice, using the voxel densities and the positions in the imaging slice of the voxels.

Optionally, calculating a density as a function of position in the imaging slice comprises adjusting one or more of the voxel densities responsive to an estimated volume of the voxel calculated from a gradient of each of the one or more gradient fields at the position of the voxel. There is further provided, in accordance with an exemplary embodiment of the invention, a method of generating an MRI image of an imaging region with reduced noise, comprising: a) applying to one or more slices of the imaging region NMR pulse sequences suitable for separating even parity and odd parity echoes, for a plurality of points in phase space; b) accumulating NMR signals for each point in phase space for each slice, and separating the accumulated signals into a set of even parity data points and a set of odd parity data points in phase space; c) identifying an expected phase of a signal part of each data point; d) subtracting from each data point a part that is expected to be substantially noise because it has phase orthogonal to the expected phase of the signal part; and e) using the sets of even parity and odd parity data points after subtracting the parts expected to be noise, to reconstruct an image with reduced noise. Optionally, the NMR pulse sequences comprise, for a point in phase space for a slice, a first NMR pulse sequence consisting of an excitation RF pulse, followed by one or more gradient pulses, followed by a train of refocusing pulses with substantially zero phase difference from the excitation pulse, and a second NMR pulse sequence consisting of an excitation pulse, followed by one or more gradient pulses, followed by a train of refocusing pulses with substantially 90 degree phase difference from the excitation pulse, and separating the signals into sets of even and odd parity data points comprises finding a sum of the accumulated signals from the first and second pulse sequences for each data point, and finding a difference of the accumulated signals from the first and second pulse sequences for each data point. Optionally, using the sets of even parity and odd parity data points to reconstruct an image comprises adding data of one parity to data of the other parity that has been flipped in each phase encoding direction, either in phase space or in position space.

Optionally, identifying an expected phase of the signal part of a data point comprises finding phases of one or more data points for a plurality of different resonance frequencies, fitting the phases to a function of resonance frequency, and identifying the value of the fitted function of resonance frequency as the expected phase for all data points in phase space for that resonance frequency.

There is further provided, according to an exemplary embodiment of the invention, a method of MRI imaging, comprising: a) applying an NMR pulse sequence to one or more imaging slices of an imaging region; b) acquiring a set of MRI data points from the one or more imaging slices; and c) reconstructing an image of the imaging region using any of the methods according to an exemplary embodiment of the invention. There is further provided, according to an exemplary embodiment of the invention, an MRI imaging system, comprising: a) a self-contained MRI probe comprising a magnet assembly that produces a static magnetic field in an imaging region, one or more sets of gradient coils that each produce a magnetic gradient field in the imaging region, and one or more RF antennas that together produce RF fields that excite nuclei in the imaging region and receive NMR signals from the excited nuclei; and b) a controller that directs the probe to produce pulse sequences of gradient fields and RF fields in the one or more slices of the imaging region, and that accumulates NMR signals to produce MRI data points in phase space, and that reconstructs an image from the MRI data points according to any of the methods according to an exemplary embodiment of the invention.

BRIEF DESCRIPTION OF THE DRAWINGS Some embodiments of the invention are herein described, by way of example only, with reference to the accompanying drawings. With specific reference now to the drawings in detail, it is stressed that the particulars shown are by way of example and for purposes of illustrative discussion of embodiments of the invention. In this regard, the description taken with the drawings makes apparent to those skilled in the art how embodiments of the invention may be practiced. In the drawings:

FIG. 1 is a flowchart showing a phase-cycling method of reconstructing an MRI image with reduced noise, according to an exemplary embodiment of the invention;

FIG. 2 schematically shows an MRI pulse sequence that could be used for the phase- cycling method of FIG. 1 ;

FIG. 3 is flowchart for a method of finding a reconstruction matrix that can be used for MRI data acquired using phase encoding in a nonlinear gradient field, according to an exemplary embodiment of the invention;

FIG. 4 schematically shows a one-dimensional density profile and a nonlinear gradient field profile, to illustrate the method of FIG. 3; FIG's. 5A and 5B schematically show voxel functions found using the method of FIG. 3 with the gradient field profile of FIG. 4;

FIG. 6 schematically shows a plot of simulated MRI data points for the density and gradient field profiles of FIG. 4; FIG. 7 schematically shows a reconstructed density profile, found for the density profile and gradient field profile of FIG. 4, using a reconstruction matrix found using the method of FIG. 3, according to an exemplary embodiment of the invention;

FIG. 8 schematically shows a two-dimensional object, to illustrate reconstructing two-dimensional MRI images using a reconstruction matrix found by the method of FIG. 3, according to an exemplary embodiment of the invention;

FIG. 9 schematically shows density profiles reconstructed in two different ways from simulated MRI data of the object of FIG. 8, according to an exemplary embodiment of the invention;

FIG. 10 schematically shows contours of constant gradient field, for each of two gradient fields that are nonlinear functions of two gradient encoding directions, to illustrate an image reconstruction method according to an exemplary embodiment of the invention; and

FIG. 11 schematically shows an MRI system which uses methods for image reconstruction according to an exemplary embodiment of the invention, including one or more of the methods of FIG's. 1, 3, and 10.

DESCRIPTION OF EMBODIMENTS OF THE INVENTION

An aspect of some embodiments of the invention concerns a method of using phase cycling to improve SNR in MRI, by separating MRI data points of odd and even parity, correcting their phase, removing a part which is likely to be substantially all noise because of its phase, and combining the even and odd parity data to form an image. The even and odd parity data is separated by running two pulse sequences for each point in phase space and each slice, for example: a first sequence consisting of an excitation pulse followed by one or more gradient pulses and a train of refocusing pulses that are 90 degrees out of phase with the excitation pulse, and a similar second sequence in which the refocusing pulses are in phase with the excitation pulse. By accumulating echoes for the first and second sequences, and taking their sum and their difference, one can obtain the even parity and odd parity data for that point in phase space. Optionally, the phase is corrected by measuring the phase of the zero-gradient or low k components of the data (in effect the k = 0 data point or a few low k data points) as a function of resonance frequency, over a range of resonance frequencies where there is a significant signal. A best fit, for example a linear fit with an offset, is found for the phase as a function of resonance frequency. The phases for all data points, including those further out in phase space, are then reduced by this fitted function. After this phase correction, the dominant points in phase space will generally all be nearly real. The remaining imaginary parts, after phase correction, are likely to be substantially all noise, and can be removed, and the data points in phase space are used to reconstruct an image in actual space. The phase-corrected odd and even parity data is combined by first flipping one of them with respect to each of the gradient encoding axes, either in phase space or in actual space, and then adding the even and odd parity data points or images together.

An aspect of some embodiments of the invention concerns a method of image reconstruction for MRI data obtained using phase encoding in one or more phase encoding directions, with a gradient field that is a nonlinear function of the position in at least one gradient encoding direction. The method uses a reconstruction matrix which takes into account the nonlinearities of the gradient fields. The method is optionally more computationally intensive than the usual FFT methods that are used with highly linear gradients in conventional MRI, but provides more accurate images when one or more gradient field is a substantially nonlinear function of position in a phase encoding direction. Optionally, the reconstruction matrix is chosen so that the associated voxel functions will match as closely as possible ideal voxel functions, for example compact functions around a set of grid points, such as delta functions or boxcar functions, subject to the constraint that the voxel functions are linear combinations of phase encoding functions associated with the nonlinear field gradients. Optionally, the method is used in one phase encoding direction and one gradient field at a time, either to reconstruct a one-dimensional image, or to reconstruct a two-dimensional image with gradient fields and voxel functions that can be treated to good approximation as separable in the two gradient encoding directions. Alternatively, the method is used in a fully two-dimensional way, to reconstruct a two- dimensional image.

An aspect of some embodiments of the invention concerns a method of calculating an image reconstruction matrix, for use with MRI data acquired using phase encoding with one or more gradient fields that are nonlinear functions of position in the phase encoding direction.

An aspect of some embodiments of the invention concerns a method of image reconstruction, for MRI data acquired using phase encoding with one or more gradient field that is a nonlinear function of position in the phase encoding direction, taking into account changes in the positions that different voxels would have if the gradient fields were linear functions of position. The method can be used if a reconstruction matrix is used to reconstruct the image, as described above, or it can be used to correct distortions in images reconstructed by conventional methods such as FFT which do not take into account the nonlinearity of the gradient fields.

Phase-cycling method for reducing noise

FIG. 1 is a flowchart showing a method of using phase-cycling to reduce noise in MRI images, according to an exemplary embodiment of the invention. The method relies on the use of pulse sequences that can be used to separate MRI data into echoes of even parity, and odd parity. Several examples of such pulse sequences, and how they are used, are described by Zur and Stokar in the paper cited above. One such example will be described here, with reference to the pulse sequences shown in FIG. 2.

At 102, a first pulse sequence is applied to the imaging region, and data is accumulated. A second pulse sequence is applied, and data is accumulated, in 104. As an example, the first pulse sequence is a CPMG type of pulse sequence. An excitation pulse 902 at frequency fi, selecting a particular imaging slice, nominally with 90 degree flip angle, is applied first. This is followed by a gradient pulse 904, followed by a refocusing pulse 906 at frequency fi at a time τ following the excitation pulse, at a phase that differs by 90 degrees from the phase of the excitation pulse. At a further time τ following refocusing pulse 906, an echo 908 occurs, and is accumulated as part of the MRI data for this value of the integrated strength of gradient pulse 906. A train of further refocusing pulses 910 is applied at frequency fi, at time intervals of 2τ, and echoes 912 are received and accumulated half-way between successive refocusing pulses. Train 910 continues for a time T, which may be adjusted to produce a desired degree of T₂ or diffusion weighting. Alternatively, T₂ or diffusion time Tj is found directly by measuring the decay rate of echoes in train 910, so the weighting does not depend on T. Train 910 is followed by another series of pulses, at a different frequency f₂ to excite a different imaging slice, including an excitation pulse 914, a gradient pulse 916, and refocusing pulses 918, 90 degrees out of phase with the excitation pulse, with echoes 920 halfway between the refocusing pulses. After a time T_R, which may be chosen to adjust the degree of Ti weighting, a second pulse sequence 922 is applied at frequency fi, with the same integrated strength of the gradient pulse, and the echoes are accumulated into an MRI data point, but this time the refocusing pulses are in phase with the excitation pulse. Again, a series of pulses, each extending for a time T, at different frequencies, is followed after a time T_R by another first type of pulse sequence at frequency fi, this time with a different integrated strength of the gradient pulse, representing a different point in phase space, generally called k-space in the case of linear field gradients. Data from the different imaging slices, from the different points in phase space, and from pulse sequences of the first or second type, is accumulated separately.

It should be noted that many variations on the foregoing description of pulse sequences are possible. For example, the time interval between the excitation pulse and the first refocusing pulse, which is also the time interval between the first refocusing pulse and the first echo, need not be the same interval τ as the interval between the other refocusing pulses in the train, but can be several times longer, as explained in the co-filed application "MRI Magnet and Coil Configurations." This may allow more time for the gradient pulse, which can result in less heating of the gradient coil for the same integrated gradient field strength. A second gradient pulse may be applied during the same longer time interval between the first refocusing pulse and the beginning of the train of refocusing pulses. Pulse sequences of the first and second type need not be applied consecutively for the same point in phase space, but, for example, pulses of the first type maybe applied at all points in phase space, before moving on to pulses sequences of the second type., or any other pattern of interleaving may be used.

Returning to FIG. 1, at 106, data from the first and second types of pulse sequences is used to obtain odd parity data and even parity data. Odd parity data is data from echoes emitted by components of the excited nuclei that have undergone an odd number of reversals in phase due to refocusing pulses, while even parity data is data from echoes emitted by components of the excited nuclei that have undergone an even number of reversals in phase sue to refocusing pulses. Although the primary contribution to echoes, at least early on in a train of refocusing pulses, will alternate between odd and even parity echoes, there will also be contributions with the other parity, due to refocusing pulses that are not exactly 180 degrees, which produce components of excited nuclei polarized in the z-direction, which are not subject to refocusing until they go back to the x-y plane.

The odd parity and even parity components of the data are optionally obtained by first adding together the data from the first and second types of pulse sequences, for a given point in phase space and a given slice, and then subtracting them.

At 108, optionally, the odd and even data is optionally separately corrected for phase errors. Such phase errors are likely to be primarily due to a shift in the timing of the echo, due to errors in the timing of the pulses, and these phase errors will be linear with resonance frequency. In addition there may be small instrumental phase errors independent of resonance frequency. In order to correct phase errors, a range of frequencies is optionally selected at which there is a substantial signal. At the dominant point in phase space, normally k = 0, the phase of the data point is measured as a function of frequency, and fit to a function of frequency, generally a linear fit with a constant offset. Alternatively, instead of only using the k = 0 component, an average is taken of the phase of a few low k components which have large signal, optionally a weighted average. Since the signal from the k = 0 point, or low k points, is expected to be real, the fitted function for the measured phase of the k = 0 points or low k points is optionally subtracted from the measured phase of all points in phase space, to yield a corrected phase for all the points in phase space., for all resonance frequencies in the imaging region. It should be noted that this procedure may not work well if there is not at least one imaging slice, and preferably a few imaging slices, for which at least one data point, typically k = 0, has a high enough SNR for its phase to be measured reasonably accurately, for example to within much better than 90 degrees. But if the phase can be measured reasonably accurately for a few such imaging slices and data points, or even for one of them, then that information can be used to reduce noise for other imaging slices and data points in phase space as well, even if they have lower SNR.

At 110, the real part of the data is taken, in order to reduce noise in the data, since the imaginary parts of phase corrected data points are likely to be only due only to noise.

This generally reduces the noise by a factor of V2^~ , and increases the SNR by the same factor.

At 112, the even and odd parity data is used to reconstruct a combined image incorporating all the data. However, one parity of data is reflected, in all the phase encoding directions, before adding it to the other parity of data, because odd echoes are reversed in phase relative to even echoes. This is done, for example, by reconstructing separate images from the even parity and odd parity data, reflecting one of the images along each of the phase encoding axes, and then adding the two images together. Alternatively, before reconstructing images, data of one of the parities is reflected in phase space along each of the phase encoding directions, and the data points of the two parities are then added together. The resulting data points are then used to reconstruct an image. The image reconstruction can be done using any known method of reconstructing an image from phase encoded MRI data. For example, in the case of an MRI system with reasonably linear gradient fields, a method such as FFT may be used. In the case of a small, self-contained MRI probe, or any other MRI system where the gradient fields are significantly nonlinear, other methods, such as those described below, may yield more accurate images. Image Reconstruction Matrix with Nonlinear Gradient Fields

FIG. 3 is a flowchart 300, describing a method for finding a reconstruction matrix which can be used to convert a vector of phase encoded MRI data points into a vector of voxel densities, taking into account nonlinearities in the gradient fields that are used to phase encode the data. The image reconstruction algorithm may be more computationally intensive then FFT, which is used when the gradient fields are highly linear, because it is much faster than multiplying an array of data points in phase space by an arbitrary reconstruction matrix, to find an array of voxel densities. However, the reconstruction matrix may yield more accurate results than FFT, and other conventional methods of image reconstruction, especially with high nonlinear gradient fields, when FFT may fail completely. Calculating the reconstruction matrix may be computationally intensive, but only has to be done once for a given set of gradient field profiles.

The method will be described first for a one-dimensional case, with only a single phase encoding direction y, and voxel function depending only y. However, the method may be used for more than one phase encoding direction, and two dimensional voxel functions, simply by replacing the single indexes for points in phase space, and grid points in configuration space, by pairs of indexes, representing two phase encoding directions, and replacing integrals over the single phase encoding direction y by integrals over both phase encoding directions. Another way of using the reconstruction matrix for two-dimensional imaging, which is much less computational intensive and which may be used when the two directions are separable to good approximation, will be described afterwards.

A one-dimensional image reconstruction may be performed if only a one- dimensional slice is being excited, or if a two-dimensional slice is being excited, but we are only interested in finding the z dependence of density for a particular ky component, for example for ky = 0, which means an average over y. This might be of interest in particular if the spin density were at least approximately separable in y and z.

At 304, phase encoding functions e_k(y) are found, for each of K points in phase space, labeled by an index k. These encoding functions are defined as exp[iφk(y)], where '_k (r) = χ ' k' \ B(y) \ 'T_{lβ (1)} is the phase at point y caused by the generally nonlinear gradient field B(y), when applied with a strength multiplied by k, for the k^m point in phase space, for a time interval T_h5. If the gradient field is superimposed on a much stronger static magnetic field, as is usually the case in MRI, then only the component of the gradient field parallel to the static magnetic field contributes to this phase. For a linear gradient field, B(y) = Gy, the phase φk(y) ⁼ γ-kGT_hs-y, and the phase encoding function e^y) = exp[iγ-kGTh_S-y]. The phase encoding functions are used to weight the spin density p(y) to produce the data points m^ for each k.

It should be understood that the spin density p(y) refers not the total density of hydrogen nuclei (or whatever species is being imaged), but the density weighted by Ti, T₂, or diffusion, or any other weighting, whatever the case may be. The weighting depends on parameters of the pulse sequence used to obtain the data points, or on analysis of the echoes that are received at different times during the pulse sequence, as is known in the art MRI. However, the reconstruction methods to be described here start with the data points mk, regardless of how they are obtained, so can be used with any type of weighting.

At 302, ideal voxel functions I_n(y) are chosen for each of N grid points indexed by n. Typically, the ideal voxel functions may be compact functions centered around grid points y_n, for example delta functions δ(y - y_n), or boxcar functions of y - y_n, or half-sin functions, or cubic splines. A goal is to find actual voxel functions f_n (y) which are as close as possible to the ideal voxel functions, given the constraint that they are linear combinations of the phase encoding functions Q^{y). The actual voxel functions depend on an N by K reconstruction matrix F_n^ which relates the data points m _k to voxel densities V _n of the reconstructed image:

rn = ∑ Fn,k ^m k (3) which may be written in vector form at V = F^#m. Putting Eq. (2) for m into Eq. (3) gives

and the quantity in parentheses is defined at the n^th voxel function for the reconstruction matrix F,

which may be expressed in matrix form at f = F»e To find a matrix F which has voxel functions that approximate the ideal voxel functions, a K by N encoding matrix Ek,_n is optionally found at 306, defined by

Ek,n = (y) ^• e_k (^r) ^{• d}y (6) which may be expressed in matrix form as E = I^H»e , where I^H is the conjugate transpose of I.

At 308, a K by K correlation matrix C of the encoding functions is optionally found, defined by

or C = e^H#e*dy in matrix form.

The reconstruction matrix F is found at 310, optionally from

F = E^{H •} C^'1 (8)

It can be shown that with this form for F, the voxel functions f _n (y) satisfy the condition that for each n,

is minimized, so the voxel functions are as close as possible to the ideal voxel functions (with this measure of the difference between them) subject to the constraint that they are expressed as linear combinations of the K phase encoding functions. Alternatively another criterion is used to choose F, for example the measure of closeness of the voxel functions to the ideal voxel functions uses absolute value instead of absolute value squared, or the measure is weighted in some way, or instead of minimizing the measure, it is only made adequately small by some criterion. But using Eq. (9) as a measure of closeness which is to be minimized has the potential advantage that F can be expressed in the simple form given by Eq. (8).

In the case where the gradient field has a linear gradient given by B(y) = Gy, and the ideal voxel functions are delta functions of y - y_n, the elements of F have a particularly simple form given by

F_n,k = exp[iγ-kGT_hs-y_n] (10)

In this special case, solving Eq. (3) may give the same result as taking a fourier transform of the data points πik to find the voxel densities V_n. This form of F is referred to herein, including in the claims, as a fourier transform matrix. In this case, finding the voxel densities may be done with a considerably faster computation by using FFT software to calculate the fourier transform, rather than by using general matrix multiplication software which does not take advantage of the special relations between the elements of F given by Eq. (10).

Optionally, at 312, the voxel functions are found using Eq. (5). In the case of a uniform field gradient, for example, the voxel functions are sine functions of γ-k_maxGT-(y - y_n), where k_max is the maximum value of k for the data points.

FIG's. 4, 5, 6, and 7 illustrate one-dimensional image reconstruction using a reconstruction matrix F found by this procedure for the nonlinear gradient field profile 402 shown in FIG. 4, in an imaging region 404 extending from y_mj_n to y_max. The density profile 406 is flat in a central portion 408 of imaging region 404, and zero elsewhere.

In the example illustrated here, the time interval T_hS that the gradient field is turned on is chosen so that φi(y) varies by 2π over the width of imaging region 404,

Using a value of T of about this magnitude has the potential advantage of avoiding aliasing within the imaging region, which could occur if larger T were used, and not spending a large part of the time sampling data from parts of k space that do not contribute very much to the reconstructed image, which could occur if much smaller T were used. The number N of voxels is 60, with grid points y_n uniformly distributed in the imaging region, and delta functions δ(y - y_n) used as the ideal voxel functions. The number K of data points is 11, going from k = -5 to k = +5. FIG's. 5 A and 5B show plots 500 of some of the voxel functions f_n(y), for n = 6, 12, 18, ... 60. The other voxel functions look similar, each one peaked near y = y_n. Although f₆₀(y) is not well localized, being peaked near both the left and right edges of the imaging region, this voxel does not contribute very much to the reconstructed density, which is very low near the edges. The other voxels resemble the sine functions that one would find for a linear gradient field, but with the width of the function depending on the local field gradient dB/dy, being wider where the gradient is smaller. FIG. 6 shows a plot 600 of the data points m^ for k = 0 to 5. The values for k = -1 to -5 are the same as for k = +1 to +5, by symmetry. FIG. 7 shows a plot 700 of the reconstructed image density, as a superposition of the voxel functions f_n(y), each one weighted by the voxel density V_n. The reconstructed image density resembles the actual density 406 shown in FIG. 4, but the edges are not as sharp because of the limited number of points in phase space used to represent the image.

In general, the method outlined in FIG. 3 for finding the reconstruction matrix F can be used for two-dimensional imaging as well, but with each one-dimensional voxel function and ideal voxel function, replaced by two-dimensional functions, of y and z, for example, and the one-dimensional array of data points and encoding functions, indexed by k, replaced by two-dimensional encoding functions, and a two-dimensional array of data points, indexed by k_y and k_z. However, for a large number of grid points and a large number of k's in each direction, the matrixes F, E, and C can be enormous, and take a long time to calculate. In practice, if the gradient fields in the two phase encoding directions, y and z, are each fairly independent of the other phase encoding coordinate, then to good approximation the two-dimensional voxel functions f_nyjnz(y,z) can each be taken to be the product of a y voxel function which depends only on y and the y grid index ny, and a z voxel function which depends only on z and the z grid index nz. Similarly, the two- dimensional encoding functions e_ky,kz(y,z) may be approximated by a product of a y encoding function e_(y) that depends on y and the index ky, and a z encoding function β(_Z) that depends only on z and the index kz. Then the two-dimensional array of data points in phase space, m, is given by

"Vfe = P^{*" d}y [^{m dz} p(y> ^z)^e{_y)ky M^e(z)fe M 02) where p(y,z) is the two-dimensional density. The two-dimensional array of voxel densities V_ny,nz is given in matrix form by

V = Fy m Fz^τ (13) where Fy and Fz are separate one-dimensional reconstruction matrixes in y and z. They are each found according to the procedure described in FIG. 3, using respectively the y gradient field profile B_(y)(y) which is taken to be independent of z, and the z gradient field profile B(_Z)(z) which is taken to be independent of y.

FIG. 8 shows a simulated oval-shaped test object 802, of uniform density, in a square imaging region extending between z_min and z_maχ, and between y_min and y_max, used to illustrate two-dimensional image reconstruction using Eq. (12). The y gradient field is assumed to have the same profile as in FIG. 4, and to be independent of z, while the z gradient field is assumed to have the same profile and to be independent of y. The same number of grid points, and the same number of k values is used in each of y and z, and the voxel functions and reconstruction matrixes in y and z have the same forms as the one- dimensional voxels shown in FIG's. 5A-5B, and as the one-dimensional reconstruction matrix used to produce the constructed image in FIG. 7. The resulting two-dimensional reconstructed image 900 is shown in FIG. 9. Like one-dimensional reconstructed image 700 in FIG. 7, reconstructed image 900 has a boundary that closely follows the boundary of the original object being imaged, but its edge is not as sharp, and it has some small oscillations near the edge, because of the limited range of ky and kz of the encoding functions used to approximate it.

Table 1 shows the values of simulated data samples nikykz calculated from Eq. (15), using the density p(y,z) shown in FIG. 6. Table 1 only lists values for ky = 0 to 5, and kz = 0 to 5, because, due to the symmetry of p(y,z) in y and z, all the πi_{ky kz}'s are real, and m_ky_;kz

Table 1. Data samples mk_y,kz kz = 5 0.4 0.4 0.2 0.3 0.2 0.3 kz = 4 0.6 0.4 0.3 0.3 0.2 0.2 kz = 3 0.7 0.3 0.4 0.3 0.3 0.2 kz = 2 0.7 0.1 0.7 0.1 0.3 0.3 kz = l 0.6 0.7 0.6 0.6 0.4 0.2 kz = 0 5.3 2.0 1.3 0.6 0.1 0.4 ky = 0 ky = l ky = 2 ky = 3 ky = 4 ky = 5

An examination of Table 1 shows that even with the limited number of data samples, most data samples mky^_z are relatively small in amplitude and contribute little to the reconstructed image, with most of the contribution to the image coming from a few larger components of m_ky^_z with ky and kz smallest in absolute value. In some embodiments of the invention, data samples taken at larger absolute values of ky and kz use lower RF power, with the nominally 90 degree excitation RF pulses, and/or the nominally 180 degree refocusing RF pulses, having actual flip angles substantially less than 90 degrees and 180 degrees respectively. Although using this option reduces the SNR for those data samples, it has very little effect on the appearance and the SNR of the reconstructed image. Table 2 lists the amplitudes of simulated data samples m_κy,κz taken with reduced RF power for larger absolute values of ky and kz, rounded to the nearest 0.05. For each (ky, kz), the flip angle α for the refocusing pulses had sin(α/2) reduced from 1 (which means α = 180 degrees) to approximately exp[-0.13(ky² + kz²)], while the nominally 90 degree excitation RF pulses, which generally contribute very little to the total RF power, were kept at 90 degrees.

Table 2. Data samples m_{κy κz} with reduced RF power kz = 5 0 0 0 0 0 0 kz = 4 0.10 0.05 0 0 0 0 kz = 3 0.30 0.10 0.10 0.05 0 0 kz = 2 0.50 0.05 0.35 0.05 0 0 kz = l 0.55 0.60 0.35 0.25 0.05 0 kz = 0 5.30 1.80 0.90 0.25 0 0 ky = 0 ky = 1 ky = 2 ky = 3 ky = 4 ky = 5

FIG. 9 shows a plot of density for a reconstructed image 902, using these values of m_Ky,κz- There is very little difference between image 902, with reduced RF power, and image 900. The total RF power in the simulated data acquisition for image 802 was 12 times lower than for image 900.

Optionally, the reduced flip angles used, particularly for the nominally 180 degree refocusing pulses which generally comprise most of the RF power, are chosen so that the RF power is approximately as low as it can be, while having little expected effect on the quality of the images. Optionally, the expected effect on the quality of the images, for a given set of flip angles, is estimated by calculating simulated images having characteristics, such as contrast and scale of structure, similar to that in actual medical images, for the application that the probe is intended to be used for, for example looking for prostate cancer, or looking for plaque in arteries.

Adjustment in positions of voxels

Once a set of voxel densities is found for a set of voxel functions, whether using one of the methods described above, or just using a conventional method such as FFT, the voxel functions will not generally have peak values at a set of points on an equally spaced a rectangular grid, but, due to the nonlinearity and non-separability of the y and z gradient fields, the voxels will typically be localized on an array of curved lines of constant y gradient field and constant z gradient field, on an imaging slice of constant resonance frequency. This may be seen in plot 1006 in FIG. 10, which shows an imaging slice with isophase curves 1002 and 1004 for the y gradient and z gradient fields. To the extent that the gradient fields are small compared to a static magnetic field, the isophase curves on the imaging slice closely follow curves where the component of the gradient field parallel to the static magnetic field is constant.

In some embodiments of the invention, in order to produce a more accurate display of the reconstructed image, the positions of each rectangular grid point labeled by (ny, nz) associated with a voxel density V_ny,_nz is adjusted to reflect the actual location of the grid point where the y and z gradient fields have a given value. The element V_{ny nz} is moved from its old assumed position (y_ny, z_nz), in which the y coordinate depends only on ny and the z coordinate depends only on nz, and is assigned to a new position (y'_ny,nz, z'_ny,nz)- The new grid points do not, in general, form a rectangular grid, with the y coordinates depending only on ny and the z coordinates depending only on nz, but the y and z coordinates of the points each depend on both ny and nz. The new grid points are related to the old grid points by solving two simultaneous equations for each grid point (ny,nz), to find the two unknown coordinates y'_ny,nz and z'_{ny nz} of the new grid point:

By Vny,nz ' ^Zny,nz ) ^~ ^_y {} '_ny ) ( 14a)

Bz Vny,nz > ^Zn'y,nz ) ⁼ ^z \^Znz ) (14b) where B_y(y) and B_z(z) are the approximate values for the y-gradient and z-gradient fields used to find the separate y-direction and z-direction reconstruction matrixes used in Eq. (12). The new grid point is the point where the actual gradient fields B_y(y,z) and B_z(y,z) have the same values as the approximate gradient fields B_y(y) and B_z(z) have at the old grid point. Another method of two-dimensional image reconstruction, used in some embodiments of the invention, is similar to the method described above using Eq. (12), but non-Cartesian coordinates Y and Z are used from the beginning, instead of y and z. Coordinates Y and Z are chosen so that the nominally y-gradient field B_y(y,z) is a function only of Y, not of Z, and the nominally z-gradient field B_z(y,z) is a function only of Z, not of Y. Then the method using Eq. (12), and described above, is carried through, but using Y and Z everywhere instead of y and z. The end result is a two-dimensional density profile depending on Y and Z. Optionally, the density profile is then converted to Cartesian coordinates y and z, and displayed.

In some embodiments of the invention, the density is corrected to take into account the varying voxel volume, due to the nonlinearity and non-separability of the gradient fields in y and z. Although this may be done automatically if a fully two-dimensional reconstruction matrix were used to find reconstruct the image, such a correction can also be done if separate one-dimensional reconstruction matrixes are used for y and z, as in Eq. (12), to correct for the non-separability of the y and z gradient fields, as well as if a method such as FFT is used, ignoring the nonlinearity of the field gradients, which may be feasible for field gradients that are only mildly nonlinear. As may be seen from FIG. 10, the volume of a voxel may depend both on the spacing between adjacent curves of constant y gradient field and constant z gradient field at the location of the voxel, and it will also, in general, depend on the angle that the curves of constant y gradient field make to the curves of constant z gradient field. If the resulting differences in voxel size are not taken into account, then the image will show artifacts, with apparently higher density in regions where the voxel size is bigger due to these effects.

Regardless of which method is used to reconstruct the two-dimensional image, it should be noted that the image will generally be more accurate if the gradient fields B_y(y,z) and B_z(y,z) are such that, within a region of an imaging slice from which substantial signals are received, there is at most one point (y,z) for a given set of gradient field values (B_y, B_z). If this is not the case, then the gradient fields will not be able to uniquely identify which point a given component of the signal comes from, and this may result in image artifacts. FIG. 10 schematically shows an example of contours 1002 of constant B_y, and contours 1004 of constant B_z, as a function of y and z in an imaging slice 1006. These are also respectively contours of constant Y and constant Z, in the second two-dimensional image reconstruction method described above. Near the corners of the imaging slice, some contours 1002 intersect some contours 1004 twice, for example at point 1008 and at point 1010. This may cause imaging artifacts, unless at least one of the intersection points, for example point 1010, is effectively outside the imaging region, for example because the RF field is too low there, or because the RF antenna is too insensitive to signals from that point, or because there are not enough nuclei to be excited there.

FIG. 12 shows an MRI system 1100 which uses any of the methods described above, singly or in combination, for reconstructing images. System 1 100 comprises an MRI probe 1102, optionally a self-contained MRI used inside the body for medical imaging, such as the prostate probe described in the co-filed applications "MRI Probe" and "MRI Magnet and Coil Configurations." Such probes are likely to use substantially nonlinear gradient fields, and/or have relatively low SNR, so would be especially likely to benefit from these image reconstruction methods, but the methods could be used in any MRI system using phase encoding, including a conventional MRI system using a bore magnet. A controller 1104 controls the pulse sequences of probe 1 102 and receives and analyzes data from probe 1102. Alternatively, these functions could be performed by separate elements of system 1100. Controller 1104 is, for example, a computer with a control interface and data interface to probe 1102. Controller 1104 is programmed, in its role as data analyzer, to reconstruct images from probe 1102, using any of the methods described above, and has a sufficient CPU, memory, and storage space to perform this task in a reasonably short time, optionally in real time.. Optionally, an output device 1 106, such as a computer monitor or a printer, linked to controller 1104, provides output images from probe 1102, optionally in real time. Alternatively, an image file, for example an image file in some standard format, is stored, and/or transferred to another location, and used to provide an image later. Optionally, an input device 1 108, such as a keyboard or mouse linked to controller 1104, is used by an operator of the system, such as a physician performing a prostate examination and/or biopsy with the system, to instruct the controller when to initiate the imaging procedure, and optionally to set one or more free parameters of the imaging procedure.

In some embodiments of the invention, system 1 100 does not include probe 1102, but controller 1104 obtains data from an MRI probe and analyzes the data, according to any of the methods described above. As used herein the term "about" refers to ± 10 % .

The terms "comprises", "comprising", "includes", "including", "having" and their conjugates mean "including but not limited to". This term encompasses the terms "consisting of and "consisting essentially of. The phrase "consisting essentially of means that the composition or method may ^ include additional ingredients and/or steps, but only if the additional ingredients and/or steps do not materially alter the basic and novel characteristics of the claimed composition or method.

As used herein, the singular form "a", "an" and "the" include plural references unless the context clearly dictates otherwise. For example, the term "a compound" or "at least one compound" may include a plurality of compounds, including mixtures thereof.

Throughout this application, various embodiments of this invention may be presented in a range format. It should be understood that the description in range format is merely for convenience and brevity and should not be construed as an inflexible limitation on the scope of the invention. Accordingly, the description of a range should be considered to have specifically disclosed all the possible subranges as well as individual numerical values within that range. For example, description of a range such as from 1 to 6 should be considered to have specifically disclosed subranges such as from 1 to 3, from 1 to 4, from 1 to 5, from 2 to 4, from 2 to 6, from 3 to 6 etc., as well as individual numbers within that range, for example, 1, 2, 3, 4, 5, and 6. This applies regardless of the breadth of the range.

Whenever a numerical range is indicated herein, it is meant to include any cited numeral (fractional or integral) within the indicated range. The phrases "ranging/ranges between" a first indicate number and a second indicate number and "ranging/ranges from" a first indicate number "to" a second indicate number are used herein interchangeably and are meant to include the first and second indicated numbers and all the fractional and integral numerals therebetween.

It is appreciated that certain features of the invention, which are, for clarity, described in the context of separate embodiments, may also be provided in combination in a single embodiment. Conversely, various features of the invention, which are, for brevity, described in the context of a single embodiment, may also be provided separately or in any suitable subcombination or as suitable in any other described embodiment of the invention. Certain features described in the context of various embodiments are not to be considered essential features of those embodiments, unless the embodiment is inoperative without those elements.

Although the invention has been described in conjunction with specific embodiments thereof, it is evident that many alternatives, modifications and variations will be apparent to those skilled in the art. Accordingly, it is intended to embrace all such alternatives, modifications and variations that fall within the spirit and broad scope of the appended claims.

All publications, patents and patent applications mentioned in this specification are herein incorporated in their entirety by reference into the specification, to the same extent as if each individual publication, patent or patent application was specifically and individually indicated to be incorporated herein by reference. In addition, citation or identification of any reference in this application shall not be construed as an admission that such reference is available as prior art to the present invention. To the extent that section headings are used, they should not be construed as necessarily limiting.

Claims

WHAT IS CLAIMED IS:

1. A method of reconstructing an image from a set of MRI data points in phase space acquired from an imaging slice, using phase encoding with one or more gradient fields that are each nonlinear as functions of position in a phase encoding direction, the method comprising: a) multiplying a vector of the data points by a reconstruction matrix which differs from a fourier transform matrix at least due to taking into account the nonlinearity of the one or more gradient fields; and b) obtaining a vector of voxel densities for the image.

2. A method according to claim 1, wherein the reconstruction matrix is the product of an encoding matrix and the inverse of a correlation matrix of phase encoding functions for the phase encoding with the one or more nonlinear gradient fields, each element of the phase encoding matrix being an integral, over at least one phase encoding direction, of a product of one of the phase encoding functions and one of a set of ideal voxel functions, each of which is compact around a different grid point in at least one of the one or more phase encoding directions.

3. A method according to claim 2, wherein the ideal voxel functions are delta functions, boxcar functions, half-sin functions, or cubic spline functions of position in at least one of the one or more phase encoding directions.

4. A method according to any of claims 1-3, wherein the MRI data points consist of K points with phase encoding differing only in characteristics of a single gradient field pulse which varies in one phase encoding direction, the voxel densities consist of densities of N voxel functions of the one phase encoding direction, the reconstruction matrix is a K by N matrix, and the vector of the voxel densities is the product of the vector of the K data points, and the reconstruction matrix for the one phase encoding direction.

5. A method according to any of claims 1-3, where the imaging slice has two phase encoding directions, and the MRI data points are phase encoded by two gradient fields, one for each of the directions.

6. A method according to claim 5, wherein the MRI data points comprise an array of Kx by Ky points in two dimensional phase space, the voxel densities comprise an array of Nx by Ny densities for an array of two-dimensional voxel functions each of which is a product of one of a set of Nx voxel functions of x, a first of the two phase encoding directions, and one of a set of Ny voxel functions of y, a second of the two phase encoding directions, and transforming the data points comprises: a) multiplying a matrix of the Kx by Ky data points by a Kx by Nx reconstruction matrix for the first of the two phase encoding directions; and b) multiplying the product by a Ky by Ny reconstruction matrix for the second of the two phase encoding directions, to obtain an Nx by Ny matrix of the voxel densities.

7. A method according to claim 6, wherein the one-dimensional voxel functions of x have a minimum of a measure of difference from a set of Nx ideal compact voxel functions localized around an array of grid points, for any set of Nx one-dimensional functions that can be expressed as linear combinations of Kx phase encoding functions of x for Kx evenly spaced amplitudes of the gradient field for x, the reconstruction matrix for x being a linear transformation from the phase encoding functions to the voxel functions of x, and wherein the one-dimensional voxel functions of y have a minimum of a measure of difference from a set of Ny ideal compact voxel functions localized around an array of grid points, for any set of Ny one-dimensional functions that can expressed as linear combinations of Ky phase encoding functions of y for Ky evenly spaced amplitudes of the gradient field for y, the reconstruction matrix for y being a linear transformation from the phase encoding functions to the voxel functions of y.

8. A method according to claim 5, wherein the MRI data points comprise K data points, each characterized by a different pair (kx,ky), the voxel functions comprise N voxel functions, each peaked at a different pair of values (x,y) of the two dimensions x and y, and transforming the data points comprises multiplying a vector of the K data points by a K by N reconstruction matrix for the two dimensions, to obtain a vector of the N voxel densities in the two dimensions.

9. A method of calculating an image reconstruction matrix for MRI data acquired from an imaging slice using phase encoding with one or more gradient fields that are each a nonlinear function of position in a phase encoding direction, the method comprising: a) choosing a set of ideal voxel functions; b) calculating a set of phase encoding functions, one for each value of each of the one or more gradient fields at which the MRI data are acquired; and c) finding a set of voxel functions, each of which has a minimum of a measure of difference from one of the ideal voxel functions, subject to the constrain that the voxel function can be expressed as a linear combination of the phase encoding functions; wherein the reconstruction matrix is the linear transformation from the phase encoding functions to the voxel functions.

10. A method according to claim 9, wherein finding the set of voxel functions comprises: a) calculating an encoding matrix, each element of which is an integral over the one or more phase encoding directions of one of the phase encoding functions and one of the ideal voxel functions; b) calculating a correlation matrix of the phase encoding functions; c) multiplying the encoding matrix by the inverse of the correlation matrix, thereby obtaining the reconstruction matrix which is a linear transformation from the phase encoding functions to the voxel functions.

1 1. A method of reconstructing an image in an imaging slice from a set of MRI data points in phase space acquired using phase encoding by one or more gradient fields, at least one of which has a field gradient that is non-uniform in magnitude or direction or both, in the imaging slice, the method comprising: a) transforming the set of data points to a set of voxel densities, each voxel being identified by a different set of values of the one of more gradient fields; b) for each voxel, determining a position in the imaging slice where the one or more gradient fields have the values in the set identifying that voxel; and c) reconstructing the image by calculating a density as a function of position in the imaging slice, using the voxel densities and the positions in the imaging slice of the voxels.

12. A method according to claim 11, wherein calculating a density as a function of position in the imaging slice comprises adjusting one or more of the voxel densities responsive to an estimated volume of the voxel calculated from a gradient of each of the one or more gradient fields at the position of the voxel.

13. A method of generating an MRI image of an imaging region with reduced noise, comprising: a) applying to one or more slices of the imaging region NMR pulse sequences suitable for separating even parity and odd parity echoes, for a plurality of points in phase space; b) accumulating NMR signals for each point in phase space for each slice, and separating the accumulated signals into a set of even parity data points and a set of odd parity data points in phase space; c) identifying an expected phase of a signal part of each data point; d) subtracting from each data point a part that is expected to be substantially noise because it has phase orthogonal to the expected phase of the signal part; and e) using the sets of even parity and odd parity data points after subtracting the parts expected to be noise, to reconstruct an image with reduced noise.

14. A method according to claim 13, wherein the NMR pulse sequences comprise, for a point in phase space for a slice, a first NMR pulse sequence consisting of an excitation RF pulse, followed by one or more gradient pulses, followed by a train of refocusing pulses with substantially zero phase difference from the excitation pulse, and a second NMR pulse sequence consisting of an excitation pulse, followed by one or more gradient pulses, followed by a train of refocusing pulses with substantially 90 degree phase difference from the excitation pulse, and separating the signals into sets of even and odd parity data points comprises finding a sum of the accumulated signals from the first and second pulse sequences for each data point, and finding a difference of the accumulated signals from the first and second pulse sequences for each data point.

15. A method according to claim 13 or claim 14, wherein using the sets of even parity and odd parity data points to reconstruct an image comprises adding data of one parity to data of the other parity that has been flipped in each phase encoding direction, either in phase space or in position space.

16. A method according to claim 13, wherein identifying an expected phase of the signal part of a data point comprises finding phases of one or more data points for a plurality of different resonance frequencies, fitting the phases to a function of resonance frequency, and identifying the value of the fitted function of resonance frequency as the expected phase for all data points in phase space for that resonance frequency.

17. A method of MRI imaging, comprising: a) applying an NMR pulse sequence to one or more imaging slices of an imaging region; b) acquiring a set of MRI data points from the one or more imaging slices; and c) reconstructing an image of the imaging region using any of the methods of claims 1-8 or 1 1-16.

18. An MRI imaging system, comprising: a) a self-contained MRI probe comprising a magnet assembly that produces a static magnetic field in an imaging region, one or more sets of gradient coils that each produce a magnetic gradient field in the imaging region, and one or more RF antennas that together produce RF fields that excite nuclei in the imaging region and receive NMR signals from the excited nuclei; and b) a controller that directs the probe to produce pulse sequences of gradient fields and RF fields in the one or more slices of the imaging region, and that accumulates NMR signals to produce MRI data points in phase space, and that reconstructs an image from the MRI data points according to any of the methods of claims 1-8 or 11-16.