V1_ReverseCorrelation/plot_orientation_tuning_shift.m

function [rcJointOrientPlot,shift]=plot_orientation_tuning_shift(w,rcJointOrient,plotNodes,plotTime,angleStep,rcOrient)

numNodes=length(rcJointOrient);
[numBins,numBins,numTimes]=size(rcJointOrient{1});
maxTime=floor(numTimes/2);
if nargin<3
  plotNodes=1:numNodes
end
if nargin<4
  plotTime=3;
end
time=plotTime+maxTime+1;
if nargin<5
  angleStep=15;
end
angle_values=[0:angleStep:180-angleStep]

%subtract mean response to ori1 in order to isolate the response to ori2, as
%described in the Methods of Felsen et al, 2002.
for node=1:numNodes,
  rcJointOrientPlot{node}=rcJointOrient{node};
  for ori_count=1:numBins
	for toff_count=1:numTimes
	  if nargin<6
		rcJointOrientPlot{node}(1:numBins,ori_count,toff_count)=...
			rcJointOrient{node}(:,ori_count,toff_count)- ...
			mean(rcJointOrient{node}(:,:,toff_count),2);
	  else
		rcJointOrientPlot{node}(1:numBins,ori_count,toff_count)=...
			rcJointOrient{node}(:,ori_count,toff_count)- ...
			rcOrient{node}(:,toff_count);
	  end
	end
  end
end
%figure(5),clf
%plot_reverse_correl_joint_orient(w,rcJointOrientPlot,angle_values,plotNodes)

%rotate orientation tuning curves so that orientation preference is approx
%central: this will allow a Gaussian to be fitted to the data. Absolute
%orientation preference is irrelevant to the following analysis of shift.
numOris=8;
cent=5; %1st node with a central orientation preference
oriVals=zeros(numNodes,numBins);
for node=1:numNodes,
  rcJointOrientPlotCentred{node}(:,:)=...
	  rcJointOrientPlot{node}(:,:,time);
  
  ori=mod(node-1,numOris)+1;
  cut=fix((cent-ori)*numBins/numOris);
  if cut>0
	range1=numBins-cut+1:numBins;
	range2=1:numBins-cut;
	%data=[data(numBins-cut+1:numBins),data(1:numBins-cut)];
	
  elseif cut<0
	cut=abs(cut);
	range1=cut+1:numBins;
	range2=1:cut;
	%data=[data(cut+1:numBins),data(1:cut)];
  end	
  if cut~=0
	%centre data for ori2
	rcJointOrientPlotCentred{node}(:,1:length(range1))=...
		rcJointOrientPlot{node}(:,range1,time);
	rcJointOrientPlotCentred{node}(:,length(range1)+1:numBins)=...
		rcJointOrientPlot{node}(:,range2,time);
	
	%centre data for ori1 by applying same shift
	rcJointOrientPlotCentred_safe{node}=rcJointOrientPlotCentred{node};
	rcJointOrientPlotCentred{node}(1:length(range1),:)=...
		rcJointOrientPlotCentred_safe{node}(range1,:);
	rcJointOrientPlotCentred{node}(length(range1)+1:numBins,:)=...
		rcJointOrientPlotCentred_safe{node}(range2,:);
  end
end

figure(6),clf
prefOri=zeros(1,numNodes);
for node=1:numNodes,
  data=mean(rcJointOrientPlotCentred{node}(:,:),1);
  data=(data-min(data)).*10000000;
  
  [A, mu, sigma, error] = fitgauss(data,angle_values);
  %mu = preferred orientation of 2nd stimulus (in degrees relative to re-centred
  %coordinates)
  if sigma>0 & error<0.1  	  
	prefOri(node)=mu; 
  else
	prefOri(node)=NaN;
  end
 
  maxsubplot(5,8,node);
  plot(10000000.*mean(rcJointOrientPlot{node}(:,:,time),1),'b');hold on,plot(data,'g');hold on,plot([1:numBins],A*exp(-0.5*((angle_values-mu)./sigma).^2)/sqrt(2*pi*sigma^2),'m');
end

figure(7),clf
ori1BinCentres=[-82.5:15:82.5];
numBinsShift=length(ori1BinCentres);
shift=ones(numNodes,numBinsShift).*NaN;
for node=plotNodes
  maxsubplot(5,8,node);
  for ori1=1:numBins, 
	data=rcJointOrientPlotCentred{node}(ori1,:);
	data=(data-min(data)).*10000000;
	[A, mu, sigma, error] = fitgauss(data,angle_values);
	%mu = peak of orientation tuning curve for a particular value of ori1 (in
    %degrees relative to re-centred coordinates)

	shiftVal=mu-prefOri(node);%tuning to ori1 - mean tuning to ori2 (in
							  %degrees)
	
	Ori1Diff=angle_values(ori1)-prefOri(node);%orientation of stimulus 1 -
                                              %mean tuning to ori2 (in degrees)
	if Ori1Diff<-90
	  Ori1Diff=90+(Ori1Diff+90);
	end
	if Ori1Diff>90
	  Ori1Diff=-90+(Ori1Diff-90);
	end
	binNum=binned_param(Ori1Diff,[-90,90],numBinsShift);

	if sigma>0 & error<0.1 & binNum>0
	  shift(node,binNum)=shiftVal;
	  plot(data,'g');hold on,plot([1:numBins],A*exp(-0.5*((angle_values-mu)./sigma).^2)/sqrt(2*pi*sigma^2),'m');
	end
  end
end
figure(8),clf
for node=plotNodes,
  maxsubplot(5,8,node);
  barh(fliplr(shift(node,:)));axis([-7.5,7.5,0,numBinsShift+1]);
end


figure(9),clf
shiftMean=nanmean(shift,1);
shiftStd=nanstd(shift,1,0);
shiftPop=sum(~isnan(shift),1)
shiftSEM=shiftStd./sqrt(shiftPop);
barh(fliplr(shiftMean));
hold on
for bin_count=1:numBinsShift
  plot([-shiftSEM(bin_count),shiftSEM(bin_count)]+shiftMean(bin_count),...
	   [1,1].*(numBinsShift+1-bin_count),'k','LineWidth',2);
end
set(gca,'FontSize',22);axis('square')
axis([-5,5,0,numBinsShift+1])
h = findobj(gca,'Type','patch');
set(h,'FaceColor','w','EdgeColor','k','LineWidth',3);
xlabel('shift (degrees)','FontSize',28)
ylabel('\theta_1 (degrees)','FontSize',28)
set(gca,'YTick',[0.5:numBinsShift/4:numBinsShift+0.5],'YTickLabel',[90:-45:-90]);
set(gcf,'PaperPosition',2.55.*[1 1 6 6]);