机器视觉_图像算法（四）——同态滤波

原理流程：

 

//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 返回功率谱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;
}

 

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);
}

 

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);

}

reference:

/*--------------------- 
https://blog.csdn.net/liujiabin076/article/details/53366678 
*/