Tomography method based on multi-grid division and pixel value dynamic constraint

Abstract

The invention provides a tomography method based on multi-grid division and pixel value dynamic constraint. The tomography method comprises three main steps of basic parameter setting, static division constraint reconstruction and dynamic division constraint reconstruction. In static division constraint reconstruction, grid blocks which are tightly connected are constructed, multiple kinds of division of a grid to be reconstructed is achieved, a constraint that pixel values are equal is applied to each grid block, reconstruction is carried out to obtain the pixel value of each grid block, weighting processing is carried out on pixel distribution under each kind of division, and pixel distribution on the grid to be reconstructed is synthesized; in the dynamic division constraint reconstruction, different connected regions are dynamically divided according to the pixel value interval and connectivity on the to-be-reconstructed grid, the constraint that the pixel values are equal is applied to each connected region, and the pixel distribution on the to-be-reconstructed grid is solved. The invention provides an imaging method for converting high-resolution grid reconstruction into reconstruction of a plurality of low-resolution grids, the anti-noise capability of a tomography system can be improved, and the method has a wide application prospect.

Description

Technical Field

[0001] The invention proposes a tomography method based on multi-grid division and dynamic constraint of pixel values, belonging to the technical field of tomography. Background Art

[0002] Tunable Diode Laser Absorption Spectroscopy (TDLAS) is a non-contact combustion diagnostic technology with the advantages of fast response speed, high sensitivity and strong reliability. It is widely used in the monitoring of parameters such as combustion field temperature, component concentration and pressure. Combining laser absorption spectroscopy with tomography technology can realize the reconstruction of the two-dimensional distribution of combustion flow field temperature and concentration. Usually, the amount of projection data in the imaging area is far less than the number of grids. Reconstructing the distribution to be measured from projection data is a typical underdetermined problem. Theoretically, there are countless solutions. Improving the underdetermination of the inverse solution problem is the key to obtaining high-quality reconstructed distribution.

[0003] Tomographic reconstruction algorithms are mainly divided into two categories, analytical methods and iterative methods. Analytical methods directly calculate the tomographic image of the object from the projection data based on mathematical theory without the need for iterative optimization. Iterative methods approximate the optimal solution through a series of step-by-step optimizations instead of directly using closed-form formulas.

[0004] The analytical algorithms mainly include Abel inverse transform and filtered back projection (FBP) algorithm, which are relatively simple to implement, highly efficient in calculation, and can provide stable reconstruction results under ideal conditions. In 2018, Zhang et al. published a paper titled "Reconstruction and simulation of temperature and CO2 concentration in an axisymmetric flame based on TDLAS" in Optik, Vol. 170, pp. 166-177. Under the assumption of axisymmetric distribution, the Fourier analysis Abel inverse transform algorithm was used to reconstruct the two-dimensional distribution of the attenuation coefficient. The reconstructed distribution was in good agreement with the computational fluid dynamics (CFD) simulation results. In 2014, Guha et al. published a paper titled "Tomographic laser absorption spectroscopy using Tikhonov regularization" in Applied Optics, Vol. 53, No. 34, pp. 8095-8103. They measured the material concentration and temperature distribution of the axisymmetric McKenna burner plume and reconstructed the radial profile of the absorption data using the Tikhonov regularized Abel inversion method, which suppressed the experimental noise amplification common in high spatial resolution reconstruction. However, the inverse Abel transform is only applicable to distributions with axisymmetric characteristics, but in most cases, it is difficult for the distribution to meet this condition, and the application scenario is relatively narrow. In 2024, Zhou et al. published a paper titled "High-Resolution TDLAS Tomography of Gaseous Temperature and H2O Concentration in Steady Flames" in Volume 73 of IEEE Transactions on Instrumentation and Measurement. They designed a reflective tunable diode laser absorption spectroscopy tomography system with an imaging time resolution of 50 milliseconds and a spatial resolution of 1.625 mm. The steady-state flame of the McKenna burner was measured. Each projection angle contained 200 laser paths, and the FBP algorithm was used to reconstruct a 100*100 pixel space, which was slightly different from the thermocouple measurement results.The FBP method is derived from the inverse Radon transform and its variants. It works well when there are enough projection data, but the reconstruction quality decreases when the available projection data is limited. Analytical algorithms are particularly suitable for situations that require real-time processing or limited hardware resources due to their speed and ease of implementation. However, when faced with non-ideal situations such as noisy data or incomplete projection data, the reconstruction results of this type of algorithm are worse than those of iterative algorithms.

[0005] Typical iterative algorithms include Algebra Reconstruction Technique (ART) and Simultaneous Algebra Reconstruction Technique (SART). These algorithms can use prior information to improve reconstruction quality and reduce artifacts in the presence of noise or incomplete data, and are more effective when reconstructing low-resolution distributions. In 2015, Peng et al. published a paper titled "Research on reconstruction algorithms for 2D temperature field based on TDLAS" in Volume 9674 of AOPC 2015: Optical and Optoelectronic Sensing and Imaging Technology. They simulated and reconstructed single-peak and double-peak temperature phantoms under different conditions with a grid of 10*10, and compared two different algorithms, ART and simulated annealing. The reconstruction quality of ART was mainly affected by ray coverage. The maximum deviation of the single-peak temperature phantom was 5.9%, while the maximum deviation of the double-peak temperature field was 25%. In 2024, Wang et al. published a paper titled "Vertical Distribution Mapping for Methane Fugitive Emissions Using Laser Path-Integral Sensing in Non-Cooperative Open Path ... In the paper, a method for mapping the vertical distribution of methane was proposed, which combined non-cooperative open-path laser path integral sensing with computer-aided tomography technology. The projection data was obtained by mobile laser scanning, and the SART algorithm was used to reconstruct a 10*10 grid. The results showed that for different methane gas bag schemes, the reconstruction error of the maximum concentration was consistently around 0.05, and the position error of the maximum concentration was between 0.01 and 0.025. In general, it is difficult to increase the projection data significantly due to limitations in hardware costs and other aspects. With the increase in the resolution requirements of the distribution to be measured, the above algorithm has poor reconstruction effect when the amount of projection data is far less than the number of unknown quantities to be calculated, and obvious artifacts will appear. When used for high-resolution reconstruction, it is necessary to make improvements on the basis of the basic iterative algorithm to alleviate the underdetermination and pathological nature of the solution problem.

[0006] There are three main ways to improve the underdetermination of the inverse problem: one is to increase the constraints of the original problem, the second is to reduce the unknowns to be solved, and the third is to solve an approximate problem of the original problem.

[0007] In terms of adding constraints to the original problem, the amount of projection data can be directly increased by increasing the number of optical paths, but this method is restricted by hardware conditions and it is difficult to achieve an order of magnitude increase; scanning optical paths can also increase the amount of projection data, but it will lead to a decrease in the system's measurement dynamic performance and measurement frame rate, which is suitable for steady-state or low-frequency changing distributions to be measured, but not suitable for objects with high-frequency changes; secondly, the prior information of the distribution to be measured can also be converted into constraints for the problem, such as the axial symmetry of the distribution, but for other distributions that do not conform to this information, the constraints need to be re-determined.

[0008] In terms of reducing the unknown quantities to be solved, one idea is to reduce the dimension of the problem by expressing the complex distribution as a linear combination of a finite number of basis functions, thereby reducing the number of parameters that need to be estimated: In 2022, Gao et al. published a paper titled "Sparse Zernike Fitting for Dynamic LAS Tomographic Images of Temperature and Water Vapor Concentration" in Volume 71 of IEEE Transactions on Instrumentation and Measurement. They proposed a sparse Zernike fitting method based on the characteristics of the smooth temperature distribution in real physical reality in tomographic imaging. They used a two-dimensional Zernike polynomial to fit the distribution to be reconstructed, and converted the calculated pixel values ​​into calculated polynomial coefficients, effectively reducing the number of unknown quantities and alleviating the underdetermination of the problem; in 2022, Zhou et al. published a paper titled "Optimized CT-TDLAS reconstruction performance evaluation of least squares with the optimized polynomial fitting" in Volume 10 of Frontiers in Physics. In the polynomial-fitting method, a two-dimensional polynomial is used to fit the distribution, and the least squares method is used to estimate the polynomial coefficients. The temperature distribution of the Bunsen burner flame is measured under the conditions of 16 horizontal and 16 vertical projection data. The SART algorithm, Tikhonov regularization algorithm and the proposed method are compared. The polynomial fitting method has better reconstruction quality. The fitting accuracy of this method depends on the type of polynomial used and the number of polynomials. When the properties of the measured distribution and the polynomial basis do not match well or the distribution structure is complex, the reconstructed distribution deviation is large. In addition, the use of globally defined basis functions may lead to poor performance of local features in the reconstructed image.By discretizing the value range of pixel values, the amount of data to be processed can be reduced, and the complexity of the problem can also be effectively reduced: In the paper "A Binary Valued Reconstruction Algorithm for Discrete TDLAS Tomography of Dynamic Flames" published by Qiu et al. in Volume 72 of IEEE Transactions on Instrumentation and Measurement in 2023, the continuously distributed gas temperature and concentration in the target combustion area are discretized into two sets of finite value pairs within the feasible range, and the problem is transformed into solving a 0-1 binary matrix. However, discretization is essentially an approximate method, which will lead to information loss, introduce additional deviations, affect the quality of the reconstructed image, and the reconstruction quality also depends on the reasonable selection of discrete values.

