Method for correcting hardening artifacts in ct images at a single voltage based on integral invariance

Abstract

The application discloses a CT image hardening artifact correction method based on integral invariance under single voltage, which comprises the following steps: collecting X-ray projection images of CT scanning under single voltage, dividing by the background value to obtain an X-ray imaging model, and writing the X-ray imaging model in a matrix form; obtaining an X-ray projection image decomposition model according to the matrix form of the X-ray imaging model; adding an integral invariance equation in the projection image decomposition model to constrain narrow energy spectrum projection obtained by projection decomposition under single voltage, and realizing CT image hardening artifact correction under single voltage; and the application can realize hardening artifact correction for the imaging condition under single voltage of a conventional CT system.

Description

Technical Field

[0001] The present invention relates to the field of CT imaging technology, and in particular to a method for correcting CT image hardening artifacts under a single voltage based on integral invariance. Background Art

[0002] The X-ray source of a conventional CT system emits a continuous, broad energy spectrum. When passing through an object, low-energy rays attenuate more rapidly than high-energy rays, resulting in a higher average energy in the X-ray beam. This leads to hardening artifacts in the reconstructed CT image. These artifacts affect not only the material but also the pixel's location, manifesting as "same object, different images" and "same image, different objects." Methods for correcting these hardening artifacts are generally categorized into hardware and software. Hardware correction methods require additional hardware or modifications to existing system hardware, increasing imaging costs. Common software correction methods include multi-energy (dual-energy) imaging correction, linearization, and statistical iterative reconstruction. Multi-energy CT imaging requires the acquisition of projection data at multiple energies, requiring multiple scans or the use of more advanced dual-energy detectors or photon-counting detectors. Multiple scans reduce imaging efficiency, and more advanced detectors require higher hardware costs. Linearization methods use polynomial fitting to determine the relationship between material thickness and projection values ​​across the continuous energy spectrum. This correction is performed on the continuous spectrum projections. This requires the production of additional samples of a specific shape to collect data for the polynomial coefficients. Statistical iterative reconstruction methods require knowledge of the X-ray energy spectrum at the current voltage, which requires additional experimental measurements and is slow to reconstruct. Therefore, existing methods for correcting CT image hardening artifacts using projection decomposition require either collecting projection data at multiple voltages (reducing imaging efficiency) or modifying hardware conditions to directly acquire projection data at multiple energies (using higher-performance hardware, increasing imaging costs). This makes them inapplicable to conventional CT systems imaging at a single voltage. Summary of the Invention

[0003] In order to solve the problems existing in the prior art, the purpose of the present invention is to provide a method for correcting hardening artifacts in CT images under a single voltage based on integral invariance. The present invention can realize hardening artifact correction for the imaging situation of a conventional CT system under a single voltage.

[0004] To achieve the above object, the present invention adopts a technical solution: a method for correcting hardening artifacts in CT images under a single voltage based on integral invariance, comprising the following steps:

[0005] Step 1: Acquire an X-ray projection image of a CT scan under a single voltage, divide it by its background value, obtain an X-ray imaging model, and write it into a matrix form;

[0006] Step 2: Obtain an X-ray projection image decomposition model according to the matrix form of the X-ray imaging model;

[0007] Step 3: Add the integral invariance equation to the projection image decomposition model to constrain the narrow energy spectrum projection obtained by projection decomposition under a single voltage, so as to realize the correction of CT image hardening artifacts under a single voltage.

[0008] As a further improvement of the present invention, the step 1 is specifically as follows:

[0009] The X-ray tube emits continuous energy spectrum X-rays, and the imaging model is written as:

[0010] I=I0∫S(E)exp[-∫ L u(x,E)dx]dE

[0011] Where: I0 represents the initial intensity of the ray, I is the intensity after ray attenuation, which is equivalent to the background grayscale and image grayscale in actual imaging, L represents the path passed by the X-ray, u(x, E) represents the attenuation coefficient of the ray with energy E at point x on the path L, and S(E) represents the normalized equivalent energy spectrum, which satisfies:

[0012] 1=∫S(E)dE

[0013] The imaging model is discretized into a summation, and both sides of the equation are divided by I0:

[0014]

[0015] Where: E r is the energy of the rth narrow energy spectrum, S(E r ) is the weighting coefficient of the rth narrow energy spectrum image, K is the type of material in the object, u rk For the kth material in E r The attenuation coefficient, d k is the length of the kth material on the path of X-rays; where S(E r )satisfy:

[0016]

[0017] Assume that the number of X-ray projection pixels at the jth scanning angle is M j (j=1, 2, ..., J), the sum of the number of projected image pixels under all scanning angles is recorded as M, then:

[0018]

[0019] The I / I0 corresponding to the mth pixel of the X-ray projection image is denoted as f m , m=1,2,…,M; then the X-ray imaging model is written in matrix form as:

[0020] F=Sexp(-UD)

[0021] Where: F=(f1,f2,……,f M ), S is the narrow energy spectrum image weighting coefficient that constitutes each voltage X-ray image, that is, s r =S(E r ), a column of U corresponds to the attenuation coefficient of a material in each narrow spectral band, and the corresponding row of D is the attenuation length of the same material at each pixel.

[0022] As a further improvement of the present invention, the step 2 is specifically as follows:

[0023] The volume expression of the kth material at the jth angle is:

[0024]

[0025] This value has nothing to do with the angle and is determined by the material distribution of the object itself, so:

[0026] T kj =T k(j+1) ,j=1,2,,J-1,k=1,2,...,K

[0027] Get the X-ray projection image decomposition model:

[0028]

[0029]

[0030] As a further improvement of the present invention, the integral invariance equation in step 3 is as follows:

[0031]

[0032] As a further improvement of the present invention, the step 3 is specifically as follows:

[0033] Set the energy spectrum value R(R) and record the maximum ray energy of the image as E max , with E max / R is the energy spectrum interval to divide the continuous energy spectrum to obtain the narrow energy spectrum interval [E' r-1 ,E' r ], r=1,2,……,R, and E'0=0, E' R =E max ; The middle energy of the rth narrow energy spectrum is recorded as E r , that is, E r =(E' r-1 +E' r ) / 2;

[0034] Denote the material type of the imaging object as K. By querying the common material attenuation coefficient database, the attenuation coefficients of each material at each energy can be obtained. Denote the attenuation coefficient of the k-th material at the X-ray energy E r as u rk , and denote U=(u rk ) RK ;

[0035] Initialize S. Initialize the assignment of S as s1 = s2 =... = s R = 1 / R;

[0036] Initialize D, D=(d km ) KM , where d km is the length of the k-th material on the X-ray penetration path corresponding to the m-th projection pixel point. Denote the maximum value of the distance between any two points on the imaging object as d max , and generate a random matrix with the same size as D as the initial value of D, and the value range of each element is between 0 and d max ;

[0037] Set the maximum number of iterations Niter, the reference relative error ε for stopping the iteration, and denote the number of iteration rounds as niter; Calculate the objective function value under the current variables denoted as h0, and initialize niter = 0;

[0038] Let niter = niter + 1, update S iteratively according to the formula , and perform row normalization on S; Update D iteratively according to the formula ; denotes the Hadamard product of matrices, that is, multiplying corresponding elements, and - denotes the Hadamard division of matrices, that is, dividing corresponding elements;

[0039] Update all d according to the integral invariance equation km ;

[0040] Calculate the objective function value under the current variables denoted as h niter ;

[0041] Calculate the absolute value of (h niter -h (niter-1) ) / h niter . If it is less than ε, the iteration ends and proceed to the next step. Otherwise, judge whether niter is less than Niter. If niter < Niter, continue the iterative update. If niter ≥ Niter, proceed to the next step;

[0042] The current UD is calculated, where each line is a narrow energy spectrum projection. The narrow energy spectrum projection is reconstructed to obtain a CT image corrected for hardening artifacts.

[0043] The beneficial effects of the present invention are:

[0044] This method adds an integral invariance equation to the projection decomposition model, constrains the narrow energy spectrum projection obtained by projection decomposition under a single voltage, realizes CT hardening artifact correction under a single voltage, and expects to obtain a CT image with hardening artifact correction. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 Directly reconstructed image and grayscale change map of red line position under conventional CT system;

[0046] Figure 2 This is the reconstructed image after hardening artifact correction and the grayscale change diagram of the red line position in an embodiment of the present invention. DETAILED DESCRIPTION

[0047] The embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0048] Example

[0049] The hardening artifact correction method is to collect projections under multiple voltages, decompose them, obtain narrow energy spectrum projections, reconstruct narrow energy spectrum images, and achieve hardening artifact correction. The following describes this type of method using the method in the literature [Wei Jiaotong, Chen Ping, Han Yan. High-contrast CT imaging by multi-voltage X-ray image decomposition [J]. Opto-Electronic Engineering, 2016, 43(08): 59-63.] as an example.

[0050] The X-ray tube emits continuous energy spectrum X-rays, and the imaging model can be written as:

[0051] I=I0∫S(E)exp[-∫ L u(x,E)dx]dE (1)

[0052] Where: I0 represents the initial intensity of the ray, I is the intensity after ray attenuation, which is equivalent to the background grayscale and image grayscale in actual imaging, L represents the path passed by the X-ray, u(x, E) represents the attenuation coefficient of the ray with energy E at point x on the path L, and S(E) represents the normalized equivalent energy spectrum, which satisfies:

[0053] 1=∫S(E)dE (2)

[0054] It is related to the X-ray initial energy spectrum, the detector scintillator, etc. In formula (1), the integral is a linear operation, and the discrete operation becomes a summation. Both sides of the equation are divided by I0:

[0055]

[0056] Where: E r is the energy of the rth narrow energy spectrum, S(E r ) is the weighting coefficient of the rth narrow energy spectrum image, K is the type of material in the object, u rk For the kth material in E r The attenuation coefficient, d k is the length of the kth material on the path of X-ray passing. r )satisfy:

[0057]

[0058] The I / I0 corresponding to the mth pixel of the X-ray image at the nth voltage is denoted as f nm , n=1,2,…,N,m=1,2,…,M. According to formula (3), the X-ray imaging of N voltages can be written in matrix form as

[0059] F=Sexp(-UD) (5)

[0060] Where: F=(f nm ) NM ,S=(s nr ) NR ,U=(u rk ) RK , D=(d km ) KM The rows of S are the weighting coefficients for the narrow energy spectrum images that constitute the X-ray images at each voltage (after dividing by the background). Each column of U corresponds to the attenuation coefficient of a material in each narrow spectrum band, and the corresponding row of D is the attenuation length of the same material at each pixel. An object consisting of a uniform mixture of multiple components is considered a single material. Each row of UD is a narrow energy spectrum projection, which is the final desired projection.

[0061] Considering the errors in actual imaging, the solution of Equation (5) is transformed into an optimization model with the minimum sum of squared errors as the objective function, as follows:

[0062]

[0063] The iterative solution formula is

[0064]

[0065]

[0066]

[0067] in: The rows of S are normalized in each iteration.

[0068] The overall solution steps are as follows:

[0069] (1) Initialize S, U, D, set the R value, and set the iteration stopping criterion: the relative difference between the objective function values ​​of two adjacent iterations is less than ε;

[0070] (2) Iteratively update S according to formula (7);

[0071] (3) Normalize the rows of S;

[0072] (4) Iteratively update the rows of U according to formula (8);

[0073] (5) Normalize the columns of U;

[0074] (6) Iteratively update D according to formula (9);

[0075] (7) When the objective function value meets the iteration stopping criterion, the iteration is stopped.

[0076] The X-ray projection image under the continuous energy spectrum is decomposed to obtain narrow energy spectrum projections, and the narrow energy spectrum CT image is reconstructed to achieve hardening artifact correction. The essence is to solve formula (5). The reason why it is necessary to collect X-ray projection data under multiple voltages is that the number of unknown quantities contained in formula (5) is less than the number of equations. To reduce the number of voltages, other constraints or equations must be added. In CT imaging scanning, the object itself is unchanged, and the distribution of each material is also unchanged. Therefore, the volume of each material in the object at different angles is unchanged (in two-dimensional CT imaging, the area is unchanged). When imaging under parallel beams, the integral of the path passed by the X-ray is the object volume (in two-dimensional CT imaging, the area). Under non-parallel beams, the integral value is close to the object volume (in two-dimensional CT imaging, the integral value is close to the object area). This integral value is independent of the angle and is determined by the object itself. Therefore, this embodiment adds an integral invariance equation to the projection decomposition to constrain the narrow energy spectrum projection obtained by projection decomposition under a single voltage, thereby achieving CT hardening artifact correction under a single voltage.

[0077] Assume that the number of X-ray projection pixels at the jth scanning angle is M j (j=1, 2, ..., J), the sum of the number of projected image pixels under all scanning angles is recorded as M, then

[0078]

[0079] The I / I0 corresponding to the mth pixel of the X-ray projection image is denoted as f m , m=1,2,…,M. According to formula (3), the X-ray imaging model is written in matrix form as

[0080] F=Sexp(-UD) (11)

[0081] Where: F=(f1,f2,……,f M ), S=(s1,s2,……,s R ), U=(u rk ) RK , D=(d km ) KM S is the weighting coefficient of the narrow energy spectrum image that constitutes each voltage X-ray image (after dividing by the background), that is,

[0082] s r =S(E r ) (12)

