2.1.1 Positron emission and annihilation

      Proton-rich isotopes may decay via positron emission, where a proton in the nucleus decays to a neutron, a positron and a neutrino. The daughter isotope has an atomic number one less than the parent. Examples of isotopes which undergo decay via positron emission are shown in table 2.1.

      As positrons travel through human tissues, they give up their kinetic energy mainly by Coulomb interactions with electrons. As the rest mass of the positron is the same as that of the electron, positrons may undergo large deviations in direction with each Coulomb interaction, and they follow a tortuous path through the tissue as they give up their kinetic energy (Fig. 2.2).

      

Table. 2.1. : Physical characteristics and production method of commonly used positron emitting radioiosotopes
Radioisotope	Half-life (min)	Maximum positron energy (MeV)	Positron range in water (FWHM in mm)	Production method
11C	20.3	0.96	0.28	cyclotron
13N	9.97	1.19	0.35	cyclotron
15O	2.03	1.70	1.22	cyclotron
18F	109.8	0.64	0.22	cyclotron
68Ga	67.8	1.89	1.35	generator
82Rb	1.26	3.35	2.60	generator

      When positrons reach thermal energies, they start to interact with electrons either by annihilation, producing two 511 keV photons which are anti-parallel in the positron's frame, or by the formation of a hydrogen-like orbiting couple called positronium. In its ground-state, positronium has two forms: ortho-positronium, where the spins of the electron and positron are parallel, and para-positronium, where the spins are anti-parallel. Para-positronium again decays by self-annihilation, generating two anti-parallel 511 keV photons. Ortho-positronium self-annihilates by the emission of three photons. Both forms are susceptible to the "pick-off" process where the positron annihilates with another electron. Free annihilation and the pick-off process are responsible for over 80% of the decay events. Variations in the momentum of the interacting particles involved in free annihilation and pick-off result in an angular uncertainty in the direction of the 511 keV photons of around 4 mrad in the observer's frame (Rickey et al 1992). In a PET tomograph of diameter 100 cm and active transaxial FOV of 60 cm, this results in a positional inaccuracy of 2-3 mm. The finite positron range and the non-collinearity of the annihilation photons give rise to an inherent positional inaccuracy not present on conventional single-photon emission techniques. However, other characteristics of PET which are discussed below more than offset this disadvantage.

      

Fig. 2.2. : Schematic representation of positron emission and annihilation

2.1.2 Positron emitting radiopharmaceuticals

      Nuclear medicine imaging involves studies of the localisation of administered radiopharmaceuticals by detection of the radiation emitted by the radionuclides. The radiopharmaceuticals can be administered to the patient by injection, orally or by inhalation. A radiopharmaceutical consists of a radionuclide with suitable radiation characteristics which is usually 'labelled' in a chemical reaction to a carrier. When labelled to a carrier, the main interest is the physiological and biological behaviour of the carrier and the radionuclide serves only as a tracer. There are several advantages of using radionuclides as tracers for in vivo studies. The physiological behaviour and the metabolism in different organs can be studied rather than the morphological structure. Important objectives such as the detection of lesions can be realised. However, since in most disease processes physiological or metabolic changes precedes changes in anatomy, nuclear medicine has a capability to provide earlier diagnosis.

      While PET was originally used as a research tool, in recent years it come to have an increasingly important clinical role. The largest area of clinical use of PET is in oncology where the most widely used tracer is ¹⁸F-fluoro-deoxy-glucose (¹⁸F-FDG). This radiopharmaceutical is relatively easy to synthesise with a high radiochemical yield. It also follows a similar metabolic pathway to glucose in vivo, except that it is not metabolised to CO₂ and water, but remains trapped within tissue. This makes it well suited to use as a glucose uptake tracer. This is of interest in oncology because proliferating cancer cells have a higher than average rate of glucose metabolism. ¹¹C-methionine is also used in oncology, where it acts as a marker for protein synthesis.

      PET has applications in cardiology, where ¹¹N-NH3 is used as a tracer for myocardial perfusion. When ¹³N-NH3 and ¹⁸F-FDG scans of the same patient are interpreted together, PET can be used to distinguish between viable and non-viable tissue in poorly perfused areas of the heart. Such information can be extremely valuable in identifying candidates for coronary by-pass surgery. PET has been used in neurology in a wide range of conditions, and in particular in severe focal epilepsy, where it may be used to compliment magnetic resonance imaging. Another reason for the importance of PET lies in the fact that, unlike conventional radiotracer techniques, it offers the possibility of quantitative measurements of biochemical and physiological processes in vivo. This is important in both research and clinical applications. For example, it has been shown that semi-quantitative measurements of FDG uptake in tumours can be useful in the grading of disease (Hoh et al 1997). By modelling the kinetics of tracers in vivo, it is also possible to obtain quantitative values of physiological parameters such as myocardial blood-flow in ml/min/g or FDG uptake in mmol/min/g providing the acquired data is an accurate measure of tracer concentration. Absolute values of myocardial blood flow can be useful in, for example, the identification of triple-vessel coronary artery disease and absolute values of FDG uptake can be useful in studies of cerebral metabolism. Some examples of radiopharmaceuticals used in PET are shown in table 2.2, together with some typical clinical and research applications.

      Prospects for the further development of several agents currently being investigated have been improved, and a number of established clinical procedures have developed an increased potential with the added application of the PET modality. These complementary roles of PET all contribute to an increased accuracy in nuclear medicine investigations.

      

Table 2.2. : Examples of positron-emitting radiopharmaceuticals used in clinical PET studies with some representative investigations where quantification is of interest
Radioisotope	Tracer compound	Physiological process	Clinical application
	methionine	protein synthesis	oncology
11C	flumazenil	bezodiazepine receptor antagonist	epilepsy
	raclopride	D2 receptor agonist	movement disorders
13N	ammonia	blood perfusion	myocardial perfusion
15O	carbon dioxide water	blood perfusion blood perfusion	brain activation studies brain activation studies
	Fluoro-deoxy-glucose	glucose metabolism	cardiology, neurology, oncology
18F	Fluoride ion Fluoro-mizonidazole	bone metabolism hypoxia	oncology oncology - response to therapy

2.1.3 Interaction of photons with matter

      When electromagnetic radiation passes through matter, some photons will be totally absorbed, some will penetrate the matter without interaction and some will be scattered in different directions, with and without loosing energy. Unlike charged particles, a well-collimated beam of photons shows a truly exponential absorption in matter. This is because photons are absorbed or scattered in a single event. That is, those collimated photons which penetrate the absorber have had no interaction, while the ones absorbed have been eliminated from the beam in a single event. For the energies of interest in nuclear medicine imaging, a photon interacts with the surrounding matter by four main interaction processes. The probability for each process depends on the photon energy, E, and the atomic number, Z, of the material.

---

---

---

---

---

2.1.4 Types of coincidence events

      Coincidence events in PET fall into 4 categories: trues, scattered, randoms and multiples. Fig. 2.5. illustrates the detection of first three of these different types of coincidences.

      

Fig. 2.5. : Types of detected coincidence events in PET

      True coincidences occur when both photons from an annihilation event are registered by detectors in coincidence, neither photon undergoes any form of interaction prior to detection, and no other event is detected within the coincidence time-window.

      A scattered coincidence is one in which at least one of the detected photons has undergone at least one Compton scattering event prior to detection. Since the direction of the photon is changed during the Compton scattering process, it is highly likely that the resulting coincidence event will be assigned to the wrong LOR. Scattered coincidences add a background to the true coincidence distribution which changes slowly with position, decreasing contrast and causing the isotope concentrations to be overestimated. They also add statistical noise to the signal. The number of scattered events depends on the volume and attenuation characteristics of the object being imaged, and on the geometry of the tomograph. Basically, we can distinguish single scattered events when one or both annihilation photons undergo one scattering and multiple scattered events where one of the two photons undergo two or more scatterings.

      Random (or false) coincidences occur when two uncorrelated photons not arising from the same annihilation event are incident on the detectors within the coincidence time window of the system. The number of random coincidences in a given LOR is closely related to the rate of single events measured by the detectors joined by that LOR and the rate of random coincidences increase roughly with the square of the activity in the FOV. As with scattered events, the number of random coincidences detected also depends on the volume and attenuation characteristics of the object being imaged, and on the geometry of the tomograph. The distribution of random coincidences is fairly uniform across the FOV, and will cause isotope concentrations to be overestimated if not corrected for. Random coincidences also add statistical noise to the data.

      A simple expression relating the number of random coincidences assigned to an LOR to the number of single events incident on the relevant detectors can be derived as follows. Let us define t as the coincidence resolving time of the system, such that any events detected with a time difference of less than t are considered to be coincident. Let r₁ be the singles rate on detector channel 1. Then in one second, the total time-window during which coincidences will be recorded is 2t r₁. If the singles rate on detector channel 2 is r₂, then the number of random coincidences R₁₂ assigned to the LOR joining detectors 1 and 2 is given by :

      

(2.11)

      This relation is true provided that the singles rate is much larger than the rate of coincidence events, and that the singles rates are small compared to the reciprocal of the coincidence resolving time t, so that dead-time effects can be ignored.

      Multiple coincidences occur when more than two photons are detected in different detectors within the coincidence resolving time. In this situation, it is not possible to determine the LOR to which the event should be assigned, and the event is rejected. Multiple events can also cause event mis-positioning.

3.1.1 Solution of the 2D problem

      The fundamental relation that links the imaged object (x,y) to its projections p_F(s) is the central-section theorem, which can be expressed in 2D as (Natterer 1986):

      

Fig. 3.2. : Illustration of the 2D central slice theorem showing the original 2D distribution f(x,y) and its 2D Fourier transform F(n_x,n_y) generally built from discrete sample points from P(n_s, F), the Fourier transform of the 1D projection data p(s,F)

      

(3.1)

      where capital letters denote Fourier transforms, and n_s denotes the coordinate on a frequency axis rotated by F relative to n_x-axis. Fig. 3.2 illustrates the various coordinates used in the derivation of the theorem. This yields the following final result:

      

, where

(3.2)

      Using this equation, (x,y) may be recovered from p_F(s) in two steps: (i) filtering: the projections are filtered by the ramp filter |n_s| in Fourier space; then (ii) backprojection: the contribution to the image volume of each of the filtered projections is computed.

3.1.2 Solution of the 3D problem

      The main difference between two- and three-dimensional image reconstruction from projections can be explained as follows (Egger 1996). Line integrals used to reconstruct the unknown source distribution (x,y), in the 2D case are characterised by two parameters: a slope and an intercept. The problem is analytically well defined and the solution exists and is unique. Whereas, in the 3D case, the image to be reconstructed, (x,y,z), depends only on three spatial variables, while line integrals are characterised by four parameters, two slopes and two intercepts. Hence, the problem becomes over-determined, since analytically the set of direct projections are sufficient to recover the source distribution in the conventional slice-by-slice manner. The additional information contained in oblique projections is analytically useless, as long as we the estimation of statistical variations is not of interest. Therefore, an infinite number of analytically equivalent reconstruction methods may be designed since the recovery procedure is no longer unique. It is worth to point out that the behaviour of these different methods may vary significantly when applied to real noisy data. Defrise (1995a) has, however, demonstrated that a recovery condition can be derived, which must be satisfied by any exact analytical algorithm.

      Again, as in the 2D case, the central section theorem is the fundamental relation that links the projection data to the image: a section of the 3D Fourier transform of the image through the origin and perpendicular to a unit vector is equal to the 2D Fourier transform of the projection of the image along . In mathematical form, this can be written as (Labbé et al 1997) :

      

(3.3)

      where capital letters denote Fourier transforms. Consequently, every measured projection samples one plane through the origin of the image in frequency space. In the case where the Fourier transform of the object to be reconstructed is known for all the frequencies, then this object is completely determined in the spatial domain. The condition that needs to be satisfied for the recovery of from is that the set of all the projections is such that for every frequency there exists at least one projection which satisfies the condition . However, these planes generally intersect and thus some frequency components are available from more than one projection. This illustrates again the over-determination of the 3D problem. A filter function is used to compensate for the fact that different frequencies are sampled more or less often by weighting every frequency component by the number of projection planes that contain .

      Fig. 3.3 summarises the reconstruction algorithms proposed so far in PET which fall as explained before within two main classes: direct analytical methods and iterative methods.

      

Fig. 3.3. : Summary of image reconstruction algorithms used in PET

3.2.1 The reprojection algorithm (3DRP)

      Three-dimensional reconstruction is complicated by the fact that all the oblique 2D projections are incompletely measured due to the finite axial extent of the tomograph. This can be compensated by forward projecting a first image estimate, reconstructed from the non-truncated direct projections. In practice, as much as 40% of the total reconstruction time is spent in estimating, and then backprojecting the unmeasured projection data. Fig. 3.4 illustrates the truncation problem during the acquisition of oblique projections for three different ring differences d.

      The algorithm combining this forward projection step with the Colsher filter is a 3D filtered backprojection algorithm (Kinahan and Rogers 1989) and is referred to here as 3DRP. Up to now, this is the most widely used reconstruction algorithm in clinical environments and consists of the following steps :

1. unmeasured oblique projection data, not accessible within the finite axial extension of the scanner, are estimated by forward-projecting through a low-statistics image reconstructed by 2D-FBP from transaxial projections;
2. the completed 2D projections are then reconstructed by the FBP technique where 2D projections are convolved with a 2D filter kernel, and then backprojected in 3D through the image volume.

      

Fig. 3.4. : Illustration of the projection truncation problem: three different cases are illustrated here during the acquisition of direct and oblique projections with a given ring difference d

3.2.2 The Fast Volume Reconstruction (FAVOR) algorithm

      The FAst VOlume Reconstruction (FAVOR) algorithm provides an exact method that does not require all 2D projections to be completely measured (Defrise et al 1992). FAVOR presents different noise propagation properties compared to other approaches, notably those using the Colsher filter referred to in the preceding section. In particular, it exhibits a low-frequency 'hump' characteristic which includes most of the useful frequency range occurring in a typical image. This can be mathematically modified to introduce a specific bias by attributing a much greater weight to the direct projections than to the oblique ones. High-frequency statistical noise is, however, not affected since it occurs at frequencies superior to those corrected by the FAVOR filter.

3.2.3 Rebinning-based algorithms

      The large amount of calculations required by truly 3D reconstruction algorithms such as those discussed in the previous sections has motivated the development of fast approximate methods for incorporating all the information measured by a volume imaging scanner into the reconstructed image. An improvement of the 3D reconstruction time requires some form of data compression to reduce the huge number of data samples. Well known approaches consist in either undersampling the data by averaging groups of contiguous samples (mashing) or compressing axially the oblique projections (span). Another approach is considered in this section: rebinning. This allows the processing of the 3D data as 2D slices and thus speeds up the reconstruction by an order of magnitude. The method used to rebin the data from the oblique sinograms into the smaller 2D data sets is the crucial step of any rebinning algorithm and has a vital influence on the accuracy and variance of the reconstructed image.

      The single-slice rebinning (SSRB) algorithm (Daube-Witherspoon and Muehllehner 1987) is based on a very simple approximation but its accuracy is acceptable only near the axis of the scanner and for small apertures. The SSRB assigns an LOR between two detector rings to the sinogram of the transaxial slice lying midway, axially between the two rings. A sinogram s_k is synthesised for each transaxial slice k=1,...,2N_rings-1 where s_k is the average of the sinogram S_ij between pairs of rings i,j located symmetrically, axially, with respect to slice k. Each slice k is then reconstructed by applying 2D FBP to the corresponding sinogram s_k.

      Lewitt et al (1994) proposed a more sophisticated rebinning approach, the multi-slice rebinning algorithm (MSRB). In this approach, a sinogram S_ij contributes to all slices k such that the LORs between diametrically opposite detectors in rings i and j intersect the slice k within a cylindrical field-of-view. The sinogram s_k is then built as a weighted average of all sinograms S_ij which satisfy that condition. After 2D reconstruction of the 2N_rings-1 sinograms s_k, each column of voxels is then deconvolved using the multiplicative iterative algorithm described by Lewitt et al (1994). MSRB is more accurate than SSRB, but suffers from instabilities in the presence of noisy data.

      The Fourier rebinning algorithm is a more sophisticated algorithm for 3D image reconstruction in which oblique rays are binned to a transaxial slice using the frequency-distance relationship of the data in Fourier space (Defrise 1995b). This algorithm rebins the 3D data into a stack of ordinary 2D direct sinograms and the image is reconstructed in 2D slice-by-slice, for instance using the standard 2D FBP algorithm (FORE+FBP) or OSEM (FORE+OSEM).

3.3.1 Algebraic reconstruction algorithms

      The Algebraic Reconstruction Technique (ART) algorithm uses a successive over-relaxation method known in numerical linear algebra (Herman and Meyer 1993). In fact, ART consists of guessing at a value for all the voxels f_j, and then modifying each element along each ray by a factor which compensates for the discrepancy between the measured ray sum p_F(s) and the calculated ray sum r_F(s):

      

(3.9)

      If the calculated ray sum is the same as the measured value, it implies that the guessed value is correct for a particular projection, however, for another projection there might be a large discrepancy. Thus the voxels of the last views while lying in the ray for the new view will be modified according to the discrepancy between the new ray and the measured value. Thus, each ray from each projection is examined and values of f falling within that ray are changed iteratively for all the projections for a certain number of iterations. Eq. (3.9) is called multiplicative ART. Another method of correcting the discrepancy between the measured and calculated projections consists of adding the difference between them. This is called the additive form of ART. It is evident that the computational effort of the ART algorithm is relatively low, but data matrices can be several hundred of Megabytes in size.

      An important modification of ART consists of setting to zero those values in the array that are clearly zero because they correspond to a ray sum that was observed as zero. This effective bounds the data and is an important boundary condition for any of the iterative techniques.

3.3.2 Statistical reconstruction algorithms

      Since PET measurements are based on a counting process, a reasonable statistical model is that the measurements have independent Poisson distributions :

      

, i = 1, .., N

(3.10)

      where are the means of the independent Poisson variables .

---

---

---

---

3.4.1 Data structure

      The hierarchy of data structures implemented in the library is presented in Fig. 3.6. This is not a hierarchy of class derivation, but of containment. Let a PETSegment be the 4D data structure by abuse of terminology. A PETSegment may be expressed in two flavours or derived classes (PETSinogram and PETViewgram), depending on the order of the coordinates View, Ring and Bin. (View is the azimuthal angle of the projection (angle), bin the radial coordinate on the projection plane (bin element), and Ring the unit in the axial direction of the scanner). Thus, the user has the capacity to decide depending on the scanner geometry what is the best and fastest way to organise the data.

      The PETSinogramOfVolume class represents the whole 3D data set before reconstruction and is always organised into an odd number of segments, where each segment contains a specific number of sinograms or viewgrams. The segment identification numbers for a sinogram with 2N+1 segments are: 0, -1, +1, -2, +2 ...., -N, +N. A well-known fact in C or C++ programming is worth repeating here: care must be taken when choosing the order of indices in multidimensional arrays: the fastest varying index last, the slowest first. Due to the large number of memory accesses this has a considerable effect. In a previous implementation (Egger 1996), it was found that using sinogram arrays indexed by [Ring][View][Bin] and an image array indexed by [plane][y][x] results in the backprojection executing more than 15% faster than when these indices are used back to front (x and y are the transaxial image coordinates). The present implementation uses viewgrams as input for the projectors described below (ring and bin are the coordinates in the inner loops).

      

Fig. 3.6. : Hierarchical view of the data structures

3.4.2 The projection/backprojection operators

      The forward projection operator allows to estimate projection data given an image. The forward projection operator implemented in our library is based on the implementation of Siddon's algorithm (Siddon 1985), calculating the length of intersection of LORs with image voxels and making use of geometrical symmetries of the image volume. Every projection element p_i is related to the image voxel values f_j according to Eq. (3.8). For a given LOR i, the length of intersection of this LOR with a given voxel ij in the image volume gives the weight a_ij of that particular voxel in the projection element. More accurate but computationally intensive strip integrals, covering the entire detector width, are often used in 2D, but are impractical for use in 3D. The image representation adopted within the library consists in taking the z-voxel size equal to half the ring spacing of the scanner to improve resolution in the z-direction. Therefore, the voxel size rather than the detector ring width must be used for the axial sampling distance of the forward projection operator.

      Symmetry properties of geometrical factors, such as lengths of intersection and voxel positions, are fully exploited to reduce the computational burden of the algorithm. Symmetry operations act on all coordinates of a LOR, its azimuthal angle F, its co-polar angle q, and its radial and axial coordinates s and z respectively. Details of these may be found in (Egger 1996). Using these symmetry properties, the lengths of intersection computed for one voxel may be re-used straight away for as many as fifteen other voxels times the number of projected sinograms for a given ring index difference.

      The backprojection operator allows to estimate the image volume given a projection data. The backprojection operator implemented in the library is based on an extension of Cho's incremental algorithm (Cho et al 1990) to 3D cylindrical geometry. This algorithm is called the 3D beamwise, incremental backprojection method (Egger et al 1998). Bi-linear interpolation between four projection elements is used to update voxel values. A beam is defined as the volume delimited by the central rays of four adjacent projection elements. Sweeping the image volume beam by beam, allows to make use of a number of beam constants.

      A searching flow algorithm is used to locate the voxels inside a beam, starting from the first voxel found. Since the beam's axial extension is equal to two voxels, a restriction to finding only one half of the actual number of voxels inside the beam can be made. Symmetry properties of the image volume may be exploited to reuse all geometrical quantities related to a particular voxel, such as the interpolation distances ds and dz, for all voxels in equivalent geometrical environments, without having to calculate them explicitly. A very fast implementation of the backprojection operator has been achieved by combining the incremental algorithm with the exploitation of the geometrical symmetries of the image volume (Egger 1996).

      In the case where voxel size is not equal to the detector size, a generalisation of the partially incremental algorithm (Egger et al 1998) to the piece-wise linear case for this choice of voxel size was recently proposed (Thielemans et al 1999). Obviously, the piecewise bilinear interpolation corresponds much better to the probability than bilinear interpolation. Moreover, it will improve resolution. The increase in CPU time for the incremental versions of piece-wise bilinear over ordinary bilinear interpolation is not considerable. When the detector size is larger than the voxel size, it has been shown that the kernel of the piecewise bilinear interpolation has a smaller spread than in bilinear interpolation, and hence the backprojected images have a higher resolution. The new piece-wise bilinear algorithm provides an efficient alternative for the Tube of Response (TOR) model. When voxel sizes are equal to detector sizes, the piecewise bilinear interpolation simply reduces to ordinary bilinear interpolation. 3.4.3 Pre-computation of the sensitivity image

      An important component of many iterative algorithms is the sensitivity image , which computes the probability that an event emitted in a given a voxel is detected or not. For current generation scanners, the number of elements can be prohibitively large. When the image and sinogram grid sizes are chosen as discussed previously, table 3.1 illustrates the predicted size of the probability matrix that can result for both the GE-Advance and the ECAT EXACT3D PET tomographs. The number of independent elements refers to the number of relevant entries of the a_ij matrix after symmetry reduction.

      Various models for the calculation of detector-pair sensitivities have been proposed in the literature. These include the solid angle model, the model based on the volume of intersection of TOR and the voxel, the length of intersection of the 'centre' LOR with the voxel model, and interpolation. The sensitivity image can be computed as the backprojection of the transition matrix A set all equal to 1 (Thielemans et al 1999). Therefore the smoothness assumption - the variation of the image intensities over a scale set by the detector sizes is very smooth - is optimally satisfied. For a given model for the detector-pair sensitivities, the approximation that the discrete backprojection can be computed by sampling the continuous backprojection, times a factor which takes into account the volume of the intersected voxel will give (nearly) the same result as the exact calculation.

      

Table 3.1. : Estimated sizes of the probability matrix data for two large PET tomographs: the GE-Advance and the EXACT3D
	GE-Advance	EXACT3D
Total number of elements	~1013	~1016
Number of non-zero elements	~1010	~1012
Number of independent elements	~109	~1011

3.5.1 COSEM : an iterative algorithm for fully 3D listmode data

      As explained in Chapter 2, listmode data acquisition of coincident annihilation photons stores events in chronological order of their registration. Therefore, the usual sinograms of PET scans can be considered as a dense representation of the registered data and the listmode as a sparse representation of the same data. In low statistics studies, the listmode representation is more compact than the sinogram representation. This is because the sinogram representation is fixed according to the geometry of the scanner (billions of LORs) and does not depend on the number of registered events.

      A new algorithm called Coincidence list OSEM (COSEM) was recently developed in a joint collaboration between the Technion Institute of Technology, Haifa, and Elscint (now Elgems) to exploit this acquisition mode successfully by obtaining better signal-to-noise rate in a timely way. It is based on patent-pending proprietary Elscint development, nevertheless the algorithm has been improved and refined with the PARAPET project (Levkovitz et al 1999). The COSEM is based on the OSEM algorithm for binned data but has several extensions that make it suitable for rotating, dual-head planar detector PET tomographs. This representation has several major advantages over standard OSEM representation :

1. the algorithm uses the exact detection points as registered in the listmode data and does not degrade positioning accuracy by binning the data;
2. the number of operations is proportional to the number of detected events;
3. the large four-dimensional sinogram representation required for true 3D reconstruction using binned data is no longer needed;
4. partitioning the ordered sets is simple and straightforward. A subset can be made of groups of events collected in the listmode data (e.g. randomly selected events, acquisition during a certain period of time, events acquired at certain angles, etc.)

      Since the data sets acquired on dual-head scanners are not binned to projection data, correction for the influences of different physical and geometrical factors should be performed either on the registered LORs or taken into account when calculating voxel detection probabilities. The application of appropriate weights to the LORs renders the correction process simpler and more computationally efficient. An exact analytical calculation of each voxel's geometrical detection probability is included in the COSEM implementation. Therefore, patient-dependent corrections are applied directly to LORs, while scanner-dependent corrections are calculated as voxel detection probabilities.

      COSEM was adapted to the Elgems coincidence gamma cameras and an evaluation on phantom and clinical studies performed in comparison to the standard manufacturer's FBP algorithm. Fig. 3.8 shows a comparison between images of the thorax phantom as reconstructed by standard filtered backprojection (left) and by the optimised COSEM algorithm (right). It has been shown that COSEM reconstructions have better signal to noise ratio and much better visibility of small tumours. Moreover, an important gain in speed has resulted by running the algorithm after the acquisition of the first events and making the reconstructed images available shortly after completion of data acquisition.

      

Fig. 3.8. : Comparison between images of the Lung and Heart phantom with 4 inserts (Duke Hospital, USA) reconstructed by standard filtered backprojection (left) and images reconstructed using COSEM optimised iterative reconstruction (right)

      The smallest observable sphere has a radius of 8 mm (Courtesy of Dr R. Levkovitz, Technion).

3.5.2 OSEM with inter-update filtering

      An enhanced version of OSEM called Inter-Update Metz Filtered OSEM (IUMF-OSEM), which uses an iterative filtering technique for enhancing image quality was developed by Jacobson et al (1999). In this technique, a Metz filter is applied to image estimates at certain intervals. This filter is known to carry out restorative image enhancement on the image estimates, and the proposed algorithm incorporates it in OSEM's updating mechanism.

      A software implementation of IUMF-OSEM was performed using the previously described object-oriented library for 3D PET reconstruction (Labbé et al 1999b). The algorithm combines the forward and back projection operators with an appropriate choice of ordered sets to fully use all symmetry properties for fast computation.

      The performance of IUMF-OSEM Vs standard OSEM and post-filtered OSEM reconstructions was tested on phantom and clinical studies. The phantom scans were acquired on the GE-Advance PET scanner. Data sets were pre-corrected for attenuation, detector efficiency, scatter, and random coincidences. Iterative filtering operations were applied at intervals of full iterations. Fig. 3.9 shows sample images of the Jaszczak phantom reconstructed with 14 subsets and 5 iterations (70 sub-iterations).

      Quantitative evaluations of the algorithm on phantom data showed that inter-update Metz filtering, with appropriately chosen filter parameters yields superior signal-to-noise ratio and contrast recovery compared to plain OSEM. Moreover, contrast enhancement was found to be higher, and noise control only slightly poorer than that obtained by post-filtered OSEM. The contrast enhancing capability of IUMF-OSEM was also demonstrated on reconstructions of the whole body phantom and on clinical studies.

      

Fig. 3.9. : Illustration of the Jaszczak phantom with rods insert (left) and a sample IUMF-OSEM reconstruction (slice 16) after 70 iterations with Metz filter power 1, radial FWHM is 5 mm and axial FWHM is 8 mm (right)

(Jacobson et al 1999)

3.5.3 The mirror descent algorithm

      A new numerical method called Ordered Subsets Mirror Descent (OSMD) for solving a non-smooth convex optimisation problem over a compact feasible set was recently introduced and applied to the problem of image reconstruction in 3D PET (Ben-Tal et al 1999). This is an accelerated version of the Mirror Descent (MD) method developed by Nemirovski and Yudin (1978). The MD method can be regarded as a non-Euclidean version of the classical subgradient method, but has a significant theoretical advantage since its rate of convergence is independent (or nearly so) of the dimension of the problem. This feature makes it an attractive choice for large-scale problems like volumetric reconstruction of PET data.

      The OSMD algorithm has been adapted to the PET reconstruction by reformulating the optimisation problem. This leads to an objective function and a constraining set with desirable properties. A proof of convergence of the OSMD method for the unconstrained case and for a constrained compact convex domain together with more details about the algorithm are described elsewhere (Ben-Tal et al 1999). A principal advantage of OSMD is that convergence to the maximum likelihood estimate is guaranteed for any choice of subsets and for any projection data. OSEM on the other hand may cycle among sub-optimal solutions if the data are noisy and/or the subsets don't satisfy certain balancing conditions described in Hudson and Larkin (1994). It is worth to point out that the preliminary results presented here come from a prototype implementation. A more comprehensive evaluation is in preparation by the authors of this algorithm. Fig. 3.10 shows slices of a Monte Carlo simulated Jaszczak phantom with rods insert for the PRT-1 tomograph reconstructed with both OSMD (left) and standard OSEM (right). The OSMD reconstruction exhibits high resolution even for the block with the smallest diameters.

      

Fig. 3.10. : One slice of a simulated Jaszczak phantom with rods insert reconstructed with OSMD (left) and standard OSEM (right) after 5 iterations with 6 subsets

(Ben-Tal et al 1999)

3.5.4 Inverse Monte Carlo-based reconstruction

      The concept of inverse Monte Carlo (IMC) was introduced in 1981 in an attempt to describe a numerical method for solving a class of inverse problems (Dunn 1981). The IMC method converts the inverse problem, through a non-iterative simulation technique, into a system of algebraic equations that can be solved by standard analytical or numerical techniques. The principal merits of IMC are that, like direct Monte Carlo, the method can be applied to complex and multivariable problems, and variance reduction procedures can be applied. The non-iterative IMC is strongly related to the variance reduction technique in direct simulation called importance sampling where the sampling process is altered by using random numbers from a modified distribution.

      Application of the Monte Carlo method for solving the inverse problem in emission computed tomography has been previously reported by other groups (Floyd et al 1986). As proposed, the method realises tomographic reconstruction for SPECT with simultaneous compensation for attenuation, scatter, and distance-dependent collimator resolution. A detection probability matrix was created by Monte Carlo solution to the photon transport equation for a simulated SPECT acquisition from a unit source activity in each source voxel. The measured projection vector equals the product of this detection probability matrix with the unknown source distribution vector. The resulting large, non-sparse system of equations was solved for the source distribution using an iterative ML-EM estimator. It is worth to note that although the technique was developed for SPECT, it is also valid for other imaging techniques like PET and transmission CT.

      A strategy for developing a 3D reconstruction algorithm for PET based on IMC taking advantage of the fact that Monte Carlo methods provide a natural framework for incorporation of detector response function, the effects of attenuation, scatter and other sources of systematic errors in the transition matrix is presented here and further discussed in future developments section (chapter 8). Our proposed application of IMC is in the context of the 3D ML-EM reconstruction algorithm (see section 3.2.2). For each iteration, the update to the source distribution is done by comparing the forward projected data of the current image estimate to the measured projections. The IMC procedure replaces the explicit forward projection step given by Eq. (3.12) in the ML-EM algorithm as by Monte Carlo generated events according to a hypothesised source distribution and attenuation map. A uniform image could be used as the initial source estimate for the first iteration. In the following iterations, voxel values are updated by the ratio of the backprojection of the measured data to the Monte Carlo generated data at that voxel. The total number of detected events is scaled to match the original measured data sets.

      A preliminary version of the IMC algorithm was implemented using the 3D reconstruction library (Labbé et al 1999b). Our approach may be distinguished from the method of Floyd et al (1986) and has been discussed as a potential application of our Monte Carlo simulator (Zaidi et al 1998), although other groups reported some related work in progress in the context of the image space reconstruction algorithm (Kudrolli et al 1999). One of the limitations of their work is the rather crude system model used which contains only the basic detector geometry. That is, the system response matrix is approximated by the length of intersection of the i^th LOR with the j^th voxel. A major advantage of our approach is that is avoids the explicit calculation and storage of the system response matrix since Monte Carlo simulation of a given source/scattering/detector system implicitly contains all the information contained in the detection probability matrix. The probability matrix is usually much too large to store in memory and has to be either stored on disk or calculated on the fly. Several groups have used symmetries and sparsness of this matrix to be able to store it in memory (Qi et al 1998).

      Obviously, IMC is a reconstruction algorithm which is based on a more realistic model of the PET acquisition process. The most dramatic improvements using IMC are expected to be apparent for low statistics studies. This is significant since clinical PET studies are usually count limited. This improvement in image noise is expected to be to due to several characteristic features of IMC: the incorporation into the EM algorithm of the knowledge that the photon counting process follows a Poisson distribution, the simultaneous solution approach, and the detailed modelling provided by the Monte Carlo method reduces the data inconsistency. The most effective clinical application of IMC could be envisaged in studies where the signal to noise ratio is too low for conventional FBP or common iterative reconstruction techniques. Specific applications include fast dynamic studies and low statistics multi-bed oncologic studies. The most obvious disadvantage of the IMC algorithm is the increased computational effort required to perform the reconstructions. However, parallel implementation may help to overcome this drawback.

3.5.5 Parallelisation and scalability

      As stated previously, the major goal of the PARAPET project is to create parallel versions of both analytic and iterative algorithms to run on distributed memory multiprocessor parallel platforms. The partners in the project use Parsytec CC series of systems (Fig. 3.11), however, the main goal is to create parallel software tools that achieve:

1. speed-up on the basis of elapsed wall-clock time;
2. scale-up in terms of the size of data sets obtained on different scanners;
3. portability and platform independence, allowing to run the software on clusters of high-end PCs and workstations and other distributed memory multiprocessor hardware architectures.

      The parallelisation library developed by our partners (sadki et al 1999) provides the programmer with a common interface for creating parallel applications for different hardware and software platforms. Both PVM (Geist et al 1994) and EPX (EPX 1996) message-passing libraries are supported now, and it is planned that MPI will be supported in the near future. The library consists of a set of C++ classes and functions that define a common parallel programming interface independent of a particular message-passing library. Therefore, the application is definitely isolated from the used low-level library. Moreover, the implemented message-passing code can be used as it is when ported to a different hardware platform (as far as the platform is supported by the library).

      

Fig. 3.11. : Photograph of the 8-node Parsytec CC system installed at Geneva University Hospital

      Parallel implementation of the reconstruction algorithms uses this portable, object-oriented message passing library and a coarse-grain parallelisation approach based on partitioning the projection space. By considering similarity with respect to data distribution and communication patterns, four groups of parallel reconstruction algorithms were identified :

1. 3DRP, FORE, ...etc.
2. OSEM, OSEM/MAP, OSMD, BI-SMART;
3. pure ART and modifications of ART based on splitting the image volume;
4. B-ART.

      Detailed analysis of data flow allowed to identify and parallelise a common kernel. For example, 3DRP, which is an analytic algorithm, shares the same strategy with OSEM, but is classified in a separate group due to its differences from the OSEM-like schemes. With current and future developments, more groups of algorithms can be identified and corresponding parallel kernels implemented. For the present implementation, a master-slave architecture was adopted. The processing time depends on data sizes and hence on the segment number (ring difference) being processed. This results in a problem of even distribution of the workload which was handled by implementing on-the-fly workload distribution scheme. In this approach, the working nodes are assigned jobs as soon as they get free, instead of using a predefined partitioning scheme. Since the size of projection data decreases with increasing ring difference, processing ring differences in ascending order seems a reasonable solution. This leads to better load balancing and allows to reduce the granularity of required processing towards the end of a sub-iteration. The achieved speed-up will of course be influenced by the different sizes of projections. The strategy adopted for data decomposition consists in partitioning the projection space across the processors (Sadki et al 1999). For instance, in the case of OSEM, this approach assigns groups of symmetry-related viewgrams to different processors, and allows several iterations of the inner loop to be executed in parallel. Data consolidation involves summing the partial computations of the update factors calculated on different processors.

      Performance evaluation was carried out on the Parsytec CC system consisting of 12 nodes (see chapter 6 for a detailed description of the system). Reconstructions were carried out for data sets of both PRT-1 and GE-Advance tomographs. Scalability issues were investigated by varying the number of nodes from 3 to 12. The efficiency of the parallelisation approach for different numbers of nodes was assessed by comparing parallel timing results with those obtained by running the code in serial mode. In the present implementation, the master node is used only for data distribution and load balancing. Efficiency is therefore defined as , where N is the total number of nodes.

      Timing results and parallelisation efficiencies together with estimated timings for an 'ideal' case of linear speed-up defined as for two different scanners are presented in Fig. 3.12. As can be seen on the illustrated charts, the parallel implementation provides a reasonable speed-up and leads to reduced reconstruction times. In particular, the implementation on 12 nodes yielded a speed-up relative to serial mode of about 7 out of a possible 11. Comparison of efficiency charts for PRT-1 and GE-Advance scanners suggests that the implementation scales well with respect to problem size. The parallelisation of the algorithms constitutes a significant step toward achieving clinically acceptable reconstruction times.

      

Fig. 3.12. : Execution time per iteration (left) and efficiency (right) of parallel implementation of OSEM algorithm on the Parsytec CC system for both PRT-1 (top) and the GE-Advance (bottom) scanners

      Both actual execution time (circles) and estimated timings for an 'ideal' case of linear speed-up (diamonds) are shown (Sadki et al 1999).

4.2.1 Random numbers generation

      Computational studies requiring the generation of random numbers are becoming increasingly common. All random number generators (RNG) are based upon specific mathematical algorithms, which are repeatable. As such, the numbers are just pseudo-random. Here, for simplicity, we shall term them just 'random' numbers. Formally, random is defined as exhibiting 'true' randomness, such as the time between 'tics' from a Geiger counter exposed to a radioactive element. Pseudo-random is defined as having the appearance of randomness, but nevertheless exhibiting a specific, repeatable pattern. Quasi-random is defined as filling the solution space sequentially (in fact, these sequences are not at all random, they are just comprehensive at a preset level of granularity). Monte Carlo methods make extensive use of random numbers to control the decision making when a physical event has a number of possible results. The RNG is always one of the most crucial subroutines in any Monte Carlo-based simulation code. A large number of generators are readily available (Marsaglia and Zaman 1994), and many of these are suitable for implementation on any computer system since today there is no significant distinction in floating point processing capabilities between a modern desktop and a mainframe computer. A typical simulation uses from 10⁷ to 10¹² random numbers, and subtle correlations between these numbers could lead to significant errors (Ferrenberg et al 1992). The largest uncertainties are typically due more to approximations arising in the formulation of the model than those caused by lack of randomness in the RNG. Mathematically speaking, the sequence of random numbers used to effect a Monte Carlo model should possess the following properties (Vattulainen et al 1995):

1. Uncorrelated sequences: the sequences of random numbers should be serially uncorrelated. Most especially, n-tuples of random numbers should be independent of one another.
2. Long period: ideally, the generator should not repeat; practically, the repetition should occur only after the generation of a very large set of random numbers.
3. Uniformity: the sequence of random numbers should be uniform, and unbiased. That is, suppose we define n-tuples and divide the n-dimensional unit hypercube into many equal sub-volumes. A sequence is uniform if in the limit of an infinite sequence, all the sub-volumes have an equal number of occurrences of random n-tuples.
4. Reproducibility: when debugging programs, it is necessary to repeat the calculations to find out how the errors occurred. The feature of reproducibility is also helpful while porting the program to a different machine
5. Speed: it is of course desirable to generate the random numbers fast.
6. Parallelisation: the generator used on vector machines should be vectorisable, with low overhead. On massively parallel architectures, the processors should not have to communicate among themselves, except perhaps during initialisation.

      Although powerful RNGs have been suggested including shift register, inversive congruentional, combinatorial and 'intelligent' methods such as those implemented in the MCNP code (Booth 1988), the most commonly used generator is the linear congruential RNG (LCRNG) (Lehmer 1949). Recently, Monte Carlo researchers have become aware of the advantages of lagged Fibonacci series (LFRNG). With extremely long periods, they are generally faster than LCRNG and have excellent statistical properties (Marsaglia and Zaman 1993). Those generators are briefly described below.

---

---

4.2.2 The source and positron emission

      In most Monte Carlo computer codes, the positron histories and thus positron range are neglected during simulations with the assumption that the positron emission point is coincident with the annihilation point. The annihilation point in the source region and the flight direction for the first of the two annihilation 511 keV photons are randomly generated. The second annihilation photon is generated in the same emission position, but with anti-parallel flight direction with respect to the first one. Both the positron range and the non-collinearity effects are accounted for by analytic models and included with the other physical phenomena affecting spatial resolution.

4.2.3 Modelling photon transport

      For radiation transport problems, the computational model includes geometry and material specifications. Every computer code contains a database of experimentally obtained quantities, known as cross sections, that determine the probability of a particle interacting with the medium through which it is transported. Every cross section is peculiar to the type and energy of the incident particle and to the kind of interaction it undergoes. These partial cross sections are summed to form the total cross section; the ratio of the partial cross section to the total cross section gives the probability of this particular interaction occurring. Cross section data for the interaction types of interest must be supplied for each material present. The model also consists of algorithms used to compute the result of interactions (changes in particle energy, direction, etc.) based on the physical principles that describe the interaction of radiation with matter and the cross section data provided. Therefore, it is extremely important to use an accurate transport model as the Monte Carlo result is only as valid as the data supplied.

      When a photon (having an energy below 1 MeV) passes through matter, any of the three interaction processes (photoelectric, incoherent scattering, coherent scattering) may occur. The probability of a photon of a given energy E undergoing absorption or scattering when traversing a layer of material Z can be expressed quantitatively in terms of a linear attenuation coefficient (cm^-1) which is dependent on the material's density, r (g.cm^-3).

      In the case of photoelectric absorption, the total photon energy is transferred to an atomic electron and the random walk is terminated. In an incoherent photon interaction, a fraction of the photon energy is transferred to the atomic electron. The direction of the scattered photon is changed to conserve the total momentum of the interaction. The Klein-Nishina expression for the differential cross section per electron for an incoherent interaction is used to sample the energy and the polar angle of the incoherently scattered photon (Klein and Nishina 1928). The coherent scattering only results in a change in the direction of the photon since the momentum change is transferred to the whole atom. The kinetic energy loss of the photon is negligible. Coherent scattering of a photon could be generated using the random number composition and rejection technique (Kahn 1956) to sample the momentum of the scattered photon and the scattering angle according to the form factor distribution. A simplified flow chart of photon transport for the energies of interest in nuclear medicine imaging is given in Fig. 4.2.

      

Fig. 4.2. : A flow diagram showing the main steps followed when simulating a photon transport history

      It is common to neglect coherent scattering in PET Monte Carlo simulation of photon transport because of its low contribution to the total cross section at 511 keV. In the following examples, the relative importance of the various processes involved in the energy range of interest (below 1 MeV) are considered for some compounds and mixtures used in nuclear medicine to justify some of the approximations made in Monte Carlo codes. Fig. 4.3 illustrates the relative strengths of the photon interactions versus energy for water, cortical bone, sodium iodide (NaI(Tl)) and bismuth germanate (BGO), respectively. For water, a moderately low-Z material, we note two distinct regions of single interaction dominance: photoelectric below and incoherent above 20 keV. The almost order of magnitude depression of the coherent contribution is some justification for the approximations discussed. The coherent contribution to the total cross section is less than 1% for energies above 250 keV. However, this contribution is in the order of 7% for high-Z materials like BGO. Therefore, efforts should be made to treat the coherent scattering process adequately for detector materials.

      

      

      

      

Fig. 4.3. : Components of photon cross sections for different tissues (H₂O and cortical bone) and detector materials (NaI(Tl) and BGO) of interest in PET imaging illustrating the relative contribution of each process

      Calculation of the distances between interactions in a medium are performed by sampling from the exponential attenuation distribution (Eq. 4.8 below). Different techniques have been proposed to improve the computation speed when sampling from the probability distributions. They are described in more detail in the following sections.

      In principle, electron transport should be included when simulating the complete electromagnetic cascade (microscopic techniques). However, the large number of interactions that may occur during electron slowing down makes it unrealistic to simulate all the physical interactions (macroscopic techniques). Secondary electrons are generally assumed to deposit all their energy at the point of interaction because of the low energies involved in nuclear medicine and therefore the short ranges of the electrons generated and their negligible Bremsstrahlung production. Therefore, electron transport has not received particular attention in nuclear imaging applications of the Monte Carlo method. However, a number of investigators considered this effect mainly for dosimetry calculations (Furhang et al 1997).

4.2.4 Analogue sampling

      Analogue Monte Carlo attempts to simulate the full statistic development of the electromagnetic cascade. If we assume that a large number of particle histories, N, are included in a batch, the individual batch estimates can be considered as drawn from a normal distribution. For a given calculation, the estimated uncertainty is proportional to the inverse of the square root of the number of histories simulated. The efficiency of a Monte Carlo calculation can therefore be defined as (Bielajew and Rogers 1989):

      

(4.5)

      where T is the calculation time required to obtain a variance estimate . For large N, should be constant as long as the calculation technique remains the same.

      As discussed earlier, the imaging system can be described in terms of pdf's. These pdf's, supplemented by additional computations, describe the evolution of the overall system, whether in space, energy, time, or even some higher dimensional phase space. The goal of the Monte Carlo method is to simulate the imaging system by random sampling from these pdf's and by performing the necessary supplementary computations needed to describe the system evolution. In essence, the physics and mathematics are replaced by random sampling of possible states from pdf's that describe the system. Thus, it is frequently necessary to sample some physical event, the probability of which is described by a known pdf. Examples include the distance to the next interaction and the energy of a scattered photon. Let x be the physical quantity to be selected and f(x) the pdf. Among the properties of the pdf is that it is integrable and non-negative. Assume that the domain of f(x) is the interval [x_min, x_max] and that it is normalised to unit area. The cumulative distribution function F(x) of the frequency function f(x) gives the probability that the random variable t is less or equal to x. It is defined as:

      

,

(4.6)

      A stochastic variable can be sampled by the use of uniformly distributed random numbers R in the range [0-1] using one of the techniques described below.

---

---

---

4.2.5 Non-analogue sampling "Variance reduction techniques"

      A direct Monte Carlo simulation using true probability functions may require an unacceptable long time to produce statistically relevant results. Photons emission is isotropic, so directional parameters may be sampled uniformly within their individual ranges. Nuclear imaging systems have a low geometrical efficiency because of the small solid angle defined by the collimator and/or the small axial aperture. Therefore, the calculation would be very ineffective in terms of required computing time (Haynor et al 1990). It is thus desirable to bias the sampling (non-analogue sampling) by introducing different types of importance sampling and other variance reduction techniques to improve the computational efficiency of the Monte Carlo method (Zubal and Harell 1992). The results obtained by non-analogue simulation are, however, biased by the variance reduction technique and a correction for this is required. A particle history weight, W, is introduced, which describes the probability of the particle following the current path. This weight is calculated for each particle history, and used in the calculation of the results. If an event occurs, the weight W is added to the counter rather than incrementing the counter by one unit. Bielajew and Rogers (1989) divided variance reduction techniques in three categories: those that concern photon transport only, those that concern electron transport only, and other more general methods. The most useful techniques in photon transport are described below. It is worth to note that concerning the Poisson nature of the activity distributions in PET, these variance reduction approximations may result in statistically relevant deviations from an otherwise unbiased simulation result since Monte Carlo simulated projection data are (as a result of variance reduction) not count but probability data.

---

---

---

---

---