[0009] In terms of solving the approximate problem of the original problem, it is generally considered to add a regularization term in the optimization so that the calculated solution is close to the solution of the original problem, such as TV regularization, Tikhonov regularization, etc. By introducing additional constraints to limit the solution space, it can effectively deal with problems such as noise and incomplete data. In 2018, Bao et al. published a paper titled "X-Rays Tomographic Reconstruction Images using Proximal Methods based on L1 Norm and TV Regularization" in Procedia Computer Science, Vol. 127, pp. 236-245. They proposed a total variation algorithm based on the L1 norm and proximal function, and applied it to the nondestructive evaluation of materials using two-dimensional reconstruction of X-ray tomographic images containing real defects. In 2021, Bao et al. published a paper titled "Relative Entropy Regularized TDLAS Tomography for Robust Temperature and Water Vapor Concentration" in IEEE Transactions on Instrumentation and Measurement, Vol. 70. In the recent paper, a relative entropy regularization method was introduced, which significantly improved the quality of tomographic temperature images compared with the existing synchronous algebraic reconstruction technology and had good robustness to TDLAS tomographic measurement noise.In 2023, Zhang et al. published a paper titled "Tikhonov regularization-based extended Kalman filter technique for robust and accurate reconstruction in diffuse optical tomography" in Journal of the Optical Society of America A, Vol. 40, No. 1, pp. 10-20. They combined Tikhonov regularization with the extended Kalman filter method to reduce the instability of the observation matrix in the tomography problem, conducted a series of numerical simulations and imaging experiments, and quantitatively compared the experimental results with the algebraic reconstruction technique and the Levenberg-Marquardt algorithm. The proposed algorithm can accurately estimate the position and size of the target, and the imaging accuracy and noise robustness are significantly improved. In 2021, Wei et al. published a paper titled "Tikhonov regularization-based extended Kalman filter technique for robust and accurate reconstruction in diffuse optical tomography" in Combustion and Flame. In the paper "Volumetric laser absorption imaging of temperature, CO and CO inlaminar flames using 3D masked Tikhonov regularization" published in Volume 224, Pages 239-247 of the Journal of Flame, the linear tomography method with Tikhonov regularization is used to extend the mid-infrared laser absorption imaging of flame temperature and substance concentration to three-dimensional space. The spatial convolution of two small-scale (less than centimeter) laminar Bunsen flames is measured and compared with the reference image based on Abel inversion of a single axisymmetric flame. The results show that the introduction of 3D regularization and masking improves the spatial resolution of gradients within the flow field, improves the overall reconstruction accuracy, and effectively reduces the number of required projection angles. The introduction of regularization terms usually makes the objective function more complicated, resulting in an increase in computational cost. Moreover, regularization is usually based on some prior knowledge or assumptions (such as smoothness, sparsity, etc.). If the actual distribution does not meet these assumptions, regularization may not bring the expected effect, or even deteriorate the reconstruction quality. The effect of regularization also depends on the choice of regularization parameters.

[0010] In addition, neural network is also one of the inverse solution methods for tomography. It can effectively use limited or incomplete projection data for reconstruction, and can even recover more accurate results from sparse sampling. After the network training is completed, the neural network can perform rapid reconstruction with good real-time performance. In 2021, Si et al. published a paper titled "Hierarchical Temperature Imaging Using Pseudoinversed Convolutional Neural Network Aided TDLAS Tomography" in Volume 70 of IEEE Transactions on Instrumentation and Measurement. They proposed a pseudo-inverse convolutional neural network for hierarchical temperature imaging. Compared with traditional convolutional neural networks, it has higher accuracy and computational efficiency. The effectiveness of the algorithm was verified by numerical simulation and experiments. In 2022, Si et al. published a paper titled "A Quality-Hierarchical Temperature Imaging Network for TDLAS Tomography" in Volume 71 of IEEE Transactions on Instrumentation and Measurement. They proposed a TDLAS tomography quality hierarchical network based on stacked long short-term memory. The network outputs two reconstructed temperature images, low-quality images and high-quality images, which correspond to different numbers of network layers and have different computational costs. Low-quality images provide more timely temperature reconstruction, which can meet the real-time dynamic monitoring of turbulence-chemical interactions. The time resolution reaches tens of thousands of frames per second. High-quality images can be stored and used for offline analysis and diagnosis, and the temperature reconstruction is more accurate. In 2024, Shi et al. published a paper titled "Multi-View Synthesis of Sparse Projection of Absorption Spectra Based on Joint GRUand U-Net" in Volume 14, Issue 9 of Applied Sciences. A sparse projection multi-view synthesis model based on U-Net was proposed. This model combines the sequence learning characteristics of the gated recurrent unit and the generalization ability of the residual network, and reconstructs the temperature field of sparse projection and synthetic projection respectively, indicating that the application of this model can significantly improve the reconstruction accuracy of high-energy combustion temperature field.Neural network methods usually require a large amount of labeled data for training. Even after sufficient training, neural networks may perform poorly on unseen distribution types or conditions and have limited generalization capabilities. Moreover, neural networks are "black box" models with poor interpretability.

[0011] Based on the above background, the present invention proposes a tomographic imaging method based on multi-grid division and dynamic constraints on pixel values, which realizes the distribution reconstruction of high-resolution grids under the condition of less projection data through static multi-grid division pixel value constraint reconstruction and dynamic multi-grid division pixel value constraint reconstruction. This method does not rely on the prior characteristics of the measured distribution, can greatly reduce the number of unknown quantities to be determined, reduce the underdetermination and pathological nature of the reconstruction problem, and realize adaptive dynamic adjustment of pixel value constraints to obtain a high-resolution pixel value distribution with fewer artifacts. This method can not only be used for laser absorption spectral tomography, but also can be used for other parameter distribution tomography that meets the path integral principle, and has broad application prospects. Summary of the invention

[0012] The present invention proposes a tomographic imaging method based on multi-grid division and dynamic constraint of pixel values, including three main steps: basic parameter setting, pixel value constraint reconstruction of static multi-grid division, and pixel value constraint reconstruction of dynamic multi-grid division. In static division constraint reconstruction, closely connected grid blocks are constructed to realize multiple divisions of the grid to be reconstructed, and pixel value equality constraints are imposed in each grid block. Reconstruction is performed to obtain the pixel value of each grid block, and the pixel distribution under each division is weighted to synthesize the pixel distribution on the grid to be reconstructed. In dynamic division constraint reconstruction, the grid to be reconstructed is dynamically divided into different connected areas according to the pixel value interval and connectivity, and pixel value equality constraints are imposed in each connected area to obtain the pixel distribution on the grid to be reconstructed. It is implemented specifically according to the following steps:

[0013] Step 1: Basic parameter settings.

[0014] Determine the number of grids for 2D distribution reconstruction The pixel set to be reconstructed B = {(1,1), (1,2), ..., (N,N)}, where (i,j) represents the pixel in the i-th row and j-th column in the grid; determine n L Projection paths and projection data Calculate the sensitivity matrix Among them, the item in the i-th row and the j-th column represents the projection path L i Located in pixels The length inside, It means rounding down, and the sensitivity matrix is ​​also expressed as W = [w (1,1) w (1,2) …w (N,N) ],in, is a column vector; the distribution of the parameter to be measured is expressed as pixel the parameter value at is expressed as α(i,j); set the scale of the divided grid block to M*M( M < N), the iterative condition error limit ε, and set the initial value of the number of iterations k to 1.

