林啟招,孫永科,何 鑫,秦 磊,邱 堅
(西南林業(yè)大學 材料科學與工程學院,云南 昆明 650224)
樹木年代學是由美國亞利桑那大學樹木年輪研究試驗室的創(chuàng)始人、天文學家Douglass在20世紀上半葉發(fā)展起來的[1]。與之相關(guān)的研究涉及氣候環(huán)境等領(lǐng)域[2]。在研究中,學者們需要測算樹木的生長輪,于是提出了基于圖像處理技術(shù)的GIS工具、R軟件、圖像處理技術(shù)等方法和手段得到生長輪圖像[3-6]。
數(shù)字圖像處理技術(shù)被用于眾多實驗教學中[7],其可以輔助人們快速測算生長輪,但仍有不足。為了有效識別橫切面微觀圖,Wang等[8]提出了自動識別木材微觀圖生長輪的方法,但未給出修復(fù)橫切面圖像上斷裂斷開生長輪界進的詳細步驟。王燕鳳等[9]提出了年輪圖像雙邊濾波增強方法,但未提及如何解決不理想的生長輪界問題。Vaz等[10]提出的生長輪自動識別與測量系統(tǒng),未考慮到生長輪界不完整的情況。Subah等[11]提出了生長輪交互式視覺分析系統(tǒng),在遇到生長輪斷裂或不完整時,需通過用戶給定起點和終點進行生長輪修復(fù),工作量大的情況下,用戶容易疲勞。Fabijańska等[12]在生長輪自動檢測研究中提出了擬合線條,具有一定弧形的輪界很可能被擬合成理想化的線段。寧霄等[13]提出的基于U-Net卷積神經(jīng)網(wǎng)絡(luò)對年輪圖像進行分割,可達到好的分割效果,但它的不足是需要進行大量的圖像訓(xùn)練學習。
對具有明顯生長輪界的木材端面掃描圖進行Canny邊緣檢測,可得到生長輪界,但所得結(jié)果經(jīng)常會出現(xiàn)邊緣斷裂,即生長輪界不連續(xù)的情況。Canny檢測出現(xiàn)邊緣斷裂是由邊緣附近的灰度值接近引起的[14],被學者們改進后加以運用[15]。木材端面掃描圖較自然界宏觀圖像而言,修復(fù)斷裂生長輪界的要求更高,現(xiàn)有的圖像修復(fù)方法無法用來修復(fù)生長輪界圖像,例如:舒彬等[16]所得的圖像修復(fù)算法可以修復(fù)自然圖像的小區(qū)域破損,劉昱等[17]提出的圖像修復(fù)算法適用于補全大面積破損的圖像。
鑒于此,基于在木材圓盤上對若干不同方向進行傳統(tǒng)人工生長輪計數(shù)思路,本文提出一種對生長輪掃描圖Canny邊緣檢測結(jié)果進行修復(fù)的算法,對作者“基于交互模式和圖像處理的針葉材生長輪測算方法”[18]的前期工作進行深入研究。在前期工作中,發(fā)現(xiàn)輪界不連續(xù)的情況包括中間斷裂和端頭斷開,如圖1所示,圖 1(a)和圖 1(b)是針葉樹翠柏斷面檢測前后對照圖,圖 1(c)和圖 1(d)是闊葉樹冬青斷面檢測前后對照圖。從上往下看,圖 1(b)的第一條輪界線是連續(xù)的、第二條端頭斷開、第三和第四條中間斷裂;圖 1(d)有多處斷裂的生長輪界線。
圖1 木材端面邊緣檢測前后對比圖
本研究試驗對象包括杉木、翠柏、北美紅杉、冬青從髓心向樹皮方向的掃描圖邊緣檢測結(jié)果,其中前3個樣品為針葉樹、第4個樣品為闊葉樹,它們的邊緣檢測圖都存在生長輪界斷裂的情況。
1.2.1 準備工作
在對生長輪界進行修復(fù)前,需導(dǎo)入生長輪掃描圖的Canny邊緣檢測結(jié)果。應(yīng)對所獲得的邊緣檢測圖進行降噪處理,例如:給定閾值,把寬度小于閾值的白色輪界區(qū)域轉(zhuǎn)為黑色,即可達到降噪的目的,作者已在前期工作[18]中進行研究并提出了解決方法。
1.2.2 算法流程
已經(jīng)去噪的生長輪掃描圖Canny邊緣檢測結(jié)果矩形圖,是本文修復(fù)算法的輸入。算法流程圖如圖 1所示。流程涉及的關(guān)鍵過程和技術(shù)在后續(xù)章節(jié)進行陳述。
圖2 生長輪界修復(fù)流程圖
1.2.3 算法描述
本文設(shè)計的算法涉及的主要過程包括:①輸入圖像;②輪界區(qū)域排序;③取一斷裂輪界區(qū)域;④找下一個與之相鄰且中心點X坐標值差小于閾值t1的斷裂輪界區(qū)域;⑤連接修復(fù);⑥延長修復(fù)。具體處理過程描述如下:
(1)輸入圖像。算法導(dǎo)入已經(jīng)去噪的生長輪掃描圖的邊緣檢測二值圖像矩形圖。該過程應(yīng)給定查找相鄰斷裂輪界所用的閾值和延長修復(fù)的閾值。
(2)輪界區(qū)域排序。對輪界區(qū)域按中心坐標X坐標值排序。區(qū)域個數(shù)如果大于0,則執(zhí)行第 3步,否則執(zhí)行第7步。
(3)取一斷裂輪界區(qū)域。取一個未形成完整輪界線的輪界區(qū)域,即輪界區(qū)域小于二值圖像寬度。如果取到,則執(zhí)行第4步,否則執(zhí)行第7步。
(4)找下一個與之相鄰的斷裂輪界區(qū)域且中心點坐標X坐標值差小于閾值t1的斷裂輪界區(qū)域。取一個滿足中心點坐標X坐標值差小于閾值t1的相鄰斷裂輪界區(qū)域。如果取到,則執(zhí)行第5步,否則執(zhí)行第6步。
(5)連接修復(fù)。以原先的左邊輪界區(qū)域最右側(cè)點為起始點,以右邊輪界區(qū)域的最左側(cè)點為終點,然后用白色像素點從起始點開始一直連接到終點。執(zhí)行完連接修復(fù)后,執(zhí)行第3步。
(6)延長修復(fù)。以輪界區(qū)域中心點為起始點,以最右側(cè)或最左側(cè)點為終點,然后用白色像素點延長這兩點間的連接線段到二值圖像側(cè)邊。執(zhí)行完連接修復(fù)后,執(zhí)行第3步。
(7)結(jié)束。結(jié)束算法的執(zhí)行,如果算法成功進行生長輪界修復(fù),則可以輸出修復(fù)后的生長輪界二值圖像。
1.2.4 圖像存儲
一幅大小為M×N像素的數(shù)字圖像f(x,y),具有M行,N列。把輸入的二值圖像存儲到長度為M×N的一維數(shù)組中,黑色像素值為0,對應(yīng)二值圖像的0,白色像素值為255,對應(yīng)二值圖像的1。數(shù)字圖像坐標如圖3所示。二值圖像存儲示例如圖4所示。
圖3 數(shù)字圖像坐標
圖4 二值圖像存儲示例
1.2.5 輪界區(qū)域排序
基于連通區(qū)域判斷[19],可以得到數(shù)字圖像上所有的輪界區(qū)域。為f(sumx, sumy, area, allPoint, leftPoint, rightPoint, midPoint),所有點X坐標值之和 sumx、區(qū)域所有點Y坐標值之和sumy、像素點個數(shù)area、區(qū)域內(nèi)所有像素點坐標allPoint、最左側(cè)點坐標leftPoint、最右側(cè)點坐標rightPoint、中心點坐標midPoint。圖5為輪界區(qū)域數(shù)據(jù)結(jié)構(gòu)體。圖6由1像素點構(gòu)成的連通區(qū)域是一個輪界區(qū)域示例,a、b、c分別對應(yīng)需被保存坐標的最左側(cè)像素點、最右側(cè)像素點、中心點。如果在X坐標方向有并列的最左側(cè)或最右側(cè)坐標,任意保存一個像素點坐標。
為了正確修復(fù)輪界,需對邊緣檢測后的二值圖像的輪界區(qū)域進行排序。從數(shù)字圖像坐標來看,各輪界區(qū)域的中心點坐標的X坐標值按升序排序,X坐標值小的靠前。fj(x,y),j為輪界區(qū)域新序號。
圖5 輪界區(qū)域結(jié)構(gòu)體
圖6 輪界區(qū)域示例
1.2.6 閾值計算
本文提出的算法共用到2個閾值,均通過計算機程序自動計算,而不是人為設(shè)置默認閾值[11],在輪界區(qū)域“找下一個與之相鄰且中心點X坐標值差小于閾值t1的斷裂輪界區(qū)域”的過程中用到2個閾值,設(shè)存在j個白色區(qū)域,有k個大白色區(qū)域?qū)挾龋ㄗ钣覀?cè)點Y坐標值-最左側(cè)Y坐標值)≥閾值t2,即t2=w/a,w為圖像寬度,a為可調(diào)參數(shù),實驗中a取3,即實驗中把寬度小于圖像寬度 1/3的白色區(qū)域當作小區(qū)域,該參數(shù)越大,篩選出來的區(qū)域越多,j個白色區(qū)域中心點Y坐標集合∏=[x1,x2, …,xj],k個大白色區(qū)域中心點X坐標集合Γ=[x1,x2, …,xk],閾值t1根據(jù)式(1)計算而來,min為求最小值函數(shù)。
1.2.7 輪界修復(fù)
用本文作者提出的“基于交互模式和圖像處理的針葉樹生長輪測算方法”[18],在生長輪邊緣檢測圖去噪的基礎(chǔ)上,對不完整輪界線區(qū)域進行修復(fù)。具體的操作根據(jù)圖2流程圖進行修復(fù)。
連接修復(fù)以原先的左邊輪界區(qū)域最右側(cè)點為起始點,以右邊輪界區(qū)域的最左側(cè)點為終點,然后用輪界像素點從起始點開始一直連接到終點;延長修復(fù)以輪界區(qū)域中心點為起始點,以最右側(cè)或最左側(cè)點為終點,然后用輪界像素點向右或向左延長這兩點間的連接線段到二值圖像側(cè)邊。
兩點連線上的點坐標計算:
其中,起始點像素點坐標為f(x1,y1),終點像素點坐標為f(x2,y2),x是大于x1且小于x2的整數(shù)。為了滿足連接線的兩側(cè)不被判斷為8鄰域連接,除了設(shè)置f(x,y)像素點為白色外,還需設(shè)置f(x,y–1),f(x,y-2),f(x,y+1)像素點為白色,如果點的坐標不在坐標體系中,則不進行設(shè)置。在具體的二值圖像中,X軸最小值xmin=0、最大值xmax=圖像高度-1,Y軸最小值ymin=0、最大值ymax=圖像寬度-1。
為了驗證算法的有效性,本文作者用C#編程語言實現(xiàn)了所述算法的可視化試驗。試驗選取了端面邊緣檢測結(jié)果圖都具有斷裂生長輪界且生長輪數(shù)分別為21、104、235、14的杉木、翠柏、北美紅杉、冬青,如表1所示。
表1 修復(fù)結(jié)果
圖7是樣品修復(fù)前后局部對照圖,框線處為被修 復(fù)部位,圖上的框線及A和B是后期添加上去,僅作本文說明用。A框線部位為2個輪界區(qū)域連接修復(fù)示例,B框線部位為輪界區(qū)域延長修復(fù)示例。圖 7(a)是1號樣品修復(fù)前后局部對照,框線處為延長修復(fù)。圖7(b)是2號樣品修復(fù)前后局部對照,A為連接修復(fù),B為延長修復(fù)。圖7(c)是3號樣品修復(fù)前后局部對照,框線為連接修復(fù)。圖7(d)是4號樣品修復(fù)前后局部對照,該樣品有10處需要修復(fù),均被延長修復(fù)或連接修復(fù)。
通過修復(fù)前后對照圖可知,本文提出的修復(fù)算法可較好地對斷裂的輪界線進行連接修復(fù)和延長修復(fù)。從修復(fù)的樹種來看,不僅可修復(fù)針葉樹端面Canny邊緣檢測結(jié)果,還可修復(fù)闊葉樹端面Canny邊緣檢測結(jié)果。從修復(fù)效果來看,修復(fù)后的生長輪界線走向與原有的輪界線走向有較好的一致性。從修復(fù)的原圖來看,針對不同的針葉樹,輸入的端面掃描圖Canny邊緣檢測輪界線粗細不同,如圖 7(c)的輪界線較其他輪界線粗,這是圖7(c)對應(yīng)的樹種晚材率高所導(dǎo)致的。
圖7 修復(fù)前后對照圖
本文提出的生長輪界修復(fù)算法,基于在木材輪盤上對若干不同方向進行傳統(tǒng)人工生長輪計數(shù)思路,它有別于獲取整個木材端面生長輪界的方法,只處理矩形區(qū)域斷裂的生長輪界線。通過C#編程語言實現(xiàn)了所述算法的可視化試驗并已集成到自行開發(fā)的“針葉樹宏觀生長輪數(shù)測算系統(tǒng)”中。試驗結(jié)果表明,該算法可以準確地修復(fù)斷裂的生長輪界,為后期獲得準確的生長輪數(shù)、計算生長輪寬度等工作做好準備。如何使生長輪界修復(fù)的區(qū)域準確沿著生長輪界伸展方向,值得繼續(xù)深入研究。