Adaptive Structure-Preserving Correction Method, Device and Imaging Equipment for CT Image Ring Artifacts

Abstract

The present invention discloses a method, apparatus and imaging device for adaptively correcting ring artifacts in CT images while preserving the structure. The steps included in the invention are as follows: Step 1, obtaining the original CT projection image of the object to be measured; Step 2, performing multi-scale decomposition on the original CT projection image to obtain initial multi-scale images with different sizes; Step 3, performing boundary-preserving filtering on each initial multi-scale image and outputting the first multi-scale image; Step 4, removing the frequency information corresponding to the artifact information in the first multi-scale image through an artifact information correction method and outputting the second multi-scale image; Step 5, sequentially reconstructing and rebuilding the second multi-scale image to obtain a CT image with ring artifacts removed. The present invention can effectively correct ring artifacts in CT images only through software algorithms without the need for special hardware support.

Description

Technical Field

[0001] The present invention relates to the technical field of CT (English full name: "Computed Tomography", Chinese full name: "Computed Tomography") image processing, and particularly to an adaptive structure-preserving correction method, device and imaging equipment for CT image ring artifacts. Background Art

[0002] An X-ray CT imaging device generally consists of a ray source subsystem, a detector subsystem, a mechanical scanning subsystem, a signal control subsystem, a shielding facility, etc. The basic process of X-ray CT imaging is to place the object to be measured on the sample stage, turn on the ray source to generate an X-ray beam, and at the same time turn on the detector to detect the X-rays after interacting with the sample to be measured; the sample to be measured rotates one week relative to the ray source and the detector to obtain CT scan data; perform certain processing on the scan data, and use a reconstruction algorithm to obtain a CT image reflecting the internal structure of the object to be measured. During the CT imaging process, there are various factors affecting the quality of CT images, such as the ray energy spectrum distribution, the stability of the ray tube current, the material of the filter, the detector response efficiency, and the reconstruction algorithm, etc.

[0003] Currently, classical CT reconstruction algorithms (such as FBP, ART) always assume that the detection efficiency of the detector is 100%, ignoring the influence of multi-energy rays and scattered photons. In fact, commonly used X-ray detectors detect X-rays indirectly, that is, convert X-rays into visible light through a scintillation crystal, then convert the visible light into an electric current through a photosensitive device, and finally integrate the electric current and convert it into a digital signal through an analog-to-digital converter. The scintillation crystal and the photosensitive device are coupled to form a detector unit, and are arranged in a certain geometric array (such as a linear array detector, a planar array detector) to form a detector subsystem. During the X-ray detection process, the detection efficiencies of different detector units for X-rays of the same energy are often inconsistent. Even for the same detector unit, its response efficiency to photons of different energies is non-linear, and the X-rays generated by the commonly used X-ray light sources in medical and industrial CT imaging devices are composed of multi-energy photons, thus making the inconsistency of the detector unit response efficiency stronger. The inconsistency of the detector unit response efficiency exists in the form of vertical stripes in the projection image and in the form of rings in the reconstructed CT image. The vertical stripe artifacts in the projection image are manifested as parallel stripes with different widths; the ring artifacts in the reconstructed CT image are manifested as a series of concentric circles or arcs.

[0004] Existing circular artifact correction methods can be divided into hardware correction methods, model-driven correction methods, and data-driven correction methods. Hardware correction methods usually control the detector to generate random small offsets during the CT scanning process. By rearranging the data, the vertical stripes in the projection image are transformed into random noise, thereby removing the artifact components. However, the random noise will cause a certain degree of decline in the spatial resolution of the image. In addition, controlling the random offset of the detector will not only increase the system complexity but also significantly increase the scanning time. Model-driven correction methods utilize the characteristics of circular artifacts in the projection image or the reconstructed image to establish their mathematical descriptions, and then design corresponding solution algorithms or image processing algorithms to remove or suppress the circular artifacts in the projection image or the reconstructed image. Wavelet-FFT filtering is a representative correction method. This method removes the vertical stripe artifacts in the projection image through wavelet decomposition and Fourier filtering. However, its filtering process cannot guarantee effective discrimination between artifact components and image boundary components, resulting in incomplete correction in some cases and even possible introduction of new artifacts. Data-driven correction methods usually use a labeled dataset composed of artifact images and ideal correction images to train a deep neural network to obtain a neural network model for artifact recognition or correction. The accuracy of this method strongly depends on the training dataset. In practical applications, due to the difficulty in obtaining the training dataset, circular artifacts are usually difficult to completely remove, and generally, it is necessary to combine with model-driven methods to improve the correction effect.

[0005] In summary, the existing correction methods have defects such as incomplete removal of circular artifacts, easy introduction of new artifacts, and difficulty in ensuring the image spatial resolution. Summary of the Invention

[0006] The purpose of the present invention is to provide a CT image circular artifact adaptive structure-preserving correction method, device, and imaging equipment, which can well retain the boundary and detail information of the image without introducing new artifacts.

