# Application of Image Reconstruction Based on Inverse Radon Transform in CT System Parameter Calibration and Imaging

Authors: Chun-li Wang; Hao-chi Zhang; Ze-yu Zeng; Jun-hui Yu; Yang Wang
Journal: Complexity
 (2021)

---

## Abstract

According to the requirements of CT system parameter calibration and imaging, using OpenCV and MATLAB software, the reverse Radon transform and the Canny edge detection algorithm in the projection edge methods can determine the position and geometry of two different media in the square tray based on given data. Besides, it can apply and analyze the shape and the absorption rate of the specified point, the accuracy, and stability of the template calibration parameters by this way, which enables designing a new template and calibrating the new parameters.


---

## Body

## 1. Introduction

CT is called computed tomography. Under the premise of not destroying the sample, using the absorption characteristics of the sample to absorb energy can be used for tomographic imaging of biological tissues, and the internal structure information of the sample can be obtained [1]. The existence of certain CT systems when installing errors affects the accuracy of imaging, which in turn affects the analysis of samples. It is necessary to calibrate the parameters of the installed CT system. Using the standardized sample template can calibrate the CT system parameters, and it can obtain the accurate imaging effect of samples with unknown structures [2].In recent years, foreign scholars have put forward innovative methods for image reconstruction. In 2021, Yuan [3] used image reconstruction technology based on micro-computed tomography (micro-CT). They obtain the distribution of pressure and water flow inside the foamed concrete in order to study the permeability characteristics of foamed concrete, so as to conduct a numerical simulation. Szczykutowicz et al. [4] proposed a proprietary deep learning image reconstruction method (DLIR) based on deep learning. Liu et al. [5, 6] considered the reconstruction of the point set of the graph, which also provided some inspiration for the innovation of image reconstruction methods. In edge monitoring, Asghari and Jalali [7] proposed to transform images by simulating the propagation of light through physical media with specific distorted diffraction characteristics, effectively revealing the transition of image intensity. Francesco et al. [8] proposed a one-dimensional antinoise edge detection algorithm to improve the accuracy of fire front prediction.In this study, according to the given data, the inverse Radon transform and Canny edge detection algorithm were used in the projection edge method to determine the positions and geometric shapes of the two different media in the square tray, so as to design a new template and calibrate new parameters.


## 2. Data Sources and Model Assumptions

The data in this paper were derived from the 2017 National College Students Mathematical Modeling Contest A. In order to facilitate the processing of the problem, the following assumptions are made: (1) the radiation source is received by the sample medium coming from the CT system light source only; (2) the CT system works are in an ideal environment, and various uncontrollable events are ignored; (3) the light has the same intensity; (4) the calibration template is spherical and elliptical on the two-dimensional plane; and (5) the new model setting differs from the initial template only in terms of geometry.


## 3. Determination of the Indicators when Uniform Solid Medium Is the Calibration Template

### 3.1. Research Ideas

According to the meaning of the question, two uniform fixed media are placed on a square tray, and the absorption strength is the same. Using MATLAB can draw the received information heat map of the calibration template composed of uniform solid media. If the rotation center of the CT system happens to be at the geometric center of the tray, it will appear that the projection map should be symmetrical. If it is not symmetrical, it means that the CT system rotation center is not at the geometric center of the tray. According to the relevant properties on the heat map and the image processing method, it can solve the position of the rotation center on the square tray and the distance between the detector units and 180 directions using by the CT system.


### 3.2. Data Processing

There is an error in the installation of the CT system, and its rotation center will have a certain offset, which will affect the imaging quality. Therefore, it is necessary to calibrate the parameters of the CT system with a sample with a known structure. Figure1 shows a schematic diagram of the rotation center offset when the object is not moving, as shown in Figure 1.Figure 1
Schematic diagram of the rotation center offset when the object is not moving.Use MATLAB to analyze the geometric information data file of the template. The center of the ellipse is the centerO of the square tray as the origin of the coordinate system, the direction of the long axis of the ellipse is the axis y, and the direction of the short axis of the ellipse is the axis x. Using it can establish a rectangular coordinate system, which is called the original coordinate system.


### 3.3. Research Methods and Result Analysis

#### 3.3.1. The Solution of Distance between Detector Units

Let the received information of the calibration template composed of a uniform solid medium beΓ. Γ reflects that the CT system light source illuminates the regular template, causing different detectors to receive the gain processing medium. This medium absorbs the attenuated ray energy. It rotates along a certain fixed center and the center is totated 180 times. Figure 2 shows the relationship between the projection data of the light source on the template and the absorption intensity Uij of each point (each slice) on the template. Assuming it is a 3×3 unit composition, give the constraints of Γ followed, that is, the ray scan schematic [1].Figure 2
Schematic diagram of ray scan.The constraint equation is(1)U11+U12+U13=E11,U11=E21,U31+U21+U11=E31,U21+U22+U23=E12,U31+U22+U13=E22,U32+U22+U12=E32,U31+U32+U33=E13,U32+U23=E23,U33+U23+U13=E33.The above constraint equation is the constraint equation of the detector after ray scanning. FromΓ, It can be seen that the CT system and its detector satisfy the above constraint equation, and its detector and light source both have 512 0≤j≤512, and its rotation angle is 0≤i≤180. Using MATLAB can make received information of the template into a heat map, and the darker the color, the larger the representative value of Eij. The image reflects the projected image of the template on the detector, which is the situation when the detector is rotated 180° and 360°, respectively, as shown in Figures 3–5 .Figure 3
Three-dimensional projection heat map from different angles.Figure 4
Projection heat map of180° (flat).Figure 5
Projection heat map of360° (flat).Analyzing the projection diagram shows that the projection diagram is divided into the projection diagram of a small ball and the projection diagram of an ellipse. The widest position of the ellipse projection is the position where the long axis of the ellipse is perpendicular to the light source and parallel to the receiver. The narrowest position is the short axis of the ellipse and the position where the light source is perpendicular to the receiver. Because the ball is regular, it can be seen from the 180° projection heat map and the 360° projection heat map that the motion projection trajectory of the ball is a sinusoidal (cosine) curve. The widthD=8mm of the projection curve of the ball at any moment is the same.According to the conclusion drawn from the heat map, the small ball is the key to solving this problem because its projection image has obvious periodicity. Because the small ball has a regular shape, the distanceL between the detectors can be solved. Let the number of detectors be N, from Figure 4 a small ball template can be obtained, and a detector of N=29 is required to complete the scan. The diameter of the ball divided by the number of detectors can be obtained by the detector spacing 0.2759 mm.


#### 3.3.2. The Solution of 180 Directions of X-Ray

According to the results of 2.3 (1) and the conclusions drawn from the projection heat map, the data files and corresponding images can be drawn. Assuming that the direction shown in Figure2 is 0°, the light source rotates counterclockwise, and the light source illuminates that the minor axis direction vertically is starting direction. When the light source rotates to the widest position and the shortest position, the light source is perpendicular to the coordinate axis. It is easy to know the scale when the light source rotates 90°. Get the initial rotation angle of CT and the angle of each rotation, the initial rotation angle φ0=29.3478° , and the algorithm θr=0.9783° of each rotation angle. So, the total rotation angle is Θt=176.0940° approximately.


#### 3.3.3. The Coordinates of the Center of Rotation

Since the CT system needs to be rotated and the center of rotation will have a certain offset, the situation is more complicated. Assume that the center of rotation is moved fromo to o′, the translation distance a in the horizontal axis direction, and the translation distance b in the vertical axis direction. At this time, the detection coordinate system is so′t.The position rotated to the vertical irradiation short axis is used as one of the reference points. In the range of 0–360°, there are two such points. These two points are the brightest positions in the heat map, and the difference is exactly half a cycle. And there is a central rotational symmetry between the two about the center of rotation [9]. The deviation is reflected by the number of receivers. That is, the result can be obtained by L=D/N converting it into a real distance. By selecting the receiver corresponding to the peak of the narrowest position of the image, the position number is recorded as N1, and the completed picture gets the receiver position number N2=512−N+1 of the special position. Then, the offset of x can be calculated:(2)δx=DN1−N22=−0.2759×291−2222=−9.5185mm.Select the receiver position number corresponding to the peak of the widest position of the image and write it asN3, and the completed picture will get the receiver position number N4=512−N+1 of the special position. Then, the offset of x can be calculated:(3)δy=DN3−N42=0.2759×277−2362=5.656mm.The coordinates of the rotation center are−9.5185mm,5.6560mm, and the related results are shown in Table 1.Table 1
Uniform solid medium is the result of various indexes of calibration template.
TitleResultsTotal rotation angleΘ≈176.0940°The center of rotation is in the position of the square pallet−9.5185mm,5.6560mmDistance between detector unitsL=0.2759mmThe CT system uses 180 directions ofX-rays:Initial rotation angle φ0/per rotation angle algorithm θrφ0=29.3478°θr=0.9783°


## 3.1. Research Ideas

According to the meaning of the question, two uniform fixed media are placed on a square tray, and the absorption strength is the same. Using MATLAB can draw the received information heat map of the calibration template composed of uniform solid media. If the rotation center of the CT system happens to be at the geometric center of the tray, it will appear that the projection map should be symmetrical. If it is not symmetrical, it means that the CT system rotation center is not at the geometric center of the tray. According to the relevant properties on the heat map and the image processing method, it can solve the position of the rotation center on the square tray and the distance between the detector units and 180 directions using by the CT system.


## 3.2. Data Processing

There is an error in the installation of the CT system, and its rotation center will have a certain offset, which will affect the imaging quality. Therefore, it is necessary to calibrate the parameters of the CT system with a sample with a known structure. Figure1 shows a schematic diagram of the rotation center offset when the object is not moving, as shown in Figure 1.Figure 1
Schematic diagram of the rotation center offset when the object is not moving.Use MATLAB to analyze the geometric information data file of the template. The center of the ellipse is the centerO of the square tray as the origin of the coordinate system, the direction of the long axis of the ellipse is the axis y, and the direction of the short axis of the ellipse is the axis x. Using it can establish a rectangular coordinate system, which is called the original coordinate system.


## 3.3. Research Methods and Result Analysis

### 3.3.1. The Solution of Distance between Detector Units

Let the received information of the calibration template composed of a uniform solid medium beΓ. Γ reflects that the CT system light source illuminates the regular template, causing different detectors to receive the gain processing medium. This medium absorbs the attenuated ray energy. It rotates along a certain fixed center and the center is totated 180 times. Figure 2 shows the relationship between the projection data of the light source on the template and the absorption intensity Uij of each point (each slice) on the template. Assuming it is a 3×3 unit composition, give the constraints of Γ followed, that is, the ray scan schematic [1].Figure 2
Schematic diagram of ray scan.The constraint equation is(1)U11+U12+U13=E11,U11=E21,U31+U21+U11=E31,U21+U22+U23=E12,U31+U22+U13=E22,U32+U22+U12=E32,U31+U32+U33=E13,U32+U23=E23,U33+U23+U13=E33.The above constraint equation is the constraint equation of the detector after ray scanning. FromΓ, It can be seen that the CT system and its detector satisfy the above constraint equation, and its detector and light source both have 512 0≤j≤512, and its rotation angle is 0≤i≤180. Using MATLAB can make received information of the template into a heat map, and the darker the color, the larger the representative value of Eij. The image reflects the projected image of the template on the detector, which is the situation when the detector is rotated 180° and 360°, respectively, as shown in Figures 3–5 .Figure 3
Three-dimensional projection heat map from different angles.Figure 4
Projection heat map of180° (flat).Figure 5
Projection heat map of360° (flat).Analyzing the projection diagram shows that the projection diagram is divided into the projection diagram of a small ball and the projection diagram of an ellipse. The widest position of the ellipse projection is the position where the long axis of the ellipse is perpendicular to the light source and parallel to the receiver. The narrowest position is the short axis of the ellipse and the position where the light source is perpendicular to the receiver. Because the ball is regular, it can be seen from the 180° projection heat map and the 360° projection heat map that the motion projection trajectory of the ball is a sinusoidal (cosine) curve. The widthD=8mm of the projection curve of the ball at any moment is the same.According to the conclusion drawn from the heat map, the small ball is the key to solving this problem because its projection image has obvious periodicity. Because the small ball has a regular shape, the distanceL between the detectors can be solved. Let the number of detectors be N, from Figure 4 a small ball template can be obtained, and a detector of N=29 is required to complete the scan. The diameter of the ball divided by the number of detectors can be obtained by the detector spacing 0.2759 mm.


### 3.3.2. The Solution of 180 Directions of X-Ray

According to the results of 2.3 (1) and the conclusions drawn from the projection heat map, the data files and corresponding images can be drawn. Assuming that the direction shown in Figure2 is 0°, the light source rotates counterclockwise, and the light source illuminates that the minor axis direction vertically is starting direction. When the light source rotates to the widest position and the shortest position, the light source is perpendicular to the coordinate axis. It is easy to know the scale when the light source rotates 90°. Get the initial rotation angle of CT and the angle of each rotation, the initial rotation angle φ0=29.3478° , and the algorithm θr=0.9783° of each rotation angle. So, the total rotation angle is Θt=176.0940° approximately.


### 3.3.3. The Coordinates of the Center of Rotation

