CN115330900B - System matrix fast iterative computation method for PET image reconstruction - Google Patents

System matrix fast iterative computation method for PET image reconstruction Download PDF

Info

Publication number
CN115330900B
CN115330900B CN202211252651.0A CN202211252651A CN115330900B CN 115330900 B CN115330900 B CN 115330900B CN 202211252651 A CN202211252651 A CN 202211252651A CN 115330900 B CN115330900 B CN 115330900B
Authority
CN
China
Prior art keywords
voxel
ray
distance
algorithm
representing
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202211252651.0A
Other languages
Chinese (zh)
Other versions
CN115330900A (en
Inventor
张弓
苏进
李学俊
王华彬
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Canada Institute Of Health Engineering Hefei Co ltd
Original Assignee
China Canada Institute Of Health Engineering Hefei Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Canada Institute Of Health Engineering Hefei Co ltd filed Critical China Canada Institute Of Health Engineering Hefei Co ltd
Priority to CN202211252651.0A priority Critical patent/CN115330900B/en
Publication of CN115330900A publication Critical patent/CN115330900A/en
Application granted granted Critical
Publication of CN115330900B publication Critical patent/CN115330900B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/424—Iterative

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Image Generation (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

The invention discloses a system matrix fast iterative computation method for PET image reconstruction, and belongs to the technical field of image processing. Firstly, an LOR line is brought into a siddon algorithm, a first voxel index coordinate passed by the LOR line is calculated, and a voxel neighborhood is obtained; secondly, determining a voxel increment through an LOR line, adding the first voxel passed by the ray with the voxel increment to obtain a voxel to be passed by the ray and searching a neighborhood of the voxel to judge whether the voxel is out of range; if the boundary is not out of range, adding the voxel and the voxel increment to obtain the next voxel, and repeating the step until the boundary is out of range; and finally, substituting the calculated neighborhood voxels into a distance algorithm to obtain the distance from the voxels to the ray, and substituting the distance into a probability formula to obtain a system matrix. The invention optimizes the orthogonal distance ray tracing algorithm and improves the imaging quality.

Description

System matrix fast iterative computation method for PET image reconstruction
Technical Field
The invention relates to the technical field of image processing, in particular to a system matrix fast iterative computation method for PET image reconstruction.
Background
Positron Emission Tomography (PET) is a medical imaging technique that has recognized clinical value for disease diagnosis and therapy monitoring. How to reconstruct high quality PET images quickly and accurately is a hotspot of current research. Usually, an iterative method is used to reconstruct the Rawdata data acquired by the PET device, and in order to accurately calculate the mapping relationship between the image space and the projection space, a system matrix is required to be solved, which reflects the mapping relationship
Figure 945922DEST_PATH_IMAGE001
The contribution of rays to image pixels directly influences the reconstruction speed and the image quality, and is the key of an iterative reconstruction algorithm.
The classical system matrix solving algorithm is the Siddon algorithm proposed by Robert l. Siddon in 1985, which finds all voxels crossed by the LOR line and their intercepts with the LOR line, representing the probability that the voxel is detected by the ray with its intercept with the ray. The increment method changes the multiplication operation in the Siddon algorithm into the addition and subtraction operation, thereby improving the calculation speed. The central pixel law gives the voxel probability value of 1 traversed by a ray, the others being 0. The area weighting algorithm uses the intersection area of the ray and the voxel as the weight of the system matrix, and the detail texture of the reconstructed image is clearer. The orthogonal distance ray tracing algorithm finds voxels around the ray by the siddon algorithm, making the ray have a width.
However, although the siddon algorithm and the incremental method and the center pixel method have high reconstruction speed, the reconstructed image is missing. The area-weighted reconstruction method has excellent quality but slow speed. The quality of the image reconstructed by the orthogonal distance ray tracing algorithm is not as good as that of the area weighting method, but the speed is higher than that of the area weighting method, but the speed is still slower than that of the siddon and the like.
Disclosure of Invention
1. Technical problem to be solved by the invention
In view of the defects of the prior art, the invention provides a system matrix fast iterative calculation method for PET image reconstruction, which optimizes an orthogonal distance ray tracing algorithm, optimizes a siddon algorithm, introduces incremental thinking, recursively calculates voxel indexes, reduces the calculation time of the siddon algorithm, optimizes the query of adjacent voxel neighborhoods by addition and subtraction operation, and optimizes a point-to-straight line distance calculation formula to improve the calculation speed of the orthogonal distance ray tracing algorithm. In addition, when the matrix probability of the system is calculated, the calculation mode of the probability value is optimized, and the imaging quality is improved.
2. Technical scheme
In order to achieve the purpose, the technical scheme provided by the invention is as follows:
the invention relates to a system matrix fast iterative computation method for PET image reconstruction, which comprises the steps of firstly bringing an LOR line into a siddon algorithm, computing a first voxel passed by the LOR line and obtaining a voxel neighborhood; secondly, determining a voxel increment through an LOR line, adding the first voxel passed by the ray with the voxel increment to obtain a voxel to be passed by the ray and searching a neighborhood of the voxel to judge whether the voxel is out of range; if the boundary is not out of range, adding the voxel and the voxel increment to obtain the next voxel, and repeating the step until the boundary is out of range; and finally, substituting the calculated neighborhood voxels into a distance algorithm, solving the distance from the voxels to the rays, and substituting the distance into a probability formula to solve a system matrix.
The method introduces incremental thinking on the basis of the siddon algorithm, directly solves the non-repetitive neighborhood of the current voxel and the last voxel, and accelerates the index calculation of the voxel coordinate; secondly, iteratively solving the distance from the neighborhood voxel to the LOR line based on the distance from the voxel to the LOR line, thereby further improving the calculation speed; finally, the probability value of the voxel completely covered by the detector is set to be a fixed value, and the probability value of other voxels is smaller as the distance from the LOR line is farther, so that the imaging quality of the algorithm is improved. A large number of evaluation experiments are carried out on the resolution prosthesis model and the line source prosthesis model, and the effectiveness of the invention is verified.
3. Advantageous effects
Compared with the prior art, the technical scheme provided by the invention has the following remarkable effects:
(1) The invention relates to a system matrix fast iterative calculation method for PET image reconstruction, which introduces increment thinking, replaces multiplication operation by addition and subtraction, simplifies the calculation of voxel coordinate index, and can directly find the non-repetitive neighborhood of the voxel coordinate index and the previous voxel;
(2) According to the rapid iterative calculation method for the system matrix of the PET image reconstruction, when the distance from a voxel to an LOR line is calculated, the distance from the last voxel to the LOR line is utilized, the calculation of the distance from the adjacent voxel to the LOR line is iteratively calculated, and the calculation speed is increased;
(3) The invention relates to a system matrix fast iterative computation method for PET image reconstruction, which optimizes the probability computation mode of the system matrix, considers the condition that a detector has a certain width and voxels are completely covered by the detector, makes the probability value of the voxels in a certain range from a central line be a constant value, and decreases the probability value larger than the range, thereby improving the imaging quality.
Drawings
FIG. 1 is a diagram of an optimization framework for an orthogonal distance ray tracing algorithm;
FIG. 2 (a) is a resolution prosthesis model;
FIG. 2 (b) is a model of a line source prosthesis;
FIG. 3 is a graph of the effect of iterative reconstruction of resolution prosthesis data using four different algorithms;
FIG. 4 is a graph of the effect of iterative reconstruction of resolution prosthesis data using four different algorithms;
FIG. 5 (a) is a gray value contrast plot of a reconstructed image using the center pixel method, the siddon algorithm, the orthogonal distance ray tracing algorithm, and the algorithm of the present invention;
FIG. 5 (b) is a schematic diagram of gray value change of a reconstructed image by a central pixel method;
FIG. 5 (c) is a schematic diagram of gray value change of a reconstructed image by the siddon algorithm;
FIG. 5 (d) is a schematic diagram showing the gray value variation of the reconstructed image by the orthogonal distance ray tracing algorithm;
FIG. 5 (e) is a schematic diagram of the gray value variation of the reconstructed image by the algorithm of the present invention;
fig. 6 is a diagram of the time it takes to compute the ray-crossing a voxel array of size 216 x 163 using different algorithms, with the abscissa representing the number of rays and the ordinate representing time. Wherein the solid line represents the center pixel method, the dashed line represents the siddon algorithm, the dotted line represents the orthogonal distance algorithm, and the dotted line represents the algorithm of the present invention.
Detailed Description
For a further understanding of the invention, reference should be made to the following detailed description taken in conjunction with the accompanying drawings and examples.
Example 1
In Positron Emission Tomography (PET) studies, iterative reconstruction of PET acquired data is typically performed using an iterative method, where a system matrix reflects a mapping relationship between an image space and a projection space, which is a key of an iterative reconstruction algorithm. The conventional orthogonal distance ray tracing algorithm is complex in calculation and low in efficiency, in order to improve the calculation speed and the imaging quality of the algorithm, the embodiment optimizes the orthogonal distance ray tracing algorithm, and provides a system matrix fast iterative calculation method for PET image reconstruction, wherein a frame diagram of the optimization algorithm is shown in FIG. 1.
This embodiment first brings the LOR line into the siddon algorithm, calculates the first voxel that the LOR line passes through and acquires the voxel neighborhood. Secondly, determining voxel increment through an LOR line, adding the first voxel passed by the ray with the voxel increment to obtain the next voxel to be passed by the ray, searching the neighborhood of the next voxel to be passed by the ray, and judging whether the voxel is out of range or not. If the boundary is not crossed, the voxel is added with the voxel increment to obtain the next voxel, and the step is repeated until the boundary is crossed. And finally, substituting the calculated neighborhood voxels into a distance algorithm to obtain the distance from the voxels to the ray, and substituting the distance into a probability formula to obtain a system matrix.
The algorithm of the embodiment can be divided into three parts, wherein the first part is to use the siddon algorithm to find the first voxel through which a ray passes; the second part is that other voxels and neighborhoods which the ray passes through are obtained on the basis of the first voxel through addition and subtraction; the third part is to utilize the distance formula provided by the invention to calculate the distance from the voxel to the LOR line, and then the calculation result is substituted into the probability value calculation method of the invention. The following will be specifically described:
finding the first voxel crossed by LOR line by siddon algorithm
For one
Figure 665617DEST_PATH_IMAGE002
The voxel array of (1), wherein
Figure 444217DEST_PATH_IMAGE003
The number of planes on the x-axis is indicated,
Figure 108416DEST_PATH_IMAGE004
indicating the number of planes on the y-axis,
Figure 630665DEST_PATH_IMAGE005
number of planes in z-axisThe R.L.Siddon algorithm firstly calculates the x-axis of the three-dimensional space coordinateN x The positions of the planes parallel to each other are calculated and then calculated on the y-axis and the z-axis respectivelyN y AndN z position of the planes parallel to each other:
Figure 204865DEST_PATH_IMAGE006
（1）
wherein the content of the first and second substances,dx、dyanddzrespectively representx、y、zThe distance between the planes in the direction is also the length, width and height of the voxel.
Let the starting point coordinate of the ray be
Figure 29733DEST_PATH_IMAGE007
The coordinate of the end point is
Figure 322174DEST_PATH_IMAGE008
After acquiring the position of the voxel array, the siddon algorithm calculates
Figure 507168DEST_PATH_IMAGE009
And
Figure 935875DEST_PATH_IMAGE010
wherein
Figure 321857DEST_PATH_IMAGE009
Represents the ratio of the distance from the first plane to the ray starting point on the x-axis to the distance of the ray in the x-axis direction,
Figure 711381DEST_PATH_IMAGE011
denotes the number one on the x-axisN x The ratio of the distance from each plane to the starting point of the ray to the distance in the x-axis direction of the ray can be obtained by the same method
Figure 575432DEST_PATH_IMAGE012
. Wherein the content of the first and second substances,
Figure 124225DEST_PATH_IMAGE013
and
Figure 540163DEST_PATH_IMAGE014
representing the point of incidence and the point of emergence, respectively, of a ray through the voxel array:
Figure 541617DEST_PATH_IMAGE015
（2）
wherein:
Figure 99373DEST_PATH_IMAGE016
for the
Figure 502673DEST_PATH_IMAGE017
And
Figure 964878DEST_PATH_IMAGE014
other values of α in between, can also be obtained by different parameter planes. By the formula (1) and the formula (2), respectively, it is possible to obtainx,y,zRange of directional plane indices
Figure 578262DEST_PATH_IMAGE018
,
Figure 49695DEST_PATH_IMAGE019
And
Figure 307501DEST_PATH_IMAGE020
i.e. the plane through which the ray will pass:
Figure 81553DEST_PATH_IMAGE021
Figure 792020DEST_PATH_IMAGE022
（3）
the same can be obtained:
Figure 67143DEST_PATH_IMAGE023
the Siddon algorithm respectively calculates the first time passing of the ray by using the plane index rangex,y,zIn the plane of
Figure 38511DEST_PATH_IMAGE024
The value of the sum of the values,αthe value represents the ratio of the distance value from the ray's intersection with a different plane to the ray's origin to the ray's distance in that direction:
Figure 108098DEST_PATH_IMAGE025
（4）
wherein the content of the first and second substances,λrepresenting a first plane through which the rays passαValue of a step of,
Figure 181227DEST_PATH_IMAGE026
denotes the number one on the x-axisi min A plane and ai max The set of ratios of the distances from the planes to the ray origin to the ray x-axis direction distances,
Figure 994462DEST_PATH_IMAGE027
、
Figure 961281DEST_PATH_IMAGE028
the correspondences represent sets in the y-axis and z-axis,
Figure 326403DEST_PATH_IMAGE029
representing the alpha value of a ray traversing a first plane in the x-axis direction,
Figure 11463DEST_PATH_IMAGE030
representing the value of alpha for a ray traversing a first plane in the y-axis direction,
Figure 238176DEST_PATH_IMAGE031
representing the alpha value of the ray through the first plane in the direction of the z-axis. Since the first voxel a ray traverses is relatively special, the ray may traverse diagonally from the voxel, and will
Figure 325081DEST_PATH_IMAGE032
Is penetrated by raysx,y,zComparing the minimum distance ratios of the planes, and if equal, expressing the minimum distance ratios as
Figure 736470DEST_PATH_IMAGE024
The values are incremented in the respective directions, let:
Figure 33459DEST_PATH_IMAGE033
Figure 188497DEST_PATH_IMAGE034
Figure 129909DEST_PATH_IMAGE035
（5）
Figure 584636DEST_PATH_IMAGE036
、
Figure 244288DEST_PATH_IMAGE037
and
Figure 62071DEST_PATH_IMAGE038
respectively represent LOR lines onx、y、zThe increment between the planes in the direction is,
Figure 123568DEST_PATH_IMAGE007
is the coordinates of the origin of the ray,
Figure 611181DEST_PATH_IMAGE008
is the end-point coordinate of the ray,
Figure 633495DEST_PATH_IMAGE039
representing the alpha value of the ray after increment in the x-axis direction,
Figure 130335DEST_PATH_IMAGE040
the alpha value of the ray after increment in the y-axis direction is shown,
Figure 46338DEST_PATH_IMAGE041
representing the alpha value of the ray after increment in the z-axis direction. Using the method proposed by r.l.siddon to calculate the index coordinates of voxels, the index coordinate of the first voxel traversed by the ray is:
Figure 95066DEST_PATH_IMAGE042
（6）
wherein the content of the first and second substances,
Figure 463730DEST_PATH_IMAGE043
as the first voxel coordinate traversed by the ray, and then by the pair of voxels
Figure 374049DEST_PATH_IMAGE043
Plus one minus one to get its 26 neighbors.
Incremental method for calculating voxel and neighborhood thereof
After the first voxel coordinate index is obtained, other voxels and their neighborhoods are obtained by an addition and subtraction method. First it is determined whether the ray is increasing or decreasing in the x, y, z direction:
Figure 410138DEST_PATH_IMAGE044
（7）
Figure 505133DEST_PATH_IMAGE045
greater than 0, indicating that the voxel traversed by the ray is indexed at
Figure 220148DEST_PATH_IMAGE046
Increasing in direction;
Figure 324370DEST_PATH_IMAGE047
less than 0, the voxel index that the ray passed is at
Figure 214966DEST_PATH_IMAGE048
Decreasing in direction. Then, solve to obtainλAfter the value:
Figure 356228DEST_PATH_IMAGE049
（8）
will be provided withλAnd
Figure 699485DEST_PATH_IMAGE050
a comparison is made to determine the next plane through which the LOR line passes:
Figure 200873DEST_PATH_IMAGE051
（9）
Figure 945975DEST_PATH_IMAGE052
（10）
Figure 382773DEST_PATH_IMAGE053
（11）
wherein
Figure 91622DEST_PATH_IMAGE043
Representing the voxels through which the ray passes,
Figure 272067DEST_PATH_IMAGE054
representing voxels
Figure 871676DEST_PATH_IMAGE043
A neighborhood that is not repeated with the last voxel. It should be noted that, in general, the indices of neighboring voxels in a ray differ by a certain amount
Figure 869588DEST_PATH_IMAGE055
But in special cases the ray will cross from the diagonal of the voxel.
Figure 921857DEST_PATH_IMAGE056
（12）
According to the formulas (8), (9), (10), (11) and (12), the voxel index coordinate to be passed by the ray next can be calculated by using the current voxel index coordinate when the ray passes through
Figure 171573DEST_PATH_IMAGE057
All voxels through which the ray passes can be traversed, wherein,
Figure 235475DEST_PATH_IMAGE058
representing the exit point of the ray through the voxel array.
Distance formula
As can be seen from the above, except that the first voxel coordinate index is obtained by complex calculation, other voxel coordinate indexes are obtained by addition and subtraction on the basis of the first voxel. Therefore, when the distance from the geometric center of the voxel to the ray is calculated, the distance from the first voxel to the ray is calculated by using a point-to-straight line distance formula, and the first voxel indexes coordinates
Figure 279655DEST_PATH_IMAGE043
Can be calculated from equation (6). Let voxel A be the first voxel traversed by the ray PQ, with an index coordinate of
Figure 943854DEST_PATH_IMAGE043
If the index coordinate of the adjacent voxel B of the voxel A is obtained by adding one to or subtracting one from the index coordinate of the voxel A, the index coordinate of the voxel B is set as
Figure 466102DEST_PATH_IMAGE059
. The distance of voxel A to ray PQ is first calculated, where the vector of AQ is
Figure 40303DEST_PATH_IMAGE060
The vector of PQ is
Figure 130750DEST_PATH_IMAGE061
The distance from voxel a to ray PQ is then formulated as follows:
Figure 157612DEST_PATH_IMAGE062
（13）
wherein the content of the first and second substances,
Figure 483551DEST_PATH_IMAGE063
，
Figure 771313DEST_PATH_IMAGE064
representing a vectorbThe die (a) of (b) is,
Figure 157295DEST_PATH_IMAGE065
representing a vector
Figure 546819DEST_PATH_IMAGE066
The mold of (4);
since voxel A and voxel B are adjacent and the vector of AQ is
Figure 676449DEST_PATH_IMAGE060
The vector for BQ can be deduced as
Figure 959663DEST_PATH_IMAGE067
The distance formula of voxel B to ray PQ is as follows:
Figure 641180DEST_PATH_IMAGE068
（14）
wherein the content of the first and second substances,
Figure 377055DEST_PATH_IMAGE069
，
Figure 310376DEST_PATH_IMAGE070
representing a vector
Figure 320533DEST_PATH_IMAGE071
Die (2).
The same can be said for other neighboring voxels. The formula (14) simplifies the point-to-straight line calculation step and increases the calculation speed.
At the time of obtaining the voxel-to-ray distance
Figure 782738DEST_PATH_IMAGE072
Then, if the distance is less thandis（disRepresenting the width of a voxel), then
Figure 396122DEST_PATH_IMAGE073
Is constant if it is greater than or equal todisThen will be
Figure 867554DEST_PATH_IMAGE074
Substituting into the probability formula:
Figure 859781DEST_PATH_IMAGE075
（15）
Figure 899413DEST_PATH_IMAGE073
denotes the firstjThe individual voxel isiThe probability of detection of a bar LOR line,Lis adjusted according to the ratio between the detector modules and the voxels.
As shown in FIG. 2 (a), it is a cross section of a disc prosthesis model with a radius of 37.5mm and hole diameters of 4.8mm,4mm,3.2mm,2.4mm,1.6mm and 1.2mm, respectively. And simulating the resolution prosthesis by using GATE software to acquire the resolution prosthesis data.
Figure 2 (b) is a model of a line source prosthesis and its position in which the gray circle represents the FOV area of the PET device, and a 5mm diameter hollow cylinder filled with a radioactive source is placed 45mm below the center of the FOV. Simulating the line source prosthesis by using GATE software, and acquiring data of the line source prosthesis.
The algorithm of the present embodiment was compared with the siddon algorithm, the center-pixel method, and the orthogonal distance ray tracing algorithm.
1) The Siddon algorithm: the voxel that the LOR line passes through is calculated first, then the intercept at which the LOR line intersects the voxel is calculated, and the intercept represents the probability that the voxel is detected by the LOR line.
2) Center Pixel Method (CPM): the voxels through which the LOR line passes are determined by the siddon algorithm and the probability value for these voxels is assigned to 1.
3) Orthogonal distance ray tracing algorithm (ORRTA): the voxels through which the LOR line passes are determined using the siddon algorithm, then the neighborhood of these voxels is found, and the distance of each voxel to the LOR line is calculated, the smaller the distance, the greater the probability that the voxel is detected by the LOR line.
In the MLEM algorithm, the resolution prosthesis data were tested using the center pixel method, the siddon algorithm, the orthogonal distance algorithm and the inventive algorithm, respectively, and the results are shown in fig. 3 (a) - (d), wherein fig. 3 (e) - (h) are partial image magnifications of fig. 3 (a) - (d). As can be seen from fig. 3, the shape of each point in the image reconstructed by the center pixel method and the siddon algorithm is irregular, while each point in the image reconstructed by the orthogonal distance ray tracing algorithm and the algorithm of the present invention is close to a circle.
The four images in the upper row of fig. 3 were evaluated using no-reference image indicators such as a Natural Image Quality Evaluator (NIQE), a Standard Deviation (SD), a Spatial Frequency (SF), an Entropy (Entropy), and the like, respectively, wherein a smaller NIQE value indicates better image quality, and a larger SD, SF, and Entropy values indicate better image quality, and the results are shown in table 1.
TABLE 1 image quality evaluation results
Method NIQE SD SF Entropy
CPM 24.1471 31.6994 8.1482 2.3048
Siddon 24.0625 30.4521 8.1426 2.2751
ORRTA 24.0520 33.8662 8.2252 2.3025
Ours 23.9745 34.1124 8.2428 2.3072
It can be concluded from fig. 3 and table 1 that the performance of our algorithm model on the resolution prosthesis exceeds that of the center pixel method, the siddon algorithm and the orthogonal distance ray tracing algorithm, which proves that the algorithm of the present invention can reconstruct an image with better quality.
In the MLEM algorithm, the line source prosthesis data were tested using the center pixel method, the siddon algorithm, the orthogonal distance algorithm, and the algorithm of the present invention, respectively, and the results are shown in fig. 4 (a) - (d), in which fig. 4 (e) - (h) are partial image enlargements of fig. 4 (a) - (d).
According to the line source prosthesis model, the line source prosthesis consists of a circle, and as can be seen from fig. 4, the images reconstructed by the central pixel method and the siddon algorithm are irregular circles, the noise around the radioactive source is more, and each point in the images reconstructed by the orthogonal distance ray tracing algorithm and the algorithm of the invention is close to the circle.
And horizontally selecting a ray passing through the line source image, and comparing gray values on the four rays, wherein a peak is a radioactive source. The results are shown in fig. 5 (a), 5 (b), 5 (c), 5 (d), and 5 (e), where zx represents the center pixel method, siddon represents the siddon algorithm, zj represents the orthogonal distance ray tracing algorithm, and outputs represents the algorithm proposed by the present invention. It can be found that the edge of the line source of the image generated by the algorithm of the invention is smoother compared with other algorithms, and the noise of the area outside the line source is relatively less, so that the algorithm of the invention can reconstruct an image with better quality.
Table 2 image quality evaluation results
Method NIQE SD SF Entropy
CPM 20.6126 4.2280 3.0470 0.8060
Siddon 20.8212 4.4116 3.0401 0.8759
ORRTA 19.0518 4.4023 2.9330 0.8754
Ours 18.5210 4.4756 3.0536 0.8947
As can be seen from Table 2, the performance of the algorithm model on the resolution prosthesis exceeds that of the center pixel method, the siddon algorithm and the orthogonal distance ray tracing algorithm, which proves that the algorithm of the invention can reconstruct an image with better quality.
In order to compare the operation speeds of the central pixel method, the siddon algorithm, the orthogonal distance ray tracing algorithm and the algorithm of the invention, the invention compares the speeds of the four algorithms by calculating the time for the ray to pass through the voxel array with the size of 216 × 216 × 163.
As shown in FIG. 6, the time consumed by the orthogonal distance ray tracing algorithm to calculate the ray to pass through the voxel is obviously longer than that consumed by the central pixel method and the siddon algorithm, while the time consumed by the algorithm of the invention is only 0.42 time of that consumed by the orthogonal distance algorithm, so that the calculation time of the orthogonal distance ray tracing algorithm is reduced by about 6/10, and the calculation speed is obviously improved.
Therefore, the algorithm provided by the invention can effectively improve the calculation speed and reduce the calculation time.
In order to verify the effectiveness of the improvement points proposed by the present invention, the present invention performed a series of ablation experiments on a resolution prosthesis model. The experimental setup was as follows: (1) The Ray model uses an orthogonal distance Ray tracing algorithm to calculate the time and imaging quality of the Ray traversing the voxel array. (2) Ray + Incre: and (3) introducing incremental thinking on the basis of a Ray model, and simplifying voxel coordinate calculation. (3) Ray + Incre + Neigh: and (3) optimizing the query of the neighborhood on the basis of (2). (4) Ray + Incre + Neigh + true: and (3) simplifying the calculation of the distance between the point and the straight line. (5) Ray + Incre + Neigh + distance + Pro: i.e. a full version of the model proposed by the present invention. The experimental results are shown in table 3, where Time represents the Time taken for 1000 rays to pass through a voxel array with a size of 216 × 216 × 163, the more the Time is, the slower the speed is, and NIQE represents the imaging quality at different stages of the ablation experiment, and the smaller the value is, the better the imaging quality is. As shown in table 3, it can be found by comparing (1) and (2) that the calculation of the voxel index by the reference increment method can be accelerated. Through the comparison of (2) and (3), the method for directly searching the non-repetitive neighborhood can be found out to reduce the time for removing the duplicate and accelerate the speed. By comparing (3) and (4), the optimization of the distance from the point to the straight line is found to greatly improve the calculation speed of the whole algorithm. Through comparison of (5) with (1), (2), (3) and (4), it can be proved that setting the probability value of the voxel completely covered by the detector to a fixed value can improve the imaging quality when the probability value is calculated.
TABLE 3 ablation test results
Number Model Time NIQE SD SF Entropy
(1) Ray 443ms 24.0520 33.8662 8.2252 2.3025
(2) Ray+Incre 401ms 24.0520 33.8662 8.2252 2.3025
(3) Ray+Incre+Neigh 376ms 24.0520 33.8662 8.2252 2.3025
(4) Ray+Incre+Neigh+tance 188ms 24.0520 33.8662 8.2252 2.3025
(5) Ray+Incre+Neigh+tance+Pro 189ms 23.9745 34.1124 8.2428 2.3072
The invention provides a new algorithm for accelerating the calculation speed and improving the imaging quality based on an orthogonal distance ray tracing algorithm. The algorithm firstly uses an increment method to calculate the coordinate index of a voxel passed by a ray, finds a neighborhood where the voxel is not repeated with the previous voxel, then calculates the distance between the voxel and the ray by using an optimized point-to-straight line distance formula, and finally calculates the probability that each voxel is detected by an LOR line through a new probability formula. Research has shown that compared with the orthogonal distance ray tracing algorithm, the new algorithm reduces the calculation time by about 60% and has certain improvement on the image quality. Therefore, the algorithm ensures the imaging quality and improves the PET imaging speed.
The present invention and its embodiments have been described above schematically, without limitation, and what is shown in the drawings is only one of the embodiments of the present invention, and the actual structure is not limited thereto. Therefore, without departing from the spirit of the present invention, a person of ordinary skill in the art should understand that the present invention shall not be limited to the embodiments and the similar structural modes without creative design.

Claims (2)

1. A system matrix fast iterative computation method for PET image reconstruction is characterized in that: firstly, an LOR line is brought into a siddon algorithm, a first voxel which is passed by the LOR line is calculated, and a voxel neighborhood is obtained; secondly, determining a voxel increment through an LOR line, adding the first voxel passed by the ray with the voxel increment to obtain a voxel to be passed by the ray and searching a neighborhood of the voxel to judge whether the voxel is out of range; if the boundary is not out of range, adding the voxel and the voxel increment to obtain the next voxel, and repeating the step until the boundary is out of range; finally, substituting the calculated neighborhood voxels into a distance algorithm, solving the distance from the voxels to the rays, and substituting the distance into a probability formula to solve a system matrix; wherein:
the specific process of solving the first voxel through which the ray passes by using the siddon algorithm is as follows:
(1) For one
Figure 433857DEST_PATH_IMAGE001
The voxel array of (1), wherein
Figure 571578DEST_PATH_IMAGE002
The number of planes on the x-axis is indicated,
Figure 423996DEST_PATH_IMAGE003
indicating the number of planes on the y-axis,
Figure 647167DEST_PATH_IMAGE004
the number of planes on the z axis is expressed, and the siddon algorithm is utilized to calculate the number of planes on the x axis in the three-dimensional space coordinateN x The positions of the planes parallel to each other are calculated and then calculated on the y-axis and the z-axis respectivelyN y AndN z the position of the planes parallel to each other;
(2) After the position of the voxel array is obtained, siddon computes the entrance and exit points of the ray through the voxel array
Figure 384179DEST_PATH_IMAGE005
And
Figure 251772DEST_PATH_IMAGE006
；
(3) Respectively obtaining the coordinate system in the coordinate system based on the step (1) and the step (2)x,y,zRange of directional plane indices
Figure 150458DEST_PATH_IMAGE007
,
Figure 860925DEST_PATH_IMAGE008
And
Figure 995103DEST_PATH_IMAGE009
；
(4) Respectively calculating the first time of ray passing through the plane index range obtained in the step (3)x,y,zThe value of α at the plane, which represents the ratio of the value of the distance from the ray origin to the ray intersection with a different plane to the distance of the ray in that direction:
Figure 107415DEST_PATH_IMAGE010
wherein the content of the first and second substances,λrepresenting a first plane through which the ray passesαA value;
Figure 52369DEST_PATH_IMAGE011
representing the second on the x-axisi min A plane and ai max The set of ratios of the distances from the planes to the ray origin to the ray x-axis direction distances,
Figure 250132DEST_PATH_IMAGE012
、
Figure 328946DEST_PATH_IMAGE013
the correspondences represent sets in the y-axis and z-axis,
Figure 420399DEST_PATH_IMAGE014
representing the first ray in the direction of the x-axisThe value of a for each of the planes,
Figure 395308DEST_PATH_IMAGE015
representing the value of alpha for a ray traversing a first plane in the y-axis direction,
Figure 345947DEST_PATH_IMAGE016
an alpha value representing the ray through a first plane in the z-axis direction;
(5) Based on a first plane traversed by the rayαValue, calculating the index coordinate of the first voxel that the ray passes through
Figure 572660DEST_PATH_IMAGE017
；
After the first voxel coordinate index is found, it is first determined whether the ray is increasing or decreasing in the x, y, z direction, and then solvedλValue of willλAnd with
Figure 659565DEST_PATH_IMAGE018
Comparing, determining the next plane to be crossed by the LOR line, and repeating the steps when the LOR line passes through
Figure 70954DEST_PATH_IMAGE019
All voxels through which the ray passes can be traversed at this time, wherein,
Figure 102364DEST_PATH_IMAGE020
an exit point representing a ray through the voxel array;
when calculating the distance from the geometric center of the voxel to the ray, firstly, the distance from the first voxel A to the ray PQ is calculated by using a point-to-straight line distance formula, wherein the vector of AQ is
Figure 522981DEST_PATH_IMAGE021
The vector of PQ is
Figure 77109DEST_PATH_IMAGE022
Distance of voxel A from ray PQThe formula is as follows:
Figure 924980DEST_PATH_IMAGE023
wherein, the first and the second end of the pipe are connected with each other,
Figure 319052DEST_PATH_IMAGE024
，
Figure 136835DEST_PATH_IMAGE025
representing a vectorbThe die (a) of (b) is,
Figure 198332DEST_PATH_IMAGE026
representing a vector
Figure 826891DEST_PATH_IMAGE027
The die of (2);
for the neighboring voxel B of voxel A, the vector of its BQ is
Figure 708259DEST_PATH_IMAGE028
The distance formula of voxel B to ray PQ is as follows:
Figure 205099DEST_PATH_IMAGE029
wherein the content of the first and second substances,
Figure 245737DEST_PATH_IMAGE030
，
Figure 435409DEST_PATH_IMAGE031
representing a vector
Figure 804074DEST_PATH_IMAGE032
The mold of (4);
the situation of other neighboring voxels can be derived in the same way;
at the time of obtaining the voxel-to-ray distance
Figure 979971DEST_PATH_IMAGE033
Then, if the distance is less thandisThen, then
Figure 484902DEST_PATH_IMAGE034
Is a constant if greater than or equal todisThen will be
Figure 845476DEST_PATH_IMAGE035
Substituting into the probability formula:
Figure 560491DEST_PATH_IMAGE036
Figure 664714DEST_PATH_IMAGE034
is shown asjThe individual voxel isiThe probability of detection of a bar LOR line,diswhich represents the width of the voxel(s),Lis adjusted according to the ratio between the detector modules and the voxels.
2. The fast iterative computation method of system matrix for PET image reconstruction as recited in claim 1, wherein: in the case of a ray crossing from a diagonal of a voxel, it willλThrough with raysx,y,zComparing the minimum distance ratios of the planes, and if equal, expressing the minimum distance ratios asαThe values are incremented in the respective directions, let:
Figure 165096DEST_PATH_IMAGE037
Figure 696572DEST_PATH_IMAGE038
Figure 774249DEST_PATH_IMAGE039
wherein, the first and the second end of the pipe are connected with each other,dx、dyanddzrespectively representx、y、zThe distance between the planes in the direction of the plane,
Figure 541217DEST_PATH_IMAGE040
、
Figure 20740DEST_PATH_IMAGE041
and
Figure 723117DEST_PATH_IMAGE042
respectively represent LOR lines onx、y、zThe increment between the planes in the direction is,
Figure 160527DEST_PATH_IMAGE043
is the coordinates of the origin of the ray,
Figure 606551DEST_PATH_IMAGE044
is the end point coordinate of the ray,
Figure 330794DEST_PATH_IMAGE045
representing the alpha value of the ray after increment in the direction of the x axis,
Figure 204072DEST_PATH_IMAGE046
representing the alpha value of the ray after increment in the y-axis direction,
Figure 990762DEST_PATH_IMAGE047
representing the alpha value of the ray after increment in the z-axis direction.
CN202211252651.0A 2022-10-13 2022-10-13 System matrix fast iterative computation method for PET image reconstruction Active CN115330900B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211252651.0A CN115330900B (en) 2022-10-13 2022-10-13 System matrix fast iterative computation method for PET image reconstruction

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211252651.0A CN115330900B (en) 2022-10-13 2022-10-13 System matrix fast iterative computation method for PET image reconstruction

