跳转至

前言

之前做的一些滤波,例如高斯,双边,均值,导向滤波均是在时域下做的滤波。而同态滤波是在频域下来做滤波,用于改善动态范围很大但是暗区细节又不清楚的图像。

原理

对于一般的图像,有这样一个模型f(x,y)=i(x,y)r(x,y)。其中f(x, y)代表原图像,i(x,y)代表照射强度,r(x,y)代表反射强度。一般的图像光照是均匀变化的,所以i应该是低频分量,而不同的物体对光的反射是具有突变型的,所以r(x,y)是高频分量。现在同时对上式两边同时取对数,在做Fourier变化,得到线性组合的频率域:

ln f(x,y)=ln i(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)/D_0)^{2n})]+Rl,巴特沃斯同态滤波器H(u,v)_2=(Rh-Rl)[1/(1+[D_0/cD(u,v)]^{2n})]+RI,指数型同态滤波器:H(u,v)_3=(Rh-Rl)exp(-c[D_0/D(u,v)])^n+RI。其中D(u,v)=[(u-u_0)^2+(v-v_0)^2]^{\frac{1}{2}},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