CN116740205A - Source track differential image reconstruction method for linear CT scanning of ray source - Google Patents
Source track differential image reconstruction method for linear CT scanning of ray source Download PDFInfo
- Publication number
- CN116740205A CN116740205A CN202310666181.0A CN202310666181A CN116740205A CN 116740205 A CN116740205 A CN 116740205A CN 202310666181 A CN202310666181 A CN 202310666181A CN 116740205 A CN116740205 A CN 116740205A
- Authority
- CN
- China
- Prior art keywords
- image
- hilbert
- dimensional
- projection
- reconstructed
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 55
- 238000002591 computed tomography Methods 0.000 title claims abstract description 46
- 230000009466 transformation Effects 0.000 claims abstract description 17
- 230000001351 cycling effect Effects 0.000 claims abstract description 6
- 230000009191 jumping Effects 0.000 claims abstract description 6
- 238000004364 calculation method Methods 0.000 claims description 23
- 230000005855 radiation Effects 0.000 claims description 21
- 239000011159 matrix material Substances 0.000 claims description 20
- 238000013519 translation Methods 0.000 claims description 18
- 238000012545 processing Methods 0.000 claims description 3
- 230000001131 transforming effect Effects 0.000 claims description 3
- 238000004422 calculation algorithm Methods 0.000 description 7
- 238000003384 imaging method Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 238000002474 experimental method Methods 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
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
- G06T11/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2200/00—Indexing scheme for image data processing or generation, in general
- G06T2200/08—Indexing scheme for image data processing or generation, in general involving all processing steps from image acquisition to 3D model generation
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Algebra (AREA)
- Computer Graphics (AREA)
- Geometry (AREA)
- Software Systems (AREA)
- Image Processing (AREA)
Abstract
A source track differential image reconstruction method facing to ray source straight line CT scanning comprises a two-dimensional method and a three-dimensional method. The method comprises the following steps: s1, initializing parameters i=1, the number of straight line scanning segments T and the rotation angle interval delta theta, and obtaining a zero space of an image to be reconstructedS2, obtaining projection of linear scanning of the ith section of ray source; s3, pre-weighting projection data; s4, differentiating the weighted projection along the direction of the source track; s5, carrying out weighted back projection on the differential projection to obtain an i-th Hilbert image; s6, carrying out limited Hilbert inverse transformation on the Hilbert image along the ith section of ray source track, and reconstructing an ith section of limited angle imageS7. let i=i+1,s8, if i is less than or equal to T, rotating the linear CT scanning system by an angle (i-1) delta theta, jumping to S2, and sequentially cycling until i is more than T to finish an imageIs performed in the reconstruction of (a). The image reconstruction method is oriented to the linear CT scanning track, can directly avoid artifacts caused by truncated projections, and reconstructs high-quality images.
Description
Technical Field
The invention belongs to the field of CT image reconstruction, and particularly relates to a source track differential image reconstruction method for linear CT scanning of a ray source.
Background
In recent years, in order to meet the higher resolution requirements, a plurality of different CT scanning modes have emerged. Liu Fenglin and Yu Haijun, etc. propose a structure of expanding imaging field of view for linear CT scanning of a radiation source, in which the detector is fixed in imaging geometry, the radiation source is linearly translated in a direction parallel to the detector row, the object to be measured is placed in the center of a turntable, and iterative image reconstruction algorithms based on TV minimization are designed [1,2]. The linear CT scanning structure of the ray source can realize the expansion of imaging vision through controlling the source translation, is easy to realize high-precision control, has great development potential in the technical field of CT scanning, and is sequentially researched along with a series of image reconstruction algorithms suitable for the scanning model.
The iterative image algorithm based on TV minimization is firstly provided, and although the iterative image algorithm can reconstruct an image without artifacts, the reconstruction time is long, the calculation cost is high, and the actual engineering requirements [1,2] are difficult to meet. Although the traditional FBP (filtered back projection) analysis type reconstruction algorithm can perform high-efficiency reconstruction, in the linear CT scanning of the enlarged imaging view field, data truncation exists in the acquired projection data, and then the reconstructed image obtained by the filtered back projection shows serious truncation artifacts [2-5]. In order to avoid truncation artifacts by adopting an FBP type reconstruction algorithm, the prior art adopts a virtual projection mode to realize equivalent global non-truncation projection, and then image reconstruction is carried out [3-5]. Although the method of converting projection into virtual projection by recombining projection can avoid truncation artifact, the reconstruction method is still not simple and direct enough.
Reference is made to:
[1] liu Fenglin, yu Haijun, li Lei, tan Chuandong a novel large field-of-view linear scanning CT system and image reconstruction method [ P ]. Chongqing city: CN111839568A,2020-10-30.
[2]H.Yu,L.Li,C.Tan,F.Liu,R.Zhou,X-ray source translation based computed tomography(STCT),Optics Express,29(2021)19743-19758.R.Clackdoyle,F.Noo,A large class of inversion formulae for the 2D Radon transform of functions of compact support,Inverse Problems,20(2004)1281.
[3] Gossypol Wen Jie, yu Haijun, chen Jie, et al source linear scan computed tomography analytical reconstruction based on derivative-hilbert transform-backprojection [ J ]. Optics, 2022,42 (11): 292-303.
[4] Li Lei, yu Haijun, tan Chuandong, duan Xiaojiao, liu Fenglin. Analytical reconstruction algorithm for radiation source translation scanning CT [ J ]. Instructions on instruments and meters, 2022,43 (02): 187-195.DOI:10.19650/j.cnki.cjsi.J2108157.
[5]Yu H,Ni S,Chen J,et al.Analytical reconstruction algorithm for multiple source-translation computed tomography(mSTCT)[J].Applied Mathematical Modelling,2023,117:251-266.
Disclosure of Invention
In order to overcome the defects in the prior art, the invention aims to provide a source track differential image reconstruction method for linear CT scanning of a ray source, which directly avoids truncation artifacts and reconstructs high-quality artifact-free images.
In order to achieve the above object, the present invention provides two-dimensional and three-dimensional image reconstruction methods. The two-dimensional method adopts the following technical scheme:
s1, initializing parameters i=1, the number of linear scanning segments T and the rotation angle interval delta theta, and reconstructing a zero space of a two-dimensional image to be reconstructed
S2, acquiring projection of linear scanning of the ith section of ray source
S3, projection of acquisitionPre-weighting to obtain weighted projection +.>
S4, weighting projectionPerforming differential operation along the direction of the translation track lambda of the ray source;
s5, carrying out weighted back projection operation on the differential projection to obtain an ith Hilbert image
S6, carrying out limited Hilbert inverse transformation on the ith section Hilbert image along the linear translation track direction of the ith section ray source, and reconstructing the ith section limited angle image
S7, executing i=i+1, and finishing superposition of images:
s8, if i is less than or equal to T, rotating the linear CT scanning system by an angle (i-1) delta theta, jumping to S2, and sequentially cycling until i is more than T to finish an imageIs performed in the reconstruction of (a).
Preferably, for the source trajectory differential two-dimensional image reconstruction method of the ray source oriented straight line CT scan, in step S3, the projectionsPre-weighting results in a weighted projection +.>The calculation formula is that
wherein ,θi Is the direction angle theta of the linear scanning track of the ray source i = (i-1) ·Δθ; lambda represents the coordinates of the source on a straight line trajectory, lambda epsilon-s, s]The method comprises the steps of carrying out a first treatment on the surface of the u represents the coordinates of the detector, and u E [ -d, d]D represents half the detector length; l represents the distance from the source to the centre of rotation in the direction perpendicular to the detector; h represents the distance from the center of rotation to the center of the detector;is a redundant weight factor.
Preferably, for the source trajectory differential two-dimensional image reconstruction method of the ray source oriented straight line CT scan, in step S4, the weighted projectionDifferential operation along the direction of the translational trajectory lambda of the source to obtain a differential projectionThe calculation formula is that
wherein ,the differential operator along the lambda direction is represented and realized by finite difference operation, the data in the projection adopts center difference, and the data of the projection boundary adopts single-side difference; lambda (lambda) * Represents crossing the point to be reconstructed->Coordinate values of the rays of (a) on the source translation trajectory,
wherein L represents a point to be reconstructedPerpendicular distance to the source trajectory, l= -xsin θ i +ycosθ i +l; h represents the vertical distance of the point x to be reconstructed from the detector, h=xsin θ i -ycosθ i +h。
Preferably, for the source trajectory differential two-dimensional image reconstruction method of the ray source oriented straight line CT scan, in step S5, the differential projectionWeighted back projection operation is carried out to obtain i-th Hilbert image +.>The calculation formula is as follows,
in the formula ,representing a point to be reconstructed; η (eta) i Representing the Hilbert image for the i th segment>A direction angle for performing a limited Hilbert inverse transform; 1/H 2 Representing the backprojection weighting factor; the back projection operation is performed in a matrix larger than the image to be reconstructed, i.e. p is added outside each edge of the matrix space of the image to be reconstructed 0 Zero in number =0.5×i, where I is the image row or column number to be reconstructed.
Preferably, for the source trajectory differential two-dimensional image reconstruction method of the linear CT scan facing the radiation source, in step S6, the pair of i-th segment Hilbert imagesThe limited Hilbert inverse transformation is carried out along the direction of the linear translation track of the ith section ray source, and the ith section limited angle image can be reconstructed>The calculation formula is that
wherein ,y1 Representing the one-dimensional Hilbert transform direction, y 1 ∈[L y +ε y ,U y -ε y], wherein [Ly ,U y ]Representing the finite interval, ε, of the Hilbert transform y Is a small positive number;representing reconstructing an i-th segment limited angle image +.>The unknown constant that needs to be calculated is determined by finding the known +.>Is used as the mean value after being subjected to the finite Hilbert transformationHilbert image of section i->A specific implementation of the limited inverse hilbert transform comprises the steps of:
1) Image of Hilbert section iDot-2 pi and rotate around the center of the image by-eta i Obtaining a rotated Hilbert image +.> wherein ηi ＝θ i Pi/2 (i=1, 2,., T), angle value positive, image counter-clockwise rotation; the angle is negative and the image rotates clockwise;
2) Calculating a weight matrix W mat I.e.
3) Searching for rotated Hilbert imagesThe values of which rows are all 0, and an array R with the elements being row numbers is obtained i ；
4) To rotating Hilbert imagesWeighting W mat And performing Hilbert transform along the matrix array direction to obtain image +.>
5) According to R obtained in step S3 i Numerical value, slave imageMiddle indexElement is extracted, the mean value is calculated along the column direction and is inverted to obtain +.>
6) Execution of
7) Will rotate the imageRotation angle eta about the center of the image i To the original state, get ∈>
8) An intermediate image region to be reconstructed is extracted,
the three-dimensional method provided by the invention adopts the following technical scheme:
s1, initializing parameters i=1, the number of straight line scanning segments T and the rotation angle interval delta theta, and reconstructing a zero space of a three-dimensional image
S2, acquiring cone beam projection of linear scanning of an ith section of ray source
S3, projecting collected cone beamsPre-weighting to obtain weighted cone beam projection +.>
S4, projecting the weighted cone beamsPerforming differential operation along the direction of the radiation source track lambda;
s5, carrying out weighted back projection operation on the differential cone beam projection to obtain an ith section of three-dimensional Hilbert image
S6, carrying out limited Hilbert inverse transformation on the ith section of three-dimensional Hilbert image layer by layer axially along the direction of linear translation track of the ith section of ray source, and reconstructing the ith section of limited angle three-dimensional image
S7, executing i=i+1, and finishing superposition of three-dimensional images:
s8, if i is less than or equal to T, enabling the linear CT scanning system to rotate by an angle (i-1) delta theta, jumping to S2, and sequentially cycling until i is more than T, and completing the three-dimensional imageIs performed in the reconstruction of (a).
Preferably, for the source trajectory differential three-dimensional image reconstruction method of the ray source-oriented straight line CT scan, in step S3, the cone beam projectionPre-weighting to obtain weighted cone beam projections +.>The calculation formula is that
wherein ,θi Is the direction angle theta of the linear scanning track of the ray source i = (i-1) ·Δθ; lambda represents a rayCoordinates of the source on a straight line trajectory, lambda E-s, s]The method comprises the steps of carrying out a first treatment on the surface of the u represents the line-wise coordinates of the panel detector, and u E [ -d, d]D represents half of the line length of the panel detector; l represents the distance from the source to the center of rotation in the direction perpendicular to the panel detector; h represents the distance from the center of rotation to the center of the panel detector;is a redundant weight factor.
Preferably, for the source trajectory differential three-dimensional image reconstruction method of the ray source oriented straight line CT scan, in step S4, the weighted cone beam projectionsDifferential operation is carried out along the lambda direction of the ray source track to obtain differential cone beam projectionThe calculation formula is that
wherein ,the differential operator along the lambda direction is represented and realized by finite difference operation, the data in the projection adopts center difference, and the data of the projection boundary adopts single-side difference; lambda (lambda) * Represents crossing the point to be reconstructed->The row-wise coordinate of the ray of (2) on the ray source trajectory λ, v representing the x-ray passing through the point to be reconstructed>The column-wise coordinates of the rays of (a) on the panel detector,
wherein L represents a point to be reconstructedPerpendicular distance to the source trajectory, l= -xsin θ i +ycosθ i +l; h represents the point to be reconstructed +.>Vertical distance to panel detector, h=xsin θ i -ycosθ i +h。
Preferably, for the source trajectory differential three-dimensional image reconstruction method of the ray source-oriented straight line CT scan, in step S5, the differential cone beam projectionsPerforming weighted back projection operation to obtain an ith section of three-dimensional Hilbert imageThe calculation formula is that
in the formula ,representing a point to be reconstructed; η (eta) i Representing +.>A direction angle for performing a limited Hilbert inverse transform; 1/H 2 Weighting factors for the backprojection; the back projection operation is performed in a matrix larger than the three-dimensional image to be reconstructed, namely, each time in the matrix space of the three-dimensional image to be reconstructedP is increased outside the edge 0 Zero in number =0.5×i, where I is the number of voxels in a certain dimension of the image to be reconstructed.
Preferably, for the source trajectory differential three-dimensional image reconstruction method of the linear CT scan facing the radiation source, in step S6, the method further comprises the step ofAxially carrying out the inverse Hilbert transformation layer by layer along the linear translation track direction of the ith section of ray source, namely, z by z k (k=1, 2,., I) inverse transforming the two-dimensional Hilbert image layer by layer, thereby reconstructing an I-th limited-angle three-dimensional image +.>The inverse transform of the two-dimensional Hilbert image is also along x j (j=1, 2,., I) one-dimensional processing is performed piece by piece, and the calculation formula is
wherein ,y1 Representing the one-dimensional Hilbert transform direction, y 1 ∈[L y +ε y ,U y -ε y], wherein [Ly ,U y ]Representing the finite interval, ε, of the Hilbert transform y Is a small positive number;representing reconstruction->The unknown constant that needs to be calculated is determined by finding the known +.>Is used as C by taking the mean value of the position points after being subjected to the Hilbert transformation yi The method comprises the steps of carrying out a first treatment on the surface of the Three-dimensional Hilbert image for the ith segment->Is subjected to a limited Hilbert inverse transformation layer by layer in the axial direction, and comprises the following steps:
1) Three-dimensional Hilbert image of the ith sectionDot-2 pi and rotate around the central axis of the three-dimensional image by-eta i Angle, obtain rotated three-dimensional Hilbert image +.> wherein ηi ＝θ i Pi/2 (i=1, 2,., T), the angle value being positive, the three-dimensional image rotating anticlockwise; the angle is negative, and the three-dimensional image rotates clockwise;
2) Calculate one-dimensional weight sequence W mat I.e.
3) for (k=1; k is less than or equal to I; k++) (i.e., into the for loop body):
3.1 From a rotated three-dimensional Hilbert imageIndexing out a kth layer two-dimensional Hilbert image
3.2 Searching for two-dimensional Hilbert imagesThe values of which rows are all 0, and an array R with the elements being row numbers is obtained i ；
3.3 For two-dimensional Hilbert imageWeighting W mat And performing Hilbert transform along the matrix array direction to obtain image +.>
3.4 According to R obtained in step S3 i Numerical value, slave imageThe elements are extracted from the middle rope, and the average value is calculated along the column direction and is inverted to obtain +.>
3.5 Execution of (1)
3.6 A two-dimensional image to be rotated)Rotation angle eta about the center of the image i To the original state, the reconstructed k layer two-dimensional +.>
4) Extracting intermediate three-dimensional image regions to be reconstructed, i.e.
The invention has the beneficial effects that:
compared with the existing reconstruction algorithm for the linear CT scanning of the ray source, the method can utilize truncated projection data to carry out accurate reconstruction, directly avoid truncation artifacts, and reconstruct a two-dimensional or three-dimensional CT image almost without artifacts.
Drawings
Fig. 1 is a schematic flow chart of the present invention.
Fig. 2 is a schematic diagram of a linear CT two-dimensional scan trajectory of a radiation source according to the present invention.
FIG. 3 is a flow chart of two-dimensional image reconstruction using a Shepp-Logan phantom as a measured object in the invention.
Fig. 4 is a two-dimensional reconstructed image of different projection numbers of a measured object using a fouild phantom in accordance with the present invention.
Fig. 5 is a schematic diagram of a ray source linear cone beam CT three-dimensional scan trajectory for which the present invention is directed.
Fig. 6 is a schematic representation of the finite inverse Hilbert transform of a three-dimensional Hilbert image of the present invention axially layer by layer along the direction of the ray source trajectory.
Fig. 7 is a flow chart of three-dimensional CT image reconstruction of the present invention.
FIG. 8 shows the image reconstruction results under different projection numbers by using a three-dimensional Shepp-Logan phantom as a measured object.
Detailed Description
Embodiments of the present invention will be described in detail below with reference to the accompanying drawings and examples.
A source track differential image reconstruction method facing to ray source straight line CT scanning is shown in figure 1, and the two-dimensional flow comprises the following steps:
s1, initializing parameters i=1, the number of linear scanning segments T and the rotation angle interval delta theta, and reconstructing a zero space of a two-dimensional image to be reconstructed
S2, acquiring projection of linear scanning of the ith section of ray source
S3, projection of acquisitionPre-weighting to obtain weighted projection +.>
S4, weighting projectionPerforming differential operation along the direction of the translation track lambda of the ray source;
s5, carrying out back projection operation on the differential projection to obtain an ith Hilbert image
S6, carrying out limited Hilbert inverse transformation on the ith section Hilbert image along the linear translation track direction of the ith section ray source, and reconstructing the ith section limited angle image
S7, executing i=i+1, and finishing superposition of images:
s8, if i is less than or equal to T, rotating the linear CT scanning system by an angle (i-1) delta theta, jumping to S2, and sequentially cycling until i is more than T to finish an imageIs performed in the reconstruction of (a).
Preferably, for the source trajectory differential two-dimensional image reconstruction method of the ray source oriented straight line CT scan, in step S3, the projectionsPre-weighting results in a weighted projection +.>The calculation formula is that
wherein ,θi Is the direction angle theta of the linear scanning track of the ray source i = (i-1) ·Δθ; lambda represents the coordinates of the source on a straight line trajectory, lambda epsilon-s, s]The method comprises the steps of carrying out a first treatment on the surface of the u represents the coordinates of the detector, and u E [ -d, d]D represents half the detector length; l represents the distance from the source to the centre of rotation in the direction perpendicular to the detector;h represents the distance from the center of rotation to the center of the detector;is a redundant weight factor.
Preferably, for the source trajectory differential two-dimensional image reconstruction method of the ray source oriented straight line CT scan, in step S4, the weighted projectionDifferential operation along the direction of the translational trajectory lambda of the source to obtain a differential projectionThe calculation formula is that
wherein ,the differential operator along the lambda direction is represented and realized by finite difference operation, the data in the projection adopts center difference, and the data of the projection boundary adopts single-side difference; lambda (lambda) * Represents crossing the point to be reconstructed->Coordinate values of the rays of (a) on the source translation trajectory,
wherein L represents a point to be reconstructedPerpendicular distance to the source trajectory, l= -xsin θ i +ycosθ i +l; h represents the vertical distance of the point x to be reconstructed from the detector, h=xsin θ i -ycosθ i +h。
Preferably, for the source trajectory differential two-dimensional image reconstruction method of the ray source oriented straight line CT scan, in step S5, the differential projectionWeighted back projection operation is carried out to obtain i-th Hilbert image +.>The calculation formula is as follows,
in the formula ,representing a point to be reconstructed; η (eta) i Representing the Hilbert image for the i th segment>A direction angle for performing a limited Hilbert inverse transform; 1/H 2 Representing the backprojection weighting factor; the back projection operation is performed in a matrix larger than the image to be reconstructed, i.e. p is added outside each edge of the matrix space of the image to be reconstructed 0 Zero in number =0.5×i, where I is the row or column number of the matrix of images to be reconstructed.
Preferably, for the source trajectory differential two-dimensional image reconstruction method of the linear CT scan facing the radiation source, in step S6, the pair of i-th segment Hilbert imagesThe limited Hilbert inverse transformation is carried out along the direction of the linear translation track of the ith section ray source, and the ith section limited angle image can be reconstructed>The calculation formula is that
wherein ,y1 Representing the one-dimensional Hilbert transform direction, y 1 ∈[L y +ε y ,U y -ε y], wherein [Ly ,U y ]Representing the finite interval, ε, of the Hilbert transform y Is a small positive number;representing reconstructing an i-th segment limited angle image +.>The unknown constant that needs to be calculated is determined by finding the known +.>Is used as the mean value after being subjected to the finite Hilbert transformationHilbert image of section i->A specific implementation of the limited inverse hilbert transform comprises the steps of:
1) Image of Hilbert section iDot-2 pi and rotate around the center of the image by-eta i Obtaining a rotated Hilbert image +.> wherein ηi ＝θ i Pi/2 (i=1, 2,., T), angle value positive, image counter-clockwise rotation; the angle is negative and the image rotates clockwise;
2) Calculating a weight matrix W mat I.e.
3) Searching for rotated Hilbert imagesThe values of which rows are all 0, and an array R with the elements being row numbers is obtained i ；
4) To rotating Hilbert imagesWeighting W mat And performing Hilbert transform along the matrix array direction to obtain image +.>
5) According to R obtained in step S3 i Numerical value, slave imageThe elements are extracted from the middle rope, and the average value is calculated along the column direction and is inverted to obtain +.>
6) Execution of
7) Will rotate the imageRotation angle eta about the center of the image i To the original state, get ∈>
8) An intermediate image region to be reconstructed is extracted,
fig. 2 shows a schematic diagram of the first two scans in a linear CT scan trajectory of a radiation source, corresponding to the 1 st and 2 nd radiation source linear scans of the present invention.
Fig. 3 shows the image reconstruction flow of the invention when the Shepp-Logan phantom is taken as a measured object, which proves the effectiveness of the method and can reconstruct an image without artifacts.
Further, using the FORBILD phantom, the number of pixels was 512×512, and the size was 8.4mm×8.4mm. The image reconstruction effect of the invention is evaluated by a numerical experiment. The projection number N of each section of ray source straight line CT scanning is set to be 50, 100, 200, 400, 600 and 800 respectively, the front projection of the FORBILD body model is calculated in a simulation mode and is reconstructed by adopting the method, and as can be known from the image reconstruction result shown in fig. 4, the method can reconstruct a high-quality image without artifacts by utilizing the projection data truncated in the ray source straight line CT scanning process, and the image reconstruction quality is increased along with the increase of the projection number.
Table 1 numerical simulation test parameters
A source track differential image reconstruction method facing to ray source straight line CT scanning is shown in figure 1, and the three-dimensional method flow comprises the following steps:
s1, initializing parameters i=1, the number of straight line scanning segments T and the rotation angle interval delta theta, and reconstructing a zero space of a three-dimensional image
S2, acquiring cone beam projection of linear scanning of an ith section of ray source
S3, projecting collected cone beamsPre-weighting to obtain weighted cone beam projection +.>
S4, projecting the weighted cone beamsPerforming differential operation along the direction of the radiation source track lambda;
s5, carrying out weighted back projection operation on the differential cone beam projection to obtain an ith section of three-dimensional Hilbert image
S6, carrying out limited Hilbert inverse transformation on the ith section of three-dimensional Hilbert image layer by layer axially along the direction of linear translation track of the ith section of ray source, and reconstructing the ith section of limited angle three-dimensional image
S7, executing i=i+1, and finishing superposition of three-dimensional images:
s8, if i is less than or equal to T, enabling the linear CT scanning system to rotate by an angle (i-1) delta theta, jumping to S2, and sequentially cycling until i is more than T, and completing the three-dimensional imageIs performed in the reconstruction of (a).
Preferably, for the source trajectory differential three-dimensional image reconstruction method of the ray source-oriented straight line CT scan, in step S3, the cone beam projectionPre-weighting to obtain weighted cone beam projections +.>The calculation formula is +.>
wherein ,θi Is the direction angle theta of the linear scanning track of the ray source i = (i-1) ·Δθ; lambda represents the coordinates of the source on a straight line trajectory, lambda epsilon-s, s]The method comprises the steps of carrying out a first treatment on the surface of the u represents the line-wise coordinates of the panel detector, and u E [ -d, d]D represents half of the line length of the panel detector; l represents the distance from the source to the center of rotation in the direction perpendicular to the panel detector; h represents the distance from the center of rotation to the center of the panel detector;is a redundant weight factor.
Preferably, for the source trajectory differential three-dimensional image reconstruction method of the ray source oriented straight line CT scan, in step S4, the weighted cone beam projectionsDifferential operation is carried out along the lambda direction of the ray source track to obtain differential cone beam projectionThe calculation formula is that
wherein ,the differential operator along the lambda direction is represented and realized by finite difference operation, the data in the projection adopts center difference, and the data of the projection boundary adopts single-side difference; lambda (lambda) * Represents crossing the point to be reconstructed->And v denotes the column coordinate on the panel detector of the radiation passing through the point x to be reconstructed,
wherein L represents a point to be reconstructedPerpendicular distance to the source trajectory, l= -xsin θ i +ycosθ i +l; h represents the point to be reconstructed +.>Vertical distance to panel detector, h=xsin θ i -ycosθ i +h。
Preferably, for the source trajectory differential three-dimensional image reconstruction method of the ray source-oriented straight line CT scan, in step S5, the differential cone beam projectionsPerforming weighted back projection operation to obtain an ith section of three-dimensional Hilbert imageThe calculation formula is that
in the formula ,representing a point to be reconstructed; η (eta) i Representing +.>A direction angle for performing a limited Hilbert inverse transform; 1/H 2 Weighting factors for the backprojection; the back projection operation is performed in a matrix larger than the three-dimensional image to be reconstructed,i.e. adding p outside each edge of the matrix space of the three-dimensional image to be reconstructed 0 Zero in number =0.5×i, where I is the number of voxels in a certain dimension of the image to be reconstructed.
Preferably, for the source trajectory differential three-dimensional image reconstruction method of the linear CT scan facing the radiation source, in step S6, the method further comprises the step ofAxially carrying out the inverse Hilbert transformation layer by layer along the linear translation track direction of the ith section of ray source, namely, z by z k (k=1, 2,., I) inverse transforming the two-dimensional Hilbert image layer by layer, thereby reconstructing an I-th limited-angle three-dimensional image +.>The inverse transform of the two-dimensional Hilbert image is also along x j (j=1, 2,., I) one-dimensional processing is performed piece by piece, and the calculation formula is
wherein ,y1 Representing the one-dimensional Hilbert transform direction, y 1 ∈[L y +ε y ,U y -ε y], wherein [Ly ,U y ]Representing the finite interval, ε, of the Hilbert transform y Is a small positive number;representing reconstruction->The unknown constant that needs to be calculated is determined by finding the known +.>Is used as +.>Three-dimensional Hilbert image for the ith segment->Is subjected to a limited Hilbert inverse transformation layer by layer in the axial direction, and comprises the following steps:
1) Three-dimensional Hilbert image of the ith sectionDot-2 pi and rotate around the central axis of the three-dimensional image by-eta i Angle, obtain rotated three-dimensional Hilbert image +.> wherein ηi ＝θ i Pi/2 (i=1, 2,., T), the angle value being positive, the three-dimensional image rotating anticlockwise; the angle is negative, and the three-dimensional image rotates clockwise;
2) Calculate one-dimensional weight sequence W mat I.e.
3) for (k=1; k is less than or equal to I; k++) (i.e., into the for loop body):
3.1 From a rotated three-dimensional Hilbert imageIndexing out a kth layer two-dimensional Hilbert image
3.2 Searching for two-dimensional Hilbert imagesThe values of which rows are all 0, and an array R with the elements being row numbers is obtained i ；
3.3 For two-dimensional Hilbert imageWeighting W mat And performing Hilbert transform along the matrix array direction to obtain image +.>
3.4 According to R obtained in step S3 i Numerical value, slave imageThe elements are extracted from the middle rope, and the average value is calculated along the column direction and is inverted to obtain +.>
3.5 Execution of (1)
3.6 A two-dimensional image to be rotated)Rotation angle eta about the center of the image i To the original state, the reconstructed k layer two-dimensional +.>
4) Extracting intermediate three-dimensional image regions to be reconstructed, i.e.
Fig. 5 shows a schematic view of a linear cone beam CT scan trajectory of a source, and fig. 6 shows a three-dimensional Hilbert image for the i-th segmentAnd (3) carrying out finite Hilbert inverse transformation axially layer by layer along the linear translation track direction of the ith section of ray source.
Fig. 7 shows a basic flow of three-dimensional CT image reconstruction according to the present invention, which can prove the effectiveness of the present method, and can reconstruct a high quality near artifact free three-dimensional CT image.
Further, a three-dimensional Shepp-Logan phantom is used, the number of pixels was 512×512, and the size was 8.4mm×8.4mm. The image reconstruction effect of the invention is evaluated by quantitative numerical experiments. The result of the reconstructed image is shown in fig. 8, the condition of different projection numbers N under each section of ray source straight line CT scan is set, the cone beam orthographic projection of the three-dimensional Shepp-Logan phantom is simulated and calculated, the reconstruction is carried out by adopting the invention, the quality of the reconstructed image is observed, N is set as 200 and 1600 respectively, and the parameter settings of the numerical experiment are shown in table 2. From the reconstruction results shown in fig. 8, it can be seen that the present invention can reconstruct a high-quality three-dimensional CT image almost free of artifacts when the number of projections is large.
Table 2 numerical simulation test parameters
The invention is not limited to the specific embodiments described above, which are intended to be illustrative only and not limiting; those skilled in the art, having the benefit of this disclosure, may make numerous forms without departing from the spirit of the invention and the scope of the claims which follow.
Claims (7)
1. A source track differential image reconstruction method facing to ray source straight line CT scanning is characterized in that the two-dimensional or three-dimensional method comprises the following steps:
s1, initializing parameters i=1, the number of linear scanning segments T and the rotation angle interval delta theta, and reconstructing a zero space of a two-dimensional image to be reconstructed
S2, acquiring projection of linear scanning of the ith section of ray sourceOr->
S3, projection of acquisitionOr->Pre-weighting to obtain weighted projection +.>Or (b)
S4, weighting projectionOr->Performing differential operation along the direction of the translation track lambda of the ray source;
s5, carrying out weighted back projection operation on the differential projection to obtain an ith Hilbert image
S6, carrying out limited Hilbert inverse transformation on the ith section Hilbert image along the linear translation track direction of the ith section ray source, and reconstructing the ith section limited angle image
S7, executing i=i+1, and finishing superposition of images:
s8, if i is less than or equal to T, rotating the linear CT scanning system by an angle (i-1) delta theta, jumping to S2, and sequentially cycling until i is more than T to finish an imageIs performed in the reconstruction of (a).
2. The method of reconstructing a source trajectory differential image for a linear CT scan of a radiation source of claim 1, wherein in step S3, said projectionOr->Pre-weighting results in a weighted projection +.>Or (b)The calculation formula is that
Or (b)
wherein ,θi Is the direction angle theta of the linear scanning track of the ray source i = (i-1) ·Δθ; lambda represents the coordinates of the source on a straight line trajectory, lambda epsilon-s, s]The method comprises the steps of carrying out a first treatment on the surface of the u represents the coordinates of the detector, and u E [ -d, d]D represents half the detector length; l represents the distance from the source to the centre of rotation in the direction perpendicular to the detector; h represents the distance from the center of rotation to the center of the detector;is a redundant weight factor.
3. The method of source trajectory differential two-dimensional image reconstruction for a linear CT scan of a radiation source according to claim 1, wherein in step S4, the weighted projectionsDifferential operation along the direction of the translational trajectory lambda of the radiation source to obtain a differential projection +.>The calculation formula is that
wherein ,the differential operator along the lambda direction is represented and realized by finite difference operation, the data in the projection adopts center difference, and the data of the projection boundary adopts single-side difference; lambda (lambda) * Representing the coordinate values of the rays passing through the point x to be reconstructed on the translational trajectory of the source,
where L represents the perpendicular distance of the point to be reconstructed x to the ray source trajectory, l= -xsin θ i +ycosθ i +l; h represents the point to be reconstructedVertical distance to detector, h=xsin θ i -ycosθ i +h。
4. Root of Chinese characterThe method of reconstructing a source trajectory differential two-dimensional image for a linear CT scan of a source of radiation of claim 1, wherein in step S5, said differential projectionsOr->Weighted back projection operation is carried out to obtain i-th Hilbert image +.>The calculation formula is as follows,
in the formula ,or->Representing a point to be reconstructed; η (eta) i Representing the Hilbert image for the i th segment>A direction angle for performing a limited Hilbert inverse transform; 1/H 2 Representing the backprojection weighting factor; the back projection operation is performed in a matrix larger than the image to be reconstructed, i.e. p is added outside each edge of the matrix space of the image to be reconstructed 0 Zero in number =0.5×i, where I is the image row or column number to be reconstructed.
5. The method of claim 1, wherein in step S6, the pair of i-th Hilbert imagesThe limited Hilbert inverse transformation is carried out along the direction of the linear translation track of the ith section ray source, and the ith section limited angle image can be reconstructed>The calculation formula is that
wherein ,y1 Representing the one-dimensional Hilbert transform direction, y 1 ∈[L y +ε y ,U y -ε y], wherein [Ly ,U y ]Representing the finite interval, ε, of the Hilbert transform y Is a small positive number;representing reconstructing an i-th segment limited angle image +.>The unknown constant that needs to be calculated is determined by finding the known +.>Is used as +.>Hilbert image of section i->A specific implementation of the limited inverse hilbert transform comprises the steps of:
1) Image of Hilbert section iDot-2 pi and windImage center rotation-eta i Obtaining a rotated Hilbert image +.> wherein ηi ＝θ i Pi/2 (i=1, 2,., T), angle value positive, image counter-clockwise rotation; the angle is negative and the image rotates clockwise;
2) Calculating a weight matrix W mat I.e.
3) Searching for rotated Hilbert imagesThe values of which rows are all 0, and an array R with the elements being row numbers is obtained i ；
4) To rotating Hilbert imagesWeighting W mat And performing Hilbert transform along the matrix array direction to obtain image +.>
5) According to R obtained in step S3 i Numerical value, slave imageThe elements are extracted from the middle cable, the average value of the elements is calculated along the column direction and is inverted to obtain C yi ；
6) Execution of
7) Will rotate the imageRotation angle eta about the center of the image i To the original state, get ∈>
8) An intermediate image region to be reconstructed is extracted,
6. the method of source trajectory differential three-dimensional image reconstruction for a linear CT scan of a radiation source according to claim 1, wherein in step S4, the weighted cone beam projectionsDifferential operation is performed along the direction of the radiation source trajectory lambda to obtain a differential cone beam projection +.>The calculation formula is that
wherein ,the differential operator along the lambda direction is represented and realized by finite difference operation, the data in the projection adopts center difference, and the data of the projection boundary adopts single-side difference; lambda (lambda) * Representing the row-wise coordinate of the radiation passing through the point x to be reconstructed on the radiation source trajectory λ, v representing the x-ray passing through the point x to be reconstructed>The column-wise coordinates of the rays of (a) on the panel detector,
where L represents the perpendicular distance of the point to be reconstructed x to the ray source trajectory, l= -xsin θ i +ycosθ i +l; h represents the point to be reconstructedVertical distance to panel detector, h=xsin θ i -ycosθ i +h。
7. The method for reconstructing a source trajectory differential three-dimensional image for linear CT scanning with a radiation source according to claim 1, wherein in step S6, said pair of i-th segment three-dimensional Hilbert imagesAxially carrying out the inverse Hilbert transformation layer by layer along the linear translation track direction of the ith section of ray source, namely, z by z k (k=1, 2,., I) inverse transforming the two-dimensional Hilbert image layer by layer, thereby reconstructing an I-th limited-angle three-dimensional image +.>The inverse transform of the two-dimensional Hilbert image is also along x j (j=1, 2,., I) one-dimensional processing is performed piece by piece, and the calculation formula is
wherein ,y1 Representing the one-dimensional Hilbert transform direction, y 1 ∈[L y +ε y ,U y -ε y], wherein [Ly ,U y ]Representing the finite interval, ε, of the Hilbert transform y Is a smallA positive number;representing reconstruction->The unknown constant that needs to be calculated is determined by finding the known +.>Is used as +.>Three-dimensional Hilbert image for the ith segment->Is subjected to a limited Hilbert inverse transformation layer by layer in the axial direction, and comprises the following steps:
1) Three-dimensional Hilbert image of the ith sectionDot-2 pi and rotate around the central axis of the three-dimensional image by-eta i Angle, obtain rotated three-dimensional Hilbert image +.> wherein ηi ＝θ i Pi/2 (i=1, 2,., T), the angle value being positive, the three-dimensional image rotating anticlockwise; the angle is negative, and the three-dimensional image rotates clockwise;
2) Calculate one-dimensional weight sequence W mat I.e.
3) for (k=1; k is less than or equal to I; k++) (i.e., into the for loop body):
3.1 From rotationThree-dimensional Hilbert imageTwo-dimensional Hilbert image of the kth layer is indexed out +.>
3.2 Searching for two-dimensional Hilbert imagesThe values of which rows are all 0, and an array R with the elements being row numbers is obtained i ；
3.3 For two-dimensional Hilbert imageWeighting W mat And performing Hilbert transform along the matrix array direction to obtain image +.>
3.4 According to R obtained in step S3 i Numerical value, slave imageThe elements are extracted from the middle rope, and the average value is calculated along the column direction and is inverted to obtain +.>
3.5 Execution of (1)
3.6 A two-dimensional image to be rotated)Rotation angle eta about the center of the image i To the original state, the reconstructed k layer two-dimensional +.>
4) Extracting intermediate three-dimensional image regions to be reconstructed, i.e.
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310302887 | 2023-03-24 | ||
| CN2023103028879 | 2023-03-24 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116740205A true CN116740205A (en) | 2023-09-12 |
Family
ID=87918017
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310666181.0A Pending CN116740205A (en) | 2023-03-24 | 2023-06-06 | Source track differential image reconstruction method for linear CT scanning of ray source |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116740205A (en) |
-
2023
- 2023-06-06 CN CN202310666181.0A patent/CN116740205A/en active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107764846B (en) | Orthogonal linear scanning CL imaging system and analysis method | |
| CN107796834B (en) | Orthogonal electronic linear scanning CL imaging system and method | |
| CN107328798B (en) | Novel ICL system and implementation method | |
| JPH11326243A (en) | Method for operating ct imaging apparatus and scan and data collection apparatus | |
| CN110057847B (en) | TR (transmitter-receiver) tomography projection rearrangement method and device | |
| CN102456227A (en) | Reconstruction method and device for CT (computerized tomography) image | |
| CN109146987B (en) | GPU-based rapid cone beam computed tomography reconstruction method | |
| US8948337B2 (en) | Computed tomography image reconstruction | |
| CN100348158C (en) | Rapid progressive three-dimensional reconstructing method of CT image from direct volume rendering | |
| CN105118039A (en) | Method and system for reconstructing cone beam CT image | |
| CN111839568B (en) | Novel large-view-field linear scanning CT system and image reconstruction method | |
| CN116630460A (en) | Detector line differential high-quality image reconstruction method for source linear scanning track | |
| Ametova et al. | Software-based compensation of instrument misalignments for X-ray computed tomography dimensional metrology | |
| CN107233105B (en) | Correction method and correction system for CT image reconstruction | |
| CN106097411B (en) | CT machine image rebuilding method and high resolution ct scanning machine | |
| CN113533392B (en) | Combined scanning CL imaging method | |
| US7272205B2 (en) | Methods, apparatus, and software to facilitate computing the elements of a forward projection matrix | |
| CN116843779A (en) | Linear scanning detector differential BPF reconstructed image sparse artifact correction method | |
| CN110264536B (en) | Method for calculating high-low resolution projection relation in parallel beam ultra-resolution reconstruction | |
| CN102599887B (en) | Optical projection tomography method based on helical scanning track | |
| CN116188615A (en) | Sparse angle CT reconstruction method based on sine domain and image domain | |
| CN116740205A (en) | Source track differential image reconstruction method for linear CT scanning of ray source | |
| CN105266839B (en) | A kind of different big visual field CT imaging methods of three sources Circular test radius | |
| CN114820927A (en) | Industrial CT three-dimensional image reconstruction method based on data rearrangement and conjugate rays | |
| Shih et al. | Fast algorithm for X-ray cone-beam microtomography |
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 |