Maximum likelihood positioning for gamma-ray imaging detectors with depth of interaction measurement

ARTICLE IN PRESS Nuclear Instruments and Methods in Physics Research A 604 (2009) 359–362

Contents lists available at ScienceDirect

Nuclear Instruments and Methods in Physics Research A journal homepage: www.elsevier.com/locate/nima

Maximum likelihood positioning for gamma-ray imaging detectors with depth of interaction measurement$ Ch.W. Lerche a,, A. Ros b, J.M. Monzo´ a, R.J. Aliaga a, N. Ferrando a, J.D. Martı´nez a, V. Herrero a, R. Esteve a, R. Gadea a, R.J. Colom a, J. Toledo a, F. Mateo a, A. Sebastia´ a, F. Sa´nchez b, J.M. Benlloch b a b

Grupo de Sistemas Digitales, ITACA, Universidad Polite´cnica de Valencia, 46022 Valencia, Spain Grupo de Fı´sica Me´dica Nuclear, IFIC, Universidad de Valencia—Consejo Superior de Investigaciones Cientı´ficas, 46980 Paterna, Spain

a r t i c l e in fo

abstract

Available online 30 January 2009

The center of gravity algorithm leads to strong artifacts for gamma-ray imaging detectors that are based on monolithic scintillation crystals and position sensitive photo-detectors. This is a consequence of using the centroids as position estimates. The fact that charge division circuits can also be used to compute the standard deviation of the scintillation light distribution opens a way out of this drawback. We studied the feasibility of maximum likelihood estimation for computing the true gamma-ray photoconversion position from the centroids and the standard deviation of the light distribution. The method was evaluated on a test detector that consists of the position sensitive photomultiplier tube H8500 and a monolithic LSO crystal (42 mm  42 mm  10 mm). Spatial resolution was measured for the centroids and the maximum likelihood estimates. The results suggest that the maximum likelihood positioning is feasible and partially removes the strong artifacts of the center of gravity algorithm. & 2009 Elsevier B.V. All rights reserved.

Keywords: Depth of interaction PET Monolithic scintillation crystal Maximum likelihood

1. Introduction Most detector designs for dedicated small animal PET scanners are based on arrays of small scintillation crystals. For these designs, the spatial resolution is mainly determined by the size of the crystal segments and an increase of the intrinsic spatial detector resolution requires smaller segments. However, the reduction of crystal dimension entails problems as light collection difficulties, inter crystal scatter, efficiency loss due to packing fractions smaller than 100%, and other. Alternatively, monolithic scintillation crystals together with an array of photomultiplier tubes (PMT) or position sensitive photomultiplier tubes (PSPMT) can be used [1,2]. Digitizing all PMT outputs corresponds to the complete sampling of the light response function (LRF) and therefore allows for highly exact 3D position estimation of the detected g-ray, but requires a very complex and expensive data acquisition system. This explains the widespread use of the center of gravity (COG) algorithm for the impact position estimation. The COG algorithm simplifies significantly the required data acquisition system and can be implemented with passive charge division circuits (CDC) [3]. However, the COG

$ This work was supported by the Spanish Ministry of Science and Innovation under Grant FPA2007-65013-C02-02 and by Generalitat Valenciana local government under Grant (GV/2007/008).  Corresponding author. Tel.: +34 96 387 7007; fax: +34 96 387 7279. E-mail address: lerche@ific.uv.es (C.W. Lerche).

0168-9002/$ - see front matter & 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.nima.2009.01.060

introduces systematic errors at the edges of the monolithic crystal which are especially severe for thick scintillation crystals.

2. Method These systematic errors are actually due to the fact that the COG computes the first order moments

m0 ðr0 Þ ¼

N;M X

f i;j ðr0 Þ

(1)

i;j

mx1 ðr0 Þ ¼

1

N;M X

m0 ðr0 Þ

i;j

xi f i;j ðr0 Þ

(2)

and

my1 ðr0 Þ ¼

N;M 1 X y f ðr Þ. m0 ðr0 Þ i;j j i;j 0