[0015] Step 2: Static partition constraint reconstruction: Use static multi-grid partitioning to divide the grid to be reconstructed into multiple grid blocks, add the constraint that the pixel values of the pixels belonging to the same grid block are equal, update the sensitivity matrix and the unknowns to be solved, and perform distribution reconstruction, and synthesize the weights of the reconstruction distributions corresponding to different partitioning methods.

[0016] The static multi-grid partitioning method adopted is to completely cover the N*N grid with non-overlapping grid blocks of size M*M (M < N) that are adjacent to each other, and satisfy:

[0017]

[0018] where q (i,j) is the set of pixels contained in an M*N grid block, uniquely determined by (i,j), is the set of M*M grid blocks, representing a covering method, uniquely determined by (i0,j0).

[0019] Traverse (i0,j0), the total number of covering methods determined by formula (1) is M 2 types, which are respectively the elements in; each covering method defines a partition of the set B of pixels to be reconstructed:

[0020]

[0021] The total number of partitions of the set B determined by formula (2) is M 2 types, and the t-th partition is also expressed as 1 ≤ t ≤ M 2 , where is the m-th element in the set S (t) represents the m-th grid block of the partition of S (t) , n t represents the number of elements in the set S (t) .

[0022] The constraint L (t) constructed according to the partition S (t) is expressed as:

[0023]

[0024] in, Indicates L (t) For grid blocks The constraints on the pixels within the same set make them belong to the same All pixels in the same pixel value have the same value.

[0025] In the constraint L (t) The simplified reconstruction model with low number of unknowns under the action is:

[0026] W (t) X (t) =P (4)

[0027] Among them, W (t) and X (t) Respectively represent the simplified sensitivity matrix and the simplified unknown quantity to be solved, satisfying:

[0028] 1≤m≤n t , W (t) The mth column vector of , T represents the transpose of the matrix, and P is the projection data.

[0029] In M 2 Under the condition that each constraint is valid separately, according to formula (4), we can get M 2 Different distributions f t , 1≤t≤M 2 .

[0030] The weighted synthesis method used for the reconstruction distribution corresponding to different partitions is to calculate M respectively. 2 The projection error of the reconstructed distribution is determined according to the projection error, and then weighted synthesis is performed:

[0031]

[0032] in, is the distribution f t , ‖·‖2 represents the vector 2-norm, and F0 is the synthetic reconstruction distribution, which is used as the initial value of the iteration to enter the dynamic partition constraint reconstruction.

[0033] Step 3: Dynamic partitioning and constrained reconstruction: In each iteration, pixels are classified according to the pixel value interval of the current reconstructed high-resolution distribution, and the grid to be reconstructed is divided into multiple connected regions according to pixel category and connectivity. Pixel value equality constraints are added to pixels belonging to the same connected region, and the distribution is solved by a simplified reconstruction model. The iteration stopping condition is determined, and the iteration is exited when the stopping condition is met, and the reconstructed distribution is output.

[0034] The pixels are classified as follows: for the reconstructed distribution F output at the k-1th iteration k-1 The minimum and maximum values ​​of all pixel values ​​in this distribution are recorded as A k-1 and B k-1 , the interval [A k-1 ,B k-1 ] is divided into h k non-overlapping intervals of equal length r j (0≤j <h k , ), h k is the number of equal-length intervals divided, which is monotonically non-decreasing as k increases; if the pixel value of a pixel belongs to interval r j , then the pixel is said to belong to the jth class.

[0035] In the kth iteration, the connected region partitioning method used is to record the maximum set of pixels belonging to the same class and closely connected in the N*N grid as a connected region, which is expressed as Among them, d k It is the total number of connected areas divided in the N*N grid, and the adjacency relationship adopts 4-adjacency relationship.

[0036] Constraints C constructed based on the divided connected regions (k) It is expressed as:

[0037]

[0038] in, Represents C (k) For connected areas The constraints on the pixels within the same set make them belong to the same All pixels in the same pixel value have the same value.

[0039] In the constraint C (k) The simplified reconstruction model with low number of unknowns under the action is:

[0040] W (k) X (k) =P (7)

[0041] Among them, W (k) and X (k) Respectively represent the simplified sensitivity matrix and the simplified unknown quantity to be solved, satisfying: 1≤m≤d k , W (k) The mth column vector of , P is the projection data.

[0042] In the constraint C(k) Under these conditions, the distribution F is calculated according to formula (7): k .

[0043] The way to determine the end of iteration is to determine the change amount Is it true? If so, output the final reconstructed distribution F k , if not, set k=k+1 and return to step 3.

[0044] Beneficial effects of the present invention: The static partitioning constrained reconstruction proposed by the present invention has nothing to do with the distribution to be measured, but only depends on the spatial distribution of the optical path and the grid to be reconstructed, and can effectively extract high-resolution information of the distribution to be measured, with a wide range of applications and strong versatility. The dynamic partitioning constrained reconstruction proposed by the present invention constructs dynamic pixel value constraints based on the reconstructed high-resolution distribution, does not require prior information, and iteratively realizes the adaptive extraction of distribution features. In addition, both the static partitioning constrained reconstruction and the dynamic partitioning constrained reconstruction proposed by the present invention greatly reduce the number of unknown quantities to be solved in the reconstruction problem, reduce the underdetermination of the problem, and do not limit the specific reconstruction solution algorithm, and have high scalability. The present invention provides an imaging method that converts high-resolution grid reconstruction into multiple low-resolution grid reconstructions, which can realize the distribution reconstruction of high-resolution grids under the condition of insufficient projection data, can improve the noise resistance of the tomography imaging system, and has broad application prospects. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 It is a flow chart of a tomography method based on multi-grid partitioning and dynamic constraint of pixel values.

[0046] Figure 2 is the parameter distribution within the set imaging area.

[0047] Figure 3 It is a weighted synthesis of the reconstruction distributions corresponding to 25 static multi-gridding constraints.

[0048] Figure 4 It is the number of unknown quantities to be solved in the iterative process and the average relative error of the reconstructed distribution.

[0049] Figure 5 is the final reconstructed distribution of the iteration. DETAILED DESCRIPTION

[0050] The present invention will be further described below in conjunction with examples.

[0051] This example uses the optical path layout of a five-angle fan-shaped laser absorption spectral tomography system to perform distribution reconstruction calculations. Among them, the fan-shaped laser source is located at the five vertices of a regular pentagon, and the detectors are located on the five sides of the regular pentagon. There are 12 detectors arranged at equal intervals on each side. The side length of the regular pentagon is 210mm. The closest distance between each detector and the vertex of the regular pentagon is 26.02mm, and the distance between two adjacent detectors is 14.36mm. Assume that the light emitted by each fan-shaped laser source can be received by 24 detectors on the opposite side and half of its left and right adjacent sides, and the five laser sources generate a total of 120 projection rays. The specific implementation steps of distribution reconstruction are as follows:

[0052] Step 1: Basic parameter settings.

[0053] Set N = 100, grid the 100mm*100mm area to be reconstructed in the center of the regular pentagon into a 100*100 pixel grid, determine the equations of 120 projection rays, calculate the distance each ray passes through each pixel, and generate the sensitivity matrix W 120×10000 . Set the absorption rate distribution equation to:

[0054]

[0055] The parameter value of each grid is determined by the coordinates of the grid center point (x p ,y p ) corresponds to α(x p ,y p ) function value, and the original parameter distribution is as shown in the attached Figure 2 As shown, it is represented as a column vector X 10000×1 , calculate the projection data P 120×1 =W 120×10000 X 10000×1 Set M = 5, that is, the grid block size is 5*5, and the iteration error limit ε = 10 -3 , number of iterations k = 1.

[0056] Step 2: Static partition constraint reconstruction.

[0057] There are 25 types of coverage determined by formula (1), namely Q (1,1) ,Q (1,2) ,…,Q (1,5) ,Q (2,1) ,Q (2,2) ,…,Q (2,5) ,…,Q (5,1) ,Q (5,2) ,…,Q (5,5) For the construction element, each cover defines a partition, a total of 25 partitions S (1) ,S (2) ,…,S (25) , corresponding to 25 constraints L(1) ,L (2) ,…,L (25) . Using the SART algorithm as the distribution reconstruction algorithm, we reconstruct under the condition that each constraint acts separately and obtain f1,f2,…,f 25 , calculate the projection error of each distribution.

