Charged-particle emission tomography.

PURPOSE: Conventional charged-particle imaging techniques - such as autoradiography - provide only two-dimensional (2D) black ex vivo images of thin tissue slices. In order to get volumetric information, images of multiple thin slices are stacked. This process is time consuming and prone to distortions, as registration of 2D images is required. We propose a direct three-dimensional (3D) autoradiography technique, which we call charged-particle emission tomography (CPET). This 3D imaging technique enables imaging of thick tissue sections, thus increasing laboratory throughput and eliminating distortions due to registration. CPET also has the potential to enable in vivo charged-particle imaging with a window chamber or an endoscope.
METHODS: Our approach to charged-particle emission tomography uses particle-processing detectors (PPDs) to estimate attributes of each detected particle. The attributes we estimate include location, direction of propagation, and/or the energy deposited in the detector. Estimated attributes are then fed into a reconstruction algorithm to reconstruct the 3D distribution of charged-particle-emitting radionuclides. Several setups to realize PPDs are designed. Reconstruction algorithms for CPET are developed.
RESULTS: Reconstruction results from simulated data showed that a PPD enables CPET if the PPD measures more attributes than just the position from each detected particle. Experiments showed that a two-foil charged-particle detector is able to measure the position and direction of incident alpha particles.
CONCLUSIONS: We proposed a new volumetric imaging technique for charged-particle-emitting radionuclides, which we have called charged-particle emission tomography (CPET). We also proposed a new class of charged-particle detectors, which we have called particle-processing detectors (PPDs). When a PPD is used to measure the direction and/or energy attributes along with the position attributes, CPET is feasible.

1. INTRODUCTION

Charged particles emitted by radioisotopes — including alpha particles, beta particles, conversion electrons, and Auger electrons —are widely used in biology, pharmacology, and radionuclide therapy. Radioactive tracers, usually beta-particle emitting, are commonly used to understand biochemical processes in biology[1-5] and to provide insight into the localization and metabolism of pharmaceuticals in drug discovery and development.[6-8] Charged particles also have desirable properties for radionuclide therapy in cancer treatment,[9-12] including localized deposition of energy, which will lead to little damage to the nearby tissues, provided that the charged particles are guided specifically to their target. Antibodies to tumor antigens as well as peptides targeting tumor-cell-surface receptors can provide this target-specific guidance. The cellular-level distribution of the charged-particle-emitting radionuclide is often valuable information, because it is useful in an accurate estimation of drug distribution or radiation dose.

When charged-particle emissions are accompanied by gamma or x-ray photons, standard imaging methods, such as positron emission tomography (PET) and single-photon emission tomography (SPECT) can be used, but the usual inherent limitations on spatial resolution will apply. The resolution of preclinical SPECT systems is about 0.5–2 mm.[13] For PET systems, the resolution for preclinical PET is about 1–2 mm.[13] Neither PET nor SPECT can provide cellular information, which requires resolution of around 10 μm. One of the fundamental factors that limit the resolution of PET is the positron annihilation range, which can be eliminated by direct positron imaging.[14] When no photons are present, direct charged-particle imaging is required.

The current technique for charged-particle imaging, autoradiography, is two-dimensional, ex vivo imaging of thin tissue slices.[8,15,16] In autoradiography, the radiopharmaceutical is introduced into the living subject, and after a suitable time for it to equilibrate, the organ or tissue of interest is removed and cut into thin slices (4- to 10-μm thick when cellular-level radionuclide distribution is required) with a cryomicrotome. One can get the three-dimensional distribution of the radioisotope by imaging multiple thin slices and coregistering them,[17] but in practice, this procedure is very difficult because of distortion of the thin slices during the transfer from the microtome to the imaging detector.

In this paper, we present a direct three-dimensional autoradiography technique, which we call charged-particle emission tomography (CPET).[18,19] With CPET, the volumetric distribution of the radionuclide can be imaged without requiring registration of 2D slices. CPET would allow imaging of thick tissue sections with high spatial resolution (perhaps rivaling that of thin-slice autoradiography), which reduces the number of tissue samples needed for imaging and hence increase laboratory throughput. Furthermore, in vivo CPET is possible with a small-animal window chamber,[20,21] clinically for superficial lesions, and potentially with an endoscope or a surgical probe.

As discussed above, a crucial component of charged-particle emission tomography is the charged-particle detector. Current charged-particle detectors, such as films,[22] phosphor imaging systems,[23] and scintillation gas detectors,[24-26] either are particle-counting or estimate only the 2-D positions of the detected particles.

The key innovation of CPET is the use of charged-particle detectors that provide information about not only the location of the particle when it interacts with the detector but also its direction and/or energy. This new class of charged-particle detectors, which we call particle-processing detectors (PPDs), detects single particles, and for each detected particle, they measure a subset of particle attributes, such as position, direction, and energy. The output of the detector is a list of attributes, and each entry of the list contains the estimated attributes of a detected particle.

