MATLAB code: Explicit finite difference method for Image Segmentation

%Main Code

clear all;

f = imread(‘blood_vessel.gif’);

f = double(f);

f0 = (f - min(min(f)))/(max(max(f))-min(min(f)));

[Nx,Ny] = size(f0); h = 1.0; dt = 0.1hh;

x = linspace(-0.5h,hNx+0.5*h,Nx+2);

y = linspace(-0.5h,hNy+0.5*h,Ny+2);

eps2 = h^2; lambda = 10.0;

oldphi(1:Nx+2,1:Ny+2)=0;

oldphi(2:Nx+1,2:Ny+1)=2*f0-1;

newphi=oldphi;

for iter = 1:25

c1 = sum(sum(f0.*(1.0+oldphi(2:Nx+1,2:Ny+1))))…

    /sum(sum(1.0+oldphi(2:Nx+1,2:Ny+1)));

c2 = sum(sum(f0.*(1.0-oldphi(2:Nx+1,2:Ny+1))))...

    /sum(sum(1.0-oldphi(2:Nx+1,2:Ny+1)));

oldphi(1,:)=oldphi(2,:);oldphi(Nx+2,:)=oldphi(Nx+1,:);

oldphi(:,1)=oldphi(:,2);oldphi(:,Ny+2)=oldphi(:,Ny+1);

for i = 2:Nx+1

for j = 2:Ny+1

newphi(i,j) = oldphi(i,j)+dt*((oldphi(i,j) …

    -oldphi(i,j)^3)/eps2+(oldphi(i-1,j) ...

    +oldphi(i+1,j)+oldphi(i,j-1)+oldphi(i,j+1) ...

    -4.0*oldphi(i,j))/h^2 ...

    -lambda*( (1.0+oldphi(i,j))*(f0(i-1,j-1)-c1)^2 ...

    -(1.0-oldphi(i,j))*(f0(i-1,j-1)-c2)^2));

end

end

oldphi=newphi;

figure(1);clf;hold on;surf(x(2:Nx+1),y(2:Ny+1),f0’);

shading interp;colormap gray; axis image;

axis([0 Nxh 0 Nyh -1 1]);axis off;

contour(x(2:Nx+1),y(2:Ny+1),oldphi(2:Nx+1,2:Ny+1)’, [0 0],’color’,’r’,’linewidth’,2);

view(180,-90);pause(0.2)

end