A multi-target automatic recognition reconstruction method based on probability density subspace

Abstract

A multi-target reconstruction method based on subspace decision optimization, acquires biological tissue surface light distribution information, grid node coordinate matrix, grid internal tetrahedron matrix, grid internal tetrahedron index matrix and system matrix as input; through four parts of reconstruction, subspace initialization, subspace decision, subspace optimization to carry out multi-target reconstruction. The method introduces the subspace idea on the multi-target reconstruction problem for the first time, breaks through the limitation of the traditional feasible region method, significantly improves the quality of multi-target reconstruction, and reduces the cost of parameter optimization. In addition, the present application does not need source number prior information, automatically judges the number of subspaces and gradually converges to the real source number, expands the application prospect of reconstruction in clinic. In addition, the reconstruction algorithm in the method can be replaced arbitrarily, and is not dependent on a specific algorithm, which provides an effective tool for three-dimensional reconstruction, especially multi-target reconstruction.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of medical imaging technology, in particular to a multi-target reconstruction method based on subspace decision optimization. BACKGROUND

[0002] Optical molecular imaging technology uses specific probes in vivo as signal sources, reflecting the molecular level changes in the lesion area of the body, and has the advantages of relatively low cost, strong specificity, high sensitivity, etc., and has the ability to detect early tumors.

[0003] In order to better determine the three-dimensional spatial position of the signal source distribution in the body (lesion), fluorescence molecular tomography and bioluminescence tomography technologies based on tomographic imaging technology have been proposed. This kind of technology has similar principles, which are to construct the process of optical signal transmission in biological tissue by means of mathematical tools, to combine the optical signals collected on the surface, and to use inverse solving algorithm to solve the internal source distribution.

[0004] Three-dimensional reconstruction is a clear inverse problem with strong ill-posedness, mainly because the collected information only accounts for a small part of all energy distribution information, and the model is not accurate in describing biological tissue, and the noise of signal collection.

[0005] In order to better solve this kind of inverse problem, the feasible region strategy is often combined with iterative solving means to reduce the space to be solved and reduce the ill-posedness. However, the current iterative reconstruction simply reduces the entire feasible region space according to certain rules, often destroying the geometric characteristics and spatial structure information of the light source; secondly, the current iterative reconstruction method ignores the multi-source situation and simply uses one feasible region, which requires a lot of cost for parameter optimization to obtain good multi-source reconstruction results; in addition, the number of current multi-source reconstruction sources is often input as a prior condition, which obviously does not conform to the clinical reality. SUMMARY

[0006] In view of the deficiencies in the above-mentioned technologies, the present application proposes a multi-target reconstruction method based on subspace decision optimization, which uses subspace instead of the traditional feasible region concept. The subspace is automatically generated in the region and dynamically changes, and converges to the true source number after a certain number of times. In addition, the subspace will make autonomous decisions during the iteration process, achieving the purpose of gradually reducing the space to be solved according to the spatial feature information. After all the subspace iterations are completed, the iteration results are processed based on the normal distribution and energy curve to obtain the final result.

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

[0008] A multi-target reconstruction method based on subspace decision optimization, comprising the following steps:

[0009] Step 1, obtaining the light distribution information of the biological tissue surface, the grid node coordinate matrix, the grid internal tetrahedron matrix, the grid internal tetrahedron index matrix and the system matrix; comprising the following steps:

[0010] 1.1, obtaining a.raw file containing a light source;

[0011] 1.2, obtaining a.grid.am file of the entire simulation model;

[0012] 1.3, obtaining a.mphtxt file of the light source to be used in the simulation model;

[0013] 1.4, obtaining the forward simulation result of the light distribution information of the biological tissue surface;

[0014] 1.5, obtaining the light distribution information of the biological tissue surface, and generating pre-run data;

[0015] Step 2, initializing the global index of the grid node to obtain the initial index column Sp containing all the grid node indexes; reconstructing the initial region formed by the grid node corresponding to the initial index column Sp to obtain the i-th original result X i , that is, the energy intensity value of the internal source;

[0016] Step 3, combining X i , the grid node coordinate matrix, the grid internal tetrahedron matrix, the grid internal tetrahedron index matrix and the system matrix to obtain the L2 norm and the cosine similarity of the current iteration, adding the two indicators and dividing by 2 to obtain the reconstruction weight value of this time

[0017] Step 4, processing as follows:

[0018] 4.1, according to the X ip belonging to the subspace feasible region index, obtaining the probability mean coordinates (X P , Y P , Z P ) of the current iteration, regarding these samples as three-dimensional random variables, and then assembling the covariance matrix M Cov of this time, obtaining the eigenvalues and eigenvectors of the covariance matrix;

