CN109507792B - 一种基于机器学习的光学显微镜自动聚焦方法 - Google Patents

Abstract

本发明提出一种基于机器学习的光学显微镜自动聚焦方法，属于医学图像处理技术领域。该方法首先将光学显微镜采集并分组后的图片用设计的原始特征和组合特征表示，并将图片与组内最清晰图片的序列差值作为图片的标签，然后采用由回归树构成的随机森林对原始特征和组合特征的重要性进行计算，结合设定的阈值多次迭代筛选出重要性较高的特征，接着用留一法和筛选出的特征将数据划分为训练集和测试集来训练梯度提升回归树，最后迭代训练得到的强回归器进行自动聚焦。

Description

一种基于机器学习的光学显微镜自动聚焦方法

技术领域

本发明属于医学图像处理技术领域，涉及一种基于机器学习的光学显微镜自动聚焦方法。

技术背景

传统的人工阅片方法给病理医生带来了繁重的劳动，而且长时间高注意力地阅片易产生视觉疲劳，大大增加了误诊的概率。近年来，随着显微镜的自动化、智能化发展，自动阅片技术开始出现并迅速发展。显微镜的自动阅片技术由自动聚焦算法拍摄清晰的镜下图像，然后进行后续病理分析，显微镜自动聚焦算法作为自动阅片技术的第一步，极大的影响了后续的病理分析过程，其速度和精度显得尤为重要。

显微视觉存在视场小、景深短、聚焦精度要求高等特点，这使得显微镜自动聚焦面临着很多亟待解决的难题。基于图像处理的自动聚焦法由于具有速度快、精度高、成本低和体积小等优点，成为了现代自动聚焦技术发展的主流。现有的基于图像处理的自动聚焦方法可以分为离焦深度法(DFD)和聚焦深度法(DFF)两类：离焦深度法(DFD)首先要对光学显微镜的成像系统进行建模，在此基础上采集模糊的图像，利用模糊图像和模型推导出该图像与焦平面的偏移量，然后驱动步进电机补偿偏移量达到聚焦平面，这种方法主要依据建模的精度，只需要采集2-3幅图像即可完成聚焦，聚焦时间比较短，但是，该方法依赖于成像系统的数学模型，不同的成像系统数学模型差异较大，数学模型建立的不准确则会导致聚焦误差，没有很好的泛化性能。聚焦深度法(DFF)是一种建立在搜索算法上的自动对焦方法，如专利文献CN105652429A公开了一种基于机器学习的显微镜细胞载玻片扫描自动聚焦方法，首先在扫描路径规划阶段，根据细胞载玻片情况确定选择扫描路径；然后在学习阶段，得到变步爬山法的步长、门限值以及低灰度值统计法的阈值；在聚焦阶段，首先使用灰度零值法判断当前位置以确定使用哪种聚焦策略，如果当前位置处于距离焦点较远处时，使用灰度零值法并使用大步长，经过试探判断焦点方向后进行移动，当位置处于焦点附近时选用较小步长，使用低灰度值统计法，经过试探判断焦点方向后进行移动，最后到达焦点位置聚焦完成。聚焦深度法(DFF)主要涉及聚焦评价函数的选择与设计、聚焦窗口的选择以及聚焦搜索策略的制定。在聚焦评价函数和聚焦窗口确定的情况下，聚焦搜索策略对聚焦性能的优劣起着决定性的作用。大量的诸如爬山法、二分法、斐波纳契搜索法、模糊控制搜索法、自适应步长法、函数曲线拟合法、离散差分方程预测法等搜索策略被用来自动聚焦，在一定程度上改善了调焦的速度和精度，但存在通用性较差，速率较慢的缺点，不适用于精度较高的显微镜自动聚焦。

发明内容

本发明针对现有技术的不足，提供一种基于机器学习的光学显微镜自动聚焦方法，将自动聚焦的调焦过程视为回归问题，通过设计特征和构造回归器，来同时获取显微镜调焦过程需要移动的步长和方向，避免了调焦过程中反复移动带来的回程误差。在保证显微镜聚焦精度的同时，大大的提高了聚焦的速度。

本发明的具体技术方案如下：

一种基于机器学习的光学显微镜自动聚焦方法，所述方法包括以下步骤：

步骤1：沿着光学显微镜的Z轴采集p张原始图片Iⁱ(c，r)，并转化为灰度图fⁱ(c，r)，其中i＝{1，2，...，p}，p＝100～300；

步骤2：重复步骤1，每次记为1组，采集q组数据共N张图片，q＝20～40；

步骤3：计算每张图片的48个原始特征

和96个组合特征

将所有图片用特征的形式表示，记为D＝{(x₁，y₁)，(x₂，y₂)，...，(x_N，y_N)}，其中x表示当前图片的特征，y表示当前图片到聚焦最清晰图片的序列差值，原始特征和组合特征的总数记为m，

和

分别表示显微镜当前位置和下一个位置图片的第j种原始特征，

表示第j种原始特征的第k 种组合特征，其中j＝{1,2, ...,24},k＝{1,2,3,4}，m＝144；

步骤4：采用自助法(bootstrapping)将所有数据集D分为n个数据子集{D₁，D₂，...，D_n}，其中未抽到的样本构成袋外数据集

步骤5：对于每个数据子集D_i，从所有特征中随机抽取w个特征(w＜m)做为回归树T_i的训练集，训练生成回归树T_i，其中i＝{1，2，...，n}，w＝1～143；

步骤6：重复步骤5，训练得到n棵回归树，用袋外数据计算均方误差，记为

其中i＝{1，2，...，n}；

步骤7：随机对

所有所述当前图片的特征x的原始特征j加入噪声干扰，再次用袋外数据计算均方误差，记为

其中i＝{1，2，...，n}，j＝{1，2，…，n}；

步骤8:计算原始特征j的重要性I(j)，其中

j＝{1，2，...，m}；

步骤9：重复步骤7-8，计算出所有特征的重要性，筛选出r个重要性较高的特征，记为 S＝{s¹,s², ...,s^r}；

步骤10:用步骤9中的特征S和留一法将数据集D划分为训练集D_train和测试集D_test；

步骤11:在训练集D_train上采用基于梯度提升的回归树(Gradient BoostingRegression Tree)方法训练得到强回归树F_M(x)，其中

F₀表示初始的弱回归树，v表示学习率，M表示迭代的总次数，m＝{1，2，...，M}，i＝1，2，...，n， j＝1，2，...，J_m，J_m表示迭代第m次时回归树的叶节点个数，γ_jm表示第m次迭代时第j个叶结点的最佳拟合值，I(x∈R_jm)表示第m次迭代当样本归结为第j个叶结点时值为1，否则为0；

步骤12：采用迭代强回归树F_M(x)进行聚焦，直到不符合迭代条件。

更进一步地，步骤3中，采用聚焦评估函数F^j作为图片原始特征，采用显微镜当前视野图片的聚焦评估函数计算值

和下一个位置的计算值

作为图片的原始特征，

聚焦评估函数F^j可以为以下24种：

原始特征F¹为autocorrelation，计算方法为

其中μ，σ²是灰度图片f(c，ｒ)均值和方差，M，N表示为灰度图片f(c，r)长宽的像素点数量,k是常数，这里设置为2；

原始特征F²为Brenner，计算方法为

原始特征F³为entropy histogram，计算方法为

其中p(k)表示图片中像素值为k的像素点的相对频率；

原始特征F⁴为firstDerivGaussian，计算方法为

其中⊙表示相关运算，

为7*7的滤波器，

其中

i，j为整数，取值为

这里设置σ为0.8；

原始特征F⁵为3×3Laplacian，计算方法为

其中⊙表示相关运算，

原始特征F⁶为5×5Laplacian，计算方法为

其中

原始特征F⁷为LoG，计算方法为

其中

为9*9的滤波器，产生的方式如下：

其中i，j为整数，取值范围如下

这里σ设置为1.2；

原始特征F⁸为M&G histogram计算方法为

其中h(k)表示图片中像素值为k的数量，μ为θ_MG+1，这里θ_MG为阈值，计算方法如下

Δ＝2(f(c，r-1)-f(c，r+1))²+2(f(c-1，r)-f(c+1，r))²

+(f(_c-1，r-1)-f(c+1，r+1))²+(f(c-1，r+1)-f(c+1，r-1))²

原始特征F⁹为M&M histogram，计算方法为

其中h(k)表示图片中像素值为k的像素点的数量，μ表示图片像素值的均值，这里设置为 128；

原始特征F¹⁰为normalized variance，计算方法为

其中μ是灰度图片f(c，r)均值；

原始特征F¹¹为3×3Prewitt，计算方法为

其中

原始特征F¹²为range histogram，计算方法为

max{h(k)}-min{h(k)}

其中h(k)表示图片中像素值为k的像素点的数量；

原始特征F¹³为3×3Roberts，计算方法为

其中

原始特征F¹⁴为3×3Scharr，计算方法为

其中

原始特征F¹⁵为3×3Sobel，计算方法为

其中

原始特征F¹⁶为3×3crossSobel，计算方法为

其中

原始特征F¹⁷为5×5Sobel，计算方法为

其中

原始特征F¹⁸为5×5crossSobel，计算方法为

其中

原始特征F¹⁹为squared gradient，计算方法为

原始特征F²⁰为threshold gradient，计算方法为

|f(c，r+1)-f(c，r)|≥threshold

只有f(c，r+1)-f(c，r)的绝对值大于等于设定的阈值才累加，这里阈值为0；

原始特征F²¹为threshold content，计算方法为

f(c，r)≥threshold

只有f(c，r)大于等于设定的阈值才累加，这里阈值为0；

原始特征F²²为variance，计算方法为

其中μ是灰度图片f(c，r)均值；

原始特征F²³为Vollath’s F4，计算方法为

原始特征F²⁴为Vollath’s F5，计算方法为

其中μ表示图片像素值的均值；

采用如下4种组合方式作为图片的组合特征：

其中j＝1，2，...，24；

其中j＝1，2，...，24；

其中j＝1，2，...，24；

其中j＝1，2，...，24。

更进一步地，步骤3中，采用当前图片到聚焦最清晰图片的序列差值作为训练集的标签。

更进一步地，步骤9中，采用随机森林计算特征重要性，结合设定的阈值迭代多次进行特征的筛选。

更进一步地，步骤10中，采用随机森林计算出来的特征重要性把数据划分训练集D_train和测试集D_test。

更进一步地，步骤11中，采用梯度提升回归树的方法对显微镜下一步要移动的步长和方向同时做出预测，采用回归的方法，拟合出来的结果正负表示显微镜下一步需要移动的方向，数值大小表示显微镜下一步需要移动的步长。

本发明能够得到以下有益效果：

本发明将自动聚焦的调焦过程视为回归问题，通过设计特征和构造回归器，来同时获取显微镜调焦过程需要移动的步长和方向，避免了调焦过程中反复移动带来的回程误差。在保证显微镜聚焦精度的同时，大大的提高了聚焦的速度。在自己采集的数据集中，设置初始位置为每组第一张图片时，平均只需要移动7.6次；精度达到了100％，任意初始化位置时最佳精度达到了99.9％，平均只需要移动6.9次。

附图说明

图1是本发明的流程图；

图2是采集的一组数据集中的部分图片以及其聚焦评估函数的曲线；

图3-图26是24种原始特征在20组数据集上的计算结果，每组中节选出100张图片；

图27是随机森林构建的过程；

图28是梯度提升树的迭代过程。

具体实施方式

为了使本发明的发明目的、技术方案和有益技术效果更加清晰，以下结合实施例，对本发明进行进一步详细说明。应当理解的是，本说明书中描述的实施例仅仅是为了解释本发明，并非为了限定本发明，实施例的具体参数设置等可因地制宜做出选择而对结果并无实质性影响。

步骤1：如图2所示，沿着光学显微镜的Z轴采集200张原始图片Iⁱ(c，r)，并转化为灰度图fⁱ(c，r)，其中i＝{1，2，...，p}；

步骤2：重复步骤1，每次记为1组(如图2所示)，采集20组数据共N张图片；

步骤3：计算每张图片的48个原始特征

和96个组合特征

和

分别表示显微镜当前位置和下一个位置图片的第j种原始特征，

表示第j种原始特征的第k种组合特征，原始特征和组合特征的具体计算方法如下：

原始特征F¹为autocorrelation(如图3所示)，横轴表示图片编号，纵轴表示每张图片对应的特征值，每种线型代表一组图片，计算方法为

其中μ，σ²是灰度图片f(c，ｒ)均值和方差，M，N表示为灰度图片f(c，r)长宽的像素点数量,k是常数，这里设置为2；

原始特征F²为Brenner(如图4所示)，计算方法为

原始特征F³为entropy histogram(如图5所示)，计算方法为

其中p(k)表示图片中像素值为k的像素点的相对频率；

原始特征F⁴为firstDerivGaussian(如图6所示)，计算方法为

其中⊙表示相关运算，

为7*7的滤波器，

其中

i，j为整数，取值为

这里设置σ为0.8；

原始特征F⁵为3×3Laplacian(如图7所示)，计算方法为

其中⊙表示相关运算，

原始特征F⁶为5×5Laplacian(如图8所示)，计算方法为

其中

原始特征F⁷为LoG(如图9所示)，计算方法为

其中

为9*9的滤波器，产生的方式如下：

其中i，j为整数，取值范围如下

这里σ设置为1.2；

原始特征F⁸为M&G histogram(如图10所示)计算方法为

其中h(k)表示图片中像素值为k的数量，μ为θ_MG+1，这里θ_MG为阈值，计算方法如下

Δ＝2(f(c，r-1)-f(c，r+1))²+2(f(c-1，r)-f(c+1，r))²

+(f(c-1，r-1)-f(c+1，r+1))²+(f(c-1，r+1)-f(c+1，r-1))²

原始特征F⁹为M&Mhistogram(如图11所示)，计算方法为

其中h((k)表示图片中像素值为k的像素点的数量，μ表示图片像素值的均值，这里设置为128；

原始特征F¹⁰为normalized variance(如图12所示)，计算方法为

其中μ是灰度图片f(c，r)均值；

原始特征F¹¹为3×3Prewitt(如图13所示)，计算方法为

其中

原始特征F¹²为range histogram(如图14所示)，计算方法为

max{h(k)}-min{h(k)}

其中h(k)表示图片中像素值为k的像素点的数量；

原始特征F¹³为3×3Roberts(如图15所示)，计算方法为

其中

原始特征F¹⁴为3×3Scharr(如图16所示)，计算方法为

其中

原始特征F¹⁵为3×3Sobel(如图17所示)，计算方法为

其中

原始特征F¹⁶为3×3crossSobel(如图18所示)，计算方法为

其中

原始特征F¹⁷为5×5Sobel(如图19所示)，计算方法为

其中

原始特征F¹⁸为5×5crossSobel(如图20所示)，计算方法为

其中

原始特征F¹⁹为squared gradient(如图21所示)，计算方法为

原始特征F²⁰为threshold gradient(如图22所示)，计算方法为

|f(c，r+1)-f(c，r)|≥threshold

只有f(c，r+1)-f(c，r)的绝对值大于等于设定的阈值才累加，这里阈值为0；

原始特征F²¹为threshold content(如图23所示)，计算方法为

f(c，r)≥threshold

只有f(c，r)大于等于设定的阈值才累加，这里阈值为0；

原始特征F²²为variance(如图24所示)，计算方法为

其中μ是灰度图片f(c，r)均值；

原始特征F²³为Vollath’s F4(如图25所示)，计算方法为

原始特征F²⁴为Vollath’s F5(如图26所示)，计算方法为

其中μ表示图片像素值的均值；

组合特征

其中j＝1，2，...，24；

组合特征

其中j＝1，2，...，24；

组合特征

其中j＝1，2，...，24；

组合特征

其中j＝1，2，...，24；

将所有图片用特征的形式表示，记为D＝{(x₁，y₁)，(x₂，y₂)，...，(x_N，y_N)}，其中x表示当前图片的特征，y表示当前图片到聚焦最清晰图片的序列差值(如图2所示聚焦评估的曲线中x轴的差值)；

步骤4：如图27所示，采用自助法(bootstrapping)将所有数据集D分为n个数据子集{D₁，D₂，…，D_n}，其中未抽到的样本构成袋外数据集

步骤5：如图27所示，对于每个数据子集D_i，从样本中随机抽取w个特征(w＜m)做为回归树T_i的训练集，训练生成回归树T_i，每棵决策树生长不受限制，也不要做剪裁，其中i＝{1，2，...，n}，m为基础特征和组合特征的总和，这里设置为144；

步骤6：重复步骤5，训练得到n棵回归树，用袋外数据计算均方误差，记为

其中i＝{1，2，...，n}；

步骤7：随机对

所有当前图片的特征x的原始特征j加入噪声干扰，再次用袋外数据计算均方误差，记为

其中i＝1{1，2，...，n}；

步骤8:计算原始特征j的重要性I(j)，计算方法如下

步骤9：重复步骤7-8，计算出所有特征的重要性，每一次将I(j)大于阈值的原始特征j加入集合S当中，最终集合S即为随机森林特征选择后的结果，记为S＝{s¹，s²，...，s^k}，sⁱ表示选择后的特征，这里阈值设置为5％；

步骤10:用步骤9中的特征S和留一法将数据集D划分为训练集D_train和测试集D_test，这里留一法表示将19组数据做为训练集D_train，余下的一组数据做为测试集D_test；

步骤11:如图28所示，用D_train＝{(x₁，y₁)，(x₁，y₁)，...，(x_n，y_n)}做为训练集，采用基于梯度提升回归树(Gradient Boosting Regression Tree)方法训练得到强回归树F_M(x)，具体过程如下

1)使用常量值γ初始化模型，

这里γ设置为

2)迭代训练模型，对于第m次迭代，其中m＝{1，2，...，M}

a)计算损失函数的负梯度在当前模型的值，将它作为残差的估计

其中L(y_i，F(x_i))为损失函数，这里设置为最小二乘(Least-square)，i＝1，2，...，n；

b)将{(x_i，r_im)}作为训练集，其中i＝1，2，...，n，拟合得到第m棵回归树，叶节点区域记为R_jm,j＝1，2，...，J_m，其中J_m为迭代第m次时叶子节点的个数；

c)对叶子区域j＝1，2，...,J_m计算最佳拟合值

d)如图28所示更新模型

其中v是学习率，I(x∈R_jm)表示第m次迭代当样本归结为第j个叶结点时值为1，否则为0；

3)更新模型,得到强学习器

上面对本专利的较佳实施方式作了详细说明，但是本专利并不限于上述实施方式，在本领域的普通技术人员所具备的知识范围内，还可以在不脱离本专利宗旨的前提下做出各种变化。

Claims (7)

1.一种基于机器学习的光学显微镜自动聚焦方法，其特征在于，包括如下步骤：

步骤1：沿着光学显微镜的Z轴采集p张原始图片Iⁱ(c，r)，并转化为灰度图fⁱ(c，r)，其中i＝{1，2，...，p}，p＝100～300；

步骤2：重复步骤1，每次记为1组，采集q组数据共N张图片，q＝20～40；

步骤3：计算每张图片的特征，特征数记为m，将所有图片用特征的形式表示，记为D＝{(x₁，y₁)，(x₂，y₂)，...，(x_N，y_N)}，其中x表示当前图片的特征，y表示当前图片到聚焦最清晰图片的序列差值；

步骤4：采用自助法将所有数据集D分为n个数据子集{D₁，D₂，...，D_n}，其中未抽到的样本构成袋外数据集

步骤5：对于每个数据子集D_i，从所有特征中随机抽取w个特征做为回归树T_i的训练集，训练生成回归树T_i，其中i＝{1，2，...，n}，w＜m，w＝1～143；

步骤6：重复步骤5，训练得到n棵回归树，用袋外数据计算均方误差，记为

其中i＝{1，2，...，n}；

步骤7：随机对

所有所述当前图片的特征x的原始特征j加入噪声干扰，再次用袋外数据计算均方误差，记为

其中i＝{1，2，...，n}，j＝{1，2，...，m}；

步骤8:计算原始特征j的重要性I(j)，其中

步骤9：重复步骤7-8，计算出所有特征的重要性，结合阈值筛，多次迭代选出重要性较高的特征；

步骤10:用步骤9中的特征和留一法将数据集D划分为训练集D_train和测试集D_test，这里留一法表示留一组数据做为测试集D_test，其余组别全部作为训练集D_train；

步骤11:在训练集D_train上采用基于梯度提升的回归树方法训练得到强回归树F_M(x)，其中

F₀表示初始的弱回归树，v表示学习率，M表示迭代的总次数，J_m表示迭代第m次时回归树的叶节点个数，γ_jm表示第m次迭代时第j个叶结点的最佳拟合值，I(x∈R_jm)表示第m次迭代当样本归结为第j个叶结点时值为1，否则为0；

步骤12：采用迭代强回归树F_M(x)进行聚焦，直到不符合迭代条件。

2.根据权利要求1所述的方法，其特征在于，步骤3中，图片的特征采用原始特征聚焦评估函数

聚焦评估函数

和组合特征

和

分别表示显微镜当前位置和下一个位置图片的第j种原始特征，

表示第j种原始特征的第k种组合特征，j＝{1,2,...,24}，k＝{1,2,3,4}。

3.根据权利要求2所述的方法，其特征在于，步骤3中，采用如下4种组合方式作为图片的组合特征：

其中j＝1，2，...，24；

其中j＝1，2，...，24；

其中j＝1，2，...，24；

其中j＝1，2，...，24。

4.根据权利要求1所述的方法，其特征在于，步骤3中，采用当前图片到聚焦最清晰图片的序列差值作为训练集的标签。

5.根据权利要求1所述的方法，其特征在于，步骤9中，采用随机森林计算特征重要性，结合设定的阈值迭代多次进行特征的筛选。

6.根据权利要求1所述的方法，其特征在于，步骤10中，采用随机森林计算出来的特征重要性把数据划分训练集D_train和测试集D_test。

7.根据权利要求1所述的方法，其特征在于，步骤11中，采用梯度提升回归树的方法对显微镜下一步要移动的步长和方向同时做出预测，采用回归的方法，拟合出来的结果正负表示显微镜下一步需要移动的方向，数值大小表示显微镜下一步需要移动的步长。

Images