#!/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.

Usage() {
    cat <<EOF

Usage: siena <input1> <input2> [options]

  -o <output-dir>       : set output directory (default output is <input1>_to_<input2>_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 <t>                : ignore from t (mm) upwards in MNI152/Talairach space
  -b <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 <mask image>       : optional user-supplied ventricle mask (default is $FSLDIR/bin/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 '<HTML><HEAD><link REL="stylesheet" TYPE="text/css" href="file:'${FSLDIR}'/doc/fsl.css"><TITLE>FSL</TITLE></HEAD><BODY><hr><TABLE BORDER=0 WIDTH="100%"><TR><TD ALIGN=CENTER><H1>SIENA Report</H1>'${thecommand}'<TD ALIGN=RIGHT><a href="'${FSLDIR}'/doc/index.html"><IMG BORDER=0 SRC="'${FSLDIR}'/doc//images/fsl-logo.jpg"></a></TR></TABLE>' > 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 "<hr><p><b>BET brain extraction results</b><p>${Ao}<br><IMG BORDER=0 SRC=\"${A}_bet.png\"><p>${Bo}<br><IMG BORDER=0 SRC=\"${B}_bet.png\">" >> 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 "<hr><p><b>FLIRT A-to-B registration results</b><p><IMG BORDER=0 SRC=\"${A}_and_${B}.gif\">" >> 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 "<hr><p><b>FLIRT standard space registration results</b><p>${Ao}<br><IMG BORDER=0 SRC=\"${A}_to_std.png\"><p>${Bo}<br><IMG BORDER=0 SRC=\"${B}_to_std.png\">" >> 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 "<hr><p><b>Field-of-view and standard space masking</b><br>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.<p>${Ao}<br><IMG BORDER=0 SRC=\"${A}_valid_mask_with_${B}.png\"><p>${Bo}<br><IMG BORDER=0 SRC=\"${B}_valid_mask_with_${A}.png\">" >> 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 "<hr><p><b>FAST tissue segmentation</b><br>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.<p>${Ao}<br><IMG BORDER=0 SRC=\"${A}_halfwayto_${B}_brain_seg.png\"><p>${Bo}<br><IMG BORDER=0 SRC=\"${B}_halfwayto_${A}_brain_seg.png\">" >> 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 "<hr><p><b>Final brain edge movement image</b><p>atrophy <IMG SRC=\".ramp2.gif\"> 0 <IMG SRC=\".ramp.gif\"> \"growth\"<p><IMG BORDER=0 SRC=\"${A}_halfwayto_${B}_render.png\"><p><b>Estimated PBVC: $pbvc_average</b>" >> report.html

if [ $do_viena = yes ] ; then
    ${FSLDIR}/bin/viena_quant ${A} ${B} ${vienadir}
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 <<EOF

<hr><p><b>SIENA Methods</b>

<p>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.

<font size=-1><em>
<P>[Smith 2001] S.M. Smith, N.&nbsp;De&nbsp;Stefano, M.&nbsp;Jenkinson, and P.M. Matthews.
<BR>&nbsp;&nbsp;&nbsp;Normalised accurate measurement of longitudinal brain change.
<BR>&nbsp;&nbsp;&nbsp;Journal of Computer Assisted Tomography, 25(3):466-475, May/June 2001.

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

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

<P>[Smith 2002b] S.M. Smith.
<BR>&nbsp;&nbsp;&nbsp;Fast robust automated brain extraction.
<BR>&nbsp;&nbsp;&nbsp;Human Brain Mapping, 17(3):143-155, November 2002.

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

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

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

EOF