[0019] 4.2, using the eigenvectors of M Cov to determine the deflection angle of the subspace and the deflection angle of the grid node coordinates, and using the eigenvalues of M Cov to determine the edge length of the subspace, and the nodes in the initial region Sp contained in the subspace are the subspace Sp_ROI obtained by updating this iteration;

[0020] 4.3, calculate the distance error between the grid node coordinates in the subspace Sp_ROI and the mean coordinate of the probability, arrange the grid nodes in the subspace Sp_ROI in ascending order according to the distance error from small to large;

[0021] 4.4, define the initial area change coefficient, update the subspace Sp_ROI, and the number of grid nodes in the subspace Sp_ROI is divided by the quotient of the number of grid nodes in the subspace Sp_ROI and β 2 The relationship is updated as the standard;

[0022] 4.5, arrange the grid nodes in the initial area Sp in descending order according to the energy intensity to obtain the subspace Sp_Descend, take the number of grid nodes in the current subspace Sp_Descend as the initial node number, update the value of β and the number of grid nodes in the subspace Sp_Descend;

[0023] 4.6, the new index column Sp is the concatenation of the Sp_ROI index column and the Sp_Descend index column, and the grid nodes corresponding to the new index column Sp form a new index area Sp; when the iteration number reaches 5, output the reconstruction result X of this time i ; otherwise, return to step 4.1;

[0024] Step 5, subspace initialization, including the following steps:

[0025] 5.1, select the Sp index area of X as the initial sample of clustering, and perform normalization operation on the energy value to obtain the vector Nodes_DensityCluster, count the length of Nodes_DensityCluster, and magnify 20 times as the amplification multiple Multiple, magnify all values of Nodes_DensityCluster by Multiple times and take the integer part, to obtain the multiplied Nodes_DensityClusterMul vector, each element of which becomes an integer;

[0026] 5.2, generate an empty matrix specimenSpace with Multiple rows and 3 columns as the clustering sample matrix; the value of each row element of Nodes_DensityClusterMul vector is used as the filling multiple, and the coordinate value of the grid node coordinate matrix corresponding to the Sp index is used as the input, and the node coordinate elements are filled in the first three columns of specimenSpace;

[0027] 5.3, take the first three columns of specimenSpace as input, set the possible number of subspaces k from 1 to 20, run the K-Means algorithm one by one, and record the sum of squared error (Sum Squared Error, SSE) of each clustering to obtain the k-SSE curve;

[0028] 5.4. Dual-line inflection point search: For the obtained k-SSE curve, define two straight lines, and calculate the compensation value D for each k value based on these two lines. offset(k) Let D offset The smallest k value is the elbow of this curve, which is also the number K of subspaces to be initialized this time;

[0029] 5.5. Reinitialize the global index of the grid nodes to obtain an initial index column Sp containing the indices of all grid nodes; reconstruct the initial region Sp formed by the grid nodes corresponding to the initial index column Sp to obtain X. i , that is, the energy intensity value of the internal light source;

[0030] 5.6. Based on the current second feasible region index Sp, the number of light sources K obtained in step 5.4, the grid node coordinate matrix, and the reconstruction result X. i Perform K-means clustering to obtain K subspace feasibility region indices SpA, SpB, etc.;

[0031] Step 6, Subspace Decision: For each subspace, execute steps 3 to 4. After every 5 iterations, return to step 4 and reinitialize the subspace.

[0032] Step 7, Subspace Optimization, includes the following steps:

[0033] 7.1 Weights for all iterations The data is sorted to obtain the frequency distribution of the weights, and then fitted using a normal distribution. The mean and standard deviation of this normal distribution are determined. The iteration counts outside the range of mean ± one standard deviation are then compared with the corresponding P-values. Err Remove, and you get the remaining X′ and P′. Err ;

[0034] 7.2, X′ ip and Weighted summation yields the result S of the number of filtering iterations. F ;

[0035] 7.3 For each subspace in the last iteration, sort the corresponding node energies in descending order to obtain the node-energy curve. Use the method in step 6.4 to find the elbow point of each node-energy curve, remove the nodes after the elbow point, and finally obtain the final result S. opt .

[0036] The beneficial effects of this invention are:

[0037] The present application firstly obtains clustering samples according to iterative reconstruction and sample resampling, so as to convert the input of clustering into the input of points, and enhance the weight of high-energy results. The number of subspaces is obtained by clustering, which is an automatic process. After the number of subspaces is determined, the subspaces are initialized and iteratively determined, so as to realize the fitting of the space to the source and reduce the space to be solved. After the reconstruction is completed, the normal distribution curve is constructed based on the weight of each reconstruction, and a part of the iteration number is filtered according to the weight. After the iteration number filtering is completed, the node energy is arranged in descending order, and the elbow point is found. The two steps realize the optimization of the subspace.

