浅析curvelet变换原理与理解 - 全文

　　Curvelet变换是基于傅里叶变换和小波变换的一种改进，其特点是有高度的各向异性，具有良好表达图形沿边缘的信息的能力，对于恢复形状的沿边缘的主要结构和抑制周边噪声有其特有优势。

　　其过程为

　　这和传统的DFT及小波变换的处理过程类似，把图表中的curvelet换成DFT和wavelet就可以了。

　　Curvelet变换是最近图像处理较新的一种多尺度几何变换算法。其发展历程在短短十年间：

　　1999年，Candès和Donoho在Ridgelet变换的基础上提出了连续曲波（Curvelet）变换——第一代Curvelet变换中的Curvelet99。

　　2002年，Strack、Candès和Donoho提出了第一代Curvelet变换中的Curvelet02。

　　2002年，Candès等人提出了第二代Curvelet变换。

　　2005年，Candès提出了两种基于第二代Curvelet变换理论的快速离散实现方法：

　　1）非均匀空间抽样的二维FFT算法（Unequally-Spaced Fast Fourier Transform，USFFT）；

　　2）Wrap算法（Wrapping-Based Transform）

　　与小波变换类似，Curvelet变换同样有其对应公式。Curvelet系数可由下式得到，即信号与小波函数内积：

　　这里j表示尺度，l表示方向，k表示位移。变换的推导以及原理是个十分复杂的过程，这需要有相当强的数学功底。

　　Curvelet变换原理

　　Curvelet 变换通过对Radon 域内的每一个投影轴作一维小波分析，并由局部脊波分析得到多尺度结构。设参量θ是常数，平移量1是变量，脊波系数R，（a，b，θ）为

　　然后对每个分块用式（1，2）进行脊波分析。上面每步均可逆，经过逆变换可得重建图像。

　　通过Radon域的多角度投影，Curvelet变换明显加强了目标边界。同时由于Curvelet变换是多尺度脊波变换的子集，它可以用投影切片定理来实现：一幅n×n的图像可表达成2n条径向线的积分结果，对每条线的投影数据计算一维小波变换，就得到脊波系数。这样的Curvelet变换冗余因子达到16J+1。由于经典算法中分块的径向抽样角度恒定，直角坐标与极坐标的转换将导致在同一块内数据的非均匀抽样，当块中心附近的抽样值满足奈奎斯特抽样定理要求时，在块的边界部分将出现欠抽样;而当块边界部分满足抽样定理要求时，块中心部分数据将会过抽样，而且块的尺寸越大，非均匀抽样的影响就越明显。因此，在原算法中尽管离散Curvelet变换的实现较简单，但是为了无失真地重构图像，保证块边界周围的抽样值足够多，它在块的中心部分实施了过抽样，因此耗时和冗余度也相应增大。

　　下面是对Curvelet变换的理解：

　　在时间域：Curvelet变换可以看做是一个椭圆以内积形式依次覆盖整个矩阵（指要变换的矩阵），这个图在很多地方都出现过。就我个人的理解详细解释下。

　　这里的椭圆就是Curvelet变换的窗口，相当于小波变换的基函数。其长轴与短轴关系为平方关系，其窗口大小视尺度j而定。与二维小波变换最大的优势在于Curvelet变换具有方向性。就是在尺度j下，椭圆做完整个矩阵的内积后（即依次覆盖完整个矩阵），可得到一组系数。然后将椭圆进行旋转适当的角度a再与矩阵做内积得到a角度的Curvelet系数。这样将椭圆进行多次旋转，每旋转一个角度就能得到尺度j一个方向的系数矩阵。在尺度j下做L次旋转，那么尺度j下就产生L个角度的系数矩阵。做完尺度j了接下来做尺度j+1下所有角度，方法同上。

　　如一幅图像上有一段弧，椭圆长轴沿着弧方向做覆盖内积得到的系数就大。这样就适合做图像边缘检测。

　　对于Curvelet变换的Matlab程序包curvlab可在网上下载。Curvlab包里有Curvelet的快速离散算法的Matlab程序和C++程序。

　　将其添加在Matlab里面：File -》 set path

　　Add Folder 添加好curvlab里相关*. M文件的路径。

　　这时在Matlab里面使用Curvelet的快速变换及反变换函数：

　　fdct_usfft（）

　　ifdct_usfft（）

　　fdct_wrapping（）

　　ifdct_wrapping（）

