濾波原理

簡(jiǎn)而言之，圖像的同態(tài)濾波是基于以入射光和反射光為基礎(chǔ)的圖像模型上的，如果把圖像函數(shù)F(x,y)表示為光照函數(shù)，即照射分量i(x,y)與反射分量r(x,y)兩個(gè)分量的乘積，那么圖像的模型可以表示為F(x,y)= i(x,y)*r(x,y)。通過對(duì)照射分量i(x,y)和反射分量r(x,y)的研究可知，照射分量一般反映灰度的恒定分量，相當(dāng)于頻域中的低頻信息，減弱入射光就可以起到縮小圖像灰度范圍的作用；而反射光與物體的邊界特性是密切相關(guān)的，相當(dāng)于頻域中的高頻信息，增強(qiáng)反射光就可以起到提高圖像對(duì)比度的作用。同態(tài)濾波器屬于頻域?yàn)V波器范疇，過濾低頻信息，放大高頻信息，達(dá)到壓縮圖像的灰度空間，并擴(kuò)展對(duì)比度的目的，其傳遞函數(shù)表現(xiàn)為在低頻部分小于1，高頻部分大于1。具體細(xì)節(jié)參考Homomorphic filtering

傳遞函數(shù)：H(u,v)

原理流程圖：

---

濾波器設(shè)計(jì)

根據(jù)上述濾波器的特性與傳遞函數(shù)，不難選擇濾波器。表達(dá)式如圖所示。

這里參數(shù)c控制上升速率，參考值2，D0參數(shù)我們通過高斯低通濾波器來說明，表示的就是橫軸截止時(shí)的值,如圖所示。

參數(shù)D0的選擇，是比較講究的。對(duì)應(yīng)高通濾波器，其決定過濾圖像能量的大小。如何自適應(yīng)取值，我們后面討論。

OpenCV實(shí)現(xiàn)

數(shù)據(jù)處理
按照上述的原理流程圖，第一步需要對(duì)數(shù)據(jù)作相應(yīng)的處理，然后進(jìn)行l(wèi)n取對(duì)數(shù)操作，這里往往容易出問題。OpenCV讀取的圖像數(shù)據(jù)f(x,y)正常是0-255范圍的，直接f(x,y)/255.0 + 1.0，然后取對(duì)數(shù)即可。此處，如果使用normalize函數(shù)，然后再加1，后面會(huì)發(fā)現(xiàn)最終結(jié)果g(x,y)是一片黑，原因是normalize對(duì)原圖的操作是a*f(x,y) + b, 根據(jù)傅里葉變換的性質(zhì)就會(huì)發(fā)現(xiàn)問題出現(xiàn)在哪里了。

DFT
OpenCV中有實(shí)現(xiàn)dft的算法，并且在Sample/c++中給出了示例，注意其中將頻譜中心移到圖像中心的操作。這里補(bǔ)充一點(diǎn)是計(jì)算圖像的功率譜，還記得上述中提及的D0計(jì)算么？ 功率譜計(jì)算公式如下：

巴特沃斯濾波器
根據(jù)求解的參數(shù)以及濾波器的傳遞函數(shù)，即可求解出濾波器，此時(shí)與DFT的結(jié)果一致，濾波器也包含實(shí)數(shù)與虛數(shù)兩部分

頻域?yàn)V波
調(diào)用mulSpectrums進(jìn)行計(jì)算，得到濾波結(jié)果,此時(shí)注意將結(jié)果的頻域中心移回到圖像的左上角。IDFT
OpenCV中有實(shí)現(xiàn)idft的算法,注意設(shè)置參數(shù)DFT_SCALE，具體查看OpenCV的官方文檔說明，接著對(duì)結(jié)果計(jì)算e指數(shù)，再分解結(jié)果的實(shí)部和虛部，計(jì)算幅值，最后幅值減1.

數(shù)據(jù)恢復(fù)
使用normalize方法，將最大值置為255，并將數(shù)據(jù)格式轉(zhuǎn)換為CV_8U，得到最終結(jié)果。

---

代碼

濾波器：

//High-Frequency-Emphasis Filters
Mat Butterworth_Homomorphic_Filter(Size sz, int D, int n, float high_h_v_TB, float low_h_v_TB,Mat& realIm)
{
    Mat single(sz.height, sz.width, CV_32F);
    cout <<sz.width <<" "<< sz.height<<endl;
    Point centre = Point(sz.height/2, sz.width/2);
    double radius;
    float upper = (high_h_v_TB * 0.1);
    float lower = (low_h_v_TB * 0.1);
    long dpow = D*D;
    float W = (upper - lower);

    for(int i = 0; i < sz.height; i++)
    {
        for(int j = 0; j < sz.width; j++)
        {
            radius = pow((float)(i - centre.x), 2) + pow((float) (j - centre.y), 2);
            float r = exp(-n*radius/dpow);
            if(radius < 0)
                single.at<float>(i,j) = upper;
            else
                single.at<float>(i,j) =W*(1 - r) + lower;       
        }
    }

    single.copyTo(realIm);
    Mat butterworth_complex;
    //make two channels to match complex
    Mat butterworth_channels[] = {Mat_<float>(single), Mat::zeros(sz, CV_32F)};
    merge(butterworth_channels, 2, butterworth_complex);

    return butterworth_complex;
}

---

DFT變換：

//DFT 返回功率譜Power
Mat Fourier_Transform(Mat frame_bw, Mat& image_complex,Mat &image_phase, Mat &image_mag)
{
    Mat frame_log;
    frame_bw.convertTo(frame_log, CV_32F);
    frame_log = frame_log/255;
    /*Take the natural log of the input (compute log(1 + Mag)*/
    frame_log += 1;
    log( frame_log, frame_log); // log(1 + Mag)

    /*2. Expand the image to an optimal size
    The performance of the DFT depends of the image size. It tends to be the fastest for image sizes that are multiple of 2, 3 or 5.
    We can use the copyMakeBorder() function to expand the borders of an image.*/
    Mat padded;
    int M = getOptimalDFTSize(frame_log.rows);
    int N = getOptimalDFTSize(frame_log.cols);
    copyMakeBorder(frame_log, padded, 0, M - frame_log.rows, 0, N - frame_log.cols, BORDER_CONSTANT, Scalar::all(0));

    /*Make place for both the complex and real values
    The result of the DFT is a complex. Then the result is 2 images (Imaginary + Real), and the frequency domains range is much larger than the spatial one. Therefore we need to store in float !
    That's why we will convert our input image "padded" to float and expand it to another channel to hold the complex values.
    Planes is an arrow of 2 matrix (planes[0] = Real part, planes[1] = Imaginary part)*/
    Mat image_planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
    /*Creates one multichannel array out of several single-channel ones.*/
    merge(image_planes, 2, image_complex);

    /*Make the DFT
    The result of thee DFT is a complex image : "image_complex"*/
    dft(image_complex, image_complex);

    /***************************/
    //Create spectrum magnitude//
    /***************************/
    /*Transform the real and complex values to magnitude
    NB: We separe Real part to Imaginary part*/
    split(image_complex, image_planes);
    //Starting with this part we have the real part of the image in planes[0] and the imaginary in planes[1]
    phase(image_planes[0], image_planes[1], image_phase);
    magnitude(image_planes[0], image_planes[1], image_mag);

    //Power
    pow(image_planes[0],2,image_planes[0]);
    pow(image_planes[1],2,image_planes[1]);

    Mat Power = image_planes[0] + image_planes[1];

    return Power;
}

---

IDFT變換

void Inv_Fourier_Transform(Mat input, Mat& inverseTransform)
{
    /*Make the IDFT*/
    Mat result;
    idft(input, result,DFT_SCALE);
    /*Take the exponential*/
    exp(result, result);

    vector<Mat> planes;
    split(result,planes);
    magnitude(planes[0],planes[1],planes[0]);
    planes[0] = planes[0] - 1.0;
    normalize(planes[0],planes[0],0,255,CV_MINMAX);

    planes[0].convertTo(inverseTransform,CV_8U);
}

---

頻譜移位

//SHIFT
void Shifting_DFT(Mat &fImage)
{
    //For visualization purposes we may also rearrange the quadrants of the result, so that the origin (0,0), corresponds to the image center.
    Mat tmp, q0, q1, q2, q3;

    /*First crop the image, if it has an odd number of rows or columns.
    Operator & bit to bit by -2 (two's complement : -2 = 111111111....10) to eliminate the first bit 2^0 (In case of odd number on row or col, we take the even number in below)*/
    fImage = fImage(Rect(0, 0, fImage.cols & -2, fImage.rows & -2));
    int cx = fImage.cols/2;
    int cy = fImage.rows/2;

    /*Rearrange the quadrants of Fourier image so that the origin is at the image center*/
    q0 = fImage(Rect(0, 0, cx, cy));
    q1 = fImage(Rect(cx, 0, cx, cy));
    q2 = fImage(Rect(0, cy, cx, cy));
    q3 = fImage(Rect(cx, cy, cx, cy));

    /*We reverse each quadrant of the frame with its other quadrant diagonally opposite*/
    /*We reverse q0 and q3*/
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);

    /*We reverse q1 and q2*/
    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2);
}

主函數(shù)：

void homomorphicFiltering(Mat src,Mat& dst)
{
    Mat img;
    Mat imgHls;
    vector<Mat> vHls;

    if(src.channels() == 3)
    {
        cvtColor(src,imgHls,CV_BGR2HSV);
        split(imgHls,vHls);
        vHls[2].copyTo(img);
    }
    else
        src.copyTo(img);

    //DFT
    //cout<<"DFT "<<endl;
    Mat img_complex,img_mag,img_phase;
    Mat fpower = Fourier_Transform(img,img_complex,img_phase,img_mag);
    Shifting_DFT(img_complex);
    Shifting_DFT(fpower);
    //int D0 = getRadius(fpower,0.15);
    int D0 = 10;
    int n = 2;
    int w = img_complex.cols;
    int h = img_complex.rows;
    //BHPF
//  Mat filter,filter_complex;
//  filter = BHPF(h,w,D0,n);
//  Mat planes[] = {Mat_<float>(filter), Mat::zeros(filter.size(), CV_32F)};
//  merge(planes,2,filter_complex);

    int rH = 150;
    int rL = 50;
    Mat filter,filter_complex;
    filter_complex = Butterworth_Homomorphic_Filter(Size(w,h),D0,n,rH,rL,filter);

    //do mulSpectrums on the full dft
    mulSpectrums(img_complex, filter_complex, filter_complex, 0);

    Shifting_DFT(filter_complex);
    //IDFT
    Mat result;
    Inv_Fourier_Transform(filter_complex,result);

    if(src.channels() == 3)
    {
        vHls.at(2)= result(Rect(0,0,src.cols,src.rows));
        merge(vHls,imgHls);
        cvtColor(imgHls, dst, CV_HSV2BGR);
    }
    else
        result.copyTo(dst);

}

總結(jié)

以上是生活随笔為你收集整理的OpenCV同态滤波的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。