CN116206003A - Cone beam CT geometric correction method suitable for irregular track - Google Patents

Abstract

The invention discloses a cone beam CT geometric correction method suitable for irregular tracks, which adopts a new optimization constraint iteration method based on a back projection model to solve geometric parameters under each projection view, most of the existing correction methods assume that the internal geometric relationship of a datum point in a phantom is accurate and known, and the problem of whether the phantom precision is suitable for a correction model is seldom and deeply studied. By the scheme of the application, the correction method is a correction method with high noise robustness, and high-precision correction can be realized even when the resolution of a system to be calibrated is twice that of a die body measurement system.

Description

Cone beam CT geometric correction method suitable for irregular track

Technical Field

The invention relates to the technical field of image processing of computer tomography, in particular to a geometric correction method for high-resolution cone beam CT with irregular scanning tracks.

Background

Cone beam CT is currently widely used in many fields, such as medical diagnosis and treatment, preclinical imaging and research, and structural analysis of materials science, among others. In general, different fields of application have different resolution requirements for cone beam CT. Cone-beam CT for clinical, preclinical and materials science has resolutions of 100 microns, 10 microns and 1 micron, respectively. As resolution increases, cone beam CT systems are expected to become more stable and more accurate. Thus, as a prerequisite for obtaining an ideal CT image, the corresponding geometrical correction step becomes more challenging.

In order to solve the correction problems encountered in different scenarios, many correction methods have been proposed in the literature over the last decades. One type is online correction, which does not depend on a particular phantom, but on a loss function of different image quality metrics such as image entropy, image gradients, or a combination of multiple metrics. Such methods are applicable where the acquisition of the geometry is regular but not repeatable. These methods are generally computationally intensive and tend to degrade image quality. Most commercial CT systems are mechanically repeatable and employ another method called off-line correction. The off-line correction method is applicable to both regular and irregular geometry acquisition. For a regular geometric acquisition system, a simple phantom, since the parameters are very few, embedding one or several pellets is sufficient to obtain all geometric parameters. In these correction cases, the model is not complex, and the manual operation is easy, so that they are very simple to use correspondingly. For irregular systems with different geometric parameters of each view, a simple phantom with low precision requirements cannot work. Because the complete projection matrix is not solvable unless some additional conditions are given, such as spatial relationships between point markers or line markers. It has been found by search that most documents assume that the geometric relationship between the markers is known. Some of these rely on specially made phantoms in which fiducial points are precisely combined in a particular pattern, such as a spiral, circle or orthogonal line. The geometry of a particular pattern is considered a known condition, so manufacturing accuracy is a key prerequisite to obtaining accurate geometric parameters. Since manufacturing errors are not easily quantified, none of the prior art methods have focused too much on this problem. In fact, if the system does not have high spatial resolution requirements, it may not be as powerful to make a dedicated phantom with low accuracy requirements. However, as resolution increases, this becomes a challenging task. Many researchers may encounter such situations: i.e. they cannot reproduce the known correction methods or achieve the same correction accuracy based on homemade phantoms. One reason for this is that the accuracy of the phantom does not meet the requirements of the corrected system, but unfortunately, it is not known at all whether this is true, since there is no suitable detection method. Measuring a roughly fabricated phantom appears to be more practical than manufacturing a high precision phantom, and therefore Li et al propose a correction method based on a projection matrix. They do not provide detailed measurement methods nor do they further study how well model measurements match correction models. In theory, the high-precision phantom and the high-robustness correction method are mutually restricted and mutually supported, and are important for correction precision, so that the best correction precision of the irregular track CBCT can be obtained under the condition that the correction model and the measurement of the phantom can be matched. Furthermore, since the fiducial points are not regularly arranged, a correction model based on a projection matrix is currently the best choice. In practical applications, the noise robustness of the model is limited, especially in the case of a low number of fiducial mark points, which requires a more reliable correction method to compensate for this deficiency. Thus, if there is a good phantom measurement method and a more reliable correction model, we can correct a high resolution CT system with irregular trajectories using a relatively coarse phantom that is hand-made. This is very helpful to researchers who want to construct high resolution cone beam CT systems without the need for tip tooling to make high precision phantoms.

Disclosure of Invention

The invention aims to: aiming at the defects of the prior art, the invention provides a cone beam CT geometric correction method suitable for irregular tracks, firstly, the internal geometric relation of a datum point in a phantom is obtained by introducing a measurement model, the measurement model is realized by a high-precision geometric correction model before expansion, and the model is suitable for cone beam CT with circular tracks; a new optimization constraint iteration method based on the back projection model is provided to solve the geometric parameters under each projection view.

The technical scheme is as follows: the invention provides a cone beam CT geometric correction method suitable for irregular tracks, which comprises the following steps of

The steps are as follows:

s1: after space coordinates of the randomly arranged pellets are obtained through calculation, an iterative model is established, and geometric correction of a single projection view is realized;

s2: two coordinate systems are arranged in the model, namely X-Y-Z and W-U-V; wherein W-U-V is referred to as the detector coordinate system, U and V are along the horizontal and vertical axes of the detector, respectively, W being unique to the detector plane; X-Y-Z and is referred to as an object coordinate system;

s3: assuming that the coordinates of the point B in the X-Y-Z and coordinate system are (X ₁ ，y ₁ ，x ₁ ) ^T The method comprises the steps of carrying out a first treatment on the surface of the The corresponding coordinates of the point B in the W-U-V coordinate system can be expressed as

Wherein: t is the translation vector and,

r is the rotation matrix of the matrix of rotation,

r is written as

Wherein eta, theta and->

Is Euler angle;

s4: the x-ray source S coordinate in the model is S (w _s ，u _s ，v _s ) ^T Point B _w The coordinates of the projection point Q of (a) are q= (0, u, v) ^T The method comprises the steps of carrying out a first treatment on the surface of the There are k points in the target coordinate system, called B _k ，K∈[1，K]Their coordinates are given by b _k ＝(x _k ，y _k ，z _k ) ^T The method comprises the steps of carrying out a first treatment on the surface of the The corresponding projection point in the detector coordinate system is Q _k Its coordinate is Q _k ＝(0，u _k ，v _k ) ^T The method comprises the steps of carrying out a first treatment on the surface of the Iterative intermediate solution of geometric parameters

B is given by equation (1) _k Converting into a detector coordinate system to obtain a corresponding point B _w-k Its coordinates are written as B _w-k ＝[R/T]B _k ；

S5: set point B _k To straight line Q _k The foot of the perpendicular line of S is H _k By establishing a straight line Q _k S and B _k H _k The equation set is composed to obtain H _k The coordinate h of (2) _k Deducing B again _k To straight line Q _k The distance of S is as follows:

ds _k ＝||b _w-k -h _k || ₂ ＝||[R|T]b _k -h _k || ₂ (2)

combining all K points, the total distance can be expressed as:

when all geometric parameters are accurate, the corresponding ds=0; the basic objective function can therefore be written as:

further, the cone beam CT geometry correction method applicable to irregular trajectories described above introduces a variance penalty term in equation (4) described above:

wherein the method comprises the steps of

After this term is added, the objective function becomes:

where λ is a weight factor with an empirical value of 10.

Further, in the cone beam CT geometry correction method applicable to irregular trajectories, in most cone beam CT systems, the x-ray source and detector remain relatively stationary during scanning; by solving equation (6), a parameter vector V for each projection view can be obtained _i ＝[s _i ，T _i ，η _i ，φ _i ，θ _i ]，i∈[1，N]Where N is the total number of projections; without noise, all s _i Should be identical; thus, all s _i Average value of (2)

Will be very close to the real coordinates of the x-ray source; usually, there are hundreds of projection views in a system that need to be corrected, so that +.>

Accuracy of (2); after the source position is determined, the objective function of 9 variables is reduced to a function of 6 variables as follows:

