CN105326475B - A kind of Bioluminescence tomography reconstruction method differentiated based on multiple light courcess - Google Patents
A kind of Bioluminescence tomography reconstruction method differentiated based on multiple light courcess Download PDFInfo
- Publication number
- CN105326475B CN105326475B CN201510589426.XA CN201510589426A CN105326475B CN 105326475 B CN105326475 B CN 105326475B CN 201510589426 A CN201510589426 A CN 201510589426A CN 105326475 B CN105326475 B CN 105326475B
- Authority
- CN
- China
- Prior art keywords
- matrix
- clustering
- node
- imaging target
- light source
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Abstract
The invention discloses a kind of Bioluminescence tomography reconstruction method differentiated based on multiple light courcess, priori is used as by the use of optical property parameter and anatomical information, with reference to multispectral fluorescence data, the accuracy of BLT light source reconstructions and the resolution ratio of reconstruction image are improved；Using adaptive area of feasible solutions iterative shrinkage strategy, reduce the pathosis of BLT Problems of Reconstruction, reduce dependence of the reconstructed results for algorithm for reconstructing；A kind of mixing clustering method of combination AP clusters and K means clusters is proposed, is avoided individually using the limitation and deficiency of respective clustering algorithm.This method has not only directly used the energy distribution information of grid node, and the information such as the affiliated tetrahedron energy of reconstruction node, centre coordinate, volume are taken into full account, the stability of clustering algorithm is improved to the full extent, is reduced algorithm and is calculated time complexity.
Description
Technical Field
The invention belongs to the field of molecular imaging, and relates to a bioluminescence tomography reconstruction method based on multi-light-source resolution.
Background
Bioluminescence tomography (hereinafter referred to as BLT) is an important imaging modality in optical molecular imaging technology, and utilizes the biochemical reaction of luciferase and a luciferin substrate in active molecules or cells for imaging; the change of the biomolecular cell level is reflected qualitatively and quantitatively through the accurate positioning of the luminous light source.
BLT reconstruction is an ill-defined problem whose solution is not unique and is susceptible to measurement noise. In the current research, the accuracy of the reconstruction result is improved mainly by setting strategies such as a light source feasible region, a self-adaptive finite element, a self-adaptive grid, an effective regularization method and the like. The Von Jinchao provides a feasible region automatic determination method according to the reconstructed light source energy density in the BLT research, so that the instability of manually calibrating a feasible region is avoided; the region selection of the automatic feasible region depends on the reconstructed light source to a great extent, and the accuracy of the light source reconstruction in each time cannot be ensured generally, so that the automatic feasible region reconstruction method is not universal; lvyujie performs self-adaptive subdivision on finite element grids in BLT research according to a parameter estimation method, so that the quantification and positioning accuracy are improved, and the robustness of a reconstruction algorithm is improved; although the self-adaptive finite element and the self-adaptive mesh have high stability, the subdivision of the mesh needs to recalculate a system matrix, the mathematical computation is complex, and the computation cost is high. Algorithms commonly used in the study of BLT reconstruction algorithms include: based on L 2 Norm reconstruction algorithm (Tikhonov regularization); a total variation reconstruction algorithm; based on L 1 Norm reconstruction algorithms (IVTCG reconstruction algorithms, homtopy algorithms, etc.), greedy thought-based reconstruction algorithms; based on L P (0<P&And lt 1) a reconstruction algorithm of the norm. However, these mathematical optimization algorithms are faced with the determination of parameters, and improper parameters will seriously affect the reconstruction results of the algorithms. In addition to single light source reconstruction, the calibration, centering, and size analysis of multiple lesion regions in BLT preclinical studies are of greater importance for guiding clinical treatment. However, reconstruction difficulty of multi-light source reconstruction is higher than that of single light source reconstruction, and a robust and efficient reconstruction method is not available in BLT research for multi-target resolution at present.Therefore, multi-source reconstruction has been a major aspect of BLT research.
Reference documents:
【1】Virostko J,Powers A C,Jansen E D.Validation of luminescent source reconstruction using single-view spectrally resolved bioluminescence images[J].Applied optics,2007,46(13):2540-2547.
【2】Frey B J,Dueck D.Clustering by passing messages between data points[J].science,2007,315(5814):972-976.
disclosure of Invention
In view of the above problems and drawbacks of the prior art, an object of the present invention is to provide a bioluminescence tomography reconstruction method based on multi-light source resolution.
In order to achieve the purpose, the invention adopts the following technical scheme:
a bioluminescence tomography reconstruction method based on multi-light source resolution specifically comprises the following steps:
step 1: obtaining anatomical structure information, optical specificity parameters and surface multispectral fluorescence data of an imaging target;
step 2: normalizing the surface multispectral fluorescence data obtained in the step (1) to obtain normalized surface fluorescence data;
and step 3: constructing a linear relation between surface fluorescence data and the distribution of an internal light source of an imaging target based on a light transmission model and a finite element method;
and 4, step 4: converting the linear relation obtained in the step 3 into an optimization problem;
and 5: performing algebraic iterative solution on the optimization problem in the step 4 by adopting a self-adaptive feasible region iterative contraction strategy to obtain an optimal solution X * And completing the reconstruction.
Further, the multi-light-source resolution-based bioluminescence tomography reconstruction method further comprises the following steps:
and 6: for the optimal solution X obtained in the step 5 * Carrying out AP center clustering, and calculating to obtain a clustering center m of each type;
and 7: taking the clustering center obtained in the step 6 as an initial clustering center of a k-means clustering algorithm, carrying out k-means clustering, and calculating to obtain a clustering center C of each type;
and step 8: the optimal solution X obtained in the step 5 is * Performing geometric 3D display to realize the display of the reconstruction result of the light source inside the imaging target; and (4) performing geometric 3D display on the clustering center C obtained in the step (7) to determine the central position of the light source in the imaging target.
Specifically, the specific process of step 3 includes:
step 3.1: obtaining an optical transmission model according to the anatomical structure information and the optical specificity parameters obtained in the step 1, wherein the optical transmission model is a diffusion approximation equation, and the transmission process of light in the imaging target is described by using the diffusion approximation equation;
step 3.2: carrying out grid dispersion on the imaging target by using a nirfast tool to obtain discrete grid data of the imaging target, wherein the discrete grid data comprises grid nodes and tetrahedrons, and the corresponding relation between the grid nodes and the tetrahedrons is obtained;
step 3.3: and (3) dispersing the diffusion approximation equation in the step (3.1) according to the discrete grid data and a finite element theory, and further constructing a linear relation between the surface fluorescence data and the distribution of the internal light source of the imaging target:
Φ(λ)＝M(λ)X
wherein phi (lambda) is surface fluorescence data with the wavelength of lambda, M (lambda) represents a system matrix of the fluorescence data under the wavelength of lambda, and X is the energy distribution of an imaging target internal light source requiring solution and is non-negative;
step 3.4: the system matrix M (λ) in step 3.3 is processed as follows:
the linear relationship in step 3.3 becomesWhereinIs the processed system matrix.
Specifically, the method for converting the obtained linear relationship into the optimization problem in step 4 is as follows:
to be in step 3.4The following treatments were carried out:
specifically, the specific process of step 5 includes:
step 5.1: setting the feasible region Z of the initial light source as the whole imaging target, and setting a region contraction factorWherein N is I For nodes contained within the whole imaged object, N F Is the number of non-0 solutions at the end of the iteration, L max The number of iterative shrinkage of the region;
step 5.2: in conjunction with the feasible region Z, the optimization problem of step 4 becomesWherein the content of the first and second substances,obtaining an optimal solution X (Z) under the feasible region Z by adopting an algebraic iterative algorithm ART for a system matrix under the feasible region Z;
step 5.3: the optimal solution X (Z) under the feasible region Z obtained according to the step 5.2 and a system matrix under the feasible regionCalculating the luminous flux phi of the imaged target surface c (lambda); calculating multi-spectral fluorescence data of a penetrating imaged objectWith luminous flux phi c A difference of (λ), i.e.Wherein i represents the number of iterations (0 < i.ltoreq.L) max )；
Step 5.4: comparing f (i) obtained in step 5.3 with min (f (1), \8230;, f (i-1)), and if f (i) is less than min (f (1), \8230;, f (i-1)), solving the global optimal solution X * (Z) replacing with X (Z), otherwise, not replacing;
step 5.5: sorting the optimal solution X (Z) obtained in the step 5.2 from large to small, and selecting a node corresponding to the previous length (Z)/beta value to construct a new feasible region Z;
step 5.6: judging whether the iteration number i reaches the maximum value L max If yes, stopping the iterative contraction of the area and outputting a global optimal solution X * (Z); if not, go to step 5.2;
step 5.7: from the global optimum solution X * (Z) obtaining an optimal solution X within the entire imaging target * If the node is located in the feasible region Z, X * ＝X * (Z); if the node is located outside the feasible region Z, X * ＝0。
Specifically, the process of step 6 includes:
step 6.1: selecting X * The node corresponding to the element other than 0 in (b) constructs a similarity matrix D, where D is a symmetric matrix D (i, k) = D (k, i), i.e., the node corresponding to the element other than 0 in (b)
WhereinWherein the coordinates of the ith node are(x i ,y i ,z i ) The coordinates of the kth node are (x) k ,y k ,z k )，X * i Is the energy of the ith node, X * k Is the energy of the kth node;
step 6.2: recording an initial reference matrix P as 0, an initial adaptive matrix A as 0, an initial representative matrix R as 0, and setting clustering iteration times;
step 6.3: calculating a representative matrix R ok1 And the adaptive matrix A ok1 The calculation formula is as follows:
wherein R is ok1 (i, k) represents the degree of suitability of the point k as a representative point of the class of the point i, A ok1 (i, k) represents the suitability of the point i as a class representative point of the point k, k 'satisfies k' ≠ k, i 'satisfies i' ≠ k;
step 6.4: calculating an updated representative matrix R and an appropriate matrix A, wherein the alternative updating process comprises the following steps:
wherein R is n+1 Representing the representative matrix R obtained in the (n + 1) th iteration n+1 ok1 The updated representative matrix; r n Representing the representative matrix R obtained in the nth iteration n ok1 The updated representative matrix; a. The n+1 Representing the representative matrix A obtained from the (n + 1) th iteration n+1 ok1 The updated representative matrix; a. The n Representing the representative matrix A obtained from the nth iteration n ok1 The updated representative matrix; r represents a damping coefficient;
step 6.5: judging whether the clustering iteration times reach the set clustering iteration times, if so, terminating the iteration, and outputting a representative matrix R and an adaptive matrix A, otherwise, turning to the step 6.3;
step 6.6: classifying each node according to the representative matrix R and the adaptive matrix A, and determining the clustering center of each class, wherein the specific method comprises the following steps:
for any node i, if arg is satisfied k max { A (i, k) + R (i, k) }, then node i and node k belong to the same class; and calculating A (m, m) + R (m, m) for any node m in each class, and selecting the node m corresponding to the maximum value in A (m, m) + R (m, m) as the clustering center.
Specifically, the process of step 7 includes:
step 7.1: calculating energy values of all tetrahedrons according to the corresponding relation between the grid nodes and the tetrahedrons obtained in the step 3, selecting the tetrahedron with the energy value being not 0 as a clustering object, and calculating the geometric center of the clustering object;
step 7.2: sequentially calculating the distance between the geometric center of the tetrahedron with each energy value being not 0 and each initial clustering center, { l 1 ,…l c Where c is the number of initial cluster centers according to the minimum distance l i Dividing the tetrahedron into ith classes;
step 7.3: recalculating the cluster center C for each class:
where p (i), a (i), c (i) are the energy, volume, and geometric centers of each type of ith non-0 tetrahedron.
Compared with the prior art, the invention has the following technical effects:
1. the invention utilizes the optical characteristic parameters and the anatomical structure information as the prior knowledge, combines the multispectral fluorescence data, and improves the accuracy of the BLT light source reconstruction and the resolution of the reconstructed image.
2. The invention adopts a self-adaptive feasible region iterative contraction strategy, reduces the ill-conditioned performance of the BLT reconstruction problem, and reduces the dependency of the reconstruction result on the reconstruction algorithm.
3. The method converts the determination problem of the reconstructed light source center into the cluster analysis problem, adopts the AP clustering and the k-means clustering algorithm, effectively utilizes the advantage that the AP clustering does not need to appoint the number of clusters in advance, and the relative scalability of the k-means clustering algorithm, and simultaneously avoids the instability of the AP clustering and the limitation that the k-means clustering must know the cluster mean value; in the hybrid clustering process, not only the energy distribution information of the grid nodes is directly used, but also information such as reconstruction energy, center coordinates, volume and the like of a tetrahedron to which the nodes belong is fully considered, so that the stability of a clustering algorithm is improved to the greatest extent, and the complexity of computing time is reduced.
Drawings
FIG. 1 is a flow chart of a method of the present invention;
FIG. 2 (a) is a simulation model and a light source distribution diagram with 10mm center distance;
FIG. 2 (b) is a diagram showing a simulation experiment model and a light source distribution diagram with a center distance of 5 mm;
FIG. 3 (a) is a graph showing the light intensity distribution of a biological surface measured by a detector with the center thereof spaced by 10 mm;
FIG. 3 (b) is a graph showing the light intensity distribution of the biological surface measured by the detector with the center thereof spaced 5 mm;
FIG. 4 (a) is a reconstructed light source distribution diagram of the present invention with 10mm light source centers apart;
FIG. 4 (b) is a diagram of the distribution of the cluster centers of the light sources in the case of 10mm distance between the centers of the light sources according to the present invention;
FIG. 5 (a) is a reconstructed light source distribution diagram of the present invention with the light source centers separated by 5 mm;
FIG. 5 (b) is a diagram of the distribution of the cluster centers of the light sources with a distance of 5mm between the centers of the light sources according to the present invention;
the embodiments of the invention will be explained and explained in further detail with reference to the figures and the detailed description.
Detailed Description
Examples
According to the above technical solution, referring to fig. 1, the bioluminescence tomography reconstruction method based on multi-light source resolution of the embodiment specifically includes the following steps:
step 1: obtaining anatomical structure information, optical specificity parameters and surface multispectral fluorescence data of an imaging target
The study was carried out on a subcutaneous lesion, the imaging target being defined as skin tissue of a certain thickness, regardless of the border, and a cylinder of radius 42mm and height 24mm was generated for simulation, see document [ 1 ].
The imaging target is subjected to CT scanning to obtain anatomical structure information of the imaging target, such as shape, size, position and other information.
The above-mentioned imaging target is fixed on an electrically controlled rotary table, an optical detection instrument is placed above the imaging target, a light source generated inside the imaging target emits fluorescence whose wavelength is determined by the biological probe material itself and whose wavelength range is 590 to 650nm, and the optical specificity parameter of the transmission of the above-mentioned fluorescence in the skin is set with reference to the document [ 1 ].
The fluorescence penetrates the imaged object and the optical detection instrument detects the multispectral fluorescence data penetrating the imaged object, which does not collect the fluorescence data of the side and bottom of the imaged object, and is therefore called single angle data collection.
The optical detection instrument adopts a CCD camera with high sensitivity.
Step 2: normalizing the surface multispectral fluorescence data obtained in the step 1
Because the fluorescence absorption and scattering properties of different spectra are different, the fluorescence intensity transmitted to the surface of the imaging target has great difference, and in order to eliminate the influence of the difference on the reconstruction result, the surface multispectral fluorescence data of each wave band is normalized:
wherein phi (lambda) is surface fluorescence data with the wavelength of lambda,normalized surface fluorescence data at wavelength λ.
And step 3: and constructing a linear relation between the surface fluorescence data and the distribution of the internal light source of the imaging target based on a light transmission model and a finite element method. The specific implementation mode is as follows:
step 3.1: and (3) obtaining a light transmission model according to the anatomical structure information and the optical specificity parameters obtained in the step (1), wherein the light transmission model is a diffusion approximation equation, and the transmission process of light in the imaging target is described by using the diffusion approximation equation.
Step 3.2: and carrying out grid discretization on the imaging target by using a nirfast tool to obtain discrete grid data of the imaging target, wherein the discrete grid data comprises grid nodes and tetrahedrons, and the corresponding relation between the grid nodes and the tetrahedrons is obtained.
Step 3.3: and (3) dispersing the diffusion approximation equation in the step (3.1) according to the discrete grid data and a finite element theory, and further constructing a linear relation between the surface fluorescence data and the distribution of the internal light source of the imaging target:
Φ(λ)＝M(λ)X
wherein M (lambda) is a system matrix representing the absorption and scattering effects of all units in the imaging target on light, and X is the energy distribution of the light source inside the imaging target requiring solution and is non-negative.
Step 3.4: since the normalization processing is performed on the surface fluorescence data under different spectra in the step 2, the system matrix M in the step 3.3 is processed as follows:
the linear relationship in step 3.3 also becomesWhereinIs the processed system matrix.
And 4, step 4: the linear relationship obtained in step 3.4 is converted into an optimization problem:
and 5: performing algebraic iterative solution on the optimization problem in the step 4 by adopting a self-adaptive feasible region iterative contraction strategy to obtain an optimal solution X * . The specific implementation method comprises the following steps:
step 5.1: setting the feasible region Z of the initial light source as the whole imaging target and setting the region shrinkage factorWherein N is I For nodes contained within the whole imaged object, N F The number of non-0 solutions at the end of the iteration is set to 1,L max The number of iterative contractions for the region was set to 50.
And step 5.2: in conjunction with the feasible region Z, the optimization problem of step 4 becomesWherein, the first and the second end of the pipe are connected with each other,for a system matrix under the feasible region Z, an algebraic iterative Algorithm (ART) is specifically adopted for solving to obtain an optimal solution X (Z) under the feasible region Z, and the optimal solution X (Z) reflects the energy of each node under the feasible region Z.
Step 5.3: the optimal solution X (Z) under the feasible region Z obtained according to the step 5.2 and a system matrix under the feasible regionCalculating the luminous flux phi of the imaged target surface c (lambda). Calculating a penetration imaging target detected by an optical detection instrumentMulti-spectral fluorescence data ofWith luminous flux phi c (λ) difference, i.e.Wherein i represents the number of iterations (0 < i.ltoreq.L) max )。
Step 5.4: comparing f (i) obtained in step 5.3 with min (f (1), \8230;, f (i-1)), and if f (i) is less than min (f (1), \8230;, f (i-1)), determining the global optimal solution X * (Z) is replaced with X (Z), otherwise, no replacement is made. .
Step 5.5: and (3) sequencing the optimal solution X (Z) obtained in the step (5.2) from large to small, and selecting nodes corresponding to the previous length (Z)/beta values to construct a new feasible region Z.
Step 5.6: judging whether the iteration number i reaches the maximum value L max If yes, stopping the iterative contraction of the area and outputting a global optimal solution X * (Z); if not, go to step 5.2.
Step 5.7: from the global optimal solution X * (Z) obtaining an optimal solution X within the entire imaging target * If the node is located in the feasible region Z, X * ＝X * (Z); if the node is outside the feasible region Z, X * ＝0。
Optimal solution X * Representing the energy distribution of the individual nodes within the whole imaged object, which is shown on the spatial distribution map, i.e. enabling the reconstruction of the light source generated inside the imaged object.
In order to accurately acquire the central position of the light source generated inside the imaging target, the method of the invention further comprises the following steps:
and 6: the optimal solution X obtained in the step 5 is * And (4) carrying out AP center clustering to obtain the clustering center of each type. Due to X * Each element corresponds to the energy of the grid node, and X is selected when the central clustering is performed * The nodes corresponding to the non-0 elements in the (1) are clustered according to the spatial distribution, and the AP clustering can automatically determine the clustering number and the clustering center, so that the method adoptsAnd the AP clustering algorithm is used for preliminarily positioning the centers of all the light sources on the basis of determining the number of the light sources.
The specific implementation method comprises the following steps:
step 6.1: selecting X * The nodes corresponding to the non-0 elements in (b) construct a similarity matrix D, where D is a symmetric matrix D (i, k) = D (k, i), i.e., D is a symmetric matrix D (i, k) = D (k, i)
WhereinWherein the coordinate of the ith node is (x) i ,y i ,z i ) The coordinates of the kth node are (x) k ,y k ,z k )，X * i Is the energy of the ith node, X * k Is the energy of the kth node.
Step 6.2: the initial reference degree matrix P is recorded as 0, the initial adaptive matrix A is recorded as 0, the initial representative matrix R is recorded as 0, and the clustering iteration times are recorded as 1000.
Step 6.3: the AP clustering algorithm delivers two types of message matrices: represents the matrix R ok1 And the adaptive matrix A ok1 Reference [ 2 ], the calculation formula is as follows:
wherein R is ok1 (i, k) represents the suitability of the point k as a representative point of the class of the point i. A. The ok1 (i, k) represents the degree of suitability of the point i as the class representative point of the point k. k 'satisfies k' ≠ k, i 'satisfies i' ≠ k.
Step 6.4: the iterative process of the AP clustering algorithm is the process of alternately updating the two message matrixes R and A, the damping coefficient R is set to be 0.8, the updated representative matrix R and the adaptive matrix A are calculated, and the alternate updating process is as follows:
wherein R is n+1 Representing the representative matrix R obtained from the (n + 1) th iteration n+1 ok1 The updated representative matrix; r is n Representing the representative matrix R obtained in the nth iteration n ok1 The updated representative matrix; a. The n+1 Represents the representative matrix A obtained by the (n + 1) th iteration n+1 ok1 The updated representative matrix; a. The n Representing the representative matrix A obtained in the nth iteration n ok1 The updated representative matrix;
step 6.5: and judging whether the clustering iteration times reach the set clustering iteration times, if so, terminating the iteration, and outputting a representative matrix R and an adaptive matrix A, otherwise, turning to the step 6.3.
Step 6.6: determining a clustering center m according to the representative matrix R and the adaptive matrix A
For any node i, if arg is satisfied k max { A (i, k) + R (i, k) }, then node i and node k belong to the same class; and calculating A (m, m) + R (m, m) for any node m in each class, and selecting the node m corresponding to the maximum value in A (m, m) + R (m, m) as the clustering center.
And 7: and (4) taking the clustering center obtained by AP clustering in the step (6) as an initial clustering center of a k-means clustering algorithm, and carrying out k-means clustering to obtain a clustering center C of each type. Since the AP clustering result has poor stability due to sparsity of the optimal solution, k-means clustering is carried out according to the geometric center of the energy non-0 tetrahedron. The concrete implementation method comprises the following steps:
step 7.1: calculating energy values of all tetrahedrons according to the corresponding relation between the grid nodes and the tetrahedrons obtained in the step 3, selecting the tetrahedron with the energy value being not 0 as a clustering object, and calculating the geometric center of the clustering object; the design can fully utilize the correlation between the grid nodes and the tetrahedron to obtain more non-zero data for clustering, and the maximum data is obtainedTo a degree reduced due to the optimal solution X * The clustering is performed by adopting tetrahedral spatial distribution, so that the reliability and stability of a clustering data source are improved, the accuracy of clustering analysis is improved, and the test error is reduced.
And 7.2: sequentially calculating the distance l between the geometric center of each tetrahedron with energy value different from 0 and each initial clustering center 1 ,…l c Where c is the number of initial cluster centers according to the minimum distance l i The tetrahedra is classified into the ith class.
Step 7.3: recalculating the cluster center C for each class:
where p (i), a (i), c (i) are the energy, volume, and geometric center of each ith non-0 tetrahedron.
And 8: utilizing matlab related tool kit and nirfast, iso2mesh related function to obtain the optimal solution X in the step 5 * Performing geometric 3D display to realize the display of the reconstruction result of the light source inside the imaging target; and (4) performing geometric 3D display on the clustering center C obtained in the step (7) to determine the central position of the light source in the imaging target.
Analysis of Experimental results
The reconstruction results and the clustering centers of the present invention are further described with reference to fig. 2, fig. 3, and fig. 4.
A mimic used in the experiment is shown in FIG. 2. It is a cylinder with a radius of 42mm and a height of 24mm, and mainly simulates the pathological changes of deep subcutaneous tissues. Two sets of mock experiments were set up with centers 10mm and 5mm apart in the experiment, where FIG. 2 (a) simulates 5 light sources with centers 10mm apart and FIG. 2 (b) simulates 5 light sources with centers 5mm apart.
The surface intensity profile required for the reconstruction process is shown in figure 3. The surface light intensity distributions of fig. 3 were all collected on the upper surface of the phantom due to the simulated skin tissue. FIG. 3 (a) is the measurement data of 5 light source surfaces with a center distance of 10 mm. FIG. 3 (b) is the measurement data of 5 light source surfaces with a center distance of 5 mm. The influence of the reduction of the distance between the centers on the reconstruction difficulty can be known by comparison.
The final multi-light source reconstruction results are shown in fig. 4 and 5. Wherein, fig. 4 shows the experimental reconstruction result and clustering result based on the present invention under the condition that the light source centers are 10mm apart. FIG. 5 shows the experimental reconstruction and clustering results based on the present invention when the light source centers are 5mm apart. Fig. 4 (a) and 5 (a) show the reconstruction results, the red region is the reconstructed light source distribution, the blue circle represents the real light source center, and the blue cross represents the clustering result center. Fig. 4 (b) and 5 (b) show final clustering result graphs, and the central point is the final clustering center. In clinical application, the lesion area can be clearly obtained from fig. 4 (a), fig. 5 (a), fig. 4 (b) and fig. 5 (b), and the area center can be accurately determined. The reconstruction method based on the invention has the advantages of small error of the position of the reconstructed light source, stable algorithm and accurate positioning of the multi-target center, and is an effective multi-light-source bioluminescence tomography reconstruction method.
Claims (6)
1. A bioluminescence tomography reconstruction method based on multi-light source resolution is characterized by comprising the following steps:
step 1: obtaining anatomical structure information, optical specificity parameters and surface multispectral fluorescence data of an imaging target;
step 2: normalizing the surface multispectral fluorescence data obtained in the step (1) to obtain normalized surface fluorescence data;
and step 3: constructing a linear relation between surface fluorescence data and the distribution of an internal light source of an imaging target based on a light transmission model and a finite element method;
and 4, step 4: converting the linear relation obtained in the step 3 into an optimization problem;
and 5: performing algebraic iterative solution on the optimization problem in the step 4 by adopting a self-adaptive feasible region iterative contraction strategy to obtain an optimal solution X * Completing reconstruction;
the bioluminescence tomography reconstruction method based on multi-light source resolution further comprises the following steps:
step 6: for the optimal solution X obtained in the step 5 * Carrying out AP center clustering, and calculating to obtain a clustering center m of each type;
and 7: taking the clustering center obtained in the step 6 as an initial clustering center of a k-means clustering algorithm, carrying out k-means clustering, and calculating to obtain a clustering center C of each type;
and step 8: the optimal solution X obtained in the step 5 is * Performing geometric 3D display to realize the display of the reconstruction result of the light source inside the imaging target; and (4) performing geometric 3D display on the clustering center C obtained in the step (7) to determine the central position of the light source in the imaging target.
2. The multi-light-source resolution-based bioluminescence tomography reconstruction method according to claim 1, wherein the specific process of the step 3 comprises:
step 3.1: obtaining an optical transmission model according to the anatomical structure information and the optical specificity parameters obtained in the step 1, wherein the optical transmission model is a diffusion approximation equation, and the transmission process of light in the imaging target is described by using the diffusion approximation equation;
step 3.2: carrying out grid dispersion on the imaging target by using a nirfast tool to obtain discrete grid data of the imaging target, wherein the discrete grid data comprises grid nodes and tetrahedrons, and the corresponding relation between the grid nodes and the tetrahedrons is obtained;
step 3.3: and (3) dispersing the diffusion approximation equation in the step (3.1) according to the discrete grid data and a finite element theory, and further constructing a linear relation between the surface fluorescence data and the distribution of the internal light source of the imaging target:
Φ(λ)＝M(λ)X
wherein phi (lambda) is surface fluorescence data with the wavelength of lambda, M (lambda) represents a system matrix of the fluorescence data under the wavelength of lambda, and X is the energy distribution of an imaging target internal light source requiring solution and is non-negative;
step 3.4: the system matrix M (λ) in step 3.3 is processed as follows:
the linear relationship in step 3.3 becomesWhereinIs the processed system matrix.
3. The multi-light-source resolution-based bioluminescence tomography reconstruction method according to claim 2, wherein the linear relationship obtained in the step 4 is converted into an optimization problem by:
to be in step 3.4The following treatments were carried out:
4. the multi-light-source resolution-based bioluminescence tomography reconstruction method according to claim 3, wherein the specific process of the step 5 comprises:
step 5.1: setting the feasible region Z of the initial light source as the whole imaging target, and setting a region contraction factorWherein N is I For nodes contained within the whole imaged object, N F Number of solutions, L, that are not 0 at the end of the iteration max The number of times of iterative shrinkage of the region;
step 5.2: in conjunction with the feasible region Z, the optimization problem of step 4 becomesWherein the content of the first and second substances,obtaining an optimal solution X (Z) under the feasible region Z by adopting an algebraic iterative algorithm ART for a system matrix under the feasible region Z;
step 5.3: the optimal solution X (Z) under the feasible region Z obtained according to the step 5.2 and a system matrix under the feasible regionCalculating the luminous flux phi of the imaged target surface c (lambda); calculating multi-spectral fluorescence data of a penetrating imaged objectWith luminous flux phi c (λ) difference, i.e.Wherein i represents the number of iterations (0 < i.ltoreq.L) max )；
Step 5.4: comparing f (i) obtained in step 5.3 with min (f (1), \8230;, f (i-1)), and if f (i) is less than min (f (1), \8230;, f (i-1)), determining the global optimal solution X * (Z) replacing with X (Z), otherwise, not replacing;
step 5.5: sorting the optimal solution X (Z) obtained in the step 5.2 from large to small, and selecting nodes corresponding to the previous length (Z)/beta values to construct a new feasible region Z;
step 5.6: judging whether the iteration number i reaches the maximum value L max If yes, stopping the iterative contraction of the area and outputting a global optimal solution X * (Z); if not, go to step 5.2;
step 5.7: from the global optimal solution X * (Z) obtaining an optimal solution X within the entire imaging target * If the node is located in the feasible region Z, X * ＝X * (Z); if the node is located outside the feasible region Z, X * ＝0。
5. The multi-light-source-resolution-based bioluminescent tomographic reconstruction method of claim 4, wherein the process of step 6 comprises:
step 6.1: selecting X * The node corresponding to the element other than 0 in (b) constructs a similarity matrix D, where D is a symmetric matrix D (i, k) = D (k, i), i.e., the node corresponding to the element other than 0 in (b)
WhereinWherein the coordinates of the ith node are (x) i ,y i ,z i ) The coordinates of the kth node are (x) k ,y k ,z k )，X * i Is the energy of the ith node, X * k Energy of the kth node;
step 6.2: recording an initial reference matrix P as 0, an initial adaptive matrix A as 0, an initial representative matrix R as 0, and setting clustering iteration times;
step 6.3: calculating a representative matrix R ok1 And the adaptive matrix A ok1 The calculation formula is as follows:
wherein R is ok1 (i, k) represents the degree of suitability of point k as a representative point of the class of point i, A ok1 (i, k) represents the suitability of the point i as a class representative point of the point k, k 'satisfies k' ≠ k, i 'satisfies i' ≠ k;
step 6.4: calculating an updated representative matrix R and an appropriate matrix A, wherein the alternative updating process comprises the following steps:
wherein R is n+1 Representing the representative matrix R obtained in the (n + 1) th iteration n+1 ok1 The updated representative matrix; r n Representing the representative matrix R obtained in the nth iteration n ok1 The updated representative matrix; a. The n+1 Representing the representative matrix A obtained from the (n + 1) th iteration n+1 ok1 The updated representative matrix; a. The n Representing the representative matrix A obtained in the nth iteration n ok1 The updated representative matrix; r represents a damping coefficient;
step 6.5: judging whether the clustering iteration times reach the set clustering iteration times, if so, terminating the iteration, and outputting a representative matrix R and an appropriate matrix A, otherwise, turning to the step 6.3;
step 6.6: classifying each node according to the representative matrix R and the adaptive matrix A, and determining the clustering center of each class, wherein the specific method comprises the following steps:
for any node i, if arg is satisfied k max { A (i, k) + R (i, k) }, then node i and node k belong to the same class; and calculating A (m, m) + R (m, m) for any node m in each class, and selecting the node m corresponding to the maximum value in A (m, m) + R (m, m) as the clustering center.
6. The multi-light-source resolution-based bioluminescent tomographic reconstruction method of claim 5, wherein the step 7 process comprises:
step 7.1: calculating energy values of all tetrahedrons according to the corresponding relation between the grid nodes and the tetrahedrons obtained in the step 3, selecting the tetrahedron with the energy value being not 0 as a clustering object, and calculating the geometric center of the clustering object;
step 7.2: sequentially calculating the distance between the geometric center of the tetrahedron with each energy value being not 0 and each initial clustering center, { l 1 ,…l c Where c is the number of initial cluster centers according to the minimum distance l i Dividing the tetrahedron into ith classes;
step 7.3: recalculating the cluster center C for each class:
where p (i), a (i), c (i) are the energy, volume, and geometric centers of each type of ith non-0 tetrahedron.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510589426.XA CN105326475B (en) | 2015-09-16 | 2015-09-16 | A kind of Bioluminescence tomography reconstruction method differentiated based on multiple light courcess |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510589426.XA CN105326475B (en) | 2015-09-16 | 2015-09-16 | A kind of Bioluminescence tomography reconstruction method differentiated based on multiple light courcess |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105326475A CN105326475A (en) | 2016-02-17 |
| CN105326475B true CN105326475B (en) | 2018-01-19 |
Family
ID=55277405
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510589426.XA Active CN105326475B (en) | 2015-09-16 | 2015-09-16 | A kind of Bioluminescence tomography reconstruction method differentiated based on multiple light courcess |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105326475B (en) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106404730B (en) * | 2016-08-30 | 2019-10-11 | 上海大学 | The description method that light based on Lattice Boltzmann method model is propagated in the medium |
| CN107184181A (en) * | 2017-05-15 | 2017-09-22 | 清华大学 | The processing method and system of Dynamic Fluorescence molecular tomographic |
| CN107374588B (en) * | 2017-08-01 | 2020-05-19 | 西北大学 | Multi-light-source fluorescent molecular tomography reconstruction method based on synchronous clustering |
| CN108095686B (en) * | 2017-11-06 | 2020-09-11 | 西北大学 | Fluorescence molecular tomography target feasible region selection method |
| CN107713995B (en) * | 2017-11-15 | 2020-09-01 | 陕西师范大学 | Bioluminescence tomography light source reconstruction method based on penalty algorithm |
| CN108444980B (en) * | 2018-01-30 | 2020-08-07 | 中国科学院上海技术物理研究所 | L IBS quantitative analysis method based on algebraic reconstruction relevance vector solution |
| CN108451508B (en) * | 2018-04-28 | 2020-05-05 | 中国科学院自动化研究所 | Biological autofluorescence three-dimensional imaging method based on multilayer perceptron |
| CN112229827B (en) * | 2020-09-07 | 2022-02-08 | 南京大学 | Real-time multispectral tomography method and device |
| CN112089434B (en) * | 2020-10-16 | 2024-05-03 | 陕西师范大学 | Multispectral bioluminescence tomography method and system |
| CN113781652B (en) * | 2021-08-16 | 2024-04-12 | 西北大学 | Multi-level probability reconstruction method based on energy density region shrinkage |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101292863A (en) * | 2007-04-25 | 2008-10-29 | 中国科学院自动化研究所 | Self-adaption finite element light source reconstruction method based on simple spectrum or mix segment measurement |
| CN101342075A (en) * | 2008-07-18 | 2009-01-14 | 北京工业大学 | Multi-optical spectrum autofluorescence dislocation imaging reconstruction method based on single view |
| CN101396262A (en) * | 2008-10-31 | 2009-04-01 | 清华大学 | Fluorescent molecule tomography rebuilding method based on linear relationship |
| CN102034266A (en) * | 2010-11-30 | 2011-04-27 | 中国科学院自动化研究所 | Rapid sparse reconstruction method and equipment for exciting tomography fluorescence imaging |
| CN102334979A (en) * | 2011-08-03 | 2012-02-01 | 中国科学院自动化研究所 | Bimodal fusion tomography method based on iterative shrinkage |
| CN103300829A (en) * | 2013-06-25 | 2013-09-18 | 中国科学院自动化研究所 | Biological autofluorescence tomography method based on iteration reweighting |
| CN103393410A (en) * | 2013-08-21 | 2013-11-20 | 西安电子科技大学 | Fluorescence molecular tomography reconstruction method based on alternative iterative operation |
| CN104586366A (en) * | 2015-01-27 | 2015-05-06 | 西安电子科技大学 | Fluorescence molecular tomography reconstruction method |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6666579B2 (en) * | 2000-12-28 | 2003-12-23 | Ge Medical Systems Global Technology Company, Llc | Method and apparatus for obtaining and displaying computed tomography images using a fluoroscopy imaging system |
| US8237786B2 (en) * | 2009-12-23 | 2012-08-07 | Applied Precision, Inc. | System and method for dense-stochastic-sampling imaging |
-
2015
- 2015-09-16 CN CN201510589426.XA patent/CN105326475B/en active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101292863A (en) * | 2007-04-25 | 2008-10-29 | 中国科学院自动化研究所 | Self-adaption finite element light source reconstruction method based on simple spectrum or mix segment measurement |
| CN101342075A (en) * | 2008-07-18 | 2009-01-14 | 北京工业大学 | Multi-optical spectrum autofluorescence dislocation imaging reconstruction method based on single view |
| CN101396262A (en) * | 2008-10-31 | 2009-04-01 | 清华大学 | Fluorescent molecule tomography rebuilding method based on linear relationship |
| CN102034266A (en) * | 2010-11-30 | 2011-04-27 | 中国科学院自动化研究所 | Rapid sparse reconstruction method and equipment for exciting tomography fluorescence imaging |
| CN102334979A (en) * | 2011-08-03 | 2012-02-01 | 中国科学院自动化研究所 | Bimodal fusion tomography method based on iterative shrinkage |
| CN103300829A (en) * | 2013-06-25 | 2013-09-18 | 中国科学院自动化研究所 | Biological autofluorescence tomography method based on iteration reweighting |
| CN103393410A (en) * | 2013-08-21 | 2013-11-20 | 西安电子科技大学 | Fluorescence molecular tomography reconstruction method based on alternative iterative operation |
| CN104586366A (en) * | 2015-01-27 | 2015-05-06 | 西安电子科技大学 | Fluorescence molecular tomography reconstruction method |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105326475A (en) | 2016-02-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105326475B (en) | A kind of Bioluminescence tomography reconstruction method differentiated based on multiple light courcess | |
| US8094149B2 (en) | Multi-spectral reconstruction method based on adaptive finite element | |
| CN101342075B (en) | Multi-optical spectrum autofluorescence dislocation imaging reconstruction method based on single view | |
| CN107924457A (en) | For the area-of-interest in lookup hematoxylin and the organization chart picture of eosin (H ＆ E) dyeing in multiplexing/super composite fluorescence organization chart picture and quantify the system and method for intra-tumor cell spaces heterogeneity | |
| CN102940482B (en) | Adaptive tomographic fluorescence imaging (TFI) reconstructing method | |
| CN102334979B (en) | Bimodal fusion tomography method based on iterative shrinkage | |
| CN102034266B (en) | Rapid sparse reconstruction method and equipment for exciting tomography fluorescence imaging | |
| WO2017004851A1 (en) | Bioluminescence tomography reconstruction algorithm based on multi-task bayesian compressed sensing method | |
| Zhang et al. | Robust reconstruction of fluorescence molecular tomography based on sparsity adaptive correntropy matching pursuit method for stem cell distribution | |
| CN102988026A (en) | Auto-fluorescence tomography re-establishing method based on multiplier method | |
| He et al. | Half thresholding pursuit algorithm for fluorescence molecular tomography | |
| CN107146261B (en) | Bioluminescence tomography quantitative reconstruction method based on magnetic resonance image prior region of interest | |
| CN103300829A (en) | Biological autofluorescence tomography method based on iteration reweighting | |
| CN105581779A (en) | Bioluminescent fault imaging reestablishment method for directly fusing structure imaging | |
| Slavine et al. | Iterative reconstruction method for light emitting sources based on the diffusion equation | |
| CN109820479B (en) | Fluorescence molecular tomography feasible region optimization method | |
| CN108095686A (en) | A kind of fluorescent molecular tomography target feasible zone choosing method | |
| CN108140236A (en) | For normalizing image processing method and device with artifact correction | |
| CN115868923A (en) | Fluorescence molecule tomography method and system based on expanded cyclic neural network | |
| Li et al. | Discover mouse gene coexpression landscapes using dictionary learning and sparse coding | |
| CN112089434B (en) | Multispectral bioluminescence tomography method and system | |
| CN113951835B (en) | Three-dimensional fluorescence microscopic imaging method based on optical fault reconstruction strategy | |
| CN113288188B (en) | Cone beam X-ray luminescence tomography method based on grouping attention residual error network | |
| Yi et al. | A permissible region extraction based on a knowledge priori for x-ray luminescence computed tomography | |
| CN107374588B (en) | Multi-light-source fluorescent molecular tomography reconstruction method based on synchronous clustering |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |