|Home | About | Journals | Submit | Contact Us | Français|
We are developing a dual panel breast-dedicated positron emission tomography (PET) system using LSO scintillators coupled to position sensitive avalanche photodiodes (PSAPD). The charge output is amplified and read using NOVA RENA-3 ASICs. This paper shows that the coincidence timing resolution of the RENA-3 ASIC can be improved using certain list-mode calibrations. We treat the calibration problem as a convex optimization problem and use the RENA-3’s analog-based timing system to correct the measured data for time dispersion effects from correlated noise, PSAPD signal delays and varying signal amplitudes. The direct solution to the optimization problem involves a matrix inversion that grows order (n3) with the number of parameters. An iterative method using single-coordinate descent to approximate the inversion grows order (n). The inversion does not need to run to convergence, since any gains at high iteration number will be low compared to noise amplification. The system calibration method is demonstrated with measured pulser data as well as with two LSO-PSAPD detectors in electronic coincidence. After applying the algorithm, the 511 keV photopeak paired coincidence time resolution from the LSO-PSAPD detectors under study improved by 57%, from the raw value of 16.3 ± 0.07 ns full-width at half-maximum (FWHM) to 6.92 ± 0.02 ns FWHM (11.52 ± 0.05 ns to 4.89 ± 0.02 ns for unpaired photons).
We are developing a positron emission tomography (PET) system for breast cancer imaging , . The system uses a two-panel geometry with each 16×9×2 cm3 panel requiring 2304 Position Sensitive Avalanche Photodiodes (PSAPDs) coupled to 8 × 8 arrays of 1 × 1 × 1 mm3 LSO (Lutetium Oxyorthosilicate) scintillation crystals for a total of 294 912 crystal elements. The use of an “edge on” detector design with the photodetector planes oriented orthogonal to the panel planes enables > 90% scintillation light collection efficiency for 1 mm resolution crystal elements, independent of interaction location, directly measured photon depth of interaction (DoI), and 3-D positioning of individual photon interactions . Simulations show that at 5 mCi whole body image dose, the system will have an event rate of 4.6 million events per second . The system will acquire list-mode data.
Each panel has 9216 readout channels. The channel count is read out using 288 RENA-3 readout ASICs, developed by NOVA R&D . The RENA-3 ASIC was chosen for its processing capabilities and density of channels—36 channels, each consisting of preamplifier, Gaussian shaper, trigger, and timestamp circuitry—and for its low power consumption of less than 6 mW per channel. Each RENA-3 chip can handle a maximum event rate of 83 kHz, though the average event rate per RENA-3 will be 16 kHz .
This RENA family has been used for gamma ray detection with ,  built using the RENA-3,  using the RENA-2 and , the RENA-1. The RENA-3 is being designed with CZT for use in small animal PET in .
In , the coincidence timing resolution of two LSO-PSAPD detectors was shown to be 15.5 ns full-width at half-maximum (FWHM), similar to the time resolution reported for the RENA-3 in  and . The timing resolution for other PMT-based breast dedicated PET imaging systems is under 10 ns FWHM: 1.2 ns , 5.2 ns , and 8.1 ns . In order to achieve sub 10 ns timing resolution with PSAPDs, additional calibration parameters were used and a calibration method was created in order to compensate for the nonoptimal timing of the RENA-3 ASIC.
We show that per crystal coincidence time calibration from list-mode data can be formulated using a maximum-likelihood framework and solved with a direct or iterative convex optimization method. We also demonstrate that additional correlated noise can be accounted for to further improve the annihilation photon pair coincidence time resolution.
The RENA-3 ASIC consists of 36 configurable channels . Each channel can be configured to accept positive or negative pulses, to act as a trigger or be read out on another channel’s trigger. On a trigger, the RENA-3 system records a timestamp for the triggered channel and amplitude data from user specified channels. The timing is obtained through two mechanisms: a leading edge discriminator and a quadrature time stamp. This section explains these two mechanisms.
The RENA-3 uses a leading edge discriminator as the trigger. This consists of a comparator with a programmable threshold value. For our application, timing is done with a negative pulse signal, so when the input signal passes below the threshold, the RENA-3 senses a trigger.
The time of the trigger event is dependent on the input signal shape. Prior to reaching the comparator, the energy pulse is shaped by a preamplifier, followed by an optional pole-zero cancellation network, differentiator, and fast shaper. The resulting pulse shape present at the comparator is seen by changing the threshold on one channel while holding another constant. Fig. 1 shows the shape of the pulse that is seen by the comparator. For this example, a square wave signal, with a rise time of 5 ns, is AC-coupled to two channels on one RENA-3 chip. The comparator threshold on one channel is held constant and the arrival time reported by this channel is a reference for the signal arrival time. The threshold voltage on the second channel is varied and the difference in its arrival time to that of the reference is plotted.
For this demonstration, the differentiator sets the dc level of the shaped pulse to 3.5 V. The negative pulse shape is visible and the last threshold to record triggers was 3.17 V, indicating that for this square wave, the shaped pulses have amplitude of around 0.33 V. The middle range has a slope of 2.6 V/µs. The 10%–90% rise time is in the 80 ns range.
Many of these factors impact the timing capabilities of the RENA-3 and much of this paper will focus on explaining how to compensate for time-walk introduced by the use of leading edge discrimination.
When the leading edge discriminator comparator threshold is passed, two different methods capture the time: a coarse digital timestamp, and a fast analog timestamp. The coarse timestamp is recorded by a counter on an FPGA when the RENA-3 reports it has seen an event. For the fast analog timestamp, when the comparator threshold is passed, the voltage of the two quadrature signals are sampled and later converted to digital values. The relationship between the signals encodes the arrival time.
The quadrature signals, U and V, are sinusoids with one lagging the other by 90° in phase. In Fig. 2, a plot of U and V signals recorded for several events shows the resulting timing circle for one channel. To find the time difference between events that occurred on two different channels, the location of each event on the timing circle is found, and the angle difference between the two measurements directly correlates to the time difference. The angle difference divided by the frequency of the signals gives the time difference. In our setup, the frequency of the signals is 1 MHz.
As seen in Fig. 2, dc offsets in the ADC values need to be taken into account before calculating the angle for timing. The skew in the circle is due to the phase difference in the two signals not being exactly 90° or the sampling occurring at slightly different times. More eccentricity can occur when the magnitude of the two quadrature signals are not identical.
Quadrature timing does not have uniform time bins, due to the phase being calculated by a uniform sampling of amplitude of the quadrature signals, rather than a direct phase measurement. With the quadrature signals maximum range of 4876 counts, the largest time bin is under 100 ps, occurring at phases of pi/8 + n * pi/4. This is a differential nonlinearity of 50 ps. Uncorrected, a 3° phase error between the signals can lead to an integral nonlinearity of 4 ns. Fig. 2 has a phase error of about 3°.
The fast timestamp is determined by the angular separation of the two events on the circle. Imperfect signals or imperfect sampling causes distortion in the circle. The radius of the circle averages 2438 counts and the standard deviation of the radius at each angle is nine counts.
The timing pickoff can be offset or distorted in a number of different ways. Some distortion is directly due to variations in the chip, such as individual channel delays and noise from the quadrature signals. Other distortion is related to the actual measurement, such as crystal location on the PSAPD and variance in amplitude of the output signal.
With identical noiseless input pulses on two channels, it is possible for the channels to report slightly different arrival times. The largest impact is likely the difference in comparator switching delays. Each comparator will have a different delay when its threshold is passed, reporting an arrival time offset by the delay amount.
Other individual channel delay variation can occur in the preamp, zero-pole cancellation network, differentiator and fast shaper that process the pulse before the comparator. These delays will vary from channel to channel and chip to chip and will also create a fixed offset to the reported photon arrival time.
Finally, any variation between channels in the threshold voltage or dc bias voltage will cause offset in measured photon arrival times. A threshold difference could be caused by variations in digital to analog converters or in reference levels. The dc bias voltages could be different due to I-R drop across the chip. If the pulse shape at the comparator is linear, these variations will affect the arrival time measurement by a fixed offset. Other pulse shapes will have varying effects.
The comparator threshold has a range of 3.75 V, with 8-bit resolution, for about 15 mV steps. The experiment with the square wave signal, described in Section II-A, showed that a one-bit difference in threshold causes an offset in reported arrival time of 5.8 ns. This is a significant difference for our application, so a hardware threshold calibration to handle the cumulative fixed offsets is insufficient to calibrate the individual channels. A post-acquisition calibration of a constant offset is necessary.
A charge injection pulse was sent to the 36 channels on the RENA-3 ASIC. The pulser data set consisted of 54 098 pulses. The same threshold DAC value was used on each channel. Histograms of the reported arrival time relative to channel 16 shows several distributions with a range greater than 25 ns. The different arrival time offsets are visible in Fig. 3. From this figure, it can be seen that an iterative method of averaging time differences should be feasible.
A second timing effect is also observed using a pulser. Pulsing two channels and plotting the U and V signals against the timing difference shows a nonconstant, nonlinear correlation between the quadrature signals and time difference. The observed correlation is shown in Fig. 4. This correlation is due to the distortion in the analog timing circle, shown in Fig. 2.
The PSAPD produces further time dispersion effects. The first is due to the variation in RC time constant across the device. The PSAPD contains a resistive sheet to determine positional information , so charge created by light that enters near a corner of the device will see an effective resistance significantly lower than that of charge created in the center of the device. Therefore, annihilation photons detected by crystals at the edge of the device will appear to arrive much sooner than those interacting in a crystal near the center, as seen in Fig. 5. This effect is also mentioned in , though no calibration procedure was proposed or reported. For our system, the time delays for the 8×8 crystal arrays on PSAPDs vary up to 20 ns . Fig. 6 shows a map of the time delay for a crystal array found in the calibration, demonstrating the time delay versus the relative position of the crystals. A constant offset in timing per crystal compensates for this effect.
The second effect on coincidence timing resolution from the detector is due to energy resolution. The energy resolution of this system is 14.4 ± 0.8% FWHM . The different amplitudes of pulses cause the leading edge discriminator to cross its threshold value at different delays. A higher amplitude signal will appear to arrive earlier than a lower amplitude signal. With the RENA-3’s fast shaper 10%–90% rise time of about 80 ns, a pulse with energy 7.2% lower would cause the system to trigger about 7.5 ns later. With the imperfect energy resolution of this system, the time walk inherent to a leading edge discriminator is significant, even when gated for the 511 keV photo peak. The delay will be even more pronounced with if lower energies are used. The delay is linearly related to pulse height if the rising edge of the pulse is approximately linear.
Many methods have been created for the timing calibration of PET scanners. Methods range from using time alignment probes and calibrating manually to using algorithms and timing histograms to calibrate offsets for each block detector –, and then for each individual crystal . Some algorithms use direct calculation of the calibration offsets , while others use iterative methods , , . Here we show that timing calibration, as a function of any number of calibration parameters, is a convex optimization problem. This allows the solution to be found using conventional methods, such as direct inversion or coordinate descent, an iterative method. The calibration parameters in this system include delay for each individual crystal and correlated noise from the quadrature timing.
We consider a list of n coincidence events (or annihilation photon pairs) detected by the system. The k,th event comprises the following information: the indices of the crystals involved ik and jk, the two arrival times and , and some additional data, which will be denoted by two column vectors and . The column vectors contain the amplitude of both pulses, and , and the quadrature time samples , and .
The timestamps measured on crystal i follows a random distribution, described by the random variable T
where τ is the true photon arrival time, δ is the per crystal delay and ε is the measurement noise. To simplify later calibration, the per crystal delay is modeled as a linear combination of the event data, including a constant bias
where ci; is the vector of calibration parameters for each crystal element.
Since the breast imaging panels are 4–8 cm apart, the noise bias effects from time of flight are negligible. The measurement noise is modeled by a zero-mean normal distribution with variance σ2.
The time difference measured between both photons in a coincidence event follows a random distribution given by
Since the time of flight, , is negligible, it is ignored. As a result, the vector of measured time difference is well described by an uncorrelated Gaussian random vector
The per-crystal delay differences are a linear combination of the per-crystal calibration factors c
where is formed using the per-event data and .
Equation (6) shows an example matrix formulation of the time calibration problem for five annihilation photon pairs in a system with four crystal elements or arrays
The goal of the time calibration is to estimate the calibration parameters per crystal element given the measured data. Maximum-likelihood (ML) estimation was used for this purpose. The log-likelihood function is given by
where is the measured time difference.
The ML calibration parameters c are such that is minimized, which is more commonly known as a “least-squares problem.” It should be noted that this approach also optimizes the coincidence time resolution of the system since
is simply the calibrated time difference.
A number of other objective functions, such as the L1 norm or a weighted sum are available to evaluate the timing error. With an initial system calibration scan with a weak source and a low randoms rate, the L2 method is sufficient. However, to adjust the calibration using data from a real scan with a high randoms rate, an alternative cost function, with a smaller penalty for large errors, could be used.
Formulating a maximum-likelihood problem amounts to solving the following least squares problem: select the appropriate vector c to minimize . This can be directly solved with the pseudo inverse, c = (A T A)−1 A T Δt.
When A becomes large, the inverse of A T A is computationally intensive, and an iterative single-coordinate descent algorithm is used. For each iteration, the parameters for a single crystal or channel are optimized while the parameters for other crystals or channels are held constant. This procedure is repeated until no further improvements in coincidence time resolution can be achieved.
The authors in  show a similar problem formation, and are, in effect, using a single-coordinate descent method, itera-tively optimizing a parameter for each coordinate, the per crystal delay time.
To calibrate the individual RENA-3 channel delays, the ASIC inter-channel delay is calibrated using a single fixed delay per channel, each an element in the vector c. One channel is used as a reference and does not have a calibration offset, so the vector c has one element fewer than the number of channels. Each row of the matrix A contains one coincidence event. The kth row represents the kth event, recorded between the channels i and j. On this row, all the coefficients are zero except the ith and jth columns, which are equal to (+1) and (−1), respectively. The size of the matrix A is the number of coincident event rows times the number of channels less the reference channel.
The quadrature timing noise is calibrated by adding information to the matrix A and vector c. A polynomial correction for the two signals, u and v, is achieved by adding correction coefficients for each of the channels into the vector c. The number of rows in the matrix A is increased and for each coincident event the recorded signal values are placed into the appropriate column. The known values, u, v. u2, v 2, and uv are included in the matrix A, and five additional calibration parameters for each channel are added to c. The constant term is included with a constant term from the delay corrections. For 10 channels, there are 59 calibration parameters.
Knowing the channel to channel variance for all combinations, the individual single timing resolution per channel can be estimated. Assuming a Gaussian noise distribution on the timing error and independence between channel measurement errors, the combined variance for two channels is the sum of the individual variance between the channels. The combined variance is directly related to the square of the FWHM of the pulser data coincident time resolution.
Two 8×8 arrays of 1 × 1 × 1 mm3 crystal elements, each coupled to a 10 × 10 mm2 PSAPD, is connected to one RENA-3 chip. The signal of the PSAPD is higher than the dynamic range of the RENA-3 chip, so attenuation is required. The negative, common, signal is split through a capacitive divider between two channels, the larger of which is saturated and used for timing information, while the smaller is used for energy information. The positive signals, which are split by the PSAPD to provide positional information, are attenuated to ground using a capacitive divider . A 22Na point source is placed between the crystal arrays, which are aligned edge-on with respect to the source (Fig. 7).
The timing calibration optimization is performed using the coordinate-descent method, treating the parameters for each 1×1×1 mm3 crystal element as a different coordinate. For each test, the calibration is performed using a subset of collected data; the coincidence time resolution is reported using the remaining data.
Both the RENA-3 and the LSO-PSAPD detectors have influence on the calibration procedure and calibration factors. In order to separate out the effects, different experiments were run. In order to determine the effect from the RENA-3, a pulser was used to generate timing data. Once the effects from the RENA-3 were known, the pulser was replaced with the LSO-PSAPD detectors. All pulser data was calibrated using the direct convex optimization method. The PSAPD detectors were calibrated using the iterative convex optimization method.
Fig. 8 shows the results from a correction for ASIC per channel delay. With the delay correction, pulser coincident time resolution from all channels to channel 16 is 4.39±0.04 ns. This is a significant improvement compared to the 25 ns of channel skew as seen in Fig. 3.
In a second calibration of the data, randoms were artificially added to the data set by arbitrarily pairing list mode data. The distribution of the resulting randoms does not change after application of the algorithm. The algorithm improves timing between coincident photons and has no nonlinear effects on randoms data.
The histograms in Fig. 9 show the results after calibration for quadrature timing noise, using the polynomial correction previously described. As with the ASIC inter-channel delay correction, the algorithm has no impact on the randoms distribution. With correction for delay and quadrature timing, pulser coincident time resolution for all channels to channel 16 is 3.37 ± 0.02 ns.
Fig. 10 shows the pulser coincidence time resolution for all pairings of 35 channels of the RENA-3. The coincidence time resolution between any two channels varies from 1.95 ± 0.01 ns at best to 5.87 ± 0.02 ns at worst. The estimated unpaired time resolution per channel is shown in Fig. 11. The best individual timing resolution is 1.38 ± 0.02 ns FWHM; the worst is 4.48 ± 0.02 ns FWHM. The estimated paired time resolutions are shown in Fig. 12, which compares well to the corresponding measured values shown in Fig. 10. The mean difference for all pairings is 2.34 ± 1.5%.
For the LSO-PSAPD detector setup, approximately 450 000 coincident events were collected. Around 50 000 coincidence events remained when both arrays were gated to ±15% around the 511 keV photopeak, shown in Fig. 13. The crystals locations were segmented manually from the acquired crystal image as described in . The flood histograms of the crystal arrays are shown in Fig. 7. The pin cushion effect in the flood histograms is a characteristic of PSAPDs due to nonlinear equivalent resistances across the sheet and is also shown in  and .
The precalibration paired photon coincidence time resolution was found to be 16.37±0.07 ns FWHM. After correction for per crystal location on the PSAPD effects, the paired coincidence time resolution improved to 8.42 ± 0.03 ns FWHM. The additional signal amplitude and quadrature timing correction improved coincidence time resolution to 6.92 ± 0.02 ns FWHM over all crystal elements in the two arrays. Timing histograms are shown in Fig. 14.
The calibration delays resulting from the per crystal PSAPD delay calibration are shown in Fig. 15. The change in RC time constant due to the capacitance of the readout system and different locations on the resistive sheet of PSAPD across the array is visible. The corner crystals have significantly less delay than center crystals and the same pattern is visible in both crystal arrays. There is also an apparent bias between the two crystal arrays. This is due the inherent delay between the two PSAPDs’ readout channels.
A closer examination, in Fig. 16, of the optimized time spectrum (6.92 ± 0.02 ns FWHM) over the two 8×8 arrays shows two side peaks outside of the center peak. The flood histograms from each of these side lobes for both PSAPD detectors are shown in Fig. 17.
Fig. 18 shows a histogram of the lines of response from the coincidence data backprojected through the space between the two detectors. The data were normalized for varying detector element efficiencies by weighting for the sensitivity of each line of response. The resulting superposition of backprojected lines converge at a point in the middle corresponding to the source position, confirming that most of the coincidence data were trues. The FWHM of a vertical profile from the histogram taken through the point source is 0.90 ± 0.05 mm.
Calibration on a data set with a high randoms rate would significantly slow the convergence of a single-coordinate descent optimization. The attempt to minimize the timing difference of randoms would counteract the minimization of the timing difference of trues. In this case and also in cases of low randoms rates, only data with near zero time difference should be optimized for each coordinate optimization.
For each crystal, or coordinate, the time difference to all others is calculated. This data could be gated to constant window or an adaptive window, so that the optimization could be calculated data with a higher percentage of trues. This could allow for patient data to be used for calibration.
For large quantities of data, only a subset of the data interacting on each crystal, or coordinate, is necessary to be used in the coordinate optimization. This would allow for less calculation per iteration.
As seen in Fig. 10, even with calibration, a significant variation in timing FWHM exists between different pairings of channels in the RENA-3 ASIC. For a global system timing window, to ensure trues are not incorrectly removed from noisy lines of response, the timing window would need to be based on the worst line of response. If a smaller timing window is used, the sensitivity weighting of the lines of response would be incorrect.
Given the proximity of the scanner to the body and heart, as well as the low uptake in the breast, a significant percentage of the events will come from out of the field of view. Even when gated for energy and timing, a high percentage of coincidence events will be randoms. Reducing the timing window in lines of response can significantly reduce randoms. A reduction in timing window of 1.4 ns will decrease the randoms on that line of response by 10%.
In order to remove more randoms, an optimum time window could be created for each coincident pair based on the channels it came from. The timing resolution for each channel would be calculated as in Fig. 11 and stored in a lookup table. When a suspected coincidence pair is analyzed, the timing resolution from the individual elements would be summed in quadrature to create the timing window for the line of response. This way randoms can be removed from the data set for better time resolution pairings, while true photon coincidences from poor time resolution pairings are not removed from the data set.
As with individual channels pairs, coincidence time resolutions can be found between individual crystal element pairs. This was not performed with this data set since the point source did not provide lines of response between all individual 1×1×1 mm3 crystal element pairs. This will be done in future work using a larger area source.
The side lobes seen in Fig. 16 are due to multiple interactions inside the array, when a single photon scatters and is absorbed into another 1×1×1 mm crystal element in the same array. The common and spatial signals seen by the RENA-3 system is the sum of the signals of the two interactions.
Since the scatter interactions occur in two different portions of the array, each is subjected to a different RC time delay. The 290 ns shaper of the RENA-3 chip effectively removes the time delay between the signals and the peak detect reads approximately the sum of the two maximums of the signals, yielding a summed energy result appearing in the 511 keV photopeak. The spatial signals will result as a weighted mean of the two interaction locations, with Anger logic resulting in a location as some combination of the interaction locations.
For timing, no shaper is used, and the RC delay of the PSAPD for the two interaction locations is significant. The leading edge discriminator determines timing from the faster event, the one with the lowest RC time delay. The lowest RC time delay is the one nearer to a corner, so this interaction will have its timing recorded.
Since the reported spatial position is between the two interactions and the timing will be from the corner, an inaccurate (too large) time shift correction will be applied. Thus, for multiple interaction events, the timing will always be overcorrected, creating the side peaks. One side peak is from a scatter interaction in one array and the other comes from a scatter interaction in the opposing array.
This is confirmed by looking at the histogram of the acquired event locations for each lobe in each array. From Fig. 17, it is apparent that for the negative lobe, the left array has visible crystals, while the right array does not, as would be expected from an average of multiple interactions. The opposite is true for the positive lobe.
A fraction of such multiple interaction events with interaction locations well-separated on the PSAPD will be removed from the data set with an appropriate timing window, since the delay across the PSAPD is on the order of 20 ns.
A number of issues can affect coincidence timing resolution when detecting LSO-PSAPD scintillation detector signals with the RENA-3 ASIC. However, many of these issues are strongly correlated with sources of timing dispersion and can be calibrated out of the system. The system calibration method used was formulated as a convex optimization problem, and by finding the solution, we are able to convert a leading edge discriminator-based ASIC with raw timing variation of tens of nanoseconds into a fast timing data acquisition system appropriate for PET. Correcting data for ASIC time channel delay, PSAPD RC time constant delay, leading edge discriminator time walk, and quadrature timing correlated noise, a paired coincidence timing of less than 7 ns is achieved. Also, with PSAPDs and per crystal time correction, the events with multiple interactions will be removed from the data set.
This work was supported by the National Institutes of Health (nih) under Grant r01 ca119056, Grant r01 ca119056-04s1 (arra), Grant r33 eb003283, Grant r01 ca120474. Asterisk indicates corresponding author.
Paul D. Reynolds, Department of Electrical Engineering, Stanford University, Stanford, CA 94305 USA (Email: paulr2/at/stanford.edu).
Peter D. Olcott, Bio-Engineering Department, Stanford University, Stanford, CA 94305 USA.
Guillem Pratx, Department of Electrical Engineering, Stanford University, Stanford, CA 94305 USA.
Frances W. Y. Lau, Department of Electrical Engineering, Stanford University, Stanford, CA 94305 USA.
Craig S. Levin, Departments of Radiology, Physics, and Electrical Engineering, Stanford University, Stanford, CA 94305 USA (Email: cslevin/at/stanford.edu).