WO2011036436A1 - X-ray imaging - Google Patents

Abstract

Using the acquisition configuration of a mammography unit and an x-ray source model, the spectrum incident upon the breast is known. Using the pixel intensity of the mammographic image and the detector model, the energy exiting the breast is known. The contribution of scattered radiation is calculated and subtracted from the acquired image. The thickness of reference material that is required to match the attenuation of the breast is therefore ascertained by numerically solving for the thickness which when traversed by the incident beam, and subsequently recorded by the image receptor results in the primary component of the pixel intensity observed in the mammographic image.

Description

X-RAY IMAGING

The invention relates to methods and apparatus for use in normalisation and quantification of x-ray projection images, for example of the human breast, including estimation and correction of the scattered radiation.

Breast cancer is one of the leading causes of cancer and mortality of women, particularly in developed countries. The earlier the diagnosis, the better the prognosis, and this has led both to greater awareness of breast cancer and to the introduction of (national) screening programmes in many countries. Currently, one of the most important technologies used to detect breast cancer is mammography - an x-ray image of the breast. It is estimated that 60 million mammograms are taken annually, a number that is rising quickly, especially as mammography becomes more regularly used in Asia.

Traditionally mammograms result in a film that is viewed on a light box. The early technology involved converting the x-ray photons emerging from the breast into light photons, which then expose a film. There has been increasing interest in computer-based analysis of mammograms, which, in the case of mammographic films, requires that they be digitised in order to be read in to the memory of a computer. Over the past 15 years, a number of digital x-ray detectors have been developed, either by converting x-ray photons to light and digitally recording the light intensity, or by converting x-rays directly into electrons, thereby allowing the incident x-ray energies to be read directly onto a computer. In both cases (a digitised film or a directly-digital mammogram) the image is represented as a large array of pixels, typically of a size 50 μηι (microns) square (for the digitised mammogram film) or between 70 and 100 μηι square for a directly-digital device. Since the present invention involves computer processing of mammographic, tomographic, or computed tomography images, we assume that the image is stored in the memory of a computer and that it is represented as an array (or, in the case of digital breast tomosynthesis arrays) of pixels.

2D Mammograms The overwhelming majority of mammograms are two-dimensional (2D) projected images of the tightly-compressed breast. In such a mammogram, the brightness recorded at each location in the image encodes the attenuation (roughly: absorption, though scattering effects also require consideration and is a subject of this invention) of the photon flux along the path from the focal point of the x-ray source to the pixel area. In fact, the brightness recorded at the image location depends on many factors, including: the specifics of the x-ray source, including the accelerating potential and the anode material; the exposure time; the way the flux is filtered before exposing the breast; the detector characteristics; and any anti-scatter devices present. In addition, and crucially, the brightness also depends upon the unknown contents of the breast that is being imaged, since tissues of different varieties attenuate the incident beam to differing extents. As a result of this complex dependency of the brightness on a multiplicity of factors, compounded by the unknown contents of the breast, most analysis - by human or computer - has been, and continues to be, qualitative rather than quantitative, that is, it relies upon human judgement, which, even for experts, shows substantial inter- and intra-observer variability. More precisely, such analysis rarely yields measurements that are either absolute and/or which enable detailed comparisons between subjects, or between images of the same subject, over time, view, or between left and right breast, or the estimation of breast density. For the same reason, computer analysis of

mammograms - for example, to detect cancer or the change of breast contents over time - is made more difficult.

A 2D mammogram is a projected view of an intrinsically 3D structure. The brightness at an image location represents the integral of the attenuation of the beam along the line from x-ray source to pixel (plus a component due to x-ray scatter). A consequence is that a region with a value that is substantially brighter than its surrounds may be the result of a single dense region (for example, a tumour) or the result of two overlapping normal regions whose sum appears abnormally bright. For this reason, the breast is tightly compressed (by a force of typically around 100- 150N) to spread tissue thinly and avoid as far as possible the problem with overlap, as well as to reduce patient dose and to maintain high signal-to-noise ratio between varying tissues. A clinician (or a computer program) aims to draw inferences from the images about the breast contents. This is intrinsically difficult, since there is huge variation in the tissue constituent of breasts across the population, hence huge variation in appearances of mammograms. Also, the signs of disease are often quite subtle, adding to the difficulty of the inference process. However, in the case of

mammography, the problem is made worse by the confounding effects referred to above: the specifics of the incident radiation; detector characteristics and the nature of the anti-scatter devices (as well as scatter itself)-

Whilst the majority of the incident x-ray radiation traverses the breast along a straight line (referred to as the primary component), a not insignificant proportion is scattered, that is the photon direction (and possibly energy) changes due to interactions with the atoms of which the breast tissue is composed. The change in direction results in a complex path between the focal spot and image receptor, and, as a result, no useful information is carried by scattered radiation, since its intensity no longer describes the attenuation properties of the tissue in the straight traversal between the focal spot and the image receptor pixel. In this way, scatter reduces image signal by introducing a low frequency blurring into the image, and hence limits the ease with which a human (and often a computer) observer can perceive and recognise a feature. To limit this effect, an anti-scatter grid is almost always incorporated between the breast and the image receptor in current clinical mammographic practice, though in tomosynthesis the movement of the focal spot introduces complications and so anti-scatter grids may have to be jettisoned for tomosynthesis. The grid is typically formed of a series of thin strips of a highly attenuating material (generally lead) orientated toward the focal spot in either a cellular or linear arrangement, and interspaced (at a greater distance than the thickness of the strips) by a low attenuating material (generally a cellulose based material). Such grids attenuate the photon flux incident upon them according to angle of incidence, and since scatter is generally of a far lower angle (measured relative to the detector surface) than that of the primary flux, it tends to be selectively rejected by the grid. To prevent the grid from being visible within the image it is necessary for it to move in a reciprocating manner in such a way that each image pixel has an attenuating strip falling directly above it for an equal time as each of the others. Unfortunately, scatter grids are not 100% efficient either in rejecting scattered photons, or in transmitting primary ones. The 'Bucky factor' quantifies the inefficiency in the transmission of the primary in terms of the multiplicative factor by which the exposure (and therefore harmful x-ray dose to the women) must be increased to obtain the same photon flux at the image receptor (and thereby magnitude of quantum noise) as that which would result in the absence of the grid. Bucky factors are of the order of 2, although values as high as 2.35 are known.

3D Mammograms: Digital Breast Tomosynthesis

A recent development in mammography, which has not yet resulted in commercially available products, but which has attracted a great deal of research interest, is digital breast tomosynthesis: a series of projection images (perhaps as many as 25) are taken at slightly different angles to the breast (in one realisation of the principle these views are 1° apart). The separate views, taken over a limited angle, are then combined in a process known as reconstruction to form a 3D representation of the breast. The depth resolution attainable depends upon the angular range over which the individual views are taken. Reconstruction is mathematically very different from tomography (as used, for example, in CT or PET) because the angular range over which tomosynthesis views are taken is very much less than 360°, and thus a considerable mathematical null space of reconstruction exists. Note that in tomosynthesis the signal-to-noise ratio of each of the (say, 15) images is considerably poorer than in the case of conventional mammography, because the total "x-ray budget" available cannot exceed that for conventional mammography. As a result the scatter component of individual tomosynthesis projections is significantly (perhaps an order of magnitude) greater than in mammography.

Spectroscopic x-ray detectors

Much research effort has been recently focused on the development of photon counting detectors. These work by processing and recording each and every incident photon as an individual entity, rather than simply integrating across all photons within the incident flux to give a total energy. Such detectors produce very low noise, low dose images, and are beginning to be used for mammography. This technology is further being developed to measure the incident x-ray spectrum at each pixel, and thus a spectral image is recorded. Specifically, each pixel of a spectral image records a vector of photon fluxes for each of a number of discrete photon energy ranges, rather than a single pixel intensity describing the flux over all photon energies present in the incident beam. One example of such technology is the pixelised silicon drift Medipix detector developed by an international collaboration hosted by CERN in conjunction with the Large Hadron Collider experiments. Such technology facilitates the addition of an "extra dimension" to the radiodensity measurement disclosed herein. For such detectors, the present invention as defined below, may be applied to each photon energy range in the vector recorded by the detector, and thus the polyenergetic values available from conventional detectors, may be replaced by a series of approximately monoenergetic radiodensity measures. As a result, the atomic and density variations in tissue (or calcification) may be quantified with far greater resolution, opening significantly greater opportunities to characterise tumour tissue, and the surrounding micro-environment. Such detectors are dealt with by the invention through a specialisation of the generalised detector model disclosed.

Summary of the Invention

The present invention provides methods, apparatus, and a computer program as defined in the appended independent claims.

The invention comprises a novel method, which may be embodied in the form of a computer program. Preferred embodiments of the invention calculate a normalised scatter compensated measure of the radiodensity per unit distance of tissue traversed between the focal point of an x-ray source and each pixel in an x-ray projection image of a mammal's breast: that is the radiodensity per unit distance encountered by the primary beam. The radiodensity may therefore be thought of as a multiplicative factor that must be applied to the equivalent thickness of reference material to glean an equal attenuation to that observed. The invention potentially overcomes the problems discussed in the preceding paragraphs describing the background to the invention, and opens the way to quantitative image analysis of mammograms, whilst at the same time producing low patient dose images (because the scatter grid may be removed) that are free from the degrading effects of scatter. The quantitative analysis facilitated by the present invention adds a "new dimension" to mammography, which traditionally only considers the relative contrast between image features; instead, as a result of processing the image according to the present invention, an absolute measure of tissue radiodensity is made available. This can be used to quantify the attenuation (and hence atomic and density characteristics) of a suspect feature, the microenvironment of a rumour, or the breast as a whole, and then may be compared to that which is known for such a lesion or a cancerous breast. For example, in the case of a suspect microcalcification, identified either by a human reader or computer, the radiodensity measurement allows the similarity of the x-ray attenuation of the feature to be compared to that of the likely calcium based compound present (typically calcium hydroxyapatite in the malignant case, and calcium oxalate dihydrate in the benign case), thereby reinforcing the diagnosis. Current imaging techniques used for diagnosis go no further than identifying a single microcalcification as a "small spot of high brightness", the extra information as to the constituent calcium salt identified as a result of processing the image using the disclosed invention being only available after invasive biopsy.

The invention significantly enhances and extends two inventions in the prior art: the hj_nt representation for mammographic image processing, developed by Highnam and Brady; and the Hounsfield unit, a standardised and accepted unit for reporting and displaying reconstructed x-ray computed tomography (CT) values.

The computation of the hj_nt representation consists of using a model of the physics of image formation to ascertain the thickness of "interesting tissue" (breast tissue that is anything other than fat) above each pixel in the mammographic image. Primarily as a result of histopathological analysis of excised breast tissue, and its comparison with its appearance in a mammographic image, a significant contribution to which comes from a study conducted by Alowami et al, this two category classification of tissues has been found to be limited due to a number of factors, the most significant being the variation between members of the population in the composition, and hence x-ray properties, of Highnam and Brady's two "constant" tissue types. In particular, the inclusion as a single constant entity of both the functional tissue, such as the glandular aggregates, and the structural tissue, such as the connective stroma, has been demonstrated by Alowami et al's findings to be counter-productive, since they conclude that it is the variation in the composition of the stroma, for example the proportion of collagen present, that is the significant factor in dictating

mammographic appearance. The h(_nt representation aims solely to measure the quantity of "interesting tissue" (which includes stroma), and the inter-population variations in stromal composition identified by Alowami et al to be of importance, violate Highnam and Brady's underlying assumption of constant x-ray properties within "interesting tissue" across the population, and as a result the stromal composition variation is lost by their technique. The disclosed invention overcomes this shortcoming by adopting a novel representation of breast x-ray properties, that of measuring radiodensity: the property of relative transparency compared to some reference material of the passage of x-ray photons through the tissue. A measure of radiodensity is not limited by assumptions as to constituent tissues, since it solely encodes the x-ray absorption properties. This approach is analogous to the

Hounsfield unit, in which attenuation is quantified relative to water, and thus acts as a reference in CT images. Breast imaging presents significant additional challenges to those encountered in imaging other parts of the body since it is fundamentally concerned with identifying differences in soft tissues, all of which have very similar x-ray attenuation properties, thus the discriminatory signal is small and the signal-to- noise ratio poor. The technique adopted for the calculation of the Hounsfield unit at any given image point consists of a comparison of the intensity with that observed in an image of a known quantity of water, which therefore ignores beam hardening and scattering effects that vary between the reference and anatomical image.

Mammographic techniques based on simple comparison with an image of a reference object are reported in the literature, though the reference material of choice varies. In order to improve accuracy and to take consideration of the limited signal-to-noise ratio in soft tissue imaging, the disclosed invention significantly enhances the standard approach of simple comparison for ascertaining radiodensity. The first enhancement in the disclosed invention comes from the computation of radiodensity per unit distance, rather than simply the radiodensity, of the tissue traversed between the focal point of an x-ray source and each pixel in an x-ray projection image. The traversal distance through the breast is estimated using a marnmography-specific ray tracing algorithm. Dividing the thickness of reference material required to produce an equal attenuation to that which is observed for the focal spot-detector pixel traversal path under consideration, by the length of said traversal path, gives the normalised measure of radiodensity per unit tissue traversal distance: a direct measurement of the attenuation properties of the tissues present, suitable for intra and inter patient comparison since it removes the dependency on breast thickness, and thus the magnitude of the compressive force used for image acquisition.

The second enhancement comes from the use of a model of the physics of image formation to calculate radiodensity, rather than relying on simple comparison of pixel intensity with those observed from calibration images. This not only allows the calculation of radiodensity per unit traversal distance, it affords better accuracy by fully considering beam hardening effects associated with variation in compressed breast thickness, and also allows the full consideration of scattered radiation and thus a suitably compensated measure of radiodensity solely arising due to the primary beam. This is similar to the Highnam and Brady's hint representation, which is also calculated using a model of the physics of image formation, however that which is disclosed herein is a considerably enhanced reworking of said method, giving far greater consideration to the underlying physical phenomena, resulting in a method devoid of the majority of Highnam and Brady's underlying simplifying assumptions as well as utilising optimal information sampling so as to facilitate a fast computation time. During testing the implementation of the disclosed invention achieved a computation of the scattered radiation in 128.71 seconds, running in a single thread of execution on a Intel Core2 6600 CPU at 2.40GHz, for the image shadow of a phantom with a base area of 10 000mm², imaged at 0.07mm pixel spacing, whilst returning an accuracy within an upper error bound of 10%. Such accuracies are of the order previously only associated with Monte Carlo photon transport modelling techniques, which require orders of magnitude more computing time (months). The disclosed invention therefore offers the unique combination of high accuracy, and fast computation time.

The x-ray beam produced by the majority of mammographic x-ray source, for example a vacuum tube, is polyenergetic, that is there is a complex spectrum of photon energies covering a considerable range of energy values. This further complicates the analysis of x-ray image formation. Often, particularly in the case of current digital breast tomosynthesis reconstruction algorithms, the beam is assumed to be monoenergetic (akin to being a laser); but while this simplifies the

mathematical analysis, it creates substantial artefacts, for example apparent beam hardening. In one embodiment of the invention, the reference material has been chosen to be an equal mixture by mass of published attenuation values for breast fibroglandular and adipose tissue (the two tissues for which the breast is often considered to approximately consist of, particularly in medical physics work); but it is understood that other choices of reference material can be made and may be preferable in certain circumstances. The particular choice of reference material used in the current embodiment has the favourable property of closely mimicking the attenuation properties of the breast across the range of mammographic photon energies (since it is based on data directly measured from breast tissue) and therefore results in a measurement with very low dependency on beam quality. Any dependency arises from the natural variation in composition of these tissues between the samples measured and those in the population at large, and the proportion of these tissues present in the breast under examination.

Part of the calculation of the radiodensity measurement at each projection image pixel involves the calculation of the contribution to the image pixel of each voxel encountered within the breast volume. In cases where multiple projection images acquired from different focal spot positions and/or detector orientations are available, most particularly in the case of digital breast tomosynthesis, this output may be used to form a set of simultaneous equations, which can be solved to yield a novel method for reconstruction of a three-dimensional voxelised radiodensity image of the breast. It will be recalled that the scatter (non-primary) component of the mammography image is a further complicating factor in analysis of mammograms (by human or by computer program). In the present invention, the effects of scatter are explicitly taken into consideration, so the measurement of radiodensity may be thought of as the attenuation encountered by the primary radiation traversing the breast. The result of the process is a direct measurement of the attenuation characteristics of the tissue present, and is thus suitable for intra and inter patient comparison. Within the consideration of scatter, the contrast limiting effects are also modelled and removed, thereby allowing the acquisition of low dose images (reduced by the bucky factor) without an anti-scatter grid present. The invention therefore facilitates the removal of the harmful dose increase associated with the scatter grid, while at the same time maintaining image signal and thereby aiding the perception and recognition of image features. It also removes the complication of fabricating the anti-scatter grid, a non-trivial process due to the high accuracy and small size of the grid components, and also the reliability issues associated with a scatter grid, in particular those' associated with the moving parts required to facilitate reciprocation. The invention will also benefit digital breast tomosynthesis, for which scatter is a more important image degrading effect and for which anti-scatter grids are difficult to use.

Embodiments of the invention will now be disclosed, by way of example only, with reference to the accompanying drawings in which: -

Fig. 1 is a flow chart (Flow Chart 1) showing a high level overview of the process for calculating the radiodensity image from a raw x-ray projection image;

Fig. 2 is a flowchart (Flow Chart 2) giving an overview of the vacuum x-ray tube model for calculating the photon beam departing the tube window for a given image detector pixel;

Fig. 3 illustrates the x-ray tube geometry in which the published x-ray spectra are measured;

Fig. 4 is a schematic illustration of a typical x-ray tube anode design of the type encountered in clinical use for mammography; Fig. 5 is a flow chart (Flow Chart 3) showing an overview of the ray tracer used to calculate the list of voxels encountered by a ray, and the distance traversed through each of them;

Fig. 6 is a flow chart (Flow Chart 4) showing an overview of the process for calculating the attenuation of a ray between the point of entry into an object, or point of photon scattering, and an image detector pixel;

Fig. 7 is a flowchart (Flow Chart 5) showing an overview of the process for calculating the scatter incident on a given image pixel arising from the traversal of a given primary ray;

Fig. 8 is a schematic overview of the process of x-ray scattering arising from a primary ray;

Fig. 9 is a flow chart (Flow Chart 6) showing an overview of the process for calculating the scatter image;

Fig. 10 is a flowchart (Flow Chart 7) showing an overview of the process for calculating the scatter incident on the pixels within the immediate neighbourhood surrounding the pixel upon which a primary beam is incident for the homogeneous equivalent breast;

Fig. 1 1 is a flow chart (Flow Chart 8) showing an overview of the image filtering applied to give consideration to the inhomogeneous mix of tissues present within the breast;

Fig. 12 is a flow chart (Flow Chart 9) showing a high level overview of the process for calculating the transfer function from the empirically sampled data;

Fig. 13 is a flow chart (Flow Chart 10) showing a procedure for a technician to calibrate an x-ray mammography machine;

Fig. 14 is a schematic illustration of a possible construction of a breast equivalent attenuator for calibrating a mammographic image detector according to an embodiment of the invention; Fig. 15 illustrates a voxel model of the breast utilised for reconstruction of 3D radiodensity volumes from a number of projection images obtaining from varying focal spot positions;

Fig. 16 is a schematic cross-section illustrating the increasing thickness of breast tissue to be traversed by the x-ray beam as the sampling angle is increased in tomosynthesis;

Fig. 17 is a graph showing the variation in radiodensity with calcification diameter, for high resolution images in which the calcification is completely sampled according to the Nyquist-Shannon sampling theorem, and low resolution (typical of current mammographic image detectors) in which under-sampling is occuring;

Fig. 18 shows photographs of a test object manufactured from tissue- equivalent resins, and an x-ray apparatus used to test the model;

Fig. 19 is a graph showing the pixel intensity of the raw images for a range of step thicknesses of the test object for a range of clinical tube potentials, exposures and target/filter combinations;

Fig. 20 is a graph of the resulting radiodensity measures following normalisation of the pixel intensity data of Fig. 19 using an embodiment of the invention;

Fig. 21 is a schematic illustration of a computer system according to an embodiment of the invention.

In order to calculate radiodensity, that is the attenuation encountered by the primary radiation traversing the breast per unit tissue traversal distance, the difference between the primary x-ray photon flux incident upon the breast; and the primary x-ray photon flux exiting the breast and being subsequently recorded by the image detector, is calculated using a detailed model of the physics of mammographic image formation. This model consists of three components which will be discussed in detail in what follows: a model of the x-ray source; a model of the image detector; and a model of the attenuation and scattering within the tissues of the breast. Flow Chart 1 in Fig. 1 provides a schematic representation of the radiodensity calculation. Since the radiodensity measure is free from the effects of scattered radiation, the scatter image calculated using the model of scattered radiation (disclosed in what follows) must be subtracted from the acquired energy absorption image (the calculation of which is disclosed in the detector model which follows) in order to ascertain the image resulting from the primary signal only. For a given projection image, using a detector that records only a single photon flux measurement at each pixel, the radiodensity is calculated using the image formation model to ascertain the thickness of reference material, per unit traversal distance of breast tissue, that leads to the primary energy absorbed by the detector for the pixel intensity observed. This function is non-linear, so a possible solution technique is the use of a numerical minimiser. However, the variation of the radiodensity measure with detector output is smooth, and therefore where accuracy permits, it may be approximated using a suitable function and thereby an analytical solution may be used.

Technological advancements such as silicon drift detectors, the technology utilised in the pixelised Medipix detector for example, have brought about the possibility of the use of x-ray spectrometry in mammographic imaging. X-ray spectrometry allows far greater specificity in the characterisation of the radiodensity of tissues, and thus in principle provides a far better metric from which to characterise important features of the breast, such as the microenvironment surrounding a tumour. In an x-ray spectrometry image a set of photon fluxes is recorded at each pixel, each describing the flux over a specific range of photon energies, rather than the sum of all the photon energies present in the beam as a whole. Several alternatives exist to apply the radiodensity measure disclosed in this invention to normalise spectrometry images; one embodiment normalises each of the photon energy ranges recorded using a single reference material; another

embodiment uses a mixture of a number of reference materials (two or possibly more depending on the exact application) for which a weighting coefficient for each material is calculated which best fits the x-ray spectrum recorded. In an identical fashion to that discussed previously in the case of measuring only a single photon flux at each pixel, the radiodensity function is non-linear, so a possible solution technique is the use of a numerical minimiser. Again also, the variation of the radiodensity measure with detector output is smooth, and therefore where accuracy permits, it may be approximated using a suitable function and thereby an analytical solution may be used.

X-Ray Source Model

The x-ray source model is used to calculate the photon beam, that is the spectrum of x-ray photons (the number of photons present in the beam for each of the discrete energy ranges considered) that emanates from the source window in any given direction. The various components of the disclosed invention which contribute to the calculation of radiodensity are equally applicable to any x-ray source technology, providing a suitable model of said source is developed to calculate the spectrum of x-ray photons emitted. The most common in current clinical use is the vacuum x-ray tube, however, flat panel parallelised sources based on pyroelectric electron emission followed by subsequent bombardment of a suitable target is another of a number of possibilities. In one embodiment of the invention, the x-ray source is taken to be a vacuum x-ray tube, and for which Flow Chart 2 in Fig. 2 provides a graphical overview of the constituent steps in the method for calculating the emitted photon spectrum.

Emission Spectral Data

A vacuum x-ray tube operates by bombarding an anode (of a carefully selected material) with electrons thermionically emitted by a cathode (the

accelerating potential between the two being what is often termed the tube voltage) in order to produce x-ray photons. The resulting spectrum is therefore governed by the emission characteristics of the anode material, which is itself governed by the electron configurations within the constituent atoms. In one embodiment of the invention, the emission spectra of the tube anode is based upon published data (e.g. K. Cranley, B. J. Gilmore, G. W. A. Fogarty, and L. Desponds, Catalogue of diagnostic X-ray spectra and other data. York: Institute of Physical Sciences in Medicine, 1997) describing the photon energy spectrum (where energies are assigned to discrete intervals, called "bins") arising from a given anode material across a range of target angles and tube potentials. It is generally accepted that the relationship between the logarithm of the measured photon count, for a given tube potential, is approximately linear with the reciprocal of the tangent of the target angle. For each photon energy bin a regression fit, or similar numerical method, can be made to establish the coefficients of this function, and it thereby provides a high accuracy interpolation technique for ascertaining spectra for target angles falling between the relatively coarse set of angles available in the published data.

Consideration of the Anode Heel Effect

The position within the anode at which x-ray photon production occurs varies for the generally thick anodes found in mammography. As a result, the incident electrons and the emitted x-ray photons are subjected to varying amounts of attenuation (filtration) by the anode material itself. Since emission is polyenergetic, and lower energies are attenuated to a greater extent than higher energies, the varying filtration leads to a variation in both magnitude and photon energy distribution of the spectrum incident upon each spatial position on the image receptor. This effect is known as the anode heel effect, and is incorporated in the invention through the calculation of an effective target angle for each image pixel. The effective target angle is the target angle in the measurement geometry of the published data for which the self filtration distance within the anode is the same as that for the ray emanating in the direction of the pixel for the mammography unit under

consideration, and therefore includes such factors as the anode angle and the angle of tilt between the x-ray tube and the image receptor.

Calculating Effective Target Angle

One of the possible methods of calculation of the effective target angle, ε , for any 3 dimensional emission direction may be calculated as follows. Previously published data (Cranley et al.) uses a basic two-dimensional tube geometry as shown in Figure 3.

Let P_M be the photon origination point, Y the point at which photon exits the anode, and p the electron penetration depth at which x-ray photon production occurs. The self filtration of the x-ray photon beam emerging on the spectrum measurement axis (90° to the electron beam), for an anode angle of θ , is therefore the distance - P_M | , which may be calculated from the electron penetration depth using trigonometry: J = - ^N

tan <

A generalised design of a mammography tube that is found in widespread clinical use is illustrated in Figure 4 (although depicted in 2D it should be noted that the x-ray emission direction is 3D).

Let P be the x-ray photon origination point, X the x-ray photon anode exit point, plane AB [n_AB,p_AB] be the anode face, and let line PC |d_PC, pj be the photon path. In this, and the following, we use the familiar notation [i, A] to denote the unit vector x that passes through the point P . For completeness, we introduce several tube design parameters: β is the tube inclination to the image receptor, Θ is the anode angle (usually zero for mammography), and φ is the angle between the emitted x-ray photon beam and the anode surface. In this tube geometry, the self filtration of the x-ray photon beam emerging in the direction PC, is the distance

|X - P| , where the position of the point X is the intersection of the photon path PC

[d_PC,pj with the target face AB [n_AB,p_AB], It is straightforward to show then that

Tube Accelerating Potential Waveform

Variations are sometimes present in the potential difference between the anode and cathode during an exposure (known as 'ripples'). Such variations have a resultant effect on the electron energy incident on the anode, and hence on x-ray photon emission. These are generally a result of the type of mains supply and rectification used. Such variations are incorporated within the present invention by calculating a time of occurrence weighted sum of a number of spectra generated across a range of constant tube potentials. Attenuation of Photon Beam for Tube Window and Beam Filtration

The beam emanating from the anode surface is suitably attenuated for absorption effects within the tube window and any additional beam filtration included by the manufacturer to suitably 'shape' the x-ray energy spectrum. The attenuation depends on the traversal distance through such filtration materials, and this is calculated according to the angle of incidence of the ray of photons. An intensity scaling can also be applied to include consideration of the inverse square law effect between the focal spot and the detector pixel position.

Empirical Calibration Corrections

Finally, two empirical calibrations are applied within the embodiment of the invention for beam intensity and for beam quality. These take account of discrepancies between the specific tube and the idealised spectra, and include such factors as manufacturing tolerances and anode wear.

The beam intensity calibration in one embodiment of the invention is obtained by scaling the entire spectrum by a constant such that the beam intensity measured in terms of air kerma matches that measured on the specific clinical machine.

The beam quality correction is obtained through the addition of filtration to the tube output of the thickness required to match the first half-value layer of the model with that measured on the clinical machine. In one embodiment of the invention, aluminium is the material in which the first half value layer is measured, and also that used for additional filtration.

Ray Tracing

In the current embodiment of the invention, a novel ray-tracing algorithm has been developed to calculate the details of the traversal path taken by each of the primary ray paths from the focal spot to each image detector pixel, and likewise those of the scattered rays. Specifically, the ray tracing algorithm calculates a list of materials/tissues encountered by the beam of x-ray photons during its traversal for which each list entry contains: the name of the material; the three-dimensional location of the point of entry of the beam; the three-dimensional location of the point of exit of the beam; and the distance traversed through the material. In one embodiment of the invention the outer surfaces of the compressed breast and the various components of the x-ray acquisition equipment that are present within the x-ray beam (for example the compression plate) are represented as vector surfaces. For example, the breast compression plate, which for certain equipment may be modelled as a cuboid, has its upper and lower surfaces modelled as vector planes defined by:

(r - a) n = 0 where r is a point on the surface (plane), a is a given point on the surface (plane), and n is a vector normal to the surface (plane).

In one embodiment of the invention the ray of photons incident upon a given pixel on the image receptor is modelled as a line in vector space defined by: r = b + A(c - b) where r is a point on the straight line joining points b and c at a distance λ from c .

The entry and exit points in three dimensional space of a ray of x-ray photons is calculated using the vector models, in the example given of the cuboid

compression plate, the intersection of the ray (line) with an outer surface of the cuboid results from solving:

(b + A_int (c - b)- a) n = 0 for /l_jnt which gives the position along the ray of the point of intersection, and the three-dimensional position, r_int , is given by: rin, = *> + ,l_int(c - b)

Modelling the Breast Tissue Volume

In one embodiment of the invention, the volume corresponding to the breast is represented as an array of small voxels whose shapes are cuboids with a chosen height and a base that spans one or more pixel areas. It is of particular note that this representation is especially useful in tomosynthesis and tomography applications, where the objective is to establish a 3D volume describing the breast consisting of a number of slices orientated parallel to the image receptor surface. In the case of projection mammography, we note that a special case may be adopted where parallelepiped voxels are used which span the entire thickness of any objects encountered by a ray (for example the breast or the compression plate) from point of entry of the photon beam to point of exit, and with a base spanning one or more pixel areas. Each voxel is assumed to consist entirely of a single type of tissue/material (that is, it is assumed homogenous), and therefore presents attenuation and scattering properties of a constant magnitude (partial volume effects are assumed negligible in this embodiment). Of course, other representations of the breast volume are also possible.

An efficient next cell calculation algorithm, depicted in Flow Chart 3 shown in Fig. 5, determines the set of voxels traversed between the points at which the primary ray enters, or the scattered ray originates, and the point of exit from the volume. This operates by stepping through the voxels, determining at each voxel if the next step is in the x, y or z direction by calculating which of the voxel walls is closest in terms of the direction of traversal, and terminating when the voxel containing the exit point is reached. This method is not only highly efficient in terms of computation; but, crucially, enables the direct formulation of the simultaneous equations that define the set of voxels that contribute to a particular pixel in the projection image (and thus is important for digital breast tomosynthesis

reconstruction).

X-Ray Breast Interaction Model: Absorption and Scattering

In the typical mammographic photon energy range consideration is given to attenuation due to photoelectric absorption, and scattering due to coherent and incoherent phenomena.

Calculating the Attenuation of a Primary Ray Traversing the Breast

Photoelectric absorption is the process by which an x-ray photon impinging on an atom transfers its entire energy to an inner shell electron, which is ejected from the atom as a result. Flow Chart 4 shown in Fig. 6 gives an overview of the process for calculating attenuation, taking into account photoelectric absorption. In the current embodiment of the invention, published data is used for the photoelectric cross-sections of the elemental constituents of any material encountered. In the molecular case, atomic values are combined using a summation weighted by the proportion of each present by mass (also gleaned from published data). A linear attenuation coefficient due to photoelectric absorption is calculated by combining the cross section data with published density data.

Calculating the Scatter Arising around a Primary Ray Traversal

It is known by those familiar with x-ray physics that both coherent (a change in direction but no energy transfer to the atom) and incoherent (a change in direction together with a portion of the energy being transferred to the atom) scattering effects occur within the mammographic energy range. A technique for the estimation of the contribution of scattered energy at each pixel within the image has been developed, and is an important aspect of the present invention. Its application enables the primary component to be established (that which has followed a straight traversal path) and thus the radiodensity measured at each pixel therefore describes solely the attenuation encountered by the primary beam.

Calculating Coherent Scatter

Coherent scattering is a result of the electric field component of the incident x-ray wave accelerating an electron, which, as it accelerates, emits radiation itself, and thus the incident wave is scattered. In the direction of travel of the incident wave, the emitted wave will be exactly out of phase, and therefore destructive interference will occur, thereby reducing its intensity. The coherent cross section for a single 'free' or isolated electron may be expressed mathematically. The atomic coherent scatter form factor, and the equivalent value in the molecular case, specifies the collective effect of interference in magnitude and phase between the scatterings of the wave front from the different electrons within an atom. It scales the scattered intensity, per electron in the system, from the free electron case to the actual value. The current embodiment of the invention utilises published data for the atomic form factor, and approximates the molecular form factor using the atomic values combined using a summation weighted by the abundance of each element present within the molecule by mass. We have discovered, via experimental investigations, that the discrepancy between the calculation of the final scatter signal at an image receptor pixel using a measured molecular form factor, and the approximated form factor, is negligible, though it is accepted that in certain specialised circumstances, for example when particularly high accuracy is required, other form factors may be employed.

Calculating Incoherent Scatter

Incoherent scattering may be thought of as an x-ray photon colliding with an outer shell electron of an atom and transferring part of its energy to it. Outer shell electrons are bound with very little energy compared to an x-ray photon so may be considered "free". The excited electron is termed a Compton electron and is ejected or moved into an excited atomic state. Due to the law of conservation of energy the photon energy is reduced. The differential Klein-Nishna (free-electron Compton) collision cross section per electron may be expressed mathematically, and in a technique akin to that in the coherent case, may be scaled by the incoherent scattering function which scales the scattered intensity from the free electron case to the actual value. The current embodiment of the invention relies upon a physical law that is well known to those versed in x-ray physics, utilises published data for the atomic incoherent scatter function, and approximates the molecular case using the atomic values combined using a summation weighted by the abundance of each element present within the molecule by mass.

Consideration of the Anti-Scatter Grid

The current embodiment of the invention includes consideration of scatter reduction devices, such as the anti-scatter grid. The likelihood of a photon successful traversing such a grid is derived from the geometry of the grid and the angle of incidence of the x-ray photon beam. Calculating the Scatter Incident on an Image Pixel Arising from a Primary Ray

In one embodiment of the invention, the scatter arising from each primary ray (that is the ray of photons following a straight traversal path linking the focal spot to an image pixel) is calculated according to the method depicted in Flow Chart 5 shown in Fig. 7. More specifically, the scattered energy contribution to the image signal at each of the image pixels within the immediate neighbourhood surrounding the pixel upon which the primary beam is incident is evaluated. The size and shape of the neighbourhood is governed by the desired accuracy of the resultant scatter estimation, together with the thickness of the compressed breast. For this reason, it is understood that many possible choices exist depending upon requirements. Each primary ray is approximated to consist of a number of short traversals, the scatter arising from the atoms within each being approximated as solely originating from an infinitesimally small point at the centre. This is illustrated in Fig. 8.

The scattered radiation is approximated as originating from the (in principle) infinitesimally small point P , and is a result of the x-ray photons traversing the material contained within distance dt . Using the molecular form factor and incoherent scatter functions together with the mathematical relations describing free electron scattering described earlier, the cross section (which when combined with density data gives a linear scattering coefficient) describing the scatter arising within the angular range required for it to be incident on each image pixel C is calculated. The attenuation and scattering occurring along the traversal path PC may then be modelled, though depending upon the accuracy required, it has been found that the effect of ignoring the secondary scattering (that is photons that are scattered for a second time or more) may be neglected, and therefore only consideration of photoelectric absorption effects is required. Through consideration of each scattering point P , and summing each of the contributions at each of the pixels, the scatter contribution for each primary ray is established. Summing the contributions from all primary rays gives the scatter signal present within the final image.

In one embodiment of the invention, the number of scatter origination points considered along the primary ray is sampled coarsely, and then the scatter arising from the unsampled intermediate points is found by interpolation using local fitting of suitable functions (for example second-order polynomials). The Nyquist-Shannon sampling theorem governs the minimum sampling frequency (since aliasing must be avoided), which in rum depends upon the variation in the scattering properties of the tissues encountered along the ray traversal.

Calculating the Scatter Image

Calculating the scatter arising around each primary ray to ascertain the scatter image is computationally intensive, and a major innovation of one embodiment of the invention is to utilise a number of high speed interpolations to streamline computation whilst maintaining sufficient accuracy to be of clinical use. The process is described graphically in Flow Chart 6 shown in Fig. 9.

Calculating the Scatter-to -Primary Ratio using a Homogenous Equivalent

Empirically, the ratio of the scattered radiation to the primary radiation at pixels located away from the edge of breast, that is those that are not significantly influenced by the boundary conditions, is approximately constant, regardless of the exact proportion of adipose to fibroglandular tissue making up the breast. One embodiment of the invention exploits this independence between tissue type and the scatter-to-primary ratio to significantly reduce the computational complexity of the scatter calculation. The scatter-to-primary ratio is calculated for each image pixel by approximating the breast based on the approximation that it consists entirely of a 50/50 adipose/fibroglandular homogenous (at an atomic level) mix of tissue. It is understood that other choices are available, and may be more optimal under certain conditions and requirements. The elemental composition and density of breast fibroglandular and adipose tissue is taken from the literature (the data reported by G. R. Hammerstein, D. W. Miller, D. R. White, M. E. Masterson, H. Q. Woodard, and J. S. Laughlin, "Absorbed radiation dose in mammography," Radiology, vol. 130, pp. 485-91, 1979 has been used, though other studies are equally applicable). The calculated scatter-to-primary ratios are then used to approximate the scatter component for the actual breast, which of course consists of an inhomogeneous mixture of tissues. The use of the homogenous breast approximation in the calculation of scatter allows many of the functions to be sampled at course intervals and then interpolation schemes utilised to ascertain the intermediate data. Since the interpolation schemes use simpler approximating functions, the computational complexity of the task is thereby reduced.

Rotational Invariance of the Scatter Arising around a Primary Ray

The amount of scatter around each primary ray in the case of the homogenous equivalent breast is approximately rotationally invariant: more specifically that is the iso-contours may be assumed circular and have offset centre points in the direction of the incident primary ray. A consequence of this is a low frequency, smoothly varying function describing the relation between scatter intensity and the angular coordinate at any given radius from the point at which the primary ray is incident. One embodiment of the invention samples this function to the extent required by the Nyquist-Shannon sampling theorem, and locally fits suitable functions (for example polynomials) between sample points, as illustrated in Flow Chart 7 shown in Fig. 10.

Spatial Interpolation of the Scatter Arising around Primary Rays

A reduction in computational complexity is obtained through sampling and interpolation to give the scatter "profiles" that arise for the primary rays incident at different spatial locations across the image receptor. An approximately two- dimensional quadratic variation across the area of the image receptor exists in scatter intensity between points at equivalent distances from the point of incidence of the associated primary ray. Exploitation of this allows the scatter "profiles" arising around each of the primary rays incident upon each pixel of the image receptor to be found with low computational execution time.

Calculating the Scatter-to-Primary Ratio Image at Low Resolution and Interpolating

Since the Fourier domain representation of the scatter image consists of only low frequency components, the Nyquist-Shannon sampling theorem enables the calculation of scatter at a considerably lower resolution than that required of the primary signal for clinical diagnosis. One embodiment of the invention calculates the scatter-to-primary ratio calculation at a low image resolution (selected subject to that allowed by the Nyquist-Shannon sampling theorem to prevent aliasing), and interpolates the result using cubic splines (though it is appreciated that other choices exist, and may be more optimal under certain conditions).

Calculating the Blurring due to Scatter at Image Discontinuities

The limitations of the approximation of a constant scatter-to-primary ratio occurs at tissue boundaries, where discontinuities in composition exist, and the scatter-to-primary ratio fluctuates significantly due to the difference in the magnitude of the scattered radiation arising from the varying tissue types within the vicinity. It is these fluctuations which result in image blurring, and thereby limit contrast and lead to inconspicuous image features. The use of the scatter-to-primary ratio calculated using the homogenous equivalent fails to include consideration of such fluctuations, and therefore the degrading effect of scatter upon the image. A further image filtering step is thus required, and is depicted graphically in Flow Chart 8 shown in Fig. 1 1.

The filter is derived by careful consideration of the nature of the scatter signal, in particular it is the additive effect of photons scattered from a large number of atoms distributed over a wide area. This results in the Fourier domain

representation of the scatter signal consisting solely of low frequency components. A low pass filter is applied to the scatter signal derived from the primary signal using the scatter-to-primary ratio calculated from the homogenous equivalent breast. In the spatial domain, the low pass filter takes the shape of the intensity profile of the scatter arising around the primary ray. A shift-invariant approximation is adopted to facilitate an optimal and therefore high speed filtering algorithm using the Fourier transform. The approximation is taken to be the shape of the scatter arising around the primary ray for which the angle of incidence upon the image receptor is perpendicular to the surface: it is accepted that other choices are possible, and may be preferable depending on circumstance and the accuracy required. The

approximation of the scatter signal may be expressed more succinctly by:

= p + s = p + r(p * k) = p + (rp)* k where is the acquired image, k is the normalised shape of the scatter arising around the primary ray for which the angle of incidence is perpendicular to the surface of the image receptor, p is the primary image, and r is scatter-to-primary ratio 'image' calculated using the homogenous equivalent breast. The application of the low pass filter should be limited to those areas of the image that are free from the effects the boundary of the breast has on the scatter signal. It is noted that discussion thus far has been concerned with the forward modelling of scatter, that is modelling of its participation in the image formation process. In the case of forward simulation the scatter signal may be found using the above relation, since the primary image is known. However, when making measurements of radiodensity from acquired images, the objective is to estimate the scatter signal from the acquired projection image, rather than from a precise knowledge of the underlying anatomy, which is unknown, and so too is the primary signal. In the reverse case, it is the acquired image, i , that is known, and a deconvolution step is therefore required, which is itself complicated by the need to apply it to the product of the scatter-to-primary ratio with the primary signal. In one embodiment of the invention the deconvolution is performed by approximating the shape (though not the actual intensity values) of the primary image by that of the acquired image, an approximation that holds well for small primary-to-scatter ratios. Under this approximation (rp)* k , the scatter image, is known, and may therefore be subtracted from the acquired image to give the scatter removed primary image.

Detector Model

In any practical realisation, a precise image detector model is determined by the particular detector technology being used. In order to model photon transport and signal production within a detector a detailed knowledge of the exact

constructional details is required, in the case of certain technologies for example the substrate thickness, material, and the exact details of any doping used are required, as well as gain factors for the electronics etc. In one embodiment of the invention, the detector characteristics are measured using an empirical procedure, rather than ascertained by modelling of the underlying physical phenomena occurring during signal production. The model of any given detector amounts to a specialisation of a technique, used in the present embodiment of measuring transfer functions which relate the x-ray photon energy absorbed from the flux of an incident beam by the 'active' (x-ray sensitive) component of the detector, to the resulting recorded pixel intensity. The polyenergetic nature of the photon beam, and thus beam hardening, is preferably given consideration because photon energy impacts upon detector response as it governs the likelihood of successful photon conversion to recordable image signal. Spectral detectors, such as the Medipix (developed by a collaboration based at CERN), record not only the overall photon flux, but the distribution of photon flux across photon energy. Such detectors have an associated energy resolution, and thus energies are recorded in "bins" sized accordingly, each bin corresponding to a range of photon energy. The model disclosed herein is applied to each of the bins separately, specifically, a separate transfer function is measured for each photon energy range for which a flux magnitude is recorded.

In the case of an amorphous selenium direct digital detector, to use a common technology as an example, the energy of the incident photon governs the number of electrons liberated from the substrate, as well as the energy the liberated electrons posses, which in turn governs the likelihood of the electron successfully traversing the substrate (together with the location of the point of production) and reaching the charge collector. The image signal is then directly proportional to the total charge, and hence the number of electrons collected. The invention utilises a set of transfer functions which relate the pixel intensity to sum of the x-ray photon energy absorbed by the "active" component of the detector, for a given beam quality (a photon energy distribution): i.e. the relationship between magnitude of incident photon flux (x-ray exposure) and pixel intensity. Each member describes a possible x-ray spectra which during clinical use may be incident upon the detector. Our extensive empirical studies show that the beam hardening fluctuations in detector response arising from variations in the incident x-ray spectra due to breast composition, can be ignored, when compared to those arising from the x-ray source beam quality configuration parameters (in the case of a vacuum x-ray tube, anode material, filtration material and thickness, and the anode-cathode accelerating potential difference); and the anatomical thickness. The current embodiment therefore calculates transfer functions describing each beam quality supported by the x-ray source in question, and for which each function is dependant on the thickness of the anatomy being imaged. In one embodiment of the invention, the detector transfer functions are sampled empirically for a set of beam qualities where said set is of a sufficient size and spacing so as to satisfy the Nyquist-Shannon sampling theorem and thus prevent aliasing. Also, for any sampled beam quality, the range of likely breast thicknesses encountered clinically are sampled empirically with a spacing sufficiently small to satisfy the Nyquist-Shannon sampling theorem and thus prevent aliasing.

Interpolation techniques using an appropriate set of basis functions, in one embodiment of the invention a second order polynomial, is used to ascertain the transfer functions applicable to beam qualities and anatomy thicknesses between sample points when required (illustrated schematically in Flow Chart 9 shown in Fig. 12).

Manufacturers of digital image detectors generally apply a "flat fielding" correction to the output images in order to ensure all pixels within the image have a uniform response. This is often achieved by regular imaging of a suitable flat fielding phantom, commonly a 40mm sheet of Polymethylmethacrylate (PMMA). The acquired image is used to calculate a multiplicative factor for each pixel so as to give a constant gain across all image pixels, and thus the image of the flat field phantom comprises of uniform pixel intensities across the spatial area of the detector. Such methods do not lead to a completely uniform field x-ray field being incident on the detector, due to such effects as the spectral variation arising from the anode heel effect; the variation in ray path length; the inverse square law; and the variation in traversal distance through the phantom. Consideration of the flat field correction is therefore required in the detector model in order to reflect the bias in the pixel gain correction introduced by the non-uniform x-ray field. In one embodiment of the invention, consideration is achieved by simulating the acquisition of the image of the flat fielding phantom, and calculating the multiplicative factors required to give a uniform pixel intensity, using the point at which the detector response was measured empirically as the reference point.

Measuring the Detector Transfer Functions

The detector transfer functions are sampled empirically using a specially developed technique which is independent of the particular detector technology. From these samples the coefficients of the function describing the response (which is usually linear in the case of digital detectors) can be fitted using well-known mathematical methods, for example regression. In order to experimentally sample the transfer function, an attenuator is required to suitably filter the x-ray beam in order that the spectrum incident upon the detector is typical of that exiting an average human breast. In the current embodiment PM A has been selected for this purpose, though it is understood that other materials may be substituted, and may be more suitable given particular requirements. A selection of attenuator thicknesses are utilised which cover the range of likely breast thicknesses encountered clinically. The incident photon beam upon the detector (in terms of both energy distribution and magnitude) is calculated using the x-ray source model disclosed within this invention, and then suitably attenuated for the photoelectric absorption of the thickness of the attenuator using the model of beam attenuation disclosed within this invention. In order to mitigate the effects of scatter, and therefore allow its exclusion from the calculation of radiation incident on the detector, a beam stop magnification technique is adopted. In this embodiment a suitable support is used to hold the attenuator as close to the x-ray source emission window as possible. Two small apertures in a highly attenuating material (of the order of 2mm in diameter and of a material such as lead) are placed above and below a given thickness of breast equivalent attenuator, aligned so as to facilitate the passing of a primary ray to a given position on the image detector. The considerable separation distance between the attenuator and the detector, together with the aligned apertures, results in a high degree of magnification, and thus only scattered radiation of a very low angle (thought to be insignificant in intensity) is included within the area of the aperture shadow upon the detector. It is therefore assumed that the pixel intensity measured within said aperture shadow is a result of only the primary radiation traversing the attenuator. Empirical readings of pixel intensity are made across a range of suitable exposures (the choice of which are dependant on the response characteristics of the exact detector in question) and the associated calculated incident x-ray energy are ascertained using the disclosed x-ray source model and assuming that all photons undergoing photoelectric absorption or scattering do not reach the detector (i.e. primary beam only), and thus the transfer function is sampled. A suitable set of images are acquired so as to sufficiently sample the range of clinically used beam qualities. In one embodiment of the invention second order polynomials are used to interpolate pixel intensity at a given exposure to provide pixel intensities for beam qualities between sample points.

In order to streamline the procedure for calibrating an x-ray mammography machine (see Flow Chart 10 shown in Fig. 13) a purpose-designed breast equivalent attenuator calibration object has been developed, which is design to sit on top of a telescopic magnification stand to hold it as close as possible to the x-ray source emission window. In one embodiment of the invention this consists of a number of PMMA thicknesses ranging between 30 and 100mm in steps of 5mm, though it is accepted that other choices of thickness and interval spacing are possible and maybe more optimal under certain conditions. Each block of PMMA of a given thickness is separated from the others by a vertical sheet of lead, although other choices of highly attenuating material may be appropriate. This structure is enclosed by a high attenuating material (for example lead) both top and bottom, and a suitably sized (for example 2mm diameter) and aligned pair of apertures are present in the top and bottom so that for each thickness of PMMA a narrow primary ray may pass through in such a way so that when in use to be incident upon a prescribed position on the image receptor. Figure 14 gives a schematic drawing describing a possible construction in which the shaded regions correspond to tissue-mimicking material, and the grey portions (horizontally and vertically) correspond to a highly x-ray attenuating material such as lead.

Applications of the Normalised Radiodensity Projection Images

Calculation of a 3D Radiodensity Volume from Multiple Projections Images

When a number of projection images acquired from a variety of orientations are available, as in the case of tomosynthesis, the set of equations describing the radiodensity at each voxel within the breast form a simultaneous set. The variety of acquisition orientations may be obtained through moving the x-ray focal spot, the image detector, or both simultaneously. Solution of the set of simultaneous equations yields a tomographically reconstructed three-dimensional volume, where each voxel describes the radiodensity encountered at the equivalent volume within the breast. This offers a novel possibility for use in breast tomosynthesis

reconstruction. Consider for example the acquisition geometry pictured in Fig. 15.

In this example, the normalised radiodensity projection image may be defined mathematically in terms of the voxels within the breast volume by:

∑∑∑T(i, j,k,x,y)x V(i,j,k) k j ; where R(x,y) denotes the normalised radiodensity projection image pixel with coordinates x, y ; V(i,j,k) is the radiodensity of the voxel i, j, k within the breast volume, and T(i, j,k, x, y) is the output of the ray tracing algorithm described previously: the distance traversed through voxel i, j, k by the primary ray incident upon pixel x, y . It is of note that T(i, j,k,x,y) is zero in the case that the ray does not pass through that breast volume element. In the case of digital breast tomosynthesis, the multiple projection images R_m(x, y) form a set of simultaneous equations describing the radiodensity at each of the voxels.

The calculation of R(x,y) may be defined mathematically as the solution of the following equation:

where P(x,y) is the observed pixel intensity in the raw projection image output by the mammographic acquisition equipment, Ό(ε) is transfer function (as a function of energy ε corresponding to the detector model, calculated according to the detector model disclosed in this invention; Α(ε) is a function describing the energy absorbed by the active component of the detector, also calculated according to the disclosed detector model; Ν(ε, ζ),α) is the number of photons of energy ε at beam quality Q and effective target angle a , calculated according to the x-ray source model disclosed in this invention; _ref (s) is the photoelectric absorption coefficient of the reference material at photon energy ε , gleaned from published data; and S(x,y) is the estimate of the scatter incident upon pixel x, y , calculated according to the model of scattered radiation disclosed in this invention. It is understood that, in addition, a function describing the transmission characteristics of an anti-scatter grid may be included within this expression whenever that is relevant. Note that the above expression relating pixel intensity P(x,y) at pixel (x, y) to the desired radiodensity R(x,y) s non-linear and so the solution is calculated using well-known numerical optimisation methods.

Experimental investigation suggests the signal-to-noise ratio in a digital mammography system is limited by the quantum noise which arises from counting a limited number of stochastic events. The noise in P(x, y) , and hence R(x,y), therefore depends upon pixel intensity, which itself is governed by incident energy. Testing using a olmogorov-Smirnov D test suggests that a Gaussian distribution provides a good fit to the noise observed in the empirical investigation (superior to that of a Poisson distribution), though it is understood that other models can be used instead. In the case of a Gaussian, as used in the current embodiment, the noise within a radiodensity image is described by the following equation:

where G„ is a stochastic Gaussian distribution noise source with 0 mean and standard deviation a(p(x,y))² which itself is a function of energy imparted to the detector. a(p(x, y)Y is established empirically for the appropriate beam quality and breast thickness.

The likelihood expression for use in digital breast tomosynthesis

reconstruction algorithms adopting the maximum likelihood approach is:

where R is the vector of normalised radiodensity projection images, and G is a Gaussian distribution with mean R_t{x,y) and standard deviation a(P,(x,y)). Again, it is understood that any other distribution can be substituted for the Gaussian and a corresponding result derived.

Numerous reconstruction algorithms exist in the literature which utilise the common aim of optimising the solution of the set of simultaneous equations by minimising the residual error between the acquired projection images, and those calculated by forward projection through the three dimensional data volume. One such example is the Simultaneous Algebraic Reconstruction Technique (SART), developed by Anderson et al (A. H. Andersen and A. C. Kak, "Simultaneous algebraic reconstruction technique (SART): a superior implementation of the art algorithm," Vltrason Imaging, vol. 6, pp. 81-94, 1984). It is appreciated that substitution of the maximum likelihood approach proposed here for any one of the algorithms in the literature which aims to optimise for a minimal residual error will yield the three dimensional radiodensity volume when the projection images are pre- processed using the radiodensity calculation algorithm disclosed herein, and that in certain cases, the use of one of the alternatives may be more optimal.

Traditional approaches to tomosynthesis have utilised a limited angular sampling range around the breast, often around ±30°. This results in a considerable null space in the simultaneous equations describing the voxels of which the breast volume comprises, and as a result resolution in the plane perpendicular to the breast compression plates is limited, with structures being "smeared out" between a number of reconstructed slices, rather than appearing sharply in only those slices which correspond to the locations in which they are physically present. The main limitation associated with increasing the sampling angle in tomosynthesis is that the thickness of breast tissue the x-ray beam is required to traverse rapidly increases with angle of acquisition, as illustrated in Fig. 16.

As the thickness of tissue increases, so too does the magnitude of the scattered radiation, and the absorbed radiation dose, which has the consequence of reducing the primary photon flux, and hence has a negative impact on the signal-to- noise ratio. Increased thicknesses of tissue are generally imaged using an x-ray beam consisting of a photon energy spectrum with a greater proportion of higher energy photons. Were this approach to be used in tomosynthesis in order to give

consideration to the larger tissue traversal distances between acquisition orientations, a systematic variation would be introduced in the appearance of the underlying anatomy between the resulting projection images. This would lead to errors in the output of the reconstruction algorithm that solves the system of simultaneous equations describing the breast voxels. The calculation of the normalised

radiodensity disclosed herein overcomes this issues, since each of the projection images are normalised so as to remove the dependency of the image pixel intensity on the beam quality employed for acquisition, and thus the resulting radiodensity values may be reconstructed to form a volume of voxels describing only the underlying tissues of the breast. Further, the removal of the effect of scattered radiation which results from the calculation of radiodensity, removes the impact of the increase of this effect resulting from increased tissue traversal thickness. The use of the radiodensity measures to normalise the projection images in tomosynthesis thus allows the variation of beam quality between projection images, and therefore the patient dose and the signal-to-noise ratio may be optimised for each projection individually, without impacting on the accuracy of the overall three dimensional reconstruction.

Quantifying the Attenuation of Tissue Features and Calcifications of the Breast for Diagnosis

One of the major advantages of using the disclosed radiodensity calculation algorithm to normalise x-ray mammograms is that the quantitative measure it provides may be used to characterise features of the breast, such as areas of tissue, the presence of which may form biomarkers indicative of malignancy. For example, the growth of a tumour often causes complex changes in the connective stroma in its surrounding, and thus an ability to characterise a lesions micro-environment is a strong tool to aid diagnosis. The disclosed radiodensity measure allows in- vivo characterisation and classification of tissue by facilitating a normalised basis upon which comparisons may be made between those values of radiodensity observed in an image of a breast, and those observed from tissue of known pathology: for example, tissue from a benign or malignant lesion, healthy tissue, the stroma surrounding an invasive carcinoma etc.

Microcalcifications, small deposits of a calcium compounds of sizes in the approximate range of 70 to 150 microns, may also be characterised using the disclosed radiodensity measure. These are a feature that is present in a significant proportion of mammograms, although they are associated with both malignancy, such as ductal carcinoma in-situ, and benign conditions, such as fat necrosis.

Histological studies suggest that microcalcifications may be divided into two types: those composed of Calcium Oxalate Dihydrate (also termed Wedellite), which are mainly found in benign ductal cysts and are rarely associated with malignancy; and those composed of Calcium phosphates, mainly Calcium Hydroxyapatite, which may occur in both benign and malignant lesions, but are seen most frequently in proliferative lesions including carcinomas. Currently computer aided detection algorithms aim to identify microcalcifications by searching for small regions of high pixel intensity. The quantitative radiodensity measure disclosed herein allows the attenuation of any such high intensity region to be reconciled against the radiodensity expected of Calcium Oxalate Dihydrate or Calcium Hydroxyapatite, and thereby provide a diagnostic classification of: likely benign; likely malignant; or an erroneous finding which is unlikely to be calcium. Figure 17 shows the variation in the radiodensity measure observed for the two Calcium compounds, and inspection reveals two significantly different relations relating the diameter of the calcification to the calculated radiodensity. In one embodiment of the invention, the calcification diagnosis algorithm takes the form of subtracting the radiodensity observed within the candidate calcification area from the radiodensity observed at the pixels in the immediate vicinity, to yield a "background compensated" radiodensity value, which is then used together with the approximate diameter measured from the area of the candidate calcification to ascertain which of the two relations shown in Fig.l 7 the observed data better fits: the closer of the two being that which is returned as the diagnosis.

The small size of microcalcifications (70 to 150 microns) compared to the resolution of the mammographic image detector (often 100 microns) results in the images of calcifications, particular the smaller ones more often associated with malignancy, being under-sampled according to the Nyquist-Shannon sampling theorem. As may be seen from inspection of Fig. 17, this impacts on the observed variation of radiodensity with calcification diameter, and therefore one embodiment of the invention utilises the "low resolution" relations plotted, where the resolution in question is matched to that of the resolution of the mammographic image detector used to acquire the image so as to mitigate this issue.

This approach to calcification diagnosis is a specific case of the more generalised technique of model predictive diagnosis proposed here. This technique uses both feedback and feed-forward inputs to a model of image formation in order to predict the appearance of a given type of lesion (or calcification, cyst or other anatomical feature; for convenience the term 'lesion' alone is used in the following) within a patient under examination, in terms of the radiodensity measure proposed herein. Through comparison between that predicted, and that observed, diagnosis is made on the basis of their similarity.

The feedback inputs to the model come from segmenting the suspected lesion in the acquired image, in order to quantify its size and identify its shape, together with any significant features within it, and the quantification of the radiodensity exhibited by the tissues in its surroundings. The term 'segmenting' means identifying the boundaries of an image feature, e.g., cyst, lesion, calcification, or a sub-feature within such a feature, e.g. a highly dense region within a lesion such as an invasive carcinoma.

The feed-forward inputs to the model are prior knowledge of the expected radiodensity of the lesion, which itself is dependent upon its elemental composition and density. These inputs are based on a hypothesis regarding the suspected lesion. The forward model is used to predict the appearance of the lesion in the tissue surroundings observed for the compressed breast thickness of the patient, and the presence or absence of the lesion is decided by a metric quantifying the similarity in radiodensity between the lesion and the surrounding tissue in the acquired patient image, with that predicted by the model. The locations of candidate regions are identified by searching for maxima in a correlation metric quantifying the similarity between the pixels corresponding to the anatomical feature, lesion, or calcification in the simulated mammographic image, with those in the acquired mammographic image compensated for the effects of the attenuation of the background tissue in the localised surroundings.

For example, in one embodiment of the invention, a cyst, a fluid-filled mass having a distinct membrane separating it from the surrounding tissue, be it a type 1 case containing fluid with an electrolyte composition similar to that of intracellular fluid (sodium-to-potassium ratio <3), or type 2 in which the fluid's content resembles that of plasma (Na/K ratio >3), is modelled using the attenuation of water (since both fluids are composed for the most part of water) in the shape of an ellipsoid, for which the major and minor diameters are measured from the acquired patient

mammographic image. The contrast in the radiodensity measure across the suspected cyst boundary is measured in both the acquired image and the predicted image computed using the image formation model, and the percentage difference between the two used as a metric to quantify the similarity between the two, from which a diagnosis may be made if the metric satisfies a given threshold or criterion.

This technique is clearly advantageous when applied to conventional projection mammography at least. Consider a 20 mm diameter cyst in a 60 mm thick breast; it will have a different radiodensity from that if it were in a 75 mm thick breast. However, with the aid of the model, one can predict the radiodensity of the cyst feature in both the 60 mm and 75 mm thicknesses, and hence inform the diagnosis.

Similarly, in another embodiment of the invention, ductal carcinoma in-situ are taken (hypothesised) to have a cylindrical shape for which the diameter may be segmented from the acquired image, and the likely radiodensity taken (hypothesised) to be that measured from population statistics within a collection of specimen radiographs of excised histology and pathology samples of various grades of in-situ carcinoma. These parameters of the hypothesised carcinoma are applied to a model to predict its appearance in a simulated mammogram. A high degree of similarity between prediction and observation not only provides diagnosis, but may also be used to grade the case (depending on the grade of the most closely matching radiodensity), thereby contributing to the selection decision of the most appropriate course of treatment. The invention may also be applied to invasive carcinomas, in one embodiment the central foci (which may be several in number in multi-focal cases), the surrounding mass, and the invasion of the surrounding tissue (for example, the alteration of the connective stroma) are segmented and their radiodensity quantified relative to the surrounding healthy tissue, and diagnosis is made on the basis of the similarity between that observed in the patient image and that predicted by the modelled conditioned using a-priori radiodensity measurements from specimen radiographs of known pathologies.

Testing

Extensive testing of the image formation model and normalisation technique has been performed using homogenous and inhomogeneous tissue equivalent phantoms (i.e. test objects) of our own design fabricated using breast tissue resins mimicking the x-ray characteristics of both adipose and fibroglandular tissue, manufactured by CIRS (Computerised Imaging Reference Systems, Inc, Norfolk, VA). A tissue equivalent phantom has also been specifically developed to investigate scattering properties, and to test the performance of our algorithmic correction against that of a conventional anti-scatter grid. Significant validation has been conducted through the comparison of simulated images using the model of image formation and those acquired using a GE Medical Systems (Milwaukee, WI) Senographe 2000D, and comparison of simulated and normalised images of our known composition phantom test objects acquired using a IMS Giotto (Bologna, Italy) Image MD, and a Hologic Lorad (Bedford, MA) Selenia. Agreement is good, with worst case accuracy bounds to have been found of around 10% on the entire model.

Figure 18 shows photographs of the test object manufactured from tissue equivalent resins (white = breast fibroglandular; grey (actually terracotta coloured) = adipose) to experimentally verify the normalisation performance of the model, and an x-ray apparatus used to test the model.

The graph of Fig. 19 shows the pixel intensity of the raw images for a selection of step sizes in the above test object. For each step size, the pixel intensity is shown for the full range of clinical tube potentials and the associated exposures for the Molybdenum target and filter combination, acquired on a Hologic Lorad

(Bedford, MA) Selenia.

Following normalisation of the data of Fig. 19 by the software

implementation of the invention disclosed herein, the graph of Fig. 20 shows the resulting radiodensity measures corresponding to the pixel intensities at the selected steps on the phantom (test object). Inspection reveals the excellent performance of the technique's normalisation.

Fig. 21 illustrates schematically a computer system embodying the invention for performing a method described above according to the invention. The software for performing the method is stored in data store 50 and executed by processor 52. Such software constitutes a computer program and may be written in any suitable language. Data corresponding to the digitised images to be analysed can also be stored in data store 50. Images and results can be displayed by the processor on display 54. The software and image data may be stored and/or distributed on any suitable computer-readable medium, such as magnetic and/or optical disc-shaped media, or solid-state computer memory, and may be fixed to the computer system or may be removable. Typically the data store 50 will be a hard drive of the computer system. The software and image data may also be distributed and/or accessed as signals sent over a communications network, such as the internet. The computer system does not necessarily have to be directly connected to the x-ray imaging system, and can perform the analysis according to a method of the invention off-line, using previously acquired image data. The computer system would receive input values regarding the image acquisition geometry, and the x-ray source and detector parameters. Other input output devices (not shown) such as keyboard, mouse, printer, connected to the processor 52, enable the user to control the operation of the computer system and output results.

Industrial Applicability As noted above, in current clinical practice image analysis is qualitative (relies upon human judgement). However, the need for quantitative analysis is increasingly recognised by clinicians, researchers, developers of computer aided diagnosis software, and others. In particular, the estimate of breast density, the study of such variations and their correlation with the occurrence of malignant lesions is expanding rapidly. Mainly due to a lack of suitable normalisation algorithms, this work is generally conducted using qualitative visual assessment, and little to no consideration is given to the effect on image appearance of the beam quality used to acquire the image, the x-ray exposure, or the scatter arising. Applications of quantitative measurements are many: from classification of tissue to use as an image quality metric.

Currently the vast majority of mammography systems utilise an anti-scatter grid consisting of cellulose interspaced focused lead strips. Though these selectively reject photons according to angle incidence, they are not 100% efficient in either passing primary photons, or rejecting scattered ones. As a result the exposure, and the hence the patient dose, has to be increased by the 'Bucky factor' generally a factor of approximately 2, though sometimes more. The software scatter correction capabilities of the invention allows for the removal of the contrast limiting blur introduced to the image by scattered radiation, and hence the scatter grid is no longer required for this purpose. The removal of the scatter grid has the highly desirable effect of reducing patient dose, which since the vast majority of mammograms are taken of health women for screening purposes must be reduced to as little as possible, particularly in light of recent announcements of the extended age range to be included with the UK breast screening programme.

The vast majority of tomographic reconstruction algorithms approximate the polyenergetic x-ray beam as monoenergetic, and thus introduce beam hardening artefacts. The contribution of scattered radiation to the signal measured at each image pixel is also generally ignored. Both theses effects are significant, and result in visible artefacts in reconstructed images. The technique proposed here overcomes these two limitations, since the radiodensity measurements have a low dependence on beam quality, and are compensated for scatter. This not only improves the quality of the reconstruction by better describing reality, it allows direct comparison of intra and inter patient images.

The normalisation technique described may be applied to any x-ray projection image for which quantified measurements are required.

Claims

Claims

A method to calculate for each pixel of an x-ray projection image of a mammal's breast a normalised quantitative measure of the x-ray attenuation characteristics of the tissue encountered along the traversal path between the point of entry of the ray of x-ray photons into the mammal's breast closest to the x-ray source, and the point of exit from the mammal's breast closest to the image detector, where said method comprises the combination of the following:

• obtaining the x-ray photon flux emitted from the x-ray source incident upon the mammal's breast;

• calculating the primary x-ray photon flux exiting the mammal's breast and being subsequently recorded by the image detector;

• calculating a single scaling factor quantifying the difference between said incident and exit primary fluxes, per unit length of the traversal path through the breast, relative to the attenuation characteristics of a reference material.

The method recited in claim 1 wherein said calculation of the primary x- ray photon flux exiting the mammal's breast and being subsequently recorded by the image detector consists of the combination of the following:

• calculating the photon flux detected at each image pixel using a

transfer function model describing the response of the image detector in terms of the input photon flux absorbed by the detector given the resultant output pixel intensity from the mammographic image;

• calculating the component of the photon flux recorded by the image detector arising from scattering phenomena occurring within the breast; • subtracting the calculated component arising from scattered radiation from the photon flux recorded by the image detector so as to yield the primary flux exiting the breast.

A method to calculate for each pixel within a spectral x-ray projection image of a mammal's breast a normalised quantitative measure of the x- ray attenuation characteristics of the tissue encountered along the traversal path between the point of entry of the ray of x-ray photons into the mammal's breast closest to the x-ray source, and the point of exit from the mammal's breast closest to the image detector, where said method comprises the combination of the following:

• obtaining the x-ray photon flux emitted from the x-ray source incident upon the mammal's breast;

• calculating the primary x-ray photon flux exiting the mammal's breast and being subsequently recorded by the image detector for each of a plurality of separate photon energy ranges for which the spectral detector quantifies the incident flux;

• calculating a plurality of weighting factors which when applied to the attenuation characteristics of the corresponding plurality of reference materials quantify the difference between said incident and exit primary fluxes across the plurality of discrete photon energy ranges for which the spectral detector quantifies the incident flux, per unit length of the traversal path through the breast calculated.

The method recited in claim 3 wherein said calculation of the primary x- ray photon flux exiting the mammal's breast and being subsequently recorded by the image detector comprises the combination of the following:

• calculating the photon flux detected at each image pixel, for each discrete photon energy range for which the spectral detector quantifies the incident flux, using a plurality of transfer function models for each recorded photon energy range, describing the response of the image detector in terms of the input photon flux in the given photon energy range absorbed by the detector given the resultant output pixel intensity from the mammographic image;

• calculating the component of the photon flux recorded by the image detector, for each discrete photon energy range for which the spectral detector quantifies the incident flux, arising from scattering phenomena occurring within the breast;

• subtracting the calculated component arising from scattered radiation from the photon flux recorded by the image detector to obtain the primary flux exiting the breast for each discrete photon energy range for which the spectral detector quantifies the incident flux.

A method comprising combining a plurality of measurements in order to calculate the photon spectrum describing the photon count in each of a plurality of energy ranges emitted from an x-ray tube, for a given direction of beam emission, said measurements comprising at least two of the following:

• a calculation of the effective target angle for the direction of emission of the beam under consideration, for the geometry of the tube of the physical x-ray unit, which yields an identical self-filtration to that present for an equal target angle in the tube geometry used for spectral measurement;

• a calculation of the spectrum emitted from the target surface for the target angle in the measured spectral data which is equal to the calculated effective target angle, and where linear interpolation of the logarithm of the measured photon count and the reciprocal of the tangent of the target angle is used to yield values between the plurality of target angles samples for which spectra have been measured;

• a calculation of the effect of the tube accelerating potential waveform using a time of occurrence weighted sum of a number of spectra generated from constant tube potentials covering the range of variations present within the accelerating potential waveform;

• a calculation of the attenuation of the beam of x-ray photons emitted from the target surface according to the length of the traversal path taken through the tube window and any added beam filtration that is present;

• a calculation of a scaling factor defining the effect of the inverse

square law between the focal spot and the image detector pixel;

• a calculation of an empirical correction factor consisting of a scaling factor which when applied to the modelled x-ray beam emitted from the tube window results in said beam having a photon flux rate which is equal to that measured experimentally on the physical x-ray unit;

• a calculation of an empirical correction factor consisting of the

addition of extra tube filtration of a thickness which when applied to the modelled x-ray beam emitted from the tube window results in said photon beam having a first half value layer equal to that measured experimentally on the physical x-ray unit.

A method for calculating the component of photon flux recorded by an image detector arising from scattering phenomena occurring within a mammal's breast during mammographic x-ray image acquisition, wherein said method comprises:

• calculating the ratio of scatter-to-primary radiation at each image pixel whilst approximating the mammal's breast by assuming the tissues of which it comprises are of a single homogenous tissue type only, where said calculation is performed at an image resolution at or above the lowest resolution allowed by the Nyquist-Shannon sampling theorem;

• interpolating the resolution of said image of the ratio of scatter-to- primary radiation to match that of the acquired mammographic image; • calculating an approximation of the scattered radiation at each image pixel in the acquired x-ray mammographic image using said calculation of the scatter-to-primary ratio and the incident photon flux corresponding to the pixel intensity recorded by the image detector in the acquired image; and

• applying a low-pass image filtering operation to the approximated scatter image.

The method recited in claim 6 wherein said calculation of the ratio of scatter-to-primary radiation within the approximated homogenous breast at each image pixel comprises:

• calculating the image resulting from the primary photon fluxes

recorded by the image detector of the beams of primary x-ray photons incident on each pixel, using a model of the x-ray source to ascertain the x-ray photon flux incident on the approximated breast, and a ray tracer to identify the traversal path through the breast and thus the resulting attenuation of the x-ray photons by photoelectric absorption and scattering phenomena to yield the primary photon flux incident upon the image detector;

• calculating the image resulting from summing the scattered photon fluxes recorded by the image detector arising from the scattering phenomena occurring along the traversal paths taken by each of the primary rays through the approximated breast, where the scattered flux contributed by each traversal is calculated for the plurality of image pixels in the neighbourhood surrounding the image pixel upon which the primary beam associated with said traversal is incident; and

• calculating the ratio of the primary flux to that of the scattered flux.

The method recited in claim 7 wherein said calculation of the scattered photon fluxes recorded by the image detector arising from the scattering phenomena occurring along the traversal paths taken by each of the primary beams through the approximated breast, comprises: • calculating for a sample subset of the primary ray traversal paths, the scattered flux incident on a plurality of image pixels in the

neighbourhood surrounding the image pixel upon which the primary ray is incident, where the size of said subset is selected to be of a sufficient size, and the members of said subset are situated at suitable locations within the image, to satisfy the Nyquist-Shannon sampling theorem; and

• interpolating between said samples to obtain the scattered radiation arising from the primary ray traversal path to all image pixels.

The method recited in claim 8 wherein said calculation of the scattered flux arising from a primary ray traversal path and being subsequently incident on the plurality of image pixels in the neighbourhood

surrounding the image pixel upon which the primary ray is incident, comprises:

• dividing the primary ray traversal path into a plurality of sub- traversals, and assuming that all scattered radiation arising along the length of each respective sub-traversal originates from a single respective point, and where taking all the sub-traversals together encompasses the full traversal path; and

• calculating for a sample subset of said sub-traversals, where the subset is of a sufficient size, and the members of the subset are situated at suitable locations within the full traversal, so as to satisfy the

Nyquist-Shannon sampling theorem, the scatter flux incident upon a one dimensional sample line of image pixels within the plurality of image pixels in the neighbourhood surrounding the image pixel upon which the primary ray in question is incident, such that a full rotation of the one dimensional sample line encompasses all pixels in the neighbourhood;

the method further comprising at least one of: • using rotational invariance in order to interpolate said scatter fluxes across the two dimensional area of the plurality of pixels on the image detector; and

• interpolating between the sample sub-traversal s using second order polynomials to yield the scattered flux from all the sub-traversals from which the full traversal path comprises.

The method recited in claim 9 wherein said application of a low-pass image filtering operation to the approximated scatter image comprises filtering said image using a normalised kernel of a substantially identical shape to that of the scattered radiation recorded by the image detector arising from the traversal path through the breast taken by a beam of primary x-ray photons, where said primary beam is of such a size that it is incident upon a single image detector pixel only.

An apparatus to enable the collection of data for characterising the pixel intensity output response of a mammographic image detector to the primary beam x-ray photon flux incident upon it at a dose equivalent to that exiting a breast during mammography, where said apparatus comprises:

• a breast equivalent attenuation object consisting of a plurality of portions of different thicknesses manufactured of a material with x- ray attenuation properties substantially equivalent to those of a mammal's breast tissue; where each portion is encased within a highly x-ray attenuating material; and where said highly attenuating material contains two small apertures, one situated above, and one situated below each portion of tissue equivalent material, aligned in such a way as to allow x-ray photons emanating from the x-ray focal spot to pass through said tissue equivalent material and subsequently expose the image detector;

• a support stand to hold the breast equivalent attenuation object up to the x-ray source emission window whilst sitting on the breast support table of a mammographic image acquisition device. A method to calculate the photon flux incident upon a pixel of a mammographic image detector from the recorded pixel intensity, for the range of x-ray photon energies which the detector considers when recording said pixel intensity, where the method comprises the following:

• acquiring calibration images of the apparatus defined in claim 1 1 such that a representative sample encompassing the range employed in clinical practice is gleaned of the incident x-ray photon fluxes arising from different anatomical thicknesses, and beam qualities;

• measuring the pixel intensities associated within the image of the shadow of each aperture of the breast equivalent attenuation object in each of the empirical calibration images;

• calculating the photon flux incident on the image detector within the image of the shadow of each aperture of the breast equivalent attenuation object in each of the images using a model of x-ray emission from the x-ray source, attenuated using the Beer-Lambert law to remove the photons photoelectrically absorbed, or scattered, within the breast equivalent attenuation object;

• interpolating between the data points for the sampled beam qualities using regression fitted second order polynomials to ascertain data points for the beam quality used to acquire the recorded image in question;

• interpolating between the data points for the sampled anatomical thicknesses using regression fitted second order polynomials to ascertain data points for the anatomical thickness present in the recorded image in question; and

• regression fitting a linear relation to the calculated interpolated data points to glean the transfer function relating pixel intensity to incident photon flux, and using said function to calculate the incident photon flux associated with the recorded pixel intensity in question. A method to calculate a three dimensional data volume describing the constituent tissue of a mammal's breast, where said volume consists of a plurality of voxels each describing the radiodensity of the equivalent volume of tissue within the breast; the method comprising:

• acquiring a plurality of x-ray projection images from a variety of focal spot positions and detector orientations;

• calculating the radiodensity image corresponding to each of the

acquired projection images according to the method recited in claim 1 or 2;

• calculating a three dimensional data volume using a method by which the radiodensity is used as a common basis between projections so as to facilitate reconstruction by the minimisation of the residual error between the radiodensity observed in each of the acquired projections images, and those calculated by forward projection through the three dimensional data volume.

A method to characterise tissue in-vivo which comprises comparing the plurality of weighting factors relating the observed attenuation

characteristics to those of the plurality of reference materials, calculated according to the method recited in claim 3 or 4, to those measured from tissue for which the underlying pathology is known.

A method to identify whether a small highly attenuating area of a mammographic image of a mammal's breast corresponds to tissue calcification, and if so to classify such area into one of a plurality of classes, wherein said method comprises a combination of the following:

• calculating the normalised radiodensity of the projection image

according to the method recited in claim 1 or 2;

• calculating the difference in radiodensity between the candidate

calcification area and the surrounding tissue;

• comparing the radiodensity values of the likely calcium compounds present, of a thickness equal to an approximate measure of the diameter of the candidate calcification area, where the radiodensity value is corrected for any under sampling due to the resolution of the mammographic image detector that acquired the x-ray image, and returning the result of which calcium compound matches the radiodensity most closely, or a negative result if no good match is found.

A computer-implemented method of identifying a feature within a mammographic image of a mammal's breast which corresponds to the physical presence of a specific anatomical feature, lesion, or calcification within the breast; the method comprising:

• inputting a hypothesis regarding the type of the anatomical feature, lesion or calcification so as to predict parameters characterising said anatomical feature, lesion or calcification;

• applying a model of image formation to calculate a prediction of the appearance of said anatomical feature, lesion or calcification in a simulated mammographic image using the parameters predicted in the previous step;

• locating and segmenting candidate regions in the mammographic image which potentially correspond to said specific anatomical feature, lesion or calcification;

• comparing the appearance of the predicted anatomical feature, lesion or calcification in the simulated mammographic image, to that observed in the acquired mammographic image, compensated for the effects of the differences in the attenuation of the background tissue in the localised surroundings, to obtain a metric of the comparison for identifying the presence of the anatomical feature, lesion, or calcification.

A method according to claim 16, wherein the parameters characterizing the anatomical feature, lesion, or calcification, comprise at least one of: x- ray attenuation characteristics, shape, size, and/or the shape and size of any internal features.

18. A method according to claim 16 or 17, wherein the hypothesis is

informed by segmenting in the acquired mammographic image the boundaries and/or internal features of the anatomical feature, lesion, or calcification, to provide the characterisation parameters.

19. A method according to claim 16 or 17 or 18 wherein the location of a candidate region, corresponding to the anatomical feature, lesion, or calcification, is identified by searching for maxima in a correlation metric quantifying the similarity between the pixels corresponding to the anatomical feature, lesion, or calcification in the simulated

mammographic image, with those in the acquired mammographic image compensated for the effects of the attenuation of the background tissue in the localised surroundings.

20. A method according to claim 16, 17, 18 or 19, wherein the pixels of the mammographic image comprise x-ray attenuation characteristics calculated according to the method of claim 1 or 2.

21. A method according to any one of claims 16 to 20, wherein the metric of the comparison comprises the percentage difference between the contrast in x-ray attenuation characteristic of the group of pixels corresponding to the suspected anatomical feature, lesion or calcification measured across the boundary of said anatomical feature, lesion or calcification, calculated for the acquired image and that simulated using the model of image formation.

22. A method according to any one of claims 16 to 21 , wherein the applying and comparing steps are repeated for each of a plurality of different specific types of anatomical feature, lesion, or calcification.

23. A computer program comprising computer-executable code that, when executed on a computer system, causes the computer system to perform a method according to any preceding method claim. A computer-readable medium storing a computer program according to claim 23.

A computer program product comprising a signal comprising a computer program according to claim 23.