Kuwahara filter (python版)

1.前言

最近刚好在学习NPR相关的知识以及论文，发现还是要从最基础的一些算法入手学习才比较好，最后的效果虽然重要，但是算法的核心思想不清楚的话，应该也很难突破现有的研究。

2.参考文章

虽然一般参考都写在最后，但是由于有些图以及算法的基础内容本人并不想过多书写，所以请多多见谅。以下列出的参考内容都具有学习kuwahara滤镜的基础知识，本文中也会结合参考博客来讲。

1. https://imagej.nih.gov/ij/plugins/kuwahara.html
2. http://blog.csdn.net/skelking/article/details/44618963
3. http://blog.csdn.net/panda1234lee/article/details/52185530
4. http://www.cnblogs.com/Imageshop/p/6219990.html

3.kuwahara 滤镜

kuwahara滤镜经常被用在非真实感渲染当中，本人是在阅读论文时发现的，但是kuwahara的原作却没有找到，如果有原论文的话希望能够有人分享出来。

关于kuwahara滤镜的一些基础介绍请先参考第二部分列表中的第二篇。关于kuwahara滤镜的原理，简单来说就是针对卷积核内部的所有像素点二次进行分割成四个区域，找出这四个区域当中数值最平稳的，然后求一下这个区域内的平均值并将这个平均值作为卷积核中心像素点的数值。

所以说从理解上并不困难，但是细节上需要解决几个问题，那就是运算速度，如果说我们首先不考虑在GPU上实现的话，那么优化点还是有很多的，本人也是结合了以上几篇参考做了一版python的复现。

接下来，具体讲一下关于kuwahara的快速算法和基础算法。首先提出一个根据算法思路最简单的版本。根据上述基础思想可知，要求四个区域中最平稳的一个区域其实在数学上也就是求一下每个区域的方差，方差越小数据越稳定，那么也就是我们要求的目标区域。求出这个区域后对该区域进行一次平均得出最后的像素数值。由此可见，kuwahara滤镜最后的效果其实本质上就是为了平滑整个图像的细节之处，虽然在此后的l0smooth效果更好，但是kuwahara滤镜的原理却让我们能更好地感受到如何解决非真实感渲染这个问题，那就是图像平滑，本身有很多细节的图像通过平滑处理后，会发现细节内容都不见了，只有图像大致的信息保留了下来，这就达到了抽象的第一步。但是，经本人测试处理一张图像速度是实在慢，当然因为本人并未在代码上做更多优化处理，只是按照这个思路完全实现，速度慢是自然而然的，下面附上这部分代码。

def kuwahara_normal_filter(image,kernel=7):    pad_size = kernel//2    height, width, channel = image.shape    out_image = image.copy()    pad_image = np.zeros((height+pad_size*2,width+pad_size*2,channel))    for c in range(channel):        pad_image[:,:,c] = np.pad(out_image[:,:,c],[pad_size,pad_size],'constant')    for h in range(height):        for w in range(width):            for c in range(channel):                #identify the area 1,2,3,4 range                #in pad image                cur_point_index = (h + pad_size,w + pad_size,c)                area = np.zeros((4,kernel//2 +1, kernel//2 +1))                area[0] = pad_image[h:(cur_point_index[0]+1),w:(cur_point_index[1]+1),c]                area[1] = pad_image[h:(cur_point_index[0]+1),cur_point_index[1]:(cur_point_index[1]+pad_size+1),c]                area[2] = pad_image[cur_point_index[0]:(cur_point_index[0]+1+pad_size),w:(cur_point_index[1]+1),c]                area[3] = pad_image[cur_point_index[0]:(cur_point_index[0]+1+pad_size),cur_point_index[1]:(cur_point_index[1]+1+pad_size),c]                std_area = [np.std(area[0]),np.std(area[1]),np.std(area[2]),np.std(area[3])]                min_std_area_index = np.argwhere(std_area==np.min(std_area))[0,0]                out_image[h,w,c] = np.sum(area[min_std_area_index])/(len(area[min_std_area_index])**2)    return out_image

接下来就是思考如何加快这部分运算，我先参考了第二部分第二篇文章的方法，那就是通过积分图来加速这个过程。为什么用积分图呢，原因在于本身kuwahara算法核心是求和求平方求方差，所以大部分的计算直白说来就是求面积的一个过程，那么用积分图来计算是非常便捷的一种方式。

积分图：对于一幅灰度的图像，积分图像中的任意一点(x,y)的值是指从图像的左上角到这个点的所构成的矩形区域内所有的点的灰度值之和。（百度百科）

说到积分图又会牵扯到积分图的加速算法，我参考的是第二部分第四篇文章的方法。具体的细节大家可以去看参考内容。但是我在此需要声明几点，首先对于积分图的制作来说，是先给图像增加一列0，作为图像矩阵的第一列。然后再给图像增加一行0，作为图像矩阵的第一行。而后再计算积分图，关于积分图的慢速算法和快速算法我也都尝试了一下。下面直接贴出代码。

慢速算法

def cal_normal_integral_image(image):    height, width ,channel = image.shape    integral_image = np.zeros((height+1,width+1,channel))    for c in range(channel):        integral_image[:,:,c] = np.vstack((np.zeros(width+1),np.hstack((np.zeros((height,1)),image[:,:,c]))))    for h in range(height):        for w in range(width):            for c in range(channel):                integral_image[h+1,w+1,c] = integral_image[h+1,w,c]  + integral_image[h,w+1,c] - integral_image[h,w,c] + integral_image[h+1,w+1,c]    return integral_image

快速算法

def cal_fast_integral_image(image):    height, width ,channel = image.shape    integral_image = np.zeros((height+1,width+1,channel))    for c in range(channel):        integral_image[:,:,c] = np.vstack((np.zeros(width+1),np.hstack((np.zeros((height,1)),image[:,:,c]))))    for h in range(height):        sum_v = [0,0,0]        for w in range(width):            for c in range(channel):                sum_v[c] += integral_image[h+1,w+1,c]                integral_image[h+1,w+1,c] = integral_image[h,w+1,c] + sum_v[c]    return integral_image  

关于快速算法和慢速算法其实差别就在于积分图的特点，稍微参考下第二篇和第四篇参考的积分图部分的图片就可以一目了然，我就不当图片搬运工了。

最后给出关于kuwahara滤镜使用积分图部分的算法的要点，此处算法关键点在于如何求出面积，求面积的话也就是说如何确定积分图中点的坐标顶点问题，不难但是非常需要谨慎小心，下标稍微算错就会导致最终结果出错。

def kuwahara_fast_filter(image,kernel=7):    pad_size = kernel//2    height, width, channel = image.shape    out_image = image.copy()    pad_image = np.zeros((height+pad_size*2,width+pad_size*2,channel))    for c in range(channel):        pad_image[:,:,c] = np.pad(out_image[:,:,c],[pad_size,pad_size],'constant')    #calculate integral image and square of integral image    integral_image = integral_utils.cal_fast_integral_image(pad_image)    square_integral_image = integral_utils.cal_fast_integral_image(pad_image**2)    for h in range(height):        for w in range(width):            for c in range(channel):                left_top_square_num = square_integral_image[h+pad_size+1,w+pad_size+1,c] - square_integral_image[h,w+pad_size+1,c] - square_integral_image[h+pad_size+1,w,c] + square_integral_image[h,w,c]                left_top_num = integral_image[h+pad_size+1,w+pad_size+1,c] - integral_image[h,w+pad_size+1,c] - integral_image[h+pad_size+1,w,c] + integral_image[h,w,c]                left_top = left_top_square_num - (left_top_num**2)/((pad_size+1)**2)                right_top_square_num = square_integral_image[h+pad_size+1,w+kernel,c] - square_integral_image[h,w+kernel,c] - square_integral_image[h+pad_size+1,w+kernel-pad_size-1,c] + square_integral_image[h,w+kernel-pad_size-1,c]                right_top_num = integral_image[h+pad_size+1,w+kernel,c] - integral_image[h,w+kernel,c] - integral_image[h+pad_size+1,w+kernel-pad_size-1,c] + integral_image[h,w+kernel-pad_size-1,c]                right_top = right_top_square_num - (right_top_num**2)/((pad_size+1)**2)                left_bottom_square_num = square_integral_image[h+kernel,w+pad_size+1,c] - square_integral_image[h+kernel,w,c] - square_integral_image[h+kernel-pad_size-1,w+pad_size+1,c] + square_integral_image[h+kernel-pad_size-1,w,c]                left_bottom_num = integral_image[h+kernel,w+pad_size+1,c] - integral_image[h+kernel,w,c] - integral_image[h+kernel-pad_size-1,w+pad_size+1,c] + integral_image[h+kernel-pad_size-1,w,c]                left_bottom = left_bottom_square_num - (left_bottom_num**2)/((pad_size+1)**2)                right_bottom_square_num = square_integral_image[h+kernel,w+kernel,c] - square_integral_image[h+kernel-pad_size-1,w+kernel,c] - square_integral_image[h+kernel,w+kernel-pad_size-1,c] + square_integral_image[h+kernel-pad_size-1,w+kernel-pad_size-1,c]                right_bottom_num = integral_image[h+kernel,w+kernel,c] - integral_image[h+kernel-pad_size-1,w+kernel,c] - integral_image[h+kernel,w+kernel-pad_size-1,c] + integral_image[h+kernel-pad_size-1,w+kernel-pad_size-1,c]                right_bottom = right_bottom_square_num - (right_bottom_num**2)/((pad_size+1)**2)                num_area = [left_top_num,right_top_num,left_bottom_num,right_bottom_num]                        std_area = [left_top,right_top,left_bottom,right_bottom]                min_std_area_index = np.argwhere(std_area==np.min(std_area))[0,0]                out_image[h,w,c] = num_area[min_std_area_index]/((pad_size+1)**2)    return out_image 

最后的最后给出实验结果：慢速kuwahara滤镜需要64.55秒计算一张400*500的图片，而快速只需要17秒左右。由此可见运算速度差别还是很大的，但是这个运算速度仍然很慢，因此GPU版本的封装是最好的方式，基本可以达到实时计算，这个以后有机会再尝试并分享出来。

4.结语

关于kuwahara滤镜的计算过程当中由于本人能力有限，可能存在不适当的代码，希望大家可以多给意见改善代码。最后附上通过kuwahara滤镜的测试图，以及代码下载的地址。

下一篇是kuwahara算法的改进算法，Generalized Kuwahara