old/mitkImageStatisticsCalculator.cpp

Go to the documentation of this file.

/*===================================================================

The Medical Imaging Interaction Toolkit (MITK)

Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.

This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.

See LICENSE.txt or http://www.mitk.org for details.

===================================================================*/


#include "mitkImageStatisticsCalculator.h"
#include "mitkImageAccessByItk.h"
#include "mitkImageCast.h"
#include "mitkExtractImageFilter.h"
#include "mitkImageTimeSelector.h"
#include "mitkITKImageImport.h"

#include <mitkExtendedStatisticsImageFilter.h>
#include <mitkExtendedLabelStatisticsImageFilter.h>

#include <itkScalarImageToHistogramGenerator.h>

#include <itkChangeInformationImageFilter.h>
#include <itkExtractImageFilter.h>
#include <itkMinimumMaximumImageCalculator.h>
#include <itkStatisticsImageFilter.h>
#include <itkLabelStatisticsImageFilter.h>
#include <itkMaskImageFilter.h>
#include <itkImageFileWriter.h>
#include <itkRescaleIntensityImageFilter.h>

#include <itkMaskImageFilter.h>

#include <itkImageRegionConstIterator.h>
#include <itkImageRegionIterator.h>

#include <itkCastImageFilter.h>
#include <itkVTKImageImport.h>
#include <itkVTKImageExport.h>
#include <itkImageDuplicator.h>

#include <vtkPoints.h>
#include <vtkImageStencil.h>
#include <vtkImageImport.h>
#include <vtkImageExport.h>
#include <vtkLassoStencilSource.h>

#include <itkFFTConvolutionImageFilter.h>
#include <itkConstantBoundaryCondition.h>
#include <itkImageDuplicator.h>

#include <itkContinuousIndex.h>
#include <itkNumericTraits.h>
#include <list>

#include <exception>

#include "itkImage.h"


//#define DEBUG_HOTSPOTSEARCH

#define _USE_MATH_DEFINES
#include <math.h>

#include "vtkLassoStencilSource.h"


namespace mitk
{
  ImageStatisticsCalculator::ImageStatisticsCalculator()
    : m_MaskingMode( MASKING_MODE_NONE ),
    m_MaskingModeChanged( false ),
    m_IgnorePixelValue(0.0),
    m_DoIgnorePixelValue(false),
    m_IgnorePixelValueChanged(false),
    m_PlanarFigureAxis (0),
    m_PlanarFigureSlice (0),
    m_PlanarFigureCoordinate0 (0),
    m_PlanarFigureCoordinate1 (0),
    m_HistogramBinSize(1.0),
    m_UseDefaultBinSize(true),
    m_UseBinSizeBasedOnVOIRegion(false),
    m_HotspotRadiusInMM(6.2035049089940),   // radius of a 1cm3 sphere in mm
    m_CalculateHotspot(false),
    m_HotspotRadiusInMMChanged(false),
    m_HotspotMustBeCompletelyInsideImage(true)
  {
    m_EmptyHistogram = HistogramType::New();
    m_EmptyHistogram->SetMeasurementVectorSize(1);
    HistogramType::SizeType histogramSize(1);
    histogramSize.Fill( 256 );
    m_EmptyHistogram->Initialize( histogramSize );

    m_EmptyStatistics.Reset();
  }

  ImageStatisticsCalculator::~ImageStatisticsCalculator()
  {
  }


  void ImageStatisticsCalculator::SetUseDefaultBinSize(bool useDefault)
  {
    m_UseDefaultBinSize = useDefault;
  }




  ImageStatisticsCalculator::Statistics::Statistics(bool withHotspotStatistics)
    :m_HotspotStatistics(withHotspotStatistics ? new Statistics(false) : nullptr)
  {
    Reset();
  }

  ImageStatisticsCalculator::Statistics::Statistics(const Statistics& other)
    :m_HotspotStatistics( nullptr)
  {
    this->SetLabel( other.GetLabel() );
    this->SetN( other.GetN() );
    this->SetMin( other.GetMin() );
    this->SetMax( other.GetMax() );
    this->SetMedian( other.GetMedian() );
    this->SetMean( other.GetMean() );
    this->SetVariance( other.GetVariance() );
    this->SetKurtosis( other.GetKurtosis() );
    this->SetSkewness( other.GetSkewness() );
    this->SetUniformity( other.GetUniformity() );
    this->SetEntropy( other.GetEntropy() );
    this->SetUPP( other.GetUPP() );
    this->SetMPP( other.GetMPP() );
    this->SetSigma( other.GetSigma() );
    this->SetRMS( other.GetRMS() );
    this->SetMaxIndex( other.GetMaxIndex() );
    this->SetMinIndex( other.GetMinIndex() );
    this->SetHotspotIndex( other.GetHotspotIndex() );

    if (other.m_HotspotStatistics)
    {
      this->m_HotspotStatistics = new Statistics(false);
      *this->m_HotspotStatistics = *other.m_HotspotStatistics;
    }
  }

  bool ImageStatisticsCalculator::Statistics::HasHotspotStatistics() const
  {
    return m_HotspotStatistics != nullptr;
  }

  void ImageStatisticsCalculator::Statistics::SetHasHotspotStatistics(bool hasHotspotStatistics)
  {
    m_HasHotspotStatistics = hasHotspotStatistics;
  }


  ImageStatisticsCalculator::Statistics::~Statistics()
  {
    delete m_HotspotStatistics;
  }

  double ImageStatisticsCalculator::Statistics::GetVariance() const
  {
    return this->Variance;
  }

  void ImageStatisticsCalculator::Statistics::SetVariance( const double value )
  {
    if( this->Variance != value )
    {
      if( value < 0.0 )
      {
        this->Variance = 0.0; // if given value is negative set variance to 0.0
      }
      else
      {
        this->Variance = value;
      }
    }
  }

  double ImageStatisticsCalculator::Statistics::GetSigma() const
  {
    return this->Sigma;
  }

  void ImageStatisticsCalculator::Statistics::SetSigma( const double value )
  {
    if( this->Sigma != value )
    {
      // for some compiler the value != value works to check for NaN but not for all
      // but we can always be sure that the standard deviation is a positive value
      if( value != value || value < 0.0 )
      {
        // if standard deviation is NaN we just assume 0.0
        this->Sigma = 0.0;
      }
      else
      {
        this->Sigma = value;
      }
    }
  }

  void ImageStatisticsCalculator::Statistics::Reset(unsigned int dimension)
  {
    SetLabel(0);
    SetN( 0 );
    SetMin( 0.0 );
    SetMax( 0.0 );
    SetMedian( 0.0 );
    SetVariance( 0.0 );
    SetMean( 0.0 );
    SetSigma( 0.0 );
    SetRMS( 0.0 );

    SetSkewness( 0.0 );
    SetKurtosis( 0.0 );
    SetUniformity( 0.0 );
    SetEntropy( 0.0 );
    SetMPP( 0.0 );
    SetUPP( 0.0 );

    vnl_vector<int> zero;
    zero.set_size(dimension);
    for(unsigned int i = 0; i < dimension; ++i)
    {
      zero[i] = 0;
    }

    SetMaxIndex(zero);
    SetMinIndex(zero);
    SetHotspotIndex(zero);

    if (m_HotspotStatistics != nullptr)
    {
      m_HotspotStatistics->Reset(dimension);
    }
  }

  const ImageStatisticsCalculator::Statistics&
    ImageStatisticsCalculator::Statistics::GetHotspotStatistics() const
  {
    if (m_HotspotStatistics)
    {
      return *m_HotspotStatistics;
    }
    else
    {
      throw std::logic_error("Object has no hostspot statistics, see HasHotspotStatistics()");
    }
  }

  ImageStatisticsCalculator::Statistics&
    ImageStatisticsCalculator::Statistics::GetHotspotStatistics()
  {
    if (m_HotspotStatistics)
    {
      return *m_HotspotStatistics;
    }
    else
    {
      throw std::logic_error("Object has no hostspot statistics, see HasHotspotStatistics()");
    }
  }

  ImageStatisticsCalculator::Statistics&
    ImageStatisticsCalculator::Statistics::operator=(ImageStatisticsCalculator::Statistics const& other)
  {
    if (this == &other)
      return *this;

    this->SetLabel( other.GetLabel() );
    this->SetN( other.GetN() );
    this->SetMin( other.GetMin() );
    this->SetMax( other.GetMax() );
    this->SetMean( other.GetMean() );
    this->SetMedian( other.GetMedian() );
    this->SetVariance( other.GetVariance() );
    this->SetSigma( other.GetSigma() );
    this->SetRMS( other.GetRMS() );
    this->SetMinIndex( other.GetMinIndex() );
    this->SetMaxIndex( other.GetMaxIndex() );
    this->SetHotspotIndex( other.GetHotspotIndex() );
    this->SetSkewness( other.GetSkewness() );
    this->SetKurtosis( other.GetKurtosis() );
    this->SetUniformity( other.GetUniformity() );
    this->SetEntropy( other.GetEntropy() );
    this->SetUPP( other.GetUPP() );
    this->SetMPP( other.GetMPP() );


    delete this->m_HotspotStatistics;
    this->m_HotspotStatistics = nullptr;

    if (other.m_HotspotStatistics)
    {
      this->m_HotspotStatistics = new Statistics(false);
      *this->m_HotspotStatistics = *other.m_HotspotStatistics;
    }

    return *this;

  }

  void ImageStatisticsCalculator::SetImage( const mitk::Image *image )
  {
    if ( m_Image != image )
    {
      m_Image = image;
      this->Modified();
      unsigned int numberOfTimeSteps = image->GetTimeSteps();

      // Initialize vectors to time-size of this image
      m_ImageHistogramVector.resize( numberOfTimeSteps );
      m_MaskedImageHistogramVector.resize( numberOfTimeSteps );
      m_PlanarFigureHistogramVector.resize( numberOfTimeSteps );

      m_ImageStatisticsVector.resize( numberOfTimeSteps );
      m_MaskedImageStatisticsVector.resize( numberOfTimeSteps );
      m_PlanarFigureStatisticsVector.resize( numberOfTimeSteps );

      m_ImageStatisticsTimeStampVector.resize( numberOfTimeSteps );
      m_MaskedImageStatisticsTimeStampVector.resize( numberOfTimeSteps );
      m_PlanarFigureStatisticsTimeStampVector.resize( numberOfTimeSteps );

      m_ImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps );
      m_MaskedImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps );
      m_PlanarFigureStatisticsCalculationTriggerVector.resize( numberOfTimeSteps );

