Digital x-ray tomosynthesis system

• US 7,177,390 B2
• Filed: 03/10/2005
• Issued: 02/13/2007
• Est. Priority Date: 03/11/2004
• Status: Expired due to Fees

• Current Assignee: Trex Enterprises Corporation
• Original Assignee: Trex Enterprises Corporation
• Inventors: Martin, Peter, Spivey, Brett
• Primary Examiner(s): Glick; Edward J.
• Assistant Examiner(s): SONG, HOON K
• Application Number:
  US11/077,666
• Publication Number:
  US 20050226369A1
• Time in Patent Office: 705 Days
• Field of Search: 378/21, 378/22, 378/23, 378/25, 378/901, 378/26
• US Class Current: 378/25
• CPC Class Codes: A61B 6/025 Tomosynthesis
• CPC Class Codes: A61B 6/4405 the apparatus being movable or portable, e.g. handheld or mounted on a trolley
• CPC Class Codes: A61B 6/4441 the rigid structure being a C-arm or U-arm
• CPC Class Codes: A61B 6/466 adapted to display 3D data
• CPC Class Codes: A61B 6/587 Alignment of source unit to detector unit
• CPC Class Codes: A61B 6/588 Setting distance between source unit and detector unit
• CPC Class Codes: G01N 2223/419 computed tomograph
• CPC Class Codes: G01N 23/044 using laminography or tomosynthesis
• CPC Class Codes: G01N 23/046 using tomography, e.g. computed tomography [CT]

Abstract

Method and device for digital x-ray tomosynthesis. Tomographic and/or three-dimensional images of an object are obtained with an x-ray source and a digital x-ray image sensor. The source, object and sensor are positioned relative to each other and attenuation data is obtained for a large number of rays of x-radiation through the object. A special algorithm is provided to convert the data into images. To calculate the images the algorithm uses iterative processes with a least squares type technique but with generalized (as opposed to specific) functions. The algorithm solves for the functions which are the images. Preferred embodiments include a system having an x-ray point source with a cone of diverging x-rays, a two-dimensional digital x-ray image sensor, two linear translation stages to independently move both the x-ray source and the digital x-ray image sensor, two rotation mechanisms to rotate the two linear translation stages, a microprocessor to control the data acquisition, and a computer programmed with a special algorithm to calculate the tomographic images. A plurality of sets of digital data (representing x-ray algorithm images of an object) are acquired by the digital x-ray image sensor, with the x-ray source and the digital x-ray image sensor located at different positions and angles relative to the object. The digital data representing the x-ray attenuation images is stored in the computer. Special mathematical algorithms then compute multiple images of the object using the acquired digital data. These images could include multiple tomographic images, a three-dimensional image, or a multiple three-dimensional images.

Claims

1. A digital x-ray tomosynthesis system comprising:

A) at least one x-ray source,

B) a two-dimensional digital x-ray image sensor,

C) an x-ray data collection-positioning mechanism for positioning, relative to each other, said at least one x-ray source, said image sensor and an object for collecting with said image sensor x-ray attenuation data representing attenuation of a large number of rays of x-radiation from said at least one source through said object to said image sensor, and

D) a computer processor programmed with an algorithm for producing tomographic or three-dimensional images of said object or portions of said object using a least squares technique utilizing generalized non-specific functions;

wherein in said algorithm a three-dimensional object being imaged is represented by a scalar function d(x_tomo,y_tomo,h) in the dimension of Hounsfield units with one Hounsfield unit equal to a 0.1 % difference in the density of water and the range of Hounsfield units is from −

1000 (air) to 1000 (bone) with 0 Hounsfields as the density of water, wherein the attenuation of x-rays directed along a line (trajectory s) is given by

$A_i<br><br>\left(s\right)=A_i<br><br>\left[{\int <br><br>}_0^s<br><br>d<br><br>\left(x_{\mathrm{tomo}},y_{\mathrm{tomo}},h\right)<br><br><br><br>d<br><br>s^{\mathrm{\prime <br><br>}}\right]$

