基于CT和MRI的多模態(tài)醫(yī)學圖像配準算法

常錦才，張笑林，武偉娜，蘇照，王雨暄，王嘉祺

( 1.華北理工大學 理學院，河北 唐山 063210；2.華北理工大學 3D建模與應用創(chuàng)新實驗室，河北 唐山 063210；3.華北理工大學 基礎醫(yī)學院，河北 唐山 063210)

引言

圖像配準提出的最初目的是為了對2幅或多幅圖像進行圖像融合，這些圖像大多數(shù)來自于同一種成像原理下的同一物體，主要是對不同視角下的圖像進行信息配準和拼接。隨著科學技術的不斷發(fā)展，各類傳感器層出不窮，人們對圖像配準的期望也不單單滿足于簡單的單模態(tài)和拼接功能，在大量應用場景和需求的催生下，針對同一物體，對不同成像原理而拍攝的圖像進行配準，整合不同模態(tài)下圖像的各類特征。多模態(tài)下的圖像配準被廣泛應用于多個領域，其中與醫(yī)學影像學相結合的圖像融合技術，已經(jīng)成為了圖像配準研究中的主流。

電子計算機斷層掃描(CT)和磁共振功能成像(MRI)是目前的醫(yī)學影像領域最為成熟，而且應用最廣的成像方式。CT和MRI都是以DICOM格式進行存儲的，這對圖像融合和配準帶來了便利。針對顱腦而言，CT成像可用于腦部腫瘤和血管的三維重建，MRI則主要是對腦內的神經(jīng)纖維進行三維重建。二者圖像的配準融合為醫(yī)學工作者和研究人員提供了極大的便利。

1圖像配準的基本框架

圖像配準的基本框架包括以下4個方面：特征空間、空間變換、最優(yōu)化搜索策略和相似性度量[1]。

特征空間的選取是圖像配準工作的重點，一般也稱為特征提取。通常情況下，特征空間的選擇隨配準算法而變化，不同的算法將具有不同的特征空間。比如，當配準算法以灰度值為基準，特征空間就可以選為圖像中像素的灰度值；而大多數(shù)配準所采用的圖像都是來自不同的圖像設備，以灰度值為標準并不適用，這類算法在配準時可以選用圖像空間位置上比較相近的點、邊緣或者是曲線、曲面等。特征空間不僅關乎于算法的運行效率和結果，還影響算法的穩(wěn)定性。大量的研究表明，理想的特征空間應該滿足：特征提取簡便快捷、特征匹配計算量小、特征數(shù)據(jù)量適宜，不受噪聲、光照度等因素影響，且對各種圖像均能使用。圖1所示為圖像配準框架。

圖1 圖像配準框架

空間變換是指配準過程中圖像變換的范圍和方式，可以劃分為全局變換、局部變換和位移場變換。其中全局變換的變換域最廣，可以覆蓋到整個圖像上，但其變換效果往往不夠精準。局部變換對全局變換進行了改進，可以對圖像的關鍵點部位進行參數(shù)上的調整，而關鍵點以外的圖像則使用插值算法進行處理。位移場變換精度最高，可以實現(xiàn)圖像上每個像素點獨立的參數(shù)變換，而這些變換則用一個巧妙的連續(xù)函數(shù)串聯(lián)起來，從而實現(xiàn)整幅圖像的變換。圖像的變換方式基于其數(shù)學原理，線性變換一般指矩陣變換，主要是對圖像進行旋轉、平移、縮放等變換，能對圖像實現(xiàn)基本的剛性配準。非線性變換的主要方式為多項式函數(shù)，能夠實現(xiàn)圖像的形變。在實際應用中，合理地使用空間變換能讓配準算法事半功倍。

最優(yōu)化搜索與相似性度量在功能上相輔相成，搜索最優(yōu)的配準參數(shù)需要以2幅圖像的相似度度量的數(shù)值為判斷依據(jù)，它也直接反映了最優(yōu)搜索策略的好與壞。實際的圖像配準算法流程看似簡單，但其背后卻包含了大量的運算過程，選擇一個高效快捷的尋優(yōu)策略是降低運算的有效途徑，最優(yōu)化搜索方法種類繁多，但在進行尋優(yōu)設計時，應該把眼光放到更快、更高效的算法上，理論數(shù)學由傳統(tǒng)的貪婪算法出發(fā)，已經(jīng)發(fā)展出多種最優(yōu)化搜索算法，包括一維搜索中的試探法、函數(shù)逼近法；無約束最優(yōu)化的梯度方法中的最速下降法、牛頓法、最小二乘法等；無約束最優(yōu)化的直接方法中的模式搜索法、Powell法等。

相似性度量是圖像配準過程中用來評判配準結果好壞的準則，相似度度量的數(shù)值直接反映了特征空間與空間變換算法的優(yōu)劣，同時它也承上啟下，它的數(shù)值也會反饋到搜索策略上，對搜索策略的選取具有重要意義。相似性的數(shù)學表達式如下所示：

(1)

公式中的xi，yi分別表示參考與浮動2幅圖像的第i個像素的灰度值；xm、ym參考與浮動2幅圖像的平均強度。針對不同的圖像，配準算法也會有相似性度量。比如相關性是一種適用于單模態(tài)圖像配準的相似性度量，它可以幫助研究人員了解浮動圖像和參考圖像之間的相似程度。相關性的數(shù)值直接反映了圖像的相似程度，若數(shù)值為1，則表示2幅圖像完全一致；若相關性的數(shù)值為0，則說明2幅圖像完全不一致，圖像可能來自不同物體，或相似性算法的選擇出現(xiàn)問題。

2多模態(tài)醫(yī)學圖像配準常用算法

2.1 最大互信息配準算法

最大互信息配準方法與傳統(tǒng)相似性度量之間的區(qū)別是，前者無視噪點的存在，也不必再對圖像進行分割。互信息不需要對參考與浮動圖像直接的強度關系進行初步假設，所以最大互信息配準方法被廣泛應用于多模態(tài)圖像配準工作，其強大的特性為配準工作帶來了諸多便利，不僅適用于多模態(tài)圖像的配準，面對破損的圖像，互信息也處理得游刃有余。

互信息用于度量一個變量包含另一個變量的信息量或者2個隨機變量的統(tǒng)計依賴性，它是信息論的基本概念，2個隨機變量A、B間的互信息可根據(jù)信息論表示如下[3]：

(2)

式中，pAB(a,b)表示隨機變量A、B的聯(lián)合概率分布，pA(a)和pB(b)分別表示A、B的邊緣概率分布。

同時，根據(jù)Shannon信息熵的定義：

H(A)=-∑pA(a)log2pA(a)

(3)

H(B)=-∑pB(b)log2pB(b)

(4)

(5)

(6)

(7)

上述公式中，H(A,B)表示條件A、B的聯(lián)合熵。H(A)、H(B)分別表示A和B的信息熵，H(A|B)表示給定B的條件下A的條件熵，H(B|A)表示給定A的條件下B的條件熵，按照琴生不等式的概念，可以證明I(A,B)是非負的，則信息熵可以將互信息表示如下：

I(A,B)=H(A)-H(A∣B)

=H(A)-H(B∣A)

=H(A)+H(B)-H(A,B)

(8)

圖2以維恩圖的形式展現(xiàn)了邊緣熵、條件熵和聯(lián)合熵之間交并差集的關系。

圖2 互信息的圖形概述

綜上所述，引入最大互信息的概念，應用到配準算法中，可以將既定的互信息測度定義為圖像配準結果的判斷依據(jù)，即配準圖像之間的互信息達到最大值。此時作為隨機向量的兩圖像系統(tǒng)的熵達到最小，也就是說2幅圖像的不確定性己經(jīng)很小了[4]。互信息的應用不僅縮減了配準所耗費的時間，也提高了算法的性能，因此，互信息被應用在許多圖像配準過程中。

從概率的角度出發(fā)，當pAB(a,b)=pA(a)pB(b)，且I(A,B)=0時，條件A和條件B相互獨立。當pAB(a,b)=pA(a)=pB(b)時，則稱條件A、B完全依賴或者條件A、B完全包含，即H(A)=H(B)=H(A,B)，此時I(A,B)最大。互信息的基本性質如表1所示。

表1 互信息的基本性質

2.2 基于Canny算子和小波變換的最大互信息算法

John F. Canny于上世紀80年代提出了基于一種理論的最優(yōu)邊緣檢測算子，稱為Canny算子，雖然其自提出至今已有30多年，但卻不影響其成為圖像配準領域之中應用最廣的邊緣檢測算法之一。Canny算法目標明確，首先要保證特征提取過程中的低錯誤率，其次要更好、更快地對位圖像中的邊緣點。由此引出了Canny算法的4個主要步驟：低通濾波降噪、有限差分法計算梯度幅值和方向、對像素點進行非極大值抑制、雙閾值算法檢測和連接邊緣。

在數(shù)學領域，小波變換的提出為許多問題的解決提供了新方向，它也是傅里葉變換和短時傅里葉變換的一個重大突破，其優(yōu)越的特性為圖像處理領域開拓了新思路。小波變換的性能也優(yōu)于傳統(tǒng)的圖像處理方法，其利用尺度函數(shù)和平移函數(shù)對函數(shù)或信號進行細化和分析，解決了許多傅里葉變換無法解決的難題。小波變換也因此得到了大量應用，尤其適用于近些年熱門的圖像處理研究。

由于小波變換的抗噪性很強，所以當Canny算子在圖像邊緣提取出來的是噪點時，小波變換無法在對應邊緣提取出同樣的噪點。因此，可以通過去除Canny算子提取的邊緣圖像中的噪聲點，得到圖像的邊緣[5]。圖3所示為邊緣提取算法的流程圖。

圖3 邊緣提取算法流程圖

2.3 基于輪廓金字塔的歸一化互信息圖像配準算法改進

對于灰度值分布不均勻，甚至分布比較極端的圖像，傳統(tǒng)算法一般很難處理。由于小波變換引入了多尺度的概念，可以利用不同尺度的圖像進行多尺度分析，從而提高了配準的精度和效率，使圖像得到良好的配準處理。

輪廓金字塔是小波變換在圖像處理領域針對序列圖像而提出的一種多尺度策略，通過對圖像進行模糊處理或下采樣，達到對序列圖像進行配準的目的，下采樣的圖片由下至上，層層遞減，堆疊起來的序列圖像就如圖金字塔一般，故而成為輪廓金字塔算法。使用CT和MRI圖像進行一次多模態(tài)醫(yī)學圖像配準，并給出輪廓提取和下采樣金字塔的圖像依據(jù)和理論概述，從圖像的角度對算法進行一個粗略的剖析。Canny算子具有優(yōu)質的邊緣檢測能力，可以對CT和MRI分別進行輪廓提取，下采樣金字塔自下而上隔行隔列地降低分辨率，使算法所需算力呈指數(shù)級下降。最后給出如圖4所示配準前后的對比圖。

圖4 配準前后的對比圖

圖4流程圖左上角的圖像為CT數(shù)據(jù)，左下角為MRI數(shù)據(jù)，從左到右分別對圖像進行了預處理、輪廓提取、下采樣和圖像配準，右側是配準結果對比

3基于3D Slicer平臺的圖像配準

3.1 圖像配準模塊簡介

3D Slicer為學術研究人員和醫(yī)療從業(yè)者提供了眾多關于針對圖像配準和圖像融合的模塊，這些模塊包含大量算法和算法測評規(guī)則。

3.1.1 General Registration(BRAINS)

General Registration BRAINS(以下簡述BRAINS)是一種基于ITK的顱腦圖像互信息配準程序，算法過程完全實現(xiàn)自動化[6]。該模塊的開發(fā)源自BRAINSFit項目，其初衷是為多模態(tài)的三維圖像進行配準。

對比ITK內的算法，BRAINS模塊功能全面，不僅包含各種類型的變換方式，且能對內部函數(shù)進行優(yōu)化。圖5為BRAINS模塊的內部選項欄。在配準前輸入?yún)⒖紙D像和浮動圖像，在Percentage of Samples一欄可以調節(jié)配準的精度，默認為0.002，可以設置參數(shù)為0.2(20%)以獲得高細節(jié)。

圖5 General Registration功能展示

3.1.2 General Registration(Elastix)

General Registration(Elastix)模塊依托于Elastix軟件，Elastix是一款基于ITK的開源軟件。該軟件由一組通常用于解決醫(yī)學圖像配準問題的算法組成。Elastix的模塊化設計允許用戶快速配置、測試和比較特定應用程序的不同注冊方法，命令行界面通過腳本支持對大量數(shù)據(jù)集的自動化處理[7]。現(xiàn)在Elastix有了SimpleElastix，使得它可以在C++、Python、Java、R、Ruby、C#和Lua等語言中使用。

3D Slicer中的Elastix去繁從簡，只保留了全自動的圖像配準功能，去除了類似于BRAINS模塊中復雜的參數(shù)設置，降低了使用門檻，配準過程迅速快捷。圖6所示為配準前后2幅圖像的對比結果，由圖9可以看出，Elastix采用的是非線性變換，且對配準圖像進行了圖像融合。

圖6 Elastix模塊配準結果前后對比圖

3.1.3 Landmark Registration

Landmark Registration的主要功能是實現(xiàn)參考圖像與浮動圖像之間的手動配準，從而得到圖像之間的變換參數(shù)，該模塊對圖像類別沒有限制，適用于多模態(tài)圖像的配準。

Landmark Registration模塊的另一個優(yōu)點是具有實時性，打開模塊后，輸入?yún)⒖紙D像和浮動圖像，選擇變換類型，然后在模塊中添加基準點(Landmark)，通過調整圖像中的基準點使圖像重合。圖7所示為Landmark Registration示意圖。

圖7 Landmark Registration示意圖

3.1.4 CheckerBoard Filter

CheckerBoard Filter模塊可用于對比配準結果。它的功能是基于參考圖像與浮動圖像，創(chuàng)造出一個新的棋盤圖像，輸出圖像將顯示2個圖像之間的相似度，浮動圖像和配準圖像之間被重采樣到相同的原點、間距和方向。模塊的使用不受圖像類別的限制，適用于檢測多模態(tài)圖像配準的結果。圖8所示為CheckerBoard Filter選項欄展示。

圖8 CheckerBoard Filter選項欄展示

圖9為冠矢軸示意圖，CheckerBoard Parameters中Check Pattern的參數(shù)表示使用者可以從水平位、矢狀位和冠狀位3個維度對圖像進行切分重組。數(shù)字代表切分后的像素塊，1代表不切分。比如2，2，1是指在前2個維度進行切分，而第三個維度不進行切分。結合醫(yī)學上3個位面的定義圖與CheckerBoard的結果，圖10所示為CheckerBoard示意圖，可以清晰地看到模塊如何根據(jù)參數(shù)進行切分。

圖9 冠矢軸示意圖

圖10 CheckerBoard示意圖

3.2 圖像配準實例

圖11所示為配準流程圖。

圖11 配準流程圖

3.2.1配準圖像準備

圖像來源于軟件自帶的樣本數(shù)據(jù)，選用CT和MRI的圖像。首先要對數(shù)據(jù)進行預處理，使用Swiss Skull Stripper模塊去除顱腦以外的部分，只保留顱腦可以使圖像的輪廓更加清晰，這樣方便在后續(xù)的手動配準進行操作。圖12所示為MRI圖像去除顱骨前后的對比圖。

圖12 顱骨去除前后對比

3.2.2相似度檢測

使用CheckerBoard Filter模塊，對去除顱骨后的CT和MRI圖像進行相似度檢測，Check Pattern參數(shù)設為4、4、4。圖13所示為2幅圖像的原圖像和切分后重組的棋盤圖像的對比圖。這幅圖像較好地展示了模塊的切分效果和2幅圖像重采樣得到的原點、間距和方向都不盡相同。

圖13 圖像配準前相似度檢測

3.2.3配準實現(xiàn)

打開Landmark Registration模塊，輸入MRI圖像作為參考圖像，輸入CT圖像作為浮動圖像。點擊Apple，進入下一步操作。在空間變換中使用仿射配準(Affine Registration)，插入基準點(如圖14所示)后，手動移動基準點，可以實時地觀測到重疊部分的圖像隨基準點的移動而移動。參考4.2.3的模塊介紹，可以看到2幅圖像在3個維度上都實現(xiàn)了配準。

圖14 插入基準點后的變換結果

配準效果滿意之后，模塊會自動保存根據(jù)調整基準點得到的變換，再打開Transform模塊，對圖像進行變換，得到配準圖像。如果沒有對圖像變換這一步，則得到的只是變換，而不是圖像。

3.2.4配準結果檢測

得到配準后的圖像后，再使用CheckerBoard模塊進行相似度檢測，2次檢測的結果如15圖所示，由圖15可以看到,配準結果較好。至此，完成了基于CT和MRI的多模態(tài)圖像配準，也得到了配準后的浮動圖像。可以看到配準后的圖像輪廓清晰可見，配準效果較好。

圖15 配準前后的棋盤圖像展示

4結論

(1)二維層面的算法配準，其原理類似于最優(yōu)化原理中的尋優(yōu)過程，重點是特征空間和最優(yōu)搜索空間的選取。醫(yī)學圖像配準工作的開展是為了后續(xù)的圖像融合和對融合后的圖像進行三維重建，所以算法還需要更加快速高效。

(2)3D Slicer 是一個極其強大的醫(yī)學影像處理開源平臺，多種功能模塊的開發(fā)為醫(yī)學圖像配準提供了更多的可能性。軟件的熟練掌握也避免低效率的底層開發(fā)。

(3)探索完成了圖像配準和融合，以及融合后三維模型的重建。直觀的三維模型，擺脫了以往手術方案制定平面化的短板，讓患者的病情一目了然，這為術前的手術方案制定和術中的穿刺導航都提供了極大的便利。