Generalized Kuwahara (python版)

1.前言

此篇文章是基于第一篇Kuwahara filter基础上，针对该算法的改进。当然并不是我自己改进，主要还是根据参考论文算法改进而来的，因为本身并没有作者的源代码，所以有可能在实现上会出错，所以还请大家多多指出错误并且给予一定的改进方案。

2.参考论文

仍然先列出来参考论文，方便大家查找。

1. M. Nagao and T. Matsuyama, ”Edge preserving smoothing”, Computer Graphics and Image Processing, vol. 9, no. 4, pp. 394–407, Apr. 1979.
2. van den Boomgaard R.: Decomposition of the Kuwahara-Nagao operator in terms of linear smoothing and mor- phological sharpening. In ISMM 2002 (2002). 2
3. Papari G., Petkov N., Campisi P.: Artistic edge and corner enhancing smoothing. IEEE Transactions on Image Pro- cessing 16, 10 (2007), 2449–2462

3.Kuwahara滤镜改进算法

如果看过我第一篇以及我第一篇引用的博客内容的朋友应该对Kuwahara滤镜的计算方式有一个基础的认识。核心的计算就是在卷积核内部找出标准差最小的一个区域然后作为目标像素点的数值。那么后人基于这个思想首先做的改进就是子窗口形状的修改，参考论文1就做了这样一件事情。就是把原先卷积核内部分成了四个长方形区域而现在则是为了从多个方向找出最平稳的数值进行最后的赋值，即子窗口的内容由长方形变成了扇形。如下图所示。

此图来源于第三篇参考文献。

在第一篇论文关于分割出子区域修改的基础上第二篇论文将高斯核引入了进来，就使得原本直接对目标区域求均值变成了和高斯核进行进行一次卷积操作。如下为第二篇参考文章的高斯核公式。

第三篇参考论文的优化则是做的最后一步，也就是如何计算最终目标像素的数值。

首先给出每个区域内部平均值和方差的计算，此处优化为加权平均和加权标准差的计算。如下公式所示。

第一个公式当中计算了原图与高斯核的卷积，即加权平均运算，因为高斯核部分约束条件是高斯核所有数的总和为1.

第二个公式则是计算加权标准差，将原图像素值平方再次与正规化后的高斯核进行卷积随后减去第一个公式中计算出的加权平均数值，得到了每个区域的加权标准差的数值。

下面一步，也就是最后的一步计算，也最为重要。

这个公式我自己的理解是，先从整体结构上分析一下，我们不考虑每个参数的具体意义，只考虑整个公式。就会发现这个貌似仍然是计算加权平均的公式，上面是数值和权重乘积，下面是权重的和。如此看来，这个就有以下理解。

我们将一个卷积核区分为了八个子区域，每个子区域的加权平均值已经算出来了，我们现在不用以前通过寻找最小标准差的方式来计算目标像素值，而是改为利用八个区域的加权平均数在进行一次加权平均作为最后的数值。那么问题来了，权重是怎么确定的呢。按照一般思路，每个区域的标准差越大也就是越不平稳，权重要降低，那么看到公式就明白权重就是上面计算出的标准差，但是直接相乘起来会导致与刚才的推算恰好相反。所以我们给刚算出的标准差一个指数，并且是一个负的指数，这就使得只要这个指数越负，那么原先标准差越小的区域就会比标准差大的区域数值大。就如同，2的3次方是8比3的3次方27小，但是2的-3次方是1/8比3的-3次方1/27大的道理一样。非常简单，但是数学含义确确实实是正确的。

下面给出一段核心代码，此处不妥的地方请大家多多包涵，并且给出建议。同时，因为划分区域已经不是长方形，那么也就不能像第一篇讲的用积分图来加速运算。此处其实需要打个问号，是否可以用类似的方法来进行快速计算呢。

def generalized_kuwahara_filter(image,kernel = 7,q = 3,div = 8):
    pad_size = kernel//2
    gauss_kernel = gaussian_kernel(kernel)
    
    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)
                div_res = divide_kernel(pad_image[:,:,c],kernel,cur_point_index,div)
                
                
                total_ms = 0
                total_s = 0
                for i in range(div):
                    cur_sum = 0
                    cur_var = 0
                    cur_weight = 0
                    cur_index = div_res[i]
                    for ci in cur_index:
                        g = gauss_kernel[ci[0]-cur_point_index[0]+pad_size,ci[1]-cur_point_index[1]+pad_size]
                        cur_val = pad_image[ci[0],ci[1],c] 
                        cur_sum += cur_val*g
                        cur_var += (cur_val**2)*g
                        cur_weight += g
                    
                    if cur_weight !=0:
                        cur_sum /= cur_weight
                        cur_var /= cur_weight
                    
                    cur_std = cur_var - cur_sum**2
                    
                    if cur_std > 1e-10:
                        cur_std = np.sqrt(cur_std)
                    else:
                        cur_std = 1e-10
                    
                    total_ms += cur_sum * np.power(cur_std,-q)
                    total_s  += np.power(cur_std,-q)
                
                res_val = 0
                if total_s > 1e-10:
                    res_val = total_ms/total_s
                    if res_val>1:
                        out_image[h,w,c] = res_val
    return out_image

4.实验结果

此次实验结果如下图所示，第一张是用kuwahara滤镜做出的结果，第二张是用参考论文3做出的结果，明显的看来在边缘以及各个区域的效果变得更加平滑了，说明此次实验成功。但是由于并没有加速计算，所以计算一张400*500的图需要花145秒左右，相当耗时，需要有加速算法才可以，这是现在并未考虑GPU加速的情况。

最后给出实验的测试代码，下载代码

----------------------------------------20181019更新-----------------------

附上WebGL的demo

----------------------------------------20190508更新-----------------------

附上WebGL版代码的github地址https://github.com/RaymondMcGuire/GPU-Based-Image-Processing-Tools