where (x_tomo,y_tomo) are the coordinates where trajectories cross the tomographic plane at z=h.

2. The system as in claim 1 wherein the attenuation of x-rays are expressed in logarithmic units.

3. The system as in claim 2 wherein logarithms of pixelated digital x-ray images acquired by the digital x-ray image sensor are represented by

m_i(x,y)=m0_i(x,y)+n_i(x,y)

where

m0_i(x,y)=w(x,y)A_i(s){circle around (× )}psƒ

(x,y)

are noise-free measurements for each N source positions (X_i,Y_i);

i=1,N, n_i(x,y) is noise, w_i(x,y) is the white field image acquired during calibration, and psƒ

(x,y) is a point spread function.

4. The system as in claim 3 wherein m0_i(x,y) is also represented by:

m0_i(x,y)=w(x,y)A

${\mathrm{m0}}_i<br><br>\left(x,y\right)=w<br><br>\left(x,y\right)<br><br>A<br><br>\left[\sec <br><br><br><br>{\mathrm{\theta <br><br>}}_i<br><br>\left(x,y\right)<br><br>{\int <br><br>}_0^{\mathrm{SID}}<br><br><br><br>d<br><br>\left[x<br><br>\left(1-\frac{h}{\mathrm{SID}}\right)+\frac{{\mathrm{hX}}_i}{\mathrm{SID}},y<br><br>\left(1-\frac{h}{\mathrm{SID}}\right)+\frac{{\mathrm{hY}}_i}{\mathrm{SID}},h\right]<br><br>d<br><br>h\right]\otimes <br><br>\mathrm{psf}<br><br>\left(x,y\right).$

5. The system as in claim 4 wherein an equation for m0_i(x,y) is inverted to express a tomographic slice image data d(x,y,h) in terms of the acquired noise free data m0(x,y) by minimizing the following generalized chi-squared function:

${\mathrm{\chi <br><br>}}^2=\underset{i}{\sum <br><br>}<br><br><br><br>\int <br><br>\int <br><br>{\left(m_i<br><br>\left(x_0,y_0\right)-{\mathrm{m0}}_i<br><br>\left(x_0,y_0\right)\right)<br><br>\left[{\left(n<br><br>\left(x_0,y_0\right)<br><br>n<br><br>\left(x_1,y_1\right)\right)}_{\mathrm{avg}}\right]}^{-1}<br><br>\left(m_i<br><br>\left(x_1,y_1\right)-\mathrm{m0}<br><br>\left(x_1,y_1\right)\right)<br><br>d<br><br>x_0<br><br>d<br><br>y_0<br><br>d<br><br>x_1<br><br>d<br><br>y_1+\mathrm{\varepsilon <br><br>}<br><br>\int <br><br>\int <br><br>d<br><br>{\left(x,y,h\right)}^2<br><br>d<br><br>x<br><br>d<br><br>y<br><br>d<br><br>h$

versus a function d(x_tomo,y_tomo,h),

where e is a regularization parameter chosen so that the first term averages 1 per measurement and m0 is implicitly a function of d(x_tomo,y_tomo,h).

6. The system as in claim 5 wherein the equation set forth in claim 4 is simplified to produce an equation:

$\mathrm{m0}<br><br>\left(x,y\right)=w_i<br><br>\left(x,y\right)<br><br>A<br><br>\left[\sec <br><br><br><br>{\mathrm{\theta <br><br>}}_i<br><br>\left(x,y\right)<br><br>{\int <br><br>}_0^{\mathrm{\infty <br><br>}}<br><br>D<br><br>\left(x+\frac{X_i<br><br>z}{\mathrm{SID}}<br><br>,y+\frac{Y_i<br><br>z}{\mathrm{SID}},z\right)<br><br>d<br><br>z\right]\otimes <br><br>\mathrm{psf}<br><br>\left(x,y\right)$

and equation set forth in claim 2 is expressed as;

