CN113066142B - Optical function imaging method combining spatial regularization and semi-blind spectrum unmixing - Google Patents

Optical function imaging method combining spatial regularization and semi-blind spectrum unmixing Download PDF

Info

Publication number
CN113066142B
CN113066142B CN202110205426.0A CN202110205426A CN113066142B CN 113066142 B CN113066142 B CN 113066142B CN 202110205426 A CN202110205426 A CN 202110205426A CN 113066142 B CN113066142 B CN 113066142B
Authority
CN
China
Prior art keywords
chromophore
spectrum
iteration
concentration distribution
spectral
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
Application number
CN202110205426.0A
Other languages
Chinese (zh)
Other versions
CN113066142A (en
Inventor
王艺涵
许梦璐
朱守平
昝克华
贺礼
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xidian University
Original Assignee
Xidian University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Xidian University filed Critical Xidian University
Priority to CN202110205426.0A priority Critical patent/CN113066142B/en
Publication of CN113066142A publication Critical patent/CN113066142A/en
Application granted granted Critical
Publication of CN113066142B publication Critical patent/CN113066142B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00Measuring for diagnostic purposes; Identification of persons
    • A61B5/0059Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence
    • A61B5/0075Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence by spectroscopy, i.e. measuring spectra, e.g. Raman spectroscopy, infrared absorption spectroscopy
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00Image generation
    • G06T2211/40Computed tomography
    • G06T2211/424Iterative

Abstract

The invention belongs to the technical field of biomedical multispectral optical function imaging, and discloses an optical function imaging method combining spatial regularization and semi-blind spectral unmixing, which inputs a multi-wavelength optical absorption coefficient map μa The method comprises the steps of carrying out a first treatment on the surface of the Initializing spectral characteristics and concentration distribution of each chromophore in target biological tissues; constructing an optimized objective function form, introducing a matched priori spectrum, and adding sparse regularization items of the concentration spatial distribution of each chromophore; inverting the objective function by using the proposed iterative solution framework, and iteratively updating unknown in-vivo spectral characteristics of the exogenous contrast agent in the biological tissue and concentration distribution images of each chromophore; judging whether an iteration stop condition is reached; otherwise, stopping iteration and simultaneously outputting an in-vivo spectrum of the exogenous contrast agent and a concentration distribution image of each chromophore. The invention can avoid the selection of the step length of the searching direction in the traditional gradient descent iteration process, and ensure the convergence of the iteration process and the non-negativity of the analysis result.

Description

Optical function imaging method combining spatial regularization and semi-blind spectrum unmixing
Technical Field
The invention belongs to the technical field of biomedical multispectral optical function imaging, and particularly relates to an optical function imaging method combining spatial regularization and semi-blind spectral unmixing.
Background
At present, cancer (malignant tumor) seriously threatens human health, and the development of the cancer is often accompanied by the change of physicochemical properties of tissues, including tissue structure, tissue composition, blood oxygen content and the like. The pathological change can cause the change of the optical property of the tissue, so that the optical function imaging technology can visualize the physiological parameter distribution in a non-invasive mode, further discover early lesions of the tissue, and achieve the purposes of reducing the death rate, improving the prognosis and improving the quality of life. The multispectral biomedical optical function imaging technology can quantitatively acquire the concentration distribution of each chromophore (namely functional parameters directly related to pathophysiology) in biological tissues by reconstructing a series of optical absorption distribution maps reflecting physiological characteristics of tissue bodies and utilizing a multispectral unmixing algorithm, so that the analysis of molecular mechanisms of disease formation and development is realized, and effective noninvasive detection and monitoring means can be provided for early diagnosis of diseases and treatment of the diseases. So the multispectral biomedical optical function imaging technology plays an important role in the biomedical research field.
The concentration distribution of each chromophore in biological tissues is analyzed by utilizing biomedical multispectral optical function imaging technology (such as multispectral diffusion optical tomography technology and multispectral photoacoustic tomography technology), namely the problem of spectral unmixing. Under a certain measurement wavelength, based on the linear condition that the optical absorption coefficient expressed by a certain pixel in the reconstructed image is a weighted sum of the molar extinction coefficient of each chromophore in the tissue body and the corresponding concentration value of each chromophore under the current measurement wavelength, the blood oxygen information and the exogenous contrast agent concentration information of the tissue region can be extracted by utilizing a multispectral unmixing algorithm, namely, the physiological and pathological information in the tissue body can be reflected qualitatively or quantitatively, and further the physiological processes of metabolism, blood flow and the like of the living body are reflected. The typical method for solving the multispectral component mixing problem is a spectrum fitting algorithm, namely a multispectral optical absorption distribution image obtained by processing the known spectrum characteristic information pixel by pixel under a least square rule frame, and the blood oxygen and exogenous contrast agent concentration distribution of a tissue vascular system is imaged by solving the least square problem. However, hemoglobin, lipids, melanin and water are the major endogenous chromophores in tissues, and their absorption spectra are known; while some exogenous contrast agents, such as indocyanine green (Indocyanine green, ICG), are only approved by the united states food and drug administration (Food and Drug Administration, FDA) as fluorescent dyes for use in humans, and have long been widely used in the detection of physiological states of various humans, the absorption spectrum of which is related to the nature of the solvent medium and the dye concentration; binding to plasma proteins results in a shift of the main peak in the absorption spectrum to a longer wavelength of about 25nm; furthermore, ICG solutions do not comply with lambert-beer's law at concentrations above 15mg/l in plasma due to the formation of aggregates with increasing concentration, and when the concentration is higher, the spectral profile exhibits two absorption peaks at λ=730 nm and λ=805 nm, respectively; at a dilution of 5mg/l, the molar extinction coefficient increases significantly at λ=805 nm, while the maximum at λ=730 nm becomes shoulder-shaped. The above problems make it impossible to accurately obtain the in-vivo spectral characteristics of the exogenous contrast agent, so that the unmixed performance of the spectral fitting algorithm is limited. In addition, the gradient descent method for optimizing and solving the least square problem usually needs to rely on an empirical method for selecting the search step length in the iterative updating process, and can not ensure the convergence and non-negativity of the analysis result.
Because the limitation of the spectrum fitting algorithm is mainly caused by the fact that the in-vivo spectrum characteristics of the exogenous contrast agent cannot be accurately acquired, blind solution mixing can be used as another method for analyzing the concentration distribution of the internal and external source chromophores in tissues, and the multispectral dataset is analyzed based on the data characteristics without the need of exact prior information about the absorption spectrum. Among the common blind unmixing algorithms are independent component analysis (Independent component analysis, ICA), vertex component analysis (Vertex component analysis, VCA), and the like. ICA uses the statistical nature of chromophore distribution to estimate the concentration of each component, and its applicable conditions require that the distribution of each chromophore within the tissue be statistically independent of each other. The VCA is based on the concept that observed data form a convex single-body in a feature space, and the absorption spectrum characteristics of each chromophore are searched through orthogonal projection, so that the concentration distribution of the chromophore is obtained, the calculation time of the chromophore is greatly reduced compared with ICA, and the VCA has important advantages when processing a large amount of data or performing real-time unmixing; however, for the multispectral optical reconstruction image with larger noise, the deviation between the absorption spectrum obtained by the VCA algorithm and the true value is larger, and the effectiveness is reduced; meanwhile, the initial projection vector of the single-body growth selected in each iteration of the algorithm has randomness, so that the VCA algorithm can have inconsistent running results. In contrast, the orthogonal projection-based object detection technique: an automatic target generation process (Automatic target generation process, ATGP), which is essentially the same as VCA, selects the determined initial projection vector to avoid randomness, and has a higher stability than VCA; in addition, the ATGP algorithm is not constrained by dimension reduction transformation, can ensure the integrity of data and prevent the characteristic information of an absorption spectrum from being destroyed, but the method is also sensitive to noise and can cause lower quantitative performance. In addition, statistical subpixel detection methods, such as adaptive matched filters (Adaptive matched filter, AMF), are proposed for multispectral optical functional imaging unmixing, and subpixel detection algorithms focus on identifying a spectral target that is significantly different from the background, as opposed to typical linear unmixing methods that attempt to resolve all of the different spectral components present in an image; however, this method has good performance only in the case of sparse distribution of target probes. Based on the analysis, in order to realize the synchronous reconstruction of unknown in-vivo spectral characteristics of exogenous contrast agents in biological tissues and concentration distribution of each chromophore, the method of the invention is required to have an optimized objective function and an efficient iterative solution framework.
Through the above analysis, the problems and defects existing in the prior art are as follows:
(1) In the prior art, the traditional spectrum blind unmixing algorithm is easily affected by noise; the use of spectral fitting algorithms is limited in cases where the in-vivo spectral characteristics of the exogenous contrast agent are unknown.
(2) Due to noise or other factors in the measurement, some conventional spectral unmixing algorithms often produce concentration negative values that are not physically significant, and cannot obtain an ideal chromophore concentration distribution image.
(3) When solving the objective function of the spectral unmixed problem, the selection of the search step length in the iterative updating process by the gradient descent method for optimizing and solving the least square problem usually needs to rely on an empirical method, and the convergence and non-negativity of the analysis result cannot be ensured.
The difficulty of solving the problems and the defects is as follows:
(1) Under the condition that the in-vivo spectrum of the exogenous contrast agent is unknown, the synchronous and accurate reconstruction of the in-vivo spectrum characteristics of the exogenous chromophore concentration distribution images and the exogenous contrast agent in each target tissue region is difficult to realize.
(2) In order to reconstruct the concentration distribution image of each chromophore with high fidelity, it is generally necessary to obtain the multi-wavelength optical absorption coefficient spectrum μ in advance, given that the spectral characteristics of each chromophore are known or have been effectively estimated a Enhancement processing is performed to reduce image noise or eliminate reconstruction artifacts due to undersampling, but in the prior art, these operations are preprocessing steps, increasing the time and computational cost of multispectral optical functional imaging.
(3) When solving the objective function of the spectral unmixed problem, in order to avoid selecting the search step length of iterative update in an empirical mode in the optimization calculation process, an algorithm for realizing the self-adaptive update of the search step length under the condition of optimizing the regularization constraint of the concentration image is invented, the solution process has convergence and non-negative property of the result, and the algorithm design has difficulty.
The meaning of solving the problems and the defects is as follows:
(1) According to the invention, by constructing an optimized solving objective function, introducing the spectrum of the inherent chromophore in the biological tissue as a priori spectrum in a multispectral unmixing reconstruction algorithm and adding spatial regularization constraint in image reconstruction, the reconstruction signal-to-noise ratio and fidelity of each internal and external chromophore concentration image and exogenous contrast agent spectrum can be improved. Avoiding the acquisition of the spectrum mu of the multi-wavelength optical absorption coefficient in advance a The enhancement processing of the method improves the algorithm efficiency, is expected to realize the concentration image calculation in near real time, and has significant significance in clinical application.
(2) The invention builds an efficient iterative solving framework, provides a searching step length self-adaptive updating calculation formula under the condition of optimizing the regularization constraint of the concentration image, and can realize the synchronization and the non-load construction of unknown in-vivo spectral characteristics of the exogenous contrast agent and the concentration images of each chromophore. The step of selecting the search step length in each iteration of the traditional optimization reconstruction algorithm is avoided while the convergence of the objective function process is guaranteed to be minimized. The invention meets the requirement of high robustness of the algorithm tool used in clinical application.
Disclosure of Invention
Aiming at the problems existing in the prior art, the invention provides an optical function imaging method combining spatial regularization and semi-blind spectrum unmixing. According to the invention, the spectrum of the inherent chromophore in the biological tissue is introduced as a priori spectrum, and the spatial regularization constraint is added in the image reconstruction, so that the effective matrix decomposition reconstruction framework is developed to carry out iterative solution on the objective function in consideration of the condition that the spectrum and the concentration of each chromophore are non-negative, thereby realizing the synchronous acquisition of the unknown in-vivo spectrum characteristics of the exogenous contrast agent in the biological tissue and the concentration distribution of each chromophore. In the iteration process, the method can avoid the selection of the step length of the search direction in the traditional gradient descent iteration process, and ensure the convergence of the iteration process and the non-negativity of the analysis result.
The invention is realized in such a way that an optical function imaging method combining spatial regularization and semi-blind spectrum unmixing is realized, the optical function imaging method combining spatial regularization and semi-blind spectrum unmixing strategies comprises the following steps:
step one, inputting a multi-wavelength optical absorption coefficient map mu a The method comprises the steps of carrying out a first treatment on the surface of the The multi-wavelength optical absorption spectrum of the biological tissue input in the first step can be directly obtained or reconstructed by using biomedical optical detection equipment (such as an endoscope, an optical mammary gland imager, a tissue biopsy microscope and the like).
Initializing the spectral characteristics and concentration distribution of each chromophore in the target biological tissue; in the initialization process, an ATGP algorithm is selected to preliminarily estimate the absorption spectrum of each chromophore in the tissue, and compared quantitatively, the ATGP algorithm is used as a pretreatment tool, so that compared with the VCA algorithm commonly used in other multispectral imaging, the method has stronger robustness; moreover, the ATGP algorithm is not constrained by dimension reduction transformation, so that the data integrity can be ensured, and the characteristic information of the absorption spectrum is prevented from being damaged. Secondly, replacing the spectrum of the corresponding chromophore analyzed from the ATGP pretreatment result by using the priori spectrum of the endogenous chromophore (hemoglobin, water, lipid, melanin and the like), thereby improving the accuracy of calculating the unknown in-vivo spectrum and concentration image of the exogenous contrast agent in the subsequent step;
Step three, constructing an optimized objective function form, introducing a matched priori spectrum, and adding sparse regularization items of the concentration spatial distribution of each chromophore; step three, introducing an accurately known priori spectrum, adding spatial constraint on concentration distribution of each chromophore by utilizing total variation regularization, and enabling the spectrum to be subjected to multi-wavelength optical absorption spectrum mu input a High quality reconstruction results are still obtained under the influence of noise or reconstruction artifacts.
Step four, inverting the objective function by using the proposed iterative solution framework, and iteratively updating the in-vivo spectral characteristics of the exogenous contrast agent in the biological tissue and the concentration distribution image of each chromophore; the step four can realize the synchronous reconstruction of unknown in-vivo spectral characteristics of the exogenous contrast agent and the concentration distribution of each chromophore, and simultaneously, the iterative updating formula provided in the step four ensures the convergence of the process of minimizing the objective function and avoids the step of selecting the searching step length in each iteration of the traditional optimization reconstruction algorithm.
Step five, judging whether an iteration stop condition is met, and returning to the step four if the iteration stop condition is not met; otherwise, stopping iteration and outputting the spectral characteristics of the exogenous contrast agent and the concentration distribution image of each chromophore. And fifthly, by setting reasonable iteration stopping criteria, the method can stably acquire high-quality reconstruction results, avoid spending a great deal of time and calculation cost, and meet the requirement of high robustness of the algorithm tool in clinical application.
Further, in the first step, a reconstructed multi-wavelength optical absorption coefficient map μ of the target tissue region is input a The method specifically comprises the following steps:
optical detection devices (e.g., optical endoscopes, optical breast imagers, tissue biopsy microscopes, etc.) use a near infrared first window wavelength range commonly used in biological tissue detectionPerforming multi-wavelength excitation and imaging on pathological tissues at a certain wavelength interval by using a light source in a second window wavelength range (600 nm-900 nm) or near infrared (1000 nm-1700 nm), wherein the image size can be set in detection equipment, for example, 35 excitation wavelengths are selected at a step length interval of 5nm in a wavelength range of 680 nm-884 nm, and an image with a size of 158×158 pixels under each wavelength is obtained; within this band, the biological tissue intrinsic chromophore reflecting the physiological change of the tumor region has unique spectral absorption characteristics; thereby obtaining the multispectral optical absorption coefficient mu of the target tissue region a And (5) a map.
In the second step, blind solution mixing pretreatment is performed by using an ATGP algorithm, a priori spectrum of a natural chromophore in a biological tissue is introduced, and the spectral characteristics and the concentration distribution of each chromophore in the target biological tissue are initialized.
Further, the specific process of the second step is as follows:
determining the total number W of chromophore components contained in biological tissue region, and mapping the optical absorption coefficient mu for multiple wavelengths a Realizing blind solution and obtaining preliminary estimation of absorption spectrum characteristics of each chromophore in tissue
Figure BDA0002950300100000061
Wherein->
Figure BDA0002950300100000062
The dimension is L×W->
Figure BDA0002950300100000063
L is the total number of wavelengths selected for measurement;
wherein L may be set to 35; since the tumor development and development process has close relation with tumor angiogenesis, hemoglobin is a main intrinsic chromophore in blood, and after fluorescent dye approved by FDA and having definite clinical application is added;
comprises three main absorbers and HbO in the near infrared spectrum 2 Hb, exogenous contrast agent such as ICG, so that the chromophore composition contained in a region of biological tissue is determinedThe total number W is 3;
calculation of absorption spectra of respective chromophores in tissue preliminarily estimated by ATGP
Figure BDA0002950300100000064
With intrinsic chromophore a priori spectral features
Figure BDA0002950300100000065
Spectral angle error between->
Figure BDA0002950300100000066
Measuring similarity between spectra; let P be the total number of intrinsic chromophore components in the biological tissue +.>
Figure BDA0002950300100000067
The dimension is L×P->
Figure BDA0002950300100000068
HbO 2 The Hb is taken as a natural chromophore in a target biological tissue area, the spectrum characteristic of the Hb is taken as a priori spectrum, and therefore P is 2;
Setting a spectral angle error threshold angle_tol, if ATGP preliminarily estimates the absorption spectrum of each chromophore in the tissue
Figure BDA0002950300100000071
And a priori spectrum->
Figure BDA0002950300100000072
The spectral angle error between them is smaller than the set threshold, i.e.>
Figure BDA0002950300100000073
The prior spectrum is used for replacing the corresponding chromophore spectrum analyzed in the pretreatment result so as to improve the accuracy of the follow-up method for calculating the unknown in-vivo spectrum and concentration of the exogenous contrast agent, and the step is finished to obtain the initialized spectrum characteristic epsilon of each chromophore in the target biological tissue;
the obtained multi-wavelength optical absorption coefficient spectrum mu a And organismsChromophore, hbO in tissue 2 Initializing spectral characteristics epsilon of Hb and ICG, and solving an initialized concentration distribution C of each chromophore in biological tissues by using a generalized inverse algorithm to complete all initialization steps:
Figure BDA0002950300100000074
further, said obtaining a preliminary estimate of the absorption spectra of these three types of chromophores
Figure BDA0002950300100000075
The method comprises the following steps:
the extracted nth target absorption spectrum t is extracted by utilizing ATGP in nth iteration n Satisfies the following formula:
Figure BDA0002950300100000076
in the method, in the process of the invention,
Figure BDA0002950300100000077
is an orthogonal projection matrix; i is an identity matrix; u (U) n-1 =[t 1 ,t 2 ,…,t n-1 ]A target absorption spectrum matrix generated for the n-1 th calculation; r is the vector of image pixels.
Further, the calculated spectral angle error formula is shown as follows:
Figure BDA0002950300100000078
Further, in the step three, an optimized objective function form is constructed, a matched priori spectrum is introduced, and sparse regularization items of the concentration space distribution of each chromophore are added;
the specific process is as follows:
in the biomedical optical function imaging technology field, the contribution of each chromophore of the tissue under a certain wavelength is distinguished according to the difference of spectral characteristics; the multi-wavelength optical absorption coefficient spectrum can be expressed as:
Figure BDA0002950300100000079
wherein mu is a An optical absorption coefficient at a certain pixel position in the target area under a certain wavelength;
Figure BDA0002950300100000081
the molar extinction coefficient of the mth chromophore at the measured wavelength; c (C) m (r) represents the concentration of the m-th chromophore at the site; n is the total number of pixels of the reconstructed optical absorption coefficient image at each wavelength, where N may be 158×158= 24964, where N may be expressed as a matrix multiplication, that is: mu (mu) a =εC;
μ a The optical absorption coefficient expressed by a certain pixel in the map is a weighted sum of the molar extinction coefficient of each chromophore in the tissue under the current measurement wavelength and the corresponding concentration content value of each chromophore, and an objective function is constructed based on the linear condition;
in the constructed form of the objective function, the fidelity term is divided into a plurality of terms to match the prior spectrum
Figure BDA0002950300100000082
In vivo spectrum to be sought- >
Figure BDA0002950300100000083
As two parts of coefficient matrix, the concentration distribution of each chromophore in tissue is separated into intrinsic chromophore concentration distribution vector +.>
Figure BDA0002950300100000084
Exogenous contrast agent concentration distribution vector +.>
Figure BDA0002950300100000085
And regularizing the term pair +.>
Figure BDA0002950300100000086
Is->
Figure BDA0002950300100000087
Realizing the constraint of space distribution;
wherein the method comprises the steps of
Figure BDA0002950300100000088
Its dimension is L×Q->
Figure BDA0002950300100000089
Q is the total number of exogenous chromophore components; i.e., Q may be set to 1;
above-mentioned
Figure BDA00029503001000000810
The dimension is P×N->
Figure BDA00029503001000000811
Figure BDA00029503001000000812
The dimension is Q x N,
Figure BDA00029503001000000813
the objective function constructed is:
Figure BDA00029503001000000814
in the formula, I TV As the total variation norms, lambda and zeta are regularization parameters; in the course of the calculation process,
Figure BDA00029503001000000815
vector and->
Figure BDA00029503001000000816
The vectors are rearranged into image matrixes respectively to meet the requirement of image total variation operation.
Further, in the fourth step, the proposed iterative solution framework can make the iterative process of minimizing the objective function converge, and avoid the step of selecting the search step in each iteration in the process of optimizing reconstruction;
in particular, proposed iterative updating of the in vivo spectrum of an exogenous contrast agent (such as ICG) in a target biological tissue
Figure BDA00029503001000000817
And the concentration profile of the respective chromophore (+.>
Figure BDA00029503001000000818
And->
Figure BDA00029503001000000819
) The updated formula of (2) is:
iterative solution to the unknown in-vivo spectral signature of the jth exogenous contrast agent, where j e [1, q ]:
Figure BDA0002950300100000091
iterative solution to the j-th exogenous contrast agent concentration profile, where j ε [1, Q ]:
Figure BDA0002950300100000092
Iterative solution to the i-th intrinsic chromophore concentration distribution, where i ε [1, P ]:
Figure BDA0002950300100000093
in the above-mentioned description of the invention,
Figure BDA0002950300100000094
a division operation element by element within two vectors representing the same dimension; />
Figure BDA0002950300100000098
Representing a Hadamard product operator; k isCurrent iteration number; lambda (lambda) k With xi k The regularization parameters selected in the kth iteration are selected, and the convergence of the updating formula is ensured by the selection of the regularization parameters; />
Figure BDA0002950300100000095
A vector representation after the concentration distribution vector of the j-th exogenous contrast agent is rearranged into an image matrix and subjected to total variation partial derivative operation; />
Figure BDA0002950300100000096
The vector representation after the partial derivative operation of the total variation of the concentration distribution vector rearrangement of the ith intrinsic chromophore into an image matrix is shown.
In the fifth step, the iteration number k is updated to k+1, whether the iteration stopping condition is reached is judged, and if the iteration stopping condition is not reached, the fourth step is returned; otherwise, stopping iteration, and outputting the finally reconstructed in-vivo spectrum of the exogenous contrast agent and the concentration distribution image of each chromophore.
Further, the specific judging mode of the iteration stop condition is as follows:
calculating the difference error of the fidelity term in the objective function in the process of updating the adjacent two iterations, setting an iteration stop criterion tol, and stopping the iteration if the error is less than or equal to tol, wherein the difference error is equal to or less than the total of the iteration stop criterion tol
Figure BDA0002950300100000097
Outputting the reconstruction result of the (a);
the tol value is set to 10 -4 The iteration stop condition is reached when the iteration times k are reached, and the synchronous acquisition of unknown in-vivo spectrum of the exogenous contrast agent and the concentration distribution image of each chromophore is finally realized;
in the process of calculating the adjacent two iterative updating, the formula of the difference error of the fidelity term in the objective function is as follows:
Figure BDA0002950300100000101
by combining all the technical schemes, the invention has the advantages and positive effects that: the invention performs an initialization operation to estimate spectral characteristics and concentration distribution of each chromophore in a target biological tissue, comprising: performing blind unmixing pretreatment on the multi-wavelength optical absorption coefficient map by adopting an automatic target generation Algorithm (ATGP), performing spectrum angle matching calculation on a spectrum unmixing result and a spectrum (i.e. a priori spectrum) of a natural chromophore (such as oxyhemoglobin, deoxyhemoglobin, water, lipid, melanin and the like) in a biological tissue, and replacing the corresponding chromophore spectrum analyzed in the pretreatment result by using the priori spectrum so as to improve the accuracy of unknown in-vivo spectrum and concentration calculation of the exogenous contrast agent in a subsequent step; then, constructing an optimized solving objective function, wherein a fidelity term of the optimized solving objective function is divided into two parts taking a matching priori spectrum and a spectrum to be solved as coefficient matrixes, and introducing spatial constraint on the concentration distribution of each chromophore by utilizing a total variation regularization Term (TV); and finally, inverting the objective function by using the proposed iteration solution framework based on non-negative matrix factorization, so that the unknown in-vivo spectrum of the exogenous contrast agent and the concentration distribution of each chromophore are synchronously acquired, the convergence of an iteration process is ensured, and the step of selecting a search step length in each iteration of the traditional optimization reconstruction algorithm is avoided.
Meanwhile, in the practical application of multispectral biomedical optical function imaging, the spectral characteristics of the exogenous contrast agent can be changed in biological tissues, the spectrum of the inherent chromophore (such as oxyhemoglobin, deoxyhemoglobin, water, lipid, melanin and the like) in the biological tissues is known, and the exogenous contrast agent can be used as a priori spectrum and introduced into an initialization step by utilizing a spectrum angle matching mode so as to improve the accuracy of unknown in-vivo spectrum and concentration calculation of the exogenous contrast agent in a subsequent step. The invention constructs an optimized solving objective function, the fidelity term of the optimized solving objective function is divided into two parts taking the matched priori spectrum and the spectrum to be solved as coefficient matrixes, and the spatial constraint on the concentration distribution of each chromophore is introduced by utilizing the total variation regularization term, so that a high-quality reconstruction result can be obtained under the condition that the input multi-wavelength optical absorption spectrum is influenced by noise. The iteration solution framework provided by the invention inverts the objective function, so that the synchronous reconstruction of unknown in-vivo spectrum of the exogenous contrast agent and concentration distribution images of each chromophore can be realized, and meanwhile, the provided iteration update formula not only ensures the convergence of the process of minimizing the objective function, but also avoids the step of selecting the search step length in each iteration of the traditional optimization reconstruction algorithm.
The results of the image reconstruction of the chromophore concentration of the target tissue for the method of the present invention are compared with those not involving an effective sparsity constraint method in table 1. The result shows that the effective space sparse regularization is beneficial to improving the reconstruction quality, and even if noise influence exists in the simulation synthesized data, the method can reconstruct the concentration distribution image with high fidelity.
TABLE 1
Figure BDA0002950300100000111
The results of the reconstruction of unknown spectra of exogenous contrast agents by the method of the present invention and the common blind unmixing algorithms (VCA and ATGP) are compared in fig. 5. The result shows that the spectrum characteristics obtained by the blind unmixing algorithm (such as VCA and ATGP) are quite different from the accurate spectrum characteristics due to the sensitivity to noise interference. In contrast, the spectrum characteristics obtained by reconstruction in the method are consistent with the preset spectrum profile, which shows that in the method step provided by the invention, the spectrum of the inherent chromophore in the biological tissue is used as the priori spectrum to be beneficial to accurately estimating the unknown spectrum characteristics of the exogenous contrast agent.
Drawings
In order to more clearly illustrate the technical solutions of the embodiments of the present application, the following description will briefly explain the drawings needed in the embodiments of the present application, and it is obvious that the drawings described below are only some embodiments of the present application, and that other drawings can be obtained according to these drawings without inventive effort for a person skilled in the art.
Fig. 1 is a schematic flow chart of an optical function imaging method combining spatial regularization and semi-blind spectral unmixing provided by an embodiment of the invention.
Fig. 2 is a schematic flow chart of an implementation of an optical functional imaging method combining spatial regularization and semi-blind spectral unmixing according to an embodiment of the present invention.
FIG. 3 is a graph showing the multi-wavelength optical absorption coefficient spectrum μ according to an embodiment of the present invention a An image schematic of 4 wavelengths.
FIG. 4 is a schematic diagram of the reconstruction results of resolving unknown spectral characteristics of exogenous contrast agents with respective chromophore concentration profile images at a noise level of 20dB provided by an embodiment of the present invention.
In fig. 4: FIGS. (a) - (c) are graphs simulating oxyhemoglobin (HbO) in a tumor microenvironment 2 ) A true concentration profile image of deoxyhemoglobin (Hb) and exogenous contrast agent (indocyanine green ICG);
FIGS. (d) through (f) illustrate the analysis of each chromophore (oxyhemoglobin HbO) by the method of the present invention 2 Deoxyhemoglobin Hb, and exogenous contrast agent ICG).
Fig. 5 is a graph of spectral characteristics of an exogenous contrast agent ICG reconstructed at a noise level of 20dB with a signal-to-noise ratio provided by an embodiment of the present invention.
Fig. 6 is a change trend chart of a difference value of an objective function in two adjacent iteration updating processes when an iteration 278 reaches an iteration stop condition after an iteration solution framework provided by the embodiment of the invention performs iteration solution.
Detailed Description
The present invention will be described in further detail with reference to the following examples in order to make the objects, technical solutions and advantages of the present invention more apparent. It should be understood that the specific embodiments described herein are for purposes of illustration only and are not intended to limit the scope of the invention.
Aiming at the problems existing in the prior art, the invention provides an optical function imaging method combining spatial regularization and semi-blind spectrum unmixing, and the invention is described in detail below with reference to the accompanying drawings.
As shown in fig. 1: the optical function imaging method combining spatial regularization and a semi-blind spectrum unmixing strategy provided by the embodiment of the invention comprises the following steps:
s101: input multi-wavelength optical absorption coefficient map mu a
S102: initializing spectral characteristics and concentration distribution of each chromophore in target biological tissues;
s103: constructing an optimized objective function form, introducing a matched priori spectrum, and adding sparse regularization items of the concentration spatial distribution of each chromophore;
s104: inverting the objective function by using the proposed iterative solution framework, and iteratively updating the in-vivo spectral characteristics of the exogenous contrast agent in the biological tissue and the concentration distribution image of each chromophore;
S105: judging whether an iteration stop condition is reached, and returning to the step S104 if the iteration stop condition is not reached; otherwise, stopping iteration and outputting the spectral characteristics of the exogenous contrast agent and the concentration distribution image of each chromophore.
One of ordinary skill in the art may also implement other steps in the optical function imaging method combining the spatial regularization and the semi-blind spectral unmixing strategy provided by the present invention, and the optical function imaging method combining the spatial regularization and the semi-blind spectral unmixing strategy provided by the present invention in fig. 1 is merely a specific embodiment.
In S101 provided by the embodiment of the invention, inputting reconstructed target tissue region multi-wavelength optical absorption coefficient map mu a The method specifically comprises the following steps:
designing a microenvironment of a blood vessel simulated tumor area, wherein the size of a simulated body image is 158 multiplied by 158 pixels; the experiment is that 35 excitation wavelengths are selected at intervals of 5nm in the wavelength range of 680nm to 884 nm; within this band, a biological tissue intrinsic chromophore (oxyhemoglobin HbO) reflecting a physiological change in a tumor region 2 Deoxyhemoglobin Hb), and exogenous contrast agents (indocyanine green ICG) have unique spectral absorption characteristics; gaussian white noise is superimposed in the simulated multispectral optical absorption profile image in such a way that a data set μ with a signal-to-noise ratio of 20dB is formed a And (5) a map. Mu, optionally 4 wavelengths a An image, as shown in fig. 3.
In the step S102 provided by the embodiment of the invention, blind solution mixing pretreatment is performed by utilizing an ATGP algorithm, and a priori spectrum of a natural chromophore in a biological tissue is introduced, so that the spectral characteristics and the concentration distribution of each chromophore in a target biological tissue are initialized;
the specific process is as follows:
(1) Determining the total number W of chromophore components contained in biological tissue region, and mapping the optical absorption coefficient mu for multiple wavelengths a Blind solution is realized, and preliminary estimation of absorption spectrum characteristics (molar extinction coefficients) of each chromophore in tissue is obtained
Figure BDA0002950300100000131
Wherein the method comprises the steps of
Figure BDA0002950300100000132
The dimension is L×W->
Figure BDA0002950300100000133
L is the total number of wavelengths selected for measurement;
wherein L is 35; since the tumor development and development process has close relation with tumor angiogenesis, hemoglobin is a main intrinsic chromophore in blood, and after fluorescent dye approved by FDA and having definite clinical application is added;
in the near infrared region of the spectrum, comprises three main absorbers (HbO) 2 Hb, exogenous contrast agent ICG), the total number of chromophore components contained within the defined biological tissue region W is 3.
Said obtaining a preliminary estimate of the absorption spectra of these three types of chromophores
Figure BDA0002950300100000134
The method comprises the following steps:
the extracted nth target absorption spectrum t is extracted by utilizing ATGP in nth iteration n Satisfies the following formula:
Figure BDA0002950300100000141
in the method, in the process of the invention,
Figure BDA0002950300100000142
is an orthogonal projection matrix; i is an identity matrix; u (U) n-1 =[t 1 ,t 2 ,…,t n-1 ]A target absorption spectrum matrix generated for the n-1 th calculation; r is the vector of image pixels.
(2) Calculation of absorption spectra of respective chromophores in tissue preliminarily estimated by ATGP
Figure BDA0002950300100000143
A priori spectral features with intrinsic chromophores>
Figure BDA0002950300100000144
Spectral angle error between->
Figure BDA0002950300100000145
Measuring similarity between spectra; let P be the total number of intrinsic chromophore components in the biological tissue +.>
Figure BDA0002950300100000146
The dimension is L×P->
Figure BDA0002950300100000147
HbO 2 And Hb as the intrinsic chromophore in the target biological tissue region, the spectral signature is a priori, so P is 2 in this embodiment.
(3) Setting a spectral angle error threshold angle_tol, if ATGP preliminarily estimates the absorption spectrum of each chromophore in the tissue
Figure BDA0002950300100000148
And a priori spectrum->
Figure BDA0002950300100000149
The spectral angle error between them is smaller than the set threshold, i.e.>
Figure BDA00029503001000001410
Then the corresponding chromophore spectrum resolved in the pretreatment result is replaced by the prior spectrum to improve the follow-upThe method is used for obtaining the initialization spectral characteristics epsilon of each chromophore in the target biological tissue by calculating the accuracy of unknown in-vivo spectrum and concentration of the exogenous contrast agent and ending the step:
the formula for calculating the spectral angle error is shown as follows:
Figure BDA00029503001000001411
(4) Thus far, a multi-wavelength optical absorption coefficient spectrum μ has been obtained a And respective chromophores (HbO) in biological tissue 2 Hb, ICG), the initialized concentration profile C of each chromophore within the biological tissue is solved using a generalized inverse algorithm to complete all initialization steps:
Figure BDA00029503001000001412
in S103 provided by the embodiment of the invention, an optimized objective function form is constructed, a matched priori spectrum is introduced, and sparse regularization items of the concentration space distribution of each chromophore are added;
the specific process is as follows:
in the field of biomedical optical functional imaging technology, the contribution of each chromophore of a tissue at a certain wavelength is distinguished according to the difference of spectral features. The multi-wavelength optical absorption coefficient spectrum can be expressed as:
Figure BDA0002950300100000151
wherein mu is a An optical absorption coefficient at a certain pixel position in the target area under a certain wavelength;
Figure BDA0002950300100000152
the molar extinction coefficient of the mth chromophore at the measured wavelength; c (C) m (r) represents the concentration of the m-th chromophore at the site; n is the optical absorption coefficient graph reconstructed at each wavelengthThe total number of pixels of the image, in this embodiment, N is 158×158= 24964, where N can be expressed as a matrix multiplication, that is: mu (mu) a =εC;
μ a The optical absorption coefficient represented by a pixel in the map is a weighted sum of the molar extinction coefficient of each chromophore in the tissue at the current measurement wavelength and the corresponding concentration content value of each chromophore, and an objective function is constructed based on the linear condition.
In the form of the objective function constructed in S103, the fidelity term is divided to match the prior spectrum
Figure BDA0002950300100000153
(HbO 2 Hb) and the in vivo spectrum to be sought +.>
Figure BDA0002950300100000154
(ICG) as two parts of the coefficient matrix, the concentration profile of each chromophore in the tissue is accordingly separated into the intrinsic chromophore concentration profile vector +.>
Figure BDA0002950300100000155
Exogenous contrast agent concentration distribution vector +.>
Figure BDA0002950300100000156
And regularizing the ++using the total variation regularization Term (TV)>
Figure BDA0002950300100000157
Is->
Figure BDA0002950300100000158
Realizing the constraint of space distribution;
wherein the method comprises the steps of
Figure BDA0002950300100000159
Its dimension is L×Q->
Figure BDA00029503001000001510
Q is the total number of exogenous chromophore components; i.e., Q is 1 in this embodiment;
above-mentioned
Figure BDA00029503001000001511
The dimension is P×N->
Figure BDA00029503001000001512
Figure BDA00029503001000001513
The dimension is Q x N,
Figure BDA00029503001000001514
the objective function constructed in S103 is:
Figure BDA00029503001000001515
in the formula, I TV For the total variation norm, λ and ζ are regularization parameters. In the course of the calculation process,
Figure BDA00029503001000001516
vector and->
Figure BDA00029503001000001517
The vectors are rearranged into image matrixes respectively to meet the requirement of image total variation operation.
In S104 provided by the embodiment of the present invention, the proposed iterative solution framework may enable the iterative process of minimizing the objective function to converge, and avoid the step of selecting the search step in each iteration in the process of optimizing reconstruction.
In particular, proposed iterative updating of the in vivo spectrum of an exogenous contrast agent (ICG) in a target biological tissue
Figure BDA00029503001000001611
And the concentration profile of the respective chromophore (+. >
Figure BDA0002950300100000161
And->
Figure BDA0002950300100000162
) The updated formula of (2) is:
iterative solution to spectral features of a j-th exogenous contrast agent, where j ε [1, Q ]:
Figure BDA0002950300100000163
/>
iterative solution to the j-th exogenous contrast agent concentration profile, where j ε [1, Q ]:
Figure BDA0002950300100000164
iterative solution to the i-th intrinsic chromophore concentration distribution, where i ε [1, P ]:
Figure BDA0002950300100000165
in the above-mentioned description of the invention,
Figure BDA0002950300100000166
a division operation element by element within two vectors representing the same dimension; />
Figure BDA00029503001000001610
Representing a Hadamard product operator; k is the current iteration number; lambda (lambda) k With xi k Regularization parameters selected in the kth iteration are selected; />
Figure BDA0002950300100000167
A vector representation after the concentration distribution vector of the j-th exogenous contrast agent is rearranged into an image matrix and subjected to total variation partial derivative operation; />
Figure BDA0002950300100000168
Representing that the concentration distribution vector of the ith intrinsic chromophore is rearranged into an image matrix and subjected to total variation partial derivative operationIs a vector representation of (c).
In S105 provided by the embodiment of the present invention, the iteration number k is updated to k+1, and whether an iteration stop condition is reached is determined, if not, the step returns to S104; otherwise, stopping iteration, and outputting the finally reconstructed in-vivo spectrum of the exogenous contrast agent and the concentration distribution image of each chromophore.
The specific judging mode of the iteration stop condition is as follows:
Calculating the difference error of the fidelity term in the objective function in the process of updating the adjacent two iterations, setting an iteration stop criterion tol, and stopping the iteration if the error is less than or equal to tol, wherein the difference error is equal to or less than the total of the iteration stop criterion tol
Figure BDA0002950300100000169
And (3) outputting the reconstruction result of the model. In the present embodiment, the tol value is set to 10 -4 The iteration stop condition is reached when the iteration number k is 278, and finally the unknown in-vivo spectrum of the exogenous contrast agent (ICG) and the chromophore (HbO) are realized 2 Hb, ICG) concentration distribution images.
In the process of calculating the adjacent two iterative updating, the formula of the difference error of the fidelity term in the objective function is as follows:
Figure BDA0002950300100000171
the invention constructs an optimized objective function form, introduces the spectrum of the inherent chromophore in the biological tissue as a priori spectrum and adds spatial regularization constraint in image reconstruction. In order to quantitatively evaluate the superiority of spatial regularization constraints in reconstructing a concentration profile image, image reconstruction results of the inventive method (Priori-NMF-TV) and those not involving the effective sparsity constraint method (Priori-NMF) were compared, and the results are shown in table 1. The parameter distance (d) is selected as a quantization index to represent the difference between the reconstructed image and the real image. The definition of the parameter d is:
Figure BDA0002950300100000172
where u and v represent the reconstructed image and the real image, respectively. The image size is mxn. The smaller the d value, the smaller the difference between the reconstructed image and the real image.
The invention can also use a hand-held or fixed-mounted biomedical optical detection device (such as an optical endoscope, a diffuse fluorescence tomography system, a diffuse optical breast tomography system, a multi-wavelength photoacoustic tomography system, a tissue biopsy microscope and the like) to excite and image in vivo tissues or tissue samples at a certain wavelength interval by utilizing a light source in a near infrared first window wavelength range (600 nm-900 nm) or a near infrared second window wavelength range (1000 nm-1700 nm), wherein the image size can be set in the detection device, for example, 35 excitation wavelengths are selected at a step length interval of 5nm in a wavelength range of 680 nm-884 nm, so as to obtain an image with a size of 158 multiplied by 158 pixels at each wavelength; within this band, the biological tissue intrinsic chromophore reflecting the physiological change of the tumor region has unique spectral absorption characteristics; thereby obtaining the multispectral optical absorption coefficient mu of the target tissue region a And (5) a map. After that, the software written based on the algorithm of the invention is operated in the biomedical optical detection system, so that the in-vivo concentration distribution image of the endogenous chromophore (hemoglobin) and the exogenous contrast agent (ICG) with biological markers which clinically reflect the physiological and pathological states of the tissues can be further obtained, the blood oxygen saturation value of the tissues can be further calculated, and important in-vivo physiological parameters can be provided for clinical diagnosis.
The technical effects of the present invention will be described in detail with reference to simulation experiments.
As shown in fig. 4-6, the experimental result of the invention shows that in the optical function imaging method combining spatial regularization and semi-blind spectrum unmixing provided by the invention, an optimized objective function is constructed, namely, the spectrum of the inherent chromophore in the biological tissue is considered as a priori spectrum, and spatial regularization constraint is added in image reconstruction, so that under the condition that the multi-wavelength optical absorption coefficient spectrum is influenced by noise, unknown in-vivo spectral characteristics of the exogenous contrast agent can be effectively estimated, the concentration distribution image of each chromophore can be reconstructed with high quality, and synchronous calculation of the unknown in-vivo spectral characteristics of the exogenous contrast agent can be realized. The results of the reconstruction of unknown spectra of exogenous contrast agents by the method of the present invention and the common blind unmixing algorithms (VCA and ATGP) are compared in fig. 5. The result shows that the spectrum characteristics obtained by the blind unmixing algorithm (such as VCA and ATGP) are quite different from the accurate spectrum characteristics due to the sensitivity to noise interference. In contrast, the spectrum characteristics obtained by reconstruction in the method are consistent with the preset spectrum profile, which shows that in the method step provided by the invention, the spectrum of the inherent chromophore in the biological tissue is used as the priori spectrum to be beneficial to accurately estimating the unknown spectrum characteristics of the exogenous contrast agent. In the iteration process, the iteration solution framework provided by the method avoids the selection of the search step length in each iteration, and ensures that the iteration process of the minimized objective function converges under the self-adaptive multiplication update.
It should be noted that the embodiments of the present invention can be realized in hardware, software, or a combination of software and hardware. The hardware portion may be implemented using dedicated logic; the software portions may be stored in a memory and executed by a suitable instruction execution system, such as a microprocessor or special purpose design hardware. Those of ordinary skill in the art will appreciate that the apparatus and methods described above may be implemented using computer executable instructions and/or embodied in processor control code, such as provided on a carrier medium such as a magnetic disk, CD or DVD-ROM, a programmable memory such as read only memory (firmware), or a data carrier such as an optical or electronic signal carrier. The device of the present invention and its modules may be implemented by hardware circuitry, such as very large scale integrated circuits or gate arrays, semiconductors such as logic chips, transistors, etc., or programmable hardware devices such as field programmable gate arrays, programmable logic devices, etc., as well as software executed by various types of processors, or by a combination of the above hardware circuitry and software, such as firmware.
The foregoing is merely illustrative of specific embodiments of the present invention, and the scope of the invention is not limited thereto, but any modifications, equivalents, improvements and alternatives falling within the spirit and principles of the present invention will be apparent to those skilled in the art within the scope of the present invention.

Claims (7)

1. An optical function imaging method combining spatial regularization and semi-blind spectral unmixing, characterized in that the optical function imaging method combining spatial regularization and semi-blind spectral unmixing comprises:
input multi-wavelength optical absorption coefficient map mu a
Initializing spectral characteristics and concentration distribution of each chromophore in target biological tissues;
constructing an optimized objective function form, introducing a matched priori spectrum, and adding sparse regularization items of the concentration spatial distribution of each chromophore;
inverting the objective function by using the proposed iterative solution framework, and iteratively updating the in-vivo spectral characteristics of the exogenous contrast agent in the biological tissue and the concentration distribution image of each chromophore;
judging whether an iteration stop condition is reached, if not, returning to perform inversion on the objective function by using the proposed iteration solution framework; otherwise, stopping iteration, and outputting spectral characteristics of the exogenous contrast agent and concentration distribution images of the chromophores;
constructing an optimized objective function form, introducing a matched priori spectrum, and adding sparse regularization items of the concentration spatial distribution of each chromophore; the specific process is as follows:
in the biomedical optical function imaging technology field, the contribution of each chromophore of the tissue under a certain wavelength is distinguished according to the difference of spectral characteristics; the multi-wavelength optical absorption coefficient spectrum is expressed as:
Figure FDA0004229497130000011
Wherein mu is a Is a multi-wavelength optical absorption coefficient map under a certain wavelength;
Figure FDA0004229497130000012
to measure the molar extinction coefficient of the mth chromophore at the wavelength; c (C) m (r) represents the concentration of the m-th chromophore at the site; n is the total number of pixels of the reconstructed optical absorption coefficient image at each wavelength, where N is 158×158= 24964, and expressed as a matrix multiplication, that is: mu (mu) a =εC;
Multi-wavelength optical absorption coefficient spectrum mu a The optical absorption coefficient represented by a certain pixel is a weighted sum of the molar extinction coefficient of each chromophore in the tissue under the current measurement wavelength and the corresponding concentration content value of each chromophore, and an objective function is constructed based on the linear condition;
in the constructed form of the objective function, the fidelity term is divided into a priori spectra
Figure FDA00042294971300000213
In vivo spectrum to be sought->
Figure FDA00042294971300000212
As two parts of coefficient matrix, the concentration distribution of each chromophore in tissue is separated into intrinsic chromophore concentration distribution vector +.>
Figure FDA0004229497130000021
Exogenous contrast agent concentration distribution vector +.>
Figure FDA00042294971300000216
And regularizing the term pair +.>
Figure FDA00042294971300000214
Is->
Figure FDA00042294971300000215
Realizing the constraint of space distribution;
wherein the method comprises the steps of
Figure FDA0004229497130000022
Its dimension is L×Q->
Figure FDA0004229497130000023
Q is the total number of exogenous chromophore components; i.e. Q is set to 1;
above-mentioned
Figure FDA0004229497130000024
The dimension is P×N->
Figure FDA0004229497130000025
The dimension is Q x N,
Figure FDA0004229497130000026
The objective function constructed is:
Figure FDA0004229497130000027
in the formula, I TV As the total variation norms, lambda and zeta are regularization parameters; in the course of the calculation process,
Figure FDA0004229497130000028
vector and->
Figure FDA0004229497130000029
The vectors are rearranged into image matrixes respectively to meet the requirement of image total variation operation; initializing concentration distribution C, initializing spectral characteristics epsilon and L to be the total number of wavelengths selected for measurement and the total number of inherent chromophore components in biological tissues P;
the proposed iterative solution framework enables the iterative process of minimizing the objective function to converge, and avoids the step of selecting the search step length in each iteration in the process of optimizing and reconstructing;
in particular, proposed iterative updating of the in vivo spectrum of an exogenous contrast agent in a target biological tissue
Figure FDA00042294971300000217
And the concentration profile of the respective chromophore, +.>
Figure FDA00042294971300000218
And->
Figure FDA00042294971300000219
The updated formula of (2) is:
iterative solution to the unknown in-vivo spectral signature of the jth exogenous contrast agent, where j e [1, q ]:
Figure FDA00042294971300000210
iterative solution to the j-th exogenous contrast agent concentration profile, where j ε [1, Q ]:
Figure FDA00042294971300000211
iterative solution to the i-th intrinsic chromophore concentration distribution, where i ε [1, P ]:
Figure FDA0004229497130000031
in the above-mentioned description of the invention,
Figure FDA0004229497130000032
a division operation element by element within two vectors representing the same dimension; />
Figure FDA0004229497130000035
Representing a Hadamard product operator; k is the current iteration number; lambda (lambda) k With xi k The regularization parameters selected in the kth iteration are selected, and the convergence of the updating formula is ensured by the selection of the regularization parameters; />
Figure FDA0004229497130000033
A vector representation after the concentration distribution vector of the j-th exogenous contrast agent is rearranged into an image matrix and subjected to total variation partial derivative operation;
Figure FDA0004229497130000034
a vector representation after partial derivative operation of rearranging the concentration distribution vector of the ith natural chromophore into an image matrix and performing total variation on the image matrix;
inputting the multi-wavelength optical absorption coefficient spectrum mu of the target tissue region obtained by direct acquisition or reconstruction of biomedical optical detection equipment a The method specifically comprises the following steps:
the optical detection device uses a light source with a near infrared first window wavelength range of 600 nm-900 nm or a near infrared second window wavelength range of 1000 nm-1700 nm in biological tissue detection to perform multi-wavelength excitation and imaging on pathological tissues at a certain wavelength interval, the image size is set in the detection device, 35 excitation wavelengths are selected at a step length interval of 5nm in a wavelength range of 680 nm-884 nm, and an image with a size of 158×158 pixels under each wavelength is obtained; within this band, the biological tissue intrinsic chromophore reflecting the physiological change of the tumor region has unique spectral absorption characteristics; thereby obtaining a multi-wavelength optical absorption coefficient spectrum mu of the target tissue region a The method comprises the steps of carrying out a first treatment on the surface of the The optical detection device is an optical endoscope, an optical mammary gland imager or a tissue biopsy microscope.
2. The optical function imaging method combining spatial regularization and semi-blind spectral unmixing according to claim 1, wherein the blind unmixing pretreatment is performed by using an ATGP algorithm, and a priori spectrum of inherent chromophores in biological tissues is introduced, so that the spectral characteristics and concentration distribution of each chromophore in target biological tissues are initialized.
3. The optical function imaging method combining spatial regularization and semi-blind spectral unmixing according to claim 1, wherein the specific process of initializing the spectral characteristics of each chromophore and its concentration distribution in the target biological tissue is as follows:
determining the total number W of chromophore components contained in biological tissue region, and mapping the optical absorption coefficient mu for multiple wavelengths a Realizing blind solution and obtaining the absorption spectrum of each chromophore in the tissue estimated preliminarily
Figure FDA0004229497130000041
Wherein->
Figure FDA0004229497130000042
The dimension is L×W->
Figure FDA0004229497130000043
L is the total number of wavelengths selected for measurement;
wherein L is set to 35; since the tumor development and development process has close relation with tumor angiogenesis, hemoglobin is a main intrinsic chromophore in blood, and after fluorescent dye approved by FDA and having definite clinical application is added;
In the near infrared region of the spectrum, three main absorbers are involved: hbO (HbO) 2 Hb, and an exogenous contrast agent ICG, so that the total number of chromophore components contained in the biological tissue region, W, is determined to be 3;
calculation of absorption spectra of respective chromophores in tissue preliminarily estimated by ATGP
Figure FDA0004229497130000044
And a priori spectrum->
Figure FDA0004229497130000045
Spectral angular error between
Figure FDA0004229497130000046
Measuring similarity between spectra; assuming P is the total number of intrinsic chromophore components in the biological tissue
Figure FDA0004229497130000047
The dimension is L×P->
Figure FDA0004229497130000048
HbO 2 The Hb is taken as a natural chromophore in a target biological tissue area, the spectrum characteristic of the Hb is taken as a priori spectrum, and therefore P is 2;
setting a spectral angle error threshold angle_tol, if ATGP preliminarily estimates the absorption spectrum of each chromophore in the tissue
Figure FDA0004229497130000049
And a priori spectrum->
Figure FDA00042294971300000410
The spectral angle error between them is smaller than the set threshold, i.e.>
Figure FDA00042294971300000411
The prior spectrum is used for replacing the spectrum of the corresponding chromophore analyzed in the result so as to improve the accuracy of the follow-up step on the calculation of the unknown in-vivo spectrum and the concentration of the exogenous contrast agent, and the step is finished to obtain the initialized spectral characteristics epsilon of each chromophore in the target biological tissue;
the obtained multi-wavelength optical absorption coefficient spectrum mu a And chromophores, hbO, in biological tissue 2 Initializing spectral characteristics epsilon of Hb and ICG, and solving an initialized concentration distribution C of each chromophore in biological tissues by using a generalized inverse algorithm to complete all initialization steps:
C=ε -1 μ a
4. The optical functional imaging method combining spatial regularization and semi-blind spectral unmixing of claim 3, wherein said obtaining a preliminary estimate of in-tissue HbO 2 Absorption spectrum of three chromophores of Hb and ICG
Figure FDA00042294971300000412
MethodThe following are provided:
the extracted nth target absorption spectrum t is extracted by utilizing ATGP in nth iteration n Satisfies the following formula:
Figure FDA0004229497130000051
in the method, in the process of the invention,
Figure FDA0004229497130000052
is an orthogonal projection matrix; i is an identity matrix; u (U) n-1 =[t 1 ,t 2 ,…,t n-1 ]A target absorption spectrum matrix generated for the n-1 th calculation; r is the vector of image pixels.
5. An optical functional imaging method combining spatial regularization and semi-blind spectral unmixing as recited in claim 3, wherein the calculated spectral angle error formula is represented by:
Figure FDA0004229497130000053
6. the optical function imaging method combining spatial regularization and semi-blind spectral unmixing according to claim 1, characterized in that the iteration number k is updated to k+1, whether an iteration stop condition is reached is judged, if not, the objective function is inverted by the proposed iteration solution frame, and the in-vivo spectral characteristics of the exogenous contrast agent in the biological tissue and the concentration distribution image of each chromophore are iteratively updated; otherwise, stopping iteration, and outputting the finally reconstructed in-vivo spectrum of the exogenous contrast agent and the concentration distribution image of each chromophore;
The specific judging mode of the iteration stop condition is as follows:
calculating the difference error of the fidelity term in the objective function in the process of updating the adjacent two iterations, setting an iteration stop criterion tol, and stopping the iteration if the error is less than or equal to tol, wherein the difference error is equal to or less than the total of the iteration stop criterion tol
Figure FDA0004229497130000054
Outputting the reconstruction result of the (a);
the tol value is set to 10 -4 When the iteration times k are the same, iteration stop conditions are reached, and finally, the synchronous acquisition of unknown in-vivo spectrums of the exogenous contrast agent and the concentration distribution images of each chromophore is realized;
in the process of calculating the adjacent two iterative updating, the formula of the difference error of the fidelity term in the objective function is as follows:
Figure FDA0004229497130000055
7. a biomedical optical function imaging terminal for implementing the optical function imaging method combining spatial regularization and semi-blind spectral unmixing of any one of claims 1-6.
CN202110205426.0A 2021-02-24 2021-02-24 Optical function imaging method combining spatial regularization and semi-blind spectrum unmixing Active CN113066142B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110205426.0A CN113066142B (en) 2021-02-24 2021-02-24 Optical function imaging method combining spatial regularization and semi-blind spectrum unmixing

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110205426.0A CN113066142B (en) 2021-02-24 2021-02-24 Optical function imaging method combining spatial regularization and semi-blind spectrum unmixing

