A method for improving the temporal resolution of PET systems based on STFT

Abstract

The present invention discloses a method for improving the time resolution of a PET system based on STFT. The method acquires waveform data by setting a point source at a specific location, then obtains short-time frequency domain amplitude information through short-time Fourier transform (STFT), then divides the training, validation, and test sets into training, validation, and test sets, and trains a residual neural network model composed of a residual module and a fully connected layer to achieve estimation of the TOF time of the PET system. The present invention uses a Monte Carlo simulation experiment to perform STFT transformation on the waveform data to obtain short-time frequency domain amplitude information, then uses the divided training set to train the residual neural network model, and tests the trained model on the validation set. Finally, the test set verifies the accuracy of TOF time estimation and the improvement of the time resolution of the PET system based on the STFT model compared to traditional CFD methods.

Description

Technical Field

[0001] The present invention belongs to the technical field of PET imaging, and in particular relates to a method for improving the temporal resolution of a PET system based on STFT. Background Art

[0002] Positron Emission Tomography (PET) is a nuclear medicine imaging technique that achieves diagnostic purposes by labeling substances essential for the metabolism of living organisms and injecting them into the body for detection. These specific drugs labeled with radionuclides are usually called tracers. Commonly used tracers include 18 F. 11 C. 15 O, etc. During metabolism in organisms, radionuclides have a chance to decay and release positrons. After a short drift, the positrons annihilate with surrounding negative electrons, producing two gamma photon pairs traveling in opposite directions with equal energy. These photon pairs are then detected by a detection system consisting of scintillation crystals and detectors, converted into electrical signals and input into a circuit. After certain signal processing, information on the concentration, location, and timing of radioactive substances in life activities can be obtained.

[0003] Furthermore, if the time difference between the arrival of two gamma photons at the detector can be determined, the exact location of the annihilation site can be calculated from this time difference information. This type of PET detection system is called a time-of-flight (TOF) PET system. However, TOF-PET systems have high time requirements. In practice, gamma photon detection is subject to the following errors: 1. Conversion depth: After entering the crystal, the gamma photon will travel a certain distance before being absorbed; 2. The crystal scintillation process: the rise and decay times of the crystal during the scintillation process; 3. Transmission time: the time it takes for the photon to exit the crystal and reach the photodetector; and 4. The single-photon dispersion time of the photodetector. This information is ultimately contained in the photodetector waveform.

[0004] Currently, the main methods for timing photodetector waveforms are the following two methods mentioned in the paper "Signal processing for picosecond resolution timing measurements" (Genat, JF et al., Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, Vol. 607, No. 2, pp. 387-393, June 11, 2009):

[0005] Leading-edge timing method; this method first sets a certain voltage threshold V th The time t1 and t2 when a pair of waveforms first exceed the threshold are taken as the photon arrival time, and the difference between the two is used to obtain the TOF time.

[0006] Constant ratio timing method; this method sets the threshold V th It is set as a percentage of the maximum value of the waveform. Similarly, the time t1 and t2 when a pair of waveforms first exceed the threshold are used as the photon arrival time, and the difference between the two is used to obtain the TOF time.

[0007] In addition to the constant ratio timing methods mentioned above, there are also the zero-crossing constant ratio timing method and the interpolation constant ratio timing method mentioned in the document "Neural network-featured timing systems for radiation detectors: performance evaluation based on bound analysis" (Ai, P et al. published in the Journal of Instrumentation, Volume 16, Issue 9, Pages 09-19 in September 2021).

[0008] Zero-crossing constant ratio timing method; this method copies each waveform into two groups, multiplies the first group of waveforms by a percentage to scale, multiplies the second group of waveforms by a negative sign and delays it, and then adds the two groups of waveforms to obtain the final waveform. The zero-crossing point of the waveform is regarded as the arrival time of the photon. Finally, a pair of waveforms are respectively subjected to zero-crossing constant ratio timing to obtain t1 and t2, and the difference between the two is obtained to obtain the TOF time.

