图像获取过程中不可避免地会收到噪声干扰，有效地去除噪声对提高图像视觉质量，提供精确的细节信息有着重要的作用^[1-2]。图像去噪算法根据处理域的不同可以将它们大致分为两类：时域去噪算法和频域去噪算法。时域去噪最经典的算法如Wiener滤波^[3]、双边带滤波(BF)^[4-5]、非局部均值滤波(NL-means)^[6]等，它们在时域定义滤波核，这类滤波器能较好的保持图像边缘，但很难保持图像的低对比度细节。频域去噪包括小波变换^[7]、短时傅里叶变换(STFT) ^[8]等方法，这一类方法对于保存图像纹理、细节有优势，但是对于图像尖锐边缘,由于Gibbs现象，处理效果欠佳。双域滤波算法将前面两种单域方法进行了结合，在抑制噪声的同时有效保护纹理细节。如Knaus ^[9]提出了双域去噪(DDID)，该算法通过联合双边带滤波将频域和时域滤波的结果进行互补获得了很好的去噪效果。张聚等人^[10]提出利用小波阈值法对小波域内的高频信号分量进行去噪，对小波域内的低频信号分量进行双边滤波处理，然后利用小波逆变换便得到去噪后的图像。

本文的主要研究对象为多模图像去噪，多模图像是指对同一目标在不同条件(不同时间、不同传感器等)下获得的两幅或者两幅以上的图像，如，对同一人体器官通过计算机断层成像(CT)和核磁共振成像(MRI)所产生的图像；对同一场景进行拍摄时带闪光灯和不带闪光的图像；用不同波段的成像设备获取的遥感图像等等。多模图像一方面可以将不同空间结构信息形成信息互补，另一方面这些图像之间又存在着空间位置上的高度相关性可以方便于对成像内容的结构进行深入分析。

目前所提出的多模联合去噪算法，如联合双边带滤波(JBF)^[11],交叉双边带滤波^[12]等，能够在时域中对多种模态的图像信息进行互补，滤波结果能够对图像中高对比度边缘进行很好的保持，但是对图像中细节纹理往往平滑过度。

针对以上提到的多模图像去噪所存在的问题，提出一种交叉双域滤波算法(CDDID)，该算法能够将时域多模联合滤波的优势与变换域小波收缩滤波的优势进行结合，并且通过时域和频域滤波交替迭代有效抑制振铃效应。相比于单模图像去噪和时域联合去噪，本文算法在保持图像高对比度边缘的同时对图像细节信息也能够进行较好地保留。

双边带滤波由两个函数构成，一个由几何空间距离决定滤波器系数，另一个由像素差值决定滤波器系数。与其他滤波方法相比，双边带滤波主要优点是容易计算，并且对图像尖锐边缘具有很好地保持，但对细节的处理效果一般。

对于一幅含噪图像y=x+η(σ²=Var[η])双边带滤波给出加权平均作为对原始图像x的估计，即

${{{\bar{y}}}_{p}}=\sum\limits_{q\in {{N}_{p}}}{{{k}_{p,q}}{{y}_{p}}}/\sum\limits_{q\in {{N}_{p}}}{{{k}_{p,q}}}$

(1)

式(1)是在特定邻域内的归一化加权平均，该邻域大小为N_p((2r+1)×(2r+1))。对于中心像素p，加权系数k_p,q由k^s_p,q和k^r_p,q两个加权因子相乘得到，即

${{k}_{p,q}}=k_{p,q}^{s}k_{p,q}^{r}$

(2)

$k_{p,q}^{s}=\exp \left\{ -{{\left| p-q \right|}^{2}}/\left( 2\sigma _{s}^{2} \right) \right\}$

(3)

$k_{p,q}^{s}=\exp \left\{ -{{\left( {{y}_{p}}-{{y}_{q}} \right)}^{2}}/\left( {{\gamma }_{r}}{{\sigma }^{2}} \right) \right\}$

(4)

式中，k^s_p,q 表示像素在空间距离决定的权重，k^r_p,q 表示像素灰度差决定的权重，σ_s和γ_r分别影响空间距离核以及灰度差核的形状。

Petschnigg等人^[11]在2004年提出的联合双边带滤波(JBF)算法是对双边带滤波^[4]算法的扩展。该算法适合于对多模图像进行去噪，在对一幅受噪声污染的图像y进行滤波时所使用的滤波权重由另一幅其他模态下较清晰的图像g来引导计算。

${{{\bar{y}}}_{p}}=\sum\limits_{q\in {{N}_{p}}}{{{k}_{p,q}}{{y}_{p}}}/\sum\limits_{q\in {{N}_{p}}}{{{k}_{p,q}}}$

(5)

核函数为

${{k}_{p,q}}=\exp \left\{ -\frac{{{\left| p-q \right|}^{2}}}{2\sigma _{s}^{2}}-\frac{{{\left( {{g}_{p}}-{{g}_{q}} \right)}^{2}}}{{{\gamma }_{r}}{{\sigma }^{2}}} \right\}$

(6)

由式(5)(6)可见，滤波核函数的灰度加权因子由引导图像g决定。该方法在对y进行去噪时由于使用了g在边缘处的信息，因此可以保留更多的图像细节。

在多模图像处理领域，由于成像设备的限制不同模态的图像可能都被噪声污染。有些情况下在对一幅图像去噪时并没有另一幅多模清晰图像来引导计算，本文提出的交叉双域滤波算能同时对两幅含噪多模图像进行滤波。该算法是一个迭代的过程，在每步迭代中，首先在时域通过交叉双边带滤波综合两种模态在边界的信息，得到保留了高对比度边缘信息的去噪图像。然后在各自模态中对所滤除的残差在频域中使用小波收缩算法获取各模态对于低对比度残量信息的估计，并且将该估计结果叠加到时域去噪结果上，由此得到一步迭代后的去噪图像。算法在时域和频域进行交替迭代，在迭代中通过滤波参数的不断调整，在逐渐获得图像更多的细节纹理的同时，对频域滤波所产生的振铃效应进行抑制，最终得到有效的滤波结果。

将原始的两幅多模噪声图像分别表示为y_p¹和y_p²(简写为y_p^α，(α=1，2))，两个图像分别对应的引导图像假设为g_p¹、g_p² (简写为g_p^α，(α=1，2))。在首次迭代时g_p¹、g_p² 分别赋值为原始噪声图像y_p¹和y_p²。在进行迭代时，分别取值为y_p^α上一步去噪的结果，此时，认为g_p^α比y_p^α具有更小的噪声。

交叉双边滤波算法对各个模态的图像在时域进行滤波时采用统一的灰度系数k_p,q ，该系数由像素在空间高斯权重k_p,q^s以及各个模态在灰度差的高斯权重k_p,q^r1、k_p,q^r2共同决定。即有

${{k}_{p,q}}=k_{p,q}^{s}k_{p,q}^{r1}k_{p,q}^{r2}$

(7)

$k_{p,q}^{s}=\exp \left\{ -{{\left| p-q \right|}^{2}}/\left( 2\sigma _{s}^{2} \right) \right\}$

(8)

$k_{p,q}^{r1}=\exp \left\{ -{{\left( g_{p}^{1}-g_{q}^{1} \right)}^{2}}/\left( {{\gamma }_{r}}\sigma _{1}^{2} \right) \right\}$

(9)

$k_{p,q}^{r2}=\exp \left\{ -{{\left( g_{p}^{2}-g_{q}^{2} \right)}^{2}}/\left( {{\gamma }_{r}}\sigma _{2}^{2} \right) \right\}$

(10)

对g_p^α、y_p^α均用式(7)中的k_p,q作为双边带滤波的核函数分别进行滤波处理，得到滤波后的图像 $\bar{g}_{p}^{a}$ 和$\bar{y}_{p}^{a}$，即

$\bar{y}_{p}^{\alpha }=\sum\limits_{q\in {{N}_{p}}}{{{k}_{p,q}}y_{q}^{\alpha }}/\sum\limits_{q\in {{N}_{p}}}{{{k}_{p,q}},\alpha =1,2}$

(11)

$\bar{g}_{p}^{\alpha }=\sum\limits_{q\in {{N}_{p}}}{{{k}_{p,q}}g_{q}^{\alpha }}/\sum\limits_{q\in {{N}_{p}}}{{{k}_{p,q}},\alpha =1,2}$

(12)

在该步骤中，滤波系数k_p,q是由两种模态共同决定的。因此在对边缘滤波时，两种模态的信息可以形成互补，只要有一种模态在该边缘上有较高的对比度，就能够有效调整滤波系数k_p,q减少对边缘的平滑。

类似双域滤波算法，可以将滤波前图像和滤波后图像的差(y_p^α－ $\bar{y}_{p}^{\alpha }$)(α=1，2)，作为图像在时域滤波后所丢失的细节信息，并称之为低对比度图像。算法通过对该低对比度图像进行小波收缩获取去噪后图像本身的细节信息然后叠加回高对比度图像，获得一步迭代的估计结果。由此在频域小波收缩前对低对比度图像做STFT计算出频域滤波系数为小波收缩作准备。STFT在离散傅里叶变换(DFT)前对信号加窗减轻边界振铃效应，窗口系数为式(7)中k^s_p,q。同时使用式(7)中灰度高斯权重k^r1_p,q和k^r2_p,q来作为DFT的权重，进一步减小边界对于频域去噪的影响。则整个过程等效为对低对比度图像作Gabor变换^[7]。所得系数G_p,f和Y_p,f的自变量为频率f，定义域为频率窗口F_p，该窗大小与N_p相同。综合以上，STFT系数及其方差为

$\sigma _{p,f}^{2}={{\sigma }_{1}}{{\sigma }_{2}}\sum\limits_{q\in {{N}_{p}}}{k_{p,q}^{2}}$

(13)

$Y_{p,f}^{\alpha }=\sum\limits_{q\in {{N}_{p}}}{{{e}^{-i2\pi \left( q-p \right)f/\left( 2r+1 \right)}}{{k}_{p,q}}\left( y_{q}^{\alpha }-\bar{y}_{p}^{\alpha } \right)}$

(14)

$G_{p,f}^{\alpha }=\sum\limits_{q\in {{N}_{p}}}{{{e}^{-i2\pi \left( q-p \right)f/\left( 2r+1 \right)}}{{k}_{p,q}}\left( g_{q}^{\alpha }-\bar{g}_{p}^{\alpha } \right)}$

(15)

式中，α=1,2。

最后一步，采用小波收缩因子K_p,f^α，α=1，2 来收缩含噪STFT系数Y_p,f^α，α=1，2。收缩因子K_p,f^α的构造利用了导引图像谱G_p,f^α，并采用欧氏距离倒数σ_p,f²/|G_p,f^α|²构造，γ_f为小波收缩因子调整参数，K_p,f^α 的效应是对高频噪声系数进行抑制，即

$K_{p,f}^{\alpha }=\exp \left\{ -{{\gamma }_{f}}\sigma _{p,f}^{2}/{{\left| G_{p,f}^{\alpha } \right|}^{2}} \right\},\alpha =1,2$

(16)

在频率窗口F_p内，由DFT零频率系数与直流分量之间的关系，可得平均所有收缩后小波系数即可得到当前点的逆DFT，产生低对比度信号去噪后的结果

$\bar{Y}_{p}^{\alpha }=\frac{1}{\left| {{F}_{p}} \right|}\sum\limits_{f\in {{F}_{p}}}{K_{p,f}^{\alpha }Y_{p,f}^{\alpha },\alpha =1,2}$

(17)

将交叉双边带滤波得到的高对比度图像与小波收缩得到的低对比度图像相加，得到一次对噪声图像的估计，即

${{y}^{\alpha }}\approx {{{\bar{y}}}^{\alpha }}+{{{\bar{Y}}}^{\alpha }},\alpha =1,2$

(18)

经过多次迭代后得到对噪声图像的最终估计。

多模图像交叉双域滤波算法步骤概括如下：

1)输入多模图像y_p^α,α=1,2。

2)g_p^α 赋值为上次迭代结果，初始值赋值为y_p^α。

3)通过式(7)计算k_p,q，通过式(11)(12)得到时域交叉滤波结果 $\bar{g}_{p}^{\alpha }$、$\bar{y}_{p}^{\alpha }$。

4)通过式(13)—(15)计算小波收缩系数。

5)通过式(16)(17)对各个模态在频域进行小波收缩。

6)通过式(18)得到各个模态的去噪估计值 ${{{\bar{y}}}^{\alpha }}+{{{\bar{Y}}}^{\alpha }}$。

7)若没有达到迭代次数，则更新g^α为g^α≈ ${{{\bar{y}}}^{\alpha }}+{{{\bar{Y}}}^{\alpha }}$ ，更新参数γ_r、γ_f ，返回到步骤2)；

否则，输出最终各模态的滤波结果 ${{{\bar{y}}}^{\alpha }}+{{{\bar{Y}}}^{\alpha }}$。

算法在步骤6)将频域处理后的结果叠加到时域中，实际上也引入了振铃效应。但是由于算法在下一次迭代时，在时域中再次使用了交叉双边带滤波，因此可以将振铃效应当做噪声进行抑制，达到了将两种滤波优势相结合的目的。从式(9)(10)以及式(16)可以看到，每次滤波的系数g_r、γ_f主要控制时域以及频域滤波的范围，随着参数的不断减小可以保留图像中更多的细节信息。但是前后两次迭代的衔接也很重要，后一次迭代中的时域滤波参数g_r过小，可能达不到对振铃效果的抑制，直接导致滤波效果较差。下面在实验部分，对参数的选取进行了进一步的讨论。

本文所使用噪声方差和算法中的参数的设置，均假设图像在处理前已经归一至[0，1]。在所有实验中，根据经验设置交叉双边带滤波的滤波半大小为r=15，设置空间滤波核大小的σ_s均设为7。根据多模图像的特点，首先设计了一组合成图像实验，来验证DDID^[9]与CDDID的去噪以及边缘保持能力。在图 1(a)(b)为原始清晰的合成图像，这两幅图片的背景灰度值均为 127，白块灰度值为 200，黑块为 50，另外两个灰块均为 135，并且两幅图在同一位置处正好是强弱边缘相对应。图 1(c)(d)分别为对图 1(a)(b)添加噪声后的图像，且噪声强度的PSNR=20 dB(σ=0.1)。使用双域滤波算法与交叉双域滤波算法的去噪结果分别如图 1(e)—(h)所示，所示结果均选取算法第3次迭代后的结果。算法中式(9)(10)中的σ₁、σ₂分别取值为和噪声强度方差相当的0.1，3次迭代参数根据经验分别设为g_r=(20,3,1),γ_f=(10,6,2)。从结果可以观察到,对图像加噪之后图像中的灰色小方块与背景色之间这样的弱边缘已经很难区分。双域滤波方法尽管能比较好地去除噪声,但是弱边缘依然比较模糊。而交叉双域滤波算法在边缘保持上体现出明显的优势，算法利用了多模图像之间空间位置的相关性与信息互补的优势,在去噪时将两幅多模图像在时域进行了交叉联合处理,合成图像1(合成图像2)中的弱边缘去噪的系数受到在合成图像2(合成图像1)中为强边缘信息的引导，最后尖锐边缘得到了较好的保持。

图 2为合成图像使用双域滤波算法与使用交叉双域滤波算法的3次迭代结果对比。从左到右分别为每次迭代后对图像中间一行进行采样获得的曲线。通过PSNR值可以看出，交叉双域滤波算法的去噪效果较双域滤波要好。图 2中的黑色曲线为原始无噪声图像，可以看到CDDID方法得到的绿色曲线与黑色曲线的重合度更高，尤其是在边缘位置具有更加尖锐的边缘。

图 3中的原始图像分别为从Brainweb数据库(http://www.bic.mni.mcgill.ca/brainweb)上获得的模拟脑部核磁共振的多模图像“T1”和“T2”。第2列为分别为对“T1”和“T2”添加高斯噪声PSNR=20 dB(σ=0.1)后的图像。本文将JBF、DDID，以及CDDID的结果进行了对比。可以发现所本文算法与时域中的联合双边带滤波方法以及单模态的双域滤波方法相比具有明显优势，本文算法更好地保持了图像的边缘信息以及纹理细节。

图 4为对Petschnigg等人^[11]所提供的闪光和非闪光测试图像进行滤波的结果。图 4第1行为闪光图像，由于原始图像噪声较小，实验中加入强度为PSNR=20 dB(σ=0.1)的高斯噪声。图 4第2行为非闪光图像，其在原始采集所得到的图像中已经包含较大噪声。由于两种模态的图像噪声强度不同，因此在实验中，按照各模态噪声的强度分别设置式(9)(10)中σ₁=0.1,σ₂=0.3。在3步迭代中，g_r的参数依次取为(100,8.7,0.7)，γ_f的参数依次取为(4,0.4,0.8)。可以看到，本文算法不但保留了图像高对比度信息，如罐子的高光区域，同时也较好地还原了图像中细节纹理，如罐子上的图案。并且在去噪后图像的PSNR值上相比其他方法有了较为明显的提高。

本文算法的参数目前需要根据经验进行设置。式(9)(10)中σ₁和σ₂在实际中一般设置为对各个模态噪声强度估计以后的方差值。在迭代中，g_r和γ_f值的设置，也对结果有一定的影响。图 5给出了不同参数组对多模图像处理结果的影响。图 5(a)为实验中准最优参数下的结果。由式(9)(10)可以看到，g_r代表对灰度值平滑的强度，如果设置较小对噪声的整体抑制较差，如图 5(b)的第1次迭代设置为g_r=1，导致滤波结果噪声整体较大。同时g_r的值在迭代时不能过小，否则对上次频域滤波的振铃效果抑制较差。如图 5(c)在第3次迭代时g_r=0.1，虽然滤波结果纹理更加清晰，但有较明显的振铃效应干扰。γ_f根据式(16)为对频域小波收缩的系数。若在迭代中，设置过小，将产生较为严重的振铃效应，在后期迭代中也很难消除，如图 5(d)在最后一次迭代中设置γ_f=0.2，振铃效应非常明显。

通过算法仿真，本文多模图像交叉双域滤波算法，能够有效利用多模图像结构互补的特点对滤波后图像的边缘进行较好的保持。另外，算法由于在时域与频域算法相结合，因此图像细节纹理也获得了较好地保留。通过与传统基于单域和单一模态的去噪算法相比，本文算法具有明显优势。

本文算法实现较为简单，并且可以推广到两个以上多模态图像去噪应用中。假设存在n(n>2)个模态的噪声图像，去噪时可将式(7)扩展为k_p,q= $ \prod\limits_{i=1}^{n}{k_{p,q}^{s}k_{p,q}^{ri}}$，其他步骤与本文算法类似。

虽然本文算法提供了一种能有效将时域联合去噪和频域去噪优势进行互补的算法，但是参数的选择同其他去噪算法类似，存在一定的主观因素。今后的研究将主要解决如何在某种合理量化评价体系下，估计全局或者局部最优的迭代参数。