(3)

where f i;j ðr0 Þ is the collected scintillation light for the photodetector element ði; jÞ at position ðxi ; yj Þ, r0 ¼ ðx0 ; y0 ; z0 Þ is the true impact position of the g-ray and ðN; MÞ the numbers of photodetector elements along x and y. All moments m0 ðr0 Þ, mx1 ðr0 Þ, and my1 ðr0 Þ depend on the depth of interaction (DOI) z0 because of f i;j ðr0 Þ. Recently, it was shown that CDC can be enhanced for simultaneously estimating the composed moment of second order

ARTICLE IN PRESS 360

C.W. Lerche et al. / Nuclear Instruments and Methods in Physics Research A 604 (2009) 359–362

m2 ðr0 Þ ¼

N;M 1 X ðx2 þ y2j Þf i;j ðr0 Þ m0 ðr0 Þ i;j i

(4)

of the LRF [4]. Thus the variance varðr0 Þ ¼ m2 ðr0 Þ  ½mx1 ðr0 Þ2  ½my1 ðr0 Þ2 of the LRF can be computed. Since this moment depends on x and y, one has to invert a coupled system of three non-linear equations. Maximum Likelihood (ML) estimation was repeatedly used for similar problems [2,5–7]. We suppose that the probability distribution of observing the measured moments m 2 mx1 ; my1 ; m2 are independent Gaussian distributions with means mx1 , my1 , var, and constant errors sx  1:5 mm, sy  1:5 mm and svar  3 mm. These errors correspond to the measured resolutions for the centroids and the standard deviation at the center of the PMT’s sensitive area [8]. The log likelihood function is then given by the sum of mean squares ln L ¼ 

 x   y 2   m1  my1 m1  mx1 2 m2  var 2 pffiffiffi ffiffiffi  p  pffiffiffi x y var 2s 2s 2s

(5)

which takes its maximum at the ML estimate r^ 0 for the g-ray impact position r^ 0 ¼ argmax ln L.

(6)

r0

In monolithic crystals, the scintillation light transport obeys the inverse square law f ðr; r0 Þ ¼

f 0 deff 4pjr  r0 j3

(7)

where f ðr; r0 Þ is the LRF of the detector and deff is the effective distance between the photo-conversion position r0 and the optical interface between crystal and PSPMT-window. deff accounts for Snell’s Law (refer to Fig. 1) and can be found by solving the quartic equation !

2

2

0 ¼ b  deff  a2 þ d a¼

n2 b n1

and

2

with qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi deff ðdeff þ tÞ2 þ r2 . b¼ deff þ t

(8)

Finally, the LRF has to be integrated over the intervals Oi;j that correspond to the pixels ði; jÞ of the PSPMT ZZ f i;j ðr0 Þ ¼ f ðr; r0 Þ dx dy. (9) Oi;j

We computed the moments mx1 , my1 , and var for a detector head consisting of a 42  42  10 mm3 monolithic LSO crystal and the PSPMT H8500 (Hamamatsu Photonics K.K., Japan) with 8  8 detector segments and 49  49 mm2 sensitive area. Since the numerical integration takes a long time, we computed the moments and gradients in advance for a set of impact positions

