X-ray image correction method and system

Abstract

The present invention discloses an X-ray image correction method and system. The X-ray image correction method includes: simulating an X-ray energy spectrum distribution with an X-ray tube voltage energy of E; dividing the X-ray energy spectrum distribution into a plurality of monoenergetic sub-energy spectra according to the energy spectrum division step; simulating the photon flux of each sub-energy spectrum X-ray passing through a detected object with different penetration lengths L; converting the photon flux of the sub-energy spectrum into a projection value P; performing polynomial fitting on the penetration length L and the weighted projection value P of each sub-energy spectrum; and correcting the actual multi-energy projection value P using the fitting curve obtained between the penetration length L and the projection value P. m , get the ideal single energy projection value P s The X-ray image correction method of the present invention overcomes the problem of experimental conditions limited by hardening correction based on Monte Carlo simulation, has the advantages of high flexibility and high simulation accuracy, and can effectively improve the uniformity of reconstructed images.

Description

Technical Field

[0001] The present invention relates to the technical field of medical image hardening correction, and in particular to an X-ray image correction method and system. Background Art

[0002] For an ideal monoenergetic ray source, when passing through a uniform object, the ray attenuation obeys Beer's theorem,

[0003] I=I0e -μL

[0004] Transformed

[0005]

[0006] Among them, I represents the intensity of the ray after penetrating the object, I0 represents the intensity of the incident ray, μ represents the linear attenuation coefficient of the object being measured for the ray, L represents the penetration length of the ray through the object, P s It is called the single energy projection value. In the case of single energy, the single energy projection value is linearly related to the thickness of the object being inspected through the ray.

[0007] In the multi-energy case, the projection value is expressed as

[0008]

[0009] Among them, I' represents the intensity of the ray after penetrating the object, represents the average intensity of the incident ray, μ Ei Represents the linear attenuation coefficient of the object under test at a ray energy of Ei. When X-rays have a continuous energy spectrum, the attenuation of the object causes the X-ray energy spectrum to harden, enhancing its penetrating power. The projection value is nonlinearly related to the length of the ray penetrating the object under test. The projection value of multi-energy X-rays is smaller than that of single-energy X-rays, and the difference increases with the penetration length. Figure 1 Function curves of projection value and penetration length for single energy and multi-energy.

[0010] Therefore, the purpose of hardening correction is to convert the projection value P of the multi-energy group X-ray into m Projection value P when corrected to the monoenergetic case s , by establishing the projection value P m The functional relationship between L and the intermediate length L can be used to calculate the projection value P of the ideal case of single energy. s , as input for CT or CBCT image reconstruction.

[0011] Traditional beam hardening calibration methods require measuring the projection value and penetration length curve for each material tested. These curves are highly dependent on experimental conditions. Whenever conditions such as the X-ray machine voltage or the material being tested change, the curves must be remeasured to complete the hardening calibration process, making this method time-consuming and labor-intensive.

[0012] In order to break through the limitation of hardening correction on experimental conditions, this invention patent is proposed. Summary of the Invention

[0013] To address the above technical issues, the present invention provides an X-ray image correction method, system, and storage medium. The X-ray image correction method performs polynomial fitting hardening correction based on the weighted projection values ​​of each sub-spectrum after energy spectrum division and the length of the object being inspected. Specifically, the following technical solutions are adopted:

[0014] An X-ray image correction method, comprising:

[0015] Simulate the X-ray energy spectrum distribution when the X-ray tube voltage energy is E;

[0016] Divide the X-ray energy spectrum distribution into several monoenergetic sub-energy spectra according to the energy spectrum division step;

[0017] Simulate the photon flux of each sub-spectrum X-ray passing through the detected object with different penetration lengths L;

[0018] Convert the photon flux of the sub-spectrum into a projection value P;

[0019] Perform polynomial fitting on the penetration length L and the weighted projection value P of each sub-spectrum;

[0020] The actual multi-energy projection value P is corrected using the fitting curve obtained by the penetration length L and the projection value P. m , get the ideal single energy projection value P s , to achieve hardening correction of X-ray images.

[0021] As an optional embodiment of the present invention, in an X-ray image correction method of the present invention, the simulation of the photon flux of each sub-spectrum X-ray passing through the detected object with different penetration lengths L is performed using Monte Carlo simulation, including:

[0022] Geometric modeling for photon flux distribution simulation: creating a first geometric body representing a photon source, creating a second geometric body representing an object to be detected, and creating a third geometric body and a fourth geometric body representing a detector plate, wherein the second, third, and fourth geometric bodies are sequentially arranged on one side of the first geometric body, and the fourth geometric body is arranged at a position Smax away from the first geometric body, wherein Smax represents a position infinitely far from the photon source, and the thickness of the second geometric body is set to represent a penetration length L of the object to be detected;

