FreesurferGoodies/transformRAStoNative - MRC CBU Imaging Wiki

Upload page content

You can upload content for the page named below. If you change the page name, you can also upload content for another page. If the page name is empty, we derive the page name from the file name.

File to load page content from
Page name
Comment
Type the missing letters from: Lodon is th capial of nglnd

location: FreesurferGoodies / transformRAStoNative

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