[0083] Each column of U corresponds to the attenuation coefficient of a material in each narrow spectral band, and the corresponding row of D is the attenuation length of the same material at each pixel. A uniform mixture of multiple components in an object is considered a single material. Each row of UD represents a narrow energy spectrum projection, and a hardening-corrected CT image is reconstructed from each narrow energy spectrum projection.

[0084] The volume (area) expression of the kth material at the jth angle is:

[0085]

[0086] This value has nothing to do with the angle and is determined by the material distribution of the object itself, so

[0087] T kj =T k(j+1) ,j=1,2,...,J-1,k=1,2,...,K (14)

[0088] Get the X-ray projection image decomposition model

[0089]

[0090]

[0091] In this model, the energy spectrum is divided into fixed segments, and each energy segment uses the intermediate energy as a reference to determine the U value. Thus, in formula (15), F and U are known quantities, and S and D are quantities to be solved. After obtaining D, UD is calculated. Each line of UD is a narrow energy spectrum projection, which is the final projection required, and a narrow energy spectrum image is obtained by reconstructing it. In the iterative solution of this model, the integral invariance constraint is first ignored, that is, in each round of iteration, S and D are first iterated using formulas (7) and (9). Here, S and F have only one row, which is a special case of n = 1 in formulas (7) and (9). Then the value of D is forced to be redistributed so that it satisfies the constraint of formula (14) as much as possible, specifically:

[0092]

[0093] In summary, the specific steps of this embodiment are:

[0094] (1) Collect the X-ray projection image of the CT scan under a single voltage and divide it by the background value to obtain F in formula (11). The angle number of the collected projection is denoted as j, j = 1, 2, ..., J. The number of pixels of the projection image at the jth angle is denoted as M j , and remember

[0095] (2) Set the energy spectrum value R(R). The maximum ray energy of the image is E max , with E max / R is the energy spectrum interval to divide the continuous energy spectrum to obtain the narrow energy spectrum interval [E' r-1 ,E' r ], r=1,2,……,R, and E'0=0, E' R =E max The middle energy of the rth narrow energy spectrum is recorded as E r , that is, E r =(E' r-1 +E' r ) / 2;

[0096] (3) The material type of the imaging object is recorded as K. The attenuation coefficient of each material at each energy can be obtained by querying the public material attenuation coefficient database. r The attenuation coefficient is denoted as u rk . Remember U=(u rk ) RK ;

[0097] (4) Initialize S. The initial value of S is s1=s2=……=s R =1 / R;

[0098] (5) Initialize D. D=(d km ) KM , where d km The length of the kth material on the path of the X-ray corresponding to the mth projection pixel. The maximum distance between any two points on the imaging object is recorded as d max Generate a random matrix with the same size as D as the initial value of D, and each element value ranges from 0 to d max between;

[0099] (6) Set the maximum number of iterations Niter, the reference relative error ε for stopping the iteration, and the number of iterations is recorded as niter. Calculate the objective function value under the current variable Denote it as h0. Initialize niter=0;

[0100] (7) Let niter = niter + 1. Iteratively update S according to formula (7) and perform row normalization on S. Iteratively update D according to formula (9);

[0101] (8) Update all d according to formula (16) km ;

[0102] (9) Calculate the objective function value under the current variables Denote it as h niter ;

[0103] (10) Calculate the absolute value of (h niter - h (niter-1) ) / h niter . If it is less than ε, end the iteration and go to the next step. Otherwise, judge whether niter is less than Niter. If niter < Niter, go to step (7). If niter ≥ Niter, go to the next step;

[0104] (11) Calculate the current UD, each row of which is a projection of a narrow energy spectrum segment, and reconstruct the CT image after hardening artifact correction from the projection of the narrow energy spectrum segment.

[0105] A certain sample is composed of two materials, titanium alloy and steel. The outside is titanium alloy and the internal circular area is steel. The directly reconstructed image and the gray level change at the red line position under the conventional CT system are as Figure 1 shown. The expected correction result of the inventive method is as Figure 2 shown. Set R = 16. In the directly reconstructed image, hardening artifacts are obvious, especially in the titanium alloy area. The gray level values of the pixels at the image edge are significantly higher than those in the middle part. After hardening artifact correction, the gray level difference in the titanium alloy area becomes smaller, which is more conducive to further quantitative analysis of the materials.

[0106] The above embodiments only represent the specific implementation manners of the present invention. The description is relatively specific and detailed, but it should not be construed as a limitation on the scope of the patent of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present invention, several deformations and improvements can still be made, and these all belong to the protection scope of the present invention.

Claims

