LD_phantom_oscill_percept_bis.m

function LD_phantom_oscill_percept_bis(subject, session, debug)

%In this version, we add multiple velocities
% subject = 'Dave';                                                                                                                                                                                                                                                     
% session = 1;                                                                                                                           
% debug = 1;


ex.version = 'vrandpercept';
global EyeData rect w xc yc eye_used 
%%%% resolution 
if debug == 1
    ex.screenWidth = 53.1;             % in cm; %laptop=27.5,office=43, %19=%T3b, miniHelm=39;
    ex.viewingDist = 53.5;             % in cm; 3Tb/office=43, miniHelm=57;
	ex.resolution = SetResolution(max(Screen('Screens')),2880, 1800 ,0); % laptop 1920,1080/ 2880, 1800 ,0
    ex.gammaCorrection = 0;       % make sure this = 1 when you're at the scanner!
else
                                                                                                                               
    ex.screenWidth = 53.1;             % in cm; % 16 in eye tracking room 425%laptop=27.5,office=43, %19=%T3b, miniHelm=39;
    ex.viewingDist = 53.5;             % in cm; %23 in eye tracking                                                                                                                          room 425 3Tb/office=43, miniHelm=57;
    ex.resolution = SetResolution(max(Screen('Screens')),1600,900,60); % ET room 1600,900,60
    ex.gammaCorrection = 1;       % make sure this = 1 when you're at the scanner!
end

%%%% keyboard
[keyboardIndices, productNames, ~] = GetKeyboardIndices ;
deviceNumber = keyboardIndices(1);
responseKeys = zeros(1,256);
responseKeys(KbName('Return'))=1; % button box 1
Screen('Preference', 'SkipSyncTests', 1);

% ex.scanNum = input('Scan number :');
ex.runNum = input('Run number :');

%%% basic naming set-up
ex.subject = subject;
ex.session = session;


%%%% set-up rand
ex.startTime = clock;
rng(sum(100*ex.startTime));
ex.rand = rng;

%%%% files and things
ex.root = pwd;
ex.date = datestr(now,30);


%%%% 2D sine wave grating inducers properties
ex.stim.spatialFreqDeg = 0.5/2; %0.286;                                          % cycles per degree of visual angle
ex.stim.contrast = 0.3 ;                                                 % in %, maybe??
ex.stim.orientation = [90]; %[90 180];                                                % in degrees
ex.stim.degFromFix = .6;                                              % in degrees of visual angle
ex.stim.gaborHDeg = 8;                                                  % in degrees of visual angle
ex.stim.gaborWDeg = 16;
%ex.stim.rectGaborWDeg = 8;
ex.stim.contrast = 0.15;%linspace(0.01,0.20,10);%[0.05, 0.10, 0.15];                                                 % in %, maybe?? %here the number of stimulus contrast levels is the number of different conditions
ex.stim.contrastMultiplicator = ex.stim.contrast/2;  % for sine wave 0.5 = 100% contrast, 0.2 = 40%
ex.stim.contrastOffset = .5; %.5 .5 0];                                  % for procedural gabor
%ex.stim.cycPerSec = [1.13*1/2,1.13*3/2]; % try multiple speeds
%ex.stim.motionRate = ex.stim.cycPerSec.*360;                                          % 1.13 cycles per second = 360 deg of phase *1.13 per sec
%ex.stim.cycles =[1, 3]; %number of cycles shifted per lap  (modified to half the number of cycles per lap)
ex.stim.cycPerSec = [1,1.4]; %drifting speed in cycles of grating per sec
ex.stim.cycles = [1 1]; %number of cycles per lap


%%%% sine wave grating timing (within block scale)
ex.initialFixation = 6;        % in seconds
ex.finalFixation = 2;          % in seconds
ex.stimDur = (ex.stim.cycles./ex.stim.cycPerSec)*2;        % in seconds. 1.77 sec refers to sine wave grating 1.77 = 2cycles/1.13cyc.sec-1 mutiplied by 2 for back and forth
%ex.stimsPerBlock = 1;      % number of back-and-forth laps of the stimulus drift
ex.blockLength =  16; %ceil(ex.stimDur*ex.stimsPerBlock);           % in seconds
ex.betweenBlocks = 2;          % in seconds
ex.flipsPerSec = 60;  % 60;         % number of phase changes we want from the visual stimulus, and thus the number of times we want to change visual stimulation on the screen
ex.flipWin = 1/ex.flipsPerSec;         % in seconds then actually in 1 sec the stimuli will change 12 times 
%e.numBlocks = 12;  % 6 for single and 6 for pair...

