$\bf 摘要$: 本文给出了王大凯等编的《图像处理中的偏微分方程方法》第 5.4.1 节的详细论述.

$\bf 关键词$: 图像滤波; 方向扩散模型; matlab 编程

1. 模型的建立

从保护图像边缘的观点出发, 我们希望扩散是沿着平行于边缘的切线方向 (即垂直于 $\n I$ 的方向) 进行. 于是得到如下 PDE: $$\bee\label{1:df} I_t=I_{\xi\xi}, \eee$$ 其中 $\xi(\perp \n I)$ 为单位矢量. 我们化简 \eqref{1:df} 如下 (若 $$\bex {\bf \eta}=\frac{\n I}{\sev{\n I}} =(\cos t,\sin t),\quad {\bf \xi}=(-\sin t,\cos t)(\perp {\bf \eta}), \eex$$ 则 $$\bex \sex{\ba{cc} \frac{\p}{\p \xi}\\ \frac{\p}{\p \eta} \ea} &=&\sex{\ba{cc} \cos t&\sin t\\ -\sin t&\cos t \ea} \sex{\ba{cc} \frac{\p}{\p x}\\ \frac{\p}{\p y} \ea}, \eex$$ 而 $$\bex \frac{\p^2}{\p \xi^2} +\frac{\p^2}{\p \eta^2} &=&\sex{\ba{cc} \frac{\p}{\p \xi}&\frac{\p}{\p \eta} \ea} \sex{\ba{cc} \frac{\p}{\p \xi}\\ \frac{\p}{\p \eta} \ea} =\sex{\ba{cc} \frac{\p}{\p x}&\frac{\p}{\p y} \ea} \sex{\ba{cc} \cos t&-\sin t\\ \sin t&\cos t \ea} \sex{\ba{cc} \cos t&\sin t\\ -\sin t&\cos t \ea} \sex{\ba{cc} \frac{\p}{\p x}\\ \frac{\p}{\p y} \ea}\\ &=& \sex{\ba{cc} \frac{\p}{\p x}&\frac{\p}{\p y} \ea} \sex{\ba{cc} \frac{\p}{\p x}\\ \frac{\p}{\p y} \ea} =\frac{\p^2}{\p x^2}+\frac{\p^2}{\p y^2} \eex$$): $$\bex I_t=I_{\xi\xi}\\ &=&\lap I-I_{\eta\eta}\\ &=&\lap I-(\sev{\n I})_\eta\quad\sex{I_\eta=\n I\cdot{\bf \eta} =\n I\cdot\frac{\n I}{\sev{\n I}}=\sev{\n I}}\\ &=&\lap I-\n \sev{\n I}\cdot {\bf \eta}\\ &=&\lap I-\frac{ \sex{I_xI_{xx}+I_yI_{xy},I_xI_{xy}+I_yI_{yy}}}{\sev{\n I}} \cdot\frac{(I_x,I_y)}{|\n I|}\\ &=&I_{xx}+I_{yy} -\frac{I_x^2I_{xx}+2I_xI_yI_{xy}+I_y^2I_{yy}}{|\n I|^2}\\ &=&\frac{I_y^2I_{xx}-2I_xI_yI_{xy}+I_x^2I_{yy}}{|\n I|^2}. \eex$$ 即有 $$\bee\label{1:df_xy} I_t=\frac{I_y^2I_{xx}-2I_xI_yI_{xy}+I_x^2I_{yy}}{|\n I|^2}. \eee$$ 注意到 $I$ 的水平集 (level set) $$\bex \chi_\lambda(I)=\sed{(x,y);\ I(x,y)\geq \lambda} \eex$$ 适合 $$\bex \lambda_1\leq \lambda_2\ra \chi_{\lambda_2}(I)\subset \chi_{\lambda_1}(I), \eex$$ 而以 ${\bf \xi}$ 为切方向的曲线之法向量 $$\bex {\bf N}=\frac{\n I}{\sev{\n I}} ={\bf \eta}. \eex$$ 因此, 由 $$\bex \kappa =-\Div{\bf N} =-\Div\frac{\n I}{\sev{\n I}} =-\frac{I_y^2I_{xx}-2I_xI_yI_{xy}+I_x^2I_{yy}}{\sev{\n I}^3} \eex$$ 知 \eqref{1:df_xy} 可进一步化简为 $$\bee\label{1:df_cuv} I_t=-\kappa \sev{\n I}. \eee$$ 这说明沿水平集切方向扩散等价于水平集做曲率运动. 再注意到曲率运动与数学形态学 (mathematical morphology) 的联系, 我们知道沿水平集切方向扩散亦等价于对图像进行中值滤波. 从中我们体会到PDE 的好处: 不必对图像的所有水平集做分解---中值集操作 (数学形态学的一个算法) ---重构这三步骤就可以实现对灰度图像的中值滤波. 如果希望削弱在边缘附近的扩散速率, 还可引进边缘停止函数 (edge stopping function): $$\bee\label{1:esf} g(r)=\frac{1}{1+\sex{\frac{r}{K}}^p},\quad \sex{p=1,2}, \eee$$ 而将 \eqref{1:df} 改为 $$\bee\label{1:df_esf} I_t=g(\sev{\n I})I_{\xi\xi}, \eee$$ 或将 \eqref{1:df_xy} 改为 $$\bee\label{1:df_xy_esf} I_t=g(\sev{\n I})\frac{I_y^2I_{xx}-2I_xI_yI_{xy}+I_x^2I_{yy}}{|\n I|^2}, \eee$$ 抑或将 \eqref{1:df_cuv} 改为 $$\bee\label{1:df_cuv_esf} I_t=-\kappa g(\sev{\n I})\sev{\n I}. \eee$$

2. 数值算法

我们采用完全显式方案求解 \eqref{1:df_xy} 或 \eqref{1:df_xy_esf}： $$\bee\label{2:scheme_df} I_{ij}^{n+1} =I^n_{ij}+\lap t\cdot Q^n_{ij}, \eee$$ 其中 $$\bex Q^n_{ij} =\frac{ \sex{D^{(0)}_yI^n_{ij}}^2 \cdot D^{(0)}_{xx}I^n_{ij} -2D^{(0)}_xI^n_{ij}\cdot D^{(0)}_yI^n_{ij}\cdot D^{(0)}_{xy}I^n_{ij} +\sex{D^{(0)}_xI^n_{ij}}^2 \cdot D^{(0)}_{yy}I^n_{ij} }{ \sex{D^{(0)}_xI^n_{ij}}^2 +\sex{D^{(0)}_yI^n_{ij}}^2 }, \eex$$ $$\bex D^{(0)}_xI^n_{ij} =\frac{I^n_{i,j+1}-I^n_{i,j-1}}{2}, \quad D^{(0)}_yI^n_{ij} =\frac{I^n_{i+1,j}-I^n_{i-1,j}}{2}, \eex$$ $$\bex D^{(0)}_{xx}I^n_{xy} &=&D^{(0)}_xI^n_{i,j+\frac{1}{2}} -D^{(0)}_xI^n_{i,j-\frac{1}{2}}\\ &=&\sex{I^n_{i,j+1}-I^n_{i,j}}-\sex{I^n_{i,j}-I^n_{i,j-1}}\\ &=&I^n_{i,j+1}+I^n_{i,j-1}-2I^n_{ij}, \eex$$ $$\bex D^{(0)}_{xy}I^n_{ij} &=&D^{(0)}_xI^n_{i+\frac{1}{2},j} -D^{(0)}_xI^n_{i-\frac{1}{2},j}\\ &=&D^{(0)}_x\frac{I^n_{i+1,j}+I^n_{ij}}{2} -D^{(0)}_x\frac{I^n_{i-1,j}+I^n_{ij}}{2}\quad \sex{\mbox{用均值近似半点值}}\\ &=&\frac{\sex{I^n_{i+1,j+1}-I^n_{i+1,j-1}}+\sex{I^n_{i,j+1}-I^n_{i,j-1}}}{4}\\ & & -\frac{\sex{I^n_{i-1,j+1}-I^n_{i-1,j-1}}+\sex{I^n_{i,j+1}-I^n_{i,j-1}}}{4}\\ &=&\frac{I^n_{i+1,j+1}+I^n_{i-1,j-1}-I^n_{i+1,j-1}-I^n_{i-1,j+1}}{4}, \eex$$ $$\bex D^{(0)}_{yy}I^n_{ij} =I^n_{i+1,j}+I^n_{i-1,j}-2I^n_{ij}, \eex$$ 或者 $$\bee\label{2:scheme_df_esf} I^{n+1}_{ij}=I^n_{ij} +\lap t\cdot g\sex{I^n_{ij}} \cdot Q^n_{ij}. \eee$$

3. 数值实验

用如下的 matlab 代码

clear all;

close all;

clc;

%%参数设定

dt=0.05;%时间步长

N=1000;%迭代次数设置

D=200;%每运行D次输出图形

nn=1;%输出图形个数初始化

I=imread(‘key.bmp‘);

I=imnoise(I,‘salt & pepper‘,0.3);

I=double(rgb2gray(I));

[ny,nx]=size(I);

J=I;%方向扩散图像初始化

K=I;%带边缘函数的方向扩散图像初始化

for n=1:N

    Dx=0.5*(J(:,[2:nx,nx])-J(:,[1,1:nx-1]));

    Dy=0.5*(J([2:ny,ny],:)-J([1,1:ny-1],:));

    Dxx=(J(:,[2:nx,nx])-J)-(J-J(:,[1,1:nx-1]));

    Dxy=0.25*((J([2:ny,ny],[2:nx,nx])-J([2:ny,ny],[1,1:nx-1]))

                -(J([1,1:ny-1],[2:nx,nx])-J([1,1:ny-1],[1,1:nx-1])));

    Dyy=(J([2:ny,ny],:)-J)-(J-J([1,1:ny-1],:));

    J=J+dt*(Dy.^2.*Dxx-2*Dx.*Dy.*Dxy+Dx.^2.*Dyy)./(eps+Dx.^2+Dy.^2);

    Dx=0.5*(K(:,[2:nx,nx])-K(:,[1,1:nx-1]));

    Dy=0.5*(K([2:ny,ny],:)-K([1,1:ny-1],:));

    Dxx=(K(:,[2:nx,nx])-K)-(K-K(:,[1,1:nx-1]));

    Dxy=0.25*((K([2:ny,ny],[2:nx,nx])-K([2:ny,ny],[1,1:nx-1]))

                -(K([1,1:ny-1],[2:nx,nx])-K([1,1:ny-1],[1,1:nx-1])));

    Dyy=(K([2:ny,ny],:)-K)-(K-K([1,1:ny-1],:));

    K=K+dt./(1+(Dx.^2+Dy.^2)/1024)

                .*(Dy.^2.*Dxx-2*Dx.*Dy.*Dxy+Dx.^2.*Dyy)./(eps+Dx.^2+Dy.^2);

    if mod(n,D)==0

        figure(nn);

        subplot(1,3,1);

        imshow(uint8(I));

        subplot(1,3,2);

        imshow(uint8(J));

        subplot(1,3,3);

        imshow(uint8(K));

        hold off;

        nn=nn+1;

    end

end

就可实现方向扩散, 见下图(每幅图由三幅小图组成, 第一幅是原图像, 第二幅是方向扩散后的图像, 而第三幅是带边缘停止函数的方向扩散图像). 从中我们可以看到带边缘停止函数的方向扩散能更好地保持边缘的锐度.