前言
之前做的一些滤波,例如高斯,双边,均值,导向滤波均是在时域下做的滤波。而同态滤波是在频域下来做滤波,用于改善动态范围很大但是暗区细节又不清楚的图像。
原理
对于一般的图像,有这样一个模型 f(x,y)=i(x,y)r(x,y)。其中f(x, y)代表原图像,i(x,y)代表照射强度,r(x,y)代表反射强度。一般的图像光照是均匀变化的,所以i应该是低频分量,而不同的物体对光的反射是具有突变型的,所以r(x,y)是高频分量。现在同时对上式两边同时取对数,在做Fourier变化,得到线性组合的频率域: lnf(x,y)=lni(x,y)+lnr(x,y) FFT(lnf(x,y))=FFT(lni(x,y))+FFT(lnr(x,y)),我们希望对低频能量进行压制,这样就降低了动态范围,而要对高频分量进行提高,这样会增强图像的对比度。因此,同态滤波器的传递函数一般在低频部分小于1,高频部分大于1。
算法过程
- 对原图像f(x,y)取对数,目的是使得图像模型中的乘法运算转化为简单的加法运算: z(x,y)=lnf(x,y)=lni(x,y)+lnr(x,y)
- 再对对数函数做傅里叶变换,将图像转换到频域
- F(z(x,y))=F[lni(x,y)]+F[lnr(x,y)]即 Z=I+R
- 选择适当的传递函数,压缩照射分量 i(x,y)的变化范围,削弱 I(u,v),增强反射分量 r(x,y)的对比度,所以总结起来就是这个滤波器需要对低频能量进行压制,以降低动态范围,同时要对高频进行提高,以增强图像对比度。所以同态滤波器的传递函数如下:
- 假设用一个同态滤波器函数 H(u,v)来处理原图像f(x,y)的对数的傅里叶变换 Z(u,v)得到: S(u,v)=H(u,v)Z(u,v)=H(u,v)I(u,v)+H(u,v)R(u,v),逆变换到时域得: s(x,y)=F−1(S(u,v)),所以 f′(x,y)=exp(s(x,y)), f′(x,y)代表滤波后的图像。整个算法的过程可以用下面的流程图来表示:
C++ opencv实现
传统的同态滤波器有高斯同态滤波 H(u,v)1=(Rh−Rl)[1−exp(−c(D(u,v)/D0)2n)]+Rl,巴特沃斯同态滤波器 H(u,v)2=(Rh−Rl)[1/(1+[D0/cD(u,v)]2n)]+RI,指数型同态滤波器: H(u,v)3=(Rh−Rl)exp(−c[D0/D(u,v)])n+RI。其中 D(u,v)=[(u−u0)2+(v−v0)2]21,Rh,RI分别为高频和低频增益,D0为截止频率,c是控制斜面锐化的常数。这里实现了第一种滤波器,也就是高斯滤波。
Mat HomoFilter(cv::Mat src){
src.convertTo(src, CV_64FC1);
int rows = src.rows;
int cols = src.cols;
int m = rows % 2 == 1 ? rows + 1 : rows;
int n = cols % 2 == 1 ? cols + 1 : cols;
copyMakeBorder(src, src, 0, m - rows, 0, n - cols, BORDER_CONSTANT, Scalar::all(0));
rows = src.rows;
cols = src.cols;
Mat dst(rows, cols, CV_64FC1);
//1. ln
for(int i = 0; i < rows; i++){
double *srcdata = src.ptr<double>(i);
double *logdata = src.ptr<double>(i);
for(int j = 0; j < cols; j++){
logdata[j] = log(srcdata[j] + 0.0001);
}
}
//2. dct
Mat mat_dct = Mat::zeros(rows, cols, CV_64FC1);
dct(src, mat_dct);
//3. 高斯同态滤波器
Mat H_u_v;
double gammaH = 1.5;
double gammaL = 0.5;
double C = 1;
double d0 = (src.rows / 2) * (src.rows / 2) + (src.cols / 2) * (src.cols / 2);
double d2 = 0;
H_u_v = Mat::zeros(rows, cols, CV_64FC1);
for(int i = 0; i < rows; i++){
double * dataH_u_v = H_u_v.ptr<double>(i);
for(int j = 0; j < cols; j++){
d2 = pow(i, 2.0) + pow(j, 2.0);
dataH_u_v[j] = (gammaH - gammaL) * (1 - exp(-C * d2 / d0)) + gammaL;
}
}
H_u_v.ptr<double>(0)[0] = 1.1;
mat_dct = mat_dct.mul(H_u_v);
//4. idct
idct(mat_dct, dst);
//exp
for(int i = 0; i < rows; i++){
double *srcdata = dst.ptr<double>(i);
double *dstdata = dst.ptr<double>(i);
for(int j = 0; j < cols; j++){
dstdata[j] = exp(srcdata[j]);
}
}
dst.convertTo(dst, CV_8UC1);
return dst;
}
主函数:
int main(){
//
Mat src = cv::imread("./1.png");
imshow("origin", src);
int originrows = src.rows;
int origincols = src.cols;
Mat dst(src.rows, src.cols, CV_8UC3);
cvtColor(src, src, COLOR_BGR2YUV);
vector <Mat> yuv;
split(src, yuv);
Mat nowY = yuv[0];
Mat newY = HomoFilter(nowY);
Mat tempY(originrows, origincols, CV_8UC1);
for(int i = 0; i < originrows; i++){
for(int j = 0; j < origincols; j++){
tempY.at<uchar>(i, j) = newY.at<uchar>(i, j);
}
}
yuv[0] = tempY;
merge(yuv, dst);
cvtColor(dst, dst, COLOR_YUV2BGR);
imshow("result", dst);
waitKey(0);
}
效果:
可以看到图片的整体亮度被拉到人眼可以接受的范围里面了,但是缺点是有些地方出现了白点。再看一组效果:
我的这个代码效果很不好,如果真的要使用同态滤波需要认真取做研究,可以参考https://github.com/lilingyu/homofilter