CN109102553B - Polar coordinate system matrix calculation method and device in two-dimensional reconstruction algorithm - Google Patents

Polar coordinate system matrix calculation method and device in two-dimensional reconstruction algorithm Download PDF

Info

Publication number
CN109102553B
CN109102553B CN201810921211.7A CN201810921211A CN109102553B CN 109102553 B CN109102553 B CN 109102553B CN 201810921211 A CN201810921211 A CN 201810921211A CN 109102553 B CN109102553 B CN 109102553B
Authority
CN
China
Prior art keywords
polar coordinate
projection
system matrix
grid
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.)
Active
Application number
CN201810921211.7A
Other languages
Chinese (zh)
Other versions
CN109102553A (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.)
Space Engineering University
Original Assignee
Space Engineering University
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 Space Engineering University filed Critical Space Engineering University
Publication of CN109102553A publication Critical patent/CN109102553A/en
Application granted granted Critical
Publication of CN109102553B publication Critical patent/CN109102553B/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
    • 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
    • G06T2210/00—Indexing scheme for image generation or computer graphics
    • G06T2210/41—Medical
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Length Measuring Devices By Optical Means (AREA)

Abstract

The invention discloses a polar coordinate system matrix calculation method and a device in a two-dimensional reconstruction algorithm, wherein the method comprises the following steps: step S100: establishing a polar coordinate projection model of a two-dimensional field to be reconstructed;step S200: calculating the overlapping area of the projection beam and the polar coordinate grid under the ith projection angle as the system matrix element value amnIn the subscript, m is the serial number of the detector element, and n is the serial number of the grid; step S300: traversing and calculating all the overlapping areas under the projection angle to obtain a system matrix A ═ a of the two-dimensional field to be reconstructedmnAnd adjusting the positions of the element values of the system matrix in the system matrix according to the change of the projection angle to obtain the system matrix under other projection angles without repeated calculation. The method can improve the precision of an algebraic iterative two-dimensional reconstruction algorithm, simultaneously avoids recalculating a system matrix when the projection angle is changed, and improves the calculation efficiency. The invention also provides a device used in the method.

Description

