CN100351057C - 基于显微图像处理的微操作工具深度信息提取方法及装置 - Google Patents

Abstract

本发明涉及一种微操作工具深度信息的提取，特别是基于显微图像处理的，用于细胞、染色体等微细生物体的深度信息的提取方法及装置，属于微操作机器人技术领域。以往通过增加额外传感装置来获得深度信息，开销大影响实际应用。本发明的技术方案：根据显微镜成象模型，通过显微图象处理，以提取显微镜点扩散参数σ；再通过神经网络建立的点扩散参数与离焦量的关系，来获取操作工具的深度信息。提取装置包括显微镜、微针、CCD摄像机、微操作机械手、计算机等。本发明的有益效果在40X显微镜条件下，深度提取标准差在2微米以内，最大有效范围是显微镜聚焦带上下各70微米；计算一个深度值的时间为0.2秒左右，在精度和效率两方面均达到了在线应用的水平。

Description

基于显微图像处理的微操作工具深度信息提取方法及装置

技术领域

本发明涉及一种微操作工具深度信息的提取，特别是基于显微图像处理的，用于生物医学实验，对象是细胞、染色体等微细生物体的深度信息的提取方法及装置，属于微操作机器人技术领域。

背景技术

微操作机器人是机器人技术向微细操作领域的延伸。机器人的基本工作方式是基于笛卡儿坐标进行操作[卢桂章，张建勋，赵新.面向生物工程实验的微操作机器人.南开大学学报，1999，32(3)：42-46]。利用微操作机器人系统进行生物医学实验，无论是驱动操作工具，还是到目标点的轨迹规划，都要求高精度；因此，操作过程中微操作工具的精确定位十分关键。微操作系统中作为视觉反馈的显微图像是一个二维图像，通过显微图像处理获取X，Y方向坐标的方法有很多；它的难点在于如何获得Z方向坐标，即深度信息。

获得深度信息的方法有两种：增加额外的传感装置和处理现有的显微图像。一种额外的传感装置是，垂直于Z-X平面或Y～Z平面增加另外的显微镜镜头和CCD摄像机系统，使得Z方向坐标可以象X，Y方向坐标一样，直接通过显微图像处理得到。另一种方式是增加微传感器，主要用于微机电系统(MEMS)的装配。面向生物医学工程的微操作机器人系统，操作对象是微细的生物体(如细胞、染色体)；操作工具是由玻璃管拉制而成的玻璃针，它的末端直径在1微米～几十微米；这就导致在操作工具末端加装传感器来获取反馈信号非常困难。此外，加装垂直于Z-X平面或Y-Z平面显微镜头必然会加大整个系统的开销，并有可能影响实际应用。

通过图像处理获取深度信息是一个受到广泛关注的研究方向[Pemtland，A.A New Sensefor Depth of Field.IEEE Trans.on Pattern Analysis and Machine Intelligence.1987，9(4)：523-531][Grossmann，P.Depth from Focus.Pattern Recognition Letters，1987，5(1)：63-69]，比较有代表性的方法是Pentland等人提出的，利用离焦状态图像信息，计算摄像机扩散参数σ，来恢复对象深度信息[2 Pemtland，A.A New Sense for Depth of Field.IEEE Trans.on PatternAnalysis and Machine Intelligence.1987，9(4)：523-531]。此方法核心之处在于：通过象素灰度梯度计算摄像机扩散参数σ。这种方法较好地解决了上述问题，最终得到的精度“可以与立体测量相比”(在18米范围内，最大误差40厘米)。此方法唯一的限制是“确保足够的高频信息来获得图像之间的变化量”。

随着微操作机器人系统研究更加广泛和深入的发展，通过显微图像处理获取深度信息受到越来越广泛的关注。我们前期的实验已经证明，显微图像具有不同于宏观图像的特点，象素灰度梯度由于对噪声敏感，不能正确反映显微图像的离焦程度[谢少荣，微系统的虚拟环境研究，天津大学博士论文，2001]，文献[Pemtland，A.A New Sense for Depth of Field.IEEETrans.on Pattern Analysis and Machine Intelligence.1987，9(4)：523-531]通过象素灰度梯度计算摄像机扩散参数σ也就无法实现。文献[Pemtland，A.A New Sense for Depth of Field.IEEETrans.on Pattern Analysis and Machine Intelligence.1987，9(4)：523-531]中，要求足够的高频信息，指的是图像内容的足够的高频信息，也就是通常说的图像内容“边缘清楚”；而对显微图像，特别是生物医学应用的显微图像，很难满足“边缘清楚”这一条件。试想，在宏观图像中，景深可以有几米，那么，至少在0.1毫米以下的宽度才能称为边缘；以生物细胞注射常用的40X显微镜为例，景深只有2微米，至少在纳米级才能成为边缘。因此，在生物医学应用的显微图像中，高频信息主要包括图像噪声，图像内容的高频信息是不存在的；通过象素灰度梯度计算摄像机扩散参数σ就无法实现。

文献[张建勋，薛大庆，卢桂章，李彬。通过显微图像特征抽取获得微操作目标纵向信息。机器人，2001，23(1)：73-77]给出了一种基于自动调焦技术实现深度信息恢复的方法。基本过程如下：1)利用Fourier变换，给出显微图像的模糊度判据；2)利用模糊度判据，获得聚焦状态最清晰的显微图像，对应聚焦状态最清晰图像的，被观察物体当前Z方向位置，定义为Z方向坐标原点；3)被观察物体沿Z方向移动后，自动记忆沿Z方向移动相对变化量，换算成Z方向坐标值。因为显微图像中高频分量是噪声，而图像内容存在于低频部分；所以，上述方法与宏观图像不同处理在于：Fourier变换后作低通滤波处理，滤掉高频分量。

