U.S. patent application number 12/099416 was filed with the patent office on 2008-09-25 for tomosynthesis imaging system and method.
This patent application is currently assigned to BRANDEIS UNIVERSITY. Invention is credited to Daniel B. Kopans, Richard Moore, Walter Phillips, Martin Stanton, Alex Stewart, Tao Wu.
Application Number | 20080232545 12/099416 |
Document ID | / |
Family ID | 34576441 |
Filed Date | 2008-09-25 |
United States Patent
Application |
20080232545 |
Kind Code |
A1 |
Wu; Tao ; et al. |
September 25, 2008 |
TOMOSYNTHESIS IMAGING SYSTEM AND METHOD
Abstract
A system for three-dimensional tomosynthesis imaging of a target
element is provided having an image acquisition element and a
processor. The image acquisition element obtains a plurality of
images of the target element from a plurality of angles and
includes a radiation source that is positionable at a plurality of
angles with respect to the target element and a radiation detector.
The radiation detector is positioned so as to detect radiation
emitted by the radiation source passing through the target element
and determine a plurality of attenuation values for radiation
passing through the target element to establish a radiation
absorbance projection image of the target element for a particular
radiation source angle. The processor is configured to apply an
iterative reconstruction algorithm to the radiation absorbance
projection images of the target element obtained from a plurality
of radiation source angles to generate a three-dimensional
reconstruction of the target element. The system can gain further
accuracy where the iterative reconstruction algorithm is applied
using cone-beam forward projection and back projection.
Inventors: |
Wu; Tao; (Woburn, MA)
; Stewart; Alex; (Waltham, MA) ; Stanton;
Martin; (Stow, MA) ; Phillips; Walter;
(Arlington, MA) ; Kopans; Daniel B.; (Waban,
MA) ; Moore; Richard; (Concord, MA) |
Correspondence
Address: |
NUTTER MCCLENNEN & FISH LLP
WORLD TRADE CENTER WEST, 155 SEAPORT BOULEVARD
BOSTON
MA
02210-2604
US
|
Assignee: |
BRANDEIS UNIVERSITY
Waltham
MA
|
Family ID: |
34576441 |
Appl. No.: |
12/099416 |
Filed: |
April 8, 2008 |
Related U.S. Patent Documents
|
|
|
|
|
|
Application
Number |
Filing Date |
Patent Number |
|
|
10776690 |
Feb 11, 2004 |
7356113 |
|
|
12099416 |
|
|
|
|
60446784 |
Feb 12, 2003 |
|
|
|
Current U.S.
Class: |
378/22 |
Current CPC
Class: |
G06T 2211/424 20130101;
A61B 6/025 20130101; A61B 6/583 20130101; G06T 11/006 20130101;
G06T 2211/436 20130101; G01N 23/044 20180201; A61B 6/502
20130101 |
Class at
Publication: |
378/22 |
International
Class: |
A61B 6/00 20060101
A61B006/00 |
Claims
1. A tomosynthesis method for creating a three-dimensional
reconstruction of a target element comprising: acquiring radiation
absorbance images of the target element through a limited plurality
of angles; and applying an iterative reconstruction algorithm to
generate the three-dimensional reconstruction of the target
element; wherein the iterative reconstruction algorithm is applied
using cone-beam forward projection and back projection.
2. A method according to claim 1, wherein the radiation absorbance
images are acquired by transmitting x-ray energy from an x-ray
source through the target element to an x-ray detector.
3. A method according to claim 2, wherein the x-ray detector is a
digital x-ray detector having a plurality of detector pixels.
4. A method according to claim 1, wherein radiation absorbance
images are acquired through a number of angles that is less than or
equal to about 100.
5. A method according to claim 1, wherein radiation absorbance
images are acquired through a range of angles that is between about
30 and 120 degrees.
6. A method according to claim 1, wherein the iterative
reconstruction algorithm is a maximum likelihood algorithm.
7. A method according to claim 3, wherein the three-dimensional
reconstruction of the target element is represented as an array of
voxels having a uniform or non-uniform size in
three-dimensions.
8. A method according to claim 7, wherein a forward projection step
is implemented by ray tracing from the x-ray source to a detector
pixel and the forward projection of the target element is obtained
by repeating the ray tracing for each detector pixel to calculate
an attenuation value for each voxel.
9. A method according to claim 8, wherein a back projection step is
implemented by locating detector pixels containing attenuation
information relating to a selected voxel and using those pixels to
update the attenuation value of the selected voxel.
10. A method according to claim 9, wherein the back projection step
includes performing a back projection for at least each voxel
corresponding to a three-dimensional reconstruction of the target
element.
11. A method according to claim 6, wherein the maximum-likelihood
estimation is implemented using an expectation-maximization
algorithm.
12. A method according to claim 1, wherein the target element is at
least a portion of a human patient.
13. A method according to claim 12, wherein the target element is a
breast of a female human patient.
14. A method according to claim 1, wherein the number of iterations
is less than or equal to about 10.
15. A system for three-dimensional tomosynthesis imaging of a
target element comprising: an image acquisition element for
obtaining a plurality of images of the target element from a
plurality of angles having: a radiation source positionable at a
plurality of angles with respect to the target element; and a
radiation detector positioned so as to detect radiation emitted by
the radiation source passing through the target element and
determine a plurality of attenuation value for radiation passing
through the target element to establish a radiation absorbance
projection image of the target element for a particular radiation
source angle; and a processor configured to apply an iterative
reconstruction algorithm to the radiation absorbance projection
images of the target element obtained from a plurality of radiation
source angles to generate a three-dimensional reconstruction of the
target element wherein the iterative reconstruction algorithm is
applied using cone-beam forward projection and back projection.
16. A system according to claim 15, wherein the radiation detector
is a digital x-ray detector having a plurality of detector
pixels.
17. A system according to claim 15, wherein radiation absorbance
projection images are acquired through a number of angles that is
less than or equal to about 100.
18. A system according to claim 15, wherein radiation absorbance
projection images are acquired through a range of angles that is
between about 30 and 120 degrees.
19. A system according to claim 15, wherein the iterative
reconstruction algorithm is a maximum likelihood algorithm.
20. A system according to claim 16, wherein the three-dimensional
reconstruction of the target element is represented as an array of
voxels having a uniform or non-uniform size in
three-dimensions.
21. A system according to claim 20, wherein a forward projection
step is implemented by ray tracing from the radiation source to a
detector pixel and the forward projection of the target element is
obtained by repeating the ray tracing for each detector pixel to
calculate an attenuation value for each voxel.
22. A system according to claim 21, wherein a back projection step
is implemented by locating detector pixels containing attenuation
information relating to a selected voxel and using those pixels to
update the attenuation value of the selected voxel.
23. A system according to claim 22, wherein the back projection
step includes performing a back projection for at least each voxel
corresponding to a three-dimensional reconstruction of the target
element.
24. A system according to claim 19, wherein the maximum-likelihood
estimation is implemented using an expectation-maximization
algorithm.
25. A computer program for three-dimensional tomosynthesis imaging
of a target element from a plurality of radiation absorbance
projection images obtained at a different angles from an image
acquisition element having a radiation source positionable at a
plurality of angles with respect to the target element and a
radiation detector positioned so as to detect radiation emitted by
the radiation source passing through the target element and
determine a plurality of attenuation value for radiation passing
through the target element to establish a radiation absorbance
projection image of the target element for a particular radiation
source angle, the computer program code being embodied in a
computer readable medium and comprising: computer program code for
applying an iterative reconstruction algorithm to the radiation
absorbance projection images of the target element obtained from a
plurality of radiation source angles to generate the
three-dimensional reconstruction of the target element wherein the
iterative reconstruction algorithm is applied using cone-beam
forward projection and back projection.
26. A computer program according to claim 25, wherein the radiation
detector is a digital x-ray detector having a plurality of detector
pixels.
27. A computer program according to claim 25, wherein radiation
absorbance projection images are acquired through a number of
angles that is less than or equal to about 100.
28. A computer program according to claim 25, wherein radiation
absorbance projection images are acquired through a range of angles
that is between about 30 and 120 degrees.
29. A computer program according to claim 25, wherein the iterative
reconstruction algorithm is a maximum likelihood algorithm.
30. A computer program according to claim 26, wherein the
three-dimensional reconstruction of the target element is
represented as an array of voxels having a uniform or non-uniform
size in three-dimensions.
31. A computer program according to claim 30, wherein a forward
projection step is implemented by ray tracing from the radiation
source to a detector pixel and the forward projection of the target
element is obtained by repeating the ray tracing for each detector
pixel to calculate an attenuation value for each voxel.
32. A computer program according to claim 31, wherein a back
projection step is implemented by locating detector pixels
containing attenuation information relating to a selected voxel and
using those pixels to update the attenuation value of the selected
voxel.
33. A computer program according to claim 32, wherein the back
projection step includes performing a back projection for at least
each voxel corresponding to a three-dimensional reconstruction of
the target element.
34. A computer program according to claim 29, wherein the
maximum-likelihood estimation is implemented using an
expectation-maximization algorithm.
Description
RELATED APPLICATIONS
[0001] This application is a continuation of U.S. application Ser.
No. 10/776,690, filed Feb. 11, 2004, issuing as U.S. Pat. No.
7,356,113 on Apr. 8, 2008, which claims priority to U.S.
Provisional Application No. 60/446,784, filed Feb. 12, 2003 and
entitled Tomosynthesis Imaging System and Method, which is
incorporated by reference herein.
BACKGROUND OF THE INVENTION
[0002] 1. Technical Field of the Invention
[0003] The present invention relates to a system and method for
imaging a target element using tomosynthesis. More specifically,
the invention relates to a system, method and computer program
product for creating a three-dimensional image of target elements
from a plurality of radiation absorbance projection images taken
from different angles.
[0004] 2. Background
[0005] Imaging of a patient's tissue has become a common screening
and/or diagnostic tool in modern medicine. One example of such
imaging is mammography, or the imaging of a patient's breast
tissue. Breast cancer remains the most common cancer among women
today, however, at this time there is no certain way to prevent
breast cancer and the best strategy for dealing with breast cancer
is early detection of the cancer so that it may be treated prior to
metastatic spread. Accordingly, it is important for patients to
have access to imaging techniques and systems that will detect very
small cancers as early in their development as possible.
[0006] Conventional mammography involves an x-ray examination of
the breast, typically using a fluorescent panel that converts
transmission x-rays from a breast into visible light photons that
expose a film. While screening using conventional mammography has
been shown to reduce breast cancer deaths by approximately 30 to
50%, this imaging technique lacks the dynamic range that would
allow it to detect small or hidden cancers, and thus permit therapy
that can improve survival rates further. In particular,
conventional mammography techniques suffer from the limitation that
three-dimensional anatomical information is projected onto a two
dimensional image. Because of this, "structure noise" such as
overlapping breast tissues makes it difficult to perceive and
characterize small lesions. This can result in a 10 to 30%
false-negative diagnosis rate, especially where the cancer is
masked by overlying dense fibroglandular tissue.
[0007] A three-dimensional approach to imaging could allow for the
separation of overlying tissue and thus improve correct diagnosis
rates for diseases such as breast cancer; however,
three-dimensional imaging has not yet been applied for this purpose
in the general population. The most widely used three-dimensional
x-ray imaging technique is computed tomography ("CT"). A CT scanner
contains a rotating frame that has one or more x-ray tubes mounted
on one side and one or more detectors on the opposite side. As the
rotating frame spins both the x-ray tube and the detector around
the patient, numerous projections of the x-ray beam attenuated by a
cross section slice of the body are acquired. These projections are
then used to reconstruct cross-sectional images of the body.
Despite the fact that CT has been found useful in detecting lesions
in the breast, it is not suitable as a technique for regular breast
imaging due to the high dose required to take a number of
projections (approximately 100 to 1,000 projections) and the low
spatial resolution (on the order of a millimeter). In addition, the
CT projections mix attenuation effects from other organs of the
body (such as those within the chest cavity) with the attenuation
of the breast, which can distort information about the breast and
causes these interposed organs to be irradiated. Still further, the
cost of CT scanning is too high to permit its use as part of an
annual exam.
[0008] A three-dimensional imaging approach called "tomosynthesis"
has also been developed. Tomosynthesis is a technique that allows
the reconstruction of tomographic planes on the basis of the
information contained in a series of projections acquired from a
series of viewpoints about the target object. They need not be
regularly spaced, numerous, or arranged in any regular geometry.
The tomosynthesis technique is promising in that it may provide
improved spatial differentiation of nearby tissues at very high
resolution comparable to projection 2D imaging, with limited
radiation. The problem of 3D reconstruction from tomosynthesis
projections has been described as intractable by those skilled in
the art.
[0009] In order for a three-dimensional imaging technique to be
successful in medical diagnosis and other applications, it should
offer: [0010] Sufficient spatial resolution and contrast resolution
to detect and characterize, for example, breast cancers; [0011]
Minimum radiation dose to a patient; [0012] Fast image acquisition;
[0013] Cost effectiveness; and [0014] 3D reconstruction that can be
performed effectively.
SUMMARY OF THE INVENTION
[0015] In one aspect, the invention provides a method that enables
to the use of tomosynthesis to efficiently provide accurate
three-dimensional imaging of a target element. This method involves
acquiring radiation absorbance images of the target element through
a limited plurality of angles and applying an iterative
reconstruction algorithm to generate the three-dimensional
reconstruction of the target element. This method can gain further
accuracy where the iterative reconstruction algorithm is applied
using cone beam forward projection and back projection.
[0016] In a further aspect of the invention, a system for
three-dimensional tomosynthesis imaging of a target element is
provided having an image acquisition element and a processor. The
image acquisition element obtains a plurality of images of the
target element from a plurality of angles and includes a radiation
source that is positionable at a plurality of angles with respect
to the target element and a radiation detector. The radiation
detector is positioned so as to detect radiation emitted by the
radiation source passing through the target element and determine a
plurality of attenuation values for radiation passing through the
target element to establish a radiation absorbance projection image
of the target element for a particular radiation source angle. The
processor is configured to apply an iterative reconstruction
algorithm to the radiation absorbance projection images of the
target element obtained from a plurality of radiation source angles
to generate a three-dimensional reconstruction of the target
element. Again, the system can gain further accuracy where the
iterative reconstruction algorithm is applied using cone-beam
forward projection and back projection.
[0017] In a still further aspect of the invention, a computer
program for three-dimensional tomosynthesis imaging of a target
element is provided. The three-dimensional images are created from
a plurality of radiation absorbance projection images obtained at
different angles from an image acquisition element having a
radiation source positionable at a plurality of angles with respect
to the target element and a radiation detector. The radiation
detector is positioned so as to detect radiation emitted by the
radiation source passing through the target element and determine a
plurality of attenuation values for radiation passing through the
target element to establish a radiation absorbance projection image
of the target element for a particular radiation source angle. The
computer program code is embodied in a computer readable medium and
includes computer program code for applying an iterative
reconstruction algorithm to the radiation absorbance projection
images of the target element obtained from a plurality of radiation
source angles to generate the three-dimensional reconstruction of
the target element wherein the iterative reconstruction algorithm
is applied using cone-beam forward projection and back
projection.
[0018] In specific embodiments of any of these aspects of the
invention, the radiation absorbance images can be acquired by
transmitting x-ray energy from an x-ray source through the target
element to an x-ray detector and the x-ray detector may have a
plurality of detector pixels. The three-dimensional reconstruction
of the target element may be represented as an array of voxels
having a uniform or non-uniform size in three-dimensions. The
forward projection step may then be implemented by ray tracing from
the x-ray source to a detector pixel and the forward projection of
the target element is obtained by repeating the ray tracing for
each detector pixel to calculate an attenuation value for each
voxel. The back projection step can be implemented by locating
detector pixels containing attenuation information relating to a
selected voxel and using those pixels to update the attenuation
value of the selected voxel. The back projection step can further
include performing a back projection for at least each voxel
corresponding to a three-dimensional reconstruction of the target
element. In the enumerated aspects of the invention or in any of
their embodiments, the iterative reconstruction algorithm may be a
maximum likelihood algorithm and the maximum likelihood estimation
can be implemented using an expectation-maximization algorithm.
[0019] The invention is particularly useful for creating
three-dimensional reconstructions of animal and more particularly
human tissue. In one preferred embodiment, the invention is
employed in mammography to create a three-dimensional
reconstruction of the breast tissue of a human female patient.
BRIEF DESCRIPTION OF THE DRAWINGS
[0020] The invention will be more fully understood from the
following detailed description taken in conjunction with the
accompanying drawings:
[0021] FIG. 1 is a diagram illustrating the geometry of a
tomosynthesis system of the invention;
[0022] FIG. 2 is a top view of the coordinate system of a
tomosynthesis system of the invention;
[0023] FIG. 3 illustrates a forward projection of the tomosynthesis
system of the invention;
[0024] FIG. 4 illustrates the path length of an x-ray beam in a
voxel in the tomosynthesis system of the invention;
[0025] FIG. 5 illustrates a projection of the path length of FIG.
4;
[0026] FIG. 6 illustrates an exception to the projection of FIG.
5;
[0027] FIG. 7 illustrates a back-projection step of the
invention;
[0028] FIG. 8A illustrates a phantom used to test a system of the
invention;
[0029] FIG. 8B illustrates a feature plate that makes up a portion
of the phantom of FIG. 8A;
[0030] FIGS. 9A, 9B, and 9C illustrate structural noise reduction
in (A) a projection of an ACR phantom; (B) a projection of a
mastectomy specimen/ACR phantom; and (C) a reconstructed ACR
phantom feature layer;
[0031] FIG. 10 is a film-screen mammogram of a patient's tissue;
and
[0032] FIGS. 11A, 11B, and 11C illustrate slices of a reconstructed
volume of the same tissue at three different depths.
DETAILED DESCRIPTION OF THE INVENTION
[0033] The systems and methods of the present invention address the
needs of the art by providing tomosynthesis apparatus and
techniques for imaging target elements that overcome the problems
of conventional three-dimensional imaging systems. The present
invention enables the use of tomosynthesis to efficiently provide
accurate three-dimensional imaging of a target element in which
overlapping sub-elements having differing attenuation
characteristics by applying a 3D reconstruction algorithm having a
novel combination of features. The algorithm can employ a cone-beam
geometry lacking in geometric simplification such as parallel-beam
based approximation methods. The algorithm can further apply the
cone-beam geometry in an iterative forward-projection and
back-projection method based on maximum-likelihood image estimation
using an estimation-maximization algorithm. The invention is
applied below to one preferred embodiment in which the system is
used for tomosynthesis mammography; however, the invention will be
useful in a variety of three-dimensional imaging situations. For
example, the invention can be applied to a variety of patient
imaging problems such as heart imaging, or imaging of the soft
tissues or bones of the hand. The imaging system of the invention
can be used for diagnoses (as is described below for tomosynthesis
mammography) or it may be used for other applications such as
three-dimensional modeling for the purpose of fitting an implant
(whether orthopedic, such as a hip or knee implant, an artificial
heart, or other type of implant) or for use in surgical navigation
systems. What follows is a description of one preferred embodiment
of the invention.
1. Tomosynthesis Mammography System
[0034] Tomosynthesis mammography is a three-dimensional breast
imaging technique. It involves acquiring projection images of a
breast at a plurality of viewpoints, typically over an arc or
linear path. Three-dimensional distribution of x-ray attenuation
coefficient of the breast volume is reconstructed from these
projections. A prototype tomosynthesis system 10 for breast imaging
is illustrated in FIG. 1. In this exemplary system, eleven
projections are acquired by moving the x-ray tube 12 over a
50.degree. arc (-25.degree. to +25.degree.) above the target
element, in this case breast tissue 18 which may be compressed by
compression paddle 16, in 5.degree. angular steps about axis of
rotation 14. Breast tissue 18 and digital detector 20 are
stationary during the image acquisition. Certain characteristics of
this exemplary embodiment of a tomosynthesis system of the
invention are described below:
Spatial resolution and contrast resolution: The tomosynthesis
system uses an amorphous-Silicon-based flat panel digital detector
20 on which a CsI crystal phosphor is grown epitaxially. It reads
out 2304.times.1800 pixels (100 .mu.m pixel pitch) via a TFT array.
The detector has a linear response over exposure levels up to 4000
mR and 12 bits of working dynamic range. Each plane of the 3D
reconstruction has about the same resolution as the detector (100
um) but the depth resolution is on the order of a millimeter. Dose:
The target/filter combination is Rh/Rh and the accelerating
potential is 25.about.33 kVp to image breasts with 3.about.8 cm
range of thickness. The total x-ray dose for acquiring 11
projections is approximately 1.5 times of that used for one
film-screen mammogram. Each projection is a low dose breast image
(approximately 1/11 of the does per projection). Patient motion:
Patient motion is reduced by fast image acquisition. Using
cone-beam x-ray geometry and area detector, a projection of the
whole breast can be recorded with one x-ray exposure at each angle.
For each projection, the exposure time is 0.1.about.0.2 s and
detector readout time is about 0.3 s. Rotation to the next angle is
performed during the detector readout. The total image acquisition
time for 11 projections is about 7 sec. Breast compression also
helps to reduce patient motion. Image acquisition geometry: The
design of the tomosynthesis system can be based on the conventional
mammography system. The MLO views have been used in most cases
since it provides the most complete coverage of the whole
breast.
2. 3D Reconstruction Algorithm
[0035] Tomosynthesis can take advantage of the high efficiency of a
digital detector in acquiring low dose breast images. Prior to the
present invention, appropriate reconstruction methods that make
good use of the low dose projections and the acquisition geometry
of the tomosynthesis system 10 have not been deployed. For an
initial evaluation, Niklason implemented a "shift-and-add" method
that is similar to backprojection [Niklason et al, 1997]. Methods
used by others [Chakraborty et al, 1984; Haaker et al, 1985;
Suryanarayanan et al, 2000] essentially did not handle the limited
statistics in low dose projection images. In theory, they were not
suitable in the case of limited number of projections and limited
angular range. Therefore, the three-dimensional information
extracted by these methods was limited, which resulted in poor
quality reconstructions.
The Maximum Likelihood (ML) algorithm is an iterative
reconstruction method [Rockmore, 1977; Shepp et al, 1982; Levitan
et al, 1987; Herbert et al, 1989; Browne et al, 1992; Manglos et
al, 1995; Pan et al 1997; Zhou et al, 1997]. It is well suited for
tomosynthesis reconstruction, which is an ill-conditioned problem
(only 11 low dose projections are available). The ML algorithm
incorporates the stochastic nature of the x-ray transmission
process so that the statistical noise in projection images is taken
into consideration in the case of low x-ray flux. It also
incorporates the information of the object into the reconstruction
in the form of constraints.
[0036] In ML reconstruction, the Likelihood function, which is the
probability of obtaining the projections Y obtained in a
measurement, given a certain model for the three-dimensional map of
attenuation coefficients u is:
L=P(Y|u) (1)
The ML solution is the 3D reconstruction that maximizes the
probability of the measured projections. Because an analytical
solution is usually intractable, an iterative algorithm is a better
choice. The incident and transmitted x-rays follow Poisson
statistics and the log-likelihood is described by:
LnL = i ( - N i - < l , u > i - Y i < l , u > i + Y i
ln N i - ln Y i ) ( 2 ) ##EQU00001##
where u is the linear attenuation coefficient; N.sub.i is the
number of incident x-ray photons to projection pixel i, before
attenuation; Y.sub.i is the number of transmitted x-ray photons to
projection pixel i, after attenuation; l.sub.ij is the path length
of beam ray i in the object (reconstruction voxel j; and
< l , u > i = i l ij u j ##EQU00002##
is the total attenuation along beam ray to pixel i.
[0037] The algorithm by Lange and Fessler [Lange and Fessler, 1995]
can be selected to solve the ML problem. At the n-th iteration, the
value of an object voxel .mu. is updated by:
u j ( n + 1 ) = u j ( n ) + u j ( n ) i l ij ( N i - < l , u ( n
) > i - Y i ) i ( l ij < l , u ( n ) > i N i - < l , u
( n ) > i ) ( 3 ) ##EQU00003##
where the notations are the same as above.
3. Implementation of the Cone-Beam Reconstruction
[0038] Cone-beam forward projection and back projection can form
the basis for iterative reconstruction according to the invention.
At the forward projection step, the projection images at 11 angles
are calculated based on the current 3D reconstruction model. At the
backprojection step, the calculated projections and the measured
projections are compared and the 3D reconstruction model is updated
according to their difference.
[0039] The forward projection to a detector pixel i at a projection
angle can be used to illustrate the whole forward projection
problem. An x-ray beam containing N.sub.i photons is incident from
the source to the center of the selected detector pixel. This beam
penetrates a series of object voxels and is sequentially attenuated
by them. The total aggregate attenuation is
<l,u.sup.(n)>.sub.i and the number of transmitted photons is
N.sub.ie.sup.-<l,u.sup.(n).sup.>.sup.i, which is the forward
projection to the pixel. This operation is repeated for all
detector pixels that form the forward projection at this angle. The
forward projections at all angles can be done in the same way
except that the "pseudo-beam" is rotated.
[0040] The 3D reconstruction model is updated at the backprojection
step. Equation 3 describes the update of a voxel j at the n-th
iteration of reconstruction. The whole image is updated by doing
the same operation on every voxel in it. At a projection angle, the
center of the voxel is projected from the source to a detector
pixel containing the attenuation information of this voxel. This
operation is repeated at other angles and totally 11 detector
pixels are found. In equation 3, the values of these 11 pixels,
both in forward projection and in measured projection are used to
update the object (reconstruction) voxel (the summation is on the
set of these 11 pixels).
[0041] 3.1 Positions of X-Ray Source, Object Voxel and Detector
Pixel
[0042] The origin of the coordinate system is at the axis of
rotation 14 as illustrated in FIG. 2. The rotation plane of the
x-ray 12 source is the YZ-plane (x=0). The detector 20 is parallel
to the XY-plane at z=21.7 cm. The distance between the source 12
and the axis of rotation 14 is D.sub.sa and the distance between
the detector 20 and the axis of rotation 14 is D.sub.da. At
projection angle .theta., the position of the x-ray source 12
is:
x.sub.s(.theta.)=0
y.sub.s(.theta.)=D.sub.sasin(.theta.) (4)
z.sub.s(.theta.)=-D.sub.sacos(.theta.)
[0043] The reconstructed object 24 is a rectangular volume,
represented by a three-dimensional array of voxels 26. The breast
volume 18 is contained in this rectangular volume 24. In a
reconstructed image, the value of a voxel is positive if it
represents breast tissue; zero if it represents the empty space out
of the breast. In the coordinate system, the position of a voxel 26
indexed by (m.sub.x, m.sub.y, m.sub.z) is:
x.sub.obj=X.sub.obj+m.sub.xd.sub.x
y.sub.obj=Y.sub.obj+m.sub.yd.sub.y (5)
z.sub.obj=Z.sub.obj+m.sub.zd.sub.z
where (X.sub.obj, Y.sub.obj, Z.sub.obj) is the position of the
center of the rectangular volume 24; d.sub.x, d.sub.y and d.sub.z
are the size of the voxel 26 in three dimensions.
[0044] The position of a detector pixel 28 indexed by (n.sub.x,
n.sub.y) is:
x.sub.p=X.sub.p+n.sub.xd'.sub.x
y.sub.p=Y.sub.p+n.sub.yd'.sub.y (6)
z.sub.p=D.sub.da
where (X.sub.p, Y.sub.p) is the position of the center of the
detector 20; d'.sub.x and d'.sub.y are the size of the pixel 28 in
X and Y dimensions.
[0045] 3.2 Forward Projection
[0046] The forward projection is implemented by ray tracing from
the x-ray source 12 to detector pixel 28. At a projection angle,
the x-ray beam to a detector pixel 12 is attenuated from the point
where the beam enters the volume 24 to the point where it goes out.
The total attenuation along the beam <l,u.sup.(n)>.sub.i is
calculated by accumulating the attenuation lu.sup.(n) by each voxel
26 on the beam line. The number of transmitted x-rays to the pixel
28 is N.sub.ie.sup.-<l,u.sup.(n).sup.>.sup.i. The forward
projection of the object 18 at this angle is obtained by repeating
this operation for all detector pixels 28. The forward projections
at other angles are calculated in the same way except the x-ray
source 12 is at a different location.
[0047] The first step of forward projection is to determine the
orientation of the x-ray beam 30 as illustrated in FIG. 3. At an
angle, the position of the x-ray source (x.sub.s, y.sub.s, z.sub.s)
12 and detector pixel (x.sub.p, y.sub.p, z.sub.p) 28 are determined
by equation 4 and 6. The orientation of the beam P.sub.(x,y,z) 30
from source 12 to the detector pixel 28 can be described by two
parameters: (1) .beta., the angle made by the beam and the
YZ-plane; (2) .alpha., the angle made by the projection of the beam
in YZ-plane and the Z-axis. These two parameters are determined
by:
.alpha.=tan.sup.-1((y.sub.p-y.sub.s)/(z.sub.p-z.sub.s))
.beta.=tan.sup.-1(x.sub.p/ {square root over
((y.sub.p-y.sub.s).sup.2+(z.sub.p-z.sub.s).sup.2)}{square root over
((y.sub.p-y.sub.s).sup.2+(z.sub.p-z.sub.s).sup.2)}) (7)
[0048] The path length |{right arrow over (P.sub.1P.sub.4)}| of the
x-ray beam 30 through a voxel 26, as illustrated in FIG. 4, is also
the distance between the centers of two successive voxels along the
beam. The position of the next voxel along the beam can be located
by shifting .DELTA.x, .DELTA.y and .DELTA.z ({right arrow over
(P.sub.1P.sub.2)}, {right arrow over (P.sub.2P.sub.3)} and {right
arrow over (P.sub.3P.sub.4)} in FIG. 4) along three dimensions from
the current voxel 26.
.DELTA.x={right arrow over (P.sub.1P.sub.2)}={right arrow over
(P.sub.1P.sub.4)}cos .beta.cos .alpha.
.DELTA.y={right arrow over (P.sub.2P.sub.3)}={right arrow over
(P.sub.1P.sub.4)}cos .beta.sin .alpha. (8)
.DELTA.z={right arrow over (P.sub.3P.sub.4)}={right arrow over
(P.sub.1P.sub.4)}sin .beta.
[0049] To calculate {right arrow over (P.sub.1P.sub.4)}, its
projection in the YZ-plane {right arrow over (P.sub.1P.sub.3)},
illustrated in FIG. 5, is calculated first:
{right arrow over (P.sub.1P.sub.3)}=d.sub.y/sin .alpha. if
.alpha.>tan.sup.-1(d.sub.y/d.sub.z); (9)
{right arrow over (P.sub.1P.sub.3)}=d.sub.x/cos .alpha. if
.alpha..ltoreq.tan.sup.-1(d.sub.y/d.sub.z)
[0050] In a similar way, the path length {right arrow over
(P.sub.1P.sub.4)} can be calculated by:
{right arrow over (P.sub.1P.sub.4)}=d.sub.x/sin .beta. if
.beta.>tan.sup.-1(d.sub.x/{right arrow over (P.sub.1P.sub.3)});
(10)
{right arrow over (P.sub.1P.sub.4)}=d.sub.x/cos .beta. if
.beta..ltoreq.tan.sup.-1(d.sub.x/{right arrow over
(P.sub.1P.sub.3)})
[0051] There are exceptions to the two cases illustrated in FIG. 5.
In a case shown in FIG. 6, the path lengths through voxel 3 and 4
cannot be described by equation 10. But the total path length of
them is equal to the path length in voxel 2. The total attenuation
by voxel 3 and 4 is equivalent to the attenuation by the shaded
area in FIG. 6, which has the same path length as voxel 2. The
equivalent attenuation is estimated by a linear interpolation of
attenuations by voxel 3 and 4. The weighting for the interpolation
is proportional to the inverse of the distance from the voxel
center to the beam line. The ratio of the weighting for voxel 3 to
that for voxel 4 is d.sub.4/d.sub.3, equivalent to r.sub.4/r.sub.3,
where d.sub.3 and d.sub.4 are the distances from the voxel center
to the beam; r.sub.3 and r.sub.4 are the distances from the voxel
center to the projection of the beam along the Y-axis.
[0052] The total attenuation along a beam to a detector pixel i is
the summation from the first voxel at the point where the beam
enters the volume to the voxel at the point where the beam goes out
of the volume. For a beam with orientation (.alpha., .beta.), the
position of the voxel at entering point is:
x.sub.0=x.sub.s+ {square root over
((y.sub.0-y.sub.s).sup.2+(z.sub.0-z.sub.s).sup.2)}{square root over
((y.sub.0-y.sub.s).sup.2+(z.sub.0-z.sub.s).sup.2)}tan .beta.
y.sub.0=y.sub.s+(z.sub.0-z.sub.s)tan .alpha.; (11)
z.sub.0=21.7-D;
where D is the thickness of the reconstruction volume. The
attenuation lu.sub.0 by the first voxel at (x.sub.0, y.sub.0,
z.sub.0) is calculated and then the tracing point is shifted
forward by (.DELTA.x, .DELTA.y, .DELTA.z) to the next voxel along
the beam, where the attenuation lu.sub.1 is calculated and added to
lu.sub.0. At the n-th step, the position being search is:
x.sub.n=x.sub.0+n.DELTA.x
y.sub.n=y.sub.0+n.DELTA.y (12)
z.sub.n=y.sub.0+n.DELTA.z
The number of steps of forward projection is V=int(D/.DELTA.z)+1.
After V steps, the total attenuation along the beam to detector
pixel i is
n = 0 M u n l n ##EQU00004##
(represented by <l,u.sup.(n)>.sub.i). The number of
transmitted x-ray photons is
N i - n = 0 N u n l n . ##EQU00005##
[0053] 3.3 Backprojection
[0054] The value of the object voxel is updated at the
backprojection step as illustrated in FIG. 7. At this step,
projection pixels containing the attenuation information of the
selected object voxel are found and used to update the value of
this voxel. At a projection angle, the position of the detector
pixel (x.sub.p, y.sub.p, z.sub.p) which contains the information of
a selected voxel is:
x.sub.p=x.sub.s+(x.sub.obj-x.sub.s)(z.sub.p-z.sub.s)/(z.sub.obj-z.sub.s)
y.sub.p=y.sub.s+(y.sub.obj-y.sub.s)(x.sub.p-x.sub.s)/(x.sub.obj-x.sub.s)
(13)
z.sub.p=21.7
where (x.sub.s, y.sub.s, z.sub.s) is the position of the x-ray
source at this angle. This operation is repeated to find detector
pixels related to this voxel at other angles. The value of this
voxel is updated by equation 3, using these detector pixels.
4. Image Reconstruction Results
[0055] 4.1 Study on an ACR Phantom/Mastectomy Specimen
[0056] A phantom 38 is composed of a piece of mastectomy specimen
40 and a feature plate 42 from an American College of Radiology
(ACR) accredited mammography phantom and placed on detector 20 as
illustrated in FIG. 8A. The feature plate 42, further illustrated
in FIG. 8B, contained nylon fibers (labeled 1 to 6 on the plate),
simulated micro-calcifications (labeled 7 to 11 on the plate) and
tumor-like masses (labeled 12 to 16 on the plate). The mastectomy
specimen 40 is a surgically removed breast tissue containing
lesions. The combination of the feature plate 42 with the
mastectomy specimen 40 makes it very hard to find features of the
ACR phantom 42. The reconstructed feature plate demonstrates how
the three-dimensional reconstruction works to improve the
visibility of features.
[0057] Ten features (fiber 1, 2, 3, 4; micro-calcification cluster
7, 8, 9 and mass 12, 13, 14) can be seen very well in a projection
of the 4 cm thick ACR phantom 42 itself (Rh/Rh, 28 kVp and 160 mAs)
as shown in FIG. 9A. With the superimposed mastectomy specimen 40,
only one feature (micro-calcification cluster 7) is visible in a
projection (Rh/Rh, 30 kVp and 140 mAs) as can be seen in FIG.
9B.
[0058] The reconstruction of the feature layer after 10 iterations
is shown in FIG. 9C. The x-ray energy and exposure are the same as
that used to create the image of FIG. 9B. More features
(micro-calcification cluster 7, 8, 9 and mass 12) can be seen in
the reconstruction. Even some low contrast features (fiber 1, 2, 3,
4) are recognizable. The number "503 059" on the label is clearer.
It is clear that the visibility of features are significantly
improved.
[0059] 4.2 3D Reconstruction of a Patient Tissue
[0060] Clinical imaging of volunteers conducted at Massachusetts
General Hospital under IRB approved protocols have been
reconstructed for comparison of conventional film-screen
mammography and to tomosynthesis mammography. As an example, a
mediolateral oblique (MLO) mammogram from a volunteer was obtained
using film-screen system (Mo/Mo, 25 kV and 330 mrad average
glandular dose). The x-ray film image is shown in FIG. 10. The
patient was found to have a non-palpable 10 mm invasive ductal
cancer with associated in situ tumor and this was proved by biopsy.
The cancer was difficult to see in the conventional screening
mammogram and was found primarily because the calcifications
associated with it drew the attention of the radiologist.
[0061] A tomosynthesis image dataset was taken with Rh/Rh
target/filter at 28 kVp and a total dose of 307 mrad. Three
reconstructed slices from the 3D reconstruction are shown in FIG.
11. Blood vessels are seen near the breast skin in FIG. 11A. A
tumor that has intraductal as well as invasive ductal cancer
elements is just out of the plane of section in FIG. 11B. The
invasive tumor mass, marked by an arrow, with associated
calcifications in the in situ portion is clearly seen in FIG. 11C,
as is a benign intramammary lymph node in the upper portion of the
image.
[0062] It is apparent from this volunteer's dataset that
overlapping structures in the conventional two-dimensional
projection images (FIG. 10) were spacially separated. A
reconstructed image provided at three different depths (FIG. 11A
illustrating a depth of Z=2 mm, FIG. 11B illustrating a depth of
Z=22 mm, and FIG. 11C illustrating a depth of Z=32 mm) makes it
easier to see the tumor and calcifications and their relative
geometry.
[0063] A person of ordinary skill in the art will appreciate
further features and advantages of the invention based on the
above-described embodiments. Accordingly, the invention is not to
be limited by what has been particularly shown and described,
except as indicated by the appended claims or those ultimately
provided in a utility application claiming priority to this
provisional application. A number of references have been referred
to in the specification by last name of the first listed author and
year of publication; those references are listed by full citation
in the Bibliography below. All publications and references cited
herein are expressly incorporated herein by reference in their
entirety, in particular, each of the references listed in the
Bibliography below is expressly incorporated for the teachings
referred to in the sections of the application above for which they
are cited.
BIBLIOGRAPHY
[0064] U.S. Pat. No. 5,872,828 to Niklason et al., entitled
"Tomosynthesis System for Breast Imaging." [0065] J. A. Browne, and
T. J. Holmes, "Developments with Maximum Likelihood X-ray Computed
Tomography," IEEE Transactions on Medical Imaging, 11(1): 40-52
(1992). [0066] D. P. Chakraborty, M. V. Yester, G. T. Barnes and A.
V. Lakshminarayanan, "Self-masking subtraction tomosynthesis,"
Radiology, 150:225-229 (1984). [0067] P. Haaker, E. Klotz, R.
Koppe, R. Linde and H. Moller, "A New Digital Tomosynthesis Method
With Less Artifacts for Angiography," Medical Physics, 12(4):
431-436 (1985). [0068] T. J. Herbert and R. M. Leahy, "A
Generalized EM Algorithm for 3-D Bayesian Reconstruction from
Poisson Data Using Gibbs Priors," IEEE Transactions on Medical
Imaging, 8(2): 194-202 (1989). [0069] K. Lange and J. A. Fessler,
"Globally Convergent Algorithm for Maximum a Posteriori
Transmission Tomography," IEEE Transactions on Image Processing, 4:
1430-1438 (1995). [0070] E. Levitan and G. T. Herman, "A Maximum A
Posteriori Probability expectation Maximization Algorithm or Image
Reconstruction in Emission Tomography," IEEE Transactions on
Medical Imaging, MI-6(3): 185-192 (1987). [0071] S. H. Manglos, G.
M. Gagne, F. D. Thomas and R. Narayanaswamy, "Transmission
Maximum-Likelihood Reconstruction with Ordered Subsets for Cone
Beam CT," Physics in Medicine and Biology, 40: 1225-1241 (1995).
[0072] L. T. Niklason, B. T. Christian, L. E. Niklason, D. B.
Kopans, D. E. Castleberry, B. H. Opsahl-Ong, C. E. Landberg, P. J.
Slanetz, A. A. Giardino, R. M. Moore, D. Albagi, M. C. Dejule, P.
A. Fitzgerald, D. F. Fobare, B. W. Giambattista, R. F. Kwasnick, J.
Liu, S. J. Lubowski, G. E. Possin, J. F. Richotte, C-Y Weinad R. F.
Wirth, "Digital Tomosynthesis in Breast Imaging," Radiology, 205:
399-406 (1997). [0073] L. T. Niklason, B. T. Christian, L. E.
Niklason, D. B. Kopans, P. J. Slanetz, D. E. Castleberry, B. H.
Opsahl-Ong, C. E. Landberg, B. W. Giambattista, "Digital Breast
Tomosynthesis: Potentially a New Method for Breast Cancer
Screening," Digital Mammography, edited by N. Karssemeijer, M.
Thijssen, J. Hendriks and L. van Erning, 13: 51-56 (Kluwer Academic
Publishers, 1998). [0074] T. Pan, B. M. W. Tsui and C. L. Byrne,
"Choice of Initial Conditions in the ML Reconstruction of Fan-Beam
Transmission with Truncated Projection Data," IEEE Transactions on
Medical Imaging, 16(4): 426-438 (1997). [0075] A. J. Rockmore and
A. Macovski, "A Maximum Likelihood Approach to Transmission Image
Reconstruction From Projections," IEEE Transactions on Nuclear
Science, 24: 1929-1935 (1977). [0076] L. A. Shepp and Y. Vardi,
"Maximum Likelihood Reconstruction for Emission Tomography," IEEE
Transactions on Medical Imaging, MI-1: 113-122 (1982). [0077] S.
Suryanarayanan, A. Karellas, S. Vedantham, S. J. Glick, C. J.
D'Orsi, S. P. Baker and R. L. Webber, "Comparison of Tomosynthesis
Methods Used with Digital Mammography," Academic Radiology,
7:1085-1097, (2000). [0078] Z. Zhou, R. M. Leahy and J. Qi,
"Approximate Maximum Likelihood Hyperparameter Estimation for Gibbs
Priors," IEEE Transactions on Image Processing, 6(6): 844-861
(1997).
* * * * *