CN113066142A - 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
CN113066142A
CN113066142A CN202110205426.0A CN202110205426A CN113066142A CN 113066142 A CN113066142 A CN 113066142A CN 202110205426 A CN202110205426 A CN 202110205426A CN 113066142 A CN113066142 A CN 113066142A
Authority
CN
China
Prior art keywords
chromophore
spectrum
iteration
spectral
optical
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.)
Granted
Application number
CN202110205426.0A
Other languages
Chinese (zh)
Other versions
CN113066142B (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

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/0059—Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence
    • A61B5/0075—Measuring 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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/424—Iterative

Landscapes

  • Health & Medical Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Pathology (AREA)
  • Medical Informatics (AREA)
  • Theoretical Computer Science (AREA)
  • Biophysics (AREA)
  • General Physics & Mathematics (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)

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 spectrum unmixing, which is used for inputting a multi-wavelength optical absorption coefficient mapμa(ii) a Initializing the spectral characteristics and concentration distribution of each chromophore in the target biological tissue; constructing an optimized objective function form, introducing a matched prior spectrum, and adding a sparse regularization term of the concentration spatial distribution of each chromophore; the proposed iterative solution framework is utilized to invert the objective function, and unknown in-vivo spectral characteristics of the external contrast agent in the biological tissue and the concentration distribution image of each chromophore are updated iteratively; judging whether an iteration stop condition is reached; otherwise, stopping iteration, and simultaneously outputting the in-vivo spectrum of the exogenous contrast agent and the concentration distribution image of each chromophore. The method can avoid the selection of the step length of the search direction in the traditional gradient descent iteration process and ensure the iterationThe convergence of the generation 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 spectrum unmixing.
Background
At present, cancer (malignant tumor) seriously threatens human health, and the development of the cancer is often accompanied with the change of physicochemical properties of tissues, including tissue structure, tissue components, blood oxygen content and the like. The pathological change can cause the change of the optical property of the tissue, so the optical function imaging technology can visualize the distribution of the physiological parameters in a non-invasive way, further discover the early pathological changes of the tissue, and achieve the aims of reducing the death rate, improving the prognosis and improving the life quality. The multispectral biomedical optical function imaging technology can quantitatively acquire the concentration distribution (namely functional parameters directly related to pathophysiology) of each chromophore in biological tissues by reconstructing a series of optical absorption distribution maps reflecting the physiological characteristics of the tissues and utilizing a multispectral unmixing algorithm, thereby realizing the analysis of molecular mechanisms of the development of disease formation and providing an effective noninvasive detection and monitoring means for early diagnosis and treatment of the disease. Therefore, the multispectral biomedical optical function imaging technology plays an important role in the field of biomedical research.
The concentration distribution of each chromophore in the biological tissue is analyzed by utilizing a biomedical multispectral optical function imaging technology (such as a multispectral diffusion optical tomography technology and a multispectral photoacoustic tomography technology), namely the spectrum unmixing problem is solved. Under a certain measuring wavelength, based on the linear condition that the optical absorption coefficient represented by a certain pixel in the reconstructed image is the weighted sum of the molar extinction coefficient of each chromophore in the tissue body at the current measuring wavelength and the corresponding concentration content value of each chromophore, 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 is reflected qualitatively or quantitatively, and the physiological processes of the organism, such as metabolism, blood flow and the like are further reflected. The typical method for solving the multi-spectral component decomposition mixing problem is a spectral fitting algorithm, namely, a multi-spectral optical absorption distribution image obtained by pixel-by-pixel processing by using known spectral characteristic information under a least square rule framework is used for imaging the blood oxygen and exogenous contrast agent concentration distribution of a tissue vascular system by solving the least square problem. However, hemoglobin, lipids, melanin and water are the main endogenous chromophores in tissues, whose absorption spectra are known; some exogenous contrast agents, such as Indocyanine green (ICG), which is the only fluorescent dye approved by the U.S. Food and Drug Administration (FDA) for human body, have been widely used for detecting various physiological states of human body for a long time, and their absorption spectra are related to the nature of solvent medium and dye concentration; binding to plasma proteins results in a shift of the main peak in the absorption spectrum by about 25nm at longer wavelengths; furthermore, ICG solutions do not follow the lambert-beer law at plasma concentrations higher than 15mg/l due to the formation of aggregates with increasing concentration, and at higher concentrations the spectral profile shows two absorption peaks, at λ 730nm and λ 805nm, respectively; at a dilution of 5mg/l, the molar extinction coefficient at λ 805nm increased significantly, while the maximum at λ 730nm became shoulder-like. The above problems make it impossible to accurately obtain the in-vivo spectral characteristics of the exogenous contrast agent, so that the spectral fitting algorithm has limited unmixing performance. In addition, the traditional gradient descent method for optimally solving the least square problem often depends on an empirical method for selecting the search step length in the iterative updating process, and the convergence and the non-negativity of the analysis result cannot be ensured.
Because the limitation of the spectrum fitting algorithm is mainly caused by the fact that the in-vivo spectral characteristics of the exogenous contrast agent cannot be accurately obtained, blind unmixing can be used as another method for analyzing the concentration distribution of the internal and external source chromophore of the tissue, a multispectral data set is analyzed based on data characteristics, and exact prior information related to absorption spectrum is not needed. Among the common blind unmixing algorithms, there are Independent Component Analysis (ICA), Vertex Component Analysis (VCA), and the like. ICA uses the statistical properties of chromophore distribution to estimate the concentrations of components, and its applicable conditions require that the distribution of chromophores within a tissue be statistically independent from each other. The VCA is based on the concept that the observation data form a convex simplex in a characteristic space, the absorption spectrum characteristics of each chromophore are searched through orthogonal projection, and then the concentration distribution of the chromophore is obtained, the calculation time of the VCA is greatly reduced compared with ICA, and the VCA has important advantages when a large amount of data are processed or real-time unmixing is carried out; however, for a multi-spectrum optical reconstruction image with large noise, the deviation between the absorption spectrum acquired by the VCA algorithm and the true value is large, and the effectiveness is reduced; meanwhile, the initial projection vector of the growth of the simplex selected in each iteration of the algorithm is random, so that inconsistency of the operation results of multiple times of the VCA algorithm may occur. In contrast, object detection techniques based on orthogonal projection: an Automatic Target Generation Process (ATGP), which is essentially the same as the VCA, selects a certain initial projection vector to avoid randomness, and has stronger stability compared with the VCA; moreover, the ATGP algorithm is not constrained by dimensionality reduction transformation, the integrity of data can be ensured, and the characteristic information of an absorption spectrum is prevented from being damaged. In addition, statistical sub-pixel detection methods, such as Adaptive Matched Filters (AMFs), are proposed for multi-spectral optical functional imaging unmixing, which focuses on identifying a spectral target that is significantly different from the background, compared to typical linear unmixing methods that attempt to resolve all the different spectral components present in the image; however, this method has good performance only in cases where the target probes are sparsely distributed. Based on the above analysis, in order to realize the synchronous reconstruction of unknown in-vivo spectral characteristics of exogenous contrast agents and the concentration distribution of each chromophore in biological tissues, the method provided by the invention is required to have an optimized objective function and an efficient iterative solution framework.
Through the above analysis, the problems and defects of the prior art are as follows:
(1) in the prior art, the traditional spectrum blind unmixing algorithm is susceptible to noise; the spectral fitting algorithm is limited in its use in cases where the 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 negative concentration values with no physical significance in the results, and an ideal chromophore concentration distribution image cannot be obtained.
(3) When an objective function of a spectral unmixing problem is solved, the traditional gradient descent method for optimizing and solving the least square problem often depends on an empirical method for selecting a search step length in an iterative updating process, and the convergence and the nonnegativity of an analytic result cannot be ensured.
The difficulty in solving the above problems and defects is:
(1) under the condition that the spectrum of the exogenous contrast agent is unknown, synchronous and accurate reconstruction of the exogenous chromophore concentration distribution images in the target tissue area and the spectrum characteristics of the exogenous contrast agent is difficult to realize.
(2) In the case where the spectral characteristics of each chromophore are known or have been effectively estimated, in order to reconstruct the concentration distribution image of each chromophore with high fidelity, it is generally necessary to previously acquire a multi-wavelength optical absorption coefficient map μaEnhancement processing is performed to reduce image noise or eliminate reconstruction artifacts caused by insufficient sampling, but in the prior art, these operations belong to preprocessing steps, which increase the time and computational cost of multispectral optical functional imaging.
(3) When an objective function of a spectral unmixing problem is solved, in order to avoid selecting a search step length of iterative update in an empirical mode in an optimization calculation process, an algorithm for realizing the self-adaptive update of the search step length is invented under the condition of optimizing the regularization constraint of a concentration image, the solving process has convergence and the result is nonnegative, and the algorithm design has difficulty.
The significance of solving the problems and the defects is as follows:
(1) the invention can improve the reconstruction signal-to-noise ratio and fidelity of each internal and external source chromophore concentration image and the exogenous contrast agent spectrum by constructing an optimized solving objective function, introducing the spectrum of the inherent chromophore in the biological tissue as a prior spectrum in a multispectral unmixing reconstruction algorithm and adding space regularization constraint in image reconstruction. Avoiding pre-acquisition of multi-wavelength optical absorption systemsNumber map muaThe enhancement processing of the method improves the algorithm efficiency, is expected to realize the near-real-time concentration image calculation, and has significant significance in clinical application.
(2) The invention constructs an efficient iterative solution framework, provides a search step self-adaptive updating calculation formula under the condition of optimizing the regularization constraint of the concentration image, and can realize the synchronization and non-weight building of unknown in-vivo spectral characteristics of the exogenous contrast agent and concentration images of various chromophores. The step of selecting the search step length in each iteration of the traditional optimization reconstruction algorithm is avoided while the convergence of the process of minimizing the objective function is ensured. The invention meets the requirement of high robustness of the used algorithm tool in clinical application.
Disclosure of Invention
Aiming at the problems in the prior art, the invention provides an optical function imaging method combining spatial regularization and semi-blind spectrum unmixing. The spectrum of the inherent chromophore in the biological tissue is introduced as a prior spectrum, space regularization constraint is added in image reconstruction, an effective matrix decomposition reconstruction framework is developed to carry out iterative solution on an objective function in consideration of the condition that the spectrum and the concentration of each chromophore are non-negative, and therefore the unknown in-vivo spectral characteristics of the exogenous contrast agent in the biological tissue and the concentration distribution of each chromophore are synchronously obtained. 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 nonnegativity of the analysis result.
The invention is realized in such a way that an optical function imaging method combining spatial regularization and semi-blind spectral unmixing comprises the following steps:
step one, inputting a multi-wavelength optical absorption coefficient map mua(ii) a The multi-wavelength optical absorption spectrum of the biological tissue input in the step one can be directly obtained or reconstructed by using biomedical optical detection equipment (such as an endoscope, an optical breast imager, a tissue biopsy microscope and the like).
Initializing the spectral characteristics and concentration distribution of each chromophore in the target biological tissue; selecting an ATGP algorithm to preliminarily estimate the absorption spectrum of each chromophore in the tissue in the initialization process, and quantitatively comparing the absorption spectra, wherein the ATGP algorithm is used as a preprocessing tool, so that the ATGP algorithm has stronger robustness compared with other VCA algorithms commonly used in multispectral imaging; moreover, the ATGP algorithm is not restricted by dimensionality reduction transformation, the integrity of data 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 prior spectrum of the endogenous chromophore (hemoglobin, water, lipid, melanin and the like) in the second step, so that the accuracy of calculating the unknown in-vivo spectrum and concentration image of the exogenous contrast agent in the subsequent step is improved;
constructing an optimized objective function form, introducing the matched prior spectrum, and adding a sparse regularization term of the spatial distribution of the concentration of each chromophore; introducing an accurately known prior spectrum, adding space constraint on the concentration distribution of each chromophore by utilizing total variation regularization, and enabling the input multi-wavelength optical absorption spectrum mu to be inaHigh quality reconstruction results are still obtained under the influence of noise or reconstruction artifacts.
Step four, the proposed iterative solution framework is utilized to invert the objective function, and the spectral characteristics of the external contrast agent in the biological tissue and the concentration distribution image of each chromophore are iteratively updated; and step four, synchronous reconstruction of unknown spectral characteristics of the exogenous contrast agent and concentration distribution of each chromophore can be realized, and meanwhile, the iterative updating formula provided by step four not only ensures the convergence of the process of minimizing the target function, but also avoids the step of selecting a search step length in each iteration of the traditional optimization reconstruction algorithm.
Step five, judging whether an iteration stopping condition is reached, and returning to the step four if the iteration stopping condition is not reached; otherwise, stopping the iteration and outputting the spectral characteristics of the exogenous contrast agent and the concentration distribution image of each chromophore. And step five, by setting a reasonable iteration stop criterion, a high-quality reconstruction result can be stably obtained, meanwhile, the cost of a large amount of time and calculation cost is avoided, and the requirement of high robustness of the used algorithm tool in clinical application is met.
Further, in the first step, the reconstructed multi-wavelength optical absorption coefficient map mu of the target tissue region is inputaThe method specifically comprises the following steps:
an optical detection device (such as an optical endoscope, an optical breast imager, a biopsy microscope and the like) uses 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) commonly used in biological tissue detection to carry out multi-wavelength excitation and imaging on lesion tissues at a certain wavelength interval, the size of an image can be set in the detection device, for example, 35 excitation wavelengths are selected at 5nm step intervals in a wavelength range of 680 nm-884 nm, and an image with the size of 158 x 158 pixels under each wavelength is obtained; in the wave band range, the inherent chromophore of the biological tissue reflecting the physiological change of the tumor area and the exogenous contrast agent have unique spectral absorption characteristics; further acquiring the multispectral optical absorption coefficient mu of the target tissue regionaAnd (4) mapping.
Further, in the second step, blind unmixing pretreatment is carried out by using an ATGP algorithm, and the spectrum prior spectrum of the chromophore inherent in the biological tissue is introduced to initialize the spectral characteristics and the concentration distribution of each chromophore in the target biological tissue.
Further, the second step specifically comprises the following steps:
determining the total number W of chromophore components contained in a region of biological tissue versus a multi-wavelength optical absorption coefficient map muaRealizing blind solution solving to obtain the initial estimation of the absorption spectrum characteristics of each chromophore in the tissue
Figure BDA0002950300100000061
Wherein
Figure BDA0002950300100000062
The dimension of the composite material is L multiplied by W,
Figure BDA0002950300100000063
l is the total number of the wavelengths selected for measurement;
wherein L can be set to 35; because the generation and development process of tumor is closely related to the generation of tumor blood vessel, and hemoglobin is the main inherent chromophore in blood, after fluorescent dye which is approved by FDA and has definite clinical application is added;
comprises three parts of main absorber and HbO in the near infrared spectrum region2Hb, exogenous contrast agent such as ICG, so that the total number W of chromophore components contained in the region of biological tissue is determined to be 3;
calculation of absorption spectra of various chromophores within tissue preliminary estimated by ATGP
Figure BDA0002950300100000064
Prior spectral characteristics of the intrinsic chromophore
Figure BDA0002950300100000065
Spectral angle error between
Figure BDA0002950300100000066
Measuring the similarity between the spectra; when P is the total number of intrinsic chromophore components in the biological tissue, then
Figure BDA0002950300100000067
The dimension of the method is L multiplied by P,
Figure BDA0002950300100000068
HbO2and Hb as an inherent chromophore in the target biological tissue region, the spectral characteristics of which are taken as prior spectra, so P is 2;
setting a spectral angle error threshold angle _ tol, and if ATGP is preliminarily estimated, setting the absorption spectrum of each chromophore in the tissue
Figure BDA0002950300100000071
With a certain prior spectrum
Figure BDA0002950300100000072
The spectral angle error therebetween is less than a set threshold, i.e.
Figure BDA0002950300100000073
Then a priori spectra are usedReplacing the spectrum of the corresponding chromophore analyzed from the pretreatment result to improve the accuracy of the subsequent method for calculating the unknown spectrum and concentration of the exogenous contrast agent, and obtaining the initialized spectral characteristics epsilon of each chromophore in the target biological tissue after the step is finished;
obtaining a multi-wavelength optical absorption coefficient map muaAnd various chromophores in biological tissues, HbO2The initialized spectral characteristics epsilon of Hb and ICG, and solving the initialized concentration distribution C of each chromophore in the biological tissue by utilizing a generalized inverse algorithm to complete all initialization steps:
Figure BDA0002950300100000074
further, the preliminary estimation of absorption spectra of the three types of chromophores is obtained
Figure BDA0002950300100000075
The method comprises the following steps:
extracting the nth target absorption spectrum t in the nth iteration by using ATGPnSatisfies the following formula:
Figure BDA0002950300100000076
in the formula (I), the compound is shown in the specification,
Figure BDA0002950300100000077
is an orthogonal projection matrix; i is an identity matrix; u shapen-1＝[t1,t2,…,tn-1]Generating a target absorption spectrum matrix for the n-1 th calculation; r is a vector of image pixels.
Further, the formula for calculating the spectral angle error is shown as follows:
Figure BDA0002950300100000078
further, in the third step, an optimized objective function form is constructed, matched prior spectra are introduced, and sparse regularization terms of the spatial distribution of the concentration of each chromophore are added;
the specific process is as follows:
in the technical field of biomedical optical function imaging, the contribution of each chromophore of a tissue under a certain wavelength is distinguished according to the difference of spectral characteristics; the multi-wavelength optical absorption coefficient map can be expressed as:
Figure BDA0002950300100000079
in the formula, muaThe optical absorption coefficient of a certain pixel position in a target area under a certain wavelength;
Figure BDA0002950300100000081
the molar extinction coefficient of the mth chromophore at the measurement wavelength; cm(r) represents the concentration of the mth chromophore at the site; n is the total number of pixels of the optical absorption coefficient image reconstructed at each wavelength, and N may be set to 158 × 158 — 24964, where the expression may be expressed as a matrix multiplication form, that is: mu.sa＝εC；
μaThe optical absorption coefficient represented by a certain pixel in the map is the 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 constructed objective function form, the fidelity term is divided to match the prior spectrum
Figure BDA0002950300100000082
And the in-vivo spectrum to be obtained
Figure BDA0002950300100000083
As two parts of the coefficient matrix, the concentration distribution of each chromophore in the tissue is correspondingly separated into intrinsic chromophore concentration distribution vectors
Figure BDA0002950300100000084
And exogenous contrast agent concentration distribution vector
Figure BDA0002950300100000085
And using the pairs of total variation regularization terms
Figure BDA0002950300100000086
And
Figure BDA0002950300100000087
implementing constraints of spatial distribution;
wherein
Figure BDA0002950300100000088
The dimension of the composite material is L multiplied by Q,
Figure BDA0002950300100000089
q is the total number of exogenous chromophore components; i.e., Q may be set to 1;
as described above
Figure BDA00029503001000000810
The dimension is P multiplied by N,
Figure BDA00029503001000000811
Figure BDA00029503001000000812
the dimension is Q multiplied by N,
Figure BDA00029503001000000813
the constructed objective function is:
Figure BDA00029503001000000814
in the formula, | · the luminance | |TVIs a full variation norm, and lambda and xi are regularization parameters; in the course of the calculation process,
Figure BDA00029503001000000815
vector sum
Figure BDA00029503001000000816
The vectors need to be 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 the step of selecting a search step length in each iteration is avoided in the optimized reconstruction process;
in particular, the proposed iterative update in the optical spectrum of an exogenous contrast agent (such as ICG) in target biological tissue
Figure BDA00029503001000000817
And the concentration distribution of the individual chromophores: (
Figure BDA00029503001000000818
And
Figure BDA00029503001000000819
) The update formula of (2) is:
iterative solution of unknown in-vivo spectral features for the j-th exogenous contrast agent, where j ∈ [1, Q ]:
Figure BDA0002950300100000091
iterative solution of the concentration distribution of the j-th exogenous contrast agent, where j ∈ [1, Q ]:
Figure BDA0002950300100000092
iterative solution of the ith intrinsic chromophore concentration distribution, where i ∈ [1, P ]:
Figure BDA0002950300100000093
in the above-mentioned formula, the compound of formula,
Figure BDA0002950300100000094
representing element-by-element division operations within two vectors of the same dimension;
Figure BDA0002950300100000098
represents the Hadamard product operator; k is the current iteration number; lambda [ alpha ]kAnd xikSelecting a regularization parameter for the kth iteration, wherein the parameter is selected to ensure the convergence of the updating formula;
Figure BDA0002950300100000095
representing the vector representation after rearranging the concentration distribution vector of the jth exogenous contrast agent into an image matrix and carrying out total variation partial derivative operation on the image matrix;
Figure BDA0002950300100000096
the vector representation is obtained by rearranging the concentration distribution vector of the i-th intrinsic chromophore into an image matrix and performing a partial derivative operation of the total variation on the image matrix.
Further, in the fifth step, the iteration number k is updated to k +1, whether an iteration stop condition is reached is judged, and if not, the fourth step is returned; otherwise, stopping iteration, and outputting the final reconstructed in-vivo spectrum of the exogenous contrast agent and the concentration distribution image of each chromophore.
Further, the specific determination manner of the iteration stop condition is as follows:
calculating the difference error of the fidelity terms in the objective function in the two adjacent iteration updating processes, setting an iteration stop criterion tol, stopping iteration if the error is less than or equal to tol, and performing iteration on the basis of the criterion tol
Figure BDA0002950300100000097
Outputting the reconstruction result;
tol value set to 10-4When the iteration times k reach the iteration stop condition, the unknown in-vivo spectrum of the exogenous contrast agent and the synchronous acquisition of the concentration distribution images of all chromophores are finally realized;
in the two adjacent iterative updating processes, the formula for calculating the difference error of the fidelity terms 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 carries out initialization operation to estimate the spectral characteristics and the concentration distribution of each chromophore in target biological tissues, and comprises the following steps: blind unmixing pretreatment is carried out on the multi-wavelength optical absorption coefficient map by adopting an automatic target generation Algorithm (ATGP), the spectrum unmixing result and the spectrum (prior spectrum) of inherent chromophores (such as oxygenated/deoxygenated hemoglobin, water, lipid, melanin and the like) in biological tissues are subjected to spectrum angle matching calculation, and the prior spectrum is used for replacing the corresponding chromophore spectrum analyzed from the pretreatment result so as to improve the accuracy of the follow-up steps on the unknown in-vivo spectrum and concentration calculation of the exogenous contrast agent; then, constructing an optimized solving objective function, wherein the fidelity term of the solving objective function is divided into two parts which take the matching prior spectrum and the spectrum to be solved as a coefficient matrix, and introducing space constraint on the concentration distribution of each chromophore by utilizing a total variation regularization Term (TV); and finally, the objective function is inverted by using the proposed iteration solving framework based on non-negative matrix factorization, so that the unknown exogenous contrast agent is synchronously obtained in the spectrum and the concentration distribution of each chromophore, the convergence of the 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, and the spectrum of the inherent chromophore (such as oxygenated/deoxygenated hemoglobin, water, lipid, melanin and the like) in the biological tissues is known and can be used as a prior spectrum and introduced into the initialization step by using a spectral angle matching mode, so that the accuracy of the subsequent steps on the unknown in-vivo spectral spectrum and concentration calculation of the exogenous contrast agent is improved. The invention constructs an optimized solving objective function, the fidelity term of the solving objective function is divided into two parts which take the matching prior spectrum and the spectrum to be solved as a coefficient matrix, and the space constraint on the concentration distribution of each chromophore is introduced by utilizing a 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 iterative solution framework provided by the invention is used for inverting the objective function, so that the unknown in-vivo spectrum of the exogenous contrast agent and the synchronous reconstruction of the concentration distribution images of all chromophores can be realized, and meanwhile, the provided iterative 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 image reconstruction of chromophore concentrations in target tissues are compared in table 1 for the method of the invention with the method not involving an effective sparsity constraint. The result shows that effective space sparse regularization is beneficial to improving reconstruction quality, and even if noise influence exists in simulation synthetic data, the method can be used for reconstructing a concentration distribution image with high fidelity.
TABLE 1
Figure BDA0002950300100000111
Fig. 5 compares the reconstruction of the unknown spectrum of the exogenous contrast agent by the method of the present invention with the common blind unmixing algorithm (VCA and ATGP). The result shows that due to the sensitivity to noise interference, the spectral characteristics acquired by using the blind unmixing algorithm (such as VCA and ATGP) are different from the accurate spectral characteristics. In contrast, the spectral features reconstructed by the method of the present invention are consistent with the preset spectral profile, which indicates that in the steps of the method of the present invention, the spectrum of the intrinsic chromophore in the biological tissue is used as the prior spectrum, which is beneficial to accurately estimate the unknown spectral features of the exogenous contrast agent.
Drawings
In order to more clearly illustrate the technical solutions of the embodiments of the present application, the drawings needed to be used in the embodiments of the present application will be briefly described below, and it is obvious that the drawings described below are only some embodiments of the present application, and it is obvious for those skilled in the art that other drawings can be obtained from the drawings without creative efforts.
Fig. 1 is a schematic flowchart of an optical functional imaging method combining spatial regularization and semi-blind spectral unmixing according to an embodiment of the present invention.
Fig. 2 is a schematic flowchart of an implementation of the optical function imaging method combining spatial regularization and semi-blind spectral unmixing according to an embodiment of the present invention.
FIG. 3 is a graph of multi-wavelength optical absorption coefficient μ provided by an embodiment of the present inventionaSchematic of the image at 4 wavelengths.
FIG. 4 is a graph illustrating the results of reconstructing an image for resolving the unknown spectral features and concentration distribution of individual chromophores in an exogenous contrast agent at a noise level of 20dB SNR according to an embodiment of the present invention.
In fig. 4: panels (a) to (c) are graphs simulating oxyhemoglobin (HbO) in the tumor microenvironment2) True concentration distribution images of deoxyhemoglobin (Hb) and exogenous contrast agent (indocyanine green ICG);
FIGS. (d) to (f) are diagrams for analyzing each chromophore (oxyhemoglobin HbO) by the method proposed by the present invention2Deoxyhemoglobin Hb, and exogenous contrast agent ICG).
Fig. 5 is a spectral feature diagram for reconstructing an exogenous contrast agent ICG at a noise level with a signal-to-noise ratio of 20dB, provided by an embodiment of the present invention.
Fig. 6 is a graph of a change trend of a difference value of an objective function in two adjacent iterative updating processes when an iteration stop condition is reached in the 278 th iteration performed by the iterative solution framework according to the embodiment of the present invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is further described in detail with reference to the following embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.
Aiming at the problems in the prior art, the invention provides an optical function imaging method combining spatial regularization and semi-blind spectral 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 the spatial regularization and the semi-blind spectrum unmixing strategy, provided by the embodiment of the invention, comprises the following steps:
s101: input multi-wavelength optical absorption coefficient map mua；
S102: initializing the spectral characteristics and concentration distribution of each chromophore in the target biological tissue;
s103: constructing an optimized objective function form, introducing a matched prior spectrum, and adding a sparse regularization term of the concentration spatial distribution of each chromophore;
s104: the proposed iterative solution framework is utilized to invert the objective function, and the in-vivo spectral characteristics of the external contrast agent in the biological tissue and the concentration distribution image of each chromophore are updated iteratively;
s105: judging whether an iteration stop condition is reached, if not, returning to the step S104; otherwise, stopping the iteration and outputting the spectral characteristics of the exogenous contrast agent and the concentration distribution image of each chromophore.
A person having ordinary skill in the art can also perform 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 shown in fig. 1 is only one specific embodiment.
In S101 provided by the embodiment of the present invention, a target tissue region multi-wavelength optical absorption coefficient map μ obtained by reconstruction is inputaThe method specifically comprises the following steps:
designing a microenvironment of a blood vessel phantom simulation tumor area, wherein the size of a phantom image is 158 multiplied by 158 pixels; in the experiment, 35 excitation wavelengths are selected at 5nm step intervals in the wavelength range of 680nm to 884 nm; in this band, a chromophore intrinsic to biological tissue (oxyhemoglobin HbO) that reflects physiological changes in the tumor region2Deoxyhemoglobin Hb) and exogenous contrast agents (indocyanine green ICG) have unique spectral absorption characteristics; in simulated multi-spectral optical absorption distributionWhite Gaussian noise is superimposed on the image, in this way a data set mu with a signal-to-noise ratio of 20dB is formedaAnd (4) mapping. Optionally 4 wavelengths of muaImage, as shown in fig. 3.
In S102 provided by the embodiment of the invention, blind unmixing pretreatment is carried out by utilizing an ATGP algorithm, and the spectrum prior spectrum of each chromophore in the biological tissue is introduced to initialize the spectral characteristics and the concentration distribution of each chromophore in the target biological tissue;
the specific process is as follows:
(1) determining the total number W of chromophore components contained in a region of biological tissue versus a multi-wavelength optical absorption coefficient map muaRealizing blind unmixing solution to obtain the initial estimation of the absorption spectrum characteristics (molar extinction coefficient) of each chromophore in the tissue
Figure BDA0002950300100000131
Wherein
Figure BDA0002950300100000132
The dimension of the composite material is L multiplied by W,
Figure BDA0002950300100000133
l is the total number of the wavelengths selected for measurement;
wherein L is 35; because the generation and development process of tumor is closely related to the generation of tumor blood vessel, and hemoglobin is the main inherent chromophore in blood, after fluorescent dye which is approved by FDA and has definite clinical application is added;
comprises three main absorbers (HbO) in the near infrared spectral region2Hb, exogenous contrast agent ICG) so that the total number W of chromophore components contained in the region of biological tissue is determined to be 3.
The described preliminary estimation of absorption spectra of these three types of chromophores
Figure BDA0002950300100000134
The method comprises the following steps:
extracting the nth target absorption spectrum t in the nth iteration by using ATGPnSatisfies the following formula:
Figure BDA0002950300100000141
in the formula (I), the compound is shown in the specification,
Figure BDA0002950300100000142
is an orthogonal projection matrix; i is an identity matrix; u shapen-1＝[t1,t2,…,tn-1]Generating a target absorption spectrum matrix for the n-1 th calculation; r is a vector of image pixels.
(2) Calculation of absorption spectra of various chromophores within tissue preliminary estimated by ATGP
Figure BDA0002950300100000143
Prior spectral characteristics of the intrinsic chromophore
Figure BDA0002950300100000144
Spectral angle error between
Figure BDA0002950300100000145
Measuring the similarity between the spectra; when P is the total number of intrinsic chromophore components in the biological tissue, then
Figure BDA0002950300100000146
The dimension of the method is L multiplied by P,
Figure BDA0002950300100000147
HbO2and Hb as the inherent chromophore in the target biological tissue region, the spectral signature of which is used as the prior spectrum, so that P is 2 in this embodiment.
(3) Setting a spectral angle error threshold angle _ tol, and if ATGP is preliminarily estimated, setting the absorption spectrum of each chromophore in the tissue
Figure BDA0002950300100000148
With a certain prior spectrum
Figure BDA0002950300100000149
The error of the spectral angle between the two is less than the set threshold valueI.e. by
Figure BDA00029503001000001410
Replacing the spectrum of the corresponding chromophore analyzed from the preprocessing result by using the prior spectrum to improve the accuracy of the subsequent method for calculating the unknown spectrum and concentration of the exogenous contrast agent, and finishing the step to obtain the initialized spectral characteristics epsilon of each chromophore in the target biological tissue:
the formula for calculating the spectral angle error is shown as follows:
Figure BDA00029503001000001411
(4) to this end, a multi-wavelength optical absorption coefficient map μ has been obtainedaAnd chromophores (HbO) within biological tissues2Hb, ICG), the initialized concentration distribution C of each chromophore in 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 prior spectrum is introduced, and a sparse regularization term of the spatial distribution of each chromophore concentration is added;
the specific process is as follows:
in the field of biomedical optical functional imaging technology, the contribution of tissue chromophores at a certain wavelength is distinguished according to the difference of spectral characteristics. The multi-wavelength optical absorption coefficient map can be expressed as:
Figure BDA0002950300100000151
in the formula, muaThe optical absorption coefficient of a certain pixel position in a target area under a certain wavelength;
Figure BDA0002950300100000152
the molar extinction coefficient of the mth chromophore at the measurement wavelength; cm(r) represents the concentration of the mth chromophore at the site; n is the total number of pixels of the optical absorption coefficient image reconstructed at each wavelength, and in this embodiment, N is 158 × 158 — 24964, which can be expressed as a matrix multiplication form, that is: mu.sa＝εC；
μaThe optical absorption coefficient represented by a certain pixel in the map is the 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 objective function form constructed in S103, the fidelity term is divided to match the prior spectrum
Figure BDA0002950300100000153
(HbO2Hb) and in vivo spectrum to be determined
Figure BDA0002950300100000154
(ICG) As two parts of the coefficient matrix, the concentration distribution of each chromophore in the tissue is separated into intrinsic chromophore concentration distribution vectors
Figure BDA0002950300100000155
And exogenous contrast agent concentration distribution vector
Figure BDA0002950300100000156
And utilizing a total variation regularization Term (TV) pair
Figure BDA0002950300100000157
And
Figure BDA0002950300100000158
implementing constraints of spatial distribution;
wherein
Figure BDA0002950300100000159
The dimension of the composite material is L multiplied by Q,
Figure BDA00029503001000001510
q is the total number of exogenous chromophore components; i.e., Q is 1 in this embodiment;
as described above
Figure BDA00029503001000001511
The dimension is P multiplied by N,
Figure BDA00029503001000001512
Figure BDA00029503001000001513
the dimension is Q multiplied by N,
Figure BDA00029503001000001514
the objective function constructed in S103 is:
Figure BDA00029503001000001515
in the formula, | · the luminance | |TVIn the case of fully-variant norms, λ and ξ are regularization parameters. In the course of the calculation process,
Figure BDA00029503001000001516
vector sum
Figure BDA00029503001000001517
The vectors need to be 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 iterative solution framework proposed in the embodiment of the present invention can converge the iterative process of minimizing the objective function, and the step of selecting a search step in each iteration is avoided in the optimized reconstruction process.
In particular, the proposed iterative update of the in-vivo spectral update of exogenous contrast agents (ICG) in target biological tissue
Figure BDA00029503001000001611
And the concentration distribution of the individual chromophores: (
Figure BDA0002950300100000161
And
Figure BDA0002950300100000162
) The update formula of (2) is:
iterative solution of the spectral features of the jth exogenous contrast agent, where j ∈ [1, Q ]:
Figure BDA0002950300100000163
iterative solution of the concentration distribution of the j-th exogenous contrast agent, where j ∈ [1, Q ]:
Figure BDA0002950300100000164
iterative solution of the ith intrinsic chromophore concentration distribution, where i ∈ [1, P ]:
Figure BDA0002950300100000165
in the above-mentioned formula, the compound of formula,
Figure BDA0002950300100000166
representing element-by-element division operations within two vectors of the same dimension;
Figure BDA00029503001000001610
represents the Hadamard product operator; k is the current iteration number; lambda [ alpha ]kAnd xikThe regularization parameter is selected during the kth iteration;
Figure BDA0002950300100000167
representing the vector representation after rearranging the concentration distribution vector of the jth exogenous contrast agent into an image matrix and carrying out total variation partial derivative operation on the image matrix;
Figure BDA0002950300100000168
the vector representation is obtained by rearranging the concentration distribution vector of the i-th intrinsic chromophore into an image matrix and performing a partial derivative operation of the total variation on the image matrix.
In S105 provided in the embodiment of the present invention, the iteration number k is updated to k +1, and it is determined whether an iteration stop condition is met, and if not, the process returns to step S104; otherwise, stopping iteration, and outputting the final reconstructed in-vivo spectrum of the exogenous contrast agent and the concentration distribution image of each chromophore.
The specific judgment mode of the iteration stop condition is as follows:
calculating the difference error of the fidelity terms in the objective function in the two adjacent iteration updating processes, setting an iteration stop criterion tol, stopping iteration if the error is less than or equal to tol, and performing iteration on the basis of the criterion tol
Figure BDA0002950300100000169
And outputting the reconstruction result. In the present embodiment, the tol value is set to 10-4When the iteration number k is 278, the iteration stop condition is reached, and finally the unknown in-vivo spectrum and each chromophore (HbO) of the exogenous contrast agent (ICG) are realized2Hb, ICG) concentration distribution image.
In the two adjacent iterative updating processes, the formula for calculating the difference error of the fidelity terms 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 prior spectrum and adds a space regularization constraint in image reconstruction. In order to quantitatively evaluate the superiority of the spatial regularization constraint in reconstructing the concentration distribution image, the image reconstruction results of the method of the present invention (Priori-NMF-TV) and the method not involving the effective sparsity constraint (Priori-NMF) were compared, and the results are shown in table 1. And selecting the parameter distance (d) as a quantization index to represent the difference between the reconstructed image and the real image. The parameter d is defined as:
Figure BDA0002950300100000172
where u and v represent the reconstructed image and the real image, respectively. The image size is M × N. The smaller the value of d, the smaller the difference between the reconstructed image and the real image.
The invention can also use a handheld or fixed biomedical optical detection device (such as an optical endoscope, a diffusion fluorescence tomography system, a diffusion optical breast tomography system, a multi-wavelength photoacoustic tomography system, a tissue biopsy microscope and the like) to perform multi-wavelength excitation and imaging on body tissues or tissue samples at a certain wavelength interval by using 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), the image size can be set in the detection device, for example, in the wavelength range of 680 nm-884 nm, 35 excitation wavelengths are selected at a step interval of 5nm to obtain images with the size of 158 × 158 pixels under each wavelength; in the wave band range, the inherent chromophore of the biological tissue reflecting the physiological change of the tumor area and the exogenous contrast agent have unique spectral absorption characteristics; further acquiring the multispectral optical absorption coefficient mu of the target tissue regionaAnd (4) mapping. After that, the software written based on the algorithm of the invention is operated in the biomedical optical detection system, the in-vivo concentration distribution image of the endogenous chromophore (hemoglobin) and the exogenous contrast agent (ICG) with the biological marker which can reflect the tissue physiology and pathological state clinically can be further directly obtained, the tissue blood oxygen saturation value can be further calculated, and the 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 to fig. 6, experimental results of the present invention indicate that, in the optical function imaging method combining spatial regularization and semi-blind spectral unmixing proposed by the present invention, an optimized objective function is constructed, that is, the spectrum of intrinsic chromophores in biological tissues is considered as a prior spectrum and spatial regularization constraints are added in image reconstruction, so that under the condition that a multi-wavelength optical absorption coefficient map is affected by noise, unknown in-vivo spectral features of exogenous contrast agents and high-quality reconstructed concentration distribution images of the chromophores can be effectively estimated, and synchronous calculation thereof can be realized. Fig. 5 compares the reconstruction of the unknown spectrum of the exogenous contrast agent by the method of the present invention with the common blind unmixing algorithm (VCA and ATGP). The result shows that due to the sensitivity to noise interference, the spectral characteristics acquired by using the blind unmixing algorithm (such as VCA and ATGP) are different from the accurate spectral characteristics. In contrast, the spectral features reconstructed by the method of the present invention are consistent with the preset spectral profile, which indicates that in the steps of the method of the present invention, the spectrum of the intrinsic chromophore in the biological tissue is used as the prior spectrum, which is beneficial to accurately estimate the unknown spectral features of the exogenous contrast agent. In the iteration process, the iteration solving framework provided by the method avoids the selection of the search step length in each iteration, so that the iteration process convergence of the minimized target function is ensured under the self-adaptive multiplication updating.
It should be noted that the embodiments of the present invention can be realized by 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 specially designed hardware. Those skilled 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 code being provided on a carrier medium such as a disk, CD-or DVD-ROM, programmable memory such as read only memory (firmware), or a data carrier such as an optical or electronic signal carrier, for example. The apparatus and its modules of the present invention may be implemented by hardware circuits such as very large scale integrated circuits or gate arrays, semiconductors such as logic chips, transistors, or programmable hardware devices such as field programmable gate arrays, programmable logic devices, etc., or by software executed by various types of processors, or by a combination of hardware circuits and software, e.g., firmware.
The above description is only for the purpose of illustrating the present invention and the appended claims are not to be construed as limiting the scope of the invention, which is intended to cover all modifications, equivalents and improvements that are within the spirit and scope of the invention as defined by the appended claims.

