CN107713995A - 一种基于惩罚算法的生物发光断层成像光源重建方法 - Google Patents
一种基于惩罚算法的生物发光断层成像光源重建方法 Download PDFInfo
- Publication number
- CN107713995A CN107713995A CN201711126597.4A CN201711126597A CN107713995A CN 107713995 A CN107713995 A CN 107713995A CN 201711126597 A CN201711126597 A CN 201711126597A CN 107713995 A CN107713995 A CN 107713995A
- Authority
- CN
- China
- Prior art keywords
- light source
- objective function
- algorithm
- fluorescence intensity
- surface fluorescence
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 50
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 26
- 238000003384 imaging method Methods 0.000 title claims abstract description 16
- 230000029918 bioluminescence Effects 0.000 title claims abstract description 15
- 238000005415 bioluminescence Methods 0.000 title claims abstract description 15
- 230000003287 optical effect Effects 0.000 claims abstract description 9
- 238000003325 tomography Methods 0.000 claims description 16
- 238000009792 diffusion process Methods 0.000 claims description 11
- 238000005457 optimization Methods 0.000 claims description 8
- 230000005540 biological transmission Effects 0.000 claims description 7
- 238000009826 distribution Methods 0.000 claims description 6
- 239000011159 matrix material Substances 0.000 claims description 6
- 230000003121 nonmonotonic effect Effects 0.000 claims description 4
- 238000010276 construction Methods 0.000 claims description 2
- QVRVXSZKCXFBTE-UHFFFAOYSA-N n-[4-(6,7-dimethoxy-3,4-dihydro-1h-isoquinolin-2-yl)butyl]-2-(2-fluoroethoxy)-5-methylbenzamide Chemical compound C1C=2C=C(OC)C(OC)=CC=2CCN1CCCCNC(=O)C1=CC(C)=CC=C1OCCF QVRVXSZKCXFBTE-UHFFFAOYSA-N 0.000 claims description 2
- 230000005855 radiation Effects 0.000 claims description 2
- 238000005516 engineering process Methods 0.000 abstract description 2
- 230000006870 function Effects 0.000 description 15
- 238000010586 diagram Methods 0.000 description 8
- 238000004590 computer program Methods 0.000 description 7
- 238000012986 modification Methods 0.000 description 5
- 230000004048 modification Effects 0.000 description 5
- 108060001084 Luciferase Proteins 0.000 description 3
- 230000001788 irregular Effects 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 238000007619 statistical method Methods 0.000 description 3
- 238000003860 storage Methods 0.000 description 3
- IGXWBGJHJZYPQS-SSDOTTSWSA-N D-Luciferin Chemical compound OC(=O)[C@H]1CSC(C=2SC3=CC=C(O)C=C3N=2)=N1 IGXWBGJHJZYPQS-SSDOTTSWSA-N 0.000 description 2
- CYCGRDQQIOGCKX-UHFFFAOYSA-N Dehydro-luciferin Natural products OC(=O)C1=CSC(C=2SC3=CC(O)=CC=C3N=2)=N1 CYCGRDQQIOGCKX-UHFFFAOYSA-N 0.000 description 2
- BJGNCJDXODQBOB-UHFFFAOYSA-N Fivefly Luciferin Natural products OC(=O)C1CSC(C=2SC3=CC(O)=CC=C3N=2)=N1 BJGNCJDXODQBOB-UHFFFAOYSA-N 0.000 description 2
- 239000005089 Luciferase Substances 0.000 description 2
- DDWFXDSYGUXRAY-UHFFFAOYSA-N Luciferin Natural products CCc1c(C)c(CC2NC(=O)C(=C2C=C)C)[nH]c1Cc3[nH]c4C(=C5/NC(CC(=O)O)C(C)C5CC(=O)O)CC(=O)c4c3C DDWFXDSYGUXRAY-UHFFFAOYSA-N 0.000 description 2
- 238000010521 absorption reaction Methods 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 239000000758 substrate Substances 0.000 description 2
- ZKHQWZAMYRWXGA-KQYNXXCUSA-J ATP(4-) Chemical compound C1=NC=2C(N)=NC=NC=2N1[C@@H]1O[C@H](COP([O-])(=O)OP([O-])(=O)OP([O-])([O-])=O)[C@@H](O)[C@H]1O ZKHQWZAMYRWXGA-KQYNXXCUSA-J 0.000 description 1
- ZKHQWZAMYRWXGA-UHFFFAOYSA-N Adenosine triphosphate Natural products C1=NC=2C(N)=NC=NC=2N1C1OC(COP(O)(=O)OP(O)(=O)OP(O)(O)=O)C(O)C1O ZKHQWZAMYRWXGA-UHFFFAOYSA-N 0.000 description 1
- JLVVSXFLKOJNIY-UHFFFAOYSA-N Magnesium ion Chemical compound [Mg+2] JLVVSXFLKOJNIY-UHFFFAOYSA-N 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- QVGXLLKOCUKJST-UHFFFAOYSA-N atomic oxygen Chemical compound [O] QVGXLLKOCUKJST-UHFFFAOYSA-N 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000002329 infrared spectrum Methods 0.000 description 1
- 229910001425 magnesium ion Inorganic materials 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000007254 oxidation reaction Methods 0.000 description 1
- 229910052760 oxygen Inorganic materials 0.000 description 1
- 239000001301 oxygen Substances 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 230000009261 transgenic effect Effects 0.000 description 1
Classifications
-
- 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/0073—Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence by tomography, i.e. reconstruction of 3D images from 2D projections
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- 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
- G06T11/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B2562/00—Details of sensors; Constructional details of sensor housings or probes; Accessories for sensors
- A61B2562/02—Details of sensors specially adapted for in-vivo measurements
- A61B2562/0233—Special features of optical sensors or probes classified in A61B5/00
- A61B2562/0238—Optical sensor arrangements for performing transmission measurements on body tissue
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Biomedical Technology (AREA)
- Molecular Biology (AREA)
- Geometry (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Radiology & Medical Imaging (AREA)
- Evolutionary Computation (AREA)
- Biophysics (AREA)
- Pathology (AREA)
- Computer Hardware Design (AREA)
- Heart & Thoracic Surgery (AREA)
- Medical Informatics (AREA)
- General Engineering & Computer Science (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Algebra (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Investigating, Analyzing Materials By Fluorescence Or Luminescence (AREA)
Abstract
本发明公开了一种基于惩罚算法的生物发光断层成像光源重建方法,涉及光学分子影像技术领域,采用了保留硬约束并惩罚软约束的方法,采用的方法能够求解非凸的目标函数,运用适当的正则化参数更新方案,自适应得到最优的正则化参数,减小了算法对所设定初始正则化参数的敏感性,运用了非凸lp范数,使得到的解更加稀疏,对光源位置的反演更加准确。
Description
技术领域
本发明涉及光学分子影像技术领域,特别是涉及一种基于惩罚算法的生物发光断层成像光源重建方法。
背景技术
生物发光断层成像(Bioluminescence tomography,BLT)是分子影像学的一个重要成像模态,具体过程是首先通过转基因技术将荧光素酶基因转录到靶细胞的DNA上以表达荧光素酶,再向生物体内注射底物荧光素,在生物体内环境中三磷酸腺苷、氧气及镁离子存在的情况下,荧光素酶会催化荧光素底物发生氧化反应,发射出近红外谱段的光子。产生的光子在生物体内经历吸收、散射等行为,其中有一部分光子逸出生物表面被高灵敏设备捕获,通过对于表面荧光强度数据的测量,能够定性,定量的反应生物体细胞的生理变化。
光在生物组织中的传播可由辐射传输方程(Radiative transport equation,RTE)进行准确描述,但由于辐射传输方程求解较为复杂,一般采用其一阶球谐近似模型,即扩散方程(Diffusion equation,DE)描述光在生物组织中的传播。求解扩散方程有三种方法:解析法、统计法和数值法。通常情况下,任意模型上的扩散方程的解析解很难得到,因此一般采用数值法或统计法。统计法中最常用的是蒙特卡洛(Monto Carlo,MC)方法,由于其本质上是一种随机统计方法,计算运行时间长。所以,目前最常用的求解扩散方程的方法是数值法中的有限元方法。
生物发光断层成像的光源重建问题是典型的逆问题。由于逆问题的病态性,以及CCD相机能探测到的生物发光信号极其有限,所以需要利用正则化方法才能获得近似解。这里常用的正则化算法包括:基于L1范数的重建方法(如不完全变量截断共轭梯度算法(Incomplete variables truncated conjugate gradient,IVTCG)),基于L2范数的重建算法(Tikhonov正则化),基于Lp(0<p<1)范数的重建方法等。然而毫无疑问,上述正则化方法均面临着正则化参数的选择的难题,并且在光源重建的定位精度和算法的稳定性等方面仍存在不足。
发明内容
本发明实施例提供了一种基于惩罚算法的生物发光断层成像光源重建方法,可以解决现有技术中存在的问题。
本发明提供了一种基于惩罚算法的生物发光断层成像光源重建方法,所述方法包括以下步骤:
步骤100,采用CCD相机采集生物体的图像,以此获得生物体的表面荧光强度数据;
步骤200,结合步骤100中获得的表面荧光强度数据及给定的光学特异性参数,根据光在生物体中传输过程的扩散近似模型,采用有限元方法构建生物体内部光源和表面荧光强度的线性关系;
步骤300,将步骤200中生物体内部光源和表面荧光强度的线性关系转化为约束优化问题;
步骤400,对于步骤300中的约束优化问题,采用惩罚算法解出最优解x*;
步骤500,利用MATLAB相关工具包,对成像结果进行展示。
本发明实施例中的一种基于惩罚算法的生物发光断层成像光源重建方法,具有以下效果:
1、采用了保留硬约束并惩罚软约束的方法,采用的方法能够求解非凸的目标函数,如lp范数。
2、运用适当的正则化参数更新方案,自适应得到最优的正则化参数,减小了算法对所设定的正则化参数的敏感性。
3、运用了lp范数,使得到的解更加稀疏。
4、对光源位置的反演更加准确。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的基于惩罚算法的生物发光断层成像光源重建方法的总体流程图;
图2(a)为包含两个相距6mm且半径均为1mm的球状光源的仿真模型及重建结果的YZ平面透视图;
图2(b)为包含两个相距6mm且半径均为1mm的球状光源的仿真模型及重建结果的XZ平面透视图;
图2(c)为包含两个相距6mm且半径均为1mm的球状光源的仿真模型及重建结果的XY平面透视图;
图3(a)为包含两个相距6mm且半径分别为1mm和2mm的球状光源的仿真模型及重建结果的YZ平面透视图;
图3(b)为包含两个相距6mm且半径分别为1mm和2mm的球状光源的仿真模型及重建结果的XZ平面透视图;
图3(c)为包含两个相距6mm且半径分别为1mm和2mm的球状光源的仿真模型及重建结果的XY平面透视图;
图4(a)为包含两个水平相距6mm,垂直相距2mm,且半径均为1mm的球状光源的仿真模型及重建结果的YZ平面透视图;
图4(b)为包含两个水平相距6mm,垂直相距2mm,且半径均为1mm的球状光源的仿真模型及重建结果的XZ平面透视图;
图4(c)为包含两个水平相距6mm,垂直相距2mm,且半径均为1mm的球状光源的仿真模型及重建结果的XY平面透视图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
参照图1,本发明提供了一种基于惩罚算法的生物发光断层成像光源重建方法,该方法包括:
步骤100,采用CCD相机采集生物体的图像,以此获得生物体的表面荧光强度数据,将包括吸收系数、散射系数和各向异性因子在内的光学特异性参数作为先验信息;
定义成像域为一个高25mm,半径20mm的均匀圆柱状仿体,成像目标内部光源发射出光子,穿透仿体并被高灵敏度的CCD相机采集到。
步骤200,结合步骤100中获得的表面荧光强度数据及给定的光学特异性参数,根据光在生物体中传输过程的扩散近似模型,采用有限元方法构建生物体内部光源和表面荧光强度的线性关系。具体方法如下:
步骤210:利用光在生物体中传播的模型以描述光在生物体中的传输过程,所采用的模型为辐射传输方程的一阶球谐近似,即扩散方程;
步骤220:利用网格化软件amira对目标离散化,将成像域剖分为多个四面体单元构成的有限元网格,得到网格化后各节点的坐标,节点和四面体之间的对应关系,节点和四面体的面之间的对应关系;
步骤230:对步骤210中的扩散方程进行离散化,结合给定的光学特异性参数构建表面荧光强度数据和内部光源分布之间的线性关系:
b=Ax
其中,b为表面荧光强度数据,A为系统矩阵,表示表面荧光强度数据和内部光源分布之间的线性关系,x是要求解的内部光源能量分布。
步骤300,将步骤200中生物体内部光源和表面荧光强度的线性关系转化为约束优化问题。具体实现方式如下:
将步骤230所得的线性关系转化为以下约束优化问题:
s.t.||Ax-b||≤σ
其中,是非凸目标函数,σ=1e-10代表噪声容限。
步骤400,对于步骤300中的约束优化问题,采用惩罚算法解出最优解x*。具体实现方式如下:
步骤410:构建步骤300中目标函数更加简化的形式:
其中
fλ,μ(x):=hλ,μ(||Ax-b||2-σ2)
步骤420:设定步骤410中目标函数的初始可行解xfeas=A+b,其中A+=AT(AAT)-1为A的伪逆,设定初始参数λ0=1,μ0=1,ε0=1,ρ=2,σ=1e-10,σ=1e-10,θ=1/ρ,τ=2,令x0=0作为算法的起始点;
步骤430:若则将xk,0的值替换为xfeas的值,用非单调近端梯度(NPG)方法解出步骤410中目标函数的近似驻点,也就是函数的一阶导数与0的距离小于εk的点,即满足:
的步骤410中目标函数的近似驻点xk,其中εk为近似程度;
步骤440:判断得出结果的残差和近似程度是否都小于设定容限,即:
max{(||Axk-b||2-σ2)+,0.01εk}≤tol
其中tol是人工设定的终止条件,若满足则xk就是目标函数的最优解x*,若不满足则执行参数更新,使λk+1=ρλk,μk+1=θμk,εk+1=θεk,xk+1=xk,k=k+1并返回步骤430;
具体的,步骤430中应用的NPG方法为:
步骤431:选择参数Lmax=108,Lmin=1,τ=2,c=10-4,M=4,k=0;
步骤432:选择初始
步骤433:求解:
步骤434:若u的目标函数值稍小于或等于前M+1次迭代中最大的目标函数值即若:
则结束搜寻并用Barzilai-Borwein(BB)方法选择下一次内部迭代时非单调近端梯度方法的参数
并使xk+1=μ,k=k+1,若大于前M+1次迭代中最大的目标函数的值,则更新Lk=τLk并返回步骤433;
步骤435:若xk,l满足的近似一阶最优条件,即:
和
得出xk=xk,l,若不满足,则返回步骤433,其中xk,l代表NPG中内迭代l次,主迭代k次得出的解,Diag(xk,l)代表第i个对角元素是的对角矩阵,|xk,l|p代表第i个元素是的矩阵。
步骤500,利用MATLAB相关工具包,对光源重建结果进行展示。
实验结果分析
结合附图2、3、4,对本发明基于惩罚算法的生物发光断层成像光源重建方法进一步进行说明:
附图2(a)、(b)、(c)中,在半径为20mm,高为24mm的圆柱状仿体中,两光源被放置在高18mm,相距6mm的地方。星型标识点为真实光源位置,不规则部分表示重建的荧光光源。两光源都是半径为1mm的球形光源,可以看出,本方法得到的重建结果对光源中心点的定位准确。
附图3(a)、(b)、(c)中,在半径为20mm,高为24mm的圆柱状仿体中,两光源被放置在高18m,相距6mm的地方,星型标识点为真实光源位置,不规则部分表示重建的荧光光源。两光源的半径分别为1mm和2mm。
附图4(a)、(b)、(c)中,在半径为20mm,高为24mm的圆柱状仿体中,两光源其中一个被放置在高16mm的地方,另一个被放置在高18mm的地方。水平距离为6mm,星型标识点为真实光源位置,不规则部分表示重建的荧光光源。
由附图2、3、4可以看出,对于不同的光源设定,基于本发明的生物发光断层成像光源重建方法所得到的重建结果具有光源中心定位准确、光源形状近似的优点。
本领域内的技术人员应明白,本发明的实施例可提供为方法、系统、或计算机程序产品。因此,本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且,本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器,使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中,使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品,该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上,使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理,从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
尽管已描述了本发明的优选实施例,但本领域内的技术人员一旦得知了基本创造性概念,则可对这些实施例作出另外的变更和修改。所以,所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。
显然,本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样,倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内,则本发明也意图包含这些改动和变型在内。
Claims (5)
1.一种基于惩罚算法的生物发光断层成像光源重建方法,其特征在于,所述方法包括以下步骤:
步骤100,采用CCD相机采集生物体的图像,以此获得生物体的表面荧光强度数据;
步骤200,结合步骤100中获得的表面荧光强度数据及给定的光学特异性参数,根据光在生物体中传输过程的扩散近似模型,采用有限元方法构建生物体内部光源和表面荧光强度的线性关系;
步骤300,将步骤200中生物体内部光源和表面荧光强度的线性关系转化为约束优化问题;
步骤400,对于步骤300中的约束优化问题,采用惩罚算法解出最优解x*;
步骤500,利用MATLAB相关工具包,对光源重建结果进行展示。
2.如权利要求1所述的基于惩罚算法的生物发光断层成像光源重建方法,其特征在于,步骤200具体包括:
步骤210:利用光在生物体中传播的模型以描述光在生物体中的传输过程,所采用的模型为辐射传输方程的一阶球谐近似,即扩散方程;
步骤220:利用网格化软件amira对目标离散化,将成像域剖分为多个四面体单元构成的有限元网格,得到网格化后各节点的坐标,节点和四面体之间的对应关系,节点和四面体的面之间的对应关系;
步骤230:对步骤210中的扩散方程进行离散化,结合给定的光学特异性参数构建表面荧光强度数据和内部光源分布之间的线性关系:
b=Ax
其中,b为表面荧光强度数据,A为系统矩阵,表示表面荧光强度数据和内部光源分布之间的线性关系,x是要求解的内部光源能量分布。
3.如权利要求2所述的基于惩罚算法的生物发光断层成像光源重建方法,其特征在于,步骤300具体包括:
将步骤230所得的线性关系转化为以下约束优化问题:
<mrow>
<munder>
<mrow>
<mi>m</mi>
<mi>i</mi>
<mi>n</mi>
</mrow>
<mi>x</mi>
</munder>
<mi>&Phi;</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
</mrow>
s.t.||Ax-b||≤σ
其中,是非凸目标函数,σ=1e-10代表噪声容限。
4.如权利要求3所述的基于惩罚算法的生物发光断层成像光源重建方法,其特征在于,步骤400具体包括:
步骤410:构建步骤300中目标函数更加简化的形式:
<mrow>
<munder>
<mrow>
<mi>m</mi>
<mi>i</mi>
<mi>n</mi>
</mrow>
<mrow>
<mi>x</mi>
<mo>&Element;</mo>
<msup>
<mi>R</mi>
<mi>n</mi>
</msup>
</mrow>
</munder>
<msub>
<mi>F</mi>
<mrow>
<mi>&lambda;</mi>
<mo>,</mo>
<mi>&mu;</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>:</mo>
<mo>=</mo>
<msub>
<mi>f</mi>
<mrow>
<mi>&lambda;</mi>
<mo>,</mo>
<mi>&mu;</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>&Phi;</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
</mrow>
其中
fλ,μ(x):=hλ,μ(||Ax-b||2-σ2)
<mrow>
<msub>
<mi>h</mi>
<mrow>
<mi>&lambda;</mi>
<mo>,</mo>
<mi>&mu;</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>s</mi>
<mo>)</mo>
</mrow>
<mo>:</mo>
<mo>=</mo>
<mi>&lambda;</mi>
<munder>
<mrow>
<mi>m</mi>
<mi>a</mi>
<mi>x</mi>
</mrow>
<mrow>
<mn>0</mn>
<mo>&le;</mo>
<mi>t</mi>
<mo>&le;</mo>
<mn>1</mn>
</mrow>
</munder>
<mo>{</mo>
<mi>s</mi>
<mi>t</mi>
<mo>-</mo>
<mfrac>
<mi>&mu;</mi>
<mi>s</mi>
</mfrac>
<msup>
<mi>t</mi>
<mn>2</mn>
</msup>
<mo>}</mo>
</mrow>
<mrow>
<mi>&Phi;</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>:</mo>
<mo>=</mo>
<msubsup>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>n</mi>
</msubsup>
<msup>
<mrow>
<mo>|</mo>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
<mo>|</mo>
</mrow>
<mfrac>
<mn>1</mn>
<mn>2</mn>
</mfrac>
</msup>
</mrow>
步骤420:设定步骤410中目标函数的初始可行解xfeas=A+b,其中A+=AT(AAT)-1为A的伪逆,设定初始参数λ0=1,μ0=1,ε0=1,ρ=2,σ=1e-10,σ=1e-10,θ=1/ρ,τ=2,令x0=0作为算法的起始点;
步骤430:若则将xk,0的值替换为xfeas的值,用非单调近端梯度方法解出步骤410中目标函数的近似驻点,也就是函数的一阶导数与0的距离小于εk的点,即满足:
<mrow>
<mi>d</mi>
<mi>i</mi>
<mi>s</mi>
<mi>t</mi>
<mrow>
<mo>(</mo>
<mn>0</mn>
<mo>,</mo>
<mo>&dtri;</mo>
<msub>
<mi>f</mi>
<mrow>
<msub>
<mi>&lambda;</mi>
<mi>k</mi>
</msub>
<mo>,</mo>
<msub>
<mi>&mu;</mi>
<mi>k</mi>
</msub>
</mrow>
</msub>
<mo>(</mo>
<msup>
<mi>x</mi>
<mi>k</mi>
</msup>
<mo>)</mo>
<mo>+</mo>
<mo>&part;</mo>
<mo>(</mo>
<mrow>
<mi>&Phi;</mi>
<mo>+</mo>
<msub>
<mi>&delta;</mi>
<msub>
<mi>s</mi>
<mn>1</mn>
</msub>
</msub>
</mrow>
<mo>)</mo>
<mo>(</mo>
<msup>
<mi>x</mi>
<mi>k</mi>
</msup>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>&le;</mo>
<msub>
<mi>&epsiv;</mi>
<mi>k</mi>
</msub>
</mrow>
的步骤410中目标函数的近似驻点xk,其中εk为近似程度;
步骤440:判断得出结果的残差和近似程度是否都小于设定容限,即:
max{(||Axk-b||2-σ2)+,0.01εk}≤tol
其中tol是人工设定的终止条件,若满足则xk就是目标函数的最优解x*,若不满足则执行参数更新,使λk+1=ρλk,μk+1=θμk,εk+1=θεk,xk+1=xk,k=k+1并返回步骤430。
5.如权利要求4所述的基于惩罚算法的生物发光断层成像光源重建方法,其特征在于,步骤430中应用的NPG方法为:
步骤431:选择参数Lmax=108,Lmin=1,τ=2,c=10-4,M=4,k=0;
步骤432:选择初始
步骤433:求解:
<mrow>
<mi>u</mi>
<mo>=</mo>
<mi>A</mi>
<mi>r</mi>
<mi>g</mi>
<munder>
<mrow>
<mi>m</mi>
<mi>i</mi>
<mi>n</mi>
</mrow>
<mi>x</mi>
</munder>
<mo>{</mo>
<mo><</mo>
<mo>&dtri;</mo>
<mi>f</mi>
<mrow>
<mo>(</mo>
<msup>
<mi>x</mi>
<mi>k</mi>
</msup>
<mo>)</mo>
</mrow>
<mo>,</mo>
<mi>x</mi>
<mo>-</mo>
<msup>
<mi>x</mi>
<mi>k</mi>
</msup>
<mo>></mo>
<mo>+</mo>
<mfrac>
<msub>
<mi>L</mi>
<mi>k</mi>
</msub>
<mn>2</mn>
</mfrac>
<mo>|</mo>
<mo>|</mo>
<mi>x</mi>
<mo>-</mo>
<msup>
<mi>x</mi>
<mi>k</mi>
</msup>
<mo>|</mo>
<msup>
<mo>|</mo>
<mn>2</mn>
</msup>
<mo>+</mo>
<mi>P</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>}</mo>
</mrow>
步骤434:若u的目标函数值小于或等于前M+1次迭代中最大的目标函数值即若:
<mrow>
<mi>F</mi>
<mrow>
<mo>(</mo>
<mi>u</mi>
<mo>)</mo>
</mrow>
<mo>&le;</mo>
<munder>
<mrow>
<mi>m</mi>
<mi>a</mi>
<mi>x</mi>
</mrow>
<mrow>
<msub>
<mrow>
<mo>|</mo>
<mrow>
<mi>k</mi>
<mo>-</mo>
<mi>M</mi>
</mrow>
<mo>|</mo>
</mrow>
<mo>+</mo>
</msub>
<mo>&le;</mo>
<mi>i</mi>
<mo>&le;</mo>
<mi>k</mi>
</mrow>
</munder>
<mi>F</mi>
<mrow>
<mo>(</mo>
<msup>
<mi>x</mi>
<mi>i</mi>
</msup>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mfrac>
<mi>c</mi>
<mn>2</mn>
</mfrac>
<mo>|</mo>
<mo>|</mo>
<mi>u</mi>
<mo>-</mo>
<msup>
<mi>x</mi>
<mi>k</mi>
</msup>
<mo>|</mo>
<msup>
<mo>|</mo>
<mn>2</mn>
</msup>
</mrow>
则结束搜寻并用Barzilai-Borwein方法选择下一次内部迭代时非单调近端梯度法的参数
<mrow>
<msubsup>
<mi>L</mi>
<mi>l</mi>
<mn>0</mn>
</msubsup>
<mo>=</mo>
<mi>m</mi>
<mi>i</mi>
<mi>n</mi>
<mo>{</mo>
<mi>m</mi>
<mi>a</mi>
<mi>x</mi>
<mo>{</mo>
<mfrac>
<mrow>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>l</mi>
</mrow>
</msup>
<mo>-</mo>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>l</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>&rsqb;</mo>
</mrow>
<mi>T</mi>
</msup>
<mo>&lsqb;</mo>
<mo>&dtri;</mo>
<msub>
<mi>f</mi>
<mrow>
<msub>
<mi>&lambda;</mi>
<mi>k</mi>
</msub>
<mo>,</mo>
<msub>
<mi>&mu;</mi>
<mi>k</mi>
</msub>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>l</mi>
</mrow>
</msup>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mo>&dtri;</mo>
<msub>
<mi>f</mi>
<mrow>
<msub>
<mi>&lambda;</mi>
<mi>k</mi>
</msub>
<mo>,</mo>
<msub>
<mi>&mu;</mi>
<mi>k</mi>
</msub>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>l</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>|</mo>
<mo>|</mo>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>l</mi>
</mrow>
</msup>
<mo>-</mo>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>l</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>|</mo>
<msup>
<mo>|</mo>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>,</mo>
<msub>
<mi>L</mi>
<mi>min</mi>
</msub>
<mo>}</mo>
<mo>,</mo>
<msub>
<mi>L</mi>
<mrow>
<mi>m</mi>
<mi>a</mi>
<mi>x</mi>
</mrow>
</msub>
<mo>}</mo>
</mrow>
并使xk+1=μ,k=k+1,若大于前M+1次迭代中最大的目标函数的值,则更新Lk=τLk并返回步骤433;
步骤435:若xk,l满足的近似一阶最优条件,即:
<mrow>
<mfrac>
<mrow>
<mo>|</mo>
<msub>
<mi>F</mi>
<mrow>
<msub>
<mi>&lambda;</mi>
<mi>k</mi>
</msub>
<mo>,</mo>
<msub>
<mi>&mu;</mi>
<mi>k</mi>
</msub>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>l</mi>
</mrow>
</msup>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msub>
<mi>F</mi>
<mrow>
<msub>
<mi>&lambda;</mi>
<mi>k</mi>
</msub>
<mo>,</mo>
<msub>
<mi>&mu;</mi>
<mi>k</mi>
</msub>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>l</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>)</mo>
</mrow>
<mo>|</mo>
</mrow>
<mrow>
<mi>max</mi>
<mo>{</mo>
<mn>1</mn>
<mo>,</mo>
<mo>|</mo>
<msub>
<mi>F</mi>
<mrow>
<msub>
<mi>&lambda;</mi>
<mi>k</mi>
</msub>
<mo>,</mo>
<msub>
<mi>&mu;</mi>
<mi>k</mi>
</msub>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>l</mi>
</mrow>
</msup>
<mo>)</mo>
</mrow>
<mo>|</mo>
<mo>}</mo>
</mrow>
</mfrac>
<mo>&le;</mo>
<mi>min</mi>
<mo>{</mo>
<msubsup>
<mi>&epsiv;</mi>
<mi>k</mi>
<mn>2</mn>
</msubsup>
<mo>,</mo>
<msup>
<mn>10</mn>
<mrow>
<mo>-</mo>
<mn>4</mn>
</mrow>
</msup>
<mo>}</mo>
</mrow>
和
<mrow>
<mo>|</mo>
<mo>|</mo>
<mi>D</mi>
<mi>i</mi>
<mi>a</mi>
<mi>g</mi>
<mrow>
<mo>(</mo>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>l</mi>
</mrow>
</msup>
<mo>)</mo>
</mrow>
<mo>&dtri;</mo>
<msub>
<mi>f</mi>
<mrow>
<mi>&lambda;</mi>
<mo>,</mo>
<mi>&mu;</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>l</mi>
</mrow>
</msup>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>p</mi>
<msup>
<mrow>
<mo>|</mo>
<msup>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>,</mo>
<mi>l</mi>
</mrow>
</msup>
<mo>|</mo>
</mrow>
<mi>p</mi>
</msup>
<mo>|</mo>
<msub>
<mo>|</mo>
<mi>&infin;</mi>
</msub>
<mo>&le;</mo>
<msqrt>
<msub>
<mi>&epsiv;</mi>
<mi>k</mi>
</msub>
</msqrt>
</mrow>
得出xk=xk,l,若不满足,则返回步骤433,其中xk,l代表NPG方法中内迭代l次,主迭代k次得出的解,Diag(xk,l)代表第i个对角元素是的对角矩阵,|xk,l|p代表第i个元素是的矩阵。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201711126597.4A CN107713995B (zh) | 2017-11-15 | 2017-11-15 | 一种基于惩罚算法的生物发光断层成像光源重建方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201711126597.4A CN107713995B (zh) | 2017-11-15 | 2017-11-15 | 一种基于惩罚算法的生物发光断层成像光源重建方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN107713995A true CN107713995A (zh) | 2018-02-23 |
CN107713995B CN107713995B (zh) | 2020-09-01 |
Family
ID=61215673
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201711126597.4A Active CN107713995B (zh) | 2017-11-15 | 2017-11-15 | 一种基于惩罚算法的生物发光断层成像光源重建方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN107713995B (zh) |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112089434A (zh) * | 2020-10-16 | 2020-12-18 | 陕西师范大学 | 一种多光谱生物发光断层成像方法和系统 |
CN113159715A (zh) * | 2021-04-06 | 2021-07-23 | 杭州远传新业科技有限公司 | 客服坐席排班方法、系统、电子装置和存储介质 |
CN113674375A (zh) * | 2021-08-19 | 2021-11-19 | 北京航空航天大学 | 一种用于半透明材质渲染的次表面散射计算方法 |
CN114220489A (zh) * | 2021-12-15 | 2022-03-22 | 北京应用物理与计算数学研究所 | 原子结构弛豫的非单调线搜索方法及装置 |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20110282181A1 (en) * | 2009-11-12 | 2011-11-17 | Ge Wang | Extended interior methods and systems for spectral, optical, and photoacoustic imaging |
CN105326475A (zh) * | 2015-09-16 | 2016-02-17 | 西北大学 | 一种基于多光源分辨的生物发光断层成像重建方法 |
CN105455780A (zh) * | 2015-11-17 | 2016-04-06 | 西北大学 | 基于双网格的有限投影的荧光分子断层成像重建方法 |
CN105581779A (zh) * | 2015-12-13 | 2016-05-18 | 北京工业大学 | 一种直接融合结构成像的生物发光断层成像重建的方法 |
CN106097441A (zh) * | 2016-06-25 | 2016-11-09 | 北京工业大学 | 基于l1范数与tv范数的复合正则化生物发光断层成像重建方法 |
CN107220961A (zh) * | 2017-06-14 | 2017-09-29 | 西北大学 | 一种基于半阈值追踪算法的荧光分子断层成像重建方法 |
-
2017
- 2017-11-15 CN CN201711126597.4A patent/CN107713995B/zh active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20110282181A1 (en) * | 2009-11-12 | 2011-11-17 | Ge Wang | Extended interior methods and systems for spectral, optical, and photoacoustic imaging |
CN105326475A (zh) * | 2015-09-16 | 2016-02-17 | 西北大学 | 一种基于多光源分辨的生物发光断层成像重建方法 |
CN105455780A (zh) * | 2015-11-17 | 2016-04-06 | 西北大学 | 基于双网格的有限投影的荧光分子断层成像重建方法 |
CN105581779A (zh) * | 2015-12-13 | 2016-05-18 | 北京工业大学 | 一种直接融合结构成像的生物发光断层成像重建的方法 |
CN106097441A (zh) * | 2016-06-25 | 2016-11-09 | 北京工业大学 | 基于l1范数与tv范数的复合正则化生物发光断层成像重建方法 |
CN107220961A (zh) * | 2017-06-14 | 2017-09-29 | 西北大学 | 一种基于半阈值追踪算法的荧光分子断层成像重建方法 |
Cited By (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112089434A (zh) * | 2020-10-16 | 2020-12-18 | 陕西师范大学 | 一种多光谱生物发光断层成像方法和系统 |
CN112089434B (zh) * | 2020-10-16 | 2024-05-03 | 陕西师范大学 | 一种多光谱生物发光断层成像方法和系统 |
CN113159715A (zh) * | 2021-04-06 | 2021-07-23 | 杭州远传新业科技有限公司 | 客服坐席排班方法、系统、电子装置和存储介质 |
CN113674375A (zh) * | 2021-08-19 | 2021-11-19 | 北京航空航天大学 | 一种用于半透明材质渲染的次表面散射计算方法 |
CN113674375B (zh) * | 2021-08-19 | 2023-10-27 | 北京航空航天大学 | 一种用于半透明材质渲染的次表面散射计算方法 |
CN114220489A (zh) * | 2021-12-15 | 2022-03-22 | 北京应用物理与计算数学研究所 | 原子结构弛豫的非单调线搜索方法及装置 |
CN114220489B (zh) * | 2021-12-15 | 2022-12-02 | 北京应用物理与计算数学研究所 | 原子结构弛豫的非单调线搜索方法及装置 |
Also Published As
Publication number | Publication date |
---|---|
CN107713995B (zh) | 2020-09-01 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN107713995B (zh) | 一种基于惩罚算法的生物发光断层成像光源重建方法 | |
CN109191564B (zh) | 基于深度学习的激发荧光断层成像三维重建方法 | |
Hazelbag et al. | Calibration of individual-based models to epidemiological data: A systematic review | |
Klose et al. | In vivo bioluminescence tomography with a blocking‐off finite‐difference method and MRI/CT coregistration | |
Liu et al. | Two schemes for quantitative photoacoustic tomography based on Monte Carlo simulation | |
CN110974166B (zh) | 基于k近邻局部连接网络的光学断层成像方法、系统 | |
CN108451508B (zh) | 基于多层感知机的生物自发荧光三维成像方法 | |
CN107358653A (zh) | 成像重建方法及装置 | |
Huang et al. | Fast and robust reconstruction method for fluorescence molecular tomography based on deep neural network | |
Cong et al. | A Born‐type approximation method for bioluminescence tomography | |
CN113409466B (zh) | 基于gcn残差连接网络的激发荧光断层成像方法 | |
WO2017004851A1 (zh) | 基于多任务贝叶斯压缩感知方法的生物发光断层成像重建算法 | |
Dlask et al. | Short-time fractal analysis of biological autoluminescence | |
CN104586366B (zh) | 一种激发荧光断层成像重建方法 | |
CN110020710A (zh) | 一种基于人工蜂群算法的射束方向及权重多目标优化方法 | |
CN102034266A (zh) | 激发荧光断层成像的快速稀疏重建方法和设备 | |
Leng et al. | Mathematical method in optical molecular imaging | |
Wang et al. | Fluorescence molecular tomography reconstruction of small targets using stacked auto-encoder neural networks | |
CN109166103A (zh) | 基于多层感知网络的激发荧光断层成像方法 | |
Murad et al. | Reconstruction and localization of tumors in breast optical imaging via convolution neural network based on batch normalization layers | |
Lou et al. | Analysis of the use of unmatched backward operators in iterative image reconstruction with application to three-dimensional optoacoustic tomography | |
CN112089434B (zh) | 一种多光谱生物发光断层成像方法和系统 | |
Kumar et al. | Monte Carlo method for bioluminescence tomography | |
Guo et al. | Sparse reconstruction for bioluminescence tomography based on the semigreedy method | |
CN107374588B (zh) | 一种基于同步聚类的多光源荧光分子断层成像重建方法 |
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 |