Title of Invention	"AN APPARATUS AND A METHOD FOR SCATTER CORRECTION IN PROJECTION RADIOGRAPHY"
Abstract	The invention relates to an apparatus for projection radiography, in particular for mammography, comprising a radiation source (2) which emits radiation (3), a detector (8) and, downstream of said detector (8), a processing unit (12) which determines in an approximately manner, on the basis of the projection data (17) supplied by the detector (8), the scatter material distribution of the object under examination (6) and reads out scatter information (20) from a data memory (19) as a function of the scatter material distribution and corrects the projection data (17) in respect of the scatter content (11) on the basis of the scatter information (20), the scatter information (20) is determined by Monte Carlo simulations which calculate the interactions of the photons with a scatter material distribution in respect of different scatter material distributions, in that the scatter information is scatter distributions (20) which describes a scatter-induced distribution over adjacent image regions of the radiation (14) emitted by the radiation source (2) and directed to a particular image region, and in that the processing unit (12) determines a scatter distribution (11) in an image region of the projection image (17) by calculating and adding, for each image region, the scatter contributions (11) of the surrounding image regions.

Title of Invention

"AN APPARATUS AND A METHOD FOR SCATTER CORRECTION IN PROJECTION RADIOGRAPHY"

Abstract

The invention relates to an apparatus for projection radiography, in particular for mammography, comprising a radiation source (2) which emits radiation (3), a detector (8) and, downstream of said detector (8), a processing unit (12) which determines in an approximately manner, on the basis of the projection data (17) supplied by the detector (8), the scatter material distribution of the object under examination (6) and reads out scatter information (20) from a data memory (19) as a function of the scatter material distribution and corrects the projection data (17) in respect of the scatter content (11) on the basis of the scatter information (20), the scatter information (20) is determined by Monte Carlo simulations which calculate the interactions of the photons with a scatter material distribution in respect of different scatter material distributions, in that the scatter information is scatter distributions (20) which describes a scatter-induced distribution over adjacent image regions of the radiation (14) emitted by the radiation source (2) and directed to a particular image region, and in that the processing unit (12) determines a scatter distribution (11) in an image region of the projection image (17) by calculating and adding, for each image region, the scatter contributions (11) of the surrounding image regions.

