A three-dimensional synthetic imaging method and device based on coded energy spectrum

Abstract

The present invention relates to a three-dimensional synthetic imaging method and device based on coded energy spectrum, including using an X-ray source to scan a phantom, and using a detector to detect the X-ray beam attenuated by the phantom, so as to reconstruct a three-dimensional digital synthetic image, wherein a K-edge coded aperture is fixed before the X-ray source, and each pixel in the K-edge coded aperture is made of K-edge material; the reconstruction process of the three-dimensional digital synthetic image includes: constructing a measurement value calculation formula adding the K-edge coded aperture, and dividing the energy spectrum interval to obtain a discrete measurement value calculation formula, and finally constructing a comprehensive measurement value calculation formula obtained by the X-ray source at each projection viewing angle, substituting the initial three-dimensional image, and comparing it with the actual measurement value obtained, performing iterative optimization, and obtaining the optimal three-dimensional image as the reconstruction result of the three-dimensional tomographic synthesis image. Compared with the prior art, the present invention can effectively perceive the energy spectrum attenuation information of an object, discover tiny lesions hidden in a two-dimensional image, and reduce unnecessary biopsies.

Description

Technical Field

[0001] The present invention relates to the field of three-dimensional tomosynthesis imaging, and in particular to a three-dimensional synthesis imaging method and device based on coded energy spectrum. Background Art

[0002] 3D tomosynthesis is an emerging tomographic imaging method that acquires a few 2D projections of an object from a limited range of angles to reconstruct tomographic images of different depths reflecting the 3D information of the object. Compared with traditional computed tomography technology, it has the advantages of low X-ray radiation dose and low economic cost. As a promising screening tool, it has been used clinically for the detection of lung cancer and breast cancer lesions.

[0003] A typical 3D tomosynthesis camera system consists of a high-voltage X-ray source, a rotating gantry, and an energy-integrating detector. It collects the X-ray intensity projected across the entire energy spectrum as measurement data, and then uses the negative logarithm of the measurement data for filtered back-projection reconstruction. However, such a system cannot obtain measurements containing energy spectrum attenuation information, which can be used to identify the material composition of an object and reduce beam hardening artifacts. An existing prototype of a spectral 3D tomosynthesis camera system uses a photon counting detector, which allows high-energy and low-energy measurements to be obtained simultaneously. However, it is currently expensive and suffers from pulse pile-up at high X-ray fluxes, which limits its use in routine clinical examinations. Summary of the invention

[0004] The purpose of the present invention is to provide a three-dimensional synthetic imaging method and device based on coded energy spectrum in order to overcome the defect of the above-mentioned prior art that the measurement value containing energy spectrum attenuation information cannot be obtained.

[0005] The purpose of the present invention can be achieved by the following technical solutions:

[0006] A three-dimensional synthetic imaging method based on coded energy spectrum, comprising: using an X-ray source to scan a phantom, and using a detector to detect the X-ray beam attenuated by the phantom, thereby reconstructing a three-dimensional digital synthetic image, wherein a K-edge coded aperture is fixed in front of the X-ray source, and each pixel in the K-edge coded aperture is made of a K-edge material, and the K-edge coded aperture includes a plurality of K-edge materials;

[0007] The reconstruction process of the three-dimensional digital composite image includes:

[0008] Construct a calculation formula for the continuous measurement value of the initial X-ray passing through each pixel of the K-edge coded aperture through the phantom to the detector;

[0009] According to the number of types of K-edge materials in the K-edge coded aperture, the continuous measurement value calculation formula is divided into energy spectrum intervals to obtain a discrete measurement value calculation formula;

[0010] According to the discrete measurement value calculation formula, a comprehensive measurement value calculation formula obtained by the X-ray source at each projection viewing angle is constructed, substituted into the initial three-dimensional image, and compared with the actual measurement values ​​obtained, iterative optimization is performed to obtain the optimal three-dimensional image as the three-dimensional tomography synthesis image reconstruction result.

[0011] Furthermore, the continuous measurement value calculation formula is:

[0012] I=∫f(E)I(E)dE