[0009] Interpolation constant ratio timing method; Since the waveform is a digital signal, the ordinary constant ratio timing method sets the time when the waveform first exceeds the threshold as the photon arrival time, which is not accurate. The interpolation constant ratio timing method performs linear interpolation between the sampling points on both sides of the threshold in the waveform, thereby improving the equivalent sampling rate and obtaining more accurate photon arrival times t1 and t2. The difference between the two is used to obtain the TOF time.

[0010] However, these methods do not utilize the potential information contained in the waveform; therefore, it is necessary to learn the potential features in the waveform through some deep learning methods to further improve the temporal resolution of the PET system. Summary of the Invention

[0011] In view of the above, the present invention provides a method for improving the coincidence time resolution of a PET system based on STFT, which can effectively improve the coincidence time resolution of the PET system.

[0012] A method for improving the temporal resolution of a PET system based on STFT comprises the following steps:

[0013] (1) Using a pair of PET detectors to detect coincidence events occurring at a point source, and then saving the detected coincidence waveform and the corresponding point source position;

[0014] (2) trimming the conforming waveforms and calculating the true TOF value corresponding to each set of conforming waveforms;

[0015] (3) Performing short-time Fourier transform (STFT) on each set of coincident waveforms to obtain their short-time frequency domain amplitude information;

[0016] (4) obtaining a large number of samples according to steps (1) to (3), wherein each set of samples includes short-time frequency domain amplitude information that conforms to the waveform and the true TOF value, and then dividing all samples into a training set and a test set;

[0017] (5) Construct a residual convolutional neural network model, which consists of multiple cascaded residual units and fully connected layers;

[0018] (6) Using the short-time frequency domain amplitude information in the training set samples as the model input and the TOF true value as the label, the above network model is trained;

[0019] (7) By inputting the short-time frequency domain amplitude information in the test set samples into the trained network model, the corresponding TOF time of the output can be directly predicted.

[0020] Furthermore, the specific implementation method of step (1) is as follows: a radioactive point source is placed on the line connecting a pair of PET detectors, the point source is moved on the line at a certain step interval, and the two detectors are used to detect the coincidence event occurring at the point source at each position to obtain multiple groups of coincidence waveforms, each group of coincidence waveforms includes two waveform sequences obtained by the two detectors in response to the paired gamma photons emitted by the same coincidence event, and the time delay between the two waveform sequences represents the PET flight time.

[0021] Furthermore, the clipping process in step (2) is to take the maximum value of each of the two waveform sequences that meet the waveform, and take the maximum value point as the benchmark, sample forward for a time of t1, and sample backward for a time of t2 to form a waveform of length t1+t2, so as to ensure that most of the waveforms obtained after clipping are rising edges of the waveform; if the point source position corresponding to any group of waveforms that meet the waveform is Δx, then the TOF true value of the waveform that meets the waveform is c is the speed of light.

[0022] Furthermore, in step (3), for the two waveform sequences f1(n) and f2(n) that conform to the waveform, short-time Fourier transform is performed on f1(n) and f2(n) using the following formula;

[0023]

[0024]

[0025] Where: F1(n, ω) and F2(n, ω) are the results of short-time Fourier transform of f1(n) and f2(n) respectively, and both are two-dimensional arrays. The horizontal axis n represents the time, the vertical axis ω represents the frequency, h() is the window function, j is the imaginary unit, and m represents the time shift during Fourier transform.

[0026] Then, the amplitudes S1(n, ω) and S2(n, ω) of F1(n, ω) and F2(n, ω) are extracted respectively, and then S1(n, ω) and S2(n, ω) are concatenated up and down to serve as the short-time frequency domain amplitude information that conforms to the waveform.

[0027] Furthermore, the residual unit is composed of three residual modules D1 to D3 connected in sequence. Each residual module is composed of a convolution layer, a batch normalization layer, and an activation function ReLU connected in sequence from input to output. The convolution kernel size of the convolution layer in D1 and D3 is 1×1, and the convolution kernel size of the convolution layer in D2 is 3×3. The output of D3 is added to the input of D1 and then processed by the activation function ReLU to serve as the output of the residual unit.

[0028] Furthermore, the output of each layer of residual units in the network model serves as the input of the next layer of residual units. The input of the first layer of residual units is the short-time frequency domain amplitude information in the training set samples, which is expanded into a one-dimensional vector after passing through multiple layers of residual units. After passing the one-dimensional vector through the fully connected layer, the TOF time corresponding to the waveform is output.

[0029] Furthermore, the process of training the network model in step (6) is as follows:

[0030] 6.1 Initialize model parameters, including the bias vector and weight matrix of each layer, learning rate, and optimizer;

[0031] 6.2 Input the short-time frequency domain amplitude information of the training set samples into the model, and the model forward propagates the output to obtain the corresponding prediction result, i.e., the TOF time, and calculate the loss function L between the prediction result and the label;

[0032] 6.3 Based on the loss function L, the optimizer is used to iteratively update the model parameters using the gradient descent method until the loss function L converges and the training is completed.

[0033] Furthermore, the loss function L adopts mean square error.

[0034] Furthermore, the optimizer adopts Adam optimizer.

[0035] This paper uses a residual neural network model based on the STFT to estimate the time of flight (TOF) of a PET system, extracting potential frequency-domain features from the waveform to facilitate better timing. This enables the PET system to achieve better temporal resolution, thereby leveraging TOF information to achieve better spatial resolution. This enables the PET system to provide more accurate information in medical imaging, further aiding clinical diagnosis. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] Figure 1 Schematic diagram of the PET system's detection of event waveforms and TOF true value data.

[0037] Figure 2 Schematic diagram of the overall structure of the network model of the present invention.

[0038] Figure 3 Schematic diagram of the structure of the residual unit in the network model of the present invention.

[0039] Figure 4 Schematic diagram comparing the TOF prediction results of the present invention and the traditional method at point source positions of -50mm, 0mm, and 50mm. DETAILED DESCRIPTION

[0040] In order to describe the present invention more specifically, the technical solution of the present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0041] The STFT-based PET system of the present invention complies with a temporal resolution improvement method, comprising the following steps:

[0042] (1) Place a radioactive point source on the line connecting a pair of PET detectors. Move the point source at a certain interval and collect waveforms from the two detectors when a coincidence event occurs. Save the point source position Δx and waveform for each set of coincidence events.

[0043] The point source is located on the line connecting the two detectors. After collecting enough coincident events each time, it moves along the line connecting the two detectors with a certain step size and continues to collect waveform data of the next batch of coincident events.

[0044] (2) For each matching event, each pair of waveforms collected is trimmed. Taking the time t0 when the waveform reaches its highest value as the reference point, waveforms with durations t1 and t2 are taken forward and backward to form a pair of waveform data with a length of t1 + t2. Based on the saved point source position Δx, the true TOF value Δt is calculated using the following formula:

[0045]

[0046] Where: c is the speed of light.

[0047] (3) Perform short-time Fourier transform (STFT) on each pair of waveform sequences f1(n) and f2(n), as follows:

[0048]

[0049] Where: f(n) is the waveform sequence, h(n) is the window function (Hanning window, etc. can be selected), F(n, ω) is a two-dimensional array, the horizontal axis is time n, and the vertical axis is ω, which is the frequency.

[0050] Then the amplitude S(n, ω) of F(n, ω) is extracted as the short-time frequency domain amplitude information of the waveform, as shown in the following formula:

[0051] S(n,ω)=|F(n,ω)|

[0052] The two two-dimensional arrays of frequency domain amplitude information S1(n, ω) and S2(n, ω) after short-time Fourier transform of f1(n) and f2(n) are concatenated up and down as the waveform short-time frequency domain amplitude information S of the event.

[0053] (4) All collected waveform short-time frequency domain amplitude information and TOF true values ​​that meet the event are randomly divided into training set, validation set and test set; the data set division requirements are: no duplication, and the training set, validation set and test set are divided in a ratio of 7:2:1.

[0054] (5) Extract the event data from the training set, use the waveform short-time frequency domain amplitude information S as the input sample, and the TOF true value Δt as the true value label to train the residual convolutional neural network. Finally, according to the performance of the model on the validation set, the TOF time estimation model of the PET system based on STFT is obtained. Specifically:

[0055] 5.1 Construct a residual neural network model, which consists of multiple layers of residual units and fully connected layers. Each residual unit is composed of multiple residual modules connected together, and each residual module is composed of a convolutional layer, an activation function, and a batch normalization layer connected together.

[0056] 5.2 Initialize the parameters of the above network model, including the bias vector and weight matrix of each residual module and fully connected layer;

[0057] 5.3 The residual neural network model is trained under the guidance of the TOF true value label of the training set, and the training status of the model is verified by the validation set. Specifically: the short-time frequency domain amplitude information of the waveform that meets the event in the training set is input into the residual neural network model to obtain the prediction result Δt′ of the residual neural network model, and the mean square error L between the output result and the true value label Δt is calculated, and the parameters of each residual module are optimized by the Adam optimizer. After each step of training, the short-time frequency domain amplitude information of the waveform that meets the event in the validation set is used as input to obtain the prediction value Δt″ of the residual neural network model and calculate its mean square error L' with the true value label; when the error of the model on the validation set no longer decreases, the final residual neural network model is obtained.

[0058] The expression of the error function L is:

[0059]

[0060] The expression of the error function L' is:

[0061]

[0062] in: is the square of the L2 norm.

[0063] (5) Take out the event data in the test set, use the waveform short-time frequency domain amplitude information as the input sample to input into the trained network model, and output the TOF time value estimated by the network model.

[0064] In the following embodiments, we use Monte Carlo simulation to obtain waveform data, which is consistent with the detection of event data. Figure 1 In the simulation, the PET detector consists of a LYSO crystal and a SiPM photodetector. The center of the line connecting the pair of PET detectors is taken as the origin, and the line is the z-axis. Point sources are placed every 1 cm and 3,000 sets of coincidence event data are collected. A total of 63,000 sets of coincidence event data are collected at 21 locations.

[0065] Each pair of acquired waveforms is cropped, with time t0, when the waveform reaches its peak, as the reference point. Waveforms of 3.5 ns and 1.5 ns are taken forward and backward, respectively, to form a pair of waveform data f1(n) and f2(n) with a length of 5 ns. Because the waveform sampling interval is 0.05 ns, each pair of waveform data is one-dimensional data with a length of 100.

[0066] Then, the short-time Fourier transform (STFT) is performed on the pair of waveform data f1(n) and f2(n), respectively, as shown in the following formula:

[0067]

[0068] Where: f(n) is the corresponding discrete waveform data, h(n) is the window function, and F(n, ω) is a two-dimensional array with the horizontal axis being time n and the vertical axis being ω, i.e., frequency.

[0069] Then, the amplitude S(n, ω) of F(n, ω) is extracted as the short-time frequency domain amplitude information of the waveform, as shown in the following formula:

[0070] S(n,ω)=|F(n,ω)|

[0071] Then, the two two-dimensional arrays of frequency domain amplitude information S1(n, ω) and S2(n, ω) after short-time Fourier transform of f1(n) and f2(n) are concatenated up and down as the waveform short-time frequency domain amplitude information S of the event.

[0072] Finally, according to the point source position Δx, the true TOF value Δt is calculated using the following formula:

[0073]

[0074] Where: c is the speed of light.

[0075] At this point, a coincidence event data consists of a two-dimensional array S representing the short-time frequency domain amplitude information of the waveform and a TOF true value Δt.

[0076] The 63,000 sets of event data collected were divided into training set, validation set and test set in a ratio of 7:2:1.

[0077] Construct residual modules; each residual module consists of a convolutional layer, an activation function, and a batch normalization layer, such as Figure 3 As shown in the figure, each residual module consists of three small modules, each consisting of a convolutional layer, a batch normalization layer, and a ReLU activation function. The convolution kernel size of the first and third small modules is 1×1, and the convolution kernel size of the second small module is 3×3. The short-time frequency domain amplitude signal S passes through the three small modules and is then added to itself. Finally, it passes through the ReLU activation function layer and serves as the output of the entire residual module.

