function vision_engine_a(Itemps,ID,ntries)
 
%executes  STRATEGIES FOR SEEING for Experiment 1 when
%generators are areas (filled regions in the plane.
%Complexity bounded by ntries as max. configuration size
%Binary matrices representing BW pictures as planes in 3D array Itemps
%where last index enumerates the images
%ID is single binary matrix that is being observed
%parameters chosen for images of approximately size 100 x 100

%opening frames
page1
pause(2)
close
page2
pause(2)
close

%1.ENTER
%display observed image
image(30.*ID)
axis('image')
title('noisy observed image')
pause
close

%start vision algorithm
[l1,l2,nalpha]=size(Itemps);

%2.ARRAYS
%set up 3D array to contain preliminary hypotheses about the observed image ID
hypotheses=zeros(l1,l2,nalpha);

%setup two 3D arrays to store horizontal and vertical deformations
%for the different hypothetical generator templates
s1s=zeros(l1,l2,nalpha);s2s=s1s;

%likelihood energies under various hypotheses 
qs=zeros(1,nalpha);
%prior energies under various hypotheses
ps=qs;
%posterior energies under various hypotheses
rs=ps;

%3.FOURIER
%use 2*r^2 Fourier coefficients
r=5;
% for use in local group deformation compute basis function
v1=[1:l1];v2=[1:l2];w=[1:r];
M1=sin(pi/l1.*kron(v1,w'));M2=sin(pi/l2.*kron(v2,w'));
%and the two matrices
w=w.^2;
N1=kron(ones(1,r),w');N2=kron(w,ones(1,r)');

%4.PREPARE
%prepare for diffeomorphism by filtering images with Gaussian
%filter G of size 16x16 and bandwidth 1.5
G=zeros(l1,l2);
v1=[42:57];v2=[52:67];
G(v1,v2)=fspecial('Gaussian',16,1.5);
%which leads to modified image IDm and template  array Itempsm
IDm=filter2(G,ID);
for k=1:nalpha
   Itempsm(:,:,k)=filter2(G,Itemps(:,:,k));
end

%5.INITIALIZE
%set so far unexplained unexplained part of modified image to be IDmu
IDmu=IDm;
%set initial attention high
maxatt=1000;

%start with lowest complexity: configuration size 1
attempt=1;
while (attempt <= ntries)&(maxatt > 30)


%6.START
%start with first template
alpha=1;

while (alpha <= nalpha)
   
   IDm=IDmu;
   
%set template to be the one associated with hypothesis (generator) alpha
Itempm=Itempsm(:,:,alpha);
%get attention for IDmu
A=filter2(ones(15),IDmu,'same');
[X2,X1]=meshgrid(1:l2,(1:l1)');
[Y,k]=max(A);
[Z,l]=max(Y);k=k(l);[k l]

'max attention is'
maxatt=Z
'centered at'
mid=[k+7,l+7];

%get attention for Itempm
A=filter2(ones(15),Itempm,'same');
[X2,X1]=meshgrid(1:l2,(1:l1)');
[Y,k]=max(A);
[Z,l]=max(Y);k=k(l);[k l]
'centered at'
mids=[k+7,l+7];

%7.TRANSLATION
%Start saccadic translation group
Itempnew=fastint(Itempm,mids(1)-mid(1),mids(2)-mid(2));
%Itempnew is translated Itempm

clf
subplot(2,2,1),image(30.*ID)
axis('image')
title('OBSERVED IMAGE')
subplot(2,2,2),image(Itempnew.*60)
axis('image')
title(' translated image Itempnew')
subplot(2,2,3),image(IDm.*60)
axis('image')
title('MODIFIED OBSERVED IMAGE IDm')
subplot(2,2,4),image(Itempm.*60)
axis('image')
title('MODIFIED TEMPLATE')

pause(1) 

%8.ROTATION
%start saccadic rotation group using rotation covariant
%center of template
Mt1=sum(sum(Itempnew.*X1))/sum(sum(Itempnew));
Mt2=sum(sum(Itempnew.*X2))/sum(sum(Itempnew));
Mt1=round(Mt1);Mt2=round(Mt2);

v1=[max(1,Mt1-10 ):min(l1,Mt1+10)];
n1=length(v1);
v2=[max(1,Mt2-10):min(l2,Mt2+10)];
n2=length(v2);
%introduce fovea sets for area around centers of attention
Itempfovea=Itempnew(v1,v2);
IDfovea=IDm(v1,v2);

subplot(2,2,1),image(IDm.*60)
axis('image')
title('MODIFIED OBSEVED IMAGE')
fovea=zeros(l1,l2);fovea(v1,v2)=30.*ones(n1,n2);
subplot(2,2,2);image(fovea)
axis('image')
title('FOVEA')
subplot(2,2,3),image(60.*IDm(v1,v2))
axis('image')
title('FOVEA PART OF MODIFIED OBSERVED IMAGE')
pause(1)
xs=1:n1;ys=1:n2;xs=kron(xs',ones(1,n2));ys=kron(ones(n1,1),ys);

%compute centers for foveas
st=sum(sum(Itempfovea));
sd=sum(sum(IDfovea));
mt1=sum(sum(xs.*Itempfovea))/st;
mt2=sum(sum(ys.*Itempfovea))/st;
md1=sum(sum(xs.*IDfovea))/sd;
md2=sum(sum(ys.*IDfovea))/sd;

X=kron(ones(1,l2),[1:l1]');
Y=kron([1:l2],ones(l1,1));

%translate further
U1x=xs-mt1;U1y=ys-mt2;

%compute moments of inertia for fovea part of
%of modified template

R1=zeros(2);R2=R1;
R1(1,1)=sum(sum(Itempfovea.*U1x.^2))/st;
R1(1,2)=sum(sum(Itempfovea.*U1x.*U1y))/st;
R1(2,1)=R1(1,2);
R1(2,2)=sum(sum(Itempfovea.*U1y.^2))/st;

%and for modified observed image

U2x=xs-md1;U2y=ys-md2;

R2(1,1)=sum(sum(IDfovea.*U2x.^2))/sd;
R2(1,2)=sum(sum(IDfovea.*U2x.*U2y))/sd;
R2(2,1)=R2(1,2);
R2(2,2)=sum(sum(IDfovea.*U2y.^2))/sd;

%get eigen elements of moment matrices
[O1,v1]=eig(R1);
[O2,v2]=eig(R2);

[v1,i1]=sort(diag(v1));
[v2,i2]=sort(diag(v2));
O1=O1(:,i1);
O2=O2(:,i2);

%make determinant =+1 for SO(2)
D=det(O1)>0;O1(:,1)=D.*O1(:,1)-(1-D).*O1(:,1);
D=det(O2)>0;O2(:,1)=D.*O2(:,1)-(1-D).*O2(:,1);
O=O1'*O2

%write in terms of cosine
c=O(1,1);
%and sine
s=O(2,1);

U1=X1-Mt1;U2=X2-Mt2;
Z1=Mt1+c.*U1+s.*U2;Z2=Mt2-s.*U1+c.*U2;
Z1=Z1-X1;Z2=Z2-X2;

%make sure that boundaries of background space haveintensity 0
Itempnew(1,:)=zeros(1,l2);Itempnew(l1,:)=zeros(1,l2);
Itempnew(:,1)=zeros(l1,1);Itempnew(:,l2)=zeros(l1,1);

Itempnew3=fastint(Itempnew,Z1,Z2);
%Itempnew's are propagated modified templates
%Istar images are propagated templates (binary)

%9. DISAMBIGUATE
%BECAUSE OF THE AMBIHUITY OF THE SO(2)-INVARIANT TURN THE ESTIMATED
%TEMPLATE 180 DEGREES
c=-c;s=-s;
U1=X1-Mt1;U2=X2-Mt2;
Z1=Mt1+c.*U1+s.*U2;Z2=Mt2-s.*U1+c.*U2;
Z1=Z1-X1;Z2=Z2-X2;

	Itempnew4=fastint(Itempnew,Z1,Z2);
   %display modified observed image and new rotated ,modified templates
   clf
subplot(2,2,1),image(IDm.*60)
title('OBSERVED (MODIFIED) IMAGE')
subplot(2,2,2),image(Itempnew3.*60)
title('NEW, ROTATED (MODIFIED) TEMPLATE')
subplot(2,2,3),image(Itempnew4.*60)
title('180 DEGREES ROTATION')
pause (1)
close

% select most likely of the two rotations

qbefore3=sum(sum((Itempnew3-IDm).^2));
qbefore4=sum(sum((Itempnew4-IDm).^2));
Itempnew=(qbefore3<qbefore4).*Itempnew3 +(qbefore3 >=qbefore4).*Itempnew4;
%Istar=(qbefore3<qbefore4).*Istar1 +(qbefore3 >=qbefore4).*Istar2;

image(Itempnew.*60)
title('SELECTED ROTATION')
pause(1)

%10.DIFF
%Start local group transformations of diffeomorphisms
[G2,G1]=gradient(Itempnew);
%Variability coefficient sigma for prior probability measure
sigma=2;b=1/sigma^2;
%and variances of random Fourier coefficients associated with
%low divergence 

%noise variance
tau=.005; c=1/tau^2;
%stepsize and number of iterations
dt=.000001;
niter=10;

%initialize Fourier coefficients
t1=zeros(r);t2=zeros(r);

%initialize local group elements to identity

s1=zeros(l1,l2);s2=s1;

gradl1=zeros(r);gradl2=gradl1;

%11.GRADIENT
%Start crude gradient search
'the displayed values rtemp are proportional to posterior energy'
pause(.3)
for t=1:niter
   
%introduce deformed modified template   
intItempnew=fastint(Itempnew,s1,s2);
%display propagated template together with the first modified one
subplot(2,2,1),image(30.*intItempnew)
axis('image')
title('propagated template')
subplot(2,2,2),image(30.*Itempnew)
axis('image')
title('first modified template')
subplot(2,2,3),image(30.*abs(intItempnew-IDm))
axis('image')
title('DIFFERENCE |prop.temp.-mod. ID|')
subplot(2,2,4),text(.02,.2,'LOCAL GROUP ITERATIONS')
pause(.8)
axis off

      intG1=fastint(G1,s1,s2);intG2=fastint(G2,s1,s2);
      s1=M1'*t1*M2;s2=M1'*t2*M2;
      
%compute prior energy gradient    
gradp1=t1.*N1.*l1;gradp2=t2.*N2.*l2;

%compute likelihood energy gradient
   for mu=1:r
       for nu=1:r
      
          gradl1(mu,nu)=sum(sum((intItempnew-IDm).*intG1.*(M1(mu,:)'*M2(nu,:))));
          gradl2(mu,nu)=sum(sum((intItempnew-IDm).*intG2.*(M1(mu,:)'*M2(nu,:))));

       end
    end
    
%and total gradient    
   grad1=b.*gradp1+c.*gradl1;
   grad2=b.*gradp2+c.*gradl2; 
   
%update Fourier coefficients
	t1=t1-grad1.*dt;
	t2=t2-grad2.*dt;
%and get deformation fields   
   s1=M1'*t1*M2;s2=M1'*t2*M2;
%get prior energy
 qafter=sum(sum((IDm-intItempnew).^2));
 T1=(t1.*N1).^2;T2=(t2.*N2).^2;
 ptemp=b*l1*sum(sum(T1))+b*l2*sum(sum(T2));
 %and add posterior energy to get posterior energy
 'posterior energy ='
 rtemp=ptemp+c*qafter
end

%deformation fields
s1s(:,:,alpha)=s1;s2s(:,:,alpha)=s2;

%display (decimated) deformation field?
u1=s1;u2=s2;
dec1=[1:ceil(l1/5)].*5;dec2=[1:ceil(l2/5)].*5;
s1=s1(dec1,dec2);s2=s2(dec1,dec2);
save l s1 s2

def(s1,s2);
pause

%likelihood energy
s1=u1;s2=u2;
qs(alpha)=c*qafter;
%prior energy
T1=(t1.*N1).^2;T2=(t2.*N2).^2;
 ps(alpha)=b*l1*sum(sum(T1))+b*l2*sum(sum(T2));

  %posterior energy
 rs(alpha)=ps(alpha)+qs(alpha)
 
 %enter deformed templates for hypothesis No. alpha
 hypotheses(:,:,alpha)=intItempnew;
 

%go to next template number if less than nalpha
alpha = alpha+1;
end

%12.CREDIBLE
%Decide on most credible template hypothesis alphastar
[Y,alphastar]=min(rs);
Imstar=hypotheses(:,:,alphastar);
'maximum belief in hypothesis No.='
alphastar
ps
qs
rs

	
subplot(1,2,1),image(60.*IDm)
title('MODIFIED ID')
axis('image')
subplot(1,2,2),image(60.*Imstar)
title(' RESTORED MODIFIED IMAGE')
axis('image')
pause(1)
close

image(zeros(l1,l2))
title('Recognized object type')
text(l1/4,l2/2,num2str(alphastar),'FontSize',128)
pause(1)
close


%13.REST
%Get unexplained part IDmu of IDm
IDmu=abs(IDmu-Imstar);
subplot(1,2,1),image(60.*IDm)
title('MODIFIED ID')
axis('image')
subplot(1,2,2),image(60.*IDmu)
title('UNEXPLAINED PART OF MODIFIED ID')
axis('image')
pause(1)

%get attention for IDmu
A=filter2(ones(15),IDmu,'same');
[X2,X1]=meshgrid(1:l2,(1:l1)');
[Y,k]=max(A);
[Z,l]=max(Y);k=k(l);[k l]

'max is'
maxatt=Z


%increse complexity of configuration if less than ntry
attempt=attempt+1;
end

if sum(sum(IDmu))>.8*sum(sum(IDm))
   'image not understood'
else
end

%14.STOP

%auxiliary functions:
function IDf
%displays image ensemble as requested by the user
load IDf                               

a = figure('Color',[0.8 0.8 0.8], ...
	'Colormap',mat0, ...
	'PointerShapeCData',mat1, ...
	'Position',[397 101 560 420], ...
	'Tag','Fig1');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'BackgroundColor',[0.752941 0.752941 0.752941], ...
	'Position',[110.25 240 114.75 45.75], ...
	'String','''Would you like to see some ID''s?''', ...
	'Style','text', ...
	'Tag','StaticText2');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','load K;number =ceil(rand(1)*16);I=IDmany(:,:,number);image(30.*I)', ...
	'Position',[75 197.25 45 15], ...
	'String','Yes', ...
	'Tag','Pushbutton1');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','delete(gcf)', ...
	'Position',[193.5 197.25 45 15], ...
	'String','No', ...
	'Tag','Pushbutton2');
b = axes('Parent',a, ...
	'Box','on', ...
	'CameraUpVector',[0 1 0], ...
	'CameraUpVectorMode','manual', ...
	'Color',[1 1 1], ...
	'ColorOrder',mat2, ...
	'Layer','top', ...
	'Position',[0.125 0.111905 0.778571 0.488095], ...
	'Tag','Axes1', ...
	'XColor',[0 0 0], ...
	'XLim',[0.5 1.5], ...
	'XLimMode','manual', ...
	'YColor',[0 0 0], ...
	'YDir','reverse', ...
	'YLim',[0.5 1.5], ...
	'YLimMode','manual', ...
	'ZColor',[0 0 0]);
c = image('Parent',b, ...
	'CData',300, ...
	'Tag','Axes1Image1', ...
	'XData',1, ...
	'YData',1);
c = text('Parent',b, ...
	'Color',[0 0 0], ...
	'HandleVisibility','callback', ...
	'HorizontalAlignment','center', ...
	'Position',[0.997696 1.61823 0], ...
	'Tag','Axes1Text4', ...
	'VerticalAlignment','cap');
set(get(c,'Parent'),'XLabel',c);
c = text('Parent',b, ...
	'Color',[0 0 0], ...
	'HandleVisibility','callback', ...
	'HorizontalAlignment','center', ...
	'Position',[0.43318 1.00739 0], ...
	'Rotation',90, ...
	'Tag','Axes1Text3', ...
	'VerticalAlignment','baseline');
set(get(c,'Parent'),'YLabel',c);
c = text('Parent',b, ...
	'Color',[0 0 0], ...
	'HandleVisibility','callback', ...
	'HorizontalAlignment','right', ...
	'Position',[0.336406 -0.32266 0], ...
	'Tag','Axes1Text2', ...
	'Visible','off');
set(get(c,'Parent'),'ZLabel',c);
c = text('Parent',b, ...
	'Color',[0 0 0], ...
	'HandleVisibility','callback', ...
	'HorizontalAlignment','center', ...
	'Position',[0.997696 0.465517 0], ...
	'Tag','Axes1Text1', ...
	'VerticalAlignment','bottom');
set(get(c,'Parent'),'Title',c);

function noi()
%interrogates the user about noise level
load noi                               

a = figure('Color',[0.8 0.8 0.8], ...
	'Colormap',mat0, ...
	'PointerShapeCData',mat1, ...
	'Position',[227 116 560 420], ...
	'Tag','Fig1');
b = uicontrol('Parent',a, ...
	'BackgroundColor',[0.752941 0.752941 0.752941], ...
	'Callback','val=get(gco,''Value'');save V val', ...
	'Max',300, ...
	'Position',[21 20 20 320], ...
	'Style','slider', ...
	'Tag','Slider1', ...
	'Value',291);
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'BackgroundColor',[0.752941 0.752941 0.752941], ...
	'Position',[202.5 128.25 133.5 93], ...
	'String','Enter noise level,expected number of pixel errors', ...
	'Style','text', ...
	'Tag','StaticText1');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'BackgroundColor',[0.752941 0.752941 0.752941], ...
	'Position',[42.75 225 46.5 20.25], ...
	'String','300', ...
	'Style','text', ...
	'Tag','StaticText2');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'BackgroundColor',[0.752941 0.752941 0.752941], ...
	'Position',[42 6.75 45 25.5], ...
	'String','0', ...
	'Style','text', ...
	'Tag','StaticText3');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'BackgroundColor',[0.752941 0.752941 0.752941], ...
	'Position',[43.5 75 47.25 25.5], ...
	'String','100', ...
	'Style','text', ...
	'Tag','StaticText5');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'BackgroundColor',[0.752941 0.752941 0.752941], ...
	'Position',[43.5 152.25 48 32.25], ...
	'String','200', ...
	'Style','text', ...
	'Tag','StaticText6');

function JD=addbin(J,errorprob)
%computes deformed image JD by adding binary noise to 
%BW image J; addition modulo 2
[l1,l2]=size(J);
errors=rand(l1,l2)<errorprob;
JD=rem(J+errors,2);

function defgridi(t1,t2)
%displays deformed grid with displacements 
%given by matrices t1 and t2
%in matrix (image) coordinates
t1=t1./5;t2=t2./5;
%above because of decimation with factor 5
[l1,l2]=size(t1);xs=1:l1;ys=1:l2;
[Y,X]=meshgrid(ys,xs);
X=X+t1;Y=Y+t2;
Y1=Y;Y=l1-X+1;X=Y1;
clf
hold on
for x=1:l1-1
	for y=1:l2-1
		plot([X(x,y),X(x+1,y)],[Y(x,y),Y(x+1,y)]);
	end
end		
for x=1:l1-1
	for y=1:l2-1
		plot([X(x,y),X(x,y+1)],[Y(x,y),Y(x,y+1)]);	
	end
end

function def(s1,s2)
%pushbuttons for display of deformation field
load def                               

a = figure('Color',[0.8 0.8 0.8], ...
	'Colormap',mat0, ...
	'PointerShapeCData',mat1, ...
	'Position',[341 -65 560 420], ...
	'Tag','Fig1');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'BackgroundColor',[0.752941 0.752941 0.752941], ...
	'Position',[120 223.5 132.75 61.5], ...
	'String',mat2, ...
	'Style','text', ...
	'Tag','StaticText1');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','load l;defgridi(s1,s2);pause;close', ...
	'Position',[121.5 172.5 131.25 15.75], ...
	'String','Display', ...
	'Tag','Pushbutton1');
b = uicontrol('Parent',a, ...
	'Units','points', ...
	'Callback','close', ...
	'Position',[121.5 138 132.75 18], ...
	'String','Do not display', ...
	'Tag','Pushbutton2');

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);