# 同态滤波

S(x,y)---->Log---->DFT---->频域滤波---->IDFT---->Exp---->T(x,y)

f(x,y)=i(x,y)r(x,y)

Inf(x,y)=lni(x,y)+lnr(x,y)

I{Inf(x,y)}=I{Ini(x,y)}+I{Inr(x,y)}

H(u,v)=(γHγL)[1ec[D2(u,v)/D20]]+γL

function I3 = homofilter(I)    %同态滤波函数
subplot(1,2,1),imshow(I);title('同态滤波之前原始图像');
I=double(rgb2gray(I));
[M,N]=size(I);
rL=0.5;
rH=5;%可根据需要效果调整参数
c=3;
d0=9;
I1=log(I+1);%取对数
FI=fft2(I1);%傅里叶变换
n1=floor(M/2);
n2=floor(N/2);
for i=1:M
for j=1:N
D(i,j)=((i-n1).^2+(j-n2).^2);
H(i,j)=(rH-rL).*(exp(c*(-D(i,j)./(d0^2))))+rL;%高斯同态滤波
end
end
I2=ifft2(H.*FI);%傅里叶逆变换
I3=real(exp(I2));
subplot(1,2,2),imshow(I3,[]);title('同态滤波增强后');  

%% 同台滤波测试函数
clc,close all,clear all;
I3 = homofilter(img);

function I3 = homofilterColor(I)
subplot(1,2,1),imshow(I);title('同态滤波之前原始图像');
% I=double(rgb2gray(I));
I=double(I);
[M,N,channals]=size(I);
rL=0.5;
rH=5;%可根据需要效果调整参数
c=3;
d0=9;
I3=I;
for ch=1:channals
I1=log(I(:,:,ch)+1);%取对数
FI=fft2(I1);%傅里叶变换
n1=floor(M/2);
n2=floor(N/2);
for i=1:M
for j=1:N
D(i,j)=((i-n1).^2+(j-n2).^2);
H(i,j)=(rH-rL).*(exp(c*(-D(i,j)./(d0^2))))+rL;%高斯同态滤波
end
end
I2=ifft2(H.*FI);%傅里叶逆变换
I3(:,:,ch)=real(exp(I2));
end

minV=min(min(min(I3)));
maxV=max(max(max(I3)));
for ch=1:channals
for i=1:M
for j=1:N
I3(i,j,ch)=255* (I3(i,j,ch)-minV)./(maxV-minV);
end
end
end

subplot(1,2,2),imshow(uint8(I3));title('同态滤波增强后');  

#include <opencv2\opencv.hpp>
#include <iostream>
#include <math.h>

using namespace cv;
using namespace std;
//同态滤波
void my_HomoFilter(Mat srcImg, Mat &dst)
{
srcImg.convertTo(srcImg, CV_64FC1);
dst.convertTo(dst, CV_64FC1);
//第一步，取对数
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] + 1);
}
}
//第二步，傅里叶变换
Mat mat_dct = Mat::zeros(srcImg.rows, srcImg.cols, CV_64FC1);
dct(srcImg, mat_dct);
namedWindow("dct", CV_WINDOW_FREERATIO);
imshow("dct", mat_dct);
//第三步，频域滤波
Mat H_u_v;
int n1 = floor(srcImg.rows / 2);
int n2 = floor(srcImg.cols / 2);
double gammaH = 5;
double gammaL = 0.5;
double C = 3;
double d0 = 9.0;
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-n1, 2.0) + pow(j-n2, 2.0);
dataH_u_v[j] = (gammaH - gammaL)*exp(C*(-d2 / (d0*d0))) + gammaL;
}
}
H_u_v.ptr<double>(0)[0] = 1.5;
mat_dct = mat_dct.mul(H_u_v);
//第四步，傅里叶逆变换
idct(mat_dct, dst);

//第五步，取指数运算
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]);
}
}
dst.convertTo(dst, CV_8UC1);
}

void main()
{
namedWindow("src:", CV_WINDOW_FREERATIO);
imshow("src:", src);
Mat dst(src.rows, src.cols, src.type());
my_HomoFilter(src, dst);
namedWindow("dst:", CV_WINDOW_FREERATIO);
imshow("dst:", dst);
waitKey();
system("pause");
}

1. http://blog.csdn.net/scottly1/article/details/42705271
2. http://blog.csdn.net/lilingyu520/article/details/46654265
3. http://blog.csdn.net/bluecol/article/details/45788803