the beneficial effects are that: compared with the prior art, the invention has the advantages that: the invention proposes an easy-to-implement geometrical correction scheme that is independent of a phantom with reference points of complex arrangement, while having good robustness to noise and measurement errors. The body model used in the embodiment of the scheme is manufactured manually, so that the requirement on processing equipment and precision is hardly met. The exact coordinates of the pellets were obtained using the NOC method proposed in the published scheme as a reference for subsequent irregular track correction. The simulation in the example shows that the global coordinate error for each pellet is on the order of 0.01 pixel. When gaussian noise with a standard deviation of 0.5 is added, the overall error remains at the same level. This result shows that the NOC method can lay a good foundation for the next step of irregular track correction. The simulation is verified in the simulation, and the result shows that the simulated annealing algorithm can effectively solve the global optimum under the condition that the initial value is far away from the true value. This means that there is no problem in solving this non-male model in practical applications. Iterative methods are generally more flexible and stable than analytical methods. The flexibility of the iterative method can be fully utilized in the study to improve the robustness to noise and measurement errors.

Drawings

FIG. 1 is a projection view geometry of the present invention;

FIG. 2 is a back projection model of the present invention;

FIG. 3 is a convergence curve of the iterative process in example 2;

in FIG. 4, (a) is the projection metric d _fpe-m (b) is a back projection metric d based on GGC method and proposed method _bpe-i Is a distribution of (3);

FIG. 5 is a corrected phantom reconstructed image based on geometric parameters obtained by GGC methods (a 1-a 4) and methods (b 1-b 4) herein;

FIG. 6 is a graph of MTF for a real image and an image obtained based on the GGC correction method and the correction method herein;

FIG. 7 is a reconstructed image based on geometric parameters obtained by the competitive GGC methods (a 1-a 4) and the methods (b 1-b 4) herein;

FIG. 8 is a first low resolution mouse femur reconstruction image in example 4;

fig. 9 is a second high resolution mouse femur reconstruction image in example 4.

Detailed Description

The technical scheme of the invention is described in detail below through the drawings, but the protection scope of the invention is not limited to the embodiments.

Firstly, after the space coordinates of the randomly arranged pellets are obtained, an iterative model is established, and the geometric correction of a single projection view is realized. The geometry of the projection view is shown in fig. 1. There are two coordinate systems, X-Y-Z and W-U-V, respectively. W-U-V is referred to as the detector coordinate system. U and V are along the horizontal and vertical axes of the detector, respectively. W is specific to the detector plane. X-Y-Z and is referred to as the object coordinate system. Assuming that the coordinates of the point B in the X-Y-Z and coordinate system are (X ₁ ，y ₁ ，x ₁ ) ^T . The corresponding coordinates of the point in the W-U-V coordinate system can be expressed as

Where T is the translation vector and R is the rotation matrix. R can be written as

Wherein eta, theta and->

Euler angles, respectively, represent yaw, pitch and roll angles, respectively. Let the x-ray source S coordinate be S (w _s ，u _s ，v _s ) ^T ，B _w The coordinates of the projection point Q are q= (0, u, v) ^T . If all 9 geometrical parameters defining this geometry +.>

Are all known, then point B _w S and Q are collinear. From this observation, the 9 geometric parameters were iteratively obtained by constructing a back projection model as shown in fig. 2 as follows. Suppose there are k points in the target coordinate system, called Bk _， K∈[1，K]Their coordinates are given by b _k ＝(x _k ，y _k ，z _k ) ^T . The corresponding projection point in the detector coordinate system is Q _k Its coordinate is Q _k ＝(0，u _k ，v _k ) ^T . If the iterative intermediate solution of the geometrical parameters is +.>

B is given by equation (1) _k Converting into a detector coordinate system to obtain a corresponding point B _w-k Its coordinates can be written as B _w-k ＝[R/T]B _k . Set point B _k To straight line Q _k The foot of the perpendicular line of S is H _k . By establishing a straight line Q _k S and B _k H _k The system of equations is composed to easily obtain H _k The coordinate h of (2) _k . Deducing B _k To straight line Q _k The distance of S is as follows:

ds _k ＝||b _w-k -h _k || ₂ ＝||[R|T]b _k -h _k || ₂ (2)

combining all K points, the total distance can be expressed as

Theoretically, when all 9 geometric parameters are accurate, the corresponding ds=0. Thus, the basic objective function can be written as:

the iterative method is more flexible than the analytical method because various penalty functions or constraints can be added to improve convergence and reduce the impact of metric errors on the final solution. In the model, the variance penalty term is first introduced as

Wherein the method comprises the steps of

After this term is added, the objective function becomes

Here λ is a weight factor with an empirical value of 10. This variance penalty is a constraint that lets all dss _k And (5) uniformly converging. When there is a measurement error, the objective function value of the optimal parameter is no longer minimum. Therefore, even if the global minimum of the objective function can be obtained, the optimal geometric parameters are not necessarily obtained. To reduce the effect of measurement errors, conventional analysis methods introduce more marker points. Based on experience with the computer vision field paper, about 28 mark points are required. One well-known calibration method proposed by Li et al is to complete the calibration of the chromatographic synthesis system with 44 points. However, in high resolution cone beam CT, fitting so many marker points in a small model is not easy due to the limitations of the field of view. Thus, a new strategy based on the flexibility of iterative methods is introduced, instead of adding marker points to reduce the impact of measurement errors. The flexibility of the iterative approach is that if any variables in the function are known, they can remain unchanged during the iteration process. In calibrating this system, the gantry of the x-ray source and detector can be rotated simultaneously by moving/rotating the object or rotating it, and the relative positions of the x-ray source and detector can be easily maintained. For example, in most cone-beam CT systems, the x-ray source and detector remain relatively stationary during the scan, and thus, by solving equation (6), a parameter vector V for each projection view can be obtained _i ＝[s _i ，T _i ，η _i ，φ _i ，θ _i ]，i∈[1，N]. Where N is the total number of projection views. Theoretically, if there is no noise, all s _i Should be identical. Thus, all s _i Average value of (2)

Will be very close to the real coordinates of the x-ray source. Usually, hundreds of projection views in a system need to be corrected, so that +.>

Is a precision of (a). After the source position is determined, the objective function of 9 variables is reduced to a function of 6 variables as follows:

thus, an approximate solution of s is used to improve the accuracy of the remaining variables. The functions F (s, T, η, Φ, θ) and F (T, η, Φ, θ) are both non-convex functions, so a simulated annealing function algorithm can be used to solve it. The simulated annealing function algorithm combines the simplex simulated annealing method with the traditional simulated annealing method, can obviously avoid sinking into local minima, and improves iteration convergence.

The specific experiments in this example are as follows

Example 1 numerical simulation of phantom measurements

In order to verify the feasibility and robustness of the measurement method, two simulation experiments are designed in the embodiment. The scan trajectory of the simulated medium cone beam CT is circular. The geometric parameters are shown in table 1. The analog detector has 2016 x 1344 pixels with a pixel size of 10 μm. Since the ideal coordinates of the pellets obtained by the measurement model were the original coordinates, multiplied by the magnification of the measurement system, both SRD and SDD were set to 8cm in order to intuitively verify the accuracy of the measurement model, and thus the magnification was 1.0. The model was simulated as a cylinder with 8 pellets arranged in a spiral around the circumference. In the first simulation, noiseless analytical projections were performed on the pellets to obtain projection data of the pellets. A total of 360 projection data is obtained over 360 °. The accuracy of the sphere estimated spatial coordinates was evaluated using the following indices:

where K is the number of pellets, d (B _{i_real} ，B _{i_est} ) Is the distance between the true and estimated coordinates of the i-th sphere. In the second simulation, the geometric parameters are the same as in the first simulation. The diameter of the sphere was set to 400 μm, simulating a real projection environment, and in addition, 360 projections were collected over 360 °.

Table 1: geometric parameters of simulation experiments

Table 2: real three-dimensional coordinates and actual three-dimensional coordinates of the pellets in the noiseless experiment. X, Y, Z represent the true values in three directions, respectively. X is X _est 、Y _est 、Z _est The measured values in three directions are respectively obtained. The unit is a pixel.

