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:\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