CN102175693B - 医用药剂中可见异物的机器视觉检测方法 - Google Patents

Abstract

本发明公开了一种医用药剂中可见异物的机器视觉检测方法，该方法通过以下步骤：(1)获取连续帧图像。(2)图像预处理消除噪声。(3)运动目标提取和分割：基于改进的三帧差分方法运动目标提取。(4)运动目标跟踪：基于SIFT匹配和最邻近匹配方法的跟踪。(5)杂质信息分析：获取杂质的大小和数目以及粗略的分析出杂质的总类。最后通过与国家标准进行判断出该瓶的是否合格。本发明的方法对于医用药剂中可见异物检测十分有效，大大提高了药剂液体中的异物的识别率，降低了药剂的误检率。

Description

医用药剂中可见异物的机器视觉检测方法

技术领域

本发明属于计算机视觉领域，涉及一种医用药剂中可见异物的机器视觉检测方法。

背景技术

透明药剂是我国医药行业重要产品之一，是临床上应用十分广泛的产品。透明药剂是医药产品的一个类别，由于具有迅速、高效、控释的特点，是医疗机构日常必须使用的药品，同时也是临床抢救药物和静脉治疗药物不可缺少的载体或溶酶，在现代临床上占据十分重要的地位。但是这些透明药剂产品的质量是否可靠。影响透明药剂产品质量的因素很多。在输液生产过程中，除了引进先进生产技术及严格的管理手段提高药液本身质量之外，其包装过程也是决定输液产品质量高低的重要因素。目前透明药剂包装的形式有：玻璃瓶、塑料瓶和塑料袋。在我国，约占全国市场90％以上的透明药剂产品还是采用的玻璃瓶包装。由于生产工艺以及封装技术的原因，透明药剂产品中可能含有玻璃碎屑、橡皮屑、纤维等异物。这些异物对人体所造成的危害非常大，如果混同药液注入人体血液后，轻则可能导致循环障碍，组织缺氧，供血不足，从而引起静脉炎、血栓、中风、水肿等疾病，重则导致血管堵塞、破裂，内脏受损，甚至危及生命。为了克服这些异物所带来的危害，必须对透明药剂产品进行细致的检测，行业中称为灯检，因为，这种检测通常是在暗室中用灯光检测的。然而能够自动检测药液中可见异物的机器，叫做灯检机。全世界上的灯检机分为两种，一种是半自动灯检机，一种是全自动灯检机。其中全自动灯检机只有德国、日本、意大利等几个少数各家可以生产。可以预见现在和未来全自动灯检机的市场会全被这些国家占据。为了与这些国家进行竞争，生产出具有自己国家知识产权的全自动灯检机是一件十分必须的事情。然而在我们国家生产时，很多问题在困扰着我们，很多技术难题需要我们攻破，由于我们采用的方法是模拟人眼的处理流程，首先使用机械装置让药瓶运动起来，以便沉浸在瓶底的异物可以飘上来，有助于我们的检测，但是这种装置在处理上很有可能会有气泡和噪声，然而在药液中的异物又十分小，所以处理起来会有以下难点：

1、由于目标图像无形状、尺寸、和纹理信息，可供处理算法利用的信息量小；

2、低信噪比，目标极易被噪声所淹没，单帧处理不能保证对目标的可靠检测。

3、实时性要求高，对算法的运算速度有较高要求。

发明内容

本发明所要解决的技术问题是提出一种医用药剂中可见异物的机器视觉检测方法，该医用药剂中可见异物的机器视觉检测方法检测精度高、准确度高。

本发明的技术解决方案如下：

一种医用药剂中可见异物的机器视觉检测方法，其特征在于，包括以下步骤：

步骤1：对待检测的药液获取连续帧图像；

步骤2：对获取的连续帧图像进行预处理以消除噪声；

步骤3：运动目标提取和分割：使用基于改进的三帧差分方法对运动目标进行提取；

步骤4：运动目标跟踪：使用基于SIFT匹配和最邻近匹配方法的对运动目标进行跟踪，为杂质信息分析提供依据；

步骤5：杂质信息分析：通过运动轨迹分析来区分杂质和气泡。

2.根据权利要求1所述的医用药剂中可见异物的机器视觉检测方法，其特征在于，所述步骤2为：采用中值滤波的方法对获得连续帧图像进行预处理操作：使用5×5的窗口作为滑动窗口，将窗口中各点的灰度值的中值来替代指定点的灰度值。

3.根据权利要求1所述的医用药剂中可见异物的机器视觉检测方法，其特征在于，所述步骤3包括运动目标提取、运动目标分割和运动物体区分三个部分：

A)运动目标提取：采用帧间差分法，通过对连续的三帧图像进行差分处理，第一帧图像与第二帧图像进行差分，第二帧图像与第三帧图像进行差分，再将第一次差分得出的结果与第二次差分得出的结果再进行一次差分；

d(x，y，t)＝|f(x，y，t)-f(x，y，t-1)|；

d(x，y，t+1)＝|f(x，y，t+1)-f(x，y，t)|；

D(x，y，t+1)＝|d(x，y，t+1)-d(x，y，t)|

其中f(x，y，t)为经过中值滤波后的第t帧图像，含有的是第t帧的运动物体；x，y分别为图像的横坐标和纵坐标；

d(x，y，t)就是第t帧图像和第t-1帧图像的差分图像，含有第t帧图像和第t-1帧图像的运动物体；

最后，D(x，y，t+1)是d(x，y，t+1)和d(x，y，t)的差分图像；而D(x，y，t+1)含有的是第t+1帧和第t-1帧的运动物体；

B)运动目标分割：

将获得的D(x，y，t+1)图像进行阈值分割，用上一时刻的前景图的零值区域作为当前时刻的背景区域，以阈值T将前景和背景分开；