This paper is organized as follows. Section 2 discusses the basic configuration, forward model, and reconstruction algorithms of CPET. Sections 1 and 4 discuss two special cases of CPET and demonstrate the feasibility of both cases with reconstruction results from simulated data. Section 5 introduces the concept of PPDs and demonstrates one PPD, which is a two-foil detector, experimentally. Section 6 summarizes our results and presents future research directions. Two supplementary materials can be downloaded online: Supplement A presents several designs of PPDs. Supplement B provides statistical analysis of a two-foil detector.

2. GENERAL DISCUSSIONS

2.A. Concepts

Two special cases of CPET are alpha emission tomography (αET, pronounced as AlphET) and beta emission tomography (BET). αET is an emission tomography technique that relies on alpha-particle-emitting radioisotopes to form tomographic images. BET is a technique to image the volumetric distribution of tracers that emit fast electrons.

In a CPET system, a PPD is placed in contact with or in close proximity to one side of a sample tissue, as shown in Fig. 1, rather than surrounding the tissue with detectors as in SPECT or PET. In CPET, neither the detector nor the tissue is moving; therefore, CPET is an extreme limited-angle tomography system. By employing PPDs, a richer dataset can be collected and used to perform 3D reconstruction of the distribution of the radionuclide that emitted the particles.

Fig. 1

Illustration of CPET systems. (a) BET system and Geant4-simulated[27,28] tracks of 400-keV electrons in a 400-μm-thick tissue. (b) αET- and Geant4-simulated tracks of alpha particles emitted from 239Pu source in a 40-μm-thick tissue. [Color figure can be viewed at wileyonlinelibrary.com]

The tissue thickness that can be imaged with CPET is limited by the type and the energy of the particles emitted from the source. For example, the positrons emitted by 18F sources have maximum energy of 633 keV, which has a continuous-slowing-down-approximation (CSDA) range[29] of 2.43 mm in soft tissue. The energy of an alpha particle is usually in the 4–8 MeV range, which corresponds to a penetration range of 26–78 μm in soft tissue. A source located deeper than the particle range in tissue cannot be imaged by CPET.

2.B. Forward model

In a CPET system, the object is the spatial density of radioactive decays per unit time, which is a three-dimensional distribution function. We denote the object as a scalar-valued function f(R), where R = (x,y,z) is a point in the 3D Euclidean space. The units of f(R) are 1/(s·mm3). The object function f(R) is also referred to as f in the following discussion.

A particle detected by the detector can be characterized by several attributes, including the 2D position r, the energy E and the direction s⊥ of the particle when it arrives at the detector plane. In our notation, s⊥ = (s,s) is a 2D vector (not a unit vector) where s and s are the direction cosines. An ideal PPD measures a subset of these attributes for each detected particle.

For the jth detection event (j = 1,. . .,J), we can organize the attributes into an ‘attribute vector’ Â, where the circumflex represents estimation. The output of a PPD is a list of J attribute vectors, where J is the total number of detected events. We will denote this list as {Â1, Â2, . . ., Â} and we will use it to introduce the following point process[30]

The mean of u(Â) over an ensemble of attribute lists (more specifically, over J and for each J over all Â) is proportional to the probability density of a particle detected with estimated attributes Â.[30]

For a given object f and a fixed exposure time τ, the mean of the list-mode data is[31]

where pr(A|R) is the probability density of a detection event with true attribute vector A given that the event originated with an emission at R; pr(Â|A) is the probability density of estimate  given that the true attribute vector is A; and S(R) is the probability that an emission at R will be detected and result in a list entry. The function S(R) is often referred to as the sensitivity function. There are three possible ways that a particle does not produce a list entry after its emission: (a) The particle is stopped in the tissue after losing all of its energy through interactions; (b) the particle escapes the tissue but it does not interact with the detector; or (c) the particle detection event is rejected based on a likelihood window[32,33] on pr(Â|A).

The integral in Eq. (2) indicates a linear system which can be expressed in a more abstract notation as

where the operator ℒ is a linear operator defined by Eq. (2). The kernel function of ℒ is

The function prf(Â|R) is the count density in attribute space provided that a particle is emitted at point R. The notation prf is used, because the function describes the system response to a point source and such a function is often referred to as a point response function.

The function pr(A|R)S(R), which describes the interactions of particles in tissue, can be approximated by Monte Carlo simulations[30,34] but analytical forms are often valuable. The analytical forms can be either derived directly or solved with Boltzmann transport equation.[30] In this paper, the function pr(A|R)S(R) is derived analytically for alpha particles and approximated with Monte Carlo simulations for beta particles.

If maximum-likelihood (ML) estimation is used to estimate the detector outputs, the function pr(Â|A) can be expressed asymptotically as a multivariate Gaussian probability distribution function with an inverse covariance matrix equal to Fisher information matrix.[35] The asymptotic properties of ML estimation is satisfied in the limit of a large number of secondaries (photoelectrons for scintillation detectors or electron-hole pairs for semiconductor detectors).

Another way to obtain the functional kernel of the operator ℒ is through experiments. The functional kernel prf(Â|R) is the mean detector response to a point source located at position R. One can scan through a grid in the object space with a point source and measure the detector response.

