A crystal level pet time correction method based on l1 regularization term

Abstract

The application provides a crystal level PET time correction method based on an L1 regular term, which changes a PET time correction method based on a response line number mode to establish a system matrix, adopts a coincidence event number mode to establish the system matrix, and realizes crystal level correction from original detector block level correction or sub-crystal level correction. Meanwhile, the application adopts a least square QR decomposition method to speed up according to the sparsity of the system matrix, and adds an L1 regular term in the method to improve the robustness of the method. Through the above mode, the application makes the PET system have better resolution, and further obtains a higher quality reconstructed image, provides more accurate information in medical imaging, and provides better help for clinical diagnosis.

Description

Technical Field

[0001] This invention belongs to the field of PET imaging technology, specifically relating to a crystal-level PET time correction method based on L1 regularization terms. Background Technology

[0002] Since the advent of PET, numerous improvements have been made to its various components, and the performance requirements for each component have been continuously increasing. Time-of-flight (TOF) PET—TOF-PET—was born in this context. For a TOF-PET system, temporal resolution is a crucial parameter, directly affecting the quality of subsequent reconstructed images. Reconstructed images are obtained from coincidence events through data transformation. For coincidence events, the size of the coincidence event window affects the number of coincidence events used for subsequent processing. When the system's temporal resolution is high, a narrower coincidence time window can be used to improve the quality of subsequent reconstructed images by reducing the impact of noisy data. A higher signal-to-noise ratio (SNR) in the reconstructed image indicates a stronger ability to detect smaller tumors and more accurately distinguish the stage of the lesion, thus providing stronger information support for doctors' diagnosis and decision-making. Simultaneously, without compromising image quality, lower doses of radionuclides can be used for image reconstruction. Reducing the dose reduces the risk of radionuclide exposure for patients and increases the imaging rate. Furthermore, reducing the dosage injected into the organism to a certain order of magnitude will allow PET technology to be extended to the treatment of conditions other than cancer, such as inflammation, cardiovascular disease, and sepsis. In conclusion, temporal resolution is closely related to the quality of the reconstructed image; therefore, temporal correction of the data is of paramount importance.

[0003] For PET coincidence time correction, methods can be categorized into direct correction and indirect correction. Direct correction primarily employs a reference detector method. This method uses a fast detector as a reference detector to record the time difference between the measured detector unit and the reference detector, and calculates the delay compensation for each detector unit's detection channel. A fast detector is formed by coupling a fast crystal BaF2 with a photomultiplier tube (PMT). The PMT circuit is designed as a constant fraction discriminator (CFD). Placing the point source Na-22 at the center of the cylindrical PET system's field of view effectively reduces the influence of the Time-of-Flight (TOF) effect. After a period of data acquisition, single-event data streams in list-mode format are collected. Each data stream contains position, energy, and time information. After data acquisition, the fast detector coincides with all crystal detector units within the detector system to obtain the corresponding time difference, which is then used as the reference time.

[0004] Indirect methods include iterative methods, special source methods, linear equation methods, and compatible equation methods, as detailed below:

[0005] Iterative methods; the iterative time correction method calculates the coincidence time difference between a detector element to be corrected and its opposite detector element within the detector's field of view, sums all coincidence time differences, and obtains the time compensation value for the corrected detector element. The coincidence time spectrum of the iterative method is based on the number of response lines; therefore, the detector element can be at the detector block level or at the crystal level as the smallest detector element. Considering the radial direction, multiple detectors form a ring, and the time difference values ​​within the corresponding sector of the detector element are averaged. Similarly, if the axial direction is considered, the values ​​on the axial slices also need to be averaged. Through multiple numerical iterations, the required compensation factor is obtained.

[0006] Specialized radiation source method: In experiments used for time correction, point sources are commonly used as radiation sources. To make the correction results robust, many researchers have proposed using specially designed radiation sources and applying them in different scenarios. Common specialized sources include: rotating line sources with fixed devices, point sources with brass cylindrical blocks for scattering, phantoms, positron emission tomography (PET) devices, and lesion data, etc. The acquired data are numerically solved using iterative methods to obtain the compensation factor.

[0007] The linear equation method involves extracting the coincidence time difference and the position information of the two detection units from the coincidence event flow to establish a system matrix and a measurement time difference vector, thereby obtaining the compensation factor. In terms of details, when data is sufficient, the system matrix can be established based on the number of response lines or the number of coincidence events. Both methods can be used for solving the problem. In earlier methods, the linear equation method based on the number of response lines was more common; however, this method is significantly limited by the amount of data and therefore makes it difficult to achieve crystal-level correction. Relatively speaking, in scenarios with low noise, the linear equation method based on the number of coincidence events is recommended. Summary of the Invention

