Combining Sketch and Tone for Pencil Drawing Production的优化过程

Combining Sketch and Tone for Pencil Drawing Production 是一个非常出色的自然图像转铅笔画的算法，具体的算法原理就不多说了，不了解的可以看一下这个博客，写得非常清楚了：
http://blog.csdn.net/bluecol/article/details/45422763
这里主要写一下我的开发过程，和用到的优化方法。
这个算法主要有三个步骤：
（1） Generate Structure Map
（2） Generate Tone Map
（3） Generate Texture Map
这里面，最费时间的是第（3）步，其中涉及到了一个超大规模稀疏线性方程组的求解，其将是第（1）步，涉及到了图像的16次卷积操作。第（2）费时相当短，不需要优化。

先看第（3）步，按照文章的实现思路，构造稀疏矩阵，然后用CG求解，对于450 * 600的图像，需要花费大约600ms，而且由于直接采用了现成的稀疏矩阵算法库，这个模块优化起来相当吃力。那有没有方法可以替代这个步骤呢？我们先观察一下生成所CG求解出来的beta有什么样的特征。用下现这两幅图作实验，其中左图是原图，右图是纹理图。

用CG求解出来的beta是这个样子：

可以在beta图上看到明显的输入的纹理图像的痕迹，因此，可以猜测，beta本质上是由输入到的纹理图像与另外一幅图合成的。再考虑到Texture map 的生成算法：
Texture Map = (Tone Map) ^beta，也就是说texture Map由于Tone Map经过一个gamma变换的衰减得到的，因此，这个“另外一幅图”肯定是和Tone Map相关的。上面右图就是Tone Map图，比较beta与Tone Map，可以很容易的得到用于合成beta的另外一幅图就是：（1 – Tone Map）：

找到了这两个用于合成的图，我们就可以借用铅笔画生成算法里将Structure Map 和 Texture Map进行合成时所采用的方法，将这两个图合成一个beta图，具体方法就是将这两个图的对应像素进行相乘，可以得到下面的beta图的一个近似，我们称为beta_:

对比，beta和beta_，可以看出，这两个非常相似，用于近似计算是可以的，下面是用这两个图分别得到的最终铅笔画，差别也很小。但是还是有肉眼可分辨的差别，原始的beta生成的结果，其纹理更加细腻，几乎看不出输入的铅笔纹理，都变成了类似一种小颗粒样的纹理。但是用beta_生成的结果，还是可以看出来明显的纹理图像的铅笔笔画痕迹的。

经过这样的处理，我们将最难优化的部分进行了简化，不用再求解超大规模稀疏线性方程组，而只用图像最基本的像素加减乘法，因此速度变得非常快，第（3）步的优化完成。当然如果想追求高品质的效果，这样的简化是不合适的。但如果想做一个可以实时处理的高速程序，这样的处理有时候不可避免。

下面看第（1）的优化。前面说过，这里最费时的是两组各8个方向的卷积运算，最原始的卷积运算的写法是下面这样：

img_y = img;    
for (int y = 0; y < height; y++, img_y += width){
    float* img_x = img_y;
    for (int x = 0; x < width; x++, img_x++){
        unsigned char* kn_y = kernel;
        conv_y = m_conv2DImg + y * m_convWidth + x;
        float sum = 0;
        for (int yy = 0; yy < kerSize; yy++, kn_y += kerSize, conv_y += m_convWidth){
            unsigned char* kn_x = kn_y;
            float* conv_x = conv_y;
            for (int xx = 0; xx < kerSize; xx++, kn_x++, conv_x++){
                sum += *conv_x * *kn_x;
            }
        }
        *img_x = sum;
    }
}

这个运算有4重for循环，可想而知其效率。经测试，采用上面的玫瑰花的图片，在9*9的卷积核下，需要180ms计算8个方向的卷积，而当卷积核增加到21*21时，时间则飙升到超过1s，这是完全无法接受的。

优化的关键有两点，一是在想办法减少内层循环的数量，二是去掉乘法运算。

先说去掉乘法运算。网上能查到的关于这个算法的实现，无论是C++的还是Matlab的，在生成卷积核时基本都采用了先生成水平方向的卷积核，然后再旋转卷积核得到其他7个。这样做最大的问题是，在旋转卷积核的操作中，有插值算法的存在，导致了卷积核中出现了非0非1的值，这样我们只能采用相乘再相加的原始方法计算卷积。其实我们可以将旋转生成的卷积核再进行一次阈值操作，或者直接自己在8个方向上画直线生成8个值为非0即1的卷积核。这样做的好处非常明显，我们可以只针对卷积核中值为1的像素进行相加操作，从而去掉了乘法。

再说减少内层循环。由于卷积核中大部分像素都是0，因此在最开始的时候，我们可以将这个卷积核放大图像的左上角，然后用一个指针数组记录下这个卷积核中非0像素对应的原图像指针，卷积操作就是将这个数组中的指针对应的值进行相加，然后再除以非0值的个数（这个数可以提前计算，然后将除法变成乘法）。卷积核每向右移动一下，这个指针数组里的所有指针都向右移动一下（+1）。这样我们就把内层的双层循环变成了单层循环。经过这样的处理，对于9 * 9的卷积核，速度几乎可以提高6 - 8倍，由180ms变成 20多ms。经过这步优化之后，代码基本上就变成了下面这样：

img_y = img;    
    for(int y = 0; y < height; y ++, img_y += width){
        float* img_x = img_y;
        for(int x = 0; x < width; x ++, img_x++){
            float sum = 0;
            for(int i = 0; i < listLen; i ++){//listLen就是指针数组的长度
                sum += *(ptmp[i]);
                ptmp[i] ++;
            }
            *img_x = sum;
        }
        for(int i = 0; i < listLen; i ++){
            ptmp[i] += kerSize;
        }
}

优化到这里之后，还能不能继续优化？毕竟内层还是有一个for循环，还是比较费时的，能不能把这个for也去掉呢？答案是肯定的。
这主要得益于这个8个方向的卷积核的非0元素都是一条直线，这样我们就可以借用积分图的思想，也就是沿直线方向计算一个一维的积分图，然后用这个一维的积分图来计算任意区间的像素之和。这样做就可以用两次遍历来完成对 x < width 和 i < listLen 的两层循环。
另外一个效率差不多的优化方法也是借用直线的性质，以水平方向的卷积核为例：

如上图，假设卷积核长度为5，我们可以按照正常方法先计算出每一行最初5个像素的和，记为sum, 卷积核向右移动一个像素之后，我们看到只有1移出来了，而6新加进来，中间的2，3，4，5都没有变化，因此上面蓝框里像素的和就可以用sum’ = sum – P_1 + P_6 来计算。利用这个方法也可以达到去除内层循环的目的。对于水平方向，采用这个方法优化之后，程序一般是下面这样子，可以看到，基本只剩下两层for了。有意思的是，采用这个方法优化之后，对于9*9大小的核，如果核内非0元素的宽度是1，则8个方向的卷积总费时与上面第一次优化之后的结果差不多，还是20ms左右。单看每个方向上的卷积时间，水平，竖直，与45度，135度这四个方向上计算速度有明显加快，但是22.5度，67.5度等四个方向计算是有很多中间量需要计算，比上面的算法略慢。但是如果采用更大的核，如21*21,则加速很明显。而且采用这种方法的好处是，卷积的计算速度与卷积核的大小没有关系，只与图像的大小相关了，速度变得很稳定。例如，采用21*21的核之后，8个方向的卷积时间只比9*9的核多了不到1ms

for(int y = 0; y < height; y++, img_y += width, conv_y += m_convWidth){
      float* conv_x = conv_y + half_size;
      float* img_x = img_y;
      //calc the first sum
      for(int i = 0; i < listLen; i++){
          ptmp[i] = conv_x + offset[i];
      }
      float sum = 0;
      for(int i = 0; i < listLen; i++){
          sum += *(ptmp[i]);
      }
      img_x[0] = sum;
      startp[0] = conv_x + offset[0];
      endp[0] = conv_x + offset[listLen - 1];
      for(int x = 1; x < width; x++){
          sum -= *(startp[0]);
          startp[0]++;
          endp[0]++;
          sum += *(endp[0]);
          img_x[x] = sum;
      }
  }

到目前为止，我们已经基本上完成了两个模块的优化，对于640*480的图，已经基本上达到实时处理了。但所谓优化无止境，如果想进一步加速的话，下面这些方法还可以偿试：
1. 8个方向的卷积采用多线程
2. 对于逐像素处理的程序块采用SSE（手机上可以用Neon）并行处理。
3. 其实在去掉求解稀疏线性方程组的模块之后，整个系统只剩下最基本的像素加减乘和卷积运算，这三种运算都可以在GPU上实现，而且速度超快。有需要的同学自己搞搞吧。