1. A method for correcting CT image hardening artifacts under a single voltage based on integral invariance, characterized in that: The following steps are involved: Step 1: Acquire an X-ray projection image of a CT scan under a single voltage, divide it by its background value, obtain an X-ray imaging model, and write it into a matrix form; The step 1 is specifically as follows: The X-ray tube emits continuous energy spectrum X-rays, and the imaging model is written as: I=I0∫S(E)exp[-∫ L u(x,E)dx]dE Where: I0 represents the initial intensity of the ray, I is the intensity after ray attenuation, which is equivalent to the background grayscale and image grayscale in actual imaging, L represents the path passed by the X-ray, u(x, E) represents the attenuation coefficient of the ray with energy E at point x on the path L, and S(E) represents the normalized equivalent energy spectrum, which satisfies: 1=∫S(E)dE The imaging model is discretized into a summation, and both sides of the equation are divided by I0: Where: E r is the energy of the rth narrow energy spectrum, S(E r ) is the weighting coefficient of the rth narrow energy spectrum image, K is the type of material in the object, u rk For the kth material in E r The attenuation coefficient, d k is the length of the kth material on the path of X-rays; where S(E r )satisfy: Assume that the number of X-ray projection pixels at the jth scanning angle is M j (j=1, 2, ..., J), the sum of the number of projected image pixels under all scanning angles is recorded as M, then: The I / I0 corresponding to the mth pixel of the X-ray projection image is denoted as f m , m=1,2,…,M; then the X-ray imaging model is written in matrix form as: F=Sexp(-UD) Where: F=(f1,f2,……,f M ), S is the narrow energy spectrum image weighting coefficient that constitutes each voltage X-ray image, that is, s r =S(E r ), a column of U corresponds to the attenuation coefficient of a material in each narrow spectral band, and the corresponding row of D is the attenuation length of the same material at each pixel; Step 2: Obtain an X-ray projection image decomposition model according to the matrix form of the X-ray imaging model; Step 3: Adding an integral invariance equation to the projection image decomposition model to constrain the narrow energy spectrum projection obtained by projection decomposition under a single voltage, thereby achieving hardening artifact correction of CT images under a single voltage; The integral invariance equation in step 3 is as follows:

2. The method for correcting CT image hardening artifacts under a single voltage based on integral invariance according to claim 1, characterized in that: The step 2 is specifically as follows: The volume expression of the kth material at the jth angle is: This value has nothing to do with the angle and is determined by the material distribution of the object itself, so: T kj =T k(j+1) ,j=1,2,…,J-1,k=1,2,…,K Get the X-ray projection image decomposition model:

3. The method for correcting CT image hardening artifacts under a single voltage based on integral invariance according to claim 1, characterized in that: The step 3 is specifically as follows: Set the energy spectrum value R(R) and record the maximum ray energy of the image as E max , with E max / R is the energy spectrum interval to divide the continuous energy spectrum to obtain the narrow energy spectrum interval [E' r-1 ,E' r ], r=1,2,……,R, and E'0=0, E' R =E max ; The middle energy of the rth narrow energy spectrum is recorded as E r , that is, E r =(E' r-1 +E' r ) / 2; The material type of the imaging object is recorded as K. The attenuation coefficient of each material at each energy can be obtained by querying the public material attenuation coefficient database. r The attenuation coefficient is denoted as u rk , remember U=(u rk ) RK ; Initialize S, and assign S to s1=s2=…=s R =1 / R; Initialize D, D=(d km ) KM , where d km is the length of the kth material on the path of the X-ray corresponding to the mth projection pixel point, and the maximum distance between any two points on the imaging object is recorded as d max , generate a random matrix with the same size as D as the initial value of D, and each element value ranges from 0 to d max between; Set the maximum number of iterations Niter, the reference relative error ε for stopping the iteration, and the number of iterations is recorded as niter; calculate the objective function value under the current variable Denoted as h0, initialize niter = 0; Let niter = niter + 1, according to the formula Iteratively update S and normalize the rows of S; according to the formula Iteratively update D; represents the Hadamard product of the matrix, that is, the corresponding elements are multiplied, and - represents the Hadamard division of the matrix, that is, the corresponding elements are divided; According to the integral invariance equation Update all d km ; Calculate the objective function value under the current variable Denoted as h niter ; Calculate (h niter -h (niter-1) ) / h niter and if its absolute value is less than ε, end the iteration and go to the next step. Otherwise, check if niter is less than Niter. If niter < Niter, continue the iterative update. If niter ≥ Niter, go to the next step; The current UD is calculated, where each line is a narrow energy spectrum projection. The narrow energy spectrum projection is reconstructed to obtain a CT image corrected for hardening artifacts.

Images