图像噪声
噪声图像是指在图像中存在随机干扰或不规则变化的图像。以下是几种常见的噪声图像类型:
-
高斯噪声图像:高斯噪声是指在图像中每个像素值上添加的服从高斯分布的随机噪声。它是最常见的图像噪声类型之一。
//对原始图像添加高斯噪声 Mat addGaussianNoise(Mat image, double mean, double stddev) { //image 原始图像 mean stddev // 创建随机数生成器 default_random_engine generator; normal_distribution<double> distribution(mean, stddev); // 添加高斯噪声 Mat noise_image = image.clone(); for (int i = 0; i < noise_image.rows; i++) { for (int j = 0; j < noise_image.cols; j++) { double noise = distribution(generator); int pixel = noise_image.at<uchar>(i, j) + noise; if (pixel < 0) { pixel = 0; } else if (pixel > 255) { pixel = 255; } noise_image.at<uchar>(i, j) = pixel; } } return noise_image; }
-
椒盐噪声图像:椒盐噪声是指在图像中随机出现亮白或暗黑的像素点,使得这些像素点的值变得与周围像素不一致。它模拟了图像传感器或信号传输中的随机错误。
//对原始图像添加椒盐噪声 void addSaltAndPepperNoise(cv::Mat &src, cv::Mat&dst, double noise_ratio, int border_size) { //src原始图像,dst结果图像,noise_ratio噪声比例,border_size边界距离 dst = src.clone(); int rows = src.rows - 2 * border_size; int cols = src.cols - 2 * border_size; int num_noise_pixels = static_cast<int>(noise_ratio * rows * cols); for (int i = 0; i < num_noise_pixels; ++i) { //随机点位 int row = rand() % rows + border_size; int col = rand() % cols + border_size; if (rand() % 2 == 0) { dst.at<uchar>(row, col) = 0; // 椒噪声 } else { dst.at<uchar>(row, col) = 255; // 盐噪声 } } }
-
瑞利噪声图像:瑞利噪声是指在图像中存在强度不均匀的干扰,常见于无线通信或医学图像中。它的幅度遵循瑞利分布。
-
指数噪声图像:指数噪声是指在图像中存在幅度随机变化的干扰,它的幅度遵循指数分布。
-
拉普拉斯噪声图像:拉普拉斯噪声是指在图像中存在有尖峰和宽峰的干扰,它的幅度遵循拉普拉斯分布。
-
泊松噪声图像:泊松噪声是指在图像中存在随机出现的光子计数的波动干扰,它的幅度遵循泊松分布。泊松噪声在低光条件下的图像中常见。
这些是常见的噪声类型,它们可以单独存在,也可以以组合形式出现在图像中。此外,还有其他一些特定于特定应用领域的噪声类型,例如相机传感器噪声、压缩噪声等。
图像平滑
图像平滑是指在图像中减少噪声或过滤掉不必要的细节,使图像更加平滑和清晰的过程。以下是几种常用的图像平滑方法:
-
均值滤波(Mean Filter):将每个像素的值替换为其邻域像素的平均值。这种方法简单有效,但可能会导致图像细节的模糊。
-
中值滤波(Median Filter):将每个像素的值替换为其邻域像素值的中值。中值滤波能够有效地去除椒盐噪声等离群点,对于保留边缘细节效果较好。
-
高斯滤波(Gaussian Filter):通过将每个像素与其邻域像素进行加权平均,使用高斯函数作为权重。高斯滤波器对于平滑图像并保持边缘信息非常有效。
-
双边滤波(Bilateral Filter):在计算加权平均值时考虑像素的空间距离和像素值之间的相似性。这种滤波方法可以在保持边缘细节的同时进行噪声抑制。
-
邻近平均滤波(Neighborhood Averaging Filter):将每个像素的值替换为其邻域像素的加权平均值,其中较远的像素具有较小的权重。该方法可用于平滑图像并减少噪声。
这些方法可以单独使用,也可以组合使用,根据具体的图像和应用需求选择适合的方法。此外,还有其他一些图像平滑方法,如双向滤波、总变差滤波等,可以根据具体情况进行选择和尝试。
中值滤波器
a)原始中值滤波器
01 原理
中值滤波是基于排序统计理论的一种能有效抑制噪声的非线性信号处理技术,中值滤波的基本原理是把数字图像或数字序列中一点的值用该点的一个邻域中各点值的中值代替,让周围的像素值接近的真实值,从而消除孤立的噪声点。方法是用某种结构的二维滑动模板,将板内像素按照像素值的大小进行排序,生成单调上升(或下降)的为二维数据序列。
02代码实现
//中值滤波器
void MedianFlitering(Mat src, Mat & dst, int ksize)
{
int R = ksize / 2;//窗口半径
dst = src.clone();
//对图像边界进行处理
cv::Mat paddedSrc;
copyMakeBorder(src, paddedSrc, R, R, R, R, BORDER_REFLECT);//对图像边界添加宽度为R个像素的镜像边框
for (int y = R; y < dst.rows - R; y++)
{
for (int x = R; x < dst.cols - R; x++)
{
vector<uchar> values;
for (int j = -R; j <= R; j++) {
for (int i = -R; i <= R; i++) {
values.push_back(paddedSrc.at<uchar>(y + j, x + i));
}
}
sort(values.begin(), values.end());//对像素值进行排序
dst.at<uchar>(y, x) = values[values.size() / 2];
}
}
}
b)直方图加速中值滤波器
01原理
直方图加速中值滤波器的基本原理是利用直方图统计和查找表的思想,通过预先计算滤波窗口内像素值的直方图,并根据直方图信息快速找到中值。
02代码实现
//用直方图对中值滤波器加速
void FastMedianFilter(const Mat& src, Mat& dst, int ksize) {
int R = ksize / 2;
int window_size = ksize * ksize;
int median_position = window_size / 2;//中间值位置
dst = src.clone();
//对图像边界进行处理
Mat paddedSrc;
copyMakeBorder(src, paddedSrc, R, R, R, R, BORDER_REPLICATE);//BORDER_REPLICATE复制边界填充
for (int y = R; y < dst.rows - R; y++) {
for (int x = R; x < dst.cols - R; x++) {
int histogram[256] = { 0 };//设置直方图
for (int j = -R; j <= R; j++) {
for (int i = -R; i <= R; i++) {
histogram[src.at<uchar>(y + j, x + i)]++;//统计在该灰度值点的个数
}
}
int count = 0;
for (int i = 0; i < 256; i++) {
count += histogram[i];//对窗口内点进行计数
if (count > median_position) {
//对边缘噪点进行处理
dst.at<uchar>(y - R, x - R) = i;//忽略填充边界
break;
}
}
}
}
}
c)实验
(1)不同窗口下三种实现方式速度对比(对椒盐噪声)
3*3(second) | 5*5(second) | 7*7(second) | |
---|---|---|---|
opencv中值滤波器 | 0.0002647 | 0.0013674 | 0.0060189 |
中值滤波 | 0.0646833 | 0.124498 | 0.179576 |
直方图加速中值滤波器 | 0.0200518 | 0.0192534 | 0.0214629 |
- 窗口越小速度越快
中值滤波器实现中,复杂度主要来自两部分:
-
对于每个像素,遍历一个大小为
ksize*ksize
的窗口并收集值。如果ksize
是K
,那么这部分的复杂度是O(K^2)
。 -
对于收集的每个像素值,对它们进行排序。在一般情况下,排序算法的时间复杂度是
O(n log n)
。在这个例子中,n
是窗口大小,也就是K^2
。
对于每个像素,总的复杂度是 O(K^2 log K^2)
。又因为对图像中的每个像素都执行这个操作,所以,如果图像的大小是 M*N
(M
行,N
列),那么整个中值滤波器的时间复杂度就是 O(M*N*K^2 log K^2)
。
直方图方法,可以将中值滤波器的复杂度降低到 O(M*N*K^2)
。
(2)讨论分析中值滤波在什么样的噪声图像下效果比较好
对比噪声类型:椒盐噪声,高斯噪声
处理椒盐噪声的效果较好
双边滤波器
双边高斯滤波器
01原理
双边滤波的核函数是空间域核与像素范围域核的综合结果:在图像的平坦区域,像素值变化很小,对应的像素范围域权重接近于1,此时空间域权重起主要作用,相当于进行高斯模糊;在图像的边缘区域,像素值变化很大,像素范围域权重变大,从而保持了边缘的信息。
02代码实现
//双边高斯滤波器
void bilateralFilterCustom(const Mat& src, Mat& dst, int d, double sigmaColor, double sigmaSpace) {
int radius = d / 2; // 计算滤波器半径
double colorCoeff = -0.5 / (sigmaColor * sigmaColor); // 颜色权重系数
double spaceCoeff = -0.5 / (sigmaSpace * sigmaSpace); // 空间权重系数
Mat srcPadded; // 定义扩展后的源图像
copyMakeBorder(src, srcPadded, radius, radius, radius, radius, BORDER_REFLECT); // 对源图像进行扩展
dst = Mat::zeros(src.size(), src.type()); // 初始化输出图像
// 遍历扩展后的源图像
for (int y = radius; y < srcPadded.rows - radius; y++) {
for (int x = radius; x < srcPadded.cols - radius; x++) {
double sum = 0.0; // 初始化滤波器加权和
double wSum = 0.0; // 初始化权重和
// 遍历滤波器窗口
for (int j = -radius; j <= radius; j++) {
for (int i = -radius; i <= radius; i++) {
double spaceDist2 = i * i + j * j; // 计算空间距离的平方
double colorDist2 = pow(srcPadded.at<uchar>(y + j, x + i) - srcPadded.at<uchar>(y, x), 2); // 计算颜色差的平方
double w = exp(spaceDist2 * spaceCoeff) * exp(colorDist2 * colorCoeff); // 计算权重
sum += srcPadded.at<uchar>(y + j, x + i) * w; // 累加加权和
wSum += w; // 累加权重和
}
}
dst.at<uchar>(y - radius, x - radius) = round(sum / wSum); // 计算输出像素值
}
}
}
直方图加速的双边滤波
01 原理
02具体实现步骤
根据给定的源图像,使用直方图加速的双边滤波器的自定义实现步骤如下:
-
计算空间权重:创建一个大小为d x d的空间权重矩阵
space_weight
,遍历每个像素位置(i, j)
,计算空间距离((i - d / 2)^2 + (j - d / 2)^2)
的负指数值,并将其赋给对应位置的空间权重值。 -
创建颜色权重查找表:创建一个大小为256 x 1的颜色权重矩阵
color_weight
,遍历每个灰度值i
,计算其平方的负指数值,并将其赋给对应位置的颜色权重值。 -
对源图像进行边界扩展:创建一个扩展后的源图像
padded_src
,使用copyMakeBorder
函数将原图像src
的边界进行镜像填充,填充宽度为d/2
。 -
初始化目标图像:创建一个与源图像大小相同的目标图像
dst
,并将其初始化为全零图像。 -
应用直方图加速双边滤波器:对于每个像素位置
(y, x)
在源图像中,遍历滤波器窗口。对于窗口内的每个位置(dy, dx)
,计算当前像素与中心像素的灰度差color_diff
,获取对应的空间权重和颜色权重。累加每个窗口位置的加权像素值到sum
中,并累加权重到wsum
中。 -
计算输出像素值:通过除以权重和
wsum
,将sum
除以wsum
得到输出像素值,并使用static_cast
将其转换为uchar
类型。将计算得到的像素值赋给目标图像dst
的相应位置。
最终,函数返回输出图像dst
。
03代码实现
//用直方图对双边滤波器进行加速
void fastBilateralFilter(const Mat& src, Mat& dst, int d, double sigma_color, double sigma_space) {
//计算空间权重
Mat space_weight(d, d, CV_32F);
for (int i = 0; i < d; i++) {
for (int j = 0; j < d; j++) {
float g = static_cast<float>(exp(-((i - d / 2) * (i - d / 2) + (j - d / 2) * (j - d / 2)) / (2 * sigma_space * sigma_space)));
space_weight.at<float>(i, j) = g;
}
}
// 创建颜色权重查找表
const int L = 256;
Mat color_weight(L, 1, CV_32F);//生成颜色直方图,颜色差距最大为255
for (int i = 0; i < L; i++) {
float g = static_cast<float>(exp(-(i * i) / (2 * sigma_color * sigma_color)));
color_weight.at<float>(i, 0) = g;
}
// 对源图像进行边界扩展
Mat padded_src;
copyMakeBorder(src, padded_src, d / 2, d / 2, d / 2, d / 2, BORDER_REFLECT);
// 初始化目标图像
dst = Mat::zeros(src.size(), src.type());
// 应用直方图加速双边滤波器
for (int y = 0; y < src.rows; y++) {
for (int x = 0; x < src.cols; x++) {
float sum = 0;
float wsum = 0;
for (int dy = -d / 2; dy <= d / 2; dy++) {
for (int dx = -d / 2; dx <= d / 2; dx++) {
int color_diff = abs(padded_src.at<uchar>(y + dy + d / 2, x + dx + d / 2) - padded_src.at<uchar>(y + d / 2, x + d / 2));//计算颜色差距color_diff,一定会落在颜色直方图内
float weight = space_weight.at<float>(dy + d / 2, dx + d / 2) * color_weight.at<float>(color_diff, 0);
sum += weight * padded_src.at<uchar>(y + dy + d / 2, x + dx + d / 2);
wsum += weight;
}
}
dst.at<uchar>(y, x) = static_cast<uchar>(sum / wsum);
}
}
}
c)问题
-
双边高斯滤波和传统高斯滤波对比
-
三种实现方式速度对比
opencv双边高斯滤波器 | 0.014423 |
双边高斯滤波器 | 1.51627 |
直方图加速双边高斯滤波器 | 0.0336481 |
在不同噪声下进行比较
噪声类型:椒盐噪声、高斯噪声
维纳滤波器
01原理
维纳滤波器,也被称为最小均方误差滤波器,是一种线性滤波器,用于恢复受到噪声影响的图像。其基本的工作原理是通过最小化输出误差的平方均值来估计原始图像。
在图像拍摄过程中由于各种原因会造成图像退化,图像退化模型如下
g
(
x
,
y
)
=
h
(
x
,
y
)
⋆
f
(
x
,
y
)
+
η
(
x
,
y
)
g(x,y)=h(x,y)⋆f(x,y)+η(x,y)
g(x,y)=h(x,y)⋆f(x,y)+η(x,y)
傅里叶变换后
G
(
u
,
v
)
=
H
(
u
,
v
)
F
(
u
,
v
)
+
N
(
u
,
v
)
G(u,v)=H(u,v)F(u,v)+N(u,v)
G(u,v)=H(u,v)F(u,v)+N(u,v)
如果退化函数和加性噪声都考虑,空域滤波器无法解决图像退化问题,逆滤波效果因为噪声的存在会变得非常差,这个时候就需要用到维纳滤波,(维纳滤波的推导写在结论中)维纳滤波公式如下:
F
(
u
,
v
)
=
[
H
(
u
,
v
)
1
∣
H
(
u
,
v
)
∣
2
+
S
η
(
u
,
v
)
/
S
f
(
u
,
v
)
∣
H
(
u
,
v
)
∣
2
]
G
(
u
,
v
)
F^(u,v)=[ H(u,v)1∣H(u,v)∣ 2 +S η(u,v)/S f (u,v)∣H(u,v)∣ 2 ]G(u,v)
F(u,v)=[H(u,v)1∣H(u,v)∣2+Sη(u,v)/Sf(u,v)∣H(u,v)∣2]G(u,v)
其中$ Sη(u,v)=∣N(u,v)∣2$ 为噪声的功率谱,这个我们可以通过用户输入的方差构造一个噪声图像
N
(
u
,
v
)
N(u,v)
N(u,v)并计算功率谱;
S f ( u , v ) = ∣ F ( u , v ) ∣ 2 Sf(u,v)=∣F(u,v)∣2 Sf(u,v)=∣F(u,v)∣2为输入图像的功率谱,
退化函数进行了简化,将退化函数置为1,因此维纳滤波公式简化为:
F
(
u
,
v
)
=
[
S
f
(
u
,
v
)
S
f
(
u
,
v
)
+
S
η
(
u
,
v
)
]
G
(
u
,
v
)
F
^
(
u
,
v
)
=
[
S
f
(
u
,
v
)
S
f
(
u
,
v
)
+
S
η
(
u
,
v
)
]
G
(
u
,
v
)
F
(
u
,
v
)
=
[
S
f
(
u
,
v
)
+
S
η
(
u
,
v
)
S
f
(
u
,
v
)
]
G
(
u
,
v
)
F ^ ( u , v ) = [ S f ( u , v ) S f ( u , v ) + S η ( u , v ) ] G ( u , v ) \hat{F}(u, v)=\left[\frac{S_{f}(u, v)}{S_{f}(u, v)+S_{\eta}(u, v)}\right] G(u, v) F ^ (u,v)=[ S f (u,v)+S η (u,v) S f (u,v) ]G(u,v)
F(u,v)=[Sf(u,v)Sf(u,v)+Sη(u,v)]G(u,v)F^(u,v)=[Sf(u,v)+Sη(u,v)Sf(u,v)]G(u,v)F(u,v)=[Sf(u,v)+Sη(u,v)Sf(u,v)]G(u,v)
维纳滤波器主要作用
- 噪声抑制: 维纳滤波器可以用于抑制图像中的噪声,例如高斯噪声、白噪声等。由于维纳滤波器可以根据噪声的功率谱和信号的功率谱动态调整滤波参数,因此能够更好地保留图像的细节信息。
- 图像去模糊: 当图像受到模糊影响时,例如运动模糊、聚焦模糊等,维纳滤波器可以用于恢复清晰的图像。维纳滤波器的去模糊效果在知道模糊函数(也称为点扩散函数)的情况下尤为有效。
- 图像复原: 在一些应用中,例如医学成像、遥感成像等,我们可能需要从受到噪声和模糊影响的图像中恢复出原始的、清晰的图像。这时,维纳滤波器是一种常用的图像复原工具。
02具体步骤
-
创建几个Mat对象,包括
pad
(经过零填充的原始图像),cpx
(双通道频域图像),dst
(滤波后的图像)。 -
获取最佳的傅里叶变换尺寸,使其为2的指数。
-
使用
copyMakeBorder
函数对原始图像进行零填充,获得填充后的图像pad
。 -
创建一个临时图像
tmpR
,将参考图像ref
调整为与pad
相同的尺寸。 -
调用
GetSpectrum
函数,传入临时图像tmpR
,获取参考图像的频谱refSpectrum
。 -
创建一个临时图像
tmpN
,生成服从标准正态分布的噪声图像。 -
调用
GetSpectrum
函数,传入临时图像tmpN
,获取噪声图像的频谱noiseSpectrum
。 -
创建一个双通道Mat数组
planes
,其中第一个通道为pad
的浮点型副本,第二个通道初始化为与pad
相同尺寸的全零图像。 -
使用
merge
函数将两个通道合并为一个复数频域图像cpx
。 -
使用
dft
函数对复数频域图像cpx
进行离散傅里叶变换。 -
使用
split
函数将复数频域图像cpx
分离为实部和虚部图像,并存储在planes
数组中。 -
计算维纳滤波因子,
-
将参考图像的频谱
refSpectrum
除以参考图像频谱与噪声图像频谱之和,得到factor
。 -
使用
multiply
函数将实部和虚部图像与维纳滤波因子factor
相乘,得到滤波后的频域图像。 -
使用
merge
函数将实部和虚部图像重新合并为一个复数频域图像cpx
。 -
使用
idft
函数对复数频域图像cpx
进行反离散傅里叶变换,并使用DFT_SCALE
和DFT_REAL_OUTPUT
标志进行缩放和提取实部。 -
使用
Rect
函数裁剪滤波后的图像dst
,以去除填充的部分,使其恢复到原始大小。 -
将滤波后的图像
dst
转换为8位无符号整数类型,并将其作为函数的返回值。
03代码实现
/维纳滤波
Mat GetSpectrum(const Mat &src)//求功率谱
{
Mat dst, cpx;//dst功率谱,cpx
Mat planes[] = { Mat_<float>(src), Mat::zeros(src.size(), CV_32F) };//分离为两个通道
merge(planes, 2, cpx);//合并通道
dft(cpx, cpx);
split(cpx, planes);//分离通道
magnitude(planes[0], planes[1], dst);//计算频域幅度值
//频谱就是频域幅度图的平方
multiply(dst, dst, dst);
return dst;
}
Mat WienerFilter(const Mat &src, const Mat &ref, int stddev)
{
Mat pad, cpx, dst;//pad是原图像0填充后的图像,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));//BORDER_CONSTANT添加恒定的像素值
//获得参考图片频谱
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));//生成正态分布随机的矩阵,高斯噪声图像
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);
//维纳滤波因子,这里退化函数设置为1,简化运算
Mat factor = refSpectrum / (refSpectrum + noiseSpectrum);
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);
// 裁剪到原始大小
dst = dst(Rect(0, 0, src.cols, src.rows));
dst.convertTo(dst, CV_8UC1);
return dst;
}
问题
(1)对模糊图像进行复原