图像增强处理之：同态滤波与Retinex算法(一)同态滤波

链接中文章，讲述了基于同态滤波和retinex算法进行图像去雾的算法，本篇对同态滤波和retinex算法进行介绍。

在生活中会得到这样的图像，动态范围很大，感兴趣的部分的灰度却很暗，范围很小，灰度层次和细节没有办法辨认，用一般的灰度线性变换法是不行的，因为扩展灰度级虽可以提高物体图像的 因为扩展灰度级虽可以提高物体图像的反差，但会使动态范围变大。而压缩灰 反差，但会使动态范围变大。而压缩灰度级，虽可以减少动态范围，但物体灰 度级，虽可以减少动态范围，但物体灰

度层次和细节就会更看不清。 度层次和细节就会更看不．这时我们需要采用同态滤波。同态滤波是一种在频域中同时将图像亮度范围进行压缩和将图像对比度进行增强的方法。

同态滤波

节选自百度百科里的一段话：同态滤波的基本原理是：将像元灰度值看作是照度和反射率两个组份的产物。由于照度相对变化很小，可以看作是图像的低频成份，而反射率则是高频成份。通过分别处理照度和反射率对像元灰度值的影响，达到揭示阴影区细节特征的目的。同态滤波处理的基本流程如下：
S(x,y)---->Log---->DFT---->频域滤波---->IDFT---->Exp---->T(x,y)
其中S(x,y)表示原始图像；T(x,y)表示处理后的图像；Log 代表对数运算；DFT 代表傅立叶变换（实际操作中运用快速傅立叶变换FFT）；IDFT 代表傅立叶逆变换（实际操作中运用快速傅立叶逆变换IFFT）；Exp
代表指数运算。



下面附一段更直观，带公式的讲解，作者:风吹夏天

对于一副图像f(x,y)可由照射分量i(x,y)和反射分量r(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)，它可用不同的可控方法影响傅里叶变换的低频和高频。如果γL和γH选定，而γL<1且γH>1，那么滤波器函数趋近于衰减低频（照射）的贡献，而增强高频反射的贡献。最终结果是同时进行动态范围的压缩和对比度的增强。



H(u,v)=(γH−γL)[1−e−c[D2(u,v)/D20]]+γL
同态滤波的目的目的：消除不均匀照度的影响而又不损失图象细节。其依据：图象的灰度由照射分量和反射分量合成。反射分量反映图象内容，随图象细节不同在空间上作快速变化。照射分量在空间上通常均具有缓慢变化的性质。照射分量的频谱落在空间低频区域，反射分量的频谱落在空间高频区。

参考1：

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;
img=imread('../image/1.jpg');
I3 = homofilter(img);

上面程序是将RGB图像灰度化之后再同态滤波的，在网上没找到直接针对彩色图像同态滤波处理的，所以我在上面程序基础上写了针对彩色图像的程序：

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('同态滤波增强后');

下面是我参照matlab程序在opencv的实现程序，针对灰度图像：

#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()
{
Mat src = imread("D:\\fcq_proMatlab\\vessel_edge_extration\\image\\homofilter.jpg",0);
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");
}
效果图如下：



参考：
http://blog.csdn.net/scottly1/article/details/42705271 http://blog.csdn.net/lilingyu520/article/details/46654265 http://blog.csdn.net/bluecol/article/details/45788803