Since the CT system needs to be rotated and the center of rotation will have a certain offset, the situation is more complicated. Assume that the center of rotation is moved fromo to o′, the translation distance a in the horizontal axis direction, and the translation distance b in the vertical axis direction. At this time, the detection coordinate system is so′t.The position rotated to the vertical irradiation short axis is used as one of the reference points. In the range of 0–360°, there are two such points. These two points are the brightest positions in the heat map, and the difference is exactly half a cycle. And there is a central rotational symmetry between the two about the center of rotation [9]. The deviation is reflected by the number of receivers. That is, the result can be obtained by L=D/N converting it into a real distance. By selecting the receiver corresponding to the peak of the narrowest position of the image, the position number is recorded as N1, and the completed picture gets the receiver position number N2=512−N+1 of the special position. Then, the offset of x can be calculated:(2)δx=DN1−N22=−0.2759×291−2222=−9.5185mm.Select the receiver position number corresponding to the peak of the widest position of the image and write it asN3, and the completed picture will get the receiver position number N4=512−N+1 of the special position. Then, the offset of x can be calculated:(3)δy=DN3−N42=0.2759×277−2362=5.656mm.The coordinates of the rotation center are−9.5185mm,5.6560mm, and the related results are shown in Table 1.Table 1
Uniform solid medium is the result of various indexes of calibration template.
TitleResultsTotal rotation angleΘ≈176.0940°The center of rotation is in the position of the square pallet−9.5185mm,5.6560mmDistance between detector unitsL=0.2759mmThe CT system uses 180 directions ofX-rays:Initial rotation angle φ0/per rotation angle algorithm θrφ0=29.3478°θr=0.9783°


## 3.3.1. The Solution of Distance between Detector Units

Let the received information of the calibration template composed of a uniform solid medium beΓ. Γ reflects that the CT system light source illuminates the regular template, causing different detectors to receive the gain processing medium. This medium absorbs the attenuated ray energy. It rotates along a certain fixed center and the center is totated 180 times. Figure 2 shows the relationship between the projection data of the light source on the template and the absorption intensity Uij of each point (each slice) on the template. Assuming it is a 3×3 unit composition, give the constraints of Γ followed, that is, the ray scan schematic [1].Figure 2
Schematic diagram of ray scan.The constraint equation is(1)U11+U12+U13=E11,U11=E21,U31+U21+U11=E31,U21+U22+U23=E12,U31+U22+U13=E22,U32+U22+U12=E32,U31+U32+U33=E13,U32+U23=E23,U33+U23+U13=E33.The above constraint equation is the constraint equation of the detector after ray scanning. FromΓ, It can be seen that the CT system and its detector satisfy the above constraint equation, and its detector and light source both have 512 0≤j≤512, and its rotation angle is 0≤i≤180. Using MATLAB can make received information of the template into a heat map, and the darker the color, the larger the representative value of Eij. The image reflects the projected image of the template on the detector, which is the situation when the detector is rotated 180° and 360°, respectively, as shown in Figures 3–5 .Figure 3
Three-dimensional projection heat map from different angles.Figure 4
Projection heat map of180° (flat).Figure 5
Projection heat map of360° (flat).Analyzing the projection diagram shows that the projection diagram is divided into the projection diagram of a small ball and the projection diagram of an ellipse. The widest position of the ellipse projection is the position where the long axis of the ellipse is perpendicular to the light source and parallel to the receiver. The narrowest position is the short axis of the ellipse and the position where the light source is perpendicular to the receiver. Because the ball is regular, it can be seen from the 180° projection heat map and the 360° projection heat map that the motion projection trajectory of the ball is a sinusoidal (cosine) curve. The widthD=8mm of the projection curve of the ball at any moment is the same.According to the conclusion drawn from the heat map, the small ball is the key to solving this problem because its projection image has obvious periodicity. Because the small ball has a regular shape, the distanceL between the detectors can be solved. Let the number of detectors be N, from Figure 4 a small ball template can be obtained, and a detector of N=29 is required to complete the scan. The diameter of the ball divided by the number of detectors can be obtained by the detector spacing 0.2759 mm.


## 3.3.2. The Solution of 180 Directions of X-Ray

According to the results of 2.3 (1) and the conclusions drawn from the projection heat map, the data files and corresponding images can be drawn. Assuming that the direction shown in Figure2 is 0°, the light source rotates counterclockwise, and the light source illuminates that the minor axis direction vertically is starting direction. When the light source rotates to the widest position and the shortest position, the light source is perpendicular to the coordinate axis. It is easy to know the scale when the light source rotates 90°. Get the initial rotation angle of CT and the angle of each rotation, the initial rotation angle φ0=29.3478° , and the algorithm θr=0.9783° of each rotation angle. So, the total rotation angle is Θt=176.0940° approximately.


## 3.3.3. The Coordinates of the Center of Rotation

