% 8th Feb 2007. I have been asked questions by Lianqing Yu, who cannot
% reproduce my figure 12. I cannot find my original code. I am therefore
% reconstructing it.

clear all
close all

nim=5000; % number of images
Pos = inline('x.*(x>0)');

x=[-63:64];
y=[-63:64];
[x2,y2]=meshgrid(x,y);
[xsize,ysize]=size(x2);

SDx = 12;
SDy = 32;
SF = 1/64;

density = 0.25;
dotsize = 5;

SFBW = 2.5;
ORBW = 30*pi/180;
SDx=sqrt(log(2))*1/SF/2/pi*(2^SFBW+1)/(2^SFBW-1)
SDy=1/SF/pi/ORBW*sqrt(log(2))

RFA = exp(-x2.^2/2/SDx^2 - y2.^2/2/SDy^2) .* cos(2*pi*SF.*x2 + pi/4);
RFB = exp(-x2.^2/2/SDx^2 - y2.^2/2/SDy^2) .* cos(2*pi*SF.*x2 - pi/4);

% Gaussian
RFA = exp(-(x2-8).^2/2/SDx^2 - y2.^2/2/SDy^2) ;
RFB = exp(-(x2+8).^2/2/SDx^2 - y2.^2/2/SDy^2) ;

[imageL,imageR] = RandomDotStereogram(xsize,ysize,0,0,density,dotsize);
figure('name','Fig 1')
subplot(1,2,1)
imagesc(imageL)
axis equal off
subplot(1,2,2)
imagesc(-imageR)
colormap gray
axis equal off

figure('name','Fig 9')
subplot(2,2,1)
imagesc(x,y,RFA)
colormap gray
ylabel('vertical position')
subplot(2,2,2)
plot(x,RFA(1,:)./max(squeeze(RFA(1,:))),'k')
axis tight
ylim([-1 1])
subplot(2,2,3)
imagesc(x,y,RFB)
ylabel('vertical position')
xlabel('horizontal position')
subplot(2,2,4)
plot(x,RFB(1,:)./max(squeeze(RFB(1,:))),'k')
xlabel('horizontal position')
axis tight
ylim([-1 1])

xdisparity = [-35 : 5 : 35];

TE = 1; TI = 2; NEAR = 3; FAR = 4; NOTCH1 = 5; NOTCH2 = 6;
CORR = 1; ANTI = 2;


response = zeros(6,2,length(xdisparity));
figure
for jim=1:nim


	% Make new big random-dot pattern
	bigxsize = xsize+max(abs(xdisparity))+1;
	bigimage = zeros(ysize,bigxsize);
	% Decide how many dots to plot:
	ndots = round(density * (bigxsize+2*dotsize) * (ysize+2*dotsize) / min(dotsize,xsize) / min(dotsize,ysize));

	halfdot = fix(dotsize/2);
	xpos = fix(rand(1,ndots)*(bigxsize+2*dotsize)) - halfdot;
	ypos = fix(rand(1,ndots)*(ysize+2*dotsize)) - halfdot;
	col = ((rand(1,ndots)>0.5)-0.5)*2;
	x1 = max(xpos,1);
	x2 = min(xpos + dotsize-1,bigxsize);
	y1 = max(ypos,1);
	y2 = min(ypos + dotsize-1,ysize);

	for j=1:ndots
		bigimage(y1(j):y2(j),x1(j):x2(j) ) = col(j);
	end

	for jd=1:length(xdisparity)
		xdisp = xdisparity(jd);
		%Pick out the right disparities from this big image:
		if xdisp>=0
			imageL = bigimage( : , [1:xsize] );
			imageR = bigimage( :  , [1:xsize] + xdisp );
		elseif xdisp<0
			imageL = bigimage( : , [1:xsize] - xdisp );
			imageR = bigimage( :  , [1:xsize]);
		end

		% Convolve images
		vLA = sum(sum(imageL.*RFA));
		vLB = sum(sum(imageL.*RFB));
		vRA = sum(sum(imageR.*RFA));
		vRB = sum(sum(imageR.*RFB));

		% Disparity-independent
		DIC = vLA^2 + vLB^2 + vRA^2 + vRB^2;

		% Now work out which
		TETI = Pos(vLA).*Pos(vRA) + Pos(-vLA)*Pos(-vRA) + Pos(vLB)*Pos(vRB)  + Pos(-vLB)*Pos(-vRB);
		response(TE,CORR,jd) = response(TE,CORR,jd) + DIC + 2 * TETI;
		response(TI,CORR,jd) = response(TI,CORR,jd) + DIC - 2 * TETI;
		TETI = Pos(vLA).*Pos(-vRA) + Pos(-vLA)*Pos(vRA) + Pos(vLB)*Pos(-vRB)  + Pos(-vLB)*Pos(vRB);
		response(TE,ANTI,jd) = response(TE,ANTI,jd) + DIC + 2 * TETI;
		response(TI,ANTI,jd) = response(TI,ANTI,jd) + DIC - 2 * TETI;

		NF = Pos(vLA)*Pos(vRB) + Pos(-vLA)*Pos(-vRB) - Pos(vLB)*Pos(vRA) - Pos(-vLB)*Pos(-vRA);
		response(NEAR,CORR,jd) = response(NEAR,CORR,jd) + DIC - 2 * NF;
		response(FAR,CORR,jd) = response(FAR,CORR,jd) + DIC + 2 * NF;
		NF = Pos(vLA)*Pos(-vRB) + Pos(-vLA)*Pos(vRB) - Pos(vLB)*Pos(-vRA) - Pos(-vLB)*Pos(vRA);
		response(NEAR,ANTI,jd) = response(NEAR,ANTI,jd) + DIC - 2 * NF;
		response(FAR,ANTI,jd) = response(FAR,ANTI,jd) + DIC + 2 * NF;

		Notch = Pos(vLA)*Pos(vRB) + Pos(-vLA)*Pos(-vRB) + Pos(vLB)*Pos(vRA) + Pos(-vLB)*Pos(-vRA);
		response(NOTCH1,CORR,jd) = response(NOTCH1,CORR,jd) + DIC + 2 * Notch;
		response(NOTCH2,CORR,jd) = response(NOTCH2,CORR,jd) + DIC - 2 * Notch;
		Notch = Pos(vLA)*Pos(-vRB) + Pos(-vLA)*Pos(vRB) + Pos(vLB)*Pos(-vRA) + Pos(-vLB)*Pos(vRA);
		response(NOTCH1,ANTI,jd) = response(NOTCH1,ANTI,jd) + DIC + 2 * Notch;
		response(NOTCH2,ANTI,jd) = response(NOTCH2,ANTI,jd) + DIC - 2 * Notch;

	end % next disparity

	if fix(jim/20)==jim/20
		subplot(2,3,1)
		hold off
		plot(xdisparity,squeeze(response(TE,CORR,:))/max2d(squeeze(response(TE,:,:))),'ko-','markerfacecol','k')
		hold on
		plot(xdisparity,squeeze(response(TE,ANTI,:))/max2d(squeeze(response(TE,:,:))),'ko:')
		title('"Tuned excitatory"')
		subplot(2,3,4)
		hold off
		plot(xdisparity,squeeze(response(TI,CORR,:))/max2d(squeeze(response(TI,:,:))),'ko-','markerfacecol','k')
		hold on
		plot(xdisparity,squeeze(response(TI,ANTI,:))/max2d(squeeze(response(TI,:,:))),'ko:')
		title('"Tuned inhibitory"')
		subplot(2,3,2)
		hold off
		plot(xdisparity,squeeze(response(NEAR,CORR,:))/max2d(squeeze(response(NEAR,:,:))),'ko-','markerfacecol','k')
		hold on
		plot(xdisparity,squeeze(response(NEAR,ANTI,:))/max2d(squeeze(response(NEAR,:,:))),'ko:')
		title('"Near"')
		subplot(2,3,5)
		hold off
		plot(xdisparity,squeeze(response(FAR,CORR,:))/max2d(squeeze(response(FAR,:,:))),'ko-','markerfacecol','k')
		hold on
		plot(xdisparity,squeeze(response(FAR,ANTI,:))/max2d(squeeze(response(FAR,:,:))),'ko:')
		title('"Far"')
		subplot(2,3,3)
		hold off
		plot(xdisparity,squeeze(response(NOTCH1,CORR,:))/max2d(squeeze(response(NOTCH1,:,:))),'ko-','markerfacecol','k')
		hold on
		plot(xdisparity,squeeze(response(NOTCH1,ANTI,:))/max2d(squeeze(response(NOTCH1,:,:))),'ko:')
		title('"Notch"')
		subplot(2,3,6)
		hold off
		plot(xdisparity,squeeze(response(NOTCH2,CORR,:))/max2d(squeeze(response(NOTCH2,:,:))),'ko-','markerfacecol','k')
		hold on
		plot(xdisparity,squeeze(response(NOTCH2,ANTI,:))/max2d(squeeze(response(NOTCH2,:,:))),'ko:')
		title('"Notch"')
		
		for j=1:6
			subplot(2,3,j)
			axis tight
		end
		drawnow
		jim
	end
end

save Fig10_GaussSDx12