阈值T＝α×δ，其中α取3，

N，M分别是图像序列的长度和宽度，AREA是背景区域的面积，

D(x，y，t+1)为t+1时刻的灰度图像；采用阈值T对D(x，y，t+1)二值化，得到t+1时刻的二值图像FD(x，y，t+1)；

【这里的思想“用上一时刻的前景图的零值区域作为当前时刻的背景区域”，所谓前景图就是二值图像，上一时刻就是时刻t，那么上一时刻的前景图就是FD(x，y，t)，上一时刻的前景图的零值区域即FD(x，y，t)＝0中(x，y)的坐标，FD(x，y，t)是t时刻的二值图像，找到这些坐标后，通过计算在D(x，y，t+1)图中这些坐标的值，进而得出一个阈值，再得出t+1时刻的二值图像FD(x，y，t+1)】

C)运动物体区分：

设定FD(x，y，t)为t时刻D(x，y，t)的二值图像，FD(x，y，t+1)和FD(x，y，t+3)均进行边长为2×2的矩形形状的开运算来消除部分噪声，然后对实施了开运算的FD(x，y，t+1)和FD(x，y，t+3)进行取交集运算：

P(x，y，t+1)＝FD(x，y，t+3)I FD(x，y，t+1)；

图像P(x，y，t+1)就是含有第t+1帧运动物体的二值图像，而且图像P(x，y，t+1)中的白色区域表示的就是的图像f(x，y，t+1)中的运动物体。

4.根据权利要求1所述的医用药剂中可见异物的机器视觉检测方法，其特征在于，所述步骤4为：

A)首先计算f(x，y，t)和f(x，y，t+1)两幅图像的运动区域SIFT特征点：

第一步：尺度空间的生成

一幅二维图像的尺度空间定义为：

L(x，y，σ)＝G(x，y，σ)*f(x，y，t)

其中G(x，y，σ)是尺度可变高斯函数，

采用加速了的SIFT算法，对图像的运动区域进行多次与G(x，y，σ)卷积操作，得到一个图像金字塔，

其中G(x，y，σ)是尺度可变高斯函数，

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

构建图像金字塔，图像金字塔共O组，每组有S层，O和S是分别设为3和5，下一组的图像由上一组图像降采样得到，而下一层图像是由上一层图像与G(x，y，σ)进行卷积获得的，用Pic(i，j)表示这个金字塔中第i组第j层图像。那么

Pic(1，1)＝f(x，y，t)

Pic(i，j)＝Pic(i，j-1)*G(x，y，σ)(1≤i≤O，2≤j＜S)，仅对运动目标区域【运动目标区域就是指前面通过分割将源输入灰度图像处理成了二值图像，这里的运动区域就是二值图像中的为白色的那些位置，就是运动物体的区域简称运动区域】，进行卷积操作，并且，Pic(i+1，j)是由Pic(i，j)降采样获得；【σ取0.5-1之间的一个定值，本发明的具体实例中取0.5】

第二步：空间极值点检测

为了寻找尺度空间的极值点，每一个采样点要和它所有的相邻点比较，看其是否比它的图像域和尺度域的相邻点大或者小，即中间的检测点和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点比较，以确保在尺度空间和二维图像空间都检测到极值点，取27个点中的最大值点或最小值点作为空间极值点；

第三步：对空间极值点进行筛选：

采用两种筛选方法中的至少一种：

第一种筛选方法：

对每一个极值点，如果满足

则将该极值点删除；

取γ＝10；

其中，Tr(H)＝F_xx+F_yy，

$\mathrm{Det}\left(H\right)=F_{\mathrm{xx}}{F}_{\mathrm{yy}}-{\left(F_{\mathrm{xy}}\right)}^2,H=\left[\begin{array}{cc}{F}_{\mathrm{xx}}& F_{\mathrm{xy}}\\ F_{\mathrm{xy}}& F_{\mathrm{yy}}\end{array}\right];$

F(x，y，σ)＝(G(x，y，kσ)-G(x，y，σ))*D(x，y，t)＝L(x，y，kσ)-L(x，y，σ)，

L(x，y，σ)＝G(x，y，σ)*D(x，y，t)；

对于F(x，y，σ)，对x求二阶导数、对y求二阶导数以及先对x求偏导再对y求偏导，并将当前极值点的像素坐标(x，y)代入，分别得到F_xx、F_yy和F_xy；其中

第二种筛选方法：

首先计算相对极值点精确位置的偏移向量【即先对F(x，y，σ)求x的二阶偏导，再将极值点的(x，y，σ)代入该二阶偏导中求得具体的偏导值作为偏移向量】，将该偏移向量加到原极值点的坐标中形成新的

再将该新的中，公式中F表示F(x，y，σ)在极值点处的值；其中向量

如果

小于预设阈值，则将当前的极值点删除；【根据Lowe在文章中的建议，该预设阈值可以设定为0.003，假设灰度已经归一化了，也就是说灰度范围是[0，1]】

第四步：给关键点分配方向：【前面提到的极值点其实是候选的关键点，经过第三步后剩下的极值点就是这里说的关键点】

首先在高斯尺度空间计算关键点的梯度模和方向：

$m(x,y)=\sqrt{{(L(x+1,y)-L(x-1,y))}^2+{(L(\mathrm{xy}+1)-L(x,y-1))}^2}$

$\mathrm{\θ}(x,y)=a\tan 2\left(\frac{L(x,y+1)-L(x,y-1)}{L(x+1,y)-L(x-1,y)}\right)$

