同态滤波

// homographicfilter.cpp : 定义控制台应用程序的入口点。
//


#include "stdafx.h"
#include <iostream>
#include <highgui.h>
#include <cv.h>


using namespace cv;
using namespace std;


void shiftDFT(cv::Mat& fImage);
void create_butterworth_pass_filter(cv::Mat &dft_Filter, int D, int n);




int _tmain(int argc, _TCHAR* argv[])
{   
cv::Mat src =  imread("1.tif", 0);
//cout<<src.rows<<" "<<src.cols<<endl;
if (src.empty())
{
cout<<"read image error!!!"<<endl;
}
int radius = 50; // pass filter parameter
int order = 3; // pass filter parameter
cv::Mat doubleImg, logImg;
src.convertTo(doubleImg, CV_32FC1, 1, 0);


//求对数
log(doubleImg, logImg);


int M = getOptimalDFTSize(src.rows);
    int N = getOptimalDFTSize(src.cols);
cv::Mat padded;
copyMakeBorder(doubleImg, padded, 0, M - src.rows, 0, N - src.cols, BORDER_CONSTANT, cv::Scalar::all(0));
cv::Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
cv::Mat complexImg;
    merge(planes, 2, complexImg);
dft(complexImg, complexImg);


cv::Mat filter = complexImg.clone();
create_butterworth_pass_filter(filter, radius, order);
shiftDFT(complexImg);
    mulSpectrums(complexImg, filter, complexImg, 0);


    shiftDFT(complexImg);
idft(complexImg, complexImg);
split(complexImg, planes);


cv::Mat imgOutput, filterOutput;
    normalize(planes[0], imgOutput, 0, 255, CV_MINMAX);
    split(filter, planes);
    normalize(planes[0], filterOutput, 0, 255, CV_MINMAX);


imwrite("dft.jpg", imgOutput);
imwrite("lpdft.jpg", filterOutput);










return 0;
}


void create_butterworth_pass_filter(cv::Mat &dft_Filter, int D, int n)
{
Mat tmp = Mat(dft_Filter.rows, dft_Filter.cols, CV_32F);


Point centre = Point(dft_Filter.rows / 2, dft_Filter.cols / 2);
double radius;


// based on the forumla in the IP notes (p. 130 of 2009/10 version)
// see also HIPR2 on-line




for(int i = 0; i < dft_Filter.rows; i++)
{   
float* ptr = tmp.ptr<float>(i);
for(int j = 0; j < dft_Filter.cols; j++)
{
radius = (double) sqrt(pow((i - centre.x), 2.0) + pow((double) (j - centre.y), 2.0));
//tmp.at<float>(i,j) = 1.2+(float)( 1 / (1 + pow((double) (D / radius), (double) (2 * n))));//high pass
//tmp.at<float>(i,j) = (float)( 1 / (1 + pow((double) (radius / D), (double) (2 * n))));//low pass
ptr[j] = (float)( 1 / (1 + pow((double) (radius / D), (double) (2 * n))));//low pass
//ptr[j] = (float)( 1 / (1 + pow((double) (D / radius), (double) (2 * n))));//high pass
}
}


cv::Mat toMerge[] = {tmp, tmp};
merge(toMerge, 2, dft_Filter);
}


void shiftDFT(cv::Mat& fImage)
{
  Mat tmp, q0, q1, q2, q3;


// first crop the image, if it has an odd number of rows or columns
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));


q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);


q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
}