[0007] To achieve the above purpose, the technical solutions adopted by the present invention are as follows:

[0008] To achieve the above purpose, the present invention provides a CT image circular artifact adaptive structure-preserving correction method, which includes:

[0009] Step 1, obtain the original CT projection image of the object to be measured;

[0010] Step 2, perform multi-scale decomposition on the original CT projection image to obtain initial multi-scale images W ψ (s, t), and this initial multi-scale image W ψ (s, t) includes a low-frequency coefficient image cA ψ (s, t), a horizontal direction detail coefficient image chD ψ(s, t), the vertical direction detail coefficient image cvD ψ (s, t) and the diagonal direction detail coefficient image cdD ψ (s, t);

[0011] Step 3, perform boundary-preserving filtering on the initial multi-scale image W ψ (s, t), and output the first multi-scale image This first multi-scale image contains the low-frequency coefficient image cA σ (s, t), the horizontal direction detail coefficient image chD ψ (s, t), the boundary-preserving filtered vertical direction detail coefficient image and the diagonal direction detail coefficient image cdD ψ (s, t), where contains boundary information and artifact information. Locate the boundary information g(s, t) through a boundary detection operator, and extract the artifact information R(s, t);

[0012] Among them, the boundary-preserving filtering method specifically includes:

[0013] Step 31, use the following formula (3-1) to filter the vertical direction detail coefficient image cvD ψ in the initial multi-scale image W ψ (s, t) to obtain the filtered vertical direction detail coefficient image

[0014]

[0015] In the formula, ψ is the multi-scale transform basis function, s and t are the scale parameter and the translation parameter respectively, represents the spatial domain filtering operator, and H(s, t) is the filter kernel;

[0016] Step 32, obtain the boundary information g(s, t) of the filtered vertical direction detail coefficient image through the boundary detection operator Edge provided by the following formula (3-2):

[0017]

[0018] Step 33, extract the artifact information R(s, t) through the artifact information extraction operator Ring provided by the following formula (3-3):

[0019]

[0020] Step 4, use the artifact information correction method provided by the following formula (4-1) to remove the first multi-scale image The frequency information corresponding to the artifact information R(s, t) in it, and output the second multi-scale image This second multi-scale image Contains the low-frequency coefficient image cA ψ (s, t), the horizontal direction detail coefficient image chD ψ (s, t), the vertical direction detail coefficient image after removing the artifact information R(s, t) And the diagonal direction detail coefficient image cdD ψ (s, t);

[0021]

[0022] Step 5, perform reconstruction and restoration on the second multi-scale image in sequence to obtain the CT image with the circular artifact removed.

[0023] Further, the method for obtaining the original projection image in Step 1 includes:[[]]

[0024] Step 11, use Equation (1-1) to represent the process of obtaining the scanning data;

[0025]

[0026] In the formula, I(u, β) is the detection intensity of the ray with equivalent energy by the detector unit u after placing the sample , is the detection efficiency of the detector unit u at the equivalent energy of after placing the sample, I0(u) is the detection intensity of the ray with equivalent energy E by the detector unit u without placing the sample, δ(u, E) is the detection efficiency of the detector unit u at the equivalent energy E without placing the sample, σ(u) is the scattered photon intensity detected by the detector unit u, E and are the equivalent energies of the ray before and after placing the sample respectively, l is the integral microelement along the straight line where the ray is located, L(u, β) represents the straight line equation from the ray source focus to the detector unit u, β is the sampling angle, y represents the coordinate of the point on the straight line, and μ(y) is the linear attenuation coefficient of the sample;

[0027] Step 12, perform dark field correction and gain correction using Equation (1-2)

[0028]

[0029] Among them, Z0 is the dark field data collected by the detector when the ray source is not exposed.

[0030] Step 13: Perform negative logarithm transformation on the scanned data corrected in Step 12 using Equation (1-3), and patch and rearrange the scanned data after the negative logarithm transformation to obtain the original CT projection data p(u, β):

[0031] p(u, β) = -lnp0(u, β) (1-3).

[0032] Furthermore, the artifacts in the original projection image appear as vertical stripe artifacts or horizontal stripe artifacts, with or without discontinuities, and the stripe gray-scale changes vary.

[0033] The present invention also provides a CT image circular artifact adaptive structure-preserving correction device, which includes:

[0034] An original CT projection image acquisition unit, which is used to acquire the original CT projection image of the object to be measured;

[0035] A multi-scale decomposition unit, which is used to perform multi-scale decomposition on the original CT projection image to obtain initial multi-scale images W ψ (s, t) of different sizes. This initial multi-scale image W ψ (s, t) includes a low-frequency coefficient image cA ψ (s, t), a horizontal direction detail coefficient image chD ψ (s, t), a vertical direction detail coefficient image cvD ψ (s, t), and a diagonal direction detail coefficient image cdD ψ (s, t);