[0008] In view of the above, the present invention provides a crystal-level PET time correction method based on L1 regularization terms, which can effectively improve the time resolution of PET systems.

[0009] A crystal-level PET time correction method based on L1 canonical terms includes the following steps:

[0010] (1) Use the PET detection system to detect coincidence events occurring at the radioactive source and determine the time resolution of the PET detection system (for comparison with the corrected time resolution of the PET detection system);

[0011] (2) Establish a PET time correction model based on the linear equation method;

[0012] (3) Based on the above PET time correction model, add L1 regularization constraints and redefine the objective function for PET time correction;

[0013] (4) Transform the above objective function into an optimization problem with constraints, and solve it;

[0014] (5) The crystal compensation value vector obtained above is added to the subsequent calculation process of the event timestamp.

[0015] Furthermore, the radiation source used in step (1) can be a point source or a cylindrical source, which is placed at the center of the PET detection system. The detected coincidence events include immediate coincidence events and delayed coincidence events, and the delayed coincidence events need to be removed from the immediate coincidence events.

[0016] Furthermore, the time resolution of the PET detection system obtained in step (1) is presented as an FWHM (Full Width at Half Maxima) plot.

[0017] Furthermore, the PET time correction model expression in step (2) is as follows:

[0018]

[0019] Where: |||2 represents the L2 norm, A is the system matrix, x is the crystal compensation value vector used for PET time correction, and b is the detection delay vector.

[0020] Furthermore, the system matrix A is an m×n matrix, where m is the number of detected coincidence events and n is the number of detection units in the PET detection system. Each row of the system matrix A corresponds to a specific coincidence event. In a row, the values ​​corresponding to two detection units i and j are 1 and -1, respectively, and the remaining values ​​are 0, indicating that the coincidence event corresponding to that row was detected by detection units i and j. The detection delay vector b is an m-dimensional vector, and the crystal compensation value vector x is an n-dimensional vector. Each row of vector b corresponds to each coincidence event, that is, it represents the detection time difference Δt between the two detection units corresponding to the coincidence event.

[0021] Furthermore, the objective function expression in step (3) is as follows:

[0022]

[0023] Where: |||2 represents the L2 norm, |||1 represents the L1 norm, A is the system matrix, x is the crystal compensation value vector used for PET time correction, b is the detection delay vector, and λ is the regularization term weight parameter.

[0024] Furthermore, in step (4), the objective function is transformed into an optimization problem with constraints, specifically as follows:

[0025] min f(x) + g(z)

[0026] stx-z = 0

[0027]

[0028] Where: |||2 represents the L2 norm, |||1 represents the L1 norm, A is the system matrix, x is the crystal compensation value vector used for PET time correction, b is the detection delay vector, λ is the regularization term weight parameter, and z is the intermediate variable.

[0029] Furthermore, the specific process of solving the optimization problem in step (4) is as follows:

[0030] x k+1 =(A T A+ρI) -1 (A T b+ρ(z k -uk ))

[0031] z k+1 =S λ / ρ (x k+1 +u k )

[0032] u k+1 =u k +x k+1 -z k+1

[0033] Where: A is the system matrix, I is the identity matrix, x is the crystal compensation value vector used for PET time correction, b is the detection delay vector, ρ is the regularization term weighting parameter, and z and u are intermediate variables. T This indicates transpose, with superscripts k and k+1 representing the iteration count, where k is a natural number, and S λ / ρ () represents the soft threshold operator.

[0034] This invention modifies the conventional PET time correction method based on the number of response lines to establish the system matrix. Instead, it uses the number of coincidence events to establish the system matrix, achieving crystal-level correction instead of the previous detector block-level or sub-crystal-level correction. Furthermore, this invention utilizes least-squares QR decomposition to accelerate the process due to the sparsity of the system matrix and adds an L1 regularization term to improve robustness. Through these methods, this invention enables PET systems to achieve better resolution, resulting in higher-quality reconstructed images. This provides more accurate information in medical imaging and better assists in clinical diagnosis. Attached Figure Description

[0035] Figure 1 This is a schematic diagram of the HITS-655K PET system used in an embodiment of the present invention.

[0036] Figure 2 This is a schematic flowchart of the crystal-level PET time correction method based on L1 regularization term of the present invention.

[0037] Figure 3 This is a schematic diagram of the steps in the PET time correction method of the present invention to solve the problem using the alternating direction multiplier method. Detailed Implementation

[0038] To describe the present invention in more detail, the technical solution of the present invention will be described in detail below with reference to the accompanying drawings and specific embodiments.

[0039] The PET detection system used in this embodiment is the HITS-655K brain PET scanner from Hamamatsu Photonics Co., Ltd., Japan. Figure 1As shown, the specific parameters of the system can be found in the reference [Watanabe M, Saito A, Isobe T, et al. Performance evaluation of a high-resolution brain PET scanner using four-layer MPPC DOI detector[J]. Physics in Medicine & Biology, 2017, 62(17)]. The scanner consists of four layers of detector modules, representing four depths of the detector. The lengths of the scintillation crystals in the four layers are 3mm, 4mm, 5mm, and 8mm from the inside out, respectively, with a total detector depth of 45mm. For the axial direction angle, the system consists of five ring detectors, with each small detector having 32 rings. For the radial direction angle, the system consists of 32 four-layer detector modules. The scintillation crystals used in the detection system are LYSO, with each crystal array consisting of 256 crystal detector units, each 32×32 in size. The size of a single crystal detector unit is 1.2mm×1.2mm, and the crystal length is related to the number of layers in its array. Therefore, the total number of crystal detector units in the scanner is 32 (radial number of crystal arrays) × 32 (axial number of crystal arrays) × 32 (number of modules) × 4 (depth number) × 5 (number of small rings) = 655360. The photosensitive device used is a Hamamatsu S10931-050P multi-pixel photon counter. The effective area of ​​a single device is 3mm × 3mm, with each pixel spaced 50um, for a total of 3600 pixels. 64 MPPCs are arranged in an 8×8 array and coupled to the crystal array to form a detector.

[0040] The flowchart of the crystal-level PET time correction method based on L1 canonical terms of this invention is as follows: Figure 2 As shown, the time correction model based on the linear equation is first established as follows:

[0041]

[0042] Where: || ||2 represents the L2 norm, A is the system matrix, x is the crystal compensation value vector, and b is the measurement time difference vector;

[0043] In the above model, the system matrix A is an m×n matrix, where m is the number of coincidence events collected and n is the number of detection units. Each row of the system matrix A represents a specific coincidence event. In each row, there will be two detection units i and j with column values ​​of 1 and -1 respectively, and the remaining column values ​​will be 0, indicating that the coincidence event is between detection unit i and detection unit j. The detection time difference vector b is an m×1 vector, and the crystal compensation value vector x is an n×1 vector. The value Δt in each row of b corresponds to each row of the system matrix A, representing the coincidence event, where the detection time difference between detection unit i and detection unit j is Δt.

[0044] To make the system matrix more intuitive, its specific form is shown below:

[0045]

[0046] Where: u i and u j Let represent the i-th and j-th detection units respectively. The corresponding vector b represents the time difference between the two detection units. In the above formula, the first event represents the detection of detection unit 1 and detection unit 15, which undergoes coincidence processing. The time difference is 1 unit time difference, which is related to the measurement of different detectors, and is generally expressed in ns or ps. The sign of the time difference is determined by subtracting the next detection unit from the previous one. Similarly, the second event represents the detection unit 1 and detection unit 16, with a time difference of -3. And so on. Each coincidence event is presented in the form of a row vector, and the vector b contains the time difference value of the coincidence event.

[0047] A PET detector typically consists of a photodetector and a crystal array coupled together, and the number of detector units, *n*, is related to the crystal region division. Taking a 16×16 crystal array as an example, the parameter *section* represents the number of region divisions. When *section* is 1, the entire 16×16 array is considered a single crystal detector unit. When *section* is 2, the 16×16 crystal is divided into 4 (the square of the *section* value) 8×8 crystal detector units of the same size, using the midline of the two sides as the dividing line. When *section* is 4, it is divided into 16 4×4 crystal detector units of the same size, and so on, with *section* values ​​of 4, 8, ... The maximum value of *section* represents the minimum number of crystal detector units; in the example given, this is 16. The number of detector units, *n*, is the square of the *section* value. However, in this invention, since the correction effect requires crystal-level correction, each crystal in the array is considered the minimum detector unit, provided storage resources allow.

[0048] The system matrix and the measurement time difference are obtained, and the linear equation is inverted to obtain the compensation factor. To improve the robustness of the results against noisy data, as follows... Figure 3 As shown, this invention introduces an L1 regularization term into the linear equation, transforming the objective function into a constrained optimization problem, and redefining the problem as follows:

[0049]