${\mathrm{\chi <br><br>}}^2=\underset{i}{\sum <br><br>}<br><br><br><br>\int <br><br>\int <br><br>{\left(m_i<br><br>\left(x_0,y_0\right)-{\mathrm{m0}}_i<br><br>\left(x_0,y_0\right)\right)<br><br>\left[{\left(n<br><br>\left(x_0,y_0\right)<br><br>n<br><br>\left(x_1,y_1\right)\right)}_{\mathrm{avg}}\right]}^{-1}<br><br>\left(m_i<br><br>\left(x_1,y_1\right)-\mathrm{m0}<br><br>\left(x_1,y_1\right)\right)<br><br>d<br><br>x_0<br><br>d<br><br>y_0<br><br>d<br><br>x_1<br><br>d<br><br>y_1+\mathrm{\varepsilon <br><br>}<br><br>\int <br><br>\int <br><br>{D<br><br>\left(x,y,z\right)}^2<br><br>d<br><br>x<br><br>d<br><br>y<br><br>d<br><br>h.$

7. The system as in claim 6 wherein a noise correlation is a noise spectrum function (NPS) using the following equation:

(n_i(x₀,y₀)n_i(x₁,y₁))_avg=∫

NPS(x₀−

x,y₀−

y,x₁−

x,y₁−

y,x,y)dxd

NPS is approximated by;

∫

NPS(x₀−

x,y₀−

y,x₁−

x,y₁−

y,x,y)dxdy=nps(x₀−

x₁,y₀−

y₁)√

{square root over (noise₁(x₀,y₀))}√

{square root over (noise_i(x₁,y₁))}

where NPS is the noise power spectrum correlation function normalized to 1 at 0 lp/mm.

8. The system as in claim 7 an “

integrated Hounsfield”

measurement M and its noise free version M0 are is defined according to the equations set forth below;

