A hybrid iterative CT imaging method and apparatus for metal artifact correction

Abstract

The application discloses a hybrid iterative CT imaging method and device for metal artifact correction, which comprises the following steps: step 1, obtaining initial values of a density image f of a first material, a density image g of a second material and a density image h of a third material; step 2, converting g and h into binary images with pixel values of 0 and value i initial values according to the initial values of value i , obtaining a segmented pixel value and a critical value set R i and τ i ; step 3, obtaining boundary pixel points of g n and h n after n iterations of g and h, forming a union set U, updating pixel points of f n after n iterations and pixel points of the union set U, keeping pixel values of g n and h n unchanged outside the union set U, performing bilateral filtering on f n , g n and h n , and obtaining f n+1 after n+1 iterations; step 4, updating R i and τ i when a preset iteration number is reached, and finally obtaining g n+1 and h n+1 after n+1 iterations; and step 5, judging whether a termination condition is reached, and returning to step 3 if not. The application can effectively remove metal artifacts in a CT reconstructed image and retain the tissue structure of an original image.

Description

Technical Field

[0001] This invention relates to the field of CT (X-ray Computed Tomography) imaging technology, and in particular to a hybrid iterative CT imaging method and apparatus for correcting metal artifacts. Background Technology

[0002] CT scans, with their advantages of high spatial resolution, short scan time, large scanning range, and ease of operation, are widely used in medical diagnosis. However, in clinical practice, it has been found that when performing CT scans on patients with implanted metallic prostheses, metal artifacts appear in the reconstructed images because the attenuation coefficient of metal is much higher than that of normal tissue in the human body. Metal artifacts severely degrade image quality and significantly affect the doctor's diagnosis, especially since some artifacts may closely resemble certain lesions, potentially leading to misdiagnosis. Therefore, researching methods for metal artifact reduction (MAR) to improve CT imaging quality is of great significance to the field of radiological diagnosis.

[0003] Based on the data format and processing method, metal artifact removal methods can generally be categorized into three types:

[0004] The first type is MAR (Mental Artifact Reduction), which is based on projection domain data inpainting. This method uses projection data from non-metallic trajectory regions—normal data unaffected by metal—to interpolate and replace the metallic trajectory regions to remove artifacts. The metallic trajectory is obtained by projecting the metal segments from the original image. However, improper methods of inpainting the projection domain data can introduce new artifacts, and the accuracy of metal segmentation and the masking of metallic trajectory regions affects the artifact removal effect. For severe and complex artifacts, the correction effect is not ideal.

[0005] The second type of method is based on iterative MAR (Mathematical Arithmetic Reconstruction), which combines a physical model with image domain data through iterative calculations to correct the original projection data and reconstruct the restored image, continuously improving the restoration effect. Currently, another team has proposed an iterative method based on prior images to remove metal artifacts. This method combines the advantages of statistical methods and methods with complete projection data to estimate and correct data bias caused by metal implants. However, the objective function of this algorithm has two regularization terms, requiring significant time for parameter tuning, and it does not consider multicolor X-ray energy spectrum information, resulting in suboptimal image quality.

[0006] The third type of method is deep learning-based MAR, which utilizes convolutional neural networks. It requires training with a large amount of data to extract fine image features and fit a complex mapping relationship from artifact-laden images to artifact-free images, thereby removing artifacts. Summary of the Invention

[0007] The purpose of this invention is to provide a hybrid iterative CT imaging method and apparatus for correcting metal artifacts to overcome or at least mitigate at least one of the aforementioned defects of the prior art.

[0008] To achieve the above objectives, the present invention provides a hybrid iterative CT imaging method for metal artifact correction, comprising:

[0009] Step 1: The original image contains three materials: material 1, material 2, and material 3. Based on the orthographic projection data of the original image, obtain the density image f = (f1, f2, ... f3) of material 1. j ,…f J ) T The density graph of the second material is g = (g1, g2, ... g). j ,…g J ) T The density graph of the third material is h = (h1, h2…h) j ,…h J ) T The initial value; where f j g j h j Let f, g, and h be the pixel values ​​of f, g, and h at the j-th pixel, respectively, where J is the total number of pixels;

[0010] Step 2, based on the respective values ​​of the second and third materials. i The initial value is the average density of the i-th material. The density images g of the second material and h of the third material are both converted to pixel values ​​of 0 and value. i The initial binary image is used to obtain a series of segmented pixel values ​​r of the i-th material. ik This forms the set of segmented pixel values ​​R for the i-th material. i =[r i1 ,r i2 ,…,r ik …,r il ], r ik The specific value selection is in value i The initial value's neighborhood is given, where l is the total number of segmented pixel values, thus obtaining a series of segmentation pixel threshold values ​​τ for the i-th material. ik This forms the set of critical values ​​τ for segmented pixels of the i-th material. i =[τ i1 ,τi2 ,…τ i(l-1) ], k is the index of the segmented pixel value and the segmented pixel threshold value, i is the index of the material type, i=1 corresponds to the first material, i=2 corresponds to the second material, i=3 corresponds to the third material;

[0011] Step 3: First, obtain the density image g of g and h after the nth iteration. n and h n The boundary pixels, and take g n and h n The union of the boundary pixels U is then used to reconstruct the density image f after the nth iteration using the E-algebraic reconstruction algorithm. n Update the pixel values ​​of all pixels and pixels within the union U, while maintaining the values ​​of g outside the union U. n pixel value g j and h n pixel value h j Unchanged, finally for f n g n and h n Bilateral filtering is performed to obtain the density image f of the first material after the (n+1)th iteration. n+1 ;

[0012] Step 4: After reaching the preset number of iterations, take the g values ​​after bilateral filtering. n and h n The average of all non-zero pixel values ​​is used as the value. i The current value of g after bilateral filtering n and h n Convert to pixel values ​​of 0 and value i The current value's binary image, then update R. i and τ i Finally, the discrete density image g of the second material after the (n+1)th iteration is obtained using the following equation (1). n+1 and the discrete density image h of the third material n+1 ;

[0013]

[0014] In the formula, T(v i ) corresponds to g n+1 g j and h n+1 h j One of them, v i For the corresponding g n g j and h n h j One of them, r i1 Let r be the value of the first segmented pixel of the i-th material.i2 r is the value of the second segmented pixel of the i-th material. il For the i-th material, the l-th segmentation pixel value, τ i1 This represents the threshold value for the first segmentation pixel of the i-th material. τ i2 This represents the threshold value for the second segmentation pixel of the i-th material. The threshold value for the (l-1)th segmentation pixel of the i-th material.

[0015] Step 5: Determine if the termination condition has been met. If not, return to step 3.

[0016] Furthermore, the first material is soft tissue, the second material is bone, and the third material is metal.

[0017] Furthermore, the initial values ​​of the density image f of the first material, the density image g of the second material, and the density image h of the third material in step 1 are obtained using a filtered back projection algorithm.

[0018] The present invention also provides a hybrid iterative CT imaging device for metal artifact correction, comprising:

[0019] The density image initialization unit is used to obtain the density image f = (f1, f2, ... f3) of the first material based on the first material, the second material, and the third material contained in the original image, and based on the orthographic projection data of the original image. j ,…f J ) T The density graph of the second material is g = (g1, g2, ... g). j ,…g J ) T The density graph of the third material is h = (h1, h2…h) j ,…h J ) T The initial value; where f j g j h j Let f, g, and h be the pixel values ​​of f, g, and h at the j-th pixel, respectively, where J is the total number of pixels;

[0020] Segmenting pixel units, which are used to divide pixels according to the values ​​of the second and third materials respectively. i The initial value is the average density of the i-th material. The density images g of the second material and h of the third material are both converted to pixel values ​​of 0 and value. i The initial binary image is used to obtain a series of segmented pixel values ​​r of the i-th material. ik This forms the set of segmented pixel values ​​R for the i-th material. i =[r i1 ,ri2 ,…,r ik …,r il ], r ik The specific value selection is in value i The initial value's neighborhood is given, where l is the total number of segmented pixel values, thus obtaining a series of segmentation pixel threshold values ​​τ for the i-th material. ik This forms the set of critical values ​​τ for segmented pixels of the i-th material. i =[τ i1 ,τ i2 ,…τ i(l-1) ], k is the index of the segmented pixel value and the segmented pixel threshold value, i is the index of the material type, i=1 corresponds to the first material, i=2 corresponds to the second material, i=3 corresponds to the third material;

