雙能計算機斷層成像（CT）技術是 CT 成像領域未來重要的發展方向。雙能 CT 重建算法主流的模型是基物質分解模型，算法的核心關鍵是求出基物質分解系數投影值。基于投影匹配的雙能 CT 投影分解算法通過建立能譜查找表，使用最小二乘法進行匹配查找得到分解系數投影值。但該方法由于查找表數據龐大，計算時間長，不利于臨床的應用。本文在該方法的基礎上，提出一種通過直線方程擬合和平面方程擬合查找表數據，快速計算分解系數投影值的改進算法。仿真實驗證明，該算法在大幅提高計算速度的同時，也能穩定地收斂到正確的解。
引用本文： 侯曉文, 盧子鵬, 滕月陽, 孝大宇, 范晟昱, 楊超然, 劉瑜伽, 康雁. 基于投影匹配的雙能計算機斷層成像投影分解加速算法. 生物醫學工程學雜志, 2018, 35(3): 376-383. doi: 10.7507/1001-5515.201605009 復制
引言
計算機斷層成像（computed tomography，CT）所用的 X 射線能譜包含不同能量的 X 射線，就如白光包含不同顏色的光。傳統 CT 成像技術把 X 射線能譜簡化為等價的單一能量的 X 射線，利用不同物質對此單一能量 X 射線的衰減信息不同得到掃描組織的形態學結構圖像。雙能 CT 成像技術不對 X 射線能譜進行簡化，能同時得到掃描組織的形態學結構圖像和化學組織成分圖像，進行化學組織成分分析，這是傳統 CT 成像技術無法企及的優點[1]。雙能 CT 成像技術的關鍵是雙能 CT 重建算法。重建算法的模型主要有雙效應分解模型和基物質分解模型。本文主要討論的是基物質分解模型。
根據重建方法的不同，雙能 CT 重建算法可以分為三大類[2]：① 投影分解重建算法，該類方法主要是通過高、低能投影數據得到分解系數投影值，然后由投影值反投影求出分解系數。② 后處理重建算法，該類方法首先使用高、低能投影數據重建出高、低能 CT 圖像，然后使用高、低能 CT 圖像求解分解系數，這種方法精確度不高。③ 迭代重建算法，通過雙能迭代，得到分解系數圖像。無論哪種方法，最終目的都是求出分解系數。李保磊、張耀軍在 2011 年提出基于投影匹配的雙能 CT 投影分解重建算法[1]，該算法屬于第一類。依據 CT 系統能譜和基物質線性衰減曲線，通過投影積分方程組建立 X 射線的高低能譜查找表，利用沿著某條 X 射線經過的路徑下實測得到的高低能投影數據，在查找表中尋找其最佳匹配點，該最佳匹配點就是對應的基物質分解系數的投影值，然后通過傳統的 CT 重建算法如濾波反投影重建算法（filter back projection，FBP）[3]重建求出對應的分解系數圖像。但實際應用中，由于查找表數據龐大，在尋找最佳匹配點的過程中需要遍歷整張查找表，計算過程十分耗時，不利于臨床應用。該算法的計算量大部分在于尋找最佳匹配點。針對這一點，本文在該算法基礎上提出一種加速算法。
1 方法
1.1 雙能 CT 重建算法的原理[4 -8 ]
雙能 CT 重建算法模型主要是雙效應分解模型和基物質分解模型：
 |
 |
u（E）代表物質在能量 E 下的線性衰減系數，式（1）對應的是雙效應分解，在一定的射線能量范圍內，物質的線性衰減可由光電效應和康普頓散射兩種作用共同組成。其中 fp（E）、fkn（E）只與能量 E 有關而與材質無關，ac、ap 分別對應的是光電效應和康普頓效應的分解系數。式（2）對應基物質分解模型，表示任何一種物質的線性衰減系數都可由兩種基物質的線性衰減系數線性表示。其中，u1（E）、u2（E）是兩種基物質在能量 E 下對應的線性衰減值，該數據可以在美國國家標準與技術研究院（National Institute of Standards and Technology，NIST）庫中查詢得到，b1、b2 是兩種基物質對應的分解系數。在本文中，我們主要是討論基物質分解模型。由式（2）可知，在某兩個單能量如 E1、E2 下，就可以得到兩個方程，其中 u（E）、u1（E）、u2（E）是已知數，而 b1、b2 是未知數，通過兩個方程組求解即可得到分解系數 b1、b2。在實際中，球管產生的 X 射線并不是單能譜射線，而是連續的混合能譜射線，因此通過式（2）無法求出基物質的分解系數。
根據基物質分解模型，在沿著某一條 X 射線的路徑下有：
 |
B1、B2 分別表示兩種基物質的分解系數投影值，根據混合能譜 X 射線條件下的朗伯比爾定律可得：
 |
其中，
代表沿著某條 X 射線路徑下高低能投影值，
為 X 射線的高低能譜數據。由式（3）可知，只要求出基物質分解系數投影值 B1、B2，再通過 FBP 即可得到基物質分解系數圖像。根據式（4）和式（5），我們可以求出重建圖像中每個像素點對應的有效原子序數 Zeff 和有效電子密度 ρe：
 |
 |
其中 ρe1 和 ρe2 是基物質的有效電子密度，Z1 和 Z2 是基物質的有效原子序數。n 的值為 3.5。綜上所述，如何計算分解系數投影值就成為雙能 CT 重建算法的難點。基于投影匹配的雙能 CT 投影分解重建算法就是為了計算分解系數投影值而提出的。
1.2 基于投影匹配的雙能 CT 投影分解重建算法[1 ]
基于投影匹配的雙能 CT 投影分解重建算法的主要思想是：根據 CT 系統能譜信息生成合理意義范圍內的一系列 B1 和 B2 數值，再根據式（3）來生成對應的高低能投影數值 PH 和 PL，建立查找表數據，然后根據給定的實測高低能投影數據
、
，匹配查找表，尋找最佳匹配點，即分解系數投影值
和
。主要包含三個步驟：
（1）設定 B1 和 B2 的范圍，
，
，根據設定的步進 d 來產生 B1（i）和 B2（j）序列。根據系統能譜數據和式（3）計算對應的高低能投影值 PH 和 PL。我們用 i 來代表查找表的行編號，j 代表查找表的列編號，則有：
 |
（2）根據實測高低能投影數據
和
，利用最小二乘法，在查找表中進行匹配比較，使得式（6）的值最小化。
 |
（3）根據最佳匹配點的坐標（i0，j0）來計算分解系數投影值：
 |
只要得到分解系數投影值，通過 FBP 就能求出基物質分解系數圖像，進而重建單能量圖像，擬合特定區域物質的線性衰減曲線以及重建有效電子密度圖像和有效原子序數圖像。該算法避開了復雜的數學公式，實現起來很簡單。但為了提高計算精度，需要把查找表的步進設定成非常小的數值，同時 B1、B2 范圍比較大，因此查找表非常龐大，查找最佳匹配點也就非常耗時。為了解決這個問題，本文提出了一種加速算法。
1.3 基于投影匹配的雙能 CT 成像投影分解加速算法
本文提出一種查找最佳匹配點的加速算法，其基本思想是：對于一對實測的高低能投影數據，通過某種方法確定其在查找表中對應的最佳匹配點所在的大致區域，然后對該區域進行遍歷、匹配和比較，找到最佳匹配點。該方法可以總結為以下兩種算法。
1.3.1 基于直線擬合查找表數據的加速投影分解算法
以 C 和 Fe 為基物質的查找表為例，把查找表中某一行數據取出來畫圖，如圖 1 所示。
從圖中可以看出，高低能投影值關于變量 j 呈現近似線性關系。根據這個規律，可以建立直線方程來擬合該行查找表數據[9]。
 |
 |
其中 PL（i，j）和 PH（i，j）表示查找表中坐標為（i，j）的數據點的低能和高能投影數值。aL、bL、aH、bH 分別是直線方程系數。式（6）用 f（j）來表示，把式（8）和式（9）代入得到：
 |
圖1
某一行 C-Fe 查找表數據畫圖顯示
Figure1.
The drawing display of a line of C-Fe lookup table data
要尋找最佳匹配點，必須最小化 f（j）的值。f（j）是關于變量 j 的一元二次函數，同時也是開口向上的凸函數，其最小值要么在極小值點取得，要么在邊界點取得。對變量 j 進行求導，并令導數等于 0，可以得到極小值點：
 |
若在極小值點取得最小值，那么最佳匹配點就在極值點附近。因此設定一個閾值 th，如圖 2 所示，遍歷區間
中每一個點，找到最佳匹配點。
圖2
某一行查找表數據一維表示
Figure2.
One-dimensional representation of a row of lookup table data
若 j0 < jmin 或者 j0 > jmax，函數 f（j）的最小值只能在邊界點 jmin 和 jmax 取得。
對于一組給定的實測高低能投影數據，上述方法只能獲取某一行數據的最佳匹配點，稱為行最佳匹配點。使用同樣的方法可以求出其他行的行最佳匹配點，最后比較所有的行最佳匹配點，找到全局最佳匹配點。該方法需要遍歷查找表的每一行，但避免了遍歷行內所有數據點，大大減少了計算量。
1.3.2 基于平面擬合查找表數據的加速投影分解算法
基于平面擬合查找表數據的加速投影分解算法用平面方程擬合查找表數據。通過建立目標函數方程來求解該高低能投影數據對應于查找表的最佳匹配點，求出對應的分解系數的投影值。圖 3 是高低能查找表三維空間顯示結果。由圖可見，高低能查找表在三維空間近似一個平面，可以通過建立平面方程[10]擬合高低能查找表。
 |
 |
式（12）和式（13）分別表示擬合低高能投影數據的平面方程，AL、BL、CL、AH、BH、CH 都是平面方程的系數，把式（12）和式（13）代入式（6）可得：
 |
圖3
C-Fe 查找表數據在三維空間顯示
Figure3.
C-Fe lookup table data are displayed in three-dimensional space
由式（14）可知，f（i，j）是關于變量 i 和 j 的開口向上的凸函數，在區間
&&
內存在極小值點。對變量 i 和變量 j 求偏導：
 |
為了簡化式子，用下面符號代表不同的方程系數：
 |
則可以得到下面兩個方程組：
 |
兩個方程組聯立求解最終得到 i0 和 j0：
 |
 |
若極小值點（i0，j0）是最小值點，最佳匹配點一定在極小值點（i0，j0）附近。因此可以設定一個閾值 th，在區間
&&
內每個數據點進行匹配計算，找到最佳匹配點。
若極小值點（i0，j0）不在設定的區間范圍內，那么最小值在查找表的邊界上取得。對應邊界區間分為四個：
 |
 |
 |
 |
在邊界區間中每個數據點進行匹配計算找到最佳匹配點，從而求出實際投影數據對應的分解系數投影值。基于平面擬合查找表數據加速投影分解算法不需要遍歷查找表的每一行，只要給出實測的高低能投影數據，就能定位對應的最佳匹配點的大致位置。同時，該方法也需要設定一個閾值，在一個特定區域遍歷和匹配計算，比較每一個數據點，才能找到最佳匹配點。
總結來說，兩種方法的思路都是先進行粗糙搜索，定位最佳匹配點的大致區域，然后進行精確搜索，找到最佳匹配點，計算得到分解系數投影值。兩種方法各有優缺點。對于列的數目遠遠多于行的數目的查找表，使用基于直線擬合查找表數據加速投影分解算法可以很快得到解。對于行列數相差不多的查找表，基于平面擬合查找表數據加速投影分解算法比較有優勢。但是為了保證準確性，基于平面擬合查找表數據加速投影分解算法的閾值 th 一般比較大，遍歷區域也就比較大，這會增加計算量。可以對這兩種方法進行混合使用，首先通過平面加速算法確定遍歷區域，然后在此遍歷區域上使用直線加速算法確定最佳匹配點，就可以起到很好的加速效果。
采用上述加速算法計算得到基物質分解系數投影值以后，通過傳統的 CT 重建算法如 FB 算法重建出對應的基物質分解系數圖像，然后就可以重建單能量圖像、有效電子密度圖像和有效原子序數圖像。
2 仿真實驗
為了敘述簡單，我們把基于投影匹配的雙能 CT 投影分解重建算法簡稱為基礎算法（basic algorithm），基于直線擬合查找表數據的加速投影分解算法簡稱為直線加速算法（linear acceleration algorithm），基于平面擬合查找表數據的加速投影分解算法簡稱為平面加速算法（plane acceleration algorithm），直線加速算法和平面加速算法混合使用則簡稱為直線-平面混合加速算法（the mixture algorithm of line-plane acceleration）。
2.1 產生仿真數據
根據已有的論文，設計體模圖像[11]，圖像的大小為 256 × 256 像素，空間分辨率為 1 mm，如圖 4 所示。
圖4
體模圖像
Figure4.
Phantom image
從 NIST 庫中可以查詢到特定物質在不同能量下的線性衰減系數。在體模圖像不同區域填充不同物質即可得到模體。本研究中仿真實驗采用的高低能峰值電壓分別是 140 kV 和 80 kV。假設探測器響應理想線性，140 kV 和 80 kV 能譜數據由開源能譜生成程序 SPECPRO 生成，如圖 5 所示。
圖5
80 kV 和 140 kV 能譜數據
Figure5.
The energy spectrum data of 80 kV and 140 kV
本實驗中采用弧形探測器，探測器數目為 525 個，旋轉角度間隔為 0.5°，一共 720 個投影角度，放射源到旋轉中心的距離是 550 mm，探測器中心到旋轉中心距離為 86.5 mm，編程工具為 Matlab，根據各個真實單能量圖像數據和 140 kV 以及 80 kV 能譜數據，生成 140 kV 和 80 kV 的混合能譜投影，如圖 6 所示。生成的高低能混合能譜投影數據就作為實測的高低能投影數據進行雙能 CT 圖像重建。
圖6
高低能投影數據仿真
Figure6.
The simulation of high and low polychromatic projections
2.2 生成高低能查找表
在本實驗中，基物質為 C 和 Fe，設定：BC-min = 0，BC-max = 30 cm，BFe-min = 0，BFe-min = 1 cm，dstep = 0.001 cm。
2.3 實驗結果
該實驗在聯想筆記本（CPU：i3；主頻：2.1 GH）上運行，使用直線加速算法（th = 50）和平面加速算法（th = 400）以及直線-平面混合加速算法對 140 kV 和 80 kV 的高低能投影數據進行雙能重建。每種加速算法運行 5 次，記錄 5 組數據，結果如表 1 所示。
表1
不同算法運行時間統計（單位：s）
Table1.
The run time statistics of different algorithms (Units: s)
基礎算法運行時間過長，只記錄了兩組數據。由表 1 可知，無論哪種加速方法，都比基礎算法快，加速效果明顯。為了保證準確性，平面加速算法設定的查找閾值比較大，遍歷計算量增加，因此運行時間相對較長。為了解決這一問題，混合使用直線加速和平面加速后，從表 1 可見，運行速度顯著提高。
下面通過實驗驗證算法的準確性。首先我們可以得到基物質分解系數圖像。如圖 7 所示，三種算法重建得到的結果幾乎相同。
圖7
C-Fe 基物質分解系數圖像
Figure7.
The decomposition coefficient image of C and Fe
根據式（2），我們可以重建各個能量下的單能量圖像，如圖 8 所示。
圖8
各個能量下的重建單能量圖像
Figure8.
The reconstructed mono-energy images of different energy
圖 8 顯示，三種算法重建的結果幾乎相同。用三種算法重建得到的單能圖像分別與真實單能圖像相減，得到圖 9。從圖 9 可見，三種加速算法得到的單能圖像與真實單能圖像的差值，除了邊界區域，其他區域的數值均接近 0，所以重建單能圖像和真實單能圖像非常接近。
圖9
重建單能量圖像與真實單能量圖像相減得到的圖像
Figure9.
The subtraction images of reconstructed mono-energy image and real mono-energy image
由式（3）和式（4）可以計算有效電子密度圖像和有效原子序數圖像，如圖 10 和圖 11 所示。
圖10
有效電子密度圖像
Figure10.
Effective electron density images
圖11
有效原子序數圖像
Figure11.
Effective atomic number image
從圖 10 可以看出三種算法重建的電子密度圖像幾乎是一樣的。圖 11 中重建的三種算法有效原子序數圖像也相差不多，但平面加速算法的中部區域卻有一些差異。這是因為平面加速算法的設定閾值 th 過小，造成了計算誤差。因為計算原子序數的公式有分母項，造成了誤差放大。加大平面加速算法的設定閾值 th，可以解決這一問題，但會增加計算時間。直線-平面混合加速算法也可以解決這個問題，而且計算時間相對于高閾值平面加速算法會降低。為了驗證這一猜想，我們用高閾值平面加速算法和直線-平面混合加速算法進行了重建，并得到相應原子序數圖，如圖 11 所示。從圖 11 可見，直線-平面混合加速算法和高閾值平面加速算法重建的效果與基礎算法、直線加速算法的效果是一樣。
選取圖像不同區域，計算不同區域內每個像素點對應的電子密度和原子序數數值，并求平均。理論上的電子密度值通過公式 （19）計算：
 |
ρ 為物質的密度，Zeff 是物質的有效原子序數，M 是摩爾質量，經過計算得到表 2。
表2
不同區域的電子密度數值計算
Table2.
The calculation of effective electron density in different region
其中相對誤差由重建值與理論值差的絕對值除以理論值，再乘以百分比得到。
單質的原子序數可以從元素周期表中得到，對于化合物的有效原子序數可以通過式（22）求得[12]：
 |
其中，Zj 是構成化合物的第 j 種元素的原子序數；m 一般在 3.0 到 3.8 的范圍取得，本文中取 3.5；aj 為構成化合物的各元素的電子數份額。經過計算得到表 3。根據表 2 和表 3 可以看出，無論是有效電子密度還是原子序數，三種算法重建值的相對誤差都在一個可接受的范圍內，進一步驗證了加速算法的有效性。
表3
不同區域的原子序數數值計算
Table3.
The calculation of effective atomic number in different region
3 結論
基于投影匹配的雙能 CT 投影分解算法避免了求解基物質分解系數投影值復雜的迭代優化過程，簡化了系統的標定過程。分解精度取決于查找表的設定步進，可以取得穩定的正確解。其缺點是求解過程需要遍歷龐大的查找表數據，計算量大，運行時間長，不利于臨床應用。本文提出的基于投影匹配的雙能 CT 成像投影分解加速算法，大大提高了算法的運行效率，對于給定實測的高低能投影數據，可以快速計算出對應的基物質分解系數的投影值，進而求出基物質分解系數圖像，重建單能量圖像、有效電子密度圖像和有效原子序數圖像。仿真實驗結果驗證了算法的有效性。同時該算法也易于并行計算，從而提高了該算法在臨床上的實用性。
引言
計算機斷層成像（computed tomography，CT）所用的 X 射線能譜包含不同能量的 X 射線，就如白光包含不同顏色的光。傳統 CT 成像技術把 X 射線能譜簡化為等價的單一能量的 X 射線，利用不同物質對此單一能量 X 射線的衰減信息不同得到掃描組織的形態學結構圖像。雙能 CT 成像技術不對 X 射線能譜進行簡化，能同時得到掃描組織的形態學結構圖像和化學組織成分圖像，進行化學組織成分分析，這是傳統 CT 成像技術無法企及的優點[1]。雙能 CT 成像技術的關鍵是雙能 CT 重建算法。重建算法的模型主要有雙效應分解模型和基物質分解模型。本文主要討論的是基物質分解模型。
根據重建方法的不同，雙能 CT 重建算法可以分為三大類[2]：① 投影分解重建算法，該類方法主要是通過高、低能投影數據得到分解系數投影值，然后由投影值反投影求出分解系數。② 后處理重建算法，該類方法首先使用高、低能投影數據重建出高、低能 CT 圖像，然后使用高、低能 CT 圖像求解分解系數，這種方法精確度不高。③ 迭代重建算法，通過雙能迭代，得到分解系數圖像。無論哪種方法，最終目的都是求出分解系數。李保磊、張耀軍在 2011 年提出基于投影匹配的雙能 CT 投影分解重建算法[1]，該算法屬于第一類。依據 CT 系統能譜和基物質線性衰減曲線，通過投影積分方程組建立 X 射線的高低能譜查找表，利用沿著某條 X 射線經過的路徑下實測得到的高低能投影數據，在查找表中尋找其最佳匹配點，該最佳匹配點就是對應的基物質分解系數的投影值，然后通過傳統的 CT 重建算法如濾波反投影重建算法（filter back projection，FBP）[3]重建求出對應的分解系數圖像。但實際應用中，由于查找表數據龐大，在尋找最佳匹配點的過程中需要遍歷整張查找表，計算過程十分耗時，不利于臨床應用。該算法的計算量大部分在于尋找最佳匹配點。針對這一點，本文在該算法基礎上提出一種加速算法。
1 方法
1.1 雙能 CT 重建算法的原理[4 -8 ]
雙能 CT 重建算法模型主要是雙效應分解模型和基物質分解模型：
|
|
u（E）代表物質在能量 E 下的線性衰減系數，式（1）對應的是雙效應分解，在一定的射線能量范圍內，物質的線性衰減可由光電效應和康普頓散射兩種作用共同組成。其中 fp（E）、fkn（E）只與能量 E 有關而與材質無關，ac、ap 分別對應的是光電效應和康普頓效應的分解系數。式（2）對應基物質分解模型，表示任何一種物質的線性衰減系數都可由兩種基物質的線性衰減系數線性表示。其中，u1（E）、u2（E）是兩種基物質在能量 E 下對應的線性衰減值，該數據可以在美國國家標準與技術研究院（National Institute of Standards and Technology，NIST）庫中查詢得到，b1、b2 是兩種基物質對應的分解系數。在本文中，我們主要是討論基物質分解模型。由式（2）可知，在某兩個單能量如 E1、E2 下，就可以得到兩個方程，其中 u（E）、u1（E）、u2（E）是已知數，而 b1、b2 是未知數，通過兩個方程組求解即可得到分解系數 b1、b2。在實際中，球管產生的 X 射線并不是單能譜射線，而是連續的混合能譜射線，因此通過式（2）無法求出基物質的分解系數。
根據基物質分解模型，在沿著某一條 X 射線的路徑下有：
|
B1、B2 分別表示兩種基物質的分解系數投影值，根據混合能譜 X 射線條件下的朗伯比爾定律可得：
|
其中，
代表沿著某條 X 射線路徑下高低能投影值，
為 X 射線的高低能譜數據。由式（3）可知，只要求出基物質分解系數投影值 B1、B2，再通過 FBP 即可得到基物質分解系數圖像。根據式（4）和式（5），我們可以求出重建圖像中每個像素點對應的有效原子序數 Zeff 和有效電子密度 ρe：
|
|
其中 ρe1 和 ρe2 是基物質的有效電子密度，Z1 和 Z2 是基物質的有效原子序數。n 的值為 3.5。綜上所述，如何計算分解系數投影值就成為雙能 CT 重建算法的難點。基于投影匹配的雙能 CT 投影分解重建算法就是為了計算分解系數投影值而提出的。
1.2 基于投影匹配的雙能 CT 投影分解重建算法[1 ]
基于投影匹配的雙能 CT 投影分解重建算法的主要思想是：根據 CT 系統能譜信息生成合理意義范圍內的一系列 B1 和 B2 數值，再根據式（3）來生成對應的高低能投影數值 PH 和 PL，建立查找表數據，然后根據給定的實測高低能投影數據
、
，匹配查找表，尋找最佳匹配點，即分解系數投影值
和
。主要包含三個步驟：
（1）設定 B1 和 B2 的范圍，
，
，根據設定的步進 d 來產生 B1（i）和 B2（j）序列。根據系統能譜數據和式（3）計算對應的高低能投影值 PH 和 PL。我們用 i 來代表查找表的行編號，j 代表查找表的列編號，則有：
|
（2）根據實測高低能投影數據
和
，利用最小二乘法，在查找表中進行匹配比較，使得式（6）的值最小化。
|
（3）根據最佳匹配點的坐標（i0，j0）來計算分解系數投影值：
|
只要得到分解系數投影值，通過 FBP 就能求出基物質分解系數圖像，進而重建單能量圖像，擬合特定區域物質的線性衰減曲線以及重建有效電子密度圖像和有效原子序數圖像。該算法避開了復雜的數學公式，實現起來很簡單。但為了提高計算精度，需要把查找表的步進設定成非常小的數值，同時 B1、B2 范圍比較大，因此查找表非常龐大，查找最佳匹配點也就非常耗時。為了解決這個問題，本文提出了一種加速算法。
1.3 基于投影匹配的雙能 CT 成像投影分解加速算法
本文提出一種查找最佳匹配點的加速算法，其基本思想是：對于一對實測的高低能投影數據，通過某種方法確定其在查找表中對應的最佳匹配點所在的大致區域，然后對該區域進行遍歷、匹配和比較，找到最佳匹配點。該方法可以總結為以下兩種算法。
1.3.1 基于直線擬合查找表數據的加速投影分解算法
以 C 和 Fe 為基物質的查找表為例，把查找表中某一行數據取出來畫圖，如圖 1 所示。
從圖中可以看出，高低能投影值關于變量 j 呈現近似線性關系。根據這個規律，可以建立直線方程來擬合該行查找表數據[9]。
|
|
其中 PL（i，j）和 PH（i，j）表示查找表中坐標為（i，j）的數據點的低能和高能投影數值。aL、bL、aH、bH 分別是直線方程系數。式（6）用 f（j）來表示，把式（8）和式（9）代入得到：
|
圖1
某一行 C-Fe 查找表數據畫圖顯示
Figure1.
The drawing display of a line of C-Fe lookup table data
要尋找最佳匹配點，必須最小化 f（j）的值。f（j）是關于變量 j 的一元二次函數，同時也是開口向上的凸函數，其最小值要么在極小值點取得，要么在邊界點取得。對變量 j 進行求導，并令導數等于 0，可以得到極小值點：
|
若在極小值點取得最小值，那么最佳匹配點就在極值點附近。因此設定一個閾值 th，如圖 2 所示，遍歷區間
中每一個點，找到最佳匹配點。
圖2
某一行查找表數據一維表示
Figure2.
One-dimensional representation of a row of lookup table data
若 j0 < jmin 或者 j0 > jmax，函數 f（j）的最小值只能在邊界點 jmin 和 jmax 取得。
對于一組給定的實測高低能投影數據，上述方法只能獲取某一行數據的最佳匹配點，稱為行最佳匹配點。使用同樣的方法可以求出其他行的行最佳匹配點，最后比較所有的行最佳匹配點，找到全局最佳匹配點。該方法需要遍歷查找表的每一行，但避免了遍歷行內所有數據點，大大減少了計算量。
1.3.2 基于平面擬合查找表數據的加速投影分解算法
基于平面擬合查找表數據的加速投影分解算法用平面方程擬合查找表數據。通過建立目標函數方程來求解該高低能投影數據對應于查找表的最佳匹配點，求出對應的分解系數的投影值。圖 3 是高低能查找表三維空間顯示結果。由圖可見，高低能查找表在三維空間近似一個平面，可以通過建立平面方程[10]擬合高低能查找表。
|
|
式（12）和式（13）分別表示擬合低高能投影數據的平面方程，AL、BL、CL、AH、BH、CH 都是平面方程的系數，把式（12）和式（13）代入式（6）可得：
|
圖3
C-Fe 查找表數據在三維空間顯示
Figure3.
C-Fe lookup table data are displayed in three-dimensional space
由式（14）可知，f（i，j）是關于變量 i 和 j 的開口向上的凸函數，在區間
&&
內存在極小值點。對變量 i 和變量 j 求偏導：
|
為了簡化式子，用下面符號代表不同的方程系數：
|
則可以得到下面兩個方程組：
|
兩個方程組聯立求解最終得到 i0 和 j0：
|
|
若極小值點（i0，j0）是最小值點，最佳匹配點一定在極小值點（i0，j0）附近。因此可以設定一個閾值 th，在區間
&&
內每個數據點進行匹配計算，找到最佳匹配點。
若極小值點（i0，j0）不在設定的區間范圍內，那么最小值在查找表的邊界上取得。對應邊界區間分為四個：
|
|
|
|
在邊界區間中每個數據點進行匹配計算找到最佳匹配點，從而求出實際投影數據對應的分解系數投影值。基于平面擬合查找表數據加速投影分解算法不需要遍歷查找表的每一行，只要給出實測的高低能投影數據，就能定位對應的最佳匹配點的大致位置。同時，該方法也需要設定一個閾值，在一個特定區域遍歷和匹配計算，比較每一個數據點，才能找到最佳匹配點。
總結來說，兩種方法的思路都是先進行粗糙搜索，定位最佳匹配點的大致區域，然后進行精確搜索，找到最佳匹配點，計算得到分解系數投影值。兩種方法各有優缺點。對于列的數目遠遠多于行的數目的查找表，使用基于直線擬合查找表數據加速投影分解算法可以很快得到解。對于行列數相差不多的查找表，基于平面擬合查找表數據加速投影分解算法比較有優勢。但是為了保證準確性，基于平面擬合查找表數據加速投影分解算法的閾值 th 一般比較大，遍歷區域也就比較大，這會增加計算量。可以對這兩種方法進行混合使用，首先通過平面加速算法確定遍歷區域，然后在此遍歷區域上使用直線加速算法確定最佳匹配點，就可以起到很好的加速效果。
采用上述加速算法計算得到基物質分解系數投影值以后，通過傳統的 CT 重建算法如 FB 算法重建出對應的基物質分解系數圖像，然后就可以重建單能量圖像、有效電子密度圖像和有效原子序數圖像。
2 仿真實驗
為了敘述簡單，我們把基于投影匹配的雙能 CT 投影分解重建算法簡稱為基礎算法（basic algorithm），基于直線擬合查找表數據的加速投影分解算法簡稱為直線加速算法（linear acceleration algorithm），基于平面擬合查找表數據的加速投影分解算法簡稱為平面加速算法（plane acceleration algorithm），直線加速算法和平面加速算法混合使用則簡稱為直線-平面混合加速算法（the mixture algorithm of line-plane acceleration）。
2.1 產生仿真數據
根據已有的論文，設計體模圖像[11]，圖像的大小為 256 × 256 像素，空間分辨率為 1 mm，如圖 4 所示。
圖4
體模圖像
Figure4.
Phantom image
從 NIST 庫中可以查詢到特定物質在不同能量下的線性衰減系數。在體模圖像不同區域填充不同物質即可得到模體。本研究中仿真實驗采用的高低能峰值電壓分別是 140 kV 和 80 kV。假設探測器響應理想線性，140 kV 和 80 kV 能譜數據由開源能譜生成程序 SPECPRO 生成，如圖 5 所示。
圖5
80 kV 和 140 kV 能譜數據
Figure5.
The energy spectrum data of 80 kV and 140 kV
本實驗中采用弧形探測器，探測器數目為 525 個，旋轉角度間隔為 0.5°，一共 720 個投影角度，放射源到旋轉中心的距離是 550 mm，探測器中心到旋轉中心距離為 86.5 mm，編程工具為 Matlab，根據各個真實單能量圖像數據和 140 kV 以及 80 kV 能譜數據，生成 140 kV 和 80 kV 的混合能譜投影，如圖 6 所示。生成的高低能混合能譜投影數據就作為實測的高低能投影數據進行雙能 CT 圖像重建。
圖6
高低能投影數據仿真
Figure6.
The simulation of high and low polychromatic projections
2.2 生成高低能查找表
在本實驗中，基物質為 C 和 Fe，設定：BC-min = 0，BC-max = 30 cm，BFe-min = 0，BFe-min = 1 cm，dstep = 0.001 cm。
2.3 實驗結果
該實驗在聯想筆記本（CPU：i3；主頻：2.1 GH）上運行，使用直線加速算法（th = 50）和平面加速算法（th = 400）以及直線-平面混合加速算法對 140 kV 和 80 kV 的高低能投影數據進行雙能重建。每種加速算法運行 5 次，記錄 5 組數據，結果如表 1 所示。
表1
不同算法運行時間統計（單位：s）
Table1.
The run time statistics of different algorithms (Units: s)
基礎算法運行時間過長，只記錄了兩組數據。由表 1 可知，無論哪種加速方法，都比基礎算法快，加速效果明顯。為了保證準確性，平面加速算法設定的查找閾值比較大，遍歷計算量增加，因此運行時間相對較長。為了解決這一問題，混合使用直線加速和平面加速后，從表 1 可見，運行速度顯著提高。
下面通過實驗驗證算法的準確性。首先我們可以得到基物質分解系數圖像。如圖 7 所示，三種算法重建得到的結果幾乎相同。
圖7
C-Fe 基物質分解系數圖像
Figure7.
The decomposition coefficient image of C and Fe
根據式（2），我們可以重建各個能量下的單能量圖像，如圖 8 所示。
圖8
各個能量下的重建單能量圖像
Figure8.
The reconstructed mono-energy images of different energy
圖 8 顯示，三種算法重建的結果幾乎相同。用三種算法重建得到的單能圖像分別與真實單能圖像相減，得到圖 9。從圖 9 可見，三種加速算法得到的單能圖像與真實單能圖像的差值，除了邊界區域，其他區域的數值均接近 0，所以重建單能圖像和真實單能圖像非常接近。
圖9
重建單能量圖像與真實單能量圖像相減得到的圖像
Figure9.
The subtraction images of reconstructed mono-energy image and real mono-energy image
由式（3）和式（4）可以計算有效電子密度圖像和有效原子序數圖像，如圖 10 和圖 11 所示。
圖10
有效電子密度圖像
Figure10.
Effective electron density images
圖11
有效原子序數圖像
Figure11.
Effective atomic number image
從圖 10 可以看出三種算法重建的電子密度圖像幾乎是一樣的。圖 11 中重建的三種算法有效原子序數圖像也相差不多，但平面加速算法的中部區域卻有一些差異。這是因為平面加速算法的設定閾值 th 過小，造成了計算誤差。因為計算原子序數的公式有分母項，造成了誤差放大。加大平面加速算法的設定閾值 th，可以解決這一問題，但會增加計算時間。直線-平面混合加速算法也可以解決這個問題，而且計算時間相對于高閾值平面加速算法會降低。為了驗證這一猜想，我們用高閾值平面加速算法和直線-平面混合加速算法進行了重建，并得到相應原子序數圖，如圖 11 所示。從圖 11 可見，直線-平面混合加速算法和高閾值平面加速算法重建的效果與基礎算法、直線加速算法的效果是一樣。
選取圖像不同區域，計算不同區域內每個像素點對應的電子密度和原子序數數值，并求平均。理論上的電子密度值通過公式 （19）計算：
|
ρ 為物質的密度，Zeff 是物質的有效原子序數，M 是摩爾質量，經過計算得到表 2。
表2
不同區域的電子密度數值計算
Table2.
The calculation of effective electron density in different region
其中相對誤差由重建值與理論值差的絕對值除以理論值，再乘以百分比得到。
單質的原子序數可以從元素周期表中得到，對于化合物的有效原子序數可以通過式（22）求得[12]：
|
其中，Zj 是構成化合物的第 j 種元素的原子序數；m 一般在 3.0 到 3.8 的范圍取得，本文中取 3.5；aj 為構成化合物的各元素的電子數份額。經過計算得到表 3。根據表 2 和表 3 可以看出，無論是有效電子密度還是原子序數，三種算法重建值的相對誤差都在一個可接受的范圍內，進一步驗證了加速算法的有效性。
表3
不同區域的原子序數數值計算
Table3.
The calculation of effective atomic number in different region
3 結論
基于投影匹配的雙能 CT 投影分解算法避免了求解基物質分解系數投影值復雜的迭代優化過程，簡化了系統的標定過程。分解精度取決于查找表的設定步進，可以取得穩定的正確解。其缺點是求解過程需要遍歷龐大的查找表數據，計算量大，運行時間長，不利于臨床應用。本文提出的基于投影匹配的雙能 CT 成像投影分解加速算法，大大提高了算法的運行效率，對于給定實測的高低能投影數據，可以快速計算出對應的基物質分解系數的投影值，進而求出基物質分解系數圖像，重建單能量圖像、有效電子密度圖像和有效原子序數圖像。仿真實驗結果驗證了算法的有效性。同時該算法也易于并行計算，從而提高了該算法在臨床上的實用性。
