evaluateContrastShrinkage.m

%% This script tests how much contrast shrinking to perform on the special 100 images to minimize the effect of negative pixel intensities
%
%    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/>.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
addpath('path/to/knkutils/imageprocessing');
addpath('path/to/knkutils/math');

% the 2 images to mix
baseFolder = 'stims/special100'; % where the base images are found
imgList = glob(fullfile(baseFolder, sprintf('*.png')));
listSz = length(imgList);

%% Phase scramble pairs of images and calculate pixels lost by shrinkage

rng('default');

coh = 25;
npairs = 50; % how many pairs to test at each shrinkage

shrinkage = 0:0.05:0.8; % percent of shrinkage at each iteration
pixelsLost = zeros(2*npairs, length(shrinkage));

% all possible image pairs
[X, Y] = meshgrid(1:listSz, 1:listSz);
idxes = logical(triu(ones(listSz), 1)); % indices above diagonal
pairsAll = [Y(idxes), X(idxes)]; % all possible pairs of NSD images

for aa = 1:length(shrinkage)
    
    fprintf('%0.02f', shrinkage(aa));
    
    % select 100 random non-self non-duplicate pairs of NSDspecial images
    p = randperm(length(pairsAll));
    pairs = pairsAll(p(1:npairs), :); % first 100 randomly permuted pairs

    for kk = 1:npairs
        fprintf('.');

            imgPaths = [imgList(pairs(kk, 1)), imgList(pairs(kk, 2))]; % the 2 images to average together

            imgs = cell(2, 1);
            imgsGrayLin = cell(2, 1);

            % get mean phase from both images
            for ii = 1:2
                imgs{ii} = imread(imgPaths{ii});
                imgs{ii} = imresize(imgs{ii}, 1080/size(imgs{ii}, 1));

                imgsGrayLin{ii} = rgb2gray(imgs{ii});
                imgsGrayLin{ii} = rgb2lin(imgsGrayLin{ii}, 'ColorSpace', 'adobe-rgb-1998', 'OutputType', 'double');
                imgsGrayLin{ii} = imadjustKeepMean(imgsGrayLin{ii}, shrinkage(aa));
            end

            magMean = sqrt(abs(fft2(imgsGrayLin{1})) .* abs(fft2(imgsGrayLin{2})));

            imgsScramb = cell(2, 1);
            for ii = 1:2
                imgsScramb{ii} = phaseScrambleGrayNotRgb(imgsGrayLin{ii}, coh, false, magMean);
                %pixelsLost(kk*2-2+ii, aa) = sum(imgsScramb{ii} < 0, 'all') / numel(imgsScramb{ii});
                pixelsLost(kk*2-2+ii, aa) = sum(imgsScramb{ii} < 0 | imgsScramb{ii} > 1, 'all') / numel(imgsScramb{ii});
            end

    end
    fprintf('\n');
end


%% Boxplot across shrinkages for current phase coherence, save output

mkdir(fullfile('output', 'contrastShrinkage'));

figure;
boxplot(pixelsLost, 'Positions', shrinkage, 'Colors', 'b', 'PlotStyle', 'compact');
yline(0, 'Color', [0.5, 0.5, 0.5]);
xlim([min(shrinkage)-0.03, max(shrinkage)+0.03]);
set(gca, 'xtick', shrinkage, 'xticklabels', shrinkage, 'box', 'off');
xlabel('Contrast shrinkage')
ylabel('Fraction of Pixels Lost');
title(sprintf('Phase Coherence %d%%', coh));
saveas(gcf, fullfile('output', 'contrastShrinkage', sprintf('pixelsLostCoh%d', coh)), 'png');


%% Examples of images scrambled with different shrinkage

rng('default');

sh = 0.3; % shrinkage

mkdir(fullfile('output', 'contrastShrinkage', sprintf('examples_sh=%0.2f', sh)));

cohs = [0, 25, 50, 100];

% select 5 random image pairs
p = randperm(length(pairsAll));
pairs = pairsAll(p(1:5), :);

