基于不可分离MRA的小波图像重建

                          作者：金国英, 罗戎蕾, 汪元美

【摘要】  本文提出了一种基于不可分离MRA的小波重建算法。算法将投影数据进行二通道的小波分解，在分解抽样中，使用不同的抽样方法。算法结果，相比通常所使用的基于可分离MRA的小波重建算法基本相同的局部重建特性，重建质量和运算速度都优于后者。

【关键词】  图像重建； 小波变换； 不可分离MRA

    Abstract: A wavelet reconstruction algorithm based on non?separable MRA has been proposed. In the algorithm  the projection data were disposed with two channel wavelet. This algorithm not only has the same local reconstruction characteristics with the algorithm based on separable MRA, but have better image quality and can be obtained by shorter working time.

    Key words: tomography;   wavelet transform;   non?separable MRA

    引言

    利用二维小波理论进行图像重建的目的为了利用小波的局部性，即使用局部数据就可以得到局部图像。这样可以减少X射线对人体的射入。通常在二维小波变换时，对二维图像进行行和列的依次滤波［1］，即将假定二维信号是关于自变量x和y可分离的，而实际二维信号种大多时不宜分开处理的。基于不可分离MRA的小波去噪算法将投影数据进行二通道的小波分解，从而直接得到小波的近似系数和细节系数，这些系数再经过逆小波变换得到最终的重建图像。对算法进行仿真实验，其结果不但具有基于可分离MRA（Mualtiresolution Analysis—多分辨率分析）的小波重建算法基本相同的局部重建特性，重建质量和运算速度都优于后者。

    1  重建算法

    首先，我们可以写双尺度方程为［2］：Φ(x)=?k∈Z2h0(k)Φ(Dx－k)，

    Ψ(x)=?k∈Z2h1(k)Ψ(Dx－k)，(1)其中，对于固定的尺度j，{Φj,k,k∈Z2}构成Vj空间的一个非正交基，{Ψij,k,i=1…(N－1)}构成空间Wj的非正交基，而对于抽样矩阵D，我们定义N=|detD|，则Φj,k(x)=|detD|－(j/2)Φ(D－jx－k)，其中，h0(k)和h1(k)分别为二维低通滤波器和二维高通滤波器，可以通过一维滤波器利用McClellan变换得到［2］：用二维零相位滤波器代替(cosω)，从而二维的频率响应为：H(ω)=?Nn=1a(n)Tn［F(ω1,ω2)］.对于梅花五点抽样方法：D=11

    1－1,（2）

    F(ω1,ω2)=12(cos(ω1)+cos(ω2)).    为了避免混叠现象，我们取h1=0(－z1,－z2),  1=－h0(－z1,－z2),（3）其中，h0、h1表示分解低通和高通分解滤波器，0、表示低通和高通重构滤波器。对于尺度函数对{Φ,Φ～}及N-1对小波函数{Ψ,Ψ～}，函数f∈L2(R2)可以分解为双正交小波序列［4］：f=?k∈Z2〈f,Φj,k〉Φ～j,k(x)

    +?N－1i=1?j≤J?k∈Z2〈f,Ψij,k〉Ψ～ij,k(x),(4)其中，〈f,Ψj,k〉是f在子空间Vj上的投影，表示f在Vj上的离散近似系数。

    从式(2)可得： |detD|=2，所以，在分辨率为j时，原图像f可分解为：近似图像Ajf［n］和细节图像Djf［n］：Ajf［n］=〈f,Φj,n〉=(f?Φ～j,0)(Djn),

    Djf［n］=〈f,Φj,n〉=(f?Ψ～j,0)(Djn),(5)其中，Φ～(x)=Φ(－x)，  Ψ～(x)=Ψ(－x),则对于滤波反投影算法［5］f=R#ΛRθf，可以写为：Ajf［n］=(f?R#ΛRθΦ～j,0)(Djn)

    =R#(Rθf?ΛRθΦ～j,0)(Djn),（6）

    Djf［n］=(f?R#ΛRθΨ～j,0)(Djn)

    =R#(Rθf?ΛRθΨ～j,0)(Djn).（7）    不可分离的二维MRA和可分离的二维MRA在形式上是相同的。在基于可分离的二维MRA算法中：Φj,0(x)=?2j(x)?2j(y)，

    ?2j(x)=2－j/2?(2－jx).    利用式(1)我们写式(5)如下，Ajf［n］=f*hj0(Djn),

    hj0［n］=?k∈Z2h0(k)h［n－Dj－1k］,

    Djf［n］=f*hj1(Djn),

    hj1［n］=?k∈Z2h1(k)hj－1［n－Dj－1k］.(8)将式(8)代入式(6)可得： Ajf［n］=∫π0(pθ?kjθ,h1)(Djn,Θ)dθ,其中，Kθ(w)=|w|×H0wcosθ,S0Spwsinθ，Hj0是hj0［n］（式 (8)）的二维傅立叶变换。

    算法可以写为以下步骤：

    (1) 对投影函数pθ(n)用kjθ,h1和kjθ,h0进行离散卷积：pθ,h0=pθ?kjθ,h1;    (2) 对结果反投影到待建图像中原来按照抽样矩阵D抽取的点上，得到小波系数：Ajf［n］，  Djf［n］；    (3) 二维小波重构：Aj－1f=Ajf?0+Djf?1.    这里要注意的是：由于滤波器组不是正交的，所以分解滤波器和重构滤波器不同的（见式 (3)）。

    2  实验仿真

    我们把Coiflet2小波用于256×256像素 Shepp?Logan头部模型的全局图像重建。对图像进行了［0,π］之间的平行投影：每个角度256条平行射线，每次旋转角度为π/180，重建的图像大小为256×256像素。重建图像如图2所示。

    对 Shepp?Logan头部模型的局部重建中，重建区域为头部模型中心半径为64个像素的图像区域，细节图像和重建图像分别如图 4和图 5所示。局部原图像和FBP算法重建的局部图像如图 6所示。

    4  算法分析

    我们在图 1的近似图像和细节图像旋转了角度，这是由于抽取的方式不同造成的，梅花五点抽取方式是按照对角线方向抽取的（如图6所示）。梅花五点抽取方式使得每层近似图像和细节图像的抽样点是上一层图像的一半。每一层近似图像和细节图像的每一维的象素点是相等的，而且等于上一层图像的1/2。基于不可分离MRA的小波重建算法的图1  一层近似图像和一层细节图像全局重构图像如图 2所示。我们进行了均方误差分析，MSE的值为0.091，比基于可分离MRA理论的小波重建方法［1］重建图像的MSE小约10％，和FBP算法的MSE基本相同。

    这是因为在滤波过程中，可分离MRA算法的卷积是先在一维下进行的，所以图像被分离为行列无关联的信号，对于不可分离MRA重建算法，卷积使用的是二维滤波，这样图像被看作是一体的，而不是行和列分离的，所以图像在重建过程中信息的损失量要小。

    另外，因为分离的小波分解的是水平和垂直方向的（如图6所示），而不可分离的小波分解的在对角线上的，在梅花五点抽样矩阵的基础上建立的不可分离MRA的重建算法，损失的图像信息不在对角线上，所以更符合人类的视觉。

    在基于不可分离MRA理论的小波多尺度重建中，对于反投影这一步的量来说是和滤波反投影法也一样的。这是因为待建图像的像素大小并没有增加，只是用于重建的两套滤波后的投影数据之间互不相同，它们分别被用来重建图像的大致轮廓和细节特征。所以m×m的图像n角度投影时，对于反投影这一步的复杂度还是O(nm2)，而且这个复杂度与重建过程中小波分解的层数无关。而对于重建过程中的滤波这一步来说，它的复杂度是和小波分解的层数是相关联的。在滤波反投影法中，滤波这一步的复杂度为O(nmlogm)，因此在多尺度重建算法中对于小波分解的每一层，滤波这一步的复杂度都是O(nmlogm)。由于最大的分解层数可以达到log22m，这就使得总的滤波的复杂度可上升到O(nmlog2m)。于是和滤波反投影法相比，基于不可分离MRA的小波多尺度重建算法会花费更多的时间用于滤波，但与基于可分离MRA的小波多尺度重建算法用于滤波的时间减少一半。把反投影和滤波结合起来之后，基于不可分离MRA的小波多尺度重建算法总的复杂度为O(nm2)。

    对于局部成像，假定我们需要重建仍旧是位于图像中心半径为ri的局部区域。在不可分离MRA重建算法中，重构滤波器的支撑区间略大于可分离MRA重建算法中的支撑区间，而且滤波器经过希尔波特变换之后也不是绝对的紧支集函数，所以需要将X射线的照射区域的半径为re也略大于可分离MRA小波重建算法，但总的算法复杂度基本相同，而且使用在滤波的时间上是可分离MRA小波重建算法的一半。而总的算法复杂度为O(nr2e)，是取决于反投影的部分。全局重建的复杂度为O(nm2)，由于re比m小了很多，所以局部重建的复杂度 O(rrr2e)小于全局重建的复杂度，进而在算法实际运行时所消耗的时间大为减少。

【】
  ［1］Farrokh Rashid?Farrokhi， K J． Ray Liu， Carlos A. Berenstein. David Walnut. Waveletbased multi?resolution local tomography ［J］. IEEE Trans Image Proces, 1997, 6(10): 1412-1430.

　　［2］Jelena Kovacevic, Martin Vetterli, Nonseparable Multidimensional Perfect Reconstruction Filter Banks and Wavelet Bases for Rn ［J］. IEEE Transactions On Information Theory, 1992, 38(3): 533-555.

　　［3］ Kak A C and Malcolm Slaney, Principles of Com?puterized Tomographic Imaging, Society of Industrial and Applied Mathematics, 2001.

　　［4］陈武凡. 小波分析及其在图像处理中的应用 ［M］. 出版社． 2002: 113-118.

　　［5］汪元美. 医学成像理论 ［M］. 浙江大学出版社, 2000, 10: 37-49.