CN113109373B - Area array industrial CT beam hardening correction method - Google Patents
Area array industrial CT beam hardening correction method Download PDFInfo
- Publication number
- CN113109373B CN113109373B CN202110402097.9A CN202110402097A CN113109373B CN 113109373 B CN113109373 B CN 113109373B CN 202110402097 A CN202110402097 A CN 202110402097A CN 113109373 B CN113109373 B CN 113109373B
- Authority
- CN
- China
- Prior art keywords
- image
- fitting
- value
- filter
- equal
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N23/00—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
- G01N23/02—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
- G01N23/04—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material
- G01N23/046—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material using tomography, e.g. computed tomography [CT]
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2223/00—Investigating materials by wave or particle radiation
- G01N2223/03—Investigating materials by wave or particle radiation by transmission
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2223/00—Investigating materials by wave or particle radiation
- G01N2223/10—Different kinds of radiation or particles
- G01N2223/101—Different kinds of radiation or particles electromagnetic radiation
- G01N2223/1016—X-ray
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2223/00—Investigating materials by wave or particle radiation
- G01N2223/30—Accessories, mechanical or electrical features
- G01N2223/303—Accessories, mechanical or electrical features calibrating, standardising
Abstract
The invention relates to a method for correcting the hardening of a beam of an area array industrial CT (computed tomography), which comprises the following steps of 1, processing N filters with different penetration thicknesses; step 2, placing each filter segment in front of a beam outlet window of the X-ray machine in a grading manner, and performing DR scanning on each filter segment by adopting the same scanning process to obtain a DR image of each filter segment; step 3, performing exponential fitting according to the gray value of the DR image of the filter plate and the corresponding penetration thickness of the filter plate in the step 2 to obtain a fitting function; step 4, establishing a beam hardening correction function; step 5, acquiring a circumferential DR image of the detected sample by using the same scanning process as in the step 2, and performing beam hardening correction on the circumferential DR image by using a beam hardening correction function to obtain a corrected DR image; and 6, reconstructing the DR image corrected in the step 5 to obtain a CT image of the detected sample. The fitting result in the method is accurate, and the beam hardening artifact of the image after the subsequent industrial CT reconstruction is greatly reduced.
Description
Technical Field
The invention relates to the field of CT detection, in particular to a method for correcting area array industrial CT beam hardening.
Background
The industrial CT detection technology is a practical nondestructive detection means developed on the X-ray detection technology, has the advantages of visual imaging, accurate quantification, positioning and qualitative determination and the like, and is widely applied to the fields of industrial nondestructive inspection, medical treatment and health care and the like. During CT detection, X-rays irradiate a detected object from all directions, the detector is utilized to obtain the initial intensity of the detected rays and the intensity of the X-rays penetrating through the object, and then the linear attenuation coefficients of substances in all parts of the object relative to the rays are calculated through a certain reconstruction algorithm, so that a tomographic image of the detected object is obtained. However, commonly used CT image reconstruction algorithms such as filtered back projection algorithm, algebraic reconstruction algorithm, etc. all assume that the radiation source is a monochromatic source. In actual industrial CT (medical CT), an X-ray machine is used as a radiation source, and X-rays emitted therefrom have a continuous spectral distribution. When the X-ray of the continuous spectrum passes through the substance to interact with the substance, the lower-energy rays are preferentially absorbed, namely the attenuation coefficient of the low-energy rays is higher than that of the high-energy rays, so that the energy spectrum of the X-ray transmitted through the substance is changed, and the change is represented by the increase of the average energy of the X-ray passing through the substance, namely the beam hardening phenomenon. The beam hardening phenomenon seriously affects the CT detection effect, and the slice with uniform density visually shows that the brightness of the reconstructed image is not uniform, the pixel value distribution on the image shows a cup-shaped or striped artifact, and the defect judgment, detection or disease diagnosis is directly affected.
Correction methods for hardening artifacts can be divided into two main categories, namely hardware and software: the hardware method is to add some correction tools on the CT system to achieve the purpose of suppressing hardening correction, such as a water bag correction method, a filter correction method, and the like, wherein the filter correction method is the most widely used method, and the principle is to place a filter at the beam outlet of an X-ray source of the CT system, and the filter is used for pre-hardening the X-ray, filtering out the low-energy part of the energy spectrum, and allowing as many high-energy rays as possible to pass through. The method is simple and easy to implement, and extra processing is not needed to be carried out on the projection data. However, since the filter filters out the low-energy spectrum part of the radiation, a large part of photons are lost by the radiation, thereby reducing the signal-to-noise ratio. But also can only inhibit the hardening artifact phenomenon to a certain extent and can not completely solve the beam hardening problem; the software method mainly utilizes a certain correction algorithm to carry out certain post-processing on the obtained X-ray projection image according to the mechanism of forming the hardening artifact, and mainly comprises a polynomial fitting method, a Monte Carlo correction method, a dual-energy method, an iterative correction method (iteration method for short) and the like. The polynomial fitting method is most widely applied, X-ray attenuation curves formed by X-rays penetrating different thicknesses are measured by setting different step test blocks, and the ray attenuation coefficients under different thicknesses are calculated by utilizing polynomial fitting interpolation, so that the method has the advantages of relatively simple implementation and satisfactory correction effect; the iteration method is not practical in engineering application due to the defects of large operation amount, low parallelizable degree and the like; the dual-energy method has a good correction effect, but needs to perform two times of scanning under the condition of known ray energy spectrum distribution, has high requirements on hardware, is complex to implement, has low efficiency, and is not wide in practical application. Although various current calibration methods achieve certain effects theoretically or practically, calibration methods which are good in calibration effect, wide in application range and easy to implement in engineering are still needed to be further researched.
At present, filter pre-hardening and software combined correction methods are generally adopted for area array industrial CT beam hardening correction, and the correction effect is limited. In addition, the area array industrial CT cannot be provided with a post-collimator, so that the influence of scattered rays is difficult to avoid, and the stepped test block beam correction method cannot be applied to the area array CT, so that further improvement is required.
Disclosure of Invention
The invention aims to solve the technical problem of the prior art and provides an area array industrial CT beam hardening correction method which improves the correction effect by improving the fitting accuracy.
The technical scheme adopted by the invention for solving the technical problems is as follows: an area array industrial CT beam hardening correction method is characterized in that: the method comprises the following steps:
step 1, processing N filters with different penetration thicknesses; the penetration thickness of each filter is smaller than or equal to the maximum penetration thickness of the detected sample; n is a positive integer;
step 2, placing each filter segment in front of a beam outlet window of the X-ray machine in a grading manner, and performing DR scanning on each filter segment by adopting the same scanning process to obtain a DR image of each filter segment; the size of the DR image of each filter is c x d;
step 3, performing exponential fitting according to the gray value of the DR image of the filter plate and the penetration thickness corresponding to the filter plate in the step 2 to obtain a fitting function q (t), wherein t is the penetration thickness, and q (t) is a gray value obtained by fitting when the penetration thickness is t;
step 4, establishing a beam hardening correction function T (e), wherein e is a gray value of each pixel point in the acquired DR image, and T (e) is an actual penetration thickness obtained by reverse calculation when a fitting function q (t) is used;
step 5, acquiring a circumferential DR image of the detected sample by using the same scanning process as in the step 2, and performing beam hardening correction on the circumferential DR image of the detected sample by using a beam hardening correction function T (e) to obtain a corrected DR image;
and 6, reconstructing the DR image corrected in the step 5 to obtain a CT image of the detected sample.
As an improvement, the specific steps in the step 3 are as follows:
3-1, sequencing the DR images of the filters obtained in the step 2 according to the sequence from small to large of the penetration thickness of the filters to obtain
Wherein the content of the first and second substances,
is penetrated by a thickness z 1 The DR image corresponding to the filter segment of (a),
is penetrated by a thickness z 2 The DR image corresponding to the filter segment of (a),
is penetrated by a thickness z n The DR image corresponding to the filter segment of (a),
is penetrated by a thickness z N DR image corresponding to the filter of (1), z 1 ＜z 2 …＜z n …<z N ；
Step 3-2, in the above
Selecting the same pixel point (x, y), wherein x is more than or equal to 1 and less than or equal to c, and y is more than or equal to 1 and less than or equal to d; and obtain
The gray value of the same pixel point (x, y);
step 3-3, selecting the front i sequenced in the step 3-1 1 The gray value of the same pixel point in the DR image is as follows:
the gray value of the same pixel point (x, y) is more than or equal to 1 and less than or equal to i 1 N or less, combined with penetration thickness pair selection
Performing exponential function fitting on the gray value of the middle pixel point (x, y) to obtain an exponential function g obtained after the 1 st fitting 1 (t)；g 1 (t) the calculation formula is as follows:
wherein t is the penetration thickness, a 1 The amplitude obtained for the 1 st fit; b 1 Attenuation coefficient, g, for the 1 st fit 1 (t) is a gray value obtained by fitting when the penetration thickness is t;
step 3-4, calculating the absolute value k of the 1 st error 1 ；k 1 The calculation formula of (2) is as follows:
wherein f is t (x, y) is the gray value of the pixel point (x, y) in the DR image corresponding to the filter with the penetration thickness t;
3-5, calculating to obtain a 1 st difference value delta f according to the following formula t (1) (x,y)；
Δf t (1) (x,y)＝f t (x,y)-g 1 (t)/R；
Wherein R is a preset constant; t ═ z 1 、z 2 、…z N ；
Step 3-6, mixing
Substituted into Δ f t (p-1) In (x, y) gives i p Gray value, where T ≧ i p ＞ i 1 And obtaining the exponential function g obtained after the p-th fitting according to the same mode in the step 3-3 to the step 3-5 p (t), p-th error absolute value k p And the p-th difference value deltaf t (p) (x，y)；
Δf t (p) (x，y)＝Δf t (p-1) (x，y)-g p (t)/R；
Wherein, a p The amplitude obtained by the p-th fitting; b p Obtaining the attenuation coefficient for the p fitting;
the initial value of p is 2;
and 3-7, sequentially adding 1 to the value of P until P is equal to P, and the number of the gray values used in each fitting meets the following conditions: i.e. i P >i P-1 …>i p >i p-1 …>i 1 And i is P N; p is a positive integer;
respectively obtaining P exponential functions g obtained after fitting according to the method in the step 3-6 1 (t)、g 2 (t)、…g P (t); p absolute values of error k 1 、k 2 、…k P And P differences Δ f t (1) (x，y)、Δf t (2) (x，y)、…Δf t (P) (x,y)；
Step 3-8, calculating P absolute error values k in the step 3-7 1 、k 2 、、、k P Extracting the times M corresponding to the absolute value of the error if the minimum value is zero; m is an element of [1, P ]]；
When M is equal to 1, then
Then M>1, then
In the scheme, the scanning process in the step 2 comprises tube voltage, tube current, amplification ratio, focus size and integration time.
Preferably, the value range of R in the steps 3-5 and 3-6 is more than 1 and less than or equal to 2.
Compared with the prior art, the invention has the advantages that: collecting DR images of filters with different penetration thicknesses in the same scanning process, selecting DR images with different numbers for multiple times to fit an exponential function, and extracting the times corresponding to the minimum value of the error absolute value after multiple times of fitting, so as to obtain the optimal polynomial exponential fit function. Therefore, the fitting method in the method can reduce errors to the greatest extent, so that the fitting result is relatively accurate, in addition, the method has high automation degree, and the beam hardening artifact of the image after subsequent industrial CT reconstruction is greatly reduced.
Drawings
FIG. 1 is a schematic view of a sample to be tested according to an embodiment of the present invention;
2(a), 2(b), 2(c), 2(d), 2(e), 2(f) are DR images corresponding to filters having penetration thicknesses of 1mm, 2mm, 4mm, 6mm, 8mm and 12mm, respectively;
FIG. 3(a) is a DR image of a sample to be tested before beam hardening correction; FIG. 3(b) is a DR image of the sample under test after beam hardening correction;
FIG. 4 is a graph showing the variation of gray scale values in FIGS. 3(a) and 3 (b);
FIG. 5(a) is a CT image obtained after reconstruction of FIG. 3 (a); fig. 5(b) is a CT image obtained by reconstructing fig. 3 (b).
Detailed Description
The invention is described in further detail below with reference to the accompanying examples.
An area array industrial CT beam hardening correction method comprises the following steps:
step 1, processing N filters with different penetration thicknesses; the penetration thickness of each filter plate is smaller than or equal to the maximum penetration thickness of the detected sample; n is a positive integer;
step 2, placing each filter segment in front of a beam outlet window of the X-ray machine in a grading manner, and performing DR scanning on each filter segment by adopting the same scanning process to obtain a DR image of each filter segment; the size of the DR image of each filter is c x d; in the embodiment, the scanning process at least comprises parameters of tube voltage, tube current, amplification ratio, focal spot size and integration time; the parameters are all process parameters adopted by the X-ray machine during DR scanning;
step 3, performing exponential fitting according to the gray value of the DR image of the filter plate and the penetration thickness corresponding to the filter plate in the step 2 to obtain a fitting function q (t), wherein t is the penetration thickness, and q (t) is a gray value obtained by fitting when the penetration thickness is t;
the method comprises the following specific steps:
3-1, sequencing the DR images of the filters obtained in the step 2 according to the sequence from small to large of the penetration thickness of the filters to obtain
Wherein the content of the first and second substances,
is penetrated by a thickness z 1 The DR image corresponding to the filter segment of (a),
is penetrated by a thickness z 2 The DR image corresponding to the filter segment of (a),
is penetrated by a thickness z n The DR image corresponding to the filter segment of (a),
is penetrated by a thickness z N DR image corresponding to the filter of (1), z 1 ＜z 2 …＜z n …<z N ；
Step 3-2, in the above
Selecting the same pixel point (x, y), wherein x is more than or equal to 1 and less than or equal to c, and y is more than or equal to 1 and less than or equal to d; and obtain
The gray value of the same pixel point (x, y);
step 3-3, selecting the front i sequenced in the step 3-1 1 The gray value of the same pixel point in the DR image is as follows:
the gray value of the same pixel point (x, y) is more than or equal to 1 and less than or equal to i 1 N or less, combined with penetration thickness pair selection
Performing exponential function fitting on the gray value of the middle pixel point (x, y) to obtain an exponential function g obtained after the 1 st fitting 1 (t)；g 1 (t) the calculation formula is as follows:
wherein t is the penetration thickness, a 1 The amplitude obtained for the 1 st fit; b 1 Attenuation coefficient, g, for the 1 st fit 1 (t) is a gray value obtained by fitting when the penetration thickness is t;
step 3-4, calculating the absolute value k of the 1 st error 1 ；k 1 The calculation formula of (2) is as follows:
wherein f is t (x, y) is the gray value of a pixel point (x, y) in the DR image corresponding to the filter with the penetration thickness of t;
3-5, calculating to obtain a 1 st difference value delta f according to the following formula t (1) (x,y)；
Δf t (1) (x,y)＝f t (x,y)-g 1 (t)/R；
Wherein R is a preset constant; t ═ z 1 、z 2 、…z N ；
The value range of R is more than 1 and less than or equal to 2; in this example, R ═ 1.5; the value of the R is determined according to an actual experiment, so that the accuracy can be finely adjusted;
step 3-6, mixing
Substituted into Δ f t (p-1) In (x, y) gives i p Gray value, where T ≧ i p ＞ i 1 And obtaining an exponential function g obtained after the p-th fitting according to the same manner in the step 3-3 to the step 3-5 p (t), p-th error absolute value k p And the p-th difference value deltaf t (p) (x,y)；
Δf t (p) (x,y)＝Δf t (p-1) (x,y)-g p (t)/R；
Wherein, a p The amplitude obtained by the p fitting; b p Obtaining the attenuation coefficient for the p fitting;
the initial value of p is 2;
and 3-7, sequentially adding 1 to the value of P until P is equal to P, and the number of the gray values used in each fitting meets the following conditions: i.e. i P >i P-1 …>i p >i p-1 …>i 1 And i is P N; p is a positive integer; preferably, the number of gray values used in each iterative fitting is increased by the same number, i p -i p-1 ＝i P -i P-1 ；
Respectively obtaining P exponential functions g obtained after fitting according to the method in the step 3-6 1 (t)、g 2 (t)、…g P (t); p absolute values of error k 1 、k 2 、…k P And P differences Δ f t (1) (x,y)、Δf t (2) (x，y)、…Δf t (P) (x，y)；
Step 3-8, calculating P absolute error values k in the step 3-7 1 、k 2 、、、k P If the absolute value of the error is the minimum value, extracting the times M corresponding to the absolute value of the error; m is an element of [1, P ]](ii) a Obtaining an optimal polynomial exponential function, and taking the optimal polynomial exponential function as a fitting function q (t);
when M is equal to 1, then
When M is>1, then
Step 4, establishing a beam hardening correction function T (e), wherein e is a gray value of each pixel point in the acquired DR image, and T (e) is an actual penetration thickness obtained by reverse calculation when a fitting function q (t) is used; namely: actual penetration thickness t ═ q -1 (e)；
Acquiring the digit s of the DR image in the step 2, and calculating by using a fitting function q (t) to obtain the gray value of a pixel point in the DR image as [0, 2 ] s ]The penetration thickness t corresponds to each other one by one;
step 5, acquiring a circumferential DR image of the detected sample by using the same scanning process as in the step 2, and performing beam hardening correction on the circumferential DR image of the detected sample by using a beam hardening correction function T (e) to obtain a corrected DR image; wherein, the beam hardening correction function T (e) is used for carrying out beam hardening correction on the circumferential DR image of the detected sample, which is the existing method and is not described herein;
and 6, reconstructing the DR image corrected in the step 5 to obtain a CT image of the detected sample. At this time, the obtained CT image of the detected sample is the beam hardening corrected image.
The method is explained in the invention by the following experiments, taking stainless steel round piece with the diameter of 40mm as the application object, as shown in figure 1; red copper is adopted as a material of a filter (not limited to red copper), the penetration thicknesses (1mm, 2mm, 4mm, 6mm, 8mm and 12mm) of 6 kinds of red copper filters are set, the 6 kinds of filters are placed in front of a beam outlet window of an X-ray machine in times, DR scanning is carried out by adopting a fixed scanning process (tube voltage 400kV, tube current 1800uA, amplification ratio 1.215, focal size 0.4mm and integration time 1200ms), and DR images of the filters with different penetration thicknesses are collected, as shown in figure 2;
and (3) carrying out beam hardening fitting correction on each point in the DR image, randomly selecting the same point (x, y) in the DR image, and showing image gray values under different filter plate penetration thicknesses in a table 1.
TABLE 1 Filter penetration thickness and DR image gray scale corresponding table
| Serial number | Filter penetration thickness (mm) | DR image gray scale value |
| 1 | 1 | 48585.8 |
| 2 | 2 | 27492 |
| 3 | 4 | 14052.6 |
| 4 | 6 | 8956.57 |
| 5 | 8 | 6164.52 |
| 6 | 12 | 3225.98 |
Selecting the DR image gray value of a filter plate with the penetration thickness of 1mm, 2mm and 4mm during the 1 st fitting; selecting DR image gray values of filters with the penetration thicknesses of 1mm, 2mm, 4mm and 6mm during fitting for the 2 nd time; selecting DR image gray values penetrating through filters with the thicknesses of 1mm, 2mm, 4mm, 6mm and 8mm during fitting for the 3 rd time; selecting DR image gray values of filter plates with the penetration thicknesses of 1mm, 2mm, 4mm, 6mm, 8mm and 12mm during fitting for the 4 th time; r is set to 1.5, and the corresponding amplitude, attenuation and absolute value of error are obtained for each fitting according to the specific method in step 3, and the specific results are shown in table 2.
TABLE 2 iterative procedure and result data
| Number of iterations | Fitting amplitude | Fitting attenuation | Absolute value of error |
| 1 | 25192 | -0.173456 | 38822.5 |
| 2 | 15137.6 | -0.252387 | 30074.1 |
| 3 | 40790.9 | -0.728681 | 8637.77 |
| 4 | 55823.3 | -1.4201 | 2033.73 |
According to the minimum value of the absolute values of the errors in table 2, that is, the absolute value of the error in the 4 th iteration is the minimum value, the corresponding fitting function is established as follows:
the circumferential DR image of the detected sample in FIG. 1 is acquired by the same process (tube voltage 400kV, tube current 1800uA, magnification ratio 1.215, focal size 0.4mm, integration time 1200ms), and the beam hardening correction function T (e) is used to perform beam hardening correction on the original DR image to obtain the corrected DR image, as shown in FIGS. 3 and 4, for CT image reconstruction, as shown in FIG. 5.
Comparing the method with a single exponential function fitting method and a segmented cubic spline interpolation fitting method; DR images are acquired by adopting filters with penetration thicknesses of 1mm, 2mm, 3mm, 4mm, 6mm, 8mm, 10mm and 12mm respectively, and gray values of the images at any point are selected as shown in Table 3. Selecting 1mm, 2mm, 4mm, 6mm, 8mm and 12mm as known quantities, calculating gray values of two points of 3mm and 10mm by adopting different fitting methods, and calculating errors between fitting results and actual measurement results.
TABLE 3 Filter penetration thickness and DR image gray scale corresponding table
| Serial number | Filter penetration thickness (mm) | DR image gray scale value |
| 1 | 1 | 48585.8 |
| 2 | 2 | 27492 |
| 3 | 3 | 18764.5 |
| 4 | 4 | 14052.6 |
| 5 | 6 | 8956.57 |
| 6 | 8 | 6164.52 |
| 7 | 10 | 4395.67 |
| 8 | 12 | 3225.98 |
Establishing an optimal fitting result p (t) by using a single exponential function fitting method according to image gray values corresponding to 1mm, 2mm, 4mm, 6mm, 8mm and 12mm, wherein the optimal fitting result p (t) is as follows:
p(t)＝67537.2e -0.3806t ；
the optimal polynomial exponential function q (t) obtained by the method of the invention is as follows:
the gray scale values of the two points of 3mm and 10mm were calculated, and the results are shown in Table 4.
TABLE 4 comparison of results of single exponential fitting, piecewise cubic spline fitting and fitting methods of the present invention
Therefore, the error between the fitting result and the actual measurement result is the minimum. In the method, the fitting result is closer to the actual measurement result by a method of performing exponential fitting for multiple times, and the problem of larger error occurs by one-time data fitting in the existing single exponential fitting method, so that the fitting result of the method is more accurate.
The foregoing is only a preferred embodiment of the present invention, and it should be noted that, for those skilled in the art, various modifications and decorations can be made without departing from the technical principle of the present invention, and these modifications and decorations should also be regarded as the protection scope of the present invention.
Claims (3)
1. A method for correcting the hardening of an area array industrial CT beam is characterized by comprising the following steps: the method comprises the following steps:
step 1, processing N filters with different penetration thicknesses; wherein, the penetration thickness of each filter is less than or equal to the maximum penetration thickness of the detected sample; n is a positive integer;
step 2, placing each filter segment in front of a beam outlet window of the X-ray machine in a grading manner, and performing DR scanning on each filter segment by adopting the same scanning process to obtain a DR image of each filter segment; the size of the DR image of each filter is c x d;
step 3, performing exponential fitting according to the gray value of the DR image of the filter plate and the penetration thickness corresponding to the filter plate in the step 2 to obtain a fitting function q (t), wherein t is the penetration thickness, and q (t) is the gray value obtained by fitting when the penetration thickness is t;
the fitting method comprises the following specific steps:
step 3-1, sequencing the DR images of the filters obtained in the step 2 according to the sequence of the penetration thicknesses of the filters from small to large to obtain
Wherein the content of the first and second substances,
is penetrated by a thickness z 1 The DR image corresponding to the filter segment of (a),
is penetrated by a thickness z 2 The DR image corresponding to the filter segment of (a),
is penetrated by a thickness z n The DR image corresponding to the filter segment of (a),
is penetrated by a thickness z N DR image corresponding to the filter of (1), z 1 ＜z 2 ...＜z n ...＜z N ；
Step 3-2, in the above
Selecting the same pixel point (x, y), wherein x is more than or equal to 1 and less than or equal to c, and y is more than or equal to 1 and less than or equal to d; and obtain
The gray value of the same pixel point (x, y);
step 3-3, selecting the front i sequenced in the step 3-1 1 The gray value of the same pixel point in the DR image is as follows:
the gray value of the same pixel point (x, y) is more than or equal to 1 and less than or equal to i 1 N or less, combined with penetration thickness pair selection
Performing exponential function fitting on gray values of the middle pixel points (x, y) to obtain an exponential function g obtained after the 1 st fitting 1 (t)；g 1 (t) the calculation formula is as follows:
wherein t is the penetration thickness, a 1 The amplitude obtained for the 1 st fit; b 1 Attenuation coefficient, g, for the 1 st fit 1 (t) is a gray value obtained by fitting when the penetration thickness is t;
step 3-4, calculating the absolute value k of the 1 st error 1 ；k 1 The calculation formula of (2) is as follows:
wherein f is t (x, y) is the gray value of the pixel point (x, y) in the DR image corresponding to the filter with the penetration thickness t;
3-5, calculating to obtain the 1 st difference value according to the following formula
Δf t (1) (x，y)＝f t (x，y)-g 1 (t)/R；
Wherein R is a preset constant; t ═ z 1 、z 2 、...z N ；
Step 3-6, mixing
Substituted into Δ f t (p-1) In (x, y) gives i p Gray value, where T ≧ i p ＞i 1 And obtaining the exponential function g obtained after the p-th fitting according to the same mode in the step 3-3 to the step 3-5 p (t) absolute value k of p-th error p And the p-th difference value deltaf t (p) (x，y)；
Δf t (p) (x，y)＝Δf t (p-1) (x，y)-g p (t)/R；
Wherein, a p The amplitude obtained by the p-th fitting; bp is an attenuation coefficient obtained by the p-th fitting;
the initial value of p is 2;
and 3-7, sequentially adding 1 to the value of P until P is equal to P, and the number of the gray values used in each fitting meets the following conditions: i.e. i P ＞i P-1 ...＞i p ＞i p-1 ...＞i 1 And i is P N; p is a positive integer;
respectively obtaining P exponential functions g obtained after fitting according to the method in the step 3-6 1 (t)、g 2 (t)、...g P (t); p absolute values of error k 1 、k 2 、...k P And P differences Δ f t (1) (x，y)、Δf t (2) (x，y)、...Δf t (P) (x，y)；
Step 3-8, calculating P absolute error values k in the step 3-7 1 、k 2 、、、k P Extracting the times M corresponding to the absolute value of the error if the minimum value is zero; m is an element of [1, P ]]；
When M is equal to 1, then
Then M > 1, then
Step 4, establishing a beam hardening correction function T (e), wherein e is the gray value of each pixel point in the acquired DR image, and T (e) is the actual penetration thickness obtained by reverse calculation when a fitting function q (t) is used;
step 5, acquiring a circumferential DR image of the detected sample by using the same scanning process as in the step 2, and performing beam hardening correction on the circumferential DR image of the detected sample by using a beam hardening correction function T (e) to obtain a corrected DR image;
and 6, reconstructing the DR image corrected in the step 5 to obtain a CT image of the detected sample.
2. The area array industrial CT beam hardening correction method of claim 1, characterized in that: the scanning process in step 2 comprises tube voltage, tube current, amplification ratio, focal spot size and integration time.
3. The area array industrial CT beam hardening correction method of claim 1, characterized in that: in the steps 3-5 and 3-6, the value range of R is more than 1 and less than or equal to 2.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110402097.9A CN113109373B (en) | 2021-04-14 | 2021-04-14 | Area array industrial CT beam hardening correction method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110402097.9A CN113109373B (en) | 2021-04-14 | 2021-04-14 | Area array industrial CT beam hardening correction method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113109373A CN113109373A (en) | 2021-07-13 |
| CN113109373B true CN113109373B (en) | 2022-09-27 |
Family
ID=76717214
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110402097.9A Active CN113109373B (en) | 2021-04-14 | 2021-04-14 | Area array industrial CT beam hardening correction method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113109373B (en) |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101639453A (en) * | 2008-07-30 | 2010-02-03 | 中国科学院过程工程研究所 | Correcting method aiming at CT ray hardening in multi-phase flow system |
| CN107563972A (en) * | 2017-08-11 | 2018-01-09 | 重庆真测科技股份有限公司 | A kind of CT data correcting methods and a kind of ladder test block |
| CN108109183A (en) * | 2016-11-25 | 2018-06-01 | 上海东软医疗科技有限公司 | Beam hardening correction method and device |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| DE602005004653T2 (en) * | 2004-03-17 | 2009-01-29 | Philips Intellectual Property & Standards Gmbh | Correction of radiation hardening and attenuation in coherent scatter computed tomography (CSCT) |
| US7391844B2 (en) * | 2005-01-14 | 2008-06-24 | General Electric Company | Method and apparatus for correcting for beam hardening in CT images |
| JP6021319B2 (en) * | 2011-12-02 | 2016-11-09 | 東芝メディカルシステムズ株式会社 | X-ray diagnostic equipment |
-
2021
- 2021-04-14 CN CN202110402097.9A patent/CN113109373B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101639453A (en) * | 2008-07-30 | 2010-02-03 | 中国科学院过程工程研究所 | Correcting method aiming at CT ray hardening in multi-phase flow system |
| CN108109183A (en) * | 2016-11-25 | 2018-06-01 | 上海东软医疗科技有限公司 | Beam hardening correction method and device |
| CN107563972A (en) * | 2017-08-11 | 2018-01-09 | 重庆真测科技股份有限公司 | A kind of CT data correcting methods and a kind of ladder test block |
Non-Patent Citations (2)
| Title |
|---|
| A Novel Beam Hardening Correction Method for Computed Tomography;Baolei Li et al;《2007 IEEE/ICME International Conference on Complex Medical En2ineerin2》;20071231;全文 * |
| 工业CT 在工件检测中X 射线硬化校正;彭光含等;《 光谱学与光谱分析》;20080630;全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113109373A (en) | 2021-07-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Van de Casteele et al. | A model-based correction method for beam hardening artefacts in X-ray microtomography | |
| Rinkel et al. | A new method for x-ray scatter correction: first assessment on a cone-beam CT experimental setup | |
| Hogan et al. | Fluorescent computer tomography: a model for correction of X-ray absorption | |
| US8041096B2 (en) | Method for creating mass density images on the basis of attenuation images captured at different energy levels | |
| Le et al. | Least squares parameter estimation methods for material decomposition with energy discriminating detectors | |
| EP2951615B1 (en) | Method and device for generating an energy-resolved x-ray image with adapted energy threshold | |
| JP5745054B2 (en) | X-ray CT apparatus, calculation apparatus, and maintenance method for X-ray CT apparatus | |
| Yanch et al. | A comparison of deconvolution and windowed subtraction techniques for scatter compensation in SPECT | |
| Pajevic et al. | Noise characteristics of 3-D and 2-D PET images | |
| US20190213715A1 (en) | System and method for controlling noise in multi-energy computed tomography images based on spatio-spectral information | |
| US7496171B2 (en) | Method for estimating the radiation scattered in a two-dimensional detector | |
| US7471759B2 (en) | Method for estimating the scattered radiation in X-ray tomography | |
| CN110636796B (en) | Beam hardening correction in X-ray dark field imaging | |
| US20160143606A1 (en) | X-ray ct apparatus | |
| US20040109528A1 (en) | Beam hardening post-processing method and X-ray CT apparatus | |
| Bornefalk et al. | Theoretical comparison of the iodine quantification accuracy of two spectral CT technologies | |
| Ahmed et al. | A review of common beam hardening correction methods for industrial X-ray computed tomography | |
| JP2022113115A (en) | Beam hardening calibration method, x-ray ct apparatus and beam hardening calibration program | |
| Van de Casteele et al. | The effect of beam hardening on resolution in X-ray microtomography | |
| CN113125476B (en) | Area array industrial CT scattering correction method | |
| CN113109373B (en) | Area array industrial CT beam hardening correction method | |
| Cai et al. | Scatter correction using beam stop array algorithm for cone-beam CT breast imaging | |
| Kim et al. | Planar cone-beam computed tomography with a flat-panel detector | |
| Mohapatra | Development and quantitative assessment of a beam hardening correction model for preclinical micro-CT | |
| Ha et al. | Linear analysis of single-shot dual-energy computed tomography with a multilayer detector |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |