Abstract:
A CE-MIRA image is acquired using a centric view order. The acquired k-space data set is motion corrected by reconstructing tracking images from successive annular segments of k-space data and detecting the location of a point spread object therein. Each k-space data segment is phase corrected to align the point spread object with a reference position.

Description:
BACKGROUND OF THE INVENTION 
     The field of the invention is nuclear magnetic resonance imaging methods and systems. More particularly, the invention relates to the correction of motion artifacts in MR images. 
     When a substance such as human tissue is subjected to a uniform magnetic field (polarizing field B 0 ), the individual magnetic moments of the spins in the tissue attempt to align with this polarizing field, but precess about it in random order at their characteristic Larmor frequency. If the substance, or tissue, is subjected to a magnetic field (excitation field B 1 ) which is in the x-y plane and which is near the Larmor frequency, the net aligned moment, M z  may be rotated, or “tipped”, into the x-y plane to produce a net transverse magnetic moment M t . A signal is emitted by the excited spins after the excitation signal B 1  is terminated, this signal may be received and processed to form an image. 
     When utilizing these signals to produce images, magnetic field gradients (G x  G y  and G z ) are employed. Typically, the region to be imaged is scanned by a sequence of measurement cycles, or “views,” in which these gradients vary according to the particular localization method being used. The resulting set of received NMR signals are digitized and processed to reconstruct the image using one of many well known reconstruction techniques. 
     Depending on the particular view order used, object motion during the acquisition of NMR image data can produce blurring or “ghosts” in the phase-encoded direction, or both. For most physiological motion each view of the NMR signal is acquired in a period short enough that the object may be considered stationary during the acquisition window. In such case the blurring and ghosting is due to the inconsistent appearance of the object from view to view. Motion that changes the appearance between views such as that produced by a patient moving, by the respiration or the cardiac cycle, or by peristalsis, is referred to hereinafter as “view-to-view motion”. Motion may also change the amplitude and phase of the NMR signal as it evolves during the pulse sequence and such motion is referred to hereinafter as “in-view motion”. 
     One method for reducing motion artifacts in NMR images is referred to in the art as “gradient moment nulling”. This method requires the addition of gradient pulses to the pulse sequence which cancel, or null, the effect on the NMR signal phase caused by spins moving in the gradients employed for position encoding. Such a solution is disclosed, for example, in U.S. Pat. No. 4,731,583 entitled “Method For Reduction of NMR Image Artifacts Due To Flowing Nuclei By Gradient Moment Nulling”. 
     A successful method for correcting MR images for motion artifacts employs navigator signals acquired during the scan. As described in U.S. Pat. No. 4,937,526, such navigator signals are acquired periodically during the scan, and the information in these signals may be used to correct the image data for patient motion. This method is also disclosed in U.S. Pat. Nos. 5,581,184; 6,118,273; 5,800,354; 6,076,006 and 5,539,312. Unfortunately, acquisition of the navigator signals increases the scan time and can disturb the steady-state magnetization established during fast scanning. 
     An automatic correction method has been proposed in International PCT application WO9801828A1 in which the entropy of the reconstructed image is examined as a focus criterion by which to iteratively adjust motion estimate. The clinical application of the autocorrection method has been made possible by the discovery of improved metrics used to evaluate the quality of the image during each iteration. A number of such metrics are disclosed in co-pending PCT patent application No. PCT/US99/08123 filed on Apr. 14, 1999 and entitled “Autocorrection of MR Images for Motion Artifacts.” One disadvantage of this correction technique is that it requires many repeated iterations and is computationally intensive. 
     The need for better motion correction methods has increased with the use of MRI systems having higher field strengths (e.g. 3 Tesla). Such systems enable the acquisition of images with smaller voxel sizes (e.g. smaller than 1 mm×1 mm×1 mm), with resulting higher image resolutions. However, such high resolution images are also easier to corrupt with motion artifacts. 
     SUMMARY OF THE INVENTION 
     The present invention is a method for correcting acquired NMR image data for subject motion, and particularly, a correction method that does not require the addition of navigator pulse sequences to the scan. The method includes acquiring a k-space NMR data set by repeating an NMR image pulse sequence to sample k-space as a series of k-space segments; reconstructing a series of tracking images from the corresponding series of k-space segments; calculating the movement of a point object in the subject during the acquisition of each k-space segment by detecting the movement of a corresponding point spread object in each tracking image; calculating a corrective phase shift for each k-space segment based on the detected movement during the acquisition; and reconstructing an image using the corrected k-space segments. The k-space segments are selected such that the location of the point object can be determined from the point spread object in their reconstructed tracking images, and in the preferred embodiment, the elliptical centric k-space view order is used for a three-dimensional acquisition, and each k-space segment is an annular region (i.e. ring) that surrounds the origin of k-space. 
     A general object of the invention is to correct acquired NMR image data for bulk (i.e. rigid body) subject motion without the need for additional, interleaved navigator pulse sequences and without the need for computationally intensive, iterative auto-correction post processing. By judiciously selecting the view order, tracking images can be reconstructed which enable segments of acquired k-space data to be separately motion corrected. A point object that produces an intense NMR signal is located in the field of view and produces a point spread object in each reconstructed tracking image. The point object may be a fiducial which is attached to or inserted in the subject being imaged, or it may be part of the subject, such as a blood vessel imaged in cross-section, containing contrast agent. Annular regions (rings) of k-space are used in the preferred embodiment since, with the elliptical centric view order, the views acquired during any time interval At will lie within such a ring. Thus, the entire acquisition can be sub-divided into a series of time intervals, corresponding to a series of concentric rings. 
     A more specific object of the invention is to motion correct a contrast-enhanced magnetic resonance angiogram (CE-MRA). The flow of contrast agent through the vasculature of interest provides intense point objects which move with the patient . A centric view order, preferably elliptical centric is used to acquire the NMR image data when the contrast agent bolus arrives in the vasculature of interest. The main advantage of using the elliptical centric view order for CE-MRA is that it can provide a high degree of venous suppression. The view order is determined by sorting the views according to their distance to the center (origin) of the k y -k z  plane. The sampling trajectory starts at the center of k-space and spirals outward. Each “ring” of k-space sampled in this spiral trajectory may form a k-space segment which can be separately motion corrected using the present invention. 
     The foregoing and other objects and advantages of the invention will appear from the following description. In the description, reference is made to the accompanying drawings which form a part hereof, and in which there is shown by way of illustration a preferred embodiment of the invention. Such embodiment does not necessarily represent the full scope of the invention, however, and reference is made therefore to the claims herein for interpreting the scope of the invention. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is a block diagram of an MRI system which employs the present invention; 
     FIG. 2 is a graphic representation of a 3DFT pulse sequence performed by the NMR system of FIG. 1 to practice a preferred embodiment of the invention; 
     FIG. 3 is a graphic representation of a conventional Cartesian k-space sampling trajectory using the pulse sequence of FIG. 2; 
     FIG. 4 is a graphic representation of a centric k-space sampling trajectory using the pulse sequence of FIG. 2; 
     FIG. 5 is a flow chart illustrating the steps employed in the preferred embodiment to motion correct an image acquired with the centric view order in FIG. 4; 
     FIG. 6 is a graphic representation of an annular region of k-space sampled during a centric acquisition and used to reconstruct a tracking image; 
     FIG. 7 is a graphic representation of point spread functions for a central disk of k-space and two subsequently-acquired rings of k-space; 
     FIG. 8 is a graphic representation showing k-space segments used in one preferred embodiment of the invention; and 
     FIGS. 9 b ,  9   d  and  9   f  are exemplary target images reconstructed from segments of k-space in FIGS. 9 a ,  9   c  and  9   e , respectively. 
    
    
     GENERAL DESCRIPTION OF THE INVENTION 
     The motion-correction method of the present invention exploits the property that during any time interval At during an elliptical centric or spiral acquisition, an annular region (ring) of raw data in the k y -k z , plane is acquired. Separate, independent tracking images are reconstructed from each ring of k-space data. Reconstructing an image solely from a ring of data introduces a point-spread function (PSF). The key property of this PSF is that any point object in the field of view being imaged displays a characteristic bull&#39;s-eye pattern in the tracking image. The bull&#39;s-eye pattern is centered at the true location of the point object. By registering the centers of the bull&#39;s-eye patterns in each of the tracking images, a motion correction can be performed in the full, NMR data set. 
     The point objects can either be externally placed fiducial markers which are fixed to the patient, or they can be small bright structures in the patient, such as contrast-enhanced vessels. The dimensions of fiducial markers can be controlled, and minimizing the size of the “point” object makes it easier to detect in tracking images reconstructed from peripheral k-space rings, (i.e., data acquired later in the scan), since object size and extent in k-space are reciprocally related. Also, the constant signal produced by a fiducial marker can make it easier to track than a vessel, which typically has a time-dependent contrast enhancement. 
     In the preferred embodiment of the invention an elliptical concentric view order such as that disclosed in U.S. Pat. No. 5,912,557 is used to sample 3D k-space as a series of annular rings of increasing diameter. As the following discussion will show, each tracking image reconstructed from such an annular ring of k-space data produces a very imperfect image, but an image which includes image objects that correspond to each strong point object in the subject. Most importantly, even though the image objects do not look like the point objects, due to the grossly undersampled k-space, one can determine the location of the point objects in the tracking image. In the case of the elliptical concentric view order, the image objects are rings which are concentric about the location of their corresponding point object positions. 
     This is illustrated in FIG. 9, where FIG. 9 a  illustrates the first k-space data segment acquired during the first 5 seconds of a scan and FIGS. 9 c  and  9   e  illustrate annular rings of k-space acquired in the following two 5 second intervals. FIG. 9 b  is an image reconstructed from the central disc of k-space data in FIG. 9 a  and because of the gross undersampling of k-space, objects are very blurred and aliased. Nevertheless, the location of point objects (in this case contrast enhanced arteries) can be determined. FIG. 9 d  illustrates the image reconstructed from the annular k-space segment of FIG. 9 c , and FIG. 9 f  illustrates the image reconstructed from the annular k-space segment in FIG. 9 e . In both cases, even though the images suffer in quality, the location of the point sources can be found. 
     Consider an annular region (i.e., a ring) of k-space with inner radius k 0 , and outer radius k 0 +Δk as shown in FIG.  6 . The effect of selectively reconstructing a tracking image from this ring of k-space data can be determined by calculating the PSF. The PSF is calculated by performing the two-dimensional 
     Fourier transform of the ring. Defining polar coordinates (k, θ) in k-space, where                k   =         k   y   2     +     k   z   2           ,                θ   =     arctan                     k   y       k   z                   (   1.1   )                 PSF        (     r   →     )       =         ∫   ring                   -   2                   π                           k   →     ·     r   →                                2        k         =       ∫     k   =     k   0           k   0     +     Δ                 k                ∫     θ   =   0       2                 π                     -   2                   π                                  kr                 cos                 θ                          θ                   k                      k                     (   1.2   )                                
     The factor of 2π in the exponential in Eq. (1.2) is required since the gyromagnetic ratio will be measured in units of linear frequency per field, e.g. 4257 Hz/Gauss. The integrals in Eq. (1.2) can be evaluated analytically to yield:                PSF        (   r   )       =           (       k   0     +     Δ                 k       )            J   1          (     2                 π                   r        (       k   0     +     Δ                 k       )         )         -       k   0            J   1          (     2                 π                   rk   0       )           r             (   1.3   )                                
     where J 1 (x) is the first order Bessel function. The full-width at half maximum (FWHM) can be calculated numerically. An approximate analytical expression for the FWHM can be determined from the Taylor expression:                  J   1          (   x   )       =       2   x     -       x   3     16     +   ⋯             (   1.4   )                                
     which yields              FWHM   ≈         1   π              8        (       k   0     +       Δ                 k     2       )           Δ                   k   3       +     4                 Δ                   k   2          k   0       +     6                 Δ                   kk   0   2       +     4        k   0   3               +   ⋯             (   1.5   )                                
     The value of the FWHM given by Eq. (1.5) provides a good estimate of the FWHM, typically underestimating its true value by about 10%. 
     To convert the parameters k 0  and Δk which appear in equations. (1.3) and (1.5) to a specified time interval during the acquisition, first the rate of k-space area traversal R is calculated:              R   =         Δ                   k   y        Δ                   k   z       TR     =     1       FOV   z          FOV   y        TR                 (   1.6   )                                
     where TR is the pulse sequence repetition time, and FOV y  and FOV z  are the fields of view in the y and z directions, respectively. If a ring of k-space is acquired during the time interval (t 0 →t 0 +Δt), then the k-space parameters k 0  and Δk which appear in the point-spread function are given by:                  k   0     =         Rt   0     π         ,     
            Δ                 k     =       -     k   0       +         k   0   2     +       R                 Δ                 t     π                     (   1.7   )                                
     FIG. 7 shows examples of the point-spread function plotted for a central disk, and two progressively larger rings of k-space, all with Δt=5 seconds: t=0-5 seconds, t=20-25 seconds, and t=30-35 seconds. The value of R is set to 1.43. cm −2 /sec. 
     The point-spread function in equation (1.3) has several important properties that will be used to analyze the motion correction method. First note that the point-spread is solely a function of the magnitude of the vector r: 
     
       
           PSF ( {right arrow over (r)} )= PSF (|| {right arrow over (r)} ||)= PSF ( r ),  r={square root over (y 2   +z   2 )}   (1.8) 
       
     
     As a consequence of equation (1.8) an image of a point source in the y-z plane will be spread out concentrically about its true location in the reconstructed tracking image. The zero-crossings of the point-spread function generate a bull&#39;s-eye pattern in the tracking image depiction of a point source. These two properties are important because it makes it feasible to track the true location of the point source in the tracking image reconstructed from a ring of k-space, simply by finding the center of the resultant bull&#39;s-eye pattern. 
     This property of the PSF is proved by calculating its convolution with a point source represented by a Dirac delta function in Cartesian or polar coordinates. The delta function provides a convenient mathematically approximation of a small, disk shaped object like a marker or an enhanced vessel in cross-section. If the point source is located at (y 0 , z 0 )=(r 0  cos φ0 , r 0 sin φ0 ):                I        (     y   ,   z     )       =         ∫       r   ′     =   0     ∞            ∫     ϕ   =   0       2                 π            δ                   (       z   ′     -     z   0       )                   δ                   (       y   ′     -     y   0       )          PSF        (       r   →                -       r   →     ′       )            r   ′             r                        ϕ           =         ∫       r   ′     =   0     ∞            ∫     ϕ   =   0       2                 π                           δ                   (       r   ′     -     r   0       )         r   ′          δ                   (       ϕ   ′     -     ϕ   0       )          PSF        (            r   →                -       r   →     ′            )            r   ′             r             ϕ           =     PSF        (            r   →                -       r   →     0            )                   (   1.9   )                                
     Another important property of the point-spread function in equation (1.3) is its functional dependence for large values of r. Since                    J   1          (   x   )       =         2     π                 x                       cos                   (     x   -       3                 π     4       )         ,       where                 x                 is        1     ,           (   1.10   )                                
     the radial dependence of the envelope of the point-spread function is                  envelope        [     PSF        (   r   )       ]       =     1     r     3   2           ,     r        1     2                 π                   (       k   0     +     Δ                 k       )                   (   1.11   )                                
     According to equations (1.9) and (1.11) the image intensity from a point source reconstructed from a ring of k-space will falloff as the inverse three-half power of the distance from the point source. In practice the larger the exponent the better, since more rapid falloff of the ringing leads to less interference between the point source and the ringing from other objects in the tracking image. 
     The functional dependence of the falloff in equation (1.11) highlights an advantage of using rings of k-space rather than of strips of k-space that would correspond to a conventional Cartesian acquisition. The equivalently normal PSF for a Cartesian acquisition is                PSF        (   y   )       =       N   z        Δ                   k   z        Δ                   k   y                       sin                   (     π                 Δ                   k   y        y     )         π                 Δ                   k   y        y                 (   1.12   )                                
     where N z  is the number of z-points (i.e., slices or sections) in the acquisition. Comparing equations (1.3) and (1.12), the envelope of the PSF from a strip of k-space falls off only as the inverse first power of y, i.e. more slowly than for the rings. This means that the pint source used for tracking can be more readily discerned from surrounding structures when the elliptical centric view order is used. 
     DESCRIPTION OF THE PREFERRED EMBODIMENT 
     Referring first to FIG. 1, there is shown the major components of a preferred MRI system which incorporates the present invention. The operation of the system is controlled from an operator console  100  which includes a keyboard and control panel  102  and a display  104 . The console  100  communicates through a link  116  with a separate computer system  107  that enables an operator to control the production and display of images on the screen  104 . The computer system  107  includes a number of modules which communicate with each other through a backplane. These include an image processor module  106 , a CPU module  108  and a memory module  113 , known in the art as a frame buffer for storing image data arrays. The computer system  107  is linked to a disk storage  111  and a tape drive  112  for storage of image data and programs, and it communicates with a separate system control  122  through a high speed serial link  115 . 
     The system control  122  includes a set of modules connected together by a backplane. These include a CPU module  119  and a pulse generator module  121  which connects to the operator console  100  through a serial link  125 . It is through this link  125  that the system control  122  receives commands from the operator which indicate the scan sequence that is to be performed. The pulse generator module  121  operates the system components to carry out the desired scan sequence. It produces data which indicates the timing, strength and shape of the RF pulses which are to be produced, and the timing of and length of the data acquisition window. The pulse generator module  121  connects to a set of gradient amplifiers  127 , to indicate the timing and shape of the gradient pulses to be produced during the scan. The pulse generator module  121  also receives patient data from a physiological acquisition controller  129  that receives signals from a number of different sensors connected to the patient, such as ECG signals from electrodes or respiratory signals from a bellows. And finally, the pulse generator module  121  connects to a scan room interface circuit  133  which receives signals from various sensors associated with the condition of the patient and the magnet system. It is also through the scan room interface circuit  133  that a patient positioning system  134  receives commands to move the patient to the desired position for the scan. 
     The gradient waveforms produced by the pulse generator module  121  are applied to a gradient amplifier system  127  comprised of G x , G y  and G z  amplifiers. Each gradient amplifier excites a corresponding gradient coil in an assembly generally designated  139  to produce the magnetic field gradients used for position encoding acquired signals. The gradient coil assembly  139  forms part of a magnet assembly  141  which includes a polarizing magnet  140  and a whole-body RF coil  152  a transceiver module  150  in the system control  122  produces pulses which are amplified by an RF amplifier  151  and coupled to the RF coil  152  by a transmit/receive switch  154 . The resulting signals radiated by the excited nuclei in the patient may be sensed by the same RF coil  152  and coupled through the transmit/receive switch  154  to a preamplifier  153 . The amplified NMR signals are demodulated, filtered, and digitized in the receiver section of the transceiver  150 . The transmit/receive switch  154  is controlled by a signal from the pulse generator module  121  to electrically connect the RF amplifier  151  to the coil  152  during the transmit mode and to connect the preamplifier  153  during the receive mode. The transmit/receive switch  154  also enables a separate RF coil (for example, a head coil or surface coil) to be used in either the transmit or receive mode. 
     The NMR signals picked up by the RF coil  152  are digitized by the transceiver module  150  and transferred to a memory module  160  in the system control  122 . When the scan is completed and an entire array of data has been acquired in the memory module  160 , an array processor  161  operates to Fourier transform the data into an array of image data. This image data is conveyed through the serial link  115  to the computer system  107  where it is stored in the disk memory  111 . In response to commands received from the operator console  100 , this image data may be archived on the tape drive  112 , or it may be further processed by the image processor  106  and conveyed to the operator console  100  and presented on the display  104 . 
     The NMR system of FIG. 1 performs a series of pulse sequences to collect sufficient NMR data to reconstruct an image. Referring particularly to FIG. 2, an exemplary pulse sequence for conducting a 3 DFT NMR scan is shown. The pulse sequence commences by the selective excitation of the entire region of interest with an RF excitation pulse  25  in the presence of a slab select G z  gradient pulse  26 . The frequency content of the excitation pulse  25  and the amplitude of the slab select G z  pulse  26  are selected to produce transverse magnetization in the region which is the subject of the 3D scan. A negative G z  pulse  27  is then produced to rephase the spins in preparation for the phase encoding and readout. 
     Phase encoding is performed along two axes: the z-axis and the y-axis. The z-axis encoding is accomplished by applying a G z  phase encoding pulse  28  and the y-axis encoding is accomplished by applying a G y  phase encoding pulse  29 . As is well-known to those skilled in the art, the magnitude of the phase encoding pulses  28  and  29  are stepped through a series of positive and negative values during the scan, but each is set to one value during each pulse sequence. As will be described in detail below, it is the order in which these phase encoding pulses  28  and  29  are stepped through their set of values which is of importance to practicing the present invention. As is well-known in the art, the magnitude of a phase encoding gradient pulse is determined by the integral of its amplitude over its duration (i.e. its area). In most pulse sequences the duration is kept constant and the phase encoding pulse magnitude is stepped through its values by changing its amplitude. 
     After phase encoding the transverse magnetization, the NMR signal  30  is read-out in the presence of a G x  read-out gradient  31 . This read-out is preceded by a negative G x  gradient pulse  32  to produce the gradient refocused NMR echo signal  30  in the usual fashion. It can be appreciated by those skilled in the art that either a full echo signal  30  as shown in FIG. 2 may be acquired, or a partial echo signal may be acquired to reduce the pulse sequence TR. In the latter case, the missing samples are estimated using a homodyne reconstruction method or the like. The 3DFT pulse sequence is then concluded by the application of a large G z  spoiler gradient pulse  33  and a G y  rewinder gradient pulse  34  to prepare the magnetization for the next pulse sequence which follows immediately. As is known to those skilled in the art, the spoiler pulse  33  dephases transverse magnetization and the rewinder pulse  34  refocuses transverse magnetization along the y-axis in preparation for the next pulse sequence. The rewinder pulse  34  is equal in magnitude, but opposite in polarity with the G y  phase encoding pulse  29 . 
     The acquisition of data in 3DFT scanning can be considered sampling of a three-dimensional “k-space”. Two of the dimensions, k y  and k z , are sampled by applying different phase encoding gradients G y  and G z  during each pulse sequence of the scan, and each acquired NMR signal contains 256 samples along a line in the k x  direction. The k y  and k z  phase encodings are defined as the gyromagnetic ratio (γ) times the area under the G y  and G z  gradient waveforms applied after the RF excitation pulse  25  and before the readout of echo signal  30 . The pulse sequence is repeated for as many repetitions as are necessary to sample all desired k y  and k z  values. For example, k y  may assume 128 different values and k z  may have 64 values. In this case the number of repetitions of the pulse sequence of FIG. 2 would be 128×64 or 8192. 
     With conventional 3DFT scanning the desired values of k y  and k z  are sampled with two nested loops. For example, the inner loop increments k y  through its 128 values and after all such samples have been made for a given value of k z , the outer loop increments k z . This process continues until all 128 values of k y  have been sampled at each of the 64 values of k z . A diagram of this prior process is shown in FIG. 3, where the vertical axis indicates the sample along k z  produced by the magnitude of the G z  phase encoding pulse, and the horizontal axis indicates the sample along k y  produced by the magnitude of the G y  phase encoding pulse. The arrow shows the temporal trajectory through (k y , k z ) space of the sampling during the scan. The sampling starts in the lower left corner with negative k y  and k z  samples, and ends at the upper right comer with positive k y  and k z  samples. In this conventional scan, therefore, (k y , k z ) space is scanned in much the same way as a television screen is scanned. 
     Centric view ordering is based on the realization that for most objects the bulk of the signal power is contained in the samples taken near the origin of (k y , k z ) space, and it is these samples which contribute most significantly to the appearance of the reconstructed image. This results from the fact that the NMR signals acquired during the scan are Fourier transformed along the k x , k y  and k z  directions to produce intensity values for an image in real (x,y,z) space. It is the nature of this transformation that the samples near the origin (k y =0, k z =0) contribute a disproportionate share to the signal power of the reconstructed image. Accordingly, it is a basic idea of centric view ordering to sample the (k y , k z ) points which contain the most signal power in as short a time and as close to the beginning of the scan as possible. This can be done by modifying the trajectory with which the (k y , k z ) space is sampled. Specifically, a square spiral (k y , k z ) trajectory, such as that shown in FIG. 4 may be used. The scan starts at or near the origin of (k y , k z ) space and progressively works its way outward in a spiral fashion. In FIG. 4, a 12×16 array of k-space samples are acquired, but this is only illustrative, since in practice many more samples are usually acquired to cover the field of view with adequate resolution. The manner in which the values of the G y  and G z  phase encoding gradients may be stepped to accomplish such a spiral trajectory scan is described in U.S. Pat. No. 5,122,747. 
     A second, and more preferred centric view order is the so-called elliptical centric view order described in U.S. Pat. No. 5,912,557. This is a view order based on the distance of the sample point from the origin of k-space. First, we assume that the field of view along y is FOV and that the field of view along z is some fraction of this, FOV/N. The k y  values which are to be sampled in the phase encoding process are +1/(2 FOV), ±3/(2 FOV), ±5/(2 FOV), . . . , ±m/(2 FOV), where 2m is the total number of k y  phase encodes desired. This assumes that the k y  origin is not sampled and the smallest nonzero k y  spatial frequencies sampled are at ±1/(2 FOV). Similarly the k z  values which are to be sampled are ±N/(2 FOV), ±N3/(2 FOV), ±N5/(2 FOV), . . . , ±N/(2 FOV), where 2 n is the total number of k z  phase encodes. All phase encodings to be sampled can then be characterized by the index (i,j), where −m≦i≦m and −n≦j≦+n and i and j assume nonzero integer values. The distance from the k-space origin to the sample point is characterized by the indices (ij) and is equal to (1/FOV) sqrt {[2¦i¦−1] 2 +N 2 [2¦j¦−1] 2 } where “sqrt” denotes the operation of taking the square root and ␣*␣ denotes absolute value of *. Once the distance is determined for all k-space sample points (i,j), the points are ranked on the basis of their respective distances with the point having the smallest distance being ranked first. In this way, the ordered list of desired phase encodings, or views, is determined. 
     The advantage of the elliptical centric view order for the motion correction described here is that k-space is sampled as a series of circular-shaped annular rings. The point spread function for circular rings is superior to the point spread function of other shaped segments of k-space in that a point source is more easily located in the reconstructed target image. 
     While most MRA clinical applications benefit from the centric view order in which the center of k-space is acquired first, there are applications in which the same k-space sampling trajectory is employed in reverse order. This “reverse centric” view order also follows a spiral-like sampling trajectory which circles the center of k-space, but in this case sampling spirals inward toward the center of k-space. 
     The characteristic of both the centric and reverse-centric view orders that is important to the present invention is that the acquisition can be divided into a series of k-space sampling segments which are encircling, ring-like segments of k-space. Each k-space sampling segment is comprised of k-space samples from one or more complete revolutions around the center of k-space. These ring-like k-space segments have point spread functions which enable the location of point objects to be determined in images reconstructed from the k-space samples. 
     Referring particularly to FIG. 5, the preferred embodiment of the invention is employed in a CE-MRA acquisition of the carotid and vertebral arteries of a patient. As indicated at process block  300 , a contrast agent is injected into the patient, and when the contrast arrives in the arteries of interest, a 3D image is acquired using the above-described pulse sequence and a centric view order as indicated at process block  304 . Contrast arrival can be detected as indicated at process block  302  in a number of different ways. In the preferred embodiment low resolution 2D images are acquired at a high frame rate and the physician initiates the image acquisition when enhanced contrast is seen in the arteries. Automatic methods can also be used to initiate the scan such as the method disclosed in U.S. Pat. No. 6,167,293. A three-dimensional array of k-space data is acquired during the scan and stored in memory. 
     As indicated in process block  306 , the first step in the image reconstruction process is to perform a fast Fourier transformation of the 3D k-space data set along the readout gradient axis. In the preferred embodiment the readout axis is the x-axis. It can be appreciated by those skilled in the art that if a partial echo acquisition is employed to acquire the k-space image data set this step is more complicated than a one-dimensional FFT. Instead, the missing k-space data are estimated using a technique such as that described by Douglas Noll, et al., “Homodyne Detection In Magnetic Resonance Imaging,”  IEEE Trans. On Medical Imaging , vol. 10, Jun. 2, 1991. Also, in that case, the homodyne reconstruction in the readout direction would be performed after the Fourier transformations along the phase encoding axes. In either case, the result is a three-dimensional array of data composed of a set of k y , k z  planes disposed along the x axis. One of these planes is selected, as indicated at process block  307 , as being representative of the bulk motion of the entire region of interest. Phase corrections will be calculated for the selected k y , k z  plane and these phase corrections are applied to corresponding k y , k z  elements in all the x-axis planes. In alternative embodiments of the invention the phase corrections can be separately calculated for each k y , k z  plane in the data set, or the data set may be divided into separately corrected groups, or slabs, of k y , k z  planes. 
     As shown in FIG. 8 the k-space defined by the k y  and k z  axes in the selected plane along the z-axis is divided into segments. These segments include a central disc of k-space samples bounded by dashed line  308  and a plurality of annular ring-shaped segments identified by dashed lines  310 ,  312  and  314 . FIG. 8 does not show the number of samples in each segment nor the total number of segments, but rather is illustrative of how k-space may be divided into segments. In the preferred embodiment the number of (k y -k z ) sample points (i.e. views) is substantially the same in each segment such that the time interval Δt required to acquire each segment is substantially the same. This time interval should be kept as short as possible and in the preferred embodiment is less then 5 seconds. 
     As indicated at process block  316 , the next step is to reconstruct an image from the segment  308 . This is done by transferring the k-space samples in segment  308  to an array in which all other sample points in the 3D array are set to zero. A 2DFT is then performed along the y and z axes and the magnitude of the transformed complex values is calculated. The resulting image forms a reference image against which the location of the subject is compared in the subsequently acquired ring segments. A loop is then entered at  318  in which the first step indicated at process block  320  is to reconstruct an image using k-space samples in the next surrounding segment. This is done with a 2DFT as described above, and a magnitude image is calculated from the result. 
     The next step as indicated at process block  322  is to determine the extent to which the tracking image is mis-registered with the reference image. This is accomplished by finding the peak of the cross-correlation function of the two images. More specifically, if the reference image R(y,z) is:          R        (     y   ,   z     )       =              ∫   ∫       central        
        segment                 (       k   y     ,     k   z       )                 2                 π                                    (         k   y        y     +       k   z        z       )                 y             z                                   
     where: 
     d(k y , k z ) is the complex k-space data, 
     i={square root over (−1)}, and 
     |*|=the modulus of the complex quantity *. 
     Similarly, the tracking image T(y,z) is given by:          T        (     y   ,   z     )       =                ∫   ∫       ring        
        segment                 (       k   y     ,     k   z       )                 2                 π                                    (         k   y        y     +       k   z        z       )                 y             z            .                            
     The degree of motion between the reference and tracking images is calculated by finding the peak of the cross-correlation function:        C   =       (       Δ                 y     ,     Δ                 z       )     =         ∫   ∫         y        -        z                 image          
        plane            R        (     y   ,   z     )            T        (       y   +     Δ                 y       ,     z   +     Δ                 z         )                          y               z     .                                
     The cross-correlation function C(Δy, Δz) is expected to peak at values (Δy, Δz) corresponding to motion that occurred between the acquisition of the data for the reference and tracking images. 
     Once the values of Δy and Δz are determined, the complex k-space data in the ring segment that produced T(y,z) is corrected with the shift theorem as indicated at process block  324 :                 (       k   y     ,     k   z       )     corrected       =            (       k   y     ,     k   z       )                     -   2                   π                                    (         k   y          Δ   y       +       k   z          Δ   z         )         .                              
     It should be noted that corrective values Δy and Δz can be either positive or negative, and in general, they will be different for each ring segment. 
     As indicated by decision block  326 , the correction procedure is repeated for each ring k-space segment until the last k-space segment is corrected. A 3D k-space data set is thus produced in which all the ring k-space segments are phase corrected to offset any motion relative to the reference location. A 3D image is then produced at process block  328  by performing a fast Fourier transformation along the two phase encoding axes (k y  and k z  in the example) of the corrected 3D k-space data set. A 3D magnitude image is produced from the resulting complex numbers at process block  330  and a 2D projection image is produced at process block  323 . Many different projection methods are known to those skilled in the art, but a maximum intensity pixel (MIP) method is preferred. 
     While the preferred embodiment of the invention is a CE-MRA scan using an elliptical centric view order, other clinical applications are also possible. The contrast enhanced blood vessels in CE-MRA provide good point sources that enable accurate motion corrections to be made. However, fiducials in the form of small (i.e. point source) phantoms that provide very intense NMR signals may be attached to the subject to track with subject motion. This enables the present invention to be used in other clinical applications where strong point source signals are not naturally present in the subject being imaged. Also, the use of fiducials which produce strong point source signals may enable view orders to be used which have point spread functions that are less desirable than the elliptical centric view order.