[0038] The main research results of the present application are as follows: (1) an unsupervised clustering method based on double straight line search combined with sample resampling is proposed, the number of subspaces is automatically generated according to the reconstruction result, and initialization is performed, and the number of subspaces will converge to the true source number in the subsequent iteration, without the need for manual input, and more in line with clinical practice; (2) the concept of subspace is introduced to jump out of the limitation of the traditional feasible region, and the space change decision is made independently for each subspace, avoiding the phenomenon that the traditional feasible region method based on a single feasible region is not good for fitting the multi-source distribution, and the number of subspaces will dynamically change to ensure the reconstruction quality of multiple targets; (3) weights are assigned to the reconstruction results of all subspaces, the weights are filtered according to the characteristics of normal distribution, and the node intensity is filtered after the filtering is completed, so as to improve the quality of the final result. The method does not need to know the number of targets in advance when reconstructing multiple targets, can dynamically divide subspaces, and gradually fit the number of real sources in the iteration process, and can better fit multiple targets. In addition, the internal reconstruction algorithm of the method can be arbitrarily replaced, and is not dependent on specific algorithms or parameters, which provides an effective tool for three-dimensional reconstruction, especially multiple target reconstruction. BRIEF DESCRIPTION OF DRAWINGS

[0039] Figure 1 The coordinate diagram of the heart, lung, liver, stomach, kidney and muscle.

[0040] Figure 2 The.grid.am file diagram of the entire simulation model;

[0041] Figure 3 The light source grid diagram;

[0042] Figure 4 The diagram of filling node coordinate elements in the first three columns of specimenSpace;

[0043] Figure 5 The light source center coordinate diagram of the test sample;

[0044] Figure 6 The experimental results of the simulation experiment;

[0045] Figure 7For simulation experiment results

[0046] Figure 8 For subspace change graph

[0047] Figure 9 For frequency filtering graph

[0048] Figure 10 For energy value filtering graph.

[0049] Specific implementation examples

[0050] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only some of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative work fall within the scope of protection of the present application.

[0051] The method of the present application is based on bioluminescence tomography technology and a physical model of light propagation in biological tissue. The intensity distribution information of the photons is obtained by collecting the photons emitted from the surface of the biological tissue. Then, the spatial position and signal intensity distribution of the internal light source target are reversely calculated by combining the anatomical structure inside the biological tissue and the optical parameters (such as absorption coefficient, scattering coefficient, etc.) of the corresponding tissues, so as to quantitatively display the early lesion change in the biological tissue. Specifically, the present application relates to a multi-target reconstruction method based on subspace decision optimization in the field of bioluminescence tomography technology. The method can be used to quickly find the accurate position of the lesion area, has the advantages of not needing the number of light sources to be input, not being dependent on specific parameters and specific algorithms, etc., and can improve the performance of the current multi-target reconstruction.

[0052] The following is the explanation of the terms appearing in the present application:

[0053] The light distribution information on the surface of the biological tissue refers to the energy distribution vector on the surface of the biological tissue.

[0054] The grid refers to a plurality of tetrahedral structures divided in the model reconstruction process. The four vertices of the tetrahedral structure are grid nodes. Adjacent nodes are connected by the edges of the tetrahedral structure. In an analysis model, a plurality of grids and a plurality of grid nodes are present.

[0055] Subspace refers to an area divided in the entire space based on certain rules.

[0056] The grid node coordinate matrix refers to the coordinates of the four vertices (nodes) of each tetrahedron, which are the coordinates of the grid nodes. The XYZ coordinates can form the grid node coordinate matrix.

[0057] Grid internal tetrahedron matrix: refers to each vertex coordinate is numbered, and four vertexes of each tetrahedron can be obtained, which is a tetrahedron matrix.

[0058] Grid internal tetrahedron index: refers to in three-dimensional reconstruction, a tetrahedron grid is used as a basic unit, the tetrahedron grid constitutes a spatial region of the entire biological tissue, and the tetrahedron is numbered to obtain a tetrahedron index,

[0059] System matrix: refers to a parameter matrix composed of mathematical relations and logical mathematical factors describing the propagation of photons in the tissue.

[0060] Global index or index: the global index refers to that all grid nodes are used for reconstruction in the initial first reconstruction, and the numbering set of the grid nodes is the global index, and each grid node number is an index;

[0061] Index region: the region composed of the grid nodes corresponding to the index is referred to as an index region.

[0062] Embodiment one: the reconstruction method of the application comprises the following steps:

[0063] Based on the digital mouse simulation model of the University of Southern California, USA, a simulation experiment is carried out. The digital mouse simulation model is obtained by extracting tissue information from imaging slice data of a real experimental animal through CT technology, and then through tomographic reconstruction means. In order to reduce the calculation complexity and save system resources, the head and tail regions of the mouse are removed, and only the trunk part is reserved, and the model size is 3.8 cm long, 2.08 cm wide and 3.5 cm high. According to the internal tissue structure, it can be divided into six organs and the remaining tissue of the body, from top to bottom, which are heart, lung, liver, stomach, kidney and muscle. As shown in Figure 1 .