Since the CT system needs to be rotated and the center of rotation will have a certain offset, the situation is more complicated. Assume that the center of rotation is moved fromo to o′, the translation distance a in the horizontal axis direction, and the translation distance b in the vertical axis direction. At this time, the detection coordinate system is so′t.The position rotated to the vertical irradiation short axis is used as one of the reference points. In the range of 0–360°, there are two such points. These two points are the brightest positions in the heat map, and the difference is exactly half a cycle. And there is a central rotational symmetry between the two about the center of rotation [9]. The deviation is reflected by the number of receivers. That is, the result can be obtained by L=D/N converting it into a real distance. By selecting the receiver corresponding to the peak of the narrowest position of the image, the position number is recorded as N1, and the completed picture gets the receiver position number N2=512−N+1 of the special position. Then, the offset of x can be calculated:(2)δx=DN1−N22=−0.2759×291−2222=−9.5185mm.Select the receiver position number corresponding to the peak of the widest position of the image and write it asN3, and the completed picture will get the receiver position number N4=512−N+1 of the special position. Then, the offset of x can be calculated:(3)δy=DN3−N42=0.2759×277−2362=5.656mm.The coordinates of the rotation center are−9.5185mm,5.6560mm, and the related results are shown in Table 1.Table 1
Uniform solid medium is the result of various indexes of calibration template.
TitleResultsTotal rotation angleΘ≈176.0940°The center of rotation is in the position of the square pallet−9.5185mm,5.6560mmDistance between detector unitsL=0.2759mmThe CT system uses 180 directions ofX-rays:Initial rotation angle φ0/per rotation angle algorithm θrφ0=29.3478°θr=0.9783°


## 4. Use the Inverse Radon Transform to Reconstruct the Graphics of the CT Receiver

### 4.1. Research Ideas

According to the received information and the projection map, the problem obtains the position, geometry, and absorptivity of the unknown medium on the square tray and specifically gives the absorptivity of 10 positions. Since the known information is the received information on the detector. It is necessary to consider using the Radon transform of CT system image reconstruction to reconstruct the image [9].


### 4.2. Data Processing

The definition of the two-dimensional Radon transform is(4)Rfs,θ=∫−∞+∞∫−∞+∞fx,y∂xcosθ+ysinθ−sdxdy.Among which,fx,y is the function of the object to be reconstructed. Rfs,θ is expressed as the Radon transformation along a certain straight line when the CT system rotates in the angle of view θ. That is, the projection along the straight line, fx,y, denoted as ps,θ. s is the distance from the origin to a certain straight line. Theoretically, the image that can be reconstructed is regarded as the inverse deduction fx,y of the Radon projection transformation function ps,θ, that is, the inverse Radon transformation problem. The formula is(5)fx,y=12π2∫02π∫−∞+∞∂ps,θ/∂sxcosθ+ysinθ−sdsdθ.Using the reception information of an unknown medium obtained by the above-mentioned CT system, we need to use the inverse Radon transform to reconstruct its image.The above geometric shape is the Shepp–Logan head model [10]. In Figure 6(a), the absorption rate of the ten positions given in the title is calibrated. The absorption rate at ten locations in the location diagram are shown in Table 2.Figure 6
The geometric shape of the received information of an unknown medium. (a) Geometry and position on the square pallet. (b) Map of receiving information of unknown medium.
(a)(b)Table 2
Absorption rate at ten locations in the location diagram.
LocationAbsorption rate10.00,18.000.000034.50,25.000.938143.50,33.000.024345.00,75.501.175448.50,55.501.034050.00,75.501.468456.00,76.501.274365.50,37.000.000079.50,18.000.000098.50,43.500.0000


### 4.3. Data Analysis

According to the calibrated CT system, the receiving information of another unknown medium is obtained. The relevant information of the unknown medium should be solved, such as the geometric shape, the position on the square tray, and the absorption rate at ten positions. The inverse Radon transform method is used to solve the problem, it can know the geometry of the unknown medium, its position in the square tray, and the absorption rate information.Figure7 shows the image reconstruction of the information received by the CT detector through the inverse Radon transformation of the receiving information of the unknown medium. The geometric image reflected is about a certain biological tissue. The ten positions given in the title are marked. In addition, the absorption rate at these ten positions is given in Table 3 (because there is no negative absorption rate, the negative numbers with small deviations are treated as 0).Figure 7
Geometry and its projection heat map. (a) Geometry and position on the square tray. (b) Projection heat map.
(a)(b)Table 3
Absorption rate at the ten positions given in the title.
LocationAbsorption rate10.00,18.000.000034.50,25.002.810343.50,33.006.676545.00,75.500.010248.50,55.500.169750.00,75.502.981756.00,76.506.448965.50,37.000.000079.50,18.006.578598.50,43.500.0321It can be seen from Table3 that 5 positions are inside the medium and 5 positions are outside the medium (the value with an absorption rate less than 0.0001 is noise).


## 4.1. Research Ideas

According to the received information and the projection map, the problem obtains the position, geometry, and absorptivity of the unknown medium on the square tray and specifically gives the absorptivity of 10 positions. Since the known information is the received information on the detector. It is necessary to consider using the Radon transform of CT system image reconstruction to reconstruct the image [9].