其中L是原图与高斯核，m是梯度模，θ是梯度方向，其中L所用的尺度为每个关键点各自所在的尺度；在以关键点为中心的领域窗口内采样，用直方图统计领域像素的梯度方向，梯度直方图范围是0°～360°；如果每10°一个柱，则共36个柱；然后将邻域内的每个采样点按梯度方向θ归入适当的柱，以梯度m作为贡献的权重；最后选择直方图的主峰值作为特征点的主方向，选取量值达到主峰值80％以上的局部峰值作为辅助方向。那么一个特征点可能会被指定具有多个方向，具有增强匹配的鲁棒性；

至此，图像的关键点已检测完毕，每个关键点有三个信息：位置、所处尺度、方向；

第五步：特征点描述子生成

【描述子生成过程：以关键点为中心取8×8的窗口。图5左部分的中央黑点为当前关键点的位置，每个小格代表关键点邻域所在尺度空间的一个像素，箭头方向代表该像素的梯度方向，箭头长度代表梯度模值，图中蓝色的圈代表高斯加权的范围(越靠近关键点的像素梯度方向信息贡献越大)。然后在每4×4的小块上计算8个方向的梯度方向直方图，绘制每个梯度方向的累加值，即可形成一个种子点，如图5右部分所示。此图中一个关键点由2×2共4个种子点组成，每个种子点有8个方向向量信息。这种邻域方向性信息联合的思想增强了算法抗噪声的能力，同时对于含有定位误差的特征匹配也提供了较好的容错性。

实际计算过程中，为了增强匹配的稳健性，对每个关键点使用4×4共16个种子点来描述，这样对于一个关键点就可以产生128个数据，即最终形成128维的SIFT特征向量。此时SIFT特征向量已经去除了尺度变化、旋转等几何变形因素的影响，再继续将特征向量的长度归一化，则可以进一步去除光照变化的影响。

128个数据的由来，本来是一个关键点由2*2个小块，4个种子构成，每个种子是8维向量，共4*8＝32维，如图5所示。现在将这个2*2个小块扩大成4*4个小块，那么就由16个种子构成了。共含有向量16*8＝128维了。这128维就是由16个种子构成的。】

以关键点为中心取8×8的窗口；然后在组成8×8窗口的4个4×4的小块上分别计算8个方向的梯度方向直方图，绘制每个梯度方向的累加值，即总共可形成4个种子点；此时经过SIFT【即尺度不变特征变换】变换的特征向量已经去除了几何变形因素的影响，再继续将特征向量的长度归一化，则进一步去除光照变化的影响；

【这里每个运动区域SIFT特征点的数据结构就是一个128维的向量和表示坐标的x和y】

B)其次对f(x，y，t)和f(x，y，t+1)两幅图像的运动区域SIFT特征点进行SIFT匹配：

待匹配的两幅图像为图像f(x，y，t)和图像f(x，y，t+1)，从图像f(x，y，t)中提取的SIFT特征点集合为：

F_t＝{a₁，a₂，L，a_m}，m为图像f(x，y，t)的特征点个数；

从图像f(x，y，t+1)中提取的SIFT特征点集合为：

F_t+1＝{b₁，b₂，L，b_n}，n为图像f(x，y，t+1)的特征点个数；

逐一对集合F_t中每个元素点计算与集合F_t+1中每个元素点的距离，这个距离是指两个向量的欧式距离，假如是I和J这两个特征点的话，而且I₁表示I的第1维数据，那么I₁₂₈表示I的第128维数据，对于特征点J也有相同的表达方式，用D_I，J表示特征点I到特征点J的距离，则

得到一个距离集合，从而得到最小距离d_min和次小【次小是倒数第二小】距离d_n-min，再判断

是否成立，Thresh是预设定的阈值，如果

成立则匹配成功，否则匹配不成功；阈值Thresh设为0.44；【这样做的目的是为了增加匹配的可靠性，假如不考虑这个Thresh的话，就会有这种可能，d_min和d_n-min相差不大，那么F_t中的改点就可以几乎和这两个点都算匹配成功，这样会大大影响这个的匹配精确度。所以设定一个阈值，它的作用就是要提高匹配的精确度，还有一个就是这个公式有问题，应该改成

才是对的，可以对两个很简单的假设进行分析，假设一d_min＝5，d_n-min＝10，和假设二d_min＝5，d_n-min＝20，很明显假设二要比假设一的匹配性更好，但是假设一的

这个比值比假设二的

大，也就意味着如果按照

来判定的话，假设一更容易满足条件，这样就违背了我的初衷，所以这个公式有问题，应该改成才是对的。而且这个精确度是由Thresh这个值决定的，这个值越小说明能匹配点的个数越少，但是精确度越高】

C)最后对f(x，y，t)和f(x，y，t+1)两幅图像的剩余运动区域进行最邻近匹配：

统计P(x，y，t)和P(x，y，t+1)中未匹配的运动物体得到两个集合：

P_t＝{a₁(x，y，s)，a₂(x，y，s)，L，a_m1(x，y，s)}，m1为图像P(x，y，t)中未匹配运动物体的数目，其中x和y表示该运动物体的质心坐标，s表示该运动物体的大小；【大小其实是面积，一个像素一个单位，其实这个面积的意思就是统计该运动区域包含的像素的数目】

P_t+1＝{b₁(x，y，s)，b₂(x，y，s)，L，b_n1(x，y，s)}，n1为图像P(x，y，t+1)中未匹配运动物体的数目，其中x和y表示该运动物体的质心坐标，s表示该运动物体的大小；

第一步骤：逐一对集合P_t中每个运动物体的质心计算与集合P_t+1中每个运动物体质心的距离，得到一个距离集合，得到最小距离d1_min；

