#!/bin/sh # siena - Structural Image Evaluation, including Normalisation, of Atrophy # # Stephen Smith, FMRIB Image Analysis Group # # Copyright (C) 1999-2007 University of Oxford # # Part of FSL - FMRIB's Software Library # http://www.fmrib.ox.ac.uk/fsl # fsl@fmrib.ox.ac.uk # # Developed at FMRIB (Oxford Centre for Functional Magnetic Resonance # Imaging of the Brain), Department of Clinical Neurology, Oxford # University, Oxford, UK # # # LICENCE # # FMRIB Software Library, Release 5.0 (c) 2012, The University of # Oxford (the "Software") # # The Software remains the property of the University of Oxford ("the # University"). # # The Software is distributed "AS IS" under this Licence solely for # non-commercial use in the hope that it will be useful, but in order # that the University as a charitable foundation protects its assets for # the benefit of its educational and research purposes, the University # makes clear that no condition is made or to be implied, nor is any # warranty given or to be implied, as to the accuracy of the Software, # or that it will be suitable for any particular purpose or for use # under any specific conditions. Furthermore, the University disclaims # all responsibility for the use which is made of the Software. It # further disclaims any liability for the outcomes arising from using # the Software. # # The Licensee agrees to indemnify the University and hold the # University harmless from and against any and all claims, damages and # liabilities asserted by third parties (including claims for # negligence) which arise directly or indirectly from the use of the # Software or the sale of any products based on the Software. # # No part of the Software may be reproduced, modified, transmitted or # transferred in any form or by any means, electronic or mechanical, # without the express permission of the University. The permission of # the University is not required if the said reproduction, modification, # transmission or transference is done without financial return, the # conditions of this Licence are imposed upon the receiver of the # product, and all original and amended source code is included in any # transmitted product. You may be held legally responsible for any # copyright infringement that is caused or encouraged by your failure to # abide by these terms and conditions. # # You are not permitted under this Licence to use this Software # commercially. Use for which any financial return is received shall be # defined as commercial use, and includes (1) integration of all or part # of the source code or the Software into a product for sale or license # by or on behalf of Licensee to third parties or (2) use of the # Software or any derivative of it for research with the final aim of # developing software products for sale or license to a third party or # (3) use of the Software or any derivative of it for research with the # final aim of developing non-software products for sale or license to a # third party, or (4) use of the Software to provide any service to an # external organisation for which payment is received. If you are # interested in using the Software commercially, please contact Isis # Innovation Limited ("Isis"), the technology transfer company of the # University, to negotiate a licence. Contact details are: # innovation@isis.ox.ac.uk quoting reference DE/9564. export LC_NUMERIC=C Usage() { cat < [options] -o : set output directory (default output is _to__siena) -d : debug (don't delete intermediate files) -B "betopts" : options to pass to BET brain extraction (inside double-quotes), e.g. -B "-f 0.3" -2 : two-class segmentation (don't segment grey and white matter separately) -t2 : T2-weighted input image (default T1-weighted) -m : use standard-space masking as well as BET -t : ignore from t (mm) upwards in MNI152/Talairach space -b : ignore from b (mm) downwards in MNI152/Talairach space (b should probably be negative) -S "sienadiffopts" : options to pass to siena_diff timepoint differencing (inside double-quotes), e.g. -S "-s -i 20" -V : run ventricle analysis (VIENA) -v : optional user-supplied ventricle mask (default is ${FSLDIR}/data/standard/MNI152_T1_2mm_VentricleMask) EOF exit 1 } [ "$2" = "" ] && Usage [ `${FSLDIR}/bin/imtest $1` = 0 ] && Usage [ `${FSLDIR}/bin/imtest $2` = 0 ] && Usage Ao=`${FSLDIR}/bin/remove_ext $1` Bo=`${FSLDIR}/bin/remove_ext $2` thecommand="siena $@" shift 2 outdir=${Ao}_to_${Bo}_siena vienadir=${outdir}/viena debug=0 betopts="" sdopts="" sdo="-m" dostd=0 stdmask=0 stdroi="" origin3=37 # `fslval ${FSLDIR}/data/standard/MNI152_T1_2mm origin3` pixdim3=2 # `fslval ${FSLDIR}/data/standard/MNI152_T1_2mm pixdim3` Vmask="${FSLDIR}/data/standard/MNI152_T1_2mm_VentricleMask" do_viena=no while [ _$1 != _ ] ; do if [ $1 = -d ] ; then debug=1 ov=-ov shift elif [ $1 = -o ] ; then outdir=$2 vienadir=${outdir}/viena shift 2 elif [ $1 = -B ] ; then betopts=$2 shift 2 elif [ $1 = -S ] ; then sdopts=$2 shift 2 elif [ $1 = -2 ] ; then sdo="$sdo -2" shift elif [ $1 = -t2 ] ; then is_t2=" -s -t 2" shift elif [ $1 = -m ] ; then stdmask=1 dostd=1 shift elif [ $1 = -t ] ; then dostd=1 stdt=`echo $2 | sed 's/-/_/g'` stdt=`echo "10 k $stdt $pixdim3 / $origin3 + p" | dc -` stdroi="$stdroi -roi 0 1000000 0 1000000 0 $stdt 0 1" shift 2 elif [ $1 = -b ] ; then dostd=1 stdb=`echo $2 | sed 's/-/_/g'` stdb=`echo "10 k $stdb $pixdim3 / $origin3 + p" | dc -` stdroi="$stdroi -roi 0 1000000 0 1000000 $stdb 1000000 0 1" shift 2 elif [ $1 = -V ] ; then do_viena=yes shift elif [ $1 = -v ] ; then Vmask=$2 shift 2 else Usage fi done # ensure full path for ventriclemask (because of cd in scripts) if [ $do_viena = yes ] ; then if [ `${FSLDIR}/bin/imtest ${Vmask}` = 0 ] ; then echo "ERROR:: cannot find image ${Vmask}" Usage fi fi sdo="${sdo}${is_t2}" mkdir -p $outdir if [ $do_viena = yes ] ; then mkdir -p $vienadir fi ${FSLDIR}/bin/imcp $Ao ${outdir}/A ${FSLDIR}/bin/imcp $Bo ${outdir}/B cd $outdir A=A B=B echo 'FSL

