專利名稱：基于插值算法的柔性體點分布模型的特征值和特征向量計算方法
技術領域：
本發明涉及柔性體點分布模型分析技術，結合圖形學、矩陣理論、線性代數、醫學圖像分析、計算機圖像處理技術和柔性體模型的參數計算方法。
背景技術：
柔性體相對于剛性體，在計算機中進行建模和分析比較復雜，而在現實生活中物體隨時間柔性變化的情況是非常普遍的，如文獻1王敏琴 韓國強 涂泳秋.基于形變模型的醫學圖像分割綜述[J].醫療衛生裝備，2009，30(2)-37-39；文獻2“Munsell，B.C.；Dalal，P.；Song Wang，″Evaluating Shape Correspondence forStatistical Shape AnalysisA Benchmark Study，″IEEE Transactions on PatternAnalysis and Machine Intelligence，vol.30，no.11，pp.2023-2039，Nov.2008，即馬塞爾，B.C；戴拉，P.；王松，”統計形狀分析中的形狀匹配估計一種校準研究”，IEEE模式分析與機器智能會刊，vol.30，no.11，2008(11)2023-2039；以及文獻3龔永義 羅笑南 賈維嘉 黃厚生.基于改進的彈簧質子模型的醫學圖像配準.計算機學報，2008(7)1224-1233。對于許多生物醫學組織，如人的心臟，雙手，骨骼，腎臟、大腦、細胞等等，這些人們可以很好地描述其外部形狀特征但卻很難給出準確的剛性模型，如文獻4陳昱 莊天戈 王合.直方圖估計互信息在非剛性圖像配準中的應用.計算機學報，2000(4)444-448；文獻5簡江濤 馮煥清 熊進.點分布模型約束的主動輪廓及其在腦MR圖像分割中的應用.中國生物醫學工程學報，2006，25(5)513-517；以及文獻6周壽軍，梁斌，陳武凡.心臟序列圖像運動估計新方法基于廣義模糊梯度矢量流場的形變曲線運動估計與跟蹤.計算機學報.2003.26(11)。點分布模型(PDM)是一種新近崛起的技術，它能很好的對研究對象進行描述，可以成功地應用于非剛性模型的分解與合成處理。近年來，柔性體的統計學建模及其與圖像解釋技術的結合催生了一些高級應用軟件，這些軟件不僅針對生物醫學，還涉及目標跟蹤、識別、以至影視等領域，如文獻7唐果 趙曉東 汪元美.三維醫學圖像分割與可視化研究.計算機學報，1998(3)204-209；文獻8段儕杰 馬竟鋒 張藝寶 侯凱 包尚聯.能量傳導模型及在醫學圖像分割中的應用[J].軟件學報，2009(5)-1106-1115；文獻9侯立華.改進Snake模型在醫學超聲圖像分割中的應用[J].計算機仿真，2008，25(8)-183-185，322；以及文獻10閔小平 王博亮 戴培山 黃曉陽.基于主動形變模型的醫學序列圖像的分割[J].系統仿真學報，2007，19(22)-5331-5335。
PDM是一種基于對象樣本分析的統計模型。從Cootes等人開始，這種基于訓練樣本上標記點的形體建模思想一直被貫徹下來，如文獻11Timothy F.Cootes，Christopher J.Taylor，“Combining point distribution models with shape models basedon finite element analysis”，Image Vision Comput.13(5)403-409(1995)，即蒂姆F.古茨，克里斯托弗J.泰勒，“有限元分析中點分布模型與形狀模型的結合”，計算機圖像與視覺，1995，13(5)403-409。最典型的想法是，使標記點自動對齊以確定圖形的位置和形變類型。對齊圖像后可以很容易的比較同組標記點在不同例圖中的位置變化，只需要計算它們的坐標就可以了。因而，PDM模型對于形狀的描述是通過輪廓點向量來完成的，它由中間形狀和形狀變化樣式組成。通過對這種分布情況建模，我們可以得到與原始訓練集相似的新形狀。
在三維空間，如果訓練集圖像中的每個圖形都由一組數量為n的標記點描述，那么它可以用一個3n維的矢量表示(如果是二維對象即為2n維矢量)。訓練集中的標記工作非常重要，通常是手動完成的。另外，在研究中我們要收集數百名例柔性體的個例以得到比較好的統計結果，如文獻12Guoyan Zheng；XiaoDong；Rajamani，K.T.；Xuan Zhang；Styner，M.；Thoranaghatte，R.U.；Nolte，L.-P.；Ballester，M.A.G.，″Accurate and Robust Reconstruction of a Surface Model of theProximal Femur From Sparse-Point Data and a Dense-Point Distribution Model forSurgical Navigation，″IEEE Transactions on Biomedical Engineering，vol.54，no.12，pp.2109-2122，Dec.2007，即張國艷，董曉，阿雅曼尼，K，T；張璇，斯蒂那，M；索拉納哈迪，R U.；諾蒂，L，P.；巴耶斯特爾，M.A.G.，“手術導航系統中基于稀疏點數據及密集點分布模型的股骨近端精確魯棒重建”，IEEE生物醫學工程會刊，vol.54，no.12，2007(12)2109-2122；文獻13Koikkalainen，J.；Tolli，T.；Lauerma，K.；Antila，K.；Mattila，E.；Lilja，M.；Lotjonen，J.，″Methods of Artificial Enlargementof the Training Set for Statistical Shape Models，″IEEE Transactions on MedicalImaging，vol.27，no.11，pp.1643-1654，Nov.2008，即考克克萊納，J.；托里，T.；勞爾瑪，K.；按提拉，K.；馬蒂拉，E.；莉亞，M.；羅提那，J.，“統計形狀模型訓練集人工擴大方法”，IEEE醫學圖像會刊，vol.27，2008(11)1643-1654。訓練集形成的概率模型維度可能非常的大(在一個心臟模型中通常可達到上萬維)，所以建立PDM的過程是非常枯燥又費時費力的。然而，在一些應用中，單獨的一個PDM模型或是有限個PDM模型都不足以很好的描述對象變化。事實上，不論是連續PDM模型還是模型插入方法都是為了滿足對象的性質，即對象時間和空間上的變化。例如在心臟建模時，需要不同時間段的模型才能很好地描述心室形狀在心臟跳動不同階段的變化。但是，在構建整個心臟運動周期的連續模型時，我們只能由已知的模型序列得到特定幾個時刻的模型，而其他任意時刻的模型就不得而知了，盡管在基于模型的圖像分割領域PDM建模出現已有十來年，但對于模型插值及相關內容的深入研究還沒有開展過。
發明內容
為了克服已有的柔性體點分布模型的不能形成實時分布的模型、模型精度差、實用性差的不足，本發明提供一種能夠形成實時分布的中間過渡模型、模型精度高、實用性良好的基于插值算法的柔性體點分布模型的特征值和特征向量計算方法。
本發明解決其技術問題所采用的技術方案是 一種基于插值算法的柔性體點分布模型的特征值和特征向量的計算方法，所述計算方法包括以下步驟 1)、載入兩個相鄰時刻Ta和Tb的柔性體分布模型Ma和Mb，其計算式為 Ma＝<ya，Φa，λa>，Mb＝<yb，Φb，λb>； 2)、通過所述柔性體分布模型Ma和Mb，分別計算平均形狀ya和yb、特征矩陣Φa和Φb、以及特征值λa和λb，設定柔性體從Ta到Tb呈連續變化， 2.1)中間模型的特征值的計算 λt≈a2λa+b2λb，(25) 這里λa和λb分別為前后兩個相鄰PDM模型的特征值向量，即 λa＝(λa1，λa2，λa3，...)，(26) λb＝(λb1，λb2，λb3，...). (27) 2.2)、中間模型的特征向量的估計 插入中間模型的協方差矩陣近似為 St≈a2Sa+b2Sb (28) 其中，協方差矩陣Sa和Sb由特征矩陣Φa和Φb近似得到，Φa和Φb從插補前已知的點分布模型中得到；設定中間模型建立前先對訓練集數據進行PCA分析以減少參數個數，并選取t種形變模式來描述整個對象的形變；由于Φ是協方差矩陣中最前面的t個特征矢量，則有 其中，Λa和Λb是由特征矢量λa和λb組合形成的對角矩陣；用奇異值分解的方法從式得到插入中間模型的特征矢量 由兩個已知的相鄰PDM模型及相關參數a、b、Λa、Λb、Φa和Φb確定正交矩陣Ut，所述正交矩陣Ut的前t行就是特征矩陣Φt。
進一步，所述步驟2)中還包括 2.3)、中間模型的協方差矩陣由下式推導得出 令 因為 顯然有 根據協方差矩陣，矩陣St由Sa、Sb和Sab計算得到，通過如下公式 用奇異值分解法，可令 通常Φab≠Φba，因為矩陣Sab并非對稱矩陣，Φab和Φba都包含t個矢量； 在插補時，點分布模型保留特征值和特征向量，中間模型的協方差矩陣通過下式計算 再進一步，所述步驟2)中還包括 2.3)、特征向量組的計算 通過式(31)來計算特征向量時，將(31)改寫為 其中是一個r×N維的矩陣，其中r是主成分數目，N是樣本維數，即標定點的數量與標記點維度的乘積； 如果Ui是矩陣G＝LΦa的主要特征向量，對應特征值為λi，則(Φaui)和λi分別為矩陣S＝ΦaL的特征向量和特征值。
由Gui＝λiui，，得 S(Φaui)＝ΦaL(Φaui)＝ΦaGui＝Φa(λiui)＝λi(Φaui)(36)。
本發明的技術構思為PDM方法是一種先驗建模方法，它建立在目標對象的一組樣本(即訓練集)已知的基礎上。這組樣本集可以提供對象的形狀及形變的統計學描述。在實際應用中，這組對象樣本的形狀描述通常是使用輪廓來進行的(如圖片上的一組像素的坐標值)。進一步的，可以在輪廓上選取一些標記點來描述目標物形狀，而這些標記點往往與目標物的某種特征相關聯。
不同樣本上的對應標記點在選取的時候遵循著“特征位置大致相同”的原則，但往往這種定位是不夠精確的，存在一定偏差。這種偏差要歸因于不同個體間的自然形變。我們期望相對于整個形體的尺度而言，這種偏差是十分微小的。PDM方法允許我們模仿這種微小的差異，并且指出哪些偏差的確是微小的，哪些相對而言更大一些。
由于訓練集中的形體是從不同的圖片集中獲得的，而不同圖片集的坐標系并不統一，因此有必要先將所有訓練集形體進行坐標歸一化，即將它們對齊。對齊的過程實際上就是對每一個樣本模型尋找適當的相似性變換的過程，包括平移、縮放和旋轉，以使得不同樣本之間盡可能的接近，另一方面，選定的變換也應該令各個樣本模型與由所有樣本得到的平均模型盡可能接近。為了方便表述，假設我們的樣本只有兩個，歸一化簡化為對齊這兩個樣本。每個樣本的形狀都由一個標記點坐標組成的向量來表示 x＝(x1，y1，z1，...，xn，yn，zn)T， (1) 其中(xi，yi，zi)是某個標記點在三維圖像中的空間坐標。如果是二維醫學圖像，則標記點用(xi，yi)表示。
圖片上標記點的獲取有多種方式，人工手動標記或采用一定得自動標記算法都是可以的。例如，Izard等人提出的一種MRI圖像標記方法，使用統一的算法規則在不同人腦圖片中進行組織結構標記，如文獻14C.Izard，B.Jedynak，and C.Stark，Automatic Landmarking of Magnetic Resonance brain Images，Proc.of theSPIE International Symposium on Medical Imaging，12-17 February 2005，San Diego，California，USA，即C.伊扎德，B.扎第納克，C.斯塔克，“大腦核磁共振圖像的自動標記”，SPIE醫學圖像網絡研討會，2005年2月12-17日，美國，加利福尼亞，圣地亞哥。得到一組經過標記的訓練樣本后，我們需要將它們對齊到一個統一的坐標系下。可以使用廣義對齊算法(GPA)來將訓練樣本對齊，并使得所有樣本形狀到平均模型的距離平方和最小，實際上也就是尋找變換Ti(平移，旋轉，縮放)，使下式最小化 Dm＝∑|m-Ti(xi)|2， (2) 其中 且xi是訓練集中的一個形狀樣本，它是一個3n維的向量表示(在三維空間中有n個標記點的情況下)。
對齊后的訓練集組成了一個三維空間上的點云，可以看作是一個概率密度函數。為了降低運算成本和內存占用，我們采用主成分分析法(PCA)尋找點云中的主元，并以前列最大的少數分量建立模型。這些被選出的分量可以描述對象主要的形變。以心動周期內左心室的模型為例，我們經過主成分分析后的用60個特征向量已經能足夠描述心室的形狀及其變化了。最后，PDM模型可以表達為 x＝x+p1b1+…+ptbt＝x+Φc(4) 其中，x是對齊后的平均形狀向量，Φ是一個3n×t的矩陣，Φ的每一列表示主成分方向上的單位向量，c是由形狀參數組成的t維列向量。可見，經過PCA分析后形狀的維度已經由3n下降到了t。
經過以上步驟我們得到了一個類似于點分布模型的統計學模型，該模型可應用于柔性體的物體建模，以及在新的圖像中定位和分割新個例。在訓練集的學習過程中我們可以得到形狀參數b的變化范圍，只要在該范圍內變化，任意b的值都可以生成合理的新樣本(仿真形體)。通常bi的變化為λj(矩陣Φ中對應第i大的特征值)或者
從樣本集中訓練出一個統計學形態模型后，就可以用PDM來解釋新的圖像個例了。為了將模型與圖像匹配，通常采用一種主動重復迭代的方法。該方法使模型不斷地變形來適應新圖像中目標對象的輪廓。模型形狀的變化受到一些統計數據的約束，也即只能在訓練樣本集中所得出的范圍內進行變化。對于形狀模型，我們還要求圖像輪廓大致位于模型點的附近。可以表示為某點的梯度與輪廓上過該點的法線方向的變化。令統計模型的輪廓沿圖像輪廓不斷移動，同時計算兩輪廓的馬氏距離最小化就可以最終得到真實的邊界位置。因此，要在特定實例中獲得目標柔性對象的形狀需要反復執行以下兩步直到收斂(1)沿模型各點的法線尋找圖像輪廓上相匹配的對應點(使得模型點與對應圖像點的馬氏距離最小)；(2)改變模型位置及形狀參數使得模型各標記點更接近圖像輪廓上找到的對應點。也就是說，模型的擬合過程實際上就是在圖像上搜尋模型點的最相似位置及不斷修正整體幾何變換T及形狀參數c的過程，使之最小化。數學上可以表示為 f＝|X-T(x+Pc)|2＝f(b，Xc，Yc，Zc，s，α，β，γ)(5) ＝|X-T(x+Pc；Xc，Yc，Zc，s，α，β，γ)|2 其中X是在當前迭代時得到的臨時模型。
這是一個尋找最小值的問題，可以通過一些非線性優化的迭代方法求解。最終我們可以確定整體變換參數T，以及個體實例的形狀參數b c＝PT(T-1(X)-x). (6) 近幾年，人們已經就基于PDM的柔性體建模方法的許多相關問題進行了深入的研究，如形體對齊、自動標記、平均形體生成、形變建模、模態分析、圖像分割等等問題，如文獻15蔣曉悅，趙榮椿，一種改進的活動輪廓圖像分割技術，中國圖象圖形學報A輯，2004，9(9)1019-1024；文獻16馬峰，唐澤圣，夏紹瑋.多尺度幾何活動曲線及MR圖像邊界提取.計算機學報.2000.23(8)；文獻17Gilliam，A.D.；Epstein，F.H.；Acton，S.T.，″Cardiac Motion Recovery via ActiveTrajectory Field Models，″IEEE Transactions on Information Technology inBiomedicine，vol.13，no.2，pp.226-235，March 2009，即吉列姆，A.D.；愛潑斯坦，F.H.；埃克森，S.T.，“基于軌跡場模型的心動復原”，IEEE生物醫學信息技術會刊，vol.13，2009(3)226-235。此外，在該思想廣泛地應用于非剛體對象建模的基礎上，并可進一步擴展為諸如活動形體模型(ASM)和活動表現建模(AAM)等以應用于更深入的圖像分析。目前，對PDM模型最成功的表達是通過PCA方法得到的，如文獻18Tobon-Gomez，C.；Butakoff，C.；Aguade，S.；Sukno，F.；Moragas，G.；Frangi，A.F.，″Automatic Construction of 3D-ASM Intensity Models by Simulating ImageAcquisitionApplication to Myocardial Gated SPECT Studies，″IEEE Transactions onMedical Imaging，vol.27，no.11，pp.1655-1667，Nov.2008，即托頓-哥麥茨，C.；巴塔考夫，C.；阿古愛德，S.；薩克諾，F.；莫拉蓋斯，G.；弗蘭基，A.F.，“基于模擬圖像采集的三維ASM密度模型自動重建應用于心臟門閘SPECT圖像研究”，IEEE醫學圖像會刊，vol.27，2008(11)1655-1667。PCA方法非常簡單易懂，如文獻19Engelen，S.，Hubert，M.，Vanden Branden，K.，A comparison of three procedures forrobust PCA in high dimensions，Workshop on Robustness for High-dimensional Data(RobHD2004)，2004，即英格倫，S.，赫伯特，M.，范登布蘭登，K.，“高維度PCA魯棒分析的三個過程對比”，高維數據魯棒性學術討論會，2004，首先把一個形體投影到經過學習PCA空間(由線性無關方向形成的正交空間)，然后從訓練樣本中提取少量與形變模式相關的參數(稱為主要成分)，再用這些參數來控制形體的變形。顯然，這種形變是沿著樣本集所確定的最大形變方向進行的。由于只研究幾個主要方向上的形變，無論是計算復雜度還是對存儲空間的要求都大大降低了。從數學角度講，首先是從訓練集中估計得到平均形狀y和形體協方差矩陣S 這樣，一個合理的形狀可以被表達為 y＝y+Φc (8) 式中向量c由對應于各個特征向量的形變參數組成。矩陣Φ則包含了一組由協方差矩陣S分解得到的特征向量，這組向量可線性組合成任意的新形狀y。根據PCA理論，c中的每個參數都控制形體在某一方向上的形變，并且，各個方向都彼此獨立，參數按照形體在對應方向上形變由大到小順序排列。本發明考慮在一個PDM模型僅包含平均形狀、特征值和特征向量而不包括原始訓練集的數據的情況下進行插值運算。即 M＝<y，Φ，λ>(9) 如上所述，單個PDM模型可以由一組柔性對象的形體獲得，例如，我們研究的醫院患者的心臟圖片(一些被觀測的柔性體個例)。然而，有時我們要觀測這些柔性體隨時間的變化情況，這樣可以得到一系列的形狀信息，例如心臟在整個活動周期內的三維變形。而醫院現有設備只能獲得若干時刻的三維圖像，即對應于某些Ti時刻的模型Mi。
對一個特定的動態形變物體，假設我們已經在離散的時間坐標(T1，T2，...)上建立起一組PDM模型(M1，M2，...)。我們可以把問題描述為給出已知的兩個PDM模型，Ta時刻的模型Ma和Tb時刻的模型Mb(模型已知意味著我們知道了平均形狀ya和yb、特征矩陣Φa和Φb、以及特征值λa和λb)，如何得到在t(Ta≤t≤Tb)時刻的中間模型Mt的相關參數。這些參數有平均形狀、特征值和特征向量，即Mt＝<yt，Φt，λt>(圖1)。這里我們假設在Ta到Tb的時間段中，一個形體從Yai連續變化到Ybi，二者都是n維的，如Yai＝(xai1，xai2，xai3，...，xain)T。這里還定義了兩個常數 b＝1-a(11) a和b反映了中間模型相對于Ma和Mb的距離。
(1)中間模型的平均形狀 根據以上的定義和假設，用線性插值(可滿足大部分的實際應用)得到Yai和Ybi之間的中間模型形狀為 yti＝ayai+bybi(12) a和b由式(10)和(11)定義，即a+b＝1。
插入的平均形狀為 所以，我們得到以下結論 定理1線性插值時，插入的PDM中間模型平均形狀是相鄰模型平均形狀的線性組合，其參數由式(13)確定。
(2)中間模型的形變 相似的，中間PDM模型的形變可以這樣算出 Δyti＝yti-yt＝(ayai+bybi)-(aya+byb) ＝a(yai-ya)+b(ybi-yb) (14) ＝aΔyai+bΔybi 從上式可以得到以下結論 定理2線性插值時，插入的PDM中間模型的形狀變化速度是相鄰模型形變速度的線性組合，其參數由式(14)確定。
(3)中間模型的協方差矩陣 中間模型的協方差矩陣可以由下式推導得出 令 因為 顯然有 最終，得到以下推論。
推論1線性插值時，中間模型的協方差矩陣可由以下表達式得到 這里Sab表示模型Ma和Mb的交變矩陣。我們知道協方差矩陣是通過對同一時間段內不同個體的采樣數據進行統計分析得出的，而交變矩陣則描述了在不同采樣時間得到的兩個對象模型的形變過程。
中間模型的估算 (1)交變矩陣 為方便我們的討論，假設模型是n維的，如 Δyai＝[da1i da2i ... dani]T (21) 則，交變矩陣為 通常情況下，樣本是正態或平均分布的，交變矩陣中的各元素趨于零，即 所以，交變矩陣對角線元素遠小于矩陣Sa和Sb中的對應元素，即 (2)特征值的估算 可以看出，交變矩陣Sab和Sba對于中間模型特征值的計算影響非常小，因此我們可以估計 λt≈a2λa+b2λb， (25) 這里λa和λb分別為前后兩個相鄰PDM模型的特征值向量，即 λa＝(λa1，λa2，λa3，...)， (26) λb＝(λb1，λb2，λb3，...).(27) 對λt中元素估計的精度主要取決于樣本集的數量和分布情況。經過計算機模擬發現，當樣本集為30個隨機均勻分布時，誤差率可維持在百分之五以下。
(3)特征向量的估計方法 當樣本集很大時，插入中間模型的協方差矩陣可以近似為 St≈a2Sa+b2Sb(28) 所以，協方差矩陣Sa和Sb可以由特征矩陣Φa和Φb近似得到，而Φa和Φb是可以從插補前已知的點分布模型中得到。假設PDM模型建立前先對訓練集數據進行PCA分析以減少參數個數，并選取t種形變模式來描述整個對象的形變。由于Φ是協方差矩陣中最前面的t個特征向量，我們知道 其中Λa和Λb是由特征向量λa和λb組合形成的對角矩陣。然后，我們可以用奇異值分解(SVD)的方法從式得到插入中間模型的特征向量 現在，已由兩個已知的相鄰PDM模型及相關參數a、b、Λa、Λb、Φa和Φb，確定了正交矩陣Ut，且Ut的前t行就是特征矩陣Φt。
(4)特征向量的精確計算方法 通過上面的公式，我們可以由兩個已知PDM模型計算得到插入中間模型的各個參數(如平均形狀、特征值、特征向量等)。然而，由于我們忽略了交變矩陣的作用，因此這種方法存在一定的誤差，尤其是當樣本集容量很小或其分布比較異常時這種誤差就更明顯。事實上我們有一種精確的方法，就是直接生成訓練集中所有相鄰兩個PDM模型的交變矩陣。根據式(20)，矩陣St可由Sa、Sb和Sab計算得到，這意味著矩陣Sab需要在構建Ma和Mb時就得到。實際上，可以通過如下公式 用奇異值分解法，可令 通常Φab≠Φba，這是因為矩陣Sab并不是一個對稱矩陣。Φab和Φba都包含t個矢量。
在插補時，由于點分布模型只保留了特征值和特征向量的信息，中間模型的協方差矩陣應通過下式計算 相似的，中間模型的特征向量可以由式(31)得到。這種方法由于不受訓練樣本分布的影響，因此是一種更加準確的解決方案。
(5)特征向量組的計算 當直接通過式(31)來計算特征向量時，我們會發現由于矩陣通常很大，計算過程對存儲空間和運算效率都有很高的要求。例如，我們實驗中有一個三維心臟模型，包含2848個標記點，通過主成分分析后選擇了90種形變模式，那么該模型的一個協方差矩陣將占據584MB的內存空間，連最簡單的加減操作可能也要用數分鐘時間。因此，為了簡化計算復雜度，我們將(31)改寫為 其中是一個r×N維的矩陣，其中r是主成分數目，N是樣本維數，即標定點的數量與標記點維度的乘積。下面我們來描述簡化后的計算方法。
定理3如果Ui是矩陣G＝LΦa的主要特征向量，對應特征值為λi，則(Φaui)和λi分別為矩陣S＝ΦaL的特征向量和特征值。
其證明相當簡單，由Gui＝λiui，，得 S(Φaui)＝ΦaL(Φaui)＝ΦaGui＝Φa(λiui)＝λi(Φaui) (36) 本發明的有益效果主要表現在能夠形成實時分布的中間過渡模型、模型精度高、實用性良好。
圖1是從時間序列中的臨近模型插值到得的瞬時PDM模型的示意圖。
圖2是典型的心動周期分為六個階段的示意圖。
圖3是心臟不同階段特征的六個PDM模型的示意圖。
圖4是插值生成新的PDM模型以進行柔性體分析的示意圖。
圖5是通過插值生成的PDM心臟模型的時間平滑變化的示意圖。
具體實施例方式 下面結合附圖對本發明作進一步描述。
參照圖1～圖5，一種基于插值算法的柔性體點分布模型的特征值和特征向量的計算方法，其特征在于所述計算方法包括以下步驟 1)、載入兩個相鄰時刻Ta和Tb的柔性體分布模型Ma和Mb，其計算式為 Ma＝<ya，Φa，λa>，Mb＝<yb，Φb，λb>； 2)、通過所述柔性體分布模型Ma和Mb，分別計算平均形狀ya和yb、特征矩陣Φa和Φb、以及特征值λa和λb，設定柔性體從Ta到Tb呈連續變化， 2.1)中間模型的特征值的計算 λt≈a2λa+b2λb，(25) 這里λa和λb分別為前后兩個相鄰PDM模型的特征值向量，即 λa＝(λa1，λa2，λa3，...)，(26) λb＝(λb1，λb2，λb3，...). (27) 2.2)、中間模型的特征向量的估計 插入中間模型的協方差矩陣近似為 St≈a2Sa+b2Sb (28) 其中，協方差矩陣Sa和Sb由特征矩陣Φa和Φb近似得到，Φa和Φb從插補前已知的點分布模型中得到；設定中間模型建立前先對訓練集數據進行PCA分析以減少參數個數，并選取t種形變模式來描述整個對象的形變；由于Φ是協方差矩陣中最前面的t個特征矢量，則有 其中，Λa和Λb是由特征矢量λa和λb組合形成的對角矩陣；用奇異值分解的方法從式得到插入中間模型的特征矢量 由兩個已知的相鄰PDM模型及相關參數a、b、Λa、Λb、Φa和Φb確定正交矩陣Ut，所述正交矩陣Ut的前t行就是特征矩陣Φt； 2.3)、中間模型的協方差矩陣由下式推導得出 令 因為 顯然有 根據協方差矩陣，矩陣St由Sa、Sb和Sab計算得到，通過如下公式 用奇異值分解法，可令 通常Φab≠Φba，因為矩陣Sab并非對稱矩陣，Φab和Φba都包含t個矢量； 在插補時，點分布模型保留特征值和特征向量，中間模型的協方差矩陣通過下式計算 2.3)、特征向量組的計算 通過式(31)來計算特征向量時，將(31)改寫為 其中是一個r×N維的矩陣，其中r是主成分數目，N是樣本維數，即標定點的數量與標記點維度的乘積； 如果Ui是矩陣G＝LΦa的主要特征向量，對應特征值為λi，則(Φaui)和λi分別為矩陣S＝ΦaL的特征向量和特征值。
由Gui＝λiui，，得 S(Φaui)＝ΦaL(Φaui)＝ΦaGui＝Φa(λiui)＝λi(Φaui) (36)。
在我們一個關于心臟圖像分析的研究項目中，可以根據心動周期定義一個規范化的時間坐標軸。例如，一個心跳周期可以分為六個生理階段心房收縮、心室收縮、心室射血、心室舒張、心室快速填充及舒張后期休息(圖2)。也有人進一步將心室射血的過程分為射血前期和射血后期。左心與右心的運動過程相似但不完全相同。主要的差別是，左心室與動脈的收縮壓比右心的大三至四倍，而且左右心泵血過程的時間方面也有顯著區別，如文獻20楊柳，汪天富，林江莉等，旋轉掃描超聲心臟圖像的中軸匹配與圖像插值，生物醫學工程學雜志，2004，21(1)28-33，41；文獻21陳強，周則明，王平安，夏德深，帶標記線左心室MR圖像的自動分割，中國圖象圖形學報A輯，2004，9(6)666-673；文獻22羅渝蘭，鄭昌瓊等，基于Snake模型的B超心臟圖像的腔體分割及算法研究，四川師范大學學報自然科學版，2002，25(1)66-69。由于整個心臟的空間形態在一個心動周期內變化很大，而單一的模型很難準確描述動態的心臟運動過程，所以我們需要建立一組模型來描述心動不同階段的心臟形態。
在gated-SPECT的圖像研究中，可參考心電圖信號進行時間規范化。我們希望通過對心臟形狀進行統計后可以得到全部的六個模型，并且這些模型能用來準確地描述心臟的特定狀態，如心臟舒張末期和心臟收縮末期。圖3在標準心電圖曲線中標出了各個階段模型出現的位置。圖中我們已經把壓力值幅度經過了修正，使得波峰處的壓力值為正常人體的收縮壓(120mm汞柱)，所以圖中也表示一條典型的左心室壓力曲線。
如果在正常的心動周期中有6個PDM模型，則式(6)中的Ti(代表某一相位所在的時刻)可以定義為
式中t是獲得gated-SPECT圖像的時間值，Hc是心率的的倒數，
表示不大于某小數的最大整數。
然而，在實際的例子中，常會產生不同數目的相位(如圖4中有8個)，這是由于心臟圖像的采樣通常有固定的頻率，如每100ms(圖4)。那么就有必要通過插值得到采樣時刻的PDM模型，以期得到更好的三維心臟形狀的分割。
柔性體點分布模型插值算法，作者研究過程中分別采用計算機仿真及真實心臟數據進行了大量的實驗。在仿真中，我們可以假設一組九維(或更高)的模型，yi＝[x1i，x2i，...，x9i]T。兩個PDM模型M1和M2的數據集都是隨機產生的。各組數據都含有100個樣本，作為生成點分布模型的訓練集，例如 為生成一個插入模型，我們從S1計算M1的特征值及特征向量，對S2也同樣。假設插入模型的左右距離為a＝0.3，b＝0.7。為了進行對比，我們還計算它們的交變矩陣Sab和相應的特征向量組。實驗結果顯示，本中所述的方法能很好的計算出插入模型的各種參數。由式(31)計算出的特征值誤差小于1.83％，由式(35)計算出的特征值誤差小于0.05％。該仿真實驗的C++程序已由所在實驗室開發完成。
用真實的心臟數據進行了測試。通過灌注SPECT成像法得到的大量圖像數據中，我們預先建立了左心室在不同時刻的五個PDM模型。每個模型由246組病例的真實圖像的訓練得到(圖像主要由巴塞羅那醫院提供)，包括一個平均形狀的VTK表達，一組特征值矢量和特征向量。所有醫學圖像的處理過程和模型的插入運算也由都C++程序實現(研究組共開發了約600MB程序代碼)。實驗時，我們在每兩個相鄰的PDM之間均勻地插入四個中間模型，得到的實驗結果令人非常滿意。圖5顯示通過插值生成的PDM模型隨時間能夠平滑地變化左心室柔性模型。
點分布模型作為一種描述柔性對象的統計學模型，已被證明在很多應用上是非常有效的，尤其是在醫學圖像分析領域。為了分析柔性體的動態變化情況，本發明首創研究了點分布模型的插值算法。研究結果表明，可由PDM的平均形狀、特征值和特征向量等參數生成任意時刻的中間模型。事實上，插入模型的平均形狀及形變就是它相鄰前后模型對應項的線性組合，即yt＝aya+byb，Δyti＝aΔyai+bΔybi。其特征值則可以采用式λt≈a2λa+b2λb計算得到。最后，為求出插入模型的協方差矩陣和特征向量，本發明公開了交變矩陣形式及其相應的計算方法。
權利要求
1、一種基于插值算法的柔性體點分布模型的特征值和特征向量的計算方法，其特征在于所述計算方法包括以下步驟
1)、載入兩個相鄰時刻Ta和Tb的柔性體分布模型Ma和Mb，其計算式為
Ma＝<ya，Φa，λa>，Mb＝<yb，Φb，λb>；
2)、所述柔性體分布模型Ma和Mb，包括計算平均形狀ya和yb、特征矩陣Φa和Φb、以及特征值λa和λb，設定柔性體從Ta到Tb呈連續變化，則
2.1)中間模型的特征值的計算
λt≈a2λa+b2λb，(25)
這里λa和λb分別為前后兩個相鄰PDM模型的特征值向量，即
λa＝(λa1，λa2，λa3，...)，(26)
λb＝(λb1，λb2，λb3，...). (27)
2.2)、中間模型的特征向量的估算
插入中間模型的協方差矩陣近似為
St≈a2Sa+b2Sb (28)
其中，協方差矩陣Sa和Sb由特征矩陣Φa和Φb近似得到，Φa和Φb從插補前已知的點分布模型中得到；設定中間模型建立前先對訓練集數據進行PCA分析以減少參數個數，并選取t種形變模式來描述整個對象的形變；由于Φ是協方差矩陣中最前面的t個特征矢量，則有
其中，Λa和Λb是由特征矢量λa和λb組合形成的對角矩陣；用奇異值分解的方法從式得到插入中間模型的特征矢量
由兩個已知的相鄰PDM模型及相關參數a、b、Λa、Λb、Φa和Φb確定正交矩陣Ut，所述正交矩陣Ut的前t行就是特征矩陣Φt。
2、如權利要求1所述的基于時變插值算法的柔性體點分布模型的特征值和特征向量的計算方法，其特征在于所述步驟2)中還包括
2.3)、中間模型的協方差矩陣由下式推導得出
令
因為
顯然有
根據協方差矩陣，矩陣St由Sa、Sb和Sab計算得到，通過如下公式
用奇異值分解法，可令
通常Φab≠Φba，因為矩陣Sab并非對稱矩陣，Φab和Φba都包含t個矢量；
在插值計算時，點分布模型保留特征值和特征向量，中間模型的協方差矩陣通過下式計算
3、如權利要求1或2所述的基于插值算法的柔性體點分布模型的特征值和特征向量的計算方法，其特征在于所述步驟2)中還包括
2.3)、特征向量組的計算
通過式(31)來計算特征向量時，將(31)改寫為
其中是一個r×N維的矩陣，其中r是主成分數目，N是樣本維數，即標定點的數量與標記點維度的乘積；
如果Ui是矩陣G＝LΦa的主要特征向量，對應特征值為λi，則(Φaui)和λi分別為矩陣S＝ΦaL的特征向量和特征值。
由Gui＝λiui，，得
S(Φaui)＝ΦaL(Φaui)＝ΦaGui＝Φa(λiui)＝λi(Φaui)(36)。
全文摘要
一種基于插值算法的柔性體點分布模型的特征值和特征向量計算方法，包括以下步驟1)載入兩個相鄰時刻Ta和Tb的柔性體分布模型Ma和Mb，其中Ma＝<ya，Φa，λa>，Mb＝<yb，Φb，λb>；2)通過所述柔性體分布模型Ma和Mb，分別計算平均形狀ya和yb、特征矩陣Φa和Φb、以及特征值λa和λb，設定柔性體從Ta到Tb呈連續變化；2.1)中間模型的特征值的計算；2.2)中間模型的特征向量的估算。本發明能夠形成實時分布的中間過渡模型、模型精度高、實用性良好。
文檔編號G06T17/00GK101639948SQ20091010218
公開日2010年2月3日 申請日期2009年8月20日 優先權日2009年8月20日
發明者陳勝勇, 杜雅慧, 秋 管, 毛國紅, 王萬良 申請人:浙江工業大學