李正偉

（成都理工大學(xué)工程技術(shù)學(xué)院，四川樂山 614007）

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

單模態(tài)醫(yī)學(xué)圖像配準(zhǔn)方法的目的是對(duì)同一成像技術(shù)(如CT)得到的圖像進(jìn)行配準(zhǔn)。文獻(xiàn)[4]提出了基于快速傅里葉變換(FFT)技術(shù)，該技術(shù)可尋求參考圖像和感測(cè)圖像之間的最佳平移、旋轉(zhuǎn)和縮放參數(shù)。這種相位相關(guān)技術(shù)的擴(kuò)展后來發(fā)展為基于Fourier-Mellin 變換(FMT)的圖像配準(zhǔn)方法[5]。文獻(xiàn)[6]提出了改進(jìn)的FMT 方法，消除了直角坐標(biāo)到對(duì)數(shù)極坐標(biāo)的轉(zhuǎn)換。評(píng)估并驗(yàn)證了該方法對(duì)噪聲的準(zhǔn)確性和魯棒性。

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

考慮到CT和MRI是兩種廣泛使用的成像方式，可以產(chǎn)生具有可比空間分辨率的橫截面圖像，因此MRI具有更高的對(duì)比度分辨率，而CT掃描則相對(duì)較低。常用的多模態(tài)配準(zhǔn)方法往往由于缺失部分?jǐn)?shù)據(jù)導(dǎo)致配準(zhǔn)無效，降維算法（即流形學(xué)習(xí)）常用于圖像和數(shù)據(jù)處理領(lǐng)域，其目的是研究嵌入高維空間的低維流形[10]。

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

1 數(shù)據(jù)降維

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

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

1.1 主成分分析(PCA)

PCA 的目的是通過正交變換將一組可能存在相關(guān)性的變量轉(zhuǎn)換為一組線性不相關(guān)的變量[12]。在PCA 中，假設(shè)數(shù)據(jù)集的方差是對(duì)數(shù)據(jù)集中存儲(chǔ)信息量的度量，則方差越大，存儲(chǔ)的信息越多。

給定矩陣X∈RD×n，其列為D維數(shù)據(jù)點(diǎn)，則矩陣Xc∈RD×n定義為中心數(shù)據(jù)點(diǎn)的矩陣（每行的樣本均值為零）。投影由Y=AXc給出，其中A為d×n矩陣。PCA 優(yōu)化目標(biāo)函數(shù)為：

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

利用Lagrange 乘數(shù)技術(shù)，將最大化問題簡(jiǎn)化為求解協(xié)方差矩陣V的特征值問題：

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

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

圖1 樣本合成數(shù)據(jù)集的主成分

1.2 Laplacian特征映射

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

Laplacian 特征映射基于Laplacian 圖，其算法包括3 個(gè)主要步驟[11]：1）構(gòu)造鄰接圖；2）選擇權(quán)值；3）尋找特征映射。每個(gè)步驟都涉及選擇適當(dāng)?shù)膶?shí)現(xiàn)方式并為多個(gè)參數(shù)確定適當(dāng)?shù)闹?。由Laplacian 特征映射所保證的保留局部性和結(jié)構(gòu)等價(jià)性是該文方法的關(guān)鍵決定因素。

2 多模態(tài)圖像配準(zhǔn)

該文提出的多模態(tài)到單模態(tài)變換的目的是將多模態(tài)輸入圖像變換成相似強(qiáng)度的坐標(biāo)系。通過提取比單個(gè)像素強(qiáng)度更多的信息來研究小圖像塊的結(jié)構(gòu)。由于來自不同模態(tài)的圖像的內(nèi)部結(jié)構(gòu)相似，因此，轉(zhuǎn)換將導(dǎo)致多模態(tài)圖像的相似結(jié)構(gòu)表示。多模態(tài)到單模態(tài)變換的圖像配準(zhǔn)如圖2 所示。

圖2 多模態(tài)到單模態(tài)變換的圖像配準(zhǔn)

2.1 構(gòu)建高維空間

在圖像配準(zhǔn)問題中，只存在參考和感測(cè)兩種圖像。在高維空間中，通過從每個(gè)輸入圖像中提取小塊來構(gòu)造點(diǎn)云。高維空間的每個(gè)維度專用于每個(gè)補(bǔ)丁的一個(gè)圖像像素?？紤]大小為s×s的每個(gè)補(bǔ)丁都比前一個(gè)補(bǔ)丁移動(dòng)了一個(gè)像素，對(duì)于大小nr×nc的圖像，高維流形在s2=D維空間中包括nr×nc=N個(gè)點(diǎn)。使用點(diǎn)的坐標(biāo)以矩陣格式表示這些點(diǎn)集，對(duì)于大小為200×200、塊大小為3×3的圖像，矩陣表示為40 000×9的矩陣，使用圖像的小塊和相應(yīng)的點(diǎn)云矩陣表示構(gòu)建高維點(diǎn)云的示例，如圖3 所示。對(duì)于大小為N像素、塊為D個(gè)像素的灰度圖像，矩陣大小為N×D。

圖3 高維空間中點(diǎn)云的矩陣表示

D的選擇取決于圖像的應(yīng)用和大小。當(dāng)選擇相對(duì)較小的D時(shí)，嵌入的流形類似于原始圖像。當(dāng)D選擇較大的值會(huì)導(dǎo)致圖像模糊。不同D值時(shí)T2-MRI圖像切片及其嵌入流形的圖像示例，如圖4所示。

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

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

2.2 流形學(xué)習(xí)

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

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

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

其中，σ∈R為唯一的參數(shù)。如果兩個(gè)節(jié)點(diǎn)i和j通過一條邊連接，則wij=1，否則wij=0。

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

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

所有邊的權(quán)重都限制在e-1/2≈0.606 和1 之間。在減少計(jì)算時(shí)間的同時(shí)，嵌入流形和配準(zhǔn)結(jié)果均有效且可靠。

最優(yōu)嵌入在保留局部信息的同時(shí)將點(diǎn)xi映射到點(diǎn)yi，即連接點(diǎn)xi和xj在分別映射到y(tǒng)i和yj之后將保持盡可能靠近。最小化以下目標(biāo)函數(shù)：

對(duì)于給定的具有wij項(xiàng)的對(duì)稱權(quán)值矩陣W，該文定義了對(duì)角度量矩陣D，使得。度量矩陣D的項(xiàng)是W的列（或行）之和。該文還將Laplacian 矩陣定義為L(zhǎng)=D-W。如果集合V包含N個(gè)點(diǎn)，則W以及D和L均為稀疏N×N對(duì)稱矩陣。

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

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

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

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

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

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

圖5 用Laplacian特征圖表示T1-MRI的結(jié)構(gòu)

算法1：Laplacian 特征映射降維（結(jié)構(gòu)表示）

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

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

步驟1：計(jì)算每?jī)蓚€(gè)數(shù)據(jù)點(diǎn)之間的距離；

步驟2：選擇k和對(duì)稱kNN 圖，構(gòu)造k近鄰構(gòu)造鄰接圖；

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

步驟4：采用熱核加權(quán)方案為每?jī)蓚€(gè)鄰居點(diǎn)之間的每個(gè)邊分配權(quán)重；

步驟5：構(gòu)造稀疏實(shí)對(duì)稱N×N矩陣W,D,L；

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

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

步驟8：對(duì)特征向量進(jìn)行排序，以遞增順序表示特征值；

步驟9：去掉第一個(gè)nc個(gè)特征向量；

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

2.3 流形對(duì)齊

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

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

圖6 兩個(gè)嵌入流形的右手規(guī)則檢查

該文提出了一種快速而直接的方法來檢查參考流形和感測(cè)流形是否符合右手規(guī)則，每?jī)煞鶎?duì)應(yīng)的特征圖像所對(duì)應(yīng)的協(xié)方差的符號(hào)是它們關(guān)聯(lián)方向的度量。如果它們有負(fù)關(guān)聯(lián)，則將-1 與感知嵌入流形的特征圖像相乘。當(dāng)使用該文方法補(bǔ)償反射誤差時(shí)，PCA計(jì)算剛體變換矩陣R，其對(duì)齊中心嵌入流形為：

其中，矩陣A和矩陣B分別為感測(cè)流形和參考流形的主方向。由于矩陣A為正交矩陣，根據(jù)PCA的性質(zhì)，旋轉(zhuǎn)矩陣簡(jiǎn)化為：

通過將兩個(gè)流形的中心移動(dòng)到原點(diǎn)，然后旋轉(zhuǎn)感測(cè)流形的主方向以對(duì)準(zhǔn)參考流形的主方向來實(shí)現(xiàn)對(duì)齊。利用PCA 對(duì)感測(cè)嵌入流形進(jìn)行剛體變換，使其與參考嵌入流形對(duì)齊。在檢查右側(cè)規(guī)則符合性后，可表示為：

其中，cs和cr分別為的感測(cè)點(diǎn)和參考點(diǎn)。

圖7 給出了流形對(duì)齊完成預(yù)期強(qiáng)度轉(zhuǎn)換的必要性。不同模態(tài)使用不同的相對(duì)強(qiáng)度水平來顯示大腦的同一部分。T1-MRI 和T2-MRI 圖像中的腦脊液、白質(zhì)和灰質(zhì)表現(xiàn)不同。通過流形學(xué)習(xí)，還沒有達(dá)到必要的強(qiáng)度調(diào)整。然而，流形學(xué)習(xí)和流形對(duì)齊在兩幅圖像上執(zhí)行強(qiáng)度變換，使得它們使用可比較的強(qiáng)度級(jí)別。

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

圖7 將多模態(tài)掃描轉(zhuǎn)換為相同強(qiáng)度的坐標(biāo)系，雖然原始的T1-MRI 和T2-MRI 存在強(qiáng)度變化，但T2-MRI 強(qiáng)度轉(zhuǎn)換后的特征圖像與T1-MRI 特征圖像使用可比較的相對(duì)強(qiáng)度水平，在圖像中產(chǎn)生必要的對(duì)比度。

2.4 圖像配準(zhǔn)

該文使用流形學(xué)習(xí)和流形對(duì)齊的多模態(tài)到單模態(tài)轉(zhuǎn)換提供了多模態(tài)掃描的結(jié)構(gòu)表示。這種新的圖像表示代替了原始圖像，適用于使用單模態(tài)配準(zhǔn)技術(shù)尋找配準(zhǔn)參數(shù)。然后，使用估計(jì)的變換矩陣對(duì)原始的多模態(tài)掃描進(jìn)行對(duì)齊。該文分別使用基于強(qiáng)度的圖像配準(zhǔn)和基于FMT[5]的圖像配準(zhǔn)來求解具有全部和部分?jǐn)?shù)據(jù)的圖像配準(zhǔn)問題?；趶?qiáng)度的配準(zhǔn)迭代利用規(guī)則步長(zhǎng)梯度下降優(yōu)化器來調(diào)整變換參數(shù)，使優(yōu)化沿著均方梯度的極值方向進(jìn)行。

3 實(shí)驗(yàn)分析

為了驗(yàn)證該文提出的強(qiáng)度變換方法的性能，使用Laplacian 特征映射和提出的流形對(duì)齊技術(shù)，在完全或部分?jǐn)?shù)據(jù)可用的情況下，來改善和促進(jìn)多模態(tài)醫(yī)學(xué)掃描配準(zhǔn)。

3.1 實(shí)驗(yàn)裝置

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

對(duì)BrainWeb 數(shù)據(jù)庫中的T1、T2 和PD 加權(quán)MR圖像進(jìn)行了實(shí)驗(yàn)。實(shí)驗(yàn)中的圖像含有3%的噪聲和20%的強(qiáng)度不均勻性。同時(shí)，還使用了CT 和Vanderbilt 數(shù)據(jù)庫中的不同MR 模式，該數(shù)據(jù)庫從RIRE 數(shù)據(jù)庫進(jìn)行收集。此外，RIRE 數(shù)據(jù)庫提供了利用對(duì)靜場(chǎng)不均勻性進(jìn)行校正而生成的校正圖像。所有圖像的大小都調(diào)整為200×200，以便節(jié)省處理時(shí)間和內(nèi)存使用。

3.2 多模態(tài)到單模態(tài)轉(zhuǎn)換

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

圖8 多模態(tài)醫(yī)學(xué)掃描的強(qiáng)度變換

3.3 全數(shù)據(jù)的多模態(tài)配準(zhǔn)

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

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

表1 每對(duì)CT-MR掃描配準(zhǔn)前后測(cè)量的MI均值和標(biāo)準(zhǔn)差

RIRE 數(shù)據(jù)庫中CT 和PD-MRI 掃描樣本的對(duì)比如圖9 所示。結(jié)果表明，即使在配準(zhǔn)過程不進(jìn)行MI值的計(jì)算和優(yōu)化，該方法仍能提高M(jìn)I 值。

圖9 RIRE數(shù)據(jù)庫中CT和PD-MRI掃描樣本的對(duì)比

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

表2 每對(duì)CT-MR模態(tài)掃描的距離誤差的平均值和標(biāo)準(zhǔn)差

在表2 中，粗體數(shù)字表示每個(gè)模態(tài)對(duì)獲得的最佳結(jié)果。除帶星號(hào)數(shù)據(jù)外，其余結(jié)果均在統(tǒng)計(jì)誤差p＜0.1%內(nèi)。CT-T1 校正掃描的配準(zhǔn)誤差均在統(tǒng)計(jì)誤差p＜0.3%內(nèi)。CT-T2 校正掃描的配準(zhǔn)誤差高于MI，但無統(tǒng)計(jì)學(xué)意義。

3.4 部分?jǐn)?shù)據(jù)的多模態(tài)配準(zhǔn)

在多模態(tài)醫(yī)學(xué)掃描與部分?jǐn)?shù)據(jù)的配準(zhǔn)問題中，從不同模態(tài)下獲取圖像時(shí)，感測(cè)圖像（部分?jǐn)?shù)據(jù)）只覆蓋部分參考圖像（全部數(shù)據(jù)），因此，需要在較大數(shù)據(jù)空間內(nèi)定位小數(shù)據(jù)。

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

圖10 部分?jǐn)?shù)據(jù)的多模態(tài)圖像配準(zhǔn)對(duì)比

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

圖11 感測(cè)圖像中的部分?jǐn)?shù)據(jù)對(duì)合成實(shí)例進(jìn)行配準(zhǔn)

4 結(jié)束語

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