
% read and return data from any images from an identified region of interest
% the region of interest is identified within SPM surfing the spm table
% 
% the code comes from spm_regions - one needs to estimate a contrast and get 
% the spm table (volume button) - from there the user moves on a given cluster
% then calls Easy_ROI in the command window. One gives a name, select box, spere or cluster.
% At variance with the VOI button, the user next select some images
% and the function returns the data, the mean, the sd and the 1st eigenvector
% - across voxels for each image (mean, sd)
% - across voxels and images (mean, sd, 1st eigenvariate)
% 
% Everything is saved in .mat with the name of the region
%
% xY.xyz           coordinate of the maximum
%   .name          name of the ROI (user choice)
%   .def           sphere, box, cluster (user choice)
%   .XYZmm         coordinate of voxels included in the analysis
%   .raw           extracted values
%   .eigen         1st eigen value across voxels for each image/subject (1 vector)
%   .mean          average across voxels for each image/subject (1 vector)
%   .sd            standard deviation
%   .overallmean   overall mean - average across voxels and images/subject (1 value)
%   .overallsd     overall standard deviation
% 
% cyril pernet june 2006


%-----------------------------------------------------------------------
% The beginning is from spm_regions.m file (Karl Friston v192 2005-06-23 13:55:27Z) 
%-----------------------------------------------------------------------
 
% get figure handles
%-----------------------------------------------------------------------
Finter = spm_figure('GetWin','Interactive');
Fgraph = spm_figure('GetWin','Graphics');
header = get(Finter,'Name');
set(Finter,'Name','ROI data extraction') % just modified the text
try
	xY;
catch
	xY = {};
end


%-Find nearest voxel [Euclidean distance] in point list in Y.mad
%-----------------------------------------------------------------------
if ~length(xSPM.XYZmm)
	spm('alert!','No suprathreshold voxels!',mfilename,0);
	Y = []; xY = [];
	return
end
try
	xyz     = xY.xyz;
catch
	[xyz,i] = spm_XYZreg('NearestXYZ',...
		      spm_XYZreg('GetCoords',hReg),xSPM.XYZmm);
	xY.xyz  = xyz;
end


% and update GUI location
%-----------------------------------------------------------------------
spm_XYZreg('SetCoords',xyz,hReg);


%-----------------------------------------------------------------------
% here add the spm_input line
%-----------------------------------------------------------------------
xY.name = spm_input('name of region','!+1','s','VOI');



%-Specify VOI
%-----------------------------------------------------------------------

if ~isfield(xY,'def')
	xY.def    = spm_input('VOI definition...','!+1','b',...
			{'sphere','box','cluster'});
end
Q       = ones(1,size(xSPM.XYZmm,2));


switch xY.def

	case 'sphere'
	%---------------------------------------------------------------
	if ~isfield(xY,'spec')
		xY.spec = spm_input('VOI radius (mm)','!+0','r',0,1,[0,Inf]);
	end
	d     = [xSPM.XYZmm(1,:) - xyz(1);
		 xSPM.XYZmm(2,:) - xyz(2);
		 xSPM.XYZmm(3,:) - xyz(3)];
	Q     = find(sum(d.^2) <= xY.spec^2);

	case 'box'
	%---------------------------------------------------------------
	if ~isfield(xY,'spec')
		xY.spec = spm_input('box dimensions [x y z] {mm}',...
			'!+0','r','0 0 0',3);
	end
	Q     = find(all(abs(xSPM.XYZmm - xyz*Q) <= xY.spec(:)*Q/2));

	case 'cluster'
	%---------------------------------------------------------------
	[x i] = spm_XYZreg('NearestXYZ',xyz,xSPM.XYZmm);
	A     = spm_clusters(xSPM.XYZ);
	Q     = find(A == A(i));
end

% voxels defined
%-----------------------------------------------------------------------
spm('Pointer','Watch')




%-----------------------------------------------------------------------
% now instead of looking at the current data and adjust, we load new images
% and compute on the selected voxels on the selected volumes
%-----------------------------------------------------------------------

% select the images

try
   P = spm_select(Inf, '.*\.img$','Select images you want to read');
catch
   P = spm_get(Inf,'.img','Select images you want to read');
end

V = spm_vol(P);


%-Get raw data 
%-----------------------------------------------------------------------
y = spm_get_data(V,xSPM.XYZ(:,Q));
xY.XYZmm = xSPM.XYZmm(:,Q);

% in case we got NaNs
N = isnan(y);
K = (sum(find(N)) > 0);

if K ==1
    col = size(y);
    col = col(2);
    for i = 1:col
        N = isnan(y(:,i));
        K = (sum(find(N)) > 0);
        if K ==1
            Badcol(i)= i;
        end
    end
    warndlg('Some voxels were outside the brain, less voxels will be used','NaNs');

    I = find(Badcol);
    y(:,I)=[];
    
    nb = size(xY.XYZmm);
    nb = nb(2);
    newnb = size(y);
    newnb = newnb(2);
    finalnb = nb - newnb;
    fprintf('%g voxels were removed from the ROI',finalnb);
    disp(' ');
    clear col Badcol I nb newnb finalnb
end

clear N K

% Compute what we need
%-----------------------------------------------------------------------


% the result for each images is given by
mean_each = mean(y,2);
sd_each = std (y,1,2); % note flag = 1, i.e. matlab divides by n

% the result across voxels and images is given by
mean_all = mean(mean_each);
sd_all = std (mean_each,1);



% compute regional response in terms of first eigenvariate
%-----------------------------------------------------------------------
% note here that we do not take into account the whitening part of
% spm_regions and there is no adjustement - the idea is that you may want to
% know the fit in many images from two conditions (con images condition1 
% then redo with con images condition2) and across subjects but the area 
% selected may not be homogenous and the 1st eigenvariate is probably better 
% than the mean - if it is homogenous then results must converge
%-----------------------------------------------------------------------


[m n]   = size(y);
if m > n
	[v s v] = svd(spm_atranspa(y));
	s       = diag(s);
	v       = v(:,1);
	u       = y*v/sqrt(s(1));
else
	[u s u] = svd(spm_atranspa(y'));
	s       = diag(s);
	u       = u(:,1);
	v       = y'*u/sqrt(s(1));
end
d       = sign(sum(v));
u       = u*d;
v       = v*d;
Y       = u*sqrt(s(1)/n);


% % set in structure
% %-----------------------------------------------------------------------
xY.raw          = y;         % extracted values
xY.eigen        = Y;         % eigenvariate - eigen values across voxels for each image
xY.mean         = mean_each; % average across voxels for each image
xY.sd           = sd_each;   % standard deviation
xY.overallmean  = mean_all;  % overall mean - average across voxels and images
xY.overallsd    = sd_all;    % overall standard deviation


% % show results along images
% %------------------------------------------------------------------------
spm_results_ui('Clear',Fgraph);
figure(Fgraph);

subplot(2,2,3)
plot(xY.eigen,'b');
hold;
plot(xY.mean,'r');
grid;
title(['1st eigenvariate (blue) and Mean (red)'],'FontSize',10);
hold off;

%subplot(2,2,1)
%plot(xY.eigen,xY.mean);
%title(['1st eigenvariate vs Mean: ' xY.name],'FontSize',10);
%axis tight square

% % show position
% %------------------------------------------------------------------------
subplot(2,2,4)
spm_dcm_display(xY,[],[],[[1 0 0];[0 1 0]]',64)


% % save
% %-----------------------------------------------------------------------

save ([xY.name], 'xY');

% clear variables from space
clear x xyz P V y Y s v u mean_each sd_each mean_all sd_all
clear xY



