Abstract:
A method for correcting the background phase in magnetic resonance phase contrast flow images includes providing a time series of velocity encoded magnetic resonance images of a patient, where the time series of velocity encoded images comprises for each time point a phase contrast image where a pixel intensity is proportional to a flow velocity, measuring a change of intensity for each pixel over the time series of phase contrast images, identifying pixels with a low measure of temporal change as stationary pixels, and calculating a correction field for the stationary pixels, where the correction field represents a background phase to be subtracted from the phase contrast image.

Description:
CROSS REFERENCE TO RELATED UNITED STATES APPLICATIONS 
     This application claims priority from “Background Phase Correction For Phase Contrast Flow Images”, U.S. Provisional Application No. 61/387,187 of Guehring, et al., filed Sep. 28, 2010, the contents of which are herein incorporated by reference in their entirety. 
    
    
     TECHNICAL FIELD 
     This disclosure is directed to methods for correcting background phase in phase contrast magnetic resonance images. 
     DISCUSSION OF THE RELATED ART 
     Phase contrast cardiovascular magnetic resonance imaging (MRI) sequences are used to measure blood flow velocity in vessels or at the heart valves. Velocity encoded MRI is based on detection of changes in the phase of proton precession. These changes are proportional to the velocity of the movement of those protons through a magnetic field with a known gradient. When using velocity encoded MRI, the result of the MRI scan is usually two sets of images for each time point in the cardiac cycle: an anatomical image and a phase contrast image where the signal intensity in each pixel is directly proportional to the velocity of a fluid flowing in the image. Note, however, that velocity encoded MRI can be performed using only the phase contrast images. The average velocity in a vessel is quantified by measuring the average signal intensity of the pixels in the cross section of the vessel, and then multiplying by a known constant. The flow is calculated by multiplying the mean velocity by the cross-sectional area of the vessel. This flow data can be used to graph flow versus time. The area under the flow versus time curve for one cardiac cycle is the stroke volume. MRI is typically used to quantify the flow over one cardiac cycle as the average of several heart beats, but it is also possible to quantify the stroke volume in real time on a beat-for-beat basis. This flow can be used to calculate cardiac output, shunt flow, and aortic or pulmonary regurgitation. 
     However, the velocity-encoded phase images should be corrected for velocity offset errors that occur during image acquisition. These errors are due to non compensated eddy-current-induced fields and concomitant gradient terms present in MR systems. As volume flow assessment in large vessels in determined by integration over time and space, relatively small velocity offsets can lead to significant errors in mean volume flow, stroke volume, and cardiac output. 
     In some cardiovascular MRI machines, a single reference region may be defined by the user to delineate stationary tissue and measure the background phase. It is assumed that the background phase in the vessel region will be identical. Therefore, the background phase of the reference region can be subtracted from the flow curve calculated from the vessel region. The challenge with this solution is that it is now known that the background is not uniform throughout the region and should not be estimated from a single reference region. 
     A more recent approach includes fitting a surface through the stationary pixels in the image. The stationary pixels are determined by calculating the variance over time of the gray levels at each pixels and choosing a predetermined percentage λ of pixels with lowest variance. The surface is then approximated by a two dimensional Taylor expansion of order k. Published results for this approach suggest that the best results were obtained for λ=25% and k=1. However, the stationary pixels are chosen arbitrarily, independent of the sequence. 
     SUMMARY OF THE INVENTION 
     Exemplary embodiments of the invention as described herein generally include methods and systems for adaptively determining stationary pixels based on the sequence, and fitting a surface through the stationary pixels. 
     According to an aspect of the invention, there is provided a method for correcting the background phase in magnetic resonance phase contrast flow images, including providing a time series of velocity encoded magnetic resonance images of a patient, where the time series of velocity encoded images comprises for each time point a phase contrast image where a pixel intensity is proportional to a flow velocity; measuring, for each pixel, a change of intensity over the time series of phase contrast images; identifying pixels with a low measure of temporal change as stationary pixels; and calculating a correction field for the stationary pixels, where the correction field represents a background phase to be subtracted from the phase contrast image. 
     According to a further aspect of the invention, the correction field is a planar surface determined by a least squares fit to the stationary points. 
     According to a further aspect of the invention, the correction field is a free-form surface. 
     According to a further aspect of the invention, the free-form surface is determined by calculating a 2D tensor of uniform 1D cubic B-splines, 
                   T   local     ⁡     (     x   ,   y     )       =       ∑     m   =   0     3     ⁢       ∑     n   =   0     3     ⁢         B   m     ⁡     (   u   )       ⁢       B   n     ⁡     (   v   )       ⁢     φ       p   +   m     ,     q   +   n                 ,         