[0064] The operation process of constructing a forward energy simulation data is described in detail as follows. The steps are as follows:

[0065] Step one: data acquisition;

[0066] I, obtain a.raw file containing a light source.

[0067] The digital mouse simulation model of the University of Southern California is a.raw file with a size of 3.8 cm long, 2.08 cm wide and 3.5 cm high, and does not contain a light source inside.

[0068] II, obtain the.grid.am file of the entire simulation model.

[0069] Amira software is a tool for meshing, which is used to mesh the.raw file obtained in I to obtain the mesh structure of the whole simulation model, and the file format is.grid.am file. The example of the whole simulation model.grid.am file obtained in this step is shown in Figure 2 .

[0070] III. Obtain the light source.mphtxt file to be used in the simulation model.

[0071] COMSOL Multiphysics software is a tool for generating standard geometric and physical meshes, which can set specific geometric shapes as needed, and the file format is.mphtxt file. The light source mesh file obtained in this step is shown in Figure 3 .

[0072] IV. Obtain the forward simulation results of the surface light distribution information of the biological tissue.

[0073] The specific luminescent probe in the biological tissue emits photons spontaneously or under excitation, and the photons undergo reflection, absorption, scattering and other physical and optical processes during the transmission from inside to outside, and form light distribution information on the surface of the biological body. The surface light distribution information of the biological tissue has a one-to-one nonlinear mapping relationship with the light source information. The.grid.am file of the whole simulation model obtained in II is used to obtain the surface light distribution information of the biological tissue by forward simulation. Mose software is used to read the mesh file obtained in II, which can convert the.grid.am file into an.off file according to the internal organization information, and each organ corresponds to an.off file. Then the Mose software can read the.off file generated by itself, simulate the photon transmission based on the Monte Carlo method, and thus perform forward energy simulation. The file obtained by forward energy simulation is in the.CW format.

[0074] V. Obtain the surface light distribution information of the biological tissue and generate the pre-data.

[0075] The surface light distribution information of the biological tissue is extracted from the forward simulation result file in the.CW format obtained in IV using the Matlab programming language, and the.grid.am mesh file generated in II is read in accordingly. According to the mesh file and the forward energy file, the surface energy distribution vector, the grid node coordinate (X, Y and Z) matrix, the grid internal tetrahedron matrix, the grid internal tetrahedron index matrix and the system matrix are generated, and these results are saved as the corresponding.mat files of Matlab, including B_Energy.mat, Nodes.mat, InsideElement.mat, InsideElementIndex.mat, G.mat and triMat.mat, a total of 6 mat files.

[0076] Step 2: Reconstruction; The .mat file obtained in Step 1 is used as input for reconstruction. The reconstruction process mainly includes the following steps:

[0077] A. Initialize the global index by generating the initial index column Sp for all nodes in the current initial region based on the Nodes.mat file;

[0078] B. For the initial index column Sp, the reconstruction algorithm proposed in He's paper "Sparse reconstruction for quantitative bioluminescence tomography based on the incomplete variables truncated conjugate gradient method" is used to reconstruct the initial region Sp formed by the grid nodes corresponding to the current initial index column Sp, thus obtaining X. i The energy intensity value of the internal light source;

[0079] C. Based on the currently obtained X i By combining the imported G.mat file and the forward simulation energy.CW file, the current L2 error rate E can be calculated. L2 Cosine similarity E cos Based on these two values, combined with the previous X i The current global probability weight value P can be calculated. i =0.5*(P L2 +P cos ),in Further obtain X i Corresponding X iP =X i *P i ; i represents the i-th iteration currently in progress, and the value of i is a positive integer;

[0080] The global probability weight value X in the current iteration ip ;

[0081] D. According to X ip Obtain the coordinates of the probability mean of the current iteration (X). P Y P Z P );

[0082] E. Treat the three-dimensional coordinates as a set of three-dimensional random variables, and assemble the covariance matrix M for this iteration. Cov Further, we obtained M Cov The eigenvalues ​​Val and eigenvectors Vec;

[0083] F. According to the Val and Vec obtained in E, the deflection angle of the subspace and the deflection angle of the grid node coordinates are determined by Vec, and the edge length in each direction of the subspace interval is determined by Val, eight vertices of the subspace region can be determined, the region surrounded by the eight vertices is divided to obtain the subspace region of the current time, and the nodes in the initial region Sp contained in the subspace are the subspace Sp_ROI obtained by the iteration update this time;

[0084] G. Record the number of grid nodes in the current subspace Sp_ROI, save it as Cut_Num, and use Mohamed A. Naser's Improved bioluminescence and fluorescence reconstruction algorithms using diffuse optical tomography, normalized data, and optimized selection of the permissible source region article to calculate the number of grid nodes in the current subspace Sp_ROI according to the formula

[0085] Calculate the current attenuation factor β, calculate the distance error between the coordinates in the subspace Sp_ROI and (X P , Y P , Z P ) in ascending order, and arrange the grid nodes in the subspace Sp_ROI in ascending order;

[0086] H. According to the relationship that Cut_Num is divided by β 2 , if it is greater than 2, the half edge length of the next subspace space is twice the length of Val; if it is less than or equal to 2 but greater than 1, the half edge length of the next subspace space is 1 times the length of Val; if it is less than or equal to 1, the half edge length of the next subspace space is 0.5 times the length of Val;

[0087] I. According to the energy intensity of the grid nodes in the initial region Sp in descending order, the subspace Sp_Descend is obtained, the current number of grid nodes is taken as the initial number of nodes, and the β value and Cut_Num are updated according to the formula

[0088] J. The new index column Sp is the index column Sp_ROI spliced with the index column Sp_Descend;

[0089] K. The new index column Sp is the Sp_ROI index column spliced with the Sp_Descend index column, and the grid nodes corresponding to the new index column Sp form a new index region Sp; when the iteration number reaches 5, the reconstruction result X of this time is output; otherwise, return to step B;

[0090] ​Step three, initialize the subspace;

[0091] L. Select the Sp index region of X as the initial sample of clustering, and normalize the energy value to obtain the vector Nodes_DensityCluster. Count the length of Nodes_DensityCluster and amplify it by 20 times as the amplification Multiple. Amplify all values of Nodes_DensityCluster by Multiple times and take the integer part to obtain the multiplied Nodes_DensityClusterMul vector. Each element of this vector becomes an integer.

[0092] M. Generate a Multiple row, 3 column empty matrix specimenSpace as the clustering sample matrix. The value of each row element of Nodes_DensityClusterMul is used as the filling multiple, and the coordinate value of the Sp index corresponding grid node coordinate matrix is used as the input to fill the node coordinate elements in the first three columns of specimenSpace. As shown in the example, the 21471th node corresponds to the coordinates (17.036, 8.751, 15.910) and needs to be amplified by 21 times. The coordinates (17.036, 8.751, 15.910) are printed 21 times in the specimenSpace matrix. Figure 4

[0093] N. Take the first three columns of specimenSpace as input, set the possible number of subspace k from 1 to 20, run the K-Means algorithm one by one, and record the Sum Squared Error (SSE) of each clustering to obtain the k-SSE curve, which is defined as follows Two straight lines:

[0094]

[0095]

[0096] According to the two straight lines as above, the compensation value D of each k value is calculated by the following formula offset(k) :

[0097]

[0098] Let D offset The minimum k value is the elbow point of this curve, that is, the number of subspace K to be initialized this time.

[0099] O. Initialize the global index of the grid nodes again to obtain the initial index column Sp containing all the grid node indexes. Reconstruct the initial region Sp formed by the grid nodes corresponding to the initial index column Sp to obtain X​i i.e. the energy intensity value of the internal light source;

[0100] P. According to the current sub-region index Sp, the number of light sources K obtained in step N, the grid node coordinate matrix, the reconstruction result X i K-means clustering is performed to obtain K sub-space feasibility region indexes SpA, SpB, etc.

[0101] Step four, sub-space decision;

[0102] Q. For each sub-space, perform the previous step two, dynamically decide the internal space change of each sub-space;

[0103] Step five, sub-space optimization;

[0104] R. The weight of all iteration numbers Sort to obtain the weight distribution frequency, and use the normal distribution to fit the distribution to determine the mean and standard deviation of the normal distribution, remove the iteration number results and the corresponding P outside the range of mean ± one standard deviation Err to obtain the remaining X' and P' Err ;

[0105] S. Weight X' ip and to obtain the filtered iteration number result S F ;

[0106] T. For each sub-space of the last iteration, arrange the corresponding node energy in descending order to obtain the node-energy curve, use the method of step N to find the elbow point of each node-energy curve, remove the nodes after the elbow point, and finally obtain the final result S opt

[0107] Example two:

[0108] In order to verify the good performance of the proposed framework, the present application sets up a double light source simulation experiment. The present application uses the digital mouse simulation model of the University of Southern California, and the model size is 3.8 cm long, 2.08 cm wide, and 3.5 cm high. In order to reduce the calculation complexity and save system resources, the head and tail regions of the mouse are removed, and only the torso part is reserved. According to its internal tissue structure, it can be divided into six organs and the remaining tissues of the body, from top to bottom, which are heart, lung, liver, stomach, kidney and muscle. The light source center coordinates of the test sample in this experiment are set at (15, 8, 15) mm and (22, 8, 15) mm as shown in Figure 5As shown, the light source radius is set to 1mm. Then, the surface light distribution information of the biological tissue is obtained by forward simulation. The experimental results show that the framework of the present application has good target positioning and morphological reconstruction performance. The reconstructed light source center coordinates of the test sample are (15.35, 8.15, 14.95) mm and (22.27, 8.11, 14.70) mm, the positioning errors are 0.39mm and 0.50mm respectively, and the DICE coefficients are 64.7% and 73.1% respectively. The experimental results of the simulation experiment are shown in Figure 6 and Figure 7 As shown. The subspace change is shown in Figure 8 The frequency filtering and energy value filtering are shown in Figure 9 and Figure 10

