图像分割基于GAC水平集方法实现图像分割matlab 源码

Posted 2021-05-24 Matlab走起

一、简介

1 原理分析




2 编程过程
2.1 准备工作

2.2 迭代运算

二、源代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
clc;
 
%读图像并将rgb图像转化为灰度图象，为了节省cpu时间，将图像变小为原来大小的1/2
 
im=imread('3.bmp');
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;      
grad_im = (J_x.^2 + J_y.^2).^0.5;      
%计算边缘函数g
kk=5;  %contrast parameter
g=1./(1+(grad_im/kk).^2); 
 
%初始化曲线为一个图，初始化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
   %%%  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);
    
   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

三、运行结果

四、备注

完整代码或者仿真咨询添加QQ1575304183