Multi-energy CT scattering correction method and device based on Boltzmann transport equation

Abstract

The invention relates to the technical field of radiation imaging, in particular to a multi-energy CT scattering correction method and device based on a Boltzmann transport equation, and the method comprises the steps: carrying out the preprocessing of to-be-corrected initial multi-energy projection data obtained after the multi-energy CT scanning of an object, and obtaining the initial correction multi-energy projection data; performing material density conversion on the initial correction multi-energy projection data to obtain material density distribution data; discretization and label marking are carried out on an actual scanning energy spectrum of an object in the multi-energy CT scanning process, and an energy spectrum matrix in a pre-constructed Boltzmann transport equation with labels is determined; calculating multi-energy scattering signal data of the object; and obtaining actual multi-energy projection data by using the multi-energy scattering signal data. Therefore, the problems that the density of the material obtained by substance decomposition is distorted due to the fact that scattering correction is not fine enough, multiple scattering estimations need to be carried out under different X-ray source conditions due to multiple energy spectrums, the processing process is tedious and complex and the like in the prior art are solved.

Description

Technical Field

[0001] The present application relates to the field of radiation imaging technology, and in particular to a multi-energy CT scatter correction method and device based on the Boltzmann transport equation. Background Art

[0002] CT (Computed Tomography) is a commonly used non-destructive testing technology that can obtain internal structure information without damaging the object being tested. Among them, ME-CBCT (Multi Energy Cone Beam Computed Tomography)

[0003] It is a type of CT that has developed rapidly in recent years. Compared with traditional single-energy CBCT, it can provide stronger material differentiation and more material structure information, and has been widely used in medical diagnosis and security inspection.

[0004] In the related technology, based on Beer's law and Compton scattering formula, combined with the physical parameter database such as the attenuation coefficient and scattering factor of the material, Monte Carlo simulation can be performed on the multi-energy CBCT scanning process to obtain the corresponding scattering signal, thereby realizing scattering correction; in view of the low frequency of the scattering signal, the frequency domain distribution of the scattering signal on the image can be modeled, and a scattering kernel can be constructed, and the scattering correction can be realized by combining the iterative algorithm; it is also possible to use the actual data or simulated data, use the scattering image as the training set, and the non-scattering image as the label, and use the deep learning technology to train the relevant network, and then obtain the corresponding scattering signal, thereby realizing scattering correction and reducing scattering artifacts.

[0005] However, in the related art, since the material decomposition counts used in spectral imaging are very sensitive to the energy spectrum, some simple descattering methods, even if they can remove most of the scattering, will cause serious distortion of the obtained material density. In addition, due to the multi-energy spectrum, multiple scattering estimates need to be performed under different X-ray source conditions. The processing process is lengthy and complicated and urgently needs to be improved. Summary of the invention

[0006] The present application provides a multi-energy CT scattering correction method and device based on the Boltzmann transport equation to solve the problems in the related art that, since material decomposition is very sensitive to the energy spectrum, some simple descattering methods can remove most of the scattering, but will cause serious distortion of the material density obtained by material decomposition, and due to the multi-energy spectrum, multiple scattering estimates need to be performed under different X-ray source conditions, and the processing process is lengthy and complicated.

[0007] The first aspect of the present application provides a multi-energy CT scatter correction method based on the Boltzmann transport equation, comprising the following steps: preprocessing the initial multi-energy projection data to be corrected obtained after the object is scanned by multi-energy computed tomography (CT) to obtain the initial corrected multi-energy projection data after the preprocessing of the initial multi-energy projection data to be corrected; performing material density conversion on the initial corrected multi-energy projection data to obtain material density distribution data after the material density conversion of the initial corrected multi-energy projection data; discretizing and labeling the actual scanning energy spectrum of the object during the multi-energy CT scanning process to determine the energy spectrum matrix in the pre-constructed labeled Boltzmann transport equation based on the actual scanning energy spectrum after discretization and labeling; calculating the multi-energy scattering signal data of the object in the multi-energy CT scanning based on the final scattering signal intensity equation in the pre-constructed labeled Boltzmann transport equation, the energy spectrum matrix and the material density distribution data; and performing scattering correction on the initial multi-energy projection data using the multi-energy scattering signal data to obtain the actual multi-energy projection data of the object after scattering correction.

[0008] Optionally, in one embodiment of the present application, before calculating the multi-energy scattering signal data of the object in the multi-energy CT scan based on the final scattering signal intensity equation in the pre-constructed labeled Boltzmann transport equation, the energy spectrum matrix and the material density distribution data, it also includes: constructing a static Boltzmann transport equation in the multi-energy CT scanning process; obtaining the photon flux integral expression and the photon flux differential expression of the photon flux and the photon source term in the static Boltzmann transport equation, the photon source term integral expression and the photon source term differential expression of the photon flux in the multi-energy CT scanning process; based on the photon flux integral expression, the photon flux differential expression, the photon source term integral expression and the photon source term differential expression, obtaining the photon flux discrete expression and the photon source term of the photon flux The invention relates to a method for obtaining a discrete expression of a photon source term of a static Boltzmann transport equation. The method comprises the following steps: generating a discrete expression of a photon source term of a static Boltzmann transport equation based on the photon flux integral expression, the photon flux differential expression, the photon source term integral expression and the photon source term differential expression, generating an initial scattering signal intensity equation of the scattering signal during the multi-energy CT scanning process based on the photon flux integral expression, the photon flux differential expression, the photon source term integral expression and the photon source term differential expression; discretizing the initial scattering signal intensity equation using the photon flux discrete expression and the photon source term discrete expression to obtain a final scattering signal intensity equation after the initial scattering signal intensity equation is discretized; determining at least one label dimension of the static Boltzmann transport equation based on the photon flux discrete expression, the photon source term discrete expression, the final scattering signal intensity equation and the at least one label dimension; and constructing a labeled Boltzmann transport equation based on the photon flux discrete expression, the photon source term discrete expression, the final scattering signal intensity equation and the at least one label dimension.

[0009] Optionally, in one embodiment of the present application, the actual scanning energy spectrum of the object during the multi-energy CT scanning process is discretized and labeled to determine the energy spectrum matrix in the pre-constructed labeled Boltzmann transport equation based on the actual scanning energy spectrum after discretization and labeling, including: judging whether the number of the actual scanning energy spectrum after the discretization and labeling is greater than half of the number of preset discrete energy groups; if the number of the actual scanning energy spectrum after the discretization and labeling is greater than half of the number of preset discrete energy groups, then determining the energy spectrum matrix in the pre-constructed labeled Boltzmann transport equation based on the preset discrete energy groups. The dimension and value of the energy spectrum matrix in the pre-constructed labeled Boltzmann transport equation are determined based on the number of actual scanning energy spectra after discretization and labeling, wherein the energy spectrum matrix is ​​a unit matrix; if the number of actual scanning energy spectra after discretization and labeling is less than or equal to half of the number of preset discrete energy groups, the dimension and value of the energy spectrum matrix in the pre-constructed labeled Boltzmann transport equation are determined based on the number of actual scanning energy spectra after discretization and labeling, wherein each column vector in the energy spectrum matrix is ​​an energy spectrum vector of each actual scanning energy spectrum after discretization and labeling.

[0010] Optionally, in one embodiment of the present application, the expression of the pre-constructed labeled Boltzmann transport equation may be, but is not limited to,:

[0011]

[0012] in, is the angular flux of photons at voxel i, with energy group g, emission direction m, label l and scattering order n; is voxel i w′ The average linear attenuation coefficient for the energy group g at position i w′ is the voxel number where the w′ path is located; w′ is a vector The path represented passes through each voxel, and the path length in each voxel exists as an accumulated variable; The intensity of the photon source at voxel i, with energy group g, emission direction m, label l and scattering order n+1; μ s；i,g′→g,m′→m is the scattering cross section of a photon with a moving direction of m′ discrete angle and energy group g′ at voxel i; is the angular flux of photons at voxel i, with energy group g′, emission direction m′, label l and scattering order n, is the position of voxel i, the energy group is g, and the emission direction is m (the m direction is given by The direction of the photon source is determined by the label l and the scattering order is k, A d is the area of ​​pixel d, θ is The angle between the vector and the normal vector of pixel d, E g is the average energy of energy group g, D(E g ) is the detector at energy E g Lower detector gain, is the spatial position vector of voxel i, is the spatial position vector of the detector pixel d, i is the discrete voxel index, which identifies the discrete spatial position; g is the energy group index, which identifies the position of different discrete energy groups; m is the discrete solid angle index, which identifies different spatial directions; l is the label index, which is used to identify photons from different energy spectra; k is the scattering order, which exists as a summation variable in the formula; N is the total number of discretized voxels; G is the number of discretized energy groups, K is the highest scattering order selected for calculation; n is the scattering order (that is, the photon is scattered n times); n+1 is the scattering order (that is, the photon is scattered n+1 times).

[0013] Optionally, in one embodiment of the present application, the expression of the final scattered signal intensity equation may be, but is not limited to,:

[0014]

[0015] In a second aspect, an embodiment of the present application provides a multi-energy CT scatter correction device based on the Boltzmann transport equation, comprising: a preprocessing module, used to preprocess the initial multi-energy projection data to be corrected obtained after the object is scanned by multi-energy CT, so as to obtain the initial corrected multi-energy projection data after the preprocessing of the initial multi-energy projection data to be corrected; a material density conversion module, used to perform material density conversion on the initial corrected multi-energy projection data, so as to obtain the material density distribution data after the material density conversion of the initial corrected multi-energy projection data; an energy spectrum discrete marking module, used to perform material density conversion on the object during the multi-energy CT scanning process; The present invention relates to a method for discretizing and labeling an actual scanning energy spectrum, so as to determine an energy spectrum matrix in a pre-constructed labeled Boltzmann transport equation based on the actual scanning energy spectrum after discretization and labeling; a calculation module, for calculating multi-energy scattering signal data of the object in a multi-energy CT scan based on a final scattering signal intensity equation in the pre-constructed labeled Boltzmann transport equation, the energy spectrum matrix and the material density distribution data; and a scattering correction module, for performing scattering correction on the initial multi-energy projection data using the multi-energy scattering signal data, so as to obtain actual multi-energy projection data of the object after scattering correction.

[0016] Optionally, in one embodiment of the present application, it also includes: a first construction module, which is used to construct a static Boltzmann transport equation in the multi-energy CT scanning process before calculating the multi-energy scattering signal data of the object in the multi-energy CT scanning based on the final scattering signal intensity equation in the pre-constructed labeled Boltzmann transport equation, the energy spectrum matrix and the material density distribution data; a first acquisition module, which is used to obtain the photon flux integral expression and the photon flux differential expression of the photon flux, the photon source term integral expression and the photon source term differential expression of the photon flux in the static Boltzmann transport equation, and the second acquisition module, which is used to obtain the photon flux discrete expression and the photon source term of the photon flux based on the photon flux integral expression, the photon flux differential expression, the photon source term integral expression and the photon source term differential expression of the photon source term. a discrete expression of the photon source term; a generating module for generating an initial scattering signal intensity equation of the scattering signal during the multi-energy CT scanning process based on the integral expression of the photon flux, the differential expression of the photon flux, the integral expression of the photon source term and the differential expression of the photon source term; a correcting module for discretizing the initial scattering signal intensity equation using the discrete expression of the photon flux and the discrete expression of the photon source term to obtain a final scattering signal intensity equation after the initial scattering signal intensity equation is discretized; a determining module for determining at least one label dimension of the static Boltzmann transport equation based on the discrete expression of the photon flux, the discrete expression of the photon source term and the final scattering signal intensity equation; and a second constructing module for constructing a labeled Boltzmann transport equation based on the discrete expression of the photon flux, the discrete expression of the photon source term, the final scattering signal intensity equation and the at least one label dimension.

[0017] Optionally, in one embodiment of the present application, the energy spectrum discrete marking module includes: a judgment unit, used to judge whether the number of actual scanning energy spectra after the discretization and labeling is greater than half of the number of preset discrete energy groups; a first determination unit, used to determine the dimension and value of the energy spectrum matrix in the pre-constructed labeled Boltzmann transport equation based on the number of preset discrete energy groups when the number of actual scanning energy spectra after the discretization and labeling is greater than half of the number of preset discrete energy groups, wherein the energy spectrum matrix is ​​a unit matrix; a second determination unit, used to determine the dimension and value of the energy spectrum matrix in the pre-constructed labeled Boltzmann transport equation based on the number of actual scanning energy spectra after the discretization and labeling when the number of actual scanning energy spectra after the discretization and labeling is less than or equal to half of the number of preset discrete energy groups, wherein each column vector in the energy spectrum matrix is ​​an energy spectrum vector of each actual scanning energy spectrum after discretization and labeling.

[0018] Optionally, in one embodiment of the present application, the expression of the pre-constructed labeled Boltzmann transport equation may be, but is not limited to,:

[0019]

[0020] in, is the angular flux of photons at voxel i, with energy group g, emission direction m, label l and scattering order n; is voxel i w′ The average linear attenuation coefficient for the energy group g at position i w′ is the voxel number where the w′ path is located; w′ is a vector The path represented passes through each voxel, and the path length in each voxel exists as an accumulated variable; The intensity of the photon source at voxel i, with energy group g, emission direction m, label l and scattering order n+1; μ s；i,g′→g,m′→m is the scattering cross section of a photon with a moving direction of m′ discrete angle and energy group g′ at voxel i; is the angular flux of photons at voxel i, with energy group g′, emission direction m′, label l and scattering order n, is the position of voxel i, the energy group is g, and the emission direction is m (the m direction is given by The direction of the photon source is determined by the label l and the scattering order is k, A d is the area of ​​pixel d, θ is The angle between the vector and the normal vector of pixel d, E g is the average energy of energy group g, D(E g ) is the detector at energy E g Lower detector gain, is the spatial position vector of voxel i, is the spatial position vector of the detector pixel d, i is the discrete voxel index, which identifies the discrete spatial position; g is the energy group index, which identifies the position of different discrete energy groups; m is the discrete solid angle index, which identifies different spatial directions; l is the label index, which is used to identify photons from different energy spectra; k is the scattering order, which exists as a summation variable in the formula; N is the total number of discretized voxels; G is the number of discretized energy groups, K is the highest scattering order selected for calculation; n is the scattering order (that is, the photon is scattered n times); n+1 is the scattering order (that is, the photon is scattered n+1 times).

[0021] Optionally, in one embodiment of the present application, the expression of the final scattered signal intensity equation may be, but is not limited to,:

[0022]

[0023] A third aspect of the present application provides an electronic device, comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the multi-energy CT scatter correction method based on the Boltzmann transport equation as described in the above embodiment.

[0024] A fourth aspect of the present application provides a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, which, when executed by a processor, implements the above-mentioned multi-energy CT scatter correction method based on the Boltzmann transport equation.

[0025] The fifth aspect of the present application provides a computer program product, including a computer program, which, when executed, implements the above-mentioned multi-energy CT scatter correction method based on the Boltzmann transport equation.

[0026] The embodiment of the present application can first pre-process the initial multi-energy projection data to be corrected obtained after the multi-energy CT scan of the object, and then obtain the initial corrected multi-energy projection number after pre-processing, and perform material density conversion to obtain the material density distribution data of the object, and then discretize and label the multi-energy spectrum data to obtain the spectrum matrix, and then combine the material distribution data and the spectrum matrix, use the pre-constructed labeled Boltzmann transport equation to calculate the multi-energy scattering signal data of the object in the multi-energy CT scan, and use the multi-energy scattering signal data to perform scattering correction on the initial multi-energy projection data to obtain the actual multi-energy projection data after scattering correction, calculate the multi-energy scattering signal data based on the pre-constructed labeled Boltzmann transport equation, and through the discretization of the phase space, convert the equation into an analytical solution, so that the program has excellent parallelism, thereby improving the calculation speed and portability. For the scattering correction requirements in multi-energy CT, an additional label dimension is introduced to realize the vectorization of the energy group space, thereby reducing the amount of calculation, and achieving a single calculation to obtain scattering signals under multiple energy spectra. This solves the problem in the related art that, since material decomposition is very sensitive to the energy spectrum, some simple descattering methods, even if they can remove most of the scattering, will cause serious distortion of the material density obtained by material decomposition, and due to the multi-energy spectrum, multiple scattering estimates need to be performed under different X-ray source conditions, resulting in a lengthy and complicated processing process.

[0027] Additional aspects and advantages of the present application will be given in part in the description below, and in part will become apparent from the description below, or will be learned through the practice of the present application. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] The above and / or additional aspects and advantages of the present application will become apparent and easily understood from the following description of the embodiments in conjunction with the accompanying drawings, in which:

[0029] Figure 1 A flowchart of a multi-energy CT scatter correction method based on the Boltzmann transport equation provided according to an embodiment of the present application;

[0030] Figure 2(a)-Figure 2(d) A schematic diagram showing a comparison between a scattering estimation result provided according to an embodiment of the present application and a Monte Carlo method using the Geant4 tool;

[0031] Figure 3 A schematic diagram of image comparison of projections collected on a multi-energy CT before and after scatter correction and fan beam standard values ​​according to an embodiment of the present application;

[0032] Figure 4 A flowchart of the working principle of a multi-energy CT scatter correction method based on the Boltzmann transport equation provided according to an embodiment of the present application;

[0033] Figure 5 It is a block diagram of a multi-energy CT scatter correction device based on the Boltzmann transport equation provided according to an embodiment of the present application;

[0034] Figure 6 It is a schematic diagram of the structure of an electronic device provided according to an embodiment of the present application. DETAILED DESCRIPTION

[0035] Embodiments of the present application are described in detail below, and examples of the embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are intended to be used to explain the present application, and should not be construed as limiting the present application.

[0036] The following describes the multi-energy CT scatter correction method and device based on the Boltzmann transport equation in an embodiment of the present application with reference to the accompanying drawings. In view of the problem mentioned in the above background technology that material decomposition is very sensitive to the energy spectrum, some simple descattering methods can remove most of the scattering, but will cause serious distortion of the material density obtained by material decomposition, and due to the multi-energy spectrum, it is necessary to perform multiple scattering estimates under different X-ray source conditions, and the processing process is lengthy and complicated. The present application provides a multi-energy CT scatter correction method based on the Boltzmann transport equation. In this method, the initial multi-energy projection data to be corrected obtained after the multi-energy CT scan of the object can be preprocessed, and then the initial corrected multi-energy projection data after preprocessing can be obtained, and the material density conversion is performed to obtain the object. The material density distribution data of the body is obtained, and then the multi-energy energy spectrum data is discretized and labeled to obtain the energy spectrum matrix. Then, the material density distribution data and the energy spectrum matrix are combined, and the pre-constructed labeled Boltzmann transport equation is used to calculate the multi-energy scattering signal data of the object in the multi-energy CT scan, and the multi-energy scattering signal data is used to perform scattering correction on the initial multi-energy projection data to obtain the actual multi-energy projection data after scattering correction. The multi-energy scattering signal data is calculated based on the pre-constructed labeled Boltzmann transport equation. Through the discretization of the phase space, the equation is converted into an analytical solution, which makes the program have superior parallelism, thereby improving the calculation speed and portability. For the scattering correction requirements in multi-energy CT, an additional label dimension is introduced to realize the vectorization of the energy group space, thereby reducing the amount of calculation, and achieving a single calculation to obtain scattering signals under multiple energy spectra. This solves the problem in the related art that, since material decomposition is very sensitive to the energy spectrum, some simple descattering methods, even if they can remove most of the scattering, will cause serious distortion of the material density obtained by material decomposition, and due to the multi-energy spectrum, multiple scattering estimates need to be performed under different X-ray source conditions, resulting in a lengthy and complicated processing process.

[0037] Specifically, Figure 1 The present invention provides a flowchart of a multi-energy CT scatter correction method based on the Boltzmann transport equation according to an embodiment of the present application.

[0038] like Figure 1 As shown, the multi-energy CT scatter correction method based on the Boltzmann transport equation includes the following steps:

[0039] In step S101, initial multi-energy projection data to be corrected obtained after multi-energy CT scanning of the object is preprocessed to obtain initial corrected multi-energy projection data after the preprocessing of the initial multi-energy projection data to be corrected.

[0040] As a possible implementation method, the embodiment of the present application can obtain multiple energy spectra of the object during the scanning process and the initial multi-energy projection data under each energy spectrum through multi-energy CT scanning, and then the embodiment of the present application preprocesses the initial multi-energy projection data, such as rough scattering correction and hardening correction, etc., which are not specifically limited in this application, to obtain the initial corrected multi-energy projection data after the initial projection data is preprocessed.

[0041] For example, the embodiment of the present application can obtain multiple energy spectra S by multi-energy CT scanning of the object. l ,l=1,2,…,L and the initial multi-energy projection data P under each energy spectrum l (N D ×1),l=1,2,3,…,L, let P=[P 1 ,P 2 ,…,P L ](N D ×L), perform rough scattering correction and hardening correction in the projection domain, and then obtain the preprocessed initial corrected multi-energy projection data P f .

[0042] For example, the present application adopts a three-dimensional rectangular coordinate system. In terms of voxel discretization, the object is discretized into 16×16×16 cube units of uniform size. In terms of three-dimensional angle, the 4π solid angle is divided into 16 longitudes and 9 latitudes at intervals of π / 8, which is a total of 114 discrete angles, according to the discretization method of longitude and latitude in geography. In terms of energy, 0~E is divided into 114 discrete angles according to the energy spectrum of the ray source. max The range is approximately evenly divided into 8 energy groups.

[0043] Furthermore, the embodiment of the present application can obtain two sets of initial multi-energy projection data of the object through multi-energy CT scanning, which can be represented by but not limited to p 1 ,p 2 , and then perform rough scatter correction and hardening correction processing on the initial multi-energy projection data, and then obtain the pre-processed initial corrected multi-energy projection data, which can be expressed as but not limited to

[0044] In step S102, material density conversion is performed on the initial corrected multi-energy projection data to obtain material density distribution data after material density conversion of the initial corrected multi-energy projection data.

[0045] For example, the embodiment of the present application can be a plurality of energy spectra S l ,l=1,2,…,L and the initial corrected multi-energy projection data P f , perform material density conversion, such as material decomposition, pre-reconstruction, etc., which is not specifically limited in this application, and then obtain the material density distribution data ρ of the object l,l=1,2,…,L, the material density distribution data can be used as input to the pre-built labeled Boltzmann transport equation.

[0046] In some embodiments, the embodiments of the present application can perform material decomposition based on the initial corrected multi-energy projection data, directly obtain multi-material density projection data, and reconstruct the material density distribution data of the object during the multi-energy CT scanning process.

[0047] For example, the embodiment of the present application can be based on the initial correction multi-energy projection data Material decomposition is performed to obtain multi-material density projection data, and then reconstruction is performed to obtain two kinds of material density distribution data, which can be expressed as but not limited to ρ 1 ,ρ 2 .

[0048] In some embodiments, the embodiments of the present application can perform material decomposition based on the initial corrected multi-energy projection data to obtain a virtual monoenergetic reconstructed image, and then convert the HU value data into material density distribution data of the object during the multi-energy CT scanning process based on the material conversion curve.

[0049] For example, the embodiment of the present application can firstly calibrate the multi-energy projection data based on the initial calibration Material decomposition is performed to obtain multi-material density projection data, and then weighted summation is performed to obtain virtual monoenergetic projection data. Then, three-dimensional reconstruction is performed to obtain a virtual monoenergetic reconstruction image. Through the material conversion curve, two kinds of material density distribution data are obtained, which can be expressed as but not limited to ρ 1 ,ρ 2 .

[0050] In some embodiments, the embodiments of the present application can perform pre-reconstruction based on the initial corrected projection data at a certain energy to obtain a three-dimensional reconstructed image, and then convert the HU value data into material density distribution data of the object during the multi-energy CT scanning process based on the material conversion curve.

[0051] For example, the embodiment of the present application can be based on the initial correction projection data at a certain energy Pre-reconstruction is performed to obtain a three-dimensional reconstructed image, and then the density distribution data of two materials are obtained through the material conversion curve, which can be expressed as but not limited to ρ 1 ,ρ 2 .

[0052] Optionally, in one embodiment of the present application, before calculating the multi-energy scattering signal data of the object in the multi-energy CT scan based on the final scattering signal intensity equation, the energy spectrum matrix and the material density distribution data in the pre-constructed labeled Boltzmann transport equation, the method further includes: constructing a static Boltzmann transport equation in the multi-energy CT scanning process; obtaining the photon flux integral expression and the photon flux differential expression of the photon flux, the photon source term integral expression and the photon source term differential expression of the photon flux and the photon source term in the static Boltzmann transport equation during the multi-energy CT scanning process; based on the photon flux integral expression, the photon flux differential expression, the photon source term integral expression and the photon source term differential expression, obtaining the photon flux discrete expression and the photon flux differential expression of the photon flux. Discrete expression of the photon source term of the photon source term; based on the integral expression of the photon flux, the differential expression of the photon flux, the integral expression of the photon source term and the differential expression of the photon source term, generate the initial scattering signal intensity equation of the scattering signal during the multi-energy CT scanning process; use the discrete expression of the photon flux and the discrete expression of the photon source term to discretize the initial scattering signal intensity equation to obtain the final scattering signal intensity equation after the initial scattering signal intensity equation is discretized; based on the discrete expression of the photon flux, the discrete expression of the photon source term and the final scattering signal intensity equation, determine at least one label dimension of the static Boltzmann transport equation; based on the discrete expression of the photon flux, the discrete expression of the photon source term, the final scattering signal intensity equation and at least one label dimension, construct a labeled Boltzmann transport equation. Among them, the expression of the pre-constructed labeled Boltzmann transport equation can be but is not limited to:

[0053]

[0054] in, is the angular flux of photons at voxel i, with energy group g, emission direction m, label l and scattering order n; is voxel i w′ The average linear attenuation coefficient for the energy group g at position i w′ is the voxel number where the w′ path is located; w′ is a vector The path represented passes through each voxel, and the path length in each voxel exists as an accumulated variable; The intensity of the photon source at voxel i, with energy group g, emission direction m, label l and scattering order n+1; μ s；i,g′→g,m′→m is the scattering cross section of a photon with a moving direction of m′ discrete angle and energy group g′ at voxel i; is the angular flux of photons at voxel i, with energy group g′, emission direction m′, label l and scattering order n, is the position of voxel i, the energy group is g, and the emission direction is m (the m direction is given by The direction of the photon source is determined by the label l and the scattering order is k, A d is the area of ​​pixel d, θ is The angle between the vector and the normal vector of pixel d, E g is the average energy of energy group g, D(E g ) is the detector at energy E g Lower detector gain, is the spatial position vector of voxel i, is the spatial position vector of the detector pixel d, i is the discrete voxel index, which identifies the discrete spatial position; g is the energy group index, which identifies the position of different discrete energy groups; m is the discrete solid angle index, which identifies different spatial directions; l is the label index, which is used to identify photons from different energy spectra; k is the scattering order, which exists as a summation variable in the formula; N is the total number of discretized voxels; G is the number of discretized energy groups, K is the highest scattering order selected for calculation; n is the scattering order (that is, the photon is scattered n times); n+1 is the scattering order (that is, the photon is scattered n+1 times).

[0055] It can be understood by those skilled in the art that the embodiments of the present application can construct the Boltzmann transport equation based on the transport and scattering problems of photons in CT scanning. Considering that the interaction time between photons and matter is very short, the relaxation time is ignored in the model, and it is assumed that the system reaches static equilibrium instantaneously. Based on the above description, the embodiments of the present application can construct the static Boltzmann transport equation in the multi-energy CT scanning process, which can be expressed as but not limited to:

[0056]

[0057] in, represents the position vector of the photon in three-dimensional space; represents the direction vector of the photon's motion in three-dimensional space; E represents the energy carried by the photon; Indicated in Position, direction of movement The flux of photons carrying energy E; Indicated in Position, direction of movement The number of newly generated photons carrying energy E; Indicated in Position, for photons with energy E, the linear attenuation coefficient of the material; Indicated in Position, the direction of movement is Photons carrying energy E' interact with matter and are scattered into new directions. The proportion of new photons carrying energy E; Emax It represents the highest energy of photons in the CT system, which is equivalent to the maximum photon energy emitted by the X-ray source.

[0058] Furthermore, the embodiments of the present application can obtain the photon flux integral expression and the photon flux differential expression of the photon flux and the photon source term in the static Boltzmann transport equation during multi-energy CT scanning, and the photon source term integral expression and the photon source term differential expression of the photon source term, and then determine the photon flux discrete expression of the photon flux, the photon source term discrete expression of the photon source term and the initial scattering signal intensity equation of the scattering signal, and use the discrete expression of the photon flux and the discrete expression of the photon source term to discretize the initial scattering signal intensity equation, and then obtain the final scattering signal intensity equation after discretization, thereby constructing the Boltzmann transport equation.

[0059] Exemplarily, the embodiment of the present application uses the iterative source method to decompose and decouple equation (1) based on the physical process in the CT scanning process. Considering the process of photon transport and scattering inside the object, the iterative source method is used to perform multi-order expansion of the photon flux and the photon source term according to the number of photon scatterings, and then obtain the photon flux integral expression and photon flux differential expression of the photon flux, the photon source term integral expression and photon source term differential expression of the photon source term, which can be expressed as but not limited to:

[0060]

[0061] Among them, φ (n) is the photon flux term consisting of photons scattered n times, S (n) is the photon scattering source term composed of photons scattered n times. This decouples the integral term and the differential term, and the following relationship exists between the photon flux and the photon source term at each order:

[0062]

[0063] In addition, it should be noted that in the embodiment of the present application, photons are emitted from the object and enter the detector to be recorded. In this stage, the photons are not scattered, and the photon flux φ received by the detector is D It can be expressed as, but not limited to:

[0064]

[0065] Among them, u and w are both integral variables. At this time, the energy deposition signal I received by the detector is scatter , that is, the initial scattered signal intensity equation of the scattered signal can be expressed as but not limited to:

[0066]

[0067] in, is the unit normal vector of the detector plane, D(E) is the energy response of the detector. In an ideal state, the detector can completely deposit energy for photons of any energy, that is, D(E) = 1.

[0068] Furthermore, the embodiment of the present application can select a suitable discretization scheme for the degenerate equation after the above integral-differential decoupling, and Discretize the object into N ordered voxels in space. In terms of energy, from 0 to E max Divide the continuous energy into G energy groups (E g ); In terms of solid angle, the 4π solid angle is discretized into M discrete space angles Then the discrete expression of photon flux and the discrete expression of photon source term are obtained, which can be expressed as but not limited to:

[0069]

[0070] Based on this, the embodiment of the present application can write equation (7) and equation (8) into matrix form, which can be expressed as but not limited to:

[0071] Ф (n) =TS (n) ,S (n-1) =CФ (n)

[0072] , (9)

[0073] After discretization, the final scattering signal intensity equation after discretization of the initial scattering signal intensity equation can be expressed as, but not limited to:

[0074]

[0075] Based on this, the embodiment of the present application can write formula (10) in matrix form, which can be expressed as but not limited to:

[0076]

[0077] Furthermore, the embodiment of the present application introduces at least one label dimension into the phase space (assuming that L labels are introduced, which is not specifically limited in the present application) based on the discrete expression of photon flux, the discrete expression of photon source term and the final scattered signal intensity equation through a discrete phase space solution method, and the phase space is transformed from Expand to The introduced label dimension has nothing to do with the actual physical process. Based on this, for different labels l, the matrices F, T, C, and D are exactly the same, that is, they can be expressed as but not limited to:

[0078]

[0079] Furthermore, the embodiment of the present application can simplify formula (12), which can be expressed as but not limited to:

[0080]

[0081] in, is a matrix independent of l.

[0082] Among them, in the actual calculation of the embodiment of the present application, it is only necessary to calculate T, C, and D once. The increase of the L dimension only increases the amount of calculation of 2K (L-1) sparse matrix multiplications, which is very small compared to the amount of calculation of T, C, and D. For example, the exponential operation contained in the solution of the T and D matrices, and the tracking of the length within the voxel, all of these calculations are much larger than the aforementioned sparse matrix multiplication. In other words, the introduction of the label dimension in the embodiment of the present application will only increase the amount of calculation slightly, but it makes it possible to calculate multiple scattering images under different energy spectra at the same time. For multiple X-ray energy spectra 1,2,…,L, remember Calculate F, T, C, and D to get Scattering images Among them, I=[I 1 ,I 2 ,…,I L ], I size is N D ×L,N D is the number of detector pixels, I l ,l=1,2,..,L, its size is N D ×1.

[0083] In particular, in the case where L≥G or it is uncertain whether a new energy spectrum scan will be used in the future, the embodiment of the present application can construct R G A set of bases in space, let’s take a set of base vectors as the unit orthogonal basis E = [E 1 ,E 2 ,…,E G ], E is a G×G diagonal matrix, so that a set of basis scattering images can be calculated Among them, I E Size is N D ×G. Therefore, for any energy spectrum The embodiment of the present application can easily obtain the corresponding scattering image In this case, since the photoelectric effect and Compton effect are taken into account in the equation, the energy of the photon will only gradually decrease and not increase. Assuming that energy group 1 is the highest energy group and energy group G is the lowest energy group, then for the base energy spectrum E i(That is, for the energy spectrum where only energy group i has values), only the photon behavior of the low-energy groups i, i+1, i+2…, G-1, G needs to be considered, without considering the behavior of the high-energy groups 1, 2,…, i-1. This enables the small computational amount introduced by increasing the label dimension to be further reduced by half, allowing the judgment condition L≥G proposed above to be further corrected to L>=G / 2.

[0084] The above analysis shows that in the embodiments of the present application, if L<G / 2, the original energy spectrum is directly discretized and labeled by energy groups, introducing L labels; if L>=G / 2, the above energy spectrum basis vector method is introduced, introducing G labels, thereby enabling the optimization of the calculation speed.

[0085] Thus, a Boltzmann transport equation with labels is constructed, which can be but is not limited to being expressed as:

[0086]

[0087] In step S103, the actual scanned energy spectrum during the multi-energy CT scan of the object is discretized and labeled, so as to determine the energy spectrum matrix in the pre-constructed Boltzmann transport equation with labels based on the discretized and labeled actual scanned energy spectrum.

[0088] As a possible implementation, the embodiments of the present application can discretize the actual scanned energy spectrum during the multi-energy CT scan of the object into multiple energy spectrum vectors, and then combine the label dimension to obtain the energy spectrum matrix.

[0089] Optionally, in an embodiment of the present application, discretizing and labeling the actual scanned energy spectrum during the multi-energy CT scan of the object to determine the energy spectrum matrix in the pre-constructed Boltzmann transport equation with labels based on the discretized and labeled actual scanned energy spectrum includes: determining whether the number of the discretized and labeled actual scanned energy spectra is greater than half of the number of preset discrete energy groups; if the number of the discretized and labeled actual scanned energy spectra is greater than half of the number of preset discrete energy groups, determining the dimension and value of the energy spectrum matrix in the pre-constructed Boltzmann transport equation with labels based on the number of preset discrete energy groups, where the energy spectrum matrix is an identity matrix; if the number of the discretized and labeled actual scanned energy spectra is less than or equal to half of the number of preset discrete energy groups, determining the dimension and value of the energy spectrum matrix in the pre-constructed Boltzmann transport equation with labels based on the number of the discretized and labeled actual scanned energy spectra, where each column vector in the energy spectrum matrix is the energy spectrum vector of each discretized and labeled actual scanned energy spectrum.

[0090] As a possible implementation method, the embodiment of the present application can determine the energy spectrum matrix in the pre-constructed labeled Boltzmann transport equation by judging whether the number of actual scanned energy spectra after discretization and labeling is greater than half of the number of certain discrete energy groups. Specifically, when the number of actual scanned energy spectra after discretization and labeling is greater than half of the number of preset discrete energy groups, the embodiment of the present application can determine the dimension and value of the energy spectrum matrix based on the number of certain discrete energy groups, wherein the energy spectrum matrix is ​​a unit matrix; when the number of actual scanned energy spectra after discretization and labeling is less than or equal to half of the number of preset discrete energy groups, the dimension and value of the energy spectrum matrix can be determined based on the number of actual scanned energy spectra after discretization and labeling, wherein each column vector in the energy spectrum matrix is ​​an energy spectrum vector after the actual scanned energy spectrum after each discretization and labeling is discretized. The number of certain discrete energy groups can be set by a person skilled in the art according to actual conditions, and this application does not impose specific restrictions.

[0091] For example, the embodiment of the present application inputs L different initial energy spectra S actually used in the multi-energy CT scanning process. l , and discretize the sampling into a certain number of scattering energy spectra, such as G certain scattering energy spectra Get the matrix If L<=G / 2, the energy spectrum matrix is ​​determined to be S, with a dimension of G×L; if L>G / 2, the energy spectrum matrix is ​​determined to be the unit matrix E, with a dimension of G×G.

[0092] In step S104, based on the final scattered signal intensity equation, energy spectrum matrix and material density distribution data in the pre-constructed labeled Boltzmann transport equation, the multi-energy scattered signal data of the object in the multi-energy CT scan is calculated. The expression of the final scattered signal intensity equation can be but is not limited to:

[0093]

[0094] It is understandable that the input of the labeled Boltzmann transport equation pre-constructed in the embodiment of the present application may include but is not limited to material density distribution data, geometric parameters and discretization parameters, etc., and the present application does not make specific restrictions. Among them, the geometric parameters may include but are not limited to the source-detection distance, detector size, detector bias, etc., and the present application does not make specific restrictions; the discretization parameters may include but are not limited to the number of voxels N, the number of discrete angles M, the number of energy groups G, etc., and the present application does not make specific restrictions.

[0095] In addition, it should be noted that in the embodiments of the present application, the calculation of the pre-constructed labeled Boltzmann transport equation is implemented by computer programming, and a GPU is used for parallel accelerated calculation, which can be specifically configured by technicians in this field according to actual conditions, and this application does not impose any specific restrictions.

[0096] In some embodiments, the embodiment of the present application can use the material density distribution data and the energy spectrum matrix as the input of the pre-constructed labeled Boltzmann transport equation, and then calculate the multi-energy scattering signal data of the object. Among them, the main content of the embodiment of the present application calculating the multi-energy scattering signal data of the object can be: first determine whether the number of actual scanning energy spectra of the object after discretization and labeling is greater than half of the number of certain discrete energy groups during the multi-energy CT scanning process, and when the number of actual scanning energy spectra after discretization and labeling is greater than half of the number of certain discrete energy groups, combine the energy spectrum matrix and the material density distribution data, and calculate the scattering signal intensity by the final scattering signal intensity equation in the pre-constructed labeled Boltzmann transport equation. The multi-energy scattering signal data can be obtained by multiplying the scattering signal intensity and the energy spectrum vector; and when the number of actual scanning energy spectra after discretization and labeling is less than or equal to half of the preset number of discrete energy groups, combine the energy spectrum matrix and the material density distribution data, and calculate the scattering signal intensity by the final scattering signal intensity equation in the pre-constructed labeled Boltzmann transport equation, that is, obtain the multi-energy scattering signal data;

[0097] For example, the embodiment of the present application inputs L different initial energy spectra actually used in the multi-energy CT scanning process. And discretize the sampling into a certain number of scattering energy spectra, such as G certain scattering energy spectra Get the matrix If L <= G / 2, the scattered signal intensity obtained according to the final scattered signal intensity equation can be expressed as but not limited to: Then the multi-energy scattering signal data can be expressed as but not limited to If L>G / 2, the scattered signal intensity obtained according to the final scattered signal intensity equation can be expressed as, but not limited to, I E , then the multi-energy scattering signal data can be expressed as but not limited to

[0098] The final scattered signal intensity equation may be expressed as, but not limited to,

[0099]

[0100] For example, the embodiment of the present application can obtain the attenuation coefficient distribution μ based on the material density distribution data. t and the linear scattering coefficient distribution μ s , and set L = 2, then input it into the Boltzmann transport equation implemented by the computer, and obtain the scattered signal intensity through a single calculation Then the multi-energy scattering signal data is

[0101] In step S105, scatter correction is performed on the initial multi-energy projection data using the multi-energy scatter signal data to obtain actual multi-energy projection data of the object after scatter correction.

[0102] As a possible implementation method, the embodiment of the present application can use the multi-energy scattering signal data to scatter-correct the initial multi-energy projection data, and then obtain the actual multi-energy projection data after scattering correction.

[0103] For example, in the present embodiment, a suitable scatter correction coefficient e may be selected first. 1 ,e 2 , and then obtain the actual multi-energy projection data after scatter correction The scatter correction process is completed. The calculation results are as follows Figure 2(a)-Figure 3 shown.

[0104] The working principle of the multi-energy CT scatter correction method based on the Boltzmann transport equation proposed in the embodiment of the present application is introduced below in conjunction with an embodiment.

[0105] Figure 4 The flowchart is a working principle of a multi-energy CT scatter correction method based on the Boltzmann transport equation according to one embodiment of the present application.

[0106] Step S401: data preprocessing.

[0107] In this embodiment of the present application, the initial multi-energy projection data (p 1 ,p 2 ,…,p n ) to perform rough scatter correction and hardening correction processing, and then obtain the pre-processed initial correction multi-energy projection data

[0108] Step S402: material density conversion.

[0109] In this embodiment of the present application, the initial correction multi-energy projection data can be Perform material density conversion, such as material decomposition, pre-reconstruction, etc., which is not specifically limited in this application, and then obtain the material density distribution data (ρ 1 ,ρ 2 ,…ρ n ).

[0110] Step S403: discretization of energy spectrum and labeling.

[0111] Among them, the embodiment of the present application can be based on actual scanning spectrum data The energy spectrum is discretized using the discretization parameters and labeled with the label dimension to obtain the energy spectrum matrix (E or S).

[0112] Step S404: Solving the pre-constructed labeled Boltzmann transport equation.

[0113] Among them, the embodiment of the present application can be based on material density distribution data (ρ 1 ,ρ 2 ,…ρ n ), energy spectrum matrix (E or S), geometric parameters and discretization parameters, using the pre-built labeled Boltzmann transport equation to solve the multi-energy scattering signal data

[0114] Step S405: scatter correction.

[0115] Among them, the embodiment of the present application can be based on the initial multi-energy projection data (p 1 ,p 2 ,…,p n ), multi-energy scattering signal data and scatter correction coefficients to perform scatter correction, and then obtain the actual multi-energy projection data after scatter correction

[0116] According to the multi-energy CT scattering correction method based on the Boltzmann transport equation proposed in the embodiment of the present application, the initial multi-energy projection data to be corrected obtained after the multi-energy CT scanning of the object can be preprocessed, and then the initial corrected multi-energy projection number after preprocessing can be obtained, and the material density conversion can be performed to obtain the material density distribution data of the object, and then the multi-energy energy spectrum data can be discretized and labeled to obtain the energy spectrum matrix, and then the material density distribution data and the energy spectrum matrix can be combined, and the multi-energy scattering signal data of the object in the multi-energy CT scanning can be calculated using the pre-constructed labeled Boltzmann transport equation, and the initial multi-energy projection data can be scatter-corrected using the multi-energy scattering signal data to obtain the actual multi-energy projection data after scattering correction, and the multi-energy scattering signal data can be calculated based on the pre-constructed labeled Boltzmann transport equation, and the equation can be solved in an analytical way through the discretization of the phase space, so that the program has excellent parallelism, thereby improving the calculation speed and portability. In order to meet the scatter correction requirements in multi-energy CT, an additional label dimension is introduced to realize the vectorization of the energy group space, thereby reducing the amount of calculation and achieving the goal of obtaining scatter signals under multiple energy spectra with a single calculation. This solves the problem in related technologies that, because material decomposition is very sensitive to the energy spectrum, some simple descattering methods can remove most of the scattering, but will cause serious distortion of the material density obtained by material decomposition, and because of the multi-energy spectrum, multiple scattering estimates need to be performed under different X-ray source conditions, and the processing process is lengthy and complicated.

[0117] Next, a multi-energy CT scatter correction device based on the Boltzmann transport equation proposed in an embodiment of the present application is described with reference to the accompanying drawings.

[0118] Figure 5 It is a block diagram of a multi-energy CT scatter correction device based on the Boltzmann transport equation provided according to an embodiment of the present application.

[0119] like Figure 5 As shown, the multi-energy CT scatter correction device 50 based on the Boltzmann transport equation includes: a preprocessing module 100 , a material density conversion module 200 , an energy spectrum discrete marking module 300 , a calculation module 400 and a scatter correction module 500 .

[0120] The preprocessing module 100 is used to preprocess the initial multi-energy projection data to be corrected obtained after the object is scanned by multi-energy CT, so as to obtain initial corrected multi-energy projection data after the preprocessing of the initial multi-energy projection data to be corrected.

[0121] The material density conversion module 200 is used to perform material density conversion on the initial corrected multi-energy projection data to obtain material density distribution data of the object.

[0122] The energy spectrum discrete labeling module 300 is used to discretize and label the actual scanning energy spectrum of the object during the multi-energy CT scanning process, so as to determine the energy spectrum matrix in the pre-constructed labeled Boltzmann transport equation based on the discretized and labeled actual scanning energy spectrum.

[0123] The calculation module 400 is used to calculate the multi-energy scattering signal data of the object in the multi-energy CT scanning based on the final scattering signal intensity equation, energy spectrum matrix and material density distribution data in the pre-constructed labeled Boltzmann transport equation.

[0124] The scatter correction module 500 is used to perform scatter correction on the initial multi-energy projection data using the multi-energy scatter signal data to obtain actual multi-energy projection data of the object after scatter correction.

[0125] Optionally, in one embodiment of the present application, it further includes: a first building module, a first acquisition module, a second acquisition module, a generation module, a correction module, a determination module and a second building module.

[0126] Among them, the first construction module is used to construct a static Boltzmann transport equation in the multi-energy CT scanning process before calculating the multi-energy scattering signal data of the object in the multi-energy CT scanning based on the final scattering signal intensity equation, energy spectrum matrix and material density distribution data in the pre-constructed labeled Boltzmann transport equation.

[0127] The first acquisition module is used to obtain the photon flux integral expression and photon flux differential expression of the photon flux and the photon source term in the static Boltzmann transport equation during the multi-energy CT scanning process, and the photon source term integral expression and photon source term differential expression of the photon flux.

[0128] The second acquisition module is used to acquire the photon flux discrete expression of the photon flux and the photon source term discrete expression of the photon source term based on the photon flux integral expression, the photon flux differential expression, the photon source term integral expression and the photon source term differential expression.

[0129] A generation module is used to generate an initial scattering signal intensity equation of a scattering signal during a multi-energy CT scanning process based on a photon flux integral expression, a photon flux differential expression, a photon source term integral expression, and a photon source term differential expression.

[0130] The correction module is used to discretize the initial scattering signal intensity equation by using the discrete expression of photon flux and the discrete expression of photon source term, so as to obtain the final scattering signal intensity equation after the initial scattering signal intensity equation is discretized.

[0131] A determination module is used to determine at least one label dimension of the static Boltzmann transport equation based on the discrete expression of the photon flux, the discrete expression of the photon source term and the final scattered signal intensity equation.

[0132] The second building module is used to build a labeled Boltzmann transport equation based on the discrete expression of photon flux, the discrete expression of photon source term, the final scattered signal intensity equation and at least one label dimension.

[0133] Optionally, in one embodiment of the present application, the energy spectrum discrete marking module 300 includes: a judgment unit, a first determination unit and a second determination unit.

[0134] The judging unit is used to judge whether the number of the actual scanned energy spectra after discretization and labeling is greater than half of the number of preset discrete energy groups.

[0135] The first determination unit is used to determine the dimension and value of the energy spectrum matrix in the pre-constructed labeled Boltzmann transport equation based on the preset number of discrete energy groups when the number of actual scanned energy spectra after discretization and labeling is greater than half of the preset number of discrete energy groups, wherein the energy spectrum matrix is ​​a unit matrix.

[0136] The second determination unit is used to determine the dimension and value of the energy spectrum matrix in the pre-constructed labeled Boltzmann transport equation based on the number of actual scanning energy spectra after discretization and labeling when the number of actual scanning energy spectra after discretization and labeling is less than or equal to half of the number of preset discrete energy groups, wherein each column vector in the energy spectrum matrix is ​​an energy spectrum vector of each actual scanning energy spectrum after discretization and labeling.

[0137] Optionally, in one embodiment of the present application, the expression of the pre-constructed labeled Boltzmann transport equation may be, but is not limited to,:

[0138]

[0139] in, is the angular flux of photons at voxel i, with energy group g, emission direction m, label l and scattering order n; is voxel i w′ The average linear attenuation coefficient for the energy group g at position i w′ is the voxel number where the w′ path is located; w′ is a vector The path represented passes through each voxel, and the path length in each voxel exists as an accumulated variable; The intensity of the photon source at voxel i, with energy group g, emission direction m, label l and scattering order n+1; μ s；i,g′→g,m′→m is the scattering cross section of a photon with a moving direction of m′ discrete angle and energy group g′ at voxel i; is the angular flux of photons at voxel i, with energy group g′, emission direction m′, label l and scattering order n, is the position of voxel i, the energy group is g, and the emission direction is m (the m direction is given by The direction of the photon source is determined by the label l and the scattering order is k, A d is the area of ​​pixel d, θ is The angle between the vector and the normal vector of pixel d, E g is the average energy of energy group g, D(E g ) is the detector at energy E g Lower detector gain, is the spatial position vector of voxel i, is the spatial position vector of the detector pixel d, i is the discrete voxel index, which identifies the discrete spatial position; g is the energy group index, which identifies the position of different discrete energy groups; m is the discrete solid angle index, which identifies different spatial directions; l is the label index, which is used to identify photons from different energy spectra; k is the scattering order, which exists as a summation variable in the formula; N is the total number of discretized voxels; G is the number of discretized energy groups, K is the highest scattering order selected for calculation; n is the scattering order (that is, the photon is scattered n times); n+1 is the scattering order (that is, the photon is scattered n+1 times).

[0140] Optionally, in one embodiment of the present application, the expression of the final scattered signal intensity equation may be, but is not limited to,:

[0141]

[0142] It should be noted that the above explanation of the embodiment of the multi-energy CT scatter correction method based on the Boltzmann transport equation is also applicable to the multi-energy CT scatter correction device based on the Boltzmann transport equation of this embodiment, which will not be repeated here.

[0143] According to the multi-energy CT scatter correction device based on the Boltzmann transport equation proposed in the embodiment of the present application, the initial multi-energy projection data to be corrected obtained after the multi-energy CT scanning of the object can be preprocessed, and then the initial corrected multi-energy projection number after preprocessing can be obtained, and the material density conversion can be performed to obtain the material density distribution data of the object, and then the multi-energy energy spectrum data can be discretized and labeled to obtain the energy spectrum matrix, and then the material density distribution data and the energy spectrum matrix can be combined, and the multi-energy scattering signal data of the object in the multi-energy CT scanning can be calculated using the pre-constructed labeled Boltzmann transport equation, and the initial multi-energy projection data can be scatter corrected using the multi-energy scattering signal data to obtain the actual multi-energy projection data after scattering correction, and the multi-energy scattering signal data can be calculated based on the pre-constructed labeled Boltzmann transport equation, and the equation can be solved in an analytical way through the discretization of the phase space, so that the program has excellent parallelism, thereby improving the calculation speed and portability. In order to meet the scatter correction requirements in multi-energy CT, an additional label dimension is introduced to realize the vectorization of the energy group space, thereby reducing the amount of calculation and achieving the goal of obtaining scatter signals under multiple energy spectra with a single calculation. This solves the problem in related technologies that, because material decomposition is very sensitive to the energy spectrum, some simple descattering methods can remove most of the scattering, but will cause serious distortion of the material density obtained by material decomposition, and because of the multi-energy spectrum, multiple scattering estimates need to be performed under different X-ray source conditions, and the processing process is lengthy and complicated.

[0144] Figure 6 This is a schematic diagram of the structure of an electronic device provided according to an embodiment of the present application. The electronic device may include:

[0145] A memory 601 , a processor 602 , and a computer program stored in the memory 601 and executable on the processor 602 .

[0146] When the processor 602 executes the program, the multi-energy CT scatter correction method based on the Boltzmann transport equation provided in the above embodiment is implemented.

[0147] Furthermore, the electronic device further comprises:

[0148] The communication interface 603 is used for communication between the memory 601 and the processor 602 .

[0149] The memory 601 is used to store computer programs that can be executed on the processor 602 .

[0150] The memory 601 may include a high-speed RAM memory, and may also include a non-volatile memory (non-volatile memory), such as at least one disk memory.

[0151] If the memory 601, the processor 602 and the communication interface 603 are implemented independently, the communication interface 603, the memory 601 and the processor 602 can be connected to each other through a bus and communicate with each other. The bus can be an Industry Standard Architecture (ISA) bus, a Peripheral Component Interconnect (PCI) bus or an Extended Industry Standard Architecture (EISA) bus. The bus can be divided into an address bus, a data bus, a control bus, etc. For ease of representation, Figure 6 Only one thick line is used in the diagram, but this does not mean that there is only one bus or only one type of bus.

[0152] Optionally, in a specific implementation, if the memory 601, the processor 602 and the communication interface 603 are integrated on a chip, the memory 601, the processor 602 and the communication interface 603 can communicate with each other through an internal interface.

[0153] The processor 602 may be a central processing unit (CPU), or an application specific integrated circuit (ASIC), or one or more integrated circuits configured to implement the embodiments of the present application.

[0154] The embodiment of the present application further provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the above-mentioned multi-energy CT scatter correction method based on the Boltzmann transport equation.

[0155] An embodiment of the present application further provides a computer program product, including a computer program, which, when executed, implements the above-mentioned multi-energy CT scatter correction method based on the Boltzmann transport equation.

[0156] In the description of this specification, the description with reference to the terms "one embodiment", "some embodiments", "example", "specific example", or "some examples" etc. means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present application. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described may be combined in any one or N embodiments or examples in a suitable manner. In addition, those skilled in the art may combine and combine the different embodiments or examples described in this specification and the features of the different embodiments or examples, without contradiction.

[0157] In addition, the terms "first" and "second" are used for descriptive purposes only and should not be understood as indicating or implying relative importance or implicitly indicating the number of technical features indicated. Therefore, a feature defined as "first" or "second" may explicitly or implicitly include at least one of the features. In the description of this application, "N" means at least two, such as two, three, etc., unless otherwise clearly and specifically defined.

[0158] Any process or method description in a flowchart or otherwise described herein may be understood to represent a module, fragment or portion of code comprising one or N executable instructions for implementing the steps of a custom logical function or process, and the scope of the preferred embodiments of the present application includes alternative implementations in which functions may not be performed in the order shown or discussed, including performing functions in a substantially simultaneous manner or in reverse order depending on the functions involved, which should be understood by technicians in the technical field to which the embodiments of the present application belong.

[0159] The logic and / or steps represented in the flowchart or otherwise described herein, for example, can be considered as an ordered list of executable instructions for implementing logical functions, and can be embodied in any computer-readable medium for use by an instruction execution system, device or apparatus (such as a computer-based system, a system including a processor, or other system that can fetch instructions from an instruction execution system, device or apparatus and execute instructions), or in combination with these instruction execution systems, devices or apparatuses. For the purpose of this specification, "computer-readable medium" can be any device that can contain, store, communicate, propagate or transmit a program for use by an instruction execution system, device or apparatus, or in combination with these instruction execution systems, devices or apparatuses. More specific examples of computer-readable media (a non-exhaustive list) include the following: an electrical connection with one or N wirings (electronic devices), a portable computer disk box (magnetic device), a random access memory (RAM), a read-only memory (ROM), an erasable and programmable read-only memory (EPROM or flash memory), a fiber optic device, and a portable compact disk read-only memory (CDROM). In addition, the computer-readable medium may even be paper or other suitable medium on which the program is printed, since the program may be obtained electronically by optically scanning the paper or other medium and then editing, interpreting or otherwise processing in a suitable manner if necessary and then storing it in a computer memory.

[0160] It should be understood that the various parts of the present application can be implemented by hardware, software, firmware or a combination thereof. In the above embodiment, the N steps or methods can be implemented by software or firmware stored in a memory and executed by a suitable instruction execution system. If implemented by hardware, as in another embodiment, it can be implemented by any one or a combination of multiple of the following technologies known in the art: a discrete logic circuit having a logic gate circuit for implementing a logic function for a data signal, a dedicated integrated circuit having a suitable combination of logic gate circuits, a programmable gate array (PGA), a field programmable gate array (FPGA), etc.

[0161] A person skilled in the art may understand that all or part of the steps in the above-mentioned embodiment method may be completed by instructing related hardware through a program, and the program may be stored in a computer-readable storage medium, which, when executed, includes one or a combination of the steps of the method embodiment.

[0162] In addition, each functional unit in each embodiment of the present application may be integrated into a processing module, or each unit may exist physically separately, or two or more units may be integrated into one module. The above-mentioned integrated module may be implemented in the form of hardware or in the form of a software functional module. If the integrated module is implemented in the form of a software functional module and sold or used as an independent product, it may also be stored in a computer-readable storage medium.

[0163] The storage medium mentioned above may be a read-only memory, a magnetic disk or an optical disk, etc. Although the embodiments of the present application have been shown and described above, it can be understood that the above embodiments are exemplary and cannot be understood as limiting the present application. A person of ordinary skill in the art may change, modify, replace and modify the above embodiments within the scope of the present application.

Claims

1. A multi-energy CT scatter correction method based on the Boltzmann transport equation, characterized in that: The following steps are involved: Preprocessing the initial multi-energy projection data to be corrected obtained after the object is scanned by multi-energy computed tomography (CT) to obtain initial corrected multi-energy projection data after the preprocessing of the initial multi-energy projection data to be corrected; Performing material density conversion on the initial corrected multi-energy projection data to obtain material density distribution data after material density conversion of the initial corrected multi-energy projection data; Discretizing and labeling the actual scanning energy spectrum of the object during the multi-energy CT scanning process, so as to determine the energy spectrum matrix in the pre-constructed labeled Boltzmann transport equation based on the discretized and labeled actual scanning energy spectrum; Calculating multi-energy scattering signal data of the object in multi-energy CT scanning based on the final scattering signal intensity equation in the pre-constructed labeled Boltzmann transport equation, the energy spectrum matrix and the material density distribution data; The multi-energy scatter signal data is used to perform scatter correction on the initial multi-energy projection data to obtain actual multi-energy projection data of the object after scatter correction.

2. The method according to claim 1, characterized in that Before calculating the multi-energy scattering signal data of the object in the multi-energy CT scan based on the final scattering signal intensity equation in the pre-constructed labeled Boltzmann transport equation, the energy spectrum matrix and the material density distribution data, the method further includes: Constructing a static Boltzmann transport equation of the object during multi-energy CT scanning; Obtaining the photon flux integral expression and the photon flux differential expression of the photon flux and the photon source term in the static Boltzmann transport equation during the multi-energy CT scanning, and the photon source term integral expression and the photon source term differential expression of the photon source term; Based on the photon flux integral expression, the photon flux differential expression, the photon source term integral expression and the photon source term differential expression, obtaining a photon flux discrete expression of the photon flux and a photon source term discrete expression of the photon source term; Based on the photon flux integral expression, the photon flux differential expression, the photon source term integral expression and the photon source term differential expression, generating an initial scattering signal intensity equation of the scattering signal during the multi-energy CT scanning process; Discretizing the initial scattered signal intensity equation by using the discrete expression of the photon flux and the discrete expression of the photon source term to obtain a final scattered signal intensity equation after the initial scattered signal intensity equation is discretized; Determining at least one label dimension of the static Boltzmann transport equation based on the discrete expression of the photon flux, the discrete expression of the photon source term, and the final scattered signal intensity equation; A labeled Boltzmann transport equation is constructed based on the discrete expression of the photon flux, the discrete expression of the photon source term, the final scattered signal intensity equation and the at least one label dimension.

3. The method according to claim 1, characterized in that The actual scanning energy spectrum of the object during the multi-energy CT scanning process is discretized and labeled to determine the energy spectrum matrix in the pre-constructed labeled Boltzmann transport equation based on the discretized and labeled actual scanning energy spectrum, including: Determining whether the number of actually scanned energy spectra after the discretization and labeling is greater than half of the number of preset discrete energy groups; If the number of the actual scanned energy spectra after the discretization and labeling is greater than half of the number of the preset discrete energy groups, the dimension and value of the energy spectrum matrix in the pre-constructed labeled Boltzmann transport equation are determined based on the number of the preset discrete energy groups, wherein the energy spectrum matrix is ​​a unit matrix; If the number of actual scanning energy spectra after discretization and labeling is less than or equal to half of the number of preset discrete energy groups, the dimension and value of the energy spectrum matrix in the pre-constructed labeled Boltzmann transport equation are determined based on the number of actual scanning energy spectra after discretization and labeling, wherein each column vector in the energy spectrum matrix is ​​an energy spectrum vector of each actual scanning energy spectrum after discretization and labeling.

4. The method according to claim 1, characterized in that: The expression of the pre-built labeled Boltzmann transport equation is: in, is the angular flux of photons at voxel i, with energy group g, emission direction m, label l and scattering order n; is voxel i w′ The average linear attenuation coefficient for the energy group g at position i w′ w ′ The voxel number where the path is located; w ′ For vector The path represented passes through each voxel, and the path length in each voxel exists as an accumulated variable; The intensity of the photon source at voxel i, with energy group g, emission direction m, label l and scattering order n+1; μ s；i,g′→g,m′→m is the position of voxel i, and the moving direction is m ′ Discrete angle and energy group is g ′ The photon scattering is the scattering cross section of photons with a moving direction of m discrete angle and energy group g; is the position of voxel i, and the energy group is g ′ , the emission direction is m ′ , the angular flux of photons with label l and scattering order n, is the position of voxel i, the energy group is g, and the emission direction is m (the m direction is given by The direction of the photon source is determined by the label l and the scattering order is k, A d is the area of ​​pixel d, θ is The angle between the vector and the normal vector of pixel d, E g is the average energy of energy group g, D(E g ) is the detector at energy E g Lower detector gain, is the spatial position vector of voxel i, is the spatial position vector of the detector pixel d, i is the discrete voxel index, which identifies the discrete spatial position; g is the energy group index, which identifies the position of different discrete energy groups; m is the discrete solid angle index, which identifies different spatial directions; l is the label index, which is used to identify photons from different energy spectra; k is the scattering order, which exists as a summation variable in the formula; N is the total number of discretized voxels; G is the number of discretized energy groups, K is the highest scattering order selected for calculation; n is the scattering order (that is, the photon is scattered n times); n+1 is the scattering order (that is, the photon is scattered n+1 times).

5. The method according to claim 1, characterized in that The final scattered signal intensity equation is expressed as:

6. A multi-energy CT scatter correction device based on the Boltzmann transport equation, characterized in that: include: A preprocessing module, used for preprocessing the initial multi-energy projection data to be corrected obtained after the object is scanned by multi-energy CT, so as to obtain initial corrected multi-energy projection data after the preprocessing of the initial multi-energy projection data to be corrected; A material density conversion module, used for performing material density conversion on the initial corrected multi-energy projection data to obtain material density distribution data after material density conversion of the initial corrected multi-energy projection data; An energy spectrum discrete labeling module, used for discretizing and labeling the actual scanning energy spectrum of the object during the multi-energy CT scanning process, so as to determine the energy spectrum matrix in the pre-constructed labeled Boltzmann transport equation based on the discretized and labeled actual scanning energy spectrum; A calculation module, configured to calculate multi-energy scattering signal data of the object in a multi-energy CT scan based on a final scattering signal intensity equation in the pre-constructed labeled Boltzmann transport equation, the energy spectrum matrix and the material density distribution data; The scatter correction module is used to perform scatter correction on the initial multi-energy projection data using the multi-energy scatter signal data to obtain actual multi-energy projection data of the object after scatter correction.

7. The device according to claim 6, characterized in that Also includes: A first construction module is used to construct a static Boltzmann transport equation in a multi-energy CT scanning process before calculating the multi-energy scattering signal data of the object in the multi-energy CT scanning based on the final scattering signal intensity equation in the pre-constructed labeled Boltzmann transport equation, the energy spectrum matrix and the material density distribution data; A first acquisition module is used to acquire the photon flux and the photon source term in the multi-energy CT scanning process, the photon flux integral expression and the photon flux differential expression of the photon flux, the photon source term integral expression and the photon source term differential expression of the photon source term in the static Boltzmann transport equation; A second acquisition module is used to acquire a photon flux discrete expression of the photon flux and a photon source term discrete expression of the photon source term based on the photon flux integral expression, the photon flux differential expression, the photon source term integral expression and the photon source term differential expression; A generating module, configured to generate an initial scattering signal intensity equation of the scattering signal during the multi-energy CT scanning process based on the photon flux integral expression, the photon flux differential expression, the photon source term integral expression and the photon source term differential expression; A correction module, used for discretizing the initial scattered signal intensity equation by using the discrete expression of the photon flux and the discrete expression of the photon source term, so as to obtain a final scattered signal intensity equation after the initial scattered signal intensity equation is discretized; A determination module, configured to determine at least one label dimension of the static Boltzmann transport equation based on the discrete expression of the photon flux, the discrete expression of the photon source term and the final scattered signal intensity equation; The second construction module is used to construct the labeled Boltzmann transport equation based on the discrete expression of the photon flux, the discrete expression of the photon source term, the final scattered signal intensity equation and the at least one label dimension.

8. An electronic device, characterized in that: include: A memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the multi-energy CT scatter correction method based on the Boltzmann transport equation as claimed in any one of claims 1 to 5.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that: The program is executed by a processor to implement the multi-energy CT scatter correction method based on the Boltzmann transport equation as described in any one of claims 1 to 5.

10. A computer program product, characterized in that It comprises a computer program, which, when executed, is used to implement the multi-energy CT scatter correction method based on the Boltzmann transport equation as claimed in any one of claims 1 to 5.

Images