Shader特效——“Kuwahara Filter”的实现 【OpenCV】【GLSL】

本文代码参考自 skelking 的一篇博客 《kuwahara filter 实现》

Kuwahara 是一种降噪低通滤波器，能够较好的保留物体的边缘。

基本思想就是 将Kuwahara 的模板以目标像素为中心分成4块邻域，然后分别计算四块邻域的方差，取方差最小的邻域计算其平均值，得到的结果作为目标像素的新值。

还有Java版本的代码：

Kuwahara Filter 作者的官网：https://imagej.nih.gov/ij/plugins/kuwahara.html

原理很简单，博客里也写的很清楚，我就不赘述了，在此我把代码结构整理得更清晰一点，并删除了一些多余的代码。

void procKuwaharaFilter(const Mat& src, Mat& dst, int kernel, int width, int height)
{
assert(src.channels() == 1
&& dst.channels() == 1);

//Integral Box 二维指针数组
int size = 4096;
double** box_sum;
double** box_sq_sum;

box_sum = new2DDoublePointer(size);
box_sq_sum = new2DDoublePointer(size);

// ----------------------------

if (kernel < 3 || kernel > 15)
{
throw runtime_error("the kernel is out of range");
}
if (kernel / 2 == 0)
{
throw runtime_error("the kernel must be odd");
}

// ---------------------------------------------------

calcIntegralBox(src, box_sum, width, height);
calcSquareIntegralBox(src, box_sq_sum, width, height);

// ---------------------------------------------------
//kuwahara滤波器
int neibourhood = kernel / 2 + 1;
float N = 1.f * neibourhood * neibourhood;

for (int h = neibourhood; h < height - neibourhood; h++)
{
for (int w = neibourhood; w < width - neibourhood; w++)
{
// 根据 Integral Box 来计算 n*n 目标区域的方差
//左上
Point lt(w - neibourhood , h - neibourhood);
Point lb(w - neibourhood , h );
Point rt(w, h - neibourhood);
Point rb(w, h );

double lt_sq_sum = calcNeibourhood(box_sq_sum, lt, rt, rb, lb);
double lt_sum = calcNeibourhood(box_sum, lt, rt, rb, lb);

double left_top = lt_sq_sum - lt_sum * lt_sum / N;
// ---------------------------------------------------
//右上
lt = Point(w, h - neibourhood);
lb = Point(w, h);
rt = Point(w + neibourhood, h - neibourhood);
rb = Point(w + neibourhood, h);

double rt_sq_sum = calcNeibourhood(box_sq_sum, lt, rt, rb, lb);
double rt_sum = calcNeibourhood(box_sum, lt, rt, rb, lb);

double right_top = rt_sq_sum - rt_sum * rt_sum / N;
// ---------------------------------------------------
//左下
lt = Point(w - neibourhood, h);
lb = Point(w - neibourhood, h + neibourhood);
rt = Point(w, h);
rb = Point(w, h + neibourhood);

double lb_sq_sum = calcNeibourhood(box_sq_sum, lt, rt, rb, lb);
double lb_sum = calcNeibourhood(box_sum, lt, rt, rb, lb);

double left_bottom = lb_sq_sum - lb_sum * lb_sum / N;
// ---------------------------------------------------
//右下
double rb_sq_sum = calcNeibourhood(box_sq_sum, lt, rt, rb, lb);
double rb_sum = calcNeibourhood(box_sum, lt, rt, rb, lb);

double right_bottom = rb_sq_sum - rb_sum * rb_sum / N;
// ---------------------------------------------------

// 比较以上四个区域的方差，选取最小方差区域
double min_var = min(min(min(left_top, right_top), left_bottom), right_bottom);

// ---------------------------------------------------
uchar kuwahara_val = 0;

// 计算矩形区域面积的均值
if (min_var == left_top)
kuwahara_val = floor(lt_sum / (neibourhood * neibourhood));
else if (min_var == right_top)
kuwahara_val = floor(rt_sum / (neibourhood * neibourhood));
else if (min_var == left_bottom)
kuwahara_val = floor(lb_sum / (neibourhood * neibourhood));
else if (min_var == right_bottom)
kuwahara_val = floor(rb_sum / (neibourhood * neibourhood));

// Clamp
if (kuwahara_val > 255)
kuwahara_val = 255;
else if (kuwahara_val < 0)
kuwahara_val = 0;

uchar* point = dst.ptr < uchar > (h);
point[w] = kuwahara_val;
}
}

// -------------------------
del2DDoublePointer(box_sum, size);
del2DDoublePointer(box_sq_sum, size);
}

效果图：



当天晚上实现的GLSL版（代码目前还比较粗糙，待优化一下再放上来）

以下是邻域大小不同的动态效果：