      for ( unsigned int t = 0; t < image->GetTimeSteps(); ++t )
      {
        m_ImageStatisticsTimeStampVector[t].Modified();
        m_ImageStatisticsCalculationTriggerVector[t] = true;
      }
    }
  }


  void ImageStatisticsCalculator::SetImageMask( const mitk::Image *imageMask )
  {
    if ( m_Image.IsNull() )
    {
      itkExceptionMacro( << "Image needs to be set first!" );
    }

    if ( m_ImageMask != imageMask )
    {
      m_ImageMask = imageMask;
      this->Modified();

      for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t )
      {
        m_MaskedImageStatisticsTimeStampVector[t].Modified();
        m_MaskedImageStatisticsCalculationTriggerVector[t] = true;
      }
    }
  }


  void ImageStatisticsCalculator::SetPlanarFigure( mitk::PlanarFigure *planarFigure )
  {
    if ( m_Image.IsNull() )
    {
      itkExceptionMacro( << "Image needs to be set first!" );
    }

    if ( m_PlanarFigure != planarFigure )
    {
      m_PlanarFigure = planarFigure;
      this->Modified();

      for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t )
      {
        m_PlanarFigureStatisticsTimeStampVector[t].Modified();
        m_PlanarFigureStatisticsCalculationTriggerVector[t] = true;
      }
    }
  }


  void ImageStatisticsCalculator::SetMaskingMode( unsigned int mode )
  {
    if ( m_MaskingMode != mode )
    {
      m_MaskingMode = mode;
      m_MaskingModeChanged = true;
      this->Modified();
    }
  }


  void ImageStatisticsCalculator::SetMaskingModeToNone()
  {
    if ( m_MaskingMode != MASKING_MODE_NONE )
    {
      m_MaskingMode = MASKING_MODE_NONE;
      m_MaskingModeChanged = true;
      this->Modified();
    }
  }


  void ImageStatisticsCalculator::SetMaskingModeToImage()
  {
    if ( m_MaskingMode != MASKING_MODE_IMAGE )
    {
      m_MaskingMode = MASKING_MODE_IMAGE;
      m_MaskingModeChanged = true;
      this->Modified();
    }
  }


  void ImageStatisticsCalculator::SetMaskingModeToPlanarFigure()
  {
    if ( m_MaskingMode != MASKING_MODE_PLANARFIGURE )
    {
      m_MaskingMode = MASKING_MODE_PLANARFIGURE;
      m_MaskingModeChanged = true;
      this->Modified();
    }
  }

  void ImageStatisticsCalculator::SetIgnorePixelValue(double value)
  {
    if ( m_IgnorePixelValue != value )
    {
      m_IgnorePixelValue = value;
      if(m_DoIgnorePixelValue)
      {
        m_IgnorePixelValueChanged = true;
      }
      this->Modified();
    }
  }

  double ImageStatisticsCalculator::GetIgnorePixelValue()
  {
    return m_IgnorePixelValue;
  }

  void ImageStatisticsCalculator::SetDoIgnorePixelValue(bool value)
  {
    if ( m_DoIgnorePixelValue != value )
    {
      m_DoIgnorePixelValue = value;
      m_IgnorePixelValueChanged = true;
      this->Modified();
    }
  }

  bool ImageStatisticsCalculator::GetDoIgnorePixelValue()
  {
    return m_DoIgnorePixelValue;
  }

  void ImageStatisticsCalculator::SetHistogramBinSize(double size)
  {
    this->m_HistogramBinSize = size;
  }
  double ImageStatisticsCalculator::GetHistogramBinSize()
  {
    return this->m_HistogramBinSize;
  }

  void ImageStatisticsCalculator::SetHotspotRadiusInMM(double value)
  {
    if ( m_HotspotRadiusInMM != value )
    {
      m_HotspotRadiusInMM = value;
      if(m_CalculateHotspot)
      {
        m_HotspotRadiusInMMChanged = true;
        //MITK_INFO <<"Hotspot radius changed, new convolution required";
      }
      this->Modified();
    }
  }

  double ImageStatisticsCalculator::GetHotspotRadiusInMM()
  {
    return m_HotspotRadiusInMM;
  }

  void ImageStatisticsCalculator::SetCalculateHotspot(bool on)
  {
    if ( m_CalculateHotspot != on )
    {
      m_CalculateHotspot = on;
      m_HotspotRadiusInMMChanged = true;
      //MITK_INFO <<"Hotspot calculation changed, new convolution required";
      this->Modified();
    }
  }

  bool ImageStatisticsCalculator::IsHotspotCalculated()
  {
    return m_CalculateHotspot;
  }

  void ImageStatisticsCalculator::SetHotspotMustBeCompletlyInsideImage(bool hotspotMustBeCompletelyInsideImage, bool warn)
  {
    m_HotspotMustBeCompletelyInsideImage = hotspotMustBeCompletelyInsideImage;
    if (!m_HotspotMustBeCompletelyInsideImage && warn)
    {
      MITK_WARN << "Hotspot calculation will extrapolate pixels at image borders. Be aware of the consequences for the hotspot location.";
    }
  }

  bool ImageStatisticsCalculator::GetHotspotMustBeCompletlyInsideImage() const
  {
    return m_HotspotMustBeCompletelyInsideImage;
  }


  /* Implementation of the min max values for setting the range of the histogram */
  template < typename TPixel, unsigned int VImageDimension >
  void  ImageStatisticsCalculator::GetMinAndMaxValue( double &min,
    double &max, int &counter, double &sigma,
    const itk::Image< TPixel, VImageDimension > *InputImage,
    itk::Image< unsigned short, VImageDimension > *MaskedImage )
  {
    typedef itk::Image< unsigned short, VImageDimension > MaskImageType;
    typedef itk::Image< TPixel, VImageDimension > ImageType;

    typedef itk::ImageRegionConstIteratorWithIndex<ImageType> Imageie;
    typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> Imageie2;

    Imageie2 labelIterator2( MaskedImage, MaskedImage->GetRequestedRegion() );
    Imageie labelIterator3( InputImage, InputImage->GetRequestedRegion() );

    max = 0;
    min = 0;
    counter = 0;
    sigma = 0;
    double SumOfSquares = 0;
    double sumSquared = 0;

    double actualPielValue = 0;
    int counterOfPixelsInROI = 0;


    for( labelIterator2.GoToBegin(); !labelIterator2.IsAtEnd(); ++labelIterator2, ++labelIterator3)
    {
      if( labelIterator2.Value()== 1.0)
      {
        counter++;

        counterOfPixelsInROI++;
        actualPielValue = labelIterator3.Value();

        sumSquared = sumSquared + actualPielValue;
        SumOfSquares = SumOfSquares + std::pow(actualPielValue,2);

        if(counterOfPixelsInROI == 1)
        {
          max = actualPielValue;
          min = actualPielValue;
        }

        if(actualPielValue >= max)
        {
          max = actualPielValue;
        }
        else if(actualPielValue <= min)
        {
          min = actualPielValue;
        }

      }
    }

    if (counter > 1)
    {
      sigma = ( SumOfSquares - std::pow( sumSquared, 2) / counter ) / ( counter-1 );
    }
    else
    {
      sigma = 0;
    }

  }


  bool ImageStatisticsCalculator::ComputeStatistics( unsigned int timeStep )
  {

    if (m_Image.IsNull() )
    {
      mitkThrow() << "Image not set!";
    }

    if (!m_Image->IsInitialized())
    {
      mitkThrow() << "Image not initialized!";
    }

    if ( m_Image->GetReferenceCount() == 1 )
    {
      // Image no longer valid; we are the only ones to still hold a reference on it
      return false;
    }

    if ( timeStep >= m_Image->GetTimeSteps() )
    {
      throw std::runtime_error( "Error: invalid time step!" );
    }

    // If a mask was set but we are the only ones to still hold a reference on
    // it, delete it.
    if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() == 1) )
    {
      m_ImageMask = nullptr;
    }

    // Check if statistics is already up-to-date
    unsigned long imageMTime = m_ImageStatisticsTimeStampVector[timeStep].GetMTime();
    unsigned long maskedImageMTime = m_MaskedImageStatisticsTimeStampVector[timeStep].GetMTime();
    unsigned long planarFigureMTime = m_PlanarFigureStatisticsTimeStampVector[timeStep].GetMTime();

    bool imageStatisticsCalculationTrigger = m_ImageStatisticsCalculationTriggerVector[timeStep];
    bool maskedImageStatisticsCalculationTrigger = m_MaskedImageStatisticsCalculationTriggerVector[timeStep];
    bool planarFigureStatisticsCalculationTrigger = m_PlanarFigureStatisticsCalculationTriggerVector[timeStep];

    if ( !m_IgnorePixelValueChanged
      && !m_HotspotRadiusInMMChanged
      && ((m_MaskingMode != MASKING_MODE_NONE) || (imageMTime > m_Image->GetMTime() && !imageStatisticsCalculationTrigger))
      && ((m_MaskingMode != MASKING_MODE_IMAGE) || (maskedImageMTime > m_ImageMask->GetMTime() && !maskedImageStatisticsCalculationTrigger))
      && ((m_MaskingMode != MASKING_MODE_PLANARFIGURE) || (planarFigureMTime > m_PlanarFigure->GetMTime() && !planarFigureStatisticsCalculationTrigger)) )
    {
      // Statistics is up to date!
      if ( m_MaskingModeChanged )
      {
        m_MaskingModeChanged = false;
      }
      else
      {
        return false;
      }
    }

    // Reset state changed flag
    m_MaskingModeChanged = false;
    m_IgnorePixelValueChanged = false;

    // Depending on masking mode, extract and/or generate the required image
    // and mask data from the user input
    this->ExtractImageAndMask( timeStep );


    StatisticsContainer *statisticsContainer;
    HistogramContainer *histogramContainer;
    switch ( m_MaskingMode )
    {
    case MASKING_MODE_NONE:
    default:
      if(!m_DoIgnorePixelValue)
      {
        statisticsContainer = &m_ImageStatisticsVector[timeStep];
        histogramContainer = &m_ImageHistogramVector[timeStep];

        m_ImageStatisticsTimeStampVector[timeStep].Modified();
        m_ImageStatisticsCalculationTriggerVector[timeStep] = false;
      }
      else
      {
        statisticsContainer = &m_MaskedImageStatisticsVector[timeStep];
        histogramContainer = &m_MaskedImageHistogramVector[timeStep];

        m_MaskedImageStatisticsTimeStampVector[timeStep].Modified();
        m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false;
      }
      break;

    case MASKING_MODE_IMAGE:
      statisticsContainer = &m_MaskedImageStatisticsVector[timeStep];
      histogramContainer = &m_MaskedImageHistogramVector[timeStep];

      m_MaskedImageStatisticsTimeStampVector[timeStep].Modified();
      m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false;
      break;

    case MASKING_MODE_PLANARFIGURE:
      statisticsContainer = &m_PlanarFigureStatisticsVector[timeStep];
      histogramContainer = &m_PlanarFigureHistogramVector[timeStep];

      m_PlanarFigureStatisticsTimeStampVector[timeStep].Modified();
      m_PlanarFigureStatisticsCalculationTriggerVector[timeStep] = false;
      break;
    }

    // Calculate statistics and histogram(s)
    if ( m_InternalImage->GetDimension() == 3 )
    {
      if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue )
      {
        AccessFixedDimensionByItk_2(
          m_InternalImage,
          InternalCalculateStatisticsUnmasked,
          3,
          statisticsContainer,
          histogramContainer );
      }
      else
      {
        AccessFixedDimensionByItk_3(
          m_InternalImage,
          InternalCalculateStatisticsMasked,
          3,
          m_InternalImageMask3D.GetPointer(),
          statisticsContainer,
          histogramContainer );
      }
    }
    else if ( m_InternalImage->GetDimension() == 2 )
    {
      if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue )
      {
        AccessFixedDimensionByItk_2(
          m_InternalImage,
          InternalCalculateStatisticsUnmasked,
          2,
          statisticsContainer,
          histogramContainer );
      }
      else
      {
        AccessFixedDimensionByItk_3(
          m_InternalImage,
          InternalCalculateStatisticsMasked,
          2,
          m_InternalImageMask2D.GetPointer(),
          statisticsContainer,
          histogramContainer );
      }
    }
    else
    {
      MITK_ERROR << "ImageStatistics: Image dimension not supported!";
    }


    // Release unused image smart pointers to free memory
    m_InternalImage = mitk::Image::ConstPointer();
    m_InternalImageMask3D = MaskImage3DType::Pointer();
    m_InternalImageMask2D = MaskImage2DType::Pointer();



    return true;
  }


  ImageStatisticsCalculator::BinFrequencyType
    ImageStatisticsCalculator::GetBinsAndFreuqencyForHistograms( unsigned int timeStep , unsigned int label ) const
  {
    const HistogramType *binsAndFrequencyToCalculate = this->GetHistogram(0);

    // ToDo: map should be created on stack not on heap
    std::map<int, double> returnedHistogramMap;

    unsigned int size = binsAndFrequencyToCalculate->Size();
    for( unsigned int bin=0; bin < size; ++bin )
    {
      double frequency = binsAndFrequencyToCalculate->GetFrequency( bin, 0 );
      //if( frequency > mitk::eps )
      {
        returnedHistogramMap.insert( std::pair<int, double>(binsAndFrequencyToCalculate->GetMeasurement( bin, 0 ), binsAndFrequencyToCalculate->GetFrequency( bin, 0 ) ) );
      }
    }

    return returnedHistogramMap;
  }

  const ImageStatisticsCalculator::HistogramType *
    ImageStatisticsCalculator::GetHistogram( unsigned int timeStep, unsigned int label ) const
  {
    if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) )
    {
      return nullptr;
    }

    switch ( m_MaskingMode )
    {
    case MASKING_MODE_NONE:
    default:
      {
        if(m_DoIgnorePixelValue)
          return m_MaskedImageHistogramVector[timeStep][label];

        return m_ImageHistogramVector[timeStep][label];
      }

    case MASKING_MODE_IMAGE:
      return m_MaskedImageHistogramVector[timeStep][label];

    case MASKING_MODE_PLANARFIGURE:
      return m_PlanarFigureHistogramVector[timeStep][label];
    }
  }

  const ImageStatisticsCalculator::HistogramContainer &
    ImageStatisticsCalculator::GetHistogramVector( unsigned int timeStep ) const
  {
    if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) )
    {
      return m_EmptyHistogramContainer;
    }

    switch ( m_MaskingMode )
    {
    case MASKING_MODE_NONE:
    default:
      {
        if(m_DoIgnorePixelValue)
          return m_MaskedImageHistogramVector[timeStep];

        return m_ImageHistogramVector[timeStep];
      }

    case MASKING_MODE_IMAGE:
      return m_MaskedImageHistogramVector[timeStep];

    case MASKING_MODE_PLANARFIGURE:
      return m_PlanarFigureHistogramVector[timeStep];
    }
  }


  const ImageStatisticsCalculator::Statistics &
    ImageStatisticsCalculator::GetStatistics( unsigned int timeStep, unsigned int label ) const
  {
    if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) )
    {
      return m_EmptyStatistics;
    }

    switch ( m_MaskingMode )
    {
    case MASKING_MODE_NONE:
    default:
      {
        if(m_DoIgnorePixelValue)
          return m_MaskedImageStatisticsVector[timeStep][label];

        return m_ImageStatisticsVector[timeStep][label];
      }
    case MASKING_MODE_IMAGE:
      return m_MaskedImageStatisticsVector[timeStep][label];

    case MASKING_MODE_PLANARFIGURE:
      return m_PlanarFigureStatisticsVector[timeStep][label];
    }
  }

  const ImageStatisticsCalculator::StatisticsContainer &
    ImageStatisticsCalculator::GetStatisticsVector( unsigned int timeStep ) const
  {
    if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) )
    {
      return m_EmptyStatisticsContainer;
    }

    switch ( m_MaskingMode )
    {
    case MASKING_MODE_NONE:
    default:
      {
        if(m_DoIgnorePixelValue)
          return m_MaskedImageStatisticsVector[timeStep];

        return m_ImageStatisticsVector[timeStep];
      }
    case MASKING_MODE_IMAGE:
      return m_MaskedImageStatisticsVector[timeStep];

    case MASKING_MODE_PLANARFIGURE:
      return m_PlanarFigureStatisticsVector[timeStep];
    }
  }



  void ImageStatisticsCalculator::ExtractImageAndMask( unsigned int timeStep )
  {
    if ( m_Image.IsNull() )
    {
      throw std::runtime_error( "Error: image empty!" );
    }

    if ( timeStep >= m_Image->GetTimeSteps() )
    {
      throw std::runtime_error( "Error: invalid time step!" );
    }

    ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New();
    imageTimeSelector->SetInput( m_Image );
    imageTimeSelector->SetTimeNr( timeStep );
    imageTimeSelector->UpdateLargestPossibleRegion();
    mitk::Image *timeSliceImage = imageTimeSelector->GetOutput();


    switch ( m_MaskingMode )
    {
    case MASKING_MODE_NONE:
      {
        m_InternalImage = timeSliceImage;
        m_InternalImageMask2D = nullptr;
        m_InternalImageMask3D = nullptr;
        if(m_DoIgnorePixelValue)
        {
          if( m_InternalImage->GetDimension() == 3 )
          {
            if(itk::ImageIOBase::USHORT != timeSliceImage->GetPixelType().GetComponentType())
              CastToItkImage( timeSliceImage, m_InternalImageMask3D );
            else
              CastToItkImage( timeSliceImage->Clone(), m_InternalImageMask3D  );
            m_InternalImageMask3D->FillBuffer(1);
          }
          if( m_InternalImage->GetDimension() == 2 )
          {
            if(itk::ImageIOBase::USHORT != timeSliceImage->GetPixelType().GetComponentType())
              CastToItkImage( timeSliceImage, m_InternalImageMask2D );
            else
              CastToItkImage( timeSliceImage->Clone(), m_InternalImageMask2D );
            m_InternalImageMask2D->FillBuffer(1);
          }
        }
        break;
      }

    case MASKING_MODE_IMAGE:
      {
        if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() > 1) )
        {
          if ( timeStep >= m_ImageMask->GetTimeSteps() )
          {
            // Use the last mask time step in case the current time step is bigger than the total
            // number of mask time steps.
            // It makes more sense setting this to the last mask time step than to 0.
            // For instance if you have a mask with 2 time steps and an image with 5:
            // If time step 0 is selected, the mask will use time step 0.
            // If time step 1 is selected, the mask will use time step 1.
            // If time step 2+ is selected, the mask will use time step 1.
            // If you have a mask with only one time step instead, this will always default to 0.
            timeStep = m_ImageMask->GetTimeSteps() - 1;
          }

          ImageTimeSelector::Pointer maskedImageTimeSelector = ImageTimeSelector::New();
          maskedImageTimeSelector->SetInput( m_ImageMask );
          maskedImageTimeSelector->SetTimeNr( timeStep );
          maskedImageTimeSelector->UpdateLargestPossibleRegion();
          mitk::Image *timeSliceMaskedImage = maskedImageTimeSelector->GetOutput();

          m_InternalImage = timeSliceImage;
          CastToItkImage( timeSliceMaskedImage, m_InternalImageMask3D );
        }
        else
        {
          throw std::runtime_error( "Error: image mask empty!" );
        }
        break;
      }
    case MASKING_MODE_PLANARFIGURE:
      {
        m_InternalImageMask2D = nullptr;

        if ( m_PlanarFigure.IsNull() )
        {
          throw std::runtime_error( "Error: planar figure empty!" );
        }
        if ( !m_PlanarFigure->IsClosed() )
        {
          throw std::runtime_error( "Masking not possible for non-closed figures" );
        }

        const BaseGeometry *imageGeometry = timeSliceImage->GetGeometry();
        if ( imageGeometry == nullptr )
        {
          throw std::runtime_error( "Image geometry invalid!" );
        }

        const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry();
        if ( planarFigurePlaneGeometry == nullptr )
        {
          throw std::runtime_error( "Planar-Figure not yet initialized!" );
        }

        const PlaneGeometry *planarFigureGeometry =
          dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry );
        if ( planarFigureGeometry == nullptr )
        {
          throw std::runtime_error( "Non-planar planar figures not supported!" );
        }

        // Find principal direction of PlanarFigure in input image
        unsigned int axis;
        if ( !this->GetPrincipalAxis( imageGeometry,
          planarFigureGeometry->GetNormal(), axis ) )
        {
          throw std::runtime_error( "Non-aligned planar figures not supported!" );
        }
        m_PlanarFigureAxis = axis;

        // Find slice number corresponding to PlanarFigure in input image
        MaskImage3DType::IndexType index;
        imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index );

        unsigned int slice = index[axis];
        m_PlanarFigureSlice = slice;


        // Extract slice with given position and direction from image
        unsigned int dimension = timeSliceImage->GetDimension();

        if (dimension != 2)
        {
          ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New();
          imageExtractor->SetInput( timeSliceImage );
          imageExtractor->SetSliceDimension( axis );
          imageExtractor->SetSliceIndex( slice );
          imageExtractor->Update();
          m_InternalImage = imageExtractor->GetOutput();
        }
        else
        {
          m_InternalImage = timeSliceImage;
        }

        // Compute mask from PlanarFigure
        AccessFixedDimensionByItk_1(
          m_InternalImage,
          InternalCalculateMaskFromPlanarFigure,
          2, axis );
      }
    }

    if(m_DoIgnorePixelValue)
    {
      if ( m_InternalImage->GetDimension() == 3 )
      {
        AccessFixedDimensionByItk_1(
          m_InternalImage,
          InternalMaskIgnoredPixels,
          3,
          m_InternalImageMask3D.GetPointer() );
      }
      else if ( m_InternalImage->GetDimension() == 2 )
      {
        AccessFixedDimensionByItk_1(
          m_InternalImage,
          InternalMaskIgnoredPixels,
          2,
          m_InternalImageMask2D.GetPointer() );
      }
    }
  }


  bool ImageStatisticsCalculator::GetPrincipalAxis(
    const BaseGeometry *geometry, Vector3D vector,
    unsigned int &axis )
  {
    vector.Normalize();
    for ( unsigned int i = 0; i < 3; ++i )
    {
      Vector3D axisVector = geometry->GetAxisVector( i );
      axisVector.Normalize();

      if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps )
      {
        axis = i;
        return true;
      }
    }

    return false;
  }


  unsigned int ImageStatisticsCalculator::calcNumberOfBins(mitk::ScalarType min, mitk::ScalarType max)
  {
    return std::ceil( ( (max - min ) / m_HistogramBinSize) );
  }


  template < typename TPixel, unsigned int VImageDimension >
  void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked(
    const itk::Image< TPixel, VImageDimension > *image,
    StatisticsContainer *statisticsContainer,
    HistogramContainer* histogramContainer )
  {
    typedef itk::Image< TPixel, VImageDimension > ImageType;
    typedef typename ImageType::IndexType IndexType;

    typedef itk::Statistics::ScalarImageToHistogramGenerator< ImageType >
      HistogramGeneratorType;

    statisticsContainer->clear();
    histogramContainer->clear();

    // Progress listening...
    typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType;
    ITKCommandType::Pointer progressListener;
    progressListener = ITKCommandType::New();
    progressListener->SetCallbackFunction( this,
      &ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate );


    // Issue 100 artificial progress events since ScalarIMageToHistogramGenerator
    // does not (yet?) support progress reporting
    this->InvokeEvent( itk::StartEvent() );
    for ( unsigned int i = 0; i < 100; ++i )
    {
      this->UnmaskedStatisticsProgressUpdate();
    }

    // Calculate statistics (separate filter)
    typedef itk::ExtendedStatisticsImageFilter< ImageType > StatisticsFilterType;
    typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New();
    statisticsFilter->SetInput( image );
    statisticsFilter->SetBinSize( 100 );
    statisticsFilter->SetCoordinateTolerance( 0.001 );
    statisticsFilter->SetDirectionTolerance( 0.001 );

    unsigned long observerTag = statisticsFilter->AddObserver( itk::ProgressEvent(), progressListener );
    try
    {
      statisticsFilter->Update();
    }
    catch (const itk::ExceptionObject& e)
    {
      mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what();
    }
    catch( const std::exception& e )
    {
      //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what();
    }

    statisticsFilter->RemoveObserver( observerTag );
    this->InvokeEvent( itk::EndEvent() );

    // Calculate minimum and maximum
    typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType;
    typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New();
    minMaxFilter->SetImage( image );
    unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener );
    minMaxFilter->Compute();
    minMaxFilter->RemoveObserver( observerTag2 );
    this->InvokeEvent( itk::EndEvent() );

    Statistics statistics;
    statistics.Reset();
    statistics.SetLabel(1);
    statistics.SetN(image->GetBufferedRegion().GetNumberOfPixels());
    statistics.SetMin(statisticsFilter->GetMinimum());
    statistics.SetMax(statisticsFilter->GetMaximum());
    statistics.SetMean(statisticsFilter->GetMean());
    statistics.SetMedian(statisticsFilter->GetMedian());
    statistics.SetVariance(statisticsFilter->GetVariance());
    statistics.SetSkewness(statisticsFilter->GetSkewness());
    statistics.SetKurtosis(statisticsFilter->GetKurtosis());
    statistics.SetUniformity( statisticsFilter->GetUniformity());
    statistics.SetEntropy( statisticsFilter->GetEntropy());
    statistics.SetUPP( statisticsFilter->GetUPP());
    statistics.SetMPP( statisticsFilter->GetMPP());
    statistics.SetSigma(statisticsFilter->GetSigma());
    statistics.SetRMS(sqrt( statistics.GetMean() * statistics.GetMean() + statistics.GetSigma() * statistics.GetSigma() ));

    statistics.GetMinIndex().set_size(image->GetImageDimension());
    statistics.GetMaxIndex().set_size(image->GetImageDimension());

    vnl_vector<int> tmpMaxIndex;
    vnl_vector<int> tmpMinIndex;

    tmpMaxIndex.set_size(image->GetImageDimension() );
    tmpMinIndex.set_size(image->GetImageDimension() );

    for (unsigned int i=0; i<statistics.GetMaxIndex().size(); i++)
    {
      tmpMaxIndex[i] = minMaxFilter->GetIndexOfMaximum()[i];
      tmpMinIndex[i] = minMaxFilter->GetIndexOfMinimum()[i];
    }

    statistics.SetMinIndex(tmpMaxIndex);
    statistics.SetMinIndex(tmpMinIndex);

    if( IsHotspotCalculated() && VImageDimension == 3 )
    {
      typedef itk::Image< unsigned short, VImageDimension > MaskImageType;
      typename MaskImageType::Pointer nullMask;
      bool isHotspotDefined(false);
      Statistics hotspotStatistics = this->CalculateHotspotStatistics(image, nullMask.GetPointer(), m_HotspotRadiusInMM, isHotspotDefined, 0 );
      if (isHotspotDefined)
      {
        statistics.SetHasHotspotStatistics(true);
        statistics.GetHotspotStatistics() = hotspotStatistics;
      }
      else
      {
        statistics.SetHasHotspotStatistics(false);
      }

      if(statistics.GetHotspotStatistics().HasHotspotStatistics() )
      {
        MITK_DEBUG << "Hotspot statistics available";
        statistics.SetHotspotIndex(hotspotStatistics.GetHotspotIndex());
      }
      else
      {
        MITK_ERROR << "No hotspot statistics available!";
      }
    }

    statisticsContainer->push_back( statistics );

    // Calculate histogram
    // calculate bin size or number of bins
    unsigned int numberOfBins = 200; // default number of bins
    if (m_UseDefaultBinSize)
    {
      m_HistogramBinSize = std::ceil( (statistics.GetMax() - statistics.GetMin() + 1)/numberOfBins );
    }
    else
    {
      numberOfBins = calcNumberOfBins(statistics.GetMin(), statistics.GetMax());
    }
    typename HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New();
    histogramGenerator->SetInput( image );
    histogramGenerator->SetMarginalScale( 100 );
    histogramGenerator->SetNumberOfBins( numberOfBins );
    histogramGenerator->SetHistogramMin( statistics.GetMin() );
    histogramGenerator->SetHistogramMax( statistics.GetMax() );
    histogramGenerator->Compute();
    histogramContainer->push_back( histogramGenerator->GetOutput() );
  }

  template < typename TPixel, unsigned int VImageDimension >
  void ImageStatisticsCalculator::InternalMaskIgnoredPixels(
    const itk::Image< TPixel, VImageDimension > *image,
    itk::Image< unsigned short, VImageDimension > *maskImage )
  {
    typedef itk::Image< TPixel, VImageDimension > ImageType;
    typedef itk::Image< unsigned short, VImageDimension > MaskImageType;

    itk::ImageRegionIterator<MaskImageType>
      itmask(maskImage, maskImage->GetLargestPossibleRegion());
    itk::ImageRegionConstIterator<ImageType>
      itimage(image, image->GetLargestPossibleRegion());

    itmask.GoToBegin();
    itimage.GoToBegin();

    while( !itmask.IsAtEnd() )
    {
      if(m_IgnorePixelValue == itimage.Get())
      {
        itmask.Set(0);
      }

      ++itmask;
      ++itimage;
    }
  }

  template < typename TPixel, unsigned int VImageDimension >
  void ImageStatisticsCalculator::InternalCalculateStatisticsMasked(
    const itk::Image< TPixel, VImageDimension > *image,
    itk::Image< unsigned short, VImageDimension > *maskImage,
    StatisticsContainer* statisticsContainer,
    HistogramContainer* histogramContainer )
  {
    typedef itk::Image< TPixel, VImageDimension > ImageType;
    typedef itk::Image< unsigned short, VImageDimension > MaskImageType;
    typedef typename ImageType::IndexType IndexType;
    typedef typename ImageType::PointType PointType;
    typedef typename ImageType::SpacingType SpacingType;
    typedef typename ImageType::Pointer ImagePointer;
    typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskImageType > LabelStatisticsFilterType;
    typedef itk::ChangeInformationImageFilter< MaskImageType > ChangeInformationFilterType;
    typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType;

    statisticsContainer->clear();
    histogramContainer->clear();

    // Make sure that mask is set
    if ( maskImage == nullptr  )
    {
      itkExceptionMacro( << "Mask image needs to be set!" );
    }

    // Make sure that spacing of mask and image are the same
    //SpacingType imageSpacing = image->GetSpacing();
    //SpacingType maskSpacing = maskImage->GetSpacing();
    //PointType zeroPoint; zeroPoint.Fill( 0.0 );
    //if ( (zeroPoint + imageSpacing).SquaredEuclideanDistanceTo( (zeroPoint + maskSpacing) ) > mitk::eps )
    //{
    //  itkExceptionMacro( << "Mask needs to have same spacing as image! (Image spacing: " << imageSpacing << "; Mask spacing: " << maskSpacing << ")" );
    //}
    // Make sure that orientation of mask and image are the same
    typedef typename ImageType::DirectionType DirectionType;
    DirectionType imageDirection = image->GetDirection();
    DirectionType maskDirection = maskImage->GetDirection();
    for( int i = 0; i < imageDirection.ColumnDimensions; ++i )
    {
      for( int j = 0; j < imageDirection.ColumnDimensions; ++j )
      {
        double differenceDirection = imageDirection[i][j] - maskDirection[i][j];
        if ( fabs( differenceDirection ) > mitk::eps )
        {
          double differenceDirection = imageDirection[i][j] - maskDirection[i][j];
          if ( fabs( differenceDirection ) > 0.001 /*mitk::eps*/ ) // TODO: temp fix (bug 17121)
          {
            itkExceptionMacro( << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")" );
          }
        }
      }
    }
    // Make sure that the voxels of mask and image are correctly "aligned", i.e., voxel boundaries are the same in both images
    PointType imageOrigin = image->GetOrigin();
    PointType maskOrigin = maskImage->GetOrigin();
    long offset[ImageType::ImageDimension];

    typedef itk::ContinuousIndex<double, VImageDimension> ContinousIndexType;
    ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex;

    image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex);
    image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex);

    for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i )
    {
      double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 );
      if ( fabs( misalignment ) > mitk::eps )
      {
        itkWarningMacro( << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << misalignment << ")" );
      }

      double indexCoordDistance = maskOriginContinousIndex[i] - imageOriginContinousIndex[i];
      offset[i] = int( indexCoordDistance + image->GetBufferedRegion().GetIndex()[i] + 0.5 );
    }

    // Adapt the origin and region (index/size) of the mask so that the origin of both are the same
    typename ChangeInformationFilterType::Pointer adaptMaskFilter;
    adaptMaskFilter = ChangeInformationFilterType::New();
    adaptMaskFilter->ChangeOriginOn();
    adaptMaskFilter->ChangeRegionOn();
    adaptMaskFilter->SetInput( maskImage );
    adaptMaskFilter->SetOutputOrigin( image->GetOrigin() );
    adaptMaskFilter->SetOutputOffset( offset );
    adaptMaskFilter->SetCoordinateTolerance( 0.001 );
    adaptMaskFilter->SetDirectionTolerance( 0.001 );


    typename MaskImageType::Pointer adaptedMaskImage;
    try
    {
      adaptMaskFilter->Update();
      adaptedMaskImage = adaptMaskFilter->GetOutput();
    }
    catch( const itk::ExceptionObject &e)
    {
      mitkThrow() << "Attempt to adapt shifted origin of the mask image failed due to ITK Exception: \n" << e.what();
    }
    catch( const std::exception& e )
    {
      //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what();
    }


    // Make sure that mask region is contained within image region
    if ( adaptedMaskImage.IsNotNull() &&
      !image->GetLargestPossibleRegion().IsInside( adaptedMaskImage->GetLargestPossibleRegion() ) )
    {
      itkWarningMacro( << "Mask region needs to be inside of image region! (Image region: "
        << image->GetLargestPossibleRegion() << "; Mask region: " << adaptedMaskImage->GetLargestPossibleRegion() << ")" );
    }


    // If mask region is smaller than image region, extract the sub-sampled region from the original image
    typename ImageType::SizeType imageSize = image->GetBufferedRegion().GetSize();
    typename ImageType::SizeType maskSize = maskImage->GetBufferedRegion().GetSize();
    bool maskSmallerImage = false;
    for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i )
    {
      if ( maskSize[i] < imageSize[i] )
      {
        maskSmallerImage = true;
      }
    }

    typename ImageType::ConstPointer adaptedImage;
    if ( maskSmallerImage )
    {
      typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New();
      extractImageFilter->SetInput( image );
      extractImageFilter->SetExtractionRegion( adaptedMaskImage->GetBufferedRegion() );
      extractImageFilter->SetCoordinateTolerance( 0.001 );
      extractImageFilter->SetDirectionTolerance( 0.001 );
      extractImageFilter->Update();
      adaptedImage = extractImageFilter->GetOutput();
    }
    else
    {
      adaptedImage = image;
    }

    // Initialize Filter
    typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType;
    typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New();
    statisticsFilter->SetInput( adaptedImage );

    try
    {
      statisticsFilter->Update();
    }
    catch( const itk::ExceptionObject& e)
    {
      mitkThrow() << "Image statistics initialization computation failed with ITK Exception: \n " << e.what();
    }
    catch( const std::exception& e )
    {
      //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what();
    }

    // Calculate bin size or number of bins
    unsigned int numberOfBins = m_HistogramBinSize; // default number of bins
    double maximum = 0.0;
    double minimum = 0.0;

    if (m_UseBinSizeBasedOnVOIRegion)
    {
      maximum = statisticsFilter->GetMaximum();
      minimum = statisticsFilter->GetMinimum();

      if (m_UseDefaultBinSize)
      {
        m_HistogramBinSize = std::ceil( static_cast<double>((statisticsFilter->GetMaximum() - statisticsFilter->GetMinimum() + 1)/numberOfBins) );
      }
      else
      {
        numberOfBins = calcNumberOfBins(statisticsFilter->GetMinimum(), statisticsFilter->GetMaximum());
      }
    }
    else
    {
      double sig = 0.0;
      int counter = 0;

      //Find the min and max values for the Roi to set the range for the histogram
      GetMinAndMaxValue( minimum, maximum, counter, sig, image, maskImage);
//      numberOfBins = maximum - minimum;
//      if(maximum - minimum <= 10)
//      {
//        numberOfBins = 100;
//      }
    }


    typename LabelStatisticsFilterType::Pointer labelStatisticsFilter = LabelStatisticsFilterType::New();
    labelStatisticsFilter->SetInput( adaptedImage );
    labelStatisticsFilter->SetLabelInput( adaptedMaskImage );
    labelStatisticsFilter->SetCoordinateTolerance( 0.001 );
    labelStatisticsFilter->SetDirectionTolerance( 0.001 );
    labelStatisticsFilter->UseHistogramsOn();
    labelStatisticsFilter->SetHistogramParameters( numberOfBins, floor(minimum), ceil(maximum) );   //statisticsFilter->GetMinimum() statisticsFilter->GetMaximum()

    // Add progress listening
    typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType;
    ITKCommandType::Pointer progressListener;
    progressListener = ITKCommandType::New();
    progressListener->SetCallbackFunction( this,
      &ImageStatisticsCalculator::MaskedStatisticsProgressUpdate );
    unsigned long observerTag = labelStatisticsFilter->AddObserver(
      itk::ProgressEvent(), progressListener );

    // Execute filter
    this->InvokeEvent( itk::StartEvent() );

    // Make sure that only the mask region is considered (otherwise, if the mask region is smaller
    // than the image region, the Update() would result in an exception).
    labelStatisticsFilter->GetOutput()->SetRequestedRegion( adaptedMaskImage->GetLargestPossibleRegion() );

    // Execute the filter
    try
    {
      labelStatisticsFilter->Update();
    }
    catch( const itk::ExceptionObject& e)
    {
      mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what();
    }
    catch( const std::exception& e )
    {
      //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what();
    }

    this->InvokeEvent( itk::EndEvent() );

    if( observerTag )
      labelStatisticsFilter->RemoveObserver( observerTag );

    // Find all relevant labels of mask (other than 0)
    std::list< int > relevantLabels = labelStatisticsFilter->GetRelevantLabels();
    unsigned int i;

    if ( labelStatisticsFilter->GetMaskingNonEmpty() )
    {
      std::list< int >::iterator it;
      for ( it = relevantLabels.begin(), i = 0;
        it != relevantLabels.end();
        ++it, ++i )
      {
        Statistics statistics; // restore previous code
        labelStatisticsFilter->GetHistogram(*it) ;
        histogramContainer->push_back( HistogramType::ConstPointer( labelStatisticsFilter->GetHistogram( (*it) ) ) );

        statistics.SetLabel (*it);
        statistics.SetN(labelStatisticsFilter->GetCount( *it ));
        statistics.SetMin(labelStatisticsFilter->GetMinimum( *it ));
        statistics.SetMax(labelStatisticsFilter->GetMaximum( *it ));
        statistics.SetMean(labelStatisticsFilter->GetMean( *it ));
        statistics.SetMedian(labelStatisticsFilter->GetMedian( *it));
        statistics.SetMedian(labelStatisticsFilter->GetMedian( *it ));
        statistics.SetVariance(labelStatisticsFilter->GetVariance( *it ));
        statistics.SetSigma(labelStatisticsFilter->GetSigma( *it ));
        statistics.SetSkewness(labelStatisticsFilter->GetSkewness( *it ));
        statistics.SetKurtosis(labelStatisticsFilter->GetKurtosis( *it ));
        statistics.SetUniformity( labelStatisticsFilter->GetUniformity( *it ));
        statistics.SetEntropy( labelStatisticsFilter->GetEntropy( *it ));
        statistics.SetUPP( labelStatisticsFilter->GetUPP( *it));
        statistics.SetMPP( labelStatisticsFilter->GetMPP( *it));
        statistics.SetRMS(sqrt( statistics.GetMean() * statistics.GetMean()
          + statistics.GetSigma() * statistics.GetSigma() ));

        // restrict image to mask area for min/max index calculation
        typedef itk::MaskImageFilter< ImageType, MaskImageType, ImageType > MaskImageFilterType;
        typename MaskImageFilterType::Pointer masker = MaskImageFilterType::New();
        bool isMinAndMaxSameValue = (statistics.GetMin() == statistics.GetMax());
        // bug 17962: following is a workaround for the case when min and max are the same, we can probably find a nicer way here
        double outsideValue = (isMinAndMaxSameValue ? (statistics.GetMax()/2) : (statistics.GetMin()+statistics.GetMax())/2);
        masker->SetOutsideValue( outsideValue );
        masker->SetInput1(adaptedImage);
        masker->SetInput2(adaptedMaskImage);
        masker->SetCoordinateTolerance( 0.001 );
        masker->SetDirectionTolerance( 0.001 );
        masker->Update();
        // get index of minimum and maximum
        typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType;
        typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New();
        minMaxFilter->SetImage( masker->GetOutput() );
        unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener );
        minMaxFilter->Compute();
        minMaxFilter->RemoveObserver( observerTag2 );
        this->InvokeEvent( itk::EndEvent() );

        typename MinMaxFilterType::IndexType tempMaxIndex = minMaxFilter->GetIndexOfMaximum();
        // bug 17962: following is a workaround for the case when min and max are the same, we can probably find a nicer way here
        typename MinMaxFilterType::IndexType tempMinIndex =
          (isMinAndMaxSameValue ? minMaxFilter->GetIndexOfMaximum() : minMaxFilter->GetIndexOfMinimum());

        // FIX BUG 14644
        //If a PlanarFigure is used for segmentation the
        //adaptedImage is a single slice (2D). Adding the
        // 3. dimension.

        vnl_vector<int> maxIndex;
        vnl_vector<int> minIndex;
        maxIndex.set_size(m_Image->GetDimension());
        minIndex.set_size(m_Image->GetDimension());

        if (m_MaskingMode == MASKING_MODE_PLANARFIGURE && m_Image->GetDimension()==3)
        {
          maxIndex[m_PlanarFigureCoordinate0] = tempMaxIndex[0];
          maxIndex[m_PlanarFigureCoordinate1] = tempMaxIndex[1];
          maxIndex[m_PlanarFigureAxis] = m_PlanarFigureSlice;

          minIndex[m_PlanarFigureCoordinate0] = tempMinIndex[0] ;
          minIndex[m_PlanarFigureCoordinate1] = tempMinIndex[1];
          minIndex[m_PlanarFigureAxis] = m_PlanarFigureSlice;
        } else
        {
          for (unsigned int i = 0; i<maxIndex.size(); i++)
          {
            maxIndex[i] = tempMaxIndex[i];
            minIndex[i] = tempMinIndex[i];
          }
        }
        // FIX END
        statistics.SetMaxIndex(maxIndex);
        statistics.SetMinIndex(minIndex);
        /*****************************************************Calculate Hotspot Statistics**********************************************/

        if(IsHotspotCalculated() && VImageDimension == 3)
        {
          bool isDefined(false);
          Statistics hotspotStatistics = CalculateHotspotStatistics(adaptedImage.GetPointer(), adaptedMaskImage.GetPointer(),GetHotspotRadiusInMM(), isDefined, *it);
          statistics.GetHotspotStatistics() = hotspotStatistics;
          if(statistics.GetHotspotStatistics().HasHotspotStatistics())
          {
            MITK_DEBUG << "Hotspot statistics available";
            statistics.SetHotspotIndex( hotspotStatistics.GetHotspotIndex() );
          }
          else
          {
            MITK_ERROR << "No hotspot statistics available!";
          }
        }
        statisticsContainer->push_back( statistics );
      }
    }
    else
    {
      histogramContainer->push_back( HistogramType::ConstPointer( m_EmptyHistogram ) );
      statisticsContainer->push_back( Statistics() );
    }
  }


  template <typename TPixel, unsigned int VImageDimension  >
  ImageStatisticsCalculator::ImageExtrema
    ImageStatisticsCalculator::CalculateExtremaWorld(
    const itk::Image<TPixel, VImageDimension> *inputImage,
    itk::Image<unsigned short, VImageDimension> *maskImage,
    double neccessaryDistanceToImageBorderInMM,
    unsigned int label)
  {
    typedef itk::Image< TPixel, VImageDimension > ImageType;
    typedef itk::Image< unsigned short, VImageDimension > MaskImageType;

    typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> MaskImageIteratorType;
    typedef itk::ImageRegionConstIteratorWithIndex<ImageType> InputImageIndexIteratorType;

    typename ImageType::SpacingType spacing = inputImage->GetSpacing();

    ImageExtrema minMax;
    minMax.Defined = false;
    minMax.MaxIndex.set_size(VImageDimension);
    minMax.MaxIndex.set_size(VImageDimension);

    typename ImageType::RegionType allowedExtremaRegion = inputImage->GetLargestPossibleRegion();

    bool keepDistanceToImageBorders( neccessaryDistanceToImageBorderInMM > 0 );
    if (keepDistanceToImageBorders)
    {
      long distanceInPixels[VImageDimension];
      for(unsigned short dimension = 0; dimension < VImageDimension; ++dimension)
      {
        // To confirm that the whole hotspot is inside the image we have to keep a specific distance to the image-borders, which is as long as
        // the radius. To get the amount of indices we divide the radius by spacing and add 0.5 because voxels are center based:
        // For example with a radius of 2.2 and a spacing of 1 two indices are enough because 2.2 / 1 + 0.5 = 2.7 => 2.
        // But with a radius of 2.7 we need 3 indices because 2.7 / 1 + 0.5 = 3.2 => 3
        distanceInPixels[dimension] = int( neccessaryDistanceToImageBorderInMM / spacing[dimension] + 0.5);
      }

      allowedExtremaRegion.ShrinkByRadius(distanceInPixels);
    }

    InputImageIndexIteratorType imageIndexIt(inputImage, allowedExtremaRegion);

    float maxValue = itk::NumericTraits<float>::min();
    float minValue = itk::NumericTraits<float>::max();

    typename ImageType::IndexType maxIndex;
    typename ImageType::IndexType minIndex;

    for(unsigned short i = 0; i < VImageDimension; ++i)
    {
      maxIndex[i] = 0;
      minIndex[i] = 0;
    }

    if (maskImage != nullptr)
    {
      MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion());
      typename ImageType::IndexType imageIndex;
      typename ImageType::PointType worldPosition;
      typename ImageType::IndexType maskIndex;

      for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt)
      {
        imageIndex = maskIndex = maskIt.GetIndex();

        if(maskIt.Get() == label)
        {
          if( allowedExtremaRegion.IsInside(imageIndex) )
          {
            imageIndexIt.SetIndex( imageIndex );
            double value = imageIndexIt.Get();
            minMax.Defined = true;

            //Calculate minimum, maximum and corresponding index-values
            if( value > maxValue )
            {
              maxIndex = imageIndexIt.GetIndex();
              maxValue = value;
            }

            if(value < minValue )
            {
              minIndex = imageIndexIt.GetIndex();
              minValue = value;
            }
          }
        }
      }
    }
    else
    {
      for(imageIndexIt.GoToBegin(); !imageIndexIt.IsAtEnd(); ++imageIndexIt)
      {
        double value = imageIndexIt.Get();
        minMax.Defined = true;

        //Calculate minimum, maximum and corresponding index-values
        if( value > maxValue )
        {
          maxIndex = imageIndexIt.GetIndex();
          maxValue = value;
        }

        if(value < minValue )
        {
          minIndex = imageIndexIt.GetIndex();
          minValue = value;
        }
      }
    }

    minMax.MaxIndex.set_size(VImageDimension);
    minMax.MinIndex.set_size(VImageDimension);

    for(unsigned int i = 0; i < minMax.MaxIndex.size(); ++i)
    {
      minMax.MaxIndex[i] = maxIndex[i];
    }

    for(unsigned int i = 0; i < minMax.MinIndex.size(); ++i)
    {
      minMax.MinIndex[i] = minIndex[i];
    }

    minMax.Max = maxValue;
    minMax.Min = minValue;

    return minMax;
  }

  template <unsigned int VImageDimension>
  itk::Size<VImageDimension>
    ImageStatisticsCalculator
    ::CalculateConvolutionKernelSize(double spacing[VImageDimension], double radiusInMM)
  {
    typedef itk::Image< float, VImageDimension > KernelImageType;
    typedef typename KernelImageType::SizeType SizeType;
    SizeType maskSize;

    for(unsigned int i = 0; i < VImageDimension; ++i)
    {
      maskSize[i] = static_cast<int>( 2 * radiusInMM / spacing[i]);

      // We always want an uneven size to have a clear center point in the convolution mask
      if(maskSize[i] % 2 == 0 )
      {
        ++maskSize[i];
      }
    }
    return maskSize;
  }

  template <unsigned int VImageDimension>
  itk::SmartPointer< itk::Image<float, VImageDimension> >
    ImageStatisticsCalculator
    ::GenerateHotspotSearchConvolutionKernel(double mmPerPixel[VImageDimension], double radiusInMM)
  {
    std::stringstream ss;
    for (unsigned int i = 0; i < VImageDimension; ++i)
    {
      ss << mmPerPixel[i];
      if (i < VImageDimension -1)
        ss << ",";
    }
    MITK_DEBUG << "Update convolution kernel for spacing (" << ss.str() << ") and radius " << radiusInMM << "mm";


    double radiusInMMSquared = radiusInMM * radiusInMM;
    typedef itk::Image< float, VImageDimension > KernelImageType;
    typename KernelImageType::Pointer convolutionKernel = KernelImageType::New();

    // Calculate size and allocate mask image
    typedef typename KernelImageType::SizeType SizeType;
    SizeType maskSize = this->CalculateConvolutionKernelSize<VImageDimension>(mmPerPixel, radiusInMM);

    Point3D convolutionMaskCenterIndex; convolutionMaskCenterIndex.Fill(0.0);
    for(unsigned int i = 0; i < VImageDimension; ++i)
    {
      convolutionMaskCenterIndex[i] = 0.5 * (double)(maskSize[i]-1);
    }

    typedef typename KernelImageType::IndexType IndexType;
    IndexType maskIndex;
    maskIndex.Fill(0);

    typedef typename KernelImageType::RegionType RegionType;
    RegionType maskRegion;
    maskRegion.SetSize(maskSize);
    maskRegion.SetIndex(maskIndex);

    convolutionKernel->SetRegions(maskRegion);
    convolutionKernel->SetSpacing(mmPerPixel);
    convolutionKernel->Allocate();

    // Fill mask image values by subsampling the image grid
    typedef itk::ImageRegionIteratorWithIndex<KernelImageType> MaskIteratorType;
    MaskIteratorType maskIt(convolutionKernel,maskRegion);

    int numberOfSubVoxelsPerDimension = 2; // per dimension!
    int numberOfSubVoxels = ::pow( static_cast<float>(numberOfSubVoxelsPerDimension), static_cast<float>(VImageDimension) );
    double subVoxelSizeInPixels = 1.0 / (double)numberOfSubVoxelsPerDimension;
    double valueOfOneSubVoxel = 1.0 / (double)numberOfSubVoxels;
    double maskValue = 0.0;
    Point3D subVoxelIndexPosition;
    double distanceSquared = 0.0;

    typedef itk::ContinuousIndex<double, VImageDimension> ContinuousIndexType;
    for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt)
    {
      ContinuousIndexType indexPoint(maskIt.GetIndex());
      Point3D voxelPosition;
      for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension)
      {
        voxelPosition[dimension] = indexPoint[dimension];
      }

      maskValue = 0.0;
      Vector3D subVoxelOffset; subVoxelOffset.Fill(0.0);
      // iterate sub-voxels by iterating all possible offsets
      for (subVoxelOffset[0] = -0.5 + subVoxelSizeInPixels / 2.0;
        subVoxelOffset[0] < +0.5;
        subVoxelOffset[0] += subVoxelSizeInPixels)
      {
        for (subVoxelOffset[1] = -0.5 + subVoxelSizeInPixels / 2.0;
          subVoxelOffset[1] < +0.5;
          subVoxelOffset[1] += subVoxelSizeInPixels)
        {
          for (subVoxelOffset[2] = -0.5 + subVoxelSizeInPixels / 2.0;
            subVoxelOffset[2] < +0.5;
            subVoxelOffset[2] += subVoxelSizeInPixels)
          {
            subVoxelIndexPosition = voxelPosition + subVoxelOffset; // this COULD be integrated into the for-loops if neccessary (add voxelPosition to initializer and end condition)
            distanceSquared =
              (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] * (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0]
            + (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] * (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1]
            + (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2] * (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2];

            if (distanceSquared <= radiusInMMSquared)
            {
              maskValue += valueOfOneSubVoxel;
            }
          }
        }
      }
      maskIt.Set( maskValue );
    }

    return convolutionKernel;
  }

  template <typename TPixel, unsigned int VImageDimension>
  itk::SmartPointer<itk::Image<TPixel, VImageDimension> >
    ImageStatisticsCalculator::GenerateConvolutionImage( const itk::Image<TPixel, VImageDimension>* inputImage )
  {
    double mmPerPixel[VImageDimension];
    for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension)
    {
      mmPerPixel[dimension] = inputImage->GetSpacing()[dimension];
    }

    // update convolution kernel
    typedef itk::Image< float, VImageDimension > KernelImageType;
    typename KernelImageType::Pointer convolutionKernel = this->GenerateHotspotSearchConvolutionKernel<VImageDimension>(mmPerPixel, m_HotspotRadiusInMM);

    // update convolution image
    typedef itk::Image< TPixel, VImageDimension > InputImageType;
    typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType;
    typedef itk::FFTConvolutionImageFilter<InputImageType,
      KernelImageType,
      ConvolutionImageType> ConvolutionFilterType;

    typename ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New();
    typedef itk::ConstantBoundaryCondition<InputImageType, InputImageType> BoundaryConditionType;
    BoundaryConditionType boundaryCondition;
    boundaryCondition.SetConstant(0.0);

    if (GetHotspotMustBeCompletlyInsideImage())
    {
      // overwrite default boundary condition
      convolutionFilter->SetBoundaryCondition(&boundaryCondition);
    }

    convolutionFilter->SetInput(inputImage);
    convolutionFilter->SetKernelImage(convolutionKernel);
    convolutionFilter->SetNormalize(true);
    MITK_DEBUG << "Update Convolution image for hotspot search";
    convolutionFilter->UpdateLargestPossibleRegion();

    typename ConvolutionImageType::Pointer convolutionImage = convolutionFilter->GetOutput();
    convolutionImage->SetSpacing( inputImage->GetSpacing() ); // only workaround because convolution filter seems to ignore spacing of input image

    m_HotspotRadiusInMMChanged = false;
    return convolutionImage;
  }

  template < typename TPixel, unsigned int VImageDimension>
  void
    ImageStatisticsCalculator
    ::FillHotspotMaskPixels( itk::Image<TPixel, VImageDimension>* maskImage,
    itk::Point<double, VImageDimension> sphereCenter,
    double sphereRadiusInMM)
  {
    typedef itk::Image< TPixel, VImageDimension > MaskImageType;
    typedef itk::ImageRegionIteratorWithIndex<MaskImageType> MaskImageIteratorType;

    MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion());

    typename MaskImageType::IndexType maskIndex;
    typename MaskImageType::PointType worldPosition;

    for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt)
    {
      maskIndex = maskIt.GetIndex();
      maskImage->TransformIndexToPhysicalPoint(maskIndex, worldPosition);
      maskIt.Set( worldPosition.EuclideanDistanceTo(sphereCenter) <= sphereRadiusInMM ? 1 : 0 );
    }
  }

  template < typename TPixel, unsigned int VImageDimension>
  ImageStatisticsCalculator::Statistics
    ImageStatisticsCalculator::CalculateHotspotStatistics(
    const itk::Image<TPixel, VImageDimension>* inputImage,
    itk::Image<unsigned short, VImageDimension>* maskImage,
    double radiusInMM,
    bool& isHotspotDefined,
    unsigned int label)
  {
    // get convolution image (updated in GenerateConvolutionImage())
    typedef itk::Image< TPixel, VImageDimension > InputImageType;
    typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType;
    typedef itk::Image< float, VImageDimension > KernelImageType;
    typedef itk::Image< unsigned short, VImageDimension > MaskImageType;

    //typename ConvolutionImageType::Pointer convolutionImage = dynamic_cast<ConvolutionImageType*>(this->GenerateConvolutionImage(inputImage));
    typename ConvolutionImageType::Pointer convolutionImage = this->GenerateConvolutionImage(inputImage);

    if (convolutionImage.IsNull())
    {
      MITK_ERROR << "Empty convolution image in CalculateHotspotStatistics(). We should never reach this state (logic error).";
      throw std::logic_error("Empty convolution image in CalculateHotspotStatistics()");
    }

    // find maximum in convolution image, given the current mask
    double requiredDistanceToBorder = m_HotspotMustBeCompletelyInsideImage ? m_HotspotRadiusInMM : -1.0;
    ImageExtrema convolutionImageInformation = CalculateExtremaWorld(convolutionImage.GetPointer(), maskImage, requiredDistanceToBorder, label);

    isHotspotDefined = convolutionImageInformation.Defined;

    if (!isHotspotDefined)
    {
      m_EmptyStatistics.Reset(VImageDimension);
      MITK_ERROR << "No origin of hotspot-sphere was calculated! Returning empty statistics";
      return m_EmptyStatistics;
    }
    else
    {
      // create a binary mask around the "hotspot" region, fill the shape of a sphere around our hotspot center
      typedef itk::ImageDuplicator< InputImageType > DuplicatorType;
      typename DuplicatorType::Pointer copyMachine = DuplicatorType::New();
      copyMachine->SetInputImage(inputImage);
      copyMachine->Update();

      typedef itk::CastImageFilter< InputImageType, MaskImageType > CastFilterType;
      typename CastFilterType::Pointer caster = CastFilterType::New();
      caster->SetInput( copyMachine->GetOutput() );
      caster->Update();
      typename MaskImageType::Pointer hotspotMaskITK = caster->GetOutput();

      typedef typename InputImageType::IndexType IndexType;
      IndexType maskCenterIndex;
      for (unsigned int d =0; d< VImageDimension;++d) maskCenterIndex[d]=convolutionImageInformation.MaxIndex[d];
      typename ConvolutionImageType::PointType maskCenter;
      inputImage->TransformIndexToPhysicalPoint(maskCenterIndex,maskCenter);

      this->FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, radiusInMM);

      // calculate statistics within the binary mask
      typedef itk::ExtendedLabelStatisticsImageFilter< InputImageType, MaskImageType> LabelStatisticsFilterType;
      typename LabelStatisticsFilterType::Pointer labelStatisticsFilter;
      labelStatisticsFilter = LabelStatisticsFilterType::New();
      labelStatisticsFilter->SetInput( inputImage );
      labelStatisticsFilter->SetLabelInput( hotspotMaskITK );
      labelStatisticsFilter->SetCoordinateTolerance( 0.001 );
      labelStatisticsFilter->SetDirectionTolerance( 0.001 );

      labelStatisticsFilter->Update();

      Statistics hotspotStatistics;
      hotspotStatistics.SetHotspotIndex(convolutionImageInformation.MaxIndex);
      hotspotStatistics.SetMean(convolutionImageInformation.Max);

      if ( labelStatisticsFilter->HasLabel( 1 ) )
      {
        hotspotStatistics.SetLabel (1);
        hotspotStatistics.SetN(labelStatisticsFilter->GetCount(1));
        hotspotStatistics.SetMin(labelStatisticsFilter->GetMinimum(1));
        hotspotStatistics.SetMax(labelStatisticsFilter->GetMaximum(1));
        hotspotStatistics.SetMedian(labelStatisticsFilter->GetMedian(1));
        hotspotStatistics.SetVariance(labelStatisticsFilter->GetVariance(1));
        hotspotStatistics.SetSigma(labelStatisticsFilter->GetSigma(1));
        hotspotStatistics.SetRMS(sqrt( hotspotStatistics.GetMean() * hotspotStatistics.GetMean()
          + hotspotStatistics.GetSigma() * hotspotStatistics.GetSigma() ));

        MITK_DEBUG << "Statistics for inside hotspot: Mean " << hotspotStatistics.GetMean()
          << ", SD " << hotspotStatistics.GetSigma()
          << ", Max " << hotspotStatistics.GetMax()
          << ", Min " << hotspotStatistics.GetMin();
      }
      else
      {
        MITK_ERROR << "Uh oh! Unable to calculate statistics for hotspot region...";
        return m_EmptyStatistics;
      }

      return hotspotStatistics;
    }
  }

  template < typename TPixel, unsigned int VImageDimension >
  void ImageStatisticsCalculator::InternalCalculateMaskFromPlanarFigure(
    const itk::Image< TPixel, VImageDimension > *image, unsigned int axis )
  {
    typedef itk::Image< TPixel, VImageDimension > ImageType;

    typedef itk::CastImageFilter< ImageType, MaskImage2DType > CastFilterType;

    // Generate mask image as new image with same header as input image and
    // initialize with 1.
    typename CastFilterType::Pointer castFilter = CastFilterType::New();
    castFilter->SetInput( image );
    castFilter->Update();
    castFilter->GetOutput()->FillBuffer( 1 );

    // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object.
    // These points are used by the vtkLassoStencilSource to create
    // a vtkImageStencil.
    const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry();
    const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 );
    const mitk::BaseGeometry *imageGeometry3D = m_Image->GetGeometry( 0 );
    // If there is a second poly line in a closed planar figure, treat it as a hole.
    PlanarFigure::PolyLineType planarFigureHolePolyline;

    if (m_PlanarFigure->GetPolyLinesSize() == 2)
      planarFigureHolePolyline = m_PlanarFigure->GetPolyLine(1);


    // Determine x- and y-dimensions depending on principal axis
    // TODO use plane geometry normal to determine that automatically, then check whether the PF is aligned with one of the three principal axis
    int i0, i1;
    switch ( axis )
    {
    case 0:
      i0 = 1;
      i1 = 2;
      break;

    case 1:
      i0 = 0;
      i1 = 2;
      break;

    case 2:
    default:
      i0 = 0;
      i1 = 1;
      break;
    }
    m_PlanarFigureCoordinate0= i0;
    m_PlanarFigureCoordinate1= i1;

    // store the polyline contour as vtkPoints object
    bool outOfBounds = false;
    vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
    typename PlanarFigure::PolyLineType::const_iterator it;
    for ( it = planarFigurePolyline.begin();
      it != planarFigurePolyline.end();
      ++it )
    {
      Point3D point3D;

      // Convert 2D point back to the local index coordinates of the selected
      // image
      // Fabian: From PlaneGeometry documentation:
      // Converts a 2D point given in mm (pt2d_mm) relative to the upper-left corner of the geometry into the corresponding world-coordinate (a 3D point in mm, pt3d_mm).
      // To convert a 2D point given in units (e.g., pixels in case of an image) into a 2D point given in mm (as required by this method), use IndexToWorld.
      planarFigurePlaneGeometry->Map( *it, point3D );

      // Polygons (partially) outside of the image bounds can not be processed
      // further due to a bug in vtkPolyDataToImageStencil
      if ( !imageGeometry3D->IsInside( point3D ) )
      {
        outOfBounds = true;
      }

      // Fabian: Why convert to index coordinates?
      imageGeometry3D->WorldToIndex( point3D, point3D );

      points->InsertNextPoint( point3D[i0], point3D[i1], 0 );
    }

    vtkSmartPointer<vtkPoints> holePoints = nullptr;

    if (!planarFigureHolePolyline.empty())
    {
      holePoints = vtkSmartPointer<vtkPoints>::New();

      Point3D point3D;
      PlanarFigure::PolyLineType::const_iterator end = planarFigureHolePolyline.end();

      for (it = planarFigureHolePolyline.begin(); it != end; ++it)
      {
        // Fabian: same as above
        planarFigurePlaneGeometry->Map(*it, point3D);
        imageGeometry3D->WorldToIndex(point3D, point3D);
        holePoints->InsertNextPoint(point3D[i0], point3D[i1], 0);
      }
    }

    // mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds
    // this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero
    double bounds[6] = {0, 0, 0, 0, 0, 0};
    points->GetBounds( bounds );
    bool extent_x = (fabs(bounds[0] - bounds[1])) < mitk::eps;
    bool extent_y = (fabs(bounds[2] - bounds[3])) < mitk::eps;
    bool extent_z = (fabs(bounds[4] - bounds[5])) < mitk::eps;

    // throw an exception if a closed planar figure is deformed, i.e. has only one non-zero extent
    if ( m_PlanarFigure->IsClosed() &&
      ((extent_x && extent_y) || (extent_x && extent_z)  || (extent_y && extent_z)))
    {
      mitkThrow() << "Figure has a zero area and cannot be used for masking.";
    }

    if ( outOfBounds )
    {
      throw std::runtime_error( "Figure at least partially outside of image bounds!" );
    }

    // create a vtkLassoStencilSource and set the points of the Polygon
    vtkSmartPointer<vtkLassoStencilSource> lassoStencil = vtkSmartPointer<vtkLassoStencilSource>::New();
    lassoStencil->SetShapeToPolygon();
    lassoStencil->SetPoints( points );

    vtkSmartPointer<vtkLassoStencilSource> holeLassoStencil = nullptr;

    if (holePoints.GetPointer() != nullptr)
    {
      holeLassoStencil = vtkSmartPointer<vtkLassoStencilSource>::New();
      holeLassoStencil->SetShapeToPolygon();
      holeLassoStencil->SetPoints(holePoints);
    }

    // Export from ITK to VTK (to use a VTK filter)
    typedef itk::VTKImageImport< MaskImage2DType > ImageImportType;
    typedef itk::VTKImageExport< MaskImage2DType > ImageExportType;

    typename ImageExportType::Pointer itkExporter = ImageExportType::New();
    itkExporter->SetInput( castFilter->GetOutput() );

    vtkSmartPointer<vtkImageImport> vtkImporter = vtkSmartPointer<vtkImageImport>::New();
    this->ConnectPipelines( itkExporter, vtkImporter );

    // Apply the generated image stencil to the input image
    vtkSmartPointer<vtkImageStencil> imageStencilFilter = vtkSmartPointer<vtkImageStencil>::New();
    imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() );
    imageStencilFilter->SetStencilConnection(lassoStencil->GetOutputPort());
    imageStencilFilter->ReverseStencilOff();
    imageStencilFilter->SetBackgroundValue( 0 );
    imageStencilFilter->Update();

    vtkSmartPointer<vtkImageStencil> holeStencilFilter = nullptr;

    if (holeLassoStencil.GetPointer() != nullptr)
    {
      holeStencilFilter = vtkSmartPointer<vtkImageStencil>::New();
      holeStencilFilter->SetInputConnection(imageStencilFilter->GetOutputPort());
      holeStencilFilter->SetStencilConnection(holeLassoStencil->GetOutputPort());
      holeStencilFilter->ReverseStencilOn();
      holeStencilFilter->SetBackgroundValue(0);
      holeStencilFilter->Update();
    }

    // Export from VTK back to ITK
    vtkSmartPointer<vtkImageExport> vtkExporter = vtkSmartPointer<vtkImageExport>::New();
    vtkExporter->SetInputConnection( holeStencilFilter.GetPointer() == nullptr
      ? imageStencilFilter->GetOutputPort()
      : holeStencilFilter->GetOutputPort());
    vtkExporter->Update();

    typename ImageImportType::Pointer itkImporter = ImageImportType::New();
    this->ConnectPipelines( vtkExporter, itkImporter );
    itkImporter->Update();

    typedef itk::ImageDuplicator< ImageImportType::OutputImageType > DuplicatorType;
    DuplicatorType::Pointer duplicator = DuplicatorType::New();
    duplicator->SetInputImage( itkImporter->GetOutput() );
    duplicator->Update();

    // Store mask
    m_InternalImageMask2D = duplicator->GetOutput();
  }


  void ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate()
  {
    // Need to throw away every second progress event to reach a final count of
    // 100 since two consecutive filters are used in this case
    static int updateCounter = 0;
    if ( updateCounter++ % 2 == 0 )
    {
      this->InvokeEvent( itk::ProgressEvent() );
    }
  }


  void ImageStatisticsCalculator::MaskedStatisticsProgressUpdate()
  {
    this->InvokeEvent( itk::ProgressEvent() );
  }

}