{{{ 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 }}}