DTI Coregistration to Anatomical Scans

Here I document my workflow for coregistration of DTI data to anatomical scans using FSL. Feel free to make annotations and comments.

I have one T2 weighted image (volume 0) and 64 diffusion weighted images (DWI) (volumes 1-64) in a 4D-NIfTI file with 128x128x96 1.5mm isovoxels.
Also a T1 image with a good brain extraction (gray matter) with 256^3 1mm isovoxels and a CSF mask.

Script

Here is my script to automate the coregistration as described above.

#!/bin/bash
# Mirko Windhoff, 2009
# $Id: dti_coregistration 213 2009-07-20 15:42:51Z mirko $
# cd /your/dti/folder/
# 5 input files needed, use relative paths, be sure, that there are no other files or directories in the working directory!
T1_REF=T1_ISO_homogen_ACPC_peeled_GM.nii.gz # brain extracted reference image (3D-NIfTI)
CSF_REF_MASK=T1_ISO_homogen_ACPC_CSF_Mask.nii.gz # csf reference image mask (3D-NIfTI)
# CSF_REF= # csf extracted reference image (3D-NIfTI), if no mask available (see below)
INPUT=data.nii.gz # T2+DWI in a 4D-NIfTI, not yet coregistered
BVECS=bvecs
BVALS=bvals
 
CSF_REF_MASK=T1_ISO_homogen_ACPC_csf_mask.nii.gz
# csf extracted reference image (3D-NIfTI), if no csf mask available:
# CSF_REF=
# echo "Creating csf mask..."
# CSF_REF_MASK=ref_csf_mask.nii.gz
# bet $CSF_REF $CSF_REF_MASK -m
# CSF reference image mask (3D-NIfTI):
INPUT=data.nii.gz # T2+DWI in a 4D-NIfTI, not yet coregistered
BVECS=bvecs
BVALS=bvals
 
timer=`date +%s`
echo "DTI Coregistration started..."
# fslsplit to get each volume in a separate NIfTI-file
echo "Splitting t2 and dwi..."
mkdir data_split
fslsplit $INPUT data_split/vol
# bet2 volume 0, to remove signal of skin (otherwise eddy_correct could take the skin as the outer border of the brain, in my case fractional intensity was set to 0.22).
echo "Betting t2..."
mkdir bet_t2
bet data_split/vol0000 bet_t2/vol0000_bet.nii.gz -f 0.2
echo "Flirting t2 to t1..."
mkdir data_coreg
flirt -in bet_t2/vol0000_bet.nii.gz -ref $T1_REF -cost normmi -omat t2_to_t1.mat
echo "Flirting dwi to t2 with correlation ratio and apply t2_to_t1 transformation..."
for filename in `imglob data_split/*.*`; do
    matfilename="${filename%.*}.mat"
    matfilename2="${filename%.*}_flirted.mat"
    flirt -in $filename -ref bet_t2/vol0000_bet.nii.gz -nosearch -omat $matfilename -paddingsize 1
    convert_xfm -omat $matfilename2 -concat t2_to_t1.mat $matfilename
    flirt -in $filename -ref $T1_REF -applyxfm -init $matfilename2 -out data_coreg/`basename "$filename"`
done
echo "Merging coregistrated t2 and dwi..."
fslmerge -t data_coreg.nii.gz `imglob data_coreg/*.*`
echo "Flirting bvecs with matlab..."
matlab -nojvm -wait -r "flirt_bvecs('data_split', '$BVECS', 'bvecs_flirt'); exit;"
echo "Dtifitting..."
dtifit --data=data_coreg.nii.gz --out=dti --mask=$CSF_REF_MASK --bvecs=bvecs_flirt --bvals=$BVALS --save_tensor
seconds=$[`date +%s`-timer]
echo "done. Calculation time $seconds seconds."
 
# echo "Removing intermediate processing files..."
# rm -rf data_split bet_t2 data_coreg

The flirt_bvecs Matlab function that is called by the script above looks like this:

function flirt_bvecs(flirt_matrices_directory, bvecs_filename, transformed_bvecs_output_filename)
  % Reads the flirt transformation matrix and applies it to the bvecs taken
  % from the bvecs file. The transformed bvecs are written to the bvecs
  % output file.
  % USAGE:
  % FLIRT_BVECS(flirt_matrix_filename, bvecs_filename, transformed_bvecs_output_filename)
  %
  % $Id: flirt_bvecs.m 219 2009-07-21 09:38:16Z mirko $
 
  % read the bvecs
  fd=fopen(bvecs_filename,'rt');
  for i=1:3
    b(i,:)=sscanf(fgetl(fd), '%f', [1 inf]);
  end;
  fclose(fd);
  % read the transformation matrices
  files=dir(fullfile(flirt_matrices_directory,'*_flirted.mat'));
  if isempty(files), error(['No files matching *_flirted.mat in ', flirt_matrices_directory, ' found.']); end;
  b_trafo=zeros(size(b));
  for i=2:size(files,1)
      fd=fopen(fullfile(flirt_matrices_directory,files(i).name),'r');
      trafo=zeros(4,4);
      trafo(1,1:4)=str2num(fgetl(fd));
      trafo(2,1:4)=str2num(fgetl(fd));
      trafo(3,1:4)=str2num(fgetl(fd));
      trafo(4,1:4)=str2num(fgetl(fd));
      fclose(fd);
      % remove translation components, bvecs are translation invariant
      trafo=trafo(1:3,1:3);
      % do some norm with the rotation matrix
      trafo=real((trafo*trafo')^(1/2)*trafo);
      % transform bvecs
      b_trafo(:,i)=trafo*b(:,i);
      b_trafo(:,i)=b_trafo(:,i)./norm(b_trafo(:,i)); % norm to length 1
  end;
  % write the transformed bvecs
  if exist(transformed_bvecs_output_filename,'file'), error([transformed_bvecs_output_filename ' already exists.']); end;
  fd=fopen(transformed_bvecs_output_filename,'wt');
  for i=1:3, fprintf(fd, '%1.5f ', b_trafo(i,:)); fprintf(fd, '\n'); end;
  fclose(fd);
end
 
 
dti_coregistration.txt ยท Last modified: 19.05.2010 22:26 by 95.208.115.158
Recent changes RSS feed Creative Commons License Valid XHTML 1.0 Valid CSS Driven by DokuWiki
Drupal Garland Theme for Dokuwiki