where (x, y) is the coordinate of pixel i, the image Ω={(x, y)|0≦x≦X, 0≦y≦Y} is represented by a sparse control point lattice φ p,q  having uniform grid spacing Δ x  and Δ y  in x and y directions, respectively, (x, y) is the coordinate of pixel i, p=└x/Δ x ┘−1, q=└y/Δ y ┘−1, u=x/Δ x −└x/Δ x ┘, and v=y/Δ y −└y/Δ y ┘, and B m  represents an m-th basis function of the B-spline.
 
     According to a further aspect of the invention, measuring a change of intensity for each pixel over the time series of phase contrast images comprises computing, for each pixel, an intensity variance over time series of phase contrast images. 
     According to a further aspect of the invention, identifying pixels with a low measure of temporal change comprises forming a histogram of intensity changes of all pixels, and determining from the histogram a threshold, where pixels whose intensity change is less than the threshold are considered as stationary pixel, and pixels whose intensity change is greater than the threshold are considered as moving pixels. 
     According to a further aspect of the invention, the background phase surface is calculated for an average phase contrast image that is an average over the time series of phase contrast images. 
     According to a further aspect of the invention, the method includes providing a time series of anatomical images, measuring, for each pixel, a change of intensity over the time series of anatomical images, normalizing the intensity change over the phase contrast images and the intensity change over the anatomical images, and selecting, for each pixel, the maximum of the intensity change over the phase contrast images and the intensity change over the anatomical images as the measure of temporal change for each pixel. 
     According to a further aspect of the invention, the method includes grouping neighboring stationary pixels together, and removing small and ill-shaped groups of pixels from the stationary pixels. 
     According to another aspect of the invention, there is provided a program storage device readable by a computer, tangibly embodying a program of instructions executable by the computer to perform the method steps for correcting the background phase in magnetic resonance phase contrast flow images. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIGS. 1(   a )-( f ) illustrate an algorithm to extract stationary pixels from a phase contrast image, according to an embodiment of the invention 
         FIG. 2  is a flowchart of an algorithm to extract stationary pixels from a phase contrast image, according to an embodiment of the invention. 
         FIG. 3  is a block diagram of an exemplary computer system for implementing a method for correcting background phase in phase contrast magnetic resonance (MR) images, according to an embodiment of the invention. 
     
    
    
     DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTS 
     Exemplary embodiments of the invention as described herein generally include systems and methods for correcting background phase in phase contrast magnetic resonance (MR) images. Accordingly, while the invention is susceptible to various modifications and alternative forms, specific embodiments thereof are shown by way of example in the drawings and will herein be described in detail. It should be understood, however, that there is no intent to limit the invention to the particular forms disclosed, but on the contrary, the invention is to cover all modifications, equivalents, and alternatives falling within the spirit and scope of the invention. 
     As used herein, the term “image” refers to multi-dimensional data composed of discrete image elements (e.g., pixels for 2-dimensional images and voxels for 3-dimensional images). The image may be, for example, a medical image of a subject collected by computer tomography, magnetic resonance imaging, ultrasound, or any other medical imaging system known to one of skill in the art. The image may also be provided from non-medical contexts, such as, for example, remote sensing systems, electron microscopy, etc. Although an image can be thought of as a function from R 3  to R or R 7 , the methods of the inventions are not limited to such images, and can be applied to images of any dimension, e.g., a 2-dimensional picture or a 3-dimensional volume. For a 2- or 3-dimensional image, the domain of the image is typically a 2- or 3-dimensional rectangular array, wherein each pixel or voxel can be addressed with reference to a set of 2 or 3 mutually orthogonal axes. The terms “digital” and “digitized” as used herein will refer to images or volumes, as appropriate, in a digital or digitized format acquired via a digital acquisition system or via conversion from an analog image. 
     An algorithm according to an embodiment of the invention to extract stationary pixels from a phase contrast image is illustrated by  FIGS. 1(   a )-( f ), and by the flowchart of  FIG. 2 .  FIG. 1(   a ) depicts an original magnitude image, with a liver  11 , lungs  12 , aortic valve  13 , and descending aorta  14 , and  FIG. 1(   b ) depicts an original phase contrast image, each selected from a sequence of velocity encoded MR images of a heart. The same reference numbers refer to the same item in  FIGS. 1(   a )-( f ). Referring to  FIG. 1(   b ), the salt and pepper regions  15  and  12  respectively correspond to outside air and the lungs, which contain air. Phase contrast MRI cannot measure velocity in air, thus air shows up as noise in these images. 
     Referring now to  FIG. 2 , given a sequence of velocity encoded MR images, an algorithm according to an embodiment of the invention begins at step  21  by determining the change of pixel intensity levels over time at each pixel in the sequence of velocity encoded MR images, illustrated in  FIG. 1(   c ). According to some embodiments of the invention, only a series of phase contrast is provided. According to other embodiments of the invention, both a time series of anatomical magnitude images and a time series of phase contrast images are provided. According to an embodiment of the invention, the change of pixel intensity levels is represented by the intensity variance of each pixel over the time series of images. According to some embodiments of the invention, only a series of phase contrast images is provided, in which case the variance is calculated over the series of phase contrast images. According to other embodiments of the invention, if both anatomical magnitude images and phase contrast images are provided, the variance is computed over time as the maximum of the variance over time of the magnitude images and the variance over time of the phase contrast images at each pixel, after a suitable normalization. 
     Next, at step  22 , the variances are analyzed to determine a threshold that separates the stationary pixels from the moving pixels. According to embodiments of the invention, a histogram is formed of the variance values. The histogram includes two peaks: one for the noisy background and air with a large variance and another peak for tissue with a smaller variance. The histogram, depicted in  FIG. 1(   d ), shows the number of pixels with a certain variance on the y-axis for all possible variances on the x-axis. One looks for the 2 peaks, separated by a valley, in that histogram. Only the first peak, closest to the origin with the largest variances, is of interest. At step  23 , a Gaussian mode is fit to first peak in the histogram, from which the mean d standard deviation (stddev) of the first peak in the histogram are computed. According to embodiments of the invention, the mean and standard deviation are used to determine a threshold. Exemplary, non-limiting thresholds definitions include: (1) Mean+n×StdDev, where n&gt;1; (2) the variance for which the value of the histogram is half the peak value; and (3) the larger of the two previous values. Exemplary, non-limiting values for n are 2 and 3. According to other embodiments of the invention, the threshold may be determined from a pre-selected percentage of the difference between a variance minimum and a variance maximum. Any pixel with a variance less than that threshold is considered a stationary pixel. All the pixels beyond the right tail of that mode are considered to be moving pixels.  FIG. 1(   e ) shows the mask formed by the stationary pixels. Notice how the descending aorta  14  and aortic valve  13  are not part of the stationary pixels. Then at step  24 , a correction field is fit through the stationary pixels of an average phase contrast image over time. This correction field, which represents the background phase, is shown in  FIG. 1(   f ). 
     In an additional optional step, neighboring stationary pixels can be grouped to together to improve spatial coherence. There may be holes, etc., in the groups, as long there exists a neighbor, such as 4- or 8-neighbor that is classified to be stationary as well. An exemplary, non-limiting algorithm for performing this grouping is a connected-component labeling algorithm. Once the groups are formed, small or very ill-shaped groups, such as very thin groups, can be removed, as they will likely not reliably contribute to the determination of the background region due to their potentially random nature and limited spatial resolution. 
     The simplest correction field to fit to the stationary pixels is a plane. This can be done using a least squares fit, which is equivalent to finding the eigenvalues and eigenvectors of the points in 3D and setting the normal of the plane to the eigenvector with smallest eigenvalue and a point on the plane to the center of mass of the points. 
     A free-form surface can also be fit to the stationary pixels to represent the correction field. A B-Spline freeform surface may be applied to approximate the bias field. In this representation, a dense 2D bias field is parameterized on a sparse control point lattice. The image is represented as Ω={(x, y)|0≦x≦X, 0≦y≦Y} and φ represents a grid of control points φ p,q  with the grid spacing being Δ x ×Δ y . This spacing between adjacent control points is uniform for each coordinate direction. A 2D tensor of uniform 1D cubic B-splines may be defined as follows: 
                 T   local     ⁡     (     x   ,   y     )       =       ∑     m   =   0     3     ⁢       ∑     n   =   0     3     ⁢         B   m     ⁡     (   u   )       ⁢       B   n     ⁡     (   v   )       ⁢     φ       p   +   m     ,     q   +   n                     