SIENA Report

'${thecommand}'
' > report.html echo "-----------------------------------------------------------------------" > report.siena echo "" >> report.siena echo " SIENA - Structural Image Evaluation, using Normalisation, of Atrophy" >> report.siena echo " part of FSL www.fmrib.ox.ac.uk/fsl" >> report.siena echo " running longitudinal atrophy measurement: siena version 2.6" >> report.siena echo " siena $@" >> report.siena echo "" >> report.siena echo "---------- extract brain --------------------------------------------" >> report.siena ${FSLDIR}/bin/bet $A ${A}_brain -s -m $betopts >> report.siena ${FSLDIR}/bin/bet $B ${B}_brain -s -m $betopts >> report.siena ${FSLDIR}/bin/fslmaths ${A}_brain -sub `$FSLDIR/bin/fslstats ${A}_brain -p 0` -mas ${A}_brain_mask ${A}_brain -odt float ${FSLDIR}/bin/fslmaths ${B}_brain -sub `$FSLDIR/bin/fslstats ${B}_brain -p 0` -mas ${B}_brain_mask ${B}_brain -odt float ${FSLDIR}/bin/overlay 0 0 $A -a ${A}_brain 1 `${FSLDIR}/bin/fslstats ${A}_brain -P 95` ${A}_brain_skull 0.9 1.1 ${A}_grot ${FSLDIR}/bin/slicer ${A}_grot -a ${A}_bet.png ${FSLDIR}/bin/overlay 0 0 $B -a ${B}_brain 1 `${FSLDIR}/bin/fslstats ${B}_brain -P 95` ${B}_brain_skull 0.9 1.1 ${B}_grot ${FSLDIR}/bin/slicer ${B}_grot -a ${B}_bet.png ${FSLDIR}/bin/imrm ${A}_grot ${B}_grot echo "

BET brain extraction results

${Ao}

${Bo}
" >> report.html echo "" >> report.siena echo "---------- register brains and skulls -------------------------------" >> report.siena echo "(do not worry about histogram warnings)" >> report.siena ${FSLDIR}/bin/siena_flirt $A $B >> report.siena 2>&1 echo "


FLIRT A-to-B registration results

" >> report.html echo "" >> report.siena echo "---------- produce valid masks --------------------------------------" >> report.siena XDIM=`${FSLDIR}/bin/fslval $A dim1` ; XDIM=`echo "$XDIM 2 - p" | dc -` YDIM=`${FSLDIR}/bin/fslval $A dim2` ; YDIM=`echo "$YDIM 2 - p" | dc -` ZDIM=`${FSLDIR}/bin/fslval $A dim3` ; ZDIM=`echo "$ZDIM 2 - p" | dc -` ${FSLDIR}/bin/fslmaths ${A}_brain_mask -mul 0 -add 1 -roi 1 $XDIM 1 $YDIM 1 $ZDIM 0 1 ${A}_valid_mask XDIM=`${FSLDIR}/bin/fslval $B dim1` ; XDIM=`echo "$XDIM 2 - p" | dc -` YDIM=`${FSLDIR}/bin/fslval $B dim2` ; YDIM=`echo "$YDIM 2 - p" | dc -` ZDIM=`${FSLDIR}/bin/fslval $B dim3` ; ZDIM=`echo "$ZDIM 2 - p" | dc -` ${FSLDIR}/bin/fslmaths ${B}_brain_mask -mul 0 -add 1 -roi 1 $XDIM 1 $YDIM 1 $ZDIM 0 1 ${B}_valid_mask ${FSLDIR}/bin/flirt -in ${B}_valid_mask -ref $A -out ${B}_valid_mask_to_${A} -applyxfm -init ${B}_to_${A}.mat -paddingsize 0 ${FSLDIR}/bin/flirt -in ${A}_valid_mask -ref $B -out ${A}_valid_mask_to_${B} -applyxfm -init ${A}_to_${B}.mat -paddingsize 0 ${FSLDIR}/bin/fslmaths ${A}_valid_mask -mul ${B}_valid_mask_to_${A} ${A}_valid_mask_with_$B ${FSLDIR}/bin/fslmaths ${B}_valid_mask -mul ${A}_valid_mask_to_${B} ${B}_valid_mask_with_$A if [ $dostd = 1 ] ; then echo "" >> report.siena echo "---------- standard space masking ----------------------------------" >> report.siena ${FSLDIR}/bin/flirt -ref ${FSLDIR}/data/standard/MNI152_T1_2mm_brain -in ${A}_brain -o ${A}_to_std -omat ${A}_to_std.mat >> report.siena ${FSLDIR}/bin/flirt -ref ${FSLDIR}/data/standard/MNI152_T1_2mm_brain -in ${B}_brain -o ${B}_to_std -omat ${B}_to_std.mat >> report.siena ${FSLDIR}/bin/convert_xfm -inverse -omat ${A}_to_std_inv.mat ${A}_to_std.mat ${FSLDIR}/bin/convert_xfm -inverse -omat ${B}_to_std_inv.mat ${B}_to_std.mat ${FSLDIR}/bin/slicer ${A}_to_std ${FSLDIR}/data/standard/MNI152_T1_2mm_brain -a ${A}_to_std.png ${FSLDIR}/bin/slicer ${B}_to_std ${FSLDIR}/data/standard/MNI152_T1_2mm_brain -a ${B}_to_std.png echo "


FLIRT standard space registration results

${Ao}

${Bo}
" >> report.html ${FSLDIR}/bin/convert_xfm -concat ${B}_to_std_inv.mat -omat ${A}_to_${B}_tmp.mat ${A}_to_std.mat RMSDIFF=`${FSLDIR}/bin/rmsdiff ${A}_to_${B}.mat ${A}_to_${B}_tmp.mat $A | sed 's/\..*$/ /g'` # last part makes it integer echo "rmsdiff for standard space transform is $RMSDIFF mm" >> report.siena if [ $RMSDIFF -ge 10 ] ; then echo "Warning! Probably failed consistency check for standard-space registrations!" echo "Warning! Probably failed consistency check for standard-space registrations!" >> report.siena fi if [ $stdmask = 1 ] ; then ${FSLDIR}/bin/flirt -in ${FSLDIR}/data/standard/MNI152_T1_2mm_brain_mask_dil -ref $A -out ${A}_stdmask -applyxfm -init ${A}_to_std_inv.mat ${FSLDIR}/bin/flirt -in ${FSLDIR}/data/standard/MNI152_T1_2mm_brain_mask_dil -ref $B -out ${B}_stdmask -applyxfm -init ${B}_to_std_inv.mat ${FSLDIR}/bin/fslmaths ${A}_stdmask -thr 0.5 -bin ${A}_stdmask ${FSLDIR}/bin/fslmaths ${B}_stdmask -thr 0.5 -bin ${B}_stdmask ${FSLDIR}/bin/fslmaths ${A}_brain_mask -mas ${A}_stdmask ${A}_brain_mask ${FSLDIR}/bin/fslmaths ${B}_brain_mask -mas ${B}_stdmask ${B}_brain_mask fi if [ "$stdroi" != "" ] ; then ${FSLDIR}/bin/fslmaths ${FSLDIR}/data/standard/MNI152_T1_2mm_brain_mask -mul 0 -add 1 $stdroi ${A}_and_${B}_stdmask ${FSLDIR}/bin/flirt -in ${A}_and_${B}_stdmask -ref $A -out ${A}_stdmask -applyxfm -init ${A}_to_std_inv.mat ${FSLDIR}/bin/flirt -in ${A}_and_${B}_stdmask -ref $B -out ${B}_stdmask -applyxfm -init ${B}_to_std_inv.mat ${FSLDIR}/bin/fslmaths ${A}_stdmask -thr 0.5 -bin ${A}_stdmask ${FSLDIR}/bin/fslmaths ${B}_stdmask -thr 0.5 -bin ${B}_stdmask ${FSLDIR}/bin/fslmaths ${A}_valid_mask_with_$B -mul ${A}_stdmask ${A}_valid_mask_with_$B ${FSLDIR}/bin/fslmaths ${B}_valid_mask_with_$A -mul ${B}_stdmask ${B}_valid_mask_with_$A fi fi ${FSLDIR}/bin/overlay 0 0 -c $A -a ${A}_valid_mask_with_$B 0.9 3 ${A}_brain_mask 0.9 1.1 ${A}_valid_mask_with_${B}_grot ${FSLDIR}/bin/slicer ${A}_valid_mask_with_${B}_grot -a ${A}_valid_mask_with_${B}.png ${FSLDIR}/bin/overlay 0 0 -c $B -a ${B}_valid_mask_with_$A 0.9 3 ${B}_brain_mask 0.9 1.1 ${B}_valid_mask_with_${A}_grot ${FSLDIR}/bin/slicer ${B}_valid_mask_with_${A}_grot -a ${B}_valid_mask_with_${A}.png ${FSLDIR}/bin/imrm ${A}_valid_mask_with_${B}_grot ${B}_valid_mask_with_${A}_grot echo "


Field-of-view and standard space masking
Red shows the common field-of-view of the two timepoint images and the standard-space-based field-of-view masking (if this was run). Blue shows the brain masks, including standard-space-based brain masking (if this was run). Green shows the intersection of the two.

${Ao}

${Bo}
" >> report.html echo "" >> report.siena echo "---------- change analysis ------------------------------------------" >> report.siena corr1=`${FSLDIR}/bin/siena_cal $A $B 1.002 $sdo $sdopts` corr2=`${FSLDIR}/bin/siena_cal $B $A 1.002 $sdo $sdopts` corr=`echo "10 k $corr1 $corr2 + 2.0 / p" | dc -` echo "corr1=$corr1 corr2=$corr2 corr=$corr" >> report.siena echo "" >> report.siena ${FSLDIR}/bin/siena_diff ${B} ${A} -c $corr $sdo $sdopts >> report.siena pbvc_backward=`grep PBVC report.siena | tail -n 1 | awk '{print $2}' | sed 's/-/_/g'` ${FSLDIR}/bin/overlay 1 0 -c ${B}_halfwayto_${A} -a ${B}_halfwayto_${A}_brain_seg 1.1 3 ${B}_halfwayto_${A}_brain_seg_grot ${FSLDIR}/bin/slicer ${B}_halfwayto_${A}_brain_seg_grot -a ${B}_halfwayto_${A}_brain_seg.png ${FSLDIR}/bin/imrm ${B}_halfwayto_${A}_brain_seg_grot echo "" >> report.siena ${FSLDIR}/bin/siena_diff ${A} ${B} -c $corr $sdo $sdopts >> report.siena ${FSLDIR}/bin/overlay 1 0 -c ${A}_halfwayto_${B} -a ${A}_halfwayto_${B}_brain_seg 1.1 3 ${A}_halfwayto_${B}_brain_seg_grot ${FSLDIR}/bin/slicer ${A}_halfwayto_${B}_brain_seg_grot -a ${A}_halfwayto_${B}_brain_seg.png ${FSLDIR}/bin/imrm ${A}_halfwayto_${B}_brain_seg_grot pbvc_forward=`grep PBVC report.siena | tail -n 1 | awk '{print $2}' | sed 's/-/_/g'` echo "


FAST tissue segmentation
These images show the tissue segmentation used to find the brain/non-brain boundary. The exact segmentation of grey matter vs. white matter is not important.

${Ao}

${Bo}
" >> report.html echo "" >> report.siena pbvc_average=`echo "10 k $pbvc_forward $pbvc_backward - 2.0 / p" | dc -` echo "finalPBVC $pbvc_average %" >> report.siena ${FSLDIR}/bin/fslmaths ${A}_to_${B}_flow -mul -1 ${A}_to_${B}_flowneg ${FSLDIR}/bin/overlay 0 0 ${A}_halfwayto_${B} -a ${A}_to_${B}_flow 0.01 1 ${A}_to_${B}_flowneg 0.01 1 ${A}_halfwayto_${B}_render ${FSLDIR}/bin/slicer ${A}_halfwayto_${B}_render -s 1 -x 0.4 gr${A}_halfwayto_${B}a.png -x 0.5 gr${A}_halfwayto_${B}b.png -x 0.6 gr${A}_halfwayto_${B}c.png -y 0.4 gr${A}_halfwayto_${B}d.png -y 0.5 gr${A}_halfwayto_${B}e.png -y 0.6 gr${A}_halfwayto_${B}f.png -z 0.4 gr${A}_halfwayto_${B}g.png -z 0.5 gr${A}_halfwayto_${B}h.png -z 0.6 gr${A}_halfwayto_${B}i.png ${FSLDIR}/bin/pngappend gr${A}_halfwayto_${B}a.png + gr${A}_halfwayto_${B}b.png + gr${A}_halfwayto_${B}c.png + gr${A}_halfwayto_${B}d.png + gr${A}_halfwayto_${B}e.png + gr${A}_halfwayto_${B}f.png + gr${A}_halfwayto_${B}g.png + gr${A}_halfwayto_${B}h.png + gr${A}_halfwayto_${B}i.png ${A}_halfwayto_${B}_render.png /bin/rm gr${A}_halfwayto_${B}?.??? /bin/cp ${FSLDIR}/etc/luts/ramp.gif .ramp.gif /bin/cp ${FSLDIR}/etc/luts/ramp2.gif .ramp2.gif echo "


Final brain edge movement image

atrophy 0 \"growth\"

Estimated PBVC: $pbvc_average" >> report.html if [ $do_viena = yes ] ; then ${FSLDIR}/bin/viena_quant ${A} ${B} ${Vmask} fi if [ $debug = 0 ] ; then $FSLDIR/bin/imrm \ ${A}_brain ${A}_brain_mask ${A}_brain_skull \ ${B}_brain ${B}_brain_mask ${B}_brain_skull \ ${A}_halfwayto_${B} ${A}_halfwayto_${B}_mask \ ${B}_halfwayto_${A} ${B}_halfwayto_${A}_mask \ ${A}_halfwayto_${B}_stdmask \ ${B}_halfwayto_${A}_stdmask \ ${A}_halfwayto_${B}_brain \ ${A}_to_std ${B}_to_std \ ${A}_halfwayto_${B}_brain_seg \ ${A}_to_${B}_flowneg \ ${B}_halfwayto_${A}_brain \ ${B}_halfwayto_${A}_brain_seg \ ${B}_to_${A}_flowneg \ ${A}_stdmask ${B}_stdmask \ ${A}_and_${B}_stdmask \ ${A}_valid_mask ${B}_valid_mask ${A}_valid_mask_to_${B} ${B}_valid_mask_to_${A} \ ${A}_valid_mask_with_$B ${B}_valid_mask_with_$A ${A}_halfwayto_${B}_valid_mask ${B}_halfwayto_${A}_valid_mask /bin/rm -f ${A}_to_${B}_tmp.mat ${A}_halfwayto_${B}_brain.vol ${B}_halfwayto_${A}_brain.vol ${B}_to_${A}.mat_avscale fi if [ $do_viena = yes ] && [ $debug = 0 ] ; then $FSLDIR/bin/imrm \ ${A}_halfwayto_sc${A}_brain \ ${A}_to_sc${A}_flow ${A}_to_sc${A}_flow_ventricles \ ${A}_to_sc${A}_edgepoints ${A}_to_sc${A}_edgepoints_ventricles \ ${A}_halfwayto_sc${A}_brain_ventricle_region_bin \ ${B}_halfwayto_sc${B}_brain \ ${B}_to_sc${B}_flow ${B}_to_sc${B}_flow_ventricles \ ${B}_to_sc${B}_edgepoints ${B}_to_sc${B}_edgepoints_ventricles \ ${B}_halfwayto_sc${B}_brain_ventricle_region_bin \ ${A}_to_${B}_edgepoints ${A}_to_${B}_edgepoints_ventricles \ ${B}_to_${A}_edgepoints ${B}_to_${A}_edgepoints_ventricles \ ${A}_halfwayto_${B}_ventricle_region_bin \ ${B}_halfwayto_${A}_ventricle_region_bin fi echo "" echo "Finished. The SIENA report can be viewed by pointing your web browser at:" echo file:`pwd`/report.html echo "Estimated percentage brain volume change (PBVC) = " echo "$pbvc_average" echo "" cat >> report.html <

SIENA Methods

Two-timepoint percentage brain volume change was estimated with SIENA [Smith 2001, Smith 2002], part of FSL [Smith 2004]. SIENA starts by extracting brain and skull images from the two-timepoint whole-head input data [Smith 2002b]. The two brain images are then aligned to each other [Jenkinson 2001, Jenkinson 2002] (using the skull images to constrain the registration scaling); both brain images are resampled into the space halfway between the two. Next, tissue-type segmentation is carried out [Zhang 2001] in order to find brain/non-brain edge points, and then perpendicular edge displacement (between the two timepoints) is estimated at these edge points. Finally, the mean edge displacement is converted into a (global) estimate of percentage brain volume change between the two timepoints.

[Smith 2001] S.M. Smith, N. De Stefano, M. Jenkinson, and P.M. Matthews.
   Normalised accurate measurement of longitudinal brain change.
   Journal of Computer Assisted Tomography, 25(3):466-475, May/June 2001.

[Smith 2002] S.M. Smith, Y. Zhang, M. Jenkinson, J. Chen, P.M. Matthews, A. Federico, and N. De Stefano.
   Accurate, robust and automated longitudinal and cross-sectional brain change analysis.
   NeuroImage, 17(1):479-489, 2002.

[Smith 2004] S.M. Smith, M. Jenkinson, M.W. Woolrich, C.F. Beckmann, T.E.J. Behrens, H. Johansen-Berg, P.R. Bannister, M. De Luca, I. Drobnjak, D.E. Flitney, R. Niazy, J. Saunders, J. Vickers, Y. Zhang, N. De Stefano, J.M. Brady, and P.M. Matthews.
   Advances in functional and structural MR image analysis and implementation as FSL.
   NeuroImage, 23(S1):208-219, 2004.

[Smith 2002b] S.M. Smith.
   Fast robust automated brain extraction.
   Human Brain Mapping, 17(3):143-155, November 2002.

[Jenkinson 2001] M. Jenkinson and S.M. Smith.
   A global optimisation method for robust affine registration of brain images.
   Medical Image Analysis, 5(2):143-156, June 2001.

[Jenkinson 2002] M. Jenkinson, P.R. Bannister, J.M. Brady, and S.M. Smith.
   Improved optimisation for the robust and accurate linear registration and motion correction of brain images.
   NeuroImage, 17(2):825-841, 2002.

[Zhang 2001] Y. Zhang, M. Brady, and S. Smith.
   Segmentation of brain MR images through a hidden Markov random field model and the expectation maximization algorithm.
   IEEE Trans. on Medical Imaging, 20(1):45-57, 2001. EOF