X-ray spectral estimation method with scatter correction

Abstract

The application discloses an X-ray energy spectrum estimation method with scatter correction, which comprises the following steps: scanning a sample to obtain original CT projection data; constructing a sample model, simulating a scatter signal of the sample model, and estimating a scatter signal caused by a detector shell; performing scatter correction on the original CT projection data to obtain scatter-corrected projection data; calculating a model thickness of the sample model under each angle; and performing energy spectrum estimation according to the scatter-corrected projection data and the model thickness under each angle. The X-ray detection energy spectrum can be more accurately and quickly estimated by using the X-ray energy spectrum estimation method with scatter correction.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of X-ray technology, and in particular to an X-ray spectrum estimation method with scatter correction. BACKGROUND

[0002] In most clinical CT scanners, the X-ray source is often polychromatic. The projection data detected by the detector is affected by the X-ray source spectrum distribution and the detector response, in addition to the effect of the scattered signal.

[0003] The following model is often used in the prior art to describe the attenuation of X-rays passing through the measured object:

[0004]

[0005] where I0(l) and I(l) are the X-ray signals before and after passing through the measured object along path l, Ω E is the energy value space, φ d (E): = φ(E)R(E) is the energy spectrum detected by the detector, E represents the energy value, φ(E) is the X-ray spectrum, R(E) is the detector response function, μ tot (x, E) is the linear attenuation coefficient, Sc l is the scattered signal intensity received by the detector pixel corresponding to path l.

[0006] As can be seen from formula (1), in addition to considering the effect of the scattered signal, accurate reconstruction of the CT image (i.e. accurate estimation of the linear attenuation coefficient of the measured object) also requires prior estimation of the accurate X-ray spectrum.

[0007] Pre-estimating the accurate X-ray spectrum distribution and applying the spectrum for projection data correction can reduce the beam hardening artifacts in CT. The methods for estimating the X-ray spectrum can be divided into the following four categories. One is to directly use a photon counting detector to measure the X-ray spectrum, but the energy response of the detector cannot be accurately known, which will lead to inaccurate estimation. Two is the model-based method, which uses empirical or semi-empirical physical models to generate the spectrum. Three is to express the pre-estimated spectrum as a linear combination of different weights of the pre-generated spectrum. Four is to estimate the X-ray spectrum in CT from the measured projection data. The fourth method is usually divided into two steps: (i) measuring the transmission data of X-rays passing through different thicknesses of known materials (such as uniform aluminum); (ii) solving the linear equation representing the attenuation process of X-rays to reconstruct the spectrum. In 2005, Sidky et al. proposed using the expectation-maximization (EM) algorithm to solve the linear integral equation for spectrum estimation.

[0008] When there is "missing data" or "incomplete data", it is complicated to solve the parameter value that maximizes the likelihood function. Rubin et al. proposed a general method for iterative computation of maximum likelihood estimation for incomplete data in 1997, which is the EM algorithm. The outline of the EM algorithm is as follows:

[0009] (i) Let and be two sample spaces, and there exists a mapping from to Let be the observed incomplete data, and be the complete data, whose probability density is f(x|θ), θ∈Θ, and Θ is the parameter value space. Let Y=T(X), and the probability density of Y is

[0010]

[0011] where The purpose of the EM algorithm is to solve the parameter θ that maximizes ln(g(y|θ)).

[0012] (ii) Solve the expectation (E step) of the likelihood function, i.e.

[0013]

[0014] i=0,1,2,... is the iteration number.

[0015] (iii) Solve the maximum point of the expectation of the log-likelihood function (M step), i.e.

[0016] (iv) Take the initial estimate

[0017] After that, repeat the E step and the M step to get an EM estimation sequence and satisfy

[0018]

[0019] Sidky et al. proposed to estimate the spectrum from the transmission data using EM in 2005, and the transmission data in the CT imaging system is described by the following formula (6)

[0020]

[0021] where Ω E is the value space of the spectrum φ(E), and φ d(E) = φ(E)R(E) is the detector spectrum, and R(E) is the detector response function. Under certain discrete conditions, equation (6) can be discretized as follows

[0022]

[0023]

[0024] where I is the sampled spectrum number, J is the index set of the X-ray path, and {w i} is φ d (E) in the selected basis expansion coefficient

[0025]

[0026] According to the research of Sidky et al., the EM algorithm is used to iteratively solve equation (7), and the following equation is obtained

[0027]

[0028] where k is the iteration number, the initial value is obtained by equation (9) after the initial spectrum is given.

[0029] The experimental results show that for the projection data without scattering, the EM algorithm is an accurate and robust spectrum estimation method. However, the scattering signal will affect the accuracy of the X-ray spectrum estimation, so it is urgent to provide an X-ray spectrum estimation scheme considering scattering. SUMMARY

[0030] The purpose of the present application is to provide an X-ray spectrum estimation method with scattering correction to overcome or at least alleviate at least one of the above-mentioned defects of the prior art.

[0031] To achieve the above purpose, the present application provides an X-ray spectrum estimation method with scattering correction, comprising:

[0032] Step 1, scanning the sample to obtain the original CT projection data;

[0033] Step 2, constructing a sample model, simulating the scattering signal of the sample model, and estimating the scattering signal caused by the detector shell;

[0034] Step 3, performing scattering correction on the original CT projection data to obtain the scattering corrected projection data

[0035] where,

[0036] denotes the transmission data after correction of the scattered signal caused by the detector housing and the object under test; I * (l j ) and denote the intensities of the X-ray signals before and after passing through the sample after correction of the scattered signal caused by the detector housing, respectively; Sc j is the scattered signal intensity received by the detector pixel corresponding to the path l j ; Ω E is the preset energy value space; φ d (E) is the energy spectrum detected by the detector, φ(E) is the X-ray energy spectrum, R(E) is the detector response function, μ tot (x, E) is the linear attenuation coefficient; l j denotes the X-ray path, j denotes the jth X-ray path, and J denotes the index set of the X-ray paths;

[0037] Step 4, calculate the model thickness of the sample model at each angle;

[0038] Step 5, perform energy spectrum estimation according to the scatter-corrected projection data and the model thickness at each angle; including:

[0039] Discretize equation (11) to obtain equations (12) and (13):

[0040]

[0041]

[0042] where I is the number of sampled energy spectra; i denotes the ith energy spectrum; is a preset energy spectrum basis; w i is the expansion coefficient of the energy spectrum φ d (E), as follows:

[0043]

[0044] Iteratively solve w i using the following equation (15):

[0045]

[0046] In equation (15), k is the iteration number, and the initial energy spectrum φ is given, and the initial value w can be determined by equation (14), and w i is obtained by iterative solution;

[0047] According to w i and equation (14), φ d (E) is obtained, and then φ d(E) = φ(E)R(E) to obtain the X-ray spectrum φ(E).

[0048] Preferably, the method further comprises:

[0049] Step 6: judging whether the estimated spectrum meets a preset condition; if yes, taking the estimated spectrum as the detected spectrum of the X-ray; otherwise, returning to step 2.

[0050] Preferably, the preset condition is that the estimated X-ray detected spectrum φ d (E) has a root mean square error with the reference spectrum less than a threshold.

[0051] Preferably, step 4 comprises:

[0052] Let the sample model be M, and let x represent a point in the sample model, and the integral value of δ(x) on each X-ray path l j is calculated as the model thickness corresponding to the angle of the X-ray path l j .

[0053] The present application has the following advantages due to the above technical solutions:

[0054] The X-ray spectrum estimation scheme with scatter correction provided by the present application simultaneously considers the influence of the scatter signals caused by the detector and the sample itself on the projection data, and uses an EM algorithm with scatter correction, so that the detected spectrum of the X-ray can be more accurately and quickly estimated. BRIEF DESCRIPTION OF DRAWINGS

[0055] Figure 1 FIG. 1 is a flowchart of the X-ray spectrum estimation method with scatter correction provided by the present application.

[0056] Figure 2 (a) is a geometric diagram of a uniform aluminum block generated by simulation data in an example of the method provided by the present application.

[0057] Figure 2 (b) is a given reference spectrum and an initial spectrum in an example of the method provided by the present application.

[0058] Figure 3 FIG. 5 is an experimental result of spectrum estimation of simulation data given in an example of the method provided by the present application.

[0059] Figure 4 FIG. 6 is a corresponding relationship diagram of the polychromatic projection value of the detected spectrum estimated from the simulation data given in an example of the method provided by the present application and the intersection length of the uniform aluminum block.

[0060] Figure 5A schematic diagram of a uniform Al block of actual data given in an example of the method provided by the present invention.

[0061] Figure 6 (a) is a schematic diagram of the estimated detection energy spectrum in an example using the method provided by the present invention.

[0062] Figure 6 (b)-(d) are schematic diagrams of multi-color projection curves of the detection energy spectrum in an example using the method provided by the present invention. DETAILED DESCRIPTION

[0063] In the accompanying drawings, the same or similar reference numerals are used to represent the same or similar elements or elements with the same or similar functions. The embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0064] In the description of the present invention, the terms "center", "longitudinal", "lateral", "front", "back", "left", "right", "vertical", "horizontal", "top", "bottom", "inside", "outside" and the like to indicate directions or positional relationships based on the directions or positional relationships shown in the accompanying drawings. They are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific direction, be constructed and operated in a specific direction. Therefore, they should not be understood as limiting the scope of protection of the present invention.

[0065] In the absence of conflict, the technical features in the various embodiments and implementations of the present invention may be combined with each other and are not limited to the embodiments or implementations to which the technical features belong.

[0066] The present invention is further described below with reference to the accompanying drawings and specific embodiments. It should be pointed out that the technical solution and design principle of the present invention are described in detail below only with an optimized technical solution, but the protection scope of the present invention is not limited to this.

[0067] The following terms are used herein, and their meanings are explained below for ease of understanding. Those skilled in the art will appreciate that the following terms may have other names, but any other names should be considered consistent with the terms listed herein without departing from their meanings.

[0068] The present invention provides an X-ray energy spectrum estimation method with scatter correction (EM-SC), which is applicable to projection data contaminated by scatter signals. Figure 1 As shown, the method includes:

[0069] Step 1: Scan the sample to obtain original CT projection data.

[0070] The sample (or phantom) is a uniform material object such as an aluminum cylinder or an aluminum cuboid. Under common conditions, the original CT projection data of the sample is scanned. Given the initial energy spectrum information and geometric parameters.

[0071] Step 2: construct a sample model, simulate the scattered signal of the sample model, and estimate the scattered signal caused by the detector shell.

[0072] After obtaining the CT projection data, CT reconstruction can be performed. This can be performed using the FDK algorithm, or other algorithms, which are not limited in this article. The sample model can be obtained by performing threshold segmentation on the results of the CT projection data reconstruction using the FDK reconstruction algorithm.

[0073] The scattering signal simulation of the sample model can estimate the scattering signal caused by the detector housing. For example, the gQMCFFD method can be used to simulate the scattering signal of the sample model, and the convolution method can be used to estimate the scattering signal caused by the detector housing. The specific scheme of using gQMCFFD to simulate the scattering signal of the model is described in detail in the invention patent No. 202111107623.5, which will not be repeated in this article. Estimating the scattering signal caused by the detector housing using the convolution method can also be achieved using existing technology. For example, the article Efficient scatter correction using asymmetric kernels[J] published by J.Star-Lack, M.Sun, A.Kaestner et al. in Physics of Medical Imaging specifically introduces a method of estimating the scattering signal caused by the detector housing using the convolution method, which will not be repeated in this article.

[0074] Step 3: Perform scatter correction on the CT projection data to obtain the scatter corrected projection data.

[0075] in,

[0076] Represents the transmission data after correcting the scattered signals caused by the detector housing and the object being measured; I * (l j )and They represent the X-ray signals after correcting the scattered signals caused by the detector housing; Sc j is the path l j The intensity of the scattered signal received by the corresponding detector pixel; Ω E is the preset energy value space; φ d (E): =φ(E)R(E) is the energy spectrum detected by the detector, φ(E) is the X-ray energy spectrum, R(E) is the detector response function, μ tot (x,E) linear attenuation coefficient; l jX-ray path, j represents the jth X-ray path, J represents the index set of X-ray paths. R(E) is a preset function related to a specific detector, and the value range is [0, 1].

[0077] Step 4, calculate the model thickness of the sample model at each angle.

[0078] This step includes: let the sample model be M, and x represents a point in the sample model, and the integral value of δ(x) at each X-ray path l j is calculated as the model thickness at the corresponding angle of each X-ray path l j .

[0079] Step 5, according to the projection data after scattering correction and the model thickness at each angle, perform energy spectrum estimation; including:

[0080] Discretize equation (11) to obtain equations (12) and (13):

[0081]

[0082]

[0083] where I is the sampled energy spectrum number; i represents the ith energy spectrum; is a preset energy spectrum basis; w i is the expansion coefficient of the energy spectrum φ d (E), as follows:

[0084]

[0085] Iteratively solve w i using the following equation (15):

[0086]

[0087] In equation (15), k is the iteration number, and after giving the initial energy spectrum , the initial value can be determined by equation (14), and w is obtained by iterative solution.

[0088] According to and equation (14), the φ d (E) after the kth iteration is obtained, and then the X-ray energy spectrum φ(E) is obtained according to φ d (E): = φ(E)R(E).

[0089] The method can further include:

[0090] Step 6, judging whether the estimated energy spectrum meets preset conditions; if yes, taking the estimated energy spectrum as a final energy spectrum; otherwise, returning to step 2.

[0091] wherein the preset conditions are that the estimated X-ray detection energy spectrum φ d (E) a root mean square error with the reference energy spectrum is less than a threshold.

[0092] The application further provides an X-ray energy spectrum estimation device with scatter correction, which is used to execute the method in any of the above embodiments or examples. The device can be installed in a CT scanner as a component, or connected with the CT scanner as an independent device.

[0093] In an example, the device comprises an input interface, a processor, and an output interface, the processor is used to execute each step in the method in any of the above embodiments or examples, the input interface is used to receive external information required by the processor for X-ray energy spectrum estimation, and the output interface is used to output the X-ray energy spectrum with scatter correction estimated by the processor.

[0094] The X-ray energy spectrum estimation scheme with scatter correction provided by the application can more accurately and quickly estimate the detection energy spectrum of X-rays.

[0095] Experimental results

[0096] This subsection gives numerical experiments of estimating energy spectrum by using the EM algorithm with scatter correction, to illustrate the accuracy and robustness of the EM algorithm with scatter correction (EM-SC). Here, the experimental results of simulated data and real data are given respectively. For simulated data, we first collect the attenuation data of X-rays along different paths through a uniform aluminum block with a volume of 80x80x80mm 3 . The specific operation is as follows: translate the aluminum block (45 degrees against the X-ray source and detector direction) to perform translation scanning. The geometric diagram of the uniform aluminum block is shown in Figure 2 (a). The collected attenuation data is contaminated with scatter signals and noise. Then, the gQMCFFD algorithm is used to estimate the scatter signals of X-rays along different paths through the uniform aluminum block, to correct the scatter signals in the projection data. Figure 2 (b) shows a given reference energy spectrum and an initial energy spectrum, and the reference energy spectrum is also called the assumed true energy spectrum.

[0097] Figure 3 The experimental results of energy spectrum estimation of simulated data are given. Figure 3 (a) is a comparison of the energy spectrum estimated by the EM algorithm after 128 iterations, the energy spectrum estimated by the EM algorithm with scatter correction (EM-SC) after 128 iterations, and the reference energy spectrum. From Figure 3(a) We can see that the estimated energy spectrum by the EM algorithm with scatter correction is highly consistent with the true energy spectrum, and the root mean square error (RMSE) between the estimated energy spectrum by the EM algorithm with scatter correction and the true energy spectrum is 0.001337, and the RMSE between the estimated energy spectrum by the EM method and the true energy spectrum is 0.00535, which shows the accuracy of the EM algorithm with scatter correction. Figure 3 (b) is the RMSE between the estimated energy spectrum by the EM algorithm and the estimated energy spectrum by the EM algorithm with scatter correction and the reference energy spectrum, and the trend of the RMSE with the number of iterations. From Figure 3 (b) we can see that the RMSE between the estimated energy spectrum by the EM method and the true energy spectrum increases with the number of iterations, which may be caused by the introduction of noise in the iteration process. This shows that the robustness of the EM algorithm is weak for data with scatter signal pollution. The RMSE between the estimated energy spectrum by the EM algorithm with scatter correction and the reference energy spectrum first decreases with the number of iterations, reaches the minimum at 128 iterations, and remains stable with the increase of the number of iterations, which shows the robustness of the EM algorithm with scatter correction.

[0098] To further illustrate the accuracy of the EM algorithm with scatter correction, we give the corresponding relationship diagram of the polychromatic projection value of the estimated detection energy spectrum of the simulated data and the intersection line length of the uniform aluminum block, as Figure 4 shown. Figure 4 (a) and (b) are the corresponding relationship between the polychromatic projection value without scatter correction and the intersection line length of the uniform aluminum block, respectively, and the corresponding relationship between the polychromatic projection value with scatter correction and the intersection line length of the uniform aluminum block, respectively, where the relationship diagram composed of points is obtained from the measured data, and the dashed line is obtained from the inverse solution of the detected detection energy spectrum. From Figure 4 (a) we can find that the scatter signal will cause the nonlinear reduction of the polychromatic projection value, and the effect will be more obvious with the increase of the intersection line length. This shows the necessity of correcting the scatter signal first, and also explains why the EM method with scatter correction is more accurate than the EM method for data with scatter. From Figure 4 (b) we find that the corresponding relationship curve between the polychromatic projection value calculated by the estimated energy spectrum by the EM algorithm with scatter correction and the intersection line length of the uniform aluminum block is highly fitted with the corresponding relationship curve between the polychromatic projection value measured and the intersection line length of the uniform aluminum block, which further illustrates the accuracy of the EM algorithm with scatter correction.

[0099] To estimate the detection energy spectrum of the actual data, this paper uses a 40x40x80mm 3 uniform Al block, and its schematic diagram is shown in Figure 5 . The estimated detection energy spectrum by the EM algorithm with scatter correction (128 iterations) and the EM method (128 iterations) is shown in Figure 6 (a). The corresponding polychromatic projection curve of the detection energy spectrum is shown in Figure 6(b)-(d) show. The results of the EM algorithm with scatter correction (EM-SC) are more accurate than the EM algorithm, but the EM algorithm with scatter correction (EM-SC) that takes into account the scatter from both the object under test and the detector housing is better than the EM algorithm with scatter correction (EM-SC0) that takes into account the scatter from only the object under test.

[0100] Finally, it should be pointed out that the above embodiments are only used to illustrate the technical solutions of the present application, and are not intended to limit the same. Those skilled in the art should understand that the technical solutions described in the foregoing embodiments can be modified, or some of the technical features can be replaced by equivalents; such modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. A method of X-ray energy spectrum estimation with scatter correction, characterized by, Comprising: Step 1, scanning the sample to obtain original CT projection data; Step 2, constructing a sample model, simulating the scattering signal of the sample model, and estimating the scattering signal caused by the detector shell; Step 3. Scatter correction is performed on the original CT projection data to obtain scatter corrected projection data wherein denotes the transmission data after correction of the scattered signal from the detector housing and the measured object; I * (l j ) and denote the intensity of the X-ray signal before and after passing through the sample, respectively, after correction of the scattered signal from the detector housing; Sc j is the scattered signal intensity received by the detector pixel corresponding to the path l j ; Ω E is the preset energy value space; φ d (E): = φ(E)R(E) is the energy spectrum detected by the detector, φ(E) is the X-ray energy spectrum, R(E) is the detector response function, μ tot (x, E) is the linear attenuation coefficient; l j denotes the X-ray path, j denotes the jth X-ray path, and J denotes the index set of the X-ray paths; Step 4, calculating the model thickness of the sample model under each angle; Step 5, performing energy spectrum estimation according to the scattering-corrected projection data and the model thickness under each angle; comprising: Discretizing formula (11) to obtain formula (12) and formula (13): Where, I is the number of sampled energy spectra; i represents the i-th energy spectrum; is the preset energy spectrum basis; w i is the energy spectrum φ d The expansion coefficient of (E) is as follows: The w is iteratively solved using the following equation (15) i : In Equation (15), k is the iteration number, and the initial energy spectrum is given After that, the initial value w can be determined by Equation (14) and is solved iteratively i ; According to w i and formula (14) φ d (E), then according to φ d (E): = φ(E) R(E) X-ray spectrum φ(E) is obtained.

2. The X-ray energy spectrum estimation method with scatter correction according to claim 1, characterized in that, Further comprising: Step 6, judging whether the estimated energy spectrum meets a preset condition; If yes, taking the estimated energy spectrum as the detection energy spectrum of X-rays; Otherwise, returning to step 2.

3. The X-ray energy spectrum estimation method with scatter correction according to claim 2, characterized in that: The preset condition is that the estimated X-ray detection energy spectrum φ d (E) a root mean square error with the reference energy spectrum is less than a threshold value.

4. The X-ray energy spectrum estimation method with scatter correction according to claim 1, characterized in that, Step 4 comprises: Let the sample model be M, and let x represents a point in the sample model, and δ(x) is calculated for each X-ray path l j The integral value of the X-ray path l j Model thickness corresponding to the angle.

Images