The centroid of the pellet projection can be extracted by a threshold segmentation method, and Gaussian noise with standard deviation of 0.5 pixel is added into the projection centroid for verifying the robustness of the measurement method. The whole simulation process is repeated for 8 times, and a final three-dimensional space coordinate is obtained by using a disclosed measurement method through a random movement model. The results of the first simulation are shown in table 2, where the first three behavior pellets have real values of spatial coordinates in three directions. The units in the table are pixels. From the fourth line to the sixth line are estimated coordinates. It can be seen that the estimated spatial coordinates of the pellets have very high accuracy, the errors in all three directions are less than 0.01 pixel, and the overall error accuracy DE is 0.013 pixel. The results of the second simulation run with 50% gaussian noise are shown in table 3. The first three rows of parameters are the estimated coordinates of the pellet in three directions. The last three rows are errors of pellets. It can be seen that most of the errors in the three directions are much smaller than 0.1 pixel, with an overall error DE of 0.039 pixel. Based on the experimental results, the phantom measurement method has good robustness to noise and great practical application potential.

Table 3: the three-dimensional coordinates of the pellets and the corresponding errors measured in the noise experiments. X is X _est 、Y _est 、Z _est The measured values in three directions are respectively obtained. X is X _error ，Y _error And Z _error Errors in three directions, respectively, are in pixels.

Example 2 Convergence of iterative method

In order to verify the convergence of the proposed geometric correction method, a simulation experiment is performed by using the phantom and arbitrarily setting a geometric parameter and a parameter value interval. The analog detector has a size of 2016 x 1344 pixels and a pixel pitch of 10 μm. The geometrical relationship between the phantom coordinate system X-Y-Z and the detector coordinate system W-U-V under a specific projection view is shown in the first row of Table 4. With the standard three-dimensional space coordinates of 8 pellets as a known condition, a non-convex objective function, equation (6), is constructed and then iteratively solved by using a simulated annealing algorithm. The simulation experiment is only used for verifying the convergence performance of the iterative algorithm, so that the measurement error is not considered in the simulation. The initial values of the experimental settings are located in the second row of table 4, which are set far enough from the true values to demonstrate the convergence performance of the algorithm. The resulting errors for the 9 geometric parameters are shown in the last line of table 4. It can be seen that the error is almost equal to 0. The convergence curve of the iterative process is shown in fig. 3. It can be seen that using the simulated annealing algorithm allows the objective function to converge effectively to a global minimum even if the initial value is far from the true value.

Table 4: true values, initial values and errors of geometric parameters of scan views

Example 3 geometric correction of Cone Beam CT System simulating non-circular trajectories

To verify the robustness of the proposed iterative correction method, this embodiment simulates a CBCT system with an elliptical trajectory. The system SDD is 8cm, the detector has 2016 x 1344 pixels, and the pixel spacing is 20 μm. Elliptic trajectory equation is

Eta and->

The rotation angles of the respective projection views are 2.0+k degrees, k ε [0, 359) (k is an integer) respectively 1 DEG and 3 deg. The phantom used was that of the above simulation experiments, containing 8 pellets with a diameter of 400 μm. The true coordinates of the pellets are shown in the first three rows of table 2. In this simulation, geometric correction was performed using the coordinates (from line 4 to line 6) measured in table 2. To further verify the noise robustness of the proposed correction method, gaussian noise with a standard deviation of 0.5 was added to the resulting measured coordinates. In order to quantitatively evaluate the correction accuracy of the algorithm, we simulated a test phantom containing 360 uniformly distributed marker points, and simulated the phantom, and we introduced two indexes of average forward projection error and average back projection error. The average forward projection error is expressed as follows:

wherein Px is _{i_m} For an ideal projection point of the spatial point xi on the mth projection view,

for which projection points are estimated, KAnd m is the number of projection views. The average backprojection error is:

here the number of the elements is the number,

is test point x _i Is described. />

Can be obtained by solving the following loss function:

wherein d (x, r) _mi ) Is the back projection ray r passing through the ith test point under the projection view from the space point x to the mth _mi Is a euclidean distance of (c). Metric index d _bpe The degree of divergence of the back-projected light is reflected, and the magnitude of the value directly affects the quality of the reconstructed image.

