MATLAB script for correlation maps: Difference between revisions

From Psych 221 Image Systems Engineering
Jump to navigation Jump to search
imported>Psych 204
Created page with '% 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 f…'
 
imported>Psych 204
No edit summary
 
(One intermediate revision by the same user not shown)
Line 1: Line 1:
% want to make a parameter map on subjectB's brain that shows the
  % 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
  % correlation in each voxel for a particular voxel chosen from subjectA's
% brain
  % brain


% load some data:
  % load some data:
clear all
  clear all
format long g
  format long g




%set directory for the two subjects to be compared
  %set directory for the two subjects to be compared
% make sure you run this file from cjimmy
  % make sure you run this file from cjimmy
if ispc
  if ispc
     subjADir='Y:\New_Localizers\dr021108-1p5mm-3mm';
     subjADir='Y:\New_Localizers\dr021108-1p5mm-3mm';
     subjBDir='Y:\New_Localizers\kw011708-1p5mm-3mm';
     subjBDir='Y:\New_Localizers\kw011708-1p5mm-3mm';
else
  else
     %    kevin subjects
     %    kevin subjects
     %    subjADir='/biac2/kgs/projects/New_Localizers/dr021108-1p5mm-3mm';
     %    subjADir='/biac2/kgs/projects/New_Localizers/dr021108-1p5mm-3mm';
Line 23: Line 23:


     % nathan's prosopagnosic subjects
     % nathan's prosopagnosic subjects
%    subjADir='/biac2/kgs/projects/Prosopagnosia/fMRIData/112909whFMRI';
  %    subjADir='/biac2/kgs/projects/Prosopagnosia/fMRIData/112909whFMRI';
% subjADir='/biac2/kgs/projects/Prosopagnosia/fMRIData/020209_mr_fmri';
  % subjADir='/biac2/kgs/projects/Prosopagnosia/fMRIData/020209_mr_fmri';
%    subjBDir='/biac2/kgs/projects/Prosopagnosia/fMRIData/52309AF';
  %    subjBDir='/biac2/kgs/projects/Prosopagnosia/fMRIData/52309AF';
%  subjADir =
  %  subjADir =
%  '/biac2/kgs/projects/Kids/fmri/localizer/adult_jc_27yo_052408';
  %  '/biac2/kgs/projects/Kids/fmri/localizer/adult_jc_27yo_052408';
subjADir = '/biac2/kgs/projects/Kids/fmri/localizer/adult_al_22yo_051108/';
    subjADir = '/biac2/kgs/projects/Kids/fmri/localizer/adult_al_22yo_051108/';
% subjBDir =
  % subjBDir =
% '/biac2/kgs/projects/Kids/fmri/localizer/adult_kw_25yo_090308';
  % '/biac2/kgs/projects/Kids/fmri/localizer/adult_kw_25yo_090308';
subjBDir = '/biac2/kgs/projects/Kids/fmri/localizer/adult_dy_25yo_041908';
  subjBDir = '/biac2/kgs/projects/Kids/fmri/localizer/adult_dy_25yo_041908';
      
      
end  
  end  


% %choose an ROI in subject A and B
  % %choose an ROI in subject A and B
% waldemar
  % roiA='lPPA';
% roiA='lPPA';
  % roiA='rFFA';
% roiA='rFFA';
  % alana
% alana
  % roiB='AllRois_Combined';
% roiB='AllRois_Combined';




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


  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];


  %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);


  % .. 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
  % get the roi data from subject 1
%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
  %then we grab the multivoxel data from the ROI
%we will probably want to set these parameters ourselves so that we know
  %this should give us mvA which is a struct with lots of useful stuff
%what happened
  %in particular
%this gives us
  %mvA.trials has information about conditions and order of stimuli etc
%mvA.glm has the betas, the design matrix and so on
  %mvA.coord has the coordinates for each voxel in the roi should be 3 x
%for a block design there will be as many betas as there are conditions
  %numVoxels
%for an event related design there will be timepoints x conditions x voxels
  %mvA.tSeries has the time series for each voxel in % signal change
%betas
  %mvA.roi had roi info
%    mvA{i} = mv_applyGlm(mvA{i});
  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
  % now do the same for subject B
%instead invoke hidden inplane with multiple rois.  then loop through the
  % set scan numbers for B
%rois when comparing to subjectA
  % scan=[1 2 3 4 5];
%go to subject b directory
  scan=[1 2 3 4];
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});


  % .. 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);
  % can get the size of the brain map in the view for a particular scan by
% which should be x,y, slice
  % doing


% when we use mv_init on our rois we end up with coordinates in the field
  sz=dataSize(hiB,1);
% mvA.roi.coords that have a particular numbering convention.  we want to
  % which should be x,y, slice
% 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
  % when we use mv_init on our rois we end up with coordinates in the field
voxelInA = 10;
  % 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


% then we compute the correlation with each voxel in subject B's roi
  % 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));
  %        vector for the correlations
coordinates=zeros(3,size(mvB.roi.coords,1));
  corrs=zeros(size(mvA.tSeries,2),size(mvB.tSeries,2));
for j = 1:size(mvA.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
     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
         %get correlation between the time series of the two voxels
Line 146: Line 144:
     end
     end
     %
     %
end
  end


% get an average correlation
  % get an average correlation


meancorr = mean(corrs,1);
  meancorr = mean(corrs,1);


% function [roiInd, coords] = roiIndices(view, coords, preserveCoords);
  % function [roiInd, coords] = roiIndices(view, coords, preserveCoords);




% in this case we think it should be
  % in this case we think it should be
[roiInd, coords] = roiIndices(hiB, mvB.roi.coords, 1);
  [roiInd, coords] = roiIndices(hiB, mvB.roi.coords, 1);


% [oldcoords  viewcoords] = roiIndices(view, ourroi, 0)
  % [oldcoords  viewcoords] = roiIndices(view, ourroi, 0)
% with the coords variable containing the mapping to the inplane.
  % with the coords variable containing the mapping to the inplane.


% then we propose to make an array of the same size as the view
  % then we propose to make an array of the same size as the view
ourmap = zeros(sz);
  ourmap = zeros(sz);


% then we will assign our correlations for each voxel to its location in
  % then we will assign our correlations for each voxel to its location in
% the matrix
  % the matrix
for i = 1:size(coords,2)
  for i = 1:size(coords,2)
    ourmap(coords(1,i),coords(2,i),coords(3,i))=meancorr(i);
      ourmap(coords(1,i),coords(2,i),coords(3,i))=meancorr(i);
end
  end


map = {ourmap,[]};
  map = {ourmap,[]};
mapName = 'cor_meanallFFA_KWall';
  mapName = 'cor_meanallFFA_KWall';
mapUnits = 'r';
  mapUnits = 'r';
co ={[],[]};
  co ={[],[]};


pth = ['Inplane/GLMs/' mapName];
  pth = ['Inplane/GLMs/' mapName];
save(pth, 'map', 'mapName', 'mapUnits');
    save(pth, 'map', 'mapName', 'mapUnits');
% save our parameter map  
  % save our parameter map  




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

Latest revision as of 09:24, 9 December 2009

  % 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.