[转帖]新小波Wedgelet（楔波）

simon21

　　1。小波的末路

　　千万别被此吓倒。因为当我看到这个小波广告牌的时候，也有些晕。实际上它刻画的是小波在二维图像处理中的致命伤。图象处理的要素被大师概括为：局部性，多分辨率，严格采样，方向性，各相异性(anisotropy)。关于前三者小波无疑是优秀的，但后两者就略显不足了。小波是点奇异的或是灰度奇异的，所以在一维信号处理中它绝对可以称得上是：powerful tool。但在二维图像处理中，小波对于图像几何特性的刻画就不那么令人振奋了。例如，一张黑色背景上面的白色实心圆，标准的cartoon PIC。(这里请注意悲剧的始作俑者：光滑的圆环边界。)如采用小波作压缩的话，你会发现必须采用多层小波分解，并且在明暗交接的边缘处很多高频分量的幅度是无法忽略或容忍的。回想起在图像处理中，小波也只有水平、垂直、对角线、低频四个方向信息，这也是远远不够的。数学上，对于光滑的二维C^2曲线(contour)，取最明显的N个小波系数逼近，它的收敛是N^(-1)而不是N^(-2)。所以，我们需要找到更多的反映图像光滑边界或边缘信息，并且找到反映图像多方向性信息的新小波。

　　2。两个方向

　　时域和频域永远是那么奇妙。Idea永远来源于这两端。在频域，一场在频谱上方向性正交分割的游戏开始了。ridgelet->curvelet->contourlet，我个人喜欢最后的。因为它是离散化的东西，更重要的是它告诉我们一个快被人遗忘的名字--方向性滤波器。于是fan filter，multirate system，polyphase domain，等的一切东西开始发热。包括完美滤波器的构造，广义的PR条件，梅花5点采样，二维提升，这些也给很多人造成了最大的麻烦--代码。于是，大家开始下载大师们提供的toolbox，这点无可厚非，只要是目的纯良。我个人仍然是桀骜不驯，从不喜欢拿来主义。但对于这个方向的难度，也深知不易，所以也就稍稍回避，有时间再战不迟，毕竟我的博士论文与小波毫无关系。在这时，另外一种新的时域基—Wedgelet，开始引起我的重视。熟悉的自适应四叉树，简单的直线分割，还有些说怪不怪的名字：dyadic square，rate-distortion (R-D) problem。毕竟我的专业不是图像处理。但觉得它的难度还行，核心程序估计不会超过3000行左右，也就开始新的挑战。

　　3。MATLAB V.S. Visual C++

　　这是至关重要的，想做Wedgelet的同学，记住C++是最好的选择。Why? 数据结构。C++的精髓实质上也就是数据结构。树这个东西，一定要用Link或Stack实现它。(但.net好像有了这个结构。)还有一个原因，速度。不怕告诉大家，在我初期写这个程序的时候，运行的时间是可怕的。但最终我优化后，也就控制在20秒左右了。

　　4。Wedgelet新基

　　什么是Wedgelet?说白了，就是在一个图像子块(dyadic square)画条线段，把它分成两个楔块，每一个楔块用唯一的灰度值表示。线的位置，两个灰度值，就近似刻画了这个子块的性质。

　　5。灰度的确定

　　我们把问题分开。假设线固定，我们怎样找出最好的灰度，让它能和子块最接近呢?大家都知道数学上是最小二乘拟和;物理上这个灰度其实就是那个被分割而成的楔块包含的所有像素的平均灰度值。

　　6。线段的噩梦

　　怎样能找到最佳的线段。每个人都知道穷举。但怎么穷举?如果盲目的对整个图像穷举，我们考虑到任何n*n大小的子块或划分，那它的工作量实际上是O(n^4)。所以，谁也不会傻到这样。

　　7。问题开始

　　任何东西，有得有失。我们把整个图像当作分解对象，我们会找到一条最优的长长的线段，把图像分成两块，同时找到最优的灰度。我们花很少的Bit就可以代表整个图像了，但后果是美女没有了，开个玩笑。就是说误差太大了，或者是PSNR太大了。但我们如果把图像弄成1*1的小块，连Wedgelet分解都不需要(每个子块就一个灰度)，还没有误差，但谁如果这样想，不就什么也没干嘛!所以我们要把PSNR和Bit率结合起来看，这就是fidelity vs. complexity 或者rate-distortion (R-D) 问题。写成数学就是A+factor*B问题，A就是误差，B就是Bit率，factor就是惩罚因子。上式就是代价函数。当factror->0，我们就只考虑误差最小，于是我们得到是整个完美的图像，但是没有进行任何的压缩。当factror->inf，我们只考虑Bit率，我们最终得到只有一个灰度级的图像，这个灰度级就是这个图像的灰度均值。所以有了factor，一切是可控的。问题就是给定factor，我们找到上面表达式的最小值。最后再说一句，A我们可以用绝对平方误差表示。B呢?实际上B我们可以用分解的块数表示。这是显然的，因为图像的Bit率是和块的多少成正比的。

　　8。解决方案

　　最优多树字典(multitree dictionary)。这个名词，我不想谈太多。我只谈它的特例，自适应四叉树。把大的子块分成四个小的子块，找到每子块的均值。比较四个小子块的代价函数之和与大子块的代价函数的大小。如果，前者小，分解进行。否则，分解中止。有人会问，这不就是自适应小波包算法吗?是的，除了代价函数小波包一般是“熵”外，其他的应该是一模一样的。最终分解的结果可以想象，在图像的平坦区(直流分量多)子块较大，在图像的纹理区(交流分量多)子块较小。

　　9。线段的解决

　　噩梦中止。其实穷举的线段应该和子块联系起来，在子块周界上离散化搜索。拿n*n子块而言，我们有4*n个离散的周界点。因此，穷举的复杂度O(n*n)。也就是说，我们先在整个图像上进行四叉树分解，然后再进行Wedgelet分解。这时，我们仍需判断分解的必要性。即若分解后的两个楔块的代价函数之和小于对应子块的代价函数，则分解进行。否则保留子块，不采用楔波代表。这里，我还要说一句，我们穷举所有直线的时候，必须判断像素和直线的关系，也就是说我们必须把属于不同楔块的像素分成两类，来获得每个楔块的均值。这里像素的中心点的位置，是判断的标准。有了直线方程，和每个像素中心点的位置，不会有人还说不能把像素分成两类了吧。呵呵，是不是有些像神经网络里面最简单的感知器，科学就是如此的相似与奇妙。

　　10。加速再加速

　　不要以为，噩梦中止后就不会有麻烦。即使这样，速度仍然慢的让你头晕，乃至痛哭流涕。于是，批处理开始了。换句话说，我们把同样大小的子块看作同一个子块，我们可以把线段的穷举一次完成，因为同样大小子块的周界有着相同的相对索引，构成它们的像素点也有着相同的相对坐标。所以，所有像素点的分类是一次完成的。这就是Key Point。

　　11。最后的话

　　这里，四叉树是可以扩展的。因为，我们分解图像的时候，可以不是一下分成四块，可以是两块。在横向或纵向的划分中找到最优的分割。如果这种分割，是正好把图像分成面积相等的两块，就是dyadic tree。如果是任意的话，就是多树字典了。关于编码，大家可以想想，这里我就不说了，但压缩绝对是Wedgelet的一个优势，甚至超过了JPEG2000。还有，小波域的Wedgelet是更有优势的东西，论文IEEE也出来了，大家可以看看。还有就是自适应的Wedgelet，采用马尔科夫预测的方法就更是神奇了，我由于专业所限，也就不深入研究了。希望大家，可以写写这个程序，会很有意思的。