CN108805975A - 一种基于改进迭代收缩阈值算法的微观3d重建方法 - Google Patents

Abstract

本发明公开了一种基于改进迭代收缩阈值算法的微观3D重建方法，包括设定相对模糊度初值，计算热扩散方程的解，求解能量方程判断是否小于能量阈值，小于能量阈值时直接得到深度信息值，不小于能量阈值时计算加速算子带入迭代公式，重新进行梯度估计，利用割线线性搜索求出的最优迭代步长带入迭代公式，再求出深度信息值。该方法能加快离焦深度恢复动态优化问题的收敛速度，提高微观3D形貌重建的精度。

Description

一种基于改进迭代收缩阈值算法的微观3D重建方法

技术领域

本发明涉及一种微观3D重建方法，尤其是涉及一种基于改进迭代收缩阈值算法的微观3D重建方法。

背景技术

在微纳米计算机视觉中，基于视觉的微纳米尺度3D形貌重建，对于较全面的理解样品特性和对操作过程进行评估具有重要意义。比较常用的微观3D重建方法主要有体积恢复方法、立体视觉深度恢复方法、聚焦深度恢复方法以及离焦深度恢复方法。与其他的3D重建方法相比，离焦深度恢复方法由于所需图片少、设备简单、操作便利等优点，近些年得到广泛关注与深入研究。离焦深度恢复最早由Pentland提出，该方法利用二维图像的离焦程度特征与景物深度之间的映射关系，反解出景物的三维深度信息。

在景物的离焦深度恢复方法中，首先必须获取景物不同程度的离焦图像，这就导致要改变摄像机参数。但是，在微纳米图像观测中，观测空间非常有限，而且使用的是具有高放大倍数的相机，所以相机的成像模型会随着摄像机参数改变而改变。限于微纳米观测中的条件限制，研究者提出了一种新的基于参数固定的单目视觉摄像机的三维离焦形貌恢复方法。这种方法利用相对模糊度及热辐射方程来求解景物的深度信息，并将深度信息的求解转化成动态优化问题来求解。

但是，在求解离焦深度恢复动态优化问题时，研究者采用经典的迭代收缩阈值算法(ISTA)来求解，由于ISTA算法是梯度下降法的延伸，迭代过程仅考虑当前点的信息进行梯度估计更新迭代点，优化过程呈“之”字形向极小值点靠近，收敛速度比较慢，而且该算法在迭代过程采用固定步长，在迭代过程靠近极小值点时，此时的固定步长可能大于实际的迭代步长，导致算法迭代效率不佳，从而使得重建的微观3D形貌精度不高。

发明内容

针对现有技术的不足，本发明的目的是提供一种基于改进迭代收缩阈值算法的微观3D重建方法，解决离焦深度恢复动态优化过程收敛效率不佳，导致重建的微观3D形貌精度不佳的问题。

本发明技术方案如下，一种基于改进迭代收缩阈值算法的微观3D重建方法，依次包括以下步骤：

步骤一、初始化基本参数，所述基本参数包括第一模糊图像和第二模糊图像；

步骤二、设定第一模糊图像和第二模糊图像的相对模糊度的初始值，并将第二模糊图像以第一模糊图像和相对模糊度表示；

步骤三、将第一模糊图像和第二模糊图像以图像的模糊化用卷积的形式表示并改写为热扩散方程组，计算热扩散方程的解；

步骤四、利用步骤三得到的热扩散方程的解计算能量方程，如果能量小于能量阈值，则由步骤二得到的相对模糊度带入深度信息表达式计算深度信息值；否则，转入步骤五；

步骤五、计算加速算子y_j+1，y_j+1表示由当前的相对模糊度开始，沿着Δσ² _j-Δσ² _j-1所构成的搜索方向进行步长为a_j所得到的相对模糊度，Δσ² _j代表当前的相对模糊度，Δσ² _j-1代表前一次的相对模糊度，y_j+1＝Δσ² _j+(n_j-1/n_j+1)(Δσ² _j-Δσ² _j-1)，a_j＝n_j-1/n_j，n₁＝α₁＝1，0＜α_k＜1；

步骤六、将加速算子y_j+1带入迭代公式重新进行梯度估计；

步骤七、利用割线线性搜索求出最优迭代步长，带入迭代公式求出此时的相对模糊度；

步骤八、将计算所得的此时的相对模糊度带入深度信息表达式，求解出深度信息值。

进一步的，所述深度信息表达式为：

其中f为摄像机焦距，v为摄像机像距，D为凸透镜半径，λ为模糊度与模糊圆之间的常数，s₂(x,y)为第二模糊图像的深度信息值，s₁(x,y)为第一模糊图像的深度信息值，Δσ²(x,y)为相对模糊度。

进一步的，所述第二模糊图像以第一模糊图像和相对模糊度表示是

其中I₂(x,y)为第二模糊图像，I₁(x,y)为第一模糊图像，Δσ²为相对模糊度；所述将第一模糊图像和第二模糊图像以图像的模糊化用卷积的形式表示并改写为热扩散方程组为：

其中为梯度算子，为微分算子，t为扩散时间。

进一步的，所述步骤四中，所述能量方程表示如下

其中I₂(x,y)为第二模糊图像，能量系数ρ＞0，l＞0，z(x,y,Δt)为热扩散方程的解，s为由相对模糊度带入深度信息表达式计算的深度信息值。

进一步的，所述步骤六中，所述迭代公式由以下方法确定：对于无约束最优化问题

E₁(s)＝∫∫(z(x,y,Δt)-I₂(x,y))²dxdy

解决上式的梯度算法为

其中，t_j表示迭代步长，而上式的二阶估计模型可以表示为

加入Tikhonov正则项得到以下迭代公式

忽略常数项，得到

Δσ² _j的计算变成解决每个元素的一维问题，它的迭代公式为

其中，T_ρ:Rⁿ→Rⁿ是收缩算子，定义为

将计算得到的加速算子带入迭代公式得到

进一步的，所述步骤七中，利用割线线性搜索求出最优迭代步长t^*的方法是：建立关于步长的一元函数

t^*＝argminψ(t)，

其中ψ(t)是关于t的一元函数，t是迭代步长，d是迭代方向，t^*是最优迭代步长，采用割线迭代格式表示为

本发明针对离焦深度恢复动态优化问题收敛效率不佳的问题，首先，在每一次迭代过程中，由当前点和前一个点的线性组合构成加速算子重新进行梯度估计，更新迭代点；其次，为了改变迭代步长固定的限制，引入割线线性搜索，动态确定每次最优迭代步长，加快算法收敛速度；最后将改进的迭代收缩阈值算法用于求解离焦深度恢复动态优化问题，加快离焦深度恢复动态优化问题的收敛速度和提高微观3D形貌重建的精度。

本发明所提供的技术方案的优点在于：

1、加快离焦深度恢复动态优化问题的收敛速度；

2、提高微观3D形貌重建的精度。

附图说明

图1为本发明方法流程示意图。

图2为第一模糊图像。

图3为第二模糊图像。

图4为真实的栅格3D形貌。

图5为四种算法深度信息收敛曲线图。

图6为ISTA算法重建的栅格3D形貌。

图7为FISTA算法重建的栅格3D形貌。

图8为MFISTA算法重建的栅格3D形貌。

图9为FL-ISTA算法重建的栅格3D形貌。

图10为ISTA算法重建的栅格3D形貌相对误差曲面。

图11为FISTA算法重建的栅格3D形貌相对误差曲面。

图12为MFISTA算法重建的栅格3D形貌相对误差曲面。

图13为FL-ISTA算法重建的栅格3D形貌相对误差曲面。

图14为四种算法均方差的收敛曲线图。

具体实施方式

为了更好地理解本发明的技术方案，以下结合附图及实施例对本发明作进一步详细说明，但不作为对本发明的限定。

请结合图1，基于改进迭代收缩阈值算法的微观3D重建方法的具体实施是这样的：