文献[张建勋，薛大庆，卢桂章，李彬。通过显微图像特征抽取获得微操作目标纵向信息。机器人，2001，23(1)：73-77]利用聚焦状态的最清晰显微图像作为定位基准建立深度信息恢复方法，存在两方面问题：1)自动调焦技术是来源于光学技术，对于光学系统，最终目的是成像，景深大意味着能在一个较大的区域成清晰的像。而微操作机器人系统需要定位信息，Z方向坐标获取基于显微镜聚焦区域给出的，则意味着定位精度不会高于景深范围。如40X显微镜系统，景深2微米，实际精度在4-5微米。提高精度的主要途径选用景深更小的显微镜，这显然是不可取的。2)被观察物体处于离焦状态，坐标是由自动记忆的Z方向位移量换算得到的；因此，定位精度依赖于Z方向位移精度，这样离焦状态Z方向定位精度进一步降低。

同时，一些研究者从研究显微镜成像角度入手，设计了各种特定的光源条件，对显微镜点扩散参数σ的获取及性质进行了分析，做了很多很好的工作[谢少荣，微系统的虚拟环境研究，天津大学博士论文，2001]；其目的在于研究显微镜成像的基本性质。但由于其主要手段是制作各种各样的点光源或等效点光源的叠加，利用点物体成像模型来辨识点扩散参数；此方法用来研究成像系统性质是足够的，但无法取得任意给定显微图像的扩散参数σ，也就无法获取与其对应的深度信息。“利用显微图像提取微小物体深度信息”这一问题，尚未得到很好解决。

发明内容

总结上述方法，基于自动调焦技术的深度信息恢复方法，仅仅利用了聚焦状态的显微图像特征恢复深度信息，而大范围的离焦状态显微图像信息并未加以利用。尽管文献[Pemtland，A.A New Sense for Depth of Field.IEEE Trans.on Pattern Analysis andMachine Intelligence.1987，9(4)：523-531]提出的通过象素灰度梯度计算摄像机扩散参数σ的方法对显微图像不适用，但通过计算摄像机扩散参数σ，获得深度信息依然是一条很好的思路。因为扩散参数σ是一个光学系统的固有参数，不会因光照等条件而改变；对于一个给定的光学系统，扩散参数σ与离焦量之间的关系是固定的。因此，本发明致力于建立恰当的方法，通过显微图像处理，获得对应的显微镜扩散参数σ；并通过实验建立起其与Z方向坐标的关系；最终，可以通过显微图像处理获得对象深度信息。将此方法实际应用于微操作过程中，提取微操作工具——微针的深度信息，取得了较好的效果。

本发明的目的是提供一种基于图像处理的微操作工具深度信息提取方法，该方法将被观察物体的深度信息提取问题转化为显微镜点扩散参数的提取问题。首先，根据显微镜成像模型，通过显微图像处理，在线提取显微镜点扩散参数；之后，通过神经网络建立点扩散参数与离焦量(Z方向深度信息)的映射关系；然后，通过上述映射关系，得到微操作工具的深度信息；最终，将此方法应用于微操作机器人系统，利用离焦状态双针互插实验验证了方法的有效性。

本发明的技术方案：

这种基于显微图像处理的微操作工具深度信息提取装置，其特点在于：倒置显微镜光源发出光经过聚光镜进入分光棱镜，将光束分成相互垂直的两束；其中一束光透过分光棱镜照亮位于XOY平面上的微针，微针通过倒置显微镜物镜成像于装在显微镜上的电荷耦合器件(CCD)摄像机，经过图像采集接口将不同离焦状态的微针图像送入计算机；另一束光经过分光棱镜反射到位于YOZ平面的纵向标尺，标尺的局部图像通过另一个显微镜物镜成像于显微摄像装置的CCD摄像机，读取实际微针离焦量；微操作机械手是由计算机控制的。

这种基于显微图像处理的微操作工具深度信息提取方法，要经过下述步骤：根据显微镜成像模型，通过显微图像处理，提取显微镜点扩散参数σ，然后依据神经网络建立起来的点扩散参数α与离焦量的关系来获得操作工具的深度信息。

下面对上述方案作一些说明：

1.显微镜成像模型

如附图1所示，根据几何光学成像原理，如果实际成像平面与聚焦成像平面有一定距离，那么物点P在像平面上所成的像不再是一个清晰的像点，而是一个模糊的圆，即焦斑。若透镜和成像平面的位置都固定不变，只改变物点的位置，则该模糊圆的直径d与物距之间的关系为：

$d\left(\mathrm{\ϵ}\right)=\frac{\mathrm{\ϵ}\·f\·D}{(u_0+\mathrm{\ϵ})(u_0-f)}---\left(1\right)$

式中α₀为成清晰像时的物距；ε为物距u₀的偏移量，又称之为离焦量，设定远离透镜方向为正，反之为负；f为透镜的焦距；D为透镜的直径。对于给定显微系统，f、D、u₀都是固定值，经标定可以获得。

由于光的衍射效应和透镜成像的非理想化，焦斑中的光强度分布可用二维高斯函数来近似：

$h(x,y)=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}{e}^{-\frac{x^2+y^2}{2{\mathrm{\σ}}^2}}---\left(2\right)$

式中h(x，y)意义为点物体的像，称为该成像系统的点扩散函数(point spread function)，σ是点扩散参数[Pemtland，A.A New Sense for Depth of Field.IEEE Trans.on Pattern Analysisand Machine Intelligence.1987，9(4)：523-531]。

卷积能很好地描述显微镜系统的成像[赵新，余斌等，基于系统辨识的显微镜点扩散参数提取方法及应用。计算机学报，2004，27(1)：140-144]。若给定的被观察物体f(u，v)，经一光学系统所成的像为g(x，y)，则有

$g(x,y)=f*h={\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}f(u,v)h(x-u,y-v)\mathrm{dudv}---\left(3\right)$

