- Email: [email protected]
Development of a Compton camera for prompt-gamma medical imaging
Radiation Physics and Chemistry (xxxx) xxxx–xxxx
Contents lists available at ScienceDirect
Radiation Physics and Chemistry journal homepage: www.elsevier.com/locate/radphyschem
Development of a Compton camera for prompt-gamma medical imaging ⁎
S. Aldawooda,b, P.G. Thirolfa, , A. Miania,c, M. Böhmerd, G. Dedesa, R. Gernhäuserd, C. Langa, S. Liprandia, L. Maierd, T. Marinšeka, M. Mayerhofera, D.R. Schaarte, I. Valencia Lozanoa, K. Parodia a
Fakultät für Physik, Department f. Medical Physics, Ludwig-Maximilians-Universität München, Am Coulombwall 1, Garching b. München, Germany Department of Physics and Astronomy, King Saud University, Riyadh, Saudi Arabia c Universitá degli Studi di Milano, Italy d Technische Universität München, Garching b. München, Germany e Delft University of Technology, Delft, The Netherlands b
A R T I C L E I N F O
A BS T RAC T
Keywords: Compton camera Prompt-γ imaging Hadron therapy Scintillation detectors
A Compton camera-based detector system for photon detection from nuclear reactions induced by proton (or heavier ion) beams is under development at LMU Munich, targeting the online range verification of the particle beam in hadron therapy via prompt-gamma imaging. The detector is designed to be capable to reconstruct the photon source origin not only from the Compton scattering kinematics of the primary photon, but also to allow for tracking of the secondary Compton-scattered electrons, thus enabling a γ-source reconstruction also from incompletely absorbed photon events. The Compton camera consists of a monolithic LaBr3:Ce scintillation crystal, read out by a multi-anode PMT acting as absorber, preceded by a stacked array of 6 double-sided silicon strip detectors as scatterers. The detector components have been characterized both under offline and online conditions. The LaBr3:Ce crystal exhibits an excellent time and energy resolution. Using intense collimated 137 Cs and 60Co sources, the monolithic scintillator was scanned on a fine 2D grid to generate a reference library of light amplitude distributions that allows for reconstructing the photon interaction position using a k-Nearest Neighbour (k-NN) algorithm. Systematic studies were performed to investigate the performance of the reconstruction algorithm, revealing an improvement of the spatial resolution with increasing photon energy to an optimum value of 3.7(1)mm at 1.33 MeV, achieved with the Categorical Average Pattern (CAP) modification of the k-NN algorithm.
1. Introduction The clinical application of tumor therapy via proton or heavy-ion beams has largely expanded over the last two decades. It exploits the favorable dose delivery properties of charged hadrons compared to conventional photon radiotherapy and major advancements in the fields of accelerator technology as well as treatment planning and diagnostic capabilities. The well-localized Bragg peak of charged particles in matter allows for a highly conformal dose deposition in the tumor, especially effective for the treatment of tumors in the vicinity of critical organs at risk. Moreover, sparing healthy tissue from unnecessary dose deposition is another benefit of particle therapy. However, while the number of clinical particle therapy facilities is still rapidly increasing (COG16article), exploiting the full benefits of the well-localized Bragg peak without applying large safety margins to account for methodological range uncertainties necessitates a reliable
⁎
monitoring and in-vivo verification of the ion-beam stopping range (Knopf and Lomax, 2013). Several experimental signatures are being investigated in view of their potential to provide precise in-vivo information on the Bragg peak position: prompt photons or delayed photons (exploiting their energy and timing properties) (Min, 2006; Golnik, 2014), secondary ions (Gwosch, 2013) or ion-induced ultrasonic shock waves (Assmann, 2015). The direct detection of promptly ( < ns) emitted γ radiation from nuclear reactions induced by the therapeutic proton or ion beam within the patient constitutes a promising option, since the distribution of prompt photons will not be blurred by physiological effects and no “wash-out” processes will occur, unlike in PET monitoring. The perspectives of “prompt-γ”-based medical imaging have been intensively studied in recent years, starting with feasibility studies to determine the correlation between the prompt γ radiation and the dose profile for mono-energetic proton beams (Min, 2006) and carbon beams (Testa, 2010). Several groups
Corresponding author. E-mail address: [email protected] (P.G. Thirolf).
http://dx.doi.org/10.1016/j.radphyschem.2017.01.024 Received 14 September 2016; Accepted 25 January 2017 0969-806X/ © 2017 Elsevier Ltd. All rights reserved.
Please cite this article as: Thirolf, P.G., Radiation Physics and Chemistry (2017), http://dx.doi.org/10.1016/j.radphyschem.2017.01.024
Radiation Physics and Chemistry (xxxx) xxxx–xxxx
S. Aldawood et al.
Fig. 2. Sketch of the layout of the LMU Compton camera prototype.
stacked array of six double sided silicon strip detectors (DSSSD), each of them with a thickness of 500 µm and an active area of 50×50 mm2. Each detector layer is segmented into 128 strips on the n and p sides, respectively, with a pitch size of 390 µm. The 256 silicon strips per detector are read out via 64-pin connectors at each of the four detector sides. Signals are transferred via AC coupling boards to the subsequent compact front-end signal processing boards, based on the GASSIPLEX (Santiard, 1994) ASIC chips. The signals processed in the six front-end boards of one detector side are then collected in a bus card and guided to a VME-based readout controller. Thus finally all 1536 DSSSD signals are fed into the combined event data stream together with the data from the absorber detector. The absorber component of the Compton camera is a monolithic 50×50×30 mm3 LaBr3:Ce3+ scintillator crystal (Gobain). This scintillator material offers favorable properties in all relevant quantities for γray detection. LaBr3:Ce3+ outperforms other commonly used scintillation materials in terms of timing and energy resolution. It exhibits a very fast decay time τ=17 ns, very high light yield (63,000 photons/ MeV), good energy resolution and reasonably high mass density (ρ=5.07 g/cm3). Similar to all lanthanum halide crystals, also LaBr3:Ce features an unavoidable amount of internal radioactivity, due to contaminations from 227Ac and its daughter products, as well as the presence of a small amount of the radioactive isotope 138La. In case of our absorber crystal, an internal activity of 2 Bq/cm3 (in total 150 Bq) was measured. The scintillation crystal is read out by a multi-anode photomultiplier tube (PMT) with 16×16 segments (each 3×3 mm2, Hamamatsu H9500). In addition to the 256 individual channels of the PMT, a sum signal can be extracted via the ‘sum dynode’ output. The 257 signals are fed into 16channel amplifier and Constant Fraction modules (mesytec MCFD-16), generating both amplified charge signals and amplitude-independent timing gates. Subsequently, the charge signals are fed into VME-based Charge-to-Digital (QDC) converter modules (mesytec MQDC-32), while time-of-flight measurements (relative to an external trigger) are enabled by sending the gate signals to Time-to-Digital (TDC) modules (mesytec MTDC-32). Finally, data acquisition of the combined data stream from the two Compton camera detector components is achieved through a frontend CPU (Power-PC RIO-3, operated with the LynxOS real-time operating system) and the MBS- and ROOT-based acquisition and analysis system MARABOU (Lutter, 2000). The scintillation detector was characterized in the laboratory with calibration sources, revealing a position independent relative energy resolution of ΔE/E =3.5(2)% and an excellent time resolution of 273(6) ps (FWHM) (Thirolf, 2014; Marinšek, 2015).
Fig. 1. Scheme of the principle of Compton camera operation in electron tracking mode.
investigate the possibilities for a range verification and in-vivo dosimetry via prompt γ radiation from nuclear reactions, either employing passively collimated imaging devices (Min et al., 2012; Roellinghoff, 2014; Perali, 2014) or electronically collimated systems based on the Compton-camera principle (Kormoll, 2011; Richard, 2011; Llosa, 2013; Polf, 2015; Thirolf, 2016). In general, a Compton camera consists of a scatterer and an absorber component. While in the conventional design of a Compton camera exclusively the photon interaction is registered, an advanced approach, pursued in the LMU Compton camera design, targets the tracking of the Compton electrons as well. Fig. 1 displays this ‘electron-tracking’ mode of operation, which is characterized by a stack of position-sensitive scatter detectors forming a tracking array (instead of a monolithic scatter detector), where the thickness of the individual detectors has to be chosen thin enough to allow the Compton scattered electrons to penetrate at least 2–3 layers without too much deflection by Moliére scattering, thus enabling their trajectory reconstruction. The primary γ ray with energy E1 interacts with one of the scatter detectors, depositing a fraction of its initial energy, while being deflected by the Compton scattering angle θ. If the remaining photon energy E2 is low enough, the scattered photon is fully absorbed by the second component of the Compton camera. The Compton scattering angle θ can be inferred by detecting the deposited energies and the interaction positions of the primary and the Compton-scattered photon in the two detector components, exploiting the Compton scattering kinematics based on energy and momentum conservation. Thus the primary photon origin can be constrained to the surface of the Compton cone, spanned by θ with its apex given by the primary interaction position in the scatter detector. From the intersection of different Compton cones, inferred from subsequent photon interactions originating from the same source, the γ-ray source position can be determined. In addition to this conventional source reconstruction scheme, the kinematical information carried by the Compton electron trajectory can be exploited to derive independent spatial information, which restricts the Compton cone to an arc segment (as indicated in Fig. 1). This mode of operation provides an increased reconstruction efficiency, since it allows to reconstruct the photon source position also from incompletely absorbed events, i.e. events, where part of the scattered photon energy escapes the absorption in the second component of the Compton camera, e.g. by Compton scattering or pair creation with subsequent single- or double escape of the 511 keV annihilation photons. 1.1. The Compton camera prototype
2. Determination of the photon interaction position in a monolithic scintillator
The geometrical layout of the LMU Compton camera prototype is schematically shown in Fig. 2. The scatter component is formed by a
The determination of the interaction position of Compton-scattered multi-MeV photons in the absorbing scintillator of the Compton 2
Radiation Physics and Chemistry (xxxx) xxxx–xxxx
S. Aldawood et al.
Fig. 3. Flowchart of the k-NN algorithm (for details see main text).
systematic studies with (collimated) 137Cs and 60Co calibration sources. The basic procedure is always to compare an unknown 2D light amplitude distribution to a set of reference light amplitude distributions previously recorded. This reference library consists of nepp ·n pos light distributions measured by scanning the scintillation detector front surface (spanning the x-y plane) with a perpendicularly impinging, tightly collimated γ source for a total of n pos irradiation positions in the x-y plane and nepp recorded photopeak events per irradiation position. Thus every reference light distribution via its x-y coordinates carries the spatial information of the photon entry position to the detector front surface. The standard version of the position reconstruction algorithm is called the “k-Nearest Neighbours” (k-NN) algorithm. The light distribution of the spatially unknown photon in the Compton camera absorber is compared to all the entries of the reference library by calculating the intensity difference DlkNN , defined as
camera is a mandatory prerequisite for reconstructing the photon source position. The determination of the photon interaction position in a monolithic scintillator is not as straightforward as in a detector based on segmented (and optically isolated) crystals. Due to internal reflections of the scintillation light, the readout segment with the highest light amplitude does not necessarily correspond to the primary photon interaction position. Moreover, scintillation light will typically spread out over several readout sensor segments, which increases the statistical variance of the light yield per segment and increases the influence of background and electronic noise. However, monolithic crystals avoid dead volume, exhibit a favorable high-energy sensitivity and facilitate the detector assembly.
2.1. The k-nearest neighbours (k-NN) algorithm Targeting the application of monolithic scintillators in PET scanner systems, position reconstruction algorithms have been developed at TU Delft to determine the entry points of 511 keV annihilation photons into the crystal front surface (van Dam, 2011). Here, the position determination is realized by comparing the 2-dimensional light amplitude distribution of the impinging photon measured by a segmented photosensor (multi-anode photomultiplier tube, (digital) silicon photomultiplier array, avalanche photodiode array) with sets of reference data containing a large number of events recorded at known entry points. This statistical method is known as k-Nearest Neighbour (kNN) algorithm, which was first introduced by Fix and Hodges in 1951 (Fix and Hodges, 1951) and later on applied to the determination of the entry point of γ rays to monolithic scintillation detectors by Maas (2006). At TU Delft extensive studies have been performed to optimize the method in view of its performance and computational as well as measurement requirements. In order to validate the applicability of the k-NN algorithm beyond the energy range of PET also for prompt-γ medical imaging with multiMeV photons and, in particular, when using the monolithic LaBr3:Ce absorbing scintillator of our Compton camera prototype, we performed
N
DlkNN =
∑ (Iunknown,i − Iref(l,i))2 , i =1
l = 1, …, (n pos·nepp )
(1)
where N is the total number of photosensor pixels (in our case corresponding to the segmentation of the multi-anode PMT), Iunknown,i is the light intensity of the i-th pixel of the unknown light distribution and Iref(l,i) is the intensity of the i-th pixel of the l-th reference light distribution. Among all the reference light distributions, the k closest matching 2D distributions, i.e. the kkNN distributions with smallest DlkNN , are selected to create the so-called k-NN histogram that contains the x and y coordinates of the kkNN distributions. This histogram is then smoothed with a moving average filter, which for each pixel considers a surrounding matrix of 5×5 pixel and assigns to the central pixel the mean value to the surroundings. The x-y coordinates of the maximum value of the smoothed k-NN histogram are then identified with the coordinates of the calculated interaction position. The top panel of Fig. 3 displays the flowchart of this reconstruction method. 3
Radiation Physics and Chemistry (xxxx) xxxx–xxxx
S. Aldawood et al.
respective reconstruction algorithm to each one of the n pos xnepp library entries. Thus successively each of the 2D light amplitude distributions, which form the reference library, is considered as “unknown”, while the rest of the library entries constitutes the reference data set. Once the reconstruction has been performed and the position of the primary photon interaction in the scintillator has been calculated, the differences of the x and y coordinates of the calculated position and the “true” position of the respective reference library entry assumed to be unknown, (Δx, Δy ) = (x,calc − xtrue , y,calc − ytrue ), are computed. These differences are then filled into the 2D ‘error histogram’, as indicated in Figs. 3 and 4. Finally, the FWHM of the error histogram (i.e. the average of the FWHM values of its projections on x and y direction) is identified with the spatial resolution achieved in the specific reconstruction scenario.
2.2. The categorical average pattern (CAP) algorithm The second method used to determine the interaction position in the absorber of the Compton camera is the Categorical Average Pattern (CAP) algorithm, a modification of the standard k-NN method (van Dam, 2011). The main difference between the standard k-NN and CAP is that in the latter the k closest matching light distributions are searched amongst the nepp reference light distributions collected independently for each irradiation position. From this subset of k CAP reference distributions an average 2D light distribution is then calculated. The difference DlCAP is then calculated between the average light distribution and the unknown one for each irradiation position according to N
DlCAP =
∑ (Iunknown,i − Iave(l,i))2 , i =1
l = 1, …, (n pos)
(2)
2.4. The reference libraries
where N is again the total number of photosensor pixels, Iunknown,i is the light intensity of the i-th pixel of the unknown event and Iave(l,i) is the intensity of the i-th pixel of the l-th average light distribution. The outcome is then a 2D plot containing n pos intensity differences, each one corresponding to an (x,y) pair of irradiation coordinates. This plot is then smoothed with a 5×5 moving average filter and the coordinates of the minimum value of the plot are assigned as the calculated (x,y) interaction position inside the scintillator crystal. In Fig. 4 the flowchart of the CAP algorithm is shown.
In the process of acquiring the light amplitude reference libraries discussed before, a sequence of five correction steps has to be applied consecutively to the raw data, which will be introduced in the following (Aldawood, 2015): (i) gain matching:Since each PMT segment is read out and processed by an individual electronics chain, gain variations may occur. Injecting two pulser signals of different amplitude into the signal processing chain, gain correction factors can be derived for each signal channel. (ii) QDC pedestal subtraction: Unavoidable electronic dark current integrated and digitized in the QDC modules (‘energy pedestal’) can be quantified in a dedicated measurement without detector input. The resulting pedestal peaks are then taken as lower energy threshold for the data acquisition process. (iii) PMT nonuniformity: Possible gain fluctuations of the multi-anode PMT segments have to be accounted for via the non-uniformity matrix provided
2.3. Determination of the spatial resolution: “Leave-One-Out” method In order to determine the spatial resolution achievable via the standard k-NN or CAP photon-interaction position reconstruction algorithms, the “Leave-One-Out”-method was applied. As schematically indicated in Figs. 3 and 4, this method is based on applying the
Fig. 4. Flowchart of the CAP algorithm (for details see main text).
4
Radiation Physics and Chemistry (xxxx) xxxx–xxxx
S. Aldawood et al.
by the manufacturer. (iv) Spatial homogeneity: Photon scattering and reflections at the side surface wrapping will induce anisotropies in the absorption of isotropic incident radiation. This can be corrected by making use of the internal radioactivity of the LaBr3:Ce scintillator material, safely assumed to be homogeneously distributed in the crystal volume. The internal activity energy spectrum is acquired and a position-dependent correction matrix is derived after gating on the energy region of interest (i.e. the energy of the collimated γ-ray source). (v) Energy gating: An energy gate in the region of the photopeak energy of interest is applied. After preparing the data acquisition as described before, the 2D detector scan can be performed. The intense calibration sources (137Cs/ 82 MBq and 60Co/20 MBq) are placed behind a Densimet collimator (ρ=18.5 g/cm3) with a diameter of 1 mm. In case of 137Cs with a photon energy of 662 keV, the collimator length was 48 mm and the source and its surrounding Pb shielding were placed on a motorized and automatic controlled x-y translation stage in front of the stationary scintillation detector. For 60Co, the higher photon energies (1.17/ 1.33 MeV) required an extended collimator length of 100 mm and a thicker Pb side shielding. Therefore, the setup was modified as indicated in Fig. 5. Now the shielded source/collimator unit was kept stationary, while the detector was moved in front of the collimator opening. In this setup, the collimator was formed by a 100×100×100 mm3 Densimet block with a central bore of 4 mm diameter, allowing to insert an exchangeable collimator tube from sintered WC (ρ=15.6 g/cm3) with an open diameter of 1 mm. In a preparatory measurement series, the position of the LaBr3:Ce crystal inside its aluminum housing was determined by scanning across the detector edges in x and y direction. The rise and fall of the detected photon count rate indicated the crystal position and served to determine the starting and ending points of the subsequent 2D scans. In Fig. 6 (top panel) a subset of the 2D reference library acquired with a collimated (1 mm opening) 137Cs (82 MBq) is shown, representing a scan of 16×16 irradiation positions with a step size of 3 mm in x and y direction. The full scan has been carried out using a fine grid of 0.5 mm step size and 103×103 irradiation positions. For each position, data were acquired during 30 s, allowing to collect at least 400 photopeak events per position. The bottom panel of Fig. 6 displays the corresponding library subset when using a 60Co calibration source (20 MBq) and after gating on the 1.33 MeV photopeak. The correlation of the light amplitude distributions with the source position is clearly visible.
Fig. 6. Subsets of the 2D light amplitude distribution libraries generated with a 137Cs (top) and 60Co source gated at 1.3bsets. 16×16 irradiation positions (step size 3 mm in x and y direction) are shown, measured with a 1 mm collimator opening.
2.5. Systematic studies of the k-NN performance Both position reconstruction algorithms, standard k-NN and CAP, were systematically studied in order to identify the optimum values of the key parameters k (number of best matching reference light distributions) and nepp (number of photopeak reference events per irradiation position). Moreover, the analysis of the achievable spatial resolution was performed for two values of the PMT segmentation, using either the full 256-fold segmentation provided by the electronical readout or combining the signals from each neighbouring four segments to realize a 64-fold segmentation a posteriori by software in the offline analysis. The study was performed for the three available photon energies from 137Cs and 60Co and for two values of the scan step size (0.5 mm, 1 mm) (Miani, 2016). Fig. 7 displays the resulting values of the spatial resolution as extracted by using the standard k-NN algorithm and a reference library with nepp = 400 photopeak events acquired per irradiation position. Clearly visible is the improved spatial resolution for the higher photon energies provided by the 60Co source compared to the curves derived at 662 keV from 137Cs. The minimum of the spatial resolution is consistently reached for values of k around 1000. The best spatial resolution is reached for a PMT granularity of 64 channels and a scan step size of 0.5 mm. Results from a similar study using the CAP algorithm are shown in
Fig. 5. Sketch of the 2D scanning system used to position the LaBr3:Ce scintillator in front of a strong, tightly collimated γ-ray source. The top inset shows photographs of the collimator block with the exchangeable sintered WC collimator tube.
5
Radiation Physics and Chemistry (xxxx) xxxx–xxxx
S. Aldawood et al.
Fig. 9. Determination of the experimental uncertainties of the spatial resolution data, obtained with the CAP algorithm at 1.33 MeV (details: see text).
Fig. 7. Spatial resolution as a function of k (number of nearest neighbours) achieved with the k-NN algorithm, using nepp =400. The 3 parameters photon energy, PMT granularity and scan pitch are distinguished by their line style, marker shape and symbol filling style, respectively.
of the number of best matching light amplitude distributions. The spatial resolution obtained with the CAP algorithm at 1.33 MeV is shown as a function of k (number of nearest neighbours) in Fig. 9. The data refer to the 0.5 mm pitch size libraries, created by using the four subsets with 100 reference events per irradiation position extracted from the full library (with 400 reference 2D light amplitude distributions per irradiation position), for both PMT granularities of 64 and 256, respectively. The index from 1 to 4 indicates the 4 library subsets. Consequently, for a fixed algorithm, PMT granularity value, pitch size and photon energy, the standard deviation of the four spatial resolution values obtained for the four data subsets was calculated and assigned as the uncertainty for this specific scenario. After having determined the optimum value of k for different irradiation energies and readout scenarios, systematic studies of the achievable spatial resolution (i.e. FWHM of the error histogram of the reconstruction algorithm) were performed as a function of the number of photopeak events per irradiation position, nepp . Fig. 10 illustrates the performance of the standard k-NN and CAP algorithms applied to different irradiation energies (662 keV, 1.17 MeV, 1.33 MeV), PMT granularities (64- and 256-fold) and scan pitch sizes (0.5 mm, 1 mm) as a function of nepp , each time applying the optimum value of k as determined before. The top panel shows the corresponding results for 662 keV from 137Cs, while the middle and bottom panels display the results for 1.17 MeV and 1.32 MeV from 60Co, respectively. Consistent trends can be inferred from the three different energy dependent data sets. In general, the spatial resolution improves with an increasing number of reference events per position, due to the larger statistics the unknown photon event can be compared to. Below nepp ≈ 50 , large fluctuations and uncertainties of the individual data sets prevent a clear identification of comparative trends. Above nepp ≈ 50 − 100 (depending on the energy), the CAP algorithm provides (within the experimental uncertainties) at least comparable and mostly even better results than the standard k-NN algorithm. Moreover, the 0.5 mm scan pitch size for a given PMT granularity results in a better spatial resolution compared to the 1 mm case. Important from a practical point of view is the finding that a lower PMT granularity (64 segments) consistently results in better spatial resolution compared to the higher (256 fold) granularity. This behaviour can most likely be attributed to the higher photon statistics per readout pixel in case of the lower granularity, even though the summing of each four pixels was not done a priori electronically in hardware, but a posteriori during the software analysis (van der Laan, 2010). The most important difference of the energy dependent data sets in Fig. 10 is the consistently improving spatial resolution as the photon energy increases.
Fig. 8. Again reference libraries with 400 photopeak events per irradiation position were used. A general feature of all curves is the optimum value of k ≈ 12 (remind that for the CAP algorithm the best matching reference library entries are searched within the ensemble of each irradiation position). While the general trends visible in the CAP data are similar to the ones discussed before for the standard k-NN algorithm, the optimum value of the spatial resolution is consistently lower for CAP compared to k-NN. The experimental uncertainties included in the data shown in Figs. 7 and 8 have been estimated as fluctuations of the spatial resolution between different realizations of the ‘leave-one-out’ procedure applied to the various parameter scenarios. In particular, the initial reference library with 400 photopeak events per irradiation position was subdivided into four subsets with 100 photopeak events each before execution of the ‘leave-one-out’ method for various values k
Fig. 8. Spatial resolution as a function of k (number of nearest neighbours) achieved with the CAP algorithm, using nepp =400. The 3 parameters photon energy, PMT granularity and scan pitch are distinguished by their line style, marker shape and symbol filling style, respectively.
6
Radiation Physics and Chemistry (xxxx) xxxx–xxxx
S. Aldawood et al.
Fig. 11. Energy dependence of the spatial resolution, comparing the two reconstruction algorithms k − NN and CAP in different scenarios. While the scan pitch size was kept constant at 0.5 mm, the PMT granularity was varied between a 64-and 256 fold segmentation of the scintillator PMT readout. The data were acquired with a 1 mm collimator and optimized values of k were used for nepp =400 photopeak events per irradiation position.
tion of the spatial resolution of its absorber component, a monolithic LaBr3:Ce scintillation crystal. The photon entry position at the crystal front surface was determined by using the k-NN reconstruction algorithm in its standard form as well as in the modified version of the Categorical Average Pattern (CAP) algorithm. The applicability of this method beyond its initial energy range of 511 keV annihilation photons (for PET applications) was verified up to 1.33 MeV, with a clear trend of improving spatial resolution as a function of the photon energy, as can be seen in Fig. 11. The plot also indicates the design goal of the Compton camera project, which aims at reaching a position reconstruction resolution in the absorbing scintillator of 3 mm at the multi-MeV photon energies to be expected from the prompt-γ spectrum of hadron-tissue interactions. According to design specification simulations (Lang, 2014), this resolution would enable an overall spatial resolution of the Compton camera of 1.5°–2° in the relevant photon energy range of 2–6 MeV (equivalent to 1.5–2 mm resolution for a distance of 50 mm between the photon source and the first scatter detector, as in the case of a small-animal irradiation scenario). In view of envisaged attempts to further improve the spatial resolution by using an even tighter γ-ray source collimation with a diameter of 0.6 mm, correspondingly reducing the photon flux by about a factor of 3, the saturation of the spatial resolution improvement for nepp ≥ 200 will help to keep the required measurement in a realistically manageable range of a few days. However, the ultimate goal is to commission the Compton camera in the relevant energy range of 4–6 MeV, not accessible via laboratory calibration sources. The energy-dependent trend of the absorber's spatial resolution established in this study nourishes confidence that the design goal can be met, taking into account that with higher photon energies also the contribution of pair creation will increase, giving rise to an even more localized interaction of the photons in the scintillation crystal. Experimental characterization in this energy range might be possible by exploiting the availability of well-collimated, highly monoenergetic photon beams from (nγ) reactions at neutron sources like the GAMS6 facility of the Institut Laue Langevin in Grenoble (Kessler, 2001), as currently being planned. Supported by DFG Cluster of Excellence MAP (Munich-Centre for Advanced Photonics) and King Saud University, Riyadh, Saudi Arabia.
Fig. 10. Spatial resolution as a function of nepp (number of photopeak events per irradiation position for Eγ = 662 keV (top), 1.17 MeV (middle) and 1.33 MeV (bottom), plotted for 3 parameters: type of reconstruction algorithm, PMT granularity and scan pitch size. According to the systematic study discussed before, the optimum value of k was chosen for each parameter value data set.
Apparently, the improvement due the the higher photon count dominates the deterioration of spatial resolution that might be expected from increased intra-crystal scattering (van der Laan, 2010; Maas, 2010). The best values of the spatial resolution determined via the CAP algorithm for nepp = 400 and 0.5 mm scan step size were 4.8(1) mm (662 keV), 3.9(1) mm (1.17 MeV) and 3.7(1) mm (1.33 MeV), respectively. 2.6. Conclusions and outlook The current study presents a key ingredient of the commissioning process of the Garching Compton camera prototype, designed for beam range monitoring in hadron therapy, namely the determina7
Radiation Physics and Chemistry (xxxx) xxxx–xxxx
S. Aldawood et al.
Llosa, G., et al., 2013. Nucl. Instrum. Methods A 718, 130. Lutter, R., et al., 2000. IEEE Trans. Nucl. Sci. 47, 280. Maas, M.C., et al., 2006. IEEE Trans. Nucl. Sci. 53, 1071. Maas, M.C., et al., 2010. Med. Phys. 37, 1904–1913. Marinšek, T., 2015. Thesis, Master., LMU Munich, unpublished. Miani, A., 2016. Thesis, Master, Univ. degli Studi di Milano, unpublished. Min, C.-H., et al., 2006. Appl. Phys. Lett. 89, 183517. Min, C.-H., Lee, H.R., Kim, C.H., Lee, S.B., 2012. Med. Phys. 39, 2100. Perali, I., et al., 2014. Phys. Med. Biol. 59, 5849. Polf, J.C., et al., 2015. Phys. Med. Biol. 60, 7085. Richard, M.-H., et al., 2011. IEEE Trans. Nucl. Sci. 58, 87. Roellinghoff, F., et al., 2014. Phys. Med. Biol. 59, 1327. Santiard, J.C., et al., 1994. CERN-ECP/94 17. Testa, M., et al., 2010. Radiat. Environ. Biophys. 49, 337. Thirolf, P.G., et al., 2016. EPJ Web of Conferences 117, 05005. Thirolf, P.G., et al., 2014. EPJ Web of Conferences 66, 11036. van Dam, H.T., et al., 2011. IEEE Trans. Nucl. Sci. 58, 2139.
References Aldawood, S., et al., 2015. Front. Oncol. 5, 270. Assmann, W., et al., 2015. Med. Phys. 42, 567. COG16article, Particle Therapy Co-operative Group,〈http://www.ptcog.ch〉 Fix, E., Hodges, J.L., 1951. Technical Report 4, USAF School of Aviation Medicine. Gobain, St., scintillation material data sheet, 〈http://www.crystals.saint-gobain.com/ uploadedFilesSG-CrystalsDocumentsBrilLanCe%20380%20Dat%20Sheet.pdf〉 Golnik, C., et al., 2014. Phys. Med. Biol. 59, 5399. Gwosch, K., et al., 2013. Phys. Med. Biol. 58, 3755. Kessler, E.G., et al., 2001. Nucl. Instrum. Methods A 457, 187. Knopf, A.-C., Lomax, A., 2013. Phys. Med. Biol. 58, R131. Kormoll, T., et al. 2012. Nucl. Instrum. Methods A, 11, 4, 2011, Erratum ibid. A 678, (13) 626-627. Krimmer, J., et al., 2015. Nucl. Instrum. Methods A 787, 98. van der Laan, D.J., et al., 2010. Phys. Med. Biol. 55, 1659–1675. Lang, C., et al., 2014. J. Instrum. 9, P01008.
8
Contents lists available at ScienceDirect
Radiation Physics and Chemistry journal homepage: www.elsevier.com/locate/radphyschem
Development of a Compton camera for prompt-gamma medical imaging ⁎
S. Aldawooda,b, P.G. Thirolfa, , A. Miania,c, M. Böhmerd, G. Dedesa, R. Gernhäuserd, C. Langa, S. Liprandia, L. Maierd, T. Marinšeka, M. Mayerhofera, D.R. Schaarte, I. Valencia Lozanoa, K. Parodia a
Fakultät für Physik, Department f. Medical Physics, Ludwig-Maximilians-Universität München, Am Coulombwall 1, Garching b. München, Germany Department of Physics and Astronomy, King Saud University, Riyadh, Saudi Arabia c Universitá degli Studi di Milano, Italy d Technische Universität München, Garching b. München, Germany e Delft University of Technology, Delft, The Netherlands b
A R T I C L E I N F O
A BS T RAC T
Keywords: Compton camera Prompt-γ imaging Hadron therapy Scintillation detectors
A Compton camera-based detector system for photon detection from nuclear reactions induced by proton (or heavier ion) beams is under development at LMU Munich, targeting the online range verification of the particle beam in hadron therapy via prompt-gamma imaging. The detector is designed to be capable to reconstruct the photon source origin not only from the Compton scattering kinematics of the primary photon, but also to allow for tracking of the secondary Compton-scattered electrons, thus enabling a γ-source reconstruction also from incompletely absorbed photon events. The Compton camera consists of a monolithic LaBr3:Ce scintillation crystal, read out by a multi-anode PMT acting as absorber, preceded by a stacked array of 6 double-sided silicon strip detectors as scatterers. The detector components have been characterized both under offline and online conditions. The LaBr3:Ce crystal exhibits an excellent time and energy resolution. Using intense collimated 137 Cs and 60Co sources, the monolithic scintillator was scanned on a fine 2D grid to generate a reference library of light amplitude distributions that allows for reconstructing the photon interaction position using a k-Nearest Neighbour (k-NN) algorithm. Systematic studies were performed to investigate the performance of the reconstruction algorithm, revealing an improvement of the spatial resolution with increasing photon energy to an optimum value of 3.7(1)mm at 1.33 MeV, achieved with the Categorical Average Pattern (CAP) modification of the k-NN algorithm.
1. Introduction The clinical application of tumor therapy via proton or heavy-ion beams has largely expanded over the last two decades. It exploits the favorable dose delivery properties of charged hadrons compared to conventional photon radiotherapy and major advancements in the fields of accelerator technology as well as treatment planning and diagnostic capabilities. The well-localized Bragg peak of charged particles in matter allows for a highly conformal dose deposition in the tumor, especially effective for the treatment of tumors in the vicinity of critical organs at risk. Moreover, sparing healthy tissue from unnecessary dose deposition is another benefit of particle therapy. However, while the number of clinical particle therapy facilities is still rapidly increasing (COG16article), exploiting the full benefits of the well-localized Bragg peak without applying large safety margins to account for methodological range uncertainties necessitates a reliable
⁎
monitoring and in-vivo verification of the ion-beam stopping range (Knopf and Lomax, 2013). Several experimental signatures are being investigated in view of their potential to provide precise in-vivo information on the Bragg peak position: prompt photons or delayed photons (exploiting their energy and timing properties) (Min, 2006; Golnik, 2014), secondary ions (Gwosch, 2013) or ion-induced ultrasonic shock waves (Assmann, 2015). The direct detection of promptly ( < ns) emitted γ radiation from nuclear reactions induced by the therapeutic proton or ion beam within the patient constitutes a promising option, since the distribution of prompt photons will not be blurred by physiological effects and no “wash-out” processes will occur, unlike in PET monitoring. The perspectives of “prompt-γ”-based medical imaging have been intensively studied in recent years, starting with feasibility studies to determine the correlation between the prompt γ radiation and the dose profile for mono-energetic proton beams (Min, 2006) and carbon beams (Testa, 2010). Several groups
Corresponding author. E-mail address: [email protected] (P.G. Thirolf).
http://dx.doi.org/10.1016/j.radphyschem.2017.01.024 Received 14 September 2016; Accepted 25 January 2017 0969-806X/ © 2017 Elsevier Ltd. All rights reserved.
Please cite this article as: Thirolf, P.G., Radiation Physics and Chemistry (2017), http://dx.doi.org/10.1016/j.radphyschem.2017.01.024
Radiation Physics and Chemistry (xxxx) xxxx–xxxx
S. Aldawood et al.
Fig. 2. Sketch of the layout of the LMU Compton camera prototype.
stacked array of six double sided silicon strip detectors (DSSSD), each of them with a thickness of 500 µm and an active area of 50×50 mm2. Each detector layer is segmented into 128 strips on the n and p sides, respectively, with a pitch size of 390 µm. The 256 silicon strips per detector are read out via 64-pin connectors at each of the four detector sides. Signals are transferred via AC coupling boards to the subsequent compact front-end signal processing boards, based on the GASSIPLEX (Santiard, 1994) ASIC chips. The signals processed in the six front-end boards of one detector side are then collected in a bus card and guided to a VME-based readout controller. Thus finally all 1536 DSSSD signals are fed into the combined event data stream together with the data from the absorber detector. The absorber component of the Compton camera is a monolithic 50×50×30 mm3 LaBr3:Ce3+ scintillator crystal (Gobain). This scintillator material offers favorable properties in all relevant quantities for γray detection. LaBr3:Ce3+ outperforms other commonly used scintillation materials in terms of timing and energy resolution. It exhibits a very fast decay time τ=17 ns, very high light yield (63,000 photons/ MeV), good energy resolution and reasonably high mass density (ρ=5.07 g/cm3). Similar to all lanthanum halide crystals, also LaBr3:Ce features an unavoidable amount of internal radioactivity, due to contaminations from 227Ac and its daughter products, as well as the presence of a small amount of the radioactive isotope 138La. In case of our absorber crystal, an internal activity of 2 Bq/cm3 (in total 150 Bq) was measured. The scintillation crystal is read out by a multi-anode photomultiplier tube (PMT) with 16×16 segments (each 3×3 mm2, Hamamatsu H9500). In addition to the 256 individual channels of the PMT, a sum signal can be extracted via the ‘sum dynode’ output. The 257 signals are fed into 16channel amplifier and Constant Fraction modules (mesytec MCFD-16), generating both amplified charge signals and amplitude-independent timing gates. Subsequently, the charge signals are fed into VME-based Charge-to-Digital (QDC) converter modules (mesytec MQDC-32), while time-of-flight measurements (relative to an external trigger) are enabled by sending the gate signals to Time-to-Digital (TDC) modules (mesytec MTDC-32). Finally, data acquisition of the combined data stream from the two Compton camera detector components is achieved through a frontend CPU (Power-PC RIO-3, operated with the LynxOS real-time operating system) and the MBS- and ROOT-based acquisition and analysis system MARABOU (Lutter, 2000). The scintillation detector was characterized in the laboratory with calibration sources, revealing a position independent relative energy resolution of ΔE/E =3.5(2)% and an excellent time resolution of 273(6) ps (FWHM) (Thirolf, 2014; Marinšek, 2015).
Fig. 1. Scheme of the principle of Compton camera operation in electron tracking mode.
investigate the possibilities for a range verification and in-vivo dosimetry via prompt γ radiation from nuclear reactions, either employing passively collimated imaging devices (Min et al., 2012; Roellinghoff, 2014; Perali, 2014) or electronically collimated systems based on the Compton-camera principle (Kormoll, 2011; Richard, 2011; Llosa, 2013; Polf, 2015; Thirolf, 2016). In general, a Compton camera consists of a scatterer and an absorber component. While in the conventional design of a Compton camera exclusively the photon interaction is registered, an advanced approach, pursued in the LMU Compton camera design, targets the tracking of the Compton electrons as well. Fig. 1 displays this ‘electron-tracking’ mode of operation, which is characterized by a stack of position-sensitive scatter detectors forming a tracking array (instead of a monolithic scatter detector), where the thickness of the individual detectors has to be chosen thin enough to allow the Compton scattered electrons to penetrate at least 2–3 layers without too much deflection by Moliére scattering, thus enabling their trajectory reconstruction. The primary γ ray with energy E1 interacts with one of the scatter detectors, depositing a fraction of its initial energy, while being deflected by the Compton scattering angle θ. If the remaining photon energy E2 is low enough, the scattered photon is fully absorbed by the second component of the Compton camera. The Compton scattering angle θ can be inferred by detecting the deposited energies and the interaction positions of the primary and the Compton-scattered photon in the two detector components, exploiting the Compton scattering kinematics based on energy and momentum conservation. Thus the primary photon origin can be constrained to the surface of the Compton cone, spanned by θ with its apex given by the primary interaction position in the scatter detector. From the intersection of different Compton cones, inferred from subsequent photon interactions originating from the same source, the γ-ray source position can be determined. In addition to this conventional source reconstruction scheme, the kinematical information carried by the Compton electron trajectory can be exploited to derive independent spatial information, which restricts the Compton cone to an arc segment (as indicated in Fig. 1). This mode of operation provides an increased reconstruction efficiency, since it allows to reconstruct the photon source position also from incompletely absorbed events, i.e. events, where part of the scattered photon energy escapes the absorption in the second component of the Compton camera, e.g. by Compton scattering or pair creation with subsequent single- or double escape of the 511 keV annihilation photons. 1.1. The Compton camera prototype
2. Determination of the photon interaction position in a monolithic scintillator
The geometrical layout of the LMU Compton camera prototype is schematically shown in Fig. 2. The scatter component is formed by a
The determination of the interaction position of Compton-scattered multi-MeV photons in the absorbing scintillator of the Compton 2
Radiation Physics and Chemistry (xxxx) xxxx–xxxx
S. Aldawood et al.
Fig. 3. Flowchart of the k-NN algorithm (for details see main text).
systematic studies with (collimated) 137Cs and 60Co calibration sources. The basic procedure is always to compare an unknown 2D light amplitude distribution to a set of reference light amplitude distributions previously recorded. This reference library consists of nepp ·n pos light distributions measured by scanning the scintillation detector front surface (spanning the x-y plane) with a perpendicularly impinging, tightly collimated γ source for a total of n pos irradiation positions in the x-y plane and nepp recorded photopeak events per irradiation position. Thus every reference light distribution via its x-y coordinates carries the spatial information of the photon entry position to the detector front surface. The standard version of the position reconstruction algorithm is called the “k-Nearest Neighbours” (k-NN) algorithm. The light distribution of the spatially unknown photon in the Compton camera absorber is compared to all the entries of the reference library by calculating the intensity difference DlkNN , defined as
camera is a mandatory prerequisite for reconstructing the photon source position. The determination of the photon interaction position in a monolithic scintillator is not as straightforward as in a detector based on segmented (and optically isolated) crystals. Due to internal reflections of the scintillation light, the readout segment with the highest light amplitude does not necessarily correspond to the primary photon interaction position. Moreover, scintillation light will typically spread out over several readout sensor segments, which increases the statistical variance of the light yield per segment and increases the influence of background and electronic noise. However, monolithic crystals avoid dead volume, exhibit a favorable high-energy sensitivity and facilitate the detector assembly.
2.1. The k-nearest neighbours (k-NN) algorithm Targeting the application of monolithic scintillators in PET scanner systems, position reconstruction algorithms have been developed at TU Delft to determine the entry points of 511 keV annihilation photons into the crystal front surface (van Dam, 2011). Here, the position determination is realized by comparing the 2-dimensional light amplitude distribution of the impinging photon measured by a segmented photosensor (multi-anode photomultiplier tube, (digital) silicon photomultiplier array, avalanche photodiode array) with sets of reference data containing a large number of events recorded at known entry points. This statistical method is known as k-Nearest Neighbour (kNN) algorithm, which was first introduced by Fix and Hodges in 1951 (Fix and Hodges, 1951) and later on applied to the determination of the entry point of γ rays to monolithic scintillation detectors by Maas (2006). At TU Delft extensive studies have been performed to optimize the method in view of its performance and computational as well as measurement requirements. In order to validate the applicability of the k-NN algorithm beyond the energy range of PET also for prompt-γ medical imaging with multiMeV photons and, in particular, when using the monolithic LaBr3:Ce absorbing scintillator of our Compton camera prototype, we performed
N
DlkNN =
∑ (Iunknown,i − Iref(l,i))2 , i =1
l = 1, …, (n pos·nepp )
(1)
where N is the total number of photosensor pixels (in our case corresponding to the segmentation of the multi-anode PMT), Iunknown,i is the light intensity of the i-th pixel of the unknown light distribution and Iref(l,i) is the intensity of the i-th pixel of the l-th reference light distribution. Among all the reference light distributions, the k closest matching 2D distributions, i.e. the kkNN distributions with smallest DlkNN , are selected to create the so-called k-NN histogram that contains the x and y coordinates of the kkNN distributions. This histogram is then smoothed with a moving average filter, which for each pixel considers a surrounding matrix of 5×5 pixel and assigns to the central pixel the mean value to the surroundings. The x-y coordinates of the maximum value of the smoothed k-NN histogram are then identified with the coordinates of the calculated interaction position. The top panel of Fig. 3 displays the flowchart of this reconstruction method. 3
Radiation Physics and Chemistry (xxxx) xxxx–xxxx
S. Aldawood et al.
respective reconstruction algorithm to each one of the n pos xnepp library entries. Thus successively each of the 2D light amplitude distributions, which form the reference library, is considered as “unknown”, while the rest of the library entries constitutes the reference data set. Once the reconstruction has been performed and the position of the primary photon interaction in the scintillator has been calculated, the differences of the x and y coordinates of the calculated position and the “true” position of the respective reference library entry assumed to be unknown, (Δx, Δy ) = (x,calc − xtrue , y,calc − ytrue ), are computed. These differences are then filled into the 2D ‘error histogram’, as indicated in Figs. 3 and 4. Finally, the FWHM of the error histogram (i.e. the average of the FWHM values of its projections on x and y direction) is identified with the spatial resolution achieved in the specific reconstruction scenario.
2.2. The categorical average pattern (CAP) algorithm The second method used to determine the interaction position in the absorber of the Compton camera is the Categorical Average Pattern (CAP) algorithm, a modification of the standard k-NN method (van Dam, 2011). The main difference between the standard k-NN and CAP is that in the latter the k closest matching light distributions are searched amongst the nepp reference light distributions collected independently for each irradiation position. From this subset of k CAP reference distributions an average 2D light distribution is then calculated. The difference DlCAP is then calculated between the average light distribution and the unknown one for each irradiation position according to N
DlCAP =
∑ (Iunknown,i − Iave(l,i))2 , i =1
l = 1, …, (n pos)
(2)
2.4. The reference libraries
where N is again the total number of photosensor pixels, Iunknown,i is the light intensity of the i-th pixel of the unknown event and Iave(l,i) is the intensity of the i-th pixel of the l-th average light distribution. The outcome is then a 2D plot containing n pos intensity differences, each one corresponding to an (x,y) pair of irradiation coordinates. This plot is then smoothed with a 5×5 moving average filter and the coordinates of the minimum value of the plot are assigned as the calculated (x,y) interaction position inside the scintillator crystal. In Fig. 4 the flowchart of the CAP algorithm is shown.
In the process of acquiring the light amplitude reference libraries discussed before, a sequence of five correction steps has to be applied consecutively to the raw data, which will be introduced in the following (Aldawood, 2015): (i) gain matching:Since each PMT segment is read out and processed by an individual electronics chain, gain variations may occur. Injecting two pulser signals of different amplitude into the signal processing chain, gain correction factors can be derived for each signal channel. (ii) QDC pedestal subtraction: Unavoidable electronic dark current integrated and digitized in the QDC modules (‘energy pedestal’) can be quantified in a dedicated measurement without detector input. The resulting pedestal peaks are then taken as lower energy threshold for the data acquisition process. (iii) PMT nonuniformity: Possible gain fluctuations of the multi-anode PMT segments have to be accounted for via the non-uniformity matrix provided
2.3. Determination of the spatial resolution: “Leave-One-Out” method In order to determine the spatial resolution achievable via the standard k-NN or CAP photon-interaction position reconstruction algorithms, the “Leave-One-Out”-method was applied. As schematically indicated in Figs. 3 and 4, this method is based on applying the
Fig. 4. Flowchart of the CAP algorithm (for details see main text).
4
Radiation Physics and Chemistry (xxxx) xxxx–xxxx
S. Aldawood et al.
by the manufacturer. (iv) Spatial homogeneity: Photon scattering and reflections at the side surface wrapping will induce anisotropies in the absorption of isotropic incident radiation. This can be corrected by making use of the internal radioactivity of the LaBr3:Ce scintillator material, safely assumed to be homogeneously distributed in the crystal volume. The internal activity energy spectrum is acquired and a position-dependent correction matrix is derived after gating on the energy region of interest (i.e. the energy of the collimated γ-ray source). (v) Energy gating: An energy gate in the region of the photopeak energy of interest is applied. After preparing the data acquisition as described before, the 2D detector scan can be performed. The intense calibration sources (137Cs/ 82 MBq and 60Co/20 MBq) are placed behind a Densimet collimator (ρ=18.5 g/cm3) with a diameter of 1 mm. In case of 137Cs with a photon energy of 662 keV, the collimator length was 48 mm and the source and its surrounding Pb shielding were placed on a motorized and automatic controlled x-y translation stage in front of the stationary scintillation detector. For 60Co, the higher photon energies (1.17/ 1.33 MeV) required an extended collimator length of 100 mm and a thicker Pb side shielding. Therefore, the setup was modified as indicated in Fig. 5. Now the shielded source/collimator unit was kept stationary, while the detector was moved in front of the collimator opening. In this setup, the collimator was formed by a 100×100×100 mm3 Densimet block with a central bore of 4 mm diameter, allowing to insert an exchangeable collimator tube from sintered WC (ρ=15.6 g/cm3) with an open diameter of 1 mm. In a preparatory measurement series, the position of the LaBr3:Ce crystal inside its aluminum housing was determined by scanning across the detector edges in x and y direction. The rise and fall of the detected photon count rate indicated the crystal position and served to determine the starting and ending points of the subsequent 2D scans. In Fig. 6 (top panel) a subset of the 2D reference library acquired with a collimated (1 mm opening) 137Cs (82 MBq) is shown, representing a scan of 16×16 irradiation positions with a step size of 3 mm in x and y direction. The full scan has been carried out using a fine grid of 0.5 mm step size and 103×103 irradiation positions. For each position, data were acquired during 30 s, allowing to collect at least 400 photopeak events per position. The bottom panel of Fig. 6 displays the corresponding library subset when using a 60Co calibration source (20 MBq) and after gating on the 1.33 MeV photopeak. The correlation of the light amplitude distributions with the source position is clearly visible.
Fig. 6. Subsets of the 2D light amplitude distribution libraries generated with a 137Cs (top) and 60Co source gated at 1.3bsets. 16×16 irradiation positions (step size 3 mm in x and y direction) are shown, measured with a 1 mm collimator opening.
2.5. Systematic studies of the k-NN performance Both position reconstruction algorithms, standard k-NN and CAP, were systematically studied in order to identify the optimum values of the key parameters k (number of best matching reference light distributions) and nepp (number of photopeak reference events per irradiation position). Moreover, the analysis of the achievable spatial resolution was performed for two values of the PMT segmentation, using either the full 256-fold segmentation provided by the electronical readout or combining the signals from each neighbouring four segments to realize a 64-fold segmentation a posteriori by software in the offline analysis. The study was performed for the three available photon energies from 137Cs and 60Co and for two values of the scan step size (0.5 mm, 1 mm) (Miani, 2016). Fig. 7 displays the resulting values of the spatial resolution as extracted by using the standard k-NN algorithm and a reference library with nepp = 400 photopeak events acquired per irradiation position. Clearly visible is the improved spatial resolution for the higher photon energies provided by the 60Co source compared to the curves derived at 662 keV from 137Cs. The minimum of the spatial resolution is consistently reached for values of k around 1000. The best spatial resolution is reached for a PMT granularity of 64 channels and a scan step size of 0.5 mm. Results from a similar study using the CAP algorithm are shown in
Fig. 5. Sketch of the 2D scanning system used to position the LaBr3:Ce scintillator in front of a strong, tightly collimated γ-ray source. The top inset shows photographs of the collimator block with the exchangeable sintered WC collimator tube.
5
Radiation Physics and Chemistry (xxxx) xxxx–xxxx
S. Aldawood et al.
Fig. 9. Determination of the experimental uncertainties of the spatial resolution data, obtained with the CAP algorithm at 1.33 MeV (details: see text).
Fig. 7. Spatial resolution as a function of k (number of nearest neighbours) achieved with the k-NN algorithm, using nepp =400. The 3 parameters photon energy, PMT granularity and scan pitch are distinguished by their line style, marker shape and symbol filling style, respectively.
of the number of best matching light amplitude distributions. The spatial resolution obtained with the CAP algorithm at 1.33 MeV is shown as a function of k (number of nearest neighbours) in Fig. 9. The data refer to the 0.5 mm pitch size libraries, created by using the four subsets with 100 reference events per irradiation position extracted from the full library (with 400 reference 2D light amplitude distributions per irradiation position), for both PMT granularities of 64 and 256, respectively. The index from 1 to 4 indicates the 4 library subsets. Consequently, for a fixed algorithm, PMT granularity value, pitch size and photon energy, the standard deviation of the four spatial resolution values obtained for the four data subsets was calculated and assigned as the uncertainty for this specific scenario. After having determined the optimum value of k for different irradiation energies and readout scenarios, systematic studies of the achievable spatial resolution (i.e. FWHM of the error histogram of the reconstruction algorithm) were performed as a function of the number of photopeak events per irradiation position, nepp . Fig. 10 illustrates the performance of the standard k-NN and CAP algorithms applied to different irradiation energies (662 keV, 1.17 MeV, 1.33 MeV), PMT granularities (64- and 256-fold) and scan pitch sizes (0.5 mm, 1 mm) as a function of nepp , each time applying the optimum value of k as determined before. The top panel shows the corresponding results for 662 keV from 137Cs, while the middle and bottom panels display the results for 1.17 MeV and 1.32 MeV from 60Co, respectively. Consistent trends can be inferred from the three different energy dependent data sets. In general, the spatial resolution improves with an increasing number of reference events per position, due to the larger statistics the unknown photon event can be compared to. Below nepp ≈ 50 , large fluctuations and uncertainties of the individual data sets prevent a clear identification of comparative trends. Above nepp ≈ 50 − 100 (depending on the energy), the CAP algorithm provides (within the experimental uncertainties) at least comparable and mostly even better results than the standard k-NN algorithm. Moreover, the 0.5 mm scan pitch size for a given PMT granularity results in a better spatial resolution compared to the 1 mm case. Important from a practical point of view is the finding that a lower PMT granularity (64 segments) consistently results in better spatial resolution compared to the higher (256 fold) granularity. This behaviour can most likely be attributed to the higher photon statistics per readout pixel in case of the lower granularity, even though the summing of each four pixels was not done a priori electronically in hardware, but a posteriori during the software analysis (van der Laan, 2010). The most important difference of the energy dependent data sets in Fig. 10 is the consistently improving spatial resolution as the photon energy increases.
Fig. 8. Again reference libraries with 400 photopeak events per irradiation position were used. A general feature of all curves is the optimum value of k ≈ 12 (remind that for the CAP algorithm the best matching reference library entries are searched within the ensemble of each irradiation position). While the general trends visible in the CAP data are similar to the ones discussed before for the standard k-NN algorithm, the optimum value of the spatial resolution is consistently lower for CAP compared to k-NN. The experimental uncertainties included in the data shown in Figs. 7 and 8 have been estimated as fluctuations of the spatial resolution between different realizations of the ‘leave-one-out’ procedure applied to the various parameter scenarios. In particular, the initial reference library with 400 photopeak events per irradiation position was subdivided into four subsets with 100 photopeak events each before execution of the ‘leave-one-out’ method for various values k
Fig. 8. Spatial resolution as a function of k (number of nearest neighbours) achieved with the CAP algorithm, using nepp =400. The 3 parameters photon energy, PMT granularity and scan pitch are distinguished by their line style, marker shape and symbol filling style, respectively.
6
Radiation Physics and Chemistry (xxxx) xxxx–xxxx
S. Aldawood et al.
Fig. 11. Energy dependence of the spatial resolution, comparing the two reconstruction algorithms k − NN and CAP in different scenarios. While the scan pitch size was kept constant at 0.5 mm, the PMT granularity was varied between a 64-and 256 fold segmentation of the scintillator PMT readout. The data were acquired with a 1 mm collimator and optimized values of k were used for nepp =400 photopeak events per irradiation position.
tion of the spatial resolution of its absorber component, a monolithic LaBr3:Ce scintillation crystal. The photon entry position at the crystal front surface was determined by using the k-NN reconstruction algorithm in its standard form as well as in the modified version of the Categorical Average Pattern (CAP) algorithm. The applicability of this method beyond its initial energy range of 511 keV annihilation photons (for PET applications) was verified up to 1.33 MeV, with a clear trend of improving spatial resolution as a function of the photon energy, as can be seen in Fig. 11. The plot also indicates the design goal of the Compton camera project, which aims at reaching a position reconstruction resolution in the absorbing scintillator of 3 mm at the multi-MeV photon energies to be expected from the prompt-γ spectrum of hadron-tissue interactions. According to design specification simulations (Lang, 2014), this resolution would enable an overall spatial resolution of the Compton camera of 1.5°–2° in the relevant photon energy range of 2–6 MeV (equivalent to 1.5–2 mm resolution for a distance of 50 mm between the photon source and the first scatter detector, as in the case of a small-animal irradiation scenario). In view of envisaged attempts to further improve the spatial resolution by using an even tighter γ-ray source collimation with a diameter of 0.6 mm, correspondingly reducing the photon flux by about a factor of 3, the saturation of the spatial resolution improvement for nepp ≥ 200 will help to keep the required measurement in a realistically manageable range of a few days. However, the ultimate goal is to commission the Compton camera in the relevant energy range of 4–6 MeV, not accessible via laboratory calibration sources. The energy-dependent trend of the absorber's spatial resolution established in this study nourishes confidence that the design goal can be met, taking into account that with higher photon energies also the contribution of pair creation will increase, giving rise to an even more localized interaction of the photons in the scintillation crystal. Experimental characterization in this energy range might be possible by exploiting the availability of well-collimated, highly monoenergetic photon beams from (nγ) reactions at neutron sources like the GAMS6 facility of the Institut Laue Langevin in Grenoble (Kessler, 2001), as currently being planned. Supported by DFG Cluster of Excellence MAP (Munich-Centre for Advanced Photonics) and King Saud University, Riyadh, Saudi Arabia.
Fig. 10. Spatial resolution as a function of nepp (number of photopeak events per irradiation position for Eγ = 662 keV (top), 1.17 MeV (middle) and 1.33 MeV (bottom), plotted for 3 parameters: type of reconstruction algorithm, PMT granularity and scan pitch size. According to the systematic study discussed before, the optimum value of k was chosen for each parameter value data set.
Apparently, the improvement due the the higher photon count dominates the deterioration of spatial resolution that might be expected from increased intra-crystal scattering (van der Laan, 2010; Maas, 2010). The best values of the spatial resolution determined via the CAP algorithm for nepp = 400 and 0.5 mm scan step size were 4.8(1) mm (662 keV), 3.9(1) mm (1.17 MeV) and 3.7(1) mm (1.33 MeV), respectively. 2.6. Conclusions and outlook The current study presents a key ingredient of the commissioning process of the Garching Compton camera prototype, designed for beam range monitoring in hadron therapy, namely the determina7
Radiation Physics and Chemistry (xxxx) xxxx–xxxx
S. Aldawood et al.
Llosa, G., et al., 2013. Nucl. Instrum. Methods A 718, 130. Lutter, R., et al., 2000. IEEE Trans. Nucl. Sci. 47, 280. Maas, M.C., et al., 2006. IEEE Trans. Nucl. Sci. 53, 1071. Maas, M.C., et al., 2010. Med. Phys. 37, 1904–1913. Marinšek, T., 2015. Thesis, Master., LMU Munich, unpublished. Miani, A., 2016. Thesis, Master, Univ. degli Studi di Milano, unpublished. Min, C.-H., et al., 2006. Appl. Phys. Lett. 89, 183517. Min, C.-H., Lee, H.R., Kim, C.H., Lee, S.B., 2012. Med. Phys. 39, 2100. Perali, I., et al., 2014. Phys. Med. Biol. 59, 5849. Polf, J.C., et al., 2015. Phys. Med. Biol. 60, 7085. Richard, M.-H., et al., 2011. IEEE Trans. Nucl. Sci. 58, 87. Roellinghoff, F., et al., 2014. Phys. Med. Biol. 59, 1327. Santiard, J.C., et al., 1994. CERN-ECP/94 17. Testa, M., et al., 2010. Radiat. Environ. Biophys. 49, 337. Thirolf, P.G., et al., 2016. EPJ Web of Conferences 117, 05005. Thirolf, P.G., et al., 2014. EPJ Web of Conferences 66, 11036. van Dam, H.T., et al., 2011. IEEE Trans. Nucl. Sci. 58, 2139.
References Aldawood, S., et al., 2015. Front. Oncol. 5, 270. Assmann, W., et al., 2015. Med. Phys. 42, 567. COG16article, Particle Therapy Co-operative Group,〈http://www.ptcog.ch〉 Fix, E., Hodges, J.L., 1951. Technical Report 4, USAF School of Aviation Medicine. Gobain, St., scintillation material data sheet, 〈http://www.crystals.saint-gobain.com/ uploadedFilesSG-CrystalsDocumentsBrilLanCe%20380%20Dat%20Sheet.pdf〉 Golnik, C., et al., 2014. Phys. Med. Biol. 59, 5399. Gwosch, K., et al., 2013. Phys. Med. Biol. 58, 3755. Kessler, E.G., et al., 2001. Nucl. Instrum. Methods A 457, 187. Knopf, A.-C., Lomax, A., 2013. Phys. Med. Biol. 58, R131. Kormoll, T., et al. 2012. Nucl. Instrum. Methods A, 11, 4, 2011, Erratum ibid. A 678, (13) 626-627. Krimmer, J., et al., 2015. Nucl. Instrum. Methods A 787, 98. van der Laan, D.J., et al., 2010. Phys. Med. Biol. 55, 1659–1675. Lang, C., et al., 2014. J. Instrum. 9, P01008.
8