Abstract:
A method of regularization of x-ray phase contrast imaging (XPCi) system measurement data includes obtaining air scan data of the XPCi system prior to the presence of an object undergoing imaging, performing Fourier analysis of the air scan data, computing air coefficients from the result of the performing step, obtaining object scan data of an object undergoing imaging on the XPCi system, regularizing the object scan data, and calculating at least one of absorption image data, differential phase image data, and dark field image data by using object coefficients. A system configured to implement the method and a non-transitory computer-readable medium are disclosed.

Description:
BACKGROUND 
       [0001]    X-ray phase contrast imaging (XPCi) creates images of an object utilizing changes in the phase of an x-ray beam as it passes through the object. XPCi improves over problems with conventional x-ray imaging, where poor image contrast can arise from small attenuation differences. Its application can be found in medicine (e.g., breast tissue tumors) and industrial applications. Because many substances induce a phase shift in the x-rays as they pass through an object, the detection of features not discernable by conventional x-ray imaging is possible with XPCi. The phase shift can be measured using an interferometer to obtain high-contrast data for low-absorption objects. 
         [0002]      FIG. 1  depicts standard XPCi imaging process  100 . One implementation of XPCi is a grating based technique, where a series of images is acquired along with a grating stepping technique in order to sample the interference patterns (i.e., fringes) for each detector pixel (x,y). To correct for system non-homogeneities, sampling is performed for an air scan  105  to create a reference of system contributions. An object scan  110  is performed to collect data for imaging. 
         [0003]    The fringes can be well approximated as a sinusoidal curve. This curve is analyzed by performing Fourier analysis  115 . The Fourier analysis computes for the air scan and the object scan the first two (complex) Fourier coefficients a(x,y,0) and a(x,y,1) (0th and 1st order, respectively)  120 ,  125 . In principle the analysis is not limited to the first two coefficients but can be done up to order N. An exemplary, conventional formulation for process  100  can be expressed as 
         [0000]      min ρ ∥           H   ρ−b ( x,y ,•)∥ 2             a ( x,y ,•)=ρ  (EQ.1)
 
         [0004]    where (x, y) is pixel position; 
         [0005]               the Fourier transform up to order N; and 
         [0006]                H  is the adjoint of          . 
         [0007]    Since           is an orthogonal operator, the adjoint equals the inverse Fourier transform, and equation 1 can be solved directly as a(x,y,•)=         b(x, y,•) 
         [0008]    The three types of contrast provided by XPCi (i.e., absorption image att(x,y)  132 , differential phase image dpc(x,y)  134 , and dark field image dci(x,y)  138 ), can be computed by a complex division  130  of the Fourier coefficients of the object scan by the Fourier coefficients of the air scan, as follows 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       p 
                        
                       
                         ( 
                         
                           x 
                           , 
                           y 
                           , 
                           n 
                         
                         ) 
                       
                     
                     = 
                     
                       
                         
                           a 
                           obj 
                         
                          
                         
                           ( 
                           
                             x 
                             , 
                             y 
                             , 
                             n 
                           
                           ) 
                         
                       
                       
                         
                           a 
                           air 
                         
                          
                         
                           ( 
                           
                             x 
                             , 
                             y 
                             , 
                             n 
                           
                           ) 
                         
                       
                     
                   
                   , 
                   
                     n 
                     = 
                     
                       0 
                        
                       
                           
                       
                        
                       … 
                        
                       
                           
                       
                        
                       N 
                     
                   
                 
               
               
                 
                   ( 
                   
                     EQ 
                     . 
                     
                         
                     
                      
                     2 
                   
                   ) 
                 
               
             
             
               
                 
                   
                     att 
                      
                     
                       ( 
                       
                         x 
                         , 
                         y 
                       
                       ) 
                     
                   
                   = 
                   
                      
                     
                       p 
                        
                       
                         ( 
                         
                           x 
                           , 
                           y 
                           , 
                           0 
                         
                         ) 
                       
                     
                      
                   
                 
               
               
                 
                   ( 
                   
                     EQ 
                     . 
                     
                         
                     
                      
                     3 
                   
                   ) 
                 
               
             
             
               
                 
                   
                     dpc 
                      
                     
                       ( 
                       
                         x 
                         , 
                         y 
                       
                       ) 
                     
                   
                   = 
                   
                      
                     
                       ( 
                       
                         p 
                          
                         
                           ( 
                           
                             x 
                             , 
                             y 
                             , 
                             1 
                           
                           ) 
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   
                     EQ 
                     . 
                     
                         
                     
                      
                     4 
                   
                   ) 
                 
               
             
             
               
                 
                   
                     dci 
                      
                     
                       ( 
                       
                         x 
                         , 
                         y 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       
                         ∑ 
                         
                           
                             n 
                             = 
                             1 
                           
                           , 
                           3 
                           , 
                           5 
                           , 
                           
                               
                           
                            
                           … 
                            
                           
                               
                           
                           , 
                           
                               
                           
                            
                           N 
                         
                       
                        
                       
                           
                       
                        
                       
                          
                         
                           p 
                            
                           
                             ( 
                             
                               x 
                               , 
                               y 
                               , 
                               n 
                             
                             ) 
                           
                         
                          
                       
                     
                     
                        
                       
                         p 
                          
                         
                           ( 
                           
                             x 
                             , 
                             y 
                             , 
                             0 
                           
                           ) 
                         
                       
                        
                     
                   
                 
               
               
                 
                   ( 
                   
                     EQ 
                     . 
                     
                         
                     
                      
                     5 
                   
                   ) 
                 
               
             
           
         
       
     
         [0009]    Subsequently, phase image  150  can be computed by integrating  140  the differential phase image  123  on a line-by-line basis. However, integration step  140  is an optional part of the process. 
         [0010]    Advanced processing strategies of XPCi images (e.g., regularization/de-noising using neighboring pixels/rows, etc.) have only been applied to the later steps of conventional process  100 . 
     
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         [0011]      FIG. 1  depicts a standard XPCi imaging process; 
           [0012]      FIG. 2  depicts an X-ray imaging system in accordance with some embodiments; 
           [0013]      FIG. 3  depicts a differential XPC imaging setup in accordance with some embodiments; 
           [0014]      FIG. 4  depicts an XPCi imaging process in accordance with some embodiments; and 
           [0015]      FIGS. 5A-5C  depict regularized Fourier analysis image contrast results for an object obtained using an XPCi system in accordance with some embodiments; 
       
    
    
     DETAILED DESCRIPTION 
       [0016]    In accordance with embodiments, regularization of XPCi imaging is applied to the beginning and central steps of the XPCi imaging process. In particular, the Fourier analysis step uses the air scan Fourier coefficients as correction factors for the analysis of the object scan. Rather than independently computing Fourier coefficients for each detector pixel, prior knowledge about the object undergoing imaging is incorporated into embodying methods to yield a global image reconstruction. For example, knowing that neighboring pixels are not completely independent from one another is sufficient prior knowledge about the object. In accordance with implementations, the choices for the regularizers can include total variation, wavelet sparsity, dictionary sparsity, etc.—each of which can act on the complex coefficients, their magnitude, and/or their phase terms. 
         [0017]    In accordance with embodiments, the pixel cost function is globalized by a regularization term that incorporates prior knowledge (e.g., model assumptions) about the object. The optimization problem can then be solved using any numerical solver for unconstrained optimization, in particular a dedicated iterative shrinkage/thresholding algorithm. 
         [0018]      FIG. 2  illustrates X-ray imaging system  10  including an X-ray source  14  that projects a beam of X-rays  16  through a subject  18  (e.g., a patient, object, sample, etc.) toward one or more detectors  20 . The detector  20  is coupled to a data acquisition system  32 . The one or more detectors  20  sense the transmitted X-rays that pass through the subject  18 , and the data acquisition system  32  converts the sensed X rays to digital signals for subsequent processing. Each detector  20  produces an electrical signal that represents the intensity of an impinging X-ray beam after it passes through the subject  18 . The operation of the X-ray source  14  may be governed by an X-ray controller  34  that provides power and timing signals to the X-ray source  14 . An image reconstructor  36  receives sampled and digitized X-ray data from the data acquisition system  32  and performs reconstructions to produce phase contrast images. The reconstructed image is applied as an input to a processor based computer  40  that stores the image in a mass storage device  42 . 
         [0019]    The computer  40  also receives commands and scanning parameters from an operator via a console  44  that has some form of operator interface, such as a keyboard, mouse, voice activated controller, or any other suitable input apparatus. An associated display  46  enables the operator to observe the reconstructed images and other data from the computer  40 . The operator-supplied commands and parameters are used by the computer  40  to provide control signals and information to the data acquisition system  32  and the X-ray controller  34 . 
         [0020]      FIG. 3  illustrates a differential XPC imaging setup  48  in which a spatially coherent X-ray beam is used to probe an object (or subject)  18 . An incoherent X-ray source  14  is provided with a blocking grating  54  to create the coherent X-ray beam. In accordance with other implementations, the spatially coherent X-ray beam may be realized by synchrotron radiation, a micro focus X-ray source, or any other suitable source. A spatially coherent X-ray beam passes a phase grating  56 , and periodic interference patterns or fringes are generated. The patterns&#39; period is typically in the order of a few microns, an interferometric technique is applied to analyze the fringes using an X-ray detector  20  (e.g., having a pixel in the order of a few 100 ∥m). Another blocking grating  60  having the same period as the fringes is placed in front of the detector  20 . 
         [0021]    During operation of the illustrated imaging setup  48 , in a series of steps, grating  56  or grating  60  are shifted by a fraction of its period in the direction orthogonal to the grating slits, and images are taken for each position. After covering the entire period, the measurements for each detector pixel may be described as the convolution of the fringes with the rectangular grating function. Using Fourier analysis, the phase, the mean value and the oscillation amplitude of the fringes are determined. During an imaging operation, in addition to the gratings  56  and  60 , the object or subject  18  is placed into the X-ray beam, and the X-rays are refracted by the object  18  and hence undergo an additional phase shift. By repeating the measurement procedure, the phase of the shifted fringes is detected and the difference of both measurements yields the phase shift due to the object  18 . The differential XPC measurement generates projections of the gradient of the cumulative phase shift due to refractive index variability of the object in a direction orthogonal to the X-ray beam and to the grating slits. 
         [0022]      FIG. 4  depicts XPCi imaging process  400  in accordance with some embodiments. Imaging process  400  applies regularized Fourier analysis to first and central steps in the imaging process. To correct for system non-homogeneities, process  400  obtains air scans, step  405 , of the system prior to the presence of the object undergoing imaging. Fourier analysis of the air scans is performed, step  410 . From the Fourier analysis, air coefficients can be computed, step  415 . After the object undergoing imaging is placed in the system, object scans can be obtained, step  420 . 
         [0023]    Fourier analysis of the object scans is regularized, step  425 . Regularization uses the prior knowledge of the object (i.e., assumptions that the object measured is not completely random). By assuming ‘non-randomness’ of the object, potential regularizations can include total variation (e.g., the object has some smoothness properties) and wavelet domain sparsity (e.g., the object contains somewhat coherent structures due to the non-randomness). In accordance with other implementations, other regularizations can be used. In accordance with embodiments, once the regularized reconstruction is formulated, the solution is a global optimization problem—i.e., no longer a pixel-by-pixel approach. 
         [0024]    Regularization is done by incorporating prior object knowledge, e.g., by adding terms to the cost function. The correction factors are used to employ the regularization in a proper way. Object prior knowledge can be exploited for the real physical object, and not for the physical object confounded by inhomogeneity factors due to the system. 
         [0025]    Fourier analysis of the regularized object scans is performed, step  430 , to obtain object coefficients. In accordance with embodiments, the regularized Fourier analysis directly results, step  430 , in the coefficients p(x,y,n) (as opposed to the conventional coefficients a obj (x,y,n)), from which absorption image  432 , differential phase image  434 , and dark field image  438  can be calculated—see Equations 3-5 (above). 
         [0026]    Equation 6 (below) includes a regularization function, but does not take system non-homogeneities into account, and is confounded by variation across the detector (i.e., regularization is not optimal). 
         [0000]      min a Σ x,y ∥           H   a ( x,y ,•)− b ( x,y ,•)∥ 2 +Σ n=0   N             n ( a (•,•, n ))  (EQ.6)
 
         [0027]    where, b(x,y,s) is measured data for x, y detector row/column; 
         [0028]    s is grating step; 
         [0029]               is Fourier basis up to order N (typically N=1); 
         [0030]    a(x,y,n) is the Fourier coefficients for object scan (order n); and 
         [0031]                n  is the regularizer function (e.g., total variation, wavelet sparsity, dictionary sparsity, etc. —which act on the complex coefficients, their magnitude and/or their phase terms). 
         [0032]    In accordance with an embodiment, correction factors from the air reference scan are inserted into data obtained from the object scan. By inserting the reference scan correction factors there are less steps needed to obtain accurate imaging data, thus less exposure (i.e., dosage) to the x-rays for the same scan obtained by conventional XPCi techniques. Alternately, for the same exposure an improved image (i.e., higher image quality) can be obtained. This approach can be expressed as Equation 7: 
         [0000]      min a Σ x,y ∥           H ( Cp )( x,y ,•)− b ( x,y ,•)∥ 2 +Σ n=0   N             n ( p (•,•, n ))  (EQ.7)
 
         [0033]    where (Cp)(x,y,n)=a air (x,y,n)p(x,y,n) is the operator to uncorrect p. 
         [0034]    Equations 6 and 7 can be solved iteratively (e.g., using an iterative shrinkage/thresholding algorithm (ISTA) or other iterative approach). 
         [0035]      FIGS. 5A-5C  depict regularized Fourier analysis image contrast results for an object obtained using an XPCi system in accordance with some embodiments. Measurements were taken with 20 steps (i.e., scans) (left column); 10 steps (middle column): and 5 steps (right column). The object, a fish, was measured under the same conditions (i.e., @35 kVp, 1 sec exp. Time, +3% noise (which was added artificially to highlight the de-noising properties of the embodying method)) using standard XPCi processing (first row); embodying XPCi processing with wavelet regularization (second row); embodying XPCi processing with total variation regularization (third row); and embodying XPCi processing with wavelet plus total variation regularization (fourth row).  FIG. 5A  depicts the absorption image contrast results.  FIG. 5B  depicts the differential phase image contrast results.  FIG. 5C  depicts the dark field image contrast results. The percentages within each image represent the root-mean-square error (RMSE) from a high-quality scan (40 steps with no noise added). As can be seen from  FIGS. 5A-5C , embodying XPCi processing in accordance with embodiments produces contrast images that have lower RMSE values than the standard approach for the same number of steps. Further, in some instances embodying XPCi processing can produce lower RMSE values in fewer steps when compared to the standard XPCi processing. Accordingly, embodying XPCi processing can lower the dosage exposure of patients and still yield acceptable, or even improved, results. 
         [0036]    In accordance with some embodiments, a computer program application stored in non-volatile memory or computer-readable medium (e.g., register memory, processor cache, RAM, ROM, hard drive, flash memory, CD ROM, magnetic media, etc.) may include code or executable instructions that when executed may instruct and/or cause a controller or processor to perform methods discussed herein such as a method for regularization of XPCi imaging by exploiting prior knowledge of the object to the beginning and central steps of the XPCi imaging process, as described above. 
         [0037]    The computer-readable medium may be a non-transitory computer-readable media including all forms and types of memory and all computer-readable media except for a transitory, propagating signal. In one implementation, the non-volatile memory or computer-readable medium may be external memory. 
         [0038]    Although specific hardware and methods have been described herein, note that any number of other configurations may be provided in accordance with embodiments of the invention. Thus, while there have been shown, described, and pointed out fundamental novel features of the invention, it will be understood that various omissions, substitutions, and changes in the form and details of the illustrated embodiments, and in their operation, may be made by those skilled in the art without departing from the spirit and scope of the invention. Substitutions of elements from one embodiment to another are also fully intended and contemplated. The invention is defined solely with regard to the claims appended hereto, and equivalents of the recitations therein.