#!/bin/bash
#   Copyright (C) 2016 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 6.0 (c) 2018, The University of
#   Oxford (the "Software")
#
#   The Software remains the property of the Oxford University Innovation
#   ("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 Oxford
#   University Innovation ("OUI"), the technology transfer company of the
#   University, to negotiate a licence. Contact details are:
#   fsl@innovation.ox.ac.uk quoting Reference Project 9564, FSL.
export LC_ALL=C

set -e

######

if [ $# -lt 1 ] ; then
    echo "Usage: `basename $0` <lesionmask> <threshold> <manualmask> <saveoutput>"
    echo " "
    echo "the script calculates measures of overlap between the lesion mask generated by BIANCA and the manual mask and lesion volumes of both"
    echo "threshold = probability threshold that will be applied to BIANCA output before calculation of the overlap measures"
    echo "saveoutput = if set to 0 it will output the measures' names and values on the screen with the following format:"
    echo "DICEsimilarityindex(SI) FalseDiscoveryRate(FDR) FalseNegativeRatio(FNR) FDRcluster-level FNRcluster-level DetectionErrorRate(DER) OutlineErrorRate(OER) MeanTotalArea(MTA) LESIONmaskVolume MANUALmaskVolume"
    echo "             if set to 1 it will save only the values (in the same order) in a .txt file"
    echo "             the output file will be saved in the same folder of the lesion mask (BIANCA output) with the name Overlap_and_Volumes_<lesionmask>_<threshold>.txt"
    exit 1
fi

lesmaskfile=$1
thresh=$2
manualmaskfile=$3
saveoutput=$4

lesmask=$(remove_ext $(basename ${lesmaskfile}))
manualmask=$(realpath $(dirname ${manualmaskfile}))/$(remove_ext $(basename ${manualmaskfile}))

# cd into the directory containing the lesion mask. It will create all the files at this level. At the end it will go back to the folder where the command was called from.
origdir=`pwd`
cd $(dirname ${lesmaskfile})

logID=`echo $(date | awk '{print $1 $2 $3}' |  sed 's/://g')`
TMPVISDIR=`mktemp -d ./biancaoverlap_${logID}_${lesmask}_XXXXXX`

# # APPLY THRESHOLD AND CALCULATES VOLUMES
$FSLDIR/bin/fslmaths $lesmask -thr $thresh -bin ${TMPVISDIR}/lesmask
lesionvol=`$FSLDIR/bin/fslstats ${TMPVISDIR}/lesmask -V | awk '{ print $2 }'`
$FSLDIR/bin/fslmaths $manualmask -bin ${TMPVISDIR}/manualmask
manualvol=`$FSLDIR/bin/fslstats ${TMPVISDIR}/manualmask -V | awk '{ print $2 }'`

# # VOXEL-WISE INDICES
tpvox=`$FSLDIR/bin/fslstats ${TMPVISDIR}/lesmask -k ${TMPVISDIR}/manualmask -V | awk '{ print $1 }'`
posvox=`$FSLDIR/bin/fslstats ${TMPVISDIR}/lesmask -V | awk '{ print $1 }'`
fpvox=`echo "$posvox - $tpvox" | bc -l`
truvox=`$FSLDIR/bin/fslstats ${TMPVISDIR}/manualmask -V | awk '{ print $1 }'`
fnvox=`echo "$truvox - $tpvox" | bc -l`

# calculation of true positive ratio, true negative ratio, Dice index on the whole image (how good is the overlap overall)
FDR=`echo "scale=5; $fpvox / $posvox" | bc -l`
FNR=`echo "scale=5; $fnvox / $truvox" | bc -l`
sumvox=`echo "$posvox + $truvox" | bc -l`
Dice=`echo "scale=5; 2 * $tpvox / $sumvox" | bc -l`
# mean total area for OER and DER (Wack et al., 2012)
MTA=`echo "scale=5; 0.5 * $sumvox" | bc -l`

# # CLUSTER LEVEL INDICES
# generate FP lesion mask (all clusters in lesionmask that do not overlap with anything in manualmask)
$FSLDIR/bin/fsl-cluster -t 0.5 -i ${TMPVISDIR}/lesmask --connectivity=26 --oindex=${TMPVISDIR}/lesmask_idx --no_table
npos=`$FSLDIR/bin/fslstats ${TMPVISDIR}/lesmask_idx -P 100`
# lesmask_idx contains cluster indices (from 1 to npos) of automatically segmented lesions (positive clusters).
$FSLDIR/bin/fslmaths ${TMPVISDIR}/manualmask -thr 0.5 -mul ${TMPVISDIR}/lesmask_idx ${TMPVISDIR}/lesmask_idx_TP
# lesmask_idx_TP contains cluster indices only for voxels in the automatic lesion mask that have an overlap with the manual mask (from 1 up to npos)
maxv=`echo $npos + 0.5 | bc -l`
npos_bins=`echo $npos + 1 | bc`
histvals_pos=`$FSLDIR/bin/fslstats ${TMPVISDIR}/lesmask_idx_TP -H $npos_bins -0.5 $maxv | sed 's/\.0*//g'`
# one histogram bin for each positive cluster (index). If there are voxels with that value, there is overlap (TP), if not, that cluster is a false positive (it was in the automatic mask, but no overlap)

$FSLDIR/bin/fslmaths ${TMPVISDIR}/lesmask -mul 0 ${TMPVISDIR}/TPauto
n=0;
fptot=0
fptot_vox=0
for val in $histvals_pos ; do
    if [ $n -gt 0 ] ; then
	    if [ $val -eq 0 ] ; then
	        # FALSE POSITIVE CLUSTER COUNT
	        fptot=`echo "$fptot + 1" | bc`
	        # FALSE POSITIVE CLUSTERS VOXELS
	        $FSLDIR/bin/fslmaths ${TMPVISDIR}/lesmask_idx -thr $n -uthr $n -bin ${TMPVISDIR}/fpcluster_temp
	        fp_nvox=`$FSLDIR/bin/fslstats ${TMPVISDIR}/fpcluster_temp -V | awk '{ print $1 }'`
	        fptot_vox=`echo "$fptot_vox + $fp_nvox" | bc -l`
 	    else
	        # TRUE POSITIVE CLUSTERS - the map will contain true positive clusters from automatic segm
 	        $FSLDIR/bin/fslmaths ${TMPVISDIR}/lesmask_idx -thr $n -uthr $n -bin -add ${TMPVISDIR}/TPauto ${TMPVISDIR}/TPauto
 	    fi
    fi
    n=`echo $n + 1 | bc`;
done

# generate FN lesion mask (all clusters in manualmask that do not have any overlap with anything in lesionmask)
$FSLDIR/bin/fsl-cluster -t 0.5 -i ${TMPVISDIR}/manualmask --connectivity=26 --oindex=${TMPVISDIR}/manualmask_idx --no_table
ntrue=`$FSLDIR/bin/fslstats ${TMPVISDIR}/manualmask_idx -P 100`
# manualmask_idx contains cluster indices (from 1 to ntrue) of manual mask (true clusters)
$FSLDIR/bin/fslmaths ${TMPVISDIR}/lesmask -mul ${TMPVISDIR}/manualmask_idx ${TMPVISDIR}/manualmask_idx_TP
# manualmask_idx_TP  contains cluster indices only for voxels in the manaul mask that have an overlap with the automatic  mask (from 1 up to ntrue)
maxv=`echo $ntrue + 0.5 | bc -l`
ntrue_bins=`echo $ntrue + 1 | bc`
histvals_true=`$FSLDIR/bin/fslstats ${TMPVISDIR}/manualmask_idx_TP -H $ntrue_bins -0.5 $maxv | sed 's/\.0*//g'`
# one histogram bin for each true cluster (index). If there are voxels with that value, there is overlap, if not, that cluster is a false negative (it was in the manual mask, but no overlap)

$FSLDIR/bin/fslmaths ${TMPVISDIR}/manualmask -mul 0 ${TMPVISDIR}/TPmanual
n=0;
fntot=0
fntot_vox=0
# number of true positive voxels comes from manual mask (ground truth)
tptot=0
for val in $histvals_true ; do
    if [ $n -gt 0 ] ; then
	    if [ $val -eq 0 ] ; then
	        # FALSE NEGATIVE CLUSTERS COUNT
 	        fntot=`echo $fntot + 1 | bc`
	        # counts +1 cluster as false negative
	        # FALSE NEGATIVE CLUSTERS VOXELS
	        $FSLDIR/bin/fslmaths ${TMPVISDIR}/manualmask_idx -thr $n -uthr $n -bin ${TMPVISDIR}/fncluster_temp
	        fn_nvox=`$FSLDIR/bin/fslstats ${TMPVISDIR}/fncluster_temp -V | awk '{ print $1 }'`
	        fntot_vox=`echo "$fntot_vox + $fn_nvox" | bc -l`
	    else
	        tptot=`echo $tptot + 1 | bc`
 	        # accumulates TP clusters -  the map will contain true positive clusters from manual segm
 	        $FSLDIR/bin/fslmaths ${TMPVISDIR}/manualmask_idx -thr $n -uthr $n -bin -add ${TMPVISDIR}/TPmanual ${TMPVISDIR}/TPmanual
 	    fi
    fi
    n=`echo $n + 1 | bc`;
done

FP_clusters=$fptot
FN_clusters=$fntot
# Detection error (DE) sum of voxels in clusters detected only by one method (Wack et al., 2012)
DE=`echo "$fptot_vox + $fntot_vox" | bc -l`
# Detection error rate (DER) - DE divided by mean total area
DER=`echo "scale=5; $DE / $MTA" | bc -l`

FDR_cluster=`echo "scale=5; $FP_clusters / $npos" | bc -l`
FNR_cluster=`echo "scale=5; $FN_clusters / $ntrue" | bc -l`

# INDICES FROM TRUE POSITIVE CLUSTERS - TO CALCULATE OUTLINE ERROR
$FSLDIR/bin/fslmaths ${TMPVISDIR}/TPauto -sub ${TMPVISDIR}/TPmanual ${TMPVISDIR}/FPvox_TP_overlap
$FSLDIR/bin/fslmaths ${TMPVISDIR}/TPmanual -sub ${TMPVISDIR}/TPauto ${TMPVISDIR}/FNvox_TP_overlap
fpvox_overlap=`$FSLDIR/bin/fslstats ${TMPVISDIR}/FPvox_TP_overlap -l 0 -V | awk '{ print $1 }'`
fnvox_overlap=`$FSLDIR/bin/fslstats ${TMPVISDIR}/FNvox_TP_overlap -l 0 -V | awk '{ print $1 }'`
# Outline error (OE): within clusters deteced by both, sum of voxels detected only by one method (Wack et al., 2012)
OE=`echo "$fpvox_overlap + $fnvox_overlap" | bc -l`
# Outline error rate (OER) - OE divided by mean total area
OER=`echo "scale=5; $OE / $MTA" | bc -l`

if [ $saveoutput -eq 0 ] ; then
    echo SI $Dice FDR $FDR FNR $FNR FDR_clusters $FDR_cluster FNR_clusters $FNR_cluster DER $DER OER $OER MTA $MTA LESION_vol $lesionvol MANUAL_vol $manualvol
else
    echo $Dice $FDR $FNR $FDR_cluster $FNR_cluster $DER $OER $MTA $lesionvol $manualvol >> Overlap_and_Volumes_${lesmask}_${thresh}.txt
fi

rm -r $TMPVISDIR
cd ${origdir}
