范忠銀 王建軍 朱青松 謝耀欽 邢 磊

(中國(guó)科學(xué)院深圳先進(jìn)技術(shù)研究院 深圳 518055)

配準(zhǔn)算法對(duì) PCA 單幅投影肺部重建的影響

范忠銀 王建軍 朱青松 謝耀欽 邢 磊

(中國(guó)科學(xué)院深圳先進(jìn)技術(shù)研究院 深圳 518055)

肺癌放射治療中，肺部腫瘤位置實(shí)成像對(duì)于臨床意義重大。在一種利用單 X 射線投影進(jìn)行成像的實(shí)時(shí)肺部 3D 成像算法中，圖像配準(zhǔn)過(guò)程引入的不準(zhǔn)確對(duì)于 PCA 模型構(gòu)建以及重建過(guò)程有重大影響。文章分析了光流法、Demons 算法、水平集算法三種配準(zhǔn)算法對(duì)重建效果的影響，并通過(guò)定性以及定量實(shí)驗(yàn)分析驗(yàn)證。結(jié)果表明，光流法配準(zhǔn)在配準(zhǔn)結(jié)果以及模型構(gòu)建方面有較好的效果。

實(shí)時(shí)肺部 3D 成像；主成分分析；彈性配準(zhǔn)

1 引 言

在現(xiàn)代肺癌的放射治療中，治療過(guò)程中對(duì)實(shí)時(shí)肺部腫瘤位置進(jìn)行精確定位十分重要[1]。肺部腫瘤定位方法可以大致分為：直接腫瘤定位、基于呼吸參照定位以及基于植入標(biāo)記物定位。現(xiàn)今的直接定位技術(shù)存在以下缺點(diǎn)：(1)方法需從前后方向獲得圖像且腫瘤邊界清晰時(shí)才能有效運(yùn)行；(2)需要經(jīng)過(guò)預(yù)處理的熒光透視圖像；(3)無(wú)法精確得到垂直于 X 射線成像平面的第三維腫瘤位置信息?；诤粑鼌⒄瘴锏亩ㄎ环椒ú恍枰~外的放射劑量來(lái)進(jìn)行定位，但是腫瘤位置和參照物內(nèi)外部關(guān)系的微小變化對(duì)定位有很大的影響[2,3]?；谥踩霕?biāo)記物的定位方法定位精度較高[4]，但這種方法需要侵入人體，且易引起氣胸[5]。

要克服這些不足，需構(gòu)建一個(gè)用于直接定位的精確肺部運(yùn)動(dòng)模型。Zhang 等[6]提出一種基于 主成分分析(Principal Component Analysis，PCA)的肺部運(yùn)動(dòng)模型，只需要較少的特征向量和 PCA 系數(shù)就可以精確表示肺部運(yùn)動(dòng)。類似的，Li 等[7,8]將序列相位與參考相位進(jìn)行彈性配準(zhǔn)得到偏移矢量場(chǎng)(Deformation Vector Field，DVF)，計(jì)算 DVF 本征模式，建立一個(gè)類似于Zhang 的 PCA 模型，重建一個(gè) 4D-CT。然后利用 PCA 模型對(duì) 4D-CT 中的參考相位圖像進(jìn)行變形，將單個(gè)的 DRRs(Digitally Reconstructed Radiographs)與 CBCT(Cone Beam Computed Tomography)投影匹配來(lái)估計(jì)原始 4D CBCT 投影數(shù)據(jù)集每個(gè)時(shí)間步長(zhǎng)的本征模式系數(shù)(即主要的系數(shù))。PCA 肺部運(yùn)動(dòng)模型隱式配準(zhǔn)且效率高，使其能精確、高效率地得到一個(gè)動(dòng)態(tài)肺部運(yùn)動(dòng)模型，且只需要少量信息。

要將運(yùn)動(dòng)模型和先驗(yàn)圖像綜合起來(lái)，最主要的挑戰(zhàn)是得到一個(gè)真實(shí)的模型，且模型具有合適數(shù)量的參數(shù)，并可以簡(jiǎn)單地在任何成像物體上歸納得出。PCA 肺部運(yùn)動(dòng)模型基于先驗(yàn)圖像(4DCT)之間的彈性配準(zhǔn)得到 DVF，這給 PCA 分析引入誤差的源頭：圖像配準(zhǔn)過(guò)程引入的不準(zhǔn)確及近似問(wèn)題。

在本文中，我們將對(duì)比幾種主流配準(zhǔn)方法在PCA 構(gòu)建模型以及重建效果的差別，并甄選出最合適的配準(zhǔn)算法。

2 基于 PCA 模型的單幅投影肺部重建方法

2.1 PCA 模型

來(lái)自于仿真治療過(guò)程的 4D-CT 數(shù)據(jù)作為先驗(yàn)知識(shí)，以此進(jìn)行三維圖像重建。在參考位相和其他位相之間，進(jìn)行彈性圖像配準(zhǔn)，產(chǎn)生一系列的 DVF。這些 DVF 能夠用 PCA 方法得到的一系列特征向量和系數(shù)來(lái)表示。通過(guò)改變 PCA 系數(shù)，能夠得到新的 DVF。把新的 DVF 應(yīng)用于參考圖像就能得到新的立體圖像。然后優(yōu)化 PCA系數(shù)(運(yùn)用構(gòu)建出的模型)使計(jì)算出來(lái)的投影數(shù)據(jù)與測(cè)量得到的立體圖像投影數(shù)據(jù)一致。將反轉(zhuǎn)的DVF 應(yīng)用于參考圖像就可以得到待重建圖像。

在 PCA 肺部運(yùn)動(dòng)模型中，DVF(相對(duì)于參考圖像是一個(gè)空間和時(shí)間函數(shù))可以通過(guò)樣本均值向量和一系列具有最大特征值的特征向量的線性組合近似估計(jì)。如下：

其中，x(t)是一個(gè)參數(shù)化的 DVF，即一個(gè)時(shí)間和空間的函數(shù)； 表示不同時(shí)間 DVF 的均值；uk表示由 PCA 得到的特征向量，是空間函數(shù)；wx(t)表示 PCA 系數(shù)，是時(shí)間的函數(shù)。因?yàn)槊恳粋€(gè)空間像素元在每一個(gè)特征向量都有一個(gè)自由度。所以，每一個(gè)特征向量定義了空間中體素運(yùn)動(dòng)的一個(gè)方向，兩個(gè)特征向量就是一個(gè)平面，三個(gè)特征向量構(gòu)成腫瘤運(yùn)動(dòng)的三維空間。這種特性使得 PCA 運(yùn)動(dòng)模型能夠準(zhǔn)確把握腫瘤運(yùn)動(dòng)的軌跡。在本文的方法中，保留具有最大特征值的 3個(gè) PCA 系數(shù)。

在得到一個(gè)參數(shù)化的 PCA 肺部運(yùn)動(dòng)模型之后，我們希望能夠得到一系列最優(yōu)化 PCA 系數(shù)，使得用新的 DVF 重建的立體圖像投影數(shù)據(jù)與測(cè)量得到的 X 射線投影數(shù)據(jù)能夠匹配。f0作為參考圖像，f是重建圖像，y是測(cè)量得到的投影圖像，P是用來(lái)計(jì)算f的投影圖像的投影矩陣。代價(jià)函數(shù)如下：

其中，

其中，U是一個(gè)矩陣，它的列是 PCA 的特征向量；w是一個(gè)包含了待優(yōu)化的 PCA 系數(shù)的向量；x是參數(shù)化的 DVF；為p階泛數(shù)。

對(duì)代價(jià)函數(shù)進(jìn)行優(yōu)化：

其中

每次迭代之后，由公式(4)得到更新的 PCA系數(shù)，從而更新 DVF，由三線性插值得到重建圖像，對(duì)應(yīng)地，必須與插值過(guò)程保持一致，從而得到正確的梯度值。

2.2 主流配準(zhǔn)方法

在 PCA 模型中，參考相位和其他相位之間配準(zhǔn)生成 DVF，這些 DVF 決定了建模的準(zhǔn)確度。配準(zhǔn)算法對(duì)于整個(gè)重建有著重大的影響。為研究圖像配準(zhǔn)方法對(duì) PCA 建模以及重建效果的影響，選取三種經(jīng)典的配準(zhǔn)方法：光流法(Optical Flow Method，OFM)、Demons 法、水平集法(Levelset)進(jìn)行比較研究。

其中，光流場(chǎng)方法起源于計(jì)算機(jī)視覺(jué)領(lǐng)域，主要用于運(yùn)動(dòng)物體檢測(cè)，由于運(yùn)動(dòng)物體的光流場(chǎng)與圖像配準(zhǔn)的位移場(chǎng)具有一定的相似性，因此光流場(chǎng)的基本原理也可以用于圖像配準(zhǔn)。主流的光流場(chǎng)方法是基于 Horn 等[9]提出的變分法，該方法基于像素灰度值的局部梯度模型和整體平滑性假設(shè)。相對(duì)于 Horn 的全局化模型，Lucas 等[10,11]提出了局部化窗口模型，其思想是將每一像素速度矢量的求解放到其鄰域中，由鄰域內(nèi)各像素的信息共同約束該點(diǎn)速度矢量。窗口化操作中，每一像素只受其周圍鄰域像素的影響，可以通過(guò)限定小的窗口來(lái)減少對(duì)圖像造成的模糊。本課題研究中采用窗口化的 Lucas 模型進(jìn)行圖像配準(zhǔn)。

Demons 算法配準(zhǔn)是目前最為流行的配準(zhǔn)方法。Demons 模型由 Thirion[12]于 1998 年提出并應(yīng)用于圖像配準(zhǔn)領(lǐng)域。將源圖像視為一個(gè)可變形網(wǎng)絡(luò)，其形變自由度可以人工設(shè)定。圖像像素對(duì)應(yīng)于網(wǎng)格的格點(diǎn)，并且都被賦予一定的極性：內(nèi)點(diǎn)或外點(diǎn)。Demons 根據(jù)格點(diǎn)的極性來(lái)控制方向，內(nèi)點(diǎn)只能進(jìn)入某個(gè)區(qū)域，外點(diǎn)只能離開某個(gè)區(qū)域。如果將源圖像放入布滿 Demons 的目標(biāo)圖像空間內(nèi)，則網(wǎng)格就會(huì)在 Demons 的控制下發(fā)生形變，最終所有內(nèi)點(diǎn)都將進(jìn)入某一區(qū)域內(nèi)，而外都將被退出到某一區(qū)域外。通過(guò)合理設(shè)置所有格點(diǎn)的極性，可以使得源圖像自動(dòng)發(fā)生形變，達(dá)到與目標(biāo)圖像相似程度最大。

水平集方法是 Sethian[13]在研究曲線以曲率相關(guān)的速度演化時(shí)提出來(lái)的，用于描述曲線或曲面的演化過(guò)程，是處理封閉界面在演化過(guò)程中幾何拓?fù)渥兓挠辛ぞ摺Ｆ渥畲蟮膬?yōu)勢(shì)在于它的穩(wěn)定性以及拓?fù)錈o(wú)關(guān)性。但 水平集方法有一個(gè)缺點(diǎn)，即計(jì)算量太大。

2.3 算法流程

本文選取了 Demons 配準(zhǔn)算法、水平集配準(zhǔn)算法和 OFM 配準(zhǔn)算法作為對(duì)比方法，來(lái)考察它們?cè)?PCA 構(gòu)建模型以及重建效果上的差別。

算法過(guò)程如下：

(1)數(shù)據(jù)準(zhǔn)備：10 個(gè)相位的肺部 CT 三維數(shù)據(jù)；

(2)選取第一個(gè)相位作為參考圖像，其他相位作為浮動(dòng)圖像，進(jìn)行彈性圖像配準(zhǔn)；

(3)PCA 提取主成分；

(4)利用主成分進(jìn)行迭代重建。

算法流程如圖 1 所示。

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

3.1 實(shí)驗(yàn)數(shù)據(jù)

本文研究的圖像數(shù)據(jù)為 4D-CT 數(shù)據(jù)。CT 掃描結(jié)束后，利用 4D 軟件對(duì)得到的大量不同位置、不同呼吸時(shí)相的數(shù)據(jù)進(jìn)行處理。以 10% 為時(shí)相間隔，把一個(gè)呼吸周期分為 10 個(gè)呼吸時(shí)相。分別命名為 0%(吸氣末)、10%(吸氣中)、20%、…、50%(呼氣末)、…、90%(吸氣末)，依次定義為 0 相位至 9 相位。

針對(duì)肺部的形變情況，我們截取了 CT 影像序列中肺腔的部分，這樣可以減少不必要的數(shù)據(jù)量，加快算法的速度。本文使用的是連續(xù) 0～127層 CT 影像信息，數(shù)據(jù)大小為 256×256×128，體素分辨率 1 mm×1 mm×1.25 mm。投影大小為 300×400，探測(cè)器像素 1 mm×1 mm。

圖 1 算法流程圖Fig. 1 The flow chart of the proposed algorithm

圖 2 是呼吸周期 10 個(gè)呼吸時(shí)相中間層(第 64層)的圖像。在圖像中顯示為黑色區(qū)域的為肺組織，可看出肺輪廓的變化情況。而毫米級(jí)的組織運(yùn)動(dòng)變化對(duì)于放療計(jì)劃都將有重大的影響，所以需要采集實(shí)時(shí)的三維體數(shù)據(jù)以根據(jù)人體組織運(yùn)動(dòng)變化制定治療計(jì)劃。

3.2 彈性圖像配準(zhǔn)

選擇呼吸周期 0 相位的圖像作為參考圖像，其他 9 個(gè)相位作為待配準(zhǔn)圖像。

圖 2 呼吸周期 10 個(gè)呼吸時(shí)相中間層(第 64 層)的圖像Fig. 2 The 64th slice of each phase in a respiratory circle

圖 3 是待配準(zhǔn)圖像為呼吸周期 50% 時(shí)相，用三種配準(zhǔn)方法得到的結(jié)果。與標(biāo)準(zhǔn)待配準(zhǔn)圖像對(duì)比，我們發(fā)現(xiàn) OFM 配準(zhǔn)的結(jié)果更加接近待配準(zhǔn)圖像，在圖像邊緣部分尤為明顯。其中，從上到下依次為從冠狀面、矢狀面和橫斷截取中間層的圖像。

表 1 為對(duì)配準(zhǔn)進(jìn)行定量分析的結(jié)果。其中，配準(zhǔn)誤差采用公式(6)計(jì)算。不同配準(zhǔn)方法得到的位置偏移不同，Demons 方法對(duì)配準(zhǔn)偏移估計(jì)過(guò)大，而 水平集方法對(duì)偏移估計(jì)較小，OFM 得到了一個(gè)較為適中的偏移，且 OFM 配準(zhǔn)誤差更小。

從圖 3 定性的說(shuō)明以及表 1 定量的計(jì)算，都可以看出：OFM 配準(zhǔn)的結(jié)果較為準(zhǔn)確。

圖 3 待配準(zhǔn)圖像為呼吸周期 50% 時(shí)相，用三種配準(zhǔn)方法得到的結(jié)果Fig. 3 The registration deformed image at 50% phase with three methods

表 1 三種配準(zhǔn)方法的定量對(duì)比Table 1 The quantitative comparison of three methods

3.3 PCA 結(jié)果

將配準(zhǔn)步驟得到的x、y、z方向上的偏移矢量進(jìn)行 PCA。

表 2 是三種算法配準(zhǔn)偏移矢量進(jìn)行 PCA 之后x、y、z方向上各個(gè)特征向量的貢獻(xiàn)率，最后一行為前三個(gè)特征向量的累計(jì)貢獻(xiàn)率?？梢钥闯觯珼emons 算法和 OFM 算法的結(jié)果 PCA 能量較為集中，前三個(gè)特征向量累計(jì)貢獻(xiàn)率也較高，從而對(duì)后續(xù)的重建更加有利。

3.4 重建結(jié)果

圖 4 為三種方法結(jié)果對(duì)比情況。其中，從上到下依次為冠狀面、矢狀面和橫斷截取中間層的圖像。從圖像的輪廓以及灰度級(jí)初略查看，三種方法的重建結(jié)果差別不大。

采用定量的方法加以考察，其中重建結(jié)果的配準(zhǔn)誤差分別為：Demons 為 0.0145727，水平集為 0.0130284，OFM 為 0.0129043?？梢缘玫?，OFM 方法得到的重建圖像與標(biāo)準(zhǔn)圖像誤差最小。

表 2 PCA 累積貢獻(xiàn)率Table 2 The accumulative score of the principal component analysis

4 討 論

本文介紹的單個(gè) X 射線圖像重建肺部三維影像算法在臨床應(yīng)用具有廣闊的潛力。本項(xiàng)研究中腫瘤定位誤差來(lái)自兩個(gè)不同的過(guò)程：(1)PCA 肺部運(yùn)動(dòng)模型；(2)投影和參考圖像之間的 2D/3D配準(zhǔn)。一個(gè)可靠的彈性配準(zhǔn)算法，特別是能夠精確地模擬較大位移運(yùn)動(dòng)的配準(zhǔn)算法，從而提高PCA 肺部運(yùn)動(dòng)模型的精確度。因此，我們采用Demons、水平集和 OFM 配準(zhǔn)算法來(lái)對(duì) PCA 肺部運(yùn)動(dòng)模型的精確度和重建效果進(jìn)行分析。綜合以上分析，在三種方法中，光流法配準(zhǔn)結(jié)果對(duì)于PCA 模型構(gòu)建和重建效果較好，而 水平集算法在配準(zhǔn)結(jié)果以及重建方面性能都較差。本文選取的配準(zhǔn)方法有限，下一步將繼續(xù)擴(kuò)大配準(zhǔn)算法的范圍來(lái)對(duì)模型構(gòu)建和重建進(jìn)行對(duì)比。

在治療過(guò)程中，病人可能在分次間發(fā)生解剖變化，如果肺部運(yùn)動(dòng)模型是基于病人模擬仿真獲得的 4D-CT 而建立，那么結(jié)果將會(huì)很不理想。在這種情況下，病人擺位過(guò)程采集的 4D CBCT建立的 PCA 肺部運(yùn)動(dòng)模型能克服這個(gè)問(wèn)題。另一個(gè)潛在的問(wèn)題是當(dāng)腫瘤的位置靠近心臟時(shí)，腫瘤的運(yùn)動(dòng)不僅會(huì)受到呼吸運(yùn)動(dòng)影響，同時(shí)也會(huì)受到無(wú)法用 PCA 運(yùn)動(dòng)模型表達(dá)的心臟運(yùn)動(dòng)影響。這是 4D-CT 的一個(gè)局限，4D-CT 與呼吸同步而非心臟運(yùn)動(dòng)。在研究中，我們假設(shè)測(cè)量的錐束投影和待重建圖像的仿真投影之間的灰度值有線性關(guān)系，但是未考慮到不同圖像模態(tài)直接的差異，用于先驗(yàn)的 4D-CT 和術(shù)中重建的 CBCT 的 X 射線能量值和光譜值均不一樣。綜合這些意見，下一步將采用病人擺位過(guò)程采集的 4D CBCT 作為先驗(yàn)圖像來(lái)改進(jìn)算法。

[1] Keall PJ, Mageras GS, Balter JM, et al. The management of respiratory motion in radiation oncology report of AAPM Task Group 76 [J]. Medical Physics, 2006, 33(10): 3874-3900.

[2] Hoisak JD, Sixel KE, Tirona R, et al. Correlation of lung tumor motion with external surrogate indicators of respiration [J]. International Journal of Radiation Oncology, Biology, Physics, 2004, 60(4): 1298-1306.

[3] Tsunashima Y, Sakae T, Shioyama Y, et al. Correlation between the respiratory waveform measured using a respiratory sensor and 3D tumor motion in gated radiotherapy [J]. International Journal of Radiation Oncology, Biology, Physics, 2004, 60(3): 951-958.

[4] Harada T, Shirato H, Ogura S, et al. Real-time tumor-tracking radiation therapy for lung carcinoma by the aid of insertion of a gold marker using bronchofiberscopy [J]. Cancer, 2002, 95(8): 1720-1727.

[5] Laurent F, Latrabe V, Vergier B, et al. CT-guided transthoracic needle biopsy of pulmonary nodules smaller than 20 mm: results with an automated 20-gauge coaxial cutting needle [J]. Clinical Radiology, 2000, 55(4): 281-287.

[6] Zhang QH, Pevsner A, Hertanto A, et al. A patientspecific respiratory model of anatomical motion for radiation treatment planning [J]. Medical Physics, 2007, 34(12): 4772-4781.

[7] Li R, Lewis JH, Jia X, et al. On a PCA-based lung motion model [J]. Physics in Medicine and Biology, 2010, 56(18): 6009-6030.

[8] Li R, Jia X, Lewis JH, et al. Real-Time volumetric image reconstruction and 3D tumor localization based on a single X-ray projection image for lung cancer radiotherapy [J]. Medical Physics, 2010, 37(6): 2822-2826.

[9] Horn BKP, Schunck BG. Determining optical flow [J]. Artificial Intelligence, 1981, 17: 185-203.

[10] Lucas B, Kanade T. An iterative image registration technique with an application to stereo vision [C] // Proceedings of the International Joint Conference on Artificial Intelligence, 1981: 674-679.

[11] Lucas BD. Generalized image matching by the method of differences [D]. Pittsburgh: Carnegie Mellon University, 1984.

[12] Thirion JP. Image matching as a diffusion process: an analogy with Maxwell’s Demons [J]. Medical Image Analysis, 1998, 2(3): 243-260.

[13] Osher S, Sethian JA. Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulation [J]. Journal of Computational Physics, 1988, 79(1): 12-49.

Deformable Image Registration Effect on PCA-based Reconstruction with Single X-ray Imaging

FAN Zhongyin WANG Jianjun ZHU Qingsong XIE Yaoqin XING Lei

(Shenzhen Institutes of Advanced Technology,Chinese Academy of Sciences,Shenzhen518055,China)

In modern lung cancer radiotherapy, it is important to have a precise knowledge of the real-time lung tumor position during the treatment delivery. For a real-time 3D lung imaging algorithm from a single X-ray projection image, the inaccuracies contributed by the image registration process affects much on the PCA modelling and construction process. We utilize 3 deformable image registration algorithms: Optical Flow method, Demons method and Levelset method to evaluate the effect. By making quantitative analysis and qualitative analysis, we get the conclusion: Optical Flow method works much better in registration and PCA modelling.

real-time 3D tumor localization; principal component analysis; deformable image registration

R 812

A

2014-04-22

： 2015-03-15

廣東省引進(jìn)科研創(chuàng)新團(tuán)隊(duì)項(xiàng)目(2011S013)

范忠銀，碩士研究生，研究方向?yàn)獒t(yī)學(xué)圖像處理、醫(yī)學(xué)影像系統(tǒng)；王建軍，碩士研究生，研究方向?yàn)橛?jì)算機(jī)控制與三維可視化；朱青松，碩士，研究方向?yàn)橛?jì)算機(jī)視覺(jué)、醫(yī)療機(jī)器人視覺(jué)；謝耀欽(通訊作者)，博士，研究員，研究方向?yàn)閳D像引導(dǎo)放射治療、醫(yī)學(xué)影像處理和分析，E-mail：yq.xie@siat.ac.cn；邢磊，博士，教授，研究方向?yàn)槟[瘤放射物理。