越來越多的醫療設備能采集人體不同特征數據并形成三維圖像。在臨床應用中, 通常需要將多幅含有不同重要信息的源圖像融合到一幅圖像中以輔助治療。而傳統的圖像融合方法主要針對二維圖像, 這些方法直接應用于三維數據會導致第三維信息的損失。本文結合近期出現的超小波變換——三維帶限剪切波變換和四組傳統融合準則, 提出新的三維磁共振圖像融合方法。最后通過四組人腦T2*和磁量圖(QSM)源數據, 將本文的方法分別與基于二維和三維的小波、雙樹復數小波的融合方法進行對比。實驗顯示基于三維變換的融合方法性能普遍優于基于二維變換的方法; 三維方法中, 剪切波克服了傳統小波變換缺乏表達方向的缺陷, 能有效地提高傳統圖像融合算法的性能, 因此, 本文的方法質量指標最高。
引用本文： 段昶, 汪學剛, 王洪, 王帥. 基于三維帶限剪切波變換的磁共振圖像融合. 生物醫學工程學雜志, 2015, 32(1): 181-186, 196. doi: 10.7507/1001-5515.20150033 復制
引言
醫學圖像融合是圖像融合的一種，已經研究了數十年，許多方法也被廣泛地應用于臨床診斷中[1-2]。融合是指將不同設備[如X線計算機斷層攝影(computer tomography，CT)、磁共振成像(magnetic resonance imaging，MRI)等]采集的重要信息提取并合并成為一幅圖像。不同設備或同類設備不同配置所生成的圖像中包含的重要信息是不同的，有些信息具有相似性，但大部分信息是互補的。CT圖像能提供人體稠密、堅硬組織的信息，而MRI圖像則主要提供軟組織的信息，如T2*能提供組織弛豫時間的對比信息，磁量圖(quantitative susceptibility mapping，QSM)可提供由多種內部磁性生物標記物或鐵、鈣、釓等對比劑引起的磁敏感對比信息[3]。通常而言，圖像融合需要先將源圖像配準，而少數情況如T2*和QSM圖像產生于同一次掃描數據，即二者已經完全配準了。
當前醫學圖像融合研究主要考慮的是二維圖像的處理，但是現在許多醫療設備都能直接采集和生成三維圖像。三維圖像中每個點的灰度值不僅與同層鄰近點互相關，也與相鄰層中的鄰近點互相關。傳統的二維融合方法直接應用于三維數據會導致第三維相關信息的丟失，因此有必要研究能直接處理三維圖像的融合方法。
融合算法可以在空域或變換域進行處理。空域中，融合圖像通常是源數據的加權平均，這樣的方法簡單易于實現，但融合圖像質量不高，邊緣等信息都有損失。變換域方法遵循以下步驟：源圖像變換到變換域，按融合準則對圖像系數進行處理，得到融合后的系數，最后將系數變回到空域，輸出即為融合圖像。這類算法中研究重點主要集中在兩點：變換的選取，融合準則的設計。許多多尺度變換都能應用于融合算法中，如離散小波變換(discrete wavelet transform，DWT)、雙樹復數小波變換(dual-tree complex wavelet transform，DTCWT)、曲波(curvelet)[4]、剪切波(shearlet)[5]等等。
剪切波變換是近幾年提出并逐步成熟的對多維數據高效表示的多尺度變換。實際上，針對小波變換缺乏對邊緣等方向性特征稀疏表示的缺點，已提出和研究了許多其他的多尺度變換。但剪切波變換是所有方法中唯一同時擁有以下優點的工具：只有一個或有限個產生函數集合，能幾乎最優地表示高維數據，對連續和離散數據采用相同方式處理，擁有緊支實現等。剪切波變換已廣泛應用于圖像處理中，如去噪[5]、緣檢測[6]、增強[7]等。剪切波也同樣適用于圖像融合。
本文提出基于三維帶限剪切波(band limited shearlet transform，BLST)的醫學圖像融合方法。融合準則方面，本文采用了最大點模值(max point’s modulus，MPM)、最大區域能量(maximum regional energy，MRE)[8]、最大絕對差值(maximum absolute difference，MAD)和最大改進拉普拉斯能量(maximum sum of modified Laplacian energy，MSMLE)[4] 4個準則。這4個準則在二維情況下是比較高效的融合準則，本文將其擴展到三維。
1 剪切波基礎理論
三維剪切波是二維剪切波的擴展，能幾乎最有效地表示三維數據中面狀(二維)和線狀(一維)奇異性特征。設：
| $ \begin{align} & {{A}_{{{2}^{j}}}}=\left(\begin{matrix} {{2}^{j}} & 0 & 0 \\ 0 & {{2}^{j/2}} & 0 \\ 0 & 0 & {{2}^{j}} \\ \end{matrix} \right), {{{\tilde{A}}}_{{{2}^{j}}}}=\left(\begin{matrix} {{2}^{j/2}} & 0 & 0 \\ 0 & {{2}^{j}} & 0 \\ 0 & 0 & {{2}^{j}} \\ \end{matrix} \right), \\ & {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{A}}}_{{{2}^{j}}}}=\left(\begin{matrix} {{2}^{j}} & 0 & 0 \\ 0 & {{2}^{j}} & 0 \\ 0 & 0 & {{2}^{j/2}} \\ \end{matrix} \right), {{S}_{R}}=\left(\begin{matrix} 1 & {{k}_{1}} & {{k}_{2}} \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{matrix} \right), \\ & {{{\tilde{S}}}_{k}}=\left(\begin{matrix} 1 & 0 & 0 \\ {{k}_{1}} & 1 & {{k}_{2}} \\ 0 & 0 & 1 \\ \end{matrix} \right), {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{S}}}_{k}}=\left(\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ {{k}_{1}} & {{k}_{2}} & 1 \\ \end{matrix} \right)\\ \end{align} $ |
其中j∈ Z。Mc=diag(c1, c2, c2), =diag(c2, c1, c2), c=diag(c2, c2, c1), c1＞0, c2＞0。三維剪切波可被定義為
| $ \begin{matrix} SH\left(\Phi, \psi, \tilde{\psi }, \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }, c \right)=\Phi \left(\Phi; {{c}_{1}} \right)\cup \Psi \left(\psi; c \right)\cup \tilde{\Psi }\left(\tilde{\psi }; c \right)\cup \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\Psi }\left(\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }; c \right), 214\Phi \left(\Phi; {{c}_{1}} \right)=\left\{ {{\Phi }_{m}}=\Phi \left(\cdot-m \right):m\in {{c}_{1}}{\boldsymbol{{\text{Z}}}^{3}} \right\}, \\ \Psi \left(\psi; c \right)=\left\{ {{\psi }_{j, k, m}}={{2}^{j}}\psi \left({{S}_{k}}{{A}_{{{2}^{j}}}}\cdot-m \right)\ge 0, \left| k \right|\le \left\lceil {{2}^{j/2}} \right\rceil, m\in {{M}_{c}}{\boldsymbol{{\text{Z}}}^{3}} \right\}, \\ \tilde{\Psi }\left(\tilde{\psi }; c \right)=\left\{ {{{\tilde{\psi }}}_{j, k, m}}={{2}^{j}}\tilde{\psi }\left({{{\tilde{S}}}_{k}}{{{\tilde{A}}}_{{{2}^{j}}}}\cdot-m \right)\ge 0, \left| k \right|\le \left\lceil {{2}^{j/2}} \right\rceil, m\in {{{\tilde{M}}}_{c}}{{\text{Z}}^{3}} \right\} \\ \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\Psi }\left(\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }; c \right)=\left\{ {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }}}_{j, k, m}}={{2}^{j}}\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }\left({{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{S}}}_{k}}{{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{A}}}_{{{2}^{j}}}}\cdot-m \right)\ge 0, \left| k \right|\le \left\lceil {{2}^{j/2}} \right\rceil, m\in {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{M}}}_{c}}{\boldsymbol{{\text{Z}}}^{3}} \right\}, j\in {{\text{N}}_{0}}, k\in {\boldsymbol{{\text{Z}}}^{2}} \\ \end{matrix}$ |
對于頻域實現的帶限剪切波，其生成函數(ξ1)= 1(ξ1)2(ξ2/ξ1) 2(ξ3/ξ1)，其中Ψ1是一個帶限小波，而Ψ2是尺度函數，表示Ψ的傅里葉變換。其頻域劃分情況由圖 1表示，關于帶限剪切波的更多知識可以參考文獻[9-10]。
圖1
金字塔三維剪切波頻域劃分
Figure1.
Partition of frequency domain for 3D pyramid adapted shearlet transform
2 融合方法
首先介紹融合準則。本文算法采用的融合準則如下：融合后的低頻系數為源數據低頻系數的均值，而高頻部分采用4個準則，分別是MPM、MRE[8]、MAD和MSMLE準則[4]，其計算公式為
| $ {{C}_{f}}=\left\{ \begin{align} & {{C}_{a}}, {{S}_{a}}\ge {{S}_{b}} \\ & {{C}_{b}}, {{S}_{a}} < {{S}_{b}} \\ \end{align} \right., $ |
其中MPM準則中St=abs(Ct)，t∈{a, b}；MRE準則中St=σt=，t∈{a, b}，Ω為一個區域，
融合算法流程圖如圖 2所示，兩幅源圖像分別用Va和Vb表示，經過三維帶限剪切波正變換后得到變換域系數Ca和Cb，通過融合準則，得到融合系數Cf，最后再通過三維帶限剪切波逆變換，得到最終的融合圖像Vf。
圖2
本文提出的三維圖像融合算法的流程圖
Figure2.
Steps of proposed 3D image fusion method
3 實驗分析
本節通過4組人腦磁共振T2*和QSM圖像數據對本文提出算法進行性能分析和比較。所有人腦的研究數據經過了學院審查委員會認可。MR掃描使用GE Signa HDxt 3.0T系統，采用8通道頭線圈。掃描序列使用3D T2*加權的多回波梯度回波序列，掃描參數為FA=20°；TR=57 ms；TE數=8；第一個TE=5.7 ms；TE間隔(ΔTE)=6.7 ms；BW=±41.67 kHz；FOV=24 cm；掃描分辨率0.57 mm×0.75 mm×2 mm。三維T2*幅度和QSM使用非線性形狀限定偶極子反演方法來重建，并且插值到128×128×128。因為在QSM生成過程中，腦組織之外的磁場被噪聲破壞，信息沒有意義，所以QSM的區域使用一個掩膜(Mask)來裁剪。該掩膜使用了一個大腦裁剪工具[11]來生成。
由于源數據為三維數據，傳統的二維算法只在x、y兩個方向上分別對每層數據進行融合，而三維融合方法則能直接對整個數據源進行融合。兩類算法的主要差別在于z軸信息的保持。對于z軸保持可通過幀間差分互信息量(inter frame difference mutual information，IFD_MI)來衡量[12-13]。本文采用的數據由于掩膜的存在，不能像文獻[12-13]那樣直接通過圖像中所有點進行計算，計算指標時只使用掩膜內的點：設上述文獻中互信息量用MIi(Dai, Dbi, Dfi)表示，Dai, Dbi, Dfi表示源圖像和融合圖像沿z軸方向的差分圖像，Dti=V ti+1 -Vti，t∈{a, b, f}，i表示序號。則本文的幀間差分互信息量用式(2)進行計算，即
| $ \text{IDF }\!\!\_\!\!\text{ MI}=\frac{1}{N-1}\sum\limits_{i=1}^{N-1}{\left(M{{I}_{i}}\left(D_{a}^{i}, D_{b}^{i}, D_{f}^{i} \right)\cdot \left(Mas{{k}^{i+1}}\cup Mas{{k}^{i}} \right)\right), } $ |
其中N為數據第三維層數量(128)，·表示點乘，Mask表示掩膜，i表示序號。
從圖 3可以看出，二維變換的幀間差分圖像存在許多與源數據不一致的區域，不同變換和準則得到的結果也有很大區別。而三維變換的結果則與源圖像非常相似，不同變換和準則的結果也幾乎沒有區別。這樣的差異也可以從表 1看出，三維變換的IFD_MI值明顯高于二維變換。而基于三維剪切波加MPM準則的方法擁有最高的IFD_MI值。由此可以得出結論，基于三維變換的方法相對基于二維變換的方法而言，具有更好的第三維信息保持能力。
圖3
二維和三維算法的幀間差分圖像
(a)源圖像幀間差分圖；(b)圖例；(c)2D DWT；(d) 2D DTCWT；(e) 3D BLST；(f) 3D DWT；(g) 3D DTCWT
Figure3. Inter frame difference image for 2D and 3D methods(a) inter frame difference image for source; (b) legend; (c) 2D DWT; (d) 2D DTCWT; (e) 3D BLST; (f) 3D DWT; (g) 3D DTCWT
表1
4組數據的IFD_MI
Table1.
IFD_MI for the four subjects
圖 4從左到右顯示了冠狀面、軸狀面、矢狀面各一層圖像的源圖像和不同方法的處理結果圖像。從圖中可以看出融合圖像中都保持了兩幅源圖像中的重要特征：相對T2*圖像，融合圖像的亮度更一致，而紋理信息更多來自QSM。不同方法的融合圖像的視覺效果幾乎一致，單純通過目測很難判定算法的優劣，還需采用性能評價指標。本文采用了互信息量(mutual information，MI)和QAB｜F[14]來評價。在文獻[14]中，QAB｜F只給出二維的情況，本文將其擴展到三維。同樣，所有的值均是基于掩膜內的數據計算得到。表 2列出4組數據的MI和QAB｜F。通過觀察可以發現，除了第二組一個指標外，指標最高的方法是基于三維帶限剪切波變換和MRE融合準則的方法。
圖4
軸狀面、冠狀面、矢狀面源圖像和不同方法的融合圖像
(a)軸狀面、冠狀面、矢狀面源圖像(左：T2*圖；右：QSM圖)；(b)3D BLST波變換融合結果(左上：MPM；右上：MRE；左下：MAD；右下：MSMLE)；(c)其他變換與最大區域能量準則的融合結果(左上：2D DWT；右上：2D DTCWT；左下：3D DWT；右下：3D DTCWT)
Figure4. Sagittal, coronal, axial source images and fused slices of different methods(a) sagittal, coronal, axial source images (Left: T2* images; Right: QSM images); (b) the fused images based on 3D BLST (Up-left: MPM; Up-right: MRE; Bottom-left: MAD; Bottom-right: MSMLE); (c) the fused images based on other transform with MRE rule (Up-left: 2D DWT; Up-right: 2D DTCWT; Bottom-left: 3D DWT; Bottom-right: 3D DTCWT)
表2
4組數據的質量指標
Table2.
Performance indices for four subjects
4 結論
現代醫療設備能提供人體結構三維圖像，而傳統的二維圖像融合算法如果直接應用于三維圖像，會導致第三維信息的損失。本文提出了基于帶限剪切波變換的三維圖像融合方法而言，并通過4組人腦MRI圖像對不同的融合方法性能進行了評估。實驗結果表明：三維圖像融合方法相對二維圖像融合方法而言，在第三維信息的保持能力上更具優勢；相對三維小波變換、三維雙樹復數小波變換而言，本文提出的融合方法具有更高的性能指標且同樣適用于其他設備(如CT、PET等)采集并已配準的醫學圖像。
引言
醫學圖像融合是圖像融合的一種，已經研究了數十年，許多方法也被廣泛地應用于臨床診斷中[1-2]。融合是指將不同設備[如X線計算機斷層攝影(computer tomography，CT)、磁共振成像(magnetic resonance imaging，MRI)等]采集的重要信息提取并合并成為一幅圖像。不同設備或同類設備不同配置所生成的圖像中包含的重要信息是不同的，有些信息具有相似性，但大部分信息是互補的。CT圖像能提供人體稠密、堅硬組織的信息，而MRI圖像則主要提供軟組織的信息，如T2*能提供組織弛豫時間的對比信息，磁量圖(quantitative susceptibility mapping，QSM)可提供由多種內部磁性生物標記物或鐵、鈣、釓等對比劑引起的磁敏感對比信息[3]。通常而言，圖像融合需要先將源圖像配準，而少數情況如T2*和QSM圖像產生于同一次掃描數據，即二者已經完全配準了。
當前醫學圖像融合研究主要考慮的是二維圖像的處理，但是現在許多醫療設備都能直接采集和生成三維圖像。三維圖像中每個點的灰度值不僅與同層鄰近點互相關，也與相鄰層中的鄰近點互相關。傳統的二維融合方法直接應用于三維數據會導致第三維相關信息的丟失，因此有必要研究能直接處理三維圖像的融合方法。
融合算法可以在空域或變換域進行處理。空域中，融合圖像通常是源數據的加權平均，這樣的方法簡單易于實現，但融合圖像質量不高，邊緣等信息都有損失。變換域方法遵循以下步驟：源圖像變換到變換域，按融合準則對圖像系數進行處理，得到融合后的系數，最后將系數變回到空域，輸出即為融合圖像。這類算法中研究重點主要集中在兩點：變換的選取，融合準則的設計。許多多尺度變換都能應用于融合算法中，如離散小波變換(discrete wavelet transform，DWT)、雙樹復數小波變換(dual-tree complex wavelet transform，DTCWT)、曲波(curvelet)[4]、剪切波(shearlet)[5]等等。
剪切波變換是近幾年提出并逐步成熟的對多維數據高效表示的多尺度變換。實際上，針對小波變換缺乏對邊緣等方向性特征稀疏表示的缺點，已提出和研究了許多其他的多尺度變換。但剪切波變換是所有方法中唯一同時擁有以下優點的工具：只有一個或有限個產生函數集合，能幾乎最優地表示高維數據，對連續和離散數據采用相同方式處理，擁有緊支實現等。剪切波變換已廣泛應用于圖像處理中，如去噪[5]、緣檢測[6]、增強[7]等。剪切波也同樣適用于圖像融合。
本文提出基于三維帶限剪切波(band limited shearlet transform，BLST)的醫學圖像融合方法。融合準則方面，本文采用了最大點模值(max point’s modulus，MPM)、最大區域能量(maximum regional energy，MRE)[8]、最大絕對差值(maximum absolute difference，MAD)和最大改進拉普拉斯能量(maximum sum of modified Laplacian energy，MSMLE)[4] 4個準則。這4個準則在二維情況下是比較高效的融合準則，本文將其擴展到三維。
1 剪切波基礎理論
三維剪切波是二維剪切波的擴展，能幾乎最有效地表示三維數據中面狀(二維)和線狀(一維)奇異性特征。設：
| $ \begin{align} & {{A}_{{{2}^{j}}}}=\left(\begin{matrix} {{2}^{j}} & 0 & 0 \\ 0 & {{2}^{j/2}} & 0 \\ 0 & 0 & {{2}^{j}} \\ \end{matrix} \right), {{{\tilde{A}}}_{{{2}^{j}}}}=\left(\begin{matrix} {{2}^{j/2}} & 0 & 0 \\ 0 & {{2}^{j}} & 0 \\ 0 & 0 & {{2}^{j}} \\ \end{matrix} \right), \\ & {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{A}}}_{{{2}^{j}}}}=\left(\begin{matrix} {{2}^{j}} & 0 & 0 \\ 0 & {{2}^{j}} & 0 \\ 0 & 0 & {{2}^{j/2}} \\ \end{matrix} \right), {{S}_{R}}=\left(\begin{matrix} 1 & {{k}_{1}} & {{k}_{2}} \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{matrix} \right), \\ & {{{\tilde{S}}}_{k}}=\left(\begin{matrix} 1 & 0 & 0 \\ {{k}_{1}} & 1 & {{k}_{2}} \\ 0 & 0 & 1 \\ \end{matrix} \right), {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{S}}}_{k}}=\left(\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ {{k}_{1}} & {{k}_{2}} & 1 \\ \end{matrix} \right)\\ \end{align} $ |
其中j∈ Z。Mc=diag(c1, c2, c2), =diag(c2, c1, c2), c=diag(c2, c2, c1), c1＞0, c2＞0。三維剪切波可被定義為
| $ \begin{matrix} SH\left(\Phi, \psi, \tilde{\psi }, \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }, c \right)=\Phi \left(\Phi; {{c}_{1}} \right)\cup \Psi \left(\psi; c \right)\cup \tilde{\Psi }\left(\tilde{\psi }; c \right)\cup \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\Psi }\left(\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }; c \right), 214\Phi \left(\Phi; {{c}_{1}} \right)=\left\{ {{\Phi }_{m}}=\Phi \left(\cdot-m \right):m\in {{c}_{1}}{\boldsymbol{{\text{Z}}}^{3}} \right\}, \\ \Psi \left(\psi; c \right)=\left\{ {{\psi }_{j, k, m}}={{2}^{j}}\psi \left({{S}_{k}}{{A}_{{{2}^{j}}}}\cdot-m \right)\ge 0, \left| k \right|\le \left\lceil {{2}^{j/2}} \right\rceil, m\in {{M}_{c}}{\boldsymbol{{\text{Z}}}^{3}} \right\}, \\ \tilde{\Psi }\left(\tilde{\psi }; c \right)=\left\{ {{{\tilde{\psi }}}_{j, k, m}}={{2}^{j}}\tilde{\psi }\left({{{\tilde{S}}}_{k}}{{{\tilde{A}}}_{{{2}^{j}}}}\cdot-m \right)\ge 0, \left| k \right|\le \left\lceil {{2}^{j/2}} \right\rceil, m\in {{{\tilde{M}}}_{c}}{{\text{Z}}^{3}} \right\} \\ \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\Psi }\left(\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }; c \right)=\left\{ {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }}}_{j, k, m}}={{2}^{j}}\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }\left({{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{S}}}_{k}}{{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{A}}}_{{{2}^{j}}}}\cdot-m \right)\ge 0, \left| k \right|\le \left\lceil {{2}^{j/2}} \right\rceil, m\in {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{M}}}_{c}}{\boldsymbol{{\text{Z}}}^{3}} \right\}, j\in {{\text{N}}_{0}}, k\in {\boldsymbol{{\text{Z}}}^{2}} \\ \end{matrix}$ |
對于頻域實現的帶限剪切波，其生成函數(ξ1)= 1(ξ1)2(ξ2/ξ1) 2(ξ3/ξ1)，其中Ψ1是一個帶限小波，而Ψ2是尺度函數，表示Ψ的傅里葉變換。其頻域劃分情況由圖 1表示，關于帶限剪切波的更多知識可以參考文獻[9-10]。
圖1
金字塔三維剪切波頻域劃分
Figure1.
Partition of frequency domain for 3D pyramid adapted shearlet transform
2 融合方法
首先介紹融合準則。本文算法采用的融合準則如下：融合后的低頻系數為源數據低頻系數的均值，而高頻部分采用4個準則，分別是MPM、MRE[8]、MAD和MSMLE準則[4]，其計算公式為
| $ {{C}_{f}}=\left\{ \begin{align} & {{C}_{a}}, {{S}_{a}}\ge {{S}_{b}} \\ & {{C}_{b}}, {{S}_{a}} < {{S}_{b}} \\ \end{align} \right., $ |
其中MPM準則中St=abs(Ct)，t∈{a, b}；MRE準則中St=σt=，t∈{a, b}，Ω為一個區域，
融合算法流程圖如圖 2所示，兩幅源圖像分別用Va和Vb表示，經過三維帶限剪切波正變換后得到變換域系數Ca和Cb，通過融合準則，得到融合系數Cf，最后再通過三維帶限剪切波逆變換，得到最終的融合圖像Vf。
圖2
本文提出的三維圖像融合算法的流程圖
Figure2.
Steps of proposed 3D image fusion method
3 實驗分析
本節通過4組人腦磁共振T2*和QSM圖像數據對本文提出算法進行性能分析和比較。所有人腦的研究數據經過了學院審查委員會認可。MR掃描使用GE Signa HDxt 3.0T系統，采用8通道頭線圈。掃描序列使用3D T2*加權的多回波梯度回波序列，掃描參數為FA=20°；TR=57 ms；TE數=8；第一個TE=5.7 ms；TE間隔(ΔTE)=6.7 ms；BW=±41.67 kHz；FOV=24 cm；掃描分辨率0.57 mm×0.75 mm×2 mm。三維T2*幅度和QSM使用非線性形狀限定偶極子反演方法來重建，并且插值到128×128×128。因為在QSM生成過程中，腦組織之外的磁場被噪聲破壞，信息沒有意義，所以QSM的區域使用一個掩膜(Mask)來裁剪。該掩膜使用了一個大腦裁剪工具[11]來生成。
由于源數據為三維數據，傳統的二維算法只在x、y兩個方向上分別對每層數據進行融合，而三維融合方法則能直接對整個數據源進行融合。兩類算法的主要差別在于z軸信息的保持。對于z軸保持可通過幀間差分互信息量(inter frame difference mutual information，IFD_MI)來衡量[12-13]。本文采用的數據由于掩膜的存在，不能像文獻[12-13]那樣直接通過圖像中所有點進行計算，計算指標時只使用掩膜內的點：設上述文獻中互信息量用MIi(Dai, Dbi, Dfi)表示，Dai, Dbi, Dfi表示源圖像和融合圖像沿z軸方向的差分圖像，Dti=V ti+1 -Vti，t∈{a, b, f}，i表示序號。則本文的幀間差分互信息量用式(2)進行計算，即
| $ \text{IDF }\!\!\_\!\!\text{ MI}=\frac{1}{N-1}\sum\limits_{i=1}^{N-1}{\left(M{{I}_{i}}\left(D_{a}^{i}, D_{b}^{i}, D_{f}^{i} \right)\cdot \left(Mas{{k}^{i+1}}\cup Mas{{k}^{i}} \right)\right), } $ |
其中N為數據第三維層數量(128)，·表示點乘，Mask表示掩膜，i表示序號。
從圖 3可以看出，二維變換的幀間差分圖像存在許多與源數據不一致的區域，不同變換和準則得到的結果也有很大區別。而三維變換的結果則與源圖像非常相似，不同變換和準則的結果也幾乎沒有區別。這樣的差異也可以從表 1看出，三維變換的IFD_MI值明顯高于二維變換。而基于三維剪切波加MPM準則的方法擁有最高的IFD_MI值。由此可以得出結論，基于三維變換的方法相對基于二維變換的方法而言，具有更好的第三維信息保持能力。
圖3
二維和三維算法的幀間差分圖像
(a)源圖像幀間差分圖；(b)圖例；(c)2D DWT；(d) 2D DTCWT；(e) 3D BLST；(f) 3D DWT；(g) 3D DTCWT
Figure3. Inter frame difference image for 2D and 3D methods(a) inter frame difference image for source; (b) legend; (c) 2D DWT; (d) 2D DTCWT; (e) 3D BLST; (f) 3D DWT; (g) 3D DTCWT
表1
4組數據的IFD_MI
Table1.
IFD_MI for the four subjects
圖 4從左到右顯示了冠狀面、軸狀面、矢狀面各一層圖像的源圖像和不同方法的處理結果圖像。從圖中可以看出融合圖像中都保持了兩幅源圖像中的重要特征：相對T2*圖像，融合圖像的亮度更一致，而紋理信息更多來自QSM。不同方法的融合圖像的視覺效果幾乎一致，單純通過目測很難判定算法的優劣，還需采用性能評價指標。本文采用了互信息量(mutual information，MI)和QAB｜F[14]來評價。在文獻[14]中，QAB｜F只給出二維的情況，本文將其擴展到三維。同樣，所有的值均是基于掩膜內的數據計算得到。表 2列出4組數據的MI和QAB｜F。通過觀察可以發現，除了第二組一個指標外，指標最高的方法是基于三維帶限剪切波變換和MRE融合準則的方法。
圖4
軸狀面、冠狀面、矢狀面源圖像和不同方法的融合圖像
(a)軸狀面、冠狀面、矢狀面源圖像(左：T2*圖；右：QSM圖)；(b)3D BLST波變換融合結果(左上：MPM；右上：MRE；左下：MAD；右下：MSMLE)；(c)其他變換與最大區域能量準則的融合結果(左上：2D DWT；右上：2D DTCWT；左下：3D DWT；右下：3D DTCWT)
Figure4. Sagittal, coronal, axial source images and fused slices of different methods(a) sagittal, coronal, axial source images (Left: T2* images; Right: QSM images); (b) the fused images based on 3D BLST (Up-left: MPM; Up-right: MRE; Bottom-left: MAD; Bottom-right: MSMLE); (c) the fused images based on other transform with MRE rule (Up-left: 2D DWT; Up-right: 2D DTCWT; Bottom-left: 3D DWT; Bottom-right: 3D DTCWT)
表2
4組數據的質量指標
Table2.
Performance indices for four subjects
4 結論
現代醫療設備能提供人體結構三維圖像，而傳統的二維圖像融合算法如果直接應用于三維圖像，會導致第三維信息的損失。本文提出了基于帶限剪切波變換的三維圖像融合方法而言，并通過4組人腦MRI圖像對不同的融合方法性能進行了評估。實驗結果表明：三維圖像融合方法相對二維圖像融合方法而言，在第三維信息的保持能力上更具優勢；相對三維小波變換、三維雙樹復數小波變換而言，本文提出的融合方法具有更高的性能指標且同樣適用于其他設備(如CT、PET等)采集并已配準的醫學圖像。
