diff --git a/performanceanalysis b/performanceanalysis new file mode 100644 index 0000000000000000000000000000000000000000..6fc22622813d0a61f579d05c00dbbec7f4371583 --- /dev/null +++ b/performanceanalysis @@ -0,0 +1,456 @@ +#!/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 +# +# Originally developed at FMRIB (Oxford Centre for Functional Magnetic Resonance +# Imaging of the Brain), Department of Clinical Neurology, Oxford +# University, Oxford, UK +# +# Modified at Department of Medical Physics and Biomedical Engineering +# Modified (2019) +# University, UCL, UK +# +# +# This software/script calculates all the following important performance metrics: +# Dice Similarity Index(SI) +# Sensitivity +# Precision +# Specificity +# Jaccardindex +# Conformity +# HausdorffDistance +# Overlapcoefficient +# Tanimotocoefficient +# Matthewscorrelationcoefficient +# Accuracy +# Voxel-level false discovery rate(FDR) +# Voxel-level false negative ratio(FNR) +# Cluster-level FDR +# Cluster-level FNR +# Detection error rate(DER) +# utline error rate(OER) +# Mean Total Area(MTA) +# Volume of segmentation +# Volume of manual mask + + + + +export LC_ALL=C +#set -e +#set -x + +###### + +if [ $# -lt 1 ] ; then + echo "Usage: `basename $0` " + echo " " + echo "For example: + performanceanalysis *lesionmask*.nii.gz 0.5 *manualmask*.nii.gz euclidean mean symmetric 1" + echo " " + echo "The script calculates measures of overlap between the lesion mask generated by any algorithm and the manual mask and lesion volumes of both" + echo " " + echo "Threshold = probability threshold that will be applied to output before calculation of the overlap measures" + echo " " + echo "Saveoutput = if set to 0 it will output the measures' names and values on the screen with the following format:" + echo " +# Dice Similarity Index(SI) +# Sensitivity +# Precision +# Specificity +# Jaccardindex +# Conformity +# HausdorffDistance +# Overlapcoefficient +# Tanimotocoefficient +# Matthewscorrelationcoefficient +# Accuracy +# Voxel-level false discovery rate(FDR) +# Voxel-level false negative ratio(FNR) +# Cluster-level FDR +# Cluster-level FNR +# Detection error rate(DER) +# utline error rate(OER) +# Mean Total Area(MTA) +# Volume of segmentation +# Volume of manual mask" + echo " " + echo "-If Saveoutput set to 1 it will save only the values (in the same order) in a .txt filethe output file will be saved in the same folder of the lesion mask with the name Overlap_and_Volumes__.txt" + exit 1 +fi + + +function progress_bar() { + bar="" + total=10 + [[ -z $1 ]] && input=0 || input=${1} + x="##" + for i in `seq 1 10`; do + if [ $i -le $input ] ;then + bar=$bar$x + else + bar="$bar " + fi + done + #pct=$((200*$input/$total % 2 + 100*$input/$total)) + pct=$(($input*10)) + echo -ne "Progress : [ ${bar} ] (${pct}%) \r" + sleep 1 + if [ $input -eq 10 ] ;then + echo -ne '\n' + fi + +} + + + +lesmaskfile=$1 +lesmaskname=`basename ${lesmaskfile} .nii.gz` +lesmaskdir=`dirname ${lesmaskfile} ` +pushd $lesmaskdir > /dev/null +lesmaskdir=`pwd` +popd > /dev/null + + +thresh=$2 + +manualmaskfile=$3 +manualmaskname=`basename ${manualmaskfile} .nii.gz` +manualmaskdir=`dirname ${manualmaskfile} ` +pushd $manualmaskdir > /dev/null +manualmaskdir=`pwd` +popd > /dev/null +manualmask=${manualmaskdir}/${manualmaskname} +lesmask=${manualmaskdir}/${lesmaskname} + +echo "Identifying the inputs required for the calculation of Hausdorff Distance" + +if [ $# -gt 3 ] ; then +ems="$4 " +else +ems="euclidean " +fi + +if [ $# -gt 4 ] ; then +ems+=" $5" +fct="$5" +else +ems+=" mean" +fct="mean" +fi + +if [ $# -gt 5 ] ; then +ems+=" $6" +symt="$6" +else +ems+=" symmetric" +symt="symmetric" +fi + +if [ $# -gt 6 ] ; then +ems+=" $7" +else +ems+=" 1.0" +fi + +if [ $# -gt 7 ] ; then +ems+=" $8" +else +ems+=" 2.0" +fi + +echo "Distance function, factor, symmetry, weight and p_norm:" $ems + + +# echo "dist:" $dist + +# echo "input arg:" $ems +# echo "fact:" $fact + +if [ $# -lt 9 ] ; then +saveoutput=1 +else +saveoutput=$9 +fi + +# euclidean mean symmetric + + +# 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 ${lesmaskdir} + +logID=`echo $(date | awk '{print $1 $2 $3}' | sed 's/://g')` +# TMPVISDIR=`mktemp -d ./biancaoverlap_${logID}_${lesmask}_XXXXXX` +RESULTSDIR=${PWD%/*} + + + + +echo "The measurement has been started..." +echo "lesmask...", $lesmask +echo "manualmask...", $manualmask + +# # APPLY THRESHOLD AND CALCULATES VOLUMES +$FSLDIR/bin/fslmaths $lesmask -thr $thresh -bin lesmask +$FSLDIR/bin/fslmaths lesmask -sub 1 -abs -bin lesmaskn + + +lesionvol=`$FSLDIR/bin/fslstats lesmask -V | awk '{ print $2 }'` +lesionvolm=`$FSLDIR/bin/fslstats lesmask -V | awk '{ print $1 }'` +$FSLDIR/bin/fslmaths $manualmask -bin manualmask +manualvol=`$FSLDIR/bin/fslstats manualmask -V | awk '{ print $2 }'` +manualvolm=`$FSLDIR/bin/fslstats manualmask -V | awk '{ print $1 }'` +$FSLDIR/bin/fslmaths manualmask -sub 1 -abs -bin manualmaskn +# $FSLDIR/bin/fslmaths lesmaskn -mul manualmaskn -bin TN + +$FSLDIR/bin/fslmaths lesmask -mul manualmask -bin intersection +$FSLDIR/bin/fslmaths lesmask -add manualmask -bin union +$FSLDIR/bin/fslmaths lesmask -sub manualmask -bin lesmask_manualmask +$FSLDIR/bin/fslmaths manualmask -sub lesmask -bin manualmask_lesmask +AND=`$FSLDIR/bin/fslstats intersection -V | awk '{print $1}'` +OR=`$FSLDIR/bin/fslstats union -V | awk '{print $1}'` +$FSLDIR/bin/fslmaths lesmask -edge edge_mask +$FSLDIR/bin/fslmaths manualmask -edge edge_manual + +lm=`$FSLDIR/bin/fslstats lesmask_manualmask -V | awk '{ print $1 }'` +ml=`$FSLDIR/bin/fslstats manualmask_lesmask -V | awk '{ print $1 }'` + + + +if [ "$lesionvolm =< $manualvolm" ] ; then + +smallerset=`echo "$lesionvolm" | bc -l` +else +smallerset=`echo "$manualvolm" | bc -l` +fi + + +# $FSLDIR/bin/fslmaths ${TMPVISDIR}/manualmask -sub 1 -abs -bin ${TMPVISDIR}/euclidean_manual +# $FSLDIR/bin/fslmaths ${TMPVISDIR}/manualmask -mul ${TMPVISDIR}/manualmask -sub 1 -add ${TMPVISDIR}/manualmask -mul 2 -sqrt -bin ${TMPVISDIR}/euclidean_manual +# $FSLDIR/bin/fslmaths ${TMPVISDIR}/lesmask -sub 1 -abs -bin ${TMPVISDIR}/euclidean_mask +# $FSLDIR/bin/fslmaths ${TMPVISDIR}/lesmask -mul ${TMPVISDIR}/lesmask -sub 1 -add ${TMPVISDIR}/lesmask -mul 2 -sqrt -bin ${TMPVISDIR}/euclidean_mask +# $FSLDIR/bin/fslmaths ${TMPVISDIR}/edge_manual -mul ${TMPVISDIR}/edge_mask -sub 1 -abs -bin ${TMPVISDIR}/non_overlap +# $FSLDIR/bin/fslmaths ${TMPVISDIR}/euclidean_manual -mul ${TMPVISDIR}/non_overlap -mul ${TMPVISDIR}/edge_mask -abs -thr 0.5 ${TMPVISDIR}/dist_mask_manual +# $FSLDIR/bin/fslmaths ${TMPVISDIR}/euclidean_mask -mul ${TMPVISDIR}/non_overlap -mul ${TMPVISDIR}/edge_manual -abs -thr 0.5 ${TMPVISDIR}/dist_manual_mask + + # hd=`python Un.py edge_mask.nii.gz edge_manual.nii.gz` +progress_bar 1 +# hd=`python universal_HD.py edge_mask.nii.gz edge_manual.nii.gz euclidean mean symmetric| awk '{ print $1 }'` +hd=`python universal_HD.py edge_mask.nii.gz edge_manual.nii.gz $ems| awk '{ print $1 }'` + # echo "python universal_HD.py edge_mask.nii.gz edge_manual.nii.gz $ems" + # hd=$(python universal_HD.py edge_mask.nii.gz edge_manual.nii.gz euclidean mean symmetric) + # 11522097 11534336 +progress_bar 2 +# # VOXEL-WISE INDICES +tpvox=`$FSLDIR/bin/fslstats lesmask -k manualmask -V | awk '{ print $1 }'` +tnvox=`$FSLDIR/bin/fslstats lesmaskn -k manualmaskn -V | awk '{ print $1 }'` +posvox=`$FSLDIR/bin/fslstats lesmask -V | awk '{ print $1 }'` + +fpvox=`echo "$posvox - $tpvox" | bc -l` +truvox=`$FSLDIR/bin/fslstats manualmask -V | awk '{ print $1 }'` + # tnvox=`$FSLDIR/bin/fslstats ${TMPVISDIR}/TN -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) +progress_bar 3 +FPR=`echo "scale=5; $fpvox / $posvox" | bc -l` + +FNR=`echo "scale=5; $fnvox / $truvox" | bc -l` + +sumvox=`echo "$posvox + $truvox" | bc -l` +progress_bar 4 +Dice=`echo "scale=5; 2 * $tpvox / $sumvox" | bc -l` + +tmp1=`echo "$tpvox + $fnvox" | bc -l` +tmp2=`echo "$tpvox + $fpvox" | bc -l` +tmp3=`echo "$tpvox + $tnvox" | bc -l` +tmp4=`echo "$tpvox + $tnvox + $fnvox + $fpvox" | bc -l` +tmp5=`echo "$tnvox + $fpvox" | bc -l` +progress_bar 5 +Sensitivity=`echo "scale=5; 100 * $tpvox / $tmp1 " | bc -l` +Precision=`echo "scale=5; 100 * $tpvox / $tmp2 " | bc -l` + +Accuracy=`echo "scale=5; 100 * $tmp3 / $tmp4 " | bc -l` +Specificity=`echo "scale=5; 100 * $tnvox / $tmp5 " | bc -l` + + +progress_bar 6 +Jaccardindex=`echo "scale=5; $AND / $OR " | bc -l` +Conformity=`echo "scale=5; 100 *( 1 - ( ( $fpvox + $fnvox ) / $tpvox ) )" | bc -l` + +HausdorffDistance=`echo "scale=5; $hd" | bc -l` + +Overlapcoefficient=`echo "scale=5; $AND / $smallerset " | bc -l` +Tanimotocoefficient=`echo "scale=5; $AND / ( $AND + $lm + $ml ) " | bc -l` + # A coefficient of +1 represents a perfect prediction, 0 no better than random prediction and −1 indicates total + # disagreement between prediction and observation + + + +Sqrt=`echo "scale=5; sqrt ( ( $tpvox + $fpvox ) * ( $tpvox + $fnvox ) * ( $tnvox + $fpvox ) * ( $tnvox + $fnvox ) )" | bc -l` + +# testing=`echo "scale=5; ( $tpvox + $fpvox ) * ( $tpvox + $fnvox ) * ( $tnvox + $fpvox ) * ( $tnvox + $fnvox ) " | bc -l` + +# echo "tpvox:" $tpvox +# echo "fpvox:" $fpvox +# echo "tnvox:" $tnvox +# echo "fnvox:" $fnvox +# echo "Sqrt:" $Sqrt +# echo "testing:" $testing +tmpx=`echo "( ( $tpvox * $tnvox ) - ( $fpvox * $fnvox ) )" | bc -l` +# echo "( ( $tpvox * $tnvox ) - ( $fpvox * $fnvox ) ):" $this +msc=`echo " $tmpx / $Sqrt " | bc -l` +# echo "thisx:" $thisx + + +Matthewscorrelationcoefficient=`echo "scale=5; $msc " | 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/cluster -t 0.5 -i lesmask --connectivity=26 --oindex=lesmask_idx --no_table +npos=`$FSLDIR/bin/fslstats lesmask_idx -P 100` + # lesmask_idx contains cluster indices (from 1 to npos) of automatically segmented lesions (positive clusters). +$FSLDIR/bin/fslmaths manualmask -thr 0.5 -mul lesmask_idx 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 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 lesmask -mul 0 TPauto +n=0; +fptot=0 +fptot_vox=0 +progress_bar 7 +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 lesmask_idx -thr $n -uthr $n -bin fpcluster_temp + fp_nvox=`$FSLDIR/bin/fslstats 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 lesmask_idx -thr $n -uthr $n -bin -add TPauto 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/cluster -t 0.5 -i manualmask --connectivity=26 --oindex=manualmask_idx --no_table + ntrue=`$FSLDIR/bin/fslstats manualmask_idx -P 100` + # manualmask_idx contains cluster indices (from 1 to ntrue) of manual mask (true clusters) + $FSLDIR/bin/fslmaths lesmask -mul manualmask_idx 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 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 manualmask -mul 0 TPmanual +n=0; +fntot=0 +fntot_vox=0 +# number of true positive voxels comes from manual mask (ground truth) +progress_bar 8 +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 manualmask_idx -thr $n -uthr $n -bin fncluster_temp + fn_nvox=`$FSLDIR/bin/fslstats 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 manualmask_idx -thr $n -uthr $n -bin -add TPmanual TPmanual + fi + fi + n=`echo $n + 1 | bc`; + done +progress_bar 9 +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` + +FPR_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 TPauto -sub TPmanual FPvox_TP_overlap +$FSLDIR/bin/fslmaths TPmanual -sub TPauto FNvox_TP_overlap +fpvox_overlap=`$FSLDIR/bin/fslstats FPvox_TP_overlap -l 0 -V | awk '{ print $1 }'` +fnvox_overlap=`$FSLDIR/bin/fslstats 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` + +Dice=`printf "%.5f" ${Dice}` +Sensitivity=`printf "%.5f" ${Sensitivity}` +Precision=`printf "%.5f" ${Precision}` +Specificity=`printf "%.5f" ${Specificity}` +Jaccardindex=`printf "%.5f" ${Jaccardindex}` +Conformity=`printf "%.5f" ${Conformity}` +HausdorffDistance=`printf "%.5f" ${HausdorffDistance}` +Accuracy=`printf "%.5f" ${Accuracy}` + +FPR=`printf "%.5f" ${FPR}` +FNR=`printf "%.5f" ${FNR}` +FPR_cluster=`printf "%.5f" ${FPR_cluster}` +FNR_cluster=`printf "%.5f" ${FNR_cluster}` +DER=`printf "%.5f" ${DER}` +OER=`printf "%.5f" ${OER}` +MTA=`printf "%.5f" ${MTA}` +Matthewscorrelationcoefficient=`printf "%.5f" ${Matthewscorrelationcoefficient}` +Tanimotocoefficient=`printf "%.5f" ${Tanimotocoefficient}` +Overlapcoefficient=`printf "%.5f" ${Overlapcoefficient}` +lesionvol=`printf "%.5f" ${lesionvol}` +manualvol=`printf "%.5f" ${manualvol}` + + + +if [ $saveoutput -eq 0 ] ; then +echo -e "Dice Similarity Index(SI): ${Dice} \nSensitivity: ${Sensitivity} \nPrecision: ${Precision} \nSpecificity: ${Specificity} \nJaccardindex: ${Jaccardindex} \nConformity: ${Conformity} \nMeanHausdorffDistance: ${HausdorffDistance} \nOverlapcoefficient: ${Overlapcoefficient} \nTanimotocoefficient: ${Tanimotocoefficient} \nMatthewscorrelationcoefficient: ${Matthewscorrelationcoefficient} \nAccuracy: ${Accuracy} \nVoxel-level false discovery rate(FDR): ${FPR} \nVoxel-level false negative ratio(FNR): ${FNR} \nCluster-level FDR: ${FPR_cluster} \nCluster-level FNR: ${FNR_cluster} \nDetection error rate(DER): ${DER} \nOutline error rate(OER): ${OER} \nMean Total Area(MTA): ${MTA} \nVolume of segmentation: ${lesionvol} \nVolume of manual mask: ${manualvol}" + +else + + +# echo SI $Dice FPR $FPR FNR $FNR FPR_clusters $FPR_cluster FNR_clusters $FNR_cluster DER $DER OER $OER MTA $MTA LESION_vol $lesionvol MANUAL_vol # $manualvol >> Overlap_and_Volumes_${lesmask}_${thresh}.txt +echo -e "The measure of overlap between the lesion mask($1) and the manual mask($3) and lesion volumes of both" >> Overlap_and_Volumes_result.txt +echo -e "" >> Overlap_and_Volumes_result.txt +echo -e "" >> Overlap_and_Volumes_result.txt +echo -e "Dice Similarity Index(SI): ${Dice} \nSensitivity: ${Sensitivity} \nPrecision: ${Precision} \nSpecificity: ${Specificity} \nJaccardindex: ${Jaccardindex} \nConformity: ${Conformity} \n${fct}${symt}HausdorffDistance: ${HausdorffDistance} \nOverlapcoefficient: ${Overlapcoefficient} \nTanimotocoefficient: ${Tanimotocoefficient} \nMatthewscorrelationcoefficient: ${Matthewscorrelationcoefficient} \nAccuracy: ${Accuracy} \nVoxel-level false discovery rate(FDR): ${FPR} \nVoxel-level false negative ratio(FNR): ${FNR} \nCluster-level FDR: ${FPR_cluster} \nCluster-level FNR: ${FNR_cluster} \nDetection error rate(DER): ${DER} \nOutline error rate(OER): ${OER} \nMean Total Area(MTA): ${MTA} \nVolume of segmentation: ${lesionvol} \nVolume of manual mask: ${manualvol}" >> Overlap_and_Volumes_result.txt + + +echo -e "${Dice} ${Sensitivity} ${Precision} ${Specificity} ${Jaccardindex} ${Conformity} ${HausdorffDistance} ${Overlapcoefficient} ${Tanimotocoefficient} ${Matthewscorrelationcoefficient} ${Accuracy} ${FPR} ${FNR} ${FPR_cluster} ${FNR_cluster} ${DER} ${OER} ${MTA} ${lesionvol} ${manualvol}" >> ${RESULTSDIR}/final_result.txt + + +# echo $Dice $FPR $FNR $FPR_cluster $FNR_cluster $DER $OER $MTA $lesionvol $manualvol >> Overlap_and_Volumes_${lesmask}_${thresh}.txt +fi +progress_bar 10 +# rm -r $TMPVISDIR +rm -rf edge_mask.nii.gz edge_manual.nii.gz fncluster_temp.nii.gz FNvox_TP_overlap.nii.gz fpcluster_temp.nii.gz FPvox_TP_overlap.nii.gz +rm -rf intersection.nii.gz lesmask.nii.gz lesmask_idx.nii.gz lesmask_idx_TP.nii.gz lesmask_manualmask.nii.gz lesmaskn.nii.gz manualmask.nii.gz +rm -rf manualmask_idx.nii.gz manualmask_idx_TP.nii.gz manualmask_lesmask.nii.gz manualmaskn.nii.gz TPauto.nii.gz TPmanual.nii.gz union.nii.gz +rm -rf performanceanalysis universal_HD.py +cd ${origdir} +echo "The measurement has been completed"