% This code is to implement the piecewise constant Chan-Vese model,but
% featured with two properties:
% 1. without reinitialization proposed by Chunming Li's paper presented at
% CVPR'2005,"Level Set Evolution Without Re-initialization: A New Variational Formulation"
% in Proceedings of CVPR'05, vol. 1, pp. 430-436.
% 2. local binary fitting also proposed by chunming li's paper, CVPR2007
% "Implicit Active Contours Driven By Local Binary Fitting Energy" in Proceedings of CVPR'07
% 3. the C-V model is incorporated with geodesic length, 2009/04/07
% (1) R. Kimmel, "Fast edge integration," in Geometric Level Set Methods in
% Imaging, Vision, and Graphics, S. Osher and N. Paragios, Eds. New
% York: Springer-Verlag, 2003, pp. 59每78.
% (2) R. Kimmel and A.M. Bruckstein. Regulariztion laplacian zero
% crossingas as optimal edge integrators, In proceedings of Image and
% Vision computing, IVCNZ01, New Zealand, Nov. 2001
% (3) R. Kimmel and A.M. Bruckstein.On edge detection edge integration
% and geometric active contours, In Proceedings of Int. Symposium on
% Mahtematical Morphology, ISMM2002, Sydney, New South Wales,
% Australia, Apr. 2002
% (4) Li Chen,Yue Zhou,YonggangWang, JieYang, GACV: Geodesic-Aided
% C每V method, Pattern Recognition 39 (2006) 1391 每 1395
% (5) Thomas Brox and Daniel Cremers, On the Statistical Interpretation of
% the Piecewise Smooth Mumford-Shah Functional, In Scale Space and Variational Methods in Computer Vision, Springer LNCS 4485,
% F. Sgallari et al. (Eds.), pp. 203-213, May 2007.
% 4. we will further improve this program.
% by Yuanquan Wang, CV center, TJUT
% 2008/06/23/
%
function PCCVnoReInitLBF(action,varargin)
if nargin<1,
action='Initialize';
end;
feval(action,varargin{:});
return;
function Initialize()
global HDmainf;
scrsz = get(0,'ScreenSize');
ww = scrsz(3);
hh = scrsz(4);
startp = ww*0.25;
endp = hh*0.25;
fw = scrsz(3)*1/2;
fh = scrsz(4)*0.5;
HDmainf = figure('Position',[startp endp fw fh],...
'Color', [0.9 0.9 0.9], ...
'NumberTitle', 'off', ...
'Name', 'Piecewise Constant Chan-Vese model without ReInit & with LBF- A demo', ...
'Units', 'pixels');
fpath = '../image/testimage/';
orgname = 'circle2';
fname1 = strcat(fpath,orgname); fname = strcat(fname1,'.bmp');
orgname = 'ring';
fname1 = strcat(fpath,orgname); fname = strcat(fname1,'.bmp');
% orgname = 'blob-like';
% fname1 = strcat(fpath,orgname); fname = strcat(fname1,'.bmp');
% orgname = 'twoObj';
% fname1 = strcat(fpath,orgname); fname = strcat(fname1,'.bmp');
%
% orgname = 'fig8';
% fname1 = strcat(fpath,orgname); fname = strcat(fname1,'.bmp');
% orgname = 'noisyNonUniform';
% fname1 = strcat(fpath,orgname); fname = strcat(fname1,'.bmp');
dot = max(find(fname == '.'));
suffix = fname(dot+1:dot+3);
if strcmp(suffix,'pgm') | strcmp(suffix,'raw')
e = rawread(fname);
else e = imread(fname); end
if isrgb(e), e = rgb2gray(e); end
if isa(e,'double')~= 1,
e = double(e);
end
maxe = max(max(e)');
mine = min(min(e)');
img = (e - mine)/(maxe - mine);
Img = e;
[ysize xsize] = size(Img);
fpos = get(HDmainf,'Position');
fw = fpos(3); fh = fpos(4);
k = 2;
xpos = 0.5*(fw - xsize*k);
ypos = 0.5*(fh - ysize*k);
HDorigPic1=subplot(1,1,1); imshow(img);
set(HDorigPic1,'Units', 'pixels','Position',[xpos ypos ysize*k xsize*k],'Units','normal');
title('Original image');
% after read input image, we should smooth it using a gaussian filter,
% also make preparation for LBF.
sigma = 0.00000001; % scale parameter in Gaussian kernel
gauss = fspecial('gaussian',round(3*sigma)*2+1,sigma); % Gaussian kernel
gaussI = convbyfft(Img,gauss);
% set the parameters for Chan-Vese model and others
mu = 1.0; % without initialization regularization
nu = 1; % length weight, maybe, we should also take the area into account
epsilon = 1;% heaviside function regularization and its derivatives
lambda1 = 0.1; % foreground weight
lambda2 = 0.1; % background weight
dt = 0.1; % timestep
iterno = 1000; % iteration number
inittype = 1; % if 1, automatically initialized using a circle located at the image center,
% if 0, manually locate the initial contour
distorconst = 1;% if 1, phi function is a distance function;
% if 0, function phi is a two constant function;
lbf = 0; % if 0, PC-CV; 1, LBF by chunming Li
sigmalbf = 3.0; % for LBF use
ged = 0; %if 0, only PC-CV; 1, PC-CV incorporated with geodesic length
% if using the geodesic length of the curve, we should calculate the
% smoothed image and its edge map.
if ged == 1,
gedk = 0.01;
sigmaged = 3.0;
gaussged = fspecial('gaussian',round(3*sigmaged)*2+1,sigmaged); % Gaussian kernel
gedI = convbyfft(BoundMirrorExpand(Img),gaussged);
[gedx,gedy] = gradient(gedI);
magGED = (gedx.^2.0 + gedy.^2.0);
% normalize this magnitude ...
fmin = min( magGED(:) );
fmax = max( magGED(:) );
magGED = ( magGED - fmin )/( fmax - fmin );
gged = 1.0./(1+ (magGED./gedk) );
figure;imshow(magGED);title('Edge map');
figure;imshow(gged);title('Edge indicator')
[gx,gy]=gradient(gged);
end
% initialize the level set function phi
c0 = 8;
if inittype == 1,
r = 10;%0.15*ysize;
xc = 0.3*xsize;
yc = 0.5*ysize;
[X,Y] = meshgrid(1:xsize, 1:ysize);
phi0 = sqrt((X-xc).^2+(Y-yc).^2) - r;
text(1,2,'Automatically initialized...','Color', 'r');
else
text(1,2,'Left click to start, right click to end','Color', 'r');
phi0 = roipoly - 0.5;
phi0 = double((phi0 < 0.0).*(bwdist(phi0 > 0.0)-0.5) - (phi0 > 0.0).*(bwdist(phi0 < 0.0)-0.5));
end
if ~distorconst,
nn = find( phi0 <0.0);
phi0(nn) = -c0;
nn = find( phi0 > 0.0);
phi0(nn) = c0;
end
hold on; contour(phi0,[0 0],'r'); hold off; axis equal;
pause(0.5);
% next, we will evolve the level set function
gaussI = Img;
tic;
gaussI = BoundMirrorExpand(gaussI);
phi = BoundMirrorExpand(phi0);
for k = 1:iterno,
phi = BoundMirrorEnsure(phi);
dirh = dirach(phi,epsilon);
kappa = curvature(phi);
ls_regul = noreinit(phi,kappa);
if ~lbf, %C-V
[c1,c2] = getcc(gaussI,phi,epsilon);
if ~ged,% just C-V
if mod(k,20)==1,
fprintf(1,'using only piecewise constant model...\n');
end
phi = phi + dt*mu*ls_regul + dt*dirh.*( nu*( kappa ) - lambda1*(gaussI - c1).^2 + lambda2*(gaussI - c2).^2);
else% PC-CV incorporated with geodesic length
if mod(k,20)==1,
fprintf(1,'using piecewise constant model plus geodesic length ...\n');
end
% using simple scheme for geodesic
[px,py]=gradient(phi);
normDP=sqrt(px.^2 + py.^2 + 1e-10);
Nx=px./normDP;
Ny=py./normDP;
geodesicKappa = gged.*kappa + (gx.*Nx + gy.*Ny) ;
phi = phi + dt*mu*ls_regul + dt*dirh.*( nu*( geodesicKappa ) - lambda1*(gaussI - c1).^2 + lambda2*(gaussI - c2).^2);
end
else% LBF
if mod(k,20)==1,
fprintf(1,'using Local Binary Fitting...\n');
end
[e1,e2,f] = getee(gaussI,phi,epsilon,sigmalbf);
phi = phi + dt*mu*ls_regul + dt*dirh.*( nu*kappa - lambda1*e1 + lambda2*e2 );
end
if ~distorconst,
nn = find(phi > c0);
phi(nn) = c0;
nn = find(phi < -c0);
phi(nn) = -c0;
end
if mod(k,10) == 0,
pause(0.5); imshow(img);
hold on; contour(BoundMirrorShrink(phi),[0 0], 'r');
title(['iteration no. ' num2str(k)]);hold off
end
end
phi=BoundMirrorShrink(phi);
fprintf(1, 'total excution time is [%s]\n', num2str(toc));
imshow(img);
hold on; contour(phi,[0 0