Polar coordinate system matrix calculation method and device in two-dimensional reconstruction algorithm
Technical Field
The invention relates to a polar coordinate system matrix calculation method and device in a two-dimensional reconstruction algorithm, and belongs to the field of optical measurement.
Background
Computer Tomography (CT) based on the principle of spectral absorption is well-established and applied to the field of medical imaging, and can provide a clear two-dimensional reconstructed image for medical diagnosis. Due to the universality of the method, the method can be combined with Tunable Diode Laser Absorption Spectroscopy (TDLAS) technology to realize two-dimensional reconstruction of the temperature and concentration of the engine combustion flow field. The CT reconstruction algorithm mainly includes a Filtered Back Projection (FBP) method and an algebraic iterative Technique (ART). The ART method has very good reconstruction accuracy and noise level compared to the FBP method, and does not require complete projection data.
The ART method is severely limited in its operation speed. The ART algorithm needs to construct a huge system matrix to participate in reconstruction operation, the volume of the system matrix can reach T bytes, and the system matrix cannot be stored in signal processing equipment. The currently adopted method is to calculate a system matrix once in the iterative process of each reconstruction operation, which greatly increases the operation amount and the operation time of the ART algorithm.
The existing tomography generally adopts a rectangular coordinate system to perform reconstruction view field meshing. The grid division mode has two problems, namely when the projection angle is changed, the system matrix must be recalculated in each iteration, and the calculation efficiency is low; secondly, when the size of the detector is considered, numerical calculation errors are introduced into the existing segmentation calculation mode, and the calculation accuracy is reduced.
Disclosure of Invention
According to one aspect of the invention, a polar coordinate system matrix calculation method in a two-dimensional reconstruction algorithm is provided. The method can improve the precision of an algebraic iterative two-dimensional reconstruction algorithm, simultaneously avoids recalculating a system matrix, and improves the calculation efficiency.
A polar coordinate system matrix calculation method in a two-dimensional reconstruction algorithm comprises the following steps:
step S100: establishing a polar coordinate projection model of the two-dimensional field to be reconstructed in an area surrounded by the point light source and the detection element array, and dividing the polar coordinate projection model at equal intervals to construct a polar coordinate grid;
step S200: calculating the overlapping area of the projection light beam generated by the point light source and the detection element and the polar coordinate grid under the ith projection angle as the system matrix element value amnIn the subscript, m is the serial number of the detector element, and n is the serial number of the grid;
step S300: at the projection angle of ith, the projection beam traverses the polar grid and the detection cellArray to obtain multiple matrix element values amnA plurality of values of elements a of said system matrixmnThe ith system matrix A is obtained by aggregationmn}；
Step S400: taking the projection angle of the (i + g) th system matrix A as { a ] according to the rotation invariance of the polar coordinate grid divisionmnEach of said system matrix element values a inmnAdjusting the position of the system to obtain the i + g system matrix;
step S500: and repeating the step S400 when g is g +1, and circulating until g is N, N is a natural positive integer to obtain N (i + g) th system matrixes, and collecting the N (i + g) th system matrixes to obtain the system matrix of the two-dimensional field to be reconstructed.
Optionally, the step S100 includes the following steps:
step S110: constructing an isosceles triangle projection area by taking the point light source as a vertex and the detection cell array as a bottom edge;
step S120: establishing a polar coordinate system by taking the circle center of the two-dimensional field to be reconstructed as an origin, and determining the geometric positions of the point light source and the detection array element in the polar coordinate system;
step S130: dividing the two-dimensional field to be reconstructed into N at equal intervals by taking the circle center of the two-dimensional field to be reconstructed as a starting point and taking delta r as an intervalrA concentric circle with said two-dimensional field of view azimuth to be reconstructed
Figure BDA0001764271590000021
As a starting point, in
Figure BDA0001764271590000022
Equally-spaced partitioning the two-dimensional field of view to be reconstructed for an interval
Figure BDA0001764271590000023
A bar ray intersecting the concentric circles to construct the polar grid.
Optionally, the step S200 includes the following steps:
step S210: polar coordinatesGrid mesh
Figure BDA0001764271590000024
Is a concentric circle rjConcentric circle rj+1Ray of radiation
Figure BDA0001764271590000025
And ray
Figure BDA0001764271590000026
Enclosing, calculating the projection beam and the concentric circle rjThe concentric circle rj+1The ray
Figure BDA0001764271590000027
And said ray
Figure BDA0001764271590000028
The coordinates of the intersection points;
step S220: calculating the projection beam and the ray according to the intersection point coordinates
Figure BDA0001764271590000029
The ray
Figure BDA00017642715900000210
The concentric circles rjOverlap area S of enclosed sector aream,i,jCalculating the projection beam and the ray according to the intersection point coordinates
Figure BDA00017642715900000211
The ray
Figure BDA00017642715900000212
The concentric circles rj+1Overlap area S of enclosed sector aream,i,j+1；
Step S230: the system matrix element value amnThe overlapping area Sm,i,j+1-said overlapping area Sm,i,jThe grid number N in the subscript is equal to Nr*(j-1)+i，NrThe number of the concentric circles.
Optionally, the two-dimensional field of view to be reconstructed is an axially symmetric region, and the projection beam at the ith projection angle in step S200 scans any symmetric region in the axially symmetric region.
Optionally, an angle difference between the ith projection angle and the (i + g) th projection angle is
Figure BDA0001764271590000031
Optionally, the two-dimensional field to be reconstructed is circular and is arranged in the isosceles triangle projection area, and the side of the two-dimensional field to be reconstructed is tangent to the waist line of the isosceles triangle projection area.
Optionally, the detection cell array is a linear array or a fan-shaped array.
The invention also provides a polar coordinate system matrix calculation device in a two-dimensional reconstruction algorithm, which comprises:
the polar coordinate grid construction module is used for establishing a polar coordinate projection model of the two-dimensional field to be reconstructed in an area surrounded by the point light source and the detection cell array, and dividing the polar coordinate projection model at equal intervals to construct a polar coordinate grid;
a system matrix element value calculation unit: used for calculating the overlapping area of the projection light beam generated by the point light source and the detection element and the polar coordinate grid under the ith projection angle as the system matrix element value amnIn the subscript, m is the serial number of the detector element, and n is the serial number of the grid;
a first system matrix calculation unit: the projection light beam traverses the polar coordinate grid and the detection element array under the ith projection angle to obtain a plurality of system matrix element values amnA plurality of values of elements a of said system matrixmnThe ith system matrix A is obtained by aggregationmn}；
A second system matrix calculation unit: the system matrix A is used for taking the projection angle of the i + g th system matrix A as { a ] according to the rotation invariance of the polar coordinate grid divisionmnEach of the system matrix elements in (1)Elemental value amnAdjusting the position of the system to obtain the i + g system matrix;
a circulation module: and the system matrix is used for returning to the second system matrix calculation unit after g is g +1, and the process is circulated until g is N, N is a natural positive integer to obtain N (i + g) th system matrices, and the N (i + g) th system matrices are collected to obtain the system matrix of the two-dimensional field to be reconstructed.
Optionally, the polar grid construction module includes:
a projection region construction module: the system is used for constructing an isosceles triangle projection area by taking the point light source as a vertex and the detection primitive array as a bottom edge;
a coordinate module: the system is used for establishing a polar coordinate system by taking the circle center of the two-dimensional field to be reconstructed as an origin, and determining the geometric positions of the point light source and the detection array element in the polar coordinate system;
a segmentation module: the method is used for dividing the two-dimensional field of view to be reconstructed into N at equal intervals by taking the circle center of the two-dimensional field of view to be reconstructed as a starting point and taking delta r as an intervalrA concentric circle with said two-dimensional field of view azimuth to be reconstructed
Figure BDA0001764271590000041
As a starting point, in
Figure BDA0001764271590000042
Equally-spaced partitioning the two-dimensional field of view to be reconstructed for an interval
Figure BDA0001764271590000043
A bar ray intersecting the concentric circles to construct the polar grid.
Optionally, the system matrix element value calculating unit includes:
an intersection point calculation module: polar grid
Figure BDA0001764271590000044
Is a concentric circle rjConcentric circle rj+1Ray of radiation
Figure BDA0001764271590000045
And ray
Figure BDA0001764271590000046
Enclosing, calculating the projection beam and the concentric circle rjThe concentric circle rj+1The ray
Figure BDA0001764271590000047
And said ray
Figure BDA0001764271590000048
The coordinates of the intersection points;
an overlap area calculation module: for calculating the projection beam and the ray according to the intersection point coordinates
Figure BDA0001764271590000049
The ray
Figure BDA00017642715900000410
The concentric circles rjOverlap area S of enclosed sector aream,i,jCalculating the projection beam and the ray according to the intersection point coordinates
Figure BDA00017642715900000411
The ray
Figure BDA00017642715900000412
The concentric circles rj+1Overlap area S of enclosed sector aream,i,j+1；
An element value calculation module: value of matrix element a for the systemmnThe overlapping area Sm,i,j+1-said overlapping area Sm,i,jThe grid number N in the subscript is equal to Nr*(j-1)+i，NrThe number of the concentric circles.
The beneficial effects of the invention include but are not limited to:
(1) according to the polar coordinate system matrix calculation method and device in the two-dimensional reconstruction algorithm, the existing method proves the advantages of the polar coordinate grid division method in reducing the system matrix and improving the reconstruction speed, but all polar coordinate grids adopt a numerical calculation method, and the error problem of the system matrix in the numerical calculation process is ignored. The invention realizes the analytic calculation method of the system matrix on the basis of the grid division of the polar coordinates, obtains the system matrix with absolute calculation precision, and further improves the reconstruction precision on the basis of improving the operation speed of an algebraic iterative algorithm.
(2) According to the polar coordinate system matrix calculation method and device in the two-dimensional reconstruction algorithm, the polar coordinate system is adopted to perform grid division on the reconstruction view field to obtain the view field grid with rotational symmetry, the reuse of the single projection system matrix is guaranteed, each element value of the system matrix is calculated by using an analytic method, and numerical errors are eliminated. The method can be used for two-dimensional reconstruction precision of medical images or combustion flow fields.
(3) The method and the device for calculating the polar coordinate system matrix in the two-dimensional reconstruction algorithm are suitable for analytic calculation of the field of view projection system matrix in the polar coordinate system, are applied to an algebraic iterative algorithm of a fan-shaped emission beam, and obtain a system matrix with absolute calculation accuracy by considering the width size of a detection element.
Drawings
FIG. 1 is a schematic block diagram of a flow chart of a polar coordinate system matrix calculation method in a two-dimensional reconstruction algorithm provided by the present invention;
FIG. 2 is a schematic flow chart of a polar coordinate system matrix calculation method in a two-dimensional reconstruction algorithm according to the present invention;
FIG. 3 is a schematic structural diagram of a polar coordinate system matrix calculation device in a two-dimensional reconstruction algorithm provided by the present invention;
FIG. 4 is a schematic view of a fan beam projection model based on polar grid division as used in the preferred embodiment of the present invention;
FIG. 5 is a schematic diagram of the system matrix element parsing computation geometry in the preferred embodiment of the present invention;
FIG. 6 is a schematic view of the overlapping of the projection beam with the polar grid in a preferred embodiment of the present invention, wherein (a) shows a first overlapping pattern, (b) shows a second overlapping pattern, (c) shows a third overlapping pattern, and (d) shows a fourth overlapping pattern;
fig. 7 is a schematic diagram of curves of residual errors obtained by an analytic calculation method and a numerical calculation method in a comparative example according to the change of the number of numerical filling points and the calculation time in the preferred embodiment of the present invention, in which a curve with a small hollow dot is the preferred embodiment of the present invention, and a curve with a small hollow diamond is the comparative example.
Illustration of the drawings:
s: a point light source;
o: a polar origin;
D1，D2: the distance from the origin of the polar coordinates to the point light source and the detection array;
dm: the mth detection element;
amn: the element of the mth row and the nth column in the coefficient matrix;
Figure BDA0001764271590000051
dividing the azimuth angle of the ith ray of the polar coordinate;
ri: dividing the jth concentric circle radius of the polar coordinate;
Figure BDA0001764271590000052
dividing the azimuth angle interval of the polar coordinates;
Δ r: dividing concentric circle intervals of polar coordinates;
L1,L2: azimuth angle
Figure BDA0001764271590000053
And
Figure BDA0001764271590000054
the corresponding ray;
Lb，Lr: detection element dmA side line of the projection beam formed with the light source;
A1，A2: side line LbAnd the concentric circle rjThe intersection point of (a);
B1，B2: side line LrAnd the concentric circle rjThe intersection point of (a);
C1，C2: concentric circles rjAnd ray L1And L2The intersection point of (a);
E1，F1: concentric circles rj+1And ray LbAnd LrThe intersection point of (a);
P1，P2: side line LbAnd ray L1And L2The intersection point of (a);
Q1，Q2: side line LrAnd ray L1And L2The intersection point of (a).
Detailed Description
The present invention will be described in detail with reference to examples, but the present invention is not limited to these examples.
Referring to fig. 1-2, the invention provides a polar coordinate system matrix calculation method in a two-dimensional reconstruction algorithm, comprising the following steps of S100: establishing a polar coordinate projection model of the two-dimensional field to be reconstructed in an area surrounded by the point light source and the detection element array, and dividing the polar coordinate projection model at equal intervals to construct a polar coordinate grid;
step S200: calculating the overlapping area of the projection light beam generated by the point light source and the detection element and the polar coordinate grid under the ith projection angle as the system matrix element value amnIn the subscript, m is the serial number of the detector element, and n is the serial number of the grid;
step S300: under the ith projection angle, the projection light beam traverses the polar coordinate grid and the detection cell array to obtain a plurality of system matrix element values amnA plurality of values of elements a of said system matrixmnThe ith system matrix A is obtained by aggregationmn}；
Step S400: taking the projection angle of the (i + g) th system matrix A as { a ] according to the rotation invariance of the polar coordinate grid divisionmnEach of said system matrix element values a inmnIn a position ofAdjusting rows to obtain the i + g system matrix;
step S500: and repeating the step S400 when g is g +1, and circulating until g is N, N is a natural positive integer to obtain N (i + g) th system matrixes, and collecting the N (i + g) th system matrixes to obtain the system matrix of the two-dimensional field to be reconstructed.
In the method, the positions of the system matrix element values in the system matrix can be adjusted according to the change of the projection angle of the system matrix under other projection angles, and the system matrix can be obtained without repeated calculation. Therefore, the calculation amount is effectively reduced, and the calculation efficiency is improved. And in the step of constructing the polar coordinate grid, the division is respectively carried out along the radius and the azimuth direction of the polar coordinate projection model.
In the conventional rectangular coordinate grid division mode, because the area of each grid is the same, the grid area is used to normalize the overlapped area and then the normalized overlapped area is used as a system matrix element value. In the invention, the grid of the polar field of view is increased along with the increase of the radius, and the normalization processing is not suitable, so the overlapping area is used as the element value of the system matrix.
In the actual tomography process, the projection angle can be changed continuously, so that the finally obtained system matrix is a three-dimensional matrix, and the number of elements can reach 106The above is very bulky.
After the polar coordinate grid is divided according to the mode, when the projection angle is at intervals
Figure BDA0001764271590000071
When scanning is performed, the overlapping mode of the projection beam and the grid is consistent at each projection angle, which is the rotation invariance of the polar grid. Therefore, the method provided by the invention can be used for calculating the system matrix of only one projection angle, and then regularly changing the positions of elements in the system matrix by utilizing the coordinate relation to obtain the system matrix of other projection angles, so that the number of the elements in the system matrix is reduced by 2 to 3 orders of magnitude, the calculation efficiency is improved, and the calculation amount is reduced.
Preferably, the step S100 includes the steps of:
as in fig. 2, step S110: constructing an isosceles triangle projection area by taking the point light source as a vertex and the detection primitive array as a bottom edge; the detection cell array comprises M detectors, and the M detectors surround into corresponding shapes according to the array shape requirement.
Step S120: establishing a polar coordinate system by taking the circle center of the two-dimensional field to be reconstructed as an origin, and determining the geometric positions of the point light source and the detection array element in the polar coordinate system;
step S130: dividing the two-dimensional field to be reconstructed into N at equal intervals by taking the circle center of the two-dimensional field to be reconstructed as a starting point and taking delta r as an intervalrA concentric circle with said two-dimensional field of view azimuth to be reconstructed
Figure BDA0001764271590000072
As a starting point, in
Figure BDA0001764271590000073
Equally-spaced partitioning the two-dimensional field of view to be reconstructed for an interval
Figure BDA0001764271590000074
A bar ray intersecting the concentric circles to construct the polar grid. The polar coordinates of each grid can be used
Figure BDA0001764271590000075
Denotes that j ∈ [1, N ] in the subscriptr],
Figure BDA0001764271590000076
There are a total of N grids, each grid having a grid width,
Figure BDA0001764271590000077
the invention provides a method for determining the overlap area of the projection beam and the polar coordinate grid by the element value contained in the system matrix to be calculated. The projection light beam is a small triangle formed by the point light source and the detector. Referring to FIGS. 4 to 5, amnFor the overlap region of the projection beam and the polar grid, M is the [1, M ] in the subscript]Numbering the detector elements, N belongs to [1, N ]]Is aPolar grid numbering.
Preferably, the two-dimensional field to be reconstructed is circular and is arranged in the isosceles triangle projection area, and the side of the two-dimensional field to be reconstructed is tangent to the waist line of the isosceles triangle projection area.
Optionally, the detection cell array is a linear array or a fan-shaped array.
Preferably, the step S200 includes the steps of:
referring to fig. 5, step S210: polar grid
Figure BDA0001764271590000081
Is a concentric circle rjConcentric circle rj+1Ray of radiation
Figure BDA0001764271590000082
And ray
Figure BDA0001764271590000083
Enclosing, calculating the projection beam and the concentric circle rjThe concentric circle rj+1The ray
Figure BDA0001764271590000084
And said ray
Figure BDA0001764271590000085
The coordinates of the intersection points;
step S220: calculating the projection beam and the ray according to the intersection point coordinates
Figure BDA0001764271590000086
The ray
Figure BDA0001764271590000087
The concentric circles rjOverlap area S of enclosed sector aream,i,jCalculating the projection beam and the ray according to the intersection point coordinates
Figure BDA0001764271590000088
The ray
Figure BDA0001764271590000089
The concentric circles rj+1Overlap area S of enclosed sector aream,i,j+1；
Step S230: the system matrix element value amnThe overlapping area Sm,i,j+1-said overlapping area Sm,i,jThe grid number N in the subscript is equal to Nr*(j-1)+i，NrThe number of the concentric circles.
The element value in the system matrix to be calculated is equal to the overlapping area of the fan-shaped projection beam formed by the point light source and each detector and the polar coordinate grid.
Preferably, the two-dimensional field of view to be reconstructed is an axially symmetric region, and the projection beam at the ith projection angle in step S200 scans any symmetric region in the axially symmetric region.
Preferably, the difference between the projection angle of the ith and the projection angle of the (i + g) th is
Figure BDA00017642715900000810
The following describes in detail a method and an apparatus for calculating a polar coordinate system matrix in a two-dimensional reconstruction algorithm according to the present invention with reference to specific examples.
Example 1 and comparative example treated subjects: a typical polar meshing example. The field of view area is a circle with a radius of 90.5mm, equally divided into 200 concentric circles at intervals in the radial direction with a pitch of 0.4525mm, and then divided into 360 equal parts in the azimuth direction in 1 ° steps, making 72000 polar grids. The distance from the light source to the circle center is 800mm, and the distance from the circle center to the detection is 200 mm. The number of detector elements is 512, and the spacing is 0.5 mm. The polar coordinate system matrix was calculated according to the above analysis method with a calculation time of 13s (Intel Core i5 CPU,4G DDR 3).
Embodiment 1 matrix calculation method for polar coordinate system in two-dimensional reconstruction algorithm
Step S100: polar coordinate projection model establishment for two-dimensional field of view to be reconstructed
The fan-shaped light beam projection model based on polar coordinate grid division is shown in fig. 4. The origin of the polar coordinates is located at the center of the field of view, and the radius of the field of view is defined as R. In practice, the effective field of view is not entirely circular, but needs to be contained within this circular field of view. The detector array structure can be linear or fan-shaped, the invention is explained by taking a linear array as an example, the detector array is set to comprise M elements, and the length of each element is d. The distance from the light source S to the detector array is D, the perpendicular line passes through the origin of the coordinate system, and the distances from the origin to the light source S and the detector array are respectively D1And D2Obviously, D ═ D1+D2。
Referring to fig. 4, the point light source S and each detector element form a projection fan passing through the two-dimensional field of view to be reconstructed. Dividing a field of view to be reconstructed into N in the direction of radius r at intervals of delta rrConcentric circles, then the azimuth angle
Figure BDA0001764271590000091
To be provided with
Figure BDA0001764271590000092
For dividing intervals into
Figure BDA00017642715900000915
The strip rays, rays and concentric circles divide the field of view into a radial grid structure, and the polar coordinates of each grid can be used
Figure BDA0001764271590000093
Denotes that j ∈ [1, N ] in the subscriptr],
Figure BDA0001764271590000094
There are a total of N grids, each grid having a grid width,
Figure BDA0001764271590000095
the system matrix element a shown in fig. 4mnExpressed as point source and m-th detector cell dmFormed projection beam and nThe area of the overlap region of the grids, where M ∈ [1, M ∈ ]]，n＝Nr*(j-1)+i,n∈[1,N]. The system matrix dimension of the resulting single projection is M × N.
Step S200: analytic calculation of overlap area of fan-shaped projection beam and polar coordinate grid
The element value of the system matrix in the invention is equal to the overlapping area of the fan-shaped projection beam and the polar coordinate grid. FIG. 5 is a geometric relationship diagram of the matrix element analytic calculation of the system according to the present invention, in which the grid is shown
Figure BDA0001764271590000096
And the projection beam dmTaking the overlapped area as an example, the specific calculation steps are as follows:
1) grid mesh
Figure BDA0001764271590000097
Can be regarded as adjacent concentric circles rj、rj+1And two adjacent rays
Figure BDA0001764271590000098
The two arcs of the grid and the origin O form two overlapped sectors, and the projection beam d is calculated firstmAnd grid sector OC1C2The overlapping area of (a);
2) solving the side line L of the projection beam by using the geometric relationbAnd the concentric circle rjRay of radiation
Figure BDA0001764271590000099
Point of intersection A1、A2、P1、P2The coordinates of (a);
3) solving the projection beam d by using geometrical relationsmSide line L ofrAnd the concentric circle rjRay of radiation
Figure BDA00017642715900000910
Point of intersection B1、B2、Q1、Q2The coordinates of (a);
4) after the intersection is obtained, the area of the overlapped region shown in the figure is calculated as follows
Figure BDA00017642715900000911
Wherein
Figure BDA00017642715900000912
And
Figure BDA00017642715900000913
can be represented by the difference between the areas of a sector and a triangle, i.e.:
Figure BDA00017642715900000914
the areas of the fan and the triangle in equation (2) can be obtained by using the polar coordinates of the intersection points, and the area of the overlapping region can be calculated.
5) Repeating (1) - (4) to calculate the area of the overlap region in consideration of the overlap region of the projection fan and the other grid fan
Figure BDA0001764271590000101
The difference between the areas of the two overlapping regions is used to obtain the value of the system matrix element
amn＝Sm,i,j+1-Sm,i,j(3)
It should be noted that the overlapping manner of the projection fan and the polar grid shown in fig. 5 is only one of many overlapping manners. Fig. 6 shows several other overlapping manners of the projection beam and the polar grid, but in any overlapping manner, the system matrix element calculation can be performed by using the analysis method in step S200. The projected beams d are shown in FIGS. 6(a) - (d), respectivelymWith other overlap patterns in 4 of the polar grid.
Step S300: acquisition of system matrix
The acquisition of the system matrix is divided into two steps, firstly the system matrix on one projection angle is acquired, and then the system matrix of other projection directions is deduced according to the grid coordinate relation. The method comprises the following specific steps:
1) traversing the detector element and the polar coordinate grid, repeating the step S200, and obtaining a system matrix A ═ a on the current projection anglemnBecause the polar coordinate grid and the detector array have good symmetry, only half of the system matrix needs to be calculated in actual operation, and the other half can be directly obtained through a symmetric relation;
2) in the polar projection model shown in fig. 4, the light source S scans in a clockwise rotation, and the projection angle increases every time
Figure BDA0001764271590000102
Under the new projection angle, the area shape and the area of the area where the same projection beam is overlapped with the polar coordinate grid are the same, only the azimuth coordinate of the grid is changed from the azimuth coordinate of the grid
Figure BDA0001764271590000103
Become to
Figure BDA0001764271590000104
The coordinates of the corresponding system matrix elements in the matrix are thus from (m, N)r(j-1) + i) to (m, N)rJ-1) + i +1), the values are unchanged.
The polar coordinate system matrix calculation method in the two-dimensional reconstruction algorithm can obtain the system matrix with absolute calculation accuracy, and eliminates errors caused by the adoption of a numerical calculation method in the prior art. Numerical calculation errors can be reduced by increasing the number of trellis fill cells, but cannot be eliminated.
Comparative example 1 numerical calculation method calculation system matrix
The object to be processed was the same as in example 1, and the conventional numerical calculation method was used. The numerical calculation method used was carried out according to the methods disclosed in the prior art documents.
The calculation results of example 1 and comparative example 1 are shown in fig. 7, and fig. 7 shows the variation of the residual error of the numerical calculation method and the analytic calculation method according to the present invention with the number of numerical fill elements. As can be seen from fig. 7, the residual error of the conventional numerical calculation method in comparative example 1 (the open circle marked curve in fig. 7 is the error of the numerical method calculation) decreases as the density of the filling elements increases, but the time consumed for the calculation (the diamond marked curve in fig. 7 is the time of the numerical method calculation) increases greatly.
The analysis method provided by the invention can calculate an absolutely accurate system matrix in a short time (in the embodiment 1, the calculation time of the analysis method is the square mark point in fig. 7), and the calculation time obtained in the embodiment 1 is about 13 seconds.
Comparative example 1 provides a method to reduce the numerical residual to 10 according to the variation shown in fig. 7-2The density of filling elements must reach 802The calculation elapsed time was 60s, which is 4 times as long as the analysis method elapsed time. The analysis method provided by the invention can calculate the result without filling points, and has high calculation efficiency and high accuracy.
Referring to fig. 3, another aspect of the present invention further provides a polar coordinate system matrix calculation apparatus in a two-dimensional reconstruction algorithm, including:
the polar coordinate grid construction module 100 is configured to establish a polar coordinate projection model of the two-dimensional field to be reconstructed in an area surrounded by the point light source and the detection cell array, and divide the polar coordinate projection model at equal intervals to construct a polar coordinate grid;
system matrix element value calculation unit 200: used for calculating the overlapping area of the projection light beam generated by the point light source and the detection element and the polar coordinate grid under the ith projection angle as the system matrix element value amnIn the subscript, m is the serial number of the detector element, and n is the serial number of the grid;
the first system matrix calculation unit 300: the projection light beam traverses the polar coordinate grid and the detection element array under the ith projection angle to obtain a plurality of system matrix element values amnA plurality of values of elements a of said system matrixmnThe ith system matrix A is obtained by aggregationmn}；
The second system matrix calculation unit 400: the system matrix A is used for taking the projection angle of the i + g th system matrix A as { a ] according to the rotation invariance of the polar coordinate grid divisionmnEach of the system matrices in (1)Value of element amnAdjusting the position of the system to obtain the i + g system matrix;
the loop module 500: and returning the obtained g-g +1 to the second system matrix calculation unit 400, and repeating the process until the g-N and N are natural positive integers to obtain N i + g system matrices, and collecting the N i + g system matrices to obtain the system matrix of the two-dimensional field to be reconstructed.
Preferably, the polar grid construction module 100 includes:
a projection region construction module: the system is used for constructing an isosceles triangle projection area by taking the point light source as a vertex and the detection primitive array as a bottom edge;
a coordinate module: the system is used for establishing a polar coordinate system by taking the circle center of the two-dimensional field to be reconstructed as an origin, and determining the geometric positions of the point light source and the detection array element in the polar coordinate system;
a segmentation module: the method is used for dividing the two-dimensional field of view to be reconstructed into N at equal intervals by taking the circle center of the two-dimensional field of view to be reconstructed as a starting point and taking delta r as an intervalrA concentric circle with said two-dimensional field of view azimuth to be reconstructed
Figure BDA0001764271590000121
As a starting point, in
Figure BDA0001764271590000122
Equally-spaced partitioning the two-dimensional field of view to be reconstructed for an interval
Figure BDA0001764271590000123
A bar ray intersecting the concentric circles to construct the polar grid.
Preferably, the system matrix element value calculation unit 200 includes:
an intersection point calculation module: polar grid
Figure BDA0001764271590000124
Is a concentric circle rjConcentric circle rj+1Ray of radiation
Figure BDA0001764271590000125
And ray
Figure BDA0001764271590000126
Enclosing, calculating the projection beam and the concentric circle rjThe concentric circle rj+1The ray
Figure BDA0001764271590000127
And said ray
Figure BDA0001764271590000128
The coordinates of the intersection points;
an overlap area calculation module: for calculating the projection beam and the ray according to the intersection point coordinates
Figure BDA0001764271590000129
The ray
Figure BDA00017642715900001210
The concentric circles rjOverlap area S of enclosed sector aream,i,jCalculating the projection beam and the ray according to the intersection point coordinates
Figure BDA00017642715900001211
The ray
Figure BDA00017642715900001212
The concentric circles rj+1Overlap area S of enclosed sector aream,i,j+1；
An element value calculation module: value of matrix element a for the systemmnThe overlapping area Sm,i,j+1-said overlapping area Sm,i,jThe grid number N in the subscript is equal to Nr*(j-1)+i，NrThe number of the concentric circles.
Although the present invention has been described with reference to a few embodiments, it should be understood that the present invention is not limited to the above embodiments, but rather, the present invention is not limited to the above embodiments.

