

%%%%% 11/12/01

clear all;


P1 = spm_get(1,'*.img','Specify background image');

num = spm_input('Number of images to dispaly',1,'1 img|2 imgs|3 imgs',[1 2 3]);




spm_figure('GetWin', 'Graphics'); 

% Clear graphics window..
spm_clf

% Display background image..
h = spm_orthviews('Image', P1,[0.05 0.05 0.9 0.9]);
                        

% Getting the 'num' no.of ROIs to display over the background image..
P2 = spm_get(num,'*.img',{'Select image(s) to display'},pwd,0);%% FORMAT P = spm_get(n,Filter,Prompt,NewWDir,CmdLine);



for i=1:num,

    %c = spm_input('Colour','+1','m','Red blobs|Yellow blobs|Green blobs|Cyan blobs|Blue blobs|Magenta blobs',[1 2 3 4 5 6],1);
	c = spm_input(['Color for image:  ' P2{i}],'+1','m','Red blobs|Yellow blobs|Green blobs|Cyan blobs|Blue blobs|Magenta blobs',[1 2 3 4 5 6],1);
    colours = [1 0 0;1 1 0;0 1 0;0 1 1;0 0 1;1 0 1]; %colors in RGB format..
    %c = spm_input(['Color for image:  ' P{i}],'+1','m','Red blobs|black blobs|dark gray |light gray |white|Magenta blobs',[1 2 3 4 5 6],1);
    %colours = [1 0 0; 0 0 0 ;.2 .2 .2 ;.8 .8 .8;1 1 1 ;1 0 1];


    % [0.6 0 0.8] for purple etc..
    spm_orthviews('AddColouredImage', h, P2{i},colours(c,:));


    % Update the display..
    spm_orthviews('Redraw');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

input('choose a voxel/ROI')

cfg = spm_get(1,'SPMcfg.mat','Choose SPMcfg File');
load(cfg)

% Don't forget to run modifySPMcfg.m in case the image files are located in 
% a directory other than 

nsess = size(Sess,2);

posmm = spm_orthviews('Pos')
% vox_no = posmm  %need to convert to voxel
%vox_no = [15 29 20];


%%%%%%%%%%% Session information  
% Sess{3}
% 
%     BFstr: 'hrf (with time derivative)'
%     DSstr: 'Variable SOA '
%       rep: 0
%       row: [1x341 double]
%       col: [7 8 9 10 11 12 13 14]
%      name: {'ikc'  'ike'  'ix'  'ik'}
%       ind: {[1 2]  [3 4]  [5 6]  [7 8]}
%        bf: {[555x2 double]  [555x2 double]  [555x2 double]  [555x2 double]}
%        sf: {[8866x1 sparse]  [8866x1 sparse]  [8866x1 sparse]  [8866x1 sparse]}
%       pst: {[1x341 double]  [1x341 double]  [1x341 double]  [1x341 double]}
%       ons: {[6x1 double]  [15x1 double]  [10x1 double]  [10x1 double]}
%        Pv: {[]  []  []  []}
%     Pname: {''  ''  ''  ''}

TR = xX.RT;

dt = xX.dt;

regressors = xX.X;


images = {VY.fname};

%%%%%   Read in entire image volumes
%  FORMAT [Y,XYZ] = spm_read_vols(V,mask)
%  V    - vector of mapped image volumes to read in (from spm_vol)
%  mask - implicit zero mask?
%  Y    - 4D matrix of image data, fourth dimension indexes images
%  XYZ  - 3xn matrix of XYZ locations returned


% Get header information etc for images.

%get number of scans
nscan = size(images,2);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%arrange filenames
for i = 1:size(images,2)
   P(i,:) = images{i};
end

%Convert mm to voxel
header = spm_vol(P(1,:))
invT = inv(header.mat);
posmm = [posmm; 1];
vox_no = invT * posmm;

fprintf('\n Reading in Image Info \n')
V  = spm_vol(P);
 
%fprintf('\n Reading in All Images Specified in SPMcfg \n')
%[Y , XYZ]  = spm_read_vols(V); 
% Need to put a loop here for the different ROIs, numbering 1 through size(P2, 1).... [TODO]
interphold = -11; 
Vroi = spm_vol(P2);
% TODO: put a loop here, ie for each roi.....
% XYZroi corresponds to the mm positions of the entire image(s), with Yroi=1 where the image overlaps with the roi and 0 otherwise.
[Yroi , XYZroi]  = spm_read_vols(Vroi{1}); % Vroi{1} is a struct, since Vroi is a cell array[of structs]
XYZroi  = [XYZroi; ones(1, size(XYZroi, 2))];% appending XYZroi with a row of ones.
Voxels  = getfield(Vroi{1}, 'mat') \ XYZroi;

% lets now get those voxels alone which form part of the specified ROI:
xyz = [];
for z = 1 : max(Voxels(3, :))
	[x y] = find(Yroi(:,:,z));
	xyz   = [xyz; [x  y  z*ones(size(x, 1), 1)]];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nScan = size(xX.X,1);

vox_ts = [];

for i = 1:nScan
   
   vox_ts(i)= spm_sample_vol(V(i), vox_no(1), vox_no(2), vox_no(3),0);   
   ROIimage_ts(i)= mean(spm_sample_vol(V(i), xyz(:, 1), xyz(:, 2), xyz(:, 3), interphold));

   
end

% for now...to be modified asap...   [TODO]
vox_ts = ROIimage_ts; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% variable_names = {F_iX0.name};
% 
% if strcmp(variable_names{1}, 'effects of interest')
%     variable_names = variable_names(2:end);
% end

variable_names = [];
nv = 0;

% extract variable list and variable data

%%%% spm regressors in "regressors";

all_vars = {};

beg = 0;

for i = 1:nsess
    
    varnames = Sess{i}.name;
    
    nvars = size(varnames,2);
                  
    var_ons = Sess{i}.ons;
    
    inds = Sess{i}.ind;
    cols = Sess{i}.col;
    row = Sess{i}.row;
    
    sessions(i).beg = row(1) + beg;
    sessions(i).fin = beg + row(end);
    
    for j = 1:nvars
    
        
        variable_names{nv+1} = varnames{j};
        nv = nv+1;
        
        variables(i).vars(j).name = varnames{j};
        
        variables(i).vars(j).ons = var_ons{j};
        
        variables(i).vars(j).column = cols(inds{j}(1));
        
        variables(i).vars(j).row =  row;
        
        variables(i).vars(j).session = i;
        
        all_vars{end+1} = [i j];
        
        %%% add the rest : regressors ....
        
        
    end
        
end





var_no = menu('Which Variable ?', variable_names);

var_inds = all_vars{var_no};

var_name = variable_names{var_no};
%var_name = var_name(12:end);



var_sess = variables(var_inds(1)).vars(var_inds(2)).session;
%var_sess = str2num(var_sess);


    
var_ons = variables(var_inds(1)).vars(var_inds(2)).ons;

var_col = variables(var_inds(1)).vars(var_inds(2)).column;

var_row = variables(var_inds(1)).vars(var_inds(2)).row;  %- sessions(var_sess).beg + 1;

var_reg = regressors(:,var_col);

var_reg = var_reg(var_row,:);


var_beg = sessions(var_sess).beg;
var_fin = sessions(var_sess).fin;



ymin_reg = min(var_reg);
ymax_reg = max(var_reg);

reg_time_series = var_reg;

% raw_time_series = Y(vox_no(1),vox_no(2),vox_no(3),var_beg:var_fin);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%EXTRACT RAW TIME SERIES FROM IMAGES%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
raw_time_series = vox_ts(var_beg:var_fin);

raw_time_series = squeeze(raw_time_series);  %FINALLY get 1xn instead of 1x1x1xn

reg_time_series = reg_time_series';
raw_time_series = raw_time_series';


ymin_reg = min(reg_time_series);
ymax_reg = max(reg_time_series);
ymin_raw = min(raw_time_series);
ymax_raw = max(raw_time_series);

time_axis_raw = ( var_row - sessions(var_sess).beg + 1)* TR;
time_axis_reg = ( var_row - sessions(var_sess).beg + 1)* TR;



%%%% raw data time series

%%% plot whole raw data and regressor time-series for a given voxel and variable
wts = figure;

warning off

vox_no = round(vox_no);

%%% raw data
%sp1 = subplot(2,1,1)
sp1 = subplot(3,1,1)
plot(time_axis_raw,raw_time_series,'r') %%%% 
title(['Raw Data Time Series for variable ' var_name ' from voxel ' num2str(vox_no(1)) ' ' num2str(vox_no(2)) ' ' num2str(vox_no(3)) ])
xlabel('TIME (sec)')
ylabel('BOLD')


% for i = 1:size(var_ons,1)
%     ons = var_ons(i);
%     hold on    
%     text(ons,ymax_raw,var_name)
%     x = [ons ons]; y = [ymin_reg-1 ymax_reg-1];
%     x = [ons ons]; y = [ymin_raw-1 ymax_raw+1]; line(x,y)
% 
% end
hold off

axis tight
axis([ -inf inf ymin_raw-1 ymax_raw+1])


%%% spm regressor
%sp2 = subplot(2,1,2)
sp2 = subplot(3,1,2)
plot(time_axis_reg,reg_time_series,'k') %%%% 
title(['variable ' var_name '  for voxel ' ...
      num2str(vox_no(1))  ' ' num2str(vox_no(2))  ' ' num2str(vox_no(3)) ])
xlabel('TIME (sec)')
ylabel('regressors')



sX   = xX.X(Sess{1}.row,Sess{1}.col);
rX    = sX(:,Sess{1}.ind{1});
q=nscan;
%spectral analyses
sp3 = subplot(3,1,3)

%axes('Position',[0.580,0.405,0.325,0.225])
gX    = abs(fft(rX)).^2;
gX    = gX*diag(1./sum(gX));
q     = size(gX,1);
Hz    = [0:(q - 1)]/(q*xX.RT);
q     = 2:fix(q/2);
plot(Hz(q),gX(q,:))
xlabel('Frequency (Hz)')
ylabel('spectral density')
%title('Frequency domain')
%grid on
axis tight





for i = 1:size(var_ons,1)
    ons = var_ons(i);
%     hold on    
    %text(ons,ymax_reg,var_name)
    
    x_ind = round(ons/TR);
%     x= ons;
%     y = raw_time_series(x_ind);
%     y = y+ymax_raw/20;
%     txt = ['x'];
%     subplot(sp1)
%     t1=text(x,y,txt);
%     set(t1,'FontWeight','bold')
%     set(t1,'FontSize',12)
%     set(t1,'Color', [0 0 1])
    
    subplot(sp1)
    hold on
    rect = get(sp1,'position');
    %tk = rect(4)/10;
    x = [ons ons];
    yc = raw_time_series(x_ind);
    tk = (ymax_raw-ymin_raw)/10;
    y = [yc-tk yc+tk];
    lin1 = line(x,y);
    set(lin1,'Color',[0 0 1],'LineWidth',2);

    subplot(sp2)
    hold on
    rect = get(sp2,'position');
    %tk = rect(4)/10;
    x = [ons ons];
    yc = reg_time_series(x_ind);
    tk = (ymax_reg-ymin_reg)/5;
    y = [yc-tk yc+tk];
    lin2 = line(x,y);
    set(lin2,'Color',[1 0 0],'LineWidth',2);
    
    
    %x = ons;
%     y = reg_time_series(x_ind);
%     y = y+ymax_reg/20;
%     txt = ['x'];
%     subplot(sp2)
%     t2=text(x,y,txt);
%     set(t2,'FontWeight','bold')
%     set(t2,'FontSize',12)
%     set(t2,'Color', [1 0 0])
    
    %x = [ons ons]; 
    %y = [ymin_reg-1 ymax_reg+1]; line(x,y)
end
hold off

axis tight
axis([ -inf inf ymin_reg-1 ymax_reg+1])



warning on



%[xX,Sess] = spm_fmri_design_show;
s=1; 
i=1;
rX    = sX(:,Sess{s}.ind{i});


%%%%%%%%%%%%% Design Specs          
% xsDes = 
% 
%                     Design: 'Variable SOA '
%            Basis_functions: 'hrf (with time derivative)'
%         Number_of_sessions: '3'
%     Conditions_per_session: '1  2  4  '
%         Interscan_interval: '1.50'
%           High_pass_Filter: '[min] Cutoff period 49 seconds'
%            Low_pass_Filter: 'none'
%     Intrinsic_correlations: 'none'
%         Global_calculation: 'mean voxel value'
%         Grand_mean_scaling: 'session specific'
%       Global_normalisation: 'None'
%             
% xX = 
% 
%          X: [908x17 double]
%         dt: 0.0577
%         RT: 1.5000
%         iH: []
%         iC: [1 2 3 4 5 6 7 8 9 10 11 12 13 14]
%         iB: [15 16 17]
%         iG: []
%     Xnames: {1x17 cell}
%          K: {[1x1 struct]  [1x1 struct]  [1x1 struct]}
%        xVi: [1x1 struct]                       



% s=uicontrol('Style','slider','Parent',figwin,...
	%	    'Position',[200 150 200 30],...
	%	    'Min',min_slide,'Max',max_slide,...
	%	    'SliderStep',[min_step max_step],...
    %                'Callback','first_movie(''Scroll'');',...
     %               'Userdata',m_struct,'String','Frame number');


