sienax

#!/bin/sh

#   sienax - Structural Image Evaluation, including Normalisation, of Atrophy (X-sectional)
#
#   Stephen Smith, FMRIB Image Analysis Group
#
#   Copyright (C) 1999-2006 University of Oxford
#
#   SHCOPYRIGHT

Usage() {
    echo ""
    echo "Usage: sienax <input_fileroot> [-d] [-B \"bet options\"] [-2] [-t2] [-t <t>] [-b <b>] [-r] [-lm <lesion_mask>] [-S \"segmentation options\"]"
    echo ""
    echo "-d              : debug (don't delete intermediate files)"
    echo "-B \"betopts\"    : options to pass to BET brain extraction (inside double-quotes), e.g. -B \"-f 0.3\""
    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 MNI152/Talairach space"
    echo "-b <b>          : ignore from b (mm) downwards in MNI152/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 "-S \"segopts\"    : options to pass to FAST segmentation (inside double-quotes), e.g. -S \"-i 20\""
    echo ""
    exit 1
}

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

shift

debug=0
regional=0
betopts=""
segopts=""
nseg=3
stdroi=""
origin3=37 # `fslval ${FSLDIR}/etc/standard/avg152T1 origin3`
pixdim3=2  # `fslval ${FSLDIR}/etc/standard/avg152T1 pixdim3`
imtype=-t1

while [ _$1 != _ ] ; do

    if [ $1 = -d ] ; then
        debug=1
        shift
    elif [ $1 = -r ] ; then
	regional=1
        shift
    elif [ $1 = -B ] ; then
        betopts=$2
        shift 2
    elif [ $1 = -S ] ; then
        segopts=$2
        shift 2
    elif [ $1 = -2 ] ; then
        nseg=2
        shift
    elif [ $1 = -t2 ] ; then
        imtype=-t2
        shift
    elif [ $1 = -t ] ; then
	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
	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 = -lm ] ; then
	lm=$2
	shift 2
    else
	Usage
    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 -m $betopts >> ${I}.sienax
${FSLDIR}/bin/fslmaths ${I}_brain -sub `$FSLDIR/bin/fslstats ${I}_brain -p 0` -mas ${I}_brain_mask ${I}_brain -odt float

echo ""                                                                        >> ${I}.sienax
echo "----------  register to stdairach 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}2std.mat >> ${I}.sienax 2>&1
${FSLDIR}/bin/avscale ${I}2std.mat ${FSLDIR}/etc/standard/avg152T1 > ${I}2std.avscale
xscale=`grep Scales ${I}2std.avscale | awk '{print $4}'`
yscale=`grep Scales ${I}2std.avscale | awk '{print $5}'`
zscale=`grep Scales ${I}2std.avscale | awk '{print $6}'`
vscale=`echo "10 k $xscale $yscale * $zscale * p"|dc -`
echo "VSCALING $vscale" >> ${I}.sienax

echo ""                                                                        >> ${I}.sienax
echo "----------  mask with std mask  ---------------------------------------" >> ${I}.sienax
${FSLDIR}/bin/convert_xfm -inverse -omat ${I}2std_inv.mat ${I}2std.mat
MASK=${FSLDIR}/etc/standard/avg152T1_brain_mask_dil
if [ "$stdroi" != "" ] ; then
    ${FSLDIR}/bin/fslmaths $MASK $stdroi ${I}_stdmaskroi
    MASK=${I}_stdmaskroi
fi
${FSLDIR}/bin/flirt -in $MASK -ref ${I}_brain -out ${I}_stdmask -applyxfm -init ${I}2std_inv.mat
${FSLDIR}/bin/fslmaths ${I}_stdmask -thr 0.5 -bin ${I}_stdmask
${FSLDIR}/bin/fslmaths ${I}_brain -mas ${I}_stdmask ${I}_stdmaskbrain

if [ $regional = 1 ] ; then
    ${FSLDIR}/bin/flirt -in ${FSLDIR}/etc/standard/avg152T1_strucseg_periph -ref ${I}_brain -out ${I}_stdmask_segperiph -applyxfm -init ${I}2std_inv.mat
    ${FSLDIR}/bin/fslmaths ${I}_stdmask_segperiph -thr 0.5 -bin ${I}_stdmask_segperiph
    ${FSLDIR}/bin/fslmaths ${FSLDIR}/etc/standard/avg152T1_strucseg -thr 4.5 -bin ${I}_tmpmask
    ${FSLDIR}/bin/flirt -in ${I}_tmpmask -ref ${I}_brain -out ${I}_stdmask_segvent -applyxfm -init ${I}2std_inv.mat
    ${FSLDIR}/bin/fslmaths ${I}_stdmask_segvent -thr 0.5 -bin ${I}_stdmask_segvent
    /bin/rm ${I}_tmpmask*
fi

