將基于經典Demons算法與加速Demons算法的彈性配準模型應用到放療患者不同分次間的錐形束CT(CBCT)圖像中,為進一步分析腫瘤及危及器官的變化提供軟件支持。通過Matlab軟件編寫三維彈性配準程序,并用此程序對宮頸癌放療患者不同分次時的兩組錐形束CT圖像進行仿真實驗驗證。結果表明經典Demons算法配準前后對比,最小均方誤(MSE)減少59.7%,相關系數(CC)提高了11.0%;加速Demons算法配準前后對比,MSE減少40.1%,CC提高7.2%。實驗驗證上述兩種基于Demons模型的彈性配準在CBCT圖像配準中取得較好的配準效果,但對細微差別處仍顯得精度不夠,且整個配準時間較長,如要應用到在線自適應放療中仍需進一步提高形變的精度并減少配準時間。
引用本文： 龐皓文, 孫小揚. 基于Demons算法的彈性配準在放療錐形束CT圖像配準中的應用. 生物醫學工程學雜志, 2014, 31(1): 103-106. doi: 10.7507/1001-5515.20140020 復制
引言
近年來圖像引導放療(image guided radiation therapy,IGRT)技術已經廣泛地運用于放療中，通過采集患者分次間錐形束CT(cone beam CT,CBCT)圖像，提高了靶區劑量精度，保護危及器官，為實時觀察靶區及危及器官解剖結構的變化提供了幫助[1]。為了進一步了解這些變化，可以通過彈性配準的方法將分次間CBCT圖像進行配準，讓臨床醫生及時了解腫瘤區域以及周圍正常器官變化，進而計算腫瘤及正常器官的實際累加劑量，及時修改治療計劃，提高放療精確度，為個體化自適應放療提供支持[2-5]。本文基于Matlab軟件編程，對放療患者CBCT圖像實施基于Demons算法的彈性配準，為進一步分析腫瘤及危及器官的變化提供軟件支持，并成功應用于臨床病例。
1 方法
1.1 Demons算法原理
Demons彈性配準算法是由Thirion[6]提出，其基本思想與Maxwell解決Gibbs效應的方法類似，通過判斷待配準圖像上各個像素點的驅動力，然后移動像素點來實現彈性配準。其驅動力來源于光流場理論方程，光流場理論的前提是假設圖像運動過程中灰度值保持不變，即
| $I\left( x,\text{ }y,\text{ }z,t \right)=\text{ }I\left( x+\Delta x,\text{ }y+\Delta y,z+\Delta z,\text{ }t+\Delta t \right),$ |
將式(1) 進行泰勒展開
| $\frac{\partial I}{\partial x}\frac{dx}{dt}+\frac{\partial I}{\partial y}\frac{dy}{dt}+\frac{\partial I}{\partial z}\frac{dz}{dt}=\frac{\partial I}{\partial t}$ |
設v=(x，y，z)，M(v)是待配準圖像灰度值，S(v)是參考圖像灰度值，式(2) 進一步簡化為
| ${{f}_{s}}\centerdot \nabla S\left( v \right)=M\left( v \right)-S\left( v \right),$ |
其中為驅動力，，可求得
| ${{f}_{s}}=\frac{M\left( v \right)-S\left( v \right)\nabla S\left( v \right)}{{{\left| \nabla S\left( v \right) \right|}^{2}}}$ |
但當?S→0時，式(4) 的定義變得不穩定，會導致較大的驅動力fs，因此在分母增加一個分量，得到
| ${{f}_{s}}=\frac{M\left( v \right)-S\left( v \right)\nabla S\left( v \right)}{{{\left| \nabla S\left( v \right) \right|}^{2}}+{{\left( M\left( v \right)-S\left( v \right) \right)}^{2}}}$ |
1.2 Demons算法的優化
1.2.1 能量函數構造
圖像彈性形變配準是尋找位移場u(v)，使配準圖像形變后M(v-u)與參考圖像S(v)在某種定義上相似。最常用的能量函數(相似性測度函數)為
| $E=\frac{1}{2}\int_{\Omega }{{{\left| M\left( v-u \right)-S\left( v \right) \right|}^{2}}dv}$ |
最小化式(6) ，利用下降法求解Euler方程得到
| $\frac{\partial u\left( v,t \right)}{\partial t}=\left[ M\left( v-u \right)-S\left( v \right) \right]{{\left. \nabla M \right|}_{v-u}}$ |
可以看到式(7) 與式(5) 相似，因此假設式(5) 即為最小化式(6) 得到的位移場，可用迭代法求出
| ${{u}_{n=+1}}\left( v \right)={{u}_{n}}\left( v \right)+\frac{M\left( v \right)-S\left( v \right)\nabla S\left( v \right)}{{{\left| \nabla S\left( v \right) \right|}^{2}}+{{\left( M\left( v \right)-S\left( v \right) \right)}^{2}}}$ |
為了在全局范圍內使該變換連續，從而保持圖像的拓撲結構，在每一次迭代后，使用高斯濾波來平滑所得到的驅動力場，如下式所示
| $\begin{align} & {{u}_{n=+1}}\left( v \right)={{G}_{\partial }}\otimes \left( {{u}_{n}}\left( v \right) \right.+ \\ & \frac{M\left( v \right)-S\left( v \right)\nabla S\left( v \right)}{{{\left| \nabla S\left( v \right) \right|}^{2}}+{{\left( M\left( v \right)-S\left( v \right) \right)}^{2}}}\left. {} \right), \\ \end{align}$ |
式中Gσ是窗寬為σ高斯低通濾波函數，n為迭代次數。
1.2.2 加速Demons算法
Rogelj等[7]基于牛頓第三定律，在圖像形變過程中增加了一個新的驅動力，將驅動浮動圖像到參考圖像的力作為正驅動力即fs，驅動參考圖像到浮動圖像的力作為負驅動力fm，即
| ${{f}_{m}}=-\frac{S\left( v \right)-M\left( v \right)\nabla M\left( v \right)}{{{\left| \nabla M\left( v \right) \right|}^{2}}+{{\left( S\left( v \right)-M\left( v \right) \right)}^{2}}},$ |
式中負號代表了負驅動力。聯合正驅動力，則總驅動力為
| $\begin{align} & f={{f}_{m}}+{{f}_{s}}=\left( M\left( v \right)-S\left( v \right) \right)\left[ \frac{\nabla S\left( v \right)}{{{\left| \nabla S\left( v \right) \right|}^{2}}+{{\left( M\left( v \right)-S\left( v \right) \right)}^{2}}} \right. \\ & +\left. \frac{\nabla M\left( v \right)}{{{\left| \nabla M\left( v \right) \right|}^{2}}+{{\left( S\left( v \right)-M\left( v \right) \right)}^{2}}} \right], \\ \end{align}$ |
Cachier等[8]為了控制每次迭代的步長，引入了歸一化因子α，式(11)改寫為
| $\begin{align} & f={{f}_{m}}+{{f}_{s}}=\left( M\left( v \right)-S\left( v \right) \right)\left[ \frac{\nabla S\left( v \right)}{{{\left| \nabla S\left( v \right) \right|}^{2}}+{{\alpha }^{2}}{{\left( M\left( v \right)-S\left( v \right) \right)}^{2}}} \right. \\ & +\left. \frac{\nabla M\left( v \right)}{{{\left| \nabla M\left( v \right) \right|}^{2}}+{{\alpha }^{2}}{{\left( S\left( v \right)-M\left( v \right) \right)}^{2}}} \right], \\ \end{align}$ |
因此，位移場可由下式進行迭代
| $\begin{align} & {{u}_{n=+1}}\left( v \right)={{G}_{\partial }}\otimes \left( {{u}_{n}}\left( v \right)+M\left( v \right)-S\left( v \right) \right) \\ & \left[ \frac{\nabla S\left( v \right)}{{{\left| \nabla S\left( v \right) \right|}^{2}}+{{\alpha }^{2}}{{\left( M\left( v \right)-S\left( v \right) \right)}^{2}}} \right.+ \\ & \left. \frac{\nabla M\left( v \right)}{{{\left| \nabla M\left( v \right) \right|}^{2}}+{{\alpha }^{2}}{{\left( S\left( v \right)-M\left( v \right) \right)}^{2}}} \right], \\ \end{align}$ |
本研究中σ=0.3，α取0.4。
1.2.3 基于Demons模型的圖像彈性配準
以M(v)與S(v)分別表示帶配準圖像和參考圖像灰度值，u(v)為位移場。步驟：① 首先初始化位移場，迭代次數n=0，u0(v)=0 ；② 設n=n+1，根據公式(9) 或(13)更新位移場un+1，更新浮動圖像為M(v-un+1(v))； ③ 判斷迭代次數大于設定值，則轉步驟④ ，否則執行步驟② ； ④ 迭代終止，M(v-un+1(v))為最終的配準圖像。
2 結果
由作者通過Matlab 7.6圖像處理工具，根據上述Demons算法原理編寫三維彈性配準程序，并用此程序對宮頸癌放療患者不同分次時的兩組CBCT圖像進行仿真實驗驗證。實驗平臺為Intel Core Duo CPU T3000(酷睿2雙核，2.0 GHz)，1G內存。患者CBCT圖像斷層大小均為512×512，層距為2.5 mm，共采集63層。實驗結果如圖 1所示，實驗參數如表 1所示。經典Demons算法配準前后對比，MSE減少59.7%，CC提高了11.0%；加速Demons算法配準前后對比，MSE減少40.1%，CC提高7.2%。經典Demons算法配準耗時975 s，加速Demons算法配準耗時513 s。
圖1
基于Demons算法的宮頸癌放療患者CBCT圖像彈性配準冠狀面圖
(a)參考圖像;(b)待配準圖像;(c)配準前差分圖;(d)Demons算法配準后差分圖;(e)加速Demons算法配準后差分圖
Figure1. Elastic registration coronal CBCT image of patients with cervical cancer based on the Demons algorithm(a) reference image; (b) pending registration image; (c) difference image before registration; (d) difference image after registration for Demons algorithm; (e) difference image after registration for accelerated Demons algorithm
表1
算法評價參數
Table1.
Algorithm evaluation parameters
3 結論
整個放療期間，由于腫瘤的退縮以及周圍正常器官的運動和充盈狀態的不一，這或許會使其實際吸收劑量發生較大改變，因此及時了解這些變化是臨床醫生十分關心的問題。彈性配準為實際觀察這些變化提供了軟件支持。本文所述兩種基于Demons算法的彈性配準對宮頸癌放療患者CBCT圖像進行實驗，從配準前后差分圖的比較及表 1兩種算法的評價參數可見，兩種配準方法均取得較好的配準效果,但耗時均較長。其中Demons算法的評價參數均好于加速Demons算法，但其配準時間是加速Demons算法的1.9倍。對于細微差別處，兩種模型仍顯得精度不夠。在放療過程中，要獲得器官真實的累加劑量，要求形變的精確性非常高。怎樣將二者結合起來，更好地模擬器官運動，為患者提供更快更精確的個體化放療是未來研究的熱點。
引言
近年來圖像引導放療(image guided radiation therapy,IGRT)技術已經廣泛地運用于放療中，通過采集患者分次間錐形束CT(cone beam CT,CBCT)圖像，提高了靶區劑量精度，保護危及器官，為實時觀察靶區及危及器官解剖結構的變化提供了幫助[1]。為了進一步了解這些變化，可以通過彈性配準的方法將分次間CBCT圖像進行配準，讓臨床醫生及時了解腫瘤區域以及周圍正常器官變化，進而計算腫瘤及正常器官的實際累加劑量，及時修改治療計劃，提高放療精確度，為個體化自適應放療提供支持[2-5]。本文基于Matlab軟件編程，對放療患者CBCT圖像實施基于Demons算法的彈性配準，為進一步分析腫瘤及危及器官的變化提供軟件支持，并成功應用于臨床病例。
1 方法
1.1 Demons算法原理
Demons彈性配準算法是由Thirion[6]提出，其基本思想與Maxwell解決Gibbs效應的方法類似，通過判斷待配準圖像上各個像素點的驅動力，然后移動像素點來實現彈性配準。其驅動力來源于光流場理論方程，光流場理論的前提是假設圖像運動過程中灰度值保持不變，即
| $I\left( x,\text{ }y,\text{ }z,t \right)=\text{ }I\left( x+\Delta x,\text{ }y+\Delta y,z+\Delta z,\text{ }t+\Delta t \right),$ |
將式(1) 進行泰勒展開
| $\frac{\partial I}{\partial x}\frac{dx}{dt}+\frac{\partial I}{\partial y}\frac{dy}{dt}+\frac{\partial I}{\partial z}\frac{dz}{dt}=\frac{\partial I}{\partial t}$ |
設v=(x，y，z)，M(v)是待配準圖像灰度值，S(v)是參考圖像灰度值，式(2) 進一步簡化為
| ${{f}_{s}}\centerdot \nabla S\left( v \right)=M\left( v \right)-S\left( v \right),$ |
其中為驅動力，，可求得
| ${{f}_{s}}=\frac{M\left( v \right)-S\left( v \right)\nabla S\left( v \right)}{{{\left| \nabla S\left( v \right) \right|}^{2}}}$ |
但當?S→0時，式(4) 的定義變得不穩定，會導致較大的驅動力fs，因此在分母增加一個分量，得到
| ${{f}_{s}}=\frac{M\left( v \right)-S\left( v \right)\nabla S\left( v \right)}{{{\left| \nabla S\left( v \right) \right|}^{2}}+{{\left( M\left( v \right)-S\left( v \right) \right)}^{2}}}$ |
1.2 Demons算法的優化
1.2.1 能量函數構造
圖像彈性形變配準是尋找位移場u(v)，使配準圖像形變后M(v-u)與參考圖像S(v)在某種定義上相似。最常用的能量函數(相似性測度函數)為
| $E=\frac{1}{2}\int_{\Omega }{{{\left| M\left( v-u \right)-S\left( v \right) \right|}^{2}}dv}$ |
最小化式(6) ，利用下降法求解Euler方程得到
| $\frac{\partial u\left( v,t \right)}{\partial t}=\left[ M\left( v-u \right)-S\left( v \right) \right]{{\left. \nabla M \right|}_{v-u}}$ |
可以看到式(7) 與式(5) 相似，因此假設式(5) 即為最小化式(6) 得到的位移場，可用迭代法求出
| ${{u}_{n=+1}}\left( v \right)={{u}_{n}}\left( v \right)+\frac{M\left( v \right)-S\left( v \right)\nabla S\left( v \right)}{{{\left| \nabla S\left( v \right) \right|}^{2}}+{{\left( M\left( v \right)-S\left( v \right) \right)}^{2}}}$ |
為了在全局范圍內使該變換連續，從而保持圖像的拓撲結構，在每一次迭代后，使用高斯濾波來平滑所得到的驅動力場，如下式所示
| $\begin{align} & {{u}_{n=+1}}\left( v \right)={{G}_{\partial }}\otimes \left( {{u}_{n}}\left( v \right) \right.+ \\ & \frac{M\left( v \right)-S\left( v \right)\nabla S\left( v \right)}{{{\left| \nabla S\left( v \right) \right|}^{2}}+{{\left( M\left( v \right)-S\left( v \right) \right)}^{2}}}\left. {} \right), \\ \end{align}$ |
式中Gσ是窗寬為σ高斯低通濾波函數，n為迭代次數。
1.2.2 加速Demons算法
Rogelj等[7]基于牛頓第三定律，在圖像形變過程中增加了一個新的驅動力，將驅動浮動圖像到參考圖像的力作為正驅動力即fs，驅動參考圖像到浮動圖像的力作為負驅動力fm，即
| ${{f}_{m}}=-\frac{S\left( v \right)-M\left( v \right)\nabla M\left( v \right)}{{{\left| \nabla M\left( v \right) \right|}^{2}}+{{\left( S\left( v \right)-M\left( v \right) \right)}^{2}}},$ |
式中負號代表了負驅動力。聯合正驅動力，則總驅動力為
| $\begin{align} & f={{f}_{m}}+{{f}_{s}}=\left( M\left( v \right)-S\left( v \right) \right)\left[ \frac{\nabla S\left( v \right)}{{{\left| \nabla S\left( v \right) \right|}^{2}}+{{\left( M\left( v \right)-S\left( v \right) \right)}^{2}}} \right. \\ & +\left. \frac{\nabla M\left( v \right)}{{{\left| \nabla M\left( v \right) \right|}^{2}}+{{\left( S\left( v \right)-M\left( v \right) \right)}^{2}}} \right], \\ \end{align}$ |
Cachier等[8]為了控制每次迭代的步長，引入了歸一化因子α，式(11)改寫為
| $\begin{align} & f={{f}_{m}}+{{f}_{s}}=\left( M\left( v \right)-S\left( v \right) \right)\left[ \frac{\nabla S\left( v \right)}{{{\left| \nabla S\left( v \right) \right|}^{2}}+{{\alpha }^{2}}{{\left( M\left( v \right)-S\left( v \right) \right)}^{2}}} \right. \\ & +\left. \frac{\nabla M\left( v \right)}{{{\left| \nabla M\left( v \right) \right|}^{2}}+{{\alpha }^{2}}{{\left( S\left( v \right)-M\left( v \right) \right)}^{2}}} \right], \\ \end{align}$ |
因此，位移場可由下式進行迭代
| $\begin{align} & {{u}_{n=+1}}\left( v \right)={{G}_{\partial }}\otimes \left( {{u}_{n}}\left( v \right)+M\left( v \right)-S\left( v \right) \right) \\ & \left[ \frac{\nabla S\left( v \right)}{{{\left| \nabla S\left( v \right) \right|}^{2}}+{{\alpha }^{2}}{{\left( M\left( v \right)-S\left( v \right) \right)}^{2}}} \right.+ \\ & \left. \frac{\nabla M\left( v \right)}{{{\left| \nabla M\left( v \right) \right|}^{2}}+{{\alpha }^{2}}{{\left( S\left( v \right)-M\left( v \right) \right)}^{2}}} \right], \\ \end{align}$ |
本研究中σ=0.3，α取0.4。
1.2.3 基于Demons模型的圖像彈性配準
以M(v)與S(v)分別表示帶配準圖像和參考圖像灰度值，u(v)為位移場。步驟：① 首先初始化位移場，迭代次數n=0，u0(v)=0 ；② 設n=n+1，根據公式(9) 或(13)更新位移場un+1，更新浮動圖像為M(v-un+1(v))； ③ 判斷迭代次數大于設定值，則轉步驟④ ，否則執行步驟② ； ④ 迭代終止，M(v-un+1(v))為最終的配準圖像。
2 結果
由作者通過Matlab 7.6圖像處理工具，根據上述Demons算法原理編寫三維彈性配準程序，并用此程序對宮頸癌放療患者不同分次時的兩組CBCT圖像進行仿真實驗驗證。實驗平臺為Intel Core Duo CPU T3000(酷睿2雙核，2.0 GHz)，1G內存。患者CBCT圖像斷層大小均為512×512，層距為2.5 mm，共采集63層。實驗結果如圖 1所示，實驗參數如表 1所示。經典Demons算法配準前后對比，MSE減少59.7%，CC提高了11.0%；加速Demons算法配準前后對比，MSE減少40.1%，CC提高7.2%。經典Demons算法配準耗時975 s，加速Demons算法配準耗時513 s。
圖1
基于Demons算法的宮頸癌放療患者CBCT圖像彈性配準冠狀面圖
(a)參考圖像;(b)待配準圖像;(c)配準前差分圖;(d)Demons算法配準后差分圖;(e)加速Demons算法配準后差分圖
Figure1. Elastic registration coronal CBCT image of patients with cervical cancer based on the Demons algorithm(a) reference image; (b) pending registration image; (c) difference image before registration; (d) difference image after registration for Demons algorithm; (e) difference image after registration for accelerated Demons algorithm
表1
算法評價參數
Table1.
Algorithm evaluation parameters
3 結論
整個放療期間，由于腫瘤的退縮以及周圍正常器官的運動和充盈狀態的不一，這或許會使其實際吸收劑量發生較大改變，因此及時了解這些變化是臨床醫生十分關心的問題。彈性配準為實際觀察這些變化提供了軟件支持。本文所述兩種基于Demons算法的彈性配準對宮頸癌放療患者CBCT圖像進行實驗，從配準前后差分圖的比較及表 1兩種算法的評價參數可見，兩種配準方法均取得較好的配準效果,但耗時均較長。其中Demons算法的評價參數均好于加速Demons算法，但其配準時間是加速Demons算法的1.9倍。對于細微差別處，兩種模型仍顯得精度不夠。在放療過程中，要獲得器官真實的累加劑量，要求形變的精確性非常高。怎樣將二者結合起來，更好地模擬器官運動，為患者提供更快更精確的個體化放療是未來研究的熱點。
