clear all
close all

info = niftiinfo('sub-1_task-objectviewing_run-01_bold.nii');
data = niftiread(info);

figure
for k=1:64
    imagesc(squeeze(data(:,:,k,60))), axis image, colormap("gray")
    title(k)
    drawnow
    pause(0.1)
end

data_slice=double(squeeze(data(:,:,36,:)));

figure, imagesc(data_slice(:,:,1),[0 2000]), axis image, colormap(gray)
figure, imagesc(sum(data_slice,3)), axis image

mask=sum(data_slice,3)>1e+5;
figure, imagesc(mask), axis image, colormap("gray")

figure
for k=1:121
    imagesc(squeeze(data_slice(:,:,k)).*mask), axis image, colormap("gray")
    title(k)
    drawnow
    pause(0.2)
end

figure
for k=2:121
    imagesc(abs(squeeze(data_slice(:,:,k)-data_slice(:,:,k-1))).*mask,[0 50])
    axis image
    title(k)
    drawnow
    pause(0.2)
end

v=var(data_slice,0,3);
figure, imagesc(v.*mask,[0 4000]), axis image, colormap("hot")

figure
for k=1:111
    imagesc(var(data_slice(:,:,k:(k+10)),0,3).*mask,[0 1000]), colormap("hot")
    axis image
    title(k)
    drawnow
    pause(0.2)
end

t=0:info.PixelDimensions(4):(info.PixelDimensions(4)*120);
t=(0:120)*info.PixelDimensions(4);
s1=squeeze(data_slice(31, 30, :));
s2=squeeze(data_slice(28, 21, :));
s3=squeeze(data_slice(33, 27, :));
figure, plot(t,s1-mean(s1)), grid, hold on, plot(t,s2-mean(s2))
plot(t, s3-mean(s3))
legend('(31,30)','(28,21)','(33,27)')


% CORRELAZIONE
data_sl_decimate=zeros(20, 32, 121);
for k=1:size(data_sl_decimate,3)
    data_sl_decimate(:,:,k)=transpose(downsample(transpose(downsample(squeeze(data_slice(:,:,k)),2)),2));
end

figure, imagesc(sum(data_sl_decimate,3))
figure, imagesc(sum(data_sl_decimate,3)>1.2e+5)
mask=sum(data_sl_decimate,3)>1.2e+5;
N=sum(mask(:));

rho=zeros(N);
for k1=1:(N-1)
    a=squeeze(data_sl_decimate(x(k1),y(k1), :));
    for k2=(k1+1):(N)
        b=squeeze(data_sl_decimate(x(k2),y(k2), :));
        rho(k1, k2)=corr(a,b);
    end
end
figure, imagesc(rho)


% rho=zeros(size(data_sl_decimate,1).*size(data_sl_decimate,2));
% for k1=1:(size(rho,1)-1)
%     a=squeeze(data_sl_decimate(rem(k1, size(data_sl_decimate,1))+1, floor(k1/size(data_sl_decimate,1))+1, :));
%     for k2=(k1+1):(size(rho,1)-1)
%         b=squeeze(data_sl_decimate(rem(k2, size(data_sl_decimate,1))+1, floor(k2/size(data_sl_decimate,1))+1, :));
%         rho(k1, k2)=corr(a,b);
%     end
% end