Method for improved breast x-ray imaging

[0002] 1. Field of the Invention

[0003] The present invention relates generally to digital xray imaging and, more particularly, to a method of digital mammography that uses dual-energy apparatus and methods for separating a single human breast x-ray image into component images, each representing a single physical substance.

[0004] 2. The Prior Art

[0005] Typical clinical mammography is performed by using x-ray films. Recently, large format two-dimensional semiconductor digital x-ray detector arrays have been available. Currently, regardless of the technology, all of the image information is contained in a single image, acquired through the use of x-rays with a single energy-spectrum. The human breast is comprised of three major substances, the lean tissue, the fat tissue, and microcalcification deposits, and each pixel of such a single image contains a mixture of all three, plus a random scatter component. The size of the contribution from each component is not known in current mammography.

[0006] It is well established that the role of the random scatter signal in an x-ray imaging is interference and distortion. The scatter blurs the image, reduces the image contrast, and degrades the image quality. The contribution of scatter in mammography is generally as high as 30% to 60% of the total image signal.

[0007] It has also been established by breast cancer research, that the cancer cells have an x-ray attenuation coefficient very close to that of glandular duct which has an average composition of typical lean tissue. Thus, only the lean tissue component provides useful information for the diagnosis of cancer. The direct value of the fat tissue information for the diagnosis is very small. At the same time, the fat tissue or adipose tissue occupies as much as about 50% of the total breast volume and is closely mingled with the lean tissue. Thus, in mammography, the fat tissue signal also acts more as an adverse factor, imposing a strong, irregular pattern to mask the useful information.

[0008] The image signal of lean tissue is used for direct determination of breast cancer or other pathological changes. The microcalcification component is an important reference because it has been found that most breast cancer patients have a large amount of microcalcification in breast tissue.

[0009] Because of the inability to separate the four basic signal components, the capability of current x-ray imaging for cancer diagnosis has been essentially limited.

[0010] The present invention is based on the dual-energy x-ray imaging method disclosed in U.S. Pat. Nos. 5,648,997 and 5,771,269 and U.S. patent application Ser. No. 09/025,926 (the Chao disclosures). In terms of the prior art, the applicability of dual-energy x-ray imaging to mammography is discussed in several journal articles. The method and apparatuses used in these articles are essentially different from the present invention. For example, Chakraborty et al. (An Energy Sensitive Cassette for Dual-Energy Mammography, 16(1) Medical Physics 7 (January/February 1989)) designed a double film cassette to use a “dual-energy subtraction method” for enhancing image contrast. Johns and Yaffe (Theoretical Optimization of Dual-Energy x-Ray Imaging with Application to Mammography, 12(3) Medical Physics 289 (May, June 1985)) conducted experiments using linear detector arrays and linearization approximations for dual energy decomposition. First, the hardware between these references and the Chao disclosures is different; the detectors of the prior art are either film cassettes, linear detector arrays, or stimulable phosphor plates. Using an intermediate digitizing method to convert the analog image data either from x-ray films or from stimulable phosphor plates into a digital format is essentially a semi-quantitative imaging system. The present invention uses two-dimensional large format integrated semiconductor detectors for direct high-accuracy quantitative imaging. Secondly, the dual-energy decomposition in the prior art is invariably based on linearization methods, while the present invention is based on directly solving a dual-energy fundamental equation system in its original form without relying on linearization. Thirdly, the present invention uses imaging hardware which combines scatter removal capability and dual-energy data acquisition in one system.

[0011] An object of the present invention is to provide a method for enhancing the diagnosis value of x-ray images for mammography.

[0012] The present invention provides a method for quantitatively separating a breast x-ray image into a number of component images: the scatter image, the fat tissue image, the lean tissue image, and the microcalcification image. The present invention uses dual-energy x-ray imaging system hardware configurations as described in U.S. Pat. Nos. 5,648,997 and 5,771,269 and U.S. patent application Ser. No. 09/025,926 (the Chao disclosures), incorporated herein by reference. The present invention does not provide any new hardware.

[0013] As to the method of the present invention, the basic methods for removing scatter and for dual-energy decomposition described in the Chao disclosures are used, and are not an aspect of the present invention. However, the method of the present invention includes significant improvements and is not a simple implementation of these patented methods. To separate the three material composition by using the basic dual-energy method, new procedures are required. The present invention provides a systematic method to achieve decomposition of a breast image into multiple component images according to physical substance.

[0014] The method of the present invention is to use the basic method for removal of scatter and the basic method for dual-energy x-ray imaging to first separate a mixed breast image into four basic image components: a scatter image, a lean tissue image, a fat tissue image, and a microcalcification image. “Microcalcification” is used interchangeably with “calcification”. These images are a first order approximation. Then the three material compositions of the human breast are taken into account. In the second order approximation, the microcalcification image, lean tissue image, and fat tissue image are separated so that each contains only a single breast component.

[0015] The method of obtaining the first order approximations includes the steps of (a) performing a calibration, as described below, to obtain a pair of numerical relationships for the front detector at the high and low energies for the microcalcification c and the soft tissue s to obtain the functions c=c(D_H,D_L) and s=s(D_H,D_L), (b) perform a calibration to obtain a pair of numerical relationships for the front detector at the high and low energies for the fat tissue f and the lean tissue g to obtain the functions f=f(D_H,D_L) and g=g(D_H,D_L), (c) illuminating the subject with x-rays of said average energy level H and of average energy L, (d) acquiring high-resolution images D_fHh(x,y) and D_fLh(x,y) from the front detection locations (x,y), where the images are composed of both primary and scatter x-rays, (e) calculating a pair of high-resolution scatter x-ray images D_fSHh(x,y) and D_fSLh(x,y), (f) calculating a pair of high-resolution primary x-ray images D_fPHh(x,y) =D_fHh(x,y)−D_fSHh(x,y) and D_fPLh(x,y)=D_fLh(x,y)−D_fSLh(x,y) (g) performing a dual-energy decomposition for the image pair D_fPHh(x,y) and D_fPLh(x,y) using the functions c=c(D_H,D_L) and s=s(D_H,D_L) to obtain two first order approximation material composition images c₁(x,y) and s₁(x,y), and (h) performing a dual-energy decomposition for the image pair D_fPHh(x,y) and D_fPLh(x,y) using the functions f=f(D_H,D_L) and g=g(D_H,D_L) to obtain two first order approximation material composition images f₁(x,y) and g₁(x,y).

