一种基于图像处理的纤维长度自动测量方法

龚 平，梁 菁

(1.梧州学院 电子与信息工程学院，广西 梧州 543002；2.梧州学院 大数据与软件工程学院，广西 梧州 543002)

长度是纤维的重要形态参数之一，也是评价纤维品质与特性的重要指标之一，直接关系其加工性能和使用价值。随着纤维测量技术的发展，结合显微镜成像并应用计算机图像处理的方法进行纤维长度的测量已经越来越受到关注。但是，由于纤维形态不一、细小、柔软且容易絮聚，给基于图像处理的纤维长度自动测量带来了巨大的困难和挑战。点击测量依靠人工用鼠标分段点击纤维的不同部位，计算机通过计算各点间线段的像素数量自动计算纤维长度和宽度，然后再进行统计，每像素代表的实际尺寸用测微尺校准，该方法测量准确，但测量速度较慢[1]。文献[2]采取人工将棉纤维单根化并进行图像采集，保证获取的图像中纤维无交叉，再通过对图像二值化、去除噪声、补断、细化和剪枝等处理，统计每根纤维的像素个数，实现测量棉纤维长度，该方法可提高棉纤维长度分布测量的准确性与可靠性，但测量速度较慢。

为提高纤维长度测量的效率，针对在显微镜下获取的无交叉的植物纤维图像，本研究提出一种自动测量长度的方法。首先对采集到的纤维图像进行背景和杂质分割，得到待测量的纤维主体；然后对纤维进行骨架化处理；最后统计像素个数并通过校准的像素参数计算纤维长度。

1 方法

1.1 提取纤维轮廓

针对采集到的无交叉的植物纤维图像，本研究自动测量纤维长度的方法流程见图1。

图1 自动测量纤维长度方法的流程

先将采集到的无交叉植物纤维图像转换为灰度图，再进行Canny边缘检测[3-4]，然后提取检测得到的各轮廓。步骤如下。

Step 1：使用高斯滤波器平滑图像，降低噪声；

Step 2：使用一阶的梯度算子计算图像中每个像素点的梯度；

Step 3：对梯度值应用非极大值抑制(Non-Maximum Suppression)，使局部灰度变化大的像素保留下来；

Step 4：应用双阈值(Double-Threshold)筛选和连接边缘。设置高阈值和低阈值的比例为3∶1。梯度值高于高阈值的像素点定义为强边缘像素，保留该像素；删除梯度值低于低阈值的像素点。对于周围无强边缘像素的像素点，从梯度值介于双阈值之间的像素点中进行选择，使得边缘尽可能的闭合。

得到图像边缘后，借助OpenCV中Imgproc.findContours()方法得到其轮廓。

1.2 剔除杂质

由于玻片间有气泡、载玻片上有污点、纤维有细小部分等原因，采集到的纤维图像会存在杂质，杂质也会被检测到轮廓(如图2所示)，因此需要对上述过程得到的各轮廓进行筛选。

a.含有气泡和污点的纤维图像;b.含有毛细纤维和污点的图像;c.对a提取的轮廓;d.对a提取的轮廓

针对每一个轮廓，本方法通过OpenCV的Imgproc.minAreaRect()方法得到包含该轮廓最小面积的斜矩形，再计算每一个矩形的面积。大量实验表明，相对于有效的待检测纤维而言，“杂质”具有较小的面积，可设置阈值移除小面积矩形对应的轮廓。此外，有效的待检测纤维应具有一定的长度，即其对应的斜矩阵应具有较大的长宽比例。基于对大量类似于图2-c，d的轮廓对应矩形面积值的数据统计，本研究方法设置阈值为Smax×p，其中Smax表示求取的最大的矩形面积值，p取值为0.05～0.15时准确率较高；设置剔除长宽比大于4∶3的矩形对应的轮廓。

对图2-c，d中轮廓剔除“杂质”后得到如图3所示的待测纤维的轮廓。

a.剔除气泡保留的纤维轮廓;b.剔除细小纤维保留的纤维轮廓

1.3 提取纤维骨架

对纤维主体进行骨架化处理前，先对其进行基本的腐蚀和膨胀处理，目的是使轮廓中的断点得以连续、使轮廓内的部分得以填充。对图3中轮廓进行腐蚀、膨胀后的效果如图4所示。

a.图3-a的腐蚀膨胀结果;b.图3-b的腐蚀膨胀结果

纤维长度自动测量的关键是细化纤维得到单像素宽度的骨架。针对不同形态的纤维，一些经典的细化算法得到的纤维骨骼会有不同程度的毛刺。本研究采用牟少敏[5]等人提出的改进的细化算法，该算法在Zhang[6]细化算法的基础上，增加了3类删除像素的条件，能够较好的保持骨骼和剔除毛刺。

设当前正在处理的像素点为p，对其八邻域编号为p1～p8，如图5所示。

图5 像素点p及八邻域像素编号

细化算法具体描述为：

算法 Thinning(I,i)

输入：待细化的二值图像I

输出：单像素宽度的骨架图像Y，遍历参数i

Step1.遍历访问图像中I所有的像素点p