[0109] From the above implementation cases, it can be seen that the method proposed by the present application has good target positioning and light source geometry recovery performance, and does not need to input the light source artificially.

[0110] For those skilled in the art, other various corresponding changes and modifications can be made to the above-described technical solutions and concepts, and all of these changes and modifications should belong to the protection scope of the claims of the present application.

[0111] The preferred embodiments of the present disclosure are described in detail above in combination with the drawings, but the present disclosure is not limited to the specific details in the above-described embodiments. Within the technical concept range of the present disclosure, various simple modifications can be made to the technical solutions of the present disclosure, and these simple modifications all belong to the protection range of the present disclosure.

[0112] In addition, it should be noted that each specific technical feature described in the above specific embodiments can be combined in any appropriate manner without contradiction. In order to avoid unnecessary repetition, the present disclosure will not further describe various possible combination manners.

[0113] In addition, any combination of various different embodiments of the present disclosure can also be made, as long as it does not deviate from the idea of the present disclosure, and it should also be considered as the disclosed content of the present disclosure.​

Claims

1. A multi-objective reconstruction method based on subspace decision optimization, characterized by comprising the following steps: Step 1: Obtain the surface light distribution information of biological tissue, the grid node coordinate matrix, the internal tetrahedral matrix of the grid, the internal tetrahedral index matrix of the grid, and the system matrix; including the following steps: 1.1 Obtain the .raw file containing the light source; 1.2 Obtain the .grid.am file of the entire simulation model; 1.3 Obtain the .mphtxt file containing the light source to be used in the simulation model; 1.

4. Forward simulation results for obtaining light distribution information on the surface of biological tissues; 1.5 Obtain light distribution information on the surface of biological tissues and generate pre-run data; Step 2: Initialize the global index of grid nodes to obtain the initial index column Sp containing the indexes of all grid nodes; The initial region formed by the grid nodes corresponding to the initial index column Sp is reconstructed to obtain the i-th original result X. i , that is, the energy intensity value of the internal light source; Step 3, combine X i The L2 norm and cosine similarity of the current iteration are obtained from the grid node coordinate matrix, the grid interior tetrahedral matrix, the grid interior tetrahedral index matrix, and the system matrix. These two indices are added together and divided by 2 to obtain the reconstruction weight value. Step 4: Follow these steps to process: 4.

1. Based on X, which belongs to the feasible region index of the subspace ip The row yields the coordinates of the probability mean (X) of the current iteration. P Y P Z P Furthermore, by treating the three coordinates as a set of three-dimensional random variables, the covariance matrix M can be assembled for this calculation. Cov This yields the eigenvalues ​​and eigenvectors of the covariance matrix. 4.2 Utilizing M Cov The eigenvectors determine the deflection angle of the subspace and the deflection angle of the grid node coordinates, using M... Cov The eigenvalues ​​determine the side length of the subspace, and the nodes in the initial region Sp contained in the subspace are the subspace Sp_ROI obtained in this iteration update; 4.3 Calculate the distance error between the coordinates of the grid nodes in the subspace Sp_ROI and the coordinates of the probability mean, and sort the grid nodes in the subspace Sp_ROI in ascending order of distance error from smallest to largest; 4.4 Define the initial region variation coefficient, update the subspace Sp_ROI by dividing the number of grid nodes in the subspace Sp_ROI by β. 2 The relationship between the quotients is used as the updated standard; 4.

5. Sort the grid nodes in the initial region Sp in descending order of energy intensity to obtain the subspace Sp_Descend. Use the number of grid nodes in the current subspace Sp_Descend as the initial number of nodes, and update the β value and the number of grid nodes in the subspace Sp_Descend. 4.6 The new index column Sp is formed by concatenating the Sp_ROI index column with the Sp_Descend index column. The grid nodes corresponding to the new index column Sp form the new index region Sp. When the iteration count reaches 5, the reconstruction result X of this iteration is output. i Otherwise, return to step 4.1; Step 5, Subspace Initialization, includes the following steps: 5.1 Select the Sp index region of X as the initial sample for clustering, and normalize the energy value to obtain the vector Nodes_DensityCluster. Count the length of Nodes_DensityCluster and magnify it by 20 times as the amplification factor Multiple. Magnify all values ​​of Nodes_DensityCluster by Multiple times and round down to obtain the multiplied Nodes_DensityClusterMul vector, where each element of the vector becomes an integer. 5.2 Generate a multiple-row, 3-column empty matrix specimenSpace as the clustering sample matrix; the value of each row element of the Nodes_DensityClusterMul vector is used as the padding factor, and the coordinate values ​​of the grid node coordinate matrix corresponding to the Sp index are used as input to fill the first three columns of specimenSpace with node coordinate elements; 5.

3. Using the first three columns of specimenSpace as input, set the number of possible subspaces k from 1 to 20, run the K-Means algorithm one by one, and record the sum of squared errors (SSE) for each clustering to obtain the k-SSE curve; 5.

4. Dual-line inflection point search: For the obtained k-SSE curve, define two straight lines, and calculate the compensation value D for each k value based on these two lines. offset(k) Let D offset The smallest value of k is the number of subspaces K to be initialized this time; 5.

5. Reinitialize the global index of the grid nodes to obtain an initial index column Sp containing the indices of all grid nodes; reconstruct the initial region Sp formed by the grid nodes corresponding to the initial index column Sp to obtain X. i , that is, the energy intensity value of the internal light source; 5.

6. Based on the current second feasible region index Sp, the number of light sources K obtained in step 5.4, the grid node coordinate matrix, and the reconstruction result X. i Perform K-means clustering to obtain K subspace feasibility region indices SpA, SpB, etc.; Step 6, Subspace Decision: For each subspace, execute steps 3 to 4. After every 5 iterations, return to step 4 and reinitialize the subspace. Step 7, Subspace Optimization, includes the following steps: 7.1 Weights for all iterations The data is sorted to obtain the frequency distribution of the weights, and then fitted using a normal distribution. The mean and standard deviation of this normal distribution are determined. The iteration counts outside the range of mean ± one standard deviation are then compared with the corresponding P-values. Err Remove, and you get the remaining X' and P'. Err ; 7.2, X' ip and Weighted summation yields the result S of the number of filtering iterations. F ; 7.3 For each subspace in the last iteration, sort the corresponding node energies in descending order to obtain the node-energy curves. Use the method in step 6.4 to find the elbow point of each node-energy curve, remove the nodes after the elbow point, and finally obtain the final result S. opt .

2. The multi-objective reconstruction method based on subspace decision optimization according to claim 1, characterized in that, The acquisition of the surface light distribution information of biological tissue, grid node coordinate matrix, grid interior tetrahedron matrix, grid interior tetrahedron index matrix, and system matrix in step 1 includes the following steps: I. Obtain the .raw file containing the light source; II. Obtain the .grid.am file of the entire simulation model; The .raw file obtained in I is meshed using Amira software to obtain the mesh structure of the entire simulation model. The file format is .grid.am. III. Obtain the .mphtxt file of the light source to be used in the simulation model; As needed, set specific geometric shapes and use COMSOL Multiphysics software to generate a light source mesh file in .mphtxt format. IV. Forward simulation results for obtaining light distribution information on the surface of biological tissues; There is a one-to-one nonlinear mapping relationship between the light distribution information on the biological surface and the light source information. The light distribution information on the surface of biological tissue is obtained from the .grid.am file of the entire simulation model obtained in II through forward simulation. The grid file obtained in II is read into the Mose software, and the .grid.am file is converted into an .off file according to its internal tissue information. Then, the Mose software can read the .off file it generates and perform photon transmission simulation based on the Monte Carlo method, thereby performing forward energy simulation. The file obtained from the forward simulation is in .CW format. V. Obtain light distribution information on the surface of biological tissues and generate pre-data; The Matlab programming language is used to extract the surface light distribution information of biological tissue from the forward simulation result file in .CW format obtained from IV, and the corresponding .grid.am grid file generated in II is read in. Based on the grid file and the forward simulation file, the surface energy distribution vector, grid node coordinate (X, Y, and Z) matrix, grid interior tetrahedral matrix, grid interior tetrahedral index matrix, and system matrix are generated. These results are saved as corresponding .mat files in Matlab, including a total of 6 .mat files: B_Energy.mat, Nodes.mat, InsideElement.mat, InsideElementIndex.mat, G.mat, and triMat.mat.

