%MakeFigure1.m
%Makes figure 1 in the commentary
%"Imaging orientation selectivity: decoding conscious perception in V1"
%Nature Neuroscience 2005;8 541-542.
%
%Written by G.M.Boynton, 4/2005

%initial parameters:
FOV = [3*3,3*3]; %Field of view (mm)
pixpermm = 25;  %Number of pixels per millimeter
voxSize = 3; %Voxel size (mm)

sz = pixpermm*FOV;  %Total size of image (in pixels)

angleList = linspace(-180,180,9);   %this determines the histogram bins
%9 numbers define 8 bins

%Fist, make fake orientation 'pinwheel' map based on
%Rojer and Schwartz' method of bandpassing random noise:

%Rojer, A.S. and E.L. Schwartz, Cat and monkey cortical columnar patterns
%modeled by bandpass-filtered 2D white noise. Biol Cybern, 1990. 62(5): p. 381-91.

%Make random noise: complex numbers where the angle is the orientation
z = exp(sqrt(-1)*rand(sz)*pi*2);

%Make a gabor filter
filtSz = 3; % 3mm
sig = .8; %  .8 mm
freq = 1; %cycles/mm (Try zero for big columns)
filtPix = ceil(filtSz*pixpermm);
[x,y] = meshgrid(linspace(-filtSz/2,filtSz/2,filtPix),linspace(-filtSz/2,filtSz/2,filtPix));
r = sqrt(x.^2+y.^2);
filt = exp(-r.^2/sig.^2).*cos(2*pi*freq*r);  %Gabor

%Convolve z with the filter
bandpassz = conv2(z,filt,'same');

%Show the pinwheel map in figure 1
figure(1)
clf
%scale it from 0 to 64
img = 64*(angle(bandpassz)+pi)/(pi*2);

%These are vectors for x and y axes, in units of millemters:
xax = linspace(.5/pixpermm,FOV(2)-.5/pixpermm,FOV(2)*pixpermm);
yax = linspace(.5/pixpermm,FOV(1)-.5/pixpermm,FOV(1)*pixpermm);

image(xax,yax,img);

set(gca,'XTick',[0:FOV(2)]);
set(gca,'YTick',[0:FOV(1)]);
colormap(hsv(64))
axis equal
axis tight
axis off
set(gca,'YDir','normal');
set(gcf,'PaperPosition',[1,1,3,3]);

hold on

%Draw the black voxel edges on the pinwheel map in figure 1.

m = FOV(1)/voxSize; %number of voxel rows
n = FOV(2)/voxSize; %number of voxel columns

%loop through the voxels
for i=1:m  %columns
    for j=1:n %rows
        xx = (j-1)*voxSize;
        yy = (i-1)*voxSize;
        patch([xx,xx+voxSize,xx+voxSize,xx,xx,xx+voxSize],[yy,yy,yy+voxSize,yy+voxSize,yy,yy],'w','LineWidth',3,'FaceColor','none');
    end
end


%On to the histograms...
%Figure 2 will show histogrrams of the amount of orientation represented in
%each bin, for each voxel.

xa = linspace(.1,.9,length(angleList));  %x-positions for bar graphs
%colors for each bars will match center of colors for each bin

tmpCol = hsv(length(angleList)*2-1);

col = tmpCol(2:2:end,:);

figure(2)
clf
axis off
axis equal
hold on

%loop through the voxels
for i=1:m  %columns
    for j=1:n %rows
        %define out the pixels inside the voxel 
        xx = (j-1)*voxSize;
        yy = (i-1)*voxSize;
        inVoxel = bandpassz(yy*pixpermm+1:(yy+voxSize)*pixpermm,xx*pixpermm+1:(xx+voxSize)*pixpermm);
        a = 180*angle(inVoxel(:))/pi;
        %calculate the pixels for each bin
        for k=1:length(angleList)-1
            nn(k) = sum(a>=angleList(k) & a<angleList(k+1));
        end
        %scale
        nn = .5*nn*(length(angleList)-1)/length(inVoxel(:));
        %draw the bins
        for k=1:length(angleList)-1
            patch(j-1+[xa(k),xa(k+1),xa(k+1),xa(k),xa(k)],i-1+[0,0,nn(k),nn(k),0],col(k,:));
        end
        %black box around the histogram:
        patch([j-1,j,j,j-1,j-1,j],[i-1,i-1,i,i,i-1,i-1],'w','LineWidth',3,'FaceColor','none');
    end
end
set(gcf,'PaperPosition',[1,1,3,3]);


%make a legend of orientation colors

figure(3)
clf
hold on
for i=1:length(angleList)-1
    a = pi*(i-1)/(length(angleList)-1);
    xc = i*2.25;
    yc = 1;

    plot([xc-cos(a),xc+cos(a)],[yc-sin(a),yc+sin(a)],'-','LineWidth',8,'Color',col(i,:));
end

axis equal
axis off