式中*为卷积操作，h(x，y)为系统的点扩散函数，它决定被观察物体离焦后的扩散程度，也就是所成显微图像g(x，y)的“模糊”程度。h(x，y)由式(2)给出，变量x，y是像空间位置信息，点扩散参数σ唯一确定被观察物体离焦后扩散程度，即显微图像的“模糊”程度。

2.显微镜点扩散参数辨识

对于一个光学成像系统，求解点扩散参数σ通常利用式(2)进行。由于得到一个理想的点物体很困难，考虑到显微图像的特殊性，选择利用式(3)的面物体成像模型来辨识σ。将式(2)代入式(3)，得到：

$g(x,y)=\frac{1}{2{\mathrm{\π\σ}}^2}{\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}f(u,v)e^{-\frac{{(x-u)}^2+{(y-v)}^2}{2{\mathrm{\σ}}^2}}\mathrm{dudv}---\left(4\right)$

给定被观察物体f(x，y)和点扩散参数σ，利用式(4)可以求出被观察物体所成的像，像与点扩散参数σ一一对应(为示区别，将计算所得的像称为虚拟像，将实际采得图像称为物体的像)。因此，当被观察物体f(x，y)的像已知时，通过与虚拟像对比，可以确定它对应的点扩散参数。基于此，给定一幅待确定σ的被观察物体的像，可以按下面步骤确定显微镜点扩散参数：

1)对于每个可能的σ，利用式(4)生成被观察物体的虚拟像；

2)将虚拟像和物体的像进行对比，求得对应点差值的平方和：

$\underset{x,y}{\mathrm{\Σ}}{(g(x,y)-\frac{1}{2{\mathrm{\π\σ}}^2}{\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}f(u,v)e^{-\frac{{(x-u)}^2+{(y-v)}^2}{2{\mathrm{\σ}}^2}}\mathrm{dudv})}^2---\left(5\right)$

3)搜索σ，寻找令虚拟像和物体的像匹配最佳的

$J=\mathrm{Min}\left\{\underset{x,y}{\mathrm{\Σ}}{(g(x,y)-\frac{1}{2{\mathrm{\π}\widehat{\mathrm{\σ}}}^2}{\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}f(u,v)e^{-\frac{{(x-u)}^2+{(y-v)}^2}{2{\widehat{\mathrm{\σ}}}^2}}\mathrm{dudv})}^2\right\}---\left(6\right)$

4)据最小二乘法原理，上一步找到的虚拟像就是物体像的“等价结果”。此时的

就是生成物体像时成像系统的点扩散参数。

3.点扩散参数的快速提取方法

直接采用上述方法求解扩散参数，需要进行大量二维卷积运算，计算效率很低。考虑到计算效率，假设采用的被观察物体是平行于X轴的无限长狭缝。由于其在X方向成像结果相同，考虑其成像时，只需计算平行于Y方向一条线上各象素点灰度值即可，我们将选取的这条线称为扫描线。同样原因，计算扫描线上各象素点灰度值时，只需在被观察物体与之对应的扫描线上作卷积即可，这样就将成像模型简化为一维的形式，则式(4)可表示为：

$g(x,y)=\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\·{\∫}_{y_{\min }}^{y_{\max }}f(x,v)e^{-\frac{{(y-v)}^2}{2{\mathrm{\σ}}^2}}\mathrm{dv}---\left(7\right)$

其中x为确定值。相应地，式(5)、(6)分别表示为：

$\underset{y}{\mathrm{\Σ}}(g(x,y)-\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\·{\∫}_{y_{\min }}^{y_{\max }}f(x,y)e^{-\frac{{(y-v)}^2}{2{\mathrm{\σ}}^2}}\mathrm{dv})---8$

$J=\mathrm{Min}\left\{\underset{y}{\mathrm{\Σ}}(g(x,y)-\frac{1}{\sqrt{2\mathrm{\π}}\widehat{\mathrm{\σ}}}\·{\∫}_{y_{\min }}^{y_{\max }}f(x,y)e^{-\frac{{(y-v)}^2}{2{\widehat{\mathrm{\σ}}}^2}}\mathrm{dv})\right\}---\left(9\right)$

实际上，只要满足被观察物体扫描线的扩散范围不超过X方向两端，则像空间中扫描线灰度分布，与对应的理想无限长狭缝任意扫描线的灰度分布相同。对于微操作工具来说，微针尖端的宽度只有1微米左右，而长度较长，因此只要合理选择扫描线位置，完全可以满足上述假设。如此，在计算扩散参数时可以只考虑一条垂直扫描线上的数据。

由于微针本身是具有一定尺度的物体，当它在显微镜下成像时，各部分并不处在同一个平面上，也就是说，图像各部分的离焦程度，也即扩散参数并不相同。由于最终目的是获得清晰的针尖，因此，扫描线要尽量选取在微针尖端，取在针尖尖端类似矩形的中部效果最好，一方面满足了上述假设，同时该处点扩散参数也较为接近，利于减小误差。因此，扫描线的位置与针尖位置是相关的，本发明首先通过图像处理的方法，获得微针的针尖位置，再通过针尖位置选择合理的扫描线位置。

另外，选取扫描线时，还存在一个“对准”的问题。即要求从每一幅离焦图像中选出来的垂直扫描线，必须与从被观察物体的清晰像中选出的那条扫描线直接对应，前者是后者生成的像。这要求在成像过程中，针尖的位置始终不变。本发明通过图像模式匹配的方法较好地解决了针尖的对准问题，从而使得选取的扫描线可以很好地支持算法的处理。

4.微针边界提取

在搜索点扩散参数的过程中，需要计算虚拟像和采样样本之间逐点灰度差的平方和，该灰度差平方和的大小决定了虚拟像和采样样本之间匹配程度的好坏。在计算过程中，虚拟像的大小和采样样本是一致的，但是求平方和时，不能将整条扫描线都计算在内，因为微针成像的关键信息都包含在微针像中，而不是外边的背景区域。因此，计算时，必须忽略背景区域的大部分点，而只考虑微针所成的像。这样求得的结果才真正反映了虚拟像和物体像之间匹配程度。

