陳彥百， 郭長勇， 程遠志

(哈爾濱工業(yè)大學 計算機科學與技術學院， 哈爾濱 150001)

引言

在骨關節(jié)炎的發(fā)病過程中，患者關節(jié)軟骨的厚度以及形態(tài)會隨著病情的發(fā)展而發(fā)生變化。通過MR圖像測定骨關節(jié)炎患者的關節(jié)軟骨厚度變化以及對其進行形態(tài)描述，已經(jīng)成為骨關節(jié)炎的診斷和治療中最廣泛使用的生物標志[1-3]。研究顯示，在一定時間段內(nèi)軟骨厚度的變化屬于局部變化，而不是整體的厚度變化[4-5]。

為了觀測關節(jié)炎患者病情變化，本文從MR圖像中跟蹤關節(jié)軟骨微小部位的厚度變化來展開預測。首先要采集不同時間掃描關節(jié)的MR圖像；然后對關節(jié)軟骨進行分割。在骨關節(jié)炎的發(fā)病過程中，關節(jié)軟骨外表面受到磨損，而在大部分發(fā)病階段，軟骨-股骨面依然保持原有正常形狀和特點，并未受到磨損[6]。因此可以使用軟骨-股骨面作為參照基準，先得到骨-軟骨交界面兩個點云(集)，然后計算出該點的關節(jié)軟骨厚度，最后計算各點的軟骨厚度差，但前提是知道兩個點云中各點的對應關系，因此研究隨即演變?yōu)樾枰蠼?個點云中各點的對應關系，也就是轉(zhuǎn)化成了一個點對應問題。

1 點對應問題的數(shù)學描述

首先給出點對應問題的數(shù)學描述：

設A、B分別為n維空間Rn的2個有限集合，A={x1,x2,…,xp}，B={y1,y2,…,yq}，A0?A。點對應問題就是找一個從A0到B的映射P(或一個從A到B的部分映射P)，使：

?xi∈A0P(xi)=yj

(1)

但其中，A0與A0中的xi所對應的yj均是未知的。

已知映射是一個特殊的二元關系，而二元關系可以用關系矩陣表示，因此，從A到B的部分映射P可以用一個p×q關系矩陣來表示，該矩陣在點對應問題中一般稱為分配矩陣，當分配矩陣為方陣時稱為配置矩陣，其中:

(2)

說明：

(1)從點對應問題的數(shù)學描述來看，找到的對應可能存在一個元素對應多個元素的情況。但在實際問題中，往往希望這種對應是一對一的，而且可能要求A(或B)中的元素都存在對應，這種對應就是從A到B單射或從B到A的單射。如果2個集合A，B中的元素一樣多，則這種對應就是從A到B的雙射。

(2)很容易證明以下結(jié)論：

① 從A到B的部分映射P總共有(q+1)p，從A到B的映射P總共有qp個。

② 如果pq，從B到A的單射P有p!/(p-q)!個。

③ 如果p=q，則從A到B的雙射P總共有p!。

雖然從一個有限集合到另一個有限集合的部分映射、映射、單射和雙射的個數(shù)都是有限的，但不要試圖通過窮舉法一一列出映射P，并驗證其是否符合要求來求解點對應問題，因為還未生成有效的生成算法將所有可能的部分映射、映射、單射窮舉出來。對于雙射，雖能將雙射的窮舉轉(zhuǎn)化成全排列的生成問題，進而利用全排列的生成算法，如：字典序法、遞增進位制數(shù)法、遞減進位制數(shù)法、鄰位對換法等完整地將所有可能的全排列無重復、無遺漏地窮舉出來，但當p或q較大時其計算量巨大，以致在允許的時間內(nèi)無法得到最終優(yōu)化解。

考察一個映射的映射圖：設A={a,b,c,d}，B={u,v,x,y,z}，f是一個從A到B的一個映射：f(a)=y，f(b)=z，f(c)=u，f(d)=x，具體如圖1所示。

圖1 函數(shù)圖

從圖論的角度看，圖1就是一個頂點集V=A∪B。進一步分析，不難得出：尋找A到B的一個(部分)映射就是從A∪B為頂點集的偶圖中尋求匹配。于是，點對應問題就轉(zhuǎn)化為偶圖的匹配問題。

定義1G=(V,E)稱為偶圖，如果G的頂點集V有一個二劃分{V1,V2}，使得G的任一條邊的2個端點在V1中，另一個在V2中，偶圖記做(V1∪V2,E)。

如果?u∈V1,v∈V2均有u,v∈E，則這個圖稱為完全偶圖，并記為Km,n，其中|V1|=m, |V2|=n。

定義2設G=(V,E)是一個圖，則：

1)G的任2條不相鄰的邊x與y稱為相互獨立的。

2)若Y?E，且Y中任2條邊都是相互獨立的，則稱Y為G的一個匹配。

3)設Y是G的一個匹配，若對G的任一匹配Y′，恒有|Y′||Y|，則稱Y是圖G的一個最大匹配。

定義3設G=(V1∪V2,E)是一個偶圖，若存在一個匹配Y使得|Y|=min{|V1|,|V2|}，則稱Y是偶圖G的完全匹配；若|V1|=|V2|，則稱Y為G的一個完美匹配。

相應地，找一個使對應點盡可能多的部分映射就是找最大匹配，找單射就是找完全匹配，找雙射就是找完美匹配。在實際問題中，最常見的就是涉及最大匹配和完全匹配的問題，當然，一般是找滿足一定(最優(yōu))條件的最大匹配和完全匹配。

偶圖的最大匹配有2種求法，分別是：網(wǎng)絡最大流算法和匈牙利算法。匈牙利算法本質(zhì)上還是網(wǎng)絡最大流算法，只是將依據(jù)偶圖匹配這個問題的特點，把最大流算法做了簡化，提高了效率。最大流算法的核心問題就是找增廣路徑，匈牙利算法也不例外。

一般情況下，找到一個偶圖的最大匹配并不意味著點對應問題已經(jīng)解決。通常一個偶圖的最大匹配并不只有一個。多數(shù)情況下，點對應問題即會轉(zhuǎn)化為一個優(yōu)化問題，因此現(xiàn)在問題已轉(zhuǎn)化為求取帶權偶圖的最大匹配，帶權偶圖的最大匹配可以用KM(Kuhn-Munkres)算法求解，也可以用文獻[7-9]等給出的圖匹配算法實現(xiàn)求解。

多數(shù)情況下，點對應問題需要求解的是滿足一定(最優(yōu))條件的完全匹配。本文正是將求解骨-軟骨交界面兩個點集之間對應問題轉(zhuǎn)化為求解偶圖的完全匹配問題。

設A，B分別表示前后兩次掃描成像后分割提取的骨-軟骨交界面的點集，A={x1,x2,…，xp}，B={y1,y2,…，yq}。將A，B作為偶圖的頂點，A中每個點與B中每個點均可能存在對應關系，因此A中每個點與B中每個點之間均可以連線作為偶圖的邊，這樣所得到的偶圖為完全偶圖G=(A∪B,E)?，F(xiàn)在需要求解滿足一定條件的完全匹配。

首先利用偶圖匹配解決完全匹配的存在性問題。

定理1Hall定理設G=(V1∪V2,E)是一個偶圖，則|V1|≤|V2|。G中存在從V1到V2的完全匹配充分必要條件是V1中任意n個頂點(n=1,2,……，|V1|)至少與V2中n個頂點相連接。

不難驗證，完全偶圖符合Hall定理條件，因此完全匹配肯定存在。在此，求解滿足一定條件的完全匹配，繼續(xù)將該問題轉(zhuǎn)化為如下的優(yōu)化問題：

(3)

(4)

(5)

pij=1 orpij=0

(6)

其中，cij表示xi與yj之間的某種度量(如距離)，約束條件(4)要求yj有且只有一個xi與之對應；約束條件(5)要求xi有且只有一個yj與之對應。滿足條件(4)、(5)、(6)的解pij可以寫成矩陣形式，該矩陣就是前面提到的配置矩陣。

上述優(yōu)化問題實質(zhì)上可以視為整數(shù)規(guī)劃或0-1規(guī)劃中分派問題或指派問題的數(shù)學模型，用Munkres分配算法來設計研究求解[10-11]。

2 點對應算法

設A，B分別為Rn的2個集合，A={x1,x2, …,xp}，B={y1,y2, …,yq}。用Munkres分配算法求解點集A與點集B之間的對應的步驟如下：

(1)計算距離矩陣D

(7)

如果p≠q，則要通過額外增加行或列的方式保證矩陣D是一個方形矩陣。一般情況下，額外增加行或列的每個元素的值取矩陣D中最大值。

(2)矩陣D中的每一行均減去該行的最小元素。

(3)在第(2)步的基礎上，再考慮矩陣每一列，對沒有零元素的列，減去該列的最小元素。

(4)用直線覆蓋零元素。用最少的直線覆蓋所有零元素的所在的行和列，如果所用最少直線條數(shù)等于矩陣的行數(shù)，則跳至第(8)步。

(5)在第(3)步得到的矩陣中，求取未被覆蓋元素的最小元素，然后將最小元素加至每一個覆蓋元素上。如果一個元素被覆蓋兩次(即有橫豎兩條直線覆蓋元素)，則該元素加兩次最小元素。

(6)考慮第(5)步得到的矩陣。矩陣的所有元素均減去最小元素(仍為第(5)步所求的未被覆蓋元素的最小元素)。

(7)再次用最少的直線覆蓋零元素。如果所用最少直線條數(shù)小于矩陣的行數(shù)，則跳至第(5)步。

(8)每行或每列僅選一個零元素。然后在矩陣中將選中的零元素改為1，其它元素都改為0，這樣就得到僅含有一個元素的0、1矩陣。

(9)考慮第(8)步所得的元素0、1的矩陣。從矩陣中刪除第(1)步所添加的空行或空列，這樣就得到了點集A和B之間的對應矩陣M，其中M是一個p×q矩陣。最后，點集A中的點xi通過公式(8)在B中找到與之對應的元素。具體公式表述如下：

(x1,x2, …,xp)=(y1,y2, …,yq)MT

(8)

本文求解點之間的對應方法可闡釋解析為如下兩步：

(1)對點集A與B實施主軸變換，得到新的點集A′與B′。

(2)對點集A′與B′，使用Munkres分配算法求解點之間的對應。

3 實驗結(jié)果與分析

為了評估研究提出的點對應算法的性能，本文采用30個豬腿骨和30名膝關節(jié)炎患者的MR圖像數(shù)據(jù)。每個豬腿骨數(shù)據(jù)有18個標記點，每個膝關節(jié)數(shù)據(jù)都有8個標記點，并均將用于算法性能評估。

3.1 MR圖像數(shù)據(jù)獲取

將30個新鮮冷凍的豬腿骨MR圖像作為實驗數(shù)據(jù)。對于每個豬腿骨，20～22根直徑為0.2 mm的竹牙簽穿過軟骨插入股骨表面后進行兩次磁共振成像，視像效果如圖2所示。在第二次成像前需對豬腿骨重新定位，然后從MR圖像中選出15個標記點，如圖3所示。標記點使用圖形用戶界面(GUI)軟件來展開選擇設計，該軟件可以將MR圖像的矢狀面、軸狀面和冠狀面直觀呈現(xiàn)在顯示器上，還具有旋轉(zhuǎn)、平移(上/下/左/右/前/后)、縮小和放大、調(diào)整窗口大小等功能。用戶可以逐片瀏覽MR切片，滾動縮小和放大標記。通過使用GUI，每根牙簽標記的中心位置得到了最終確定。單人間隔一周重復6次結(jié)果取平均，得到每根牙簽標記的中心位置的坐標，誤差為0.012±0.002 mm。此外，研究選擇的標記經(jīng)過放射科醫(yī)生的公正確認。

圖2豬腿骨

Fig.2Porcineknee

圖3 帶有標記點的豬腿骨MR圖像

30名患者的MR圖像數(shù)據(jù)獲取。30名患者信息如下：11名男性，19名女性，年齡范圍為32～67歲，平均年齡為52.6歲。每名患者2次掃描成像，中間間隔12個月。使用GUI軟件，結(jié)合文獻[9]方法確定8個標記點的坐標，實際位置即如圖4所示。為了使標記點的坐標更趨精確，單人花費2周時間重復6次，結(jié)果取平均作為標記點的坐標，誤差為0.012±0.002 mm，研究選擇的標記也經(jīng)過放射科醫(yī)生的公正確認。

圖4 帶有標記點的膝關節(jié)點MR圖像

成像條件：1.5T掃描機(GE Medical Systems，Milwaukee，WI)，三維脂肪抑制快速梯度回波(SPGR)成像序列；成像參數(shù)：TR/TE=9.0/4.2；flipangle=30°；FOV=320×320 mm；sectionthickness=2 mm；matrix=256×256；scantime=4:30。

3.2 關節(jié)軟骨的分割

關節(jié)軟骨的分割主要采用Lynch等提出的一種半自動分割方法[12]，下面則針對該方法概述為如下3個步驟，具體內(nèi)容可參閱文獻[12]。

(1)手動初始化。在每張MR圖像切片的軟骨末端標記兩個點，標點位置則詳見圖5(a)；在關節(jié)軟骨的中心MR圖像切片上標記六個點。

(2)檢測軟骨。計算機優(yōu)化這些點的位置，將其外推到相鄰切片中，詳情可見圖5(b)，并自動調(diào)整點使其位于在每個切片的軟骨線(cartilage sheet)上，操作效果即如圖5(c)所示。接下來，計算機能提供覆蓋股骨髁的最佳估計值軟骨的3D程度的運算支持，圖像輸出則如圖5(d)所示。在這個過程中，可以人工糾正位置出錯的點。

(3)檢測軟骨內(nèi)、外交界面。通過人工校正，計算機能調(diào)整軟骨線至軟骨的內(nèi)邊緣(骨與軟骨的交界面)和外邊緣(軟骨與滑膜的交界面)，運行結(jié)果可如圖5(e)所示。

圖5 關節(jié)軟骨的分割

3.3 算法性能的評估

下面將利用上面獲取的數(shù)據(jù)，通過與其它三種求解點對應方法的比較研究來評估本次提出算法的性能。這三種方法分別是：Besl等提出的ICP算法中確定點之間對應的最鄰近方法(一個點集中的點與另一個點集中到其距離最近的點相對應[13])；Sanroma等提出的一種用來解決點對應問題的圖匹配算法[8]；Guo等提出的基于稀疏特性來解決點對應問題的圖匹配算法[9]。

對于每個豬腿骨數(shù)據(jù)，將其第一次掃描成像的MR圖像中提取出的18個標記點作為參考點，與第二次掃描成像的MR圖像中提取出所有點分別用4種方法進行匹配，這樣可以得到18對對應點，然后利用公式(9)計算這些對應點間的歐氏距離之和，然后對30個豬腿骨數(shù)據(jù)得到的距離之和取平均，可得運算結(jié)果如表1所示。

(9)

其中，(xAi,yAi,zAi)(i=1,2,…,18)表示參考點的18個標記點坐標，(xBi,yBi,zBi)(i=1,2,…,18)表示與這18個標記點對應的點坐標。

對于每個膝關節(jié)數(shù)據(jù)進行相似處理，最終可得結(jié)果如表1所示。由表1可以看出，本次研究給出方法產(chǎn)生的匹配誤差是最低的。

表1 4種方法的匹配誤差

4 結(jié)束語

為了觀測關節(jié)炎患者病情變化，本文從MR圖像中跟蹤關節(jié)軟骨微小部位的厚度變化來實現(xiàn)研究觀測。研究過程可表述如下：通過對關節(jié)軟骨進行不同時間的掃描，得到MR圖像，再對關節(jié)軟骨進行分割并計算軟骨厚度，而后則是計算各對應點之間的軟骨厚度差，但前提是必須知道2個軟骨點云中各點的對應關系。為此，研究提出了一種解決方法，首先利用一種代數(shù)方法使2個點云在空間對齊，然后使用Munkres分配算法求得點之間的對應關系，仿真結(jié)果表明研究提出的方法配準誤差最低。