function [rhos,analysis,ys]=vision_engine_b(ID,Itemps,templs)
%executes inference engine in Strategies for Seeing, Experiment 2
%where generators are arcs
%ID=oberved image
%Itemps =matrix of templates extended periodically to length tmax
%templs =vector of template lengths
%to run load workspace testb

%1.ENTER.
%make image binary by thresholding
I=ID>200;
%regularize image
I=bwmorph(I,'close');
I=bwmorph(I,'fill');
I=bwmorph(I,'spur');
I=bwmorph(I,'clean');
I=double(I);

%opening frames
page1
pause(1)
close
page2
pause(1)
close
page3b
pause(1)
close


%2. SACCADIC SEARCH FOR FOVEAS
%find the ncomp topological components and their sizes
labels=bwlabel(I);
ncomp=max(max(labels));
attention=zeros(1,ncomp);
for i=1:ncomp
   attention(i)=sum(sum(labels==i));
end 
'number of topological components in observed image'
ncomp


compnumber=1;
[maxatt,i]=max(attention);
while maxatt>100
   fovea=labels==i;
  
   %find boundary B for fovea and chaincode BC
   B=bwperim(fovea);
   %regularize boundary
   B=bwmorph(B,'spur');
   B=bwmorph(B,'spur');
B=bwmorph(B,'spur');
B=bwmorph(B,'spur');

 B=double(B);
 
 %display observed image and current fovea
 image(.1.*ID)
axis('image')
title('OBSERVED IMAGE')
pause
close
image(30.*B)
axis('image')
title('CURRENT FOVEA')
pause
close



 %3. PREPARE 
 %find starting pixel
 
 v=any(B);
 [Y,colindex]=max(v);
 [Z,rowindex]=max(B(:,colindex));
 start=[rowindex,colindex];
 %introduce window W of 8-neighborhood
 W=[2 3 4;1 0 5; 8 7 6];
 %and vector of offsets
 off=[-1 0 1];
 %also 2-column matrix with x and y increments for given code
 incr=[0 -1;-1 -1;-1 0;-1 1;0 1;1 1;1 0; 1 -1];
 
 %4. INITIALIZE
 %initializes two column coordinate matrix
 xs=[];
 %put frame of 0's around B
 [l1,l2]=size(B);
 Bframe=zeros(l1+2,l2+2);
 Bframe(2:l1+1,2:l2+1)=B;

 %remember frame!
 start=start+[1 1];
 Bframe(start(1),start(2))=0;

 %initialize for chain coding
 Neigh=Bframe(start(1)+off,start(2)+off);
 chaincode=Neigh.*W;
 code=max(max(chaincode));
 BC=code;
 next=start+incr(code,:);
 Bframe(next(1),next(2))=0;
 xs=[xs;next];
 
 Neigh=Bframe(next(1)+off,next(2)+off);
 chaincode=Neigh.*W;
 code=max(max(chaincode));
 BC=[BC;code];
 Bframe(next(1),next(2))=0;
 next=next+incr(code,:);
 xs=[xs;next];
 
 Neigh=Bframe(next(1)+off,next(2)+off);
 chaincode=Neigh.*W;
 code=max(max(chaincode));
 BC=[BC;code];
 Bframe(next(1),next(2))=0;
 next=next+incr(code,:);
 xs=[xs;next];
 
 'chaincoding boundary...'
 
 %5. START CHAINCODING LOOP
 t=3;
 %start iterating along boundary
 %use sup norm distance for termination criterion
 

 while max(abs(next-start))>1
   Neigh=Bframe(next(1)+off,next(2)+off);
      chaincode=Neigh.*W;
      code=max(max(chaincode));
      BC=[BC;code];
            Bframe(next(1),next(2))=0;
      next=next+incr(code,:);
      xs=[xs;next];
      t=t+1;
 end      
   Bframe(start(1),start(2))=1;
   Neigh=Bframe(next(1)+off,next(2)+off);
      chaincode=Neigh.*W;
      code=max(max(chaincode));
      BC=[BC;code];
      
 
      
%6. COMPUTE CURVATURE     
      %compute curvature at equidistant points of ID
      %separated by distance dl
      %computes curvature function rhos
%for coordinate matrix ys. 
l=length(xs(:,1));
%decimate boundary data by a factor d
d=5;
m=floor(l/d);
%extend perimeter by one period
xs=[xs;xs(1,:)];
xs=[xs;xs];
xs=xs(d.*[1:m],:);

%6. REGULARIZE BOUNDARY
%regularize discrete boundary by cyclic convolution
l=length(xs(:,1));
a=.2;p=floor(l/4);v=a.^[0:p-1];w=a.^[p:-1:1];
u=[v,zeros(1,l-2*p),w];u=u./sum(u);
uhat=fft(u);x1shat=fft(xs(:,1)');x2shat=fft(xs(:,2)');
%execute convolution by Fourier transform
xs(:,1)=real(ifft(uhat.*x1shat))';xs(:,2)=real(ifft(uhat.*x2shat))';
%extend boundary data one period
xs=[xs;xs];

'computing curvature...'

%7. GET ARCLENGTH, CURVATURES
%compute arclengths
lengths=zeros(1,m);
for t=1:m
   lengths(t)=norm(xs(t+1,:)-xs(t,:));
end

cuml=cumsum(lengths);
%now get discretized curvatures rhos
dphis=zeros(1,m);rhos=zeros(1,m);

for t=1:m
 
         a=xs(mod(t+3,m)+1,:)-xs(t,:);a=a./norm(a);
         b=xs(t,:)-xs(mod(t-5,m)+1,:);b=b./norm(b);
         %change to image (matrix) coordinates
         a=[a(2), -a(1)];b=[b(2),-b(1)];
  
      
      x=a(1)*b(1)+a(2)*b(2);
      y=a(1)*b(2)-a(2)*b(1);
      r=sqrt(x^2+y^2);
      %dphis(t)=asin(y/r);
      dphis(t)=(x>=0)*asin(y/r)+(x<0)*(pi-asin(y/r));
      rhos(t)=2*dphis(t)/(norm(xs(mod(t+3,m)+1,:)-xs(t,:))+norm(xs(t,:)-xs(mod(t-5,m)+1,:)));
   
   end
   %arc element length=dl
   
   dl=5;
   n=floor(cuml(m)/dl);
   %interpolate curvature function to equi-distant points
   rhos=spline([0,cuml(1:m-1)],rhos,dl.*[0:n-1]); 
   ys=zeros(n,2);
   y1=spline([0,cuml(1:m-1)],xs(1:m,1)',dl.*[0:n-1]); 
   y2=spline([0,cuml(1:m-1)],xs(1:m,2)',dl.*[0:n-1]); 
   ys(:,1)=y1';ys(:,2)=y2';

'dynamic programming, backwards...'
pause(.1)
'Please wait a while...'


  %8. INITIALIZE FOR DYNAMIC PROGRAMMING     
      %prepare for  pattern analysis by dynamic programming
      % for state variables given as 
%2-vectors (alpha,t); alpha for arctype and t for arclength
%alpha=1,2...nalpha;t=1,2,..tmax. %jumpc=cost for jumping from arctype to arctype

jumpc=.2;
%elastc=cost of elastic change with fixed arctype
elastc=.005;
[nalpha,tmax]=size(Itemps);
alphas=1:nalpha;ts=1:tmax;
M=zeros(nalpha,tmax);

%result is 2-column matrix analysis,first column alpha values,
%second columnt t values
ID=rhos;
l=length(ID);
analysis=zeros(l,2);
u=ones(1,nalpha);v=ones(1,tmax);
%cost starting at t with (alphavalue,tvalue)
cost=zeros(l,nalpha,tmax);

%initialize dynamic programming (backwards)
for alpha =1:nalpha
   for t=1:tmax
      diffs=mod(t-1-ts,templs(alpha));
     
      M=(ID(l)-Itemps(alpha,t))^2+jumpc.*kron((alphas~=alpha)',v)+...
         elastc.*kron((alphas==alpha)',diffs.*(templs(alpha)-diffs));
      [x,m]=min(M);
      [y,n]=min(x);
      cost(l,alpha,t)=y;
   end
end

%9. BACKWARD LOOP OF DYNAMIC PROGRAMMING
%start iterating dynamic programming,backward loop,t1 arc length parameter
for t1=l-1:-1:1
   for alpha =1:nalpha
      for t=1:tmax
         
                  diffs=mod(t-1-ts,templs(alpha));
                  cost1=reshape(cost(t1+1,:,:),[nalpha tmax]);
      	M=cost1+(ID(t1)-Itemps(alpha,t))^2+jumpc.*kron((alphas~=alpha)',v)+...
        	 elastc.*kron((alphas==alpha)',diffs.*(templs(alpha)-diffs));
      	[x,m]=min(M);
      	[y,n]=min(x);
      	cost(t1,alpha,t)=y;
   	end
   end
   ['arc element no. ',num2str(t1)]
end
c1=reshape(cost(1,:,:),[nalpha tmax]);
[U,m]=min(cost1);
[V,n]=min(U);
analysis(1,1)=m(n);analysis(1,2)=n;
alpha=analysis(1,1);t=analysis(1,2);
cost(1,m(n),n)

%10. FORWARD LOOP OF DYNAMIC PROGRAMMING
'dynamic programming, forward'
pause(.1)

%now iterate forwards
for t1=2:l
           
                  diffs=mod(n+1-ts,templs(alpha));
         cost1=reshape(cost(t1-1,:,:),[nalpha tmax]);

      	M=cost1+(ID(t1-1)-Itemps(alpha,t))^2+jumpc.*kron((alphas~=alpha)',v)+...
            elastc.*kron((alphas==alpha)',diffs.*(templs(alpha)-diffs));
                	[x,m]=min(M);
         [y,n]=min(x);
     analysis(t1,1)=m(n);alpha=m(n);
   analysis(t1,2)=n;t=n;

end
    
%11. DISPLAY RESULT OF PATTERN ANALYSIS 
   %save and print current analysis on screen
   %eval(['analysis',num2str(compnumber),'= analysis'])
   
   %display analysis
   'COLOR CODE:'
   'red=1'
   'green=2'
   'blue=3'
   'yellow=4'  
   'magenta=5'

   color_b(ys,analysis)
   title('MAP estimates of arcelements')
      axis('image')
   pause
   close
   
   %display deformation along arcs
   plot(analysis(:,2))
   axis('image')
   title('Estimated deformation of arcs')
   
   pause
   close
   
%12. GO TO NEXT FOVEA IF ITS ATTENTION CALLS FOR IT   
   attention(i)=0;
   [maxatt,i]=max(attention);
   compnumber=compnumber+1;
end 

%Auxiliary functions

function page1()
load page1                             

a = figure('Color',[0.8 0.8 0.8], ...
	'Colormap',mat0, ...
	'PointerShapeCData',mat1, ...
	'Position',[120 120 560 420]);
b = axes('Parent',a, ...
	'Box','on', ...
	'CameraUpVector',[0 1 0], ...
	'CameraUpVectorMode','manual', ...
	'Color',[1 1 1], ...
	'ColorOrder',mat2, ...
	'Layer','top', ...
	'Visible','off', ...
	'XColor',[0 0 0], ...
	'XLim',[0.5 800.5], ...
	'XLimMode','manual', ...
	'YColor',[0 0 0], ...
	'YDir','reverse', ...
	'YLim',[0.5 600.5], ...
	'YLimMode','manual', ...
	'ZColor',[0 0 0]);
c = image('Parent',b, ...
	'CData',mat3, ...
	'XData',[1 800], ...
	'YData',[1 600]);
c = text('Parent',b, ...
	'Color',[0 0 0], ...
	'HandleVisibility','callback', ...
	'HorizontalAlignment','center', ...
	'Position',[399.576 -11.8167 0], ...
	'VerticalAlignment','bottom');
set(get(c,'Parent'),'Title',c);

function page2()
load page2                             

a = figure('Color',[0.8 0.8 0.8], ...
	'Colormap',mat0, ...
	'PointerShapeCData',mat1, ...
	'Position',[120 120 560 420]);
b = axes('Parent',a, ...
	'Box','on', ...
	'CameraUpVector',[0 1 0], ...
	'Color',[1 1 1], ...
	'ColorOrder',mat2, ...
	'Layer','top', ...
	'Visible','off', ...
	'XColor',[0 0 0], ...
	'XLim',[0.5 800.5], ...
	'XLimMode','manual', ...
	'YColor',[0 0 0], ...
	'YDir','reverse', ...
	'YLim',[0.5 600.5], ...
	'YLimMode','manual', ...
	'ZColor',[0 0 0]);
c = image('Parent',b, ...
	'CData',mat3, ...
	'XData',[1 800], ...
	'YData',[1 600]);
c = text('Parent',b, ...
	'Color',[0 0 0], ...
	'HandleVisibility','callback', ...
	'HorizontalAlignment','center', ...
	'Position',[399.576 -11.8167 0], ...
	'VerticalAlignment','bottom');
set(get(c,'Parent'),'Title',c);


function color_b(ys,analysis)
%displays result of vision_engine_b with colorcoded arc elements
colors='rgbym';
n=length(ys(:,1));ys=[ys;ys(1,:)];
for i=1:n
   v=[i i+1];
   plot(ys(v,1)',ys(v,2)',colors(analysis(i,1)'))
   hold on
end
hold off


function fig = page3b()

load page3b

h0 = figure('Color',[0.8 0.8 0.8], ...
	'Colormap',mat0, ...
	'PointerShapeCData',mat1, ...
	'Position',[197 121 560 406], ...
	'Tag','Fig1');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[0.7607843137254901 0.749019607843137 0.647058823529412], ...
	'FontSize',16, ...
	'ListboxTop',0, ...
	'Position',[29.25 189.75 338.25 48.75], ...
	'String','STRATEGY WHEN GENERATORS ARE ARCS', ...
	'Style','text', ...
	'Tag','StaticText1');
if nargout > 0, fig = h0; end