r0 lying on a three dimensional regular grid of interpolation nodes px;y;z (43  43  11 points). For the calibration of the position reconstruction we measured the mean centroids and mean variance at each of the 81 test position and assumed that these mean values represent approximately the centroids and variances at 5 mm DOI and the corresponding positions. We then computed the differences between these values and the model prediction for the corresponding position at an interaction depth of 5 mm. Bilinear interpolation was used to add these differences to the interpolation nodes px;y;z . The resulting values were stored in look-uptables which allows for fast trilinear interpolation of the expected detector response for any point 21 mmpx0 ; y0 p21 mm and 0 mmpzp10 mm during the ML-estimation. The maximization of the log-likelihood function was done using the Polak–Ribiere method [9]. This conjugate gradient algorithm does not require knowledge of the second derivatives but constructs the sequence of search direction only by line minimizations and gradient evaluations. The algorithm was implemented using NVIDIA’s Compute Unified Device Architecture (CUDA) [10]. The CUDA C language extension provides access to the NVIDIA Graphics Processing Unit (GPU) for general purpose scientific computing and takes advantage of the massively multithreaded single instruction multiple data (SIMD) architecture of these devices [11]. Two identical g-ray detectors were used to evaluate the performance of the proposed DOI estimation method. Each module consists of a single monolithic LSO scintillator block with spatial dimensions 42  42  10 mm3 . The five crystal surfaces that were not coupled to the PSPMT were fine ground and covered with AY103 black epoxy resin. The LSO block was coupled by optical grease (rhodorsil Paˆte 7, Silitech AG, Bern, Switzerland) to the entrance window of H8500. Energy, both centroids, and the variance of the signal distributions of coincidence events were measured at 81 different positions (see Fig. 2). The test detector was placed with translational stages such that the collimated beam impinges at these positions. 1 92 000 coincidence events were registered before moving to the next position. All events whose energy lies outside of the full width at tenth maximum interval are discarded. Collimation of the g-ray beam was achieved by placing the 22 Na source as close as possible to the test detector and by accepting only events inside a small circular spot (diameter 1.4 cm) at the center of the coincidence detector’s sensitive area. The normalized moments for all events were

19

1

2

3

. . .

9

14.25 r0

Scintillator

9.5 4.75 0

a

d b deff n1 n2

PMT-Window

−4.75 −9.5 −14.25

t

robservation r Fig. 1. Refraction of the scintillation light path (thick, green line). For interpretation of the references, to color in this figure legend the reader is referred to web version of this article.

−19 y [mm]

73

81

−9.5 19 9.5 0 x −19 [mm] −14.25 −4.75 4.75 14.25

Fig. 2. Location of the 81 test positions. The sequence of measurements starts at point 1 ð19; 19Þ mm and follows the arrows.

ARTICLE IN PRESS C.W. Lerche et al. / Nuclear Instruments and Methods in Physics Research A 604 (2009) 359–362

computed and the true impact position was estimated by maximizing ln L.

3. Results Figs. 3(a) and (b) show density plots of the centroids and the ML-estimates, respectively. In the plot that shows the centroids, only the 5  5 centric beam positions can be identified. A significant improvement is observed for the density plot of the ML-estimates where almost all 9  9 beam positions can be identified. Figs. 3(c)–(f) show the spatial resolutions for centroids and ML-estimates for all 81 beam positions. These resolutions were measured by fitting the 1D-histograms of both centroids and both ML-estimates to Gaussian curves and determining their full width half maximum values. The measured resolutions were corrected for the finite source diameter of  0:92 mm [8] by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dxcorr ¼ dx2meas  0:922 (10) qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dy2meas  0:922 .

(11)

4. Discussion and conclusions The presented method is sufficiently fast to allow real-time position reconstruction for the Albira small animal PET up to activities of 200 mCi [13]. The Albira PET scanner uses eight g-ray detectors as described above and has an octagonal scanner geometry. Its sensitivity was reported to be 4 kcps/MBq. Fig. 3(b) shows that the quality of the position reconstruction is different

20

20

10

10

0

-10

-20

-20 -10

0 x [mm]

10

20

-20

20

15

15

δy [mm]

20

10

-10

0 x [mm]

10

20

10 5

5 0

20

40 position#

60

80

20

20

15

15

δy [mm]

δx [mm]

0

-10

-20

δx [mm]

This was repeated for all 81 beam positions. Table 1 gives the mean, maximum, minimum, and standard deviation values for the resolutions that were measured at all 81 beam positions. Note that also the best achieved spatial resolution (minimum value) was better in the case of ML positioning. Execution time was measured for two implementations (see Table 2). The multithreaded implementation on the NVIDIA GPU 80 outperformed the single-threaded implementation on the AMD Athlon 64  2 3800þ by a factor of 12. Execution times reported for the cMiCE PET are very similar to our results for the GPU [12]. Fig. 4 shows the distribution of necessary iterations for all reconstructed events.

