% load stimulus images
% Natural images 1-300
% Black image 301
% Gray  image 302
% While image 303
load('ImageData.mat');
imSize = 256; % 256x256 pixels

% load Raster spikes and Receptive fields
load('DATA.mat');

% example cell and image in Fig.1
CellID = 44;
ImageID = 228;

figure;
set(gcf,'Position', [120 100 120+1000 100+850 ]);

subplot(2,2,1)
% Generate stimulus image
buf = ImageData(:,ImageID);
buf = reshape(buf,imSize,imSize);
imagesc(buf);
hold on;
set(gca,'YDir','normal');
colormap gray;
caxis([-0.5 0.5])

% Generate RF outline with 3 standard deviation
RF = RF_parameter(CellID);
i = zeros(2,1e3);
for j=1:size(i,2)
    alpha = 2*pi*(j-1)/(size(i,2)-1);
    i(:,j) = 3*sqrtm(RF.sigma)*[cos(alpha);sin(alpha)]+RF.mu;
end
plot(i(1,:),i(2,:),'b','LineWidth',2);

% Get Raster
data = squeeze(Raster_data(CellID,:,:,:));
Raster = squeeze(data(ImageID,:,:));

% Get all the trials for this cell
NTrial = nnz(sum(Raster,1));
Raster = Raster(:,1:NTrial)';

subplot(2,2,[2 4])
imagesc(Raster)
set(gca,'YDir','normal');
xlabel('Time')
ylabel('Trial')