第二步骤：由于根据物体运动的连续性，不可能在相邻两帧的时间内，物体的位置会发生太大的变化，于是设定一个最小相邻距离d，只有当d1_min＜d成立的最小距离，才可以用于匹配；【这个是与药瓶中异物的移动最大速度和获取连续图像的时间间隔有关的，设定d是60。】

第三步骤：通过了距离判断，找到一对匹配运动物体；进一步判断是否满足

$\frac{\max \{a_i\left(s\right),b_j\left(s\right)\}}{\min \{a_i\left(s\right),b_j\left(s\right)\}}<S_{\mathrm{num}},$

假定当前匹配的是P_t中的第i个元素和P_t+1中的第j个元素，a_i(s)表示P_t中的第i个元素的面积，b_j(s)表示P_t+1中的第j个元素的面积，S_num为预设的比值，S_num设为2；

第四步骤：匹配完了一对之后，将这两个元素分别从P_t和P_t+1中删除，然后再返回第一步骤进行循环，直到P_t和P_t+1中有一个里面没有元素了，就表示匹配完成。

所述步骤5的杂质信息分析过程如下：基于步骤4中的匹配结果记录运动物体的运动轨迹，再通过两条原则来判断有关运动物体的性质：

(1)根据杂质密度大于液体密度判定轨迹向下的运动物体是杂质；

(2)根据气泡密度小于液体密度判定轨迹向上的运动物体是气泡。

【最后对气泡的数目和大小进行统计分析，与国家标准进行比较最后判断出当前药液是否合格。】

有益效果：

与现有技术相比，本发明的有益效果如下：

首先，在运动目标提取和分割中本发明采用基于改进的三帧差分方法对运动目标进行提取，可以大大提高对运动目标提取的速度和准确率。

其次，在运动物体跟踪方面采用基于SIFT匹配和最邻近匹配方法的对运动目标进行跟踪。应用SIFT运算时，并不是对所有区域进行SIFT运算，前面已经提到，只是对运动区域进行SIFT运算。通过这样的加速大大提高了跟踪的运算速度。而且通过这两种匹配的结合可以大大提高跟踪的准确性。

附图说明

图1为本发明所涉及方法的总的流程图；

图2为两种匹配结合的算法流程图；

图3为输入算法处理的一组含有运动物体的图片；(a)第13帧图像，(b)第15帧图像；

图4为提取出运动物体之后的一组二值黑白图片；(a)第10帧运动物体图像，(b)第15帧运动物体图像；

图5为由关键点邻域梯度信息生成特征向量示意图。

具体实施方式

以下将结合附图和具体实施例对本发明做进一步详细说明：

实施例1：

本发明的药液中可见异物检测和跟踪的方法，具体包括以下步骤：

(1)对待检测的药液获取连续帧图像，以供后面需要。

(2)对图像进行预处理消除部分噪声。

(3)运动目标提取和分割：使用基于改进的三帧差分方法对运动目标进行提取。

(4)运动目标跟踪：使用基于SIFT匹配和最邻近匹配方法的对运动目标进行跟踪，以便为后面的杂质信息分析提供依据。

(5)杂质信息分析：获取杂质的大小和数目以及粗略的分析出杂质的总类。

进一步，在上述基于药液中可见异物检测和跟踪的方法中，步骤(1)是从用机械装置将药瓶旋转后，药瓶停止但是药液还是在运动，里面的异物也仍然在运动，这是通过摄像机获取药液的连续帧图像，为后续处理提供数据。

步骤(2)是对获得的图像进行预处理操作，使用中值滤波的方法。中值滤波的方法是通过采用一个含有奇数个点的滑动窗口，将窗口中各点的灰度值的中值来替代指定点(一般为窗口的中心点)的灰度值。对于一幅图像来说完成中值滤波可以选择不同的窗口形状和尺寸。经过实验，这里我们发现5×5的窗口对后面的处理对后面的处理最有帮助。

步骤(3)是提取运动目标和分割，由于药液中的可见异物是运动的，我们对药液中可见异物的提取也就是对药液中运动物体提取。一般有两类方法，一种是帧间差分法，一种是背景建模法。由于在这样的环境下，图像与建模出来的背景差别不是很大，所以不能使用背景建模的方法。我们采用帧间差分法，通过对连续的三帧图像进行差分处理，第一帧图像与第二帧图像进行差分，第二帧图像与第三帧图像进行差分，这样完成之后，再将第一次差分得出的结果与第二次差分得出的结果再进行一次差分。

d(x，y，t)＝|f(x，y，t)-f(x，y，t-1)|

d(x，y，t+1)＝|f(x，y，t+1)-f(x，y，t)|

D(x，y，t+1)＝|d(x，y，t+1)-d(x，y，t)|

其中f(x，y，t)为经过5×5窗口的中值滤波后的第t帧图像，当然它分别含有的是第t帧的运动物体。d(x，y，t)就是第t帧和第t-1帧的差分图像。它含有第t帧和第t-1帧的运动物体。最后，D(x，y，t+1)是d(x，y，t+1)和d(x，y，t)的差分图像。而D(x，y，t+1)含有的是第t+1帧和第t-1帧的运动物体。

阈值分割：

将获得的D(x，y，t+1)图像进行阈值分割，将前景和背景分开。阈值的选取可通过计算当前图像背景区域的均方差值δ来实时估计。由于一般相邻时刻背景区域变化不大，可用上一时刻的前景图的零值区域(即背景区域)作为当前时刻的背景区域。

其具体计算公式如下：

D(x，y，t+1)为t+1时刻的灰度图像，我们设定FD(x，y，t+1)为t+1时刻的二值图像，N，M分别是图像序列的长度和宽度，AREA是背景区域的面积，T为设定的阈值。

当FD(x，y，t)＝0时的处理：