步骤一、初始化基本参数，包括摄像机参数焦距f，像距v，凸透镜半径D，模糊度与模糊圆之间的常数λ，第一模糊图像I₁，第一模糊图像I₂，能量阈值ε，能量系数ρ,l，初始迭代步长t₀，初始化深度信息值为S₀。

步骤二、根据已知条件设定相对模糊度Δσ²初始值。第一模糊图像和第二模糊图像为I₁(x,y)和I₂(x,y)，模糊扩散系数分别为σ₁和σ₂，那么相对模糊度表示为Δσ²＝σ₂ ²-σ₁ ²，并且其中，η_i(i＝1,2)表示模糊圆半径；相对模糊度可以由以下公式计算得到

因此，最终的深度信息为

步骤三、计算热扩散方程的解z(x,y,Δt)。图像的模糊化用卷积的形式表示为

I₂(x,y)＝I₁(x,y)*H(x,y) (2)

其中，H(x,y)代表点扩散函数，而“*”代表卷积。根据点扩散原理，点扩散函数可由二维高斯函数来近似表示，由于热辐射方程具有各向同性，所以式(2)又可以表示为

其中，“c”表示扩散系数， 是拉普拉斯算子

如果深度图是理想平面，那么“c”是一个常数，否则此时可以表示为如下形式

其中和分别代表梯度算子和微分算子

两张不同程度的模糊图像I₁(x,y)和I₂(x,y)，同时利用所求的相对模糊度Δσ²，写成如下形式

所以公式(4)可以写成如下形式

因此，由方程组(5)计算得到热扩散方程的解z(x,y,Δt)。

步骤四、根据求得的热扩散方程的解z(x,y,Δt)来计算损失的能量，损失的能量方程表示如下

上式第三项是关于深度图的平滑约束，第四项是保证深度图的有界性，实际上能量系数ρ＞0，l＞0，这是为了尽量少增加能量损失值，s为由相对模糊度带入深度信息表达式计算的深度信息值。如果所求能量值小于能量阈值ε，则停止，此时的深度信息值即为所求。否则，转入步骤五。

步骤五、计算加速算子y_j+1。加速算子由当前点和前一个点的线性组合构成，并且表达如下

y_j+1＝Δσ² _j+a_j(Δσ² _j-Δσ² _j-1)

其中，Δσ² _j,Δσ² _j-1分别代表当前相对模糊度和前一次相对模糊度，Δσ² _j-Δσ² _j-1表示搜索方向，a_j表示由当前相对模糊度Δσ² _j开始，沿着Δσ² _j-Δσ² _j-1所构成的搜索方向进行迭代所需要的步长因子，而y_j+1表示由当前相对模糊度Δσ² _j开始，沿着Δσ² _j-Δσ² _j-1所构成的搜索方向进行步长为a_j所得到的相对模糊度。a_j＝n_j-1/n_j，n₁＝α₁＝1，0＜α_k＜1。

步骤六、将计算得到的加速算子y_j+1带入迭代公式重新进行梯度估计。

对于离焦深度恢复动态优化问题，采用迭代收缩阈值算法进行求解。首先，对于无约束最优化问题

E₁(s)＝∫∫(z(x,y,Δt)-I₂(x,y))²dxdy

解决上式的梯度算法为

其中，t_j表示迭代步长，而上式的二阶估计模型可以表示为

当加入Tikhonov正则项，即

可以得到以下迭代公式

忽略常数项，可以得到

由于Tikhonov正则项是可分的，因此Δσ² _j的计算可以变成解决每个元素的一维问题，它的迭代公式为

其中，T_ρ:Rⁿ→Rⁿ是收缩算子，定义为

将计算得到的加速算子带入式(6)，可以得到

步骤七、利用割线线性搜索求出最优迭代步长t^*，带入迭代公式求出此时的Δσ² _j+1。为了求解最优迭代步长，建立关于步长的一元函数

t^*＝argminψ(t)

其中ψ(t)是关于t的一元函数，t是迭代步长，d是迭代方向，t^*是最优迭代步长，因此可以采用一维线性搜索方法来求极小值

