MATLAB script for the first approach

From Psych 221 Image Systems Engineering
Revision as of 08:53, 9 December 2009 by imported>Psych 204
Jump to navigation Jump to search
   clear all
   format long g
   %set directory for the two subjects to be compared
   if ispc
   subjADir='Y:\New_Localizers\dr021108-1p5mm-3mm';
   subjBDir='Y:\New_Localizers\kw011708-1p5mm-3mm';
   else
   %     subjADir='/biac2/kgs/projects/New_Localizers/dr021108-1p5mm-3mm';
   % %     subjBDir='/biac2/kgs/projects/New_Localizers/kw011708-1p5mm-3mm';
   %     subjBDir='/biac2/kgs/projects/New_Localizers/kgs020408-1p5mm-3mm';
   %     
  home = '/biac2/kgs/projects/SameDifferent/'
   %    subjADir='/biac2/kgs/projects/Prosopagnosia/fMRIData/112909whFMRI';    
   %    subjBDir='/biac2/kgs/projects/Prosopagnosia/fMRIData/52309AF';    
   % subjADir = '/biac2/kgs/projects/Kids/fmri/localizer/adult_al_22yo_051108';
   % subjBDir ='/biac2/kgs/projects/Kids/fmri/localizer/adult_jc_27yo_052408';
   subjADir ='/biac2/kgs/projects/Kids/fmri/localizer/adult_kw_25yo_090308';

subjBDir = '/biac2/kgs/projects/Kids/fmri/localizer/adult_dy_25yo_041908';

end

% %choose an ROI in subject A % roiA='lh_LOfaces_event_loc'; % roiA='lh_PPA_event';

% roiA='lPPA'; % roiA='rSTSfaces';

% get all the ROIs from subject A % cd(subjADir); % cd Inplane/ROIs

%finds all the file names and chops off the .mat

% rois=dir('*.mat'); % for i=1:size(rois,1) % roisA(i) = cellstr(strrep(rois(i).name,'.mat',)); % end

% golijeh kid data rois roisA={'lFFA_MBvAC_p3','rFFA_MBvAC_p3','lLO_ACvT_p3','rLO_ACvT_p3','lMT_p4_al', 'rMT_p4_al',...

   'lPPA_IOvAC_p3','rPPA_IOvAC_p3'};

% ,'lSTS_MBvAC_p3', 'rSTS_MBvAC_p3' %load ROIs from subject B % roiB={'lFFA','lPPA','rPPA','pSTSfaces'};

% get all the ROIs from subject A % cd(subjBDir); % cd Inplane/ROIs

%finds all the file names and chops off the .mat

% rois=dir('*.mat'); % for i=1:size(rois,1) % roisB(i) = cellstr(strrep(rois(i).name,'.mat',)); % end

% golijeh kid data rois roisB={'lFFA_MBvAC_p3','rFFA_MBvAC_p3','lLO_ACvT_p3','rLO_ACvT_p3','lMT_p4_al','rMT_p4',...

   'lPPA_IOvAC_p3','rPPA_IOvAC_p3'};

% ,'lpSTS_MBvACIO_p3', 'rSTS_MBvAC_p3' %want to get the motion corrected data type which is usually 3 dt = 3; %'MotionComp_RefScan5'; %then pick the scans from the motion corrected data that were (in this %case) taken from the event related adaptation experiment for %dr021108-1p5mm-3mm this is scans 1:8 % for waldemar scan went % mtloc,eccbias1,eccbias2,4xretinotopy,objloc1,objloc2 % so to get our scans in order we have % scan = [8 9 1 2 3];

% golijeh kid scans scan = [1 2 3 4]; %find the betas over time for all the ROIs loaded. each voxel represented %as a vector with the B for each category/stimulus

%go to path


% .. for subject A brain' %go to subjects data directory cd(subjADir); %initialize our hidden inplane set to the right data type,scan, and roi hiA = initHiddenInplane(dt, scan, roisA); %then we grab the multivoxel data from the ROI %this should give us mvA which is a struct with lots of useful stuff %in particular %mvA.trials has information about conditions and order of stimuli etc %mvA.coord has the coordinates for each voxel in the roi should be 3 x %numVoxels %mvA.tSeries has the time series for each voxel in % signal change %mvA.roi had roi info

for i = 1:size(roisA,2)

   mvA{i} = mv_init(hiA, roisA{i}, scan, dt);

% mvA{i} = mv_applyGlm(mvA{i}); end %then we apply our glm to each of the voxels in the roi %we will probably want to set these parameters ourselves so that we know %what happened %this gives us %mvA.glm has the betas, the design matrix and so on %for a block design there will be as many betas as there are conditions %for an event related design there will be timepoints x conditions x voxels %betas


% scan=[1 2 3 4 5]; % scan=[1 2];

% kevin scans scan=[1 2 3 4];

% .. for all ROIs in subj B %instead invoke hidden inplane with multiple rois. then loop through the %rois when comparing to subjectA %go to subject b directory cd(subjBDir); dt=3; %load data with all rois hiB = initHiddenInplane(dt,scan,roisB);

for i = 1:size(roisB,2);

   mvB{i} = mv_init(hiB, roisB{i}, scan, dt);

% mvB{i} = mv_applyGlm(mvB{i}); end


%find the correlation coefficient between a voxel in subject A and each %voxel in each ROI in subj B

%get the number of voxels by getting the size and then the number of %columns

%trim the size of the tSeries down (if necessary) % for i = 1:nROIsInB % mvB{i}.tSeries(1249:1296,:) = []; % end

% as a figure, lets just make one plot comparing all the ROIS % size should be numROIs subjA x numROIs in subjB % each entry is the % of voxels in subjA ROI that were most highly % correlated with the ROI in subjB % make matrix

ROImix = zeros(size(mvA,2),size(mvB,2));

% for each ROI in subject A for m = 1:size(mvA,2) % for each voxel in subject A

   for i = 1:size(mvA{m}.tSeries, 2) %number of voxels in ROI in A

% vector for the correlations

     corrs=[];

% for each ROI in subject B

     for j = 1:size(mvB,2) %number of ROIs

% for each voxel in subject B

         for k = 1:size(mvB{j}.tSeries, 2); %number of voxels in the ROI in B
             %get correlation between the time series of the two voxels
             r=corrcoef(mvA{m}.tSeries(:,i), mvB{j}.tSeries(:,k)); %records the correlations in a

% 2D matrix with each ROI's set of correlations in a row % so corrs is j rois by k voxels, where k is the largest number of voxels % in any of the rois

             corrs(j,k) = r(2);
             
             
             %for glm betas correlation

% r=corrcoef(mvA{m}.glm.betas(:,i),mvB{j}.glm.betas(:,k)); % corrs(j,k) = r;

         end
     end   

% then for each voxel in ROIA find the maximum correlation in Brain b

     maxCorr(i) = max(max(corrs)); % get the max correlation from all ROIs

% find the name of the ROI

     [row,col] = find(corrs == max(max(corrs))); % finds the row that the correlation was in and then records the name from that ROI

% if there is more than one ROI with the same max correlation pick the % second.... probably doesn't happen much.

     if length(row) > 1
         row = row(2);
     end
     
     maxCorrLocnum(i)=row;

% get the name of the roi

     maxCorrLoc(i) = cellstr(mvB{row}.roi.name); %record from where that coef came from
   end


   % make a histogram showing how many times each roi in subject b was chosen
   % as a match
   % get names of chosen rois
    a=unique(maxCorrLoc);
   %  get count of each from maxCorrLoc
   roipicks=[];
   for i = 1:length(a)
    roipicks=[roipicks,size(find(strcmp(cellstr(maxCorrLoc),a(i))),2)];
   end
   % bar plot

% figure(gcf+1); % bar(roipicks); % set(gca,'XTickLabel',a); % Title(mvA{m}.roi.name);


   %get a sense for the kind of correlations we are getting

% figure(gcf+1); hist(maxCorr);


% now for this ROI in A we want to get the %of times each ROI in B was % picked % for each ROI in subjB, find the number of times it was picked % normalized by number of voxels in ROI from subjA

   for z = 1:size(mvB,2)
       ROImix(m,z) = size(find(maxCorrLocnum==z),2);%/size(mvA{m}.tSeries, 2);
   end
   

end


% normalize ROImix denom=sum(ROImix,2); rdenom = repmat(denom,1,size(ROImix,2)); nROImix = ROImix./rdenom; % plot figure figure(gcf+1); imagesc(nROImix); set(gca,'YTick',1:size(roisA,2)); set(gca,'XTick',1:size(roisB,2)); set(gca,'YTickLabel',roisA); set(gca,'XTickLabel',roisB); %set title and axis labels colorbar;