基于多轮廓模型的红枣体积和表面积在线测量方法

吴明清，弋晓康，罗华平，李传峰，唐晓燕，陈坤杰

(1.南京农业大学工学院，南京210031；2.塔里木大学机械电气化工程学院，阿拉尔843300；3.江苏省农业农村厅绿色食品办公室，南京210036)

0 引 言

新疆因光照丰富，降雨量少，昼夜温差大的独特环境而出产优质红枣[1]。近年来，新疆红枣产量随种植面积逐年扩大，红枣已经成为新疆特色农业经济的支柱产业，红枣产业的快速发展对红枣采后加工的技术和设备提出了更高的要求[2]。分级是红枣采收后贮藏、加工及销售过程中的重要环节，同时也是提高红枣质量和商品化价值的重要手段。目前红枣主要的分级方法有人工、机械和机器视觉技术。人工分级的效率和精度都比较低；机械分级克服了人工分级缺点，极大的提高了分级效率和质量，但机械分级只能根据红枣的长短径进行分级，分级指标比较单一，不能够真实反映红枣的形状[3]。形状是影响红枣等级和价格的重要因素。机器视觉技术利用图像传感器，获取红枣的大小、形状、颜色、纹理和表面缺陷等信息，实现红枣的分级[4-15]。新疆干制红枣的品种有骏枣和灰枣，其中骏枣长径、短径和长度各不相同。由于二维图像只能提供红枣的两维信息，而且还因视角变化，造成几何形状和投影面积的差异[16]，会降低红枣的分级精度。

体积和表面积都是农产品重要的物理特征，准确获取体积和表面积对农产品的大小分级有重要的意义[17]。但传统类球体体积的测量多采用排水法测量，表面积采用削皮或切片的方法进行测量[18]。效率低、无法实现实时测量，无法用于农产品的在线分级。近年来，已有不少文献研究不规则农产品的体积和表面积的测量方法。李景彬等[10]将红枣近似成类圆柱体，提取红枣二维图像特征参数得出类圆柱体体积，建立图像体积与质量的线性回归方程。Clayton等[19]分别以球体和椭球体为模型估算4种苹果的表面积，采用削皮测量苹果的表面积，两种模型分别低估了实际表面积15%和18%。Uyar 等[20]通过三维激光扫描仪获取农产品三维数字图像的体积和表面积，分别计算鸡蛋和球体的体积和表面积，其中球体体积和表面积的计算误差小于1%，鸡蛋测体积计算误差在1%～3.15%，鸡蛋测表面积计算误差小于1%。熊妮娜等[21]应用三维激光扫描仪系统，通过扫描获取单株树木的三维空间的点云数据，将树冠近似为多个圆台体，求他们的体积之和来计算树冠体积。龚爱平等[22]利用三维线框模型图像处理方法测量不规则球类形状的农产品的体积和表面积，用所提及方法、排水法和扫描图像处理法对不规则的柑橘、苹果和梨进行体积和表面积测量，结果显示3 种方法的相关系数大于0.98；在测量精度方面，当图像数量大于10幅时，测量精度大于98%，测量数量小于7幅时则精度下降明显，测量系统能在5 s内完成被测物体的体积和表面积的测量。但是，目前国内针对红枣体积和表面积检测的研究几乎处于空白，根据体积进行红枣分级的研究还没有相关报道。

本文搭建红枣图像采集的简易装置，利用图像处理技术进行红枣图像的实时采集和特征提取，再通过特征参数构建红枣的三维模型并计算模型的体积和表面积，为将来开发基于图像技术和体积参数的红枣分级设备提供理论和技术基础。

1 测量原理

1.1 二维轮廓与三维轮廓空间的映射推导

三维轮廓模型是表示物体几何形体的模型，是计算机中对物体的形状可视化表现的一种方式，模型由对象轮廓的边缘点按映射关系转换得到，其特点是存储量小、速度快、操作灵活和处理时间较短[22]。要建立被测物体的三维轮廓模型，需要对旋转过程的物体规定时间内按等时间间隔测量采集N 副图像，完成等间隔角度图像的采集，对采集的到图像进行处理，提取图像边缘轮廓点的坐标，按以下原理进行映射转换。

在二维空间，假设向量p，q 围绕原点旋转，角度为θ，顺时针为正，逆时针为负。图1 展示向量p，q 绕原点旋转，得到新向量p'，q'，构造矩阵如式(1)所示：

图1二维向量原点旋转图Fig.1 Two-dimensional vector origin rotation diagram

图像轮廓点的坐标点围绕着坐标系z 轴进行旋转，同一点旋转前后的关系式如(2)所示：

式中（x，y）为旋转前图像坐标的位置；θ 为图像的旋转角，(°)；（x'，y'）为旋转后图像坐标的位置。设坐标系(xoy)内建立空间直角坐标系xyz，空间直角坐标系xyz 转动后形成新的直角坐标系x'y'z'，如图2 所示，转换公式如(3)所示。

图2 三维旋转坐标系Fig.2 Three-dimensional rotational coordinates

式中（x，y，z）T为空间直角坐标系中任意一点；（x'，y'，z'）T为空间直角坐标系x'y'z'中的坐标，(△x，△y，△z)T为空间直角坐标系xyz 到空间直角坐标系x'y'z'的平移参数。其中Rx(φ)，Ry(α)，Rz(θ)为空间直角坐标系xyz到空间直角坐标系x'y'z'的旋转矩阵，φ,α,θ 表示红枣轮廓将其围绕x,y,z,轴旋转角度(°)，如(4)、(5)、(6)所示：

在公式(4)中，假设当坐标点围绕x 轴进行旋转，物体不发生位移，(△x，△y，△z)T=(0，0，0)T，其中(x，y，z)T表示坐标点在平面直角坐标系zoy 上的坐标，(x'，y'，z')T表示轮廓上的像素点在平面直角坐标系z'o'y'上的坐标。

1.2 多轮廓体积和表面积的测量原理

多轮廓三维点云模型如图3 所示，两个相邻轮廓的间角为△θ，轮廓个数用N，轮廓的个数N=180°/△θ，假设轮廓个数N=20，则轮廓在z轴x心上的投影如图3a、b所示，ri和ri+1为投影面轮廓边缘点和轮廓轴心点之间的距离，投影面轮廓边缘相邻两点A，P 和轴心点Z0构成三角形A=△Z0PA，相邻轮廓点和轴心点Z0形成的投影面积可以用式(7)表示：

式中S 为轮廓点的投影面积，cm2，Ai，Ai+1，Ai+2…，Ai+n为投影面轮廓边缘相邻两点和轴心Z0点组成三角形的面积，cm2；n 为轮廓点形成三角的个数。在图3c 中，在Z 轴方向，相邻轮廓投影面形成一个椭圆台，多轮廓模型的体积可以用式(8)表示。

式中V 为多轮廓模型的近似体积，cm3；ab 为多轮廓模型纵轴的长度，cm；Si，Si+1为第i，i+1 层投影面的面积,cm2；△h为相邻投影面的高度，cm。

在图3b中，轮廓投影面的边缘点的周长可以用式(9)表示。

式中C为投影面积的周长；Li，Li+1，Li+2…，Li+n为轮廓边缘相邻两点的距离，cm，n 为边缘轮廓点的个数。在图3c 中，在z 轴方向，相邻轮廓投影面形成一个椭圆台，多轮廓模型的近似表面积可以用式(10)表示。

式中P为多轮廓模型的近似表面积，cm2；Ci，Ci+1为第i，i+1层投影面轮廓周长，cm，△h 为相邻投影面的高度，近似等于椭圆台的母线的长度，cm。

图3体积和表面积计算示意图Fig.3 Schematic diagram of volume and surface area calculation

1.3 图像的分割和点云处理

本文在VS2010 开发平台下采用VC++编程语言，结合PCL(Point Cloud Library1.6 )点云库和OpenCV 跨平台计算机视觉函数库(Open Source Computer Vision Library3.0)对图像进行处理。图像分割的目是将目标从图像中精确的分割出来，对提取的结果进行滤波去噪，提高测量的精度[23]。原始图像如4a 所示。相机拍摄过程受到白色背景，灯光的影响，图像存在噪音，需要对灰度化图像进行滤波处理，本文应用加权平均法将彩色图像转换为单色图像方法进行灰度化[24]，计算公式为Gray=0.299R+0.587G+0.114B，式中R、G、B 分别代表彩色图像3 个基本单色光红色、绿色、蓝色的强度等级，其取值范围均为0～255，灰度化效果如图4b 所示。均值滤波法[25]是一种把相邻像素的相应分量值的平均值作为中心点像素相应分量值的方法，本文采用均值滤波(5×5算子)的方法对红枣图像进行滤波操作，滤波效果如图4c 所示。图像的二值化就是将图像的像素点灰度值设置为0 或255，将整幅图呈现出明显的黑白效果，本文采用Otsu 算法，该算法进行二值化具有简单，时效性高，能够兼顾边缘细节和分割效果的优点[26]，二值化效果如图4d 所示。二值化分割后，图像的边缘存在不连续现象，影响边缘轮廓的提取，应用开运算法使图像边缘光滑，应用腐蚀算法增加连通区域，效果如图4e f 所示。轮廓边缘坐标的提取方法分为两步：1)首先应用Opencv 中的findContours()函数对图像轮廓进行检查，检索模式为RETR_EXTERNAL，该模式只检索最外层轮廓；模式标识符选CHAIN_APP_ROX_NONE，其作用获取轮廓边缘每一个像素的坐标位置；2)应用Opencv中的drawContours()函数绘制外轮廓，函数中参数thickness 来控制轮廓线的粗细，一般设置为1，效果如图4g所示。轮廓边缘像素点是二维坐标，需要将二维坐标转换为三维坐标[26]；具体操作如下，构建一个三维结构体CvPoint3D32f(float x，float y，float z)，将二维坐标存储到三维结构体中；再将三维结构体CvPoint3D32f 的数据类型转化为PointXYZ 点云类型，存储为.pcd 文件格式[27]，效果如图4h 所示。由于图像采集过程，不同间隔角度的轮廓坐标都有一定的差异，点云的坐标系不同，需要对点云进行平移到统一的坐标系中，即采用getMin-Max3D 最小包围盒函数计算轮廓坐标中心(x，y，z)，轮廓点坐标中心(x，y，z)与原点坐标(0，0，0)位置相减就实现轮廓平移，效果如图4i 所示。轮廓以Z 轴为旋转中心按二维轮廓与三维轮廓空间的映射关系转换为多轮廓模型[28]，效果如图4j 所示，模型中点云十分密集，存在冗余，无法有效运算，降低了轮廓模型体积和表面积的计算速度，需要对多轮廓点云模型进行下采样操作。本文采用体素栅格[29]的点云下采样方法，设置素网格立方体的大小为0.1 cm×0.1 cm×0.1 cm，如图4k 所示。

图4图像分割和处理过程Fig.4 Image segmentation and processing

2 试验材料与装置

2.1 试验样品

试验样本为阿克苏骏枣，购买于2018年10月，采用排水法测量红枣的体积，将红枣按体积大小划分成5 个等级，如表1 所示，每个等级取20 个红枣，共计100 个红枣样本备用。

2.2 试验仪器和设备

电子数显游标卡尺（中国桂林量具刃具有限责任公司，精度：0.01mm），直嘴溢水烧杯，量筒（量程：100 mL，精度：1 mL），锡纸，剪刀，记号笔，密封袋，标签纸，镊子等。

2.3 试验平台的设计

自主开发一套机器视觉图像采集装置，如图5 所示，试验装置由工业相机，计算机，旋转台，光源，暗箱等组成。相机采用维视MV-EM510C 型工业相机及4.0 mm 焦距镜头，像素尺寸3.45μm×2.2 μm，帧率15 帧/s，曝光时间为30～5 000 000 μs，以太网与计算机相连，相机通过云台夹固定在暗箱边框，水平朝向旋转台中心，距离被测物体约300 mm；照明光源采用两个LED 灯；旋转台采用宝康隆NA250 电动旋转台，转速60 s/圈，旋转角度360°，旋转台半径为125 mm，在旋转台的中心固定一根长度为160 mm 的钢针，钢针表面涂为白色。

图5图像采集装置Fig.5 Image acquisition device

2.4 红枣图像的采集方法

用户界面用MFC(Microsoft Foundation Class)进行界面设计，该用户操作采用软触发，定时器控制，在红枣的背后放置白板，减少复杂环境的干扰。镜头到红枣的距离150 mm，主光轴高度300 mm，在图像采集时，首先将红枣样本果蒂部位刺插在针尖上，启动旋转按钮，调试光源和焦距，当样本清晰呈现在相机镜头正前方，触发相机。30 s内分别得到样本的图像的帧数为12，15，20，30，45；间隔角度分别为15°，12°，9°，6°，4°；采集角度范围为0°～180°，连续采集图像，图像的存储格式为JPG。

3 结果与分析

3.1 排水法和图像法的实际值测量及标定

本文采用排水法测量红枣体积，操作过程如下：将蒸馏水倒入直嘴溢口烧杯，直至直口嘴溢口烧杯不漏水为止；用镊子夹住红枣浸入烧杯中使水溢出到100 mL量筒中，读出量筒中水的体积，即红枣的体积，每个样本测量3 次，取平均值；再用平均值减去镊子的体积即为红枣的体积。本文图像的处理的方法是用锡纸包裹红枣表面，用剪刀剪切锡纸包裹的多余部分，将锡纸剥离压平放置相机镜头的黑色A4 纸上，采集俯视图如图6所示，图中白色区域为红枣表面锡纸图像。在图像采集过程中，相机实际成像与理想成像之间存在非线性光学变形，为了减少系统误差，提高系统测量的精度，需要对系统进行标定。本文参照[30-31]提出的标定方法，提取图像的特征向量为像素单位，将其转换为图像的实际值。在系统硬件参数、物距确定的情况下，图像的物理尺寸和像素尺寸比值固定不变[32]，如式(11)所示。

图6红枣表面锡纸图像采集Fig.6 Tin paper on red jujube surface image acquisition

式中A为比例系数，Lr为物理尺寸，Lp为像素尺寸；实际面积=0.005 6×像素尺寸，实际体积=0.000 42×像素尺寸，应用Matlab中的bwarea函数统计白色区域面积。

3.2 不同等级红枣的体积和表面积统计分析

将红枣按体积大小划分为5 个等级，每个等级20 个红枣，计算各等级的最小值、平均值、最大值、标准偏差和变异系数。不同等级红枣的实际体积和表面积统计分析结果如表1所示，体积的标准偏差在1.09～2.05之间，变异系数在6%～12%之间；表面积的标准偏差在2.56～2.88 之间，变异系数在4%～10%之间。

表1 不同等级红枣实际值统计分析Table 1 Different grades of red jujube actual value statistical analysis

3.3 球体体积与表面积测量

球体的图像如图7a所示；提取球体的边缘轮廓如图7b所示；将边缘轮廓像素点坐标位置转换为三维点云，如图7c所示；将球体边缘轮廓的中心偏移到原点位置，如图7d所示；按映射关系转换为多轮廓模型，如图7e所示；将多轮廓模型进行下采样，图7f 所示；对多轮廓模型进行垂直Z轴中心点进行投影，如图7g所示；提取投影面的边缘点，如图7h所示。

图7 球体图像和点云处理过程Fig.7 Process processing of spherical image and point cloud

球体样本的的直径为4.2，3.4，3.0，2.8和2.4 cm，按等间隔角度采集图像，间隔角度设置为15°、12°、9°、6°和4°，对应采样图像帧数12、15、20、30和45，采集范围为0°～180°。

3.4 球体投影面高度、间角和直径的误差分析

在表2 中，设球体直径为3.4 cm 固定值，轮廓间隔角度为15°、12°、9°、6°和4°和投影高度为0.1，0.2，0.3，0.4和0.5 cm。为分析投影高度和轮廓间角度对多轮廓球体模型体积和表面积的计算结果有影响。根据表2，找出多轮廓模型实际值和测量值之间相对误差的最大值和最小值，来设定合适的投影高度和轮廓间角度，当轮廓间角度为4°，投影高度为0.1 cm 时，体积的最小的相对误差为6.0%；当轮廓间角为15°时，表面积最小相对误差为1.0%。在本试验中，多轮廓球体模型的体积随轮廓间角和投影高度的相对误差的增大而增大，表面积的相对误差随轮廓间角和投影高度的增大而减小；随球体模型轮廓间角的增大，耗时和内存随之减少，间隔角度4°、6°、9°、12°和15°相对应的消耗时间为180，120，100，60和30 s。

表2 不同投影高度和间隔角度对球体模型的测量误差Table 2 Measurement error of spherical model with different projection height and contour angle

为分析不同直径对球体模型体积和表面积的计算结果的影响，在表3 中，轮廓间角和投影高度为固定值，分别为9°和0.1、0.4 cm，直径为4.2，3.4，3.0，2.8和2.4 cm。观察表3，多轮廓体积模型体积和表面积的相对误差随直径增大变化不明显，但球体直径越小，误差越大；计算多轮廓球体模型体积和表面积实际值和测量值之间相对误差的均值，分别为9.1%和4.34%。

3.5 不同等级红枣的误差分析

通过对多轮廓球体模型的分析，分别确定红枣模型的投影高度分别为0.1 cm 和0.5 cm，轮廓间角为9°，每个等级红枣数量为20 个。计算不同等级红枣模型的体积和表面积。通过观察表4，多轮廓红枣模型体积和表面积的均方根误差随着等级变化不明显；体积平均相对误差随等级增大而增大，表面积的平均相对误差随着等级变化不明显；不同等级红枣体积模型的均方根误差和平均相对误差的均值分别为2.45 cm3和10.02%；表面积均方根误差和平均相对误差的均值为3.65 cm2和7.09%。

表3 不同直径球体的测量误差Table 3 Measurement errors of spheres with different diameter

表4 不同等级红枣体积和表面积的测量误差Table 4 Measurement error of volume and surface area of different grades of red jujube

4 结论

本文搭建图像采集的简易装置，按不同角度采集红枣的图像，编写图像处理软件，提取红枣二维图像轮廓构建多轮廓三维模型，分析在不同轮廓角度，不同投影高度和不同直径下球体和红枣三维模型的体积和表面积。

1）当球体直径为3.4 cm 时，球体的轮廓间角和轮廓投影高度分别为15°，12°，9°，6°，4°和0.1，0.2，0.3，0.4，0.5 cm 时，球体体积的测量值随轮廓间角和投影高度增大增大，而表面积的测量值轮廓间角和投影高度增大减小，体积和表面积的最大相对误差分别为6.0%和1.0%。

2）当球体轮廓间角为9°，球体体积和表面积的投影高度为0.1 cm 和0.4 cm，直径为4.2，3.4，3.0，2.8，24 cm时，不同球体的体积和表面积随直径变化不明显，但直径越小，体积和表面积的误差越大，体积和表面积的相对误差的均值为9.1%和4.34%。

3）当红枣轮廓间角为9°，体积和表面积的投影高度为0.1 cm和0.5 cm时，红枣的体积测量平均相对误差随等级的增大而增大，表面积随等级变化不明显，体积的均方根误差和平均相对误差的均值为2.45 cm3和10.2%；表面积的的均方根误差和平均相对误差的均值为3.65 cm2和7.09%。

本文采用多轮廓模型测量红枣的体积和表面积，判断红枣过程有一定偏差，需要改进算法来提高测量的精度。上述结论对测量红枣的体积的表面积的测量提供一种无损的测量方法，为红枣三维特征的分级装备开发提供技术参考。