图像滤波常用算法实现及原理解析
导读
图像滤波是一种非常重要的图像处理技术，本文详细介绍了四种常见的图像滤波算法，并附上源码，包括自适应中值滤波、高斯滤波、双边滤波和导向滤波。
前言
本文介绍四种常见的图像滤波算法，并附上源码。图像滤波是一种非常重要的图像处理技术，现在大火的卷积神经网络其实也是滤波的一种，都是用卷积核去提取图像的特征模式。不过，传统的滤波，使用的卷积核是固定的参数，是由经验非常丰富的人去手动设计的，也称为手工特征。而卷积神经网络的卷积核参数初始时未知的，根据不同的任务由数据和神经网络反向传播算法去学习得到的参数，更能适应于不同的任务。
目录
- 自适应中值滤波
- 高斯滤波
- 双边滤波
- 导向滤波
自适应中值滤波
中值滤波器是一种常用的非线性滤波器，其基本原理是：选择待处理像素的一个邻域中各像素值的中值来代替待处理的像素。主要功能使某像素的灰度值与周围领域内的像素比较接近，从而消除一些孤立的噪声点，所以中值滤波器能够很好的消除椒盐噪声。不仅如此，中值滤波器在消除噪声的同时，还能有效的保护图像的边界信息，不会对图像造成很大的模糊（相比于均值滤波器）。
中值滤波器的效果受滤波窗口尺寸的影响较大，在消除噪声和保护图像的细节存在着矛盾：滤波窗口较小，则能很好的保护图像中的某些细节，但对噪声的过滤效果就不是很好，因为实际中的噪声不可能只占一个像素位置；反之，窗口尺寸较大有较好的噪声过滤效果，但是会对图像造成一定的模糊。另外，根据中值滤波器原理，如果在滤波窗口内的噪声点的个数大于整个窗口内非噪声像素的个数，则中值滤波就不能很好的过滤掉噪声。
自适应中值滤波器
常规的中值滤波器，在噪声的密度不是很大的情况下，效果不错。但是当噪声出现的概率较高时，常规的中值滤波的效果就不是很好了。有一个选择就是增大滤波器的窗口大小，这虽然在一定程度上能解决上述的问题，但是会给图像造成较大的模糊。
常规的中值滤波器的窗口尺寸是固定大小不变的，就不能同时兼顾去噪和保护图像的细节。这时就要寻求一种改变，根据预先设定好的条件，在滤波的过程中，动态的改变滤波器的窗口尺寸大小，这就是自适应中值滤波器 Adaptive Median Filter。在滤波的过程中，自适应中值滤波器会根据预先设定好的条件，改变滤波窗口的尺寸大小，同时还会根据一定的条件判断当前像素是不是噪声，如果是则用邻域中值替换掉当前像素；不是，则不作改变。
自适应中值滤波器有三个目的：
- 滤除椒盐噪声
- 平滑其他非脉冲噪声
- 尽可能的保护图像中细节信息，避免图像边缘的细化或者粗化。
自适应中值滤波算法描述
自适应滤波器不但能够滤除概率较大的椒盐噪声，而且能够更好的保护图像的细节，这是常规的中值滤波器做不到的。自适应的中值滤波器也需要一个矩形的窗口 ，和常规中值滤波器不同的是这个窗口的大小会在滤波处理的过程中进行改变（增大）。需要注意的是，滤波器的输出是一个像素值，该值用来替换点 处的像素值，点 是滤波窗口的中心位置。
在描述自适应中值滤波器时需要用到如下的符号：
- 窗口中的最小灰度值
- 窗口中的最大灰度值
- 窗口中的灰度值的中值
- 表示坐标 处的灰度值
- 允许的最大窗口尺寸
自适应中值滤波器有两个处理过程，分别记为：和。
A :
如果A1 > 0 且 A2 < 0，跳转到 B；
否则，增大窗口的尺寸 如果增大后窗口的尺寸 ，则重复A过程。否则，输出 𝑍𝑚𝑒𝑑
B:
如果B1 > 0 且 B2 <0则输出 ,否则输出
自适应中值滤波原理说明
过程A的目的是确定当前窗口内得到中值 𝑍𝑚𝑒𝑑 是否是噪声。如果 𝑍𝑚𝑖𝑛<𝑍𝑚𝑒𝑑<𝑍𝑚𝑎𝑥，则𝑍𝑚𝑒𝑑不是噪声，这时转到过程B测试当前窗口的中心位置的像素 𝑍𝑥𝑦 是否是一个噪声点。如果 𝑍𝑚𝑖𝑛<𝑍𝑥𝑦<𝑍𝑚𝑎𝑥，则 𝑍𝑥𝑦 不是一个噪声，此时输出𝑍𝑥𝑦 ;如果不满足上述条件，则可判定 𝑍𝑥𝑦 是噪声，这是输出中值 𝑍𝑚𝑒𝑑 （在A中已经判断出 𝑍𝑚𝑒𝑑 不是噪声）。
如果在过程A中，得到则 𝑍𝑚𝑒𝑑 不符合条件 𝑍𝑚𝑖𝑛<𝑍𝑚𝑒𝑑<𝑍𝑚𝑎𝑥 ，则可判断得到的中值 𝑍𝑚𝑒𝑑 是一个噪声。在这种情况下，需要增大滤波器的窗口尺寸，在一个更大的范围内寻找一个非噪声点的中值，直到找到一个非噪声的中值，跳转到B；或者，窗口的尺寸达到了最大值，这时返回找到的中值，退出。
从上面分析可知，噪声出现的概率较低，自适应中值滤波器可以较快的得出结果，不需要去增加窗口的尺寸；反之，噪声的出现的概率较高，则需要增大滤波器的窗口尺寸，这也符合种中值滤波器的特点：噪声点比较多时，需要更大的滤波器窗口尺寸。
算法实现
有了算法的详细描述，借助于OpenCV对图像的读写，自适应中值滤波器实现起来也不是很困难。首先定义滤波器最小的窗口尺寸以及最大的窗口尺寸。要进行滤波处理，首先要扩展图像的边界，以便对图像的边界像素进行处理。copyMakeBorder根据选择的BorderTypes使用不同的值扩充图像的边界像素，具体可参考OpenCV的文档信息。下面就是遍历图像的像素，对每个像素进行滤波处理。需要注意一点，不论滤波器多么的复杂，其每次的滤波过程，都是值返回一个值，来替换掉当前窗口的中心的像素值。函数adpativeProcess就是对当前像素的滤波过程，其代码如下：
uchar adaptiveProcess(const Mat &im, int row,int col,int kernelSize,int maxSize)
{
vector<uchar> pixels;
for (int a = -kernelSize / 2; a <= kernelSize / 2; a++)
for (int b = -kernelSize / 2; b <= kernelSize / 2; b++)
{
pixels.push_back(im.at<uchar>(row + a, col + b));
}
sort(pixels.begin(), pixels.end());
auto min = pixels[0];
auto max = pixels[kernelSize * kernelSize - 1];
auto med = pixels[kernelSize * kernelSize / 2];
auto zxy = im.at<uchar>(row, col);
if (med > min && med < max)
{
// to B
if (zxy > min && zxy < max)
return zxy;
else
return med;
}
else
{
kernelSize += 2;
if (kernelSize <= maxSize)
return adpativeProcess(im, row, col, kernelSize, maxSize); // 增大窗口尺寸，继续A过程。
else
return med;
}
}有了上面这个函数，剩下的只需要对全部像素做一个遍历即可，更为完整的代码，请见我的Github地址：
https://github.com/zhangqizky/common-image-filteringgithub.com
高斯滤波
高斯滤波也是一种非常常见的滤波方法，其核的形式为：
其中是图像中的点的坐标, 为标准差，高斯模板就是利用这个函数来计算的，x和y都是代表，以核中心点为坐标原点的坐标值。这里想说一下 的作用，当 比较小的时候，生成的高斯模板中心的系数比较大，而周围的系数比较小，这样对图像的平滑效果不明显。而当 比较大时，生成的模板的各个系数相差就不是很大，比较类似于均值模板，对图像的平滑效果比较明显。
高斯滤波没有特别多可说的，最主要的作用是滤除高斯噪声，即符合正态分布的噪声。
实现的方式有两种，第一种是按照公式暴力实现，代码如下:
//O(m * n * ksize^2)
void GaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma)
{
CV_Assert(src.channels() || src.channels() == 3); //只处理3通道或单通道的图片
double **GaussianTemplate = new double *[ksize];
for(int i = 0; i < ksize; i++){
GaussianTemplate[i] = new double [ksize];
}
generateGaussianTemplate(GaussianTemplate, ksize, sigma);
//padding
int border = ksize / 2;
copyMakeBorder(src, dst, border, border, border, border, BORDER_CONSTANT);
int channels = dst.channels();
int rows = dst.rows - border;
int cols = dst.cols - border;
for(int i = border; i < rows; i++){
for(int j = border; j< cols; j++){
double sum[3] = {0};
for(int a = -border; a <= border; a++){
for(int b = -border; b <= border; b++){
if(channels == 1){
sum[0] += GaussianTemplate[border+a][border+b] * dst.at<uchar>(i+a, j+b);
}else if(channels == 3){
Vec3b rgb = dst.at<Vec3b>(i+a, j+b);
auto k = GaussianTemplate[border+a][border+b];
sum[0] += k * rgb[0];
sum[1] += k * rgb[1];
sum[2] += k * rgb[2];
}
}
}
for(int k = 0; k < channels; k++){
if(sum[k] < 0) sum[k] = 0;
else if(sum[k] > 255) sum[k] = 255;
}
if(channels == 1){
dst.at<uchar >(i, j) = static_cast<uchar >(sum[0]);
}else if(channels == 3){
Vec3b rgb = {static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2])};
dst.at<Vec3b>(i, j) = rgb;
}
}
}
for(int i = 0; i < ksize; i++)
delete[] GaussianTemplate[i];
delete[] GaussianTemplate;
}当核比较大时，高斯滤波会比较费时，此时可以使用分离X和Y通道的形式来实现可分离的高斯滤波。为什么可以这么做？因为高斯函数中没有这样的耦合项，即x和y是相对独立的，此时就可以将两个维度分离开来。
//分离实现高斯滤波
//O(m*n*k)
void separateGaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma){
assert(src.channels()==1 || src.channels() == 3); //只处理单通道或者三通道图像
//生成一维的
double *matrix = new double[ksize];
double sum = 0;
int origin = ksize / 2;
for(int i = 0; i < ksize; i++){
double g = exp(-(i-origin) * (i-origin) / (2 * sigma * sigma));
sum += g;
matrix[i] = g;
}
for(int i = 0; i < ksize; i++) matrix[i] /= sum;
int border = ksize / 2;
copyMakeBorder(src, dst, border, border, border, border, BORDER_CONSTANT);
int channels = dst.channels();
int rows = dst.rows - border;
int cols = dst.cols - border;
//水平方向
for(int i = border; i < rows; i++){
for(int j = border; j < cols; j++){
double sum[3] = {0};
for(int k = -border; k<=border; k++){
if(channels == 1){
sum[0] += matrix[border + k] * dst.at<uchar>(i, j+k);
}else if(channels == 3){
Vec3b rgb = dst.at<Vec3b>(i, j+k);
sum[0] += matrix[border+k] * rgb[0];
sum[1] += matrix[border+k] * rgb[1];
sum[2] += matrix[border+k] * rgb[2];
}
}
for(int k = 0; k < channels; k++){
if(sum[k] < 0) sum[k] = 0;
else if(sum[k] > 255) sum[k] = 255;
}
if(channels == 1)
dst.at<Vec3b>(i, j) = static_cast<uchar>(sum[0]);
else if(channels == 3){
Vec3b rgb = {static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2])};
dst.at<Vec3b>(i, j) = rgb;
}
}
}
//竖直方向
for(int i = border; i < rows; i++){
for(int j = border; j < cols; j++){
double sum[3] = {0};
for(int k = -border; k<=border; k++){
if(channels == 1){
sum[0] += matrix[border + k] * dst.at<uchar>(i+k, j);
}else if(channels == 3){
Vec3b rgb = dst.at<Vec3b>(i+k, j);
sum[0] += matrix[border+k] * rgb[0];
sum[1] += matrix[border+k] * rgb[1];
sum[2] += matrix[border+k] * rgb[2];
}
}
for(int k = 0; k < channels; k++){
if(sum[k] < 0) sum[k] = 0;
else if(sum[k] > 255) sum[k] = 255;
}
if(channels == 1)
dst.at<Vec3b>(i, j) = static_cast<uchar>(sum[0]);
else if(channels == 3){
Vec3b rgb = {static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2])};
dst.at<Vec3b>(i, j) = rgb;
}
}
}
delete [] matrix;
}同样，完整的代码请见：
https://github.com/zhangqizky/common-image-filtering/tree/maingithub.com
双边滤波
双边滤波是一种非线性滤波方法，是结合了图像的邻近度和像素值相似度的一种折中，在滤除噪声的同时可以保留原图的边缘信息。整个双边滤波是由两个函数构成：一个函数是由空间距离决定的滤波器系数，另外一个诗由像素差值决定的滤波器系数。整个双边滤波的公式如下：
其中权重系数 取决于定义域核：
和值域核
的乘积。其中定义域核影响的是空间位置，如果把图像看成一个二维函数，那么定义域就是图像的坐标，值域就是该坐标处对应的像素值。定义域核就是普通的高斯核，全局使用一个就可以。但值域核是需要对每个像素点滑动进行计算的。
那么如何理解双边滤波呢
高斯滤波的滤波核的意义是，滤波后的像素值等于窗口内的像素值的加权平均值，权值系数是符合高斯分布，距离该点越近，权值越大。但是没有考虑像素值与当前点的差距。现在加上值域核，意义就在，滤波后当前点的像素值还会受到领域内像素值与自身的像素值差异的影响，不仅仅是距离来决定。这样，在平缓的区域里，由于像素值差异非常小，则值域的权重趋向于1，所以双边滤波就近似为高斯滤波。而在边缘区域中，由于像素值的差异比较大，则值域核趋向于0，权重下降，即当前像素受到领域内像素影响比较小，从而保留了边缘信息。
双边滤波的代码
opencv中提供了bilateralFilter()函数来实现双边滤波操作，其原型如下：
void cv::bilateralFilter(InputArray src,
OutputArray dst,
int d,
double sigmaColor,
double sigmaSpace,
int borderType = BORDER_DEFAULT
)
- InputArray src: 输入图像，可以是Mat类型，图像必须是8位整型或浮点型单通道、三通道的图像。
- OutputArray dst: 输出图像，和原图像有相同的尺寸和类型。
- int d: 表示在过滤过程中每个像素邻域的直径范围。如果这个值是非正数，则函数会从第五个参数sigmaSpace计算该值。
- double sigmaColor: 颜色空间过滤器的值，这个参数的值月大，表明该像素邻域内有越宽广的颜色会被混合到一起，产生较大的半相等颜色区域。（这个参数可以理解为值域核的 和 )
- double sigmaSpace: 坐标空间中滤波器的sigma值，如果该值较大，则意味着越远的像素将相互影响，从而使更大的区域中足够相似的颜色获取相同的颜色。当d>0时，d指定了邻域大小且与sigmaSpace无关，否则d正比于sigmaSpace. （这个参数可以理解为空间域核的 和 ）
- int borderType=BORDER_DEFAULT: 用于推断图像外部像素的某种边界模式，有默认值BORDER_DEFAULT.
具体代码如下：
#include <iostream>
#include <opencv2/opencv.hpp>
using namespace std;
using namespace cv;
//定义全局变量
const int g_ndMaxValue = 100;
const int g_nsigmaColorMaxValue = 200;
const int g_nsigmaSpaceMaxValue = 200;
int g_ndValue;
int g_nsigmaColorValue;
int g_nsigmaSpaceValue;
Mat g_srcImage;
Mat g_dstImage;
//定义回调函数
void on_bilateralFilterTrackbar(int, void*);
int main()
{
g_srcImage = imread("lena.jpg");
//判断图像是否加载成功
if(g_srcImage.empty())
{
cout << "图像加载失败!" << endl;
return -1;
}
else
cout << "图像加载成功!" << endl << endl;
namedWindow("原图像", WINDOW_AUTOSIZE);
imshow("原图像", g_srcImage);
//定义输出图像窗口属性和轨迹条属性
namedWindow("双边滤波图像", WINDOW_AUTOSIZE);
g_ndValue = 10;
g_nsigmaColorValue = 10;
g_nsigmaSpaceValue = 10;
char dName[20];
sprintf(dName, "邻域直径 %d", g_ndMaxValue);
char sigmaColorName[20];
sprintf(sigmaColorName, "sigmaColor %d", g_nsigmaColorMaxValue);
char sigmaSpaceName[20];
sprintf(sigmaSpaceName, "sigmaSpace %d", g_nsigmaSpaceMaxValue);
//创建轨迹条
createTrackbar(dName, "双边滤波图像", &g_ndValue, g_ndMaxValue, on_bilateralFilterTrackbar);
on_bilateralFilterTrackbar(g_ndValue, 0);
createTrackbar(sigmaColorName, "双边滤波图像", &g_nsigmaColorValue,
g_nsigmaColorMaxValue, on_bilateralFilterTrackbar);
on_bilateralFilterTrackbar(g_nsigmaColorValue, 0);
createTrackbar(sigmaSpaceName, "双边滤波图像", &g_nsigmaSpaceValue,
g_nsigmaSpaceMaxValue, on_bilateralFilterTrackbar);
on_bilateralFilterTrackbar(g_nsigmaSpaceValue, 0);
waitKey(0);
return 0;
}
void on_bilateralFilterTrackbar(int, void*)
{
bilateralFilter(g_srcImage, g_dstImage, g_ndValue, g_nsigmaColorValue, g_nsigmaSpaceValue);
imshow("双边滤波图像", g_dstImage);
}导向滤波
需要有高斯滤波和双边滤波的相关知识背景才能更好的理解导向滤波。在导向滤波中，首先利用了局部线性模型。这个模型认为某函数上一点与其近邻部分的点成线性关系，一个复杂的函数就可以用很多局部的线性函数来表示，当需要求该函数上某一点的值时，只需要计算所有包含该点的线性函数的值并取平均值即可。这种模型，在表示非解析函数上，非常有用。
同理，我们可以认为图像是一个二维函数，并且假设该函数的输出与输入在一个二维窗口内满足线性关系，如下：
其中，是输出像素的值，是输入图像的值，和是像素索引，和是当窗口中心位于k时该线性函数的系数。其实，输入图像不一定是待滤波的图像本身，也可以是其他图像即引导图像，这也是为何称为引导滤波的原因。对上式两边取梯度，可以得到:
即当输入图像有梯度时，输出也有类似的梯度，现在可以解释为什么引导滤波有边缘保持特性了。下一步是求出线性函数的系数，也就是线性回归，即希望拟合函数的输出值与真实值之间的差距最小，也就是让下式最小:
这里只能是待滤波图像，并不像那样可以是其他图像。同时，之前的系数用于防止求得的过大，也是调节滤波器滤波效果的重要参数(相当于L2正则化的权重惩罚)。接下来利用最小二乘法的原理令 和 得到2个二元一次方程，求解得到：
其中 是 在窗口的平均值， 是 在窗口 的方差， 是窗口 中的像素个数， 是待滤波图像在窗口 中的均值。在计算每个窗口的线性系数时，我们可以发现一个像素会被多个窗口包含，也就是说，每个像素都由多个线性函数所描述。因此，如之前所说，要具体求某一点的输出值时，只需将所有包含该点的线性函数值平均即可，如下:
这里， 是所有包含像素 的窗口， 是其中心位置。
当把引导滤波用作边缘保持滤波器时，往往有 ，如果 ，显然是为最小值的解，从上式可以看出，这时的滤波器没有任何作用，将输入原封不动的输出。如果 ，在像素强度变化小的区域（或单色区域），有近似于（或等于0，而近似于（或等于） ，即做了一个加权均值滤波；而在变化大的区域，近似于1，近似于0，对图像的滤波效果很弱，有助于保持边缘。而 的作用就是界定什么是变化大，什么是变化小。在窗口大小不变的情况下，随着的增大，滤波效果越明显。
在滤波效果上，引导滤波和双边滤波差不多，然后在一些细节上，引导滤波较好(在PS的磨皮美白中，经过实践，效果更好)。引导滤波最大的优势在于，可以写出时间复杂度与窗口大小无关的算法，因此在使用大窗口处理图片时，其效率更高。
同样，OpenCV中也有导向滤波的接口。具体代码如下：
void cv::ximgproc::guidedFilter ( InputArray guide,
InputArray src,
OutputArray dst,
int radius,
double eps,
int dDepth = -1
)
guide | 引导图像，3通道，如果大于3通道则只有前三个会被用到 |
|---|---|
src | 待滤波图像 |
dst | 输出图像 |
radius | 导向滤波的窗口 |
eps | 正则化参数 |
dDepth | 可选，图像的深度参数 |
这边有个基于scipy实现的python代码，可以参考一下：
import numpy as np
import scipy as sp
import scipy.ndimage
def box(img, r):
""" O(1) box filter
img - >= 2d image
r - radius of box filter
"""
(rows, cols) = img.shape[:2]
imDst = np.zeros_like(img)
tile = [1] * img.ndim
tile[0] = r
imCum = np.cumsum(img, 0)
imDst[0:r+1, :, ...] = imCum[r:2*r+1, :, ...]
imDst[r+1:rows-r, :, ...] = imCum[2*r+1:rows, :, ...] - imCum[0:rows-2*r-1, :, ...]
imDst[rows-r:rows, :, ...] = np.tile(imCum[rows-1:rows, :, ...], tile) - imCum[rows-2*r-1:rows-r-1, :, ...]
tile = [1] * img.ndim
tile[1] = r
imCum = np.cumsum(imDst, 1)
imDst[:, 0:r+1, ...] = imCum[:, r:2*r+1, ...]
imDst[:, r+1:cols-r, ...] = imCum[:, 2*r+1 : cols, ...] - imCum[:, 0 : cols-2*r-1, ...]
imDst[:, cols-r: cols, ...] = np.tile(imCum[:, cols-1:cols, ...], tile) - imCum[:, cols-2*r-1 : cols-r-1, ...]
return imDst
def _gf_color(I, p, r, eps, s=None):
""" Color guided filter
I - guide image (rgb)
p - filtering input (single channel)
r - window radius
eps - regularization (roughly, variance of non-edge noise)
s - subsampling factor for fast guided filter
"""
fullI = I
fullP = p
if s is not None:
I = sp.ndimage.zoom(fullI, [1/s, 1/s, 1], order=1)
p = sp.ndimage.zoom(fullP, [1/s, 1/s], order=1)
r = round(r / s)
h, w = p.shape[:2]
N = box(np.ones((h, w)), r)
mI_r = box(I[:,:,0], r) / N
mI_g = box(I[:,:,1], r) / N
mI_b = box(I[:,:,2], r) / N
mP = box(p, r) / N
# mean of I * p
mIp_r = box(I[:,:,0]*p, r) / N
mIp_g = box(I[:,:,1]*p, r) / N
mIp_b = box(I[:,:,2]*p, r) / N
# per-patch covariance of (I, p)
covIp_r = mIp_r - mI_r * mP
covIp_g = mIp_g - mI_g * mP
covIp_b = mIp_b - mI_b * mP
# symmetric covariance matrix of I in each patch:
# rr rg rb
# rg gg gb
# rb gb bb
var_I_rr = box(I[:,:,0] * I[:,:,0], r) / N - mI_r * mI_r;
var_I_rg = box(I[:,:,0] * I[:,:,1], r) / N - mI_r * mI_g;
var_I_rb = box(I[:,:,0] * I[:,:,2], r) / N - mI_r * mI_b;
var_I_gg = box(I[:,:,1] * I[:,:,1], r) / N - mI_g * mI_g;
var_I_gb = box(I[:,:,1] * I[:,:,2], r) / N - mI_g * mI_b;
var_I_bb = box(I[:,:,2] * I[:,:,2], r) / N - mI_b * mI_b;
a = np.zeros((h, w, 3))
for i in range(h):
for j in range(w):
sig = np.array([
[var_I_rr[i,j], var_I_rg[i,j], var_I_rb[i,j]],
[var_I_rg[i,j], var_I_gg[i,j], var_I_gb[i,j]],
[var_I_rb[i,j], var_I_gb[i,j], var_I_bb[i,j]]
])
covIp = np.array([covIp_r[i,j], covIp_g[i,j], covIp_b[i,j]])
a[i,j,:] = np.linalg.solve(sig + eps * np.eye(3), covIp)
b = mP - a[:,:,0] * mI_r - a[:,:,1] * mI_g - a[:,:,2] * mI_b
meanA = box(a, r) / N[...,np.newaxis]
meanB = box(b, r) / N
if s is not None:
meanA = sp.ndimage.zoom(meanA, [s, s, 1], order=1)
meanB = sp.ndimage.zoom(meanB, [s, s], order=1)
q = np.sum(meanA * fullI, axis=2) + meanB
return q
def _gf_gray(I, p, r, eps, s=None):
""" grayscale (fast) guided filter
I - guide image (1 channel)
p - filter input (1 channel)
r - window raidus
eps - regularization (roughly, allowable variance of non-edge noise)
s - subsampling factor for fast guided filter
"""
if s is not None:
Isub = sp.ndimage.zoom(I, 1/s, order=1)
Psub = sp.ndimage.zoom(p, 1/s, order=1)
r = round(r / s)
else:
Isub = I
Psub = p
(rows, cols) = Isub.shape
N = box(np.ones([rows, cols]), r)
meanI = box(Isub, r) / N
meanP = box(Psub, r) / N
corrI = box(Isub * Isub, r) / N
corrIp = box(Isub * Psub, r) / N
varI = corrI - meanI * meanI
covIp = corrIp - meanI * meanP
a = covIp / (varI + eps)
b = meanP - a * meanI
meanA = box(a, r) / N
meanB = box(b, r) / N
if s is not None:
meanA = sp.ndimage.zoom(meanA, s, order=1)
meanB = sp.ndimage.zoom(meanB, s, order=1)
q = meanA * I + meanB
return q
def _gf_colorgray(I, p, r, eps, s=None):
""" automatically choose color or gray guided filter based on I's shape """
if I.ndim == 2 or I.shape[2] == 1:
return _gf_gray(I, p, r, eps, s)
elif I.ndim == 3 and I.shape[2] == 3:
return _gf_color(I, p, r, eps, s)
else:
print("Invalid guide dimensions:", I.shape)
def guided_filter(I, p, r, eps, s=None):
""" run a guided filter per-channel on filtering input p
I - guide image (1 or 3 channel)
p - filter input (n channel)
r - window raidus
eps - regularization (roughly, allowable variance of non-edge noise)
s - subsampling factor for fast guided filter
"""
if p.ndim == 2:
p3 = p[:,:,np.newaxis]
out = np.zeros_like(p3)
for ch in range(p3.shape[2]):
out[:,:,ch] = _gf_colorgray(I, p3[:,:,ch], r, eps, s)
return np.squeeze(out) if p.ndim == 2 else out
def test_gf():
import imageio
cat = imageio.imread('cat.bmp').astype(np.float32) / 255
tulips = imageio.imread('tulips.bmp').astype(np.float32) / 255
r = 8
eps = 0.05
cat_smoothed = guided_filter(cat, cat, r, eps)
cat_detail = cat / cat_smoothed
print(cat_detail.shape)
cat_smoothed_s4 = guided_filter(cat, cat, r, eps, s=4)
imageio.imwrite('cat_smoothed.png', cat_smoothed)
imageio.imwrite('cat_smoothed_s4.png', cat_smoothed_s4)
imageio.imwrite('cat_smoothed_detailed.png',cat_detail)
tulips_smoothed4s = np.zeros_like(tulips)
tulips_detailed = np.zeros_like(tulips)
for i in range(3):
tulips_smoothed4s[:,:,i] = guided_filter(tulips, tulips[:,:,i], r, eps, s=4)
tulips_detailed = tulips / tulips_smoothed4s
imageio.imwrite('tulips_detailed.png',tulips_detailed)
imageio.imwrite('tulips_smoothed4s.png', tulips_smoothed4s)
tulips_smoothed = np.zeros_like(tulips)
for i in range(3):
tulips_smoothed[:,:,i] = guided_filter(tulips, tulips[:,:,i], r, eps)
imageio.imwrite('tulips_smoothed.png', tulips_smoothed)
if __name__ == '__main__':
test_gf()一副图像，经过mask是图像本身的导向滤波之后，得到一张细节图和一张滤波图。下面从左到右分别是原图，细节图和滤波图。其实这是现在很多low-level领域的预处理步骤。拿到细节图之后可以用卷积神经网络做下面的处理。
