CN100447816C - Three-dimensional Analytical Simulation Method of CT Projection Data - Google Patents

Abstract

The invention relates to a method for CT projection image reconstruction. A three-dimensional analytical simulation method for CT projection data disclosed by the present invention includes: 1) determining a model definition method; 2) determining a scanning mode definition method; 3) setting an analytical simulation program in a computer; 4) by 1), 2), and 3) Generate CT projection data p_jk. Since the present invention can be used to change the model definition and scan parameter definition in an interactive manner, various models and scan parameters can be defined as required without rewriting the simulation projection program, thereby obtaining analytical simulation projection data conveniently and quickly .

Description

Three-dimensional Analytical Simulation Method of CT Projection Data

technical field

The invention relates to a method for CT projection image reconstruction, in particular to a simulation method for CT projection image reconstruction.

Background technique

CT technology is widely used in the fields of industrial non-destructive testing (NDT), medical diagnosis and scientific research. CT projection data must be used in the research process of CT reconstruction algorithm.

There are two sources of CT projection data: 1) collected from an actual CT system; 2) generated by computer simulation.

Due to the high cost of actual CT systems, most research institutions do not have actual CT systems. In addition, in order to study the performance of the reconstruction algorithm, the projection data of a specific model is required. When collecting the projection data in the actual CT system, the corresponding model must be customized, which is very expensive, and due to the constraints of actual manufacturing capacity, some models cannot be manufactured. Based on the above reasons, computer simulation calculation projection data is widely used in the research process of CT reconstruction algorithm.

There are two types of computer simulation methods for computing projection data: 1) analytical simulation methods; 2) numerical simulation methods. The analytical simulation method obtains the simulated projection data through direct calculation of mathematical expressions, and the required model definition is directly written in the computer program code. If the model needs to be changed, the simulated projection calculation program must be rewritten, but the simulated projection data is ideal. Contains discretized noise. The numerical simulation method treats the model as a collection of pixels. The simulation process is to calculate each pixel and combine the results to obtain the simulated projection data. The model can be defined by each pixel. Changing the model does not require changing the simulation projection calculation program, so it is very convenient to use. However, the numerical simulation method has two shortcomings: 1) the calculation speed is slow; 2) discretization noise is introduced in the projection process. The way to reduce discretization noise is to reduce the pixel size, but this will lead to a cubic increase in computing time. For the study of CT algorithms, discretization noise is not expected in the simulated projection data, and the actual selection of the simulated projection calculation method needs to be based on the calculation time, the tolerance to discretized noise, and the convenience of model modification. In this process, the appropriate calculation method is selected. Under current conditions, the calculation speed of the analytical simulation method is usually tens to hundreds of times that of the numerical simulation method.

Contents of the invention

(1) Technical problems to be solved

The purpose of the present invention is to provide a three-dimensional analytical simulation method of CT projection data with fast calculation speed and easy model modification, so that the simulated projection data of various models without discretization noise can be obtained conveniently and quickly.

(2) Technical solutions

In order to achieve the above object, the present invention takes the following solutions: including: 1) determine the model definition method; 2) determine the scan mode definition method; 3) set the analytical simulation program in the computer; 4) generate by 1), 2), 3) CT projection data p ^jk .

Wherein, the model definition method includes: A defines a model by a combination of a plurality of basic geometric bodies with respective attenuation coefficients specified; B assumes that a model is composed of i basic geometric bodies _Mi , and the attenuation coefficient of each basic geometric body is _mi ; C carries on the model definition through the computer interactive interface.

Wherein, the scanning mode definition method includes:

a The complete simulated projection process consists of j projections;

探测器中心位置

b specifies for each projection the center position of the ray source for this projection

Detector center position

相对探测器中心位置

的几何关系；c specifies the number k of detection units of the detector and the center of each detector unit

relative detector center position

geometric relationship;

相对射线源中心位置的几何关系及各条射线与探测器交点相对探测器单元中心

的几何关系、各条射线对探测器单元的模拟投影数值的贡献权重w_jkl，权重应是归一化的，即

d specifies the number of rays received by each detector unit l, the starting point of each ray

Relative to the center position of the ray source The geometric relationship of each ray and the intersection point of the detector Relative to the center of the detector unit