4.3.1 Detector modelling

      The critical component of emission tomography is the scintillation detector. Monte Carlo simulation of detector responses and efficiencies is one of the areas which has received considerable attention (Andreo 1991, Zaidi 1999a). An early contribution to the field providing a detailed description of the techniques used was due to Zerby (1963). Many detector modelling applications were developed in the PET field including the pioneering work of Derenzo (1981) who simulated arrays of detectors of different materials and sizes to study the effect of the inter-crystal septa and later on to optimise the optical coupling between BGO crystals and PMT's (Derenzo and Riles 1982) by taking into account the reflection and scattering along the detection system. The search for an appropriate detector for this imaging modality was conducted in a comparative study of several crystals including BGO, CsF and NaI(Tl) (Bottigli et al 1985), BaF₂ used in time-of-flight PET (Bice et al 1990), and liquid Xenon (Lopes et al 1995). Binkley (1994) modelled the impulse response of a PMT, front-end amplifier, and constant fraction discriminator to evaluate the effects of front-end bandwidth and constant fraction delay and fraction for timing system optimisations of BGO scintillation detectors.

      The penetration of annihilation photons into the detector material before interaction is a statistical process which leads to significant displacement and anisotropy of the point-spread function. Compensation for crystal penetration is thus an important issue to recover the spatial resolution in PET (Huesman et al 1989). Comanor et al (1996) investigated algorithms to identify and correct for detector Compton scatter in hypothetical PET modules with 3x3x30 mm BGO crystals coupled to individual photosensors. The true crystal of first interaction was determined by the simulation for eventual comparison with the crystal identified by a given algorithm. They reported a mis-identification fraction of 12% if the detector has good energy and position resolution when using position of interaction to identify forward scatter.

      Numerous strategies have been proposed for constructing detector modules that measure the depth of interaction (DOI), but most of them proved impractical to implement or provided insufficient DOI measurement resolution. Two important questions can be addressed through Monte Carlo simulation: (i) what fraction of events will be mis-identified because of noise fluctuations in the photomultiplier tubes (PMT's) or photodiode array and (ii) how will the DOI measurement resolution affect the reconstructed resolution of a PET camera? The position-dependent light distribution has been used to measure the 511 keV photon interaction position in the crystal on an event by event basis to reduce radial elongation (DeVol et al 1993). Different geometrical modifications were also simulated leading to a proposal of a 2.2x5x30 mm BGO crystal, for which a 2.2 mm FWHM light distribution is predicted, which should yield a PET detector module with DOI measurement resolution of 3.6 mm FWHM. A test module with one 3x3x30 mm BGO crystal, one 3 mm square PIN photodiode, and one PMT operated at -20º C with an amplifier peaking time of 4 µs, and a measured DOI resolution of 5 to 8 mm FWHM has been proposed by Moses and Derenzo (1994). Simulations predicted that this virtually eliminates radial elongation in a 60 cm diameter BGO tomograph. The performance of a single detector element must be extrapolated using Monte Carlo simulations to predict the performance of a multi-element module or a complete PET camera.

      The Triumph PET group has developed a simulation tool to model position encoding multicrystal detectors for PET that treats the interactions of energetic photons in a scintillator, the geometry of the multicrystal array, as well as the propagation and detection of individual scintillation photons (Tsang et al 1995). Design studies of a whole-body PET tomograph with the capacity to correct for the parallax error induced by the DOI of annihilation photons were also performed (Moisan et al 1995). The experimental energy, depth, and transverse position resolutions of BGO block detectors were used as main inputs to the simulations to avoid extensive light transport in position encoding blocks. An improved model for energy resolution which includes the non-proportionality of the scintillation response of BGO and the statistical noise from photoelectron amplification in the PMT's was also proposed (Vozza et al 1997). Simulation studies have also been carried out to investigate the feasibility of using a triangular detection module for PET with neural networks to reconstruct the coordinates of the photon absorption point and thus recover the DOI information (Delorme et al 1996). Another exciting application is the use of a PET imaging system for monitoring the dose delivered by proton and gamma-ray radiotherapy beams (Litzenberg et al 1997). By measuring the amount and ratio of the beam-induced positron-emitting activity, the dose distribution and tissue composition may be determined.

4.3.2 Imaging systems and collimators design

      Monte Carlo techniques were extensively used to analyse the performance of new collimators design for planar gamma camera (Webb et al 1992), SPECT (Kimiaei et al 1997) and PET imaging (Digby et al 1990). Practical guidance could be offered for understanding trade-offs that must be considered for clinical imaging. Selective comparisons among different collimators could also be presented for illustrative and teaching purposes. In particular, Monte Carlo techniques were used to determine the effects of crystals with straight and pointed tips and septa on spatial resolution and efficiency (Thompson et al 1986), to compare the singles to true coincident events ratios in well collimated single, multi-slice and open collimator 3D configurations (Thompson 1988), to evaluate tungsten inter-plane septa of different thicknesses and geometries (Digby et al 1990) and to assess the effect of collimation on the scatter fraction (Thompson 1986).

      The design of SPECT (Bradshaw et al 1985) and PET (Keller and Lupton 1983) systems using the Monte Carlo method has received considerable attention and a large number of applications were the result of such investigations. The Monte Carlo method has also been used in the design of single-slice (Lupton and Keller 1983) and multi-slice PET scanners (Dahlbom et al 1989, Michel et al 1991). A VME bus-based microcomputer system has been used to implement a model for simulation of the flux of gamma rays in cylindrical PET detector systems (Stearns et al 1988). The program is capable of tracing over one million photons per hour and has been used to explore some of the effects of 'opening up' planar detector geometries into volumetric imagers. Rogers et al (1988) compared some of the performance parameters of a tomograph based on large area NaI(Tl) detectors to similar parameters of conventional small crystal machines. Michel et al (1991) used the GEANT package from CERN (Brun et al 1994) to study the response function and the scatter fraction in two PET scanners with and without interplane septa. The simulation of a large multiplane PET camera named HISPET (Del Guerra and Nelson 1988) and a planar imaging system made of two matrices, each one consisting of 400 (2x2x30 mm³) crystals of YAP:Ce using the EGS4 system have also been reported (Bollini et al 1997). Thompson (1990) investigated the effects of detector material and structure on PET spatial resolution and efficiency in terms of the number of interactions and tangential component of the mean square distance between the centroid and the point of first interaction.

      Several researchers used Monte Carlo simulation methods to study potential designs of dedicated small animal positron tomographs (Tzanakos and Pavlopoulos 1992, Pavlopoulos and Tzanakos 1996). An important conclusion drawn from these studies is that unlike human imaging where both sensitivity and spatial resolution limitations significantly effect the quantitative imaging performance of a tomograph, the imaging performance of dedicated animal tomographs is almost solely based upon its spatial resolution limitations (Miyaoka 1991). Recently, a conceptual design for a PET camera developed to image the human brain and small animals has been presented (Moses et al 1997). The authors performed a Monte Carlo simulation to predict the spatial resolution for a single plane PET camera with 3 mm LSO crystals. They concluded that the detector modules must be able to measure the DOI on an event by event basis in order to eliminate radial elongation artefacts, and that such depth information can be incorporated into the reconstruction algorithm in an artefact free way with a simple rebinning method.

4.3.3 Image reconstruction algorithms

      Monte Carlo simulations have been shown to be very useful for validation and comparative evaluation of image reconstruction techniques since it is possible to obtain a reference image to which reconstructed images should be compared. Fig. 4.5 shows transaxial slices of a simulated cylindrical phantom containing six cold spheres with diameters ranging from 9.5 to 31.8 mm reconstructed using the widely-used reprojection algorithm of Kinahan and Rogers (1989). According to visual inspection, the quality of the reconstructions is superior when the projections are obtained with good statistics. As a result of the enhanced quality, the small cold lesions are clearly visible when compared to those of the low count studies. This type of studies is useful to assess contrast recovery and to quantitatively predict the effect of different image reconstruction techniques on human accuracy for lesions detection.

      The ability to theoretically model the propagation of photon noise through emission computed tomography reconstruction algorithms is crucial in evaluating the reconstructed image quality as a function of parameters of the algorithm. Wilson et al (1994) used a Monte Carlo approach to study the noise properties of the ML-EM algorithm and to test the predictions of the theory. The ML-EM statistical properties were calculated from sample averages of a large number of images with different noise realisations. The agreement between the more exact form of the theoretical formulation and the Monte Carlo formulation was better than 10% in most cases examined, and for many situations the agreement was within the expected error of the Monte Carlo experiments. The same methodology was also followed to analyse a MAP-EM algorithm incorporating an independent gamma prior, and a one-step-late (OSL) version of a MAP-EM algorithm incorporating a multivariate Gaussian prior, for which familiar smoothing priors are special cases (Wang and Gindi 1997).

      

Fig. 4.5. : Effect of limited statistics on the quality of reconstructed images and lesion detection

      The Jaszczak phantom with spheres insert consists of a water-filled cylindrical phantom containing uniform activity in which 6 cold spheres with diameters ranging from 9.5 to 31.8 mm are included. The reference image (top left) and reconstructed slices generated with 1 Mcounts (top right), 10 Mcounts (bottom left), and 50 Mcounts (bottom right) are shown.

      The search for unified reconstruction algorithms led to the development of inverse Monte Carlo (IMC) reconstruction techniques (Floyd et al 1986). The authors used IMC to perform tomographic reconstruction for SPECT with simultaneous compensation for attenuation, scatter, and distance-dependent collimator resolution (see section 3.5.3).

4.3.4 Attenuation and scatter correction techniques

      The presence of scatter and attenuation in the images limits the accuracy of activity quantification in emission computed tomography (Zaidi 1996b). With no corrections, the uncertainty could be as high as 50-100%. Monte Carlo calculations have been found to be powerful tools to quantify and correct for photon attenuation and scattering in nuclear medicine imaging since the user has the ability to separate the detected photons into their components: primary events, scatter events, contribution of down-scatter events, .. etc. Monte Carlo modelling thus allows a detailed investigation of the spatial and energy distribution of Compton scatter which would be difficult to perform using present experimental techniques, even with very good energy resolution detectors (Lowry and Cooper 1987).

      Compton scattering effects in water on profiles of activity have been simulated (Logan and Bernstein 1983). The accuracy of experimental methodologies used for scatter fraction and scatter spatial distribution determination were evaluated using the Monte Carlo method (Acchiappati et al 1989). Barney et al (1991) developed an analytical simulation for single- and multiple-scattered annihilation photons in PET. The multiple-scatter model showed good agreement with Monte Carlo simulation of total object scatter. The authors also proposed a scatter correction method which uses the analytical simulation and exploits the inherent smoothness of the scatter distribution to account for three-dimensional effects in scatter distribution and object shape. Scatter components in PET divided into primaries, object scatter, gantry scatter and mixed scatter and their effects on the degradation of reconstructed images were also investigated (Adam et al 1996). Quantification of those components for a small animal PET prototype were also reported (Ziegler and Kuebler 1993). A Monte Carlo study of the acceptance to scattered events in a depth encoding large aperture camera made of position encoding blocks modified to have DOI resolution through a variation in the photopeak pulse height was performed (Moisan et al 1996). It was reported that the poorer discrimination of object scatters with depth sensitive blocks does not lead to a dramatic increase of the scatter fraction.

      Levin et al (1995) developed a correction method that uses the 3D reconstructed image volume as the source intensity distribution for a photon-tracking Monte Carlo simulation. The history of each annihilation photon's interactions in the scattering medium is followed, and the sinograms for the scattered and unscattered photon pairs are generated in a simulated 3D PET acquisition. The calculated scatter contribution is used to correct the original data set. Monte Carlo techniques were also used to estimate 'best possible' weighting functions for different energy-based scatter correction schemes and to examine the optimal number of energy windows for both NaI(Tl) and BGO scintillators (Haynor et al 1997).

4.3.5 Dosimetry and treatment planning

      There is no doubt that the area where early Monte Carlo calculations in the field have been performed is dosimetry modelling and computations (Andreo 1991). The approach adopted by the Medical Internal Radiation Dose (MIRD) committee of the Society of Nuclear Medicine was published as different pamphlets in a series of supplements to the Journal of Nuclear Medicine. Some of these pamphlets made extensive use of Monte Carlo calculations to derive specific absorbed fractions for electron and photon sources uniformly distributed in organs of mathematical phantoms. In particular, calculation of absorbed fractions for positron emitters relevant to neurologic studies were reported by Bice et al (1985). Interest in Monte Carlo-based dose calculations with b-emitters has been revived with the application of radiolabelled monoclonal antibodies to radioimmunotherapy (RIT).

      Currently, the preferred strategy is personalised human dosimetry with radiolabelled antibodies and this approach may become routinely employed clinically. The dose distribution pattern is often calculated by generalising a point source dose distribution (Akabani et al 1997), but direct calculation by Monte Carlo techniques is also frequently reported because it allows media of inhomogeneous density to be considered (Tagesson et al 1996, Furhang et al 1996). The development of a 3D treatment planner based on SPECT/PET imaging is an area of considerable research interest and several dose calculation algorithms have been developed. Fig. 4.6 lists the essential steps required in developing a 3D treatment planning program for RIT. Projection data acquired from an emission tomographic imaging system are processed to reconstruct transverse section images which yields a count density map of source regions in the body. This count density is converted to an activity map using the sensitivity derived from a calibration phantom. In the final step, this activity distribution is converted to a dose rate or dose map by convolving the activity distribution with dose point kernels or by direct Monte Carlo calculations. To elaborate a treatment plan for an individual patient, prospective dose estimates can be made by using a tracer activity of radiolabelled antibody to obtain biodistribution information prior to administration of a larger therapeutic activity. The clinical implementability of treatment planning algorithms will depend to a significant extent on the time required to generate absorbed dose estimates for a particular patient.

      

Fig 4.6. : Diagram showing the essential steps required in developing a three-dimensional internal dosimetry program based on quantitative positron emission tomography

4.3.6 Pharmacokinetic modelling

      A few applications of Monte Carlo techniques have been reported in the field of pharmacokinetic modelling. Casciari et al (1995) developed a compartmental model of [F-18] fluoromisonidazole transport and metabolism to compute the volume average kappa in tissue regions from [F-18] fluoromisonidazole PET time-activity data and characterised it using Monte Carlo simulations and PET time-activity data. This model was able to accurately determine kappa for a variety of computer generated time-activity curves including those for hypothetical heterogeneous and poorly perfused tissue regions. Compartmental models allow also the in vivo analysis of radioligand binding to receptor sites in the human brain. Benzodiazepine receptor binding was studied using a three-compartmental model (Millet et al 1996). The validity of the results of the coefficient of variation of each parameter were verified with statistical results provided by Monte Carlo simulation. A simplified approach involving linear-regression straight-line parameter fitting of dynamic scan data was developed for both specific and non-specific models (Thie et al 1997). Monte-Carlo simulations were used to evaluate parameter standard deviations, due to data noise, and much smaller noise-induced biases. The authors reported good agreement between regression and traditional methods.

5.1.1 Object-oriented programming concepts

      OOP is a method of programming using reusable components (called objects). An object is a software bundle of variables and related methods (Fig. 5.1(a)). Software objects are often used to model real world objects found in everyday life. Software objects interact and communicate with each other using messages. A class is a blueprint or prototype that defines the variables and the methods common to all objects of a certain kind. OOP is based on three major concepts: inheritance, encapsulation, and polymorphism. Inheritance is one way that object-oriented systems minimise duplication of code and effort. It allows the programmer to define a hierarchy of objects, where each object takes on the attributes and behaviours of its ancestors (Fig. 5.1(b)). This hierarchy is an inverted tree structure (with the root at the top and the leaves at the bottom) where each lower level of the tree defines more specific attributes about a particular class of object. For instance, the programmer can build a tree of windows, menus, or user objects (classes) with the top-most one implementing the most generic attributes and functions and the bottom-most ones implementing the most specific. With this structure, attributes and functionality can be defined at their highest level of reuse and be inherited by the rest of the objects below them. This process is usually called sub-classing, as the newly created class is a subordinate (or descendant) of another.

      

Fig. 5.1. : (a) Symbolic representation of an object: a group of related functions or methods and a data structure that serves those functions. (b) Inheritance hierarchy: each new class is the accumulation of all the class definitions in its inheritance chain

      In this example, class B inherit both from A, its superclass, and the root class. For instance, members of the B class will have methods and instance variables defined in all three classes B, A and root.

      The benefit of inheritance is that attributes and behaviours do not have to be separately defined for all objects in the system. This allows to increase productivity, promote consistency, reduce the chances for errors, and simplify system maintenance. Additionally, when standard functionality must change for a particular class of object, the programmer needs only to change it in one place and all descendant objects will automatically inherit the change with no additional coding. Encapsulation is a technique of combining data and code into an object class. For example, the programmer can define a function in a window, menu, or user object, and that function is then a part of that class. This encapsulates the data and the processing routines into one package. The main benefit of encapsulation is that it allows to hide information or shield complexity. When the functions that operate on an object are embedded in the object class itself, and only interact with that object through this defined interface, other system components are insulated from the details of these operations. Encapsulation has far-reaching implications in code maintenance and quality. When used properly, simple changes to systems will no longer cause other pieces of the system to break unexpectedly. Polymorphism is the ability of different objects to respond to the same message differently. It is a powerful concept and enables to define functions at a high level and let the details of the operation be handled by the object responsible for them. It also allows the programmer to exploit logical similarities between different classes of objects. Polymorphism is only possible with encapsulation, and it becomes even more powerful when combined with inheritance. Within an inheritance tree, any subclass of the class structure can modify or extend the way it behaves when it receives standard function calls from other system components without necessarily affecting those calling components.

      The major milestones for OOP were Simula, Smalltalk, and C11. Recently Java has become a popular OOP language (Rodgers 1996). There are, however, a number of other object-oriented languages with different strengths and focuses. The Smalltalk style of OOP quickly caught on in the Lisp community and led to Lisp extensions such as Flavors and CLOS. Eiffel and BETA are examples of statically typed languages. Self is an example of a prototype-based language. Objective C (NextStep 1992a) is a Smalltalk-like extension of C which can also be used as an extension to C++ when using the NextStep development environment (NextStep 1992b). At first glance, this may seem superfluous since C++ is itself an object-oriented extension of C. But C++ was designed primarily as a better C, and not necessarily as a full featured object-oriented language. It lacks some of the possibilities for object-oriented design that dynamic typing and dynamic binding bring to Objective C. At the same time, it has useful language features not found in Objective C. In addition to programming languages, there are a number of languages for analysis, design, and data modelling. For example, the unified modelling language integrates notations from several popular analysis and design methods. The main characteristics -objects and classes, message passing, abstraction, encapsulation, inheritance, and polymorphism- are the essence of OOP. Each contributes significantly to the overall utility of the paradigm, and each allows the paradigm to address one of the many software engineering concerns that motivated it. While different programming languages implement them in various combinations and to varying degrees, any language that implements them all is considered object-oriented.

5.1.2 Organising object-oriented projects

      A complete project management strategy includes strategies for creating, organising, and maintaining source files, building the application from its sources, running and debugging the application, revising the source files to fix bugs, and installing the finished application or preparing it for others to install. Object-oriented programming helps restructure the programming task in ways that benefit collaboration. It helps in eliminating the need to collaborate on low-level implementation details, while providing structures that facilitate collaboration at a higher level. When programs are designed at a high level of abstraction, the division of labour is more easily conceived. It can match the division of the program on logical lines; the way a project is organised can grow out of its design. The connections between the various components of an object-oriented program are worked out early in the design process. They can be well-defined, at least for the initial phase of development, before implementation begins. During implementation, only this interface needs to be coordinated, and most of that falls naturally out of the design. Since each class encapsulates its implementation and has its own name space, there is no need to coordinate implementation details.

      Object-oriented design (OOD) is the first step in the construction of an object-oriented system. OOP and object-oriented testing complete the process. OOD transforms the object-oriented analysis model of the real world into an implementation-specific model that can be realised in software. Unlike conventional software design methods, OOD results in a design that achieves a number of different levels of modularity. Conventional approaches to software design apply a distinct notation and set of heuristics to map the analysis model into a design model. OOD is different from traditional design in that the same concepts and terminology can be used all the way from analysis through design to implementation. The OOD process can be described as a pyramid composed of four layers. The foundation layer (subsystem design) contains a representation of each of the subsystems that enable the software to achieve its customer-defined requirements and to implement the technical infrastructure that supports customer requirements. The next layer up, (class and object) contains the class hierarchies that enable the system to be created using generalisations and increasingly more targeted specialisations. The third layer (message) contains the details that enable each object to communicate with its collaborators. This layer establishes the external and internal interfaces for the system. The top layer (responsibilities) contains the data structure and algorithmic design for all attributes and operations for each object. An additional foundation layer, which focuses on the design of domain objects, also exists. Domain objects play a key role in building the infrastructure for the object-oriented system by providing support for human-computer interface activities, task and data management.

      The three challenges of OOD are to identify the objects and classes needed to implement the software, describe the relationships between the identified objects and classes, and define the behaviour of the objects by describing the function implementations of each class. When the important objects have been found, the details of the class relationships and object behaviour are developed through a set of interconnected models that capture four dimensions of the software architecture: logical static, logical dynamic, physical static, and physical dynamic (Booch 1991). The logical static model describes the relationships of instantiation (class membership), association, aggregation, and generalisation. This is called the system object model. The generalisation relationship implies inheritance of attributes and methods. Graphical notations derived from the entity relationship notation for data modelling are frequently used to define the object model.

      There are a number of different OOD methods. Although each differs from its counterparts, all conform to the design pyramid and all approach the design process through two levels of abstractions: design of the system and subsystems, and design of individual objects. An example of an OOD method is the Jacobson method (Booch 1991) which comprises the following steps:

• Consider adaptations to make the idealised analysis model fit the real world environment.
• Create blocks as the primary design object.
• Create an interaction diagram that shows how stimuli are passed between blocks.
• Organise blocks into subsystems.
• Review the design work.

5.2.1 A simulation model for 3D PET

      The simulation model assumes a cylindrical array of detector crystals and known spatial distributions of annihilation sources and scatter phantoms. The positron histories and thus positron range are neglected during simulations with the assumption that the positron emission point is coincident with the annihilation point. Both the positron range and the non-collinearity effects are accounted for by analytic models and included with the other physical phenomena affecting spatial resolution. The annihilation point in the source region and the flight direction for the first of the two 511 keV annihilation photons are randomly generated. The second annihilation photon is generated in the same emission position, but with anti-parallel flight direction with respect to the first one. Pairs of annihilation photons are generated within the source objects in or out of the scanner's FOV; they are tracked until termination of their history either by interacting with scatter or detector objects, or by escaping the positron tomograph geometry and FOV. Fig. 5.2. illustrates the detection of these different types of coincidences. The Marsaglia algorithm (Marsaglia and Zaman 1993) was used to generate uniformly distributed pseudo-random numbers. This sequence of 24 bit pseudo-random numbers has a period of about 2¹⁴⁴ and has passed stringent statistical tests for randomness and independence.

      

Fig. 5.2. : Schematic of a simulation model for 3D positron tomography

      Detection of true, random, object scattered and detector scattered coincidences are indicated along the corresponding LORs.

      Dash lines indicate photon paths. Different steps are followed when tracing the photon in both the phantom and the detector volume. Both annihilation photons are tracked independently until a photoelectric interaction occurs in the detector volume or they escape from the geometrical acceptance of the imaging system. Photoelectric absorption as well as coherent and incoherent scatterings are taken into account to simulate photon interaction with scatter and detector objects. Interaction cross sections and scattering distributions are computed from parametrizations implemented in the GEANT simulation package (Brun et al 1994). The pathlength of the interacting photon is randomly generated according to the exponential attenuation based on the interaction length. The total cross section at the energy of the interacting photon determines the interaction length of the exponential distribution. The relative ratios of the cross sections for photoelectric effect, incoherent and coherent scattering to the total cross section are used to choose randomly which process occurs at each interaction vertex. A random number R uniformly distributed between 0 and 1 is sampled. The interaction type will then be (see section 2.1.3.5):

1. photoelectric absorption if R m £ m_photo;
2. incoherent scattering if m_photo £R µ < (m_photo+m_incoh);
3. else coherent scattering.

      In the case of photoelectric absorption, the total photon energy is transferred to an atomic electron and the random walk is terminated. The photoelectron propagation and the characteristic X-rays (with low energy with respect to 511 keV in the considered media) are not simulated since they have a negligible effect in the PET detection process. In an incoherent photon interaction, a fraction of the photon energy is transferred to the atomic electron. The direction of the scattered photon is changed to conserve the total momentum of the interaction. The Klein-Nishina expression for the differential cross section per electron for an incoherent interaction is used to sample the energy and the polar angle of the incoherently scattered photon. Fig. 5.3 shows energy and angular distribution of Compton-scattered photons having an initial energy of 511 keV in water as calculated by the Monte Carlo simulations. It illustrates clearly that small angle scattering is more likely than large angle scattering. The coherent scattering results only in a change in the direction of the scattered photon since the momentum change is transferred to the whole atom. In this latter case, the kinetic energy loss of the photon is negligible. The random number composition and rejection technique (Kahn 1956) is used to sample the momentum of the scattered photon and the scattering angle according to the form factor distributions. Coherent scatter distributions are sharply forward-peaked and vary considerably with atomic number and energy (Kaplan et al 1998).

      

Fig. 5.3.a : Energy distributions of Compton single-scattered photons having an initial energy of 511 keV in water as scattering medium

      

Fig. 5.3.b : Angular distributions of Compton single-scattered photons having an initial energy of 511 keV in water as scattering medium

      The coordinates after the photon travels a distance d are calculated using direction cosines. If the photon is deflected by a polar angle q and an azimuthal angle F, into a new direction (U_k+1, V_k+1, W_k+1) relative to the incoming direction (U_k, V_k, W_k), the coordinates after travelling the distance d are given by

      

(5.1)

      where

      

(5.2)

      In the case where W_k is close to unity, the calculations are then reduced to

      

(5.3)

      Fig. 5.4 shows the basic steps involved in tracing the photon in the phantom.

      

Fig. 5.4. : Flow diagram showing the steps followed when tracking the photon path in the phantom without applying variance reduction techniques

      A 3D PET scanner is simulated as a number of detection rings, each ring consisting of a number of scintillator crystals. The detection block, typical of current generation PET scanners is simulated by grouping crystals in matrices. Photon history within a crystal, across crystals within a block and across blocks is tracked. Crystals are considered as adjacent in the transaxial plane and in the axial direction. Two external shields are simulated as tungsten rings located at the two axial edges of the tomograph, partly shielding from radiation coming from outside the scanner FOV.

      In the detector blocks, at each interaction vertex, the energy released is recorded. Tracking is stopped either by a photoelectric absorption, escape of the photon from the block volume, or by a Compton scattering leaving less than 5 keV to the scattered photon. The energies of all interaction vertices are summed to yield the total absorbed energy in the scintillation detector. This total energy is assumed to be converted to scintillation light using a Gaussian random smearing to account for the combined energy resolution, DE/E , of the scintillator and its photomultipliers. Photon fractional loss in not taken into account. A knowledge of the energy resolution of the detector block is required since the FWHM must be known for a specific energy, e.g. at 511 keV. The energy resolution is then a function of 1/ .

      

(5.4)

      The energy resolution of the detector blocks depends on the PET tomograph being simulated. For example, for the ECAT 953B scanner, it was estimated as 23% at 511 keV (Spinks et al 1992). The position blurring step then calculates the mean detection coordinates (X,Y,Z) of each incident photon. This is done by computing the centroid of all interaction vertices, each weighted by the ratio of its individual energy to the total energy. A lower cut-off energy is set to model the energy acceptance window used in PET scanners to reject scattered radiation. When a photon falls below this user selectable energy value, its history is terminated. The mean X and Y coordinates of each photon are smeared to account for spatial resolution and position miscoding in the detector block.

      The spatial resolution in PET is a three-dimensional parameter. For simulation purposes, we need to understand the effect of the spatial resolution limitation projected in any given direction. There are different independent components that contribute to the overall spatial resolution. It is assumed that these blurring factors are isotropic in both transaxial and axial directions. The contributions to spatial resolution on measured projection data near the centre of a tomograph with adequate sampling are the crystal size, detector penetration (depth of interaction) and inter-crystal scatter, crystal decoding process, annihilation photon accolinearity, and positron range, and finally source size and motion. Therefore, the measured (intrinsic) resolution is a convolution of the various resolution components. If the different resolution components are assumed to be Gaussian in shape and are described by FWHM, then the combined resolution is the squared sum of the individual resolution components. All factors related to detector block effect are combined in a single component, s_block. For each detected coincidence, a convolution with a Gaussian having an appropriate width (s_gauss) is performed:

      

(5.5)

      where:

      s_range: accounts for the loss of resolution due to the positron range;

      s_non-col: accounts for the loss of resolution due to the annihilation photons non-collinearity effect. This value is estimated in the centre, where this effect is maximum, and assumed constant over the whole FOV;

      s_block: accounts for the loss of resolution due to the block blurring effect.

      s_block is approximated by the following formula (Moses et al 1998), where d is the detector size and b is a factor due to the crystal decoding process (0 if the crystal of interaction is individually coupled to a photodetector and approximately 2.2 mm otherwise). This factor while empirically determined for BGO detector modules, has been found to be necessary to 'predict' the performance of several new high-resolution PET cameras (Moses et al 1998). The same paper suggests an empirical formula for calculating the linear fluctuations due to photons non-collinearity in any given direction which can be modelled as a Gaussian distribution (around 180° with an FWHM of about 0.3°) centred at the origin and expressed as s_non-col = 0.0022 x D_s, where D_s is the scanner's diameter (in mm). The intrinsic spatial resolution calculation for the ECAT 953B scanner when using ¹⁸F is estimated as 3.96 mm according to:

      s_range = 0.5 mm for ¹⁸F;

      s_non-col = 0022 x 760 = 1.67 mm;

      s_block = 3.56 mm.

      

Fig. 5.5. : Flow diagram showing the main steps followed when tracking the photon history in the scintillation crystals volume

      The incident energy E_i , the total imparted energy E_n, and the imparted energy at each interaction site in the crystal E' are shown.

      The transport of light produced during scintillation in the crystals and the block-associated electronic chain (photomultipliers, pre-amplifiers and amplifiers) were not simulated. Random coincidences are not simulated but could be analytically modelled after simulation using the singles rate. Interaction within scatter or detector objects can be switched on and off. In case interaction within detector objects is switched off, any photon impinging on a detector is assumed to deposit all its energy in the detector crystal. When it is switched on, photon pairs are recorded once two photons resulting from one annihilation event have passed the energy window set for discrimination. Radial samples are assumed to be equidistant, although ring curvature can be taken into account for sampling. Fig. 5.5 shows the different steps followed when tracing the photon in the scintillation crystals volume.

      To find the number of detector combinations to cover a transaxial FOV of radius R_fov, the fan-beam angle a is determined, with which the cylindrical FOV will project itself onto the periphery of the detector ring seen from the side of one of the detectors of the ring (Fig. 5.6). Thus, the fraction of detectors which a given detector must be able to record coincidences within the FOV is:

      

(5.6)

      where

      

(5.7)

      

Fig. 5.6. : Illustration of the scanner geometry and relationship between the tomograph's ring radius R_s, the field-of-view radius R_fov and the size of the fan-beam angle a

      If the number of detectors per detector ring is denoted by N_d, the number of actual detectors, N to cover the FOV is:

      

(5.8)

      The number of coincidence combinations per detector ring is then

      

(5.9)

5.2.2 Software design

      Unlike procedural programming languages which separate data from operations on data defined by procedures and functions, object-oriented programming languages consist of a collection of interacting high-level units, the objects, that combine both data and operations on data. This renders objects not much different from ordinary physical objects. This resemblance to real things gives objects much of their power and appeal. They can not only model components of real systems, but equally as well fulfil assigned roles as components in software systems. This programming paradigm appears to be quite well appropriate to Monte Carlo modelling as there is a direct correspondence between objects of the application domain - source and scatter volumes, detectors - and objects of the computational model - source, scatter and detector objects -.

      The first version of the Monte Carlo simulator, Eidolon, was written in Objective-C (NextStep 1992a), an object-oriented programming language based on ANSI C, and run under the NextStep 3.3 development environment (NextStep 1992b) on an Hewlett-Packard 712/60 workstation with 48 Mb RAM memory. In order to ease the job of incrementally adding capabilities to the Monte Carlo simulator, a modular design featuring dynamically loadable program elements or bundles was adopted. The basic building block is a model element object class which allows elements to be browsed, inspected, adjusted, created and destroyed through a graphical inspector. This was then used to implement simple parametric source, detector and scatter object classes and sinogram and image object classes to view and save the generated data. A controller object oversees the simulation process (Fig. 5.7).

      

Fig. 5.7. : A flow diagram showing the main components of the Monte Carlo simulator

      Eidolon, like any object-oriented program, uses objects built from a variety of classes. In Objective-C, objects are defined by their class. Another characteristic of Objective-C is that many decisions are postponed from compile time to run time. Therefore, the Objective-C language depends on a run-time system to execute compiled code. The application is organised in a convenient way by dividing the program into different compartments called subprojects. Each subproject could be compiled and debugged separately. The software is structured as follows: the graphical interface (not used on Parsytec-CC system), objective-C classes, categories and protocols, and procedural programming C routines. The graphical interface as well as the link between graphical objects and objective-C code was designed using the NextStep Interface Builder (nib) development tool.

      Class definitions are additive; each new class is based on another class through which it inherits methods and instance variables. The new class simply adds to or modifies what it inherits. Inheritance links all classes in a hierarchical tree with a single class, the Object class, as its root. Eidolon takes advantage of this property by organising a hierarchical tree such that methods used by different objects need not be implemented twice. Basic classes are used to handle gamma pair creation, photon interaction and detection methods as well as other matrix handling routines, 3D shape definitions and geometrical calculations. Categories are extensively used to add methods to an existing class. This simplifies the management of large classes and allows to beneficiate of incremental compilation. Protocols are also used to capture similarities among classes that are not hierarchically related. A protocol is a list of method declarations, not associated with a particular class definition. Thus, any class or many classes can implement them by adopting the protocol.

      In the version implemented under NextStep environment, Eidolon exploits the dynamic loading property of objective-C, thus classes and categories are loaded while the program is running. The new code is incorporated into the program and treated identically to classes and categories loaded at the beginning. For example, objects adjustable by the user through the graphical interface for a particular simulation study, such as source, scatterer or scanner parameters are dynamically loaded into the kernel at run-time. These are referred to as dynamically loadable modules or bundles. Unfortunately, dynamic loading is not yet supported on gcc/PowerPC-AIX4.1.4. Thus bundles are handled as the other classes in the version running on the Parsytec system.

5.2.3 Graphical user interface

      The polymorphism of object-oriented programs yields simpler programming interfaces, since the same names and conventions can be reused in any number of different classes. The result is less to learn, a greater shared understanding of how the whole system works, and a simpler path to cooperation and collaboration. A user interface must meet the needs of both novice and experienced users:

• For the novice or infrequent user, it must be simple and easy both to learn and to remember. It shouldn't require any relearning after an extended absence from the computer.
• For the more experienced user, it must be fast and efficient. Nothing in the user interface should get in the way or divert the user's attention from the task at hand.

      The challenge is to accommodate both these goals in ways that don't conflict to combine simplicity with efficiency. A graphical (window-based) user interface is well suited to this task. Because graphical objects can be endowed with recognisable features of real objects, users can borrow on their everyday experience when they approach the computer. Graphical buttons work like we would expect real buttons to work, windows behave much like separate tablets or sheets of paper, sliders and other graphical objects act like their physical counterparts off-screen. The computer becomes less an entirely new world with its own rules, and more an extension of the more familiar world away from the computer screen. This not only makes the user interface easier to learn and remember, it also permits operations to be simpler and more straightforward. Picking an option is as easy as flicking a switch. Resizing a window is as direct and simple as pulling on a tab. The same attributes of the user interface that provide simplicity for novice users can also result in efficiency for more expert users.

      

      The reference image (top right) and the sinograms (bottom right) may be viewed as they are generated and are periodically updated.

      The graphical user interface designed for Eidolon allows the user to select scanner parameters such as the number of detector rings, detector material and sizes, energy discrimination thresholds and detector energy resolution (Fig. 5.8). It also allows one to choose a set of simple 3D shapes, such as parallelepiped, ellipsoid or cylindroid for both the annihilation sources and the scattering media, as well as their respective activity concentrations and chemical compositions. One may view the reference source image and the resulting simulated sinograms, as they are generated. The reference image and sinogram displays are periodically updated.

5.2.4 Monte Carlo code features

      The Monte Carlo simulator allows to simulate various configurations of volume-imaging PET scanners and to obtain detailed information on the different processes occurring within the phantom and detectors. For example, energy pulse-height distribution, point-spread function (PSF) and the scatter fraction could be obtained. It also allows to separate the detected photons into their components: primary events, scatter events, contribution of down-scatter events, ..etc., and to perform detailed investigations of the spatial and energy distribution of Compton scatter which would be difficult to perform using present experimental techniques. The PSF is useful for evaluating the characteristics of an imaging system. Fig. 5.9 shows the simulated PSF for the ECAT 953B scanner which clearly shows that either scattered and primary events displaced relative to the position of the origin of the emission site will be acquired. The effect of this is the broadening of the PSF tails.

      

Fig. 5.9. : A simulated point-spread function is shown together with the scattered part for a point source on axis

      For a detector system with ideal resolution, photons will be detected as orthogonal to the crystal surface (solid line). In reality, primary photons will interact with the detector material after travelling a certain distance resulting in position blurring and degradation of the spatial resolution (dashed line). Due to the finite energy resolution, scattered photons will also be detected (solid line). The resulting total PSF (primary + scatter) is also shown (dotted line).

      The simulator is an efficient tool that can be used to generate data sets in a controllable manner in order to assess different reconstruction algorithms. As the 'best' algorithm can only be selected with respect to a certain task, different 'basic' performance measures can be used. Image degrading effects are illustrated using simulated projections of the digitised 3D Hoffman brain phantom (Hoffman et al 1990). A slice of this phantom is shown in Fig. 5.10. The ratio between the activity in white, grey matter and ventricles was chosen as 1:4:0, respectively.

      The projections of this phantom at different levels of fidelity are generated. The strengths of the image degrading factors are characteristic of a ¹⁸F-FDG brain study and an energy window of [380-850] keV. Fig. 5.11.E-H shows the effects of different aspects of image degradation on FBP reconstructions. The loss of resolution caused by detector blurring (FWHM = 4 mm) on projection data and FBP reconstructions is shown in Fig. 5.11.B and 5.11.F.

      

Fig. 5.10. : Transaxial slice of the digital 3D Hoffman brain phantom

      

Fig. 5.11. : Simulation of the 3D Hoffman brain phantom in different imaging situations showing 2D projections (top row) and their filtered backprojection reconstructions with the reprojection algorithm (bottom row)

      In case A and E, a 2D projection and reconstructed image neglecting attenuation, scatter and detector blurring. In B and F, effects of detector blurring are included, while in C and G, effects of detector blurring, attenuation and scatter are included in the simulation. Finally, in D and H, effects of detector blurring, attenuation and scatter are included and appropriate corrections for attenuation and scatter applied.

      The scatter components in 3D PET can be treated as the sum of object and detector scatters. The object scatter consists in photons which have interacted within the object, mainly by Compton scattering (99.7% at 511 keV for H₂O). The detector scatter consists in photons which have interacted within detector blocks by Compton scattering (53% at 511 keV for BGO) or by Rayleigh scattering (6% at 511 keV for BGO).

      In a real measurement, scatter components in the detectors and in the FOV are indistinguishable from other secondary effects. Monte Carlo simulation of energy distributions provide insight into the scattering processes arising in 3D positron tomography. The detected events can thus be distinguished and separated in an efficient way into scattered and unscattered coincidences. The sinograms resulting from the simulation of the Hoffman phantom have been separated into unscattered and scattered coincidences and are illustrated in Fig. 5.12.

      

Fig. 5.12. : Representative simulated sinograms of the Hoffman brain phantom taking into account the effects of detector blurring, attenuation and scatter showing only unscattered events (left), scattered events (middle) and total detected events (right)

      Fig 5.13 shows energy pulse-height distribution obtained by simulation of a line source in the centre of a water-filled cylindrical phantom. The scattered events in the energy pulse-height distribution have been separated according to the order of scattering (coherent and incoherent in the phantom). It is clear from viewing Fig. 5.13 that events from some scattered photons will not be rejected by the [380-850 keV] energy discrimination due to the limited energy resolution.

      

Fig. 5.13. : The energy distribution due to scattered photons resulting from the simulation of a line source placed in the centre of a 20 cm diameter water-filled cylinder is separated into different contributions (total scattering or different orders of photon scattering)

      Energy resolution is proportional to the inverse square root of the deposited energy and is simulated by convolving the deposited energy with a Gaussian function whose FWHM is 23% for 511 keV photons.

      The simulation package was used to obtain information on the different processes occurring within the phantom. Fig. 5.14 shows the distribution of the number of scattering interactions of annihilation photons within the phantom. Approximately, 32% of the events undergo one or more scattering interactions before they exit the phantom, with the mean number of interactions being equal to 0.46. As shown, a significant portion of the photons have interacted and lost energy in the phantom, in particular a photon that is produced near the rim of the phantom has a higher probability of escaping without interaction.

      

Fig. 5.14. : Distribution of the number of scattering interactions within the phantom

      Approximately 32% of the annihilation photons undergo one or more scattering interactions within the phantom for this particular geometry and energy window setting.

      The detector scatter is characteristic of the detection system and depends on the energy discrimination threshold. Due to the high density of detector materials, detector scatter contribution has a narrow distribution around the detector of first interaction. For this reason, the contribution of detector scatter has negligible effect on the overall response function and has been ignored in the determination of the scatter fraction. Fig. 5.15 shows the distribution of the number of scattering interactions within the BGO crystals. As shown, about 50% of the incident photons undergo one or more scattering interactions. Fig. 5.16 shows the distribution of the distance travelled in the BGO crystal per interaction. The energy dependence of the attenuation coefficient as well as the geometry of the crystal must be taken into account in understanding this distribution. It should be noticed here that the shape of the distribution is determined by the energy of the photons (since the attenuation coefficient is a function of the energy) and the geometry of the detection system.

      

Fig. 5.15. : Distribution of the number of scattering interactions within the detector blocks

      Approximately 55% of the annihilation photons undergo one or more scattering interactions within the phantom.

      

Fig. 5.16. : Distribution of the distance travelled by annihilation photons inside the crystal between interactions

      The mean distance is about 1.04 cm.

      The lower level threshold (LLT) can be easily changed and its effect on the scatter components studied in an effective way. Fig. 5.17 shows transaxial profiles in the projection space for a simulated line source in a 20 cm diameter water-filled cylinder as a function of the LLT. In each case, the cylinder is centred on the origin of the graph. It is clear that the boundary (at -10 cm, +10 cm) of the object has no influence on the profile. However, the LLT set-up greatly influences the amount of scatter in the projection data.

      

Fig. 5.17. : Sum of one-dimensional transaxial projections resulting from the simulation of a line source placed in a 20 cm diameter cylinder filled with water as a function of the lower energy discrimination threshold

      This illustrates the compromise that should be attained between increasing the lower level threshold to reduce scattered events and the variance in the reconstructed images resulting from limited statistics.

5.2.5 Validation of the simulation model

      In order to assess the accuracy in modelling the physical response of PET tomographs in 3D mode, physical parameters (intrinsic spatial resolution, energy spectra and scatter fraction), representative of the performance of PET scanners were estimated to verify the simulated distribution against phantom measurements and published data. We have chosen to compare our simulations with measurements (Spinks et al 1992) and published data from the Triumph PET group (Moisan et al 1996). Some measurements were performed on the ECAT 953B (16 rings of 384 detectors each with a ring radius of 38 cm) and others on the GE-Advance (18 rings of 672 detectors each with a ring radius of 46.35 cm) PET scanners operated in 3D mode. It should be pointed out that our simulation does not include backscatter from equipment in the laboratory and the housing of the detector blocks. Nor does it simulate the natural background radiation always present in the scanning rooms. It is therefore expected that the shape of the measured and simulated energy spectra will differ in the low-energy region.

---

---

---

---

---

5.3.1 Object modelling

      Object modelling is fundamental for performing photon and electron transport efficiently by means of a Monte Carlo method. It consists of a description of the geometry and material characteristics for an object (Ogawa et al 1997). The material characteristics of interest include density and energy-dependent cross sections. The modelling includes simple geometry (SG), shape-based (SB), and voxel-based (VB) approaches. The three approaches use a piecewise uniform distribution of object characteristics to model an object. With the SG model, an object is composed of a simple combination of primitives such as cylinders and spheres. The Utah phantom (Fig. 5.29), which was designed with a high degree of inhomogeneity both transaxially and axially in order to compare and test scatter correction techniques in 3D PET can be accurately represented using this approach. The outer compartment is generally used to provide activity from outside the field-of-view. The SB approach represents the boundaries of shapes by mathematical equations. Regular shapes such as sphere, cylinder, rectangular solid, etc. have been used to approximate irregularly-shaped regions. The VB approach discretises an object into tiny cubes (voxels) with uniform characteristics. An object is thus represented by a union of voxels of the same size.

      

Fig. 5.29. : Diagram of the Utah phantom

      Extensions of SG and SB models such as the solid geometry-based (SGB) approach (Wang et al 1992) includes more primitives (ellipsoids, elliptic cylinders, tapered elliptic cylinders, rectangular solids, and their subsets: half, quarter, and eighth) and uses an inclusion tree data structure to provide relationships between primitives. These extensions provide simple irregular shape modelling. To allow anthropomorphic modelling the composite model (Wang et al 1993), which is an extension to the SGB approach adds to the primitives a voxelised rectangular solid primitive. Combinatorial approaches to solid modelling, which describe complex structures as set-theoretic combinations of simple objects, are limited in their ease of use and place unrealistic constraints on the geometric relations between objects such as excluding common boundaries. An approach to volume-based solid modelling has been developed, which is based upon topologically consistent definitions of boundary, interior, and exterior of a region (Li and Williamson 1992). From these definitions, union, intersection, and difference routines have been developed that allow involuted and deeply nested structures to be described as set-theoretic combinations of ellipsoids, elliptic cylinders, prisms, cones, and planes that accommodate shared boundaries.

      An octree-based method (OCT) which describes an object by using several sizes of cubic regions was proposed by Ogawa and Maeda (1995) to increase the calculation speed in photon transport since the number of voxels is much smaller than that of the VB approach. The 'octree string' is generated from a set of serial cross sections automatically. The same author developed a modelling method called the maximum rectangular region (MRR) method (Ogawa et al 1997). In this approach, a MRR for a given voxel is selected within a homogeneous, irregularly shaped region from a set of cross sections. The search is performed by checking the six sides of a box (MRR) including the voxel of interest. With the MRR representation of the object, high speed calculation of photon transport can be accomplished because an object can be described by means of fewer regions than in the VB or the OCT representation methods. Fig. 5.30 illustrates the calculation time required for the VB, OCT and MRR approaches for different image matrix sizes.

      

Fig. 5.30. : Calculation times for the VB, OCT, and MRR representation of an object

(Ogawa et al 1997)

5.3.2 Anthropomorphic phantoms

      Modelling of imaging and other medical applications is best done with phantom models that match the gross parameters of an individual patient. Computerised anthropomorphic phantoms can either be defined by mathematical (analytical) functions, or digital volume arrays. The mathematical specifications for phantoms that are available assume a specific age, height, and weight. People, however, exhibit a variety of shapes and sizes. In the first MIRD pamphlets, several organs including the skeletal system were represented schematically using geometric forms (cylinders, cones and ellipsoids) (Snyder et al 1978). The representation of internal organs with this mathematical phantom is very crude since the simple equations can only capture the most general description of the organ's position and geometry. A version of this phantom has been updated to include female organs.

      Mathematical phantoms are still evolving and are being constantly improved. The heterogeneity of the body has been taken into account by including soft tissues, bone and lungs with different compositions and densities. For certain organs such as the stomach and the bladder, a distinction should be made between the organ contents and the organ wall. The Mathematical CArdiac Torso (MCAT) phantom is an anthropomorphic phantom, developed at the University of North Carolina at Chapel Hill, that has been used in emission computed tomography imaging research (LaCroix 1997). Using mathematical formulae, the size, shape and configurations of the major thoracic structures and organs such as the heart, liver, breasts, and rib cage are realistically modelled for imaging purposes. Though anatomically less realistic than phantoms derived from CT or MR images of patients, the MCAT phantom has the advantage that it can be easily modified to simulate a wide variety of patient anatomies. In addition, the MCAT phantom simulates a dynamic, beating heart including changes in myocardial wall thickness, changes in chamber volumes, apical movement and heart rotation during the cardiac cycle. The phantom consists of two physical models: a 3D distribution of attenuation coefficients and a 3D distribution of radionuclide uptake for the various thoracic organs. The 3D attenuation coefficient phantom classifies all thoracic tissues into one of 5 types: muscle (vasculature and other soft tissues), lung, fat (such as in the breasts), trabecular bone and cortical bone. The MCAT phantom has become a valuable tool in imaging studies where reasonably realistic, but anatomically variable patient data needs to be simulated. The graphic in Fig. 5.31 illustrates the MCAT phantom. The graphic is a volume-rendering of anterior and posterior views with some sections removed for visualisation purposes.

      

Fig. 5.31. : Surface rendered images of the 3D MCAT phantom developed at Chapel Hill (a) Anterior view with outer body surface and ribs removed to show the various organs modelled. (b) posterior view with the rib cage present

(LaCroix 1997)

      Some calculations make use of more accurate representations of individuals based on volumetric scans, such as CT, MRI, and PET. As an improvement to the mathematical anthropomorphic phantoms, a new family of phantoms was constructed from CT data. The human phantoms present advantages towards the location and shape of the organs, in particular the hard bone and bone marrow. A physical brain phantom has also been developed to simulate the activity distributions found in the human brain in the cerebral blood flow and metabolism studies currently employed in PET (Hoffman et al 1990). The phantom utilises thin layers of Lucite to provide apparent relative concentrations of 4, 1 and 0 for grey matter, white matter and ventricles, respectively, in the brain. A clinically realistic source distribution simulating brain imaging was thus created in digital format. Zubal et al (1994) developed a typical anthropomorphic VB adult phantom by manual segmentation of CT transverse slices of a living human male performed by medical experts. A computerised 3D volume array modelling all major internal structures of the body was then created. Each voxel of the volume contains an index number designating it as belonging to a given organ or internal structure. These indexes can then be used to assign a value, corresponding to, e.g. density or activity. Two versions of the phantom exist, representing either the complete human torso with an isotropic voxel resolution of 1.5 mm, or a dedicated head phantom with 0.5 mm voxel size. The dedicated brain phantom was created from the high resolution MRI scans of a human volunteer. This volume array represents a high resolution model of the human anatomy and serves as a VB anthropomorphic phantom.

      The VB model allows to envision some calculations requiring the use of more accurate representations of individuals based on volumetric scans, but obviously at the expense of increased computing time. Ray-tracing approaches based on an implementation of Siddon's algorithm are used to compute the length of intersection of each ray with every voxel. By considering pixel boundaries as the intersection of orthogonal sets of equally spaced planes, rather than independently, the computational complexity is greatly reduced. The points of intersection of a line with a set of parallel planes are particularly simple to determine (Fig. 5.32). Once one such point has been found, all the others may be found by recursion. A detailed description of the algorithm may be found elsewhere (Siddon 1985).

      

Fig. 5.32. : Principle of the ray-tracing algorithm

      The lengths of intersection of a line with an array of pixels may be computed independently for every pixel (left), or, more efficiently, by considering pixel boundaries as the intersection of orthogonal sets of equally spaced planes (middle and right).

      The current version of Eidolon is capable of modelling photon transport through media in which the material properties vary in three dimensions. Both SB and VB models are supported within the simulator. A detailed description of the way to specify the source and scattering medium for both models is given elsewhere (Zaidi et al 1997b).

6.1.1 Parallel computing platforms

      The theoretical peak performance of a computer is determined by counting the number of floating-point additions and multiplications that can be completed during period of time, usually the cycle time of the machine (Dongarra 1998). Today's large machines measure their speed in GFlops (10⁹ operations/s). Each CPU generally contains 2 multiply and 2 add pipes. When all of these can be employed simultaneously as for instance in a dot product or a vector update operation, 4 Flops/cycle can be attained. The Cray T90 has a cycle time of 2.2 ns and a maximum number of 32 processors, thus the peak performance is :

      .

      Easily portable Monte Carlo user codes are generally used for timing benchmark purposes on different computers. According to Van der Steen and Dongarra (1998), the classification of high-performance computers is based on the way instructions and data streams are arranged and comprises four main architectural classes.

---

---

---

---

6.1.2 Parallelisation strategies for Monte Carlo codes

      Sequential programs make the most effective use of the available processing power: they alone guarantee maximum use of the CPU. In parallel programs, communication management introduces an unavoidable overhead, resulting in less efficient use of the overall CPU power. Moreover, according to Amdahl's law (Amdahl 1967), parallelisation efficiency is decreased by a factor representing the fraction of operations that must be executed in sequential order. When this fraction reaches one, we are confronted with a wholly unparallelisable code, and the speed-up is zero no matter how many processors are used. The efficiency of parallel programs is furthermore reduced by a factor equal to the fraction of processor idle time, which is highly dependent on the software parallelisation techniques used by the programmer.

      In photon transport simulation, scalar or serial Monte Carlo codes track the history of one particle at a time, and the total calculation time is the sum of the time consumed in each particle history. Many Monte Carlo applications have characteristics that make them easy to map onto computers having multiple processors. Some of these parallel implementations require little or no inter-processor communication, and are typically easy to code on a parallel computer. Others require frequent communication and synchronisation among processors and in general are more difficult to write and debug. A common way to parallelise Monte Carlo is to put identical "clones" on the various processors; only the random numbers are different. It is therefore important for the sequences on the different processors to be uncorrelated so that each processor does not end up simulating the same data (De Matteis and Pagnutti 1995). That is, given an initial segment of the sequence on one process, and the random number sequences on other processes, we should not be able to predict the next element of the sequence on the first process. For example, it should not happen that if random numbers of large magnitude are obtained on one process, large numbers are more likely to be obtained on another. Furthermore, in developing any parallel Monte Carlo code, it is important to be able to reproduce exactly Monte Carlo runs in order to trace program execution.

      The basic principles of parallel and vector processing are illustrated in Fig. 6.1. In history-based parallel processing, each particle (p₁, p₂, .., p_m) is assigned to one process which tracks its complete factual history (e₁, e₂, .., e_n). In event-based vector processing, a process treats only part of each particle history (e₁, e₂, .., e_n) and particles (p₁, p₂, .., p_m) "flow" from process to process. The first approach seems to be best suited to the physics of the problem. By spreading out the work among many processors, a speed-up that approaches the number of processors being used could be attained. This can be achieved through the use of parallel processing environments including arrays of transputers (Askew et al 1988, Ma 1994), vector parallel supercomputers (Miura 1987, Smith et al 1993), massively parallel computers (Miura and Babb 1987), or a cluster of workstations in a local area network using a parallel computing simulator such as PVM (Kirkby and Delpy 1997). While vector processors were designed for single instructions operating on multiple data (SIMD), parallel computers are for multiple instructions operating independently on multiple data (MIMD). In this later case, memories are distributed on each processor and the communications between processors are carried out by passing messages. Processors in MIMD machines are subject to asynchronous and unbalanced external effects and are thus, for all practical purposes, impossible to keep aligned in time in a predictable way. If the assumption is made that random number generation from a single generator will occur across processors in a certain predictable order, then that assumption will quite likely be wrong. A number of techniques have been developed that guarantee reproducibility in multiprocessor settings and with various types of Monte Carlo problems (Bhavsar and Isaac 1987). Although PVM (Geist et al 1994) or MPI (Snir et al 1996) can be used to provide a machine independent message passing interface, uniform assignment of particles to processes is not suited for workstation clusters because all the processors in the cluster may not have the same computation power in a heterogeneous environment, and the processing speed is bound by the slowest process. Moreover, the programmer needs to handle the different formats (big-endian and little-endian) computers use to keep data in memory.

      

Fig. 6.1. : Comparison between history-based scalar processing, event-based vector processing and history-based parallel processing

      In history-based scalar processing, one particle history is tracked at a time. In history-based parallel processing, each particle (p₁, p₂, .., p_m) is assigned to one process which tracks its complete history (e₁, e₂, .., e_n). In event-based vector processing, a process treats only part of each particle history (e₁, e₂, .., e_n) and particles (p₁, p₂, .., p_m) 'flow' from process to process.

      During the last decade, investigations were carried out to run the EGS4 photon-electron Monte Carlo transport code on vector machines (Brown and Martin 1984, Martin and Brown 1987) and a speed up of about 8 was reported with the vectorized code (Miura 1987). The same code was also implemented on a multiple-transputer system (Ma 1994) and a parallel computer (Miura and Babb 1987). Different approaches for parallel execution of this code including large-grain data flow techniques were reported (Babb and Storc 1988). Last but not least, Smith et al (1993) described a vectorized code for simulation of SPECT imaging that uses an event-based algorithm in which photon history computations are performed within DO loops. The indices of the DO loops range over the number of photon histories, and take advantage of the vector processing unit of the Stellar GS1000 computer for pipelined computations.

6.1.3 The Parsytec-CC system

      The Parsytec CC system* ²  is an autonomous unit at the card rack level. The CC card rack subsystem provides the system with its infrastructure including power supply and cooling. The system is configured as a standard 19" rack mountable unit which accepts the various 6U plug-in modules. The Parsytec CC system is a distributed memory, message passing parallel computer and is globally classified into the MIMD category of parallel computers. Its architecture is based on the mainstream Motorola MPC 604 processor running at 133 MHz with 512KB L2-cache. The modules are connected together at 1 Gbits/s with high speed (HS) link technology according to the IEEE 1355 standard, allowing data transfer at up to 75 Mbytes/s. The communication controller is integrated in the processor nodes through the PCI bus. PCI is also the standard which is used in the I/O modules and customer-specific I/O functions. All the nodes are directly connected to the same router which implements an active hardware 8 by 8 crossbar switch for up to 8 connections using the HS link. A schematic representation of the CC node is given in Fig. 6.2. It combines a memory controller (MPC 104) with a MPC 604 PowerPC processor and its local memory through an on-board PCI bus. The system board uses the MPC 105 chip to provide memory control, DRAM refresh and memory decoding for banks of DRAM and/or Flash. The CPU bus speed is limited to 66 MHz while the PCI bus speed is 33 MHz at maximum.

      

Fig. 6.2. : Schematic layout of the Parsytec CC node, showing the main PowerPC 604 CPU, its local memory through an on-board PCI bus, and the MPC 105 chip for memory control

      The software is based on IBM's AIX 4.1 UNIX operating system together with Parsytec's parallel programming environment Embedded PARIX (EPX 1996). Thus, it combines a standard UNIX environment (compilers, tools, libraries) with an advanced software programming development environment. The system was integrated to the local area network using standard Ethernet. Currently a CC node has a peak performance of 266 MFlops. The peak performance of the 8-node CC system installed at Geneva University Hospital is therefore 2.1 GFlops. A schematic view is shown in Fig. 6.3. To develop parallel applications using EPX, data streams and function tasks are allocated to a network of nodes. The data handling between processors requires just a few system calls. Standard routines for synchronous communication such as send and receive are available as well as asynchronous system calls. The full set of EPX calls establishes the EPX application programming interface (API). The destination for any message transfer is defined through a virtual channel that ends at any user defined process. Virtual channels are user defined and managed by EPX. The actual message delivery system software utilises the router.

      

Fig. 6.3. : Schematic view of the Parsytec CC8 system installed at Geneva University Hospital and integrated to the local area network (LAN)

      It consists mainly of a router, a networked entry node, one I/O node and six computing nodes housed in the same 19" rack.

6.1.4. Porting the software on the parallel platform

      Although, the NextStep development environment provided a great tool allowing to reduce the amount of time necessary for writing applications, the portability issue still remained hard to solve since these applications were tied to machines running NextStep. Fortunately, the GNU C compiler (gcc) available from the Free Software Foundation (FSF, Boston, MA 02111-1307 USA) comes with an Objective-C compiler since version 2.7.1. The current distributions of gcc (version 2.8.1) includes an Objective-C compiler and a run-time library. Hopefully, this makes it possible to port Eidolon on most of the current platforms and operating systems. Eventually, a problem was encountered during installation of gcc on the Parsytec CC system because Objective-C with gcc was broken for the PowerPC-AIX 4.1 architecture. Though we were able to get gcc and the Objective-C runtime library compiled properly, there was a problem apparent during run-time which was due to the fact that the Objective-C constructors were not gathered and executed upon application start-up as they were supposed to. This functionality is provided through the collect2.c program within gcc. Thanks to Scott Christley from the GNUstep project* ³ , we were able to hand patch gcc in order to make it work. A general-purpose class library of non-graphical Objective-C objects designed in the Smalltalk tradition was also installed and method calls changed accordingly. The libobjects library is to GNU's Objective-C what libg++ is to GNU's C++.

      The parallelised code consists of the following processes :

1. a master 'control' process to generate tasks;
2. a number of slave 'simulation' processes to execute tasks and to generate results;
3. an 'analysis' process to collect data and to analyse results.

      A simplified flow chart of the parallel simulator is given in Fig. 6.4. Information about the number and identity of nodes on which to perform the simulation, interaction within scatter medium and/or detector blocks (scatter-free, or considering scattering and attenuation) as well as details of scanner geometry, source and scatterer parameters are specified within series of ASCII files. The control process reads those files from the disk, performs the simulation, then writes the results. To avoid possible file corruption when two processes are editing the same file, a small delay between execution of the program on the selected processors has been arranged. Each processor writes a file named 'phantom' followed by its identification number (1 to 12). The random nature of the Monte Carlo method results in the fact that each processor could be writing data to its files at any time, making it impossible to safely add the resulting data files from many processors for analysis. A separate program is called to carry on the analysis phase and is launched only when all processors finished the calculations. The EPX partition management software (NRM) was used to define sets of partitions consisting of one processor each. Each individual processor run its own simulation with different initial configurations. One designated processor (entry node) serves as a controller which supplies other processors with some parameters that govern the simulation processes. Besides running its own sequence, the controller assesses the status of the whole simulation.

      

Fig. 6.4. : Flow chart of the parallel simulation processes

6.1.5 Distribution of random number sequences

      For many computational science applications, such as Monte Carlo simulation, it is crucial that the generators have good randomness properties. This is particularly true for large-scale simulations done on high-performance parallel computers.

      Many approaches to vectorized and parallel random number generation have been described in the literature (Sarno et al 1990, Srinavasan et al 1997). We can distinguish three general approaches to the generation of random numbers on parallel computers: centralised, replicated, and distributed. In the centralised approach, a sequential generator is encapsulated in a task from which other tasks request random numbers. This avoids the problem of generating multiple independent random sequences, but is unlikely to provide good performance. Furthermore, it makes reproducibility hard to achieve: the response to a request depends on when it arrives at the generator, and hence the result computed by a program can vary from one run to the next. In the replicated approach, multiple instances of the same generator are created (for example, one per task). Each generator uses either the same seed or a unique seed derived, for example, from a task identifier. Clearly, sequences generated in this fashion are not guaranteed to be independent and, indeed, can suffer from serious correlation problems. However, the approach has the advantages of efficiency and ease of implementation and should be used when appropriate. In the distributed approach, responsibility for generating a single sequence is partitioned among many generators, which can then be parcelled out to different tasks. The generators are all derived from a single generator; hence, the analysis of the statistical properties of the distributed generator is simplified.

      A Monte Carlo particle history is a Markov chain because the next interaction or movement of a particle is always determined by the current state of the particle. The histories of two particles become identical only when the same random number sequence is used to sample the next state. The multiplicative lagged-Fibonacci random number generator (MLFRNG) can be parameterized through its initial values because of the tremendous number of different cycles. Different streams are produced by assigning each stream a different cycle. An elegant seeding algorithm that accomplishes this is required. The MLFRNG is initialised through the function start_random_number (int seed_a, int seed_b). Supplying two seeds to start_random_number is therefore required once at program start-up before requesting any random numbers. The correspondence between pairs of seeds and generated sequences of pseudo-random numbers is many-to-one. If we choose the seeds carefully, then we can ensure that each random sequence starts out in a different cycle, and so two sequences will not overlap. Thus the seeds are parameterized (that is, sequence i gets a seed from cycle i, the sequence number being the parameter that determines its cycle). The question now becomes, how to initialise separate cycles to ensure that the seed tables on each processor are random and uncorrelated? This problem is addressed by Mascagni et al (1995) where they describe a canonical form for initialising Fibonacci generators. This canonical form is determined by l and k, but is independent of m. In general, the canonical form for initialising Fibonacci generators requires word (l -1) to be set to all zero bits and the least significant bits of all words in the register to be set to zero, with the exception of one or two characteristic bits that depend on l and k. Using this methodology, one can actually distribute random number seeds to different processors for doing the same calculation and ensure that particles for a specific process will not start at the same position of the sequence.

6.1.6 Performance and timing results

      Besides the fidelity of the simulated PET data and the ease with which software models can be prepared, the remaining important consideration in performing simulation of PET imaging systems is the computational time required to generate adequate data sets. The time needed to perform a simulation study depends on the complexity of the chosen sets of source, scatter and detector objects, and on selected interactions. Timing of a single MPC 604 processor was carried out against the HP9000/712 station. The same simulation benchmark was used in the timing where a line source was simulated in the centre of a water-filled cylindrical phantom for the ECAT 953B PET scanner operated in 3D mode (16 rings of 384 detectors each with a ring radius of 38 cm) in scatter-free imaging, and when considering attenuation and scattering within the phantom and detector blocks. The resulting computing time in minutes and the ratios between 4, 8, 12 and one single node to track 10 million annihilation pair photon histories are given in Table 6.1. The time required by slave processors to write data sets on disk is not included in the time quoted in table 6.1. The time required to perform the I/O operations and summation of data sets is negligible compared to that required for large simulation studies. A linear scaling of the computing time with the number of processors has been achieved.

      

Table 6.1. : Comparison of the computing times for different simulation studies
	HP 9000 712/60	Parsytec CC1	Parsytec CC4	Parsytec CC8	Parsytec CC12
Scatter-free imaging	72.33	25.88	6.38 (4)	3.7 (7)	2.2 (11.2)
Detector-scatter imaging	1425.65	436.55	110.26 (4)	54.6 (8.2)	36.7 (11.9)
Object-scatter imaging	2256.2	493.35	122.1 (4)	63.97 (7.7)	42.5 (11.6)
Full-scatter imaging	2540.45	657	170.5 (3.9)	83.66 (7.8)	55.7 (11.8)

      Computing times required to track 10 million annihilation pairs on both the HP 9000 712/60 workstation and the Parsytec CC system having 1, 4, 8 and 12 computing nodes are given in minutes. The speed-up achieved is also reported (in bold).

      There is no theoretical limit on the number of processors to be used. Upgrade of the system is reasonable to obtain better performance. The techniques developed in this work are also suitable for the implementation of the Monte Carlo code on other parallel computer systems. The package was also successfully ported on other platforms including Sun SPARC workstations and the Fujitsu AP3000 parallel system consisting of 80 nodes with 60 UltraSPARC 300 MHz single processors, 16 UltraSPARC 200 MHz single processors and 4 UltraSPARC 200 MHz dual processors (installed at Imperial College parallel computing centre, London).

6.2.1 Photon cross section libraries and parametrizations

      Accurate Monte Carlo simulations rely on detailed understanding and modelling of radiation transport and on the availability of reliable, physically consistent databases (Zaidi 1999a). As discussed and historically reviewed in some detail in Hubbell (1999), there exist many compilations of photon cross section data. The discrepancies and envelope of uncertainty of available interaction data have been examined from time to time, including the effects of molecular and ionic chemical binding, particularly in the vicinity of absorption edges. Calculations of photon interaction data are generally in terms of atomic cross sections, in units of cm²/atom, customarily in units of barns/atom (or b/atom) where 1 barn = 10^-24 cm².

      The Storm and Israel (1970) photon cross section data have been used extensively in medical physics. This is a 1970 compilation of data for elements 1 through 100 and energies 1 keV to 100 MeV, and contains mass attenuation coefficients, mass energy-transfer coefficients, and mass-energy absorption coefficients, presented in units of barns/atom. The medical physics community makes extensive use of these coefficients in different applications including Monte Carlo modelling. Table 6.2 lists Monte Carlo codes used to simulate nuclear medicine imaging systems together with the corresponding photon cross section libraries. In the MCNP Monte Carlo code (Briesmeister 1997), the photon interaction tables for Z=84, 85, 87, 88, 89, 91, and 93 are based on the compilation of Storm and Israel from 1 keV to 15 MeV. For all other elements from Z=1 through Z=94 the photon interaction tables are based on evaluated data from Evaluated Nuclear Data Files (ENDF) (Hubbell et al 1975) from 1 keV to 100 MeV. Data above 15 MeV for the Storm and Israel data and above 100 MeV for the ENDF data come from adaptation of the Livermore Evaluated Photon Data Library (EPDL89) (Cullen et al 1989) and go up to 100 GeV. The original EGS4 system (Nelson et al 1985) also uses compilation by Storm and Israel for the photoelectric and pair production cross sections. Recently, cross section data for this code, based on the PHOTX library (Trubey et al 1989) was created by Sakamoto (1993). ITS (Halbleib et al 1992) includes cross sections for bound electrons, the effect of which is ignored in the default EGS4 package. For the photon energy range over 1 keV to 50 MeV, of most interest to medical physicists, a particular attention is called to a recently assembled electronic database, including values of energy absorption coefficients, developed by Boone and Chavez (1996). In addition to coefficient data, other useful data such as the density, atomic weight, K, L₁, L₂, L₃, M, and N edges, and numerous characteristic emission energies are output from the program, depending on a single input variable.

      

Table 6.2. : Common Monte Carlo software packages used for simulating nuclear medical imaging systems and corresponding photon cross section libraries and parametrizations
Monte Carlo code	Photon cross section library
EGS4 (Nelson et al 1985)	Storm & Israel (Storm and Israel 1970) PHOTX (Sakamoto 1993)
ITS including TIGER CYLTRAN and ACCEPT (Halbleib et al 1992)	XCOM (Berger and Hubbell 1987)
MCNP (Briesmeister 1997)	Storm & Israel (Storm and Israel 1970) ENDF (Hubbell et al 1975) EPDL89 (Cullen et al 1989)
GEANT (Brun et al 1994)	Parametric model (Brun et al 1994)
SIMSET (Harrison et al 1993)	Parametric model (Picard et al 1993) EPDL94 (Perkins 1994), since 1997 (Kaplan 1998)
SIMIND (Ljungberg and Strand 1989)	XCOM (Berger and Hubbell 1987)
SIMSPECT (Yanch et al 1993)	Storm & Israel (Storm and Israel 1970) ENDF (Hubbell et al 1975) EPDL89 (Cullen et al 1989)
MCMATV (Smith et al 1993)	NBS (Hubbell 1969) and polynomial fitting
PETSIM (Thompson et al 1992)	Parametric model (Picard et al 1993)
Triumph's simulator (Tsang et al 1995)	GEANT's parametric model (Brun et al 1994)
EIDOLON (Zaidi et al 1998, 1999b)	GEANT's parametric model (Brun et al 1994) EPDL97 (Cullen et al 1997)

      The Lawrence Livermore National Laboratory (LLNL) houses the world's most extensive nuclear and atomic cross section database, which parametrizes the interactions of photons, electrons/positrons, neutrons, protons, and other heavy charged particles. A key feature of the LLNL database is that it is the only exhaustive interaction cross section compilation available. The Evaluated Photon Data Library, EPDL97 (Cullen et al 1997), compilation is the result of a long collaboration between LLNL and the National Institute of Standards and Technology (NIST), USA.

      The most recent EPDL97 attenuation coefficients were compared for some human tissues (water, cortical bone, brain tissue, inflated lung and air) and detector materials (NaI(Tl), BGO, BaF₂, LSO:Ce and LuAP:Ce) of interest in PET imaging to earlier libraries (PHOTX, XCOM) and parametrizations implemented in Monte Carlo simulation packages (GEANT, PETSIM) over the energy interval from 1 to 1000 keV. Although XCOM, PHOTX and EPDL97 are treated in this work as independent databases, it is recognised that they are more or less closely related. In particular, XCOM and PHOTX were both produced at NIST and EPDL97 is a result of a long and fruitful collaboration between LLNL and NIST. However, significant differences between the different libraries were reported for low energies (Cullen et al 1997) and the cross section data are sensitive to the type of interpolation used for intermediate energies. Appropriate conversions between barns and cm^-1 were made for each element. The macroscopic cross section, s (cm^-1) for a compound or a mixture is given by:

      

(6.1)

      where N_a is Avogadro's number (6.0221367 10²³ atoms/mol), r the density of the material, Z_i, A_i, w_i, and s_i are the atomic number, the atomic mass, the fraction by weight, and the atomic cross section of the i^th component of the medium, respectively. Human tissue composition data are taken from ICRU 44 report (ICRU 1989).

      For libraries where the mass attenuation coefficients are directly available for individual elements, the linear attenuation coefficients for mixtures and compounds are obtained according to:

      

(6.2)

---

---

---

---

---

6.2.2 Comparative evaluation

      A comparison between an up-to-date source of cross section data developed by LLNL in collaboration with NIST, USA, the EPDL97 (Cullen et al 1997), with other more familiar photon interaction databases: XCOM (Berger and Hubbell 1987) and PHOTX (Trubey et al 1989), and parametrizations implemented in public domain simulation packages, GEANT (Brun et al 1994) and PETSIM (Picard et al 1993) in interval from 1 to 1000 keV was performed for some human tissues and detector materials of interest in PET imaging.

      To illustrate differences between the cross section libraries, linear attenuation coefficients for water are plotted in Fig. 6.5 as a function of energy, along with the individual contributions from photoelectric absorption, incoherent scatter and coherent scatter. EPDL97 data points are given with 1 keV steps. The corresponding coefficients in different libraries appear to be qualitatively similar in this display format.

      In order to investigate differences between the EPDL97 database and the other libraries and parametrizations, percent differences between the individual coefficients are shown in Fig. 6.6 and Fig. 6.7 for water and NaI(Tl), respectively. The percent difference was calculated as , where m_c is the coefficient being compared and m_EPDL97 are the EPDL97 coefficients. The comparisons shown in Fig. 6.6 and Fig. 6.7 are continuous and were calculated at 1 keV intervals between 1 and 1000 keV for both XCOM and PHOTX, between 10 and 1000 keV for GEANT, and between 15 and 511 keV for PETSIM. It can be seen that the differences can be as large as 20%. For PHOTX, the strong variations and oscillations visible in the plots as high frequency quasiperiodic noise is thought to be due to the discrete nature of the source data and the crude interpolation used between the known coefficients which was chosen mainly to mimic the standard method used by most Monte Carlo developers.

      

Fig. 6.5. : Total and partial linear attenuation coefficients versus photon energy for water (H₂O) computed using different photon cross section libraries and parametrizations plotted on a log scale

      The different processes : (a) photoelectric absorption, (b) incoherent scattering and (c) coherent scattering contributing to (d) the total linear attenuation coefficient are shown. The coefficients appear to be qualitatively similar in this display format.

      The parametric model used in PETSIM has an apparent problem in modelling the incoherent coefficients for low energies (<100 keV). This can be explained by the fact that XCOM uses a combination of the Klein-Nishina formula and non-relativistic Hartree-Fock incoherent scattering functions (which is about equal to 1 above 100 keV) to calculate the incoherent scattering, whereas only the Klein-Nishina relationship is used in PETSIM (Picard et al 1993). The discrepancy is even greater for high atomic number materials since the electrons can no longer be considered free. However, this should not affect the accuracy of photon transport simulations within the phantoms since the low energy cut-off is generally higher than 200 keV. But, it might have implications in multi-spectral imaging where acquisitions in lower energy windows are simulated and used for example in scatter correction based on a multi-energy window setting.

      

Fig. 6.6. : Comparisons (percent differences) between the different libraries and the EPDL97 database for water (H₂O)

      The coefficients shown are: (a) total, (b) photoelectric, (c) incoherent and (d) coherent. The comparisons were calculated at energies given to the keV resolution between 1 and 1000 keV for both XCOM and PHOTX, between 10 and 1000 keV for GEANT, and between 15 and 511 keV for PETSIM.

      

Fig. 6.7. : Comparisons (percent differences) between the different libraries and the EPDL97 database for sodium iodide (NaI(Tl))

      The coefficients shown are: (a) total, (b) photoelectric, (c) incoherent and (d) coherent. The comparisons were calculated at energies given to the keV resolution between 1 and 1000 keV for both XCOM and PHOTX, between 10 and 1000 keV for GEANT, and between 15 and 511 keV for PETSIM.

      The average positive difference between the GEANT and the EPDL97 coefficients for water are 0.3%, 9.61%, 0.31%, and 0.85%, for the total, photoelectric, incoherent and coherent attenuation coefficients, respectively. Larger differences can be seen near the K or L edges for some of the compounds. A good agreement is observed between EPDL97 and XCOM for all compounds except the low energy region for the coherent attenuation coefficients. With some exceptions, most of the data points fall within the ±3% region.

      

Table 6.4. : Global estimates of relative percent differences between PHOTX and XCOM libraries and EPDL97 coefficients. Mean relative errors with standard deviations over the energy range 250 to 511 keV are shown
	XCOM				PHOTX
	Total	Photoelectric	Incoherent	Coherent	Total	Photoelectric	Incoherent	Coherent
Water	0.24±0.17	1.74±0.78	0.24±0.18	0.40±0.29	0.28±0.22	7.67±7.29	0.26±0.20	3.74±3.20
Bone (Cort)	0.21±0.16	1.40±0.58	0.20±0.15	0.29±0.17	0.59±0.29	7.53±7.28	0.50±0.19	4.0±3.08
Lung (Infl)	0.23±0.18	1.72±0.74	0.23±0.18	0.42±0.29	0.63±0.22	7.09±7.37	0.61±0.20	4.07±3.09
Brain tissue	0.23±0.17	4.28±0.69	0.24±0.18	0.75±0.36	0.64±0.22	7.73±7.58	0.61±0.20	3.87±3.22
Air	0.29±0.21	1.58±0.68	0.29±0.21	0.33±0.24	0.60±0.38	7.83±7.66	0.73±0.31	4.15±3.22
NaI(Tl)	0.44±0.39	1.44±0.70	0.12±0.07	1.58±0.79	2.83±3.07	6.37±6.11	0.32±0.12	4.85±2.06
BGO	0.73±0.45	1.30±0.56	0.11±0.08	2.25±1.32	3.52±3.46	5.01±4.58	0.22±0.11	3.83±1.50
BaF2	0.42±0.35	1.33±0.61	0.12±0.06	1.56±0.89	2.92±3.15	6.32±6.03	0.32±0.12	4.78±1.94
LSO:Ce	0.50±0.32	0.91±0.44	0.13±0.08	1.98±1.16	5.14±3.92	7.68±5.85	1.54±0.13	5.50±1.90
LuAP:Ce	0.43±0.29	0.83±0.42	0.09±0.06	1.98±1.16	4.96±3.85	7.69±5.89	1.44±0.12	5.45±1.89

      

Table 6.5. : Global estimates of relative percent differences between GEANT and PETSIM parametrizations and EPDL97 coefficients
	GEANT				PETSIM
	Total	Photoelectric	Incoherent	Coherent	Total	Photoelectric	Incoherent
Water	0.30 ± 0.13	9.61 ± 6.82	0.31 ± 0.13	0.85 ± 0.44	0.11 ± 0.08	5.12 ± 3.79	0.29 ± 0.11
Bone (Cort)	0.32 ± 0.16	3.10 ± 1.57	0.33 ± 0.16	1.65 ± 0.80	0.42 ± 0.17	4.21 ± 3.15	0.50 ± 0.14
Lung (Infl)	0.39 ± 0.20	8.37 ± 5.93	0.39 ± 0.20	0.86 ± 0.48	0.25 ± 0.08	4.71 ± 2.75	0.14 ± 0.11
Brain tissue	0.39 ± 0.20	8.18 ± 5.81	0.40 ± 0.20	0.86 ± 0.49	0.73 ± 0.08	4.91 ± 3.34	0.34 ± 0.11
Air	0.33 ± 0.24	9.73 ± 6.90	0.33 ± 0.24	1.24 ± 0.97	0.24 ± 0.16	4.87 ± 3.43	0.34 ± 0.22
NaI(Tl)	0.98 ± 0.68	1.73 ± 1.58	0.51 ± 0.33	2.36 ± 1.33	5.41 ± 0.46	2.29 ± 1.99	1.96 ± 0.57
BGO	1.01 ± 0.72	1.96 ± 1.02	1.60 ± 0.46	5.03 ± 2.68	1.92 ± 0.78	1.76 ± 1.31	13.94 ± 0.72
BaF2	0.95 ± 0.71	1.75 ± 1.51	0.51 ± 0.32	3.32 ± 1.77	5.03 ± 0.35	2.28 ± 1.73	2.20 ± 0.55
LSO:Ce	1.33 ± 0.78	1.38 ± 0.96	1.81 ± 0.47	5.09 ± 2.45	3.40 ± 0.36	1.83 ± 1.42	7.46 ± 0.61
LuAP:Ce	1.30 ± 0.76	1.38 ± 0.96	1.74 ± 0.45	5.06 ± 2.43	9.24 ± 1.02	1.83 ± 1.37	5.66 ± 0.50

      Mean relative errors with standard deviations over the energy range 250 to 511 keV are shown.

      The mean relative errors for water linear attenuation coefficients in the energy range over 250 to 511 keV between XCOM and EPDL97 are (1.74±0.78)%, (0.24±0.18)%, (0.40±0.29)% and (0.24±0.17)% for photoelectric effect, incoherent scattering, coherent scattering and total linear attenuation coefficients, respectively. They are (5.12±3.79)%, (0.29±0.11)% and (0.11±0.08)% between PETSIM and EPDL97 for photoelectric effect, incoherent scattering and total linear attenuation coefficients, respectively. While the overall average difference was small, there was a larger difference (>50%) between cross sections for some compounds. The overall percent differences are summarised in tables 6.4 and 6.5 for the different materials investigated. 6.2.3 Effect on simulation of PET tomographs

      The scatter fraction is of great importance for quantitative estimation of the scattering contribution. The ratios of scatter fractions estimated using the different cross section libraries to those estimated using the EPDL97 library as a function of the lower level threshold (LLT) are shown in Fig. 6.8. There is little difference between estimated scatter fractions. The effect of the cross section library can hardly be observed when calculating this quantity which involves only the resulting numbers of unscattered and scattered events. However, it can be seen that PETSIM slightly underestimates the scatter fraction as compared to other libraries.

      A detailed investigation of the error in scanner input parameters on the accuracy of the simulated data particularly the energy and spatial resolution of the tomograph was not helpful in drawing a general conclusion regarding the sensitivity to a particular parameter. That is, in the absence of an accurate model of the variations of the crystal energy response within a single block (Vozza et al 1997), changing the energy resolution parameter alone revealed small variations of the scatter fraction similarly to the results reported here when changing the photon cross section library. Determination of the sources of systematic errors is not an easy task.

      

Fig. 6.8. : Ratios between scatter fractions estimated using the different photon cross section libraries and those estimated using the EPDL97 library for the ECAT 953B as a function of the lower level threshold

      

Fig. 6.9. : Comparison between simulated scatter profiles of a 20 cm diameter uniformly water-filled right circular cylinder generated using the different photon cross section libraries: 380 keV LLT (left) and 450 keV LLT (right)

      Scatter profiles are summed over all angles.

      

Fig. 6.10. : Comparison between simulated profiles of an axial line source centred in a 20 cm diameter water-filled cylinder generated using the different photon cross section libraries: (a) 300 keV LLT and (b) 400 keV LLT

      Scatter profiles for a uniformly distributed source in a cylindrical phantom are illustrated in Fig. 6.9 for a LLT of 380 keV and 450 keV. Raising or lowering the LLT setting changes the shape of the scatter profile; however, even with a higher LLT, the FWHM of the scatter response function is still significant. The LLT setting affects the accuracy of a correction method based on estimating the scatter profile from the scatter 'tails'. Energy settings of 380 keV place the energy threshold at the lower energy tail of the photopeak for detector blocks with 23% energy resolution. A 511 keV photon can scatter 30° and still possess 450 keV of energy. PETSIM and GEANT slightly overestimate scatter components for both LLT settings as compared to EPDL97. It has been shown that using different libraries can lead to differences in the scatter distributions of approximately ±10%. In Figs. 6.10a and 6.10b, the difference between summed profiles of the line source simulated using the EPDL97 and PETSIM photon cross section libraries are shown for two different settings of the LLT. The profile demonstrates how well defined the scatter portion of the curve is outside the object. In this logarithmic display format, the differences are just barely visible on the figure. However, a detailed analysis revealed that using different libraries can lead to differences in the scatter distributions of approximately ±6%. This is in agreement with the results obtained from Fig. 6.9.

      

Fig. 6.11. : (a) Integral profile through a 2D projection (single angular view) of a simulated Hoffman brain phantom before and after scatter correction, (b) Integral profile through a 2D projection generated using the different photon cross section libraries, and (c) same as (b) after applying scatter correction

      Fig. 6.11(a) shows profiles through a simulated 2D projection of the 3D brain phantom without and with scatter correction using the convolution-subtraction method (Bailey and Meikle 1994). Comparisons between the profiles estimated using different libraries before and after applying scatter correction are illustrated in Figs. 6.11(b) and 6.11(c), respectively. An example of a representative slice from the reference image of the Hoffman brain phantom is shown in Fig. 6.12 together with reconstructed images using the reprojection algorithm (Kinahan and Rogers 1989) in a comparison of the 3D images obtained without (Fig. 6.13) and with scatter compensation (Fig. 6.14). Profiles through the images are also shown. Results using this approach for scatter correction in cerebral 3D PET scans (homogeneous attenuating region) have proven to be accurate (Bailey and Meikle 1994). Although the data presented here were obtained from low count studies, there is a clear difference in the resulting profiles on reconstructed images when using different libraries for generating data sets. The data sets were obtained with the same statistics and random number seeds, so the differences are due mainly to the cross section libraries used to generate them. This might strongly influence the parameters derived for modelling the scatter response function and the conclusions drawn when evaluating scatter correction techniques using Monte Carlo simulated source distributions.

      

      

      

      

      

Fig. 6.12. : Selected slice from the 3D reference image of the Hoffman brain phantom and horizontal profile drawn at the centre of the slice

      

      

Fig. 6.13. : Reconstructed images of the Hoffman brain phantom and horizontal profiles for non-scatter corrected data sets generated with different photon cross section libraries: (A) EPDL97, (B) GEANT, (C) XCOM, (D) PHOTX, and (E) PETSIM

      

Fig. 6.14. : Reconstructed images of the Hoffman brain phantom and horizontal profiles for scatter corrected data sets using the convolution-subtraction method generated with different photon cross section libraries: (A) EPDL97, (B) GEANT, (C) XCOM, (D) PHOTX, and (E) PETSIM

      Accurate Monte Carlo modelling is still the gold standard for the simulation of PET data and is an indispensable tool to develop and evaluate scatter models. The recent interest in the development of high-resolution PET scanners and multiple-processor parallel processing systems has created new opportunities for Monte Carlo simulation in PET imaging (Zaidi 1999a).This motivates the question: what is the ultimate accuracy that can be achieved using the Monte Carlo method? The accuracy of simulations is limited by several factors such as accuracy in modelling of annihilation photon non-collinearity, positron range, the photon transport model and associated physical databases used, and the obvious complexity in modelling the human body. Because the overall accuracy is a combination of all these components, when designing a Monte Carlo package that simulates a PET imaging system it is quite important to understand the effect of each component.

      Data on the scattering and absorption of photons are fundamental for all Monte Carlo calculations since the accuracy of the simulation depends on the accuracy in the probability functions, i.e. the cross section data. An aspect that deserves attention is the fact that there exists a difference in the physical cross section data used in available Monte Carlo codes. This is of special importance when comparing results from different codes. The use of different cross section data and approximations will usually yield different results and the accuracy of the results is generally hard to quantify. However, given the results presented in this paper, it is expected that photon transport in non-uniformly attenuating medium will exhibit noticeable differences when using different cross section libraries. Therefore, the reported errors should have a significant effect on most, if not all realistic simulations of nuclear medicine applications. In the first version of Eidolon, interaction cross sections and scattering distributions were computed from parametrizations implemented in the GEANT simulation package (Brun et al 1994). It has been demonstrated in this work that different photon cross section libraries show quite large variations as compared to the most recent EPDL97 nuclear data files. Both XCOM and PHOTX use form factor corrections for incoherent cross sections, but do not use anomalous scattering factors, so their low energy coherent cross sections is orders of magnitude too large.

      In the last few years, the LLNL photon and electron data bases have been greatly improved in terms of the detail included as well as the accuracy of the data. At the same time there has been an enormous increase in available inexpensive computer power. The combination of improved photon and electron data bases and increased available inexpensive computer power allows us today to calculate results in greater detail, with greater accuracy, using accurate methods such as Monte Carlo modelling. The cross section values produced by the LLNL (in collaboration with NIST) are thought to be the most up-to-date and accurate coefficients available, and the LLNL database is now a national and international standard (Boone and Chavez 1996). EPDL97 is already a standard in the nuclear reactor industry, since its adoption by the Cross Section Evaluation Working Group (CSEWG) for use with ENDF/B-VI. In our opinion, it would be desirable that the EPDL97 is adopted as an International standard not because it is perfect, but rather because it could serve as a standard against which everyone can test, and based on user feedback the authors of this library would then maintain EPDL97 up-to-date to reflect the latest measurements and calculations. Therefore, we highly recommend that EPDL97 cross section data serve as a standard and supersede earlier cross section libraries used in Monte Carlo codes to simulate medical imaging systems. This will help to eliminate potential differences between the results obtained with different codes. The results presented in this section have shown that the effect of using different photon cross section libraries is noticeable on both generated projection data and reconstructed images (Figs 6.9 - 6.14). The data sets were obtained with the same statistics and same random number seeds. So, the differences are due mainly to the cross section libraries used to generate them. The library was carefully designed in the form of look-up tables for efficiency in data storage, access, and use by the Monte Carlo software (Zaidi 1999b). This work presented just a direct comparison of the photon interaction data. It doesn't address the companion Evaluated Electron Data Library (EEDL), and Evaluated Atomic Data Library (EADL), which are needed to complete the picture of photon transport and energy deposition. Those libraries are the only complete and consistent set of data libraries available to treat the complete coupled photon-electron transport problem. These libraries include details that were not previously available, or not considered when performing calculations using traditional photon interaction data.

      Together with the optimisation of the computing time performances of the Monte Carlo software, photon transport in 3D PET could be efficiently modelled to better understand scatter correction techniques. The accuracy in correction methods can thus be evaluated in an unbiased way since systematic errors can be controlled. The improved accuracy and speed-up due to parallelisation makes Monte Carlo simulation more attractive as a tool for modelling photon transport in 3D PET. These issues are further investigated in the following chapter of this thesis.

7.3.1 Scatter distribution functions

      It has been shown (Fig. 7.5) that the projections of a line source placed in a uniform water-filled cylinder are dominated in the wings by the object scatter and in the peak by the unscattered photons. As expected, the object scatter is fairly well described by the monoexponential functions. The detector scatter, on the other hand, is closely located around the primary peak. This observation seems reasonable, since the point of interaction of a photon scattered in the crystal is closer to the centroid of interaction than that of a photon scattered in the object, leading to a smaller deviation of the scattered photons. Since the detector scatter has no significant contribution to the wings of the scatter distribution function, it can be added to the unscattered component.

      This discussion remains valid, when the line source is shifted out of the symmetry centre. The amplitude of the short side of the projection, compared to that of the symmetrical case is increased, since the path length of the photons through the phantom becomes shorter, whereas the amplitude of the long side is decreased, due to a longer pathway through the attenuating medium. The wings of the scatter distribution functions are still well parametrised by mono-exponentials, but with different parameters for the left and the right side. The data presented reveal that the amount of multiple scatter is at least 10-15% of the scattered photons. The rationale of several scatter correction algorithms is to consider the multiple scatter by rescaling the single scatter component. This approach seems to be questionable in view of the results obtained by the Monte Carlo calculations. Obviously, a rescaling does not change the shape of the different scatter components.

      Fig. 7.6 shows scatter distribution functions when the line source is located at the centre of the FOV (left) and displaced 5 cm radially off-centre (right). The projections for the line source at the centre are symmetrical and were summed over all projection angles. This is not the case when the source is moved off-centre. In this case, a profile of a single angular view was used. The scatter function kernel was assumed to be mono-exponential with exponent decreasing linearly as a function of the radial position. A linear regression analysis on the wings of the log plot of the projection data was used to estimate the slope, a (in units of cm^-1). The data in the extreme bins of the projection were excluded to avoid edge effects.

      

Fig. 7.6. : Illustration of the scatter response functions when the line source is located at the centre of a water-filled cylindrical phantom (left) and displaced 5 cm radially (right)

      The projections for the line source placed at the centre of the FOV are symmetrical and were summed over all projection angles.

      The scatter function was examined with respect to three parameters: (i) variation in lower level threshold setting; (ii) variation in radial position of the line source within the cylindrical phantom; and (iii) variation in the diameter of the cylindrical phantom. For the case when the line source was located at the centre of the cylinder, table 7.2 shows a progressive increase in the slope of the scatter function from 0.037 cm^-1 to 0.099 cm^-1 with increasing lower level threshold. This is consistent with the fact that the scatter tails are less pronounced when the lower level threshold is increased. For comparison purposes, the measured values of the srf slopes for the ECAT 953B are also shown. Note that when the line source was located at the centre of the cylinder and the LLT was varied, the source is in the same location in the sinograms at all projections. This allows to estimate a and its standard deviation from all 2D radial measurements. When the line source was translated radially along the x-axis of the scanner, only the projection corresponding to the first view (F = 0) was considered. The effect of increasing the axial acceptance angle of the sinograms up to dr = 12 was negligible. A good agreement between measured and simulated data was observed except for a LLT setting of 250 keV.

      

Table 7.2. : Measured and simulated values of a (in cm-1), the slope of the scatter response function exp(-ar) for different lower level energy window settings obtained by linear regression analysis
Lower energy window setting (keV)
250
300
350
380
400

      Table 7.3 shows the variation of a as a function of the radial displacement. The slope of the scatter function is seen to increase on one side of the peak while decreasing on the other side. A variation by a factor of over two is also observed. As the line source is moved radially, the slope decreases (i.e. an increase in scatter) outside the object up to some limit, while the slope increases inside the object (i.e. a decrease in scatter), which is consistent with the observations made by Bailey and Meikle (1994). By observing the profiles of the line source at different radial locations, it was noticed that many scattered events are still recorded even when the line source is placed just inside the edge of the cylinder (9.5 cm).

      

Table 7.3. : Simulated and measured variation of a (in cm-1), the slope of the scatter function exp(-ar) for either side of the log profile at different radial locations
Radial distance from centre (cm)
0
5
7.5
9.5

      All simulations were made with a lower energy window threshold of 380 keV.

      On the basis of the results presented in the preceding tables, the effect of the size of the cylinder on the spatial variation of scatter estimation was further investigated. The variation of a as a function of the radius of the cylindrical phantom is illustrated in table 7.4. It can be seen that the slope increases with increasing the size of the cylinder on both sides of the peak. This parameter changes by a factor of around two when increasing the radius of the cylinder from 5 cm to 20 cm.

      

Table 7.4. : Variation of a (in cm-1), the slope of the scatter function exp(-ar) for different sizes of the water-filled cylinder serving as scattering medium
Radius of the cylinder (cm)	slope of the srf (cm-1)
5	0.069
8	0.074
10	0.084
12	0.090
15	0.101
18	0.119
20	0.129

      All simulations were made with a lower energy window threshold of 380 keV.

7.3.2 Scatter fractions

      The variation of the scatter fraction (SF) was investigated in the same way as the scatter distribution function. However, in addition to the line source geometry, the variation of the scatter fraction in a uniform cylindrical phantom as a function of its size and for three lower energy thresholds (250, 380 and 450 keV) was also studied. The scatter fraction is estimated directly from the results of the Monte Carlo simulation as the ratio between scattered and total coincidences. Fig. 7.7 shows the SF versus the radius R of the cylindrical phantom. Obviously, SF(R) increases with increasing R. To fit the data, a simple function was derived . As reported by Adam et al (1996), this function fits the scatter fraction of a uniform activity distribution in cylindrical phantoms (Fig. 7.7), as well as centred (Fig. 7.8) and a 5 cm off-centred line sources (Fig. 7.9) reasonably well, although some crude simplifying assumptions have been made for the derivation of this function. This simple formula for the scatter fraction was derived by assuming that the probability of a photon not to be scattered is proportional to (where x is the pathlength of the photon). The same assumption is valid for the second photon of the pair, leading to a probability for both photons not to be scattered proportional to . In the case of a centred line source, the total length through the cylinder is always twice the radius R. Thus, the probability of two coincident photons not to be scattered is proportional to . The probability of two coincident photons to be scattered could then be obtained since the two probabilities are complementary. Moreover, the percentage of Compton interactions in the object consist of more than 99.7% at 511 keV for H₂O which renders the number of interactions by photoelectric absorption negligible. The second argument is that for the majority of detected scattered events, only one of the two annihilation photons undergoes a single Compton interaction. The single scatter distribution was reported to consist of more than 75% of detected scattered events (Ollinger 1996). The fit function is therefore and the normalisation constraint that the limit for should be 1 is verified. This simple functional dependence for the scatter fraction was also assumed for homogeneously distributed activities.

      

Fig. 7.7. : Monte Carlo calculations of the variation of the scatter fraction as a function of the radius R of the cylindrical phantoms for a homogeneous activity distribution for three different LLT settings

      The fitted curves are also shown.

      With increasing the lower energy discrimination level, the scatter fraction decreases. An increase of the energy threshold reduces the number of registered scattered coincidences, but it also decreases the sensitivity leading to a reduced number of true coincidences. As expected (and confirmed by Monte Carlo simulations), the largest effect is achieved for multiple scattered photons and photons that are scattered under a large angle.

      

Fig. 7.8. : Monte Carlo calculations of the variation of the scatter fraction as a function of the radius R of the cylindrical phantom for a central line source for three different LLT settings

      The fitted curves are also shown.

      

Fig. 7.9. : Monte Carlo calculations of the variation of the scatter fraction as a function of the radius R of the cylindrical phantom for a line source 5 cm off-axis for three different LLT settings

      The fitted curves are also shown.

      It was observed that, for the same phantom size, the scatter fraction for a line source is always larger than that for a homogeneously distributed activity. The higher scatter fraction for a line source centred in the cylinder compared to a homogeneous activity distribution in the same cylinder can be explained as follows. For a line source, the two annihilation photons must always travel the maximum path length through the phantom, which is not the case for a homogeneously distributed activity. The activity distribution in cerebral scans can be considered as a mixture of both components. Consequently, its scatter fraction is between the two extremes. The scatter fraction for a non-uniformly attenuating medium in the thorax region (oncology clinical study described in section 7.6), compared to a cylindrical phantom with a comparable (similar) cross-section (R=15 cm), is larger than that for a homogeneously distributed activity, but smaller than that for a line source.

      Better fitting was achieved using third order polynomials, however, preference was given to the simple parametrization presented above for both simplicity and reproducibility purposes. Since the function fits the object scatter fraction for uniform activity distributions and line sources (centred and off-centre) in cylindrical shaped objects reasonably well (Figs. 7.7-7.9), it allows an estimation of the amount of object scatter for simple activity distributions. According to the definition presented in table 7.1, unscattered and detector scatter components are not indistinguishable until they enter the crystal volume since they are both not scattered in the object. This means that the detector scatter is correlated with the primaries and not with the object scatter, which was already assumed in other studies (Bentourkia et al 1995). Hence, the detector scatter becomes relevant for small objects, such as small animals, scanned on high-resolution scanners, when the object scatter fraction is very small (Ziegler and Kuebler 1993).

7.4.1 Energy window-based approaches

      Multiple energy window methods have been originally developed for SPECT and have been in use for more than 15 years (Zaidi 1996b). The development of 3D acquisition mode and improvements in the detector energy resolution have allowed the implementation of scatter correction based on the analysis of energy spectra. Several groups investigated the potential of acquiring data in two (Grootoonk et al 1996, Bendriem et al 1993), three (Shao et al 1994) and multiple (Bentourkia et al 1995) energy windows to develop corrections for scattering in 3D PET.

      Two variants of dual-energy window techniques have been proposed: methods estimating the scatter component in the photopeak window from the events recorded in a lower energy window placed just below the photopeak and methods estimating the unscattered component in the photopeak window from the unscattered counts recorded in a high energy window in the upper portion of the photopeak. The dual energy window (DEW) technique (Grootoonk et al 1996) belongs to the former while the estimation of trues method (ETM) (Bendriem et al 1993) belongs to the latter.

      In the practical implementation of the dual-energy window method on the ECAT 953B scanner, detected coincidence events are assigned to the upper energy (UE) window when both photons deposit energy between 380 keV and 850 keV, or to the lower energy (LE) window when one or both photons deposit within 200 keV and 380 keV (Grootoonk et al 1996). Both energy windows are assumed to contain object scattered and unscattered events. The coincidence events in both energy windows in each sinogram can be expressed as follows:

      

(7.5a) (7.5b)

      where _unsc and _sc refer to the unscattered and scattered components, respectively. These two equations contain four unknown parameters and so additional data are required to solve the system of equations and extract the unscattered components. Let us define R_sc as the ratio of scattered events and R_unsc as the ratio of unscattered events detected in both windows:

      

,

(7.6)

      Thus Eq. (7.5a and 7.5b) can be re-arranged to give the following expression for scattered events in the upper window:

      

(7.7)

      UE_unsc can then be obtained from equation (7.7) by subtracting UE_sc from UE. The scaling parameters R_sc and R_unsc are derived from measurements of the ratios of counts from line sources due to unscattered (measurements in air) and scattered events (measurements in a head-sized phantom) in the two energy windows, respectively.

      The estimation of trues method (Bendriem et al 1993) consists in acquiring data simultaneously into two energy windows: a high window with a lower energy threshold higher than 511 keV and a regular acquisition window including the higher window. Therefore, both windows have the same upper energy threshold (UET) value. The hypothesis of this method is that the number of unscattered coincidences recorded in a given energy range depends on the energy settings of the window and the angle of incidence of the annihilation photons on the detector face. Hence, the unscattered component in the high energy window can be related to the unscattered coincidences in the standard wider window through a function of the energy settings, the radial position in the sinogram for a given LOR and the axial opening for a given radial position. This calibrating function is in principle independent of the source distribution. The unscattered component in the wide energy window can thus be calculated and subsequently subtracted from the data recorded in the regular window to produce a scattered sinogram. The unscattered component in the regular window is then obtained by smoothing that sinogram and subtracting it from the data recorded in the standard window.

      The triple energy window method (Shao et al 1994) was suggested as an extension of the DEW technique. Coincidence events are recorded in three windows: two overlapping windows having the same UET settings (450 keV) located below the photopeak window and a regular window centred on the photopeak and adjacent to the low windows. A calibrating function that accounts for the distribution of scattered coincidences at low energies is obtained by calculating the ratio of the coincidence events recorded in both low energy windows for the scanned object and for a homogeneous uniform cylinder. The scatter component in the standard acquisition window is then estimated from the calibrating function and the narrower low energy window.

      The multispectral method is based on the acquisition of data in a very large number (typically 256) of windows of the same energy width (16x16 energy values for the two coincident photons). The spatial distribution of scattered and unscattered components in each window can be adjusted using simple mono-exponential functions (Bentourkia et al 1995). The method presents the advantage of increased sensitivity due to smaller rejection of detected coincidence events (~10%) since the position of the scattered photons in the detectors can be restored. The statistical noise in each window and the necessity to accommodate the acquisition hardware with the required electronic boards to allow acquisition in multiple windows remain the major obstacles for implementation of the method on commercial PET scanners.

7.4.2 Convolution-deconvolution based approaches

      Techniques based on convolution or deconvolution estimate the distribution of scatter from the standard photopeak data. The scatter fraction (SF) which gives an indication about the expected amount of scatter and the scatter response function (srf) which defines the spatial distribution of scatter in the photopeak data are usually the two parameters required for estimation of scatter component and need to be determined a priori. A pure additive model of the imaging system is generally assumed where the recorded data are composed of an unscattered and a scattered component plus a noise term due to statistical fluctuations, and can be written in the following form:

      

(7.8)

      where p_o are the observed data, p_u and p_s are the unscattered and scattered components respectively, and is the noise term. The problem to be addressed consists in estimating the unscattered distribution (p_u) from the measured data (p_o) contaminated by scatter, or alternatively estimate the scattered component (p_s) and then derive p_u. The proposed methods differ in the way the scatter function is defined.

      The convolution-subtraction (CVS) technique developed for 3D PET (Bailey and Meikle 1994) operates directly on projection data (pre-reconstruction correction) and is based on the assumption that the convolution of the source distribution with the srf gives an estimate of the scatter component. Two assumptions have been made: the stationary and nonstationary assumptions. In the stationary assumption, the scatter is assumed to be analytically defined and not dependent on the object, activity distribution, etc. The nonstationary assumption overcomes this problem by taking into consideration the dependence of scatter upon source locations, object size, detector angle, etc. Using the stationary assumption, the scatter component can be related to the unscattered component by the convolution relation:

      

(7.9)

      where denotes the convolution operator. The scatter estimate (p_s) is then subtracted from the measured data (p_o) after scaling by an appropriate scatter fraction (SF). The process can be repeated iteratively with each step using the previous estimate of the scatter-free distribution as input to the scatter estimation (Bailey and Meikle 1994). The observed data are used as a first approximation to the unscattered distribution (p_u) and the process is repeated iteratively with each step using the previous estimate of the scatter-free distribution as input to the scatter estimation:

      

(7.10)

      where the '' indicates that the parameter is a rough estimate of the scatter and n is the iteration number. The rationale is that with each iteration, the input to the scatter estimation step more closely approximates p_u. It was also suggested to use a damping factor to prevent oscillation in the result.

      The convolution-subtraction approach can also be applied to reconstructed images (post-reconstruction). In this case, the scatter estimates are reconstructed and then subtracted from the non-corrected reconstructed images of the acquired data (Lercher and Wienhard 1994).

      Scattered photons degrade the point-spread function (PSF) of the PET camera; the long tails of the PSF are mainly due to scatter. Thus, deconvolution methods which correct the images for the PSF will also, implicitly correct for scatter. In general, the PSF will act as a low-pass filter. Deconvolution will restore the high-frequency contents of the signal emitted by the object, at least as long as they have not yet been completely removed by PSF. Not only the high-frequencies in the signal are restored, but also the high-frequency noise is amplified, which in turn can degrade the image again. Therefore, the restoration filter is often combined with a low-pass filter that balances that image improvement by deconvolution and its degradation due to amplification of noise. Well known examples of these filters are the Wiener and Metz filters. Links et al (1992) studied the use of two-dimensional Fourier filtering to simultaneously increase quantitative recovery and reduce noise. The filter is based on the inversion of the scanner's measured transfer function, coupled with high-frequency roll-off. In phantom studies, they found improvements in both 'hot' and 'cold' sphere quantification. Fourier-based image restoration filtering is thus capable of improving both accuracy and precision in PET.

      The curve fitting approach is based on the hypothesis that detected events assigned to LORs outside of the source object must have scattered and that the scatter distribution corresponds to a low-frequency component that is relatively insensitive to the source distribution. Estimation of the unscattered component can thus be performed in three successive steps: (i) fitting the activity outside the source object with an analytical function (e.g. Gaussian), (ii) interpolating the fit inside the object, and (iii) subtracting the scatter component from the observed data (Cherry and Huang 1995). The accuracy of this class of scatter correction methods depends on how accurately the scatter components can be estimated. The appropriate choice of a set of parameters which should be optimised for each PET scanner and for different radioactivity distributions is the dominant factor.

7.4.3 Approaches based on direct calculation of scatter distribution

      This class of methods assume that the distribution of scattered events can be estimated accurately from the information contained in the emission and transmission data. For the majority of detected scattered events, only one of the two annihilation photons undergoes a single Compton interaction. The rationale of most of these methods is that the overall scatter distribution can be computed from the single scatter distribution (~75% of detected scattered events) and that this latter can be scaled to model the distribution of multiple-scattered events (Ollinger 1996). The multiple scatter distribution is generally modelled as an integral transformation of the single scatter distribution. The same approach can also be applied to scatter estimation in transmission imaging. Monte Carlo simulation studies of various phantom geometries demonstrated the accuracy of this method for fully 3D PET imaging by direct comparison of analytic calculations with Monte Carlo estimates (Barney et al 1991).

      The model-based scatter correction method developed by Ollinger (1996) uses a transmission scan, an emission scan, the physics of Compton scatter, and a mathematical model of the scanner in a forward calculation of the number of events for which one photon has undergone a single Compton interaction. A single-scatter simulation technique for scatter correction where the mean scatter contribution to the net true coincidence data is estimated by simulating radiation transport through the object was also suggested and validated using human and chest phantom studies (Watson et al 1997). However, the above methods do not correct for scatter from outside the field-of-view. This effect can be directly taken into account by acquiring short, auxiliary scans adjacent to the axial volume being investigated. This implicitly assumes that the distribution of scatter from outside the FOV has the same shape as that of scatter from inside the FOV. These extra data are naturally available in whole-body imaging. However, this method is impractical for isotopes with a short half-life or rapid uptake relative to the scanning interval.

      The true scatter component (being an experimentally impossible measurement) can be accurately estimated using rigorous Monte Carlo simulations. Given a known radioactive source distribution and density of the object, Monte Carlo techniques allow detected events to be classified into unscattered and scattered events and the scatter component to be determined in an independent way. However, the source and scattering geometry is generally not known in clinical studies. Levin et al (1995) used filtered-backprojection reconstructions as an estimate of the true source distribution. This input image is then treated as a 3D source intensity distribution for a photon-tracking simulation. The number of counts in each pixel of the image is assumed to represent the isotope concentration at that location. The image volume planes are then stacked and placed at the desired position in the simulated scanner geometry assuming a common axis. The program then follows the history of each photon and its interactions in the scattering medium and traces escaping photons in the block detectors in a simulated 3D PET acquisition. The distributions of scattered and total events are calculated and sorted into their respective sinograms. The scatter component is equal to the difference between measured data and the scaled and smoothed scatter component. To reduce the calculation time, coarser sampling of the image volume over regions equivalent to 4x4 pixels was adopted assuming that the Compton scatter distribution varies slowly over the object. For obvious reasons, the implemented method do not correct for scatter from outside the field-of-view and further refinements of the technique were required to take this effect into account as described in section 7.5.1.

7.4.4 Iterative reconstruction-based scatter correction approaches

      Improvement in computer speed and recent advances in acceleration of reconstruction algorithms (Hudson and Larkin 1994) have led to renewed interest in algorithms which incorporate scatter as part of the emission model. Iterative reconstruction-based scatter compensation has received considerable attention during the last decade. Earlier studies (Frey et al 1993, Beekman et al 1997b) indicate that when scatter is modelled in iterative reconstruction, higher contrast and lower noise can be obtained than when the technique of scatter subtraction is used. In addition, model-based iterative reconstruction does not require patient-dependent parametrization. Development of scatter models which can be incorporated in statistical reconstruction (e.g. ML-EM) of PET continues to be appealing, however, implementation must be efficient to be clinically applicable.

      One of the requirements of this method is to compute the scatter response function at each point in the attenuator with respect to each projection view. An efficient algorithm for scatter estimation was recently described in which the spatial scatter distribution is implemented as a spatially invariant convolution for points of constant depth in tissue (Hutton and Baccarne 1998). The scatter estimate is weighted by a space-dependent build-up factor based on the measured attenuation in the body. This approach was validated by Monte Carlo simulation studies of a realistic thorax phantom.

      The normal approach for implementation of a scatter model is to incorporate the scatter estimation directly in the transition matrix (a_ij), although efficiency has been improved by utilising a dual matrix approach (Kamphuis et al 1998) in which scatter is incorporated in the forward projection step only. In this case, a_ij is considerably larger than is necessary if only attenuation and geometric factors are included and computation is therefore slow since scatter is essentially recalculated and added each iteration. In the case of constant pre-calculated (or measured) scatter, for example, as in multi-energy window methods, this scatter estimate either can be subtracted directly from projections (p_i) prior to reconstruction, or alternatively can be introduced in the denominator of the ML-EM equation (Eq. 7.11). This latter approach results in better noise properties than direct subtraction.

      

(7.11)

      where S_i is the scatter estimated on all projections. Hutton and Baccarne (1998) proposed that S can be estimated once, based on an estimate of the activity distribution obtained from Eq. (3.12) in chapter 3, with the final scatter-corrected result obtained by proceeding with Eq. (7.11). It is worth to point out that most of the research performed in this field is related to SPECT imaging. Further development and validation of this class of algorithms in 3D PET are still required.

7.5.1 Monte Carlo-based scatter correction (MCBSC)

      The possibility of using Monte Carlo simulation to calculate and correct for Compton scatter in 3D PET was first suggested by Levin et al (1995). As proposed, the method uses as an estimate of the true distribution, the non-scatter corrected reconstructed 3D image volume for a simulated 3D PET acquisition, assuming that the number of events in each voxel accurately represents the number of annihilation photons originating at the corresponding location within the simulated body. Monte Carlo simulations requires as input the voxelised source distribution as well as the corresponding distribution of attenuation coefficients. The contours and attenuation coefficients of the different organs and tissues are obtained from a segmented transmission image derived from a short transmission scan.

      The sinograms of simulated scattered and unscattered events are created and used to estimate the scatter distribution for any given plane. For obvious reasons, the Monte Carlo generated data sets have a lower number of events than the measured data sets. Scaling of the simulated scatter distribution is therefore required to normalise the statistics of the simulated and the measured data. This scale factor was calculated as the ratio between the uncorrected emission sinograms and the total simulated emission sinogram (unscattered + scattered). Correction for scatter is performed by subtracting the estimated scatter distribution, , multiplied by the calculated scale factor from the measured emission data set according to:

      

(7.12)

      where k is the ratio of total counts in the acquired data to the total counts in the simulated data. To reduce computation time, the authors proposed to use coarse sampling of the input emission and segmented transmission images. That is, a 128x128 image volume is resampled to 32x32 by grouping 16 pixels (4x4) in a single pixel. As proposed by Levin et al (1995), the reconstructed emission image used as input to the simulations was not pre-corrected for scatter. This results in an overestimation of the amount of scatter and blurring in the emission image due to miss-positioning of events. This standard method is refereed as MCBSC1 in this section. The authors claim that the negative systematic errors introduced by the inclusion of scatter in the input image are less relevant than the statistical noise in the data sets. Fig. 7.10 illustrates the basic principles of the standard method (MCBSC1).

      

Fig. 7.10. : Flow-chart illustrating the basic principles of the original Monte Carlo-based scatter correction method (MCBSC1)

      The limitations of the technique are also mentioned (in bold).

      The contribution of scatter from out-of-FOV activity is particularly significant in 3D whole body scanning. To take this effect into account, activity from outside the FOV should be incorporated in the calculation of the scatter distribution. For this purpose, data from several bed positions should be made available as input to the Monte Carlo simulation. The effect of pre-correcting the projection data for scatter using an approximate method (e.g. CVS) was also investigated in this work. The data sets were pre-corrected for scatter with the CVS method and the reconstructed images are then used as input to the Monte Carlo simulator. This approach seems reasonable from the point-of-view of the physics of PET imaging and is refereed in what follows as MCBSC2 which can be distinguished from MCBSC1 by the modifications mentioned above. Fig. 7.11 illustrates the basic principles of the modified version of the algorithm (MCBSC2).

      

Fig. 7.11. : Flow-chart illustrating the principles of the modified Monte Carlo-based scatter correction method (MCBSC2)

      The modifications and improvements brought to the technique are outlined (in bold).

      The effects of pre-correcting the first image estimate for scatter on the estimation of scatter distribution is illustrated in Fig. 7.12 for a clinical brain study. A representative plane of the Monte Carlo simulated scatter sinograms using as input reconstructed images not corrected for scatter (left), and the same plane of Monte Carlo simulated sinograms using as input reconstructed images corrected for scatter (middle) are shown. Normalised differences between both sinogram are also shown (right). Integral profiles of the simulated low statistics scatter sinogram planes showing small differences between the two distributions are also illustrated in Fig. 7.13.

      

Fig. 7.12. : Clinical cerebral study showing a representative plane of the Monte Carlo simulated scatter sinograms using as input reconstructed images not corrected for scatter (left), and the same plane of Monte Carlo simulated sinograms using as input reconstructed images corrected for scatter (middle)

      Normalised differences between both sinograms are also shown (right).

      

Fig. 7.13. : Integral profiles of the simulated scatter sinogram planes shown in Fig. 7.12 illustrating both non-corrected (solid line) and scatter corrected (dotted line) data

7.5.2 Statistical reconstruction-based scatter correction (SRBSC)

      Reconstruction-based scatter compensation is a technique in which the scatter response function is modelled during the reconstruction process. Several research groups devised different algorithms belonging to this class of methods (Floyd et al 1986, Bowsher and Floyd 1991, Beekman et al 1997b). In this section, a new technique for scatter correction in 3D PET called statistical reconstruction-based scatter correction (SRBSC) is proposed (Zaidi 2000a). The principle of the method is based on the hypothesis that the image corresponding to scattered events in the projection data consist of almost low-frequency components of activity distribution and that the low-frequency components will converge faster than the high-frequency ones in successive iterations of statistical reconstruction methods such as ML-EM. Other investigators independently reported a related method for scatter correction in SPECT imaging (Liu et al 1999).

      Iterative algorithms possess a non-uniform convergence property. Low-frequency components of the image tend to be recovered earlier in iterative reconstruction than high-frequency components. There has been some evidence presented by other researchers supporting this hypothesis (Pan and Yagle 1991). The study of convergence properties of the ML-EM algorithm by Fourier analysis revealed clearly the non-uniform frequency response of the EM algorithm (Tanaka 1987). Moreover, preliminary investigations of inverse Monte Carlo-based reconstruction indicate that the recovery of spatial frequency information is achieved at different numbers of iterations for different spatial frequencies: higher spatial frequencies appear at higher iterations while the lower frequencies (smooth structures) are well defined at early iterations (Floyd et al 1986).

      The distribution of scattered photons in nuclear imaging (SPECT and PET) has been studied extensively (Floyd et al 1984, Frey and Tsui 1990, Ljungberg and Strand 1990, Barney et al 1991, Bentourkia et al 1995). In terms of the frequency response, the scatter components of PET projection data tend to be dominated by low-frequency information, though there is some higher frequency information present. The SRBSC approach takes advantage of this by estimating the scatter component from forward projection of images reconstructed in early iterations of OSEM.

      Following the additive imaging model presented in section 7.4.2 and neglecting statistical noise, the recorded projection data can be divided in two components: unscattered (p_u) and scattered (p_s). The total response function of the scanner can therefore be divided in two response kernels corresponding to the scattered and unscattered components, srf and urf, respectively.

      

(7.13)

      where fwdproject is the forward projection operator. Within the limits of our assumptions, the activity distribution, f, can be roughly divided in two parts:

      

(7.14)

      where f_L denotes the low-frequency image and f_H the high-frequency one. The range of low-frequency components is assumed to be equivalent to the frequency range corresponding to the scatter component p_s as detailed in the following derivation of the method.

      It has been shown in the preceding section that the scatter response kernels extend a wide region outside the source activity, i.e. contain low-frequency components while the unscatter response kernels are limited to a small region corresponding to the location of the source distribution. Thus, the scatter components correspond to the low-frequency range in the source distribution (fig. 7.5). Based on the assumption that the high-frequency components will be smeared, i.e. filtered by the scatter response kernels, the scatter components in the projection data can be approximated as follows:

      

(7.15)

      The proposed SRBSC method for scatter correction exploits the properties mentioned above. Fig. 7.14 shows a flow-chart of the general principles of the method. The basic steps followed when applying the method consist of the following:

• estimate the low-frequency components by only one OSEM iteration;
• obtain the scatter components by forward projection of the estimated image convolved with Monte Carlo simulated scatter response kernel, srf;
• subtract the estimated scatter components from measured projections;
• reconstruct the image using any available reconstruction algorithm (analytic, iterative) using scatter corrected projections.

      

Fig. 7.14. : Flow-chart illustrating the general principles of statistical reconstruction-based scatter correction (SRBSC) technique

      The scatter distribution is modelled as the convolution of the forward projected data with the scatter response function on the 2D projections. As described in the CVS method (Bailey and Meikle 1994), the correction method is implemented on 2D projections corresponding to polar angle q » 0. These are formed from acquired sinograms with ring differences equal to 0 or unity (dr=0, ±1). The estimated scatter distribution is then used to correct any oblique sinogram in the 3D data set coplanar with a direct plane on the axial axis (Z) of the scanner. It is assumed that the scatter distribution do not vary significantly over a small range of q. The immediate consequence is that the computation time is significantly decreased since much fewer scatter projections need to be estimated.

      To confirm the fast convergence properties of the low-frequency components in statistical reconstructions, a clinical oncology study was reconstructed with 1, 2, 3 and 10 iterations of OSEM. This study was reconstructed by pure OSEM without inter-update or post-filtering. Fig. 7.15 illustrates a slice reconstructed with different numbers of iterations. It can be noticed that the quality of the images is degraded when the number of iterations is increased up to 10 where the ML-EM solution starts to diverge. This is a well known phenomena in iterative reconstruction and may be explained by the fact that OSEM is based on convergence of the reprojected data of estimated images to the measured noisy projections. Therefore, the reconstructions will converge to noisy images and may diverge from the solution image when we increase the number of iterations.

      

Fig. 7.15. : OSEM reconstructions of an PRT-1 oncology study at the level of the thorax using different numbers of iterations: they are from top left, clockwise: one full iteration, 2 iterations, 3 iterations, and 10 iterations, respectively

      The scatter components converge very rapidly in the first iterations of OSEM. Fig. 7.16 shows the profiles of scatter distribution in one sinogram plane as estimated by the SRBSC method after 1, 2 and 3 full iterations of OSEM. Small but no appreciable differences are visible on the curves. Similar results have been obtained on other projection planes. This indicates that the estimated scatter components converge in the first iterations of OSEM and has some interesting practical consequences, the most important being that one iteration of OSEM is sufficient to assess the distribution of low-frequency scatter. This property makes the method very fast which renders its implementation in clinical routine viable.

      

Fig. 7.16. : Profile of the scatter component in a projection plane of the clinical oncology study estimated by SRBSC with one (solid line), two (dotted line) and three full iterations (dashed line) of OSEM

7.6.1 Phantom simulation studies

      Previously reported Monte Carlo validations of scatter correction techniques and related parameters have often been made by using simple source and attenuating medium geometries and compositions. In addition to these evaluations, there is also need to investigate more clinically realistic source distributions to validate and compare scatter correction techniques. A comparison of the relative performance of different methods are presented in this section using Monte Carlo simulated data of a digital brain phantom (Hoffman et al 1990). Evaluations have also been performed for the Utah phantom to characterise the differences between the methods for more standard geometries and allow determination of whether performance in the standard geometry predicts performance in the clinically realistic source distribution.

      Different regions of interest (ROIs) were defined and the average number of events within each ROI computed for both ideal, non-corrected and scatter corrected images. With a few exceptions, most papers dealing with the evaluation of scatter correction techniques compare relative concentrations within different compartments of a given phantom with the background compartment serving as a reference. This approach is criticised because it and does not necessarily reflect the accuracy of the correction procedure and might bias the evaluation procedure. Hence, the evaluation reported here was performed in absolute terms. Contrast recovery can be viewed as an example of parameter estimation. It is defined as:

      

(7.16)

      where N_BG is the count density in the ROI corresponding to the background and N_cold is the count density in the ROI corresponding to the cold cylinder. Obviously, the ideal value is 100%. The activity recovery coefficients, defined as the percentage ratio of events in a ROI in the corrected images to events in the same ROI in the ideal images, were calculated according to:

      

(7.17)

      where N_ROI is the count density in the ROI. The signal-to-noise ratio was defined as the mean number of events divided by the standard deviation of pixel intensities in a ROI defined within the background region A of the Utah phantom (Fig. 7.17).

      

(7.18)

---

---

7.6.2 Experimental phantom studies

      Experimental phantom studies were performed on the ECAT 953B scanner operated in 3D mode. Both the Hoffman 3D brain and Utah phantoms were scanned and analysed according to the same protocol used for the simulated data sets.

---

---

7.6.3 Clinical studies

      The different scatter correction algorithms were also tested on clinical data obtained on the PRT-1 scanner from the Geneva University Hospital (Townsend et al 1993). A cerebral and whole-body clinical study were selected from the database and used for clinical evaluation of the scatter correction methods.

---

---