diff --git a/NOTES b/NOTES
new file mode 100644
index 0000000000000000000000000000000000000000..f3a6ed8f2e4aa3d7d677316553b1f41a0a94a03b
--- /dev/null
+++ b/NOTES
@@ -0,0 +1,3 @@
+
+siena_diff - cannot compile with optimisation turned on (gcc bug)
+
diff --git a/doc/index.html b/doc/index.html
new file mode 100644
index 0000000000000000000000000000000000000000..8319c3e4e826126d8e658953ff5d5d2828d81f7c
--- /dev/null
+++ b/doc/index.html
@@ -0,0 +1,396 @@
+<!-- {{{ start -->
+
+<HTML><HEAD>
+<TITLE>SIENA - Structural Brain Change Analysis</TITLE>
+</HEAD><BODY BACKGROUND="../images/fsl-bg.jpg">
+<hr><TABLE BORDER=0 WIDTH="100%"><TR>
+<TD ALIGN=CENTER><H1>SIENA Structural Brain Change Analysis<br>User
+Guide<br></H1>(SIENA - Structural Image Evaluation, using
+Normalisation, of Atrophy - Version 2.2)
+<TD ALIGN=RIGHT><a href="../index.html"><IMG BORDER=0 SRC="../images/fsl-logo.jpg"></a>
+</TR></TABLE>
+
+<!-- }}} -->
+<!-- {{{ Introduction -->
+
+<hr><H2>Introduction</H2>
+
+<p>SIENA is a package for both single-time-point ("cross-sectional")
+and two-time-point ("longitudinal") analysis of brain change, in
+particular, the estimation of atrophy (volumetric loss of brain
+tissue). SIENA has already been used in many clinical studies.
+
+<p><b>siena</b> estimates percentage brain
+volume change (PBVC) betweem two input images, taken of the same
+subject, at different points in time. It calls a series of image
+analysis programs (supplied with <a href="../index.html">FSL</a>) to
+strip the non-brain tissue from the two images, register the two
+brains (under the constraint that the skulls are used to hold the
+scaling constant during the registration) and analyse the brain change
+between the two time points.
+
+<p><b>sienax</b> estimages total brain tissue
+volume, from a single image, after registration to standard
+(Talairach) space. It calls a series of FSL programs: It first strips
+non-brain tissue, and then uses the brain and skull images to estimate
+the scaling between the subject's image and Talairach space. It then
+runs tissue segmentation to estimate the volume of brain tissue, and
+multiplies this by the estimated scaling factor, to reduce
+head-size-related variability between subjects.
+
+<p><b>Contributors</b>: There have been many contributions of various
+kinds from members of the FMRIB analysis group and collaborators
+mentioned on the <a href="../index.html">FSL page</a>.
+
+<p>For more detail on SIENA and updated journal references, see the <A
+HREF="http://www.fmrib.ox.ac.uk/analysis/research/siena/">SIENA web
+page</A>. If you use SIENA in your research, please quote the journal
+references listed there.
+
+<!-- }}} -->
+<!-- {{{ FSL Tools used -->
+
+<hr><H2>FSL Tools used</H2>
+
+This section lists the generic FSL programs that SIENA uses.
+
+<p><a href="../bet/index.html">bet</a> - Brain Extraction Tool. This
+automatically removes all non-brain tissue from the image. It can
+optionally output the binary brain mask that was derived during this
+process, and output an estimate of the external surface of the skull,
+for use as a scaling constraint in later registration.
+
+<p><b>pairreg</b>, a script supplied with <a
+href="../flirt/index.html">FLIRT</a> - FMRIB's Linear Image
+Registration Tool. This script calls FLIRT with a special optimisation
+schedule, to register two brain images whilst at the same time using
+two skull images to hold the scaling constant (in case the brain has
+shrunk over time, or the scanner calibration has changed). The script
+first calls FLIRT to register the brains as fully as possible. This
+registration is then applied to the skull images, but only the scaling
+and skew are allowed to change. This is then applied to the brain
+images, and a final pass optimally rotates and translates the brains
+to get the best final registration.
+
+<p><a href="../fast/index.html">fast</a> - FMRIB's Automated
+Segmentation Tool. This program automatically segments a brain-only
+image into different tissue types (normally background, grey matter,
+white matter, CSF and other). It also corrects for bias field. It is
+used in various ways in the SIENA scripts. Note that both <b>siena</b>
+and <b>sienax</b> allow you to choose between segmentation of grey
+matter and white matter as separate classes or a single class. It is
+important to choose the right option here, depending on whether there
+is or is not reasonable grey-white contrast in the image.
+
+<!-- }}} -->
+<!-- {{{ Two-Time-Point Estimation -->
+
+<hr><H2>Two-Time-Point Estimation</H2>
+
+<h3>Usage</h3>
+
+<p>The script <a href="siena">siena</a> (see <a
+href="siena_usage">usage</a>) is run simply by typing
+
+<p><b>siena <input1_fileroot> <input2_fileroot></b>
+
+<p>where the two input fileroots are analyze images without the .hdr
+or .img extensions. Note that the input fileroot must not contain
+directory names - i.e. all must be done within a single directory.
+
+<p>Other options are:
+
+<p><b>-d</b> : debug (don't delete intermediate files)
+
+<p><b>-f <BET threshold></b> : Threshold for BET brain
+extraction (default 0.5) - reduce this to make brain estimates larger
+and vice versa
+
+<p><b>-2</b> : two-class segmentation (don't segment grey and white
+matter separately) - use this if there is poor grey/white contrast
+
+<p><b>-m</b> : use Talairach-space masking as well as BET (e.g. if it
+is proving hard to get reliable brain segmentation from BET, for
+example if eyes are hard to segment out) - register to Talairach space
+in order to use a pre-defined Talairach-space brain mask
+
+<p><b>-t <t></b>: ignore from t (mm) upwards in Talairach space
+- if you need to ignore the top part of the head (e.g. if some
+subjects have the top missing and you need consistency across
+subjects)
+
+<p><b>-b <b></b>: ignore from b (mm) downwards in Talairach space; b should probably be -ve
+
+
+<h3>What the script does</h3>
+
+<p><b>siena</b> carries out the following steps:
+
+<p>Run <b>bet</b> on the two input images, producing as output, for
+each input: extracted brain, binary brain mask and skull image. If you
+need to call BET with a different threshold than the default of 0.5,
+use <b>-f <threshold></b>.
+
+<p>Run <b>siena_flirt</b>, a separate script, to register the two
+brain images. This first calls the FLIRT-based registration script
+<b>pairreg</b> (which uses the brain and skull images to carry out
+constrained registration). It then deconstructs the final transform
+into two half-way transforms which take the two brain images into a
+space halfway between the two, so that they both suffer the same
+amount of interpolation-related blurring. Finally the script produces
+a multi-slice gif picture showing the registration quality, with one
+transformed image as the background and edges from the other
+transformed image superimposed in red.
+
+<p>The final step is to carry out change analysis on the registered
+brain images. This is done using the program <b>siena_diff</b>. (In
+order to improve slightly the accuracy of the <b>siena_diff</b>
+program, a self-calibration script <b>siena_cal</b>, described later,
+is run before this.) <b>siena_diff</b> carries out the following
+steps:
+
+<UL>
+
+<LI>Transforms original whole head images and brain masks for each
+time point into the space halfway between them, using the two halfway
+transforms previously generated.
+
+<LI>Combines the two aligned masks using logical OR (if either is 1
+then the output is 1).
+
+<LI>The combined mask is used to mask the two aligned head
+images, resulting in aligned brain images.
+
+<LI>The change between the two aligned brain images is now estimated,
+using the following method (note that options given to the
+<b>siena</b> script are passed on to <b>siena_diff</b>): Apply tissue
+segmentation to the first brain image. At all points which are
+reported as boundaries between brain and non-brain (including internal
+brain-CSF boundaries), compute the distance that the brain surface has
+moved between the two time points. This motion of the brain edge
+(perpendicular to the local edge) is calulated on the basis of
+sub-voxel correlation (matching) of two 1D vectors; these are taken
+from the 3D images, a fixed distance either side of the surface point,
+and perpendicular to it, and are differentiated before correlation,
+allowing some variation in the two original images. Compute mean
+perpendicular surface motion and convert to PBVC.
+
+<LI>To make this conversion between mean perpendicular edge motion and
+PBVC, it is necessary to assume a certain relationship between real
+brain surface area, number of estimated edge points and real brain
+volume. This number can be estimated for general images, but will vary
+according to slice thickness, image sequence type, etc, causing small
+scaling errors in the final PBVC. In order to correct for this,
+self-calibration is applied, in which <b>siena</b> calls
+<b>siena_cal</b>. This script runs <b>siena_diff</b> on one of the
+input images relative to a scaled version of itself, with the scaling
+pre-determined (and therefore known). Thus the final PBVC is known in
+advance and the estimated value can be compared with this to get a
+correction factor for the current image. This is done for both input
+images and the average taken, to give a correction factor to be fed
+into <b>siena_diff</b>.
+
+</UL>
+
+Note that all output is in the same directory as the input, so this
+must be writable by the user. The output files are (assuming the input
+images are called "A" and "B"):
+
+<UL>
+
+<LI>A_to_B.siena the output information from the <b>siena</b> script.
+
+<LI>A_halfwayto_B_render.hdr a colour rendered image of edge motion
+superimposed on the halfway A image.
+
+<LI>B_regto_A.gif a gif image showing the results of the registration,
+using one transformed image as the background and the other as the
+coloured edges foreground.
+
+<LI>A_to_B.mat the transformation taking A to B, using the brain and
+skull images.
+
+<LI>B_to_A.mat the transformation taking B to A, using the brain and
+skull images.
+
+<LI>A_halfwayto_B.mat and B_halfwayto_A.mat the transformations taking
+the images to the halfway positions.
+
+</UL>
+
+<!-- }}} -->
+<!-- {{{ Single-Time-Point Estimation -->
+
+<hr><H2>Single-Time-Point Estimation</H2>
+
+<h3>Usage</h3>
+
+<p>The script <a href="sienax">sienax</a> (see <a
+href="sienax_usage">usage</a>) is run simply by typing
+
+<p><b>sienax <input_fileroot></b>
+
+<p>where the input fileroot is an analyze image without the .hdr or
+.img extension. Note that the input fileroot must not contain
+directory names - i.e. all must be done within a single directory.
+
+
+<p>Other options are:
+
+<p><b>-d</b> : debug (don't delete intermediate files)
+
+<p><b>-f <BET threshold></b> : Threshold for BET brain
+extraction (default 0.5) - reduce this to make brain estimates larger
+and vice versa
+
+<p><b>-2</b> : two-class segmentation (don't segment grey and white
+matter separately) - use this if there is poor grey/white contrast
+
+<p><b>-t <t></b>: ignore from t (mm) upwards in Talairach space
+- if you need to ignore the top part of the head (e.g. if some
+subjects have the top missing and you need consistency across
+subjects)
+
+<p><b>-b <b></b>: ignore from b (mm) downwards in Talairach space; b should probably be -ve
+
+
+
+<h3>What the script does</h3>
+
+<p><b>sienax</b> carries out the following steps:
+
+<UL>
+
+ <LI>Run <b>bet</b> on the single input image, outputting the
+ extracted brain, and the skull image. If you need to call BET with a
+ different threshold than the default of 0.5, use <b>-f
+ <threshold></b>.
+
+ <LI>Run <b>pairreg</b> (which uses the brain and skull images to
+ carry out constrained registration); the MNI152 standard brain is
+ the target (reference), using brain and skull images derived from
+ the MNI152. Thus, as with two-time-point atrophy, the brain is
+ registered (this time to the standard brain), again using the skull
+ as the scaling constraint. Thus brain tissue volume (estimated
+ below) will be relative to a "normalised" skull size. (Ignore the
+ "WARNING: had difficulty finding robust limits in histogram"
+ message; this appears because FLIRT isn't too happy with the unusual
+ histograms of skull images, but is nothing to worry about in this
+ context.) Note that all later steps are in fact carried out on the
+ original (but stripped) input image, not the registered input image;
+ this is so that the original image does not need to be resampled
+ (which introduces blurring). Instead, to make use of the
+ normalisation described above, the brain volume (estimated by the
+ segmentation step described below) is scaled by a scaling factor
+ derived from the normalising transform, before being reported as the
+ final normalised brain volume.
+
+ <LI>A standard brain image mask, (derived from the MNI152 and
+ slightly dilated) is transformed into the original image space (by
+ inverting the normalising transform found above) and applied to the
+ brain image. This helps ensure that the original brain extraction
+ does not include artefacts such as eyeballs.
+
+ <LI>Segmentation is now run on the masked brain using
+ <b>fast</b>. If there is reasonable grey-white contrast, grey
+ matter and white matter volumes are reported separately, as well as
+ total brain volume (this is the default behaviour). Otherwise
+ (i.e. if <b>sienax</b> was called with the <b>-2</b> option), just
+ brain/CSF/background segmentation is carried out, and only brain
+ volume is reported. Before reporting, all volumes are scaled by the
+ normalising scaling factor, as described above, so that all
+ subjects' volumes are reported relative to a normalised skull size.
+
+</UL>
+
+Note that all output is in the same directory as the input, so this
+must be writable by the user. The output files are (assuming the input
+image is called "A"):
+
+<UL>
+
+<LI>A.sienax the output information from the <b>sienax</b> script.
+
+<LI>A_render.hdr a colour rendered image showing the segmentation
+output superimposed on top of the original image.
+
+<LI>A2tal.mat the transformation that takes the input image into
+standard space.
+
+</UL>
+
+<!-- }}} -->
+<!-- {{{ Manual Segmentation Correction -->
+
+<!-- <hr><H2>Manual Segmentation Correction</H2> -->
+
+<!-- <p>If you want to manually edit the segmentation image, to correct for -->
+<!-- any segmentation errors, you should follow the steps below. The -->
+<!-- following assumes the use of MEDx to edit images. The main steps -->
+<!-- described are for correcting <b>sienax</b>, followed by a description -->
+<!-- of how to adapt this procedure for <b>siena</b>. -->
+
+<!-- <UL> -->
+
+<!-- <LI>Load the ...seg.hdr image into MEDx. -->
+
+<!-- <LI>In lightbox mode, make sure the dimensionality is 3D. -->
+
+<!-- <LI>In "look-at-voxel-values" mode rather than the "graphics-edit" -->
+<!-- mode (controlled by the button to the right of the 3D/2D button), -->
+<!-- check the intensities of the classes of interest - for example, -->
+<!-- 0=background, 1=CSF, 2=grey, 3=white. -->
+
+<!-- <LI>Using <b>Toolbox->Measurement->Statistics</b>, set the lower and -->
+<!-- upper thresholds to give the current count of G+W (eg lower=1.5, -->
+<!-- upper=3.5) and measure the <b>Area</b> (in fact the volume). This is -->
+<!-- slightly different (and less accurate) than that reported by fast -->
+<!-- as it does not directly take into account partial volume effect. -->
+
+<!-- <LI>In lightbox mode, change dimensionality from 3D to 2D. -->
+
+<!-- <LI>Using the pencil tool from the toolbox, outline all regions of -->
+<!-- one particular type that you wish to change to a different type. -->
+
+<!-- <LI>Use <b>Folder->Graphic->Select All->2D Graphics</b> to highlight -->
+<!-- all graphics that you've just drawn. Use -->
+<!-- <b>Folder->Graphic->Grouping->Derive 3D Contour</b> to convert these -->
+<!-- into one 3D object. -->
+
+<!-- <LI>Bring up the image calculator and enter the intensity of the -->
+<!-- class that you wish to set all highlighted voxels to (eg 2). Press -->
+<!-- <b>OK</b>. Delete the 3D contour with -->
+<!-- <b>Folder->Graphic->Delete</b>, and reset the display range using -->
+<!-- <b>Folder->Display->Display Range->Compute</b>. -->
+
+<!-- <LI>Now repeat the last 4 steps if you wish to also change other -->
+<!-- areas to a different class type to that which you have just set. -->
+
+<!-- <LI>Now re-run the measurement of "<b>Area</b>" using -->
+<!-- statistics. Subtract the two volumes, multiply by the VSCALING value -->
+<!-- (from the .sienax output file) and amend the BRAIN volume reported -->
+<!-- at the end of this file. -->
+
+<!-- </UL> -->
+
+<!-- If you want to amend the segmentation carried out within <b>siena</b>, -->
+<!-- carry out the above (without bothering with any of the statistics -->
+<!-- steps) and save the image, overwriting the original segmentation -->
+<!-- output. Now if you re-run <b>siena</b> the edited segmentation result -->
+<!-- will be used (as <b>siena</b> doesn't re-run segmentation if a -->
+<!-- segmentation output image is already in existence). -->
+
+<!-- }}} -->
+<!-- {{{ end -->
+
+<p><HR><FONT SIZE=1>Copyright © 2000, University of
+Oxford. Written by <A
+HREF="http://www.fmrib.ox.ac.uk/~steve/index.html">S. Smith</A>.</FONT>
+
+<br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>
+<br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>
+<br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>
+
+</BODY></HTML>
+
+<!-- }}} -->
diff --git a/doc/siena b/doc/siena
new file mode 100755
index 0000000000000000000000000000000000000000..2b1800cc1d11cf3b74518ce2ace5dea7c68c6017
--- /dev/null
+++ b/doc/siena
@@ -0,0 +1,265 @@
+#!/bin/sh
+
+# siena - Structural Image Evaluation, including Normalisation, of Atrophy
+#
+# Stephen Smith, FMRIB Image Analysis Group
+#
+# Copyright (C) 1999-2002 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 3.0 (c) 2002, 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/1112.
+
+Usage() {
+ echo ""
+ echo "Usage: siena <input1_fileroot> <input2_fileroot> [-d] [-f <BET threshold>] [-2] [-t2] [-m] [-t <t>] [-b <b>] [siena_diff options]"
+ echo ""
+ echo "-d : debug (don't delete intermediate files)"
+ echo "-f <f> : Threshold for BET brain extraction (default 0.5)"
+ echo "-2 : two-class segmentation (don't segment grey and white matter separately)"
+ echo "-t2 : T2-weighted input image (default T1-weighted)"
+ echo "-m : use Talairach-space masking as well as BET"
+ echo "-t <t> : ignore from t (mm) upwards in Talairach space"
+ echo "-b <b> : ignore from b (mm) downwards in Talairach space; b should probably be -ve"
+ echo ""
+ exit
+}
+
+[ "$1" = "" ] && Usage
+[ "$2" = "" ] && Usage
+[ ! -f $1.hdr ] && Usage
+[ ! -f $2.hdr ] && Usage
+A=$1
+B=$2
+
+echo "-----------------------------------------------------------------------" > ${A}_to_${B}.siena
+echo "" >> ${A}_to_${B}.siena
+echo " SIENA - Structural Image Evaluation, using Normalisation, of Atrophy" >> ${A}_to_${B}.siena
+echo " part of FSL www.fmrib.ox.ac.uk/fsl" >> ${A}_to_${B}.siena
+echo " running longitudinal atrophy measurement: siena version 2.2" >> ${A}_to_${B}.siena
+echo " siena $@" >> ${A}_to_${B}.siena
+echo "" >> ${A}_to_${B}.siena
+
+shift 2
+
+if [ _$FSLDIR = _ ] ; then
+ FSLDIR=/usr/local/fsl
+ export FSLDIR
+fi
+
+debug=0
+betf=0.5
+sdo="-m"
+dotal=0
+talmask=0
+talroi=""
+origin3=37 # `avwval ${FSLDIR}/etc/standard/avg152T1 origin3`
+pixdim3=2 # `avwval ${FSLDIR}/etc/standard/avg152T1 pixdim3`
+
+for opts in $@ ; do
+
+ if [ $opts = -d ] ; then
+ debug=1
+ ov=-ov
+ shift
+ fi
+
+ if [ $opts = -f ] ; then
+ betf=$2
+ shift 2
+ fi
+
+ if [ $opts = -2 ] ; then
+ sdo="$sdo -2"
+ shift
+ fi
+
+ if [ $opts = -t2 ] ; then
+ is_t2=" -s -t2"
+ shift
+ fi
+
+ if [ $opts = -m ] ; then
+ talmask=1
+ dotal=1
+ shift
+ fi
+
+ if [ $opts = -t ] ; then
+ dotal=1
+ talt=`echo $2 | sed 's/-/_/g'`
+ talt=`echo "10 k $talt $pixdim3 / $origin3 + p" | dc -`
+ talroi="$talroi -roi 0 1000000 0 1000000 0 $talt 0 1"
+ shift 2
+ fi
+
+ if [ $opts = -b ] ; then
+ dotal=1
+ talb=`echo $2 | sed 's/-/_/g'`
+ talb=`echo "10 k $talb $pixdim3 / $origin3 + p" | dc -`
+ talroi="$talroi -roi 0 1000000 0 1000000 $talb 1000000 0 1"
+ shift 2
+ fi
+
+done
+
+sdo="${sdo}${is_t2}"
+
+echo "---------- extract brain --------------------------------------------" >> ${A}_to_${B}.siena
+${FSLDIR}/bin/bet $A ${A}_brain -s -m -f $betf >> ${A}_to_${B}.siena
+${FSLDIR}/bin/bet $B ${B}_brain -s -m -f $betf >> ${A}_to_${B}.siena
+
+echo "" >> ${A}_to_${B}.siena
+echo "---------- register brains and skulls -------------------------------" >> ${A}_to_${B}.siena
+echo "(do not worry about histogram warnings)" >> ${A}_to_${B}.siena
+${FSLDIR}/bin/siena_flirt $A $B >> ${A}_to_${B}.siena 2>&1
+
+echo "" >> ${A}_to_${B}.siena
+echo "---------- produce valid masks --------------------------------------" >> ${A}_to_${B}.siena
+XDIM=`${FSLDIR}/bin/avwval $A dim1` ; XDIM=`echo "$XDIM 2 - p" | dc -`
+YDIM=`${FSLDIR}/bin/avwval $A dim2` ; YDIM=`echo "$YDIM 2 - p" | dc -`
+ZDIM=`${FSLDIR}/bin/avwval $A dim3` ; ZDIM=`echo "$ZDIM 2 - p" | dc -`
+${FSLDIR}/bin/avwmaths ${A}_brain_mask -mul 0 -add 1 -roi 1 $XDIM 1 $YDIM 1 $ZDIM 0 1 ${A}_valid_mask
+XDIM=`${FSLDIR}/bin/avwval $B dim1` ; XDIM=`echo "$XDIM 2 - p" | dc -`
+YDIM=`${FSLDIR}/bin/avwval $B dim2` ; YDIM=`echo "$YDIM 2 - p" | dc -`
+ZDIM=`${FSLDIR}/bin/avwval $B dim3` ; ZDIM=`echo "$ZDIM 2 - p" | dc -`
+${FSLDIR}/bin/avwmaths ${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/avwmaths ${A}_valid_mask -mul ${B}_valid_mask_to_${A} ${A}_valid_mask_with_$B
+${FSLDIR}/bin/avwmaths ${B}_valid_mask -mul ${A}_valid_mask_to_${B} ${B}_valid_mask_with_$A
+
+if [ $dotal = 1 ] ; then
+ echo "" >> ${A}_to_${B}.siena
+ echo "---------- Talairach space masking ----------------------------------" >> ${A}_to_${B}.siena
+ ${FSLDIR}/bin/flirt -ref ${FSLDIR}/etc/standard/avg152T1_brain -in ${A}_brain -omat ${A}_to_tal.mat >> ${A}_to_${B}.siena
+ ${FSLDIR}/bin/flirt -ref ${FSLDIR}/etc/standard/avg152T1_brain -in ${B}_brain -omat ${B}_to_tal.mat >> ${A}_to_${B}.siena
+ ${FSLDIR}/bin/convert_xfm -matonly -inverse -omat ${A}_to_tal_inv.mat ${A}_to_tal.mat
+ ${FSLDIR}/bin/convert_xfm -matonly -inverse -omat ${B}_to_tal_inv.mat ${B}_to_tal.mat
+
+ ${FSLDIR}/bin/convert_xfm -matonly -concat ${B}_to_tal_inv.mat -omat ${A}_to_${B}_tmp.mat ${A}_to_tal.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 Talairaching is $RMSDIFF mm" >> ${A}_to_${B}.siena
+ if [ $RMSDIFF -ge 10 ] ; then
+ echo "Warning! Probably failed consistency check for Talairach registrations!"
+ echo "Warning! Probably failed consistency check for Talairach registrations!" >> ${A}_to_${B}.siena
+ fi
+
+ if [ $talmask = 1 ] ; then
+ ${FSLDIR}/bin/flirt -in ${FSLDIR}/etc/standard/avg152T1_brain_mask_dil2 -ref $A -out ${A}_talmask -applyxfm -init ${A}_to_tal_inv.mat
+ ${FSLDIR}/bin/flirt -in ${FSLDIR}/etc/standard/avg152T1_brain_mask_dil2 -ref $B -out ${B}_talmask -applyxfm -init ${B}_to_tal_inv.mat
+ ${FSLDIR}/bin/avwmaths ${A}_brain_mask -mas ${A}_talmask ${A}_brain_mask
+ ${FSLDIR}/bin/avwmaths ${B}_brain_mask -mas ${B}_talmask ${B}_brain_mask
+ fi
+
+ if [ "$talroi" != "" ] ; then
+ ${FSLDIR}/bin/avwmaths ${FSLDIR}/etc/standard/avg152T1_brain_mask -mul 0 -add 1 $talroi ${A}_and_${B}_talmask
+ ${FSLDIR}/bin/flirt -in ${A}_and_${B}_talmask -ref $A -out ${A}_talmask -applyxfm -init ${A}_to_tal_inv.mat
+ ${FSLDIR}/bin/flirt -in ${A}_and_${B}_talmask -ref $B -out ${B}_talmask -applyxfm -init ${B}_to_tal_inv.mat
+ ${FSLDIR}/bin/avwmaths ${A}_valid_mask_with_$B -mul ${A}_talmask ${A}_valid_mask_with_$B
+ ${FSLDIR}/bin/avwmaths ${B}_valid_mask_with_$A -mul ${B}_talmask ${B}_valid_mask_with_$A
+ fi
+fi
+
+echo "" >> ${A}_to_${B}.siena
+echo "---------- change analysis ------------------------------------------" >> ${A}_to_${B}.siena
+corr1=`${FSLDIR}/bin/siena_cal $A $B 1.002 $sdo $@`
+corr2=`${FSLDIR}/bin/siena_cal $B $A 1.002 $sdo $@`
+corr=`echo "10 k $corr1 $corr2 + 2.0 / p" | dc -`
+echo "corr1=$corr1 corr2=$corr2 corr=$corr" >> ${A}_to_${B}.siena
+
+echo "" >> ${A}_to_${B}.siena
+${FSLDIR}/bin/siena_diff ${B} ${A} -c $corr $sdo $@ >> ${A}_to_${B}.siena
+pbvc_backward=`grep PBVC ${A}_to_${B}.siena | tail -1 | awk '{print $2}' | sed 's/-/_/g'`
+
+echo "" >> ${A}_to_${B}.siena
+${FSLDIR}/bin/siena_diff ${A} ${B} -c $corr $sdo $@ >> ${A}_to_${B}.siena
+pbvc_forward=`grep PBVC ${A}_to_${B}.siena | tail -1 | awk '{print $2}' | sed 's/-/_/g'`
+
+echo "" >> ${A}_to_${B}.siena
+pbvc_average=`echo "10 k $pbvc_forward $pbvc_backward - 2.0 / p" | dc -`
+echo "finalPBVC $pbvc_average %" >> ${A}_to_${B}.siena
+
+${FSLDIR}/bin/avwmaths ${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
+
+if [ $debug = 0 ] ; then
+ /bin/rm -f ${A}_brain.hdr ${A}_brain.img ${A}_brain_mask.hdr ${A}_brain_mask.img ${A}_brain_skull.hdr ${A}_brain_skull.img \
+ ${B}_brain.hdr ${B}_brain.img ${B}_brain_mask.hdr ${B}_brain_mask.img ${B}_brain_skull.hdr ${B}_brain_skull.img \
+ ${A}_halfwayto_${B}.hdr ${A}_halfwayto_${B}.img ${A}_halfwayto_${B}_mask.hdr ${A}_halfwayto_${B}_mask.img \
+ ${B}_halfwayto_${A}.hdr ${B}_halfwayto_${A}.img ${B}_halfwayto_${A}_mask.hdr ${B}_halfwayto_${A}_mask.img \
+ ${A}_halfwayto_${B}_talmask.hdr ${A}_halfwayto_${B}_talmask.img \
+ ${B}_halfwayto_${A}_talmask.hdr ${B}_halfwayto_${A}_talmask.img \
+ ${A}_halfwayto_${B}_brain.hdr ${A}_halfwayto_${B}_brain.img \
+ ${A}_halfwayto_${B}_brain_seg.hdr ${A}_halfwayto_${B}_brain_seg.img \
+ ${A}_to_${B}_flowneg.hdr ${A}_to_${B}_flowneg.img \
+ ${B}_halfwayto_${A}_brain.hdr ${B}_halfwayto_${A}_brain.img \
+ ${B}_halfwayto_${A}_brain_seg.hdr ${B}_halfwayto_${A}_brain_seg.img \
+ ${B}_to_${A}_flowneg.hdr ${B}_to_${A}_flowneg.img \
+ ${B}_to_${A}.mat_avscale \
+ ${A}_halfwayto_${B}_brain.vol ${B}_halfwayto_${A}_brain.vol \
+ ${A}_talmask.hdr ${A}_talmask.img ${B}_talmask.hdr ${B}_talmask.img \
+ ${A}_and_${B}_talmask.hdr ${A}_and_${B}_talmask.img \
+ ${A}_to_tal_inv.mat ${B}_to_tal_inv.mat ${A}_to_${B}_tmp.mat \
+ ${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.???
+fi
+
+echo "$pbvc_average"
diff --git a/doc/siena_usage b/doc/siena_usage
new file mode 100755
index 0000000000000000000000000000000000000000..76d8a1590e54e146e972f854769cd2aa0ea3653c
--- /dev/null
+++ b/doc/siena_usage
@@ -0,0 +1,11 @@
+
+Usage: siena <input1_fileroot> <input2_fileroot> [-d] [-f <BET threshold>] [-2] [-t2] [-m] [-t <t>] [-b <b>] [siena_diff options]
+
+-d : debug (don't delete intermediate files)
+-f <f> : Threshold for BET brain extraction (default 0.5)
+-2 : two-class segmentation (don't segment grey and white matter separately)
+-t2 : T2-weighted input image (default T1-weighted)
+-m : use Talairach-space masking as well as BET
+-t <t> : ignore from t (mm) upwards in Talairach space
+-b <b> : ignore from b (mm) downwards in Talairach space; b should probably be -ve
+
diff --git a/doc/sienax b/doc/sienax
new file mode 100755
index 0000000000000000000000000000000000000000..44bd675e3672de1a70ade58fc7c1bc5f7487654f
--- /dev/null
+++ b/doc/sienax
@@ -0,0 +1,269 @@
+#!/bin/sh
+
+# sienax - Structural Image Evaluation, including Normalisation, of Atrophy (X-sectional)
+#
+# Stephen Smith, FMRIB Image Analysis Group
+#
+# Copyright (C) 1999-2002 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 3.0 (c) 2002, 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/1112.
+
+Usage() {
+ echo ""
+ echo "Usage: sienax <input_fileroot> [-d] [-f <BET threshold>] [-2] [-t2] [-t <t>] [-b <b>] [-r] [-lm <lesion_mask>] [segmentation options]"
+ echo ""
+ echo "-d : debug (don't delete intermediate files)"
+ echo "-f <f> : Threshold for BET brain extraction (default f=0.5)"
+ echo "-2 : two-class segmentation (don't segment grey and white matter separately)"
+ echo "-t2 : T2-weighted input image (default T1-weighted)"
+ echo "-t <t> : ignore from t (mm) upwards in Talairach space"
+ echo "-b <b> : ignore from b (mm) downwards in Talairach space. b should probably be -ve"
+ echo "-r : regional - use standard-space masks to give peripheral cortex GM volume (3-class segmentation only) and ventricular CSF volume"
+ echo "-lm <mask> : use lesion (or lesion+CSF) mask to remove incorrectly labelled \"grey matter\" voxels"
+ echo ""
+ exit
+}
+
+[ _$1 = _ ] && Usage
+[ ! -f $1.hdr ] && Usage
+I=$1
+
+echo "-----------------------------------------------------------------------" > ${I}.sienax
+echo "" >> ${I}.sienax
+echo " SIENA - Structural Image Evaluation, using Normalisation, of Atrophy" >> ${I}.sienax
+echo " part of FSL www.fmrib.ox.ac.uk/fsl" >> ${I}.sienax
+echo " running cross-sectional atrophy measurement: sienax version 2.2" >> ${I}.sienax
+echo " sienax $@" >> ${I}.sienax
+echo "" >> ${I}.sienax
+
+shift
+
+if [ _$FSLDIR = _ ] ; then
+ FSLDIR=/usr/local/fsl
+ export FSLDIR
+fi
+
+debug=0
+regional=0
+betf=0.5
+nseg=3
+talroi=""
+origin3=37 # `avwval ${FSLDIR}/etc/standard/avg152T1 origin3`
+pixdim3=2 # `avwval ${FSLDIR}/etc/standard/avg152T1 pixdim3`
+imtype=-t1
+
+for opts in $@ ; do
+
+ if [ $opts = -d ] ; then
+ debug=1
+ shift
+ fi
+
+ if [ $opts = -r ] ; then
+ regional=1
+ shift
+ fi
+
+ if [ $opts = -f ] ; then
+ betf=$2
+ shift 2
+ fi
+
+ if [ $opts = -2 ] ; then
+ nseg=2
+ shift
+ fi
+
+ if [ $opts = -t2 ] ; then
+ imtype=-t2
+ shift
+ fi
+
+ if [ $opts = -t ] ; then
+ talt=`echo $2 | sed 's/-/_/g'`
+ talt=`echo "10 k $talt $pixdim3 / $origin3 + p" | dc -`
+ talroi="$talroi -roi 0 1000000 0 1000000 0 $talt 0 1"
+ shift 2
+ fi
+
+ if [ $opts = -b ] ; then
+ talb=`echo $2 | sed 's/-/_/g'`
+ talb=`echo "10 k $talb $pixdim3 / $origin3 + p" | dc -`
+ talroi="$talroi -roi 0 1000000 0 1000000 $talb 1000000 0 1"
+ shift 2
+ fi
+
+ if [ $opts = -lm ] ; then
+ lm=$2
+ shift 2
+ fi
+
+done
+
+if [ $regional = 1 ] ; then
+ if [ $nseg = 2 ] ; then
+ echo "Can't do regional analysis with 2-class segmentation"
+ exit
+ fi
+fi
+
+echo "---------- extract brain --------------------------------------------" >> ${I}.sienax
+${FSLDIR}/bin/bet $I ${I}_brain -s -f $betf >> ${I}.sienax
+
+echo "" >> ${I}.sienax
+echo "---------- register to talairach space using brain and skull --------" >> ${I}.sienax
+echo "(do not worry about histogram warnings)" >> ${I}.sienax
+${FSLDIR}/bin/pairreg ${FSLDIR}/etc/standard/avg152T1_brain ${I}_brain ${FSLDIR}/etc/standard/avg152T1_skull ${I}_brain_skull ${I}2tal.mat >> ${I}.sienax 2>&1
+${FSLDIR}/bin/avscale ${I}2tal.mat ${FSLDIR}/etc/standard/avg152T1 > ${I}2tal.avscale
+xscale=`grep Scales ${I}2tal.avscale | awk '{print $4}'`
+yscale=`grep Scales ${I}2tal.avscale | awk '{print $5}'`
+zscale=`grep Scales ${I}2tal.avscale | awk '{print $6}'`
+vscale=`echo "10 k $xscale $yscale * $zscale * p"|dc -`
+echo "VSCALING $vscale" >> ${I}.sienax
+
+echo "" >> ${I}.sienax
+echo "---------- mask with talairach mask ---------------------------------" >> ${I}.sienax
+${FSLDIR}/bin/convert_xfm -matonly -inverse -omat ${I}2tal_inv.mat ${I}2tal.mat
+MASK=${FSLDIR}/etc/standard/avg152T1_brain_mask_dil
+if [ "$talroi" != "" ] ; then
+ ${FSLDIR}/bin/avwmaths $MASK $talroi ${I}_talmaskroi
+ MASK=${I}_talmaskroi
+fi
+${FSLDIR}/bin/flirt -in $MASK -ref ${I}_brain -out ${I}_talmask -applyxfm -init ${I}2tal_inv.mat
+${FSLDIR}/bin/avwmaths ${I}_brain -mask ${I}_talmask ${I}_talmaskbrain
+
+if [ $regional = 1 ] ; then
+ ${FSLDIR}/bin/flirt -in ${FSLDIR}/etc/standard/avg152T1_strucseg_periph -ref ${I}_brain -out ${I}_talmask_segperiph -applyxfm -init ${I}2tal_inv.mat
+ ${FSLDIR}/bin/avwmaths ${FSLDIR}/etc/standard/avg152T1_strucseg -thr 4.5 -bin ${I}_tmpmask
+ ${FSLDIR}/bin/flirt -in ${I}_tmpmask -ref ${I}_brain -out ${I}_talmask_segvent -applyxfm -init ${I}2tal_inv.mat
+ /bin/rm ${I}_tmpmask*
+fi
+
+echo "" >> ${I}.sienax
+echo "---------- segment tissue into types --------------------------------" >> ${I}.sienax
+if [ $nseg = 2 ] ; then
+ ${FSLDIR}/bin/fast -c 2 $imtype -e -ov $@ ${I}_talmaskbrain >> ${I}.sienax 2>&1
+ echo "" >> ${I}.sienax
+ echo "---------- convert brain volume into normalised volume --------------" >> ${I}.sienax
+ S=`${FSLDIR}/bin/avwstats ${I}_talmaskbrain_pve_1 -m -v`
+ xa=`echo $S | awk '{print $1}'`
+ xb=`echo $S | awk '{print $3}'`
+ brain=`echo "2 k $xa $xb * $vscale * p" | dc -`
+else
+ ${FSLDIR}/bin/fast $imtype -e -ov $@ ${I}_talmaskbrain >> ${I}.sienax 2>&1
+
+ if [ _$lm != _ ] ; then
+ ${FSLDIR}/bin/avwmaths_32R $lm -bin -mul -1 -add 1 -mul ${I}_talmaskbrain_pve_1 ${I}_talmaskbrain_pve_1
+ fi
+
+ echo "" >> ${I}.sienax
+ echo "---------- convert brain volume into normalised volume --------------" >> ${I}.sienax
+ if [ $regional = 1 ] ; then
+ ${FSLDIR}/bin/avwmaths_32R ${I}_talmaskbrain_pve_1 -mas ${I}_talmask_segperiph ${I}_talmaskbrain_pve_1_segperiph
+ S=`${FSLDIR}/bin/avwstats ${I}_talmaskbrain_pve_1_segperiph -m -v`
+ xa=`echo $S | awk '{print $1}'`
+ xb=`echo $S | awk '{print $3}'`
+ xg=`echo "2 k $xa $xb * $vscale * p" | dc -`
+ echo "" >> ${I}.sienax
+ echo "pgrey $xg" >> ${I}.sienax
+
+ ${FSLDIR}/bin/avwmaths_32R ${I}_talmaskbrain_pve_0 -mas ${I}_talmask_segvent ${I}_talmaskbrain_pve_0_segvent
+ S=`${FSLDIR}/bin/avwstats ${I}_talmaskbrain_pve_0_segvent -m -v`
+ xa=`echo $S | awk '{print $1}'`
+ xb=`echo $S | awk '{print $3}'`
+ xg=`echo "2 k $xa $xb * $vscale * p" | dc -`
+ echo "" >> ${I}.sienax
+ echo "vcsf $xg" >> ${I}.sienax
+ fi
+ S=`${FSLDIR}/bin/avwstats ${I}_talmaskbrain_pve_1 -m -v`
+ xa=`echo $S | awk '{print $1}'`
+ xb=`echo $S | awk '{print $3}'`
+ grey=`echo "2 k $xa $xb * $vscale * p" | dc -`
+ S=`${FSLDIR}/bin/avwstats ${I}_talmaskbrain_pve_2 -m -v`
+ xa=`echo $S | awk '{print $1}'`
+ xb=`echo $S | awk '{print $3}'`
+ white=`echo "2 k $xa $xb * $vscale * p" | dc -`
+ brain=`echo "2 k $white $grey + p" | dc -`
+ echo "" >> ${I}.sienax
+ echo "GREY $grey" >> ${I}.sienax
+ echo "WHITE $white" >> ${I}.sienax
+fi
+
+echo "" >> ${I}.sienax
+echo "BRAIN $brain" >> ${I}.sienax
+echo "" >> ${I}.sienax
+
+${FSLDIR}/bin/overlay 1 1 -c ${I} -a ${I}_talmaskbrain_seg 1.9 5 ${I}_render
+
+if [ $regional = 1 ] ; then
+ ${FSLDIR}/bin/overlay 1 1 -c ${I} -a ${I}_talmaskbrain_pve_1_segperiph 0.3 0.7 ${I}_periph_render
+ ${FSLDIR}/bin/overlay 1 1 -c ${I} -a ${I}_talmaskbrain_pve_0_segvent 0.3 0.7 ${I}_vent_render
+fi
+
+if [ $debug = 0 ] ; then
+ /bin/rm -f ${I}_brain.hdr ${I}_brain.img ${I}_brain_skull.hdr ${I}_brain_skull.img \
+ ${I}_talmask* ${I}2tal.avscale ${I}2tal_inv.mat
+fi
+
+echo "$brain"
+
diff --git a/doc/sienax_usage b/doc/sienax_usage
new file mode 100755
index 0000000000000000000000000000000000000000..8c2c12c7042d63aa1e94bf9c433c5b28538ca2c9
--- /dev/null
+++ b/doc/sienax_usage
@@ -0,0 +1,12 @@
+
+Usage: sienax <input_fileroot> [-d] [-f <BET threshold>] [-2] [-t2] [-t <t>] [-b <b>] [-r] [-lm <lesion_mask>] [segmentation options]
+
+-d : debug (don't delete intermediate files)
+-f <f> : Threshold for BET brain extraction (default f=0.5)
+-2 : two-class segmentation (don't segment grey and white matter separately)
+-t2 : T2-weighted input image (default T1-weighted)
+-t <t> : ignore from t (mm) upwards in Talairach space
+-b <b> : ignore from b (mm) downwards in Talairach space. b should probably be -ve
+-r : regional - use standard-space masks to give peripheral cortex GM volume (3-class segmentation only) and ventricular CSF volume
+-lm <mask> : use lesion (or lesion+CSF) mask to remove incorrectly labelled "grey matter" voxels
+
diff --git a/lesions_0.1 b/lesions_0.1
new file mode 100755
index 0000000000000000000000000000000000000000..71ca63e636c441c8f65bbcad99fd8351e49637c6
--- /dev/null
+++ b/lesions_0.1
@@ -0,0 +1,90 @@
+#!/bin/sh
+
+#### process options ###################################################################
+
+if [ _$2 = _ ] ; then
+ echo "Usage: lesions <t1_root> <t2_root> [-f]"
+ echo "-f : second image is FLAIR not t2"
+ exit
+fi
+
+nclasses=3
+t2segopts="-t2 -a"
+t2lesions=$2_brain_pve_0
+
+if [ _$3 = _-f ] ; then
+ echo "processing second input as FLAIR image"
+ nclasses=4
+ t2segopts=""
+ t2lesions=$2_brain_pve_3
+fi
+
+
+
+#### T1 ################################################################################
+
+# brain-extract t1 and run segmentation; pve_0 should be CSF only
+${FSLDIR}/bin/bet $1 $1_brain -m
+${FSLDIR}/bin/fast -e -ov $1_brain
+
+# register t1 to standard space and invert transform
+${FSLDIR}/bin/flirt -ref ${FSLDIR}/etc/standard/avg152T1_brain -in $1_brain -omat $1_brain_2_tal.mat
+${FSLDIR}/bin/convert_xfm -inverse -matonly -omat tal_2_$1_brain.mat $1_brain_2_tal.mat
+
+# dilate (6mm) standard space structural segmentation, bring into t1 space and re-binarise
+${FSLDIR}/bin/avwmaths ${FSLDIR}/etc/standard/avg152T1_strucseg -thr 3.5 -bin -dil -dil -dil strucseg_2_$1_brain
+${FSLDIR}/bin/flirt -in strucseg_2_$1_brain -ref $1_brain -out strucseg_2_$1_brain -applyxfm -init tal_2_$1_brain.mat
+${FSLDIR}/bin/avwmaths strucseg_2_$1_brain -bin strucseg_2_$1_brain
+
+# produce t1-derived CSF mask and also dilate in-plane: 5 voxels outside central mask, 1 voxel inside mask
+${FSLDIR}/bin/avwmaths $1_brain_pve_0 -thr 0.3 -bin $1_csf
+${FSLDIR}/bin/avwmaths $1_csf -dil2 -dil2 -dil2 -dil2 -dil2 -sub strucseg_2_$1_brain -thr 0.5 -bin $1_grot
+${FSLDIR}/bin/avwmaths $1_csf -dil2 -mul strucseg_2_$1_brain -add $1_grot -bin $1_csf_dil
+
+
+
+#### T2 ################################################################################
+
+# register t2 to t1 and invert transform
+${FSLDIR}/bin/flirt -in $2 -ref $1 -omat $2_to_$1.mat
+${FSLDIR}/bin/convert_xfm -inverse -matonly -omat $1_to_$2.mat $2_to_$1.mat
+
+# transform t1 brain mask into t2 space, dilate slightly and apply to t2 to get t2 brain
+${FSLDIR}/bin/flirt -in $1_brain_mask -out $2_brain_mask -ref $2 -applyxfm -init $1_to_$2.mat
+#${FSLDIR}/bin/avwmaths $2_brain_mask -dil $2_brain_mask
+${FSLDIR}/bin/avwmaths $2_brain_mask $2_brain_mask
+${FSLDIR}/bin/avwmaths $2 -mas $2_brain_mask $2_brain
+
+# run segmentation on t2; pve_0 should be CSF+lesions; transform into t1 space, remask with t1 brain mask
+${FSLDIR}/bin/fast $t2segopts -c $nclasses -e -ov $2_brain
+${FSLDIR}/bin/flirt -in $t2lesions -ref $1 -out $1_lesion+CSF -applyxfm -init $2_to_$1.mat
+${FSLDIR}/bin/avwmaths $1_lesion+CSF -mas $1_brain_mask $1_lesion+CSF
+
+
+
+#### combine to remove CSF from CSF+lesions ############################################
+
+# combine masks and use to mask out CSF from t2-derived lesion+CSF probability; output lesion volume
+${FSLDIR}/bin/avwmaths_32R $1_csf -mul -1 -add 1 -mul $1_lesion+CSF -thr 0.3 -bin $1_lesions+ucsf
+${FSLDIR}/bin/avwmaths_32R $1_csf_dil -mul -1 -add 1 -mul $1_lesion+CSF $1_lesions
+echo `${FSLDIR}/bin/avwstats $1_lesions -m -v | awk '{print "2 k " $1 " " $3 " * 1000 / p" }' | dc -` > $1_lesions.txt
+echo "$1 `cat $1_lesions.txt`"
+
+# create display output
+${FSLDIR}/bin/flirt -in $2 -ref $1 -out $2_grot -applyxfm -init $2_to_$1.mat
+${FSLDIR}/bin/slicer $2_grot -A 400 $2
+${FSLDIR}/bin/flirt -in $2_brain_seg -ref $1 -out $2_grot2 -applyxfm -init $2_to_$1.mat
+${FSLDIR}/bin/slicer $2_grot2 -i 0 $nclasses -A 400 $2_brain_seg
+${FSLDIR}/bin/overlay 0 1 $1 -a $1_lesions 0.5 1 $1_lesions_render
+${FSLDIR}/bin/slicer $1_lesions_render -A 400 $1_lesions_render
+${FSLDIR}/bin/slicer $1_brain_seg strucseg_2_$1_brain -A 400 $1_brain_seg
+${FSLDIR}/bin/slicer $1 -A 400 $1
+${FSLDIR}/bin/slicer $1_lesions+ucsf -i 0 1 -A 400 $1_lesions+ucsf
+${FSLDIR}/bin/convert -colors 200 +append $2 $2_brain_seg $1_lesions_render $1_brain_seg $1 $1_lesions+ucsf $1_lesions.gif
+
+
+
+#### clean up #########################################################################
+
+/bin/rm -f $1 $2 $1_brain_seg $2_brain_seg $1_lesions_render* $1_brain_delesioned_seg $1_lesions_render* $1_lesions+ucsf $1_grot* $2_grot*
+
diff --git a/lesions_0.2 b/lesions_0.2
new file mode 100755
index 0000000000000000000000000000000000000000..a97e8d3cbcd2d6f7d88742399788cf83718030f6
--- /dev/null
+++ b/lesions_0.2
@@ -0,0 +1,102 @@
+#!/bin/sh
+
+
+
+#### process options ###################################################################
+
+if [ _$2 = _ ] ; then
+ echo "Usage: lesions <t1_root> <t2_root> [-f] [-lt <thr>]"
+ echo "-f : second image is FLAIR not T2"
+ echo "-lt <thr> : set lesion threshold (default 0.5)"
+ exit
+fi
+
+A=$1
+B=$2
+shift 2
+
+lesion_thr=0.5
+nclasses=3
+t2segopts="-t2 -a"
+t2lesions=${B}_brain_pve_0
+
+for opts in $@ ; do
+
+ if [ $opts = -f ] ; then
+ echo "processing second input as FLAIR image"
+ nclasses=4
+ t2segopts=""
+ t2lesions=${B}_brain_pve_3
+ shift
+ fi
+
+ if [ $opts = -lt ] ; then
+ lesion_thr=$2
+ shift 2
+ fi
+
+done
+
+
+
+#### T1 ################################################################################
+
+# brain-extract t1 and run segmentation; pve_0 should be CSF only
+${FSLDIR}/bin/bet $A ${A}_brain -m
+${FSLDIR}/bin/fast -e -ov ${A}_brain
+
+# register t1 to standard space and invert transform
+${FSLDIR}/bin/flirt -ref ${FSLDIR}/etc/standard/avg152T1_brain -in ${A}_brain -omat ${A}_brain_2_tal.mat
+${FSLDIR}/bin/convert_xfm -inverse -matonly -omat tal_2_${A}_brain.mat ${A}_brain_2_tal.mat
+
+# dilate (6mm) standard space ventricles+deep-grey, bring into t1 space and re-binarise to make "CENTRAL MASK"
+${FSLDIR}/bin/avwmaths ${FSLDIR}/etc/standard/avg152T1_strucseg -thr 3.5 -bin -dil -dil -dil strucseg_2_${A}_brain
+${FSLDIR}/bin/flirt -in strucseg_2_${A}_brain -ref ${A}_brain -out strucseg_2_${A}_brain -applyxfm -init tal_2_${A}_brain.mat
+${FSLDIR}/bin/avwmaths strucseg_2_${A}_brain -bin strucseg_2_${A}_brain
+
+# produce t1-derived CSF mask and also dilate in-plane: 5 voxels outside CENTRAL MASK, 1 voxel inside mask
+${FSLDIR}/bin/avwmaths ${A}_brain_pve_0 -thr 0.3 -bin ${A}_csf
+${FSLDIR}/bin/avwmaths ${A}_csf -dil2 -dil2 -dil2 -dil2 -dil2 -sub strucseg_2_${A}_brain -thr 0.5 -bin ${A}_grot
+${FSLDIR}/bin/avwmaths ${A}_csf -dil2 -mul strucseg_2_${A}_brain -add ${A}_grot -bin ${A}_csf_dil
+
+
+
+#### T2 or FLAIR #######################################################################
+
+# register t2 to t1 and invert transform
+${FSLDIR}/bin/flirt -in $B -ref $A -omat ${B}_to_${A}.mat
+${FSLDIR}/bin/convert_xfm -inverse -matonly -omat ${A}_to_${B}.mat ${B}_to_${A}.mat
+
+# transform t1 brain mask into t2 space and apply to t2 to get t2_brain
+${FSLDIR}/bin/flirt -in ${A}_brain_mask -out ${B}_brain_mask -ref $B -applyxfm -init ${A}_to_${B}.mat
+${FSLDIR}/bin/avwmaths $B -mas ${B}_brain_mask ${B}_brain
+
+# run segmentation on t2; $t2lesions should be CSF+lesions; transform into t1 space, remask with t1 brain mask
+${FSLDIR}/bin/fast $t2segopts -c $nclasses -e -ov ${B}_brain
+${FSLDIR}/bin/flirt -in $t2lesions -ref $A -out ${A}_lesion+CSF -applyxfm -init ${B}_to_${A}.mat
+${FSLDIR}/bin/avwmaths ${A}_lesion+CSF -mas ${A}_brain_mask ${A}_lesion+CSF
+
+
+
+#### combine to remove CSF from CSF+lesions ############################################
+
+# combine masks and use to mask out CSF from t2-derived lesion+CSF probability; output lesion volume
+${FSLDIR}/bin/avwmaths_32R ${A}_csf_dil -mul -1 -add 1 -mul ${A}_lesion+CSF ${A}_lesions
+${FSLDIR}/bin/avwmaths_32R ${A}_lesions -thr $lesion_thr ${A}_lesions_thr
+echo `${FSLDIR}/bin/avwstats ${A}_lesions_thr -m -v | awk '{print "2 k " $1 " " $3 " * 1000 / p" }' | dc -` > ${A}_lesions.txt
+echo "$A `cat ${A}_lesions.txt`"
+
+# create display output
+ ${FSLDIR}/bin/slicer $A -A 400 $A
+ ${FSLDIR}/bin/slicer ${A}_brain_seg strucseg_2_${A}_brain -A 400 ${A}_brain_seg
+${FSLDIR}/bin/flirt -in $B -ref $A -out ${B}_grot -applyxfm -init ${B}_to_${A}.mat ; ${FSLDIR}/bin/slicer ${B}_grot -A 400 $B
+${FSLDIR}/bin/flirt -in ${B}_brain_seg -ref $A -out ${B}_grot -applyxfm -init ${B}_to_${A}.mat ; ${FSLDIR}/bin/slicer ${B}_grot -i 0 $nclasses -A 400 ${B}_brain_seg
+${FSLDIR}/bin/overlay 0 1 $A -a ${A}_lesions $lesion_thr 1 ${A}_lesions_render ; ${FSLDIR}/bin/slicer ${A}_lesions_render -A 400 ${A}_lesions_render
+${FSLDIR}/bin/convert -colors 200 +append $A ${A}_brain_seg $B ${B}_brain_seg ${A}_lesions_render ${A}_lesions.gif
+
+
+
+#### clean up #########################################################################
+
+/bin/rm -f $A $B ${A}_brain_seg ${B}_brain_seg ${A}_lesions_render* ${A}_brain_delesioned_seg ${A}_lesions_render* ${A}_grot* ${B}_grot*
+
diff --git a/lesions_0.3 b/lesions_0.3
new file mode 100644
index 0000000000000000000000000000000000000000..c415decfe3617d186ce178d15295de88f97a6aaf
--- /dev/null
+++ b/lesions_0.3
@@ -0,0 +1,102 @@
+#!/bin/sh
+
+
+
+#### process options ###################################################################
+
+if [ _$2 = _ ] ; then
+ echo "Usage: lesions <t1_root> <t2_root> [-f] [-lt <thr>]"
+ echo "-f : second image is FLAIR not T2"
+ echo "-lt <thr> : set lesion threshold (default 0.5)"
+ exit
+fi
+
+A=$1
+B=$2
+shift 2
+
+lesion_thr=0.5
+nclasses=2
+t2segopts="-t2 -a"
+t2lesions=${B}_brain_pve_0
+
+for opts in $@ ; do
+
+ if [ $opts = -f ] ; then
+ echo "processing second input as FLAIR image"
+ nclasses=4
+ t2segopts=""
+ t2lesions=${B}_brain_pve_3
+ shift
+ fi
+
+ if [ $opts = -lt ] ; then
+ lesion_thr=$2
+ shift 2
+ fi
+
+done
+
+echo "working" > ${A}_lesions.txt
+
+#### T1 ################################################################################
+
+# brain-extract t1 and run segmentation; pve_0 should be CSF only
+${FSLDIR}/bin/bet $A ${A}_brain -m
+${FSLDIR}/bin/fast -e -ov ${A}_brain
+
+# register t1 to standard space and invert transform
+${FSLDIR}/bin/flirt -ref ${FSLDIR}/etc/standard/avg152T1_brain -in ${A}_brain -omat ${A}_brain_2_tal.mat
+${FSLDIR}/bin/convert_xfm -inverse -matonly -omat tal_2_${A}_brain.mat ${A}_brain_2_tal.mat
+
+# dilate (6mm) standard space ventricles+deep-grey, bring into t1 space and re-binarise to make "CENTRAL MASK"
+${FSLDIR}/bin/avwmaths ${FSLDIR}/etc/standard/avg152T1_strucseg -thr 3.5 -bin -dil -dil -dil strucseg_2_${A}_brain
+${FSLDIR}/bin/flirt -in strucseg_2_${A}_brain -ref ${A}_brain -out strucseg_2_${A}_brain -applyxfm -init tal_2_${A}_brain.mat
+${FSLDIR}/bin/avwmaths strucseg_2_${A}_brain -bin strucseg_2_${A}_brain
+
+# produce t1-derived CSF mask and also dilate in-plane: 5 voxels outside CENTRAL MASK, 1 voxel inside mask
+${FSLDIR}/bin/avwmaths ${A}_brain_pve_0 -thr 0.3 -bin ${A}_csf
+${FSLDIR}/bin/avwmaths ${A}_csf -dil2 -dil2 -dil2 -dil2 -dil2 -sub strucseg_2_${A}_brain -thr 0.5 -bin ${A}_grot
+${FSLDIR}/bin/avwmaths ${A}_csf -dil2 -mul strucseg_2_${A}_brain -add ${A}_grot -bin ${A}_csf_dil
+
+
+
+#### T2 or FLAIR #######################################################################
+
+# register t2 to t1 and invert transform
+${FSLDIR}/bin/flirt -in $B -ref $A -omat ${B}_to_${A}.mat
+${FSLDIR}/bin/convert_xfm -inverse -matonly -omat ${A}_to_${B}.mat ${B}_to_${A}.mat
+
+# transform t1 brain mask into t2 space and apply to t2 to get t2_brain
+${FSLDIR}/bin/flirt -in ${A}_brain_mask -out ${B}_brain_mask -ref $B -applyxfm -init ${A}_to_${B}.mat
+${FSLDIR}/bin/avwmaths $B -mas ${B}_brain_mask ${B}_brain
+
+# run segmentation on t2; $t2lesions should be CSF+lesions; transform into t1 space, remask with t1 brain mask
+${FSLDIR}/bin/fast $t2segopts -c $nclasses -e -ov ${B}_brain
+${FSLDIR}/bin/flirt -in $t2lesions -ref $A -out ${A}_lesion+CSF -applyxfm -init ${B}_to_${A}.mat
+${FSLDIR}/bin/avwmaths ${A}_lesion+CSF -mas ${A}_brain_mask ${A}_lesion+CSF
+
+
+
+#### combine to remove CSF from CSF+lesions ############################################
+
+# combine masks and use to mask out CSF from t2-derived lesion+CSF probability; output lesion volume
+${FSLDIR}/bin/avwmaths_32R ${A}_csf_dil -mul -1 -add 1 -mul ${A}_lesion+CSF ${A}_lesions
+${FSLDIR}/bin/avwmaths_32R ${A}_lesions -thr $lesion_thr ${A}_lesions_thr
+echo `${FSLDIR}/bin/avwstats ${A}_lesions_thr -m -v | awk '{print "2 k " $1 " " $3 " * 1000 / p" }' | dc -` > ${A}_lesions.txt
+echo "$A `cat ${A}_lesions.txt`"
+
+# create display output
+ ${FSLDIR}/bin/slicer $A -A 400 $A
+ ${FSLDIR}/bin/slicer ${A}_brain_seg strucseg_2_${A}_brain -A 400 ${A}_brain_seg
+${FSLDIR}/bin/flirt -in $B -ref $A -out ${B}_grot -applyxfm -init ${B}_to_${A}.mat ; ${FSLDIR}/bin/slicer ${B}_grot -A 400 $B
+${FSLDIR}/bin/flirt -in ${B}_brain_seg -ref $A -out ${B}_grot -applyxfm -init ${B}_to_${A}.mat ; ${FSLDIR}/bin/slicer ${B}_grot -i 0 $nclasses -A 400 ${B}_brain_seg
+${FSLDIR}/bin/overlay 0 1 $A -a ${A}_lesions $lesion_thr 1 ${A}_lesions_render ; ${FSLDIR}/bin/slicer ${A}_lesions_render -A 400 ${A}_lesions_render
+${FSLDIR}/bin/convert -colors 200 +append $A ${A}_brain_seg $B ${B}_brain_seg ${A}_lesions_render ${A}_lesions.gif
+
+
+
+#### clean up #########################################################################
+
+/bin/rm -f $A $B ${A}_brain_seg ${B}_brain_seg ${A}_lesions_render* ${A}_brain_delesioned_seg ${A}_lesions_render* ${A}_grot* ${B}_grot*
+
diff --git a/makedoc b/makedoc
new file mode 100755
index 0000000000000000000000000000000000000000..1dc41bb04c25dfee3141550fc2227f5cbdc2f424
--- /dev/null
+++ b/makedoc
@@ -0,0 +1,8 @@
+#!/bin/sh
+
+/bin/cp -f siena sienax doc
+
+[ -f ${FSLDIR}/etc/insertcopyright ] && ${FSLDIR}/etc/insertcopyright doc/siena doc/sienax
+
+sh ./siena > doc/siena_usage 2>&1
+sh ./sienax > doc/sienax_usage 2>&1
diff --git a/siena b/siena
new file mode 100755
index 0000000000000000000000000000000000000000..8db468baec60783332aa5262bcd153d5a923d8c5
--- /dev/null
+++ b/siena
@@ -0,0 +1,205 @@
+#!/bin/sh
+
+# siena - Structural Image Evaluation, including Normalisation, of Atrophy
+#
+# Stephen Smith, FMRIB Image Analysis Group
+#
+# Copyright (C) 1999-2002 University of Oxford
+#
+# SHCOPYRIGHT
+
+Usage() {
+ echo ""
+ echo "Usage: siena <input1_fileroot> <input2_fileroot> [-d] [-f <BET threshold>] [-2] [-t2] [-m] [-t <t>] [-b <b>] [siena_diff options]"
+ echo ""
+ echo "-d : debug (don't delete intermediate files)"
+ echo "-f <f> : Threshold for BET brain extraction (default 0.5)"
+ echo "-2 : two-class segmentation (don't segment grey and white matter separately)"
+ echo "-t2 : T2-weighted input image (default T1-weighted)"
+ echo "-m : use Talairach-space masking as well as BET"
+ echo "-t <t> : ignore from t (mm) upwards in Talairach space"
+ echo "-b <b> : ignore from b (mm) downwards in Talairach space; b should probably be -ve"
+ echo ""
+ exit
+}
+
+[ "$1" = "" ] && Usage
+[ "$2" = "" ] && Usage
+[ ! -f $1.hdr ] && Usage
+[ ! -f $2.hdr ] && Usage
+A=$1
+B=$2
+
+echo "-----------------------------------------------------------------------" > ${A}_to_${B}.siena
+echo "" >> ${A}_to_${B}.siena
+echo " SIENA - Structural Image Evaluation, using Normalisation, of Atrophy" >> ${A}_to_${B}.siena
+echo " part of FSL www.fmrib.ox.ac.uk/fsl" >> ${A}_to_${B}.siena
+echo " running longitudinal atrophy measurement: siena version 2.2" >> ${A}_to_${B}.siena
+echo " siena $@" >> ${A}_to_${B}.siena
+echo "" >> ${A}_to_${B}.siena
+
+shift 2
+
+if [ _$FSLDIR = _ ] ; then
+ FSLDIR=/usr/local/fsl
+ export FSLDIR
+fi
+
+debug=0
+betf=0.5
+sdo="-m"
+dotal=0
+talmask=0
+talroi=""
+origin3=37 # `avwval ${FSLDIR}/etc/standard/avg152T1 origin3`
+pixdim3=2 # `avwval ${FSLDIR}/etc/standard/avg152T1 pixdim3`
+
+for opts in $@ ; do
+
+ if [ $opts = -d ] ; then
+ debug=1
+ ov=-ov
+ shift
+ fi
+
+ if [ $opts = -f ] ; then
+ betf=$2
+ shift 2
+ fi
+
+ if [ $opts = -2 ] ; then
+ sdo="$sdo -2"
+ shift
+ fi
+
+ if [ $opts = -t2 ] ; then
+ is_t2=" -s -t2"
+ shift
+ fi
+
+ if [ $opts = -m ] ; then
+ talmask=1
+ dotal=1
+ shift
+ fi
+
+ if [ $opts = -t ] ; then
+ dotal=1
+ talt=`echo $2 | sed 's/-/_/g'`
+ talt=`echo "10 k $talt $pixdim3 / $origin3 + p" | dc -`
+ talroi="$talroi -roi 0 1000000 0 1000000 0 $talt 0 1"
+ shift 2
+ fi
+
+ if [ $opts = -b ] ; then
+ dotal=1
+ talb=`echo $2 | sed 's/-/_/g'`
+ talb=`echo "10 k $talb $pixdim3 / $origin3 + p" | dc -`
+ talroi="$talroi -roi 0 1000000 0 1000000 $talb 1000000 0 1"
+ shift 2
+ fi
+
+done
+
+sdo="${sdo}${is_t2}"
+
+echo "---------- extract brain --------------------------------------------" >> ${A}_to_${B}.siena
+${FSLDIR}/bin/bet $A ${A}_brain -s -m -f $betf >> ${A}_to_${B}.siena
+${FSLDIR}/bin/bet $B ${B}_brain -s -m -f $betf >> ${A}_to_${B}.siena
+
+echo "" >> ${A}_to_${B}.siena
+echo "---------- register brains and skulls -------------------------------" >> ${A}_to_${B}.siena
+echo "(do not worry about histogram warnings)" >> ${A}_to_${B}.siena
+${FSLDIR}/bin/siena_flirt $A $B >> ${A}_to_${B}.siena 2>&1
+
+echo "" >> ${A}_to_${B}.siena
+echo "---------- produce valid masks --------------------------------------" >> ${A}_to_${B}.siena
+XDIM=`${FSLDIR}/bin/avwval $A dim1` ; XDIM=`echo "$XDIM 2 - p" | dc -`
+YDIM=`${FSLDIR}/bin/avwval $A dim2` ; YDIM=`echo "$YDIM 2 - p" | dc -`
+ZDIM=`${FSLDIR}/bin/avwval $A dim3` ; ZDIM=`echo "$ZDIM 2 - p" | dc -`
+${FSLDIR}/bin/avwmaths ${A}_brain_mask -mul 0 -add 1 -roi 1 $XDIM 1 $YDIM 1 $ZDIM 0 1 ${A}_valid_mask
+XDIM=`${FSLDIR}/bin/avwval $B dim1` ; XDIM=`echo "$XDIM 2 - p" | dc -`
+YDIM=`${FSLDIR}/bin/avwval $B dim2` ; YDIM=`echo "$YDIM 2 - p" | dc -`
+ZDIM=`${FSLDIR}/bin/avwval $B dim3` ; ZDIM=`echo "$ZDIM 2 - p" | dc -`
+${FSLDIR}/bin/avwmaths ${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/avwmaths ${A}_valid_mask -mul ${B}_valid_mask_to_${A} ${A}_valid_mask_with_$B
+${FSLDIR}/bin/avwmaths ${B}_valid_mask -mul ${A}_valid_mask_to_${B} ${B}_valid_mask_with_$A
+
+if [ $dotal = 1 ] ; then
+ echo "" >> ${A}_to_${B}.siena
+ echo "---------- Talairach space masking ----------------------------------" >> ${A}_to_${B}.siena
+ ${FSLDIR}/bin/flirt -ref ${FSLDIR}/etc/standard/avg152T1_brain -in ${A}_brain -omat ${A}_to_tal.mat >> ${A}_to_${B}.siena
+ ${FSLDIR}/bin/flirt -ref ${FSLDIR}/etc/standard/avg152T1_brain -in ${B}_brain -omat ${B}_to_tal.mat >> ${A}_to_${B}.siena
+ ${FSLDIR}/bin/convert_xfm -matonly -inverse -omat ${A}_to_tal_inv.mat ${A}_to_tal.mat
+ ${FSLDIR}/bin/convert_xfm -matonly -inverse -omat ${B}_to_tal_inv.mat ${B}_to_tal.mat
+
+ ${FSLDIR}/bin/convert_xfm -matonly -concat ${B}_to_tal_inv.mat -omat ${A}_to_${B}_tmp.mat ${A}_to_tal.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 Talairaching is $RMSDIFF mm" >> ${A}_to_${B}.siena
+ if [ $RMSDIFF -ge 10 ] ; then
+ echo "Warning! Probably failed consistency check for Talairach registrations!"
+ echo "Warning! Probably failed consistency check for Talairach registrations!" >> ${A}_to_${B}.siena
+ fi
+
+ if [ $talmask = 1 ] ; then
+ ${FSLDIR}/bin/flirt -in ${FSLDIR}/etc/standard/avg152T1_brain_mask_dil2 -ref $A -out ${A}_talmask -applyxfm -init ${A}_to_tal_inv.mat
+ ${FSLDIR}/bin/flirt -in ${FSLDIR}/etc/standard/avg152T1_brain_mask_dil2 -ref $B -out ${B}_talmask -applyxfm -init ${B}_to_tal_inv.mat
+ ${FSLDIR}/bin/avwmaths ${A}_brain_mask -mas ${A}_talmask ${A}_brain_mask
+ ${FSLDIR}/bin/avwmaths ${B}_brain_mask -mas ${B}_talmask ${B}_brain_mask
+ fi
+
+ if [ "$talroi" != "" ] ; then
+ ${FSLDIR}/bin/avwmaths ${FSLDIR}/etc/standard/avg152T1_brain_mask -mul 0 -add 1 $talroi ${A}_and_${B}_talmask
+ ${FSLDIR}/bin/flirt -in ${A}_and_${B}_talmask -ref $A -out ${A}_talmask -applyxfm -init ${A}_to_tal_inv.mat
+ ${FSLDIR}/bin/flirt -in ${A}_and_${B}_talmask -ref $B -out ${B}_talmask -applyxfm -init ${B}_to_tal_inv.mat
+ ${FSLDIR}/bin/avwmaths ${A}_valid_mask_with_$B -mul ${A}_talmask ${A}_valid_mask_with_$B
+ ${FSLDIR}/bin/avwmaths ${B}_valid_mask_with_$A -mul ${B}_talmask ${B}_valid_mask_with_$A
+ fi
+fi
+
+echo "" >> ${A}_to_${B}.siena
+echo "---------- change analysis ------------------------------------------" >> ${A}_to_${B}.siena
+corr1=`${FSLDIR}/bin/siena_cal $A $B 1.002 $sdo $@`
+corr2=`${FSLDIR}/bin/siena_cal $B $A 1.002 $sdo $@`
+corr=`echo "10 k $corr1 $corr2 + 2.0 / p" | dc -`
+echo "corr1=$corr1 corr2=$corr2 corr=$corr" >> ${A}_to_${B}.siena
+
+echo "" >> ${A}_to_${B}.siena
+${FSLDIR}/bin/siena_diff ${B} ${A} -c $corr $sdo $@ >> ${A}_to_${B}.siena
+pbvc_backward=`grep PBVC ${A}_to_${B}.siena | tail -1 | awk '{print $2}' | sed 's/-/_/g'`
+
+echo "" >> ${A}_to_${B}.siena
+${FSLDIR}/bin/siena_diff ${A} ${B} -c $corr $sdo $@ >> ${A}_to_${B}.siena
+pbvc_forward=`grep PBVC ${A}_to_${B}.siena | tail -1 | awk '{print $2}' | sed 's/-/_/g'`
+
+echo "" >> ${A}_to_${B}.siena
+pbvc_average=`echo "10 k $pbvc_forward $pbvc_backward - 2.0 / p" | dc -`
+echo "finalPBVC $pbvc_average %" >> ${A}_to_${B}.siena
+
+${FSLDIR}/bin/avwmaths ${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
+
+if [ $debug = 0 ] ; then
+ /bin/rm -f ${A}_brain.hdr ${A}_brain.img ${A}_brain_mask.hdr ${A}_brain_mask.img ${A}_brain_skull.hdr ${A}_brain_skull.img \
+ ${B}_brain.hdr ${B}_brain.img ${B}_brain_mask.hdr ${B}_brain_mask.img ${B}_brain_skull.hdr ${B}_brain_skull.img \
+ ${A}_halfwayto_${B}.hdr ${A}_halfwayto_${B}.img ${A}_halfwayto_${B}_mask.hdr ${A}_halfwayto_${B}_mask.img \
+ ${B}_halfwayto_${A}.hdr ${B}_halfwayto_${A}.img ${B}_halfwayto_${A}_mask.hdr ${B}_halfwayto_${A}_mask.img \
+ ${A}_halfwayto_${B}_talmask.hdr ${A}_halfwayto_${B}_talmask.img \
+ ${B}_halfwayto_${A}_talmask.hdr ${B}_halfwayto_${A}_talmask.img \
+ ${A}_halfwayto_${B}_brain.hdr ${A}_halfwayto_${B}_brain.img \
+ ${A}_halfwayto_${B}_brain_seg.hdr ${A}_halfwayto_${B}_brain_seg.img \
+ ${A}_to_${B}_flowneg.hdr ${A}_to_${B}_flowneg.img \
+ ${B}_halfwayto_${A}_brain.hdr ${B}_halfwayto_${A}_brain.img \
+ ${B}_halfwayto_${A}_brain_seg.hdr ${B}_halfwayto_${A}_brain_seg.img \
+ ${B}_to_${A}_flowneg.hdr ${B}_to_${A}_flowneg.img \
+ ${B}_to_${A}.mat_avscale \
+ ${A}_halfwayto_${B}_brain.vol ${B}_halfwayto_${A}_brain.vol \
+ ${A}_talmask.hdr ${A}_talmask.img ${B}_talmask.hdr ${B}_talmask.img \
+ ${A}_and_${B}_talmask.hdr ${A}_and_${B}_talmask.img \
+ ${A}_to_tal_inv.mat ${B}_to_tal_inv.mat ${A}_to_${B}_tmp.mat \
+ ${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.???
+fi
+
+echo "$pbvc_average"
diff --git a/siena_cal b/siena_cal
new file mode 100755
index 0000000000000000000000000000000000000000..540b5d3797490baf3c7e6d55bf99e561a63571e7
--- /dev/null
+++ b/siena_cal
@@ -0,0 +1,85 @@
+#!/bin/sh
+
+# siena_cal - Self-calibration part of SIENA
+#
+# Stephen Smith, FMRIB Image Analysis Group
+#
+# Copyright (C) 1999-2001 University of Oxford
+#
+# SHCOPYRIGHT
+
+Usage() {
+ echo "Usage: siena_cal <input1_fileroot> <input2_fileroot> <scale> [siena_diff options]"
+ exit 1
+}
+
+[ "$1" = "" ] && Usage
+[ "$2" = "" ] && Usage
+[ ! -f $1.hdr ] && Usage
+
+if [ _$FSLDIR = _ ] ; then
+ FSLDIR=/usr/local/fsl
+ export FSLDIR
+fi
+
+input=$1
+inputB=$2
+m=$3
+M=`echo "5 k 1.0 $m / p" | dc -`
+shift 3
+
+o1=`${FSLDIR}/bin/avwval ${input} dim1`
+o2=`${FSLDIR}/bin/avwval ${input} dim2`
+o3=`${FSLDIR}/bin/avwval ${input} dim3`
+
+v1=`${FSLDIR}/bin/avwval ${input} pixdim1 | sed "s/-//g"`
+v2=`${FSLDIR}/bin/avwval ${input} pixdim2 | sed "s/-//g"`
+v3=`${FSLDIR}/bin/avwval ${input} pixdim3 | sed "s/-//g"`
+
+t1=`echo "5 k $m 1.0 - $o1 * $v1 * _0.5 * p" | dc -`
+t2=`echo "5 k $m 1.0 - $o2 * $v2 * _0.5 * p" | dc -`
+t3=`echo "5 k $m 1.0 - $o3 * $v3 * _0.5 * p" | dc -`
+
+echo "$m 0.0 0.0 $t1" > ${input}_halfwayto_sc${input}.mat
+echo "0.0 $m 0.0 $t2" >> ${input}_halfwayto_sc${input}.mat
+echo "0.0 0.0 $m $t3" >> ${input}_halfwayto_sc${input}.mat
+echo "0.0 0.0 0.0 1.0" >> ${input}_halfwayto_sc${input}.mat
+
+t1=`echo "5 k $M 1.0 - $o1 * $v1 * _0.5 * p" | dc -`
+t2=`echo "5 k $M 1.0 - $o2 * $v2 * _0.5 * p" | dc -`
+t3=`echo "5 k $M 1.0 - $o3 * $v3 * _0.5 * p" | dc -`
+
+echo "$M 0.0 0.0 $t1" > sc${input}_halfwayto_${input}.mat
+echo "0.0 $M 0.0 $t2" >> sc${input}_halfwayto_${input}.mat
+echo "0.0 0.0 $M $t3" >> sc${input}_halfwayto_${input}.mat
+echo "0.0 0.0 0.0 1.0" >> sc${input}_halfwayto_${input}.mat
+
+/bin/cp -f ${input}.hdr sc${input}.hdr
+/bin/cp -f ${input}.img sc${input}.img
+/bin/cp -f ${input}_brain_mask.hdr sc${input}_brain_mask.hdr
+/bin/cp -f ${input}_brain_mask.img sc${input}_brain_mask.img
+/bin/cp -f ${input}_valid_mask_with_${inputB}.hdr ${input}_valid_mask_with_sc${input}.hdr
+/bin/cp -f ${input}_valid_mask_with_${inputB}.img ${input}_valid_mask_with_sc${input}.img
+
+${FSLDIR}/bin/siena_diff ${input} sc${input} $@ > ${input}_to_sc${input}.siena
+
+vm=`echo "10 k $m 6 ^ 1.0 - 100.0 * p" | dc -`
+cal=`grep PBVC ${input}_to_sc${input}.siena | awk '{print $2}' | sed "s/-//g"`
+corr=`echo "10 k $vm $cal / p" | dc -`
+
+#echo "$m $vm $cal $corr"
+echo $corr
+
+/bin/rm -f sc${input}.hdr sc${input}.img sc${input}_brain_mask.hdr sc${input}_brain_mask.img \
+ ${input}_halfwayto_sc${input}.hdr ${input}_halfwayto_sc${input}.img \
+ ${input}_halfwayto_sc${input}_brain.hdr ${input}_halfwayto_sc${input}_brain.img \
+ ${input}_halfwayto_sc${input}_brain_seg.hdr ${input}_halfwayto_sc${input}_brain_seg.img \
+ ${input}_halfwayto_sc${input}_mask.hdr ${input}_halfwayto_sc${input}_mask.img \
+ ${input}_halfwayto_sc${input}_talmask.hdr ${input}_halfwayto_sc${input}_talmask.img \
+ sc${input}_halfwayto_${input}.hdr sc${input}_halfwayto_${input}.img \
+ sc${input}_halfwayto_${input}_mask.hdr sc${input}_halfwayto_${input}_mask.img \
+ ${input}_halfwayto_sc${input}.mat sc${input}_halfwayto_${input}.mat \
+ ${input}_to_sc${input}_flow.hdr ${input}_to_sc${input}_flow.img \
+ ${input}_to_sc${input}.siena ${input}_halfwayto_sc${input}_brain.vol \
+ ${input}_valid_mask_with_sc${input}.hdr ${input}_valid_mask_with_sc${input}.img \
+ ${input}_halfwayto_sc${input}_valid_mask.hdr ${input}_halfwayto_sc${input}_valid_mask.img
diff --git a/siena_flirt b/siena_flirt
new file mode 100755
index 0000000000000000000000000000000000000000..6719eea39b66e3e4811ef28533d3e1f849d143b5
--- /dev/null
+++ b/siena_flirt
@@ -0,0 +1,50 @@
+#!/bin/sh
+
+# siena_flirt - Registration part of SIENA
+#
+# Stephen Smith, FMRIB Image Analysis Group
+#
+# Copyright (C) 1999-2001 University of Oxford
+#
+# SHCOPYRIGHT
+
+Usage() {
+ echo "Usage: siena_flirt <input1_fileroot> <input2_fileroot>"
+ exit 1
+}
+
+[ "$1" = "" ] && Usage
+[ "$2" = "" ] && Usage
+[ ! -f $1.hdr ] && Usage
+[ ! -f $2.hdr ] && Usage
+
+if [ _$FSLDIR = _ ] ; then
+ FSLDIR=/usr/local/fsl
+ export FSLDIR
+fi
+
+# register using brains and skulls
+${FSLDIR}/bin/pairreg $1_brain $2_brain $1_brain_skull $2_brain_skull $2_to_$1.mat
+${FSLDIR}/bin/pairreg $2_brain $1_brain $2_brain_skull $1_brain_skull $1_to_$2.mat
+
+# replace both transforms with "average" (reduces error level AND makes system symmetric)
+F=$1_to_$2.mat
+B=$2_to_$1.mat
+${FSLDIR}/bin/convert_xfm -matonly -concat $B -omat tmp_${F}_then_${B} $F
+${FSLDIR}/bin/avscale tmp_${F}_then_${B} $1 > tmp_${F}_then_${B}.avscale
+${FSLDIR}/bin/extracttxt Backward tmp_${F}_then_${B}.avscale 4 1 > tmp_${F}_then_${B}_halfback
+${FSLDIR}/bin/convert_xfm -matonly -concat tmp_${F}_then_${B}_halfback -omat $F $F
+${FSLDIR}/bin/convert_xfm -inverse -matonly -omat $B $F
+/bin/rm tmp_${F}_then_${B} tmp_${F}_then_${B}.avscale tmp_${F}_then_${B}_halfback
+
+# replace the .mat matrix that takes 2->1 with one that takes 2->halfway and one that takes 1->halfway
+${FSLDIR}/bin/avscale $2_to_$1.mat $1_brain > $2_to_$1.mat_avscale
+${FSLDIR}/bin/extracttxt Forward $2_to_$1.mat_avscale 4 1 > $2_halfwayto_$1.mat
+${FSLDIR}/bin/extracttxt Backward $2_to_$1.mat_avscale 4 1 > $1_halfwayto_$2.mat
+
+# produce picture for registration evaluation
+${FSLDIR}/bin/flirt -out $1_halfwayto_$2 -applyxfm -init $1_halfwayto_$2.mat -ref $1 -in $1_brain
+${FSLDIR}/bin/flirt -out $2_halfwayto_$1 -applyxfm -init $2_halfwayto_$1.mat -ref $1 -in $2_brain
+${FSLDIR}/bin/slicer $2_halfwayto_$1 $1_halfwayto_$2 -s 1 -x 0.4 gr$1$2a -x 0.5 gr$1$2b -x 0.6 gr$1$2c -y 0.4 gr$1$2d -y 0.5 gr$1$2e -y 0.6 gr$1$2f -z 0.4 gr$1$2g -z 0.5 gr$1$2h -z 0.6 gr$1$2i
+${FSLDIR}/bin/convert -colors 100 +append gr$1$2a gr$1$2b gr$1$2c gr$1$2d gr$1$2e gr$1$2f gr$1$2g gr$1$2h gr$1$2i $2_regto_$1.gif
+/bin/rm gr$1$2a gr$1$2b gr$1$2c gr$1$2d gr$1$2e gr$1$2f gr$1$2g gr$1$2h gr$1$2i $2_halfwayto_$1.img $2_halfwayto_$1.hdr $1_halfwayto_$2.img $1_halfwayto_$2.hdr
diff --git a/siena_flow2tal b/siena_flow2tal
new file mode 100755
index 0000000000000000000000000000000000000000..81a8d4ead84edfd0b28f81fb1857597362054368
--- /dev/null
+++ b/siena_flow2tal
@@ -0,0 +1,61 @@
+#!/bin/sh
+
+# SIENAL - Local analysis of SIENA flow field
+#
+# Stephen Smith, FMRIB Image Analysis Group
+#
+# Copyright (C) 2002-2003 University of Oxford
+#
+# SHCOPYRIGHT
+
+Usage() {
+ echo "Usage: sienal <fileroot2> <fileroot2> [-d]"
+ echo "-d : debug (don't delete intermediate files)"
+ exit
+}
+
+[ "$1" = "" ] && Usage
+[ "$2" = "" ] && Usage
+[ ! -f $1.hdr ] && Usage
+[ ! -f $2.hdr ] && Usage
+A=$1
+B=$2
+shift 2
+
+if [ _$FSLDIR = _ ] ; then
+ FSLDIR=/usr/local/fsl
+ export FSLDIR
+fi
+
+debug=0
+
+for opts in $@ ; do
+
+ if [ $opts = -d ] ; then
+ debug=1
+ fi
+
+done
+
+# register to MNI152 if not already done, and create halfway2tal transform
+if [ ! -f ${A}_to_tal.mat ] ; then
+ if [ ! -f ${A}_brain.hdr ] ; then
+ ${FSLDIR}/bin/bet $A ${A}_brain
+ fi
+ ${FSLDIR}/bin/flirt -ref ${FSLDIR}/etc/standard/avg152T1_brain -in ${A}_brain -omat ${A}_to_tal.mat
+fi
+${FSLDIR}/bin/convert_xfm -omat ${A}_halfwayto_${B}_inv.mat -inverse ${A}_halfwayto_${B}.mat
+${FSLDIR}/bin/convert_xfm -omat ${A}_halfwayto_${B}_to_tal.mat -concat ${A}_to_tal.mat ${A}_halfwayto_${B}_inv.mat
+
+# dilate flow image, -> standard space -> mask with ss edge image -> blur -> remask
+${FSLDIR}/bin/avwmaths ${A}_to_${B}_flow -dil -dil -dil -dil -dil -dil -dil -dil -dil -dil -dil -dil -dil -dil ${A}_to_${B}_flow_dil
+${FSLDIR}/bin/flirt -in ${A}_to_${B}_flow_dil -ref ${FSLDIR}/etc/standard/avg152T1_brain -out ${A}_to_${B}_flow_to_tal -applyxfm -init ${A}_halfwayto_${B}_to_tal.mat
+${FSLDIR}/bin/avwmaths ${A}_to_${B}_flow_to_tal -mas ${FSLDIR}/etc/standard/avg152T1_edges ${A}_to_${B}_flow_to_tal
+${FSLDIR}/bin/ip ${A}_to_${B}_flow_to_tal ${A}_to_${B}_flow_to_tal 0 -s 10
+${FSLDIR}/bin/avwmaths ${A}_to_${B}_flow_to_tal -mas ${FSLDIR}/etc/standard/avg152T1_edges ${A}_to_${B}_flow_to_tal
+
+# cleanup
+if [ $debug = 0 ] ; then
+ /bin/rm -f ${A}_halfwayto_${B}_inv.mat ${A}_brain.hdr ${A}_brain.img ${A}_to_${B}_flow_dil.hdr ${A}_to_${B}_flow_dil.img
+fi
+
diff --git a/sienax b/sienax
new file mode 100755
index 0000000000000000000000000000000000000000..70440ebb33962ee6484c0a8ef5ae648b83a955e0
--- /dev/null
+++ b/sienax
@@ -0,0 +1,209 @@
+#!/bin/sh
+
+# sienax - Structural Image Evaluation, including Normalisation, of Atrophy (X-sectional)
+#
+# Stephen Smith, FMRIB Image Analysis Group
+#
+# Copyright (C) 1999-2002 University of Oxford
+#
+# SHCOPYRIGHT
+
+Usage() {
+ echo ""
+ echo "Usage: sienax <input_fileroot> [-d] [-f <BET threshold>] [-2] [-t2] [-t <t>] [-b <b>] [-r] [-lm <lesion_mask>] [segmentation options]"
+ echo ""
+ echo "-d : debug (don't delete intermediate files)"
+ echo "-f <f> : Threshold for BET brain extraction (default f=0.5)"
+ echo "-2 : two-class segmentation (don't segment grey and white matter separately)"
+ echo "-t2 : T2-weighted input image (default T1-weighted)"
+ echo "-t <t> : ignore from t (mm) upwards in Talairach space"
+ echo "-b <b> : ignore from b (mm) downwards in Talairach space. b should probably be -ve"
+ echo "-r : regional - use standard-space masks to give peripheral cortex GM volume (3-class segmentation only) and ventricular CSF volume"
+ echo "-lm <mask> : use lesion (or lesion+CSF) mask to remove incorrectly labelled \"grey matter\" voxels"
+ echo ""
+ exit
+}
+
+[ _$1 = _ ] && Usage
+[ ! -f $1.hdr ] && Usage
+I=$1
+
+echo "-----------------------------------------------------------------------" > ${I}.sienax
+echo "" >> ${I}.sienax
+echo " SIENA - Structural Image Evaluation, using Normalisation, of Atrophy" >> ${I}.sienax
+echo " part of FSL www.fmrib.ox.ac.uk/fsl" >> ${I}.sienax
+echo " running cross-sectional atrophy measurement: sienax version 2.2" >> ${I}.sienax
+echo " sienax $@" >> ${I}.sienax
+echo "" >> ${I}.sienax
+
+shift
+
+if [ _$FSLDIR = _ ] ; then
+ FSLDIR=/usr/local/fsl
+ export FSLDIR
+fi
+
+debug=0
+regional=0
+betf=0.5
+nseg=3
+talroi=""
+origin3=37 # `avwval ${FSLDIR}/etc/standard/avg152T1 origin3`
+pixdim3=2 # `avwval ${FSLDIR}/etc/standard/avg152T1 pixdim3`
+imtype=-t1
+
+for opts in $@ ; do
+
+ if [ $opts = -d ] ; then
+ debug=1
+ shift
+ fi
+
+ if [ $opts = -r ] ; then
+ regional=1
+ shift
+ fi
+
+ if [ $opts = -f ] ; then
+ betf=$2
+ shift 2
+ fi
+
+ if [ $opts = -2 ] ; then
+ nseg=2
+ shift
+ fi
+
+ if [ $opts = -t2 ] ; then
+ imtype=-t2
+ shift
+ fi
+
+ if [ $opts = -t ] ; then
+ talt=`echo $2 | sed 's/-/_/g'`
+ talt=`echo "10 k $talt $pixdim3 / $origin3 + p" | dc -`
+ talroi="$talroi -roi 0 1000000 0 1000000 0 $talt 0 1"
+ shift 2
+ fi
+
+ if [ $opts = -b ] ; then
+ talb=`echo $2 | sed 's/-/_/g'`
+ talb=`echo "10 k $talb $pixdim3 / $origin3 + p" | dc -`
+ talroi="$talroi -roi 0 1000000 0 1000000 $talb 1000000 0 1"
+ shift 2
+ fi
+
+ if [ $opts = -lm ] ; then
+ lm=$2
+ shift 2
+ fi
+
+done
+
+if [ $regional = 1 ] ; then
+ if [ $nseg = 2 ] ; then
+ echo "Can't do regional analysis with 2-class segmentation"
+ exit
+ fi
+fi
+
+echo "---------- extract brain --------------------------------------------" >> ${I}.sienax
+${FSLDIR}/bin/bet $I ${I}_brain -s -f $betf >> ${I}.sienax
+
+echo "" >> ${I}.sienax
+echo "---------- register to talairach space using brain and skull --------" >> ${I}.sienax
+echo "(do not worry about histogram warnings)" >> ${I}.sienax
+${FSLDIR}/bin/pairreg ${FSLDIR}/etc/standard/avg152T1_brain ${I}_brain ${FSLDIR}/etc/standard/avg152T1_skull ${I}_brain_skull ${I}2tal.mat >> ${I}.sienax 2>&1
+${FSLDIR}/bin/avscale ${I}2tal.mat ${FSLDIR}/etc/standard/avg152T1 > ${I}2tal.avscale
+xscale=`grep Scales ${I}2tal.avscale | awk '{print $4}'`
+yscale=`grep Scales ${I}2tal.avscale | awk '{print $5}'`
+zscale=`grep Scales ${I}2tal.avscale | awk '{print $6}'`
+vscale=`echo "10 k $xscale $yscale * $zscale * p"|dc -`
+echo "VSCALING $vscale" >> ${I}.sienax
+
+echo "" >> ${I}.sienax
+echo "---------- mask with talairach mask ---------------------------------" >> ${I}.sienax
+${FSLDIR}/bin/convert_xfm -matonly -inverse -omat ${I}2tal_inv.mat ${I}2tal.mat
+MASK=${FSLDIR}/etc/standard/avg152T1_brain_mask_dil
+if [ "$talroi" != "" ] ; then
+ ${FSLDIR}/bin/avwmaths $MASK $talroi ${I}_talmaskroi
+ MASK=${I}_talmaskroi
+fi
+${FSLDIR}/bin/flirt -in $MASK -ref ${I}_brain -out ${I}_talmask -applyxfm -init ${I}2tal_inv.mat
+${FSLDIR}/bin/avwmaths ${I}_brain -mask ${I}_talmask ${I}_talmaskbrain
+
+if [ $regional = 1 ] ; then
+ ${FSLDIR}/bin/flirt -in ${FSLDIR}/etc/standard/avg152T1_strucseg_periph -ref ${I}_brain -out ${I}_talmask_segperiph -applyxfm -init ${I}2tal_inv.mat
+ ${FSLDIR}/bin/avwmaths ${FSLDIR}/etc/standard/avg152T1_strucseg -thr 4.5 -bin ${I}_tmpmask
+ ${FSLDIR}/bin/flirt -in ${I}_tmpmask -ref ${I}_brain -out ${I}_talmask_segvent -applyxfm -init ${I}2tal_inv.mat
+ /bin/rm ${I}_tmpmask*
+fi
+
+echo "" >> ${I}.sienax
+echo "---------- segment tissue into types --------------------------------" >> ${I}.sienax
+if [ $nseg = 2 ] ; then
+ ${FSLDIR}/bin/fast -c 2 $imtype -e -ov $@ ${I}_talmaskbrain >> ${I}.sienax 2>&1
+ echo "" >> ${I}.sienax
+ echo "---------- convert brain volume into normalised volume --------------" >> ${I}.sienax
+ S=`${FSLDIR}/bin/avwstats ${I}_talmaskbrain_pve_1 -m -v`
+ xa=`echo $S | awk '{print $1}'`
+ xb=`echo $S | awk '{print $3}'`
+ brain=`echo "2 k $xa $xb * $vscale * p" | dc -`
+else
+ ${FSLDIR}/bin/fast $imtype -e -ov $@ ${I}_talmaskbrain >> ${I}.sienax 2>&1
+
+ if [ _$lm != _ ] ; then
+ ${FSLDIR}/bin/avwmaths_32R $lm -bin -mul -1 -add 1 -mul ${I}_talmaskbrain_pve_1 ${I}_talmaskbrain_pve_1
+ fi
+
+ echo "" >> ${I}.sienax
+ echo "---------- convert brain volume into normalised volume --------------" >> ${I}.sienax
+ if [ $regional = 1 ] ; then
+ ${FSLDIR}/bin/avwmaths_32R ${I}_talmaskbrain_pve_1 -mas ${I}_talmask_segperiph ${I}_talmaskbrain_pve_1_segperiph
+ S=`${FSLDIR}/bin/avwstats ${I}_talmaskbrain_pve_1_segperiph -m -v`
+ xa=`echo $S | awk '{print $1}'`
+ xb=`echo $S | awk '{print $3}'`
+ xg=`echo "2 k $xa $xb * $vscale * p" | dc -`
+ echo "" >> ${I}.sienax
+ echo "pgrey $xg" >> ${I}.sienax
+
+ ${FSLDIR}/bin/avwmaths_32R ${I}_talmaskbrain_pve_0 -mas ${I}_talmask_segvent ${I}_talmaskbrain_pve_0_segvent
+ S=`${FSLDIR}/bin/avwstats ${I}_talmaskbrain_pve_0_segvent -m -v`
+ xa=`echo $S | awk '{print $1}'`
+ xb=`echo $S | awk '{print $3}'`
+ xg=`echo "2 k $xa $xb * $vscale * p" | dc -`
+ echo "" >> ${I}.sienax
+ echo "vcsf $xg" >> ${I}.sienax
+ fi
+ S=`${FSLDIR}/bin/avwstats ${I}_talmaskbrain_pve_1 -m -v`
+ xa=`echo $S | awk '{print $1}'`
+ xb=`echo $S | awk '{print $3}'`
+ grey=`echo "2 k $xa $xb * $vscale * p" | dc -`
+ S=`${FSLDIR}/bin/avwstats ${I}_talmaskbrain_pve_2 -m -v`
+ xa=`echo $S | awk '{print $1}'`
+ xb=`echo $S | awk '{print $3}'`
+ white=`echo "2 k $xa $xb * $vscale * p" | dc -`
+ brain=`echo "2 k $white $grey + p" | dc -`
+ echo "" >> ${I}.sienax
+ echo "GREY $grey" >> ${I}.sienax
+ echo "WHITE $white" >> ${I}.sienax
+fi
+
+echo "" >> ${I}.sienax
+echo "BRAIN $brain" >> ${I}.sienax
+echo "" >> ${I}.sienax
+
+${FSLDIR}/bin/overlay 1 1 -c ${I} -a ${I}_talmaskbrain_seg 1.9 5 ${I}_render
+
+if [ $regional = 1 ] ; then
+ ${FSLDIR}/bin/overlay 1 1 -c ${I} -a ${I}_talmaskbrain_pve_1_segperiph 0.3 0.7 ${I}_periph_render
+ ${FSLDIR}/bin/overlay 1 1 -c ${I} -a ${I}_talmaskbrain_pve_0_segvent 0.3 0.7 ${I}_vent_render
+fi
+
+if [ $debug = 0 ] ; then
+ /bin/rm -f ${I}_brain.hdr ${I}_brain.img ${I}_brain_skull.hdr ${I}_brain_skull.img \
+ ${I}_talmask* ${I}2tal.avscale ${I}2tal_inv.mat
+fi
+
+echo "$brain"
+