為抑制因射線源掃描時隨機抖動造成的幾何偽影，并且實現射線源掃描軌跡靈活，滿足任務驅動掃描成像要求，本文提出了一種自由軌跡錐束計算機斷層掃描（CBCT）成像的方法。本方法提出二維平面的幾何標定方法，基于該方法可以實現幾何標定模體和成像物體同時成像，進而通過在線標定方法求出幾何參數，再結合幾何參數應用于交替方向乘子法（ADMM）進行圖像迭代重建。實驗表明，本方法可以獲得較高質量的重建圖像，圖像對比度較高，特征邊緣清晰。模擬實驗結果均方根誤差（RMSE）較小，結構相似性（SSIM）值均在 0.99 以上。實際實驗結果有著較低的圖像信息熵（IE）值和較高的對比度噪聲比（CNR）值。該方法為 CBCT 實現軌跡自由以及獲得高質量重建圖像提供了一定的實用價值。
引用本文： 蔡江澤, 段曉曼, 齊宏亮, 陳宇思, 馬健暉, 周凌宏, 徐圓. 幾何標定模體與成像物同步掃描的自由軌跡錐束計算機斷層掃描成像方法. 生物醫學工程學雜志, 2021, 38(5): 951-959. doi: 10.7507/1001-5515.202101066 復制
引言
錐束計算機斷層掃描（cone-beam computed tomography，CBCT）系統具有輻射劑量低、射線利用率高、重建圖像空間分辨率高等優點[1]，現已被應用于口腔[2]、乳腺[3]和圖像引導放射治療（image guided radiation therapy，IGRT）成像[4]等方面。最常見的 CBCT 掃描軌跡為圓軌跡[5]，圓軌跡滿足濾波反投影解析求解逆問題的需求，但同時限制了進一步降低輻射劑量的可能，另外在某些特殊掃描任務時圓軌道會存在截斷偽影等成像弊端，同時實際設備研制中由于機械安裝和加工等原因有時并非是理想的圓軌跡。所以一些學者提出基于 CBCT 系統的非圓掃描軌跡，特別是針對機械臂計算機斷層掃描（computed tomography，CT）系統和 C 型臂系統，如球面正弦波軌跡[6]、三弧軌跡[5]、雙橢圓和交叉頂點軌跡[7]。與理想的圓掃描軌跡相比，符合 Orlov 條件的非圓掃描軌跡具有視場大、投影數據完整等優點，從而可以消除或者降低錐束偽影[5, 8-9]。
當前，一些靈活的非圓軌跡被設計提出，主要用于解決數字減影血管造影（digital subtraction angiography，DSA）三維（three-dimensional，3D）成像的視野（field of view，FOV）擴增及稀疏采樣問題，例如非圓軌跡被用于擴展軸向 FOV[10]和橢圓 FOV[11]，提高 3D 采樣和數據完備性[7, 12-13]，避免產生錐束偽影。另外，術前需要使用 CBCT 查看某特定的解剖結構任務如骨折部位、出血部位時，Ouadah 等[14]提出一種基于對象的個性化掃描軌跡成像，并在臨床機器人 C 型臂上實現了成像及圖像重建，該方法不僅可以降低劑量，而且提高了成像質量。因此，自由軌跡 CBCT 作為當今 CBCT 研究熱點之一，將在未來臨床中產生積極的影響。
基于此，本研究提出一種自由軌跡 CBCT 成像的方法，本方法將幾何標定模體與成像物體同時成像，提出基于平面標定模體的在線幾何標定方法求出幾何參數，再將幾何參數用于 ADMM 算法進行圖像迭代重建，實現掃描軌跡自由。標定模體上鋼珠點的遮擋會導致成像物體部分投影數據缺失，解析重建對投影完備性要求較高，若使用解析重建會導致重建圖像出現嚴重的條形偽影。而迭代重建對投影完備性要求較低，可以加入重建圖像先驗信息抑制噪聲和偽影，提高圖像質量。因此本研究采用迭代重建。文中使用的 ADMM 是一種最優化算法，結合了對偶上升法的可分解性和乘子法優良的收斂特性，可以快速解決大規模的分布式凸優化問題[15-17]，已被應用到 CT、磁共振（magnetic resonance，MR）、正電子發射斷層掃描（positron emission tomography，PET）重建等領域，本文的方法也適合其他迭代優化算法。
1 方法
對于 CBCT 的幾何標定，常規離線標定所用的 3D 圓筒形幾何標定模體是封閉式的，不易擺放，很難和成像物體同時出現在成像的 FOV 里。實驗使用的幾何標定模體為二維平面標定板（如圖 1 所示），板上嵌有鋼珠點且是開放式的，可任意擺放，更適用于同步成像任務。為了能用平面標定板上的鋼珠點求出系統的幾何參數，成像時至少要用到 2 個平面標定板。在實際 CBCT 系統中，射線源和探測器繞成像物體在 360° 范圍內掃描，對于每塊標定板，在某些掃描角度，板上的鋼珠點在投影圖像上會發生重疊。當一個板上大部分鋼珠點都在投影圖像上重疊，就無法利用該板進行幾何參數計算。故至少使用 3 塊標定板進行幾何參數求解。本實驗將三塊標定板兩兩之間呈 60° 等角度擺放。
圖1
二維平面標定板實物圖
Figure1.
Physical drawing of two-dimensional plane calibration plate
自由軌跡 CBCT 成像方法的流程圖見圖 2。首先，將幾何標定模體放置在成像物體附近并進行同步成像；對獲得的投影圖像進行鋼珠點分割并定位其質心坐標；然后，利用獲得的投影圖像上鋼珠點的質心坐標和鋼珠點在標定板坐標系下的坐標，根據基于平面標定模體的在線幾何標定方法求解得到幾何參數；最后，ADMM 算法結合幾何參數進行圖像迭代重建。
圖2
自由軌跡 CBCT 成像方法的流程圖
Figure2.
Flowchart of free trajectory CBCT imaging method
1.1 基于平面模體的在線幾何標定
下面講述如何借助于幾何標定板求解出系統的幾何參數。首先構建幾何映射模型，目的是將成像物體上的體素點與探測器上的像素點之間建立空間上的對應關系。坐標轉換示意圖如圖 3 所示。
圖3
坐標轉換示意圖
Figure3.
Schematic diagram of coordinate transformation
模體坐標系上的坐標（xw，yw，zw）經過旋轉平移后會變成射線源坐標系上的坐標（xs，ys，zs），該變化可用式（1）表示；R 為旋轉矩陣，t 為平移矩陣。若已知射線源到探測器的垂直距離（source to image distance，SID），則可以得到射線源坐標系上的坐標（xs，ys，zs）在探測器上的映射坐標（u，v），見式（2）。
 |
 |
考慮到射線源坐標系和探測器坐標系是按照平行的關系建立的，故（u，v）轉換到探測器坐標系上的坐標（ui，vi）只需加上射線源坐標系原點和探測器坐標系原點之間的平移量 u0 和 v0，見式（3）。將式（1）、（2）、（3）整合成式（4）。
 |
 |
實際中，（xw，yw，zw）為平面模體上鋼珠點的坐標，（ui，vi）為鋼珠點的投影坐標，一個鋼珠點可以形成一個等式，多個鋼珠點便可以組成方程組，用奇異值分解求解可得到幾何參數 R、t、SID、u0、v0。
1.2 幾何參數用于重建
對于基于 TV 正則化的 CBCT 圖像重建，目標函數定義為：
 |
在式（5）中，x 是待重建圖像，y 是投影數據，A 是系統矩陣。標量 μ 是平衡因子，用于調節保真項和正則項之間的權重，其取值與圖像內容的復雜性等因素相關。TV 具有良好的偽影和噪聲抑制特性，同時也能較好地保留圖像特征邊緣[18]。在目標函數中引入變量 z，通過 ADMM 方法求解得下列目標子函數[19]：
 |
 |
 |
其中，ρ > 0 為對應的増廣拉格朗日懲罰參數。交替最小化關于 x、z 的目標子函數以同時達到最優。
用共軛梯度算法求解式（6），更新圖像 x：
 |
用收縮算子精確求解式（7）：
 |
 |
在求解式（9）時，會涉及到多次的前向、反向投影的運算，需要用到射線源、探測器上投影點在成像物體坐標系下的坐標。假設射線源在射線源坐標系的坐標為（0，0，0），則根據式（1），射線源在成像物體坐標下的坐標（xs，ys，zs）見式（12）；對于待重建的體素點（x，y，z），結合式（1），其在探測器上的投影點在成像物體坐標系下的坐標（xd，yd，zd）見式（13）。這樣就將幾何參數應用于 ADMM 算法中進行圖像重建。
 |
 |
1.3 評價指標
重建圖像質量評價指標定義如下：
1）均方根誤差（root mean square errors，RMSE）：
 |
f（i，j）為求得的重建圖像在（i，j）處的值，g（i，j）為準確值。RMSE 的值越小，表明兩幅圖像差別越小。
2）結構相似性（structural similarity，SSIM）：SSIM 是一種衡量兩幅圖像相似度的指標，其公式可參見文獻[20]。SSIM 的范圍為 ? 1 到 1，SSIM 值越接近 1 說明兩幅圖像越相似。
3）圖像信息熵（image information entropy，IE）：
 |
B 為最大灰度值，H（b）和 h（b）依次為直方圖、歸一化直方圖。IE 越大，幾何偽影越明顯。
4）對比度噪聲比（contrast noise ratio，CNR）：
 |
A 和 B 代表兩個包含不同目標的感興趣區域（region of interest，ROI），E 和 σ 分別是 ROI 的灰度均值和方差。background 是指不包含主要目標物的 ROI，一般來說，CNR 越大，圖像質量越好。
2 結果
2.1 仿真體模實驗
本實驗的程序在 CUDA 編程環境和配備 2 個 NVIDIA GTX TITAN Z GPU 卡的硬件系統下開發。實驗使用數字仿真頭頸體模，利用 Siddon 的光線追蹤算法來獲取仿真投影數據。設定重建分辨率為 512 × 512 × 100，體素大小為 0.24 mm × 0.24 mm × 0.24 mm，探測器分辨率為 512 × 384，像素大小為 0.781 25 mm × 0.781 25 mm。μ 和 ρ 經手動調節后分別取 2 000、10，迭代次數為 110。實驗包含 3 種掃描軌跡（見圖 4）：① 圓；② 球面正弦波，即在圓軌跡的基礎上，射線源在旋轉軸方向上按 z = 2 × sin θ 移動，θ 為投影角度；③ 橢圓。圓和球面正弦波軌跡按等投影角間隔取 360 張投影，橢圓取 270 張投影。根據不同的掃描軌跡在 matlab 上生成相應的幾何參數 R、t、SID、u0、v0，對仿真投影進行重建。圖 5 表明，想要重建出正確的圖像需要精確匹配的幾何參數。選取重建圖像的第 81、91、93 層進行觀察。從圖 6 可以看到，3 種掃描軌跡重建圖像在特定層上與準確值相比，肉眼上已看不出任何差別。對比圖 6 紅框對應局部放大圖，重建圖像對比度均較高，特征邊緣較為清晰。
圖4
三種掃描軌跡示意圖
Figure4.
Schematic diagram of three scanning trajectories
圖5
不同軌跡的投影用圓軌跡幾何參數重建的結果
Figure5.
The projection of different trajectories reconstructed by geometric parameters of circular trajectory
圖6
仿真實驗中不同軌跡重建圖像的比較
Figure6.
Comparison of reconstructed images of different trajectories in simulation experiments
圖 7 是重建圖像在第 91 層的輪廓線對比，可以看出，3 種掃描軌跡的重建結果與準確值相差不大，圖像衰減值無明顯差異，進一步說明利用精確掃描幾何結構的幾何參數可以得到較高質量的重建圖像
圖7
仿真實驗中不同軌跡重建圖像的輪廓線比較
Figure7.
Profile comparison of reconstructed images with different trajectories in simulation experiments
同時，采用 SSIM[20]及 RMSE 兩種量化指標評價圖像整體質量，如表 1 所示。結果顯示，3 種掃描軌跡的重建圖像 RMSE 值均較小，表示與準確值相差很小，而且與準確值相比在結構上達到了 99% 的相似。
表1
仿真實驗不同軌跡重建結果定量指標分析
Table1.
Quantitative index analysis of different trajectory reconstruction results in simulation experiment
2.2 真實體模實驗
真實體模實驗采用口腔模體（見圖 8b），探測器分辨率為 192 0 × 153 6，像素大小為 0.127 mm × 0.127 mm，每隔 1° 采集一次投影，共 360 張投影。本實驗用實驗室搭建的機器人 CBCT 系統（見圖 8a）進行圓軌跡掃描。實際在機器人 CBCT 系統的安裝環節，因機械加工或安裝精度原因，實際的幾何結構和理想幾何結構出現偏差，往往導致重建圖像出現幾何偽影，因此需要進行幾何偽影校正。Rougee 等[21]提出直接線性變換（direct linear transformation，DLT）幾何標定方法，被視為 CBCT 系統幾何標定方法的金標準[22]。實驗采用 3 種方法進行重建：① 幾何標定前直接用理想幾何結構重建；② 用 DLT 進行幾何標定后重建；③ 用本文方法標定得到的幾何參數進行重建。設定重建分辨率為 512 × 512 × 512，體素大小為 0.3 mm × 0.3 mm × 0.3 mm，μ 和 ρ 均分別取 2 000、10，迭代次數為 110。選取重建圖像的第 180、219、285 層進行比較，見圖 9。可以看到，未進行幾何標定的重建圖像有很明顯的幾何偽影，出現重影模糊，結構扭曲；而經過 DLT 幾何標定后和本方法的重建圖像極為相似，均沒有明顯的幾何偽影。
圖8
本研究使用的實驗裝置
a. 機器人 CBCT 系統；b. 口腔模體
Figure8. Experimental setup used in this studya. robot CBCT system; b. dental phantom
圖9
實際實驗重建結果比較
Figure9.
Comparison of actual experimental reconstruction results
圖 10 顯示了重建結果第 180 層圖像的輪廓線比較。可以看到，在高衰減區域邊緣，未進行幾何標定得到的重建圖像未表現出高對比，而后兩種方法的重建圖像衰減值相近，且本方法的重建圖像更為真實。也可以看出，本方法在標定精度上是和金標準大體一致的。
圖10
實際實驗中重建圖像的輪廓線對比
Figure10.
Profile comparison of reconstructed images in actual experiments
使用 IE[22-23]和 CNR[24]對圖像進行定量分析，見表 2。一般來說，IE 越大，圖像幾何偽影越明顯，CNR 越大，圖像質量越好。可以看出，總體而言未進行幾何標定的重建圖像質量較差，IE 高，CNR 低；兩種方法的重建圖像質量高，沒有明顯的幾何偽影，IE 降低，CNR 升高，而且兩種方法的重建結果相差不大。
表2
實際實驗重建結果定量指標分析
Table2.
Quantitative index analysis of actual experimental reconstruction results
3 討論與結論
表 1 說明 3 種掃描軌跡的重建圖像與準確值非常接近。對比這三種軌跡的 RMSE、SSIM 值，發現正弦波軌跡和橢圓軌跡的重建質量比圓軌跡稍差。對于沿旋轉軸上成像物任一層平面，在圓軌跡掃描時，該平面在掃描時受到均勻射線的照射；而球面正弦波軌跡，由于射線源在旋轉軸方向上有 z = 2 × sin θ 的移動，導致該平面在 360° 掃描時受到的射線照射不均勻，每個角度都在改變；橢圓軌跡則是因為掃描時每個角度下射線源到旋轉中心的距離是變化的，導致該平面在 360° 掃描時受到的射線照射也是不均勻的。上述現象可能使球面正弦波軌跡和橢圓軌跡的重建質量比圓軌跡的重建質量更差。真實體模實驗表明本文提出的成像方法幾何標定精度高，與 DLT 差異不大，能有效抑制幾何偽影，得到高質量的重建圖像。圖 9 紅框對應的局部放大圖，可以看到本方法獲得的重建圖像組織結構、特征邊緣較為清晰，與標定前的重建圖像相比沒有明顯的幾何偽影。我們觀察表 2 中重建圖像第 219 層不同的方法的 IE 值，發現 DLT、本方法的 IE 值反而比標定前的值要大，與 1.3 節評價指標描述相悖。可以發現用 IE 值評判幾何偽影校正好壞的策略不一定準確。而且并沒有文獻直接表明 IE 可以作為幾何偽影的評價，其結果受所選區域影響較大。本文選取 IE 作為其中的評價指標之一，目的是測試相關指標和幾何校正效果的關系。從測試情況發現，IE 和幾何校正效果之間并非呈線性正相關。
本文提出一種幾何標定模體與成像物同步掃描的自由軌跡 CBCT 成像方法，該方法可以實現實時精確的幾何標定，有效抑制幾何偽影，獲得高質量圖像。該方法可以歸納為三步。第一，幾何標定模體和成像物同步成像；第二，利用基于平面標定模體的在線幾何標定方法求解得到幾何參數；第三，利用幾何參數進行 ADMM 迭代重建。該方法為 CBCT 實現軌跡自由以及獲得高質量重建圖像提供了一定的實用價值。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
錐束計算機斷層掃描（cone-beam computed tomography，CBCT）系統具有輻射劑量低、射線利用率高、重建圖像空間分辨率高等優點[1]，現已被應用于口腔[2]、乳腺[3]和圖像引導放射治療（image guided radiation therapy，IGRT）成像[4]等方面。最常見的 CBCT 掃描軌跡為圓軌跡[5]，圓軌跡滿足濾波反投影解析求解逆問題的需求，但同時限制了進一步降低輻射劑量的可能，另外在某些特殊掃描任務時圓軌道會存在截斷偽影等成像弊端，同時實際設備研制中由于機械安裝和加工等原因有時并非是理想的圓軌跡。所以一些學者提出基于 CBCT 系統的非圓掃描軌跡，特別是針對機械臂計算機斷層掃描（computed tomography，CT）系統和 C 型臂系統，如球面正弦波軌跡[6]、三弧軌跡[5]、雙橢圓和交叉頂點軌跡[7]。與理想的圓掃描軌跡相比，符合 Orlov 條件的非圓掃描軌跡具有視場大、投影數據完整等優點，從而可以消除或者降低錐束偽影[5, 8-9]。
當前，一些靈活的非圓軌跡被設計提出，主要用于解決數字減影血管造影（digital subtraction angiography，DSA）三維（three-dimensional，3D）成像的視野（field of view，FOV）擴增及稀疏采樣問題，例如非圓軌跡被用于擴展軸向 FOV[10]和橢圓 FOV[11]，提高 3D 采樣和數據完備性[7, 12-13]，避免產生錐束偽影。另外，術前需要使用 CBCT 查看某特定的解剖結構任務如骨折部位、出血部位時，Ouadah 等[14]提出一種基于對象的個性化掃描軌跡成像，并在臨床機器人 C 型臂上實現了成像及圖像重建，該方法不僅可以降低劑量，而且提高了成像質量。因此，自由軌跡 CBCT 作為當今 CBCT 研究熱點之一，將在未來臨床中產生積極的影響。
基于此，本研究提出一種自由軌跡 CBCT 成像的方法，本方法將幾何標定模體與成像物體同時成像，提出基于平面標定模體的在線幾何標定方法求出幾何參數，再將幾何參數用于 ADMM 算法進行圖像迭代重建，實現掃描軌跡自由。標定模體上鋼珠點的遮擋會導致成像物體部分投影數據缺失，解析重建對投影完備性要求較高，若使用解析重建會導致重建圖像出現嚴重的條形偽影。而迭代重建對投影完備性要求較低，可以加入重建圖像先驗信息抑制噪聲和偽影，提高圖像質量。因此本研究采用迭代重建。文中使用的 ADMM 是一種最優化算法，結合了對偶上升法的可分解性和乘子法優良的收斂特性，可以快速解決大規模的分布式凸優化問題[15-17]，已被應用到 CT、磁共振（magnetic resonance，MR）、正電子發射斷層掃描（positron emission tomography，PET）重建等領域，本文的方法也適合其他迭代優化算法。
1 方法
對于 CBCT 的幾何標定，常規離線標定所用的 3D 圓筒形幾何標定模體是封閉式的，不易擺放，很難和成像物體同時出現在成像的 FOV 里。實驗使用的幾何標定模體為二維平面標定板（如圖 1 所示），板上嵌有鋼珠點且是開放式的，可任意擺放，更適用于同步成像任務。為了能用平面標定板上的鋼珠點求出系統的幾何參數，成像時至少要用到 2 個平面標定板。在實際 CBCT 系統中，射線源和探測器繞成像物體在 360° 范圍內掃描，對于每塊標定板，在某些掃描角度，板上的鋼珠點在投影圖像上會發生重疊。當一個板上大部分鋼珠點都在投影圖像上重疊，就無法利用該板進行幾何參數計算。故至少使用 3 塊標定板進行幾何參數求解。本實驗將三塊標定板兩兩之間呈 60° 等角度擺放。
圖1
二維平面標定板實物圖
Figure1.
Physical drawing of two-dimensional plane calibration plate
自由軌跡 CBCT 成像方法的流程圖見圖 2。首先，將幾何標定模體放置在成像物體附近并進行同步成像；對獲得的投影圖像進行鋼珠點分割并定位其質心坐標；然后，利用獲得的投影圖像上鋼珠點的質心坐標和鋼珠點在標定板坐標系下的坐標，根據基于平面標定模體的在線幾何標定方法求解得到幾何參數；最后，ADMM 算法結合幾何參數進行圖像迭代重建。
圖2
自由軌跡 CBCT 成像方法的流程圖
Figure2.
Flowchart of free trajectory CBCT imaging method
1.1 基于平面模體的在線幾何標定
下面講述如何借助于幾何標定板求解出系統的幾何參數。首先構建幾何映射模型，目的是將成像物體上的體素點與探測器上的像素點之間建立空間上的對應關系。坐標轉換示意圖如圖 3 所示。
圖3
坐標轉換示意圖
Figure3.
Schematic diagram of coordinate transformation
模體坐標系上的坐標（xw，yw，zw）經過旋轉平移后會變成射線源坐標系上的坐標（xs，ys，zs），該變化可用式（1）表示；R 為旋轉矩陣，t 為平移矩陣。若已知射線源到探測器的垂直距離（source to image distance，SID），則可以得到射線源坐標系上的坐標（xs，ys，zs）在探測器上的映射坐標（u，v），見式（2）。
|
|
考慮到射線源坐標系和探測器坐標系是按照平行的關系建立的，故（u，v）轉換到探測器坐標系上的坐標（ui，vi）只需加上射線源坐標系原點和探測器坐標系原點之間的平移量 u0 和 v0，見式（3）。將式（1）、（2）、（3）整合成式（4）。
|
|
實際中，（xw，yw，zw）為平面模體上鋼珠點的坐標，（ui，vi）為鋼珠點的投影坐標，一個鋼珠點可以形成一個等式，多個鋼珠點便可以組成方程組，用奇異值分解求解可得到幾何參數 R、t、SID、u0、v0。
1.2 幾何參數用于重建
對于基于 TV 正則化的 CBCT 圖像重建，目標函數定義為：
|
在式（5）中，x 是待重建圖像，y 是投影數據，A 是系統矩陣。標量 μ 是平衡因子，用于調節保真項和正則項之間的權重，其取值與圖像內容的復雜性等因素相關。TV 具有良好的偽影和噪聲抑制特性，同時也能較好地保留圖像特征邊緣[18]。在目標函數中引入變量 z，通過 ADMM 方法求解得下列目標子函數[19]：
|
|
|
其中，ρ > 0 為對應的増廣拉格朗日懲罰參數。交替最小化關于 x、z 的目標子函數以同時達到最優。
用共軛梯度算法求解式（6），更新圖像 x：
|
用收縮算子精確求解式（7）：
|
|
在求解式（9）時，會涉及到多次的前向、反向投影的運算，需要用到射線源、探測器上投影點在成像物體坐標系下的坐標。假設射線源在射線源坐標系的坐標為（0，0，0），則根據式（1），射線源在成像物體坐標下的坐標（xs，ys，zs）見式（12）；對于待重建的體素點（x，y，z），結合式（1），其在探測器上的投影點在成像物體坐標系下的坐標（xd，yd，zd）見式（13）。這樣就將幾何參數應用于 ADMM 算法中進行圖像重建。
|
|
1.3 評價指標
重建圖像質量評價指標定義如下：
1）均方根誤差（root mean square errors，RMSE）：
|
f（i，j）為求得的重建圖像在（i，j）處的值，g（i，j）為準確值。RMSE 的值越小，表明兩幅圖像差別越小。
2）結構相似性（structural similarity，SSIM）：SSIM 是一種衡量兩幅圖像相似度的指標，其公式可參見文獻[20]。SSIM 的范圍為 ? 1 到 1，SSIM 值越接近 1 說明兩幅圖像越相似。
3）圖像信息熵（image information entropy，IE）：
|
B 為最大灰度值，H（b）和 h（b）依次為直方圖、歸一化直方圖。IE 越大，幾何偽影越明顯。
4）對比度噪聲比（contrast noise ratio，CNR）：
|
A 和 B 代表兩個包含不同目標的感興趣區域（region of interest，ROI），E 和 σ 分別是 ROI 的灰度均值和方差。background 是指不包含主要目標物的 ROI，一般來說，CNR 越大，圖像質量越好。
2 結果
2.1 仿真體模實驗
本實驗的程序在 CUDA 編程環境和配備 2 個 NVIDIA GTX TITAN Z GPU 卡的硬件系統下開發。實驗使用數字仿真頭頸體模，利用 Siddon 的光線追蹤算法來獲取仿真投影數據。設定重建分辨率為 512 × 512 × 100，體素大小為 0.24 mm × 0.24 mm × 0.24 mm，探測器分辨率為 512 × 384，像素大小為 0.781 25 mm × 0.781 25 mm。μ 和 ρ 經手動調節后分別取 2 000、10，迭代次數為 110。實驗包含 3 種掃描軌跡（見圖 4）：① 圓；② 球面正弦波，即在圓軌跡的基礎上，射線源在旋轉軸方向上按 z = 2 × sin θ 移動，θ 為投影角度；③ 橢圓。圓和球面正弦波軌跡按等投影角間隔取 360 張投影，橢圓取 270 張投影。根據不同的掃描軌跡在 matlab 上生成相應的幾何參數 R、t、SID、u0、v0，對仿真投影進行重建。圖 5 表明，想要重建出正確的圖像需要精確匹配的幾何參數。選取重建圖像的第 81、91、93 層進行觀察。從圖 6 可以看到，3 種掃描軌跡重建圖像在特定層上與準確值相比，肉眼上已看不出任何差別。對比圖 6 紅框對應局部放大圖，重建圖像對比度均較高，特征邊緣較為清晰。
圖4
三種掃描軌跡示意圖
Figure4.
Schematic diagram of three scanning trajectories
圖5
不同軌跡的投影用圓軌跡幾何參數重建的結果
Figure5.
The projection of different trajectories reconstructed by geometric parameters of circular trajectory
圖6
仿真實驗中不同軌跡重建圖像的比較
Figure6.
Comparison of reconstructed images of different trajectories in simulation experiments
圖 7 是重建圖像在第 91 層的輪廓線對比，可以看出，3 種掃描軌跡的重建結果與準確值相差不大，圖像衰減值無明顯差異，進一步說明利用精確掃描幾何結構的幾何參數可以得到較高質量的重建圖像
圖7
仿真實驗中不同軌跡重建圖像的輪廓線比較
Figure7.
Profile comparison of reconstructed images with different trajectories in simulation experiments
同時，采用 SSIM[20]及 RMSE 兩種量化指標評價圖像整體質量，如表 1 所示。結果顯示，3 種掃描軌跡的重建圖像 RMSE 值均較小，表示與準確值相差很小，而且與準確值相比在結構上達到了 99% 的相似。
表1
仿真實驗不同軌跡重建結果定量指標分析
Table1.
Quantitative index analysis of different trajectory reconstruction results in simulation experiment
2.2 真實體模實驗
真實體模實驗采用口腔模體（見圖 8b），探測器分辨率為 192 0 × 153 6，像素大小為 0.127 mm × 0.127 mm，每隔 1° 采集一次投影，共 360 張投影。本實驗用實驗室搭建的機器人 CBCT 系統（見圖 8a）進行圓軌跡掃描。實際在機器人 CBCT 系統的安裝環節，因機械加工或安裝精度原因，實際的幾何結構和理想幾何結構出現偏差，往往導致重建圖像出現幾何偽影，因此需要進行幾何偽影校正。Rougee 等[21]提出直接線性變換（direct linear transformation，DLT）幾何標定方法，被視為 CBCT 系統幾何標定方法的金標準[22]。實驗采用 3 種方法進行重建：① 幾何標定前直接用理想幾何結構重建；② 用 DLT 進行幾何標定后重建；③ 用本文方法標定得到的幾何參數進行重建。設定重建分辨率為 512 × 512 × 512，體素大小為 0.3 mm × 0.3 mm × 0.3 mm，μ 和 ρ 均分別取 2 000、10，迭代次數為 110。選取重建圖像的第 180、219、285 層進行比較，見圖 9。可以看到，未進行幾何標定的重建圖像有很明顯的幾何偽影，出現重影模糊，結構扭曲；而經過 DLT 幾何標定后和本方法的重建圖像極為相似，均沒有明顯的幾何偽影。
圖8
本研究使用的實驗裝置
a. 機器人 CBCT 系統；b. 口腔模體
Figure8. Experimental setup used in this studya. robot CBCT system; b. dental phantom
圖9
實際實驗重建結果比較
Figure9.
Comparison of actual experimental reconstruction results
圖 10 顯示了重建結果第 180 層圖像的輪廓線比較。可以看到，在高衰減區域邊緣，未進行幾何標定得到的重建圖像未表現出高對比，而后兩種方法的重建圖像衰減值相近，且本方法的重建圖像更為真實。也可以看出，本方法在標定精度上是和金標準大體一致的。
圖10
實際實驗中重建圖像的輪廓線對比
Figure10.
Profile comparison of reconstructed images in actual experiments
使用 IE[22-23]和 CNR[24]對圖像進行定量分析，見表 2。一般來說，IE 越大，圖像幾何偽影越明顯，CNR 越大，圖像質量越好。可以看出，總體而言未進行幾何標定的重建圖像質量較差，IE 高，CNR 低；兩種方法的重建圖像質量高，沒有明顯的幾何偽影，IE 降低，CNR 升高，而且兩種方法的重建結果相差不大。
表2
實際實驗重建結果定量指標分析
Table2.
Quantitative index analysis of actual experimental reconstruction results
3 討論與結論
表 1 說明 3 種掃描軌跡的重建圖像與準確值非常接近。對比這三種軌跡的 RMSE、SSIM 值，發現正弦波軌跡和橢圓軌跡的重建質量比圓軌跡稍差。對于沿旋轉軸上成像物任一層平面，在圓軌跡掃描時，該平面在掃描時受到均勻射線的照射；而球面正弦波軌跡，由于射線源在旋轉軸方向上有 z = 2 × sin θ 的移動，導致該平面在 360° 掃描時受到的射線照射不均勻，每個角度都在改變；橢圓軌跡則是因為掃描時每個角度下射線源到旋轉中心的距離是變化的，導致該平面在 360° 掃描時受到的射線照射也是不均勻的。上述現象可能使球面正弦波軌跡和橢圓軌跡的重建質量比圓軌跡的重建質量更差。真實體模實驗表明本文提出的成像方法幾何標定精度高，與 DLT 差異不大，能有效抑制幾何偽影，得到高質量的重建圖像。圖 9 紅框對應的局部放大圖，可以看到本方法獲得的重建圖像組織結構、特征邊緣較為清晰，與標定前的重建圖像相比沒有明顯的幾何偽影。我們觀察表 2 中重建圖像第 219 層不同的方法的 IE 值，發現 DLT、本方法的 IE 值反而比標定前的值要大，與 1.3 節評價指標描述相悖。可以發現用 IE 值評判幾何偽影校正好壞的策略不一定準確。而且并沒有文獻直接表明 IE 可以作為幾何偽影的評價，其結果受所選區域影響較大。本文選取 IE 作為其中的評價指標之一，目的是測試相關指標和幾何校正效果的關系。從測試情況發現，IE 和幾何校正效果之間并非呈線性正相關。
本文提出一種幾何標定模體與成像物同步掃描的自由軌跡 CBCT 成像方法，該方法可以實現實時精確的幾何標定，有效抑制幾何偽影，獲得高質量圖像。該方法可以歸納為三步。第一，幾何標定模體和成像物同步成像；第二，利用基于平面標定模體的在線幾何標定方法求解得到幾何參數；第三，利用幾何參數進行 ADMM 迭代重建。該方法為 CBCT 實現軌跡自由以及獲得高質量重建圖像提供了一定的實用價值。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
