Novel crystal timing calibration method based on total variation

A novel crystal timing calibration method based on total variation (TV), abbreviated as ‘TV merge’, has been developed for a high-resolution positron emission tomography (PET) system. The proposed method was developed for a system with a large number of crystals, it can provide timing calibration at the crystal level. In the proposed method, the timing calibration process was formulated as a linear problem. To robustly optimize the timing resolution, a TV constraint was added to the linear equation. Moreover, to solve the computer memory problem associated with the calculation of the timing calibration factors for systems with a large number of crystals, the merge component was used for obtaining the crystal level timing calibration values. Compared with other conventional methods, the data measured from a standard cylindrical phantom filled with a radioisotope solution was sufficient for performing a high-precision crystal-level timing calibration. In this paper, both simulation and experimental studies were performed to demonstrate the effectiveness and robustness of the TV merge method. We compare the timing resolutions of a 22Na point source, which was located in the field of view (FOV) of the brain PET system, with various calibration techniques. After implementing the TV merge method, the timing resolution improved from 3.34 ns at full width at half maximum (FWHM) to 2.31 ns FWHM.

Keywords: PET, timing calibration, total variation constraint, merge component (Some figures may appear in colour only in the online journal)

Introduction
Timing resolution is one of the most important performance parameters of positron emission tomography (PET) systems. A good timing resolution implies that a narrow coincidence window can be used for detecting the coincident pairs of annihilation photons, reducing the impact of random event counts. In addition, a good timing resolution implies a more precise measurement of the apparent time of the photons annihilation, which yields more precise estimation of the time-of-flight (TOF) (Kadrmas et al 2009) in PET scanning, significantly affecting the quality of the reconstructed images (Conti 2011). Therefore, a PET system should have a good timing resolution. However, the PET system's timing resolution is easily affected by the timing delays, which are usually caused by variations in the scintillation responses, by the timing of pick-off circuits and by circuit delays across different channels. These factors tend to negatively affect the overall timing resolution of a PET system. Therefore, a high-accuracy and crystal-level timing calibration method is required to compensate for these delays. In the crystal-level timing calibration, every crystal has a unique timing calibration; by contrast, in the detector-level calibration, every detector block has a unique calibration, and all crystals belonging to the same detector block have the same timing calibration. Since the time delays are different from crystals to crystals and the differences are unknown, a high-accuracy crystal-level timing calibration is a challenge issue.
A variety of methods have been developed for the timing calibration of PET scanners, including those that use a dedicated fast photo-multiplier tube (PMT) combined with a radioactive source as a reference signal for measuring the time delays for each measured event (Thompson et al 2005, Hancock andThompson 2010). Some researchers have proposed using specially designed sources, such as a rotating line source (Thompson et al 2005), a point source with a scatter cylinder (Perkins et al 2005), and an annihilation target (Li et al 2016), for timing calibration. For all varieties of the specially designed sources, the true location of positron emission can be calculated or determined before or during timing calibration. Hence, the difference between the true temporal difference and the measured value can be used as the timing calibration value. However, the necessity of a specially designed source limits the applicability of these methods. To solve this problem, some researchers proposed using the scan data of patients (or phantoms) along with the corresponding reconstructed images, to calculate the timing calibration factors (Werner and Karp 2013). Others have proposed to scan a centred cylindrical phantom and to use an iterative method for estimating the timing calibration values (Park et al 2008). Because these approaches require no additional detection elements or special radioactivity sources, it is easier to obtain data. However, in the latter approach results are easily affected by noise in the data (or require preprocessing of the data to obtain good reconstructed images). Therefore, a novel, high-accuracy and robust timing calibration method is required (Reynolds et al 2011).
We have already developed a segmented-detector-unit timing calibration method using an iterative shrinkage threshold (IST) with an L1 norm constraint in the governing linear equation (Yu et al 2015). The method provides a sub-detector-based timing calibration procedure for a PET system. Each segmented detector unit has a dedicated timing calibration value and all crystals in the same detector share the same timing calibration value. Compared with the detector-level timing calibration, this method offers a more robust and accurate timing calibration. However, because using a large number of crystals is computationally demanding, the IST method does not yield a true crystal-level timing calibration, and is quite time-consuming when used for calibrating a PET system with a large number of crystals. Therefore, developing a fast, crystal-level and high-accuracy timing calibration method remains a significant challenge.
Toward solving this problem, we have proposed a novel crystal-level timing calibration method based on a total variation (TV) constraint, henceforth called 'the TV merge method'. The proposed method formulates the timing calibration problem as a linear problem and then introduces a TV constraint to guarantee a highly accurate estimation. In addition, a merge component is used during the calculation process of the timing calibration factors, to achieve timing calibration at the crystal level and to make the computation less demanding. Both simulation and experimental data were used for demonstrating the effectiveness and robustness of the TV merge method. As a reference method, we used the least squares method (Park et al 2008), which is widely used as a software timing calibration method.

