基于減法直方圖算法的土壤圖像陰影檢測

曾紹華，王 琪，佘春燕，王 帥，羅達璐

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%，本文算法對低于此陰影占比土壤圖像陰影檢測的適應性尚需進一步實驗驗證，對低陰影占比土壤圖像陰影檢測可能尚需要進一步研究.