The geometric relationship of , the contribution weight w _jkl of each ray to the analog projection value of the detector unit, the weight should be normalized, that is

e Define the scanning mode through the computer interactive interface.

Wherein, the analysis simulation program includes:

a) Divide the set physical factors to be simulated into two types: the influence on the ray position and the influence on the projection value. Adjust the start and end coordinates of the ray according to the effect on the ray position, and adjust the final simulated projection value according to the effect on the projection value.

b) Convert the coordinates of the starting point and end point of the ray to the local coordinate system of the geometry through coordinate transformation, and calculate the intersection coordinates of the ray and the geometry in the local coordinate system of the geometry.

c) For the direct superposition combination method, the length of the intersection line is calculated through the coordinates in the local coordinate system of the geometry, and then the attenuation coefficient is calculated.

d) For the alternative combination method, the coordinates of the intersection points are converted to the system coordinate system, and the distance from the two intersection points to the starting point of the ray is calculated in the system coordinate system, and then the attenuation coefficient is calculated.

Wherein, the CT projection data p _jk includes:

A) The simulated projection value p _jkl of a ray is a line segment ${\stackrel{r}{r}}_{\mathrm{jkl}}={\stackrel{r}{d}}_{\mathrm{jkl}}-{\stackrel{r}{\mathrm{the\; s}}}_{\mathrm{jkl}}$ The integral of all attenuation coefficients m on , that is:

${\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; jkl</span>}}<span\; class="notranslate">\; =</span>{\stackrel{<span\; class="notranslate">\; `</span>}{\mathrm{<span\; class="notranslate">\; o</span>}}}_{{\stackrel{\mathrm{<span\; class="notranslate">\; r</span>}}{\mathrm{<span\; class="notranslate">\; r</span>}}}_{\mathrm{<span\; class="notranslate">\; jkl</span>}}}\mathrm{<span\; class="notranslate">\; du</span>}$

B) The simulated projection value p _jk of each detector unit is:

C) The calculation method of p _jkl is related to the attenuation coefficient of the model. The attenuation coefficient of the model is obtained by combining the attenuation coefficients of each geometry in a certain way. Different combination methods must use different calculation methods to calculate p _jkl

位于组成此模型的m个几何体内，则模型该点的衰减系数为

这种组合方式下p_jkl的计算方法为：Wherein, the combination method includes: the direct accumulation of the attenuation coefficients of each geometric body, that is: a point in the model

Located within the m geometries that make up the model, the attenuation coefficient of this point in the model is

The calculation method of p _jkl in this combination is:

First calculate the line segment The intersection line length L _ijkl with each geometric body that makes up the model, and then calculate the simulated projection value p _jkl :

位于组成此模型的m个几何体内，则模型该点的衰减系数为 $m{r}_r=m_m,$ 这种组合方式下p_jkl的计算方法为：设线段

与组成模型的i个几何体中的N个几何体有两个交点，和第n个几何体的两个交点位于第n层，到

的距离分别是a_jkln和b_jkln，设a_jkln＜b_jkln。Wherein, the combination method includes: the replacement combination of the attenuation coefficients of each geometric body, that is: a point in the model

Located within the m geometries that make up the model, the attenuation coefficient of this point in the model is $m{r}_r=m_m,$ The calculation method of p _jkl in this combination is as follows: set the line segment

There are two intersections with N geometries among the i geometries that make up the model, and the two intersections with the nth geometry are located on the nth layer, to

The distances are a _jkln and b _jkln respectively, let a _jkln <b _jkln .

a _jkl0 = 0, $b_{\mathrm{jkl}0}=\left|{\stackrel{r}{r}}_{\mathrm{jkl}}\right|,$ m ₀ =0. Layer 0 is a line segment The line segments of each layer may be divided into shorter line segments by higher-level intersection points, and the attenuation coefficient value of each point is the attenuation coefficient value of the geometry corresponding to the highest layer line segment to which it belongs. From this, the simulated projection value p _jkl can be calculated.

(3) Beneficial effects

1) Since it can be used to change the model definition and scan parameter definition in an interactive way, various models and scan parameters can be defined as required without rewriting the simulation projection program, so as to obtain analytical simulation projection data conveniently and quickly . 2) The calculation speed is fast, and the simulated projection data of various models without discretization noise can be obtained.

Description of drawings

Fig. 1 is the simulation projection value calculation flow chart of a ray when the replacement combination of each geometry attenuation coefficient of the present invention;

Fig. 2 is a schematic diagram of the analytical simulation projection process of the present invention;

Fig. 3 is a schematic diagram of the model definition interactive interface of the present invention;

Fig. 4 is a schematic diagram of an interactive interface for setting simulated projection parameters in the present invention.

Fig. 5 is a calculation flow chart of the analytical simulation projection program of the present invention.

1) specific implementation

The following examples are used to illustrate the present invention, but are not intended to limit the scope of the present invention.

When the present invention is implemented, it includes steps: 1) determine the model definition method; 2) determine the scan mode definition method; 3) set the analytical simulation program in the computer; 4) generate CT by step 1), step 2), and step 3) Projection data p _jk .

As shown in FIG. 2 , an exemplary embodiment of the present invention includes a computer interaction interface and a method for processing user input by an analysis simulation program. Two coordinate systems are used in the embodiment: the system coordinate system and the geometry local coordinate system. The center of the model is located at the origin of the system coordinate system, and the positions of the ray source and the detector are defined in the system coordinate system. The center of the geometry is located at the origin of the local coordinate system of the geometry, and the symmetry axis coincides with the coordinate axis.

The interactive interface of model definition is shown in Figure 3. The user determines the coordinates (x _i , y _i , z _i ) of the center of the geometry that compose the model in the system coordinate system and the inclination angle a _i of the symmetry axis of the geometry in the system coordinate system , f _i , and at the same time determine the geometric attenuation coefficient m _i and its combination.

The interactive interface for setting simulation projection parameters is shown in Figure 4. The user determines the scanning mode and the physical factors to be simulated.

The process of parsing the simulation program is shown in Figure 5, and the processing method for user input is as follows:

1) The physical influence factors to be simulated are decomposed into two categories: the influence on the ray position (such as various offsets, etc.) and the influence on the projection value on the ray (such as noise, etc.).

和终点坐标 2) Calculate the starting point coordinates of each ray according to the setting of scanning parameters and the setting of physical factors affecting the ray position

and end point coordinates

3) Transform the starting point coordinates and the ending point coordinates to obtain the starting point coordinates of each ray in the local coordinate system of each geometry and end point coordinates

和

4) Calculate the intersection coordinates of the ray and the geometry in the local coordinate system of each geometry: there are no two intersections (no intersection or only one intersection); or there are two intersections

and

If the attenuation coefficients of the geometry are combined in a direct accumulation manner, the length L _ijkl of the intersection line between the ray and each geometry is calculated according to the distance formula between two points, and then the projection value p _jkl of each ray can be calculated. The calculation method of p _jkl is:

First calculate the line segment The intersection line length L _ijkl with each geometric body that makes up the model, and then calculate the simulated projection value p _jkl :

和

再计算交点到射线起点的距离a_ijkl和b_ijkl，然后即可计算得到各条射线的投影值p_jkl。If the attenuation coefficients of the geometry are combined by replacement, the coordinates of the intersection point are transformed to obtain their coordinates in the system coordinate system

and

Then calculate the distances a _ijkl and b _ijkl from the intersection point to the starting point of the ray, and then calculate the projection value p _jkl of each ray.

The corresponding calculation formula is:

${\mathrm{<span\; class="notranslate">\; a</span>}}_{\mathrm{<span\; class="notranslate">\; ijkl</span>}}<span\; class="notranslate">\; =</span>\sqrt{{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{{\mathrm{<span\; class="notranslate">\; a</span>}}_{\mathrm{<span\; class="notranslate">\; ijkl</span>}}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; jkl</span>}}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{{\mathrm{<span\; class="notranslate">\; a</span>}}_{\mathrm{<span\; class="notranslate">\; ijkl</span>}}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; jkl</span>}}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{{\mathrm{<span\; class="notranslate">\; a</span>}}_{\mathrm{<span\; class="notranslate">\; ijkl</span>}}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; jkl</span>}}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; ,</span>$

${\mathrm{<span\; class="notranslate">\; b</span>}}_{\mathrm{<span\; class="notranslate">\; ijkl</span>}}<span\; class="notranslate">\; =</span>\sqrt{{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{{\mathrm{<span\; class="notranslate">\; b</span>}}_{\mathrm{<span\; class="notranslate">\; ijkl</span>}}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; jkl</span>}}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{{\mathrm{<span\; class="notranslate">\; b</span>}}_{\mathrm{<span\; class="notranslate">\; ijkl</span>}}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; jkl</span>}}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{{\mathrm{<span\; class="notranslate">\; b</span>}}_{\mathrm{<span\; class="notranslate">\; ijkl</span>}}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; jkl</span>}}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; .</span>$

The contribution weight w _jkl of each ray to the analog projection value of the detector unit can be selected as:

${\mathrm{<span\; class="notranslate">\; w</span>}}_{\mathrm{<span\; class="notranslate">\; jkl</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}$

As shown in Figure 1, under the combination of the above replacement methods, the calculation method of p _jkl is:

与组成模型的i个几何体中的N个几何体有两个交点，和第n个几何体的两个交点位于第n层，到

的距离分别是a_jkln和b_jkln，设a_jkln＜b_jkln。a_jkl0＝0， $b_{\mathrm{jkl}0}=\left|{\stackrel{r}{r}}_{\mathrm{jkl}}\right|,$ m₀＝0。第0层为线段

各层的线段都可能被更高层的交点分割为更短的线段，各点的衰减系数值为其所属的最高层线段对应的几何体的衰减系数值。由此按图1所示的方法计算模拟投影数值p_jkl。将所有的交点按到

的距离从小到大的顺序排列，然后以此查找各位于最高层的线段。使用一个变量p_jkl保存模拟投影值，初值置为0。使用变量S_c、E_c保存当前正在处理的线段起点、终点位置，变量n_c保存当前处理的线段所在层。line segment

There are two intersections with N geometries among the i geometries that make up the model, and the two intersections with the nth geometry are located on the nth layer, to

The distances are a _jkln and b _jkln respectively, let a _jkln <b _jkln . a _jkl0 = 0, $b_{\mathrm{jkl}0}=\left|{\stackrel{r}{r}}_{\mathrm{jkl}}\right|,$ m ₀ =0. Layer 0 is a line segment

The line segments of each layer may be divided into shorter line segments by higher-level intersection points, and the attenuation coefficient value of each point is the attenuation coefficient value of the geometry corresponding to the highest layer line segment to which it belongs. From this, the simulated projection value p _jkl is calculated according to the method shown in FIG. 1 . Press all intersection points to

The distances are arranged in descending order, and then use this to find the line segments at the highest level. Use a variable p _jkl to save the simulated projection value, and set the initial value to 0. Use variables S _c and E _c to save the starting point and end position of the line segment currently being processed, and variable n _c to save the layer where the line segment currently being processed is located.

较近的点)，则先计算E_c，E_c为所有高于当前层n_c且大于S_c的a_jkln类点和当前层对应的b_jkln类点(即第n层中距离

较远的点)b_jlnc中最小的点，然后重新计算当前层n_c，将n_c置为E_c所在层；如果当前处理的线段起点S_c为b_jkln类点，则先重新计算当前层n_c，将n_c置为满足对应的a_jkln类点小于S_c、对应的b_jkln类点小于E_c(此时为前一次处理的线段终点)且本身层数小于n_c的所有层中的最大值，然后计算当前处理的线段终点E_c，E_c为所有高于当前层n_c且大于S_c的a_jkln类点和当前层对应的b_jkln类点(即第n层中距离

较远的点)b_jklnc中最小的点。随后查找当前层n_c对应的衰减系数m_nc，然后将(E_c-S_c)′m_nc累加到p_jkl中，这样就完成了一次线段处理。随后将当前处理线段起点置为前一次处理的线段终点，然后重复以上的线段处理过程，直到起点到达所有点中的最大值b_jkl0时即完成计算。S _c starts from 0, find the corresponding E _c and n _c according to the following method: If the starting point S _c of the line segment currently processed is a _jkln type point (that is, the distance in the nth layer

closer points), then calculate E _c first, E _c is all a _jkln class points higher than the current layer n _c and greater than S _c and b _jkln class points corresponding to the current layer (that is, the middle distance of the nth layer

The farther point) the smallest point in b _jlnc , then recalculate the current layer n _c , and set n _c as the layer where E _c is located; if the starting point S _c of the line segment currently processed is a b _jkln class point, then recalculate the current layer first n _c , set n _c to satisfy that the corresponding a _jkln class point is smaller than S _c , the corresponding b _jkln class point is smaller than E _c (at this time, it is the end point of the line segment of the previous processing) and the number of layers is less than n _c in all layers , and then calculate the currently processed end point E _c of the line segment, E _c is all a _jkln class points higher than the current layer n _c and greater than S _c and b _jkln class points corresponding to the current layer (that is, the nth layer middle distance

the farther point) b the smallest point in _jklnc . Then find the attenuation coefficient m _nc corresponding to the current layer n _c , and then add (E _c -S _c )′m _nc to p _jkl , thus completing a line segment processing. Then set the starting point of the current processing line segment as the end point of the previous processing line segment, and then repeat the above line segment processing process until the starting point reaches the maximum value b _jkl0 of all points to complete the calculation.

Claims (5)

1. A three-dimensional analytical simulation method for CT projection data, characterized in that it comprises: 1) Determine the model definition method; 2) Determine the scan mode definition method, the scan mode definition method includes: a The complete simulated projection process consists of j projections; b specifies for each projection the center position of the ray source for this projection Detector center position c specifies the number k of detection units of the detector and the center of each detector unit

relative detector center position

geometric relationship; d specifies the number of rays received by each detector unit l, the starting point of each ray

Relative to the center position of the ray source

The geometric relationship of each ray and the intersection point of the detector

Relative to the center of the detector unit The geometric relationship of , the contribution weight w _jkl of each ray to the analog projection value of the detector unit, the weight should be normalized, that is

e definition of the scanning method through the computer interactive interface; 3) set analysis simulation program in computer; 4) Generate CT projection data p _jk from 1), 2), and 3). 2. The method for three-dimensional analysis and simulation of CT projection data as claimed in claim 1, characterized in that said model definition method comprises: A defining a model through the combination of a plurality of basic geometric bodies with respective attenuation coefficients specified; B a model consisting of i The attenuation _coefficient of each basic geometry is _mi ; C defines the model through the computer interface. 3. The method for three-dimensional analysis and simulation of CT projection data according to claim 1, wherein the analysis and simulation program includes: a) Divide the set physical factors to be simulated into two categories: influence on the position of the ray and influence on the projection value, adjust the coordinates of the starting point and end point of the ray according to the influence on the position of the ray, and adjust the final simulated projection value according to the influence on the projection value ; b) Convert the coordinates of the starting point and end point of the ray to the local coordinate system of the geometry through coordinate transformation, and calculate the intersection coordinates of the ray and the geometry in the local coordinate system of the geometry; c) For the direct superposition combination method, the length of the intersection line is calculated by the coordinates in the local coordinate system of the geometry, and then the attenuation coefficient is calculated; d) For the alternative combination method, the coordinates of the intersection point are converted to the system coordinate system, and the distance from the two intersection points to the starting point of the ray is calculated in the system coordinate system, and then the attenuation coefficient is calculated. 4. The method for three-dimensional analysis and simulation of CT projection data according to claim 2, characterized in that the combination method includes: direct accumulation of the attenuation coefficients of each geometric body, that is, a point in the model

Located within the m geometries that make up the model, the attenuation coefficient of this point in the model is

The calculation method of p _jkl in this combination is: First calculate the line segment The intersection line length L _ijkl with each geometric body that makes up the model, and then calculate the simulated projection value p _jkl :

5. The method for three-dimensional analysis and simulation of CT projection data according to claim 2, characterized in that the combination method includes: the replacement combination of the attenuation coefficients of each geometric body, that is, a point in the model

Located within the m geometries that make up the model, the attenuation coefficient of this point in the model is The calculation method of p _jkl in this combination is: line segment

There are two intersections with N geometries among the i geometries that make up the model, and the two intersections with the nth geometry are located on the nth layer, to

The distances are a _jkln and b _jkln respectively, a _jkln <b _jkln , a _jkl0 =0,

m ₀ ＝0, the 0th layer is a line segment

The line segments of each layer are divided into shorter line segments by the intersection points of higher layers, and the attenuation coefficient value of each point is the attenuation coefficient value of the geometry corresponding to the highest layer line segment to which it belongs, from which the simulated projection value p _jkl is calculated.

Images