A fluorescence tomographic reconstruction imaging method based on a physical model deep network

Abstract

The application discloses a fluorescence tomography reconstruction imaging method based on a physical model depth network, acquires tissue surface light intensity information through a CCD camera, performs mesh partitioning on a phantom by using Amira software, respectively partitions two meshes with different precisions, respectively solves corresponding system matrices and adjacent matrices of mesh nodes by using the mesh after partitioning and optical parameters of each region in the phantom through a finite element method, sets a transfer matrix corresponding to the power density of different meshes, then transmits surface light intensity data into a depth network, solves in a coarse mesh, and demarcates a feasible area for solving a fine mesh according to the solving result of the coarse mesh, and solves in the feasible area, so that the norm of the obtained solution is small, the result is relatively sparse, the error of the norm and the result forward and the real surface light intensity is considered during solving, the position of a fluorescent group in the tissue can be positioned more accurately, and the application has important significance for studying tumors and early diagnosis of tumors.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of medical imaging, and particularly relates to a fluorescence tomography reconstruction imaging method based on a physical model deep network. BACKGROUND

[0002] Fluorescence tomography reconstruction refers to calculating the position and power of a fluorescence light source according to the light intensity information transmitted to the surface inside the tissue combined with the optical parameters of the tissue, including bioluminescence tomography, fluorescence molecular tomography, X-ray luminescence computed tomography, and Cherenkov luminescence tomography. These imaging technologies have the advantages of high sensitivity and high specificity, and can be applied to early diagnosis of cancer and the like.

[0003] Fluorescence tomography reconstruction inside the tissue is an ill-posed and underdetermined problem, and needs to solve the result with the minimum norm, and the interference has a greater impact on the solution result. The traditional method ([1] Zhang X, Badea C T, and Johnson GA. Three-dimensional reconstruction in free-space whole-body fluorescence tomography of mice using optically reconstructed surface and atlas anatomy. Journal of Biomedical Optics. 2009, 14(6). 0640101-0640109; [2] Han D, Tian J, et al. A fast reconstruction algorithm for fluorescence molecular tomography with sparsity regularization. Optics Express. 2010, 18(8). 8630-8646.) needs to demarcate a feasible region for solving, reduces the node number for solving and is convenient for solving, but needs manual operation. The method using deep learning (Yuan G, Kun W, et al. Nonmodel-based bioluminescence tomography using a machine-learning reconstruction strategy. Optica. 2018, 5(11).) needs to make a large amount of data set for training. SUMMARY

[0004] In order to overcome the above-mentioned prior art defects, the purpose of the present application is to provide a fluorescence tomography reconstruction imaging method based on a physical model deep network, which ensures that the norm of the obtained solution is small and the result is relatively sparse; the norm and the error of the result forward and the real surface light intensity are considered during solving; the position of the fluorescent group in the tissue can be more accurately located, which has important significance for studying tumors and their early diagnosis.

[0005] In order to achieve the above-mentioned purpose, the technical scheme adopted by the present application is:

[0006] A fluorescence tomography reconstruction imaging method based on a physical model deep network, comprising the following steps:

[0007] 1) Obtain the tissue surface light intensity data through a CCD camera, divide the phantom into a coarse grid and a fine grid by using Amira software; use the divided grid and the optical parameters of each region in the phantom to solve the coarse grid system matrix A1 and the fine grid system matrix A2 by the finite element method, and solve the coarse grid adjacency matrix W1 and the fine grid adjacency matrix W2 of each node of the grid; match the surface light intensity data to the two grids to form the coarse grid surface light intensity b1 and the fine grid surface light intensity b2 corresponding to the coarse grid and the fine grid;

[0008] 2) Set the coarse grid transfer matrix B1 and the fine grid transfer matrix B2 of the power density, the transfer matrix transfers the power density to the adjacent node position, and the column and the row of the transfer matrix should be 1;

[0009] 3) Use the coarse grid surface light intensity b1, the coarse grid system matrix A1 and the coarse grid adjacency matrix W1 to solve the coarse grid system matrix A1 as a physical model constraint in the deep network;

[0010] 4) Use the result obtained in step 3) to divide the feasible region in the fine grid, and solve by using the fine grid surface light intensity b2, the fine grid system matrix A2 and the fine grid adjacency matrix W2.

[0011] The step 1) is specifically:

[0012] 1.1) Divide the phantom into a grid by using Amira software;

[0013] 1.2) Discretize the diffusion equation according to the nodes of the grid to solve the system matrix;

[0014] The diffusion equation is:

[0015]

[0016] Wherein, D(r) is the diffusion coefficient, S(r) is the light source, Φ(r) is the light flux generated by the light source, μ a is the absorption coefficient, and q is the refractive index;

[0017] Discretization is:

[0018] KΦ = F

[0019] Where K is the stiffness matrix, Φ is the light flux matrix, and F is the light source matrix.

[0020] Then:

[0021] Φ = K -1 F

[0022] The system matrix A is the result of setting 0 to the corresponding rows of the internal nodes of Φ.

[0023] 1.3) Obtain the adjacency matrix W of all nodes in the grid.

[0024] The step 3) is specifically:

[0025] 3.1) The surface light intensity b1 of the coarse grid is transmitted into the deep network, and the output x of the deep network is the light source power density distribution, and the objective function is:

[0026] loss = mse(b1, A1x) + mse(b1, A1B1x) 2 x)

[0027] Adjust the internal parameters of the network by the mean square error of the forward fluorescence distribution and the real surface;

[0028] 3.2) Solve multiple times, record all solution results, and take the non-zero nodes after multiplying the sum of all results by the coarse grid adjacency matrix W1 as the feasible region;

[0029] 3.3) Continue to solve in the feasible region, and the objective function is:

[0030] loss = mse(b1, A1x) + mse(b1, A1B1x) + mse(b1, A1B1x) + mse(b1, A1B1x) 2 x) + mse(b1, A1B1x) 3 x)

[0031] The iteration of the feasible region is the non-zero node after multiplying the solution result x by the coarse grid adjacency matrix W1; record the loss value and result of all iteration results, when the program cannot produce a feasible region when solving the all-zero solution, select the result with the smallest loss value from the previous results to produce a feasible region, if all the loss values of the previous results are greater than the loss value of the all-zero solution, use the initial feasible region to continue solving; when the number of iterations is reached, select the result with the smallest loss value in the record as the solution result of this time;

[0032] 3.4) continue to solve the residual error of the solution in the same way, stop iteration when the value of the solution is less than a certain value; accumulate all the solutions and discard the nodes with smaller power density as the final coarse grid solution.

[0033] The input of the deep network is the measured tissue surface light intensity, and the output of the deep network is the light source distribution; the forward result corresponding to the given light source distribution is calculated through a system matrix, and the difference between the forward result and the measured tissue surface light intensity is taken as a loss function and transmitted to the neural network, and the parameters are learned through iteration.

[0034] The step 4) is specifically:

[0035] 4.1) according to the solution result of the coarse grid, taking the fine grid nodes within 1mm range of the non-0 value node coordinates as the solution feasible region of the fine grid;

[0036] 4.2) the fine grid surface light intensity b2 is transmitted into the deep network, and the output x of the deep network is the light source power density distribution, and the objective function is:

[0037] loss=mse(b2,A2x)+mse(b2,A2B2x)+mse(b2,A2B2 2 x)+mse(b2,A2B2 3 x)

[0038] The loss value and result of all iteration results are recorded, when the program cannot generate a feasible region when the all-zero solution is obtained, the result with the minimum loss is selected from the previous results to generate a feasible region, if the loss value of all previous results is greater than the loss value of the all-zero solution, the initial feasible region is used for further solving; when the number of iterations is reached, the result with the minimum loss in the record is selected as the solution result of this time;

[0039] 4.3) continue to solve the residual error of the solution in the same way, stop iteration when the value of the solution is less than a certain value; accumulate all the solutions and discard the nodes with smaller power density as the final solution.

[0040] The beneficial effects of the present application are:

[0041] Since the present application adopts a physical model as a constraint, the step of preparing a training set required by a data-driven deep network is reduced, a physical drive is added to the deep network, and the interpretability is increased;

[0042] Since the present application adopts the idea of gradually reducing the feasible region for solving, it is not necessary to manually define the feasible region, and the accuracy and robustness of the solution are increased. BRIEF DESCRIPTION OF DRAWINGS

[0043] Figure 1Flow chart of the method of the present application.

[0044] Figure 2 Deep network structure diagram of the embodiment of the present application.

[0045] Figure 3 Flow chart of the feasible region solving of the embodiment of the present application.

[0046] Figure 4 Result diagram of the embodiment of the present application. DETAILED DESCRIPTION

[0047] In order to make the objectives, technical solutions and advantages of the present application clearer, the present application is further described in detail below with reference to embodiments. The embodiments described herein are only used to explain the present application and do not limit the present application.

[0048] Referring to Figure 1 A fluorescence tomographic reconstruction imaging method based on a physical model deep network, comprising the following steps:

[0049] 1) Obtain the tissue surface light intensity data by a CCD camera, divide the phantom into a coarse grid and a fine grid by Amira software, and divide the grid into a coarse grid and a fine grid; use the divided grid and the optical parameters of each region in the phantom to solve the coarse grid system matrix A1 and the fine grid system matrix A2 by a finite element method, and solve the coarse grid adjacency matrix W1 and the fine grid adjacency matrix W2 of each node of the grid; match the surface light intensity data to the two grids to form the coarse grid surface light intensity b1 and the fine grid surface light intensity b2 corresponding to the coarse grid and the fine grid; specifically:

[0050] 1.1) Divide the phantom into a grid by Amira software;

[0051] 1.2) Discretize the diffusion equation according to the nodes of the grid to solve the system matrix;

[0052] The diffusion equation is:

[0053]

[0054] Wherein, D(r) is the diffusion coefficient, S(r) is the light source, Φ(r) is the light flux generated by the light source, μ a is the absorption coefficient, and q is the refractive index;

[0055] Discretization is:

[0056] KΦ=F

[0057] Wherein, K is the stiffness matrix, Φ is the light flux matrix, and F is the light source matrix;

[0058] Then:

[0059] Φ=K-1 F

[0060] System matrix A is the result of Φ internal node corresponding row 0;

[0061] 1.3) find all the nodes in the grid adjacent matrix W;

[0062] 2) set the power density of the coarse grid transfer matrix B1, fine grid transfer matrix B2, transfer matrix will be the power density to its adjacent node position, the column and of the transfer matrix should be 1;

[0063] In order to set up conveniently, directly use the intermediate variable stiffness matrix F in the calculation of system matrix to construct, that is, let each column of transfer matrix B be the corresponding column of light source matrix F divided by the sum of the column;

[0064] 3) use coarse grid surface light intensity b1 and coarse grid system matrix A1 and coarse grid adjacent matrix W1, and take coarse grid system matrix A1 as the physical model constraint in the deep network to solve; Specifically:

[0065] 3.1) pass the coarse grid surface light intensity b1 into the deep network, and the output x of the deep network is the light source power density distribution, and the objective function is:

[0066] loss = mse(b1, A1x) + mse(b1, A1B1x) 2 x)

[0067] Adjust the internal parameters of the deep network by the mean square error of the forward fluorescence distribution and the real surface;

[0068] The deep network uses full connection network, the learning rate is 0.0005, the iteration number is 50 times; The network layer number is 3 layers, the length of the input vector is the surface node number, and the length of the output vector is the internal node number;

[0069] 3.2) solve 5 times, record all the results, and take the non-zero nodes after multiplying the coarse grid adjacent matrix W1 as the feasible region;

[0070] 3.3) refer to Figure 3 , continue to solve in the feasible region, and the objective function is:

[0071] loss = mse(b1, A1x) + mse(b1, A1B1x) + mse(b1, A1B1x) + mse(b1, A1B1x) + mse(b1, A1B1x) 2 x) 3 x)

[0072] The iteration of the feasible region is the non-zero node after multiplying the result x of the solution by the coarse grid adjacency matrix W1; the loss value and result of all iterations are recorded, when the program cannot generate a feasible region with the all-zero solution, the result with the minimum loss is selected from the previous results to generate a feasible region, if the loss values of all previous results are greater than the loss value of the all-zero solution, the initial feasible region is used to continue solving; when the number of iterations is reached, the result with the minimum loss in the record is selected as the solution result of this time;

[0073] Referring to Figure 2 , the input of the deep network is the measured tissue surface light intensity, and the output of the deep network is the light source distribution; the forward result corresponding to the given light source distribution can be calculated through the system matrix, and the difference between the forward result and the measured tissue surface light intensity is taken as the loss function to the neural network, and the parameters are learned through iteration;

[0074] 3.4) The residual error of the solution result is solved in the same way as 3.3), and the iteration is stopped when the sum of the solution result is less than 5% of the initial solution result; all solution results are accumulated and the nodes with a value less than 30% of the maximum value are discarded as the final coarse grid solution result;

[0075] 4) Using the result obtained in step 3), the feasible region is divided in the fine grid, and the fine grid surface light intensity b2, the fine grid system matrix A2 and the fine grid adjacency matrix W2 are used to solve; specifically:

[0076] 4.1) According to the solution result of the coarse grid, the fine grid nodes within 1mm range of the non-zero node coordinates are taken as the solution feasible region of the fine grid;

[0077] 4.2) The fine grid surface light intensity b2 is input into the deep network, and the output x of the deep network is the light source power density distribution, and the objective function is:

[0078] loss=mse(b2,A2x)+mse(b2,A2B2x)+mse(b2,A2B2 2 x)+mse(b2,A2B2 3 x)

[0079] The loss value and result of all iterations are recorded, when the program cannot generate a feasible region with the all-zero solution, the result with the minimum loss is selected from the previous results to generate a feasible region, if the loss values of all previous results are greater than the loss value of the all-zero solution, the initial feasible region is used to continue solving; when the number of iterations is reached, the result with the minimum loss in the record is selected as the solution result of this time;

[0080] The deep network of the embodiment uses a full connection network, the learning rate is 0.001, the iteration number is 100 times; the network layer number is 3 layers, the length of the input vector is the surface node number, and the length of the output vector is the feasible region node number;

[0081] 4.3) The residual of the solution result is solved in the same way as 4.2), and the iteration is stopped when the sum of the solution result is less than 5% of the initial solution result; all solution results are accumulated and the nodes with a value less than 30% of the maximum value are discarded as the final fine grid solution result.

[0082] Referring to Figure 4 , Figure 4 The results obtained by the method of the embodiment in mouse tissue have an average error of 0.22 mm.

Claims