Publications (2)

Publication Number Publication Date
CN113066142A CN113066142A (en) 2021-07-02
CN113066142B true CN113066142B (en) 2023-07-07

Family

ID=76558948

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110205426.0A Active CN113066142B (en) 2021-02-24 2021-02-24 Optical function imaging method combining spatial regularization and semi-blind spectrum unmixing

Country Status (1)

Country Link
CN (1) CN113066142B (en)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110596017B (en) * 2019-09-12 2022-03-08 生态环境部南京环境科学研究所 Hyperspectral image soil heavy metal concentration assessment method based on space weight constraint and variational self-coding feature extraction
CN113793646B (en) * 2021-09-29 2023-11-28 上海交通大学 Spectral image unmixing method based on weighted non-negative matrix factorization and application thereof

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101697008A (en) * 2009-10-20 2010-04-21 北京航空航天大学 Hyperspectral unmixing method for estimating regularized parameter automatically
CN102314685A (en) * 2011-07-23 2012-01-11 北京航空航天大学 Hyperspectral image sparse unmixing method based on random projection
CN104027132A (en) * 2014-06-09 2014-09-10 苏州大学 Device and method based on multispectral photoacoustic tomography
CN104952050A (en) * 2015-07-07 2015-09-30 西安电子科技大学 Self-adaptive hyperspectral image unmixing method based on region segmentation

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11510600B2 (en) * 2012-01-04 2022-11-29 The Trustees Of Dartmouth College Method and apparatus for quantitative and depth resolved hyperspectral fluorescence and reflectance imaging for surgical guidance
US10438258B2 (en) * 2016-03-21 2019-10-08 The Procter & Gamble Company Method and apparatus for generating graphical chromophore maps

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101697008A (en) * 2009-10-20 2010-04-21 北京航空航天大学 Hyperspectral unmixing method for estimating regularized parameter automatically
CN102314685A (en) * 2011-07-23 2012-01-11 北京航空航天大学 Hyperspectral image sparse unmixing method based on random projection
CN104027132A (en) * 2014-06-09 2014-09-10 苏州大学 Device and method based on multispectral photoacoustic tomography
CN104952050A (en) * 2015-07-07 2015-09-30 西安电子科技大学 Self-adaptive hyperspectral image unmixing method based on region segmentation

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
A regularization modification to linear spectral unmixing algorithm;Ye Zhang et al;《2012 IEEE International Geoscience and Remote Sensing Symposium》;20121231;全文 *
基于深度非负矩阵分解与聚类空间处理的高光谱解混研究;祝伟;《中国优秀硕士学位论文全文数据库电子期刊 基础科学辑》;20210115;第2021年卷(第1期);全文 *
高光谱遥感图像光谱解混方法研究及其应用;甘玉泉;《中国博士学位论文全文数据库电子期刊 工程科技II辑》;20190615;第2019年卷(第6期);全文 *