A CPET system has translational symmetry and rotational symmetry if the following conditions are satisfied: (a) the detector is large compared to the size of the object, and (b) the medium where radioactive sources are located is uniform. The symmetries can be used to simplify the characterization of the imaging system.

2.C. Reconstruction algorithms

The end goal of a CPET system is to solve the inverse problem to reconstruct the 3D distribution of the radioactive tracers, f(R), given a list of detected events {Â1, Â2, . . ., Â}. The inverse problem can be solved implicitly with iterative algorithms by minimizing some scalar-valued functional Q(f, u) that depends on the object f and the list-mode data u. Two iterative algorithms, the Landweber algorithm and the list-mode maximum likelihood expectation maximization (LMMLEM) algorithm, are introduced and used in BET and αET, respectively. On one hand, we chose Landweber algorithm for BET, because there is no analytical models to describe the complex beta-particle interactions with matter. On the other hand, MLEM is implemented for αET, because alpha particles approximately travel in straight lines which is easy to describe analytically and provides the necessary probability function required to implement the MLEM algorithm.[36]

The feasibility of CPET is demonstrated with reconstruction results from simulated data. Using Geant4-simulated data,[27,28] we reconstruct charged-particle-emitting objects for BET and αET in the following two sections, respectively.

3. BET: BETA EMISSION TOMOGRAPHY

3.A. Introduction

In BET, molecules or cells of interest are labeled so that they emit fast electrons (i.e., beta particles, conversion electrons, or Auger electrons); by imaging the fast electrons, the distribution of the molecules or cells of interest can be reconstructed.

Beta particles have broad energy spectra, while conversion electrons and Auger electrons have discrete energy spectra. Due to its small mass, a fast electron tends to have a tortuous path while traversing matter, as shown in Fig. 1(a). The mean energy loss of a fast electron in a material is a function of many variables, including the path length of the particle, the initial energy of the particle and the material through which the particle propagates. Despite the tortuous tracks of beta particles in tissue, the directions and the residual energies of the detected beta particles both carry information about where the particles are emitted. The particle attributes, including position, direction and/or energy, enable 3D reconstruction of BET.

The detector employed by BET can be any type of PPD that is sensitive to fast electrons. The attributes measured can be (a) position and direction or (b) position, direction, and energy of each detected particle.

In the following discussion, for demonstration purposes, we are going to focus on a 18F source and one specific detector that measures direction and position of each detected particle. The isotope 18F, which emits positrons isotropically with broad energy spectra, is widely used in PET imaging. Note, in PET, the annihilation photons are detected; while in BET, the positrons are detected before annihilation.

3.B. System configuration

We introduce a two-foil detector[37] that is sensitive to beta particles and measures the position and direction of each detected particle. A proposed setup is illustrated in Fig. 2. The two-foil detector consists of two layers of ultrathin phosphor foil, an image intensifier, a high-numerical-aperture lens system, and a CMOS sensor.

Fig. 2

Illustration of a BET system. A beta particle emitted within the tissue propagates towards the detector (dashed box). A two-foil detector placed on the surface of the tissue measures position (x̂, ŷ) and direction (ŝ, ŝ) of the beta particles leaving the tissue.

When a charged particle enters the detector and interacts with both phosphor foils, light is generated at both interaction locations (one interaction location for one foil). The light produced by the phosphors is amplified by the image intensifier, collected with the high NA lens system, and imaged by the CMOS sensor. By processing signals collected on the CMOS sensor, one can estimate the two interaction locations, and hence the position and direction of the incident particle. More details about a two-foil detector are provided in Section 5.

A two-foil detector with the following properties is considered: At each phosphor foil, 500 photons are emitted upon interaction of a beta particle with the phosphor foil; the quantum efficiency of the image intensifier is 80%; z1 = 220 μm and z2 = 20 μm (z1 and z2 are illustrated in Fig. 2). In such a system, the number of photoelectrons produced is N = 200, and the distance between the two foils is d = z1–z2 = 200 μm. According to a Fisher information analysis presented in Supplement B, the uncertainties of the detector are σ(tan θ) = σ(tan θ) = 0.0676 and σ(x) = σ(y) = 13.5 μm.

We focus on the potential advantages that particle-processing detectors bring in to charged-particle imaging. Therefore, we studied the system performance based on a nearly perfect detector. The assumption that 200 photoelectrons will be produced is optimistic with current technology. As presented in Ref. 37, when a 1-MeV electron penetrates a layer of 3.5-μm-thick P43 phosphor, the energy deposited in the phosphor follows a Landau distribution[38] with peak located at 1.5 keV. The most probable energy deposition, 1.5 keV, corresponds to about 53 photons produced in P43 phosphor.[37] We assumed that the number of optical photons generated in each phosphor is 500, which can be achieved by using a thicker phosphor or phosphor with higher light output efficiency. The quantum efficiency of the image intensifier, which is assumed to be 80%, is higher than the maximum quantum efficiency of a currently available image intensifier in the market, which is about 50%.[39] Furthermore, the scattering of beta particles within phosphor layer 1[37] is not considered. However, we believe that a particle-processing detector (a two-foil detector or other detectors discussed in Supplement A) that measures the position, direction, and/or energy of an incident beta particle to relatively high resolution will be feasible as the technology improves.

