專利名稱：基于插值算法的柔性體點(diǎn)分布模型的特征值和特征向量計(jì)算方法
技術(shù)領(lǐng)域：
本發(fā)明涉及柔性體點(diǎn)分布模型分析技術(shù)，結(jié)合圖形學(xué)、矩陣?yán)碚摗⒕€性代數(shù)、醫(yī)學(xué)圖像分析、計(jì)算機(jī)圖像處理技術(shù)和柔性體模型的參數(shù)計(jì)算方法。

背景技術(shù)：
柔性體相對于剛性體，在計(jì)算機(jī)中進(jìn)行建模和分析比較復(fù)雜，而在現(xiàn)實(shí)生活中物體隨時(shí)間柔性變化的情況是非常普遍的，如文獻(xiàn)1王敏琴 韓國強(qiáng) 涂泳秋.基于形變模型的醫(yī)學(xué)圖像分割綜述[J].醫(yī)療衛(wèi)生裝備，2009，30(2)-37-39；文獻(xiàn)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.；王松，”統(tǒng)計(jì)形狀分析中的形狀匹配估計(jì)一種校準(zhǔn)研究”，IEEE模式分析與機(jī)器智能會(huì)刊，vol.30，no.11，2008(11)2023-2039；以及文獻(xiàn)3龔永義 羅笑南 賈維嘉 黃厚生.基于改進(jìn)的彈簧質(zhì)子模型的醫(yī)學(xué)圖像配準(zhǔn).計(jì)算機(jī)學(xué)報(bào)，2008(7)1224-1233。對于許多生物醫(yī)學(xué)組織，如人的心臟，雙手，骨骼，腎臟、大腦、細(xì)胞等等，這些人們可以很好地描述其外部形狀特征但卻很難給出準(zhǔn)確的剛性模型，如文獻(xiàn)4陳昱 莊天戈 王合.直方圖估計(jì)互信息在非剛性圖像配準(zhǔn)中的應(yīng)用.計(jì)算機(jī)學(xué)報(bào)，2000(4)444-448；文獻(xiàn)5簡江濤 馮煥清 熊進(jìn).點(diǎn)分布模型約束的主動(dòng)輪廓及其在腦MR圖像分割中的應(yīng)用.中國生物醫(yī)學(xué)工程學(xué)報(bào)，2006，25(5)513-517；以及文獻(xiàn)6周壽軍，梁斌，陳武凡.心臟序列圖像運(yùn)動(dòng)估計(jì)新方法基于廣義模糊梯度矢量流場的形變曲線運(yùn)動(dòng)估計(jì)與跟蹤.計(jì)算機(jī)學(xué)報(bào).2003.26(11)。點(diǎn)分布模型(PDM)是一種新近崛起的技術(shù)，它能很好的對研究對象進(jìn)行描述，可以成功地應(yīng)用于非剛性模型的分解與合成處理。近年來，柔性體的統(tǒng)計(jì)學(xué)建模及其與圖像解釋技術(shù)的結(jié)合催生了一些高級(jí)應(yīng)用軟件，這些軟件不僅針對生物醫(yī)學(xué)，還涉及目標(biāo)跟蹤、識(shí)別、以至影視等領(lǐng)域，如文獻(xiàn)7唐果 趙曉東 汪元美.三維醫(yī)學(xué)圖像分割與可視化研究.計(jì)算機(jī)學(xué)報(bào)，1998(3)204-209；文獻(xiàn)8段儕杰 馬竟鋒 張藝寶 侯凱 包尚聯(lián).能量傳導(dǎo)模型及在醫(yī)學(xué)圖像分割中的應(yīng)用[J].軟件學(xué)報(bào)，2009(5)-1106-1115；文獻(xiàn)9侯立華.改進(jìn)Snake模型在醫(yī)學(xué)超聲圖像分割中的應(yīng)用[J].計(jì)算機(jī)仿真，2008，25(8)-183-185，322；以及文獻(xiàn)10閔小平 王博亮 戴培山 黃曉陽.基于主動(dòng)形變模型的醫(yī)學(xué)序列圖像的分割[J].系統(tǒng)仿真學(xué)報(bào)，2007，19(22)-5331-5335。
PDM是一種基于對象樣本分析的統(tǒng)計(jì)模型。從Cootes等人開始，這種基于訓(xùn)練樣本上標(biāo)記點(diǎn)的形體建模思想一直被貫徹下來，如文獻(xiàn)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.泰勒，“有限元分析中點(diǎn)分布模型與形狀模型的結(jié)合”，計(jì)算機(jī)圖像與視覺，1995，13(5)403-409。最典型的想法是，使標(biāo)記點(diǎn)自動(dòng)對齊以確定圖形的位置和形變類型。對齊圖像后可以很容易的比較同組標(biāo)記點(diǎn)在不同例圖中的位置變化，只需要計(jì)算它們的坐標(biāo)就可以了。因而，PDM模型對于形狀的描述是通過輪廓點(diǎn)向量來完成的，它由中間形狀和形狀變化樣式組成。通過對這種分布情況建模，我們可以得到與原始訓(xùn)練集相似的新形狀。
在三維空間，如果訓(xùn)練集圖像中的每個(gè)圖形都由一組數(shù)量為n的標(biāo)記點(diǎn)描述，那么它可以用一個(gè)3n維的矢量表示(如果是二維對象即為2n維矢量)。訓(xùn)練集中的標(biāo)記工作非常重要，通常是手動(dòng)完成的。另外，在研究中我們要收集數(shù)百名例柔性體的個(gè)例以得到比較好的統(tǒng)計(jì)結(jié)果，如文獻(xiàn)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.，“手術(shù)導(dǎo)航系統(tǒng)中基于稀疏點(diǎn)數(shù)據(jù)及密集點(diǎn)分布模型的股骨近端精確魯棒重建”，IEEE生物醫(yī)學(xué)工程會(huì)刊，vol.54，no.12，2007(12)2109-2122；文獻(xiàn)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.，“統(tǒng)計(jì)形狀模型訓(xùn)練集人工擴(kuò)大方法”，IEEE醫(yī)學(xué)圖像會(huì)刊，vol.27，2008(11)1643-1654。訓(xùn)練集形成的概率模型維度可能非常的大(在一個(gè)心臟模型中通?？蛇_(dá)到上萬維)，所以建立PDM的過程是非?？菰镉仲M(fèi)時(shí)費(fèi)力的。然而，在一些應(yīng)用中，單獨(dú)的一個(gè)PDM模型或是有限個(gè)PDM模型都不足以很好的描述對象變化。事實(shí)上，不論是連續(xù)PDM模型還是模型插入方法都是為了滿足對象的性質(zhì)，即對象時(shí)間和空間上的變化。例如在心臟建模時(shí)，需要不同時(shí)間段的模型才能很好地描述心室形狀在心臟跳動(dòng)不同階段的變化。但是，在構(gòu)建整個(gè)心臟運(yùn)動(dòng)周期的連續(xù)模型時(shí)，我們只能由已知的模型序列得到特定幾個(gè)時(shí)刻的模型，而其他任意時(shí)刻的模型就不得而知了，盡管在基于模型的圖像分割領(lǐng)域PDM建模出現(xiàn)已有十來年，但對于模型插值及相關(guān)內(nèi)容的深入研究還沒有開展過。


發(fā)明內(nèi)容
為了克服已有的柔性體點(diǎn)分布模型的不能形成實(shí)時(shí)分布的模型、模型精度差、實(shí)用性差的不足，本發(fā)明提供一種能夠形成實(shí)時(shí)分布的中間過渡模型、模型精度高、實(shí)用性良好的基于插值算法的柔性體點(diǎn)分布模型的特征值和特征向量計(jì)算方法。
本發(fā)明解決其技術(shù)問題所采用的技術(shù)方案是 一種基于插值算法的柔性體點(diǎn)分布模型的特征值和特征向量的計(jì)算方法，所述計(jì)算方法包括以下步驟 1)、載入兩個(gè)相鄰時(shí)刻Ta和Tb的柔性體分布模型Ma和Mb，其計(jì)算式為 Ma＝<ya，Φa，λa>，Mb＝<yb，Φb，λb>； 2)、通過所述柔性體分布模型Ma和Mb，分別計(jì)算平均形狀ya和yb、特征矩陣Φa和Φb、以及特征值λa和λb，設(shè)定柔性體從Ta到Tb呈連續(xù)變化， 2.1)中間模型的特征值的計(jì)算 λt≈a2λa+b2λb，(25) 這里λa和λb分別為前后兩個(gè)相鄰PDM模型的特征值向量，即 λa＝(λa1，λa2，λa3，...)，(26) λb＝(λb1，λb2，λb3，...). (27) 2.2)、中間模型的特征向量的估計(jì) 插入中間模型的協(xié)方差矩陣近似為 St≈a2Sa+b2Sb (28) 其中，協(xié)方差矩陣Sa和Sb由特征矩陣Φa和Φb近似得到，Φa和Φb從插補(bǔ)前已知的點(diǎn)分布模型中得到；設(shè)定中間模型建立前先對訓(xùn)練集數(shù)據(jù)進(jìn)行PCA分析以減少參數(shù)個(gè)數(shù)，并選取t種形變模式來描述整個(gè)對象的形變；由于Φ是協(xié)方差矩陣中最前面的t個(gè)特征矢量，則有 其中，Λa和Λb是由特征矢量λa和λb組合形成的對角矩陣；用奇異值分解的方法從式得到插入中間模型的特征矢量 由兩個(gè)已知的相鄰PDM模型及相關(guān)參數(shù)a、b、Λa、Λb、Φa和Φb確定正交矩陣Ut，所述正交矩陣Ut的前t行就是特征矩陣Φt。
進(jìn)一步，所述步驟2)中還包括 2.3)、中間模型的協(xié)方差矩陣由下式推導(dǎo)得出 令 因?yàn)? 顯然有 根據(jù)協(xié)方差矩陣，矩陣St由Sa、Sb和Sab計(jì)算得到，通過如下公式 用奇異值分解法，可令 通常Φab≠Φba，因?yàn)榫仃嘢ab并非對稱矩陣，Φab和Φba都包含t個(gè)矢量； 在插補(bǔ)時(shí)，點(diǎn)分布模型保留特征值和特征向量，中間模型的協(xié)方差矩陣通過下式計(jì)算 再進(jìn)一步，所述步驟2)中還包括 2.3)、特征向量組的計(jì)算 通過式(31)來計(jì)算特征向量時(shí)，將(31)改寫為 其中是一個(gè)r×N維的矩陣，其中r是主成分?jǐn)?shù)目，N是樣本維數(shù)，即標(biāo)定點(diǎn)的數(shù)量與標(biāo)記點(diǎn)維度的乘積； 如果Ui是矩陣G＝LΦa的主要特征向量，對應(yīng)特征值為λi，則(Φaui)和λi分別為矩陣S＝ΦaL的特征向量和特征值。
由Gui＝λiui，，得 S(Φaui)＝ΦaL(Φaui)＝ΦaGui＝Φa(λiui)＝λi(Φaui)(36)。
本發(fā)明的技術(shù)構(gòu)思為PDM方法是一種先驗(yàn)建模方法，它建立在目標(biāo)對象的一組樣本(即訓(xùn)練集)已知的基礎(chǔ)上。這組樣本集可以提供對象的形狀及形變的統(tǒng)計(jì)學(xué)描述。在實(shí)際應(yīng)用中，這組對象樣本的形狀描述通常是使用輪廓來進(jìn)行的(如圖片上的一組像素的坐標(biāo)值)。進(jìn)一步的，可以在輪廓上選取一些標(biāo)記點(diǎn)來描述目標(biāo)物形狀，而這些標(biāo)記點(diǎn)往往與目標(biāo)物的某種特征相關(guān)聯(lián)。
不同樣本上的對應(yīng)標(biāo)記點(diǎn)在選取的時(shí)候遵循著“特征位置大致相同”的原則，但往往這種定位是不夠精確的，存在一定偏差。這種偏差要?dú)w因于不同個(gè)體間的自然形變。我們期望相對于整個(gè)形體的尺度而言，這種偏差是十分微小的。PDM方法允許我們模仿這種微小的差異，并且指出哪些偏差的確是微小的，哪些相對而言更大一些。
由于訓(xùn)練集中的形體是從不同的圖片集中獲得的，而不同圖片集的坐標(biāo)系并不統(tǒng)一，因此有必要先將所有訓(xùn)練集形體進(jìn)行坐標(biāo)歸一化，即將它們對齊。對齊的過程實(shí)際上就是對每一個(gè)樣本模型尋找適當(dāng)?shù)南嗨菩宰儞Q的過程，包括平移、縮放和旋轉(zhuǎn)，以使得不同樣本之間盡可能的接近，另一方面，選定的變換也應(yīng)該令各個(gè)樣本模型與由所有樣本得到的平均模型盡可能接近。為了方便表述，假設(shè)我們的樣本只有兩個(gè)，歸一化簡化為對齊這兩個(gè)樣本。每個(gè)樣本的形狀都由一個(gè)標(biāo)記點(diǎn)坐標(biāo)組成的向量來表示 x＝(x1，y1，z1，...，xn，yn，zn)T， (1) 其中(xi，yi，zi)是某個(gè)標(biāo)記點(diǎn)在三維圖像中的空間坐標(biāo)。如果是二維醫(yī)學(xué)圖像，則標(biāo)記點(diǎn)用(xi，yi)表示。
圖片上標(biāo)記點(diǎn)的獲取有多種方式，人工手動(dòng)標(biāo)記或采用一定得自動(dòng)標(biāo)記算法都是可以的。例如，Izard等人提出的一種MRI圖像標(biāo)記方法，使用統(tǒng)一的算法規(guī)則在不同人腦圖片中進(jìn)行組織結(jié)構(gòu)標(biāo)記，如文獻(xiàn)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.斯塔克，“大腦核磁共振圖像的自動(dòng)標(biāo)記”，SPIE醫(yī)學(xué)圖像網(wǎng)絡(luò)研討會(huì)，2005年2月12-17日，美國，加利福尼亞，圣地亞哥。得到一組經(jīng)過標(biāo)記的訓(xùn)練樣本后，我們需要將它們對齊到一個(gè)統(tǒng)一的坐標(biāo)系下?？梢允褂脧V義對齊算法(GPA)來將訓(xùn)練樣本對齊，并使得所有樣本形狀到平均模型的距離平方和最小，實(shí)際上也就是尋找變換Ti(平移，旋轉(zhuǎn)，縮放)，使下式最小化 Dm＝∑|m-Ti(xi)|2， (2) 其中 且xi是訓(xùn)練集中的一個(gè)形狀樣本，它是一個(gè)3n維的向量表示(在三維空間中有n個(gè)標(biāo)記點(diǎn)的情況下)。
對齊后的訓(xùn)練集組成了一個(gè)三維空間上的點(diǎn)云，可以看作是一個(gè)概率密度函數(shù)。為了降低運(yùn)算成本和內(nèi)存占用，我們采用主成分分析法(PCA)尋找點(diǎn)云中的主元，并以前列最大的少數(shù)分量建立模型。這些被選出的分量可以描述對象主要的形變。以心動(dòng)周期內(nèi)左心室的模型為例，我們經(jīng)過主成分分析后的用60個(gè)特征向量已經(jīng)能足夠描述心室的形狀及其變化了。最后，PDM模型可以表達(dá)為 x＝x+p1b1+…+ptbt＝x+Φc(4) 其中，x是對齊后的平均形狀向量，Φ是一個(gè)3n×t的矩陣，Φ的每一列表示主成分方向上的單位向量，c是由形狀參數(shù)組成的t維列向量?？梢?，經(jīng)過PCA分析后形狀的維度已經(jīng)由3n下降到了t。
經(jīng)過以上步驟我們得到了一個(gè)類似于點(diǎn)分布模型的統(tǒng)計(jì)學(xué)模型，該模型可應(yīng)用于柔性體的物體建模，以及在新的圖像中定位和分割新個(gè)例。在訓(xùn)練集的學(xué)習(xí)過程中我們可以得到形狀參數(shù)b的變化范圍，只要在該范圍內(nèi)變化，任意b的值都可以生成合理的新樣本(仿真形體)。通常bi的變化為λj(矩陣Φ中對應(yīng)第i大的特征值)或者
從樣本集中訓(xùn)練出一個(gè)統(tǒng)計(jì)學(xué)形態(tài)模型后，就可以用PDM來解釋新的圖像個(gè)例了。為了將模型與圖像匹配，通常采用一種主動(dòng)重復(fù)迭代的方法。該方法使模型不斷地變形來適應(yīng)新圖像中目標(biāo)對象的輪廓。模型形狀的變化受到一些統(tǒng)計(jì)數(shù)據(jù)的約束，也即只能在訓(xùn)練樣本集中所得出的范圍內(nèi)進(jìn)行變化。對于形狀模型，我們還要求圖像輪廓大致位于模型點(diǎn)的附近?？梢员硎緸槟滁c(diǎn)的梯度與輪廓上過該點(diǎn)的法線方向的變化。令統(tǒng)計(jì)模型的輪廓沿圖像輪廓不斷移動(dòng)，同時(shí)計(jì)算兩輪廓的馬氏距離最小化就可以最終得到真實(shí)的邊界位置。因此，要在特定實(shí)例中獲得目標(biāo)柔性對象的形狀需要反復(fù)執(zhí)行以下兩步直到收斂(1)沿模型各點(diǎn)的法線尋找圖像輪廓上相匹配的對應(yīng)點(diǎn)(使得模型點(diǎn)與對應(yīng)圖像點(diǎn)的馬氏距離最小)；(2)改變模型位置及形狀參數(shù)使得模型各標(biāo)記點(diǎn)更接近圖像輪廓上找到的對應(yīng)點(diǎn)。也就是說，模型的擬合過程實(shí)際上就是在圖像上搜尋模型點(diǎn)的最相似位置及不斷修正整體幾何變換T及形狀參數(shù)c的過程，使之最小化。數(shù)學(xué)上可以表示為 f＝|X-T(x+Pc)|2＝f(b，Xc，Yc，Zc，s，α，β，γ)(5) ＝|X-T(x+Pc；Xc，Yc，Zc，s，α，β，γ)|2 其中X是在當(dāng)前迭代時(shí)得到的臨時(shí)模型。
這是一個(gè)尋找最小值的問題，可以通過一些非線性優(yōu)化的迭代方法求解。最終我們可以確定整體變換參數(shù)T，以及個(gè)體實(shí)例的形狀參數(shù)b c＝PT(T-1(X)-x). (6) 近幾年，人們已經(jīng)就基于PDM的柔性體建模方法的許多相關(guān)問題進(jìn)行了深入的研究，如形體對齊、自動(dòng)標(biāo)記、平均形體生成、形變建模、模態(tài)分析、圖像分割等等問題，如文獻(xiàn)15蔣曉悅，趙榮椿，一種改進(jìn)的活動(dòng)輪廓圖像分割技術(shù)，中國圖象圖形學(xué)報(bào)A輯，2004，9(9)1019-1024；文獻(xiàn)16馬峰，唐澤圣，夏紹瑋.多尺度幾何活動(dòng)曲線及MR圖像邊界提取.計(jì)算機(jī)學(xué)報(bào).2000.23(8)；文獻(xiàn)17Gilliam，A.D.；Epstein，F(xiàn).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(xiàn).H.；?？松?，S.T.，“基于軌跡場模型的心動(dòng)復(fù)原”，IEEE生物醫(yī)學(xué)信息技術(shù)會(huì)刊，vol.13，2009(3)226-235。此外，在該思想廣泛地應(yīng)用于非剛體對象建模的基礎(chǔ)上，并可進(jìn)一步擴(kuò)展為諸如活動(dòng)形體模型(ASM)和活動(dòng)表現(xiàn)建模(AAM)等以應(yīng)用于更深入的圖像分析。目前，對PDM模型最成功的表達(dá)是通過PCA方法得到的，如文獻(xiàn)18Tobon-Gomez，C.；Butakoff，C.；Aguade，S.；Sukno，F(xiàn).；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(xiàn).；莫拉蓋斯，G.；弗蘭基，A.F.，“基于模擬圖像采集的三維ASM密度模型自動(dòng)重建應(yīng)用于心臟門閘SPECT圖像研究”，IEEE醫(yī)學(xué)圖像會(huì)刊，vol.27，2008(11)1655-1667。PCA方法非常簡單易懂，如文獻(xiàn)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魯棒分析的三個(gè)過程對比”，高維數(shù)據(jù)魯棒性學(xué)術(shù)討論會(huì)，2004，首先把一個(gè)形體投影到經(jīng)過學(xué)習(xí)PCA空間(由線性無關(guān)方向形成的正交空間)，然后從訓(xùn)練樣本中提取少量與形變模式相關(guān)的參數(shù)(稱為主要成分)，再用這些參數(shù)來控制形體的變形。顯然，這種形變是沿著樣本集所確定的最大形變方向進(jìn)行的。由于只研究幾個(gè)主要方向上的形變，無論是計(jì)算復(fù)雜度還是對存儲(chǔ)空間的要求都大大降低了。從數(shù)學(xué)角度講，首先是從訓(xùn)練集中估計(jì)得到平均形狀y和形體協(xié)方差矩陣S 這樣，一個(gè)合理的形狀可以被表達(dá)為 y＝y(tǒng)+Φc (8) 式中向量c由對應(yīng)于各個(gè)特征向量的形變參數(shù)組成。矩陣Φ則包含了一組由協(xié)方差矩陣S分解得到的特征向量，這組向量可線性組合成任意的新形狀y。根據(jù)PCA理論，c中的每個(gè)參數(shù)都控制形體在某一方向上的形變，并且，各個(gè)方向都彼此獨(dú)立，參數(shù)按照形體在對應(yīng)方向上形變由大到小順序排列。本發(fā)明考慮在一個(gè)PDM模型僅包含平均形狀、特征值和特征向量而不包括原始訓(xùn)練集的數(shù)據(jù)的情況下進(jìn)行插值運(yùn)算。即 M＝<y，Φ，λ>(9) 如上所述，單個(gè)PDM模型可以由一組柔性對象的形體獲得，例如，我們研究的醫(yī)院患者的心臟圖片(一些被觀測的柔性體個(gè)例)。然而，有時(shí)我們要觀測這些柔性體隨時(shí)間的變化情況，這樣可以得到一系列的形狀信息，例如心臟在整個(gè)活動(dòng)周期內(nèi)的三維變形。而醫(yī)院現(xiàn)有設(shè)備只能獲得若干時(shí)刻的三維圖像，即對應(yīng)于某些Ti時(shí)刻的模型Mi。
對一個(gè)特定的動(dòng)態(tài)形變物體，假設(shè)我們已經(jīng)在離散的時(shí)間坐標(biāo)(T1，T2，...)上建立起一組PDM模型(M1，M2，...)。我們可以把問題描述為給出已知的兩個(gè)PDM模型，Ta時(shí)刻的模型Ma和Tb時(shí)刻的模型Mb(模型已知意味著我們知道了平均形狀ya和yb、特征矩陣Φa和Φb、以及特征值λa和λb)，如何得到在t(Ta≤t≤Tb)時(shí)刻的中間模型Mt的相關(guān)參數(shù)。這些參數(shù)有平均形狀、特征值和特征向量，即Mt＝<yt，Φt，λt>(圖1)。這里我們假設(shè)在Ta到Tb的時(shí)間段中，一個(gè)形體從Yai連續(xù)變化到Y(jié)bi，二者都是n維的，如Yai＝(xai1，xai2，xai3，...，xain)T。這里還定義了兩個(gè)常數(shù) b＝1-a(11) a和b反映了中間模型相對于Ma和Mb的距離。
(1)中間模型的平均形狀 根據(jù)以上的定義和假設(shè)，用線性插值(可滿足大部分的實(shí)際應(yīng)用)得到Y(jié)ai和Ybi之間的中間模型形狀為 yti＝ayai+bybi(12) a和b由式(10)和(11)定義，即a+b＝1。
插入的平均形狀為 所以，我們得到以下結(jié)論 定理1線性插值時(shí)，插入的PDM中間模型平均形狀是相鄰模型平均形狀的線性組合，其參數(shù)由式(13)確定。
(2)中間模型的形變 相似的，中間PDM模型的形變可以這樣算出 Δyti＝y(tǒng)ti-yt＝(ayai+bybi)-(aya+byb) ＝a(yai-ya)+b(ybi-yb) (14) ＝aΔyai+bΔybi 從上式可以得到以下結(jié)論 定理2線性插值時(shí)，插入的PDM中間模型的形狀變化速度是相鄰模型形變速度的線性組合，其參數(shù)由式(14)確定。
(3)中間模型的協(xié)方差矩陣 中間模型的協(xié)方差矩陣可以由下式推導(dǎo)得出 令 因?yàn)? 顯然有 最終，得到以下推論。
推論1線性插值時(shí)，中間模型的協(xié)方差矩陣可由以下表達(dá)式得到 這里Sab表示模型Ma和Mb的交變矩陣。我們知道協(xié)方差矩陣是通過對同一時(shí)間段內(nèi)不同個(gè)體的采樣數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析得出的，而交變矩陣則描述了在不同采樣時(shí)間得到的兩個(gè)對象模型的形變過程。
中間模型的估算 (1)交變矩陣 為方便我們的討論，假設(shè)模型是n維的，如 Δyai＝[da1i da2i ... dani]T (21) 則，交變矩陣為 通常情況下，樣本是正態(tài)或平均分布的，交變矩陣中的各元素趨于零，即 所以，交變矩陣對角線元素遠(yuǎn)小于矩陣Sa和Sb中的對應(yīng)元素，即 (2)特征值的估算 可以看出，交變矩陣Sab和Sba對于中間模型特征值的計(jì)算影響非常小，因此我們可以估計(jì) λt≈a2λa+b2λb， (25) 這里λa和λb分別為前后兩個(gè)相鄰PDM模型的特征值向量，即 λa＝(λa1，λa2，λa3，...)， (26) λb＝(λb1，λb2，λb3，...).(27) 對λt中元素估計(jì)的精度主要取決于樣本集的數(shù)量和分布情況。經(jīng)過計(jì)算機(jī)模擬發(fā)現(xiàn)，當(dāng)樣本集為30個(gè)隨機(jī)均勻分布時(shí)，誤差率可維持在百分之五以下。
(3)特征向量的估計(jì)方法 當(dāng)樣本集很大時(shí)，插入中間模型的協(xié)方差矩陣可以近似為 St≈a2Sa+b2Sb(28) 所以，協(xié)方差矩陣Sa和Sb可以由特征矩陣Φa和Φb近似得到，而Φa和Φb是可以從插補(bǔ)前已知的點(diǎn)分布模型中得到。假設(shè)PDM模型建立前先對訓(xùn)練集數(shù)據(jù)進(jìn)行PCA分析以減少參數(shù)個(gè)數(shù)，并選取t種形變模式來描述整個(gè)對象的形變。由于Φ是協(xié)方差矩陣中最前面的t個(gè)特征向量，我們知道 其中Λa和Λb是由特征向量λa和λb組合形成的對角矩陣。然后，我們可以用奇異值分解(SVD)的方法從式得到插入中間模型的特征向量 現(xiàn)在，已由兩個(gè)已知的相鄰PDM模型及相關(guān)參數(shù)a、b、Λa、Λb、Φa和Φb，確定了正交矩陣Ut，且Ut的前t行就是特征矩陣Φt。
(4)特征向量的精確計(jì)算方法 通過上面的公式，我們可以由兩個(gè)已知PDM模型計(jì)算得到插入中間模型的各個(gè)參數(shù)(如平均形狀、特征值、特征向量等)。然而，由于我們忽略了交變矩陣的作用，因此這種方法存在一定的誤差，尤其是當(dāng)樣本集容量很小或其分布比較異常時(shí)這種誤差就更明顯。事實(shí)上我們有一種精確的方法，就是直接生成訓(xùn)練集中所有相鄰兩個(gè)PDM模型的交變矩陣。根據(jù)式(20)，矩陣St可由Sa、Sb和Sab計(jì)算得到，這意味著矩陣Sab需要在構(gòu)建Ma和Mb時(shí)就得到。實(shí)際上，可以通過如下公式 用奇異值分解法，可令 通常Φab≠Φba，這是因?yàn)榫仃嘢ab并不是一個(gè)對稱矩陣。Φab和Φba都包含t個(gè)矢量。
在插補(bǔ)時(shí)，由于點(diǎn)分布模型只保留了特征值和特征向量的信息，中間模型的協(xié)方差矩陣應(yīng)通過下式計(jì)算 相似的，中間模型的特征向量可以由式(31)得到。這種方法由于不受訓(xùn)練樣本分布的影響，因此是一種更加準(zhǔn)確的解決方案。
(5)特征向量組的計(jì)算 當(dāng)直接通過式(31)來計(jì)算特征向量時(shí)，我們會(huì)發(fā)現(xiàn)由于矩陣通常很大，計(jì)算過程對存儲(chǔ)空間和運(yùn)算效率都有很高的要求。例如，我們實(shí)驗(yàn)中有一個(gè)三維心臟模型，包含2848個(gè)標(biāo)記點(diǎn)，通過主成分分析后選擇了90種形變模式，那么該模型的一個(gè)協(xié)方差矩陣將占據(jù)584MB的內(nèi)存空間，連最簡單的加減操作可能也要用數(shù)分鐘時(shí)間。因此，為了簡化計(jì)算復(fù)雜度，我們將(31)改寫為 其中是一個(gè)r×N維的矩陣，其中r是主成分?jǐn)?shù)目，N是樣本維數(shù)，即標(biāo)定點(diǎn)的數(shù)量與標(biāo)記點(diǎn)維度的乘積。下面我們來描述簡化后的計(jì)算方法。
定理3如果Ui是矩陣G＝LΦa的主要特征向量，對應(yīng)特征值為λi，則(Φaui)和λi分別為矩陣S＝ΦaL的特征向量和特征值。
其證明相當(dāng)簡單，由Gui＝λiui，，得 S(Φaui)＝ΦaL(Φaui)＝ΦaGui＝Φa(λiui)＝λi(Φaui) (36) 本發(fā)明的有益效果主要表現(xiàn)在能夠形成實(shí)時(shí)分布的中間過渡模型、模型精度高、實(shí)用性良好。



圖1是從時(shí)間序列中的臨近模型插值到得的瞬時(shí)PDM模型的示意圖。
圖2是典型的心動(dòng)周期分為六個(gè)階段的示意圖。
圖3是心臟不同階段特征的六個(gè)PDM模型的示意圖。
圖4是插值生成新的PDM模型以進(jìn)行柔性體分析的示意圖。
圖5是通過插值生成的PDM心臟模型的時(shí)間平滑變化的示意圖。

具體實(shí)施例方式 下面結(jié)合附圖對本發(fā)明作進(jìn)一步描述。
參照圖1～圖5，一種基于插值算法的柔性體點(diǎn)分布模型的特征值和特征向量的計(jì)算方法，其特征在于所述計(jì)算方法包括以下步驟 1)、載入兩個(gè)相鄰時(shí)刻Ta和Tb的柔性體分布模型Ma和Mb，其計(jì)算式為 Ma＝<ya，Φa，λa>，Mb＝<yb，Φb，λb>； 2)、通過所述柔性體分布模型Ma和Mb，分別計(jì)算平均形狀ya和yb、特征矩陣Φa和Φb、以及特征值λa和λb，設(shè)定柔性體從Ta到Tb呈連續(xù)變化， 2.1)中間模型的特征值的計(jì)算 λt≈a2λa+b2λb，(25) 這里λa和λb分別為前后兩個(gè)相鄰PDM模型的特征值向量，即 λa＝(λa1，λa2，λa3，...)，(26) λb＝(λb1，λb2，λb3，...). (27) 2.2)、中間模型的特征向量的估計(jì) 插入中間模型的協(xié)方差矩陣近似為 St≈a2Sa+b2Sb (28) 其中，協(xié)方差矩陣Sa和Sb由特征矩陣Φa和Φb近似得到，Φa和Φb從插補(bǔ)前已知的點(diǎn)分布模型中得到；設(shè)定中間模型建立前先對訓(xùn)練集數(shù)據(jù)進(jìn)行PCA分析以減少參數(shù)個(gè)數(shù)，并選取t種形變模式來描述整個(gè)對象的形變；由于Φ是協(xié)方差矩陣中最前面的t個(gè)特征矢量，則有 其中，Λa和Λb是由特征矢量λa和λb組合形成的對角矩陣；用奇異值分解的方法從式得到插入中間模型的特征矢量 由兩個(gè)已知的相鄰PDM模型及相關(guān)參數(shù)a、b、Λa、Λb、Φa和Φb確定正交矩陣Ut，所述正交矩陣Ut的前t行就是特征矩陣Φt； 2.3)、中間模型的協(xié)方差矩陣由下式推導(dǎo)得出 令 因?yàn)? 顯然有 根據(jù)協(xié)方差矩陣，矩陣St由Sa、Sb和Sab計(jì)算得到，通過如下公式 用奇異值分解法，可令 通常Φab≠Φba，因?yàn)榫仃嘢ab并非對稱矩陣，Φab和Φba都包含t個(gè)矢量； 在插補(bǔ)時(shí)，點(diǎn)分布模型保留特征值和特征向量，中間模型的協(xié)方差矩陣通過下式計(jì)算 2.3)、特征向量組的計(jì)算 通過式(31)來計(jì)算特征向量時(shí)，將(31)改寫為 其中是一個(gè)r×N維的矩陣，其中r是主成分?jǐn)?shù)目，N是樣本維數(shù)，即標(biāo)定點(diǎn)的數(shù)量與標(biāo)記點(diǎn)維度的乘積； 如果Ui是矩陣G＝LΦa的主要特征向量，對應(yīng)特征值為λi，則(Φaui)和λi分別為矩陣S＝ΦaL的特征向量和特征值。
由Gui＝λiui，，得 S(Φaui)＝ΦaL(Φaui)＝ΦaGui＝Φa(λiui)＝λi(Φaui) (36)。
在我們一個(gè)關(guān)于心臟圖像分析的研究項(xiàng)目中，可以根據(jù)心動(dòng)周期定義一個(gè)規(guī)范化的時(shí)間坐標(biāo)軸。例如，一個(gè)心跳周期可以分為六個(gè)生理階段心房收縮、心室收縮、心室射血、心室舒張、心室快速填充及舒張后期休息(圖2)。也有人進(jìn)一步將心室射血的過程分為射血前期和射血后期。左心與右心的運(yùn)動(dòng)過程相似但不完全相同。主要的差別是，左心室與動(dòng)脈的收縮壓比右心的大三至四倍，而且左右心泵血過程的時(shí)間方面也有顯著區(qū)別，如文獻(xiàn)20楊柳，汪天富，林江莉等，旋轉(zhuǎn)掃描超聲心臟圖像的中軸匹配與圖像插值，生物醫(yī)學(xué)工程學(xué)雜志，2004，21(1)28-33，41；文獻(xiàn)21陳強(qiáng)，周則明，王平安，夏德深，帶標(biāo)記線左心室MR圖像的自動(dòng)分割，中國圖象圖形學(xué)報(bào)A輯，2004，9(6)666-673；文獻(xiàn)22羅渝蘭，鄭昌瓊等，基于Snake模型的B超心臟圖像的腔體分割及算法研究，四川師范大學(xué)學(xué)報(bào)自然科學(xué)版，2002，25(1)66-69。由于整個(gè)心臟的空間形態(tài)在一個(gè)心動(dòng)周期內(nèi)變化很大，而單一的模型很難準(zhǔn)確描述動(dòng)態(tài)的心臟運(yùn)動(dòng)過程，所以我們需要建立一組模型來描述心動(dòng)不同階段的心臟形態(tài)。
在gated-SPECT的圖像研究中，可參考心電圖信號(hào)進(jìn)行時(shí)間規(guī)范化。我們希望通過對心臟形狀進(jìn)行統(tǒng)計(jì)后可以得到全部的六個(gè)模型，并且這些模型能用來準(zhǔn)確地描述心臟的特定狀態(tài)，如心臟舒張末期和心臟收縮末期。圖3在標(biāo)準(zhǔn)心電圖曲線中標(biāo)出了各個(gè)階段模型出現(xiàn)的位置。圖中我們已經(jīng)把壓力值幅度經(jīng)過了修正，使得波峰處的壓力值為正常人體的收縮壓(120mm汞柱)，所以圖中也表示一條典型的左心室壓力曲線。
如果在正常的心動(dòng)周期中有6個(gè)PDM模型，則式(6)中的Ti(代表某一相位所在的時(shí)刻)可以定義為
式中t是獲得gated-SPECT圖像的時(shí)間值，Hc是心率的的倒數(shù)，

表示不大于某小數(shù)的最大整數(shù)。
然而，在實(shí)際的例子中，常會(huì)產(chǎn)生不同數(shù)目的相位(如圖4中有8個(gè))，這是由于心臟圖像的采樣通常有固定的頻率，如每100ms(圖4)。那么就有必要通過插值得到采樣時(shí)刻的PDM模型，以期得到更好的三維心臟形狀的分割。
柔性體點(diǎn)分布模型插值算法，作者研究過程中分別采用計(jì)算機(jī)仿真及真實(shí)心臟數(shù)據(jù)進(jìn)行了大量的實(shí)驗(yàn)。在仿真中，我們可以假設(shè)一組九維(或更高)的模型，yi＝[x1i，x2i，...，x9i]T。兩個(gè)PDM模型M1和M2的數(shù)據(jù)集都是隨機(jī)產(chǎn)生的。各組數(shù)據(jù)都含有100個(gè)樣本，作為生成點(diǎn)分布模型的訓(xùn)練集，例如 為生成一個(gè)插入模型，我們從S1計(jì)算M1的特征值及特征向量，對S2也同樣。假設(shè)插入模型的左右距離為a＝0.3，b＝0.7。為了進(jìn)行對比，我們還計(jì)算它們的交變矩陣Sab和相應(yīng)的特征向量組。實(shí)驗(yàn)結(jié)果顯示，本中所述的方法能很好的計(jì)算出插入模型的各種參數(shù)。由式(31)計(jì)算出的特征值誤差小于1.83％，由式(35)計(jì)算出的特征值誤差小于0.05％。該仿真實(shí)驗(yàn)的C++程序已由所在實(shí)驗(yàn)室開發(fā)完成。
用真實(shí)的心臟數(shù)據(jù)進(jìn)行了測試。通過灌注SPECT成像法得到的大量圖像數(shù)據(jù)中，我們預(yù)先建立了左心室在不同時(shí)刻的五個(gè)PDM模型。每個(gè)模型由246組病例的真實(shí)圖像的訓(xùn)練得到(圖像主要由巴塞羅那醫(yī)院提供)，包括一個(gè)平均形狀的VTK表達(dá)，一組特征值矢量和特征向量。所有醫(yī)學(xué)圖像的處理過程和模型的插入運(yùn)算也由都C++程序?qū)崿F(xiàn)(研究組共開發(fā)了約600MB程序代碼)。實(shí)驗(yàn)時(shí)，我們在每兩個(gè)相鄰的PDM之間均勻地插入四個(gè)中間模型，得到的實(shí)驗(yàn)結(jié)果令人非常滿意。圖5顯示通過插值生成的PDM模型隨時(shí)間能夠平滑地變化左心室柔性模型。
點(diǎn)分布模型作為一種描述柔性對象的統(tǒng)計(jì)學(xué)模型，已被證明在很多應(yīng)用上是非常有效的，尤其是在醫(yī)學(xué)圖像分析領(lǐng)域。為了分析柔性體的動(dòng)態(tài)變化情況，本發(fā)明首創(chuàng)研究了點(diǎn)分布模型的插值算法。研究結(jié)果表明，可由PDM的平均形狀、特征值和特征向量等參數(shù)生成任意時(shí)刻的中間模型。事實(shí)上，插入模型的平均形狀及形變就是它相鄰前后模型對應(yīng)項(xiàng)的線性組合，即yt＝aya+byb，Δyti＝aΔyai+bΔybi。其特征值則可以采用式λt≈a2λa+b2λb計(jì)算得到。最后，為求出插入模型的協(xié)方差矩陣和特征向量，本發(fā)明公開了交變矩陣形式及其相應(yīng)的計(jì)算方法。
權(quán)利要求
1、一種基于插值算法的柔性體點(diǎn)分布模型的特征值和特征向量的計(jì)算方法，其特征在于所述計(jì)算方法包括以下步驟
1)、載入兩個(gè)相鄰時(shí)刻Ta和Tb的柔性體分布模型Ma和Mb，其計(jì)算式為
Ma＝<ya，Φa，λa>，Mb＝<yb，Φb，λb>；
2)、所述柔性體分布模型Ma和Mb，包括計(jì)算平均形狀ya和yb、特征矩陣Φa和Φb、以及特征值λa和λb，設(shè)定柔性體從Ta到Tb呈連續(xù)變化，則
2.1)中間模型的特征值的計(jì)算
λt≈a2λa+b2λb，(25)
這里λa和λb分別為前后兩個(gè)相鄰PDM模型的特征值向量，即
λa＝(λa1，λa2，λa3，...)，(26)
λb＝(λb1，λb2，λb3，...). (27)
2.2)、中間模型的特征向量的估算
插入中間模型的協(xié)方差矩陣近似為
St≈a2Sa+b2Sb (28)
其中，協(xié)方差矩陣Sa和Sb由特征矩陣Φa和Φb近似得到，Φa和Φb從插補(bǔ)前已知的點(diǎn)分布模型中得到；設(shè)定中間模型建立前先對訓(xùn)練集數(shù)據(jù)進(jìn)行PCA分析以減少參數(shù)個(gè)數(shù)，并選取t種形變模式來描述整個(gè)對象的形變；由于Φ是協(xié)方差矩陣中最前面的t個(gè)特征矢量，則有
其中，Λa和Λb是由特征矢量λa和λb組合形成的對角矩陣；用奇異值分解的方法從式得到插入中間模型的特征矢量
由兩個(gè)已知的相鄰PDM模型及相關(guān)參數(shù)a、b、Λa、Λb、Φa和Φb確定正交矩陣Ut，所述正交矩陣Ut的前t行就是特征矩陣Φt。
2、如權(quán)利要求1所述的基于時(shí)變插值算法的柔性體點(diǎn)分布模型的特征值和特征向量的計(jì)算方法，其特征在于所述步驟2)中還包括
2.3)、中間模型的協(xié)方差矩陣由下式推導(dǎo)得出
令
因?yàn)?br> 顯然有
根據(jù)協(xié)方差矩陣，矩陣St由Sa、Sb和Sab計(jì)算得到，通過如下公式
用奇異值分解法，可令
通常Φab≠Φba，因?yàn)榫仃嘢ab并非對稱矩陣，Φab和Φba都包含t個(gè)矢量；
在插值計(jì)算時(shí)，點(diǎn)分布模型保留特征值和特征向量，中間模型的協(xié)方差矩陣通過下式計(jì)算
3、如權(quán)利要求1或2所述的基于插值算法的柔性體點(diǎn)分布模型的特征值和特征向量的計(jì)算方法，其特征在于所述步驟2)中還包括
2.3)、特征向量組的計(jì)算
通過式(31)來計(jì)算特征向量時(shí)，將(31)改寫為
其中是一個(gè)r×N維的矩陣，其中r是主成分?jǐn)?shù)目，N是樣本維數(shù)，即標(biāo)定點(diǎn)的數(shù)量與標(biāo)記點(diǎn)維度的乘積；
如果Ui是矩陣G＝LΦa的主要特征向量，對應(yīng)特征值為λi，則(Φaui)和λi分別為矩陣S＝ΦaL的特征向量和特征值。
由Gui＝λiui，，得
S(Φaui)＝ΦaL(Φaui)＝ΦaGui＝Φa(λiui)＝λi(Φaui)(36)。
全文摘要
一種基于插值算法的柔性體點(diǎn)分布模型的特征值和特征向量計(jì)算方法，包括以下步驟1)載入兩個(gè)相鄰時(shí)刻Ta和Tb的柔性體分布模型Ma和Mb，其中Ma＝<ya，Φa，λa>，Mb＝<yb，Φb，λb>；2)通過所述柔性體分布模型Ma和Mb，分別計(jì)算平均形狀ya和yb、特征矩陣Φa和Φb、以及特征值λa和λb，設(shè)定柔性體從Ta到Tb呈連續(xù)變化；2.1)中間模型的特征值的計(jì)算；2.2)中間模型的特征向量的估算。本發(fā)明能夠形成實(shí)時(shí)分布的中間過渡模型、模型精度高、實(shí)用性良好。
文檔編號(hào)G06T17/00GK101639948SQ20091010218
公開日2010年2月3日 申請日期2009年8月20日 優(yōu)先權(quán)日2009年8月20日
發(fā)明者陳勝勇, 杜雅慧, 秋 管, 毛國紅, 王萬良 申請人:浙江工業(yè)大學(xué)