Iterative reconstruction scheme for phase contrast tomography

Abstract

A method of tomographic imaging the real part of the refractive index of the internal structure of an image volume and a corresponding apparatus for performing the method. The method includes a two-step process of first retrieving phase projections of an image volume from projective intensity measurements through the image volume. The second step of the method includes iterative image reconstruction of an image of the image volume using phase projections taken at multiple projection angles. The first step of intensity measurements and subsequent phase retrieval can use one of many possible phase retrieval methods including the Bronnikov phase retrieval method. The second step of iterative image reconstruction can be performed using the total variation minimization method, the algebraic reconstruction technique method, or any other iterative image reconstruction method.

Claims

1. A method of computer aided imaging, comprising: obtaining intensity data representing an intensity of radiation at a detector array, wherein a first dataset of the intensity data represents the intensity of the radiation at the detector array after propagating a first distance past an exit plane of an image volume, a second dataset of the intensity data represents the intensity of the radiation at the detector array after propagating a second distance past the exit plane of the image volume, and the first distance is greater than the second distance; determining, using a phase retrieval method and the intensity data, a phase of the radiation at the exit plane of the image volume; and reconstructing an image of a real part of an index of refraction of the image volume using an iterative reconstruction algorithm, including performing an affine projection of an image vector onto an affine space determined by a plurality of Radon transforms corresponding to the intensity data, and the image vector represents the real part of the index of refraction of the image volume.

2. The method according to claim 1, wherein the radiation is X-ray radiation.

3. The method according to claim 1, wherein the radiation is a collimated beam having a plurality of rays, each of which is substantially parallel with all other of the plurality of rays of the radiation.

4. The method according to claim 1, wherein the radiation is a cone beam having a plurality of rays, each of which substantially divergent with respect to all other of the plurality of rays of the radiation.

5. The method according to claim 1, wherein the radiation is a fan beam.

6. The method according to claim 1, wherein the phase retrieval method is one of a transport intensity equation algorithm, a contrast transfer function algorithm, a mixed transport-intensity-equation/contrast-transfer-function algorithm, a Bronnikov algorithm, a modified Bronnikov algorithm, a phase-attenuation duality algorithm, a single material algorithm, a two materials algorithm, a Fourier method with Born approximation method algorithm, and a Fourier method with Rytov approximation algorithm.

7. The method according to claim 1, wherein the reconstruction algorithm is one of an algebraic reconstruction technique algorithm, a total variation minimization method algorithm, a simultaneous iterative reconstruction algorithm, an iterative least squares technique algorithm, and a simultaneous algebraic reconstruction technique algorithm.

8. An apparatus for determining an image, comprising: a processing circuit configured to obtain intensity data representing an intensity of radiation at a detector array, wherein a first dataset of the intensity data represents the intensity of the radiation at the detector array after propagating a first distance past an exit plane of an image volume, a second dataset of the intensity data represents the intensity of the radiation at the detector array after propagating a second distance past the exit plane of the image volume, and the first distance is greater than the second distance; determine, using a phase retrieval method and the intensity data, a phase of the radiation at the exit plane of the image volume; and reconstruct an image of a real part of an index of refraction of the image volume using an iterative reconstruction algorithm, including performing an affine projection of an image vector onto an affine space determined by a plurality of Radon transforms corresponding to the intensity data, and the image vector represents the real part of the index of refraction of the image volume.

9. The apparatus according to claim 8, wherein the processing circuitry is further configured to reconstruct the image of the real part of the index of refraction using the iterative reconstruction algorithm, wherein the iterative reconstruction algorithm is a total variation minimization algorithm.

10. The apparatus according to claim 8, wherein the processing circuitry is further configured to reconstruct the image of the real part of the index of refraction using the iterative reconstruction algorithm, wherein the iterative reconstruction algorithm is one of an algebraic reconstruction technique algorithm, a simultaneous iterative reconstruction algorithm, an iterative least squares technique algorithm, and a simultaneous algebraic reconstruction technique algorithm.

