1GAC_fenge_gac_图像分割_

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This program impliment general GAC model： % ut=m(div(|▽u|-▽u/|▽u|)+s(u)(div(g▽u/|▽u|)+cg) % by level_set method with upwind scheme.It will call functions: % createimage(),which put the current zero_level_set (Curve) onto the % original image, and reinitialize the embeding functin u(). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear; clc; %读图像并将rgb图像转化为灰度图象，为了节省cpu时间，将图像变小为原来大小的1/2 im=imread('3.bmp'); %im=rgb2gray(im); im=double(im(:,:,1)); sigma=1.5;% scale parameter in Gaussian kernel for smoothing. G=fspecial('gaussian',15,sigma); im=conv2(im,G,'same'); % smooth image by Gaussiin convolution %im = imresize( im, 0.5 ); figure(1);imagesc(im);colormap(gray);hold on; [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; grad_im = (J_x.^2 + J_y.^2).^0.5; %计算边缘函数g kk=5; %contrast parameter g=1./(1+(grad_im/kk).^2); %%%% Initialize U u = zeros(nrow,ncol); p=2;w=8; for i=1:nrow for j=1:ncol if ((i==w|i==nrow-w)&(j>=w&j<=(ncol-w)))|((i>=w&i<=(nrow-w))&(j==w|j==ncol-w)) u(i,j)=0;%在曲线上 elseif i>w&i<(nrow-w)&j>w&j<(ncol-w) u(i,j)=-2*p;%在曲线内 else u(i,j)=2*p;%在曲线外 end end end %将当前曲线加入到原图像中,然后写入文件 newim=createimage(im,u,0); figure(2);imagesc(newim,[0,255]);colormap(gray);hold on; %选定迭代步长 delta_t=5; %选定常数速度 c=0.15; %迭代开始 for iterations=1:1000 %%% compute Roe upwind gradient 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); Termq_e =u_x_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); Termq_w = u_x_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); Termq_s = u_y_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); Termq_n =u_y_n-u_y_n./(u_G_n+eps); divgn = Term_e - Term_w + Term_s - Term_n; divgnq = Termq_e - Termq_w + Termq_s - Termq_n; m=0.0001; s=p./(pi.*(p^2+u.^2)); term = divgn+c*g; termp=s.*term; termq=m*divgnq; u=u+delta_t*(termp+termq); if(mod(iterations,100)==0) newim=createimage(im,u,0 ); figure(3);imagesc(newim,[0,255]);colormap(gray);hold on; end end