基于Morse-Smale拓撲特征的文物碎片拼接算法

袁潔 周明全 耿國華 張雨禾

三維破碎物體的匹配與拼接是計算機圖形學、模式識別、可視化技術等眾多領域里一個頗具挑戰的問題.20世紀末,隨著計算機技術的快速發展,將計算機技術引入到文物碎片的虛擬拼接過程中,提高了文物復原的效率,對文物保護和復原至今有著重要的意義.

文物數字化虛擬拼接技術,根據復原過程中是否有專家參與可分為自動復原方法和交互式復原方法.其中,現有的剛體匹配算法主要結合破損文物斷裂區域的幾何特征,作為判斷相鄰碎片之間相似性的主要依據,從而實現破損文物的完整拼合[1].

根據文物的厚度信息,自動化虛擬復原技術可大致分為兩類:薄壁類和非薄壁類文物碎片.早期薄壁類破損文物的虛擬拼接,主要以提取斷裂面輪廓線,并將其映射成為二維平面曲線,實現碎片拼合[1?6].例如Cooper等[1]通過結合碎片輪廓線及其法向,以最大似然函數為基礎進行自底向上的搜索,完成碎片模型的自動拼合.樊少榮等[2]區分三角網格曲面模型的外表面與斷裂面從而正確提取出內輪廓線和外輪廓線,實現碎片的精確拼合.Oxholm等[3]將斷裂面輪廓線上頂點的曲率、撓率以及顏色信息進行結合,組成特征屬性串,采用最長公共子串的方法實現輪廓線間的匹配.Willi等[4]將貝葉斯分析應用在半自動拼合方法中,該算法對斷裂面光滑平整且軸對稱的三維模型有很好的匹配效果,其局限性是只能作用于有限形狀的三維模型.Zhang等[5]基于模板匹配的思想確定顱骨碎片的位置關系,再結合相鄰碎片輪廓曲線實現碎片拼接.Huang等[6]通過計算模型表面頂點的積分不變量值提取碎塊表面尖銳的邊緣線,基于向前搜索技術和表面一致性的約束方法實現碎片的兩兩拼合.空間輪廓曲線匹配方法在幾何特征的數值化計算時具有簡單高效的特點;但由于其采樣點數量有限,因此這類方法對噪聲較敏感并且易忽略厚度信息的部分文物碎片模型效果較差.

非薄壁類剛體匹配對文物重建、數字化遺產保護有著重要的意義[7].Papaioannou等[8]采用Z緩沖方法,獲取三維模型斷裂面投影,計算當前位置斷裂區域的“位置誤差”,通過最小化誤差獲得最優匹配,實現虛擬復原,該方法適用于雕塑、紀念碑等大體積破碎模型的修復.李姬俊男等[9]基于相鄰斷裂面間的凹凸互補性,通過提取約束性的特征簇完成碎片間的兩兩拼合.Sahner等[10]通過計算斷裂區域頂點的坐標和法向,采用層次聚類方法進行破損文物碎塊之間的拼合.由于斷裂面中包含了較多的特征信息,能夠準確地反映模型斷裂區域的鄰接關系,因此,針對非薄壁類文物碎片,基于空間曲面的碎片拼合,具有一定優勢,但針對斷裂面部位受損較嚴重的碎塊,匹配的正確性無法保證.

因此針對破損文物斷裂部位邊緣受損而引起的輪廓線不能充分表示斷裂面幾何信息的問題,本文基于文獻[11]中的特征線提取方法,提出一種基于Morse-Smale的斷裂面拓撲幾何特征的破碎剛體自動拼接方法.本文算法包括斷裂面拓撲圖的生成、四邊形描述符的定義、匹配集的確定以及碎片拼合4個主要部分,步驟如圖1所示.本文算法的核心思想是采用曲面四邊形描述符表示斷裂面幾何特征信息.在傳統的基于空間曲線的碎片拼合方法中,針對斷裂部位邊緣受損而提取的輪廓線不精確,錯誤的曲線匹配導致斷裂面滲透;傳統的基于空間曲面的碎片拼合方法,并沒有考慮斷裂面特征點之間的拓撲關系,可能出現斷裂面的錯誤配準.因此,為了提高碎片拼合的效果,基于Morse-Smale復形[12]本身的四邊形性質及其可控性和魯棒性,本文借助能有效融合斷裂面上凹凸信息的四邊形曲面,更好地反映空間曲面特征,并且可精確地確定鄰接碎片的匹配關系,同時采用凹凸互補閾值方法篩選出最優匹配集,將計算量較大的空間點對應關系轉換成少數的四邊形匹配,有效地提高了剛體轉換的效率,最后采用窮舉搜索的方法實現破損碎片精確拼接.

圖1 本文方法步驟Fig.1 Procedures for the method proposed in the paper

1 拓撲圖提取及相關術語

本文通過曲面四邊形表示斷裂面上的幾何特征,從而實現碎片拼合,因此,本節提取能有效融合斷裂面幾何特征的拓撲圖.以三角網格模型為研究對象,首先,依據頂點的顯著度指標函數,提取關鍵點并分類;然后,采用最大角度法構建特征線;最后,對提取的斷裂面Morse-Smale復形進行簡化.

1.1 Morse-Smale理論

Morse理論主要用來研究n維空間內的d維流行M(M?Rn)的形狀,通過Morse函數,將二維流行表面分割為一系列Morse單元,而Morse-Smale復形為這些莫斯單元的集合.其主要思想是將形狀的拓撲研究同映射函數集合特性的定量分析結合起來,通過定義特征點并判斷其類型,按照一定的尋徑方法利用特征線將其連接起來,形成關鍵網絡或轉換成其他數據模型,從而表達流形表面的拓撲結構與形態特征,是對特征提取和幾何形狀的分析,保存其拓撲性質及形態的強大工具.

拓撲特征是對圖像中區域結構形狀的總體描述,其特點是不受旋轉與平移等變形影響,描述能力強、計算速度快的優點,可有效、準確地表示物體模型或圖像的局部形狀.常見的圖像拓撲特征有基于骨架結構表示的目標檢測算法和基于形態圖表示的三維識別算法.

基于Morse函數的特征提取算法在由二維圖像擴展到三模型表面時,依據三角網格頂點的顯著度度量指標(曲率值、法向量改變值等頂點屬性)提取的極值點、升降線與Morse-Smale復形做為模型表面的顯著特征,避免了沒有實際意義的非顯著特征提取,保證了Morse-Smale復形的拓撲完整性.

1.2 特征點提取

針對三維模型表面的網格數據,特征提取是指采用指標函數來提取出具有明顯特征的點集合.本文采用曲度函數f(v)作為三角網格上各頂點的顯著度指標函數[13],定義如下:

其中,v為三角網格上任意一個頂點;k2max(v)和k2min(v)分別為頂點v的最大和最小主曲率.

采用曲度函數作為指標函數是為了將曲率絕對值最大的特征點提取出來.斷裂面上曲率絕對值最大的點即為凹點、凸點和鞍點,均是能夠表示物體模型表面幾何特征的點.

1.3 特征線構建

基于上節提取的斷裂面特征點,本節通過最大角度法[14]構建特征線,構造能完整表示斷裂面幾何關系的拓撲圖.

首先,對提取的特征點進行判定,并對其進行分類.定義vi(i=1,2,···,n)為兵馬俑碎塊斷裂面的三角網格上的任意頂點,其中,n為三角網格數,其所對應的顯著度指標函數為f(vi).若其一階鄰域點的顯著度函數值均小于該頂點的顯著度函數值,則定義該頂點為極大值點;同理,若其一階鄰域點的顯著度函數值均大于該頂點的顯著度函數值,該頂點則被定義為極小值點.若目標點的鄰域中有某段連續區域內的頂點,其指標函數值均大于目標點的函數值,則稱該段連續區域能取得的最大范圍為上升區間;若其指標函數值均小于目標點的函數值,則稱該段連續區域能取得的最大范圍為下降區間.定義鞍點是同時連接兩個上升區間和兩個下降區間的特征點.如圖2所示,分別對應斷裂面上曲率絕對值最大的頂點,即特征點.

圖2 頂點鄰域關系圖Fig.2 Diagrams of vertex neighborhoods

基于以上分析,本文采用最大角度法構建特征線,構建規則如下:對任意一個鞍點都存在2條沿著指標函數值下降至極小值點和2條沿著指標函數值上升到極大值點的特征線,計算鞍點與其鄰域頂點的邊,選擇上升或下降最快的點作為下一個路徑點位,構成特征線.以升弧為例,以鞍點為出發點,方向為鞍點上升區間中指標函數值最大的鄰點所指示的方向,直至到達極大值點.設vi為三角網格上的任意頂點,其方向為則有:

式中,vj為vi的一階鄰域點.降弧構造方式與升弧類似,沿著最大角度反方向搜索.

1.4 拓撲簡化

針對所提取的斷裂面拓撲圖,本節采用注銷和移動操作[15]進行簡化,從而獲取準確的純四邊形拓撲圖.

Morse-Smale復形中的特征點以及升弧、降弧構成初始特征集合.一般情況,特征線會將模型表面劃分成若干四邊形區域,所有四邊形都是由兩個鞍點,一個極大值點和一個極小值點構成.在現實情況中,存在一些特殊問題,如圖3(a)深色弧線部分生成的環狀,對某個鞍點在上升空間中尋找極大值點時,兩條升弧同時連接到一個極大值點,使得拓撲圖中出現環狀回路,以及中心區域可能出現孤立極值點,如圖3(c).

圖3 環狀和孤立臨界點Fig.3 Diagrams of cyclic and solitary critical points

綜上所述,針對拓撲圖中環狀回路的問題,本文采用注銷的方法消除多余或者不規則的關鍵點.注銷是一種雙重邊緣收縮的方法,即去除鞍點與極值點之間連成的兩條重復特征線,其中包含所有的弧線以及這個鞍點和兩個胞腔,結果如圖3(b).若鞍點v所連接的兩條降弧和鞍點x的極小值點重復,而特征線均為升弧,因此將這些升弧合并為一條升弧,且刪除重復的兩條降線,簡化結果如圖3(d).

針對初始特征集合,某些鞍點在特征線構建過程中,搜索下一個路徑點位時可能存在微小誤差,線段間誤差積累,導致構建的初始莫爾斯復形中存在一些次要特征線,如圖4所示,鞍點的下一個最優路徑點位可能是沿著圖4(c)中虛線所示的方向,而實際在構建過程中是沿著圖4(a)或4(b)中三角形邊方向.

圖4 路徑點位選擇Fig.4 Paths selected by points

文獻[15]中將這些次要特征線及其所連接的臨界點一并刪除,但是針對模型表面的幾何特征,刪除誤差積累導致的次要特征線可能會造成模型表面特征的不完整性,因此本文基于此方法進行改進,簡化規則如下:1)通過定義特征線的顯著度,確定其重要性程度;2)設定閾值δ0,當特征線的顯著度低于該閾值時,定義為次要特征線;3)針對次要特征線進行移動操作從而尋找最準確的特征線.

一般地,升弧上的指標函數值與其所在區域的指標函數差值越大,則該升弧越重要,顯著度越高;降弧上的函數值與其所在區域的函數值差值越小,則該降弧越重要,顯著度越高[16].因此,本文定義特征線的顯著度為該特征線與其鄰接區域的平均指標函數值之差.

基于以上分析,每個鞍點連接2條升弧和2條降弧,本文將其合并為一條升弧R和一條降弧r,假設降弧r,其所對應的上升域為D(m1)和D(m2);升弧R,對應的下降域為A(m1)和A(m2).則升弧R相對于A(mi)的顯著度為

其中,i=1,2.

最大角度法構建特征線過程中,特征線方向沿著三角網格上三角形邊前進.因此,每條特征線都是由若干個特征點pi(i=1,2,···,n)連接而成的,即升弧R是由若干首尾相連的線段組成,如圖5所示.

圖5 特征線構成Fig.5 Feature line composition

因此,R上的積分可近似表示為

式中,pj1,pj2是特征線R上的第j條線段上的2個端點,區域A(Mi)上的積分則可近似表示為

式中,pk1,pk2以及pk3分別為區域A(Mi)中第k個三角形的頂點;Ak為該三角形的面積.

通常,對任意的特征線都有兩個鄰接區域,因此會相應的有兩個顯著度.假設特征線R的2個鄰域分別為A(M1)和A(M2),相對應的2個顯著度為δ(R,A(M1))和δ(R,A(M2)),取較小值為該特征線的整體顯著度:

同理,降弧r的兩個鄰接區域為D(m1)和D(m2),則r的顯著度為

其中,i=1,2.則該條降弧的整體顯著度為

對任意特征線,若其特征線顯著度低于所給定的閾值,則將該特征線進行移動操作,即從該特征線的起始鞍點開始,重新搜索下一個最優的路徑點位,從而獲取該區域內最準確并且顯著度最高的特征線,構造斷裂面上的純四邊形拓撲圖.

2 特征描述符定義及搜索

為了將拓撲圖中曲面四邊形構造成能夠有效表示斷裂面幾何信息的四邊形描述符,首先,定義基準點和0值面,采用頂點高度差的方法對曲面四邊形進行編碼,獲得四邊形描述符;然后,通過凹凸互補誤差篩選出最優四邊形匹配集.

2.1 四邊形定義

斷裂面上的Morse-Smale復形,由多個四邊形曲面構成,這些曲面四邊形將斷裂面進行剖分,四邊形本身具有各向異性的特點,使得四邊形網格更易與模型特征對應,因此本文將斷裂面上的曲面四邊形作為特征描述符,以相鄰碎片斷裂面間的凹凸特性為基準進行碎塊拼合.設待匹配的兩個斷裂面分別為?A和?B,定義?A為源面,?B為目標面,對于源面?A上的Morse-Smale復形中的任意曲面四邊形,假設兩個鞍點分別為s1和s2,極大值為M1,極小值為m1.

如圖6所示,設以極大值點M1所在的切平面l0為0值面,定義極大值點M1為基準點,鞍點s1、s2,極小值點m1為目標點.用h1表示鞍點s1到該切平面的垂直距離,極小值點m1到該切平面的垂直距離為h2,鞍點s2的對應距離為h3,極大值點的對應距離為h0,將該高度差值作為每個頂點的屬性值,并按照一定規則進行編碼.

圖6 頂點高度差值Fig.6 Vertex height difference

編碼規則:本文以源面?A上任意曲面四邊形為例,以四邊形中極大值點開始,順時針編碼,如圖6所示的四邊形中,極大值點M1為基準點,則定義hA0=0;依次鞍點s1、極小值點m1、鞍點s2所對應的屬性值為hA1,hA2,hA3,則該四邊形的一個屬性編碼為.

類似地,目標面?B在選取基準點時,根據匹配碎塊斷裂面的凹凸互補性,基準點選擇特征四邊形中的極小值點m1,以其切平面作為0值面,然后求得其余三個目標點到0直面的垂直距離,以同一編碼規則獲得任意特征四邊形的屬性編碼.

2.2 最優匹配集搜索

基于本文定義的四邊形屬性,從破損的兵馬俑斷裂面上提取出特征四邊形集,如圖7所示,但由于特征四邊形的數目較多,斷裂面?A和?B之間的匹配關系會比較復雜.因此,本文設置凹凸互補閾值條件,減少誤匹配四邊形集合個數,求取少量的最優匹配集合.

圖7 斷裂面拓撲圖Fig.7 Topological diagram of fracture surface

假設?A和?B為兩個待匹配的斷裂面,和分別為上的初始特征四邊形集合.設sAi和sBi是待匹配的兩個斷裂面?A和?B上相對應的可匹配特征四邊形,則定義凹凸互補誤差為

假設四邊形sA1為斷裂面?A上的任意特征曲面四邊形,屬性編碼如下:

相對應的,斷裂面?B上對應的特征四邊形sB1為

則凹凸互補誤差公式為

其中,x1=hA0?hB0,x2=hA1?hB2,x3=hA3?hB3,x4=hA2?hB1.當ε越接近0,表示該對特征四邊形的匹配度越高.因此,設定凹凸互補閾值ε0＞0,對四邊形sA1,若在斷裂面?B上滿足ε(sA1,sBi)＜ε0的特征四邊形存在m個,則取最優匹配對 (sA1,sB1)=, 并將該最優匹配對加入到最優匹配集s=;否則,刪除該匹配對.本文通過凹凸互補閾值方法約束特征四邊形,篩選出匹配度較高的匹配集,提高斷裂面的匹配精度并且降低剛體轉換的復雜性.圖8是初始特征匹配到最優特征集匹配的結果對比圖.

圖8 匹配集優化結果Fig.8 Optimization results of matching set

3 剛體變換

對特征四邊形進行初步篩選后,確定了待匹配斷裂面上的最優特征集合S′A和,以及由此構成的最優匹配集S=, 但由于兩個碎片之間的方向與位置均不同,為了能將碎塊完整的拼合,需根據兩個四邊形匹配集計算出對齊碎片的旋轉矩陣R和平移矩陣T[17].

基于特征四邊形區域的質心坐標,采用四元組方法[18]求解出最優匹配集中的3組質心不共線特征四邊形匹配對的旋轉矩陣R和平移矩陣T,并采用窮舉搜索法[19]篩選出匹配正確率較高的剛體變換,實現碎片的精確拼合,算法步驟如下:

步驟1.設任意特征四邊形S的4個頂點為{p1,p2,p3,p4},點pi的三維坐標為(xi,yi,zi),定義S的質心為.

步驟2.根據三角形理論,從匹配集中篩選出所有滿足該條件且質心不共線的3組特征四邊形對,構成集合,其中,,設M中的元素個數為e個.

步驟3.從上述集合M中選取任意一組特征四邊形對,采用四元組法計算出斷裂面?A和?B之間的剛體變換矩陣R和T,遍歷集合M中所有的特征對元素,求得剛體變換數組D.

步驟4.計算D[t],t=1,···,e的匹配正確率;c為正確匹配對個數,初值設為0;集合S中任意的特征四邊形對經過D[t]變換之后,若滿足,則表明該特征四邊形對匹配正確,c=c+1,否則為錯誤匹配;遍歷集合S可求出D[t]的匹配正確率q=c/m;遍歷整個數組D[e],計算所有特征對剛體變換的正確匹配率.

步驟5.選取匹配正確率最高的剛體變換,完成待拼接碎片的精確對齊.

4 實驗結果與分析

為檢驗本文算法的有效性,以秦始皇陵1號坑出土的部分兵馬俑碎片網格模型作為實驗數據.本文實驗采用Visual C++(Visual studio 2010)和 OpenGL編程,在 Intel Core i7-3770 CPU/3.4GHz,4GB內存的PC機上實現.實驗閾值均是通過經驗值以及數據統計獲取的,δ0=0.02,ε0=5.0mm,σ1=2.0mm,碎片在斷裂處均存在邊界受損,幾何特征缺失現象.

4.1 本文算法結果

圖9～11所示為G10-26號俑、G10-19號俑和G10-23號俑采用本文算法在部分鄰接碎片上的實驗結果,可以看出,本文算法在邊緣受損的碎片模型上可取得較好的拼合結果.圖12～14所示為本文算法在一組整俑上的拼接結果,實驗結果表明,本文算法在斷裂部位邊緣受損的整俑模型上,如圖12中矩形框標注的有小部分缺失的碎片,仍可取得較好的拼合結果.本文實驗數據采用的兵馬俑模型均為陶俑模型,碎片斷裂面具有一定厚度,碎片數量較多并且網格數據總量較大,因此在鄰接碎片的搜索時復雜度較高.實驗數據如表1所示,鄰接碎片在斷裂面提取的特征四邊形隨碎片網格數增加,碎片的數量決定了鄰接碎片的斷裂面個數,其對特征四邊形個數也有直接的影響;由表1中數據可得,本文算法可取得較高的匹配精度.

圖9 G10-26號俑部分鄰接碎片拼接結果Fig.9 Reassembly results of some adjacent fragments of Warriors G10-26

圖10 G10-19號俑部分鄰接碎片拼接結果Fig.10 Reassembly results of some adjacent fragments of Warriors G10-19

圖11 G10-23號俑部分鄰接碎片拼接結果Fig.11 Reassembly results of some adjacent fragments of Warriors G10-23

圖12 G10-18號俑部分鄰接碎片拼接結果Fig.12 Reassembly results of some adjacent fragments of Warriors G10-18

本文將圖7中的碎片裁剪為薄壁模型后,斷裂面上的特征四邊形提取結果示意圖如圖15,有效特征四邊形定義為四邊均處于碎片模型斷裂面上的四邊形,如圖15中矩形框標注所示.由于薄壁模型斷裂面上的特征點數量較少,提取的有效特征四邊形數量不足從而無法完整表示斷裂面的幾何結構,造成特征匹配不準確的問題,從而影響后續鄰接碎片位置的判斷.因此,本文算法針對薄壁類文物碎片模型,拼合效果較差.

表1 兵馬俑碎片實驗數據Table 1 Experimental datas of the Terracotta Army fragments

圖13 G10-36號俑部分鄰接碎片拼接結果Fig.13 Reassembly results of some adjacent fragments of Warriors G10-36

4.2 算法運行時間

本文算法針對成組的兵馬俑碎片拼合時間如表2所示,其中,斷裂面曲面四邊形特征提取為T1,本文每次提取2個碎片的斷裂面四邊形特征;特征匹配集獲取為T2,同理,每次獲取2個碎片的特征匹配集;T3為碎片拼合.實驗數據表明,隨著碎片數據的增加,運行時間也會隨之增加,這是因為斷裂面上四邊形特征描述符提取為本文算法最耗時的部分,需要遍歷斷裂面網格上的所有頂點.

圖14 G10-22號俑部分鄰接碎片拼接結果Fig.14 Reassembly results of some adjacent fragments of Warriors G10-22

圖15 薄壁碎片斷裂面特征提取圖Fig.15 Feature extraction in thin fragments

4.3 對比實驗結果

圖16所示分別為采用基于形態圖方法[20]及基于斷裂面幾何特征方法[7]的碎片拼合結果,選取3組碎片模型為實驗數據,分別為:G10-18號俑,G10-19號俑及G10-36號俑部分鄰接碎片.該3組鄰接碎片具有不同的特點:G10-18號俑碎片模型斷裂部位邊緣受損較少;G10-19號俑碎片模型斷裂部位邊緣受損較嚴重且斷裂面存在缺損;G10-36號俑鄰接碎片斷裂部位面積較大,碎片存在缺失.

表2 本文算法運行時間Table 2 Execute times of the proposed algorithm

若兩個斷裂面可實現正確匹配,那么當兩個碎片完成精確配準之后,相鄰的兩個斷裂面之間不應該發生明顯的碰撞,存在縫隙過大以及相互重疊交錯[21].圖16(a)所示為G10-18號俑碎片采用文獻[20]中提出的基于形態圖方法的拼合結果,將碎片模型轉換成點云數據,可看出,當碎片的斷裂處存在缺損時,基于斷裂面出輪廓線的拼合方法會出現相互滲透的現象,如圖16(a)中矩形框中標注所示,難以達到較好的拼合結果;圖16(b)中所示為采用文獻[7]的基于斷裂面特征簇的拼合方法對G10-19號俑部分碎片的拼合結果,由于碎片斷裂面存在缺損,則導致幾何信息無法有效表示斷裂面結構,造成斷裂面相互滲透現象,如圖16(b)中矩形框標注所示;當碎片存在小部分的缺損時,則會造成不準確的位置匹配,出現碎片表面紋飾信息不連續的現象,如圖16(c)所示.

實驗結果表明,相較于其他拓撲特征的拼合方法和基于斷裂面幾何特征的拼合方法,本文算法通過融合斷裂面幾何特征的四邊形曲面,增強了特征描述符的約束性,能夠準確有效地判斷碎片的鄰接位置,針對斷裂面邊緣受損的碎片具有更好的拼合效果.

5 結語

針對破損文物斷裂部位邊緣受損而引起的輪廓線不能充分表示斷裂面幾何特征的問題,本文提出了一種基于Morse-Smale的斷裂面拓撲特征的破損文物自動拼接算法.該算法將斷裂面上特征點的曲度函數作為特征檢測算法,根據斷裂面幾何特征提取特征點;采用最大角度法對特征點進行特征線構建,并計算特征線重要度簡化拓撲圖,得到斷裂面的純四邊形拓撲圖;將拓撲圖中四邊形曲面構造成為能有效表示斷裂面集合特征的描述符,彌補了傳統方法反映局部特征的不足,提高了匹配的精度和效率,同時擴展了匹配算法的使用范圍.本文算法需要較為精確地識別出碎片的斷裂面,解決薄壁類破損碎片斷裂面特征平滑,無法取得精確地拼合結果的難題.

由于本文算法采用了拓撲圖四邊形為特征描述符,因此在提取特征描述符時較耗時,尋找一種更加魯棒、高效地特征提取算法以減少計算量,將是我們下一步工作的重點.

圖16 傳統算法實驗結果Fig.16 Experimental results of traditional methods