空域滤波算法（均值，中值，高斯，双边） C++ +OPENCV

#include<bits/stdc++.h>
#include<opencv2/opencv.hpp>
#include<opencv2/highgui/highgui.hpp>
#include<time.h>
using namespace std;
using namespace cv;
typedef Mat mat;
typedef unsigned int uint;
typedef unsigned char uchar;
typedef long long ll;
const int maxn=1005;
const int N=1000;
ll a[maxn][maxn];
vector<double> val_mask;
int f[100];
/*********各种函数*******/
/*
mat img(3,3,CV_8UC1,Scalar(1));构造
mat img=(Mat_<uchar>(3,3)<<1,2,3,4,5,6,7,8,9);构造
cvtColor(src,dst,COLOR_BGR2GRAY);BGR变灰度
img.convertTo(dst,CV_32F)变成float数据
Scalar(127,0,255);设置三BGR通道的数值
addWeighted(src1,k1,src2,1-k1,0.0,dst);简单图像融合
blur(src,dst,Size(5,5),Point(-1,-1));系统均值滤波
GaussianBlur(src,dst,Size(5,5),11,11);系统高斯
medianBlur(src,dst,5);系统中值
filter2D(src,dst,-1,kernel,Point(-1,-1),0);卷积
*/
/**********高斯函数********/
double gaussian(double x,double deta)//高斯函数
{
    return exp(-0.5*(x*x)/(deta*deta));
}
/*********加椒盐噪声******/
mat add_salt_noise(mat pc,double k)
{
    int n=pc.rows,m=pc.cols;
    for(int i=0;i<n;i++)
        for(int j=0;j<m;j++)
        {
            double s=rand()%N/(double)(N);
            if(s<=k) pc.at<uchar>(i,j)=0;
        }
    return pc;
}
/********加高斯噪声*******/

double gaussnoise(double mu, double sigma)
{
    //定义小值，numeric_limits<double>::min()是函数,返回编译器允许的double型数最小值
    const double epsilon = std::numeric_limits<double>::min();
    static double z0, z1;
    static bool flag = false;
    flag = !flag;
    //flag为假构造高斯随机变量x
    if(!flag) return z1 * sigma + mu;
    //构造随机变量
    double u1, u2;
    do
    {
        u1 = rand() * (1.0 / RAND_MAX);
        u2 = rand() * (1.0 / RAND_MAX);
    } while (u1 <= epsilon);
    //flag为真构造高斯随机变量x
    z0 = sqrt(-2.0 * log(u1)) * cos(2 * CV_PI * u2);
    z1 = sqrt(-2.0 * log(u1)) * sin(2 * CV_PI * u2);
    return z0 * sigma + mu;
}
mat add_gauss_noise(mat pc,double mu,double sigma,double k)
{
    int n=pc.rows,m=pc.cols;
    for(int i=0;i<n;i++)
        for(int j=0;j<m;j++)
        {
            int val=pc.at<uchar>(i,j)+gaussnoise(mu,sigma)*k;
            val=min(val,255);
            val=max(val,0);
            pc.at<uchar>(i,j)=(uchar)val;
        }
    return pc;
}
/********两矩阵卷积********/
mat matrix_convolve(mat pc,mat mode)
{
    int n=pc.rows,m=pc.cols;
    int step=mode.rows/2;
    int modesum=0;
    for(int i=0;i<mode.rows;i++)
        for(int j=0;j<mode.cols;j++)
            modesum+=mode.at<int>(i,j);//计算算子和
    mat ans=pc.clone();
    for(int i=0;i<n;i++)
        for(int j=0;j<m;j++)
            if(i-step>=0&&i+step<n&&j-step>=0&&j+step<m)//将模板在图片游离
            {
                int sum=0;
                for(int r=i-step;r<=i+step;r++)
                    for(int c=j-step;c<=j+step;c++)
                        sum+=(int)pc.at<uchar>(r,c)*mode.at<int>(r-i+step,c-j+step);//对齐进行线性乘积
                ans.at<uchar>(i,j)=(uchar)(sum/modesum);//乘积和除算子和
            }
    return ans;
}
/*********均值滤波********/
mat get_ave(mat pc,int step)
{
    mat mode(step,step,CV_32SC1,Scalar(1));//填充一个全为1的宽度为step的mat
    return matrix_convolve(pc,mode);//复杂度O(pc.rows*pc.cols*step*step)
}
/*********均值滤波优化算法******/
mat fast_get_ave(mat pc,int step)
{
    int x=step/2,sum=step*step;
    int n=pc.rows,m=pc.cols;
    for(int i=0;i<n;i++)
        for(int j=0;j<m;j++)
        {
            a[i+1][j+1]=(ll)pc.at<uchar>(i,j)+a[i+1][j]+a[i][j+1]-a[i][j];//处理二维前缀和
            if(i-step+1>=0&&j-step+1>=0)
                pc.at<uchar>(i-x,j-x)=(uchar)((a[i+1][j+1]-a[i+1-step][j+1]-a[i+1][j+1-step]+a[i+1-step][j+1-step])/sum);
        }
    return pc;
}
/******高斯滤波（加权均值滤波）***/
mat get_gauss_ave(mat pc,int step)
{
    mat gs3=(Mat_<int>(3,3)<<1,2,1,2,4,2,1,2,1);
    mat gs5=(Mat_<int>(5,5)<<1,4,7,4,1,4,16,26,16,4,7,26,41,26,7,4,16,26,16,4,1,4,7,4,1);
    if(step==3) return matrix_convolve(pc,gs3);
    else if(step==5) return matrix_convolve(pc,gs5);
}
/*********中值滤波********/
mat get_mid(mat pc,int step)//pc表示待处理图片,step表示模板宽度
{
    step=step/2;
    mat ans=pc.clone();
    int n=pc.rows,m=pc.cols;
    for(int i=0;i<n;i++)
        for(int j=0;j<m;j++)
            if(i-step>=0&&i+step<n&&j-step>=0&&j+step<m)//模板游离
            {
                int k=0;
                for(int r=i-step;r<=i+step;r++)
                    for(int c=j-step;c<=j+step;c++)
                        f[k++]=pc.at<uchar>(r,c);//模板内所有值装到数组里
                sort(f,f+k);//排序
                ans.at<uchar>(i,j)=f[k/2];//取中值即可
            }
    return ans;
}
/*********双边滤波*****/
mat get_space_Mask(int step,double space_sigma)//空间域模板
{
    mat mask(step,step,CV_64FC1);
    int pos=step/2;
    for(int i=0;i<step;i++)
        for(int j=0;j<step;j++)
        {
            float x=abs(pos-i);
            float y=abs(pos-j);
            cout<<gaussian(sqrt(x*x+y*y),space_sigma)<<endl;
            mask.at<double>(i,j)=gaussian(sqrt(x*x+y*y),space_sigma);
        }
    cout<<mask<<endl;
    return mask;
}
void get_val_Mask(double val_sigma)//值域模板
{
    val_mask.clear();
    for(int i=0;i<256;i++)
        val_mask.push_back(gaussian(i,val_sigma));
}
mat double_filter(mat pc,int step,double spacesigma,double valsigma)
{
    mat space_mask=get_space_Mask(step,spacesigma);
    mat dst=pc.clone();
    mat new_mask(step,step,CV_64FC1);
    get_val_Mask(valsigma);
    step=step/2;
    int n=pc.rows,m=pc.cols;
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<m;j++)
        {
            if(i-step>=0&&i+step<n&&j-step>=0&&j+step<m)//将模板在图片游离
            {
                double sum=0;
                for(int r=i-step;r<=i+step;r++)
                    for(int c=j-step;c<=j+step;c++)
                    {
                        new_mask.at<double>(r-i+step,c-j+step)=space_mask.at<double>(r-i+step,c-j+step)*val_mask[abs((int)pc.at<uchar>(i,j)-(int)pc.at<uchar>(r,c))];
                        sum+=space_mask.at<double>(r-i+step,c-j+step);
                    }
                new_mask=new_mask/sum;
                sum=0;
                double all=0;
                for(int r=i-step;r<=i+step;r++)
                    for(int c=j-step;c<=j+step;c++)
                    {
                        sum+=new_mask.at<double>(r-i+step,c-j+step);
                        all+=new_mask.at<double>(r-i+step,c-j+step)*(double)pc.at<uchar>(r,c);
                    }
                all=all/sum;
                int c=(int)all;
                dst.at<uchar>(i,j)=c;
            }
        }
    }
    return dst;
}
/*********去除旁瓣效应****/
mat bounary_get_mid(mat pc,int step)//pc表示待处理图片,step表示模板宽度
{
    step=step/2;
    mat ans=pc.clone();
    int n=pc.rows,m=pc.cols;
    for(int i=0;i<n;i++)
        for(int j=0;j<m;j++)
            {
                int k=0;
                for(int r=i-step;r<=i+step;r++)
                    for(int c=j-step;c<=j+step;c++)
                    {
                        int x=r,y=c;
                        if(x<0) x=-x;
                        else if(x>=n) x=2*n-x-2;
                        if(y<0) y=-y;
                        else if(y>=m) y=2*m-y-2;//对超出边界的进行变换即可
                        f[k++]=pc.at<uchar>(x,y);//模板内所有值装到数组里
                    }
                sort(f,f+k);//排序
                ans.at<uchar>(i,j)=f[k/2];//取中值即可
            }
    return ans;
}
int main()
{
//    srand(time(NULL));
    mat src = imread("E://PC//6.jpg",0);
    imshow("src",src);
    mat ans;
//    mat pc=add_salt_noise(src,0.1);imshow("add_salt",pc);//加入椒盐噪声
    mat pc=add_gauss_noise(src,4,2,2);//435
    imshow("add_gauss_noise",pc);//加入高斯噪声
//    imshow("mid_filter_3*3",get_mid(pc,3));//中值滤波
//    medianBlur(pc,ans,3);//系统中值滤波
//    imshow("mid_function_3*3",ans);
//    int T=clock();
//    ans=get_ave(src,10);//均值滤波
//    cout<<"10*10 ordinary ave_filter need time:"<<clock()-T<<" ms"<<endl;
//    T=clock();
//    ans=fast_get_ave(src,10);
//    cout<<"10*10 fast ave_filter need:"<<clock()-T<<" ms"<<endl;
//    blur(pc,ans,Size(3,3),Point(-1,-1));//系统均值滤波
//    imshow("ave_function_3*3",ans);
//    imshow("gauss_filter_3*3",get_gauss_ave(src,3));//高斯滤波
    GaussianBlur(pc,ans,Size(3,3),1,3);//系统高斯滤波
    imshow("gauss_function_3*3",ans);
    imshow("double_filter_3*3",double_filter(pc,3,2,15));//双边滤波//3 15 15
//    imshow("mid",get_mid(pc,5));
//    imshow("remove_bounary_noise_mid",bounary_get_mid(pc,5));
    waitKey();
    return 0;
}