Natural movies data and code

naturalMovies/code natural movies/BC_getNaturalMovies.m

@@ -0,0 +1,123 @@
+function [movieSignal]=BC_getNaturalMovies(pixelSize)
+%[naturalMovie_output]=BC_getNaturalMovies(loadData)
+%this function shows the prinicpal calculation for how to preprocess the
+%natural movies from the "CatCam" database. The movie intensity values are
+%scaled to have the same mean light intensity as white noise stimulus and 
+%converted to WeberContrast values. Further, the movie dimensions are
+%adapted to match the dimension of the white noise STA filter to later 
+%perform the prediction with the linear-nonlinear model (LN-model)computed 
+%with white noise stimulus.
+%Inputs:
+%               pixelSize = pixel size from the white noise stimulus
+%
+%Outpus:        movieSignal = preprocessed movie signal, ready for the
+%               LN-model analysis
+%
+%Code by HS Feb 2021
+
+%-------------------------------------------------------------------------------
+%% Constants
+%movie dimensions of the original tiff files
+NxMovie=320; %x pixel tiff movies
+NyMovie=240; %y pixel tiff movies
+%spatiotemporal white noise dimensions x-y (they are all the same for all
+%recorded cells)
+monitorXPixels=800; %x pixels of OLED monitor; 
+monitorYPixels=600; %y pixels of OLED monitor; 
+NxWhiteNoise=ceil(monitorXPixels/pixelSize); %number of squares shown X dimension
+NyWhiteNoise=ceil(monitorYPixels/pixelSize); %number of squares shown Y dimension
+
+%-------------------------------------------------------------------------------
+%% 1. Preprocess the movie contrast values
+%the movies can be downloaded under the following
+%link:https://zenodo.org/record/46481. Each movie frame is presented as an 
+%individual .tiff file. See documentation.data.pdf for
+%further instructions, also which tiff-filename to take per movie.
+%Ones the files are downloaded they have to be prepared as follow:
+
+%to show how the movie signal is preprocessed, I generated a random
+%example signal of the same dimension as the natural movies (320x240)
+nbrFrames=998;
+genRandSignal=rand(NyMovie,NxMovie,nbrFrames)*255; %a movie signal (0-255 intensity values, as movie itself)) of 240 (y-dim) 320 (x-dim) and 998 (time frames); thats where you have to put the original movie signal at 320x240xnumber of frames 
+% %%code to load the tiff files in 1 to 998
+% loadTiff='Y:\FromPeopleToPeople\Helene2Tim\tiffMovie4';
+% indx=1;
+% ContentStackDir_tif=dir([loadTiff '\']); %what is inside the folder
+% for l=1:1000 %nbr of files in the folder
+%     if numel(strfind(ContentStackDir_tif(l).name,'.tif'))>0
+%         %tif files containes one value per pixel -> it is an intensity
+%         %value for grayscale. So no color information. It is a in uint16
+%         %imfinfo([ StackDir ImageStack.name '\' ContentStackDir(l).name ])
+%         %gives you the infos.
+%         %16-bit intensity images is usually [0, 65535].
+%         %http://de.mathworks.com/help/matlab/creating_plots/image-types.html
+%         [ImageStack_tif.imagesOri{indx},cmp]=imread([loadTiff '\' ContentStackDir_tif(l).name ]);%read tiff file
+%         indx=indx+1;
+%     end
+% end
+% genRandSignal=cell2mat(ImageStack_tif.imagesOri);
+oneDIamge = reshape(genRandSignal,[],1);
+doubleIamge=double(oneDIamge);
+normImage=doubleIamge/255; %get intensity values between 0-1
+clear oneDIamge doubleIamge
+%scale the light intensity to have the same mean as the spatiotemporal white noise stimulus
+%which is: 0.5 (see also Method of the mansucript by Schreyer & Gollisch,
+%2021)
+a=mean(normImage);
+dummy = normImage*(0.5/a);
+%correct contrast values if needed (this part is done in stimulator program)
+dummy(dummy<0) = 0;
+dummy(dummy>1) = 1;
+
+%change the vlaues to weber contrast +/1
+meanInt=mean(dummy);
+movie_WebCont=(dummy-meanInt)/meanInt;
+%reshape signal to have again a 3 D structure
+movie_WebCont3D = reshape(movie_WebCont,[NyMovie,NxMovie,nbrFrames]);
+clear movie_WebCont dummy oneDIamge
+
+%------------------------------------------------------------------------------------------
+%% 2. cut the natural movie to have the same dimension as the spatiotemporal white noise
+%the movie has a dimension of 320x240, to present the movie on the retina,
+%the movie was upsampled by 3 (960x720) and cut out to the OLED resolution
+%of 800x600 (cut is from the upper left ->see Documentation.Data.pdf).
+%The spatiotemporal white noise was shown at 89x67 checker size resolution.
+%And the filter (STA) is at that resolution. We want now to have the same
+%resolution for the movies (89x67) so that we can easily convolve the STA with the
+%movie.
+%To get the dimensions, there are multiple ways. The most intuitive way and easiest to understand is:
+%First upsample the image by 3 and cut it to 800x600 OLED resolution, then make
+%a mean per 9x9 pixels (for all cells with movie the checker size of the
+%white noise stimulus was 9x9 -> see also Documentation_Data.pdf). Then you 
+%get the same resolution as for the white noise stimulus (89x67).
+%Alternative way (maybe faster from the coding side and what I show here) is to directly make the
+%average over 3x3 pixels (dont upsample it first by 3) and then cut it to
+%89x67 resolution to match the white noise STA (filter).
+
+%%make average per 3x3 movie pixels
+NyMovie2=ceil(NyMovie/3);
+NxMovie2=ceil(NxMovie/3);
+indx=1;%x dimension
+indx2=1;%y dimension
+movie_WebContSmall=nan(NyMovie2,NxMovie2,size(movie_WebCont3D,3));
+for z=1:size(movie_WebCont3D,3) %over movie frames
+    for x=1:3:size(movie_WebCont3D,2)-2 %over x-axis (first three)
+        for y=1:3:size(movie_WebCont3D,1)-2 %do the corresponding y-axis
+            tempx=movie_WebCont3D(y:y+2,x:x+2,z);
+            movie_WebContSmall(indx2,indx,z)=mean(tempx(:)); 
+            indx2=indx2+1;%y dimension
+        end
+        indx=indx+1;
+        indx2=1;
+    end
+    indx=1;
+    indx2=1;
+end
+
+%cut the image to a 89x67 resolution (see Documentation_data.pdf) the cut
+%is done from the upper left corner
+movieSignal=movie_WebContSmall(size(movie_WebContSmall,1)-(NyWhiteNoise)+1:end...
+    ,1:NxWhiteNoise,:); %starts from y pixel 14 (the first 13 are cut out in stimulator program)
+
+end
+

naturalMovies/code natural movies/BC_mainNaturalMovies.m

@@ -8,7 +8,7 @@ clear;
 close all;
 
 %% Preprocessing of the natural movie contrast signal
-loadData='Z:\DataAvailable\naturalMovies\'; %path of the BC files for natural movies (has to be adapted to your own path)
+loadData='Z:\naturalMovies\'; %path of the BC files for natural movies (has to be adapted to your own path)
 filename='16d200016Comp.mat'; %BC file to analyse: for example here we use file 16d200016Comp from cell 5, recorded in retina 2 on 161219
 Hz=25; % given in the excel file LogInfo_NaturalMovies.xlsx for the corresponding file
 movienbr=4; %given in the excel file LogInfo_NaturalMovies.xlsx for the corresponding file

naturalMovies/code natural movies/BC_predNaturalMovies.m

@@ -1,13 +1,14 @@
-function [naturalMovie_output]=BC_predNaturalMovies(loadData,filename,Hz,STA_output,NL_output)
-%[naturalMovie_output]=BC_predNaturalMovies(loadData,filename,pixelSize,Hz,STA_output,NL_output)
+function [naturalMovie_output]=BC_predNaturalMovies(loadData,filename,movieSignal,STA_output,NL_output,Hz)
+%[naturalMovie_output]=BC_predNaturalMovies(movieSignal,STA_output,NL_output)
 %this function shows the prinicpal calculation for the prediction of the 
 %response to the natural movie, from the linear-nonlinear model (LN-model)
 %computed with white noise stimulus.
 %Inputs:
 %               loadData = the folder path of the natural movie files
 %               filename = the filename of the BC of interest
-%               movieData
-%               Hz = stimulus update rate
+%               movieSignal = the preprocessed natural movie signal the
+%               output of the BC_getNaturalMovies.m. Important that the
+%               movie signal has the same dimensions as the STA_output
 %               STA_output =  output of the BC_STA function; containing the
 %               STA (filter) of the white noise stimulus  
 %               NL_output =  output of the BC_NL function; containing the
@@ -16,27 +17,33 @@ function [naturalMovie_output]=BC_predNaturalMovies(loadData,filename,Hz,STA_out
 %Outpus:        naturalMovie_output = structure containing the output of the
 %               prediction
 %
-%
-
-%-------------------------------------------------------------------------------
-%% 1. Get the movie contrast values
-%the movies can be downloaded under the following
-%link:https://zenodo.org/record/46481. See documentation.data.pdf for
-%further instructions, also which frames to take.
-%they have to be prepared as follow:
+% code by HS last modiefied Feb 2021
 
+%% Constants
+findBreak=500; %the movie is repeated multiple times. Between the movies there is 
+%a break/gray screen of >0.5s. To find the start of each movie trial, you find the gray screen between the movies
 
+%-------------------------------------------------------------------------------
+%% 1. cut the natural movie signal to the size of the filter
+%in the BC_STA.m function, the STA was cut. Here the natural movie signal
+%has to be cut to the same dimensions.
+movieSignal_zoom=movieSignal(STA_output.rangeY_zoom,STA_output.rangeX_zoom,:); 
+%make the signal 2D for prediction
+movieSignal2D_zoom=reshape(movieSignal_zoom,size(movieSignal_zoom,1)*...
+    size(movieSignal_zoom,2),size(movieSignal_zoom,3));
 
 %-------------------------------------------------------------------------------
 %% 2. compute the prediction
 %% 2.1 convolution
-%convolve the new stimulus sequence (not used for prediciton) with the
+%normalize the filter
+linearFiler=normSTA(STA_output.STA2D_zoom);
+%convolve the natural movie stimulus sequence  with the
 %normalized STA
-stim2pred=filter2(STA_norm,stimPred,'valid');
+stim2pred=filter2(linearFiler,movieSignal2D_zoom,'valid');
 
 %% 2.2 compute the prediction through linear interpolation
-NL_x=bins.NL_xbin;
-NL_y=bins.NL_ybin;
+NL_x=NL_output.bins.NL_xbin;
+NL_y=NL_output.bins.NL_ybin;
 %example of how to make the prediction through linear interpolation
 predMemP_y=nan(1,size(stim2pred,2));
 for kk=1:size(stim2pred,2)
@@ -51,36 +58,74 @@ for kk=1:size(stim2pred,2)
 end
 
 %-------------------------------------------------------------------------------
-%% 3. evaluate the prediction performance
+%% 4. compute the average response to the move
+%% 4.1 load movie response
+data=load([loadData,filename]); %load the file
+
+%sort the different pulses for the different trials
+diffttls=diff(data.ttls_ms);
+nbrTrials=find(diffttls>findBreak); %to find out which pulse correspond to a new trial (gray break)
+nbrFrames=diff(nbrTrials);
+if isempty(nbrFrames)
+    nbrFrames=nbrTrials;
+end
+%I get the membrane potential per trial & frames per trial
+reOrderFrames=nan(nbrFrames(1),length(nbrTrials));
+membranePerTrial_origRes=cell(1,length(nbrTrials));
+for kk=1:length(nbrTrials)
+    if kk==1 %first trial
+        reOrderFrames(:,kk)=data.ttls_ms(1:nbrTrials(kk));
+        membranePerTrial_origRes{kk}=data.membPot(reOrderFrames(1,kk):reOrderFrames(end,kk)-1); %-1 because the last frame is the start of the gray screen
+    else
+        reOrderFrames(:,kk)=data.ttls_ms(nbrTrials(kk-1)+1:nbrTrials(kk)); %end of previous trial +1 (starts at 999)
+        membranePerTrial_origRes{kk}=data.membPot(reOrderFrames(1,kk):reOrderFrames(end,kk)-1);
+    end
+end
+%I only include trials here that here fully finished until the break!
+
+%------------------------------------------------------------------------------------------
+%% 4.2 Natural Movie Build average per frame
+%make one value for the membrane potential between two stimulus signal (two
+%frames), where the same image was shown
+memPperTrial_Movie=avgPerFrameNatMovies(data.membPot,...
+    reOrderFrames,nbrFrames(1),Hz);
+
+memP_AVGMovie=mean(memPperTrial_Movie,2); 
+
+%------------------------------------------------------------------------------------------
+%% 5 compute the pearson correlation coefficient
 %pearson correlation
-corrcoeff_temp=corrcoef(voltPred(STA_output.STA_fnbr:end),predMemP_y); %the first value is at the length of the filter
+corrcoeff_temp=corrcoef(memP_AVGMovie(STA_output.STA_fnbr:end),predMemP_y); 
 pearCorrSqrt=corrcoeff_temp(2)^2;
+%not with this method the prediction value is a bit lower than what we report in the
+%manuscript. To get the same values as in the mansucript, you have to
+%cut 3sigma around the RF and perform the SVD (see also Manuscript and
+%BC_STA, BC_NL)
 
 %-------------------------------------------------------------------------------
-%% 4. plot the output
+%% 6. plot the output
 figure;
-%STA
-subplot(2,2,1)
-[~,indx_Zoom]=max(abs(STA3D_norm(:))); %max intensity value of STA
-[~,~,posZ]=ind2sub(size(STA3D_norm),indx_Zoom); %get the frame position
-imagesc(STA3D_norm(:,:,posZ))
+%plot the first frame of the move
+subplot(2,1,1)
+imagesc(movieSignal(:,:,1))
 colormap('gray'); %how the sitmulus was shown
 axis equal; axis tight;
-title('STA best frame')
-%NL
-subplot(2,2,2)
-plot(NL_xAxis,NL_yAxis,'.','color','k') %plot non-binned signal
-hold on;
-plot(bins.NL_xbin,bins.NL_ybin,'o','markerfacecolor','r','markeredgecolor','r') %plot binned signal
-title('NL')
-axis tight
-xlabel('generator signal')
-ylabel('voltage (mV)')
-%prediction
-subplot(2,2,3:4)
-plot(voltPred(STA_output.STA_fnbr:end),'k') %original membrane potential
+xticks([]); yticks([]);
+
+%plot the prediction trace
+subplot(2,1,2)
+%only for visual comparison we normalize the predicted trace to have the
+%same mean and std as the response to the movie.
+x=predMemP_y;
+m2=mean(memP_AVGMovie);
+s2=std(memP_AVGMovie);
+s1=std(x);
+m1=mean(x);
+corrPredMemPred=m2+(x-m1)*(s2/s1);
+
+plot(memP_AVGMovie(STA_output.STA_fnbr:end),'k') %original membrane potential
 hold on;
-plot(predMemP_y,'r')%predicted membrane potential
+plot(corrPredMemPred,'r')%predicted membrane potential
 title(['pCorr:' mat2str(round(pearCorrSqrt,2)) ' Filename: ' STA_output.filename])
 xlabel('stimulus frames')
 ylabel('voltage (ms)')
@@ -88,8 +133,8 @@ legend('original response','predicted response')
 
 %-------------------------------------------------------------------------------
 %% 5. add variables into the output
-pred_output.pearCorrSqrt=pearCorrSqrt;
-
+naturalMovie_output.pearCorrSqrt=pearCorrSqrt;
+naturalMovie_output.memP_AVGMovie=memP_AVGMovie;
 end
 
 function [IndS,IndL]= getIndex(x,valConv,Points)
@@ -112,4 +157,15 @@ elseif isempty(IndL)|| length(IndL)<Points %if you are on the right border
         IndS=find(x<valConv,Points*2-length(IndL));
     end
 end
-end
+end
+
+% %make SVD of the filter
+% [sp,sigma,temp] = svd(STA_output.STA2D_zoom,'econ');
+% sp=sp(:,1);
+% imagesc(reshape(sp,size(STA_output.rangeY_zoom,2),size(STA_output.rangeX_zoom,2)));
+% temp=temp(:,1);
+% %linear filters
+% linearFiler1=normSTA(sp);
+% linearFiler2=normSTA(temp);
+% stim2pred=filter2(linearFiler1,movieSignal2D_zoom,'valid');
+% stim2pred=filter2(linearFiler2',stim2pred,'valid');