[0023] For the physical modeling of photon flux distribution simulation, select the corresponding materials for the second, third, and fourth geometric body models. For the source modeling, set the source particle type to photon source, the photon source to a fixed point source, and set the coordinates of the point. Set the energy of the photon source to monoenergetic, and the size to the size of the several monoenergetic sub-energy spectra into which the X-ray energy spectrum distribution is divided. Set the direction parameter of the photon source to anisotropy, and the photon source emits in a single direction.

[0024] Photon flux distribution simulation counting modeling, the counting particle type is photon, the counting grid element number is the geometric model grid element number used for counting, and the number of photons emitted by the monoenergetic beam of the photon source passing through the object is counted.

[0025] As an optional embodiment of the present invention, an X-ray image correction method of the present invention saves the photon flux of each sub-energy spectrum X-ray passing through the detected object with different penetration lengths L using Monte Carlo simulation as single-energy fluxes, constructs a flux database containing each single-energy flux, and when performing hardening correction of the X-ray image, calls the corresponding single-energy flux in the flux database according to the several monoenergetic sub-energy spectra divided into, and quickly calculates the fitting curve.

[0026] As an optional embodiment of the present invention, in an X-ray image correction method of the present invention, converting the photon flux of the sub-spectrum into a projection value P includes:

[0027] The Monte Carlo simulation of each monoenergetic X-ray passing through the object to be detected obtains the photon flux, and the photon flux simulated by the monoenergetic beam is converted into the detector plate response. The detector plate response is proportional to the ray intensity, and the projection value is calculated after accumulation and normalization. The detector plate response calculation formula is:

[0028]

[0029] Where i represents the energy group, t represents the thickness, μ i Indicates the mass attenuation coefficient of the detected object corresponding to the i-th energy group, μ en，i represents the mass energy absorption coefficient of CsI corresponding to the i-th energy group, Φ i,0 represents the photon flux of the i-th energy group through the thickness t in Monte Carlo simulation, represents the average energy of the i-th energy group, I i,t It represents the energy absorbed by the CsI detector when the i-th energy group passes through a water layer with a thickness of t, I i,t Proportional to the detector plate response, the I of the i energy group i,t Perform accumulation and normalization to obtain I t , which represents the energy absorbed by the CsI detector when the weighted multi-energy X-ray beam passes through the object with a thickness of t. The calculation formula is as follows:

[0030]

[0031] The projection value is calculated using the flat panel brightness, and the formula is as follows:

[0032] Where I0 represents the incident ray intensity.

[0033] As an optional embodiment of the present invention, in an X-ray image correction method of the present invention, performing polynomial fitting on the penetration length L and the weighted projection value P of each sub-spectrum includes:

[0034] Using the projection value P as the independent variable and the material penetration length L as the dependent variable, an N-order polynomial least squares fitting is performed to obtain the functional relationship between L and P as follows:

[0035]

[0036] Among them, N is not less than 4.

[0037] As an optional embodiment of the present invention, in an X-ray image correction method of the present invention, the actual multi-energy projection value P is corrected by using the fitting curve obtained by the penetration length L and the projection value P. m , get the ideal single energy projection value P s , calculated using the following formula:

[0038] Ps = kf-1(Pm), where k is the equivalent attenuation coefficient;

[0039] As input for polynomial hardening correction of images such as CT or CBCT.

[0040] As an optional embodiment of the present invention, in an X-ray image correction method of the present invention, the simulated X-ray energy spectrum distribution with a tube voltage energy of the X-ray tube being E is simulated using Monte Carlo simulation, comprising:

[0041] Geometric modeling: creating a fifth geometric body representing the electron source, a sixth geometric body representing the tungsten target, a seventh geometric body representing the aluminum filter plate, and an eighth geometric body representing the detector plate. The sixth geometric body has an inclined target surface with a target angle θ. The fifth geometric body is located vertically above the inclined target surface of the sixth geometric body. The seventh and eighth geometric bodies are sequentially arranged horizontally to the right of the inclined target surface of the sixth geometric body.

[0042] Physical modeling: For the sixth, seventh, and eighth geometric bodies used in this simulation experiment, select the corresponding materials for material modeling. For source modeling, set the source particle type to electron, the electron source to a fixed point source and set the coordinate value of the point. Set the energy of the electron source to monoenergetic, and the electron source to emit in one direction.

[0043] Counting modeling, the counting particle type is photons, the counting gate element number is the geometric model gate element number used for counting, the photon flux is counted, and an auxiliary energy counting card is added to set the energy group for the current count.

[0044] As an optional embodiment of the present invention, in an X-ray image correction method of the present invention, the energy spectrum division step size in the step of dividing the X-ray energy spectrum distribution into a plurality of monoenergetic sub-energy spectra according to the energy spectrum division step size is determined by an equal step division method;

[0045] Optionally, the energy spectrum step size is in the range of 0.1keV-10keV.

[0046] The present invention also provides a system for implementing the X-ray image correction method, comprising:

[0047] The geometric and physical modeling module establishes the geometric model, material model and source type of X-ray irradiation of the object with different penetration lengths L;

[0048] Monte Carlo simulation module, which performs Monte Carlo particle transport based on physical modeling and calculates the projection value P by counting statistics;

[0049] Polynomial fitting module, which performs polynomial fitting on different penetration lengths L and projection values ​​P;

[0050] The image correction module uses the fitting curve obtained from the penetration length L and the projection value P to correct the actual multi-energy projection value P m , get the ideal single energy projection value P s , and thus applied to the back-projection reconstruction of X-ray images.

[0051] The present invention also provides a computer-readable storage medium storing a computer-executable program, wherein when the computer-executable program is executed, the X-ray image correction method is implemented.

[0052] Compared with the prior art, the present invention has the following beneficial effects:

[0053] The present invention provides an X-ray image correction method. First, the X-ray energy spectrum distribution of the X-ray tube is simulated. Then, the energy spectrum division step is determined to divide the energy spectrum into a number of mono-energetic X-ray beams. The Monte Carlo method is used to simulate the photon flux of the mono-energetic X-ray beams after the division step passes through the object of different thickness L and reaches the detector. Then, the photon flux calculated for each mono-energetic beam is converted into the detector flat panel response and weighted and normalized to obtain the multi-energy projection value P. m Finally, a polynomial is used to fit the functional relationship between the two. Several sets of projection values ​​close to the origin are taken to fit a straight line with the penetration length L as the independent variable. The slope k of the straight line is the equivalent attenuation coefficient, which is obtained from the equivalent monoenergetic projection P s Filtered back projection reconstructs the image. Therefore, the X-ray image correction method of the present invention overcomes the problem of experimental conditions limited by hardening correction based on Monte Carlo simulation, has the advantages of high flexibility and high simulation accuracy, and can effectively improve the uniformity of the reconstructed image. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 Comparison of the function curves of projection value and penetration length in single energy and multi-energy in the background technology of the present invention;

[0055] Figure 2 A flowchart of an X-ray image correction method according to an embodiment of the present invention;

[0056] Figure 3 The geometric modeling of X-ray energy spectrum distribution simulation in the embodiment of the present invention;

[0057] Figure 4 Geometric modeling for simulating photon flux distribution in an embodiment of the present invention;

[0058] Figure 5 Polynomial curve of projection value and penetration length at 100kV tube voltage in the embodiment of the present invention;

[0059] Figure 6 Polynomial curve of projection value and penetration length at 120kV tube voltage according to the embodiment of the present invention;

[0060] Figure 7 Grayscale value curve of row 250 in the uniformity detection diagram before and after hardening correction under 100 kV conditions in the embodiment of the present invention;

[0061] Figure 8 Gray value curve of row 250 in the uniformity detection diagram before and after hardening correction under 120kV conditions of the embodiment of the present invention;

[0062] Figure 9 Cross-sectional view of image uniformity of CBCT reconstructed at 100 kV according to an embodiment of the present invention;

[0063] Figure 10 A cross-sectional view of image uniformity of CBCT reconstructed under 120 kV according to an embodiment of the present invention;

[0064] Figure 11 The embodiment of the present invention selects a distribution example diagram of a region of interest on a CBCT image. DETAILED DESCRIPTION

[0065] To make the purpose, technical solutions and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of them.

[0066] Therefore, the following detailed description of the embodiments of the present invention is not intended to limit the scope of the claimed invention, but merely represents some embodiments of the present invention. All other embodiments derived by persons of ordinary skill in the art based on the embodiments of the present invention without creative effort shall fall within the scope of protection of the present invention.

[0067] It should be noted that, in the absence of conflict, the embodiments of the present invention and the features and technical solutions therein may be combined with each other.

[0068] It should be noted that similar reference numerals and letters denote similar items in the following drawings, and therefore, once an item is defined in one drawing, it does not need to be further defined or explained in subsequent drawings.

[0069] In the description of the present invention, it should be noted that the terms "upper" and "lower" and the like indicate orientations or positional relationships based on the orientations or positional relationships shown in the accompanying drawings, or the orientations or positional relationships in which the inventive product is typically placed when in use, or the orientations or positional relationships commonly understood by those skilled in the art. Such terms are intended solely to facilitate the description of the present invention and simplify the description, and are not intended to indicate or imply that the device or element referred to must have a specific orientation, be constructed, or operate in a specific orientation. Therefore, they should not be construed as limitations on the present invention. Furthermore, the terms "first" and "second" and the like are used solely for distinction and should not be construed as indicating or implying relative importance.

[0070] See also Figure 2 As shown, an X-ray image correction method of this embodiment includes:

[0071] Simulate the X-ray energy spectrum distribution when the X-ray tube voltage energy is E;

[0072] Divide the X-ray energy spectrum distribution into several monoenergetic sub-energy spectra according to the energy spectrum division step;

[0073] Simulate the photon flux of each sub-spectrum X-ray passing through the detected object with different penetration lengths L;

[0074] Convert the photon flux of the sub-spectrum into a projection value P;

[0075] Perform polynomial fitting on the penetration length L and the weighted projection value P of each sub-spectrum;

[0076] The actual multi-energy projection value P is corrected using the fitting curve obtained by the penetration length L and the projection value P. m , get the ideal single energy projection value P s , to achieve hardening correction of X-ray images.

[0077] In this embodiment, an X-ray image correction method first simulates the X-ray energy spectrum distribution of the X-ray tube, then determines the energy spectrum division step size to divide the energy spectrum into a number of mono-energetic X-ray beams, and uses the Monte Carlo method to simulate the photon flux of the mono-energetic X-ray beams after the division step size passing through the object of different thickness L to reach the detector. Then, the photon flux calculated for each mono-energetic beam is converted into the detector panel response and weighted and normalized to obtain the multi-energy projection value P m Finally, a polynomial is used to fit the functional relationship between the two. Several sets of projection values ​​close to the origin are taken to fit a straight line with the penetration length L as the independent variable. The slope k of the straight line is the equivalent attenuation coefficient, which is obtained from the equivalent monoenergetic projection P s Filtered back projection reconstructs the image. Therefore, the X-ray image correction method of this embodiment overcomes the problem of experimental conditions limited by hardening correction based on Monte Carlo simulation, has the advantages of high flexibility and high simulation accuracy, and can effectively improve the uniformity of the reconstructed image.

[0078] As an optional implementation manner of this embodiment, in an X-ray image correction method of this embodiment, the simulation of the photon flux of each sub-spectrum X-ray passing through the detected object with different penetration lengths L is performed using Monte Carlo simulation, including:

[0079] Geometric modeling for photon flux distribution simulation: creating a first geometric body representing a photon source, creating a second geometric body representing an object to be detected, and creating a third geometric body and a fourth geometric body representing a detector plate, wherein the second, third, and fourth geometric bodies are sequentially arranged on one side of the first geometric body, and the fourth geometric body is arranged at a position Smax away from the first geometric body, wherein Smax represents a position infinitely far from the photon source, and the thickness of the second geometric body is set to represent a penetration length L of the object to be detected;

[0080] For the physical modeling of photon flux distribution simulation, select the corresponding materials for the second, third, and fourth geometric body models. For the source modeling, set the source particle type to photon source, the photon source to a fixed point source, and set the coordinates of the point. Set the energy of the photon source to monoenergetic, and the size to the size of the several monoenergetic sub-energy spectra into which the X-ray energy spectrum distribution is divided. Set the direction parameter of the photon source to anisotropy, and the photon source emits in a single direction.

[0081] Photon flux distribution simulation counting modeling, the counting particle type is photon, the counting grid element number is the geometric model grid element number used for counting, and the number of photons emitted by the monoenergetic beam of the photon source passing through the object is counted.

[0082] Specifically, see Figure 4 As shown, the geometric modeling of the photon flux distribution simulation of this embodiment is performed by selecting the Sphere tool to create a small sphere with r = 2 cm (i.e., the first geometric body) to represent the photon source 500, and creating three cylinders (i.e., the second geometric body, the third geometric body, and the fourth geometric body). One cylinder represents the water layer 600, and the remaining two cylinders represent the carbon plate 700 and the cesium iodide detector plate 400. The outermost cesium iodide detector plate 400 is placed at a distance Smax from the photon source 500. The distance between the cesium detector plate 400 and the photon source 500 should be as far as possible, more than 10 times the distance to the other geometric bodies, and regarded as infinite. The purpose is to ignore the influence of non-normal incident photons scattered on the CsI photosensitive material. Figure 4 In this embodiment, the thickness of the water layer 600 is in the range of 0 to 20 cm.

[0083] Photon flux distribution simulation geometric modeling is used to establish a photon-sensing scene, simulating the physical path of the photon source passing through the water layer, the detector plate surface, and irradiating the CsI photosensitive material. Among them, the photon source modeling is used to simulate the photon source, the water layer modeling is used to simulate the irradiated objects of different thicknesses, the detector plate surface modeling is used to represent the outer surface of the detector plate, and the CsI modeling is used to represent the photon-sensitive material.

[0084] In the physical modeling, the corresponding materials are selected for the water layer and the detector plate geometric model. In the source modeling, the source particle type is set to photon, the source is a fixed point source, and the coordinate value of the point is set; the energy of the source particle is set to monoenergetic, and the size is the size of the calculated X-ray energy spectrum distribution divided into several monoenergetic sub-energy spectra; the direction parameter of the source is set to anisotropy, the reference direction is (0, 0, -1), the source particle is emitted in a single direction, and the cosine value of the angle between the emission direction and the reference direction is 1.

[0085] In the counting modeling, the counting particle type is photon, the counting element number is the geometric model element number used for counting, and the statistical sub-energy spectrum is used to count photons passing through the object. The T6 type of counting is selected. The T6 type is the counting type of the Monte Carlo simulation software, which represents the energy deposition of statistical photons.

[0086] As an optional implementation of this embodiment, an X-ray image correction method of this embodiment saves the photon flux of each sub-energy spectrum X-ray passing through the detected object with different penetration lengths L using Monte Carlo simulation as a single energy flux, constructs a flux database containing each single energy flux, and when performing hardening correction of the X-ray image, calls the corresponding single energy flux in the flux database according to the several monoenergetic sub-energy spectra divided, and quickly calculates the fitting curve.

[0087] Thus, this embodiment proposes an X-ray image correction method that divides a monoenergetic energy spectrum into sub-spectrums. This method uses Monte Carlo simulation of monoenergetic fluxes and saves them as a flux database for future use. During hardening correction, the "energy spectrum distribution + sub-spectrum division + database" approach allows for rapid calculation of fitting curves for application to X-ray image hardening correction. This rapid calculation of fitting curves for a given energy spectrum offers the advantages of high computational speed and adaptability. In contrast, existing techniques require Monte Carlo simulations to generate fitting curves for each energy level, a time-consuming process.

[0088] For example, if the existing Monte Carlo simulation results in a 170kV fitting curve, another Monte Carlo simulation is required to determine the fitting curve for 120kV. This method directly derives the fitting curves for X-rays at different energies by dividing the Monte Carlo simulation results for a single energy, eliminating the time-consuming Monte Carlo simulation.

[0089] In an X-ray image correction method of this embodiment, converting the photon flux of the sub-spectrum into a projection value P includes:

[0090] The Monte Carlo simulation of each monoenergetic X-ray passing through the object to be detected obtains the photon flux, and the photon flux simulated by the monoenergetic beam is converted into the detector plate response. The detector plate response is proportional to the ray intensity, and the projection value is calculated after accumulation and normalization. The detector plate response calculation formula is:

[0091]

[0092] Where i represents the energy group, t represents the thickness, μ i Indicates the mass attenuation coefficient of the detected object corresponding to the i-th energy group, μ en，i represents the mass energy absorption coefficient of CsI corresponding to the i-th energy group, Φ i，0 represents the photon flux of the i-th energy group through the thickness t in Monte Carlo simulation, represents the average energy of the i-th energy group, I i，tIt represents the energy absorbed by the CsI detector when the i-th energy group passes through a water layer with a thickness of t, I i，t Proportional to the detector plate response, the I of the i energy group i，t Perform accumulation and normalization to obtain I t , which represents the energy absorbed by the CsI detector when the weighted multi-energy X-ray beam passes through the object with a thickness of t. The calculation formula is as follows:

[0093]

[0094] The projection value is calculated using the flat panel brightness, and the formula is as follows:

[0095]

[0096] In an X-ray image correction method of this embodiment, performing polynomial fitting on the penetration length L and the weighted projection values ​​P of each sub-spectrum includes:

[0097] Using the projection value P as the independent variable and the material penetration length L as the dependent variable, an N-order polynomial least squares fitting is performed to obtain the functional relationship between L and P as follows:

[0098]

[0099] Among them, N is not less than 4, a i is the coefficient of the polynomial fit, representing the i-th power term P i The degree of contribution to the penetration length L.

[0100] In an X-ray image correction method of this embodiment, the actual multi-energy projection value P is corrected by using the fitting curve obtained by the penetration length L and the projection value P. m , get the ideal single energy projection value P s , calculated using the following formula:

[0101] Ps=kf-1(Pm)

[0102] As input for polynomial hardening correction of images such as CT or CBCT.

[0103] As an optional implementation manner of this embodiment, in an X-ray image correction method of this embodiment, the simulated X-ray energy spectrum distribution with a tube voltage energy of the X-ray tube being E is simulated using Monte Carlo simulation, including:

[0104] Geometric modeling: creating a fifth geometric body representing the electron source, a sixth geometric body representing the tungsten target, a seventh geometric body representing the aluminum filter plate, and an eighth geometric body representing the detector plate. The sixth geometric body has an inclined target surface with a target angle θ. The fifth geometric body is located vertically above the inclined target surface of the sixth geometric body. The seventh and eighth geometric bodies are sequentially arranged horizontally to the right of the inclined target surface of the sixth geometric body.

[0105] Physical modeling: For the sixth, seventh, and eighth geometric bodies used in this simulation experiment, select the corresponding materials for material modeling. For source modeling, set the source particle type to electron, the electron source to a fixed point source and set the coordinate value of the point. Set the energy of the electron source to monoenergetic, and the electron source to emit in one direction.

[0106] Counting modeling, the counting particle type is photon, the counting gate element number is the geometric model gate element number used for counting, select t4 type counting, t4 type is the counting type of Monte Carlo simulation software, t4 means statistical photon flux, add auxiliary energy counting card to set energy group for the current counting.

[0107] Specifically, see Figure 3 As shown, in this embodiment, geometric modeling first uses the Sphere tool to create a small sphere with an r value of 2 cm (i.e., the fifth geometric body) to represent the electron source 100. A tungsten target 200 with a cylindrical base and a cone as its top (i.e., the sixth geometric body) is created, where the cone's height h and base radius are r, tanθ = h / r, and θ is the target angle. Finally, two rectangular blocks (i.e., the seventh and eighth geometric bodies) are created, one representing the aluminum filter plate 300 and the other representing the cesium iodide detector plate 400. In this embodiment, the thickness of the aluminum filter plate 300 is 4.5 mm, and the distance S1 between the cesium iodide detector plate 400 and the tungsten target 200 is 750 mm.

[0108] Physical modeling begins with material modeling, selecting appropriate materials for the tungsten target, aluminum filter, and cesium iodide detector used in this simulation scenario. In source modeling, set the source particle type to electron, the source to a fixed point source, and the coordinates of the point; set the source particle energy to monoenergetic; and set the source direction parameters to anisotropic, with a reference direction of (0, 0, -1). The source particle emits in a single direction, and the cosine of the angle between the emission direction and the reference direction is 1.

[0109] In counting modeling, the counting particle type is photon, the counting element number is the geometric model element number used for counting, select t4 type counting, add an auxiliary energy counting card to set the energy group for the current counting.

[0110] In an X-ray image correction method of this embodiment, the energy spectrum division step size in dividing the X-ray energy spectrum distribution into a plurality of monoenergetic sub-energy spectra according to the energy spectrum division step size is determined by an equal step size division method.

[0111] Optionally, the energy spectrum step size is in the range of 0.1keV-10keV.

[0112] This embodiment also provides a system for implementing the X-ray image correction method, including:

[0113] The geometric and physical modeling module establishes the geometric model, material model and source model of X-ray irradiation of the object with different penetration lengths L;

[0114] Monte Carlo simulation module, which performs Monte Carlo particle transport based on physical modeling and calculates the projection value P by counting statistics;

[0115] Polynomial fitting module, which performs polynomial fitting on different penetration lengths L and projection values ​​P;

[0116] The image correction module uses the fitting curve obtained from the penetration length L and the projection value P to correct the actual multi-energy projection value P m , get the ideal single energy projection value P s , and thus applied to the back-projection reconstruction of X-ray images.

[0117] To verify the method of the present invention, the least squares method was used to fit a fourth-order polynomial with the multi-energy projection value P as the independent variable and the penetration length L as the dependent variable, and the L=f(P) curve was obtained. The cases where the 100kV and 120kV X-ray energy spectrum was divided into 0.5keV steps were studied respectively. The fitting curves obtained are shown in Fig. Figure 5 and Figure 6 The method was applied to CBCT reconstruction using projection images, and the image uniformity of the reconstructed CBCT images was evaluated. Uniformity evaluates the ability of a CBCT image to reproduce specific grayscale values ​​within the imaging area, with smaller values ​​being better.

[0118] In the module corresponding to the uniformity test, the grayscale value of the 250th row in the cross-sectional image area is taken for uniformity comparison, and a curve of the relationship between the grayscale value and the pixel position is drawn, as shown in Figure 7 、 Figure 8 As shown in the figure, the two curves show that the hardening effect of multi-spectral X-rays causes a cupping artifact in the image. This artifact is darker in areas with longer penetration lengths, and the brightness gradually decreases as the penetration length increases. The brightness of the center of the cylinder is significantly lower than that of the cylinder's edges. After hardening correction, the cupping artifact is significantly improved, and the cylinder's brightness values ​​are more uniform, which is reflected in the flatter curve after correction.

[0119] Select 5 different locations in the image area and draw circles with a diameter of 1 cm to obtain 5 areas. The distribution of the 5 areas is as follows: Figure 11As shown, the diameter of each of the five circles is 1 cm. Region E is the central ROI, and regions AD are the uniformly distributed outer ROIs. According to standard requirements for CBCT image homogeneity, the maximum absolute difference between the reference HU value of the central ROI and the average HU value of the four outer ROIs must be ≤ 30 HU.

[0120] like Figure 9 and Figure 10 The uniformity of CBCT images reconstructed at 100kV and 120kV was compared and compared with the standard for CBCT image uniformity. As can be seen from the table, the uniformity of CBCT images reconstructed after hardening correction is superior to that of uncorrected CBCT images and meets the standard for uniformity, indicating that hardening correction improves CBCT image quality.

[0121] Table 1 Comparison of uniformity of reconstructed CBCT images before and after hardening correction at 100 kV

[0122]

[0123] Table 2 Comparison of uniformity of reconstructed CBCT before and after hardening correction at 120 kV

[0124]

[0125]

[0126] This embodiment further provides a computer-readable storage medium storing a computer-executable program. When the computer-executable program is executed, the X-ray image correction method is implemented.

[0127] The computer-readable storage medium described in this embodiment may include a data signal propagated in baseband or as part of a carrier wave, which carries a readable program code. This propagated data signal may take a variety of forms, including but not limited to electromagnetic signals, optical signals, or any suitable combination of the above. The computer-readable storage medium may also be any readable medium other than a readable storage medium, which may send, propagate, or transmit a program for use by or in conjunction with an instruction execution system, device, or component. The program code contained on the computer-readable storage medium may be transmitted using any appropriate medium, including but not limited to wireless, wired, optical cable, RF, etc., or any suitable combination of the above.

[0128] This embodiment further provides an electronic device, including a processor and a memory, wherein the memory is used to store a computer executable program. When the computer program is executed by the processor, the processor executes the X-ray image correction method.

[0129] The electronic device is implemented as a general-purpose computing device. The processor may be one or multiple processors operating in concert. The present invention also does not exclude distributed processing, meaning the processors may be dispersed across different physical devices. The electronic device of the present invention is not limited to a single entity but may also be the sum of multiple physical devices.

[0130] The memory stores a computer executable program, typically a machine-readable code, which can be executed by the processor to enable the electronic device to perform the method of the present invention, or at least some of the steps in the method.

[0131] The memory includes a volatile memory, such as a random access memory unit (RAM) and / or a cache memory unit, and may also be a non-volatile memory, such as a read-only memory unit (ROM).

[0132] It should be understood that the electronic devices of the present invention may also include elements or components not shown in the above examples. For example, some electronic devices also include display units such as screens, and some electronic devices also include human-computer interaction elements such as buttons and keyboards. As long as the electronic device can execute a computer-readable program stored in its memory to implement the method of the present invention or at least some of the steps of the method, it can be considered an electronic device covered by the present invention.

[0133] Through the above description of the implementation mode, it is easy for those skilled in the art to understand that the present invention can be implemented by hardware capable of executing a specific computer program, such as the system of the present invention, and the electronic processing unit, server, client, mobile phone, control unit, processor, etc. contained in the system. The present invention can also be implemented by computer software that executes the method of the present invention, such as control software executed by a microprocessor, an electronic control unit, a client, a server, etc. However, it should be noted that the computer software that executes the method of the present invention is not limited to being executed by one or a specific hardware entity, and it can also be implemented in a distributed manner by unspecified specific hardware. For computer software, the software product can be stored in a computer-readable storage medium (which can be a CD-ROM, a USB flash drive, a mobile hard disk, etc.), or it can be distributed and stored on a network, as long as it enables an electronic device to execute the method according to the present invention.

[0134] The above embodiments are only used to illustrate the present invention and are not intended to limit the technical solutions described in the present invention. Although this specification has described the present invention in detail with reference to the above embodiments, the present invention is not limited to the above specific implementation methods. Therefore, any modification or equivalent replacement of the present invention; and all technical solutions and improvements thereof that do not depart from the spirit and scope of the invention are included in the scope of the claims of the present invention.

Claims

1. An X-ray image correction method, characterized in that: include: Simulate the X-ray energy spectrum distribution when the X-ray tube voltage energy is E; Divide the X-ray energy spectrum distribution into several monoenergetic sub-energy spectra according to the energy spectrum division step; Simulate the photon flux of each sub-spectrum X-ray passing through the detected object with different penetration lengths L; Convert the photon flux of the sub-spectrum into a projection value P; Perform polynomial fitting on the penetration length L and the weighted projection value P of each sub-spectrum; The actual multi-energy projection value P is corrected using the fitting curve obtained by the penetration length L and the projection value P. m , get the ideal single energy projection value P s , realize hardening correction of X-ray images; The simulation of the photon flux of each sub-spectrum X-ray passing through the detected object with different penetration lengths L is performed using Monte Carlo simulation, including: Geometric modeling for photon flux distribution simulation: creating a first geometric body representing a photon source, creating a second geometric body representing an object to be detected, and creating a third geometric body and a fourth geometric body representing a detector plate, wherein the second, third, and fourth geometric bodies are sequentially arranged on one side of the first geometric body, and the fourth geometric body is arranged at a position Smax away from the first geometric body, wherein Smax represents a position infinitely far from the photon source, and the thickness of the second geometric body is set to represent a penetration length L of the object to be detected; For the physical modeling of photon flux distribution simulation, select the corresponding materials for the second, third, and fourth geometric body models. For the source modeling, set the source particle type to photon source, the photon source to a fixed point source, and set the coordinates of the point. Set the energy of the photon source to monoenergetic, and the size to the size of the several monoenergetic sub-energy spectra into which the X-ray energy spectrum distribution is divided. Set the direction parameter of the photon source to anisotropy, and the photon source emits in a single direction. Photon flux distribution simulation counting modeling, the counting particle type is photon, the counting grid element number is the geometric model grid element number used for counting, and the number of photons emitted by the monoenergetic beam of the photon source passing through the object is counted; The converting of the photon flux of the sub-energy spectrum into a projection value P comprises: The Monte Carlo simulation of each monoenergetic X-ray passing through the object to be detected obtains the photon flux, and the photon flux simulated by the monoenergetic beam is converted into the detector plate response. The detector plate response is proportional to the ray intensity, and the projection value is calculated after accumulation and normalization. The detector plate response calculation formula is: Where i represents the energy group, t represents the thickness, μ i Indicates the mass attenuation coefficient of the detected object corresponding to the i-th energy group, μ en，i represents the mass energy absorption coefficient of CsI corresponding to the i-th energy group, Φ i，0 represents the photon flux of the i-th energy group through the thickness t in Monte Carlo simulation, represents the average energy of the i-th energy group, I i，t It represents the energy absorbed by the CsI detector when the i-th energy group passes through a water layer with a thickness of t, I i，t Proportional to the detector plate response, the I of the i energy group i，t Perform accumulation and normalization to obtain I t , which represents the energy absorbed by the CsI detector when the weighted multi-energy X-ray beam passes through the object with a thickness of t. The calculation formula is as follows: The projection value is calculated using the flat panel brightness, and the formula is as follows: Where I0 represents the incident ray intensity.

