Image Reconstruction Method, System, Electronic Device and Storage Medium of Energy Spectrum CT

Abstract

The present invention discloses an image reconstruction method, system, electronic device and storage medium for spectral CT. The image reconstruction method includes: acquiring spectral projection data and inputting the spectral projection data into an image reconstruction model to obtain a spectral image to be updated; extracting similar image patches from the spectral image to be updated to construct a hypercomplex matrix; updating the image reconstruction model based on the hypercomplex matrix, and obtaining an updated spectral image based on the updated image reconstruction model. The image reconstruction method of the present invention uses a hypercomplex matrix to represent multi-dimensional spectral images, encodes the inherent low-rank property of spectral images in a holistic manner in the hypercomplex space, and better preserves the integrity of the multi-dimensional data structure; fully characterizes the cross-space non-local self-similarity and the global correlation along the energy spectrum contained in all energy segment images, suppresses quantum noise and weakens the strip artifacts caused by photon starvation.

Description

Technical Field

[0001] The present invention relates to the field of medical imaging technology, and in particular to an image reconstruction method, system, electronic equipment and storage medium for spectral CT. Background Art

[0002] Spectral CT (Computed Tomography) has broad clinical application prospects due to its energy resolution and material differentiation capabilities. Reconstructing high-quality spectral images is one of the keys to achieving accurate material decomposition. However, due to the small number of photons received in each energy band (only 1 / 2 to 1 / 8 of conventional CT), the signal-to-noise ratio of a single energy band image is low, especially in high attenuation areas such as bones and shoulders, where severe photon starvation occurs, resulting in obvious strip artifacts in the reconstructed image (especially low-energy band images). Therefore, spectral CT reconstruction places high demands on the denoising and artifact removal capabilities of the image reconstruction algorithm.

[0003] Studies have shown that the use of prior information inherent in spectral images can effectively improve the quality of reconstruction. Piece-Wise Constant is one of the most commonly used prior knowledge in the field of CT reconstruction. Based on this, the existing technology has developed regularization terms such as Total Variation (TV) and structural tensor TV for iterative reconstruction of spectral CT. However, these methods only consider the spatial prior within a single energy band image, while ignoring the correlation between different energy bands. Obviously, images of different energy bands have highly similar morphological structures and texture features, and only differ in CT values. In order to make full use of the strong correlation between energy bands, low-rank-based prior information is used in spectral CT reconstruction, and methods such as prior rank, intensity and sparsity model (PRISM), tensor PRISM, tensor-based dictionary learning algorithm, self-similarity-assisted reconstruction model in the atlas tensor, and fourth-order non-local tensor decomposition model have been proposed. However, some of these methods ignore the similarities of image patches at different spatial locations, while others cause the loss of multi-dimensional data structure information by expanding tensors into matrices during low-rank regularization. Summary of the invention

[0004] The technical problem to be solved by the present invention is to overcome the defect that the tensor-based spectral CT image reconstruction method in the prior art causes the loss of multidimensional data structure information, and to provide an image reconstruction method, system, electronic device and storage medium for spectral CT.

[0005] The present invention solves the above technical problems through the following technical solutions:

[0006] The present invention provides an image reconstruction method for spectral CT, the image reconstruction method comprising:

[0007] Acquiring energy spectrum projection data, and inputting the energy spectrum projection data into an image reconstruction model to obtain an energy spectrum image to be updated;

[0008] Extracting similar image blocks from the energy spectrum image to be updated to construct a hypercomplex matrix;

[0009] The image reconstruction model is updated based on the hypercomplex matrix, and an updated energy spectrum image is obtained based on the updated image reconstruction model.

[0010] Preferably, the step of updating the image reconstruction model based on the hypercomplex matrix and obtaining an updated energy spectrum image based on the updated image reconstruction model comprises:

[0011] Iteratively updating the image reconstruction model, that is, updating the image reconstruction model in each iteration based on the hypercomplex matrix of the energy spectrum image reconstructed last time;

[0012] After the step of iteratively updating the image reconstruction model, the image reconstruction method further comprises:

[0013] When the iterative updating process satisfies the preset convergence condition, the iterative reconstruction is stopped and the current energy spectrum image is determined as the final energy spectrum image.

[0014] Preferably, the step of updating the image reconstruction model based on the hypercomplex matrix comprises:

[0015] updating a regularization term of the image reconstruction model based on the hypercomplex matrix;

[0016] The regularization term is the sum of several sub-regularization terms, and the sub-regularization terms include the sum of the ranks of the hypercomplex matrix.

[0017] Preferably, the similar image blocks are image blocks with similar structures and containing different energy segments.

[0018] Preferably, the hypercomplex matrix comprises a quaternion matrix or an octonion matrix;

[0019] The step of extracting similar image blocks from the energy spectrum image to be updated to construct a hypercomplex matrix includes:

[0020] Cutting the energy spectrum image to be updated into a plurality of first three-dimensional image blocks;

[0021] Searching for a set of first three-dimensional image blocks of the most similar preset number for each of the first three-dimensional image blocks;

[0022] Stretching and stacking each group of the most similar first three-dimensional image blocks to obtain a second three-dimensional image block;

[0023] The energy segments of the second three-dimensional image block are substituted into the real part and / or the imaginary part of the hyper-complex matrix respectively to obtain a plurality of hyper-complex matrices.

[0024] Preferably, the iterative update process satisfies a preset convergence condition including that the number of iterative updates reaches a preset threshold number of iterations; and / or,

[0025] The steps of obtaining an updated energy spectrum image based on the updated image reconstruction model include:

[0026] An alternating minimization algorithm is used to solve the minimization problem of the image reconstruction model to obtain an updated energy spectrum image.

[0027] Preferably, after the step of updating the regularization term of the image reconstruction model based on the hypercomplex matrix, the image reconstruction method further comprises:

[0028] The rank function of the regularization term is approximated by using a substitute function of the rank function; the rank function is used to characterize the sum of the ranks of the hypercomplex matrix.

[0029] The present invention also provides an image reconstruction system for spectral CT, the image reconstruction system comprising:

[0030] An energy spectrum projection data acquisition module, used for acquiring energy spectrum projection data, and inputting the energy spectrum projection data into an image reconstruction model to obtain an energy spectrum image to be updated;

[0031] A hypercomplex matrix construction module, used for extracting similar image blocks from the energy spectrum image to be updated to construct a hypercomplex matrix;

[0032] A reconstruction model updating module is used to update the image reconstruction model based on the hypercomplex matrix, and obtain an updated energy spectrum image based on the updated image reconstruction model.

[0033] Preferably, the reconstruction model updating module is further used for iteratively updating the image reconstruction model, that is, in each iteration, updating the image reconstruction model based on the hypercomplex matrix of the energy spectrum image obtained by the previous reconstruction;

[0034] The image reconstruction system further comprises:

[0035] The energy spectrum image determination module is used to stop the iterative reconstruction and determine the current energy spectrum image as the final energy spectrum image when the iterative update process meets the preset convergence condition.

[0036] Preferably, the reconstruction model updating module is specifically used to update the regularization term of the image reconstruction model based on the hypercomplex matrix;

[0037] The regularization term is the sum of several sub-regularization terms, and the sub-regularization terms include the sum of the ranks of the hypercomplex matrix.

[0038] Preferably, the similar image blocks are image blocks with similar structures and containing different energy segments.

[0039] Preferably, the hypercomplex matrix comprises a quaternion matrix or an octonion matrix;

[0040] The hypercomplex matrix construction module is specifically used to cut the energy spectrum image to be updated into a plurality of first three-dimensional image blocks;

[0041] The hypercomplex matrix construction module is specifically used to search for a group of first three-dimensional image blocks of a preset number that are most similar to each of the first three-dimensional image blocks;

[0042] The hypercomplex matrix construction module is specifically used to stretch and stack each group of the most similar first three-dimensional image blocks to obtain a second three-dimensional image block;

[0043] The hypercomplex matrix construction module is specifically used to substitute the energy segments of the second three-dimensional image block into the real part and / or the imaginary part of the hypercomplex matrix to obtain a plurality of hypercomplex matrices.

[0044] Preferably, the preset convergence condition in the energy spectrum image determination module includes that the number of iteration updates reaches a preset iteration number threshold; and / or,

[0045] The reconstruction model updating module is specifically used to solve the minimization problem of the image reconstruction model by adopting an alternating minimization algorithm to obtain an updated energy spectrum image.

[0046] Preferably, the reconstruction model updating module is also used to approximate the rank function of the regularization term using a substitute function of the rank function; the rank function is used to characterize the sum of the ranks of the hypercomplex matrix.

[0047] The present invention also 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 implements the above-mentioned image reconstruction method for spectral CT when executing the computer program.

[0048] The present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the above-mentioned image reconstruction method of spectral CT.

[0049] The positive and progressive effects of the present invention are:

[0050] The image reconstruction method of energy spectrum CT provided by the present invention provides extracting similar image blocks from the energy spectrum image to be updated to construct a hypercomplex matrix, updating the image reconstruction model based on the hypercomplex matrix, and obtaining an updated energy spectrum image based on the updated image reconstruction model. In this iterative process, a hypercomplex matrix is ​​used to represent the multi-dimensional energy spectrum image, and the inherent low rank of the energy spectrum image is encoded in a holistic manner in the hypercomplex space, which makes up for the defect that the reconstruction model based on tensor representation will destroy its inherent structure when the tensor is expanded into a matrix in all directions, and better preserves the integrity of the multi-dimensional data structure; fully describes the non-local self-similarity across space and the global correlation along the energy spectrum contained in all energy band images, suppresses quantum noise and weakens the stripe artifacts caused by photon starvation, and achieves higher quality reconstruction. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] Figure 1 This is a first flow chart of the image reconstruction method of spectral CT in Example 1 of the present invention.

[0052] Figure 2 This is a second flow chart of the image reconstruction method of spectral CT in Example 1 of the present invention.

[0053] Figure 3 This is a comparison diagram of the effects of reconstructing the XCAT digital phantom using the image reconstruction method of the spectral CT of the prior art and the image reconstruction method of the present embodiment in Example 1 of the present invention.

[0054] Figure 4 This is a comparison diagram of the effects of reconstructing the Gammex real phantom using the image reconstruction method of the spectral CT of the prior art and the image reconstruction method of the present embodiment in Example 1 of the present invention.

[0055] Figure 5 This is a first structural schematic diagram of the image reconstruction system of the energy spectrum CT in Example 2 of the present invention.

[0056] Figure 6 This is a second structural schematic diagram of the image reconstruction system of the energy spectrum CT in Example 2 of the present invention.

[0057] Figure 7 Schematic diagram of the structure of an electronic device in Embodiment 3 of the present invention. DETAILED DESCRIPTION

[0058] The present invention is further described below by way of examples, but the present invention is not limited to the scope of the examples.

[0059] Example 1

[0060] Please refer to Figure 1, which is a first flow chart of the image reconstruction method of spectral CT in this embodiment. Specifically, Figure 1 As shown, the image reconstruction method includes:

[0061] S101, obtaining energy spectrum projection data, and inputting the energy spectrum projection data into the image reconstruction model to obtain the energy spectrum image to be updated; specifically, the image reconstruction model in this embodiment includes a data fidelity term and a regularization term. The role of the fidelity term is to ensure that an energy spectrum image consistent with the energy spectrum projection data is obtained. The regularization term is equivalent to adding an additional constraint in the process of minimizing the fidelity term, that is, not only to satisfy the fidelity term is small enough, but also to be completed under the constraint that the regularization term cannot be too large. The energy spectrum (Spectral) of this embodiment includes dual-energy (Dual-Energy) and multi-energy (Multi-Energy). That is, the energy spectrum of this embodiment can be "dual energy" based on technologies such as dual sources, double-layer detectors, split filters (Split Filter), fast kV switching, and temporally sequential scanning (Temporally Sequential Scanning), or "multi-energy" based on technologies such as continuous multiple different kV scans, photon counting detectors, etc. Furthermore, the multiple energy segments of the energy spectrum CT of this embodiment can be two energy spectra of dual-energy CT, or multiple energy spectra or multiple energy bins (Energy Bin) of multi-energy CT.

[0062] S102: extract similar image blocks from the energy spectrum image to be updated to construct a hypercomplex matrix.

[0063] S103 . Update the image reconstruction model based on the hypercomplex matrix, and obtain an updated energy spectrum image based on the updated image reconstruction model.

[0064] In this embodiment, similar image blocks are image blocks with similar structures that contain different energy bands. Specifically, in the same CT image, there is a strong similarity between image blocks located at different spatial positions, which is called non-local self-similarity across space. On the other hand, images of different energy bands of energy spectrum CT have highly similar morphological structures and texture features, and this similarity is called global correlation along the energy spectrum. These two prior features make the energy spectrum image have good low rank. Based on this low rank, a high-quality energy spectrum image can be reconstructed from projection data with strong noise.

[0065] In an optional embodiment, the hypercomplex matrix includes a quaternion matrix; specifically, the quaternion matrix is ​​used to represent the multi-dimensional (spatial dimension and energy dimension) energy spectrum image as a whole, rather than expanding the energy spectrum image into a real matrix or real vector along a certain dimensional direction, thereby avoiding the loss of useful structural information along other dimensional directions.

[0066] Quaternions are a generalization of complex numbers, consisting of a real part and three imaginary parts, for example:

[0067]

[0068] in, is a quaternion, q 0 ,q 1 ,q 2 ,q 3 are all real numbers, i, j, k are three imaginary units, Represents a quaternion space.

[0069] Similar to a real matrix, a quaternion matrix is ​​a matrix composed of quaternions, denoted as:

[0070]

[0071] in, is a quaternion matrix, All are real matrices.

[0072] When using X 0 ,X 1 ,X 2 ,X 3 When representing different energy bands of the energy spectrum image, the quaternion matrix It can be used to represent an energy spectrum image with four energy bands. 0 If set to zero, It can be used to represent the energy spectrum image with three energy bands. 2 and X 3 Set to zero (at this time degenerates into a complex matrix), then It can be used to represent energy spectrum images with two energy bands. This representation method is consistent with the multi-channel characteristics of energy spectrum images and ensures the integrity of the multi-dimensional data structure.

[0073] Furthermore, quaternion representation can be extended to octonion representation.

[0074] Similar to quaternions, octonions are a further generalization of complex numbers, consisting of one real part and seven imaginary parts, for example:

[0075]

[0076] in, is an octonion, q 0 ,q 1 ,q 2 ,q 3 ,q 4 ,q 5 ,q6 ,q 7 are all real numbers, i, j, k, l, m, n, o are seven imaginary units, Represents the octonion space. The octonion matrix is ​​a matrix composed of octonions, recorded as:

[0077]

[0078] in, is an octonion matrix, All are real matrices.

[0079] When using X 0 ,X 1 ,X 2 ,X 3 ,X 4 ,X 5 ,X 6 ,X 7 When representing different energy bands of the energy spectrum image, the octonion matrix It can be used to represent an energy spectrum image with eight energy bands. 0 If set to zero, It can be used to represent an energy spectrum image with seven energy bands.

[0080] Therefore, it is only necessary to replace the quaternion representation with the octonion representation, and other technical contents such as the objective function, optimization strategy, algorithm framework, and iteration scheme described below in this embodiment do not need to be changed, so that the image reconstruction method of energy spectrum CT based on octonion representation can be directly obtained to achieve high-quality reconstruction of seven-segment images or eight-segment images.

[0081] Similarly, those skilled in the art should know that based on the technical content disclosed in this embodiment, quaternion representation and octonion representation can be further extended to hypercomplex representation, and those skilled in the art can obtain an image reconstruction method of spectral CT based on hypercomplex representation. Therefore, the image reconstruction method of spectral CT based on hypercomplex representation is within the protection scope of the present invention.

[0082] For the convenience of explanation, this embodiment takes the reconstruction of a single energy spectrum image with three energy bands as an example. The energy spectrum image to be reconstructed is represented as h and w represent the height and width of the image, respectively. The general form of the spectral CT iterative reconstruction model can be written as:

[0083]

[0084] in, S = h × w × 3 is a vectorized energy spectrum image, that is, a column vector formed by concatenating the vectorized images of each energy band. f(x) is a data fidelity term, g(x) is a regularization term, and κ is a regularization parameter. In this embodiment, f(x) is defined as follows:

[0085]

[0086] in, is the projection data, is the system matrix, is the weighting matrix. The value of W is not restricted and can usually be taken as the inverse matrix of the covariance matrix of the observation value y.

[0087] Please refer to Figure 2 , which is a second flow chart of the image reconstruction method of spectral CT in this embodiment. Specifically, Figure 2 As shown, in this embodiment, step S102 includes:

[0088] S1021, cutting the energy spectrum image to be updated into a plurality of first three-dimensional image blocks; specifically, taking a quaternion matrix as an example, cutting the energy spectrum image with a size of h×w×3 into a plurality of first three-dimensional image blocks; Cut into L size 3D image blocks l=1,…,L.

[0089] S1022, for each first three-dimensional image block, respectively search for a group of first three-dimensional image blocks of a preset number that are most similar; specifically, taking a quaternion matrix as an example, search for the first three-dimensional image blocks that are most similar to the first three-dimensional image block from all L 3D image blocks. The most similar n-1 3D image blocks, a total of n 3D image blocks (i.e. and its n-1 similar image patches).

[0090] S1023, stretching and stacking each group of the most similar first three-dimensional image blocks to obtain a second three-dimensional image block; specifically, taking the quaternion matrix as an example, for the n 3D image blocks obtained in step S1022, stretching each energy segment of each 3D image block into a column vector to obtain n 2D image blocks of size m×3; stacking the n 2D image blocks obtained above into a 3D image block of size m×n×3

[0091] S1024, respectively substitute the energy segments of the second three-dimensional image block into the real part and / or imaginary part of the hypercomplex matrix to obtain a plurality of hypercomplex matrices; specifically, taking the quaternion matrix as an example, substitute the energy segments of the second three-dimensional image block into the real part and / or imaginary part of the hypercomplex matrix to obtain a plurality of hypercomplex matrices; The three energy segments are substituted into the three imaginary parts of the quaternion matrix to obtain a quaternion matrix of size m×n After repeating steps S1022 and S1023, L quaternion matrices of size m×n can be obtained.

[0092] In this embodiment, step S103 includes:

[0093] The image reconstruction model is updated iteratively, that is, in each iteration, the image reconstruction model is updated based on the hypercomplex matrix of the energy spectrum image obtained by the previous reconstruction.

[0094] Specifically, in an optional implementation, step S103 includes:

[0095] S1031, updating the regularization term of the image reconstruction model based on the hypercomplex matrix;

[0096] The regularization term is the sum of several sub-regularization terms, and the sub-regularization term includes the sum of the ranks of the hypercomplex matrix. Specifically, in this embodiment, g(x) is defined as follows:

[0097]

[0098] Among them, g r (x), r = 1, ..., R is R sub-regular terms, one of which is g 1 (x) is L quaternion matrix The sum of the ranks of the other sub-regular terms g 2 (x),…,g R (x) is not limited and can be any possible regularization term. It is worth noting that directly defining “one of the sub-regularization terms” as g 1 (x), but it can actually be g 2 (x),…,g R In this embodiment, for the sake of illustration, it can be directly assumed that g 2 (x) = ... = g R (x) = 0, then g(x) = g 1 (x).

[0099] Taking the above quaternion matrix as an example, since the quaternion matrix Different columns are similar image blocks and different imaginary parts are similar energy segments, so It should have obvious low rank. To this end, we can l=1,…,L performs low-rank regularization to suppress noise, which is exactly the sub-regularization term g 1 The function of (x).

[0100] In an optional implementation, after step S1031, step S103 further includes:

[0101] S1032. Approximate the rank function of the regularization term using a substitute function of the rank function; the rank function is used to characterize the sum of the ranks of the hypercomplex matrix.

[0102] The rank of is calculated as follows:

[0103]

[0104]

[0105] in, is the quaternion matrix The singular value decomposition of diag(∑ l ) is the singular value matrix ∑ l The column vector composed of the main diagonal elements of is the singular value vector, ‖·‖ 0 yes Norm.

[0106] Since optimization problems involving rank functions are difficult to solve directly, a substitute function of the rank function is often used to approximate the rank function before solving it. In this embodiment, the weighted Schatten-τ norm can be used as a substitute function for the rank function, as follows:

[0107]

[0108] Among them, 0<τ<1, i=1,…,I is the quaternion matrix The singular values ​​of i ,i=1,…,I is a monotonically non-decreasing weight, that is, it satisfies ω 1 ≤ω 2 ≤…≤ω I The optimization problem after using the rank function's substitute function to approximate the rank function is as follows:

[0109]

[0110] Assume g 2 (x) = ... = g R (x) = 0, then we get:

[0111]

[0112] Furthermore, in addition to the weighted Schatten-τ norm, nuclear norm, weighted nuclear norm, capped nuclear norm, truncated nuclear norm, Schatten-τ norm, truncated Schatten-τ norm, smoothly clipped absolute deviation (SCAD) penalty, minimax concave penalty (MCP), logarithm penalty, log-determinant penalty, Geman penalty, Laplace penalty, etc. can also be used as alternative functions to the rank function.

[0113] The present invention does not limit the weight ω of the weighted Schatten-τ norm i , i = 1, ..., I, the weight of the weighted Schatten-τ norm only needs to be monotonically non-decreasing, that is, it only needs to satisfy ω 1 ≤ω 2 ≤…≤ω I This means that the quaternion matrix The singular values ​​of The larger the value, the more weight ω is applied to it. i The smaller it is, the less it is compressed. This is expected, because larger singular values ​​encode texture information and structural details of the image, so they should be preserved more.

[0114] In an optional implementation, after step S1032, step S103 further includes:

[0115] S1033, using an alternating minimization algorithm to solve the minimization problem of the image reconstruction model to obtain an updated energy spectrum image;

[0116] After step S1033, the image reconstruction method of spectral CT further includes:

[0117] S104, judging whether the iterative updating process satisfies the preset convergence condition, if so, executing step S105, if not, taking the current energy spectrum image obtained in this iteration as the energy spectrum image to be updated, and executing step S102.

[0118] Specifically, the convergence condition preset in step S104 may include one of the following conditions:

[0119] The objective function of the image reconstruction model reaches a minimum value;

[0120] The number of iteration updates reaches the preset iteration threshold.

[0121] S105 , stopping the iterative reconstruction and determining the current energy spectrum image as the final energy spectrum image.

[0122] Taking the above quaternion matrix representation as an example, the augmented Lagrangian function of the optimization problem to be solved is:

[0123]

[0124] in, is an auxiliary variable, is the Lagrange multiplier, μ is the penalty parameter. Then the framework of Alternating Direction Method of Multipliers (ADMM) can be used to solve it, and the iterative scheme is as follows:

[0125]

[0126] Among them, (a) is the fidelity subproblem, (b) is the denoising subproblem, (c) is the Lagrange multiplier update formula, (d) is the penalty parameter update formula, and the superscript k represents the kth iteration. After several iterations, the reconstructed images of all energy bands can be obtained.

[0127] Further, in solving the optimization problem min x When f(x)+κg(x), in addition to ADMM, you can also use alternating minimization algorithms such as Fast ADMM (FADMM), Iterative Shrinkage-Thresholding Algorithm (ISTA), Fast ISTA (FISTA), Forward-Backward Splitting (FBS), Generalized FBS (GFBS), Douglas-Rachford Splitting (DRS), Primal-Dual Splitting (PDS), Half Quadratic Splitting (HQS), and Split Bregman.

[0128] The following takes the reconstruction of XCAT (a 4D computer phantom) digital phantom and Gammex (a CT phantom) real phantom as examples to illustrate the effectiveness of the proposed algorithm. According to the imaging principle of photon counting CT, the XACT phantom is projected to generate projection data of three energy bands [20,57]keV, [57,82]keV and [82,120]keV. The traditional filtered backprojection (Filtered Backprojection, FBP) algorithm, the fourth-order nonlocal tensor decomposition (Fourth-Order Nonlocal Tensor Decomposition, FONT) algorithm based on tensor representation and the image reconstruction method of spectral CT proposed in this embodiment are used for reconstruction, and the following is obtained: Figure 3 The results are shown in Figure 2. The first row is the reconstructed image of FBP, the second row is the reconstructed image of FONT, and the third row is the reconstructed image of the proposed image reconstruction method. The first to third columns correspond to the three energy bands of [20,57]keV, [57,82]keV, and [82,120]keV, respectively, and the display windows are [0.18,0.32]cm -1 、[0.15,0.23]cm -1 and [0.13,0.21]cm -1 .

[0129] The Gammex phantom was scanned with a photon counting CT machine at a tube current of 300 mA, including three energy bands: [20,55]keV, [55,80]keV, and [80,140]keV. The image reconstruction was performed using the traditional FBP algorithm, the FONT algorithm based on tensor representation, and the energy spectrum CT image reconstruction method proposed in this embodiment, respectively, to obtain the following: Figure 4 The results are shown in Figure 1. The first row is the reconstructed image of FBP, the second row is the reconstructed image of FONT, and the third row is the reconstructed image of the proposed image reconstruction method. The first to third columns correspond to the three energy ranges of [20,55]keV, [55,80]keV, and [80,140]keV, and the display windows are [840,1340]HU, [820,1320]HU, and [800,1300]HU, respectively.

[0130] from Figure 3 and Figure 4It can be seen that compared with the traditional FBP algorithm, the proposed image reconstruction method significantly suppresses quantum noise, effectively weakens the stripe artifacts caused by photon starvation, and fully improves the reconstruction quality of energy spectrum images (especially low-energy images). On the other hand, compared with the FONT algorithm based on tensor representation, the proposed image reconstruction method uses quaternion representation, that is, the quaternion matrix is ​​used to completely represent the entire 3D image block that should have low rank. The singular value decomposition and subsequent singular value shrinkage operation (Shrinkage Operation) of the quaternion matrix can be directly performed without expanding the quaternion matrix into a 2D real matrix. This makes up for the defect that the FONT algorithm destroys its internal structure when expanding the tensor into a matrix in all directions, better preserves the integrity of the multidimensional data structure, avoids the loss of useful structural information, and makes the reconstructed image higher in quality.

[0131] Example 2

[0132] Please refer to Figure 5 , which is a first structural schematic diagram of the image reconstruction system of the energy spectrum CT in this embodiment. Specifically, Figure 5 As shown, the image reconstruction system comprises:

[0133] The energy spectrum projection data acquisition module 1 is used to acquire energy spectrum projection data and input the energy spectrum projection data into the image reconstruction model to obtain the energy spectrum image to be updated; specifically, the image reconstruction model in this embodiment includes a data fidelity term and a regularization term. The role of the fidelity term is to ensure that an energy spectrum image consistent with the energy spectrum projection data is obtained. The regularization term is equivalent to adding an additional constraint in the process of minimizing the fidelity term, that is, not only the fidelity term must be small enough, but also the regularization term must not be too large. The energy spectrum of this embodiment includes dual energy and multi-energy. That is, the energy spectrum of this embodiment can be "dual energy" based on technologies such as dual sources, dual-layer detectors, split filters, fast kV switching, and time-series scanning, or "multi-energy" based on technologies such as continuous multiple different kV scans and photon counting detectors. Furthermore, the multiple energy segments of the energy spectrum CT in this embodiment can be two energy spectra of dual-energy CT, or multiple energy spectra or multiple energy bins of multi-energy CT.

[0134] A hypercomplex matrix construction module 2, used for extracting similar image blocks from the energy spectrum image to be updated to construct a hypercomplex matrix;

[0135] The reconstruction model updating module 3 is used to update the image reconstruction model based on the hypercomplex matrix, and obtain an updated energy spectrum image based on the updated image reconstruction model.

[0136] In this embodiment, similar image blocks are image blocks with similar structures that contain different energy bands. Specifically, in the same CT image, there is a strong similarity between image blocks located at different spatial positions, which is called non-local self-similarity across space. On the other hand, images of different energy bands of energy spectrum CT have highly similar morphological structures and texture features, and this similarity is called global correlation along the energy spectrum. These two prior features make the energy spectrum image have good low rank. Based on this low rank, a high-quality energy spectrum image can be reconstructed from projection data with strong noise.

[0137] In an optional embodiment, the hypercomplex matrix includes a quaternion matrix; specifically, the quaternion matrix is ​​used to represent the multi-dimensional (spatial dimension and energy dimension) energy spectrum image as a whole, rather than expanding the energy spectrum image into a real matrix or real vector along a certain dimensional direction, thereby avoiding the loss of useful structural information along other dimensional directions.

[0138] Quaternions are a generalization of complex numbers, consisting of a real part and three imaginary parts, for example:

[0139]

[0140] in, is a quaternion, q 0 ,q 1 ,q 2 ,q 3 are all real numbers, i, j, k are three imaginary units, Represents a quaternion space.

[0141] Similar to a real matrix, a quaternion matrix is ​​a matrix composed of quaternions, denoted as:

[0142]

[0143] in, is a quaternion matrix, All are real matrices.

[0144] When using X 0 ,X 1 ,X 2 ,X 3 When representing different energy bands of the energy spectrum image, the quaternion matrix It can be used to represent an energy spectrum image with four energy bands. 0 If set to zero, It can be used to represent the energy spectrum image with three energy bands. 2 and X 3 Set to zero (at this time degenerates into a complex matrix), then It can be used to represent energy spectrum images with two energy bands. This representation method is consistent with the multi-channel characteristics of energy spectrum images and ensures the integrity of the multi-dimensional data structure.

[0145] Furthermore, quaternion representation can be extended to octonion representation.

[0146] Similar to quaternions, octonions are a further generalization of complex numbers, consisting of one real part and seven imaginary parts, for example:

[0147]

[0148] in, is an octonion, q 0 ,q 1 ,q 2 ,q 3 ,q 4 ,q 5 ,q 6 ,q 7 are all real numbers, i, j, k, l, m, n, o are seven imaginary units, Represents the octonion space. The octonion matrix is ​​a matrix composed of octonions, recorded as:

[0149]

[0150] in, is an octonion matrix, All are real matrices.

[0151] When using X 0 ,X 1 ,X 2 ,X 3 ,X 4 ,X 5 ,X 6 ,X 7 When representing different energy bands of the energy spectrum image, the octonion matrix It can be used to represent an energy spectrum image with eight energy bands. 0 If set to zero, It can be used to represent an energy spectrum image with seven energy bands.

[0152] Therefore, it is only necessary to replace the quaternion representation with the octonion representation, and other technical contents such as the objective function, optimization strategy, algorithm framework, and iteration scheme described below in this embodiment do not need to be changed, so that the image reconstruction method of energy spectrum CT based on octonion representation can be directly obtained to achieve high-quality reconstruction of seven-segment images or eight-segment images.

[0153] Similarly, those skilled in the art should know that based on the technical content disclosed in this embodiment, quaternion representation and octonion representation can be further extended to hypercomplex representation, and those skilled in the art can obtain an image reconstruction method of spectral CT based on hypercomplex representation. Therefore, the image reconstruction method of spectral CT based on hypercomplex representation is within the protection scope of the present invention.

[0154] For the convenience of explanation, this embodiment takes the reconstruction of a single energy spectrum image with three energy bands as an example. The energy spectrum image to be reconstructed is represented as h and w represent the height and width of the image, respectively. The general form of the spectral CT iterative reconstruction model can be written as:

[0155]

[0156] in, S = h × w × 3 is a vectorized energy spectrum image, that is, a column vector formed by concatenating the vectorized images of each energy band. f(x) is a data fidelity term, g(x) is a regularization term, and κ is a regularization parameter. In this embodiment, f(x) is defined as follows:

[0157]

[0158] in, is the projection data, is the system matrix, is the weighting matrix. The value of W is not restricted and can usually be taken as the inverse matrix of the covariance matrix of the observation value y.

[0159] In an optional implementation, the hypercomplex matrix construction module 2 is specifically used to cut the energy spectrum image to be updated into a plurality of first three-dimensional image blocks; specifically, taking the quaternion matrix as an example, the energy spectrum image x with a size of h×w×3 is cut into L first three-dimensional image blocks with a size of 3D image blocks l=1,…,L.

[0160] The hypercomplex matrix construction module 2 is specifically used to search for a set of first three-dimensional image blocks of the most similar preset number for each first three-dimensional image block; specifically, taking the quaternion matrix as an example, searching for the first three-dimensional image blocks with the same preset number from all L 3D image blocks. The most similar n-1 3D image blocks, a total of n 3D image blocks (i.e. and its n-1 similar image patches).

[0161] The hypercomplex matrix construction module 2 is specifically used to stretch and stack each group of the most similar first three-dimensional image blocks to obtain a second three-dimensional image block; specifically, taking the quaternion matrix as an example, for the above-mentioned n 3D image blocks, each energy segment of each 3D image block is stretched into a column vector to obtain n 2D image blocks of size m×3; the above-mentioned n 2D image blocks are stacked into a 3D image block of size m×n×3

[0162] The hypercomplex matrix construction module 2 is specifically used to substitute the energy segments of the second three-dimensional image block into the real part and / or imaginary part of the hypercomplex matrix to obtain a plurality of hypercomplex matrices; specifically, taking the quaternion matrix as an example, the 3D image block The three energy segments are substituted into the three imaginary parts of the quaternion matrix to obtain a quaternion matrix of size m×n Repeatedly calling the hypercomplex matrix construction module 2 can obtain L quaternion matrices of size m×n

[0163] Please refer to Figure 6 , which is a second structural schematic diagram of the image reconstruction system of the energy spectrum CT in this embodiment. Specifically, Figure 6 As shown, the reconstruction model updating module 3 is used to iteratively update the image reconstruction model, that is, in each iteration, the image reconstruction model is updated based on the hypercomplex matrix of the energy spectrum image obtained by the previous reconstruction;

[0164] In an optional implementation, the reconstruction model updating module 3 is specifically used to update the regularization term of the image reconstruction model based on the hypercomplex matrix; the regularization term is the sum of several sub-regularization terms, and the sub-regularization term includes the sum of the ranks of the hypercomplex matrix. In this implementation, g(x) is defined as follows:

[0165]

[0166] Among them, g r (x), r = 1, ..., R is R sub-regular terms, one of which is g 1 (x) is L quaternion matrix l=1,…,the sum of the ranks of L, and other sub-regular terms g 2 (x),…,g R (x) is not limited and can be any possible regularization term. It is worth noting that directly defining “one of the sub-regularization terms” as g 1 (x), but it can actually be g 2 (x),…,g R In this embodiment, for the sake of illustration, it can be directly assumed that g 2 (x) = ... = g R(x) = 0, then g(x) = g 1 (x).

[0167] Taking the above quaternion matrix as an example, since the quaternion matrix Different columns are similar image blocks and different imaginary parts are similar energy segments, so It should have obvious low rank. To this end, we can Low-rank regularization is performed to suppress noise, which is exactly the sub-regularization term g 1 The function of (x).

[0168] In an optional implementation, the reconstruction model updating module 3 is further used to approximate the rank function of the regularization term using a substitute function of the rank function; the rank function is used to characterize the sum of the ranks of the hypercomplex matrix.

[0169] The rank of is calculated as follows:

[0170]

[0171]

[0172] in, is the quaternion matrix The singular value decomposition of diag(∑ l ) is the singular value matrix ∑ l The column vector composed of the main diagonal elements of is the singular value vector, ‖·‖ 0 yes Norm.

[0173] Since optimization problems involving rank functions are difficult to solve directly, a substitute function of the rank function is often used to approximate the rank function before solving it. In this embodiment, the weighted Schatten-τ norm can be used as a substitute function for the rank function, as follows:

[0174]

[0175] Among them, 0<τ<1, i=1,…,I is the quaternion matrix The singular values ​​of i ,i=1,…,I is a monotonically non-decreasing weight, that is, it satisfies ω 1 ≤ω 2 ≤…≤ω I The optimization problem after using the rank function's substitute function to approximate the rank function is as follows:

[0176]

[0177] Assume g 2 (x) = ... = gR (x) = 0, then we get:

[0178]

[0179] Furthermore, in addition to the weighted Schatten-τ norm, the nuclear norm, weighted nuclear norm, capped nuclear norm, truncated nuclear norm, Schatten-τ norm, truncated Schatten-τ norm, smooth clipping absolute deviation penalty, minimax concave penalty, Logarithm penalty, Log-Determinant penalty, Geman penalty, Laplace penalty, etc. can also be used as alternative functions to the rank function.

[0180] The present invention does not limit the weight ω of the weighted Schatten-τ norm i , i = 1, ..., I, the weight of the weighted Schatten-τ norm only needs to be monotonically non-decreasing, that is, it only needs to satisfy ω 1 ≤ω 2 ≤…≤ω I This means that the quaternion matrix The singular values ​​of The larger the value, the more weight ω is applied to it. i The smaller it is, the less it is compressed. This is expected, because larger singular values ​​encode texture information and structural details of the image, so they should be preserved more.

[0181] In an optional implementation, the reconstruction model updating module 3 is further configured to adopt an alternating minimization algorithm to solve the minimization problem of the image reconstruction model to obtain an updated energy spectrum image;

[0182] The image reconstruction system of the energy spectrum CT further includes:

[0183] The energy spectrum image determination module 4 is used to stop the iterative reconstruction and determine the current energy spectrum image as the final energy spectrum image when the iterative updating process meets the preset convergence condition.

[0184] Specifically, the convergence condition preset in the energy spectrum image determination module 4 may include one of the following conditions:

[0185] The objective function of the image reconstruction model reaches a minimum value;

[0186] The number of iteration updates reaches the preset iteration threshold.

[0187] Taking the above quaternion matrix representation as an example, the augmented Lagrangian function of the optimization problem to be solved is:

[0188]

[0189] in, is an auxiliary variable, is the Lagrange multiplier, μ is the penalty parameter. Then the ADMM framework can be used to solve it, and the iterative scheme is as follows:

[0190]

[0191] Among them, (a) is the fidelity subproblem, (b) is the denoising subproblem, (c) is the Lagrange multiplier update formula, (d) is the penalty parameter update formula, and the superscript k represents the kth iteration. After several iterations, the reconstructed images of all energy bands can be obtained.

[0192] Further, in solving the optimization problem min x When f(x)+κg(x), in addition to ADMM, you can also use alternating minimization algorithms such as FADMM, ISTA, FISTA, FBS, GFBS, DRS, PDS, HQS, and split Bregman.

[0193] The reconstruction of the XCAT digital phantom and the Gammex real phantom are used as examples to illustrate the effectiveness of the proposed algorithm. According to the imaging principle of photon counting CT, the XACT phantom is projected to generate projection data of three energy bands: [20,57]keV, [57,82]keV, and [82,120]keV. The traditional FBP algorithm, the FONT algorithm based on tensor representation, and the image reconstruction method of spectral CT proposed in this embodiment are used for reconstruction, and the following are obtained: Figure 3 The results are shown in Figure 2. The first row is the reconstructed image of FBP, the second row is the reconstructed image of FONT, and the third row is the reconstructed image of the proposed image reconstruction method. The first to third columns correspond to the three energy bands of [20,57]keV, [57,82]keV, and [82,120]keV, respectively, and the display windows are [0.18,0.32]cm -1 、[0.15,0.23]cm -1 and [0.13,0.21]cm -1 .

[0194] The Gammex phantom was scanned with a photon counting CT machine at a tube current of 300 mA, including three energy bands: [20,55]keV, [55,80]keV, and [80,140]keV. The image reconstruction was performed using the traditional FBP algorithm, the FONT algorithm based on tensor representation, and the energy spectrum CT image reconstruction method proposed in this embodiment, respectively, to obtain the following: Figure 4The results are shown in Figure 1. The first row is the reconstructed image of FBP, the second row is the reconstructed image of FONT, and the third row is the reconstructed image of the proposed image reconstruction method. The first to third columns correspond to the three energy ranges of [20,55]keV, [55,80]keV, and [80,140]keV, and the display windows are [840,1340]HU, [820,1320]HU, and [800,1300]HU, respectively.

[0195] from Figure 3 and Figure 4 It can be seen that compared with the traditional FBP algorithm, the proposed image reconstruction method significantly suppresses quantum noise, effectively weakens the stripe artifacts caused by photon starvation, and fully improves the reconstruction quality of energy spectrum images (especially low-energy images). On the other hand, compared with the FONT algorithm based on tensor representation, the proposed image reconstruction method uses quaternion representation, that is, the quaternion matrix is ​​used to completely represent the entire 3D image block that should have low rank. The singular value decomposition and subsequent singular value compression operations can be directly performed on the quaternion matrix without expanding the quaternion matrix into a 2D real matrix. This makes up for the defect that the FONT algorithm destroys its internal structure when expanding the tensor into a matrix in all directions, better preserves the integrity of the multidimensional data structure, avoids the loss of useful structural information, and makes the reconstructed image higher in quality.

[0196] Example 3

[0197] Figure 7 The present invention is a schematic diagram of the structure of an electronic device provided in Embodiment 3 of the present invention. The electronic device includes a memory, a processor, and a computer program stored in the memory and executable on the processor, and the image reconstruction method of spectral CT of Embodiment 1 is implemented when the processor executes the computer program. Figure 7 The electronic device 30 shown is only an example and should not bring any limitation to the functions and scope of use of the embodiments of the present invention.

[0198] like Figure 7 As shown, the electronic device 30 may be in the form of a general-purpose computing device, for example, it may be a server device. The components of the electronic device 30 may include, but are not limited to: at least one processor 31, at least one memory 32, and a bus 33 connecting different system components (including the memory 32 and the processor 31).

[0199] The bus 33 includes a data bus, an address bus, and a control bus.

[0200] The memory 32 may include a volatile memory, such as a random access memory (RAM) 321 and / or a cache memory 322 , and may further include a read-only memory (ROM) 323 .

[0201] The memory 32 may also include a program / utility 325 having a set (at least one) of program modules 324, such program modules 324 including but not limited to: an operating system, one or more application programs, other program modules, and program data, each of which or some combination may include an implementation of a network environment.

