function Local_Chromatic_Integration_Stimulus(varargin)
%
%===================================================================================
% Parameter             Default                     Usage
%===================================================================================
%
% selecttextfile        false           Browse prompt to select the parameters text file
% stimduration          600             Stimulus presentation duration in refresh rate unit (600 = 10 sec).
% maxcontrast           0.2             Maximum Weber contrast value (ranges from -1 to 1)
% mincontrast           -0.2            Minimum Weber contrast value (ranges from -1 to 1)
% contrastdiff          0.02            Contrast steps from mincontrast to maxcontrast e.g. -20 : 2 : 20
% stixelwidth           20              Width of each local stimulus in pixels
% stixelheight          20              Height of each local stimulus in pixels
% gapwidth              2               Number of locations to be avoided around each locally selected stimulus area
% gapheight             2               Same as gapwidth but for the y-axis of each selected location
% seedlocation          -1000           Seed for random number generator for locations’ list (should be negative)
% seedonoroffpixels     -10000          Seed for random number generator of binary contrasts (not used here)
% seedcolorcontrast     -2000           Seed for random number generator for color contrast list (should be negative)
% redmeanintensity      0               Mean intensity for red gun of the screen
% greenmeanintensity	0.5             Mean intensity for green gun of the screen
% bluemeanintensity     0.5             Mean intensity for blue gun of the screen
% screensize            864 x 480       Screen resolution, default value is the resolution of the lightcrafter projector
% refreshrate           60              Screen refresh rate in Hz
% drawcircle            true            Option to display stimulus as circle or squares
% drawannulus           false           Option to display the stimulus as annulus (not available for this stimulus)
% fullscreen            false           Option to display the stimulus in full-screen mode
% help                  false           Option to check the list of available parameters
% lmargin               0               Margins from left side of the screen (not available for this stimulus)
% rmagin                0               Margins from right side of the screen (not available for this stimulus)
% tmargin               0               Margins from top of the screen (not available for this stimulus)
% bmargin               0               Margins from bottom of the screen (not available for this stimulus)
% coneisolating         false           Option to activate opsin-isolation (excluded for simplicity)
%
%===================================================================================


para            =   read_stimulus_parameters(varargin{:});
if para.help, return; end


offcontrasts    =   fliplr(0:-para.contrastdiff:para.mincontrast);
oncontrasts     =   0:para.contrastdiff:para.maxcontrast;

greencontrasts  =   [offcontrasts,oncontrasts];
bluecontrasts   =   [oncontrasts,offcontrasts];

% convert to weber contrast (mean+(contrast*mean))
greencontrasts  =   para.greenmeanintensity + (greencontrasts * para.greenmeanintensity);
bluecontrasts   =   para.bluemeanintensity + (bluecontrasts * para.bluemeanintensity);

% first get the corrdinates of all the locations
[f, v]          =   draw_locations(para);

% The darwing of the screen goes here
monitorsize =   get(0,'ScreenSize');
scpos       =   [monitorsize(3)/2-(para.screensize(1)/2), ...
    monitorsize(4)/2-(para.screensize(2)/2), para.screensize];
% make an screen like figure
fh          =   figure('Menu','none','ToolBar','none','Position',scpos,'Color',0.5*[1 1 1]);
fh.Name     =   'Example of Local Chromatic Integration Stimulus, from Khani and Gollisch (2021)';
fh.Colormap =   gray;
ah          =   axes('Units','Normalize','Position',[0 0 1 1]);
                axis(ah,[0, para.screensize(1),0, para.screensize(2)]);
                axis(ah,'off');

% inital variables
framecounter=   0;
seedlocation    =   para.seedlocation;      % to randomize locations
seedcolor       =   para.seedcolorcontrast; % to randomize contrast order
stimdur         =   ceil((para.stimduration ./ para.refreshrate)/ 0.5); % start by an estimate of one cycle

while ishandle(fh)
    
    if framecounter == 2,   tic;    end     % start the timer on the second frame
    frameMod        =   mod(framecounter,stimdur);
    
    if frameMod     ==  0
        [contrastlist, seedlocation, seedcolor] = select_sparse_locations(seedlocation, seedcolor, greencontrasts, bluecontrasts, para);
    end
    if framecounter ==  0
        p   =     patch('Faces',f,'Vertices',v,'FaceVertexcData',contrastlist,'Facecolor','Flat','edgecolor','none');
    else
        p.FaceVertexCData            =          contrastlist;
    end
     
    % this is to draw the frames relatively accurately.
    drawnow;
    java.lang.Thread.sleep(1/para.refreshrate*1e3);
    %pause(1/para.refreshrate);
    if framecounter == 2
       t            =  toc; % take the time of one trial and adjust the time based on that.
       stimdur      =  ceil((para.stimduration ./ para.refreshrate)/ t); 
    end
    
    framecounter    =   framecounter+1;
end

end

%--------------------------------------------------------------------------------------------------%
%----------                               sub-functions                                  ----------%
%--------------------------------------------------------------------------------------------------%

function varargout = fisher_Yates_shuffle_order(seed,inputVec,varargin)
%
%%% fisher_Yates_shuffle_order %%%
%
%
% This function generate psudorandom permution similar to randperm in MATLAB
% but works with psudorandom number generator ran1. It also gives back the
% most recent seed value to continue the permuation in case of repeated trials.
% note that the direction of permutation is along x-axix or for columns of
% MATALB not for the rows.
%
%
% ===============================Inputs====================================
%
%   seed : seed value for random number generation.
%   inputVec : input vector used for permutation.
%
%================================Output====================================
%
%   testOrder : vector of permuted indices for the inputVec.
%   newSeed : the recent seed that used in the ran1 function.
%   outputVec : the permuted input vector along x-axis
%
% Note that the permution algorithem is based on Fisher-Yates shuffle
% algorithem identical to what is used in the stimulus program.
% for more info check :
% https://en.wikipedia.org/wiki/Fisher%E2%80%93Yates_shuffle
%
% written by Mohammad, 01.02.2016

newSeed         =   seed;
testOrder       =   zeros(1,size(inputVec,2));
testOrder(1)    =   1;

for i   =   2:length(inputVec)-1
    
    [randVal,newSeed] = ran1(newSeed);
    
    j               =   ceil(i*randVal);    % based on Fischer-Yates algorithem
    testOrder(i)    =   testOrder(j);
    testOrder(j)    =   i;
    
end

testOrder           =   [testOrder(end),testOrder(1:end-1)]+1;   % to match MATLAB indexing
varargout{1}        =   testOrder;
varargout{2}        =   newSeed;
for j   =   1:size(inputVec,1)
    varargout{3}(j,:)   =   inputVec(j,testOrder);
end

end

%--------------------------------------------------------------------------------------------------%


function [f, v] = draw_locations(para)

Nx      =   ceil(para.screensize(1)/para.stixelwidth);
Ny      =   ceil(para.screensize(2)/para.stixelheight);

if para.drawcircle,     numpoints = 100;        else,    numpoints = 5;     end

locs    =   nan(Nx*Ny,numpoints,2);
iter    =   1;
for ii  =   1 : para.stixelheight : para.screensize(2)
    for jj  =   1 : para.stixelwidth : para.screensize(1)
        if para.drawcircle
            locs(iter,1:numpoints,1)    =   jj + ((para.stixelwidth/2).*cos(linspace(0, 2*pi,numpoints))) - para.stixelwidth/2;
            locs(iter,1:numpoints,2)    =   ii + ((para.stixelheight/2).*sin(linspace(0, 2*pi,numpoints))) - para.stixelheight/2;
        else % little hack to draw squares
            locs(iter,1:numpoints,1)    =   [jj, jj+para.stixelwidth, jj+para.stixelwidth, jj, jj]-para.stixelwidth;
            locs(iter,1:numpoints,2)    =   [ii, ii, ii+para.stixelheight, ii+para.stixelheight,ii]-para.stixelheight;
        end
        iter                            =   iter+1;
    end
end
x       =   squeeze(locs(:,:,1))';
y       =   squeeze(locs(:,:,2))';

f       =   reshape(1:Nx*Ny*numpoints,numpoints,Nx*Ny)';    % faces for patch fuction
v       =   [x(:),y(:)];                                    % vertices for patch function
%c      =   repmat(0.5*[1 1 1],Nx*Ny,1);                    % colors for patch function
end

%--------------------------------------------------------------------------------------------------%

function [stimMat, seedlocationout, seedcolorout] = select_sparse_locations(seedlocation, seedcolor, greencontrasts, bluecontrasts, para)

numcontrasts    =   size(greencontrasts,2);
Nx              =   ceil(para.screensize(1)/para.stixelwidth);
Ny              =   ceil(para.screensize(2)/para.stixelheight);

stimMat         =   zeros(Ny*Nx,3)+0.5;         % set all locations to gray
locationlist    =   (0:Nx*Ny-1)';
neighbourlist   =   get_neighbourlist(Nx, Ny, para.gapwidth, para.gapheight);

% make an estimate of all possible locations
possibleXYlocs  =   zeros(floor(numel(locationlist)/min(cellfun(@numel,neighbourlist))),1);

%outside loop is faster
[randvals, seedlocationout]   =   ran1(seedlocation, numel(possibleXYlocs));
randIdx         =   1;

availablelist   =   locationlist;
currXYlocs      =   possibleXYlocs;
idx             =   1;

while numel(availablelist)>0
    
    randval     =   randvals(randIdx);
    randIdx     =   randIdx + 1;
    stixelInd   =   floor(randval * size(availablelist,1)) + 1;
    toadd       =   availablelist(stixelInd);
    toremove    =   neighbourlist{Ny * mod(toadd,Nx) + floor(toadd/Nx)+1};
    %using the mex function for speed
    % availablelist(builtin('_ismemberoneoutput',availablelist,sort(toremove)))=[];
    % this is for older matlab, for new matlab use the line below
    availablelist(ismembc(availablelist,sort(toremove)))    =   [];
    
    currXYlocs(idx)     =   toadd+1;
    idx         =   idx + 1;
end

contlist        =   cell(1,ceil((idx-1)/length(1:numcontrasts)));
for jj  =    1:ceil((idx-1) / length(1 : numcontrasts))
    [contlist{jj}, seedcolor] = fisher_Yates_shuffle_order(seedcolor, 1:numcontrasts);
end
seedcolorout    =   seedcolor;
contindices     =   cell2mat(contlist)';
%contperframe = contindices(1:idx-1);
stimMat(currXYlocs(1:idx-1),1)      =   0;
stimMat(currXYlocs(1:idx-1),2)      =   greencontrasts (contindices(1:idx-1));
stimMat(currXYlocs(1:idx-1),3)      =   bluecontrasts (contindices(1:idx-1));

end

%--------------------------------------------------------------------------------------------------%

function neighbourlist = get_neighbourlist(nx, ny, gapx, gapy)

neighbourlist       =       cell(nx*ny,1);
idx                 =       1;
for x   =   0: nx-1
    for y   =   0: ny-1
        localneigbour       =       [];
        for currX           = x-gapx : x+gapx
            for currY       = y-gapy : y+gapy
                if (currX  >= 0 && currY >= 0 )
                    localneigbour   =   [localneigbour, nx * currY + currX];    %#ok
                end
            end
        end
        neighbourlist{idx}      =    localneigbour;
        idx                     =    idx+1;
    end
end

end

%--------------------------------------------------------------------------------------------------%

function paraout = read_stimulus_parameters(varargin)

% first parse the user inputs
p = inputParser();      % check the user options.
p.addParameter('selecttextfile', false,  @(x) islogical(x) || (isnumeric(x) && ismember(x,[0,1])));
p.addParameter('stimduration', 30, @isnumeric);
p.addParameter('contrastdiff', 0.02, @isnumeric);
p.addParameter('mincontrast',-0.2, @isnumeric);
p.addParameter('maxcontrast', 0.2, @isnumeric);

p.addParameter('stixelwidth', 20, @isnumeric);
p.addParameter('stixelheight', 20, @isnumeric);
p.addParameter('gapwidth', 2, @isnumeric);
p.addParameter('gapheight', 2, @isnumeric);

p.addParameter('seedlocation', -1000, @isnumeric);
p.addParameter('seedonoroffpixels', -10000, @isnumeric);
p.addParameter('seedcolorcontrast', -2000, @isnumeric);

p.addParameter('redmeanintensity', 0, @isnumeric);
p.addParameter('greenmeanintensity', 0.5, @isnumeric);
p.addParameter('bluemeanintensity', 0.5, @isnumeric);
p.addParameter('lmargin',0, @isnumeric);
p.addParameter('rmargin', 0, @isnumeric);
p.addParameter('bmargin', 0, @isnumeric);
p.addParameter('tmargin', 0, @isnumeric);
p.addParameter('coneisolating', true,  @(x) islogical(x) || (isnumeric(x) && ismember(x,[0,1])));
p.addParameter('drawcircle', true, @(x) islogical(x) || (isnumeric(x) && ismember(x,[0,1])));
p.addParameter('drawannulus', false, @(x) islogical(x) || (isnumeric(x) && ismember(x,[0,1])));
p.addParameter('screensize', [864 480], @isnumeric);
p.addParameter('refreshrate', 60, @isnumeric);
p.addParameter('fullscreen', false,  @(x) islogical(x) || (isnumeric(x) && ismember(x,[0,1])));
p.addParameter('help', false,  @(x) islogical(x) || (isnumeric(x) && ismember(x,[0,1])));

p.parse(varargin{:});
% defualt parameters
defpara     =    p.Results;

if defpara.help
    help_info_for_parameters(defpara);
    paraout     =   defpara;
    return;
end

if defpara.selecttextfile
    % now read the text files
    [stimfile, stimpath] = uigetfile('*.txt','Select a chromatic integration stimulus parameter file','chromatic_integration.txt');
    
    fid     =   fopen([stimpath,filesep,stimfile]);
    tline   =   fgetl(fid);
    while ischar(tline)
        tline   =   fgetl(fid);
        if tline == -1, break; end
        fn      =   extractBefore(tline,' = ');
        if isempty(fn), continue; end
        val     =   extractAfter(tline,' = ');
        % to convert to double
        if ~isnan(str2double(val))
            val =   str2double(val);
        end
        % to convert to logical
        if strcmp(val,'true'),        val = true;       end
        if strcmp(val,'false'),       val = false;       end
        % store the values in a structure
        stimpara.(fn)   =   val;
    end
    fclose(fid);
    
    % compare the text file to the defualt values and fill the missing parameters
    fn      =       fieldnames(defpara);
    for ii  =   1:numel(fn)
        if isfield(stimpara,fn{ii})
            paraout.(fn{ii})    =       stimpara.(fn{ii});
        else
            paraout.(fn{ii})    =       defpara.(fn{ii});
        end
    end
else
    paraout     =       defpara;
end

if paraout.fullscreen
    monitorsize         =   get(0,'ScreenSize');
    paraout.screensize  =   monitorsize(3:4) ;
end

end

%--------------------------------------------------------------------------------------------------%

function help_info_for_parameters(defpara)

fn          =   fieldnames(defpara);
fprintf(['\n\n',repmat('==',1,50),'\r\n']);
fprintf([repmat(' ',1,35),'List of Stimulus Parameters\r\n']);
fprintf([repmat('==',1,50),'\r\n']);
maxtxtlen   =   max(cellfun(@numel,fn))+10;
for ii      =  1:numel(fn)
    g       =  repmat(' ',1,maxtxtlen - numel(fn{ii}));
   fprintf(['\t-- \t%s',g,':',repmat(' ',1,10),'%s\n'],fn{ii},num2str(defpara.(fn{ii}))); 
end
fprintf(['\n',repmat('==',1,50),'\r\n']);

end