[0016] The method of obtaining the second order approximations corrects for microcalcification effects and includes the steps of (a) identifying all microcalcification points c₁(x_k,y_k) and all non-microcalcification points c₁(x_i,y_i) in the image c₁(x,y), (b) constructing a background image B(x,y) where points B(x_i,y_i)=c₁(x_i,y_i) and where points B(x_k,y_k) are interpolated from points c₁(x_j,y_j) around points C₁(x_k,y_k), (c) subtracting the background image B(x,y) from the image c₁(x,y) to obtain a second order approximation microcalcification image c₂(x,y), (d) identifying all zero points c₂(x_o,y_o) and non-zero points c₂(x_m,y_m) in the image c₂(x,y); (e) constructing a second order approximation fat tissue image f₂(x,y) where points f₂(x_o,y_o)=f₁(x_o,y_o) and where points f₂(x_m,y_m) are interpolated from points f₁(x_n,y_n) around points f₁(x_m,y_m); and (f) constructing a second order approximation lean tissue image g₂(x,y) where points g₂(x_o,y_o)=g₁(x_o,y_o) and where points g₂(x_m,y_m) are interpolated from points g₁(x_n,y_n) around points g₁(x_m, y_m).

[0017] Another object of the present invention is to provide an improved dual-energy calibration method so that the decomposed two materials are the actual breast tissue components, instead of merely equivalent materials. Currently, dual-energy calibration is performed by measuring x-ray attenuation curves using attenuation plates of a pair of different materials having known values of thickness. Generally, the materials used for calibration cannot be the same as those actually present in the human body, because the actual substances in the human body are too complex to be made as quantitative materials. For example, it is a common practice to use aluminum as a representation of human bone material and to use Lucite as a representation of human average soft tissue. Thus, the dual-energy results can only provide an equivalent aluminum quantity and an equivalent Lucite quantity in terms of x-ray attenuation in the human body. By using the standard calibration method and the equivalent decomposition, the second order approximation of the present invention cannot be achieved. Thus, the improved dual-energy calibration method is part of the present invention for a high accuracy decomposition of the human breast into pure single-component images.

[0018] The calibration method of the present invention includes the steps of (1) determining a dual-energy equation system for two known materials u and v as D_H=D_DH[u,v], D_L=D_DL[u,v], (2) conducting functional decomposition for energy-dependent attenuation coefficient functions μ_u(E)=u_p×μ_p(E)+u_q×μ_q(E) and μ_v(E)=v_p×μ_p(E)+v_q×μ_q(E) using standard least-square data fitting methods to obtain constants u_p, u_q, v_p, and v_q, (3) calculating u×(u_p+v_p)=p and v×(u_q+v_q)=q for each coordinate pair (u,v) to obtain dual-energy equation system D_H=D_DH[p,q], D_L=D_DL[p,q], and (4) solving the equation system D_H=D_DH[p,q], D_L=D_DL[p,q] for the materials p and q as a function of variable pair (D_H,D_L) through numerical inversion to obtain the equation system p=p[D_H,D_L] and q=q[D_H,D_L]. This calibration is performed for the component pair microcalcification and the soft tissue, which is an average combination of fat tissue and lean tissue, and is then performed for the component pair fat tissue and lean tissue.

[0019] Other objects of the present invention will become apparent in light of the following drawings and detailed description of the invention.

[0020] For a fuller understanding of the nature and object of the present invention, reference is made to the accompanying drawings, wherein:

[0021] FIG. 1 is a schematic presentation of the hardware system of the present invention;

[0022] FIG. 2 is a flow diagram of the method for performing a decomposition of a breast image into first order component images;

[0023] FIG. 3 is a flow diagram of the method for obtaining a relatively pure microcalcification image;

[0024] FIG. 4 is a flow diagram of the method for obtaining a microcalcification image based on the information provided by the fat and lean tissue images; and

[0025] FIG. 5 is a flow diagram of the method for further correction of the fat tissue and lean tissue images.

[0026] Hardware

[0027] The dual-energy breast x-ray imaging system of the present invention has two hardware configurations as described in the Chao disclosures, one described in U.S. Pat. Nos. 5,648,997 and 5,771,269 and the other described in U.S. patent application Ser. No. 09/025,926, all of which are incorporated herein by reference. The present invention utilizes these hardware configurations without change. The common features of the two configurations, as shown in FIG. 1, include a dual energy x-ray source 14, a front detector 16, a beam selector 18, and a rear detector 26. The subject under examination 12 is located between the x-ray source 14 and the front detector 16. The x-ray source 14 emits two consecutive pulses, a high-energy pulse at an average energy level H followed by a low-energy pulse at an average energy level L. Each pulse has a single energy spectrum, which is composed of bremsstrahlung radiation and discrete line emissions. A single energy x-ray spectrum generally cannot be represented by a single energy value, but has an distributed energy range. With the emission of each x-ray pulse, the front detector records a composite high-spatial-resolution image, the rear detector records one or two low-spatial-resolution images, depending upon the specific beam selector used.

[0028] A portion of the x-rays 32 passes through the subject 12 directly to the front detector assembly 16 without a change in their direction of propagation. These x-rays 32 are called the primary x-rays and convey true information about the subject 12. The remainder of the x-rays 34 are randomly scattered as a result of interaction with the material of the subject 12. These x-rays 34 are called scatter and cause a distortion of the true information.