Full Text	Description Apparatus and methods for scatter correction in projection radiography, particularly mammography The invention relates to an apparatus for projection radiography comprising a radiation source, a detector and, downstream of said detector, a processing unit which uses the projection data supplied by the detector to approximately determine the scatter material distribution of the object under examination and which reads out scatter information from a data memory as a function of the scatter material distribution and corrects the projection data in respect of the scatter component on the basis of said scatter information. The invention further relates to a method with scatter correction for projection radiography and a method for obtaining scatter information. An apparatus and methods of this kind are known from US 6,104,777 A. The scatter produced in the object under examination (breast), the intensity of which may almost attain the order of magnitude of the direct, unscattered, image-producing primary radiation, results in image quality impairment by reducing the contrast, increasing the noise, and ultimately in respect of the image post-processing methods used for differentiating between various types of tissue in the images produced, in particular between glandular and fatty tissue in the breast. For differentiation in respect of two tissue types, techniques using a single energy spectrum, i.e. a single x-ray tube voltage, or the dual energy method using two voltage values are known in mammography. In both cases scatter compensation is required; with the dual energy method this is also because the amount of scatter is different in the two energy spectra. To reduce scatter, mechanical measures have already been proposed. The use of slit collimators requires mechanical displacement of the slit collimators over the breast measuring field and is therefore time-consuming. Anti-scatter grids not only reduce the scatter, but also the imaging- producing primary radiation. Arguments in favor of dispensing with anti-scatter grids have been ongoing for years. For compression thicknesses of less than 4-5 cm, the dose could even be reduced or the SNR (= signal-to-noise ratio) increased if the grid is removed. On the other hand there are applications in which it is not technically possible to use a grid, e.g. in tomosynthesis. A large number of computer correction methods have already been proposed. Methods of this kind are known e.g. from M. DARBOUX, J.M. DINTEN: Physical model based scatter correction in mammography. In: Proc. SPIE, Vol. 3032, 1997, pages 405 to 410 and from J.M. DINTEN and J.M. VOLLE: Physical model based restoration of mammographies. In. Proc. SPIE, Vol. 3336, 1998, pages 641 to 650 and in US 6,104,777 A. These are convolution/deconvolution methods in which a scatter intensity distribution is approximated as a convolution of the primary radiation distribution using suitable convolution kernels . In the cited documents an analytical model is thus proposed with which the physical scatter process in the scatter object (breast) is explicitly computed as an integral transformation. However, this explicit analytical representation only describes first-order scatter, not multiple scatter. The intensity distribution of multiply scattered photons is assumed to be a spatially constant background over the detector surface and must be estimated from previously determined tables. The analytical model for calculating just the first-order scatter contribution requires 4-dimensional numerical integrations (3 space coordinates + energy spectrum) for each detector pixel, i.e. it is compute-intensive. Approximations are therefore required in order to reduce the computational complexity. Because of the high computational complexity it is proposed to perform the calculations in advance and tabulate the results. In addition, W. KALENDER: Monte Carlo calculations of x-ray scatter data for diagnostic radiology. In: Phys. Med. Biol., 1981, Vol. 26, No. 5, pages 835 to 849 describes the use of Monte Carlo methods for simulating radiation propagation in radiography. Proceeding from this prior art, the object of the invention is therefore to specify an apparatus and methods enabling improved scatter correction compared to the prior art to be performed. This object is achieved by an apparatus and the methods having the features set forth in the independent claims. Advantageous further developments and embodiments are set forth in claims dependent thereon. With the apparatus and method, the projection images supplied by a detector are analyzed in a processing unit. It is first attempted to approximately determine the scatter material distribution, in mammography typically the proportions of glandular and fatty tissue, of an object under examination. In another processing step, scatter information depending on the scatter material distribution is read out of a data memory. This scatter information can then be used to correct the projection images in respect of the scatter content of said projection images, it being essential that the scatter information read out of the data memory has been determined in advance by a Monte Carlo simulation which takes the multiple interaction of the photons with the object under examination into account. The basis for the solution described here is optimally correct physical modeling. In contrast to the prior art, modeling is possible which takes a much larger number of details into account, namely in the following respect: the occurrence of multiple scatter and polychromasia and the geometrical relationships, in particular the peculiarities of the scatter distribution at the edges of the objects, can be simulated. During the scatter correction itself, merely table access is required, possibly with subsequent interpolation, and the calculation of the scatter distribution in the detector plane is reduced to 2-dimensional integrations over the detector plane. Despite the relatively simple implementation of the scatter correction, the procedure is not limited to special cases and requires no drastic simplifications or approximations, such as simplified acquisition geometry, monochromatized radiation, simplifications of the physical model or a Taylor development by orders of approximation or similar. In a preferred embodiment, the scatter material distribution is specifically determined for different regions of the projection image. To perform scatter correction in an image region, the scatter contributions of the surrounding contributions which depend on the specific scatter material distribution are then determined and corrected accordingly. In this way it is possible to take local scatter variations into account. In another preferred embodiment, in the region of the edges of the object under examination, scatter information which takes the particular geometrical relationships in the region of the object edge into account is used for scatter correction. The scatter information is preferably obtained under the assumption that the scatter material distribution is homogeneous along the radiation direction. Particularly in the context of mammography, such an assumption only results in slight deviations from the actual scatter distribution. The specific scatter information associated with an image area can also be calculated under the assumption that the object under examination is also homogeneously structured in the direction perpendicular to the beam, thereby simplifying the calculation of the scatter information. However, if particularly high accuracy in calculating the scatter information is required, any inhomogeneity perpendicular to the beam direction can be taken into account. In a preferred embodiment, the scatter material distribution is determined by analyzing the ratio of incident radiation intensity to the unscattered primary radiation in an image region, the values for the primary radiation being ascertained by means of a scatter correction based on the scatter information associated with a characteristic homogeneous scatter material distribution. The processing steps carried out by the processing unit can also be executed iteratively, the calculated primary radiation components being used to simplify the approximate calculation of the scatter components and thus arrive at improved values for the primary radiation. The scatter correction does not generally need to be performed at full detector resolution. It may occasionally be sufficient to perform scatter correction at selected grid points and interpolate between the determined scatter correction values at the selected grid points. Further advantages and embodiments of the invention will emerge from the following description in which examples of the invention are specifically explained with reference to the accompanying drawings in which: Figure 1 shows the configuration of a mammography machine in which a breast under examination is compressed between two compression plates and irradiated with x-radiation; Figure 2 illustrates a simplified breast structure assumed for calculating the scatter correction; Figure 3 shows a flowchart of a method performed for scatter correction; Figure 4 shows the breast tissue distribution assumed for calculating a simple scatter beam spread function; and Figure 5 shows the breast structure assumed for calculating a precise scatter beam spread function. Figure 1 shows the configuration of a mammography machine 1 in which x-radiation 3 is produced with the aid of a radiation source 2. The divergence of the x-radiation 3 is limited if necessary using a collimator 4 which is shown in Figure 1 as a single beam diaphragm. However, the collimator 4 can also be contrived such that a plurality of virtually parallel x-ray beams is produced. Such a collimator 4 can be implemented e.g. as an iris diaphragm. The mammography machine 1 additionally has compression plates 5 between which a breast 6 is compressed. The x-radiation 3 passes through the compression plates 5 and the breast 6 and generally crosses an air gap 7 before the x-radiation 3 is incident on an x-ray detector 8 comprising a plurality of individual detector elements 9, the so-called detector pixels. The portion of x-radiation 3 passing through the breast 6 without interacting with the breast 6 is known as the primary radiation 10. On the other hand, the portions of x-radiation 3 incident on the x-ray detector 8 after at least one scattering within the breast 6 are referred to as secondary radiation 11. It should be noted that the term scatter is to be understood as any kind of interaction between the x-radiation 3 and the material of the breast 6 causing a change in the propagation direction of the photons of the x-radiation 3. Since, as mentioned in the introduction, the secondary radiation 11 may considerably distort the structure of the breast 6 imaged by the primary radiation 10, it is advantageous if the secondary radiation 11 can be removed from the projection images of the breast 6 captured by the x- ray detector 8. For this purpose, a processing unit 12 connected downstream of the x-ray detector 8 performs a scatter correction. In order to be able to perform the scatter correction, model assumptions are made concerning the structure of the breast 6. In particular it is assumed that the tissue structure of the breast 6 which is essentially composed of glandular and fatty tissue can be described by a homogeneous tissue distribution along the propagation direction of the x-radiation 3. Figure 2 accordingly shows different regions 13, 14 and 15 in the breast 6 whose different shadings are designed to illustrate different amounts of fatty and glandular tissue along the propagation direction of the x-radiation 3. In the context of projection radiography, this constitutes a simplification which does not result in major deviations from the actual scatter distribution. On the basis of this model assumption, the scatter correction can now be performed, the sequence of which is shown in Figure 3. After image capture 16, a projection image 17 is present which reproduces the primary radiation 10 and secondary radiation incident on the x-ray detector 8. The projection image 17 undergoes data reduction 18 in which different breast regions 13, 14 and 15 are each assigned specific tissue distributions. In addition, information relating to the geometrical relationships, in particular the edges of the breast 6, are obtained. With the aid of the information obtained in data reduction 18 concerning the physical constitution of the breast 6, a scatter beam spread function (SBSF) 20 assignable to the particular breast region 13, 14 and 15 can be looked up in a breast SBSF atlas 19. Using the SBSFs 20 and an estimate for the primary radiation 10, a scatter correction 21 can then be performed. The correction values generated as part of scatter correction 21 can be directly applied to the projection images 17 if the scatter correction has been calculated for each of the detector pixels 9 of the x-ray detector 8. Because of the minimal scatter variation across the x-ray detector 8, it may be sufficient to perform the scatter correction for selected detector regions. These can be individual grid points or groups of detector pixels 9. The scatter correction for the detector pixels 9 for which no scatter correction has yet been determined can be determined by an interpolation 22 which produces a correction image 23 having the same resolution as the projection image 17. By combination 24 of the projection image 17 and the correction image 23 there is finally produced a finished structure image 25 which preferably contains exclusively the structure of the breast 6 imaged by the primary radiation 10. The requirements for the radiation correction described here and the associated processing steps will now be described in detail below: Requirements: It is firstly assumed that the sensitivity spectrum N(E) critical for imaging is known: the radiation of the x-ray tubes is polychromatic, the energy spectrum Qu(E) of the photons emitted as bremsstrahlung at the anode depends on the high voltage U applied with which the electrons are accelerated from the cathode to the anode; the maximum photon energy is then Emax(U) = U(keV/kV)=eU; however, it is not just the image spectrum that is critical for imaging, but also the transparency of spectral filters W(E) used and the spectral response sensitivity ΗD(E) of the detector 8. The resulting (normalized) spectral distribution is defined by: It is secondly assumed that - for a given resulting spectral distribution NU(E) and given breast layer thickness H which is defined by the spacing of the compression plates 5 - the attenuation of the detector signal (of primary x-radiation, without scatter) is present in pre-calculated form as a function of the proportion of glandular tissue or fatty tissue(if necessary validated by measurements), i.e. the following function is given in tabular form: it being assumed that the compressed breast 6 completely fills out the layer thickness H between the compression plates 5. As shown in Figure 4 this condition is no longer met in the region of a few cm near a breast tip 26 and outside in the region of unattenuated radiation. As will be explained in detail below, these image field regions must be dealt with separately as part of a pre-correction, e.g. by suitable extrapolation of the tissue layer thickness H to 0. For mathematical reasons the logarithmic attenuation signal is more useful than the non-logarithmic attenuation function F in equation (# 2): The function fH is monotonic and continuous and consequently invertible, e.g. by inverse interpolation. It can therefore be assumed that the inverse function is also available in tabular form. It is thirdly assumed that the so-called breast SBSF atlas 19 is available, for the method described here is based on knowledge of the relevant SBSFs 20 (scatter beam spread functions). An SBSF 20 describes in each case the spatial intensity distribution of the scatter on the x-ray detector 8 implemented as a flat-panel detector for a thin x-ray beam which penetrates the scatter object (breast) according to Figure 1 at a predefined location. The SBSF 20 depends on capture parameters and on object parameters. Capture parameters are, for example, the tube voltage which affects the photon emission spectrum which, moreover, is also dependent on the anode material, the pre-filtering, the air gap, the SID (source-image distance), the collimation, the spectral response sensitivity of the x-ray detector 8 and the presence or absence of an anti-scatter grid. An object parameter is on the one hand the layer thickness H of the breast 6 and, on the other, the different proportion of glandular or fatty tissue along the propagation direction of the x-radiation 3. It is assumed that the SBSFs 20 are available for the most important capture and object parameters arising, i.e. that there exists a set of tables created in advance, the so- called breast SBSF atlas 19, which can be used to determine with sufficient accuracy the associated SBSF 20 for the specific capture conditions for each proportion of fatty and glandular tissue (scatter material distribution) along an x- ray beam, e.g. by interpolation in the breast SBSF atlas 19 or by semi-empirical conversions for parameters on which the SBSF is only weakly dependent or for which functional dependencies are known, such as in the case of the SID. The breast SBSF atlas 19 is created in advance by means of Monte Carlo simulation calculation. Monte Carlo simulation permits adequate modeling of the physical processes of absorption and multiple scattering (predominantly coherent scattering in the lower frequency range in mammography) during passage through the scatter object, in particular the breast 6, taking account of the capture conditions (anode material, filter, voltage, air gap, SID, field size, and possibly anti-scatter grid). This is the major advantage of the Monte Carlo methods over analytical simulation models which are generally limited to single scattering and in which in most cases various simplifications and approximation are also introduced in order to reduce the cost/complexity. The calculation of scatter distributions on the basis of a Monte Carlo simulation will be familiar to the average person skilled in the art and as such is not part of the subject matter of the application. Description of the individual steps: The scatter correction is subdivided into the following individual steps which can be repeated in an iterative cycle: 0. Empty image calibration and determination of the effective attenuation signal (even a simple general scatter pre-correction being recommended); 1. Determination of the proportion of glandular and fatty tissue; 2. Estimation of the scatter distribution (more accurate SBSF model); 3. Estimation of the primary radiation distribution (scatter correction); 4. Iterative repetition from step 1. or end. Steps 0. and 1. must be performed for each measuring beam, i.e. for each pixel (j,k), the term pixel being used in the following both for the detector pixels 9 and for detector regions comprising a plurality of detector pixels. Step 0: Io calibration and attenuation signal with pre- correction If Io(j,k) is the empty image which is identical to the measured intensity distribution in the beam path without scatter object, I(j,k) the measured intensity distribution with scatter object (breast), then the effective attenuation signal for total radiation, i.e. the superposition of primary and secondary (=scattered) radiation, is given by: In general it will be advisable in respect of step 1. to carry out even here a pre-correction of the scatter background which shall be denoted by S(0) . Methods for estimating S(0) are appended below. S(0) can be location- dependent, but is constant in the simplest case. The pre- correction already provides an estimate of the primary attenuation signal (normalized primary intensity) Step 1: Estimation of specific tissue proportions If P(j, k) is initially assumed to represent only primary radiation without scatter, with equation (#4) and (#3) this yields for the glandular tissue component: and for the glandular tissue weight per unit area [g/cm2] : and for the fatty tissue weight per unit area: As the abovementioned assumption does not strictly apply, an iterative procedure is required. This will be explained in greater detail in connection with remarks concerning step 4. Step 2: Optimally correct estimation of the scatter distribution over the entire projection image This step involves several sub-steps: 2.1 Look-up in the breast SBSF atlas Generation of the SBSF atlas 19 will now be described in further detail. In step 1, a(j, k) was calculated for each beam to which a pixel (j,k) is assigned. For the calculated value of α(j,k) and H as well as further parameters such as air gap, spectrum and other parameters, the associated SBSF 20 is then generally determined from the breast SBSF atlas 19 by interpolation: SBSF((λx, λv); a; H; air gap, voltage, filter, detector,...) SBSF is a two-dimensional function or rather a two- dimensional field (data array) depending on the row and column coordinates on the x-ray detector 8. Each SBSF 20 is focused on a center, namely the particular beam or rather the relevant pixel with the coordinates (0,0) and reduces as a function of distance from the beam center. The distance from the center in both coordinate directions is characterized by an index pair (λX,λy) . The SBSF 20 is a kind of point or line image function, the beam corresponding to the point or line in reality. To characterize the interpolation we employ the notation: This applies to any pixels (j,k) and therefore the total scatter distribution is described by equation (#9). 2.3 Low-pass filtering Because of the multiple scatter processes producing it in the body, the scatter distribution is relatively smooth and therefore exhibits a low-frequency Fourier spectrum. In order to eliminate high-frequency error components induced by the preceding steps, 2-dimensional smoothing is advisable. Step 3: Scatter correction Initially the available data is actually uncorrected, i.e. measurement-based data containing the superposition of primary radiation 10 (direct, unscattered radiation) and secondary radiation 11 (= scatter). After normalization according to equation (#5a) we have: P initially unknown but wanted (normalized) primary- radiation 10 S unknown secondary radiation 11, but estimated (normalized) using the proposed model. Normalization should be understood as division by the intensity distribution without scatter object. Equation (#9) directly yields a subtractive scatter correction: for estimating the primary radiation distribution. Another correction which is recommended in cases of a relatively large amount of secondary radiation 11 is multiplicative scatter correction: Note that the corrections in equation (#11) and equation (#12) are only approximate and do not provide identical results. For S/T«l, however, (#11) becomes (#12). Step 4: Iteration In equation (#11) and (#12) the scatter radiation term S, which for its part must be calculated by equation (#9), appears on the right-hand side ; however, equation (#9) is defined by means of the (unknown) primary radiation P which for its part appears on the left-hand side of equation (#11) and (#12) and is only to be calculated by one of these equations. P therefore appears both on the left- and right- hand side of equation (#11) and (#12). Such implicit equations must be solved iteratively. We write for S in equation S_(#9) : Iteration is performed for the subtractive method as follows: Start of iteration with pre-correction which will be described in greater detail below: For the multiplicative method, the iteration is performed as follows: Start of iteration with pre-correction which will be described in greater detail below: The sequence of iterations is aborted if the result between step n and n+1 only varies slightly. In many cases even one cycle suffices (n=l). SNR improvement by statistical estimation: ML and Bayesian methods Interestingly the multiplicative correction method (#15b) can be derived from a statistical estimation approach according to the maximum likelihood principle (ML). Although in the relevant technical literature a simple convolutional model is used for the scatter operator S_(P) in equation (#13a), for example, in A. H. BAYDUSH, C. E. FLOYD: Improved image quality in digital mammography with image processing. In: Med. Phys., Vol. 27, July 2000, pages 1503 to 1508, ML can basically be applied independently of the specific scatter model, particularly also in the case of the scatter model described here. The feature of a method based on the ML principle is that although the SNR (= signal-to-noise ratio) is usually improved after a few iterations, if the iterations are continued, the noise increases uncontrollably and the SNR deteriorates again. In order to counteract this runaway of the ML algorithm, Bayesian estimation methods are recommended, resulting in algorithms which differ from equation (#15b) in having a stabilizing additional term on the right-hand side. The effect of this additional term on convergence rate, SNR and the compromise between noise and local resolution can be controlled by parameters. Pre-corrections In the previous comments concerning steps 1 and 2.1, equations (#6) and (#7), it was assumed that the compressed breast 6 completely fills out the layer thickness H between the compression plates 5 and that the function fH-1 can be evaluated. As shown in Figure 4 this condition is no longer fulfilled in the region of a few cm near a breast tip 26 and outside in the region of unattenuated x-radiation 3. These image field regions must be dealt with separately as part of a pre-correction. In the region of unattenuated x-radiation 3 outside the breast 6, the effective attenuation signal according to equation (#5a) must theoretically be = 1, but is generally > 1 because of the presence of scatter. The difference in the image region outside the breast 6. From the normal image region of the fully compressed breast 6 to the region near the breast tip 26 a suitable extrapolation of the tissue layer thickness from H to 0 must be performed. In this image region, H must therefore generally be assumed to be variable in equations (#2), (#6) and (#7) . If necessary, segmentation into 3 image regions can also be performed as described in K. NYKANEN, S. SILTANEN: X-ray scattering in full field digital mammography. In Med. Phys. , Vol. 30(7), July 2003, pages 1864 to 1873. In the normal image region with constant tissue layer thickness H a scatter pre-correction can look like this: as there has not yet been any evaluation of the tissue proportions (glandular/fatty tissue), initially 100 % fat can be assumed. Although because of the lower density of fat (0.92 compared to 0.97 g/cm3 for glandular tissue) the scatter is underestimated, for a Oth-order correction this estimation is significantly better than no correction at all. α=0 is inserted in equation (#7) and the subsequent equations, making the scatter kernel SBSF location- independent, in particular independent of the pixel index (j.k), and equation (#9) is reduced to a genuine convolution. Creating the breast SBSF atlas Of interest in the SBSF concept is the distribution of the scatter produced in the scatter body in the detector plane when, as shown in Figure 4, the (unscattered) primary radiation (i.e. a mini cone beam 27) is focused on a single detector pixel 9. If this is done consecutively for each detector pixel 9 and all the associated SBSFs 20 are summed, the total scatter distribution is obtained for the case where the entire detector surface is illuminated - and not only individual detector pixels 9. As already described above in connection with the third requirement and step 2, the breast SBSF atlas 19 of the scatter beam spread functions (SBSF) comprises the scatter intensity distributions normalized to the intensity of the primary radiation 10 in the detector pixel 9 (assuming that the mini cone beam 27 is focused on just one pixel 9) as a function of a plurality of different parameter configurations: SBSF(( λx,λy); a;H; air gap, voltage, filter, detector,... ) (#17) and also contains the dependency of the x-ray energy spectrum on the tube voltage, pre-filtering, radiation-sensitive detector material, e.g. the type of scintillation crystal, and the dependency on the presence or absence of an anti- scatter grid and where applicable the dependence on the type of anti-scatter grid as well as the dependence on other parameters. The creation of an SBSF series will now be explained: First the parameters characterizing the relevant mammography machine 1 are defined: SID, air gap, anode material of the x- ray tubes (and associated emission spectra), detector material, pre-filter material (e.g. compression plates), and other parameters. Then comes the compression thickness H, the voltage, the spectral filters used and other variables, the voltage and if necessary the spectral filter (thickness) generally being modified as a function of the compression thickness H in order to optimize image quality. For this parameter configuration, the parameter a describing the tissue composition according to equation (#2a) is varied between 0 (fat only) and 1 (glandular tissue only): the calculation using the tried and tested Monte Carlo method produces a set of different SBSFs 20, each α-value being assigned an SBSF 20. The tissue thickness H is varied between > 0 and up to approximately 10 cm and another set of SBSFs 20 is again calculated for each H . The voltage and the spectral filters can also be varied, the variation being linked to H or also independent of H . However, in the latter case multiple variations are possible. In addition, the calculation can be continued for all the parameter combinations. For calculating the SBSFs 20, simplifications can be performed which are well justified: • Disregarding the beam divergence of the x-radiation 3 due to the cone beam geometry by assuming approximately parallel beam geometry; this is justified in that generally SID >> H ; this is achieved in that the SBSF 20 remains location- and pixel-independent for an identical beam configuration; by identical configuration is meant that, for each pixel, the material distribution is the same along the mini cone beam 27 and in the lateral neighborhood. • To improve the statistics for the Monte Carlo method and reduce the computational complexity, pixels approximately an order of magnitude larger (e.g. 1 xl mm2 or 2 x 2 mm2) than the actual detector pixels 9 ( calculate the SBSFs 20; this is justified by the low- frequency Fourier spectrum of the spatial scatter distribution. • The succession of fatty and glandular tissue is replaced by a mixture; although the scatter depends (for the same weight per unit area and path length) on whether the denser tissue is nearer the x-ray detector 8 or nearer the radiation source 2, according to J.M. DINTEN and J.M. Voile: Physical model based restoration of mammographies. In Proc. SPIE, Vol. 3336, 1998, 641-650, the differences occurring under mammographic conditions can be disregarded. Advantages The solution proposed here has the following advantages: If required, the method can be incorporated in existing mammography machines without mechanical reconstruction. Moreover it is a method which on the one hand shares the adequacy of physical modeling using the Monte Carlo method, but on the other hand - because all the time-consuming calculations are carried out in advance where possible and the necessary data is stored in tables - ultimately involves relatively low computational complexity for the scatter correction. The modeling accuracy of the scatter correction described here is essentially greater than that of the known (analytical) physical models, as a number of simplifying assumptions and approximations can be dispensed with. The possibilities of the scatter correction proposed here go far beyond the possibilities of the long known convolution/deconvolution methods. Disregarding the specific technical embodiment of the method and looking at it from a mathematical standpoint, the method can be regarded in the mathematical sense as a generalization of the long known convolution/deconvolution method. On the one hand, by using approximations and dispensing with accuracy, it can be categorized as method of this type and then shares its advantages, such as the possibility of using the so-called FFT (= Fast Fourier Transformation). On the other hand, however, the method described here can also be extended in terms of SNR improvement, e.g. by extending the iterative multiplicative algorithm in the direction of statistical Bayesian estimation. In this context it should be re-emphasized that only pre- calculation of the SBSFs 20 enables the method described here to be performed in full generality. Examples Example 1: In this example, scatter correction is performed, as described above with equations (#5)-(#9) and (#13)-(#15), using homogeneous location-dependent scatter beam spread functions 20 (= SBSF). For creating said scatter beam spread functions 20 it is assumed by way of simplification that the tissue distribution characterized by the proportion α(j, k) of glandular tissue along the beam leading from the source to the detector pixel continues in a constantly homogeneous manner according to Figure 4 at right angles to the beam, i.e. parallel to the compression plates 5. It is therefore assumed with respect to the scatter contribution of the beam in the pixel (j, k) that the tissue composition in the lateral neighborhood of the beam does not vary abruptly. Although this is no longer relevant at the edge of the breast, special treatment could be provided there. Note, however, that the actual location-dependent inhomogeneity of the tissue composition is allowed for by a specifically different amount of glandular tissue α(j',k') for each pixel (j', k') and a specific scatter contribution dependent thereon. The SBSFs 20 are therefore generally different for each pixel. Example 1a: In this example la the method is essentially performed as in example 1. However, the following simplifications are made: For each pre-specified layer thickness and the other parameters such as voltage and pre-filtering, a common SBSF 20 is used for all the pixels. In this case the SBSF 20 is therefore selected on a location-independent basis. The selection can be made, for example, by suitable averaging over the tissue compositions present. ∆S in equation (#7) and (#9) then becomes independent of the pixel index (j,k); the double index (j,k) can - similarly as in equations (#16a) to (#16c) - be omitted. The important feature is that the integral in equation (#9) becomes a genuine convolution which can be efficiently executed by FFT (=Fast Fourier Transformation). Example 1b: In this example lb the method is likewise performed essentially as in example 1. In this case, however, a uniform convolution kernel (for all the layer thicknesses) is used for the scatter calculation. The fact that for a small layer thickness relatively less scatter is produced than for a large layer thickness must be taken into account by means of scaling factors which are a function of the layer thickness and other parameters such as voltage and filtering. Approximately the same computational complexity is necessary for example 1b as for example la. On the other hand, much less memory space is required for storing the breast SBSF atlas 19 in this example. Notes on examples la and lb: In general the simplified examples la and lb share the characteristic that the convolutional models for the scatter can be inverted using the Fourier transformation. This is known as deconvolution. The examples described here differ from the conventional deconvolution methods in using one or more scatter beam spread functions 20 obtained in advance by Monte Carlo simulation. With regard to performing deconvolution, reference is made to a publication by J.A. SEIBERT and J.M. BOONE: X-ray scatter removal by deconvolution. In Med.Phys., Vol. 15, 1988, pages 567 to 575. Reference is also made to the more recent publication P. ABBOTT et al: Image deconvolution as an aid to mammographic artifact identification I: basic techniques. In: Proc.SPIE, Vol. 3661, 1999, pages 698 to 709 which deals with deconvolution using regularization techniques for noise suppression. Another deconvolution method with thickness- dependent convolution is described in D.G. TROTTER et al: Thickness-Dependent Scatter-Correction Algorithm for Digital Mammography. In: Proc. SPIE, Vol. 4682, 2002, pages 469 to 478. In this method an iteration with relaxation is performed. Example 2 In this example the method is essentially performed as in example 1, but employing scatter beam spread functions 20 which have been calculated for an inhomogeneous medium. Figure 5, for example, illustrates the case where a breast region 28 has a different composition from that of a surrounding breast region 29. This enables it to be taken into account that the SBSF 20 depends not only on the tissue composition along the mini cone beam 27 supposedly focused on the detector pixel 27 but also on the tissue composition in the lateral neighborhood into which photons are scattered and can be further scattered again in the direction of the pixel. However, the effective extent of the lateral neighborhood is not very large because of the average free path length mammography energy range between about 20 and 40 keV. It would therefore suffice to assume the tissue composition to be homogeneous in a lateral half space, but generally different from the mini cone beam 27. The allowance for inhomogeneous SBSFs 20 with differences between beam and neighborhood might be relevant particularly at the breast edge . This example therefore constitutes a generalization of the above-described examples 1, la and lb, as in this case the SBSFs 20 depend not only on a tissue parameter a, but also on a surrounding tissue parameter y to be newly introduced. In this case the breast SBSF atlas 19 would therefore have an additional dimension. For the sake of clarity, the following table compares the different characteristics of examples 1, la, lb and 2: The method described here can also be applied to so-called dual energy methods which will be known to the average person skilled in the art. With the so-called dual energy method, which is used primarily in mammography or in bone densitometry, images are recorded simultaneously using two different energy spectra. The recordings using different energy spectra are provided by two different voltages and if possible also different spectral filtering so that the spectral regions effectively corresponding to the two measurements overlap one another as little as possible. By means of a computational process which is essentially based on the solution of a generally nonlinear system of two equations assigned to the two spectra, finer tissue differentiation can be achieved compared to a recording using one energy spectrum. For computation to be successful, the scatter components must be eliminated as much as possible, as otherwise the artifacts induced by the scatter components are in some circumstances stronger than the actual tissue image. Because of the differences in scatter for the two spectra, effective scatter correction is therefore critically important for the quality of the dual energy method. The proposed scatter correction method can also be used in this context. The geometrical parameters are identical for the two recordings, but the spectrally dependent parameters are different. The correction must be carried out for each of the two recordings according to the described formula, the only difference being that different SBSFs 20 must be used according to the different spectra. WE CLAIM 1. An apparatus for projection radiography, in particular for mammography, comprising a radiation source (2) which emits radiation (3), a detector (8) and, downstream of said detector (8), a processing unit (12) which determines in an approximately manner, on the basis of the projection data (17) supplied by the detector (8), the scatter material distribution of the object under examination (6) and reads out scatter information (20) from a data memory (19) as a function of the scatter material distribution and corrects the projection data (17) in respect of the scatter content (11) on the basis of the scatter information (20), characterized in that - the scatter information (20) is determined by Monte Carlo simulations which calculate the interactions of the photons with a scatter material distribution in respect of different scatter material distributions, in that - the scatter information is scatter distributions (20) which describes a scatter-induced distribution over adjacent image regions of the radiation (14) emitted by the radiation source (2) and directed to a particular image region, and in that - the processing unit (12) determines a scatter distribution (11) in an image region of the projection image (17) by calculating and adding, for each image region, the scatter contributions (11) of the surrounding image regions. 2. The apparatus as claimed in claim 1, wherein the scatter distributions (20) are scalable using the intensity of the unscattered primary radiation (10) incident on the detector (8). 3. The apparatus as claimed in claim 1 or 2, wherein the processing unit (12) evaluates scatter information (20) specific to the particular scatter material distribution for different image regions of a projection image (17). 4. The apparatus as claimed in one of claims 1 to 3, wherein the processing unit (12) determines a scatter distribution (11) in an image region of the projection image (17) by convolving the primary radiation distribution with a scatter distribution (20). 5. The apparatus as claimed in one of claims 1 to 4, wherein the processing unit (12) determines the unscattered primary radiation (10) by solving the implicit equation P + S (P) = T, where P is the distribution of the unscattered primary radiation (10), S (P) the secondary radiation distribution (11) dependent on the unscattered primary radiation (10) and T the measured total radiation distribution in the projection images (17). 6. The apparatus as claimed in one of claims 1 to 5, wherein processing unit (12), for a first approximate determination of the scatter material distribution of the object under examination (6), estimates the amount of scatter (11) contained in the projection image (17) on the basis of scatter information associated with a typical scatter material distribution. 7. The apparatus as claimed in one of claims 1 to 6, wherein the processing unit (12) performs the processing steps iteratively as claimed in one of claims 1 to 6. 8. The apparatus as claimed in one of claims 1 to 7, wherein the scatter information (20) is calculated under the assumption of a homogeneous scatter material distribution in the beam direction. 9. The apparatus as claimed in one of claims 1 to 8, wherein scatter information (20) allowing for the outer contour (26) of the object under examination (6) is stored in the data memory (19). 10.The apparatus as claimed in one of claims 1 to 9, wherein scatter information (20) calculated under the assumption of a homogeneous scatter material distribution perpendicular to the beam direction is stored in the data memory (19). 11.The apparatus as claimed in one of claims 1 to 10, wherein scatter information (20) determined allowing for an inhomogeneous scatter material distribution perpendicular to the beam direction is stored in the data memory (19). 12.The apparatus as claimed in one of claims 1 to 11, wherein scatter information depending on parameters of the radiation source (2) is stored in the data memory (19). 13.The apparatus as claimed in one of claims 1 to 12, wherein the processing unit (12) determines the scatter components (12) at selected grid points and works out the correction values for individual detector elements (9) by interpolation between the grid points. 14. The apparatus as claimed in one of claims 1 to 13, wherein the object under examination (6) is compressible in a compression device (5) and that the processing unit (12) uses the physical configuration of the surfaces of the compression device (5) facing the object under examination to determine the path length of the radiation (3) through the object under examination. ABSTRACT TITLE: "AN APPARATUS AND A METHOD FOR SCATTER CORRECTION IN PROJECTION RADIOGRAPHY" The invention relates to an apparatus for projection radiography, in particular for mammography, comprising a radiation source (2) which emits radiation (3), a detector (8) and, downstream of said detector (8), a processing unit (12) which determines in an approximately manner, on the basis of the projection data (17) supplied by the detector (8), the scatter material distribution of the object under examination (6) and reads out scatter information (20) from a data memory (19) as a function of the scatter material distribution and corrects the projection data (17) in respect of the scatter content (11) on the basis of the scatter information (20), the scatter information (20) is determined by Monte Carlo simulations which calculate the interactions of the photons with a scatter material distribution in respect of different scatter material distributions, in that the scatter information is scatter distributions (20) which describes a scatter-induced distribution over adjacent image regions of the radiation (14) emitted by the radiation source (2) and directed to a particular image region, and in that the processing unit (12) determines a scatter distribution (11) in an image region of the projection image (17) by calculating and adding, for each image region, the scatter contributions (11) of the surrounding image regions.

Full Text

Description
Apparatus and methods for scatter correction in projection
radiography, particularly mammography
The invention relates to an apparatus for projection
radiography comprising a radiation source, a detector and,
downstream of said detector, a processing unit which uses the
projection data supplied by the detector to approximately
determine the scatter material distribution of the object
under examination and which reads out scatter information from
a data memory as a function of the scatter material
distribution and corrects the projection data in respect of
the scatter component on the basis of said scatter
information.
The invention further relates to a method with scatter
correction for projection radiography and a method for
obtaining scatter information.
An apparatus and methods of this kind are known from US
6,104,777 A.
The scatter produced in the object under examination (breast),
the intensity of which may almost attain the order of
magnitude of the direct, unscattered, image-producing primary
radiation, results in image quality impairment by reducing the
contrast, increasing the noise, and ultimately in respect of
the image post-processing methods used for differentiating
between various types of tissue in the images produced, in
particular between glandular and fatty tissue in the breast.
For differentiation in respect of two tissue types, techniques
using a single energy spectrum, i.e. a single x-ray tube

voltage, or the dual energy method using two voltage values
are known in mammography. In both cases scatter compensation
is required; with the dual energy method this is also because
the amount of scatter is different in the two energy spectra.
To reduce scatter, mechanical measures have already been
proposed. The use of slit collimators requires mechanical
displacement of the slit collimators over the breast
measuring field and is therefore time-consuming. Anti-scatter
grids not only reduce the scatter, but also the imaging-
producing primary radiation. Arguments in favor of dispensing
with anti-scatter grids have been ongoing for years. For
compression thicknesses of less than 4-5 cm, the dose could
even be reduced or the SNR (= signal-to-noise ratio)
increased if the grid is removed. On the other hand there are
applications in which it is not technically possible to use a
grid, e.g. in tomosynthesis.
A large number of computer correction methods have already
been proposed. Methods of this kind are known e.g. from M.
DARBOUX, J.M. DINTEN: Physical model based scatter correction
in mammography. In: Proc. SPIE, Vol. 3032, 1997, pages 405 to
410 and from J.M. DINTEN and J.M. VOLLE: Physical model based
restoration of mammographies. In. Proc. SPIE, Vol. 3336,
1998, pages 641 to 650 and in US 6,104,777 A. These are
convolution/deconvolution methods in which a scatter
intensity distribution is approximated as a convolution of
the primary radiation distribution using suitable convolution
kernels . In the cited documents an analytical model is thus
proposed with which the physical scatter process in the
scatter object (breast) is explicitly computed as an integral
transformation. However, this explicit analytical
representation only describes first-order scatter, not

multiple scatter. The intensity distribution of multiply
scattered photons is assumed to be a spatially constant
background over the detector surface and must be estimated
from previously determined tables. The analytical model for
calculating just the first-order scatter contribution
requires 4-dimensional numerical integrations (3 space
coordinates + energy spectrum) for each detector pixel, i.e.
it is compute-intensive. Approximations are therefore
required in order to reduce the computational complexity.
Because of the high computational complexity it is proposed
to perform the calculations in advance and tabulate the
results.
In addition, W. KALENDER: Monte Carlo calculations of x-ray
scatter data for diagnostic radiology. In: Phys. Med. Biol.,
1981, Vol. 26, No. 5, pages 835 to 849 describes the use of
Monte Carlo methods for simulating radiation propagation in
radiography.
Proceeding from this prior art, the object of the invention
is therefore to specify an apparatus and methods enabling
improved scatter correction compared to the prior art to be
performed.
This object is achieved by an apparatus and the methods
having the features set forth in the independent claims.
Advantageous further developments and embodiments are set
forth in claims dependent thereon.
With the apparatus and method, the projection images supplied
by a detector are analyzed in a processing unit. It is first
attempted to approximately determine the scatter material
distribution, in mammography typically the proportions of

glandular and fatty tissue, of an object under examination.
In another processing step, scatter information depending on
the scatter material distribution is read out of a data
memory. This scatter information can then be used to correct
the projection images in respect of the scatter content of
said projection images, it being essential that the scatter
information read out of the data memory has been determined
in advance by a Monte Carlo simulation which takes the
multiple interaction of the photons with the object under
examination into account.
The basis for the solution described here is optimally
correct physical modeling. In contrast to the prior art,
modeling is possible which takes a much larger number of
details into account, namely in the following respect: the
occurrence of multiple scatter and polychromasia and the
geometrical relationships, in particular the peculiarities of
the scatter distribution at the edges of the objects, can be
simulated. During the scatter correction itself, merely table
access is required, possibly with subsequent interpolation,
and the calculation of the scatter distribution in the
detector plane is reduced to 2-dimensional integrations over
the detector plane. Despite the relatively simple
implementation of the scatter correction, the procedure is
not limited to special cases and requires no drastic
simplifications or approximations, such as simplified
acquisition geometry, monochromatized radiation,
simplifications of the physical model or a Taylor development
by orders of approximation or similar.
In a preferred embodiment, the scatter material distribution
is specifically determined for different regions of the
projection image. To perform scatter correction in an image

region, the scatter contributions of the surrounding
contributions which depend on the specific scatter material
distribution are then determined and corrected accordingly.
In this way it is possible to take local scatter variations
into account.
In another preferred embodiment, in the region of the edges
of the object under examination, scatter information which
takes the particular geometrical relationships in the region
of the object edge into account is used for scatter
correction.
The scatter information is preferably obtained under the
assumption that the scatter material distribution is
homogeneous along the radiation direction. Particularly in
the context of mammography, such an assumption only results
in slight deviations from the actual scatter distribution.
The specific scatter information associated with an image
area can also be calculated under the assumption that the
object under examination is also homogeneously structured in
the direction perpendicular to the beam, thereby simplifying
the calculation of the scatter information.
However, if particularly high accuracy in calculating the
scatter information is required, any inhomogeneity
perpendicular to the beam direction can be taken into
account.
In a preferred embodiment, the scatter material distribution
is determined by analyzing the ratio of incident radiation
intensity to the unscattered primary radiation in an image
region, the values for the primary radiation being

ascertained by means of a scatter correction based on the
scatter information associated with a characteristic
homogeneous scatter material distribution.
The processing steps carried out by the processing unit can
also be executed iteratively, the calculated primary
radiation components being used to simplify the approximate
calculation of the scatter components and thus arrive at
improved values for the primary radiation.
The scatter correction does not generally need to be
performed at full detector resolution. It may occasionally be
sufficient to perform scatter correction at selected grid
points and interpolate between the determined scatter
correction values at the selected grid points.
Further advantages and embodiments of the invention will
emerge from the following description in which examples of
the invention are specifically explained with reference to
the accompanying drawings in which:
Figure 1 shows the configuration of a mammography machine
in which a breast under examination is compressed
between two compression plates and irradiated
with x-radiation;
Figure 2 illustrates a simplified breast structure assumed
for calculating the scatter correction;
Figure 3 shows a flowchart of a method performed for
scatter correction;

Figure 4 shows the breast tissue distribution assumed for
calculating a simple scatter beam spread
function; and
Figure 5 shows the breast structure assumed for
calculating a precise scatter beam spread
function.
Figure 1 shows the configuration of a mammography machine 1
in which x-radiation 3 is produced with the aid of a
radiation source 2. The divergence of the x-radiation 3 is
limited if necessary using a collimator 4 which is shown in
Figure 1 as a single beam diaphragm. However, the collimator
4 can also be contrived such that a plurality of virtually
parallel x-ray beams is produced. Such a collimator 4 can be
implemented e.g. as an iris diaphragm.
The mammography machine 1 additionally has compression plates
5 between which a breast 6 is compressed. The x-radiation 3
passes through the compression plates 5 and the breast 6 and
generally crosses an air gap 7 before the x-radiation 3 is
incident on an x-ray detector 8 comprising a plurality of
individual detector elements 9, the so-called detector
pixels.
The portion of x-radiation 3 passing through the breast 6
without interacting with the breast 6 is known as the primary
radiation 10. On the other hand, the portions of x-radiation
3 incident on the x-ray detector 8 after at least one
scattering within the breast 6 are referred to as secondary
radiation 11.

It should be noted that the term scatter is to be understood
as any kind of interaction between the x-radiation 3 and the
material of the breast 6 causing a change in the propagation
direction of the photons of the x-radiation 3.
Since, as mentioned in the introduction, the secondary
radiation 11 may considerably distort the structure of the
breast 6 imaged by the primary radiation 10, it is
advantageous if the secondary radiation 11 can be removed
from the projection images of the breast 6 captured by the x-
ray detector 8. For this purpose, a processing unit 12
connected downstream of the x-ray detector 8 performs a
scatter correction. In order to be able to perform the
scatter correction, model assumptions are made concerning the
structure of the breast 6. In particular it is assumed that
the tissue structure of the breast 6 which is essentially
composed of glandular and fatty tissue can be described by a
homogeneous tissue distribution along the propagation
direction of the x-radiation 3. Figure 2 accordingly shows
different regions 13, 14 and 15 in the breast 6 whose
different shadings are designed to illustrate different
amounts of fatty and glandular tissue along the propagation
direction of the x-radiation 3. In the context of projection
radiography, this constitutes a simplification which does not
result in major deviations from the actual scatter
distribution.
On the basis of this model assumption, the scatter correction
can now be performed, the sequence of which is shown in
Figure 3.
After image capture 16, a projection image 17 is present
which reproduces the primary radiation 10 and secondary

radiation incident on the x-ray detector 8. The projection
image 17 undergoes data reduction 18 in which different
breast regions 13, 14 and 15 are each assigned specific
tissue distributions. In addition, information relating to
the geometrical relationships, in particular the edges of the
breast 6, are obtained. With the aid of the information
obtained in data reduction 18 concerning the physical
constitution of the breast 6, a scatter beam spread function
(SBSF) 20 assignable to the particular breast region 13, 14
and 15 can be looked up in a breast SBSF atlas 19. Using the
SBSFs 20 and an estimate for the primary radiation 10, a
scatter correction 21 can then be performed. The correction
values generated as part of scatter correction 21 can be
directly applied to the projection images 17 if the scatter
correction has been calculated for each of the detector
pixels 9 of the x-ray detector 8. Because of the minimal
scatter variation across the x-ray detector 8, it may be
sufficient to perform the scatter correction for selected
detector regions. These can be individual grid points or
groups of detector pixels 9. The scatter correction for the
detector pixels 9 for which no scatter correction has yet
been determined can be determined by an interpolation 22
which produces a correction image 23 having the same
resolution as the projection image 17. By combination 24 of
the projection image 17 and the correction image 23 there is
finally produced a finished structure image 25 which
preferably contains exclusively the structure of the breast 6
imaged by the primary radiation 10.
The requirements for the radiation correction described here
and the associated processing steps will now be described in
detail below:

Requirements:
It is firstly assumed that the sensitivity spectrum N(E)
critical for imaging is known:
the radiation of the x-ray tubes is polychromatic, the energy
spectrum Qu(E) of the photons emitted as bremsstrahlung at the
anode depends on the high voltage U applied with which the
electrons are accelerated from the cathode to the anode; the
maximum photon energy is then Emax(U) = U(keV/kV)=eU; however, it
is not just the image spectrum that is critical for imaging,
but also the transparency of spectral filters W(E) used and
the spectral response sensitivity ΗD(E) of the detector 8. The
resulting (normalized) spectral distribution is defined by:

It is secondly assumed that - for a given resulting spectral
distribution NU(E) and given breast layer thickness H which
is defined by the spacing of the compression plates 5 - the
attenuation of the detector signal (of primary x-radiation,
without scatter) is present in pre-calculated form as a
function of the proportion of glandular tissue or fatty

tissue(if necessary validated by measurements), i.e. the
following function is given in tabular form:

it being assumed that the compressed breast 6 completely
fills out the layer thickness H between the compression
plates 5. As shown in Figure 4 this condition is no longer
met in the region of a few cm near a breast tip 26 and
outside in the region of unattenuated radiation. As will be
explained in detail below, these image field regions must be
dealt with separately as part of a pre-correction, e.g. by
suitable extrapolation of the tissue layer thickness H to 0.

For mathematical reasons the logarithmic attenuation signal
is more useful than the non-logarithmic attenuation function
F in equation (# 2):

The function fH is monotonic and continuous and consequently
invertible, e.g. by inverse interpolation. It can therefore
be assumed that the inverse function

is also available in tabular form.
It is thirdly assumed that the so-called breast SBSF atlas 19
is available, for the method described here is based on
knowledge of the relevant SBSFs 20 (scatter beam spread
functions). An SBSF 20 describes in each case the spatial
intensity distribution of the scatter on the x-ray detector 8
implemented as a flat-panel detector for a thin x-ray beam
which penetrates the scatter object (breast) according to
Figure 1 at a predefined location. The SBSF 20 depends on
capture parameters and on object parameters.
Capture parameters are, for example, the tube voltage which
affects the photon emission spectrum which, moreover, is also
dependent on the anode material, the pre-filtering, the air
gap, the SID (source-image distance), the collimation, the
spectral response sensitivity of the x-ray detector 8 and the
presence or absence of an anti-scatter grid.

An object parameter is on the one hand the layer thickness H
of the breast 6 and, on the other, the different proportion
of glandular or fatty tissue along the propagation direction
of the x-radiation 3.
It is assumed that the SBSFs 20 are available for the most
important capture and object parameters arising, i.e. that
there exists a set of tables created in advance, the so-
called breast SBSF atlas 19, which can be used to determine
with sufficient accuracy the associated SBSF 20 for the
specific capture conditions for each proportion of fatty and
glandular tissue (scatter material distribution) along an x-
ray beam, e.g. by interpolation in the breast SBSF atlas 19
or by semi-empirical conversions for parameters on which the
SBSF is only weakly dependent or for which functional
dependencies are known, such as in the case of the SID.
The breast SBSF atlas 19 is created in advance by means of
Monte Carlo simulation calculation. Monte Carlo simulation
permits adequate modeling of the physical processes of
absorption and multiple scattering (predominantly coherent
scattering in the lower frequency range in mammography)
during passage through the scatter object, in particular the
breast 6, taking account of the capture conditions (anode
material, filter, voltage, air gap, SID, field size, and
possibly anti-scatter grid). This is the major advantage of
the Monte Carlo methods over analytical simulation models
which are generally limited to single scattering and in which
in most cases various simplifications and approximation are
also introduced in order to reduce the cost/complexity. The
calculation of scatter distributions on the basis of a Monte
Carlo simulation will be familiar to the average person

skilled in the art and as such is not part of the subject
matter of the application.
Description of the individual steps:
The scatter correction is subdivided into the following
individual steps which can be repeated in an iterative cycle:
0. Empty image calibration and determination of the
effective attenuation signal (even a simple general
scatter pre-correction being recommended);
1. Determination of the proportion of glandular and fatty
tissue;
2. Estimation of the scatter distribution (more accurate
SBSF model);
3. Estimation of the primary radiation distribution
(scatter correction);
4. Iterative repetition from step 1. or end.
Steps 0. and 1. must be performed for each measuring beam,
i.e. for each pixel (j,k), the term pixel being used in the
following both for the detector pixels 9 and for detector
regions comprising a plurality of detector pixels.
Step 0: Io calibration and attenuation signal with pre-
correction
If Io(j,k) is the empty image which is identical to the
measured intensity distribution in the beam path without
scatter object, I(j,k) the measured intensity distribution
with scatter object (breast), then the effective attenuation
signal for total radiation, i.e. the superposition of primary
and secondary (=scattered) radiation, is given by:

In general it will be advisable in respect of step 1. to
carry out even here a pre-correction of the scatter
background which shall be denoted by S(0) . Methods for
estimating S(0) are appended below. S(0) can be location-
dependent, but is constant in the simplest case. The pre-
correction already provides an estimate of the primary
attenuation signal (normalized primary intensity)

Step 1: Estimation of specific tissue proportions
If P(j, k) is initially assumed to represent only primary
radiation without scatter, with equation (#4) and (#3) this
yields for the glandular tissue component:

and for the glandular tissue weight per unit area [g/cm2] :

and for the fatty tissue weight per unit area:

As the abovementioned assumption does not strictly apply, an
iterative procedure is required. This will be explained in
greater detail in connection with remarks concerning step 4.

Step 2: Optimally correct estimation of the scatter
distribution over the entire projection image
This step involves several sub-steps:
2.1 Look-up in the breast SBSF atlas
Generation of the SBSF atlas 19 will now be described in
further detail.
In step 1, a(j, k) was calculated for each beam to which a
pixel (j,k) is assigned. For the calculated value of α(j,k)
and H as well as further parameters such as air gap,
spectrum and other parameters, the associated SBSF 20 is then
generally determined from the breast SBSF atlas 19 by
interpolation:
SBSF((λx, λv); a; H; air gap, voltage, filter, detector,...)
SBSF is a two-dimensional function or rather a two-
dimensional field (data array) depending on the row and
column coordinates on the x-ray detector 8. Each SBSF 20 is
focused on a center, namely the particular beam or rather the
relevant pixel with the coordinates (0,0) and reduces as a
function of distance from the beam center. The distance from
the center in both coordinate directions is characterized by
an index pair (λX,λy) . The SBSF 20 is a kind of point or line
image function, the beam corresponding to the point or line
in reality.
To characterize the interpolation we employ the notation:

This applies to any pixels (j,k) and therefore the total
scatter distribution is described by equation (#9).
2.3 Low-pass filtering
Because of the multiple scatter processes producing it in the
body, the scatter distribution is relatively smooth and
therefore exhibits a low-frequency Fourier spectrum. In order
to eliminate high-frequency error components induced by the
preceding steps, 2-dimensional smoothing is advisable.
Step 3: Scatter correction
Initially the available data is actually uncorrected, i.e.
measurement-based data containing the superposition of
primary radiation 10 (direct, unscattered radiation) and
secondary radiation 11 (= scatter).
After normalization according to equation (#5a) we have:

P initially unknown but wanted (normalized) primary-
radiation 10
S unknown secondary radiation 11, but estimated
(normalized) using the proposed model.
Normalization should be understood as division by the
intensity distribution without scatter object.
Equation (#9) directly yields a subtractive scatter
correction:

for estimating the primary radiation distribution.
Another correction which is recommended in cases of a
relatively large amount of secondary radiation 11 is
multiplicative scatter correction:

Note that the corrections in equation (#11) and equation
(#12) are only approximate and do not provide identical
results. For S/T«l, however, (#11) becomes (#12).
Step 4: Iteration
In equation (#11) and (#12) the scatter radiation term S,
which for its part must be calculated by equation (#9),
appears on the right-hand side ; however, equation (#9) is
defined by means of the (unknown) primary radiation P which
for its part appears on the left-hand side of equation (#11)
and (#12) and is only to be calculated by one of these

equations. P therefore appears both on the left- and right-
hand side of equation (#11) and (#12). Such implicit
equations must be solved iteratively. We write for S in
equation S_(#9) :

Iteration is performed for the subtractive method as follows:
Start of iteration with pre-correction which will be
described in greater detail below:

For the multiplicative method, the iteration is performed as
follows:
Start of iteration with pre-correction which will be
described in greater detail below:

The sequence of iterations is aborted if the result between
step n and n+1 only varies slightly. In many cases even one
cycle suffices (n=l).
SNR improvement by statistical estimation: ML and Bayesian
methods
Interestingly the multiplicative correction method (#15b)
can be derived from a statistical estimation approach
according to the maximum likelihood principle (ML). Although
in the relevant technical literature a simple convolutional
model is used for the scatter operator S_(P) in equation
(#13a), for example, in A. H. BAYDUSH, C. E. FLOYD: Improved
image quality in digital mammography with image processing.
In: Med. Phys., Vol. 27, July 2000, pages 1503 to 1508, ML
can basically be applied independently of the specific
scatter model, particularly also in the case of the scatter
model described here.
The feature of a method based on the ML principle is that
although the SNR (= signal-to-noise ratio) is usually
improved after a few iterations, if the iterations are
continued, the noise increases uncontrollably and the SNR
deteriorates again. In order to counteract this runaway of
the ML algorithm, Bayesian estimation methods are
recommended, resulting in algorithms which differ from
equation (#15b) in having a stabilizing additional term on
the right-hand side. The effect of this additional term on
convergence rate, SNR and the compromise between noise and
local resolution can be controlled by parameters.
Pre-corrections

In the previous comments concerning steps 1 and 2.1,
equations (#6) and (#7), it was assumed that the compressed
breast 6 completely fills out the layer thickness H between
the compression plates 5 and that the function fH-1 can be
evaluated. As shown in Figure 4 this condition is no longer
fulfilled in the region of a few cm near a breast tip 26 and
outside in the region of unattenuated x-radiation 3. These
image field regions must be dealt with separately as part of
a pre-correction. In the region of unattenuated x-radiation 3
outside the breast 6, the effective attenuation signal
according to equation (#5a) must theoretically be = 1, but is
generally > 1 because of the presence of scatter. The
difference

in the image region outside the breast 6.
From the normal image region of the fully compressed breast
6 to the region near the breast tip 26 a suitable
extrapolation of the tissue layer thickness from H to 0
must be performed. In this image region, H must therefore
generally be assumed to be variable in equations (#2), (#6)
and (#7) .
If necessary, segmentation into 3 image regions can also be
performed as described in K. NYKANEN, S. SILTANEN: X-ray

scattering in full field digital mammography. In Med. Phys. ,
Vol. 30(7), July 2003, pages 1864 to 1873.
In the normal image region with constant tissue layer
thickness H a scatter pre-correction can look like this: as
there has not yet been any evaluation of the tissue
proportions (glandular/fatty tissue), initially 100 % fat
can be assumed. Although because of the lower density of fat
(0.92 compared to 0.97 g/cm3 for glandular tissue) the
scatter is underestimated, for a Oth-order correction this
estimation is significantly better than no correction at
all. α=0 is inserted in equation (#7) and the subsequent
equations, making the scatter kernel SBSF location-
independent, in particular independent of the pixel index
(j.k), and equation (#9) is reduced to a genuine
convolution.

Creating the breast SBSF atlas
Of interest in the SBSF concept is the distribution of the
scatter produced in the scatter body in the detector plane
when, as shown in Figure 4, the (unscattered) primary
radiation (i.e. a mini cone beam 27) is focused on a single
detector pixel 9. If this is done consecutively for each
detector pixel 9 and all the associated SBSFs 20 are summed,
the total scatter distribution is obtained for the case where
the entire detector surface is illuminated - and not only
individual detector pixels 9.
As already described above in connection with the third
requirement and step 2, the breast SBSF atlas 19 of the
scatter beam spread functions (SBSF) comprises the scatter
intensity distributions normalized to the intensity of the
primary radiation 10 in the detector pixel 9 (assuming that
the mini cone beam 27 is focused on just one pixel 9) as a
function of a plurality of different parameter
configurations:
SBSF(( λx,λy); a;H; air gap, voltage, filter, detector,... ) (#17)
and also contains the dependency of the x-ray energy spectrum
on the tube voltage, pre-filtering, radiation-sensitive
detector material, e.g. the type of scintillation crystal,
and the dependency on the presence or absence of an anti-
scatter grid and where applicable the dependence on the type
of anti-scatter grid as well as the dependence on other
parameters.
The creation of an SBSF series will now be explained:

First the parameters characterizing the relevant mammography
machine 1 are defined: SID, air gap, anode material of the x-
ray tubes (and associated emission spectra), detector
material, pre-filter material (e.g. compression plates), and
other parameters. Then comes the compression thickness H,
the voltage, the spectral filters used and other variables,
the voltage and if necessary the spectral filter (thickness)
generally being modified as a function of the compression
thickness H in order to optimize image quality.
For this parameter configuration, the parameter a describing
the tissue composition according to equation (#2a) is varied
between 0 (fat only) and 1 (glandular tissue only): the
calculation using the tried and tested Monte Carlo method
produces a set of different SBSFs 20, each α-value being
assigned an SBSF 20.
The tissue thickness H is varied between > 0 and up to
approximately 10 cm and another set of SBSFs 20 is again
calculated for each H . The voltage and the spectral filters
can also be varied, the variation being linked to H or also
independent of H . However, in the latter case multiple
variations are possible. In addition, the calculation can be
continued for all the parameter combinations.
For calculating the SBSFs 20, simplifications can be
performed which are well justified:
• Disregarding the beam divergence of the x-radiation 3 due
to the cone beam geometry by assuming approximately
parallel beam geometry; this is justified in that generally
SID >> H ; this is achieved in that the SBSF 20 remains

location- and pixel-independent for an identical beam
configuration; by identical configuration is meant that,
for each pixel, the material distribution is the same along
the mini cone beam 27 and in the lateral neighborhood.
• To improve the statistics for the Monte Carlo method and
reduce the computational complexity, pixels approximately
an order of magnitude larger (e.g. 1 xl mm2 or 2 x 2 mm2)
than the actual detector pixels 9 ( calculate the SBSFs 20; this is justified by the low-
frequency Fourier spectrum of the spatial scatter
distribution.
• The succession of fatty and glandular tissue is replaced by
a mixture; although the scatter depends (for the same
weight per unit area and path length) on whether the denser
tissue is nearer the x-ray detector 8 or nearer the
radiation source 2, according to J.M. DINTEN and J.M.
Voile: Physical model based restoration of mammographies.
In Proc. SPIE, Vol. 3336, 1998, 641-650, the differences
occurring under mammographic conditions can be disregarded.
Advantages
The solution proposed here has the following advantages:
If required, the method can be incorporated in existing
mammography machines without mechanical reconstruction.
Moreover it is a method which on the one hand shares the
adequacy of physical modeling using the Monte Carlo method,
but on the other hand - because all the time-consuming
calculations are carried out in advance where possible and

the necessary data is stored in tables - ultimately involves
relatively low computational complexity for the scatter
correction.
The modeling accuracy of the scatter correction described here
is essentially greater than that of the known (analytical)
physical models, as a number of simplifying assumptions and
approximations can be dispensed with.
The possibilities of the scatter correction proposed here go
far beyond the possibilities of the long known
convolution/deconvolution methods. Disregarding the specific
technical embodiment of the method and looking at it from a
mathematical standpoint, the method can be regarded in the
mathematical sense as a generalization of the long known
convolution/deconvolution method. On the one hand, by using
approximations and dispensing with accuracy, it can be
categorized as method of this type and then shares its
advantages, such as the possibility of using the so-called FFT
(= Fast Fourier Transformation). On the other hand, however,
the method described here can also be extended in terms of SNR
improvement, e.g. by extending the iterative multiplicative
algorithm in the direction of statistical Bayesian estimation.
In this context it should be re-emphasized that only pre-
calculation of the SBSFs 20 enables the method described here
to be performed in full generality.
Examples
Example 1:

In this example, scatter correction is performed, as described
above with equations (#5)-(#9) and (#13)-(#15), using
homogeneous location-dependent scatter beam spread functions
20 (= SBSF). For creating said scatter beam spread functions
20 it is assumed by way of simplification that the tissue
distribution characterized by the proportion α(j, k) of
glandular tissue along the beam leading from the source to
the detector pixel continues in a constantly homogeneous
manner according to Figure 4 at right angles to the beam,
i.e. parallel to the compression plates 5. It is therefore
assumed with respect to the scatter contribution of the beam
in the pixel (j, k) that the tissue composition in the lateral
neighborhood of the beam does not vary abruptly. Although
this is no longer relevant at the edge of the breast, special
treatment could be provided there.
Note, however, that the actual location-dependent
inhomogeneity of the tissue composition is allowed for by a
specifically different amount of glandular tissue α(j',k') for
each pixel (j', k') and a specific scatter contribution
dependent thereon. The SBSFs 20 are therefore generally
different for each pixel.
Example 1a:
In this example la the method is essentially performed as in
example 1.
However, the following simplifications are made:
For each pre-specified layer thickness and the other
parameters such as voltage and pre-filtering, a common SBSF
20 is used for all the pixels. In this case the SBSF 20 is

therefore selected on a location-independent basis. The
selection can be made, for example, by suitable averaging
over the tissue compositions present. ∆S in equation (#7)
and (#9) then becomes independent of the pixel index (j,k);
the double index (j,k) can - similarly as in equations (#16a)
to (#16c) - be omitted.
The important feature is that the integral in equation (#9)
becomes a genuine convolution which can be efficiently
executed by FFT (=Fast Fourier Transformation).
Example 1b:
In this example lb the method is likewise performed
essentially as in example 1.
In this case, however, a uniform convolution kernel (for all
the layer thicknesses) is used for the scatter calculation.
The fact that for a small layer thickness relatively less
scatter is produced than for a large layer thickness must be
taken into account by means of scaling factors which are a
function of the layer thickness and other parameters such as
voltage and filtering.
Approximately the same computational complexity is necessary
for example 1b as for example la. On the other hand, much
less memory space is required for storing the breast SBSF
atlas 19 in this example.
Notes on examples la and lb:
In general the simplified examples la and lb share the
characteristic that the convolutional models for the scatter

can be inverted using the Fourier transformation. This is
known as deconvolution. The examples described here differ
from the conventional deconvolution methods in using one or
more scatter beam spread functions 20 obtained in advance by
Monte Carlo simulation.
With regard to performing deconvolution, reference is made to
a publication by J.A. SEIBERT and J.M. BOONE: X-ray scatter
removal by deconvolution. In Med.Phys., Vol. 15, 1988, pages
567 to 575. Reference is also made to the more recent
publication P. ABBOTT et al: Image deconvolution as an aid to
mammographic artifact identification I: basic techniques. In:
Proc.SPIE, Vol. 3661, 1999, pages 698 to 709 which deals with
deconvolution using regularization techniques for noise
suppression. Another deconvolution method with thickness-
dependent convolution is described in D.G. TROTTER et al:
Thickness-Dependent Scatter-Correction Algorithm for Digital
Mammography. In: Proc. SPIE, Vol. 4682, 2002, pages 469 to
478. In this method an iteration with relaxation is
performed.
Example 2
In this example the method is essentially performed as in
example 1, but employing scatter beam spread functions 20
which have been calculated for an inhomogeneous medium.
Figure 5, for example, illustrates the case where a breast
region 28 has a different composition from that of a
surrounding breast region 29.
This enables it to be taken into account that the SBSF 20
depends not only on the tissue composition along the mini

cone beam 27 supposedly focused on the detector pixel 27 but
also on the tissue composition in the lateral neighborhood
into which photons are scattered and can be further scattered
again in the direction of the pixel. However, the effective
extent of the lateral neighborhood is not very large because
of the average free path length mammography energy range between about 20 and 40 keV. It
would therefore suffice to assume the tissue composition to
be homogeneous in a lateral half space, but generally
different from the mini cone beam 27. The allowance for
inhomogeneous SBSFs 20 with differences between beam and
neighborhood might be relevant particularly at the breast
edge .
This example therefore constitutes a generalization of the
above-described examples 1, la and lb, as in this case the
SBSFs 20 depend not only on a tissue parameter a, but also
on a surrounding tissue parameter y to be newly introduced.
In this case the breast SBSF atlas 19 would therefore have an
additional dimension.
For the sake of clarity, the following table compares the
different characteristics of examples 1, la, lb and 2:

The method described here can also be applied to so-called
dual energy methods which will be known to the average person
skilled in the art. With the so-called dual energy method,
which is used primarily in mammography or in bone
densitometry, images are recorded simultaneously using two
different energy spectra. The recordings using different
energy spectra are provided by two different voltages and if
possible also different spectral filtering so that the
spectral regions effectively corresponding to the two
measurements overlap one another as little as possible. By
means of a computational process which is essentially based
on the solution of a generally nonlinear system of two
equations assigned to the two spectra, finer tissue
differentiation can be achieved compared to a recording using
one energy spectrum. For computation to be successful, the
scatter components must be eliminated as much as possible, as
otherwise the artifacts induced by the scatter components are
in some circumstances stronger than the actual tissue image.
Because of the differences in scatter for the two spectra,
effective scatter correction is therefore critically
important for the quality of the dual energy method.
The proposed scatter correction method can also be used in
this context. The geometrical parameters are identical for
the two recordings, but the spectrally dependent parameters
are different.
The correction must be carried out for each of the two
recordings according to the described formula, the only
difference being that different SBSFs 20 must be used
according to the different spectra.

WE CLAIM
1. An apparatus for projection radiography, in particular for mammography,
comprising a radiation source (2) which emits radiation (3), a detector (8)
and, downstream of said detector (8), a processing unit (12) which
determines in an approximately manner, on the basis of the projection
data (17) supplied by the detector (8), the scatter material distribution of
the object under examination (6) and reads out scatter information (20)
from a data memory (19) as a function of the scatter material distribution
and corrects the projection data (17) in respect of the scatter content (11)
on the basis of the scatter information (20),
characterized in that
- the scatter information (20) is determined by Monte Carlo simulations
which calculate the interactions of the photons with a scatter material
distribution in respect of different scatter material distributions, in that
- the scatter information is scatter distributions (20) which describes a
scatter-induced distribution over adjacent image regions of the radiation
(14) emitted by the radiation source (2) and directed to a particular image
region, and in that
- the processing unit (12) determines a scatter distribution (11) in an image
region of the projection image (17) by calculating and adding, for each
image region, the scatter contributions (11) of the surrounding image
regions.

2. The apparatus as claimed in claim 1, wherein the scatter distributions (20)
are scalable using the intensity of the unscattered primary radiation (10)
incident on the detector (8).
3. The apparatus as claimed in claim 1 or 2, wherein the processing unit (12)
evaluates scatter information (20) specific to the particular scatter
material distribution for different image regions of a projection image
(17).
4. The apparatus as claimed in one of claims 1 to 3, wherein the processing
unit (12) determines a scatter distribution (11) in an image region of the
projection image (17) by convolving the primary radiation distribution with
a scatter distribution (20).
5. The apparatus as claimed in one of claims 1 to 4, wherein the processing
unit (12) determines the unscattered primary radiation (10) by solving the
implicit equation P + S (P) = T, where P is the distribution of the
unscattered primary radiation (10), S (P) the secondary radiation
distribution (11) dependent on the unscattered primary radiation (10) and
T the measured total radiation distribution in the projection images (17).
6. The apparatus as claimed in one of claims 1 to 5, wherein processing unit
(12), for a first approximate determination of the scatter material
distribution of the object under examination (6), estimates the amount of
scatter (11) contained in the projection image (17) on the basis of scatter
information associated with a typical scatter material distribution.

7. The apparatus as claimed in one of claims 1 to 6, wherein the processing
unit (12) performs the processing steps iteratively as claimed in one of
claims 1 to 6.
8. The apparatus as claimed in one of claims 1 to 7, wherein the scatter
information (20) is calculated under the assumption of a homogeneous
scatter material distribution in the beam direction.
9. The apparatus as claimed in one of claims 1 to 8, wherein scatter
information (20) allowing for the outer contour (26) of the object under
examination (6) is stored in the data memory (19).
10.The apparatus as claimed in one of claims 1 to 9, wherein scatter
information (20) calculated under the assumption of a homogeneous
scatter material distribution perpendicular to the beam direction is stored
in the data memory (19).
11.The apparatus as claimed in one of claims 1 to 10, wherein scatter
information (20) determined allowing for an inhomogeneous scatter
material distribution perpendicular to the beam direction is stored in the
data memory (19).
12.The apparatus as claimed in one of claims 1 to 11, wherein scatter
information depending on parameters of the radiation source (2) is stored
in the data memory (19).

13.The apparatus as claimed in one of claims 1 to 12, wherein the processing
unit (12) determines the scatter components (12) at selected grid points
and works out the correction values for individual detector elements (9)
by interpolation between the grid points.
14. The apparatus as claimed in one of claims 1 to 13, wherein the object
under examination (6) is compressible in a compression device (5) and
that the processing unit (12) uses the physical configuration of the
surfaces of the compression device (5) facing the object under
examination to determine the path length of the radiation (3) through the
object under examination.

ABSTRACT

TITLE: "AN APPARATUS AND A METHOD FOR SCATTER CORRECTION IN
PROJECTION RADIOGRAPHY"
The invention relates to an apparatus for projection radiography, in particular for
mammography, comprising a radiation source (2) which emits radiation (3), a
detector (8) and, downstream of said detector (8), a processing unit (12) which
determines in an approximately manner, on the basis of the projection data (17)
supplied by the detector (8), the scatter material distribution of the object under
examination (6) and reads out scatter information (20) from a data memory (19)
as a function of the scatter material distribution and corrects the projection data
(17) in respect of the scatter content (11) on the basis of the scatter information
(20), the scatter information (20) is determined by Monte Carlo simulations
which calculate the interactions of the photons with a scatter material
distribution in respect of different scatter material distributions, in that the
scatter information is scatter distributions (20) which describes a scatter-induced
distribution over adjacent image regions of the radiation (14) emitted by the
radiation source (2) and directed to a particular image region, and in that the
processing unit (12) determines a scatter distribution (11) in an image region of
the projection image (17) by calculating and adding, for each image region, the
scatter contributions (11) of the surrounding image regions.

"AN APPARATUS AND A METHOD FOR SCATTER CORRECTION IN PROJECTION RADIOGRAPHY"

Documents:

Inventors:

PCT Conventions: