MATLAB script for correlation maps

From Psych 221 Image Systems Engineering
Jump to navigation Jump to search
  % want to make a parameter map on subjectB's brain that shows the
  % correlation in each voxel for a particular voxel chosen from subjectA's
  % brain
  % load some data:
  clear all
  format long g


  %set directory for the two subjects to be compared
  % make sure you run this file from cjimmy
  if ispc
   subjADir='Y:\New_Localizers\dr021108-1p5mm-3mm';
   subjBDir='Y:\New_Localizers\kw011708-1p5mm-3mm';
  else
   %     kevin subjects
   %     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';
   % location of code
   home = '/biac2/kgs/projects/SameDifferent/'
   % nathan's prosopagnosic subjects
  %     subjADir='/biac2/kgs/projects/Prosopagnosia/fMRIData/112909whFMRI';
  % subjADir='/biac2/kgs/projects/Prosopagnosia/fMRIData/020209_mr_fmri';
  %     subjBDir='/biac2/kgs/projects/Prosopagnosia/fMRIData/52309AF';
  %  subjADir =
  %  '/biac2/kgs/projects/Kids/fmri/localizer/adult_jc_27yo_052408';
   subjADir = '/biac2/kgs/projects/Kids/fmri/localizer/adult_al_22yo_051108/';
  % subjBDir =
  % '/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 and B
  % roiA='lPPA';
  % roiA='rFFA';
  % alana
  % roiB='AllRois_Combined';


  % roiA='lPPA_IOvAC_p3';
  roiA='lFFA_MBvAC_p3';
  roiB = 'allROIS';


  %want to get the motion corrected data type which is usually 3 but might be
  %anything depending on the experimenter and experiments.
  dt = 3; %'MotionComp_RefScan5';
  %then pick the scans we will use.
  % mtloc,eccbias1,eccbias2,4xretinotopy,objloc1,objloc2
  % so to get our scans in order we have
  % scan = [8 9 1 2 3];
  scan = [1 2 3 4];


  % .. for subject A brain
  %go to subjects data directory
  cd(subjADir);
  %initialize our hidden inplane set to the right data type,scan(s), and roi
  % hiA = initHiddenInplane(dt, scan, roisA);
  hiA = initHiddenInplane(dt, scan);


  % get the roi data from subject 1
  %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
  mvA = mv_init(hiA, roiA, scan, dt,1);   
  %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
  %     mvA{i} = mv_applyGlm(mvA{i});   


  % now do the same for subject B
  % set scan numbers for B
  % scan=[1 2 3 4 5];
  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);
  % get dataTpe for subjectB
  dt=3;
  %load data with all rois
  hiB = initHiddenInplane(dt,scan,roiB);
  mvB = mv_init(hiB, roiB, scan, dt, 1);
  %     mvB{i} = mv_applyGlm(mvB{i});



  % can get the size of the brain map in the view for a particular scan by
  % doing
  sz=dataSize(hiB,1);
  % which should be x,y, slice
  % when we use mv_init on our rois we end up with coordinates in the field
  % mvA.roi.coords that have a particular numbering convention.  we want to
  % transform that into the space that the view uses which we hope is done
  % with the function
  % so we pick an arbitrary voxel from subject A's roi
  voxelInA = 10;
  % then we compute the correlation with each voxel in subject B's roi


  %         vector for the correlations
  corrs=zeros(size(mvA.tSeries,2),size(mvB.tSeries,2));
  coordinates=zeros(3,size(mvB.roi.coords,1));
  for j = 1:size(mvA.tSeries, 2)
   for k = 1:size(mvB.tSeries, 2); %number of voxels in the ROI in B
       %get correlation between the time series of the two voxels
       r=corrcoef(mvA.tSeries(:,j), mvB.tSeries(:,k)); %records the correlations in a
       % record correlations
       corrs(j,k) = r(2);
       %     in order to get the same number of timeseries as anatomical
       %     coordinates must set preserveCoord variable in mv_init to be 1
       % some time series will be duplicated but we don't care
       %     coordinates(:,k) = mvB.roi.coords(:,k)
   end
   %
  end
  % get an average correlation
  meancorr = mean(corrs,1);
  % function [roiInd, coords] = roiIndices(view, coords, preserveCoords);


  % in this case we think it should be
  [roiInd, coords] = roiIndices(hiB, mvB.roi.coords, 1);
  % [oldcoords  viewcoords] = roiIndices(view, ourroi, 0)
  % with the coords variable containing the mapping to the inplane.
  % then we propose to make an array of the same size as the view
  ourmap = zeros(sz);
  % then we will assign our correlations for each voxel to its location in
  % the matrix
  for i = 1:size(coords,2)
      ourmap(coords(1,i),coords(2,i),coords(3,i))=meancorr(i);
  end
  map = {ourmap,[]};
  mapName = 'cor_meanallFFA_KWall';
  mapUnits = 'r';
  co ={[],[]};
  pth = ['Inplane/GLMs/' mapName];
   save(pth, 'map', 'mapName', 'mapUnits');
  % save our parameter map 


  % -> load parameter map in mrvista and see what comes out.