% EE684_fmri.m  - Program to Read and Manipulate fMRI data
%
% x_dim = x-dimension of brain
% y_dim = y-dimension of brain
% z_dim = number of brain slices
% t_dim = t dimension of brain
% filename = file from which brain time courses are read
% btc(x,y,s,:) = time course for pixel(x,y) in slice s
%
% 'steve' data set consists of 31 epochs with 64 images in each epoch
% (64x64x1x1984). 
%
% Activation can be found at:
%   tcx= [34,35,30,33,32,33,31,29,35,32,10]
%   tcy= [50,50,48,50,50,49,48,48,15,32,10]
%
% Credit goes to Stephen Laconte for collecting and sharing this data set


clear all;
close all;

% Read anatomical image
anat_filename= 'steve_anat.sdt';
x_dima= 256;                  % brain fMRI image dimension in x
y_dima= 256;                  % brain fMRI image dimension in y
z_dim= 1;                   % number of brain fMRI slices
t_dim= 1;                % number of time points collected     

fid = fopen(anat_filename,'rb','b');
anat_sdt = fread(fid,'float32');
fclose(fid);
im_anat= reshape(anat_sdt,x_dima, y_dima)';
figure(1);
imagesc(im_anat);
colormap(gray);

% Read functional MRI data
fmri_filename= 'steve_fct.sdt';
x_dim= 64;                  % brain fMRI image dimension in x
y_dim= 64;                  % brain fMRI image dimension in y
z_dim= 1;                   % number of brain fMRI slices
t_dim= 1984;                % number of time points collected     

fid = fopen(fmri_filename,'rb','b');
brain_sdt = fread(fid,'float32');
fclose(fid);

btc = zeros(x_dim,y_dim,z_dim,t_dim);

index=1; 
for t=1:t_dim,
    for slice=1:z_dim,
        for x=1:x_dim,
            row = (y_dim*(index-1)+1):(y_dim*(index));
            btc(x,:,slice,t) = brain_sdt(row);
            index = index + 1;
        end;
    end;
end;

% animate fMRI images in data set
figure(2)
implay(squeeze(btc),10);


% These are the coordinates of sample points to test
tcy= [34,35,30,33,32,33,31,29,35,32,10];
tcx= [50,50,48,50,50,49,48,48,15,32,10];


Nt= 1024;               % number of points in timecourse to consider
t= (1:Nt);
ind= 1;
sig= squeeze(btc(tcx(ind)+1,tcy(ind)+1,1,t));
dsig= SpectralSubtraction(sig, 0.00002, 1);     % denoise signal

figure(3)
%plot( t, sig, 'b-', t, dsig, 'r-')
 subplot(4,1,1), plot(t, sig)
 subplot(4,1,2), plot(t, dsig)
 subplot(4,1,3), plot( t, sig, 'b-', t, dsig, 'r-')
 subplot(4,1,4), plot(t, sig-dsig)


% Detection functional activation by cross-correlation 
fmri_xcor= zeros(x_dim,y_dim);      % array to contain correlation values
ref= ones(64,1);
ref(1:32)= -1;
sref= repmat(ref,Nt/64);            % reference signal

for x=1:x_dim,
for y=1:x_dim,
    x1= squeeze(btc(x,y,1,:));
    %x1= SpectralSubtraction(x1, 0.00002, 1);     % denoise signal
    fmri_xcor(x,y)= max(abs(xcorr(x1, ref, 32)));
%     figure(4)
%     plot(sref);
%     cor1= xcorr(x1, ref, 64);
%     plot(cor1);
%     [Y,I]= max(cor1);
%     plot(I:Nt,dsig(I:end),'b', I:Nt,3500*sref(1:Nt-I+1),'r');
%  
end
end

figure(5);
activ_level= prctile (reshape(abs(fmri_xcor),x_dim*y_dim,1),99)
imagesc(wthresh(fmri_xcor,'h',activ_level));
colorbar

% Here, we try Fourier interpolation for better results
size_factor= x_dima/x_dim;
factiv= zeros(x_dima,y_dima);
factiv(1:x_dim, 1:y_dim)= fftshift(fft2(abs(fmri_xcor)));

factiv_interp= abs(ifft2(factiv)); 
activ_level= prctile (reshape(factiv_interp,x_dima*y_dima,1),99);
figure(8);
cmap2 = parula;
cmap1 = gray;
cmap = [cmap1;cmap2];
colormap(cmap)

imagesc(im_anat*100+wthresh(factiv_interp,'h',activ_level)*10);
colorbar