Claims (10)

1. An optical function imaging method combining spatial regularization and semi-blind spectral unmixing, the optical function imaging method combining spatial regularization and semi-blind spectral unmixing comprising:
input multi-wavelength optical absorption coefficient map mua；
Initializing the spectral characteristics and concentration distribution of each chromophore in the target biological tissue;
constructing an optimized objective function form, introducing a matched prior spectrum, and adding a sparse regularization term of the concentration spatial distribution of each chromophore;
the proposed iterative solution framework is utilized to invert the objective function, and the in-vivo spectral characteristics of the external contrast agent in the biological tissue and the concentration distribution image of each chromophore are updated iteratively;
judging whether an iteration stopping condition is met, if not, returning to utilize the proposed iteration solving framework to invert the target function; otherwise, stopping the iteration and outputting the spectral characteristics of the exogenous contrast agent and the concentration distribution image of each chromophore.
2. The method of claim 1, wherein the target tissue region multi-wavelength optical absorption coefficient map μ directly obtained or reconstructed by the biomedical optical inspection device is inputaThe method specifically comprises the following steps:
the optical detection equipment uses a light source commonly used in biological tissue detection, wherein the wavelength of the light source is within a near-infrared first window wavelength range of 600 nm-900 nm or within a near-infrared second window wavelength range of 1000 nm-1700 nm to carry out multi-wavelength excitation and imaging on lesion tissues at a certain wavelength interval, and the size of an image can be detectedSetting in the equipment, for example, in the wavelength range of 680nm to 884nm, selecting 35 excitation wavelengths at the step interval of 5nm to obtain images with the size of 158 × 158 pixels under each wavelength; in the wave band range, the inherent chromophore of the biological tissue reflecting the physiological change of the tumor area and the exogenous contrast agent have unique spectral absorption characteristics; further acquiring the multispectral optical absorption coefficient mu of the target tissue regionaA map; the optical detection device is an optical endoscope, an optical breast imager or a biopsy microscope.
3. The method for optical functional imaging combining spatial regularization and semi-blind spectral unmixing as defined in claim 1 wherein the spectral features and concentration distributions of individual chromophores within a target biological tissue are initialized using the ATGP algorithm for blind unmixing pre-treatment and introducing the inherent chromophore prior spectra in the biological tissue.
4. The method for optically functional imaging with combined spatial regularization and semi-blind spectral unmixing as claimed in claim 1 wherein the specific process of initializing the spectral features and concentration distributions of individual chromophores within a target biological tissue is as follows:
determining the total number W of chromophore components contained in a region of biological tissue versus a multi-wavelength optical absorption coefficient map muaRealizing blind solution solving to obtain the initial estimation of the absorption spectrum characteristics of each chromophore in the tissue
Figure FDA0002950300090000021
Wherein
Figure FDA0002950300090000022
The dimension of the composite material is L multiplied by W,
Figure FDA0002950300090000023
l is the total number of the wavelengths selected for measurement;
wherein L can be set to 35; because the generation and development process of tumor is closely related to the generation of tumor blood vessel, and hemoglobin is the main inherent chromophore in blood, after fluorescent dye which is approved by FDA and has definite clinical application is added;
comprises three parts of main absorber and HbO in the near infrared spectrum region2Hb, exogenous contrast agent ICG, so that the total number W of chromophore components contained in the region of biological tissue is determined to be 3;
calculation of absorption spectra of various chromophores within tissue preliminary estimated by ATGP
Figure FDA0002950300090000024
Prior spectral characteristics of the intrinsic chromophore
Figure FDA0002950300090000025
Spectral angle error between
Figure FDA0002950300090000026
Measuring the similarity between the spectra; when P is the total number of intrinsic chromophore components in the biological tissue, then
Figure FDA0002950300090000027
The dimension of the method is L multiplied by P,
Figure FDA0002950300090000028
HbO2and Hb as an inherent chromophore in the target biological tissue region, the spectral characteristics of which are taken as prior spectra, so P is 2;
setting a spectral angle error threshold angle _ tol, and if ATGP is preliminarily estimated, setting the absorption spectrum of each chromophore in the tissue
Figure FDA0002950300090000029
With a certain prior spectrum
Figure FDA00029503000900000210
The spectral angle error therebetween is less than a set threshold, i.e.
Figure FDA00029503000900000211
Replacing the spectrum of the corresponding chromophore analyzed from the preprocessing result by using the prior spectrum so as to improve the accuracy of the subsequent steps in calculating the unknown spectrum and concentration of the exogenous contrast agent, and finishing the step to obtain the initialized spectral characteristics epsilon of each chromophore in the target biological tissue;
obtaining a multi-wavelength optical absorption coefficient map muaAnd various chromophores in biological tissues, HbO2The initialized spectral characteristics epsilon of Hb and ICG, and solving the initialized concentration distribution C of each chromophore in the biological tissue by utilizing a generalized inverse algorithm to complete all initialization steps:
Figure FDA00029503000900000212
5. the method of claim 4, wherein the preliminary estimate of absorption spectra of the three types of chromophores is obtained
Figure FDA0002950300090000031
The method comprises the following steps:
extracting the nth target absorption spectrum t in the nth iteration by using ATGPnSatisfies the following formula:
Figure FDA0002950300090000032
in the formula (I), the compound is shown in the specification,
Figure FDA0002950300090000033
is an orthogonal projection matrix; i is an identity matrix; u shapen-1＝[t1,t2,…,tn-1]Generating a target absorption spectrum matrix for the n-1 th calculation; r is a vector of image pixels.
6. The method of claim 4 in combination with spatially regularizing and semi-blind spectral unmixing optical function imaging, wherein said computed spectral angle error is formulated as follows:
Figure FDA0002950300090000034
7. the optical functional imaging method combining spatial regularization and semi-blind spectral unmixing of claim 1, wherein an optimized objective function form is constructed, a matched prior spectrum is introduced, and a sparse regularization term for the spatial distribution of chromophore concentrations is added; the specific process is as follows:
in the technical field of biomedical optical function imaging, the contribution of each chromophore of a tissue under a certain wavelength is distinguished according to the difference of spectral characteristics; the multi-wavelength optical absorption coefficient map is represented as:
Figure FDA0002950300090000035
in the formula, muaThe optical absorption coefficient of a certain pixel position in a target area under a certain wavelength;
Figure FDA0002950300090000036
the molar extinction coefficient of the mth chromophore at the measurement wavelength; cm(r) represents the concentration of the mth chromophore at the site; n is the total number of pixels of the optical absorption coefficient image reconstructed at each wavelength, and N may be set to 158 × 158 — 24964, where the expression may be expressed as a matrix multiplication form, that is: mu.sa＝εC；
μaThe optical absorption coefficient represented by a certain pixel in the map is the 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 constructed objective function form, the fidelity term is divided to match the prior spectrum
Figure FDA0002950300090000041
And the in-vivo spectrum to be obtained
Figure FDA0002950300090000042
As two parts of the coefficient matrix, the concentration distribution of each chromophore in the tissue is correspondingly separated into intrinsic chromophore concentration distribution vectors
Figure FDA0002950300090000043
And exogenous contrast agent concentration distribution vector
Figure FDA0002950300090000044
And using the pairs of total variation regularization terms
Figure FDA0002950300090000045
And
Figure FDA0002950300090000046
implementing constraints of spatial distribution;
wherein
Figure FDA0002950300090000047
The dimension of the composite material is L multiplied by Q,
Figure FDA0002950300090000048
q is the total number of exogenous chromophore components; i.e., Q may be set to 1;
as described above
Figure FDA0002950300090000049
The dimension is P multiplied by N,
Figure FDA00029503000900000410
Figure FDA00029503000900000411
the dimension is Q multiplied by N,
Figure FDA00029503000900000412
the constructed objective function is:
Figure FDA00029503000900000413
in the formula, | · the luminance | |TVIs a full variation norm, and lambda and xi are regularization parameters; in the course of the calculation process,
Figure FDA00029503000900000414
vector sum
Figure FDA00029503000900000415
The vectors need to be rearranged into image matrixes respectively to meet the requirement of image total variation operation.
8. The optical functional imaging method combining spatial regularization and semi-blind spectral unmixing as defined in claim 1 wherein the proposed iterative solution framework enables the iterative process of minimizing the objective function to converge and circumvents the step of selecting search steps in each iteration in the optimized reconstruction process;
in particular, the proposed iterative update of the in-vivo spectra of exogenous contrast agents in target biological tissue
Figure FDA00029503000900000416
And the concentration distribution of each chromophore,
Figure FDA00029503000900000417
and
Figure FDA00029503000900000418
the update formula of (2) is:
iterative solution of unknown in-vivo spectral features for the j-th exogenous contrast agent, where j ∈ [1, Q ]:
Figure FDA00029503000900000419
iterative solution of the concentration distribution of the j-th exogenous contrast agent, where j ∈ [1, Q ]:
Figure FDA00029503000900000420
iterative solution of the ith intrinsic chromophore concentration distribution, where i ∈ [1, P ]:
Figure FDA0002950300090000051
in the above-mentioned formula, the compound of formula,
Figure FDA0002950300090000052
representing element-by-element division operations within two vectors of the same dimension;
Figure FDA0002950300090000057
represents the Hadamard product operator; k is the current iteration number; lambda [ alpha ]kAnd xikSelecting a regularization parameter for the kth iteration, wherein the parameter is selected to ensure the convergence of the updating formula;
Figure FDA0002950300090000053
representing the vector representation after rearranging the concentration distribution vector of the jth exogenous contrast agent into an image matrix and carrying out total variation partial derivative operation on the image matrix;
Figure FDA0002950300090000054
the vector representation is obtained by rearranging the concentration distribution vector of the i-th intrinsic chromophore into an image matrix and performing a partial derivative operation of the total variation on the image matrix.
9. The optical function imaging method combining spatial regularization and semi-blind spectral unmixing as claimed in claim 1, wherein the iteration number k is updated to k +1, determining whether an iteration stop condition is reached, if not, returning to utilize the proposed iteration solution framework to invert the objective function, and iteratively updating the in-vivo spectral features of the exogenous contrast agent in the biological tissue and the concentration distribution images of each chromophore; otherwise, stopping iteration, and outputting the final reconstructed in-vivo spectrum of the exogenous contrast agent and the concentration distribution image of each chromophore;
the specific judgment mode of the iteration stop condition is as follows:
calculating the difference error of the fidelity terms in the objective function in the two adjacent iteration updating processes, setting an iteration stop criterion tol, stopping iteration if the error is less than or equal to tol, and performing iteration on the basis of the criterion tol
Figure FDA0002950300090000055
Outputting the reconstruction result;
tol value set to 10-4When the iteration times k are the same, the iteration stopping condition is achieved, and finally the unknown in-vivo spectrum of the exogenous contrast agent and the synchronous acquisition of the concentration distribution images of all chromophores are realized;
in the two adjacent iterative updating processes, the formula for calculating the difference error of the fidelity terms in the objective function is as follows:
Figure FDA0002950300090000056
10. a biomedical optical function imaging terminal, characterized in that the biomedical optical function imaging terminal is used for implementing the optical function imaging method combining spatial regularization and semi-blind spectral unmixing according to any one of claims 1 to 9.
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 true CN113066142A (en) 2021-07-02
CN113066142B 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)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110596017A (en) * 2019-09-12 2019-12-20 生态环境部南京环境科学研究所 Hyperspectral image soil heavy metal concentration assessment method based on space weight constraint and variational self-coding feature extraction
CN113793646A (en) * 2021-09-29 2021-12-14 上海交通大学 Spectral image unmixing method based on weighted nonnegative matrix decomposition and application thereof