3.C. Forward model

A Monte Carlo simulation toolkit —Geant4[27,28] —is used to simulate the propagation of beta particles in tissue. From the simulated data, the distribution pr(A|R)S(R) can be obtained. Point sources are simulated at different locations R in the 3D object space. As the system has translational symmetry on xy-plane, enough information about the system can be acquired by simulating a point source, whose location is a function of just depth. In other words, the system operator can be built by scanning a point source through one line parallel to z-axis in the 3D object space instead of scanning through a 3D grid, which significantly reduces the amount of simulation work. The simulation results in a set of attributes lists, each of which corresponds to one depth.

3.D. Reconstruction algorithm

The forward model of the system is given by Eq. (3). The inverse problem to find the object f given list-mode data u(Â) can be solved by a Landweber iterative algorithm. With positivity constraints, the algorithm takes the form[30]

where k is the iteration index, C is a constant, ℒ† is the adjoint operator of the system operator ℒ and 𝓟+ is an operator that sets negative values to zero and makes no change to the positive values. The constant C can be tuned to ensure the convergence of the algorithm. If we denote μ as the maximum eigenvalue of ℒ†ℒ operator, then the algorithm converges to a solution that satisfies ℒ†u – ℒ†ℒf̂ = 0 when |C|2μ<2.[30] If the initial estimate contains no null functions, the algorithm converges to a solution with no null functions, because the operator ℒ† eliminates null components.

The functional form of ℒ†u is

where prf(Â|R) is considered as a function of R and the form of prf(Â|R) is given by Eq. (4).

The operator ℒ†ℒ transforms a function in the 3D-object space to another function in the 3D object space: where k(R,R′) is the functional kernel of ℒ†ℒ given by

The reconstruction algorithm is developed on graphics processing units (GPU) with the CUDA programming language.

3.E. Reconstruction results with simulated data

To acquire a set of data as an input to the reconstruction algorithm, radioactive objects were simulated in a slab of tissue. The radioactivity was confined in a volume of 0.6 ×0.6 ×0.1 mm3 centered at (0,0, −0.05) mm. Three objects were placed at different depths. Beta particles were emitted isotropically in all directions. A planar detector that measures position and direction of each particle was assumed to be placed on the surface of the tissue.

The objects included two bars centered at z = 85 μm and z = 15 μm respectively and a ‘CGRI’ pattern centered at z = 50 μm. The thickness of the bars was 10 μm and the thickness of the ‘CGRI’ pattern was 20 μm. Within the boundaries of each structure, 18F sources are distributed uniformly. The xy-projection view and yz-projection view of the 3D objects are shown in the upper left corner of Fig. 3(a).

Fig. 3

Two projection views of the 3D reconstructed volume (a) versus the 2D autoradiograph (b). Two projection views of the true emission pattern are shown in the upper left corner.

With 106 simulated events, 4.0×105 beta particles were collected by the detector. The output of the simulation was a list of particle attributes. The particle–tissue interaction was simulated in Geant4. To simulate the detector estimation uncertainty, we added Gaussian random variables to the simulated attributes. The uncertainty of the estimates was assumed to be σ(tan θ) = σ(tan θ) = 0.0676 for direction and σ(x) = σ(y) = 13.5 μm for position as discussed in Section 3.B.

The list-mode calibration and detector data were binned into histograms. In the data space, the size of the bins in each dimension were defined as: Δx = Δy = 4 μm, Δθ = 2° and Δφ = 2°.

For the reconstruction, we considered a discretized field of view with 256 ×256 ×20 voxels. Each voxel measures 4 μm ×4 μm ×5 μm. With very good approximation, the system approximately exhibits lateral shift invariance which simplifies the calculations if we represent the object by its discrete Fourier transform with respect to variables (x,y). The iterative steps were calculated in Fourier space. Every 10 iterations, the estimated object was transformed from Fourier space back to object space and positivity constraints were enforced in object space. Reconstruction results were obtained with 1000 iterations. The top projected view and the side projected view of the reconstructed 3D volume are shown in Fig. 3(a).

For comparison, a conventional beta autoradiograph was calculated from the same detector data. In conventional autoradiography, the direction of a particle is not measured and the output image is an integration of the radioactivity over time at the detector surface. Therefore, an autoradiograph is calculated by ignoring the direction information and binning the detected particles into 4 μm ×4 μm pixels. The autoradiograph is shown in Fig. 3(b).

3D BET reconstruction from simulated data was achieved. In comparison to the autoradiograph, beta emission tomography shows improved resolution. Furthermore, the yz-projection of the reconstructed volume shows that BET reveals some depth information.

4. αET : ALPHA EMISSION TOMOGRAPHY

4.A. Introduction

The distinctive physical properties of alpha particles enable αET. Unlike beta particles, alpha particles have discrete energy spectra, with highly monoenergetic emission lines associated with particular nuclear transitions. In low-atomic-number materials such as water or tissue, alpha particles interact with matter primarily through Coulomb forces between their positive charges and the negative charges of the orbital electrons within the absorber atoms. At any given time, the particle is interacting with many electrons, so the net effect is to decrease its velocity continuously until the particle is stopped.