11. The apparatus according to claim 8, wherein the processing circuitry is further configured to determine the phase of the radiation at the exit plane the image volume, wherein the phase retrieval method is one of a transport intensity equation algorithm, a contrast transfer function algorithm, a mixed transport-intensity-equation/contrast-transfer-function algorithm, a Bronnikov algorithm, a modified Bronnikov algorithm, a phase-attenuation duality algorithm, a single material algorithm, a two materials algorithm, a Fourier method with Born approximation method algorithm, and a Fourier method with Rytov approximation algorithm.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) A more complete appreciation of the invention and many of the attendant advantages thereof will be readily obtained as the same becomes better understood by reference to the following detailed description when considered in connection with the accompanying drawings, wherein:

(2) FIG. 1 shows the arrangement of a projective measurement of an image object onto a detector plan;

(3) FIG. 2 shows a flowchart diagram of a two-step phase contrast tomography method for retrieving phase projections and iterative image reconstruction of the real part of the index of refraction;

(4) FIG. 3 shows the real part of the index of refraction of a an image object;

(5) FIG. 4 shows a diagram of a projective measurement of a an image object;

(6) FIG. 5 shows a reconstruction of the real part of the index of refraction of an image object using the Bronnikov reconstruction formula with 256 projections;

(7) FIG. 6 shows a reconstruction of the real part of the index of refraction of an image object using the Bronnikov reconstruction formula with 32 projections;

(8) FIG. 7 shows a reconstruction of the real part of the index of refraction of an image object using the two-step phase retrieval iterative image reconstruction method with 32 projections;

(9) FIG. 8 shows a flowchart diagram of a generalized two-step phase contrast tomography method for retrieving phase projections and iterative image reconstruction of the real part of the index of refraction; and

(10) FIG. 9 shows an embodiment of computer hardware for performing a two-step phase contrast tomography method.

DETAILED DESCRIPTION

(11) In one embodiment, the present disclosure provides a method for computer aided imaging comprising: obtaining intensity data that indicates an intensity of radiation at a detector array, wherein the radiation has propagated through an image volume, the intensity data is obtained for the radiation propagating in a plurality of projection angles through the image volume, and the intensity data indicates the radiation intensity at a plurality of distances from the image volume; determining a phase of the radiation for each of the plurality of projection angles at a first boundary immediately adjacent to the image volume, wherein the phase is determined using a phase retrieval method; and reconstructing an image of a real part of an index of refraction of the image volume using an iterative reconstruction algorithm, where the iterative reconstruction algorithm includes the step of an performing an affine projection of an image vector onto an radon transform affine space, where the radon transform affine space is determined by a value of the intensity data and a radon transform for a propagation ray of the radiation that corresponds to the value of the intensity data, and the image vector is an approximation of the image of the real part of the index of refraction of the image volume.

(12) In one aspect, the present disclosure provides the radiation is X-ray radiation.

(13) In one aspect, the present disclosure provides that the radiation beam is a colimated beam with parallel rays.

(14) In one aspect, the present disclosure provides that the radiation beam is a divergent beam configured in a cone beam.

(15) In one aspect, the present disclosure provides that the radiation beam is a divergent beam configured in a fan beam.

(16) In one aspect, the present disclosure provides that the phase retrieval method is a method taken from the group consisting of the transport intensity equation algorithm, the contrast transfer function algorithm, the mixed transport-intensity-equation/contrast-transfer-function algorithm, the Bronnikov algorithm, the modified Bronnikov algorithm, the phase-attenuation duality algorithm, the single material algorithm, the two materials algorithm, the Fourier method with Born approximation method algorithm, and the Fourier method with Rytov approximation algorithm.

(17) In one aspect, the present disclosure provides that the reconstruction algorithm is taken from the group consisting of an algebraic reconstruction technique algorithm, a total variation minimization method algorithm, a simultaneous iterative reconstruction algorithm, an iterative least squares technique algorithm, and a simultaneous algebraic reconstruction technique algorithm.

(18) In one aspect, the present disclosure provides that the phase and amplitude of the radiation at the first boundary of the image volume is determined using differential phase contrast computed tomography, where the plurality of rays traverse a series of two diffraction gratings that are separated by a predetermined fraction of the Talbot length.

