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<dst_height; y++)
    {
      float *center_pix = src_buf + ((y+1) * src_stride + 1) * 4;
      dst_offset = dst_stride * y;
      for (x=0; x<dst_width; x++)
        {
          for (c=0; c<3; c++) /* do each color component individually */
            {
              float  metric_reference[AXES];
              int    axis;
              int    direction;
              float  sum;
              int    count;

              for (axis = 0; axis < AXES; axis++)
                { /* initialize original metrics for the horizontal, vertical
                     and 2 diagonal metrics */
                  float *before_pix  = center_pix + offsets[axis];
                  float *after_pix   = center_pix + offsets[SYMMETRY(axis)];

                  metric_reference[axis] =
                    GEN_METRIC (before_pix[c], center_pix[c], after_pix[c]);
                }

              sum   = center_pix[c];
              count = 1;

              /* try smearing in data from all neighbours */
              for (direction = 0; direction < NEIGHBOURS; direction++)
                {
                  float *pix   = center_pix + offsets[direction];
                  float  value = pix[c] * 0.5 + center_pix[c] * 0.5;
                  int    axis;
                  int    valid;

                  /* check if the non-smoothing operating check is true if
                   * smearing from this direction for any of the axes */
                  valid = 1; /* assume it will be valid */
                  for (axis = 0; axis < AXES; axis++)
                    {
                      float *before_pix = center_pix + offsets[axis];
                      float *after_pix  = center_pix + offsets[SYMMETRY(axis)];
                      float  metric_new =
                             GEN_METRIC (before_pix[c], value, after_pix[c]);

                      if (BAIL_CONDITION(metric_new, metric_reference[axis]))
                        {
                          valid = 0; /* mark as not a valid smoothing, and .. */
                          break;     /* .. break out of loop */
                        }
                    }
                  if (valid) /* we were still smooth in all axes */
                    {        /* add up contribution to final result  */
                      sum += value;
                      count ++;
                    }
                }
              dst_buf[dst_offset*4+c] = sum / count;
            }
          dst_buf[dst_offset*4+3] = center_pix[3]; /* copy alpha unmodified */
          dst_offset++;
          center_pix += 4;
        }
    }
}

   这个代码看上去比较混乱,没办法,大型软件没有哪一个代码看上去能让人省心的,而且也不怎么讲究效率,我测试了一个3K*2K的彩色图,在GIMP里大概在4S左右处理完成,属于很慢的了,看这个代码,也知道大概有width * height * 3 * 8 * 4 * Iter个循环,计算量确实是相当的大。

  我试着尝试优化这个算法。

  优化的第一步是弄明白算法的原理,在GIMP的UI界面上可当鼠标停留在Noise Reduction菜单上时,会出现Anisotroic smoothing operation字样,所以初步分析他是属于各项异性扩散类的算法。稍微分析下代码,也确实是。明显这属于一个领域滤波器,对每一个像素,求取其3*3领域内的累加值,但是3*3领域的权重并不是平均分布或者是高斯分布,而是和领域的值有关的,如果领域的值是一个边缘点,他将不参与到累加中,权重为0,否则权重为1。

  具体一点,对于领域里的任何一点,我们先求取其和中心点的平均值,对应 float value = pix[c] * 0.5 + center_pix[c] * 0.5; 这条语句,然后计算这个值在2个45度对角线及水平和垂直方向的梯度(4个梯度值)是否比中心点在四个方向的梯度都小,如果都小,说明这个领域点不属于边缘点,可以往这个方向扩散,把他计入到统计值中,如果有任何一个方向小了,则不参与最终的计算。

  上面的过程可以看成是标准的各项异性扩散的特殊在特殊处理,他具有各项异性扩散的特性,也具有一些特殊性。

  下一步,稍微分析下最简单的优化方法。第一,我们知道,在大循环里一般不建议嵌套入小的循环,这样是很低效的。我们观察到上面代码里的

      for (c=0; c<3; c++) /* do each color component individually */

  这个语句主要是为了方便表达3通道的处理的方便,但是其实三通道之间的处理时没有任何联系的,对于这样的算法,很明显,我们可以一次性当然处理R G B R G B R G B ,而不需要像GIMP这个代码这样按照 RRR  GGG  BBB这样的顺序来写,GIMP这种写法浪费了很多CPU的CACHE,毕竟R和G和B在内存类分布本来就是连续的。这样就减少了一个小循环。

  第二个优化的点是,对于普通的图像数据,我们可以考虑不用浮点数来处理,毕竟上述计算里只有*0.5这样的浮点操作,我们考虑将原先的图像数据放大一定的倍数,然后用整形来玩,在处理完后,在缩小到原来的范围,比如使用short类型应该就足够了,我把数据放大16倍或者32倍,甚至8倍应该都能获得足够的精度。

  第三个优化点,程序中是使用的Pow来判断梯度的大小的,其实可以不用,直接使用绝对值的结果和Pow是完全一样的,而绝对值的计算量比pow要小很多,对于整数则更为如此(还可以不考虑pow数据类型的改变,比如short的绝对值还是short类型,但是其pow可能就需要用int来表示了,这在SIMD优化会产生不同的结果)。

  第四个优化点是 for (axis = 0; axis < AXES; axis++)这个小循环我们应该把它直接展开。

  第五点,我们还可以考虑我在其他文章里提到的支持Inplace操作的方式,这样noise_reduction这个函数的输入和输出就可以是同一个内存。

  第六点还有小点上的算法改进,比如一些中间计算没必要重复进行,有些可以提到外部来等。

  综合上面的描述,我整理除了一个优化的C语言版本的程序,如下所示:

void IM_AnisotropicDiffusion3X3(short *Src, short *Dest, int Width, int Height, int SplitPos, int Stride)
{
    int Channel = Stride / Width;

    short *RowCopy = (short *)malloc((Width + 2) * 3 * Channel * sizeof(short));
    short *First = RowCopy;
    short *Second = RowCopy + (Width + 2) * Channel;
    short *Third = RowCopy + (Width + 2) * 2 * Channel;
    memcpy(Second, Src, Channel * sizeof(short));
    memcpy(Second + Channel, Src, Width * Channel * sizeof(short));                                                    //    拷贝数据到中间位置
    memcpy(Second + (Width + 1) * Channel, Src + (Width - 1) * Channel, Channel * sizeof(short));

    memcpy(First, Second, (Width + 2) * Channel * sizeof(short));                                                    //    第一行和第二行一样

    memcpy(Third, Src + Stride, Channel * sizeof(short));                                                            //    拷贝第二行数据
    memcpy(Third + Channel, Src + Stride, Width * Channel* sizeof(short));
    memcpy(Third + (Width + 1) * Channel, Src + Stride + (Width - 1) * Channel, Channel* sizeof(short));

    for (int Y = 0; Y < Height; Y++)
    {
        short *LinePD = Dest + Y * Stride;
        if (Y != 0)
        {
            short *Temp = First; First = Second; Second = Third; Third = Temp;
        }
        if (Y == Height - 1)
        {
            memcpy(Third, Second, (Width + 2) * Channel * sizeof(short));
        }
        else
        {
            memcpy(Third, Src + (Y + 1) * Stride, Channel * sizeof(short));
            memcpy(Third + Channel, Src + (Y + 1) * Stride, Width * Channel * sizeof(short));                            //    由于备份了前面一行的数据,这里即使Src和Dest相同也是没有问题的
            memcpy(Third + (Width + 1) * Channel, Src + (Y + 1) * Stride + (Width - 1) * Channel, Channel * sizeof(short));
        }
        for (int X = 0; X < SplitPos * Channel; X++)
        {
            short LT = First[X], T = First[X + Channel], RT = First[X + 2 * Channel];
            short L = Second[X], C = Second[X + Channel], R = Second[X + 2 * Channel];
            short LB = Third[X], B = Third[X + Channel], RB = Third[X + 2 * Channel];
            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);
            
            int Sum = C_C, Amount = 2;

            short LT_C = LT + C;
            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;
            }
            short T_C = T + C;
            if ((IM_Abs(T_C - LT_RB) < Dist1) && (IM_Abs(T_C - T_B) < Dist2) && (IM_Abs(T_C - RT_LB) < Dist3) && (IM_Abs(T_C - L_R) < Dist4))
            {
                Sum += T_C;
                Amount += 2;
            }
            short RT_C = RT + C;
            if ((IM_Abs(RT_C - LT_RB) < Dist1) && (IM_Abs(RT_C - T_B) < Dist2) && (IM_Abs(RT_C - RT_LB) < Dist3) && (IM_Abs(RT_C - L_R) < Dist4))
            {
                Sum += RT_C;
                Amount += 2;
            }
            short L_C = L + C;
            if ((IM_Abs(L_C - LT_RB) < Dist1) && (IM_Abs(L_C - T_B) < Dist2) && (IM_Abs(L_C - RT_LB) < Dist3) && (IM_Abs(L_C - L_R) < Dist4))
            {
                Sum += L_C;
                Amount += 2;
            }
            short R_C = R + C;
            if ((IM_Abs(R_C - LT_RB) < Dist1) && (IM_Abs(R_C - T_B) < Dist2) && (IM_Abs(R_C - RT_LB) < Dist3) && (IM_Abs(R_C - L_R) < Dist4))
            {
                Sum += R_C;
                Amount += 2;
            }
            short LB_C = LB + C;
            if ((IM_Abs(LB_C - LT_RB) < Dist1) && (IM_Abs(LB_C - T_B) < Dist2) && (IM_Abs(LB_C - RT_LB) < Dist3) && (IM_Abs(LB_C - L_R) < Dist4))
            {
                Sum += LB_C;
                Amount += 2;
            }
            short B_C = B + C;
            if ((IM_Abs(B_C - LT_RB) < Dist1) && (IM_Abs(B_C - T_B) < Dist2) && (IM_Abs(B_C - RT_LB) < Dist3) && (IM_Abs(B_C - L_R) < Dist4))
            {
                Sum += B_C;
                Amount += 2;
            }
            short RB_C = RB + C;
            if ((IM_Abs(RB_C - LT_RB) < Dist1) && (IM_Abs(RB_C - T_B) < Dist2) && (IM_Abs(RB_C - RT_LB) < Dist3) && (IM_Abs(RB_C - L_R) < Dist4))
            {
                Sum += RB_C;
                Amount += 2;
            }
            LinePD[X] = Sum / Amount;
        }
    }
    free(RowCopy);
}

  调用函数

int IM_ReduceNoise(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int SplitPos,  int Strength)
{
    int Channel = Stride / Width;
    if ((Src == NULL) || (Dest == NULL))                        return IM_STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))                            return IM_STATUS_INVALIDPARAMETER;
    if ((Channel != 1) && (Channel != 3))                         return IM_STATUS_INVALIDPARAMETER;

    Strength = IM_ClampI(Strength, 1, 10);
    SplitPos = IM_ClampI(SplitPos, 0, Width);
    int Status = IM_STATUS_OK;
    short *Temp = (short *)malloc(Height * Stride * sizeof(short));
    if (Temp == NULL)    return IM_STATUS_OUTOFMEMORY;
    for (int Y = 0; Y < Height * Stride; Y++)
    {
        Temp[Y] = Src[Y] << 3;
    }
    for (int Y = 0; Y < Strength; Y++)
    {
        IM_AnisotropicDiffusion3X3(Temp, Temp, Width, Height, SplitPos, Stride);
    }
    for (int Y = 0; Y < Height * Stride; Y++)
    {
        Dest[Y] = Temp[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指令没有提供整数的除法计算方法,还只能转换到浮点后,再次转换回来。

  这样用SSE处理后,还是同一幅测试图像,在同一台PC上速度能提升到400ms(4次迭代),比之前的普通的C语言提高了约17倍的速度。

  在现代CPU中,具有AVX2指令集已经是很普遍的了,单AVX2能同时处理256字节的数据,比SSE还要多一倍,我也常使用AVX2进行优化处理,速度能达到250ms,相当于普通C语言的28倍之多(但是AVX编程里有很多坑,这些坑都拜AVX不是完全的按照SSE的线性扩展导致的,这个后续有时间我单独提出)。

  经过测试,1080P的图像使用4次迭代大约需要80ms,3次迭代55ms,2次迭代月40ms,也就是说前面的一些方法和缩小所使用的时间几乎可以忽略。

  选了几幅有特点的图进行了去燥测试,其中分界线左侧的位处理的效果,右侧为未处理的。

 

  但是,这个算法也还是不是很好,他对于图像容易出现轻微的油画效果,对于一些细节特别丰富的图像非常明显,比如下图:

  这个应该是不太可以接受的,也许可以通过修改部分权重的规则来改变这个现象。这个属于后期研究的问题了。

     另外,在GIMP里也提供了这个算法的OPENCL实现,有兴趣的可以源代码里找一找,不晓得速度怎么样。

  本文Demo下载地址:  https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar,见其中的Denoise -> Anisotroic Diffusion 菜单。

  写博不易,欢迎土豪打赏赞助。

 

 

posted @ 2019-11-18 08:45  Imageshop  阅读(4397)  评论(3编辑  收藏  举报