Materials and methods
The proposed TV merge method consists of three components: (1) a crystal grouping component, (2) a TV component, and (3) a merge component. Schematic of the TV merge method is shown in figure 1. The crystal grouping and the merge components are closely related, because both are used for solving the memory problem. The TV component is entirely independent. Therefore, we first describe the TV component.

TV component
The basic timing calibration process can be formulated as a linear equation where x is an × n 1 vector of timing calibration factors for compensating for temporal delays, with n being the number of the detector units or the number of crystals. ∆T is the peak location of the time histogram acquired at the channel of concern. When placing a cylinder phantom at the centre of a scanner, the peak location of the time histogram at each line of response (LOR) is zero, and its width reflects the size of the cylinder and the temporal resolution of the system. A non-zero ∆T signifies the mis-calibration of the time calibration factors. In equation (1), ∆T is an × m 1 vector, where m is the number of LORs. A is an × m n system matrix that defines the system's LORs where the () un indexes of the detection unit. For high-accuracy timing calibration estimation, we propose to add a TV constraint to equation (1). The TV-based model of equation (1) can be expressed as an optimization problem: where ( ) x TV represents the TV constraint and is defined as where is the Frobenius (or Euclidean) norm, the sum of the discrete gradient of x in every smallest unit i. The TV constraint is a widely used regularization approach capable of properly maintaining edges and removing noise from data (Rudin et al 1992, Beck and Teboulle 2009a, Strong and Chan 2003.
The values of the time delays of the crystals are different, even though the crystals belong to the same crystal array or detector block. However, the differences of the time delay between different crystals in the same detector block are smaller than those between crystals in different detector blocks. In the vector of timing calibration factors x, all of the timing calibration factors of crystals are arranged according to their detectors. Because the values of the time delays of different detectors are different, groups of timing calibrations factors corresponding to different detectors are clearly distinguishable in this vector. However, noise in the data is likely to weaken the boundaries and reduce the accuracy of the timing calibration estimation. Thus, the purpose of using the TV constraint is to maintain the group boundaries, marking different crystal arrays, and to remove the noise in the dataset.
Solving equation (2) directly is difficult. Therefore, we introduce a new auxiliary variable u, which is defined by = u Dx i i . The main purpose of this process is to transfer D x i out of the norm . The optimization problem in equation (2) is clearly equivalent to Several solvers have been developed for dealing with this TV optimization problem, including the second-order cone program (SOCP) (Alizadeh andGoldfarb 2003, Goldfarb andYin 2005), the iterative shrinkage/thresholding (IST) (Bioucas-Dias andFigueiredo 2007, Beck andTeboulle 2009b) algorithms, and the augmented Lagrangian method (ALM) (Wu andTai 2010, Wu et al 2011). In this work, we used the ALM and we transformed equation (3) into a minimization problem of the augmented Lagrangian function, denoted by ( ) L x u , i . The corresponding augmented Lagrangian problem of the optimization function in equation (3) can be expressed as where the first term in equation (4) is the regularization term, the second and the fourth terms are linear, and the third and the fifth terms are quadratic. The linear and quadratic parts ensure the accuracy and robustness of the TV constraint. Here, θ i , β ∆ i and λ are the multipliers or weights of these four penalty parts. In general, β i and λ are constants, whereas θ i and ∆ are updated in each iteration. According to the classical augmented Lagrangian method, exact minimization is required in each iteration, which is time-consuming and not necessary. We therefore suggest using the alternating direction method (ADM) for solving equation (4). The ADM, which was originally proposed to deal with parabolic and elliptic differential equations (Peaceman and Rachford 1955), was embedded here to solve equation (4) efficiently. More specifically, at the kth iteration we first fixed u at its currently estimated value u k , and then obtained x can be achieved by solving the following sub-problem: When fixing u i for calculating x, the terms containing u i k are constant and are omitted. The gradient of ( ) Q x with respect to x can be easily obtained from equation (5), as follows: , we can obtain the exact minimizer of ( ) Q x directly. However, this calcul ation is too costly to implement numerically, especially when the size of the matrix is large. To obtain the values of + x k 1 more efficiently, the steepest descent method with a proper step length was used iteratively by applying the recurrence formulation where α k is the step length. Each iteration of the steepest descent method demands updating the gradient to update the estimated value of + x k 1 . Therefore, the step length should be chosen carefully to obtain an accurate solution. In this paper, α k was determined by Armijo-like nonmonotone line search (Beck and Teboulle 2009b).
Similar to the method for obtaining + x k 1 , we fix + x k 1 and calculate the solution of + u k 1 in the next step which can be expressed as follows: When the constant β > 0, the minimization in equation (8) can be easily calculated using the 2D-shrinkage formula (Bioucas-Dias and Figueiredo 2007). Therefore, we can obtain This algorithm has several parameters. First, β i and λ are the quadratic penalty parameters. The values of these two parameters determine the accuracy of the estimated values relative to X Yu et al Phys. Med. Biol. 61 (2016) 7833 the true values. In each iteration, these parameters are constant. The values may vary in the range from 2 4 to 2 13 , according to the exact level of noise and required accuracy. A detailed description of the optimal selection of β i and λ is given in Beck and Teboulle (2009b) and Wu and Tai (2010). In this paper, λ was 2 8 , and β i was 2 4 . These values were selected after performing many tests and comparing the results for determining the optimal values. The Lagrangian multipliers θ i and ∆ were initialized to 0 and were updated during each iteration. The exact update formulas of the two Lagrangian multipliers were proposed by Hestenes (1969), as follows: After a number of iterations, the final timing calibration estimations are obtained. The program terminates after a certain number of iterations (300 in this paper), or when the relative stopping criterion (based on empirical estimations) is reached:

Crystal grouping
Both the crystal grouping and merge components were designed to achieve a real crystal level calibration and to solve the memory problem in the crystal level estimation. Schematic of these components is shown in figure 2. In this paper, we first provide an introduction to crystal grouping, which consists of two components, the × N N grouping component and the / × × N N 1 1 grouping component. The former component is designed to obtain the final timing calibration factors that are used for compensating for the temporal shifts or delays in the PET system. The latter component is the special grouping during the process of estimating the timing calibration factors. In the TV merge method, the value of N in both components were the same.
In the × N N components, the crystal array in each detector block was grouped into × N N units. Here, N is a divisor of the number of crystals. For example, the valid values of N are 2, 4, 8, 16, or 32 for a crystal array. However, memory capacity significantly affects on the final value of N. As an example, in the case of a system with 640 detector blocks with × 32 32 crystals each, 16 GB of memory are sufficient only for N = 4 (Yu et al 2015). It is difficult to obtain × × N N n crystal grouping timing calibration factors directly for all detector blocks (n is the number of detector blocks) when the value of N is relatively large. The / × × N N 1 1 crystal grouping and merge components were used to solve this limitation. The main purpose of the / × × N N 1 1 crystal grouping is to reduce the dimension of the system matrix and the timing calibration sequence. In the proposed method, the crystal array in each detector block is first grouped or projected onto × N 1 units by row and × N 1 units by column, respectively. For example, if N is 32, then every × 32 32 crystal array is projected onto × 32 1 units by row and × 1 32 units by column. The mathematical formulation of the proposed timing calibration in equations (1) and (2) is the 'crystal grouping-free for mulation'. It means that all types of crystal groupings can be used for timing calibration based on equations (1) and (2). For example, when the × N 1 crystal grouping is used, the dimension of the timing calibration factors parameter x is × × N n 1 (while for × N N crystal grouping, it becomes × × N N n). Because the ∆T is the peak location of the time histogram X Yu et al Phys. Med. Biol. 61 (2016) 7833 acquired at the channel of concern, when the × N 1 crystal grouping is used, the corre sponding time histogram for calculating ∆T changes from the histogram of all crystals ( × N N crystal grouping) to that of the sum over all crystals in the same row. If every crystal coincides with all other crystals, the size of ∆T is ( ) ( )/ × × × × × − N n N n 1 1 1 2 (theoretically, for all possible LORs, the number of m is ( )/ = − m n n 1 2). Then, the size and the formulation of the system matrix can also be determined. After that, we can use the solution process from equations (3)-(13) to calculate the timing calibration estimations for the × N 1 crystal grouping. However, because we aim to obtain a true crystal-level calibration, and the estimations using the × N 1/ × N 1 crystal grouping only contain the calibration information in row and column directions, we need a method for merging them into the final × N N timing calibration. To accomplish this, we proposed the merge component.

Merge component
The merge component is designed to solve the computer memory problem when the number of the crystal groupings becomes large. The basic idea behind using the merge component is that we use the row and column timing calibration estimations with the × N 1/ × N 1 crystal grouping to obtain the final × × N N n timing calibration factors, instead of directly calculating the final × × N N n timing calibration factors from the measured data. The merge component is used to reduce the memory and time associated with the calculation of the timing calibration factors. However, it is difficult to provide mathematical justification of the exactness of the proposed merge component method. Therefore, we designed some experiments in the Experiments and Results section, and used the empirical results obtained in these experiments to justify the proposed merge component.
The merge method assumes that the × × N n 1 estimation contains the exact timing calibration information (values and distributions) of the crystal array its rows, while the × × N n 1 estimation contains this information in its columns. These two estimations can be merged into the final × × N N n crystal grouping of timing calibration factors using merge component. We emphasize that both × × N n 1 and × × N n 1 estimations can be directly calculated by using the TV algorithm in section 2.1. For this, the number of the detector units should be set  . The process of the back-projection can be expressed as  ). The final × × N N n timing calibration factors are obtained by averaging the maps A and B.
The main purpose using the merge component is to solve the computation and memory problems caused by an enormous number of crystals. However, the assumption that the × × N n 1 and × × N n 1 estimations contain all of the exact timing calibration information is not perfect. Therefore, the merge component introduces additional estimation errors in the estimation of the final timing calibration factors. However, we assume that it only insignificantly affects the final timing calibration. This assumption is validated in the following section.

Experimental setup
2.4.1. Description of the system. To evaluate the effectiveness of the proposed TV merge method, a high-resolution brain PET scanner (Hamamatsu HITS-655000) was used. The brain PET scanner has been described in detail in previous papers (Yu et al 2015), thus, below we provide only a brief summary of it. The scanner consisted of 32 detector modules, which were arranged in a 430 mm-diameter ring. The transaxial field of view (FOV) was 330 mm, whereas the axial FOV was 202 mm. The detector module consisted of four detector layers. Each detector layer had five detector units, which were lined up axially, and front-end circuits including ASICs for multi-pixel photon counters (MPPCs: Hamamatsu S10931-050P) and signal processing circuits. Each detector unit consisted of a × 32 32 LYSO scintillator array, with a 1.2 mm pitch, and an 8 × 8 array of MPPCs. The intrinsic timing resolutions of the detector module in coincidence with a BaF 2 crystal coupled to a fast PMT for detector units from the first to the forth layers, were 850 ps, 862 ps, 983 ps, and 930 ps FWHM. The overall number of detector rings in the system was 168, and the overall number of crystals was 655 360. Individual events detected by individual crystals were sent to data processing servers and were recorded as single-event data pattern. Coincidence detections were performed by software using the single-event data in the data processing servers.

Real scanning setup.
The real scanning experimental data used for estimating the timing calibration factors were measured using a cylindrical phantom with the diameter of 140 mm and height of 205 mm. The phantom was filled with an 18 F-FDG solution and was located at the centre of the scanner's FOV. The initial activity was 282.3 MBq, and the overall scanning time was 10 h. A Na 22 point source (473 kBq) with the diameter of 0.25 mm, located at the centre of the FOV, was also measured, for calculating the timing resolution, to evaluate the timing calibration performance of different methods. The scanning time was 5 min. Schematics of the measurement setup and data processing are shown in figure 3.

Experiments and results
In this section, we demonstrate the effectiveness of the TV merge method by employing both simulated and real scanning data. We use quantitative analysis for the simulation, because the true values are known. Regarding the real scanning data, the timing resolutions for different methods are provided, to demonstrate the effectiveness of the TV merge method.

Simulation studies
The accuracy of the TV merge method depends on two factors. The first factor is the improvement in the estimation accuracy using the TV constraint, while the second factor is the influence incurred by the merge component. Therefore, one simulation and two comparisons were performed to evaluate the effectiveness of the proposed method using data created in MATLAB, before performing real scanning experiments. To evaluate the influence incurred by the merge component, we need a method that is exactly the same as the TV merge method but is not feature the merge component. Therefore, we developed such a method and named it 'the TV method'. The mathematical procedure underlying the TV method was introduced in section 2.1. Because the TV method does not contain the merge component, the × N N crystal grouping is directly used in this method and it cannot be used for obtaining a crystal-level calibration. The simulation was used to compare the estimation accuracy of the least squares method, the TV method, and the TV merge method. The least squares and the TV merge methods were compared for evaluating the improvement introduced by the TV constraint. The TV and the TV merge methods were compared for evaluating the influence incurred by the merge method. In the least squares method, the linear equation can be formulated as a maximumlikelihood problem − ∆ Ax T min 2 2 . This can be directly solved using the pseudo inverse, . When A is large, an iterative single-coordinate descent algorithm is used (Barzilai and Borwein 1988).
In the simulation, the data y are simulated based on the linear model in equation (1) of the PET timing calibration process, and can be expressed as where A is the system matrix. The overall number of the simulated detector units was 40, and each of the individual detector units coincided with all others. Therefore, the size of the system matrix was 40 columns by 780 (=40(40 − 1)/2) rows. The vector x true contained the true timing calibration factors. The value of x true was selected randomly, with a uniform distribution between −12.0 ns to 12.0 ns (the minimal interval was 0.01 ns). The quantity ε represented the 1% Gaussian noise. To evaluate the accuracy of the estimation, the estimation errors (E r ) of the least squares method, the TV method, and the TV merge method, are where x true contains the true values of timing calibration factors and x est is the estimation of x true . For the least squares methods, E r was 0.1023. For the TV merge method, E r was 0.0848, corresponding to an improvement of 17% over the value for the least squares method. Thus, the TV merge method yielded a more accurate estimation compared with that of the least squares method. For the TV method, E r was 0.0846, and the difference caused by the merge method was 0.24%. It is a relatively small difference and is considered to fall within an acceptable range. Therefore, the TV merge method yielded a more accurate estimation of the timing calibration factors than the least squares method, in terms of the estimation error E r .

Real scanning studies
3.2.1. Comparison between the least squares and the TV merge methods. The conventional least squares method and the TV merge methods were compared in terms of the accuracy of the timing calibration, by evaluating the timing resolutions. Because the system matrix is large, the iterative single-coordinate descent algorithm was used in the least squares method. And the least squares method was only performed in the detector level. For a fair comparison, the TV merge method was also used in the detector level. Therefore, the size of the crystal grouping was 1 for both the least squares and TV merge methods. In this experiment, only the data set of ring unit 2 (located at the centre of the overall system) was used for calculating the timing calibration factors.
The histograms of the point source, compensated for the time delays at each channel, without calibration and using timing calibration factors calculated using these two methods, are shown in figure 4. The timing resolutions were extracted by performing Gaussian fits. The FWHM values were then obtained by evaluating the fit functions. The timing resolutions without calibration, for the least squares and TV merge methods were 3.29 ns, 3.04 ns, and 2.89 ns, respectively, representing an improvement of 12.16% (from the setup without calibration to the one that uses the TV merge method) and 4.93% (from the setup that uses the LS method to the one that uses the TV merge method). The calculation time for the least squares method was 268 s and the calculation time for the TV merge method was 41 s. These results demonstrate that a superior and fast timing calibration can be obtained using the TV merge method. However, compared with the results of the simulated data, the improvement was relatively small. The noise in real data is more complex than that in the simulated data, therefore, the more complex noise reduces the estimation accuracy.

Evaluation of the estimation variance incurred by the merge component.
This experiment was designed to evaluate the influence incurred by the merge component. Although we showed that this contribution is not significant in the simulation studies, real data are usually more complex than simulated data. It is important to check whether the merge component would incur a serious estimation error when used on a real scanning data set.
Because true calibration factors are not known in a real system, we have to evaluate the influence incurred by the merge component by comparing the results of the TV merge method and a method without the merge component. Therefore, the TV method, which was mentioned above, was used in this experiment. The mathematical foundation of this TV method was introduced in the section that described the TV component. Because the merge component is not used in the TV method and the computer memory problem is un-solved, the TV method can only be used for calculating the timing calibration factors when the crystal grouping is small. Because the maximal allowable value of N was 4, which was limited by the computer memory, only the × 2 2 and × 4 4 crystal groupings were used in this study for determining the magnitude of the potential estimation error incurred by merge component.
The timing resolutions of the TV and TV merge methods in the × 2 2 and × 4 4 crystal groupings are listed in table 1. The timing resolutions of the TV and TV merge methods were 2.367 ns and 2.378 ns, respectively, for the × 2 2 crystal grouping, and 2.321 ns and 2.350 ns for the × 4 4 crystal grouping. The difference of the TV merge to TV in the two groupings were 0.011 ns ( × 2 2) and 0.029 ns ( × 4 4), and the Var ratios ( / ( ) = Var ratio variance time resolution TV ) were 0.46% and 1.16%, respectively. Therefore, the estimation result, by contrast, is almost not affected by the merge component.