$\stackrel{\&OverBar;}{D(t+1)}=\frac{1}{\mathrm{AREA}}\underset{x=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{M}{\mathrm{\&Sigma;}}}D(x,y,t+1)$

$\sqrt{\frac{\underset{x=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{M}{\mathrm{\&Sigma;}}}{[D(x,y,t+1)-\stackrel{\&OverBar;}{D(t+1)}]}^2}{\mathrm{AREA}}}$

T＝α×δ

α值的设定与场景的复杂度有关，一般的场景取3，如果场景复杂的话，可以选取大一点4～8。就我们的场景，我们把α定为3。这样就可以将运动物体提取出来，但是提取出来的不够彻底，是因为它是对含有两帧运动物体图像的D(x，y，t+1)进行处理的，下面还应该做的就是将这两帧的物体区分开来。

运动物体区分：

具体方法就是通过两幅含有两帧运动物体的图像求他们的交集，由上面的运算我们可以知道D(x，y，t+1)含有的是第t+1帧和第t-1帧的运动物体，同理也可以知道D(x，y，t+3)含有的是第t+3帧和第t+1帧的运动物体，这样一来如果对D(x，y，t+1)和D(x，y，t+3)求一次“交集”，就可以将含有第t+1帧，这里当然不是指直接对D(x，y，t+1)和D(x，y，t+3)进行操作，而是对他们的二值图像FD(x，y，t+1)和FD(x，y，t+3)进行操作。

FD(x，y，t+1)首先进行边长为2×2的矩形形状的开运算来消除部分噪声，FD(x，y，t+3)也进行这样操作，最后在对结束了开运算的FD(x，y，t+1)和FD(x，y，t+3)进行取交集运算。

P(x，y，t+1)＝FD(x，y，t+3)I FD(x，y，t+1)

则目前P(x，y，t+1)就是含有第t+1帧运动物体的图像，而且它已经是一张二值图像了。而且图像P(x，y，t+1)的白色区域表示的就是的f(x，y，t+1)运动物体。

接下来，对P(x，y，t+1)进行统计计算，统计出里面运动物体的数目，每个运动物体的大小和质心。并对他们进行编号，以供后面处理。

步骤(4)是运动目标跟踪使用基于SIFT匹配和最邻近匹配方法的对运动目标进行跟踪。

首先介绍一下SIFT匹配算法：

SIFT算法是一种提取局部特征的算法，在尺度空间寻找极值点，提取位置，尺度，旋转不变量。一幅图像通过该算法的运算可以求出很多SIFT特征点。这种特征变换方法使图像的局部特征，对旋转、尺度缩放和亮度变化保持不变，对视角变化、仿射变换和噪声也保持一定的稳定性，这也正是我们选择这种算法的原因，正因为它有着对旋转特性保持不变的特点。

第一步：尺度空间的生成

一幅二维图像的尺度空间定义为：

L(x，y，σ)＝G(x，y，σ)*f(x，y，t)

其中G(x，y，σ)是尺度可变高斯函数，

为了有效的在尺度空间检测到稳定的关键点，采用高斯差分尺度空间(DOG scale-space)。利用不同尺度的高斯差分核与图像卷积生成。

F(x，y，σ)＝(G(x，y，kσ)-G(x，y，σ))*I(x，y)＝L(x，y，kσ)-L(x，y，σ)

构建图像金字塔：图像金字塔共O组，每组有S层(这里的O和S是由程序设定的)，下一组的图像由上一组图像降采样得到了。而下一层图像是由上一层图像与G(x，y，σ)进行卷积获得的。不过这里我们采用加速的SIFT运算，就是只是对图像的运动物体区域进行与G(x，y，σ)的卷积运算。在这里我们用Pic(i，j)表示这个金字塔中第i组第j层图像。那么

Pic(1，1)＝f(x，y，t)

Pic(i，j)＝Pic(i，j-1)*G(x，y，σ)(1≤i≤O，2≤j＜S)

并且，Pic(i+1，j)是由Pic(i，j)降采样获得。

第二步：空间极值点检测

为了寻找尺度空间的极值点，每一个采样点要和它所有的相邻点比较，看其是否比它的图像域和尺度域的相邻点大或者小。如图3所示，中间的检测点和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点比较，以确保在尺度空间和二维图像空间都检测到极值点。

第三步：精确确定极值点位置

通过拟和三维二次函数以精确确定关键点的位置和尺度(达到亚像素精度)，同时去除不稳定的边缘响应点，以增强匹配稳定性、提高抗噪声能力。

边缘响应的去除：

一个定义不好的高斯差分算子的极值在横跨边缘的地方有较大的主曲率，而在垂直边缘的方向有较小的主曲率。主曲率通过一个2x2的Hessian矩阵H求出：

$H=\left[\begin{array}{cc}{F}_{\mathrm{xx}}& F_{\mathrm{xy}}\\ F_{\mathrm{xy}}& F_{\mathrm{yy}}\end{array}\right]$

导数由采样点相邻差估计得到。

F的主曲率和H的特征值成正比，令α为最大特征值，β为最小的特征值，则

Tr(H)＝F_xx+F_yy

Det(H)＝F_xxF_yy-(F_xy)²

的值在两个特征值相等的时候最小，随着r的增大而增大，因此，为了检测主曲率是否在某域值r下，只需检测

$\frac{\mathrm{Tr}{\left(H\right)}^2}{\mathrm{Det}\left(H\right)}<\frac{{(\mathrm{\&gamma;}+1)}^2}{\mathrm{\&gamma;}}$

通过实验参考Lowe的文章，我们取γ＝10。

第四步：给关键点分配方向

利用关键点邻域像素的梯度方向分布特性为每个关键点指定方向参数，使算子具备旋转不变性。首先在高斯尺度空间计算特征点的梯度模和方向公式：

