具有全數據或部分數據的多模態醫學圖像配準

李正偉

（成都理工大學工程技術學院，四川樂山 614007）

圖像分析涉及到圖像分割、重構、分類和檢測等多種任務[1]。通過比較或融合兩個或多個圖像可以尋求重要信息，圖像之間不可避免的錯位使得圖像配準成為研究重點。圖像配準的目的是在同一場景的兩個或多個圖像之間找到最佳的空間對齊，涉及圖像配準的圖像可以使用相同的(單模態)或不同的(多模態)成像技術捕獲[2]。在計算機斷層掃描(CT)和磁共振成像(MRI)的情況下，在軸平面上拍攝器官的完整橫截面圖像。因此，將新掃描與原始掃描進行比較將需要具有部分重疊的圖像配準。由于部分重疊的兩幅圖像之間缺乏一一對應關系，因此使用最常用的配準方法很難獲得精確的配準[3]。

單模態醫學圖像配準方法的目的是對同一成像技術(如CT)得到的圖像進行配準。文獻[4]提出了基于快速傅里葉變換(FFT)技術，該技術可尋求參考圖像和感測圖像之間的最佳平移、旋轉和縮放參數。這種相位相關技術的擴展后來發展為基于Fourier-Mellin 變換(FMT)的圖像配準方法[5]。文獻[6]提出了改進的FMT 方法，消除了直角坐標到對數極坐標的轉換。評估并驗證了該方法對噪聲的準確性和魯棒性。

多模態圖像配準算法側重于尋找使用各種成像方式(如CT 和MRI)生成圖像之間的對應關系，并從不同醫學成像方式的融合中提供密集的視覺信息。文獻[7]實現了一種基于圖像強度全局和局部變化的通用圖像配準方法。該算法在全數據的多模態圖像和部分數據的單模態圖像上進行了實驗，但由于采用彈性模型估計變形場而存在特殊的缺陷。基于信息論的度量，如互信息(MI)[8]及其變化，已廣泛用于基于強度的模態間配準，但其存在收斂速度慢、精度低和對實現決策敏感等缺陷。文獻[9]提出了標準化互信息(NMI)度量，當邊緣熵之和的增長速度快于聯合熵時，可以處理由于失準而導致的MI 值的增長。

考慮到CT和MRI是兩種廣泛使用的成像方式，可以產生具有可比空間分辨率的橫截面圖像，因此MRI具有更高的對比度分辨率，而CT掃描則相對較低。常用的多模態配準方法往往由于缺失部分數據導致配準無效，降維算法（即流形學習）常用于圖像和數據處理領域，其目的是研究嵌入高維空間的低維流形[10]。

該文改進了文獻[11]中利用Laplacian 特征映射來表示多模態圖像配準的方法，實現了多模態到單模態的轉換。在多模態醫學圖像上直接應用單模態配準技術以及恢復強放縮、旋轉和平移，并且對完全重疊和部分重疊兩種類型的圖像都適用。此外，還提出了一種快速且易于使用的對齊技術，以此補償由主向量反射引起的隨機符號歧義。對齊技術有助于更復雜的優化算法以及基于強度度量的應用。提出的強度變換解決了部分重疊的多模態醫學圖像配準問題。從實驗的角度出發，提出了一套定性和定量的分析方法來檢驗系統的性能和準確性，其中臨床輸入圖像同時受到多個自由度（平移、旋轉、放縮）的約束，在不同的多模態配準中獲得了較低的平均絕對誤差。

1 數據降維

將高維數據映射到低維空間的變換降維技術可分為線性和非線性方法。該文引入主成分分析(PCA)[12]和Laplacian 特征映射[11]分別用于線性和非線性技術。光譜嵌入可作為尋找與特定矩陣的上（或下）特征值相關聯的一組特征向量。

考慮集合x1,…,xn∈M的流形M表面上有n個采樣數據點，其中xi∈RD。尋找一組對應點y1,…,yn，其中yi∈RD(d＜＜D)，用于保留X=[x1,...,xn]的特定相關信息。

1.1 主成分分析(PCA)

PCA 的目的是通過正交變換將一組可能存在相關性的變量轉換為一組線性不相關的變量[12]。在PCA 中，假設數據集的方差是對數據集中存儲信息量的度量，則方差越大，存儲的信息越多。

給定矩陣X∈RD×n，其列為D維數據點，則矩陣Xc∈RD×n定義為中心數據點的矩陣（每行的樣本均值為零）。投影由Y=AXc給出，其中A為d×n矩陣。PCA 優化目標函數為：

其中，A為d×D正交矩陣。Y的協方差矩陣為：

利用Lagrange 乘數技術，將最大化問題簡化為求解協方差矩陣V的特征值問題：

根據線性代數，任何MMT形式的矩陣都是對稱矩陣。因此，V為對稱矩陣，由此得出結論：根據譜定理，V為正交可對角化矩陣，并且只有實的非負特征值。正交可對角化矩陣V具有正交的非零實特征向量。

矩陣V具有最大數量D個特征值和大小為D×1 的D個特征向量。設λ1≥λ2≥···≥λD≥0（按降序排列），相應的正交特征向量為。矩陣V的特征向量是數據點的主成分(或主方向)，數據點沿主分量(或主方向)方差最大。圖1 給出了樣本合成數據集的主要組成部分。這些主成分構成一組標準基向量，將高維數據點投影到低維空間，同時保留大部分信息。所需的線性變換矩陣A由與前d個最大特征值相關的前d個主成分組成[13]。

圖1 樣本合成數據集的主成分

1.2 Laplacian特征映射

Laplacian 特征映射的目的是找到高維流形的低維嵌入作為結構表示。在保留局部信息的同時，提取結構信息并將其嵌入到低維空間中。結構表示得益于保留局部性和結構等價性。保留局部性是指同一流形的兩個相似塊映射到新的坐標后，其結構表示是相似的。結構等價性是指不同流形的兩個相似塊的結構表示相似。

Laplacian 特征映射基于Laplacian 圖，其算法包括3 個主要步驟[11]：1）構造鄰接圖；2）選擇權值；3）尋找特征映射。每個步驟都涉及選擇適當的實現方式并為多個參數確定適當的值。由Laplacian 特征映射所保證的保留局部性和結構等價性是該文方法的關鍵決定因素。

2 多模態圖像配準

該文提出的多模態到單模態變換的目的是將多模態輸入圖像變換成相似強度的坐標系。通過提取比單個像素強度更多的信息來研究小圖像塊的結構。由于來自不同模態的圖像的內部結構相似，因此，轉換將導致多模態圖像的相似結構表示。多模態到單模態變換的圖像配準如圖2 所示。

圖2 多模態到單模態變換的圖像配準

2.1 構建高維空間

在圖像配準問題中，只存在參考和感測兩種圖像。在高維空間中，通過從每個輸入圖像中提取小塊來構造點云。高維空間的每個維度專用于每個補丁的一個圖像像素。考慮大小為s×s的每個補丁都比前一個補丁移動了一個像素，對于大小nr×nc的圖像，高維流形在s2=D維空間中包括nr×nc=N個點。使用點的坐標以矩陣格式表示這些點集，對于大小為200×200、塊大小為3×3的圖像，矩陣表示為40 000×9的矩陣，使用圖像的小塊和相應的點云矩陣表示構建高維點云的示例，如圖3 所示。對于大小為N像素、塊為D個像素的灰度圖像，矩陣大小為N×D。

圖3 高維空間中點云的矩陣表示

D的選擇取決于圖像的應用和大小。當選擇相對較小的D時，嵌入的流形類似于原始圖像。當D選擇較大的值會導致圖像模糊。不同D值時T2-MRI圖像切片及其嵌入流形的圖像示例，如圖4所示。

圖4 不同D值時T2-MRI圖像切片及其嵌入流形的圖像

圖4(a)給出了原始T2-MRI 圖像的切片。利用Laplacian 特征映射對圖像進行了研究。圖4(b)～(d)說明了嵌入流形后參數D的選定值對渲染圖像清晰度的影響。隨著D的增加，嵌入流形的清晰度降低。從圖4(d)可以看出，使用15×15 的補丁會導致大腦結構內部細節的丟失。

2.2 流形學習

Laplacian 特征映射降維需要構造一個無向鄰接圖G=(V,E)，圖中包含高維空間RD中的一組點xi∈V和一組指定鄰域連通度的邊E。對于給定的點集，k近鄰通過在高維空間中找到每個點的k近鄰并指定該點與其每個近鄰之間的邊來構造圖[14]。通過構造有向圖，當節點i屬于節點j的k近鄰集合，而節點j不在節點i的k近鄰集合時，在Laplacian 特征映射方法中，通過在有單向邊的地方添加更多邊（稱為對稱kNN 圖）或刪除單向邊（稱為相互kNN 圖）來生成無向圖[15]。

為了構造鄰接圖，該文將k近鄰與對稱kNN 圖相結合。在降維過程中，參數j的選擇只是內存速度與確保具有單個連接之間的權衡。由于醫學掃描使用平滑的強度級別范圍，從而獲得平滑的高維點云，因此，對于200×200 像素左右的醫學掃描，能夠獲得k=20 的單個連接分量。如果輸入圖像有噪點，可能導致多個連接組件的圖，因此使用低通濾波器或增加k。

選擇每條邊的權重決定了連接的能力，該文借鑒文獻[16]中圖像的Laplacian 算子與流形上的Laplacian Beltrami(L-B)算子的可比性。為了得到L-B算子的最佳近似，該文使用以下熱核加權方式：

其中，σ∈R為唯一的參數。如果兩個節點i和j通過一條邊連接，則wij=1，否則wij=0。

相鄰點之間的距離是局部信息的度量方法。因此，通過選擇合適的參數σ值來保持各向異性指數加權方案。在式（4）的指數加權方案中，鄰近點的權重高于遠處點的權重。當σ與||xi-xj||相比相對較大時，所有連接的權值都接近于1，從而產生無權鄰接圖。相反，相對較小的σ將導致不顯著的邊緣，進而使得特征值的收斂失敗。該文提出構造與σ相關的權重矩陣S(σ2)為：

對于多個σ值，繪制S曲線與對數刻度(σ2)的關系圖。S曲線在σ=0 和σ=∞處有兩個漸近線，并用于選擇σ，其中圖形在總范圍的上半部分呈線性。為每幅圖像計算S曲線增加了額外的處理時間。該文選擇熱核的方差等于所有邊的最大平方距離：

所有邊的權重都限制在e-1/2≈0.606 和1 之間。在減少計算時間的同時，嵌入流形和配準結果均有效且可靠。

最優嵌入在保留局部信息的同時將點xi映射到點yi，即連接點xi和xj在分別映射到yi和yj之后將保持盡可能靠近。最小化以下目標函數：

對于給定的具有wij項的對稱權值矩陣W，該文定義了對角度量矩陣D，使得。度量矩陣D的項是W的列（或行）之和。該文還將Laplacian 矩陣定義為L=D-W。如果集合V包含N個點，則W以及D和L均為稀疏N×N對稱矩陣。

通過添加正交約束yTD1=0 來消除平凡解，并約束yTDy=1 來消除嵌入中的任意比例因子，將式(7)的最小化問題簡化為：

由廣義特征值問題的最小特征值解得到式（8）的解向量y：

式（9）為L-B 算子的特征分解或譜分解，將矩陣用其特征值和特征向量表示，特征值的集合為L的譜。L-B 算子特征值和特征向量具有等距不變量性質，即如果流形沒有拉伸（例如彎曲成額外的維度），光譜值不會改變。因此，兩個方向不同的流形將具有相似的譜表示。式（4）和（7）表明，L為實對稱且半正定矩陣。因此，本征函數1 和本征值λ=0是式（9）的平凡解。特征值為0 的多重性與圖的連通分量的數目相關，對于具有多個連通分量的圖，L為塊對角矩陣，其中每個塊表示連通分量的Laplacian矩陣。

考慮到Laplacian 矩陣的性質，需要去掉與特征值為0 的所有特征向量，并使用下一個具有最小非零特征值的d個特征向量來將流形嵌入到低d維歐幾里得空間中。嵌入流形的形式為N×d矩陣Y=[y1,…,yd]，其中第i行顯示第i點的嵌入坐標。

該文使用ARPACK 軟件求解廣義特征值問題的d個最小特征值，ARPACK 代表ARnoldi 包可用于逼近幾個特征值和對應的特征向量。由于矩陣L為對稱矩陣，因此，將Arnoldi 過程簡化為Lanczos 過程的變體，即隱式重新啟動Lanczos 方法(IRLM)[17]。由特征解算器進行的計算近似使得最小特征值不完全等于零。與大于0 的最小特征值相對應的第一個d特征向量是在低維空間嵌入流形的正交基。此時，尺寸為N×1 的d個特征向量在d維空間中定位嵌入流形的N個點。

輸入圖像的結構表示是通過重新排序特征向量到原始輸入圖像的大小。每個特征向量表示具有特定派生度量的輸入圖像。特征向量的指數越高，結構表示的導數較高。使用特征圖像1、6、11 和18 的人腦T1-MRI 的結構表示，如圖5 所示。Laplacian 特征映射降維的必要步驟，如算法1 所示。

圖5 用Laplacian特征圖表示T1-MRI的結構

算法1：Laplacian 特征映射降維（結構表示）

輸入：高維空間中N個點的集合X=(x1,…,xN)T∈RN×D

輸出：低維空間中N個點的集合Y=(y1,…,yN)T∈RN×d

步驟1：計算每兩個數據點之間的距離；

步驟2：選擇k和對稱kNN 圖，構造k近鄰構造鄰接圖；

步驟3：對于任意的xi,xj∈V，定義σ2=max(||xi-xj||2)；

步驟4：采用熱核加權方案為每兩個鄰居點之間的每個邊分配權重；

步驟5：構造稀疏實對稱N×N矩陣W,D,L；

步驟6：從L中找出連接組件的數量nc；

步驟7：求解(nc+d) 個特征值和相應特征向量的廣義特征值問題Ly=λDy；

步驟8：對特征向量進行排序，以遞增順序表示特征值；

步驟9：去掉第一個nc個特征向量；

步驟10：返回剩余的d個特征向量。

2.3 流形對齊

Lanczos 算法以隨機初始向量開始初始化，雖然Laplacian 特征映射保留了局部信息，但不能確保嵌入流形在低維空間的對齊。因此，需要進行流形對齊[18]。

PCA 提供由平移和旋轉組成剛體變換的參數。由于任何特征向量都是有效的特征向量，因此兩個嵌入流形之間的映射可能需要旋轉。采用PCA 的剛體變換無法處理這種隨機符號模糊問題。如果兩個流形的主成分遵循由向量的叉積在垂直主向量之間定義的右手或左手規則，PCA 將使流形具有足夠的精度確保對齊。兩個嵌入流形的右手規則檢查，如圖6 所示。其中，參考流形和感測流形的前3 個主方向分別遵循左手和右手規則。

圖6 兩個嵌入流形的右手規則檢查

該文提出了一種快速而直接的方法來檢查參考流形和感測流形是否符合右手規則，每兩幅對應的特征圖像所對應的協方差的符號是它們關聯方向的度量。如果它們有負關聯，則將-1 與感知嵌入流形的特征圖像相乘。當使用該文方法補償反射誤差時，PCA計算剛體變換矩陣R，其對齊中心嵌入流形為：

其中，矩陣A和矩陣B分別為感測流形和參考流形的主方向。由于矩陣A為正交矩陣，根據PCA的性質，旋轉矩陣簡化為：

通過將兩個流形的中心移動到原點，然后旋轉感測流形的主方向以對準參考流形的主方向來實現對齊。利用PCA 對感測嵌入流形進行剛體變換，使其與參考嵌入流形對齊。在檢查右側規則符合性后，可表示為：

其中，cs和cr分別為的感測點和參考點。

圖7 給出了流形對齊完成預期強度轉換的必要性。不同模態使用不同的相對強度水平來顯示大腦的同一部分。T1-MRI 和T2-MRI 圖像中的腦脊液、白質和灰質表現不同。通過流形學習，還沒有達到必要的強度調整。然而，流形學習和流形對齊在兩幅圖像上執行強度變換，使得它們使用可比較的強度級別。

圖7 T1-MRI特征圖像1和T2-MRI強度變換特征圖像1之間的顏色相似性

圖7 將多模態掃描轉換為相同強度的坐標系，雖然原始的T1-MRI 和T2-MRI 存在強度變化，但T2-MRI 強度轉換后的特征圖像與T1-MRI 特征圖像使用可比較的相對強度水平，在圖像中產生必要的對比度。

2.4 圖像配準

該文使用流形學習和流形對齊的多模態到單模態轉換提供了多模態掃描的結構表示。這種新的圖像表示代替了原始圖像，適用于使用單模態配準技術尋找配準參數。然后，使用估計的變換矩陣對原始的多模態掃描進行對齊。該文分別使用基于強度的圖像配準和基于FMT[5]的圖像配準來求解具有全部和部分數據的圖像配準問題。基于強度的配準迭代利用規則步長梯度下降優化器來調整變換參數，使優化沿著均方梯度的極值方向進行。

3 實驗分析

為了驗證該文提出的強度變換方法的性能，使用Laplacian 特征映射和提出的流形對齊技術，在完全或部分數據可用的情況下，來改善和促進多模態醫學掃描配準。

3.1 實驗裝置

該文實驗在Intel(R) Xeon(R) E3-1245 CPU，3.50 GHz，32 Gb 內存的個人計算機上運行Matlab R2013a 環境進行。為了驗證所提出方法的有效性，利用Elastix toolbox 進行基于NMI[9]的配準。在實驗中，NMI 和MI 度量使用50 個具有直方圖的所有像素進行計算。

對BrainWeb 數據庫中的T1、T2 和PD 加權MR圖像進行了實驗。實驗中的圖像含有3%的噪聲和20%的強度不均勻性。同時，還使用了CT 和Vanderbilt 數據庫中的不同MR 模式，該數據庫從RIRE 數據庫進行收集。此外，RIRE 數據庫提供了利用對靜場不均勻性進行校正而生成的校正圖像。所有圖像的大小都調整為200×200，以便節省處理時間和內存使用。

3.2 多模態到單模態轉換

為了直觀地驗證所提出的強度變換性能，使用與參考特征圖像中相似的相對強度級別來重新繪制感測特征圖像。使用來自不同數據集的圖像來檢驗該方法的有效性，該文分別研究了每幅圖像的流形，使用大小為3×3 的塊得到了9 維空間中包含40 000個點的流形，構造了k=20 近鄰的Laplacian 圖，并將流形嵌入到三維空間中。以T1 加權MRI 為參考，將其余部分與之對齊。多模態到單模態轉換的結果如圖8 所示。其中，第一行為人腦的原始MRI 和CT 掃描，第二行流形學習后每種模式的第一幅特征圖像；第三行流形對齊后的相同特征圖像。T2、PD-MRI和CT 的特征圖像的強度映射與T1-MRI 的特征圖像配準。在經過多種流形對齊之后，所有模態都使用了與T1-MRI 特征圖像1 類似的強度映射。

圖8 多模態醫學掃描的強度變換

3.3 全數據的多模態配準

為了驗證該文方法在全數據多模態圖像配準中的有效性，利用RIRE 數據集進行了兩組實驗。由于目前的方法只適用于二維圖像，因此獲得了切片對齊。首先，使用相同的參數D=9，k=10 和d=3 分別對每幅圖像的高維流形進行研究。然后，采用流形對齊技術將多模態掃描轉換為單模態強度坐標系。其次，利用各模態的第一幅特征圖像估計剛體配準參數。為此，該文采用了正則梯度下降優化器和均方誤差度量作為單模態配準的方法。最后，在原始圖像集上引入了配準參數。

使用所提出的方法檢查了用多種MR 模式配準CT 掃描的性能。由于沒有RIRE 測試數據集的真實信息，還需測量配準前后圖像對之間的相似性度量MI。MI量化了參考圖像和感測圖像之間的統計相關性程度。當前數據集包含10 名患者的CT、T1、T2 和PD-MRI掃描。其中，只有6名患者有經過校正的MRI。在統計誤差p＜0.5%內，每對CT-MR 掃描配準前后測量的MI 的均值和標準差（單位：mm），如表1 所示。

表1 每對CT-MR掃描配準前后測量的MI均值和標準差

RIRE 數據庫中CT 和PD-MRI 掃描樣本的對比如圖9 所示。結果表明，即使在配準過程不進行MI值的計算和優化，該方法仍能提高MI 值。

圖9 RIRE數據庫中CT和PD-MRI掃描樣本的對比

為了進一步定量地評價該文方法的性能，利用RIRE 訓練數據集，在[-π/4,π/4]范圍內隨機生成旋轉角度，并沿每個軸進行平移以此約束大腦區域保持在坐標系內，進而創建了一個CT（感測圖像）的非對齊切片。采用3 種不同的方法對失真圖像進行配準：1）文獻[11]中提出的基于Laplacian 圖的配準方法；2）文獻[8]中提出的MI 度量的多模態配準方法；3）該文方法。利用5 個隨機重定位點的平均絕對誤差(MAE)計算配準誤差。由于RIRE 訓練數據集只有患者的多模態圖像，該文對每對模態重復實驗30 次。距離誤差的平均值和標準差結果（單位：mm），如表2所示。

表2 每對CT-MR模態掃描的距離誤差的平均值和標準差

在表2 中，粗體數字表示每個模態對獲得的最佳結果。除帶星號數據外，其余結果均在統計誤差p＜0.1%內。CT-T1 校正掃描的配準誤差均在統計誤差p＜0.3%內。CT-T2 校正掃描的配準誤差高于MI，但無統計學意義。

3.4 部分數據的多模態配準

在多模態醫學掃描與部分數據的配準問題中，從不同模態下獲取圖像時，感測圖像（部分數據）只覆蓋部分參考圖像（全部數據），因此，需要在較大數據空間內定位小數據。

利用隨機參數對遙感圖像進行旋轉、平移和裁剪，生成局部數據的合成實例。在實驗中，在感測圖像中加入了額外的放縮。使用所提出的方法對兩幅圖像進行研究，然后將每幅圖像的特征圖像傳遞給FMT 算法。FMT 算法將相應地計算轉換參數(平移、旋轉和放縮)，利用估計的變換參數對原始的全圖像和局部圖像進行聯合配準。原始圖像和配準圖像對比結果如圖10 所示。

圖10 部分數據的多模態圖像配準對比

將患者曾經掃描的T1-MRI 與覆蓋患者大腦部分的新捕獲的T2-MRI（感測）進行配準。兩幅圖像均在軸向采集，在進行新的掃描（T2）時，患者的成像位置與T1 不同。參考圖像和感測圖像以及多種方法的配準結果，如圖11 所示。

圖11 感測圖像中的部分數據對合成實例進行配準

4 結束語

該文將所提出的多模態到單模態的變換作為預處理步驟，而不必考慮圖像的重疊，這有助于恢復平移、旋轉和放縮。在醫學圖像配準策略的背景下，該方法可作為參數化方法，但不需要對參數進行精確的匹配和賦值。為進一步調整高維流形的光滑性提供了自由度，最終將對特征值問題的收斂速度產生積極的影響。為了實現精確的結構表示，利用所提出的多模態到單模態變換，生成多個特征圖像作為結構表示的輸出。雖然該文只考慮了人腦圖像的剛性變換（除了考慮仿射變換的部分數據的配準外）來評估所提出的變換性能，但對其他圖像的配準以及非剛性配準都可以借助此方法。