基于减法直方图算法的土壤图像阴影检测

曾绍华，王 琪，佘春燕，王 帅，罗达璐

1. 重庆师范大学 计算机与信息科学学院，重庆 401331； 2. 重庆市数字农业服务工程技术研究中心，重庆 401331；3. 重庆市农业技术推广总站，重庆 401147； 4. 重庆市沙坪坝区农业农村委员会，重庆 400030

野外自然环境下机器视觉采集的土壤图像，由于光照的遮挡以及土壤自身表面凹凸不平，导致图像中土壤区域存在阴影，阴影会对土壤识别产生干扰．为了消除阴影过渡区域对后续土种识别子图切割、 土种识别特征提取的影响，需对阴影进行检测分割.

阈值分割是常用的图像分割算法[1-5]．Otsu[6]是一种经典的阈值分割算法，它采用类间方差最大化的优化准则来获取分割阈值，它在类间数据量均衡时分割效果好，在类间数据量差异较大时不能达到好的分割效果．TSAI[7]提出色度和亮度比的Otsu阴影分割方法，常用于街区阴影分割，并将颜色较暗建筑物区域分割为阴影．瑚敏君等[8]利用HSI颜色空间数据定义一个测度(S-I)/(S+I)，根据建筑物阴影检测的经验值指定分割阈值进行阴影检测．MENG Y等[9]提出了减法直方图的阈值分割方法，通过亮度和比率特征分别计算其保留直方图，再分别求它们的最大梯度区间，选择最大梯度区间的F中位值点作为亮度和比率特征的分割阈值点，分割背景部分(非阴影)和目标部分(阴影)； 它解决了“阴影、 非阴影数据量差异较大”和“较暗建筑物区域与阴影”的错分问题．但其减法直方图方法的比率特征描述不适宜土壤图像的阴影和非阴影的区分描述，算法中超参数获取花费了大量的时间，保留直方图的计算区间存在冗余，导致时间复杂度较大．因此，土壤图像阴影检测的减法直方图算法尚需改进.

1 原始减法直方图算法[9]

1.1 亮度和比率特征的定义

亮度(I)和比率(α)特征的定义如下[10-12]：

I(i，j)=1-(R*(i，j)+G*(i，j)+B*(i，j))/(255*N)

(1)

α(i，j)=B*(i，j)/max(R*(i，j)，G*(i，j))

(2)

其中，R*(i，j)，G*(i，j)，B*(i，j)分别表示图像第i行第j列像素RGB三通道的灰度值．N为通道数3.

1.2 原始减法直方图算法

原始减法直方图算法描述如算法1所示．在步骤Step4-3中，最大梯度区间的中位点值：

算法1 原始减法直方图算法

F=(1-C1)×V+C1×P

(3)

V表示ek中最小的非零波谷值，P表示ek最右边峰值，参数C1∈[0，1]文献[9]取的经验值0.4.

1.3 原始减法直方图算法检测土壤图像阴影存在的问题

从包含阴影和非阴影的土壤彩色图像中随机抽取60张构成20组图像，运用公式(1)、 (2)将图像分别转换为I和α的灰度图.

图1显示，土壤图像转换为α灰度图后，阴影与非阴影的区分度较低．其余19组土壤图像有相同的结果．对原始减法直方图算法分析，原始算法保留直方图的计算区间为64等份值域，阈值的搜索存在冗余，有效缩短阈值的搜索区间可降低时间复杂度； 原始算法中参数C1在[0，1]区间内按一定步长进行实验搜索获得经验值0.4，对于土壤图像阴影检测不具科学性，鲁棒性差.

图1 No.15组土壤图像的I和α灰度图

2 改进的减法直方图算法

2.1 重构比率α'特征

为了提升土壤图像阴影和非阴影的区分度，重构α为α'：

α'(i，j)=

(4)

其中，τ为拉伸因子，τ∈[2，3].

2.2 寻找阴影分割阈值搜索区间

原始减法直方图算法中保留直方图的计算区间为64等份值域，存在冗余.

分别对I和α'特征的64等分直方图B0和R'0高斯平滑获得B和R'.

2.2.1 计算相邻频数比r

计算B和R'相邻频数比r

rH(i)=log(H(i+1)/H(i))

s.t.H={B，R'}

(5)

其中，H(i)，H(i+1)分别表示直方图在i和i+1位置的频数.

2.2.2 平滑相邻频数比r

由于高斯平滑处理后的直方图仍可能出现如图2红圈所示的伪谷，为了消除其对rH的影响，对rH进行平滑[13-14]处理，具体平滑方法为

图2 No.15组第3幅图I特征的64等分直方图B0经高斯平滑后的结果B

1) if (rH(i)<0&&rH(i+1)>0&&rH(i+2)<0)

rH(i+1)=-0.1

2) if (rH(i)<0&&rH(i+1)>0&&rH(i+2)>0&&rH(i+3)<0)

rH(i+1)=-0.1，rH(i+2)=-0.1

3) if (rH(i)>0&&rH(i+1)<0&&rH(i+2)>0)

rH(i+1)=0.1

4) if (rH(i)>0&&rH(i+1)<0&&rH(i+2)<0&&rH(i+3)>0)

rH(i+1)=0.1，rH(i+2)=0.1

2.2.3 求平滑后相邻频数比r直方图的波谷向量

平滑后相邻频数比r直方图中，如果rH(i)<0&&rH(i+1)>0，则i+1为波谷点．顺序提取所有波谷点位置构成向量φH，并在向量φH的第一个元素前插入0作为第一个峰的位置起点，在最后一个元素后插入63作为最后一个峰的终点．即φH为平滑后相邻频数比r直方图的完整波谷点向量，其任意相邻点构成区间上平滑后r直方图的波峰.

2.2.4 获得阴影分割阈值的搜索区间

平滑后r直方图的每一个波峰对应于B和R'上的一个真正的波峰．取φH中任意相邻点构成区间范围，分别计算B和R'直方图在该区间的面积：

s.t.H={B，R'}

(6)

其中，φH(j)、φH(j+1)为φH中相邻2个元素的值，H(k)表示B或R'第k个位置的对应频数；w表示直方图64等份中1等份的长度.

2.2.5 获取阴影分割阈值的搜索区间

根据上述算法思想，获取I和α'特征的非阴影和阴影分割阈值搜索区间算法如下所示.

2.3 自适应获取阴影检测阈值

2.3.1 重构保留直方图

原始减法直方图算法理想状态保留直方图(如图3a)，其最大梯度差的梯度区间中位点为F．实际状态保留直方图(如图3b)，以最大梯度差的梯度区间中位点为F，其结果不准确，影响分割效果．为解决上述缺陷，重构保留直方图．运用式(7)对保留直方图ek进行拉伸，增大保留直方图最大梯度区间的梯度差.

图3 保留直方图

(7)

其中，ω为拉伸因子，ω≥0.1，size为特征等份数64，tanh为双曲正切函数.

2.3.2 获取高保留率点F'

将e'k连续5等份分为一组，第1组为[e'k(0)，e'k(4)]，第2组为[e'k(5)，e'k(9)]，…，第13组为[e'k(60)，e'k(63)]，第13组最后一位用0补齐； 求第1组到第13组每一组中e'k(i)的最大值； 在13组中剔除最大值为0的组，从剩余组中找到最大值最小的那一组，其非0最小值所在e'k中的位置即为非零波谷点v'．保留直方图e'k中v'右边的最大值位置为k'．则高保留率点F'为

(8)

其中，MidT表示保留直方图中第i个等份的中位值点，e'k(i)表示保留直方图中第i等份，i在[v'+1，k']范围内.

2.3.3 获取阴影检测阈值

图4 I与α'阈值位置交点

2.3.4 阴影检测阈值搜索算法

根据上述算法思想获得的阴影检测阈值搜索算法如算法3所示.

算法3 阴影检测阈值搜索算法

3 阴影检测阈值搜索加速算法

3.1 阴影检测阈值搜索算法分析

3.2 改进的阴影检测阈值搜索算法思想

算法3阴影检测阈值搜索算法获得分割阈值TB和TR'为2条折线(如图4)的交点．通常情况下，2条折线有1个交点，一条折线的一端仅跨越另一条折线1次，连接每条折线端点的2条直线也相交．因此，将折线相交近似为直线相交，逐步逼近折线交点，可减少搜索次数．其方法是：首次取端点构成直线，获取交点，再用端点与交点的中位点在折线上的位置点替代端点重构直线，获取新的交点，这样反复直到搜索到折线交点.

3.2.1 迭代获取阴影检测阈值

图5 迭代获取阴影检测阈值

(9)

tB∈[tBleft，tBright]，tR'∈[tR'left，tR'right]

将①与③和④进行比较，②与③和④进行比较，若相等停止迭代； 若不相等则重构①②和③④连成2条直线，更新2条直线交点(TB，TR')，更新①②③④点，直到迭代停止．迭代停止时的交点(TB，TR')即为最后的阴影检测阈值.

3.2.2 阴影检测阈值搜索加速算法

3.3 阴影检测阈值搜索加速算法迭代次数及收敛性证明

实际上，设I左右两主峰间距离间隔为disB，α'当前左右两主峰间距离间隔为disR'； 每迭代一次disB减少一半，disR'也减少一半； 迭代总次数为n=┌log2disB┐+┌log2disR'┐(向上取整)．通常，迭代4～8次到达折线交点，获得分割阈值．上述结论也证明加速算法是收敛的.

4 实验结果与分析

4.1 实验样本获取

根据重庆土壤分类与代码[DB50/T796-2017][15]，对重庆璧山区分布的4土属(暗紫泥、 红棕紫泥、 灰棕紫泥、 黄棕紫泥)34土种，在野外自然环境下用土锹锹出耕层0～20 cm的土壤，拍摄土壤(心土)自然断口，使土壤区域位于图像中心位置，并占全图面积50%以上．从已采集的共计342张土壤图像中，随机抽取60张在土壤区域存在部分阴影的土壤彩色图像，用文献[16]分割出土壤区域图像组成20组实验图像样本.

4.2 实验设计

为了验证本文算法的有效性设计实验如下：

1) 阴影检测精度对比实验(实验1)：应用实验样本，进行本文算法2+算法3和本文算法2+算法4，文献[9]原始减法直方图算法(包括原始算法eSH和加速算法iSH算法)、 文献[7]算法、 文献[8]算法的对比实验，检测土壤彩色图像阴影，检验本文算法阴影检测精度，验证本文算法2+算法4计算的实际迭代次数与是否在理论(证明)最大迭代次数内的一致性问题.

算法2 寻找I和α'特征的非阴影和阴影分割阈值搜索区间算法

算法4 阴影检测阈值搜索加速算法

2) 阴影检测效率对比实验(实验2)：应用实验样本，进行本文算法2+算法3和本文算法2+算法4，文献[9]原始减法直方图算法(包括原始算法eSH和加速算法iSH算法)、 文献[7]算法、 文献[8]算法的对比实验，检测土壤彩色图像阴影，检验本文算法阴影检测效率.

3) 全图为阴影或非阴影检测实验(实验3)：取土壤图像子块全图为阴影或非阴影(包括被人工用非阴影部分替换阴影部分像素的人工合成子图)，验证本文算法2对于全图为阴影或非阴影图像检测的有效性.

4.3 实验环境

本文实验在Intel(R) Xeon(R) Silver 4114CPU @ 2.20GHz 2.19 GHz (2处理器)，内存64GB，显卡NVIDIA TITAN V，Windows 10专业工作站版64位，Visual Studio 2017，OpenCV 3.4.1环境下完成(未使用多CPU、 GPU并行加速).

4.4 仿真实验结果及分析

4.4.1 阴影检测对比实验(实验1和实验2)结果

1) 阴影检测实验图像结果

对20组实验图像样本进行本文算法2+算法3和本文算法2+算法4、 文献[9]原始减法直方图算法(包括原始算法eSH和加速算法iSH算法)、 文献[7]算法、 文献[8]算法的阴影检测对比实验．随机选取No.3和No.15组实验样本，阴影检测结果如图6所示.

注：本文算法实验参数：τ=2.5，ω=1(对参数τ在[2，3]范围内以0.1为步长进行搜索实验，τ=2.5时分割效果最好； 对参数ω在[0，3]范围内以0.1为步长搜索实验，ω≥0.1均有效)．No.3组中文献[9]eSH算法C1实验参数分别为：0.6，0.65，0.5； No.3组中文献[9]iSH算法C1实验参数分别为：0.6，0.45，0.4； No.15组中文献[9]eSH算法C1实验参数分别为：0.5，0.7，0.8； No.15组中文献[9]iSH算法C1实验参数分别为：0.5，0.55，0.7．图6 实验图像(No.3和No.15组)阴影检测图像结果

2) 阴影检测实验数据结果

上述实验的阴影检测——分割阴影与非阴影精度(用亮度标准差评估)数据如表1所示.

表1 实验样本分割后阴影与非阴影的亮度标准差及迭代次数

续表1

上述实验的阴影检测执行10次，平均时间花销如表2所示.

表2 阴影检测实验的时间花销 s

4.4.2 阴影检测实验结果分析

1) 重构α'特征的有效性分析

阴影检测实验图像结果图6(No.3组和No.15组样本)b，c，f，g显示：本文算法2+算法3和算法2+算法4使用重构的α'特征阴影检测，相比于文献[9]的eSH算法和iSH算法使用α特征阴影检测，本文改进的算法分割效果更好．表1显示：本文算法2+算法3和算法2+算法4 对20组样本阴影检测，分割的阴影与非阴影的亮度标准差均值分别为23.987 1，24.542 8和24.713 7，24.200 2； 文献[9]eSH算法和iSH算法对20组样本阴影检测，分割的阴影与非阴影的亮度标准差均值分别为15.996 0，51.773 1和21.091 9，52.367 9．本文改进算法分割的非阴影亮度标准差值远小于文献[9]eSH算法和iSH算法分割的非阴影亮度标准差值，且其分割的阴影与非阴影的亮度标准差均值也远小于文献[9]的算法，说明文献[9]算法阴影分割存在不足，没有从非阴影中分割出来，本文改进算法的阴影检测精度更高．上述结果证明，重构α'特征对土壤图像阴影检测是有效的.

2) 阴影检测精度分析

阴影检测实验图像结果图6(No.3组和No.15组样本)d，e，f，g显示：文本算法2+算法3和算法2+算法4对20组样本阴影检测，相比于文献[7]算法和文献[8]算法，本文改进算法分割效果更好．表1显示：文献[7]算法和文献[8]对20组样本阴影检测，分割的阴影与非阴影的亮度标准差均值分别为18.036 2，41.665 5和20.851 7，33.930 6．本文改进算法分割的非阴影亮度标准差值远小于文献[7-9]的对比算法分割的非阴影亮度标准差值，且其分割的阴影与非阴影的亮度标准差均值之和也都远小于对比算法，说明对比算法阴影分割不足，部分阴影没有从非阴影中分割出来，本文改进算法阴影检测精度更高.

3) 本文算法2+算法4的实际迭代次数结果分析

表1显示：本文算法2+算法4对20组样本阴影检测的阈值搜索迭代次数为4～6次，本文3.3节理论证明阴影检测阈值搜索最大迭代次数为6次，实验结果证明算法计算的实际迭代次数在理论(证明)最大迭代次数内，与理论证明结果是一致的.

4) 2条分割阈值折线交点数的分析与处理

