# Patent

## Title
Pre-reconstruction calibration, data correction, and material decomposition method and apparatus for photon-counting spectrally-resolving X-ray detectors and X-ray imaging

## Abstract

An apparatus and method of processing X-ray projection data obtained using photon-counting detectors and having multiple spectral components. The processing of the projection data includes correcting for nonlinear detector response, where the detector response model includes: pileup, ballistic deficit effects, polar effects, and characteristic X-ray escape. The processing of the projection data also includes a material decomposition mapping the projection data from spectral components into material components corresponding to high-Z and low-Z materials. The material decomposition includes a noise balancing process where the allocation of spectral components between a high-energy and a low-energy combination of spectral components is adjusted such that both high- and low-energy components have signal-to-noise ratios of similar magnitude. For computed tomography (CT) applications, material decomposition can be followed by image reconstruction and then image post-processing and presentation. For non-CT applications, material decomposition can be followed by image post-processing and presentation.

## Background
BACKGROUND
Field
This disclosure relates to data processing of X-ray projection data obtained using spectrally-resolving photon-counting X-ray detectors for both computed tomography (CT) and non-CT applications, and more particularly relates to the data processing steps of calibration, correcting for detector artifacts, measurement artifacts, and performing material decomposition of X-ray projection data.
Description of the Related Art
Computed tomography (CT) systems and methods are widely used, particularly for medical imaging and diagnosis. CT systems generally create images of one or more sectional slices through a subject's body. A radiation source, such as an X-ray source, irradiates the body from one side. A collimator, generally adjacent to the X-ray source, limits the angular extent of the X-ray beam, so that radiation impinging on the body is substantially confined to a planar region defining a cross-sectional slice of the body. At least one detector (and generally many more than one detector) on the opposite side of the body receives radiation transmitted through the body substantially in the plane of the slice. The attenuation of the radiation that has passed through the body is measured by processing electrical signals received from the detector.
Conventionally, energy-integrating detectors have been used to measure CT projection data. Now, recent technological developments are making photon-counting detectors a feasible alternative to conventional energy-integrating detectors. Photon-counting detectors have many advantages including their capacity for preforming spectral CT. To obtain the spectral nature of the transmitted X-ray data, the photon-counting detectors split the X-ray beam into its component energies or spectrum bins and count a number of photons in each of the bins. Since spectral CT involves the detection of transmitted X-rays at two or more energy levels, spectral CT generally includes dual-energy CT by definition.
Many clinical applications can benefit from spectral CT technology, which can provide improvement in material differentiation and beam hardening correction. Further, semiconductor-based photon-counting detectors are a promising candidate for spectral CT, which is capable of providing better spectral information compared with conventional spectral CT technology (e.g., dual-source, kVp-switching, etc.).
Photon-counting detectors are configured to acquire the spectral nature of the X-ray source. To obtain the spectral nature of the transmitted X-ray data, the photon-counting detector counts a number of photons in each of a plurality of energy bins. The use of the spectral nature of the X-ray source in CT is often referred to as spectral CT. Since spectral CT involves the detection of transmitted X-rays at two or more energy levels, spectral CT generally includes dual-energy CT by definition.
Semiconductor based photon-counting detectors used in spectral CT can detect incident photons and measure photon energy for every event. However, due to factors such as interaction depth and ballistic deficit, the measured photon energy cannot be related to incident photon energy uniquely. Furthermore, at high flux, pulse-pileup may also cause a loss in photon count and a distortion in photon energy. Accordingly, accurate image reconstruction can be achieved by efficiently estimating parameters of a response function of the photon-counting detectors.
There are several effects that can cause the detected spectrum to deviate from the X-ray spectrum incident on the photon-counting detectors, including: pileup (i.e., multiple detection events occurring within the detector response time), ballistic deficit effects, polar effects, characteristic X-ray escape, and space-charge effects.
Regarding pileup and ballistic deficit, due to the dead time ( ̃100 ns), which is determined by the type of semiconductor (e.g. CZT or Cd Te), the semiconductor thickness and readout circuitry, pulse pileup at high X-ray flux ( ̃108 cps/mm2) can be very severe, and the measured spectral signals can be distorted. The distorted spectral signal can cause artifacts in the reconstructed images. Furthermore, the dead time is not a constant for a given readout circuit due to the location of the pulse formation within the detector cell. However, if the pileup effect can be corrected in the detector model, then image quality can be improved.
Regarding the polar effect, when X-ray radiation is incident on a detector element at an oblique angle rather than normal incidence, then X-rays will enter the detector element through multiple faces of the detector element. The pileup and ballistic deficit will depend on which face the X-rays enter through. Thus, a detector response model benefits from including the differences in the pileup and ballistic deficit due to oblique X-rays illuminating multiple faces of the detector (i.e., the polar effect).
Regarding characteristic X-ray escape, when high energy photons impinge on a detector, the inner shell electrons from atoms of the detector are ejected from the atom as “photoelectrons.” After the ionization or excitation, the atom is in an excited state with a vacancy (hole) in the inner electron shell. Outer shell electrons then fall into the created holes, thereby emitting photons with energy equal to the energy difference between the two states. Since each element has a unique set of energy levels, each element emits a pattern of X-rays characteristic of the element, termed “characteristic X-rays.” The intensity of the X-rays increases with the concentration of the corresponding element.
In many materials such as Cadmium Telluride (CdTe) or Cadmium Zinc Telluride (CZT) or the like, the characteristic X-rays primarily involve K-shell (closest shell to the nucleus of an atom) electrons. If the characteristic X-rays escape from the detector, the detector signal is incorrect and the loss of energy incurred manifests itself as errors in the output spectrum of the detectors. Thus, the measured spectral signal can be distorted and may cause artifacts in the reconstructed image.
Uncorrected, each of the discussed measurement/detector artifacts distorts the detected spectrum relative to the incident spectrum ultimately degrading the quality of reconstructed images and the material decomposition derived from the data.
One advantage of spectral CT, and spectral X-ray imaging in general, is that materials having atoms with different atomic number Z also have different spectral profiles for attenuation. Thus, by measuring the attenuation at multiple X-ray energies, materials can be distinguished and the attenuation can be attributed to a particular atom (i.e., effective Z). This attribution enables spectral projection data to be mapped from the spectral domain to the material domain using a material decomposition. In some instances, this material decomposition is performed using a dual-energy analysis method.
The dual-energy analysis method can be used because the attenuation of X-rays in biological materials is dominated by two physical processes (i.e., photoelectric absorption and Compton scattering). Thus, the attenuation coefficient as a function of energy can be approximated by the decomposition

μ(E,x,y)=μPE(E,x,y)+μC(E,x,y),

where μPE(E, x, y) is the photoelectric attenuation and μC(E, x, y) is the Compton attenuation. This attenuation coefficient can be rearranged instead into a decomposition of a high-Z material (i.e., material 1) and a low-Z material (i.e., material 2) to become

μ(E,x,y)≈μ1(E)c1(x,y)+μ2(E)c2(x,y),

wherein c1,2(x,y) are spatial functions describing how much the imaged object located at position (x,y) is represented by materials 1 and 2, respectively.




BRIEF DESCRIPTION OF THE DRAWINGS
A more complete understanding of this disclosure is provided by reference to the following detailed description when considered in connection with the accompanying drawings, wherein:
FIG. 1 shows a flow diagram of an implementation of a computed tomography data processing method;
FIG. 2 shows a flow diagram of an implementation of a calibration parameter estimation method;
FIG. 3 shows a flow diagram of an implementation of a method for calculating a detector response model;
FIG. 4 shows an implementation of X-ray radiation entering a detector element at an oblique angle;
FIG. 5 shows a schematic of an implementation of a reference detector in an X-ray imaging apparatus;
FIG. 6 shows a flow diagram of an implementation of a basis material decomposition method;
FIG. 7 shows a flow diagram of an implementation of a projection length calculation method;
FIG. 8 shows a flow diagram of an implementation of an iterative detector response correction and material decomposition method; and
FIG. 9 shows a schematic of an implementation of a computed tomography scanner.



DETAILED DESCRIPTION
Referring now to the drawings, wherein like reference numerals designate identical or corresponding parts throughout the several views, FIG. 1 shows a flow diagram of a method for reconstructing an image of an object OBJ based on a series of projection measurements of the object OBJ performed at different projection directions (i.e., computed tomography (CT) using projective measurements). The data processing is a product of three inputs, including: calibration values 106, projection data 104, and reference measurement values 108. The projection data have multiple spectral components, making it compatible with material decomposition based on the different spectral absorption characteristics of high-Z and low-Z materials (i.e., greater contribution of photoelectric absorption for bone (high-Z) than for water (low-Z)). In addition to being applicable to CT applications as shown in FIG. 1, processes 110 and 120 are also applicable to non-CT applications involving projective measurements, including radiography, mammography, and tomosynthesis.
The first process 110 of the image reconstruction method 100 corrects the projection data for nonlinearities, loss mechanisms, and other aspects of the detectors and the measurement process. This can include applying various calibration values 106 and reference measurement values 108 to the projection data.
Next, the method 100 proceeds to process 120 wherein the spectral projection data is decomposed from spectral components into material components while still in the projection domain. While it is possible to first reconstruct images of the object OBJ using spectral components and then perform material decomposition in the image domain (swapping the order of process 120 with process 130 contrary to the order illustrated in FIG. 1), this alternative order of the processing steps will not be considered herein. Although processes 110 and 120 are conceptually distinct, in practice (as is discussed later) the execution of processes 110 and 120 can overlap or blend together in certain implementations of method 100.
After process 120, the method 100 proceeds to process 130 wherein multiple images are reconstructed using an image reconstruction process (e.g., an inverse Radon transformation). The image reconstruction can be performed using a back-projection method, a filtered back-projection, a Fourier-transform-based image reconstruction method, an iterative image reconstruction method (e.g., algebraic reconstruction technique or the like), a matrix-inversion image reconstruction method, or a statistical image reconstruction method. For non-CT applications (e.g., radiography, mammography, and tomosynthesis) the process 130 is omitted, and the non-CT application can proceed directly from process 120 to process 140.
After process 130, the method 100 proceeds to process 140 wherein post-processing steps are performed on the data, including: volume rendering, smoothing, filtering, and various methods for combining the material images to convey physical concept (e.g., maps of the attenuation, density, or effective Z density).
Finally, in step 150 of method 100 the image is presented to a user. The image presentation can be performed by displaying the image on a digital screen (e.g., LCD monitor), by printing the image on a suitable medium (e.g., paper or an X-ray film), or by storing the image on a computer-readable medium.
The discussion herein is focused primarily on process 110 and process 120. These processes are applicable to both CT and non-CT applications.
The projection data correction process 110 is based on a detector response function, which is a modeled by the detected energy Sout(E) derived from the energy spectrum incident on the detector Sin(E), wherein an implementation of the detector response function is given by

Sout(E)=ne−nτ∫dE0R0(E,E0)Sin(E0)+n2e−nτ∫∫dE0dE1R1(E,E0,E1)Sin(E0)Sin(E1)

wherein R0 is the linear response function, R1 is the quadratic response function representing first-order pileup, and T is the dead time of the detector. Each of R0, R1, and τ can depend on the detector element and the incident angle of the X-ray radiation. The incident spectrum is given by

Sin(Ei)=Sair(E)exp[−μ1(E)L1−μ2(E)L2],

wherein μ1 and μ2 are the attenuation coefficients of the basis materials for the material decomposition, L1 and L2 are the projection lengths. The X-ray flux n for each detector is given by

n=nair∫dE0Sin(E0)exp[−μ1(E0)L1−μ2(E0)L2],

wherein nair=A·Iref is the calculated flux based on a reference intensity measurement Iref of an X-ray source and A is a calibration value. The factor nair represents the X-ray flux at the detector for projective measurements taken in the absence of an imaged object OBJ (i.e., in the presence of air as an image object OBJ). A unique value of Iref is measured by a detector near the X-ray source for each projection angle.

The number of counts in a given energy bin is calculated by the expression

Nm=T∫dEwm(E)SDet(E),

wherein T is the integration time and wm(E) is the spectral function of the mth energy bin of the photon counting detectors. For example, the spectral function could be a square function, which is defined as





  
    
      w
      m
    
    ⁡
    
      (
      E
      )
    
  
  =
  
    {
    
      
        
          
            1
          
          
            
              
                E
                
                  m
                  ⁢
                  
                    
                  
                  ⁢
                  1
                
              
              <
              E
              <
              
                E
                
                  m
                  ⁢
                  
                    
                  
                  ⁢
                  2
                
              
            
          
        
        
          
            0
          
          
            otherwise
          
        
      
      .
    
  




The values expressed by the symbols R0, R1, and τ, μ1, μ2, A, Iref, Sair, L1 and L2 can be categorized into four types: (1) constant values (fixed prior to calibrations or projective measurements), (2) calibration values (determined by calibrations performed prior to projection measurements on imaged object OBJ), (3) reference values (calibrations during the projection measurements on imaged object OBJ), and (4) result values (values being solved for when the detector model is solved).
The constants are the attenuation coefficients μ1, and μ2, which are determined by established models for the chosen high-Z (e.g., bone) and the low-Z (e.g., water) materials. There is only one reference value—the reference intensity Iref. The results values are the projection lengths L1 and L2. Thus, all other values/functions (i.e., R0, R1, τ, A, and Sair) are calibration values.
FIG. 2 shows as implementation of a method 200 of determining R0, R1, and τ, according to U.S. patent application Ser. No. 14/479,955 incorporated herein by reference in its entirety. In U.S. patent application Ser. No. 13/866,965, incorporated herein by reference in its entirety, the method 200 is applied to determining calibration constants. In U.S. patent application Ser. No. 14/535,396, incorporated herein by reference in its entirety, the ballistic deficit model is expanded to include polar effects. By including polar effects, the detector model is generalized to include variations in the detector response when the X-ray radiation impinges the detectors at angles deviating from normal incidents.
Models accounting for pileup and the ballistic effect describe how the recorded spectrum Sout differs from the incident spectrum Sin due to multiple detection events occurring within a given detection window (i.e., the response time of the detector), resulting in a larger detection signal resulting from the combination of photoelectrons coming from multiple detection events being attributed to a single detection event. Therefore, pileup shifts the recorded spectrum SM towards higher energies because higher energies arising from multiple detections (i.e. the pileup of detection events within a detection time window) are tallied as high-energy single detection events. Thus, the detected spectrum Sout deviates from the incident spectrum Sin.
In contrast, characteristic X-ray escape causes the recorded spectrum Sout to shift towards lower energies relative to the incident spectrum Sin because some of the energy absorbed into the semiconductor in a detection event is reemitted as a characteristic X-ray. Thus, rather than all of the absorbed X-ray energy being applied to creating photoelectrons, some of the energy is reemitted at a characteristic energy corresponding to the difference between an empty inner electron shell (e.g., the K-, L-, or M-shell) and an outer shell. The characteristic X-ray escape is often referred to as K-escape due to the characteristic X-ray energy typically corresponding to the K-edge. The detector model for K-escape is described in U.S. patent application Ser. No. 14/190,170, incorporated herein by reference in its entirety.
FIG. 3 shows an implementation of a process 310 that models the combined effects of the energy spectrum distortions arising from characteristic pileup, ballistic deficit, and polar effects. The first step 322 of process 310 accounts for the ballistic deficit, and step 324 of process 310 accounts for the spectral shifts arising from the polar effect.
In one implementation, without considering K-escape, the probability that energy E is deposited at depth z0 is,




  
    
      S
      dep
    
    ⁡
    
      (
      
        E
        ,
        
          z
          0
        
      
      )
    
  
  =
  
    
      
        S
        
          i
          ⁢
          
            
          
          ⁢
          n
        
      
      ⁡
      
        (
        E
        )
      
    
    ⁢
    
      e
      
        
          -
          
            
              μ
              CZT
            
            ⁡
            
              (
              E
              )
            
          
        
        ⁢
        
          
            z
            0
          
          /
          sin
        
        ⁢
        
          
        
        ⁢
        β
      
    
    ⁢
    
      
        
          
            μ
            CZT
          
          ⁡
          
            (
            E
            )
          
        
        
          sin
          ⁢
          
            
          
          ⁢
          β
        
      
      .
    
  




Returning to FIG. 2, the diagram illustrates an apparatus for determining a detector pileup model for each photon-counting detector in a spectral CT scanner. In particular, FIG. 2 illustrates a detector pileup model 201 that receives an incident spectrum Sin(E), a count rate, and a parameter vector a. Based on the received values, the detector pileup model 201 generates a simulated measured spectrum Sout(E; a), wherein a includes R0, R1, and τ. The process of determining the simulated measured spectrum Sout(E; a) will be described in more detail below.
As shown in FIG. 2, the model parameter estimation cost function 203 compares the simulated measured spectrum Sout(E; a) with an actual measured spectrum SDet(E), and updates a parameter vector a so as to minimize a predetermined cost function. The updated parameter vector a, which includes R0, R1, and τ, is fed back to the detector pileup model to generate a new simulated measured spectrum Sout(E; a). This process continues for a predetermined number of iterations or until the change in the parameter vector a falls below a predetermined threshold. A respective optimal parameter a is found for each photon-counting detector in a scanner.
The input spectrum Sin(E) can be determined by calculation (all vendors have models to calculate the output from their tubes) or measurements (by using a gold-standard spectroscopic detector, e.g., a high-purity germanium spectrometer) for each PCD in a scanner. The measured spectrum SDet is the output spectrum from each PCD corresponding to each incident spectrum. Where different X-ray sources are used for calibration measurements and for projection measurements of the imaged object OBJ, then Sin(E) and Sair(E) will be different. However, where the same X-ray source is used for both calibration and projection measurements, then Sin(E) is Sair(E). Whether or not the X-ray source for projection measurements is also used for calibrations, Sair(E) is determined prior to the projection measurements of the imaged object OBJ. Sair(E) may be determined in a similar manner to Sin(E) is determined. One of ordinary skill in the art will recognize many ways of obtaining Sair(E).
The parameter vector a includes a dead time value τ, a time threshold T to determine whether, e.g., double photon events are peak pileup events or tail pileup events (although this threshold applies to determining whether a peak- or tail-pileup event occurs at any pileup order), and individual probabilities of different number of quasi-coincident photons χ0, χ1p, χ1t, χ2p, χ2t, etc. For example, χ0 is the probability of single photon events, χ1p is the probability of peak double pileup events, χ1t is the probability of tail double pileup events, χ2p is the probability of triple peak pileup events, etc. Note that, in practice, the contribution of higher order spectra drops quickly with increasing order and, in most embodiments, those spectra can be ignored in calculating the summed spectrum. Note that the sum of individual probabilities will be equal to or less than 1.
For the no pileup case, the first component spectrum is calculated as follows. The component spectrum becomes:

S0(E)=e−nτd∫∫dz0dE0χ0Sdep(E0,z0)δ(E−νp(tTOF0;z0,E0))

wherein the integration runs over the whole volume with the energy condition determined by δ(•) the Dirac delta function. The time of flight tTOF0 and the pre-amplifier voltage νp are described in U.S. patent application Ser. No. 13/866,965.

For peak pileup, the spectrum is




  
    
      S
      1
      p
    
    ⁡
    
      (
      E
      )
    
  
  =
  
    n
    ⁢
    
      
    
    ⁢
    
      e
      
        
          -
          n
        
        ⁢
        
          
        
        ⁢
        
          τ
          d
        
      
    
    ⁢
    
      ∫
      
        ∫
        
          ∫
          
            ∫
            
              ∫
              
                
                  d
                  
                    z
                    0
                  
                
                ⁢
                
                  d
                  
                    E
                    0
                  
                
                ⁢
                
                  d
                  
                    t
                    1
                  
                
                ⁢
                
                  d
                  
                    z
                    1
                  
                
                ⁢
                
                  d
                  
                    E
                    1
                  
                
                ⁢
                
                  χ
                  1
                  p
                
                ⁢
                
                  
                    S
                    dep
                  
                  ⁡
                  
                    (
                    
                      
                        E
                        0
                      
                      ,
                      
                        z
                        0
                      
                    
                    )
                  
                
                ⁢
                
                  
                    S
                    dep
                  
                  ⁡
                  
                    (
                    
                      
                        E
                        1
                      
                      ,
                      
                        z
                        1
                      
                    
                    )
                  
                
                ×
                
                  
                    δ
                    (
                    
                      E
                      -
                      
                        
                          
                            v
                            p
                          
                          ⁡
                          
                            (
                            
                              
                                t
                                TOF
                                0
                              
                              ;
                              
                                z
                                0
                              
                            
                            )
                          
                        
                        ⁢
                        
                          e
                          
                            -
                            
                              
                                
                                  t
                                  max
                                
                                -
                                
                                  t
                                  TOF
                                  0
                                
                              
                              
                                τ
                                p
                              
                            
                          
                        
                      
                      -
                      
                        
                          
                            v
                            p
                          
                          ⁡
                          
                            (
                            
                              
                                t
                                TOF
                                1
                              
                              ;
                              
                                z
                                1
                              
                            
                            )
                          
                        
                        ⁢
                        
                          e
                          
                            -
                            
                              
                                
                                  t
                                  max
                                
                                -
                                
                                  t
                                  TOF
                                  1
                                
                                -
                                
                                  t
                                  1
                                
                              
                              
                                τ
                                p
                              
                            
                          
                        
                      
                    
                    )
                  
                  .
                
              
            
          
        
      
    
  




Similar spectra can be calculated for tail pileup and higher order pileup. The total output spectrum then becomes the sum of each of the individual pileup components,

Sout=S0(E)+S1p(E)+S1t0(E)+S1t1(E)+S2p(E)+ . . . .

The ballistic deficit detector model from U.S. patent application Ser. No. 13/866,965 can be improved by incorporating a polar effect model from U.S. patent application Ser. No. 14/535,396. The polar effect model generalizes the ballistic deficit model to include the case where oblique X-rays enter a detector element through faces other than the front face of the detector, as shown in FIG. 4. The polar effect is included by modifying the ballistic deficit expressions for expressions for S0(E), S1p(E), S1t0(E), etc. to include the effect that some X-rays deviating from normal are incident on different faces of the detector. Further, X-rays incident on different detector faces have different propagation lengths through the detector material and different migration lengths for photoelectrons from their point of origin to the anode of the detector. These differences result in different linear and nonlinear detector response functions.
For example, FIG. 4 shows an example of a detector element with incident X-rays deviating from normal incidence. The X-rays are incident on both the top surface and the side surface. Those X-rays incident and absorbed on the side surface have, on average, shorter propagation lengths in the semiconductor of the detector element and a short distance between the origin of the photoelectrons and the detector anode.
Therefore the linear and nonlinear detector response and the corresponding models can change depending on the X-ray angle incidence. Details of how the detector response model can be modified to include the polar effect are described in U.S. patent application Ser. No. 14/535,396, which describes how these polar angle corrections can be applied to the expressions discussed for detector response based on the ballistic deficit.
In one implementation, the calibrations also include a method for reference calibration, as discussed in U.S. patent application Ser. No. 14/103,137, incorporated herein by reference in its entirety. A reference calibration is utilized to establish a mapping between a reference detector signal and a true count incident on a photon-counting detector. The true count rate/flux is then utilized to determine projection data from the detector, which is in turn utilized to generate images.
FIG. 5 illustrates a CT scanner that utilizes a reference detector to monitor a variation of the X-ray source during a scan. The CT scanner generally includes an X-ray source, a bowtie filter, a reference detector, and an arrangement of photon counting detectors (PCDs) shown as PCD1 through PCDN, where PCDN is the Nth PCD. The imaged object OBJ is positioned between the X-ray source and the PCDs. As illustrated in FIG. 5, the reference detector can be located between the bowtie filter and the object in one implementation. The reference detector can also be located between the X-ray source and the bowtie filter, as also illustrated in FIG. 5. In these and other implementations, the reference detector can be installed so as to rotate with the bowtie filter and the X-ray source. Further, multiple reference detectors can be provided, and outputs thereof can be used to create an average reference signal that is used for reference calibration. Also, multiple reference calibrations can be obtained by using multiple reference detectors, where the reference calibrations can be averaged or statistically changed to obtain a summary reference calibration.
According to aspects of implementations of this disclosure, an additional filter can be additionally/optionally included. Further, a phantom can be additionally/optionally included, as illustrated in FIG. 5. Also, in one implementation, the bowtie filter can be removed or not used.
In one implementation, the CT scanner can use a non-spectral reference detector. The reference detector is an energy-integrating detector, and only measures X-ray intensity variations. The reference detector does not measure spectrum variation.
For a given incident spectrum, the following relationship exists: nPCDair∝Iref.
A signal from an energy-integrating reference detector is Iref=T·nrefair·∫dE·E·Sref(E), where T is the integration time, nrefair is the count rate at the reference detector, E is the energy factor, and Sref (E) is the normalized spectrum at the reference detector. From FIG. 5, given a same spectrum, nPCDair∝nrefair since they only differ by geometrical factors (i.e., different distances to the X-ray tube) and attenuation paths to the X-ray tube (different bowtie paths, filter paths, etc.). Therefore, nPCDair∝Iref.
A constant A is determined such that nPCDair=A·Iref, where Iref is the reference detector signal (arbitrary unit), and nPCDair is the true count rate (i.e., the true count rate across all energies) on the PCD without an object. nPCD is then determined from nPCDair and a basis material thickness {L1, L2} according to

nPCD=nPCDair·∫dE·Sair(E)e−μ1(E)L1−μ2(E)L2,

where Sair(E) is the normalized X-ray source spectrum in the absence of the object OBJ (i.e., the object OBJ is replaced by air), which spectrum Sair(E) is known to the manufacturer.

Having determined each of the calibration values R0, R1, τ, A, and Sair during the calibration measurements and analysis, the next step of method 100 is to perform the projection data correction process 110. In process 110 the correction factor is given by

NmCorr.=NmRaw−NmNonlin. 

where NmCorr. is the corrected count value of the mth energy bin of the PCD, NmRaw is the raw count value recorded from the detector, and NmNonlin. is the calculated count from the nonlinear detector response. In order to calculate NmNonlin. the values for the projections lengths L1 and L2 have to be assumed.

The nonlinear count value NmNonlin. is calculated according to

NmNonlin.=∫dEwm(E)SNonlin.(E),

where, in one implementation, the nonlinear spectrum correction is given by the first order pileup

SNonlin.(E)≡S1,out(E)n2e−nτ∫∫(dE0dE1R1(E,E0,E1)Sin(E0)Sin(E1).

In another implementation, the nonlinear spectrum correction includes higher-order pileup terms.
In one implementation, the nonlinear correction can be expressed as

SCorr.(E)=SRaw(E)−S1,out(E).

In one implementation, the energy spectrum distortion due to the linear term in the detector response model can also be corrected. For example, in the case that the detector counts are arranged into five energy bins, the linear response of the detector can be expressed as a matrix equation given by

{right arrow over (N)}Det=R{right arrow over (N)}ln,

where




  
    
      
        R
        _
      
      =
      
        [
        
          
            
              
                R
                
                  1
                  ,
                  1
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  1
                  ,
                  2
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  1
                  ,
                  3
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  1
                  ,
                  4
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  1
                  ,
                  5
                
                
                  (
                  0
                  )
                
              
            
          
          
            
              
                R
                
                  2
                  ,
                  1
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  2
                  ,
                  2
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  2
                  ,
                  3
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  2
                  ,
                  4
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  2
                  ,
                  5
                
                
                  (
                  0
                  )
                
              
            
          
          
            
              
                R
                
                  3
                  ,
                  1
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  3
                  ,
                  2
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  3
                  ,
                  3
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  3
                  ,
                  4
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  3
                  ,
                  5
                
                
                  (
                  0
                  )
                
              
            
          
          
            
              
                R
                
                  4
                  ,
                  1
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  4
                  ,
                  2
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  4
                  ,
                  3
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  4
                  ,
                  4
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  4
                  ,
                  5
                
                
                  (
                  0
                  )
                
              
            
          
          
            
              
                R
                
                  5
                  ,
                  1
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  5
                  ,
                  2
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  5
                  ,
                  3
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  5
                  ,
                  4
                
                
                  (
                  0
                  )
                
              
            
            
              
                R
                
                  5
                  ,
                  5
                
                
                  (
                  0
                  )
                
              
            
          
        
        ]
      
    
    ,
    
      
        
          N
          →
        
        Det
      
      =
      
        [
        
          
            
              
                N
                1
                
                  (
                  Det
                  )
                
              
            
          
          
            
              
                N
                2
                
                  (
                  Det
                  )
                
              
            
          
          
            
              
                N
                3
                
                  (
                  Det
                  )
                
              
            
          
          
            
              
                N
                4
                
                  (
                  Det
                  )
                
              
            
          
          
            
              
                N
                5
                
                  (
                  Det
                  )
                
              
            
          
        
        ]
      
    
    ,
    and
  




  
    
      
        N
        →
      
      
        I
        ⁢
        
          
        
        ⁢
        n
      
    
    =
    
      
        [
        
          
            
              
                N
                1
                
                  (
                  In
                  )
                
              
            
          
          
            
              
                N
                2
                
                  (
                  In
                  )
                
              
            
          
          
            
              
                N
                3
                
                  (
                  
                    I
                    ⁢
                    
                      
                    
                    ⁢
                    n
                  
                  )
                
              
            
          
          
            
              
                N
                4
                
                  (
                  In
                  )
                
              
            
          
          
            
              
                N
                5
                
                  (
                  In
                  )
                
              
            
          
        
        ]
      
      .
    
  



The linear response corrected counts {right arrow over (N)}ln can be obtained by inverting the matrix R to solve the matrix equation, as given by

{right arrow over (N)}ln=R−1{right arrow over (N)}Det.

In one implementation, the first correction step is correcting for the nonlinear detector response to obtain the nonlinear corrected counts,

{right arrow over (N)}Det={right arrow over (N)}Corr=[N1Corr.N2Corr.N3Corr.N4Corr.N5Corr.]T,

where T denotes the transpose. The second correction step is correcting for the linear detector response to obtain

{right arrow over (N)}ln=R−1{right arrow over (N)}Det.

After correcting for the detector response, the next step is the basis material decomposition 120. As shown in FIG. 6, the material decomposition process 120 begins with the step of apportioning the counts from the multiple energy bins into two energy components corresponding to a high spectrum SH(E) and a low spectrum SL(E). The noise balancing process of apportioning the counts into high- and low-energy components is described in U.S. patent application Ser. No. 13/906,110, incorporated herein by reference in its entirety. Noise balancing enables the delineation between the high spectrum and the low spectrum to be different for different detectors in order to prevent the case where the signal to noise ratio for either the high- or low-energy components counts became significantly lower than the other energy component (i.e., the signal-to-noise ratio of the high- and low-energy components become unbalanced) and degrades the image quality of the reconstructed image.
In FIG. 6, the noise balancing in step 610 results in partitioning the counts from the energy bins into high- and low-energy components according to



  
    
      
        N
        H
      
      =
      
        
          ∑
          m
        
        ⁢
        
          
            a
            m
            
              (
              H
              )
            
          
          ⁢
          
            N
            m
            
              Corr
              .
            
          
        
      
    
    ,
    and
  




  
    
      
        N
        L
      
      =
      
        
          ∑
          m
        
        ⁢
        
          
            a
            m
            
              (
              L
              )
            
          
          ⁢
          
            N
            m
            
              Corr
              .
            
          
        
      
    
    ,
  



where am(H)+am(L)=1 and the values am(H) and am(L) are determined by noise balancing.

In an alternative implementation, the values am(H) and am(L) are predetermined values that are not determined by a noise-balancing process.
The next step 620 is the projection length calculation, and includes the process of decomposing the high- and low-energy components projection data into projection lengths of high- and low-Z materials, as described in Yu Zou and Michael Silver, “Analysis of fast kV-switching in dual energy CT using a pre-reconstruction decomposition technique,” Proc. SPIE Vol. 6913, 691313-1 (2008), incorporated herein by reference in its entirety. One implementation of the material decomposition method of process 620 is also discussed in U.S. patent application Ser. No. 10/0,189,212, incorporated herein by reference in its entirety.
One implementation of the process 620 for material decomposition is shown in FIG. 7. Similar to the detector counts for the high- and low-energy components, the high and low spectra can be given by



  
    
      
        
          S
          H
        
        ⁡
        
          (
          E
          )
        
      
      =
      
        
          ∑
          m
        
        ⁢
        
          
            a
            m
            
              (
              H
              )
            
          
          ⁢
          
            
              S
              air
            
            ⁡
            
              (
              E
              )
            
          
        
      
    
    ,
    and
  




  
    
      
        
          S
          L
        
        ⁡
        
          (
          E
          )
        
      
      =
      
        
          ∑
          m
        
        ⁢
        
          
            a
            m
            
              (
              L
              )
            
          
          ⁢
          
            
              S
              air
            
            ⁡
            
              (
              E
              )
            
          
        
      
    
    ,
  



where SH(E) and SL(E) are respectively the detected high- and low-energy spectra in the absence of the object OBJ (i.e., the object OBJ is air), and where SH(E) and SL(E) have been normalized such that

∫dESH(E)=∫dESL(E)=1.

By taking the natural logarithm of the detector counts, the log-projection data are obtained as

gH(l)=−ln(NH/NHair) and

gL(l)=−ln(NL/NLair),

where the path length l signifies the trajectory corresponding to the path integral/Radon transform of the X-rays incident on a given detector element.

In one implementation, as shown in FIG. 7, L1 and L2 are found using perturbation theory by treating the variations around the mean of the attenuation coefficients μ1(E) and μ2(E) as perturbations. First, the mean attenuation for the high- and low-energy spectra are given by

μ1,2H,L=∫SH,L(E)μ1,2(E)dE, 

and the variations around the mean are given by

Δμ1,2H,L(E)=μ1,2(E)−μ1,2H,L.

Thus, the log-projection data can be expressed as

gH(l)=−ln∫SH(E)exp[−μ1HL1(l)−Δμ1H(E)L1(l)−μ2HL2(l)−Δμ2H(E)L2(l)]dE 

gL(l)=−ln∫SL(E)exp[−μ1LL1(l)−Δμ1L(E)L1(l)−μ2LL2(l)−Δμ2L(E)L2(l)]dE. 

Simplifying these expressions, the log-projection data can be written as

gH(l)=μ1HL1(l)+μ2HL2(l)−gH(BH)(L1(l),L2(l))

gL(l)=μ1LL1(l)+μ2LL2(l)−gL(BH)(L1(l),L2(l))

where

gH,L(BH)(L1(l),L2(l))≡ln∫SH,L(E)exp[−L1(l)Δμ1H,L(E)−L2(l)Δμ2H,L(E)]dE 

is the beam-hardening perturbation.

The first step 712 of process 620 is initializing n=0 and also initializing the values L1 and L2 to their respective zeroth order perturbation values by solving the matrix equation




  
    (
    
      
        
          
            g
            H
          
        
      
      
        
          
            g
            L
          
        
      
    
    )
  
  =
  
    
      (
      
        
          
            
              
                μ
                _
              
              1
              H
            
          
          
            
              
                μ
                _
              
              2
              H
            
          
        
        
          
            
              
                μ
                _
              
              1
              L
            
          
          
            
              
                μ
                _
              
              2
              L
            
          
        
      
      )
    
    ⁢
    
      (
      
        
          
            
              L
              1
            
          
        
        
          
            
              L
              2
            
          
        
      
      )
    
  




to obtain





  
    (
    
      
        
          
            L
            1
          
        
      
      
        
          
            L
            2
          
        
      
    
    )
  
  =
  
    
      
        D
        
          -
          1
        
      
      ⁡
      
        (
        
          
            
              
                
                  μ
                  _
                
                2
                L
              
            
            
              
                -
                
                  
                    μ
                    _
                  
                  2
                  H
                
              
            
          
          
            
              
                -
                
                  
                    μ
                    _
                  
                  1
                  L
                
              
            
            
              
                
                  μ
                  _
                
                1
                H
              
            
          
        
        )
      
    
    ⁢
    
      (
      
        
          
            
              g
              H
            
          
        
        
          
            
              g
              L
            
          
        
      
      )
    
  




where D is the determinant D=μ1Hμ2L−μ1Lμ2H.

The second step 714 of process 620 and the first step in the iterative loop is updating the beam-hardening perturbation values using the nth order perturbation in the equation

gH,L(BH)(L1(l),L2(l))≡ln∫SH,L(E)exp[−L1(l)Δμ1H,L(E)−L2(l)Δμ2H,L(E)]dE. 

The third step 716 of process 620 is to update the values of L1 and L2 by solving for the n+1th perturbation by solving the matrix equation



  
    
      (
      
        
          
            
              
                g
                H
              
              +
              
                
                  g
                  H
                  
                    (
                    BH
                    )
                  
                
                ⁡
                
                  (
                  
                    
                      L
                      1
                    
                    ,
                    
                      L
                      2
                    
                  
                  )
                
              
            
          
        
        
          
            
              
                g
                L
              
              +
              
                
                  g
                  L
                  
                    (
                    BH
                    )
                  
                
                ⁡
                
                  (
                  
                    
                      L
                      1
                    
                    ,
                    
                      L
                      2
                    
                  
                  )
                
              
            
          
        
      
      )
    
    =
    
      
        (
        
          
            
              
                
                  μ
                  _
                
                1
                H
              
            
            
              
                
                  μ
                  _
                
                2
                H
              
            
          
          
            
              
                
                  μ
                  _
                
                1
                L
              
            
            
              
                
                  μ
                  _
                
                2
                L
              
            
          
        
        )
      
      ⁢
      
        (
        
          
            
              
                L
                1
              
            
          
          
            
              
                L
                2
              
            
          
        
        )
      
    
  




  
    
      to
      ⁢
      
        
      
      ⁢
      
        obtain
        ⁢
        
          

        
        (
        
          
            
              
                L
                1
                n
              
            
          
          
            
              
                L
                s
                n
              
            
          
        
        )
      
    
    =
    
      
        
          D
          
            -
            1
          
        
        ⁡
        
          (
          
            
              
                
                  
                    μ
                    _
                  
                  2
                  L
                
              
              
                
                  -
                  
                    
                      μ
                      _
                    
                    2
                    H
                  
                
              
            
            
              
                
                  -
                  
                    
                      μ
                      _
                    
                    1
                    L
                  
                
              
              
                
                  
                    μ
                    _
                  
                  1
                  H
                
              
            
          
          )
        
      
      ⁢
      
        
          (
          
            
              
                
                  
                    g
                    H
                  
                  +
                  
                    
                      g
                      H
                      
                        (
                        BH
                        )
                      
                    
                    ⁡
                    
                      (
                      
                        
                          L
                          1
                          
                            n
                            -
                            1
                          
                        
                        ,
                        
                          L
                          2
                          
                            n
                            -
                            1
                          
                        
                      
                      )
                    
                  
                
              
            
            
              
                
                  
                    g
                    L
                  
                  +
                  
                    
                      g
                      L
                      
                        (
                        BH
                        )
                      
                    
                    ⁡
                    
                      (
                      
                        
                          L
                          1
                          
                            n
                            -
                            1
                          
                        
                        ,
                        
                          L
                          2
                          
                            n
                            -
                            1
                          
                        
                      
                      )
                    
                  
                
              
            
          
          )
        
        .
      
    
  



After step 718, step 720 of process 620 inquiries whether stopping criteria have been satisfied. In one implementation, the stopping criteria are satisfied when the values L1 and L2 have adequately converged according to some predetermined criteria, such as whether the difference between each current and previous values of L1 and L2 are less than a predefined threshold, or whether a maximum number of iterations has been reached. If stopping criteria have not been satisfied, then the loop variable n is incremented at step 718 and the loop begins again starting from step 714.
The process 620 of process 120 provides new values of the projection lengths L1 and L2 different from the original projection length values used to calculate the nonlinear correction in process 110 shown in FIG. 1. An iterative process of feeding back the projection lengths calculated in the material decomposition process 120 back to the beginning of process 110 and repeating processes 110 and 120 multiple times will result in convergence to an optimal value for the projection lengths.
Therefore, in one implementation, as shown in FIG. 8, the projection lengths L1 and L2 from the material decomposition step 850 will be fed back to the projection data correction step 820 and steps 820 through 850 are repeated again and again until the projection lengths converge such that the projection lengths remain essentially unchanged from the beginning of step 820 through the end of step 850.
FIG. 8 shows a single iterative method 800 replacing both the projection data correction process 110 and the basis material decomposition process 120. The method 800 begins at step 810 with an initial guess of the projection lengths L1 and L2 for a given PCD and a given projection measurement angle. The loop begins at step 820 by calculating the count rate

n=nair∫dE0Sair(E0)exp[−μ1(E0)L1−μ2(E0)L2],

and the nonlinear correction term (e.g., SNonlin.(E)), which in the implementation shown in FIG. 8 includes only the first order pileup (i.e., two detection events within the detection window) given by

S1,out(E)=∫∫dE0dE1R1(E,E0,E1)Sin(E0)Sin(E1),

wherein

Sin(E)=Sair(E)exp[−μ1(E)L1−μ2(E)L2].

Next at step 830 of method 800, the corrected detector spectrum is calculated correcting for pileup according to

SCorr.(E)=SRaw(E)−SNonlin.(E).

The method 800 then proceeds to step 840, wherein noise balancing is performed by mapping the detector counts into high- and low-energy components in preparation for material decomposition. Material decomposition is then performed in step 850, wherein new values for the projection lengths L1 and L2 are calculated.

Finally, at step 860, an inquiry is made into whether the stopping criteria have been satisfied. The stopping criteria can include, e.g., conditions regarding whether the projection lengths L1 and L2 have converged, and whether the maximum number of loop iterations have been reached.
The material decomposition step 850 has a loop (referred to herein as the small loop)—given by process 620 as shown in FIG. 6—with its own stopping criteria. The small-loop stopping criteria of the step 850 may depend on the current iteration number of the large loop (i.e., the large being the loop that includes steps 820 through 860). For example, in one implementation, the step 850 loop can have a lower number of maximum iterations and more relaxed convergence threshold when the large-loop iteration variable is small, preventing the small loop in step 850 from spending an excessive time optimizing the projection lengths when the value of the projection lengths are at best rough approximations based on the bigger picture of the status of the large loop.
After obtaining the projections lengths by either exiting step 860 of method 800 or by completing process 120 in method 100, the projection lengths are prepared for further processing.
In a CT process, the projection lengths are processed by reconstructing an image from the projection lengths, as illustrated in process 130 of method 100 shown in FIG. 1. The image reconstruction process can be performed using filtered back-projection, iterative image reconstruction methods, or stochastic image reconstruction methods. After reconstructing an image of the object OBJ, the CT process then proceeds to the post-processing and image presentation steps as discussed in relation to FIG. 1.
In a non-CT application, such as radiography where the projection image is the final result, the data processing method proceeds directly from the material decomposition process 120 (or alternatively step 860 of method 800) to the post-processing 140 and image presentation 150 steps illustrated in FIG. 1.
FIG. 9 shows a computed tomography (CT) scanner having both energy-integrating detectors arranged in a third-generation geometry and photon-counting detectors arranged in a fourth-generation geometry. Illustrated in FIG. 9 is an implementation for placing the photon-counting detectors (PCDs) in a predetermined fourth-generation geometry in combination with a detector unit 903 in a predetermined third-generation geometry in a CT scanner system. The diagram illustrates relative positions among an object OBJ to be scanned resting on a table 916, an X-ray source 912, a collimator/filter 914, an X-ray detector 903, and the photon-counting detectors PCD1 through PCDN in a gantry 940. Also shown in FIG. 9 is circuitry and hardware for acquiring, storing, processing, and distributing X-ray projection data. The circuitry and hardware include: a processor 970, a network controller 980, a memory 978, and a data acquisition system 976.
In one implementation, the CT scanner includes PCDs but does not include an energy-integrating detector unit 903
In general, the photon-counting detectors PCD1 through PCDN each output a photon count for each of a predetermined number of energy bins. In addition to the photon-counting detectors PCD1 through PCDN arranged in the fourth-generation geometry, the implementation shown in FIG. 9 includes a detector unit 903 having energy-integrating detectors arranged in a conventional third-generation geometry. The detector elements in the detector unit 903 can be more densely placed along the detector unit surface than the photon-counting detectors.
In one implementation, the photon-counting detectors are sparsely placed around the object OBJ in a predetermined geometry such as a circle. For example, the photon-counting detectors PCD1 through PCDN are fixedly placed on a predetermined circular component 920 in the gantry 940. In one implementation, the photon-counting detectors PCD1 through PCDN are fixedly placed on the circular component 920 at predetermined equidistant positions. In an alternative implementation, the photon-counting detectors PCD1 through PCDN are fixedly placed on the circular component 910 at predetermined non-equidistant positions. The circular component 920 remains stationary with respect to the object OBJ and does not rotate during the data acquisition.
Both the X-ray source 912, collimator 914 (e.g., a bow tie filter), and the detector unit 903 rotate around the object OBJ while the photon-counting detectors PCD1 through PCDN are stationary with respect to the object OBJ. In one implementation, the X-ray source 912 and collimator 914 are mounted on a first rotating portion 910 such as the annular frame in the gantry 940 so that the X-ray source 912 projects X-ray radiation with a predetermined source fan beam angle θA towards the object OBJ while the X-ray source 912 rotates around the object OBJ inside the sparsely placed photon-counting detectors PCD1 through PCDN. Furthermore, an additional detector unit 903 is mounted on a second rotating portion 930 in the third-generation geometry. The rotating portion 930 mounts the detector unit 903 at a diametrically opposed position from the X-ray source 912 across the object OBJ and rotates outside the stationary circular component 920, on which the photon-counting detectors PCD1 through PCDN are fixedly placed in a predetermined sparse manner.
In one implementation, the X-ray source 912 optionally travels a helical path relative to the object OBJ, which is moved in a predetermined direction that is perpendicular to the rotational plane of the rotating portion 910.
As the X-ray source 912 and the detector unit 903 rotate around the object OBJ, the photon-counting detectors PCDs and the detector unit 903 respectively detect the transmitted X-ray radiation during data acquisition. The photon-counting detectors PCD1 through PCDN intermittently detect the X-ray radiation that has been transmitted through the object OBJ and individually output a count value representing a number of photons, for each of the predetermined energy bins. On the other hand, the detector elements in the detector unit 903 continuously detect the X-ray radiation that has been transmitted through the object OBJ and output the detected signals as the detector unit 903 rotates. In one implementation, the detector unit 903 has densely placed energy-integrating detectors in predetermined channel and segment directions on the detector unit surface.
In one implementation, the X-ray source 912, the photon-counting detectors and the detector unit 903 collectively form three predetermined circular paths that differ in radius. The photon-counting detectors are sparsely placed along a first circular path around the object OBJ while at least one X-ray source 912 rotates along a second circular path around the object OBJ. Further, the detector unit 903 travels along a third circular path. The above exemplary embodiment illustrates that the third circular path is the largest and outside the first and second circular paths around the object OBJ. Although not illustrated, an alternative embodiment optionally changes the relative relation of the first and second circular paths so that the second circular path for the X-ray source 912 is larger and outside the first circular path of the sparsely placed photon-counting detectors PCD1 through PCDN around the object OBJ. Furthermore, in another alternative embodiment, the X-ray source 912 also optionally travels on the same third circular path as the detector unit 903.
There are other alternative embodiments for placing the photon-counting detectors in a predetermined fourth-generation geometry in combination with the detector unit in a predetermined third-generation geometry in the CT scanner. Several alternative embodiments of the X-ray CT Scanner as described in U.S. patent application Ser. No. 13/0,291,097, herein incorporated by reference in its entirety.
In one implementation, the X-ray source 912 is optionally a single energy source in certain embodiments. In another implementation, the X-ray source 912, which is configured to perform a kV-switching function for emitting X-ray radiation at a predetermined high-level energy and at a predetermined low-level energy. In still another alternative embodiment, the X-ray source 912 is a single source emitting a broad spectrum of X-ray energies. In still another embodiment, the X-ray source 912 is more than a single X-ray emitter and each emitter can emit X-rays separately and emits a different spectrum of X-ray energies.
The detector unit 903 can use energy integrating detectors such as scintillation elements with photo-multiplier tubes or avalanche photo-diodes to detect the resultant scintillation photons from scintillation events resulting from the X-ray radiation interacting with the scintillator elements. The scintillator elements can be crystalline (e.g., NaI(Tl), CsI(Tl), CsI(Na), CsI(pure), CsF, KI(Tl), LiI(Eu), BaF2, CaF2(Eu), ZnS(Ag), CaWO4, CdWO4, YAG(Ce), Y3Al5O12(Ce), GSO, LSO, LaCl3(Ce), LaBr3(Ce), LYSO, BGO, LaCl3(Ce), LaBr3(Ce), C14H10, C14H12, and C10H8), an organic liquid (e.g., an organic solvent with a fluor such as p-terphenyl (C18H14), PBD (C20H14N2O), butyl PBD (C24H22N2O), or PPO (C15H11NO)), a plastic (e.g., a flour suspended in a solid polymer matrix), or other know scintillator.
The photon-counting detectors can use a direct X-ray radiation detectors based on semiconductors, such as cadmium telluride (CdTe), cadmium zinc telluride (CZT), silicon (Si), mercuric iodide (HgI2), and gallium arsenide (GaAs), which have much faster time response than indirect detectors. The fast time response of direct detectors enables them to resolve individual X-ray detection events with only limited pile-up even at the high X-ray fluxes typical of clinical X-ray imaging applications. The amount of energy of the X-ray detected is proportional to the signal generated at the direct detector, and the energies of detection events can be binned into a discrete number of corresponding energy bins yielding spectrally resolved X-ray projection measurements.
The CT scanner also includes a data channel that routes projection measurement results from the photon counting detectors and the detector unit 903 to a data acquisition system 976, a processor 970, memory 978, network controller 980. The data acquisition system 976 controls the acquisition, digitization, and routing of projection data from the detectors. The data acquisition system 976 also includes radiography control circuitry to control the rotation of the annular rotating frames 910 and 930. In one implementation data acquisition system 976 will also control the movement of the bed 916, the operation of the X-ray source 912, and the operation of the X-ray detectors 903. The data acquisition system 976 can be a centralized system or alternatively it can be a distributed system. In an implementation, the data acquisition system 976 is integrated with the processor 970. The processor 970 performs functions including reconstructing images from the projection data, pre-reconstruction processing of the projection data, and post-reconstruction processing of the image data. The pre-reconstruction processing of the projection data can include correcting for detector calibrations, detector nonlinearities, polar effects, noise balancing, and material decomposition. Post-reconstruction processing can include filtering and smoothing the image, volume rendering processing, and image difference processing as needed. The image reconstruction process can be performed using filtered back projection, iterative image reconstruction methods, or stochastic image reconstruction methods. Both the processor 970 and the data acquisition system 976 can make use of the memory 976 to store, e.g., projection data, reconstructed images, calibration data and parameters, and computer programs.
The processor 970 can include a CPU that can be implemented as discrete logic gates, as an Application Specific Integrated Circuit (ASIC), a Field Programmable Gate Array (FPGA) or other Complex Programmable Logic Device (CPLD). An FPGA or CPLD implementation may be coded in VHDL, Verilog, or any other hardware description language and the code may be stored in an electronic memory directly within the FPGA or CPLD, or as a separate electronic memory. Further, the memory may be non-volatile, such as ROM, EPROM, EEPROM or FLASH memory. The memory can also be volatile, such as static or dynamic RAM, and a processor, such as a microcontroller or microprocessor, may be provided to manage the electronic memory as well as the interaction between the FPGA or CPLD and the memory.
Alternatively, the CPU in the reconstruction processor may execute a computer program including a set of computer-readable instructions that perform the functions described herein, the program being stored in any of the above-described non-transitory electronic memories and/or a hard disk drive, CD, DVD, FLASH drive or any other known storage media. Further, the computer-readable instructions may be provided as a utility application, background daemon, or component of an operating system, or combination thereof, executing in conjunction with a processor, such as a Xenon processor from Intel of America or an Opteron processor from AMD of America and an operating system, such as Microsoft VISTA, UNIX, Solaris, LINUX, Apple, MAC-OS and other operating systems known to those skilled in the art. Further, CPU can be implemented as multiple processors cooperatively working in parallel to perform the instructions.
In one implementation, the reconstructed images can be displayed on a display. The display can be an LCD display, CRT display, plasma display, OLED, LED or any other display known in the art.
The memory 978 can be a hard disk drive, CD-ROM drive, DVD drive, FLASH drive, RAM, ROM or any other electronic storage known in the art.
The network controller 980, such as an Intel Ethernet PRO network interface card from Intel Corporation of America, can interface between the various parts of the CT scanner. Additionally, the network controller 980 can also interface with an external network. As can be appreciated, the external network can be a public network, such as the Internet, or a private network such as an LAN or WAN network, or any combination thereof and can also include PSTN or ISDN sub-networks. The external network can also be wired, such as an Ethernet network, or can be wireless such as a cellular network including EDGE, 3G and 4G wireless cellular systems. The wireless network can also be WiFi, Bluetooth, or any other wireless form of communication that is known.
While certain implementations have been described, these implementations have been presented by way of example only, and are not intended to limit the teachings of this disclosure. Indeed, the novel methods, apparatuses and systems described herein may be embodied in a variety of other forms; furthermore, various omissions, substitutions and changes in the form of the methods, apparatuses and systems described herein may be made without departing from the spirit of this disclosure.

## Claims

1. An apparatus to process projection measurements, the apparatus comprising:
processing circuitry configured to
obtain a plurality of datasets, wherein each dataset represents radiation detected by a detector having a plurality of detector elements, the radiation having been transmitted through an object, wherein each dataset corresponds to a respective energy bin having an energy-detection spectrum that is narrower than a bandwidth of an energy spectrum of the transmitted radiation;
calculate a nonlinear detector response of each dataset and subtract the nonlinear detector response from the dataset to obtain a corresponding linear-response dataset, wherein the nonlinear detector response being a response of the detector in which an output of the detector varies nonlinearly with respect to an input of the detector; and
perform a material decomposition of the plurality of linear-response datasets into a plurality of projection lengths of a first material and a second material.

2. The apparatus according claim 1, wherein the processing circuitry is further configured to obtain the plurality of datasets, which represent the detection of X-ray radiation.

3. The apparatus according claim 2, wherein the processing circuitry is further configured to:
calculate the nonlinear detector response using a detector response model, which includes at least one of pileup effects and ballistic deficit effects, polar effects, and characteristics X-ray escape effects.

4. The apparatus according claim 3, wherein the processing circuitry is further configured to:
calculate the nonlinear detector response of each respective detector element of the detector using a corresponding nonlinear response function including a dead time of the corresponding detector element τ, a quadratic response function of the corresponding detector element R1, an X-ray flux at the corresponding detector element, and an incident spectrum that is a spectrum of the energies of the X-rays incident on the corresponding detector element.

5. The apparatus according claim 4, wherein the processing circuitry is further configured to calculate the nonlinear detector response using the X-ray flux at a detector element of the detector that is determined using a reference intensity representing an X-ray intensity at an X-ray source that is related to the X-ray flux at the detector element using a corresponding intensity calibration value A, and the X-ray flux at the detector element and the incident spectrum are each calculated using the plurality of projection lengths, which includes a projection length of the first material and a projection length of the second material.

6. The apparatus according claim 5, wherein the processing circuitry is further configured to:
calibrate the plurality of detector elements by determining calibration values, including the detector response model, using calibration data that includes the dead time τ of each respective detector element of the detector, the quadratic response function R1 of each respective detector element, and the corresponding intensity calibration value A of each respective detector element;
wherein adjusting the calibration values changes a spectrum of the detector-response model, and
the calibration of the detector is performed in an iterative calibration loop by iteratively adjusting the calibration values until either a value returned by a cost function measuring the difference between a spectrum of the plurality of data sets and the spectrum of the detector-response model falls below a predefined threshold or a maximum number of calibration loop iterations is reached.

7. The apparatus according claim 2, wherein the processing circuitry is further configured to:
perform a material decomposition of the plurality of linear-response datasets using noise balancing to
combine the plurality of linear-response datasets into a low-energy projection-data component having a first X-ray energy spectrum, and
combine the plurality of linear-response datasets into a high-energy projection-data component having a second X-ray energy spectrum,
wherein an energy expectation value of the first X-ray energy spectrum is less than an energy expectation value of the second X-ray energy spectrum.

8. The apparatus according claim 7, wherein the processing circuitry is further configured to:
perform the material decomposition of the plurality of linear-response datasets into the plurality of projection lengths by calculating for each detector element of the plurality of detector elements a projection length of the first material and a projection length of the second material, wherein
the calculation of the plurality of projection lengths includes calculating perturbations of beam hardening corrections for variations of attenuation values of the first and second materials around mean attenuation values.

9. The apparatus according claim 2, further comprising:
a polychromatic X-ray source radiating X-rays;
a photon counting spectrally discriminating X-ray detector configured to detect the X-rays that are radiated from the X-ray source, where the X-ray detections are arranged into energy bins according to the detected energy of each X-ray detection, and each energy bin is assigned to a corresponding dataset of the plurality of datasets.

10. The apparatus according claim 2, wherein the processing circuitry is further configured to:
reconstruct a plurality of images using the plurality projection lengths of the first material and using the plurality projection lengths of the second material, where the reconstruction method is one of a filtered back-projection method, an iterative image-reconstruction method, and a statistical image-reconstruction method.

11. The apparatus according claim 2, wherein the processing circuitry is further configured to:
solve for the plurality of projection lengths of the first material and plurality of projection lengths of the first material by performing in an iterative loop the steps of
calculating the nonlinear detector response and subtracting the nonlinear detector response from each respective dataset to obtain the corresponding linear-response dataset, and
then performing the material decomposition on each linear-response dataset;

wherein the plurality of projection lengths of the first material and of the second material resulting from a previous iteration of material decomposition are used as the plurality of first and second material projection lengths for calculating the nonlinear detector response of a next iteration, and the iterative loop stops when a predefined stopping criteria is satisfied.

12. A computed tomography apparatus, comprising:
a radiation source configured to emit radiation toward an image object;
a detector having a plurality of detector elements configured to detect the radiation having been transmitted through the image object and to generate projection data representing an intensity of the radiation at the detector;
a rotation mount configured to rotate the radiation source around the image object, wherein the radiation source is fixedly connected to the rotation mount;
processing circuitry configured to
obtain a plurality of datasets of the projection data, wherein each dataset represents radiation detected by the detector, the radiation having been transmitted through the image object, wherein each dataset corresponds to a respective energy bin having an energy detection spectrum that is narrower than a bandwidth of an energy spectrum of the transmitted radiation;
calculate a nonlinear detector response of each dataset and subtract the nonlinear detector response from the dataset to obtain a corresponding linear-response dataset, wherein the nonlinear detector response being a response of the detector in which an output of the detector varies nonlinearly with respect to an input of the detector; and

perform a material decomposition of the plurality of linear-response datasets into a plurality of projection lengths of a first and a second material.

13. A method of processing X-ray projection data, comprising:
obtaining a plurality of datasets, wherein each dataset represents radiation detected by a detector having a plurality of detector elements, the radiation having been transmitted through an object, and each dataset corresponds to a respective energy bin having an energy detection spectrum that is narrower than a bandwidth of an energy spectrum of the transmitted radiation;
calculating a nonlinear detector response of each dataset and subtracting the nonlinear detector response from the dataset to obtain a corresponding linear-response dataset, wherein the nonlinear detector response being a response of the detector in which an output of the detector varies nonlinearly with respect to an input of the detector;
performing a material decomposition of the plurality of linear-response datasets into a plurality of projection lengths of a first and a second material.

14. The method according to claim 13, wherein the step of calculating a nonlinear detector response uses a detector response model for calculating at least one of the nonlinear detector response including pileup of detection events, ballistic deficit effects, polar effects, and characteristics X-ray escape effects.

15. The method according to claim 13, wherein the step of calculating a nonlinear detector response is calculated using a respective nonlinear response function of each corresponding detector element of the plurality of detector elements, the nonlinear response function including a dead time of the corresponding detector element τ, a quadratic response function of the corresponding detector element R1, an X-ray flux at the corresponding detector element, and an incident spectrum that is a spectrum of the energies of the X-rays incident on the corresponding detector element.

16. The method according to claim 15, further comprising:
calibrating the plurality of detector elements by determining calibration values, including the detector response model, using calibration data that includes the dead time τ of each respective detector element, the quadratic response function R1 of each respective detector element, and the respective intensity calibration value A of each respective detector element;
wherein adjusting the calibration values changes a spectrum of the detector-response model, and
the calibration of the plurality of detector elements is performed in an iterative calibration loop by iteratively adjusting the calibration values until either a value returned by a cost function measuring the difference between a spectrum of the plurality of data sets and the spectrum of the detector-response model falls below a predefined threshold or a maximum number of calibration loop iterations is reached.

17. A non-transitory computer-readable medium storing executable instructions, wherein the instructions, when executed by processing circuitry, cause the processing circuitry to perform the method according to claim 13.