3. The multi-objective reconstruction method based on subspace decision optimization according to claim 2, characterized in that, The .mat file obtained in step 1 is used as input for subspace reconstruction, which includes the following steps: I. Initialize the global index: Based on the Nodes.mat file, generate the initial index column Sp for all nodes in the current initial region; II. For the initial index column Sp, reconstruct the initial region Sp formed by the grid nodes corresponding to the current initial index column Sp to obtain X. i The energy intensity value of the internal light source; III. Based on the currently obtained X i By combining the input G.mat file and the forward simulation .CW file, the L2 error rate E for the current iteration is calculated. L2 Cosine similarity E cos Based on these two values, combined with the previous X i The current global probability weight value P can be calculated. i =0.5*(P L2 +P cos ),in Further obtain X i Corresponding X iP =X i *P i ; i represents the i-th iteration currently in progress, and the value of i is a positive integer; Calculate the global probability weight value X for the current iteration. ip ; IV. According to X ip Obtain the coordinates of the probability mean of the current iteration (X). P Y P Z P ); V. Treat the three-dimensional coordinates as a set of three-dimensional random variables, and assemble the covariance matrix M for this step. Cov Further, we obtained M Cov The eigenvalues ​​Val and eigenvectors Vec; VI. Based on Val and Vec, use Vec to determine the deflection angle of the subspace and the deflection angle of the grid node coordinates, and use Val to determine the side length of each direction of the subspace interval. This allows us to determine the eight vertices of the subspace region, divide the region enclosed by the eight vertices, and obtain the current subspace region. The nodes in the initial region Sp contained in the subspace are the subspace Sp_ROI obtained in this iteration update. VII. Record the number of grid nodes in the current subspace Sp_ROI and save it as Cut_Num. Then, use the formula β = (number of grid nodes in Sp_ROI / number of grid nodes in Sp_ROI) to calculate the number of nodes. Calculate the current attenuation factor β, and calculate the coordinates in the subspace Sp_ROI and (X) P Y P Z P The distance error between the nodes in the Sp_ROI subspace is sorted in ascending order. VIII. Based on Cut_Num divided by β 2 The relationship is as follows: if the value is greater than 2, the half-side length of the next subspace is twice the Val length; if the value is less than or equal to 2 but greater than 1, the half-side length of the next subspace is once the Val length; if the value is less than or equal to 1, the half-side length of the next subspace is 0.5 times the Val length. IX. The subspace Sp_Descend is obtained by sorting the grid nodes in the initial region Sp in descending order of energy intensity. The current number of grid nodes is used as the initial number of nodes, according to the formula... Update the β value and the Cut_Num value; X, The new index column Sp is formed by concatenating the index column Sp_ROI with the index column Sp_Descend; XI. The new index column Sp is formed by concatenating the Sp_ROI index column with the Sp_Descend index column. The grid nodes corresponding to the new index column Sp form the new index region Sp. When the number of iterations reaches 5, output the reconstruction result X of this step; otherwise, return to step II.

4. The multi-objective reconstruction method based on subspace decision optimization according to claim 1, characterized in that, The initialization of the subspace includes the following steps: I. Select the Sp index region of X as the initial sample for clustering, and normalize the energy value to obtain the vector Nodes_DensityCluster. Calculate the length of Nodes_DensityCluster and multiply it by 20 as the amplification factor Multiple. Multiply all values ​​of Nodes_DensityCluster by Multiple and round down to obtain the multiplied Nodes_DensityClusterMul vector, where each element of the vector becomes an integer. II. Generate a multiple-row, 3-column empty matrix specimenSpace as the clustering sample matrix; the value of each row element of the Nodes_DensityClusterMul vector is used as the padding factor, and the coordinate values ​​of the grid node coordinate matrix corresponding to the Sp index are used as input to fill the first three columns of specimenSpace with node coordinate elements. III. Using the first three columns of specimenSpace as input, set the number of possible subspaces k from 1 to 20, run the K-Means algorithm one by one, and record the sum of squared errors (SSE) for each clustering to obtain the k-SSE curve. Define the following two lines: Based on the two lines above, the compensation value D for each k value is calculated using the following formula. offset(k) : Let D offset The smallest k value is the elbow of this curve, which is also the number K of subspaces to be initialized this time; IV. Reinitialize the global index of the grid nodes to obtain an initial index column Sp containing the indices of all grid nodes; reconstruct the initial region Sp formed by the grid nodes corresponding to the initial index column Sp to obtain X. i , that is, the energy intensity value of the internal light source; V. Based on the current second feasible region index Sp, the number of light sources K obtained in step N, the grid node coordinate matrix, and the reconstruction result X. i Perform K-means clustering to obtain K subspace feasibility region indices SpA and SpB.

5. The multi-objective reconstruction method based on subspace decision optimization according to claim 1, characterized in that, The subspace optimization includes the following steps: I. Weights for all iterations The data is sorted to obtain the frequency distribution of the weights, and then fitted using a normal distribution. The mean and standard deviation of this normal distribution are determined. The iteration counts outside the range of mean ± one standard deviation are then compared with the corresponding P-values. Err Remove, and you get the remaining X' and P'. Err ; II. X' ip and Weighted summation yields the result S of the number of filtering iterations. F ; III. For each subspace in the last iteration, sort the corresponding node energies in descending order to obtain the node-energy curves. Use the method in step 6.4 to find the elbows of each node-energy curve, remove the nodes after the elbows, and finally obtain the final result S. opt .

Images