土壤图像可分为全图为阴影、 全图为非阴影和包含部分阴影三类图像．算法2可检测出全图为阴影或非阴影图像，不需继续阴影检测．阴影检测是对包含部分阴影图像检测，是一个典型的二分类问题，通常2条分割阈值折线交点数为1．由于受到土壤杂质等因素影响，极少数出现2条分割阈值折线交点数大于等于2的多个交点情形，选择其亮度标准差均值最小的交点分割阈值，实验结果证明是有效的.

5) 阴影检测时间效率分析

表2显示：对20组样本执行10次阴影检测任务，本文算法2+算法3和算法2+算法4的平均时间花销分别为42.687±2.405、 14.050±1.540； 文献[9]eSH算法和iSH算法分别为61.016±3.171、 4.718±0.570； 文献[7]算法为1.897±0.169； 文献[8]算法为1.142±0.129．上述数据显示：本文原始算法——算法2+算法3比文献[9]原始算法——eSH算法快近1/3，文献[9]加速算法——iSH算法比本文加速算法——算法2+算法4快近3倍是因为iSH算法完全依赖于数据集计算的阴影检测经验值，在经验值附近搜索降低运行时间开销，而本文算法2+算法4是不需要任何先验知识的完全自适应单张图像阴影检测算法，是不需要任何人为参数的更先进的自适应算法．文献[7]算法是自定义单测度的Otsu算法和文献[8]是利用HSI颜色空间数据定义一个测度(S-I)/(S+I)，根据建筑物阴影检测的经验值指定分割阈值进行阴影检测，尽管它们的耗时都非常低，但它们主要是用于建筑物阴影检测的，不适用于土壤图像阴影检测.

4.4.3 全图为阴影或非阴影检测实验(实验3)结果与分析

对20组土壤图像实验样本中取全图为阴影或非阴影的子块(如图7)．在阴影和包含少量阴影的非阴影区域各取80*50像素大小的子块，对A子块中阴影像素用同子块中非阴影像素替换，人工合成非阴影子块A和阴影子块B作为检测样本，用算法2检测A、 B子块的I和α'特征直方图.

图7显示：非阴影子块A和阴影子块B的I和α'特征直方图呈现单峰结构．如果检测为单峰结构说明土壤图像为全图非阴影土壤图像或全阴影土壤图像，无需进一步分割，可直接进入下一步切割为用于土种识别的子图，进行子图归一化和土种识别操作.

图7 全图为阴影或非阴影检测实验图像结果

20组土壤图像实验样本中，对其他所有图像生成的非阴影和阴影子块检测，具有相同的检测结果.

5 结 论

本文对土壤图像阴影检测的典型二分类问题和基于I、α特征阴影检测的原始减法直方图算法进行研究，依据典型二分类问题分类特征的双峰特性，改进减法直方图算法用于土壤图像阴影检测.

1) 改进减法直方图算法将α特征重构为α'特征，提升该特征的双峰特性； 改进减法直方图算法对获得的保留直方图进行拉伸，增大保留直方图的最大梯度差，使搜索保留直方图的最大梯度求解阈值变得更直观.

2) 改进减法直方图算法从左右峰值点开始搜索分割阈值，减少了搜索范围，比改进前减少了近1/3计算时间开销．理论证明改进减法直方图加速算法，最坏情形计算保留直方图和F'k的次数为12次，即减法直方图的上确界为12； 通常情况下，迭代总次数为n=┌log2disB┐+┌log2disR'┐．它也证明改进减法直方图算法是收敛的.

3) 改进减法直方图加速算法改进了文献[9] iSH加速算法收敛速度对先验知识的依赖，成为更先进的不需要任何先验知识的自适应算法.

4) 通过重构α'特征，相关算法对比和直方图单峰检测实验，结果显示：算法是有效的.

本文20组实验土壤图像样本中，最小的阴影占比为2.386%，本文算法对低于此阴影占比土壤图像阴影检测的适应性尚需进一步实验验证，对低阴影占比土壤图像阴影检测可能尚需要进一步研究.