%%% Code for exploring the ideas in the Kamitani & Tong paper
%%% by Rajeev Raizada, Oct.2008
%%%
%%% Links to the relevant papers:
%%% The Kamitani & Tong paper:
% http://www.psy.vanderbilt.edu/tonglab/publications/Kamitani&Tong_NN2005.pdf
%%% Geoff Boynton's commentary, which has the figure
%%% which this code is based on:
% http://faculty.washington.edu/gboynton/publications/boynton-natureneuro05.pdf
%%% Rojer & Schwartz's paper, which showed how to make
%%% patterns which look like V1 by filtering noise:
% http://eslab.bu.edu/publications/articles/1990/rojer1990cat.pdf

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% First of all, we want to make a grid of random orientations.
% Then, we're going to blur those orientations
% by averaging neighbouring orientations together.
% This will make a pattern which looks very much like a real
% V1 orientation map.
%
% However, averaging together orientations is not quite as straightforward
% as just averaging together a bunch of scalar numbers,
% because orientations loop back on themselves.
% (A scalar is a number which has just one dimension.
%  A vector has two or more dimensions, so it has a direction
%  as well as magnitude).
% The average of a 1-degree orientation and a 179 degree-orientation
% is *not* equal to (1+179)/2 = 90 degrees.
% Both 1deg and 179deg are very close to horizontal,
% and the average of them is exactly horizontal, namely 0 degrees.
% So, we have to treat the orientations as vectors
% when we average them together, not just as scalar numbers.


