PET/CT成像呼吸運動B樣條校正

潘李鵬，賀建峰，封 碩，崔 銳，馬 磊，相 艷，易三莉，張 俊

昆明理工大學 信息工程與自動化學院，昆明 650500

1 引言

在PET/CT掃描成像過程中，病人的各種運動往往會造成偽影的產生，尤其是呼吸心跳等運動產的偽影，會嚴重影響PET/CT圖像的質量。它將會干擾醫生的診斷、制定治療計劃和治療效果評估。針對呼吸運動校正，目前常用方法之一有呼吸門控方法[1]。它把呼吸運動劃分成不同的時間相位（時相），運用特定外部設備跟蹤呼吸運動的位置狀態，再將不同呼吸周期中具有相同時相的數據疊加后再重建。這樣各時相的呼吸偽影便大大減輕。一般門控方法需要額外的硬件（如紅外監測裝置、電磁裝置、彈性傳感帶等）[2-4]來監測呼吸運動的位置狀態。

呼吸門控除了用上面提到的硬件跟蹤方法外，也采用基于掃描數據的數據驅動（data-driven）門控方法[5-9]。這種方法不需要硬件跟蹤設備，而是基于已采集的數據通過算法估算出具有相同位置狀態的時相。但是目前現有的數據驅動方法基本上是在圖像空間內進行操作來進行門控校正的，即在已重建的圖像序列上，運用不同的算法分析處理圖像序列估算出運動時相。由于是在圖像序列上估算運動參數，這類方法需要重建大量圖像序列，運算量非常大，而且其精度仍需進一步改善[5]。例如，數據驅動方法中常用的質點中心方法（Centre of Mass），根據每幅圖像序列不同狀態的質心點（例如圖像灰度值質心點）差異，估算出呼吸時相進行門控[7，9]。這種方法在呼吸運動速度過快，超出掃描范圍或掃描數據采集不足的情況下會影響估算結果的精確性。

利用硬件或數據驅動的門控方法雖然能夠降低運動成像偽影，但會導致圖像的信噪比大大降低。這是因為它只提取了掃描總數據的10%～20%進行門控處理，而大部分有效數據被白白浪費了[10]。為了提高信噪比，只有利用更多的掃描數據進行門控。這可以將掃描時間加大到正常時間的5～10倍；或者成比例地增加放射性示蹤劑的含量，但勢必會使患者接受更多的輻射劑量，是不可取的。

為此，研究人員提出用圖像配準的方法利用所有掃描采集數據進行圖像校正。這種方法是先估算出已門控時相配準參數，然后再根據配準參數將所有已門控時相數據全部配準到參考相位以提高圖像校正信噪比[11-14]。例如光流法（optical flow algorithm），在已重建所有圖像序列的基礎上，運用算法估算出每兩幅圖像之間每個像素的運動狀態來確定其配準運動參數，然后進行圖像配準[15-16]。它需要重建圖像并對所有序列圖像里的像素進行處理，其算法的復雜性和運算量不應忽視[5]。

為此，本文提出一種基于運動特征的呼吸運動校正方法，通過時間幀對應提取與PET圖像相對應的CT圖像運動特征，來配準校正呼吸運動圖像。

2 方法

PET圖像因為其固有的缺陷和不足使得圖像質量較差噪聲也較高，這就造成了PET圖像的結合邊緣比較模糊。這樣的圖像中所提取的運動信息的可靠度較低，難以滿足臨床運用對精度的需求。在PET/CT的一次掃描中就可以同時獲取PET和CT的圖像，因此本文提出通過對CT圖像進行運動特征信息提取，然后使用這些信息對PET圖像進行呼吸運動校正。本文所使用的運動特征信息是一系列的控制點，通過這些控制點對PET圖像進行B樣條變換就可以達到校正目的。

2.1 基于B樣條的FFD配準方法

在本文中進行圖像配準的目的是通過對不同時間幀狀態下不同呼吸相的肺部CT圖像進行對比和校正，得出用于校正PET圖像所需的運動特征信息。而這種配準的實質是尋找一最優變換T：(x，y，z)→(x′，y′，z′)，此式映射在一個動態圖像序列I(x，y，z，t)在時刻t上任意一 點對應到時刻t0圖像I(x′，y′，z′，t0) 上點的關系。分析人體肺部呼吸運動特性可知，肺部呼吸運動主要表現形式為非剛性形變運動。

基于B樣條[17-18]的配準已被用于心臟圖像運動分析和乳房MRI圖像配準等應用上，證明其可以被用于人體非剛性組織的運動配準。這里選擇基于B樣條的FFD（Free form Deformation）方式來對肺部影像進行處理。本文首次提出使用基于B樣條的配準來提取CT圖像的運動信息，用于PET圖像的校正。FFD是一種用于對立體模型進行任意自由變形的一種技術，其可以對任意平面、角度、參數的表面進行變形。這種變換可以是全局性的也可以是局部性的。FFD的基本思想是使用由一系列的控制所組成的網格來控制定義一個物體的形變。

為了定義一個基于FFD的樣條，定義圖像空間為Ω={(x，y)|0≤x＜X，o≤y＜Y}。Φ表示一個在均勻空間σ內nx×ny網格控制點Φi，j。然后FFD可以被寫成熟悉的一維三次B樣條張量積：

其中i=(x/nx)-1，j=(y/ny)-1，u=x/nx-(x/n)，v=y/ny-(y/n)，并且Bl表示第l個B樣條基函數：

相對于 thin-plate樣條[19]或者 elastic-body 樣條[20]，B樣條擁有更好的局部控制特性，能夠有效地計算大量的控制點。

控制點集Φ作為B樣條FFD和非剛性變換角度的參數。非剛性變換角度可以通過控制點Φ構成網格的分辨率來表示。當控制點網格間距比較大也就是網格比較稀疏時，可用于全局非剛性變換，當網格間距比較小時也就是網格比較密集時可用于局部變換。同時，網格的間距越小，計算的復雜度就越高。因此在選擇網格密度時，需要權衡多方因素，在計算時間合理的前提下找到變換效果最好的網格間距。

為了能夠更好地對呼吸運動進行校正，引入塔形多層網絡結構校正的方法。首先使用分辨率較低的控制點網格對圖像進行整體校正，隨后使用分辨率高的網格對圖像進行校正：初始數據第一次經過分辨率32×32的控制點網格校正后得到第一次中間圖像數據，隨后對該數據用分辨率64×64的控制點網格進行校正并得到第二個中間圖像數據，最后用分辨率128×128的控制點網格進行最后一次校正，最終得到所需要的校正圖像。

2.2 基于CT圖像的運動特征提取

首先分析呼吸運動特性，在呼吸運動期間，肺部的運動變化呈現出周期性的特點。也就是說在一個周期內，呼吸運動相位均有差異，而在周期間這一運動過程是不斷重復的過程，類似正弦曲線的波形如圖1所示。其中橫軸t表示呼吸周期，每一個重復的弧形曲線表示一個呼吸周期，曲線上升段表示吸氣狀態此時肺部膨脹，曲線下降端表示呼吸狀態，此時肺部縮小。在每個周期的起始和末尾肺部都達到自然態，而曲線在一個呼吸周期的中點也就是達到波峰的時候肺部體積達到最大，在此處圖像面積也最大。

圖1 呼吸運動曲線圖

為了能更好達到校正的目的，正確有效地提取運動信息至關重要。在運用CT圖像提取運動信息對PET圖像進行校正之前，要解決PET圖像和CT圖像之前的時間周期匹配問題。目前CT的成像速度要大大快于PET圖像。這就造成了CT圖像每一幀的圖像要多于PET圖像，如圖2所示。本文使用時間軸對稱采樣的方法來解決這一問題。提取CT圖像中時間片時刻和PET圖像完全相同的圖像，這些圖像的位置和幾何結構與同時刻的PET圖像完全相同。事實上，這一過程可以視為對CT圖像周期數據的一個稀疏化過程。對這些CT圖像進行運動特征提取后，就可用于PET圖像呼吸運動校正。

圖2 經過周期同步后的CT和PET圖像

運動特征提取是找尋PET連續圖像中各個不同時間幀的對應圖像之間的差異，用于對這些圖像做一些補償處理使得這些圖像都能夠向基準圖像靠攏，這樣就可以去除由于運動所造成圖像偽影和圖像質量降級。通過尋找圖像間差異方式實現運動特征提取。

選取基于B樣條的FFD配準方法來獲取呼吸特性，那么這種呼吸特性的表現形式為描述呼吸周期內各個呼吸相與基準相差異的一組控制點Φi。這一過程可由以下公式表示：

其中MR表示一個B樣條配準的過程，fi表示在一個呼吸周期內的某一幀，ft表示用于配準的基準幀。根據公式（2）對周期內每一幀圖像均做配準特征提取處理后得：

最終得到Φ={Φ1，Φ2，Φ3，…}為呼吸運動一周期內的運動特性。

為了獲取運動特征信息，首先需要選定用于差異特征提取的基準幀圖像，也就是基準呼吸相。所有的圖像都將把基準幀圖像做為參考，做配準來提取運動特征（控制點集Φ）。本文使用NCAT生成的肺部體模[21]，作為研究對象，設置一個呼吸周期有11個呼吸相，也就是說用11個圖像序列來模擬一個呼吸過程。

2.3 選取基準圖像幀

先確定一個基準幀，在配準過程中這一基準幀可以是動態變化的。用B樣條獲取特征，當獲取到運動信息后通過這些控制點來對圖像做變換使其向基準靠攏。如圖3所示，本文選取第6幀作為初始基準幀，經過一次配準后圖像會作為新的基準幀與相鄰的幀來做配準。其中SF表示初始基準幀，在獲得第5幀的運動特征信息后，對第5幀圖像做B樣條變換得到第5幀圖像的動態基準幀DSF5，用于和第4幀配準，其他幀也可以以此類推，最終得到所有的運動特征信息（控制點集Φ）。

2.4 校正效果定量分析方法

在運動校正后需要對校正結果進行評估，可以通過視覺觀察圖像質量評價校正效果。但是這樣往往難以發現不同校正方式和參數在校正結果上細微的差異，而且人眼有的時候難以做到準確判斷圖像質量。本文提出使用RMCR（Respiration Motion Correction Rate）呼吸運動校正率來對校正結果進行定量分析。

圖3 動態基準方法

在校正過程中，所有的圖像切片都是與基準圖像做配準和變換的，如果對校正后圖像和基準圖像做差就可以得到一個配準后圖像與基準圖像的面積差DAIT，如公式（3）所示。其中IMST表示基準圖像，IMA表示校正后圖像。同理，對未校正圖像和基準圖像做面積差可得配準前圖像差DBIT，如公式（4）所示，其中IMB表示配準前圖像的面積。

定義RMCR未校正前圖像差和校正后圖像差的差與未校正前圖像差的比值，它描述了在校正的過程中有多少差異被校正，從理想狀態來說，校正后的圖像應該與基準圖像相同或者非常接近那么RMCR的數值就應該較大，證明大量的差異已經被校正。RMCR如公式（5）所示。因為在實際應用中要校正的圖像都是多幀組成，而每個時間幀中又包含很多圖像切片，因此可以取所有切片的平均RMCR作為對整個校正過程的效果評估，如公式（6）所示，其中M表示一個呼吸周期內的幀數，而N表示一個時間幀內所包含的圖像切片數。

3 實驗結果

實驗校正的具體步驟如圖4所示，首先對CT圖像進行同步處理。由于CT圖像與PET圖像存在采集時間上的差異，造成兩種圖像存在周期性的差異。由于CT采集圖像速度快于PET圖像，在一個相同時長的周期中CT圖像可以分割出的相位多于PET圖像。那么根據PET圖像的相位時間點來匹配CT圖像，去除CT圖像多余PET圖像中的相位圖像以實現周期同步的目的。第二步，對處理后的CT圖像進行運動特征提取。在這一過程中，首先對同一周期不同時間相位的CT圖像進行和基準幀的配準，可以得到每一個相位的CT圖像要變換成基準圖像所需的控制點集Φ。第三步，就是對所有的PET圖像進行B樣條變換，其中每一個相位的PET圖像都使用與其相對應的CT圖像的運動信息（控制點信息Φ）進行變換。經過這一系列的變換后，可以得到所需的經過運動校正的PET圖像。本文設計了兩種實驗：幾何位移形變測試和像素體模測試。

圖4 基于運動信息呼吸校正方法流程圖

3.1 幾何位移形變測試

肺部的運動是一個非常復雜的過程，簡單的使用線性位移和規則形變難以模擬肺部的復雜運動。因此，本文設計了一個結合了線性位移和形變位移的運動模型來模擬肺部運動，如圖5所示。128像素×128像素的二值化圖像中繪制一個50×50的正方形，然后加入膨脹縮小和線性位移信息。圖像融合后的一些重要信息都被膨脹信息所遮擋使得圖像效果下降。而本文的校正方法使用網格尺寸為32×32的網格進行校正，如圖6所示。校正后圖像遮擋的問題得到了解決，也沒有出現由于位移而造成的運動模糊。但是依然存在圖像角落模糊的情況，這一問題可以通過提高控制點網格分辨率的方法改善。

圖5 位移形變模型

圖6 位移形變模型

3.2 基于NCAT模擬數據的運動校正實驗結果

簡單的幾何圖像測試還不能說明本文方法的實際執行情況，因此采用更為接近真實的數據，基于NCAT體模呼吸模型對提出的方法進行驗證。NCAT可以被當作人體模型應用于PET/CT掃描仿真平臺的GATE軟件包中[22]。

NCAT是一個人體全身的模擬，由于其數據源自于真實的人體數據，所以包含較多的細節使得處理過程會很長。而本文的主要研究對象為人體的肺部運功，這里運用NCAT生成一個包含運動信息的呼吸周期圖像。在這一周期中一共有11幀圖像（其中1～5幀為吸氣7～11為呼氣），每一幀包含著31個圖像切片。每一幀圖像的各個切片都包含著相鄰幀圖像對應切片圖像的差異運動信息，就可以運用本文提出的方法（如圖4所示）進行呼吸校正。

在獲取可以用于模擬呼吸過程的運動體模信息后，就可以將該模型引入到基于GATE的PET仿真環境中，用于模擬PET掃描過程獲取包含呼吸運動的PET仿真圖像。圖7所示為一個只包含人體胸腔的NCAT體模在經過GATE仿真并還原后的圖像。

圖7 含人體胸腔的NCAT在經過GATE仿真并還原后的圖像

圖8給出了校正結果，11幀圖像的第15個切片在校正前和校正后的比較。可以發現所有的圖像都向基準幀第6幀靠攏，這樣也就達到了校正的目的。

圖8 校正實驗結果

從實驗結果來看，校正后的圖像和未校正圖像有著明顯的差異，所有圖像的形狀都接近本文所選定的初始基準第6幀圖像（作為基準幀其本身不需要校正；由于校正的時候使用了差值的方法使得在對周圍的像素起到了平均化的效果，使得圖像出現了一些花紋，這些隨著圖像分辨率的提高這一問題就會減弱）。

3.3 結果討論

根據三網格層迭代校正的實驗數據，對其進行定量分析，用本文提出的RMCR方法如公式（6），本文的實驗數據中M為11，N為31，計算這次實驗的平均呼吸運動校正率。由于DAIT和DBIT在之前校正程序中已經作為參數存在，所以可以直接在進行校正的同時獲取；可通過簡單的迭代程序實現來計算最終整體的RMCR。

經過計算本實驗的RMCR值為0.882，也就是說有10幀圖像中的每一切片和基準幀圖像切片的差異有88.2%被消除。實驗圖像為128像素×128像素，每個圖像的面積約為150～600像素，校正后大約還剩余11.8%的差異也就是17.7～70.8像素的差異（面積大的差異像素多，面積小的差異像素少），放在整體128像素×128像素的圖像中這個誤差為0.108%～0.43%之間，效果還是好的。

可以發現第5幀圖像在校正前和校正之后沒有什么差別，這是因為第5幀圖像和基準圖像十分接近（其DBIT＜10，兩幅圖像的面積差小于10個像素），在誤差允許范圍之內，沒有必要對其進行校正，對其校正可能會使圖像質量下降。

迭代校正次數與圖像質量的關系也十分密切，不是處理次數越多效果就越好，需要根據實際圖像的不同來確定迭代處理的次數，如果迭代次數過多也會降低圖像質量。本文使用的圖像數據為128像素×128像素，對其進行3次迭代校正效果較好，如果迭代4次就會出現過度校正現象。

另外，在實驗的過程中發現網格分辨率對校正圖像質量有很大的影響，對于簡單的幾何圖像，當網格分辨率較低的時候對圖像質量影響有限，但是對于復雜的不規則圖像做非剛性運動校正時，網格分辨率低會嚴重影響校正效果。將NCAT圖像分別用2×2，4×4，8×8以及32×32，64×64，128×128這幾種分辨率進行3次迭代，其效果差異十分明顯，在低網格分辨率下的非剛性運動校正出現了畸變，而在高分辨率下效果比較理想。

4 結論

針對PET/CT成像中的呼吸運動影響圖像質量問題，本文提出了基于CT圖像提取呼吸運動特征的B樣條校正方法。通過基于位移形變剛性運功實驗和基于NCAT呼吸模型的非剛性運動校正實驗，證明本文方法在控制點分辨率大小以及迭代次數合理的情況下，能夠明顯改善PET/CT呼吸運動圖像質量。

[1]許全盛，袁克虹，于麗娟，等.PET/CT圖像呼吸運動偽影校正研究進展[J].中國生物醫學工程學報，2009，28（4）：573-580.

[2]Schaerer1 J，Fassi A，Riboldi M，et al.Multi-dimensional respiratory motion tracking from markerless optical surface imaging based on deformable mesh registration[J].Physics in Medicine and Biology，2012，57：357-373.

[3]Chang G，Chang T，Pan T，et al.Implementation of an automated respiratory amplitude gating technique for PET/CT：clinical evaluation[J].Journal of Nuclear Medicine，2010，51（1）：16-24.

[4]Chamberland M，Wassenaar R，Spencer B，et al.Performance evaluation of real-time motion tracking using positron emission fiducial markers[J].Medical Physics，2011，38（2）：810-819.

[5]Kesnera A，Kuntner C.A new fast and fully automated software based algorithm for extracting respiratory signal from raw PET data and its comparison to other methods[J].Medical Physics，2010，37（10）：5550-5559.

[6]Kesner A，Bundschuh R.Respiratory gated PET derived in a fully automated manner from raw PET data[J].IEEE Transactions on Nuclear Science，2009，56（3）：677-686.

[7]Buther F，Dawood M，Stegger L，et al.List mode-driven cardiac and respiratory gating in PET[J].The Journal of Nuclear Medicine，2009，50（5）：674-681.

[8]He J，G.Keefe O，Gong S，et al.A novel method for respiratory motion gated with geometrical sensitivity of the scanner in 3D PET[J].IEEE Transactions on Nuclear Science，2008，55（5）：2557-2565.

[9]Schleyer P，Doherty M，Marsden P.Extension of a datadriven gating technique to 3D，whole body PET studies[J].Physics in Medicine and Biology，2011，56：3953-3965.

[10]Nehmeh S，Erdi Y.Respiratory motion in positron emission tomography/computed tomography：a review[J].Seminars in Nuclear Medicine，2008，38（3）：167-176.

[11]Hahn D，Daum V，Hornegger J.Automatic parameter selection for multimodal image registration[J].IEEE Transactions on Medical Imaging，2010，29（5）：1140-1155.

[12]Loeckx D，Slagmolen P，Maes F，et al.Nonrigid image registration using conditional mutual information[J].IEEE Transactions on Medical Imaging，2010，29（1）：19-29.

[13]Vandemeulebrouckea J，Rit S，Kybic J，et al.Spatiotemporal motion estimation for respiratory-correlated imaging of the lungs[J].Medical Physics，2011，38（1）：166-178.

[14]Huang X，Ren J，Guiraudon G，et al.Rapid dynamic image registration of the beating heart for diagnosis and surgical navigation[J].IEEE Transactionson MedicalImaging，2009，28（11）：1802-1814.

[15]Dawood M，Lang N，Jiang X，et al.Respiratory motion correction in 3-D PET data with advanced optical flow algorithms[J].IEEE Transactionson MedicalImaging，2008，27（8）：1164-1175.

[16]Dawood M，Brune C.A continuity equation based optical flow method for cardiac motion correction in 3D PET data[C]//Proceedings of the 5th International Workshop on Medical Imaging and Augmented Reality，2010，6326：88-97.

[17]Duijster A，De Backer S，Scheunders P.Multicomponent image restoration，an experimental study[C]//Kamel M.Proceedings of the 4th International Conference on Image Analysis and Recognition，MontrealCanada，August 22-24，2007，4633：58-68.

[18]Jiang M，Wang G，Skinner M W，et al.Blind deblurring of spiral CT images[J].IEEE Transactions on Medical Imaging，2003，22（7）：837-845.

[19]Khare A.A new method for deblurring and denoising of medical images using complex wavelet transform[C]//Proceedings of IEEE Engineering in Medicine and Biology 27th Annual Conference，2005，7：1897-1900.

[20]Herraiz J L，Samuel S，Vicente E，et al.Optimal and robust PET data sinogram restoration based on the response of thesystem[C]//ProceedingsofIEEE NuclearScience Symposium Conference，2006，6：3404-3407.

[21]Segars W P.Development and application of the new dynamic NURBS-based Cardiac-Torso（NCAT）phantom[D].The University of North Carolina at Chapel Hill，2001.

[22]Jan S，Santim G，Strul D，et al.GATE：a simulation toolkit for PET and SPECT[J].Physics in Medicine and Biology，2004，49（19）：4543-4561.