In this embodiment, a well-known analysis correction algorithm GGC method is used as a reference comparison method, which indicates that larger noise can reduce correction accuracy, and then the quality of the reconstructed image can be reduced, and the method proposed in the application can reduce the influence of noise. And scanning and reconstructing the calibration phantom and one volume data by using two groups of geometric parameters respectively obtained by a comparison method and the proposed method. The reconstruction method used is the general FDK algorithm provided by astm a Toolbox. The resolution of the reconstructed phantom image was quantitatively analyzed using MTF (modulation transfer function). D for both methods _fpe-m And d _bpe-i The distribution of (2) is shown in figure 4. It can be seen that the forward projection error and the back projection error d of the proposed method _fpe-m And d _bpe-i The value distribution of (2) is very concentrated and the average forward projection error and the average back projection error obtained by the method of the application are much smaller than those obtained by the comparison method. This isSome quantitative analysis results prove that the proposed iterative correction model can better suppress noise and measurement errors, and obtain more accurate geometric parameters. Since the back projection error is directly related to the quality of the reconstructed image, it can be inferred that the reconstructed image obtained based on the geometric parameters obtained by the method has better quality and higher resolution. The reconstructed images of the corrected phantom obtained based on the two methods are shown in fig. 5, wherein the first to third column images are axial position images, and the fourth column image is sagittal position image. The 2 nd and 4 th line images are enlarged images of the square mark region in the 1 st and 3 rd line images, respectively. (a) The line images are obtained based on the obtained geometrical parameters of the contrast method and (b) the line images are reconstructed based on the geometrical parameters obtained by the proposed method. It is apparent that the edges of the pictures in row (b) are more clear. Since tomographic images of the sphere in the phantom look like planar discs, the MTF curve is easily obtained by extracting the boundaries of these discs, fitting the Edge Spread Function (ESF) with a gaussian error function. Three MTF curves are shown in fig. 6, from reconstructed images of the real image, the geometric parameters obtained based on the method of the present application and the comparative GGC method, respectively. It can be seen that the MTF curve of the reconstructed image obtained based on the correction method of the present application is closer to the real image, and this result also accords with the visual observation result. These results all show that the correction method of the present application can obtain more accurate geometric parameters in the case of larger noise. The result of the comparison of the reconstructed images based on the volume data of the two sets of geometric parameters is shown in fig. 7, (a) the line images are reconstructed based on the geometric parameters obtained by the comparison-based GGC method, and (b) the line images are obtained based on the geometric parameters obtained by the proposed method. These images contain many fine structural features. From the magnified region of the image, it can be seen that the image reconstructed by the geometric parameters based on the correction method of the present application has a clearer boundary. (marked with arrows).

Example 4 validation in high resolution Micro-CT based on gantry Structure

Currently there are generally two types of scanning devices for Micro-CT. The first is based on sample rotation, in which the sample is rotated while the source and detector remain stationary. The device generally has better spatial resolution and is commonly used for scanning of ex vivo samples. Thus, it is commonly referred to as ex vivo Micro-CT. The second is based on gantry rotation, in which the sample is stationary and the source and detector mounted on the gantry are rotating. The advantage of this device is that it can be used to scan living small animals. Such devices are commonly referred to as in vivo Micro-CT. The disadvantage of this setup is that it has a lower spatial resolution compared to ex vivo Micro-CT. One of the main reasons is that in ex vivo Micro-CT, the circularity and eccentricity of the gantry is more difficult to guarantee due to the fact that the gantry equipped with the x-ray source and detector is too heavy, especially in case of high resolution scanning. In general, for Micro-CT ex vivo with a resolution on the order of microns, the circularity of the trajectory is relatively easy to guarantee, so that the ex vivo Micro-CT can be used to measure the hand-made phantom. In this experiment, the body of the phantom was made of polymethyl methacrylate, containing 8 metal pellets of about 400 μm diameter helically embedded in the body. We performed 8 scans of the phantom using ex vivo Micro-CT with a resolution of 9 μm. During each scan, the position of the phantom is randomly changed to uniformly obtain 360 projections over a 360 range. The three-dimensional coordinates of the 8 groups of balls obtained through measurement can be subjected to a registration procedure to obtain final three-dimensional coordinates of the balls. The mean and variance of the three-dimensional coordinates of the pellets are shown in table 5, and it can be seen from the table that the variances of the three components of the coordinates of the pellets are all at the sub-pixel level, which means that in a practical system, the measurement results are stable and the measurement accuracy is high.

Table 5: mean and variance of three-dimensional coordinates of 8 pellets in real experiment

After obtaining the coordinate estimation of the sphere in the model, we calibrate a high-resolution Micro-CT based on gantry rotation (living body) by using the calibration model of the application to verify the accuracy of the coordinate estimation and the performance of the proposed correction method in practical application. The Micro-CT was equipped with an x-ray tube with a focus of 5 μm and a CCD camera with dimensions 4008 x 2672 pixels with a pixel pitch of 9 μm. The x-ray source and the CCD camera are arranged at two ends of the portal frame, the distance between the x-ray source and the CCD camera is 20cm, and the magnification is set to be about 2. The operating voltage was 60kV and the operating current was 130uA. CCD cameras have two binding modes, 1×1 and 2×2. In the case of binning2×2, the resolution of the system is about 9 μm, which is essentially the same as the in vitro Micro-CT resolution used to measure the calibration phantom. The first experiment below is based on this binding mode. The calibration phantom was placed on a scanning bed to obtain 400 projections over a 360 range. The Micro-CT orbit is no longer a standard circle due to the large load and high resolution. To demonstrate this we use the NOC method. The NOC method is a high-precision calibration method suitable for geometric calibration of circular scanning tracks. The classical GGC method used in the simulation experiment is also used for intuitively displaying the measurement precision of the small ball, and meanwhile, the method can also be used as a comparison method for comparing the application performance in the actual scene. After the geometrical parameters are obtained by the three correction methods respectively, 720-angle mouse femur projection data are acquired for reconstructing a tomographic image. Assuming that the track is round, we use the traditional FDK reconstruction method, and the other two cases use the general FDK algorithm provided by ASTRA Toolbox to reconstruct the image. In a second experiment, the binding mode was set to 1 x 1, where the Micro CT resolution was 4.5 μm, which is doubled over the ex vivo Micro-CT resolution used to measure the coordinates of the beads, where the measurement errors of the beads would be significantly amplified, and these increased errors would have a greater negative impact on the calibration accuracy, so that the 1 x 1 mode could be used to further test the performance of the geometry correction method in the extreme case of more severe noise or measurement errors. In the second experiment, we still adopted the classical GGC method as a comparison method, the samples tested were mouse femur samples, and the total of 720 projections, and the reconstruction method was consistent with the previous experiment. The reconstructed image of the femur of the mice from the first experiment is shown in FIG. 8 with a resolution of 9. Mu.m. In the figures, the images of columns (a), (b) and (c) are obtained based on the geometric parameters obtained by the NOC, GGC and the method of the present application, respectively. Each column contains an axial image and a sagittal image. It can be seen from the first column that the image is very blurred and has many ghosts and artifacts. Since the NOC method is a high-precision geometric correction method suitable for circular scan trajectories, the only reason for these artifacts and blurring is that the trajectory is not a standard circle. The second column can be used for seeing better image quality, no obvious artifact and combining the conclusion obtained by the first column of images, the system track is further proved to be not round, the measured small ball coordinates are accurate, and the GGC method can be guaranteed to achieve a relatively good correction result on the basis of 8 small balls. By comparing the images of the second and third columns, it can be found that the image quality of the two columns is similar. This result shows that the correction method of the present application has performance comparable to the comparative GGC method when the resolution of the phantom measurement system is similar to that of the system to be corrected. In the magnified areas selected from the reconstructed images of columns 2 and 3, some subtle differences in detail can be seen in the two images, as can the sharper edges of the image from the magnified areas of the third column. This seems to indicate that the correction method of the present application has a higher calibration accuracy, whereas the reason for the insignificant performance improvement may be that the measurement of the pellets is sufficiently accurate that the calibration accuracy of the two methods is comparable. In a second experiment, the measured pellets were used for calibration of Micro-CT with a resolution of 4.5 μm, where the measurement error was at least twice as great as before. The reconstructed image of the femur of the mice from the second experiment is shown in fig. 9. The first and second line images are the reconstruction results of the geometrical parameters obtained by the contrast-based GGC method and the correction method of the present application, respectively, from which it can be seen intuitively that these images have a higher resolution than the images obtained in the first experiment, and that the differences between the first and second line images are very pronounced, i.e. the picture of the first line has a lot of artifacts, but the second line hardly sees any artifacts. Furthermore, the details of the second line of pictures are very clear. The resolution of the system to be corrected is known to be much higher than that of the phantom measurement system, so that the noise is increased when the centroid of the projection of the small sphere obtained in the system to be corrected is extracted, and the measurement error of the small sphere is at least 2 times that of the original sphere, but the geometric parameters which are accurate enough can still be obtained by using the proposed correction method to generate a high-quality reconstructed image. These results demonstrate that the proposed iterative correction model has better ability to suppress noise and measurement errors than the classical approach.

