Diff for "FreesurferGoodies/transformRAStoNative" - MRC CBU Imaging Wiki
location: Diff for "FreesurferGoodies/transformRAStoNative"
Differences between revisions 3 and 4
Revision 3 as of 2012-07-05 10:30:27
Size: 15052
Editor: IanCharest
Comment:
Revision 4 as of 2013-03-07 21:23:59
Size: 15052
Editor: localhost
Comment: converted to 1.6 markup
No differences found!

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

CbuImaging: FreesurferGoodies/transformRAStoNative (last edited 2013-03-07 21:23:59 by localhost)