[0202] The processor 31 executes various functional applications and data processing by running the computer programs stored in the memory 32, such as the image reconstruction method of spectral CT in embodiment 1 of the present invention.

[0203] The electronic device 30 may also communicate with one or more external devices 34 (e.g., keyboards, pointing devices, etc.). Such communication may be performed via an input / output (I / O) interface 35. Furthermore, the model-generated device 30 may also communicate with one or more networks (e.g., a local area network (LAN), a wide area network (WAN), and / or a public network, such as the Internet) via a network adapter 36. As shown, the network adapter 36 communicates with other modules of the model-generated device 30 via a bus 33. It should be understood that, although not shown in the figure, other hardware and / or software modules may be used in conjunction with the model-generated device 30, including but not limited to: microcode, device drivers, redundant processors, external disk drive arrays, RAID (RAID) systems, tape drives, and data backup storage systems.

[0204] It should be noted that although several units / modules or sub-units / modules of the electronic device are mentioned in the above detailed description, this division is merely exemplary and not mandatory. In fact, according to an embodiment of the present invention, the features and functions of two or more units / modules described above can be embodied in one unit / module. Conversely, the features and functions of one unit / module described above can be further divided into multiple units / modules to be embodied.

[0205] Example 4

[0206] This embodiment provides a computer-readable storage medium on which a computer program is stored. When the computer program is executed by a processor, the image reconstruction method of spectral CT of Embodiment 1 is implemented.

[0207] The readable storage medium may include but is not limited to: a portable disk, a hard disk, a random access memory, a read-only memory, an erasable programmable read-only memory, an optical storage device, a magnetic storage device or any suitable combination of the above.

[0208] In a possible implementation manner, the present invention may also be implemented in the form of a program product, which includes a program code. When the program product is run on a terminal device, the program code is used to enable the terminal device to execute the image reconstruction method of the energy spectrum CT of Example 1.

[0209] The program code for executing the present invention may be written in any combination of one or more programming languages, and may be executed entirely on a user device, partially on a user device, as an independent software package, partially on a user device and partially on a remote device, or entirely on a remote device.

[0210] Although the specific embodiments of the present invention are described above, it should be understood by those skilled in the art that this is only for illustration and the protection scope of the present invention is defined by the appended claims. Those skilled in the art may make various changes or modifications to these embodiments without departing from the principles and essence of the present invention, but these changes and modifications all fall within the protection scope of the present invention.

Claims

1. A method for image reconstruction of spectral CT. It is characterized in that The image reconstruction method comprises: Acquiring energy spectrum projection data, and inputting the energy spectrum projection data into an image reconstruction model to obtain an energy spectrum image to be updated; the image reconstruction model includes a regularization term; Extracting similar image blocks from the energy spectrum image to be updated to construct a hypercomplex matrix; the similar image blocks are image blocks with similar structures and containing different energy bands; The regularization term of the image reconstruction model is updated based on the hypercomplex matrix, and an updated energy spectrum image is obtained based on the updated image reconstruction model.

2. The image reconstruction method according to claim 1, It is characterized in that The steps of updating the image reconstruction model based on the hypercomplex matrix and obtaining an updated energy spectrum image based on the updated image reconstruction model include: Iteratively updating the image reconstruction model, that is, updating the image reconstruction model in each iteration based on the hypercomplex matrix of the energy spectrum image reconstructed last time; After the step of iteratively updating the image reconstruction model, the image reconstruction method further comprises: When the iterative updating process satisfies the preset convergence condition, the iterative reconstruction is stopped and the current energy spectrum image is determined as the final energy spectrum image.

3. The image reconstruction method according to claim 1, It is characterized in that The regularization term is the sum of several sub-regularization terms, and the sub-regularization terms include the sum of the ranks of the hypercomplex matrix.

4. The image reconstruction method according to claim 1, It is characterized in that The hypercomplex matrix includes a quaternion matrix or an octonion matrix; The step of extracting similar image blocks from the energy spectrum image to be updated to construct a hypercomplex matrix includes: Cutting the energy spectrum image to be updated into a plurality of first three-dimensional image blocks; Searching for a set of first three-dimensional image blocks of the most similar preset number for each of the first three-dimensional image blocks; Stretching and stacking each group of the most similar first three-dimensional image blocks to obtain a second three-dimensional image block; The energy segments of the second three-dimensional image block are substituted into the real part and / or the imaginary part of the hyper-complex matrix respectively to obtain a plurality of hyper-complex matrices.

5. The image reconstruction method according to claim 2, It is characterized in that The iterative update process satisfies a preset convergence condition including that the number of iterative updates reaches a preset threshold of the number of iterations; and / or, The steps of obtaining an updated energy spectrum image based on the updated image reconstruction model include: An alternating minimization algorithm is used to solve the minimization problem of the image reconstruction model to obtain an updated energy spectrum image.

6. The image reconstruction method according to claim 3, It is characterized in that After the step of updating the regularization term of the image reconstruction model based on the hypercomplex matrix, the image reconstruction method further includes: The rank function of the regularization term is approximated by using a substitute function of the rank function; the rank function is used to characterize the sum of the ranks of the hypercomplex matrix.

7. An image reconstruction system for spectral CT, It is characterized in that The image reconstruction system comprises: An energy spectrum projection data acquisition module, used for acquiring energy spectrum projection data, and inputting the energy spectrum projection data into an image reconstruction model to obtain an energy spectrum image to be updated; the image reconstruction model includes a regularization term; A hypercomplex matrix construction module, used for extracting similar image blocks from the energy spectrum image to be updated to construct a hypercomplex matrix; the similar image blocks are image blocks with similar structures and containing different energy bands; A reconstruction model updating module is used to update the regularization term of the image reconstruction model based on the hypercomplex matrix, and obtain an updated energy spectrum image based on the updated image reconstruction model.

8. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, It is characterized in that When the processor executes the computer program, the image reconstruction method for spectral CT according to any one of claims 1 to 6 is implemented.

9. A computer-readable storage medium having a computer program stored thereon, It is characterized in that When the computer program is executed by a processor, the image reconstruction method of spectral CT according to any one of claims 1 to 6 is implemented.

Images