Comparison between the TV merge method across different crystal groupings.
The effectiveness of the timing calibration at different crystal levels was also evaluated. Different crystal groupings were used to determine whether a smaller crystal grouping results in higher accuracy of the timing calibration, and to delineate the optimal crystal grouping for the timing calibration. The calibrated timing resolutions for different crystal groupings are listed in table 2.
The trend of the timing resolutions for different crystal groupings is shown in figure 5. In principle, smaller crystal groupings improve the accuracy of the timing calibration. However, the improvement in the timing resolution from × 16 16 to × 32 32 grouping was not significant. The difference between the × 16 16 and × 32 32 crystal groupings was marginal (only 0.006 ns, within the range of the statistical error).

The improvement introduced by the TV merge method.
The effectiveness of the TV merge method was evaluated. The timing resolutions for the detector-level least squares calibration and the crystal grouping TV merge method were calculated and are shown in figure 6. The timing resolutions were 3.04 ns and 2.312 ns, respectively. The relative improvement of the TV merge method over the least squares method was ~0.7 ns (~22.9%). Thus, the TV merge method significantly improved the accuracy of the timing calibration. The reconstructed images of the cylinder phantom with detector-level least squares calibration and crystal grouping timing calibration are shown in figure 7. The used coincidence time window were 7 ns and 5 ns for least squares and TV merge, respectively (twice of the timing resolution). The  reconstruction algorithm was list-mode DRAMA (Tanaka and Kudo 2010). To evaluate the improvement in the reconstructed images introduced by using the TV merge calibration, the standard deviation (Stdev) and the ratio of the Stdev to mean (Stdev/m) of the reconstructed images were used for quantification. The Stdev and Stdev/m of the image reconstructed using the least squares method were 142.8 and 0.105, respectively, while the Stdev and Stdev/m of the image reconstructed using the TV merge method were 129.2 and 0.097. The corresponding improvements were 9.5% and 7.6%. Therefore, the TV merge method significantly improves the quality of the reconstructed images.

Discussion
In light of these results, we conclude that the proposed TV merge method is capable of providing a crystal-level, fast, and high-accuracy timing calibration for PET systems. Compared with all possible crystal groupings, the crystal grouping appears to be optimal for PET systems.
However, more work remains to be done to develop the optimal method for timing calibration. Although timing calibration at the crystal level can yield a better calibration than that performed at the detector unit level, the count problem still prevents making full use of the TV merge method. The count problem arises because each value of the time difference vector is obtained by finding the time difference value at the peak position in the histogram for each LOR. If LORs do not have sufficient counts, the time difference histogram of LORs will be dominated by noise. Therefore, it is hard to accurately determine the exact peak position, which results in a significant variance in the final timing calibration estimation. This problem can be solved by using long scanning times. For example, the overall scanning time for obtaining the measurement data set in this paper was 10 h. However, long scanning times are quite inconvenient in practice. Therefore, an ideal timing calibration method should be able to extract highly accurate timing calibration factors from short-term data. Therefore, in our future work we will focus on developing a timing calibration method for working with short-term scanning data.

Conclusions
A novel crystal timing calibration method based on TV has been developed and evaluated using a high-resolution brain PET system. The method provides crystal-level, fast, highaccuracy timing calibration for a PET system with a large number of crystals. The proposed method utilizes a TV constraint and a merge component. The TV constraint improves the accuracy and robustness of the timing calibration estimation and avoids the potential error in the estimation. The merge component successfully solves the memory problem and produces a timing calibration at the crystal level for a system with a large number of crystals. Therefore, the proposed method provides a fast, true crystal-level, high-accuracy and robust timing calibration for PET systems.