何智翔 丁曉青
(智能技術與系統(tǒng)國家重點實驗室,清華信息科學與技術國家實驗室,清華大學 電子工程系,北京100084)
在醫(yī)學和高空目標觀測等領域,存在大量的不同光質(zhì)多模圖像之間的目標檢測問題。因為在這些領域中,存在大量不同感光成像設備生成的圖像,我們稱這些圖像為異質(zhì)多模(multimodal或者 multisensor)圖像[1,2]。
在有關異質(zhì)多模圖像的研究中,Ma Jun和Zeng等[3,4]利用 sobel算子生成原圖像的邊緣圖,然后利用邊緣圖的互相關進行匹配。參考圖像和觀測圖像雖然光質(zhì)不同,但是處于同一視角并且是同一場景,背景比較簡單。類似的方法也出現(xiàn)在了文獻[5,6]中,也是邊緣圖被用來進行配準;不同的是,參考圖像和觀測圖像不是同一視角;但是因為規(guī)定圖像間只有有限的平移、旋轉、尺度變換,所以通過簡單的仿射變換就可以消除。Cheng和Zheng等[7]則不使用邊緣圖,而是利用SVM從提取到的邊緣中抽取物體的輪廓,然后使用hausdorff距離進行輪廓匹配。Manjunath和Huang[7,8]也使用輪廓,但他們使用其他方法來提取輪廓并匹配它們。此外,在文獻[7-9]的研究中,它們的參考圖像和觀測圖像雖然是同一場景,但背景并不嚴格一致,視角也不完全相同,但是相差不多。此外,在最近的研究中,如 Michael和Krotosky等[2,10],嘗試利用三維場景數(shù)據(jù)對異質(zhì)多模圖像進行配準,不是簡單的圖匹配。這些研究,大多都是針對異質(zhì)多模圖像的配準,通過配準來檢測目標。
本文研究的主要對象是不同視角的異質(zhì)多模圖像,參考圖像和觀測圖像之間除了光質(zhì)不同外,還具有以下2個主要的困難:首先就是參考圖像和觀測圖像的視角不同,仿射變換無法完全消除其影響;其次就是參考圖像和觀測圖像包含復雜的背景,我們采用的實驗數(shù)據(jù)是可見光俯視參考圖像和紅外光側視觀測圖像,如圖1所示。其中三維物體的高度信息導致了場景不一致且紅外光觀測圖的清晰度低,受到了噪聲強烈的干擾。
圖1 本文研究所使用的數(shù)據(jù)Fig.1 Image data used in this research
就我們所知,現(xiàn)有的這些研究[3-9]的實驗數(shù)據(jù)都具有不太復雜的背景(噪聲少,場景比較相似等)以及相近的視角,所以他們的方法可能不能很好地處理這種情況。
不同于前面提到的研究[3-9],本文不使用基于邊緣的特征。因為依賴于邊緣的方法,容易受到圖像中噪聲和視角不同所帶來的影響,特別是對于本文研究的低清晰度紅外觀測圖。而本文采用的區(qū)域特征對于視角變換可能引起的形狀畸變和噪聲的干擾更加魯棒,相對比較穩(wěn)定。
最終,本文在貝葉斯最小風險理論的指導下,利用馬爾可夫隨機場(Markov Random Field或MRF)模型來描述在異質(zhì)多模圖像中用MSER[11]的方法提取得到的地面穩(wěn)定區(qū)域特征之間的空間約束關系及統(tǒng)計特性,通過配準異質(zhì)多模圖像實現(xiàn)最終的目標檢測。
馬爾可夫隨機場是建立在馬爾可夫隨機鏈的基礎上的,在MRF中,當前位置的狀態(tài)只與相鄰位置上的狀態(tài)有關。MRF模型提供了不確定性描述與先驗知識聯(lián)系的紐帶,并利用觀測圖像,根據(jù)統(tǒng)計決策和估計理論中的最優(yōu)準則確定目標函數(shù)。
利用MRF模型,我們能夠建立基于圖像中的地面穩(wěn)定區(qū)域特征的圖模型,從而可以依賴于圖模型的配準進而在觀測圖像中檢測甚至定位參考圖像中的目標。
通過比較圖1中的參考圖像和觀測圖像,我們可以得到關于異質(zhì)多模圖像的一些特點。首先,異質(zhì)多模圖像中的同一物體區(qū)域由于圖像模式不同,會呈現(xiàn)出不同的紋理特征;而由于同一場景中的不同物體存在各自的變化規(guī)律,因此也很難通過類似于光照變化的調(diào)整方法來消除這種差異。其次,雖然物體的紋理不一致,但是它們的大致形狀是相似的,特別是在同一視角之下。
正是由于形狀是異質(zhì)多模圖像間最穩(wěn)定的特征,現(xiàn)有的一些基于紋理的特征并不適用于異質(zhì)多模圖像的目標檢測,而需要使用與物體形狀相關的特征。以HOG特征為例,如圖1中的廣場,因為紋理特征不一致,在局部區(qū)域上計算HOG特征,會發(fā)現(xiàn)兩者的梯度強度和方向都不一致,從而導致檢測失敗。
在文獻[11]中,Matas通過將注水算法引入連通域提取的過程,從而得到了一種有效的提取穩(wěn)定區(qū)域的方法,稱為最大穩(wěn)定極值區(qū)域提取(Maximally Stable Extremal Regions,簡 寫 為MSER)。由于MSER選取的是那些形狀在一定步長的灰度變化下不敏感的區(qū)域,因此它能夠有效地克服噪聲的干擾,提取在不同圖像模式下形狀足夠穩(wěn)定的物體區(qū)域,比如湖泊、廣場、水田等。因此,MSER可以被用來提取異質(zhì)多模圖像中的穩(wěn)定物體區(qū)域來描述物體的形狀。
在最后的目標檢測過程中,需要尋找一組從觀測圖像中提取的區(qū)域配準圖模型。我們采用RANSAC[12]的思想來尋找這個最佳配準。
簡單說來,RANSAC算法就是一種從樣本中準確擬合數(shù)學模型的算法,包括去除噪聲點(野值)和留下有效值等等,采用隨機抽樣驗證的方法。它利用已知模型,從樣本集中隨機選取若干個點,并擬合模型參數(shù),然后根據(jù)參數(shù)將樣本集中的點添加入候選集合,直到誤差小于閾值或者候選集合中的點滿足一定的數(shù)目要求。采用RANSAC算法可以從可能的與模板匹配的物體區(qū)域集合中快速尋找和圖模型的最佳匹配,而不必使用窮舉算法去遍歷每一個可能的組合與圖模型的配準結果。
常用的利用邊緣特征進行匹配的算法是Hausdorff距離[7,13]。在前面介紹的利用邊緣特征對多模圖像進行匹配的研究中,Cheng和Zheng[7]使用該算法來衡量多模圖像中邊緣的相似程度。
Hausdorff距離屬于一種邊緣匹配測度,它以2個邊緣之間的最大的最小距離作為度量,能夠較好地克服噪聲的干擾;但是對模板形狀的要求比較高。正是因為基于邊緣特征的方法對于形狀不一致比較敏感,本文最終選擇了利用 MSER提取穩(wěn)定區(qū)域的方法,通過計算二值連通域的互相關值來進行匹配。盡管利用 MSER在異質(zhì)多模圖像中提取的同一物體的穩(wěn)定區(qū)域也并不一定能保證形狀完全一致,但是在利用給定的攝像機外部參數(shù)進行視角變換的情況下,足夠保證模板與對應區(qū)域的正確匹配。
在第4節(jié),我們將給出本文的利用區(qū)域特征匹配的方法和文獻[7]中采用Hausdorff距離進行邊緣特征匹配的方法的實驗結果對比。
我們的目標檢測方法的核心是利用 MSER提取物體區(qū)域,然后在MRF模型的基礎上進行匹配。前面已經(jīng)提到過我們?yōu)槭裁催x擇區(qū)域特征而不是邊緣特征來描述物體的形狀,在這里我們將給出更加詳細的解釋,主要是以下3點理由:首先是觀測圖像的清晰度很低,噪聲的干擾非常嚴重,物體的形狀很難依靠邊緣特征進行準確的描述;其次是缺乏精確的視角變換參數(shù),不能通過變換使同一物體的形狀完全一致;最后因為使用MSER的方法來提取不同模式圖像中的同一物體,使用區(qū)域特征更加適合。
由于觀測圖像與參考圖像視角的差異極大,尤其是參考圖像為俯視圖而觀測圖像為側視圖,因此,圖像中物體的高度信息會對檢測產(chǎn)生極大的干擾。為了降低物體的高度信息對檢測的影響,我們選擇那些高度信息可以忽略的地面物體作為穩(wěn)定區(qū)域特征,并利用這些地面區(qū)域特征作為對圖像的描述。
于是在貝葉斯最小風險理論下,地面穩(wěn)定區(qū)域看做是圖像上的前景RFG,而其他區(qū)域則是背景RBG。于是圖像I=RFG∪RBG且RFG∩RBG=?。地面穩(wěn)定區(qū)域在圖像中的位置L=(l1,l2,…,ln),那么前景圖像的條件概率分布為
背景圖像的條件分布概率為
最優(yōu)的配準結果為
(3)式直接求解是困難的。
但是在MRF模型的框架下,圖像中的相鄰地面物體區(qū)域之間存在聯(lián)系,于是可以用一個統(tǒng)一的模型描述地面物體區(qū)域的結構約束和地面物體區(qū)域的統(tǒng)計特性。
似然概率為
由(5)式可知,似然概率可以轉化為下式:
由于拋開空間位置約束,各個穩(wěn)定區(qū)域可以認為是相互獨立的,所以(7)式成立
令
于是
其物理意義在于,最優(yōu)的圖模型匹配結果要求在觀測圖像中匹配上的地面穩(wěn)定區(qū)域為參考圖像上的前景的概率越大越好,為背景的概率越低越好,且滿足空間約束的概率要高。
基于該模型,我們的方法的基本步驟如下。
首先進行視角變換,因為多模圖像間最穩(wěn)定的特征是與形狀相關的特征,而參考圖像和觀測圖像之間的視角差異太大,必須將它們變換到相近的視角上,這樣才能利用形狀來進行匹配。
接下來,MSER被用來提取參考圖像和觀測圖像中的物體穩(wěn)定區(qū)域。之后,參考圖像中的在地面上的物體穩(wěn)定區(qū)域被選為模板(前景圖像),并構建MRF模型。
最后,通過計算穩(wěn)定區(qū)域之間的近似程度,建立每個模板在觀測圖像上的候選匹配區(qū)域集合,并按照(9)式,在滿足空間位置約束的情況下選取最優(yōu)的匹配。
在本文的研究中,由于攝像機的視角參數(shù)已知,所以我們直接利用給定的攝像機視角參數(shù)將觀測圖像變換到俯視的視角上。不變換到側視圖中是因為視角變換的誤差在側視圖中很難通過簡單的旋轉和縮放消除,如圖2所示。
圖2 參考圖像變換與觀測圖像變換的比較Fig.2 The comparison between the transformation results of reference image and sensed image
從圖2還可以看到,視角變換確實對存在高度信息的物體形狀產(chǎn)生了很大的影響,因此,我們抽取合適的地面物體區(qū)域作為模板。
視角變換后,我們將利用MSER提取圖像中的穩(wěn)定區(qū)域。由于MSER算法提取穩(wěn)定區(qū)域的特點,它只能夠提取灰度小于某個閾值或者大于某個閾值的區(qū)域。為了能夠充分利用MSER,我們引入了圖像灰度變換。
本文不僅僅在原圖像I中提取了穩(wěn)定區(qū)域,還在反色圖像Irev以及利用公式(10)和(11)進行灰度變換后的圖像Itra中提取了穩(wěn)定區(qū)域。圖3顯示了灰度變換后的參考圖像,對于觀測圖像也是如此。I(x,y)是圖像上位于(x,y)處的灰度值,b為圖像的寬度,h為圖像的高度。
圖3 灰度變換后的圖像Fig.3 The gray transformation of an image
從參考圖像中得到的這些地面穩(wěn)定物體區(qū)域將被看做前景圖像用來建立圖模型以描述目標,而從觀測圖像中提取到的穩(wěn)定區(qū)域則被作為待匹配的物體區(qū)域,如圖4所示。從每個模板中提取特征,加上對應的幾何約束,就構成了參數(shù)化的圖模型描述。
圖4 參考圖像和觀測圖像中提取的MSERFig.4 The stable regions extracted from the images
在匹配模型之前,我們需要計算從觀測圖像中提取的穩(wěn)定區(qū)域與每個模板之間的距離dcorr來衡量它們之間的形狀相似程度。該相似程度反映了公式(9)中概率的大小。在本文中,我們使用最簡單的互相關方法(CC),如(12)式所示。其中,T為模板對應的二值圖,I為觀測圖像中提取的穩(wěn)定區(qū)域對應的二值圖。ST是T中前景像素的面積,SI則是I中前景像素的面積。
在得到第i個模板與觀測圖像中提取到的第j個穩(wěn)定區(qū)域的距離dj,corri以后,如果
那么第j個區(qū)域將被加入與第i個模板對應的第i個候選匹配區(qū)域集合中。
圖5是隨機選取的2個模板的匹配結果示意圖,其中,灰白色區(qū)域表示的是模板的候選匹配區(qū)域集合在觀測圖上的位置。在這2個模板的匹配結果中,有91.06%的正確匹配區(qū)域包含在該模板對應的候選區(qū)域集合中,而只有71.64%的正確匹配區(qū)域的dcorr最小。按照直接模板匹配的方法,最終的檢測率將很難達到80%。
圖5 觀測圖像和與模板匹配的穩(wěn)定區(qū)域Fig.5 Matching result of the codewords
最后我們將利用RANSAC的思想對圖模型進行配準,并實現(xiàn)最終的目標檢測。因為我們的模型不僅已知,而且參數(shù)確定,同時進行了2次隨機抽樣過程,所以我們的配準方法稱為模型確定二 次隨機抽樣一致性算法MDQRANSAC(model determinate quadratic RANdom sample consensus)。
對于MDQRANSAC算法,模型參數(shù)不需要通過隨機選取的數(shù)據(jù)去擬合,它是用來判斷一組隨機選取的數(shù)據(jù)是否滿足約束的。而2次隨機抽樣過程分別是對模板集合的隨機抽樣和對模板的候選正確匹配區(qū)域集合的隨機抽樣。
我們首先隨機選取2個模板,然后從選取的模板所對應的候選正確匹配區(qū)域集合中隨機選取區(qū)域。這樣可以選出2個區(qū)域,分別對應2個模板。然后我們利用模板的幾何位置,計算選取的穩(wěn)定區(qū)域的幾何位置誤差epos。
由于模板的幾何位置采用極坐標表示,因此θi和ρi是第i個模板的位置,而θji和ρji是與第i個模板對應的穩(wěn)定區(qū)域j的位置。幾何位置誤差ej,posi反映了 MRF模型中的結構約束的大小。如果ej,pos大于一個i事先確定的閾值tposi,區(qū)域j將被認為不滿足MRF模型,需要重新選取i和j加入初始候選集合Scorr。
接下來根據(jù)MRF模型可以得到計算候選集合與圖模型的歸一化誤差enormal的公式(14)
其中Ncorr是Scorr中區(qū)域的個數(shù),α是一個保證歸一化誤差最大值為1的約束量,顯然,α與tpos有關。我們使用下面的公式計算α,Ntemplate是圖模型中包含的模板個數(shù)。
在(14)式中,之所以使用 exp(-4.0(1-dj,corri))作為ei,pos的權值,是因為當 dj,corri= 0時,該模板匹配到了正確的區(qū)域。但是在實際檢測中,正確的匹配區(qū)域對應的epos不一定最小,導致enormal也不一定最小。因此,需要在epos之前乘以一個權值exp(-4.0(1-dj,corri)),使得已經(jīng)確定的正確匹配區(qū)域得到enormal最小。
接下來與原RANSAC算法相同,隨機選擇模板再隨機選擇對應區(qū)域,不斷有區(qū)域被加入Scorr。如果enormal<ebest,那么Sbest將被更新為當前的Scorr;如果ebest小于一個閾值tpos或者到達循環(huán)次數(shù)的上限,算法停止;如果不規(guī)定循環(huán)次數(shù)上限且沒有ebest<tpos,那么算法將退化為窮舉算法,遍歷所有的對應。
當算法停止后,此時在Sbest中得到的將是與圖模型的距離最近的配準結果;如果Sbest不存在,就認為觀測圖像中的目標不存在。另外,在MDQRANSAC中,并沒有對滿足模型參數(shù)的區(qū)域個數(shù)的要求,只要存在穩(wěn)定區(qū)域被添加進候選集合,算法就會計算enormal。
顯然,我們提出的算法也可以被應用在其他模型確定且存在多個對應,要求尋找最佳匹配的研究中,只是模型參數(shù)、各種誤差約束的計算方法等,根據(jù)應用的不同會有所改變。
在我們的實驗中,參考圖像為2幅可見光俯視圖,而觀測圖像則為4組紅外光側視圖序列,總共1 573幅圖像,每2組序列對應1幅參考圖像。其中的一幅參考圖像和對應的1幅觀測圖像如圖1所示。此外,數(shù)據(jù)中還包括與每幅觀測圖像相對應的攝像機的視角參數(shù)。
根據(jù)我們的方法,給出了最終的目標檢測結果。每個參考圖像都利用MSER提取了5個模板用來建立圖模型,平均每個模板對應10個左右的穩(wěn)定區(qū)域。如果遍歷所有可能的組合,大約需要循環(huán)6 631 000次;但在MDQRANSAC算法中,規(guī)定了循環(huán)次數(shù)的上限為1000次。Intel Xeon在3.20GHz CPU,5.98GB內(nèi)存的計算機上使用單線程處理,MDQRANSAC算法的平均處理時間為15ms。
首先給出使用我們的方法與直接利用模板進行區(qū)域匹配的方法的比較。直接利用模板計算相關的方法將首先計算每個模板與觀測圖像的互相關矩陣,然后選取觀測圖像上與模板相關度比較大的幾個對應位置,最后利用MDQRANSAC按照幾何位置關系搜索最終匹配結果,如表1所示。
表1 本文的方法與直接模板匹配的比較Table 1 The comparison of the method used in this study and model matching
可以看到,利用MSER提取物體穩(wěn)定區(qū)域的方法克服了噪聲的影響,獲得了物體準確的形狀描述,提高了最終結果的識別率。而在原觀測圖像上直接進行模板匹配的方法受到噪聲的干擾以及灰度變化的影響,識別率較差。這從一個方面說明了利用 MSER提取穩(wěn)定區(qū)域方法的有效性,能夠提高區(qū)域正確匹配的準確率。
為了說明區(qū)域特征的優(yōu)點,接下來給出的是采用區(qū)域特征的實驗結果和采用邊緣特征的實驗結果對比,如表2所示。其中邊緣特征的結果是用觀測圖像的邊緣圖和模板的邊緣,通過計算Hausdorff距離得到。為了保證實驗條件的一致,區(qū)域特征的方法也直接利用模板在觀測圖上計算互相關得到。最終的結果包含了從參考圖上提取的全部模板的匹配結果。
表2 區(qū)域特征匹配與Hausdorff的方法比較Table 2 The comparison of the region feature matchingand the Hausdorff
從實驗結果來看,基于區(qū)域特征的方法有效克服了存在的噪聲以及視角變換不準確帶來的誤差。而基于邊緣特征的方法則沒有取得很好的實驗結果。因為觀測圖像的清晰度很低,存在大量噪聲邊緣,同時形狀的不完全一致也影響了匹配的性能。序列2基于邊緣的結果要好于基于區(qū)域的結果是因為存在一個模板的匹配準確率只有5.56%,拖累了本文方法最終的匹配結果。
表3最后給出采用將觀測圖像進行視角變換與將參考圖像進行視角變換的檢測結果對比。
表3 不同圖像視角變換的比較結果Table 3 The comparison of different view transformations
從實驗結果來看,對觀測圖像進行視角變換得到的結果的識別率更高。因為變換到俯視的視角上,模板與對應的穩(wěn)定區(qū)域之間只存在簡單的旋轉和尺度誤差,對區(qū)域特征的匹配結果影響較?。欢儞Q到側視視角上,模板與穩(wěn)定區(qū)域之間存在仿射畸變,對最終匹配結果的影響較大。其中,序列2對參考圖像進行變換后,由于視角誤差,變換后的區(qū)域的朝向與觀測圖像明顯不一致,導致了最終的識別率很低,如圖6所示?;野咨珵槟0逶谟^測圖像中的匹配位置。
圖6 序列2的匹配結果Fig.6 Matching result of Group 2
而序列3和序列4由于攝像機的視角參數(shù)相對比較準確,因此得到了較高的檢測率(80%以上)。
圖7是使用我們的方法得到的一個檢測結果(灰白色區(qū)域)。從參考圖像中提取得到的模板總共有5個,但是最終的檢測結果中,僅僅依靠其中3個,就保證了觀測圖像與圖模型的距離小于閾值。
本文通過匹配地面物體區(qū)域,將三維目標的匹配問題轉化為了二維的圖匹配問題,為解決不同視角下具有復雜背景的異質(zhì)多模圖像的目標檢測,提供了一個新的方法。
圖7 檢測結果Fig.7 Target detection result
從實驗結果對比來看,本文基于物體區(qū)域特征匹配的方法確實獲得了比邊緣特征更加優(yōu)異的檢測性能。同時采用MSER來提取穩(wěn)定區(qū)域,也有效克服了噪聲的干擾,獲得了穩(wěn)定的區(qū)域提取結果,提高了最終結果的準確率。
另外,從表1和表2的結果來看,采用MRF模型檢測目標,擺脫了復雜背景的干擾,比直接使用模板進行匹配的方法也獲得了更好的性能。
本文的方法也能夠被用在其他不同視角且具有復雜背景的異質(zhì)光圖像或同質(zhì)光圖像的目標檢測問題中,只要參考圖像和觀測圖像中包含能夠被穩(wěn)定提取的物體區(qū)域。
當然,本文仍然存在一些不足,主要就是在計算匹配區(qū)域的距離中,我們使用了最簡單的互相關的方法。盡管該方法在實驗數(shù)據(jù)集上獲得了很好的匹配結果,但是我們?nèi)匀恍枰粋€更加魯棒的形狀匹配算法。另外,本文是對地面物體區(qū)域進行配準,沒有考慮高度信息無法忽略的三維目標。
我們的下一步工作是消除地面目標區(qū)域對結果的影響,將我們的方法完全擴展到三維目標的檢測和匹配上,同時改進我們的區(qū)域匹配算法。
[1]Zitova B,F(xiàn)lusser B.Image registration methods:a survey[J].Image and Vision Computing,2003,21:977-1000.
[2]Hild M,Umeda G.Image Registration in Stereo Based Multimodal Imaging Systems[C]//In Proceedings of the 4th International Symposium on Image and Signal Processing and Analysis.Los Alamitos:IEEE Computer Society,2005:70-75.
[3]馬俊,曹治國.基于邊緣信息的紅外與可見光圖像匹配技術[J].計算機與數(shù)字工程,2006,34(12):30-32.
[4]曾文峰.紅外與可見光圖像融合中的快速配準方法[J].紅外與激光工程,2002,31(2):158-160.
[5]陶冰潔,王敬儒,張啟衡.采用仿射變換的紅外與可見光圖像配準方法[J].光電工程,2004,31(11):39-41.
[6]Li H H,Zhou Y T.Automatic visual/IR image registration[J].Optical Engineering,1996,35(2):391-400.
[7]Li H,Manjunath B S,Mitra S K.A Contour Based Approach to Multisensor Image Registration[J].IEEE Trans of Image Processing,1995,4(3):320-334.
[8]Huang X S,Chen Z.A wavelet-based Multisensor image registration algorithm[C]//In Proceedings of 6th international conference on Signal Processing.Los Alamitos:IEEE Computer Society,2002:773-776.
[9]Cheng H,Zheng S,Yu Q Z,et al.Matching of SAR Images and Optical Images Based on Edge Feature Extracted via SVM[C]//In Proceedings of the 7th international Conference on Signal Processing.Los Alamitos:IEEE Computer Society,2004:930-933.
[10]Krotosky S J,Trivedi M M.Mutual information based registration of multimodal stereo videos for person tracking[J].Computer Vision and Image Understanding,2007,106(2/3):270-287.
[11]Matas J,Chum O,Urban M,et al.Robust widebaseline stereo from maximally stable extremal regions[J].Image and Vision Computing,2004,22(10):761-767.
[12]Fischler M A,Bolles R C.Random Sample Consensus:A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography[J].Communications of the ACM,1981,24(6):381-395.
[13]Huttenlocher D P,Klanderman G A,Rucklidge W J.Comparing Images Using the Hausdorff Distance[J].IEEE Trans.On Pattern Analysis and Machine Intelligence,1993,15(9):850-863.