[0013]

[0014]

[0015] Where I is the measured value of the detector, E is the X-ray photon energy, f(E) is the transmittance of each pixel in the K-edge coded aperture, I(E) is the intensity attenuation value of the initial X-ray passing through the phantom, and t f is the thickness of the K-edge material on one pixel in the K-edge coding aperture, μ f (E) is the linear attenuation coefficient of the K-edge material on a pixel in the K-edge coded aperture, ψ is the projection angle, D(E) is the detector energy spectrum response function, I0(E) is the initial intensity of the X-ray, is the linear attenuation coefficient of the model, is the projection path of the X-ray beam reaching the detector.

[0016] Furthermore, the discrete measurement value calculation formula is:

[0017]

[0018] Where K = F + 1, F is the number of types of K-edge materials in the K-edge coding aperture, E k is the Kth X-ray photon energy interval.

[0019] Furthermore, the process of obtaining the comprehensive measurement value calculation formula specifically includes:

[0020] The component information of the model to be scanned is loaded into the discrete measurement value calculation formula, and then the measurement value calculation formula on the jth projection path at the pth projection viewing angle is constructed. Finally, by transforming it into a matrix form, a measurement value calculation formula obtained from P projection viewing angles is constructed, which is the comprehensive measurement value calculation formula.

[0021] Furthermore, the component information of the phantom to be scanned includes fat and glands, and the calculation expression of the linear attenuation coefficient of the phantom is:

[0022]

[0023] Where τ(E) represents the effective attenuation coefficient of the component, Indicates that the component is at point The specific gravity of the subscripts 1 and 2 correspond to the components fat and glandular, respectively;

[0024] The calculation formula for the measured value on the jth projection path at the pth projection viewing angle is:

[0025]

[0026] In the formula, is the measured value on the jth projection path at the pth projection angle, H p is the projection matrix of the p-th view, is the intersection of the jth projection path and the i-th pixel in the phantom, ρ i，l is the proportion of the lth material in the ith element, is the effective attenuation coefficient of the lth component in the kth energy spectrum interval, and N is the number of pixels into which the object to be reconstructed is discretized.

[0027] Furthermore, the matrix form is:

[0028] y p =(W⊙exp(-H p XT))·1 K

[0029] In the formula, To use the size of M x ×M y The measured value obtained by the detector is is the ray intensity in each energy spectrum interval without the model, where W j，k =f j (E k )D(E k )I0(E K ), Each column represents the weight of the lth basic material. represents the spectrum attenuation of L materials in a total of K energy spectrum intervals, ⊙ represents point-by-point multiplication; 1 K Represents a vector of length K with all values ​​1;

[0030] The calculation formula for the measurement value obtained from P projection viewing angles is:

[0031]

[0032] In the formula, represents the Kronecker product.

[0033] Furthermore, the iterative optimization process is specifically as follows:

[0034] Introduce three-dimensional total variation regularization and construct the optimization problem of three-dimensional image X;

[0035] Introduce the intermediate variable Z, substitute it into the optimization problem of the three-dimensional image X, and update X and Z alternately until the preset stop condition is met. The X obtained at this time is used as X * The three-dimensional tomosynthesis image is obtained by X * The product of the effective attenuation coefficient T of the base material is obtained.

[0036] Furthermore, the optimization problem of the three-dimensional image X is expressed as:

[0037]

[0038] In the formula, is the three-dimensional tensor form of X, λ is the weight parameter used to weigh the regularization term and the data term, is the calculation formula of the comprehensive measurement value of the three-dimensional image X, and y is the actual detection value of the detector;

[0039] The expression of the optimization problem after introducing the intermediate variable Z is:

[0040]

[0041] Furthermore, the Gauss-Newton method is used to update the three-dimensional image X, and the iterative shearing algorithm is used to update the intermediate variable Z.

[0042] The present invention also provides a three-dimensional synthetic camera device based on coded energy spectrum, comprising a memory and a processor, wherein the memory stores a computer program, and the processor calls the computer program to execute the steps of the method described above.

[0043] Compared with the prior art, the present invention has the following advantages:

[0044] (1) The present invention proposes a three-dimensional synthetic imaging method based on coded energy spectrum, which uses a coded aperture (K-edge coded aperture) made of multiple K-edge materials to modulate the complex cone beam X-ray energy spectrum emitted by high-voltage X-rays, wherein pixels of different K-edge materials in the K-edge coded aperture have different modulation effects; after the energy spectrum is modulated, the X-ray beam can effectively perceive the energy spectrum attenuation information of the object during the measurement process, so that the energy spectrum tomographic image of the object can be reconstructed from the obtained measurement data.

[0045] (2) The method of the present invention is used for three-dimensional synthesis of the breast, and can obtain three-dimensional tomographic images that reflect the structure and composition information of the breast. These images can be used to detect tiny lesions hidden in two-dimensional images; eliminate false positives and false negatives caused by two-dimensional projection images; and use the structure and composition information of the mass reflected in the energy spectrum three-dimensional tomographic images to judge disease information and reduce unnecessary biopsies.

[0046] (3) The present invention provides a corresponding reconstruction algorithm. First, the transmittance of each pixel in the K-edge coded aperture is added to the X-ray measurement value representation. Taking into account the characteristic that the linear attenuation coefficient of the K-edge material increases sharply, the entire X-ray energy spectrum is divided into multiple intervals. The measurement value representation is realized in a discretized form, which is conducive to improving the matching accuracy in the subsequent iterative process. Then, the measurement value representation is converted into a matrix form to realize the comprehensive representation of each projection viewing angle and projection path, and finally achieve comparison optimization with the actual detection value.

[0047] (4) In the iterative optimization process of obtaining the three-dimensional image X, the present invention adds a three-dimensional total variation (3DTV) regularization to improve the quality of the reconstructed image, introduces an intermediate variable Z, and iteratively updates X and Z, thereby further improving the matching accuracy of the obtained three-dimensional image X and the convergence speed in the iterative process. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 This is the energy spectrum attenuation curve of silver element, which has a steep increase at its K-edge energy of 25.49KeV;

[0049] Figure 2 The energy spectrum curve of the 50kV high voltage X-ray tube and the spectrum lines attenuated by the filter composed of the elements Ag and Pd are used to extract quasi-monoenergetic X-rays from a wide energy spectrum X-ray beam;

[0050] Figure 3 The schematic diagram of the coded spectrum digital breast 3D tomosynthesis system. The right side shows the K-edge coded aperture, where different color blocks represent pixels made of different K-edge materials and the corresponding filtered X-ray beams.

[0051] Figure 4 This is a 2D projection of the phantom used in the experiment at an angle of 0°, which contains two materials: fat and glands (soft tissue);

[0052] Figure 5 Energy spectrum transmission curves of the six K-edge materials used in the experiment to make K-edge coded apertures, including selenium, strontium, niobium, ruthenium, silver, and indium;

[0053] Figure 6 It is a three-dimensional tomosynthesis image of the two basic materials contained in the phantom, where Z = 2 to 10 represents the tomographic image of the even-numbered layers of the phantom; for ease of comparison, Figure 6 The true value corresponding to each image and the reconstructed image of the traditional digital breast 3D tomosynthesis without the K-edge coded aperture are also provided. It can be seen that the proposed coded energy spectrum digital breast 3D tomosynthesis photography technology can effectively realize material decomposition breast 3D tomosynthesis photography.

[0054] In the figure, attenuation is the attenuation value, energy is the energy, intensity is the intensity, original is the original ray spectrum, Filtered by Ag is the ray spectrum after filtering by adding silver material, Filtered by Pd is the ray spectrum after filtering by adding palladium material, Monoenergetic is the monoenergetic spectrum, and Spectrum is the energy spectrum. DETAILED DESCRIPTION

[0055] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Generally, the components of the embodiments of the present invention described and shown in the drawings here can be arranged and designed in various different configurations.

[0056] Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the invention claimed for protection, but merely represents selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0057] It should be noted that similar reference numerals and letters denote similar items in the following drawings, and therefore, once an item is defined in one drawing, further definition and explanation thereof is not required in subsequent drawings.

[0058] Example 1

[0059] This embodiment provides a three-dimensional synthetic imaging method based on coded energy spectrum, which mainly relies on the K-edge coded aperture described below, which is a pixelated coded aperture composed of multiple K-edge materials. It is first described below.

[0060] 1. K-edge coded aperture

[0061] As the energy of X-rays increases, the attenuation coefficient of matter to X-rays gradually decreases, but it resonates and absorbs X-rays of certain specific photon energies, causing the corresponding X-ray attenuation coefficient to increase sharply. This attenuation is caused by the K-layer electrons of the atom, so it is called the K-edge attenuation characteristic, and the corresponding energy is called the K-edge energy. We call this material a K-edge material. Different K-edge materials have different K-edge energies, and therefore have different energy spectrum modulation effects on X-rays. Figure 1 The energy spectrum attenuation curve of the silver element is shown, which has a steep increase at its K-edge energy of 25.49KeV; materials with adjacent atomic numbers have similar K-edge energies, which can be used to extract quasi-monoenergetic X-rays from a wide energy spectrum X-ray beam. A variety of K-edge materials can be used to make a multi-spectrum modulated K-edge coded aperture that the present invention relies on, in which each pixel is made of a certain thickness of K-edge material to modulate the energy spectrum of the projected X-ray beam. Figure 2 As shown, it is combined with the traditional digital breast 3D synthetic ray, which is the 3D synthetic imaging method based on coded energy spectrum proposed by the present invention, in which the K-edge coded aperture is fixed in front of the X-ray source and rotates synchronously with it during the scanning process without additional movement. The cone beam X-ray emitted by the ray source is modulated by the K-edge coded aperture and then illuminates the phantom. The projected X-ray beam carries the energy spectrum attenuation information of the phantom to the detector, and the measurement data containing the energy spectrum attenuation information can be obtained.

[0062] Since the above projection process is nonlinear, it is impossible to directly reconstruct the multi-spectrum three-dimensional digital composite image using the existing algorithm. Therefore, this embodiment also proposes a corresponding reconstruction algorithm, and its technical solution is as follows.

[0063] 2. 3D Reconstruction Algorithm

[0064] According to the Beer-Lambert law, when a monoenergetic X-ray with an initial intensity of I0(E) passes through the phantom, its intensity is attenuated to Where D(E) is the detector energy spectrum response function, is the linear attenuation coefficient of the object, is the projection path of the ray beam reaching the detector. Figure 2 In the coded energy spectrum digital breast 3D synthesis camera system shown in FIG, a K-edge coded aperture is used to modulate the energy spectrum of the X-ray beam before the X-ray source. The transmittance of one pixel can be expressed as Among them, t f and μ f(E) are the thickness and linear attenuation coefficient of the K-edge material on the pixel, and ψ is the projection angle. Therefore, the measurement value obtained by the integrating detector can be expressed as I = ∫f(E)I(E)dE.

[0065] As mentioned above, the linear attenuation coefficient of a K-edge material is f There is a sudden increase. When F kinds of K-edge filter materials are used for encoding, the entire X-ray energy spectrum is divided into K = F + 1 intervals: [E min , E1], [E1, E2], …, [E F , E max ], where E min and E max are the minimum and maximum X-ray source photon energies, respectively. Since the energy spectrum defined by similar atomic numbers is very narrow, the linear attenuation coefficient of the phantom and K-edge material in the kth interval is and μ f (E) can be respectively represented by their effective attenuation coefficients and μ f (E k ) approximation. Thus, the measured value I can be expressed as:

[0066]

[0067] Material decomposition is one of the important applications of spectral X-ray tomography. For 3D synthetic tomography used for breast scanning, the scanned phantom components usually include fat and glands, so its effective attenuation coefficient can be divided into a linear combination of these two components:

[0068]

[0069] Where τ(E) represents the effective attenuation coefficient of the component, Indicates that the component is at point In order to reconstruct using the inverse imaging process, this embodiment discretizes the object into N pixels, so that the measurement on the jth projection path at the pth projection viewing angle can be expressed as:

[0070]

[0071] Among them, H p is the projection matrix of the pth view, where the elements represents the intersection of the jth projection path and the i-th pixel in the phantom, ρ i，l represents the proportion of the lth material in the i-th element, represents the effective attenuation coefficient of the lth component in the kth energy spectrum interval. For simplicity, this embodiment can further rewrite (3) into a matrix form:

[0072] y p =(W⊙exp(-H p XT))·1 K (4)

[0073] in, Indicates that the utilization size is M x ×M y The measurement value obtained by the detector; Represents the ray intensity of each energy spectrum interval in the empty space without an object, where W j，k =f j (E k )D(E k )I0(E K ); Each column represents the weight of the lth basic material. represents the spectrum attenuation of L materials in a total of K energy spectrum intervals; ⊙ represents point-by-point multiplication; 1 K represents a vector of length K and all values ​​are 1. Further, the measurements obtained from a total of P projection viewpoints can be expressed as:

[0074]

[0075] in represents the Kronecker product.

[0076] According to the above forward projection model, the reconstruction of the tomographic image X reflecting the basic material information in the phantom can be expressed as an unconstrained optimization problem Considering the similarity between multi-layer tomographic images, the present invention introduces a three-dimensional total variation (3DTV) regularization to each basic material tomographic image to improve the quality of the reconstructed image. It is defined as:

[0077]

[0078] in Used to adapt to the resolution in different modes. Using formula (6) as the regularization term for reconstruction, the optimization problem can be obtained:

[0079]

[0080] in, Represents the three-dimensional tensor form of X, and λ is used to weigh the regularization term and the data term. In order to solve the reconstruction problem, an intermediate variable Z is introduced, and X and Z are updated alternately and iteratively:

[0081]

[0082] Where t represents the number of iterations. The Gauss-Newton method can be used to update X, while the iterative shearing algorithm can be used to update Z. Repeat multiple iterations until the preset stop condition (number of iterations, relative error, or error reduction) is met, and the resulting X is X * . The energy spectrum tomography image can be obtained by X * The product of the effective attenuation coefficient T of the base material is obtained.

[0083] This embodiment also provides a three-dimensional synthesis camera device based on coded energy spectrum, including a memory and a processor, wherein the memory stores a computer program, and the processor calls the computer program to execute the steps of the three-dimensional synthesis camera method based on coded energy spectrum as described above.

[0084] The preferred specific embodiments of the present invention are described in detail above. It should be understood that a person skilled in the art can make many modifications and changes based on the concept of the present invention without creative work. Therefore, any technical solution that can be obtained by a person skilled in the art through logical analysis, reasoning or limited experiments based on the concept of the present invention on the basis of the prior art should be within the scope of protection determined by the claims.

Claims

1. A three-dimensional synthetic imaging method based on coded energy spectrum, comprising: An X-ray source is used to scan a phantom, and a detector detects the X-ray beam attenuated by the phantom, thereby reconstructing a three-dimensional digital synthetic image, characterized in that a K-edge coded aperture is fixed in front of the X-ray source, each pixel in the K-edge coded aperture is made of a K-edge material, and the K-edge coded aperture includes a variety of K-edge materials; The reconstruction process of the three-dimensional digital composite image includes: Construct a calculation formula for the continuous measurement value of the initial X-ray passing through each pixel of the K-edge coded aperture through the phantom to the detector; According to the number of types of K-edge materials in the K-edge coded aperture, the continuous measurement value calculation formula is divided into energy spectrum intervals to obtain a discrete measurement value calculation formula; According to the discrete measurement value calculation formula, a comprehensive measurement value calculation formula obtained by the X-ray source at each projection viewing angle is constructed, substituted into the initial three-dimensional image, and compared with the actual measurement value, iterative optimization is performed to obtain the optimal three-dimensional image as the three-dimensional tomosynthesis image reconstruction result; The process of obtaining the comprehensive measurement value calculation formula specifically includes: Loading the component information of the phantom to be scanned into the discrete measurement value calculation formula, then constructing the measurement value calculation formula on the j-th projection path at the p-th projection viewing angle, and finally constructing the measurement value calculation formula obtained from P projection viewing angles by transforming into a matrix form, which is the comprehensive measurement value calculation formula; The component information of the phantom to be scanned includes fat and glands, and the calculation expression of the linear attenuation coefficient of the phantom is: Where τ(E) represents the effective attenuation coefficient of the component, Indicates that the component is at point The specific gravity of the subscripts 1 and 2 correspond to the components fat and glandular, respectively; The calculation formula for the measured value on the jth projection path at the pth projection viewing angle is: In the formula, is the measured value on the jth projection path at the pth projection angle, H p is the projection matrix of the p-th view, is the intersection of the jth projection path and the i-th pixel in the phantom, ρ i,l is the proportion of the lth material in the ith element, is the effective attenuation coefficient of the lth component in the kth energy spectrum interval, and N is the number of pixels into which the model to be reconstructed is discretized; The matrix form is: the p =(W⊙exp(-H p XT))·1 L In the formula, To use the size of M x ×M y The measured value obtained by the detector is is the ray intensity in each energy spectrum interval without the phantom, where Each column represents the weight of the lth basic material. represents the spectrum attenuation of L materials in a total of K energy spectrum intervals, ⊙ represents point-by-point multiplication; 1 K Represents a vector with a length of K and all element values ​​1; The calculation formula for the measurement value obtained from P projection viewing angles is: In the formula, represents the Kronecker product.

2. The three-dimensional synthetic imaging method based on coded energy spectrum according to claim 1, characterized in that: The continuous measurement value calculation formula is: I=∫f(E)I(E)dE Where I is the measured value of the detector, E is the energy of the X-ray photon, f(S) is the transmittance of each pixel in the K-edge coded aperture, I(E) is the intensity attenuation value of the initial X-ray passing through the phantom, and t f is the thickness of the K-edge material on one pixel in the K-edge coding aperture, μ f (E) is the linear attenuation coefficient of the K-edge material on a pixel in the K-edge coded aperture, ψ is the projection angle, D(E) is the detector energy spectrum response function, I0(E) is the initial intensity of the X-ray, is the linear attenuation coefficient of the model, is the projection path of the X-ray beam reaching the detector.

3. The three-dimensional synthetic imaging method based on coded energy spectrum according to claim 2, characterized in that: The discrete measurement value calculation formula is: Where K = F + 1, F is the number of types of K-edge materials in the K-edge coding aperture, E k is the Kth X-ray photon energy interval.

4. The three-dimensional synthetic imaging method based on coded energy spectrum according to claim 1, characterized in that: The iterative optimization process is specifically as follows: Introduce three-dimensional total variation regularization and construct the optimization problem of three-dimensional image X; Introduce the intermediate variable Z, substitute it into the optimization problem of the three-dimensional image X, and update X and Z alternately until the preset stop condition is met. The X obtained at this time is used as X * The three-dimensional tomosynthesis image is obtained by X * The product of the effective attenuation coefficient T of the base material is obtained.

5. The three-dimensional synthetic imaging method based on coded energy spectrum according to claim 4, characterized in that: The expression of the optimization problem of the three-dimensional image X is: In the formula, is the three-dimensional tensor form of X, λ is the weight parameter used to weigh the regularization term and the data term, is the calculation formula of the comprehensive measurement value of the three-dimensional image X, and y is the actual detection value of the detector; The expression of the optimization problem after introducing the intermediate variable Z is:

6. The three-dimensional synthetic imaging method based on coded energy spectrum according to claim 5, characterized in that: The Gauss-Newton method is used to update the three-dimensional image X, and the iterative shearing algorithm is used to update the intermediate variable Z.

7. A three-dimensional synthetic camera device based on coded energy spectrum, characterized in that: The method comprises a memory and a processor, wherein the memory stores a computer program, and the processor calls the computer program to execute the steps of any one of the methods according to claims 1 to 6.

Images