An image processing method and apparatus for a digital x-ray virtual filter

Abstract

The application discloses an image processing method and device for a digital X-ray virtual filter, which comprises the following steps: firstly, obtaining initial parameters of each pixel point in an imaging image of a digital X-ray device according to the imaging image; then, performing scattering effect compensation on the initial parameters of each pixel point in the imaging image; and finally, taking the imaging image after the scattering effect compensation as an output image of the digital X-ray device. The compensation value of the scattering effect compensation on the initial parameters of each pixel point is obtained according to a ray interference compensation mathematical model. Since the ray interference compensation mathematical model is obtained by evaluating and modeling the scattering signals of the digital X-ray device, the influence of scattered rays on the digital X-ray image can be eliminated to the maximum extent.

Description

Technical Field

[0001] The present invention relates to the technical field of digital X-ray imaging, and in particular to an image processing method and device for a digital X-ray virtual grid. Background Art

[0002] Digital X-ray imaging (DR) is an advanced X-ray imaging technology that combines computer digital image processing with X-ray radiation technology. For X-ray images, the thicker and denser the object, the higher the radiation energy required to penetrate the object and achieve effective image contrast. However, as the radiation energy increases and the density of the object increases, the radiation is more likely to scatter. When scattered X-rays hit the detector, they are superimposed on the original direct X-rays and form an image, resulting in a nonlinear increase in image grayscale and localized fog, which affects the true image contrast. Therefore, eliminating the effects of image scattering is crucial.

[0003] Currently, the primary method for suppressing X-ray scattered radiation is physical grids. This approach utilizes metal plates made of specialized materials. The specially structured radiation-attenuating material filters out scattered radiation that changes direction after penetrating the human body, preventing it from reaching the detector for imaging. Physical grids are optional and require a de-gridding image processing program with appropriate specifications. While physical grids suppress most scattered radiation, they also absorb some of the original radiation. Therefore, to achieve the same image grayscale level, the imaging dose must be increased, which in turn increases the patient's radiation exposure. On the other hand, for imaging areas where scattering is less significant, grids are unnecessary. This requires frequent grid installation and removal. Therefore, ease of use and maintenance are often a key consideration in X-ray equipment development. Sensors are often incorporated to monitor the grid's loading status, allowing physicians to easily monitor its status and prevent misoperation. These factors increase equipment development and maintenance costs. Furthermore, when using physical grids, attention must be paid to specifications such as focal length, front and back faces, and centering line, increasing the operator's workload. Compared to physical grids, digital grids can more easily eliminate the effects of image scattering. Currently, there are related image processing methods for digital grids, but most of them only infer relevant information based on the captured image and cannot fully eliminate the effects of scattered rays. Summary of the Invention

[0004] The technical problem to be solved by this application is how to eliminate the influence of scattered rays on digital X-ray images to the greatest extent possible.

[0005] According to the first aspect, an embodiment provides an image processing method for a digital X-ray virtual grid, comprising:

[0006] Acquiring initial parameters of each pixel in the imaging image according to the imaging image of the digital X-ray device; the initial parameters include pixel grayscale value, pixel depth value and / or pixel size;

[0007] Scattering effect compensation is performed on the initial parameters of each pixel point of the imaging image, and the imaging image after the scattering effect compensation is used as the output image of the digital X-ray device; wherein the compensation value for the scattering effect compensation of the initial parameters of each pixel point is obtained based on a preset ray interference compensation mathematical model, and the ray interference compensation mathematical model is obtained by evaluating and modeling the scattering signal of the digital X-ray device.

[0008] According to the second aspect, an embodiment provides a computer-readable storage medium, on which a program is stored. The program can be executed by a processor to implement the method according to the first aspect.

[0009] According to a third aspect, an embodiment provides an image processing apparatus for a digital X-ray virtual grid, configured to apply the image processing method according to the first aspect to post-process an image acquired by a digital X-ray apparatus, the image processing apparatus comprising:

[0010] An initial parameter acquisition unit is used to acquire initial parameters of each pixel in the imaging image according to the imaging image of the digital X-ray device; the initial parameters include pixel grayscale value, pixel depth value and / or pixel size;

[0011] a parameter compensation unit, configured to perform scattering effect compensation on the initial parameters of each pixel of the imaging image; wherein a compensation value for the scattering effect compensation of the initial parameters of each pixel is obtained based on a preset ray interference compensation mathematical model, and the ray interference compensation mathematical model is obtained by evaluating and modeling the scattering signal of the digital X-ray device;

[0012] An image output unit is configured to use the imaging image after the scattering effect compensation as an output image of the digital X-ray device.

[0013] According to the image processing method of the above embodiment, since the ray interference compensation mathematical model used therein is obtained by evaluating and modeling the scattered signals of the digital X-ray device, the influence of scattered rays on the digital X-ray image can be eliminated to the greatest extent. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1 1 is a flow chart of an image processing method according to an embodiment;

[0015] Figure 2A schematic diagram of the core scattering kernel superposition principle in one embodiment;

[0016] Figure 3 Schematic diagram of the Monte Carlo simulation method for simulating pencil beam X-ray scattering distribution;

[0017] Figure 4 Schematic diagram of asymmetric shift-variant scattering kernel function parameters in one embodiment;

[0018] Figure 5 A schematic diagram of an experimental measurement device for automatically optimizing scattering model parameters in one embodiment;

[0019] Figure 6 A flowchart of a deconvolution algorithm in one embodiment;

[0020] Figure 7 1 is a flowchart of a Gaussian pyramid-based multi-scale residual image BLR iterative algorithm in one embodiment;

[0021] Figure 8 A comparison diagram of the final scattering suppression effects of the shifted scattering kernel and the traditional scattering kernel function in one embodiment;

[0022] Figure 9 A comparison diagram of the effects of an LR-based deconvolution algorithm and a traditional deconvolution algorithm in one embodiment;

[0023] Figure 10 FIG. 4 is a structural block diagram of an image processing device in another embodiment. DETAILED DESCRIPTION

[0024] The present invention will be further described in detail below by means of specific embodiments in conjunction with the accompanying drawings. Similar elements in different embodiments are numbered with associated similar elements. In the following embodiments, many detailed descriptions are provided to enable the present application to be better understood. However, those skilled in the art will readily appreciate that some of the features may be omitted in different circumstances, or may be replaced by other elements, materials, or methods. In some cases, some operations related to the present application are not shown or described in the specification. This is to avoid the core portion of the present application being overwhelmed by excessive descriptions, and for those skilled in the art, it is not necessary to describe these related operations in detail. They will fully understand the related operations based on the description in the specification and the general technical knowledge in the art.

[0025] In addition, the features, operations, or characteristics described in the specification may be combined in any appropriate manner to form various embodiments. Furthermore, the steps or actions in the method description may be reordered or adjusted in a manner readily apparent to those skilled in the art. Therefore, the various sequences in the specification and drawings are provided solely for the purpose of clearly describing a particular embodiment and are not intended to be mandatory, unless otherwise specified.

[0026] Component numbers used herein, such as "first" and "second," are used solely to distinguish the components being described and do not convey any sequential or technical meaning. References to "connection" and "coupling" herein, unless otherwise specified, include both direct and indirect connections (couplings).

[0027] Traditional methods for achieving scattered radiation suppression in imaging system design include narrow beam scanning, increasing air gaps, and configuring grid hardware. These methods have some effectiveness in specific system designs, but none can completely eliminate scattered radiation interference. For example, in chest and abdomen scanning modalities, beam width reduction is difficult. Due to the large size of this part of the body, scattered radiation interference is most pronounced. Increasing air gaps is a common clinical technique, but digital X-ray devices are difficult to apply in specialized scenarios (for example, mobile DR systems used for bedside monitoring cannot employ this conventional technique). Grids are a mature and widely used technology, but their requirements for mounting location and specific source-to-detector distance make them inconvenient to use in portable DR imaging systems.

[0028] To address the limitation of digital X-ray devices that they are only suitable for specific situations, an image processing method is proposed in an embodiment of the present application. This method performs image post-processing solely by compensating for scattering effects. This method can achieve an effect approximately equivalent to that of a real hardware grid while reducing the imaging dose by 1 / 3. This frees digital X-ray devices from application environment restrictions and makes them convenient, universal, and effective.

[0029] Example 1:

[0030] Please refer to Figure 1 , is a flow chart of an image processing method in one embodiment, the image processing method is used to post-process images acquired by a digital X-ray device, and specifically includes:

[0031] Step 101: Acquire an original image.

[0032] Initial parameters of each pixel in the imaging image are acquired according to the imaging image of the digital X-ray device, wherein the initial parameters include pixel grayscale value, pixel depth value and / or pixel size.

[0033] Step 102: Scattering effect compensation.

[0034] Scattering effect compensation is performed on the initial parameters of each pixel point of the imaging image, wherein the compensation value for the scattering effect compensation of the initial parameters of each pixel point is obtained based on a preset ray interference compensation mathematical model, and the ray interference compensation mathematical model is obtained by evaluating and modeling the scattering signal of the digital X-ray device.

[0035] Step 103: output the result image.

[0036] The imaging image after scattering effect compensation is used as the output image of the digital X-ray device.

[0037] In one embodiment, the mathematical model for ray interference compensation is:

[0038] ;

[0039] ;

[0040] ;

[0041] Where x and y are the pixel coordinates in the image, is the position vector of a single pixel in the projected image, x c and y c are projection coordinates, and For The scattered ray function at the center describes the coordinates; I p is the primary ray function, which is used to represent the parameter value of each pixel point generated by the primary ray signal; s is a scattered ray function, used to represent the parameter value of each pixel point generated by the scattered ray signal; the initial parameter value of each pixel point in the imaging image is the sum of the parameter values ​​generated by the primary ray signal and the scattered ray signal for the same pixel point.

[0042] In one embodiment, the image processing method further includes:

[0043] According to the scattered ray function I s and primary ray function I p The relationship between the ray interference compensation mathematical model is solved, and the scattered ray function I s and primary ray function I p The relationship function is the product of the scattered signal intensity function and the signal distribution function, and its expression is:

[0044] ;

[0045] Where h is the scattered ray function I s and primary ray function Ip The relationship function, is the scattered signal intensity function of the small-angle scattered signal, is the scattered signal intensity function of the large-angle scattered signal.

[0046] The scattered signal intensity function expression for small-angle scattered signals is:

[0047] ;

[0048] Among them, α is the preset fitting intensity coefficient; is the scattering effect threshold cutoff function; I0 is the air projection value, and ;

[0049] The distribution function form of the scattered signal intensity function for large-angle scattered signals is:

[0050] ;

[0051] Wherein, β is the peak width coefficient obtained by fitting, and d is used to control the distribution range of the scattered ray signal; is the two-dimensional Zernike orthogonal basis function Linear superposition polynomials, c i are the polynomial coefficients.

[0052] The measured ray function I and the scattered ray function I s And the primary ray function I p Expand into vector form respectively , , , then the scattered ray function I s and primary ray function I p The following relationship exists:

[0053] ;

[0054] Among them, S jk is the element of the scattering matrix S to be obtained, which represents all For a certain pixel j composition share.

[0055] Each column element in the scattering matrix S is not a co-shift variable The discrete numerical expansion form of the measurement ray function I The expression is:

[0056] ;

[0057] The imaging matrix H = diag{1,…,1} + S is the scattering matrix plus a diagonal matrix with all elements equal to 1. The scattering matrix is ​​a block diagonal matrix with all elements non-negative. The values ​​of a large number of elements away from the main diagonal gradually decrease and tend to zero.

[0058] In one embodiment, a meta-heuristic iterative optimization algorithm of generalized pattern search is applied to calculate the corresponding scattered ray function I for the image data obtained from at least two test phantom experiments. s The error function is obtained, and the ray interference compensation mathematical model is optimized by a multi-objective function optimization method based on the wolf pack model.

[0059] The following is a theoretical explanation of the modeling process of the mathematical model for ray interference compensation (for the convenience of description, the mathematical model for ray interference compensation in the following content is collectively referred to as the scattering kernel model):

[0060] Please refer to Figure 2 , which is a schematic diagram of the core scattering and superposition principle in one embodiment, including primary rays, scattered rays and irradiated objects. The observation signal obtained in the general X-ray imaging process (i.e., the detector measurement signal) can be regarded as the primary ray I p Signal and scattered rays I s Signal superposition. For each pixel in the imaging image, there is:

[0061] I(x,y)= I s (x,y) + I p (x,y);

