% ROI analysis
% 
% Xu Cui, 28 February 2005
%
% ROI stands for region of interest. The region is predefined usually from
% glm contrast. For example, you find visual cortex is activated when a
% flash is presented through glm analysis, you may then do more analysis on
% this particular region.
% 
% ROI analysis usually means a plot of BOLD signal in this region against
% time. In some cases, you may also want to do correlation analysis. I will
% explain how to do these analysis below. I will also give an example.
%
% Before you do any ROI analysis, you should have the following data in hand:
% 1. Preprocessed image files which contains the BOLD signals
%   The file names are usually of the format swaf*.img and swaf*.hdr.
% 2. The onset time of each event of interest.
% 3. The coordinates of the voxels in the region of interest, or a mask
% image.
% 4. If you want to find the correlation between two regions, you need
% another set of 3.
% 5. If you want to find the correlation between BOLD signal and the
% stimulus property (such as the intensity of flash), you need data of this
% property.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The procedure:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% I. BOLD signal against time
%
% 1. Load the image files of a single subject into a 4-D matrix. The first
% three dimensions are spatial dimensions and the last one, temporal.
P = spm_get('Files', pwd, 'swaf*.img'); % read the file names. You should make sure the order is correct.
V = spm_vol(P);
signal = spm_read_vols(V); % signal is a 4-D matrix

% 2. Average the signals across voxels in the region of interest
% If you don't have a list of ROI coordinates but have the mask image file of that ROI, convert it to matrix coordinate:
coord = mask2coord('mymask.img'); % mymask.img contains my region of interest. coord is a Nx3 matrix.
meansignalROI = zeros(size(signal,4), 1);
tmp0 = zeros(size(coord,1),1);
for i=1:size(signal,4)
		tmp = signal(:,:,:,i);
        for j=1:size(coord,1)
	        tmp0(j) = tmp(coord(j,1),coord(j,2),coord(j,3));
	    end
		meansignalROI(i) = mean(tmp0);
        tmp0 = 0*tmp0;
end
% meansignalROI is now a vector, with 1 temporal dimension

% 3. Detrend the signal, also calculate the baseline.
detrendsignal = detrend(meansignalROI);
baseline = mean(meansignalROI);
% I used linear detrend above. If you find the trend of the signal is not
% linear (by plotting the raw data), you should use nonlinear detrending.
% You may use moving average to get the nonlinear trend and subtract it
% from the raw signal.

% 4. Cut the signal at the onset times of the event of interest, and align
% the resulted pieces. 
% I assume you already have the event onset time (in seconds) called 
% eventOnset, which is a vector. Also I assume you are only interested in 
% BOLD change within a certain time window after onset, called duration.
% (Typically, duration = 20 seconds is enough). I also assume the temporal
% resolution of BOLD signal is 2s (which is true in HNL).
for j=1:length(eventOnset)
    for i=1:duration/2
        s0 = eventOnset(j)/2 + 1 + i - 1;
        eventResponse(j,i) = detrendsignal(floor(s0)) + (detrendsignal(ceil(s0))-detrendsignal(floor(s0))).*(s0 - floor(s0));
    end
end
eventResponse = eventResponse/baseline;
% now eventResponse contains the relative BOLD change for every event. Each
% row is BOLD signal over time for a certain event. Note that the code has
% taken care of the events whose onset is between two adjacent sampling by
% interpolation.

% 5. Average the signal across events.
meanEventResponse = mean(eventResponse); 
% now meanEventResponse is a vector.

% 6. Repeat step 1 to 5 for other types of events and other subjects.
% now you get meanEventResponse for every subject. Assume you save it in a matrix mER, where the ith row is the meanEventResponse for the ith subject.

% 7. Average the signals across subjects; also calculate the standard error.
meanAll = mean(mER);
seAll = std(mER)/sqrt(size(mER,1));

% 8. Plot the result.
errorbar([0:2:duration-1],meanAll, seAll); % with x-axis time of unit second.
% you may also want to plot the BOLD for other types of events for
% comparison.

% 9. If you want to compare BOLD of two types of events (e.g. visual vs
% auditory), you do a two-sample t test on every time point. For comparison
% across more than two types of event, do repeated ANOVA.

% II. correlation between two regions
%
% In step 5 in the above procedure, you get the BOLD change for a
% particular region (denoted by R1). You repeat the procedure for another
% region R2. Now you have two vectors, s1 and s2.
corr = corrcoef(s1,s2); % corr is the correlation between this two regions
% You may also want to know if there are any delay of the activation of one
% region compared to the other region. You do a cross-correlogram, which is
% a plot of correlation against time shift. Below is an example to
% calculate the correlation when the shift is 5 scan (10 seconds).
shift = 5;
corr = corrcoef(s1(1:end-shift),s2(shift+1:end));
% you then repeat for other shift size and plot correlation against shift.
% If the peak is at shift=0, then there is no delay between two regions;
% otherwise there is.

% In some case you may not want to calculate the correlation between two
% averaged signals. You may want to calculate the correlation between raw
% signals, and average the correlations. I myself prefer the latter.

% III. correlation between BOLD and stimulus property.
% In step 4 of the above procedure, you already have the BOLD signal for
% every event. You also have the property of the stimulus (e.g. flash
% intensity). You want to calculate the correlation between them. 
%
% You want to summarize the activation for a single event into one
% number. You could use the peak activation. I recommend to use the mean
% BOLD around the peak time (which are usually at time 4, 6, and 8
% seconds.)
meanPeak = mean(eventResponse(:,[3 4 5]),2);
plot(stimulusProperty, meanPeak);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Note: the data are made up.
%
% We have preprocessed data at /mnt/nfs/shared/unix_gateway/roi/data
% We have two regions of interest, they are at 
% /mnt/nfs/shared/unix_gateway/roi/region1.img
% /mnt/nfs/shared/unix_gateway/roi/region2.img
% We have one type of event - flash, the onset of which is eventOnset, the
% intensity is intensity.
% 

os = system_dependent('getos'); % which operating system are you using?

% onset time of event flash
eventOnset =[0          60         120         180         240         300         360         420         480         540 ...
         600         660         720         780         840         900         960        1020        1080        1140]; % in seconds
duration = 20; % seconds
% use xjView to display where are the regions
xjView region1.img
xjView region2.img
% read mask file and get the roi coordinates
if ~isempty(findstr(os, 'Linux')) % if you are in linux
    coords1 = mask2coord('/mnt/nfs/shared/unix_gateway/roi/region1.img');   % coordinates of region1
    coords2 = mask2coord('/mnt/nfs/shared/unix_gateway/roi/region2.img');   % coordinates of region2   
elseif ~isempty(findstr(os, 'Windows')) % if you are in windows
    coords1 = mask2coord('W:\roi\region1.img');   % coordinates of region1
    coords2 = mask2coord('W:\roi\region2.img');   % coordinates of region2   
end
intensity = [    1.7643    4.0658    0.0493    0.6945    1.0138    0.9936    3.0190    1.3609    0.9941    0.0764    3.7339    2.2255 ...
    4.6591    2.3300    2.0932    4.2311    2.6258    1.0132    3.3607    4.1906]; % intensity of flash
% load preprocessed imaging data
if ~isempty(findstr(os, 'Linux')) % if you are in linux
    P = spm_get('Files', '/mnt/nfs/shared/unix_gateway/roi/data', 'swaf*.img'); % read the file names. You should make sure the order is correct.
elseif ~isempty(findstr(os, 'Windows')) % if you are in windows
    P = spm_get('Files', 'W:\roi\data', 'swaf*.img'); % read the file names. You should make sure the order is correct.
end
V = spm_vol(P);
signal = spm_read_vols(V); % signal is a 4-D matrix

% now average across voxels in region 2
coord = coords2;
meansignalROI = zeros(size(signal,4), 1);
tmp0 = zeros(size(coord,1),1);
for i=1:size(signal,4)
	tmp = signal(:,:,:,i);
    for j=1:size(coord,1)
        tmp0(j) = tmp(coord(j,1),coord(j,2),coord(j,3));
    end
	meansignalROI(i) = mean(tmp0);
    tmp0 = 0*tmp0;
end

% Detrend the signal, also calculate the baseline.
detrendsignal = detrend(meansignalROI);
baseline = mean(meansignalROI);

% cut and past/align
for j=1:length(eventOnset)
    for i=1:duration/2
        s0 = eventOnset(j)/2 + 1 + i - 1;
        eventResponse(j,i) = detrendsignal(floor(s0)) + (detrendsignal(ceil(s0))-detrendsignal(floor(s0))).*(s0 - floor(s0));
    end
end
eventResponse = eventResponse/baseline;

% Average the signal across events.
meanEventResponse = mean(eventResponse); 
% calculate the standard error
stdEventResponse = std(eventResponse)/sqrt(size(eventResponse,1)); 

% plot BOLD vs time for region 2
figure;errorbar([0:2:duration-1],meanEventResponse,stdEventResponse)
xlabel('time in second')
ylabel('BOLD')

% Assume BOLD signal in regions 2 also respond a second type of event -
% beep. And assume the eventResponse for beep is:
eventResponseBeep =[
   -0.0001    0.0003    0.0005   -0.0004    0.0028    0.0012    0.0011    0.0002   -0.0003   -0.0009
   -0.0001    0.0007   -0.0005    0.0001    0.0021    0.0022    0.0012    0.0002   -0.0003   -0.0003
    0.0002    0.0003    0.0012   -0.0002    0.0017    0.0026    0.0015    0.0008    0.0001   -0.0005
   -0.0006   -0.0003    0.0008   -0.0000    0.0015    0.0014   -0.0002    0.0003   -0.0010    0.0000
    0.0002   -0.0007   -0.0005   -0.0005    0.0004    0.0004   -0.0002    0.0000   -0.0007   -0.0002
    0.0003    0.0006    0.0001    0.0005    0.0021    0.0012    0.0006   -0.0003    0.0004    0.0001
   -0.0004   -0.0001    0.0002    0.0002    0.0014    0.0025    0.0010    0.0008   -0.0005    0.0003
    0.0000    0.0014    0.0000    0.0007    0.0014    0.0014    0.0003    0.0006   -0.0001   -0.0001
   -0.0002   -0.0000    0.0002   -0.0005    0.0016    0.0016    0.0006   -0.0010   -0.0001   -0.0005
    0.0004   -0.0003   -0.0003   -0.0001    0.0015    0.0014    0.0013    0.0003    0.0003    0.0005
    0.0001   -0.0003    0.0002    0.0007    0.0018    0.0017    0.0010    0.0001   -0.0005   -0.0010
    0.0004    0.0007   -0.0000    0.0006    0.0001    0.0010    0.0005    0.0004    0.0004    0.0001
   -0.0007   -0.0003    0.0004   -0.0003    0.0010    0.0018    0.0006    0.0006   -0.0006    0.0001
    0.0001    0.0002    0.0005    0.0007    0.0013    0.0019    0.0010   -0.0001   -0.0007    0.0001
   -0.0006    0.0012    0.0002    0.0017    0.0021    0.0010    0.0021    0.0001   -0.0000   -0.0001
    0.0001    0.0003    0.0001    0.0008    0.0017    0.0015    0.0014   -0.0006   -0.0003   -0.0003
   -0.0003    0.0001    0.0008    0.0008    0.0025    0.0014    0.0013    0.0004   -0.0005   -0.0001
    0.0003    0.0006    0.0001    0.0002    0.0007    0.0012    0.0010   -0.0000    0.0001   -0.0005
   -0.0004    0.0005    0.0001    0.0013    0.0017    0.0027    0.0004    0.0009   -0.0003   -0.0001
    0.0000   -0.0001   -0.0000   -0.0001    0.0013    0.0014    0.0003    0.0006    0.0005    0.0003
    0.0001    0.0002    0.0007    0.0012    0.0004    0.0012    0.0001   -0.0003   -0.0001   -0.0005
   -0.0001   -0.0001   -0.0003    0.0005    0.0017    0.0019    0.0002   -0.0003    0.0003    0.0005
    0.0008   -0.0000    0.0002    0.0009    0.0023    0.0012    0.0013    0.0006    0.0009   -0.0003
    0.0006   -0.0004    0.0006   -0.0003    0.0012    0.0010    0.0008   -0.0006   -0.0010   -0.0010
   -0.0003   -0.0001   -0.0007    0.0004    0.0016    0.0020    0.0010    0.0003   -0.0006   -0.0002
   -0.0010   -0.0001   -0.0005   -0.0000    0.0013    0.0022    0.0006   -0.0001    0.0007    0.0006
   -0.0011    0.0004    0.0003    0.0004    0.0019    0.0017    0.0004    0.0014   -0.0006    0.0010
   -0.0004   -0.0001    0.0003    0.0002    0.0025    0.0014    0.0008    0.0003   -0.0000   -0.0003
   -0.0011   -0.0000    0.0001    0.0011    0.0007    0.0014    0.0007    0.0009    0.0001    0.0005
   -0.0001    0.0007    0.0003    0.0007    0.0018    0.0012    0.0009   -0.0000    0.0004   -0.0006];
