在患者胸部正電子發射成像(PET)過程中,呼吸運動會造成成像質量下降。本文使用三層網格B-樣條彈性配準的方法修正重組和降噪后的正弦圖來進行呼吸運動校正。使用GATE仿真NCAT呼吸運動體模數據進行實驗,結果顯示圖像質量得到了很大地提高,運動偽影和結構信息得到了修復,表明本文的方法是有效的。
引用本文： 封碩, 崔銳, 賀建峰, 馬磊. 基于正弦圖的正電子發射成像呼吸運動彈性配準校正. 生物醫學工程學雜志, 2014, 31(5): 1005-1010. doi: 10.7507/1001-5515.20140189 復制
引言
在正電子發射成像(positron emission tomography,PET)圖像采集的過程中，人的呼吸會帶動相關的組織器官偏離原來的位置，從而會導致成像質量下降，這是長期困擾核醫學成像研究者的重要問題。與呼吸靜止態相比較，因為運動使組織位移，放射源在原位置上停留的時間變少，相同的采集時間內采集的放射源發出的準確計數就會變少，而且放射源漂移到的位置也會產生多余的計數。這些因素作用在最后的重建圖像中，產生的效果就是圖像邊緣模糊并帶有偽影，對比度下降，增加診斷難度，病灶定位不準確增加治療難度[1]。
其它的醫學成像手段雖然也會受到呼吸運動的影響，例如計算機斷層掃描(computed tomography,CT)全身掃描，數據采集時間只需3~5 s，可以通過人為指導患者呼吸就可以消除呼吸運動的影響，不會對成像質量產生影響。但是PET的數據采集需要幾分鐘，這期間要經過幾十個呼吸周期[1]，已經不能精確地人為控制患者呼吸來達到消除呼吸運動的目的，需要新的技術來校正呼吸運動對成像質量的影響。
PET經歷采集原始數據、重建、重建后圖像處理三個過程。現有的消除運動偽影的方法可以分為：基于原始數據的方法、基于圖像重建過程的方法、基于重建后圖像的方法。
基于原始數據的方法,其校正的數據源有表模式(list mode)數據或直方圖模式(histogram mode)記錄的數據。表模式以γ光子入射位置的坐標和計數時刻為信息，結合以時間刻度記錄的呼吸信號，以表格的形式記錄下來。直方圖模式將入射坐標轉換為響應線，時間信息轉換為運動狀態，將事件進行統計，便于圖像重建。Lamare等[2]通過使用仿射配準和彈性配準獲取不同運動狀態間的位移轉換信息，修正表模式數據的響應線位置，實現了呼吸運動的校正。Chamberland等[3]借助運動追蹤儀器,通過追蹤系統監測標記物在數據采集過程中的三維空間位移來求解位移轉換參數，可以精確測量物體表層形變，適合校正剛性變換和仿射變換，求得的位移轉換參數也可以校正表模式數據。Liu等[4]通過尋找體內運動和體外運動信號之間的關系，利用體外信號求取體內運動信息，用求得的運動信息修正采集的事件。直方圖模式信息更為簡潔，包含重建所必需的信息，對直方圖數據進行運動校正，雖然不能以單個事件為單位進行處理，但是運算量更小，速度更快。
基于圖像重建的方法把對圖像的運動信息加入到重建過程中，通過構造加入了運動信息參數的重建模型來達到校正運動偽影的目的。Qiao等[5]通過將配準所得的位移信息變化為插值矩陣，再用插值矩陣修正系統矩陣，最后使用迭代算法使用修正后的系統矩陣重建圖像，能夠達到較好的運動校正效果。Fin等[6]將位移信息加入到系統矩陣中，求得多個系統矩陣的平均值用于計算重建圖像。
基于重建后圖像的校正方法，主要有圖像配準和圖像恢復兩種。圖像配準的方法主要使用仿射變換配準和彈性配準兩類。圖像恢復的方法是將呼吸運動當成高頻振動造成的圖像模糊來處理，可以不用采集呼吸運動信號，直接對PET全過程采集數據的重建圖像進行處理，常用的有盲去卷積方法[7]。基于小波分解的反卷積方法也表現出優良的恢復效果[8]。
本文提出用三層網格B-樣條彈性配準的方法并通過計算求得彈性形變對采集數據的影響，是一種基于直方圖模式原始數據的校正方法，其原理是根據建立的彈性形變模型，通過樣條擬合的方法對偏離PET探測區域的數據進行猜想，從而修正運動造成的計數錯誤。
1 正弦圖配準校正方法
人體呼吸運動涉及多個器官和組織，包括肺、膈肌、心臟、肝臟等器官的形變運動[9]，是彈性的形變運動，通過縮放、旋轉和平移并不能消除重建圖中呼吸運動帶來的形變，只有彈性的圖像配準方法才適合重建后圖像的配準。觀察呼吸運動中胸腔橫切面在PET視場內的軌跡，我們發現，水平方向的形變較小，垂直方向和矢狀方向形變較大[10],而且胸腔近似橢圓形，胸腔的擴張運動也是近似以視場中心為基準且沒有旋轉的。
1.1 呼吸信號采集
在原來三維空間采集的基礎上，給采集的原始數據加上時間維度，變為四維采集。四維采集反映了三維空間影像隨時間變化的運動信息。可以從標注的時間軸上劃分不同的呼吸運動狀態[11]。目前廣泛應用的方法是,用實時位置監控系統跟蹤胸部的紅外反射標記的運動，通過識別標記軌跡來劃分呼吸運動狀態[1]。本實驗將一個完整呼吸周期劃分為9個時相，每個時相設為0.5 s。如圖 1所示，可以考慮同一時相內的采集數據，不存在呼吸運動影響；不同呼吸周期相同時相的數據，對應相同的呼吸運動狀態。
圖1
一個呼吸周期的時相劃分
Figure1.
Frames of one respiratory circle
由圖 1知，根據運動狀態采集的設置，第5個時相(g5)對應吸氣末狀態，這時候肺部充滿空氣，胸腔飽滿，重建后的圖像細節飽滿。所以，以g5作為配準的基準幀，進行B-樣條配準。我們需要處理不同時相的采集數據使它們對應相同的運動狀態，就可以保持原有計數總量，不會影響重建圖像的質量，這就是校正呼吸運動的目的。
1.2 正弦圖
PET獲得的數據以探測器位置和符合事件的計數來存儲，但這種方式并不利于原始數據參與圖像重建。通常我們用正弦圖來存儲獲取的數據，一幅二維圖像的正弦圖是一個二維數組，分別以投影角度θ和探測器對所確定直線到視場中心的徑向距離r為橫、縱坐標，以符合探測的計數值作為一個單元的取值。這種從像空間轉化到正弦圖的變換，稱為雷登變換。如圖 2所示，像空間中的一個點在正弦圖上是一條正弦曲線。重建結果圖中的一個像素的像素值是由與該像素空間位置重合的實體產生的符合計數確定，反映了該處放射性核素的分布濃度。
圖2
像空間與正弦圖關系示例
Figure2.
Relationship between image and sinogram
公式(1)表示PET圖像采集的原理，S為采集到的正弦圖，共有i個元素，代表不同θ、不同r的i對探測器組合。I為真實的人體內放射性核素濃度的分布，共有j個元素，ε為噪聲。A為系統矩陣，描述了I和S之間的雷登變換關系，有i行j列。當使用PET采集的時候，根據呼吸信號劃分k個數據集。S和I之間的關系為
| ${{S}_{k}}=A*{{I}_{k}}+{{\varepsilon }_{k}}$ |
k代表呼吸信號采集時劃分的時相個數，本文以g5為基準時相。圖像重建的過程就是Sk?Ik的過程，即式(1)的逆求解過程。然而對應不同時相的Sk求解得到的Ik也對應不同時相的運動狀態。圖像重建所期望的結果是對應相同運動狀態的。
正弦圖的每個單元的值是多個像素區域的放射性濃度的累加，將二維放射性濃度空間分布信息壓縮到特定方向的一維投影，借助濾波手段進一步平滑噪聲，使正弦圖達到配準的基本要求。這樣就可以在重建前使用原始數據進行配準求得轉換參數，也能達到運動校正的目的。
從圖 3中我們可以看到正弦圖的配準控制點(a)在像空間(b)中的分布，(a)中的每一個控制點都代表(b)中一條直線穿過的像素點的放射性濃度的積分。本文利用正弦圖中控制點網在重建后圖像像空間中的這種分布特性，提出針對以視場中心為基準的環形形變使用正弦圖配準。
重建后圖像像空間的彈性形變反映到正弦圖上是多個像素彈性形變的某一角度分量的累加，在呼吸運動的過程中不存在圖像旋轉，所以形變后正弦圖在θ軸方向幾乎無變化。生成的形變位移場描述的是r軸方向和矢狀方向的形變。
圖3
配準控制點在正弦圖和像空間中的意義
(a) 配準自動生成的控制點網;(b)控制點網在像空間中的分布
Figure3. Control points in sinogram and its distribution in image(a) control points in sinogram; (b) control points’ distribution in image
1.3 三層網格B-樣條正弦圖彈性配準
B-樣條圖像配準方法，對彈性圖像的配準效果精度高，適用性強，是目前廣泛使用的彈性圖像配準方法[12]。本文擬采用此方法對PET呼吸成像中的正弦圖數據集Sk(r,θ,z)進行配準校正，求得不同時相的位置轉換參數。這里引入自由形變模型(free form deform，FFD),認為待配準圖像是嵌入自由形變空間的浮動圖，圖像隨空間改變[13]。本文中除基準時相外的其他時相正弦圖集都是待配準浮動圖。因為運動模型不存在旋轉，默認θ無需配準。按照時相k劃分配準浮動圖集合，選取某一投影角度θ0，得到配準集合：
| $S(r,{{\theta }_{0}},z)=[{{S}_{1}}(r,{{\theta }_{0}},z),{{S}_{2}}(r,{{\theta }_{0}},z)\ldots {{S}_{k}}(r,{{\theta }_{0}},z)]$ |
以二維圖像單層網格B-樣條圖像配準為例，S5表示基圖，Sk表示浮動圖，則S5中一點r(x,y)與浮動圖中對應點,有關系為
| ${{S}_{5}}\left( r \right)={{S}_{k}}\left( r+u\left( x,y \right) \right)\text{ },$ |
其中u(x,y)表示基準圖和浮動圖之間的形變位移。于是衡量校正基準圖與浮動圖的相似性評價函數為
| $E\left( u \right)=\int \Omega {{({{S}_{5}}\left( r,u \right)-{{S}_{k}}\left( r \right))}^{2}}dr\text{ },$ |
式中Ω表示像空間，對圖像配準的問題就變成最小化評價函數argmin(E(u))問題[13]。
如果帶入完整的形變空間計算，計算量太大，通常用控制點網φi,j描述形變位移，再使用B-樣條插值的方法求得非控制點網的形變位移。如果控制點網格間距為dx和dy，用x表示取不大于x的整數值，i=x/dx-1,j=y/dy-1,s、t表示點r(x，y)在控制點網中的相對位置，s=(x/dx)-x/dx，t=(y/dy)-y/dy，用r周圍4×4個控制點為基點構造三階B-樣條曲面，控制點網的局部形變場表示為
| $L\left( r \right)=\sum\limits_{m=0}^{3}{\sum\limits_{n=0}^{3}{{{B}_{m}}\left( s \right)}}{{B}_{n}}\left( t \right){{\varphi }_{i+m,j+n}}$ |
B樣條基函數：
B0(s)=(1-s)2/6,
B1(s)=(3s3-6s2+4)/6,
B2(s)=(-3s3+3s2+3s+1)/6,
B3=(s)=s3/6,
由此得到整幅圖像的形變校正參數Tk→5(r,θ0,z)。
如果浮動圖的形變很大，會發生迭代次數超限誤差仍然沒有收斂的情況，是因為自動生成的控制點網格過密，不能很好地擬合全局形變。本文使用8×2、16×4、32×8三層網格的方法，先用稀疏網格進行全局配準，再逐步增加控制點密度，校準局部配準精度[9]，可以加快配準速度，達到更高的配準精度。完整的配準流程如圖 4所示。
圖4
配準流程圖
Figure4.
Registration flowchart
使用配準所得形變校正參數Tk→5(r,θ0,z)校正Sk(r,θ0,z),使所有數據對應g5運動狀態，得到單一角度θ0，第k時相的校正后的正弦圖數據集S′k(r,θ0,z),即
| $S\prime k(r,{{\theta }_{0}},z)={{T}_{k\to 5}}(r,{{\theta }_{0}},z)\cdot {{S}_{k}}(r,{{\theta }_{0}},z)\text{ },$ |
| $\begin{align} & S\prime (r,{{\theta }_{0}},z)=T(r,{{\theta }_{0}},z)\cdot S(r,{{\theta }_{0}},z)= \\ & [S\prime 1(r,{{\theta }_{0}},z),S\prime 2(r,{{\theta }_{0}},z)\ldots S\prime k(r,{{\theta }_{0}},z)]\text{ }, \\ \end{align}$ |
| $S\prime ={{[S\prime (r,{{\theta }_{0}},z)S\prime (r,{{\theta }_{1}},z)\ldots S\prime (r,{{\theta }_{n-1}},z)]}^{T}},$ |
式(8)中n代表θ的個數，使用S′四維正弦圖數據集重建得到的圖像就是經過運動校正的圖像。
1.4 校正評價方法
本文使用體模圖作為參考圖，采用均方誤差(mean square error,MSE)、峰值信噪比(peak signal noise ratio,PSNR)、平均結構相似度(mean structural similarity,MSSIM)三種方法來評價成像質量。MSE評價方法計算了標準圖像和待評價圖像對應像素點的灰度值誤差的均方值，該數值越小，說明兩幅圖像之間的差異越小。PSNR是計算標準圖像與待評價圖像之間的MSE相對于最大灰度值的對數，PSNR數值越大，表明圖像失真越小。MSSIM評價方法是模擬人類視覺系統模型，從圖像組成的角度來解釋圖像的結構信息，獨立于灰度和對比度，能夠反映場景中物體結構的屬性[14]。
2 實驗驗證
2.1 實驗設備及耗時
本文使用的數據是用GATE(Geant4 application for tomographic emission)[15]和NCAT(a digital NURBS based four-dimensional cardiac-torso phantom)[15]仿真呼吸成像過程。GATE是基于Geant4物理函數驅動的蒙特卡羅仿真軟件，是目前國際核醫學成像領域進行仿真實驗的重要輔助工具。用GATE 仿真GE Discovery ST PET/CT機型成像過程。使用NCAT的四維胸腔像素體模，可以準確地模擬真實人體在PET數據采集過程中的組織運動[16]。表 1列出實驗用到的計算機配置和軟件，能夠完成本實驗。使用的正弦圖分辨率128×180×32，重建圖像分辨率128×128×32。
表1
實驗環境
Table1.
Experiment environment
校正方法中計算量最大的部分為配準的迭代優化過程，使用Matlab編寫的配準代碼，設置最大迭代次數38次耗時28 364.71 s完成全部數據校正，完成某一角度θ的配準耗時平均157.52 s。經多次實驗證明，關于迭代次數和三層網格的大小設置取決于運動形變的幅度和PET空間分辨率，迭代次數和配準網格大小設置不合適會嚴重影響校正效果。
2.2 校正結果
圖 5中顯示降噪后的正弦圖,位于第15切片，未校正正弦圖是沒有劃分運動狀態進行配準的正弦圖累計圖。與運動校正正弦圖相比較，箭頭所指的位置存在明顯的運動模糊，這是由于運動造成的鄰近的響應線穿過的區域改變，有效計數分散到相鄰區域，淹沒了原來區域的特征。運動校正圖為對采集的四維正弦圖使用三次B-樣條進行配準校正后的圖像。因為本文使用配準的方法進行運動校正，對受運動影響的響應線計數進行近似插值，沒有像基于原始數據運動校正的方法一樣剔除錯誤數據，所以計數量改變較小；圖像邊界計較清晰，說明配準校正了運動帶來的模糊。
圖5
未校正和運動校正圖像對應的正弦圖比較
Figure5.
Uncorrected and motion-corrected sinograms
2.3 結果評價
針對圖 5，表 2以體模圖經過雷登變換計算的正弦圖為標準圖像，分別對未校正正弦圖、運動校正正弦圖進行評價分析。我們認為無運動正弦圖的質量是運動校正所能達到的極限質量。從表 2的評價結果來看，校正后的正弦圖與未校正正弦圖相比較，MSE減少331.443 3,相當于未校正的45.38%，說明配準校正修復了大量的計數信息，同時也反映了運動對圖像的計數值改變是巨大的；PSNR提升2.626 3，說明配準校正修復了大部分運動模糊；MSSIM提升0.0529，整幅圖像綜合質量已有很大的提升。本文計算的MSSIM為待評價圖像與體模正弦圖圖像的MSSIM之比。
表2
評價結果
Table2.
Result of evaluation
將圖 5所示正弦圖進行重建，我們所采用的迭代重建方法是最大似然期望最大化(maximum likehood expectation maximuzation,MLEM)算法[17]。重建結果如圖 6所示，這兩幅圖分別對應圖 5所示正弦圖的重建圖，使用MLEM迭代算法，15次迭代的重建結果。圖中箭頭所指的位置是運動形變嚴重的部位，可以看到未校正圖和無運動圖相比較有較大區別，箭頭所指區域的空間位移明顯。經過正弦圖配準校正后，形變得到較好的修正。從重建結果圖來看，經過正弦圖配準校正，重建圖像的運動偽影得到修正，邊界更加清晰。
圖6
未校正和校正后正弦圖的重建圖像比較
Figure6.
Reconstruction images of uncorrected and corrected sinograms
3 結論
本文使用B-樣條方法對PET呼吸成像過程中的正弦圖進行呼吸運動校正。實驗結果表明基于正弦圖的彈性配準校正呼吸運動的方法能夠提高重建圖像的對比度，修正模糊邊緣，削減偽影，得到較好的校正結果。我們也采用MSE、PSNR、MSSIM等方法評價實驗結果，顯示圖像質量有明顯改善，逼近靜態采集的成像效果，能夠達到期望的校正效果。
引言
在正電子發射成像(positron emission tomography,PET)圖像采集的過程中，人的呼吸會帶動相關的組織器官偏離原來的位置，從而會導致成像質量下降，這是長期困擾核醫學成像研究者的重要問題。與呼吸靜止態相比較，因為運動使組織位移，放射源在原位置上停留的時間變少，相同的采集時間內采集的放射源發出的準確計數就會變少，而且放射源漂移到的位置也會產生多余的計數。這些因素作用在最后的重建圖像中，產生的效果就是圖像邊緣模糊并帶有偽影，對比度下降，增加診斷難度，病灶定位不準確增加治療難度[1]。
其它的醫學成像手段雖然也會受到呼吸運動的影響，例如計算機斷層掃描(computed tomography,CT)全身掃描，數據采集時間只需3~5 s，可以通過人為指導患者呼吸就可以消除呼吸運動的影響，不會對成像質量產生影響。但是PET的數據采集需要幾分鐘，這期間要經過幾十個呼吸周期[1]，已經不能精確地人為控制患者呼吸來達到消除呼吸運動的目的，需要新的技術來校正呼吸運動對成像質量的影響。
PET經歷采集原始數據、重建、重建后圖像處理三個過程。現有的消除運動偽影的方法可以分為：基于原始數據的方法、基于圖像重建過程的方法、基于重建后圖像的方法。
基于原始數據的方法,其校正的數據源有表模式(list mode)數據或直方圖模式(histogram mode)記錄的數據。表模式以γ光子入射位置的坐標和計數時刻為信息，結合以時間刻度記錄的呼吸信號，以表格的形式記錄下來。直方圖模式將入射坐標轉換為響應線，時間信息轉換為運動狀態，將事件進行統計，便于圖像重建。Lamare等[2]通過使用仿射配準和彈性配準獲取不同運動狀態間的位移轉換信息，修正表模式數據的響應線位置，實現了呼吸運動的校正。Chamberland等[3]借助運動追蹤儀器,通過追蹤系統監測標記物在數據采集過程中的三維空間位移來求解位移轉換參數，可以精確測量物體表層形變，適合校正剛性變換和仿射變換，求得的位移轉換參數也可以校正表模式數據。Liu等[4]通過尋找體內運動和體外運動信號之間的關系，利用體外信號求取體內運動信息，用求得的運動信息修正采集的事件。直方圖模式信息更為簡潔，包含重建所必需的信息，對直方圖數據進行運動校正，雖然不能以單個事件為單位進行處理，但是運算量更小，速度更快。
基于圖像重建的方法把對圖像的運動信息加入到重建過程中，通過構造加入了運動信息參數的重建模型來達到校正運動偽影的目的。Qiao等[5]通過將配準所得的位移信息變化為插值矩陣，再用插值矩陣修正系統矩陣，最后使用迭代算法使用修正后的系統矩陣重建圖像，能夠達到較好的運動校正效果。Fin等[6]將位移信息加入到系統矩陣中，求得多個系統矩陣的平均值用于計算重建圖像。
基于重建后圖像的校正方法，主要有圖像配準和圖像恢復兩種。圖像配準的方法主要使用仿射變換配準和彈性配準兩類。圖像恢復的方法是將呼吸運動當成高頻振動造成的圖像模糊來處理，可以不用采集呼吸運動信號，直接對PET全過程采集數據的重建圖像進行處理，常用的有盲去卷積方法[7]。基于小波分解的反卷積方法也表現出優良的恢復效果[8]。
本文提出用三層網格B-樣條彈性配準的方法并通過計算求得彈性形變對采集數據的影響，是一種基于直方圖模式原始數據的校正方法，其原理是根據建立的彈性形變模型，通過樣條擬合的方法對偏離PET探測區域的數據進行猜想，從而修正運動造成的計數錯誤。
1 正弦圖配準校正方法
人體呼吸運動涉及多個器官和組織，包括肺、膈肌、心臟、肝臟等器官的形變運動[9]，是彈性的形變運動，通過縮放、旋轉和平移并不能消除重建圖中呼吸運動帶來的形變，只有彈性的圖像配準方法才適合重建后圖像的配準。觀察呼吸運動中胸腔橫切面在PET視場內的軌跡，我們發現，水平方向的形變較小，垂直方向和矢狀方向形變較大[10],而且胸腔近似橢圓形，胸腔的擴張運動也是近似以視場中心為基準且沒有旋轉的。
1.1 呼吸信號采集
在原來三維空間采集的基礎上，給采集的原始數據加上時間維度，變為四維采集。四維采集反映了三維空間影像隨時間變化的運動信息。可以從標注的時間軸上劃分不同的呼吸運動狀態[11]。目前廣泛應用的方法是,用實時位置監控系統跟蹤胸部的紅外反射標記的運動，通過識別標記軌跡來劃分呼吸運動狀態[1]。本實驗將一個完整呼吸周期劃分為9個時相，每個時相設為0.5 s。如圖 1所示，可以考慮同一時相內的采集數據，不存在呼吸運動影響；不同呼吸周期相同時相的數據，對應相同的呼吸運動狀態。
圖1
一個呼吸周期的時相劃分
Figure1.
Frames of one respiratory circle
由圖 1知，根據運動狀態采集的設置，第5個時相(g5)對應吸氣末狀態，這時候肺部充滿空氣，胸腔飽滿，重建后的圖像細節飽滿。所以，以g5作為配準的基準幀，進行B-樣條配準。我們需要處理不同時相的采集數據使它們對應相同的運動狀態，就可以保持原有計數總量，不會影響重建圖像的質量，這就是校正呼吸運動的目的。
1.2 正弦圖
PET獲得的數據以探測器位置和符合事件的計數來存儲，但這種方式并不利于原始數據參與圖像重建。通常我們用正弦圖來存儲獲取的數據，一幅二維圖像的正弦圖是一個二維數組，分別以投影角度θ和探測器對所確定直線到視場中心的徑向距離r為橫、縱坐標，以符合探測的計數值作為一個單元的取值。這種從像空間轉化到正弦圖的變換，稱為雷登變換。如圖 2所示，像空間中的一個點在正弦圖上是一條正弦曲線。重建結果圖中的一個像素的像素值是由與該像素空間位置重合的實體產生的符合計數確定，反映了該處放射性核素的分布濃度。
圖2
像空間與正弦圖關系示例
Figure2.
Relationship between image and sinogram
公式(1)表示PET圖像采集的原理，S為采集到的正弦圖，共有i個元素，代表不同θ、不同r的i對探測器組合。I為真實的人體內放射性核素濃度的分布，共有j個元素，ε為噪聲。A為系統矩陣，描述了I和S之間的雷登變換關系，有i行j列。當使用PET采集的時候，根據呼吸信號劃分k個數據集。S和I之間的關系為
| ${{S}_{k}}=A*{{I}_{k}}+{{\varepsilon }_{k}}$ |
k代表呼吸信號采集時劃分的時相個數，本文以g5為基準時相。圖像重建的過程就是Sk?Ik的過程，即式(1)的逆求解過程。然而對應不同時相的Sk求解得到的Ik也對應不同時相的運動狀態。圖像重建所期望的結果是對應相同運動狀態的。
正弦圖的每個單元的值是多個像素區域的放射性濃度的累加，將二維放射性濃度空間分布信息壓縮到特定方向的一維投影，借助濾波手段進一步平滑噪聲，使正弦圖達到配準的基本要求。這樣就可以在重建前使用原始數據進行配準求得轉換參數，也能達到運動校正的目的。
從圖 3中我們可以看到正弦圖的配準控制點(a)在像空間(b)中的分布，(a)中的每一個控制點都代表(b)中一條直線穿過的像素點的放射性濃度的積分。本文利用正弦圖中控制點網在重建后圖像像空間中的這種分布特性，提出針對以視場中心為基準的環形形變使用正弦圖配準。
重建后圖像像空間的彈性形變反映到正弦圖上是多個像素彈性形變的某一角度分量的累加，在呼吸運動的過程中不存在圖像旋轉，所以形變后正弦圖在θ軸方向幾乎無變化。生成的形變位移場描述的是r軸方向和矢狀方向的形變。
圖3
配準控制點在正弦圖和像空間中的意義
(a) 配準自動生成的控制點網;(b)控制點網在像空間中的分布
Figure3. Control points in sinogram and its distribution in image(a) control points in sinogram; (b) control points’ distribution in image
1.3 三層網格B-樣條正弦圖彈性配準
B-樣條圖像配準方法，對彈性圖像的配準效果精度高，適用性強，是目前廣泛使用的彈性圖像配準方法[12]。本文擬采用此方法對PET呼吸成像中的正弦圖數據集Sk(r,θ,z)進行配準校正，求得不同時相的位置轉換參數。這里引入自由形變模型(free form deform，FFD),認為待配準圖像是嵌入自由形變空間的浮動圖，圖像隨空間改變[13]。本文中除基準時相外的其他時相正弦圖集都是待配準浮動圖。因為運動模型不存在旋轉，默認θ無需配準。按照時相k劃分配準浮動圖集合，選取某一投影角度θ0，得到配準集合：
| $S(r,{{\theta }_{0}},z)=[{{S}_{1}}(r,{{\theta }_{0}},z),{{S}_{2}}(r,{{\theta }_{0}},z)\ldots {{S}_{k}}(r,{{\theta }_{0}},z)]$ |
以二維圖像單層網格B-樣條圖像配準為例，S5表示基圖，Sk表示浮動圖，則S5中一點r(x,y)與浮動圖中對應點,有關系為
| ${{S}_{5}}\left( r \right)={{S}_{k}}\left( r+u\left( x,y \right) \right)\text{ },$ |
其中u(x,y)表示基準圖和浮動圖之間的形變位移。于是衡量校正基準圖與浮動圖的相似性評價函數為
| $E\left( u \right)=\int \Omega {{({{S}_{5}}\left( r,u \right)-{{S}_{k}}\left( r \right))}^{2}}dr\text{ },$ |
式中Ω表示像空間，對圖像配準的問題就變成最小化評價函數argmin(E(u))問題[13]。
如果帶入完整的形變空間計算，計算量太大，通常用控制點網φi,j描述形變位移，再使用B-樣條插值的方法求得非控制點網的形變位移。如果控制點網格間距為dx和dy，用x表示取不大于x的整數值，i=x/dx-1,j=y/dy-1,s、t表示點r(x，y)在控制點網中的相對位置，s=(x/dx)-x/dx，t=(y/dy)-y/dy，用r周圍4×4個控制點為基點構造三階B-樣條曲面，控制點網的局部形變場表示為
| $L\left( r \right)=\sum\limits_{m=0}^{3}{\sum\limits_{n=0}^{3}{{{B}_{m}}\left( s \right)}}{{B}_{n}}\left( t \right){{\varphi }_{i+m,j+n}}$ |
B樣條基函數：
B0(s)=(1-s)2/6,
B1(s)=(3s3-6s2+4)/6,
B2(s)=(-3s3+3s2+3s+1)/6,
B3=(s)=s3/6,
由此得到整幅圖像的形變校正參數Tk→5(r,θ0,z)。
如果浮動圖的形變很大，會發生迭代次數超限誤差仍然沒有收斂的情況，是因為自動生成的控制點網格過密，不能很好地擬合全局形變。本文使用8×2、16×4、32×8三層網格的方法，先用稀疏網格進行全局配準，再逐步增加控制點密度，校準局部配準精度[9]，可以加快配準速度，達到更高的配準精度。完整的配準流程如圖 4所示。
圖4
配準流程圖
Figure4.
Registration flowchart
使用配準所得形變校正參數Tk→5(r,θ0,z)校正Sk(r,θ0,z),使所有數據對應g5運動狀態，得到單一角度θ0，第k時相的校正后的正弦圖數據集S′k(r,θ0,z),即
| $S\prime k(r,{{\theta }_{0}},z)={{T}_{k\to 5}}(r,{{\theta }_{0}},z)\cdot {{S}_{k}}(r,{{\theta }_{0}},z)\text{ },$ |
| $\begin{align} & S\prime (r,{{\theta }_{0}},z)=T(r,{{\theta }_{0}},z)\cdot S(r,{{\theta }_{0}},z)= \\ & [S\prime 1(r,{{\theta }_{0}},z),S\prime 2(r,{{\theta }_{0}},z)\ldots S\prime k(r,{{\theta }_{0}},z)]\text{ }, \\ \end{align}$ |
| $S\prime ={{[S\prime (r,{{\theta }_{0}},z)S\prime (r,{{\theta }_{1}},z)\ldots S\prime (r,{{\theta }_{n-1}},z)]}^{T}},$ |
式(8)中n代表θ的個數，使用S′四維正弦圖數據集重建得到的圖像就是經過運動校正的圖像。
1.4 校正評價方法
本文使用體模圖作為參考圖，采用均方誤差(mean square error,MSE)、峰值信噪比(peak signal noise ratio,PSNR)、平均結構相似度(mean structural similarity,MSSIM)三種方法來評價成像質量。MSE評價方法計算了標準圖像和待評價圖像對應像素點的灰度值誤差的均方值，該數值越小，說明兩幅圖像之間的差異越小。PSNR是計算標準圖像與待評價圖像之間的MSE相對于最大灰度值的對數，PSNR數值越大，表明圖像失真越小。MSSIM評價方法是模擬人類視覺系統模型，從圖像組成的角度來解釋圖像的結構信息，獨立于灰度和對比度，能夠反映場景中物體結構的屬性[14]。
2 實驗驗證
2.1 實驗設備及耗時
本文使用的數據是用GATE(Geant4 application for tomographic emission)[15]和NCAT(a digital NURBS based four-dimensional cardiac-torso phantom)[15]仿真呼吸成像過程。GATE是基于Geant4物理函數驅動的蒙特卡羅仿真軟件，是目前國際核醫學成像領域進行仿真實驗的重要輔助工具。用GATE 仿真GE Discovery ST PET/CT機型成像過程。使用NCAT的四維胸腔像素體模，可以準確地模擬真實人體在PET數據采集過程中的組織運動[16]。表 1列出實驗用到的計算機配置和軟件，能夠完成本實驗。使用的正弦圖分辨率128×180×32，重建圖像分辨率128×128×32。
表1
實驗環境
Table1.
Experiment environment
校正方法中計算量最大的部分為配準的迭代優化過程，使用Matlab編寫的配準代碼，設置最大迭代次數38次耗時28 364.71 s完成全部數據校正，完成某一角度θ的配準耗時平均157.52 s。經多次實驗證明，關于迭代次數和三層網格的大小設置取決于運動形變的幅度和PET空間分辨率，迭代次數和配準網格大小設置不合適會嚴重影響校正效果。
2.2 校正結果
圖 5中顯示降噪后的正弦圖,位于第15切片，未校正正弦圖是沒有劃分運動狀態進行配準的正弦圖累計圖。與運動校正正弦圖相比較，箭頭所指的位置存在明顯的運動模糊，這是由于運動造成的鄰近的響應線穿過的區域改變，有效計數分散到相鄰區域，淹沒了原來區域的特征。運動校正圖為對采集的四維正弦圖使用三次B-樣條進行配準校正后的圖像。因為本文使用配準的方法進行運動校正，對受運動影響的響應線計數進行近似插值，沒有像基于原始數據運動校正的方法一樣剔除錯誤數據，所以計數量改變較小；圖像邊界計較清晰，說明配準校正了運動帶來的模糊。
圖5
未校正和運動校正圖像對應的正弦圖比較
Figure5.
Uncorrected and motion-corrected sinograms
2.3 結果評價
針對圖 5，表 2以體模圖經過雷登變換計算的正弦圖為標準圖像，分別對未校正正弦圖、運動校正正弦圖進行評價分析。我們認為無運動正弦圖的質量是運動校正所能達到的極限質量。從表 2的評價結果來看，校正后的正弦圖與未校正正弦圖相比較，MSE減少331.443 3,相當于未校正的45.38%，說明配準校正修復了大量的計數信息，同時也反映了運動對圖像的計數值改變是巨大的；PSNR提升2.626 3，說明配準校正修復了大部分運動模糊；MSSIM提升0.0529，整幅圖像綜合質量已有很大的提升。本文計算的MSSIM為待評價圖像與體模正弦圖圖像的MSSIM之比。
表2
評價結果
Table2.
Result of evaluation
將圖 5所示正弦圖進行重建，我們所采用的迭代重建方法是最大似然期望最大化(maximum likehood expectation maximuzation,MLEM)算法[17]。重建結果如圖 6所示，這兩幅圖分別對應圖 5所示正弦圖的重建圖，使用MLEM迭代算法，15次迭代的重建結果。圖中箭頭所指的位置是運動形變嚴重的部位，可以看到未校正圖和無運動圖相比較有較大區別，箭頭所指區域的空間位移明顯。經過正弦圖配準校正后，形變得到較好的修正。從重建結果圖來看，經過正弦圖配準校正，重建圖像的運動偽影得到修正，邊界更加清晰。
圖6
未校正和校正后正弦圖的重建圖像比較
Figure6.
Reconstruction images of uncorrected and corrected sinograms
3 結論
本文使用B-樣條方法對PET呼吸成像過程中的正弦圖進行呼吸運動校正。實驗結果表明基于正弦圖的彈性配準校正呼吸運動的方法能夠提高重建圖像的對比度，修正模糊邊緣，削減偽影，得到較好的校正結果。我們也采用MSE、PSNR、MSSIM等方法評價實驗結果，顯示圖像質量有明顯改善，逼近靜態采集的成像效果，能夠達到期望的校正效果。