Step2.如果Ip为0则Yp=0

否则，进行Step3-Step9

Step3.计算b=p1+p2+p3+p4+p5+p6+p7+p8

Step4.考察八邻域中相邻像素的关系：

t1=(p2 == 0 && p3 == 1)?1:0 t2=(p3 == 0 && p4 == 1)?1:0

t3=(p4 == 0 && p5 == 1)?1:0 t4=(p5 == 0 && p6 == 1)?1:0

t5=(p6 == 0 && p7 == 1)?1:0 t6=(p7 == 0 && p8 == 1)?1:0

t7=(p8 == 0 && p1 == 1)?1:0 t8=(p1 == 0 && p2 == 1)?1:0

Step5.计算a=t1+t2+t3+t4+t5+t6+t7+t8

Step6.计算 m1 = i==1?p2×p4×p6:p2×p4×p8

m2 = i==1?p4×p6×p8:p2×p6×p8

Step7.计算mp=p2+p3×2+p4×4+p5×8+p6×16+p7×32+p8×64+p1×128

Step8.进行布尔计算 flag=(mp==65||mp==5||mp==20||mp==80||mp==13||mp==22||mp==52||mp==133||mp==141||mp==54)

Step9.如果(b >= 2 && b <= 6)&&(a == 1||flag)&&(m1 == 0 && m2 == 0)为true，则Yp=0

将腐蚀、膨胀后的结果进行二值化，使用上述算法进行纤维细化得到如图6所示的纤维骨架。

a.含有气泡各污点的纤维图像;b.含有细小纤维和污点的图像

1.4 追踪像素链

细化后得到了纤维的单像素宽度的骨架，有p、p1～p8位置上的像素值取值非0即1。

将骨架图像中的像素分为三类进行处理：

(1)非骨架区域的点，有p+p1+p2+…+p8=0；

(2)骨架的端点(起点或终点)，有p+p1+p2+…+p8=2；

(3)骨架中的衔接点，有p+p1+p2+…+p8≥3(对于不存在毛刺或分叉的骨架，有p+p1+p2+…+p8=3)。

(4)在骨架图像中首先寻找骨架的端点，记为s，然后将其八邻域中访问值为1的像素记为q，同时将已访问过的s记为0；再以q作为新的起点，记为s，重复上述操作，直至访问到骨架的另一个端点。此过程中记录像素链中当前起点像素和后继像素的相邻关系，即记录出现图5所示p1-p、p3-p、p5-p或p7-p相邻的次数k1，记录像素链中出现p2-p、p6-p相邻的次数k2，记录像素链中出现p4-p、p8-p相邻的次数k3。若该条纤维骨架的像素个数n，有n=k1+k2+k3-1。

a.含有气泡各污点的纤维图像;b.含有细小纤维和污点的图像

1.5 换算纤维长度

不同的图像采集系统观测到的纤维图像尺寸会有不同。在计算纤维长度时需要先行标定当前环境下的单位长度。在特定的物镜倍数下，采集标尺图像(如图8-a,b)和待测纤维图像，使二者的图像尺寸一致。为减少误差，本方法并不直接标定单个像素的长度，而是记录1 mm对应在水平方向、垂直方向各占多少像素，再按比例计算得到纤维的长度。

a.水平方向1 mm标尺成像;b.垂直方向1 mm标尺成像

假设图像宽度为w，高度为h，1 mm在水平方向长度为tw个像素，在垂直方向为th个像素。对于相同标准下采集的尺寸为w*h的待测纤维图像，其中的纤维长度len通过式(1)进行计算。

(1)

2 测试与试验

测试不同形态的无交叉纤维图像，验证本方法的鲁棒性。试验结果表明，本方法对背景不复杂的无交叉纤维能以较高的准确率得到纤维的骨架，从而计算得到纤维长度(如图9-a,b所示)；但存在有一定程度的漏检与误检(如图9-c所示)，结合人工删除误检的纤维骨架、点击增加漏检的纤维路径，效果会更好。

a.准确提取无交叉的机械阔叶木桨纤维;b.准确提取无交叉的棉纤维;c.漏检右上角有效纤维

经典的人工点击纤维路径结合计算机逐段计算测量的人机交互方法具有较高的准确性和可靠性，但速度较慢，对人工判断有较大的依赖。本研究自动测量方法计算与人机交互逐段测量方法的部分计算结果对比见表1。

表1 自动测量与人机交互逐段测量结果的对比 mm

不考虑经典人机交互测量的误差，从表1中可以看出，本研究自动测量的方法在部分结果上与人机交互测量结果平均有0.015 06的差异，此差异主要来自于对纤维骨架像素的选择。从时间效率方面来看，该方法在速度上明显优于人机交互逐段测量的方法。

4 总结

本研究提出并实现了一种基于图像处理技术的无交叉纤维自动测量长度的方法。测试结果表明，该方法能够针对不同形态的无交叉纤维图像实现长度的自动测量。后续改进的思路：一方面尝试自适应的阈值算法，以提高纤维主体提取的准确性;另一方面改进纤维骨架的细化算法，以便更有效地处理单像素骨架中的毛刺。