meanEventResponseBeep = mean(eventResponseBeep);
stdEventResponseBeep = std(eventResponseBeep)/sqrt(size(eventResponseBeep,1)); 
% plot flash and beep together
figure;errorbar([0:2:duration-1],meanEventResponse,stdEventResponse)
hold on;errorbar([0:2:duration-1],meanEventResponseBeep,stdEventResponseBeep,'r')
xlabel('time in second')
ylabel('BOLD')
legend('','flash','','beep')

% Do two sample t test on each time point
for ii=1:10
    [h,p(ii)]=ttest2(eventResponseBeep(:,ii),eventResponse(:,ii));
end
% we find p<0.05 for the 5th to 7th time point.

% We want to know if the activity of region 1 is correlated with the
% intensity of  flash
coord = coords1;
meansignalROI = zeros(size(signal,4), 1);
tmp0 = zeros(size(coord,1),1);
for i=1:size(signal,4)
		tmp = signal(:,:,:,i);
        for j=1:size(coord,1)
	        tmp0(j) = tmp(coord(j,1),coord(j,2),coord(j,3));
	    end
		meansignalROI(i) = mean(tmp0);
        tmp0 = 0*tmp0;
end
detrendsignal = detrend(meansignalROI);
baseline = mean(meansignalROI);
for j=1:length(eventOnset)
    for i=1:duration/2
        s0 = eventOnset(j)/2 + 1 + i - 1;
        eventResponse1(j,i) = detrendsignal(floor(s0)) + (detrendsignal(ceil(s0))-detrendsignal(floor(s0))).*(s0 - floor(s0));
    end
end
eventResponse1 = eventResponse1/baseline;
meanEventResponse1 = mean(eventResponse1); 
stdEventResponse1 = std(eventResponse1)/sqrt(size(eventResponse1,1)); 
figure;errorbar([0:2:duration-1],meanEventResponse1,stdEventResponse1)
xlabel('time in second')
ylabel('BOLD')
% now find the correlation
meanPeak = mean(eventResponse1(:,[3 4 5]),2);
plot(intensity, meanPeak,'.');
corr = corrcoef(intensity, meanPeak);
disp(corr(1,2))
% we find the correlation is 0.98, which means the BOLD signal in region 1
% and flash intensity are highly correlated. But you will not get such nice
% correlation in reality -- this data is made up.

% We want to know the correlation between region 1 and region 2.
shift = [-5:5]; % in scan number
c = corrlag(meanEventResponse, meanEventResponse1,shift);
figure;plot(shift, c,'o-');
xlabel('shift in scan number')
% we see the correlation peaks at shift = 2, which means region 1 leads
% reagion 2 by 2 scans (4 seconds).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% I said I made the data up. How? See below:
eventOnset1 = [1:30:603];
eventOnset2 = [1:30:603]+2;
intensity = [    1.7643    4.0658    0.0493    0.6945    1.0138    0.9936    3.0190    1.3609    0.9941    0.0764    3.7339    2.2255 ...
    4.6591    2.3300    2.0932    4.2311    2.6258    1.0132    3.3607    4.1906    0.0982];
for ii=1:size(P,1)
    hrf = spm_hrf(2);
    hrf = [hrf; zeros(60,1)];
    jj = ii-eventOnset1;
    pos = find(jj>=0);
    jj = jj(pos(end))+1;
    coords1 = mask2coord('/mnt/nfs/shared/unix_gateway/roi/region1.img');
    intensity1 = 500+randn(size(coords1,1),1)+hrf(jj)*intensity(floor((ii-1)/30)+1);
    coords2 = mask2coord('/mnt/nfs/shared/unix_gateway/roi/region2.img');
    jj = ii-eventOnset2;
    pos = find(jj>=0);
    if isempty(pos)
        jj = 1;
    else
        jj = jj(pos(end))+1;
    end
    
    intensity2 = 500+randn(size(coords2,1),1)+hrf(jj);    
    mni2mask(cor2mni([coords1;coords2]), P(ii,:), [intensity1;intensity2], P(ii,:));
end