[0029] There are two types of beam selectors: type I and type II beam selectors. The type I beam selector blocks all scatter, allowing only primary x-rays to be detected by the rear detector array at a number of selected locations (i,j). From a pair of dual-energy images recorded by the rear detector at the selected locations (i,j), the scatter images at the two energy levels for the front detector can be calculated. Then, a pair of high-spatial-resolution primary images at the front detector can be calculated. By using the dual-energy data decomposition method as described in U.S. Pat. Nos. 5,648,997 and 5,771,269, two material composition images can be calculated.

[0030] The type II beam selector blocks primary x-rays from a number of locations on the rear detector and allows both primary x-rays and scatter to other locations of the rear detector. From these two images, the dual-energy primary x-ray images at selected locations (i,j) on the rear detector can be calculated. Once these images are determined, the process is the same as that for the type I beam selector configuration.

[0031] Method

[0032] The method of the present invention comprises a number of parts. The first is a preliminary decomposition of a mixed-component breast x-ray image into five component images: the high-energy scatter image data, the low-energy scatter image data, the fat tissue image, the lean tissue image, and the microcalcification image. This part uses the scatter removal method and the dual-energy decomposition method as described in the Chao disclosures. However, the calibration method is new, and is an aspect of the present invention. The new calibration method allows the dual-energy decomposition to provide true breast tissue component images, instead of aluminum equivalent or Lucite equivalent material images. As a result, the five component images are a reasonably good representation of human breast information. However, important further improvements can still be made, because the human breast is composed of three tissue components, a simple dual-energy decomposition still has some limitations. Some of the component images obtained in the first part of the procedure still contain a small portion of mixed signals, and are not a pure single-component image. These preliminary image data in five image components is denoted as a first order approximation.

[0033] The second part of the method of the present invention provides more accurate single-component images for use in medical diagnosis. The result is the second order approximation.

[0034] First Order Approximation

[0035] The steps of the first order approximation are shown in the flow diagram of FIG. 2. The first step is to perform two calibrations for different material pairs under the exactly same system conditions.

[0036] The first calibration is to obtain a pair of numerical relationships for the front detector at the high energy level H and the low energy value L for two materials c and S, where c denotes the material having an x-ray attenuation coefficient as a function of x-ray energy μ_c(E) the same as that of microcalcification deposits, and s denotes the material having an x-ray attenuation coefficient as a function of x-ray energy μ_s(E) that same as that of an average tissue substance. An average breast soft tissue composition can be chosen as having 50% fat tissue and 50% lean tissue. When using other values of fat and lean tissue composition, the present invention holds equally true. The following is a description of how to perform dual-energy calibration for true tissue components.

[0037] Dual-energy calibration methods are described in detail in the Chao disclosures. In the present invention, it is required that the calibration be conducted regarding the actual material composition that exists in the human breasts. For conducting calibration measurements, according to the methods described in the Chao disclosures, a number of material plates of the two composition materials should be used in the calibration measurements. Currently, such standard tissue samples have not been available. The present invention provides an improved dual-energy calibration method that can use another pair of materials for calibration measurements to obtain the accurate calibration data for the desired different material pair. The scientific foundation of the improved calibration method is described below:

[0038] Assume that it is desired to obtain numerical equations for the calibration data material pair of (p,q), but the direct material (p,q) data is not available due to certain practical limitations. On the other hand, another material pair (u,v) is available for use in the calibration measurements. Also, all of the x-ray attenuation coefficients as x-ray energy functions μ_u(E), μ_v(E), μ_p(E) and μ_q(E) are known. In mammography, for one type of calibration, such a (p,q) material pair can be, for example, glandular tissue and adipose tissue in the human breast. The (u,v) material pair can be, for example, a pair of commercially available materials, such polyethylene (CH₂) sheet and water (H₂O). For another type of calibration, such a (p,q) material pair can be calcification and soft tissue in the human breast. Then the (u,v) material pair can be, for example, a pair of commercially available materials, such as aluminum (Al) and a plastic material, Lucite™. The dual-energy equation system is written as $\begin{array}{cc}\begin{array}{c}{D}_{\mathrm{fPHh}}\left(x,y\right)=\int [{\Phi }_{0H}\left(E\right)\times \exp (-({\mu }_u\left(E\right)\times u\left(x,y\right)+\\ {\mu }_v\left(E\right)\times v\left(x,y\right))]\times S_f\left(E\right)dE\end{array}\mathrm{and}& \text{(1a)}\\ \begin{array}{c}{D}_{\mathrm{fPLh}}\left(x,y\right)=\int [{\Phi }_{0L}\left(E\right)\times \exp (-({\mu }_u\left(E\right)\times u\left(x,y\right)+\\ {\mu }_v\left(E\right)\times v\left(x,y\right))]\times S_f\left(E\right)dE\end{array}& \text{(1b)}\end{array}$

[0039] where Φ_0H(E) and Φ_0L(E) are the energy spectra of the x-ray source 14 at the higher energy level H and the lower energy level L, respectively. The projection mass densities u(x,y) and v(x,y) are in units of gram/centimeter²(g/cm²). The mass attenuation coefficients μ_u(E) and μ_v(E) are expressed in units of centimeter²/gram (cm²/g). All of these values are known. S_f(E) is the x-ray spectral sensitivity of the front detector.

[0040] Conduct functional decomposition for the energy-dependent attenuation coefficient function μ_u(E)

[0041] μ_u(E)=u_p×μ_p(E)+u_q×μ_q(E)  (2a)

[0042] where the constants up and uq can be obtained through standard least-square data fitting methods. Similarly

[0043] μ_v(E)=v_p×μ_p(E)+v_q×μ_q(E)  (2b)

[0044] where the constant values u_pand u_qcan be obtained through standard least-square data fitting methods.

[0045] It should be noted that such functional decomposition cannot be generally true from a mathematics point of view. However, it has been established through extensive analysis of experimental data that for substances composed of chemical elements with low to medium atomic number Z, such as chemical elements carbon, hydrogen, oxygen, nitrogen, up to calcium, as in the human body, and as long as the x-ray energy is within the medical diagnostic energy range, such as between 20 keV and 150 keV, a third x-ray mass attenuation function can always be written as a sum of two other x-ray mass attenuation functions with good approximation. The energy-dependent attenuation function in the exponent of the exp( ) of equations (1a) and (1b) is denoted UF(E): $\begin{array}{cc}\begin{array}{c}-\mathrm{UF}\left(E\right)={\mu }_u\left(E\right)\times u\left(x,y\right)+{\mu }_v\left(E\right)\times v\left(x,y\right)\\ =\left[u\left(x,y\right)\times u_p+v\left(x,y\right)\times v_p\right]\times {\mu }_p\left(E\right)+\\ \left[u\left(x,y\right)\times u_q+v\left(x,y\right)\times v_q\right]\times {\mu }_q\left(E\right)\end{array}& \left(3\right)\end{array}$

[0046] When

[0047] u(x,y)×u_p+v(x,y)×v_p≡p(x,y)  (4a)

[0048] and

[0049] u(x,y)×u_q+v(x,y)×v_q−q(x,y)  (4b)

[0050] where the symbol “≡” represents definition, then

[0051] −UF(E)=μ_p(E)×p(x,y)+μ_q(E)×q(x,y)  (5)

[0052] Therefore, the procedure for the improved calibration includes the following steps:

[0053] (1) Conduct calibration measurements by using the (u,v) material pair as described in the Chao disclosures. As a result, this numerical equation system is obtained:

[0054] D_H=D_DH[u,v]  (6a)

[0055] D_L=D_DL[u,v]  (6b)

[0056] (2) Conduct functional decomposition for the energy dependent attenuation coefficient functions μ_u(E) and μ_v(E)

[0057] μ_u(E)=u_p×μ_p(E)+u_q×μ_q(E)  (7a)

[0058] μ_v(E)=v_p×μ_p(E)+v_q×μ_q(E)  (7b)

[0059] by using standard least-square data fitting methods, and obtain the constants u_p, u_q, v_p, and v_q.

[0060] (3) For each data pair (u,v) in equations (6a),(6b), calculate

[0061] u×u_p+v×v_p=p  (8a)

[0062] u×u_q+v×v_q=q  (8b)

[0063] and obtain a new dual-energy numerical equation system

[0064] D_H=D_DH[p,q]  (9a)

[0065] D_L=D_DL[p,q]  (9b)

[0066] (4) Solve the equation system D_H=D_DH[p,q] and D_L=D_DL[p,q] for the variable pair (p,q) as a function of the variable pair (D_H,D_L) through numeric inversion and obtain

[0067] p=p[D_H,D_L]  (10a)

[0068] q=q[D_H,D_L]  (10b)

[0069] For performing two types of dual-energy decomposition, two corresponding types of calibration are conducted. The first type of calibration is for the material pair (c,s) for microcalcification and soft tissue in the human breast. By using the above-described calibration procedures, a pair of quantitative explicit numerical functions

[0070] D_H=D_DH[c,s]  (11a)

[0071] D_L=D_DL[c,s]  (11b)

[0072] are established, where the notations D_H[ ] and D_L[ ] represent functional relationships with c and s as variables. The equation system is solved by, for example, numerical inversion to obtain a numerical equation system

[0073] c=c[D_H,D_L]  (12a)

[0074] s=s[D_H,D_L]  (12b)

[0075] This numerical equation system is applicable to all normalized pixels. From the measured signal pair (D_H,D_L) for any pixel, the material composition data pair c and s for that pixel can be readily found.

[0076] The second type of calibration is to obtain a pair of numerical relationships for the front detector at the high energy level H and the low energy value L for two materials f and g, where f denotes the material having an x-ray attenuation coefficient μ_f(E) the same as that of fat tissue, and g denotes the material having an x-ray attenuation coefficient μ_g(E) the same as that of lean tissue. Thus, the quantitative explicit functions

[0077] D_H=D_DH[f,g]  (13a)

[0078] D_L=D_DL[f,g]  (13b)

[0079] are established. The equation systems is solved by, for example, numerical inversion to obtain a numerical equation system

[0080] f=f[D_H,D_L]  (14a)

[0081] g=g[D_H,D_L]  (14b)

[0082] This numerical equation system is applicable to all normalized pixels. From the measured signal pair (D_H,D_L), the material composition data pair (f,g) can be readily found.

[0083] To avoid misunderstanding, two points need to be mentioned. First, for better accuracy, the calibration materials should be chosen to have a composition as close as possible to the actual tissue compositions in terms of x-ray interaction. This should be clear based on general error theory for experimental measurements. For example, when μ_u(E) is close to μ_p(E), and μ_v(E) is close to μ_q(E), the error due to transformation of (7a) and (7b) is smaller.

[0084] The second point is that there exists a number of trivial alternatives of the order of the calculation steps. For example, after step (1), conducting calibration measurements by using the (u,v) material pair and establishing the numerical relationships (6a) and (6b), instead of performing a material pair transformation (8a) and (8b) to obtain the numerical relationships (9a) and (9b), one can directly conduct the inversion for the (6a) and (6b) to obtain

[0085] u=u[D_H,D_L]  (15a)

[0086] v=v[D_H,D_L]  (15b)

[0087] Then, applying the relationships (15a) and (15b) to the image data, D_fPHh(x,y) and D_fPLh(x,y), two material images u(x,y) and v(x,y) can be obtained. After that, by using the equations (8a) and (8b)

[0088] p(x,y)=u(x,y)×u_p+v(x,y)×v_p  (16a)

[0089] q(x,y)=u(x,y)×u_q+v(x,y)×v_q  (16b)

[0090] the entire images u(x,y) and v(x,y) are transformed, point by point, into two images p(x,y) and q(x,y). This is only a trivial variation of the order for mathematical calculations without any essential change of the method. Therefore, it is contemplated that the procedures as described in equations (1a),(1b) through (12a),(12b) include these and other similar trivial alternatives.

[0091] The second step of the method of the present invention is to acquire a pair of high-spatial-resolution composite images of the front detector D_fHh(x,y) and D_fLh(x,y), where the subscript f denotes an image of the front detector 16, subscript H denotes a high-energy-level image, subscript L denotes a low-energy-level image, subscript h denotes a high-resolution image, and (x,y) denotes the two-dimensional Cartesian coordinates of a detector cell of the front detector 16. Each of the images in image data pair D_fHh(x,y) and D_fLh(x,y) contains a scatter component and a primary x-ray component. The primary x-ray component contains a fat tissue component, a lean tissue component, and a microcalcification component.

[0092] The third step is to separate each image of the image air D_fHh(x,y) and D_fLh(x,y) into the scatter component and primary x-ray component. This is done through a number of data acquisition and calculation steps, as described in the above-identified U.S. patents, to calculate a pair of high-spatial-resolution scatter-only images at the front detector D_fSHh(x,y) and D_fSLh(x,y), where the subscript S denotes that this is a scatter-only image. This step is the same as the corresponding step in the previously-identified patents.

[0093] The fourth step is to calculate the high-spatial-resolution primary x-ray image pair D_fPHh(x,y) and D_fPLh(x,y) from the equations

[0094] D_fPHh(x,y)=D_fHh(x,y)−D_fSHh(x,y)  (17a)

[0095] D_fPLh(x,y)=D_fLh(x,y)−D_fSLh(x,y)  (17b)

[0096] where the subscript P denotes that this image is composed only of primary x-rays. This step is also the same as described in the U.S. patents.

[0097] The fifth step is to perform a dual-energy decomposition for the image pair D_fPHh(x,y) and D_fPLh(x,y) by using the calibrated numerical equation system (12a),(12b). As a result, two material composition images c₁(x,y) and s₁(x,y) are obtained. In this image pair, c₁(x,y) is basically a microcalcification image and s₁(x,y) is basically an average soft tissue image of the breast.

[0098] The sixth and final step is to perform a dual-energy decomposition for the image pair D_fPHh(x,y) and D_fPLh(x,y) by using the calibrated numerical equation system (14a),(14b). As a result, two material composition images f₁(x,y) and g₁(x,y) are obtained. In this image pair, f₁(x,y) is basically a fat tissue image and g₁(x,y) is basically a lean tissue image of the breast. However, both f₁(x,y) and g₁(x,y) contain microcalcification information, which has not been explicitly specified and appropriately quantified.

[0099] At this point, usable results have been obtained by the decomposition of a single breast image into five component images: the scatter image D_fs(x,y), the microcalcification image c₁(x,y), the average mixed tissue image s₁(x,y), the fat tissue image b₁(x,y), and the lean tissue image g₁(x,y). Using these component images, information needed to make a diagnosis is significantly enhanced over that of the prior art.

[0100] However, the power of dual-energy decomposition has not yet been fully exploited. The above results may be considered as a good first approximation for breast imaging data. The reason that these component images can only be considered as first approximation is because that, in performing the above procedures of dual-energy decomposition, the breast is assumed to be composed of only two materials. A real breast is composed of three materials.

[0101] Second Order Approximation

[0102] The next part of the present specification describes a method of second order approximation for taking into account the three material components of a human breast. First, the scientific principles are explained, and then the steps of the method are described.

[0103] In a human breast, there are basically three significant material compositions: adipose tissue, which can be represented by a fat tissue, glandular duct tissue, which can be represented by a lean tissue, and microcalcification deposits, which contain calcium. Decomposition of a system into three independent material components is generally beyond the capability of dual-energy imaging method. If there existed a triple-energy decomposition method, performing data acquisition at three different energy levels, then three material compositions could be independently separated. However, there currently is no triple-energy decomposition method.

[0104] The following explains how to use the dual-energy method to approximately separate the three material components of a human breast. This is possible because there are special conditions of the human breast that can be utilized to achieve a good approximation. The first special condition is that the x-ray attenuation of fat tissue and lean tissue in the breast tissue are very close when compared to that of the microcalcification. The second condition is that microcalcification deposits in the human breast are generally limited to existing in small, isolated regions, and in many cases, to be within only few single image pixels.

[0105] The dual-energy imaging equation of the breast can be represented by the equation pair $\begin{array}{cc}\begin{array}{c}{D}_{\mathrm{fPHf}}\left(x,y\right)=\int [{\Phi }_{0H}\left(E\right)\times \exp (-({\mu }_c\left(E\right)\times c\left(x,y\right)+\\ {\mu }_f\left(E\right)\times f\left(x,y\right)+{\mu }_g\left(E\right)\times g\left(x,y\right))]\times \\ S_f\left(E\right)dE\end{array}\mathrm{and}& \text{(18a)}\\ \begin{array}{c}{D}_{\mathrm{fPLh}}\left(x,y\right)=\int [{\Phi }_{0L}\left(E\right)\times \exp (-({\mu }_c\left(E\right)\times c\left(x,y\right)+\\ {\mu }_f\left(E\right)\times f\left(x,y\right)+{\mu }_g\left(E\right)\times g\left(x,y\right))]\times \\ S_f\left(E\right)dE\end{array}& \text{(18b)}\end{array}$

[0106] The projection mass densities c(x,y), f(x,y), and g(x,y) of the subject are in units of gram/centimeter²(g/cm²). The mass attenuation coefficients μ_c(E), μ_f(E) , and μ_g(E are known and expressed in units of centimeter²/gram (cm²/g). The energy-dependent exponent of the exponential function exp( ) in equation system (18a), (18b) is denoted as AF(E), and is $\begin{array}{cc}\begin{array}{c}-\mathrm{AF}\left(E\right)={\mu }_c\left(E\right)\times c\left(x,y\right)+{\mu }_f\left(E\right)\times f\left(x,y\right)+\\ {\mu }_g\left(E\right)\times g\left(x,y\right)\end{array}& \left(19\right)\end{array}$

[0107] There are three independent mass attenuation functions μ_c(E), μ_f(E), and μ_g(E) connected by three unknown parameters for material density values c(x,y), f(x,y), and g(x,y). Dual-energy measurements only provide two measured signals D_fPHh(x,y) and D_fPLh(x,y), from which only two unknown parameters can be determined.

[0108] The basic method used for quantitative evaluation of each tissue component is a method of function decomposition. One energy-dependent mass attenuation function must be expressed as a sum of the other two functions. For example, to perform a basic dual-energy decomposition for fat tissue and lean tissue, the calcium function μ_c(E) is decomposed as

[0109] μ_c(E)=R_cf×μ_f(E)+R_cg×μ_g(E)  (20)

[0110] where R_cfand R_cgare constants. Note that for the human body system, where the x-ray energy is within medical diagnostic energy range, the third x-ray mass attenuation function can always be written as a sum of other two x-ray mass attenuation functions with good approximation. The exact values of R_cfand R_cgcan be calculated by using standard least-square data fitting methods. For example, within an x-ray energy range of 16 keV and 100 keV, by applying least-square data fitting methods to equation (20), it can be found that R_cf{tilde over (=)}−34 and R_cg{tilde over (=)}+33. Accurate values of the R_cfand R_cgare dependent on the x-ray energy range used in the hardware system.

[0111] Some general observations can be made regarding the parameter pair R_cfand R_cg. First, R_cgis always positive, and R_cfis always negative. Secondly, the absolute values of R_cfand R_cgare always rather large and are always close to each other. For example, for the diagnostic energy range, the absolute value of R_cfand R_cgcan be assumed to be between 10 to 50. The underlying physical reason for these two characteristic features comes from the fact that calcium has a much larger attenuation coefficient at lower energy than both fat tissue and lean tissue, and at the low energy, the lean tissue has a slightly larger attenuation coefficient than that of fat. These characteristic features produce very uncommon results in the decomposed component x-ray images: when the image is decomposed into fat tissue and lean tissue component images, each microcalcification point produces a positive high-intensity attenuation change in the lean tissue image at that point, each microcalcification point produces a negative high-intensity attenuation change in the fat tissue image at that point, and these positive and negative points have an exact quantitative relationship. For example, one microgram of microcalcification produces a positive attenuation change corresponding to approximately 30 micrograms of lean tissue intensity change, and one microgram of microcalcification produces a negative attenuation point corresponding to approximately 30 micrograms of fat tissue intensity change. These characteristic features are utilized in the present invention for quantitatively separating the superimposed microcalcification image from the lean tissue image and from the fat tissue image.

[0112] AF(E) can be again written in a form that contains only two basis functions $\begin{array}{cc}\begin{array}{c}-\mathrm{AF}\left(E\right)={\mu }_f\left(E\right)\times \left[f\left(x,y\right)+R_{\mathrm{cf}}\times c\left(x,y\right)\right]+\\ {\mu }_g\left(E\right)\times \left[g\left(x,y\right)+R_{\mathrm{cg}}\times c\left(x,y\right)\right]\end{array}& \left(21\right)\end{array}$

[0113] or

[0114] −AF(E)=μ_f(E)×fc(x,y)+μ_f(E)×gc(x,y)  (22)

[0115] where

[0116] fc(x,y)=f(x,y)+R_cf×c(x,y)  (23a)

[0117] gc(x,y)=g(x,y)+R_cg×c(x,y)  (23b)

[0118] Because AF(E) contains only two independent functions μ_f(E) and μ_g(E), by using dual energy decomposition procedures as described in the Chao disclosures, two unique energy-independent parameters can be obtained through solving equation system (18a) and (18b), the solution pair is fc(x,y) and gc(x,y). Note that the solution pair fc(x,y) and gc(x,y) are exactly the solution pair of the first approximation. That is

[0119] f₁(x,y)=fc(x,y)=f(x,y)+R_cf×c(x,y)  (24a)

[0120] g₁(x,y)=gc(x,y)=g(x,y)+R_cg×c(x,y)  (24b)

[0121] After performing dual-energy decomposition according to the procedures of the first order approximation, the third material component is projected and superimposed onto the two basis component images as expressed in (24a) and (24b), as R_cf×c(x,y) on the fat tissue image and as R_cg×c(x,y) on the lean tissue image. Generally speaking, the portion of superimposed third component image cannot be separated from the two component images because what can be calculated by using dual-energy decomposition are only fc(x,y) and gc(x,y).

[0122] The present invention utilizes certain characteristic specific conditions for human breast imaging to approximately separate the third material composition. It is known that, in the human breast, the microcalcification deposits are limited to existing in small regions. For a small region, the soft tissue density can be assumed to change approximately smoothly. The following is a outline of the basic steps of separating microcalcification image spots from the soft tissue image. First, identify the microcalcification spots on the image. The microcalcification spots must be (a) a positive signal/negative signal pair, each having approximately equal absolute values of amplitude, but having opposite signs, and (b) superimposed at exactly the same corresponding locations on the lean tissue image and the fat tissue image. Second, interpolated an image intensity from nearby surrounding soft tissue pixels to replace the image intensity of these suspected microcalcification spots. Third, subtract the interpolated signal intensity from the original signal intensity gc(x,y) and fc(x,y) at the identified spot to determine if the intensity change is consistent with that predicted by decomposition equation (20) for the microcalcification function. If the microcalcification intensity on the two soft tissue images are quantitatively consistent with that predicted by the equation (20), the spot meets the preliminary test for a microcalcification spot. In addition to this preliminary test, there is a further test, to be described below.

[0123] The following is a description of the scientific foundation of how to obtain pure microcalcification single component image through the second decomposition.

[0124] In the second decomposition at the accuracy level of first order approximation, the two basis functions are chosen to be μ_c(E) and μ_sa(E), where μ_c(E) is the energy-dependent attenuation function for microcalcification and μ_sa(E) is the attenuation function for an average soft tissue composition. Two images c₁(x,y) and s₁(x,y) are obtained as a result of the second decomposition at an accuracy level of first order approximation. In the first order approximation, the breast was assumed to have only two material compositions with the mass attenuation coefficient functions μ_c(E) and μ_sa(E). According to the actual breast material composition, the energy-dependent exponent in the equations (18a) and (18b) is denoted as BF(E): $\begin{array}{cc}\begin{array}{c}-\mathrm{BF}\left(E\right)={\mu }_c\left(E\right)\times c\left(x,y\right)+{\mu }_f\left(E\right)\times f\left(x,y\right)+\\ {\mu }_g\left(E\right)\times g\left(x,y\right)\end{array}\mathrm{or}& \left(25\right)\\ \begin{array}{c}-\mathrm{BF}\left(E\right)={\mu }_c\left(E\right)\times c\left(x,y\right)+{\mu }_{\mathrm{sa}}\left(E\right)\times s_a\left(x,y\right)+\\ {\mathrm{\delta \mu }}_s\left(E\right)\times \delta s\left(x,y\right)\end{array}& \left(26\right)\end{array}$

[0125] where the average tissue mass attenuation coefficient function μ_sa(E)=50% μ_f(E)+50% μ_g(E). δμ_s(E)×δs(x,y) is the local deviation at the pixel (x,y) from the averaged soft tissue attenuation μ_sa(E)×s_a(x,y). For example, when the tissue composition at (x,y) is exactly 50% fat tissue and 50% lean tissue, then δμ_s(E)=0. Otherwise, for performing dual-energy decomposition according to the basis function μ_c(E) and μ_sa(E), the energy dependent function δμ_s(E) must be decomposed into two basis functions μ_c(E) and μ_sa(E) as

[0126] δμ_s(E)=R_sc×μ_c(E)+R_ss×μ_sa(E)  (27)

[0127] where R_scand R_ssare two constants, the exact values of which can be calculated by using standard least-square data fitting methods, when the local tissue composition is known. Then, $\begin{array}{cc}\begin{array}{c}-\mathrm{BF}\left(E\right)={\mu }_{\mathrm{sa}}\left(E\right)\times \left[s_a\left(x,y\right)+R_{\mathrm{ss}}\times \delta s\left(x,y\right)\right]+\\ {\mu }_c\left(E\right)\times \left[c\left(x,y\right)+R_{\mathrm{sc}}\times \delta s\left(x,y\right)\right]\end{array}& \left(28\right)\end{array}$

[0128] or

[0129] −BF(E)=μ_c(E)×cs(x,y)+μ_sa(E)×ss(x,y)  (29)

[0130] where

[0131] cs(x,y)=c(x,y)+R_sc×δs(x,y)  (30a)

[0132] and

[0133] ss(x,y)=s_a(x,y)+R_ss×δs(x,y)  (30b)

[0134] Some general observations can be made regarding the parameter pair R_ssand R_sc. First, R_ssis not much different from 1.0 and R_scis much smaller than 1.0. The underlying physical reason stems from the fact that the energy-dependent behavior of the function δμ_s(E) is not much different from that of the average tissue function μ_sa(E), and very much different from that of μ_c(E).

[0135] Note that the dual-energy decomposition results obtained at the accuracy level of the first order approximation c₁(x,y) is exactly cs(x,y) and that s₁(x,y) is exactly ss(x,y). That is,

[0136] c₁(x,y)=cs(x,y)=c(x,y)+R_sc×δs(x,y)  (31a)

[0137] s₁(x,y)=ss(x,y)=s_a(x,y)+R_ss×δs(x,y)  (31b)

[0138] Therefore, because of the local deviation of the actual tissue composition from the assumed average tissue composition, the microcalcification image c₁(x,y) obtained through the first order approximation is composed of two images. One is the true microcalcification image c(x,y), and the other image is produced by local deviation of the tissue composition from the assumed average composition. The exact signal intensity of the superimposed image is R_sc×δs(x,y). Generally speaking, such superimposed images due to the existence of a third material cannot be separated by using the dual-energy method. However, based on utilizing some specific conditions for the human breast, the present invention provides approximation methods for removing the effects of the third material R_sc×δs(x,y) on the true microcalcification image c(x,y).

[0139] For separating the third material effects on the microcalcification image, two characteristic feature differences between the image R_sc×δs(x,y) and that of image c(x,y) are utilized: (a) the local deviation of soft tissue composition from the average composition R_sc×δs(x,y) has a smooth spatial distribution, in a way similar to that of the soft tissue spatial distribution sa(x,y), while the microcalcification image has an abrupt intensity change at certain small regions; and (b) the amplitude of deviation image R_sc×δs(x,y) is small, because, as mentioned above, R_scis small.

[0140] FIG. 3 and the following describe the procedure for obtaining a pure microcalcification image at an accuracy level of the second order approximation from the second decomposition.

[0141] (1) In the microcalcification image c₁(x,y), first identify all the microcalcification points (x_k,y_k). The microcalcification points are prominently different from other image points because the microcalcification image c₁(x,y) is composed essentially of high-intensity isolated image points presented on a smoothly changing background. The microcalcification points are identified by pixels with a sudden intensity change at isolated points. Sometimes, microcalcification can take the form of clusters, so the points are not always individually isolated. The present invention contemplates that such clusters containing a small number of microcalcification points with a signal intensity much higher than those of adjacent points exist.

[0142] (2) Construct a background image B(x,y) for c₁(x,y) from pixels (x_i,y_i)≠(x_k,y_k). That is, background image B(x_i,y_i)=c₁(x_i,y_i) for x_i≠x_kand y_i≠y_kfor all the non-microcalcification points. The signals in B(x_i,y_i) have excluded the microcalcification points, so they must be produced by local tissue composition deviation, hence are microcalcification background signals.

[0143] (3) In the background image B(x,y) for those pixels (x_k,y_k) that have been identified as microcalcification points, the signal intensity is assumed to be equal to an averaged signal intensity value interpolated from the signal intensity of the surrounding pixels. That is,

[0144] B(x_k,y_k)=averaged[c₁(x_j,y_j)]  (32)

[0145] for j around k. Because the tissue composition has a relatively smooth change in comparison with that of the microcalcification deposits, the averaged c₁(x_j,y_j) will be very close to the actual tissue composition variation at the point (x_k,y_k).

[0146] (4) Subtract the tissue background image B(x,y) from the first approximation microcalcification image c₁(x,y) to obtain the second order approximation microcalcification image

[0147] c₂(x,y)=c₁(x,y)−B(x,y)  (33)

[0148] In the image c₂(x,y), the image signal intensity should accurately represent quantitative microcalcification density without being affected by tissue composition variations.

[0149] The steps for construction of a microcalcification image c_fg(x,y) based on the information in the first order approximation fat and lean tissue images f₁(x,y) and g₁(x,y) are described in FIG. 4 and are as follows:

[0150] (1) Identify those points (x_m,y_m) where an abrupt intensity change in relation to those of adjacent points takes place simultaneously in the fat tissue image f₁(x,y) and the lean tissue image g₁(x,y). The points (x_m,y_m) are suspected calcification sites as given by the tissue images f₁(x,y) and g₁(x,y).

[0151] (2) Verify whether the abrupt signal change at each point (x_m,y_m) is consistent with that predicted. Only when the predicted signal change on the fat tissue image f₁(x,y) and lean tissue image g₁(x,y) can both be quantitatively confirmed is the new microcalcification image signal c_fg(x,y) constructed. To do this, first calculate a pair of interpolated signal intensities f_1a(x_m,y_m) and g_1a(x_m,ym), so that f_1a(x_m,y_m)=averaged[f₁(x_n,y_n)] for n around m, and g_1a(x_m,y_m)=averaged[g₁(x_n,yn)] for n around m. Then calculate

[0152] Δf=−[f₁(x_m,y_m)−f_1a(x_m,y_m)]  (34)

[0153] and

[0154] Δg=g_1a(x_m,y_m)−g₁(x_m,y_m)  (35)

[0155] When these two quantities are close to being within a certain range, for example, between 20 to 40, the microcalcification at the pixel (x_m,y_m) is confirmed. When Δf and Δg are not quantitatively consistent, the most probable reasons are that (a) the point selected according to the criteria of abrupt intensity change is not a calcification deposit, but actually represents some peculiar soft tissue structure, or (b) the point could contain calcification, but the signal intensity of the calcification is too low and submerged in noise. In both cases, the selected point cannot be taken as a calcification deposit, and should be excluded from the calcification image according to the information provided by the fat tissue image f₁(x,y) and lean image g₁(x,y). Thus, the calcification image c_fg(x,y) may not be entirely identical to the calcification image c₂(x,y). Even with these possible slight differences between the two calcification images obtained from different data analysis methods, both c₂(x,y) and c_fg(x,y) provide their own reference data. Both images are useful.

[0156] (3) According to the results of step (2), let

[0157] c_fg(x_m,y_m)=−[f₂(x_m,y_m)−f₁(x_m,y_m)]×R_cf  (36)

[0158] or

[0159] c_fg(x_m,y_m)=[g₂(x_m,y_m)−g₁(x_m,y_m)]×R_cg  (37)

[0160] Either value of equations (36) or (37) can be used as the calcification image intensity, because these values are close to each other. Choosing one as the final result is largely a matter of operator preference and depends on the operator's judgement. If a further comparison is needed, the value of equation (36) may be better than that of equation (37) because, generally, the signal intensity of fat tissue in the breast has a smoother distribution than that of lean tissue. Consequently, the interpolated and averaged value may be more accurate.

[0161] The steps for obtaining a second order approximation of the fat tissue image f₂(x,y) and lean tissue image g₂(x,y) in which the microcalcification image is removed are described in FIG. 5 and are as follows:

[0162] (1) For each pixel of the second order microcalcification image c_fg(x,y), identify those points c_fg(x_o,y_o) whose intensity is zero and those points c_fg(x_m,y_m) whose intensity is not zero.

[0163] (2) Construct images f₂(x,y) and g₂(x,y), where there are no changes between the first order and second order images at the zero intensity points, that is, f₂(x_o,y_o)=f₁(x_o,y_o) and g₂(x_o,y_o)=g₁(x_o,y_o)

[0164] (3) Set the remainder of the points f₂(x_m,y_m) and g₂(x_m,y_m) to the average of the signal intensity of the points in the vicinity, that is, f₂(x_m,y_m)=averaged[f₁(x_n,y_n)] for n around m and g₂(x_m,y_m)=averaged[g₁(x_n,y_n)] for n around m. After going over all pixels (x,y), two corrected tissue images f₂(x,y) and g₂(x,y) are obtained with microcalcification effects corrected.

[0165] Thus it has been shown and described a method of digital mammography which satisfies the objects set forth above.

[0166] Since certain changes may be made in the present disclosure without departing from the scope of the present invention, it is intended that all matter described in the foregoing specification and shown in the accompanying drawings be interpreted as illustrative and not in a limiting sense.