　　在Matlab中Curvelet变换后返回的是一个cell矩阵。这个矩阵装着各个尺度、各个方向的系数值。

　　C= fdct_usfft（***）；

　　C{j}表示一个cell矩阵，装着尺度j上所有方向的系数。

　　C{j}{l}就表示是一个二维矩阵，表示尺度j，方向l上的所有系数，

　　一般尺度越高，其对应的都是高频系数。

　　curvelet matlab示例代码理解

　　1. fdct_wrapping

　　function C = fdct_wrapping（x， is_real， finest， nbscales， nbangles_coarse）

　　% fdct_wrapping.m - Fast Discrete Curvelet Transform via wedge wrapping - Version 1.0

　　%

　　% Inputs

　　% x M-by-N matrix 输入为MxN的矩阵

　　%

　　% Optional Inputs

　　% is_real Type of the transform 转化的类型

　　% 0： complex-valued curvelets 复数值的曲波变化

　　% 1： real-valued curvelets 实数值的曲波变化

　　% ［default set to 0］ 默认设置为0

　　% finest Chooses one of two possibilities for the coefficients at the

　　% finest level： 选择一种表示方式计算最优级的系数

　　% 1： curvelets 曲波变化

　　% 2： wavelets 小波变化

　　% ［default set to 2］ 默认设置为2

　　% nbscales number of scales including the coarsest wavelet level

　　% 包含最粗小波级在内的伸缩数

　　% ［default set to ceil（log2（min（M，N）） - 3）］

　　% nbangles_coarse

　　% number of angles at the 2nd coarsest level， minimum 8，

　　% 第二粗糙级的角度数，最小为8

　　% must be a multiple of 4. ［default set to 16］

　　% 必须为4的倍数，默认为16

　　% Outputs

　　% C Cell array of curvelet coefficients.

　　% C{j}{l}（k1，k2） is the coefficient at

　　% - scale j： integer， from finest to coarsest scale，

　　% 从最佳尺度到最粗尺度

　　% - angle l： integer， starts at the top-left corner and

　　% increases clockwise，

　　% 从左上角开始顺时针增长

　　% - position k1，k2： both integers， size varies with j

　　% and l.

　　% If is_real is 1， there are two types of curvelets，

　　% ‘cosine’ and ‘sine’。 For a given scale j， the ‘cosine’

　　% coefficients are stored in the first two quadrants （low

　　% values of l）， the ‘sine’ coefficients in the last two

　　% quadrants （high values of l）。

　　%

　　% See also ifdct_wrapping.m， fdct_wrapping_param.m

　　%

　　% By Laurent Demanet， 200412345678910111213141516171819202122232425262728293031323334353637383940414243

　　2. DCT基本示例代码理解

　　% fdct_wrapping_demo_basic.m -- Displays a curvelet both in the spatial and frequency domains.

　　m = 1024;

　　n = 1024;

　　X = zeros（m，n）;

　　%forward curvelet transform

　　disp（‘Take curvelet transform： fdct_wrapping’）;

　　tic; C = fdct_wrapping（X，0，2，8，64）; toc; %tic toc配合使用测量程序运行时间

　　%specify one curvelet

　　s = 7; %从1开始增大，空间域变细，频率域变粗

　　w = 10;%从1（左上角）开始增大，空间域顺时针旋转，与笛卡尔corona相对应

　　［A，B］ = size（C{s}{w}）;%尺度为s，方向为w，的矩阵大小

　　a = ceil（（A+1）/2）;

　　b = ceil（（B+1）/5）;

　　C{s}{w}（a，b） = 1; %该尺度、方向中心位置元素设置为1

　　%adjoint curvelet transform

　　disp（‘Take adjoint curvelet transform： ifdct_wrapping’）;

　　tic; Y = ifdct_wrapping（C，0）; toc;%进行反曲波变化，得到空间域图像

　　%display the curvelet

　　F = ifftshift（fft2（fftshift（Y）））;

　　subplot（1，2，1）; colormap gray; imagesc（real（Y））; axis（‘image’）; 。。。

　　title（‘a curvelet： spatial viewpoint’）;

　　subplot（1，2，2）; colormap gray; imagesc（abs（F））; axis（‘image’）; 。。。

　　title（‘a curvelet： frequency viewpoint’）;

　　%get parameters

　　［SX，SY，FX，FY，NX，NY］ = fdct_wrapping_param（C）;