Claims (10)

1. A polar coordinate system matrix calculation method in a two-dimensional reconstruction algorithm is characterized by comprising the following steps:
step S100: establishing a polar coordinate projection model of a two-dimensional field to be reconstructed in an area surrounded by a point light source and a detection element array, and dividing the polar coordinate projection model at equal intervals to construct a polar coordinate grid;
step S200: calculating the overlapping area of the projection light beam generated by the point light source and the detection element and the polar coordinate grid under the ith projection angle as the system matrix element value amnIn the subscript, m is the serial number of the detector element, and n is the serial number of the grid;
step S300: under the ith projection angle, the projection light beam traverses the polar coordinate grid and the detection cell array to obtain a plurality of system matrix element values amnA plurality of values of elements a of said system matrixmnThe ith system matrix A is obtained by aggregationmn}；
Step S400: taking the projection angle of the (i + g) th system matrix A as { a ] according to the rotation invariance of the polar coordinate grid divisionmnEach of said system matrix element values a inmnAdjusting the position of the system to obtain an i + g system matrix;
step S500: and repeating the step S400 when g is g +1, and circulating until g is N, N is a natural positive integer to obtain N (i + g) th system matrixes, and collecting the N (i + g) th system matrixes to obtain the system matrix of the two-dimensional field to be reconstructed.
2. The method for calculating the matrix of the polar coordinate system in the two-dimensional reconstruction algorithm according to claim 1, wherein the step S100 comprises the steps of:
step S110: constructing an isosceles triangle projection area by taking the point light source as a vertex and the detection cell array as a bottom edge;
step S120: establishing a polar coordinate system by taking the circle center of the two-dimensional field to be reconstructed as an origin, and determining the geometric positions of the point light source and the detection primitive array in the polar coordinate system;
step S130: dividing the two-dimensional field to be reconstructed into N at equal intervals by taking the circle center of the two-dimensional field to be reconstructed as a starting point and taking delta r as an intervalrA concentric circle with said two-dimensional field of view azimuth to be reconstructed
Figure FDA0002325587630000011
As a starting point, in
Figure FDA0002325587630000012
Equally-spaced partitioning the two-dimensional field of view to be reconstructed for an interval
Figure FDA0002325587630000013
A bar ray intersecting the concentric circles to construct the polar grid.
3. The method for calculating the matrix of the polar coordinate system in the two-dimensional reconstruction algorithm according to claim 1, wherein the step S200 comprises the steps of:
step S210: polar grid
Figure FDA0002325587630000021
Is a concentric circle rjConcentric circle rj+1Ray of radiation
Figure FDA0002325587630000022
And ray
Figure FDA0002325587630000023
Enclosing, calculating the projection beam and the concentric circle rjThe concentric circle rj+1The ray
Figure FDA0002325587630000024
And said ray
Figure FDA0002325587630000025
The coordinates of the intersection points;
step S220: calculating the projection beam and the ray according to the intersection point coordinates
Figure FDA0002325587630000026
The ray
Figure FDA0002325587630000027
The concentric circles rjOverlap area S of enclosed sector aream,i,jCalculating the projection beam and the ray according to the intersection point coordinates
Figure FDA0002325587630000028
The ray
Figure FDA0002325587630000029
The concentric circles rj+1Overlap area S of enclosed sector aream,i,j+1；
Step S230: the system matrix element value amnThe overlapping area Sm,i,j+1-said overlapping area Sm,i,jThe grid number N in the subscript is equal to Nr*(j-1)+i，NrThe number of the concentric circles.
4. The method for calculating a matrix of a polar coordinate system in a two-dimensional reconstruction algorithm according to claim 1, wherein the two-dimensional field of view to be reconstructed is an axially symmetric region, and the projection beam at the ith projection angle in step S200 scans any one of the axially symmetric regions.
5. The method for calculating the matrix of the polar coordinate system in the two-dimensional reconstruction algorithm according to claim 2, wherein the angle difference between the ith projection angle and the (i + g) th projection angle is
Figure FDA00023255876300000210
6. The method for calculating the matrix of the polar coordinate system in the two-dimensional reconstruction algorithm according to claim 2, wherein the two-dimensional field of view to be reconstructed is a circle and is arranged in the projection area of the isosceles triangle, and the side of the two-dimensional field of view to be reconstructed is tangent to the waist line of the projection area of the isosceles triangle.
7. The method for calculating the matrix of the polar coordinate system in the two-dimensional reconstruction algorithm of claim 1, wherein the detection cell array is a linear array or a fan-shaped array.
8. A polar coordinate system matrix calculation device in a two-dimensional reconstruction algorithm is characterized by comprising:
the polar coordinate grid construction module is used for establishing a polar coordinate projection model of a two-dimensional field to be reconstructed in an area surrounded by the point light source and the detection element array, and dividing the polar coordinate projection model at equal intervals to construct a polar coordinate grid;
a system matrix element value calculation unit: used for calculating the overlapping area of the projection light beam generated by the point light source and the detection element and the polar coordinate grid under the ith projection angle as the system matrix element value amnIn the subscript, m is the serial number of the detector element, and n is the serial number of the grid;
a first system matrix calculation unit: the projection light beam traverses the polar coordinate grid and the detection element array under the ith projection angle to obtain a plurality of system matrix element values amnA plurality of values of elements a of said system matrixmnThe ith system matrix A is obtained by aggregationmn}；
A second system matrix calculation unit: the system matrix A is used for taking the projection angle of the i + g th system matrix A as { a ] according to the rotation invariance of the polar coordinate grid divisionmnEach of said system matrix element values a inmnAdjusting the position of the system to obtain an i + g system matrix;
a circulation module: and the system matrix is used for returning to the second system matrix calculation unit after g is g +1, and the process is circulated until g is N, N is a natural positive integer to obtain N (i + g) th system matrices, and the N (i + g) th system matrices are collected to obtain the system matrix of the two-dimensional field to be reconstructed.
9. The apparatus for calculating a polar coordinate system matrix in a two-dimensional reconstruction algorithm according to claim 8, wherein the polar grid construction module comprises:
a projection region construction module: the system is used for constructing an isosceles triangle projection area by taking the point light source as a vertex and the detection primitive array as a bottom edge;
a coordinate module: the system is used for establishing a polar coordinate system by taking the circle center of the two-dimensional field to be reconstructed as an origin, and determining the geometric positions of the point light source and the detection cell array in the polar coordinate system;
a segmentation module: the method is used for dividing the two-dimensional field of view to be reconstructed into N at equal intervals by taking the circle center of the two-dimensional field of view to be reconstructed as a starting point and taking delta r as an intervalrA concentric circle with said two-dimensional field of view azimuth to be reconstructed
Figure FDA0002325587630000031
As a starting point, in
Figure FDA0002325587630000032
Equally-spaced partitioning the two-dimensional field of view to be reconstructed for an interval
Figure FDA0002325587630000033
A bar ray intersecting the concentric circles to construct the polar grid.
10. The apparatus for calculating a matrix of a polar coordinate system in a two-dimensional reconstruction algorithm according to claim 8, wherein the system matrix element value calculating unit includes:
an intersection point calculation module: polar coordinatesGrid mesh
Figure FDA0002325587630000034
Is a concentric circle rjConcentric circle rj+1Ray of radiation
Figure FDA0002325587630000035
And ray
Figure FDA0002325587630000036
Enclosing, calculating the projection beam and the concentric circle rjThe concentric circle rj+1The ray
Figure FDA0002325587630000037
And said ray
Figure FDA0002325587630000038
The coordinates of the intersection points;
an overlap area calculation module: for calculating the projection beam and the ray according to the intersection point coordinates
Figure FDA0002325587630000039
The ray
Figure FDA00023255876300000310
The concentric circles rjOverlap area S of enclosed sector aream,i,jCalculating the projection beam and the ray according to the intersection point coordinates
Figure FDA0002325587630000041
The ray
Figure FDA0002325587630000042
The concentric circles rj+1Overlap area S of enclosed sector aream,i,j+1；
An element value calculation module: value of matrix element a for the systemmnThe overlapping area Sm,i,j+1-said overlapping area Sm,i,jThe grid number N in the subscript is equal to Nr*(j-1)+i，NrThe number of the concentric circles.
CN201810921211.7A 2018-06-27 2018-08-14 Polar coordinate system matrix calculation method and device in two-dimensional reconstruction algorithm Active CN109102553B (en)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201810677221 2018-06-27
CN2018106772210 2018-06-27

