SSE图像算法优化系列三十：GIMP中的Noise Reduction算法原理及快速实现。
GIMP源代码链接：https://gitlab.gnome.org/GNOME/gimp/-/archive/master/gimp-master.zip
GEGL相关代码链接：https://gitlab.gnome.org/GNOME/gegl/-/archive/master/gegl-master.zip
最近因为要研究下色温算法，顺便下载了最新的GIMP软件，色温算法倒是找到了（有空单独来讲下），也顺便看看GIMP都有些什么更新，嗯，更新还是蛮多的，界面UI上有很多改动，有些已经改的面目全非了。随便瞄了一下Enhance菜单，发现里面有一个Nosie Reduction算法，试了下，还有点效果。于是在github上下载了GIMP的源代码，可是在源代码里搜索相关的关键词确没有发现任何的相关代码，后来才发现很多东西都有个GEGL关键词，结果一百度，原来他是一个单独的软件包，于是有下载了GEGL的源代码，终于在gegl-master\operations\common\里面看到了noise-reduction.c文件。
其核心的代码如下：
复制代码
static void
noise_reduction (float *src_buf, /* source buffer, one pixel to the left
and up from the starting pixel */
int src_stride, /* stridewidth of buffer in pixels */
float *dst_buf, /* destination buffer */
int dst_width, /* width to render */
int dst_height, /* height to render */
int dst_stride) /* stride of target buffer */
{
int c;
int x,y;
int dst_offset;
#define NEIGHBOURS 8
#define AXES (NEIGHBOURS/2)
#define POW2(a) ((a)*(a))
/* core code/formulas to be tweaked for the tuning the implementation */
#define GEN_METRIC(before, center, after) \
POW2((center) * 2 - (before) - (after))
/* Condition used to bail diffusion from a direction */
#define BAIL_CONDITION(new,original) ((new) > (original))
#define SYMMETRY(a) (NEIGHBOURS - (a) - 1) /* point-symmetric neighbour pixel */
#define O(u,v) (((u)+((v) * src_stride)) * 4)
int offsets[NEIGHBOURS] = { /* array of the relative distance i float
* pointers to each of neighbours
* in source buffer, allows quick referencing.
*/
O( -1, -1), O(0, -1), O(1, -1),
O( -1, 0), O(1, 0),
O( -1, 1), O(0, 1), O(1, 1)};
#undef O
dst_offset = 0;
for (y=0; y> 3;
}
free(Temp);
return IM_STATUS_OK;
}
复制代码
是不是看起来比上面的GIMP得要舒服些，而且中间也大概只要原始图像2倍的一个临时内存了。在速度和内存占用方面都前进了很多。
我测试前面提到的那副3K*2K的图像，耗时要7S多，但是我测试表面GIMP用了多核的，如果论单核，我这里的速度要比他快2倍多。
很明显，这个速度是不可以接受的，我们需要继续优化。
我还是老套路，使用SIMD指令做处理，看到上面的代码，其实真的觉得好容易改成SIMD的。
short LT_RB = LT + RB, RT_LB = RT + LB;
short T_B = T + B, L_R = L + R, C_C = C + C;
short Dist1 = IM_Abs(C_C - LT_RB), Dist2 = IM_Abs(C_C - T_B);
short Dist3 = IM_Abs(C_C - RT_LB), Dist4 = IM_Abs(C_C - L_R);
这些加减绝对值都有完全对应的SSE指令。 _mm_add_epi16、 _mm_sub_epi16、_mm_abs_epi16，基本上就是照着写。
稍微复杂一点就是这里：
if ((IM_Abs(LT_C - LT_RB) < Dist1) && (IM_Abs(LT_C - T_B) < Dist2) && (IM_Abs(LT_C - RT_LB) < Dist3) && (IM_Abs(LT_C - L_R) < Dist4))
{
Sum += LT_C;
Amount += 2;
}
在C语言里，这里判断会进行短路计算，即如果前一个条件已经不满足了，后续的计算就不会进行。但是在SIMD指令里，是没有这样的机制的。我们只能全部计算，然后在通过某一种条件组合。
在合理，要实现符合条件就进行累加，不符合条件就不做处理的需求，我们需要稍作修改，即不符合条件不是不做处理，而是加0，加0对结果没有影响的。主要借助下面的_mm_blendv_epi8来实现。
复制代码
__m128i LT_C = _mm_add_epi16(LT, C);
Flag1 = _mm_cmplt_epi16(_mm_abs_epi16(_mm_sub_epi16(LT_C, LT_RB)), Dist1); // 只能全部都计算，但还是能提速
Flag2 = _mm_cmplt_epi16(_mm_abs_epi16(_mm_sub_epi16(LT_C, T_B)), Dist2);
Flag3 = _mm_cmplt_epi16(_mm_abs_epi16(_mm_sub_epi16(LT_C, RT_LB)), Dist3);
Flag4 = _mm_cmplt_epi16(_mm_abs_epi16(_mm_sub_epi16(LT_C, L_R)), Dist4);
Flag = _mm_and_si128(_mm_and_si128(Flag1, Flag2), _mm_and_si128(Flag3, Flag4));
Sum = _mm_adds_epu16(Sum, _mm_blendv_epi8(Zero, LT_C, Flag));
Amount = _mm_adds_epu16(Amount, _mm_blendv_epi8(Zero, Two, Flag));
复制代码
注意到我们这里用到了_mm_adds_epu16，无符号的16位加法，这是因为我需要尽量的提速，因此需要减少类型转换的次数。同时，我们看到在统计累加值时，我们并没有求平均值，而是直接用的累加值，这样理论上最大的累加值就是 255 * n * (8 + 1) * 2 < 65535， 这样n最大能取15，但是15不是个好数据，在放大和缩小都不能用移位来实现，所以我最后取得放大系数为8。
另外，在最后还有个16位的整数的除法问题，这个没有办法，SSE