数字图像处理图像均值滤波的MATLAB程序

1. 滤波器介绍

滤波器作为图像处理课程的重要内容，大致可分为两类，空域滤波器和频率域滤波器。本文主要介绍常用的四种滤波器：中值滤波器、均值滤波器、高斯滤波器、双边滤波器，并基于opencv做出实现。空域的滤波器一般可以通过模板对原图像进行卷积。

注意：空域滤波器和频率域滤波器对比

1）空间域指图像本身，空域变换直接对图像中的像素进行操作。

2）图像变换是将图像从空间域变换到某变换域（如 傅立叶变换中的频率域）的数学变换，在变换域 中进行处理，然后通过反变换把处理结果返回到空间域。

3）图像在空域上具有很强的相关性，借助于正交变 换可使在空域的复杂计算转换到频域后得到简化

4）借助于频域特性的分析，将更有利于获得图像的 各种特性和进行特殊处理

 

2、理论知识：

图像的空域滤波无非两种情况，线性滤波和非线性滤波。

      滤波的意思就是对原图像的每个像素周围一定范围内的像素进行运算，运算的范围就称为掩膜。而运算就分两种了，如果运算只是对各像素灰度值进行简单处理（如乘一个权值）最后求和，就称为线性滤波；而如果对像素灰度值的运算比较复杂，而不是最后求和的简单运算，则是非线性滤波；如求一个像素周围3x3范围内最大值、最小值、中值、均值等操作都不是简单的加权，都属于非线性滤波。

常见的线性滤波有：均值滤波、高斯滤波、盒子滤波、拉普拉斯滤波等等，通常线性滤波器之间只是模版系数不同。

非线性滤波利用原始图像跟模版之间的一种逻辑关系得到结果，如最值滤波器，中值滤波器和双边滤波器等。

1、线性滤波

线性滤波器表达公式：

，其中均值滤波器和高斯滤波器属于线性滤波器，首先看这两种滤波器

均值滤波器：

模板：

 

从待处理图像首元素开始用模板对原始图像进行卷积，均值滤波直观地理解就是用相邻元素灰度值的平均值代替该元素的灰度值。

高斯滤波器：

高斯滤波一般针对的是高斯噪声，能够很好的抑制图像输入时随机引入的噪声，将像素点跟邻域像素看作是一种高斯分布的关系，它的操作是将图像和一个高斯核进行卷积操作： 

 

模板：通过高斯内核函数产生的

高斯内核函数：

例如3*3的高斯内核模板：

中值滤波：同样是空间域的滤波，主题思想是取相邻像素的点，然后对相邻像素的点进行排序，取中点的灰度值作为该像素点的灰度值。

统计排序滤波器，对椒盐噪声有很好的抑制

详细请参考：数字图像处理之椒盐噪声和中值滤波

中值滤波将窗口函数里面的所有像素进行排序取得中位数来代表该窗口中心的像素值，对椒盐噪声和脉冲噪声的抑制效果特别好，同时又能保留边缘细节，用公式表示是： 

双边滤波(Bilateral filter)也是一种非线性的滤波方法，是结合图像的空间邻近度和像素值相似度的一种折衷处理，同时考虑空域信息和灰度相似性，达到保边去噪的目的，具有简单，非迭代、局部的特点，它比高斯滤波多了一个高斯方差σd，用公式表示就是： 

w(x,y)为加权系数，取决于定义域核和值域核的乘积。

注意：

1）均值模糊无法克服边缘像素信息丢失的缺陷，原因是均值滤波是基于平均权重的。

2）高斯模糊部分克服了该缺陷，但无法完全避免，因为没有考虑像素值的不同。

3）高斯双边模糊-是边缘保留额滤波方法，避免了边缘信息丢失，保留了图像轮廓不变。

3. 实验

 

结论：从滤波的结果可以看出各种滤波算法对图像的作用非常不同，有些变化非常大，有些甚至跟原图一样。在实际应用时，应根据噪声的特点、期望的图像和边缘特征等来选择合适的滤波器，这样才能发挥图像滤波的最大优点。

 

 

4. C++实现

4.1均值滤波

static void exchange(int& a, int& b)
 {    
     int t = 0;
     t = a;
     a = b;
     b = t;
 }
  
 static void bubble_sort(int* K, int lenth)
 {
     for (int i = 0; i < lenth; i++)
         for (int j = i + 1; j < lenth; j++)
         {
             if (K[i]>K[j])
                 exchange(K[i], K[j]);
         }
 }
 ///产生二维的高斯内核
 static cv::Mat generate_gassian_kernel(double u, double sigma, cv::Size size)
 {
     int width = size.width;
     int height = size.height;
     cv::Mat gassian_kernel(cv::Size(width, height), CV_64FC1);
     double sum = 0;
     double sum_sum = 0;
     for (int i = 0; i < width; i++)
         for (int j = 0; j < height; j++)
         {
             sum = 1.0 / 2.0 / CV_PI / sigma / sigma * exp(-1.0 * ((i - width / 2)*(i - width / 2) + (j - width / 2)*(j - width / 2)) / 2.0 / sigma / sigma);
             sum_sum += sum;
             gassian_kernel.ptr<double>(i)[j] = sum;
         }
     for (int i = 0; i < width; i++)
         for (int j = 0; j < height; j++)
         {
             gassian_kernel.ptr<double>(i)[j] /= sum_sum;
         }
     return gassian_kernel;
 }
 ///均值滤波
 void lmt_main_blur(cv::Mat& img_in, cv::Mat& img_out, int kernel_size)
 {
     img_out = img_in.clone();
     cv::Mat mat1;
     cv::copyMakeBorder(img_in, mat1, kernel_size, kernel_size, kernel_size, kernel_size, cv::BORDER_REPLICATE);
  
     int cols = mat1.cols;
     int rows = mat1.rows;
     int channels = img_out.channels();
     const uchar* const pt = mat1.ptr<uchar>(0);
     uchar* pt_out = img_out.ptr<uchar>(0);
  
     for (int i = kernel_size; i < rows - kernel_size; i++)
     {
         for (int j = kernel_size; j < cols - kernel_size; j++)
         {
             if (channels == 1)
             {
                 long long int sum_pixel = 0;
                 for (int m = -1 * kernel_size; m < kernel_size; m++)
                     for (int n = -1 * kernel_size; n < kernel_size; n++)
                     {
                         sum_pixel += pt[(i + m)*cols + (j + n)];
                     }
                 img_out.ptr<uchar>(i - kernel_size)[j - kernel_size] = (double)sum_pixel / (kernel_size*kernel_size * 4);
             }
             else if (channels == 3)
             {
                 long long int sum_pixel = 0;
                 long long int sum_pixel1 = 0;
                 long long int sum_pixel2 = 0;
                 for (int m = -1 * kernel_size; m < kernel_size; m++)
                     for (int n = -1 * kernel_size; n < kernel_size; n++)
                     {
                         sum_pixel += pt[((i + m)*cols + (j + n))*channels + 0];
                         sum_pixel1 += pt[((i + m)*cols + (j + n))*channels + 1];
                         sum_pixel2 += pt[((i + m)*cols + (j + n))*channels + 2];
                     }
                 img_out.ptr<uchar>(i - kernel_size)[(j - kernel_size)*channels + 0] = (double)sum_pixel / (double)(kernel_size*kernel_size * 4);
                 img_out.ptr<uchar>(i - kernel_size)[(j - kernel_size)*channels + 1] = (double)sum_pixel1 / (double)(kernel_size*kernel_size * 4);
                 img_out.ptr<uchar>(i - kernel_size)[(j - kernel_size)*channels + 2] = (double)sum_pixel2 / (double)(kernel_size*kernel_size * 4);
             }
         }
     }
  
 }
 ///中值滤波
 void lmt_median_blur(cv::Mat& img_in, cv::Mat& img_out, int kernel_size)
 {
     img_out = img_in.clone();
     cv::Mat mat1;
     cv::copyMakeBorder(img_in, mat1, kernel_size, kernel_size, kernel_size, kernel_size, cv::BORDER_REPLICATE);
  
     int cols = mat1.cols;
     int rows = mat1.rows;
     int channels = img_out.channels();
  
     cv::Mat mat[3];
     cv::Mat mat_out[3];
     cv::split(mat1, mat);
     cv::split(img_out, mat_out);
     for (int k = 0; k < 3; k++)
     {
         const uchar* const pt = mat[k].ptr<uchar>(0);
         uchar* pt_out = mat_out[k].ptr<uchar>(0);
         for (int i = kernel_size; i < rows - kernel_size; i++)
         {
             for (int j = kernel_size; j < cols - kernel_size; j++)
             {
                 long long int sum_pixel = 0;
                 int* K = new int[kernel_size*kernel_size * 4];
                 int ker_num = 0;
                 for (int m = -1 * kernel_size; m < kernel_size; m++)
                     for (int n = -1 * kernel_size; n < kernel_size; n++)
                     {
                         K[ker_num] = pt[(i + m)*cols + (j + n)];
                         ker_num++;
                     }
                 bubble_sort(K, ker_num);
                 mat_out[k].ptr<uchar>(i - kernel_size)[j - kernel_size] = K[ker_num / 2];
             }
         }
     }
     cv::merge(mat_out, 3, img_out);
 }
 ///高斯滤波
 void lmt_gaussian_blur(cv::Mat& img_src, cv::Mat& img_dst, cv::Size kernel_size)
 {
     img_dst = cv::Mat(cv::Size(img_src.cols, img_src.rows), img_src.type());
     int cols = img_src.cols;
     int rows = img_src.rows;
     int channels = img_src.channels();
     cv::Mat gassian_kernel = generate_gassian_kernel(0, 1, kernel_size);
     int width = kernel_size.width / 2;
     int height = kernel_size.height / 2;
     for (int i = height; i < rows - height; i++)
     {
         for (int j = width; j < cols - width; j++)
         {
             for (int k = 0; k < channels; k++)
             {
                 double sum = 0.0;
                 for (int m = -height; m <= height; m++)
                 {
                     for (int n = -width; n <= width; n++)
                     {
                         sum += (double)(img_src.ptr<uchar>(i + m)[(j + n)*channels + k]) * gassian_kernel.ptr<double>(height + m)[width + n];
                     }
                 }
                 if (sum > 255.0)
                     sum = 255;
                 if (sum < 0.0)
                     sum = 0;
                 img_dst.ptr<uchar>(i)[j*channels + k] = (uchar)sum;
             }
         }
     }
  
     
 }
 ///双边滤波
 void lmt_bilateral_filter(cv::Mat& img_in, cv::Mat& img_out, const int r, double sigma_d, double sigma_r)
 {
     int i, j, m, n, k;
     int nx = img_in.cols, ny = img_in.rows, m_nChannels = img_in.channels();
     const int w_filter = 2 * r + 1; // 滤波器边长  
  
     double gaussian_d_coeff = -0.5 / (sigma_d * sigma_d);
     double gaussian_r_coeff = -0.5 / (sigma_r * sigma_r);
     double  **d_metrix = new double *[w_filter];
     for (int i = 0; i < w_filter; ++i)
         d_metrix[i] = new double[w_filter];
     
     double r_metrix[256];  // similarity weight  
     img_out = cv::Mat(img_in.size(),img_in.type());
     uchar* m_imgData = img_in.ptr<uchar>(0);
     uchar* m_img_outData = img_out.ptr<uchar>(0);
     // copy the original image  
     double* img_tmp = new double[m_nChannels * nx * ny];
     for (i = 0; i < ny; i++)
         for (j = 0; j < nx; j++)
             for (k = 0; k < m_nChannels; k++)
             {
                 img_tmp[i * m_nChannels * nx + m_nChannels * j + k] = m_imgData[i * m_nChannels * nx + m_nChannels * j + k];
             }
  
     // compute spatial weight  
     for (i = -r; i <= r; i++)
         for (j = -r; j <= r; j++)
         {
             int x = j + r;
             int y = i + r;
  
             d_metrix[y][x] = exp((i * i + j * j) * gaussian_d_coeff);
         }
  
     // compute similarity weight  
     for (i = 0; i < 256; i++)
     {
         r_metrix[i] = exp(i * i * gaussian_r_coeff);
     }
  
     // bilateral filter  
     for (i = 0; i < ny; i++)
         for (j = 0; j < nx; j++)
         {
             for (k = 0; k < m_nChannels; k++)
             {
                 double weight_sum, pixcel_sum;
                 weight_sum = pixcel_sum = 0.0;
  
                 for (m = -r; m <= r; m++)
                     for (n = -r; n <= r; n++)
                     {
                         if (m*m + n*n > r*r) continue;
  
                         int x_tmp = j + n;
                         int y_tmp = i + m;
  
                         x_tmp = x_tmp < 0 ? 0 : x_tmp;
                         x_tmp = x_tmp > nx - 1 ? nx - 1 : x_tmp;   // 边界处理，replicate  
                         y_tmp = y_tmp < 0 ? 0 : y_tmp;
                         y_tmp = y_tmp > ny - 1 ? ny - 1 : y_tmp;
  
                         int pixcel_dif = (int)abs(img_tmp[y_tmp * m_nChannels * nx + m_nChannels * x_tmp + k] - img_tmp[i * m_nChannels * nx + m_nChannels * j + k]);
                         double weight_tmp = d_metrix[m + r][n + r] * r_metrix[pixcel_dif];  // 复合权重  
  
                         pixcel_sum += img_tmp[y_tmp * m_nChannels * nx + m_nChannels * x_tmp + k] * weight_tmp;
                         weight_sum += weight_tmp;
                     }
  
                 pixcel_sum = pixcel_sum / weight_sum;
                 m_img_outData[i * m_nChannels * nx + m_nChannels * j + k] = (uchar)pixcel_sum;
  
             } // 一个通道  
  
         } // END ALL LOOP  
     for (i = 0; i < w_filter; i++)
         delete[] d_metrix[i];
     delete[] d_metrix;
 }

 

