Methods for determining position and energy in scintillation detectors
Patent Information
- Application Number
- DE502020011108
- Authority / Receiving Office
- DE · DE
- Patent Type
- Patents
- Current Assignee / Owner
- Priority Date
- 2019-10-15
- Filing Date
- 2020-09-26
- Publication Date
- 2025-06-12
- Estimated Expiration
- 2040-09-26
AI Technical Summary
Existing methods for determining the energy and position of scintillation events in scintillation detectors suffer from significant statistical uncertainties and measurement errors due to Poisson fluctuations and production tolerances, leading to inaccuracies and increased computational time.
An iteration-free maximum likelihood (ML) algorithm is employed to determine the energy and position of scintillation events, which accounts for Poisson statistics and is tolerant to missing data and production tolerances, allowing for real-time processing with minimal hardware requirements.
The ML-based algorithm significantly reduces computational time and maintains high precision in determining energy and position, achieving low dead time and high spatial resolution in scintillation detectors, even at low light levels.
Description
[0001] The invention relates to a method for position and energy determination in scintillation detectors.
[0002] Scintillation detectors are elementary components of a wide variety of particle detectors used in particle and neutrino physics, nuclear medicine (e.g., positron emission tomography (PET), Compton cameras, and single photon computed tomography (SPECT)), radiological imaging, and radiation protection. Scintillation detectors are primarily used to detect particles that can trigger scintillation events, such as gamma photons, α-particles, or β-particles. These can be elementary particles, such as leptons, gamma-ray photons, or X-ray photons, or particles composed of elementary particles, such as mesons, baryons, or ions. A scintillation detector always consists of a scintillator and a photodetector. Scintillators can be in monocrystalline form (e.g. BGO, LSO, etc.), polycrystalline form (e.g. Ultra Fast Ceramics), liquid form (e.g.Xenon) or gaseous form (e.g. high pressure xenon). Solid-state scintillators can be in the form of continuous crystals or as fully or partially segmented crystals. In fully segmented scintillators, the individual scintillator segments are also called scintillator pixels. The individual scintillator segments are normally partially or completely optically separated from one another, e.g. by semi-transparent, opaque or reflective layers. A scintillator consisting of several scintillator segments is called a scintillator array or scintillator matrix. Photomultiplier tubes (PMTs), multi-channel plates (MCPs), avalanche photodiodes (APDs) and silicon photomultipliers (SiPMs) are used as photodetectors. The SiPMs can be implemented using either analog technology (aSiPMs) or digital technology (dSiPMs).
[0003] In many applications of scintillation detectors, it is necessary to determine not only the energy of the particle but also the arrival time of the particle and ideally the three-dimensional position, but at least the two-dimensional position of the photoconversion within the scintillator volume.
[0004] In two-dimensional positioning, the photoconversion position is determined in a plane parallel to the photosensitive surface of the photodetectors. These two coordinates are referred to as the x-coordinate and y-coordinate. The additional third coordinate in three-dimensional photoconversion positioning is generally referred to as the depth of interaction and is referred to as the z-coordinate.
[0005] In most cases, the particle's arrival time is measured by the analysis electronics downstream of the photodetectors, e.g., by threshold discriminators or constant fraction discriminators, or combinations of both. With dSiPMs, the arrival time of individual photons can be measured directly in the photodetector without downstream analysis electronics and made available directly for measurement data processing.
[0006] Position-sensitive photodetectors are required to determine the photoconversion position. Position-sensitive photomultiplier tubes (PSPMTs), micro-channel plates (MCPs), APD arrays, and SiPM arrays are used for this purpose. The latter consist of matrices of normally independent, individual SiPMs or APDs, which are combined into a module through electronic integration. PSPMTs are typically implemented with segmented anodes, a common photocathode, and focusing dynodes, which is why the individual anode elements do not operate independently of each other. Both the anode segments of PSPMTs and MCPs, as well as the individual SiPMs and APDs of SiPM arrays and APD arrays, are referred to as photodetector pixels.
[0007] An incoming gamma photon, also called a primary gamma photon, interacts with the scintillator via the photoelectric effect, pair production, or Compton effect. In the photoelectric effect, the energy of the primary gamma photon is completely transferred to an electron in the scintillator, which then excites the scintillator material. In the Compton effect, only a portion of the energy of the primary gamma photon is transferred to an electron in the scintillator, which then excites the scintillator material. The gamma photon retains the remaining energy and can interact with the scintillator again via the photoelectric effect or Compton effect. This process repeats until an interaction via the photoelectric effect occurs, in which the gamma photon is annihilated, or the gamma photon leaves the scintillator without further interaction. The latter event is called a Compton escape. An event with multiple interactions is called a Compton cascade.In pair production, a positron and an electron are created and the energy of the primary gamma photon is completely transferred to these two particles.
[0008] In contrast to gamma photons, the range of electrons, namely photoelectrons and Compton electrons, and positrons in the scintillator is very short (<_ 150 µm at electron energies of 511 keV). The energy released by the electron or positron to the scintillator excites the scintillator's scintillation centers, which then decay and emit scintillation light isotropically within a short time interval. The number of scintillation photons is approximately proportional to the energy released by the particle, for example a gamma photon, during the interaction. In complete Compton cascades, i.e. the particle is completely converted into scintillation light in the scintillator and no Compton escape occurs, the total number of scintillation photons is therefore approximately proportional to the energy of the primary particle, for example the gamma photon.The process by which energy is transferred from the gamma photon via the photo or Compton electron to the scintillator crystal and converted into scintillation photons is called photoconversion.
[0009] An event in which the energy of exactly one primary gamma photon or particle is converted into scintillation light in one or more photoconversions shall be referred to as a scintillation event.
[0010] The scintillation light emitted in the short time interval, or a signal proportional to it from the photodetector, e.g., voltage, current, or charge, is integrated either directly by the photodetector or by downstream electronics over a defined, constant time interval. The integration is initiated by triggering electronics that compare the rapidly increasing scintillation light intensity at the beginning of the scintillation pulse with a threshold value (threshold discriminator) and start the integration when the threshold is exceeded. This threshold value is chosen to be sufficiently large compared to the thermal noise of the photodetector or other noise sources, thus avoiding continuous triggering of the integration due to noise signals. In digital SiPMs, the integration can be performed directly in the SiPM and by counting the active microcells (also called single avalanche photodiodes (SAPDs)) of the SiPM.
[0011] Due to the isotropic emission of the scintillation light, the scintillator would have to be completely enclosed with photodetectors for complete detection of the scintillation light. For economic reasons and due to technical feasibility, normally only one side of the scintillator is optically coupled to a photodetector. The remaining side surfaces are coated with a reflector, which reflects the scintillation photons so that they reach the photodetector after one or more internal reflections. In both continuous scintillators and segmented scintillators or scintillator arrays, a characteristic scintillation light distribution results in the plane of the optically coupled photodetector, position-sensitive photodetector, or photodetector array due to the isotropic emission of the scintillation photons and the internal reflections at the remaining scintillator surfaces.This scintillation light distribution exhibits a maximum in the xy plane, i.e., the plane parallel to the sensitive surface of the photodetector, at the photoconversion position. This position will be discussed in more detail later (Fig. x PK , y PK ). The greater the distance of the photodetector pixel from the position ( x PK , y PK ) in the xy plane, the smaller the amount of scintillation light detectable by the photodetector pixel becomes. At a sufficiently large distance from ( x PK ,y PK ) the amount of scintillation light can also approach zero.
[0012] From this scintillation light distribution, the energy and photoconversion position of the gamma photon can be determined. By inserting a light guide, which in the simplest case consists of a layer of material transparent to the scintillation light, the scintillation light distribution can be easily varied and adapted to the photodetector pixel size for optimal photoconversion position determination. The most widely used algorithm for determining the energy and photoconversion position is the determination of the expected value, also called the Anger method after its inventor Hal Anger. For a photodetector array with N x Photodetector pixels in x-direction and N y Photodetector pixels in the y-direction can be used to determine the energy 〈 E 〉 Anger and the two coordinates 〈 X 〉 anger and < Y 〉 Anger the photo conversion position according to E Anger = ∑ i x N x ∑ i y N y q i x , i y X Anger = 1 E Anger ∑ i x N x ∑ i y N y x i x ⋅ q i x , i y Y Anger = 1 E Anger ∑ i x N x ∑ i y N y y i y ⋅ q i x , i y determine, where i x = 1, ... , N x the photodetector pixel index in x-direction, i y = 1, ... , N y the photodetector pixel index in y-direction, x ix the center of gravity of the i x -th photodetector pixel index in x-direction, y iy the center of gravity of the i y -th photodetector pixel index in y-direction and q i x , i y which are generated by the photodetector pixels ( i x , i y ) are signals detected, which are proportional to the total area of the photodetector pixels at the position ( x ix , y iy ) and the amount of scintillator light integrated in a defined time interval. Depending on the photodetector used and its electrical output circuit, the q i x , i y These can be analog or digital values, the number of active micro-cells or single avalanche photodiodes (SPAD), charges, voltages, or currents. The method for energy and photoconversion position determination described by formulas 1 - 3 can be improved according to (Chen-Yi & Goertzen, 2013 [1]) by q ix,iy through weighted signals w ix, i y q ix,iy replaced, whereby the w ix,iy are to be determined individually.
[0013] However, the determination of energy and photoconversion position using equations 1 - 3 or [1] has two significant disadvantages. First, it does not take into account that the detection of photons is a Poisson process, which is why the signals q ix,iy follow a Poisson statistic, and the statistical measurement error of q ix,iy proportional to q i x , i y This results in large statistical uncertainties for 〈 X 〉 Anger and < Y〉 anger . Second, photodetector arrays often consist of individual photodetector pixels that operate completely independently of each other. In particular, the triggering electronics for the temporal integration of the signals of each individual photodetector pixel operate independently of all other photodetector pixels in the photodetector array. The same can also apply to MCPs and PSPMTs. In conjunction with the scintillation light distribution, the detectable light quantity for photodetector pixel positions that are far from ( x PK ,y PK ) away, assume a value which is close to the threshold value for the photodetector pixel ( i x , i y ). Due to the Poisson statistics of the signals q ix,iy This leads to the effect that the integration of the signal for photodetector pixels ( i x , i y ) with a large margin of ( x PK , y PK ) is triggered randomly, depending on whether the signal, which is subject to Poisson fluctuations, is above or below the threshold value. If the integration of the signal proportional to the scintillation light is not started, the signal for this photodetector pixel is q ix,iy = 0. This means that the number of photodetector pixels increases with q ix,iy > 0 differs from scintillation event to scintillation event and results in the positions calculated according to formulas 2 and 3 (〈 X 〉 Anger ,〈 Y 〉 Anger ) can have significant statistical errors (Lerche, et al., 2016 [2]). In (Schug, et al., 2015 [3]) this problem is circumvented by using the signals q ix,iy the photodetector pixels without signal, i.e. those for which the integration was not triggered because the integrated scintillation light quantity is below the threshold due to the Poisson fluctuation, are replaced by a value which is calculated from the signals with q ix,iy > 0 for the same scintillation event. However, the extrapolated signal does not correspond to the actual scintillation light quantity for the corresponding photodetector pixel and the signal for a maximum of one photodetector pixel can be q ix,iy = 0 can be extrapolated.
[0014] An alternative method to Equations 1 - 3 for determining the energy and position of scintillation events is the determination of the maximum likelihood (ML) estimate as described in (DeWitt et al., 2010 [4]), (Johnson-Williams et al., 2010 [5]), (Wang et al., 2016 [6]) and [2]. In [4], [5], and [6], iterative ML algorithms for determining the 2D or 3D position of the scintillation event in continuous scintillators are described, which are suitable for implementation in field programmable gate arrays (FPGAs). In the implementations according to [4] and [5], the Poisson distributions of the photodetector pixel signals are approximated by Gaussian distributions. In [2] an iterative ML implementation for determining the 2D position and energy of the gamma photon or other particle in segmented scintillators is described.For this last implementation, it could be shown that the problem described above, in which photodetector pixels with . q ix,iy = 0, can be effectively solved, since ML-based algorithms allow the position and energy of the scintillation event to be determined exclusively from the photodetector pixels with signals q ix,iy > 0, i.e., with incomplete data. Furthermore, using ML-based methods, it is possible to specifically raise the threshold values of individual photodetector pixels. This ensures that fewer photodetector pixels measure signals with values above the corresponding threshold values per scintillation event, and therefore fewer integrations are triggered overall, which reduces the overall dead time of the scintillation detector. If the integration of the signal was triggered by a photodetector element, no further integration can be triggered until the end of the integration and any subsequent data processing steps. During this time, no further scintillation event can be detected. This time is called the detector dead time. In addition, fewer signals need to be transmitted, because signal values q ix,iy = 0 are not forwarded to the data acquisition unit because they contain no information, thus reducing the amount of data to be transmitted. Raising the threshold using the standard procedure according to Equations 1 - 3 is not possible without reducing the precession of the position and energy values.
[0015] The most important disadvantage of the Anger method is the problem already described above, in which, due to Poisson fluctuations, photodetector pixels with q ix,iy = 0, which can lead to significant mispositioning [2]. In addition, the signals q ix,iy are subject to further measurement errors in addition to the Poisson fluctuations. The causes of these additional measurement errors are tolerances in the production of PSPMTs, MCPs, SiPM arrays or APD arrays, which are primarily noticeable through different intrinsic signal amplifications and thus through different signal strengths with essentially the same amount of scintillation light. Further possible production tolerances occur in the positioning of the individual photodetector pixels, ie the photodetector pixel positions ( x ix ,y iy ), which are included as weighting factors in Equations 1-3, are subject to errors due to production tolerances. Production tolerances also occur in the production of continuous scintillators and segmented scintillators, for example, small variations in the scintillator segment size and in the amount of scintillator light per unit particle energy (light yield), variations in the reflectivity of the scintillator surfaces, variations in the transparency of optical couplings, etc. Due to these and other sources of error, the energy and position of the scintillation event determined using Equations 1-3 are subject to errors and require a correction of the energy 〈 E 〉 Anger and the position (〈 X 〉 Anger , 〈 Y 〉 Anger ) following their determination. The correction values must be determined independently for each PSPMT, MCP, SiPM array, APD array, continuous scintillator, and segmented scintillator using a calibration measurement and must be repeated regularly due to aging effects of the components.
[0016] When using ML-based methods to determine the energy and position of the scintillation event, look-up tables (LUTs) containing energy and position reference values are required. Required calibration data can be integrated into these reference value tables by creating a separate reference value table for each scintillation detector. However, the required reference value tables in all previous ML-based methods are so large that it is not possible to store the reference value tables for all scintillation detectors of a complete PET or SPECT scanner in quickly accessible memory (e.g., QDR RAM, UltaRAM, BRAM, and flip-flops in FPGAs and cache in CPUs and GPUs).The reference tables must therefore be stored in external SDRAM or DRAM modules, which, however, has a very detrimental effect on the overall calculation time of the energy and position of the scintillation event due to the significantly lower read speed of these memory types.
[0017] In addition, all previously known ML-based methods for determining the energy and position of the scintillation event are formulated iteratively. This means that the algorithm must be run several times before the final result is obtained. Based on a predefined termination condition, which in most cases evaluates whether a desired precession of the result has been achieved, a decision is made for each individual scintillation event as to whether further iterations are required. As a result of such a definition, iteration numbers depend on the individual scintillation event, which adversely affects FPGA implementability and the overall computation time. Alternatively, an average optimal iteration number can be defined in advance. This improves FPGA implementability but leads to many calculations with suboptimal iteration numbers.All previously known ML-based methods [4], [5], [6] and [2] are based on an iterative formulation of the algorithm and are therefore too slow to process all scintillation events of a typical PET or SPECT scanner in a reasonable time and with adequate hardware expenditure. In a state-of-the-art human whole-body PET scanner, depending on the organ being examined and the radiopharmaceutical used, between 2 and 4 million coincident scintillation events per second can occur, for which the energies and positions must be determined. In a dedicated organ-specific PET scanner, e.g. for the breast or head, this rate can be twice as high. Due to the required high data processing speed, Anger-based methods are therefore preferred in human PET scanners.
[0018] Furthermore, state-of-the-art PET scanners initially select coincident scintillation events. Scintillation events for which no coincident scintillation events are detected, so-called single events, are ignored. This significantly reduces the amount of data to be processed. However, it complicates other necessary corrections, such as random coincidence correction and scatter correction. These are more precise and easier to determine when all coincidence events and single events are processed. In a state-of-the-art human whole-body PET scanner, depending on the organ being examined and the radiopharmaceutical used, between 40 and 80 million single scintillation events per second can occur. Depending on the organ being examined and the radiopharmaceutical used, the energies and positions for these events must be determined to enable single-based random coincidence correction and scatter correction.
[0019] The ML-based algorithm for continuous scintillators presented in [5] can process up to 117,000 scintillation events per second per FPGA. Consequently, a coincidence processing platform would require up to 4x10 6 < / 117,000 = 44 FPGAs of the type mentioned in this study. A single-processing platform would require up to 80x10 6 < / 117,000 = 684 FPGAs. Therefore, this implementation would not allow for a cost-effective data processing platform.
[0020] The ML-based algorithm for continuous scintillators presented in [4] can process up to 360,000 scintillation events per second per FPGA. Consequently, a coincidence processing platform would require up to 4x10 6 < / 360,000 = 11 FPGAs of the type mentioned in this study. A single-processing platform would require up to 80x10 6 < / 360,000 ≈ 223 FPGAs. Therefore, this implementation would not allow for the construction of a cost-effective data processing platform.
[0021] The ML-based algorithm for continuous scintillators presented in [6] can process up to 15x10 6< scintillation events per second and FPGA. Consequently, a coincidence processing platform would require up to 4x10 6< / 15x10 6< = 1 FPGAs of the type mentioned in this study. For a single-processing platform, up to 80x10 6< / 15x10 6< ≈ 6 FPGAs would be required. With this implementation, it would therefore be possible to build a cost-effective data processing platform. However, in this implementation, the projections of the signals are first q ix,iy on the x- and y-axis according to formulas 4 and 5 Q i x = ∑ i x N y q i x , i y Q i y = ∑ i x N x q i x , i y In order to use the calculated Q ix and Q iy the energy and position of the scintillation events can be determined with sufficient precision, it must be repeated for all q ix,iy > 0. Therefore, the individual threshold of the photodetector pixels must be set low enough to trigger integrations for all photodetector pixels. As described above, this leads to a significantly increased detector dead time.
[0022] The ML-based algorithm for segmented scintillators presented in [2] can process up to 840,000 scintillation events per second in a multi-CPU (central processing unit) system with 40 threads. Consequently, a coincidence processing platform would require up to 4x10 6 < / 1840,000 = 5 data processing systems of the type mentioned in this study (see (Schug, et al., 2016 [7]) and (Goldschmidt, et al., 2015 [8]). For a single-processing platform, up to 80x10 6 < / 840,000 = 95 data processing systems would be required. An FPGA implementation of this variant has not been proposed and, as in the other ML-based methods mentioned, proves to be very difficult due to the high demand for fast memory access. With this implementation, it would therefore not be possible to build a cost-effective data processing platform.
[0023] The publication by Ch. W. Lerche et al., "Fast circuit topology for spatial signal distribution analysis," 17th REAL-TIME CONFERENCE - IEEE-NPSS TECHNICAL COMMITTEE ON COMPUTER APPLICATIONS IN NUCLEAR AND PLASMA SCIENCE; 2010, PISCATAWAY, NJ, USA, January 1, 2010 (2010-01-01), pages 1-8, DOI: 10.1109 / RTC.2010.5750391, discloses a method for position and energy determination in scintillation detectors, in which the photoconversion energy and the photoconversion position of particles that trigger scintillation events are determined from the distribution of the released scintillation light from one or more scintillation events, as scanned by a photodetector, in an iteration-free process.
[0024] The object of the invention is to overcome the disadvantages of the prior art. In particular, an accurate and rapid method for position and energy determination in scintillation detectors for medical and molecular imaging is to be enabled in order to enable PET cameras, SPECT cameras, Compton cameras, and scintigraphy cameras with high spatial resolution, low dead time, and reasonable hardware performance requirements for the data processing unit. The method should be tolerant to missing data and to the Poisson fluctuations in the signal that typically occur in scintillation detectors. Tolerant here means that the positioning error due to the missing data is so small that it does not lead to artifacts or increased image noise in the image reconstructed from the data. The method should take into account that the detection of photons is a Poisson process, which is why the signals q ix,iy follow a Poisson statistic and the statistical measurement error of q ix,iy proportional to q i x , i y The implementation of the method in CPUs or FPGAs should be so resource-efficient that all scintillation events detected in a PET, SPECT, or scintigraphy examination can be positioned in real time with only a few CPUs and / or FPGAs while data is being acquired. At low light levels, large statistical uncertainties for individual photodetector pixels should not affect the precession of the position and energy values. Raising the threshold value of individual photodetector pixels should not result in a reduction in the precession of the position and energy values. Measurement errors and incorrect positioning of incoming particles should be reduced. Production tolerances in the manufacture of sensors should lead to lower inaccuracies in the determination of position and energy values. Subsequent correction of the energy and position as in the Anger method (〈 E 〉 Anger , and (< X 〉 Anger , 〈 Y 〉) Anger )) should no longer be necessary. Computation times for determining energy and position should be minimized, and energy and position should be determined using an iteration-free method.
[0025] Starting from the preamble of claim 1, the object is achieved according to the invention with the features of the characterizing part of claim 1.
[0026] The method according to the invention overcomes the disadvantages of the prior art. In particular, an accurate and rapid method for position and energy determination in scintillation detectors for medical and molecular imaging is provided, enabling PET cameras, SPECT cameras, Compton cameras, and scintillation cameras to be provided with high spatial resolution, low dead time, and reasonable hardware performance requirements for the data processing unit. The method is tolerant of missing data and the Poisson fluctuations in the signal that typically occur in scintillation detectors. The method takes into account that the detection of photons is a Poisson process, which is why the signals q ix,iy follow a Poisson statistic, and the statistical measurement error of q ix.iy proportional to q i x , i y The method allows implementation in CPUs or FPGAs, which is so resource-efficient that all scintillation events detected in a PET, SPECT, or scintigraphy examination can be positioned in real time with just a few CPUs and / or FPGAs, even during data acquisition. With low light levels, large statistical uncertainties for individual photodetector pixels are prevented from affecting the precession of the position and energy values. Higher threshold values for the individual photodetector pixels can be realized without leading to a reduction in the precession of the position and energy values. Measurement errors and incorrect positioning of incoming particles that can trigger scintillation events are reduced. Particles that can trigger scintillation events include gamma or X-ray photons, α particles, or β particles.In principle, these can be elementary particles, such as leptons or gamma or X-ray photons, or particles composed of elementary particles, such as mesons, baryons, or ions. These are referred to as particles below. Production tolerances during sensor production lead to lower inaccuracies in the determination of position and energy values. Additional energy correction is not required. Computation times for evaluating the measurement results are minimized. The energy and position are determined using an iteration-free method.
[0027] Advantageous further developments of the invention are specified in the subclaims.
[0028] In the following, the invention is described in its general form, without this being to be interpreted in a restrictive manner.
[0029] According to the invention, a method for position and energy determination in scintillation detectors is provided, in which an iteration-free algorithm for determining the energy and position of the scintillation event is determined according to equations (6), (7), and (8). According to the invention, a scintillation event can be triggered by particles that can trigger scintillation events, for example, gamma photons, X-ray photons, α-particles, or β-particles. In principle, these can be elementary particles, such as leptons or photons, or particles composed of elementary particles, such as mesons, baryons, or ions. L ˜ m j q n 1 ⋯ q n t = ∑ i = 1 t log 2 ˜ μ m j , n i ⋅ q n i − max q n i ⋅ ∑ i = 1 t μ m j , n i m ML = argmax m j ∈ m 1 ⋯ m p L ˜ m j q n 1 ⋯ q n t E ML = norm m ML ∑ i = 1 t q n i / ∑ i = 1 t μ m ML , n i
[0030] In equations (6) - (8) the following mean: : log-likelihood for the signals { q n 1 ,···, q nt } and a scintillation event in the crystal with the index m j n i ∈ { n 1 ,···, n t }: Indices of the photodetector pixels for which q ni > q th and t ≤ N applies. q th denotes the threshold value set for the photodetector pixels and t the number of photodetector pixels with a signal above the threshold value q th . m j ∈ (m 1 ,···,mp}: Indices of the scintillator segments that are taken into account for the calculation. m ML : Index of the scintillator segment in which scintillation most likely occurred. q ni : Signals from the photodetector array for the photodetector pixels n i . E ML : Most probable total energy of the scintillation event as determined by ML algorithm. log 2 ˜ : Approximation for the logarithm to base 2 norm mML : Calibration factor for correct calculation of energy for scintillator segment m ML . µ mj,ni : Probabilities that a scintillation photon, which is in scintillator segment m j emitted in the photodetector pixel n i is detected. µ mML,ni : How µ mj,ni , however for m j = m ML ;
[0031] The photoconversion energy and photoconversion position of particles triggering scintillation events are calculated from the distribution of the released scintillation light from one or more scintillation events scanned by a photodetector in an iteration-free procedure according to formulas (6), (7) and (8).
[0032] The algorithm used in the invention according to equations (6) to (8) eliminates iteration, which reduces computing time and leads to high data processing rates. The Poisson statistics underlying the scintillation signals are taken into account, and the position of the photodetector pixels in the xy plane can be chosen arbitrarily and does not have to lie on a Cartesian grid. The method is robust to incomplete data, which is why photodetector pixels without a signal do not pose a problem, and thus a sufficiently short scintillation detector dead time can be achieved while simultaneously maintaining high precession of the determined energy and position values of the scintillation events.
[0033] In equations (6) - (8) m ML the index of the scintillator segment in which the scintillation most likely occurred, E ML the most probable total energy of the scintillation event, q ni the signals from the photodetector array, where only for { n 1 ,···, n t } Photodetector pixel signals with q ni > q th present and t ≤ N applies. t be different for each scintillation event and the size of t can be influenced by the threshold setting for the photodetector pixels. Small values of t Between 5 and 20 are advantageous for a photodetector pixel size of (3 - 5 mm) 2< and a scintillator segment size of (1 - 3 mm) 2< for fast calculation according to equations (6) - (8). Since the order of the photodetector pixel indices is irrelevant for the calculation with ML-based methods, the numbering of the photodetector pixels does not have to reflect their geometric arrangement. The photodetector pixels can be arranged arbitrarily in the plane of the photodetector array for the algorithm, and in particular, for a Cartesian arrangement, N x = N y In equations (6) - (8) m j ∈ { m 1 ,···, m p } are the indices of the scintillator segments that are taken into account for the calculation according to equations (6) - (8). In the majority of cases, the scintillation event takes place entirely in only one scintillator segment, also called scintillator pixel. The further this scintillator segment is from the photodetector pixels for which a signal q ni > 0, the less likely this scintillator segment is to emit scintillation light. Consequently, a ranking of the scintillator segments according to their distance d from the center of the scintillation light distribution in the xy plane. The center of the scintillation light distribution is determined by the position of the photodetector pixel with maximum signal q ni These rankings can even be determined in advance for each photodetector pixel and stored in a look-up table (LUT) with a size of N ⋅ M rel ⋅ log 2 M Bits are stored. M rel Number of relevant scintillator segments. This can be in the range 1 ≤ M rel ≤ M can be freely chosen. Larger values of M rel lead to more precise results but longer processing times. In equations (6) - (8) log 2 ˜ an approximation for the logarithm to base 2. The approximation of the logarithm can be done, for example, as in (Gutierrez & Valls, 2010 [9]), where a very low accuracy with a mean relative error of 2%, a mean absolute error of 0.11 and its own maximum absolute error of 0.17 is sufficient. In equations (6) - (8), norm mML a calibration factor for the correct calculation of the energy and µ mj,ni the probabilities that a scintillation photon, which is in scintillator segment m j emitted in the photodetector pixel n i is detected.
[0034] The probabilities µ m,n are determined in advance by measurement, simulation or calculation and stored in a look-up table (LUT) with a size of M · N · P Bit. P denotes the required precession of the probability values, M the total number of scintillator segments used in the scintillation detector, and N the total number of photodetector pixels used in the scintillation detector. P depends on the detector type and should be greater than 8 bits. The calibration factors norm m must be determined in advance by measurement, simulation or calculation and stored in a look-up table (LUT) with a size of M · P Bit can be stored. µ m,n and norm m can be calculated as follows from the measured light distributions averaged over several scintillation events Î m,n can be determined according to equations (9) and (10). Î m,n the average light intensity for the photodetector pixel n when the scintillation in the scintillator segment m takes place. μ m , n = I ^ m , n / ∑ n = 1 N I ^ m , n norm m = ∑ n = 1 N I ^ m , n / max I ^ m , n
[0035] In equations (9) and (10) the following mean: N = N x · N y : Total number of photodetector pixels in the scintillation detector, where N x the number of photodetector pixels in the x-direction and N y is the number of photodetector pixels in the y-direction.
[0036] norm m : Calibration factors for correct calculation of energy.
[0037] Î m,n : average light intensity for the photodetector pixel n when the scintillation takes place in the scintillator segment m.
[0038] The required LUTs can be stored in external Dynamic Random Access Memory (DRAM), Synchronous Random Access Memory (SRAM), Quad Data Rate (QDR), SRAM or memory modules with comparable performance.
[0039] For the complete determination of E ML and m ML starting from a set of t signals { q n 1 ,···, q nt } from a scintillation detector with M scintillator segments and a photodetector array with N Photodetector pixels, where for the signals [ q n 1 ,···, q nt } applies: q ni > 0 ∀ i ∈ 1, ... , t ≤ N , the following calculation steps are required: 1. Identify the photodetector pixel index n max with the maximum signal q n max . Are there multiple pixels with maximum signal q nmax You can either choose only one or continue with both. n max as well as q nmαx are buffered in registers of the FPGAs or CPUs. 2. From the LUT, in which the scintillator segment indices are arranged in descending order of their distance d from the position of the photodetector pixel n max are stored the 1 ≤ p ≤ M most relevant scintillator segment indices { m 1 , ···, m p } and buffered in memory cells of the FPGAs or CPUs. 3. From the LUT, in which the detection probabilities µ m,n are stored, the { m 1 ,···, m p } × { q n 1 ···, q nt } relevant probabilities µ m j , n i with i ∈ 1, ... , t and j ∈ 1, ... ,pread out and temporarily stored in memory cells of the FPGAs or CPUs. 4. The approximated logarithms log 2 ˜ μ m j , n i are determined with the q ni and summed according to equation (6) and buffered in memory cells of the FPGAs or CPUs. 5. The probabilities µ mj,ni are summed according to equation (6) and the sum is multiplied by max( q ni ) and buffered in memory cells of the FPGAs or CPUs. 6. The scintillator segment index m ML , for which the sum ∑ i = 1 t log 2 ˜ μ m j , n i . q n i − max q n i ⋅ ∑ i = 1 t μ m j , n i is largest is identified and buffered in a memory cell of the FPGAs or CPUs. 7. The probabilities µ mML,ni are summed according to equation (8) and the result is buffered in a memory cell of the FPGAs or CPUs. 8. From the LUT, in which the calibration factors norm m for the correct calculation of energy, norm mML read out and with the sum of the photodetector pixel signals ∑ i = 1 t q n i multiplied and divided by the sum of the probabilities µmML,ni divided. For implementations in FPGAs, it makes sense to outsource the division to the image reconstruction processor, since division in FPGAs requires a lot of resources. The additional data volume required to transfer the dividend and divisor instead of just the quotient is negligible.
[0040] The calculation according to steps 1 - 8 is not iterative. Divisions are not absolutely necessary and multiplications are minimized. The required memory space with very fast access (e.g. cache in CPU, flip-flops, UltraRAM, BRAM, or comparable in FPGAs) is minimized according to the invention to such an extent that all required data can be accommodated in commercially available FPGAs and CPUs. In addition, the required data transfer of data that cannot be accommodated in the CPU cache or in FPGA flip-flops or in FPGA, BRAM or FPGA UltraRAM is minimized. For the calculation of the logarithm, a very fast, approximate implementation can be chosen, since high precision is not required for the estimation of E ML and m ML is required. The calculation according to steps 1 - 8 is much more robust and precise than the implementation of the Anger method (equations (1) - (3)). The calculation according to steps 1 - 8 is significantly faster compared to all cited, alternative ML-based methods. In particular, with state-of-the-art high-end CPUs, the execution of steps 1 - 8 for 5 million scintillation events is possible in one second, which is why only 16 threads are necessary for the previously mentioned 80 x 10 6< single scintillation events. With an FPGA implementation of the calculation steps 1 - 8, processing of the 80 x 10 6< single scintillation events is possible with only 4 high-end FPGAs.
[0041] The formulation of the ML-based algorithm allows in particular an effective use of the parallelization possibilities in CPUs (duplication) and the parallelization possibilities in FPGAs (duplication and pipeline), as in Figur 3 Described. A parallelization of the multiplications and the calculation of the log 2 ˜ μ m j , n i is of fundamental importance for a sufficiently fast and accurate calculation.
[0042] The described ML-based algorithm can also be used with continuous scintillators by dividing (quantizing) the three-dimensional, continuous scintillator volume into a finite number of subvolumes. For example, let the three-dimensional, continuous scintillator volume be of dimensions H × B × T, the height can be M H Intervals of length H / M H , the width in M B Intervals of length B / M B and the depth in M T Intervals of length T / M T These three-dimensional intervals are then treated as individual scintillation segments. The calculation is identical to that for actually segmented scintillators.
[0043] The figures show detectors and units for determining the energy and position of particles in scintillation detectors in schematic form: It shows: Fig. 1: A scintillation detector Fig. 2: A signal generation in a scintillation detector Fig. 3: A parallelized unit for calculating the approximate logarithms and for multiplication Fig. 4: Unit for determining the ML estimate for the energy and position of the scintillation event Fig. 5: Unit for determining the ML estimate for the energy and position when using FPGAs and CPUs simultaneously Fig. 6: Unit for determining the ML estimate for the energy and position when using CPUs
[0044] Figur 1 shows a typical structure of a scintillation detector with multiple layers of segmented scintillators (1), (2). One to four layers of segmented scintillators are possible. The bottom layer of segmented scintillators (2) is coupled to the photodetector array (4) via a light guide (3), which in this simple case consists of a plane-parallel material layer transparent to the scintillation light. The photodetector array can be a PSMPT, an MCP, a SiPM array, or an APD array. Typical thicknesses of the light guide are 0.1 mm - 2 cm, depending on the detector size and the granularity of the scintillator and the photodetector array.
[0045] In Figur 2 Identical device components have the same reference numerals as in the previous figures. This shows a representation of a single-layer segmented scintillator. The functionality of multi-layer scintillation detectors for three-dimensional photoconversion position determination for measuring depth of interaction with multi-layer segmented scintillators is analogous. The scintillation light (5) from a single scintillator segment is distributed (6) via the light guide (3) over the entire sensitive area of the photodetector array (4). Depending on the threshold setting q th the photodetector pixels of the photodetector array (4), an integration is then triggered and the signals q ix,iy > q th (7) are provided by the photodetector or the downstream electronics. Signals with q ix,iy < q th are not used to calculate energy and position.
[0046] Figur 3 shows a combined, parallelized multiplication unit consisting of several individual multiplication units in pipeline operation (8). In the memory unit (9), e.g., QDR, DRAM, SRAM, etc., the detection probabilities µ m , n permanently stored. In this representation, the column address of a single detection probability µ m , n the index n of the photodetector pixel and the row address the index m of the scintillation segment. Implementations with different assignments are also possible. The relevant detection probabilities µ mj,ni are stored in buffers (11) (e.g. UltaRAM, BRAM, FlipFlops, Cache, etc.) and in a unit (12) for storing and calculating the log 2 ˜ μ m j , n i The photodetector pixel signals are read in via a data interface (13) and p copies of the t signals {q n 1 ,···, q nt } are stored in buffers (14), with each memory unit (15) containing exactly one signal value. At the outputs (16), the values for ∑ i = 1 t log 2 ˜ μ m j , n i ⋅ q n i be read out.
[0047] In Figur 4 In this figure, the same device components have the same reference numerals as in the previous figures. It shows a determination of the ML estimate for the energy and position of the scintillation event. (13) denotes a data slice speed over which the t signals { q n 1 ,···, q nt } are received. (17) is a unit for determining the photodetector pixel index n i with the maximum photodetector pixel value q ni and the maximum photodetector pixel value q ni . With (20) a memory unit, e.g. QDR, DRAM, SRAM, is used to permanently store the indices of the p most relevant scintillator segment indices { m 1 ,···, m p } for the photodetector pixel index n i The reference numeral (19) denotes an optional unit for determining the log 2 ˜ μ m j , n i Values, if the log 2 ˜ μ m j , n i Values not in a combined, parallelized unit according to Figur 3 be determined (20): Unit for determining the sums ∑ i = 1 t μ m j , n i . (18) is a combined parallelized multiplication unit, as shown in Figur 3 to calculate equation (6). (21) is a unit for determining the scintillator segment index m ML with the highest likelihood. If several scintillator segments have the same likelihood, one of the scintillator segment indices with the highest likelihood is selected. (22) denotes a memory unit, e.g., QDR, DRAM, SRAM, for permanently storing the calibration factors. norm m . (23) is a unit for calculating the sum of the probabilities µ mML,ni . (24) is a unit for calculating the sum of the photodetector pixel values q ni . (25) represents a unit for calculating the quotient of the sum of the photodetector pixel values q ni and the sum of the probabilities µ mML,ni .dar. (26) is an output of the scintillator segment index m ML with highest likelihood. (29) is an output of the most probable energy E ML Alternatively, the sum of the probabilities µ mML,ni via the output (28) and the sum of the photodetector pixel values q ni be output via the output (27) and the division (25) is moved to a downstream CPU (in the case of FPGA-based implementation).
[0048] Figur 5 shows the implementation of the determination of the ML estimate for the energy and position of the scintillation event from the scintillation detector (30) using an FPGA unit (31) and a CPU unit (32). The CPU unit (32) is required for further calculations, such as the confidence search and image reconstruction.
[0049] Figur 6 shows the implementation of the ML estimate for the energy and position of the scintillation product from the scintillation detector (30) using only a CPU unit (32). The CPU unit (32) is required for further calculations, such as coincidence search and image reconstruction.
[0050] The invention can be used, for example, for a scintillation detector for PET or SPECT or scintigraphy or Compton cameras, consisting of a single-layer, segmented scintillator, a light guide and a photodetector array (PSPMT, MCP, APD array, SiPM array) and an electronics with FPGA and memory, wherein in the FPGA the calculation steps 1 - 8 and the multiplication unit are implemented as in Figuren 3 und 4 are implemented.
[0051] Likewise, the use in a scintillation detector for PET or SPECT or scintigraphy or Compton cameras, consisting of a multi-layer, segmented scintillator, a light guide and a photodetector array (PSPMT, MCP, APD array, SiPM array) and an electronics with FPGA and memory, where in the FPGA the calculation steps 1-8 and the multiplication unit are implemented as in Figuren 3 und 4 implemented are possible.
[0052] Furthermore, the invention can be used in a scintillation detector for PET or SPECT or scintigraphy or Compton cameras, consisting of a single-layer, continuous scintillator and a photodetector array (PSPMT, MCP, APD array, SiPM array) and an electronics with FPGA and memory, wherein in the FPGA the calculation steps 1 - 8 and the multiplication unit are implemented as in Figur 3 are implemented.
[0053] A further application of the invention is in a scintillation detector for PET or SPECT or scintigraphy or Compton cameras, consisting of a single-layer, continuous scintillator, a light guide and a photodetector array (PSPMT, MCP, APD array, SiPM array) and an electronics with FPGA and memory, wherein in the FPGA the calculation steps 1 - 8 and the multiplication unit are implemented as in Figures 3 and 4 implemented are possible.
[0054] In the last four applications, the invention can be applied in a way in which all calculation steps 1 - 8 are implemented in a CPU and not in an FPGA.
[0055] In the last five applications, the invention can be applied in a manner in which photodetector arrays are mounted on more than one side of the scintillator. In segmented scintillators, the top and bottom sides can be aligned as shown in Figures 1 and 2 can be used to read the scintillation light with photodetectors. With continuous scintillators, all six sides can be used to read the scintillation light with photodetectors. Example:
[0056] Crucial to enabling sufficiently high processing rates for the single scintillation events and coincidence scintillation events occurring in a typical PET scanner is the use of an iteration-free algorithm, as this enables efficient implementation in FPGAs and the use of the resulting parallelization possibilities (e.g., processing pipelines and duplication of processing instances). The use of an ML-based algorithm is preferable because it takes into account the Poisson statistics underlying the scintillation signals and allows the positions of the photodetector pixels in the xy plane to be chosen arbitrarily, and not necessarily on a Cartesian grid as with the Anger method.
[0057] The use of an ML-based algorithm is also preferable because ML-based algorithms are robust against incomplete data, which is why photodetector pixels without signal are not a problem, and thus a sufficiently short scintillation detector dead time can be achieved while simultaneously maintaining high precession of the determined energy and position values of the scintillation events. In order to achieve short dead times, it is also advantageous to use segmented scintillators for large scintillation detectors, since then the scintillation light cannot spread throughout the entire detector volume and the photodetector pixels with values q ix,iy > q thare limited in number and location in the xz plane. This allows for individual operation of the photodetector pixels, allowing multiple independent scintillation events to be read out in a scintillation detector. This significantly reduces the dead time of the entire scintillation detector. For a scintillation detector consisting of a photodetector array (PSPMT, SiPM array, APD array) with N = N x · N y Photodetector pixels and a single or multi-layer segmented scintillator with M = ∑ l M l scintillator segments and M l = M l,x · M l,y Scintillator segments capable l An iteration-free ML-based algorithm for determining the energy and position of the scintillation event can be specified as follows: L ˜ m j q n 1 ⋯ q n t = ∑ i = 1 t lo g 2 ˜ μ m j , n i ⋅ q n i − max q n i ⋅ ∑ i = 1 t μ m j , n i m ML = argmax m j ∈ m 1 ⋯ m p L ˜ m j q n 1 ⋯ q n t E ML = nor m m ML ∑ i = 1 t q n i / ∑ i = 1 t μ m ML , n i
[0058] Here, M l , x the number of scintillator segments in the x-direction in the position land M l,y the number of scintillator segments in the y-direction in the position l For single-layer scintillation detectors (only one layer of scintillator segments), the layer index is omitted l . Cited literature:
[0059] [1] Chen-Yi, L. & Goertzen, A., 2013. Improved event positioning in a gamma ray detector using an iterative position-weighted centre-of-gravity algorithm. Physics in Medicine & Biology, 58(14), p. 189. [2] Lerche, C. W. et al., 2016. Maximum likelihood positioning and energy correction for scintillation detectors. Physics in Medicine & Biology, 61(4), p. 1650. [3] Schug, D. et al., 2015. Data Processing for a High Resolution Preclinical PET Detector Based on Philips DPC Digital SiPMs. IEEE TRANSACTIONS ON NUCLEAR SCIENCE, 62(3), p. 669. [4] DeWitt, D. et al., 2010. Design of an FPGA-based algorithm for real-time solutions of statistics-based positioning. IEEE transactions on nuclear science, 57(1), pp. 71-77. [5] Johnson-Williams, N. et al., 2010. Design of a Real Time FPGA-Based Three Dimensional Positioning Algorithm. IEEE Transactions on Nuclear Science, 58(1), pp. 26-33. [6] Wang, Y. et al., 2016.An FPGA-Based Real-Time Maximum Likelihood 3D Position Estimation for a Continuous Crystal PET Detector. IEEE Transactions on Nuclear Science, 63(1), pp. 37-43. [7] Schug, D. et al., 2016. Initial PET performance evaluation of a preclinical insert for PET / MRI with digital SiPM technology. Physics in Medicine & Biology, Volume 61, p. 2851-2878. [8] Goldschmidt, B. et al., 2015. Software-based real-time acquisition and processing of PET detector raw data. IEEE transactions on biomedical engineering , 63(2), pp. 316-327. [9] Gutierrez, R. & Valls, J., 2010. Low cost hardware implementation of logarithm approximation. IEEE Transactions on Very Large Scale Integration (VLSI) Systems, 19(12), pp. 2326-2330.
Claims
1. A method for position and energy determination in scintillation detectors (30), in which the photoconversion energy and photoconversion position of particles triggering scintillation events are calculated from a distribution of the scintillation light released by a scintillation event or multiple scintillation events sampled by a photodetector in an iteration-free process, characterized in that the photoconversion energy and photoconversion position of particles are calculated according to equations (6), (7) and (8), L ¯ m j q n 1 ⋯ q n t = ∑ i = 1 t log 2 ˜ μ m j , n i ⋅ q n i − max q n i ⋅ ∑ i = 1 t μ m j , n i m ML = argmax m j ∈ m 1 ⋯ m p L ˜ m j q n 1 ⋯ q n t E ML = norm m ML ∑ i = 1 t q n i / ∑ i = 1 t μ m ML , n i wherein L̃: log-likelihood for the signals {qn1, ..., qnt} and a scintillation result in the crystal of the index mj, N = Nx · Ny: total number of photodetector pixels in the scintillation detector, wherein Nx is the number of photodetector pixels in the x-direction and Ny is the number of photodetector pixels in the y-direction, ni ∈ {n1, ..., nt}: indices of the photodetector pixels for which qni > qth and t ≤ N, qth denoting the threshold value set for the photodetector pixels, and t denoting the number of the photodetector pixels with a signal above the threshold value qth, nmax: photodetector pixel indices with the maximum signal qnmax, M = ΣlMl: total number of the scintillator segments in the scintillation detector, wherein Ml = Ml,x · Ml,y is the number of the scintillator segments in the location l of a multi-layer scintillation detector, Ml,x is the number of the scintillator segments in the x-direction in the location 1 and Ml,y is the number of the scintillator segments in the y-direction in the location l, the location index 1 being omitted for single-layer scintillation detectors (30) (only one layer of scintillator segments), Mrel: number of relevant scintillator segments in the range of 1 ≤ Mrel ≤ M, which are selected freely, mj ∈ {m1, ..., mp}: indices of the relevant scintillator segments taken into account in the calculation according to equations (6), (7) and (8), wherein 1 ≤ p ≤ Mrel, mML: index of the scintillator segment in which the scintillation has likely occurred, qni: signals from the photodetector array (4) for the photodetector pixel ni, EML: most probable total energy of the scintillation result upon determination by means of the algorithm according to equations (6), (7) and (8), log 2 ˜ : approximation for the logarithm to base 2, normm: calibration factors for correct calculation of the energy, normmML: calibration factor for correct calculation of the energy for scintillator segment mML, µmj,ni: probabilities that a scintillation photon emitted by scintillation segment mj is detected in the photodetector pixel ni, µmML,ni: like µmj,ni but for mj = mML, wherein this algorithm is implemented in an FPGA (31) or CPU (32).
2. The method according to claim 1, characterized in that one or more photodetector pixel indices nmax of the photodetector with the maximum signals are identified, and in that the 1 ≤ p ≤ Mrel most relevant scintillator segment indices {m1, ···,mp} are read from value tables storing the scintillator segment indices in descending order based on their distance d from the position of a photodetector pixel nmax, and buffered in storage cells of FPGAs (31) or CPUs (32) and used for calculating the log-likelihood according to equation (6).
3. The method according to claim 2, characterized in that the {m1, ..., mp} × {qn1, ···,qnt} most relevant probabilities µmj,ni with i ∈ 1, ..., t and j ∈ 1, ..., p are read from the value table storing the detection probabilities µm,n, and buffered in registers of FPGAs (31) or CPUs (32), and used for calculating the log-likelihood according to equation (6).
4. The method according to any one of claims 1 - 3, characterized in that the approximated logarithms log 2 ˜ μ m j , n i are determined, buffered in storage cells of FPGAs (31) or CPUs (32), multiplied with the qni, and summed up according to equation (6).
5. The method according to any one of claims 1 - 4, characterized in that the probabilities µmj,ni according to equation (6) are summed up, and the sum is multiplied with max(qni), and the results are buffered in storage cells of FPGAs (31) or CPUs (32).
6. The method according to any one of claims 1 - 5, characterized in that the scintillator segment index mML for which the sum of ∑ i = 1 t log 2 ˜ μ m j , n i ⋅ q n i − max q n i ⋅ ∑ i = 1 t μ m j , n i is greatest is identified, and the value is buffered in storage cells of FPGAs (31) or CPUs (32).
7. The method according to any one of claims 1 - 6, characterized in that the probabilities µmML,ni according to equation (6) are summed up, and the sum is buffered in storage cells of FPGAs (31) or CPUs (32).
8. The method according to any one of claims 1 - 7, characterized in that normmML is read from a value table storing the calibration factors normm for correct calculation of the energy, and multiplied with the sum of photodetector pixel signals ∑ i = 1 t q n i , and divided by the sum of the probabilities µmML,ni according to equation (8).
9. The method according to claim 8, characterized in that in the case of implementations in FPGAs (31), the division is outsourced to an image reconstruction computer.
10. The method according to any one of claims 1 to 9, characterized in that the scintillation events are triggered by particles, consisting of the group of gamma photons, α-particles, β-particles, leptons, X-ray photons or particles composed of elementary particles, such as mesons, baryons or ions.