2. The X-ray image correction method according to claim 1, characterized in that: The photon fluxes of each sub-energy spectrum X-ray passing through the inspected object with different penetration lengths L simulated by Monte Carlo simulation are saved as monoenergetic fluxes respectively, and a flux database containing each monoenergetic flux is constructed. When performing hardening correction of X-ray images, the corresponding monoenergetic fluxes in the flux database are called according to the several monoenergetic sub-energy spectra, and the fitting curve is quickly calculated.

3. The X-ray image correction method according to claim 1, characterized in that: The polynomial fitting of the penetration length L and the weighted projection value P of each sub-spectrum includes: Using the projection value P as the independent variable and the material penetration length L as the dependent variable, an N-order polynomial least squares fitting is performed to obtain the functional relationship between L and P as follows: Among them, N is not less than 4.

4. The X-ray image correction method according to claim 3, characterized in that: The actual multi-energy projection value P is corrected by using the fitting curve obtained by obtaining the penetration length L and the projection value P. m , get the ideal single energy projection value P s , calculated using the following formula: Ps = kf-1(Pm), where k is the equivalent attenuation coefficient; As input for polynomial hardening correction of CT or CBCT images.

5. The X-ray image correction method according to claim 1, characterized in that: The simulated X-ray energy spectrum distribution with the tube voltage energy of the X-ray tube being E is simulated using Monte Carlo simulation, including: Geometric modeling: creating a fifth geometric body representing the electron source, a sixth geometric body representing the tungsten target, a seventh geometric body representing the aluminum filter plate, and an eighth geometric body representing the detector plate. The sixth geometric body has an inclined target surface with a target angle θ. The fifth geometric body is located vertically above the inclined target surface of the sixth geometric body. The seventh and eighth geometric bodies are sequentially arranged horizontally to the right of the inclined target surface of the sixth geometric body. Physical modeling: For the sixth, seventh, and eighth geometric bodies used in this simulation experiment, select the corresponding materials for material modeling. For source modeling, set the source particle type to electron, the electron source to a fixed point source and set the coordinate value of the point. Set the energy of the electron source to monoenergetic, and the electron source to emit in one direction. Counting modeling, the counting particle type is photons, the counting gate element number is the geometric model gate element number used for counting, the photon flux is counted, and an auxiliary energy counting card is added to set the energy group for the current count.

6. The X-ray image correction method according to claim 1, characterized in that: The energy spectrum division step size in dividing the X-ray energy spectrum distribution into a plurality of monoenergetic sub-energy spectra according to the energy spectrum division step size is determined by using an equal step size division method; The energy spectrum step size is in the range of 0.1keV-10keV.

7. A system for implementing the X-ray image correction method according to any one of claims 1 to 6, characterized in that: include: The geometric and physical modeling module establishes the geometric model, material model and source model of X-ray irradiation of the object with different penetration lengths L; Monte Carlo simulation module, which performs Monte Carlo particle transport based on physical modeling and calculates the projection value P by counting statistics; Polynomial fitting module, which performs polynomial fitting on different penetration lengths L and projection values ​​P; The image correction module uses the fitting curve obtained from the penetration length L and the projection value P to correct the actual multi-energy projection value P m , get the ideal single energy projection value P s , and thus applied to the back-projection reconstruction of X-ray images.

8. A computer-readable storage medium, characterized in that A computer executable program is stored, and when the computer executable program is executed, an X-ray image correction method according to any one of claims 1 to 6 is implemented.

Images