%%%% conditions & layout (across blocks scale)

ex.conds = {'RectMinbgSp8Vel3','RectMeanbgSp8Vel3'...%,
    ... %
    }; 
ex.numConds = length(ex.conds);
% with line of code below we will have 1 condition per block, randomized. we might need to change that
% later, to have the conditions randomized within each block
ex.repsPerRun = 1;              % repetitions of each condition per run
ex.numBlocks = ex.numConds*ex.repsPerRun;

ex.condShuffle = [];
for i =1:ex.repsPerRun
    ex.condShuffle = [ex.condShuffle, Shuffle([1:ex.numConds])];
end


ex.trialFlips = (0:ex.flipWin:ex.blockLength(1));%+ex.betweenBlocks
ex.trialFlips = ex.trialFlips(1:end-1);
%%%% fixation 
ex.fixSizeDeg =  .25;            % in degrees, the size of the biggest white dot in the fixation

%%%% screen
ex.backgroundColor = [127 127 127];%[108.3760 108.3760 108.3760];%;  % color based on minimum gating luminance 
ex.fontSize = 26;

%% %%%%%%%%%%%%%%%%%
   % timing model  %
   %%%%%%%%%%%%%%%%%
ex.onSecs = ones(1,ex.blockLength(1));
ex.longFormBlocks = Expand(ex.onSecs,ex.flipsPerSec,1); %1 when block, 0 when between block
ex.longFormFlicker = repmat(ones(1,1),1,length(ex.longFormBlocks)); %1 all the way to ensure flip at every time selected
length(ex.longFormBlocks)

% %% create the timing model of stimulus conditions for this particular run
% clear i
for i =1:ex.numConds
    conditions(i).name = ex.conds(i);
    conditions(i).startTimes = [];
end

%%
%%%%%%%%%%%%%%%
% open screen %
%%%%%%%%%%%%%%%
HideCursor;
Priority(9);
%Priority(0);
%%%% open screen
screen=max(Screen('Screens'));
if debug
    [w, rect]=Screen('OpenWindow',screen,ex.backgroundColor,[],[],[],[],[],kPsychNeed32BPCFloat);
    %[w, rect]=Screen('OpenWindow',screen,ex.backgroundColor,[100 100 600 400],[],[],[],[],kPsychNeed32BPCFloat); %might need to switch 900 and 600 by 1600 and 1200 for room 425
    
else    
    [w, rect]=Screen('OpenWindow',screen,ex.backgroundColor,[],[],[],[],[],kPsychNeed32BPCFloat); %might need to switch 900 and 600 by 1600 and 1200 for room 425
