# 【图像分割】GAC水平集方法图像分割【Matlab 090期】

1 原理分析

2 编程过程
2.1 准备工作

2.2 迭代运算

## 二、源代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
clc;

%读图像并将rgb图像转化为灰度图象，为了节省cpu时间，将图像变小为原来大小的1/2

im=rgb2gray(im);
im=double(im);
im = imresize( im, 0.5 );
figure(1);imshow(uint8(im));
[nrow,ncol]=size(im);

%为了计算函数g，先对图像作预平滑。
J= gauss( im,3,2 );   %guassian平滑

%求图像梯度模值
J_x = (J(:,[2:ncol ncol])-J(:,[1 1:ncol-1]))/2;
J_y = (J([2:nrow nrow],:)-J([1 1:nrow-1],:))/2;
%计算边缘函数g
kk=5;  %contrast parameter

%初始化曲线为一个图，初始化u为带符号距离函数
curvIndex=zeros(5*nrow*ncol,2);
curvImag=zeros([nrow,ncol]);
num=0;
for i = 2 : nrow - 1
for j = 2 : ncol - 1
if  ((i==8|i==nrow-8)&(j>=8&j<=(ncol-8)))|((i>=8&i<=(nrow-8))&(j==8|j==ncol-8))
num=num+1;
InitCurvImag(i,j)=255;
curvIndex(num,1)=i;
curvIndex(num,2)=j;
for k=1:num-1              %check is it the last point of the curve
if curvIndex(k,1)==i & curvIndex(k,2)==j
num=num-1;break
end
end
end
end
end

u = zeros(nrow,ncol);
dist=zeros(1,num);
for j=1:ncol
for i=1:nrow
for k=1:num
dist(k)=sqrt((i-curvIndex(k,1)).^2+(j-curvIndex(k,2)).^2);
end
u(i,j)=min(dist);
if i>8&i<(nrow-8)&j>8&j<(ncol-8)
u(i,j)=-u(i,j);
end
end
end

%将当前曲线加入到原图像中,然后写入文件
newim=createimage(im,u,0);
figure(2);imshow(uint8(newim));

%选定迭代步长
delta_t=20;
%选定常数速度
c=0.2;

%迭代开始
for iterations=1:4000
u_x_e = u(:,[2:ncol,ncol])-u;
u_y_e = (u([2:nrow,nrow],:)+u([2:nrow,nrow],[2:ncol,ncol])-u([1 1:nrow-1],:)-u([1 1:nrow-1],[2:ncol ncol]))/4;
u_G_e = sqrt(u_x_e .^2+u_y_e .^2);
g_e = 0.5*(g(:,[2:ncol,ncol])+g);
Term_e = g_e.*u_x_e./(u_G_e+eps);

u_x_w  = u-u(:,[1 1:ncol-1]);
u_y_w = (u([2:nrow,nrow],:)+u([2:nrow,nrow],[1 1:ncol-1])-u([1,1:nrow-1],:)-u([1,1:nrow-1],[1 1:ncol-1]))/4;
u_G_w = sqrt(u_x_w.^2+u_y_w.^2);
g_w = 0.5*(g(:,[1 1:ncol-1])+g);
Term_w = g_w.*u_x_w./(u_G_w+eps);

u_y_s = u([2:nrow,nrow],:)-u;
u_x_s = (u(:,[2:ncol,ncol])+u([2:nrow,nrow],[2:ncol,ncol])-u(:,[1 1:ncol-1])-u([2:nrow,nrow],[1 1:ncol-1]))/4;
u_G_s = sqrt(u_y_s.^2+ u_x_s.^2);
g_s = 0.5*(g([2:nrow,nrow],:)+g);
Term_s = g_s.*u_y_s./(u_G_s+eps);

u_y_n = u-u([1 1:nrow-1],:);
u_x_n = (u(:,[2:ncol,ncol])+u([1 1:nrow-1],[2:ncol,ncol])-u(:,[1 1:ncol-1])-u([1 1:nrow-1],[1 1:ncol-1]))/4;
u_G_n = sqrt(u_y_n.^2+u_x_n.^2);
g_n = 0.5*(g([1,1:nrow-1],:)+g);
Term_n = g_n.*u_y_n./(u_G_n+eps);
function Ig=gauss(I,ks,sigma2)
%private function: gauss (by Guy Gilboa):
% Ig=gauss(I,ks,sigma2)
% ks - kernel size (odd number)
% sigma2 - variance of Gaussian

[Ny,Nx]=size(I);
hks=(ks-1)/2;  % half kernel size
if (Ny<ks)   % 1d convolutin
x=(-hks:hks);
flt=exp(-(x.^2)/(2*sigma2));  % 1D gaussian
flt=flt/sum(sum(flt));  % normalize
% expand
x0=mean(I(:,1:hks)); xn=mean(I(:,Nx-hks+1:Nx));
eI=[x0*ones(Ny,ks) I xn*ones(Ny,ks)];
Ig=conv(eI,flt);
Ig=Ig(:,ks+hks+1:Nx+ks+hks);  % truncate tails of convolution
else
%% 2-d convolution
x=ones(ks,1)*(-hks:hks); y=x';
flt=exp(-(x.^2+y.^2)/(2*sigma2));  % 2D gaussian
flt=flt/sum(sum(flt));  % normalize
% expand
if (hks>1)
xL=mean(I(:,1:hks)')'; xR=mean(I(:,Nx-hks+1:Nx)')';
else
xL=I(:,1); xR=I(:,Nx);
end
eI=[xL*ones(1,hks) I xR*ones(1,hks)];
if (hks>1)
xU=mean(eI(1:hks,:)); xD=mean(eI(Ny-hks+1:Ny,:));
else
xU=eI(1,:); xD=eI(Ny,:);
end


04-07

12-20
09-24
08-11 2352
07-03 1万+
10-16 1632
09-14