Publications (2)

Publication Number Publication Date
CN115330900A CN115330900A (en) 2022-11-11
CN115330900B true CN115330900B (en) 2022-12-27

Family

ID=83914994

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211252651.0A Active CN115330900B (en) 2022-10-13 2022-10-13 System matrix fast iterative computation method for PET image reconstruction

Country Status (1)

Country Link
CN (1) CN115330900B (en)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104116519A (en) * 2014-07-30 2014-10-29 西安电子科技大学 Rotating plate double-panel PET system based on polygon view symmetry and imaging method thereof
CN104968275A (en) * 2013-01-31 2015-10-07 株式会社东芝 System optics in back projection and/or forward projection for model-based iterative reconstruction
CN105931280A (en) * 2016-03-29 2016-09-07 中北大学 GPU-based fast three-dimensional CT (computed tomography) iterative reconstruction method
CN109949411A (en) * 2019-03-22 2019-06-28 电子科技大学 A kind of image rebuilding method based on three-dimensional weighted filtering back projection and statistics iteration

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8314796B2 (en) * 2006-02-24 2012-11-20 The Board Of Trustees Of The Leland Stanford Junior University Method of reconstructing a tomographic image using a graphics processing unit
WO2021114306A1 (en) * 2019-12-14 2021-06-17 深圳先进技术研究院 High-precision pet reconstruction method and apparatus based on gpu acceleration

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104968275A (en) * 2013-01-31 2015-10-07 株式会社东芝 System optics in back projection and/or forward projection for model-based iterative reconstruction
CN104116519A (en) * 2014-07-30 2014-10-29 西安电子科技大学 Rotating plate double-panel PET system based on polygon view symmetry and imaging method thereof
CN105931280A (en) * 2016-03-29 2016-09-07 中北大学 GPU-based fast three-dimensional CT (computed tomography) iterative reconstruction method
CN109949411A (en) * 2019-03-22 2019-06-28 电子科技大学 A kind of image rebuilding method based on three-dimensional weighted filtering back projection and statistics iteration

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Fast calculation of the exact radiological path for a three-dimensional CT array;Robert L. Siddon;《Medical Physics》;19850331;第252-255页 *
四平板PET的系统矩阵计算与快速重建;程健;《中国优秀硕士学位论文全文数据库 信息科技辑》;20200215;正文第4-6、9-17页 *