Except at their very end, the tracks tend to be quite straight because the particle is not significantly deflected by one encounter, and interactions are statistically uniform in all directions. The distance an alpha particle traveled is therefore characterized by the energy deposited in a given absorber material.[40] Hence, the path length is a function of the residual energy of the particle. When an alpha particle is detected on a plane at location (x,y) with energy E, the source would be restricted to a spherical shell centered at (x,y) with radius determined by E, if homogeneous medium is assumed. The subscript d denotes detector.

4.B. System configuration

The detector used for αET can be any PPD that is sensitive to alpha particles and the attributes measured can be (a) position and energy, (b) position and direction, or (3) position, direction, and energy. In this paper, we consider one specific setup as an example.

The αET imaging system considered here includes a hybrid semiconductor pixel detector which measures the position and energy of each detected particle. Our requirement for the detector is that it provides accurate position information as well as good energy resolution. Semiconductor detectors allow for good energy resolution because the average energy necessary to create an electron-hole pair is smaller than that needed for other types of charged particle detectors. The setup is illustrated in Fig. 4. Alpha particles emitted by a radioactive isotope inside the tissue are detected by the silicon sensor and produce measurable signals in the detector application-specific integrated circuits (ASICs).

Fig. 4

Illustration of an αET system. Alpha particles emitted within the tissue propagate along straight lines. A silicon pixel detector (dashed box) placed on the surface of the tissue measures position (x̂, ŷ) and energy Ê of the alpha particles leaving the tissue.

The hybrid semiconductor pixel detector (ModuPIX; WidePIX company) considered has 256 ×256 pixels, each 55 μm ×55 μm. The detecting area of the silicon sensor is about 14 mm ×14 mm. The large number of electron-hole pairs produced upon interaction of an alpha particle with the detector material gives rise to pixel intensities that approximate a 2D Gaussian function of spatial variables. By fitting the signals to a Gaussian function, a subpixel spatial resolution of about 750 nm for equivalent 10 MeV alpha particles was achieved.[41] With bias voltage at 100 V, the energy resolution was reported at about 50 keV FWHM for 5.5 MeV alpha particles.[42]

4.C. Forward model

Alpha particles travel through tissue in a nearly straight path with energy decreasing continuously. The energy E of an alpha particle as a function of the distance ℓ the particle traveled has been described by National Institute of Standards and Technology (NIST).[29] The count density on true underlying attribute vectors when a particle is emitted from position R is

where R = (x,y,z), which is the emission position; and , which represents the distance from the emission point of an alpha particle to the position where it enters the detector. We set the plane z = 0 at the detector–tissue boundary.

The probability of a particle emitted at position R being detected, which is also referred to as the sensitivity function, is

where ℓ0 is the range of an alpha particle in a given material; rect(t) is a rectangular function defined by

The particle range ℓ0 depends on the energy of the particle and the material the particle travels through. For a 5.15 MeV alpha particle (239Pu), ℓ0 is about 39.4 μm in water (approximate soft tissue).

If there is an accurate model[43] to describe the detector signals and maximum-likelihood estimation[30,44] is used to estimate the particle attributes, then the Fisher information matrix[35] can be used to derive the measurement uncertainty pr(Â|A). Here we make an assumption that the three estimated attributes (x̂, ŷ, Ê) are not correlated. With this assumption, the covariance matrix is diagonal. Denote (σ, σ, σ) as the diagonal elements of the covariance matrix,

where = (σ,σ,σ).

Combining Eqs. (4), (9), and (12), we can calculate prf(Â|R), the functional kernel of ℒ.

4.D. LMMLEM algorithm

The maximum-likelihood expectation maximization (MLEM) algorithm is an iterative algorithm to solve an inverse problem. The name of MLEM algorithm comes from its property that it maximizes the likelihood (ML) for a Poisson data model and its derivation by alternating expectation (E) and maximization (M) steps. When list-mode data are used as an input, list-mode MLEM or LMMLEM algorithm can be used for the reconstruction.

Consider a discrete object to be estimated from list-mode measurement data. Denote the discrete object as f = {f1, f2, …, f}, where f is the expected number of alpha particles emitted from voxel n (n = 1, 2, …, N and N is the total number of voxels) per unit time given by

where V is the volume of the n voxel in object space. The units of f are 1/s. The probability density of measuring attributes  when the detected particle is emitted from voxel n is

where S is the voxel sensitivity defined as the probability of a particle emitted from voxel n being detected.

The list-mode EM algorithm takes the form:[45-48]

where k is the iteration index, τ is the exposure time, {Â1, Â2, … Â} is the list-mode data, and J is the number of particles detected.

MLEM is a multiplicative algorithm and preserves positivity. If the initial estimate f(0) is non-negative, estimates from all subsequent iterations remain non-negative. Furthermore, the list-mode data from αET are Poisson point processes, as radioactive decays are independent; therefore, MLEM converges to an estimate that maximizes the log-likelihood for a given set of αET data.

4.E. Reconstruction result with simulated data

Simulated data are generated as an input to the reconstruction algorithm. The simulated objects consist of two layers: Chinese characters ‘Jun’ and ‘Ding’ located at z = 6 μm and z = 16 μm, respectively, each 4 μm thick. Within the boundaries of each character, 239Pu sources, which emit 5.15 MeV alpha particles, are distributed uniformly. The number of particles simulated is 106. Two slices of the true decay patterns are shown in Figs. 5(a) and 5(b). As a side note, the object is continuous, but the decay patterns are displayed on the same voxel grid as the reconstruction results.

Fig. 5

Two slices of the true decay pattern (a), (b), two slices of the object reconstructed (c), (d), and a conventional autoradiography (e) from 106 events. The top images (a, c) show one slice located at z = 5.5 μm, and the bottom images (b), (d) show one slice located at z = 16.5 μm. Each slice is 1 μm thick. (Adapted from Ref.[49].)

The alpha-particle propagation in a slab of tissue is simulated in Geant4. The output of the Geant4 simulation is a list of attributes containing information about the position and the residual energy of each particle that enters the detector. The uncertainty of the detector is assumed to be σ = σ = 320 nm for the position measurements and σ/E = 1% for the energy measurement, as discussed in Section 4.B. Gaussian random variables are added to the simulated attributes to simulate the effect of the detector estimation uncertainty.

For reconstruction, we discretize a slab of 1 mm × 1 mm × 50 μm tissue with 1 μm3 cubic voxels. As the position-measurement uncertainties of the detector, σ and σ, are relatively small compared to the voxel size, we approximate the detector response G(x; x̂, σ) with a Dirac-delta function δ(x − x̂). The probability density of measuring  = (x̂, ŷ, Ê) when a particle is emitted from voxel n can be approximated as

where

and

The function G(ℓ; ℓ̂, σ) is a one-dimensional Gaussian function with mean ℓ̂ and standard deviation σ, both of which depend on Ê.

The residual, which is defined as the distance between f̂ and f̂(−1), was calculated for each iteration and used as an indicator for whether the total number of iterations was sufficient. When the residual is getting close to 0, the estimated source distribution gradually converges to f̂, where f̂ maximizes the likelihood for a Poisson data model. Tests we performed showed that the residual does not change significantly after about 100 steps. Reconstruction results with 200 iterations are shown in Figs. 5(c) and 5(d), where the central part of two slices in the 3D reconstructed volume are displayed. The reconstructed objects closely resemble the true decay patterns shown in Figs. 5(a) and 5(b).

In order to compare αET with the conventional imaging method (alpha autoradiography), an autoradiograph is calculated from the same simulated detector data. In autoradiography, the energy of an alpha particle is not measured and the output is an integration of the radioactivity over time at the detector surface. The autoradiograph is calculated by ignoring the energy information and binning the detected particles into 1 μm × 1 μm pixels. The conventional alpha autoradiograph is shown in Fig. 5(e). In comparison to αET, the autoradiograph has worse xy-planar resolution. Furthermore, the autoradiograph in Fig. 5(e) does not show features related with the second object ‘Ding’ [shape shown in Fig. 5(b)]. In other words, information about objects located deeper in tissue is lost when shallower objects are present.

Figures 5(c) and 5(d) shows only 2D slices (on xy-plane) of the 3D reconstructed volume. To further check the reconstructed volume in z-direction, integrations of the true emission pattern and the reconstructed volume in xy-plane are calculated and shown in Fig. 6. The result in Fig. 6(b) shows that the depth information of the radionuclide distribution is well preserved, as radioactivity from one slice does not leak significantly through other slices.

Fig. 6

The integration over xy-plane of the true emission pattern (a) vs. that of the reconstructed objects (b).

5. PARTICLE-PROCESSING DETECTORS

5.A. Concepts

Inspired by research on list-mode data[45,50] and photon-processing detectors,[51,53] a new detector concept — particle-processing detector (PPD) — is introduced. PPDs detect single particles, and for each detected particle, they measure a set of attributes, which may include the interaction position, the propagation direction, and the residual energy of the particle.

Charged-particle detectors are often based on scintillators[54,55] or semiconductors.[56,57] Charged particles impinging on a detector interact with the detector along their tracks. Upon interacting with charged particles, PPDs receive energy from the particle and generate secondaries, more specifically, scintillator-based PPDs generate photons and semiconductor-based PPDs generate electron-hole pairs.[44] The secondaries are produced close to the interaction locations, then propagate to a plane where the signals are collected.

The signals produced by the secondaries may contain information about the interaction locations and the energy deposited. From the interaction locations, information about the impinging position and propagation direction of the charged particle may be extracted. If the particle is stopped in the detector, the energy deposited is the residual energy of the particle.

The particle attributes are estimated from the raw detector signals, e.g., pixel values on a CMOS sensor or a semiconductor detector. The signal processing requires an accurate model specific to each hardware setup. In the signal-processing step, maximum-likelihood estimation (MLE) provides a rigorous way to estimate attributes from raw signals. MLE is a valuable tool for signal processing, because the estimates are asymptotically unbiased, efficient, normally distributed, and consistent.[44]

Several designs of PPDs are described in Supplement A. In the following sections, we demonstrate the feasibility of PPD by presenting experimental results obtained with a two-foil directional charged-particle detector.[37]

5.B. Experimental results: a two-foil detector

5.B.1. System configuration

A two-foil detector[37] aims at measuring position and direction of each detected charged particle. A two-foil setup is introduced in Section 3.B and illustrated in the dashed box in Fig. 2.

A charged particle passes through the phosphor and deposits some of its kinetic energy there. The phosphor produces visible light in the process. The particle interacts with the first phosphor layer which produces light and then passes on to the second layer which also produces light. The light generated at the first layer spreads out more than that generated at the second phosphor layer, and this fact is used to determine at which layer each flash of light was generated. From the image of the light produced by each particle on the light sensor, the two interaction positions can be estimated. If the first phosphor foil is thin enough that a charged particle passing through it does not get significantly deflected, the direction of the particle can be estimated from the two interaction positions and the distance between the two foils.

The ultrathin phosphor foil consists of a layer of 3.5-μm-thick P43 phosphor powder coated on a clear 3-μm-thick Mylar foil (Applied Scintillation Technologies). The two phosphor foils are stacked with an airgap between them with the two phosphor sides facing each other. The air-gap separation between the two phosphor foils is about 1 mm. The image intensifier (ProxiVision), which amplifies the excited light signals to achieve sufficient sensitivity, has an S20 photocathode and a one-stage Micro-Channel Plate (MCP). We couple two 50-mm F/1.2 camera lenses to form an image of the intensifier’s output screen at unity magnification on an ultrafast CCD camera (Phantom V1210, Vision Research Inc.). This camera has a 1280 × 800 pixel array with 28 μm pitch and can operate at a speed of 10,000 fps. The camera comes with an FPGA card which enables realtime signal processing. The light-capturing system is similar to that of the iQID camera.[55,58]

The experimental setup is shown in Fig. 7. The entire system is enclosed in a light-tight box so that the ambient-light background in the box is negligible. We used the setup to measure the position and direction of alpha particles.

Fig. 7

Experimental setup of a two-foil detector (cfr. Fig. 2 for a schematic diagram.) The detector includes: (a) scintillator, (b) image intensifier, (c) lens system with high numerical aperture, and (d) CMOS camera. [Color figure can be viewed at wileyonlinelibrary.com]

5.B.2. Maximum-likelihood estimation

Signals produced when an alpha particle interacts with the two-foil detector are shown in Fig. 9. In the following discussions, we denote the phosphor foil near the source as the first foil and the phosphor foil near the imaging intensifier as the second foil. In each figure, there is a cluster of pixels with large values (black), which is the signal produced by the second foil. The signal produced by the first foil forms a cluster of discrete specks, where each speck is a photoelectron produced in the photo-cathode of the image intensifier. We use maximum-likelihood estimation (MLE)[44] to estimate the positions (x1,y1,x2,y2), where (x1,y1) is the interaction location in the first foil and (x2,y2) is the interaction location in the second foil.

Fig. 9

Signals of one alpha-particle-detection event are shown in Fig. 9, with (a) corresponding to the θ0 = 45° collimator and (b) corresponding to the θ0 = 0° collimator.

If we denote the raw signals on the CMOS camera as g and the true attributes on a detected particle as A, the probability density function for the signals conditioned on a particular attribute vector, which is denoted as pr(g|A), is referred to as a likelihood function. As discussed in Supplement B, if we assume the pixel values are independent and Gaussian distributed, the likelihood function is

where g is the value of the m pixel on the CMOS camera; ḡ(A) is the mean of g; and σ(A) is the standard deviation of g. The functions ḡ(A) and σ(A) are functions of the true underlying attributes A.

An ML algorithm[30,58] was developed to estimate the two interaction locations of each detected particle. The frames with event signals are selected with a threshold on pixel values. A cropped image of size 41×41 pixels centered at the maximum pixel value is fed to a contracting grid algorithm[59] to estimate r̂2, where the circumflex accent indicates estimation results. The whole frame and r̂2 are then fed to another contracting grid algorithm to estimate r̂1.

An early approach performed simultaneous estimation of (x1,y1,x2,y2) by (1) estimating (x2,y2), (2) using the estimated (x2,y2) as an input to estimate (x1,y1), (3) using the estimated (x1,y1) as an input to estimate (x2,y2), and (4) repeating (2) and (3) for a few iterations. However, in the experiment, we discovered that the estimation results of step (3) were not very different from the results of step (1). Therefore, we decided to do the estimations with only the first two steps. Even though the parameters are not simultaneously estimated, the results are not much different from a simultaneous estimation. This is probably due to the fact that the signals (the peak pixel value) produced by the second foil (the one close to the input surface of the image intensifier) is many times stronger than the signals produced by the first foil.

5.B.3. ḡ(A) and σ(A)

The mean and standard deviation of the detector signals are measured with one-foil setups. The final pixel values are sum of the signals produced by each phosphor layer, hence the mean of the final pixel value g is

where ḡ is the mean signals produced by the i phosphor.

