本章主要分为两大内容：
第一，针对只存在噪声影响下的图像复原；
第二，针对存在退化函数和噪声的两种作用下的图像复原；
5.1图像复原的概念以及模型
1、概念
与图像增强类似，图像复原的目的也是改善给定的图像，但是图像增强是一个主观的过程，而图像复原是一个客观的过程。
复原技术是面向退化模型的，并且采用相反的过程进行处理，以便恢复出原图像 。
1
2
2、图像退化/复原模型
1、 退化复原模型如下：
2、上述模型在空间域中表示为：
五角星表示为空间域中卷积，空间域中在卷积=频率域中在乘积，故在频率域中在乘积为
3、在频率域中表示为：
5.2只有加性噪声下在图像复原
5.2.1噪声模型
首先噪声的估计是取一小段直方图，来对比看是哪一种噪声在直方图，然后就根据公式计算该噪声的均布直方图；
加性噪声，加性噪声和图像信号强度不相关，这类噪声可以看着理想无噪声图像f和噪声的和。
1.高斯噪声
也称为正态噪声，源于电子电路噪声污染、低照明度和高温带来的传感器污染
概率密度函数（PDF）：
其中，z表示灰度值，z一杠表示灰度值均值， 点他表示z的标准差，其平方为z的方差；
高斯函数曲线如下所示：
2、瑞丽噪声
适用于近似歪斜在直方图十分适用；有助于在深度成像中表征噪声现象
其概率密度公式如下：
瑞丽密度曲线：
注意距离原点在位移和密度在基本形状向右变形
3、爱尔兰（伽马）噪声
应用于激光成像中
概率密度函数：
此公式被称为爱尔兰密度
4、指数噪声
应用于激光成像中
是伽马噪声中b=1在特殊情况
概率密度函数为：
5、均匀噪声
仿真中使用许多随机数生成器在基础是非常有用的
概率密度函数为：
6、脉冲（椒盐）噪声
概率密度函数为：
注意：
如果b>a，则灰度级b在图像中表现为一个亮点，反之，灰度级a表现为一个暗点；
如果Pa或Pb为零，则脉冲称为单级脉冲；
如果Pa\ Pb都不为零，尤其是两者近似相等的时候，脉冲噪声值将类似于在图像上随机分布在胡椒和盐粉微粒，故双击脉冲噪声又称为椒盐噪声；
正脉冲为数字图像在灰度最大值255为白点，负脉冲为灰度最小值0为黑点
6、周期噪声（本章唯一的空间相关噪声）
一幅图像的周期噪声是在图像获取期间由电力或机电干扰产生的；
周期噪声可以通过频率域滤波进行显著在减少；
5.2.2噪声参数估计
首先噪声的估计是取一小段直方图，来对比看是哪一种噪声在直方图，然后就根据公式计算该噪声的均布直方图；
1、高斯噪声的以c++实现：
#include<opencv2/opencv.hpp>
#include<iostream>
#include <limits>
using namespace std;
using namespace cv;
//高斯噪声实现
//生成高斯噪声
double generateGaussianNoise(double mu, double sigma)
{
//定义小值
const double epsilon = numeric_limits<double>::min();//是函数，返回编译器允许的 double 型数 最小值。
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);//伪随机数生成函数 rand 所能返回的最大数值。
u2 = rand() * (1.0 / RAND_MAX);
} while (u1 <= epsilon);
//flag为真构造高斯随机变量
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 addGaussianNoise(Mat &srcImag)
{
Mat dstImage = srcImag.clone();
int channels = dstImage.channels();
int rowsNumber = dstImage.rows;
int colsNumber = dstImage.cols*channels;
//判断图像的连续性
if (dstImage.isContinuous())
{
colsNumber *= rowsNumber;
rowsNumber = 1;
}
for (int i = 0; i < rowsNumber; i++)
{
for (int j = 0; j < colsNumber; j++)
{
//添加高斯噪声
int val = dstImage.ptr<uchar>(i)[j] + generateGaussianNoise(0, 2.235) * 32;
if (val < 0)
val = 0;
if (val > 255)
val = 255;
dstImage.ptr<uchar>(i)[j] = (uchar)val;
}
}
return dstImage;
}
int main()
{
Mat src = imread("C:/Users/征途/Desktop/vs-cpp/第五章/01.jpg");
if (!src.data)
{
cout << "输入有误" << endl;
return -1;
}
imshow("src", src);
Mat dst = addGaussianNoise(src);
imshow("dst", dst);
waitKey(0);
return 0;
}
2、椒盐噪声实现代码
#include<opencv2/opencv.hpp>
#include<iostream>
#include <limits>
using namespace std;
using namespace cv;
//椒盐噪声
Mat addSaltNoise(const Mat srcImage, int n)//n为个数
{
Mat dstImage = srcImage.clone();
for (int k = 0; k < n; k++)
{
//随机取值行列
int i = rand() % dstImage.rows;
int j = rand() % dstImage.cols;
//图像通道判定
if (dstImage.channels() == 1)
{
dstImage.at<uchar>(i, j) = 255; //盐噪声
}
else
{
dstImage.at<Vec3b>(i, j)[0] = 255;
dstImage.at<Vec3b>(i, j)[1] = 255;
dstImage.at<Vec3b>(i, j)[2] = 255;
}
}
for (int k = 0; k < n; k++)
{
//随机取值行列
int i = rand() % dstImage.rows;
int j = rand() % dstImage.cols;
//图像通道判定
if (dstImage.channels() == 1)
{
dstImage.at<uchar>(i, j) = 0; //椒噪声
}
else
{
dstImage.at<Vec3b>(i, j)[0] = 0;
dstImage.at<Vec3b>(i, j)[1] = 0;
dstImage.at<Vec3b>(i, j)[2] = 0;
}
}
return dstImage;
}
int main()
{
Mat src = imread("C:/Users/征途/Desktop/vs-cpp/第五章/01.jpg");
if (!src.data)
{
cout << "输入有误" << endl;
return -1;
}
imshow("src", src);
Mat dst = addGaussianNoise(src);
Mat dst2 = addSaltNoise(src, 5000);
imshow("dst", dst);
imshow("dst2", dst2);
waitKey(0);
return 0;
}
5.2.3只存在噪声的复原–空间滤波
(1) 如果不考虑退化函数，图像退化模型就简化为图像噪声模型
只存在加性噪声情况下：图像增强问题成为单纯的图像去噪问题，可以通过空间域滤波等众多方法解决
空间域中常见的滤波器分为均值/统计排序/自适应滤波器.
1.均值滤波器:
算术/几何/谐波/逆谐波均值滤波器.
2.统计排序滤波器:
中值/最大值/最小值/中点滤波器.
3.自适应滤波器:
自适应均值/自适应中值滤波器
1、均值滤波器
（1）、算数均值滤波器实现加性噪声下的图像复原
计算速度快，但是有点模糊 ，传统意义上的均值滤波
令S(xy)表示中心在点（x,y）处，大小为mxn的矩形子图像窗口的一组坐标；
算术均值滤波器在S(xy)区域中计算被污染图像g(xy)的平均值
#include<opencv2/opencv.hpp>
#include<iostream>
#include <limits>
#include<math.h>
using namespace std;
using namespace cv;
//算数均值滤波
void ArithAverFilter(Mat &src, Mat &dst)
{
if (src.channels() == 3)//彩色图像
{
for (int i = 1; i < src.rows; i++)
{
for (int j = 1; j < src.cols; j++)
{
if ((i - 1 >= 0) && (j - 1) >= 0 && (i + 1) < src.rows && (j + 1) < src.cols)
{
//边缘不进行处理
dst.at<Vec3b>(i, j)[0] = (src.at<Vec3b>(i, j)[0] + src.at<Vec3b>(i - 1, j - 1)[0] + src.at<Vec3b>(i - 1, j)[0] + src.at<Vec3b>(i, j - 1)[0] +
src.at<Vec3b>(i - 1, j + 1)[0] + src.at<Vec3b>(i + 1, j - 1)[0] + src.at<Vec3b>(i + 1, j + 1)[0] + src.at<Vec3b>(i, j + 1)[0] +
src.at<Vec3b>(i + 1, j)[0]) / 9;
dst.at<Vec3b>(i, j)[1] = (src.at<Vec3b>(i, j)[1] + src.at<Vec3b>(i - 1, j - 1)[1] + src.at<Vec3b>(i - 1, j)[1] + src.at<Vec3b>(i, j - 1)[1] +
src.at<Vec3b>(i - 1, j + 1)[1] + src.at<Vec3b>(i + 1, j - 1)[1] + src.at<Vec3b>(i + 1, j + 1)[1] + src.at<Vec3b>(i, j + 1)[1] +
src.at<Vec3b>(i + 1, j)[1]) / 9;
dst.at<Vec3b>(i, j)[2] = (src.at<Vec3b>(i, j)[2] + src.at<Vec3b>(i - 1, j - 1)[2] + src.at<Vec3b>(i - 1, j)[2] + src.at<Vec3b>(i, j - 1)[2] +
src.at<Vec3b>(i - 1, j + 1)[2] + src.at<Vec3b>(i + 1, j - 1)[2] + src.at<Vec3b>(i + 1, j + 1)[2] + src.at<Vec3b>(i, j + 1)[2] +
src.at<Vec3b>(i + 1, j)[2]) / 9;
}
else
{
//边缘赋值
dst.at<Vec3b>(i, j)[0] = src.at<Vec3b>(i, j)[0];
dst.at<Vec3b>(i, j)[1] = src.at<Vec3b>(i, j)[1];
dst.at<Vec3b>(i, j)[2] = src.at<Vec3b>(i, j)[2];
}
}
}
}
if (src.channels() == 1)
{
//灰度图像
for (int i = 0; i < src.rows; i++)
{
for (int j = 0; j < src.cols; j++)
{
if ((i - 1 >= 0) && (j - 1) >= 0 && (i + 1) < src.rows && (j + 1) < src.cols)
{
//边缘不进行处理
dst.at<uchar>(i, j) = (src.at<uchar>(i, j) + src.at<uchar>(i - 1, j - 1) + src.at<uchar>(i - 1, j) + src.at<uchar>(i, j - 1) +
src.at<uchar>(i - 1, j + 1) + src.at<uchar>(i + 1, j - 1) + src.at<uchar>(i + 1, j + 1) + src.at<uchar>(i, j + 1) +
src.at<uchar>(i + 1, j)) / 9;
}
else
{
//边缘赋值
dst.at<uchar>(i, j) = src.at<uchar>(i, j);
}
}
}
}
imshow("arithAverFilter", dst);
}
int main()
{
Mat src = imread("C:/Users/征途/Desktop/vs-cpp/第五章/02.jpg");
if (!src.data)
{
cout << "输入有误" << endl;
return -1;
}
Mat dst;
imshow("src", src);
ArithAverFilter(src,dst);
imshow("dst", dst);
//imshow("dst2", dst2);
waitKey(0);
return 0;
}
（2）、几何均值滤波器实现加性噪声下图像复原
图像模糊程度低
对子窗口的元素相乘（需要判断是否为零），并对乘积求1/m*n次幂
实现代码：
#include<opencv2/opencv.hpp>
#include<iostream>
#include <limits>
#include<math.h>
#include<vector>
using namespace std;
using namespace cv;
void GeoAverFliter(const Mat &src, Mat &dst)
{
Mat _dst(src.size(), CV_32FC1);
double power = 1.0 / 9;
cout << "power:" << power << endl;
double geo = 1;
for (int i = 0; i < src.rows; i++)
{
for (int j = 0; j < src.cols; j++)
{
if ((i - 1) > 0 && (i + 1) < src.rows && (j - 1) > 0 && (j + 1) < src.cols)
{
//下面应用公式
if (src.at<uchar>(i, j) != 0) geo = geo * src.at<uchar>(i, j);
if (src.at<uchar>(i, j+1) != 0) geo = geo * src.at<uchar>(i, j+1);
if (src.at<uchar>(i, j-1) != 0) geo = geo * src.at<uchar>(i, j-1);
if (src.at<uchar>(i+1, j) != 0) geo = geo * src.at<uchar>(i+1, j);
if (src.at<uchar>(i-1, j) != 0) geo = geo * src.at<uchar>(i-1, j);
if (src.at<uchar>(i+1, j+1) != 0) geo = geo * src.at<uchar>(i+1, j+1);
if (src.at<uchar>(i-1, j-1) != 0) geo = geo * src.at<uchar>(i-1, j-1);
if (src.at<uchar>(i+1, j-1) != 0) geo = geo * src.at<uchar>(i+1, j-1);
if (src.at<uchar>(i-1, j+1) != 0) geo = geo * src.at<uchar>(i-1, j+1);
_dst.at<uchar>(i, j) = pow(geo, power);
}
else
{
dst.at<float>(i, j) = src.at<uchar>(i, j);
}
}
}
_dst.convertTo(dst, CV_8UC1);
}
int main()
{
Mat src = imread("C:/Users/征途/Desktop/vs-cpp/第五章/02.jpg");
if (!src.data)
{
cout << "输入有误" << endl;
return -1;
}
Mat dst1,dst2;
imshow("src", src);
ArithAverFilter(src,dst1);
imshow("dst1", dst1);
GeoAverFliter(src, dst2);
imshow("dst2", dst2);
waitKey(0);
return 0;
}
（3）、谐波均值滤波
对于盐粒噪声效果好，不适合胡椒噪声，善于处理像高斯噪声那样的其他噪声；
//谐波均值滤波器
void xieboFliter(const Mat &src, Mat &dst)
{
Mat dst4 (src.size(), CV_32FC1);
double mn = 9.0;
double he = 0;
for (int i = 0; i < src.rows; i++)
{
for (int j = 0; j < src.cols; j++)
{
if ((i - 1) > 0 && (i + 1) < src.rows && (j - 1) > 0 && (j + 1) < src.cols)
{
he = (1.0 / src.at<uchar>(i, j)) + (1.0 / src.at<uchar>(i, j + 1)) + (1.0 / src.at<uchar>(i, j - 1)) + (1.0 / src.at<uchar>(i + 1, j)) + (1.0 / src.at<uchar>(i - 1, j))
+ (1.0 / src.at<uchar>(i + 1, j + 1)) + (1.0 / src.at<uchar>(i + 1, j - 1)) + (1.0 / src.at<uchar>(i - 1, j + 1)) + (1.0 / src.at<uchar>(i - 1, j - 1));
dst4 = mn / he;
}
}
}
dst4.convertTo(dst, CV_8UC1);
}
2、统计排序滤波器
先对卷积空间的像素点排序,并分别取中值/最大值/最小值/中点.
（1）中值滤波器
使用像素领域中的灰度级的中值替代该像素点的值
针对单极脉冲和双极脉冲脉冲，即椒盐噪声特别有效。
（2）最大值和最小值滤波器
最大值滤波器对发现图像中的亮点非常有效，应用于胡椒噪声；
最小值滤波器对发现图像中的最暗点非常有效，应用于盐粒噪声；
#include<opencv2/opencv.hpp>
#include<iostream>
#include <limits>
#include<math.h>
#include<vector>
using namespace std;
using namespace cv;
//uchar max(int *arr)
uchar max1(uchar n1, uchar n2, uchar n3, uchar n4, uchar n5,
uchar n6, uchar n7, uchar n8, uchar n9)
{
int arr[9] = { n1 ,n2, n3, n4, n5, n6, n7, n8, n9 };
int max = arr[0];
for (int i = 1; i < 9; i++)
{
max = (max > arr[i]) ? max : arr[i];
}
return max;
}
void maxfilter( Mat &src, Mat &dst)
{
Mat dst0(src.size(), CV_32FC1);
for (int i = 0; i < src.rows; i++)
{
for (int j = 0; j < src.cols; j++)
{
if ((i - 1) > 0 && (i + 1) < src.rows && (j - 1) > 0 && (j + 1) < src.cols)
{
//int arr[9] = { src.at<uchar>(i, j), src.at<uchar>(i, j + 1), src.at<uchar>(i, j - 1), src.at<uchar>(i + 1, j), src.at<uchar>(i - 1, j), src.at<uchar>(i + 1, j + 1), src.at<uchar>(i + 1, j - 1), src.at<uchar>(i - 1, j + 1), src.at<uchar>(i - 1, j - 1) };
//dst0.at<uchar>(i, j) = max1(&arr[9]);
dst0.at<uchar>(i, j) = max1(src.at<uchar>(i, j), src.at<uchar>(i, j + 1), src.at<uchar>(i, j - 1), src.at<uchar>(i + 1, j), src.at<uchar>(i - 1, j), src.at<uchar>(i + 1, j + 1), src.at<uchar>(i + 1, j - 1), src.at<uchar>(i - 1, j + 1), src.at<uchar>(i - 1, j - 1));
}
else
{
dst0.at<uchar>(i, j) = src.at<uchar>(i, j);
}
}
}
dst0.convertTo(dst, CV_8UC1);
}
int main()
{
Mat src, dst1;
src = imread("C:/Users/征途/Desktop/vs-cpp/第五章/05.jpg");
if (!src.data)
{
cout << "输入有误" << endl;
return -1;
}
cvtColor(src, src, CV_BGR2GRAY);
imshow("s", src);
maxfilter(src, dst1);
imshow("dst1", dst1);
waitKey(0);
return 0;
}
（3）中点滤波器
适用于混合噪声，如高斯噪声和椒盐噪声混合的情况。
计算滤波器包围区域中最大值最小值中间的中点；
对 随机分布噪声（高斯噪声、均匀噪声）非常有效
#include<opencv2/opencv.hpp>
#include<iostream>
#include <limits>
#include<math.h>
#include<vector>
using namespace std;
using namespace cv;
int mysort(int *arr)
//uchar mysort(uchar n1, uchar n2, uchar n3, uchar n4, uchar n5,
//uchar n6, uchar n7, uchar n8, uchar n9)
{
bool flag = true;
for (int i = 0; i < 9 ; i++)
{
while (flag = false)
{
for (int j = 0; j < 9 - i - 1; j++)
{
if (arr[j] > arr[j + 1])
{
int temp = arr[j];
arr[j] = arr[j + 1];
arr[j + 1] = temp;
}
}
}
}
int a = (arr[0] + arr[8]) / 2;
return a;//返回中值
}
//单通道
void myzhongdainBlur(Mat &src, Mat&dst)
{
Mat dst0(src.size(), src.type());
for (int i = 0; i < src.rows; i++)
{
for (int j = 0; j < src.cols; j++)
{
if ((i - 1) > 0 && (i + 1) < src.rows && (j - 1) > 0 && (j + 1) < src.cols)
{
int arr[9] = { src.at<uchar>(i, j), src.at<uchar>(i, j + 1), src.at<uchar>(i, j - 1), src.at<uchar>(i + 1, j), src.at<uchar>(i - 1, j), src.at<uchar>(i + 1, j + 1), src.at<uchar>(i + 1, j - 1), src.at<uchar>(i - 1, j + 1), src.at<uchar>(i - 1, j - 1) };
dst0.at<uchar>(i, j) = mysort(&arr[9]);
}
else
{
dst0.at<uchar>(i, j) = src.at<uchar>(i, j);
}
}
}
dst0.convertTo(dst, CV_8UC1, 1, 0);
}
int main()
{
Mat src, dst1;
src = imread("C:/Users/征途/Desktop/vs-cpp/第五章/05.jpg");
if (!src.data)
{
cout << "输入有误" << endl;
return -1;
}
cvtColor(src, src, CV_BGR2GRAY);
imshow("s", src);
myzhongdainBlur(src, dst1);
namedWindow("dst1", WINDOW_NORMAL);
imshow("dst1", dst1);
cvWaitKey(0);
return 0;
}
（4）修正的阿尔法均值滤波器
假设在领域S(xy)中去除g(s，t)最低灰度值的d/2和最高灰度值的d/2，令剩下的gr(s,t)代表剩下的mn-d个像素。由这些剩余像素的平均值形成的滤波器称为修正阿尔法滤波器：
d属于0到mn-1
如果d=0,则转化为算数均值滤波器；
如果d=mn-1，则转化为中值滤波器；``
3、自适应滤波器
（1）、自适应局部降低噪声滤波器
自适应滤波器的行为变化基于由m*n矩形窗口Sxy定义的区域内图像的统计特性，它的性能要明显优于前面介绍的滤波器，代价是滤波器的复杂度。
5.2.4利用频率域滤波消除周期噪声
主要有带阻滤波器、带通滤波器、陷波滤波器
消除对象是周期性噪声
1、带阻滤波器
带阻滤波器为Hbr(u,v)
主要有理想低通高通滤波器、巴特沃斯低通高通滤波器、高斯低通高通滤波器
这几个滤波器前面已经实现过，在此不再赘述
2、带通滤波器
带通滤波器与带阻滤波器正好相反，为Hbp（u，v）;
传递函数为：
Hbp(u,v) = 1-Hbr(u,v)
3、陷波滤波器
陷波滤波器是更有用的选择性滤波器；
陷波滤波器拒绝（或通过）事先定义的关于频率矩形中心的一个领域的频率。
陷波滤波器属于零相移滤波器，必须关于原点对称，因此一个中心位于陷波滤波器（u0，v0)的陷波在位置（-u0,v0)必须有一个对应的陷波
1、陷波带阻滤波器
陷波带阻滤波器可以用中心已被平移到陷波滤波器中心的高通滤波器的乘积构造。一般形式为：
上式中其中，Hk(u,v);H-k(u,v)是高通滤波器，他们的中心分别位于（u0,v0）,(-u0,-v0)中，这些中心是根据频率矩形的中心(M/2,N/2）来确定的。对于每个滤波器，距离的计算是由下式计算：
对于n阶巴特沃斯带阻滤波器，它包含三个陷波对：
2、陷波带通滤波器
5.3、即存在退化函数又有噪声干扰情况下的图像复原
在图像拍摄过程中由于各种原因会造成图像退化，图像退化模型如下：
其中，⋆ \star⋆为卷积符号，f ( x , y ) f(x,y)f(x,y)为输入图像，g ( x , y ) g(x,y)g(x,y)为退化图像，h ( x , y ) h(x,y)h(x,y)为退化函数，η ( x , y ) \eta(x,y)η(x,y)为加性噪声，将上式进行傅里叶变换有：
1、估计退化函数
估计退化函数的三种方法：1、观察法 2、试验法、 3、数字建模法
1、图像观察法
假设提供了一幅退化图像，而没有退化函数H HH的知识，那么估计该函数的一个方法就是收集图像自身的信息。例如，如果图像是模糊的，那么观察包含简单结构的一小部分图像，像某一物体和背景的一部分。为了减少观察时的噪声影响，可以寻找最强信号内容区。
其中
2、试验估计法
如果可以使用与获取退化图像的设备相似的装置，理论上得到一个准确的退化估计是可能的。与退化图像类似的图像可以通过各种设置得到，退化这些图像使其尽可能接近希望复原的图像。利用相同的系统设置(退化函数),由成像一个脉冲(小亮点)得到退化的冲激响应。
实际上就是利用产生图像相似的装置，进行设置参数从而得到一个图像，这个图像类似于退化图像，然后由这个退化图像得到原图的退化函数；
其中，G是观察图像的傅里叶变换，A是一个常量，表示冲激强度。
3、建模型估计
由于退化模型可解决图像复原问题，因此多年来一直在应用。在某些情况下，模型要把引起退化的环境因素考虑在内。
建模估计中可以通过运动数学模型将退化函数构造为：
大气湍流模型：
2、逆滤波实现图像复原
逆滤波是由退化函数直接图像复原的最初手段，也是最简单的手段；
在该方法中，用退化函数除退化图像的傅立叶变换来计算原始图像的傅立叶变换估计：
传递函数为：
不考虑噪声：
故可以反变换求：
考虑噪声影响：
需要注意的是：
3、最小均方误差滤波(维纳滤波)
维纳滤波是一种建立在最小化统计准则的基础上的复原方法，在平均意义上，它可以看成是最优的。
维纳滤波是在频域中处理图像的一种算法，是一种非常经典的图像增强算法，不仅可以进行图像降噪，还可以消除由于运动等原因带来的图像模糊。
该方法建立在退化函数和噪声统计的基础上，考虑了两者的共同作用，目标是为了找到未污染图像f的一个估计值f，使他们之间的均方误差最小。
维纳滤波综合了退化函数和噪声统计特征两个方面进行复原处理，在认为图像和噪声是随机过程的基础上，以恢复图像和原图像的均方误差最小为准则。
其中，上面的E是参数的期望值。
假定：噪声和图像不相关；其中一个有零均值(一般假设噪声是零均值)；估计的灰度级是退化图像灰度级的线性函数；维纳滤波的公式如下：
维纳滤波：一个复数量与其共轭的乘积等于该复数量幅度的平方，这个结果就是维纳滤波。
维纳滤波中没有退化函数H为0的情况，除非对于相同的u值和v值，整个分母都为0；
下面的例子中我们还对退化函数进行了简化，将退化函数置为1，因此维纳滤波公式简化为：
代码实现中：退化函数为1时：
实现代码：
#include<iostream>
#include<opencv2/opencv.hpp>
#include<string>
using namespace cv;
using namespace std;
//获得图像频谱的函数写法
Mat GetSpectrum(const Mat &src)
{
Mat dst, cpx;
Mat planes[] = { Mat_<float>(src), Mat::zeros(src.size(), CV_32F) };
merge(planes, 2, cpx);//合并通道
dft(cpx, cpx);//傅里叶变换
split(cpx, planes);//分离通道，实部re=planes[0]虚部im = planes[1]
magnitude(planes[0], planes[1], dst);//计算二维矢量的幅值函数，//求幅度图
//频谱就是频域幅度图的平方
multiply(dst, dst, dst);//频谱就是幅度图的平方
return dst;
}
//添加噪声
double generateGaussianNoise(double mu, double sigma)
{
//定义小值
const double epsilon = numeric_limits<double>::min();//是函数，返回编译器允许的 double 型数 最小值。
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);//伪随机数生成函数 rand 所能返回的最大数值。
u2 = rand() * (1.0 / RAND_MAX);
} while (u1 <= epsilon);
//flag为真构造高斯随机变量
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 addGaussianNoise(Mat &srcImag)
{
Mat dstImage = srcImag.clone();
int channels = dstImage.channels();
int rowsNumber = dstImage.rows;
int colsNumber = dstImage.cols*channels;
//判断图像的连续性
if (dstImage.isContinuous())
{
colsNumber *= rowsNumber;
rowsNumber = 1;
}
for (int i = 0; i < rowsNumber; i++)
{
for (int j = 0; j < colsNumber; j++)
{
//添加高斯噪声
int val = dstImage.ptr<uchar>(i)[j] + generateGaussianNoise(0, 2.235) * 32;
if (val < 0)
val = 0;
if (val > 255)
val = 255;
dstImage.ptr<uchar>(i)[j] = (uchar)val;
}
}
return dstImage;
}
//维纳滤波的实现
Mat wienerfilter(const Mat &src, const Mat &ref, int stddev)
{
//这些图片是过程中会用到的，pad是原图像0填充后的图像，cpx是双通道频域图，mag是频域幅值图，dst是滤波后的图像
Mat pad, cpx, dst;
//获取傅里叶变化最佳图片尺寸，为2的指数
int m = getOptimalDFTSize(src.rows);
int n = getOptimalDFTSize(src.cols);
//对原始图片用0进行填充获得最佳尺寸图片
copyMakeBorder(src, pad, 0, m - src.rows, 0, n - src.cols, BORDER_CONSTANT, Scalar::all(0));
//获得参考图片频谱
Mat tmpR(pad.rows, pad.cols, CV_8U);
//resize(ref, tmpR, tmpR.size()); //重新规划大小
Mat refSpectrum = GetSpectrum(tmpR);//频谱
//获得噪声频谱
Mat tmpN(pad.rows, pad.cols, CV_32F);
//randn(tmpN, Scalar::all(0), Scalar::all(stddev));//randn函数常用来产生高斯白噪声信号,生成m×n随机矩阵，其元素服从均值为0，方差为stddev的标准正态分布。
addGaussianNoise(tmpN);//生成高斯噪声
Mat noiseSpectrum = GetSpectrum(tmpN);
//对src进行傅里叶变换
Mat planes[] = { Mat_<float>(pad), Mat::zeros(pad.size(), CV_32F) };
merge(planes, 2, cpx);//合并两通道然后输出
dft(cpx, cpx);
split(cpx, planes);
//维纳滤波因子
Mat factor = refSpectrum / (refSpectrum + noiseSpectrum);
//cpx = factor * cpx;
multiply(planes[0], factor, planes[0]);
multiply(planes[1], factor, planes[1]);
// 重新合并实部planes[0]和虚部planes[1]
merge(planes, 2, cpx);
//进行反傅里叶变换
idft(cpx, dst, DFT_SCALE | DFT_REAL_OUTPUT);//idft傅里叶反变换
dst.convertTo(ref, CV_8UC1);
return ref;
}
int main()
{
Mat src, dst;
src = imread("C:/Users/征途/Desktop/vs-cpp/第五章/05.jpg");
if (!src.data)
{
cout << "输入错误" << endl;
}
imshow("src", src);
wienerfilter(src, dst, 1);
imshow("dst", dst);
waitKey(0);
return 0;
}
文章知识点与官方知识档案匹配，可进一步学习相关知识
————————————————
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
原文链接：https://blog.csdn.net/wangjia2575525474/article/details/117281214
章目录
前言
一、原理
1.概率密度函数及示意图灰度分布表示曲线
2.高斯随机数的产生
二、代码
前言
数字图像处理c++ opencv（VS2019 opencv4.53）持续更新
一、原理
1.概率密度函数及示意图灰度分布表示曲线
2.高斯随机数的产生
Box-Muller 算法：先得到服从均匀分布的随机数，再将服从均匀分布的随机数转变为服从正态分布。
设U1和U2是服从均匀分布的随机数，则通过下列公式
或者
得到正态分布随机数Z，然后通过均值和标准差计算对应的高斯随机数:
二、代码
代码如下：
#include<iostream>
#include<opencv2/opencv.hpp>
#include <math.h>
#include <time.h>
#define PI 3.14159
using namespace cv;
using namespace std;
void Add_GaussianNoise(Mat img_input, Mat& img_output, double mean, double sigma, int k);
int main()
{
Mat image, image_gray, image_output, image_output2; //定义输入图像，灰度图像，输出图像,
image = imread("lena.png"); //读取图像；
if (image.empty())
{
cout << "读取错误" << endl;
return -1;
}
imshow("image", image);
//转换为灰度图像
cvtColor(image, image_gray, COLOR_BGR2GRAY);
imshow("image_gray", image_gray);
//1、利用opencv添加高斯噪声
Mat noise = Mat::zeros(image_gray.size(), image_gray.type());
randn(noise, 1, 15); //randn(a,b,c)产生高斯噪声，其中a为输出矩阵，b为均值，c为方差
image_output = Mat::zeros(image_gray.size(), image_gray.type());
add(image_gray, noise, image_output, Mat(), -1); //将灰度图与噪声矩阵相加得到噪声图像
imshow("image_output", image_output);
//2、自己实现高斯噪声的添加
Add_GaussianNoise(image_gray, image_output2, 1, 15, 1);//输入图像，输出图像，均值，方差，系数
imshow("image_output2", image_output2);
waitKey(0); //暂停，保持图像显示，等待按键结束
return 0;
}
//产生高斯随机数
double generateGaussianNoise(double mu, double sigma)
{
static double V1, V2, S;
static int phase = 0;
double X;
double U1, U2;
if (phase == 0)
{
U1 = (double)rand() / RAND_MAX; //随机数1
U2 = (double)rand() / RAND_MAX; //随机数2
V1 = sqrt((-2) * log(U1));
V2 = 2 * PI * U2;
S = V1 * cos(V2); //通过算法实现高斯随机数：
}
return mu + sigma * S;
}
//噪声图像
void Add_GaussianNoise(Mat img_input, Mat& img_output, double mean, double sigma, int k)
{
img_output.create(img_input.rows, img_input.cols, img_input.type());
for (int x = 0; x < img_input.rows; x++)
{
for (int y = 0; y < img_input.cols; y++)
{
double temp = saturate_cast<uchar>(img_input.at<uchar>(x, y) + k * generateGaussianNoise(mean, sigma));
img_output.at<uchar>(x, y) = temp;
}
}
}
结果：
文章知识点与官方知识档案匹配，可进
————————————————
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
原文链接：https://blog.csdn.net/qq_44785013/article/details/123853483
扫一扫
369
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
