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

$\bf 关键词$: 图像分割; 活动轮廓模型; matlab 编程

1 模型的建立

在图像中, 对象与背景的区别有时表现为平均灰度的明显不同. 由于这类图像既没有明显的边缘 ($\sev{\n I}$ 大), 也缺乏明显的纹理 (texture, 灰度变化有一定的规律, 并形成一定的 patten), 故测地线活动轮廓 (geodesic active contour, GAC, 或 snake) 模型不能成功地实现图像分割.

为此, T. Chan 与 L. Vese 提出了如下的能量泛函: $$\bee\label{1:E} E(c_1,c_2,C)\equiv \mu\oint_C\rd s +\lambda_1\iint_{\Omega_1} (I-c_1)^2\rd x\rd y +\lambda_2\iint_{\Omega_2} (I-c_2)^2\rd x\rd y, \eee$$ 其中

(1) $\mu>0$, $\lambda_1$, $\lambda_2>0$ 是各比分的权重;

(2) 曲线 $C$ 是简单的, 并将这个图像区域 $\Omega$ 分为内外两部分, 分别记为 $\Omega_1$, $\Omega_2$;\

(3) $c_1$, $c_2$ 是两常数, 表分片常数图像.

此模型的特点是

(1)    一方面使边界曲线最短;

(2)    另一方面使背景 $(\Omega_2)$ 与对象 $(\Omega_1)$ 区别开来 (\eqref{1:E} 中的后两项). 这称作无边缘活动轮廓 (active contour without edge, C-V, 或 GAR) 模型.

引入 $\delta$ 函数, 将 $E$ 写成 $$\bee\label{1:Eu}\bea E(c_1,c_2,u) &=\mu\iint_\Omega \delta(u)\sev{\n u}\rd x\rd y\\ &\quad+\lambda_1\iint_{\Omega_1}(I-c_1)^2\sez{1-H(u)}\rd x\rd y +\lambda_2\iint_{\Omega_2}(I-c_2)^2H(u)\rd x\rd y, \eea\eee$$ 其中

(1) $\delta(\cdot)$ 为 Dirac $\delta$ 函数, 其有光滑逼近元 $$\bex \delta_\ve(x)=\frac{1}{\pi}\cdot\frac{\ve}{x^2+\ve^2}\quad \sex{0<\ve\ll 1}; \eex$$

(2) $\sev{\n u}$ 是变换之 Jacobian;

(3) $H(\cdot)$ 是 Heaviside 函数: $$\bex H(z)=\left\{\ba{ll} 1,&z\geq 0,\\ 0,&z<0, \ea\right. \eex$$ 其有光滑逼近元 $$\bex H_\ve(z)=\frac{1}{2}+\frac{1}{\pi}\arctan \frac{z}{\ve}, \eex$$ 且 $$\bex H‘(z)=\delta(z),\quad\mbox{ 在 }\mathcal{D}‘(\bbR)\mbox{ 上}. \eex$$   

(4) 在函数 $u$ 固定的条件下, 由 $$\bex \iint_{\Omega_1} (I-c_1)^2\rd x\rd y =c_1^2\iint_{\Omega_1}\rd x\rd y -2c_1\iint_{\Omega_1}I\rd x\rd y +\iint_{\Omega_1}I^2\rd x\rd y, \eex$$ $$\bex \iint_{\Omega_2} (I-c_2)^2\rd x\rd y =c_2^2\iint_{\Omega_2}\rd x\rd y -2c_2\iint_{\Omega_2}I\rd x\rd y +\iint_{\Omega_2}I^2\rd x\rd y, \eex$$ 及二次函数的性质知当且仅当 $$\bee\label{1:c1c2} c_1=\frac{\iint_{\Omega_1}I\rd x\rd y}{\iint_{\Omega_1}\rd x\rd y},\quad c_2=\frac{\iint_{\Omega_2}I\rd x\rd y}{\iint_{\Omega_2}\rd x\rd y} \eee$$ 时, $E$ 取得最小.

现在, 写出 \eqref{1:Eu} 的变分 (variation) 如下 $$\bex \frac{\delta E}{\delta u} =-\mu \delta(u)\Div\sex{\frac{\n u}{\sev{\n u}}} -\lambda_1\delta(u)(I-c_1)^2 +\lambda_2\delta(u)(I-c_2)^2. \eex$$

从而梯度下降流为 $$\bee\label{1:gdf} \frac{\p u}{\p t} =-\frac{\delta E}{\delta u} =\mu \delta(u)\Div\sex{\frac{\n u}{\sev{\n u}}} +\lambda_1\delta(u)(I-c_1)^2 -\lambda_2\delta(u)(I-c_2)^2. \eee$$

2.数值算法

为计算方便, 须离散化, 而用 $\delta_\ve$ 替代 $\delta$. 这样, 用半隐式 (semi-implicit) 方案: $$\bee\label{2:scheme} u^{n+1}_{ij} =u^n_{ij} +\tau \delta_\ve(u^n_{ij}) \sez{\mu Q(u^{n+1}_{ij})+(I_{ij}-c^n_1)^2-(I_{ij}-c^n_2)^2}, \eee$$ 其中已选 $\lambda_1=\lambda_2=1$, $Q(u^{n+1}_{ij})$ 采用向前向后差分相结合的格式: $$\bex Q(u_{ij}) &=D^{(+)}_x\sex{g_{ij}D^{(-)}_x u_{ij}} +D^{(+)}_y\sex{g_{ij}D^{(-)}_yu_{ij}}\quad\sex{g_{ij}=\frac{1}{\sev{\n u}_{ij}}}\\ &=g_{i,j+1}D^{(-)}_xu_{i,j+1} -g_{i,j}D^{(-)}_xu_{ij} +g_{i+1,j}D^{(-)}_yu_{i+1,j} -g_{ij}D^{(-)}_yu_{ij}\\ &=g_{i,j+1}\sex{u_{i,j+1}-u_{ij}} -g_{ij}\sex{u_{ij}-u_{i,j-1}}\\ &  +g_{i+1,j}\sex{u_{i+1,j}-u_{ij}} -g_{ij}\sex{u_{ij}-u_{i-1,j}}\\ &=\sez{g_{i,j+1}u_{i,j+1}+g_{ij}u_{i,j-1} +g_{i+1,j}u_{i+1,j}+g_{ij}u_{i-1,j}}\\ & -\sex{g_{i,j+1}+g_{ij}+g_{i+1,j}+g_{i,j}}u_{ij}, \eex$$ 且为保证计算的稳定性, 我们采用半隐式方案求解: $$\bee\label{2:Q}\bea Q(u^{n+1}_{ij}) &=\sez{g^n_{i,j+1}u^n_{i,j+1}+g^n_{ij}u^n_{i,j-1} +g^n_{i+1,j}u^n_{i+1,j}+g^n_{ij}u^n_{i-1,j}}\\ & +\sex{g^n_{i,j+1}+g^n_{ij}+g^n_{i+1,j}+g^n_{i,j}}u^{n+1}_{ij},  \eea\eee$$ 其中 $g^n_{ij}$ 的选取如下 (依赖于其对应的 $u_{ij}$, 一个方向上用向前或向后差分, 另一个方向上用中心差分):

(1)$u^n_{i,j+1}$ 的系数 $$\bex C_1\equiv g_{i,j+1} =\frac{1}{\sqrt{\sex{ u_{i,j+1}-u_{ij} }^2+\sex{\frac{ u_{i+1,j+1}-u_{i-1,j+1} }{2}}^2}}; \eex$$

(2)$u^n_{i,j-1}$ 的系数 $$\bex C_1\equiv g_{i,j+1} =\frac{1}{\sqrt{\sex{ u_{ij}-u_{i,j-1} }^2+\sex{\frac{ u_{i+1,j-1}-u_{i-1,j-1} }{2}}^2}}; \eex$$

(3)$u^n_{i+1,j}$ 的系数 $$\bex C_1\equiv g_{i,j+1} =\frac{1}{\sqrt{\sex{ u_{i+1,j}-u_{ij} }^2+\sex{\frac{ u_{i+1,j+1}-u_{i+1,j-1} }{2}}^2}}; \eex$$

(4)$u^n_{i-1,j}$ 的系数 $$\bex C_1\equiv g_{i,j+1} =\frac{1}{\sqrt{\sex{ u_{ij}-u_{i-1,j} }^2+\sex{\frac{ u_{i-1,j+1}-u_{i-1,j-1} }{2}}^2}}. \eex$$ 把 \eqref{2:Q} 代入 \eqref{2:scheme}, 得 $$\bex & \sez{1+\mu\tau \delta_\ve(u^n_{ij})} u^{n+1}_{ij}=u^n_{ij} +\tau \delta_\ve(u^n_{ij})\cdot\\ & \quad\cdot \sez{ \mu\sex{ C_1u^n_{i,j+1} +C_2u^n_{i,j-1} +C_3u^n_{i+1,j} +C_4u^n_{i-1,j} } +(I_{ij}-c^n_1)^2 -(I_{ij}-c^n_2)^2 }, \eex$$ 其中 $c^n_1$, $c^n_2$ 由 \eqref{1:c1c2} 离散化而来: $$\bex c^n_1 =\frac{ \sum_{ij}I_{ij}\sez{1-H_\ve\sex{u^n_{ij}}} }{ \sum_{ij}\sez{1-H_\ve\sex{u^n_{ij}}} },\quad c^n_2 =\frac{ \sum_{ij}I_{ij}H_\ve\sex{u^n_{ij}} }{ \sum_{ij}H_\ve\sex{u^n_{ij}} }. \eex$$

3. 数值试验

选取 $\ve=1.0$, $\mu=250$, $\tau=0.1$, 则通过如下的 matlab 代码: clear all;

close all;

clc;

Img=imread(‘brain.bmp‘);

Img=double(RGB2gray(Img));

figure(1);

imshow(uint8(Img));

[ny,nx]=size(Img);%获取图像大小

%%将初始曲线设为圆

%初始圆的圆心

c_i=floor(ny/2);

c_j=floor(nx/2);

%初始圆的半径

r=c_i/3;

%%初始化u为距离函数

u=zeros([ny,nx]);

for i=1:ny

    for j=1:nx

        u(i,j)=r-sqrt((i-c_i).^2+(j-c_j).^2);

    end

end

%%将初始圆形曲线叠加在原始图片上

figure(2);

imshow(uint8(Img));

hold on;

[c,h]=contour(u,[0 0],‘r‘);

%%初始化参数

epsilon=1.0;%Ｈeaviside函数参数设置

mu=250;%弧长的权重

dt=0.1;%时间步长

nn=0;%输出图像个数初始化

%%迭代计算开始

for n=1:1000

    %%计算正则化的Heaviside函数

    H_u=0.5+1/pi.*atan(u/epsilon);

    %%由当前u计算出参数c1,c2

    c1=sum(sum((1-H_u).*Img))/sum(sum(1-H_u));

    c2=sum(sum(H_u.*Img))/sum(sum(H_u));

    %%用当前c1,c2更新u

    delta_epsilon=1/pi*epsilon/(pi^2+epsilon^2);%delta函数的正则化

    m=dt*delta_epsilon;%临时参数m存储时间步长与delta函数的乘积

    %%计算四邻点的权重

    C1=1./sqrt(eps+(u(:,[2:nx,nx])-u).^2

        +0.25*(u([2:ny,ny],:)-u([1,1:ny-1],:)).^2);

    C2=1./sqrt(eps+(u-u(:,[1,1:nx-1])).^2

        +0.25*(u([2:ny,ny],[1,1:nx-1])-u([1,1:ny-1],[1,1:nx-1])).^2);

    C3=1./sqrt(eps+(u([2:ny,ny],:)-u).^2

        +0.25*(u([2:ny,ny],[2:nx,nx])-u([2:ny,ny],[1,1:nx-1])).^2);

    C4=1./sqrt(eps+(u-u([1,1:ny-1],:)).^2

        +0.25*(u([1,1:ny-1],[2:nx,nx])-u([1,1:ny-1],[1,1:nx-1])).^2);

    C=1+mu*m.*(C1+C2+C3+C4);

    u=(u+m*

        (mu*(C1.*u(:,[2:nx,nx])+C2.*u(:,[1,1:nx-1])

            +C3.*u([2:ny,ny],:)+C4.*u([1,1:ny-1],:))

        +(Img-c1).^2-(Img-c2).^2)

       )./C;

    %%每运行两百次显示曲线和分片常数曲线

    if mod(n,200)==0

        nn=nn+1;

        f=Img;

        f(u>0)=c1;

        f(u<0)=c2;

        figure(nn+2);

        subplot(1,2,1);imshow(uint8(f));

        subplot(1,2,2);imshow(uint8(Img));

        hold on;

        [c,h]=contour(u,[0,0],‘r‘);

        hold off;

    end

end

就可实现对人的大脑切片的分割, 见下图:

注记:

1.为书写方便, matlab 程序中有的部分已经分断开来, 比如 $C_1,\cdots,C_4$ 及 $u$, 您在书写的时候可不能分段哦.

2.从上述程序可以看到在 matlab 中多用矩阵运算比单纯的迭代运算起来快得多.

4. 致谢

The author would like to thank Dr. S.J. Li and Dr. P. Li for suggesting the book by D.K. Wang and etc. And special thank goes to Dr. P. Li for exploiting this matlab procedure.