生長(zhǎng)輪掃描圖Canny邊界檢測(cè)結(jié)果修復(fù)算法

林啟招，孫永科，何 鑫，秦 磊，邱 堅(jiān)

（西南林業(yè)大學(xué) 材料科學(xué)與工程學(xué)院，云南 昆明 650224）

樹木年代學(xué)是由美國(guó)亞利桑那大學(xué)樹木年輪研究試驗(yàn)室的創(chuàng)始人、天文學(xué)家Douglass在20世紀(jì)上半葉發(fā)展起來的[1]。與之相關(guān)的研究涉及氣候環(huán)境等領(lǐng)域[2]。在研究中，學(xué)者們需要測(cè)算樹木的生長(zhǎng)輪，于是提出了基于圖像處理技術(shù)的GIS工具、R軟件、圖像處理技術(shù)等方法和手段得到生長(zhǎng)輪圖像[3-6]。

數(shù)字圖像處理技術(shù)被用于眾多實(shí)驗(yàn)教學(xué)中[7]，其可以輔助人們快速測(cè)算生長(zhǎng)輪，但仍有不足。為了有效識(shí)別橫切面微觀圖，Wang等[8]提出了自動(dòng)識(shí)別木材微觀圖生長(zhǎng)輪的方法，但未給出修復(fù)橫切面圖像上斷裂斷開生長(zhǎng)輪界進(jìn)的詳細(xì)步驟。王燕鳳等[9]提出了年輪圖像雙邊濾波增強(qiáng)方法，但未提及如何解決不理想的生長(zhǎng)輪界問題。Vaz等[10]提出的生長(zhǎng)輪自動(dòng)識(shí)別與測(cè)量系統(tǒng)，未考慮到生長(zhǎng)輪界不完整的情況。Subah等[11]提出了生長(zhǎng)輪交互式視覺分析系統(tǒng)，在遇到生長(zhǎng)輪斷裂或不完整時(shí)，需通過用戶給定起點(diǎn)和終點(diǎn)進(jìn)行生長(zhǎng)輪修復(fù)，工作量大的情況下，用戶容易疲勞。Fabijańska等[12]在生長(zhǎng)輪自動(dòng)檢測(cè)研究中提出了擬合線條，具有一定弧形的輪界很可能被擬合成理想化的線段。寧霄等[13]提出的基于U-Net卷積神經(jīng)網(wǎng)絡(luò)對(duì)年輪圖像進(jìn)行分割，可達(dá)到好的分割效果，但它的不足是需要進(jìn)行大量的圖像訓(xùn)練學(xué)習(xí)。

對(duì)具有明顯生長(zhǎng)輪界的木材端面掃描圖進(jìn)行Canny邊緣檢測(cè)，可得到生長(zhǎng)輪界，但所得結(jié)果經(jīng)常會(huì)出現(xiàn)邊緣斷裂，即生長(zhǎng)輪界不連續(xù)的情況。Canny檢測(cè)出現(xiàn)邊緣斷裂是由邊緣附近的灰度值接近引起的[14]，被學(xué)者們改進(jìn)后加以運(yùn)用[15]。木材端面掃描圖較自然界宏觀圖像而言，修復(fù)斷裂生長(zhǎng)輪界的要求更高，現(xiàn)有的圖像修復(fù)方法無法用來修復(fù)生長(zhǎng)輪界圖像，例如：舒彬等[16]所得的圖像修復(fù)算法可以修復(fù)自然圖像的小區(qū)域破損，劉昱等[17]提出的圖像修復(fù)算法適用于補(bǔ)全大面積破損的圖像。

鑒于此，基于在木材圓盤上對(duì)若干不同方向進(jìn)行傳統(tǒng)人工生長(zhǎng)輪計(jì)數(shù)思路，本文提出一種對(duì)生長(zhǎng)輪掃描圖Canny邊緣檢測(cè)結(jié)果進(jìn)行修復(fù)的算法，對(duì)作者“基于交互模式和圖像處理的針葉材生長(zhǎng)輪測(cè)算方法”[18]的前期工作進(jìn)行深入研究。在前期工作中，發(fā)現(xiàn)輪界不連續(xù)的情況包括中間斷裂和端頭斷開，如圖1所示，圖 1(a)和圖 1(b)是針葉樹翠柏?cái)嗝鏅z測(cè)前后對(duì)照?qǐng)D，圖 1(c)和圖 1(d)是闊葉樹冬青斷面檢測(cè)前后對(duì)照?qǐng)D。從上往下看，圖 1(b)的第一條輪界線是連續(xù)的、第二條端頭斷開、第三和第四條中間斷裂；圖 1(d)有多處斷裂的生長(zhǎng)輪界線。

圖1 木材端面邊緣檢測(cè)前后對(duì)比圖

1 材料與方法

1.1 材料

本研究試驗(yàn)對(duì)象包括杉木、翠柏、北美紅杉、冬青從髓心向樹皮方向的掃描圖邊緣檢測(cè)結(jié)果，其中前3個(gè)樣品為針葉樹、第4個(gè)樣品為闊葉樹，它們的邊緣檢測(cè)圖都存在生長(zhǎng)輪界斷裂的情況。

1.2 方法

1.2.1 準(zhǔn)備工作

在對(duì)生長(zhǎng)輪界進(jìn)行修復(fù)前，需導(dǎo)入生長(zhǎng)輪掃描圖的Canny邊緣檢測(cè)結(jié)果。應(yīng)對(duì)所獲得的邊緣檢測(cè)圖進(jìn)行降噪處理，例如：給定閾值，把寬度小于閾值的白色輪界區(qū)域轉(zhuǎn)為黑色，即可達(dá)到降噪的目的，作者已在前期工作[18]中進(jìn)行研究并提出了解決方法。

1.2.2 算法流程

已經(jīng)去噪的生長(zhǎng)輪掃描圖Canny邊緣檢測(cè)結(jié)果矩形圖，是本文修復(fù)算法的輸入。算法流程圖如圖 1所示。流程涉及的關(guān)鍵過程和技術(shù)在后續(xù)章節(jié)進(jìn)行陳述。

圖2 生長(zhǎng)輪界修復(fù)流程圖

1.2.3 算法描述

本文設(shè)計(jì)的算法涉及的主要過程包括：①輸入圖像；②輪界區(qū)域排序；③取一斷裂輪界區(qū)域；④找下一個(gè)與之相鄰且中心點(diǎn)X坐標(biāo)值差小于閾值t1的斷裂輪界區(qū)域；⑤連接修復(fù)；⑥延長(zhǎng)修復(fù)。具體處理過程描述如下：

（1）輸入圖像。算法導(dǎo)入已經(jīng)去噪的生長(zhǎng)輪掃描圖的邊緣檢測(cè)二值圖像矩形圖。該過程應(yīng)給定查找相鄰斷裂輪界所用的閾值和延長(zhǎng)修復(fù)的閾值。

（2）輪界區(qū)域排序。對(duì)輪界區(qū)域按中心坐標(biāo)X坐標(biāo)值排序。區(qū)域個(gè)數(shù)如果大于0，則執(zhí)行第 3步，否則執(zhí)行第7步。

（3）取一斷裂輪界區(qū)域。取一個(gè)未形成完整輪界線的輪界區(qū)域，即輪界區(qū)域小于二值圖像寬度。如果取到，則執(zhí)行第4步，否則執(zhí)行第7步。

（4）找下一個(gè)與之相鄰的斷裂輪界區(qū)域且中心點(diǎn)坐標(biāo)X坐標(biāo)值差小于閾值t1的斷裂輪界區(qū)域。取一個(gè)滿足中心點(diǎn)坐標(biāo)X坐標(biāo)值差小于閾值t1的相鄰斷裂輪界區(qū)域。如果取到，則執(zhí)行第5步，否則執(zhí)行第6步。

（5）連接修復(fù)。以原先的左邊輪界區(qū)域最右側(cè)點(diǎn)為起始點(diǎn)，以右邊輪界區(qū)域的最左側(cè)點(diǎn)為終點(diǎn)，然后用白色像素點(diǎn)從起始點(diǎn)開始一直連接到終點(diǎn)。執(zhí)行完連接修復(fù)后，執(zhí)行第3步。

（6）延長(zhǎng)修復(fù)。以輪界區(qū)域中心點(diǎn)為起始點(diǎn)，以最右側(cè)或最左側(cè)點(diǎn)為終點(diǎn)，然后用白色像素點(diǎn)延長(zhǎng)這兩點(diǎn)間的連接線段到二值圖像側(cè)邊。執(zhí)行完連接修復(fù)后，執(zhí)行第3步。

（7）結(jié)束。結(jié)束算法的執(zhí)行，如果算法成功進(jìn)行生長(zhǎng)輪界修復(fù)，則可以輸出修復(fù)后的生長(zhǎng)輪界二值圖像。

1.2.4 圖像存儲(chǔ)

一幅大小為M×N像素的數(shù)字圖像f(x,y)，具有M行，N列。把輸入的二值圖像存儲(chǔ)到長(zhǎng)度為M×N的一維數(shù)組中，黑色像素值為0，對(duì)應(yīng)二值圖像的0，白色像素值為255，對(duì)應(yīng)二值圖像的1。數(shù)字圖像坐標(biāo)如圖3所示。二值圖像存儲(chǔ)示例如圖4所示。

圖3 數(shù)字圖像坐標(biāo)

圖4 二值圖像存儲(chǔ)示例

1.2.5 輪界區(qū)域排序

基于連通區(qū)域判斷[19]，可以得到數(shù)字圖像上所有的輪界區(qū)域。為f(sumx, sumy, area, allPoint, leftPoint, rightPoint, midPoint)，所有點(diǎn)X坐標(biāo)值之和 sumx、區(qū)域所有點(diǎn)Y坐標(biāo)值之和sumy、像素點(diǎn)個(gè)數(shù)area、區(qū)域內(nèi)所有像素點(diǎn)坐標(biāo)allPoint、最左側(cè)點(diǎn)坐標(biāo)leftPoint、最右側(cè)點(diǎn)坐標(biāo)rightPoint、中心點(diǎn)坐標(biāo)midPoint。圖5為輪界區(qū)域數(shù)據(jù)結(jié)構(gòu)體。圖6由1像素點(diǎn)構(gòu)成的連通區(qū)域是一個(gè)輪界區(qū)域示例，a、b、c分別對(duì)應(yīng)需被保存坐標(biāo)的最左側(cè)像素點(diǎn)、最右側(cè)像素點(diǎn)、中心點(diǎn)。如果在X坐標(biāo)方向有并列的最左側(cè)或最右側(cè)坐標(biāo)，任意保存一個(gè)像素點(diǎn)坐標(biāo)。

為了正確修復(fù)輪界，需對(duì)邊緣檢測(cè)后的二值圖像的輪界區(qū)域進(jìn)行排序。從數(shù)字圖像坐標(biāo)來看，各輪界區(qū)域的中心點(diǎn)坐標(biāo)的X坐標(biāo)值按升序排序，X坐標(biāo)值小的靠前。fj(x,y)，j為輪界區(qū)域新序號(hào)。

圖5 輪界區(qū)域結(jié)構(gòu)體

圖6 輪界區(qū)域示例

1.2.6 閾值計(jì)算

本文提出的算法共用到2個(gè)閾值，均通過計(jì)算機(jī)程序自動(dòng)計(jì)算，而不是人為設(shè)置默認(rèn)閾值[11]，在輪界區(qū)域“找下一個(gè)與之相鄰且中心點(diǎn)X坐標(biāo)值差小于閾值t1的斷裂輪界區(qū)域”的過程中用到2個(gè)閾值，設(shè)存在j個(gè)白色區(qū)域，有k個(gè)大白色區(qū)域?qū)挾龋ㄗ钣覀?cè)點(diǎn)Y坐標(biāo)值-最左側(cè)Y坐標(biāo)值）≥閾值t2，即t2=w/a，w為圖像寬度，a為可調(diào)參數(shù)，實(shí)驗(yàn)中a取3，即實(shí)驗(yàn)中把寬度小于圖像寬度 1/3的白色區(qū)域當(dāng)作小區(qū)域，該參數(shù)越大，篩選出來的區(qū)域越多，j個(gè)白色區(qū)域中心點(diǎn)Y坐標(biāo)集合∏=[x1,x2, …,xj]，k個(gè)大白色區(qū)域中心點(diǎn)X坐標(biāo)集合Γ=[x1,x2, …,xk]，閾值t1根據(jù)式（1）計(jì)算而來，min為求最小值函數(shù)。

1.2.7 輪界修復(fù)

用本文作者提出的“基于交互模式和圖像處理的針葉樹生長(zhǎng)輪測(cè)算方法”[18]，在生長(zhǎng)輪邊緣檢測(cè)圖去噪的基礎(chǔ)上，對(duì)不完整輪界線區(qū)域進(jìn)行修復(fù)。具體的操作根據(jù)圖2流程圖進(jìn)行修復(fù)。

連接修復(fù)以原先的左邊輪界區(qū)域最右側(cè)點(diǎn)為起始點(diǎn)，以右邊輪界區(qū)域的最左側(cè)點(diǎn)為終點(diǎn)，然后用輪界像素點(diǎn)從起始點(diǎn)開始一直連接到終點(diǎn)；延長(zhǎng)修復(fù)以輪界區(qū)域中心點(diǎn)為起始點(diǎn)，以最右側(cè)或最左側(cè)點(diǎn)為終點(diǎn)，然后用輪界像素點(diǎn)向右或向左延長(zhǎng)這兩點(diǎn)間的連接線段到二值圖像側(cè)邊。

兩點(diǎn)連線上的點(diǎn)坐標(biāo)計(jì)算：

其中，起始點(diǎn)像素點(diǎn)坐標(biāo)為f(x1,y1)，終點(diǎn)像素點(diǎn)坐標(biāo)為f(x2,y2)，x是大于x1且小于x2的整數(shù)。為了滿足連接線的兩側(cè)不被判斷為8鄰域連接，除了設(shè)置f(x,y)像素點(diǎn)為白色外，還需設(shè)置f(x,y–1)，f(x,y-2)，f(x,y+1)像素點(diǎn)為白色，如果點(diǎn)的坐標(biāo)不在坐標(biāo)體系中，則不進(jìn)行設(shè)置。在具體的二值圖像中，X軸最小值xmin=0、最大值xmax=圖像高度-1，Y軸最小值ymin=0、最大值ymax=圖像寬度-1。

2 結(jié)果與分析

為了驗(yàn)證算法的有效性，本文作者用C#編程語言實(shí)現(xiàn)了所述算法的可視化試驗(yàn)。試驗(yàn)選取了端面邊緣檢測(cè)結(jié)果圖都具有斷裂生長(zhǎng)輪界且生長(zhǎng)輪數(shù)分別為21、104、235、14的杉木、翠柏、北美紅杉、冬青，如表1所示。

表1 修復(fù)結(jié)果

圖7是樣品修復(fù)前后局部對(duì)照?qǐng)D，框線處為被修 復(fù)部位，圖上的框線及A和B是后期添加上去，僅作本文說明用。A框線部位為2個(gè)輪界區(qū)域連接修復(fù)示例，B框線部位為輪界區(qū)域延長(zhǎng)修復(fù)示例。圖 7(a)是1號(hào)樣品修復(fù)前后局部對(duì)照，框線處為延長(zhǎng)修復(fù)。圖7(b)是2號(hào)樣品修復(fù)前后局部對(duì)照，A為連接修復(fù)，B為延長(zhǎng)修復(fù)。圖7(c)是3號(hào)樣品修復(fù)前后局部對(duì)照，框線為連接修復(fù)。圖7(d)是4號(hào)樣品修復(fù)前后局部對(duì)照，該樣品有10處需要修復(fù)，均被延長(zhǎng)修復(fù)或連接修復(fù)。

通過修復(fù)前后對(duì)照?qǐng)D可知，本文提出的修復(fù)算法可較好地對(duì)斷裂的輪界線進(jìn)行連接修復(fù)和延長(zhǎng)修復(fù)。從修復(fù)的樹種來看，不僅可修復(fù)針葉樹端面Canny邊緣檢測(cè)結(jié)果，還可修復(fù)闊葉樹端面Canny邊緣檢測(cè)結(jié)果。從修復(fù)效果來看，修復(fù)后的生長(zhǎng)輪界線走向與原有的輪界線走向有較好的一致性。從修復(fù)的原圖來看，針對(duì)不同的針葉樹，輸入的端面掃描圖Canny邊緣檢測(cè)輪界線粗細(xì)不同，如圖 7(c)的輪界線較其他輪界線粗，這是圖7(c)對(duì)應(yīng)的樹種晚材率高所導(dǎo)致的。

圖7 修復(fù)前后對(duì)照?qǐng)D

3 結(jié)語

本文提出的生長(zhǎng)輪界修復(fù)算法，基于在木材輪盤上對(duì)若干不同方向進(jìn)行傳統(tǒng)人工生長(zhǎng)輪計(jì)數(shù)思路，它有別于獲取整個(gè)木材端面生長(zhǎng)輪界的方法，只處理矩形區(qū)域斷裂的生長(zhǎng)輪界線。通過C#編程語言實(shí)現(xiàn)了所述算法的可視化試驗(yàn)并已集成到自行開發(fā)的“針葉樹宏觀生長(zhǎng)輪數(shù)測(cè)算系統(tǒng)”中。試驗(yàn)結(jié)果表明，該算法可以準(zhǔn)確地修復(fù)斷裂的生長(zhǎng)輪界，為后期獲得準(zhǔn)確的生長(zhǎng)輪數(shù)、計(jì)算生長(zhǎng)輪寬度等工作做好準(zhǔn)備。如何使生長(zhǎng)輪界修復(fù)的區(qū)域準(zhǔn)確沿著生長(zhǎng)輪界伸展方向，值得繼續(xù)深入研究。