Also Published As

Publication number Publication date
CN113066142A (en) 2021-07-02

Similar Documents

Publication Publication Date Title
Wang et al. 3D conditional generative adversarial networks for high-quality PET image estimation at low dose
Guven et al. Diffuse optical tomography with a priori anatomical information
Delannoy et al. SegSRGAN: Super-resolution and segmentation using generative adversarial networks—Application to neonatal brain MRI
Bench et al. Toward accurate quantitative photoacoustic imaging: learning vascular blood oxygen saturation in three dimensions
CN113066142B (en) Optical function imaging method combining spatial regularization and semi-blind spectrum unmixing
Hyde et al. Data specific spatially varying regularization for multimodal fluorescence molecular tomography
CN108903942B (en) Method for identifying spatial difference by using complex fMRI spatial source phase
Schellenberg et al. Semantic segmentation of multispectral photoacoustic images using deep learning
Del Amor et al. Automatic segmentation of epidermis and hair follicles in optical coherence tomography images of normal skin by convolutional neural networks
Dehner et al. Deep-learning-based electrical noise removal enables high spectral optoacoustic contrast in deep tissue
Yamashita et al. Multi-subject and multi-task experimental validation of the hierarchical Bayesian diffuse optical tomography algorithm
Yin et al. Exploring the complementarity of THz pulse imaging and DCE-MRIs: Toward a unified multi-channel classification and a deep learning framework
CN115349137A (en) Correction of flow projection artifacts in OCTA volumes using neural networks
Combalia et al. Digitally stained confocal microscopy through deep learning
Ram et al. Three-dimensional segmentation of the ex-vivo anterior lamina cribrosa from second-harmonic imaging microscopy
Attyé et al. TractLearn: A geodesic learning framework for quantitative analysis of brain bundles
Wu et al. System-level optimization in spectroscopic photoacoustic imaging of prostate cancer
Aja-Fernández et al. Validation of deep learning techniques for quality augmentation in diffusion MRI for clinical studies
Lutzweiler et al. High-throughput sparsity-based inversion scheme for optoacoustic tomography
Lin et al. Rapid measurement of epidermal thickness in OCT images of skin
Geremia et al. Classification forests for semantic segmentation of brain lesions in multi-channel MRI
Edjlali et al. Lq− Lp optimization for multigrid fluorescence tomography of small animals using simplified spherical harmonics
D’Acunto et al. 3D image reconstruction using Radon transform
O’Kelly et al. Evaluating online filtering algorithms to enhance dynamic multispectral optoacoustic tomography
Jüstel et al. Spotlight on nerves: portable multispectral optoacoustic imaging of peripheral nerve vascularization and morphology

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant