function labelSizes=transformRAStoNative
%
% based on the spm8_RAS2SPMcoords (Linda Henriksson 2009)
% based on spm8_label2ROI.m (authors: Linda H & Tom)
%
% also using the resizeROI tools developed by Niko Kriegeskorte.
%
% Ian Charest - MRC-CBU - March 2012
functionPath='/imaging/ic01/iRSA_fMRI/functions/';addpath(genpath(functionPath));
mridataPath='/imaging/ic01/iRSA_fMRI/mridata_sessionSpecific/';
% add the freesurfer matlab path in the path
freesurferRoot='/imaging/local/linux/freesurfer_4.3.0/freesurfer/matlab';addpath(genpath(freesurferRoot));
fsPath='/imaging/ic01/iRSA_fMRI/subjects/';
subjectNames={...
'CBXXXXXX','CBUYYYYYY',...
};
ROIidentifier={'hIT'...
};
hemi={'lh'...
'rh'};
hemiIdentifier={'left'...
'right'};
nSubjects=length(subjectNames);
% generate multiple ROI sizes
nVoxels=round(logspace(log10(20),log10(800),20));
labelSizes=nan(2,nSubjects);
for subI=1:nSubjects
%transform the RAS coordinates back to spm in the subject's native
%space.
thisSubject = subjectNames{subI};
%define the path in which the register.dat file is in the freesurfer
%subject folder.
thisRegistrationPath = [fsPath,thisSubject,'/new_t_maps'];
thisStats = [fsPath,thisSubject,'/new_t_maps/spmT_0004.hdr']; %all vs baseline contrast t-map
thisSubjectMaskLocation=[mridataPath,thisSubject,'/masks'];
if exist(thisSubjectMaskLocation)~=7
mkdir(thisSubjectMaskLocation)
end
for hemiI=1:2
% automatically labeled hIT in freesurfer
thisLabel=[fsPath,thisSubject,'/label/',hemi{hemiI},'.',ROIidentifier{1},'.label'];
% get from RAS (freesurfer) to SPM
vox=spm8_RAS2SPMcoords(thisLabel,thisRegistrationPath,thisStats);
% get the all vs. baseline t-map (the one not fiddled about to match freesurfer subject's anat)
labelSizes(hemiI,subI)=size(vox,2);
thisSPMt_map=[mridataPath,thisSubject,'/faceLocStats/spmT_0001.hdr'];
P = spm_vol(thisSPMt_map);
Y = spm_read_vols(P);
hITvox = vox';
% with one subject, we have a problem of having 0 as index on 2nd
% dim.
hITy=hITvox(:,2);
hITvox=hITvox(logical(hITy~=0),:);
hITINDs=sub2ind(size(Y),hITvox(:,1),hITvox(:,2),hITvox(:,3));
wholeMaskhIT=zeros(size(Y));
wholeMaskhIT(hITINDs)=1;
newRoiMapMetadataStruct = P;
newRoiMapMetadataStruct.fname = fullfile([thisSubjectMaskLocation,filesep,hemiIdentifier{hemiI},'_',ROIidentifier{1},'_automaticLabel.img']);
spm_write_vol(newRoiMapMetadataStruct, wholeMaskhIT);
newRoiMapMetadataStruct = P;
for voxelsI=1:length(nVoxels)
% for hIT
newhITRoi = resizeRoi(single(hITvox),Y, nVoxels(voxelsI),wholeMaskhIT);
newhITRoimask = roi2mask(newhITRoi,size(Y));
% update metaData structure specifying hIT
newRoiMapMetadataStruct.fname = fullfile([thisSubjectMaskLocation,filesep,hemiIdentifier{hemiI},ROIidentifier{1},'_',num2str(nVoxels(voxelsI)),'.img']);
newRoiMapMetadataStruct.descrip = ['binary ROI mask for ' ROIidentifier{1}];
newRoiMapMetadataStruct.dim = size(newhITRoimask );
spm_write_vol(newRoiMapMetadataStruct,newhITRoimask );
clear newhITRoi newhITRoimask
end
end
end
%store how many voxels were in the whole hIT masks.
save('/imaging/ic01/iRSA_fMRI/mridata_sessionSpecific/Details/iRSA_labelSizes_hIT.mat','labelSizes')
end
function newRoi=resizeRoi(roi, map, nVox, mask)
% USAGE
% newRoi=resizeRoi(roi, map, nVox[, mask])
%
% FUNCTION
% redefine a roi to make it conform to the size
% given by nVox. the new roi is defined by a region
% growing process, which (1) is seeded at the voxel
% that has the maximal statistical parameter within
% the passed statistical map and (2) is prioritized
% by the map's values.
%
% ARGUMENTS
% roi (r)egion (o)f (i)nterest
% a matrix of voxel positions
% each row contains ONE-BASED coordinates (x, y, z) of a voxel.
%
% map a 3D statistical-parameter map
% the map must match the volume, relative to which
% the roi-voxel coords are specified in roi.
%
% nVox number of voxels the resized roi is to have
%
% [mask] optional binary mask. if present, the region growing is
% restricted to the nonzero entries of it.
% PARAMETERS
%maxNlayers=2; %defines how many complete layers are maximally added
if ~exist('mask','var') || (exist('mask','var') && isempty(mask))
mask=ones(size(map));
end
map(~mask)=min(map(:));
% DEFINE THE VOLUME
vol=zeros(size(map));
% FIND THE SEED (A MAXIMAL MAP VALUE IN ROI&mask)
mapINDs=sub2ind(size(vol),roi(:,1),roi(:,2),roi(:,3)); %single indices to MAP specifying voxels in the roi
roimap=map(mapINDs); %column vector of statistical-map subset for the roi
[roimax,roimax_roimapIND]=max(roimap); %the maximal statistical map value in the roi and its index within roimap
seed_mapIND=mapINDs(roimax_roimapIND); %seed index within map
newRoi=[];
if nVox==0; return; end;
if mask(seed_mapIND)
vol(seed_mapIND)=1;
[x,y,z]=ind2sub(size(vol),seed_mapIND);
newRoi=[newRoi;[x,y,z]];
end
if nVox==1 || isempty(newRoi)
return;
end
% GROW THE REGION
for i=2:nVox
% DEFINE THE FRINGE
cFringe=vol;
[ivolx,ivoly,ivolz]=ind2sub(size(vol),find(vol));
superset=[ivolx-1,ivoly,ivolz;
ivolx+1,ivoly,ivolz;
ivolx,ivoly-1,ivolz;
ivolx,ivoly+1,ivolz;
ivolx,ivoly,ivolz-1;
ivolx,ivoly,ivolz+1];
% exclude out-of-volume voxels
outgrowths = superset(:,1)<1 | superset(:,2)<1 | superset(:,3)<1 | ...
superset(:,1)>size(vol,1) | superset(:,2)>size(vol,2) | superset(:,3)>size(vol,3);
superset(find(outgrowths),:)=[];
% draw the layer (excluding multiply defined voxels)
cFringe(sub2ind(size(vol),superset(:,1),superset(:,2),superset(:,3)))=1;
cFringe=cFringe&mask;
cFringe=cFringe-vol;
if size(find(cFringe),1)==0
break; % exit the loop (possible cause of empty fringe: the whole volume is full)
end
% FIND A MAXIMAL-MAP-VALUE FRINGE VOXEL...
mapINDs=find(cFringe); %single indices to MAP specifying voxels in the fringe
fringemap=map(mapINDs); %column vector of statistical-map subset for the fringe
[fringemax,fringemax_fringemapIND]=max(fringemap); %the maximal statistical map value in the roi and its index within roimap
fringemax_mapIND=mapINDs(fringemax_fringemapIND); %seed index within map
% ...INCLUDE IT
vol(fringemax_mapIND)=1;
[x,y,z]=ind2sub(size(vol),fringemax_mapIND);
newRoi=[newRoi;[x,y,z]];
end
end
% roi2mask is a function with two arguments:
% ARGUMENTS:
% roi: this is a 3xn matrix where each row contains the coordinates for
% a point inside the roi
% volSize_vox: this is a 1x3 vector containing the dimensions of the
% scanned volume. E.g., [32 64 64]
% RETURNS:
% mask: a volume of size volSize_vox which is all 0s, except for the
% points indicated by roi, which are 1s.
function mask=roi2mask(roi,volSize_vox)
roi_INDs=sub2ind(volSize_vox,roi(:,1),roi(:,2),roi(:,3)); %single indices to MAP specifying voxels in the roi
mask=false(volSize_vox);
mask(roi_INDs)=true;
end
function [vox mm prob]=spm8_RAS2SPMcoords(rasdata,fspath,forighdr)
% usage: [vox,mm]=spm8_togglebetweencoords(rasdata,fspath,forighdr)
% e.g. [vox mm]=spm8_RAS2SPMcoords([4 -68 4],'/mnt/vision/fssubjects/Linda_Henriksson/functional/0000','/opt/dicom_images/Linda_Henriksson/Linda_Henriksson/0000/epi8/meanalinda8_170.hdr')
% e.g. [vox mm]=spm8_RAS2SPMcoords('/mnt/vision/fssubjects/Linda_Henriksson/label/phaseproject/rh_V1.label','/mnt/vision/fssubjects/Linda_Henriksson/functional/0000','/opt/dicom_images/Linda_Henriksson/Linda_Henriksson/0000/epi8/meanalinda8_170.hdr')
%
% Transfer coordinates from freesurfer labeled area OR individual RAS
% coordinates to functional spm-coordinates using aras2fcoords.m.
%
% rasdata = RAS coordinates voxels or a Freesurfer label file
% e.g. one voxel: [4 -68 4]
% OR
% more than one voxels in rows: [4 -68 4; 1 -67 3] or in columns [4 1; -68 -67; 4 3]
% OR
% '/mnt/vision/fssubjects/Linda_Henriksson/label/phaseproject/rh_V1.label'
%
% fspath = path to freesurfer functional folder, where files register.dat
% and v2r.dat ((transformation between voxels and ras-coordinates - defined
% by user) can be found,
% e.g. '/mnt/vision/fssubjects/Linda_Henriksson/functional/0000'
%
% forighdr = functional data hdr file with SPM.mat coordinate
% transformations,
% e.g. '/opt/dicom_images/Linda_Henriksson/Linda_Henriksson/0000/epi8/meanalinda8_170.hdr'
%
% based on spm8_label2ROI.m (authors: Linda H & Tom)
% Linda H, November 2009
% Ian Charest. Added the probability as an output (using Hinds accurate
% probabilistic estimate of V1 we get a label of probabilities).
if ischar(rasdata)
disp(sprintf('\n Input is a label file'))
[label r a s t]=textread(rasdata,'%d %f %f %f %f','headerlines',2);
vox=[]; mm=[]; prob=[];
for i=1:length(r)
[fvox,fmm]=spm8_aras2fcoords([r(i) a(i) s(i)]',fspath,forighdr);
vox=[vox,fvox]; mm = [mm,fmm];
prob=[prob t(i)];
end
% round voxels to integer values
vox = round(vox);
% exclude overlapping voxels
[nonsense,ind]=unique(vox','rows');
vox = vox(:,ind);
mm = mm(:,ind);
prob = prob(ind);
disp(sprintf('\n output: %d voxels \n',size(vox,2)))
else
disp('\n Input is RAS coordinate vector')
vox=[]; mm=[];
[r,c]=size(rasdata);
if r==3 || c ==3
if r ==3 && c ==3
disp('\n input: 3 voxel(s) \n function assumes that different colums are coordinates for different voxels')
elseif r ==3 && c ~=3
disp(sprintf('\n input: %d voxel(s)',c))
elseif r ~=3 && c==3
disp(sprintf('\n input: %d voxel(s)',r))
rasdata = rasdata';
end
else
error(sprintf('\n INTERRUPTED, because: \n RAS coordinates input must be a x-by-3 (or 3-by-x) [R A S] vector,\n'))
end
for i=1:(size(rasdata,2))
[fvox,fmm]=spm8_aras2fcoords([rasdata(1,i) rasdata(2,i) rasdata(3,i)]',fspath,forighdr);
vox=[vox,fvox]; mm = [mm,fmm];
end
% round voxels to integer values
vox = round(vox);
% exclude overlapping voxels
[nonsense,ind]=unique(vox','rows');
vox = vox(:,ind);
mm = mm(:,ind);
disp(sprintf('\n output: %d voxel(s) \n',size(vox,2)))
end
end
function [fvoxelcoords,fmmcoords]=spm8_aras2fcoords(rascoords,fspath,forigmat)
% function [fvoxelcoords,fmmcoords] = aras2fcoords(rascoords,fspath,forigmat)
%
% Earlier version: ras2fmm, Important modification: register.dat included
%
% Coordinate transforms from Freesurfer RAS (Right-Anterior-Superior)
% coordinates to original SPM functional data voxel and mm coordinates.
%
% Inputs:
% rascoords = [R A S]' coordinates of a point in freesurfer surface,
% eg. [54 -68 20]'
% fspath = path to freesurfer data with register.dat and v2r.dat info,
% eg. '/remote/ws0_opt2/fssubjects/Linda_Henriksson/functional/2957_phaseadapt'
% forigmat = data mat file with SPM.mat coordinate transformations,
% eg. '/opt/dicom_images/Linda_Henriksson/Linda_Henriksson/2957/epi8/meanalinda8_272.mat'
%
% Outputs:
% fvoxelcoords = [X Y Z] SPM voxel coordinates for the input point
% fmmcoords = [X Y Z] SPM mm coordinates for the input point
%
% Linda Henriksson, May 2008
%
% lh updated for spm8, June 2009
%
% lh rounding error of vox -> mm corrected, Nov 2011
cd(fspath)
[h,l]=size(rascoords);
% check that all inputs are available
if h*l~=3
error(sprintf('\n TASK INTERRUPTED, because: \n RAS coordinates input must be a 1-by-3 (or 3-by-1) [R A S] vector,\n'))
elseif l==3
rascoords=rascoords';
end
if ~exist('register.dat')
error(sprintf('\n TASK INTERRUPTED, because: \n there is no register.dat file in path \n %s \n please coregister the functional data to freesurfer anatomical in SPM and create register.dat file in FREESURFER with command like \n eg. tkregister2 --mov meanalinda8_170.img --s Linda_Henriksson --regheader --noedit --reg register.dat \n',fspath))
end
if ~exist('v2r.dat')
error(sprintf('\n TASK INTERRUPTED, because: \n there is no v2r.dat file in path \n %s \n please create a vox2ras_tkr matrix for the functional data in FREESURFER with command \n mri_info meana*.img --vox2ras-tkr --o v2r.dat',fspath))
end
if ~exist('vox2ras_0to1.m')
error(sprintf('\n TASK INTERRUPTED, because: \n You need to have freesurfer matlab functions in the path!! \n eg. addpath /opt2/freesurfer/matlab \n'))
end
% 1. coordinate transform
regmat=textread('register.dat','',4,'headerlines',4);
% register.dat contains a matrix with rotation and translation info between
% anatomical RAS and functional "RAS" coordinates
% 2. coordinate transform
vox2ras_tkr=load('v2r.dat');
% v2r.dat contains a matrix with rotation and translation info between functional
% data "RAS" coordinates and voxel coordinates
% | -dC 0 0 +Nc/2 |
% | 0 0 +dS -Ns/2 | , voxel sizes dC, dR, dS
% | 0 -dR 0 +Nr/2 | and Nc = width, Nr = height, Ns =depth
% | 0 0 0 1 |
% eg. with imaging parameters FOV = 180 mm x 180 mm, slice thickness = 2.8
% and number of slices = 26 the matrix looks like:
% | -2.8125 0 0 180/2 |
% | 0 0 2.8 -26*2.8/2 |
% | 0 -2.8125 0 180/2 |
% | 0 0 0 1 |
% with the info how the origo of RAS coordinates is translated to the origo
% of voxel coordinates and how the data is rotated to match correct points.
% The definition of voxel coordinates is different in FREESURFER and SPM,
% therefore we use this freesurfer matlab function to convert "0-based vox2ras matrix to 1-based matrix"
vox2ras_tkr=vox2ras_0to1(vox2ras_tkr);
% and then combine the two transformation matrices
regmat(:,5)=[];
mat1=inv(vox2ras_tkr)*regmat;
% and get the functional voxel coordinates for rascoords
fvoxelcoords=mat1(1:3,1:3)*rascoords+mat1(1:3,4);
% 3. coordinate transform
% load coordinate transformations between functional mm and voxel coordinates in SPM.mat space
% ---SPM8
P = spm_vol(forigmat);
M = P.mat;
% SPM2: load(sprintf('%s',forigmat));
% and get the functional mm coordinates for rascoords
% rounding of vox coordinates added (11/2011) to give more accurate mm coordinates
fmmcoords=M(1:3,1:3)*round(fvoxelcoords)+M(1:3,4);
% -----
end
location: FreesurferGoodies / transformRAStoNative