for kk = 1:size(pairs, 1)
        imgPaths = [imgList(pairs(kk, 1)), imgList(pairs(kk, 2))]; % the 2 images to average together

        imgs = cell(2, 1);
        imgsGrayLin = cell(2, 1);

        % get mean phase from both images
        for ii = 1:2
            imgs{ii} = imread(imgPaths{ii});
            imgs{ii} = imresize(imgs{ii}, 1080/size(imgs{ii}, 1));
            
            imgName = regexp(imgPaths{ii}, 'shared\d{4}', 'match');
            imwrite(imgs{ii}, fullfile('output', 'contrastShrinkage', sprintf('examples_sh=%0.2f', sh), sprintf('%s_orig.png', imgName{1})));

            imgsGrayLin{ii} = rgb2gray(imgs{ii});
            imgsGrayLin{ii} = rgb2lin(imgsGrayLin{ii}, 'ColorSpace', 'adobe-rgb-1998', 'OutputType', 'double');
            imgsGrayLin{ii} = imadjustKeepMean(imgsGrayLin{ii}, sh);
        end

        magMean = sqrt(abs(fft2(imgsGrayLin{1})) .* abs(fft2(imgsGrayLin{2})));

        for cc = 1:length(cohs)
            
            for ii = 1:2
                imgScramb = phaseScrambleGrayNotRgb(imgsGrayLin{ii}, cohs(cc), false, magMean);
                pixLost = sum(imgScramb < 0 | imgScramb > 1, 'all') / numel(imgScramb);
                imgOut = lin2rgb(imgScramb, 'ColorSpace', 'adobe-rgb-1998', 'OutputType', 'uint8');
                
                imgName = regexp(imgPaths{ii}, 'shared\d{4}', 'match');
                imwrite(imgOut, fullfile('output', 'contrastShrinkage', sprintf('examples_sh=%0.2f', sh), sprintf('%s_coh=%d_pixLost=%0.2f.png', imgName{1}, cohs(cc), pixLost)));
    
            end
            
        end

end


%% Functions

function imout = imadjustKeepMean(im, shrinkage)
    assert(shrinkage >= 0 && shrinkage < 1, 'Shrinkage must be on the range [0, 1)');
    assert(ismatrix(im), 'img should be 2D (grayscale) between 0 and 1');
    
    inmean = mean(im, 'all');

    outlims = [inmean*shrinkage, 1 - (1-inmean)*shrinkage]; % shrink proportionally toward mean, preserving mean
    
    imout = imadjust(im, [0, 1], outlims);
    %imout = imout - mean(imout, 'all') + mean(im, 'all'); % keep the same mean as input
    
end

% Version of phaseScrambleGray that doesn't convert to rgb at the end so its pixel distribution can be examined
function [imgScramb, phaseChange] = phaseScrambleGrayNotRgb(imgLin, coh, perm, inputMag)

    % request linearized grayscale input

    if nargin < 3 || isempty(perm), perm = true; end % don't use permuted so that both images are scrambled with the same uniform distribution
    
    %assert(isa(img, 'uint8'), 'Expects image in uint8');
    assert(coh >= 0 & coh <= 100, 'pct must be given as percent of contrast within [0, 100]');
    
    %imgGray = rgb2gray(img); % convert to grayscale
    %imgLin = rgb2lin(imgGray, 'ColorSpace', 'adobe-rgb-1998', 'OutputType', 'double'); % undo gamma to yield linear rgb (so that pixel intensity scales roughly with intensity)
    
    % FFT on linearized grayscale img
    imgLinFft = fft2(imgLin);
    if nargin < 4 || isempty(inputMag)
        mag = abs(imgLinFft);   % magnitude (from image or from input)
    else
        mag = inputMag;
    end
    ph = angle(imgLinFft);  % phase

    % Get random phase
    if perm
        ph2 = angle(generatepermutedphase(ph)); % phase permuted from input image
    else
        ph2 = angle(generaterandomphaseFromInput(ph)); % phase randomly reassigned from uniform distribution, but dc and nyquist coefficients kept the same
        %ph2 = angle(generaterandomphase(length(ph)));
        %ph2(1) = ph(1);
    end                         
    
    phdiff = circulardiff(ph2, ph, 2*pi);         % angle needed to reach random phase

    % Apply mixture between real and phase scrambled image, ifft
    imgScramb = real(ifft2(mag .* exp(1i * (ph + (1 - coh/100)*phdiff))));
    
    phFinal = angle(fft2(imgScramb)); % (redundantly) obtain final phase distribution
    phaseChange = mean(abs(sort(phFinal(:)) - sort(ph(:)))); % sort initial and final phase distributions, calculate average difference between distributions
    %phaseChange = mean(abs(circulardiff(sort(phFinal(:)), sort(ph(:)), 2*pi))); % using circulardiff? probably unnecessary cuz phases r sorted
    
    %imgOut = lin2rgb(imgScramb, 'ColorSpace', 'adobe-rgb-1998', 'OutputType', 'uint8'); % gamma correct grayscale image, convert back to uint8
    
end