[0062] Where, I is the detector measurement signal of the digital X-ray imaging device, I s is the initial ray signal, I p It is a scattered ray signal, which can be understood as the detector imaging signal of the digital X-ray imaging device is the superposition of the initial ray signal and the scattered ray signal. The initial ray signal and the scattered ray signal are not independent of each other, but are correlated. The focus of the mathematical modeling of the X-ray imaging process is the estimation method of the scattered signal. In the embodiment of the present application, a scattering kernel superposition model is used to understand the scattered signal as the superposition result of multiple primary ray signals after convolution of the kernel function. The traditional scattering kernel model believes that the kernel function is in a centrally symmetric linear shift-invariant form. This assumption is derived from the method used for scatter correction in CT imaging and is applicable to the geometry of imaging systems using arc detectors. However, for DR systems using flat-panel detectors, especially when the ray cone angle is large, this assumption fails. This can be observed in Monte Carlo simulations, please refer to Figure 3Figure 1 shows a schematic diagram of a Monte Carlo simulation method for simulating pencil beam X-ray scattering distribution. After a wide cone-angle ray passes through a 20 cm deep water phantom (water cube), the resulting scattered kernel signal distribution deviates from the primary principal ray position and is no longer centrosymmetric. This deviation varies with varying water phantom depths and ray cone angles, so the scattered kernel signal is no longer linearly shift-invariant. Based on these observations, this application proposes a generalized form of the asymmetric shift-variant scattering kernel function. This adds a non-centrosymmetric component to the centrosymmetric component. This component is spatially shift-variant (i.e., each pixel in the image has a different component). This non-centrosymmetric component is expressed using Zernike basis function polynomials, achieving a complex function representation with a small number of parameters. This modeling approach more closely approximates actual physical processes, resulting in a more uniform scattering suppression effect. The imaging matrix contains all information about the primary ray and the scattered rays, establishing a direct link between the observed image (actual image) and the ideal scatter-free image (scatter-suppressed image). This computational method based on the imaging matrix overcomes the problem of using fast Fourier transforms (FFTs) to accelerate the shift-variant convolution kernel computation. Because this matrix is ​​a semi-sparse block diagonal matrix, the convolution process of the shift-variant convolution kernel can use parallel acceleration algorithms for conventional matrix operations such as GEMM, SpGEMM, SpAMM, etc. Therefore, the imaging matrix is ​​the core calculation method of this application, which enables the specific implementation of the linear shift-variant scattering kernel function modeling.

[0063] In one embodiment of the present application, a method for automatically optimizing model parameters is also disclosed. Compared with the symmetric shift-invariant scattering kernel, the number of parameters in the asymmetric shift-variant scattering kernel function increases exponentially. Due to the problems of a sharp increase in the number of measurement experiments and the improvement of the precision requirements of experimental instruments, the method of directly measuring the distribution of asymmetric shift-variant scattering kernels has become cumbersome and difficult. In one embodiment of the present application, a method combining Monte Carlo simulation and a small number of low-cost imaging experiments is proposed to automatically optimize the parameters. Among them, the parameter fitting method based on Monte Carlo simulation is to use existing Monte Carlo simulation software to simulate the physical process of pencil beam low-energy X-rays penetrating a water phantom. The energy deposition signal distribution of the secondary scattered rays generated in the process is recorded in the detector. For a specific imaging system, the ray source energy spectrum and the detector response function are calibrated. This makes the simulation results closer to the imaging system. The parameters of the shift-variant scattering kernel in the model are preliminarily fitted based on the simulation data, and the initial values ​​of the parameters are automatically established.

[0064] In one embodiment of the present application, a parameter optimization strategy based on experimental measurement data is employed. After initial parameter values ​​are established, imaging experiments are conducted within the target system. The scattered signal distribution during X-ray imaging of the phantom is measured using a continuous scanning method using a sliding lead block ray blocker array. The primary ray signal is measured using a scanning method using a sliding pencil beam limiter. These two types of data are used to further optimize the shifted scattering kernel parameters to achieve optimal results for the specific system.

[0065] In addition, the scattering kernel superposition model completes the forward modeling from the scattering-free primary ray signal to the actual imaging with scattering. Its essence is a convolution process, which belongs to the image degradation phenomenon. However, obtaining a scattering-free image from the actual imaging transformation is a corresponding inverse problem. Its solution is a deconvolution algorithm, which belongs to the image restoration problem. In one embodiment of the present application, a stable and effective deconvolution algorithm is disclosed, that is, a deconvolution solution method for a shift-varying convolution kernel is given. The deconvolution solution algorithms given in existing solutions are mostly semi-empirical methods, and their convergence and uniqueness of the convergence results have not been rigorously derived. In fact, whether it converges or not and the convergence results are related to the manual selection of the free parameters and the SPR of the imaging object, and the performance is unstable and the scattering suppression effect is limited. In addition, its output result is often accompanied by noise amplification. In the present invention, based on the LR deconvolution algorithm, a stable and strictly convergent deconvolution solution algorithm for a shift-varying convolution kernel is given. It uses a bilateral regularization term to suppress the noise of the deconvolution output image while retaining edge details. In addition, a Gaussian pyramid structure is introduced to perform deconvolution operations in residual images in different frequency domains, which effectively improves the convergence speed of the algorithm and suppresses stripe artifacts.

[0066] In one embodiment of the present application, a parallel acceleration scheme for a deconvolution algorithm is also provided. The Gaussian pyramid deconvolution algorithm proposed in the above scheme is executed sequentially layer by layer, which results in a long computation time when providing a high-quality deconvolution image. In another embodiment of the present application, another parallel acceleration method is provided in which each layer is executed sequentially. This ensures that the computation time does not increase linearly with the number of pyramid layers.

[0067] like Figure 2 As shown in the figure, it is assumed that the scattering kernel model corresponding to the scattered ray signal is obtained by convolving the primary ray signal with the kernel function related to the digital X-ray imaging device. Further simplification, under the condition of nearly parallel beam geometry imaging (i.e., small cone angle), the scattering kernel model is a linear spatial shift invariant (specially shift invariant) form, and the initial ray signal I s It can be expressed as:

[0068] ;

[0069] Initial ray signal I sThe convolution form is simplified as follows:

[0070] ;

[0071] ;

[0072] ;

[0073] in, is the position vector of a single pixel in the projected image, and For The scattering kernel function described by the center is the coordinate, and (x, y) is the pixel coordinate in the image. s and I p Solve the scattering kernel function I p , the product of the scattered signal intensity function and the signal distribution function is expressed as:

[0074] ;

[0075] Through the attenuation analysis of small-angle scattering signals, the following signal intensity function form can be obtained:

[0076] Assume I is the measurement signal, I s is the scattered signal, I p is the primary signal, and the obtained signal is I p (The real image after filtering out the scattered signal), h is I s and I p The relationship between them is expressed as:

[0077] ;

[0078] Among them, α is the intensity coefficient obtained by fitting, and is the scattering effect threshold cutoff function.

[0079] When the grayscale value of the projected image is close enough to the grayscale value of the exposed image, the scattering effect is not considered, that is, the following conditions must be met:

[0080] ; ;

[0081] Where I0 is the air projection value (the signal received by the detector in the absence of any object). ∈ is the cutoff threshold parameter, which can be adjusted between 0.9 and 1 according to the specific situation.

[0082] Please refer to Figure 4, is a schematic diagram of the parameters of an asymmetric shift-variant scattering kernel function in one embodiment. Conventional spatial shift-invariant signal distribution functions are typically centrosymmetric. However, considering the large-angle scattering signal shifting toward the central ray under large cone-angle geometry, the signal distribution no longer has a centrosymmetric distribution form. Therefore, in one embodiment of the present application, an improved scattering kernel function is proposed to address this issue. The distribution function form is:

[0083] ;

[0084] In the distribution function expression above, the first term, a two-dimensional exponential decay function, is a symmetric term representing the basis corresponding to the ideal central ray. β is the peak width coefficient obtained by fitting, and d is used to control the distribution range of the scattered signal. The second term describes the asymmetric effect introduced by deviations from the central ray and is expressed as:

[0085] ;

[0086] in, The essence of is the two-dimensional Zernike orthogonal basis function Linear superposition polynomials, c i are the polynomial coefficients.

[0087] The complete Zernike basis function system can describe any continuously varying two-dimensional distribution function. Monte Carlo simulations show that large cone angles correspond to asymmetric trends in the scattering kernel, while the distribution around the central ray is centrosymmetric. Therefore, in practical applications, using a selected N=5 basis functions is sufficient. See Table 1 below for an expression of the five Zernike basis functions used in one embodiment:

[0088] Table 1

[0089]

[0090] Considering that the distribution function at different positions rotates, the rotation angle is introduced Since the domain of the Zernike basis function is only the unit circle (i.e. 0<ρ<1), a range normalization factor d=3β is added.

[0091] In the spatially shifted scattering kernel function, the initial ray signals for different cone angle geometries all have different non-centrosymmetric scattering signal components, which vary with different pixel coordinates in the projected image. Therefore, the Zernike polynomial coefficients in the scattering kernel function are functions of the pixel coordinates. In one embodiment of the present application, this variation trend is assumed to be continuous, and a bivariate polynomial is used to model the variation trend of the Zernike polynomial coefficients with the ray cone angle, namely:

[0092] ;

[0093] In one embodiment, setting d=2, we can obtain:

[0094] ;

[0095] As above The expressions are all in ideal continuous form. The following derives the matrix linear operation expressions in discrete form.

[0096] Any linear operation can be converted into the form of matrix multiplication. The spatial shift-invariant convolution kernel can be reconstructed into a block diagonal matrix (Toeplitz matrix), which is a special circulant matrix. The convolution matrix corresponding to the spatial shift-variant convolution kernel is also in the form of a block diagonal matrix. However, since the convolution kernel corresponding to each row is different, the row elements are also different and are no longer in the form of a circulant matrix. The measurement ray function I and the scattered ray function I are converted into s And the primary ray function I p Expand into vector form respectively , , , then the scattered ray function I s and primary ray function I p The following relationship exists:

[0097] ;

[0098] Among them, S jk is the element of the scattering matrix S to be obtained, which represents all the initial signals For the scattered signal in a certain pixel j The composition of each column in the scattering matrix S is not the same shift variable scattering kernel The discrete numerical expansion form of the measured signal is The expression is:

[0099] ;

[0100] The imaging matrix (unfiltered scattered ray signals) is H = diag{1,…,1} + S, which is the scattering matrix plus a diagonal matrix with all 1s. This matrix is ​​a block diagonal matrix with all nonnegative elements. The values ​​of the elements away from the main diagonal gradually decrease and approach zero. From a numerical computational perspective, the imaging matrix is ​​a semi-sparse diagonal matrix, and parallel acceleration methods for semi-sparse diagonal matrix operations can be applied here. Convolution operations on spatially shift-invariant convolution kernels can be accelerated using fast Fourier transforms, but this strategy cannot be used for shift-variant convolution kernels. Conventional loop nesting approaches to parallel convolution acceleration have certain limitations. In particular, when accelerating using a graphics computing unit (GPU), the shared memory and thread registers are limited, restricting the sampling capacity of the convolution kernel and ultimately affecting the accuracy of the convolution calculation. (This section refers to the limitations of non-matrix acceleration.) Therefore, using the imaging matrix approach to perform purely linear operations enables the application of conventional matrix operation parallel acceleration algorithms (GEMM and SpGEMM).

[0101] In one embodiment, a step-by-step fitting approach is employed to accurately fit the asymmetric components of the large-cone-angle scattering kernel and mitigate the impact of large small-angle scattering signals on fitting accuracy. The idea is to calculate the residuals between the Monte Carlo simulation data and the main terms of the fitting function in steps, which are then used to refine the correction terms in the subsequent fitting model function, thereby automatically adjusting the parameters of the scattering kernel model.

[0102] Since the scattering kernel function is correlated with the geometry of the imaging system, X-ray energy, metal filtration, and the depth of the ray penetration into the imaging object. In the prior art, there are many defects in the experimental measurement method. In the embodiment of the present application, a method of simulating the distribution of scattered ray signals by the Monte Carlo method is adopted. Although there is a certain deviation between the Monte Carlo simulation results constructed due to the modeling deviation and the actual imaging system measurement results. However, considering the accuracy of the endogenous physical process in the Monte Carlo simulation mechanism, the obtained scattering signal distribution is a good starting point for the application of scattering kernel modeling of the real system. In order to ensure that the Monte Carlo simulation results can reflect the characteristics of a specific imaging system to a certain extent, the key parameters therein (including the ray source energy spectrum and the detector response parameters) need to be calibrated. In one embodiment of the present application, the energy deposition response function curve of the flat-panel detector is directly measured using monoenergetic synchrotron X-rays with high energy resolution at different energy levels. The specific method includes:

[0103] Step 201: Adjust the field size of the synchrotron radiation source to within 10 cm x 10 cm, and adjust the crystal monochromator to achieve the optimal energy resolution.

[0104] Step 202: Measure the detector dark field and complete the offset correction. Control the flat panel detector to face the X-ray at 90 degrees, and continuously move the detector to measure the bright field image to complete the gain correction.

[0105] Step 203: Change the synchrotron radiation beam energy (30 keV-120 keV) and the incident angle (0°, 30°, 45°, 60°, 70°, 80°), and calculate the beam flux based on the beam control ionization chamber data.

[0106] Step 204: Process the effective area of ​​the image captured by the detector and convert it to obtain an effective output value.

[0107] Step 205: normalize all data to obtain the deposition energy response curves of the detector at different energies and different incident angles.

[0108] In one embodiment, a method for calibrating an energy spectrum of a radiation source includes:

[0109] Step 301: Calculate the initial value of the energy spectrum corresponding to the energy level and the metal filtering state using a semi-empirical model.

[0110] Step 302: Measure the attenuation ratio of the X-ray source after penetrating aluminum ladders and acrylic plates of different thicknesses.

[0111] Step 303: Obtain the linear attenuation coefficients of the metal aluminum and the acrylic plate.

[0112] Step 304: Based on the attenuation measurement data of the aluminum ladder and the initial energy spectrum value obtained in step 301, the optimal energy spectrum estimate is calculated using the truncated singular value decomposition (TSVD) algorithm weighted by the initial energy spectrum. The equivalent energy spectrum is then converted based on the detector response function.

[0113] Step 305: Calculate the attenuation ratios of different acrylic plates based on the equivalent energy spectra, and calculate the deviations from the values ​​measured in step 302.

[0114] Step 306: If the deviation in step 305 is greater than the set threshold, update the parameters and return to step 304; if the deviation in step 5 is less than the set threshold, end.

[0115] In one embodiment, the Monte Carlo simulation method for scattering distribution is to use existing Monte Carlo software to simulate the scattering distribution formed by pencil beam rays penetrating water phantoms of different thicknesses. The virtual detector end records the photon flux of different energies and different incident angles. The primary particles and the secondary particles formed by Compton scattering and Rayleigh scattering are recorded as and Based on the obtained detector response function, the incident photon flux was converted to the detector output signal. The water phantom thickness ranged from 0 cm to 50 cm, with 1 cm intervals. The cone angle ranged from 0° to 15°, with 1° intervals.

[0116] The parameters of the scattering kernel function model are fitted below. The specific fitting methods include:

[0117] Step 401: performing data fitting on the linear shift-invariant parameters of the scattering kernel function model.

[0118] Define the optimization loss function for:

[0119] ;

[0120] The nonlinear multivariate function fitting was completed using the Simplex Downhill algorithm to obtain the intensity coefficient α and the peak width coefficient β.

[0121] Step 402: Fit the non-centrosymmetric parameters of the linear shift scattering kernel function to the data. This step uses the least squares method to fit the coefficients. Continuous Zernike functions are orthogonal within the unit circle domain, but this is not true for discrete Zernike functions. Here, the Gram-Schmidt orthogonalization method is used to use the basis vector Z k (k=1,…,N) Construct mutually orthogonal discrete Zernike basis vectors V k (k=1,…,N), the formula is:

[0122] ;

[0123] Among them, G is the transformation matrix in the form of a lower triangular matrix.

[0124] Orthogonal basis vectors V k The calculation method is:

[0125] ;

[0126] Transformation matrix G element g k,j The calculation method is:

[0127] ;

[0128] Among them, j=1,…,k-1, and when j=k, g k,j =1.

[0129] The residual data can be expressed as a linear superposition of discrete Zernike basis vectors, namely:

[0130] ;

[0131] The following task is used to calculate the superposition coefficient b i(i=1,…,N), specifically including:

[0132] Construct residual data baseline:

[0133] .

[0134] Then, after excluding the interference of the mean baseline, the following relationship holds:

[0135] .

[0136] Therefore, the optimization loss function used in this step is Defined as:

[0137] ;

[0138] Use the least squares principle to calculate the optimization loss function About the optimization variable b i The gradient vector is:

[0139] ,

[0140] The coefficient b can be solved i The expression is:

[0141] ;

[0142] The polynomial coefficients of the original Zernike function are then calculated as:

[0143] .

[0144] Step 403: Fitting to obtain the variation trend of the shift-variant scattering kernel parameters with different pixel coordinates. In the linear shift-variant scattering kernel model, the Zernike function coefficient c i (x c ,y c ) is a function of the projected image coordinates. In step 402, the Zernike polynomial coefficients corresponding to different non-centrosymmetric signals in discrete coordinates have been obtained. Based on this, the following bivariate polynomial coefficients describing the variation trend of these coefficients are fitted.

[0145] Using optimized loss function Defined as:

[0146] ;

[0147] According to the least squares principle, We get N sets of linear equations:

[0148] ;

[0149] in, .

[0150] vector ; and the vector ; Directly by finding the matrix singular value decomposition method:

[0151] , get the matrix pseudo inverse, and finally get the optimization coefficient:

[0152] .

[0153] In one embodiment of the present application, a parameter optimization method for a scattering kernel function model is also disclosed. The parameter optimization method optimizes model parameters based on actual measurement data of a specific imaging system, thereby improving the uniformity of scattering suppression of images in the system.

[0154] Please refer to Figure 5 This is a schematic diagram of the experimental measurement setup used to automatically optimize scattering model parameters in one embodiment. The dense scattering signal is directly measured by continuously sliding the rear ray blocker. The dense distribution of the same primary ray signal is directly measured by moving the pencil beam limiter to multiple positions. This is combined with an adaptive grid search algorithm without gradient information to continuously optimize the scattering kernel function model parameters. This forms an equivalent experimental method for indirectly measuring the true linear shift-varying scattering kernel. This also enhances the model's generalization performance in real data and improves the stability of the scattering suppression algorithm.

[0155] First, define the following intermediate function:

[0156] ;

[0157] ;

[0158] ;

[0159] ;

[0160] ;

[0161] ;

[0162] In one embodiment, multiple phantoms are used and moved to different positions. Specifically, the following steps are performed: measuring the chest phantom, abdomen phantom, upper limb phantom, and lower limb phantom respectively, obtaining corresponding signals from these four phantoms, and calculating their corresponding error functions. Ultimately, these four functions need to be optimized uniformly to extract a multi-objective optimization problem, which can be expressed as:

[0163] ;

[0164] In one embodiment, a multi-objective function optimization method based on a wolf pack model is applied, the core of which is a meta-heuristic iterative optimization algorithm for generalized pattern search. Rotate with iteration to narrow the search dimension.

[0165] In one embodiment, the Lucy-Richardson (LR) iterative algorithm is used to solve the deconvolution problem proposed in the mathematical modeling above, and certain improvements and innovations are made to address the specific problem solved in this application. The LR algorithm is deduced in the context of X-ray imaging, specifically including the following:

[0166] The process of detecting X-ray photons by a flat-panel detector is a random process that satisfies the Poisson distribution. Assume that the detector has a total of K pixels. Its expected value is:

[0167] ;

[0168] Signals are generated by capturing X-ray photons The expression is:

[0169] ;

[0170] Generate signal The posterior probability expression is:

[0171] ;

[0172] Taking the maximum likelihood estimation for the above posterior probability expression, it can be rewritten as:

[0173] ;

[0174] Without changing the maximum point of the likelihood function, the constant term is omitted. , we can get:

[0175] .

[0176] Substitute the imaging modeling expression , then:

[0177] ;

[0178] The gradient descent numerical iteration method is used to calculate the extreme value of the likelihood function. In a single-step iteration, the input image vector and the output image vector are respectively , , and the iterative equation is written as:

[0179] ;

[0180] Among them, λ is the step size parameter used to adjust the convergence speed.

[0181] Then, calculate the likelihood function gradient components separately:

[0182]

[0183] ;

[0184] Notice The elements can be written as the transposed matrix of element, then the likelihood function gradient vector expression is:

[0185] ;

[0186] in is a unit vector with the same size as the image vector, let , bringing in the gradient descent iterative equation, we can finally get the LR algorithm iterative equation:

[0187] ;

[0188] Among them, only is a matrix multiplication operation, and all other multiplication and division operations are element-wise operations. For clarity, the component form is written below:

[0189] ;

[0190] To facilitate writing the subscript of the quantity, as well as Respectively represent , The weight.

[0191] The LR algorithm iterative equation is usually used iteratively, and the number of iterations required for convergence is generally 10-20 times.

[0192] The error function used to judge convergence is usually:

[0193] ;

[0194] when When it is less than the preset convergence threshold End the iteration.

[0195] The LR algorithm is a product operation, which naturally guarantees the non-negativity of the convergence result and the convergence result is stable. However, when it processes large-scale images, the overall efficiency of the algorithm is low due to the long number of convergence iterations. Moreover, since the overall spectrum distribution of the signal in the X-ray image is relatively wide, the convergence speed of high-frequency and low-frequency signals is different. In areas with obvious edges, the Gibbs oscillation phenomenon (Gibbs Phenomenon) is prone to occur, resulting in obvious stripe artifacts in the deconvolution image. In addition, since the influence of Gaussian noise is not considered in the derivation process, the convergence effect of the algorithm is also related to the initial image effect. The invention has made a number of improvements to these problems, which are described in detail below.

[0196] Generally, improvement strategies for the LR algorithm focus on modifying the maximum likelihood function. The maximum likelihood function based on the Poisson process itself does not focus on the spatial relationship between image signals. No prior assumptions are made about multi-scale feature information and noise correlation information. By adding a regularization term to the likelihood function, the algorithm can converge to the prior assumptions preset in the regularization term. Typical regularization terms are the Tikhonov-Miller (TM) regularization term (i.e., the gradient 2-norm regularization term) and the total variation (TV) regularization term (i.e., the gradient 1-norm regularization term). They constrain local signal changes by simply limiting the overall gradient of the image to achieve the effect of noise reduction. However, this also leads to more serious problems with local smoothing of the image. Therefore, the regularization term should have the effect of preserving adaptive edge information. The following introduces a bilateral regularization constrained LR algorithm (Bilateral Lucy-Richardson, BLR). Its essence is derived from the idea of ​​bilateral filtering. The bilateral regularization expression is:

[0197] ;

[0198] Among them, Ω is the neighborhood of pixel i, and its size controls the scope of the regularization term. The function It is a scope control item.

[0199] And the function As a smoothing control term, it imposes a certain penalty effect on the signal fluctuation in the neighborhood. Because the regularization term needs to be minimized, the likelihood function is adjusted to:

[0200] ;

[0201] Then solving the deconvolution problem is transformed into a minimum value solving optimization problem, and its minimum value solving optimization formula is:

[0202] ;

[0203] Repeat the gradient descent method to find the loss function The minimum point of :

[0204] ;

[0205] Calculating the loss function Gradient vector:

[0206] ;

[0207] make:

[0208] ;

[0209] Substituting the gradient descent iterative equation into the final BLR algorithm iterative equation is:

[0210] ;

[0211] The gradient of the bilateral regularization function can be calculated independently according to the following formula:

[0212] ;

[0213] in, It is a function of the neighborhood pixel j of the deconvolution input image pixel i, and its j-th component is obtained by the following transformation:

[0214] ;

[0215] and is an image translation operator, which operates by shifting the image Follow the vector Perform translation. Vector The starting point is the coordinate of pixel position i, and the end point is the coordinate of pixel position j.

[0216] Compared with the TV or TM regularization term constrained LR algorithm, the convergence result of the BLR algorithm can retain more details. As mentioned above, the convergence speed of high and low frequency signals in the LR algorithm is different, and the situation is the same for the BLR algorithm. In order to speed up the overall convergence speed of the algorithm, the invention adopts a strategy of completing the BLR deconvolution operation in parallel by constructing a multi-scale pyramid hierarchical level. In order to adapt to the simultaneous completion of deconvolution iterative operations in image information of different scales, the invention abandons the use of the original image for calculation and uses the residual image instead. During the iterative process of the residual image, the low-power signal of the detail to be restored is no longer interfered by the strong gradient signal at the edge of the original image, so that the stripe artifacts are suppressed.

[0217] The details are explained below.

[0218] Please refer to Figure 6, which is a flow chart of the deconvolution algorithm in one embodiment. The original X-ray image is downsampled N-1 times by a factor of 2 to construct an N-layer Gaussian pyramid. The N-layer Gaussian pyramid is constructed using Gaussian filtering and then subtraction.

[0219] In one embodiment, for simplicity of explanation, N=3 is used as an example. The input images after downsampling are:

[0220] ;

[0221] ;

[0222] ;

[0223] The residual image of the corresponding level is calculated as follows:

[0224]

[0225] Then the BLR algorithm is used to obtain the deconvolution image of each layer of residual image:

[0226] ;

[0227] The final output of the complete deconvolution image expression is:

[0228] ;

[0229] Please refer to Figure 7 , is a flowchart of a Gaussian pyramid-based multi-scale residual image BLR iterative algorithm in one embodiment. Within the framework of the Gaussian pyramid, the BLR algorithm can be further optimized. In the aforementioned multi-level parallel BLR calculation process, inter-layer information is not fully exploited, so the algorithm flow is further modified to address this issue. The upsampled image after deconvolution of the previous layer is used as a guide image. To fully exploit its edge information, the invention draws on the idea of ​​an image-guided denoising algorithm to modify the BLR algorithm into a joint bilateral LR algorithm (JBLR).

[0230] The regularization term in the minimum value solution optimization formula is modified to the guide function regularization term, and its expression is:

[0231] ;

[0232] The smoothing control items are:

[0233] ;

[0234] Need to optimize , whose expression is:

[0235] ;

[0236] It can be foreseen that the iterative equations to be solved are:

[0237] ;

[0238] It should be noted that the image is guided in the iterative equation It is also constantly updated.

[0239] When the iteration converges, Converges to .

[0240] In addition, the regularization function gradient calculation also has the following changes:

[0241] ;

[0242] Following the principle of incorporating inter-layer information from the Gaussian pyramid, the residual image calculation method has also undergone some modifications. The subtracted image is no longer a Gaussian-smoothed image of the layer itself, but rather a degraded image obtained by upsampling the deconvolved image from the previous layer. Furthermore, to eliminate the smoothing effect of bilinear interpolation during the upsampling process, a simplified BLR algorithm based on Gaussian kernel convolution is added after upsampling.

[0243] Still taking N=3 as an example, the calculation methods of the residual images of each layer are listed.

[0244] To distinguish different types of convolution kernels in the BLR algorithm, the symbol BLR(·;α) is used to represent the BLR algorithm based on the convolution kernel α. The same is true for JBLR(·;α). The downsampled images of the input image are:

[0245] ;

[0246] ;

[0247] ;

[0248] The residual image of the third layer is:

[0249] ;

[0250] The deconvolution image of the third layer is:

[0251] ;

[0252] The upsampled image of the third layer is:

[0253] ;

[0254] The residual image of the second layer is:

[0255] ;

[0256] The deconvolution image of the second layer is:

[0257] ;

[0258] The upsampled image of the second layer is:

[0259] ;

[0260] The residual image of the first layer is:

[0261] ;

[0262] The deconvolution image of the first layer is:

[0263] ;

[0264] The upsampled image of the first layer is:

[0265] ;

[0266] Before the first layer outputs the final scatter-suppressed X-ray image (i.e. The convolution kernel used in the BLR step is a specially designed Gaussian scattering kernel, K. Its function is to compensate for the scattering effects of the flat-panel detector package and the internal scintillator. This completes the calculation of the deconvolution algorithm.

[0267] There are alternative methods for the Monte Carlo simulation involved in the embodiments of this application, including:

[0268] 1) Detector Response Function Calibration Method: Scintillator-related parameters are generally confidential to flat-panel detector manufacturers and difficult to obtain. However, if the chemical composition, doping elements, and geometric parameters of the detector scintillator are known, Monte Carlo simulation of the detector's response function to monoenergetic radiation of varying energies can be considered as an alternative.

[0269] 2) Radiation source energy spectrum calibration method: Radiation source energy spectrum measurement can be replaced with a GeTd gamma-ray spectrometer. This is generally used for characterization and analysis of industrial X-ray tubes. However, if used for clinical X-ray spectrum measurement below 120 kV, noise and low-energy interference signals may be problematic. Further signal calibration is required. Given the parameters of the X-ray tube anode target, internal metal filter, and electron energy spread, the Monte Carlo method can be used to simulate the process of electron beam bombardment of the anode target to generate X-rays. The phase space of the emitted photons is recorded for use in subsequent simulations.

[0270] The embodiments of the present application involve experimental instruments, and two special engineering devices are used in the model parameter optimization method, namely a lead block ray blocker array and a pencil beam limiter. There is also a lead block ray blocker array for independently measuring scattered signals. In the absence of this array tooling, a separate ray blocker can be used instead, and the same effect can be achieved in the end. However, the experimental efficiency will definitely be reduced exponentially. In addition, the pencil beam limiter is used to independently measure the primary signal, and the direction of its light hole needs to be adjusted and aligned with the ray cone angle. Tooling that does not have this function can achieve an approximate effect, but in the case of large cone angles, the larger ray penumbra will reduce the experimental accuracy. In addition, a fan-shaped beam limiter can be used as an alternative to increase the experimental efficiency. However, this will lead to a decrease in measurement accuracy.

[0271] The alternative modeling and algorithm solutions involved in the embodiments of this application include:

[0272] (1) There are many alternatives to the mathematical modeling of the scattering kernel function. In the embodiments of this application, only one feasible solution that has been verified to have practical effects is provided. Symmetric components can be replaced by the superposition of multiple Gaussian distributions, which can improve computational efficiency as a decomposable convolution kernel. Asymmetric components can be replaced by orthogonal complete real function bases such as two-dimensional B-spline functions and Bessel functions (cylindrical harmonic functions).

[0273] (2) In the optimization algorithm, there are many options for the nonlinear optimization method based on non-gradient information mentioned in the main text, which will not be described here.

[0274] (3) The image degradation problem based on the imaging matrix proposed in the embodiment has a variety of common image restoration algorithms, such as Weiner filtering, regularized inverse filtering, Landweber algorithm, Tiknonov-Miller algorithm, and FISTA (Fast Iterative soft-threshold) algorithm. Their effects vary with the adjustment of the free parameters, and similar results should be obtained in the end.

[0275] (4) The hardware involved in applying the above method specifically includes: the deconvolution algorithm can be ported to any personal computer, workstation or cloud server. The serial version can be simply completed by using the single-thread operation of a single central processing unit (CPU). The parallel version can be completed by programming an independent graphics processing unit (GPU) or a multi-threaded CPU. In addition, the embedded computing processing unit integrated in the DR system can also implement the algorithm proposed in the present invention. Possible solutions include Raspberry Pi, Field-programmable gate array (FPGA) and Nvidia's Jetson embedded system equipped with an independent GPU unit.

[0276] Please refer to Figure 8 , is a comparison diagram of the final scattering suppression effect of the shifted scattering kernel and the traditional scattering kernel function in an embodiment, FIG082 is a DR image of a human-shaped phantom without scattering suppression, FIG083 is a DR image of a child human-shaped phantom with scattering suppression (shifted scattering kernel), and FIG084 is a DR image of a child human-shaped phantom with scattering suppression (traditional scattering kernel).

[0277] Please refer to Figure 9 This is a comparison chart of the effects of the LR-based deconvolution algorithm and the traditional deconvolution algorithm in an embodiment. Figure 091 is a scatter-suppressed virtual DR image of the head (LR deconvolution), Figure 092 is a scatter-suppressed virtual DR image of the head (traditional deconvolution), Figure 093 is a scatter-free virtual DR image of the head, and Figure 094 is a scatter-containing virtual DR image of the head.

[0278] Please refer to Figure 10 , is a structural block diagram of an image processing device in another embodiment, which is used to process an image acquired by digital X-rays using the image processing method described above, and specifically includes an initial parameter acquisition unit 10, a parameter compensation unit 20, and an image output unit 30. The initial parameter acquisition unit 10 is used to acquire the initial parameters of each pixel in the imaging image based on the imaging image of the digital X-ray device. In one embodiment, the initial parameters include pixel grayscale value, pixel depth value and / or pixel size. The parameter compensation unit 20 is used to compensate the initial parameters of each pixel in the imaging image for scattering effects, wherein the compensation value for the scattering effect compensation of the initial parameters of each pixel is obtained based on a preset ray interference compensation mathematical model, and the ray interference compensation mathematical model is obtained by evaluating and modeling the scattering signal of the digital X-ray device. The image output unit 30 is used to use the imaging image after scattering effect compensation as the output image of the digital X-ray device.

[0279] The image processing method disclosed in one embodiment of this application proposes an asymmetric shifted scattering kernel that is closer to real physical processes, resulting in more accurate scattering signal estimation, a more uniform descattering effect, and no bright or dark fluctuations in the output image. This method improves the stability of scatter suppression using grid technology. Unlike traditional semi-empirical deconvolution algorithms in this type of solution, the LR algorithm used in one embodiment of this application exhibits strict convergence stability and is also effective for large-volume chest and abdominal DR images with the most severe scattering. It not only simultaneously suppresses noise artifacts in scatter-suppressed images, but also adaptively adjusts the bilateral filtering regularization term to achieve simultaneous noise suppression. Furthermore, an automatic parameter optimization method is proposed in one embodiment, enabling precise optimization of scattering model parameters using a limited amount of specialized experimental data. Parameter setting through experimental methods eliminates the need for tedious repetition (existing techniques require 30-40 imaging cycles, which are susceptible to fluctuation errors in low-intensity radiation detection). This method measures relevant signals under high radiation intensity conditions, improving measurement data accuracy while reducing the number of required measurements (parameter setting can be completed with as few as four imaging cycles).

[0280] The image processing method disclosed in the embodiments of the present application first obtains initial parameters for each pixel in an image captured by a digital X-ray device; then, scattering effect compensation is performed on these initial parameters for each pixel in the image; and finally, the image after scattering effect compensation is used as the output image of the digital X-ray device. The compensation value for scattering effect compensation of the initial parameters for each pixel is obtained based on a mathematical model for ray interference compensation. Because the mathematical model for ray interference compensation is obtained by evaluating and modeling the scattering signal of the digital X-ray device, it can minimize the impact of scattered rays on the digital X-ray image.

[0281] Those skilled in the art will appreciate that all or part of the functions of the various methods in the above embodiments can be implemented by hardware or by computer program. When all or part of the functions in the above embodiments are implemented by computer program, the program can be stored in a computer-readable storage medium, and the storage medium can include: read-only memory, random access memory, disk, optical disk, hard disk, etc., and the program is executed by a computer to implement the above functions. For example, the program is stored in the memory of the device, and when the program in the memory is executed by the processor, all or part of the above functions can be implemented. In addition, when all or part of the functions in the above embodiments are implemented by computer program, the program can also be stored in a storage medium such as a server, another computer, disk, optical disk, flash disk or mobile hard disk, and saved in the memory of the local device by downloading or copying, or the system of the local device is updated. When the program in the memory is executed by the processor, all or part of the functions in the above embodiments can be implemented.

[0282] The above examples are used to illustrate the present invention, which are only used to help understand the present invention and are not intended to limit the present invention. Those skilled in the art can make several simple deductions, modifications or substitutions based on the concept of the present invention.

Claims

1. An image processing method for a digital X-ray virtual grid, characterized in that: include: Acquiring initial parameters of each pixel in the imaging image according to the imaging image of the digital X-ray device; The initial parameters include pixel grayscale value, pixel depth value and / or pixel size; performing scattering effect compensation on the initial parameters of each pixel of the imaging image, and using the imaging image after the scattering effect compensation as the output image of the digital X-ray device; wherein the compensation value for the scattering effect compensation of the initial parameters of each pixel is obtained based on a preset ray interference compensation mathematical model, and the ray interference compensation mathematical model is obtained by evaluating and modeling the scattering signal of the digital X-ray device; The mathematical model of ray interference compensation is: ； Where x and y are the pixel coordinates in the image, and ⊗ is the convolution operator; , is the position vector of a single pixel in the projected image, x c and y c is the projection coordinate; , for The scattered ray function centered at describes the coordinates; is the primary ray function I s , used to represent the parameter value of each pixel generated by the primary ray signal; is the scattered ray function I p , used to represent the parameter value of each pixel generated by the scattered ray signal; is the scattering kernel function; The initial parameter value of each pixel in the imaged image is the sum of the parameter values ​​generated for the same pixel by the primary ray signal and the scattered ray signal; According to the scattered ray function I s and primary ray function I p The relationship between the ray interference compensation mathematical model is solved, and the scattered ray function I s and primary ray function I p The relationship function is the product of the scattered signal intensity function and the signal distribution function, and its expression is: ； in, is the scattered ray function I s and primary ray function I p The relationship function, is the scattered signal intensity function of the small-angle scattered signal, is the scattered signal intensity function of the large-angle scattered signal; The scattered signal intensity function expression for small-angle scattered signals is: ； in, is the preset strength coefficient obtained by fitting; is the scattering effect threshold cutoff function; is the scattered signal intensity function expression of the small-angle scattered signal; is the air projection value, and ,in, is the cutoff threshold parameter, The value of is between 0.9 and 1; The distribution function form of the scattered signal intensity function for large-angle scattered signals is: ； in, is the peak width coefficient obtained by fitting, d is the distribution range used to control the scattered ray signal; is the distribution function of the scattered signal intensity function of the large-angle scattered signal; is the two-dimensional Zernike orthogonal basis function The linear superposition polynomial of are the polynomial coefficients; N is a natural number, indicating the number of selected Zernike basis functions.

2. The image processing method according to claim 1, wherein: Also includes: The measured ray function I and the scattered ray function I s And the primary ray function I p Expand into vector form respectively , , , then the scattered ray function I s and primary ray function I p The following relationship exists: ； Among them, S jk is the element of the scattering matrix S to be obtained, which represents all For a certain pixel j The composition of K is the maximum number of columns in the scattering matrix.

3. The image processing method according to claim 2, wherein: Also includes: The meta-heuristic iterative optimization algorithm of generalized pattern search is applied to calculate the corresponding scattered ray function I for the image data obtained from at least two test phantom experiments. s The error function is obtained, and the ray interference compensation mathematical model is optimized by a multi-objective function optimization method based on the wolf pack model.

4. A computer-readable storage medium, characterized in that The medium stores a program, which can be executed by a processor to implement the image processing method according to any one of claims 1 to 3.

5. An image processing device for a digital X-ray virtual grid, characterized in that: The image processing device is configured to post-process an image acquired by a digital X-ray device using the image processing method according to any one of claims 1 to 3, the image processing device comprising: An initial parameter acquisition unit is used to acquire initial parameters of each pixel in the imaging image according to the imaging image of the digital X-ray device; the initial parameters include pixel grayscale value, pixel depth value and / or pixel size; a parameter compensation unit, configured to perform scattering effect compensation on the initial parameters of each pixel of the imaging image; wherein a compensation value for the scattering effect compensation of the initial parameters of each pixel is obtained based on a preset ray interference compensation mathematical model, and the ray interference compensation mathematical model is obtained by evaluating and modeling the scattering signal of the digital X-ray device; An image output unit is configured to use the imaging image after the scattering effect compensation as an output image of the digital X-ray device.

Images