[0058] The following table shows the set S for all partitions. (t) The number of elements n t , the projection error of the reconstructed distribution β t , where n t Numerically equal to the number of unknowns after adding constraints.

[0059]

[0060] In the 25 partitioning cases, the added constraints reduced the number of original unknowns from 10,000 to 400, 420, and 441, a reduction of 96.00%, 95.80%, and 95.59%, respectively.

[0061] The distribution of the synthesis is shown in the attached Figure 3 As shown, the projection error of the distribution is 24.06%. Compared with the original distribution, the average relative error is 48.88%, and the image similarity SSIM is 0.8558.

[0062] Step 3: Dynamic partition constraint reconstruction.

[0063] Use the distribution synthesized in step 2 as the initial value of the distribution iteration and set h k =min(10+3k,50), using SART algorithm as the distribution reconstruction algorithm. The number of connected regions divided during the iteration is d k (also the number of unknown quantities to be determined), the average relative error between the reconstructed distribution and the original distribution is shown in the attached Figure 4 The added pixel value constraint reduces the number of original unknowns from 10,000 to about 1,100, a reduction of about 89%, greatly reducing the underdetermination of the problem to be solved.

[0064] After 50 iterations, the final reconstructed distribution is as shown in the attached figure. Figure 5 As shown, compared with the original distribution, the average relative error of the reconstructed distribution is 13.18%, and the image similarity SSIM is 0.9878, which shows that the present invention can reconstruct a high-resolution parameter distribution with a certain quality under the condition of relatively low projection data.

[0065] The above description of the present invention and its embodiments is not limited thereto, and the accompanying drawings are only one embodiment of the present invention. Without departing from the purpose of the present invention, any structure or embodiment similar to the technical solution designed without creativity shall fall within the protection scope of the present invention.

Claims

1. A tomography method based on multi-grid partitioning and dynamic pixel value constraint, which includes three main steps: basic parameter setting, static partitioning constraint reconstruction, and dynamic partitioning constraint reconstruction; in the static partitioning constraint reconstruction, closely adjacent grid blocks are constructed to achieve multiple partitions of the grid to be reconstructed, the constraint of equal pixel values is imposed within each grid block, reconstruction is performed to obtain the pixel values of each grid block, the pixel distributions under each partition are weighted, and the pixel distribution on the grid to be reconstructed is synthesized; in the dynamic partitioning constraint reconstruction, according to the pixel value range and connectivity on the grid to be reconstructed, it is dynamically divided into different connected regions, the constraint of equal pixel values is imposed within each connected region, and the pixel distribution on the grid to be reconstructed is obtained by solving, and it is specifically implemented according to the following steps: Step 1. Basic parameter setting: Determine the number of grids for 2D distribution reconstruction The pixel set to be reconstructed B = {(1,1), (1,2), ..., (N,N)}, where: (i,j) represents the pixel in the i-th row and j-th column of the grid; determine n L Projection paths and projection data Calculate the sensitivity matrix Among them, the item in the i-th row and the j-th column represents the projection path L i Located in pixels The length inside, It means rounding down, and the sensitivity matrix is ​​also expressed as W = [w (1,1) w (1,2) …w (N,N) ],in, is a column vector; the distribution of the parameter to be measured is expressed as Pixel The parameter value at is expressed as α(i,j); set the grid block size The iteration condition error limit ε, the number of iterations k is set to 1; Step 2. Static partitioning constraint reconstruction: Construct closely adjacent grid blocks to achieve multiple partitions of the grid to be reconstructed, impose the constraint of equal pixel values within each grid block, perform reconstruction to obtain the pixel values of each grid block, weight the pixel distributions under each partition, and synthesize the pixel distribution on the grid to be reconstructed; Step 3. Dynamic partitioning constraint reconstruction: In each iteration, according to the pixel value range and connectivity on the grid to be reconstructed, the grid to be reconstructed is dynamically divided into different connected regions, the constraint of equal pixel values is imposed within each connected region, and the pixel distribution on the grid to be reconstructed is obtained by solving.

