0引言为了保证缫丝的质量，在蚕茧缫丝前需要进行选茧工作，以剔除影响缫丝质量的次、下茧。根据国家标准(中华人民共和国国家质量监督检验检疫总局和中国国家标准化管理委员会，2015)中的规定，下茧主要包括双宫、薄皮、黄斑茧、口茧、柴印茧、畸形茧等几大类。其中黄斑茧、口茧、柴印茧这3种下茧在形状上与正常蚕茧的形状并无太大区别，其瑕疵特征主要体现在蚕茧表面图像信息上，因此需要对蚕茧表面图像信息进行特征提取并识别。国内目前的主要研究方向是通过图像处理技术实现蚕茧质量的检测(唐科梦等，2014)，如刘莫尘等人(2018a)提出的可根据图像中的颜色和面积特征对方格簇蚕茧进行定位分割，并采用FCM(fuzzy C-means)及HSV(hue, saturation, value)模型对方格蔟中的黄斑茧进行检测并剔除(刘莫尘等，2018b)；陈浩等人(2016)提出的辅助选茧方法中采用Matlab对蚕茧表面污迹面积进行计算；赵明岩等人(2016)设计的自动选茧装置通过分析图像中蚕茧的形状、大小、颜色实现对蚕茧品质的判别。由于蚕茧的形状为类椭球体(周志宇等，2006)，无法通过普通相机的一次图像采集得到其完整的表面信息(吴晓军和郭亚雪，2015)，需要对同一蚕茧的不同角度采集图像后再进行特征识别。当污迹横跨分布在两幅图像上时，识别算法无法整体分析该污迹的特征，对蚕茧表面瑕疵信息的识别准确率有一定影响。针对上述问题，目前常用的方法是图像拼接(Wu等，2012)，如姚立健等人(2015)提出基于SIFT(scale invariant feature transform)算子的苹果表面图像拼接方法；汪军等人(2017)针对刑侦领域柱体表面指纹等遗留物证的图像拼接；王玉亮等人(2011)基于SIFT特征的眼底图像拼接等。拼接后的图像包含了非平面物体表面的完整信息，既方便算法的识别，也方便检验人员的复检。但蚕茧的表面是椭球面，其两侧边缘信息在成像面上的投影存在着信息压缩，无法直接对两幅具有重合区域的不同角度下的蚕茧图像进行拼接，需要先对图像进行展开。现有的图像展开方法多是针对半径一致的柱面进行展开，如王伟等人(2014)通过对装置的标定获取圆柱的半径后对柱面上二维条形码的畸变校正，丁超等人(2018)对深孔内壁表面图像的展开，Dai(2010)和Cao等人(2003)对卷曲文档表面图像的展开复原，但单一半径的展开方法在蚕茧这一椭球面的展开应用上效果不佳。因此，本文基于柱面展开理论，结合蚕茧的形状特征，建立蚕茧数学模型的等效阶梯柱面模型，通过对等效阶梯柱面模型表面图像的展开，等效实现蚕茧表面图像的展开。1蚕茧表面图像采集装置设计为获取蚕茧的尺寸信息，使用游标卡尺对某一批次的蚕茧进行实际测量，测得该批次蚕茧的直径区间为16—20 mm，长度区间为28—37 mm，表面遍布茧丝形成的不规则褶皱并附着有少量浮丝，摩擦力较大。根据蚕茧的形状特征，设计了蚕茧表面图像采集装置，如图 1所示。滚轴1和滚轴2通过3个齿轮联动，使其沿同一方向以相同的速度转动。当两根滚轴同步转动时，放置在两滚轴之间的蚕茧在摩擦力的作用下沿与滚轴相反的方向转动，装置上方的相机采集蚕茧表面的图像，相机距离滚轴轴心的距离设置为200 mm。 图1 蚕茧表面图像采集装置 Silkworm cocoon surface image acquisition deviceFig 1((a) schematic; (b)physical picture) 2蚕茧表面图像展开数学模型2.1蚕茧数学模型2.1.1一般蚕茧模型对采集到的蚕茧图像进行预处理：灰度化、滤波、二值化、边缘提取(Cope和Rockett，2000)。根据Fitzgibbon等人(1999)提出的直接最小二乘拟合算法对提取得到的边缘进行椭圆拟合，得到蚕茧图像的长短轴以及旋转角，如图 2所示。 图2 蚕茧图像拟合 Cocoon image fittingFig 2((a) original image; (b) elliptic fitting image) 通过蚕茧的椭圆拟合结果图，可以发现蚕茧轮廓近似为标准椭圆曲线，因此本文假设理想的蚕茧模型为标准椭圆曲线绕其长轴旋转而成的椭球体，如图 3所示。 图3 理想蚕茧模型 The ideal silkworm cocoon modelFig 3定义图 3中椭圆短轴$CD$的长度为蚕茧的直径$2b_{c}$；椭圆长轴$AB$的长度为蚕茧的长度$2a_{c}$。则该椭圆绕长轴$AB$旋转一圈得到的椭球体方程即为理想蚕茧模型方程，即 1 $\frac{{x^{2}}}{{a^{2}_{c}}}+ \frac{{y^{2}+z^{2}}}{{b^{2}_{c}}}=1$ 2.1.2蚕茧等效阶梯柱面模型为了对蚕茧表面图像进行展开，参考周平等人(2010)关于鸡蛋体积与表面积的测试方法，根据微积分的思想原理，提取蚕茧模型某一微区段进行分析，如图 4所示，图中$EFHG$所在平面垂直于短轴的水平面，$EF$、$GH$两个水平面之间的距离为Δ$h$，点$E$′′、$F$′′分别为点$E$、$F$在$GH$水平面上的正投影，取线段$GE$′′的中点为点$G$′，线段$HF$′′的中点为点$H$′，即 图4 蚕茧模型切分示意图 Silkworm cocoon model segmentationFig 4 2 $G′H′= \frac{{1}}{{2}}(GH+EF)$ 当Δ$h$足够小时，椭球$EFHG$部分近似等于圆柱$E′F′H′G′$。在图像中，Δ$h$最小取值为一个像素值。依据上述模型，可以将理想蚕茧模型沿轴线方向切分为$n$段，从上到下依次为第$1, 2, …, n-1, n$段，每一段的高度Δ$h=a_{c}/n$，则理想蚕茧模型等效为一系列圆柱的叠加，由此得到蚕茧的等效阶梯柱面模型如图 5所示。 图5 蚕茧等效阶梯柱面模型转换示意图 Cocoon transformation into equivalent step cylinder modelFig 5此时，蚕茧的表面积$S$可近似转换为阶梯柱面模型的表面积，即 3 $S=Δh· \sum\limits^n_{k=1}2 \mathsf{π}r_{k}+ \mathsf{π}(r^{2}_{1}+r^{2}_{n})$ 式中，$r_{k}$为第$k$段圆柱的半径；$r_{1}$为第$1$段圆柱的半径；$r_{n}$为第$n$段圆柱的半径。假设图 4中的$EFHG$段为蚕茧模型的第$k$段，取$G$点坐标为$(X_{k}, Y_{k})$，$E$点坐标为$(X_{k-1}, Y_{k-1})$，由式(2)得$r_{k}$为 4 $r_{k}=Y_{k-1}+ \frac{{Y_{k}-Y_{k-1}}}{{2}}= \frac{{Y_{k-1}+Y_{k}}}{{2}}$ 因为$E$、$F$、$G$、$H$皆为蚕茧表面上的点，因此其坐标均满足椭圆的一般方程，即 5 $\frac{{ X^{2}_{k}}}{{a^{2}_{c}}}+ \frac{{Y^{2}_{k}}}{{b^{2}_{c}}}=1$ 圆柱从第$1$段到第$n$段自上而下依次排列，则第$k$段的$G$点坐标$X_{k}$为 6 $X_{k}=a_{c}-kΔh$ 联立式(4)—(6)可得$r_{k}$与$X_{k}$的关系为 7 $r_{k}= \frac{{b_{c}}}{{2a_{c}}} (\sqrt{a^{2}_{c}-X^{2}_{k-1}}- \sqrt{a^{2}_{c}-X^{2}_{k}})$ 2.2等效阶梯图像展开对于一个半径为$r$，高为$h$的回转圆柱面，将其沿任意一条素线展开，可以得到一个边长为$2$π$r$，高为$h$的矩形，如图 6所示，该矩形包含了圆柱表面的全部信息。 图6 柱面展开模型 Cylinder unfold modelFig 6设圆柱上$A$点在投影展开面的映射点$A′$为展开图的起点，则对于圆柱面上任意一点$P(x_{p}, y_{p}, z_{p})$在展开图上的坐标$P′(u, v)$为 8 $u=-x_{p}$ 9 $v=l_{PD}=r·{\rm arctan} \frac{{z_{p}}}{{y_{p}}}$ 式中，$l_{PD}$为$PD$两点之间的弧长。因为蚕茧的直径(18 mm)较小，忽略相机小孔成像模型引入的误差(Bors等，1997)，假设相机获得的是蚕茧的平行投影图像。根据图像中蚕茧的像素尺寸，建立阶梯柱面模型坐标系$O_{\rm W}-X_{\rm W}Y_{\rm W}Z_{\rm W}$，成像面坐标系$o-xy$如图 7所示，柱面模型第$k$段圆柱上的圆弧$PC$在成像面上的对应线段为$P″C$。 图7 柱面模型—成像面坐标系 Cylindrical model-imaging plane coordinate systemFig 7设$P″$的坐标为$(x, y)$，$r_{k}$为第$k$段的圆柱在成像面$o-xy$中的半径值，$α$为$O_{W}P$与$Y_{W}$轴负半轴的夹角，即 10 $y=-r_{k}·{\rm cos}α$ 11 $x=-x_{p}$ 12 $α={\rm arctan} \frac{{z_{p}}}{{y_{p}}}$ 联立式(8)—式(12)可得单一圆柱面展开图与原图的坐标点映射关系为 13 $x=u$ 14 $y=r_{k}{\rm cos} \frac{{v}}{{r_{k}}}$ 在图像坐标系$o-xy$中，图像原点为图像左上角的顶点，设原图图像中心坐标为$(x_{c}, y_{c})$，图 7中$o-xy$坐标系沿$y$轴的负方向平移$y_{c}$距离得到的坐标系即为原图像坐标系，则式(14)在原图像坐标系中表达式变换为 15 $y=y_{c}-r_{k}{\rm cos} \frac{{v}}{{r_{k}}}$ 这里采用的等效阶梯圆柱模型，由于圆柱的半径$r_{k}$随着$X$值的变化而变化，等效阶梯圆柱展开图中每一阶梯圆柱的展开起点不一致。设展开图的中心坐标为$(u_{c}, v_{c})$，对于第$k$段圆柱，其在展开图上的起点坐标$v_{k}$为 16 $v_{k}=v_{c}- \frac{{1}}{{2}} \mathsf{π}r_{k}$ 此时，第$k$段圆柱在展开图上的弧长$l$为 17 $l=v-v_{k}$ 则式(10)转换为如下方程 18 $y=y_{c}-r_{k}{\rm cos} \frac{{v-v_{k}}}{{r_{k}}}$ 综上，得到阶梯柱面模型展开图坐标$(u, v)$与原图坐标$(x, y)$的映射关系为 19 $\left\{ \begin{align} &x=u \\ & y=y_{c}-r_{k}{\rm cos} \frac{{v-v_{k}}}{{r_{k} }} \end{align} \right.$ 3等效阶梯柱面模型展开结果3.1实验材料为检验映射公式的有效性，根据蚕茧的平均参数，通过3D打印制作了一个长为33 mm，直径为18 mm的蚕茧等效阶梯柱面模型，该模型阶梯高度为1 mm。模型上贴有大小为1 mm×1 mm的黑白相间棋盘格。实验相机为锐视一声一视—Q4，图像尺寸为800×640像素，棋盘格在图像中的高度Δ$h$=14像素。实验的软件环境为：Windows10 64位系统下的Visual Studio 2013+OpenCV 2.4.11，硬件环境为：8GB内存，Core i5-7300HQ 2.50 GHz CPU。实验开始前运用Matlab的相机标定工具箱对相机的内外参数进行标定，标定所得参数用以消除在后续实验中镜头引入的畸变(Zhang，2000)。3.2实验步骤对采集到的图像进行预处理后，根据椭圆拟合所得的旋转角，对其进行旋转，将阶梯模型的长轴调整为竖直状态，随后以拟合得到的外接矩形信息作为输入进行ROI区域提取，得到阶梯柱面表面图像如图 8所示。 图8 阶梯模型表面图像 Step model surface imageFig 8创建长为$bc$π，高为$a_{c}$的空白展开图，采用线性插值算法进行反向插值，公式为 20 $ \begin{array}{c}f(u, v)=(1-d)f(x, y)+\\d·f(x, y+1)\end{array} $ 21 $d=y-\left\lfloor {y} \right\rfloor $ 式中，$f(u, v)$为展开图中$(u, v)$点的像素值；$f(x, y)$为原图中$(x, y)$点的像素值，$\left\lfloor {\;} \right\rfloor $表示向下取整函数，$d$为$y$与不大于$y$的最大整数的差值，即$y$的小数部分。从上到下依次遍历展开图的像素点，因为$x=u$，所以将$x$代入式(7)即可得到该像素点所在的“半径”值，将$v$代入本文提出映射公式(19)求得该像素点在原图中的对应坐标$y$，将原图中对应像素点的三通道像素值依据式(20)(21)分别赋给展开图，得到的阶梯柱面的变径展开图如图 9所示。 图9 阶梯柱面模型变径展开图 Stepped cylinder model variable diameter unfold diagramFig 9以等效阶梯柱面模型的最大半径段的半径值作为展开半径，对图 8进行柱面展开，得到该模型的等径展开图如图 10所示。 图10 阶梯柱面模型等径展开图 Equal-diameter unfold of the stepped cylinder modelFig 10由于实验所用的方格均为1 mm×1 mm的小方格，在理想状态下，展开后图像中每个方格的大小应该相同，因此以图像中心的白色小格坐标及大小为基准，即可计算得到理想展开图中所有方格角点的理论值。实际上，方格粘贴时存在上下方格对齐误差，故而本文采用图像中心列的小方格作为参考基准。在等径展开图与变径展开图中分别选取32个位置对应的方格，即图 11中用红点标记的方格。 图11 标注后的展开图 Unfolded image after markingFig 11((a)marked equal-diameter unfold image; (b) marked variable-diameter unfold image) 3.3实验结果分析通过对比两种展开方法方格左上角的坐标值与理论值之间的差值，检验本文展开方法的效果。由于展开图只进行了原图$y$轴方向上的到展开图中$v$轴方向上的映射变换，在$x$轴上为平行映射，因此，这里主要对展开图中的$v$坐标值进行分析比较。以方格左上角角点编号为横坐标，角点$v$轴坐标值为纵坐标，将32个方格按坐标值从小到大依次排列，并用线段连接相邻两点，得到的坐标值折线图如图 12中的实线所示。图 12中的红色五角星所在位置为中心基准方格左上角坐标的$v$值。 图12 $v$轴坐标值折线图 Value line chart of $v$-axis coordinateFig 12两种不同展开方法相对于理论坐标值的误差折线如图 13所示。 图13 误差折线图 Error line chartFig 13结合图 12和图 13两幅折线图可以看到，相较于以最大值展开的等径展开坐标折线，由本文提出的方法所得到的坐标值折线与理论值的折线拟合程度更高，误差范围更小。为了进一步分析误差折线，将角点的$u$轴坐标值折线加入分析，两种展开算法的$u$轴坐标值折线图如图 14所示，图中两虚线之间为图像的低纬度区域。 图14 $u$轴坐标值折线图 Value line chart of $u$-axis coordinateFig 14综合角点所在的横纵坐标可以发现，两种展开算法不仅在$v$轴方向上的中心区域具有优良的展开效果，在$u$轴方向上的中心区域即椭球的低纬度区域也有着良好的展开效果，但是在图像边缘及高纬度区域，本文提出的变径展开效果要优于等径展开的效果。3.4蚕茧图像展开实例蚕茧表面图像预处理过程如图 15所示，获取蚕茧表面图像15(a)后，对图像进行灰度化、滤波、二值化处理之后得到图 15(b)，通过Canny算法提取轮廓的边缘，得到图 15(c)，对边缘曲线进行椭圆拟合，得到图 15(d)；根据拟合所得的椭圆外接矩形，对图 15(a)进行ROI区域提取，并进行旋转，得到椭圆长轴保持竖直状态且去除了多余背景区域的图 15(e)。 图15 蚕茧表面图像展开实例 Example of silkworm cocoon surface image unfoldFig 15((a) original image; (b)threshold; (c) edge extraction; (d) elliptic fitting; (e) region of interest extraction) 图像获取装置的电机启动，待蚕茧旋转一定角度后电机停转，相机获得与上一幅图像有重叠区域的另一个角度蚕茧图像，执行同样预处理的步骤，得到该蚕茧另一个角度图像的ROI区域，采用本文展开方法对图像进行展开，分别得到两幅蚕茧表面图像的展开图如图 16所示。 图16 同一蚕茧不同角度下的表面图像展开图 Unfolding of the surface image of the same silkworm cocoon at different anglesFig 16对两幅图像进行手动拼接，得到图 17，图中蚕茧的纹理在拼接图像重叠处及重叠区域边缘的连续性保持较好，表明展开图的拼接效果良好。 图17 蚕茧展开拼接图 Unfolded picture stitchingFig 174结论本文通过对蚕茧边缘轮廓的椭圆拟合，得到蚕茧的数学模型，将该数学模型转化为等效阶梯柱面模型，建立等效阶梯柱面模型表面图像展开的数学模型，从而实现对蚕茧表面图像的展开。利用等效阶梯柱面模型和实际蚕茧图像对本文方法进行验证，实验结果表明，本文所用阶梯柱面展开方法相对于采用单一半径值的展开方法具有更好的展开效果，可为后续蚕茧表面图像的拼接与识别提供良好的环境基础。但与此同时，在坐标折线图中可以观察到实验中存在着除粘贴误差之外的误差，该误差主要是由小孔成像模型引入，因此将在后续的研究中引入这一误差源，改进展开方法以消除这一误差。