$m(x,y)=\sqrt{{(L(x+1,y)-L(x-1,y))}^2+{(L(x,y+1)-L(x,y)-1)}^2}$

$\mathrm{\&theta;}(x,y)=a\tan 2\left(\frac{L(x,y+1)-L(x,y-1)}{L(x+1,y)-L(x-1,y)}\right)$

其中L是原图与高斯核的卷积，m是梯度模，θ是梯度方向，其中L所用的尺度为每个关键点各自所在的尺度。在实际计算时，我们在以关键点为中心的领域窗口内采样，用直方图统计领域像素的梯度方向。梯度直方图范围是0°～360°。如果每10°一个柱，则共36个柱。然后将邻域内的每个采样点按梯度方向θ归入适当的柱，以梯度m作为贡献的权重。最后选择直方图的主峰值作为特征点的主方向，选取量值达到主峰值80％以上的局部峰值作为辅助方向。那么一个特征点可能会被指定具有多个方向，具有增强匹配的鲁棒性。

至此，图像的关键点已检测完毕，每个关键点有三个信息：位置、所处尺度、方向。

第五步：特征点描述子生成

描述自生成过程：以关键点为中心取8×8的窗口。图a左部分的中央黑点为当前关键点的位置，每个小格代表关键点邻域所在尺度空间的一个像素，箭头方向代表该像素的梯度方向，箭头长度代表梯度模值，图中蓝色的圈代表高斯加权的范围(越靠近关键点的像素梯度方向信息贡献越大)。然后在每4×4的小块上计算8个方向的梯度方向直方图，绘制每个梯度方向的累加值，即可形成一个种子点，如图a右部分所示。此图中一个关键点由2×2共4个种子点组成，每个种子点有8个方向向量信息。这种邻域方向性信息联合的思想增强了算法抗噪声的能力，同时对于含有定位误差的特征匹配也提供了较好的容错性。

实际计算过程中，为了增强匹配的稳健性，对每个关键点使用4×4共16个种子点来描述，这样对于一个关键点就可以产生128个数据，即最终形成128维的SIFT特征向量。此时SIFT特征向量已经去除了尺度变化、旋转等几何变形因素的影响，再继续将特征向量的长度归一化，则可以进一步去除光照变化的影响。

第六步：两幅图像特征点的匹配

经过上述过程，两幅输入图像就可以分别获得SIFT特征点，此时每个特征点都已经确定了它的确切位置，尺度和方向。我们就可以利用局部特征描述子进行特征点的匹配。步骤如下：

待匹配的两幅图像为图像A和图像B，从图像A中提取的SIFT特征点集合为：

F_a＝{a₁，a₂，L，a_m}，m为图像A的特征点个数。

从图像B中提取的SIFT特征点集合为：

F_b＝{b₁，b₂，L，b_n}，n为图像B的特征点个数。

逐一对集合F_a中每个元素点计算与集合F_b中每个元素点的距离，得到一个距离集合，对距离集合进行排序，得到最小距离d_min和次小距离d_n-min。然后再判断一下：

(Thresh是我们设定的一个阈值)。

式子是否成立，如果成立则匹配成功，否则匹配不成功。降低这个阈值，SIFT匹配的点会减少，但是会增加匹配的稳定性。这里这个阈值我们设为0.44。

其次介绍基于SIFT匹配和最邻近匹配的跟踪算法：

我们这里的跟踪实质是指匹配，由于对P(x，y，t)的统计我们已经得出了运动物体的一些基本数据，如数目，大小，质心等等，而且给所有的运动物体都设定了标号。等到下一帧处理时，我们处理完P(x，y，t+1)，就又会获得这些数据。要对运动进行跟踪，就是找到两帧得出物体的匹配情况，有三种基本情况需要考虑，第一种就是第t帧的运动物体在第t+1帧还在(运动物体还在界面内)，第二种就是第t帧的运动物体在第t+1帧已经不在了(运动物体离开了可以观察的界面)，第三种就是第t帧没有的运动物体在第t+1帧出现了(运动物体的突然出现)。对于这三种情况可以一起考虑。

具体算法如下：

首先，为了进行SIFT算法的加速，通过对两帧图像f(x，y，t)和f(x，y，t+1)的运动物体区域(由于步骤(3)已经得出了两帧图像的运动物体数目，以及运动物体的质心和大小)进行SIFT算法计算，而不是对两帧图像完全进行SIFT算法计算，可以找出两幅图像运动区域的匹配点，进而可以对该运动区域进行匹配。

其次，对剩下没有匹配到的物体再使用最邻近匹配算法，当然这算法是有限制的，就是根据物体运动的连续性，不可能在相邻两帧的时间内，物体的位置会发生太大的变化，所以这里有一个最邻近匹配距离，来表示两个物体匹配的最短距离，如果剩下的两个物体的距离超过了这个距离的话，那么也不会对他们进行匹配。

其中最邻近匹配的思路是：

统计P(x，y，t)和P(x，y，t+1)中未匹配的运动物体得到两个集合：

P_t＝{a₁(x，y，s)，a₂(x，y，s)，L，a_m(x，y，s)}，m为图像P(x，y，t)中未匹配运动物体的数目，其中x和y表示该运动物体的质心坐标，s表示该运动物体的大小。

P_t+1＝{b₁(x，y，s)，b₂(x，y，s)，L，b_n(x，y，s)}，n为图像P(x，y，t+1)中未匹配运动物体的数目，其中x和y表示该运动物体的质心坐标，s表示该运动物体的大小。

第一步：逐一对集合P_t中每个运动物体的质心计算与集合P_t+1中每个运动物体质心的距离，得到一个距离集合，对距离集合进行排序，得到最小距离d_min。