[0036] A boundary-preserving filtering unit, which is used to perform boundary-preserving filtering on each initial multi-scale image W ψ (s, t) and output a first multi-scale image This first multi-scale image includes a low-frequency coefficient image cA ψ (s, t), a horizontal direction detail coefficient image chD ψ (s, t), a boundary-preserving filtered vertical direction detail coefficient image and a diagonal direction detail coefficient image cdD ψ (s, t), where includes boundary information and artifact information. The boundary information g(s, t) is located through a boundary detection operator, and the artifact information R(s, t) is extracted;

[0037] Among them, the boundary-preserving filtering unit specifically includes:

[0038] A filtering subunit, which is used to use the following formula (3-1) to process the vertical direction detail coefficient image cvD in the initial multi-scale image W ψ (s, t) ψFilter (s, t) to obtain the filtered vertical direction detail coefficient image

[0039]

[0040] Where ψ is the multi-scale transform basis function, s and t are the scale parameter and the translation parameter respectively, represents the spatial domain filtering operator, and H(s, t) is the filter kernel;

[0041] The boundary detection sub-unit, which is used to obtain the boundary information g(s, t) of the filtered vertical direction detail coefficient image through the boundary detection operator Edge provided by the following formula (3-2):

[0042]

[0043] The artifact information extraction sub-unit, which is used to extract the artifact information R(s, t) through the artifact information extraction operator Ring provided by the following formula (3-3):

[0044]

[0045] The artifact correction unit, which is used to remove the frequency information corresponding to the artifact information R(s, t) in the first multi-scale image through the artifact information correction method provided by the following formula (4-1), and output the second multi-scale image This second multi-scale image contains the low-frequency coefficient image cA ψ (s, t), the horizontal direction detail coefficient image chD ψ (s, t), the vertical direction detail coefficient image after removing the artifact information R(s, t) and the diagonal direction detail coefficient image cdD ψ (s, t);

[0046]

[0047] The reconstruction unit, which is used to reconstruct and rebuild the second multi-scale image in sequence to obtain the CT image without ring artifacts.

[0048] Furthermore, the original CT projection image acquisition unit specifically includes:

[0049] The first data correction sub-unit, which represents the scanning data acquisition process using formula (1-1);

[0050]

[0051] ​​where I(u, β) is the detection intensity of the ray with equivalent energy by detector unit u after placing the sample, is the detection efficiency of detector unit u at the equivalent energy after placing the sample, I0(u) is the detection intensity of the ray with equivalent energy E by detector unit u without placing the sample, δ(u, E) is the detection efficiency of detector unit u at the equivalent energy E without placing the sample, σ(u) is the intensity of scattered photons detected by detector unit u, E and are the equivalent energies of the ray before and after placing the sample respectively, l is the integral microelement along the straight line where the ray is located, L(u, β) represents the straight line equation from the focus of the ray source to detector unit u, β is the sampling angle, y represents the coordinate of the point on the straight line, and μ(y) is the linear attenuation coefficient of the sample; Dark field correction and gain correction are performed using Equation (1-2);

[0052] where Z0 is the dark field data collected by the detector when the ray source is not exposed.

[0053]

[0054] The second data correction subunit uses Equation (1-3) to perform negative logarithm transformation on the corrected scan data, and patches and rearranges the scan data after negative logarithm transformation to obtain the original CT projection data p(u, β):

[0055] p(u, β) = -lnp0(u, β) (1-3).

[0056] Furthermore, the artifacts in the original projection image appear as vertical stripe artifacts or horizontal stripe artifacts, with or without discontinuity, and the stripe gray levels vary.

[0057] The present invention also provides an imaging device, which includes:

[0058] One or more X-ray sources;

[0059] One or more X-ray detectors;

[0060] One or more mechanical control systems;

[0061] One or more data workstations;

[0062] One or more data workstations;

[0063] Multiple application programs;

[0064] and one or more programs, wherein the one or more programs are stored in the workstation memory, and when the one or more programs are executed by the application program on the workstation, the imaging device is caused to execute the CT image circular artifact adaptive structure-preserving correction method or correction device as described above.

[0065] Due to the above technical solutions adopted by the present invention, it has the following advantages:

[0066] 1. The present invention first concentrates the stripe artifacts in the low-frequency coefficient image of the multi-scale decomposition coefficients and removes them by using a one-dimensional filtering method;

[0067] 2. The filtering method adopted by the present invention can well retain the detailed structure information of the image while removing the circular artifacts of the image;

[0068] 3. The present invention does not require parameter adjustment and can adaptively correct the circular artifacts of samples with different structural characteristics, such as samples containing smooth structures, texture structures, and piecewise constant structures. Description of the Drawings

[0069] Figure 1 It is a schematic diagram of the original CT projection image of a measured object containing a smooth structure after preprocessing;