echo ""                                                                        >> ${I}.sienax
echo "----------  segment tissue into types  --------------------------------" >> ${I}.sienax
if [ $nseg = 2 ] ; then
    ${FSLDIR}/bin/fast1 -c 2 $imtype -e -ov $segopts ${I}_stdmaskbrain >> ${I}.sienax 2>&1
    echo "" >> ${I}.sienax
    echo "----------  convert brain volume into normalised volume  --------------" >> ${I}.sienax
    S=`${FSLDIR}/bin/fslstats ${I}_stdmaskbrain_pve_1 -m -v`
    xa=`echo $S | awk '{print $1}'`
    xb=`echo $S | awk '{print $3}'`
    ubrain=`echo "2 k $xa $xb * 1 / p" | dc -`
    nbrain=`echo "2 k $xa $xb * $vscale * 1 / p" | dc -`
else
    if [ _$lm != _ ] ; then
	${FSLDIR}/bin/fslmaths $lm -bin -mul -1 -add 1 -mul ${I}_stdmaskbrain ${I}_stdmaskbrain -odt float
    fi

    ${FSLDIR}/bin/fast1 $imtype -e -ov $segopts ${I}_stdmaskbrain >> ${I}.sienax 2>&1

    if [ _$lm != _ ] ; then
	${FSLDIR}/bin/fslmaths $lm -bin -max ${I}_stdmaskbrain_pve_2 ${I}_stdmaskbrain_pve_2 -odt float
	${FSLDIR}/bin/fslmaths $lm -bin -mul 3 -max ${I}_stdmaskbrain_seg ${I}_stdmaskbrain_seg -odt int
    fi

    echo "" >> ${I}.sienax
    echo "----------  convert brain volume into normalised volume  --------------" >> ${I}.sienax
    echo "" >> ${I}.sienax
    echo "tissue             volume    unnormalised-volume" >> ${I}.sienax
    if [ $regional = 1 ] ; then
	${FSLDIR}/bin/fslmaths ${I}_stdmaskbrain_pve_1 -mas ${I}_stdmask_segperiph ${I}_stdmaskbrain_pve_1_segperiph -odt float
	S=`${FSLDIR}/bin/fslstats ${I}_stdmaskbrain_pve_1_segperiph -m -v`
	xa=`echo $S | awk '{print $1}'`
	xb=`echo $S | awk '{print $3}'`
	uxg=`echo "2 k $xa $xb * 1 / p" | dc -`
	xg=`echo "2 k $xa $xb * $vscale * 1 / p" | dc -`
	echo "pgrey              $xg $uxg (peripheral grey)" >> ${I}.sienax

	${FSLDIR}/bin/fslmaths ${I}_stdmaskbrain_pve_0 -mas ${I}_stdmask_segvent ${I}_stdmaskbrain_pve_0_segvent -odt float
	S=`${FSLDIR}/bin/fslstats ${I}_stdmaskbrain_pve_0_segvent -m -v`
	xa=`echo $S | awk '{print $1}'`
	xb=`echo $S | awk '{print $3}'`
	uxg=`echo "2 k $xa $xb * 1 / p" | dc -`
	xg=`echo "2 k $xa $xb * $vscale * 1 / p" | dc -`
	echo "vcsf               $xg $uxg (ventricular CSF)" >> ${I}.sienax
    fi
    S=`${FSLDIR}/bin/fslstats ${I}_stdmaskbrain_pve_1 -m -v`
    xa=`echo $S | awk '{print $1}'`
    xb=`echo $S | awk '{print $3}'`
    ugrey=`echo "2 k $xa $xb * 1 / p" | dc -`
    ngrey=`echo "2 k $xa $xb * $vscale * 1 / p" | dc -`
    echo "GREY               $ngrey $ugrey" >> ${I}.sienax
    S=`${FSLDIR}/bin/fslstats ${I}_stdmaskbrain_pve_2 -m -v`
    xa=`echo $S | awk '{print $1}'`
    xb=`echo $S | awk '{print $3}'`
    uwhite=`echo "2 k $xa $xb * 1 / p" | dc -`
    nwhite=`echo "2 k $xa $xb * $vscale * 1 / p" | dc -`
    echo "WHITE              $nwhite $uwhite" >> ${I}.sienax

    ubrain=`echo "2 k $uwhite $ugrey + 1 / p" | dc -`
    nbrain=`echo "2 k $nwhite $ngrey + 1 / p" | dc -`
fi

echo "BRAIN              $nbrain $ubrain" >> ${I}.sienax

${FSLDIR}/bin/overlay 1 1 -c ${I} -a ${I}_stdmaskbrain_seg 1.9 5 ${I}_render

if [ $regional = 1 ] ; then
    ${FSLDIR}/bin/overlay 1 1 -c ${I} -a ${I}_stdmaskbrain_pve_1_segperiph 0.3 0.7 ${I}_periph_render
    ${FSLDIR}/bin/overlay 1 1 -c ${I} -a ${I}_stdmaskbrain_pve_0_segvent   0.3 0.7 ${I}_vent_render
fi

if [ $debug = 0 ] ; then
    /bin/rm -f `$FSLDIR/bin/imglob ${I}_brain* ${I}_stdmask*`
    /bin/rm -f ${I}2std.avscale ${I}2std_inv.mat
fi

echo "$nbrain"