1. A method of fluorescence tomographic reconstruction imaging based on a physical model deep network, characterized in that, The method comprises the following steps: 1) Obtain the light intensity data of the tissue surface by CCD camera, and divide the phantom into a coarse grid and a fine grid by Amira software; use the divided grid and the optical parameters of each region in the phantom to solve the coarse grid system matrix , the fine grid system matrix , the coarse grid adjacency matrix of each node of the grid , and the fine grid adjacency matrix of each node of the grid by finite element method; The surface light intensity data is matched to two grids to form a coarse grid surface light intensity corresponding to the coarse grid and a fine grid surface light intensity corresponding to the fine grid ;​ 2) Set the coarse grid transfer matrix for power density , the fine grid transfer matrix The transfer matrix transfers the power density to its neighboring node locations, the column and row sum of the transfer matrix should be 1; 3) Using coarse mesh surface light intensity and coarse mesh system matrix and coarse mesh adjacency matrix and coarse mesh transfer matrix , coarse mesh system matrix solved in deep network as physical model constraints; 4) Using the results obtained in step 3), the feasible region is divided into fine grids, and the surface light intensity , fine grid system matrix and fine grid adjacency matrix and fine grid transfer matrix are solved.

2. The method of claim 1, wherein, The step 1) is specifically: 1.1) meshing the phantom by using Amira software; 1.2) discretizing the diffusion equation according to the nodes of the mesh, and obtaining a system matrix; The diffusion equation is: wherein, is a diffusion coefficient, is a light source, is a luminous flux generated by the light source, is an absorption coefficient, is a refractive index; The discretization is: wherein, is the stiffness matrix, is the light flux matrix, is the light source matrix; Then: System matrix For The result of setting to zero the corresponding row of an internal node 1.3) Find the adjacency matrix of all nodes in the grid .

3. The method of claim 1, wherein, The step 3) is specifically: 3.1) the raw grid surface light intensity to the incoming deep network, the output of the deep network for the light source power density distribution, the objective function is: Adjusting the internal parameters of the network by the mean square error of the fluorescence distribution forward and the real surface; 3.2) Solve many times, record all the results, and multiply the sum of all results by the non-zero nodes after the coarse grid adjacency matrix as the feasible region; 3.2) Solve many times, record all the results, and multiply the sum of all results by the non-zero nodes after the coarse grid adjacency matrix as the feasible region; 3.3) continuing to solve in the feasible region, and the objective function is: The iteration of feasible region is the result of solving Multiply coarse grid adjacency matrix non-zero nodes after the iteration; record the loss value and result of all iteration results, when the program cannot generate a feasible region for the all-zero solution, select the result with the smallest loss value from the previous results to generate a feasible region, if the loss value of all previous results is greater than the loss value of the all-zero solution, continue to solve using the initial feasible region; When the number of iterations is reached, the result with the minimum loss in the record is selected as the solving result of this time; 3.4) the residual of the solving result is solved in the same way, and the iteration is stopped when the value of the solving result is less than a certain value; all solving results are accumulated and the nodes with a value less than 30% of the maximum value are discarded as the final coarse grid solving result.

4. The method of claim 3, wherein, The input of the deep network is the measured tissue surface light intensity, and the output of the deep network is the light source distribution; the tissue surface light intensity corresponding to the given light source distribution is calculated by the system matrix, and the difference between the calculated tissue surface light intensity and the measured tissue surface light intensity is taken as the loss function and transmitted to the neural network, and the parameters are learned by iteration.

5. The method of claim 1, wherein, The step 4) is specifically: 4.1) according to the solving result of the coarse grid, taking the fine grid nodes within 1mm range of the non-zero node coordinates as the solving feasible region of the fine grid; 4.2) the fine grid surface light intensity the incoming deep network, the output of the deep network for the light source power density distribution, the objective function is: Record the loss value and result of all iteration results, when the program cannot produce a feasible region for the all-zero solution, select the result with the minimum loss in the previous results to produce a feasible region, if the loss value of all previous results is greater than the loss value of the all-zero solution, continue to solve by using the initial feasible region; When the number of iterations is reached, the result with the minimum loss in the record is selected as the solving result of this time; 4.3) the residual of the solving result is solved in the same way, and the iteration is stopped when the value of the solving result is less than a certain value; all solving results are accumulated and the nodes with a value less than 30% of the maximum value are discarded as the final coarse grid solving result.

Images