[0021] The first iteration unit is used to first acquire the density image g of g and h after the nth iteration. n and h n The boundary pixels, and take g n and h n The union of the boundary pixels U is then used to reconstruct the density image f after the nth iteration using the E-algebraic reconstruction algorithm. n Update the pixel values ​​of all pixels and pixels within the union U, while maintaining the values ​​of g outside the union U. n pixel value g j and h n pixel value h j Unchanged, finally for f n g n and h n Bilateral filtering is performed to obtain the density image f of the first material after the (n+1)th iteration. n+1 ;

[0022] The second iteration unit is used to take the bilaterally filtered g after reaching a preset number of iterations. n and h n The average of all non-zero pixel values ​​is used as the value. i The current value of g after bilateral filtering n and h n Convert to pixel values ​​of 0 and value i The current value's binary image, then update R. i and τ i Finally, the discrete density image g of the second material after the (n+1)th iteration is obtained using the following equation (1). n+1 and the discrete density image h of the third material n+1 :

[0023]

[0024] In the formula, T(v i ) corresponds to g n+1 g j and h n+1 h j One of them, v i For the corresponding g n g j and h n h j One of them, r i1 Let r be the value of the first segmented pixel of the i-th material. i2 r is the value of the second segmented pixel of the i-th material. il For the i-th material, the l-th segmentation pixel value, τ i1 This represents the threshold value for the first segmentation pixel of the i-th material. τ i2 This represents the threshold value for the second segmentation pixel of the i-th material. τ i(l-1) The threshold value for the (l-1)th segmentation pixel of the i-th material.

[0025] The iteration termination judgment unit is used to determine whether the termination condition has been met. If not, the first iteration unit and the second iteration unit execute the iteration operation in sequence.

[0026] Furthermore, the first material is soft tissue, the second material is bone, and the third material is metal.

[0027] Furthermore, the initial values ​​of the density image f of the first material, the density image g of the second material, and the density image h of the third material in the density image initialization unit are obtained using a filtered back projection algorithm.

[0028] The present invention has the following advantages due to the adoption of the above technical solutions:

[0029] This invention takes into account multicolor X-ray energy spectrum information and energy spectrum CT reconstruction algorithms, and proposes a single-energy spectrum hybrid iterative image reconstruction algorithm by combining the idea of ​​Discrete Algebraic Reconstruction Technique (DART). The entire iterative process does not require supervision, effectively removes metal artifacts in CT reconstructed images, and preserves the original image tissue structure. Attached Figure Description

[0030] Figure 1 The figures show the phantoms used in the numerical experiments. In the figure, (a) is the numerical phantom, (b) is the bone material, (c) is the water material, and (d) is the titanium material.

[0031] Figure 2The diagram shows the image results obtained by the method of this invention under noise-free conditions using FBP (uncorrected image), linear interpolation (LIMAR), beam hardening correction (BCMAR), normalized linear interpolation (NMAR), etc. The display window (center and width) of the first four images is set to (200, 600) HU, and the display window of the last image is set to (50, 300) HU.

[0032] Figure 3 The diagram shows the image results obtained by the method of this invention under noisy conditions using FBP (uncorrected image), linear interpolation (LIMAR), beam hardening correction (BCMAR), and normalized linear interpolation (NMAR). Detailed Implementation

[0033] The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0034] The hybrid iterative CT imaging method for metal artifact correction provided in this invention includes:

[0035] Step 1: The original image contains three materials: a first material, a second material, and a third material. Based on the orthographic projection data p of the original image... a Obtain the density image of the first material f = (f1, f2, ... f) j , ...f J ) T The density graph of the second material is g = (g1, g2, ... g2). j , ...g J ) T The density graph of the third material is h = (h1, h2…h) j , ...h J ) T The initial value; where f j g j h j Let f, g, and h be the pixel values ​​of f, g, and h at the j-th pixel, respectively, where J is the total number of pixels.

[0036] Step 2: Calculate the value of the second and third materials respectively. i initial value, value i The initial value is the average density of the i-th material, for example, the value of titanium. i 4.5 g / cm 3 value of bones i 1.8 g / cm 3 Soft tissue value i 1g / cm 3 The value of metallic iron i 7.9 g / cm3 .

[0037] Then, based on the respective values ​​of the second and third materials... i The initial value is the average density of the i-th material. The density images g of the second material and h of the third material are both converted to pixel values ​​of 0 and value. i The initial binary image can be used to obtain a series of segmented pixel values ​​r of the i-th material. ik This forms the set of segmented pixel values ​​R for the i-th material. i =[r i1 r i2 , ..., r ik …, r il ], r ik The specific value selection is in value i The initial value's neighborhood, k is the index of the segmented pixel value, l is the total number of segmented pixel values, and i is the material type index, i=1 corresponds to the first material, i=2 corresponds to the second material, i=3 corresponds to the third material, for example: l=5, r i1 =value i -0.2, r i2 =value i -0.1, r i3 =value i r i4 =value i +0.1, r i5 =value i +0.2.

[0038] According to r ik The specific value, using A series of segmentation pixel threshold values ​​τ of the i-th material can be obtained. ik This forms the set of critical values ​​τ for segmented pixels of the i-th material. i =[τ i1 , τ i2 ,…τ i(l-1) ], k is also the index of the segmentation pixel threshold.

[0039] Step 3: First, obtain the density image g of g and h after the nth iteration. n and h n The boundary pixels, and take g n and h n The union of the boundary pixels U. Then, the extended algebraic reconstruction technique (E-ART) is used to reconstruct the density image f after the nth iteration.n Update the pixel values ​​of all pixels and pixels within the union U, while maintaining the values ​​of g outside the union U. n pixel value g j and h n pixel value h j Keeping it unchanged, finally use a bilateral filter to adjust f. n g n and h n Bilateral filtering is performed to obtain the density image f of the first material after the (n+1)th iteration. n+1 The bilateral filter, used to smooth images and reduce noise while preserving edges, is a non-linear bilateral filter. The kernel function of the bilateral filter is a combination of the spatial domain kernel and the pixel range domain kernel. A bilateral filter with a radius of 1, a variance of 1 for the spatial domain kernel, and a variance of 0.1 for the range kernel is applied to the image.

[0040] Step 4: During the entire iterative process (from Step 3 to Step 5), after reaching the preset number of iterations (e.g., 3 times), the bilaterally filtered g... n and h n Convert to a binary image and obtain the value. i Update R to its current value. i and τ i Finally, the discrete density image g of the second material after the (n+1)th iteration is obtained using the following equation (1). n+1 and the discrete density image h of the third material n+1 ;

[0041]

[0042] In the formula, T(v i ) corresponds to g j and h j One of them, v i Let r be the current pixel value of the i-th material. i1 Let r be the value of the first segmented pixel of the i-th material. i2 r is the value of the second segmented pixel of the i-th material. il For the i-th material, the l-th segmentation pixel value, τ i1 This represents the threshold value for the first segmentation pixel of the i-th material. τ i2 This represents the threshold value for the second segmentation pixel of the i-th material. τ i(l-1) The threshold value for the (l-1)th segmentation pixel of the i-th material.

[0043] Step 5: Determine if the termination condition has been met. If not, return to step 3. The termination condition can be either executing a fixed number of iterations or the total projection error.

[0044] In one embodiment, the first material is soft tissue, the second material is bone, and the third material is metal.

[0045] In one embodiment, the initial values ​​of the density images f of the first material, g of the second material, and h of the third material in step 1 are obtained using a filtered back projection (FBP) algorithm. For example, if the first material is soft tissue, the second material is bone, and the third material is metal, then the filtered back projection algorithm specifically includes:

[0046] Before iteration, an initial density image is obtained using a filtered back-projection algorithm. Then, threshold segmentation is used to obtain the metal density image h from the original image, and the metal portion is removed from the original image. Next, the image with the metal portion removed is divided into three categories—air, soft tissue, and bone—using the Otsu's algorithm with multiple thresholds, resulting in the soft tissue density image f and the bone density image g, respectively. At this point, the grayscale values ​​of the three images correspond to the attenuation coefficients. However, this embodiment aims to use the density of the substances for iterative reconstruction. Therefore, dividing each image by the average mass attenuation coefficient of each substance yields the initial values ​​for the soft tissue density image f, the bone density image g, and the metal density image h.

[0047] In addition to this method, the initial values ​​of the density images f of the first material, g of the second material, and h of the third material can also be obtained by using the algebraic reconstruction algorithm (ART), the joint iterative reconstruction algorithm (SIRT), or the joint algebraic reconstruction algorithm (SART).

[0048] In one embodiment, the nonlinear model is linearized by performing a first-order Taylor expansion of the projection data with respect to the current reconstructed image, and finally solved iteratively using the ART algorithm. Considering that the specific applications of the embodiment often involve three base materials: bone, soft tissue, and metal, in order to decompose these three materials, the extended algebraic iterative reconstruction algorithm used in step 3 is characterized by the following equation (2):

[0049]

[0050] In the formula, p aLet I0 and p be the projection data of the a-th X-ray. a I represents the light intensity before and after the X-rays pass through the sample, respectively. a It is characterized by the following equation (3) or equation (4), where S(E) is the energy spectrum information of X-rays at incident energy E, S(E) = I0(E) / I0, I0(E) is the initial intensity of X-rays at energy level E, and μ(x,E) is the attenuation coefficient of the material at pixel x in the image at incident energy E, characterized by the following equation (5), l a Let be the path of the a-th X-ray, where a = 1, 2, ..., L, and L represents the total number of X-rays.

[0051]

[0052]

[0053] In the formula, μ(x) is the linear attenuation coefficient of the material at pixel x in the image.

[0054]

[0055] Where f(x), g(x), and h(x) are respectively equivalent to the density image f of the first material, the density image g of the second material, and the density image h of the third material in step 1. θ(E) and ω(E) are the mass attenuation coefficients for the first, second, and third materials, respectively. The specific values ​​of these three coefficients can be found in the X-ray mass attenuation coefficient table published by the National Institute of Standards and Technology (NIST).

[0056] Therefore, the base material information is incorporated, and the energy sampling interval is set to 1. The interval can also be represented by δ. For simplicity, here we directly set the effective energy range of the X-ray spectrum to p. a Subintervals, projected data p a The possible discrete expression is as follows (6):

[0057]

[0058] Among them, S m , θ m ω m They are S(E), θ(E) and ω(E) are the sampled values ​​of the m-th subinterval, and S(E) is the sampled value of ... The discrete forms of θ(E) and ω(E) are respectively (S1, S2, ... S M ) T , (θ1, θ2, ... θ) M ) T (ω1, ω2, ... ω)M ) T D a =(d a1 d a2 , ...d aJ ) T Let d be the projection vector. aj This represents the contribution of the j-th pixel value to the projection along the a-th X-ray path, and is generally equal to the length of the intersection line between the a-th X-ray path and the pixel corresponding to the j-th pixel value.

[0059] Then, the density images of the three materials obtained in the nth iteration, i.e., the current iteration image (f) n g n h n A first-order Taylor expansion of the nonlinear model at point ) transforms it into the following linear equation (7):

[0060]

[0061] Wherein, i.e., the current iteration image (f) n g n h n The projection data obtained after the nth iteration by forward projection along the a-th X-ray path.

[0062] In the projection data p a and the current iteration image (f) n g n h n Given the above equations, the orthogonal projection points on the hyperplane corresponding to the above equations are calculated using an algebraic reconstruction algorithm. The specific iterative form is as follows: (8)

[0063]

[0064] This algorithm can reconstruct high-quality images without noise by utilizing energy spectrum information. However, the reconstruction problem of energy spectrum CT is pathological; during reconstruction, noise in the projection data is further amplified, leading to severe degradation of the final reconstructed image.

[0065] Discrete algebraic iterative reconstruction algorithm is a type of discrete iterative algorithm that uses the ideas of discrete CT in its reconstruction process. During iteration, it only considers updating boundary points, and through adaptive adjustment of the boundaries, the reconstructed image boundaries continuously approach the true boundaries. This significantly reduces the number of unknowns and greatly improves the ill-conditioned nature of the problem. Therefore, it can be inferred that the discrete algebraic iterative reconstruction algorithm has positive implications for energy spectrum CT image reconstruction.

[0066] The algorithm of this invention is based on the idea of ​​discrete algebraic iterative reconstruction algorithm. It utilizes the characteristic of shared boundaries between different base material images. When updating the base material image, the updated point is the union of all image boundary points. Through iterative feedback and correction of the projection, it can effectively reduce metal artifacts and suppress noise. Because it considers energy spectrum information, the reconstructed image is better than conventional methods, further improving image resolution and material decomposition accuracy.

[0067] This invention also provides a hybrid iterative CT imaging device for metal artifact correction, comprising a density image initialization unit, a pixel segmentation unit, a first iteration unit, a second iteration unit, and an iteration termination determination unit, wherein:

[0068] The density image initialization unit is used to obtain the density image f = (f1, f2, ... f3) of the first material based on the first material, the second material, and the third material contained in the original image, and based on the orthographic projection data of the original image. j , ...f J ) T The density graph of the second material is g = (g1, g2, ... g2). j , ...g J ) T The density graph of the third material is h = (h1, h2…h) j , ...h J ) T The initial value; where f j g j h j Let f, g, and h be the pixel values ​​of f, g, and h at the j-th pixel, respectively, where J is the total number of pixels.

[0069] The pixel segmentation unit is used to divide the pixel based on the values ​​of the second and third materials respectively. i The initial value is the average density of the i-th material. The density images g of the second material and h of the third material are both converted to pixel values ​​of 0 and value. i The initial binary image is used to obtain a series of segmented pixel values ​​r of the i-th material. ik This forms the set of segmented pixel values ​​R for the i-th material. i =[r i1 r i2 , ..., r ik …, r il ], r ik The specific value selection is in value i The initial value's neighborhood is given, where l is the total number of segmented pixel values, thus obtaining a series of segmentation pixel threshold values ​​τ for the i-th material. ik This forms the set of critical values ​​τ for segmented pixels of the i-th material. i =[τi1 , τ i2 ,…τ i(l-1 [], k is the index of the segmented pixel value and the segmented pixel threshold value, i is the index of the material type, i=1 corresponds to the first material, i=2 corresponds to the second material, and i=3 corresponds to the third material.

[0070] The first iteration unit is used to first obtain the density image g of g and h after the nth iteration. n and h n The boundary pixels, and take g n and h n The union of the boundary pixels U is then used to reconstruct the density image f after the nth iteration using the E-algebraic reconstruction algorithm. n Update the pixel values ​​of all pixels and pixels within the union U, while maintaining the values ​​of g outside the union U. n pixel value g j and h n pixel value h j Unchanged, finally for f n g n and h n Bilateral filtering is performed to obtain the density image f of the first material after the (n+1)th iteration. n+1 .

[0071] The second iteration unit is used to take the bilaterally filtered g after reaching a preset number of iterations. n and h n The average of all non-zero pixel values ​​is used as the value. i The current value of g after bilateral filtering n and h n Convert to pixel values ​​of 0 and value i The current value's binary image, then update R. i and τ i Finally, the discrete density image g of the second material after the (n+1)th iteration is obtained using the above equation (1). n+1 and the discrete density image h of the third material n+1 .

[0072] The iteration termination judgment unit is used to determine whether the termination condition has been met. If not, the first iteration unit and the second iteration unit execute the iteration operation in sequence.

[0073] In one embodiment, the initial values ​​of the density image f of the first material, the density image g of the second material, and the density image h of the third material in the density image initialization unit are obtained using a filtered back projection algorithm.

[0074] This invention compares and simulates traditional linear interpolation methods (LIMAR), normalized interpolation methods (NMAR), and beam hardening correction MAR algorithms (BCMAR) using energy spectrum information. Specific experimental data are as follows:

[0075] The phantom used in the numerical experiments of this invention is a titanium-added Forbild head phantom, such as... Figure 1 As shown. The sample is composed of three substances: bone, water, and titanium. The density of bone is 1.8, water was selected with three densities: 1.0450, 1.0500, and 1.0600, and titanium has a density of 4.5. In the simulation, a parallel beam setup was used to acquire 360 ​​equally spaced projections at 180-degree intervals, with a field of view diameter of 303 mm. 10 tests were conducted. 6 The noise-free and Poisson noise corresponding to the photon emission flux. The reconstructed image size is 620*620.

[0076] Figure 2 This image compares the proposed method with linear interpolation (LIMAR), normalized linear interpolation (NMAR), beam hardening correction (BCMAR) under noise-free conditions. From left to right, the images are: uncorrected image (FBP result), LIMAR result, NMAR result, BCMAR result, and the result of the proposed method. The display window (center and width) for the first four images is set to (200, 600) HU, and the display window for the last image is set to (50, 300) HU.

[0077] Figure 3 This paper compares the proposed method with linear interpolation (LIMAR), normalized linear interpolation (NMAR), beam hardening correction (BCMAR) under noisy conditions. From left to right, the images are the uncorrected image (FBP result), LIMAR result, NMAR result, BCMAR result, and the result of the proposed method. The display window (center and width) of the first four images is set to (200, 600) HU, and the display window of the last image is set to (50, 300) HU.

[0078]

[0079] Table 1 Comparison of SSIM and PSNR

[0080] Reconstruction results as follows Figure 2 and Figure 3As shown in the figure. From the contour and reconstruction results, this invention found that the proposed algorithm has better spatial resolution and noise suppression than other methods. For the method mentioned in this invention, the structural similarity (SSIM) and peak signal-to-noise ratio (PSNR) of the de-metal artifact image and the standard image were calculated respectively, and the results are shown in Table 1. The data in the table show that the structural similarity and PSNR results of the proposed method are significantly higher than other methods.

[0081] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them. Those skilled in the art should understand that modifications can be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A hybrid iterative CT imaging method for metal artifact correction, characterized in that, include: Step 1: The original image contains three materials: material 1, material 2, and material 3. Based on the orthographic projection data of the original image, obtain the density image f = (f1, f2, ... f3) of material 1. j ,…f J ) T The density graph of the second material is g = (g1, g2, ... g). j ,…g J ) T The density graph of the third material is h = (h1, h2…h) j ,…h J ) T The initial value; where f j g j h j Let f, g, and h be the pixel values ​​of f, g, and h at the j-th pixel, respectively, where J is the total number of pixels; Step 2, based on the respective values ​​of the second and third materials. i The initial value is the average density of the i-th material. The density images g of the second material and h of the third material are both converted to pixel values ​​of 0 and value. i The initial binary image is used to obtain a series of segmented pixel values ​​r of the i-th material. ik This forms the set of segmented pixel values ​​R for the i-th material. i =[r i1 ,r i2 ,…,r ik …,r il ], r ik The specific value selection is in value i The initial value's neighborhood is given, where l is the total number of segmented pixel values, thus obtaining a series of segmentation pixel threshold values ​​τ for the i-th material. ik This forms the set of critical values ​​τ for segmented pixels of the i-th material. i =[τ i1 ,τ i2 ,… i(l-1) ], k is the index of the segmented pixel value and the segmented pixel threshold value, i is the index of the material type, i=1 corresponds to the first material, i=2 corresponds to the second material, i=3 corresponds to the third material; Step 3: First, obtain the density image g of g and h after the nth iteration. n and h n The boundary pixels, and take g n and h n The union of the boundary pixels U is then used to reconstruct the density image f after the nth iteration using the E-algebraic reconstruction algorithm. n Update the pixel values ​​of all pixels and pixels within the union U, while maintaining the values ​​of g outside the union U. n pixel value g j and h n pixel value h j Unchanged, finally for f n g n and h n Bilateral filtering is performed to obtain the density image f of the first material after the (n+1)th iteration. n+1 ; Step 4: After reaching the preset number of iterations, take the g values ​​after bilateral filtering. n and h n The average of all non-zero pixel values ​​is used as the value. i The current value of g after bilateral filtering n and h n Convert to pixel values ​​of 0 and value i The current value's binary image, then update R. i and τ i Finally, the discrete density image g of the second material after the (n+1)th iteration is obtained using the following equation (1). n+1 and the discrete density image h of the third material n+1 ; In the formula, T(v i ) corresponds to g n+1 g j and h n+1 h j One of them, v i For the corresponding g n g j and h n h j One of them, r i1 Let r be the value of the first segmented pixel of the i-th material. i2 r is the value of the second segmented pixel of the i-th material. il For the i-th material, the l-th segmentation pixel value, τ i1 This represents the threshold value for the first segmentation pixel of the i-th material. τ i2 This represents the threshold value for the second segmentation pixel of the i-th material. τ i(l-1) The threshold value for the (l-1)th segmentation pixel of the i-th material. Step 5: Determine if the termination condition has been met. If not, return to step 3.

2. The hybrid iterative CT imaging method for metal artifact correction as described in claim 1, characterized in that, The first material is soft tissue, the second material is bone, and the third material is metal.

3. The hybrid iterative CT imaging method for metal artifact correction as described in claim 1 or 2, characterized in that, The initial values ​​of the density images f of the first material, g of the second material, and h of the third material in step 1 are obtained using a filtered back projection algorithm.

4. A hybrid iterative CT imaging device for metal artifact correction, characterized in that, include: The density image initialization unit is used to obtain the density image f = (f1, f2, ... f3) of the first material based on the first material, the second material, and the third material contained in the original image, and based on the orthographic projection data of the original image. j ,…f J ) T The density graph of the second material is g = (g1, g2, ... g). j ,…g J ) T The density graph of the third material is h = (h1, h2…h) j ,…h J ) T The initial value; where f j g j h j Let f, g, and h be the pixel values ​​of f, g, and h at the j-th pixel, respectively, where J is the total number of pixels; Segmenting pixel units, which are used to divide pixels according to the values ​​of the second and third materials respectively. i The initial value is the average density of the i-th material. The density images g of the second material and h of the third material are both converted to pixel values ​​of 0 and value. i The initial binary image is used to obtain a series of segmented pixel values ​​r of the i-th material. ik This forms the set of segmented pixel values ​​R for the i-th material. i =[r i1 ,r i2 ,…,r ik …,r il ], r ik The specific value selection is in value i The initial value's neighborhood is given, where l is the total number of segmented pixel values, thus obtaining a series of segmentation pixel threshold values ​​τ for the i-th material. ik This forms the set of critical values ​​τ for segmented pixels of the i-th material. i =[τ i1 ,τ i2 ,…τ i(l-1) ], k is the index of the segmented pixel value and the segmented pixel threshold value, i is the index of the material type, i=1 corresponds to the first material, i=2 corresponds to the second material, i=3 corresponds to the third material; The first iteration unit is used to first acquire the density image g of g and h after the nth iteration. n and h n The boundary pixels, and take g n and h n The union of the boundary pixels U is then used to reconstruct the density image f after the nth iteration using the E-algebraic reconstruction algorithm. n Update the pixel values ​​of all pixels and pixels within the union U, while maintaining the values ​​of g outside the union U. n pixel value g j and h n pixel value h j Unchanged, finally for f n g n and h n Bilateral filtering is performed to obtain the density image f of the first material after the (n+1)th iteration. n+1 ; The second iteration unit is used to take the bilaterally filtered g after reaching a preset number of iterations. n and h n The average of all non-zero pixel values ​​is used as the value. i The current value of g after bilateral filtering n and h n Convert to pixel values ​​of 0 and value i The current value's binary image, then update R. i and τ i Finally, the discrete density image g of the second material after the (n+1)th iteration is obtained using the following equation (1). n+1 and the discrete density image h of the third material n+1 ; In the formula, T(v i ) corresponds to g n+1 g j and h n+1 h j One of them, v i For the corresponding g n g j and h n h j One of them, r i1 Let r be the value of the first segmented pixel of the i-th material. i2 r is the value of the second segmented pixel of the i-th material. il For the i-th material, the l-th segmentation pixel value, τ i1 This represents the threshold value for the first segmentation pixel of the i-th material. τ i2 This represents the threshold value for the second segmentation pixel of the i-th material. τ i(l-1) The threshold value for the (l-1)th segmentation pixel of the i-th material. The iteration termination judgment unit is used to determine whether the termination condition has been met. If not, the first iteration unit and the second iteration unit execute the iteration operation in sequence.

5. The hybrid iterative CT imaging apparatus for metal artifact correction as described in claim 4, characterized in that, The first material is soft tissue, the second material is bone, and the third material is metal.

6. The hybrid iterative CT imaging apparatus for metal artifact correction as described in claim 4 or 5, characterized in that, The initial values ​​of the density image f of the first material, the density image g of the second material, and the density image h of the third material in the density image initialization unit are obtained using a filtered back projection algorithm.

Images