采用割线迭代格式表示为

将利用割线线性搜索计算得到的最优迭代步长t^*带入式(7)，此时可以写成以下形式

步骤八、将步骤七计算所得的Δσ² _j+1带入深度信息表达式(1)，求解出深度信息值s₂。

使用迭代收缩阈值算法(ISTA)、快速迭代收缩阈值算法(FISTA)、单调快速迭代收缩阈值算法(MFISTA)和FL-ISTA算法对标准500nm尺度栅格进行深度信息重建，标准纳米尺度栅格高500nm、宽1500nm。

如图2和图3所示为标准栅格的两幅离焦图像，由于使用摄像机参数固定的离焦深度恢复方法用于微观3D重建时，摄像机参数都是不允许变化的，通过改变物距来获得两张深度信息变化的离焦图像，图2为深度变化前的离焦图像，图3为深度变化后的离焦图像。

如图4所示为标准栅格的真实3D形貌，高500nm，宽1500nm。

如图5所示为四种算法所得深度信息值收敛曲线图。由图5可以看出，在迭代初期，四种算法都能以较快的速度收敛，属于正常现象。但是由于在离焦深度恢复动态优化问题中，FL-ISTA算法在ISTA算法的基础上，利用当前点和前一个点的线性组合构成加速算子，重新进行梯度估计，加快迭代速度；引入割线性搜索寻找最优迭代步长，提高了算法效率。由60-120次的迭代过程中明显可以看出，FL-ISTA算法得到的深度信息值下降最快，并且趋于收敛时得到的深度信息值相对较小，此时更加接近真实的标准500nm栅格尺度，使得恢复出来的3D形貌尺度更精确。

如图6-9所示分别为ISTA、FISTA、MFISTA和FL-ISTA算法处理标准栅格的3D形貌重建结果图。由图5到图8可以看出，四种算法重建的3D形貌，尽管在边缘处误差稍微大，但重建的3D形貌大致符合标准500nm尺度栅格的整体形貌。从重建的3D形貌明显看出，前三种算法重建的3D形貌中深度信息值相差不大，分别为600nm、590nm、590nm左右，而FL-ISTA算法中深度信息值明显下降很多，深度信息值540nm左右，更接近标准500nm栅格尺度，使得恢复出来的3D形貌更精确。并且，在相同的条件下，四种算法对标准纳米尺度栅格进行深度重建实验，经多次试验结果表明，ISTA、FISTA、MFISTA算法深度重建平均运行时间分别为105s、105s、100s左右，而本文所提FL-ISTA算法的平均运行时间仅为30s左右，一定程度上提高了算法的效率。

如图10-13所示分别为ISTA、FISTA、MFISTA和FL-ISTA算法处理标准栅格的3D形貌相对误差结果图。为了进一步检验本文算法的精度，首先，利用公式(8)算了估计深度与真实深度值s之间的相对误差曲面；然后，利用公式(9)误差的均方差，图14给出200次迭代均方差的收敛曲线图，实验结果由多次试验综合得出；最后，因为已知标准栅格的准确高度为500nm，使用公式(10)计算平均误差。

如图14所示为200次迭代过程中均方差的收敛曲线图。可以看出，ISTA、FISTA和MFISTA算法均方差分别为0.05、0.048、0.045，而FL-ISTA算法均方差为0.041，与ISTA算法相比，均方差下降了18百分点，与FISTA和MFISTA算法相比，均方差也较小，因此，FL-ISTA算法重建的微观3D形貌尺度误差相对稳定。因为对于像细胞等样品的观测与测量，重建的3D形貌尺度与真实的3D形貌尺度误差过大或过小，会导致对样品尺度的估计不准确，进而可能造成操作过程中对样品的破坏，从而造成严重后果。ISTA、FISTA和MFISTA算法重建3D形貌平均误差分别为161nm、160nm、158nm，而FL-ISTA算法平均误差仅为96nm，与ISTA算法相比，平均误差下降了40百分点，与FISTA和MFISTA算法相比，平均误差也较小，FL-ISTA算法重建3D形貌尺度与真实3D形貌尺度相差不大，以便于显微镜下对样品的操作。

Claims (6)

1.一种基于改进迭代收缩阈值算法的微观3D重建方法，其特征在于，依次包括以下步骤：

步骤一、初始化基本参数，所述基本参数包括第一模糊图像和第二模糊图像；

步骤二、设定第一模糊图像和第二模糊图像的相对模糊度的初始值，并将第二模糊图像以第一模糊图像和相对模糊度表示；

步骤三、将第一模糊图像和第二模糊图像以图像的模糊化用卷积的形式表示并改写为热扩散方程组，计算热扩散方程的解；

步骤四、利用步骤三得到的热扩散方程的解计算能量方程，如果能量小于能量阈值，则由步骤二得到的相对模糊度带入深度信息表达式计算深度信息值；否则，转入步骤五；

步骤五、计算加速算子y_j+1，y_j+1表示由当前的相对模糊度开始，沿着Δσ² _j-Δσ² _j-1所构成的搜索方向进行步长为a_j所得到的相对模糊度，Δσ² _j代表当前的相对模糊度，Δσ² _j-1代表前一次的相对模糊度，y_j+1＝Δσ² _j+(n_j-1/n_j+1)(Δσ² _j-Δσ² _j-1)，a_j＝n_j-1/n_j，

步骤六、将加速算子y_j+1带入迭代公式重新进行梯度估计；

步骤七、利用割线线性搜索求出最优迭代步长，带入迭代公式求出此时的相对模糊度；

步骤八、将计算所得的此时的相对模糊度带入深度信息表达式，求解出深度信息值。

2.根据权利要求1所述的基于改进迭代收缩阈值算法的微观3D重建方法，其特征在于，所述深度信息表达式为：

，

其中f为摄像机焦距，v为摄像机像距，D为凸透镜半径，λ为模糊度与模糊圆之间的常数，s₂(x,y)为第二模糊图像的深度信息值，s₁(x,y)为第一模糊图像的深度信息值，Δσ²(x,y)为相对模糊度。

3.根据权利要求1所述的基于改进迭代收缩阈值算法的微观3D重建方法，其特征在于，所述第二模糊图像以第一模糊图像和相对模糊度表示是

其中I₂(x,y)为第二模糊图像，I₁(x,y)为第一模糊图像，Δσ²为相对模糊度；所述将第一模糊图像和第二模糊图像以图像的模糊化用卷积的形式表示并改写为热扩散方程组为：

其中为梯度算子，为微分算子，t为扩散时间。

4.根据权利要求1所述的基于改进迭代收缩阈值算法的微观3D重建方法，其特征在于，所述步骤四中，所述能量方程表示如下

其中I₂(x,y)为第二模糊图像，能量系数ρ＞0，l＞0，z(x,y,Δt)为热扩散方程的解，s为由相对模糊度带入深度信息表达式计算的深度信息值。

5.根据权利要求1所述的基于改进迭代收缩阈值算法的微观3D重建方法，其特征在于，所述步骤六中，所述迭代公式由以下方法确定：对于无约束最优化问题

E₁(s)＝∫∫(z(x,y,Δt)-I₂(x,y))²dxdy

解决上式的梯度算法为

其中，t_j表示迭代步长，而上式的二阶估计模型可以表示为

加入Tikhonov正则项得到以下迭代公式

忽略常数项，得到

Δσ² _j的计算变成解决每个元素的一维问题，它的迭代公式为

其中，T_ρ:Rⁿ→Rⁿ是收缩算子，定义为

将计算得到的加速算子带入迭代公式得到

6.根据权利要求1所述的基于改进迭代收缩阈值算法的微观3D重建方法，其特征在于，所述步骤七中，利用割线线性搜索求出最优迭代步长t^*的方法是：建立关于步长的一元函数

t^*＝argminψ(t)，

其中ψ(t)是关于t的一元函数，t是迭代步长，d是迭代方向，t^*是最优迭代步长，采用割线迭代格式表示为