第二步：由于根据物体运动的连续性，不可能在相邻两帧的时间内，物体的位置会发生太大的变化，于是设定一个最小相邻距离d，只有当d_min＜d成立的最小距离，才可以用于匹配。

第三步：通过了距离判断，找到一对匹配运动物体。这里还要考虑的一个条件就是他们的大小，既然能匹配上就说明他们是同一个物体，他们的大小应该相差不大，然而物体的大小在图像中就是表示面积，所以还应该满足下面的式子：

$\frac{\max \{a_i\left(s\right),b_j\left(s\right)\}}{\min \{a_i\left(s\right),b_j\left(s\right)\}}<S_{\mathrm{num}}$

假定当前匹配的是P_t中的第i个元素和P_t+1中的第j个元素，a_i(s)表示P_t中的第i个元素的面积，b_j(s)表示P_t+1中的第j个元素的面积，S_num表示的是我们设定的规定比值，我们设为2。

第四步：匹配完了一对之后，将这两个元素分别从P_t和P_t+1中删除，然后在进行第一步运算，如果P_t和P_t+1中有一个里面没有元素了，就表示匹配完成了，进行下一步的轨迹判断。

通过以上两种方法的匹配结合，可以把我们考虑的三种情况中的第一种情况给检测出来，至于第二种和第三种情况的运动物体是无法匹配上的，但是对于第二种情况物体消失我们已经记录了它在观察界面的运动轨迹，可以对它进行判断了，对于第三种情况的物体新出现，我们会在它出现时记录它的相关信息，等到下次进行匹配时，就可以匹配上了。

步骤(5)是杂质信息分析。

通过以上的匹配，我们可以跟踪到运动物体，可以记录运动物体的运动轨迹。

我们通过两条原则来判断有关运动物体的性质。

(1)由于杂质密度一般比液体大，所以我们判定轨迹向下的运动物体是杂质。

(2)由于气泡密度一般比液体小，所以我们判定轨迹向上的运动物体是气泡，不是杂质。最终通过运动轨迹分析来区分杂质和气泡，最后对气泡的数目和大小进行统计分析，与国家标准进行比较最后判断出它是否合格。

Claims (1)

1.一种医用药剂中可见异物的机器视觉检测方法，其特征在于，包括以下步骤： 

步骤1：对待检测的药液获取连续帧图像； 

步骤2：对获取的连续帧图像进行预处理以消除噪声： 

采用中值滤波的方法对获得的连续帧图像进行预处理操作：使用5×5的窗口作为滑动窗口，将窗口中各点的灰度值的中值来替代指定点的灰度值； 

步骤3：运动目标提取和分割：使用基于改进的三帧差分方法对运动目标进行提取； 

步骤4：运动目标跟踪：使用基于SIFT匹配和最邻近匹配方法对运动目标进行跟踪，为杂质信息分析提供依据； 

步骤5：杂质信息分析：通过运动轨迹分析来区分杂质和气泡； 

所述步骤3包括运动目标提取、运动目标分割和运动物体区分三个部分： 

A)运动目标提取：采用帧间差分法，通过对连续的三帧图像进行差分处理，第一帧图像与第二帧图像进行差分，第二帧图像与第三帧图像进行差分，再将第一次差分得出的结果与第二次差分得出的结果再进行一次差分； 

d(x，y，t)＝|f(x，y，t)-f(x，y，t-1)|； 

d(x，y，t+1)＝|f(x，y，t+1)-f(x，y，t)|； 

D(x，y，t+1)＝|d(x，y，t+1)-d(x，y，t)| 

其中f(x，y，t)为经过中值滤波后的第t帧图像，含有的是第t帧的运动物体；x，y分别为图像的横坐标和纵坐标； 

d(x，y，t)就是第t帧图像和第t-1帧图像的差分图像，含有第t帧图像和第t-1帧图像的运动物体； 

最后，D(x，y，t+1)是d(x，y，t+1)和d(x，y，t)的差分图像；而D(x，y，t+1)含有的是第t+1帧和第t-1帧的运动物体； 

B)运动目标分割： 

将获得的D(x，y，t+1)图像进行阈值分割，用上一时刻的前景图的零值区域作为当前时刻的背景区域，以阈值T将前景和背景分开； 

阈值T＝α×δ，其中α取3，

N，M分别是图像序列的长度和宽度，AREA是背景区域的面积，D(x，y，t+1)为t+1时刻的灰度图像；采用阈值T对D(x，y，t+1)二值化，得到t+1时刻的二值图像FD(x，y，t+1)； 

C)运动物体区分： 

设定FD(x，y，t)为t时刻D(x，y，t)的二值图像，FD(x，y，t+1)和FD(x，y，t+3)均进行边长为2×2的矩形形状的开运算来消除部分噪声，然后对实施了开运算的FD(x，y，t+1)和FD(x，y，t+3)进行取交集运算： 

P(x，y，t+1)＝FD(x，y，t+3)∩FD(x，y，t+1)； 

图像P(x，y，t+1)就是含有第t+1帧运动物体的二值图像，而且图像P(x，y，t+1)中的白色区域表示的就是图像f(x，y，t+1)中的运动物体； 

所述步骤4为： 

A)首先计算f(x，y，t)和f(x，y，t+1)两幅图像的运动区域SIFT特征点： 

第一步：尺度空间的生成 

一幅二维图像的尺度空间定义为： 

L(x，y，σ)＝G(x，y，σ)*f(x，y，t) 

其中G(x，y，σ)是尺度可变高斯函数，

采用加速了的SIFT算法，对图像的运动区域进行多次与G(x，y，σ)卷积操作，得到一个图像金字塔， 