[0078] Construct a residual neural network model; the model consists of three layers of residual modules and a fully connected layer, such as Figure 2 As shown, the residual modules 1, 2, and 3 have the same structure as Figure 3 As shown in the figure, the short-time frequency domain amplitude signal S is expanded into a one-dimensional vector after passing through three residual modules, and the one-dimensional vector is passed through a linear layer to obtain the final TOF time prediction value Δt′.

[0079] The residual neural network model is trained under the guidance of the true value labels of the training set, and the training status of the model is verified by the validation set. Specifically, the short-time frequency domain amplitude information of each waveform that meets the event in the training set is input into the model to obtain the model's prediction result Δt', calculate the mean square error L between the prediction result and the true value label Δt, and optimize the parameters of each residual module and linear layer through the Adam optimizer. After each step of training, the short-time frequency domain amplitude information of the waveform that meets the event in the validation set is used as input to obtain the model's prediction value Δt" and calculate its mean square error L' with the true value label; when the error of the model on the validation set no longer decreases, the final STFT-based residual neural network model is obtained, where the mean square error L and L' are calculated as follows:

[0080]

[0081]

[0082] in: is the square of the L2 norm.

[0083] Finally, the event data in the test set is taken out and input into the trained residual neural network model, and the TOF time value predicted by the model is output; the time difference distribution spectrum of the TOF time value at each position is statistically predicted, where the time difference distribution spectrum at -50mm, 0mm, and 50mm is as follows Figure 4 As shown, the comparison method is the constant fraction timing method (CFD). It can be clearly seen that the time difference distribution spectrum predicted by the model of the present invention is significantly more concentrated, and its full width at half maximum (FWHM) is lower.

[0084] The full width at half maximum of the time difference distribution between the STFT-based model of the present invention and the traditional CFD method from -100 mm to 100 mm is shown in Table 1:

[0085] Table 1

[0086] Position / mm CFD method FWHM / ps STFT method FWHM / ps Relative improvement -100 99.82 60.73 39.16％ -90 105.38 75.79 28.09％ -80 108.66 77.65 28.54％ -70 103.86 84.49 18.64％ -60 108.70 82.88 23.75％ -50 106.90 83.56 21.83％ -40 108.85 82.16 24.52％ -30 112.60 84.97 24.54％ -20 114.47 77.05 32.69％ -10 113.69 85.88 24.46％ 0 117.59 85.56 27.24％ 10 111.01 78.15 29.60％ 20 123.31 84.48 31.49％ 30 104.15 76.16 26.87％ 40 114.79 85.46 25.55％ 50 105.94 82.25 22.36％ 60 107.43 79.37 26.12％ 70 110.32 80.60 26.94％ 80 111.26 79.09 28.91％ 90 114.50 72.49 36.68％ 100 107.30 53.46 50.18％

[0087] From the above experimental results, we can see that the STFT-based TOF time estimation method for the PET system of the present invention effectively improves the half-maximum full width of the PET system at various positions and improves the time resolution of the PET system.

[0088] The above description of the embodiments is intended to facilitate understanding and application of the present invention by those skilled in the art. It is apparent that those skilled in the art can readily make various modifications to the above embodiments and apply the general principles described herein to other embodiments without requiring creative effort. Therefore, the present invention is not limited to the above embodiments. Any improvements or modifications made by those skilled in the art based on the disclosure of the present invention should fall within the scope of protection of the present invention.

Claims

1. A method for improving the temporal resolution of a PET system based on STFT, comprising the following steps: (1) Using a pair of PET detectors to detect coincidence events occurring at a point source, and then saving the detected coincidence waveform and the corresponding point source position; (2) trimming the conforming waveforms and calculating the true TOF value corresponding to each set of conforming waveforms, where TOF stands for time of flight; (3) Performing STFT transformation on each set of coincident waveforms to obtain their short-time frequency domain amplitude information. STFT transformation is short-time Fourier transform. (4) obtaining multiple groups of samples according to steps (1) to (3), each group of samples including short-time frequency domain amplitude information that conforms to the waveform and the true TOF value, and then dividing all samples into a training set and a test set; (5) Construct a residual convolutional neural network model, which consists of multiple cascaded residual units and fully connected layers; (6) Using the short-time frequency domain amplitude information in the training set samples as the model input and the TOF true value as the label, the above network model is trained; (7) By inputting the short-time frequency domain amplitude information in the test set samples into the trained network model, the corresponding TOF time of the output can be directly predicted.

