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上实现,而且速度超快。有需要的同学自己搞搞吧。