CN110517213B - 一种显微镜用的基于拉普拉斯金字塔的实时景深延拓方法 - Google Patents

Abstract

本发明公开了一种显微镜用的基于拉普拉斯金字塔的实时景深延拓方法，包括：筛选模块，根据景深变化判据筛选输入图像；配准模块，平衡亮度并估计图像间的形变，进行修正；融合模块，实时输入图像和前一次融合图像分别建立拉普拉斯金字塔，采用最大对比度原则融合金字塔各层，最后还原融合后的金字塔为最终融合图像。本发明针对实时处理需求，提出了一套完善的景深延拓方案，可用于实时处理显微观测图像，并获得大景深、高质量的景深延拓图像。借助上述方法，用户只需手动旋转细调焦旋钮，无需关心旋转速度的快慢，就可自动得到一幅超景深显微图像。上述方法对于景深变化具有较强鲁棒性，不易受景深变化程度或方向的改变而影响最终融合质量。

Description

一种显微镜用的基于拉普拉斯金字塔的实时景深延拓方法

技术领域

本发明涉及图像处理技术领域，具体涉及一种显微镜用的基于拉普拉斯金字塔的实时景深延拓方法。

背景技术

景深延拓技术是指将同一个观测对象，在不同聚焦层下拍摄的图像融合成一幅大景深图像的技术，该技术在显微数码成像领域有非常重大的意义。

目前获取大景深图像的主要方法是通过步进电机控制显微镜的Z轴调节旋钮，将不同层的图像成像到图像传感器，在定位以后，一张张保存下来，然后对所有不同层的图像进行后继的独立融合，最后得到一幅大景深的图像。

例如，公开号为CN 108319010A的专利说明书公开了一种超景深图像采集系统，采用一个含有能电动升降调景深的显微镜，所述的显微镜含有底座、立柱、调节主机升降的齿轮齿条传动机构和装在立柱上的步进电机，载物台及光源。此系统还包括装在显微镜底盘内的控制器和外置的计算机。其工作流程是；在载物台上放置待观察采集的物体，显微镜从上至下在计算机通过控制器控制步进电机等距运行，每停顿一次计算机通过CCD摄像头采集一幅局部清晰图像，然后计算机将采集的N幅局部清晰图像融合，则得到一幅突破景深限制的从上至下都清晰的图像。

另一种方法是通过步进电机精确控制Z轴调焦机械，对普通显微镜Z轴进行等步距扫描，边扫描边融合。这种方法也需要对传统的显微镜进行改装，对专业应用有一定的意义，但对于普通用户，会增加成本，意义不大。

除此以外，常见的景深延拓算法用于实时处理有很多困难。原因可以归结于以下几点，一是配准困难，实时处理时不能轻易的确定用于配准的参考图像，而且在焦深跨度大的情况下，并不一定能找到足够可信的特征点对；二是融合质量问题，局部融合产生的图像撕裂，或者连续融合时的不断加权都会影响最后的融合图像质量，这些运算需要大量的计算开支，也加剧了实时处理的难度。

公开号为CN 108020509A的专利说明书公开了一种光学投影层析成像的方法，拍摄生物样品若干张不同角度不同焦面的图像；利用拉普拉斯金字塔算法将同一个角度不同焦面的图像融合成长景深图片；利用滤波反投影算法将生物样品不同角度的长景深图片重建。该技术方案仍无法解决配准问题。

发明内容

针对本领域存在的不足之处，本发明提供了一种显微镜用的基于拉普拉斯金字塔的实时景深延拓方法，可实现显微镜图像的实时景深延拓。

一种显微镜用的基于拉普拉斯金字塔的实时景深延拓方法，重复进行实时图像处理流程，直至完成实时景深延拓；

所述实时图像处理流程包括将实时图像依次传入筛选模块、配准模块和融合模块，得到融合图像和输出图像；

所述筛选模块用于判断当前实时图像与前一帧传入配准模块的实时图像的景深是否发生变化，若发生变化，则将当前实时图像作为输入图像传入所述配准模块，并作为参考图像用于对比下一帧实时图像的景深是否发生变化；若未发生变化，则直接进入下一帧实时图像的判断和筛选；

所述配准模块用于平衡输入图像和前一次实时图像处理流程所得的融合图像的亮度，再通过建立输入图像和融合图像的图像金字塔，在不同尺度上计算两幅图像的形变关系，并将融合图像根据形变关系进行配准；

所述融合模块根据输入图像和配准后的融合图像建立拉普拉斯金字塔，并根据最大对比度原则将两组拉普拉斯金字塔逐层融合得到一组融合金字塔，融合金字塔以带有锐化的方式还原，得到输出图像，以未锐化的方式还原，得到用于下一次融合的融合图像。

本发明方法包括以下三大模块：筛选模块，根据景深变化判据筛选输入图像；配准模块，平衡亮度并估计图像间的形变，进行修正；融合模块，实时输入图像和前一次融合图像分别建立拉普拉斯金字塔，采用最大对比度原则融合金字塔各层，最后还原融合后的金字塔为最终融合图像。

所述筛选模块从实时图像中筛选出图像进入融合队列，基于尽量让更多位于不同焦面深度的图像参与景深延拓且尽可能节约运算资源这两个前提，本发明将这个过程分为静止状态和焦面深度改变两种状态。后者由于焦面深度改变，图像原本成像清晰的部分变得模糊，原本不在景深范围内的部分变得清晰，而前者像素的改变主要来自于噪声。

为了区分焦面深度变化和噪声引起的图像变化，作为优选，所述筛选模块判断景深是否发生变化的方法包括步骤：

(1-A)将实时图像和参考图像灰度化后相减取绝对值获得差异图像；

(1-B)将差异图像按式(I)进行阈值化处理：

其中，T为阈值，为差异图像的像素均值和像素方差的线性组合；

(1-C)阈值化处理后的差异图像经压缩得到8*8大小的压缩差异图像，所述压缩差异图像的像素值是原差异图像对应区域的像素平均值，若压缩差异图像的像素值均不小于1，判定实时图像和参考图像的景深发生了变化；否则，判定实时图像和参考图像的景深未发生变化。

进一步优选，所述阈值T按式(II)计算：

T＝ε+nσ (II)，

其中，ε为差异图像的像素均值，σ为差异图像的像素方差，n为常数。经测试，在一优选例中，n为3。

作为优选，所述筛选模块实时感应显微镜的Z轴状态，若显微镜的Z轴处于静止状态，所述筛选模块停止将未发生景深变化的实时图像传入配准模块，直至实时图像的景深发生变化。

焦面深度的改变会影响图像的亮度，为了减小亮度变化对景深延拓结果的影响，需要平衡输入图像和前一次融合得到的融合图像的整体亮度。此外，本发明将输入图像作为参考图像，把前一次融合得到的融合图像向参考图像进行配准。

作为优选，所述配准模块平衡输入图像和前一次实时图像处理流程所得的融合图像的亮度的方法包括步骤：

(2-A)对输入图像和融合图像进行像素采样，并分别累加每幅图像的采样像素的R值、G值和B值；

(2-B)将融合图像的各通道的采样像素之和与输入图像的同一通道的采样像素之和的比值作为平衡输入图像亮度的系数，所述系数与输入图像的同一通道的像素值相乘，完成输入图像的亮度平衡。

在一优选例中，步骤(2-A)中的像素采样大小为：水平方向和竖直方向的步长分别为被采样图像的宽和高的百分之一。

步骤(2-B)中，以R通道为例，平衡输入图像亮度的计算公式为：

其中，

为输入图像平衡亮度后的各像素R值，

为输入图像平衡亮度前的各像素R值，

为输入图像平衡亮度前的采样像素的R值之和，

为融合图像的采样像素的R值之和。G通道、B通道的计算方法与R通道相同。

作为优选，所述配准模块进行配准的方法包括步骤：建立输入图像和融合图像的图像金字塔，遵循由粗及细的原则，从尺度最小的一层开始计算形变系数，使得图像差值D为最小值，并将所得形变系数作为初始值传递到相邻尺度更大的下一层，如此循环，直至获得原始分辨率下的形变系数，并将融合图像根据所得形变系数进行配准；

所述图像差值D按式(III)计算：

D＝∑(I₁(x,y)-I₂(x,y；p))² (III)，

其中，I₁为输入图像，I₂为融合图像，p为形变系数。

进一步优选，所述形变系数以平移量t、缩放量s和旋转量r描述；

除尺度最小的一层外，其它层的形变系数由相邻尺度较小一层的形变系数计算得到，按式(IV)～(VI)根据m+1层的形变系数计算第m层形变系数：

t_m＝2×t_m+1 (IV)，

r_m＝r_m+1 (V)，

s_m＝s_m+1 (VI)；

其中，t_m为第m层的平移量，t_m+1为第m+1层的平移量，r_m为第m层的旋转量，r_m+1为第m+1层的旋转量，s_m为第m层的缩放量，s_m+1为第m+1层的缩放量。

作为优选，所述根据最大对比度原则将两组拉普拉斯金字塔逐层融合为一组金字塔的计算公式如(VII)所示：

其中，

为融合图像拉普拉斯金字塔的第k层中坐标为(x,y)的像素值，

为输入图像拉普拉斯金字塔的第k层中坐标为(x,y)的像素值，P_O(x,y；k)为融合金字塔的第k层中坐标为(x,y)的像素值。

为了获得更好的显示效果，本发明在输出融合结果前会对结果进行适当锐化。但经过了锐化的图像在下一次融合过程中再次锐化会加剧噪声的影响，所以本发明仍保留未锐化的融合结果参与下一次融合。

以锐化的方式还原融合金字塔，是指对尺度较大的金字塔层添加增益系数，以起到增强高频信号的作用，作为优选，所述融合金字塔以带有锐化的方式还原的计算公式如式(VIII)～(X)所示：

P₀₁(k)＝f(k)×P_o(k)+b(k)×P_o1(k+1) (VIII)，

其中，P_o(k)是融合金字塔的第k层，P_o1(k+1)是被还原的融合金字塔第k+1层在每两行、每两列之间插值后再与建立拉普拉斯金字塔时使用的高斯滤波器卷积得到的和P_o(k)大小相同的图像，P_o1(k)是最后得到的被还原的输出金字塔第k层，f(k)为前景系数，b(k)为背景系数，K为还原后的输出金字塔的总层数，fg、bg为常数。在一优选例中，fg为1.25，bg为0.93。

本发明与现有技术相比，主要优点包括：本发明针对实时处理需求，提出了一套完善的景深延拓方案，可以用于实时处理显微观测图像，并获得大景深、高质量的景深延拓图像。借助本发明方法，用户只需手动旋转细调焦旋钮，无需关心旋转速度的快慢，就可自动得到一幅超景深显微图像。本发明方法对于景深变化具有较强鲁棒性，不易受景深变化程度或方向的改变而影响最终融合质量。

附图说明

图1为本发明的显微镜用的基于拉普拉斯金字塔的实时景深延拓方法的融合流程示意图；

图2为实施例的筛选模块的判断和筛选流程示意图；

图3为实施例的配准模块的配准流程示意图；

图4为实施例的融合模块的融合流程示意图。

具体实施方式

下面结合附图及具体实施例，进一步阐述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。下列实施例中未注明具体条件的操作方法，通常按照常规条件，或按照制造厂商所建议的条件。

本实施例的显微镜用的基于拉普拉斯金字塔的实时景深延拓方法进行一次融合的流程如图1所示，在任务开始时，系统首先会读入一帧实时图像，作为筛选模块和配准模块的参考图像。随后的实时图像首先会进入筛选模块，与参考图像对比，区分像素变化是由噪声还是焦面变化引起的。如果未检测到焦面变化引起的像素变化，实时图像被舍弃，系统继续等待比较下一张图像；如果焦面发生变化，实时图像会被加入到融合队列，等待后续步骤的处理，同时也将更新筛选模块中的参考图像为当前实时图像。

图2是筛选模块的一个具体实例，输入的实时图像和参照图像会转化为灰度图进行处理。由于随机噪声可以以高斯噪声描述，将差异图像像素均值和方差的线性组合作为阈值能够有效区分噪声和焦面深度变化产生的像素变化，这一过程的具体步骤如下：

(1-A)将实时图像转化为灰度图，和参考图像的灰度图相减并取绝对值，获得差异图像；

(1-B)计算差异图像中像素均值和方差，根据均值和方差按下式计算阈值，将差异图像二值化，阈值T是差异图像像素均值和方差的线性组合：

T＝ε+nσ

其中，ε为差异图像的像素均值，σ为差异图像的像素方差，n为3。

将二值化的差异图像压缩到8*8大小，压缩差异图像的每个像素值是原差异图像对应区域的平均值。当压缩差异图像的像素不全小于1时，就认为两幅图像之间景深发生了变化，将实时图像传入下一个步骤，同时当前图像也会用作下一次筛选过程的参考图像；否则，继续下一次筛选过程。

随后图像会被传入配准模块，以实时输入图像I₁作为参考图像，对前一次的融合结果I₂进行配准。在配准之前，I₁会根据I₂平衡亮度，这一过程的具体步骤如下：

(2-A)以一定的步长，在R通道、G通道和B通道分别统计I₁和I₂的像素和，本实施例中水平方向和竖直方向的步长选为图像宽和高的百分之一。

(2-B)根据I₁和I₂在R通道、G通道和B通道抽样像素和的比值，平衡I₁的整体亮度。

然后再以I₁作为参考图像，对I₂进行配准。配准的方式利用了图像金字塔，从低分辨率图像开始计算形变系数，进行由粗及细的计算方法，以加速配准的过程。在本实施例中，形变系数以平移量T、缩放量S和旋转量R描述。

图3是经过平衡亮度后配准模块的一个具体实例，处理图像的分辨率为4096*3080，具体步骤如下：

(3-A)对I₁和I₂建立图像金字塔，直至尺度最小的一层宽或高小于32，本实施例中图像金字塔分为8层，尺度最小的一层尺寸是32*24，出于显示效果原因，图3只表示了前4层；

(3-B)在融合图像I₂金字塔尺度最小的一层上施加一定的形变p，再计算其与作为参照图像的实时图像I₁金字塔对应层的图像差值D，不断改变形变p中的平移量t、缩放量s和旋转量r，直到图像差值D达到极小值；

图像差值D按下式计算：D＝∑(I₁(x,y)-I₂(x,y；p))²；

在较小尺度层获得的形变系数p会用于计算尺度较大一层的初始形变系数，从第m+1层得到的形变系数推算第m层的初始型变系数推算的方式如下：

t_m＝2×t_m+1，

r_m＝r_m+1，

s_m＝s_m+1；

(3-C)重复步骤(3-B)，直到在原始分辨率下得到形变系数，对上一次融合结果I₂进行修正，然后将配准后的两幅图像传输给下一步。

最后图像会被传递到融合模块，在融合模块中，实时输入图像I₁和经过配准的前一次融合图像I₂通过建立拉普拉斯金字塔来实现融合，融合金字塔可以以加强高频信息的方式还原并输出，或者以不失真的方式还原参与下一次的融合过程。

图4是拉普拉斯金字塔融合模块中，对于一个通道处理的实例，出于显示原因，拉普拉斯金字塔也只显示了4层，其具体步骤如下：

(4-A)对I₁和I₂建立拉普拉斯金字塔，直至尺度最小的一层宽或高小于32，本例中拉普拉斯金字塔分8层；

(4-B)以最大对比度原则融合金字塔，以拉普拉斯金字塔第k层各像素为例，最大对比度原则可以表述为如下公式：

(4-C)在需要实时显示融合结果，或者输出融合图像时，融合金字塔需要以锐化的方式还原，以第k层拉普拉斯金字塔的还原过程为例，具体计算方式为：

P_o1(k)＝f(k)×P_o(k)+b(k)×P_o1(k+1)，

其中，P_o1(k)为还原后的输出金字塔的第k层，P_o(k)为融合金字塔的第k层，P_o1(k+1)是从还原后的第k+1层金字塔，通过插值和卷积，放大至和P_o(k)大小相同的图像，f(k)为锐化时添加的前景系数，b(k)为锐化时添加的背景系数，f(k)和b(k)可以由当前的金字塔层与总层数的关系计算得到，K为还原后的输出金字塔的总层数，fg为1.25，bg为0.93。

(4-D)用于继续融合的图像需要以不失真的方式从融合金字塔还原，这一方式中f(k)和b(k)都为常数1。

至此完成一次实时图像处理，重复进行上述实时图像处理流程，直至任务结束，完成一次完整的实时景深延拓过程。

本实施例方法可应用于与显微相机配套的软件，用于实时景深延拓时，将实时图像流传输到软件的输入接口，手动调整显微镜物距，即可生成景深实时延拓图像流。

此外应理解，在阅读了本发明的上述描述内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。

Claims (9)

1.一种显微镜用的基于拉普拉斯金字塔的实时景深延拓方法，其特征在于，重复进行实时图像处理流程，直至完成实时景深延拓；

所述实时图像处理流程包括将实时图像依次传入筛选模块、配准模块和融合模块，得到融合图像和输出图像；

所述筛选模块用于判断当前实时图像与前一帧传入配准模块的实时图像的景深是否发生变化，若发生变化，则将当前实时图像作为输入图像传入所述配准模块，并作为参考图像用于对比下一帧实时图像的景深是否发生变化；若未发生变化，则直接进入下一帧实时图像的判断和筛选；

所述配准模块用于平衡输入图像和前一次实时图像处理流程所得的融合图像的亮度，再通过建立输入图像和融合图像的图像金字塔，在不同尺度上计算两幅图像的形变关系，并将融合图像根据形变关系进行配准；

所述融合模块根据输入图像和配准后的融合图像建立拉普拉斯金字塔，并根据最大对比度原则将两组拉普拉斯金字塔逐层融合得到一组融合金字塔，融合金字塔以带有锐化的方式还原，得到输出图像，以未锐化的方式还原，得到用于下一次融合的融合图像。

2.根据权利要求1所述的显微镜用的基于拉普拉斯金字塔的实时景深延拓方法，其特征在于，所述筛选模块判断景深是否发生变化的方法包括步骤：

(1-A)将实时图像和参考图像灰度化后相减取绝对值获得差异图像；

(1-B)将差异图像按式(I)进行阈值化处理：

其中，T为阈值，为差异图像的像素均值和像素方差的线性组合；

(1-C)阈值化处理后的差异图像经压缩得到8*8大小的压缩差异图像，所述压缩差异图像的像素值是原差异图像对应区域的像素平均值，若压缩差异图像的像素值均不小于1，判定实时图像和参考图像的景深发生了变化；否则，判定实时图像和参考图像的景深未发生变化。

3.根据权利要求2所述的显微镜用的基于拉普拉斯金字塔的实时景深延拓方法，其特征在于，所述阈值T按式(II)计算：

T＝ε+nσ (II)，

其中，ε为差异图像的像素均值，σ为差异图像的像素方差，n为常数。

4.根据权利要求1所述的显微镜用的基于拉普拉斯金字塔的实时景深延拓方法，其特征在于，所述筛选模块实时感应显微镜的Z轴状态，若显微镜的Z轴处于静止状态，所述筛选模块停止将未发生景深变化的实时图像传入配准模块，直至实时图像的景深发生变化。

5.根据权利要求1所述的显微镜用的基于拉普拉斯金字塔的实时景深延拓方法，其特征在于，所述配准模块平衡输入图像和前一次实时图像处理流程所得的融合图像的亮度的方法包括步骤：

(2-A)对输入图像和融合图像进行像素采样，并分别累加每幅图像的采样像素的R值、G值和B值；

(2-B)将融合图像的各通道的采样像素之和与输入图像的同一通道的采样像素之和的比值作为平衡输入图像亮度的系数，所述系数与输入图像的同一通道的像素值相乘，完成输入图像的亮度平衡。

6.根据权利要求1所述的显微镜用的基于拉普拉斯金字塔的实时景深延拓方法，其特征在于，所述配准模块进行配准的方法包括步骤：建立输入图像和融合图像的图像金字塔，遵循由粗及细的原则，从尺度最小的一层开始计算形变系数，使得图像差值D为最小值，并将所得形变系数作为初始值传递到相邻尺度更大的下一层，如此循环，直至获得原始分辨率下的形变系数，并将融合图像根据所得形变系数进行配准；

所述图像差值D按式(III)计算：

D＝∑(I₁(x，y)-I₂(x，y；p))² (III)，

其中，I₁为输入图像，I₂为融合图像，p为形变系数。

7.根据权利要求6所述的显微镜用的基于拉普拉斯金字塔的实时景深延拓方法，其特征在于，所述形变系数以平移量t、缩放量s和旋转量r描述；

除尺度最小的一层外，其它层的形变系数由相邻尺度较小一层的形变系数计算得到，按式(IV)～(VI)根据m+1层的形变系数计算第m层形变系数：

t_m＝2×t_m+1 (IV)，

r_m＝r_m+1 (V)，

s_m＝s_m+1 (VI)；

其中，t_m为第m层的平移量，t_m+1为第m+1层的平移量，r_m为第m层的旋转量，r_m+1为第m+1层的旋转量，s_m为第m层的缩放量，s_m+1为第m+1层的缩放量。

8.根据权利要求1所述的显微镜用的基于拉普拉斯金字塔的实时景深延拓方法，其特征在于，所述根据最大对比度原则将两组拉普拉斯金字塔逐层融合为一组金字塔的计算公式如(VII)所示：

其中，

为融合图像拉普拉斯金字塔的第k层中坐标为(x，y)的像素值，

为输入图像拉普拉斯金字塔的第k层中坐标为(x，y)的像素值，P_O(x，y；k)为融合金字塔的第k层中坐标为(x，y)的像素值。

9.根据权利要求1所述的显微镜用的基于拉普拉斯金字塔的实时景深延拓方法，其特征在于，所述融合金字塔以带有锐化的方式还原的计算公式如式(VIII)～(X)所示：

P_O1(k)＝f(k)×P_O(k)+b(k)×P_O1(k+1) (VIII)，

其中，P_O(k)是融合金字塔的第k层，P_O1(k+1)是被还原的融合金字塔第k+1层在每两行、每两列之间插值后再与建立拉普拉斯金字塔时使用的高斯滤波器卷积得到的和P_O(k)大小相同的图像，P_O1(k)是最后得到的被还原的输出金字塔第k层，f(k)为前景系数，b(k)为背景系数，K为还原后的输出金字塔的总层数，fg、bg为常数。

Images