% First, let's make a grid of random orientation angles 
% between 0 and 2*pi ( that's 360 degrees, in radians).

% We'll call the orientations "theta",
% as people usually use that greek letter for orientation
grid_size = 300;
theta = 2*pi*rand(grid_size);

% Turn these angles into vectors with x and y coords
theta_x = cos(theta);
theta_y = sin(theta);

% Make a Gaussian, for blurring the angles together
% Matlab has a useful built-in function for making 2D grids: "meshgrid"
gaussian_grid_step_size = 0.1;
% A Gaussian is pretty much zero by the time
% it gets to 5 standard deviations from the mean
[x_grid,y_grid] = meshgrid( -5 : gaussian_grid_step_size : 5 );

gaussian = exp( -( x_grid.^2 + y_grid.^2 ) );

% Locally average the vector coords with a Gaussian blur
% The 'valid' part at the end tells Matlab to throw away
% the parts of the output where the Gaussian is hanging off an edge
blurred_theta_x = conv2(theta_x,gaussian,'valid'); 
blurred_theta_y = conv2(theta_y,gaussian,'valid');

% Turn these vectors back into angles
tan_of_blurred_theta = blurred_theta_y./(blurred_theta_x); % tan = sin/cos
blurred_theta = atan( tan_of_blurred_theta );

% Let's plot these orientation fields, with the axis scaled
% to be in units of pi
blurred_theta_pi_units = blurred_theta / pi;
theta_pi_units = theta / pi;

figure(1);
clf;
imagesc(theta_pi_units);
colormap hsv; % hsv is a colormap where the max and min are both red
colorbar;
caxis([-0.5 0.5]);
title('Random orientations, before blurring');

figure(2);
clf;
imagesc(blurred_theta_pi_units);
colormap hsv; 
colorbar;
caxis([-0.5 0.5])
title('Random orientations, after blurring');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Now we're going to slice the space up into voxels,
%%%% and see what's in each voxel
    
%%% Slice the space up into voxels, and histogram each of these
num_voxels_per_dim = 4;
voxel_size = size(blurred_theta,1) / num_voxels_per_dim;

voxel_contents = zeros(num_voxels_per_dim,num_voxels_per_dim,voxel_size^2);

for x = 1:num_voxels_per_dim,
    for y = 1:num_voxels_per_dim,
        
        this_voxel = blurred_theta( (x-1)*voxel_size + [1:voxel_size], ...
                                    (y-1)*voxel_size + [1:voxel_size] );
        voxel_contents(x,y,:) = this_voxel(:);
    end;
end;

%%%% Put thick black lines showing the voxels on Fig 2
figure(2);
hold on;
for x = 1:(num_voxels_per_dim-1),
    for y = 1:(num_voxels_per_dim-1),
        
        plot(x*voxel_size*[1 1],[0 num_voxels_per_dim*voxel_size],'k-','LineWidth',3);
        plot([0 num_voxels_per_dim*voxel_size],y*voxel_size*[1 1],'k-','LineWidth',3);

    end;
end;

%%% Plot these histograms with the right colors
%%% We'll histogram the voxels into 8 bins
num_orientation_bins = 8;
bin_floors = (pi * [0:(num_orientation_bins-1)]/num_orientation_bins) - pi/2;
bin_ceilings = bin_floors + pi/num_orientation_bins;
bin_centers = (bin_floors + bin_ceilings) / 2;

%%% We'll set the activation of each voxel for each orientation
%%% to be the number of pixels in each orientation-bin in each voxel
%%% First, initialise this to be all zeros
voxel_activations = zeros(num_voxels_per_dim,num_voxels_per_dim,num_orientation_bins);

figure(3);
clf;
colormap hsv;
for x = 1:num_voxels_per_dim,
    for y = 1:num_voxels_per_dim,
        
        subplot(num_voxels_per_dim,num_voxels_per_dim,(x-1)*num_voxels_per_dim+y);
    
        freq = hist(voxel_contents(x,y,:),bin_centers);
        %%% We'll set the activation of each voxel for each orientation
        %%% to be the number of pixels in each orientation-bin in each voxel
        voxel_activations(x,y,:) = freq;

        bar_h = bar(bin_centers,freq);
        
        %%% Set the bar colours using a trick from here:
        %%% http://www.mathworks.com/support/solutions/data/1-1T6UHS.html
        bar_children = get(bar_h,'Children');
        set(bar_children,'CData',bin_centers);
        title([ num2str(x) ',' num2str(y) ]);
        caxis([-pi/2 pi/2]);
        axis('off');
        
    end;
end;

%%%% Set the y-scaling on the histograms to go up to the maximum activation
max_activation = max(voxel_activations(:));
figure(3);
for x = 1:num_voxels_per_dim,
    for y = 1:num_voxels_per_dim,
        subplot(num_voxels_per_dim,num_voxels_per_dim,(x-1)*num_voxels_per_dim+y);
        axis([-pi/2 pi/2 0 max_activation]);
    end;
end;


%%%%% Let's draw the orientations corresponding to each histogram bin
figure(4);
clf;
hold on;
hsv_color_mat = colormap(hsv); %%% Colormap matrices have 64 rows
for this_bin = 1:num_orientation_bins,
    this_orientation = bin_centers(this_bin);
    hsv_mat_row = round( 64 * (this_orientation/pi + 0.5) );
    plot(cos(this_orientation)*[-1 1],sin(this_orientation)*[-1 1], ...
        'Color',hsv_color_mat(hsv_mat_row,:),'LineWidth',5);
end;
axis('off');
axis('equal');
title('Which colour stands for which orientation','FontSize',20);
   
%%% Now let's present one orientation at a time,
%%% and look at the activation across the voxels.
%%% We'll show the orientation for each bin

figure(5);
clf;
hsv_with_black_first_row = hsv_color_mat;
hsv_with_black_first_row(1,:) = [ 0 0 0 ];
colormap( hsv_with_black_first_row );  %%% Make the 1st row black

for this_bin = 1:num_orientation_bins,
    
    orientation_floor = bin_floors(this_bin);
    orientation_ceiling = bin_ceilings(this_bin);
    
    activated_pixels =  ( blurred_theta >= orientation_floor ) .* ...
                        ( blurred_theta < orientation_ceiling );

    subplot(3,3,this_bin),
        hsv_mat_row = round( 64 * (bin_centers(this_bin)/pi + 0.5) );
        image( hsv_mat_row * activated_pixels);
    axis('equal');
    axis('off');
    %%%% Put thick black lines showing the voxels onto each subplot
    hold on;
    for x = 1:(num_voxels_per_dim-1),
        for y = 1:(num_voxels_per_dim-1),

            plot(x*voxel_size*[1 1],[0 num_voxels_per_dim*voxel_size],'w-','LineWidth',1);
            plot([0 num_voxels_per_dim*voxel_size],y*voxel_size*[1 1],'w-','LineWidth',1);

        end;
    end;
end;


%%%%% Let's show the patterns of voxel activation for each orientation
%%%%% which we calculated above
figure(6);
clf;
colormap gray;
for this_bin = 1:num_orientation_bins, 
    subplot(3,3,this_bin),
    imagesc(voxel_activations(:,:,this_bin));
    caxis([0 max_activation]);
    axis('equal');
    axis('off');
end;
subplot(3,3,2),
title('Spatial patterns of voxel activation for each orientation','FontSize',18); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Now let's do some basic classification of those spatial patterns,
%%%%% just using correlation.

% Suppose we are given a pattern of voxel fMRI activation.
% Can we figure out which orientation gave rise to it?
% We can do this by calculating the correlation of that pattern
% with each of our original measured activation patterns,
% and see which one it matches best with.

% First, we'll turn each activation-pattern-grid into a vector
activation_vectors = zeros(num_voxels_per_dim^2, num_orientation_bins);

for this_bin = 1:num_orientation_bins,
    this_activation_grid = voxel_activations(:,:,this_bin);
    %%% Putting (:) after a matrix stacks it up into a single vector
    activation_vectors(:,this_bin) = this_activation_grid(:);
end;

%%%% Let's look at the correlation matrix of these vectors
figure(7);
clf;
colormap gray;
imagesc( corr(activation_vectors) );
colorbar;
title('Correlations of each orientation''s fMRI-pattern with the others','FontSize',18); 

%%%%% Each pattern is perfectly correlated with itself, of course.
%%%%% Without any noise, we will get perfect classification.
%%%%% But fMRI has lots of noise in it.
%%%%% So, let's add in some noise, and see how the correlations end up.

noise_factor = 0.5;
vec_length = num_voxels_per_dim^2;
noisy_vectors = zeros(vec_length,num_orientation_bins);

for this_bin = 1:num_orientation_bins,
    noise_multiplier = 1 + noise_factor*randn(vec_length,1);
    this_vector = activation_vectors(:,this_bin);
    this_noisy_voxel_vec = this_vector .* noise_multiplier;
    noisy_vectors(:,this_bin) = this_noisy_voxel_vec;
end;

%%%% Let's look at the correlation matrix of the noisy vectors
%%%% with the original activation vectors
figure(8);
clf;
colormap gray;
noisy_corrs = corr(activation_vectors,noisy_vectors);
imagesc( noisy_corrs );
colorbar;
title('Correlations of the noisy activation vectors with the originals','FontSize',18); 

%%%% For each noisy vector, figure out which original
%%%% orientation-activation-vector it has the highest correlation with
best_matching_orientations = zeros(1,num_orientation_bins);

for this_bin = 1:num_orientation_bins,
    corrs_of_this_noisy_vector = noisy_corrs(:,this_bin);
    [ max_corr, best_match_ind  ] = max( corrs_of_this_noisy_vector );
    best_matching_orientations( this_bin ) = best_match_ind;
end;

%%%% Now let's do this 100 times, and see how often the best match 
%%%% is correct, for our particualr noise-factor
noise_factor = 1;
num_reps = 100;
best_matches_rec = zeros(num_reps,num_orientation_bins);
for rep_num = 1:num_reps,
    for this_bin = 1:num_orientation_bins,
        noise_multiplier = 1 + noise_factor*randn(vec_length,1);
        this_vector = activation_vectors(:,this_bin);
        this_noisy_voxel_vec = this_vector .* noise_multiplier;
        noisy_vectors(:,this_bin) = this_noisy_voxel_vec;
    end;
    noisy_corrs = corr(activation_vectors,noisy_vectors);
    for this_bin = 1:num_orientation_bins,
        corrs_of_this_noisy_vector = noisy_corrs(:,this_bin);
        [ max_corr, best_match_ind  ] = max( corrs_of_this_noisy_vector );
        best_matching_orientations( this_bin ) = best_match_ind;
    end;
    best_matches_rec(rep_num,:) = best_matching_orientations;
end;
        
%%%% Let's make polar plots of how often each orietnation
%%%% gave the best match, like in Fig.2 of Kamitani & Tong
figure(9);
clf;
for this_bin = 1:num_orientation_bins,
    %%%% First, let's count up the best matches for each orientation
    best_matches_for_this_orientation = best_matches_rec(:,this_bin);
    best_matches_count = hist(best_matches_for_this_orientation, ...
                              [1:num_orientation_bins]-0.5);
    %%%% Now let's plot this as a polar plot 
    %%%% The orientation bins only go from -pi/2 to pi/2,
    %%%% but the polar plot goes from -pi to pi, 
    %%%% so we'll plot it twice back-to-back.
    %%%% Actually, we'll do it three times as a hack,
    %%%% as it makes the plot lines join up into a circle and look nice
    subplot(3,3,this_bin),
        h=polar([bin_centers pi+bin_centers bin_centers],repmat(best_matches_count,1,3));
        set(h,'LineWidth',3);
end;