The present invention is directed to imaging using a spiral scan and more particularly to such imaging that provides rapid and accurate off-resonance correction.
Spiral scanning is extensively used in magnetic resonance (MR) imaging and has a number of desirable properties, such as short scan time, resistance to motion artifacts, capability to achieve short echo time, and excellent gradient efficiency. One major limitation, however, is image blurring due to off-resonance effects from main field inhomogeneity and tissue-induced susceptibility variations. This problem is more pronounced with long readout times and at high field strength.
Many image deblurring methods have been developed to correct for off-resonance effects in spiral imaging.
The linear off-resonance correction method is a widely used method that requires little computation and is robust at low signal-noise ratio (SNR) regions. However, this technique assumes a linear field map and therefore residual blurring can manifest at the locations where the field value deviates from the linear variation. Image domain deconvolution approximates the off-resonance phase as a separable quadratic function, which allows image deblurring to be carried out as rapid one-dimensional deconvolutions. A drawback of this method is that it can become ineffective when the deviation of the off-resonance phase from the quadratic function cannot be ignored. Conjugate phase reconstruction is an effective but computationally expensive off-resonance correction method. Fast alternatives, including time-segmented reconstruction, frequency-segmented reconstruction, and multifrequency interpolation, and variants have been developed. A fundamental assumption for conjugate phase reconstruction and its alternatives is that the magnetic field varies slowly in space. Alternative image reconstruction methods have been developed for off-resonance correction when this assumption is violated.
All of the previously described pixelwise off-resonance correction methods require a reliable knowledge of the field map that indicates the resonant frequency of spins in each image voxel. Typically a field map is derived from two images acquired by the same spiral scan at different echo times. Field maps acquired by this method can be inaccurate, since the field map itself is susceptible to the off-resonance effect. Another problem associated with this field map acquisition method is that often only a low resolution field map can be acquired in the scan time available. Partial volume effects may cause a low resolution field map to be not accurate enough for a pixelwise off-resonance correction in certain cases. For field map based off-resonance correction methods, the image blurring can remain or even be exaggerated at regions where the estimated off-resonance frequency value deviates from the actual off-resonance frequency.
Automatic off-resonance correction methods are alternative image deblurring methods for non-Cartesian imaging. These methods can estimate a field map from the imaging data set itself and therefore perform off-resonance correction without acquiring a field map. Based on the principle that image is real when reconstructed on resonance, Noll et al (Noll D C, Pauly J M, Meyer C H, Nishimura D G, Macovski A. Deblurring for Non-2DFT MRI. Magn Reson Med 1992; 25: 319-333; see also U.S. Pat. No. 5,311,132) first proposed an automatic correction method that determines a field map by minimizing an objective function that corresponds to the imaginary components of the image. Noll's method is prone to estimation errors when the target field map resolution is high, causing image artifacts in certain cases. To address this problem, Man et al (Man L C, Pauly J M, Macovski A. Improved automatic off-resonance correction without a field map in spiral imaging Magn Reson Med 1997; 37: 906-913) proposed to perform Noll's automatic correction method in two stages: first to estimate a low resolution field map and then to estimate a high resolution field map based on the searching frequency constraint from the low resolution field map. A modified version of this two-stage method is to first estimate a linear map using a fast automatic method and then employ it as a frequency constraint for a second stage automatic off-resonance correction. A common problem for all of these automatic methods, however, is that they are computationally inefficient, which can become impractical when processing a large volume data set.
It is therefore an object of the invention to provide a computationally efficient method.
To achieve the above and other objects, the present invention is directed to a new off-resonance correction algorithm and related method and system for spiral imaging. The present invention is similar to the previously reported two-stage automatic off-resonance correction method, but has some major differences. In the first stage, we acquire a low resolution map using two single shot spirals instead of estimating it using an automatic method. This modification saves computation, and avoids potential estimation errors that can occur when using an automatic method with inappropriate parameters. In the second stage, we developed a modified version of automatic off-resonance correction with a frequency constraint from the first stage map. The second stage automatic deblurring is more computationally efficient than previous automatic methods and is insensitive to the variation of parameters required in an automatic method. Phantom and volunteer experiments show that the proposed method can provide improved off-resonance correction for spiral imaging compared to the existing map based off-resonance correction methods.
The present invention provides improved off-resonance correction in spiral imaging. In our method, a low resolution field map is acquired by two single shot spirals and a computationally efficient automatic method is performed for off-resonance correction with a frequency constraint from the acquired low resolution field map. Two different strategies to implement this frequency constraint, a linear map constraint and a polynomial map constraint, were investigated. Compared to the previously reported map based off-resonance correction methods, our experiments indicate that the proposed method can provide off-resonance correction with improved accuracy in spiral imaging.
The present invention provides an efficient off-resonance correction method and related system for spiral imaging with improved accuracy. We use two single shot spirals to collect a low resolution field map and then perform a computationally efficient automatic method for off-resonance correction with a frequency constraint from the acquired low resolution field map. The present invention can be used for online reconstruction and can provide more accurate off-resonance correction in spiral imaging than field map based off-resonance correction methods.
A preferred embodiment of the invention will be disclosed with reference to the drawings, in which:
A preferred embodiment of the present invention will be set forth in detail with reference to the drawings, in which like reference numerals refer to like elements or steps throughout.
The MR signal, ignoring relaxation, from an object with spin density m(r) is given by
s(t)=∫m(r)exp{−i[2πk(t)·r+Δω(r)t]}dr [1]
where r is the spatial position, k(t) is the k-space trajectory, and Δω(r) is the spatial distribution of the off-resonance.
In the common situation where the spatial variation of the off-resonance frequency is relatively slow, conjugate phase reconstruction is an effective way to remove the image blurring induced by the off-resonance effects in non-Cartesian imaging. Conjugate phase reconstruction can be expressed as following:
m(r)=∫s(t)W(t)exp{i[2πk(t)·r+Δω(r)]}dt, [2]
where W(t) is the density compensation accounting for the non-uniform sampling in the k-space for non-Cartesian data acquisition. Note that a priori knowledge of the off-resonance map Δω(r) is necessary in conjugate phase reconstruction in order to achieve accurate off-resonance correction. Conjugate phase reconstruction is computationally expensive. Fast alternatives, including time-segmented reconstruction, frequency-segmented reconstruction, and multifrequency interpolation, can achieve comparable accuracy of image reconstruction but with significantly reduced computational cost.
Automatic off-resonance correction methods can be performed without a priori knowledge of the off-resonance map. In the original automatic off-resonance correction method proposed by Noll et al, the time signal is demodulated by a set of equally spaced constant frequencies, {Δωi, i=1, 2 . . . }, and complex images are reconstructed from these demodulated signals, which can be expressed as follows:
m(r;Δωi)=∫s(t)exp(iΔωit)W(t)exp[i2πk(t)·r]dt, [3]
where m(r; Δωi) represents the complex image reconstructed at frequency Δωi. Eq. [3] can be performed by a fast image reconstruction method, such as gridding. After calculating the complex images, an objective function is then calculated pixel by pixel for automatic evaluation of the field map. The objective function is typically chosen as the sum of the absolute value of the imaginary image raised to a certain power within a summation window around the point of interest. The local off-resonance frequency is chosen as the constant frequency which minimizes the objective function calculated at that location. The final deblurred image can be either formed by piecing up the pixels from the complex images or reconstructed using a fast conjugate phase reconstruction method with a smoothed version of the estimated field map.
The automatic method can be better understood by the concept of the point spread function (PSF). As shown in the appendix, the complex image m(r;Δωi) can also be equivalently expressed by the following formula:
m(r; Δωi)=∫mtrue(r′)psf(r|r′; Δωi)dr′ [4]
where mtrue(r′) is the real-valued physical image, and psf (r|r′;Δωi) is the PSF evaluated at location r given that the input impulse is at location r′ and the demodulation frequency is Δωi. Eq. [4] indicates that the complex images m(r;Δωi) can be represented as the sum of product of the physical image and the PSFs with the location of input impulse varying throughout FOV. From the mathematical expression of psf(r|r′;Δωi), we can tell that when the demodulation frequency Δωi is tuned to the exact off-resonance frequency value at the point of interest, the PSF with the input impulse located at the same location will become a real-valued impulse located at that point. The PSFs with input impulse located at nearby locations can be approximated as real impulses when the spatial variation of the field is slow. The PSFs with input impulse located far away from the point of interest may have a significant imaginary component; however, they will decay away at the location of the point of interest when the field varies smoothly. Consequently, the image value at the point of interest calculated from Eq. [4] approximately equals the value of the real-valued physical image at that point. Conversely, if the demodulation frequency Δωi is tuned away from the exact off-resonance frequency at the point of interest, the PSFs at the point of interest and also at the neighboring locations will spread out and have both real and imaginary energy, which results in the image value at the point of interest calculated from Eq. [4] becoming complex. Consequently, the criterion that the complex image m(r;Δωi) should have minimal imaginary energy can be employed to determine whether or not the demodulation frequency Act), is close to the actual off-resonance frequency.
The imaginary energy of the PSFs can decrease if the demodulation frequency is tuned far away from the actual off-resonance frequency, which can induce spurious minima in the objective function in automatic methods. Increasing the size of the summation window used for objective function calculation may reduce the possibility of encountering spurious minima; however, an increased size of the summation window may also decrease the accuracy of the estimated field map. To address this problem, Man et al proposed performing automatic off-resonance correction in two stages. The first stage is to use a large summation window to estimate a low resolution field map using an automatic method. The second stage is to use a small summation window to estimate a high resolution field map with the searching frequency constrained by the low resolution field map. Since spurious minima usually occur when the demodulation frequency is tuned far away from the exact off-resonance frequency, by constraining the searching range around the low resolution field map, the possibility of encountering spurious minima of the objective function is reduced. The present inventors (Chen W, Meyer C R Fast automatic linear off-resonance correction method for spiral imaging Magn Reson Med 2006; 56: 457-462) previously proposed a modified version of this two-stage automatic off-resonance correction method where a linear map is estimated using fast automatic method and then is used as frequency constraint for the second stage automatic deblurring. One major limitation of all of these automatic methods; however, is computational inefficiency, which makes them impractical in processing large datasets, such as 3D volume and multi-coil data sets.
The algorithm according to the preferred embodiment is compared to previously described two-stage automatic off-resonance correction methods but has some important differences. In the previously described two-stage automatic off-resonance methods, a low resolution field map or a linear map is estimated by the automatic method itself. A problem associated with these methods is that an inappropriate choice of a parameter value required in automatic method, such as the amount of incidental phase to be removed, can lead to estimation error in the first stage map. A low resolution field map can also be acquired efficiently by collecting two single shot spirals with different echo times during the approach to steady state. In the present algorithm, we employ this method to acquire a low resolution field map. This strategy saves computation time and the accuracy of acquired low resolution map is sufficient to provide a frequency constraint for the second stage automatic deblurring.
In the second stage, we perform pixelwise automatic off-resonance correction with searching frequencies constrained by the first stage map. The constraint from the first stage map allows us to use a small summation window in the second stage automatic deblurring to achieve high resolution off-resonance correction. Unlike previous two-stage automatic off-resonance correction methods, linear off-resonance correction is incorporated into image reconstruction in the present algorithm, which can be expressed as following:
m(r;Δω′i)=ωs′(t)exp(iΔω′it)W′(t)exp(i2π(k′x(t)x+k′y(t)y)dt [5]
with
s′(t)≡s(t)·exp(i2πf0t)
k′x(t)≡kx(t)+αt,
k′y(t)≡ky(t)+βt
where f0, α and β are the center frequency, x and y gradient of the extra acquired low resolution map, respectively, W′(t) is the density compensation function after considering the gradient of field inhomogeneity, and Δω′i(i=1, 2 . . . ) are a series of constant frequencies used to search for the non-linear components of the off-resonance frequency. Note that Δω′i(i=1, 2 . . . ) in Eq. [5] defines a searching range of only the non-linear components of the off-resonance frequencies. Conversely, the previous automatic methods search for the actual off-resonance frequency directly. This modification can significantly improve the computational efficiency of our algorithm, which will be discussed in more detail below.
After calculating images in Eq. [5], the objective function is then calculated as:
ε(r;Δω′1)=∫H(r)|Im{m(r;Δω′i)}|αdr, [6]
where H(r) and α represent the summation window and a constant power, respectively. A small summation window can be used here to achieve high resolution off-resonance correction because the possibility of encountering a spurious minimum of objective function is reduced with the frequency constraint from the acquired linear map. The map of the non-linear component of off-resonance frequency is determined as:
Δω′L(r)(r), such that L(r)=arg miniε(r;Δω′i). [7]
After acquiring Δω′L(r)(r), we smooth it using a small window and the final image is reconstructed using Eq. [5] after replacing Δω′i with smoothed Δω′L(r)(r).
When the field inhomogeneity has a significant linear gradient, the PSFs in Eq. [4] can be skewed, which may reduce the effectiveness of automatic methods. One advantage of the present method is that this effect is inherently corrected in Eq. [5]. The proposed method also allows us to use a recursive method to calculate the objective function, which is impossible to implement in the previous two-stage automatic methods.
In certain cases, when the field inhomogeneity is highly non-linear, we can estimate a higher order polynomial map from the acquired low resolution field map and use it to constrain the second stage automatic off-resonance correction. A smooth polynomial approximation to the low resolution field map can be estimated rapidly by minimizing the following weighted least square norm:
∥|m(r)|(Δωp(r)−Δω(r))∥2, [8]
where |m(r)| represents the magnitude of the low resolution image which is reconstructed from the map data sets; Δωp(r) and Δω(r) represent the polynomial approximation and the low resolution map, respectively.
When using a polynomial map as a frequency constraint for the second stage automatic deblurring, the image reconstruction can be expressed as follows:
m(r;Δω″i)=∫s(t)exp{j[(Δω″i+Δωp(r)]t}W(t)exp[i2πk(t)·r]dt, [9]
where k(t) and W(t) are same as those defined in Eq. [2], and {Δω″i, i=1, 2 . . . } is a series of constant frequencies used to search for the deviation of the actual off-resonance frequency from the acquired polynomial map. Note that similar to Eq [5], in Eq. [9] we search for the deviation of the actual off-resonance frequency from the first stage map instead of the actual off-resonance frequency itself. For Eq. [9], calculating images at {Δω″i, i=1, 2 . . . } can be achieved by using multifrequency interpolation or the Taylor approximation method given a spatially varying polynomial map Δωp(r). After calculating the images in Eq. [9], the objective function can be calculated as described by Eq. [6] and the deviation of the actual off-resonance frequency from the first stage polynomial map can be determined as described in Eq. [7]. For polynomial map constrained automatic off-resonance correction, the objective function calculation can also be performed using the recursive method described in the inventors' previous paper cited above to improve the computational efficiency.
Multifrequency interpolation requires reconstructing a series of base images at different demodulation frequencies. To reduce the computational cost, we can interchange the order of demodulation and gridding for base image reconstruction. A strategy we employed is to use Delaunay triangulation to generate a sampling time mask on Cartesian grids and then to perform demodulation directly on the k-space data after gridding. When generating the time mask, the sampling time on the corners of Cartesian grids is set as total readout time. When implementing this strategy in multifrequency interpolation, gridding still needs to be performed for data in a small region close to the origin in k-space. To reduce the computation time in off-resonance correction based on a polynomial map constraint, the time mask can be calculated offline and loaded during online reconstruction for a given k-space trajectory. For off-resonance correction based on a linear map constraint, the time mask needs to be calculated online since the k-space trajectory is warped by the linear gradient. However, in our experiments, we usually found that the time mask generated offline can still be used for off-resonance correction based on a linear map constraint.
As discussed in previous works, the incidental phase caused by sources other than off-resonance effects must be removed when performing automatic off-resonance correction. We use the method described in our previous paper cited above to remove the incidental phase. In our experiments, we found a benign property of the present off-resonance correction method is that the second stage automatic deblurring is insensitive to variations of parameters used to determine the amount of incidental phase to be removed. A practical note is that when incorporating linear off-resonance correction in automatic deblurring, as described in Eq. [5], the incidental phase should be estimated from the demodulated signal s′(t) rather than the raw signal s(t).
Compared to the previous automatic methods, we have developed the following strategies to reduce the computational cost of the proposed algorithm. First, since the first stage field map is acquired from extra data instead of by automatic methods, there is little computational cost for first stage field map estimation. Second, in the present method, we use the automatic method to search for the deviation of the actual off-resonance frequency map from the first stage map. Conversely, in previous methods, the automatic method was used to search for the actual off-resonance frequency directly.
Various hardware can be used. For example, as shown in
When implementing the proposed off-resonance correction method, we typically use the following parameters: the size of summation window was 12×12; the power a of the objective function was 1; the first 1.2 ms of the readout was used to estimate the incidental phase (14); the searching range of the deviation of the actual map from the first stage map was from −50 Hz to 50 Hz; and the increment of searching frequencies was 10 Hz.
We have performed field map based linear off-resonance correction, multifrequency interpolation, and the proposed method on the same data sets, and compared the results with each other. For some data sets, a high resolution field map was acquired by repeating two multi-shot spiral scans with echo delay 1 ms, which was then used to support map-based multifrequency interpolation off-resonance correction.
The images were overall less blurry after applying linear off-resonance correction. However, residual image blurring remained in many data sets. In some data sets, image blurring was exaggerated in local regions after applying linear off-resonance correction indicating strong non-linearity of the field inhomogeneity in theses regions. Multifrequency interpolation based on an extra acquired field map can provide high quality off-resonance correction for some of the data sets. However, for other data sets, residual image blurring remained after off-resonance correction using mutlifrequency interpolation based on either low resolution or high resolution field map, indicating the inaccuracy of the acquired field maps in local regions. The proposed off-resonance correction method was the most robust among these methods in providing high quality off-resonance correction. The linear map constrained automatic method was employed in most of our in vivo applications. The polynomial map constrained automatic method was found to provide more accurate off-resonance correction than the linear map constrained automatic method in certain applications, particularly phantom studies, due to strong non-linear field inhomogeneity.
The present off-resonance correction method is computationally efficient. The image reconstruction system of our 1.5T Avanto scanner uses an AMD Opetron™. Our online reconstruction program takes about 1 to 2 seconds and 2 to 3 seconds when using linear and polynomial map constraint, respectively, to process a single coil data set with the parameters specified in the method section.
When implementing automatic off-resonance correction, we need to specify the value of parameters including the size of the summation window, the power of the objective function, and the amount of the incidental phase to be removed. Incorrect choices of the value of these parameters may lead to image artifacts. The original automatic off-resonance correction method proposed by Noll et al is sensitive to variations of these parameters. A moderate change of the value of these parameters may induce spurious minima. A desirable property of our proposed algorithm is that it is insensitive to variation in these parameters. The value of these parameters was kept constant when applying our algorithm for our regular spiral scans. When reducing the number of spiral interleaves to less than five, we found we needed to increase the size of the summation window when implementing our method.
We have also implemented our method for regional deblurring when only certain sections of images are of clinical interest. Performing regional deblurring can further reduce the computational cost. The linear map constraint can have improved accuracy when applying our algorithm for regional deblurring since the linear approximation only applies to a local region.
As previously indicated by Man et al, off-resonance correction can reduce image blurring artifacts and increase spatial resolution in spiral imaging without losing SNR. We did not observe visible SNR change when applying our proposed algorithm for off-resonance correction.
In the present algorithm, we used two additional single shot spirals to acquire a low resolution field map. These single shot spirals are also used to set up the approach to steady state. A low resolution field map acquisition can be incorporated into imaging data acquisition by using a variable density spiral scan. Further study should be investigated whether the proposed method can be incorporated into this variable density spiral imaging method. We recently reported a rapid automatic method to estimate a linear map in our previous paper cited above. This method can be combined with the proposed method to perform rapid pixelwise automatic off-resonance correction without acquiring a low resolution field map.
Our experimental data sets were all acquired from 1.5T scanners. Image blurring for spiral imaging at 3T can be more severe since the off-resonance effects are usually stronger at 3T. One likely challenge in applying our algorithm at 3T is that the acquired low resolution field map may become less accurate, since the single shot spirals are susceptible to off-resonance effects. If this occurs, we may apply a rapid automatic center frequency correction on the two single shot spirals before calculating the low resolution map.
One limitation of the present method is that it assumes the field is slowly varying in space. Though this assumption is valid in common applications, it can be violated in certain regions in the body, such as near the sinus. The present algorithm is inapplicable at those regions.
Automatic off-resonance correction methods use the minimum of the imaginary component of the image as a focusing criterion. Flow can also induce phase in MR images. As discussed by Man et al, this is not a problem when applying automatic off-resonance correction for spiral imaging since it is resistant to motion artifacts. For spiral imaging off isocenter, concomitant gradients may also induce image phase and can cause image blurring artifacts.
While a preferred embodiment has been set forth in detail above, those skilled in the art who have reviewed the present disclosure will readily appreciate that other embodiments can be realized within the scope of the present invention. For example, numerical limitations are illustrative rather than limiting, as are disclosures of specific technologies and sources. Therefore, the present invention should be construed as limited only by the appended claims.
Automatic off-resonance correction involves image reconstruction at a set of constant frequencies {Δωi}. The signal is first demodulated by the constant frequencies Δωi, and then followed by a rapid image reconstruction method, such as gridding, to form an image. Here we derive how the reconstructed image is related to the PSF function.
The image reconstructed at demodulation frequency Δωi can be expressed as follows:
m(r,Δωi)=∫s(t)exp(iΔωit)W(t)exp[i2πk(t)·r]dt [A.1]
Substituting Eq. [1] into Eq. [A.1] and rearranging the order of the integrals, we have:
m(r,Δωi)=∫mtrue(r′)dr′∫W(t)exp{−i2πk(t)·(r′−r)−[Δω(r′)−Δωi]t}dt
Given that
where psf(r|r′;Δωi) is the PSF evaluated at location r given the input impulse is at r′ and the demodulation frequency is Δωi, we have
m(r,Δωi)=∫mtrue(r′)psf(r|r′;Δωi)dr′ [A.2]
The present application claims the benefit of U.S. Provisional Patent Application No. 60/903,021, filed Feb. 23, 2007, whose disclosure is hereby incorporated by reference in its entirety into the present disclosure.
The work leading to the present invention was supported by NIH Grant No. HL079110. The government has certain rights in the invention.
Number | Name | Date | Kind |
---|---|---|---|
5311132 | Noll et al. | May 1994 | A |
6396269 | Hajnal et al. | May 2002 | B1 |
7005854 | Mitchell et al. | Feb 2006 | B2 |
7285955 | Roberts et al. | Oct 2007 | B2 |
7486839 | Moriguchi et al. | Feb 2009 | B2 |
20040243169 | Ricci et al. | Dec 2004 | A1 |
20050017717 | Duerk et al. | Jan 2005 | A1 |
20050040821 | Hargreaves et al. | Feb 2005 | A1 |
Number | Date | Country | |
---|---|---|---|
60903021 | Feb 2007 | US |