${\mathrm{M0}}_i<br><br>\left(x,y\right)=\int <br><br>D<br><br>\left(x+\mathrm{\alpha <br><br>}<br><br><br><br>z,y+\mathrm{\beta <br><br>z},z\right)<br><br>d<br><br>z=\frac{1}{\sec <br><br><br><br>{\mathrm{\theta <br><br>}}_i<br><br>\left(x,y\right)}<br><br>A^{-1}[\left(\frac{{\mathrm{m0}}_i<br><br>\left(x,y\right)}{w_i<br><br>\left(x,y\right)}\right)\otimes <br><br>\left({\mathrm{psf}}^{-1}<br><br>\left(x,y\right)\right)<br><br>M_i<br><br>\left(x,y\right)=\frac{1}{\sec <br><br><br><br>{\mathrm{\theta <br><br>}}_i<br><br>\left(x,y\right)}<br><br>A<br><br>\left[\left(\frac{m_i<br><br>\left(x,y\right)}{w_i<br><br>\left(x,y\right)}\right)\otimes <br><br>\left[{\mathrm{psf}}^{-1}<br><br>\left(x,y\right)\right]\right]$

where Mo_i is a transform of raw measurements m0.

9. The system as in claim 8 wherein a first order Taylor expansion is used to bring a slowly varying white field term w_i(x,y) outside of the convolution, where the difference of M and M0, i.e.

$M_i$

( x , y ) - M0 i ⁡

( x , y ) = ( m i ⁡

( x , y ) - m0 i ⁡

( x , y ) ) ⊗

( psf - 1 ⁡

( x , y ) ) sec ⁢

⁢

θ

i ⁡

( x , y ) ⁢

w i ⁡

( x , y ) ⁢

A ′

⁡

( M i ⁡

( x , y ) ⁢

sec ⁢

⁢

θ

i ⁡

( x , y ) )

is a transformed noise term and

where

$A^{\mathrm{\prime <br><br>}}<br><br>\left(v\right)=\frac{d<br><br>}{d<br><br>v}<br><br>A<br><br>\left(v\right).$

10. The system as in claim 9 wherein w_i(x,y)A′

(M_i(x,y)secθ

_i(x,y)) is the derivative of a signal versus thickness defining “

dsignal”

which can be calculated or measured and utilized to calculate;

m_i(x,y)−

m0_i(x,y)=[(M0_i(x,y)−

M_i(x,y))dsignal_i(x,y)secθ

_i(x,y)]{circle around (×

)}psƒ

(x,y)

so that the optimization function set for the in claim 2 is expressed as;

${\mathrm{\chi <br><br>}}^2={\mathrm{\Sigma <br><br>}<br><br>\left[\left(\left(M_i-{\mathrm{M0}}_i\right)<br><br>\frac{{\mathrm{dsignal}}_i<br><br>\sec <br><br><br><br>{\mathrm{\theta <br><br>}}_i}{\sqrt{{\mathrm{noise}}_i}}\right)\otimes <br><br>\mathrm{psf}\otimes <br><br>\left({\mathrm{nps}}^{-1}\right)\otimes <br><br><br><br>\mathrm{mtf}\otimes <br><br>\left[\left(M_i-{\mathrm{M0}}_i\right)<br><br>\frac{{\mathrm{dsignal}}_i<br><br>\sec <br><br><br><br>{\mathrm{\theta <br><br>}}_i}{\sqrt{{\mathrm{noise}}_i}}\right]\right]}_{x=0<br><br>y=0}+\mathrm{\varepsilon <br><br>}<br><br>\int <br><br>\int <br><br>{D<br><br>\left(x,y,z\right)}^2<br><br>d<br><br>x<br><br>d<br><br>y<br><br>d<br><br>z.$

11. The system as in claim 10 wherein the convolution in the middle of the last equation of claim 10 defines a detector quantum efficiency

dqe=psƒ

(nps<sup>−

1</sup>){circle around (×

)}psƒ

that defines a ratio

$\frac{\mathrm{noise}}{{\mathrm{dsignal}}^2}=\mathrm{\Delta <br><br>}<br><br><br><br>\mathrm{t2}$

which represents a “

thickness noise”

squared and Δ

t2 is treated as a function only of thickness and can be calculated or measured.

12. The system as in claim 11 wherein a noise function

${\sigma }_{}$

i ⁡

( x , y ) = noise i ⁡

( x , y ) dsignal i ⁡

( x , y ) ⁢

sec ⁢

⁢

θ

i ⁡

( x , y ) - 1 = sec ⁢

⁢

θ

i ⁡

( x , y ) - 1 ⁢

Δ

⁢

⁢

t2 ⁡

( M i ⁡

( x , y ) ⁢

sec ⁢

⁢

θ

i ⁡

( x , y ) )

and the last equation set forth in claim 10 is expressed as

${\mathrm{\chi <br><br>}}^2=\underset{i}{\sum <br><br>}<br><br><br><br>{\left[\left(\frac{M_i-\int <br><br>D<br><br>\left(x+{\mathrm{\alpha <br><br>}}_i<br><br>z,y+{\mathrm{\beta <br><br>}}_i<br><br>z,z\right)}{{\mathrm{\sigma <br><br>}}_i<br><br>\left(x,y\right)}\right)\otimes <br><br><br><br><br><br>\mathrm{dqe}\otimes <br><br>\left(\frac{M_i-\int <br><br>D<br><br>\left(x+{\mathrm{\alpha <br><br>}}_i<br><br>z,y+{\mathrm{\beta <br><br>}}_i<br><br>z,z\right)}{{\mathrm{\sigma <br><br>}}_i<br><br>\left(x,y\right)}\right)\right]}_{x=0,y=0}+<br><br><br><br>\mathrm{\varepsilon <br><br>}<br><br>\int <br><br>\int <br><br>{D<br><br>\left(x,y,z\right)}^2<br><br>d<br><br>x<br><br>d<br><br>y<br><br>d<br><br>z$

which is solved by calculating the derivative of χ

² versus D(x₀,y₀,z₀) and calculating where this derivative equals 0 to produce;

$0=\underset{i}{\sum <br><br>}<br><br><br><br>{\left[\left(\left(\frac{\int <br><br>D<br><br>\left(x+{\mathrm{\alpha <br><br>}}_i<br><br>z,y+{\mathrm{\beta <br><br>}}_i<br><br>z,z\right)<br><br>d<br><br>z}{{\mathrm{\sigma <br><br>}}_i<br><br>\left(x,y\right)}\right)\otimes <br><br><br><br><br><br>\mathrm{dqe}\right)<br><br>\frac{1}{{\mathrm{\sigma <br><br>}}_i<br><br>\left(x,y\right)}\right]}_{x=x_o-{\mathrm{\alpha <br><br>}}_i<br><br>z,y=y_o-{\mathrm{\beta <br><br>}}_i<br><br>z}+\mathrm{\varepsilon <br><br>}<br><br><br><br>D<br><br>\left(x_o,y_o,z\right).$

13. The system as in claim 12 where D is defined from a generator that is discrete, i.e.:

$D<br><br>\left(x,y,z\right)=\underset{j}{\sum <br><br>}<br><br><br><br>G<br><br>\left(x-{\mathrm{\alpha <br><br>}}_j<br><br>z,y-{\mathrm{\beta <br><br>}}_j<br><br>z,z\right)$

so that the equation set forth in claim 11 becomes;

$0=\underset{i}{\sum <br><br>}<br><br><br><br>{\left[\left(\left(\frac{\int <br><br>\underset{j}{\sum <br><br>}<br><br><br><br>G<br><br>\left(x-{\mathrm{\alpha <br><br>}}_j<br><br>z+{\mathrm{\alpha <br><br>}}_i<br><br>z,y-{\mathrm{\beta <br><br>}}_j<br><br>z+{\mathrm{\beta <br><br>}}_i<br><br>z\right)<br><br>d<br><br>z-M_i}{{\mathrm{\sigma <br><br>}}_i<br><br>\left(x,y\right)}\right)\otimes <br><br><br><br><br><br>\mathrm{dge}\right)<br><br>\frac{1}{{\mathrm{\sigma <br><br>}}_i<br><br>\left(x,y\right)}\right]}_{x=x_o-{\mathrm{\alpha <br><br>}}_i<br><br>z,y=y_o-{\mathrm{\beta <br><br>}}_i<br><br>z}+\mathrm{\varepsilon <br><br>}<br><br>\underset{i}{\sum <br><br>}<br><br><br><br>G<br><br>\left(x_o-{\mathrm{\alpha <br><br>}}_i<br><br>z_o,y_o-{\mathrm{\beta <br><br>}}_i<br><br>z_o\right).$

14. The system as in claim 13 wherein the i-summation are identically assumed to be equal to zero, so that:

$0=\left(\left(\frac{\int <br><br>\underset{j}{\sum <br><br>}<br><br><br><br>G<br><br>\left(x-{\mathrm{\alpha <br><br>}}_j<br><br>z+{\mathrm{\alpha <br><br>}}_i<br><br>z,y-{\mathrm{\beta <br><br>}}_j<br><br>z+{\mathrm{\beta <br><br>}}_i<br><br>z\right)<br><br><br><br>d<br><br>z-M_i}{{\mathrm{\sigma <br><br>}}_i<br><br>\left(x,y\right)}\right)\otimes <br><br>\mathrm{dge}\right)<br><br>\frac{1}{{\mathrm{\sigma <br><br>}}_i<br><br>\left(x,\right)}+\mathrm{\varepsilon <br><br>}<br><br><br><br>G_i<br><br>\left(x,y\right).$

15. The system as in claim 14 wherein the equation in claim 14 is simplified by defining a function C as

$C_{\mathrm{ij}}$

( x , y ) = ∫

0 Z max ⁢

δ

⁡

[ x + ( α

i - α

j ) ⁢

z , y + ( β

i - β

j ) ⁢

z ] ⁢

⁢

ⅆ

z

wherein z_max in claim 14 equation is determined by the height of an object that is imaged, and the integral excludes absorption below the detector or above z_max in the solution so that the last equation of claim 14 is expressed as;

$\mathrm{dge}\otimes <br><br>\left(\frac{M_i}{{\mathrm{\sigma <br><br>}}_i}\right)=\mathrm{dge}\otimes <br><br>\left(\frac{C_{\mathrm{ij}}\otimes <br><br>G_j}{{\mathrm{\sigma <br><br>}}_i}\right)+{\mathrm{\varepsilon <br><br>\sigma <br><br>}}_i<br><br>G_j<br><br><br><br>\left(\mathrm{sum}<br><br><br><br>\mathrm{over}<br><br><br><br>j\right)$

which is solved for G, then D is solved from G and d is solved from D.

16. The system as in claim 15 wherein s in the equation set forth in claim 15 is assumed to be constant, and the equation is inverted using Fourier transforms wherein the convolutions become products

ℑ

(dqe)ℑ

(M_i)=ℑ

(dqe)ℑ

(C_ij)ℑ

(G_j)+ε

σ

_i²ℑ

(G_i)

where ℑ

denotes a Fourier transform so that the equation is then expressed as

ℑ

(G₁)=(ℑ

(dqe)ℑ

(C_ij)+ε

σ

_i²δ

_ij)<sup>−

1</sup>ℑ

(dqe)ℑ

(M_i)

and the inverses are performed individually over each spatial frequency.

17. The system as in claim 16 wherein inverse matrices are computed once and stored as a look-up table for improved computational speed and stored data is required for each noise level and each value of z_max.

18. The system as in claim 15 wherein a non-constant s is assumed and multiple iterations are utilized with a constant σ

_trial, solving equation the last equation in claim 12 using a Fourier transform method, calculating an error term, then iterating by reconstructing for the error term but using successively different values for the trial σ

value σ

_trial, wherein an the error term is calculated as follows

$\mathrm{error}=\mathrm{dge}\otimes <br><br>\left(\frac{M_i}{{\mathrm{\sigma <br><br>}}_i}\right)-\mathrm{dge}\otimes <br><br>\left(\frac{C_{\mathrm{ij}}\otimes <br><br>G_j}{{\mathrm{\sigma <br><br>}}_i}\right)+\mathrm{\varepsilon <br><br>}<br><br><br><br>{\mathrm{\sigma <br><br>}}_i<br><br>G_j$

the error term is inverted using a constant σ

dqe{circle around (×

)}(C_ij{circle around (×

)}Δ

G_i)+ε

σ

_trial²Δ

G_i=error,

then a reconstructed generator coefficients are updated

G_i→

G_i+Δ

G_i

then the process is repeated with a different value of σ

_trial so that eventually, the “

error”

becomes sufficiently small that solution is obtained.

19. The system as in claim 18 wherein boundary conditions are handled by reflection conditions with tapering and a minimization function is:

${\mathrm{\chi <br><br>}}^2=\underset{i}{\sum <br><br>}<br><br><br><br>[{\left(\frac{C_{\mathrm{ij}}\otimes <br><br>G_j}{{\mathrm{\sigma <br><br>}}_i}-\frac{M_i}{{\mathrm{\sigma <br><br>}}_i}\otimes <br><br>\mathrm{dge}\otimes <br><br>\left(\frac{C_{\mathrm{ik}}\otimes <br><br>G_k}{{\mathrm{\sigma <br><br>}}_i}-\frac{M_i}{{\mathrm{\sigma <br><br>}}_i}\right)\right]}_{x=0,y=0}+\mathrm{\varepsilon <br><br>}<br><br>\left(G_k\otimes <br><br>C_{\mathrm{kj}}\otimes <br><br>G_j\right)<br><br>{❘<br><br>}_{x=0,y=0}$

so the first term equal (# of tube positions)×

(# of pixels) when ε

is set correctly.