專利名稱：能夠減少物體運動造成的成像偽影的生物醫(yī)學(xué)圖像配準(zhǔn)方法和計算機(jī)程序產(chǎn)品的制作方法
技術(shù)領(lǐng)域：
本發(fā)明涉及圖像配準(zhǔn)的方法，該方法由以下步驟組成 a)提供同一物體的至少第一幅和第二幅數(shù)字或數(shù)字化圖像，或者至少第一組和第二組該物體的截面圖像，這些圖像由像素或體素的二維或三維陣列形成； b)在第一幅圖像或第一組圖像中，通過選取一定數(shù)量的像素或體素，將這些像素或體素設(shè)為標(biāo)志或特征，來定義一定數(shù)量的標(biāo)志，即所謂的特征，并生成要跟蹤的所述特征的列表； c)通過對選定為特征的每個像素或體素確定從第一幅到第二幅圖像或從第一組到第二組圖像的光流矢量，從第一幅到第二幅圖像或從第一組到第二組圖像，跟蹤選定為特征的每個像素或體素的位置； d)通過對第二幅或第二組圖像的像素或體素施加逆光流，來對第一幅和第二幅圖像或第一組和第二組圖像進(jìn)行配準(zhǔn)；
背景技術(shù)：
這種為時間序列的部分的圖像的配準(zhǔn)方法可以從例如文獻(xiàn)US6,553,152中獲悉。
文獻(xiàn)US6,553,152提供了用于兩幅圖像的配準(zhǔn)方法，該方法是通過由操作員對標(biāo)志的目視識別以及對兩幅圖像上的對應(yīng)像素進(jìn)行標(biāo)記來完成的。
這種方法使得標(biāo)志的選擇完全依賴操作者的技能，并且當(dāng)圖像中不存在顯著的或意義明確的可識別結(jié)構(gòu)時，該方法在實踐中非常難于實現(xiàn)。考慮核磁共振成像(MRI)成像或超聲圖像或X射線圖像之類的診斷圖像可知，根據(jù)這些圖像的結(jié)構(gòu)特征進(jìn)行對于標(biāo)志的目視識別是極端困難的，并且可能造成重大的錯誤。
許多不同的算法試圖識別輪廓，如Fischer，B.；Modersitzki，J.Curvature Based Registration with Applications to MRMammography；Springer-Verlag Berlin HeidelbergICCS 2002，LNCS 2331，2002；pp202-206所公開的算法，或者試圖識別被測量過兩次或多次或者甚至用多種醫(yī)療器械進(jìn)行過研究的器官或部位的放射線圖像內(nèi)的標(biāo)志，例如Shi J，Tomasi C.Good Features to Track.1994 IEEE Conference onComputer Vision and Pattern Recognition(CVPR’94)1994；593-600、Tommasini T，F(xiàn)usiello A，Trucco E，Roberto V.Making Good FeaturesTrack Better Proc.IEEE Int.Conf.on Computer Vision and PatternRecognition(CVPR’98)1998；178-183、以及Ruiter，N.V.；Muller，T.O.；Stotzka，R.；Kaiser，W.A.Elastic registration of x-ray mammograms andthree-dimensional magnetic resonance imaging data.J Digit Imaging 2001，14，52-55中所公開的算法。
對在不同時刻拍攝的同一物體的平面圖像進(jìn)行配準(zhǔn)以計算物體隨時間的運動，這種配準(zhǔn)稱為二維配準(zhǔn)。實現(xiàn)相同功能但在例如MRI或計算斷層(CT)掃描儀的截面圖像組的三維圖像組上作用的算法，被稱為三維配準(zhǔn)。在這些圖像組內(nèi)，標(biāo)志或物體的運動可能在任何方向上發(fā)生。
根據(jù)要比較的不同成像方式的數(shù)量，配準(zhǔn)算法可以分為單模式配準(zhǔn)算法和多模式配準(zhǔn)算法。
比較兩個MRI序列就是一種單模式配準(zhǔn)。
通常，單模式配準(zhǔn)算法設(shè)法識別兩幅圖像內(nèi)的一定標(biāo)志或復(fù)雜形狀，大多數(shù)情況下處理的是二維圖像。在第二步則比較兩幅圖像內(nèi)對應(yīng)標(biāo)志的位置，并計算其運動矢量。這個信息用來將第二幅圖像的像素移回到新的位置，以消除運動偽影。
“剛性配準(zhǔn)”將二維或三維的圖像/體積作為單個單元整體移動到一定方向?！胺莿傂耘錅?zhǔn)”則僅僅將一定面積/體積的單個或多個像素/體素移往不同的方向。
考慮很柔軟的乳房組織這一例子可知，需要非剛性配準(zhǔn)來消除運動偽影。
由于乳房不具有任何像骨頭一樣的成形結(jié)構(gòu)，這樣特殊的一致性使得它的直接相鄰的部分能夠在不同的方向上移動。因此，顯而易見的是，任何配準(zhǔn)算法或任何其他現(xiàn)有的方法都必須在整個乳房組織各處識別和散布其標(biāo)志。
有一組算法定義了必須尋找的一定數(shù)量的標(biāo)志。它們將按其效價(valence)分類。對于乳房MR圖像而言，由于MRI中脂肪/空氣之間的較高的對比度，這些方法常常會在乳房的外邊界處和胸腔結(jié)構(gòu)處找到許多標(biāo)志(Huwer，S.；Rahmel，J.；Wangenheim，A.v.Data-DrivenRegistration for Lacal Deformations.Pattern Recognition Letters 1996，17，951-957)。
使用壓縮板或者至少通過俯臥壓在雙乳上，通常得到皮膚的所有部分的良好固定，并使其與壓縮板直接接觸(圖1)。相反地，雙乳的內(nèi)部區(qū)域因心跳和呼吸而繼續(xù)輕微移動。由于與外部皮膚/空氣邊界比較而言脂肪與軟組織間的對比度低得多，因此這些算法中乳房中心內(nèi)的有效標(biāo)志數(shù)保持較少。
標(biāo)志選取的一般問題存在于如下的事實中，即選取過少的標(biāo)志可能導(dǎo)致不充分或不精確的配準(zhǔn)。但是，選取額外的標(biāo)志并不一定能保證精確的配準(zhǔn)，這是因為觀察人員或者任何程序可能不得不包括具有較低的重分配確定性(certainty of being reallocated)的標(biāo)志，舉例而言，可以由乳房中心的不同組織類型的低對比度導(dǎo)致該較低的重分配確定性。增加標(biāo)志數(shù)總是會顯著增大配準(zhǔn)的計算復(fù)雜度。
一種特別適用于所謂的非剛性配準(zhǔn)的特征選擇和特征跟蹤算法是所謂的Lucas&Kanade算法及其特定的金字塔式實現(xiàn)，該算法及其實現(xiàn)在Lucas BD，Kanade T.An Iterative Image Registration Techniquewith an Application to Stereo Vision，IJCAI81；674-679 1981和BouguetJY.Pyramidal Implementation of the Lucas Kanade Feature Tracker，TechReport as Part of Open Computer Vision Library Documentation，Intel Corporation enclosed hereinafter as appendix 1中進(jìn)行了詳細(xì)的披露。Lucas&Kanade技術(shù)及其金字塔式特征跟蹤實現(xiàn)特別適合用于自動識別圖像中的可靠標(biāo)志，以及適用于在不同時刻拍攝的圖像之間跟蹤這些標(biāo)志。選定特征的偏移或位移就確定為位移或所謂的光流矢量。這種技術(shù)得到了廣泛的應(yīng)用，例如，用于機(jī)器人學(xué)中的計算機(jī)視覺，并且成為一些基本課程中所教的數(shù)字圖像處理和計算機(jī)視覺領(lǐng)域內(nèi)的技術(shù)人員的常識的部分。
但是，當(dāng)考慮將Lucas & Kanade技術(shù)應(yīng)用到數(shù)字的或數(shù)字化的診斷圖像中，更特別的是應(yīng)用到乳房(診斷)這樣的實際用途時，所述的技術(shù)會識別太多特征，當(dāng)然也包括位于興趣區(qū)之外的太多的特征，并且在譬如包圍乳房的空氣區(qū)易受噪聲的干擾。
從下面的例子中可以更詳細(xì)地看到上述問題。
對比度增強(qiáng)(CE)核磁共振乳腺成像(MRM)是一種有益且重要的研究，對于侵略性腫瘤具有極高的靈敏度。它表現(xiàn)出與乳房軟組織無關(guān)的較高的負(fù)預(yù)測值。近來的多中心研究表明，靈敏度和特異性隨圖像分析技術(shù)的不同而不同。利用1.5T掃描儀，分別計算出96％、93％或86％的靈敏度及30％、50％或70％的特異性。因此，MRM的主要缺點仍是較高的假陽性(false positive)診斷百分率。
通常，MR成像在1.0或1.5特斯拉成像儀上完成，并且使用了廠家的雙乳線圈(在美國則經(jīng)常使用單乳線圈設(shè)備)，而患者處于俯臥姿勢。該治療方案至少包含任意方向上的動態(tài)對比度增強(qiáng)T1加權(quán)序列，且采樣時間在30s到180s之間。在快速濃注釓噴酸二葡甲胺鹽或任何其他的順磁性MR造影劑之前，至少進(jìn)行一次測量，注射之后則要進(jìn)行多次測量。
為了計算乳房組織的一個區(qū)域的實際造影劑攝取量，必須將施加造影劑之前的每個體素的信號強(qiáng)度從施加了造影劑之后的信號強(qiáng)度中減去。由于呼吸和心跳造成雙乳的最小限度的運動以及MR圖像的切片厚度帶來的部分體積效應(yīng)，不同時刻拍攝的同一圖像位置的體素并不確切示出同一塊乳房組織。因此，減影圖像并不完全漆黑一片(除腫瘤外)。微小運動效應(yīng)可以在乳房的外邊界處最清楚地展現(xiàn)出來。由于脂肪組織信號強(qiáng)度較高，而周圍空氣的信號強(qiáng)度約為零，非常小的運動將脂肪置于空氣的之前圖像位置處。從脂肪的高信號強(qiáng)度中減去空氣的零信號強(qiáng)度，不恰當(dāng)?shù)啬７铝薈M應(yīng)用所引起的高信號增量。結(jié)果，減影圖像在成像乳房的至少部分處示出表示造影劑大攝取量的粗白線(圖2乳房的白色外邊界利用符號進(jìn)行標(biāo)記)。
考慮笛卡爾系統(tǒng)定義的三維或體積圖像采集，除了沿x和y方向上的任意運動之外，人們總能發(fā)現(xiàn)沿z方向的某種運動(朝胸部方向)。這是由研究過程中患者的松弛造成的。
為了減少偽影，通常借助某種壓縮方法(因廠家而異)將乳房固定在乳房線圈內(nèi)。但是，盡管這樣，最小限度的運動總會發(fā)生。因此，縱然對乳房加以固定，其圖像在x、y和z方向總會顯現(xiàn)輕微的運動偽影，它們在減影圖像中看來是明亮的像素。如果沒有運動偽影，那么減影圖像除了腫瘤組織所在的區(qū)域外應(yīng)該完全呈黑色。
在乳房外周圍空氣的噪聲信號中可檢測到一定數(shù)量的標(biāo)志。除了必須對考慮的每個特征進(jìn)行跟蹤并帶來不必要的計算負(fù)擔(dān)這一事實外，由于噪聲是個隨機(jī)現(xiàn)象，噪聲中發(fā)現(xiàn)的所有特征或標(biāo)志在第二幅或隨后的圖像中將沒有對應(yīng)的標(biāo)志。但是事實上，在第二幅或隨后的圖像或圖像組中搜索與之前在噪聲中發(fā)現(xiàn)的標(biāo)志相對應(yīng)的標(biāo)志的算法，將通過某種方式在第二幅或第二組圖像中重分配一定量的標(biāo)志。這就會導(dǎo)致錯誤的圖像配準(zhǔn)，并對結(jié)果產(chǎn)生影響。


發(fā)明內(nèi)容
本發(fā)明的目的在于提供一種配準(zhǔn)圖像尤其是三維結(jié)構(gòu)的生物醫(yī)學(xué)診斷圖像的方法，它通過跟蹤三維運動的能力，通過減少數(shù)字或數(shù)字化的圖像中興趣區(qū)之外的噪聲所引起的偽影，以及通過與相鄰組織類型的高或低信號強(qiáng)度無關(guān)地在整個乳房區(qū)各處擴(kuò)散標(biāo)志，來有效克服已知配準(zhǔn)方法的缺陷。
依照本發(fā)明，圖像配準(zhǔn)方法由以下步驟組成 a)提供通過對于組織或組織區(qū)或解剖區(qū)進(jìn)行MRI、超聲或X射線掃描而獲取的第一幅數(shù)字或數(shù)字化的圖像或第一組截面圖像；所述圖像由像素或體素的二維或三維陣列構(gòu)成；提供通過在第二時刻對相同組織或組織區(qū)或解剖區(qū)進(jìn)行MRI、超聲或X射線掃描而獲取的相同解剖區(qū)的至少第二幅數(shù)字或數(shù)字化的圖像或第二組截面圖像，在獲取圖像時可選的在所述組織或組織區(qū)或所述解剖區(qū)存在對比媒質(zhì)； b)在第一幅圖像或第一組圖像中通過選取一定數(shù)量的像素或體素，將該像素或體素設(shè)為標(biāo)志或特征，來定義一定數(shù)量的標(biāo)志，即所謂的特征，并生成要跟蹤的所述特征的列表； c)通過對每個選定為特征的像素或體素確定從第一幅到第二幅圖像或從第一組到第二組圖像的光流矢量，從第一幅到第二幅圖像或從第一組到第二組圖像跟蹤選定為特征的每個像素或體素的位置； d)通過對第二幅或第二組圖像的像素或體素施加逆光流，來對第一幅和第二幅圖像或第一組和第二組圖像進(jìn)行配準(zhǔn)； e)在步驟d中配準(zhǔn)之后，從第二幅圖像的第二圖像數(shù)據(jù)減去第一幅圖像的圖像數(shù)據(jù)； f)顯示通過從第二幅圖像的第二圖像數(shù)據(jù)減去第一幅圖像的圖像數(shù)據(jù)所得到的圖像數(shù)據(jù)。
步驟b“特征選擇”可以如下更詳細(xì)地描述 B1)在第一幅圖像或第一組截面圖像的每個像素或體素的周圍定義像素或體素鄰域，所述的像素或體素鄰域包括有限數(shù)量的像素或體素；對于單幅圖像，選用二維鄰域；對于截面圖像組，則選用三維鄰域； B2)對于每個像素或體素，確定所述像素或體素鄰域的所有像素或體素的平均信號強(qiáng)度值； B3)定義平均信號強(qiáng)度值閾值； B4)比較步驟B2中確定的每個像素或體素鄰域的平均信號強(qiáng)度，并將所述平均信號強(qiáng)度值與平均信號強(qiáng)度值閾值比較； B5)如果步驟B4中所述鄰域的平均信號強(qiáng)度值大于閾值，則將所述像素或體素定義為要跟蹤的特征，并將其添加到要跟蹤的特征的列表中。
所述平均信號強(qiáng)度閾值是可調(diào)的，可以發(fā)生變化。經(jīng)驗的或?qū)嶒灥臉?biāo)準(zhǔn)可以用來確定這種閾值，例如，在一組測試或樣本圖像上執(zhí)行本方法并比較取不同閾值時獲得的結(jié)果，可以確定該閾值。
盡管按照本發(fā)明的方法，可以應(yīng)用任何用于選擇可被定義成要跟蹤的圖像特征的像素或體素的方法，使用根據(jù)著名的Lucas & Kanade方法的自動特征選擇方法卻能得到最好的結(jié)果，Lucas & Kanade方法在Lucas BD，Kanade T.An Iterative Image Registration Technique withan Application to Stereo Vision，IJCAI81；674-679 1981中進(jìn)行了更為詳細(xì)的描述。
利用位移或偏移矢量，即所謂的光流，來跟蹤特征從第一幅圖像到第二幅圖像的偏移量，也是已知的方法。在這種情況下，本發(fā)明也可應(yīng)用到利用了標(biāo)志或特征的一些跟蹤方法中。尤其是如果考慮了診斷成像期間受限乳房的微小移動和Lucas & Kanade特征選擇方法的使用這種特殊情況，本發(fā)明優(yōu)選地采用所謂的“Lucas & Kanade特征跟蹤器的金字塔式實現(xiàn)”(見附錄1)更好。這種實現(xiàn)為普通技術(shù)人員所熟悉，并在一些大學(xué)課程中被講授。Lucas & Kanade特征跟蹤器的這種金字塔式實現(xiàn)的精確而詳細(xì)的描述在Bouguet J-Y.PyramidalImplementation of the Lucas Kanade Feature Tracker，Tech Report as Partof Open Computer Vision Library Documentation，Intel Corporation、以及Shi J，Tomasi C.Good Features to Track.1994 IEEE Conference onComputer Vision and Pattern Recognition(CVPR’94)1994；593-600、以及Tommasini T，F(xiàn)usiello A，Trucco E，Roberto V.Making Good FeaturesTrack Better.Proc.IEEE Int.Conf.on Computer Vision and PatternRecognition(CVPR’98)1998；178-183中進(jìn)行了公開。
一旦對所有特征定義了光流矢量，即可通過對第二幅圖像施加逆光流，以便將第二幅圖像的每個特征及其周圍像素移回到所述特征在第一幅圖像中的位置，可以實現(xiàn)兩幅圖像的配準(zhǔn)。
前面已經(jīng)提到，上述方法可以用來配準(zhǔn)兩幅診斷圖像或三維圖像體積，這些圖像核磁共振成像(MRI)或計算X射線斷層成像(CT)中用相同或不同類型的掃描儀獲取。
上述公開的配準(zhǔn)方法還與用于MRI或超聲圖像之類的對比媒質(zhì)增強(qiáng)診斷成像尤其是對比度增強(qiáng)核磁共振乳腺成像(MRM)的方法結(jié)合地提供。在該方法中，為了計算乳房組織的一個區(qū)域的實際造影劑攝取量，必須將造影劑施加之前的每個體素的信號強(qiáng)度從造影劑施加之后的信號強(qiáng)度中減去。因此，將某時刻拍攝的、成像組織中不存在造影劑的第一幅圖像與第二或稍后時刻拍攝的、造影劑存在且灌注在成像組織中的第二幅圖像配準(zhǔn)，有助于在很大程度上成功地抑制偽影，如果不進(jìn)行配準(zhǔn)，這些偽影將因為患者的運動而出現(xiàn)在所述的減影圖像數(shù)據(jù)中。這個簡單的思想是，把所述的兩幅圖像相減，將使得與造影劑匯集的組織區(qū)相對應(yīng)的像素或體素的強(qiáng)度水平很高。所有其他不存在造影劑的區(qū)看起來將是一些低強(qiáng)度像素或體素，且具體地呈暗灰色或黑色。
進(jìn)一步的改進(jìn)是從屬權(quán)利要求的主題。



從下文對本方法的更詳細(xì)的描述以及所附繪圖和圖表可更清楚地看到本發(fā)明的特征及派生優(yōu)點，其中 圖1示出了乳房MRI中使用的乳房固定方法的示意性例子。
圖2示出減運算過程以及不執(zhí)行配準(zhǔn)步驟時出現(xiàn)的運動偽影左上和右上部分是同一MR掃描儀位置的兩幅乳房截面圖。測量第一幅圖像后約2.5分鐘施加造影劑，施加造影劑之后測量第二幅圖像。兩幅圖像都只是整個三維截面圖像組的典型截面。
圖3-6示意性地說明了由像素陣列形成的二維圖像，其中特征選擇標(biāo)準(zhǔn)由示例給出。
圖7-10示意性地說明了特征跟蹤是如何實現(xiàn)的。
圖11通過示出對比媒質(zhì)施加前后不同時刻拍攝的兩幅乳房圖像、執(zhí)行所述圖像的特征跟蹤和配準(zhǔn)之后對兩幅圖像相減得到的圖像，并將這些圖像與本圖左側(cè)的圖2中的例子進(jìn)行對比，說明了本發(fā)明的結(jié)果。圖2和圖11的所有圖像均取自同一MR試驗，并且沒有進(jìn)行任何除配準(zhǔn)之外的進(jìn)一步的處理。

具體實施例方式 雖然本發(fā)明的一個相關(guān)方面是將Lucas&Kanade算法從二維推廣到三維，下文的圖例只關(guān)于二維圖像以便理解起來更簡單。所述方法到三維圖像的推廣由用數(shù)學(xué)公式表示的一般性描述中給出。
工作于三維從而工作于表示三維體積數(shù)據(jù)的截面圖像組中的自動特征選擇，實質(zhì)上是對Lucas&Kanade的二維特征檢測算法的推廣，并在下面進(jìn)行更詳細(xì)的解釋 將很短時差內(nèi)記錄的同一患者的兩個圖像體積分別定義為I和J。
在第一圖像體積I內(nèi)，必須識別一定數(shù)量的已選取體素，它們將被認(rèn)為兩個體積間要跟蹤的特征。
令I(lǐng)為原始/第一圖像體積，J為下一體積，要求在其中尋找對應(yīng)的特征。
第一步驟由在第一圖像體積內(nèi)計算最小特征值來構(gòu)成。其如下實現(xiàn) 將第一圖像體積轉(zhuǎn)換成三個梯度體積Ix，Iy，Iz(描述圖像體積的笛卡兒坐標(biāo)系統(tǒng)的x，y和z方向)，其中x，y和z為每個單個體素的坐標(biāo)。
梯度體積是根據(jù)每個體素相對于其包含在所謂窗口內(nèi)的相鄰體素的強(qiáng)度梯度來計算的，所述窗口一般在x方向上的大小為(2x+1)、y方向上的大小為(2y+1)以及在z方向上的大小為(2z+1)，其中x，y，z＝1，2，3，4，…，n個體素。
考慮中心在目標(biāo)體素的3×3×3體素陣列的鄰域的大小，則所述三個梯度體積被定義為 a)X方向上的梯度體積 b)Y方向上的梯度體積 c)Z方向上的梯度體積 作為下一步驟，對每個體素，使用如上表示的所有三個之前計算出的梯度體積，來構(gòu)造梯度矩陣G。如下生成梯度矩陣 其中 對于每個體素的梯度矩陣，如下計算該梯度矩陣的最小特征值λm 定義Graphics Gems IV，Paul S.Heckbert，1994，Academic Press，S193-98， c＝m200·m020e＝m110·m101 f＝m101·m101 p＝-m200-m020-m002 q＝c+(m200+m020)m002-d-e-f r＝(e-c)m002+dm200-2(m110·m011·m101)+f·m020 其中 注由于G是個中心矩矩陣，故
梯度矩陣G的特征值計算為 然后定義最小特征值λm的閾值。
之后，應(yīng)用下述標(biāo)準(zhǔn)來考慮體素是否代表可跟蹤特征 1.如果λm大于閾值，則將該體素的實際位置保存在要跟蹤特征的列表中。
2.如果λm小于所有最小特征值λm中最大值的某個百分比，則將其從要跟蹤特征的列表中刪除。
3.如果在實際特征周圍的限定區(qū)域(它是可調(diào)的，如3×3×3)中存在另一特征，則將具有較小最小特征值λm的特征從要跟蹤特征的列表中刪除。
4.如果特征位置周圍的三維鄰域的平均信號強(qiáng)度值(例如1)小于或等于取值小得接近零(如10)的平均信號強(qiáng)度值閾值，則認(rèn)為該特征位于乳房外的周圍空氣中，并將其從要跟蹤特征的列表中去除。
5.最后，保留可限定最大數(shù)量的特征。如果要跟蹤特征的列表中存在更多的特征，則從列表中刪除那些具有較小最小特征值λm的特征。
對于二維情況，這些步驟在圖3-6中進(jìn)行了可視的解釋。這種情況下，像素陣列用小方格陣列表示。對于與具有滿足第一條標(biāo)準(zhǔn)的最小特征值的選定特征的位置對應(yīng)的不同像素進(jìn)行標(biāo)記。在圖3中，這些特征和位置通過賦予對應(yīng)像素不同外觀來突出。在圖4中，5×5鄰域是通過在這種5×5像素陣列中的內(nèi)切圓來突出。
圖1中給出了用于固定患者乳房的典型裝置。通常，MR成像在1.0或1.5特斯拉成像儀上完成，并且使用了廠家的雙乳線圈(在美國則經(jīng)常使用單乳線圈設(shè)備)，而患者處于俯臥姿勢。該治療方案至少包含任意方向上的動態(tài)對比度增強(qiáng)T1加權(quán)序列，且采樣時間在30s到180s之間(各個方案有所不同，此外還可能進(jìn)行脂肪抑制)。在以0.1mmol/kg-0.2mmol/kg體重的劑量來快速濃注釓噴酸二葡甲胺鹽或任何其他的順磁性MR造影劑之前，至少進(jìn)行一次測量，注射之后則要進(jìn)行多次測量。
在造影劑注入肘前靜脈并因而到達(dá)乳房組織之前的時刻t0，至少要獲取第一MRI圖像體積。而在緊接著造影劑注射后的時刻t1至少要獲取另一圖像體積。通常，獲取的是時間延遲體積的序列(例如，在施加造影劑之后在某行對同一體積采集次數(shù)多達(dá)六次)。
為了簡化當(dāng)前圖表，從截面圖像的三維體積中取出的兩幅圖像具有一定的(切片)厚度。因此，討論體素陣列而非像素更為準(zhǔn)確。每個體素具有與對應(yīng)信號強(qiáng)度相關(guān)的預(yù)定灰度級。在圖3-10所示例子中，網(wǎng)格代表體素陣列。網(wǎng)格的每個單元就是形成圖像的體素。
以圖3為例，其中區(qū)分了作為標(biāo)志的四種體素，它們按照對應(yīng)梯度矩陣的最小特征值進(jìn)行排序。用“102”表示的體素是具有最大特征值的體素，用“202”表示的體素是具有較大特征值的體素，用“302”表示的體素是具有較小特征值的體素，而用“402”表示的體素是具有最小特征值的體素。
在圖4中，對于每個體素，定義了以選定為特征的像素為中心、直徑為5像素/體素的歐氏距離的鄰域。所述鄰域的選取在圖4中用圓圈表示出來。
在圖4中，第一種情況用502’表示的較粗圓圈突顯出來。該處有兩個選定為特征的像素位于上述用像素的歐氏距離所定義的同一鄰域內(nèi)。這兩個像素用202’和302’表示，按照上述定義，選定為特征的像素202’的梯度矩陣的特征值大于第二個像素302’。因此，該第二個像素302’就被去除。
圖5給出了進(jìn)一層的步驟，圖3中選定為特征的像素302’不再出現(xiàn)。此處在圖5中，像素202”設(shè)有表示其鄰域的圓圈。此處也出現(xiàn)了與圖4中相似的情況，即在像素202”的鄰域內(nèi)出現(xiàn)像素302”，該像素302”的梯度矩陣的特征值小于像素202”的梯度矩陣特征值。
因此，在這種情況下，現(xiàn)在將選定為特征的像素302”從要跟蹤特征的列表中去除。
最后，在圖6的像素陣列中，選定為特征的原始像素只有一些得到保留?，F(xiàn)在在單個鄰域內(nèi)再也找不到一個以上的特征。因此，可跟蹤特征在整個圖像/體積上得到擴(kuò)散。
根據(jù)本發(fā)明的進(jìn)一層的步驟，像素的強(qiáng)度要進(jìn)一步考察，并將那些強(qiáng)度小于預(yù)定值且選定為特征的像素也從對應(yīng)于要跟蹤特征的像素的列表中去除。近似為零的很低的閾值被用于表示周圍空氣中發(fā)現(xiàn)的特征。這個步驟沒有在圖中示出，但它對于普通技術(shù)人員是顯而易見的。
一旦確定了被認(rèn)為是二維或三維圖像的特征的像素或體素，隨后就得進(jìn)行特征跟蹤。根據(jù)本發(fā)明的優(yōu)選實施例，這是通過使用所謂的“Lucas & Kanade特征跟蹤器的金字塔式實現(xiàn)”來進(jìn)行的。Jean YvesBouguet在其由Intel公司微處理器研究實驗室出版的文章“PyramidalImplementation of the Lucas Kanade Feature Tracker，Description of theAlgorithm”中給出了詳細(xì)的描述。
在上述出版物中公開了其二維例子的理論在下文中經(jīng)過進(jìn)一步改動從而適用于三維情況。下面的步驟必須對表示要跟蹤特征的每個體素進(jìn)行，各體素按照上述方法和算法進(jìn)行了區(qū)分。
定義u為體積I中的一個點，v為體積J中的對應(yīng)點。在第一步，必須由體積I和J構(gòu)造兩個金字塔，且IL和JL代表第L(0，…，m)級的體積，塔基(L＝0)表示原始體積。每個下一層體積的大小通過將直接相鄰的一定數(shù)量的像素合并成一個均值而縮減。其余層的數(shù)目由體積的大小來決定。
下一步，在最高層初始化所謂的全局金字塔運動矢量g 隨后對所有的級L計算最終的運動矢量dL。
這是按照如下步驟實現(xiàn)的 a)在每個體積IL內(nèi)，點uL的位置計算為 式中p為實際的體積位置 b)在體積的金字塔表示的每級L，沿笛卡爾坐標(biāo)x、y、z的每個方向計算梯度。
對于笛卡爾坐標(biāo)x，按照下述等式計算體積梯度 該等式相應(yīng)地也適用于其他兩個笛卡爾坐標(biāo)y和z c)隨后用上述體積梯度按照下式計算梯度矩陣 此處，值ω定義了影響體素的區(qū)域尺寸或鄰域。
矩陣元素以類似于按照前文定義得到的矩陣元素的方式來獲得，即 m110＝IΔx(ux，uy，uz)·IΔy(ux，uy，uz) m101＝IΔx(ux，uy，uz)·IΔz(ux，uy，uz) m001＝IΔy(ux，uy，uz)·IΔy(ux，uy，uz) d)對于每個級L，初始化迭代矢量
e)其后對于k取從1到最大計數(shù)K或者直到

的最小偏移量，計算如下的體積差 根據(jù)下式計算失配矢量 之后如下定義光流矢量 借助上面定義的光流矢量可將迭代運動矢量計算為 并且使其等于級L的最終光流 上述步驟對每個級重復(fù)進(jìn)行，直到最后一級L＝0。
這樣對于下一級L-1，全局光流由上述計算的值來計算，為 f)結(jié)果，最終的光流矢量為 d＝g0+d0 g)現(xiàn)在，與體積I中的點U對應(yīng)的體積J中的點v的坐標(biāo)可計算為 v＝u+d 圖7-10試圖用圖形示意性地說明上述跟蹤方法的效果，當(dāng)然，這只是真實情況的簡化表示。
像在圖3-6中的例子一樣，在圖7-10中不同時刻T0和T1得到的兩幅圖像I和J被表示成二維像素陣列。
對應(yīng)于要跟蹤特征并按上面公開的方法選取的像素20在圖像I的二維像素陣列中指出。成像物體(此例中為乳房)不發(fā)生運動時，圖像J中相同位置的像素20”將對應(yīng)于像素20。像素20”的某像素鄰域用灰色陰影方框表示，對應(yīng)于中心在像素20”的7×7像素陣列，并用數(shù)字10表示。
圖像I的特征20的新位置用20’指示。這意味著像素20的位置處的乳房組織隨著時間推移移動到像素20’的位置。當(dāng)然，圖中給出的所述像素20’僅僅出于解釋跟蹤的實現(xiàn)方式的目的。實際情況中該像素是未知的，否則就不必進(jìn)行跟蹤了。
前面用數(shù)學(xué)公式描述的、跟蹤對應(yīng)于選定特征的像素的一般方法在下文中通過圖中的實例加以描述，并且為簡便起見限于二維情況。不過，記住跟蹤的一般數(shù)學(xué)描述的普通技術(shù)人員將能理解和體會三維情況下跟蹤發(fā)生的方式。
跟蹤開始時，假定圖像I中的像素20位置不變，即具有和圖像J中的像素20”相同的位置。該像素20”在下文中就稱為跟蹤像素，而選定為特征的圖像I的像素20就稱為特征像素。應(yīng)用本算法時，首先計算對應(yīng)于跟蹤像素20”的鄰域、在本例情況下即中心在圖像J中的像素20”處的所述7×7像素陣列(該值可調(diào))的某區(qū)域的逆梯度矩陣。這個鄰域?qū)?yīng)于所謂的跟蹤區(qū)域。通過應(yīng)用已改編從而適用于三維圖像體積并將稍后詳加解釋的Lucas & Kanade的特征跟蹤算法的金字塔式實現(xiàn)，就可找到跟蹤特征在圖像J中的新位置。
圖像J中跟蹤特征在20’的新位置決定了光流矢量。這在第一個迭代步驟中定義和計算。比較圖7和8可見，跟蹤區(qū)域以使其中心更靠近相同組織區(qū)的新位置20’的方式發(fā)生了位移。
通過重復(fù)計算表示圖像I中特征的像素20與圖像J中跟蹤區(qū)域中心的新識別位置之間的失配矢量，進(jìn)一步迭代地進(jìn)行跟蹤。這就導(dǎo)致新的光流矢量定義以及所述像素鄰域或跟蹤區(qū)域的新位移，如圖9所不。
確定圖像I中特征像素20和跟蹤區(qū)域的中心之間的失配矢量的步驟一直執(zhí)行，直到失配矢量達(dá)到某個可調(diào)節(jié)的最小值。
這種情形如圖10所示。參見圖7-9可知，像素20’的位置保持不變以便突出這一事實，即移動跟蹤區(qū)域，使其最終以新正確位置20’為中心，該新正確位置20’和圖像I中的像素20表示相同的乳房部位。
上述例子僅僅用了由一個像素代表的一個特征來實現(xiàn)。上述步驟必須針對圖像中的每個選定特征進(jìn)行。
跟蹤所有的特征到其新位置后，整個圖像/體積中、不過僅限于代表乳房組織的區(qū)域內(nèi)的不同位置處存在一定量的光流矢量。所述算法相對于其他配準(zhǔn)算法的最大進(jìn)步在于，它將光流矢量遍布圖像/體積地展開。
其后，通過對圖像/體積J的每個特征及其周圍鄰域施加逆光流矢量，來實現(xiàn)變形(morphing)。這樣得到新的虛擬圖像/圖像體積J’。
圖11展示了所述方法處理實際診斷圖像的效果。左側(cè)第一幅圖像反映的是未施加造影劑在時間碼t0時進(jìn)行的乳房試驗的典型切片。左側(cè)第二幅圖像是在施加造影劑后的時間碼t1時且在相同的切片位置處獲取的切片。左側(cè)第三幅圖像顯示了從時間碼t1時間碼t0時的每個像素的信號強(qiáng)度的逐像素相減所得到的結(jié)果，即所謂的減影圖像。在時間間隔t1-t0內(nèi)，如果不存在乳房的相對移動，所得圖像在沒有造影劑的地方將包含較寬部分的信號強(qiáng)度為零且呈黑色的像素。時間間隔t1-t0內(nèi)的運動效果是，由于圖像/體積I和J的相同位置的像素不對應(yīng)于相同的成像組織區(qū)域，而產(chǎn)生偽影。在本例中，減影圖像顯示有個乳房包含對比度增強(qiáng)腫瘤，如小箭頭所示。大多數(shù)類型的惡性和良性乳房腫瘤均表現(xiàn)出由血管增生造成的附加造影劑攝取量。圖中雙乳頂部的另外白色結(jié)構(gòu)屬于運動偽影。雙乳中央的其余白色結(jié)構(gòu)是由正常乳房軟組織的微少造影劑攝取量和運動結(jié)合所造成的。
圖11的右側(cè)展示了應(yīng)用所述配準(zhǔn)方法得到的相應(yīng)圖像。時間碼t0時的圖像沒有被本算法改變，因而仍保持原樣。但是經(jīng)過三維特征跟蹤后，計算出了新的虛擬圖像體積。右側(cè)時間碼t1’時的第二幅圖像展示了實行配準(zhǔn)后相同切片位置處的對應(yīng)虛擬圖像。右側(cè)最后一幅圖像顯示了時間碼t1和時間碼t0時的每個像素的信號強(qiáng)度的逐像素相減所得到的結(jié)果，即新的(虛擬的)減影圖像。仍然有個乳房示出了對比度增強(qiáng)腫瘤，如小箭頭所示。乳房內(nèi)代表正常乳房軟組織的區(qū)域仍然反映出適度的造影劑攝取量，并且由于運動偽影的減少而不那么亮。由于運動偽影補(bǔ)償，雙乳頂部較寬部分的白色結(jié)構(gòu)消失了。
圖11的例子是一種簡化的情況，這是因為僅僅考慮了兩個圖像的序列，且它們分別在施加造影劑前和后獲取。
實際上，可以在造影劑注射前后每隔固定時間間隔來采集一系列截面圖像(例如，超過50幅)，并可在每對連續(xù)圖像之間計算減影圖像。
除了在一個系列內(nèi)采集一幅以上的截面圖像外，也可一個系列接一個系列地采集兩個系列以上(如5個甚至更多)的圖像。這就得到一些減影圖像/體積。當(dāng)然，本公開算法可以應(yīng)用于所有的系列組合，而不管這些系列是在造影劑施加前或施加后獲取的，目的是減小其間發(fā)生的任何運動。
在下面的附錄中，給出了Lucas & Kanade跟蹤算法的金字塔式實現(xiàn)的詳細(xì)理論。
依照本發(fā)明的上述方法可以通過被編碼成能由計算機(jī)運行的程序的算法來實現(xiàn)。
該程序可以保存在例如磁性、光磁或光學(xué)基片的便攜式或靜態(tài)存儲設(shè)備上，如一張或多張軟盤、一張或多張磁光盤、一張或多張CD-ROM或DVD-ROM，或者也可以保存在移動或固定硬盤上。另外，該程序也可以保存在譬如筆式驅(qū)動或類似設(shè)備之類的固態(tài)電子設(shè)備上。
編碼程序的偽代碼例子如下 //build gradient volumes in X，Y and Z direction； buildGradientVolumes() { for all positions in volume do { X-GradientVolume@position＝ GreyValue(originalVolume@position(x-1，y，z))- GreyValue(originalVolume@position(x+1，y，z))； Y-GradientVolume@position＝ GreyValue(originalVolume@position(x，y-1，z))- GreyValue(originalVolume@position(x，y+1，z))； Z-GradientVolume@position＝ GreyValue(originalVolume@position(x，y，z-1))- GreyValue(originalVolume@position(x，y，z+1))； } } ////////////////////////////////////////////////////////////////////////////////////////////////// ////// //selection of features ////////////////////////////////////////////////////////////////////////////////////////////////// ////// for all Voxel in First Volume do { //build matrix with X-，Y-and Z-gradients； gradientMatrix＝buildGradientMatrix(position(x，y，z)， i_SearchRegionSize)； calcMnimumEigenvalue of gradientMatrix； if(MnimumEigenvalue＞i_GlobalMnimumQuality) { addToFeatureList(x，y，z，MnEigenValue)； } } for all Features in FeatureList from highest to lowest eigenvalue to { if(Feature-＞MnEigenValue＜i_percentvalue of OverallMaximumEigenvalue){ deleteFeatureFromList()； } if(distance to another Feature＜i_MnDistance){ deleteFeatureFromList()； } } if(Number of Features＞i_MaximumNumberOfFeatures) { deleteFromList(FeaturesWithSmallerMnEigenValue)； } ////////////////////////////////////////////////////////////////////////////////////////////////// ////// //tracking of selected features ////////////////////////////////////////////////////////////////////////////////////////////////// ////// for each Feature in FeatureList do { //build matrix from X-，Y-and Z-gradients of the original Volume at feature position GradientMatrix＝buildGradientMatrix(featurePosition， i_TrackingRegionSize)； invertedGradientMatrix＝invert(GradientMatrix)； do { if(FeatureMovementVector＞i_MaximumAllowedDisplacement) { featureTracked＝false； break； } for all Offsets in i_TrackingRegionSize do { diff＝ GreyValue(OriginalVolume@OriginalFeaturePosition+Offset) - GreyValue(MovedVolume@TrackedFeaturePosition+Offset)； vectorX＝vectorX+diff* GreyValue(X-GradientVolume@actualFeaturePosition)； vectorY＝vectorY+diff* GreyValue(Y-GradientVolume@actualFeaturePosition)； vectorZ＝vectorZ+diff* GreyValue(Z-GradientVolume@actualFeaturePosition)； } mismatchVector＝mathVector(vectorX，vectorY，vectorZ)； currentMovementVector＝invertedGradientMatrix* mismatchVector； FeatureMovementVector＝FeatureMovementVector+ currentMovementVector； if(currentMovementVector＜i_MnimumDisplacement){ featureTracked＝true； break； } if(NumberOfIterations＞i_MaximumAllowedIterations){ featureTracked＝false； break； } NumberOfIterations＝NumberOfIterations+1； } While(true) } ////////////////////////////////////////////////////////////////////////////////////////////////// ////// //calculation of gradient matrix function buildGradientMatrix(position，regionSize) { for all Offsets in regionSize do { Ix＝GreyValue(X-GradientVolume@positon+Offset)； Iy＝GreyValue(Y-GradientVolume@positon+Offset)； Iz＝GreyValue(Z-GradientVolume@positon+Offset)； Ixx＝Ixx+Ix*Ix； Iyy＝Iyy+Iy*Iy； Izz＝Izz+Iz*Iz； Ixy＝Ixy+Ix*Iy； Ixz＝Ixz+Ix*Iz； Iyz＝Iyz+Iy*Iz； } return gradientMatrix(Ixx Ixy Ixz Ixy Iyy Iyz Ixz Iyz Izz)； }
權(quán)利要求
1.用于減少物體運動造成的成像偽影的生物醫(yī)學(xué)圖像配準(zhǔn)方法，其包括以下步驟
a)提供同一物體的至少第一幅和第二幅數(shù)字或數(shù)字化圖像，或者至少第一組和第二組該物體的截面圖像，所述圖像由像素或體素的二維或三維陣列形成；
b)通過選取一定數(shù)量的被設(shè)為標(biāo)志或特征的像素或體素，在第一幅圖像或第一組圖像中定義一定數(shù)量的標(biāo)志，即所謂的特征，并生成要跟蹤的所述特征的列表；
c)通過對被選定為特征的每個像素或體素確定從第一幅到第二幅圖像或從第一組到第二組圖像的光流矢量，從第一幅到第二幅圖像或從第一組到第二組圖像跟蹤選定為特征的每個像素或體素的位置；
d)通過對第二幅或第二組圖像的像素或體素施加逆光流，來對第一幅和第二幅圖像、或?qū)Φ谝唤M和第二組圖像進(jìn)行配準(zhǔn)；
其特征在于
執(zhí)行自動特征選擇步驟，該步驟包括
B1)在第一幅圖像或第一組截面圖像的每個像素或體素周圍定義像素或體素鄰域，所述的像素或體素鄰域包括有限數(shù)量的像素或體素；對于單幅圖像，選用二維鄰域；對于截面圖像組，則選用三維鄰域；
B2)對于每個像素或體素，確定所述像素或體素鄰域的所有像素或體素的平均信號強(qiáng)度值；
B3)定義平均信號強(qiáng)度值閾值；
B4)比較在步驟B2中為每個像素或體素鄰域確定的平均信號強(qiáng)度，并將所述平均信號強(qiáng)度值與平均信號強(qiáng)度值閾值進(jìn)行比較；
B5)如果步驟B4中所述鄰域的平均信號強(qiáng)度值大于閾值，則將所述像素或體素定義為要跟蹤特征，并將其添加到要跟蹤特征列表中。
2.根據(jù)權(quán)利要求1所述的方法，其特征在于，借助實驗，通過在一組測試或樣本圖像上執(zhí)行所述方法，并比較取不同閾值時獲得的結(jié)果，來確定所述平均信號強(qiáng)度閾值。
3.根據(jù)權(quán)利要求1或2所述的方法，其特征在于，通過所謂的Lucas & Kanade自動特征跟蹤方法或算法，來自動地跟蹤被選作在第二幅圖像中要跟蹤的特征或標(biāo)志的圖像的像素或體素。
4.根據(jù)權(quán)利要求3所述的方法，其特征在于，將發(fā)明出來僅作用于二維圖像的Lucas & Kanade算法的原始實現(xiàn)推廣到作用于三維圖像體積而非二維圖像，以適用于生物醫(yī)學(xué)截面圖像體積，其包含如下步驟
a)提供由三維體素陣列構(gòu)成的第一三維圖像體積I；
b)對于每個目標(biāo)體素定義體素鄰域，該鄰域包圍所述目標(biāo)體素，并具有中心在所述目標(biāo)體素的體素陣列，該體素陣列具有可調(diào)節(jié)的大小，例如3×3×3；
c)定義笛卡爾坐標(biāo)系統(tǒng)，以及相對于所述笛卡爾坐標(biāo)系統(tǒng)的每個軸，計算所謂的梯度體積，由第一圖像體積的每個體素相對于其鄰域體素的強(qiáng)度梯度來形成所述梯度體積，且定義為
a)X方向的梯度體積
b)Y方向的梯度體積
c)Z方向的梯度體積
d)對于每個目標(biāo)體素，計算所謂的梯度矩陣，所述梯度矩陣定義為
其中
e)對于三維圖像體積的每個目標(biāo)體素的每個梯度矩陣，應(yīng)用下式來計算最小特征值λm
定義
c＝m200·m020e＝m110·m101 f＝m101·m101 p＝-m200-m020-m002
q＝c+(m200+m020)m002-d-e-f
r＝(e-c)m002+dm200-2(m110·m011·m101)+f·m020
其中
注意由于G是中心矩矩陣，所以
如下，計算所述梯度矩陣的最小特征值
λm＝minimum(λ1，λ2，λ3)
f)定義最小特征值λm的閾值；
g)將對應(yīng)梯度矩陣的最小特征值λm滿足一定標(biāo)準(zhǔn)的每個體素選擇作為圖像I中要跟蹤的特征，這些標(biāo)準(zhǔn)如下
i)λm大于所述閾值；
ii)λm不小于所有最小特征值λm中的最大值的某個百分比；
iii)如果在保持為選定特征的目標(biāo)體素周圍的已定義體素鄰域內(nèi)存在另一個被選作特征的體素，則僅將其梯度矩陣具有較大的最小特征值的一個體素選為特征，而將另一個體素從可跟蹤特征列表中刪除；
iv)選作特征的體素周圍的三維鄰域的平均信號強(qiáng)度值大于可調(diào)節(jié)平均信號強(qiáng)度值閾值；
5.根據(jù)權(quán)利要求4所述的方法，其特征在于，選作特征的體素的鄰域的平均信號強(qiáng)度的閾值被設(shè)為大約10。
6.根據(jù)上述一項或多項權(quán)利要求所述的方法，其特征在于，對于第一圖像體積I中被選作特征的每個體素，相對于所述特征在第二較晚時刻得到的相同物體的另一圖像J的體素陣列中的位置，對所述特征進(jìn)行跟蹤，且通過所謂的Lucas & Kanade特征跟蹤器的金字塔式實現(xiàn)來執(zhí)行所述跟蹤。
7.根據(jù)權(quán)利要求6所述的方法，其特征在于，所謂的Lucas &Kanade特征跟蹤器的金字塔式實現(xiàn)包含如下步驟
將點u定義為對應(yīng)于圖像I的三維體素陣列的圖像體積I中的點，將v定義為對應(yīng)于圖像J的三維體素陣列的圖像體積J中的對應(yīng)點，這些點具有圖像I中被選作特征的體素的坐標(biāo)；
-由圖像體積I和J構(gòu)造兩個理想金字塔，且IL和JL表示第L(0，…，m)級體積，塔基(L＝0)表示原始體積；
-每個下一層體積的大小通過將在直接鄰域中的一定數(shù)量的像素合并成一個均值而縮減；
-定義所謂的全局金字塔運動矢量g，并在最高層Lm上對其進(jìn)行初始化
i)將點u的位置定義為uL，并通過以下等式來計算所述位置，其中L具有金字塔的實際層或級的值
式中p為實際的體積位置
ii)根據(jù)體積梯度的以下等式，計算笛卡爾坐標(biāo)x、y和z的每個方向上的體積梯度
iii)使用梯度體積，按照下式來計算梯度矩陣
其中
m110＝IΔx(ux，uy，uz)·IΔy(ux，uy，uz)
m101＝IΔx(ux，uy，uz)·IΔz(ux，uy，uz)
m001＝IΔy(ux，uy，uz)·IΔy(ux，uy，uz)
此處，值ω定義了影響表示跟蹤特征的目標(biāo)體素的區(qū)域尺寸或體素鄰域。
iv)對于每個級L，初始化如下定義的迭代矢量
v)對于k取從1到最大計數(shù)K或者直到
的最小偏移，計算如下的體積差
vi)根據(jù)下式計算失配矢量
vii)按照下式確定光流矢量
式中，G-1為在步驟iii)中確定的梯度矩陣G的逆矩陣
viii)使用上面定義的光流矢量如下計算迭代運動矢量
并且使其等于級L的最終光流
ix)對各級重復(fù)步驟i)-viii)，直到最后級L＝0，得到下式定義的最終的光流矢量
d＝g0+d0
x)通過下式，確定與代表體積I中要跟蹤特征的點u對應(yīng)的體積J中的點v的坐標(biāo)
v＝u+d
xi)對于對應(yīng)于圖像I中選定特征的每個體素，重復(fù)執(zhí)行上述方法。
8.根據(jù)權(quán)利要求7所述的方法，其特征在于，兩幅圖像I和J的配準(zhǔn)通過如下方式實現(xiàn)對圖像J中識別為與表示某特征的體素相對應(yīng)的每個點施加逆光流矢量，其中該特征與圖像I中對應(yīng)于某選定特征的體素的點U相對應(yīng)；或者對與圖像I中識別為選定特征的體素對應(yīng)的每個點u施加光流矢量。
9.根據(jù)前面的一項或多項權(quán)利要求所述的方法，其特征在于，它與用于對比媒質(zhì)增強(qiáng)診斷成像，尤其是對比媒質(zhì)增強(qiáng)MRI或超聲成像的方法結(jié)合提供，該方法包含下列步驟
a)提供由MRI或超聲獲取的同一物體的至少第一幅和第二幅數(shù)字或數(shù)字化圖像，或者至少第一組和第二組該物體的截面圖像，所述圖像由像素或體素的二維或三維陣列形成，其中對于相同組織或組織區(qū)或者相同的解剖區(qū)的掃描，是在第二時刻或者任何隨后時刻、且在所述組織或組織區(qū)或者所述解剖區(qū)中存在對比媒質(zhì)時進(jìn)行的；
b)在第一幅圖像或第一組圖像中，通過選取一定數(shù)量的像素或體素，將該像素或體素設(shè)置作為標(biāo)記或特征，來定義一定數(shù)量的標(biāo)志，即所謂的特征，并生成要跟蹤的所述特征的列表；
c)通過對選定為特征的每個像素或體素確定從第一幅到第二幅圖像或從第一組到第二組圖像的光流矢量，來從第一幅到第二幅圖像或從第一組到第二組圖像跟蹤選定為特征的每個像素或體素的位置；
d)通過對第二幅或第二組圖像的像素或體素施加逆光流，來對第一幅和第二幅圖像或第一組和第二組圖像進(jìn)行配準(zhǔn)。
10.根據(jù)權(quán)利要求1-9所述的方法，其特征在于，在所述特征選擇步驟之后以及執(zhí)行所述特征跟蹤步驟之前，進(jìn)行進(jìn)一步的自動特征選擇步驟，該步驟包括
B1)在第一幅圖像或第一組截面圖像的每個像素或體素周圍定義像素或體素鄰域，所述的像素或體素鄰域包括有限數(shù)量的像素或體素；對于單幅圖像，選用二維鄰域；對于截面圖像組，則選用三維鄰域；
B2)對于每個像素或體素，確定所述像素或體素的鄰域的所有像素或體素的平均信號強(qiáng)度值；
B3)定義平均信號強(qiáng)度值閾值；
B4)比較在步驟B2中確定的每個像素或體素鄰域的平均信號強(qiáng)度值，并將所述平均信號強(qiáng)度值與所述平均信號強(qiáng)度值閾值進(jìn)行比較；
B5)如果步驟B4中所述鄰域的平均信號強(qiáng)度值大于閾值，則將所述像素或體素定義為要跟蹤的特征，并將其添加到要跟蹤的特征的列表中。
11.用于配準(zhǔn)兩個數(shù)字或數(shù)字化圖像/圖像體積的計算機(jī)程序，在該計算機(jī)程序中對根據(jù)前面的一項或多項權(quán)利要求的所述方法進(jìn)行編碼，該計算機(jī)程序的特征在于以下偽代碼
//build gradient volumes in X，Y and Z direction；
buildGradientVolumes()
{
for all positions in volume do
{
X-GradientVolume@position＝
GreyValue(originalVolume@position(x-1，y，z))-
GreyValue(originalVolume@position(x+1，y，z))；
Y-GradientVolume@position＝
GreyValue(originalVolume@position(x，y-1，z))-
GreyValue(originalVolume@position(x，y+1，z))；
Z-GradientVolume@position＝
GreyValue(originalVolume@position(x，y，z-1))-
GreyValue(originalVolume@position(x，y，z+1))；
}
}
//////////////////////////////////////////////////////////////////////////////////////////////////
//////
// selection of features
//////////////////////////////////////////////////////////////////////////////////////////////////
//////
for all Voxel in First Volume do
{
//build matrix with X-，Y-and Z-gradients；
gradientMatrix ＝ buildGradientMatrix(position(x，y，z)，
i_SearchRegionSize)；
calcMnimumEigenvalue of gradientMatrix；
if(MnimumEigenvalue＞i_GlobalMnimumQuality)
{
addToFeatureList(x，y，z，MnEigenValue)；
}
}
for all Features in FeatureList from highest to lowest eigenvalue to
{
if(Feature->MnEigenValue ＜i_percentvalueof
OverallMaximumEigenvalue){
deleteFeatureFromList()；
}
if(distance to another Feature＜i_MnDistance){
deleteFeatureFromList()；
}
}
if(Number of Features＞i_MaximumNumberOfFeatures)
{
deleteFromList(FeaturesWithSmallerMnEigenValue)；
}
//////////////////////////////////////////////////////////////////////////////////////////////////
//////
//tracking of selected features
//////////////////////////////////////////////////////////////////////////////////////////////////
//////
for each Feature in FeatureList do
{
//build matrix from X-，Y-and Z-gradients of the original Volume
at feature position
GradientMatrix＝buildGradientMatrix(featurePosition，
i_TrackingRegionSize)；
invertedGradientMatrix＝invert(GradientMatrix)；
do
{
if(FeatureMovementVector＞i_MaximumAllowedDisplacement)
{
featureTracked＝false；
break；
}
for all Offsets in i_TrackingRegionSize do
{
diff＝
GreyValue(OriginalVolume@OriginalFeaturePosition+Offset)
-
GreyValue(MovedVolume@TrackedFeaturePosition+Offset)；
vectorX＝vectorX+diff*
GreyValue(x-GradientVolume@actualFeaturePosition)；
vectorY＝vectorY+diff*
GreyValue(Y-GradientVolume@actualFeaturePosition)；
vectorZ＝vectorZ+diff*
GreyValue(Z-GradientVolume@actualFeaturePosition)；
}
mismatchVector＝mathVector(vectorX，vectorY，vectorZ)；
currentMovementVector＝invertedGradientMatrix*
mismatchVector；
FeatureMovementVector＝FeatureMovementVector+
currentMovementVector；
if(currentMovementVector＜i_MnimumDisplacement){
featureTracked＝true；
break；
}
if(NumberOflterations＞i_MaximumAllowedIterations){
featureTracked＝false；
break；
}
NumberOflterations＝NumberOflterations+1；
}
While(true)
}
//////////////////////////////////////////////////////////////////////////////////////////////////
//////
//calculation of gradient matrix
function buildGradientMatrix(position，regionSize)
{
for all Offsets in regionSize do
{
Ix＝GreyValue(X-GradientVolume@positon+Offset)；
Iy＝GreyValue(Y-GradientVolume@positon+Offset)；
Iz＝GreyValue(Z-GradientVolume@positon+Offset)；
Ixx＝Ixx+Ix*Ix；
Iyy＝Iyy+Iy*Iy；
Izz＝Izz+Iz*Iz；
Ixy＝Ixy+Ix*Iy；
Ixz＝Ixz+Ix*Iz；
Iyz＝Iyz+Iy*Iz；
}
return gradientMatrix(Ixx Ixy Ixz
Ixy Iyy Iyz
Ixz Iyz Izz)；
}
全文摘要
公開了一種用于減少由物體運動造成的成像偽影的生物醫(yī)學(xué)圖像的配準(zhǔn)方法，在該方法中，在第一幅和第二幅待配準(zhǔn)圖像中自動判定并跟蹤一定數(shù)量的特征或標(biāo)志，以確定上述第一幅和第二幅圖像之間的光流矢量。隨后通過對第二幅圖像施加逆光流進(jìn)行配準(zhǔn)。上述自動特征選擇步驟是通過對每個像素確定其定義的鄰域內(nèi)的平均信號強(qiáng)度，隨后將該平均信號強(qiáng)度與預(yù)先確定的閾值進(jìn)行比較來實現(xiàn)的。如果所述鄰域的平均信號強(qiáng)度值大于預(yù)先確定的閾值，那么對應(yīng)的像素就被定義為特征，并添加到要跟蹤的特征的列表中。
文檔編號G06T3/00GK101133431SQ200680003997
公開日2008年2月27日 申請日期2006年1月31日 優(yōu)先權(quán)日2005年2月3日
發(fā)明者托尼·維爾納·沃姆維格, 海納·法貝爾, 迪爾克·邁耶 申請人:布拉科成像S.P.A.公司