## 4.2. Data Processing

The definition of the two-dimensional Radon transform is(4)Rfs,θ=∫−∞+∞∫−∞+∞fx,y∂xcosθ+ysinθ−sdxdy.Among which,fx,y is the function of the object to be reconstructed. Rfs,θ is expressed as the Radon transformation along a certain straight line when the CT system rotates in the angle of view θ. That is, the projection along the straight line, fx,y, denoted as ps,θ. s is the distance from the origin to a certain straight line. Theoretically, the image that can be reconstructed is regarded as the inverse deduction fx,y of the Radon projection transformation function ps,θ, that is, the inverse Radon transformation problem. The formula is(5)fx,y=12π2∫02π∫−∞+∞∂ps,θ/∂sxcosθ+ysinθ−sdsdθ.Using the reception information of an unknown medium obtained by the above-mentioned CT system, we need to use the inverse Radon transform to reconstruct its image.The above geometric shape is the Shepp–Logan head model [10]. In Figure 6(a), the absorption rate of the ten positions given in the title is calibrated. The absorption rate at ten locations in the location diagram are shown in Table 2.Figure 6
The geometric shape of the received information of an unknown medium. (a) Geometry and position on the square pallet. (b) Map of receiving information of unknown medium.
(a)(b)Table 2
Absorption rate at ten locations in the location diagram.
LocationAbsorption rate10.00,18.000.000034.50,25.000.938143.50,33.000.024345.00,75.501.175448.50,55.501.034050.00,75.501.468456.00,76.501.274365.50,37.000.000079.50,18.000.000098.50,43.500.0000


## 4.3. Data Analysis

According to the calibrated CT system, the receiving information of another unknown medium is obtained. The relevant information of the unknown medium should be solved, such as the geometric shape, the position on the square tray, and the absorption rate at ten positions. The inverse Radon transform method is used to solve the problem, it can know the geometry of the unknown medium, its position in the square tray, and the absorption rate information.Figure7 shows the image reconstruction of the information received by the CT detector through the inverse Radon transformation of the receiving information of the unknown medium. The geometric image reflected is about a certain biological tissue. The ten positions given in the title are marked. In addition, the absorption rate at these ten positions is given in Table 3 (because there is no negative absorption rate, the negative numbers with small deviations are treated as 0).Figure 7
Geometry and its projection heat map. (a) Geometry and position on the square tray. (b) Projection heat map.
(a)(b)Table 3
Absorption rate at the ten positions given in the title.
LocationAbsorption rate10.00,18.000.000034.50,25.002.810343.50,33.006.676545.00,75.500.010248.50,55.500.169750.00,75.502.981756.00,76.506.448965.50,37.000.000079.50,18.006.578598.50,43.500.0321It can be seen from Table3 that 5 positions are inside the medium and 5 positions are outside the medium (the value with an absorption rate less than 0.0001 is noise).


## 5. Analysis of Accuracy and Stability of Parameter Calibration

### 5.1. Research Ideas

According to the solution of problem 2, it can be concluded that the most important thing that affects the accuracy and stability of the parameter calibration is the ellipse. Therefore, the focus is on the ellipse and the center of rotation. According to the Canny edge detection algorithm, the center of rotation is solved from another direction and then judged error [11]. Then, design a new template, perform the above operations again, solve the error, and compare with the previous error in order to verify the accuracy and stability of the new template.


### 5.2. Data Processing