As described above, although the present invention has been shown and described with reference to certain preferred embodiments, it is not to be construed as limiting the invention itself. Various changes in form and details may be made therein without departing from the spirit and scope of the invention as defined by the appended claims.

Claims (3)

1. The cone beam CT geometric correction method suitable for the irregular track is characterized by comprising the following steps of:

s1: after space coordinates of the randomly arranged pellets are obtained through calculation, an iterative model is established, and geometric correction of a single projection view is realized;

s2: two coordinate systems are arranged in the model, namely X-Y-Z and W-U-V; wherein W-U-V is referred to as the detector coordinate system, U and V are along the horizontal and vertical axes of the detector, respectively, W being unique to the detector plane; X-Y-Z and is referred to as an object coordinate system;

s3: assuming that the coordinates of the point B in the X-Y-Z and coordinate system are (X ₁ ,y ₁ ,x ₁ ) ^T The method comprises the steps of carrying out a first treatment on the surface of the The corresponding coordinates of the point B in the W-U-V coordinate system can be expressed as

Wherein: t is the translation vector and,

r is the rotation matrix of the matrix of rotation,

r is written as

Wherein eta, theta and->

Euler angles, respectively represent yaw, pitch and roll angles;

s4: the x-ray source S coordinate in the model is S (w _s ,u _s ,v _s ) ^T Point B _w The coordinates of the projection point Q of (a) are q= (0, u, v) ^T The method comprises the steps of carrying out a first treatment on the surface of the There are k points in the target coordinate system, called B _k ,K∈[1,K]Their coordinates are given by b _k ＝(x _k ,y _k ,z _k ) ^T The method comprises the steps of carrying out a first treatment on the surface of the The corresponding projection point in the detector coordinate system is Q _k Its coordinate is Q _k ＝(0,u _k ,v _k ) ^T The method comprises the steps of carrying out a first treatment on the surface of the Iterative intermediate solution of geometric parameters

B is given by equation (1) _k Converting into a detector coordinate system to obtain a corresponding point B _w-k Its coordinates are written as B _w-k ＝[R|T]B _k ；

S5: set point B _k To straight line Q _k The foot of the perpendicular line of S is H _k By establishing a straight line Q _k S and B _k H _k The equation set is composed to obtain H _k The coordinate h of (2) _k Deducing B again _k To straight line Q _k The distance of S is as follows:

ds _k ＝||b _w-k -h _k || ₂ ＝||[R|T]b _k -h _k || ₂ (2)

combining all K points, the total distance can be expressed as:

when all geometric parameters are accurate, the corresponding ds=0; the basic objective function can therefore be written as:

2. the cone beam CT geometry correction method for irregular trajectories according to claim 1, wherein:

introducing a variance penalty term in equation (4):

wherein the method comprises the steps of

After this term is added, the objective function becomes:

where λ is a weight factor with an empirical value of 10.

3. The cone beam CT geometry correction method for irregular trajectories according to claim 1, wherein: in most cone-beam CT systems, the x-ray source and detector remain relatively stationary during scanning; by solving equation (6), a parameter vector V for each projection view can be obtained _i ＝[s _i ,T _i ,η _i ,φ _i ,θ _i ],Y∈[1,N]N is the total number of projection views; without noise, all s _i Should be identical; thus all s _i Average value of (2)

Will be very close to the real coordinates of the x-ray source; usually, hundreds of projection views in a system need to be corrected, so that +.>

Accuracy of (2); after the source position is determined, the objective function of 9 variables is reduced to a function of 6 variables as follows:

/>

Images