2. A tomographic imaging method based on multi-grid division and pixel value dynamic constraint according to claim 1, characterized in that: During static partitioning constraint reconstruction, the grid to be reconstructed is divided into multiple grid blocks by static multi-grid partitioning, the constraint of equal pixel values is added to the pixels belonging to the same grid block, the sensitivity matrix and the unknowns to be solved are updated and distribution reconstruction is performed, and the reconstruction distributions corresponding to different partitioning methods are weighted and synthesized: The adopted static multi-grid partitioning method is to completely cover the N*N grid with non-overlapping grid blocks of size M*M (M < N) that are adjacent to each other, and satisfy: Among them, q (i,j) is a set of pixels contained in an M*N grid block, uniquely determined by (i, j). It is a set of M*M grid blocks, representing a covering method, which is uniquely determined by (i0, j0); Traversing (i0, j0), the total number of coverage methods determined by formula (1) is M 2 Species, respectively Elements in; each coverage method Both define a partition of set B: The partitions of set B determined by formula (2) are M 2 The tth division is expressed as in, For the set S (t) The mth element in represents S (t) The mth grid block divided, n t Represents the set S (t) The total number of elements in According to the division S (t) The constructed constraint L (t) It is expressed as: in, Indicates L (t) For grid blocks The constraints on the pixels within the same set make them belong to the same The pixel values ​​of all pixels are the same; In the constraint L (t) The simplified reconstruction model with low number of unknowns under the action is: W (t) X (t) =P (4) Among them, W (t) and X (t) Respectively represent the simplified sensitivity matrix and the simplified unknown quantity to be solved, satisfying: W (t) The mth column vector of , T represents the transpose of the matrix, and P is the projection data; In M 2 Under the condition that each constraint is valid separately, according to formula (4), we can get M 2 Different distributions f t , 1≤t≤M 2 ; The weighted synthesis method used for the reconstruction distribution corresponding to different partitions is to calculate M respectively. 2 The projection error of the reconstructed distribution is determined according to the projection error, and then weighted synthesis is performed: in, For the distribution f t , ‖·‖2 represents the vector 2-norm, and F0 is the synthetic reconstruction distribution, which is used as the initial value of the iteration to enter the dynamic partition constraint reconstruction.

3. The tomography method based on multi-grid division and pixel value dynamic constraint according to claim 1, characterized in that: In the dynamic partitioning constraint reconstruction, in each iteration, the pixels are classified according to the pixel values on the current grid to be reconstructed, the grid to be reconstructed is divided into multiple connected regions by pixel classification and connectivity, the constraint of equal pixel values is added to the pixels belonging to the same connected region, the reconstruction model is simplified, and the pixel distribution on the grid to be reconstructed is obtained by solving: The pixel classification method is as follows: for the reconstructed distribution F output at the k-1th iteration k-1 The minimum and maximum values ​​of all pixel values ​​in this distribution are recorded as A k-1 and B k-1 , the interval [A k-1 ,B k-1 ] is divided into h k non-overlapping intervals of equal length h k is the number of equal-length intervals divided, which is monotonically non-decreasing as k increases; if the pixel value of a pixel belongs to interval r j , then the pixel is said to belong to the jth class; In the kth iteration, the connected region partitioning method used is to record the maximum set of pixels belonging to the same class and closely connected in the N*N grid as a connected region, which is expressed as Among them, d k is the total number of connected regions divided in the N*N grid, and the adjacency relationship adopts 4-adjacency relationship; Constraints C constructed based on the divided connected regions (k) It is expressed as: in, Represents C (k) For connected areas The constraints on the pixels within the same set make them belong to the same The pixel values ​​of all pixels are the same; In the constraint C (k) The simplified reconstruction model with low number of unknowns under the action is: W (k) X (k) =P (7) Among them, W (k) and X (k) Respectively represent the simplified sensitivity matrix and the simplified unknown quantity to be solved, satisfying: W (k) The mth column vector of , P is the projection data; In the constraint C (k) Under these conditions, the distribution F is calculated according to formula (7): k ; The way to determine the end of iteration is to determine the change amount Is it true? If so, output the final reconstructed distribution F k , if it does not hold, set k=k+1 and re-reconstruct the dynamic partitioning constraints.

Images