小波去噪算法的简易实现及其扩展（小波锐化、高斯拉普拉斯金字塔去噪及锐化）之一。
早年就接触过小波的概念，那个时候看什么小波十讲这类的，看的可真谓云里雾里，一大堆数学公式，头大的要死。做去噪的时候也看很多人说小波去噪算法效果不错，不过网络上有的都是matlab代码，而matlab的小波包里的函数是已经写好的内嵌函数，是无法看到代码的。因此，一直以来，也从未想过自己动手写个小波去噪之类的效果。
偶尔翻阅了一下GIMP软件的菜单，再次看到了在其Filters-->Enhance菜单下有个wavelet-decompose菜单，点击一下，发现原图像是没有任何增强的效果的，但是在其图层界面里增加了一些列的图层，如下图所示：
后面搜索一些参考资料，大概明白了他的意识这个做分解是为后续的增强做铺垫的，因为他分解成了多个层后，可以单独对每个层进行一些特别的处理，GIMP官方的文档对其说明如下：
This filter decomposes the active layer or selection into several layers, named “scales”, each of them containing a particular set of details. Finest details are in first layers and they become larger until you get to the last one, at bottom. This last layer is called “residual” and holds what is left after all detail layers have been removed; it represents the global contrast and colors of the image.
Each of scale layers are set to combine using the Grain Merge layer mode. This means that pixels that have a 50% value will not affect the final result. So, painting a wavelet scale with neutral gray (R:50% G:50% B:50%) will erase details.
Wavelet-decompose is a wonderful filter for skin smoothing and retouching, removing blemishes, wrinkles, spots from your photos. It can be used also for sharpening and local contrast enhancement and for removing stains, colors, tones. All this is well explained in tutorials mentioned above.
这个帮助最后提到这个分解可以用于皮肤光滑、磨皮、移除瑕疵、斑点等，或者做锐化以及局部增强等等功能。
似乎很是强大。
仔细看看GIMP分解后的图，我们发现他将图像分解为了多个图层，图层的数量取决用户界面的参数，比如选择5层，他实际上是生成了6个图层，额外增加了一个特殊的Residual（残余）层，我们试着尝试解析他的代码。
在GIMP的源代码里搜索wavelet，可以发现gimp-master\plug-ins\common这个目录下有个wavelet-decompose.c文件，再打开这个文件，稍微分析下这个代码，发现其中需要一个非常核心的函数：wavelet_blur，这个函数确没有在gimp-master这个文件夹里，而是在gegl-master这里。wavelet_blur函数又涉及到一个wavelet-blur-1d的文件。
得益于早年我翻译和抽取过很多GIMP的函数，以及自己对图像处理本身算法的了解，虽然GIMP的代码写的很晦涩，但是拼接多年的经验，还是成功的把这个代码抽取出来。下面简要的分析下：
在wavelet-decompose.c里有一段核心的东西如下：
1 for (id = 0 ; id < wavelet_params.scales; id++) 2 { 3 GimpLayer *blur; 4 GimpLayer *tmp; 5 gchar scale_name[20]; 6 7 gimp_progress_update ((gdouble) id / (gdouble) wavelet_params.scales); 8 9 scale_layers[id] = new_scale; 10 11 g_snprintf (scale_name, sizeof (scale_name), _("Scale %d"), id + 1); 12 gimp_item_set_name (GIMP_ITEM (new_scale), scale_name); 13 14 tmp = gimp_layer_copy (new_scale); 15 gimp_image_insert_layer (image, tmp, parent, 16 gimp_image_get_item_position (image, 17 GIMP_ITEM (new_scale))); 18 wavelet_blur (GIMP_DRAWABLE (tmp), pow (2.0, id)); 19 20 blur = gimp_layer_copy (tmp); 21 gimp_image_insert_layer (image, blur, parent, 22 gimp_image_get_item_position (image, 23 GIMP_ITEM (tmp))); 24 25 gimp_layer_set_mode (tmp, grain_extract_mode); 26 new_scale = gimp_image_merge_down (image, tmp, 27 GIMP_EXPAND_AS_NECESSARY); 28 scale_layers[id] = new_scale; 29 30 gimp_item_set_visible (GIMP_ITEM (new_scale), FALSE); 31 32 new_scale = blur; 33 } 34 35 gimp_item_set_name (GIMP_ITEM (new_scale), _("Residual"));
明显这个循环就是要生成各个图层的内容的。
第14行从new_scale层拷贝数据到tmp，然后第18行进行一个wavelet_blur得到Blur图像 ，注意那个模糊的最后一个参数，是二的整数次幂的变化，即随着层数的增加，由1->2->4->8->16，依次类推，至于这个模糊后续再继续详解。
第25行设置tmp层的混合模式为grain_extract， 第26行执行图层向下混合并将数据保存到new_scale中，这个时候就是相当于把Blur图像和tmp图像进行grain_extract混合，这个混合模式PS中是没有的，我们可以在GIMP的代码gimpoperationlayermode-blend.c中找到其代码：
void gimp_operation_layer_mode_blend_grain_extract (GeglOperation *operation, const gfloat *in, const gfloat *layer, gfloat *comp, gint samples) { while (samples--) { if (in[ALPHA] != 0.0f && layer[ALPHA] != 0.0f) { gint c; for (c = 0; c < 3; c++) comp[c] = in[c] - layer[c] + 0.5f; } comp[ALPHA] = layer[ALPHA]; comp += 4; layer += 4; in += 4; } }
即两者相减然后加上0.5，注意这里的数据范围都是[0,1]。
第28句就是把new_scale赋值给当前层，第32句的复制又把刚刚模糊后的数据赋值给new_scale。
当下一次循环开始的时候，新的new_scale实际上已经是上一次模糊后的值了，这个必须得到重视。
第35句则是把最后一次模糊后的值直接添加一个新的层中，并把该层命名为Residual。、
整个的过程其实就是这么简单，我们可以看到除了最后一层外，其他的层其实都是那上一次模糊后的值减去这次模糊后的值，所以他们相间后就得到了不同尺度的细节信息。
下面我们来看看这个函数中最为核心的wavelet_blur是怎么回事，在wavelet-blur.c中，并没有给出什么具体的代码实现，只有这样一段函数：
static void attach (GeglOperation *operation) { GeglNode *gegl = operation->node; GeglNode *input = gegl_node_get_input_proxy (gegl, "input"); GeglNode *output = gegl_node_get_output_proxy (gegl, "output"); GeglNode *vblur = gegl_node_new_child (gegl, "operation", "gegl:wavelet-blur-1d", "orientation", 1, NULL); GeglNode *hblur = gegl_node_new_child (gegl, "operation", "gegl:wavelet-blur-1d", "orientation", 0, NULL); gegl_node_link_many (input, hblur, vblur, output, NULL); gegl_operation_meta_redirect (operation, "radius", hblur, "radius"); gegl_operation_meta_redirect (operation, "radius", vblur, "radius"); gegl_operation_meta_watch_nodes (operation, hblur, vblur, NULL); }
这些都是一些大型软件喜欢用的东西，看的是晕头转向，核心的就是两个gegl_operation_meta_redirect调用，其实一看后面的参数也就知道了，先水平模糊，然后在垂直模糊，我们直接调转到对应的真正描述算法的部分去，即wavelet-blur-1d.c文件中。
打开wavelet-blur-1d.c文件，可以快速的看到有wav_hor_blur以及wav_ver_blur2个函数名，很明显，这个验证了我们前面的猜测。两个函数的函数体的内容基本完全相同。我们以wav_hor_blur为例：
static void wav_hor_blur (GeglBuffer *src, GeglBuffer *dst, const GeglRectangle *dst_rect, gint radius, const Babl *format) { gint x, y; GeglRectangle write_rect = {dst_rect->x, dst_rect->y, dst_rect->width, 1}; GeglRectangle read_rect = {dst_rect->x - radius, dst_rect->y, dst_rect->width + 2 * radius, 1}; gfloat *src_buf = gegl_malloc (read_rect.width * sizeof (gfloat) * 3); gfloat *dst_buf = gegl_malloc (write_rect.width * sizeof (gfloat) * 3); for (y = 0; y < dst_rect->height; y++) { gint offset = 0; read_rect.y = dst_rect->y + y; write_rect.y = dst_rect->y + y; gegl_buffer_get (src, &read_rect, 1.0, format, src_buf, GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_CLAMP); for (x = 0; x < dst_rect->width; x++) { wav_get_mean_pixel_1D (src_buf + offset, dst_buf + offset, radius); offset += 3; } gegl_buffer_set (dst, &write_rect, 0, format, dst_buf, GEGL_AUTO_ROWSTRIDE); } gegl_free (src_buf); gegl_free (dst_buf); }
整体其实没啥，前面的一堆就是分配内存，搞定读写的区域范围，核心的算法部分又是调用wav_get_mean_pixel_1D 函数，所以又只能跳转到wav_get_mean_pixel_1D 这个函数中：
static inline void wav_get_mean_pixel_1D (gfloat *src, gfloat *dst, gint radius) { gint i, offset; gdouble weights[3] = {0.25, 0.5, 0.25}; gdouble acc[3] = {0.0, }; for (i = 0; i < 3; i++) { offset = i * radius * 3; acc[0] += src[offset] * weights[i]; acc[1] += src[offset + 1] * weights[i]; acc[2] += src[offset + 2] * weights[i]; } dst[0] = acc[0]; dst[1] = acc[1]; dst[2] = acc[2]; }
这个函数就非常明朗了，要干啥也一清二楚。仔细看代码，发现原来他只是一个3个像素求加权的过程，中心点的权重占了一半，左右2个像素的权重各占1/4,。
但是注意这里的坐标偏移即offset变量，他是使用的i * radius * 3，乘以3是因为3通道，而乘以Radius则，表示他的取样并不是相邻一个像素的左右取样，而是相邻Radius个像素。
注意，因为在wav_hor_blur函数中，对src源数据指针已经有了一个向左radius的偏移，所以这里的i=0时的坐标依旧在中心点的左侧，即下述代码解决了找个问题：
GeglRectangle read_rect = {dst_rect->x - radius, dst_rect->y,
dst_rect->width + 2 * radius, 1};
说到这里，基本上这个分解的过程就已经描述完成了，下面我们来做些总结。
第一，GIMP里的这个水平和垂直方向的模糊，虽然他是水平和垂直的可分离的卷积，其实可以直接整合成一个卷积核的模糊，因为他这是3*3的，这个计算量很小，这个新的卷积和，可以用matlab计算如下：
这样有利于算法的进一步加速。
第二、前面讲的grain_extract模式的计算是in[c] - layer[c] + 0.5f; 但是注意，正在的数据应该不需要加上这个0.5f，Gimp加上这个只是为了最终显示的这个结果方便，不然这个计算结果很多是小于0的，就直接弄成黑色了。
第三、还是前面那个模糊，我们要特别注意在每次迭代的时候，虽然卷积核是一样的，但是随着层数的增加，取样的位置在越来越远离中心点，我们用下面一副图说明这个问题：
中线点是黑色的那个点，每次都参与卷积，红色的8个点是半径为1是的取样位置，绿色的8个点是半径为2时的位置，蓝色的为半径为4时的取样位置，黄色的是半径为8时的结果给，青色的是半径为16时的取样位置。注意，每次的取样图也是不一样的，这种也叫做Dilated convolution。
第四：和传统的小波分解获得的梯级结果不同（如下图所示）， GIMP这个考虑到了图层的一些显示方便，以及实际的可操作性，其生产的每层结果大小都是和原图一样的，而这个操作也是上述模糊为什么每次的半径都要扩大一倍的意思，因为原本是需要每层的大小都是上一层的一半，然后在执行半径为1的模糊，现在图层大小不变，因此就扩展取样点的位置，而不改变取样点的数量，这也是GIMP这个小波的分解的精髓所在。
我们可以在网络中找到一些使用该插件进行图像增强处理的例子，比如在https://patdavid.net/2011/12/getting-around-in-gimp-skin-retouching/这个链接中，提供了一些Skin Retouching的操作过程和结果，有兴趣的朋友可以自行试验下。
可以看到细节方面还是增强的很是细腻了。
当然，GIMP这个软件的框架太大了，他的代码更多的是实现效果，而不是考虑速度，而且GIMP也只提供了分解的过程，后续如何利用他以及如何增强需要用户自己出创作，因此，后续我还将进一步描述这个算法如何进行优化，以及如何进行一些简单的增强应用。
如果想时刻关注本人的最新文章，也可关注公众号：
