#!/bin/bash

VERSION="0.0"

if [[ ! -d "$ANTSPATH" ]];
  then
    echo We can\'t find the ANTs path -- does not seem to exist.  Please \(re\)define \$ANTSPATH in your environment.
    exit 1
  fi

function Usage {
    cat <<USAGE

`basename $0` creates feature images for tumor segmentation using random decision forests.
   See <reference goes here>

Usage:

`basename $0` -d imageDimension
              -a inputImage(s)
              -t symmetricTemplate(s)
              -c clusterCenters
              -n imageNames
              -r neighgorhoodRadius
              -s smoothSigma
              -f differencePair
              -x maskImage
              -o outputPrefix

Example:

  bash $0 -d 3 -a t1.nii.gz -t symmetricTemplate.nii.gz -o output

Required arguments:

     -d:  image dimension                       2 or 3 (for 2- or 3-dimensional image)
     -a:  input image                           Anatomical image(s) comprised of flair, t2, t1, or t1 contrast.
                                                It is assumed that the following preprocessing steps have
                                                already been performed:
                                                   * brain extraction
                                                   * bias correction
                                                   * intensity normalized (i.e. histogram matching and intensity truncation)
     -c:  cluster centers                       Array describing the intensity centers of the intensity normalized images.
                                                Need one for each input image.  Should be of the form:  e.g. 0.14x0.57x0.37x0.83x0.95
                                                (for 7 classes: csf, gm, wm, edema, tumor, core, enhancement).  Note that either the cluster
                                                centers are specified (for testing) or the truth labels (for training) but not both. If these
                                                aren't specified, the number of classes needs to be specified ('-b') option
     -b:  number of clusters                    If -c is specified, this option is not needed.
     -g:  truth labels                          Truth labels.  Note that either the cluster centers are specified (for testing)
                                                or the truth labels (for training) but not both.
     -t:  symmetric anatomical templates        Symmetric templates.  Need to be specified in the same order as
                                                the input anatomical images.
     -z:  0,1                                   enable otb texture features
     -x:  mask image                            Mask image defining the region of interest.
     -o:  output prefix                         The following images are created:
                                                  * ${OUTPUT_PREFIX}N4Corrected.${OUTPUT_SUFFIX}
                                                  * ${OUTPUT_PREFIX}Segmentation.${OUTPUT_SUFFIX}
                                                  * ${OUTPUT_PREFIX}SegmentationPosteriors.${OUTPUT_SUFFIX}

Optional arguments:

     -r   radius                                Neighborhood voxel radius for local statistics image (default = 2)
     -s   smoothingSigma                        Standard deviation of gaussian smoothing kernel (in physical space) for
                                                the following images:
                                                  * symmetric template difference
                                                  * contralateral difference
     -p:  Brain segmentation priors             Tissue *probability* priors. Specified using c-style formatting, e.g.
                                                -p labelsPriors%02d.nii.gz.
     -f:  difference pair                       pair of indices \"i.e. 3x1\" to create difference image \"image[3] - image[1]\"
     -l:  tumor core label                      used to create distance feature map (default = 5)
     -n   imageNames                            used in the naming of the images (otherwise, labeled IMAGE0, IMAGE1, etc)
     -u   symmetric template mask

USAGE
    exit 1
}

echoParameters() {
    cat <<PARAMETERS

    Using createFeatureImages with the following arguments:
      image dimension         = ${DIMENSION}
      anatomical image        = ${ANATOMICAL_IMAGES[@]}
      symmetric templates     = ${SYMMETRIC_TEMPLATES[@]}
      symmetric template mask = ${SYMMETRIC_TEMPLATE_MASK}
      cluster centers         = ${CLUSTER_CENTERS[@]}
      image names             = ${IMAGE_NAMES[@]}
      radii                   = ${RADII[@]}
      smoothing sigma         = ${SMOOTHING_SIGMA}
      priors                  = ${SEGMENTATION_PRIOR}
      difference pairs        = ${DIFFERENCE_PAIRS[@]}
      output prefix           = ${OUTPUT_PREFIX}

PARAMETERS
}


#    local  myresult='some value'
#    echo "$myresult"

# Echos a command to both stdout and stderr, then runs it
function logCmd() {
  cmd="$*"
  echo "BEGIN >>>>>>>>>>>>>>>>>>>>"
  echo $cmd
  logCmdOutput=$( $cmd | tee /dev/tty )
  echo "END   <<<<<<<<<<<<<<<<<<<<"
  echo
  echo
}

################################################################################
#
# Main routine
#
################################################################################

HOSTNAME=`hostname`
DATE=`date`

CURRENT_DIR=`pwd`/
OUTPUT_DIR=${CURRENT_DIR}/tmp$RANDOM/
OUTPUT_PREFIX=${OUTPUT_DIR}/tmp
OUTPUT_SUFFIX="nii.gz"
LABELFILE=Code/dfltlabels.txt

DIMENSION=3

ANATOMICAL_IMAGES=()
NORMALIZED_IMAGES=()
RESCALED_IMAGES=()
RAWIMAGES_IMAGES=()
SYMMETRIC_TEMPLATE=()
SYMMETRIC_TEMPLATE_MASK=""
CLUSTER_CENTERS=()
NUMBER_OF_LABELS=7
IMAGE_NAMES=()
DIFFERENCE_PAIRS=()

SEGMENTATION_PRIOR=""

MASK_IMAGE=""

RADII=()
SMOOTHING_SIGMA=0
CORE_LABEL=5
CENTROIDSLICE=1
ENABLEOTB=0

################################################################################
#
# Programs and their parameters
#
################################################################################

if [[ $# -lt 3 ]] ; then
  Usage >&2
  exit 1
else
  while getopts "a:b:c:d:e:f:g:h:l:n:o:p:r:s:t:u:x:z:" OPT
    do
      case $OPT in
          a) #anatomical image
       ANATOMICAL_IMAGES[${#ANATOMICAL_IMAGES[@]}]=$OPTARG
       ;;
          b)
       NUMBER_OF_LABELS=$OPTARG
       ;;
          c) # cluster centers
       CLUSTER_CENTERS[${#CLUSTER_CENTERS[@]}]=$OPTARG
       ;;
          d) #dimensions
       DIMENSION=$OPTARG
       if [[ ${DIMENSION} -gt 4 || ${DIMENSION} -lt 2 ]];
         then
           echo " Error:  ImageDimension must be 2, 3, or 4 "
           exit 1
         fi
       ;;
          e)
       CENTROIDSLICE=$OPTARG
       ;;
          f)
       DIFFERENCE_PAIRS[${#DIFFERENCE_PAIRS[@]}]=$OPTARG
       ;;
          g)
       TRUTH_LABELS=$OPTARG
       ;;
          h) #help
       Usage >&2
       exit 0
       ;;
          l)
       CORE_LABEL=$OPTARG
       ;;
          n)
       IMAGE_NAMES[${#IMAGE_NAMES[@]}]=$OPTARG
       ;;
          o) #output prefix
       OUTPUT_PREFIX=$OPTARG
       ;;
          p) #brain segmentation label prior image
       SEGMENTATION_PRIOR=$OPTARG
       ;;
          r)
       RADII[${#RADII[@]}]=$OPTARG
       ;;
          s)
       SMOOTHING_SIGMA=$OPTARG
       ;;
          t)
       SYMMETRIC_TEMPLATES[${#SYMMETRIC_TEMPLATES[@]}]=$OPTARG
       ;;
          u)
       SYMMETRIC_TEMPLATE_MASK=$OPTARG
       ;;
          x)
       MASK_IMAGE=$OPTARG
       ;;
          z)
       ENABLEOTB=$OPTARG
       ;;
          *) # getopts issues an error message
       echo "ERROR:  unrecognized option -$OPT $OPTARG"
       exit 1
       ;;
      esac
  done
fi

################################################################################
#
# Preliminaries:
#  1. Check existence of inputs
#  2. Figure out output directory and mkdir if necessary
#
################################################################################

for (( i = 0; i < ${#ANATOMICAL_IMAGES[@]}; i++ ))
  do
    if [[ ! -f ${ANATOMICAL_IMAGES[$i]} ]];
      then
        echo "The specified image \"${ANATOMICAL_IMAGES[$i]}\" does not exist."
      fi
    if [[ -z ${IMAGE_NAMES[$i]} ]];
      then
        IMAGE_NAMES[$i]=IMAGE${i}
      fi
  done

for (( i = 0; i < ${#SYMMETRIC_TEMPLATES[@]}; i++ ))
  do
    if [[ ! -f ${SYMMETRIC_TEMPLATES[$i]} ]];
      then
        echo "The specified template image \"${SYMMETRIC_TEMPLATES[$i]}\" does not exist."
        exit 1
      fi
  done

if [[ ! -f ${MASK_IMAGE} ]];
  then
    echo "The specified mask image \"${MASK_IMAGE}\" does not exist."
    exit 1
  fi
echo centroid $CENTROIDSLICE
echo OTB $ENABLEOTB

#OUTPUT_IMAGE=${OUTPUT_PREFIX}MASKERODE.${OUTPUT_SUFFIX}
#if [[ ! -f ${OUTPUT_IMAGE} ]];
#  then
#    logCmd c3d ${MASK_IMAGE} -erode 1 3x3x0vox -o $OUTPUT_IMAGE
#  fi
#MASK_IMAGE=$OUTPUT_IMAGE
#
#PNG_IMAGE=${OUTPUT_PREFIX}MASKERODE.png
#if [[ ! -f ${PNG_IMAGE} ]];
#  then
#    logCmd c3d $OUTPUT_IMAGE -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc $PNG_IMAGE
#  fi


if [[ ${#ANATOMICAL_IMAGES[@]} -lt 3 ]];
  then
      echo "Expecting Images: (Pre) Art Ven Del "
      exit 1
  fi
#if [[ ${#ANATOMICAL_IMAGES[@]} -ne ${#SYMMETRIC_TEMPLATES[@]} ]];
#  then
#      echo "The number of symmetric templates does not match the number of anatomical images."
#      exit 1
#  fi


if [[ ${#CLUSTER_CENTERS[@]} -gt 0 ]];
  then
   if [[ ${#ANATOMICAL_IMAGES[@]} -ne ${#CLUSTER_CENTERS[@]} ]]
     then
         echo "The number of cluster center arrays does not match the number of anatomical images."
         exit 1
     fi

    CLUSTERS=( `echo ${CLUSTER_CENTERS[0]} | tr 'x' ' '` )
    NUMBER_OF_LABELS=${#CLUSTERS[@]}
    for (( i = 1; i < ${#CLUSTER_CENTERS[@]}; i++ ))
      do
        CLUSTERS=( `echo ${CLUSTER_CENTERS[i]} | tr 'x' ' '` )
        if [[ ${#CLUSTERS[@]} -ne NUMBER_OF_LABELS ]];
          then
            echo "The number of labels is not equal across the cluster center arrays."
            exit 1
          fi
      done
  fi

OUTPUT_DIR=${OUTPUT_PREFIX%\/*}
if [[ ! -d $OUTPUT_DIR ]];
  then
    echo "The output directory \"$OUTPUT_DIR\" does not exist. Making it."
    mkdir -p $OUTPUT_DIR
  fi

echoParameters >&2

echo "---------------------  Running `basename $0` on $HOSTNAME  ---------------------"

time_start=`date +%s`

################################################################################
#
# Normalize images -> truncate intensity and denoise
#
################################################################################

for (( i = 0; i < ${#ANATOMICAL_IMAGES[@]}; i++ ))
  do
    OUTPUT_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_DENOISE.${OUTPUT_SUFFIX}
    RAWIMAGE_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_RAWIMAGE.${OUTPUT_SUFFIX}
    RESCALE_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_RESCALE.${OUTPUT_SUFFIX}
    GRADIENT_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_GRADIENT.${OUTPUT_SUFFIX}

    NORMALIZED_IMAGES[${#NORMALIZED_IMAGES[@]}]=$OUTPUT_IMAGE
    RESCALED_IMAGES[${#RESCALED_IMAGES[@]}]=$RESCALE_IMAGE
    RAWIMAGES_IMAGES[${#RAWIMAGES_IMAGES[@]}]=$RAWIMAGE_IMAGE
    if [[ ! -f ${OUTPUT_IMAGE} ]];
      then
        logCmd ${ANTSPATH}ImageMath 3 $RAWIMAGE_IMAGE TruncateImageIntensity ${ANATOMICAL_IMAGES[$i]} 0.01 0.99 200
        #logCmd ${ANTSPATH}ImageMath 3 $OUTPUT_IMAGE TruncateImageIntensity ${ANATOMICAL_IMAGES[$i]} 0.01 0.99 200
#         logCmd ${ANTSPATH}N4BiasFieldCorrection -d 3 -c[20x20x20x10,0] -x $MASK_IMAGE -b [200] -s 2 -i $OUTPUT_IMAGE -o $OUTPUT_IMAGE
        #logCmd ${ANTSPATH}ImageMath 3 $RAWIMAGE_IMAGE m $MASK_IMAGE $OUTPUT_IMAGE
        logCmd ${ANTSPATH}TotalVariationDenoisingImage $RAWIMAGE_IMAGE $OUTPUT_IMAGE .1 10 4
#        logCmd ${ANTSPATH}ImageMath 3 $OUTPUT_IMAGE HistogramMatch $OUTPUT_IMAGE ${SYMMETRIC_TEMPLATES[$i]}  200 12
      fi

    # create png
    PNG_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}.png
    if [[ ! -f ${PNG_IMAGE} ]];
      then
        logCmd c3d ${ANATOMICAL_IMAGES[$i]} -slice z $CENTROIDSLICE -clip 0 200 -pad 1x0x0vox 0x0x0vox 200 -color-map grey -type uchar -omc $PNG_IMAGE
      fi

    # create png
    PNG_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_DENOISE.png
    if [[ ! -f ${PNG_IMAGE} ]];
      then
        logCmd c3d $OUTPUT_IMAGE  -slice z $CENTROIDSLICE -clip 0 200 -pad 1x0x0vox 0x0x0vox 200 -color-map grey -type uchar -omc $PNG_IMAGE
      fi

    # compute image gradient
    if [[ ! -f ${GRADIENT_IMAGE} ]];
      then
        logCmd ${ANTSPATH}GradientMagnitudeImageFilter $OUTPUT_IMAGE $GRADIENT_IMAGE
      fi
    # create png
    PNG_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_GRADIENT.png
    if [[ ! -f ${PNG_IMAGE} ]];
      then
        logCmd c3d $GRADIENT_IMAGE -slice z $CENTROIDSLICE -clip 0 25 -pad 1x0x0vox 0x0x0vox 25 -color-map jet -type uchar -omc $PNG_IMAGE
      fi

    if [[ ! -f ${RESCALE_IMAGE} ]];
      then
        logCmd ${ANTSPATH}ImageMath 3 $RESCALE_IMAGE RescaleImage $OUTPUT_IMAGE 0 1
      fi
  done

# expecting tri-phasic + precontrast
if [[ ${#ANATOMICAL_IMAGES[@]} -eq 4 ]];
  then
      echo "Expecting Images: Pre Art Ven Del "
      # multicomponent gmm imagg
      VECTOR_ATROPOS_IMAGE=${OUTPUT_PREFIX}ALL_ATROPOS_GMM.${OUTPUT_SUFFIX}
      VECTOR_ATROPOS_IMAGE_POSTERIORS=${OUTPUT_PREFIX}ALL_ATROPOS_GMM_POSTERIORS%d.${OUTPUT_SUFFIX}
      if [[ ! -f ${VECTOR_ATROPOS_IMAGE} ]];
         then
              SEG_BASE="${ANTSPATH}/Atropos -d ${DIMENSION} -a ${NORMALIZED_IMAGES[0]} -a ${NORMALIZED_IMAGES[1]} -a ${NORMALIZED_IMAGES[2]} -a ${NORMALIZED_IMAGES[3]} "
              SEG_0="-i KMeans[${NUMBER_OF_LABELS},${CLUSTER_CENTERS[i]}] -p Socrates[1] -x $MASK_IMAGE"
              if [[ ${#CLUSTER_CENTERS[@]} -eq 0 ]];
                then
                  SEG_0="-i KMeans[${NUMBER_OF_LABELS}] -p Socrates[1] -x $MASK_IMAGE"
                fi
              SEG_1="-c [3,0] -k Gaussian -m [0.1,1x1]"
              if [[ $DIMENSION -eq 3 ]]
                then
                  SEG_1="-c [3,0] -k Gaussian -m [0.1,1x1x1]"
                fi
              SEG_2="-o [${VECTOR_ATROPOS_IMAGE},${VECTOR_ATROPOS_IMAGE_POSTERIORS}]"
      
              logCmd $SEG_BASE $SEG_0 $SEG_1 $SEG_2
         fi
      PNG_IMAGE=${OUTPUT_PREFIX}ALL_ATROPOS_GMM.png
      if [[ ! -f ${PNG_IMAGE} ]];
        then
          logCmd c3d ${VECTOR_ATROPOS_IMAGE}   -slice z $CENTROIDSLICE -dup -oli ${LABELFILE} 1.0   -type uchar -omc $PNG_IMAGE
        fi
      
      # combine all images into multi-component image for PK Modeling
      VECTOR_IMAGE=${OUTPUT_PREFIX}_DENOISE.${OUTPUT_SUFFIX}
      if [[ ! -f ${VECTOR_IMAGE} ]];
         then
           logCmd c3d ${NORMALIZED_IMAGES[0]} ${NORMALIZED_IMAGES[1]} ${NORMALIZED_IMAGES[2]} ${NORMALIZED_IMAGES[3]} -omc $VECTOR_IMAGE
         fi
      # run PK Model
      MODEL_IMAGE=${OUTPUT_PREFIX}DENOISE
      if [[ ! -f ${MODEL_IMAGE}ArtDelDeriv.nii.gz ]];
         then
           logCmd ${ANTSPATH}FitCTContrastTransportModel $VECTOR_IMAGE  $MASK_IMAGE $MODEL_IMAGE 
         fi
      # create png
      PNG_IMAGE=${MODEL_IMAGE}ArtDelDeriv.png
      if [[ ! -f ${PNG_IMAGE} ]];
        then
           logCmd c3d ${MODEL_IMAGE}PreArtDeriv.nii.gz    -slice z $CENTROIDSLICE -clip -inf inf  -color-map grey -type uchar -omc ${MODEL_IMAGE}PreArtDeriv.png
           logCmd c3d ${MODEL_IMAGE}ArtVenDeriv.nii.gz    -slice z $CENTROIDSLICE -clip -inf inf  -color-map grey -type uchar -omc ${MODEL_IMAGE}ArtVenDeriv.png
           logCmd c3d ${MODEL_IMAGE}VenDelDeriv.nii.gz    -slice z $CENTROIDSLICE -clip -inf inf  -color-map grey -type uchar -omc ${MODEL_IMAGE}VenDelDeriv.png
           logCmd c3d ${MODEL_IMAGE}AreaUnderCurve.nii.gz -slice z $CENTROIDSLICE -clip -inf inf  -color-map grey -type uchar -omc ${MODEL_IMAGE}AreaUnderCurve.png
           logCmd c3d ${MODEL_IMAGE}ArtDelDeriv.nii.gz    -slice z $CENTROIDSLICE -clip -inf inf  -color-map grey -type uchar -omc $PNG_IMAGE
        fi
  fi


################################################################################
#
# Calculate statistics images
#
################################################################################

for (( i = 0; i < ${#NORMALIZED_IMAGES[@]}; i++ ))
  do

    for (( j = 0; j < ${#RADII[@]}; j++ ))
      do
        # mean image
        OUTPUT_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_MEAN_RADIUS_${RADII[$j]}.${OUTPUT_SUFFIX}
        if [[ ! -f ${OUTPUT_IMAGE} ]];
          then
            logCmd ${ANTSPATH}/ImageMath ${DIMENSION} $OUTPUT_IMAGE NeighborhoodStats ${NORMALIZED_IMAGES[$i]} 0 ${RADII[$j]}
          fi
        # create png
        PNG_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_MEAN_RADIUS_${RADII[$j]}.png
        if [[ ! -f ${PNG_IMAGE} ]];
          then
             logCmd c3d $OUTPUT_IMAGE -slice z $CENTROIDSLICE -clip 0 200 -pad 1x0x0vox 0x0x0vox 200  -color-map grey -type uchar -omc $PNG_IMAGE
          fi


        # standard deviation image
        OUTPUT_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_SIGMA_RADIUS_${RADII[$j]}.${OUTPUT_SUFFIX}
        if [[ ! -f ${OUTPUT_IMAGE} ]];
          then
            logCmd ${ANTSPATH}/ImageMath ${DIMENSION} $OUTPUT_IMAGE NeighborhoodStats ${NORMALIZED_IMAGES[$i]} 4 ${RADII[$j]}
          fi
        PNG_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_SIGMA_RADIUS_${RADII[$j]}.png
        if [[ ! -f ${PNG_IMAGE} ]];
          then
             logCmd c3d $OUTPUT_IMAGE -slice z $CENTROIDSLICE -clip 0 25 -pad 1x0x0vox 0x0x0vox 25  -color-map grey -type uchar -omc $PNG_IMAGE
          fi

        # skewness image
        OUTPUT_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_SKEWNESS_RADIUS_${RADII[$j]}.${OUTPUT_SUFFIX}
        if [[ ! -f ${OUTPUT_IMAGE} ]];
          then
            logCmd ${ANTSPATH}/ImageMath ${DIMENSION} $OUTPUT_IMAGE NeighborhoodStats ${RAWIMAGES_IMAGES[$i]} 5 ${RADII[$j]}
          fi
        # create png
        PNG_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_SKEWNESS_RADIUS_${RADII[$j]}.png
        if [[ ! -f ${PNG_IMAGE} ]];
          then
                logCmd c3d $OUTPUT_IMAGE -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc $PNG_IMAGE
          fi

        # entropy image
        OUTPUT_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ENTROPY_RADIUS_${RADII[$j]}.${OUTPUT_SUFFIX}
        if [[ ! -f ${OUTPUT_IMAGE} ]];
          then
            logCmd ${ANTSPATH}/ImageMath ${DIMENSION} $OUTPUT_IMAGE NeighborhoodStats ${RAWIMAGES_IMAGES[$i]} 7 ${RADII[$j]}
          fi
        # create png
        PNG_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ENTROPY_RADIUS_${RADII[$j]}.png
        if [[ ! -f ${PNG_IMAGE} ]];
          then
             logCmd c3d $OUTPUT_IMAGE -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc $PNG_IMAGE
          fi

        ## TODO - otb taking too long...
        # otb cooccurence matrix images
        OTBBASEIMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_OTB
        # otb should write haralick image last...
        if [[ (! -f ${OTBBASEIMAGE}HaralickCorrelation_${RADII[$j]}.${OUTPUT_SUFFIX}) && ${ENABLEOTB} -eq 1 ]];
          then
            logCmd ${ANTSPATH}/otbScalarImageToTexturesFilter ${RESCALED_IMAGES[$i]} $OTBBASEIMAGE 20 ${RADII[$j]}
          fi
        # create png
        PNG_IMAGE=${OTBBASEIMAGE}HaralickCorrelation_${RADII[$j]}.png
        if [[ ! -f ${PNG_IMAGE} && ${ENABLEOTB} -eq 1 ]];
          then
            logCmd c3d ${OTBBASEIMAGE}ClusterProminence_${RADII[$j]}.nii.gz       -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc ${OTBBASEIMAGE}ClusterProminence_${RADII[$j]}.png
            logCmd c3d ${OTBBASEIMAGE}Energy_${RADII[$j]}.nii.gz                  -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc ${OTBBASEIMAGE}Energy_${RADII[$j]}.png
            logCmd c3d ${OTBBASEIMAGE}Inertia_${RADII[$j]}.nii.gz                 -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc ${OTBBASEIMAGE}Inertia_${RADII[$j]}.png
            logCmd c3d ${OTBBASEIMAGE}ClusterShade_${RADII[$j]}.nii.gz            -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc ${OTBBASEIMAGE}ClusterShade_${RADII[$j]}.png
            logCmd c3d ${OTBBASEIMAGE}Entropy_${RADII[$j]}.nii.gz                 -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc ${OTBBASEIMAGE}Entropy_${RADII[$j]}.png
            logCmd c3d ${OTBBASEIMAGE}InverseDifferenceMoment_${RADII[$j]}.nii.gz -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc ${OTBBASEIMAGE}InverseDifferenceMoment_${RADII[$j]}.png
            logCmd c3d ${OTBBASEIMAGE}Correlation_${RADII[$j]}.nii.gz             -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc ${OTBBASEIMAGE}Correlation_${RADII[$j]}.png
            logCmd c3d ${OTBBASEIMAGE}HaralickCorrelation_${RADII[$j]}.nii.gz     -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc ${OTBBASEIMAGE}HaralickCorrelation_${RADII[$j]}.png
          fi

      done
  done

################################################################################
#
# Calculate normalized distance image (inside mask values between [0,1])
#
################################################################################

OUTPUT_DISTANCE_IMAGE=${OUTPUT_PREFIX}NORMALIZED_DISTANCE.${OUTPUT_SUFFIX}
if [[ ! -f ${OUTPUT_DISTANCE_IMAGE} ]];
  then
    logCmd ${ANTSPATH}ImageMath ${DIMENSION} $OUTPUT_DISTANCE_IMAGE MaurerDistance $MASK_IMAGE 0
    logCmd ${ANTSPATH}ImageMath ${DIMENSION} $OUTPUT_DISTANCE_IMAGE m $MASK_IMAGE $OUTPUT_DISTANCE_IMAGE
    logCmd ${ANTSPATH}ImageMath ${DIMENSION} $OUTPUT_DISTANCE_IMAGE Normalize $OUTPUT_DISTANCE_IMAGE
  fi
# create png
PNG_IMAGE=${OUTPUT_PREFIX}NORMALIZED_DISTANCE.png
if [[ ! -f ${PNG_IMAGE} ]];
  then
    logCmd c3d $OUTPUT_DISTANCE_IMAGE -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc $PNG_IMAGE
  fi


## ################################################################################
## #
## # Compute transform between anatomical image and symmetric template
## #
## ################################################################################
## 
## OUTPUT_REGISTRATION_PREFIX=${OUTPUT_PREFIX}ANTs_REGISTRATION
## 
## MESH_SIZE='10x10'
## MESH_SIZE2='0x0'
## if [[ $DIMENSION -eq 3 ]];
##   then
##     MESH_SIZE=${MESH_SIZE}x11
##     MESH_SIZE2=${MESH_SIZE2}x0
##   fi
## 
## REG_BASE="${ANTSPATH}/antsRegistration -d ${DIMENSION} -w [0.01,0.995] -o ${OUTPUT_REGISTRATION_PREFIX}"
## REG_LEV0="-r [${NORMALIZED_IMAGES[0]},${SYMMETRIC_TEMPLATES[0]},1]"
## # REG_LEV1="-t Rigid[0.1] -m MI[${NORMALIZED_IMAGES[0]},${SYMMETRIC_TEMPLATES[0]},1,32,Regular,0.25] -s 2x1x0 -f 4x2x1 -c [500x250x100,1e-8,15]"
## # REG_LEV2="-t Affine[0.1] -m MI[${NORMALIZED_IMAGES[0]},${SYMMETRIC_TEMPLATES[0]},1,32,Regular,0.25] -s 2x1x0 -f 4x2x1 -c [500x250x100,1e-8,15]"
## # REG_LEV3="-t BSplineSyN[0.1,${MESH_SIZE},${MESH_SIZE2}] -m CC[${NORMALIZED_IMAGES[0]},${SYMMETRIC_TEMPLATES[0]},1,4] -s 2x1x0 -f 4x2x1 -c [70x50x10,1e-8,15]"
## REG_LEV1="-t Rigid[0.1] -m MI[${NORMALIZED_IMAGES[0]},${SYMMETRIC_TEMPLATES[0]},1,32,Regular,0.25] -s 2x1x0 -f 8x4x2 -c [500x250x100,1e-8,15]"
## REG_LEV2="-t Affine[0.1] -m MI[${NORMALIZED_IMAGES[0]},${SYMMETRIC_TEMPLATES[0]},1,32,Regular,0.25] -s 2x1x0 -f 8x4x2 -c [500x250x100,1e-8,15]"
## REG_LEV3="-t BSplineSyN[0.1,${MESH_SIZE},${MESH_SIZE2}] -m CC[${NORMALIZED_IMAGES[0]},${SYMMETRIC_TEMPLATES[0]},1,4] -s 2x1x0x0 -f 6x4x2x1 -c [70x50x10x0,1e-8,15]"
## 
## AFFINE=${OUTPUT_REGISTRATION_PREFIX}0GenericAffine.mat
## WARP=${OUTPUT_REGISTRATION_PREFIX}1Warp.nii.gz
## INVERSE_WARP=${OUTPUT_REGISTRATION_PREFIX}1InverseWarp.nii.gz
## 
## if [[ ! -f $WARP ]];
##   then
##     logCmd $REG_BASE $REG_LEV0 $REG_LEV1 $REG_LEV2 $REG_LEV3
##   fi
## 
## # if template mask is specified, we warp it
## OUTPUT_IMAGE=${OUTPUT_PREFIX}SYMMETRIC_TEMPLATE_MASK_WARPED.${OUTPUT_SUFFIX}
## if [[ ! -f $OUTPUT_IMAGE && -f $SYMMETRIC_TEMPLATE_MASK ]];
##   then
##     logCmd ${ANTSPATH}/antsApplyTransforms -d ${DIMENSION} -n MultiLabel -r ${NORMALIZED_IMAGES[0]} -i ${SYMMETRIC_TEMPLATE_MASK} -o ${OUTPUT_IMAGE} -t $WARP -t $AFFINE
## 
##     # if the template mask is specified then we also calculate the distance
##     # image as a quasi-coordinate system and see how it is deformed with the
##     # warp.  This warped coordinate system is a better estimate of the
##     # mask of the individual subject.
## 
##     OUTPUT_IMAGE=${OUTPUT_PREFIX}SYMMETRIC_TEMPLATE_MASK_WARPED_DISTANCE.${OUTPUT_SUFFIX}
##     logCmd ${ANTSPATH}ImageMath ${DIMENSION} $OUTPUT_IMAGE MaurerDistance $SYMMETRIC_TEMPLATE_MASK 0
##     logCmd ${ANTSPATH}ImageMath ${DIMENSION} $OUTPUT_IMAGE m $SYMMETRIC_TEMPLATE_MASK $OUTPUT_IMAGE
##     logCmd ${ANTSPATH}ImageMath ${DIMENSION} $OUTPUT_IMAGE Normalize $OUTPUT_IMAGE
##     logCmd ${ANTSPATH}/antsApplyTransforms -d ${DIMENSION} -n Linear -r ${NORMALIZED_IMAGES[0]} -i ${OUTPUT_IMAGE} -o ${OUTPUT_IMAGE} -t $WARP -t $AFFINE
##   fi
## 
## # log jacobian image
## OUTPUT_IMAGE=${OUTPUT_PREFIX}LOG_JACOBIAN.${OUTPUT_SUFFIX}
## if [[ ! -f ${OUTPUT_IMAGE} ]];
##   then
##     logCmd ${ANTSPATH}/ANTSJacobian ${DIMENSION} $WARP ${OUTPUT_REGISTRATION_PREFIX} 1
##     logCmd mv ${OUTPUT_REGISTRATION_PREFIX}logjacobian.nii.gz $OUTPUT_IMAGE
##   fi
## 
## for (( i = 0; i < ${#NORMALIZED_IMAGES[@]}; i++ ))
##   do
## 
##     # symmetric template difference images
##     OUTPUT_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_SYMMETRIC_TEMPLATE_DIFFERENCE.${OUTPUT_SUFFIX}
##     if [[ ! -f ${OUTPUT_IMAGE} ]];
##       then
##         logCmd ${ANTSPATH}/antsApplyTransforms -d ${DIMENSION} -n BSpline -r ${NORMALIZED_IMAGES[$i]} -i ${SYMMETRIC_TEMPLATES[$i]} -o ${OUTPUT_IMAGE} -t $WARP -t $AFFINE
##         logCmd ${ANTSPATH}/ImageMath ${DIMENSION} $OUTPUT_IMAGE - ${NORMALIZED_IMAGES[$i]} $OUTPUT_IMAGE
##         logCmd ${ANTSPATH}/SmoothImage ${DIMENSION} $OUTPUT_IMAGE $SMOOTHING_SIGMA $OUTPUT_IMAGE 1
##       fi
## 
##     # contralateral difference images
##     OUTPUT_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_CONTRALATERAL_DIFFERENCE.${OUTPUT_SUFFIX}
##     if [[ ! -f ${OUTPUT_IMAGE} ]];
##       then
##         logCmd ${ANTSPATH}/antsApplyTransforms -d ${DIMENSION} -n BSpline -i ${NORMALIZED_IMAGES[$i]} -r ${SYMMETRIC_TEMPLATES[$i]} -o ${OUTPUT_IMAGE} -t [$AFFINE,1] -t $INVERSE_WARP
##         if [[ $DIMENSION -eq 3 ]];
##           then
##             logCmd ${ANTSPATH}/PermuteFlipImageOrientationAxes ${DIMENSION} $OUTPUT_IMAGE $OUTPUT_IMAGE 0 1 2 1 0 0
##           else
##             logCmd ${ANTSPATH}/PermuteFlipImageOrientationAxes ${DIMENSION} $OUTPUT_IMAGE $OUTPUT_IMAGE 0 1 1 0
##           fi
##         logCmd ${ANTSPATH}/CopyImageHeaderInformation ${DIMENSION} ${SYMMETRIC_TEMPLATES[$i]} $OUTPUT_IMAGE 1 1 1
##         logCmd ${ANTSPATH}/antsApplyTransforms -d ${DIMENSION} -n BSpline -i $OUTPUT_IMAGE -r ${NORMALIZED_IMAGES[$i]} -o ${OUTPUT_IMAGE} -t $WARP -t $AFFINE
##         logCmd ${ANTSPATH}/ImageMath ${DIMENSION} $OUTPUT_IMAGE - ${NORMALIZED_IMAGES[$i]} $OUTPUT_IMAGE
##         logCmd ${ANTSPATH}/SmoothImage 3 $OUTPUT_IMAGE $SMOOTHING_SIGMA $OUTPUT_IMAGE 1
##       fi
## 
##   done
## 
################################################################################
#
# Construct GMM or MAP-MRF probability images for each anatomical image
# Also create geometric feature images for each atropos labeled output
#
################################################################################

for (( i = 0; i < ${#NORMALIZED_IMAGES[@]}; i++ ))
  do

    if [[ -n ${SEGMENTATION_PRIOR} ]];
      then

        OUTPUT_ATROPOS_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_MAP_MRF.${OUTPUT_SUFFIX}
        OUTPUT_ATROPOS_IMAGE_POSTERIORS=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_MAP_MRF_POSTERIORS%d.${OUTPUT_SUFFIX}

        if [[ ! -f ${OUTPUT_ATROPOS_IMAGE} ]];
          then

            COMMAND_LINE_LABELS=''
            for (( j = 1; j <= ${NUMBER_OF_LABELS}; j++ ))
              do
                COMMAND_LINE_LABELS="${COMMAND_LINE_LABELS} -y $j"
              done

            ATROPOSN4CMD="bash ${ANTSPATH}/antsAtroposN4.sh \
              -d ${DIMENSION} \
              -b Socrates[0] \
              -a ${NORMALIZED_IMAGES[$i]} \
              -x ${MASK_IMAGE} \
              -m 1 \
              -n 5 \
              -c ${NUMBER_OF_LABELS} \
              ${COMMAND_LINE_LABELS} \
              -p ${SEGMENTATION_PRIOR} \
              -w 0.0 \
              -o ${OUTPUT_PREFIX}${IMAGE_NAMES[$i]} \
              -k 0 \
              -t [10x10x10x10,0] \
              -s ${OUTPUT_SUFFIX}"
            echo "bias correction needed ? " $ATROPOSN4CMD
            # $ATROPOSN4CMD

            # f=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}Segmentation.${OUTPUT_SUFFIX}
            # newfile=${f/Segmentation/_ATROPOS_MAP_MRF};
            # logCmd mv $f $newfile

            # for f in `ls ${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}SegmentationPosteriors*.nii.gz`;
            #   do
            #     newfile=${f/SegmentationPosteriors/_ATROPOS_MAP_MRF_POSTERIORS};
            #     logCmd mv $f $newfile
            #   done
            echo "Run ATROPOS only..." 
            SEG_BASE="${ANTSPATH}/Atropos -d ${DIMENSION} -a ${NORMALIZED_IMAGES[$i]}"
            SEG_0="-i PriorProbabilityImages[${NUMBER_OF_LABELS},${SEGMENTATION_PRIOR},0.0] -p Socrates[0] -x $MASK_IMAGE"
            SEG_1="-c [5,0] -k Gaussian -m [0.1,1x1x1]"
            SEG_2="-o [${OUTPUT_ATROPOS_IMAGE},${OUTPUT_ATROPOS_IMAGE_POSTERIORS}]"

            logCmd $SEG_BASE $SEG_0 $SEG_1 $SEG_2
          fi

        OUTPUT_ATROPOS_FEATURES_PREFIX=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_MAP_MRF_
        OUTPUT_ATROPOS_DISTANCE_IMAGE=${OUTPUT_ATROPOS_FEATURES_PREFIX}LABEL${CORE_LABEL}_DISTANCE.${OUTPUT_SUFFIX}

        if [[ ! -f ${OUTPUT_ATROPOS_FEATURES_PREFIX}ECCENTRICITY.nii.gz ]];
          then
            logCmd ${ANTSPATH}/GetConnectedComponentsFeatureImages ${DIMENSION} ${OUTPUT_ATROPOS_IMAGE} ${OUTPUT_ATROPOS_FEATURES_PREFIX}
            logCmd ${ANTSPATH}/ThresholdImage ${DIMENSION} ${OUTPUT_ATROPOS_IMAGE} ${OUTPUT_ATROPOS_DISTANCE_IMAGE} ${CORE_LABEL} ${CORE_LABEL} 1 0
            logCmd ${ANTSPATH}/ImageMath ${DIMENSION} ${OUTPUT_ATROPOS_DISTANCE_IMAGE} MaurerDistance ${OUTPUT_ATROPOS_DISTANCE_IMAGE} 1
          fi

      else

        OUTPUT_ATROPOS_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_GMM.${OUTPUT_SUFFIX}
        OUTPUT_ATROPOS_IMAGE_POSTERIORS=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_GMM_POSTERIORS%d.${OUTPUT_SUFFIX}

        SEG_BASE="${ANTSPATH}/Atropos -d ${DIMENSION} -a ${NORMALIZED_IMAGES[$i]}"
        SEG_0="-i KMeans[${NUMBER_OF_LABELS},${CLUSTER_CENTERS[i]}] -p Socrates[1] -x $MASK_IMAGE"
        if [[ ${#CLUSTER_CENTERS[@]} -eq 0 ]];
          then
            SEG_0="-i KMeans[${NUMBER_OF_LABELS}] -p Socrates[1] -x $MASK_IMAGE"
          fi
        SEG_1="-c [0,0] -k Gaussian -m [0.1,1x1]"
        if [[ $DIMENSION -eq 3 ]]
          then
            SEG_1="-c [0,0] -k Gaussian -m [0.1,1x1x1]"
          fi
        SEG_2="-o [${OUTPUT_ATROPOS_IMAGE},${OUTPUT_ATROPOS_IMAGE_POSTERIORS}]"

        if [[ ! -f ${OUTPUT_ATROPOS_IMAGE} ]];
          then
            logCmd echo "-c [0,0] impose no spatial regularization ??"
            logCmd $SEG_BASE $SEG_0 $SEG_1 $SEG_2
          fi
        #echo $SEG_BASE $SEG_0 $SEG_1 $SEG_2

        OUTPUT_ATROPOS_FEATURES_PREFIX=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_GMM_
        OUTPUT_ATROPOS_DISTANCE_IMAGE=${OUTPUT_ATROPOS_FEATURES_PREFIX}LABEL${CORE_LABEL}_DISTANCE.${OUTPUT_SUFFIX}
        OUTPUT_ATROPOS_CANNY_IMAGE=${OUTPUT_ATROPOS_FEATURES_PREFIX}CANNY.${OUTPUT_SUFFIX}

        #echo ${ANTSPATH}/GetConnectedComponentsFeatureImages ${DIMENSION} ${OUTPUT_ATROPOS_IMAGE} ${OUTPUT_ATROPOS_FEATURES_PREFIX}
        #echo ${ANTSPATH}CannyEdgeDetectionImageFilter ${OUTPUT_ATROPOS_IMAGE} ${OUTPUT_ATROPOS_CANNY_IMAGE}
        #echo ${ANTSPATH}/ImageMath ${DIMENSION} ${OUTPUT_ATROPOS_CANNY_IMAGE} MaurerDistance ${OUTPUT_ATROPOS_CANNY_IMAGE} 1
        if [[ ! -f ${OUTPUT_ATROPOS_FEATURES_PREFIX}ECCENTRICITY.nii.gz ]];
          then
            logCmd ${ANTSPATH}CannyEdgeDetectionImageFilter ${OUTPUT_ATROPOS_IMAGE} ${OUTPUT_ATROPOS_CANNY_IMAGE}
            logCmd ${ANTSPATH}/GetConnectedComponentsFeatureImages ${DIMENSION} ${OUTPUT_ATROPOS_IMAGE} ${OUTPUT_ATROPOS_FEATURES_PREFIX}
            logCmd echo ${ANTSPATH}/ThresholdImage ${DIMENSION} ${OUTPUT_ATROPOS_IMAGE} ${OUTPUT_ATROPOS_DISTANCE_IMAGE} ${CORE_LABEL} ${CORE_LABEL} 1 0
            logCmd c3d ${OUTPUT_ATROPOS_IMAGE} -thresh ${CORE_LABEL} ${CORE_LABEL} 1 0 ${OUTPUT_DISTANCE_IMAGE} -thresh .17 inf 1 0 -multiply -comp -thresh 1 1 1 0  -o ${OUTPUT_ATROPOS_DISTANCE_IMAGE} 
            #logCmd c3d ${OUTPUT_ATROPOS_IMAGE} -thresh ${CORE_LABEL} ${CORE_LABEL} 1 0 -comp -thresh 1 1 1 0 -holefill 1 0 ${OUTPUT_ATROPOS_DISTANCE_IMAGE} 
            logCmd ${ANTSPATH}/ImageMath ${DIMENSION} ${OUTPUT_ATROPOS_DISTANCE_IMAGE} MaurerDistance ${OUTPUT_ATROPOS_DISTANCE_IMAGE} 1
            logCmd ${ANTSPATH}/ImageMath ${DIMENSION} ${OUTPUT_ATROPOS_CANNY_IMAGE} MaurerDistance ${OUTPUT_ATROPOS_CANNY_IMAGE} 1
          fi
        # create png
        PNG_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_GMM_CANNY.png
        if [[ ! -f ${PNG_IMAGE} ]];
          then
            # display segmentation
            logCmd c3d ${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_GMM.${OUTPUT_SUFFIX}                              -slice z $CENTROIDSLICE -dup -oli ${LABELFILE} 1.0   -type uchar -omc ${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_GMM.png
            # display images
            logCmd c3d ${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_GMM_VOLUME_TO_SURFACE_AREA_RATIO.${OUTPUT_SUFFIX} -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc ${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_GMM_VOLUME_TO_SURFACE_AREA_RATIO.png
            logCmd c3d ${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_GMM_PHYSICAL_VOLUME.${OUTPUT_SUFFIX}              -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc ${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_GMM_PHYSICAL_VOLUME.png
            logCmd c3d ${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_GMM_ECCENTRICITY.${OUTPUT_SUFFIX}                 -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc ${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_GMM_ECCENTRICITY.png
            logCmd c3d ${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_GMM_ELONGATION.${OUTPUT_SUFFIX}                   -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc ${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_GMM_ELONGATION.png
            logCmd c3d ${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_GMM_CANNY.${OUTPUT_SUFFIX}                        -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc $PNG_IMAGE
          fi

        for (( j = 1; j <= ${NUMBER_OF_LABELS}; j++ ))
          do
            COMMAND_LINE_LABELS="${COMMAND_LINE_LABELS} -y $j"
            # create png
            PNG_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_GMM_POSTERIORS$j.png
            if [[ ! -f ${PNG_IMAGE} ]];
              then
                 logCmd c3d ${OUTPUT_PREFIX}${IMAGE_NAMES[$i]}_ATROPOS_GMM_POSTERIORS$j.${OUTPUT_SUFFIX}                 -slice z $CENTROIDSLICE -clip 2% 98%  -color-map grey -type uchar -omc $PNG_IMAGE
              fi

          done

      fi
  done

################################################################################
#
# Create difference pairs images
#
################################################################################

for (( i = 0; i < ${#DIFFERENCE_PAIRS[@]}; i++ ))
  do
    DIFFERENCE_PAIR_ARRAY=( ${DIFFERENCE_PAIRS//x/ } )

    if [[ ${#DIFFERENCE_PAIR_ARRAY[@]} -eq 2 ]];
      then
        MINUEND_INDEX=${DIFFERENCE_PAIR_ARRAY[0]}
        SUBTRAHEND_INDEX=${DIFFERENCE_PAIR_ARRAY[1]}

        OUTPUT_DIFFERENCE_PAIR_IMAGE=${OUTPUT_PREFIX}${IMAGE_NAMES[${MINUEND_INDEX}]}_${IMAGE_NAMES[${SUBTRAHEND_INDEX}]}_DIFFERENCE.nii.gz

        if [[ ! -f ${OUTPUT_DIFFERENCE_PAIR_IMAGE} ]];
          then
            logCmd ${ANTSPATH}/ImageMath ${DIMENSION} ${OUTPUT_DIFFERENCE_PAIR_IMAGE} - ${NORMALIZED_IMAGES[${MINUEND_INDEX}]} ${NORMALIZED_IMAGES[${SUBTRAHEND_INDEX}]}
          fi
      fi
  done

################################################################################
#
# End of main routine
#
################################################################################

time_end=`date +%s`
time_elapsed=$((time_end - time_start))

echo
echo "--------------------------------------------------------------------------------------"
echo " Done with creating feature images"
echo " Script executed in $time_elapsed seconds"
echo " $(( time_elapsed / 3600 ))h $(( time_elapsed %3600 / 60 ))m $(( time_elapsed % 60 ))s"
echo "--------------------------------------------------------------------------------------"