Curvelet（共7篇）
Curvelet 篇1
1 引 言
由于小波变换在空域和频域上都具有良好的局域特性,近年来他在图像消噪中的运用越来越广泛,但是小波分析主要反映奇异点的位置和特性,而二维图像的边缘有许多曲线和直线,使得小波变换在处理图像时具有一定的局限性。为了克服这种局限性,EJ.Candes提出了Curvelet变换[3],Curvelet变换是一种具有方向性的多尺度变换,他能够有效描述沿直线的奇异特性,因此在对图像进行处理时能够比小波变换更好地保护图像中的线性特征。
2 Curvelet变换
Curvelet变换的分解和重建过程如图1所示:
分解过程首先是对函数进行子带分解,即将f用滤波器P0分解为(P0f,Δ1f,Δ2f,…),然后是对不同尺度的子带函数进行平滑分块,即定义平滑窗口ωQ,使其位于方形区域Q=[k1/2s,(k1+1)/2s)×[k2/2s,(k2+1)/2s),将窗口和每个子带函数相乘,实现平滑分块。接下来对平滑分割后的各个方块进行重正规化处理,即对每一个二进方块Q,定义(TQf)(x1,x2)=2sf(2sx1-k1,2sx2-k2),用该算子将位于方形区域Q的函数正规化到单位尺度。对f进行平滑分块和重正规化后,可以得到:
undefined
Curvelet变换的核心是Ridgelet变换[1],他的定义为:
undefined
可以看出,Ridgelet变换是沿着直线xcos θ+ysin θ=t的一维小波变换,他是一种具有方向性的多尺度变换,因此能更能够有效描述沿直线的奇异性。也就是说,他在处理图像线性区域时能够起到比小波变换更好的效果。
Ridgelet变换通过Radon变换和小波变换实现,对于平面(x,y)在R2上的函数f(x,y),他的Radon变换就是该函数在各个角度上直线的投影,即:
undefined
式(3)中δ是Delta-Dirac函数。而f(x,y)的Ridgelet变换系数可以通过对他的Radon变换系数进行小波变换来计算:
undefined
根据傅里叶定理:
undefined
其中F是f(x,y)的傅里叶变换。
因此,f(x,y)的Radon变换可以由f(x,y)的二维傅里叶变换在径向做一维傅里叶反变换得到,即首先对f(x,y)做二维傅里叶变换,再将直角坐标向极坐标转换,最后在极坐标方向实现一维傅里叶反变换,得到f(x,y)的Radon变换。
Curvelet重建是分解过程的逆变换,首先是进行Ridgelet反变换,再进行重正规化,然后是平滑集成,最后进行子带重建。
3 基于Curvelet变换的图像消噪
小波萎缩阈值消噪算法是现在运用最广泛的小波消噪方法,他主要是根据信号和噪声在小波变换后的不同特性进行消噪。在小波变换下,噪声的平均幅值与尺度因子2j成反比;平均模极大值个数与2j成反比。 即噪声的能量随着尺度的增加而迅速减小。而在小波变换下图像信号的平均幅值不会随着尺度的增加而明显减小;而且,噪声在不同尺度上的小波变换是高度不相关的。信号的小波变换一般具有很强的相关性,相邻尺度上的局部极大值几乎出现在相同的位置上,并且有相同的符号。根据这些特点,可以选择合适的阈值,将小于阈值的系数置零,保留大于阈值的小波系数,然后经过阈值函数映射得到估计系数,最后对估计系数进行反变换,就可以实现消噪和重建。
同小波变换一样,Curvelet变换系数也具有相同的特点,可以通过阈值化处理去除噪声。算法步骤如下:
(1) 对含噪声图像进行Curvelet变换,得到Curvelet变换系数;
(2) 对图像的Curvelet变换系数进行阈值操作,若系数大于阈值δ则保留,若小于δ则将其置零;
(3) 对处理后的Curvelet变换系数进行Curvelet反变换,得到消噪后的图像。
4 实验结果
这里选取512×512的Lena图像进行实验,将方差为σ=0.1和σ=0.078的高斯白噪声n加入到图像中。结果由式(5)给出:
undefined
图像的峰值性噪比由以下公式给出:
undefined
其中均方差:
undefined
实验中对2种噪声情况下的典型小波硬阈值消噪方法和Curvelet消噪方法进行比较(见表1),采用undefined作为阈值,发现无论是在峰值性噪比,还是对图像边缘的处理效果,Curvelet消噪方法都比典型的小波消噪法有一定的提高。
5 结 语
将Curvelet变换运用于图像消噪中,得到比典型的小波硬阈值消噪算法更好的效果。可以看出在Lena图像的帽檐和肩膀等线性特征较明显的部分,Curvelet消噪明显要强于小波消噪,而且经过Curvelet消噪后的图像, PSNR要高于经过小波消噪后的图像。
参考文献
[1]Candes E J.Ridgelets:Theory and Applications PhD Thesis,Stanford University,1998.
[2]Candes E J,Donoho D L.Curvelets-A Surprisingly EffectiveNonadaptive Representation for Objects with Edges.Curvesand Surfcaces,Vanderbilt University Press,Nashvielle TN,2000:105-120.
[3]Starck J L,Candes E,Donoho D.The Curvelet Transformfor Image Denoising[J].IEEE Trans.Image Processing,2002(11):670-684.
[4]倪林,Y.Miao.一种更适合图像处理的多尺度变换———cur-velet变换[J].计算机工程与应用,2004,40(28):21-26.
[5]李晖晖,郭雷,刘航.基于二代curvelet变换的图像融合研究[J].光学学报,2006(5):19-24.
Curvelet 篇2
近些年来，多传感器图像融合理论研究日渐深入，应用领域也日益广泛，在医学[1]、机器视觉、环境保护和遥感[2]等领域都取得了比较广泛的关注和应用。图像融合是指将配准后的图像采用一定的方法融合在一起的技术，由于待融合的多幅图像之间具有信息的冗余性和互补性，所以融合过程需要在降低冗余信息影响的同时尽量提取互补信息。
早期的图像融合尝试均是在空间域上实现的，运用的主要算法有加权平均法、主成分分析法、Brovery变换融合算法。随后，陈武、曹喆等人将小波变换应用到图像融合领域，提出基于小波系数方向的子带系数选择策略[3,4]。小波分析兼具多分辨率和时频局部化特性，成为图像融合的重要技术手段。多尺度分析工具[5,6]不断的发展，本文将真正意义上的二维信号分析工具Curvele引入图像融合，根据图像的特点，提出一种基于区域自适应特性的Curvelet域融合算法，对Curvelet分解后的低频和高频子带进行融合操作，在低、高频子带采用不同的融合规则。实验表明，本文算法优于传统的图像融合算法。
1 Curvelet变换
1999年，Candes和Donoho等人在Ridgelet变换的基础上提出了连续曲波（Curvelet一代）变换，然而其变换过程数学实现较为复杂。于是Candes等人在2002年提出第二代Curvelet变换，并于2005年提出两种基于第二代Curvelet变换理论的快速离散实现方法[7]。二代Curvelet较第一代变换更加简单，冗余度更低，运算更为迅速，实现也更加方便。与传统小波变换相比，Curvelet变换更适于刻画图像的几何特征，如曲线、直线等，采用“楔形基”来逼近C2（二阶连续可微）的奇异点，充分考虑奇异点的几何形状，并具有任意角度的方向性（各向异性），更适合对图像的处理与应用，此外，Curvelet变换对图像的几何特征（曲线、直线）具有更好的系数表达能力，采用少数较大的Curvelet变换系数进行表示，克服了小波变换传播重要特征到多个尺度上的缺陷，变换后能量更加集中，更利于跟踪和分析图像中的重要特征[8]。
Curvelet变换有两种不同的实现算法，这两种方法的主要区别在于不同尺度和方向下空间网格的选择方法不同。下面简单介绍两种方法[9]:
(1)USFFT算法
对f[t1,t2]∈L2(R2)做二维FFT变换得到Fourier采样序列F[n1,n2]，其中
针对频域中不同尺度和方向参数(j,l)，对F[n1,n2]进行插值得到F[n1,n2-n1tanθ1]，其中(n1,n2)∈Pj;
将F[n1,n2-n1tanθ1]与拟合窗口相乘之后得到如下结果，即
对序列Fj,l[n1,n2]做二维FFT逆变换，得到Curvelet系数CD(j,l,k)。
(2)Wrapping算法
对f[t1,t2]∈L2(R2)做二维FFT变换得到Fourier采样序列F[n1,n2]，其中
针对频域中不同尺度和方向参数(j,l)，用拟合窗乘以F[n1,n2-n1tanθ1];
围绕原点Wrap局部化F,
对序列Fj,l[n1,n2]做二维FFT逆变换，得到Curvelet系数CD(j,l,k)。
对比以上两种方法，都是先通过FFT变换到频域，再在频域中进行局部化，然后对局部化后的结果做二维FFT逆变换得到离散Curvelet变换系数。
2 基于区域特性的Curvelet图像融合算法
2.1 图像融合流程
基于Curvelet变换的图像融合的实现流程如图1所示。算法的实现步骤表述如下：
(1)Curvelet分解。将配准好的源图像A和B进行Curvelet变换，得到相应的Curvelet系数集合。分解尺度和方向为默认值。分解后得到的矩阵C为一个胞矩阵，其中C{j}{l}就表示是一个二维矩阵，代表尺度j、方向l上的所有系数，高尺度对应高频系数。
(2）图像融合。对于分解后的低频子带和所有高频子带，使用基于区域特性的融合规则进行判别和融合处理，得到各尺度上融合后的Curvelet系数。
(3)Curvelet重构。重构是分解的逆过程，对融合后的Curvelet系数进行Curvelet逆变换，得到重构的融合图像，该图像包含原有多幅图像中的信息。
2.2 融合规则研究
融合规则是图像融合算法的核心，其对最终的融合效果有着决定性作用。分解后的低频和高频子带具有不同的物理信息，低频子带表征图像的近似部分，而高频子带表征图像的细节部分。通常对于多分辨率融合，在高频和低频上采用不同的融合规则。总的来说，现有的融合规则可分为两大类：基于像素的融合规则和基于区域的融合规则。
2.2.1 融合规则分类
(1）基于像素的融合规则
Burt最早提出基于像素选取的融合规则，基于像素点绝对值最大来选取最终的融合像素值。Petrovic等人提出考虑分解层内各子带图像相关性的像素选取融合规则。近几年，融合规则层出不穷，陈武等人针对小波处理图像边缘的不足，根据人类视觉系统对局部对比度敏感的特性，采用了基于区域对比度的像素选取融合规则。基于像素的融合规则具有实现简单、融合速度快的优点，但其仅仅以单个像素作为融合对象，并未考虑图像像素间的相关性，因此融合效果较差，适用场合也非常有限。
(2）基于区域的融合规则
考虑到图像像素之间的相关性，Burt等人提出基于区域特性选择的加权平均融合规则，将系数值的融合选择与其所在的局部区域联系起来。Qingping Li等人则利用区域频率划分区域进行融合实验[10]。Curvelet提出后，玄立超等人提出了基于Curvelet变换的区域能量高频融合策略[11]。基于区域的融合规则就是将某位置邻域的能量、梯度、方差等特征作为一种测度来指导该位置处的系数选取，邻域大小可以是3×3、5×5等。由于基于窗口区域的融合规则考虑了图像相邻像素间的相关性，因此减少了融合像素的错误选取，从而提高了融合效果。
结合现有融合规则的优点和不足，本文提出以下融合规则对高、低频分别处理。对低频采用基于区域方差显著性的加权融合规则，充分考虑到邻域像素间的关联性，并有效保留细节和边缘。高频采用区域能量匹配融合规则，抑制噪声的同时能更好反映图像区域特征。
2.2.2 低频子带融合规则
低频子带是包含了原图像的主体信息，决定了图像的大致轮廓。这里采取基于区域方差显著性的加权融合规则对低频系数进行处理，并将匹配度与自适应阈值比较来确定融合策略。实现过程如下：
(1）计算区域方差显著性
式中：G(X,p)为以点p为中心的局部区域Q的区域方差显著性；w(q)为权值，权值通过行和列的高斯分布加权相加得到；C(X,q)为区域Q形成的矩阵；为图像X的低频系数矩阵内区域Q的平均值。
计算区域方差匹配度：
式中：M(p)表示图像A和B低频系数在p点的局部区域方差匹配度，数值在0～1之间变化，其值越小表明两幅图像低频系数在p点的相关程度越低。
融合策略的确定。设T为匹配度阈值，不同于传统的经验选取[8,9]，此处T为自适应阈值，其定义为：
式中：G(A,p)和G(B,p)分别为源图像A,B区域方差显著性。当M(p)<T时，采用选择融合策略：
当时，采用加权平均融合策略：
其中，
式中：C(A,p),C(B,p)分别为图像A,B的低频系数在p点的值。
该策略是基于邻域像素间的关联性，并且是基于区域方差，可以有效地保留细节和边缘，因此采用该策略得到的融合图像将比较清晰，细节比较丰富。
2.2.3 高频子带融合规则
不同尺度及方向下的高频系数矩阵，其中包含了图像的细节特征，由于人眼对于单个像素的灰度取值并不敏感，图像清晰度是由区域内像素共同体现的。因此传统的基于绝对值取大的高频融合规则是比较片面，忽略了像素之间的关联性，并且忽略了图像之间的相关性。所以这里引入基于局部能量匹配的融合规则，可以有效抑制噪声，并且能更好地反映图像的区域特征。实现过程如下：
(1）计算区域能量显著性
式中，Ej,l(m,n)表示尺度j和方向l下高频系数的局部区域能量；Q为以点P(m,n)为中心选取的局部区域。
(2）计算区域能量匹配度。Mj,l(m,n)的取值为0～1，取值越小表明相关程度越低。
(3）融合策略的选取
设T1为能量匹配度阈值，一般取0.5～1，这里根据经验取0.85。当Mj,l(k1,k2)<T时，采用选择融合策略：
当时，采用加权融合策略：
其中，
该策略基于区域空间能量，可以有效保留图像的边缘信息，能得到细节比较出众的融合图像。
3 图像融合效果评价指标
对于不同融合目的的图像，需要采用不同的评价标准，从而对融合结果进行客观正确的评价[12]。例如，提高信息量时，对于融合图像的信息量是否增加，可以根据互信息、交叉熵、熵和标准差等指标来评价；提高图像清晰度时，往往是要求在保持原有信息不丢失的情况下，增强图像的细节信息和纹理特征，评价这种目的的融合时可选用平均梯度、空间频率等指标来评价。需要注意的是，客观评价法离不开主观评价，因此应该将两者结合起来进行综合评价。主观评价主要是由人眼来观察区分结果的好坏，但人眼分辨力有限，所以需要引入客观评价指标。用来评价融合结果的客观指标[13]如下：
(1）信息熵
信息熵可以从概率分布的角度来衡量图像的丰富程度，图像信息熵越大，表明图像所包含的信息量越大。熵值的大小可以反映图像对细节的表达能力。其定义如下：
式中：H代表信息熵；p(l)表示图像中像素灰度级为l的出现概率，即所有灰度为l的像素点数N1与图像中所有像素点数N之比。
(2）交叉熵
交叉熵用来反映两幅图像的信息差异。通过对融合图像和源图像交叉熵的计算，就可以得到两幅图像所包含信息量的差异。一般来说，交叉熵越小，表明融合图像从源图像中提取的信息越多，信息差异越小。其定义如下：
式中PA(l)和PF(l)分别为源图像A和融合图像F的灰度级概率分布。通常把两幅源图像分别和融合图像的交叉熵的的均方根定义为均方根交叉熵RCE。即：
通过RCE可以表示融合图像F和图像A、B之间的联合差异。一般情况下均方根交叉熵越小，融合图像继承的信息量越多，信息差异越小，融合效果也就越好。
(4）相关系数
图像的相关系数反应了两幅图像的相关度。相关系数越大，其融合效果越好。
其定义如下：
式中：Iai,j,Ibi,j为两幅图像在(i,j)点的灰度值；eIa和eIb为两幅图像的均值。把两幅源图像分别和融合图像的相关系数的均方根定义为均方根相关系数RC，即：
一般来说，均方根相关系数越大，融合图像从源图像中获得的信息量越多，融合效果也就越好。
(5）空间频率
空间频率可以反映出融合图像的细节表达能力，所以通常作为度量图像清晰度的指标。空间频率定义为：
式中：RF和CF分别为行频和列频，f(i,j)为(i,j)处的灰度值。
一般来讲，空间频率越大，图像的层次越多，融合图像就越清晰。
4 实验结果对比分析
选用两组不同领域图像作为源图像进行融合实验，第一组是医学图像，大小为256×256；第二组是多聚焦图像，大小为256×256。第三组是多波段图像，大小为512×512。实验中高、低频子块大小定义为3*3，分别采用基于低频加权平均、高频绝对值最大的离散小波变换（方法1）、基于低频加权平均和高频区域频率的离散小波变换（方法2）、基于文献[11]的Curvelet变换（方法3）、基于本文算法实现。
4.1 医学图像融合结果分析
医学图像融合结果分析如图2，表1所示。
图2和表1分析如下：
(1）基于小波变换的融合算法，基于二代Curvelet变换的融合算法的比较：
主观上看，基于小波变换的融合图像（c)(d）明显在边缘的处理上不够清晰，会发现有方块效应，使得融合图像没有任何意义；而基于二代Curvelet变换的融合图像（e)(f）就很好的保留了图像的边缘信息，融合图像更加自然。
(2）传统的融合算法和本文融合方法的比较：
表1中，医学图像融合是为了获得更多信息，将图3(f）与（c）、（d）、（e）对比来看，（f）均方根交叉熵更小，即本文融合算法融合结果从源图像继承了更多信息。整体来看，本文算法融合结果也比其他融合结果更加清晰。综上所述，医学图像的融合实验中，基于Curvelet变换的本文算法优于传统Curvelet、小波变换的融合算法。
4.2 多聚焦图像融合结果分析
多聚焦图像融合结果分析如图3，表2所示。
图3和表2分析如下：
(1）基于小波变换的融合算法，基于二代Curvelet变换的融合算法的比较：
主观来看，基于小波变换的融合图像（c）、（d）在小钟表的边缘有虚影，使得融合图像比较模糊，边缘处理不够细致。而基于Curvelet变换的图像就没有这一缺点。
(2）融合方法之间的对比：
表2中，基于Curvelet变换的融合图像包含有更丰富的信息，并且与源图像相关性更高。对比交叉熵和相关系数可以看出，本文融合算法比文献[11]融合算法从源图像中继承更多的信息的同时保留了更高的相关性。
综上所述，多聚焦图像融合实验中，本文对于高、低频的融合规则比传统融合规则可获得更好的融合效果。
4.3 多波段遥感图像融合结果分析
多波段遥感图像融合结果分析如图4，表3所示。
图4和表3分析如下：
(1）基于小波变换的融合算法，基于二代Curvelet变换的融合算法的比较：主观来看，基于Curvelet变化的融合图像（f)(g）相对基于小波变换的融合图像（d）、（f）来说，图像对比度更大，边缘细节更清晰。
(2）融合方法之间的比较：表3中，从信息熵、交叉熵、相关性来看，本文融合方法都优于传统融合算法，从源图像中继承了更多的信息，但图像频率略有降低。
综上所述，多波段遥感图像的融合实验中，基于Curvelet变换的融合算法优于小波变换融合算法，本文算法中的高、低频系数基于区域系数相关性法也比传统的模值取大取平均方法更适用于多波段图像的融合。
总体来说，基于二代Curvelet变换的本文融合算法在不同领域的图像融合中效果比较理想。
5 结语
本文分析了传统小波变换在边缘处理上的不足，使用了多分辨率分析中具有多尺度多方向性的二代Curvelet变换，同时提出针对高、低频系数的融合规则，并对不同领域的图像做了大量的融合实验。实验结果表明，本文融合算法能较好的保持图像目标信息，同时边缘也比较清晰。但是，该算法也存在不足，Curvelet变化在图像的细节特征方面有天然的弱势；故结合多种小波[1]、基于区域分割[14]或借助其他技术[15]使得融合图像在保留目标信息的同时具有更好的细节特征和边缘轮廓是未来图像融合的研究重点。
摘要：为克服小波变换在二维或更高维度空间分析中的缺陷,提高图像融合质量,提出基于二代Curvelet变换的图像融合改进算法。引入可以有效分析图像中的曲线奇异性,能更加合理处理图像边缘信息的Curvelet变换对图像进行分解,对图像分解后的低频部分采用自适应阈值的区域方差高斯加权融合方法,增加图像像素之间的关联,并有效保留细节和边缘。对高频部分采用区域能量融合方法来降低噪声,增强图像的细节。采用该算法对多组不同图像进行融合实验,并用信息熵、交叉熵、相关系数、空间频率等对融合图像进行客观评价。实验结果表明,该算法优于传统的融合规则和算法,能在保持更好清晰度的同时获得更丰富的图像内容。
Curvelet 篇3
遥感影像融合是指把来自不同传感器或同一传感器不同尺度的同一场景的两幅或多幅影像, 采用一定的算法生成一组新的信息或合成影像, 以提高影像的清晰度和可识别性, 获得单一影像所不能提供的特征信息[1]。遥感影像融合在航空、航天、军事侦察、灾害预报等很多军事及民用领域有着举足轻重的地位。它们为资源调查、环境监测等提供了丰富而又宝贵的资料[2], 从而构成了用于全球变化研究、环境监测、资源调查和灾害防治等多方面应用。
根据信息抽象程度以及融合应用层次的不同, 遥感信息融合可划分为像素级、特征级和决策级三个层次。像素级融合是遥感信息融合各层次中最成熟的一级, 像素级融合能保持尽可能多的原始数据, 可提供其他融合层次所不能提供的细微信息, 具有最高的精度。目前光谱域主要的融合算法有代数运算法[3]、Brovey变换[4]、IHS变换法[5]、主成分变换法[3]、基于平滑滤波的融合[6]、以及现在热门的小波分析法[7]等。
小波变换是目前广泛使用的非自适应方法, 二维可分小波是一维小波的扩展而成, 各向同性的性质导致方向选择性差, 不能有效地捕捉轮廓信息, 是一种“非稀疏”图像表示法, 因此有必要寻求比小波更有效的方法。为了解决这个问题, 一些新的高维函数的最优表示方法Ridgelet、Curvelet等应运而生, 这些变换所采用的基的支撑区间表现出更高的方向敏感性, 即具有“各向异性 (Anisotropy) ”, 能获得图像更“稀疏”表示, 更好地反映图像的细节信息, 因此它们比小波能更好地表现边缘特征, 更适合于进行图像融合。本文介绍了Curvelet变换以及利用Curvelet变换进行自适应遥感图像融合, 最后对各种融合方法进行分析比较。
1 Curvelet变换
1.1 Curvelet变换
Candes和Donoho于1999年提出了Curvelet变换, 即一代Curvelet变换, 它是在Ridgelet变换的基础上发展而来, 与单尺度的Ridgelet相比, Curvelet是在所有可能的尺度s≥0上进行的Ridgelet变换。可以把一代Curvelet变换看成一种多尺度的Ridgelet变换。一代Curvelet变换是在连续域中定义的, 它的数字实现比较复杂, 需要经过子带分解、平滑分块、正规化和Ridgelet分析等一系列步骤, 而且Curvelet金字塔的分解也带来了巨大的数据冗余量, 因此E.J.Candes等又提出了实现更简单、更便于理解的快速Curvelet变换算法, 即二代Curvelet (Fast Curvelet transform) 。
二代Curvelet与一代Curvelet在构造上已经完全不同, 一代Curvelet的构造思想是通过足够小的分块将曲线近似到每个分块中的直线来看待, 然后利用局部的Ridgelet分析其特征。而二代Curvelet和Ridgelet理论并没有关系, 实现过程也无需用到Ridgelet, 二者之间的相同点仅在于紧支撑、框架等抽象的数学意义[8,9,10]。Curvelet变换有连续的Curvelet变换和离散的Curvelet变换, 它们对频域窗口的表示是不同的, 如图1所示, 其中常用的是离散的Curvelet变换。
1.2 实现方法及系数分析
在文献[11]中, 作者提供了两种快速离散的Curvelet变换方法:基于USFFT (Unequispased FFT) 变换法和基于Wrapping变换法。由于USFFT算法采用装配技术, 因此计算量比Wrapping方法大, 速度不如Wrapping方法[12,13], 因此, 本文中采用的是基于Wrapping的变换方法, 该方法的实现过程是:
(1) 对于在笛卡尔坐标下的二维函数
(2) 对于每一尺度j和每一个角度1, 得到乘积
(3) 围绕原点Wrap局部化
(4) 对局部化的
经过Curvelet变换后得到C{j}{1} (k1, k2) 结构的系数中, j表示尺度, 1表示方向, (k1, k2) 表示j层上第1个方向的矩阵坐标。512×512的图像的系数格式如表1所示。
从表1中可以看出, 大小为512×512的图像经过五层Curvelet分解, 被划分为六个尺度层, 最内层是Coarse尺度层 (C{1}) , 即图像的低频系数, 大小为21×21;最外层称为Finest尺度层 (C{6}) , 为高频系数, 大小和原始图像一样;中间层称为Detail尺度层 (C{2}-C{5}) , 为中、高频系数。其中每层被分解为四个大方向, 每个大方向被分解为N个小方向, 本例中, N=4, 8, 8, 16。Detail尺度层中各个方向子带的大小为m×n, 各方向子带的数目为k。
2 基于二代Curvelet变换的自适应遥感图像融合算法
2.1 融合规则的选取
在对多光谱图像和全色图像进行二代Curvelet分解后, 分别得到不同尺度的系数, 其中, 我们对低频系数子图像和高频系数子图像选用不同的融合规则进行融合。
(1) 低频信息融合规则
低频信息融合规则采用基于区域能量的加权系数的自适应选取。首先求出以 (m, n) 为中心位置的局部能量EJA (m, n) 和EJB (m, n) , 后按照式 (4) 可计算出融合图像低频系数的选择。
C (m, n) =EJA (m, n) EJA (m, n) +EJB (m, n) CJA (m, n) +
EJB (m, n) EJA (m, n) +EJB (m, n) CJB (m, n) (4)
按照式 (4) , 随着 (m, n) 中心位置的变化, 区域能量也随着变化, 加权系数也随着变化, 其中加权系数可以根据区域能量自适应的变化。体现了自适应的特征。
由于局部区域能量较大的中心像素代表了原始图像的明显特征, 由式 (4) 若求得图像A的某一区域能量较大, 则对应的加权系数也会较大, 若区域能量较小, 则对应的加权系数也会较小, 符合原始图像本身特征对加权系数影响的特点。所以这种自适应的融合规则, 是有效的、可行的。
(2) 高频信息融合规则
Curvelet变换中的高频细节信息中包含着丰富的图像特征细节信息。其中高频信息中绝对值较大的系数对应着一些突变, 如图像的边缘、纹理等重要特征信息。因此高频信息的融合规则选择基于区域特征的算法有利于提取原始图像中的特征信息。高频分量选择以空间频率为度量标准结合自适应因子进行基于区域的自适应融合。
分别计算两幅图像Curvelet变换分解系数中对应尺度下相同方向的以 (i, j) 为中心的局部区域的空间频率, 得到对应尺度下相同方向的融合后高频系数, 即:
式中D
2.2 基于Curvelet变换的自适应遥感图像融合方法
Curvelet变换具有空间和频率的局域性和方向性, 是“各向异性 (Anisotropy) ”的, 能获得图像的“稀疏”表示, 更好地反映图像的细节信息, 将Curvelet变换应用到图像融合中, 可以获得比小波变换效果更好的融合图像。针对Curvelet的优点, 本文提出了基于二代Curvelet变换的全色和多光谱图像的融合方法, 其步骤为:
(1) 首先, 利用经过重采样的多光谱图像对全色图像进行直方图匹配;
(2) 分别对多光谱图像和新的全色图像进行Curvelet变换, 得到各自的不同尺度层的系数;
(3) 对于Coarse尺度层系数 (即低频系数) , 选用基于区域能量的加权系数的自适应融合规则进行融合, 得到融合图像Coarse尺度层系数;
(4) 对于Detail尺度层系数, 选用以空间频率为度量标准结合自适应因子进行融合, 得到融合图像的Detail尺度层系数;
(5) Finest尺度层系数体现了图像的细节、边缘特征, 因此, 我们通过对全色图像的Finest尺度层系数进行加权 (系数扩大1.5倍) , 增强图像的边缘, 然后作为融合图像的Finest尺度层系数。
(6) 最后, 对得到的融合图像的系数进行Curvelet逆变换, 得到融合图像。
其过程如图2所示。
3 实验结果及其分析
本文对渭河流域的西安市一段的经过配准的多组IKONOS1m分辨率全色图像和4m分辨率多光谱图像进行融合实验, 本文的所有实验均在Matlab下完成。图3是其中一组实验图像, 图像大小为512×512像素, 256灰度级。
图3 (a) 为1m分辨率的IKONOS全色图像, 图3 (b) 是经过重采样的IKONOS多光谱图像。
本文用以下六种融合方法进行融合: (1) IHS变换融合方法; (2) PCA变换的融合方法; (3) 加权平均的融合方法; (4) 高通滤波变换法融合方法; (5) 小波变换的融合; (6) 本文方法。实验结果图像如图4所示。
本文对得到的融合图像从以下二个方面对六种融合方法得到的融合影像的质量做了评价。
3.1 空间质量评价准则
(a) 均值
对一幅图像来说, 均值指的是像素的灰度平均值, 对人眼为平均亮度。如果均值适中 (灰度值在128附近) , 则视觉效果良好。
(b) 信息熵
图像信息熵是衡量图像信息丰富程度的一个重要指标, 其定义为:
通过对图像信息熵的比较可以对比出图像的细节表现能力。熵的大小, 反映了图像携带的信息量的多少, 熵值越大, 说明融合图像携带的信息量越大。如果图像中所有灰度级出现概率越趋于相等, 则包含的信息量越趋于最大。
(c) 清晰度
图像的清晰度采用平均梯度法来衡量, 设平均梯度为:
在图像中, 某一方向的灰度级越大, 平均梯度值就越大, 否则值越小。平均梯度反映了图像对微小细节反差表达的能力何纹理变化特征, 可用来评价图像的模糊程度。一般来说,
3.2 光谱质量评价标准
(a) 偏差指数
偏差指数是用来比较融合图像和低分辨率多光谱图像偏离的程度的。融合图像的偏差指数定义为融合后图像的亮度分量I'与原始多光谱图像亮度分量了差值的绝对值与原始多光谱图像亮度分量I的比值, 定义为:
其中图像的大小为MxN。偏差指数的大小, 反映了融合结果保持的光谱程度, 偏差指数越大, 说明融合图像的光谱失真越大, 融合效果差。
(b) 光谱扭曲程度
它直接反映了多光谱图像的光谱失真程度, 光谱扭曲定义为:
其中V (i, j) 和V' (i, j) 分别为原始图像和融合图像在 (i, j) 点上的灰度值。扭曲程度的大小, 同样也是用来反映融合结果保持的光谱程度, 扭曲程度越小, 说明融合图像的光谱失真越小, 融合效果越好。
(c) 相关系数
融合图像F与原始图像A的相关系数能反映两幅图像的光谱特征的相似程度, 其定义如下:
式中,
图像的融合目的是为了最大限度地保留两幅图像中的重要信息, 因此对多光谱和高分辨率遥感图像来说, 既要保留原多光谱图像中的地物光谱信息, 又要引入高分辨率图像中的细节纹理信息, 为了客观评价图像融合结果, 需要综合考虑这两方面的因素, 实验取偏差指数、光谱扭曲程度、相关系数、标准差、熵、清晰度、等客观指标来比较各种融合方法的性能和优缺点, 结果如表2所示。
从表2中可以看出, IHS变换法的和其它方法比较起来, 偏差指数和光谱扭曲程度较小, 而相关系数较大, 说明IHS变换的融合图像和原始多光谱图像偏离较小, 保持的多光谱信息较大, 但很好地保持了空间信息。利用直接平均融合法得到的融合图像空间分辨率提高不大, 存在明显的光谱变异。利用高通滤波融合法可以得到增强的融合图像, 但融合图像的光谱还是产生了比较严重的失真现象。采用PCA融合的影像, 保留了原图像的高频信息, 目标的细部特征也更加清晰, 但PCA融合后地物边界有点发虚, 影像几何结构信息有一定的丢失。基于curverlet变换的融合方法得到的遥感图像的整体上在保持多光谱图像的光谱特征方面常用的遥感图像融合方法, 熵值最大, 说明它保持的多光谱信息丰富, 清晰度和空间频率最高, 说明很好地保持了空间信息, 同时图像更加清晰。
4 结束语
本文利用基于Curvelet变换的自适应多传感图像融合算法对高空间分辨率全色影像和低空间分辨率多光谱影像进行融合。对全色图像和多光谱图像图像做Curvelet变换, 分解成低频系数和高频系数, 分解后的低频系数选取区域能量的加权系数自适应融合规则, 高低频系数采用以空间频率为度量标准结合自适应因子进行自适应融合规则, 再对得到的低频和高频系数进行Curvelet反变换得到融合图像。与 IHS法、PCA法、拉氏法、小波法等传统融合方法相比较, 本文的方法不仅很好地保留了光谱信息, 而且较大地提高了多光谱影像的空间结构与细节信息, 具有良好的融合效果。
摘要：提出一种基于Curvelet变换的自适应遥感图像融合算法。首先对全色图像和多光谱图像图像做Curvelet变换, 分解成低频系数和高频系数, 分解后的低频系数选取区域能量的加权系数自适应融合规则, 高低频系数采用以空间频率为度量标准结合自适应因子进行自适应融合规则, 再对得到的低频和高频系数进行Curvelet反变换得到融合图像。实验结果表明该算法是一种有效可行遥感图像融合算法。
Curvelet 篇4
关键词：主成分分析,二代Curvelet变换,融合
图像融合是通过一些算法对各种不同传感器获取的图像信息进行综合利用。借助信息的互补性,融合后的图像可以满足更多应用需求,比原始图像更利于提取和分辨不同地物信息[1]。传统的遥感影像融合算法如PCA变换、小波变换等在进行遥感影像融合时都有一定的局限性。PCA变换会失去原有的物理特性。小波变换仅仅反映信号的点奇异性,也就是仅能反映奇异“点”的位置和特征,但在高维特征信息表达上存在困难[2]。
单一的融合方法不能取得令人满意的效果,所以要将多种融合方法的优势结合起来充分发挥其各自的优势。Curvelet变换是在小波变换基础上发展起来的一种新的多尺度变换,它的结构元素包括尺度和位置参数,比小波变换还多了方位参数,使得curvelet变换具有良好的方位特性。因此,curvelet变换对图像的边缘,如曲线、直线等几何特征的表达更加优于小波[3]。本文综合分析了PCA和curvelet变换的优缺点,提出了一种基于PCA和二代curvelet变换的新融合算法。为了验证该算法的有效性,对遥感图像进行仿真实验并对结果进行定量分析。结果表明这种方法不但增强了图像的空间信息,还保证了更多的光谱特性信息,可以获得好的融合效果,可用于图像分类和特征信息提取等方面。
1 Curvelet和PCA变换相结合的融合方法
1.1 Curvelet变换
Curvelet变换是由Candes和Donoho于1999年提出的[3],它的基础源自脊波理论[4,5,6]。该理论综合了脊波适于表达线特征和小波适于表达点特征,充分利用多尺度分析,所以适用于一系列图像处理问题并在应用中得到了较好的结果[7]。但是,第一代curvelet的数字实现比较复杂,需要子带分解、平滑分块、正规化和脊波分析等一系列步骤,而且curvelet金字塔的分解也带来了巨大的数据冗余量,因此,Donoho等人又于2002年提出了实现更简单、更便于理解的快速curvelet变换算法,即第二代curvelet[8]。
第二代curvelet与第一代curvelet在构造上已经完全不同,第一代curvelet的构造思想是通过足够小的分块将曲线近似到每个分块中的直线来看待,然后利用局部的脊波分析其特性。而第二代的curvelet和脊波理论并没有关系,实现过程也无需用到脊波,二者之间的相同点仅在于紧支撑,框架等抽象的数学意义。
本文提出的遥感影像融合算法采用基于Wrap的算法实现,即Wrap＿based curvelet,它的基本结构与离散curvelet变换的数字实现方式相一致,有点不同的是在每个尺度和角度下空间格网转变为曲波更为简单。我们采用规则矩形来代替倾斜格网,以同样的方式定义“笛卡尔”曲波为b:
式(1)中的S-θlTb由(k12-j,k22-j/2)代替,θ∈(-4π,4π)或θ∈(43π,45π)。
围绕原点Wrap,是Wrap＿based curvelet的核心思想点。见图1,意思是对任意区域,在具体实现时通过周期化技术映射到原点的仿射区域中,这种映射应该是一一对应的。将左上角的椭圆在原点位置重新装配映射完成后变为方形区域,这样才能用二维数组表示,才能利用二维傅立叶反变换。否则,在不装配成方形的的情况下,也没有办法对左上角的斜四边形进行二维傅立叶反变换。而一旦装配完后,计算量就大大增加了,除非实现零装配,但零装配的边界特性是非常之差的[9]。
1.2 PCA变换
PCA变换又称为主成分分析。它是建立在图像统计特征基础上的多维线性变换,具有方差信息浓缩,数据压缩的作用。变换后的第一主成分包含了总信息量的绝大部分,一般在80%以上,并且第一主成分相当于原来各波段的加权和,反映了地物总的辐射强度[10]。
PCA变换广泛应用于图像压缩,去除随机噪声,细部特征的增强和分析等。对低分辨率图像进行PCA变换并得到各成分变量,然后对高分辨率图像进行光谱拉伸,使它具有和低分辨率图像的第一主成分一致的均值和方差,用它替代低分辨图像变换后的第一主成分。最后进行PCA反变换到RGB空间,即融合图像。这种方法的优点在于它适用于多光谱图像的所有波段,但其有不足之处。由于在PCA融合算法中只是用高分辨图像来简单替代低分辨率图像的第一主成分。故会有低分辨率图像第一主成分分量中一些反映光谱特性的信息损失,因而使得融合图像的光谱畸变严重。
1.3 Curvelet和PCA变换相结合的融合方法
本文利用curvelet变换和PCA变换的优势,提出了一种遥感图像融合的新算法。具体流程如图2和图3所示。
2 实验
为了验证本文方法的有效性和适用性,选取Erdas软件实例中的dmtm.img和spots.img两幅影像作为数据源。图4中dmtm.img是由近红外波段、红波段和绿波段假彩色合成的低分辨率多光谱图像,spots.img图像是高分辨率全色波段图像。分别对实验区采用二维小波变换、PCA变换及本文提出的方法进行图像融合,结果如图4所示。在融合之前,已经对高分辨率图像和多光谱图像进行了几何校正等预处理操作。
3 实验结果分析
3.1 质量评价指标[1]
3.1.1 标准差
标准差是衡量一幅图像信息丰富程度的重要指标。
MN表示图像f的大小,f(i,j)表示图像f中像素点(i,j)的灰度值,f-表示对图像中所有像素值求平均,公式为
3.1.2 信息熵
图像信息熵是对图像信息量的测量,熵值越大,图像包含的信息量越丰富,图像质量就越好。
根据Shannon信息熵理论,一幅8 bit图像的熵为:
在这个公式中,Pi表示第i个像素值的概率。
后来Chavez(1984年美国)提出了多波段图像熵的表达方法最佳指数OIF。
stdi表示i波段标准差;Ri,j表示i和j波段的相关系数。
3.1.3 清晰度
图像清晰度用平均梯度衡量。
式(6)中
3.1.4 光谱扭曲度
Costantn等人用偏差指数来反映融合图像和原始图像的光谱信息匹配度。
g(i,j)是融合后图像g第(i,j)个像素的灰度值。
f(i,j)是原始多光谱图像f第(i,j)个像素的灰度值。
3.2 结果分析
从表1中我们可以看出,这三种方法各自不同波段统计出的标准差的值相差不大。比较它们的信息熵发现,本文提出的融合方法的OIF指数比PCA和小波变换融合后的值高,说明本文方法融合后的图像包含的信息量更大,信息对比更强而冗余小。在清晰度方面,某些波段PCA方法得到的值最高,本文方法次之,小波变换方法清晰度最差。清晰度反映了图像的变化率的程度。清晰度值越大,图像越清晰。这说明本文方法能较多地保留图像的细节纹理信息和空间结构特征,与PCA方法相比还达不到最优。最后比较光谱扭曲度参数,发现某些波段PCA方法值最高,本文方法次之,小波变换方法最小。光谱扭曲度值越小,说明融合图像增强分辨率和保留图像的光谱信息能力更强。因此本文方法的光谱扭曲度表现居中,能较好地保留图像的光谱信息但还有进一步改进的空间。比较上述所有参数,PCA变换和小波变换的融合方法都在某个方面表现较差,综合起来考虑,本文提出的PCA和curvelet变换相结合的融合方法可增强图像空间信息,保证更多的光谱特性信息,获得较好的融合效果。
4 结论
本文在分析二代curvelet变换和PCA变换特点的基础上,提出了一种curvelet变换和PCA变换相结合的遥感图像融合方法,将这两种方法各自的优势结合起来,克服其在遥感图像融合中的自身缺陷。实验结果表明,本文提出的方法可以增强多光谱影像的空间细节纹理信息,在保持数据光谱特性方面具有突出优势。但是,本文方法得到的融合结果仍存在部分细节纹理信息丢失的现象。由于curvelet变换对融合结果的影响主要体现在分解层数以及系数调整等方面,怎样更好地解决这些问题来增强图像的清晰度,减少光谱扭曲度,尽可能多地保留光谱信息,是本文需要继续研究的方向。
参考文献
[1]赵英时.遥感应用分析原理与方法.北京:科学出版社,2003
[2]张强,郭宝龙.基于Curvelet变换的多光谱图像与全色波段图像融合.系统工程与电子技术,2006;12(28):1186—1187
[3] Candes E J,Donoho D L.Curvelet-A surprisingly effective nonadaptive representation for objects with edges.In:Curve and SurfaceFitting:Saint-Malo 1999.Cohen A.,Rabut C.,Schumaker L L,eds.Nashville,TN:Van-derbilt Univ Press,1999
[4] Candes E J.Ridgelets:theory and applications.Department of Statis-tics,Standford University,1998;23—38
[5] Candes E J.Monoscale ridgelets for the Representation of image withedges.Department of Statistics,Stanford Univesity,Stanford,CA,1999;1—26
[6] Donoho D L,Flesia A G.Digital ridgelet transform based on trueridge functions.In:Stoeckler J,Welland G V,eds.Beyond Wavelets.Pittsburgh,PA,USA:Academic Press,2002:1—33
[7] Cheng Zhili,Wang Hongxia,Luo Yong.Wavelet’s theory and applica-tion.Beijing:Science Press,2004
[8] Candes E,Demanet L,Donoho D,et al.Fast discrete curvelet trans-forms.Multiscale Mode ling&Simulation,2006;5(3):861—899
[9] Candes E J,Demanet L,Donoho D L,et al.Fast discrete curvelet yr-ansforms.Applied and Computational Mathematics.California Insti-tute of Technology,2005:1—43
Curvelet 篇5
近年来, 小波理论在信号处理中的应用得到了很大的发展, 这主要得益于其对信号的时、频局域分析能力及其对一维有界变差函数类的最优逼近性能[1]。小波在表示具有点状奇异性的目标时, 体现出良好的性能, 而在表示线和超平面奇异性时并不是最优的。其主要原因是二维可分离小波基只有有限方向, 即水平、垂直、对角。方向性的缺乏使小波变换不能充分利用图像本身的几何正则性。为了克服这种局限, 多尺度几何分析应运而生, 出现了一系列多尺度几何变换的方法, 比如Ridgelet[2]变换、Curvelet变换[3]、Beamlet[4]变换等。尤其是Curvelet变换, 其基的支撑区间满足各向异性尺度关系 (width≈length2) , 可以很好地逼近图像中的奇异曲线。Candes和Donoho于1999年提出了第一代Curvelet变换理论, 它是由Ridgelet理论衍生而来的。由于其数字实现比较复杂, Candes等人又提出了实现更简单、更便于理解的快速Curvelet变换算法, 即第二代Curvelet变换。
与第一代Curvelet变换相比, 第二代Curvelet变换在实现过程中没有引入Ridgelet变换, 而是在频域中直接给出了Curvelet基的具体表示形式。
文献[5]中将Curvelet变换应用于图像去噪中并获得了良好的效果, 该方法在去除噪声的同时往往会“过扼杀”Curvelet系数, 导致图像的部分细节丢失, 边缘出现模糊。针对这一问题, 本文提出了一种基于Curvelet变换的阈值补偿去噪算法。该算法设置了两个阈值T和T0, T0是对“过扼杀”系数的补偿;对于Curvelet系数绝对值大于T的系数, 采用软阈值补偿硬阈值。试验表明, 该方法具有较高的输出信噪比和边缘保持效果。
1 Curvelet变换
在二维空间R2中, 定义空间位置变量为x, 频率域变量为ω, 频率域下极坐标变量为r和θ。令W (r) 和V (t) 为平滑、非负、实值的半径窗和角度窗, 其支撑区间分别为r∈ (1/2, 2) 和t∈[-1, 1], 且满足以下两个容许条件:
(1)
(2)
对于所有j≥j0的尺度, 定义频率窗:
(3)
其中⎣j/2」是j/2的整数部分。Uj为极坐标下的一种楔形窗, Uj的支撑区间是受W和V支撑区间限制获得的楔形区域, 楔形区域符合各向异性尺度的特性。令母曲波 (mother Curvelet) 为φj (x) , 其傅立叶变换, 则在2-j尺度上的Curvelet均可由ϕj通过旋转和平移获得。引入均匀的旋转角度序列:θl=2π·2-⎣j/2」·l, l=0, 1, …, 0≤ϑl<2π和平移参数k= (k1, k2) ∈Z2。
综合以上概念, 定义在尺度2-j, 方向θl, (k1, k2) 平移参数处的Curvelet为:
ϕj, l, k (x) =ϕj (Rθl (x-x) ) (4)
其中:
x=R (k1·2-j, k2·2-j/2) (5)
Rθl由θl旋转获得。
Curvelet变换便可以表示为:
(6)
2 基于阈值的图像去噪方法
基于小波变换的阈值去噪算法是现在运用最广泛的小波去噪方法, 它主要根据信号和噪声在小波变换后的不同特性进行去噪。小波变换能将信号的能量集中到少数小波系数上, 而白噪声在变换后仍是白噪声, 而且有相同的幅值。相对而言, 信号的小波系数必然大于那些能量分散的且幅值较小的噪声小波系数。根据这些特性, 可以选取合适的阈值, 对小波系数进行阈值处理得到估计系数, 然后进行反变换, 就可以实现图像去噪和重建。
利用Curvelet进行阈值去噪的基本思想与小波变换去噪方法一致。假设X (x, y) 表示无噪图像, N (x, y) 表示高斯白噪声, 则加性噪声图像Y (x, y) 可以表示为:Y (x, y) =X (x, y) +N (x, y) 。对噪声图像Y (x, y) 进行Curvelet变换生成Curvelet系数CD (j, l, k) , 估计阈值T, 将Curvelet系数CD (j, l, k) 应用阈值函数处理生成, 最后进行Curvelet反变换重建图像 (x, y) 。
常用的阈值处理方法包括硬阈值函数和软阈值函数。硬阈值函数是将一部分系数置为0, 保留剩下的部分;软阈值函数则是将一部分系数置为0, 其余的进行萎缩。为了表述方便, 将CD (j, l, k) 记作CD, 记作。其表达式如下:
硬阈值函数:
(7)
软阈值函数:
(8)
其中, sgn符号函数, 即
在此基础上, 人们提出了软硬阈值折中法:
(9)
其中, μ∈[0, 1], 为常数。当μ取0和1时, 分别对应硬阈值法和软阈值法。
3 改进的Curvelet去噪方法
3.1 阈值补偿去噪方法
用软阈值法去噪后的图像比较柔和, 但往往残留不少轻微噪声。使用Curvelet变换去噪时, 残留的噪声表现为长条状, 而长条状噪声严重影响图像的视觉效果, 因此采用硬阈值法更为适宜[6]。而硬阈值法处理Curvelet系数CD时, 由于在±T处不连续, 重构图像容易产生伪吉布斯现象;与此同时, 硬阈值往往还会造成“过扼杀”Curvelet系数, 导致图像的部分细节丢失。为克服这种缺点, 本文将软硬阈值结合起来, 构造一类新的阈值函数:
(10)
其中, -1<α<1;T0=βT, 0<β<1;λ=T0/T。
传统的软硬阈值折中算法在|CD|<T时, 都将Curvelet系数置为零, 这很可能导致有用信息被作为噪声滤除, 产生“过扼杀”现象。采用本文的算法, 引入一个新的阈值T0, 对在T0和T之间的Curvelet系数进行萎缩处理, 既防止将有用的信息滤除, 又可以保持边缘细节, 是对硬阈值函数的补偿。
Curvelet系数|CD|≥T时, 当α趋于0时, 为硬阈值法;当α趋于1时, 为软阈值法。新的阈值函数兼有软硬阈值法的优点, 只要选取合适的α, 就能得到较好的去噪效果。在试验过程中, α取负值时去噪效果更好。改进后的阈值函数仍然保留软硬阈值的优点, 它增强了硬阈值对于PSNR的贡献, 其系数大于1, 提高了重构图像的客观评价标准;又结合软阈值使图像比较柔和, 其系数小于0, 对硬阈值起均衡补偿的作用, 同时还提高了重构图像的主观评价标准, 所以称改进的阈值方法为阈值补偿法。
3.2 算法描述
综上所述, 本文的去噪算法可以归纳为:
(1) 选择合适的分解层数, 将含噪图像Y (x, y) 进行Curvelet正变换, 得到相应的Curvelet分解系数CD;
(2) 对噪声N (x, y) 进行阈值估计, 得到T和T0;
(3) 对Curvelet系数CD进行阈值处理, 得到;
(4) 对经阈值处理过的Curvelet系数进行Curvelet反变换, 得到去噪后的重建图像。
4 试验与分析
为了验证本文去噪算法的有效性, 进行实验, 选取PSNR和主观视觉效果作为去噪效果的度量标准。实验采用256×256×8bit的灰度图像Lena、Barbara、Peppers和Cell作为测试图像, 噪声是均值为0的加性高斯白噪声。对比了四种去噪算法:小波去噪算法、硬阈值算法、软硬阈值折中算法、阈值补偿算法。并且对阈值补偿算法和硬阈值算法的性能作了比较。
实验中阈值T=3k (在最细尺度上T=4k) 。k=Mean (abs (A-Mean (A) ) ) /0.6745, A是在不同尺度上噪声的Curvelet系数矩阵, Mean是对A取中值。实验给出在不同去噪算法下的PSNR值 (如表1所示) 和Lena图像在不同噪声条件下各算法的PSNR值 (如表2所示) 。其中软硬阈值折中算法μ=0.01;阈值补偿算法α=-0.01, β=0.91;性能提高是阈值补偿算法和硬阈值算法相比。
可以看出, 相同噪声条件下小波去噪算法的PSNR比基于Curvelet变换的去噪算法要低一些;和硬阈值相比, 折中算法PSNR有一些降低;相对于折中算法降低的值而言, 阈值补偿算法的PSNR则有一定提高。同一图像在不同的噪声条件下, 阈值补偿算法得到的PSNR是最高的。试验表明基于Curvelet变换的阈值去噪算法要优于传统的小波去噪算法, 本文的阈值补偿算法要优于前两种Curvelet阈值去噪算法。
为了进一步说明去噪效果, 给出了Lena图像在σ=20去噪后的视觉效果图。图1为原始图像, 图2为含噪图像, 图3为小波去噪图像, 图4、5、6分别为基于Curvelet变换的三种去噪算法后的重构图像。可以看出, 图3的主观视觉效果比Curvelet去噪图像效果要差;图4有较明显的伪影;与图4相比, 由于图5采用了折中算法, 图像就比较柔和;图6有效地抑制了伪影的同时, 还较好地保持了图像的边缘特征。
对于一幅大小为N×N的图像, 小波变换的时间复杂度为O (N2) , Curvelet变换的时间复杂度为O (N2logN) ;与硬阈值算法相比, 补偿阈值算法的时间复杂度与其大小一样。在时间复杂度增加不大的情况下, Curvelet变换去噪优于小波变换;阈值补偿算法优于硬阈值, 而时间复杂度却没有增加。从时间复杂度上可以看出, 基于Curvelet变换的阈值补偿去噪算法是简单有效的。
5 结束语
本文提出了一种基于Curvelet变换的阈值补偿图像去噪方法, 与传统的阈值函数相比, 它结合了软硬阈值函数的优点。即克服了硬阈值法“过扼杀”Curvelet系数的缺点, 又有效地抑制了伪吉布斯现象。试验表明, 新的阈值函数可以有效去除高斯白噪声, 重构图像不但在视觉效果上有所改善, PSNR也得到了提高。
参考文献
[1]焦李成, 谭山.图像的多尺度几何分析:回顾和展望[J].电子学报, 2003, 31 (12A) :1975-1981.
[2]Candes E J.Ridgelets:Theory and Applications[D].USA:Departmentof Statistics.Stanford University, 1998:23-28.
[3]Candes E J, Demanet L, Donoho D L.Fast Discrete Curvelet Transforms[R].Applied and Computational Mathematics.California Institute ofTechnology, 2005:1-43.
[4]Donoho D L, Xiaoming H.Beamlets and multiscale image analysis[R].2001, Technical repotr, Stanford Univ.
[5]Starck J L, Candes E J, Donoho D L.The curvelet transform for imagedenoising[J].IEEE Transaction on Image Processing, 2002, 11 (6) :131-141.
Curvelet 篇6
关键词：超声图像,斑点噪声,去噪,curvelet变换,自适应阈值
0 引言
超声诊断相较于其它的医学影像诊断, 具有简单、安全、有效、重复性强、价格低等特点, 因此超声诊断已经成为医学临床诊断的重要手段之一。在超声成像中, 当人体组织的结构尺寸与入射超声波波长相近或小于波长时, 超声束发生散射, 相位不同的散射回波相互干涉产生斑点噪声, 以至于图像对比度低, 给超声图像的后期处理带来很大的困难[1]。由于在超声诊断中主要依靠丰富经验的医生靠肉眼来对超声图像进行判断, 然而B超图像中的斑点噪声不仅降低超声图像中进行病灶分割和匹配的速度与准确率, 而且严重影响医生对正常和病变组织的识别能力。
由于超声图像噪声主要是以斑点噪声的形式出现, 斑点噪声是一种与图像信号无关的乘性噪声, 其存在使得超声图像信噪比低, 图像变得模糊不清。斑点噪声主要分布在高频部分, 而超声图像的一些重要特征信息如边缘、纹理等细节信息也是位于高频部分, 传统滤除高频部分噪声的方法大多针对的是加性噪声, 会对图像的相关细节信息造成破坏, 并不适用于超声图像去噪。因此针对超声图像去噪具有的特殊性与难度, 在进行相关去噪的同时应尽可能地保留图像边缘信息, 在不降低图像的分辨率的前提下, 最大限度的抑制斑点噪声是超声图像去噪的目标。
现有去除斑点噪声的算法大体分为三类:一是空间域局部统计滤波算法 (Lee滤波和Kuan滤波) [2], 该类算法主要是基于中心像素及周围像素的统计关系进行的, 窗口尺寸越大就越平滑, 在一些具有复杂结构的图像中很难在平滑和保留细节之间寻找到平衡;二是各向异性扩散滤波算法 (PM模型) [3], PM模型主要运用图像梯度的单调减函数来表示扩散系数, 而由于含噪图像中梯度具有很大的不稳定性, 会随着图像的平滑程度的增加而下降, 因此在去噪效果上并不理想;三是基于多尺度变换的滤波算法 (小波) [4], 二维小波变换的基是各向同性的, 无法很好地表达边缘的信息, 这就使得传统小波变换在处理二维图像时表现出一定的局限性。Candes等在1998年提出一种新的多尺度变换方法, 即Ridgelet变换, 对于具有超平面奇异性的高维多变量函数具有良好的逼近能力, 能够有效地处理高维信号的奇异性, 较好地实现对此类信号的逼近。随后Candes和Donoho等又在Ridgelet变换的基础上提出了curvelet变换, 相对于小波变换这一多尺度分析工具, curvelet变换的最大特点是具有高度的各向异性, 因此具有更强的表达图像中沿边缘信息的能力, 其各向异性特征非常有利于图像边缘的高效表示, 因此利用curvelet变换对超声图像进行相关去噪具有独特的优势。
1 curvelet去噪算法
Candes和Donoho在Ridgelet变换的基础上提出了curvelet变换并且构造了curvelet的紧框架, 对于具有光滑奇异性曲线的目标函数, curvelet提供了高效、稳定以及接近最优的表示, 相对于小波变换而言其最大特点是具有高度的各向异性, 是一种具有方向性、带通、多分辨的函数分析方法。curvelet变换直接以边缘为基本表示元素, 并且是各向异性的, 具有很强的方向性, 非常有利于图像边缘的高效表示[5]。从概念上讲, curvelet变换是多尺度金字塔的, 它在每个尺度有很多方向和位置元素, 而这些元素的几何多尺度特性使它与传统的像小波这样的多尺度表示法隔离开, 新框架直接从频域进行多尺度分析, 不再依赖这些几何特性。新的curvelet框架直接从频域进行分析, 不再通过脊波变换实现, 被称为第二代curvelet变换。
由文献[6]可知, curvelet变换可表示为:
并且其离散curvelet系数可表示如下:
该离散curvelet系数cD ( j, l, k) 可以被划分成Coarse层、Detail层和Fine层三个部分。从频率分布上看, Coarse层是由低频系数组成, 包含了图像的概貌;Detail层是由中高频系数组成, 主要包含的是边缘特征;Fine层是由高频系数组成, 体现的是图像的细节、边缘特征。由curvelet系数的分布和超声图像斑点噪声分布在高频部分的特点可以知道, 基于curvelet变换去噪的重点在于如何选取一个合适的阈值对Detail层和Fine层的系数进行处理, 尽可能去除噪声同时又保证超声图像的细节信息不被破坏。
2 基于自适应模糊阈值的curvelet去噪算法
阈值的确定是本文去噪算法的关键问题, 合适的阈值应能够兼顾平滑 (去除噪声) 和拟合 (与原始图像的近似程度) 两方面的要求。当选取的阈值过大时, 虽然能有效地去除噪声, 但会产生过多的零系数, 这样便破坏了原始图像的奇异性结构, 造成伪影和模糊;当选取的阈值过小时, 虽然会更好地和原始信号接近, 但去噪又变得不彻底。传统curvelet去噪算法的阈值确定, 往往需要未受噪声干扰的图像信息, 或者估计噪声的方差, 这些信息可以通过多次实验分析得到, 但在实际应用中就有一定的局限性。因此根据图像和噪声经过curvelet分解后, 在不同尺度和方向的系数的不同特点, 本文提出一种基于自适应模糊阈值的超声图像去噪算法。
本文借鉴Shark[7]提出的模糊小波阈值去噪的方法, 并结合局部方差确定模糊区域, 再确定其阈值进行去噪。局部方差大则表示是在图像的边界、纹理等突变点上, 信号对其有主要的影响, 局部方差小则表示在图像的平滑区域, 噪声对其有主要的影响[8]。
模糊性是指事物在形态及类属方面有不明确性, 其根源在于各类相似的事物之间的变换是连续的, 是由一系列的渐进过渡形成的, 不是断续的, 无联系的。这些过渡形式互相渗透, 互相贯通, 使得两个不同但相似的概念间没有明确的分界线。
所谓给定了论域U上的一个模糊子集A是指:对任何u∈U, 都指定了一个数uA (u) ∈[0, 1]与之对应, 称uA (u) 为u对A的隶属度。这意味构造一个映射:
这个映射称为其模糊子集的隶属度函数。
首先, 利用Donoho提出的多尺度统一阈值去噪方法, 得到阈值
其中N 2为第m层的信号长度, 噪声方差可以用中值估计法得到即
其中y为该层的curvelet系数, 并构建一个阈值区域[CTm, Tm]其中C∈ (0, 1) 。然后根据文献[9]求出其各层系数的边缘标准差
其中为含噪信号的curvelet系数的方差估计
其次, 构建方差的模糊区域根据对噪声占主体的curvelet系数应当尽量遏制, 并对信号占主体部分的系数应当适当保留这一去噪要求。当局部方差较大时, 则认为该点属于图像的边界或纹理上, 信号对其有主要的影响, 相应的隶属度d (s) 2较小;而当局部方差较小时, 认为该点属于图像的平滑区域, 噪声对其有主要的影响, 相应的隶属度d (s) 2较大。故本文提出采用梯形分布的隶属度函数, 其表达式为:
其中s2为局部方差, a, b为常量。
局部方差将采用文献[9]的方法, 其表达式为:
其中Zi, j 是以 (i, j) 为中心的移动窗, L为奇数。
最后根据隶属度, 可确定其自适应阈值为:
当隶属度d (s ) 2较小时其局部方差越大, 此时在图像的边界、纹理等突变点上, 信号对其有主要的影响, 相应的阈值较小;而当隶属度d (s) 2较大时其局部方差小, 此时图像的平滑区域, 噪声对其有主要的影响, 相应的阈值较大。这样, 通过隶属度确定的自适应阈值可以尽可能的滤除噪声, 同时又能保留较多的细节信息。
本文的去噪算法过程如下:首先对含噪的超声图像进行curvelet变换可得到curvelet系数cD ( j, l, k) ;由于噪声主要分布在系数的Detail和Fine层, 我们在Detail层和Fine层构建阈值区域结合局部方差确定其隶属度函数, 确定各层各系数的阈值;接着使用软阈值处理规则对各系数进行处理确定最后的curvelet系数;最后将处理后的curvelet系数进行curvelet反变换, 得到去噪后的超声图像。
3 实验结果及分析
为验证算法效果, 本文使用临床超声图像进行仿真实验。在实验测试中, 我们采用某进口高档B超仪器获取B超图像作为参考标准, 如图1 (a) 所示;以较低档的B超仪器获取图像进行斑点噪声处理, 如图1 (b) 所示;下面我们分别使用中值滤波、小波、基于curvelet的硬阈值处理和本文提出的基于curvelet的自适应阈值算法对超声图像进行去噪, 从而对各算法进行比较, 测试结果如图1 (c) 至图1 (f) 所示。
从图中1 (c) 可以看出, 中值滤波的去噪效果并不明显, 中值滤波的性能受滤波窗口的影响较大, 当滤波窗口较小时能较好地保护图像的边缘信息, 但去噪能力受到抑制, 而当滤波窗口增大时虽能较好地去除噪声, 但会破坏图像的边缘信息;小波去噪的原理就是在小波分解结构中保留低频分量, 对高频分量进行量化处理, 在新生成的保留低频分量和量化后高频分量的分解结构图的基础上, 再利用小波重构算法进行图像重构。但是二维小波变换的基是各向同性的, 变换系数的局部模极大值只能反映出这个小波系数出现位置是“过”边缘的, 而无法表达“沿”边缘的信息。由图1 (d) 可看出小波算法造成图像边缘模糊, 损失大量细节, 这就使得传统小波在图像去噪时具有一定的局限性。curvelet变换可以较好地解决小波只能恢复含水平和垂直方向的噪声图像这一局限性, 但传统的curvelet的去噪算法如硬阈值法、软阈值法、块阈值法的关键问题在于选择合适的阈值, 其决定着图像的去噪效果。从图1 (e) 中可以看出, 传统的curvelet的去噪算法虽较好地保留下了图像的边缘, 但去噪效果也并不显著;而从图1 (f) 可以看出本文提出的基于curvelet变换的自适应阈值去噪算法, 能够解决传统curvelet去噪算法在阈值确定上往往需要未受噪声干扰的图像信息, 或者估计噪声的方差这一局限性, 具有更好地实际应用意义。
为了更好地定性与定量评估本文算法的性能, 我们将低档B超仪器获取的超声图像 (图1 (b) ) 加入不同程度的斑点噪声, 并用四种算法对其进行去噪, 并通过定义峰值信噪比 (Peak Signal to Noise Ratio, PSNR) 对四种去噪算法进行比较。
其中, L表示图像量化级数, PSNR的值越大, 就说明算法的去噪效果越好。
比较结果如下表所示, VAR表示的是加入斑点噪声的方差, VAR=0表示低档B超仪器获取的超声图像 (图1 (b) ) , 而VAR=0.2, VAR=0.4则分别是在VAR=0的图像上加入不同噪声方差而获得的图像。
从表1可以看到, 不同的去噪算法都能够较好的提高带噪声超声图像的PSNR, 而随着噪声的增加, 去噪后的图像PSNR也逐渐减小。由于小波造成图像边缘模糊, 因此其PSNR在四种算法中较低, 去噪效果较差。相比于中值滤波、小波和传统的curvelet算法, 本文提出的基于curvelet的自适应阈值去噪算法的PSNR有一定程度的提高, 也说明了本文算法在超声图像的去噪效果上具有较明显的优势。
5 总结
为改善传统医学超声图像去噪算法容易丢失边缘和纹理等细节信息这些缺点, 本文提出了一种基于curvelet变换的自适应阈值的超声图像去噪声算法。该算法利用curvelet系数局部方差在超声图像纹理与平滑度的差异, 结合模糊数学中的模糊区域和隶属度函数的概念确定各curvelet系数的自适应阈值, 从而实现对超声图像的去噪。经实验表明本文算法在医学超声图像斑点噪声的抑制上有良好的效果, 且具有较好的主观视觉, 可以改善低档B超仪器在获取超声图像时带有较大斑点噪声的缺点, 具有实际的应用意义。
参考文献
[1]乔蕊, 侯燕.医学超声图像有效去噪方法研究与仿真[J].计算机仿真, 2013, 30 (5) :391-394.
[2]Zhong H, Zhang J, Liu G.Robust polarimetric SAR despeckling based on nonlocal means and distributed Lee filter[J].IEEE Trans Geosci Remot Sens, 2014, 52 (7) :4168-4210.
[3]Ham B, Min D, Sohn K.Revisiting the relationship between adaptive smoothing and anisotropic diffusion with modified filters[J].IEEE Trans Imag Process, 2013, 22 (3) :1096-1107.
[4]Quellec G, Lamard M, Cazuguel G, et al.Fast wavelet-based image characterization for highly adaptive image retrieval[J].IEEE Trans Imag Process, 2012, 21 (4) :1613-1623.
[5]Cand SE, Demanet L, Donoho D, et al.Fast discrete curvelet transforms[J].Multiscal Model Simulat, 2006, 5 (3) :861-899.
[6]Wang HZ, Qian LY, Zhao JT.An image denoising method based on fast discrete curvelet transform and total variation[C].Int Conf Signal Process, 2010, 1040-1043.
[7]Shark LK, Yu C.Denoising by optimal fuzzy thresholding in wavelet domain[J].Electronics lett, 2000, 36 (6) :581-582.
[8]关新平, 刘冬, 唐英干, 等.基于局部方差的模糊小波阈值图像去噪[J].系统工程与电子技术, 2006, 28 (5) :650-653.
Curvelet 篇7
1 快速全变差图像复原算法(FTVd G)
设给定的退化图像f(x,y),模糊算子H,噪声n,原始图像u(x,y)(待求图像),梯度算子Di则全变分图像恢复的基本模型
式中,∑||Diμ||是离散TV,该模型就是通常TV/L2模型这个模型也称为非线性去卷积模型。基于这个变分模型,Wotao Yin等人提出一种交替最小化图像恢复算法——FTVd G算法[4],在去模糊方面具有很好的效果。该算法来自一种新的半二次模型,它适用于各向同性和各向异性的TV范数,每一步迭代计算的复杂性在于三次快速傅里叶变换,具有稳健的收敛性。根据最优化的变量分裂和惩罚技术,可以引入一个附加变量,在每个像素点i产生一个值wi来改变非可微Diu二范数,且在Diu和wi加个惩罚项,从而得到双变量优化模型:
式(2)实际上是两个最优化问题其中,固定其中一个,就能求另一个,可以交替进行,得到一系列迭代值。FTVd G算法的思想主要是由于原TV项计算起来比较复杂,从而引入一个新的TV项容易计算的变量,把一个最优化问题变成两个最优化问题来求解,且两个保真项有解析解,可以用快速傅里叶计算,所以计算的速度比一般的TV恢复算法要快的多。但是它对噪声比较敏感,且在每次迭代中会放大原来的噪声,容易产生一些新的人工噪声。当噪声水平加大时,在恢复过程中噪声会被放大,因此复原效果就比较差。
2 基于Curvelet阈值收缩的TV图像复原算法(Cu TVd G)
2.1 CuTVdG图像复原模型
基于稀疏性约束能很好的实现去噪的特性,我们可以利用图像信号稀疏特性来压制和消除每次迭代产生的噪声,所以我们提出一个新的图像复原模型(TV-l1)模型:
式(3)模型由两个正则项和一个保真项组成,第一项是个TV正则项,代表图像边缘信息的先验知识;第二项是稀疏性正则项,代表图像在某个变换域(Ψ)中的稀疏性的先验知识,主要作用是压制噪声,本文的Ψ选用Curvelet变换,因为Curvelet变换比小波变换具有更好的稀疏性,而且Curvelet变换对奇异性信号具有很好的表现力;第三项是图像的保真项。实际上这个新的模型是在TV去模糊复原模型再加上一个基于稀疏性约束的去噪模型,这样在保真图像时,又能去除噪声。整个模型由两个正则项和一个保真项组成,这样能很好的复原带有噪声的图像。
根据FTVd G算法的思想,我们可以把TV正则项进行分裂,引入一个新的变量w,把一个最优化问题变成两个最优化问题,然后交替求解。这样式(3)又变成下面这个双变量变分模型:
实际上,模型(4)可以分为两个最优化子问题,即分别关于w,u的最小化问题,固定其中一个,来求另一个。
交替迭代w,u,算出最终u的解,第一个最小化问题本质上是个TV问题,可以按照TV的算法进行求解,第二个最小化问题是个基于稀疏性的图像复原问题,它的前两项是个二次项,如果单独求可以得到解析解,但后面有个l1范数,所以必须把它与前两项结合起来求解,这种情况可以用算子分裂法求解。
2.2 Cu TVd G图像复原模型的求解算法
2.2.1 算子分裂法
算子分裂法是凸分析中的经典算法,通过分解策略将原先复杂的问题转化较为简单子问题的迭代求解,从而降低了问题的难度,在求解复杂优化问题时得到了广泛的应用。算子分裂法采取分而治之的策略,它将一个复杂问题转化较为简单子问题的迭代求解,每一个子问题只包含一个算子。
令基于稀疏性约束的图像复原变分模型为:
对于约束z=Ψx,引入拉格朗日乘子p,得到模型(7)的拉格朗日乘子函数:
根据凸共轭的定义,则可以得到变分问题的对偶模型:
式中p为对偶变量,我们可以看出对偶模型中保真项中含有不可分离项Ψp,如果保真项可微,那么就可以对原变分问题的进行求解。令算子
对算子A和B我们将采用Dual Douglas-Rachford[5]分裂法求解,得到如下的迭代公式:
2.1.2 CuTVdG模型数值计算
变分模型(5)是一个TV求解问题,我们可以对此进行求导,得到一个Euler_Lagrange方程,最终的计算方法如下式所示:
式中,||Diμ||为梯度的模,实际上可以把上式看做一个阈值收缩过程。
对于变分模型(6)是个比较复杂的问题,可以借助Dual Douglas-Rachford算子分裂法进行求解则有
式中的Ψ表示Curvelet变换[6],根据公式(10),分别对uk+1和zk+1进行变分求解,再对pk+1进行更新处理。
令,其中D(1),D(2)分别表示对x轴方向和y轴方向进行差分的差分算子矩阵,wi=((v1)i,(v2)i)∈R2,表示式(11)在点i所求的离散值,则有:
其中,F,F-1,Ψ分别表示傅里叶变换、逆傅里叶变换和Curvelet变换,*表示复共轭,°表示矩阵对应位置的元素相乘,除法也是如此,式中,令Ψ-1Ψ=ηI且η>0。
这里的是个阈值收缩函数,式(15)的作用就进行Curvelet软阈值收缩处理,对应的是逐个元素的软阈值收缩。p的迭代格式也比较好求,它的解如下
综合上述各个子问题的数值求解,可以得出新的TV耦合模型的数值计算算法。整个算法实际上是一个交替迭代的过程,最后算法收敛停止迭代。
对于u的计算我们在傅里叶域求解,这样能大大提高计算的速度。同时把稀疏表示锁定在Curvelet域,因为Curvelet变换具有各向异性,是完备的,它的基本表示元素是边缘,而且能非常准确地表示图像的边缘方向,方向性极强。使用Curvelet阈值收缩法又能很好的压制图像的噪声,消除在去模糊中带来的放大的噪声。
3 实验与分析
本文选用Lena图作为实验对象,对原始图像分别加入尺寸为,标准差为10的高斯模糊,再分别加入标准差是0.0255白噪声,将从视觉效果、峰值信噪比(PSNR)、结构相似度指数(SSIM)、迭代次数(Iteration)和相当误差(Re Err)方面进行对比。
FTVd G算法在低PSNR时就停止迭代,说明FTVd G算法压制噪声能力低;而Cu TVd G算法为了得到更高的PSNR继续迭代,这样做能很好的压制和消除噪声,所以尽管多迭代几次,但是能得到更好的复原效果。
4 结论
本文给出一个基于全变差(TV)正则化和稀疏性约束的耦合图像复原(TV-l1)模型。该模型是在全变差图像模型的基础上加上图像Curvelet变换下l1稀疏性约束,达到图像边缘结构和纹理特征保持的图像复原。同时针对该最优化模型的求解问题,基于算子分裂法原理,设计了一种多步迭代的数值算法。实验证明本文算法复原图像的视觉质量优于快速TV复原算法(FTVd G)的复原结果。
参考文献
[1]Rudin L,Osher S,Fatemi E.Nonliner total variation noise removal algorithms[J].Phys.D,1992,60:259-268.
[2]Jamylle Lauri Carter.Dual methods for TV-based image restoration[D].UCLA CAM,2002,2-13.
[3]Aujol J F,Gilboa G.Implementation and parameter selection for BV-Hilbert space regularizations[R].UCLA CAM Report,2004:4-66.
[4]Wang Y,Yang J,Yin W.A new alternating minimization algorithm for total variation image reconstruction[J].SIAM J.Imaging Sci,2008(1):248-272.
[5]Combettes P L,Pesquet J C.A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[J].IEEE J.SelectedTopics Signal Processing,Vol.1,2007:564-574.
[6]Candes E J,Demanet L,Donoho D,et al.Fast discrete curvelet transforms[R].Pasadena,California:California institute of Technology,2005:1-43.
[7]Daubechies I,Defrise M,DeMol C.An iterative thresholding algorithm for linear inverse problems[J].Comm,Pure Appl,Math.,2004(11):1413-1457.
【Curvelet】推荐阅读：