CN117237205A - Image processing method and device for digital X-ray virtual grid - Google Patents

Abstract

The application discloses an image processing method and device for a digital X-ray virtual grid, which comprises the steps of firstly, acquiring initial parameters of each pixel point in an imaging image of a digital X-ray device according to the imaging image; then, carrying out scattering effect compensation on initial parameters of each pixel point of the imaging image; and finally, taking the imaging image subjected to scattering effect compensation as an output image of the digital X-ray device. The compensation value for performing scattering effect compensation on the initial parameter of each pixel point is obtained according to a ray interference compensation mathematical model. The radiation interference compensation mathematical model is obtained by evaluating and modeling the scattered signals of the digital X-ray device, so that the influence of the scattered rays on the digital X-ray image can be eliminated to the greatest extent.

Description

Image processing method and device for digital X-ray virtual grid

Technical Field

The application relates to the technical field of digital X-ray imaging, in particular to an image processing method and device for a digital X-ray virtual grid.

Background

Digital Radiography (DR) is an advanced radiography technique formed by a combination of computerized digital image processing techniques and X-ray radiation techniques. For X-ray images, as the object is taken thicker and the density is higher, the higher the energy of the radiation needs to penetrate the object to create an effective image contrast. However, when the radiation energy is higher and the density of the object to be photographed is also higher, the radiation is more likely to scatter. The X-ray scattered on the detector is overlapped with the original straight ray to form an image together, so that the gray level of the image is increased in a nonlinear manner, local fog is formed, the true contrast of the image is affected, and therefore, how to eliminate the scattering effect of the image is important.

At present, a physical grid method is mainly adopted for suppressing X-scattered rays, and a metal plate made of special materials is adopted in the method, wherein the radiation attenuation materials with special structures can filter out scattered rays which change in the direction after penetrating through a human body, and the scattered rays are prevented from reaching a detector to further image, so that the purpose of eliminating the scattered rays is achieved. The physical grid is used as a selecting accessory and needs to be matched with a grid pattern removing image processing program of corresponding specification parameters; the physical grid can inhibit most of scattered rays, but can absorb part of original straight rays at the same time, so that shooting dosage must be increased in order to achieve the same image gray level, and the radiation level of a patient is increased; on the other hand, for some shooting parts with little scattering influence, the grids are not needed to be used, so that the grids are required to be frequently assembled and disassembled, and therefore, in the development of X-ray equipment, convenience in use and convenience in maintenance are often required to be considered, a relevant sensor needs to be added to sense the loading state of the grids, so that doctors can conveniently observe the use state of the grids to avoid misoperation, and the factors increase the development and maintenance cost of the equipment; meanwhile, the physical grid is used, the specification parameters such as focal length, front and back sides, centering line and the like are required to be paid attention to, and the operation intensity of doctors is increased. Compared with a physical grid, the digital grid can eliminate scattering influence of images more conveniently. There are also related image processing methods of digital grids, but most of them are processed by estimating related information only from a photographed image, and cannot eliminate the influence of scattered rays to the maximum extent.

Disclosure of Invention

The application aims to solve the technical problem of eliminating the influence of scattered rays on a digital X-ray image to the greatest extent.

According to a first aspect, in an embodiment there is provided an image processing method for a digital X-ray virtual grid, comprising:

acquiring initial parameters of each pixel point in an imaging image of a digital X-ray device according to the imaging image; the initial parameters include pixel gray values, pixel depth values, and/or pixel sizes;

performing scattering effect compensation on the initial parameters of each pixel point of the imaging image, and taking the imaging image subjected to the scattering effect compensation as an output image of the digital X-ray device; the compensation value for performing scattering effect compensation on the initial parameter of each pixel point is obtained according to a preset ray interference compensation mathematical model, and the ray interference compensation mathematical model is obtained by evaluating and modeling a scattering signal of the digital X-ray device.

According to a second aspect, an embodiment provides a computer readable storage medium having stored thereon a program executable by a processor to implement the method according to the first aspect.

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

an initial parameter obtaining unit, configured to obtain initial parameters of each pixel point in an imaging image of a digital X-ray device according to the imaging image; the initial parameters include pixel gray values, pixel depth values, and/or pixel sizes;

the parameter compensation unit is used for carrying out scattering effect compensation on the initial parameters of each pixel point of the imaging image; the compensation value for performing scattering effect compensation on the initial parameter of each pixel point is obtained according to a preset ray interference compensation mathematical model, and the ray interference compensation mathematical model is obtained by evaluating and modeling a scattering signal of the digital X-ray device;

and the image output unit is used for taking the imaging image subjected to the scattering effect compensation as an output image of the digital X-ray device.

According to the image processing method of the embodiment, the applied ray interference compensation mathematical model is obtained by evaluating and modeling the scattering signal of the digital X-ray device, so that the influence of the scattered rays on the digital X-ray image can be eliminated to the greatest extent.

Drawings

FIG. 1 is a flow chart of an image processing method in an embodiment;

FIG. 2 is a schematic diagram of core scattering kernel stacking principle in one embodiment;

FIG. 3 is a schematic diagram of a method for simulating pencil beam X-ray scattering distribution using a Monte Carlo simulation method;

FIG. 4 is a diagram of asymmetric mobile scattering kernel function parameters in one embodiment;

FIG. 5 is a schematic diagram of an experimental measurement apparatus for automatically optimizing scattering model parameters in one embodiment;

FIG. 6 is a diagram of a deconvolution algorithm flow architecture in one embodiment;

FIG. 7 is a flowchart of a multi-scale residual image BLR iterative algorithm based on Gaussian pyramid in one embodiment;

FIG. 8 is a graph showing the final scattering suppression effect of a shifted scattering kernel versus a conventional scattering kernel function in one embodiment;

FIG. 9 is a graph showing the effectiveness of the LR based deconvolution algorithm versus the conventional deconvolution algorithm in one embodiment;

fig. 10 is a block diagram showing the structure of an image processing apparatus in another embodiment.

Detailed Description

The application will be described in further detail below with reference to the drawings by means of specific embodiments. Wherein like elements in different embodiments are numbered alike in association. In the following embodiments, numerous specific details are set forth in order to provide a better understanding of the present application. However, one skilled in the art will readily recognize that some of the features may be omitted, or replaced by other elements, materials, or methods in different situations. In some instances, related operations of the present application have not been shown or described in the specification in order to avoid obscuring the core portions of the present application, and may be unnecessary to persons skilled in the art from a detailed description of the related operations, which may be presented in the description and general knowledge of one skilled in the art.

Furthermore, the described features, operations, or characteristics of the description may be combined in any suitable manner in various embodiments. Also, various steps or acts in the method descriptions may be interchanged or modified in a manner apparent to those of ordinary skill in the art. Thus, the various orders in the description and drawings are for clarity of description of only certain embodiments, and are not meant to be required orders unless otherwise indicated.

The numbering of the components itself, e.g. "first", "second", etc., is used herein merely to distinguish between the described objects and does not have any sequential or technical meaning. The term "coupled" as used herein includes both direct and indirect coupling (coupling), unless otherwise indicated.

Methods for realizing scattered radiation suppression at the design level of the traditional imaging system include methods such as narrow beam scanning, air gap increase, grid hardware configuration and the like. These methods have some effect in certain system designs, but none of them completely eliminates scattered radiation interference. For example, in a chest-abdomen scanning modality, the beam width is difficult to reduce. The body part has larger volume, so that the scattered ray interference effect is most obvious. Increasing the air distance is a routine skill used clinically, but digital X-ray devices are difficult to apply in special use scenarios (e.g., mobile DR systems applied to patient-side monitoring are not able to employ this routine skill). Grid is currently a well-established technology in widespread use, but its requirement for mounting locations and specific board supply distances (source to detector distance) makes it inconvenient to use in portable DR imaging systems.

In order to solve the limitation that the digital X-ray device is only suitable for specific situations, in the embodiment of the application, an image processing method is provided, the image post-processing is only carried out in a scattering effect compensation mode, the effect approximately equivalent to a real hardware grid can be achieved under the condition that the imaging dosage is reduced by 1/3, and the digital X-ray device is further enabled to relieve the application environment limitation and has the characteristics of convenience, universality and effectiveness.

Embodiment one:

referring to fig. 1, a flowchart of an image processing method in an embodiment is shown, where the image processing method is used for post-processing an image acquired by a digital X-ray device, and specifically includes:

step 101, an original image is acquired.

An initial parameter of each pixel point in an imaging image of a digital X-ray device is obtained according to the imaging image, wherein the initial parameter comprises a pixel gray value, a pixel depth value and/or a pixel size.

Step 102, scattering effect compensation.

And carrying out scattering effect compensation on the initial parameters of each pixel point of the imaging image, wherein a compensation value for carrying out scattering effect compensation on the initial parameters of each pixel point is obtained according to a preset ray interference compensation mathematical model, and the ray interference compensation mathematical model is obtained by evaluating and modeling the scattering signals of the digital X-ray device.

And step 103, outputting a result image.

And taking the imaging image subjected to scattering effect compensation as an output image of the digital X-ray device.

In one embodiment, the ray disturbance compensation mathematical model is:

；

；

；

where x and y are the pixel coordinates in the image,for projection of an image a single pixel position vector, x _c And y _c For projecting coordinates +.>To->Describing coordinates for a scattered ray function of the center; i _p A parameter value representing each pixel point generated from the primary ray signal as a primary ray function; i _s A parameter value representing each pixel point generated by the scattered ray signal as a scattered ray function; the initial parameter value of each pixel point in the imaging image is the sum of the parameter values generated for the same pixel point due to the primary ray signal and the scattered ray signal.

In one embodiment, the image processing method further includes:

according to the scattered ray function I _s And primary ray function I _p Solving the ray interference compensation mathematical model, and scattering ray function I _s And primary ray function I _p Is the product of the scattered signal intensity function and the signal distribution function, and is expressed as:

；

where h is the scattered ray function I _s And primary ray function I _p Is a function of the relation of (c),scattering signal intensity function for small angle scattering signal, < >>As a function of scattered signal strength for large angle scattered signals.

The scatter signal intensity function for a small angle scatter signal is expressed as:

；

wherein alpha is a preset fitting obtained intensity coefficient;is a scattering effect threshold cut-off function; i ₀ Is the air projection value, and ∈>；

The distribution function of the scattering signal intensity function for the large angle scattering signal is as follows:

；

wherein, beta is a peak width coefficient obtained by fitting, and d is a distribution range for controlling scattered ray signals;is a two-dimensional Zernike orthogonal basis function +.>Linear superposition polynomial of c _i Is a polynomial coefficient.

To measure ray function I and scatter raysFunction I _s Primary ray function I _p Respectively spread into vector form，/>，/>Scattered ray function I _s And primary ray function I _p There is a relationship between:

；

wherein S is _jk As elements of the scattering matrix S to be obtained, which represents allFor a certain pixel j +.>Is a component part of (a).

Each column of elements in the scattering matrix S is not the same shiftBased on this measurement radiation function I +.>The expression of (2) is:

；

wherein, the imaging matrix h=diag {1, …,1} +s, that is, the scattering matrix is a block diagonal matrix with one element being 1, all elements are non-negative, and the values of a large number of elements far from the main diagonal gradually decrease and tend to be zero.

In one embodiment, a meta-heuristic iterative optimization algorithm is applied to the generalized pattern search to test at least two test motifsThe acquired image data respectively calculate the corresponding scattered ray functions I _s And optimizing the ray interference compensation mathematical model by a multi-objective function optimization method based on a wolf-swarm model.

The following theory describes the modeling process of the ray interference compensation mathematical model (for convenience of description, the ray interference compensation mathematical model is collectively referred to as a scattering kernel model in the following description):

referring to FIG. 2, a schematic diagram of core scattering nuclear superposition principle in an embodiment is shown, wherein the observation signal (i.e. detector measurement signal) obtained by a general X-ray imaging process including primary rays, scattered rays and an irradiated object can be regarded as primary rays I _p Signal and scattered radiation I _s Superposition of signals. For each pixel point in the imaged image there is:

I(x,y)= I _s (x,y) + I _p (x,y)；

wherein I is a detector measurement signal of the digital X imaging device, I _s For initial ray signal, I _p For scattered radiation signals, it is understood that the detector imaging signal of the digital X-ray imaging device is a superposition of the primary radiation signal and the scattered radiation signal, which are not independent of each other but have a correlation. The focus in the mathematical modeling process of the X-ray imaging process is on a method for estimating scattered signals. In the embodiment of the application, a scattering kernel superposition model is used, and a scattering signal is understood as a superposition result of a plurality of primary ray signals after the convolution of a kernel function. The core function is considered to be in a form of central symmetry and linear movement in the traditional scattering core model. This assumption stems from the method used for scatter correction in CT imaging, which applies to imaging system geometries that use arcuate detectors. But this assumption fails for DR systems using flat panel detectors, especially where the cone angle of the radiation is large. This can be observed in Monte Carlo simulation, please refer to FIG. 3, which is a schematic diagram of a method for simulating pencil beam X-ray scattering distribution by Monte Carlo simulation, in which the obtained scattering nuclear signal distribution is deviated from the primary main ray position after the large cone angle rays pass through the water phantom (water cube) with the depth of 20cm, and is not centrosymmetric any more Form of the application. The degree of deviation is different when the depth of the water phantom and the ray cone angle are changed, so that the scattering nuclear signal is not in a linear motion-invariant form. Based on the above observations, the present application proposes a generalized version of the asymmetric shift-and-scatter kernel function. The non-centrosymmetric component is added on the basis of the centrosymmetric component. And the component is a spatially shifted version (i.e., the corresponding component is different for each pixel of the image). The non-centrosymmetric component is expressed by a zernike basis function polynomial, implementing a complex functional representation of a small number of parameters. The modeling method is more approximate to the actual physical process, so that the scattering inhibition effect is more uniform. Wherein the imaging matrix contains all the information of the primary rays and the scattered rays, and establishes a direct connection between the observed image (actual imaging) and the ideal non-scattered image (scattering-suppressed image). The calculation method based on the imaging matrix form can solve the problem that Fast Fourier Transform (FFT) is not utilized to accelerate shift convolution kernel operation. Because the matrix is a semi-sparse block diagonal matrix, the convolution process of the shift-change convolution kernel can use a parallel acceleration algorithm of conventional matrix operation, such as GEMM, spGEMM, spAMM, and the like. The imaging matrix is the core calculation method of the application, which enables the modeling of the linear-motion scattering kernel function to be implemented.

The embodiment of the application also discloses a method for automatically optimizing the model parameters, and compared with a symmetric shift-invariant scattering kernel, the number of parameters in an asymmetric shift-variant scattering kernel function is multiplied. Because of the problems of the rapid increase of the number of measurement experiments, the improvement of the precision requirement of an experimental instrument and the like, the method for directly measuring the asymmetric shift scattering kernel distribution becomes complicated and difficult. In one embodiment of the application, a method combining Monte Carlo simulation and a small amount of low-cost imaging experiments is provided, and parameter adjustment is automatically optimized. The parameter fitting method based on Monte Carlo simulation is to simulate the physical process of penetrating the water mold body by using the existing Monte Carlo simulation software. The energy deposition signal distribution of the secondary scattered rays generated during the recording process in the detector. The source spectrum and detector response functions are calibrated for a particular imaging system. So that the simulation results are closer to the imaging system. And performing preliminary fitting on the parameters of the mobile scattering cores in the model according to the simulation data, and automatically completing establishment of initial values of the parameters.

In an embodiment of the present application, an imaging experiment is performed in the target system after the initial value of the parameter is established based on the parameter optimization strategy of the experimental measurement data. The scatter signal distribution in phantom X-ray imaging is measured separately by a method of continuous scanning with a sliding lead block ray blocker array. The primary radiation signals are measured separately by means of a sliding pencil beam limiter scanning method. And further optimizing the parameters of the shift scattering kernel by using the two data so as to achieve the optimal effect under a specific system.

In addition, the scatter kernel superposition model performs forward modeling from non-scattered primary ray signals to actual imaging with scatter. The essence is convolution process, which belongs to image degradation phenomenon. However, the conversion of actual imaging to a non-scattering image is a problem of the inverse. The solving method is deconvolution algorithm, which belongs to the image restoration problem. In an embodiment of the present application, a stable and effective deconvolution algorithm is specifically disclosed, that is, a deconvolution solution method for providing a shift-transform convolution kernel. Most of deconvolution solving algorithms given in the existing schemes are semi-empirical methods, and the convergence and uniqueness of the convergence result are not strictly deduced. In fact, whether it converges or not and the convergence result are both related to the manual selection of free parameters therein and to the SPR of the imaged object, exhibit instability and limited scatter suppression. In addition, the output result is often accompanied by noise amplification. The application provides a stable and strictly convergent deconvolution algorithm of a shift convolution kernel based on an LR deconvolution algorithm. Which suppresses deconvolution output image noise while preserving edge detail using bilateral regularization terms. In addition, a Gaussian pyramid structure is introduced, deconvolution operation is carried out on residual images in different frequency domains, so that the algorithm convergence speed is effectively improved, and the streak artifacts are restrained.

The application also provides a parallel acceleration scheme of the deconvolution algorithm, wherein the calculation process of the deconvolution algorithm with the Gaussian pyramid structure is sequentially executed layer by layer, and the calculation time is long under the condition of giving high-quality deconvolution images. In one embodiment of the present application, another parallel acceleration method is provided in which the layers are sequentially executed. So that the calculation time length does not linearly increase with the deepening of the pyramid layers.

As shown in fig. 2, it is assumed that the scattering kernel model corresponding to the scattered ray signal is obtained by convolving the primary ray signal with a kernel function associated with the digital X-ray imaging apparatus. Further simplified, the scattering kernel model is in the form of linear spatial shift invariance (specially shift invariant) under nearly parallel beam geometry imaging conditions (i.e., smaller cone angle), initial ray signal I _s Can be expressed as:

；

initial ray signal I _s Is identified in simplified form as:

；

；

；

wherein,for projection image single pixel position vector, but +.>To->Coordinates are described for the central scattering kernel function, and (x, y) is the pixel coordinates in the image. According to I _s And I _p Is used for solving the relation of the scattering kernel function I _p The scatter signal strength function and the signal distribution function product are expressed as:

；

By attenuation analysis of the small angle scattering signal, the following signal strength function form can be obtained:

setting I to be a measurement signal, I _s Is a scattered signal, I _p Is the primary signal, and is obtained as I _p (real image after filtering out scattered signals), h is I _s And I _p The relation between the expression is:

；

where α is the intensity coefficient resulting from the fitting, anIs a scattering effect threshold cutoff function.

When the gray level of the projection image is close to that of the space exposure image, the scattering effect is not considered, namely the condition needs to be satisfied:

；/>；

wherein I is ₀ Is the air projection value (signal received by the detector without any object). E is a cutoff threshold parameter, which can be adjusted between 0.9 and 1 as the case may be.

Referring to fig. 4, a schematic diagram of parameters of an asymmetric shift-change scattering kernel function in an embodiment is shown, where a conventional spatial shift-invariant signal distribution function is generally in a central symmetric form. However, considering the situation that the large angle scattering signal deviates to the central ray under the large cone angle geometry, the signal distribution no longer has a central symmetry distribution form, so the distribution function form of the scattering kernel function improved for this problem in an embodiment of the application is as follows:

；

The first term in the expression of the distribution function form is a symmetric term, and represents the corresponding basic part of the ideal central ray. Beta is the peak width coefficient obtained by fitting, and d is the distribution range for controlling the scattering signal. And the second term is used to describe the asymmetric effect introduced by the off-center ray case, expressed as:

；

wherein,is essentially a two-dimensional Zernike (Zernike) orthogonal basis function +.>Linear superposition polynomial of c _i Is a polynomial coefficient.

The complete zernike basis function system may describe any continuously varying two-dimensional distribution function. The Meng Ka simulation observation shows that the large cone angle corresponds to the asymmetric trend in the scattering core, and the surrounding central rays are distributed in a central symmetry mode. In practical use, it is therefore sufficient to use a selected n=5-term basis function, the expression of which refers to the following table one, which is a 5-term zernike basis function expression applied in one embodiment:

list one

Considering that the distribution function at different positions has rotation, the rotation angle is introducedIs described. And the zernike basis function domain is only a unit circle range (i.e., 0<ρ<1) The range normalization factor d=3β is added.

In the spatially-shifted scattering kernel function, the original ray signals of different cone angle geometries all have different non-centrosymmetric scattering signal components, which vary with different pixel coordinates of the projection image. The zernike polynomial coefficients in the scattering kernel function are a function of the pixel coordinates. In an embodiment of the present application, it is assumed that the trend has continuity, and further, a binary polynomial (bivariate polynomial) is used to model the trend of the zernike polynomial coefficient along with the change of the ray cone angle, namely:

；

In one embodiment, setting d=2, one can obtain:

；

as aboveThe expressions are all desirably continuous forms, and matrix linear operation expressions in discrete forms are derived below.

Any linear operation can be converted into a form of matrix multiplication. The spatially invariant convolution kernel can be reconstructed as a block diagonal matrix (Toeplitz matrix) and is a special cyclic matrix. The convolution matrix corresponding to the convolution kernel of the space shift 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 circular matrix. Will measure ray function I, scatter ray function I _s Primary ray function I _p Respectively spread into vector form，/>，/>Scattered ray function I _s And primary ray function I _p There is a relationship between:

；

wherein S is _jk As elements of the scattering matrix S to be obtained, which represent all the initial signalsFor a scattering signal in a pixel j>Is a component part of (a). Each column element in the scattering matrix S is not a different shift scattering kernel +.>Is based on the discrete value expansion of the measurement signal +.>The expression of (2) is:

；

wherein, the imaging matrix (without filtering scattered ray signals) h=diag {1, …,1} +s, i.e. the scattering matrix plus a diagonal matrix with one element being 1, the matrix is a block diagonal matrix, all elements are non-negative, and a large number of element values far from the main diagonal gradually decrease and tend to be zero. From the perspective of numerical computation, the imaging matrix is a semi-sparse diagonal matrix, and parallel acceleration methods for the semi-sparse diagonal matrix operation can be applied here. Convolution operations on spatially invariant convolution kernels may be accelerated using a fast fourier transform, but the shifted convolution kernels cannot use this strategy. There is a limit to the parallel convolution acceleration method using the conventional loop nesting manner. In particular, in the acceleration process using a graphics computing unit (GPU), the shared memory and the thread registers are limited, and the convolution kernel sampling limit is limited, so that the convolution calculation accuracy is finally affected. (this is to say that acceleration is limited in non-matrix form). Pure linear operation is performed by means of an imaging matrix, so that a conventional matrix-operation parallel acceleration algorithm (GEMM, spGEMM) can be applied.

In one embodiment, to accurately fit the asymmetric components in the large cone angle scattering kernel, the influence of the larger value of the small angle scattering signal on the fitting accuracy is avoided, and a step fitting method is adopted. The thought is to calculate Meng Ka the residual errors of the simulation data and the main terms of the fitting function step by step, and the residual errors are used for finer correction terms in the subsequent fitting model function so as to realize automatic adjustment of the parameters of the scattering kernel model.

Since the scattering kernel function is related to the imaging system geometry, X-ray energy, metal filtering, and the depth of penetration of the radiation into the imaged object. In the prior art, a method for simulating scattered ray signal distribution by a Monte Carlo method is adopted in the embodiment of the application. Although the Meng Ka simulation results constructed due to modeling bias deviate somewhat from the actual imaging system measurement results. However, considering the accuracy of the endogenous physical processes in the Monte Carlo simulation mechanism, the obtained scattering signal distribution is a good starting point for the application of the real system scattering kernel modeling. To ensure that the Monte Carlo simulation results reflect to some extent the characteristics of a particular imaging system, the critical parameters (including the source spectrum and detector response parameters) need to be calibrated. In an embodiment of the application, the energy deposition response function curve measurement of the flat panel detector is directly measured by utilizing the single-energy synchrotron radiation X-rays with different energy levels and energy resolutions, and the specific method comprises the following steps:

Step 201: the size of the field of the synchrotron radiation light source is adjusted to be within 10cm x and 10cm, and the crystal dichroic mirror is adjusted to enable the energy resolution to reach the optimal state.

Step 202: and measuring the dark field of the detector, and finishing offset correction. The flat panel detector is controlled to be 90 degrees and opposite to X-rays, and the detector is continuously moved to measure the bright field image, so that gain correction is completed.

Step 203: the synchrotron radiation beam energy (30 keV-120 keV) and the angle of incidence (0 °, 30 °, 45 °, 60 °, 70 °, 80 °) are varied and the beam flux is calculated from the beam control ionization chamber data.

Step 204: and processing an effective area of the image acquired by the detector, and converting the effective area to obtain an effective output value.

Step 205: and carrying out normalization processing on all the data to obtain a deposition energy response curve of the detector under different energies and different incidence angles.

In one embodiment, a method for calibrating a radiation source energy spectrum includes:

step 301: and calculating the initial value of the energy spectrum under the corresponding energy level and metal filtering state through a semi-empirical model.

Step 302: and measuring the attenuation proportion of the X-ray source after penetrating through aluminum ladders with different thicknesses and acrylic plates.

Step 303: and obtaining the linear attenuation coefficient of the metal aluminum and acrylic plate.

Step 304: and calculating an energy spectrum optimal estimated value through an initial energy spectrum weighted truncated singular value decomposition (truncated singular value decomposition, TSVD) algorithm according to the attenuation measurement data of the aluminum ladder and the energy spectrum initial value obtained in the step 301. And converting the equivalent energy spectrum according to the response function of the detector.

Step 305: and calculating the attenuation ratio of different acrylic plates according to the equivalent energy spectrum, and calculating the deviation between the attenuation ratio and the measured value in the step 302.

Step 306: if the deviation in step 305 is greater than the set threshold, the update parameters are returned to step 304; and if the deviation in the step 5 is smaller than the set threshold value, ending.

In one embodiment, the Monte Carlo simulation scattering distribution method uses the existing Monte Carlo software to simulate the scattering distribution formed by penetrating pencil beam rays through water mold bodies with different thicknesses. The virtual detector end records photon fluxes of different energies and different incidence angles. Recording primary particles and secondary particles formed by Compton scattering and Rayleigh scattering, respectively, denoted asAnd->. And converting the incident photon flux into a detector output signal according to the obtained detector response function. The thickness range of the water mold body is 0 cm-50 cm, and the simulation is carried out at intervals of 1 cm.The cone angle ranges from 0 degrees to 15 degrees, and simulation is performed at intervals of 1 degree.

The following is a fitting method for parameters of a scattering kernel function model, and the specific fitting method comprises the following steps:

step 401: data fitting is performed for the linear motion invariant part parameters in the scattering kernel function model.

Defining an optimized loss functionThe method comprises the following steps:

；

and (3) completing the nonlinear multi-element function fitting by using a Simplex optimization algorithm (Simplex downlink) to obtain an intensity coefficient alpha and a peak width coefficient beta.

Step 402: and carrying out data fitting on the non-central symmetry part parameters in the linear shift scattering kernel function. This step uses the least squares idea for coefficient fitting. The unit circles of the continuous zernike function are orthogonal to each other in the definition domain, but the discrete zernike function is not. The basis vector Z is used here by the Schmidt (Gram-Schmidt orthogonalization) method _k (k=1, …, N) construction of mutually orthogonal discrete zernike basis vectors V _k (k=1, …, N) of the formula:

；

where G is a transform matrix in the form of a lower triangle (lower triangular matrix).

Orthogonal basis vector V _k The calculation method comprises the following steps:

；

transform matrix G element G _k,j The calculation method comprises the following steps:

；

where j=1, …, k-1, and when j=k, g _k,j =1。

The residual data may be represented as a linear superposition of discrete zernike basis vectors, i.e.:

；

the following tasks are used to calculate the superposition coefficient b _i (i=1, …, N), specifically including:

constructing a residual data base line:

。

then, after the mean baseline interference is eliminated, the following relationship holds:

。

therefore, the loss function is optimized for this stepThe definition is as follows:

；

calculation of optimized loss function using least squares principleWith respect to the optimization variable b _i The gradient vector of (2) is:

，

can solve coefficient b _i The expression is:

；

the original zernike polynomial coefficients were then calculated as:

。

step 403: fitting to obtain the variation trend of the mobile scattering kernel parameter along with different pixel coordinates. In the linear shift scattering kernel model, the zernike function coefficient c _i (x _c ,y _c ) As a function of projected image coordinates. The zernike polynomial coefficients corresponding to the different non-centrosymmetric signals in the discrete coordinates have been obtained in step 402. Binary polynomial coefficients describing the trend of the coefficient change are fitted based on this fitting.

Using an optimised loss functionThe definition is as follows:

；

according to the least squares principle, byObtaining N groups of linear equations:

；/>

wherein,。

(Vector)the method comprises the steps of carrying out a first treatment on the surface of the And vector quantityThe method comprises the steps of carrying out a first treatment on the surface of the The method directly comprises the following steps of:

obtaining a matrix pseudo-inverse, and finally obtaining an optimization coefficient as follows:

。

the application also discloses a parameter optimization method of the scattering kernel function model, and the parameter optimization method optimizes model parameters according to the real measurement data of a specific imaging system, so that the scattering inhibition uniformity of the imaging image in the system is improved.

Referring to fig. 5, a schematic diagram of an experimental measurement apparatus for automatically optimizing parameters of a scattering model according to an embodiment is shown, in which a dense scattering signal is directly measured by means of a continuously sliding back ray stopper, and the same dense distribution signal of primary rays is directly measured by moving a pencil beam limiter to a plurality of positions. And then, continuously optimizing the scattering kernel function model parameters by combining a self-adaptive grid search algorithm without gradient information. An equivalent experimental method for indirectly measuring the true linear shift scattering core is formed. Meanwhile, generalization performance of the model in real data is enhanced, and stability of a scattering suppression algorithm is improved.

The following intermediate functions are first defined:

；

；

；

；

；

；

in one embodiment, multiple mold bodies are used to move different positions, including: respectively measuring chest mould body, abdomen mould body, upper limb mould body and lower limb mould body, obtaining the corresponding signals of the four mould bodies, respectively calculating the corresponding error functions, and finally, refining a multi-objective optimization problem by uniformly optimizing the four functions, wherein the multi-objective optimization problem is expressed as follows:

；

in one embodiment, a multi-objective function optimization method based on a wolf-swarm model is applied, and the kernel is a meta-heuristic iterative optimization algorithm of generalized mode search. Marking the hierarchyAnd the search dimension is reduced by iteration rotation.

In one embodiment, a Lewye-Richardson (LR) iterative algorithm is applied to solve the deconvolution problem presented in the mathematical modeling above, and certain improvements and innovations are made for solving the specific problem presented in the present application. The following derives the LR algorithm in the context of X-ray imaging, including in particular:

the process of detecting X-ray photons by the flat panel detector is a random process that satisfies poisson distribution. The detector is provided with K pixels in total. It is at the expected value:

；

in which case signals are generated by capturing X-ray photons The expression is:

；

generating a signalPosterior of (2)The probability expression is:

；

taking a maximum estimate (Maximum Likelihood Estimation) of the posterior probability expression, one can rewrite:

；

constant term is omitted without changing the maximum point of likelihood functionThe method can obtain:

。

bringing into imaging modeling expressionsThe following steps are:

；

and calculating the extreme points of the likelihood function by adopting a gradient descent (Gradient Descendent) numerical iteration method. The input image vector and the output image vector in a single step iteration are respectively，/>And the iterative equation is written as:

；

wherein λ is a step size parameter used to adjust the convergence rate.

Then, likelihood function gradient components are calculated separately:

；/>

note thatElements can be written as being transposed of their matrix>Is->The element, likelihood function gradient vector expression is:

；

wherein the method comprises the steps ofIs a unit vector, the size of which is the same as that of the image vector, let +.>The LR algorithm iterative equation can be finally obtained by bringing the gradient descent iterative equation into the gradient descent iterative equation:

；

wherein onlyThe matrix multiplication operation is realized, and other multiplication and division operations are element operations. For clarity, the component forms are written below:

；

for convenience of writing component subscripts, whereAnd +.>Respectively represent->，/>Is a component of (a).

The LR algorithm iterative equation is typically used iteratively, and the number of iterations required for convergence is typically 10-20.

The error function for judging convergence is typically:

；

when (when)Less than the preset convergence threshold->The iteration is ended.

The LR algorithm is a product operation, and naturally ensures non-negativity of a convergence result, and the convergence result is stable. However, when processing larger-scale images, the overall working efficiency of the algorithm is lower due to the longer convergence iteration times. And because the whole frequency spectrum of the signals in the X-ray image is distributed widely, the convergence rates of the high-frequency signals and the low-frequency signals are different. In the region where the edges are noticeable, gibbs Phenomenon (Gibbs Phenomenon) is liable to occur such that noticeable streak artifacts appear in the deconvoluted image. In addition, the algorithm convergence effect is also related to the initial image effect because the influence of Gaussian noise is not considered in the derivation process. The invention is directed to various improvements to these problems, and will be described in detail below.

The improvement strategy for LR algorithms is generally focused on modification of the maximum likelihood function. The poisson process-based maximum likelihood function itself does not concern the spatial relationship between the image signals. No prior assumption is made of the multi-scale feature information and the noise correlation information. By adding the regular term in the likelihood function, the algorithm convergence junction effect can be directed to a priori assumptions preset in the regular term. Typical regularization terms are Tikhonov-Miller (TM) regularization terms (i.e., gradient 2-norm regularization terms), total Variation (TV) regularization terms (i.e., gradient 1-norm regularization terms), which constrain local signal variations by a simple method that limits the overall gradient of the image to achieve noise reduction. However, with the consequent problem of more local smoothing of the image. There should be a preserving effect of adaptive edge information in the regularization term. A Bilateral regularization LR algorithm (BLR) was introduced below. The essence of the method is derived from the idea of Bilateral filtering (Bilatial Filter), and the Bilateral regular expression is as follows:

；

Wherein Ω is a neighborhood of pixel i, the size of which controls the range of action of the regularization term, functionIs an action range control item.

While the functionTo smooth the control terms, a certain penalty effect is given to signal fluctuation in the neighborhood. Since the regularization term needs to be minimized, the likelihood function form adjusts to:

；

the deconvolution problem is solved and converted into a minimum value solution optimization problem, and the minimum value solution optimization formula is as follows:

；

repeated gradient descent method for finding loss functionIs the minimum point of (2):

；

calculating a loss functionGradient vector:

；

and (3) making:

；

the iteration equation with gradient descent can finally obtain the iteration equation of the BLR algorithm as follows:

；

the bilateral regular function gradient can be independently calculated according to the following formula, and specifically:

；

wherein,is a function of the deconvolution input image pixel i neighborhood pixel j, whose jth component is transformed by:

；

whileIs an image shift operator which operates to shift the image +.>According to vector->Translation is performed. Vector->The start point is the i pixel position coordinate and the end point is the j pixel position coordinate.

The convergence result of the BLR algorithm can retain more detail than the TV or TM regularization term constrained LR algorithm. The above-mentioned problem of different convergence rates of the high and low frequency signals in the LR algorithm is the same in the BLR algorithm. In order to accelerate the overall convergence speed of the algorithm, the invention adopts a strategy of completing BLR deconvolution operation in parallel by constructing a multi-scale pyramid hierarchical layer. In order to adapt to the simultaneous completion of deconvolution iterative operation in different scale image information, the invention gives up to use the original image for calculation, but uses the residual image instead. In the iteration process of the residual image, the detail low-power signal to be recovered is not interfered by the strong gradient signal of the original image edge, so that the streak artifact is restrained.

The following is a detailed description.

Referring to FIG. 6, a deconvolution algorithm flow structure diagram is shown in an embodiment, in which an N-layer Gao Sijin sub-tower is constructed after performing N-1 times 2 times downsampling on an original X-ray image. And constructing an N-layer Gaussian pyramid by using a Gaussian filtering and difference solving method.

In one embodiment, n=3 is taken as an example for simplicity of description. The downsampled input images are respectively:

；

；

；

the residual image for the corresponding hierarchy is calculated as follows:

then, respectively solving a deconvolution image of each layer of residual images by using a BLR algorithm as follows:

；

the final output complete deconvolution image expression is:

；

referring to fig. 7, a flowchart of an iterative algorithm of a multi-scale residual image BLR based on gaussian pyramid in one embodiment is shown. The BLR algorithm may be further optimized in the framework of a gaussian pyramid. In the multi-level parallel BLR calculation process described above, the inter-layer information is not fully mined, so the algorithm flow is further modified for this. The image which is up-sampled after deconvolution of the upper layer is used as a guide image, and in order to fully mine the edge information of the image, the BLR algorithm is modified into a Joint Bilateral (JBLR) LR algorithm by referring to the idea of an image guide noise reduction algorithm.

Modifying a regular term in a minimum solution optimization formula into a guide function regular term, wherein the expression is as follows:

；/>

the smoothing control items are as follows:

；

then it is necessary to optimizeThe expression is:

；

it is envisioned that solving the iterative equation is:

；

it should be noted that the image is guided in an iterative equationAnd is updated continuously.

After the convergence of the iteration(s),convergence to->。

In addition, the regular function gradient calculation also varies as follows:

；

the method for calculating the residual image is also changed to a certain extent according to the thought of combining Gaussian pyramid interlayer information. The subtracted image is no longer a gaussian smoothed image of the layer, but the upper layer outputs a degraded image after up-sampling of the deconvoluted image. And to eliminate the smoothing effect of bilinear interpolation during upsampling, a simplified BLR algorithm based on gaussian kernel convolution is added after upsampling.

Also taking n=3 as an example, a calculation method of each layer residual image is exemplified.

To distinguish between different types of convolution kernels in the BLR algorithm, the BLR algorithm based on the convolution kernel α is denoted by the symbol BLR (. Cndot.; α), and the JBLR (. Cndot.; α) is understood in the same way. The downsampled images of the input image are respectively:

；

；

；

the residual image of the third layer is:

；

the deconvoluted image of the third layer is:

；

the upsampled image of the third layer is:

；

The residual image of the second layer is:

；/>

the deconvoluted image of the second layer is:

；

the upsampled image of the second layer is:

；

the residual image of the first layer is:

；

the deconvoluted image of the first layer is:

；

the upsampled image of the first layer is:

；

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. The function of this is to compensate for scattering effects in the flat panel detector package and the internal scintillator. The calculation of the deconvolution algorithm part is completed.

The Monte Carlo simulation involved in the embodiment of the application can have alternative methods, which specifically comprise:

1) The method for calibrating the response function of the detector comprises the following steps: the scintillator-related specific parameters are typically core secrets of the flat panel detector manufacturer and are difficult to obtain. However, given knowledge of the detector scintillator chemistry, doping elements, and constituent structural geometry, it is possible to consider as an alternative to simulate the response function of the detector to different energy monoenergetic rays by the Monte Carlo method.

2) The method for calibrating the energy spectrum of the radiation source comprises the following steps: the source energy spectrum measurement may be replaced with a GeTd gamma ray spectrometer. It is commonly used in industrial X-ray tube characterization analysis. If the device is used for clinical X-ray energy spectrum measurement below 120kV, the device has the problems of noise and low-energy interference signals. Further signal correction is required. Under the condition that the anode target material of the X-ray tube, built-in metal filtration and electron scattering (electron scattering) parameters are known, a Monte Carlo method can be used for simulating the process of generating X rays by bombarding the anode target with electron beams. The outgoing photon phase space is recorded for the subsequent simulation process.

The embodiment of the application relates to an experimental instrument, and two special engineering devices, namely a lead block ray blocker array and a pen-shaped beam limiter, are used in the model parameter optimization method. There is also an array of lead stop ray blockers for independent measurement of the scatter signal. Without this array tooling, a separate ray blocker may be used instead, ultimately achieving the same result. But must be multiplied to reduce the efficiency of the experiment. Also, pencil beam limiters are used to independently measure the primary signal, with the aperture orientation thereof being aligned as required to adjust with the cone angle of the ray. The tool without the function can achieve the approximate effect, but under the condition of a large taper angle, the experimental precision is reduced by a larger radial penumbra. In addition, a fan-shaped beam limiter can be used as a substitute to achieve the purpose of increasing experimental efficiency. But this may lead to a decrease in measurement accuracy.

Alternative schemes for modeling and algorithms involved in embodiments of the present application include:

(1) There are various alternatives to the mathematical modeling of the scattering kernel function, and only one of the possible alternatives that is verified to have a practical effect is provided in the embodiments of the present application. The symmetrical components can be replaced by a multi-Gaussian distribution superposition form, and the symmetrical components can be taken as a decomposable convolution kernel to improve the calculation efficiency. The asymmetric components can be replaced by orthographically complete real function bases such as two-dimensional B-spline functions, bessel functions (column harmonics), and the like.

(2) In the optimization algorithm, many options are available for the nonlinear optimization method based on the non-gradient information, which are not described in detail herein.

(3) The imaging matrix-based image degradation problem presented in the embodiments has a variety of general image restoration algorithms, such as Weiner filtering, regularized inverse filtering (inverse filtering), landweber algorithm, tiknonov-Miller algorithm, and FISTA (Fast Iterative soft-threshold) algorithm. The effect varies with the adjustment of the free parameters thereof, and similar results should be obtained finally.

(4) The hardware related to the application of the method specifically comprises the following steps: the deconvolution algorithm may be migrated on any personal computer, workstation or cloud server. The serial version may be accomplished simply with a single threaded operation of a single block Central Processing Unit (CPU). Parallel versions may be accomplished by independent Graphics Processors (GPUs) or multithreaded CPU programming. In addition, the embedded computing processing unit integrated by the DR system can also realize the algorithm proposed by the invention, and possible schemes include Raspberry Pi, field programmable gate array (Field-programmable gate array, FPGA) and embedded system Jetson of Nvidia company carrying independent GPU units.

Referring to fig. 8, a final scattering-suppressing effect comparison of the function of the transition scattering kernel and the conventional scattering kernel is shown as a non-scattering-suppressing human phantom DR image, a scattering-suppressing child phantom DR image (transition scattering kernel) is shown as 082, and a scattering-suppressing child phantom DR image (conventional scattering kernel) is shown as 083 in one embodiment.

Please refer to fig. 9, which is a diagram illustrating the comparison between the effects of the LR-based deconvolution algorithm and the conventional deconvolution algorithm, fig. 091 is a scatter-suppressed-head virtual DR image (LR deconvolution), fig. 092 is a scatter-suppressed-head virtual DR image (conventional deconvolution), fig. 093 is a non-scattered-head virtual DR image, and fig. 094 is a scatter-containing-head virtual DR image.

Referring to fig. 10, a block diagram of an image processing apparatus according to another embodiment, which is configured to process an image acquired by digital X-rays by applying the image processing method as described above, specifically includes an initial parameter acquiring unit 10, a parameter compensating unit 20, and an image output unit 30. The initial parameter acquiring unit 10 is configured to acquire initial parameters of each pixel point in an imaging image of a digital X-ray apparatus according to the imaging image. In one embodiment, the initial parameters include pixel gray values, pixel depth values, and/or pixel sizes. The parameter compensation unit 20 is configured to perform scattering effect compensation on an initial parameter of each pixel of the imaging image, where a compensation value for performing scattering effect compensation on the initial parameter of each pixel is obtained according to a preset mathematical model for ray interference compensation, and the mathematical model for ray interference compensation is obtained by evaluating and modeling a scattering signal of the digital X-ray device. The image output unit 30 is configured to take the imaging image subjected to the scattering effect compensation as an output image of the digital X-ray apparatus.

In the image processing method disclosed by the embodiment of the application, the proposed asymmetric shift scattering kernel is closer to the actual physical process, so that the scattering signal estimation is more accurate, the final scattering removal effect is more uniform, and the bright and dark fluctuation phenomenon does not exist in the output image. The method improves the stability of scattering suppression of the grid technology, is different from the traditional semi-empirical deconvolution algorithm in the scheme, has strict convergence stability in the LR algorithm used in the embodiment of the application, has good effect on the large-volume chest and abdomen DR image with the most serious scattering, can simultaneously suppress the noise method phenomenon in the scattering suppression image, and also adapts to the bilateral filtering regular term to achieve the effect of suppressing the noise while suppressing the scattering signal. In addition, an automatic parameter optimization method is provided in an embodiment, and accurate optimization of scattering model parameters can be completed by using a few special experimental data. When parameter setting is carried out in an experimental mode, tedious repetition is not needed (30-40 imaging is needed in the prior art, fluctuation errors of low-intensity radiation detection are easily received by repeated experiments), the application tries to measure related signals under the condition of high radiation intensity, and reduces the needed measurement times (the parameter setting can be completed by at least 4 imaging) while improving the measurement data precision.

According to the image processing method disclosed by the embodiment of the application, initial parameters of each pixel point in an imaging image are acquired according to the imaging image of a digital X-ray device; then, carrying out scattering effect compensation on initial parameters of each pixel point of the imaging image; and finally, taking the imaging image subjected to scattering effect compensation as an output image of the digital X-ray device. The compensation value for performing scattering effect compensation on the initial parameter of each pixel point is obtained according to a ray interference compensation mathematical model. The radiation interference compensation mathematical model is obtained by evaluating and modeling the scattering signal of the digital X-ray device, so that the influence of the scattered radiation on the digital X-ray image can be eliminated to the greatest extent.

Those skilled in the art will appreciate that all or part of the functions of the various methods in the above embodiments may be implemented by hardware, or may be implemented by a computer program. When all or part of the functions in the above embodiments are implemented by means of a computer program, the program may be stored in a computer readable storage medium, and the storage medium may include: read-only memory, random access memory, magnetic disk, optical disk, hard disk, etc., and the program is executed by a computer to realize the above-mentioned 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 functions described above can be realized. In addition, when all or part of the functions in the above embodiments are implemented by means of a computer program, the program may be stored in a storage medium such as a server, another computer, a magnetic disk, an optical disk, a flash disk, or a removable hard disk, and the program in the above embodiments may be implemented by downloading or copying the program into a memory of a local device or updating a version of a system of the local device, and when the program in the memory is executed by a processor.

The foregoing description of the invention has been presented for purposes of illustration and description, and is not intended to be limiting. Several simple deductions, modifications or substitutions may also be made by a person skilled in the art to which the invention pertains, based on the idea of the invention.

Claims (10)

1. An image processing method for a digital X-ray virtual grid, comprising:

acquiring initial parameters of each pixel point in an imaging image of a digital X-ray device according to the imaging image; the initial parameters include pixel gray values, pixel depth values, and/or pixel sizes;

performing scattering effect compensation on the initial parameters of each pixel point of the imaging image, and taking the imaging image subjected to the scattering effect compensation as an output image of the digital X-ray device; the compensation value for performing scattering effect compensation on the initial parameter of each pixel point is obtained according to a preset ray interference compensation mathematical model, and the ray interference compensation mathematical model is obtained by evaluating and modeling a scattering signal of the digital X-ray device.

2. The image processing method of claim 1, wherein the ray interference compensation mathematical model is:

；

；

；

Where x and y are the pixel coordinates in the image,for projection of an image a single pixel position vector, x _c And y _c For projecting coordinates +.>To->Describing coordinates for a scattered ray function of the center; i _p A parameter value representing each pixel point generated from the primary ray signal as a primary ray function; i _s A parameter value representing each pixel point generated by the scattered ray signal as a scattered ray function; the initial parameter value of each pixel point in the imaging image is the sum of the parameter values generated for the same pixel point due to the primary ray signal and the scattered ray signal.

3. The image processing method according to claim 2, further comprising:

according to the scattered ray function I _s And primary ray function I _p Solving the ray interference compensation mathematical model, and scattering ray function I _s And primary ray function I _p Is the product of the scattered signal intensity function and the signal distribution function, and is expressed as:

；

where h is the scattered ray function I _s And primary ray function I _p Is a function of the relation of (c),scattering signal intensity function for small angle scattering signal, < >>As a function of scattered signal strength for large angle scattered signals.

4. The image processing method according to claim 3, further comprising:

the scatter signal intensity function for a small angle scatter signal is expressed as:

；

wherein,the intensity coefficient is obtained by preset fitting; />Is a scattering effect threshold cut-off function; i ₀ Is the air projection value, and ∈> 。

5. The image processing method as claimed in claim 4, further comprising:

the distribution function of the scattering signal intensity function for the large angle scattering signal is as follows:

；

wherein,to fit the resulting peak width coefficients, d is the distribution range for controlling the scattered ray signal；/>Is a two-dimensional Zernike orthogonal basis function +.>Is a polynomial coefficient.

6. The image processing method according to claim 5, further comprising:

will measure ray function I, scatter ray function I _s Primary ray function I _p Respectively spread into vector form，，/>Scattered ray function I _s And primary ray function I _p There is a relationship between:

；

wherein S is _jk As elements of the scattering matrix S to be obtained, which represents allFor a certain pixel j +.>Is a component part of (a).

7. The image processing method according to claim 6, further comprising:

Each column of elements in the scattering matrix S is not the same shiftBased on this measurement radiation function I +.>The expression of (2) is:

；

wherein, the imaging matrix h=diag {1, …,1} +s, that is, the scattering matrix is a block diagonal matrix with one element being 1, all elements are non-negative, and the values of a large number of elements far from the main diagonal gradually decrease and tend to be zero.

8. The image processing method according to claim 7, further comprising:

a meta heuristic iterative optimization algorithm of generalized mode search is applied to respectively calculate corresponding scattered ray functions I of image data obtained by at least two test die body experiments _s And optimizing the ray interference compensation mathematical model by a multi-objective function optimization method based on a wolf-swarm model.

9. A computer-readable storage medium, characterized in that the medium has stored thereon a program executable by a processor to implement the image processing method according to any one of claims 1 to 8.

10. An image processing device for a digital X-ray virtual grid, characterized by applying an image processing method according to any of claims 1 to 8 for post-processing an image acquired by a digital X-ray device, the image processing device comprising:

An initial parameter obtaining unit, configured to obtain initial parameters of each pixel point in an imaging image of a digital X-ray device according to the imaging image; the initial parameters include pixel gray values, pixel depth values, and/or pixel sizes;

the parameter compensation unit is used for carrying out scattering effect compensation on the initial parameters of each pixel point of the imaging image; the compensation value for performing scattering effect compensation on the initial parameter of each pixel point is obtained according to a preset ray interference compensation mathematical model, and the ray interference compensation mathematical model is obtained by evaluating and modeling a scattering signal of the digital X-ray device;

and the image output unit is used for taking the imaging image subjected to the scattering effect compensation as an output image of the digital X-ray device.