To analyze the accuracy and stability of the parameter calibration in question 1, it is mainly to determine the position of the rotation center and its error analysis. Since this is a two-dimensional plane template, we need to find the abscissaδx and ordinate δy of the rotation center. And the problem is that the coordinates of the rotation center obtained in 1 are compared in order to judge their stability and accuracy. The flow chart for calculating the center of rotation is shown in Figure 8.Figure 8
Flowchart for calculating the center of rotation.To solve the center of rotation, the following process will be used to solve. The abscissa and ordinate of the center of rotation are calculated using the centroid shift method (used in question 1), and the abscissa and ordinate of the center of rotation are calculated using the projection edge method.The inverse Radon transform [12, 13] of the Γ data is used to obtain the reconstructed image. In order to better extract the edge of the image and use the edge detection algorithm, we perform the denoising process after grayscale processing in Figure 9.Figure 9
Comparison of the reconstructed image and denoised image. (a) Reconstructed image. (b) Denoised image.
(a)(b)However, the above four operators are traditional edge detection operators such as Log, Sobel, Prewitt, Roberts, and most of them are local gradient operators, as shown in Figure10. These gradient operators are sensitive to noise, so processing images is not very practical.Figure 10
The edge operator extracted by the edge detection method.In 1986, Canny took the lead in proposing the optimal stepped edge detection operator [5], which is the optimal result for continuous jump edge processing interfered by white noise [14]. Canny operator extraction effect of the projection image boundary is better, and it is closer to the real boundary contour.OpenCV can not only perform Canny operator extraction but also find the centroid of the graph, which is the theoretical rotation center position:(6)px,py=287.466pixels,215.869pixels.The position of this point is expressed in the original coordinate system, whereδx¯, δy¯ are axis x and axis y coordinates of the theoretical rotation center, Nt represents the total number of detectors, and ij represents the picture pixels. The calculation formula is(7)δx¯=Nti×i2−pxD,δy¯=Ntj×j2−pyD.The theoretical rotation center coordinates are−9.1764mm,5.1533mm; the rotation center solved from problem 2 is −9.5185mm,5.6560mm. To analyze the stability and accuracy of the two, it is only necessary to analyze the absolute value of the difference between the two in the two directions of the coordinate system:(8)Δδx¯=0.3421mm,Δδy¯=0.5027mm.Calculate its error:Δδx¯/δx¯=3.7280% and Δδy¯/δy¯=9.7549%.Although the error is not very large, the error in they axis direction is close to 10%. Therefore, we need to correct the plane template used in Question 2 so that the error between the theoretical and actual rotation centers is as small as possible.Since the projection heat map of the small circle is a relatively standard chord function, the sphere should be a better template. The geometric shape of the small circle remains unchanged, and the ellipse becomes a standard circle of 15 mm. Assume that the center of rotation is the rotation required in the second problem center, i.e.,−9.5185mm,5.6560mm, as shown in Figure 11. Comparison of the reconstructed image and denoised image is given in Figure 12.Figure 11
Schematic diagram of the new template.Figure 12
Comparison of reconstructed image with a new template and reduced noise image. (a) Reconstructed image. (b) Denoising image.
(a)(b)The following will use the projection edge method to solve the theoretical rotation center of the new template. Firstly, it uses the Radon transform to solve the absorption intensity matrix, and we use the inverse Radon transform on the absorption intensity matrix to obtain the reconstructed image [15]. According to the comparison of the rotation center error of the two, judging the improved accuracy. In order to better extract the edges of the image, we perform the noise reduction after grayscale processing on Figure 11.The following will use OpenCV to process the edge image as is shown in Figure13. The software can not only perform Canny operator extraction [16] but also find out the centroid coordinates px,py=205.597,149.419pixels of the image.Figure 13
OpenCV’s Canny operator extracts the position of the theoretical rotation center of the new template.The position of this point is expressed in the original rectangular coordinate system, whereδx¯′, δy¯′ are the x axis and y axis coordinates of the theoretical rotation center, Nt represents the total number of receivers, and i′j′ represents the picture pixels. The calculation formula is(9)δx¯′=Nti′×i′2−px′D,δy¯′=Ntj′×j′2−py′D.The theoretical rotation center is calculated as−9.1764mm,5.1533mm.From problem 2, we solved the center of rotation:−9.5185mm,5.6560mm. To analyze the stability and accuracy of the two, it needs to analyze the absolute value of the difference between the two in the two directions of the coordinate system:(10)Δδx¯′=1.4818mm,Δδy¯′=0.0000mm.The error is calculated as(11)Δδx¯δx¯=16.1479%,Δδy¯δy¯=0.Although the error on the axis coordinate is a bit larger, the error on the axis is 0. Therefore, changing the two-dimensional error to a one-dimensional error will reduce the impact of the actual error. In summary, change the ellipse to radius, and the standard circle is feasible.


## 5.1. Research Ideas

According to the solution of problem 2, it can be concluded that the most important thing that affects the accuracy and stability of the parameter calibration is the ellipse. Therefore, the focus is on the ellipse and the center of rotation. According to the Canny edge detection algorithm, the center of rotation is solved from another direction and then judged error [11]. Then, design a new template, perform the above operations again, solve the error, and compare with the previous error in order to verify the accuracy and stability of the new template.


## 5.2. Data Processing

To analyze the accuracy and stability of the parameter calibration in question 1, it is mainly to determine the position of the rotation center and its error analysis. Since this is a two-dimensional plane template, we need to find the abscissaδx and ordinate δy of the rotation center. And the problem is that the coordinates of the rotation center obtained in 1 are compared in order to judge their stability and accuracy. The flow chart for calculating the center of rotation is shown in Figure 8.Figure 8
Flowchart for calculating the center of rotation.To solve the center of rotation, the following process will be used to solve. The abscissa and ordinate of the center of rotation are calculated using the centroid shift method (used in question 1), and the abscissa and ordinate of the center of rotation are calculated using the projection edge method.The inverse Radon transform [12, 13] of the Γ data is used to obtain the reconstructed image. In order to better extract the edge of the image and use the edge detection algorithm, we perform the denoising process after grayscale processing in Figure 9.Figure 9
Comparison of the reconstructed image and denoised image. (a) Reconstructed image. (b) Denoised image.
(a)(b)However, the above four operators are traditional edge detection operators such as Log, Sobel, Prewitt, Roberts, and most of them are local gradient operators, as shown in Figure10. These gradient operators are sensitive to noise, so processing images is not very practical.Figure 10
The edge operator extracted by the edge detection method.In 1986, Canny took the lead in proposing the optimal stepped edge detection operator [5], which is the optimal result for continuous jump edge processing interfered by white noise [14]. Canny operator extraction effect of the projection image boundary is better, and it is closer to the real boundary contour.OpenCV can not only perform Canny operator extraction but also find the centroid of the graph, which is the theoretical rotation center position:(6)px,py=287.466pixels,215.869pixels.The position of this point is expressed in the original coordinate system, whereδx¯, δy¯ are axis x and axis y coordinates of the theoretical rotation center, Nt represents the total number of detectors, and ij represents the picture pixels. The calculation formula is(7)δx¯=Nti×i2−pxD,δy¯=Ntj×j2−pyD.The theoretical rotation center coordinates are−9.1764mm,5.1533mm; the rotation center solved from problem 2 is −9.5185mm,5.6560mm. To analyze the stability and accuracy of the two, it is only necessary to analyze the absolute value of the difference between the two in the two directions of the coordinate system:(8)Δδx¯=0.3421mm,Δδy¯=0.5027mm.Calculate its error:Δδx¯/δx¯=3.7280% and Δδy¯/δy¯=9.7549%.Although the error is not very large, the error in they axis direction is close to 10%. Therefore, we need to correct the plane template used in Question 2 so that the error between the theoretical and actual rotation centers is as small as possible.Since the projection heat map of the small circle is a relatively standard chord function, the sphere should be a better template. The geometric shape of the small circle remains unchanged, and the ellipse becomes a standard circle of 15 mm. Assume that the center of rotation is the rotation required in the second problem center, i.e.,−9.5185mm,5.6560mm, as shown in Figure 11. Comparison of the reconstructed image and denoised image is given in Figure 12.Figure 11
Schematic diagram of the new template.Figure 12
Comparison of reconstructed image with a new template and reduced noise image. (a) Reconstructed image. (b) Denoising image.
(a)(b)The following will use the projection edge method to solve the theoretical rotation center of the new template. Firstly, it uses the Radon transform to solve the absorption intensity matrix, and we use the inverse Radon transform on the absorption intensity matrix to obtain the reconstructed image [15]. According to the comparison of the rotation center error of the two, judging the improved accuracy. In order to better extract the edges of the image, we perform the noise reduction after grayscale processing on Figure 11.The following will use OpenCV to process the edge image as is shown in Figure13. The software can not only perform Canny operator extraction [16] but also find out the centroid coordinates px,py=205.597,149.419pixels of the image.Figure 13
OpenCV’s Canny operator extracts the position of the theoretical rotation center of the new template.The position of this point is expressed in the original rectangular coordinate system, whereδx¯′, δy¯′ are the x axis and y axis coordinates of the theoretical rotation center, Nt represents the total number of receivers, and i′j′ represents the picture pixels. The calculation formula is(9)δx¯′=Nti′×i′2−px′D,δy¯′=Ntj′×j′2−py′D.The theoretical rotation center is calculated as−9.1764mm,5.1533mm.From problem 2, we solved the center of rotation:−9.5185mm,5.6560mm. To analyze the stability and accuracy of the two, it needs to analyze the absolute value of the difference between the two in the two directions of the coordinate system:(10)Δδx¯′=1.4818mm,Δδy¯′=0.0000mm.The error is calculated as(11)Δδx¯δx¯=16.1479%,Δδy¯δy¯=0.Although the error on the axis coordinate is a bit larger, the error on the axis is 0. Therefore, changing the two-dimensional error to a one-dimensional error will reduce the impact of the actual error. In summary, change the ellipse to radius, and the standard circle is feasible.


## 6. Conclusion

The accuracy of the solution obtained by this model and the restoration of the reconstructed image are high, and the mechanism analysis and image processing are combined, which can clearly calibrate the parameters of the CT system. Image processing can be used to analyze the accuracy and stability of parameter calibration. Due to the pursuit of verifying the accuracy and stability of the calibration parameters of the CT system, the Canny operator and other image processing methods used are more complicated. For some simpler CT systems, the parameters can be calibrated quickly and effectively using the template to perform the system parameter calibration.



---

*Source: 5360716-2021-08-31.xml*