划分显微图像中微针所在的区域，基本上是一个简化了的边界提取问题：不需要准确标记出微针的边界，只需要确保划分的区域完全将微针包括在内即可。经过多次试验，本发明给出了一个基于中值滤波的方法，实现了从扫描线中提取微针所在的区域。

5.点扩散参数的搜索方法

黄金分割搜索算法是搜索点扩散参数时采用的方法。该算法作为一种一维搜索算法，要求被搜索函数在给定搜索区间内具有单峰特性。根据对式(9)的分析可知，当点扩散参数作为自变量时，评估值函数总体上是单峰的。因此，点扩散参数的搜索问题可以由黄金分割搜索算法解出。

另一方面，该算法的搜索效率较高。从分析的结果来看，在[1，100]区间内(该区间是点扩散参数的合理取值范围)进行搜索，当搜索精度为0.005时，算法迭代10次左右即能给出结果；当精度为0.1，只需要迭代5次。实际上，当扩散参数小于1/6时，本像素点对其他像素点的影响可以忽略不计。因此，搜索精度不必设的很高，这样可以进一步提高搜索效率。

6.虚拟图像数据的计算算法

根据式(3)可知，在计算虚拟图像数据时，需要通过积分来实现卷积：首先获取积分限，并将离散的被积函数(被观察物体)通过插值的方法转化为连续形式，最终完成积分。结果表明，通过此种方法，搜索一个点扩散参数大致需要10-20秒。

实际上，由于获得的采样图像本身是离散的，没有必要使用连续卷积，这并不会提高虚拟数据的精度，却在很大程度上增加了计算量。因此，本发明使用离散卷积来代替连续卷积，并且根据傅立叶变换的性质，在计算过程中，采用“傅立叶变换-乘积-逆傅立叶变换”的方法来实现离散卷积。经过上述改进，在其他条件不变的情况下，搜索一个点扩散参数仅需0.1秒。

7.通过神经网络建立点扩散参数与离焦量的关系

本发明借助当前研究非线性系统的主要工具——人工神经网络来建立点扩散参数与离焦量的映射关系。

BP网络是应用最广的一类人工神经网络。它是一个多层的、具有非线性可导传递函数的神经网络。神经元是神经网络最基本的组成部分，一般地，一个有R个输入的神经元模型如附图13所示。其中P为输入向量，w为权向量，b为阈值，f为传递函数，a为神经元输出。所有输入P通过权重w进行加权求和后加上阈值b再经传递函数f的作用后即为该神经元的输出a，即：

$a=f(\underset{i}{\mathrm{\Σ}}{P}_i{w}_i+b)---\left(10\right)$

传递函数可以是任何可微的函数，常用的有Sigmoid型传递函数与线性传递函数。附图14是一个两级BP网络拓扑结构，由输入层、输出层以及一个隐含层互连而成。训练神经网络，即将网络的输入和输出反复作用于网络，不断调整其权重和阈值，以使网络误差函数值达到最小，从而实现输入、输出间的非线性映射。

理论上已经证明，一个带阈值的、隐含层采用S型传递函数，输出层采用线性传递函数的二级BP网络可以逼近任何具有有限个不连续点的函数。实际中采用的就是上述二级BP网络，从结果观察，此种网络可以很好的描述点扩散参数与深度信息之间的关系。

本发明的有益效果：本发明提出并实现了一种基于显微图像处理的微操作工具深度信息提取方法。此方法在40X显微镜条件下，深度提取标准差在2微米以内，最大有效范围是显微镜聚焦带上下各70微米；在普通微机上，计算一个深度值的时间为0.2秒左右，在精度和效率两方面均达到了在线应用的水平。

进一步，将上述方法应用于微操作机器人系统，设计并实现了双针互插实验，验证了方法的有效性，同时也将微操作机器人的操作空间进行了拓展，使得操作在聚焦状态和离焦状态下均可进行。

附图说明

图1.显微镜成像模型

图2.深度信息提取处理程序模块划分

图3.深度信息提取处理流程

图4.点扩散参数辩识程序流程图

图5.微针显微图像及标尺图像

图6.点扩散参数拟合结果

图7.点扩散参数拟合标准差

图8.微针点扩散参数与离焦量的关系

图9.深度信息提取结果验证

图10.微操作机器人系统双针互插实验

图11.显微物镜扩散函数测试装置结构示意图

图12.多组标准狭缝板局部放大图

图13.神经元模型

图14.二级BP神经网络拓扑结构

其中：1.倒置显微镜光源中的灯泡  2.显微镜光源中的聚光镜  3.分光棱镜  4.微操作工具(微针)  5.显微镜上的物镜  6.显微镜上的电荷耦合器件(CCD)摄像机  7.纵向标尺  8.采集纵向位置信息的显微物镜  9.采集纵向位置信息的电荷耦合器件(CCD)10.微操作机械手  11.计算机

具体实施方式

下面结合附图对本发明作具体说明：

这种基于显微图像处理的微操作工具深度信息提取装置，其特点在于：倒置显微镜光源1发出光经过聚光镜2进入分光棱镜3，将光束分成相互垂直的两束；其中一束光透过分光棱镜照亮位于XOY平面上的微针4，微针通过倒置显微镜物镜5成像于装在显微镜上的CCD摄像机6，经过图像采集接口将不同离焦状态的微针图像送入计算机11；另一束光经过分光棱镜反射到位于YOZ平面的纵向标尺7，标尺的局部图像通过另一个显微镜物镜8成像于显微摄像装置的CCD摄像机9，读取实际微针离焦量；微操作机械手10是由计算机11控制的。

这种基于显微图像处理的微操作工具深度信息提取方法，其特点在于要经过下述步骤：