2. The method for improving temporal resolution of a PET system according to claim 1, wherein: The specific implementation method of step (1) is as follows: a radioactive point source is placed on the line connecting a pair of PET detectors, the point source is moved on the line at a certain step interval, and the coincidence events occurring at the point source are detected by the two detectors at each position to obtain multiple groups of coincidence waveforms. Each group of coincidence waveforms includes two waveform sequences obtained by the two detectors in response to the paired gamma photons emitted by the same coincidence event. The time delay between the two waveform sequences represents the PET flight time.

3. The method for improving temporal resolution of a PET system according to claim 1, wherein: The clipping process in step (2) is to take the maximum value of each of the two waveform sequences that meet the waveform, and take the maximum value point as the reference, sample forward for a time of t1, and sample backward for a time of t2 to form a waveform of length t1+t2, so as to ensure that most of the waveforms obtained after clipping are rising edges of the waveform; if the point source position corresponding to any group of waveforms that meet the waveform is Δx, then the TOF true value of the waveform that meets the waveform is Δt and c is the speed of light.

4. The method for improving temporal resolution of a PET system according to claim 3, wherein: In step (3), for the two waveform sequences f1(n) and f2(n) that conform to the waveform, STFT transformation is performed on f1(n) and f2(n) using the following formula; Where: F1(n,ω) and F2(n,ω) are the results of STFT transformation of f1(n) and f2(n) respectively, and both are two-dimensional arrays. The horizontal axis n represents the time, the vertical axis ω represents the frequency, h() is the window function, j is the imaginary unit, and m represents the time shift during Fourier transform. Then, the amplitudes S1(n,ω) and S2(n,ω) of F1(n,ω) and F2(n,ω) are extracted respectively, and then S1(n,ω) and S2(n,ω) are concatenated up and down to obtain the short-time frequency domain amplitude information that conforms to the waveform.

5. The method for improving temporal resolution of a PET system according to claim 1, wherein: The residual unit is composed of three residual modules D1, D2, and D3 connected in sequence. Each residual module is composed of a convolution layer, a batch normalization layer, and an activation function ReLU connected in sequence from input to output. The convolution kernel size of the convolution layer in D1 and D3 is 1×1, and the convolution kernel size of the convolution layer in D2 is 3×3. The output of D3 is added to the input of D1 and then processed by the activation function ReLU to serve as the output of the residual unit.

6. The method for improving temporal resolution of a PET system according to claim 1, wherein: The output of each layer of residual units in the network model serves as the input of the next layer of residual units. The input of the first layer of residual units is the short-time frequency domain amplitude information in the training set samples, which is expanded into a one-dimensional vector after passing through multiple layers of residual units. After passing the one-dimensional vector through the fully connected layer, the TOF time corresponding to the waveform is output.

7. The method for improving temporal resolution of a PET system according to claim 1, wherein: The process of training the network model in step (6) is as follows: 6.1 Initialize model parameters, including the bias vector and weight matrix of each layer, learning rate, and optimizer; 6.2 Input the short-time frequency domain amplitude information of the training set samples into the model, and the model forward propagates the output to obtain the corresponding prediction result, i.e., the TOF time, and calculate the loss function L between the prediction result and the label; 6.3 Based on the loss function L, the optimizer is used to iteratively update the model parameters using the gradient descent method until the loss function L converges and the training is completed.

8. The method for improving temporal resolution of a PET system according to claim 7, wherein: The loss function L adopts mean square error.

9. The method for improving temporal resolution of a PET system according to claim 7, wherein: The optimizer adopts Adam optimizer.

Images