y [mm]

y [mm]

dycorr ¼

361

10

0

20

40 position#

60

80

0

20

40 position#

60

80

10 5

5 0

20

40 position#

60

80

Fig. 3. Density plots of centroids and ML-estimates and spatial resolutions for the 81 beam positions. Error bars represent standard deviations. (a) density plot of measured centroids; (b) density plot of ML-estimates; (c) measured x centroid resolutions; (d) measured x ML-estimate resolutions; (e) measured y centroid resolutions; (f) measured y ML-estimate resolutions.

ARTICLE IN PRESS 362

C.W. Lerche et al. / Nuclear Instruments and Methods in Physics Research A 604 (2009) 359–362

Table 1 Mean, std. dev., max. and min. values for resolutions of the centroids and reconstructed positions. Mean

Std. dev.

Min.

Max.

dx (ML-estimate) (mm) dy (ML-estimate) (mm) dx (Centroid) (mm) dy (Centroid) (mm)

2.1 1.9 3.4 3.3

0.9 0.9 3.2 3.1

1.1 1.1 1.4 1.3

6 4.8 20.9 19.9

12 Fraction of Events [%]

Measure

14

10 8 6 4 2

Table 2 Execution times of the algorithm for different processing units.

0 0

Processor

Processed single events per second

AMD Athlon 64  2 3800þ (single-threaded routine) NVIDIA GeForce 8800 GT (multi-threaded routine)

20  103

5

10

15 Iterations/Event

20

25

Fig. 4. Distribution of necessary iterations. The mean value is 13 iterations.

240  103

for the four corners and sides. This is probably due to an offset of the crystal, i.e., the center of the crystal does not exactly align with the center of the PSPMT’s window. This leads to offsets in the integration limits of the response model Eq. (9) and, therefore, to differences in the reconstruction quality. Nevertheless, a significant improvement of the spatial resolution near the crystal borders was achieved. Almost all 9  9 beam positions can be identified. The measured resolutions for the reconstructed positions are very similar to the values reported in [8], where impact position estimation was done by polynomial interpolation. The positioning scheme for the cMiCE PET performs slightly better at the crystal borders [12]. However, the LSO crystals used for the cMiCE scanner are at least 20% thinner and 19% larger than the crystals for the present study. For lower values of the ratio between crystal thickness and side length, the COG artifact will be smaller. Further improvement of the presented reconstruction

scheme can be achieved by a more sophisticated calibration of the LRF-model and a better LRF-model itself, e.g., by taking into account possible residual reflections at the black epoxy surfaces. References [1] [2] [3] [4] [5] [6] [7] [8]

[9] [10] [11] [12] [13]

G. Muehllehner, et al., IEEE Trans. Nucl. Sci. NS-35 (1988) 670. J. Joung, et al., Nucl. Instr. and Meth. A 489 (2002) 584. S. Siegel, et al., IEEE Trans. Nucl. Sci. NS-43 (1996) 1634. C.W. Lerche, et al., IEEE Trans. Nucl. Sci. NS-52 (2005) 560. D. Gagnon, et al., IEEE Trans. Med. Image 12 (1993) 101. L. Furenlid, et al., in: Real Time Conference, 14th IEEE-NPSS, 2005, 4pp. S. Moehrs, et al., Phys. Med. Biol. 51 (2006) 1113. C.W. Lerche, Depth of interaction enhanced gamma-ray imaging for medical applications, Ph.D. Thesis, University of Valencia, 2006, arXiv:physics/ 0611011. W.H. Press, et al., Numerical Recipes in C, second ed., Cambridge University Press, Cambridge, 1992. NVIDIA Corp., Nvidia cuda programming guide 2.0, hhttp://www.nvidia.comi, 2008. J.D. Owens, et al., Comput. Gr. Forum 26 (2007) 80. R. Miyaoka, et al., IEEE Trans. Nucl. Sci. NS-54 (2007) 1561. J. Benlloch, et al., Rev. de Fı´s. Me´d. 8 (2007) 315.