1)获取采样图像及标尺图像

2)读取并处理清晰图像，获得扫描线位置，保存需要数据

3)按照顺序，依次读取并处理采样图像，保存需要数据

4)求取采样图像扫描线上的微针边界

5)计算采样图像的点扩散参数

6)判断所有采样图像是否已处理完，是则继续，否则转3)

7)读取所有采样图像的标尺信息，获得离焦量

8)通过神经网络建立点扩散参数与离焦的映射关系，保留结果

9)获取需要进行深度信息提取的模糊图像

10)读取并处理该图像，保存需要数据

11)求取模糊图像扫描线上的微针边界

12)计算模糊图像的点扩散参数

13)利用神经网络结果，求取模糊图像中物体的离焦量

上面所述的提取方法，步骤5)包括下列步骤：

1)随机设置图像的点扩散参数

2)根据成像模型，计算虚拟图像数据

3)计算评估值

4)根据评估值变化，重新获取点扩散参数

5)判断参数精度是否达到要求，是则继续，否则转2)

6)保留使评估值较小的扩散参数

7)随机设置y方向上偏移量

8)加上偏移量，重新计算虚拟图像数据

9)计算评估值

10)根据评估值变化，重新获取偏移量

11)判断精度是否达到要求，是则继续，否则转8)

12)保留使评估值较小的偏移量

13)判断偏移量是否改变，是则继续，否则保留结果，搜索结束

14)和偏移前比较，判断评估值是否减小，是则继续，否则保留结果搜索结束

15)保留偏移量，重新搜索点扩散参数，转1)

显微镜成像模型，由于光的衍射效应和透镜成像的非理想化，成像系统的点扩散函数可用二维高斯函数来近似：

$h(x,y)=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}{e}^{-\frac{x^2+y^2}{2{\mathrm{\σ}}^2}}---\left(2\right)$

式中h(x，y)意义为点物体的像，称为点扩散函数，σ是点扩散参数；

卷积能很好地描述显微镜系统的成像[6]。若给定的被观察物体f(u，v)，经一光学系统所成的像为g(x，y)，则有

$g(x,y)=f*h={\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}f(u,v)h(x-u,y-v)\mathrm{dudv}---\left(3\right)$

式中*为卷积操作，h(x，y)为系统的点扩散函数，它决定被观察物体离焦后的扩散程度，也就是所成显微图像g(x，y)的“模糊”程度。h(x，y)由式(2)给出，变量x，y是像空间位置信息，点扩散参数σ唯一确定被观察物体离焦后扩散程度，即显微图像的“模糊”程度。

点扩散参数σ通常利用式(2)进行求解，由于得到一个理想的点物体很困难，考虑到显微图像的特殊性，选择利用式(3)的面物体成像模型来辨识σ。将式(2)代入式(3)，得到：

$g(x,y)=\frac{1}{2{\mathrm{\π\σ}}^2}{\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}f(u,v)e^{-\frac{{(x-u)}^2+{(y-v)}^2}{2{\mathrm{\σ}}^2}}\mathrm{dudv}---\left(4\right)$

给定被观察物体f(x，y)和点扩散参数σ，利用式(4)可以求出被观察物体所成的像，像与点扩散参数σ一一对应，为示区别，将计算所得的像称为虚拟像，将实际采得图像称为物体的像：因此，当被观察物体f(x，y)的像已知时，通过与虚拟像对比，可以确定它对应的点扩散参数；基于此，给定一幅待确定σ的被观察物体的像，可以按下面步骤确定显微镜点扩散参数；

1)对于每个可能的σ，利用式(4)生成被观察物体的虚拟像；

2)将虚拟像和物体的像进行对比，求得对应点差值的平方和：

$\underset{x,y}{\mathrm{\Σ}}{(g(x,y)-\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}{\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}f(u,v)e^{-\frac{{(x-u)}^2+{(y-v)}^2}{2{\mathrm{\σ}}^2}}\mathrm{dudv})}^2---\left(5\right)$

3)搜索σ，寻找令虚拟像和物体的像匹配最佳的

$J=\mathrm{Min}\left\{\underset{x,y}{\mathrm{\Σ}}{(g(x,y)-\frac{1}{2\mathrm{\π}{\widehat{\mathrm{\σ}}}^2}{\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}f(u,v)e^{-\frac{{(x-u)}^2+{(y-v)}^2}{2{\widehat{\mathrm{\σ}}}^2}}\mathrm{dudv})}^2\right\}---\left(6\right)$

4)据最小二乘法原理，上一步找到的虚拟像就是物体像的“等价结果”，此时的

就是生成物体像时成像系统的点扩散参数。

如果微操作工具是细长型微针，那么点扩散参数σ可以通过以下步骤快速确定：

1)对于每个可能的σ，利用式(7)生成被观察物体的虚拟像；

$g(x,y)=\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\·{\∫}_{y_{\min }}^{y_{\max }}f(x,v)e^{\frac{{(y-v)}^2}{2{\mathrm{\σ}}^2}}\mathrm{dv}---\left(7\right)$

2)将虚拟像和物体的像进行对比，求得对应点差值的平方和：

$\underset{y}{\mathrm{\Σ}}(g(x,y)-\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\·{\∫}_{y_{\min }}^{y_{\max }}f(x,v)e^{\frac{{(y-v)}^2}{2{\mathrm{\σ}}^2}}\mathrm{dv})---\left(8\right)$

3)搜索σ，寻找令虚拟像和物体的像匹配最佳的

$J=\mathrm{Min}\left\{\underset{y}{\mathrm{\Σ}}(g(x,y)-\frac{1}{\sqrt{2\mathrm{\π}}\widehat{\mathrm{\σ}}}\·{\∫}_{y_{\min }}^{y_{\max }}f(x,v)e^{\frac{{(y-v)}^2}{2{\widehat{\mathrm{\σ}}}^2}}\mathrm{dv})\right\}---\left(9\right)$

4)根据最小二乘法原理，上一步找到的虚拟像就是物体像的“等价结果”，此时的

就是生成物体像时成像系统的点扩散参数σ。

扫描线要选取在微针尖端，取在针尖尖端类似矩形的中部效果最好，并且要求每一幅离焦图像中选出来的垂直扫描线必须与从被观察物体的清晰像中选出的那条扫描线直接“对应”，也就是要求在成像过程中，针尖的位置始终不变。

点扩散参数σ在搜索的过程中，计算虚拟像和采样样本之间逐点灰度差的平方和时，必须忽略背景区域的大部分点，而只考虑微针所成的像，微针成像的关键信息都包含在微针像中。

微针成像所在的区域，不需要准确标记出微针的边界，只需要确保划分的区域完全将微针包括在内即可，经过多次试验，本发明给出了一个基于中值滤波的方法，实现了从扫描线中提取微针所在的区域。

点扩散参数σ搜索，可以采用黄金分割搜索算法，该算法作为一种一维搜索算法，要求被搜索函数在给定搜索区间内具有单峰特性，根据对式(9)的分析可知，当点扩散参数作为自变量时，评估值函数总体上是单峰的，因此，点扩散参数的搜索问题可以由黄金分割搜索算法解出，该算法的搜索效率较高。

在计算虚拟图像数据时，可以使用离散卷积代替连续卷积，在计算过程中，采用“傅立叶变换一乘积-逆傅立叶变换”来实现离散卷积。

点扩散参数σ与离焦量的映射关系可借助当前研究非线性系统的主要工具——人工神经网络BP来建立，本发明采用的是带阈值的、隐含层采用S型传递函数、输出层采用线性传递函数的两层BP神经网络。

实施例

1.微操作工具深度信息提取

深度信息提取的处理程序模块划分及处理流程见附图2、3。下面对具体步骤进行详细描述。

1)采样显微图像获取

采集的显微图像包括微针图像及侧面标尺两部分。微操作系统的显微镜物镜放大倍数是40倍；在微操作机器人的控制下，操作臂移动的运动步长为2微米。采集的图像包括近离焦和远离焦两种，采集的起点和终点选择在微针在图像中不能分辨的时刻。微针近离焦和远离焦的程度在70微米左右，超出这个范围后，微针就和背景融合在一起，难以分辨了。附图5给出的微针图像及对应标尺图像是一组图像数据中比较有代表性的几幅，该组共有100幅显微图像。图下数字为根据标尺读取的离焦量，此处只显示了近离焦的情况，离焦量为负值，远离焦的情况与之类似，离焦量为正值。

2)点扩散参数辨识

根据本发明介绍的点扩散参数的快速提取方法，及针对微针图像所作的特殊处理，可以准确、迅速地求出微针显微图像的点扩散参数。附图4给出了点扩散参数辨识的程序流程图；附图6显示了部分图像的拟合结果，与附图5一一对应。图中由“+”描出的线表示采样图像扫描线上的灰度值(取值范围为[0，255])，由“.”描出的线是根据最优扩散参数σ计算得到的虚拟像。

用同样的方法可以求出所有图像的点扩散参数，附图7给出了各幅图像采样扫描线与虚拟像之间的拟合标准差(灰度取值范围为[0，1])，从图中可以观察到，所有的标准差都在4％以下，拟合结果达到了预定的要求。

3)点扩散参数和离焦量的映射关系建立

根据通过辨识方法求得的点扩散参数，加上从标尺图像读出的离焦量，可以拟合出两者的映射关系。点扩散参数和离焦量的关系分成两段，一段描述微针在近离焦下的变化趋势，另一段描述其在远离焦下的变化趋势。

本发明通过二层BP神经网络获取点扩散参数与离焦量的关系，BP网络的基本情况如下：网络输入层只有一个节点，是图像的点扩散参数；隐含层节点的个数可配置，实际选取的是30个节点，采用的是S型传递函数；输出层也只有一个节点，采用线性传递函数，输出的结果就是微针的离焦量。下表分别列举了近离焦与远离焦情况下，30个隐含层节点的权重与域值。

	近离焦权重	近离焦域值	远离焦权重	远离焦域值
					1	-0.980000	98.000000	-0.980000	98.000000
2	0.980000	-95.117647	-0.980000	95.117647
					3	0.980000	-92.235294	-0.980000	92.235294
4	0.980000	-89.352941	-0.980000	89.352941
					5	-0.980000	86.470588	-0.980000	86.470588
6	0.980000	-83.588235	-0.980000	83.588235
					7	0.980000	-80.705882	0.980000	-80.705882
8	0.980000	-77.823529	0.980000	-77.823529
					9	-0.980000	74.941176	0.980000	-74.941176
10	-0.980000	72.058824	0.980000	-72.058824

11	-0.980000	69.176471	-0.980000	69.176471
					12	-0.980000	66.294118	0.980000	-66.294118
13	0.980000	-63.411765	0.980000	-63.411765
					14	0.980000	-60.529412	0.980000	-60.529412
15	-0.980000	57.647059	-0.980000	57.647059
					16	0.980000	-54.764706	0.980030	-54.764706
17	-0.980000	51.882353	0.986107	-51.882212
					18	0.979947	-49.000002	5.323811	-233.107915
19	2.498959	-96.439338	-4.028866	-46.229822
					20	-0.846112	43.238769	-0.623575	43.108117
21	-5.874143	-40.529506	-4.140310	70.527853
					22	-1.243301	-37.783145	-2.900080	119.154198
23	0.109782	-34.745844	1.107216	-36.352941
					24	-1.886418	15.147839	-0.932524	20.563595
25	-3.973573	147.667518	-16.114888	293.220243
					26	-0.070673	2.070272	3.458688	-142.172420
27	-0.354659	6.042087	4.451287	-87.171526
					28	-0.459049	11.282148	-0.261895	10.823594
29	-0.814922	4.490244	-0.364621	4.838316
					30	-5.519722	-14.393637	26.630394	-111.474262

                            表1神经网络训练结果

附图8显示了神经网络的训练结果，其中，横坐标为显微镜的点扩散参数，单位为象素数，纵坐标为对应的离焦量，单位为微米。从图中可以看到，当微针处于近离焦状态(ε＜0)，且离焦量不是很大时，点扩散参数与离焦量的映射关系有一个突变；据分析，这是由于微针成像本身形成的。通过对其他采样图像组的研究可知，此种情况下映射关系的突变在各组采样图像中都存在，而且是一致的。

4)深度信息提取及结果检验

当点扩散参数与离焦量的函数关系建立起来后，对于一幅离焦图像，计算其扩散参数，并带入上述关系，就可以获得该图像中微针的离焦量，即深度信息。

为了验证深度提取的正确性，将另外一组采样图像作为检验样本，将其点扩散参数带入上述关系中，计算得到离焦量ε′，将ε′与由标尺观察到的离焦量ε进行对比，结果如附图9所示，左图显示了原始采样结果和深度提取结果的对比，图中，“.”描述了原始采样的结果，“+”描述了深度提取的结果；右图显示了ε和ε′的差异。

在上述验证中，共使用离焦量在±70μm内的检验样本67个，结果表明，离焦量的平均误差为1.012μm，80％以上的深度偏差在误差范围(±2μm)内。进一步，可以得到深度提取均方差为3.2279，标准差为1.797。

以同样的方法处理其他采样图像组，得到的结果如下表所示。由此可知，上述深度信息提取方法是有效的。

2μm范围内百分比	样本	82.5％	79.4％	80.6％	90.6％
						标准差	1.797	2.002	1.994	1.5689

                         表2深度提取验证结果

2.微操作机器人系统双针互插实验

为了进一步直观验证深度提取的结果，本发明还设计了微操作机器人系统双针互插实验。在实验中，共使用两根微针，分处于屏幕两侧，与x轴平行，如附图10(a)所示。左侧微针针尖铰粗，针内径约为20μm；右侧微针针尖较细，针尖宽度2μm。实验目标即将右侧微针插入左侧微针中，步骤如下：

1)调整两针针尖位置，使之在屏幕中保持清晰，保存右针清晰图像，作为成像时的被观察物体。

2)随机在z方向上移动右针，使之离焦变模糊(此时，左针依然保持清晰)。

3)对右针进行深度信息提取，获得其离焦量。显示在对话框上的信息为：计算得到的点扩散参数是23.65；深度信息提取结果为46.8μm。

4)将此离焦量作用于左针，使之在z方向上移动该距离，此时，左右两针都是模糊的，且在理论上处于同一水平面。

5)在离焦状态下，完成双针互插操作。

6)调节显微镜，观察结果。

根据上述步骤完成实验，如附图10所示。上述实验可在离焦±70μm范围内顺利完成。

Claims (11)

1.一种基于显微图像处理的微操作工具深度信息提取装置，其特征在于：提取装置包括：倒置显微镜光源(1)、在其光线方向上依次设置的聚光镜(2)和用于将光线分成相互垂直的透射光和反射光两束光的分光棱镜(3)，在分光棱镜(3)的透射光一侧依次设置有位于XOY平面上的微针(4)、用于将微针(4)成像的倒置显微镜物镜(5)和用于摄录图像的电荷耦合器件CCD摄像机(6)，以及通过图像采集接口接收电荷耦合器件CCD摄像机(6)摄录的图像的计算机(11)；在分光棱镜(3)的反射光一侧依次设置有位于YOZ平面的纵向标尺(7)和用于将标尺(7)局部成像的显微镜物镜(8)以及用于摄录图像并将图像传入计算机的电荷耦合器件CCD摄像机(9)；微操作机械手(10)由计算机(11)控制。

2.一种基于显微图像处理的微操作工具深度信息提取方法，其特征在于要经过下述步骤：

1)获取采样图像及标尺图像；

2)读取并处理清晰图像，获得扫描线位置，保存需要数据；

3)按照顺序，依次读取并处理采样图像，保存需要数据；

4)求取采样图像扫描线上的微针边界；

5)计算采样图像的点扩散参数；

6)判断所有采样图像是否已处理完，是则继续，否则转3)；

7)读取所有采样图像的标尺信息，获得离焦量；

8)通过神经网络建立点扩散参数与离焦量的映射关系，保留结果；

9)获取需要进行深度信息提取的模糊图像；

10)读取并处理该模糊图像，保存需要数据；

11)求取模糊图像扫描线上的微针边界；

12)计算模糊图像的点扩散参数；

13)利用神经网络结果，求取模糊图像中物体的离焦量。

3.根据权利要求2所述的提取方法，其特征在于：步骤5)包括下列步骤：

1)随机设置图像的点扩散参数；

2)根据成像模型，计算虚拟图像数据；

3)计算评估值；

4)根据评估值变化，重新获取点扩散参数；

5)判断参数精度是否达到要求，是则继续，否则转2)；

6)保留使评估值较小的扩散参数；

7)随机设置y方向上偏移量；

8)加上偏移量，重新计算虚拟图像数据；

9)计算评估值；

10)根据评估值变化，重新获取偏移量；

11)判断精度是否达到要求，是则继续，否则转8)；

12)保留使评估值较小的偏移量；

13)判断偏移量是否改变，是则继续，否则保留结果，搜索结束；

14)和偏移前比较，判断评估值是否减小，是则继续，否则保留结果搜索结束；

15)保留偏移量，重新搜索点扩散参数，转1)。

4.根据权利要求3所述的提取方法，其特征在于：显微镜成像模型描述如下：由于光的衍射效应和透镜成像的非理想化，成像系统的点扩散函数可用二维高斯函数来近似：

$h(x,y)=\frac{1}{2{\mathrm{\π\σ}}^2}{e}^{-\frac{x^2+y^2}{{2\mathrm{\σ}}^2}}---\left(2\right)$

式中h(x，y)意义为点物体的像，称为点扩散函数，σ是点扩散参数；

卷积能很好地描述显微镜系统的成像，若给定的被观察物体f(u，v)，经一光学系统所成的像为g(x，y)，则有

$g(x,y)=f*h={\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}f(u,v)h(x-u,y-v)\mathrm{dudv}---\left(3\right)$

式中*为卷积操作，h(x，y)为系统的点扩散函数，它决定被观察物体离焦后的扩散程度，也就是所成显微图像g(x，y)的“模糊”程度；h(x，y)由式(2)给出，变量x，y是像空间位置信息，点扩散参数σ唯一确定被观察物体离焦后扩散程度，即显微图像的“模糊”程度；点扩散参数σ通常利用式(2)进行求解，由于得到一个理想的点物体很困难，考虑到显微图像的特殊性，选择利用式(3)的面物体成像模型来辨识σ，将式(2)代入式(3)，得到：

$g(x,y)=\frac{1}{2{\mathrm{\π\σ}}^2}{\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}f(u,v)e^{-\frac{{(x-u)}^2+{(y-v)}^2}{{2\mathrm{\σ}}^2}}\mathrm{dudv}---\left(4\right)$

给定被观察物体f(x，y)和点扩散参数σ，利用式(4)可以求出被观察物体所成的像，像与点扩散参数σ一一对应，为示区别，将计算所得的像称为虚拟像，将实际采得图像称为物体的像；因此，当被观察物体f(x，y)的像已知时，通过与虚拟像对比，可以确定它对应的点扩散参数；基于此，给定一幅待确定σ的被观察物体的像，可以按下面步骤确定显微镜点扩散参数：

1)对于每个可能的σ，利用式(4)生成被观察物体的虚拟像；

2)将虚拟像和物体的像进行对比，求得对应点差值的平方和：

${\underset{x,y}{\mathrm{\Σ}}(g(x,y)-\frac{1}{2{\mathrm{\π\σ}}^2}{\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}f(u,v)e^{-\frac{{(x-u)}^2+{(y-v)}^2}{{2\mathrm{\σ}}^2}}\mathrm{dudv})}^2---\left(5\right)$

3)搜索σ，寻找令虚拟像和物体的像匹配最佳的

$J=\mathrm{Min}\left\{{\underset{x,y}{\mathrm{\Σ}}(g(x,y)-\frac{1}{2{\mathrm{\π}\widehat{\mathrm{\σ}}}^2}{\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}f(u,v)e^{-\frac{{(x-u)}^2+{(y-v)}^2}{{2\widehat{\mathrm{\σ}}}^2}}\mathrm{dudv})}^2\right\}---\left(6\right)$

4)据最小二乘法原理，上一步找到的虚拟像就是物体像的“等价结果”，此时的

就是生成物体像时成像系统的点扩散参数。

5.根据权利要求4所述的提取方法，其特征在于：对于微操作工具是细长型微针的情况，点扩散参数σ通过以下步骤快速确定：

1)对于每个可能的σ，利用式(7)生成被观察物体的虚拟像；

$g(x,y)=\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\·{\∫}_{y_{\min }}^{y_{\max }}f(x,v)e^{-\frac{{(y-v)}^2}{2{\mathrm{\σ}}^2}}\mathrm{dv}---\left(7\right)$

2)将虚拟像和物体的像进行对比，求得对应点差值的平方和：

$\underset{y}{\mathrm{\Σ}}(g(x,y)=\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\·{\∫}_{y_{\min }}^{y_{\max }}f(x,v)e^{-\frac{{(y-v)}^2}{2{\mathrm{\σ}}^2}}\mathrm{dv})---\left(8\right)$

3)搜索σ，寻找令虚拟像和物体的像匹配最佳的

$J=\mathrm{Min}\{\underset{y}{\mathrm{\Σ}}(g(x,y)-\frac{1}{\sqrt{2\mathrm{\π}}\widehat{\mathrm{\σ}}}\·{\∫}_{y_{\min }}^{y_{\max }}f(x,v)e^{-\frac{{(y-v)}^2}{2{\widehat{\mathrm{\σ}}}^2}}\mathrm{dv})\}---\left(9\right)$

4)根据最小二乘法原理，上一步找到的虚拟像就是物体像的“等价结果”，此时的

就是生成物体像时成像系统的点扩散参数σ。

6.根据权利要求5所述的提取方法，其特征在于：扫描线要选取在微针尖端，并且要求每一幅离焦图像中选出来的垂直扫描线必须与从被观察物体的清晰像中选出的那条扫描线直接“对应”，也就是要求在成像过程中，针尖的位置始终不变。

7.根据权利要求5所述的提取方法，其特征在于：点扩散参数σ在搜索的过程中，计算虚拟像和采样样本之间逐点灰度差的平方和时，必须忽略背景区域的大部分点，而只考虑微针所成的像，微针成像的关键信息都包含在微针像中。

8.根据权利要求7所述的提取方法，其特征在于：通过基于中值滤波的方法从扫描线中提取微针所在的区域。

9.根据权利要求5所述的提取方法，其特征在于：搜索点扩散参数σ时采用黄金分割搜索算法。

10.根据权利要求4所述的提取方法，其特征在于：在计算虚拟图像数据时，可以使用离散卷积代替连续卷积，在计算过程中，采用“傅立叶变换-乘积-逆傅立叶变换”来实现离散卷积。

11.根据权利要求2所述的提取方法，其特征在于：点扩散参数σ与离焦量的映射关系借助人工神经网络来建立，采用的是带阈值的、隐含层采用S型传递函数、输出层采用线性传递函数的两层BP神经网络。

Images