[0050] Where: || ||1 represents the L1 norm, and λ is the regularization term weight parameter;

[0051] Solving the above problem is equivalent to solving a constraint problem:

[0052] min f(x) + g(z)

[0053] stx-z = 0

[0054] in: g(z) = λ||z||1, where z is an intermediate variable.

[0055] This invention uses the Alternating Direction Multiplier Method (ADMM) to solve the above problem. The constrained optimization problem is formally rewritten using ADMM, resulting in the following form:

[0056] x k+1 :=(A T A+ρI) -1 (A T b+ρ(z k -u k ))

[0057] z k+1 :=S λ / ρ (x k+1 +u k )

[0058] u k+1 :=u k +x k+1 -z k+1

[0059] Where: ρ is the weight parameter, I is the identity matrix, and S... λ / ρ This is a soft threshold operation.

[0060] The x-step update requires solving the ridge regression calculation. Unlike the previous inversion method, given the sparsity of the system matrix A, the least squares QR decomposition method (LSQR) can be used to solve it, which can improve the calculation speed.

[0061] The soft thresholding operation involved in the z-step update, which is used to solve the L1 regularization term, is specifically represented as follows:

[0062]

[0063] The following embodiment verifies the effectiveness of the method using data collected from experiments with a radioactive point source on the HITS-655K system. A total of 5,330,535 coincidence events were collected for calibration, as shown in Table 1:

[0064] Table 1

[0065]

[0066] Based on the above experiments, this invention, by using a sparse matrix for storage, achieves crystal-level correction, whereas previous methods could only achieve sub-crystal-level correction. Therefore, it can be applied to large-scale PET detector systems.

[0067] To verify the effectiveness of the L1 regularization term, we conducted experiments on a small detector single-ring, and Table 2 shows the experimental results. It can be observed that the FWHM obtained by the method based on the L1 regularization term is better than that obtained by the original LS method.

[0068] Table 2

[0069]

[0070] Based on the above experiments, and by comparing with traditional methods, it can be seen that the PET system time correction method based on L1 regularization term of this invention effectively improves the accuracy and enhances the time resolution of the PET system.

[0071] The above description of the embodiments is provided to enable those skilled in the art to understand and apply the present invention. 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 creative effort. Therefore, the present invention is not limited to the above embodiments, and any improvements and modifications made to the present invention by those skilled in the art based on the disclosure thereof should be within the scope of protection of the present invention.

Claims

1. A crystal-level PET time correction method based on L1 canonical terms, comprising the following steps: (1) Use the PET detection system to detect coincidence events occurring at the radioactive source and determine the time resolution of the PET detection system; (2) A PET time correction model based on the linear equation method is established, and its expression is as follows: in: Describing the L2 norm, A For the system matrix, x This is the crystal compensation value vector used for PET time correction. b To detect the delay vector; (3) Based on the above PET time correction model, add L1 regularization constraints and redefine the objective function for PET time correction, the expression of which is as follows: in: Describing the L1 norm, λ These are the regularization term weight parameters; (4) The objective function above is transformed into an optimization problem with constraints, as shown below: ， The optimization problem is then solved, and the specific process is as follows: in: I It is the identity matrix. ρ For regularization term weight parameters, z and u As an intermediate variable, T Indicates transpose, superscript k and k +1 indicates the number of iterations. k For natural numbers, S λ / ρ ( ) represents the soft threshold operator; (5) The crystal compensation value vector obtained above is added to the subsequent calculation process of the event timestamp.

2. The crystal-level PET time correction method according to claim 1, characterized in that: The radiation source used in step (1) is a point source or a cylindrical source, which is placed at the center of the PET detection system. The detected coincidence events include immediate coincidence events and delayed coincidence events. Delayed coincidence events need to be removed from immediate coincidence events.

3. The crystal-level PET time correction method according to claim 1, characterized in that: The time resolution of the PET detection system obtained in step (1) is presented in the form of an FWHM diagram.

4. The crystal-level PET time correction method according to claim 1, characterized in that: The system matrix A for m × n A matrix of size, where m The number of matching events detected. n The number of detection units in the PET detection system, system matrix A Each row corresponds to a specific coincidence event, and there are two probe units in a row. i and j The corresponding values ​​are 1 and -1, and the rest are 0, indicating that the matching event for that row was detected by the detection unit. i and j Detected; The detection delay vector b for m dimensional vector, the crystal compensation value vector x for n dimensional vector, vector b Each row of values ​​in the table corresponds to a coincidence event, representing the detection time difference Δ between the two detection units corresponding to the coincidence event. t .

Images