其中G(x，y，σ)是尺度可变高斯函数， 

构建图像金字塔，图像金字塔共O组，每组有S层，O和S是分别设为3和5，下一组 的图像由上一组图像降采样得到，而下一层图像是由上一层图像与G(x，y，σ)进行卷积获得的，用Pic(i，j)表示这个金字塔中第i组第j层图像；那么 

Pic(1，1)＝f(x，y，t) 

Pic(i，j)＝Pic(i，j-1)*G(x，y，σ)(1≤i≤O，2≤j＜S)，仅对运动目标区域，进行卷积操作，并且，Pic(i+1，j)是由Pic(i，j)降采样获得； 

第二步：空间极值点检测 

为了寻找尺度空间的极值点，每一个采样点要和它所有的相邻点比较，看其是否比它的图像域和尺度域的相邻点大或者小，即中间的检测点和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点比较，以确保在尺度空间和二维图像空间都检测到极值点，取27个点中的最大值点或最小值点作为空间极值点； 

第三步：对空间极值点进行筛选： 

首先计算相对极值点精确位置的偏移向量

将该偏移向量加到原极值点的坐标中形成新的

再将该新的代入中，公式中F表示F(x，y，σ)在极值点处的值；其中向量

如果

小于预设阈值，则将当前的极值点删除； 

第四步：给关键点分配方向： 

首先在高斯尺度空间计算关键点的梯度模和方向： 

其中L是原图与高斯核，m是梯度模，θ是梯度方向，其中L所用的尺度为每个关键点各自所在的尺度；在以关键点为中心的领域窗口内采样，用直方图统计领域像素的梯度方向，梯度直方图范围是0°～360°；如果每10°一个柱，则共36个柱；然后将邻域内的每个采样点按 梯度方向θ归入适当的柱，以梯度m作为贡献的权重；最后选择直方图的主峰值作为特征点的主方向，选取量值达到主峰值80％以上的局部峰值作为辅助方向 ,那么一个特征点可能会被指定具有多个方向，具有增强匹配的鲁棒性； 

至此，图像的关键点已检测完毕，每个关键点有三个信息：位置、所处尺度、方向； 

第五步：特征点描述子生成 

以关键点为中心取8×8的窗口；然后在组成8×8窗口的4个4×4的小块上分别计算8个方向的梯度方向直方图，绘制每个梯度方向的累加值，即总共可形成4个种子点；此时经过SIFT变换的特征向量已经去除了几何变形因素的影响，再继续将特征向量的长度归一化，则进一步去除光照变化的影响； 

B)其次对f(x，y，t)和f(x，y，t+1)两幅图像的运动区域SIFT特征点进行SIFT匹配： 

待匹配的两幅图像为图像f(x，y，t)和图像f(x，y，t+1)，从图像f(x，y，t)中提取的SIFT特征点集合为： 

F_t＝{a₁，a₂，…，a_m}，m为图像f(x，y，t)的特征点个数； 

从图像f(x，y，t+1)中提取的SIFT特征点集合为： 

F_t+1＝{b₁，b₂，…，b_n}，n为图像f(x，y，t+1)的特征点个数； 

逐一对集合F_t中每个元素点计算与集合F_t+1中每个元素点的距离，这个距离是指两个向量的欧式距离，假如是I和J这两个特征点的话，而且I₁表示I的第1维数据，那么I₁₂₈表示I的第128维数据，对于特征点J也有相同的表达方式，用D_I，J表示特征点I到特征点J的距离，则

得到一个距离集合，从而得到最小距离d_min和次小距离d_n-min，再判断

是否成立，Thresh是预设定的阈值，如果成立则匹配成功，否则匹配不成功；阈值Thresh设为0.44； 

C)最后对f(x，y，t)和f(x，y，t+1)两幅图像的剩余运动区域进行最邻近匹配： 

统计P(x，y，t)和P(x，y，t+1)中未匹配的运动物体得到两个集合： 

P_t＝{a₁(x，y，s)，a₂(x，y，s)，…，a_m1(x，y，s)}，m1为图像P(x，y，t)中未匹配运动物体的数目，其中x和y表示该运动物体的质心坐标，s表示该运动物体的大小； 

P_t+1＝{b₁(x，y，s)，b₂(x，y，s)，…，b_n1(x，y，s)}，n1为图像P(x，y，t+1)中未匹配运动物体的数目，其中x和y表示该运动物体的质心坐标，s表示该运动物体的大小； 

第一步骤：逐一对集合P_t中每个运动物体的质心计算与集合P_t+1中每个运动物体质心的距离，得到一个距离集合，得到最小距离d1_min； 

第二步骤：由于根据物体运动的连续性，不可能在相邻两帧的时间内，物体的位置会发生太大的变化，于是设定一个最小相邻距离d，只有当d1_min＜d成立的最小距离，才可以用于匹配； 

第三步骤：通过距离判断，找到一对匹配运动物体；进一步判断是否满足 

假定当前匹配的是P_t中的第i个元素和P_t+1中的第j个元素，a_i(s)表示P_t中的第i个元素的面积，b_j(s)表示P_t+1中的第j个元素的面积，S_num为预设的比值，S_num设为2； 

第四步骤：匹配完了一对之后，将这两个元素分别从P_t和P_t+1中删除，然后再返回第一步骤进行循环，直到P_t和P_t+1中有一个里面没有元素了，就表示匹配完成； 

所述步骤5的杂质信息分析过程如下：基于步骤4中的匹配结果记录运动物体的运动轨迹，再通过两条原则来判断有关运动物体的性质： 

(1)根据杂质密度大于液体密度判定轨迹向下的运动物体是杂质； 

(2)根据气泡密度小于液体密度判定轨迹向上的运动物体是气泡。 

Images