Publications (2)

Publication Number Publication Date
CN109102553A CN109102553A (en) 2018-12-28
CN109102553B true CN109102553B (en) 2020-05-05

Family

ID=64849486

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810921211.7A Active CN109102553B (en) 2018-06-27 2018-08-14 Polar coordinate system matrix calculation method and device in two-dimensional reconstruction algorithm

Country Status (1)

Country Link
CN (1) CN109102553B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110490925A (en) * 2019-07-19 2019-11-22 东北大学 A kind of method and apparatus that reconstruction region divides
CN111737381B (en) * 2020-05-11 2023-08-25 江苏北斗卫星应用产业研究院有限公司 Regional land parcel overlapping identification and overlapping area calculation method based on space-time big data
CN115512083B (en) * 2022-09-20 2023-04-11 广西壮族自治区地图院 Multi-inclination-angle numerical control sand table self-adaptive projection method

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1920950A (en) * 2006-09-25 2007-02-28 北京理工大学 Characteristic waveform decomposition and reconfiguration method based on Haar wavelet exaltation
CN107553476A (en) * 2017-10-12 2018-01-09 无锡信捷电气股份有限公司 Bipolar coordinate control system and implementation method based on CAD auxiliary

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN100447816C (en) * 2005-12-31 2008-12-31 清华大学 3D analysis and analog method for CT projection data
CN102289784A (en) * 2011-07-29 2011-12-21 中山大学 Method for simplifying digital geometric image based on point cloud model
CN104013403B (en) * 2014-05-23 2016-02-03 浙江理工大学 A kind of three-dimensional cardiac MR imaging method based on resolution of tensor sparse constraint
CN106462987B (en) * 2014-06-13 2019-11-19 美国西门子医疗解决公司 The improvement image reconstruction for volume based on data for projection collection
CN104020508A (en) * 2014-06-16 2014-09-03 中国矿业大学(北京) Direct ray tracking algorithm for geological radar chromatography detecting
CN106097291B (en) * 2016-06-08 2018-10-30 天津大学 A kind of ART flames slice reconstructing method based on gradient total variation

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1920950A (en) * 2006-09-25 2007-02-28 北京理工大学 Characteristic waveform decomposition and reconfiguration method based on Haar wavelet exaltation
CN107553476A (en) * 2017-10-12 2018-01-09 无锡信捷电气股份有限公司 Bipolar coordinate control system and implementation method based on CAD auxiliary

