#!/bin/sh

#   lesions - lesion analysis from T1&T2 or T1&FLAIR - not yet completed!!!!
#
#   John Vickers, Mark Jenkinson & Stephen Smith, FMRIB Image Analysis Group
#
#   Copyright (C) 2003-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 5.0 (c) 2012, The University of
#   Oxford (the "Software")
#   
#   The Software remains the property of the University of Oxford ("the
#   University").
#   
#   The Software is distributed "AS IS" under this Licence solely for
#   non-commercial use in the hope that it will be useful, but in order
#   that the University as a charitable foundation protects its assets for
#   the benefit of its educational and research purposes, the University
#   makes clear that no condition is made or to be implied, nor is any
#   warranty given or to be implied, as to the accuracy of the Software,
#   or that it will be suitable for any particular purpose or for use
#   under any specific conditions. Furthermore, the University disclaims
#   all responsibility for the use which is made of the Software. It
#   further disclaims any liability for the outcomes arising from using
#   the Software.
#   
#   The Licensee agrees to indemnify the University and hold the
#   University harmless from and against any and all claims, damages and
#   liabilities asserted by third parties (including claims for
#   negligence) which arise directly or indirectly from the use of the
#   Software or the sale of any products based on the Software.
#   
#   No part of the Software may be reproduced, modified, transmitted or
#   transferred in any form or by any means, electronic or mechanical,
#   without the express permission of the University. The permission of
#   the University is not required if the said reproduction, modification,
#   transmission or transference is done without financial return, the
#   conditions of this Licence are imposed upon the receiver of the
#   product, and all original and amended source code is included in any
#   transmitted product. You may be held legally responsible for any
#   copyright infringement that is caused or encouraged by your failure to
#   abide by these terms and conditions.
#   
#   You are not permitted under this Licence to use this Software
#   commercially. Use for which any financial return is received shall be
#   defined as commercial use, and includes (1) integration of all or part
#   of the source code or the Software into a product for sale or license
#   by or on behalf of Licensee to third parties or (2) use of the
#   Software or any derivative of it for research with the final aim of
#   developing software products for sale or license to a third party or
#   (3) use of the Software or any derivative of it for research with the
#   final aim of developing non-software products for sale or license to a
#   third party, or (4) use of the Software to provide any service to an
#   external organisation for which payment is received. If you are
#   interested in using the Software commercially, please contact Isis
#   Innovation Limited ("Isis"), the technology transfer company of the
#   University, to negotiate a licence. Contact details are:
#   innovation@isis.ox.ac.uk quoting reference DE/9564.

#### process options #############################################################################

if [ _$2 = _ ] ; then
    echo "Usage:  lesions <t1_root> <t2_root> [-f] [-hypo <hypoparam>]"
    echo "-f        :  second image is FLAIR not T2"
    exit
fi

A=$1
B=$2
shift 2
hypoparam=1.0;

nclasses=2
t2segopts="-t2"
# (no apriori on t2)
t2lesions=${B}_brain_pve_0

for opts in $@ ; do

    if [ $opts = -f ] ; then
        echo "processing second input as FLAIR image"
	nclasses=3
	t2segopts="-t2"
	t2lesions=${B}_brain_pve_2
	shift
    fi

    if [ $opts = -hypo ] ; then
        hypoparam=$2
	shift 2
    fi


done

##############################################################################################################################
#calculate minimum dimensom of t1 image for isotropic dilation################################################################



sizex=`fslval $A pixdim1` 
sizey=`fslval $A pixdim2`
sizez=`fslval $A pixdim3`

sizex=`echo sqrt \( $sizex ^ 2 \) | bc -l`
sizey=`echo sqrt \( $sizey ^ 2 \) | bc -l`
sizez=`echo sqrt \( $sizez ^ 2 \) | bc -l`

min=$sizex
num=`echo $min/$sizey|bc`
if [ $num -ge 1 ]; then
 min=$sizey;
fi
num=`echo $min/$sizez|bc`
if [ $num -ge 1 ]; then
 min=$sizez;
fi
num=`echo $min/1.0|bc`
if [ $num -ge 1 ]; then
 min=1.0;
fi


###############################################################################################################################

###############################################################################################################################

#### T1 #######################################################################################################################
# brain-extract t1 and run segmentation; pve_0 should be CSF only##############################################################
${FSLDIR}/bin/bet $A ${A}_brain -m
echo "Starting T1 Segmentation"
######First Line for installed fast, Second line for newer fast binary installed in ~/FastTwo/fast#############################
${FSLDIR}/bin/fast3.4 -e -ov -A ${A}_brain
echo "Finished T1 Segmentation"

#register t1 to standard space and invert transform###########################################################################
${FSLDIR}/bin/flirt -ref ${FSLDIR}/data/standard/MNI152_T1_2mm -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" originally 3 dils
${FSLDIR}/bin/fslmaths ${FSLDIR}/data/standard/MNI152_T1_2mm -thr 3.5 -bin -dilM -dilM -dilM 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/fslmaths 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/flirt -in ${A}_brain_pve_0 -ref ${A}_brain_pve_0 -out ${A}_brain_pve_0_iso -applyisoxfm $min
echo "Created T1 CSF"
##########original 5 dil outside mask 1 in#####################################################################################
#${FSLDIR}/bin/fslmaths ${A}_brain_pve_0_iso -dilM -dilM -dilM -dilM -dilM -sub strucseg_2_${A}_brain ${A}_grot
${FSLDIR}/bin/fslmaths ${A}_brain_pve_0_iso -kernel sphere 5.1 -dilF ${A}_grot 
${FSLDIR}/bin/fslmaths ${A}_grot -sub strucseg_2_${A}_brain ${A}_grot
#1 dil
#${FSLDIR}/bin/fslmaths ${A}_brain_pve_0_iso -dilM -mul strucseg_2_${A}_brain -add ${A}_grot ${A}_csf_dil
${FSLDIR}/bin/fslmaths ${A}_brain_pve_0_iso -kernel sphere 1.1 -dilF ${A}_csf_dil 
${FSLDIR}/bin/fslmaths ${A}_csf_dil -mul strucseg_2_${A}_brain -add ${A}_grot ${A}_csf_dil

/bin/rm -f ${A}_grot*
echo "Dilated T1 CSF"

#### T2 or FLAIR ##############################################################################################################

# register t2 to t1 and invert transform#######################################################################################
echo "Starting Registration of T1 & T2"
${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/fslmaths $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#################
echo "Starting T2 Segmentation"
######First Line for installed fast, Second line for newer fast binary installed in ~/FastTwo/fast#############################
${FSLDIR}/bin/fast3.4 $t2segopts -c $nclasses -e -ov ${B}_brain

echo "Finished T2 Segmentation"

#${FSLDIR}/bin/flirt -in $t2lesions -ref $A -out ${A}_lesion+CSF -applyxfm -init ${B}_to_${A}.mat
#${FSLDIR}/bin/fslmaths ${A}_lesion+CSF -mas ${A}_brain_mask ${A}_lesion+CSF
#${FSLDIR}/bin/flirt -in ${A}_lesion+CSF -ref ${A}_lesion+CSF -out ${A}_lesion+CSF_iso -applyisoxfm $min
${FSLDIR}/bin/flirt -in $t2lesions -ref $A -out ${A}_lesion+CSF_iso -applyisoxfm $min -init ${B}_to_${A}.mat


#### combine to remove CSF from CSF+lesions ##################################################################################
echo "Combining T1 and T2"
# combine masks and use to mask out CSF from t2-derived lesion+CSF probability; output lesion volume##########################
${FSLDIR}/bin/fslmaths ${A}_csf_dil -mul -1 -add 1 -mul ${A}_lesion+CSF_iso ${A}_continuouslesions_iso -odt float
#/bin/rm -f ${A}_csf*



###hypointense################################################################################################################
echo "Finding T1 hypointense lesions"

meanden0=`${FSLDIR}/bin/fslstats ${A}_brain_pve_0 -m`
${FSLDIR}/bin/fslmaths ${A}_brain_pve_0 -mul ${A}_brain ${A}_multemp -odt float
meannum0=`${FSLDIR}/bin/fslstats ${A}_multemp -m`
mean0=`echo $meannum0/$meanden0|bc -l`

meanden1=`${FSLDIR}/bin/fslstats ${A}_brain_pve_1 -m`
${FSLDIR}/bin/fslmaths ${A}_brain_pve_1 -mul ${A}_brain ${A}_multemp -odt float
meannum1=`${FSLDIR}/bin/fslstats ${A}_multemp -m`
mean1=`echo $meannum1/$meanden1|bc -l`

meanden2=`${FSLDIR}/bin/fslstats ${A}_brain_pve_2 -m`
${FSLDIR}/bin/fslmaths ${A}_brain_pve_2 -mul ${A}_brain ${A}_multemp -odt float
meannum2=`${FSLDIR}/bin/fslstats ${A}_multemp -m`
mean2=`echo $meannum2/$meanden2|bc -l`

#denom=`echo $mean0+$mean2|bc`
num2=`echo $mean1|bc -l`
num1=`echo $hypoparam/1.0|bc`
if [ $num1 -lt 1 ]; then
 num2=`echo $mean0-$hypoparam*$mean0+$mean1*$hypoparam|bc -l`;
fi
if [ $num1 -ge 1 ]; then
  num2=`echo $hypoparam*$mean2-$hypoparam*$mean1+2.0*$mean1-$mean2|bc -l`;
fi

echo "Mean of CSF   $mean0"
echo "Mean of Grey  $mean1"
echo "Mean of White $mean2"
echo "Hypointense upper threshold $num2"


${FSLDIR}/bin/flirt -in ${A}_continuouslesions_iso -ref $A -out ${A}_lesions -applyxfm
${FSLDIR}/bin/fslmaths $A -uthr $num2 -bin ${A}_hypo
${FSLDIR}/bin/fslmaths ${A}_lesions -mas ${A}_hypo ${A}_hypolesions

echo `${FSLDIR}/bin/fslstats ${A}_lesions -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}_brain_* ${B}_brain_* ${A}_lesions_render*  ${A}_lesions_render* ${A}_grot* ${B}_grot*
#/bin/rm -f ${A}_lesions_thr*
#/bin/rm -f ${A}_lesions*
#/bin/rm -f ${A}_csf* 
####JV#############################################################################################################################