<code - script_hemisphere_comparison.m>

%
%   Plot the distribution of beta values, in each hemisphere, for the main
%   effect covariate, as well as for the peak drift covariate, separating
%   between whole-fusiform mask and ROI only.
%   _______________________________________________
%   by Marcelo G Mattar (10/21/2013)
%   mattar@sas.upenn.edu

% Parameters
datasets = {'MDS_Multi'}; % Choose between 'MDS_RGI' and 'MDS_Multi'
masksize = 800;
% Data location
dataDir = '~/Data/FaceSpace/';

% Subjects to include
subjectList_RGI = {...
    % [list of subject ID strings redacted]
    };
subjectList_Multi = {...
    % [list of subject ID strings redacted]
    };
filename_endings = {'000' '005' '010' '015' '020' '025' '030' '035' '040' '045' '050' '055' '060' '065' '070' '075' '080' '085' '090' '095' '100'};
coefficient_values = 0:0.05:1;
numCoef = length(filename_endings);

for d = 1:length(datasets)
    % Determine which dataset we'll work on
    dataset = datasets{d};
    switch dataset
        case 'MDS_RGI'
            subjectList = subjectList_RGI;
            dataset_lowercase = '-mds_rgi';
        case 'MDS_Multi'
            subjectList = subjectList_Multi;
            dataset_lowercase = '-mds_multi';
    end
    numSubjects = length(subjectList);
    
    % Load masks
    mask_lh = load_mgh([dataDir 'ROIs/binaryMask_' dataset '_' num2str(masksize) '_lh.mgh']);
    mask_rh = load_mgh([dataDir 'ROIs/binaryMask_' dataset '_' num2str(masksize) '_rh.mgh']);
    fusiform_lh = load_mgh([dataDir 'ROIs/fusiform_lh.mgh']);
    fusiform_rh = load_mgh([dataDir 'ROIs/fusiform_rh.mgh']);
    mask_lh = logical(mask_lh);
    mask_rh = logical(mask_rh);
    fusiform_lh = logical(fusiform_lh);
    fusiform_rh = logical(fusiform_rh);
    
    % Initialize variables
    DriftBetas_lh = zeros(length(mask_lh),numSubjects,numCoef);
    DriftBetas_rh = zeros(length(mask_rh),numSubjects,numCoef);
    averageBetasROI_lh = zeros(numSubjects,length(filename_endings));
    averageBetasROI_rh = zeros(numSubjects,length(filename_endings));
    
    % Load data
    for k=1:numCoef
        DriftBetas_lh(:,:,k) = squeeze(load_mgh([dataDir 'group_glms/' dataset '/Drift_' filename_endings{k} '/lh.Drift_' filename_endings{k} '.fwhm10.mgh']));
        DriftBetas_rh(:,:,k) = squeeze(load_mgh([dataDir 'group_glms/' dataset '/Drift_' filename_endings{k} '/rh.Drift_' filename_endings{k} '.fwhm10.mgh']));
        averageBetasROI_lh(:,k) = mean(DriftBetas_lh(mask_lh,:,k))';
        averageBetasROI_rh(:,k) = mean(DriftBetas_rh(mask_rh,:,k))';
    end
    
    % Calculate summary distributions for beta_peak
    groupBetas_lh = squeeze(mean(DriftBetas_lh,2));
    groupBetas_rh = squeeze(mean(DriftBetas_rh,2));
    betaPeak_lh = max(groupBetas_lh,[],2);
    betaPeak_rh = max(groupBetas_rh,[],2);
    betaPeak_mask_lh = betaPeak_lh(mask_lh,:);
    betaPeak_mask_rh = betaPeak_rh(mask_rh,:);
    betaPeak_fusiform_lh = betaPeak_lh(fusiform_lh,:);
    betaPeak_fusiform_rh = betaPeak_rh(fusiform_rh,:);
    
    % Calculate summary distributions for main effect
    mainEffect_lh = squeeze(load_mgh([dataDir 'group_glms/' dataset '/DriftMainEffect/lh.DriftMainEffect.fwhm10.mgh']));
    mainEffect_rh = squeeze(load_mgh([dataDir 'group_glms/' dataset '/DriftMainEffect/rh.DriftMainEffect.fwhm10.mgh']));
    groupMainEffect_lh = mean(mainEffect_lh,2);
    groupMainEffect_rh = mean(mainEffect_rh,2);
    mainEffect_mask_lh = groupMainEffect_lh(mask_lh,:);
    mainEffect_mask_rh = groupMainEffect_rh(mask_rh,:);
    mainEffect_fusiform_lh = groupMainEffect_lh(fusiform_lh,:);
    mainEffect_fusiform_rh = groupMainEffect_rh(fusiform_rh,:);
    averageMainEffectROI_lh = mean(mainEffect_lh(mask_lh,:));
    averageMainEffectROI_rh = mean(mainEffect_rh(mask_rh,:));
    
    % Plot everything in a histogram
    h1 = figure(1);
    screen_size = get(0, 'ScreenSize');
    set(h1,'position',[0 0 screen_size(3)/2 screen_size(4)/2]);
    
    axis_limits = [-10/1000 15/1000 0 1000];
    subplot(221);
    [h,x] = hist(mainEffect_fusiform_lh,linspace(axis_limits(1),axis_limits(2),100));
    bar(x,h,'blue');
    hold on;
    [h,x] = hist(mainEffect_mask_lh,linspace(axis_limits(1),axis_limits(2),100));
    bar(x,h,'green');
    axis(axis_limits);
    ylabel('Main Effect','FontSize', 20);
    subplot(222);
    [h,x] = hist(mainEffect_fusiform_rh,linspace(axis_limits(1),axis_limits(2),100));
    bar(x,h,'blue');
    hold on;
    [h,x] = hist(mainEffect_mask_rh,linspace(axis_limits(1),axis_limits(2),100));
    bar(x,h,'green');
    axis(axis_limits);
    legend({'Fusiform gyrus','ROI'})
    
    axis_limits = [-1/10000 5/10000 0 600];
    subplot(223);
    [h,x] = hist(betaPeak_fusiform_lh,linspace(axis_limits(1),axis_limits(2),100));
    bar(x,h,'blue');
    hold on;
    [h,x] = hist(betaPeak_mask_lh,linspace(axis_limits(1),axis_limits(2),100));
    bar(x,h,'green');
    axis(axis_limits);
    ylabel('Beta peak','FontSize', 20);
    xlabel('LH','FontSize', 20);
    subplot(224);
    [h,x] = hist(betaPeak_fusiform_rh,linspace(axis_limits(1),axis_limits(2),100));
    bar(x,h,'blue');
    hold on;
    [h,x] = hist(betaPeak_mask_rh,linspace(axis_limits(1),axis_limits(2),100));
    bar(x,h,'green');
    axis(axis_limits);
    xlabel('RH','FontSize', 20);
    
end

</code>