(19) In one embodiment, the present disclosure provides an apparatus comprising: circuitry configured to intensity data that indicates an intensity of radiation at a detector array, wherein the radiation has propagated through an image volume, the intensity data is obtained for the radiation propagating in a plurality of projection angles through the image volume, and the intensity data indicates the radiation intensity at a plurality of distances from the image volume; circuitry configured to determine a phase of the radiation for each of the plurality of projection angles at a first boundary immediately adjacent to the image volume, wherein the phase is determined using a phase retrieval method; and circuitry configured to reconstruct an image of a real part of an index of refraction of the image volume using an iterative reconstruction algorithm, where the iterative reconstruction algorithm includes the step of an performing an affine projection of an image vector onto an radon transform affine space, where the radon transform affine space is determined by a value of the intensity data and a radon transform for a propagation ray of the radiation that corresponds to the value of the intensity data, and the image vector is an approximation of the image of the real part of the index of refraction of the image volume.

(20) It is well known that diffraction patterns can be categorized into three regions: the near-field region, the Fresnel region, and the Fraunhofer region. Assuming that a mixed phase-and-amplitude object having characteristic length D is positioned at z=0, a uniform wavefront of wavelength propagating in the positive z direction will first pass through the object. After passing through the object, the wavefront evolves as it propagates first through the near-field region (e.g., z<D.sup.2/), then through the Fresnel region (e.g., zD.sup.2/), and finally through the Fraunhofer region (e.g., z>D.sup.2/). Immediately after passing through the object at z=0, only the absorption profile of the object is observed in the intensity profile of the wave. As the wavefront propagates into the near-field region (i.e., z<D.sup.2/, where D is a characteristic length of the image object and is the wavelength of the radiation) interference fringes (also known as diffraction rings) begin to form in regions where there is a sharp change in the phase (i.e., where the second derivative of the phase is non-zero). In the near field region, these interference fringes are sometimes referred to as edge enhancement, and the intensity profile in the near field is given by

(21) $I_{}^{z=d}\left(x^{},y^{}\right)I_{}^{z=0}\left(x^{},y^{}\right)\left(1-\frac{d}{2}{}^2\left(x^{},y^{}\right)\right),$
where (x, y) is the phase, and d is the distance between the image object and the image plane. Thus, in the near field the phase can be found as

(22) $\left(x^{},y^{}\right)-\frac{2}{d}{}^{-2}\left(\frac{I_{}^{z=d}\left(x^{},y^{}\right)}{I_{}^{z=0}\left(x^{},y^{}\right)}-1\right).$
Recognizing that the inverse Laplacian can be expressed in terms of Fourier transforms as

(23) ${}^{-2}=-\frac{1}{4{}^2}{}_2^{-1}\frac{1}{v_{x^{}}^2+v_{y^{}}^2}{}_2,$
where the Laplacian is taken in only the x and y directions, .sub.2.sup.1 and .sub.2 are the inverse two-dimensional Fourier transform and the two-dimensional Fourier transform, respectively, and v.sub.x, and v.sub.y, are the spatial frequencies corresponding respectively to x and y. Therefore, in the near field, where z<<D.sup.2/, the phase of a wavefront can be retrieved from a series of intensity measurements at z=0 and z=d using the expression

(24) $\left(x^{},y^{},{}_i\right)=\frac{1}{2}{}_2^{-1}\left\{\frac{{\hat{g}}_{}\left(v_{x^{}},v_{y^{}}\right)}{v_{x^{}}^2+v_{y^{}}^2}\right\},$
where .sub.(v.sub.x,v.sub.y,)=.sub.2 {g.sub.(x, y)}, and

(25) $g_{}\left(x^{},y^{}\right)=\frac{1}{d}\left(\frac{I_{}^{z=d}\left(x^{},y^{}\right)}{I_{}^{z=0}\left(x^{},y^{}\right)}-1\right).$
The above method of phase retrieval is often referred to as the Bronnikov phase retrieval method. It is also noted that the above expression for the phase retrieval is very similar to the expression inside the integral for the Bronnikov reconstruction formula. Also, it is noted that the outer integral of the Bronnikov reconstruction formula is very similar to a back projection operator, which is a method for solving the inverse Radon transformation problem.

(26) Referring now to the drawings, wherein like reference numerals designate identical or corresponding parts throughout the several views, FIG. 1 shows the geometry of a projective measurement for phase-contrast CT imaging. An image object 212 having a real part of the index of refraction given by the function (X,Y,Z) is positioned with an image volume. A highly coherent X-ray beam propagates along the z-axis, which is at an angle relative to the Z-axis. Multiple projection images are taken at several different projection angles . In each projection image the X-rays pass through the image object and continue to propagate until the X-rays impinge a detector plane 220. The detector plane 220 will include an array of pixels to detect the X-ray intensity at discrete locations across the detector plane 220. These measurements give rise to a map of the intensity I.sub..sup.z=d(x, y) across the detector, which is recorded into computer memory storage and later processed to create a CT image of the image object 212. In FIG. 2 the detector 220 is shown to be a distance d from the closest boundary of the image volume.

