Quantifying image quality in
diagnostic radiology using
simulation of the imaging
system and model observers
Gustaf Ullman Radiation Physics, Department of Medicine and Health Faculty of Health Sciences Linköping University, Sweden Linköping 2008©Gustaf Ullman, 2008
Cover picture/illustration: An oil painting by Gustaf Ullman representing a chest radiograph Published articles and figures have been reprinted with the permission of the copyright holder. Printed in Sweden by LiU‐Tryck, Linköping, Sweden, 2008 ISBN 978‐91‐7393‐952‐2 ISSN 0345‐0082
Don’t worry about saving these songs! And if one of our instruments breaks, it doesn’t matter We have fallen into the place where everything is Music. The strumming and the flute notes rise into the atmosphere, and even if the whole world’s harp should burn up, there would still be hidden instruments playing. So the candle flickers and goes out. We have a piece of flint and a spark. This singing art is sea foam. The graceful movements come from a pearl somewhere on the ocean floor. Poems reach up like spindrift and the edge of driftwood along the beach, wanting! They derive from a slow and powerful root that we can’t see Stop the words now. Open the window in the center of your chest, and let the spirits fly in and out. Jalal al‐Din Rumi
1. INTRODUCTION ... 1 1.1. Radiation protection in diagnostic radiology ... 1 1.2. Optimisation of diagnostic radiology ... 2 1.3. Optimisation using a Monte Carlo based computational model ... 2 2. OBJECTIVE ... 5 3. MONTE CARLO BASED COMPUTATIONAL MODEL OF THE IMAGING SYSTEM... 7 3.1. Introduction... 7 3.2. Computational model of the x‐ray imaging systems ... 9 3.2.1. Model of the imaging system ... 9 3.2.2. Monte Carlo simulation of photon transport... 14 3.2.3. Scoring quantities... 18 3.2.4. Calculated quantities ... 19 3.3. Calculation of images from the high‐resolution phantom ... 20 3.4. Uncertainties... 22 3.4.1. Stochastic uncertainties ... 22 3.4.2. Systematic uncertainties... 22 4. ASSESSMENT OF IMAGE QUALITY ... 25 4.1. Introduction... 25 4.2. Image quality assessment as developed in this work... 26 4.2.1. The task... 26 4.2.2. Model of the imaging system and patient... 27 4.2.3. Observers... 29 4.2.4. Figures of merit ... 30 5. RESULTS AND DISCUSSION ... 41 5.1. Ideal observer with a simplified patient‐model ... 41
5.2. Low resolution voxel phantom ... 43 5.3. High resolution voxel phantom... 44 5.4. Ideal observer with simple anatomical background... 46 5.5. Correlation to human observers ... 49 5.6. Model observers with complex anatomical background ... 52 6. SUMMARY AND CONCLUSIONS... 59 7. FUTURE WORK ... 61 8. ACKNOWLEDGEMENTS ... 63 9. REFERENCES... 65
ABSTRACT
Accurate measures of both clinical image quality and patient radiation risk are needed for successful optimisation of medical imaging with ionising radiation. Optimisation in diagnostic radiology means finding the image acquisition technique that maximises the perceived information content and minimises the radiation risk or keeps it at a reasonably low level. The assessment of image quality depends on the diagnostic task and may in addition to system and quantum noise also be hampered by overlying projected anatomy.
The main objective of this thesis is to develop methods for assessment of image quality in simulations of projection radiography. In this thesis, image quality is quantified by modelling the whole x‐ray imaging system including the x‐ray tube, patient, anti‐scatter device, image detector and the observer. This is accomplished by using Monte Carlo (MC) simulation methods that allow simultaneous estimates of measures of image quality and patient dose. Measures of image quality include the signal‐to‐noise‐ratio, SNR, of pathologic lesions and radiation risk is estimated by using organ doses to calculate the effective dose. Based on high‐resolution anthropomorphic phantoms, synthetic radiographs were calculated and used for assessing image quality with model‐observers (Laguerre‐Gauss (LG) Hotelling observer) that mimic real, human observers. Breast and particularly chest imaging were selected as study cases as these are particularly challenging for the radiologists.
In chest imaging the optimal tube voltage in detecting lung lesions was investigated in terms of their SNR and the contrast of the lesions relative to the ribs. It was found that the choice of tube voltage depends on whether SNR of the lesion or the interfering projected anatomy (i.e. the ribs) is most important for detection. The Laguerre‐Gauss (LG) Hotelling observer is influenced by the projected anatomical background and includes this into its figure‐of‐merit, SNRhot,LG. The LG‐observer was found to be a better model of the radiologist than the ideal observer that only includes the quantum noise in its analysis. The measures of image quality derived from our model are found to correlate relatively well with the radiologist’s assessment of image quality. Therefore MC simulations can be a valuable and an efficient tool in the search for dose‐ efficient imaging systems and image acquisition schemes.
LIST OF PAPERS
This thesis is based on the following papers
I. Gustaf Ullman, Michael Sandborg, David R Dance, Martin Yaffe, Gudrun Alm Carlsson. A search for optimal x‐ray spectra in iodine contrast media mammography. Phys. Med. Biol. 50, 3143–3152 (2005)* II. Gustaf Ullman, Michael Sandborg, David R Dance, Roger Hunt, and
Gudrun Alm Carlsson. Distributions of scatter to primary ratios and signal to noise ratios per pixel in digital chest imaging. Radiat Prot Dosim, 114, no 1‐3, 355‐358 (2005)*
III. Gustaf Ullman, Michael Sandborg, David R Dance, Roger A Hunt and Gudrun Alm Carlsson. Towards optimization in digital chest radiography using Monte Carlo modelling. Phys Med Biol 51, 2729‐ 2743 (2006)*
IV. Michael Sandborg, Anders Tingberg, Gustaf Ullman, David R Dance and Gudrun Alm Carlsson. Comparison of clinical and physical measures of image quality in chest and pelvis computed radiography at different tube voltages. Med. Phys. 33(11) 4169‐4175 (2006)*
V. Gustaf Ullman, Alexandr Malusek, Michael Sandborg, David R. Dance and Gudrun Alm Carlsson. Calculation of images from an anthropomorphic chest phantom using Monte Carlo methods. Proc of SPIE 6142, (2006)*
VI. Gustaf Ullman, Magnus Båth, Gudrun Alm Carlsson, David R Dance, Markku Tapiovaara, and Michael Sandborg. Development of a Monte Carlo based model for optimization using the Laguerre‐Gauss Hotelling observer. (To be submitted to Med Phys)
Other peer reviewed papers by the author not included in the thesis
1. Gustaf Ullman, Michael Sandborg, David R Dance, Roger Hunt, and Gudrun Alm Carlsson. The influence of patient thickness, tube voltage and image detector on patient dose and detail signal to noise ratio in digital chest imaging. Radiat Prot Dosim, 114, no 1‐3, 294‐297, 2005 2. Markus Håkansson, Magnus Båth, Sara Börjesson, Susanne Kheddache,
Gustaf Ullman, Lars Gunnar Månsson. Nodule detection in digital chest radiography: effect of nodule location. Radiat Prot Dosim 114, no 1‐3, 92‐96, 2005
3. R A Hunt, D R Dance, P R Bakic, A D A Maidment, M Sandborg, G Ullman and G Alm Carlsson. Calculation of the properties of digital mammograms using a computer simulation. Radiat Prot Dosim 114, no 1‐3, 395‐398, 2005
4. D R Dance, R A Hunt, P R Bakic, A D A Maidment, M Sandborg, G Ullman and G Alm Carlsson. Breast dosimetry using a high‐resolution voxel phantom. Radiat Prot Dosim 114, no 1‐3, 359‐363, 2005
5. Roger A Hunt, David R Dance, Marc Pachoud, Gudrun Alm Carlsson, Michael Sandborg, Gustaf Ullman and Francis R Verdun. Monte Carlo simulation of a mammographic test phantom. Radiat Prot Dosim, 114, no 1‐3, 432‐435, 2005.
Conference presentations
1. Ullman G, Sandborg M, Dance D R, Skarpathiotakis M, Yaffe MJ, Alm Carlsson G. (2002) A search for optimal x‐ray energy spectra in digital iodine subtraction mammography using Monte Carlo simulation of the imaging chain. Digital Mammography IWDM 2002: Proceedings of the Workshop, Bremen, Germany, June 2002. Ed. Peitgen H‐O (Springer‐ Verlag, Berlin) pp152‐154, 2002
2. M. Båth, M. Håkansson, S. Börjesson, S. Kheddache, C. Hoeschen, O. Tischenko, F. O. Bochud, F. R. Verdun, G. Ullman, L. G. Månsson. Investigation of components affecting the detection of lung nodules in digital chest radiography. Accepted for presentation at Medical Imaging, 12‐17 February 2005, San Diego, USA. Proc. SPIE 5749, 231‐242, 2005.
Internal reports (not reviewed) 1. Gustaf Ullman, Michael Sandborg, Roger Hunt and David R Dance. Implementation of simulation of pathologies in chest and breast imaging Report no 94, ISRN ULI‐RAD‐R‐‐94—SE, 2003
2. Gustaf Ullman, Michael Sandborg and Gudrun Alm Carlsson. Validation of a voxel‐phantom based Monte Carlo model and calibration of digital systems. Report no 95, ISRN ULI‐RAD‐R‐‐95—SE, 2003
3. Gustaf Ullman, M Sandborg, D R Dance, R Hunt and G Alm Carlsson Optimisation of chest radiology by computer modelling of image quality measures and patient effective dose Report no 97, ISRN ULI‐ RAD‐R‐‐97—SE, 2004
4. Gustaf Ullman, M Sandborg, Anders Tingberg, D R Dance, Roger Hunt and G Alm Carlsson Comparison of clinical and physical measures of image quality in chest PA and pelvis AP views at varying tube voltages Report no 98, ISRN ULI‐RAD‐R‐‐98—SE, 2004
5. Gustaf Ullman, M Sandborg, D R Dance, M Båth, M Håkansson, S Börjesson, R Hunt and G Alm Carlsson On the extent of quantum noise limitation in digital chest radiography Report no 99, ISRN ULI‐RAD‐R‐‐ 99—SE, 2004
6. Gustaf Ullman, Michael Sandborg, David R Dance, Roger Hunt and Gudrun Alm Carlsson Distributions of scatter‐to‐primary and signal‐to‐ noise ratios per pixel in digital chest imaging Report no 100, ISRN ULI‐ RAD‐R‐‐100—SE, 2004
ABBREVIATIONS
AGD Average glandular dose ALARA As low as reasonable achievable APR Apical pulmonary region AUC Area under the ROC curve BKE Background known exactly BV Background varying C Contrast CC Cranio‐caudal C/CB Nodule‐to‐bone contrast Crel Relative contrast DQE Detective quantum efficiency E Effective dose FN False negative FOM Figure of merit FP False positive Ht Equivalent dose HIL Hilar region Kc, air Collision air kerma LG Laguerre‐Gauss LAT Lateral pulmonary region LME Lower mediastinal region LNT Linear non‐threshold hypothesis MC Monte Carlo MTF Modulation transfer function NPS Noise power spectrum PA Posterior Anterior RET Retrocardial region ROC Receiver operating characteristics SKE Signal known exactly SNR Signal‐to‐noise ratio SNRhot, LG Laguerre‐Gauss Hotelling observer signal‐to‐noise ratio SNRI Ideal observer signal‐to‐noise ratio SNRp Signal‐to‐noise ratio per pixel TN True negativeTP True positive UME Upper mediastinal region VGA Visual grading analysis VGAS VGA score Energy imparted per unit area from primary photons p A ε s A ε p Energy imparted per unit area from scattered photons Mean energy imparted per primary photon λ Mean squared energy imparted per primary photon 2 p λ s Mean energy imparted per scattered photon λ Mean squared energy imparted per scattered photon 2 s λ
1. INTRODUCTION
1.1. Radiation protection in diagnostic radiology
Diagnostic x‐ray examinations can support the radiologist with valuable information that can be utilised to give a patient an accurate diagnosis, and subsequently a successful treatment. However, imaging with ionising radiation is also associated with a small risk for cancer induction or genetic detriment. When x‐ray photons are scattered or absorbed inside the cells of the human body, ionisations occur that can alter molecular structures and thus make harm to the cell. The most important damage to the cell is damage in the DNA since this may induce mutations. Ultimately, the damage may lead to that the cell is killed, and if enough cells are killed, the function of the tissue or organ will be deteriorated. This type of acute harm due to large radiation exposures is referred to as a deterministic effect. However, at the relatively low radiation exposures in diagnostic radiology, the damages caused by ionising radiation are often rather easily repaired. Yet, sometimes the damage on the DNA is more complex. This can cause mutations or chromosome aberrations, which in turn may lead to a modified cell but with retained reproduction capacity. In some cases, such modified cells can result in a cancer. In the case where the harmful effects of ionising radiation are only known statistically, it is referred to as a stochastic effect. The risk related to stochastic effects to a human from exposure from ionising radiation is often quantified with the effective dose, E (ICRP 1991, ICRP 2007).
According to the linear non‐threshold (LNT) hypothesis, there is a linear relation between the effective dose and risk for cancer induction (ICRP 2005) and means that the collective dose can be used as a measure of the harm to the population. The collective dose from medical radiography is according to the Swedish radiation protection authority (Andersson et al 2007) 8000 manSv per year or 0.9 mSv on average per capita, and contributes the largest fraction of the total dose to the population from man‐made sources. Diagnostic radiology is invaluable for the health care but due to the radiation risks, radiation protection of the patient becomes an important issue. Three different principles are used for radiation protection (ICRP 2007). The first principle is justification. Ionising radiation should only be used in those situations where it brings more good than harm. The second principle is
optimisation. It means that, in those cases where the use of ionising radiation is
justified, doses should be kept as low as reasonable achievable. This is often referred to as the ALARA (As Low As Reasonably Achievable) principle. The third principle is dose limits to the individual. However, this principle is more applicable for personnel rather than for patients in diagnostic radiology. 1.2. Optimisation of diagnostic radiology Optimisation means to balance the diagnostic information (image quality) and patient dose so as to maximize the ratio between the two; either to keep the information constant and minimize the dose or to increase information at constant dose. The dose to the patient undergoing an x‐ray examination has, in digital systems, a close relation to the quantum noise in the image. The quantum noise depends on the number of photons incident on the image detector and is approximately described with a compound poisson distribution, which takes the energy absorption properties of the detector into account. If we use too few photons, the image will be noisy and it will make it difficult or even impossible for the radiologist to give a correct diagnosis. It may also take longer time for the radiologist to give a diagnosis using a noisy image. Yet, above a certain dose level, the quantum noise may become negligible in comparison to the noise naturally present in the projected anatomy (Hoeschen et al 2005). There will therefore be limited benefit to increase the dose above this level.
How to make the trade off between the dose to the patient and the image quality is a complex subject. A key aspect for the optimisation of diagnostic radiology is to understand the relative importance of the quantum noise in the image and the structures in the projected anatomy that act as noise. Several authors including Kundel et al (1985), Samei et al (1999), Burgess et al (2001) and Håkansson et al (2005b) have acknowledged the importance of projected anatomy in relation to quantum noise. The consensus from these studies is that at normal exposures, the projected anatomy is the most important factor in hampering the detection of subtle nodules in chest radiographs and mammograms.
1.3. Optimisation using a Monte Carlo based computational model
One method that has been utilised to search in a systematic way for the optimal imaging parameters in diagnostic radiology is to use a model of the imaging system, including the patient and observer, and to simulate the photon transport through the imaging system using the Monte Carlo method.
With this method it is possible to simultaneously calculate the dose to the patient and measures of image quality.
However, the physical measures of image quality derived from simulations must in some sense give us information on the usefulness of the image for a radiologist to solve a specific clinical task. Our physical measures of image quality must therefore correlate to clinical measures of image quality. Two methods for assessment of clinical image quality are given attention in this work, receiver operating characteristics (ROC) (Metz 1986) and visual grading analysis (VGA) (Tingberg 2000). A challenge in this work has been to develop a model, which includes realistic measures of image quality that takes the projected anatomy into account.
2. OBJECTIVE
While patient doses are relatively straightforward to calculate, image quality assessment is a more complex task and crucial for the optimisation process. The main objective of this thesis is therefore to further develop methods for assessment of image quality in x‐ray projection radiography. The main method is Monte Carlo photon transport simulation (Monte Carlo model) through the whole x‐ray imaging system including a model of the image observer. As study cases, chest posterior‐anterior (PA) and mammography cranio‐caudal (CC) projections are used as these are particularly challenging for the radiologist. The specific objectives are: • To study how physical measures influencing image quality are distributed over the image plane (paper II) • To develop methods for calculating physical image quality measures from simulated radiographs and search for correlations between these measures and measures of clinical image quality (papers III and IV) • To develop patient models of higher realism and finer anatomical structures for calculation of synthetic x‐ray images to be used for image quality analysis (papers V and VI)
• To complete our model of the imaging system by including a more realistic model observer that can be used to directly make any task‐ related clinical image quality assessment from synthetic images calculated by the model (papers V and VI) • To use our model of the imaging system towards optimisation of image quality and patient dose (paper I and III)
3. MONTE CARLO BASED
COMPUTATIONAL MODEL OF THE
IMAGING SYSTEM
3.1. Introduction
The Monte Carlo method relies on taking random samples from known distributions and is particularly useful for studying complex problems with many degrees of freedom. One of the first applications of the method was in Los Alamos, USA, during the Second World War where it was used to simulate neutron diffusion. Today, Monte Carlo methods are employed in widely diverse fields, from the evaluation of shares on the stock market (Glasserman 2003) to the calculation of energy levels of molecules with quantum Monte Carlo (Ceperley and Alder 1986).
In radiation physics, the Monte Carlo method is employed for simulating radiation transport, mathematically described by the Bolzmann equation. There are several general‐purpose computer codes available for the study of radiation transport, for example, MCNP (Monte Carlo N‐Particle transport) (Briesmeister 2000) developed in Los Alamos and designed to transport neutrons, electrons and photons; EGSnrc (Electron Gamma Shower) (Kawrakow and Rogers 2003, Nelson et al 1985), initially developed in Stanford, which transports photons and electrons; PENELOPE (PENetration and Energy Loss Of Positrons and Electrons) (Baro et al 1995) developed at University of Barcelona, and used to transports electrons, positrons and photons.
In diagnostic radiology, one of the most common applications of the Monte Carlo method is in patient dosimetry. There are several Monte Carlo computer codes that are used to estimate the effective dose. Jones and Wall (1985) used the Monte Carlo method to compute organ doses using a mathematical representation (Cristy 1980) of a human anatomy. Zankl and Wittman (2001) have developed a family of more realistic, segmented anthropomorphic voxel phantoms for organ dosimetry for external photon beams. Zankl and Petoussi‐ Henss (2002) calculated conversion factors based on the VIP man (Spitzer and Whitlock 1998) anthropomorphic model. The user‐friendly Monte Carlo computer program PCXMC by Servomaa and Tapiovaara (1998) calculates
organ and effective doses based on either measured air kerma‐area product or entrance air kerma values.
There are also Monte Carlo codes developed for optimisation in diagnostic radiology. Such codes rely on the fact that they are able to estimate both organ or effective doses and measures of image quality. The main application of the Monte Carlo method is for estimating the negative effect of scattered photons reaching the image detector. Chan and Doi (1985) used the Monte Carlo method to characterise scattered radiation in x‐ray imaging. Chan et al (1985) also investigated the performance of anti‐scatter grids in screen‐film imaging whereas Sandborg et al (1994a) did task‐dependent, anti‐scatter grid optimisation for digital imaging. More recently McVey et al (2003) did an optimisation study of lumbar spine radiography and Lazos et al (2003) have developed a software package for mammography. The Lazos model also includes a realistic model of the breast (Bliznakova et al 2003). Son et al (2006) have developed software that calculates images from the visual human (VIP man)(Xu et al 2000). They have used the EGSnrc code as a basis of the model, used model observers and calculated effective dose.
In this work we have used an in‐house Monte Carlo code VOXMAN adapted for conditions usually encountered in diagnostic radiology. It originates from Dance and Day (1984) and Persliden (1986) who independently developed computer programs to estimate scattered radiation in the image plane in mammography and conventional radiography, respectively. A few years later, Dance et al (1992) and Sandborg et al (1994b) merged the codes and did further validation of the computer programs. McVey et al (2003) replaced the simple homogeneous water or tissue phantoms, used in the earlier versions of the code, by a voxelised anthropomorphic male phantom developed by Zubal et al (1994). This step enabled more realistic organ dosimetry and made it possible to describe how measures of physical image quality vary in the image plane behind the patient.
The main focus of this thesis is on chest imaging. Therefore we have mainly used the VOXMAN model, adapted to simulate chest radiography. In paper I we used the version of the computer program dedicated for mammography. This computer program was further developed by Hunt et al (2005) to incorporate an anthropomorphic model of the breast developed by Bakic et al (2002). A brief description of the VOXMAN model is presented below.
3.2. Computational model of the x‐ray imaging systems
The Monte Carlo based computational method used in this thesis models the x‐ray imaging system and simulates photon transport from the source through patient, anti‐scatter grid and into the image detector. The computational method consists of the following components:
A model of the imaging system. This comprises different sources of input data and the imaging geometry. Input data includes x‐ray spectrum, patient‐based voxel phantom, anti‐scatter grid, table or chest support couch and image detector.
Monte Carlo simulation of photon transport through the imaging system. The model uses different variance reduction techniques, briefly described below, to increase the efficiency of the model.
Scoring variables such as organ and effective doses and calculation of different measures of image quality such as contrast and signal‐to‐noise ratio of nodule lesions or anatomical structures within the patient model.
3.2.1. Model of the imaging system
The input data files are described below including geometry, x‐ray spectra, voxel phantom, anti‐scatter grid and image detector.
3.2.1.1 Imaging geometry
The imaging geometries for the chest and mammography models are shown in figures 3.1 and 3.2, respectively. Examples of specific imaging configurations are listed in table 3.1. Substantial variations of the imaging system configuration were employed particularly in papers I, III and IV and to some extent also in papers II, V and VI.
Figure 3.1. The simulated imaging geometry used in chest PA radiography including an x‐ray source, voxel phantom, anti‐scatter grid and image detector.
Figure 3.2. The simulated imaging geometry used in Cranio‐caudal (CC) mammography. Notations are a) focus‐detector distance; thickness of b) breast, c) compression plate, d) adipose layer, e) contrasting detail, f) breast support, g) anti‐scatter grid and h) image detector.
Table 3.1. Examples of imaging system configurations for chest and breast
imaging system component Chest PA Breast CC imaging used in this work. Typical values Focus‐detector distance (cm) 180 65 Tube voltage (kV) 90‐150 ) Cu Cu Patient PA or breast thickness (cm) 0‐28 es 10‐25 ocalcifications Compression plate ‐ exiglas erial arbon fibre / Al / grid ratio (mg/cm2) 20‐55 Total filtration (mm 0.1‐0.5 mm 0.3 mm 25μm Rh 2 2‐8 Typical diagnostic tasks and size of details Nodul mm Micr and soft tissue masses 3 mm pl Grid interspace mat C Carbon fibre Lamella frequency (cm‐1) 40 / 12 60 / 5
Image detector material BaFCl, CsI CsI
Image detector thickness 100 100
.2.1.2 X‐ray spectra
was calculated with a computer program based on a
t
the VOXMAN model, the relative fractions of Bremsstrahlung and
3
The x‐ray spectrum
spectral model by Birch and Marshall (1979). The program calculates Bremsstrahlung and characteristic x‐rays from a tungs en or molybdenum anode target and allows the user to select appropriate thicknesses of added filtration of aluminum, copper or molybdenum. In paper I, tungsten, molybdenum and rhodium anode target spectra were instead calculated with MCNP4C Monte Carlo code since the Birch and Marshall program did not include a rhodium target or a rhodium filter.
In
characteristic x‐rays were computed and a random number selected from which of the distributions the photon emerged. If a Bremsstrahlung photon was selected, the photon energy was chosen using rejection sampling (Sandborg et al 1994b).
3.2.1.3
was
Low‐resolution chest phantom
anthropomorphic voxel phantoms were Voxel phantoms
tom a A Mammography phan
In the mammography model, simple representation of the female breast used (Ullman et al 2005). The breast was assumed to be a cylinder with semicircular cross‐section and made of a homogeneous mixture of glandular and adipose tissue in the central region surrounded by an adipose layer. The tissue compositions were taken from Hammerstein et al (1979). The density of glandular tissue was 1.04 g cm−3 and for adipose tissue 0.93 g cm−3. The
glandularity of the central part of the breast model was for the main part set to 50%, but was allowed to vary between 10‐90% to represent both dense and fatty breasts.
B
In the chest model, three different
used as a model of the patient. The main phantom was the one developed by Zubal et al (1994) and used in papers II, III and IV. The Zubal phantom (displayed in figure 3.3) was segmented into organs such as lungs, heart and bone marrow. It therefore allows for calculation of organ and effective doses. The female‐specific organs: breast, ovaries and uterus were added manually to the male body to enable effective dose to be calculated (McVey et al 2003). However, the phantom has relatively large voxels (3x3x4 mm3) and is
therefore not suitable for calculating realistic images (see figure 3.4 below). In addition, the lungs are comparably small since the phantom was based on a CT scan where the male patient was imaged with non‐inflated lungs and in a non‐upright position, contrary to the typical chest PA imaging situation. igure 3.3. Volume‐rendered representation of the Zubal phantom (left) and F the Kyoto Kaguku (PBU‐X‐21) phantom (right). An outline of the lungs, trachea, heart and breast are shown in the left image.
C High‐resolution chest phantoms
s (voxel size
Two high‐resolution voxel phantom 0.97x0.97x0.6 mm) were
he manufacturer of the Kyoto Kagaku phantom claims that it is composed of created from CT scans of two different anthropomorphic thorax phantoms: the Alderson phantom and the Kyoto Kagaku PBU‐X‐21 phantom. The Alderson phantom was used in paper V and did not include smaller vessels. The more recent Kyoto Kagaku phantom was more realistic since it contained a more human‐like distribution of small and medium‐sized vessels. Simulated x‐ray images of the Zubal, Alderson and Kyoto Kagaku phantoms are shown in figure 3.4 demonstrating an increasing realism from left to right.
T
materials with linear attenuation coefficients (μ values) resembling those of human tissues. However, they failed to provide us with details of the atomic compositions of the tissue substitute materials, which complicate a more rigorous comparison with real x‐ray images of their chest phantom (see paper VI). A more detailed description of the segmentation of the Kyoto Kagaku phantom is given in Malusek (2008). The Alderson and Kyoto Kagaku phantoms are used mainly for simulation of synthetic images with high resolution but are not yet segmented into organs and tissue types and can therefore not be used for direct calculation of effective dose. Figure 3.4. Projection images of the three chest voxel‐phantoms used in this .2.1.4 Anti‐scatter grid
as simulated by specifying the lamella thickness, thesis. To the left the low‐resolution Zubal phantom, central the Alderson
phantom and to the right the Kyoto Kagaku phantom.
3
The anti‐scatter grid w
interspace material and thickness as well as cover thickness and grid ratio. Typical grids are listed in table 3.1. In the Monte Carlo program the focused grid was simulated by an analytical transmission formula developed by Day
and Dance (1983). Scattered photons generated in the grid itself was simulated by a separate Monte Carlo simulation in a parallel grid (Sandborg et al 1994b). 3.2.1.5 Image detector ge detector was simulated in a separate Monte Carlo he image detector thickness was specified in terms of a surface density in mg .2.2. Monte Carlo simulation of photon transport tion is described by
(
n n n n)
n = r ,E ,Ω ,w The response of the imamodel of a semi‐infinite layer of the image detector material. This model is included as a sub‐routine to the main VOXMAN program. Energy imparted per unit area was assumed proportional to the image detector signal. The image detector model does not include the transport of secondary electrons as the kerma approximation was assumed. The transport of light photons was also neglected. In papers V and VI, the detector response was calculated separately with MCNP4C and was used for the calculation of primary projections (Malusek 2008).
T cm-2
(see table 3.1). In paper I, the detector material was needle crystals of CsI; in paper II, III and IV an unstructured mixture of BaFCl grains simulating a computed radiography (CR) fluorescent screen and finally in paper VI, a Gd2O2S indirect flat panel (DR) fluorescent screen was employed. 3 The physical state of the photon after the n:th interac α (3.1)
here rn is the position, En is the energy, Ωn is the solid angle and wn is the
.2.2.1 Photon interaction, cross‐sections and material compositions w
statistical weight of the photon. Photon interaction types in the energy range of diagnostic radiology (10‐150 keV) are: coherent scattering, incoherent scattering and photoelectric effect. The photon interactions are described by the differential cross sections for these events based on the atomic composition of the materials and tissue types in the geometry. A flow chart of the main steps in the Monte Carlo program is given in figure 3.6. A central part of the Monte Carlo method is the utilisation of a random number generator. In this work we have used the random number generators embedded in UNIX or LINUX operating systems.
3
) , ( ) cos 1 ( 2 2 2 2 Z x F r d dσcoh = e + θ Ω (3.2)
where re is the classical electron radius, x is defined by sin(θ/2)
hc E x= where h is Planck’s constant and c the speed of light. F is the atomic form factor, θ the scattering angle and Z the atomic number.
For incoherent scattering the differential cross‐section is given by the Klein‐ Nishina relation times the incoherent scattering function S(x,Z): ) , ( sin 2 2 2 2 Z x S E E E E E E r d d incoh e ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ′− + ′ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ′ = Ω θ σ (3.3)
Here, E is the incident photon energy and E´ is the scattered photon energy given by the Compton relation ) cos 1 ( 1+κ − θ = ′ E E 2 /m c E e = κ (3.4) where , me is the electron rest mass. For the photoelectric effect, it is assumed that the photon is locally absorbed in interactions with atoms of low atomic number (such as carbon and oxygen). Elements with high atomic numbers such as those in the grid (lead) and detector materials (e.g. barium, gadolinium, cesium, iodine) may emit (high‐ energy) characteristic x‐rays if vacancies are created in the K‐ or L‐shells.
To describe these photon interactions we have used tabulated cross‐sections from the XCOM library by Berger and Hubbell (1987). Cross‐sections for compounds were computed based on the relative weight of individual elements. The atomic form factors, F(x,Z), were given by Hubbell and Överbö (1979) and the incoherent scattering functions, S(x,Z), were given by Hubbell et
al (1975).
In the voxel phantom model of the patient, each organ is identified with one of four tissue types, with different densities: average soft tissue (1.03 g cm−3),
healthy lung (0.26 g cm−3), cortical bone (1.49 g cm−3) and bone spongiosa (1.18
g cm−3). The tissue densities and compositions were taken from ICRU 46 (1992)
3.2.2.2 Primary photons Primary transmission was calculated by first sampling the initial energy from the pre‐calculated energy spectrum, and use Siddon’s algorithm (Siddon 1985) to calculate the radiological path‐length from the focus to the detector n N n n d L
∑
= = 1 μ (3.5)The radiological path‐length can be used to calculate the contribution to the energy imparted per unit area from primary photons as dE E f e E r E s E L p A ( , ) ) ( 2 , ξ ξ ε Ω −
∫
= E sΩ, ( , ) (3.6)where is the source intensity,f E ξ is the detector absorption efficiency function depending on the photon energy E and cosine of incidence angle,
θ cos = ) , ( , at the detector surface; r is the focus detector distance. ξ Equation 3.6 is calculated by two separate methods. In the original version of the VOXMAN program, this calculation was embedded inside the Monte Carlo code. It is then calculated by first sampling the photon energy E, from the x‐ray spectrum, subsequently the optical path L from the source to a point in the detector is calculated. Finally the photon is transported in a semi‐infinite layer corresponding to the detector thickness in order to calculate f E ξ
) , (
. However, this did not allow for the calculation of high‐resolution images, since we were forced to run the Monte Carlo simulation for all these points. In papers I, III and IV, the energy imparted to the image detector per unit area from primary photons was therefore calculated to a very limited number of points (1‐15) in the detector plane corresponding to those points where the projection of the contrasting detail or lesion was located. In paper II, V and VI the Monte Carlo simulations where performed for 40 x 40 points in the detector. This is rather time‐consuming and is only feasible to perform with high precision with a fast, modern computer.
In papers V and VI, a different method was implemented for the calculation of primary projections. The detector absorption efficiency function f Eξ was calculated separately with MCNP4C (Malusek 2008) and the projections were calculated analytically averaged over the energy spectrum. This allowed for a separate calculation of primary projections with high resolution.
3.2.2.3 Scattered photons and variance reduction techniques
The simulation of scattered photons is time consuming. Therefore, different variance reduction techniques, described briefly below, were used to increase the efficiency. Monte Carlo methods that do not employ any variance reducing techniques are often referred to as analogue Monte Carlo methods. An algorithm for sampling the free path of the scattered photon, referred to as the Coleman’s algorithm is also briefly described below.
A Coleman’s algorithm
The free path of the scattered photon is sampled using an algorithm described by Coleman (1968). The sampling of the free path consists of several steps. First, the distance to the first interaction point is sampled for a homogeneous medium with the linear attenuation coefficient, μmax, corresponding to the
material with highest attenuation (e.g. bone). The sampling is performed by testing whether a sampled random number from a uniform distribution in the interval [0,1] is less than the quotient μ/μmax, where μ is the attenuation
coefficient of the material at the interaction point. If yes then the new point is accepted and the algorithm ends. If no then the sampling of the distance to the first interaction in the homogenous medium is repeated until the sampled random number is less than μ/μmax.
B Collision density estimator
Analogue Monte Carlo methods are inefficient in estimating scattered photons in the image plane due to the low probability that a scattered photon will pass a given small target area in the image detector. Therefore in the VOXMAN code, the collision density estimator (Persliden and Alm Carlsson 1986) is used. The contribution to the energy imparted per unit area at a given point of interest in the image detector is obtained from each interaction point in the phantom. The contribution εs∗ is derived through s n N n n n s w T , 1 ) (α λ ε
∑
= ∗= (3.7)where λn,s is the contribution from the n:th interaction and T(α) is the
probability for the photon of state αn to be scattered to the point of interest; wn
is the photon weight. In the collision density estimator, incoherent and coherent scattering are treated separately. The radiological path‐length from the interaction point to the point of interest in the detector is calculated with Siddon’s algorithm as in the case of primary photons.
C Analytical averaging of survival and Russian roulette The main purpose of the Monte Carlo model is to achieve an accurate estimate of scattered photons generated in the patient and emerging towards the image detector. If photons are absorbed in the patient they will not contribute to this estimate. Therefore, a technique known as analytical averaging of survival is used which does not allow photons to interact by the photoelectric effect. All interactions in the phantom are therefore constrained to be either coherent or incoherent scatterings. The new statistical weight, wn+1, for the photon after the n:th interaction is calculated from the cross‐sections for photoelectric, τ(E) and
scattering processes, σ(E) to correct for the bias which this method would otherwise introduce.
( )
( ) ( )
E E E w wn n σ τ σ (3.8) + = +1 For high photon energies, E, the ratio wn+1/wn is less than, but close to 1 and the statistical weight is only slightly reduced at each interaction. However, as the photon energy is reduced the relative importance of photoelectric cross‐ sections increases, and the number of interactions before a scattered photon escapes from the phantom geometry may be large. Hence, the statistical weight of the photon may eventually be low and so its contribution to the estimated image detector signal. Therefore, an unbiased procedure called Russian roulette (Salvat et al 2003) is used. Once the weight is less than 0.05 a random number is selected and in 95% of the cases the photon history is terminated; in the other 5% of the cases the photon history continues with a twenty times higher (100%/5%) statistical weight wn+1 compared to the original weight. Photon histories are also terminated once the photon is scattered out of the boundaries of the phantom. 3.2.3. Scoring quantities 3.2.3.1 Energy imparted to the image detectorEstimates of the mean energy imparted per unit surface area of the image detector from scattered photons,
ε
Asand primary photonsp
A
ε
are computed. The total energy imparted is calculated as the sum of primary and scatter contributions Ap Ast
A ε ε
ε = +
p
. In order to estimate the signal‐to‐noise ratio and variance in the image detector signal, the first and second moments of energy imparted per incident primary (
λ
andλ
2p) and scattered (λ
sand ) photon at the image detector was calculated (Dick and Motz 1981, Sandborg and Alm Carlsson 1992).2
s
λ
3.2.4. Calculated quantities 3.2.4.1. Contrast
A contrasting detail (e.g. corresponding to a lesion) is added to the model with a thickness and location specified by the user. The contrasting detail is not added directly into the voxel phantom but in an artificial way in a subroutine inside the VOXMAN model. The contrast of this detail is calculated as 1 1 2 1 1 1 p s p p p C ε ε ε ε ε + ⋅ − = (3.9)
where εp1 is the mean energy imparted to the detector per unit area from
primary photons with the nodule present, εp2 the mean energy imparted per
unit area from primary photons with the nodule absent and εs is the mean
energy imparted per unit area from scattered photons. 3.2.4.2. Signal‐to‐noise ratios The program calculates two types of signal‐to‐noise ratios. The signal‐to‐noise ratio per pixel, SNRp is calculated as 2 2 s s p p p A p N N SNR λ λ ε ⋅ + ⋅ = (3.10)
where N is the number of photons incident on a pixel. The indices p and s stands for contributions from primary and scattered photons, respectively. The quantities λ and λ2 are mean and mean squared values of the energy imparted
to a specified pixel per incident photon.
Given the location and thickness of a specified lesion (detail), the computer program calculates the signal‐to‐noise ratio for this detail with a projection area corresponding to one pixel. It is here called the SNRMC and is given by 2 2 1 1 2 2 1 1 s s p p p p p p MC N N N N SNR λ λ λ λ ⋅ + ⋅ ⋅ − ⋅ = (3.11) where the index n=1 refers to a pixel in the image behind the nodule, and n=2 refers to the same pixel with the nodule absent.
3.2.4.1 Air collision kerma
The air collision kerma, Kc,air is given by
( )
(
( )
)
∫
⋅ ⋅ = E Φ E μ E /ρ dE Kc,air E en air air en/ ) ( , (3.12)where μ ρ is the mass energy absorption coefficient for air and is the differential photon fluence with respect to energy. ) ( E E Φ t tH w E=
∑
3.2.4.2 Effective doseThe effective dose is the tissue‐weighted sum of the equivalent doses in all specified tissues or organs of the body calculated according to ICRP 60 (ICRP 1991).
(3.13)
where wt is the tissue or organ weighting factor and Ht the equivalent dose for
that tissue or organ. It is recognized that a new ICRP publication 103 (ICRP 2007) has recently been adopted and uses slightly different values of the tissue weighting factors in the calculation of effective dose. A detailed analysis of the effect on the figures of merit due to this change, for example signal‐to‐noise ratio per effective dose, SNR2/E, has not been performed here. The absolute values of SNR2/E may change, but the main conclusions on for example the appropriate tube voltage for chest PA radiography is unlikely to be affected by the change of weighting factors, particularly since the weighting factors for the lungs are the same in ICRP 60 as in ICRP 103.
3.3. Calculation of images from the high‐resolution phantom
The scatter projection, SNRp and other quantities calculated with the MC method are rescaled (i.e. from 40 x 40 points) to fit the number of pixels in the primary image. In paper V this number is 1760 x 1760, and in paper VI the number of pixels is 2688 x 2688. The rescaling is performed using a bilinear interpolation function in MATLAB. The interpolated scatter projections are added to the primary image to give the estimate of the mean energy imparted per unit area of the detector for the i:th pixel, s,
p , t ,i Ai Ai A ε ε ε = + . Noise is subsequently added to the image. In paper V, white gaussian noise is added with the standard deviation
i p, p , SNR i A i ε σ = .
The white noise is generated by adding sampled random numbers for each pixel i from the appropriate distribution to the calculated image. In paper VI, correlated noise is added with a method similar to the one used in Båth et al (2005c). In order to add correlated noise, knowledge of the noise power spectrum (NPS) for the clinical system is needed. The NPS is then normalized to correspond to unit variance. A random phase is added to the square root of the normalized NPS with the constraint that the random phase image φ( vu, ) should have the symmetry φ(u,v)= −φ(− −u, v). By taking an inverse Fourier transform of this spectrum a real and correlated noise image is created. The vector δ is a multivariate random variable with mean corresponding to a null vector and covariance matrix corresponding to the measured NPS. The noise fluctuation for each pixel is rescaled with the relation
εˆ δ εˆ i A i i A, σ ˆ , δε = ⋅δε . It is assumed that the NPS was invariant under a logarithmic transformation. The total energy imparted to the detector per unit area including noise fluctuations becomes Ai t i A i A , , t , ε δε ε = + t ,i A ε c b a gi = ln( Ai)+ t , ε . The pixel value in the i:th pixel is calculated by taking a logarithmic transformation of using the relation (3.14)
where the parameters a, b and c are calculated with non‐linear regression to make the best fit to the real phantom images. The method to add primary, scattered and noise images is illustrated in figure 3.5.
+
=
+
Figure 3.5. The method for calculating images illustrated in a cutout in the retrocardial region (see further figure 4.1.). From the left: primary projection, scatter projection, noise image and total image (to the right).
3.4. Uncertainties
3.4.1. Stochastic uncertainties
The choice of the number of photon histories used in the simulation is a trade off between computer time and statistical precision. If the statistical precision of the simulation is doubled, the computer time is increased by a factor of 4. For instance, the computer time for a typical simulation (for calculating the scatter contribution to a synthetic image) on the computer Alpha (AMD Opteron processor 250, 2.4 GHz; 6.26 GB RAM) in 40 x 40 points of interest, at the tube voltage 141 kV, with a precision of 1% (one standard deviation) takes approximately 17 hours. The statistical uncertainty has to be kept low when we are calculating images, since if the statistical uncertainty is too high, the scatter projection often contains artifacts when it is interpolated to a higher resolution.
3.4.2. Systematic uncertainties
There are several sources of uncertainty that affect the results. These include uncertainties in the cross sections, but also uncertainties in estimation of the different parameters in the input files. In an internal report, Ullman et al (2003) studied the effects of uncertainties in x‐ray spectrum half value layer (HVL), field size, grid lamella thickness and detector thickness. The conclusion from this study is that the systematic uncertainty due to these factors in the estimated εs/εp behind the grid is approximately 11%. Later, analysis of variance (ANOVA) and regression analysis was used to analyse similar uncertainties and the systematic uncertainty in εs/εp behind the grid was estimated to be approximately 9%. However, these relatively large uncertainties only apply in those cases where we attempt to mimic a real x‐ray imaging system and compare measured quantities from that system with our calculated quantities. In the cases where we only use simulations to study relative differences between alternative acquisition schemes for an imaging system, the stochastic uncertainty is more relevant.
Figure 3.6. Flow chart describing the most important steps in the Monte Carlo program Yes
Set-up geometry, voxel phantom, cross-sections, image detector and grid
Select photon energy Calculate contributions from primary photons analytically to point of interest
(Siddon’s algorithm) Start calculating contributions
from scattered photons Select direction of motion
of incident photon
Calculate free path with Coleman’s algorithm.
Calculate the contribution to collision density estimator from scattered photons.
Select type of interaction and assign weight Store energy imparted to the
phantom
Sample new direction of motion and continue photon history
Terminate by Russian
roulette?
Select new path length (Coleman’s algorithm) Yes Is the next interaction within the phantom? No Was this the
last history? Calculate scoring
4. ASSESSMENT OF IMAGE QUALITY
4.1. Introduction
Image quality assessment means quantifying the usefulness of an image to solve a specific diagnostic task. It is preferred if this image quality assessment is objective. There are, according to Barrett (1990), four criteria that are essential for objective assessment of image quality.
A. The task
Image quality can only be described in relation to a well‐defined task. This often means detection of an object in a structured or homogeneous background. A summary of different tasks to be solved in chest radiography is given in ICRU 70 (ICRU 2003). Among those, the search for malignant nodules at different positions in the lung is a common case treated in the literature (Samei et al 1999). In mammography it is common to search for calcifications or masses (Burgess et al 2001). Several different types of tasks are discussed in the literature (Barrett and Myers 2004). B. Image and object properties We need to understand the physical and statistical properties of the imaging system as well as the object being imaged. For example, the radiologist has an internal model of both the human anatomy as well as a large “data bank” of common pathologies in order be able to distinguish a malignant nodule from normal anatomy. The radiologist also needs to understand some of the physics behind the imaging technology in order to recognize some of the artefacts that may be present in the image. C. Observer
An observer that can perform the task is needed. This can be a human or a model observer. It is the human observer (radiologist) that is the final decision maker. Therefore, measures of image quality should always take the human observer into account. Yet, clinical trials involving human observers are costly and time demanding. Model observers may then be used to give insights into how image quality depends on the physical and technical image acquisition parameters.
D. Figure of merit
The figure of merit (FOM) is a number that tells us how well the observer performs the task. The FOM depends on the detection task; a commonly used FOM is the signal to noise ratio (SNR) or AUC, the area under the receiver operating characteristic curve (ROC). In some cases, image quality is described by physical properties derived from the image rather than by the performance of an observer on a specific task, for example, the modulation transfer function (MTF) or the noise power spectrum (NPS). We will refer to this as physical
image quality.
4.2. Image quality assessment as developed in this work
In this work, different tasks and figures of merit have been used, changing along with an improved modeling of the imaging system including the patient and improved model observers. In some cases, figures of merit based on pure physical measures of image quality were used and correlated to measures of the performance of human observers in clinical trials. A summary of the development is given below.
4.2.1. The task
Two main types of tasks are used in this work. The first task is to search for a known signal (e.g. a nodule) in a known background, referred to as the SKE/BKE (signal known exactly/background known exactly) task. The second task is to search for a known signal (e.g. a nodule) in a varying background, referred to as the SKE/BV (signal known exactly/background varying) task.
The SKE/BKE task was addressed in paper I, which was devoted to the optimization of tube voltage and filtration in iodine subtraction mammography. The task was to detect a blood vessel filled with iodine contrast of thickness 6 mg cm‐2 against a homogeneous background.
In paper IV the task was to compare two images from an x‐ray chest phantom and see in which image the structures corresponding to the image criteria were most clearly visible (the VGA study). In some sense, the background can be considered as known in this study since the same anthropomorphic phantom (the Alderson chest phantom) is used for all comparisons. Thus, the observer may be able to remember the background. In addition, the phantom used in this study does not contain as fine and complex structural details as are present in real phantoms (see further figure 5.4).