Citations (6)

* 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
US20160278678A1 (en) * 2012-01-04 2016-09-29 The Trustees Of Dartmouth College Method and apparatus for quantitative and depth resolved hyperspectral fluorescence and reflectance imaging for surgical guidance
US20170270349A1 (en) * 2016-03-21 2017-09-21 Xerox Corporation Method and apparatus for generating graphical chromophore maps

Patent Citations (6)

* 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
US20160278678A1 (en) * 2012-01-04 2016-09-29 The Trustees Of Dartmouth College Method and apparatus for quantitative and depth resolved hyperspectral fluorescence and reflectance imaging for surgical guidance
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
US20170270349A1 (en) * 2016-03-21 2017-09-21 Xerox Corporation Method and apparatus for generating graphical chromophore maps

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
YE ZHANG ET AL: "A regularization modification to linear spectral unmixing algorithm", 《2012 IEEE INTERNATIONAL GEOSCIENCE AND REMOTE SENSING SYMPOSIUM》 *
甘玉泉: "高光谱遥感图像光谱解混方法研究及其应用", 《中国博士学位论文全文数据库电子期刊 工程科技II辑》 *
祝伟: "基于深度非负矩阵分解与聚类空间处理的高光谱解混研究", 《中国优秀硕士学位论文全文数据库电子期刊 基础科学辑》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110596017A (en) * 2019-09-12 2019-12-20 生态环境部南京环境科学研究所 Hyperspectral image soil heavy metal concentration assessment method based on space weight constraint and variational self-coding feature extraction
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
CN113793646A (en) * 2021-09-29 2021-12-14 上海交通大学 Spectral image unmixing method based on weighted nonnegative matrix decomposition and application thereof
WO2023051661A1 (en) * 2021-09-29 2023-04-06 上海交通大学 Spectral image de-mixing method based on weighted non-negative matrix factorization, and application thereof
CN113793646B (en) * 2021-09-29 2023-11-28 上海交通大学 Spectral image unmixing method based on weighted non-negative matrix factorization and application thereof

