• 同态滤波(Homomorphic filtering)是一种广泛用于信号和图像处理的技术，将原本的信号经由非线性映射，转换到可以使用线性滤波器的不同域，做完运算后再映射回原始域。同态的性质就是保持相关的属性不变，而同态滤波的好处是将原本复杂的运算转为效能相同但相对简单的运算。
• 同态滤波利用去除乘性噪声(multiplicative noise)，可以同时增加对比度以及标准化亮度，借此达到图像增强的目的

Homomorphic filtering

  照射-反射模型可用于开发一种频率域处理过程，该过程通过同时压缩灰度级范围和增强对比度来改善一幅图像的表现，这种方法被称为同态滤波。一幅图像 f(x,y)f(x,y)f(x,y) 可以表示为其照射分量 i(x,y)i(x,y)i(x,y) 和反射分量 r(x,y)r(x,y)r(x,y) 的乘积，即

(1)f(x,y)=i(x,y)r(x,y)\begin{aligned} f(x,y)=i(x,y)r(x,y)\tag{1} \end{aligned} f(x,y)=i(x,y)r(x,y)​(1)

同态滤波的滤波器函数有如下形式

(2)H(u,v)=(γH−γL)[1−e−c(D(u,v)2/D02)]+γL\begin{aligned} H(u,v)=(\gamma_H-\gamma_L)[1-e^{-c(D(u,v)^2/D_0^2)}]+\gamma_L \tag{2} \end{aligned} H(u,v)=(γH​−γL​)[1−e−c(D(u,v)2/D02​)]+γL​​(2)

其中，γH>1，γL<1\gamma_H>1\text{，}\gamma_L<1γH​>1，γL​<1，常数 ccc 控制函数坡度的锐利度。上式中，D0D_0D0​ 是一个正常数，D(u,v)D(u,v)D(u,v) 则是频率域中点 (u,v)(u,v)(u,v) 与频率矩形中心的距离，即

(3)D(u,v)=[(u−P/2)2+(v−Q/2)2]1/2\begin{aligned} D(u,v)=[(u-P/2)^2+(v-Q/2)^2]^{1/2}\tag{3} \end{aligned} D(u,v)=[(u−P/2)2+(v−Q/2)2]1/2​(3)

输入图像为 f(x,y)f(x,y)f(x,y)，同态滤波结果为 g(x,y)g(x,y)g(x,y) 基本步骤如下：
STEP1: 对数变换 z(x,y)=ln⁡f(x,y)=ln⁡i(x,y)+ln⁡r(x,y)z(x,y)=\ln f(x,y)=\ln i(x,y)+\ln r(x,y)z(x,y)=lnf(x,y)=lni(x,y)+lnr(x,y)[1]

STEP2: 傅里叶变换 Z(u,v)=F{z(x,y)}=F{ln⁡f(x,y)}=F{ln⁡i(x,y)}+F{ln⁡r(x,y)}Z(u,v)=\mathfrak{F}\{z(x,y)\}=\mathfrak{F}\{\ln f(x,y)\}=\mathfrak{F}\{\ln i(x,y)\}+\mathfrak{F}\{\ln r(x,y)\}Z(u,v)=F{z(x,y)}=F{lnf(x,y)}=F{lni(x,y)}+F{lnr(x,y)}

STEP3: 频域相乘 S(u,v)=H(u,v)Z(u,v)S(u,v)=H(u,v)Z(u,v)S(u,v)=H(u,v)Z(u,v)

STEP4: 傅里叶逆变换 s(x,y)=F−1S(u,v)s(x,y)=\mathfrak{F}^{-1}{S(u,v)}s(x,y)=F−1S(u,v)

STEP5: 指数变换 g(x,y)=es(x,y)g(x,y)=e^{s(x,y)}g(x,y)=es(x,y)

MATLAB实现

同态滤波函数homo_f如下：输入参数如式(2)(2)(2)中所述。

function Output=homo_f(I,rL,rH,c,D0)
    [M,N,C]=size(I);                % 图像矩阵维度
    P = 2*M; Q = 2*N;               % 图像填充后大小
    Output=log(double(I)+1);        % 对数变换
    D=dftuv(P,Q);                   % D(u,v)
    H=(rH-rL).*(1-exp(-c.*(D./D0^2)))+rL; % 滤波函数
    H=ifftshift(H);                 %对H做反中心化
    % 循环计算
    for i=1:C
        FI=fft2(Output(:,:,i),P,Q); %傅里叶变换 
        Ou1=ifft2(H.*FI);           %傅里叶逆变换
        Ou2 =Ou1(1:M, 1:N);         %恢复图像大小
        Ou3=exp(Ou2)-1;             %取指数
        Output(:,:,i)=Ou3;
    end
    Output=uint8(Output);
% D(u,v)计算函数，输入参数为图像大小
function D=dftuv(M,N)
    u=0:M-1;v=0:N-1;
    idx=find(u>M/2);u(idx)=u(idx)-M;
    idy=find(v>N/2);v(idy)=v(idy)-N;
    [U,V]=meshgrid(v,u);D=(U.^2+V.^2);
end
end
复制代码

Reference

• [1] Gonzalez R. Digital Image Processing.
• [2] Gonzalez R, Woods R, Eddins S. Digital Image Processing Using MATLAB.

1. 为了避免处理ln⁡(0)\ln(0)ln(0)的情形，处理时每个元素加1，在滤波结束时再加1 ↩︎