(27) In FIG. 1, the detector plane 220 is a two-dimensional plane. In an alternative embodiment, the detector can be one-dimensional detector array 224 and the X-ray beam can be a long and narrow pencil beam. These one-dimensional projections are disadvantageous for phase contrast imaging because one-dimensional projections only provide interference fringes along a single dimension (i.e., the long axis of the detector and beam). Therefore, in order to obtain a three-dimensional phase-contrast image while using one-dimensional projection measurements, for a given projection angle each image pixel in the two-dimensional pixel array should be measured twice, once at each of two orthogonal orientations of the pencil beam (e.g., one measurement where =0 and one measurement where =90, where is the angle of the detector relative to the y-axis of the image plane).

(28) In the interest of simplicity and concreteness, only the case of parallel X-ray beams is discussed here. However, one of ordinary skill in the art would recognize that cone X-ray beams and fan X-ray beams can also be used. Also, X-ray beams that are not perfectly collimated can be used, so long as the X-ray beam has a sufficiently high degree of spatial coherence.

(29) Also, this technique is not limited to X-ray beams. Any type of radiation or other type of wavefront with a sufficiently high degree of spatial coherence can be used to image the real part of the index of refraction of an image object 212 using the proposed imaging method presented here.

(30) Based on the abovementioned similarities between the inner integrals of the Bronnikov reconstruction formula and the near-field phase retrieval and the similarities between the outer integral of the Bronnikov reconstruction formula and the inverse Radon transform via back projections, the proposed algorithm applies a two-step process that is analogous to first performing the inner integrals of the Bronnikov reconstruction formula as the first step and then performing the outer integral of the Bronnikov reconstruction formula Bronnikov reconstruction formula as the second step.

(31) The first step in the two-step process is to retrieve the phase of a radiation wavefront using a phase retrieval method such as the Bronnikov phase retrieval method. This phase retrieval is performed by making a series of intensity measurements at various distances from the image object 212 and then using the Bronnikov phase retrieval method or any other related method for phase retrieval. The Bronnikov phase retrieval method assumes the intensity measurements were performed in the near field, whereas other phase retrieval methods, such as the contrast transfer method, relax this assumption while imposing other assumptions. Depending on the particulars of a given imaging problem, one phase retrieval method may be favored over the rest because the underlying assumptions of the phase retrieval method conform more closely to the given parameters of the imaging problem. A discussion of the underlying assumptions of several phase retrieval methods can be found in A. Buvall, et al., X-ray phase contrast imaging suitable for tomography, Opt. Express 19, 10359 (2011), incorporated herein by reference in its entirety. Specifically, A. Buvall, et al. discusses the Bronnikov phase retrieval method method and six other phase retrieval method including the assumptions for each phase retrieval method. Similarly, M. Langer, et al., Quantitative comparison or direct phase retrieval algorithms in in-line phase tomography, Med. Phys., 35, 4556 (2008), incorporated herein by reference in its entirety, discusses the Bronnikov phase retrieval method and three other phase retrieval method including the assumptions for each phase retrieval method.

(32) The second step in the proposed two-step phase-contrast CT method uses an iterative reconstruction algorithm to obtain an image of the real part of the index of refraction for the image object 212. There are many iterative CT reconstruction algorithms and any of these algorithms can be used depending on such factors as a priori knowledge about the image object and knowledge of the projective imaging system. For example, the total variation (TV) minimization algorithm is popular for imaging the soft matter commonly encountered in various organs. This algorithm assumes a priori knowledge that the tomographic images are relatively constant over extended volumes (e.g., the region occupied by an organ), and rapid variations in the image occur primarily at the boundaries of internal structures (e.g., at the boundaries between organs). Therefore the TV-minimization algorithm seeks to impose conditions that reinforce this structure by minimizing the l.sub.1-norm on the gradient of the tomographic images. The steps of the TV-minimization algorithm are mostly identical to the ART algorithm discussed above, except that the regularization in the TV-minimization algorithm minimizes the l.sub.1-norm on the gradient of the tomographic images. Details of the TV-minimization algorithm can be found in E. Sidky, et al., Accurate image reconstruction from few-views and limited-angle data in divergent beam CT, J Xray Sci Technol, 14, 119 (2008), incorporated herein in its entirety by reference.

