#!/bin/sh

#   sienax - Structural Image Evaluation, including Normalisation, of Atrophy (X-sectional)
#
#   Stephen Smith, FMRIB Image Analysis Group
#
#   Copyright (C) 1999-2004 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.2beta (c) 2004, 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
}

if [ _$FSLDIR = _ ] ; then
    FSLDIR=/usr/local/fsl
    export FSLDIR
fi

[ _$1 = _ ] && Usage
I=`${FSLDIR}/bin/remove_ext $1`;
[ `${FSLDIR}/bin/imtest $I` = 0 ] && Usage

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.3"        >> ${I}.sienax
echo " sienax $@"                                                              >> ${I}.sienax
echo ""                                                                        >> ${I}.sienax

shift

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 `$FSLDIR/bin/imglob ${I}_brain* ${I}_talmask*`
    /bin/rm -f ${I}2tal.avscale ${I}2tal_inv.mat
fi

echo "$brain"