5. Opencv API实现

opencv相关函数简介：

双边滤波函数：bilateralFilter(InputArray src, OutputArray dst, int d, double sigmaColor, double sigmaSpace,int borderType=BORDER_DEFAULT )

   src待滤波图像

   dst滤波后图像

   d滤波器半径

   sigmaColor滤波器值域的sigma

   sigmaSpace滤波器空间域的sigma

   borderType边缘填充方式 BORDER_REPLICATE BORDER_REFLECT BORDER_DEFAULT BORDER_REFLECT_101BORDER_TRANSPARENT BORDER_ISOLATED

 

均值滤波函数：blur(InputArray src, OutputArray dst, Size ksize, Point anchor=Point(-1,-1), intborderType=BORDER_DEFAULT );

   src待滤波图像

   dst滤波后图像

   ksize 均值滤波器的大小

Piont（-1，-1）指中心

   anchor均值滤波器的锚点也就是模板移动点

   borderType边缘填充方式 BORDER_REPLICATE BORDER_REFLECT BORDER_DEFAULT BORDER_REFLECT_101BORDER_TRANSPARENT BORDER_ISOLATED

 

高斯滤波函数：GaussianBlur(InputArray src, OutputArray dst, Size ksize, double sigmaX, double sigmaY=0,int borderType=BORDER_DEFAULT );

   src待滤波图像

   dst滤波后图像

   ksize 高斯滤波器的大小Size（x,y）,x和y必须是整数且是奇数。

   sigmaX 高斯滤波器的x方向的滤波器高斯sigma

   sigmaY 高斯滤波器的y方向的滤波器高斯sigma

   borderType边缘填充方式 BORDER_REPLICATE BORDER_REFLECT BORDER_DEFAULT BORDER_REFLECT_101BORDER_TRANSPARENT BORDER_ISOLATED

 

中值滤波函数：medianBlur(InputArray src, OutputArray dst, int ksize );

    src待滤波图像

    dst滤波后图像

    ksize 中值滤波器的大小

函数演示：

void bilateral_filter_show(void)
 {
     cv::Mat mat1 = cv::imread("F:\\CVlibrary\\obama.jpg", CV_LOAD_IMAGE_GRAYSCALE); //灰度图加载进来，BGR->HSV 然后取H参数
     if (mat1.empty())
         return;
     cv::imshow("原图像", mat1); 
     cv::Mat src = cv::imread("F:\\CVlibrary\\obama.jpg");
     cv::imshow("原始彩色图像", src);
     std::cout << "channel = " << mat1.channels() << std::endl;
     
     cv::Mat mat3;
     cv::bilateralFilter(src, mat3, 5, 50, 50,cv::BORDER_DEFAULT);
     cv::imshow("opencv给出的双边滤波器", mat3);
     cv::Mat mat4;
     cv::blur(src, mat4, cv::Size(3, 3));
     cv::imshow("均值滤波", mat4);
     cv::Mat mat5;
     cv::GaussianBlur(src, mat5, cv::Size(5, 5), 1,1);
     cv::imshow("高斯滤波器", mat5);
     cv::Mat mat6;
     cv::medianBlur(src, mat6, 3);
     cv::imshow("中值滤波", mat6); 
     cv::Mat mat7;
     lmt_gaussian_blur(src, mat7, cv::Size(5, 5));
     cv::imshow("my gaussian image",mat7);
  
     cv::waitKey(0);
 }

高斯、中值、均值、双边滤波的效果

#include "cv.h"
#include "highgui.h"
#include <iostream>

using namespace std;
using namespace cv;

int main(int argc, char* argv[])
{
        Mat src = imread("misaka.jpg");
        Mat dst;

        //参数是按顺序写的

        //高斯滤波
        //src:输入图像
        //dst:输出图像
        //Size(5,5)模板大小，为奇数
        //x方向方差
        //Y方向方差
        GaussianBlur(src,dst,Size(5,5),0,0);
        imwrite("gauss.jpg",dst);
        
        //中值滤波
        //src:输入图像
        //dst::输出图像
        //模板宽度，为奇数
        medianBlur(src,dst,3);
        imwrite("med.jpg",dst);
        
        //均值滤波
        //src:输入图像
        //dst:输出图像
        //模板大小
        //Point(-1,-1):被平滑点位置，为负值取核中心
        blur(src,dst,Size(3,3),Point(-1,-1));
        imwrite("mean.jpg",dst);

        //双边滤波
        //src:输入图像
        //dst:输入图像
        //滤波模板半径
        //颜色空间标准差
        //坐标空间标准差
        bilateralFilter(src,dst,5,10.0,2.0);//这里滤波没什么效果，不明白
        imwrite("bil.jpg",dst);

        waitKey();

        return 0;
}