where (x, y) is the coordinate of pixel i, p=└x/Δ x ┘−1, q=└y/Δ y ┘−1, u=x/Δ x −└x/Δ x ┘, and v=y/Δ y −└y/Δ y ┘. B m  represents the m-th basis function of the B-spline. The basis functions of cubic B-splines have limited support. Therefore changing a control point in the grid affects only a 4×4 region around that control point.
 
     T local (x, y) is computed for all pixels in the image, including pixel positions belonging to moving and stationary. The input to compute T local (x, y) are the phase values for all stationary pixels. The output of T local (x, y) is a smoothing correction field for every pixel in the image. For a pixel location (x 1 , y 1 ) belonging to a stationary pixel, T local (x 1 , y 1 ) is a smoothing approximation of the background phase value at location (x 1 , y 1 ). For a pixel location (x 2 , y 2 ) belonging to a moving pixel, the phase value is initially unknown. Thus, T local (x 2 , y 2 ) is computed and treated as an estimation of background phase value of (x 2 , y 2 ). 
     So, the effect of computing T local (x, y) is to obtain a dense pixel-wise phase image for both stationary and moving pixels, given the phase values of the stationary pixels as the input. The properties of a bspline ensure that T local (x, y) is smooth across the image region. For a pixel (x, y) that is stationary, T local (x, y) will not in general equal the original phase value. This is reasonable, as the original phase values of stationary pixels include noise and phase changes coming from real tissue, while the background phase should be smooth and slowly changing. By computing a smoothing approximation of stationary phase values, the embedded components corresponding to background phase can be extracted. 
     A method according to an embodiment of the invention should produce good results for the following two reasons. First, the stationary pixels are not chosen arbitrarily using a fixed threshold, but rather using a threshold that depends on the images themselves. Second a free form surface fitting algorithm according to an embodiment of the invention has already been used successfully for cardiac perfusion surface coil correction in other settings. 
     It is to be understood that embodiments of the present invention can be implemented in various forms of hardware, software, firmware, special purpose processes, or a combination thereof. In one embodiment, the present invention can be implemented in software as an application program tangible embodied on a computer readable program storage device. The application program can be uploaded to, and executed by, a machine comprising any suitable architecture. 
       FIG. 3  is a block diagram of an exemplary computer system for implementing a method for correcting background phase in phase contrast magnetic resonance (MR) images, according to an embodiment of the invention. Referring now to  FIG. 3 , a computer system  31  for implementing the present invention can comprise, inter alia, a central processing unit (CPU)  32 , a memory  33  and an input/output (I/O) interface  34 . The computer system  31  is generally coupled through the I/O interface  34  to a display  35  and various input devices  36  such as a mouse and a keyboard. The support circuits can include circuits such as cache, power supplies, clock circuits, and a communication bus. The memory  33  can include random access memory (RAM), read only memory (ROM), disk drive, tape drive, etc., or a combinations thereof. The present invention can be implemented as a routine  37  that is stored in memory  33  and executed by the CPU  32  to process the signal from the signal source  38 . As such, the computer system  31  is a general purpose computer system that becomes a specific purpose computer system when executing the routine  37  of the present invention. 
     The computer system  31  also includes an operating system and micro instruction code. The various processes and functions described herein can either be part of the micro instruction code or part of the application program (or combination thereof) which is executed via the operating system. In addition, various other peripheral devices can be connected to the computer platform such as an additional data storage device and a printing device. 
     It is to be further understood that, because some of the constituent system components and method steps depicted in the accompanying figures can be implemented in software, the actual connections between the systems components (or the process steps) may differ depending upon the manner in which the present invention is programmed. Given the teachings of the present invention provided herein, one of ordinary skill in the related art will be able to contemplate these and similar implementations or configurations of the present invention. 
     While the present invention has been described in detail with reference to exemplary embodiments, those skilled in the art will appreciate that various modifications and substitutions can be made thereto without departing from the spirit and scope of the invention as set forth in the appended claims.