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所得的线性关系转化为以下约束优化问题:s.t.||Ax-b||≤σ
其中, 是非凸目标函数,σ=1e-10代表噪声容限。
4.如权利要求3所述的基于惩罚算法的生物发光断层成像光源重建方法,其特征在于,步骤400具体包括:步骤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的值,用非单调近端梯度方法解出步骤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。
5.如权利要求4所述的基于惩罚算法的生物发光断层成像光源重建方法,其特征在于,步骤430中应用的NPG方法为:步骤431:选择参数Lmax=108,Lmin=1,τ=2,c=10-4,M=4,k=0;
步骤432:选择初始
步骤433:求解:
步骤434:若u的目标函数值小于或等于前M+1次迭代中最大的目标函数值即若:则结束搜寻并用Barzilai-Borwein方法选择下一次内部迭代时非单调近端梯度法的参数并使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个元素是的矩阵。