李登旺 王洪君* 萬(wàn)洪林 尹 勇 孫偉峰（山東大學(xué)信息科學(xué)與工程學(xué)院，濟(jì)南 5000）

2（山東師范大學(xué)物理與電子科學(xué)學(xué)院，濟(jì)南 250014）

3（山東省腫瘤醫(yī)院，濟(jì)南 250117）

4（中國(guó)石油大學(xué)信息與控制工程學(xué)院，東營(yíng) 257061）

引言

近年來(lái)，醫(yī)學(xué)影像形變配準(zhǔn)及其應(yīng)用的研究正在受到醫(yī)學(xué)界和工程界的高度重視。如何提高形變配準(zhǔn)算法的精度和速度，是圖像引導(dǎo)放射治療（image guided radiation therapy）和自適應(yīng)放療（adaptive radiation therapy）領(lǐng)域的研究熱點(diǎn)。目前，多數(shù)配準(zhǔn)算法復(fù)雜，耗時(shí)巨大，有些還存在需要進(jìn)行手工標(biāo)記、不能實(shí)現(xiàn)自動(dòng)配準(zhǔn)等諸多問(wèn)題，束縛了它們?cè)谂R床上的應(yīng)用［1］。因此，研究能夠有效描述形變特點(diǎn)的模型，提高形變配準(zhǔn)的精度和速度，以達(dá)到臨床應(yīng)用的要求，是當(dāng)前醫(yī)學(xué)界和工程界一個(gè)富有挑戰(zhàn)性的課題，具有非常重要的理論意義和臨床價(jià)值。

例如，隨著放射治療技術(shù)的迅速發(fā)展，臨床上已采用了以三維適形放射治療（3D-CRT）和調(diào)強(qiáng)放療（IMRT）為代表的精確放療技術(shù)。目前，精確放療通常的做法是先CT定位，由醫(yī)生一次性制定放療計(jì)劃并依照此計(jì)劃實(shí)施逐次照射。然而，由于在治療過(guò)程中會(huì)不可避免的產(chǎn)生器官運(yùn)動(dòng)、變形和各種誤差（如擺位誤差等），這會(huì)使腫瘤（靶區(qū)）和危及器官偏離射野，從而導(dǎo)致腫瘤區(qū)域欠劑量照射和危及器官的過(guò)劑量照射。為解決這些問(wèn)題，科學(xué)家將放射治療機(jī)與成像設(shè)備結(jié)合在一起，在治療時(shí)采集有關(guān)的圖像信息，確定治療靶區(qū)和重要結(jié)構(gòu)的位置及相應(yīng)的運(yùn)動(dòng)軌跡，并在必要時(shí)進(jìn)行位置和劑量分布的校正，這稱為圖像引導(dǎo)放射治療［2］。

近年發(fā)展起來(lái)的錐形束CT（cone beam CT，CBCT）可以通過(guò)多個(gè)途徑確定和跟蹤靶區(qū)并引導(dǎo)放療，大大提高了IGRT的精度，成為目前IGRT開發(fā)和應(yīng)用的熱點(diǎn)。IGRT系統(tǒng)可以在直線加速器治療以前，通過(guò)kV級(jí)CBCT系統(tǒng)得到病人的CBCT三維重建圖像，與制定計(jì)劃時(shí)的Planning CT圖像（PCT）進(jìn)行配準(zhǔn)操作，得到CBCT重建圖像與計(jì)劃治療區(qū)的偏移情況，指導(dǎo)擺位信息［3］。由于該系統(tǒng)要求配準(zhǔn)速度要快，因此，目前系統(tǒng)中指導(dǎo)擺位的自動(dòng)配準(zhǔn)方法多是基于骨形特征的配準(zhǔn)［2，4］，其精度依賴于骨形特征提取的準(zhǔn)確性，但是，由于器官運(yùn)動(dòng)、變形以及治療過(guò)程等方面因素會(huì)使腫瘤相對(duì)骨形的位置發(fā)生改變，同時(shí)，對(duì)于盆腔內(nèi)部等不含骨組織的區(qū)域，無(wú)法進(jìn)行骨形的配準(zhǔn)。PCT與CBCT的形變配準(zhǔn)研究是當(dāng)前關(guān)于自適應(yīng)放療的熱點(diǎn)問(wèn)題，國(guó)際上多位學(xué)者已經(jīng)進(jìn)行了理論探索?；贐樣條的形變配準(zhǔn)方法［5］，基于FEM的形變配準(zhǔn)方法［6］，基于Demons的形變配準(zhǔn)方法［7］，基于光流學(xué)模型［8］的形變配準(zhǔn)方法分別被用于PCT與CBCT的形變配準(zhǔn)研究，用來(lái)指導(dǎo)基于CBCT的圖像引導(dǎo)放射系統(tǒng)中的自適應(yīng)放療。由于計(jì)劃CT與CBCT配準(zhǔn)時(shí)，全局形變和局部形變同時(shí)存在，而且對(duì)于不同病例的形變差別較大，盡管上述各種方法在恢復(fù)形變時(shí)分別具有各自的局部?jī)?yōu)勢(shì)，但用于自適應(yīng)放療恢復(fù)形變時(shí)，仍存在局限性，在精度、速度和魯棒性上仍需要改進(jìn)。

為了改進(jìn)基于CBCT的圖像引導(dǎo)放射系統(tǒng)中的形變配準(zhǔn)算法，提出一種基于正交小波變換的計(jì)劃CT與錐形束CT圖像的形變配準(zhǔn)方法。我們采用正交小波基函數(shù)由粗到精地表示形變域，低分辨率下的小波系數(shù)表示全局形變特征，高分辨率下的小波系數(shù)表示局部形變特征，這些小波系數(shù)由粗糙尺度到精細(xì)尺度逐步通過(guò)極小化形變能量函數(shù)估計(jì)得到，形變能量函數(shù)由Navier偏微分方程來(lái)設(shè)計(jì)。為了使自適應(yīng)放療的形變配準(zhǔn)算法具有較高的魯棒性，形變能量函數(shù)的外力部分采用歸一化互信息作為約束。仿真實(shí)驗(yàn)結(jié)果及臨床應(yīng)用均驗(yàn)證了本算法的有效性。

1 Navier偏微分方程

在形變配準(zhǔn)中，Navier偏微分方程可以作為描述發(fā)生形變的各向同性物質(zhì)的平衡狀態(tài)方程［9］?；诒“鍢訔l，B樣條等數(shù)學(xué)徑向基函數(shù)的配準(zhǔn)方法是通過(guò)插值和擬合理論來(lái)接近形變域，而基于Navier偏微分方程的方法可以通過(guò)內(nèi)力使物質(zhì)形成各種自由變形，并且通過(guò)約束力來(lái)保證恢復(fù)形變，可以得到形變域的數(shù)學(xué)物理解。因此Navier偏微分方程是描述非線性形變配準(zhǔn)的有效方法。Navier偏微分方程基于連續(xù)力學(xué)模型，模型中器官或者組織的形變是由內(nèi)力和外力共同作用而成，當(dāng)形變過(guò)程完成時(shí)內(nèi)力與外力達(dá)到平衡。達(dá)到平衡狀態(tài)后的Navier偏微分方程表達(dá)為

其中，θ是三維膨脹系數(shù)，表達(dá)為

式中，X=（x1，x2，x3）T是發(fā)生形變組織的三維空間坐標(biāo)系，F(xiàn)=（F1，F(xiàn)2，F(xiàn)3）T是作用于形變組織各處的外力。形變常量μ和λ取決于發(fā)生形變組織的自身性質(zhì)。u=（u1，u2，u3）T是配準(zhǔn)算法所需恢復(fù)參考對(duì)象與目標(biāo)對(duì)象之間的形變域。在發(fā)生形變組織或器官的邊緣，外力近似為零，對(duì)應(yīng)圖像的輪廓或者邊緣；而對(duì)于發(fā)生形變的其它部位，均需要外力和內(nèi)力的平衡。對(duì)于給定的形變域表示方法和受力，物體的形變域完全由外力和形變常數(shù)來(lái)決定。為了利用Navier偏微分方程來(lái)解決彈性形變問(wèn)題，需要適當(dāng)選擇外力和內(nèi)力模型，并且用合適的方法來(lái)表達(dá)形變域，這對(duì)于能否準(zhǔn)確恢復(fù)形變域具有決定作用。

平衡方程式（1）的前兩項(xiàng)，是造成組織或者器官形變的內(nèi)力。將式（2）帶入式（1），得到如下的內(nèi)力方程：

式中，λ和μ是Lame系數(shù)，它們的取值由彈性形變物質(zhì)的自身性質(zhì)決定。在不影響算法精度的前提下，為了簡(jiǎn)化計(jì)算的復(fù)雜度，選擇λ=0和μ=1。

2 利用正交小波基函數(shù)實(shí)現(xiàn)形變配準(zhǔn)

2.1 正交小波基函數(shù)表示形變域

臨床上涉及到的病例種類較多，因此形變特征較豐富，時(shí)常是全局形變與局部形變共存，同一組織或器官的不同部位大形變區(qū)域與小形變區(qū)域共存。對(duì)形變域進(jìn)行由粗及精描述，即對(duì)全局形變和局部形變分別描述，可以更加準(zhǔn)確可靠地恢復(fù)形變域。同時(shí)，由粗及精的策略可以使配準(zhǔn)過(guò)程避免陷入局部極值，提高配準(zhǔn)的魯棒性。傳統(tǒng)由粗及精策略利用塔式數(shù)據(jù)結(jié)構(gòu)［10］。小波（Wavelet）變換作為信號(hào)的一種多分辨率表示方法，自然適合對(duì)形變域由粗及精地描述［11］。

對(duì)于形變配準(zhǔn)問(wèn)題，將形變域表示為小波系數(shù)的函數(shù)，則參考圖像、目標(biāo)圖像、形變域有如下數(shù)學(xué)表達(dá)：

式中，（x，y，z）是參考圖像的空間坐標(biāo)系，（x′，y′，z′）是目標(biāo)圖像的空間坐標(biāo)系。形變域u=（u1，u2，u3）是小波系數(shù)c的函數(shù)，作為目標(biāo)參數(shù)，它可以由極小化能量函數(shù)得到。如前所述，將參考圖像和目標(biāo)圖像之間的形變域通過(guò)小波變換表示。由于正交小波變換具有能量集中特性，在表示形變特征時(shí)具有優(yōu)勢(shì)。因此，對(duì)于三維醫(yī)學(xué)圖像配準(zhǔn)，形變域u（X），在三維子空間N=（Nx，Ny，Nz）中，用三維離散正交小波分解表示可得

式（5）中各個(gè)參數(shù)的含義與普通的三維正交小波分解含義相同［12］。表示形變域的正交小波基函數(shù)在三維網(wǎng)格N=（Nx，Ny，Nz）中，以2j為空間間隔進(jìn)行平移變換。基函數(shù)采用的是一維小波基函數(shù)的三維張量積形式，并且由小波系數(shù)加權(quán)。把式（5）中對(duì)形變域參數(shù)的表示形式帶入式（3），可得到Navier偏微分方程內(nèi)力部分的線性組合形式。這一形式與各向同性的靜態(tài)Navier偏微分方程的Wavelet-Galerkin離散矩陣成線性正比，運(yùn)算上可以優(yōu)化［13］。

2.2 歸一化互信息形成外力約束

為了使Navier偏微分方程處于平衡狀態(tài)，從而恢復(fù)形變域，需為式（3）中的平衡方程設(shè)計(jì)合適的外力。外力的設(shè)計(jì)需要利用配準(zhǔn)數(shù)據(jù)對(duì)的特征。僅僅利用局部特征如曲面引導(dǎo)、對(duì)應(yīng)點(diǎn)特征引導(dǎo)或其它局部相似特征引導(dǎo)等，雖然可以獲得較快的配準(zhǔn)速度，但由于全局信息的丟失，恢復(fù)形變將會(huì)受限。因此，為了能夠?qū)θ中巫兒途植啃巫兺瑫r(shí)恢復(fù)，選取適當(dāng)?shù)奶卣鱽?lái)產(chǎn)生外力，成為是否可以準(zhǔn)確恢復(fù)形變域的關(guān)鍵所在。

在配準(zhǔn)過(guò)程中，外力的作用是使待配準(zhǔn)圖像數(shù)據(jù)對(duì)的相似或者相同的區(qū)域達(dá)到空間上的對(duì)應(yīng)。因此，可以利用參考圖像和目標(biāo)圖像的相似性測(cè)度來(lái)模擬外力?；隗w素相似性，特別是互信息的醫(yī)學(xué)圖像配準(zhǔn)方法精度高，魯棒性好，易實(shí)現(xiàn)全自動(dòng)配準(zhǔn)，并且可以實(shí)現(xiàn)多模態(tài)圖像的準(zhǔn)確配準(zhǔn)［14］。同時(shí)，互信息不僅包含圖像局部的信息，而且包含全局的信息，這對(duì)于恢復(fù)全局形變和局部形變具有重要意義。由于歸一化互信息與待配準(zhǔn)圖像數(shù)據(jù)對(duì)的重疊程度無(wú)關(guān)，相對(duì)于基于互信息的配準(zhǔn)，魯棒性更好，也是近年來(lái)醫(yī)學(xué)圖像配準(zhǔn)時(shí)普遍采用的相似性測(cè)度準(zhǔn)則?；跉w一化互信息的相似性測(cè)度定義為：

式中，H（X）是圖像X的香農(nóng)熵。因此，我們定義外力與參考圖像和目標(biāo)圖像的相似性函數(shù)S（u）成比例。為了估計(jì)目標(biāo)函數(shù)中的小波系數(shù)，需要對(duì)模擬平衡狀態(tài)的能量函數(shù)進(jìn)行極小化，因此，在這里，定義外力是歸一化互信息的倒數(shù)：

2.3 形變配準(zhǔn)算法

至此，得到了用于構(gòu)成Navier偏微分方程來(lái)解決彈性形變問(wèn)題的內(nèi)力和外力，其中形變域用正交小波基函數(shù)來(lái)表示，內(nèi)力和外力的平衡最終形成了關(guān)于形變域u的能量函數(shù)。因此，最終的能量函數(shù)可以表達(dá)為：

式中，w是加權(quán)常數(shù)，實(shí)驗(yàn)中，采用w為常數(shù)1。因此，形變域是關(guān)于小波系數(shù)的函數(shù)，即能量函數(shù)也是關(guān)于小波系數(shù)的函數(shù)。在配準(zhǔn)過(guò)程中，代表全局形變特征的的大尺度小波系數(shù)首先被估計(jì)得到，然后估計(jì)對(duì)應(yīng)局部形變特征的小尺度小波系數(shù)。為了減少需要被估計(jì)的小波參數(shù)的數(shù)量，根據(jù)子帶和尺度的不同我們把小波系數(shù)進(jìn)行分類，首先估計(jì)大尺度的小波系數(shù)，然后順序估計(jì)同一尺度上的其他子帶系數(shù)，最后估計(jì)小尺度上系數(shù)。從另一個(gè)角度看，兩幅圖像之間的形變由不同尺度和不同子帶的小波系數(shù)來(lái)描述。我們采用遞進(jìn)策略尋找這些小波系數(shù)，先估計(jì)大尺度上的小波系數(shù)，然后再估計(jì)小尺度上的小波系數(shù)。圖（1）是求解二維形變域時(shí)估計(jì)小波系數(shù)的順序。

形變域經(jīng)過(guò)小波變換以后，從理論上講，需要估計(jì)的小波系數(shù)和待配準(zhǔn)目標(biāo)圖像的體素個(gè)數(shù)是相同的。假設(shè)三維的圖像大小是N=（Nx，Ny，Nz）=（512，512，512），優(yōu)化算法一次需要估計(jì)的小波系數(shù)個(gè)數(shù)是512×512×512=134 217 128。如果估計(jì)全部的小波系數(shù)，這樣可以精確恢復(fù)形變域，但是需要耗費(fèi)大量的時(shí)間，這對(duì)于需要實(shí)時(shí)實(shí)現(xiàn)的自適應(yīng)放療系統(tǒng)是不能接受的。考慮到大多數(shù)的器官和組織產(chǎn)生的形變具有近似光滑性，我們用正交小波基函數(shù)表達(dá)形變時(shí)，只需要用一些低尺度的小波系數(shù)即可以較精確地恢復(fù)形變域，這對(duì)于配準(zhǔn)的實(shí)時(shí)實(shí)現(xiàn)具有重要的意義。采用正交小波變換7層分解，僅使用第6和第7尺度上的小波系數(shù)來(lái)估計(jì)形變域，實(shí)驗(yàn)表明，配準(zhǔn)結(jié)果具有較高的精度，而數(shù)據(jù)運(yùn)算量卻大大減少，對(duì)于N=（Nx，Ny，Nz）=（512，512，512）圖像大小，僅需要優(yōu)化算法總共估計(jì)8×8×8=512個(gè)參數(shù)，同時(shí)，子帶和低尺度上有較多的零系數(shù)，可以有效提高配準(zhǔn)速度。

圖1 二維圖示小波系數(shù)的估計(jì)次序，次序是：1-2-3-4-22-33-44Fig.1 Example for the progressive order in which registration parameters were estimated for 2D，from 1-2-3-4-22-33-44

為了簡(jiǎn)化在極小化能量函數(shù)時(shí)的復(fù)雜度，我們采用的正交小波基函數(shù)應(yīng)當(dāng)滿足三次正交性質(zhì)，即這些小波基函數(shù)自身相互正交，同時(shí)分別與自身的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)相互正交。滿足這一性質(zhì)可以使極小化能量函數(shù)的內(nèi)力部分得到運(yùn)算的簡(jiǎn)化［13］。采用接近三次正交的不完全正交三階樣條小波來(lái)進(jìn)行表示形變域。采用Marquardt-Levenberg（M-L）優(yōu)化算法對(duì)小波系數(shù)進(jìn)行估計(jì)［13］。詳細(xì)的配準(zhǔn)過(guò)程采用如下的步驟：

步驟1：輸入?yún)⒖紙D像R和目標(biāo)圖像T，去掉圖像中的床位投影。

步驟2：把形變域三維小波分解后，把小波系數(shù)的初始值置為0。

步驟3：使用M-L優(yōu)化算法遞歸的極小化式（8），尋找合適的小波系數(shù)，以小波系數(shù)的尺度和子帶作為循環(huán)參數(shù)進(jìn)行遞歸。具體方法采用步驟4和步驟5。

步驟4：根據(jù)小波分解塔式結(jié)構(gòu)，首先估計(jì)最粗尺度J上的小波系數(shù)，從代表概貌的子帶1開始，然后估計(jì)子帶2，一直到子帶8；接下來(lái)估計(jì)J-1尺度上的小波系數(shù)，從子帶2開始，然后估計(jì)子帶3，一直到子帶8。

步驟5：重復(fù)步驟4，一直到金字塔剩余的各個(gè)尺度，估計(jì)每一個(gè)尺度上子帶2到子帶8的小波系數(shù)，根據(jù)所需的精度控制在某個(gè)尺度上收斂。

步驟6：根據(jù)估計(jì)到的小波系數(shù)得到參考圖像與目標(biāo)圖像之間的形變域。

3 實(shí)驗(yàn)及結(jié)果分析

實(shí)驗(yàn)數(shù)據(jù)來(lái)自山東省腫瘤醫(yī)院Varian’s CBCT系統(tǒng)，選取20組計(jì)劃CT（planning CT，PCT）和治療時(shí)現(xiàn)場(chǎng)采集的對(duì)應(yīng)CBCT病例數(shù)據(jù)進(jìn)行配準(zhǔn)實(shí)驗(yàn)。CBCT與普通CT成像區(qū)別是一個(gè)角度直接形成三維的錐形束投影體數(shù)據(jù)，一周的照射即可快速形成三維的斷層成像，但是這種成像的特點(diǎn)會(huì)導(dǎo)致獲得的CBCT圖像總是含有偽影，同時(shí)含有較大的噪聲，如圖2所示。

圖2 CT圖像。（a）計(jì)劃CT；（b）CBCTFig.2 CT images.（a）planning CT；（b）cone beam CT image

為了定量驗(yàn)證算法的有效性，CBCT圖像可以認(rèn)為是由普通CT經(jīng)過(guò)灰度均衡變換（Gray_Trans）并加入乘性斑點(diǎn)噪聲（Noise）形成的，即：

式中，ImageCBCT表示變換后形成的模擬CBCT圖像，ImageCT表示原始的CT圖像，經(jīng)過(guò)式（9）的變換后形成CBCT，與原始的CT形成具有金標(biāo)準(zhǔn)性質(zhì)的配準(zhǔn)對(duì)。對(duì)配準(zhǔn)對(duì)中的CBCT圖像進(jìn)行定量的變形，變形方法采用ITK軟件包中的B樣條矢量變形方法?？紤]到實(shí)際病例中形變的多樣性，采用了3種由小到大的定量變形，如圖3所示。

同時(shí)，臨床治療時(shí)，CBCT與CBCT的配準(zhǔn)可以對(duì)腫瘤的變化或者病變器官進(jìn)行跟蹤，PCT與PCT的配準(zhǔn)可以比較治療前后的變化或者不同病例之間病變的差別。因此，對(duì)PCT和CBCT也進(jìn)行了同模態(tài)的形變配準(zhǔn)，同樣，對(duì)PCT或者CBCT分別進(jìn)行圖3中由小到大3種變形，然后再與自身進(jìn)行配準(zhǔn)。

通過(guò)計(jì)算已知的形變域與形變配準(zhǔn)算法獲得的形變域差值進(jìn)行定量分析配準(zhǔn)算法的性能，定義

式中，N是像素個(gè)數(shù)，Wi是由形變配準(zhǔn)算法計(jì)算出的第i個(gè)形變矢量的峰值，而Ki是第i個(gè)已知形變矢量的峰值。由于基于B樣條的自由形變配準(zhǔn)方法被學(xué)者認(rèn)為在配準(zhǔn)醫(yī)學(xué)圖像時(shí)具有優(yōu)勢(shì)［1］，因此，本研究的方法與基于B樣條的自由形變配準(zhǔn)方法進(jìn)行了比較。表1是分別使用本算法與基于B樣條的自由形變配準(zhǔn)算法配準(zhǔn)后計(jì)算出的形變域差值的結(jié)果。為了增加實(shí)驗(yàn)結(jié)果的可靠性，每個(gè)結(jié)果均取20組配準(zhǔn)數(shù)據(jù)對(duì)的均值。

圖3 3種已知的形變域。（a）～（c）從左到右的形變域分別由小到大，分別記為D1，D2，D3；（d）形變域的峰值Fig.3 Three level of deformation fields.From（a）～（c），the deformation is from small to large and we index them using D1，D2，D3；（d）meta-image magnitude for the deformation

表1 基于小波基函數(shù)和B樣條配準(zhǔn)方法得到的形變域差值比較Tab.1 Deformation difference comparison between the wavelet based and B-spline based methods

由表1可以看出，由本方法獲得的形變域差總小于基于B樣條的自由形變配準(zhǔn)方法得到的形變域差，這表明本方法配準(zhǔn)時(shí)可以獲得更高的精度，同時(shí)具有穩(wěn)定的魯棒性，可以更準(zhǔn)確地恢復(fù)形變域。同時(shí)，研究發(fā)現(xiàn)，對(duì)于較小的形變1，例如圖3中形變1，兩種方法均可準(zhǔn)確恢復(fù)，但在具有較大形變時(shí)，例如圖3中形變3，基于B樣條的自由形變配準(zhǔn)方法恢復(fù)性能較差，這種情況下，如果首先進(jìn)行全局的仿射變換，然后再進(jìn)行基于B樣條的配準(zhǔn)，會(huì)更好的恢復(fù)形變，達(dá)到和基于小波基函數(shù)配準(zhǔn)方法接近的性能。所提出的方法，即使在具有較大形變的情況下，也可以準(zhǔn)確快速地恢復(fù)形變域。

最后，采集了3組計(jì)劃CT和日常治療時(shí)相應(yīng)的CBCT圖像，包括形變較大的胸部，前列腺圖像和形變較小的頭頸部圖像，并把它們分別用基于本文的配準(zhǔn)方法和基于B樣條的自由形變配準(zhǔn)方法進(jìn)行了配準(zhǔn)。由于這種情況下沒(méi)有金標(biāo)準(zhǔn)，采用了配準(zhǔn)后的棋盤圖像來(lái)對(duì)比結(jié)果，并且臨床醫(yī)師和物理師對(duì)結(jié)果進(jìn)行了主觀評(píng)價(jià)。

圖4是計(jì)劃CT與CBCT配準(zhǔn)實(shí)驗(yàn)結(jié)果比較圖，由圖4可以看出（可從計(jì)劃CT與CBCT棋盤圖像銜接的位置觀測(cè)），基于B樣條的自由形變配準(zhǔn)方法在恢復(fù)形變性能上要比本文方法弱，特別是對(duì)于形變較大和較多的部位效果比較明顯。由于胸部和前列腺部位的形變相比頭頸部位要大，并且全局形變和局部形變特征均存在，本研究提出的方法更能表現(xiàn)出可以由粗到精的恢復(fù)全局和局部形變的良好性能。對(duì)于頭部圖像，形變本身較小，兩種方法均可以較準(zhǔn)確的恢復(fù)形變，但是觀測(cè)棋盤銜接處，本研究提出的方法在恢復(fù)頭頸部局部形變時(shí)同樣表現(xiàn)了更好的性能。實(shí)驗(yàn)結(jié)果經(jīng)過(guò)幾位臨床醫(yī)師和物理師的評(píng)價(jià)和分析，認(rèn)為本方法在恢復(fù)形變上，具有更好的性能，可以指導(dǎo)臨床的應(yīng)用。

在速度上，對(duì)于較小的形變時(shí)兩種方法均可準(zhǔn)確快速恢復(fù)形變，且平均速度近似相同，對(duì)于實(shí)驗(yàn)中采用的512×512×64的數(shù)據(jù)，平均時(shí)間均在4min左右。對(duì)于較大的形變，基于B樣條的配準(zhǔn)方法平均時(shí)間在8min左右，且在實(shí)驗(yàn)數(shù)據(jù)中發(fā)生誤配的概率在50％左右，而基于小波的配準(zhǔn)方法平均時(shí)間在6min左右，且在實(shí)驗(yàn)數(shù)據(jù)中均可以準(zhǔn)確配準(zhǔn)。因此，相對(duì)于基于B樣條的自由形變配準(zhǔn)方法，本方法具有速度和魯棒性上的優(yōu)勢(shì)。

4 臨床應(yīng)用

計(jì)劃CT與CBCT的配準(zhǔn)最直接的應(yīng)用就是在IGRT系統(tǒng)放射治療時(shí)，配準(zhǔn)結(jié)果的形變域參數(shù)指導(dǎo)治療靶區(qū)的射線進(jìn)行自適應(yīng)調(diào)整，使靶區(qū)的劑量攝取合理，而對(duì)周圍正常組織器官達(dá)到合理限度的傷害，這需要把算法移植入機(jī)器，然后控制機(jī)器運(yùn)動(dòng)，臨床實(shí)際應(yīng)用仍在討論當(dāng)中。另一個(gè)方面的重要應(yīng)用是計(jì)劃CT與CBCT的準(zhǔn)確配準(zhǔn)可以自動(dòng)準(zhǔn)確提取CBCT的組織器官，代替耗費(fèi)大量物理師時(shí)間的手動(dòng)分割。由于放射治療時(shí)，計(jì)劃CT中的器官在治療之前由物理師進(jìn)行人工分割，為了在治療時(shí)或者是在分次治療中計(jì)算劑量在靶區(qū)和正常器官上的分布，以便于劑量更加合理的分布，可以利用形變配準(zhǔn)算法對(duì)治療時(shí)或者是分次治療中的CBCT圖像進(jìn)行器官和靶區(qū)的自動(dòng)分割。

圖4 胸部病例，前列腺病例，頭頸部病例計(jì)劃CT與CBCT配準(zhǔn)后結(jié)果棋盤圖像（在每組棋盤圖像中，上行是本文提出的方法的配準(zhǔn)結(jié)果，下行是用基于B樣條的自由形變配準(zhǔn)方法得到的配準(zhǔn)結(jié)果）。（a）胸部配準(zhǔn)結(jié)果；（b）前列腺部位配準(zhǔn)結(jié)果；（c）頭頸部位病例配準(zhǔn)結(jié)果Fig.4 Checkboards after planning CT and CBCT deformable registration（in each group，the upper row is with wavelet based method，the bottom row is with FFD based method）.（a）registration result for chest images；（b）registration result for prostate images；（c）registration result for head and neck（H＆N）

利用配準(zhǔn)算法進(jìn)行器官自動(dòng)分割的具體方法是：治療前，物理師在制定放療計(jì)劃時(shí)，對(duì)CT圖像上器官和靶區(qū)的輪廓進(jìn)行勾勒。由于計(jì)劃 CT與CBCT可以準(zhǔn)確配準(zhǔn)，這些在計(jì)劃CT上的器官和靶區(qū)根據(jù)計(jì)劃CT與CBCT之間的形變域變化，把計(jì)劃CT的器官輪廓重新映射到CBCT數(shù)據(jù)上。對(duì)臨床20個(gè)病例的CBCT圖像中肝臟數(shù)據(jù)進(jìn)行了自動(dòng)分割，下圖是隨機(jī)抽取了其中3個(gè)病例圖像中的8片數(shù)據(jù)，進(jìn)行肝臟的分割顯示。

如圖5所示，經(jīng)過(guò)山東省腫瘤醫(yī)院多位醫(yī)師評(píng)價(jià)，由配準(zhǔn)算法自動(dòng)分割得到的CBCT圖像肝臟部分與人工分割得到的CBCT肝臟部分基本相同，經(jīng)過(guò)計(jì)算，平均重疊區(qū)域在97％以上，可以用于制定計(jì)劃系統(tǒng)。這主要得益于配準(zhǔn)算法可以精確恢復(fù)形變，把在計(jì)劃CT上的人工分割輪廓經(jīng)過(guò)形變域映射到治療中的CBCT圖像上，從而由原先的計(jì)劃CT肝臟邊界輪廓轉(zhuǎn)化為CBCT上的肝臟邊界輪廓。其它器官以及靶區(qū)的分割類似，這些器官和靶區(qū)的自動(dòng)分割對(duì)于在治療過(guò)程中快速合理計(jì)算放射劑量的分布具有重要的作用。

圖5 由配準(zhǔn)算法自動(dòng)分割CBCT病例圖像中肝臟結(jié)果（每個(gè)病例隨機(jī)抽取了其中8片的分割結(jié)果。（a）～（c）分別為其中3個(gè)病例Fig.5 The liver automatic segmentation result using deformable registration algorithm for CBCT images（8 slicers are selected randomly from each case）.（a）～（c）are the results of three patients，respectively

5 結(jié)論和討論

本研究提出了改進(jìn)的基于CBCT的圖像引導(dǎo)放射系統(tǒng)中計(jì)劃CT與日常治療中的CBCT圖像的配準(zhǔn)算法，采用基于正交小波基函數(shù)來(lái)表述形變域，這些不同尺度的小波基函數(shù)分別表示了不同的形變特征，粗尺度的小波系數(shù)代表了全局的形變，精細(xì)尺度的小波系數(shù)代表了局部的形變，不同子帶的小波系數(shù)也表征特有的形變細(xì)節(jié)。同時(shí)利用Navier偏微分方程來(lái)設(shè)計(jì)形變能量函數(shù)，方程的自由膨脹部分作為外力，內(nèi)力部分利用互信息作為距離函數(shù)來(lái)平衡方程的外力。表示形變域的小波部分通過(guò)優(yōu)化算法由粗尺度到精細(xì)尺度逐步估計(jì)得到。實(shí)驗(yàn)表明，相對(duì)于傳統(tǒng)的基于B樣條的自由形變配準(zhǔn)方法，本文提出的配準(zhǔn)方法具有精度和速度上的優(yōu)勢(shì)。同時(shí)，本方法可用于其它用途和其它種類的醫(yī)學(xué)圖像配準(zhǔn)。

［1］Holden M.A review of geometric transformations for nonrigid body registration［J］.IEEE Transactions on Medical Imaging，2008，27：111 -128.

［2］Sorcini B，Tilikidis A.Clinical application of image guided radiotherapy，IGRT（on the varian OBI platform）［J］.Cancer Radiothérapie，2006，10：252 -257.

［3］Jean P，Ali H，Chen J，et al.Low-dose megavoltage conebeam CT for radiation therapy［J］.International Journal of Radiation Oncology，Biology，Physics，2005，61：552 -560.

［4］Cheng B，Yang Yong，Li Fang，et al.Performance characteristics and quality assurance aspects of kilovoltage conebeam CT on medical linear accelerator［J］.Medical dosimetry：official journal of the American Association of Medical Dosimetrists，2007，32：80-85.

［5］Paquin D，Levy D，Xing Lei.Multiscale registration of planning CT and daily cone beam CT images for adaptive radiation therapy［J］.Medical Physics，2009，36：4-11.

［6］Brock KK，Dawson LA，Sharpe MB，et al.Feasibility of a novel deformable image registration technique to facilitate classification，targeting，and monitoring of tumor and normal tissue［J］.International Journal of Radiation Oncology，Biology，Physics，2006，64：1245-1254.

［7］Nithiananthan S，Brock KK，Daly MJ，etal.Demons deformable registration for CBCT-guided procedures in the head and neck：convergence and accuracy［J］.Medical Physics，2009，36：4755-4764.

［8］Noe K，Senneville BD，Vindelev U，et al.Acceleration and validation of optical flow based deformable registration for image guided radiotherapy［J］.Acta Oncologica，2008，47：1286 -1293.

［9］Bajcsy R，Kovacic S.Multiresolution elastic matching［J］.Computer Vision，Graphics，and Image Processing，1989，46：1-12.

［10］Pluim JPW，Maintz JBA，Viergever MA.Mutual information matching in multiresolution contexts［J］.Image and Vision Computing，2001，19：45-52.

［11］Mallat SG.A wavelet tour of signal processing［M］.San Diego：Academic Press，1998.

［12］彭玉華.小波變換與工程應(yīng)用［M］北京：科學(xué)出版社，1999.

［13］Gefen S，Tretiak O，Nissanov J.Elastic 3-D alignment of rat brain histological images［J］.IEEE Transactions on Medical Imaging，2003，22：1480-1489.

［14］Pluim JPW，Maintz JBA，Viergever MA.Mutual-informationbased registration of medical images：a survey［J］.IEEE Transactions on Medical Imaging，2003，22：986-1004.