We measured ḡ1 and ḡ2 separately with only one phosphor foil at a time. In the experiment, we used a collimated 239Pu source. The collimator is a hole with 1-mm diameter and 1-cm height. The camera was operated at a high frame rate, so that the probability of having more than one event in any given frame is negligible. The frames with single event are select and the ‘shift and add’ algorithm[60] is used to calculate the mean detector responses.

Figure 8(a) shows the experimental mean detector response of the first phosphor foil; and Fig. 8(b) is the measured PSF of the second phosphor foil. Each result is obtained with 250 events. The experimental mean detector response takes into account many sources of blur including the light propagation from the phosphor to the image intensifier, the optics in the image intensifier, the lens system, and the camera pixels.

Fig. 8

Mean detector response of one phosphor foil located (a) about 1 mm away from the image intensifier and (b) on the input surface of the image intensifier. The PSF is calculated from 250 alpha events with ‘shift and add’ algorithm.[60]

As discussed in Supplement B, when the CMOS readout noise is small compared to the total variance of the signals, the variance of the signals can be calculated by

where C is a constant related to the mean and variance of the gain provided by the image intensifier.

5.B.4. Results

We examine the capability of a two-foil detector in measuring the incident direction of alpha particles. To select particles that travel in certain directions, plastic slit collimators (3D printed) are placed between the source and the phosphor foils. The coordinate system is defined so that the z-axis is perpendicular to the image plane and the particle directions are limited in y-direction.

The direction of a particle can be described by two angles θ and θ, where θ is the angle between z-axis and the xz-plane projection of the direction, and θ is the angle between z-axis and the yz-plane projection of the direction. Two slits are used in the experiment, both of which are 250 μm wide and 1 cm thick. Therefore, θ, is limited by the slits to θ0±1.43°, where θ0 = 45° and θ0 = 0° respectively for the two slits.

In the experiment, a 6 kBq 239Pu disk source, which emits 5.15 MeV alpha particles, is placed about 1 cm away from the phosphor screen. The disk source is on a plane that is parallel with the phosphor screen. The diameter of the disk is about 3 mm. With the slit in position, the measured activity is 3.75 Bq. The camera is operated at 100 fps; on average, there are less than four events recorded per 100 frames.

Two signals each from one alpha particle detection event are shown in Fig. 9, with (a) corresponding to the θ0 = 45° collimator and (b) corresponding to the θ0 = 0° collimator. For a particle incident at θ0 = 45°, ŷ1 and ŷ2 are separated; for a particle incident at θ0 = 0°, ŷ1 and ŷ2 are about the same.

In Fig. 10, histograms of θ̂ are shown in (a, b). Based on its definition, θ̂ relates to ŷ1 and ŷ2 as tan for each event, where d is the distance between the two foils. The uncertainty (expressed as standard deviation) of the angular measurement is 14.5° for alpha particles incident at θ0 = 45°, and 13.7° for particles incident at θ0 = 0°.

Fig. 10

Histograms of estimated angle, θ, (a) when a 45° collimator is used, versus (b) when a 0° collimator is used. [Color figure can be viewed at wileyonlinelibrary.com]

A two-foil detector measures both position and direction. If the residual energy is of interest, one can simply change the second phosphor foil with a thick layer of scintillator that stops the particle. However, the performance of the detector is limited by the quantum efficiency of the image intensifier, especially when the incident particle is a beta particle.

6. CONCLUSIONS AND FUTURE WORK

In this paper, we proposed a new approach for charged-particle imaging, which we called charged particle emission tomography (CPET). CPET techniques, such as beta emission tomography (BET) and alpha emission tomography (αET), reconstruct the three-dimensional distribution of charged-particle-emitting radionuclides. A mathematical forward model of CPET is described. Reconstruction results from simulated data are presented for BET with position and direction information, and for αET with position and energy information.

A crucial component of a CPET system is a particle-processing detector (PPD). A PPD processes signals from each event and measures a subset of particle attributes, such as position, direction, and energy. Preliminary results are obtained from a two-foil charged-particle detector. Experiments show that the detector is able to measure the position and direction of incident alpha particles. For future work, the phosphor foil can be changed into thin monolithic scintillation crystals, which exhibits better light uniformity.

Our simulation results show that a PPD that measures more attributes than the position from each detected particle enables CPET. For BET, in addition to position, direction information provides more detail about where the beta particles originated. For αET, in addition to position, energy information makes a three-dimensional reconstruction of an alpha-radioactive distribution possible. For future work, other CPET modalities, including (a) BET with position, direction, and energy information, (b) αET with position and direction information, and (c) αET with position, direction, and energy information can be considered.

Although objective assessment of image quality should be used to evaluate the performance of CPET, we relied on visual inspection of reconstructed data to show the potential of CPET when compared to conventional techniques, such as autoradiography. A task-based image quality assessment of CPET system can be carried out in the future.

Another avenue for future work includes studying the resolution limits of CPET and how it is related to the energy spectra of the particles emitted and the depth of the sources.

Data S1

Particle-processing detector designs.

Data S2

Statistical analysis of two-foil detectors.