-
3-D剪切波系统是由各向异性膨胀矩阵和剪切矩阵组成[10],它们分别控制尺度和不同尺度的方向,设各向异性膨胀矩阵为:
$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{A}}_{\alpha ,{2^j}}} = \left[ {\begin{array}{*{20}{c}} {{2^j}}&0&0\\ 0&{{2^{\alpha j/2}}}&0\\ 0&0&{{2^{\alpha j/2}}} \end{array}} \right]}\\ {{{\mathit{\boldsymbol{\tilde A}}}_{\alpha ,{2^j}}} = \left[ {\begin{array}{*{20}{c}} {{2^{\alpha j/2}}}&0&0\\ 0&{{2^j}}&0\\ 0&0&{{2^{\alpha j/2}}} \end{array}} \right]}\\ {{{\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over A} }}}_{\alpha ,{2^j}}} = \left[ {\begin{array}{*{20}{c}} {{2^{\alpha j/2}}}&0&0\\ 0&{{2^{\alpha j/2}}}&0\\ 0&0&{{2^j}} \end{array}} \right]} \end{array}} \right. $
(1) 剪切矩阵为:
$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{S}}_k} = \left[ {\begin{array}{*{20}{c}} 1&{{k_1}}&{{k_2}}\\ 0&1&0\\ 0&0&1 \end{array}} \right]}\\ {{{\mathit{\boldsymbol{\tilde S}}}_k} = \left[ {\begin{array}{*{20}{c}} 1&0&0\\ {{k_1}}&1&{{k_2}}\\ 0&0&1 \end{array}} \right]}\\ {{{\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over S} }}}_k} = \left[ {\begin{array}{*{20}{c}} 1&0&0\\ 0&1&0\\ {{k_1}}&{{k_2}}&1 \end{array}} \right]} \end{array}} \right. $
(2) 式中, j∈Z为尺度参量,α∈(0, 2)为各向异性程度,k=(k1, k2)∈Z2为剪切方向,Z为整数集。移位采样晶格为下式:
$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{M}}_c} = {\mathop{\rm diag}\nolimits} \left( {{c_1},{c_2},{c_2}} \right)}\\ {{{\mathit{\boldsymbol{\tilde M}}}_c} = {\mathop{\rm diag}\nolimits} \left( {{c_2},{c_1},{c_2}} \right)}\\ {{{\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over M} }}}_c} = {\mathop{\rm diag}\nolimits} \left( {{c_2},{c_1},{c_2}} \right)} \end{array}} \right. $
(3) 式中,c1>0, c2>0为平移量。
对于ϕ,ψ, $\widetilde \psi $, $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \psi }$∈L2(R2), 其中ϕ为剪切波尺度函数,ψ, $\widetilde \psi $,$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \psi }$为剪切波函数,R表示实数集,L2(R2)表示实数域R上2维均方可积函数集合,对每一个尺度j, 有α=(αj)j,αj∈(0, 2),且c=(c1, c2)∈(R+)2,R+为正实数,m为位置参量,通用剪切系统SH(ϕ, ψ, $\widetilde \psi $, $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \psi }$; α, c)被定义为:
$ \begin{array}{*{20}{c}} {{\rm{SH}}(\phi ,\psi ,\tilde \psi ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \psi } ;\alpha ,c) = \varPhi (\phi ;{c_1}) \cup }\\ {\varPsi (\psi ;\alpha ,c) \cup \tilde \varPsi (\tilde \psi ;\alpha ,c) \cup \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \varPsi } (\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \psi } ;\alpha ,c)} \end{array} $
(4) 式中,
$ \varPhi (\phi ;{c_1}) = \{ {\phi _m} = \phi ( \cdot - {c_1}m),m \in {{\bf{Z}}^2})\} $
(5) $ \begin{array}{*{20}{c}} {\varPsi (\psi ;\alpha ,c) = \left\{ {{\psi _{j,k,m}} = {2^{\frac{{{\alpha _j} + 1}}{4}j}}\psi \left( {{\mathit{\boldsymbol{S}}_k}^{\rm{T}}{\mathit{\boldsymbol{A}}_{{\alpha _j},{2^j}}} \cdot - {\mathit{\boldsymbol{M}}_c}m} \right),} \right.}\\ {\left. {\left( {j \ge 0,|k| \le \left\lceil {{2^{j\left( {{\alpha _j} - 1} \right)/2}}} \right\rceil ,m \in {{\bf{Z}}^2}} \right)} \right\}} \end{array} $
(6) $ \begin{array}{*{20}{c}} {\tilde \varPsi (\tilde \psi ;\alpha ,c) = \left\{ {{{\tilde \psi }_{j,k,m}} = {2^{\frac{{{\alpha _j} + 1}}{4}j}}{\rm{ }}\tilde \psi \left( {{\mathit{\boldsymbol{S}}_k}^{\rm{T}}{\rm{ }}{{\mathit{\boldsymbol{\tilde A}}}_{{\alpha _j},{2^j}}} \cdot - {{\mathit{\boldsymbol{\tilde M}}}_c}m} \right),} \right.}\\ {\left. {\left( {j \ge 0,|k| \le \left\lceil {{2^{j\left( {{\alpha _j} - 1} \right)/2}}} \right\rceil ,m \in {{\bf{Z}}^2}} \right)} \right\}} \end{array} $
(7) $ \begin{array}{*{20}{c}} {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \varPsi } (\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \psi } ;\alpha ,c) = \left\{ {{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \psi } }_{j,k,m}} = {2^{\frac{{{\alpha _j} + 1}}{4}j}}\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over \psi } \left( {{\mathit{\boldsymbol{S}}_k}^{\rm{T}}{{\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over A} }}}_{{\alpha _j},{2^j}}} \cdot - {{\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over M} }}}_c}m} \right)} \right.,}\\ {\left. {\left( {j \ge 0,|k| \le \left\lceil {{2^{j\left( {{\alpha _j} - 1} \right)/2}}} \right\rceil ,m \in {{\bf{Z}}^2}} \right)} \right\}} \end{array} $
(8) 对于ShearLab 3D中$\widehat \psi $(ξ)是一种不可分剪切波[11],在L2(R2)空间是一个紧支撑的剪切波框架, 能最优表达高维数据的空间特征其生成函数为:
$ \begin{array}{*{20}{l}} {\hat \psi (\xi ) = \left( {P\left( {\frac{{{\xi _1}}}{2},{\xi _2}} \right){{\hat \psi }_1}({\xi _1}){{\hat \phi }_1}({\xi _2})} \right) \cdot }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {P\left( {\frac{{{\xi _1}}}{2},{\xi _3}} \right){{\hat \phi }_1}({\xi _3})} \right)} \end{array} $
(9) 式中,函数ϕ1和ψ1是紧密支持的函数,2维扇形滤波函数P都满足一定条件[11]。
-
首先从ADNI(Alzheimer’s Disease Neuroimaging Initiative)数据库中下载AD患者的PET和MRI脑功能影像,分别在统计参量图(statistical parametric mapping, SPM)中对源图像经过头动校正、配准、归一化处理得到两个待处理对象图像A, B。然后利用ShearLab 3D分解分别得到低频子带Cl, A(x, y, z)和Cl, B(x, y, z)高频子带Ch, 0, A(x, y, z),Ch, 0, B(x, y, z), 其中高频子带以方差为阈值获得强边缘子带为中频子带Cm, A(x, y, z)和留下的高频子带Ch, 1, A(x, y, z), Ch, 1, B(x, y, z)。本文中提出将低中高频融合规则扩展为3维,以加权局部能量与拉普拉斯算子加权和融合低频系数并以锐化矩阵作为权重参量捕捉图像中灰度跳变区域,以绝对值为活动度量融合中频子带,以3个3维扩展的底层特征加权融合高频系数。最后,经过ShearLab 3D逆变换得到融合结果PET/ MRI图像F,流程如图 1所示。
-
ShearLab 3D分解生成的低频子带不但包含图像绝大多数的能量,也含有图像的纹理特征。因此低频子带融合规则的选择会影响融合结果的好坏,常用的加权平均法使融合结果轮廓不完整清晰、对比度降低。像素取大法的融合结果中MRI图像的脑沟回等细节特征不足。
本文中考虑到每个体素点的空间信息将2维局部拉普拉斯能量加权[15]扩展到3维形式。设EWLE, s(x, y, z)为低频子带位置(x, y, z)的加权局部能量(weight local energy,WLE),s∈{A, B}分别表示待融合的PET和MRI图像。为了突出了MRI图像边缘信息,本文中将W设置为(2r+1)×(2r+1)×(2r+1)的锐化矩阵,r是W的半径,即:
$ \begin{array}{*{20}{c}} {{E_{{\rm{WLE}},\mathit{\boldsymbol{s}}}}(x,y,z) = \sum\limits_{m = - r}^r {\sum\limits_{n = - r}^r {\sum\limits_{p = - r}^r \mathit{\boldsymbol{W}} } } (m,n,p) \times }\\ {{C_{1,s}}{{(x + m,y + n,z + p)}^2}} \end{array} $
(10) 取r=1,W为3×3×3的矩阵,是以拉普拉斯变形模板M为基础,考虑到系数之间的相关性并以中心点到邻域的欧式距离为依据,将矩阵W的系数中距离为1的设置为-2,将距离为$\sqrt 2 $的设置为1,矩阵系数中心选取锐化程度居中的-4。
$ \mathit{\boldsymbol{M}} = \left[ {\begin{array}{*{20}{l}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 1}&{ - 2}&{{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 1}\\ { - 2}&{{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 4}&{ - 2}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 1}&{ - 2}&{{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 1} \end{array}} \right] $
(11) $ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{W}}_{3 \times 3 \times 3}} = \\ \left[ {\begin{array}{*{20}{c}} 1&1&1\\ 1&{ - 2}&1\\ 1&1&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} 1&{ - 2}&1\\ { - 2}&{ - 4}&{ - 2}\\ 1&{ - 2}&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} 1&1&1\\ 1&{ - 2}&1\\ 1&1&1 \end{array}} \right] \end{array} $
(12) 设LWSTML, s(x, y, z)是低频子带位置(x, y, z)处二十六邻域的改进拉普拉斯算子的加权和(weight sum of twenty six-neighborhood based modified laplacian,WSTML)。LTML, s[15]为改进的拉普拉斯算子(twenty six-neighborhood based modified laplacian,TML),将相邻点的权重设为欧氏距离。为了充分利用相邻体素间的关系, 将LTML, s扩展到3维,即:
$ \begin{array}{*{20}{c}} {{L_{{\rm{WSTML}},\mathit{\boldsymbol{s}}}}(x,y,z) = \sum\limits_{m = - r}^r {\sum\limits_{n = - r}^r {\sum\limits_{p = - r}^r \mathit{\boldsymbol{W}} } } (m,n,p) \times }\\ {{L_{{\rm{TML}},\mathit{\boldsymbol{s}}}}{{(x + m,y + n,z + p)}^2}} \end{array} $
(13) $ \begin{array}{*{20}{c}} {{L_{{\rm{TML}},\mathit{\boldsymbol{s}}}}(x,y,z) = \sum\limits_{s = - r}^r {\sum\limits_{t = - r}^r {\sum\limits_{l = - r}^r {\left\{ {{{[x - (x + i)]}^2} + } \right.} } } }\\ {{{\left. {{{[y - (y + g)]}^2} + {{[z - (z + e)]}^2}} \right\}}^{ - 1/2}} \times }\\ {\left| {{C_{1,\mathit{\boldsymbol{s}}}}(x,y,z) - {C_{1,\mathit{\boldsymbol{s}}}}(x + i,y + q,z + e)} \right|} \end{array} $
(14) 式中, i,q,e都为整数且不能同时为0;s, W, r同上。
位置(x, y, z)的低频融合系数Fl(x, y, z)为下式,选取局部拉普拉斯能量大的系数为融合图像的低频系数:
$ \begin{array}{*{20}{l}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {F_1}(x,y,z) = }\\ {\left\{ {\begin{array}{*{20}{l}} {{C_{1,\mathit{\boldsymbol{A}}}}(x,y,z),\left( {{E_{{\rm{WLE}},\mathit{\boldsymbol{A}}}} + {L_{{\rm{WSTML}},\mathit{\boldsymbol{A}}}} \ge {E_{{\rm{WLE}},\mathit{\boldsymbol{B}}}} + {L_{{\rm{WSTML}},\mathit{\boldsymbol{B}}}}} \right)}\\ {{C_{1,\mathit{\boldsymbol{B}}}}(x,y,z),({\rm{ otherwise }})} \end{array}} \right.} \end{array} $
(15) -
ShearLab 3D分解后高频子带包含图像的边缘和纹理结构,为了更好地融合图像的边缘信息在高频子带中提取包含源图像强边缘信息中频子带。由于方差可以判断图像的纹理、边缘、平滑区域,其中纹理区方差最大,边缘次之,平滑最小[16]。因此以方差为阈值得到中频子带,且方差的选取与图像特点有关。
设μs与$\overline {{C_{{\rm{h}}, 0, \mathit{\boldsymbol{s}}}}} $分别ShearLab 3D分解后高频子带的方差和均值,s∈{A, B},m×n×p为高频子带图像大小。根据PET图像边缘模糊和MRI图像结构清晰边缘分界明显的特点,并通过大量实验手动设置阈值范围得到融合图像的客观评价指标,最终选取主客观评价最优的阈值0.03 < μA < 0.04,0.2 < μB < 0.3。
$ {\mu _\mathit{\boldsymbol{s}}} = \frac{{\sum\limits_{x = 0}^m {\sum\limits_{y = 0}^n {\sum\limits_{z = 0}^p {{{\left[ {{C_{{\rm{h}},\mathit{\boldsymbol{s}}}}(x,y,z) - \overline {{C_{{\rm{h}},0,\mathit{\boldsymbol{s}}}}} } \right]}^2}} } } }}{{m \times n \times p}} $
(16) 中频子带采取基于绝对值的活动度量[17]保留中频子带中的边缘信息,并充分考虑到空间结构信息将其扩展到3维情况。
设Cm, A(x, y, z)与Cm, B(x, y, z)为PET, MRI图像的中频子带,以下式二值化处理:
$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{G}}(x,y,z) = }\\ {\left\{ {\begin{array}{*{20}{l}} {1,\left( {\left| {{C_{{\rm{m}},\mathit{\boldsymbol{A}}}}(x,y,z)} \right| > \left| {{C_{{\rm{m}},\mathit{\boldsymbol{B}}}}(x,y,z)} \right|} \right)}\\ {0,({\rm{ otherwise }})} \end{array}} \right.} \end{array} $
(17) 式中, G(x, y, z)为二值化矩阵,当PET图像的中频子带系数绝对值大于MRI图像的中频子带系数绝对值时, G(x, y, z)=1,否则为0。
d(x, y, z)为中频子带融合的决策图,g(x, y, z)为活动指标,通过以下等式计算:
$ \mathit{\boldsymbol{d}}(x,y,z) = \left( {\begin{array}{*{20}{l}} {1,\left( {g(x,y,z){\rm{ > }}\frac{{M \times N \times P}}{2}} \right)}\\ {0,({\rm{otherwise}})} \end{array}} \right. $
(18) $ \mathit{\boldsymbol{g}}(x,y,z) = \{ ({x_0},{y_0},{z_0}) \in \varOmega |\mathit{\boldsymbol{G}}(x,y,z)|\} $
(19) 式中, Ω是以(x, y, z)为中心的M×N×P的滑动窗口。
中频子带的融合系数为Fm(x, y, z),当决策图d(x, y, z)=1时, PET图像的中频系数为融合图像的中频子带系数。
$ {F_{\rm{m}}}(x,y,z) = \left\{ {\begin{array}{*{20}{l}} {{C_{{\rm{m}},\mathit{\boldsymbol{A}}}}(x,y,z),(\mathit{\boldsymbol{d}}(x,y,z) = 1)}\\ {{C_{{\rm{m}},\mathit{\boldsymbol{B}}}}(x,y,z),({\rm{ otherwise }})} \end{array}} \right. $
(20) -
高频子带包含了MRI源图像脑沟回等细节特征,对AD病情的诊断非常重要。因此高频子带融合规则应有效的保留细节信息,局部相位一致的大部分特征为边缘状和角状特征[18]。局部对比度、局部能量可以反映对比度和亮度信息。本文中采取这3种底层特征加权[19]的融合规则,考虑到每个体素点的空间信息将其扩展到3维形式。
设PPC, s(x, y, z)为高频子带(x, y, z)处的相位一致(phase congruency,PC)的点,s∈(A, B),即:
$ {P_{{\rm{PC}},\mathit{\boldsymbol{s}}}}(x,y,z) = \frac{{\sum\limits_k {{E_{{\theta _k}}}} (x,y,z)}}{{\varepsilon + \sum\limits_n {\sum\limits_k {{A_{n,{\theta _k}}}} } (x,y,z)}} $
(21) 式中, θk是k的方向角,An, θk为第n个傅里叶分量的幅值和角度,ε=0.001是正数去除图像中的DC分量,Eθk(x, y, z)可见参考文献[19]。
SLSCM, s(x, y, z)为(x, y, z)处的局部对比度(local measure of sharpness change, LSCM),(2M+1)×(2N+1)×(2P+1)是邻域大小,SSCM, s(x, y, z)为高频子带(x, y, z)处的对比度(measure of sharpness change, SCM)的值。
$ \begin{array}{*{20}{l}} {{S_{{\rm{LSCM,}}\mathit{\boldsymbol{s}}}}(x,y,z) = \sum\limits_{a = - M}^M {\sum\limits_{b = - N}^N {\sum\limits_{t = - P}^P {{S_{{\rm{SCM,}}\mathit{\boldsymbol{s}}}}} } } (x + a,}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} y + b,z + t)} \end{array} $
(22) $ \begin{array}{l} {S_{{\rm{SCM}},\mathit{\boldsymbol{s}}}}(x,y,z) = \sum\limits_{\left( {{x_0},{y_0},{z_0}} \right) \in {\varOmega _0}} {\left[ {{C_{{\rm{h}},1,\mathit{\boldsymbol{s}}}}(x,y,z) - } \right.} \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left. {{C_{{\rm{h}},1,\mathit{\boldsymbol{s}}}}\left( {{x_0},{y_0},{z_0}} \right)} \right] \end{array} $
(23) 式中, 取M=N=P=1,Ω0是(x, y, z)的3×3×3的局部区域,(x0, y0, z0)表示Ω0内的体素点。
ELE, s(x, y, z)为高频子带(x, y, z)处的局部能量(local energy, LE),且M=N=P=1。
$ \begin{array}{*{20}{c}} {{E_{{\rm{LE}},\mathit{\boldsymbol{s}}}}(x,y,z) = }\\ {\sum\limits_{a = - M}^M {\sum\limits_{b = - N}^N {\sum\limits_{t = - P}^P {{{\left[ {{C_{{\rm{h}},1,\mathit{\boldsymbol{s}}}}(x + a,y + b,z + t)} \right]}^2}} } } } \end{array} $
(24) 用计算出来的PPC, s,SLSCM, s,ELE, s组成活跃度测量Ns; α1,β1和γ1是用来调整PPC, s,SLSCM, s,ELE, s的参量,分别为1, 2, 2[19] :
$ \begin{array}{*{20}{c}} {{N_\mathit{\boldsymbol{s}}}(x,y,z) = {{\left[ {{P_{{\rm{PC}},\mathit{\boldsymbol{s}}}}(x,y,z)} \right]}^{{\alpha _1}}} \cdot }\\ {{{\left[ {{S_{{\rm{LSCM}},\mathit{\boldsymbol{s}}}}(x,y,z)} \right]}^{{\beta _1}}} \cdot {{\left[ {{E_{{\rm{LE}},\mathit{\boldsymbol{s}}}}(x,y,z)} \right]}^{{\gamma _1}}}} \end{array} $
(25) d(x, y, z)为高频子带融合的决策图,通过(18)式计算,其中g(x, y, z)由等下式计算:
$ \begin{array}{*{20}{c}} {g(x,y,z) = \left\{ {\left( {{x_0},{y_0},{z_0}} \right) \in {\varOmega _1}\mid {N_\mathit{\boldsymbol{s}}} \ge } \right.}\\ {\left. {\max \left( {{N_\mathit{\boldsymbol{A}}}\left( {{x_0},{y_0},{z_0}} \right),{N_\mathit{\boldsymbol{B}}}\left( {{x_0},{y_0},{z_0}} \right)} \right)} \right\}} \end{array} $
(26) 式中, Ω1是以(x, y, z)为中心的M×N的滑动窗口。
高频子带的融合系数Fh(x, y, z)与(20)式相同,当决策图d(x, y, z)=1时, PET图像系数为融合图像的高频系数。
-
本文中所有实验均采用Window 8操作系统,使用MATLAB 2018a平台进行仿真实验, 且PET和MRI数据均为同一AD患者同一时期的脑影像。首先在空域和变换域中采取相同的融合规则进行比较,然后在ShearLab 3D域中将本文中的融合规则与现有的PET/MRI算法[17]进行比较。选取信息熵(entropy, IE)、平均梯度(average gradient, AG)、空间频率(spatial frequency, SF)、边缘强度(edge intensity, EI)、综合熵(comprehensive entropy, CE)作为客观评价标准, 它们的值越大融合效果越好。
-
为了验证ShearLab 3D相较于空域算法和其它变换域算法的优越性,在加权平均、像素取大、小波域和ShearLab 3D域内融合PET/MRI图像,其中小波域和ShearLab 3D域低频系数使用加权平均法,高频子带系数使用体素绝对值取大法,且都为一层分解。图 2所示为PET和MRI图像和不同算法产生的结果图像。图 2a是AD患者注射18F-FDG的PET图像,其中左上是冠状面,右上是矢状面,左下是轴状面; 图 2b表示AD患者MRI脑图像; 图 2c是空域加权平均结果; 图 2d为空域像素取大结果; 图 2e为变换域小波变换结果; 图 2f为变换域ShearLab 3D的融合结果。
从主观角度分析,空域的加权平均相对于其它方法融合图像对比度略有降低导致边缘轮廓模糊。像素取大法由于PET和MRI呈现原理不同导致融合结果丢失MRI图像纹理信息,变换域的小波变换融合图像的边缘轮廓纹理等信息都受到了块效应的影响。ShearLab 3D融合图像从很大程度上弥补了它们的不足,可以观察到脑萎缩程度及相应组织的代谢情况。
从客观角度分析,观察表 1、图 2c和图 2e,其各项评价指标数值普遍较低,尤其是图 2c中的AG, EI, SF指数过低,表示融合图像纹理特征不够清晰,灰度变化率低对比度差。图 2e中IE和CE较低融合图像信息丢失严重。图 2d中的IE, EI, CE都较好,但融合图像中只含有极少部分MRI图像的细节纹理,缺乏实用性。图 2f的AG最好,其它指标与最优相比相差甚少,相对于其它算法,ShearLab 3D轮廓清晰、细节明显更具有实用性。
-
在ShearLab 3D域将本文中提出的融合规则与参考文献[17]中提出的PET/MRI的融合规则进行对比。ShearLab 3D为一层分解,图 3a是参考文献[17]中融合规则的结果,图 3b是本文中融合算法结果。从视觉角度看, 两幅图像均保留了原图像的主要信息,但是图 3a中的脑沟回处边缘轮廓不明显且对比度较低,无法判断脑萎缩状况,降低了医疗诊断的实用性。本文中算法在脑沟回处边缘轮清晰,并准确地判断海马体葡萄糖代谢降低萎缩严重,右侧额叶与右侧颞叶葡萄糖代谢降低发生病变。
Figure 3. PET/MRI images with different fusion rules in ShearLab 3D domain a—reference [17] fusion results b—algorithm fusion results in this paper
从客观指标分析观察表 2,本文中算法融合结果的IE值略低于参考文献[17]中的算法,但SF, AG, CE的评价指数都高于对比算法且EI有明显提高,表明融合结果的体素活跃程度高,对微小细节的表示度高,图像清晰。综合主客观评价, 本文中算法的融合结果边缘纹理特征丰富,清晰水平更好。
基于ShearLab 3D变换的3维PET/MRI图像融合
3-D PET/MRI image fusion based on ShearLab 3D transform
-
摘要: 为了解决阿尔茨海默病3维正电子发射断层成像(PET)与核磁共振成像(MRI)相同位置强度不同问题,同时保留MRI大脑皮质、脑回沟、海马体等的萎缩情况, 首先将两幅源图像在统计参量图(SPM)中预处理,再利用3维剪切波系统(ShearLab 3D)最优表达高维数据的能力对图像分解,生成低高频子带,并以方差作为阈值将高频子带分为中高频子带。低频子带使用3维扩展的加权局部能量与改进拉普拉斯算子的加权和加权的融合规则,并引入锐化矩阵为权重参量,使融合图像边缘清晰; 中频子带以绝对值为活动度量增强图像的边缘信息; 高频子带以3个3维底层特征加权融合规则增强图像的细节特征。最后,利用ShearLab 3D逆变换获得PET/MRI图像。结果表明,ShearLab 3D变换的融合结果整体优于空域和小波变换;ShearLab 3D方法中将不同融合规则对比分析,该算法融合结果的平均梯度、空间频率、边缘强度、综合熵分别提高了11.09%,22.58%,152.68%,0.58%,解决了边缘模糊、细节不清晰的问题。该研究为PET/MRI图像融合提供了参考。
-
关键词:
- 图像处理 /
- 图像融合 /
- PET/MRI /
- ShearLab 3D软件 /
- 拉普拉斯算子的加权和
Abstract: In order to solve the problem that the different intensity of the same position of the 3-D positron emission tomography(PET) and the magnetic resonance imaging(MRI) image of Alzheimer's disease and retain the MRI atrophy of the cerebral cortex, cerebral sulcus, hippocampus, etc. Two images were firstly pre-processed in SPM to obtain two images. Then, using ShearLab 3D transform to process the advantages of high-dimensional data to decompose to obtain low and high-pass subbands. The high frequency subband was divided into intermediate and high frequency subband with the variance as threshold. The fusion principle of low-pass subbands was based on the method of three-dimensional extended weighted local energy and weighted sum of modified Laplacian based on 26 neighborhoods. The sharpening operator was introduced as a weight parameter to make the edge of the fused image clear. Intermediate subband enhances the edge information with absolute value activity. The high-pass subbands were combined with three three-dimensional low-level visual features to enhance the detailed features of the image. Finally, PET/MRI fusion images were obtained using ShearLab 3D inverse transform. The results show that the fusion result of ShearLab 3D transform is better than the spatial algorithm and wavelet transform as a whole. In the ShearLab 3D method, different fusion rules are compared and analyzed. The average gradient, spatial frequency, edge strength, and comprehensive entropy of the fusion result of this algorithm were improved by 11.09%, 22.58%, 152.68%, and 0.58%, respectively. It solves the problems of blurred edges and unclear details in fusion image. This study provides a reference for PET/MRI image fusion. -
Figure 3. PET/MRI images with different fusion rules in ShearLab 3D domain a—reference [17] fusion results b—algorithm fusion results in this paper
Table 1. Objective evaluation of PET/MRI images in different domains
-