1. 频域中图像平滑流程
在一副图像中,平滑的图像信号在频谱中贡献低频分量,而高频部分往往是图像的细节和边界,噪声也具有丰富的高频分量。所以对图像进行低通滤波可以消减噪声,但也可能使边界和图像细节变得模糊。
频域中图像滤波的处理流层如下图:
其中DFT为离散傅里叶变换,IDFT为离散傅里叶反变换,
H
(
m
,
n
)
H(m,n)
H(m,n)为滤波器的传递函数。滤波器的输出为:
2. 理想低通滤波器
二维理想低通滤波器(Ideal Low Pass Filter,ILPF )的传递函数定义为:(连续形式给出)
其中
D
(
u
,
v
)
=
(
u
2
+
v
2
)
1
2
,
D
0
D(u,v)=(u^2+v^2)^\frac{1}{2},D0
D(u,v)=(u2+v2)21,D0为截止频率。
H
(
u
,
v
)
H(u,v)
H(u,v)的三维示意图如下:
二维剖面如下:
传递函数
H
(
u
,
v
)
H(u,v)
H(u,v)表示以
D
0
D0
D0为半径的圆域内的分量无损的通过,其他分量则完全的被滤除掉。
对于数字图像,
D
(
u
,
v
)
D(u,v)
D(u,v)可以用离散的
D
(
m
,
n
)
D(m,n)
D(m,n)代替,
D
(
m
,
n
)
D(m,n)
D(m,n)同样表示为某频率分量离零频率分量的距离。因为图像进行DFT后零频率分量在频谱图的四个角,为了观察和操作方便我们需要把零频率分量移到频谱中心,操作完后再移回原位。
下图是一副数字图像的低通滤波示意图:(零频率分量已经移到频谱中心,中心最亮的点既是零频率分量)
对其进行低通滤波,相当于圆内的频率通过,圆外滤除。
理想低通滤波器关键代码:
Mat ImageSmoothing::LowFilter(Mat input, double D0)//理想低通滤波
{//D0为低通滤波器的截止频率
if (input.channels() != 1) {
std::cout << "LowFilter只能处理单通道图像,要处理多通道图像请用:LowFilter3D" << std::endl;
exit(0);
}
Mat src = input.clone();
Mat srcDft = DFT::dft2(src);//傅里叶变换
srcDft = DFT::shift(srcDft);//将零频率分量移到频谱中心
int M = srcDft.rows;
int N = srcDft.cols;
int m1 = M / 2;
int n1 = N / 2;
int channels = srcDft.channels();
for (int m = 0; m < M; m++) {
double* current = srcDft.ptr<double>(m);
for (int n = 0; n < N; n++) {
if (sqrt((m - m1) * (m - m1) + (n - n1) * (n - n1)) > D0) {
current[n*channels] = 0.0;
current[n * channels+1] = 0.0;
}
}
}
srcDft = DFT::shift(srcDft);//还原
Mat dst;
idft(srcDft, dst, DFT_SCALE | DFT_REAL_OUTPUT);
dst.convertTo(dst, CV_8U);
return dst;
}
其中 DFT::shift与 DFT::dft2定义如下:
Mat DFT::shift(Mat input)
{
Mat src = input.clone();
src = src(Range(0, src.rows & -2), Range(0, src.cols & -2));
int y0 = src.rows / 2;
int x0 = src.cols / 2;
Mat q0(src, Rect(0, 0, x0, y0)); //左上角
Mat q1(src, Rect(x0, 0, x0, y0)); //右上角
Mat q2(src, Rect(0, y0, x0, y0)); //左下角
Mat q3(src, Rect(x0, y0, x0, y0)); //右下角
swapArea(q0, q3);
swapArea(q1, q2);
return src;
}
Mat DFT::dft2(Mat input)
{
if (input.channels() != 1) {
std::cout << "只处理单通道的图像" << std::endl;
exit(0);
}
Mat src = input.clone();
//扩充图像尺寸为DFT的最佳尺寸
Mat padded;
int newRows = getOptimalDFTSize(src.rows);
int newCols = getOptimalDFTSize(src.cols);
copyMakeBorder(src, padded, newRows - src.rows, 0, newCols - src.cols, 0, BORDER_CONSTANT, Scalar(0));
//将planes融合合并成一个多通道数组
Mat plans[] = { Mat_<double>(padded),Mat::zeros(padded.size(),CV_64F) };
Mat mergeArray;
merge(plans, 2, mergeArray);
//傅里叶变换
dft(mergeArray, mergeArray);
return mergeArray;
}
上面的代码只能实现单通道图像,对于多通道图像,可以使用上面的LowFilter函数对其每一个通道分别处理即可,代码如下:
Mat ImageSmoothing::LowFilter3D(Mat input, double D0)
{
Mat src = input.clone();
int channels = src.channels();
Mat *dst=new Mat[channels];
split(src, dst);
for (int i = 0; i < channels; i++) {
dst[i] = LowFilter(dst[i],D0);
}
Mat dstImg;
merge(dst, channels, dstImg);
return dstImg;
}
理想低通滤波效果:(取
D
0
=
100
D0=100
D0=100)
理想低通滤波器由于是锐截止,处理后的图像会出现不应有的亮环,也就是“振铃效应”,而且图像也会变得模糊一些。
3. 巴特沃斯低通滤波器
巴特沃斯低通滤波器(Butterworth Low Pass Filter,BLPF)的传递函数如下:
其中
D
(
u
,
v
)
=
(
u
2
+
v
2
)
1
2
D(u,v)=(u^2+v^2)^\frac{1}{2}
D(u,v)=(u2+v2)21,
D
0
D_0
D0为截止频率,
n
n
n为阶数,取正整数,控制衰减速度。
H
(
u
,
v
)
H(u,v)
H(u,v)的三维示意图与二维剖面图如下:
这种滤波器不是锐截止,无“振铃”现象,提高了图像的细节清晰度。
巴特沃斯低通滤波器关键代码:
Mat ImageSmoothing::ButterworthLPF(Mat input, double D0, int n0)
{
//D0为巴特沃斯滤波器的截止频率
//n0控制衰减速度
if (input.channels() != 1) {
std::cout << "ButterworthLPF只能处理单通道图像,要处理多通道图像请用:ButterworthLPF3D" << std::endl;
exit(0);
}
Mat src = input.clone();
Mat srcDft = DFT::dft2(src);
srcDft = DFT::shift(srcDft);//将零频率分量移到频谱中心
int M = srcDft.rows;
int N = srcDft.cols;
int m1 = M / 2;
int n1 = N / 2;
int channels = srcDft.channels();
for (int m = 0; m < M; m++) {
double* current = srcDft.ptr<double>(m);
for (int n = 0; n < N; n++) {
double D= sqrt((m - m1) * (m - m1) + (n - n1) * (n - n1));
double H = 1 / (1 + (sqrt(2) - 1) *pow(D / D0,2*n0));
current[n * channels] *= H;
current[n * channels + 1] *= H;
}
}
srcDft = DFT::shift(srcDft);//还原
Mat dst;
idft(srcDft, dst, DFT_SCALE | DFT_REAL_OUTPUT);
dst.convertTo(dst, CV_8U);
return dst;
}
同样的,对于多通道图像,可以使用上面的ButterworthLPF函数对其每一个通道分别处理即可,代码如下:
Mat ImageSmoothing::ButterworthLPF3D(Mat input, double D0, int n0)
{
Mat src = input.clone();
int channels = src.channels();
Mat* dst = new Mat[channels];
split(src, dst);
for (int i = 0; i < channels; i++) {
dst[i] = ButterworthLPF(dst[i], D0,n0);
}
Mat dstImg;
merge(dst, channels, dstImg);
return dstImg;
}
巴特沃斯低通滤波效果:(取
D
0
=
100
,
n
=
2
D0=100,n=2
D0=100,n=2)
4. 指数低通滤波器
指数滤波器的传递函数为:
其中
D
(
u
,
v
)
=
(
u
2
+
v
2
)
1
2
D(u,v)=(u^2+v^2)^\frac{1}{2}
D(u,v)=(u2+v2)21,
D
0
D_0
D0为截止频率,
n
n
n为阶数,取正整数,控制衰减速度。
H
(
u
,
v
)
H(u,v)
H(u,v)的三维示意图与二维剖面图如下:
指数低通滤波器关键代码:
Mat ImageSmoothing::expLPF(Mat input, double D0, int n0)
{//指数低通滤波器
if (input.channels() != 1) {
std::cout << "expLPF只能处理单通道图像,要处理多通道图像请用:expLPF3D" << std::endl;
exit(0);
}
Mat src = input.clone();
Mat srcDft = DFT::dft2(src);
srcDft = DFT::shift(srcDft);//将零频率分量移到频谱中心
int M = srcDft.rows;
int N = srcDft.cols;
int m1 = M / 2;
int n1 = N / 2;
int channels = srcDft.channels();
for (int m = 0; m < M; m++) {
double* current = srcDft.ptr<double>(m);
for (int n = 0; n < N; n++) {
double D = sqrt((m - m1) * (m - m1) + (n - n1) * (n - n1));
double H = exp((log(1/sqrt(2)))*pow(-D/D0,n0));
current[n * channels] *= H;
current[n * channels + 1] *= H;
}
}
srcDft = DFT::shift(srcDft);//还原
Mat dst;
idft(srcDft, dst, DFT_SCALE | DFT_REAL_OUTPUT);
dst.convertTo(dst, CV_8U);
return dst;
}
同样的,对于多通道图像,可以使用上面的expLPF函数对其每一个通道分别处理即可,代码如下:
Mat ImageSmoothing::expLPF3D(Mat input, double D0, int n0)
{
Mat src = input.clone();
int channels = src.channels();
Mat* dst = new Mat[channels];
split(src, dst);
for (int i = 0; i < channels; i++) {
dst[i] = expLPF(dst[i], D0, n0);
}
Mat dstImg;
merge(dst, channels, dstImg);
return dstImg;
}
指数低通滤波效果:(取
D
0
=
100
,
n
=
2
D0=100,n=2
D0=100,n=2)
5. 梯形低通滤波器
梯形滤波器的传递函数为:
其中
D
(
u
,
v
)
=
(
u
2
+
v
2
)
1
2
D(u,v)=(u^2+v^2)^\frac{1}{2}
D(u,v)=(u2+v2)21,
D
0
D_0
D0为截止频率,
D
0
D_0
D0和
D
1
D_1
D1按图像给定的特点预先给定,且
D
0
<
D
1
D_0<D_1
D0<D1。
H
(
u
,
v
)
H(u,v)
H(u,v)的三维示意图与二维剖面图如下:
梯形低通滤波器关键代码:
Mat ImageSmoothing::trapezoidLPF(Mat input, double D0, double D1)
{//梯形低通滤波
if (input.channels() != 1) {
std::cout << "trapezoidLPF只能处理单通道图像,要处理多通道图像请用:trapezoidLPF3D" << std::endl;
exit(0);
}
if (D0>=D1) {
std::cout << "输入参数有误,参数必须是D0<D1" << std::endl;
exit(0);
}
Mat src = input.clone();
Mat srcDft = DFT::dft2(src);
srcDft = DFT::shift(srcDft);//将零频率分量移到频谱中心
int M = srcDft.rows;
int N = srcDft.cols;
int m1 = M / 2;
int n1 = N / 2;
int channels = srcDft.channels();
for (int m = 0; m < M; m++) {
double* current = srcDft.ptr<double>(m);
for (int n = 0; n < N; n++) {
double D = sqrt((m - m1) * (m - m1) + (n - n1) * (n - n1));
double H;
if (D < D0) {
H = 1;
}
else if (D > D1) {
H = 0;
}
else {
H = (D - D1) / (D0 - D1);
}
current[n * channels] *= H;
current[n * channels + 1] *= H;
}
}
srcDft = DFT::shift(srcDft);//还原
Mat dst;
idft(srcDft, dst, DFT_SCALE | DFT_REAL_OUTPUT);
dst.convertTo(dst, CV_8U);
return dst;
}
同样的,对于多通道图像,可以使用trapezoidLPF函数对其每一个通道分别处理即可,代码如下:
Mat ImageSmoothing::trapezoidLPF3D(Mat input, double D0, double D1)
{
Mat src = input.clone();
int channels = src.channels();
Mat* dst = new Mat[channels];
split(src, dst);
for (int i = 0; i < channels; i++) {
dst[i] = trapezoidLPF(dst[i], D0, D1);
}
Mat dstImg;
merge(dst, channels, dstImg);
return dstImg;
}
梯形低通滤波效果:(取
D
0
=
100
,
D
1
=
150
D_0=100,D_1=150
D0=100,D1=150)