IIR递归高斯滤波总结
2016-04-22 15:41
309 查看
本文主要总结递归高斯算法的基本原理,分析,实现,指令优化实现等
原始高斯滤波算法计算量很大,网上资料很多,不赘述了,自行search
二维高斯滤波在图像处理中非常非常重要,最常见的是高斯模糊,做模糊滤镜,其本质是空域的低通滤波器,在本人另外的两篇博客同态滤波和MSRCR中有提及。对于同态滤波,需要对对数变换后的结果进行高通滤波,常规是变换到频域进行处理,也可以不进行频域变换,直接进行空域高斯滤波实现低通滤波,然后再用减法实现高通滤波。对于MSRCR算法,每个颜色通道需要进行三个强度的高斯滤波,复杂度可想而知。而且常规的高斯滤波需要考虑滤波窗大小,可以自己参考opencv里面的高斯滤波,可以自己搜一下复杂度分析,对于视频等实时处理,原始算法显然没什么实用性。
后面有人提出先进行水平一维滤波,然后再进行竖直方向一维滤波,对于一个点的滤波,计算数据量由NxN减少为2*N,但还是和窗口大小N有关
有看到一种递归减半滤波平面大小的递归实现,原理不太清楚,可参考如下资料,文中fastFilter函数的定义:
http://blog.csdn.net/smallstones/article/details/41787079
然后是本文主要介绍的IIR递归高斯滤波,这个滤波器能近似接近高斯滤波器,也是分成两次一维滤波,关键是不需要设定滤波窗口大小,复杂度跟窗口大小无关,对于一维滤波先进行一次正向滤波,然后进行一次逆向滤波,所以对于大小为NXN大小的二维矩阵,计算复杂度为O(NXN)
计算公式参考:
http://www.cnblogs.com/ImageVision/archive/2012/06/11/2545555.html
实现代码:
GIMP中有代码:http://gimp.sourcearchive.com/documentation/2.6.1/contrast-retinex_8c-source.html
[html] view
plain copy
typedef struct
{
gint scale;
gint nscales;
gint scales_mode;
gfloat cvar;
} RetinexParams;
/*
* Calculate the coefficients for the recursive filter algorithm
* Fast Computation of gaussian blurring.
*/
static void
compute_coefs3 (gauss3_coefs *c, gfloat sigma)
{
/*
* Papers: "Recursive Implementation of the gaussian filter.",
* Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995.
* formula: 11b computation of q
* 8c computation of b0..b1
* 10 alpha is normalization constant B
*/
gfloat q, q2, q3;
q = 0;
if (sigma >= 2.5)
{
q = 0.98711 * sigma - 0.96330;
}
else if ((sigma >= 0.5) && (sigma < 2.5))
{
q = 3.97156 - 4.14554 * (gfloat) sqrt ((double) 1 - 0.26891 * sigma);
}
else
{
q = 0.1147705018520355224609375;
}
q2 = q * q;
q3 = q * q2;
c->b[0] = (1.57825+(2.44413*q)+(1.4281 *q2)+(0.422205*q3));
c->b[1] = ( (2.44413*q)+(2.85619*q2)+(1.26661 *q3));
c->b[2] = ( -((1.4281*q2)+(1.26661 *q3)));
c->b[3] = ( (0.422205*q3));
c->B = 1.0-((c->b[1]+c->b[2]+c->b[3])/c->b[0]);
c->sigma = sigma;
c->N = 3;
/*
g_printerr ("q %f\n", q);
g_printerr ("q2 %f\n", q2);
g_printerr ("q3 %f\n", q3);
g_printerr ("c->b[0] %f\n", c->b[0]);
g_printerr ("c->b[1] %f\n", c->b[1]);
g_printerr ("c->b[2] %f\n", c->b[2]);
g_printerr ("c->b[3] %f\n", c->b[3]);
g_printerr ("c->B %f\n", c->B);
g_printerr ("c->sigma %f\n", c->sigma);
g_printerr ("c->N %d\n", c->N);
*/
}
static void
gausssmooth (gfloat *in, gfloat *out, gint size, gint rowstride, gauss3_coefs *c)
{
/*
* Papers: "Recursive Implementation of the gaussian filter.",
* Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995.
* formula: 9a forward filter
* 9b backward filter
* fig7 algorithm
*/
gint i,n, bufsize;
gfloat *w1,*w2;
/* forward pass */
bufsize = size+3;
size -= 1;
w1 = (gfloat *) g_try_malloc (bufsize * sizeof (gfloat));
w2 = (gfloat *) g_try_malloc (bufsize * sizeof (gfloat));
w1[0] = in[0];
w1[1] = in[0];
w1[2] = in[0];
for ( i = 0 , n=3; i <= size ; i++, n++)
{
w1
= (gfloat)(c->B*in[i*rowstride] +
((c->b[1]*w1[n-1] +
c->b[2]*w1[n-2] +
c->b[3]*w1[n-3] ) / c->b[0]));
}
/* backward pass */
w2[size+1]= w1[size+3];
w2[size+2]= w1[size+3];
w2[size+3]= w1[size+3];
for (i = size, n = i; i >= 0; i--, n--)
{
w2
= out[i * rowstride] = (gfloat)(c->B*w1
+
((c->b[1]*w2[n+1] +
c->b[2]*w2[n+2] +
c->b[3]*w2[n+3] ) / c->b[0]));
}
g_free (w1);
g_free (w2);
}
[html] view
plain copy
[cpp] view
plain copy
/*
@function: 1channel 2D gauss filter
*/
void gausssmooth1ch(float *src, float *dst, int width, int height, gauss3_coefs *c)
{
for (int y=0; y<height; y++) {
gausssmooth(src+y*width, dst+y*width, width, 1, c);
}
for (int x=0; x<width; x++) {
gausssmooth(dst+x, dst+x, height, width, c);
}
}
加速实现:
intel有详细的介绍和SSE实现:https://software.intel.com/en-us/articles/iir-gaussian-blur-filter-implementation-using-intel-advanced-vector-extensions
总而言之,这个算法非常非常实用
原始高斯滤波算法计算量很大,网上资料很多,不赘述了,自行search
二维高斯滤波在图像处理中非常非常重要,最常见的是高斯模糊,做模糊滤镜,其本质是空域的低通滤波器,在本人另外的两篇博客同态滤波和MSRCR中有提及。对于同态滤波,需要对对数变换后的结果进行高通滤波,常规是变换到频域进行处理,也可以不进行频域变换,直接进行空域高斯滤波实现低通滤波,然后再用减法实现高通滤波。对于MSRCR算法,每个颜色通道需要进行三个强度的高斯滤波,复杂度可想而知。而且常规的高斯滤波需要考虑滤波窗大小,可以自己参考opencv里面的高斯滤波,可以自己搜一下复杂度分析,对于视频等实时处理,原始算法显然没什么实用性。
后面有人提出先进行水平一维滤波,然后再进行竖直方向一维滤波,对于一个点的滤波,计算数据量由NxN减少为2*N,但还是和窗口大小N有关
有看到一种递归减半滤波平面大小的递归实现,原理不太清楚,可参考如下资料,文中fastFilter函数的定义:
http://blog.csdn.net/smallstones/article/details/41787079
然后是本文主要介绍的IIR递归高斯滤波,这个滤波器能近似接近高斯滤波器,也是分成两次一维滤波,关键是不需要设定滤波窗口大小,复杂度跟窗口大小无关,对于一维滤波先进行一次正向滤波,然后进行一次逆向滤波,所以对于大小为NXN大小的二维矩阵,计算复杂度为O(NXN)
计算公式参考:
http://www.cnblogs.com/ImageVision/archive/2012/06/11/2545555.html
实现代码:
GIMP中有代码:http://gimp.sourcearchive.com/documentation/2.6.1/contrast-retinex_8c-source.html
[html] view
plain copy
typedef struct
{
gint scale;
gint nscales;
gint scales_mode;
gfloat cvar;
} RetinexParams;
/*
* Calculate the coefficients for the recursive filter algorithm
* Fast Computation of gaussian blurring.
*/
static void
compute_coefs3 (gauss3_coefs *c, gfloat sigma)
{
/*
* Papers: "Recursive Implementation of the gaussian filter.",
* Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995.
* formula: 11b computation of q
* 8c computation of b0..b1
* 10 alpha is normalization constant B
*/
gfloat q, q2, q3;
q = 0;
if (sigma >= 2.5)
{
q = 0.98711 * sigma - 0.96330;
}
else if ((sigma >= 0.5) && (sigma < 2.5))
{
q = 3.97156 - 4.14554 * (gfloat) sqrt ((double) 1 - 0.26891 * sigma);
}
else
{
q = 0.1147705018520355224609375;
}
q2 = q * q;
q3 = q * q2;
c->b[0] = (1.57825+(2.44413*q)+(1.4281 *q2)+(0.422205*q3));
c->b[1] = ( (2.44413*q)+(2.85619*q2)+(1.26661 *q3));
c->b[2] = ( -((1.4281*q2)+(1.26661 *q3)));
c->b[3] = ( (0.422205*q3));
c->B = 1.0-((c->b[1]+c->b[2]+c->b[3])/c->b[0]);
c->sigma = sigma;
c->N = 3;
/*
g_printerr ("q %f\n", q);
g_printerr ("q2 %f\n", q2);
g_printerr ("q3 %f\n", q3);
g_printerr ("c->b[0] %f\n", c->b[0]);
g_printerr ("c->b[1] %f\n", c->b[1]);
g_printerr ("c->b[2] %f\n", c->b[2]);
g_printerr ("c->b[3] %f\n", c->b[3]);
g_printerr ("c->B %f\n", c->B);
g_printerr ("c->sigma %f\n", c->sigma);
g_printerr ("c->N %d\n", c->N);
*/
}
static void
gausssmooth (gfloat *in, gfloat *out, gint size, gint rowstride, gauss3_coefs *c)
{
/*
* Papers: "Recursive Implementation of the gaussian filter.",
* Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995.
* formula: 9a forward filter
* 9b backward filter
* fig7 algorithm
*/
gint i,n, bufsize;
gfloat *w1,*w2;
/* forward pass */
bufsize = size+3;
size -= 1;
w1 = (gfloat *) g_try_malloc (bufsize * sizeof (gfloat));
w2 = (gfloat *) g_try_malloc (bufsize * sizeof (gfloat));
w1[0] = in[0];
w1[1] = in[0];
w1[2] = in[0];
for ( i = 0 , n=3; i <= size ; i++, n++)
{
w1
= (gfloat)(c->B*in[i*rowstride] +
((c->b[1]*w1[n-1] +
c->b[2]*w1[n-2] +
c->b[3]*w1[n-3] ) / c->b[0]));
}
/* backward pass */
w2[size+1]= w1[size+3];
w2[size+2]= w1[size+3];
w2[size+3]= w1[size+3];
for (i = size, n = i; i >= 0; i--, n--)
{
w2
= out[i * rowstride] = (gfloat)(c->B*w1
+
((c->b[1]*w2[n+1] +
c->b[2]*w2[n+2] +
c->b[3]*w2[n+3] ) / c->b[0]));
}
g_free (w1);
g_free (w2);
}
[html] view
plain copy
[cpp] view
plain copy
/*
@function: 1channel 2D gauss filter
*/
void gausssmooth1ch(float *src, float *dst, int width, int height, gauss3_coefs *c)
{
for (int y=0; y<height; y++) {
gausssmooth(src+y*width, dst+y*width, width, 1, c);
}
for (int x=0; x<width; x++) {
gausssmooth(dst+x, dst+x, height, width, c);
}
}
加速实现:
intel有详细的介绍和SSE实现:https://software.intel.com/en-us/articles/iir-gaussian-blur-filter-implementation-using-intel-advanced-vector-extensions
总而言之,这个算法非常非常实用
相关文章推荐
- 使用递归高斯滤波器实现快速高斯模糊
- Recursive implementation of the Gaussian filter[翻译]
- github(其他类似github)下载到本地配置
- Android中的LayoutInflater
- 删除字符串中的空格
- hbase的thrift接口
- OS-->Git操作演练(项目实用命令)
- Xcode插件指南,让你的开发更便捷
- <<程序员面试宝典>>读书笔记 3
- Bootstrap 3学习笔记 -栅格
- 实验2 作业调度
- hibernate--多对多单向关联 (重点!!!)
- 一张图理解订阅者模式
- 用calibre脚本抓取ipdf的epub 3.0.1规范
- 进程的状态与转换
- Java的热部署(后期完善)
- smtp协议的基本命令
- 什么是盐加密 为什么使用盐加密密码
- servlet详解(第一篇)
- webix前端界面框架