(33) One advantage of the TV minimization algorithm is that it has been demonstrated to produce high quality tomographic images with much fewer projection images compared with filtered back projection image reconstructions.

(34) As an alternative embodiment, one of ordinary skill in the art will recognize that there are other regularization methods besides those already mentioned of imposing non-negativity on the tomographic images and the regularization method that minimizes the l.sub.1-norm of the gradient of the tomographic images. Thus, one of ordinary skill in the art will recognize that other regularization methods can be used in the iterative tomographic image reconstruction.

(35) FIG. 2 shows a flowchart of one embodiment of the proposed phase-contrast CT imaging method. In this embodiment the method is subdivided into two partsthe phase projection loop 310 and the iterative tomographic image reconstruction loop 330.

(36) The phase projection loop 310 acquires phase projections of the image object 212 at projection angles =.sub.i where i=1 . . . N.sub. for a total of N.sub. unique phase projections.

(37) The first step 314 in the phase projection loop 310 increments the loop variable i and sets the current projection angle to the angle .sub.i.

(38) The second step 316 in the phase projection loop 310 consists of measuring the intensity projections onto detector plane 220 located at z=0 and onto the detector plane 220 when it is located at z=d. If the image object 212 is a pure phase object (i.e., the image object 212 does not absorb X-rays), then the intensity measured at z=0 will always be identical to the X-ray beam intensity in the absence of an image object. In this case the intensity measurements simplify to a single measurement of I.sub..sup.z=d(x, y), and the measurement of I.sub..sup.z=0(x, y) can be obtained from a previous calibration measurement taken in the absence of the image object 212.

(39) The third step 318 in the phase projection loop 310 consists of calculating the quantity .sub.(v.sub.x,v.sub.y,) from the measured intensities I.sub..sup.z=d(x, y) and I.sub..sup.z=0(x, y), where

(40) $g_{}\left(x^{},y^{}\right)=\frac{1}{d}{}_2\left\{\frac{I_{}^{z=d}\left(x^{},y^{}\right)}{I_{}^{z=0}\left(x^{},y^{}\right)}-1\right\}.$

(41) The fourth step 320 in the phase projection loop 310 consists of calculating the phase projection (x, y, .sub.i) from .sub.(v.sub.x,v.sub.y,), where

(42) $\left(x^{},y^{},{}_i\right)=\frac{1}{2}{}_2^{-1}\left\{\frac{{\hat{g}}_{}\left(v_{x^{}},v_{y^{}}\right)}{v_{x^{}}^2+v_{y^{}}^2}\right\}.$

(43) The phase projection loop 310 is repeated N.sub. times in order to obtain N.sub. unique phase projections at different projection angles.

(44) After completing the phase projection loop 310, the iterative tomographic image reconstruction loop 330 is performed. All of the phase projection values (x,y,.sub.i) are arranged into a vector {right arrow over (g)} with components g.sub.n and a matrix M with components M.sub.n,m is defined such that each row vector {right arrow over (M)}.sub.n is the discretized radon transform corresponding to the projection value g.sub.n.

(45) The first step 332 of the image reconstruction loop 330 is to initialize a tomographic image vector {right arrow over ()}.sub.start and the intermediate tomographic image vector {right arrow over ()}.sub.Inter.sup.n=0 to zero. Alternatively, these vectors can be initialized to any other value including using the filtered back projection method, the Bronnikov analytical expression given above, or some other method to estimate an image of the image object 212.

(46) The second step 334 of the image reconstruction loop 330 performs affine projections onto the rows of the discretized radon transform matrix M, where the input vector for the current affine projection {right arrow over ()}.sub.i results from the projection of {right arrow over ()}.sub.i-1 onto the affine space corresponding to {right arrow over (M)}.sub.i-1 and g.sub.i-1, where
{right arrow over ()}.sub.i={right arrow over ()}.sub.i-1{right arrow over (M)}.sub.i-1(g.sub.i-1{right arrow over (M)}.sub.i-1.Math.{right arrow over ()}.sub.i-1/{right arrow over (M)}.sub.i-1.Math.{right arrow over (M)}.sub.i-1).