Also Published As

Publication number Publication date
CN115330900A (en) 2022-11-11

Similar Documents

Publication Publication Date Title
Zhang et al. Cross-based local stereo matching using orthogonal integral images
CN104700451B (en) Point cloud registering method based on iterative closest point algorithm
CN107016672B (en) Reconstruction method and device of medical scanning image and medical imaging system
CN111311704A (en) Image reconstruction method and device, computer equipment and storage medium
CN102222318A (en) Grid artifact detection and suppression method and device
CN110060316B (en) Ring artifact correction method for multi-region segmentation in CT reconstruction
CN105931280A (en) GPU-based fast three-dimensional CT (computed tomography) iterative reconstruction method
CN110415280B (en) Remote sensing image and building vector registration method and system under multitask CNN model
CN111105452A (en) High-low resolution fusion stereo matching method based on binocular vision
CN114612318B (en) Three-dimensional modeling method, system and equipment based on cultural relic CT image contour line
CN115294527A (en) Subway tunnel damage detection method based on computer vision
CN110335322B (en) Road recognition method and road recognition device based on image
US20220030132A1 (en) Method and device for correcting a scanned image and image scanning system
CN115330900B (en) System matrix fast iterative computation method for PET image reconstruction
CN110268441A (en) The method for obtaining the 3D model data of multiple components of object
CN106611173B (en) Crystal center position diagram generation method
CN105844690A (en) Rapid three-dimensional CT iteration reconstruction system based on GPU
WO2023065780A1 (en) Convolution algorithm-based image processing method and chip
CN105701847A (en) Algebraic reconstruction method of improved weight coefficient matrix
CN105222730B (en) A kind of industry CT physical dimension measuring method based on image restoration
CN114663416A (en) Computer vision-based copper bar cold shut defect detection method and system
CN114943761A (en) Method and device for extracting center of light stripe of central line structure of FPGA (field programmable Gate array)
Kornilov et al. Deep Neural Networks for Ring Artifacts Segmentation and Corrections in Fragments of CT Images
CN110728669B (en) Video mosaic detection method
CN113729747A (en) Spherical metal marked cone beam CT metal artifact removing system and method

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant