前言

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

原理

\quad 对于一般的图像,有这样一个模型 f ( x , y ) = i ( x , y ) r ( x , y ) f(x,y)=i(x,y)r(x,y) f(x,y)=i(x,y)r(x,y)。其中f(x, y)代表原图像,i(x,y)代表照射强度,r(x,y)代表反射强度。一般的图像光照是均匀变化的,所以i应该是低频分量,而不同的物体对光的反射是具有突变型的,所以r(x,y)是高频分量。现在同时对上式两边同时取对数,在做Fourier变化,得到线性组合的频率域: l n f ( x , y ) = l n i ( x , y ) + l n r ( x , y ) ln f(x,y)=ln i(x,y)+lnr(x,y) lnf(x,y)=lni(x,y)+lnr(x,y) F F T ( l n f ( x , y ) ) = F F T ( l n i ( x , y ) ) + F F T ( l n r ( x , y ) ) FFT(lnf(x,y))=FFT(lni(x,y))+FFT(lnr(x,y)) FFT(lnf(x,y))=FFT(lni(x,y))+FFT(lnr(x,y)),我们希望对低频能量进行压制,这样就降低了动态范围,而要对高频分量进行提高,这样会增强图像的对比度。因此,同态滤波器的传递函数一般在低频部分小于1,高频部分大于1。

算法过程

  • 对原图像f(x,y)取对数,目的是使得图像模型中的乘法运算转化为简单的加法运算: z ( x , y ) = l n f ( x , y ) = l n i ( x , y ) + l n r ( x , y ) z(x,y)=lnf(x,y)=lni(x,y)+lnr(x,y) z(x,y)=lnf(x,y)=lni(x,y)+lnr(x,y)
  • 再对对数函数做傅里叶变换,将图像转换到频域
  • F ( z ( x , y ) ) = F [ l n i ( x , y ) ] + F [ l n r ( x , y ) ] F(z(x,y))=F[lni(x,y)]+F[lnr(x,y)] F(z(x,y))=F[lni(x,y)]+F[lnr(x,y)] Z = I + R Z=I+R Z=I+R
  • 选择适当的传递函数,压缩照射分量 i ( x , y ) i(x,y) i(x,y)的变化范围,削弱 I ( u , v ) I(u,v) I(u,v),增强反射分量 r ( x , y ) r(x,y) r(x,y)的对比度,所以总结起来就是这个滤波器需要对低频能量进行压制,以降低动态范围,同时要对高频进行提高,以增强图像对比度。所以同态滤波器的传递函数如下:
  • 假设用一个同态滤波器函数 H ( u , v ) H(u,v) H(u,v)来处理原图像f(x,y)的对数的傅里叶变换 Z ( u , v ) Z(u,v) 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(u,v)=H(u,v)Z(u,v)=H(u,v)I(u,v)+H(u,v)R(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 ) ) s(x,y)=F^{-1}(S(u,v)) s(x,y)=F1(S(u,v)),所以 f ( x , y ) = e x p ( s ( x , y ) ) f'(x,y)=exp(s(x,y)) f(x,y)=exp(s(x,y)) f ( x , y ) f'(x,y) f(x,y)代表滤波后的图像。整个算法的过程可以用下面的流程图来表示:

C++ opencv实现

传统的同态滤波器有高斯同态滤波 H ( u , v ) 1 = ( R h R l ) [ 1 e x p ( c ( D ( u , v ) / D 0 ) 2 n ) ] + R l H(u,v)_1=(Rh-Rl)[1-exp(-c(D(u,v)/D_0)^{2n})]+Rl H(u,v)1=(RhRl)[1exp(c(D(u,v)/D0)2n)]+Rl,巴特沃斯同态滤波器 H ( u , v ) 2 = ( R h R l ) [ 1 / ( 1 + [ D 0 / c D ( u , v ) ] 2 n ) ] + R I H(u,v)_2=(Rh-Rl)[1/(1+[D_0/cD(u,v)]^{2n})]+RI H(u,v)2=(RhRl)[1/(1+[D0/cD(u,v)]2n)]+RI,指数型同态滤波器: H ( u , v ) 3 = ( R h R l ) e x p ( c [ D 0 / D ( u , v ) ] ) n + R I H(u,v)_3=(Rh-Rl)exp(-c[D_0/D(u,v)])^n+RI H(u,v)3=(RhRl)exp(c[D0/D(u,v)])n+RI。其中 D ( u , v ) = [ ( u u 0 ) 2 + ( v v 0 ) 2 ] 1 2 D(u,v)=[(u-u_0)^2+(v-v_0)^2]^{\frac{1}{2}} D(u,v)=[(uu0)2+(vv0)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