writeXmlCcepVisualDetPhase.m

%% This code is adapated from writeXmlCcepVisualDeterministic to:
% - generate phase scrambled instances of required images on the fly
% - reference the phase scrambled images using absolute paths instead of entity names
%
%    If this code is used in a publication, please cite the manuscript:
%    "H Huang, KN Kay, NM Gregg, G Ojeda Valencia, M In, C Kapeller, Y Shu, GA Worrell, KJ Miller, and D Hermes.
%    Single pulse electrical stimulation in white matter modulates iEEG visual responses in human early visual cortex. (Under Review)"
%
%    A preprint is available currently at doi: https://doi.org/10.1101/2025.05.05.652264.
%
%    The dataset corresponding to this code and manuscript is in BIDS format (version 1.10.0) on OpenNeuro (ds006485),
%    and it will be made publicly available upon manuscript acceptance.
%
%    Copyright (C) 2025 Harvey Huang
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <https://www.gnu.org/licenses/>.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%% Path, Define images to use. Check and configure this for each subject

mkdir('output/subjectParadigms');

% Subject and names of the two stimuli to use, without path info. Only indicate up to the nsd number here.
sub = 'ECoG'; % dummy subject
img1 = 'shared0369_nsd28069'; % right hand
img2 = 'shared0646_nsd47071'; % left hand
imgRestName = 'rest127';
coh = [0, 25, 50, 100];

% SEED - change this for the different runs of the experiment
seed = 1011; % 818, 1011, 218, 503, nan for ordered; replaced previous 1009 with 1011 to avoid consecutive blanks

% Other parameters that shouldn't need changing
isis = [nan, 0, 0.1, 0.2]; % nan is sham stim, 0 means simultaneous, other values indicate interval in s from electrical to visual.
nTrials = 6; % number of trials per ordinary condition
halfDurRest = 0.5; % half-duration of the rest interval (not accounting for jitter)
halfDurImg = 0.5; % half-duration of the image duration
maxLagFrames = 24; % maximum number of frames to add to end of rest as random jitter

assert(any(isis == 0), 'Error: Need to have simultaneous isi in order to program blank visual condition');

%% Prepare data

assert(length(coh) == 4, 'Experiment is designed for 4 phase-scrambled coherence levels');
if ~isnan(seed)
    rng(seed);
end

coh = sort(coh); % sort from low to high

% Contains necessary task info for all stim x noise x ISI conditions
conditions = [];
for kk = 1:4 % ISI
    
    % removed previous condition of having sham at half the number of trials, bc often only 2 or 3 runs can be completed
    
    for ii = 1:2 % image
        for jj = 1:4 % noise levels
            c = struct();
            c.isi = isis(kk);
            c.img = sprintf('Img%d', ii);
            c.coh = coh(jj);
            c.groupID = ii*4 + jj - 3; % 2, 3, 4, 5 for img1 and 6, 7, 8, 9 for img2

            % for 0 coherence trials, only need half as many because the 2 images are equivalent
            if c.coh == 0
                c.nTrials = ceil(nTrials/2);
            else
                c.nTrials = nTrials;
            end

            conditions = [conditions; c];
        end
    end

    % only do a stim at 0 condition, dont need to have sham or any other ISIs
    if isis(kk) ~= 0, continue; end

    % add blank stim-only trials
    c = struct();
    c.isi = isis(kk);
    c.img = 'Blank';
    c.coh = nan;
    c.groupID = 10; % last groupID
    c.nTrials = 2*nTrials; % double the number with stim only because only 1 ISI condition, to get reliable CCEPs
    conditions = [conditions; c];
    
end

% Create trials structure that repeats each condition for n Trials in deterministic order
trials = [];
for ii = 1:length(conditions)
    for jj = 1:conditions(ii).nTrials % repeated nTrials time
        c = conditions(ii);
        lagS = (randi(maxLagFrames + 1) - 1)/60; % inclusive between 0 and maxLagFrames, expressed in seconds
        c.lag = floor(lagS*1e3)/1e3; % floor to nearest millisecond

        % for blank, the same rest image will be used, so no need to configure ISI or image number
        if strcmp(c.img, 'Blank')
            c.id = '';
        else
            c.id = sprintf('isi%03dms-%d', 1e3*c.isi, jj); % unique IDs that match the filenames saved
        end

        trials = [trials; c];
    end
end
trials = rmfield(trials, 'nTrials');

% randomize order of trials
if ~isnan(seed)
    trials = trials(randperm(length(trials)));
end

%% Write xml file

% Open file, write header stuff
mkdir(fullfile('output/subjectParadigms', sub));
fid = fopen(sprintf('output/subjectParadigms/%s/ccepVisualPhase_seed%d.xml', sub, seed), 'w');
fprintf(fid, '<?xml version="1.0" encoding="utf-8"?>\n\n');

% Necessary stuff before main task area
fprintf(fid, '<Data xmlns="xsd"\n');
fprintf(fid, '\txmlns:ns="http://www.w3.org/2001/XMLSchema-instance"\n');
fprintf(fid, '\tns:schemaLocation="xsd C:\\_SVN\\ECoGmeToo\\trunk\\sources\\ParadigmTools\\src\\main\\Paradigm\\ParadigmSchema.xsd">\n\n');
fprintf(fid, '\t<Paradigm BaseFolder="mediaPhaseScrambled/" TaskOrder="Deterministic">\n\n');

% Pre-run rest
fprintf(fid, '\t<Task ns:type="SingleTask" ID="ST_PreRest" DurationSeconds="5">\n');
fprintf(fid, '\t\t<Stimulus ns:type="PictureStimulus" FileName="%s_OffNostim.png"/>\n', imgRestName);
fprintf(fid, '\t</Task>\n\n');

% Main task
fprintf(fid, '\t<!-- Main task area with all conditions in determinstic sequential order -->\n');
for ii = 1:length(trials)
%for ii = 100:110
    switch trials(ii).img
        case 'Img1'
            pathStr = sprintf('%s_coh%d', img1, trials(ii).coh); % pathStr contains image type and coherence, for backwards compatibility to contrast-adjusted imgs
        case 'Img2'
            pathStr = sprintf('%s_coh%d', img2, trials(ii).coh);
        case 'Blank'
            pathStr = imgRestName; % use rest image
        otherwise
            poopydoopy;
    end
        
    if isnan(trials(ii).isi) % sham tasks
        fprintf(fid, '\t<!-- Trial #%d: Sham, %s coherence %d -->\n', ii, trials(ii).img, trials(ii).coh);
        fprintf(fid, '\t<Task ns:type="MultiTask" ID="T%d_MT_Sham_%s_Coh%d" Label="ShamStim" TaskOrder="Deterministic" SampleSize="4">\n', ...
                  ii, trials(ii).img, trials(ii).coh);

        % Rest, diode on
        fprintf(fid, '\t\t<Task ns:type="SingleTask" ID="T%d_Sham_%s_Coh%d-RESTON" DurationSeconds="%0.03f" Group="1">\n', ii, trials(ii).img, trials(ii).coh, halfDurRest);
        fprintf(fid, '\t\t\t<Stimulus ns:type="PictureStimulus" FileName="%s_OnNostim.png"/>\n', imgRestName);
        fprintf(fid, '\t\t</Task>\n');

        % Rest, diode off
        fprintf(fid, '\t\t<Task ns:type="SingleTask" ID="T%d_Sham_%s_Coh%d-RESTOFF" DurationSeconds="%0.03f" Group="1">\n', ii, trials(ii).img, trials(ii).coh, halfDurRest + trials(ii).lag); % variably lagged
        fprintf(fid, '\t\t\t<Stimulus ns:type="PictureStimulus" FileName="%s_OffNostim.png"/>\n', imgRestName);
        fprintf(fid, '\t\t</Task>\n');

        % Image with diode on, no stim
        fprintf(fid, '\t\t<Task ns:type="SingleTask" ID="T%d_Sham_%s_Coh%d-IMGON" DurationSeconds="%0.03f" Group="%d">\n', ii, trials(ii).img, trials(ii).coh, halfDurImg, trials(ii).groupID);
        fprintf(fid, '\t\t\t<Stimulus ns:type="PictureStimulus" FileName="%s_%s_OnNostim.png"/>\n', pathStr, trials(ii).id);
        fprintf(fid, '\t\t</Task>\n');

        % Image with diode off, no stim
        fprintf(fid, '\t\t<Task ns:type="SingleTask" ID="T%d_Sham_%s_Coh%d-IMGOFF" DurationSeconds="%0.03f" Group="%d">\n', ii, trials(ii).img, trials(ii).coh, halfDurImg, trials(ii).groupID);
        fprintf(fid, '\t\t\t<Stimulus ns:type="PictureStimulus" FileName="%s_%s_OffNostim.png"/>\n', pathStr, trials(ii).id);
        fprintf(fid, '\t\t</Task>\n');

    elseif trials(ii).isi == 0 % instant stim tasks
        fprintf(fid, '\t<!-- Trial #%d: ISI=0, %s coherence %d -->\n', ii, trials(ii).img, trials(ii).coh);
        fprintf(fid, '\t<Task ns:type="MultiTask" ID="T%d_MT_Instant_%s_Coh%d" Label="InstantStim" TaskOrder="Deterministic" SampleSize="5">\n', ...
                  ii, trials(ii).img, trials(ii).coh);

        % Rest, diode on
        fprintf(fid, '\t\t<Task ns:type="SingleTask" ID="T%d_Instant_%s_Coh%d-RESTON" DurationSeconds="%0.03f" Group="1">\n', ii, trials(ii).img, trials(ii).coh, halfDurRest);
        fprintf(fid, '\t\t\t<Stimulus ns:type="PictureStimulus" FileName="%s_OnNostim.png"/>\n', imgRestName);
        fprintf(fid, '\t\t</Task>\n');

        % Rest, diode off
        fprintf(fid, '\t\t<Task ns:type="SingleTask" ID="T%d_Instant_%s_Coh%d-RESTOFF" DurationSeconds="%0.03f" Group="1">\n', ii, trials(ii).img, trials(ii).coh, halfDurRest + trials(ii).lag); % variably lagged
        fprintf(fid, '\t\t\t<Stimulus ns:type="PictureStimulus" FileName="%s_OffNostim.png"/>\n', imgRestName);
        fprintf(fid, '\t\t</Task>\n');

        % deal with double underscore problem
        if isempty(trials(ii).id) % no ID case, as in for Blank visual stimulus trials
            fnameStringCombined = pathStr;
        else
            fnameStringCombined = sprintf('%s_%s', pathStr, trials(ii).id);
        end

        % Image with diode on and electrical stimulation (0.1s)
        fprintf(fid, '\t\t<Task ns:type="SingleTask" ID="T%d_Instant_%s_Coh%d-IMGONSTIM" DurationSeconds="0.100" Group="%d">\n', ii, trials(ii).img, trials(ii).coh, trials(ii).groupID);
        fprintf(fid, '\t\t\t<Stimulus ns:type="PictureStimulus" FileName="%s_OnStim.png"/>\n', fnameStringCombined);
        fprintf(fid, '\t\t</Task>\n');

        % Image with diode on, no stim (halfDurImg - 0.1s)
        fprintf(fid, '\t\t<Task ns:type="SingleTask" ID="T%d_Instant_%s_Coh%d-IMGON" DurationSeconds="%0.03f" Group="%d">\n', ii, trials(ii).img, trials(ii).coh, halfDurImg - 0.1, trials(ii).groupID);
        fprintf(fid, '\t\t\t<Stimulus ns:type="PictureStimulus" FileName="%s_OnNostim.png"/>\n', fnameStringCombined);
        fprintf(fid, '\t\t</Task>\n');

        % Image with diode off, no stim
        fprintf(fid, '\t\t<Task ns:type="SingleTask" ID="T%d_Instant_%s_Coh%d-IMGOFF" DurationSeconds="%0.03f" Group="%d">\n', ii, trials(ii).img, trials(ii).coh, halfDurImg, trials(ii).groupID);
        fprintf(fid, '\t\t\t<Stimulus ns:type="PictureStimulus" FileName="%s_OffNostim.png"/>\n', fnameStringCombined);
        fprintf(fid, '\t\t</Task>\n');

    else % ISI > 0 tasks
        assert(trials(ii).isi > 0 && trials(ii).isi < halfDurRest, 'ISI must be positive and shorter than half the rest duration');
        fprintf(fid, '\t<!-- Trial #%d: ISI=%0.0fms, %s coherence %d -->\n', ii, trials(ii).isi*1e3, trials(ii).img, trials(ii).coh);
        fprintf(fid, '\t<Task ns:type="MultiTask" ID="T%d_MT_ISI%0.0fms_%s_Coh%d" Label="ISI%0.0fmsStim" TaskOrder="Deterministic" SampleSize="5">\n', ...
                  ii, trials(ii).isi*1e3, trials(ii).img, trials(ii).coh, trials(ii).isi*1e3);

        % Rest, diode on
        fprintf(fid, '\t\t<Task ns:type="SingleTask" ID="T%d_ISI%0.0fms_%s_Coh%d-RESTON" DurationSeconds="%0.03f" Group="1">\n', ii, trials(ii).isi*1e3, trials(ii).img, trials(ii).coh, halfDurRest);
        fprintf(fid, '\t\t\t<Stimulus ns:type="PictureStimulus" FileName="%s_OnNostim.png"/>\n', imgRestName);
        fprintf(fid, '\t\t</Task>\n');

        % Rest, diode off, halfDurRest minus ISI
        fprintf(fid, '\t\t<Task ns:type="SingleTask" ID="T%d_ISI%0.0fms_%s_Coh%d-RESTOFF" DurationSeconds="%0.03f" Group="1">\n', ii, trials(ii).isi*1e3, trials(ii).img, trials(ii).coh, halfDurRest - trials(ii).isi + trials(ii).lag); % variably lagged
        fprintf(fid, '\t\t\t<Stimulus ns:type="PictureStimulus" FileName="%s_OffNostim.png"/>\n', imgRestName);
        fprintf(fid, '\t\t</Task>\n');

        % Rest, diode off, stim for ISI duration
        fprintf(fid, '\t\t<Task ns:type="SingleTask" ID="T%d_ISI%0.0fms_%s_Coh%d-RESTOFFSTIM" DurationSeconds="%0.03f" Group="1">\n', ii, trials(ii).isi*1e3, trials(ii).img, trials(ii).coh, trials(ii).isi);
        fprintf(fid, '\t\t\t<Stimulus ns:type="PictureStimulus" FileName="%s_OffStim.png"/>\n', imgRestName);
        fprintf(fid, '\t\t</Task>\n');

        % Image with diode on, no stim
        fprintf(fid, '\t\t<Task ns:type="SingleTask" ID="T%d_ISI%0.0fms_%s_Coh%d-IMGON" DurationSeconds="%0.03f" Group="%d">\n', ii, trials(ii).isi*1e3, trials(ii).img, trials(ii).coh, halfDurImg, trials(ii).groupID);
        fprintf(fid, '\t\t\t<Stimulus ns:type="PictureStimulus" FileName="%s_%s_OnNostim.png"/>\n', pathStr, trials(ii).id);
        fprintf(fid, '\t\t</Task>\n');

        % Image with diode off, no stim
        fprintf(fid, '\t\t<Task ns:type="SingleTask" ID="T%d_ISI%0.0fms_%s_Coh%d-IMGOFF" DurationSeconds="%0.03f" Group="%d">\n', ii, trials(ii).isi*1e3, trials(ii).img, trials(ii).coh, halfDurImg, trials(ii).groupID);
        fprintf(fid, '\t\t\t<Stimulus ns:type="PictureStimulus" FileName="%s_%s_OffNostim.png"/>\n', pathStr, trials(ii).id);
        fprintf(fid, '\t\t</Task>\n');

    end

    fprintf(fid, '\t</Task>\n\n'); % close task for current condition
    
end

% Post-run rest
fprintf(fid, '\t<!-- End on a rest with photodiode on so as to time the end of the final visual stimulus -->\n');
fprintf(fid, '\t<Task ns:type="SingleTask" ID="ST_PostRestOn" DurationSeconds="1" Group="1">\n');
fprintf(fid, '\t\t<Stimulus ns:type="PictureStimulus" FileName="%s_OnNostim.png"/>\n', imgRestName);
fprintf(fid, '\t</Task>\n\n');

% End screen
fprintf(fid, '\t<Task ns:type="SingleTask" ID="ST_PostRestText" DurationSeconds="4">\n');
fprintf(fid, '\t\t<Stimulus ns:type="TextStimulus" Caption="Thank You!"/>\n');
fprintf(fid, '\t</Task>\n\n');

% End paradigm and close file
fprintf(fid, '\t</Paradigm>\n</Data>');
fclose(fid);

fprintf('Total duration of run = %0.1f s\n', length(trials)*2*(halfDurRest + halfDurImg) + sum([trials.lag]) + 10); % all trials + lag + 5s pre-rest + 5s post-rest
fprintf('Saved to output/subjectParadigms/%s/ccepVisualPhase_seed%d.xml\n', sub, seed);

return % to prevent execution of next section

%% Save all instances of phase-scrambled images to output

baseFolder = 'data/derivatives/visual_stimuli/originals'; % where the base images are found

shrinkage = 0.3; % shrink by 30% contrast

mkdir(sprintf('output/subjectParadigms/%s/mediaPhaseScrambled', sub));
addpath('path/to/knkutils/imageprocessing');
addpath('path/to/knkutils/math');

% Put in the rest files
imgBlank = 127*ones([1080, 1080, 3], 'uint8'); % color of each pixel = [127, 127, 127];

imgBlank = addFixation(imgBlank, 0.7, 0.5, 'r');
imgBlank = padToSize(imgBlank, [1080, 1920], [127, 127, 127]);

imgBlankOnStim = addPhotodiode(addPhotodiode(imgBlank, 70, 'w', 'bottom_right'), 70, 'w', [1765, 1011]); % 16 left of bottom right corners
imgBlankOnNostim = addPhotodiode(addPhotodiode(imgBlank, 70, 'w', 'bottom_right'), 70, 'k', [1765, 1011]);
imgBlankOffStim = addPhotodiode(addPhotodiode(imgBlank, 70, 'k', 'bottom_right'), 70, 'w', [1765, 1011]);
imgBlankOffNostim = addPhotodiode(addPhotodiode(imgBlank, 70, 'k', 'bottom_right'), 70, 'k', [1765, 1011]);

%imshow(imgBlank); axis equal
imwrite(imgBlankOnStim, sprintf('output/subjectParadigms/%s/mediaPhaseScrambled/%s_OnStim.png', sub, imgRestName));
imwrite(imgBlankOnNostim, sprintf('output/subjectParadigms/%s/mediaPhaseScrambled/%s_OnNostim.png', sub, imgRestName));
imwrite(imgBlankOffStim, sprintf('output/subjectParadigms/%s/mediaPhaseScrambled/%s_OffStim.png', sub, imgRestName));
imwrite(imgBlankOffNostim, sprintf('output/subjectParadigms/%s/mediaPhaseScrambled/%s_OffNostim.png', sub, imgRestName));

% get mean phase from both images
img1load = imread(fullfile(baseFolder, sprintf('%s.png', img1)));
img2load = imread(fullfile(baseFolder, sprintf('%s.png', img2)));
img1load = imresize(img1load, 1080/size(img1load, 1));
img2load = imresize(img2load, 1080/size(img2load, 1));
img1load = rgb2gray(img1load);
img2load = rgb2gray(img2load);
img1load = rgb2lin(img1load, 'ColorSpace', 'adobe-rgb-1998', 'OutputType', 'double');
img2load = rgb2lin(img2load, 'ColorSpace', 'adobe-rgb-1998', 'OutputType', 'double');

% reduce contrast of images to calculate for magmean
img1load = imadjustKeepMean(img1load, shrinkage);
img2load = imadjustKeepMean(img2load, shrinkage);

% calculate mean magnitude spectrum
spec1 = abs(fft2(img1load)); spec2 = abs(fft2(img2load));
magMean = sqrt(spec1.*spec2); % magnitude spectrum appears to be log-normal distributed. Therefore geometrically average
r = corrcoef(spec1(:), spec2(:));
sim = atanh(r(1, 2)); % fisher transform

for ii = 1:length(trials)
    
    fprintf('Saving image for trial %d: %s coherence %d, ID: %s\n', ii, trials(ii).img, trials(ii).coh, trials(ii).id);
    
    switch trials(ii).img
        case 'Img1'
            img = imread(fullfile(baseFolder, sprintf('%s.png', img1)));
            name = img1; % for writing
        case 'Img2'
            img = imread(fullfile(baseFolder, sprintf('%s.png', img2)));
            name = img2;
        case 'Blank'
            fprintf('Nothing saved\n');
            continue; % rest images already exist, no need to write anything
    end
    
    imgLarge = imresize(img, 1080/size(img, 1));
    
    % phase scramble, using mean magnitude spectrum. Using uniform noise to guarantee 0-coh cond are identical, adding contrast shrinkage to image
    imgOut = phaseScrambleGray(imgLarge, trials(ii).coh, false, magMean, shrinkage); 
    
    imgOut = addFixation(imgOut, 100*0.2/8.4*(323/1080), 0.5);
    imgOut = padToSize(imgOut, [1080, 1920], [127, 127, 127]);
    
    imgOutOnNostim = addPhotodiode(addPhotodiode(imgOut, 70, 'w', 'bottom_right'), 70, 'k', [1765, 1011]); % white image pd, black stim pd
    imgOutOffNostim = addPhotodiode(addPhotodiode(imgOut, 70, 'k', 'bottom_right'), 70, 'k', [1765, 1011]); % black image pd, black stim pd
    
    imwrite(imgOutOnNostim, sprintf('output/subjectParadigms/%s/mediaPhaseScrambled/%s_coh%d_%s_OnNostim.png', sub, name, trials(ii).coh, trials(ii).id));
    imwrite(imgOutOffNostim, sprintf('output/subjectParadigms/%s/mediaPhaseScrambled/%s_coh%d_%s_OffNostim.png', sub, name, trials(ii).coh, trials(ii).id));
    
    if trials(ii).isi == 0 % instant stim case requires the Onstim condition
        imgOutOnStim = addPhotodiode(addPhotodiode(imgOut, 70, 'w', 'bottom_right'), 70, 'w', [1765, 1011]); % white image pd, white stim pd
        imwrite(imgOutOnStim, sprintf('output/subjectParadigms/%s/mediaPhaseScrambled/%s_coh%d_%s_OnStim.png', sub, name, trials(ii).coh, trials(ii).id));
    end
    
end

fprintf('\nSimilarity of images 1 and 2 = %0.3f\n', sim);

% write metadata about shrinkage
fid = fopen(sprintf('output/subjectParadigms/%s/mediaPhaseScrambled/metadata.txt', sub), 'w');
fprintf(fid, 'img1\t%s\nimg2\t%s\nshrinkage\t%0.02f', img1, img2, shrinkage);
fclose(fid);