(47) The third step 336 of the image reconstruction loop 330 performs the chosen regularization method (e.g., the TV minimization algorithm) on the current intermediate tomographic image vector {right arrow over ()}.sub.Inter.sup.n, where {right arrow over ()}.sub.Inter.sup.n={right arrow over ()}.sub.i=Ndata. The value of the start tomographic image vector {right arrow over ()}.sub.Start for the next loop is then set to the result of the regularization method.

(48) The loop then repeats until the image estimate {right arrow over ()}.sub.Inter.sup.n satisfies a predefined convergence criteria, e.g., when difference between the current tomographic image {right arrow over ()}.sub.Inter.sup.n and previous tomographic image {right arrow over ()}.sub.Inter.sup.n-1 reduces below a predefined threshold.

(49) The final step 350 of the two-step phase contrast imaging method is setting the final tomographic image {right arrow over ()}.sub.Out to be equal to the intermediate tomographic image {right arrow over ()}.sub.Inter.sup.n from the final iteration of the image reconstruction loop 330.

(50) FIGS. 3-7 show an example of how the proposed method using the TV minimization algorithm compares with the Bronnikov reconstruction formula for image reconstruction, and especially how the two methods compare for the case where there are only a few projection views. The few-view scenario is important because fewer views can result in a lower dosage of radiation to the patient.

(51) FIG. 3 shows an example of an image object, where grey indicates an index of refraction of 1, white indicates an index of refraction greater than 1 and black indicates an index of refraction less than one. The index of refraction is the real part of the index of refraction, unless explicitly stated otherwise. In FIGS. 3-7, the imaginary part of the index of refraction is assumed to be zero. Thus, this example considers only the case of a pure phase object.

(52) FIG. 4 shows the arrangement for a projection image. The X-rays propagate from the source through the image object until the X-rays are absorbed by the detector array at z=d. The intensity measurement is registered on the detector array, and the intensity measurement is then read and recorded into a computer memory where it is stored for later processing.

(53) FIG. 5 shows the reconstructed image of the index of refraction using the Bronnikov reconstruction formula with 256 projections. Comparing FIG. 5 with FIG. 3, it can be seen that the reconstructed image is consistent with the image object. However, as the number of projections decreases the reconstructed image becomes poorer.

(54) FIG. 6 shows the reconstructed image of the index of refraction using the Bronnikov reconstruction formula with 32 projections. Streaks are present in the reconstructed image, similar to those streaks found in more conventional absorption-based CT for the case of filtered back projections using only a few projection views.

(55) FIG. 7 shows the reconstructed image created by using the proposed method with the TV minimization algorithm and only 32 projections. Comparing FIG. 7 with FIG. 3 and with FIG. 5, it is observed that, with far fewer projections, the disclosed method achieves similar image quality to a Bronnikov reconstruction formula that uses almost 10 times the number of projections.

(56) In FIGS. 3-7, the Bronnikov phase retrieval method (also called the Bronnikov phase retrieval algorithm) has been used for the phase retrieval step in the two-step phase contrast CT method. In alternative embodiments, other phase retrieval methods can be used. In many imaging scenarios, these alternative phase retrieval methods may be preferable to the Bronnikov phase retrieval method depending on the parameters of the imaging scenario. Conventional phase retrieval methods include the transport intensity equation algorithm, the contrast transfer function algorithm, the mixed transport-intensity-equation/contrast-transfer-function algorithm, the Bronnikov algorithm, the modified Bronnikov algorithm, the phase-attenuation duality algorithm, the single material algorithm, the two materials algorithm, the Fourier method with Born approximation method algorithm, and the Fourier method with Rytov approximation algorithm. In many of these phase retrieval methods, the intensity profile has to be measured at multiple distances from the image object 212, where the image positions are different than the z=d and z=0 image positions used in the Bronnikov phase retrieval algorithm.

(57) FIG. 8 shows a flowchart for an alternative embodiment of the disclosed inline phase contrast CT method, where the phase retrieval step 914 is generalized to be any of the above-mentioned phase retrieval methods, and the intensity measurement step 912 is generalized to include the measurements required in order to perform the selected phase retrieval method.

(58) Alternative embodiments of the proposed method can use other regularization terms/methods and/or variations on the above mentioned iterative CT reconstruction methods. Some examples of iterative reconstruction methods that can be applied in place of the TV-minimization algorithm include the algebraic reconstruction technique algorithm, the simultaneous iterative reconstruction algorithm, the iterative least squares technique algorithm, and the simultaneous algebraic reconstruction technique algorithm.

(59) In an alternative embodiment, the proposed method can be modified to work for differential phase contrast CT. As described above, in differential phase contrast CT the phase information is derived from measurements of interference patterns created when the radiation beam propagates through multiple diffraction gratings that have been positioned in the X-ray beam path after the image object 212. In this case the differential phase contrast CT measurements and data processing provide the phase information directly, thus eliminating the necessity of a separate phase retrieval step. Therefore, the disclosed two-step phase contrast CT method reduces to a one-step phase contrast CT method when applied to differential phase contrast imaging. In the differential phase contrast adaptation of the proposed phase contrast CT method, the phase projections are obtained directly from the differential phase contrast measurements and data processing. Given the phase projections, the index of refraction image of the image object 212 is created by the iterative image reconstruction loop 330.

(60) Next, a hardware description according to exemplary embodiments is described with reference to FIG. 9. In FIG. 9, the computational hardware includes a CPU 1000 which performs the disclosed two-step phase contrast CT method described above. The process data and instructions may be stored in memory 1002. These processes and instructions may also be stored on a storage medium disk 1004 such as a hard drive (HDD) or portable storage medium or may be stored remotely. Further, the claimed advancements are not limited by the form of the computer-readable media on which the instructions of the inventive process are stored. For example, the instructions may be stored on CDs, DVDs, in FLASH memory, RAM, ROM, PROM, EPROM, EEPROM, hard disk or any other information processing device with which the CT processing and reconstruction hardware communicates, such as a server or computer.

(61) Further, the claimed advancements may be provided as a utility application, background daemon, or component of an operating system, or combination thereof, executing in conjunction with CPU 1000 and an operating system such as Microsoft Windows 7, UNIX, Solaris, LINUX, Apple MAC-OS and other systems known to those skilled in the art.

(62) CPU 1000 may be a Xenon or Core processor from Intel of America or an Opteron processor from AMD of America, or may be other processor types that would be recognized by one of ordinary skill in the art. Alternatively, the CPU 1000 may be implemented on an FPGA, ASIC, PLD or using discrete logic circuits, as one of ordinary skill in the art would recognize. Further, CPU 1000 may be implemented as multiple processors cooperatively working in parallel to perform the instructions of the inventive processes described above.

(63) The device in FIG. 9 also includes a network controller 1006, such as an Intel Ethernet PRO network interface card from Intel Corporation of America, for interfacing with a network 1030. As can be appreciated, the network 1030 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 network 1030 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 network 1030 can also be WiFi, Bluetooth, or any other wireless form of communication that is known.

(64) The device further includes a display controller 1008, such as a NVIDIA GeForce GTX or Quadro graphics adaptor from NVIDIA Corporation of America for interfacing with display 1010, such as a Hewlett Packard HPL2445w LCD monitor. A general purpose I/O interface 1012 interfaces with a keyboard and/or mouse 1014 as well as a touch screen panel 1016 on or separate from display 1010. General purpose I/O interface also connects to a variety of peripherals 1018 including printers and scanners, such as an OfficeJet or DeskJet from Hewlett Packard.

(65) A sound controller 1020 is also provided in the device, such as Sound Blaster X-Fi Titanium from Creative, to interface with speakers/microphone 1022 thereby providing sounds and/or music.

(66) The general purpose storage controller 1024 connects the storage medium disk 1004 with communication bus 1026, which may be an ISA, EISA, VESA, PCI, or similar, for interconnecting all of the components of the device. A description of the general features and functionality of the display 1010, keyboard and/or mouse 1014, as well as the display controller 1008, storage controller 1024, network controller 1006, sound controller 1020, and general purpose I/O interface 1012 is omitted herein for brevity as these features are known.

(67) While certain embodiments have been described, these embodiments have been presented by way of example only, and are not intended to limit the scope of the inventions. 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 the inventions. The accompanying claims and their equivalents are intended to cover such forms or modifications as would fall within the scope and spirit of the inventions.