Also Published As

Publication number Publication date
CN113066142B (en) 2023-07-07

Similar Documents

Publication Publication Date Title
Gröhl et al. Deep learning for biomedical photoacoustic imaging: A review
Bench et al. Toward accurate quantitative photoacoustic imaging: learning vascular blood oxygen saturation in three dimensions
CN110946553B (en) Hyperspectral image-based in-vivo tissue optical parameter measurement device and method
Chlis et al. A sparse deep learning approach for automatic segmentation of human vasculature in multispectral optoacoustic tomography
CN111103275B (en) PAT prior information assisted dynamic FMT reconstruction method based on CNN and adaptive EKF
CN113066142B (en) Optical function imaging method combining spatial regularization and semi-blind spectrum unmixing
Fuduli et al. Melanoma detection using color and texture features in computer vision systems
Yang et al. Quantitative photoacoustic blood oxygenation imaging using deep residual and recurrent neural network
Yin et al. Exploring the complementarity of THz pulse imaging and DCE-MRIs: Toward a unified multi-channel classification and a deep learning framework
Pugazhenthi et al. Skin disease detection and classification
Dehner et al. Deep-learning-based electrical noise removal enables high spectral optoacoustic contrast in deep tissue
Yin et al. Pattern identification of biomedical images with time series: Contrasting THz pulse imaging with DCE-MRIs
Ram et al. Three-dimensional segmentation of the ex-vivo anterior lamina cribrosa from second-harmonic imaging microscopy
Attye et al. TractLearn: A geodesic learning framework for quantitative analysis of brain bundles
Combalia et al. Digitally stained confocal microscopy through deep learning
Wu et al. System-level optimization in spectroscopic photoacoustic imaging of prostate cancer
Zhang et al. Explainable liver tumor delineation in surgical specimens using hyperspectral imaging and deep learning
Jüstel et al. Spotlight on nerves: portable multispectral optoacoustic imaging of peripheral nerve vascularization and morphology
Liu et al. RFARN: Retinal vessel segmentation based on reverse fusion attention residual network
Mu et al. Cardiac transmembrane potential imaging with GCN based iterative soft threshold network
Mao et al. Single generative networks for stain normalization and quality enhancement of histological images in digital pathology
US11948296B2 (en) Method and system for assessing fibrosis in a tissue sample
Cruz-Guerrero et al. Multi and hyperspectral image unmixing with spatial coherence by extended blind end-member and abundance extraction
Cano et al. Deep learning assisted classification of spectral photoacoustic imaging of carotid plaques
Vemuri et al. Hyperspectral camera selection for interventional health-care

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