% Load images (n=303)
% Natural images: 1-300
% Black: image 301
% Gray:  image 302
% White: image 303
imagedata = h5read('responses2naturalimages.h5','/images');

% Load spike rasters and receptive field information
spikedata = h5read('responses2naturalimages.h5','/spikes');
mudata = h5read('responses2naturalimages.h5','/mu'); % Center coordinates
sigmadata = h5read('responses2naturalimages.h5','/sigma'); % Sigma matrices (convariance matrices)

% Prepare figure
figure;
set(gcf,'Position', [100 100 100+880 100+300]);

% IDs for sample cell and image of Fig. 1 in manuscript
CellID = 44;
ImageID = 228;

% Plot sample image
subplot(1,2,1)
thisimage = squeeze(imagedata(:,:,ImageID));
imagesc(thisimage);
hold on;
set(gca,'YDir','normal');
axis([0, 256, 0, 256] );
axis equal;
xlabel('Pixels')
ylabel('Pixels')
colormap gray;
caxis([-1 1])

% Plot RF outline at 3 sigma
mu = mudata(:,CellID);
sigma = sigmadata(:,:,CellID);
i = zeros(2,1e3);
for j=1:size(i,2)
    alpha = 2*pi*(j-1)/(size(i,2)-1);
    i(:,j) = 3*sqrtm(sigma)*[cos(alpha);sin(alpha)]+mu;
end
plot(i(1,:),i(2,:),'b','LineWidth',2);

% Get spike rasters for sample cell and sample image
data = squeeze(spikedata(:,:,:,CellID));
Raster = squeeze(data(:,:,ImageID));

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

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