% TEST_CCLENGTH   Demo for perimeter estimation based on chain codes
%
% (c)2010, Cris Luengo, Centre for Image Analysis, Uppsala, Sweden.

function test_cclength

R = 80;
N = 16;
phi = linspace(0,pi/2,41);
perim = zeros(length(phi),4);
for ii=1:length(phi)
   for jj=1:N
      offset = rand(2,1);
      x = (xx+offset(1))*cos(phi(ii)) - (yy+offset(2))*sin(phi(ii));
      y = (xx+offset(1))*sin(phi(ii)) + (yy+offset(2))*cos(phi(ii));
      img = abs(x)<R & abs(y)<R;
      cc = dip_imagechaincode(+img,2,1);
      perim(ii,1) = perim(ii,1) + cclength(cc.chain,'count');
      perim(ii,2) = perim(ii,2) + cclength(cc.chain,'simple');
      perim(ii,3) = perim(ii,3) + cclength(cc.chain,'weighted');
      perim(ii,4) = perim(ii,4) + cclength(cc.chain,'optimal');
   end
end
perim = perim/N;
figure
plot(phi,perim)
hold on
plot([0,pi/2],8*R*[1,1],'k:')
set(gca,'xlim',[0,pi/2],'xtick',0:pi/8:pi/2,'xticklabel',0:22.5:90)
legend({'pixel count','simple','weighted','corner count'},'Location','SouthEast')
xlabel('angle')
ylabel('estimated perimeter')
title('Estimating line length using chain codes')
print -dpng -r0 cclength1.png


rad = [2,3,5,10,20,30,50,100,200]';
sz = [10,10]+2*rad(end);
trueperim = 2*pi*rad;
err = zeros(length(rad),4);
N = 16;
for ii=1:length(rad)
   for jj=1:N
      img = ( (xx(sz)+rand)^2 + (yy(sz)+rand)^2 ) < rad(ii)^2;
      cc = dip_imagechaincode(+img,2,1);
      err(ii,1) =  err(ii,1) + abs(trueperim(ii)-cclength(cc.chain,'count'));
      err(ii,2) =  err(ii,2) + abs(trueperim(ii)-cclength(cc.chain,'simple'));
      err(ii,3) =  err(ii,3) + abs(trueperim(ii)-cclength(cc.chain,'weighted'));
      err(ii,4) =  err(ii,4) + abs(trueperim(ii)-cclength(cc.chain,'optimal'));
   end
end
err = err*(100/N)./repmat(trueperim,[1,4]);
figure
hold all
plot(rad,err);
set(gca,'xscale','log','yscale','log','xlim',[rad(1)-2,rad(end)+2])
legend({'pixel count','simple','weighted','corner count'},'Location','SouthWest')
xlabel('radius')
ylabel('absolute perimeter error (%)')
title('Estimating circle perimeter using chain codes')
print -dpng -r0 cclength2.png


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function p = cclength(cc,method)

switch method

case 'count'
   % The simplest method, just counts border pixels
   p = length(cc);

case 'simple'
   % Diagonal pixels have a length of sqrt(2)
   even = mod(cc,2)==0;
   p = sum(even) + sum(~even)*sqrt(2);

case 'weighted'
   % Straight and diagonal pixels have weights that
   % are optimal for a line under a random orientation.
   % Proffitt and Rosen, 1979.
   even = mod(cc,2)==0;
   p = sum(even)*0.948 + sum(~even)*1.340;

case {'optimal','best'}
   % Same as 'weighted', but also counting corners.
   % Vossepoel and Smeulders, 1982.
   even = mod(cc,2)==0;
   corner = cc~=[cc(end),cc(1:end-1)];
   p = sum(even)*0.980 + sum(~even)*1.406 - sum(corner)*0.091;
   
end

% Account for the fact we're measuring the length of the chain,
% which joins the pixel centers, and therefore represents a
% slightly smaller object than that we discritized. This
% correction assumes the object is a circle.
p = p + pi;