Also Published As

Publication number Publication date
CN109102553A (en) 2018-12-28

Similar Documents

Publication Publication Date Title
CN109102553B (en) Polar coordinate system matrix calculation method and device in two-dimensional reconstruction algorithm
US8204173B2 (en) System and method for image reconstruction by using multi-sheet surface rebinning
US7672424B2 (en) Image reconstruction with voxel dependent interpolation
CN100561520C (en) Rebuild the method and apparatus of (n+1) dimension image function from the Radon data
US7684539B2 (en) X-ray tomograph
CN102044081B (en) Reconstruction of 3d image datasets from X-ray cone-beam data
FI116324B (en) Method and arrangement for producing a tomography image
US20030202637A1 (en) True 3D cone-beam imaging method and apparatus
JP2000107168A (en) Three dimensional computerized tomography and tomograph
CN103679768B (en) Method for reconstructing CT image data with weighted backprojection, including computation unit and CT system for the method
US5708691A (en) X-ray computed tomographic imaging device and x-ray computed tomographic method
US5341460A (en) Method and apparatus for producing a three-dimensional computerized tomography image of an object with improved conversion of cone beam data to radon data
CN104834015A (en) Pulse flame four-dimensional measurement method based on deflectometric chromatography
US11501474B2 (en) Collimators for medical imaging systems and image reconstruction methods thereof
US6560308B1 (en) Method and system for approximating missing data in cone beam x-ray CT reconstruction
US20080013676A1 (en) Computerized Tomography Method With Helical Relative Movement and Conical Beam
Butler et al. Maximum a posteriori estimation for SPECT using regularization techniques on massively parallel computers
CN104050631A (en) Low-dose CT image reconstruction method
US20200029917A1 (en) Interior ct image generation method
CN101176117B (en) Reconstruction method for helical cone-beam CT
JPH10201751A (en) Method and device for eliminating boundary error in conical beam image
US7272205B2 (en) Methods, apparatus, and software to facilitate computing the elements of a forward projection matrix
Cho et al. Iterative image reconstruction using inverse Fourier rebinning for fully 3-D PET
Spencer et al. An automated peak fitting procedure for processing protein diffraction data from a linear position-sensitive detector
Gullberg et al. An iterative reconstruction algorithm for single photon emission computed tomography with cone beam geometry

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