[0070] Figure 2 is Figure 1 a schematic diagram of the final CT image directly reconstructed from the original CT projection image in, and the image still has circular artifacts;

[0071] Figure 3 In, a is a schematic diagram of the decomposition coefficient containing the low-frequency characteristics of the projection image in the multi-scale decomposition coefficient image used in step 2;

[0072] Figure 3 In, b is a schematic diagram of the decomposition coefficient containing the stripe artifact characteristics in the multi-scale decomposition coefficient image used in step 2;

[0073] Figure 4 It is a schematic diagram of the second multi-scale image reconstructed from the decomposition coefficients after filtering in step 5;

[0074] Figure 5 It is a schematic diagram of the difference between the second multi-scale image and the original CT projection image in step 4;

[0075] Figure 6 It is the final CT image in step 5;

[0076] Figure 7 It is a schematic diagram of the original CT projection image of a measured object containing a texture structure after preprocessing;

[0077] Figure 8 For Figure 7 Schematic diagram of the final CT image after directly reconstructing the original CT projection image, and the image still has ring artifacts;

[0078] Figure 9 In [reference], a is a schematic diagram of the decomposition coefficient containing the low-frequency features of the projection image in the multi-scale decomposition coefficient image used in step 2;

[0079] Figure 9 In [reference], b is a schematic diagram of the decomposition coefficient containing the stripe artifact features in the multi-scale decomposition coefficient image used in step 2;

[0080] Figure 10 Schematic diagram of the second multi-scale image reconstructed from the decomposition coefficients after filtering in step 5;

[0081] Figure 11 Schematic diagram of the difference between the second multi-scale image and the original CT projection image in step 4;

[0082] Figure 12 Is the final CT image in step 5. Specific implementation mode

[0083] The present invention will be described in detail below with reference to the drawings and embodiments.

[0084] The devices used in the embodiments of the present invention are all commonly used devices in the art without special regulations; the methods used in the present invention are all commonly used methods in the art without special regulations.

[0085] The CT image ring artifact adaptive structure-preserving correction method provided by the embodiments of the present invention includes:

[0086] Step 1, obtain the original CT projection image of the object to be measured.

[0087] The object to be measured can be understood as a sample containing a smooth structure or a texture structure, or containing both structures or multiple structures at the same time. The purpose of the present invention is to correct the ring artifacts in the CT images of samples containing various structures. Figure 1 And Figure 7 Show the characteristic shape of the ring artifact in the projection image, Figure 1 In [reference], the object to be measured contains a smooth structure, while, Figure 7 In [reference], the object to be measured contains a texture structure. Among them, the method for obtaining the original projection image includes:

[0088] Step 11, the CT data scanning formula can be expressed as:

[0089]

[0090] where I(u, β) is the detection intensity of the detector unit u for the rays with equivalent energy after placing the sample, for the rays with equivalent energy , after placing the sample, is the detection efficiency of the detector unit u at the equivalent energy of , I0(u) is the detection intensity of the detector unit u for the rays with equivalent energy E without placing the sample, δ(u, E) is the detection efficiency of the detector unit u at the equivalent energy E without placing the sample, σ(u) is the intensity of scattered photons detected by the detector unit u, E and are the equivalent energies of the rays before and after placing the sample respectively, l is the integral element along the straight line where the ray is located, L(u, β) represents the straight line equation from the ray source focus to the detector unit u, β is the sampling angle, y represents the coordinate of the point on the straight line, and μ(y) is the linear attenuation coefficient of the sample.

[0091] Step 12, perform dark field correction and gain correction on the scanned data using Equation (1-2).

[0092]

[0093] where Z0 is the dark field data collected by the detector when the ray source is not exposed.

[0094] Step 13, perform negative logarithm transformation on the scanned data corrected in Step 12, and patch and rearrange the scanned data after negative logarithm transformation to obtain the original CT projection data p(u, β).

[0095] p(u, β) = -lnp0(u, β) (1-3)

[0096] Due to the inconsistency of the response efficiency δ(u, E) of the detector unit u, when the data collected at different angles β are stacked, obvious vertical stripe-like artifacts can be observed on the image of the original projection data p(u, β) through computer visualization. When directly using the classical reconstruction algorithm, obvious ring-like artifacts will appear on the reconstructed image.

[0097] When the abscissa of the original CT projection image represents the detector unit position u and the ordinate represents the scanning angle β, the artifacts at this time appear as vertical stripe structures in the original CT projection image. When the ordinate of the original CT projection image represents the detector unit position u and the abscissa represents the scanning angle β, the artifacts at this time appear as horizontal stripe structures in the original CT projection image.

[0098] Of course, according to the different actual scanning modes, Step 11, Step 12 or Step 13 can also be selectively executed to obtain a CT projection image with a "sine" structure.

[0099] Such as Figure 2 andFigure 8 As shown, both are schematic diagrams of the final CT images directly reconstructed from the original CT projection images. Obviously, these images have severe circular artifacts. The prior art cannot effectively remove these circular artifacts, especially when the stripes in the projection images (such as Figure 1 shown) are discontinuous.

[0100] Regarding the circular artifacts, the embodiments of the present invention set the following steps:

[0101] Step 2: Perform multi-scale decomposition on the original CT projection image to obtain initial multi-scale images with different sizes.

[0102] Among them, as shown in a and b of Figure 3 and Figure 9 , the method of multi-scale decomposition (as shown in the following formula) can be the wavelet transform method, and the obtained initial multi-scale image is the basis function coefficient image. Of course, existing methods such as the Laplacian pyramid decomposition method, Curvelet multi-scale decomposition method, or Shearlet multi-scale decomposition method can also be selected.

[0103] In one embodiment, multi-scale decomposition is performed on the original projection data p(u, β) using multi-scale transformation:

[0104]

[0105] In the formula, ψ s , t (u, β) is the multi-scale transformation basis function, and s and t are the scale parameter and the translation parameter respectively. cA ψ (s, t), chD ψ (s, t), cvD ψ (s, t) and cdD ψ (s, t) are the low-frequency coefficient image, horizontal direction detail coefficient image, vertical direction detail coefficient image, and diagonal direction detail coefficient image of the multi-scale decomposition respectively. The initial multi-scale image W ψ (s, t), chD ψ (s, t), cvD ψ (s, t) and cdD ψ (s, t) together form the initial multi-scale image W ψ (s, t).

[0106] Step 3: Perform boundary-preserving filtering on the initial multi-scale image to obtain the first multi-scale image, in which the low-frequency coefficient image and the boundary information are both retained.

[0107] Among them, the boundary-preserving filtering method specifically includes:

[0108] Step 31: Use the following formula (3-1) to process the initial multi-scale image Wψ Vertical direction detail coefficient image cvD in (s, t) ψ Filter (s, t) to obtain the filtered vertical direction detail coefficient image

[0109]

[0110] Wherein, Represents a spatial domain filtering operator, H(x) is the filter kernel, and this filter kernel makes Retain both the image boundary information and the artifact information at the same time. From cA ψ (s, t), chD ψ (s, t), And cdD ψ (s, t) constitute the first multi-scale image

[0111] Step 32, obtain the boundary information g(s, t) of the filtered vertical direction detail coefficient image through the boundary detection operator Edge provided by the following formula (3-2): :[[]]END]]

[0112]

[0113] Step 33, since the artifact information presents a piecewise constant feature in , extract the artifact information R(s, t) through the artifact information extraction operator Ring provided by the following formula (3-3), and the image low-frequency coefficient image is retained in cA ψ (s, t):

[0114]

[0115] The boundary-preserving filtering method adopted in step 3 includes, but is not limited to, spatial domain filtering and frequency domain filtering, as well as filtering methods based on optimization models and deep learning, and is also not limited to one-dimensional filtering and two-dimensional filtering, which significantly improves the artifact removal effect. In the existing methods, non-boundary-preserving filtering is usually adopted, and non-boundary-preserving filtering will introduce new artifacts.

[0116] Step 4, remove the frequency information corresponding to the artifact information R(s, t) in the first multi-scale image through the artifact information correction method provided by the following formula (4-1), and output the second multi-scale image This second multi-scale image Contains the low-frequency coefficient image cA ψ (s, t), the horizontal direction detail coefficient image chD ψ (s, t), the vertical direction detail coefficient image after removing the artifact information R(s, t) and diagonal detail coefficient image cdD ψ (s, t);

[0117]

[0118] It should be noted that "artifact component correction" is a method designed specifically for ring artifact characteristics. This embodiment uses simple artifact component correction to accurately remove the frequency information corresponding to the ring artifact, so the correction effect is ideal. This method is particularly suitable for situations where there are discontinuous vertical stripe artifacts in the original projection image, such as Figure 5 and Figure 11 As shown, the existing method cannot distinguish the frequency information of the ring artifact and the normal image, so when removing the ring artifact, the frequency information of a part of the normal image is also removed, resulting in loss of image information.

[0119] The multi-scale image in the above embodiment refers to a basis function coefficient image obtained by using a multi-scale decomposition method. The basis function includes but is not limited to basis functions in existing literature, and also includes basis functions constructed according to relevant principles.

[0120] Step 5: Transform the second multi-scale image Reconstruct and obtain the projection image with ring artifacts removed Then Reconstruction is performed to obtain the final CT image with ring artifacts removed, such as Figure 6 and Figure 12 As shown. The “reconstruction” method is the inverse process of the multi-scale decomposition process provided in step 2, and the second multi-scale image is obtained, as shown in Figure 4 and Figure 10 shown.

[0121] The method of the present invention does not depend on the geometric shape and internal structural characteristics of the scanned sample and does not require other auxiliary hardware. In addition, the present invention can be applied to correct other imaging systems based on multi-detector units, such as satellite remote sensing images, spectral radiation measurements, synchrotron radiation tomography images, etc.

[0122] The method of the present invention demonstrates excellent correction effects on samples with diverse structural features, such as smooth, textured, and piecewise constant structures. It is also applicable to a variety of CT scanning systems, including parallel-beam, fan-beam, and cone-beam systems. Furthermore, the invention can also be used to correct image artifacts in other imaging systems, such as satellite remote sensing images, spectral radiometry, and synchrotron radiation tomography images.

[0123] An embodiment of the present invention further provides a device for adaptively correcting ring artifacts in CT images, which includes:

[0124] An original CT projection image acquisition unit, which is used to acquire the original CT projection image of the object to be measured;

[0125] A multi-scale decomposition unit, which is used to perform multi-scale decomposition on the original CT projection image to obtain initial multi-scale images W ψ (s, t), and the initial multi-scale image W ψ (s, t) includes a low-frequency coefficient image cA ψ (s, t), a horizontal direction detail coefficient image chD ψ (s, t), a vertical direction detail coefficient image cvD ψ (s, t), and a diagonal direction detail coefficient image cdD ψ (s, t);

[0126] A boundary-preserving filtering unit, which is used to perform boundary-preserving filtering on each initial multi-scale image W ψ (s, t) and output a first multi-scale image This first multi-scale image includes a low-frequency coefficient image cA ψ (s, t), a horizontal direction detail coefficient image chD ψ (s, t), the boundary-preserving filtered vertical direction detail coefficient image and a diagonal direction detail coefficient image cdD ψ (s, t), where it contains boundary information and artifact information, locates the boundary information g(s, t) through a boundary detection operator, and extracts the artifact information R(s, t);

[0127] Among them, the boundary-preserving filtering unit specifically includes:

[0128] A filtering subunit, which is used to filter the vertical direction detail coefficient image cvD in the initial multi-scale image W ψ (s, t) by using the following formula (3-1) to obtain the filtered vertical direction detail coefficient image ψ (s, t)

[0129]

[0130] In the formula, ψ is a multi-scale transform basis function, s and t are respectively a scale parameter and a translation parameter, represents a spatial domain filtering operator, and H(s, t) is a filtering kernel;

[0131] A boundary detection subunit, which is used to obtain the boundary information g(s, t) of the filtered vertical direction detail coefficient image through the boundary detection operator Edge provided by the following formula (3-2): of:

[0132]

[0133] An artifact information extraction subunit, which is used to extract the artifact information R(s, t) through the artifact information extraction operator Ring provided by the following formula (3-3):

[0134]

[0135] An artifact correction unit, which is used to remove the frequency information corresponding to the artifact information R(s, t) in the first multi-scale image through the artifact information correction method provided by the following formula (4-1), and output a second multi-scale image in which the artifact information R(s, t) is located, and output a second multi-scale image This second multi-scale image includes a low-frequency coefficient image cA ψ (s, t), a horizontal direction detail coefficient image chD ψ (s, t), and a vertical direction detail coefficient image after removing the artifact information R(s, t) and a diagonal direction detail coefficient image cdD ψ (s, t);

[0136]

[0137] A reconstruction unit, which is used to reconstruct and rebuild the second multi-scale image in sequence to obtain a CT image without ring artifacts.

[0138] In one embodiment, the original CT projection image acquisition unit specifically includes:

[0139] A first data correction subunit, which represents the scanning data acquisition process using formula (1-1);

[0140]

[0141] In the formula, I(u, β) is the detection intensity of the ray with equivalent energy by the detector unit u after placing the sample, is the detection efficiency of the detector unit u at the equivalent energy after placing the sample, is the detection efficiency of the detector unit u at the equivalent energy when placing the sample, I0(u) is the detection intensity of the ray with equivalent energy E by the detector unit u without placing the sample, δ(u, E) is the detection efficiency of the detector unit u at the equivalent energy E without placing the sample, σ(u) is the scattered photon intensity detected by the detector unit u, E and ​They are the equivalent energies of the rays before and after placing the sample, \(l\) is the integral microelement along the straight line where the ray is located, \(L(u, \beta)\) represents the straight-line equation from the focal point of the ray source to the detector unit \(u\), \(\beta\) is the sampling angle, \(y\) represents the coordinate of the point on the straight line, and \(\mu(y)\) is the linear attenuation coefficient of the sample;

[0142] Perform dark-field correction and gain correction using Equation (1-2);

[0143]

[0144] Among them, \(Z_0\) is the dark-field data collected by the detector when the ray source is not exposed.

[0145] The second data correction subunit, which performs negative logarithm transformation on the corrected scan data using Equation (1-3), and patches and rearranges the scan data after negative logarithm transformation to obtain the original CT projection data \(p(u, \beta)\):

[0146] \(p(u, \beta)=-\ln p_0(u, \beta)\) (1-3).

[0147] In one embodiment, the artifacts of the original projection image appear as vertical stripe artifacts or horizontal stripe artifacts, with or without discontinuities, and the stripe gray-scale changes vary.

[0148] An embodiment of the present invention also provides an imaging device, which includes:

[0149] One or more X-ray sources;

[0150] One or more X-ray detectors;

[0151] One or more mechanical control systems;

[0152] One or more data workstations;

[0153] Multiple application programs;

[0154] And one or more programs, where the one or more programs are stored in the workstation memory, and when the one or more programs are executed by the application programs on the workstation, the imaging device is enabled to execute the CT image circular artifact adaptive structure-preserving correction method or correction device as described in the above embodiments.

[0155] The present invention can be applied to the correction of circular artifacts in medical and industrial X-ray CT imaging devices, and can also be used to correct other imaging systems based on multiple detector units, such as satellite remote sensing images, spectral radiometry, synchrotron radiation tomography CT images, etc.

[0156] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than limiting it. Those of ordinary skill in the art should understand that the technical solutions described in the foregoing embodiments can be modified, or some of the technical features can be equivalently replaced; these modifications or replacements 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. An adaptive structure-preserving correction method for CT image circular artifacts, characterized in that Including: Step 1: Obtain the original CT projection image of the object to be measured; Step 2: Perform multi-scale decomposition on the original CT projection image to obtain initial multi-scale images Wψ(s, t) of different sizes. The initial multi-scale image Wψ(s, t) includes a low-frequency coefficient image cAψ(s, t), a horizontal direction detail coefficient image chDψ(s, t), a vertical direction detail coefficient image cvDψ(s, t), and a diagonal direction detail coefficient image cdDψ(s, t); Step 3, perform boundary-preserving filtering on the initial multi-scale image \(W_{\psi}(s,t)\) and output the first multi-scale image This first multi-scale image includes the low-frequency coefficient image \(cA_{\psi}(s,t)\), the horizontal direction detail coefficient image \(chD_{\psi}(s,t)\), and the boundary-preserving filtered vertical direction detail coefficient image and the diagonal direction detail coefficient image \(cdD_{\psi}(s,t)\), where it contains boundary information and artifact information. Locate the boundary information \(g(s,t)\) through a boundary detection operator and extract the artifact information \(R(s,t)\); Among them, the boundary-preserving filtering method specifically includes: Step 31, filter the vertical direction detail coefficient image cvDψ(s, t) in the initial multi-scale image Wψ(s, t) using the following formula (3-1) to obtain the filtered vertical direction detail coefficient image where ψ is the multi-scale transform basis function, s and t are the scale parameter and the translation parameter respectively, denotes the spatial domain filtering operator, and H(s, t) is the filter kernel; Step 32: Obtain the boundary information g(s, t) of the filtered vertical direction detail coefficient image by using the edge detection operator Edge provided by the following formula (3-2): ​ Step 33: Use the artifact information extraction operator Ring provided by the following formula (3-3) to extract artifact information R(s, t): Step 4: Remove the frequency information corresponding to the artifact information R(s, t) in the first multi-scale image through the artifact information correction method provided by the following formula (4-1), and output a second multi-scale image The second multi-scale image contains a low-frequency coefficient image cAψ(s, t), a horizontal direction detail coefficient image chDψ(s, t), and a vertical direction detail coefficient image after removing the artifact information R(s, t), and a diagonal direction detail coefficient image cdDψ(s, t); Step 5, the second multi-scale image is reconstructed and rebuilt in sequence to obtain a CT image with ring artifacts removed.

2. The CT image annular artifact adaptive structure-preserving correction method according to claim 1, wherein The method for obtaining the original projection image in Step 1 includes: Step 11: Use formula (1-1) to represent the process of obtaining scanning data; where I(u, β) is the detection intensity of the detector unit u for the rays with equivalent energy after placing the sample of the rays, is the detection efficiency of the detector unit u at the equivalent energy of after placing the sample, I0(u) is the detection intensity of the detector unit u for the rays with equivalent energy E without placing the sample, δ(u, E) is the detection efficiency of the detector unit u at the equivalent energy E without placing the sample, σ(u) is the intensity of the scattered photons detected by the detector unit u, E and are the equivalent energies of the rays before and after placing the sample respectively, l is the integral element along the straight line where the rays are located, L(u, β) represents the straight line equation from the focus of the ray source to the detector unit u, β is the sampling angle, y represents the coordinate of the point on the straight line, and μ(y) is the linear attenuation coefficient of the sample; Step 12: Perform dark field correction and gain correction using formula (1-2) Among them, Z0 is the dark field data collected by the detector when the X-ray source is not exposed; Step 13: Perform negative logarithm transformation on the scanning data corrected in Step 12 using formula (1-3), and patch and rearrange the scanning data after negative logarithm transformation to obtain the original CT projection data p(u, β): p(u, β) = -lnp0(u, β) (1-3).

3. The CT image annular artifact adaptive structure-preserving correction method according to claim 1 or 2, characterized in that, The artifacts in the original projection image appear as vertical stripe artifacts or horizontal stripe artifacts, with or without breaks, and the gray levels of the stripes vary.

4. An adaptive structure-preserving correction device for CT image circular artifacts, characterized in that, Including: An original CT projection image acquisition unit, which is used to obtain the original CT projection image of the object to be measured; A multi-scale decomposition unit, which is used to perform multi-scale decomposition on the original CT projection image to obtain initial multi-scale images Wψ(s, t) of different sizes. The initial multi-scale image Wψ(s, t) includes a low-frequency coefficient image cAψ(s, t), a horizontal direction detail coefficient image chDψ(s, t), a vertical direction detail coefficient image cvDψ(s, t), and a diagonal direction detail coefficient image cdDψ(s, t); A boundary-preserving filtering unit, which is used to perform boundary-preserving filtering on each initial multi-scale image Wψ(s, t) and output a first multi-scale image This first multi-scale image Includes a low-frequency coefficient image cAψ(s, t), a horizontal direction detail coefficient image chDψ(s, t), a boundary-preserving filtered vertical direction detail coefficient image And a diagonal direction detail coefficient image cdDψ(s, t), where Contains boundary information and artifact information, locates the boundary information g(s, t) through a boundary detection operator, and extracts the artifact information R(s, t); Among them, the boundary-preserving filtering unit specifically includes: A filtering sub-unit, which is used to filter the vertical direction detail coefficient image cvDψ(s, t) in the initial multi-scale image Wψ(s, t) by using the following formula (3-1) to obtain the filtered vertical direction detail coefficient image Where ψ is the multi-scale transform basis function, s and t are the scale parameter and the translation parameter respectively, represents the spatial domain filtering operator, and H(s, t) is the filter kernel; A boundary detection subunit, which is used to obtain the boundary information g(s, t) of the filtered vertical direction detail coefficient image through the boundary detection operator Edge provided by the following formula (3-2): ​ An artifact information extraction subunit, which is used to extract artifact information R(s, t) using the artifact information extraction operator Ring provided by the following formula (3-3): An artifact correction unit, which is used to remove the frequency information corresponding to the artifact information R(s, t) in the first multi-scale image through the artifact information correction method provided by the following formula (4-1), and outputs a second multi-scale image The second multi-scale image contains a low-frequency coefficient image cAψ(s, t), a horizontal direction detail coefficient image chDψ(s, t), and a vertical direction detail coefficient image after removing the artifact information R(s, t), and a diagonal direction detail coefficient image cdDψ(s, t); Reconstruction unit, which is used to reconstruct and then reconstruct the second multi-scale image in sequence to obtain a CT image with ring artifacts removed. Reconstruct and then reconstruct in sequence to obtain a CT image with ring artifacts removed.

5. The CT image annular artifact adaptive structure-preserving correction device according to claim 4, wherein The original CT projection image acquisition unit specifically includes: A first data correction subunit, which uses formula (1-1) to represent the process of obtaining scanning data; where I(u, β) is the detection intensity of the detector unit u for the rays with equivalent energy after placing the sample of the rays, is the detection efficiency of the detector unit u at the equivalent energy of after placing the sample, I0(u) is the detection intensity of the detector unit u for the rays with equivalent energy E without placing the sample, δ(u, E) is the detection efficiency of the detector unit u at the equivalent energy E without placing the sample, σ(u) is the intensity of the scattered photons detected by the detector unit u, E and are the equivalent energies of the rays before and after placing the sample respectively, l is the integral element along the straight line where the ray is located, L(u, β) represents the straight line equation from the ray source focus to the detector unit u, β is the sampling angle, y represents the coordinate of the point on the straight line, and μ(y) is the linear attenuation coefficient of the sample; Perform dark field correction and gain correction using formula (1-2); Among them, Z0 is the dark field data collected by the detector when the X-ray source is not exposed; A second data correction subunit, which performs negative logarithm transformation on the corrected scanning data using formula (1-3), and patches and rearranges the scanning data after negative logarithm transformation to obtain the original CT projection data p(u, β): p(u, β) = -lnp0(u, β) (1-3).

6. The CT image annular artifact adaptive structure-preserving correction device according to claim 4 or 5, characterized in that The artifacts in the original projection image appear as vertical stripe artifacts or horizontal stripe artifacts, with or without breaks, and the gray levels of the stripes vary.

7. An imaging device, characterized in that, Including: One or more X-ray sources; One or more X-ray detectors; One or more mechanical control systems; One or more data workstations; Multiple application programs; and one or more programs, wherein the one or more programs are stored in the workstation memory, and when the one or more programs are executed by the application program on the workstation, the imaging device is caused to execute the CT image annular artifact adaptive structure-preserving correction method or correction device according to any one of claims 1-6.

Images