opencv 同态濾波實作 homofilter

同态濾波的原理很簡單，詳細介紹請自行search，本文主要專注于opencv下homofilter的實作，并給出代碼和運作結果，此前的參考了一些網上opencv的homofilter實作代碼，然并卵啊，花了點時間自己做。

算法基本流程如下：

1. 原始輸入圖像: 通常是亮度分量，Y分量，色度分量不用處理

2. ln：對數變化，原理上是說将乘性信号分解為加信信号，具體原理自己搜一下

3. 變換到頻域：dft，fft，dct都行，我用的是dct，dct變換後的結果是實數，節省空間，低頻分量集中在左上角，高頻分量集中在右下角

4. high pass filter：關鍵部分，homofilter觀點認為光照不均勻的部分集中在低頻分量，将低頻分量濾除可以消去光照不均的影響，進而提升暗處和高亮處物體的可見度，但整體亮度變化區間會變小。具體濾波器的設計不多講了，這裡值得注意的是dct結果的原點處的值是原始圖像的能量總和，是不能直接濾除的，否則你得到的是一張黑色的圖檔，什麼也看不見，因為總能量被減小後，圖像的整理亮度下降。是以濾波器原點值設為1，這樣保證濾波後的圖檔和原始圖檔能量一緻，既平均亮度一緻，如果你想提升整體亮度，也可以将濾波原點值設為大于1。

5. 頻域變換到空域：idft, idct...

6. exp：指數運算，ln的反變換

其中3，4，5隻是完成高通濾波，也可以直接在空域用滑動視窗濾波實作，嘗試了7x7的gaussianblur，建構一個高通濾波，效果不好

代碼如下：

void my_HomoFilter(Mat srcImg, Mat &dst)
{
    
    srcImg.convertTo(srcImg, CV_64FC1);
    dst.convertTo(dst, CV_64FC1);
    //1. ln
    for (int i = 0; i < srcImg.rows; i++)
    {
        double* srcdata = srcImg.ptr<double>(i);
        double* logdata = srcImg.ptr<double>(i);
        for (int j = 0; j < srcImg.cols; j++)
        {
            logdata[j] = log(srcdata[j]+0.0001);
        }
    }
    
    //spectrum
    //2. dct
    Mat mat_dct = Mat::zeros(srcImg.rows, srcImg.cols, CV_64FC1);
    dct(srcImg, mat_dct);
    imshow("dct", mat_dct);
    
    //3. linear filter
    Mat H_u_v;
    double gammaH = 1.5;
    double gammaL = 0.5;
    double C = 1;
    double d0 = (srcImg.rows/2)*(srcImg.rows/2) + (srcImg.cols/2)*(srcImg.cols/2);
    double d2 = 0;
    H_u_v = Mat::zeros(srcImg.rows, srcImg.cols, CV_64FC1);
    
    double totalWeight = 0.0;
    for (int i = 0; i < srcImg.rows; i++)
    {
        double * dataH_u_v = H_u_v.ptr<double>(i);
        for (int j = 0; j < srcImg.cols; j++)
        {
            d2 = pow((i), 2.0) + pow((j), 2.0);
            dataH_u_v[j] = 	(gammaH - gammaL)*(1 - exp(-C*d2/d0)) + gammaL;
            totalWeight += dataH_u_v[j];
        }
    }
    H_u_v.ptr<double>(0)[0] = 1.1;
    
    //H_u_v = Mat::ones(srcImg.rows, srcImg.cols, CV_64FC1);
    imshow("H_u_v", H_u_v);


    //imshow("before filter", mat_dct);

    mat_dct = mat_dct.mul(H_u_v);
    //Mat tmp = mat_dct.mul(H_u_v);
    //tmp.copyTo(mat_dct);
    //4. idct
    idct(mat_dct, dst);
    
#if 0
    //spatial high high pass filter
    Mat tmp = Mat::zeros(srcImg.rows, srcImg.cols, CV_64FC1);
    GaussianBlur(srcImg, tmp, Size(9, 9), 1.5, 1.5);
    const double alpha = 0.5;
    
    for (int i = 0; i < srcImg.rows; i++)
    {
        double* srcdata = srcImg.ptr<double>(i);
        double* blurdata = tmp.ptr<double>(i);
        double* dstdata = dst.ptr<double>(i);
        for (int j = 0; j < srcImg.cols; j++)
        {
            dstdata[j] = (1+alpha)*srcdata[j] - alpha*blurdata[j];
            //dstdata[j] = blurdata[j];
            
        }
    }
     
    
    
#endif
    //5. exp
    for (int i = 0; i < srcImg.rows; i++)
    {
        double* srcdata = dst.ptr<double>(i);
        double* dstdata = dst.ptr<double>(i);
        for (int j = 0; j < srcImg.cols; j++)
        {
            dstdata[j] = exp(srcdata[j]);
        }
    }
    
    //imshow("dst", dst);
    dst.convertTo(dst, CV_8UC1);

}
           

主程式代碼：

Mat src_mat = imread("/Users/lilingyu1/Documents/video_enhance/data/dark_1136x1136.jpg", 0);
    imshow("src:", src_mat);
    Mat dst_mat(src_mat.rows, src_mat.cols, src_mat.type());
    my_HomoFilter(src_mat, dst_mat);
   
    
    imshow("dst:", dst_mat);
    cvWaitKey(0);
           

運作結果如下，可以看出整體亮度範圍變小，暗處的細節顯現，缺點是涉及頻域變換，計算複雜，經過濾波不可避免會丢失細節