end
Screen(w, 'TextSize', ex.fontSize);
Screen('BlendFunction', w, GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

%gamma correction, file prepared for room 425
if ex.gammaCorrection 
  %load gamma correction file
  load('/Users/tonglab/Desktop/monitor_calibration/425dell_22-12-09/phase2_photometry_22-12-09.mat');
  Screen('LoadNormalizedGammaTable', w, inverseCLUT);
end
%%%% timing optimization
frameInt = Screen('GetFlipInterval',w);
slack = frameInt/2;
frameRate =  1/frameInt;%Screen('NominalFrameRate',w);
% %%% Still red dot for 1 sec before trial starts
% ex.stillDotPhase = zeros(1,ex.flipsPerSec*ex.trialFixation);

flipTimes = [0:frameInt*frameRate/ex.flipsPerSec:ex.blockLength(1)]; %multiply frameInt by 60/12 = 5 to flip the image every 5 frames
ex.stim.flipTimes = flipTimes(1:length(flipTimes)-1);
% ex.stim.tempPhase1 = nan(ex.numConds,ex.repsPerRun);
% ex.stim.tempPhase2 = nan(ex.numConds,ex.repsPerRun);
% ex.stim.oscillation1 = nan(ex.numConds,ex.repsPerRun,length(ex.stim.flipTimes));
% ex.stim.oscillation2 = nan(ex.numConds,ex.repsPerRun,length(ex.stim.flipTimes));
% ex.stim.phases = nan(ex.numConds,ex.repsPerRun,length(ex.stim.flipTimes));
% 
% clear c r
% for c =1:ex.numConds
%     for r = 1:ex.repsPerRun
%         ex.stim.tempPhase1(c,r) = pi;
%        % ex.stim.tempPhase2(c,r) = rand(1,1)*2*pi;
%         ex.stim.oscillation1(c,r,:) = cos(2*pi*(1/ex.stimDur(1))*ex.stim.flipTimes+ex.stim.tempPhase1(c,r));
%         %ex.stim.oscillation2(c,r,:) = cos(2*pi*(1/ex.stimDur(2))*ex.stim.flipTimes+ex.stim.tempPhase2(c,r));
%         ex.stim.spatialPhase = 90;
%        % ex.stim.phases(c,r,:) = (ex.stim.oscillation1(c,r,:).*180*ex.stim.cycles(1)+ ex.stim.oscillation2(c,r,:).*180*ex.stim.cycles(2))/2 + ex.stim.spatialPhase; %./ex.stimDur-2*pi*flipTimes./ex.stimDur make it oscillatory
%     end
% end

ex.stim.tempPhase1 = pi;
ex.stim.spatialPhase = 90;
ex.stim.oscillation1 = cos(2*pi*(1/ex.stimDur(1))*ex.stim.flipTimes+ex.stim.tempPhase1);
ex.stim.phases = ex.stim.oscillation1.*180*ex.stim.cycles(1) + ex.stim.spatialPhase;
%reasonning behind the calculation of ex.stim.phases:
%GOAL: render the back and forth of grating drifts oscillatory in time instead
%of linear to smooth the signal phase shifts at the time of drift direction
%reversal
%1) The length (x) of grating shift after each flip is the variable that will have to be oscillatory
%over time

%2)oscillatory signal formula is x(t) = A*cos(w*t+phase) (x is the oscillatory signal, w is the angular
%frequency in rad/s, t is the time in s, phase is a phase shift in rad at time t (here 0),
%A is the amplitude of the signal

%3)Here the angular frequency is calculated based on the idea of one
%oscillation per stimulus. Considering a stimulus as 1 cycle forward and 1 cycle backwards(definition might change based on how many cycles are desired per stimulus), one oscillation will represent 2 cycles =ex.stim.cycPerSec*ex.stimDur)
% thus the temporal frequency of the oscillation should be 1/ex.stimDur
% (this is w)

%%4) We multiply w by t. Here the resolution of the signal is given by the number of flips per second, so the duration is given by flipTimes.
%Finally, we multiply by 2*pi for the units to be in radian.

%5)Then we subtract the phase at each flip at time t which here is 0.

%6)The amplitude of cos() varies between -1 and +1. Here we want our
%displacement (grating shift) in degrees with a displacement of 360? (of the spatial grating) every half temporal oscillation so we multiply by 180 for the range to be [-180;+180].

%7) we multiply by the number of cycles desired to drift over one lap
%(ex.stim.cycles).

%%%% scale the stim params for the screen
ex.ppd = pi* rect(3) / (atan(ex.screenWidth/ex.viewingDist/2)) / 360;
ex.fixSize = round(ex.fixSizeDeg*ex.ppd);
ex.gaborHeight = round(ex.stim.gaborHDeg*ex.ppd);                 % in pixels, the size of our objects
ex.gaborWidth = round(ex.stim.gaborWDeg*ex.ppd);                 % in pixels, the size of our objects 
ex.rawGaborHeight = ex.gaborHeight*3;
ex.rawGaborWidth = ex.gaborWidth*2;

%% Create only one big sinewave grating image saved for each repetition and each condition

ex.rectSWave = nan(ex.numConds, ex.repsPerRun,ex.rawGaborHeight,ex.rawGaborWidth);
ex.rectSWaveID = nan(ex.numConds, ex.repsPerRun);
clear c r
for c =1:ex.numConds 
    for r = 1:ex.repsPerRun %only save the first image of each trial, that we will move during the trial
        phase = ex.stim.phases(1);
        ex.rectSWave(c,r,:,:) = makeSineGrating(ex.rawGaborHeight,ex.rawGaborWidth,ex.stim.spatialFreqDeg,...
            ex.stim.orientation,phase,ex.stim.contrastOffset(1),ex.stim.contrastMultiplicator,...
            ex.ppd);
        ex.rectSWaveID(c,r) = Screen('MakeTexture', w, squeeze(ex.rectSWave(c,r,:,:)));
    end
end
%% Sine wave gratings locations (in the task loop since it changes)
xc = rect(3)/2; % rect and center, with the flexibility to resize & shift center - change vars to zero if not used.
yc = rect(4)/2; %+e.vertOffset;

%% create drifting red dots position
clear flipTimes
flipTimes = [0:frameInt*frameRate/ex.flipsPerSec:ex.blockLength(1)]; %multiply frameInt by 60/12 = 5 to flip the image every 5 frames
flipTimes = flipTimes(1:length(flipTimes)-1);
ex.stimDriftPosDeg = nan(ex.numConds,ex.repsPerRun,length(ex.stim.oscillation1(1,:)));
ex.stimDriftPos = nan(ex.numConds,ex.repsPerRun,length(ex.stim.oscillation1(1,:)));
ex.stimLongDriftPos = nan(ex.numConds,ex.repsPerRun,length(flipTimes));
clear c r
for c = 1:ex.numConds 
    for r = 1:ex.repsPerRun
%         ex.stimDriftPosDeg(c,r,:) = (ex.stim.oscillation1(c,r,:).*ex.stim.cycles(1).*1/(2*ex.stim.spatialFreqDeg)+ ex.stim.oscillation2(c,r,:).*ex.stim.cycles(2).*1/(2*ex.stim.spatialFreqDeg))/2;
        ex.stimDriftPosDeg(c,r,:) = ex.stim.oscillation1.*ex.stim.cycles(1).*1/(2*ex.stim.spatialFreqDeg);
        ex.stimFixSpatialPhase = 0; %-(1/(8*ex.stim.spatialFreqDeg))*ex.ppd;%(1/(4*ex.stim.spatialFreqDeg))*ex.ppd;
        ex.stimDriftPos(c,r,:) = ex.stimDriftPosDeg(c,r,:).*ex.ppd +ex.stimFixSpatialPhase;
%         ex.stimStillDotPhase = ex.stimDriftPos(c,r,1);
        ex.stimLongDriftPos(c,r,:) = [squeeze(ex.stimDriftPos(c,r,:))']; %drift for 1.25 (5/4) periods of oscilation
    end
end



% figure();
% subplot(4,1,1)
% plot(squeeze(ex.stimDriftPosDeg(1,1,:)))
% subplot(4,1,2)
% plot(squeeze(ex.stimDriftPosDeg(1,2,:)))
% subplot(4,1,3)
% plot(squeeze(ex.stimDriftPosDeg(1,3,:)))
% subplot(4,1,4)
% plot(squeeze(ex.stimDriftPosDeg(1,4,:)))
% ylabel('Stimulus position (dva)')
% xlabel('Flip #')


%% Create rectangular masks for the gratings
gray1 = repmat(mean(squeeze(ex.rectSWave(1,1,1,:))), [1,3]);
gray2 = repmat(min(min(squeeze(ex.rectSWave(1,1,:,:)),[],1)), [1,3]);

%phantom control condition
coaperture=Screen('OpenOffscreenwindow', w, gray1);
Screen('FillRect',coaperture, [255 255 255 0], [xc-(1/2)*ex.gaborWidth yc-(3/2)*ex.gaborHeight xc+ex.gaborWidth/2 yc-ex.gaborHeight/2]);
Screen('FillRect',coaperture, [255 255 255 0], [xc-(1/2)*ex.gaborWidth yc+ex.gaborHeight/2 xc+ex.gaborWidth/2 yc+(3/2)*ex.gaborHeight]);

%phantom condition
phaperture=Screen('OpenOffscreenwindow', w, gray2);
Screen('FillRect',phaperture, [255 255 255 0], [xc-(1/2)*ex.gaborWidth yc-(3/2)*ex.gaborHeight xc+ex.gaborWidth/2 yc-ex.gaborHeight/2]);
Screen('FillRect',phaperture, [255 255 255 0], [xc-(1/2)*ex.gaborWidth yc+ex.gaborHeight/2 xc+ex.gaborWidth/2 yc+(3/2)*ex.gaborHeight]);


%% %%%% initial window - wait for backtick
DrawFormattedText(w,'Can you see the visual phantom within the gap at the center of the screen? \n\n Press Space to start'... % :  '...
    ,'center', 'center',[0 0 0]);
Screen(w, 'Flip', 0);
%WaitSecs(2);
KbTriggerWait(KbName('Space'), deviceNumber);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                         experiment                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%% START task TASK/FLIPPING
KbQueueCreate(deviceNumber,responseKeys);
n = 1;
blockCnt = 1;
cnt = 0;
cMin = 1;
cMean = 2;
onOffs = [diff(ex.longFormBlocks) 0];
ex.flipTime = nan(length(ex.trialFlips),length(ex.condShuffle));
%%% initial fixation
if n == 1 && blockCnt == 1 %for first block
    ex.tasktstart = clock;
    ex.startRun = GetSecs();
    Screen('FillRect', w, gray1);
    Screen(w, 'Flip', 0);
    WaitSecs(ex.initialFixation);
end
%%% Launch the task
c = 1;
%flip through the block 
while(1) % n <= length(ex.trialFlips)
    KbQueueStart();
    ex.longFormBlocks(n)
    if n == 1 && cnt == 0
        thisCond = ex.condShuffle(c);
        condName = conditions(thisCond).name{:};
    end
    %%%% draw sine wave grating stimulus %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if contains(condName, 'Minbg')
        Screen('FillRect', w, gray2);
        
        xOffset = ex.stimLongDriftPos(cMin,1,n)-ex.stimLongDriftPos(cMin,1,1); %baseline correct the position since every image already hase a spatial phase shift in the sinewave
        ex.rectLRect =  CenterRectOnPoint([0 0 ex.rawGaborWidth ex.rawGaborHeight],xc+xOffset,yc);
        % stim
        Screen('DrawTexture', w, ex.rectSWaveID(cMin,1),[],ex.rectLRect);
        
        Screen('DrawTexture',w,phaperture);
        
    elseif contains(condName, 'RectMeanbg')
        Screen('FillRect', w, gray1);
        xOffset = ex.stimLongDriftPos(cMean,1,n)-ex.stimLongDriftPos(cMean,1,1); %baseline correct the position since every image already hase a spatial phase shift in the sinewave
        ex.rectLRect =  CenterRectOnPoint([0 0 ex.rawGaborWidth ex.rawGaborHeight],xc+xOffset,yc);
        % stim
        Screen('DrawTexture', w, ex.rectSWaveID(cMean,1),[],ex.rectLRect);
        
        Screen('DrawTexture',w,coaperture);
    end
    
    %%%%%%%%%%% FLIP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if n == 1
        [VBLT, ex.startTrial, FlipT, missed] = Screen(w, 'Flip', 0);%[VBLTimestamp StimulusOnsetTime FlipTimestamp Missed] = Screen('Flip', windowPtr [, when] [, dontclear]...
        flipTimes = ex.startTrial;
        
    else
        [VBLT,flipTime, FlipT, missed] = Screen(w, 'Flip',ex.startTrial + ex.trialFlips(n) - slack); %,   %%% ex.flipTime(n,c)
        flipTimes = [flipTimes, flipTime];
    end
    KbQueueStop();
    [pressed, firstPress]= KbQueueCheck();
    if (pressed && ismember(find(firstPress,1), [KbName('Return') KbName('ENTER')]))
        n = 0;
        c = c+1;
        cnt = 0;
    end
    n = n+1;
    if n >= length(ex.longFormBlocks)
        n = 1;
        cnt = cnt+1;
    end

    KbQueueFlush();
    if c > ex.numBlocks
        break;
    end
end

%     if n == 1382%360%421%420%359%1382%1561%
%         break;
%     end

%%%%%%%%%%%%%%%%%%
% done! wrap up  %
%%%%%%%%%%%%%%%%%%

ex.runTime = GetSecs - ex.startRun;


savedir = fullfile(ex.root,'data',sprintf('s%s_%s/',subject,ex.version));
if ~exist(savedir); mkdir(savedir); end
savename = fullfile(savedir, strcat(sprintf('/s%s_smooth_pursuit_%s_date%s_fix',subject,ex.version,num2str(ex.date)), '.mat'));
%save(savename,'ex');
save(savename,'ex','-v7.3')

ShowCursor;
Screen('Close');
Screen('CloseAll');
fclose all;