diff --git a/Modules/DiffusionImaging/DiffusionCore/PythonRequirements.txt b/Modules/DiffusionImaging/DiffusionCore/PythonRequirements.txt
index 4df3c9b702..034237bfe9 100644
--- a/Modules/DiffusionImaging/DiffusionCore/PythonRequirements.txt
+++ b/Modules/DiffusionImaging/DiffusionCore/PythonRequirements.txt
@@ -1,6 +1,8 @@
 numpy
 SimpleITK
 dipy
 sklearn
 https://github.com/MIC-DKFZ/batchgenerators/archive/master.zip
-https://github.com/MIC-DKFZ/TractSeg/archive/v1.8.zip
+https://github.com/MIC-DKFZ/TractSeg/archive/v1.9.zip
+torch
+torchvision
diff --git a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp
index a4246ca1e1..46db14c13f 100755
--- a/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp
+++ b/Modules/DiffusionImaging/FiberTracking/Fiberfox/itkTractsToDWIImageFilter.cpp
@@ -1,1749 +1,1750 @@
 /*===================================================================
 
 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 "itkTractsToDWIImageFilter.h"
 #include <boost/progress.hpp>
 #include <vtkSmartPointer.h>
 #include <vtkPolyData.h>
 #include <vtkCellArray.h>
 #include <vtkPoints.h>
 #include <vtkPolyLine.h>
 #include <itkImageRegionIteratorWithIndex.h>
 #include <itkResampleImageFilter.h>
 #include <itkNearestNeighborInterpolateImageFunction.h>
 #include <itkBSplineInterpolateImageFunction.h>
 #include <itkCastImageFilter.h>
 #include <itkImageFileWriter.h>
 #include <itkRescaleIntensityImageFilter.h>
 #include <itkWindowedSincInterpolateImageFunction.h>
 #include <itkResampleDwiImageFilter.h>
 #include <itkKspaceImageFilter.h>
 #include <itkDftImageFilter.h>
 #include <itkAddImageFilter.h>
 #include <itkConstantPadImageFilter.h>
 #include <itkCropImageFilter.h>
 #include <mitkAstroStickModel.h>
 #include <vtkTransform.h>
 #include <iostream>
 #include <fstream>
 #include <exception>
 #include <itkImageDuplicator.h>
 #include <itksys/SystemTools.hxx>
 #include <mitkIOUtil.h>
 #include <itkDiffusionTensor3DReconstructionImageFilter.h>
 #include <itkDiffusionTensor3D.h>
 #include <itkTractDensityImageFilter.h>
 #include <mitkLexicalCast.h>
 #include <itkTractsToVectorImageFilter.h>
 #include <itkInvertIntensityImageFilter.h>
 #include <itkShiftScaleImageFilter.h>
 #include <itkExtractImageFilter.h>
 #include <itkResampleDwiImageFilter.h>
 #include <boost/algorithm/string/replace.hpp>
 #include <omp.h>
 #include <cmath>
 #include <thread>
 
 namespace itk
 {
 
 template< class PixelType >
 TractsToDWIImageFilter< PixelType >::TractsToDWIImageFilter()
   : m_StatusText("")
   , m_UseConstantRandSeed(false)
   , m_RandGen(itk::Statistics::MersenneTwisterRandomVariateGenerator::New())
 {
   m_DoubleInterpolator = itk::LinearInterpolateImageFunction< ItkDoubleImgType, float >::New();
   m_NullDir.Fill(0);
 }
 
 template< class PixelType >
 TractsToDWIImageFilter< PixelType >::~TractsToDWIImageFilter()
 {
 
 }
 
 template< class PixelType >
 TractsToDWIImageFilter< PixelType >::DoubleDwiType::Pointer TractsToDWIImageFilter< PixelType >::
 SimulateKspaceAcquisition( std::vector< DoubleDwiType::Pointer >& compartment_images )
 {
   unsigned int numFiberCompartments = m_Parameters.m_FiberModelList.size();
   // create slice object
   ImageRegion<2> sliceRegion;
   sliceRegion.SetSize(0, m_WorkingImageRegion.GetSize()[0]);
   sliceRegion.SetSize(1, m_WorkingImageRegion.GetSize()[1]);
   Vector< double, 2 > sliceSpacing;
   sliceSpacing[0] = m_WorkingSpacing[0];
   sliceSpacing[1] = m_WorkingSpacing[1];
 
   DoubleDwiType::PixelType nullPix; nullPix.SetSize(compartment_images.at(0)->GetVectorLength()); nullPix.Fill(0.0);
   auto magnitudeDwiImage = DoubleDwiType::New();
   magnitudeDwiImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing );
   magnitudeDwiImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin );
   magnitudeDwiImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
   magnitudeDwiImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
   magnitudeDwiImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
   magnitudeDwiImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
   magnitudeDwiImage->SetVectorLength( compartment_images.at(0)->GetVectorLength() );
   magnitudeDwiImage->Allocate();
   magnitudeDwiImage->FillBuffer(nullPix);
 
   m_PhaseImage = DoubleDwiType::New();
   m_PhaseImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing );
   m_PhaseImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin );
   m_PhaseImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
   m_PhaseImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
   m_PhaseImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
   m_PhaseImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
   m_PhaseImage->SetVectorLength( compartment_images.at(0)->GetVectorLength() );
   m_PhaseImage->Allocate();
   m_PhaseImage->FillBuffer(nullPix);
 
   m_KspaceImage = DoubleDwiType::New();
   m_KspaceImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing );
   m_KspaceImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin );
   m_KspaceImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
   m_KspaceImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
   m_KspaceImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
   m_KspaceImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
   m_KspaceImage->SetVectorLength( m_Parameters.m_SignalGen.m_NumberOfCoils );
   m_KspaceImage->Allocate();
   m_KspaceImage->FillBuffer(nullPix);
 
   // calculate coil positions
   double a = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(0)*m_Parameters.m_SignalGen.m_ImageSpacing[0];
   double b = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1)*m_Parameters.m_SignalGen.m_ImageSpacing[1];
   double c = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(2)*m_Parameters.m_SignalGen.m_ImageSpacing[2];
   double diagonal = sqrt(a*a+b*b)/1000;   // image diagonal in m
 
   m_CoilPointset = mitk::PointSet::New();
   std::vector< itk::Vector<double, 3> > coilPositions;
   itk::Vector<double, 3> pos; pos.Fill(0.0); pos[1] = -diagonal/2;
   itk::Vector<double, 3> center;
   center[0] = a/2-m_Parameters.m_SignalGen.m_ImageSpacing[0]/2;
   center[1] = b/2-m_Parameters.m_SignalGen.m_ImageSpacing[2]/2;
   center[2] = c/2-m_Parameters.m_SignalGen.m_ImageSpacing[1]/2;
   for (unsigned int c=0; c<m_Parameters.m_SignalGen.m_NumberOfCoils; c++)
   {
     coilPositions.push_back(pos);
     m_CoilPointset->InsertPoint(c, pos*1000 + m_Parameters.m_SignalGen.m_ImageOrigin.GetVectorFromOrigin() + center );
 
     double rz = 360.0/m_Parameters.m_SignalGen.m_NumberOfCoils * itk::Math::pi/180;
     vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity();
     rotZ[0][0] = cos(rz);
     rotZ[1][1] = rotZ[0][0];
     rotZ[0][1] = -sin(rz);
     rotZ[1][0] = -rotZ[0][1];
 
     pos.SetVnlVector(rotZ*pos.GetVnlVector());
   }
 
   auto num_slices = compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2);
   auto num_gradient_volumes = static_cast<int>(compartment_images.at(0)->GetVectorLength());
   auto max_threads = omp_get_max_threads();
   int out_threads = Math::ceil(std::sqrt(max_threads));
   int in_threads = Math::floor(std::sqrt(max_threads));
   if (out_threads > num_gradient_volumes)
   {
     out_threads = num_gradient_volumes;
     in_threads = Math::floor(static_cast<float>(max_threads/out_threads));
   }
   PrintToLog("Parallel volumes: " + boost::lexical_cast<std::string>(out_threads), false, true, true);
   PrintToLog("Threads per slice: " + boost::lexical_cast<std::string>(in_threads), false, true, true);
 
   std::list< std::tuple<unsigned int, unsigned int, unsigned int> > spikes;
   if (m_Parameters.m_Misc.m_DoAddSpikes)
     for (unsigned int i=0; i<m_Parameters.m_SignalGen.m_Spikes; i++)
     {
       // gradient volume and slice index
       auto spike = std::tuple<unsigned int, unsigned int, unsigned int>(
             m_RandGen->GetIntegerVariate()%num_gradient_volumes,
             m_RandGen->GetIntegerVariate()%num_slices,
             m_RandGen->GetIntegerVariate()%m_Parameters.m_SignalGen.m_NumberOfCoils);
       spikes.push_back(spike);
     }
 
   PrintToLog("0%   10   20   30   40   50   60   70   80   90   100%", false, true, false);
   PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false);
   unsigned long lastTick = 0;
 
   boost::progress_display disp(static_cast<unsigned long>(num_gradient_volumes)*compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2));
 
 #pragma omp parallel for num_threads(out_threads)
   for (int g=0; g<num_gradient_volumes; g++)
   {
     if (this->GetAbortGenerateData())
       continue;
 
     std::list< std::tuple<unsigned int, unsigned int> > spikeSlice;
 #pragma omp critical
     {
 
       for (auto spike : spikes)
         if (std::get<0>(spike) == static_cast<unsigned int>(g))
           spikeSlice.push_back(std::tuple<unsigned int, unsigned int>(std::get<1>(spike), std::get<2>(spike)));
     }
 
     for (unsigned int z=0; z<num_slices; z++)
     {
       std::vector< Float2DImageType::Pointer > compartment_slices;
       std::vector< float > t2Vector;
       std::vector< float > t1Vector;
 
       for (unsigned int i=0; i<compartment_images.size(); i++)
       {
         DiffusionSignalModel<double>* signalModel;
         if (i<numFiberCompartments)
           signalModel = m_Parameters.m_FiberModelList.at(i);
         else
           signalModel = m_Parameters.m_NonFiberModelList.at(i-numFiberCompartments);
 
         auto slice = Float2DImageType::New();
         slice->SetLargestPossibleRegion( sliceRegion );
         slice->SetBufferedRegion( sliceRegion );
         slice->SetRequestedRegion( sliceRegion );
         slice->SetSpacing(sliceSpacing);
         slice->Allocate();
         slice->FillBuffer(0.0);
 
         // extract slice from channel g
         for (unsigned int y=0; y<compartment_images.at(0)->GetLargestPossibleRegion().GetSize(1); y++)
           for (unsigned int x=0; x<compartment_images.at(0)->GetLargestPossibleRegion().GetSize(0); x++)
           {
             Float2DImageType::IndexType index2D; index2D[0]=x; index2D[1]=y;
             DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z;
 
             slice->SetPixel(index2D, compartment_images.at(i)->GetPixel(index3D)[g]);
           }
 
         compartment_slices.push_back(slice);
         t2Vector.push_back(signalModel->GetT2());
         t1Vector.push_back(signalModel->GetT1());
       }
 
       if (this->GetAbortGenerateData())
         continue;
 
       for (unsigned int c=0; c<m_Parameters.m_SignalGen.m_NumberOfCoils; c++)
       {
         int numSpikes = 0;
         for (auto ss : spikeSlice)
           if (std::get<0>(ss) == z && std::get<1>(ss) == c)
             ++numSpikes;
 
         // create k-sapce (inverse fourier transform slices)
         auto idft = itk::KspaceImageFilter< Float2DImageType::PixelType >::New();
         idft->SetCompartmentImages(compartment_slices);
         idft->SetT2(t2Vector);
         idft->SetT1(t1Vector);
         if (m_UseConstantRandSeed)
         {
           int linear_seed = g + num_gradient_volumes*z + num_gradient_volumes*compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)*c;
           idft->SetRandSeed(linear_seed);
         }
         idft->SetParameters(&m_Parameters);
         idft->SetZ((float)z-(float)( compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)
                                      -compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)%2 ) / 2.0);
         idft->SetZidx(z);
         idft->SetCoilPosition(coilPositions.at(c));
         idft->SetFiberBundle(m_FiberBundle);
         idft->SetTranslation(m_Translations.at(g));
         idft->SetRotationMatrix(m_RotationsInv.at(g));
         idft->SetDiffusionGradientDirection(m_Parameters.m_SignalGen.GetGradientDirection(g)*m_Parameters.m_SignalGen.GetBvalue()/1000.0);
         idft->SetSpikesPerSlice(numSpikes);
         idft->SetNumberOfThreads(in_threads);
         idft->Update();
 
 #pragma omp critical
         if (numSpikes>0)
         {
           m_SpikeLog += "Volume " + boost::lexical_cast<std::string>(g) + " Coil " + boost::lexical_cast<std::string>(c) + "\n";
           m_SpikeLog += idft->GetSpikeLog();
         }
 
         Complex2DImageType::Pointer fSlice;
         fSlice = idft->GetOutput();
 
         // fourier transform slice
         Complex2DImageType::Pointer newSlice;
         auto dft = itk::DftImageFilter< Float2DImageType::PixelType >::New();
         dft->SetInput(fSlice);
         dft->SetParameters(m_Parameters);
         dft->SetNumberOfThreads(in_threads);
         dft->Update();
         newSlice = dft->GetOutput();
 
         // put slice back into channel g
         for (unsigned int y=0; y<fSlice->GetLargestPossibleRegion().GetSize(1); y++)
           for (unsigned int x=0; x<fSlice->GetLargestPossibleRegion().GetSize(0); x++)
           {
             DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z;
             Complex2DImageType::IndexType index2D; index2D[0]=x; index2D[1]=y;
 
             Complex2DImageType::PixelType cPix = newSlice->GetPixel(index2D);
             double magn = sqrt(cPix.real()*cPix.real()+cPix.imag()*cPix.imag());
             double phase = 0;
             if (cPix.real()!=0)
               phase = atan( cPix.imag()/cPix.real() );
 
             DoubleDwiType::PixelType real_pix = m_OutputImagesReal.at(c)->GetPixel(index3D);
             real_pix[g] = cPix.real();
             m_OutputImagesReal.at(c)->SetPixel(index3D, real_pix);
 
             DoubleDwiType::PixelType imag_pix = m_OutputImagesImag.at(c)->GetPixel(index3D);
             imag_pix[g] = cPix.imag();
             m_OutputImagesImag.at(c)->SetPixel(index3D, imag_pix);
 
             DoubleDwiType::PixelType dwiPix = magnitudeDwiImage->GetPixel(index3D);
             DoubleDwiType::PixelType phasePix = m_PhaseImage->GetPixel(index3D);
 
             if (m_Parameters.m_SignalGen.m_NumberOfCoils>1)
             {
               dwiPix[g] += magn*magn;
               phasePix[g] += phase*phase;
             }
             else
             {
               dwiPix[g] = magn;
               phasePix[g] = phase;
             }
 
             //#pragma omp critical
             {
               magnitudeDwiImage->SetPixel(index3D, dwiPix);
               m_PhaseImage->SetPixel(index3D, phasePix);
 
               // k-space image
               if (g==0)
               {
                 DoubleDwiType::PixelType kspacePix = m_KspaceImage->GetPixel(index3D);
                 kspacePix[c] = idft->GetKSpaceImage()->GetPixel(index2D);
                 m_KspaceImage->SetPixel(index3D, kspacePix);
               }
             }
           }
       }
 
       if (m_Parameters.m_SignalGen.m_NumberOfCoils>1)
       {
         for (int y=0; y<static_cast<int>(magnitudeDwiImage->GetLargestPossibleRegion().GetSize(1)); y++)
           for (int x=0; x<static_cast<int>(magnitudeDwiImage->GetLargestPossibleRegion().GetSize(0)); x++)
           {
             DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z;
             DoubleDwiType::PixelType magPix = magnitudeDwiImage->GetPixel(index3D);
             magPix[g] = sqrt(magPix[g]/m_Parameters.m_SignalGen.m_NumberOfCoils);
 
             DoubleDwiType::PixelType phasePix = m_PhaseImage->GetPixel(index3D);
             phasePix[g] = sqrt(phasePix[g]/m_Parameters.m_SignalGen.m_NumberOfCoils);
 
             //#pragma omp critical
             {
               magnitudeDwiImage->SetPixel(index3D, magPix);
               m_PhaseImage->SetPixel(index3D, phasePix);
             }
           }
       }
 
       ++disp;
       unsigned long newTick = 50*disp.count()/disp.expected_count();
       for (unsigned long tick = 0; tick<(newTick-lastTick); tick++)
         PrintToLog("*", false, false, false);
       lastTick = newTick;
     }
   }
 
 
   PrintToLog("\n", false);
   return magnitudeDwiImage;
 }
 
 template< class PixelType >
 TractsToDWIImageFilter< PixelType >::ItkDoubleImgType::Pointer TractsToDWIImageFilter< PixelType >::
 NormalizeInsideMask(ItkDoubleImgType::Pointer image)
 {
   double max = itk::NumericTraits< double >::min();
   double min = itk::NumericTraits< double >::max();
 
   itk::ImageRegionIterator< ItkDoubleImgType > it(image, image->GetLargestPossibleRegion());
   while(!it.IsAtEnd())
   {
     if (m_Parameters.m_SignalGen.m_MaskImage.IsNotNull() && m_Parameters.m_SignalGen.m_MaskImage->GetPixel(it.GetIndex())<=0)
     {
       it.Set(0.0);
       ++it;
       continue;
     }
 
     if (it.Get()>max)
       max = it.Get();
     if (it.Get()<min)
       min = it.Get();
     ++it;
   }
   auto scaler = itk::ShiftScaleImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New();
   scaler->SetInput(image);
   scaler->SetShift(-min);
   scaler->SetScale(1.0/(max-min));
   scaler->Update();
   return scaler->GetOutput();
 }
 
 template< class PixelType >
 void TractsToDWIImageFilter< PixelType >::CheckVolumeFractionImages()
 {
   m_UseRelativeNonFiberVolumeFractions = false;
 
   // check for fiber volume fraction maps
   unsigned int fibVolImages = 0;
   for (std::size_t i=0; i<m_Parameters.m_FiberModelList.size(); i++)
   {
     if (m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage().IsNotNull())
     {
       PrintToLog("Using volume fraction map for fiber compartment " + boost::lexical_cast<std::string>(i+1), false);
       fibVolImages++;
     }
   }
 
   // check for non-fiber volume fraction maps
   unsigned int nonfibVolImages = 0;
   for (std::size_t i=0; i<m_Parameters.m_NonFiberModelList.size(); i++)
   {
     if (m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage().IsNotNull())
     {
       PrintToLog("Using volume fraction map for non-fiber compartment " + boost::lexical_cast<std::string>(i+1), false);
       nonfibVolImages++;
     }
   }
 
   // not all fiber compartments are using volume fraction maps
   // --> non-fiber volume fractions are assumed to be relative to the
   // non-fiber volume and not absolute voxel-volume fractions.
   // this means if two non-fiber compartments are used but only one of them
   // has an associated volume fraction map, the repesctive other volume fraction map
   // can be determined as inverse (1-val) of the present volume fraction map-
   if ( fibVolImages<m_Parameters.m_FiberModelList.size() && nonfibVolImages==1 && m_Parameters.m_NonFiberModelList.size()==2)
   {
     PrintToLog("Calculating missing non-fiber volume fraction image by inverting the other.\n"
                "Assuming non-fiber volume fraction images to contain values relative to the"
                " remaining non-fiber volume, not absolute values.", false);
 
     auto inverter = itk::InvertIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New();
     inverter->SetMaximum(1.0);
     if ( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage().IsNull()
          && m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage().IsNotNull() )
     {
       // m_Parameters.m_NonFiberModelList[1]->SetVolumeFractionImage(
       // NormalizeInsideMask( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage() ) );
       inverter->SetInput( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage() );
       inverter->Update();
       m_Parameters.m_NonFiberModelList[0]->SetVolumeFractionImage(inverter->GetOutput());
     }
     else if ( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage().IsNull()
               && m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage().IsNotNull() )
     {
       // m_Parameters.m_NonFiberModelList[0]->SetVolumeFractionImage(
       // NormalizeInsideMask( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage() ) );
       inverter->SetInput( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage() );
       inverter->Update();
       m_Parameters.m_NonFiberModelList[1]->SetVolumeFractionImage(inverter->GetOutput());
     }
     else
     {
       itkExceptionMacro("Something went wrong in automatically calculating the missing non-fiber volume fraction image!"
                         " Did you use two non fiber compartments but only one volume fraction image?"
                         " Then it should work and this error is really strange.");
     }
     m_UseRelativeNonFiberVolumeFractions = true;
 
     nonfibVolImages++;
   }
 
   // Up to two fiber compartments are allowed without volume fraction maps since the volume fractions can then be determined automatically
   if (m_Parameters.m_FiberModelList.size()>2
       && fibVolImages!=m_Parameters.m_FiberModelList.size())
     itkExceptionMacro("More than two fiber compartment selected but no corresponding volume fraction maps set!");
 
   // One non-fiber compartment is allowed without volume fraction map since the volume fraction can then be determined automatically
   if (m_Parameters.m_NonFiberModelList.size()>1
       && nonfibVolImages!=m_Parameters.m_NonFiberModelList.size())
     itkExceptionMacro("More than one non-fiber compartment selected but no volume fraction maps set!");
 
   if (fibVolImages<m_Parameters.m_FiberModelList.size() && nonfibVolImages>0)
   {
     PrintToLog("Not all fiber compartments are using an associated volume fraction image.\n"
                "Assuming non-fiber volume fraction images to contain values relative to the"
                " remaining non-fiber volume, not absolute values.", false);
     m_UseRelativeNonFiberVolumeFractions = true;
 
 //    mitk::LocaleSwitch localeSwitch("C");
     //        itk::ImageFileWriter<ItkDoubleImgType>::Pointer wr = itk::ImageFileWriter<ItkDoubleImgType>::New();
     //        wr->SetInput(m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage());
     //        wr->SetFileName("/local/volumefraction.nrrd");
     //        wr->Update();
   }
 
   // initialize the images that store the output volume fraction of each compartment
   m_VolumeFractions.clear();
   for (std::size_t i=0; i<m_Parameters.m_FiberModelList.size()+m_Parameters.m_NonFiberModelList.size(); i++)
   {
     auto doubleImg = ItkDoubleImgType::New();
     doubleImg->SetSpacing( m_WorkingSpacing );
     doubleImg->SetOrigin( m_WorkingOrigin );
     doubleImg->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
     doubleImg->SetLargestPossibleRegion( m_WorkingImageRegion );
     doubleImg->SetBufferedRegion( m_WorkingImageRegion );
     doubleImg->SetRequestedRegion( m_WorkingImageRegion );
     doubleImg->Allocate();
     doubleImg->FillBuffer(0);
     m_VolumeFractions.push_back(doubleImg);
   }
 }
 
 template< class PixelType >
 void TractsToDWIImageFilter< PixelType >::InitializeData()
 {
   m_Rotations.clear();
   m_Translations.clear();
   m_MotionLog = "";
   m_SpikeLog = "";
 
   // initialize output dwi image
   m_Parameters.m_SignalGen.m_CroppedRegion = m_Parameters.m_SignalGen.m_ImageRegion;
   if (m_Parameters.m_Misc.m_DoAddAliasing)
     m_Parameters.m_SignalGen.m_CroppedRegion.SetSize( 1, m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(1)
                                                       *m_Parameters.m_SignalGen.m_CroppingFactor);
 
   itk::Point<double,3> shiftedOrigin = m_Parameters.m_SignalGen.m_ImageOrigin;
   shiftedOrigin[1] += (m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1)
                        -m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(1))*m_Parameters.m_SignalGen.m_ImageSpacing[1]/2;
 
   m_OutputImage = OutputImageType::New();
   m_OutputImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing );
   m_OutputImage->SetOrigin( shiftedOrigin );
   m_OutputImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
   m_OutputImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
   m_OutputImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
   m_OutputImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
   m_OutputImage->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() );
   m_OutputImage->Allocate();
   typename OutputImageType::PixelType temp;
   temp.SetSize(m_Parameters.m_SignalGen.GetNumVolumes());
   temp.Fill(0.0);
   m_OutputImage->FillBuffer(temp);
 
   PrintToLog("Output image spacing: [" + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_ImageSpacing[0]) + "," + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_ImageSpacing[1]) + "," + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_ImageSpacing[2]) + "]", false);
   PrintToLog("Output image size: [" + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(0)) + "," + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(1)) + "," + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(2)) + "]", false);
 
   // images containing real and imaginary part of the dMRI signal for each coil
   m_OutputImagesReal.clear();
   m_OutputImagesImag.clear();
   for (unsigned int i=0; i<m_Parameters.m_SignalGen.m_NumberOfCoils; ++i)
   {
     typename DoubleDwiType::Pointer outputImageReal = DoubleDwiType::New();
     outputImageReal->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing );
     outputImageReal->SetOrigin( shiftedOrigin );
     outputImageReal->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
     outputImageReal->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
     outputImageReal->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
     outputImageReal->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
     outputImageReal->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() );
     outputImageReal->Allocate();
     outputImageReal->FillBuffer(temp);
     m_OutputImagesReal.push_back(outputImageReal);
 
     typename DoubleDwiType::Pointer outputImageImag = DoubleDwiType::New();
     outputImageImag->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing );
     outputImageImag->SetOrigin( shiftedOrigin );
     outputImageImag->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
     outputImageImag->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
     outputImageImag->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
     outputImageImag->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
     outputImageImag->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() );
     outputImageImag->Allocate();
     outputImageImag->FillBuffer(temp);
     m_OutputImagesImag.push_back(outputImageImag);
   }
 
   // Apply in-plane upsampling for Gibbs ringing artifact
   double upsampling = 1;
   if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging && m_Parameters.m_SignalGen.m_ZeroRinging==0)
     upsampling = 2;
   m_WorkingSpacing = m_Parameters.m_SignalGen.m_ImageSpacing;
   m_WorkingSpacing[0] /= upsampling;
   m_WorkingSpacing[1] /= upsampling;
   m_WorkingImageRegion = m_Parameters.m_SignalGen.m_ImageRegion;
   m_WorkingImageRegion.SetSize(0, m_Parameters.m_SignalGen.m_ImageRegion.GetSize()[0]*upsampling);
   m_WorkingImageRegion.SetSize(1, m_Parameters.m_SignalGen.m_ImageRegion.GetSize()[1]*upsampling);
   m_WorkingOrigin = m_Parameters.m_SignalGen.m_ImageOrigin;
   m_WorkingOrigin[0] -= m_Parameters.m_SignalGen.m_ImageSpacing[0]/2;
   m_WorkingOrigin[0] += m_WorkingSpacing[0]/2;
   m_WorkingOrigin[1] -= m_Parameters.m_SignalGen.m_ImageSpacing[1]/2;
   m_WorkingOrigin[1] += m_WorkingSpacing[1]/2;
   m_WorkingOrigin[2] -= m_Parameters.m_SignalGen.m_ImageSpacing[2]/2;
   m_WorkingOrigin[2] += m_WorkingSpacing[2]/2;
   m_VoxelVolume = m_WorkingSpacing[0]*m_WorkingSpacing[1]*m_WorkingSpacing[2];
 
   PrintToLog("Working image spacing: [" + boost::lexical_cast<std::string>(m_WorkingSpacing[0]) + "," + boost::lexical_cast<std::string>(m_WorkingSpacing[1]) + "," + boost::lexical_cast<std::string>(m_WorkingSpacing[2]) + "]", false);
   PrintToLog("Working image size: [" + boost::lexical_cast<std::string>(m_WorkingImageRegion.GetSize(0)) + "," + boost::lexical_cast<std::string>(m_WorkingImageRegion.GetSize(1)) + "," + boost::lexical_cast<std::string>(m_WorkingImageRegion.GetSize(2)) + "]", false);
 
   // generate double images to store the individual compartment signals
   m_CompartmentImages.clear();
   int numFiberCompartments = m_Parameters.m_FiberModelList.size();
   int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size();
   for (int i=0; i<numFiberCompartments+numNonFiberCompartments; i++)
   {
     auto doubleDwi = DoubleDwiType::New();
     doubleDwi->SetSpacing( m_WorkingSpacing );
     doubleDwi->SetOrigin( m_WorkingOrigin );
     doubleDwi->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
     doubleDwi->SetLargestPossibleRegion( m_WorkingImageRegion );
     doubleDwi->SetBufferedRegion( m_WorkingImageRegion );
     doubleDwi->SetRequestedRegion( m_WorkingImageRegion );
     doubleDwi->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() );
     doubleDwi->Allocate();
     DoubleDwiType::PixelType pix;
     pix.SetSize(m_Parameters.m_SignalGen.GetNumVolumes());
     pix.Fill(0.0);
     doubleDwi->FillBuffer(pix);
     m_CompartmentImages.push_back(doubleDwi);
   }
 
   if (m_FiberBundle.IsNull() && m_InputImage.IsNotNull())
   {
     m_CompartmentImages.clear();
 
     m_Parameters.m_SignalGen.m_DoAddMotion = false;
     m_Parameters.m_SignalGen.m_DoSimulateRelaxation = false;
 
     PrintToLog("Simulating acquisition for input diffusion-weighted image.", false);
     auto caster = itk::CastImageFilter< OutputImageType, DoubleDwiType >::New();
     caster->SetInput(m_InputImage);
     caster->Update();
 
     if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging && m_Parameters.m_SignalGen.m_ZeroRinging==0)
     {
       PrintToLog("Upsampling input diffusion-weighted image for Gibbs ringing simulation.", false);
 
       auto resampler = itk::ResampleDwiImageFilter< double >::New();
       resampler->SetInput(caster->GetOutput());
       itk::Vector< double, 3 > samplingFactor;
       samplingFactor[0] = upsampling;
       samplingFactor[1] = upsampling;
       samplingFactor[2] = 1;
       resampler->SetSamplingFactor(samplingFactor);
       resampler->SetInterpolation(itk::ResampleDwiImageFilter< double >::Interpolate_WindowedSinc);
       resampler->Update();
       m_CompartmentImages.push_back(resampler->GetOutput());
     }
     else
       m_CompartmentImages.push_back(caster->GetOutput());
 
     VectorType translation; translation.Fill(0.0);
     MatrixType rotation; rotation.SetIdentity();
     for (unsigned int g=0; g<m_Parameters.m_SignalGen.GetNumVolumes(); g++)
     {
       m_Rotations.push_back(rotation);
       m_RotationsInv.push_back(rotation);
       m_Translations.push_back(translation);
     }
   }
 
   // resample mask image and frequency map to fit upsampled geometry
   if (m_Parameters.m_SignalGen.m_MaskImage.IsNotNull() && m_Parameters.m_SignalGen.m_MaskImage->GetLargestPossibleRegion()!=m_WorkingImageRegion)
   {
     PrintToLog("Resampling tissue mask", false);
     // rescale mask image (otherwise there are problems with the resampling)
     auto rescaler = itk::RescaleIntensityImageFilter<ItkUcharImgType,ItkUcharImgType>::New();
     rescaler->SetInput(0,m_Parameters.m_SignalGen.m_MaskImage);
     rescaler->SetOutputMaximum(100);
     rescaler->SetOutputMinimum(0);
     rescaler->Update();
 
     // resample mask image
     auto resampler = itk::ResampleImageFilter<ItkUcharImgType, ItkUcharImgType>::New();
     resampler->SetInput(rescaler->GetOutput());
     resampler->SetSize(m_WorkingImageRegion.GetSize());
     resampler->SetOutputSpacing(m_WorkingSpacing);
     resampler->SetOutputOrigin(m_WorkingOrigin);
     resampler->SetOutputDirection(m_Parameters.m_SignalGen.m_ImageDirection);
     resampler->SetOutputStartIndex ( m_WorkingImageRegion.GetIndex() );
     auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction<ItkUcharImgType, double>::New();
     resampler->SetInterpolator(nn_interpolator);
     resampler->Update();
     m_Parameters.m_SignalGen.m_MaskImage = resampler->GetOutput();
   }
   // resample frequency map
   if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull() && m_Parameters.m_SignalGen.m_FrequencyMap->GetLargestPossibleRegion()!=m_WorkingImageRegion)
   {
     PrintToLog("Resampling frequency map", false);
     auto resampler = itk::ResampleImageFilter<ItkFloatImgType, ItkFloatImgType>::New();
     resampler->SetInput(m_Parameters.m_SignalGen.m_FrequencyMap);
     resampler->SetSize(m_WorkingImageRegion.GetSize());
     resampler->SetOutputSpacing(m_WorkingSpacing);
     resampler->SetOutputOrigin(m_WorkingOrigin);
     resampler->SetOutputDirection(m_Parameters.m_SignalGen.m_ImageDirection);
     resampler->SetOutputStartIndex ( m_WorkingImageRegion.GetIndex() );
     auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction<ItkFloatImgType, double>::New();
     resampler->SetInterpolator(nn_interpolator);
     resampler->Update();
     m_Parameters.m_SignalGen.m_FrequencyMap = resampler->GetOutput();
   }
 
   m_MaskImageSet = true;
   if (m_Parameters.m_SignalGen.m_MaskImage.IsNull())
   {
     // no input tissue mask is set -> create default
     PrintToLog("No tissue mask set", false);
     m_Parameters.m_SignalGen.m_MaskImage = ItkUcharImgType::New();
     m_Parameters.m_SignalGen.m_MaskImage->SetSpacing( m_WorkingSpacing );
     m_Parameters.m_SignalGen.m_MaskImage->SetOrigin( m_WorkingOrigin );
     m_Parameters.m_SignalGen.m_MaskImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
     m_Parameters.m_SignalGen.m_MaskImage->SetLargestPossibleRegion( m_WorkingImageRegion );
     m_Parameters.m_SignalGen.m_MaskImage->SetBufferedRegion( m_WorkingImageRegion );
     m_Parameters.m_SignalGen.m_MaskImage->SetRequestedRegion( m_WorkingImageRegion );
     m_Parameters.m_SignalGen.m_MaskImage->Allocate();
     m_Parameters.m_SignalGen.m_MaskImage->FillBuffer(100);
     m_MaskImageSet = false;
   }
   else
   {
     PrintToLog("Using tissue mask", false);
   }
 
   if (m_Parameters.m_SignalGen.m_DoAddMotion)
   {
     if (m_Parameters.m_SignalGen.m_DoRandomizeMotion)
     {
       PrintToLog("Random motion artifacts:", false);
       PrintToLog("Maximum rotation: +/-" + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_Rotation) + "°", false);
       PrintToLog("Maximum translation: +/-" + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_Translation) + "mm", false);
     }
     else
     {
       PrintToLog("Linear motion artifacts:", false);
       PrintToLog("Maximum rotation: " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_Rotation) + "°", false);
       PrintToLog("Maximum translation: " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_Translation) + "mm", false);
     }
   }
   if ( m_Parameters.m_SignalGen.m_MotionVolumes.empty() )
   {
     // no motion in first volume
     m_Parameters.m_SignalGen.m_MotionVolumes.push_back(false);
 
     // motion in all other volumes
     while ( m_Parameters.m_SignalGen.m_MotionVolumes.size() < m_Parameters.m_SignalGen.GetNumVolumes() )
     {
       m_Parameters.m_SignalGen.m_MotionVolumes.push_back(true);
     }
   }
   // we need to know for every volume if there is motion. if this information is missing, then set corresponding fal to false
   while ( m_Parameters.m_SignalGen.m_MotionVolumes.size()<m_Parameters.m_SignalGen.GetNumVolumes() )
   {
     m_Parameters.m_SignalGen.m_MotionVolumes.push_back(false);
   }
 
   m_NumMotionVolumes = 0;
   for (unsigned int i=0; i<m_Parameters.m_SignalGen.GetNumVolumes(); ++i)
   {
     if (m_Parameters.m_SignalGen.m_MotionVolumes[i])
       ++m_NumMotionVolumes;
   }
   m_MotionCounter = 0;
 
   // creat image to hold transformed mask (motion artifact)
   m_TransformedMaskImage = ItkUcharImgType::New();
   auto duplicator = itk::ImageDuplicator<ItkUcharImgType>::New();
   duplicator->SetInputImage(m_Parameters.m_SignalGen.m_MaskImage);
   duplicator->Update();
   m_TransformedMaskImage = duplicator->GetOutput();
 
   // second upsampling needed for motion artifacts
   ImageRegion<3>      upsampledImageRegion = m_WorkingImageRegion;
   VectorType    upsampledSpacing = m_WorkingSpacing;
   upsampledSpacing[0] /= 4;
   upsampledSpacing[1] /= 4;
   upsampledSpacing[2] /= 4;
   upsampledImageRegion.SetSize(0, m_WorkingImageRegion.GetSize()[0]*4);
   upsampledImageRegion.SetSize(1, m_WorkingImageRegion.GetSize()[1]*4);
   upsampledImageRegion.SetSize(2, m_WorkingImageRegion.GetSize()[2]*4);
   itk::Point<double,3> upsampledOrigin = m_WorkingOrigin;
   upsampledOrigin[0] -= m_WorkingSpacing[0]/2;
   upsampledOrigin[0] += upsampledSpacing[0]/2;
   upsampledOrigin[1] -= m_WorkingSpacing[1]/2;
   upsampledOrigin[1] += upsampledSpacing[1]/2;
   upsampledOrigin[2] -= m_WorkingSpacing[2]/2;
   upsampledOrigin[2] += upsampledSpacing[2]/2;
 
   m_UpsampledMaskImage = ItkUcharImgType::New();
   auto upsampler = itk::ResampleImageFilter<ItkUcharImgType, ItkUcharImgType>::New();
   upsampler->SetInput(m_Parameters.m_SignalGen.m_MaskImage);
   upsampler->SetOutputParametersFromImage(m_Parameters.m_SignalGen.m_MaskImage);
   upsampler->SetSize(upsampledImageRegion.GetSize());
   upsampler->SetOutputSpacing(upsampledSpacing);
   upsampler->SetOutputOrigin(upsampledOrigin);
 
   auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction<ItkUcharImgType, double>::New();
   upsampler->SetInterpolator(nn_interpolator);
   upsampler->Update();
   m_UpsampledMaskImage = upsampler->GetOutput();
 }
 
 template< class PixelType >
 void TractsToDWIImageFilter< PixelType >::InitializeFiberData()
 {
   m_mmRadius = m_Parameters.m_SignalGen.m_AxonRadius/1000;
 
   auto caster = itk::CastImageFilter< itk::Image<unsigned char, 3>, itk::Image<double, 3> >::New();
   caster->SetInput(m_TransformedMaskImage);
   caster->Update();
 
   vtkSmartPointer<vtkFloatArray> weights = m_FiberBundle->GetFiberWeights();
   float mean_weight = 0;
   for (int i=0; i<weights->GetSize(); i++)
     mean_weight += weights->GetValue(i);
   mean_weight /= weights->GetSize();
 
   if (mean_weight>0.000001)
     for (int i=0; i<weights->GetSize(); i++)
       m_FiberBundle->SetFiberWeight(i, weights->GetValue(i)/mean_weight);
   else
     PrintToLog("\nWarning: streamlines have VERY low weights. Average weight: " + boost::lexical_cast<std::string>(mean_weight) + ". Possible source of calculation errors.", false, true, true);
 
 
   auto density_calculator = itk::TractDensityImageFilter< itk::Image<double, 3> >::New();
   density_calculator->SetFiberBundle(m_FiberBundle);
   density_calculator->SetInputImage(caster->GetOutput());
   density_calculator->SetBinaryOutput(false);
   density_calculator->SetUseImageGeometry(true);
   density_calculator->SetOutputAbsoluteValues(true);
   density_calculator->Update();
   double max_density = density_calculator->GetMaxDensity();
 
   double voxel_volume = m_WorkingSpacing[0]*m_WorkingSpacing[1]*m_WorkingSpacing[2];
   if (m_mmRadius>0)
   {
     std::stringstream stream;
     stream << std::fixed << setprecision(2) << itk::Math::pi*m_mmRadius*m_mmRadius*max_density;
     std::string s = stream.str();
     PrintToLog("\nMax. fiber volume: " + s + "mm².", false, true, true);
 
     {
       double full_radius = 1000*std::sqrt(voxel_volume/(max_density*itk::Math::pi));
       std::stringstream stream;
       stream << std::fixed << setprecision(2) << full_radius;
       std::string s = stream.str();
       PrintToLog("\nA full fiber voxel corresponds to a fiber radius of ~" + s + "µm, given the current fiber configuration.", false, true, true);
     }
   }
   else
   {
     m_mmRadius = std::sqrt(voxel_volume/(max_density*itk::Math::pi));
     std::stringstream stream;
     stream << std::fixed << setprecision(2) << m_mmRadius*1000;
     std::string s = stream.str();
     PrintToLog("\nSetting fiber radius to " + s + "µm to obtain full voxel.", false, true, true);
   }
 
   // a second fiber bundle is needed to store the transformed version of the m_FiberBundleWorkingCopy
   m_FiberBundleTransformed = m_FiberBundle->GetDeepCopy();
 }
 
 
 template< class PixelType >
 bool TractsToDWIImageFilter< PixelType >::PrepareLogFile()
 {
   if(m_Logfile.is_open())
     m_Logfile.close();
 
   std::string filePath;
   std::string fileName;
 
   // Get directory name:
   if (m_Parameters.m_Misc.m_OutputPath.size() > 0)
   {
     filePath = m_Parameters.m_Misc.m_OutputPath;
     if( *(--(filePath.cend())) != '/')
     {
       filePath.push_back('/');
     }
   }
   else
     return false;
 
   // Get file name:
   if( ! m_Parameters.m_Misc.m_ResultNode->GetName().empty() )
   {
     fileName = m_Parameters.m_Misc.m_ResultNode->GetName();
   }
   else
   {
     fileName = "";
   }
 
   if( ! m_Parameters.m_Misc.m_OutputPrefix.empty() )
   {
     fileName = m_Parameters.m_Misc.m_OutputPrefix + fileName;
   }
 
   try
   {
     m_Logfile.open( ( filePath + '/' + fileName + ".log" ).c_str() );
   }
   catch (const std::ios_base::failure &fail)
   {
     MITK_ERROR << "itkTractsToDWIImageFilter.cpp: Exception " << fail.what()
                << " while trying to open file" << filePath << '/' << fileName << ".log";
     return false;
   }
 
   if ( m_Logfile.is_open() )
   {
     PrintToLog( "Logfile: " + filePath + '/' + fileName + ".log", false );
     return true;
   }
   else
     return false;
 }
 
 
 template< class PixelType >
 void TractsToDWIImageFilter< PixelType >::GenerateData()
 {
   PrintToLog("\n**********************************************", false);
   // prepare logfile
   PrepareLogFile();
   PrintToLog("Starting Fiberfox dMRI simulation");
 
   m_TimeProbe.Start();
 
   // check input data
   if (m_FiberBundle.IsNull() && m_InputImage.IsNull())
     itkExceptionMacro("Input fiber bundle and input diffusion-weighted image is nullptr!");
 
   if (m_Parameters.m_FiberModelList.empty() && m_InputImage.IsNull())
     itkExceptionMacro("No diffusion model for fiber compartments defined and input diffusion-weighted"
                       " image is nullptr! At least one fiber compartment is necessary to simulate diffusion.");
 
   if (m_Parameters.m_NonFiberModelList.empty() && m_InputImage.IsNull())
     itkExceptionMacro("No diffusion model for non-fiber compartments defined and input diffusion-weighted"
                       " image is nullptr! At least one non-fiber compartment is necessary to simulate diffusion.");
 
   if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume)  // no partial volume? remove all but first fiber compartment
     while (m_Parameters.m_FiberModelList.size()>1)
       m_Parameters.m_FiberModelList.pop_back();
 
   if (!m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition)  // No upsampling of input image needed if no k-space simulation is performed
     m_Parameters.m_SignalGen.m_DoAddGibbsRinging = false;
 
   if (m_UseConstantRandSeed)  // always generate the same random numbers?
     m_RandGen->SetSeed(0);
   else
     m_RandGen->SetSeed();
 
   InitializeData();
   if ( m_FiberBundle.IsNotNull() )    // if no fiber bundle is found, we directly proceed to the k-space acquisition simulation
   {
     CheckVolumeFractionImages();
     InitializeFiberData();
 
     int numFiberCompartments = m_Parameters.m_FiberModelList.size();
     int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size();
 
     double maxVolume = 0;
     unsigned long lastTick = 0;
     int signalModelSeed = m_RandGen->GetIntegerVariate();
 
     PrintToLog("\n", false, false);
     PrintToLog("Generating " + boost::lexical_cast<std::string>(numFiberCompartments+numNonFiberCompartments)
                + "-compartment diffusion-weighted signal.");
 
     std::vector< int > bVals = m_Parameters.m_SignalGen.GetBvalues();
     PrintToLog("b-values: ", false, false, true);
     for (auto v : bVals)
       PrintToLog(boost::lexical_cast<std::string>(v) + " ", false, false, true);
     PrintToLog("\nVolumes: " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.GetNumVolumes()), false, true, true);
 
     PrintToLog("\n", false, false, true);
     PrintToLog("\n", false, false, true);
 
     unsigned int image_size_x = m_WorkingImageRegion.GetSize(0);
     unsigned int region_size_y = m_WorkingImageRegion.GetSize(1);
     unsigned int num_gradients = m_Parameters.m_SignalGen.GetNumVolumes();
     int numFibers = m_FiberBundle->GetNumFibers();
     boost::progress_display disp(numFibers*num_gradients);
 
     if (m_FiberBundle->GetMeanFiberLength()<5.0)
       omp_set_num_threads(2);
 
     PrintToLog("0%   10   20   30   40   50   60   70   80   90   100%", false, true, false);
     PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false);
 
     for (unsigned int g=0; g<num_gradients; ++g)
     {
       // move fibers
       SimulateMotion(g);
 
       // Set signal model random generator seeds to get same configuration in each voxel
       for (std::size_t i=0; i<m_Parameters.m_FiberModelList.size(); i++)
         m_Parameters.m_FiberModelList.at(i)->SetSeed(signalModelSeed);
       for (std::size_t i=0; i<m_Parameters.m_NonFiberModelList.size(); i++)
         m_Parameters.m_NonFiberModelList.at(i)->SetSeed(signalModelSeed);
 
       // storing voxel-wise intra-axonal volume in mm³
       auto intraAxonalVolumeImage = ItkDoubleImgType::New();
       intraAxonalVolumeImage->SetSpacing( m_WorkingSpacing );
       intraAxonalVolumeImage->SetOrigin( m_WorkingOrigin );
       intraAxonalVolumeImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
       intraAxonalVolumeImage->SetLargestPossibleRegion( m_WorkingImageRegion );
       intraAxonalVolumeImage->SetBufferedRegion( m_WorkingImageRegion );
       intraAxonalVolumeImage->SetRequestedRegion( m_WorkingImageRegion );
       intraAxonalVolumeImage->Allocate();
       intraAxonalVolumeImage->FillBuffer(0);
       maxVolume = 0;
       double* intraAxBuffer = intraAxonalVolumeImage->GetBufferPointer();
 
       if (this->GetAbortGenerateData())
         continue;
 
       vtkPolyData* fiberPolyData = m_FiberBundleTransformed->GetFiberPolyData();
       // generate fiber signal (if there are any fiber models present)
       if (!m_Parameters.m_FiberModelList.empty())
       {
         std::vector< double* > buffers;
         for (unsigned int i=0; i<m_CompartmentImages.size(); ++i)
           buffers.push_back(m_CompartmentImages.at(i)->GetBufferPointer());
 #pragma omp parallel for
         for( int i=0; i<numFibers; ++i )
         {
           if (this->GetAbortGenerateData())
             continue;
 
           float fiberWeight = m_FiberBundleTransformed->GetFiberWeight(i);
 
           int numPoints = -1;
           std::vector< itk::Vector<double, 3> > points_copy;
 #pragma omp critical
           {
             vtkCell* cell = fiberPolyData->GetCell(i);
             numPoints = cell->GetNumberOfPoints();
             vtkPoints* points = cell->GetPoints();
             for (int j=0; j<numPoints; j++)
               points_copy.push_back(GetItkVector(points->GetPoint(j)));
           }
 
           if (numPoints<2)
             continue;
 
           double seg_volume = fiberWeight*itk::Math::pi*m_mmRadius*m_mmRadius;
           for( int j=0; j<numPoints - 1; ++j)
           {
             if (this->GetAbortGenerateData())
             {
               j=numPoints;
               continue;
             }
 
             itk::Vector<double> v = points_copy.at(j);
 
             itk::Vector<double, 3> dir = points_copy.at(j+1)-v;
             if ( dir.GetSquaredNorm()<0.0001 || dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2] )
               continue;
             dir.Normalize();
 
             itk::Point<float, 3> startVertex = points_copy.at(j);
             itk::Index<3> startIndex;
             itk::ContinuousIndex<float, 3> startIndexCont;
             m_TransformedMaskImage->TransformPhysicalPointToIndex(startVertex, startIndex);
             m_TransformedMaskImage->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont);
 
             itk::Point<float, 3> endVertex = points_copy.at(j+1);
             itk::Index<3> endIndex;
             itk::ContinuousIndex<float, 3> endIndexCont;
             m_TransformedMaskImage->TransformPhysicalPointToIndex(endVertex, endIndex);
             m_TransformedMaskImage->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont);
 
             std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(m_WorkingSpacing, startIndex, endIndex, startIndexCont, endIndexCont);
 
             // generate signal for each fiber compartment
             for (int k=0; k<numFiberCompartments; ++k)
             {
               double signal_add = m_Parameters.m_FiberModelList[k]->SimulateMeasurement(g, dir)*seg_volume;
               for (std::pair< itk::Index<3>, double > seg : segments)
               {
                 if (!m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(seg.first) || m_TransformedMaskImage->GetPixel(seg.first)<=0)
                   continue;
 
                 double seg_signal = seg.second*signal_add;
 
                 unsigned int linear_index = g + num_gradients*seg.first[0] + num_gradients*image_size_x*seg.first[1] + num_gradients*image_size_x*region_size_y*seg.first[2];
 
                 // update dMRI volume
 #pragma omp atomic
                 buffers[k][linear_index] += seg_signal;
 
                 // update fiber volume image
                 if (k==0)
                 {
                   linear_index = seg.first[0] + image_size_x*seg.first[1] + image_size_x*region_size_y*seg.first[2];
 #pragma omp atomic
                   intraAxBuffer[linear_index] += seg.second*seg_volume;
                   double vol = intraAxBuffer[linear_index];
                   if (vol>maxVolume) { maxVolume = vol; }
                 }
               }
 
             }
           }
 
 #pragma omp critical
           {
             // progress report
             ++disp;
             unsigned long newTick = 50*disp.count()/disp.expected_count();
             for (unsigned int tick = 0; tick<(newTick-lastTick); ++tick)
               PrintToLog("*", false, false, false);
             lastTick = newTick;
           }
         }
       }
 
       // axon radius not manually defined --> set fullest voxel (maxVolume) to full fiber voxel
       double density_correctiony_global = 1.0;
       if (m_Parameters.m_SignalGen.m_AxonRadius<0.0001)
         density_correctiony_global = m_VoxelVolume/maxVolume;
 
       // generate non-fiber signal
       ImageRegionIterator<ItkUcharImgType> it3(m_TransformedMaskImage, m_TransformedMaskImage->GetLargestPossibleRegion());
       while(!it3.IsAtEnd())
       {
         if (it3.Get()>0)
         {
           DoubleDwiType::IndexType index = it3.GetIndex();
           double iAxVolume = intraAxonalVolumeImage->GetPixel(index);
 
           // get non-transformed point (remove headmotion tranformation)
           // this point lives in the volume fraction image space
           itk::Point<float, 3> volume_fraction_point;
           if ( m_Parameters.m_SignalGen.m_DoAddMotion )
             volume_fraction_point = GetMovedPoint(index, false);
           else
             m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, volume_fraction_point);
 
           if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume)
           {
             if (iAxVolume>0.0001) // scale fiber compartment to voxel
             {
               DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index);
               pix[g] *= m_VoxelVolume/iAxVolume;
               m_CompartmentImages.at(0)->SetPixel(index, pix);
 
               if (g==0)
                 m_VolumeFractions.at(0)->SetPixel(index, 1);
             }
             else
             {
               DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index);
               pix[g] = 0;
               m_CompartmentImages.at(0)->SetPixel(index, pix);
               SimulateExtraAxonalSignal(index, volume_fraction_point, 0, g);
             }
           }
           else
           {
             // manually defined axon radius and voxel overflow --> rescale to voxel volume
             if ( m_Parameters.m_SignalGen.m_AxonRadius>=0.0001 && iAxVolume>m_VoxelVolume )
             {
               for (int i=0; i<numFiberCompartments; ++i)
               {
                 DoubleDwiType::PixelType pix = m_CompartmentImages.at(i)->GetPixel(index);
                 pix[g] *= m_VoxelVolume/iAxVolume;
                 m_CompartmentImages.at(i)->SetPixel(index, pix);
               }
               iAxVolume = m_VoxelVolume;
             }
 
             // if volume fraction image is set use it, otherwise use global scaling factor
             double density_correction_voxel = density_correctiony_global;
             if ( m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()!=nullptr && iAxVolume>0.0001 )
             {
               m_DoubleInterpolator->SetInputImage(m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage());
               double volume_fraction = mitk::imv::GetImageValue<double>(volume_fraction_point, true, m_DoubleInterpolator);
               if (volume_fraction<0)
                 mitkThrow() << "Volume fraction image (index 1) contains negative values (intra-axonal compartment)!";
               density_correction_voxel = m_VoxelVolume*volume_fraction/iAxVolume; // remove iAxVolume sclaing and scale to volume_fraction
             }
             else if (m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()!=nullptr)
               density_correction_voxel = 0.0;
 
             // adjust intra-axonal compartment volume by density correction factor
             DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index);
             pix[g] *= density_correction_voxel;
             m_CompartmentImages.at(0)->SetPixel(index, pix);
 
             // normalize remaining fiber volume fractions (they are rescaled in SimulateExtraAxonalSignal)
             if (iAxVolume>0.0001)
             {
               for (int i=1; i<numFiberCompartments; i++)
               {
                 DoubleDwiType::PixelType pix = m_CompartmentImages.at(i)->GetPixel(index);
                 pix[g] /= iAxVolume;
                 m_CompartmentImages.at(i)->SetPixel(index, pix);
               }
             }
             else
             {
               for (int i=1; i<numFiberCompartments; i++)
               {
                 DoubleDwiType::PixelType pix = m_CompartmentImages.at(i)->GetPixel(index);
                 pix[g] = 0;
                 m_CompartmentImages.at(i)->SetPixel(index, pix);
               }
             }
 
             iAxVolume = density_correction_voxel*iAxVolume; // new intra-axonal volume = old intra-axonal volume * correction factor
 
             // simulate other compartments
             SimulateExtraAxonalSignal(index, volume_fraction_point, iAxVolume, g);
           }
         }
         ++it3;
       }
     }
 
     PrintToLog("\n", false);
   }
 
   if (this->GetAbortGenerateData())
   {
     PrintToLog("\n", false, false);
     PrintToLog("Simulation aborted");
     return;
   }
 
   DoubleDwiType::Pointer doubleOutImage;
   double signalScale = m_Parameters.m_SignalGen.m_SignalScale;
   if ( m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition ) // do k-space stuff
   {
     PrintToLog("\n", false, false);
     PrintToLog("Simulating k-space acquisition using "
                +boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_NumberOfCoils)
                +" coil(s)");
 
     switch (m_Parameters.m_SignalGen.m_AcquisitionType)
     {
     case SignalGenerationParameters::SingleShotEpi:
     {
       PrintToLog("Acquisition type: single shot EPI", false);
       break;
     }
     case SignalGenerationParameters::ConventionalSpinEcho:
     {
       PrintToLog("Acquisition type: conventional spin echo (one RF pulse per line) with cartesian k-space trajectory", false);
       break;
     }
     case SignalGenerationParameters::FastSpinEcho:
     {
       PrintToLog("Acquisition type: fast spin echo (one RF pulse per ETL lines) with cartesian k-space trajectory (ETL: " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_EchoTrainLength) + ")", false);
+      PrintToLog("Effective TE: " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_tEcho*m_Parameters.m_SignalGen.m_EchoTrainLength/2), false);
       break;
     }
     default:
     {
       PrintToLog("Acquisition type: single shot EPI", false);
       break;
     }
     }
     if(m_Parameters.m_SignalGen.m_tInv>0)
       PrintToLog("Using inversion pulse with TI " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_tInv) + "ms", false);
 
     if (m_Parameters.m_SignalGen.m_DoSimulateRelaxation)
       PrintToLog("Simulating signal relaxation", false);
     if (m_Parameters.m_SignalGen.m_NoiseVariance>0 && m_Parameters.m_Misc.m_DoAddNoise)
       PrintToLog("Simulating complex Gaussian noise: " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_NoiseVariance), false);
     if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull() && m_Parameters.m_Misc.m_DoAddDistortions)
       PrintToLog("Simulating distortions", false);
     if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging)
     {
       if (m_Parameters.m_SignalGen.m_ZeroRinging > 0)
         PrintToLog("Simulating ringing artifacts by zeroing " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_ZeroRinging) + "% of k-space frequencies", false);
       else
         PrintToLog("Simulating ringing artifacts by cropping high resolution inputs during k-space simulation", false);
     }
     if (m_Parameters.m_Misc.m_DoAddEddyCurrents && m_Parameters.m_SignalGen.m_EddyStrength>0)
       PrintToLog("Simulating eddy currents: " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_EddyStrength), false);
     if (m_Parameters.m_Misc.m_DoAddSpikes && m_Parameters.m_SignalGen.m_Spikes>0)
       PrintToLog("Simulating spikes: " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_Spikes), false);
     if (m_Parameters.m_Misc.m_DoAddAliasing && m_Parameters.m_SignalGen.m_CroppingFactor<1.0)
       PrintToLog("Simulating aliasing: " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_CroppingFactor), false);
     if (m_Parameters.m_Misc.m_DoAddGhosts && m_Parameters.m_SignalGen.m_KspaceLineOffset>0)
       PrintToLog("Simulating ghosts: " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_KspaceLineOffset), false);
 
     doubleOutImage = SimulateKspaceAcquisition(m_CompartmentImages);
     signalScale = 1; // already scaled in SimulateKspaceAcquisition()
   }
   else    // don't do k-space stuff, just sum compartments
   {
     PrintToLog("Summing compartments");
     doubleOutImage = m_CompartmentImages.at(0);
 
     for (unsigned int i=1; i<m_CompartmentImages.size(); i++)
     {
       auto adder = itk::AddImageFilter< DoubleDwiType, DoubleDwiType, DoubleDwiType>::New();
       adder->SetInput1(doubleOutImage);
       adder->SetInput2(m_CompartmentImages.at(i));
       adder->Update();
       doubleOutImage = adder->GetOutput();
     }
   }
   if (this->GetAbortGenerateData())
   {
     PrintToLog("\n", false, false);
     PrintToLog("Simulation aborted");
     return;
   }
 
   PrintToLog("Finalizing image");
   if (m_Parameters.m_SignalGen.m_DoAddDrift && m_Parameters.m_SignalGen.m_Drift>0.0)
     PrintToLog("Adding signal drift: " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_Drift), false);
   if (signalScale>1)
     PrintToLog("Scaling signal", false);
   if (m_Parameters.m_NoiseModel)
     PrintToLog("Adding noise: " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_NoiseVariance), false);
   ImageRegionIterator<OutputImageType> it4 (m_OutputImage, m_OutputImage->GetLargestPossibleRegion());
   DoubleDwiType::PixelType signal; signal.SetSize(m_Parameters.m_SignalGen.GetNumVolumes());
   boost::progress_display disp2(m_OutputImage->GetLargestPossibleRegion().GetNumberOfPixels());
 
   PrintToLog("0%   10   20   30   40   50   60   70   80   90   100%", false, true, false);
   PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false);
   int lastTick = 0;
 
   while(!it4.IsAtEnd())
   {
     if (this->GetAbortGenerateData())
     {
       PrintToLog("\n", false, false);
       PrintToLog("Simulation aborted");
       return;
     }
 
     ++disp2;
     unsigned long newTick = 50*disp2.count()/disp2.expected_count();
     for (unsigned long tick = 0; tick<(newTick-lastTick); tick++)
       PrintToLog("*", false, false, false);
     lastTick = newTick;
 
     typename OutputImageType::IndexType index = it4.GetIndex();
     signal = doubleOutImage->GetPixel(index)*signalScale;
 
     for (unsigned int i=0; i<signal.Size(); i++)
     {
       if (m_Parameters.m_SignalGen.m_DoAddDrift)
       {
         double df = -m_Parameters.m_SignalGen.m_Drift/(signal.Size()*signal.Size());
         signal[i] += signal[i]*df*i*i;
       }
     }
 
     if (m_Parameters.m_NoiseModel)
       m_Parameters.m_NoiseModel->AddNoise(signal);
 
     for (unsigned int i=0; i<signal.Size(); i++)
     {
       if (signal[i]>0)
         signal[i] = floor(signal[i]+0.5);
       else
         signal[i] = ceil(signal[i]-0.5);
     }
     it4.Set(signal);
     ++it4;
   }
   this->SetNthOutput(0, m_OutputImage);
 
   PrintToLog("\n", false);
   PrintToLog("Finished simulation");
   m_TimeProbe.Stop();
 
   if (m_Parameters.m_SignalGen.m_DoAddMotion)
   {
     PrintToLog("\nHead motion log:", false);
     PrintToLog(m_MotionLog, false, false);
   }
 
   if (m_Parameters.m_Misc.m_DoAddSpikes && m_Parameters.m_SignalGen.m_Spikes>0)
   {
     PrintToLog("\nSpike log:", false);
     PrintToLog(m_SpikeLog, false, false);
   }
 
   if (m_Logfile.is_open())
     m_Logfile.close();
 }
 
 
 template< class PixelType >
 void TractsToDWIImageFilter< PixelType >::PrintToLog(std::string m, bool addTime, bool linebreak, bool stdOut)
 {
   // timestamp
   if (addTime)
   {
     if ( m_Logfile.is_open() )
       m_Logfile << this->GetTime() << " > ";
     m_StatusText += this->GetTime() + " > ";
     if (stdOut)
       std::cout << this->GetTime() << " > ";
   }
 
   // message
   if (m_Logfile.is_open())
     m_Logfile << m;
   m_StatusText += m;
   if (stdOut)
     std::cout << m;
 
   // new line
   if (linebreak)
   {
     if (m_Logfile.is_open())
       m_Logfile << "\n";
     m_StatusText += "\n";
     if (stdOut)
       std::cout << "\n";
   }
   if ( m_Logfile.is_open() )
     m_Logfile.flush();
 }
 
 template< class PixelType >
 void TractsToDWIImageFilter< PixelType >::SimulateMotion(int g)
 {
   if ( m_Parameters.m_SignalGen.m_DoAddMotion &&
        m_Parameters.m_SignalGen.m_DoRandomizeMotion &&
        g>0 &&
        m_Parameters.m_SignalGen.m_MotionVolumes[g-1])
   {
     // The last volume was randomly moved, so we have to reset to fiberbundle and the mask.
     // Without motion or with linear motion, we keep the last position --> no reset.
     m_FiberBundleTransformed = m_FiberBundle->GetDeepCopy();
 
     if (m_MaskImageSet)
     {
       auto duplicator = itk::ImageDuplicator<ItkUcharImgType>::New();
       duplicator->SetInputImage(m_Parameters.m_SignalGen.m_MaskImage);
       duplicator->Update();
       m_TransformedMaskImage = duplicator->GetOutput();
     }
   }
 
   VectorType rotation;
   VectorType translation;
 
   // is motion artifact enabled?
   // is the current volume g affected by motion?
   if ( m_Parameters.m_SignalGen.m_DoAddMotion
        && m_Parameters.m_SignalGen.m_MotionVolumes[g]
        && g<static_cast<int>(m_Parameters.m_SignalGen.GetNumVolumes()) )
   {
     // adjust motion transforms
     if ( m_Parameters.m_SignalGen.m_DoRandomizeMotion )
     {
       // randomly
       rotation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[0]*2)
           -m_Parameters.m_SignalGen.m_Rotation[0];
       rotation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[1]*2)
           -m_Parameters.m_SignalGen.m_Rotation[1];
       rotation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[2]*2)
           -m_Parameters.m_SignalGen.m_Rotation[2];
 
       translation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[0]*2)
           -m_Parameters.m_SignalGen.m_Translation[0];
       translation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[1]*2)
           -m_Parameters.m_SignalGen.m_Translation[1];
       translation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[2]*2)
           -m_Parameters.m_SignalGen.m_Translation[2];
 
       m_FiberBundleTransformed->TransformFibers(rotation[0], rotation[1], rotation[2], translation[0], translation[1], translation[2]);
 
     }
     else
     {
       // linearly
       rotation = m_Parameters.m_SignalGen.m_Rotation / m_NumMotionVolumes;
       translation = m_Parameters.m_SignalGen.m_Translation / m_NumMotionVolumes;
       m_MotionCounter++;
 
       m_FiberBundleTransformed->TransformFibers(rotation[0], rotation[1], rotation[2], translation[0], translation[1], translation[2]);
 
       rotation *= m_MotionCounter;
       translation *= m_MotionCounter;
     }
     MatrixType rotationMatrix = mitk::imv::GetRotationMatrixItk(rotation[0], rotation[1], rotation[2]);
     MatrixType rotationMatrixInv = mitk::imv::GetRotationMatrixItk(-rotation[0], -rotation[1], -rotation[2]);
 
     m_Rotations.push_back(rotationMatrix);
     m_RotationsInv.push_back(rotationMatrixInv);
     m_Translations.push_back(translation);
 
     // move mask image accoring to new transform
     if (m_MaskImageSet)
     {
       ImageRegionIterator<ItkUcharImgType> maskIt(m_UpsampledMaskImage, m_UpsampledMaskImage->GetLargestPossibleRegion());
       m_TransformedMaskImage->FillBuffer(0);
 
       while(!maskIt.IsAtEnd())
       {
         if (maskIt.Get()<=0)
         {
           ++maskIt;
           continue;
         }
 
         DoubleDwiType::IndexType index = maskIt.GetIndex();
         m_TransformedMaskImage->TransformPhysicalPointToIndex(GetMovedPoint(index, true), index);
 
         if (m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(index))
           m_TransformedMaskImage->SetPixel(index, 100);
         ++maskIt;
       }
     }
   }
   else
   {
     if (m_Parameters.m_SignalGen.m_DoAddMotion && !m_Parameters.m_SignalGen.m_DoRandomizeMotion && g>0)
     {
       rotation = m_Parameters.m_SignalGen.m_Rotation / m_NumMotionVolumes;
       rotation *= m_MotionCounter;
       m_Rotations.push_back(m_Rotations.back());
       m_RotationsInv.push_back(m_RotationsInv.back());
       m_Translations.push_back(m_Translations.back());
     }
     else
     {
       rotation.Fill(0.0);
       VectorType translation; translation.Fill(0.0);
       MatrixType rotation_matrix; rotation_matrix.SetIdentity();
 
       m_Rotations.push_back(rotation_matrix);
       m_RotationsInv.push_back(rotation_matrix);
       m_Translations.push_back(translation);
     }
   }
 
   if (m_Parameters.m_SignalGen.m_DoAddMotion)
   {
     m_MotionLog += boost::lexical_cast<std::string>(g) + " rotation: " + boost::lexical_cast<std::string>(rotation[0])
         + "," + boost::lexical_cast<std::string>(rotation[1])
         + "," + boost::lexical_cast<std::string>(rotation[2]) + ";";
 
     m_MotionLog += " translation: " + boost::lexical_cast<std::string>(m_Translations.back()[0])
         + "," + boost::lexical_cast<std::string>(m_Translations.back()[1])
         + "," + boost::lexical_cast<std::string>(m_Translations.back()[2]) + "\n";
   }
 }
 
 template< class PixelType >
 itk::Point<float, 3> TractsToDWIImageFilter< PixelType >::GetMovedPoint(itk::Index<3>& index, bool forward)
 {
   itk::Point<float, 3> transformed_point;
   float tx = m_Translations.back()[0];
   float ty = m_Translations.back()[1];
   float tz = m_Translations.back()[2];
 
   if (forward)
   {
     m_UpsampledMaskImage->TransformIndexToPhysicalPoint(index, transformed_point);
     m_FiberBundle->TransformPoint<>(transformed_point, m_Rotations.back(), tx, ty, tz);
   }
   else
   {
     tx *= -1; ty *= -1; tz *= -1;
     m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, transformed_point);
     m_FiberBundle->TransformPoint<>(transformed_point, m_RotationsInv.back(), tx, ty, tz);
   }
   return transformed_point;
 }
 
 template< class PixelType >
 void TractsToDWIImageFilter< PixelType >::
 SimulateExtraAxonalSignal(ItkUcharImgType::IndexType& index, itk::Point<float, 3>& volume_fraction_point, double intraAxonalVolume, int g)
 {
   int numFiberCompartments = m_Parameters.m_FiberModelList.size();
   int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size();
 
   if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume)
   {
     // simulate signal for largest non-fiber compartment
     int max_compartment_index = 0;
     double max_fraction = 0;
     if (numNonFiberCompartments>1)
     {
       for (int i=0; i<numNonFiberCompartments; ++i)
       {
         m_DoubleInterpolator->SetInputImage(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage());
         double compartment_fraction = mitk::imv::GetImageValue<double>(volume_fraction_point, true, m_DoubleInterpolator);
         if (compartment_fraction<0)
           mitkThrow() << "Volume fraction image (index " << i << ") contains values less than zero!";
 
         if (compartment_fraction>max_fraction)
         {
           max_fraction = compartment_fraction;
           max_compartment_index = i;
         }
       }
     }
 
     DoubleDwiType::Pointer doubleDwi = m_CompartmentImages.at(max_compartment_index+numFiberCompartments);
     DoubleDwiType::PixelType pix = doubleDwi->GetPixel(index);
     pix[g] += m_Parameters.m_NonFiberModelList[max_compartment_index]->SimulateMeasurement(g, m_NullDir)*m_VoxelVolume;
     doubleDwi->SetPixel(index, pix);
 
     if (g==0)
       m_VolumeFractions.at(max_compartment_index+numFiberCompartments)->SetPixel(index, 1);
   }
   else
   {
     std::vector< double > fractions;
     if (g==0)
       m_VolumeFractions.at(0)->SetPixel(index, intraAxonalVolume/m_VoxelVolume);
 
     double extraAxonalVolume = m_VoxelVolume-intraAxonalVolume;    // non-fiber volume
     if (extraAxonalVolume<0)
     {
       if (extraAxonalVolume<-0.001)
         MITK_ERROR << "Corrupted intra-axonal signal voxel detected. Fiber volume larger voxel volume! " << m_VoxelVolume << "<" << intraAxonalVolume;
       extraAxonalVolume = 0;
     }
     double interAxonalVolume = 0;
     if (numFiberCompartments>1)
       interAxonalVolume = extraAxonalVolume * intraAxonalVolume/m_VoxelVolume;   // inter-axonal fraction of non fiber compartment
     double nonFiberVolume = extraAxonalVolume - interAxonalVolume;        // rest of compartment
     if (nonFiberVolume<0)
     {
       if (nonFiberVolume<-0.001)
         MITK_ERROR << "Corrupted signal voxel detected. Fiber volume larger voxel volume!";
       nonFiberVolume = 0;
       interAxonalVolume = extraAxonalVolume;
     }
 
     double compartmentSum = intraAxonalVolume;
     fractions.push_back(intraAxonalVolume/m_VoxelVolume);
 
     // rescale extra-axonal fiber signal
     for (int i=1; i<numFiberCompartments; ++i)
     {
       if (m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage()!=nullptr)
       {
         m_DoubleInterpolator->SetInputImage(m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage());
         interAxonalVolume = mitk::imv::GetImageValue<double>(volume_fraction_point, true, m_DoubleInterpolator)*m_VoxelVolume;
         if (interAxonalVolume<0)
           mitkThrow() << "Volume fraction image (index " << i+1 << ") contains negative values!";
       }
 
       DoubleDwiType::PixelType pix = m_CompartmentImages.at(i)->GetPixel(index);
       pix[g] *= interAxonalVolume;
       m_CompartmentImages.at(i)->SetPixel(index, pix);
 
       compartmentSum += interAxonalVolume;
       fractions.push_back(interAxonalVolume/m_VoxelVolume);
       if (g==0)
         m_VolumeFractions.at(i)->SetPixel(index, interAxonalVolume/m_VoxelVolume);
     }
 
     for (int i=0; i<numNonFiberCompartments; ++i)
     {
       double volume = nonFiberVolume;
       if (m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()!=nullptr)
       {
         m_DoubleInterpolator->SetInputImage(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage());
         volume = mitk::imv::GetImageValue<double>(volume_fraction_point, true, m_DoubleInterpolator)*m_VoxelVolume;
         if (volume<0)
           mitkThrow() << "Volume fraction image (index " << numFiberCompartments+i+1 << ") contains negative values (non-fiber compartment)!";
 
         if (m_UseRelativeNonFiberVolumeFractions)
           volume *= nonFiberVolume/m_VoxelVolume;
       }
 
       DoubleDwiType::PixelType pix = m_CompartmentImages.at(i+numFiberCompartments)->GetPixel(index);
       pix[g] += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement(g, m_NullDir)*volume;
       m_CompartmentImages.at(i+numFiberCompartments)->SetPixel(index, pix);
 
       compartmentSum += volume;
       fractions.push_back(volume/m_VoxelVolume);
       if (g==0)
         m_VolumeFractions.at(i+numFiberCompartments)->SetPixel(index, volume/m_VoxelVolume);
     }
 
     if (compartmentSum/m_VoxelVolume>1.05)
     {
       MITK_ERROR << "Compartments do not sum to 1 in voxel " << index << " (" << compartmentSum/m_VoxelVolume << ")";
       for (auto val : fractions)
         MITK_ERROR << val;
     }
   }
 }
 
 template< class PixelType >
 itk::Vector<double, 3> TractsToDWIImageFilter< PixelType >::GetItkVector(double point[3])
 {
   itk::Vector<double, 3> itkVector;
   itkVector[0] = point[0];
   itkVector[1] = point[1];
   itkVector[2] = point[2];
   return itkVector;
 }
 
 template< class PixelType >
 vnl_vector_fixed<double, 3> TractsToDWIImageFilter< PixelType >::GetVnlVector(double point[3])
 {
   vnl_vector_fixed<double, 3> vnlVector;
   vnlVector[0] = point[0];
   vnlVector[1] = point[1];
   vnlVector[2] = point[2];
   return vnlVector;
 }
 
 template< class PixelType >
 vnl_vector_fixed<double, 3> TractsToDWIImageFilter< PixelType >::GetVnlVector(Vector<float,3>& vector)
 {
   vnl_vector_fixed<double, 3> vnlVector;
   vnlVector[0] = vector[0];
   vnlVector[1] = vector[1];
   vnlVector[2] = vector[2];
   return vnlVector;
 }
 
 template< class PixelType >
 double TractsToDWIImageFilter< PixelType >::RoundToNearest(double num) {
   return (num > 0.0) ? floor(num + 0.5) : ceil(num - 0.5);
 }
 
 template< class PixelType >
 std::string TractsToDWIImageFilter< PixelType >::GetTime()
 {
   m_TimeProbe.Stop();
   unsigned long total = RoundToNearest(m_TimeProbe.GetTotal());
   unsigned long hours = total/3600;
   unsigned long minutes = (total%3600)/60;
   unsigned long seconds = total%60;
   std::string out = "";
   out.append(boost::lexical_cast<std::string>(hours));
   out.append(":");
   out.append(boost::lexical_cast<std::string>(minutes));
   out.append(":");
   out.append(boost::lexical_cast<std::string>(seconds));
   m_TimeProbe.Start();
   return out;
 }
 
 }
diff --git a/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp
index e94fee0b22..4eab432c56 100644
--- a/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp
+++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp
@@ -1,1143 +1,1143 @@
 
 /*===================================================================
 
 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 <mitkTestFixture.h>
 #include <mitkTestingMacros.h>
 #include <mitkPlanarPolygon.h>
 
 #include <mitkIOUtil.h>
 
 #include <mitkPlanarFigureMaskGenerator.h>
 #include <mitkImageMaskGenerator.h>
 #include <mitkImageStatisticsConstants.h>
 
 /**
  * \brief Test class for mitkImageStatisticsCalculator
  *
  * This test covers:
  * - instantiation of an ImageStatisticsCalculator class
  * - correctness of statistics when using PlanarFigures for masking
  */
 class mitkImageStatisticsCalculatorTestSuite : public mitk::TestFixture
 {
   CPPUNIT_TEST_SUITE(mitkImageStatisticsCalculatorTestSuite);
   MITK_TEST(TestUninitializedImage);
   MITK_TEST(TestCase1);
   MITK_TEST(TestCase2);
   MITK_TEST(TestCase3);
   MITK_TEST(TestCase4);
   MITK_TEST(TestCase5);
   MITK_TEST(TestCase6);
   MITK_TEST(TestCase7);
   MITK_TEST(TestCase8);
   MITK_TEST(TestCase9);
   MITK_TEST(TestCase10);
   MITK_TEST(TestCase11);
   MITK_TEST(TestCase12);
   MITK_TEST(TestPic3DCroppedNoMask);
   MITK_TEST(TestPic3DCroppedBinMask);
   MITK_TEST(TestPic3DCroppedMultilabelMask);
   MITK_TEST(TestPic3DCroppedPlanarFigure);
   MITK_TEST(TestUS4DCroppedNoMaskTimeStep1);
   MITK_TEST(TestUS4DCroppedBinMaskTimeStep1);
   MITK_TEST(TestUS4DCroppedMultilabelMaskTimeStep1);
   MITK_TEST(TestUS4DCroppedPlanarFigureTimeStep1);
   MITK_TEST(TestUS4DCroppedAllTimesteps);
   MITK_TEST(TestUS4DCropped3DMask);
   CPPUNIT_TEST_SUITE_END();
 
 public:
 
   void TestUninitializedImage();
 
   void TestCase1();
   void TestCase2();
   void TestCase3();
   void TestCase4();
   void TestCase5();
   void TestCase6();
   void TestCase7();
   void TestCase8();
   void TestCase9();
   void TestCase10();
   void TestCase11();
   void TestCase12();
   
   void TestPic3DCroppedNoMask();
   void TestPic3DCroppedBinMask();
   void TestPic3DCroppedMultilabelMask();
   void TestPic3DCroppedPlanarFigure();
 
   void TestUS4DCroppedNoMaskTimeStep1();
   void TestUS4DCroppedBinMaskTimeStep1();
   void TestUS4DCroppedMultilabelMaskTimeStep1();
   void TestUS4DCroppedPlanarFigureTimeStep1();
   void TestUS4DCroppedAllTimesteps();
   void TestUS4DCropped3DMask();
 private:
 	mitk::Image::ConstPointer m_TestImage;
 
 	mitk::Image::ConstPointer m_Pic3DCroppedImage;
 	mitk::Image::Pointer m_Pic3DCroppedBinMask;
 	mitk::Image::Pointer m_Pic3DCroppedMultilabelMask;
 	mitk::PlanarFigure::Pointer m_Pic3DCroppedPlanarFigure;
 
 	mitk::Image::ConstPointer m_US4DCroppedImage;
 	mitk::Image::Pointer m_US4DCroppedBinMask;
 	mitk::Image::Pointer m_US4DCroppedMultilabelMask;
 	mitk::Image::Pointer m_US4DCropped3DBinMask;
 	mitk::PlanarFigure::Pointer m_US4DCroppedPlanarFigure;
 
 	mitk::PlaneGeometry::Pointer m_Geometry;
 
 	// creates a polygon given a geometry and a vector of 2d points
 	mitk::PlanarPolygon::Pointer GeneratePlanarPolygon(mitk::PlaneGeometry::Pointer geometry, std::vector <mitk::Point2D> points);
 
 	// universal function to calculate statistics
 	const mitk::ImageStatisticsContainer::Pointer ComputeStatistics(mitk::Image::ConstPointer image,
 		mitk::MaskGenerator::Pointer maskGen = nullptr,
 		mitk::MaskGenerator::Pointer secondardMaskGen = nullptr,
 		unsigned short label = 1);
 
 	void VerifyStatistics(mitk::ImageStatisticsContainer::ImageStatisticsObject stats,
 		mitk::ImageStatisticsContainer::RealType testMean, mitk::ImageStatisticsContainer::RealType testSD, mitk::ImageStatisticsContainer::RealType testMedian = 0);
 
 	// T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy)
 	void VerifyStatistics(mitk::ImageStatisticsContainer::ImageStatisticsObject stats,
 		mitk::ImageStatisticsContainer::VoxelCountType N,
 		mitk::ImageStatisticsContainer::RealType mean,
 		mitk::ImageStatisticsContainer::RealType MPP,
 		mitk::ImageStatisticsContainer::RealType skewness,
 		mitk::ImageStatisticsContainer::RealType kurtosis,
 		mitk::ImageStatisticsContainer::RealType variance,
 		mitk::ImageStatisticsContainer::RealType stdev,
 		mitk::ImageStatisticsContainer::RealType min,
 		mitk::ImageStatisticsContainer::RealType max,
 		mitk::ImageStatisticsContainer::RealType RMS,
 		mitk::ImageStatisticsContainer::IndexType minIndex,
 		mitk::ImageStatisticsContainer::IndexType maxIndex);
 };
 
 void mitkImageStatisticsCalculatorTestSuite::TestUninitializedImage()
 {
 	/*****************************
 	* loading uninitialized image to datastorage
 	******************************/
 	MITK_INFO << std::endl << "Test uninitialized image: -----------------------------------------------------------------------------------";
 	mitk::Image::Pointer image = mitk::Image::New();
 	mitk::DataNode::Pointer node = mitk::DataNode::New();
 	node->SetData(image);
 
 	mitk::ImageStatisticsCalculator::Pointer is = mitk::ImageStatisticsCalculator::New();
 	CPPUNIT_ASSERT_THROW(is->GetStatistics(), mitk::Exception);
 }
 
 void mitkImageStatisticsCalculatorTestSuite::TestCase1()
 {
 	/*****************************
 	 * one whole white pixel
 	 * -> mean of 255 expected
 	 ******************************/
 	MITK_INFO << std::endl << "Test case 1:-----------------------------------------------------------------------------------";
 
 	std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd");
 	m_TestImage = mitk::IOUtil::Load<mitk::Image>(filename);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading an mitk::Image", m_TestImage.IsNotNull());
 
 	m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0);
 	CPPUNIT_ASSERT_MESSAGE("Failed getting image geometry", m_Geometry.IsNotNull());
 
 	mitk::Point2D pnt1; pnt1[0] = 10.5; pnt1[1] = 3.5;
 	mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5;
 	mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5;
 	mitk::Point2D pnt4; pnt4[0] = 10.5; pnt4[1] = 4.5;
 	std::vector<mitk::Point2D> points{ pnt1,pnt2,pnt3,pnt4 };
         auto figure = GeneratePlanarPolygon(m_Geometry, points);
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 
 	mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New();
 	planFigMaskGen->SetInputImage(m_TestImage);
         planFigMaskGen->SetPlanarFigure(figure.GetPointer());
 
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer()));
 	auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0);
 
 	this->VerifyStatistics(statisticsObjectTimestep0, 255.0, 0.0, 255.0);
 }
 
 void mitkImageStatisticsCalculatorTestSuite::TestCase2()
 {
 	/*****************************
 	 * half pixel in x-direction (white)
 	 * -> mean of 255 expected
 	 ******************************/
 	MITK_INFO << std::endl << "Test case 2:-----------------------------------------------------------------------------------";
 
 	std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd");
 	m_TestImage = mitk::IOUtil::Load<mitk::Image>(filename);
 	m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0);
 
 	mitk::Point2D pnt1; pnt1[0] = 10.0; pnt1[1] = 3.5;
 	mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5;
 	mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5;
 	mitk::Point2D pnt4; pnt4[0] = 10.0; pnt4[1] = 4.5;
 	std::vector<mitk::Point2D> points{ pnt1,pnt2,pnt3,pnt4 };
         auto figure = GeneratePlanarPolygon(m_Geometry, points);
 
 	mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New();
 	planFigMaskGen->SetInputImage(m_TestImage);
         planFigMaskGen->SetPlanarFigure(figure.GetPointer());
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer()));
 	auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0);
 
 	this->VerifyStatistics(statisticsObjectTimestep0, 255.0, 0.0, 255.0);
 }
 
 void mitkImageStatisticsCalculatorTestSuite::TestCase3()
 {
 	/*****************************
 	 * half pixel in diagonal-direction (white)
 	 * -> mean of 255 expected
 	 ******************************/
 	MITK_INFO << std::endl << "Test case 3:-----------------------------------------------------------------------------------";
 
 	std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd");
 	m_TestImage = mitk::IOUtil::Load<mitk::Image>(filename);
 	m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0);
 
 	mitk::Point2D pnt1; pnt1[0] = 10.5; pnt1[1] = 3.5;
 	mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5;
 	mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5;
 	std::vector<mitk::Point2D> points{ pnt1,pnt2,pnt3 };
         auto figure = GeneratePlanarPolygon(m_Geometry, points);
 
 	mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New();
 	planFigMaskGen->SetInputImage(m_TestImage);
         planFigMaskGen->SetPlanarFigure(figure.GetPointer());
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer()));
 	auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0);
 
 	this->VerifyStatistics(statisticsObjectTimestep0, 255.0, 0.0, 255.0);
 }
 
 void mitkImageStatisticsCalculatorTestSuite::TestCase4()
 {
 	/*****************************
 	 * one pixel (white) + 2 half pixels (white) + 1 half pixel (black)
 	 * -> mean of 191.25 expected
 	 ******************************/
 	MITK_INFO << std::endl << "Test case 4:-----------------------------------------------------------------------------------";
 
 	std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd");
 	m_TestImage = mitk::IOUtil::Load<mitk::Image>(filename);
 	m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0);
 
 	mitk::Point2D pnt1; pnt1[0] = 1.1; pnt1[1] = 1.1;
 	mitk::Point2D pnt2; pnt2[0] = 2.0; pnt2[1] = 2.0;
 	mitk::Point2D pnt3; pnt3[0] = 3.0; pnt3[1] = 1.0;
 	mitk::Point2D pnt4; pnt4[0] = 2.0; pnt4[1] = 0.0;
 	std::vector<mitk::Point2D> points{ pnt1,pnt2,pnt3,pnt4 };
         auto figure = GeneratePlanarPolygon(m_Geometry, points);
 
 	mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New();
 	planFigMaskGen->SetInputImage(m_TestImage);
         planFigMaskGen->SetPlanarFigure(figure.GetPointer());
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer()));
 	auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0);
 
 	this->VerifyStatistics(statisticsObjectTimestep0, 191.25, 110.41823898251593, 253.72499847412109);
 }
 
 void mitkImageStatisticsCalculatorTestSuite::TestCase5()
 {
 	/*****************************
 	 * whole pixel (white) + half pixel (gray) in x-direction
 	 * -> mean of 191.5 expected
 	 ******************************/
 	MITK_INFO << std::endl << "Test case 5:-----------------------------------------------------------------------------------";
 
 	std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd");
 	m_TestImage = mitk::IOUtil::Load<mitk::Image>(filename);
 	m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0);
 
 	mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5;
 	mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5;
 	mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5;
 	mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.5;
 	std::vector<mitk::Point2D> points{ pnt1,pnt2,pnt3,pnt4 };
         auto figure = GeneratePlanarPolygon(m_Geometry, points);
 
 	mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New();
 	planFigMaskGen->SetInputImage(m_TestImage);
         planFigMaskGen->SetPlanarFigure(figure.GetPointer());
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer()));
 	auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0);
 
 	this->VerifyStatistics(statisticsObjectTimestep0, 191.50, 63.50, 128.63499999046327);
 }
 
 void mitkImageStatisticsCalculatorTestSuite::TestCase6()
 {
 	/*****************************
 	 * quarter pixel (black) + whole pixel (white) + half pixel (gray) in x-direction
 	 * -> mean of 191.5 expected
 	 ******************************/
 	MITK_INFO << std::endl << "Test case 6:-----------------------------------------------------------------------------------";
 
 	std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd");
 	m_TestImage = mitk::IOUtil::Load<mitk::Image>(filename);
 	m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0);
 
 	mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5;
 	mitk::Point2D pnt2; pnt2[0] = 9.25; pnt2[1] = 3.5;
 	mitk::Point2D pnt3; pnt3[0] = 9.25; pnt3[1] = 4.5;
 	mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.5;
 	std::vector<mitk::Point2D> points{ pnt1,pnt2,pnt3,pnt4 };
         auto figure = GeneratePlanarPolygon(m_Geometry, points);
 
 	mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New();
 	planFigMaskGen->SetInputImage(m_TestImage);
         planFigMaskGen->SetPlanarFigure(figure.GetPointer());
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer()));
 	auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0);
 
 	this->VerifyStatistics(statisticsObjectTimestep0, 191.5, 63.50, 128.63499999046327);
 }
 
 void mitkImageStatisticsCalculatorTestSuite::TestCase7()
 {
 	/*****************************
 	 * half pixel (black) + whole pixel (white) + half pixel (gray) in x-direction
 	 * -> mean of 127.66 expected
 	 ******************************/
 	MITK_INFO << std::endl << "Test case 7:-----------------------------------------------------------------------------------";
 
 	std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd");
 	m_TestImage = mitk::IOUtil::Load<mitk::Image>(filename);
 	m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0);
 
 	mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5;
 	mitk::Point2D pnt2; pnt2[0] = 9.0; pnt2[1] = 3.5;
 	mitk::Point2D pnt3; pnt3[0] = 9.0; pnt3[1] = 4.0;
 	mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.0;
 	std::vector<mitk::Point2D> points{ pnt1,pnt2,pnt3,pnt4 };
         auto figure = GeneratePlanarPolygon(m_Geometry, points);
 
 	mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New();
 	planFigMaskGen->SetInputImage(m_TestImage);
         planFigMaskGen->SetPlanarFigure(figure.GetPointer());
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer()));
 	auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0);
 
 	this->VerifyStatistics(statisticsObjectTimestep0, 127.666666666666667, 104.10358089689113, 128.7750015258789);
 }
 
 void mitkImageStatisticsCalculatorTestSuite::TestCase8()
 {
 	/*****************************
 	 * whole pixel (gray)
 	 * -> mean of 128 expected
 	 ******************************/
 	MITK_INFO << std::endl << "Test case 8:-----------------------------------------------------------------------------------";
 
 	std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd");
 	m_TestImage = mitk::IOUtil::Load<mitk::Image>(filename);
 	m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0);
 
 	mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5;
 	mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 11.5;
 	mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 11.5;
 	mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5;
 	std::vector<mitk::Point2D> points{ pnt1,pnt2,pnt3,pnt4 };
         auto figure = GeneratePlanarPolygon(m_Geometry, points);
 
 	mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New();
 	planFigMaskGen->SetInputImage(m_TestImage);
         planFigMaskGen->SetPlanarFigure(figure.GetPointer());
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer()));
 	auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0);
 
 	this->VerifyStatistics(statisticsObjectTimestep0, 128.0, 0.0, 128.0);
 }
 
 void mitkImageStatisticsCalculatorTestSuite::TestCase9()
 {
 	/*****************************
 	 * whole pixel (gray) + half pixel (white) in y-direction
 	 * -> mean of 191.5 expected
 	 ******************************/
 	MITK_INFO << std::endl << "Test case 9:-----------------------------------------------------------------------------------";
 
 	std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd");
 	m_TestImage = mitk::IOUtil::Load<mitk::Image>(filename);
 	m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0);
 
 	mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5;
 	mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 12.0;
 	mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 12.0;
 	mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5;
 	std::vector<mitk::Point2D> points{ pnt1,pnt2,pnt3,pnt4 };
         auto figure = GeneratePlanarPolygon(m_Geometry, points);
 
 	mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New();
 	planFigMaskGen->SetInputImage(m_TestImage);
         planFigMaskGen->SetPlanarFigure(figure.GetPointer());
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer()));
 	auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0);
 
 	this->VerifyStatistics(statisticsObjectTimestep0, 191.5, 63.50, 128.63499999046327);
 }
 
 void mitkImageStatisticsCalculatorTestSuite::TestCase10()
 {
 	/*****************************
 	 * 2 whole pixel (white) + 2 whole pixel (black) in y-direction
 	 * -> mean of 127.66 expected
 	 ******************************/
 	MITK_INFO << std::endl << "Test case 10:-----------------------------------------------------------------------------------";
 
 	std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd");
 	m_TestImage = mitk::IOUtil::Load<mitk::Image>(filename);
 	m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0);
 
 	mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5;
 	mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 13.5;
 	mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 13.5;
 	mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5;
 	std::vector<mitk::Point2D> points{ pnt1,pnt2,pnt3,pnt4 };
         auto figure = GeneratePlanarPolygon(m_Geometry, points);
 
 	mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New();
 	planFigMaskGen->SetInputImage(m_TestImage);
         planFigMaskGen->SetPlanarFigure(figure.GetPointer());
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer()));
 	auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0);
 
 	this->VerifyStatistics(statisticsObjectTimestep0, 127.666666666666667, 104.10358089689113, 128.7750015258789);
 }
 
 void mitkImageStatisticsCalculatorTestSuite::TestCase11()
 {
 	/*****************************
 	 * 9 whole pixels (white) + 3 half pixels (white)
 	 * + 3 whole pixel (black) [ + 3 slightly less than half pixels (black)]
 	 * -> mean of 204.0 expected
 	 ******************************/
 	MITK_INFO << std::endl << "Test case 11:-----------------------------------------------------------------------------------";
 	std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd");
 	m_TestImage = mitk::IOUtil::Load<mitk::Image>(filename);
 	m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0);
 
 	mitk::Point2D pnt1; pnt1[0] = 0.5; pnt1[1] = 0.5;
 	mitk::Point2D pnt2; pnt2[0] = 3.5; pnt2[1] = 3.5;
 	mitk::Point2D pnt3; pnt3[0] = 8.4999; pnt3[1] = 3.5;
 	mitk::Point2D pnt4; pnt4[0] = 5.4999; pnt4[1] = 0.5;
 	std::vector<mitk::Point2D> points{ pnt1,pnt2,pnt3,pnt4 };
         auto figure = GeneratePlanarPolygon(m_Geometry, points);
 
 	mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New();
 	planFigMaskGen->SetInputImage(m_TestImage);
         planFigMaskGen->SetPlanarFigure(figure.GetPointer());
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer()));
 	auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0);
 
 	this->VerifyStatistics(statisticsObjectTimestep0, 204.0, 102.00, 253.724998474121083);
 }
 
 void mitkImageStatisticsCalculatorTestSuite::TestCase12()
 {
 	/*****************************
 	 * half pixel (white) + whole pixel (white) + half pixel (black)
 	 * -> mean of 212.66 expected
 	 ******************************/
 	MITK_INFO << std::endl << "Test case 12:-----------------------------------------------------------------------------------";
 
 	std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd");
 	m_TestImage = mitk::IOUtil::Load<mitk::Image>(filename);
 	m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0);
 
 	mitk::Point2D pnt1; pnt1[0] = 9.5; pnt1[1] = 0.5;
 	mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 2.5;
 	mitk::Point2D pnt3; pnt3[0] = 11.5; pnt3[1] = 2.5;
 	std::vector<mitk::Point2D> points{ pnt1,pnt2,pnt3 };
         auto figure = GeneratePlanarPolygon(m_Geometry, points);
 
 	mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New();
 	planFigMaskGen->SetInputImage(m_TestImage);
         planFigMaskGen->SetPlanarFigure(figure.GetPointer());
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer()));
 	auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0);
 
 	this->VerifyStatistics(statisticsObjectTimestep0, 212.666666666666667, 59.8683741404609923, 254.36499786376954);
 }
 
 // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy)
 void mitkImageStatisticsCalculatorTestSuite::TestPic3DCroppedNoMask()
 {
 	MITK_INFO << std::endl << "Test Pic3D cropped without mask:-----------------------------------------------------------------------------------";
 
 	std::string Pic3DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3D_cropped.nrrd");
 	m_Pic3DCroppedImage = mitk::IOUtil::Load<mitk::Image>(Pic3DCroppedFile);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading Pic3D_cropped", m_Pic3DCroppedImage.IsNotNull());
 	//calculated ground truth via script
 	mitk::ImageStatisticsContainer::VoxelCountType expected_N = 27;
 	mitk::ImageStatisticsContainer::RealType expected_mean = -564.1481481481481481;
 	mitk::ImageStatisticsContainer::RealType expected_MPP = 113.66666666666667;
         //mitk::ImageStatisticsContainer::RealType expected_median = -825;
 	mitk::ImageStatisticsContainer::RealType expected_skewness = 0.7120461106763573;
 	mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1.8794464383714844;
 	mitk::ImageStatisticsContainer::RealType expected_variance = 140541.38545953357;
 	mitk::ImageStatisticsContainer::RealType expected_standarddev = 374.88849736892911;
 	mitk::ImageStatisticsContainer::RealType expected_min = -927;
 	mitk::ImageStatisticsContainer::RealType expected_max = 147;
 	mitk::ImageStatisticsContainer::RealType expected_RMS = 677.35110431630551;
 	mitk::ImageStatisticsContainer::IndexType expected_minIndex;
 	expected_minIndex.set_size(3);
 	expected_minIndex[0] = 2;
 	expected_minIndex[1] = 1;
 	expected_minIndex[2] = 1;
 
 	mitk::ImageStatisticsContainer::IndexType expected_maxIndex;
 	expected_maxIndex.set_size(3);
 	expected_maxIndex[0] = 0;
 	expected_maxIndex[1] = 1;
 	expected_maxIndex[2] = 2;
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_Pic3DCroppedImage));
 	auto statisticsObject = statisticsContainer->GetStatisticsForTimeStep(0);
 
 	VerifyStatistics(statisticsObject,
 		expected_N,
 		expected_mean,
 		expected_MPP,
 		expected_skewness,
 		expected_kurtosis,
 		expected_variance,
 		expected_standarddev,
 		expected_min,
 		expected_max,
 		expected_RMS,
 		expected_minIndex,
 		expected_maxIndex);
 }
 
 // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy)
 void mitkImageStatisticsCalculatorTestSuite::TestPic3DCroppedBinMask()
 {
 	MITK_INFO << std::endl << "Test Pic3D cropped binary mask:-----------------------------------------------------------------------------------";
 
 	std::string Pic3DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3D_cropped.nrrd");
 	m_Pic3DCroppedImage = mitk::IOUtil::Load<mitk::Image>(Pic3DCroppedFile);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading Pic3D_cropped", m_Pic3DCroppedImage.IsNotNull());
 
 	std::string Pic3DCroppedBinMaskFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3D_croppedBinMask.nrrd");
 	m_Pic3DCroppedBinMask = mitk::IOUtil::Load<mitk::Image>(Pic3DCroppedBinMaskFile);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading Pic3D binary mask", m_Pic3DCroppedBinMask.IsNotNull());
 	//calculated ground truth via script
 	mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1.0765697398089618;
 	mitk::ImageStatisticsContainer::RealType expected_MPP = -nan("");
 	mitk::ImageStatisticsContainer::RealType expected_max = -22;
 	mitk::ImageStatisticsContainer::RealType expected_mean = -464;
 	mitk::ImageStatisticsContainer::RealType expected_min = -846;
 	mitk::ImageStatisticsContainer::VoxelCountType expected_N = 4;
 	mitk::ImageStatisticsContainer::RealType expected_RMS = 595.42631785973322;
 	mitk::ImageStatisticsContainer::RealType expected_skewness = 0.0544059290851858;
 	mitk::ImageStatisticsContainer::RealType expected_standarddev = 373.14407405183323;
 	mitk::ImageStatisticsContainer::RealType expected_variance = 139236.50;
 	mitk::ImageStatisticsContainer::IndexType expected_minIndex;
 	expected_minIndex.set_size(3);
 	expected_minIndex[0] = 1;
 	expected_minIndex[1] = 0;
 	expected_minIndex[2] = 0;
 
 	mitk::ImageStatisticsContainer::IndexType expected_maxIndex;
 	expected_maxIndex.set_size(3);
 	expected_maxIndex[0] = 0;
 	expected_maxIndex[1] = 0;
 	expected_maxIndex[2] = 1;
 
 	mitk::ImageMaskGenerator::Pointer imgMaskGen = mitk::ImageMaskGenerator::New();
 	imgMaskGen->SetImageMask(m_Pic3DCroppedBinMask);
 	imgMaskGen->SetInputImage(m_Pic3DCroppedImage);
 	imgMaskGen->SetTimeStep(0);
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_Pic3DCroppedImage, imgMaskGen.GetPointer(), nullptr, 1));
 	auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0);
 
 	VerifyStatistics(statisticsObjectTimestep0,
 		expected_N,
 		expected_mean,
 		expected_MPP,
 		expected_skewness,
 		expected_kurtosis,
 		expected_variance,
 		expected_standarddev,
 		expected_min,
 		expected_max,
 		expected_RMS,
 		expected_minIndex,
 		expected_maxIndex);
 }
 
 // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy)
 void mitkImageStatisticsCalculatorTestSuite::TestPic3DCroppedMultilabelMask()
 {
 	MITK_INFO << std::endl << "Test Pic3D cropped multilabel mask:-----------------------------------------------------------------------------------";
 
 	std::string Pic3DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3D_cropped.nrrd");
 	m_Pic3DCroppedImage = mitk::IOUtil::Load<mitk::Image>(Pic3DCroppedFile);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading Pic3D_cropped", m_Pic3DCroppedImage.IsNotNull());
 
 	std::string Pic3DCroppedMultilabelMaskFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3D_croppedMultilabelMask.nrrd");
 	m_Pic3DCroppedMultilabelMask = mitk::IOUtil::Load<mitk::Image>(Pic3DCroppedMultilabelMaskFile);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading Pic3D multilabel mask", m_Pic3DCroppedMultilabelMask.IsNotNull());
 	//calculated ground truth via script
 	mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1.5;
 	mitk::ImageStatisticsContainer::RealType expected_MPP = -nan("");
 	mitk::ImageStatisticsContainer::RealType expected_max = -22;
 	mitk::ImageStatisticsContainer::RealType expected_mean = -586.33333333333333;
 	mitk::ImageStatisticsContainer::RealType expected_min = -916;
 	mitk::ImageStatisticsContainer::VoxelCountType expected_N = 3;
 	mitk::ImageStatisticsContainer::RealType expected_RMS = 710.3006405741163;
 	mitk::ImageStatisticsContainer::RealType expected_skewness = 0.6774469597523700;
 	mitk::ImageStatisticsContainer::RealType expected_standarddev = 400.92421007245525;
 	mitk::ImageStatisticsContainer::RealType expected_variance = 160740.22222222222;
 	mitk::ImageStatisticsContainer::IndexType expected_minIndex;
 	expected_minIndex.set_size(3);
 	expected_minIndex[0] = 2;
 	expected_minIndex[1] = 0;
 	expected_minIndex[2] = 1;
 
 	mitk::ImageStatisticsContainer::IndexType expected_maxIndex;
 	expected_maxIndex.set_size(3);
 	expected_maxIndex[0] = 0;
 	expected_maxIndex[1] = 0;
 	expected_maxIndex[2] = 1;
 
 	mitk::ImageMaskGenerator::Pointer imgMaskGen = mitk::ImageMaskGenerator::New();
 	imgMaskGen->SetImageMask(m_Pic3DCroppedMultilabelMask);
 	imgMaskGen->SetInputImage(m_Pic3DCroppedImage);
 	imgMaskGen->SetTimeStep(0);
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_Pic3DCroppedImage, imgMaskGen.GetPointer(), nullptr, 2));
 	auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0);
 
 	VerifyStatistics(statisticsObjectTimestep0,
 		expected_N,
 		expected_mean,
 		expected_MPP,
 		expected_skewness,
 		expected_kurtosis,
 		expected_variance,
 		expected_standarddev,
 		expected_min,
 		expected_max,
 		expected_RMS,
 		expected_minIndex,
 		expected_maxIndex);
 }
 
 // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy)
 void mitkImageStatisticsCalculatorTestSuite::TestPic3DCroppedPlanarFigure()
 {
 	MITK_INFO << std::endl << "Test Pic3D cropped planar figure:-----------------------------------------------------------------------------------";
 
 	std::string Pic3DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3D_cropped.nrrd");
 	m_Pic3DCroppedImage = mitk::IOUtil::Load<mitk::Image>(Pic3DCroppedFile);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading Pic3D_cropped", m_Pic3DCroppedImage.IsNotNull());
 
 	std::string Pic3DCroppedPlanarFigureFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3D_croppedPF.pf");
 	m_Pic3DCroppedPlanarFigure = mitk::IOUtil::Load<mitk::PlanarFigure>(Pic3DCroppedPlanarFigureFile);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading Pic3D planar figure", m_Pic3DCroppedPlanarFigure.IsNotNull());
 	//calculated ground truth via script
 	mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1;
 	mitk::ImageStatisticsContainer::RealType expected_MPP = -nan("");
 	mitk::ImageStatisticsContainer::RealType expected_max = -67;
 	mitk::ImageStatisticsContainer::RealType expected_mean = -446;
 	mitk::ImageStatisticsContainer::RealType expected_min = -825;
 	mitk::ImageStatisticsContainer::VoxelCountType expected_N = 2;
 	mitk::ImageStatisticsContainer::RealType expected_RMS = 585.28369189650243;
 	mitk::ImageStatisticsContainer::RealType expected_skewness = 0;
 	mitk::ImageStatisticsContainer::RealType expected_standarddev = 379;
 	mitk::ImageStatisticsContainer::RealType expected_variance = 143641;
 	mitk::ImageStatisticsContainer::IndexType expected_minIndex;
 	expected_minIndex.set_size(3);
 	expected_minIndex[0] = 1;
 	expected_minIndex[1] = 1;
 	expected_minIndex[2] = 1;
 
 	mitk::ImageStatisticsContainer::IndexType expected_maxIndex;
 	expected_maxIndex.set_size(3);
 	expected_maxIndex[0] = 0;
 	expected_maxIndex[1] = 1;
 	expected_maxIndex[2] = 1;
 
 	mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New();
 	pfMaskGen->SetInputImage(m_Pic3DCroppedImage);
 	pfMaskGen->SetPlanarFigure(m_Pic3DCroppedPlanarFigure);
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_Pic3DCroppedImage, pfMaskGen.GetPointer()));
 	auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0);
 
 	VerifyStatistics(statisticsObjectTimestep0,
 		expected_N,
 		expected_mean,
 		expected_MPP,
 		expected_skewness,
 		expected_kurtosis,
 		expected_variance,
 		expected_standarddev,
 		expected_min,
 		expected_max,
 		expected_RMS,
 		expected_minIndex,
 		expected_maxIndex);
 }
 
 // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy)
 void mitkImageStatisticsCalculatorTestSuite::TestUS4DCroppedNoMaskTimeStep1()
 {
 	MITK_INFO << std::endl << "Test US4D cropped without mask timestep 1:-----------------------------------------------------------------------------------";
 
         std::string US4DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_cropped.nrrd");
 	m_US4DCroppedImage = mitk::IOUtil::Load<mitk::Image>(US4DCroppedFile);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading US4D_cropped", m_US4DCroppedImage.IsNotNull());
 	//calculated ground truth via script
 	mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1.5398359155908228;
 	mitk::ImageStatisticsContainer::RealType expected_MPP = 157.74074074074073;
 	mitk::ImageStatisticsContainer::RealType expected_max = 199;
 	mitk::ImageStatisticsContainer::RealType expected_mean = 157.74074074074073;
 	mitk::ImageStatisticsContainer::RealType expected_min = 101;
 	mitk::ImageStatisticsContainer::VoxelCountType expected_N = 27;
 	mitk::ImageStatisticsContainer::RealType expected_RMS = 160.991718213494823;
 	mitk::ImageStatisticsContainer::RealType expected_skewness = 0.0347280313508018;
 	mitk::ImageStatisticsContainer::RealType expected_standarddev = 32.189936997387058;
 	mitk::ImageStatisticsContainer::RealType expected_variance = 1036.19204389574722;
 	mitk::ImageStatisticsContainer::IndexType expected_minIndex;
 	expected_minIndex.set_size(3);
 	expected_minIndex[0] = 0;
 	expected_minIndex[1] = 2;
 	expected_minIndex[2] = 0;
 
 	mitk::ImageStatisticsContainer::IndexType expected_maxIndex;
 	expected_maxIndex.set_size(3);
 	expected_maxIndex[0] = 0;
 	expected_maxIndex[1] = 0;
 	expected_maxIndex[2] = 1;
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer=mitk::ImageStatisticsContainer::New();
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_US4DCroppedImage));
 	auto statisticsObjectTimestep1 = statisticsContainer->GetStatisticsForTimeStep(1);
 
 	VerifyStatistics(statisticsObjectTimestep1,
 		expected_N,
 		expected_mean,
 		expected_MPP,
 		expected_skewness,
 		expected_kurtosis,
 		expected_variance,
 		expected_standarddev,
 		expected_min,
 		expected_max,
 		expected_RMS,
 		expected_minIndex,
 		expected_maxIndex);
 }
 
 // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy)
 void mitkImageStatisticsCalculatorTestSuite::TestUS4DCroppedBinMaskTimeStep1()
 {
 	MITK_INFO << std::endl << "Test US4D cropped with binary mask timestep 1:-----------------------------------------------------------------------------------";
 
 	std::string US4DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_cropped.nrrd");
 	m_US4DCroppedImage = mitk::IOUtil::Load<mitk::Image>(US4DCroppedFile);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading US4D_cropped", m_US4DCroppedImage.IsNotNull());
 
 	std::string US4DCroppedBinMaskFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_croppedBinMask.nrrd");
 	m_US4DCroppedBinMask = mitk::IOUtil::Load<mitk::Image>(US4DCroppedBinMaskFile);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading US4D binary mask", m_US4DCroppedBinMask.IsNotNull());
 	//calculated ground truth via script
 	mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1.5863739712889191;
 	mitk::ImageStatisticsContainer::RealType expected_MPP = 166.75;
 	mitk::ImageStatisticsContainer::RealType expected_max = 199;
 	mitk::ImageStatisticsContainer::RealType expected_mean = 166.75;
 	mitk::ImageStatisticsContainer::RealType expected_min = 120;
 	mitk::ImageStatisticsContainer::VoxelCountType expected_N = 4;
 	mitk::ImageStatisticsContainer::RealType expected_RMS = 169.70636405273669;
 	mitk::ImageStatisticsContainer::RealType expected_skewness = -0.4285540263894276;
 	mitk::ImageStatisticsContainer::RealType expected_standarddev = 31.538666744172936;
 	mitk::ImageStatisticsContainer::RealType expected_variance = 994.6874999999999;
 	mitk::ImageStatisticsContainer::IndexType expected_minIndex;
 	expected_minIndex.set_size(3);
 	expected_minIndex[0] = 0;
 	expected_minIndex[1] = 0;
 	expected_minIndex[2] = 2;
 
 	mitk::ImageStatisticsContainer::IndexType expected_maxIndex;
 	expected_maxIndex.set_size(3);
 	expected_maxIndex[0] = 1;
 	expected_maxIndex[1] = 1;
 	expected_maxIndex[2] = 1;
 
 	mitk::ImageMaskGenerator::Pointer imgMask1 = mitk::ImageMaskGenerator::New();
 	imgMask1->SetInputImage(m_US4DCroppedImage);
 	imgMask1->SetImageMask(m_US4DCroppedBinMask);
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer=mitk::ImageStatisticsContainer::New();
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_US4DCroppedImage, imgMask1.GetPointer(), nullptr, 1));
 	auto statisticsObjectTimestep1 = statisticsContainer->GetStatisticsForTimeStep(1);
 
 	VerifyStatistics(statisticsObjectTimestep1,
 		expected_N,
 		expected_mean,
 		expected_MPP,
 		expected_skewness,
 		expected_kurtosis,
 		expected_variance,
 		expected_standarddev,
 		expected_min,
 		expected_max,
 		expected_RMS,
 		expected_minIndex,
 		expected_maxIndex);
 }
 
 // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy)
 void mitkImageStatisticsCalculatorTestSuite::TestUS4DCroppedMultilabelMaskTimeStep1()
 {
 	MITK_INFO << std::endl << "Test US4D cropped with mulitlabel mask timestep 1:-----------------------------------------------------------------------------------";
 
 	std::string US4DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_cropped.nrrd");
 	m_US4DCroppedImage = mitk::IOUtil::Load<mitk::Image>(US4DCroppedFile);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading US4D_cropped", m_US4DCroppedImage.IsNotNull());
 
 	std::string US4DCroppedMultilabelMaskFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_croppedMultilabelMask.nrrd");
 	m_US4DCroppedMultilabelMask = mitk::IOUtil::Load<mitk::Image>(US4DCroppedMultilabelMaskFile);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading US4D multilabel mask", m_US4DCroppedMultilabelMask.IsNotNull());
 	//calculated ground truth via script
 	mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1.0432484564918287;
 	mitk::ImageStatisticsContainer::RealType expected_MPP = 159.75;
 	mitk::ImageStatisticsContainer::RealType expected_max = 199;
 	mitk::ImageStatisticsContainer::RealType expected_mean = 159.75;
 	mitk::ImageStatisticsContainer::RealType expected_min = 120;
 	mitk::ImageStatisticsContainer::VoxelCountType expected_N = 4;
 	mitk::ImageStatisticsContainer::RealType expected_RMS = 163.74446555532802;
 	mitk::ImageStatisticsContainer::RealType expected_skewness = -0.004329226115093;
 	mitk::ImageStatisticsContainer::RealType expected_standarddev = 35.947009611371016;
 	mitk::ImageStatisticsContainer::RealType expected_variance = 1292.187500000000227;
 	mitk::ImageStatisticsContainer::IndexType expected_minIndex;
 	expected_minIndex.set_size(3);
 	expected_minIndex[0] = 0;
 	expected_minIndex[1] = 0;
 	expected_minIndex[2] = 2;
 
 	mitk::ImageStatisticsContainer::IndexType expected_maxIndex;
 	expected_maxIndex.set_size(3);
 	expected_maxIndex[0] = 0;
 	expected_maxIndex[1] = 0;
 	expected_maxIndex[2] = 1;
 
 	mitk::ImageMaskGenerator::Pointer imgMask1 = mitk::ImageMaskGenerator::New();
 	imgMask1->SetInputImage(m_US4DCroppedImage);
 	imgMask1->SetImageMask(m_US4DCroppedMultilabelMask);
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_US4DCroppedImage, imgMask1.GetPointer(), nullptr, 1));
 	auto statisticsObjectTimestep1 = statisticsContainer->GetStatisticsForTimeStep(1);
 
 	VerifyStatistics(statisticsObjectTimestep1,
 		expected_N,
 		expected_mean,
 		expected_MPP,
 		expected_skewness,
 		expected_kurtosis,
 		expected_variance,
 		expected_standarddev,
 		expected_min,
 		expected_max,
 		expected_RMS,
 		expected_minIndex,
 		expected_maxIndex);
 }
 
 // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy)
 void mitkImageStatisticsCalculatorTestSuite::TestUS4DCroppedPlanarFigureTimeStep1()
 {
 	MITK_INFO << std::endl << "Test US4D cropped planar figure timestep 1:-----------------------------------------------------------------------------------";
 
 	std::string US4DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_cropped.nrrd");
 	m_US4DCroppedImage = mitk::IOUtil::Load<mitk::Image>(US4DCroppedFile);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading US4D_cropped", m_US4DCroppedImage.IsNotNull());
 
 	std::string US4DCroppedPlanarFigureFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_croppedPF.pf");
 	m_US4DCroppedPlanarFigure = mitk::IOUtil::Load<mitk::PlanarFigure>(US4DCroppedPlanarFigureFile);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading US4D planar figure", m_US4DCroppedPlanarFigure.IsNotNull());
 	//calculated ground truth via script
 	mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1;
 	mitk::ImageStatisticsContainer::RealType expected_MPP = 172.5;
 	mitk::ImageStatisticsContainer::RealType expected_max = 197;
 	mitk::ImageStatisticsContainer::RealType expected_mean = 172.5;
 	mitk::ImageStatisticsContainer::RealType expected_min = 148;
 	mitk::ImageStatisticsContainer::VoxelCountType expected_N = 2;
 	mitk::ImageStatisticsContainer::RealType expected_RMS = 174.23116827938679;
 	mitk::ImageStatisticsContainer::RealType expected_skewness = 0;
 	mitk::ImageStatisticsContainer::RealType expected_standarddev = 24.5;
 	mitk::ImageStatisticsContainer::RealType expected_variance = 600.25;
 	mitk::ImageStatisticsContainer::IndexType expected_minIndex;
 	expected_minIndex.set_size(3);
 	expected_minIndex[0] = 2;
 	expected_minIndex[1] = 2;
 	expected_minIndex[2] = 2;
 
 	mitk::ImageStatisticsContainer::IndexType expected_maxIndex;
 	expected_maxIndex.set_size(3);
 	expected_maxIndex[0] = 2;
 	expected_maxIndex[1] = 2;
 	expected_maxIndex[2] = 1;
 
 	mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New();
 	pfMaskGen->SetInputImage(m_US4DCroppedImage);
 	pfMaskGen->SetPlanarFigure(m_US4DCroppedPlanarFigure);
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer;
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_US4DCroppedImage, pfMaskGen.GetPointer()));
 	auto statisticsObjectTimestep1 = statisticsContainer->GetStatisticsForTimeStep(1);
 
 	VerifyStatistics(statisticsObjectTimestep1,
 		expected_N,
 		expected_mean,
 		expected_MPP,
 		expected_skewness,
 		expected_kurtosis,
 		expected_variance,
 		expected_standarddev,
 		expected_min,
 		expected_max,
 		expected_RMS,
 		expected_minIndex,
 		expected_maxIndex);
 }
 
 void mitkImageStatisticsCalculatorTestSuite::TestUS4DCroppedAllTimesteps()
 {
 	MITK_INFO << std::endl << "Test US4D cropped all timesteps:-----------------------------------------------------------------------------------";
 
 	std::string US4DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_cropped.nrrd");
 	m_US4DCroppedImage = mitk::IOUtil::Load<mitk::Image>(US4DCroppedFile);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading US4D_cropped", m_US4DCroppedImage.IsNotNull());
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer=mitk::ImageStatisticsContainer::New();
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_US4DCroppedImage));
 	for (int i = 0; i < 4; i++)
 	{
 		CPPUNIT_ASSERT_MESSAGE("Error computing statistics for multiple timestep", statisticsContainer->TimeStepExists(i));
 	}
 }
 
 void mitkImageStatisticsCalculatorTestSuite::TestUS4DCropped3DMask()
 {
 	MITK_INFO << std::endl << "Test US4D cropped with 3D binary Mask:-----------------------------------------------------------------------------------";
 
 	std::string US4DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_cropped.nrrd");
 	m_US4DCroppedImage = mitk::IOUtil::Load<mitk::Image>(US4DCroppedFile);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading US4D_cropped", m_US4DCroppedImage.IsNotNull());
 
 	std::string US4DCropped3DBinMaskFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_cropped3DBinMask.nrrd");
 	m_US4DCropped3DBinMask = mitk::IOUtil::Load<mitk::Image>(US4DCropped3DBinMaskFile);
 	CPPUNIT_ASSERT_MESSAGE("Failed loading Pic3D binary mask", m_US4DCropped3DBinMask.IsNotNull());
 	//calculated ground truth via script
 	mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1;
 	mitk::ImageStatisticsContainer::RealType expected_MPP = 198;
 	mitk::ImageStatisticsContainer::RealType expected_max = 199;
 	mitk::ImageStatisticsContainer::RealType expected_mean = 198;
 	mitk::ImageStatisticsContainer::RealType expected_min = 197;
 	mitk::ImageStatisticsContainer::VoxelCountType expected_N = 2;
 	mitk::ImageStatisticsContainer::RealType expected_RMS = 198.00252523642217;
 	mitk::ImageStatisticsContainer::RealType expected_skewness = 0;
 	mitk::ImageStatisticsContainer::RealType expected_standarddev = 1;
 	mitk::ImageStatisticsContainer::RealType expected_variance = 1;
 	mitk::ImageStatisticsContainer::IndexType expected_minIndex;
 	expected_minIndex.set_size(3);
 	expected_minIndex[0] = 1;
 	expected_minIndex[1] = 2;
 	expected_minIndex[2] = 1;
 
 	mitk::ImageStatisticsContainer::IndexType expected_maxIndex;
 	expected_maxIndex.set_size(3);
 	expected_maxIndex[0] = 1;
 	expected_maxIndex[1] = 1;
 	expected_maxIndex[2] = 1;
 
 	mitk::ImageMaskGenerator::Pointer imgMask1 = mitk::ImageMaskGenerator::New();
 	imgMask1->SetInputImage(m_US4DCroppedImage);
 	imgMask1->SetImageMask(m_US4DCropped3DBinMask);
 
 	mitk::ImageStatisticsContainer::Pointer statisticsContainer = mitk::ImageStatisticsContainer::New();
 	CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_US4DCroppedImage, imgMask1.GetPointer(), nullptr, 1));
 	auto statisticsObjectTimestep1 = statisticsContainer->GetStatisticsForTimeStep(1);
 
 	VerifyStatistics(statisticsObjectTimestep1,
 		expected_N,
 		expected_mean,
 		expected_MPP,
 		expected_skewness,
 		expected_kurtosis,
 		expected_variance,
 		expected_standarddev,
 		expected_min,
 		expected_max,
 		expected_RMS,
 		expected_minIndex,
 		expected_maxIndex);
 }
 
 mitk::PlanarPolygon::Pointer mitkImageStatisticsCalculatorTestSuite::GeneratePlanarPolygon(mitk::PlaneGeometry::Pointer geometry, std::vector <mitk::Point2D> points)
 {
 	mitk::PlanarPolygon::Pointer figure = mitk::PlanarPolygon::New();
 	figure->SetPlaneGeometry(geometry);
 	figure->PlaceFigure(points[0]);
         for (unsigned int i = 1; i < points.size(); i++)
 	{
 		figure->SetControlPoint(i, points[i], true);
 	}
 	return figure;
 }
 
 const mitk::ImageStatisticsContainer::Pointer
 mitkImageStatisticsCalculatorTestSuite::ComputeStatistics(mitk::Image::ConstPointer image,
 	mitk::MaskGenerator::Pointer maskGen,
 	mitk::MaskGenerator::Pointer secondardMaskGen,
 	unsigned short label)
 {
 	mitk::ImageStatisticsCalculator::Pointer imgStatCalc = mitk::ImageStatisticsCalculator::New();
 	imgStatCalc->SetInputImage(image);
 
 	if (maskGen.IsNotNull())
 	{
 		imgStatCalc->SetMask(maskGen.GetPointer());
 		if (secondardMaskGen.IsNotNull())
 		{
 			imgStatCalc->SetSecondaryMask(secondardMaskGen.GetPointer());
 		}
 	}
 
 	return imgStatCalc->GetStatistics(label);
 }
 
 void mitkImageStatisticsCalculatorTestSuite::VerifyStatistics(mitk::ImageStatisticsContainer::ImageStatisticsObject stats,
 	mitk::ImageStatisticsContainer::RealType testMean, mitk::ImageStatisticsContainer::RealType testSD, mitk::ImageStatisticsContainer::RealType testMedian)
 {
-	mitk::ImageStatisticsContainer::RealType meanObject;
-	mitk::ImageStatisticsContainer::RealType standardDeviationObject;
-	mitk::ImageStatisticsContainer::RealType medianObject;
+    mitk::ImageStatisticsContainer::RealType meanObject = 0;
+    mitk::ImageStatisticsContainer::RealType standardDeviationObject = 0;
+    mitk::ImageStatisticsContainer::RealType medianObject = 0;
 	CPPUNIT_ASSERT_NO_THROW(meanObject = stats.GetValueConverted<mitk::ImageStatisticsContainer::RealType>(mitk::ImageStatisticsConstants::MEAN()));
 	CPPUNIT_ASSERT_NO_THROW(standardDeviationObject = stats.GetValueConverted<mitk::ImageStatisticsContainer::RealType>(mitk::ImageStatisticsConstants::STANDARDDEVIATION()));
 	CPPUNIT_ASSERT_NO_THROW(medianObject = stats.GetValueConverted<mitk::ImageStatisticsContainer::RealType>(mitk::ImageStatisticsConstants::MEDIAN()));
 	CPPUNIT_ASSERT_MESSAGE("Calculated mean grayvalue is not equal to the desired value.", std::abs(meanObject - testMean) < mitk::eps);
 	CPPUNIT_ASSERT_MESSAGE("Calculated grayvalue sd is not equal to the desired value.", std::abs(standardDeviationObject - testSD) < mitk::eps);
 	CPPUNIT_ASSERT_MESSAGE("Calculated median grayvalue is not equal to the desired value.", std::abs(medianObject - testMedian) < mitk::eps);
 }
 
 // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy)
 void mitkImageStatisticsCalculatorTestSuite::VerifyStatistics(mitk::ImageStatisticsContainer::ImageStatisticsObject stats,
 	mitk::ImageStatisticsContainer::VoxelCountType N,
 	mitk::ImageStatisticsContainer::RealType mean,
 	mitk::ImageStatisticsContainer::RealType MPP,
 	mitk::ImageStatisticsContainer::RealType skewness,
 	mitk::ImageStatisticsContainer::RealType kurtosis,
 	mitk::ImageStatisticsContainer::RealType variance,
 	mitk::ImageStatisticsContainer::RealType stdev,
 	mitk::ImageStatisticsContainer::RealType min,
 	mitk::ImageStatisticsContainer::RealType max,
 	mitk::ImageStatisticsContainer::RealType RMS,
 	mitk::ImageStatisticsContainer::IndexType minIndex,
 	mitk::ImageStatisticsContainer::IndexType maxIndex)
 {
 	mitk::ImageStatisticsContainer::VoxelCountType numberOfVoxelsObject;
-	mitk::ImageStatisticsContainer::RealType meanObject;
+    mitk::ImageStatisticsContainer::RealType meanObject = 0;
 	mitk::ImageStatisticsContainer::RealType mppObject;
-	mitk::ImageStatisticsContainer::RealType skewnessObject;
-	mitk::ImageStatisticsContainer::RealType kurtosisObject;
-	mitk::ImageStatisticsContainer::RealType varianceObject;
-	mitk::ImageStatisticsContainer::RealType standardDeviationObject;
-	mitk::ImageStatisticsContainer::RealType minObject;
-	mitk::ImageStatisticsContainer::RealType maxObject;
-	mitk::ImageStatisticsContainer::RealType rmsObject;
+    mitk::ImageStatisticsContainer::RealType skewnessObject = 0;
+    mitk::ImageStatisticsContainer::RealType kurtosisObject = 0;
+    mitk::ImageStatisticsContainer::RealType varianceObject = 0;
+    mitk::ImageStatisticsContainer::RealType standardDeviationObject = 0;
+    mitk::ImageStatisticsContainer::RealType minObject = 0;
+    mitk::ImageStatisticsContainer::RealType maxObject = 0;
+    mitk::ImageStatisticsContainer::RealType rmsObject = 0;
 	mitk::ImageStatisticsContainer::IndexType minIndexObject;
 	mitk::ImageStatisticsContainer::IndexType maxIndexObject;
 	CPPUNIT_ASSERT_NO_THROW(numberOfVoxelsObject = stats.GetValueConverted<mitk::ImageStatisticsContainer::VoxelCountType>(mitk::ImageStatisticsConstants::NUMBEROFVOXELS()));
 	CPPUNIT_ASSERT_NO_THROW(meanObject = stats.GetValueConverted<mitk::ImageStatisticsContainer::RealType>(mitk::ImageStatisticsConstants::MEAN()));
 	CPPUNIT_ASSERT_NO_THROW(mppObject = stats.GetValueConverted<mitk::ImageStatisticsContainer::RealType>(mitk::ImageStatisticsConstants::MPP()));
 	CPPUNIT_ASSERT_NO_THROW(skewnessObject = stats.GetValueConverted<mitk::ImageStatisticsContainer::RealType>(mitk::ImageStatisticsConstants::SKEWNESS()));
 	CPPUNIT_ASSERT_NO_THROW(kurtosisObject = stats.GetValueConverted<mitk::ImageStatisticsContainer::RealType>(mitk::ImageStatisticsConstants::KURTOSIS()));
 	CPPUNIT_ASSERT_NO_THROW(varianceObject = stats.GetValueConverted<mitk::ImageStatisticsContainer::RealType>(mitk::ImageStatisticsConstants::VARIANCE()));
 	CPPUNIT_ASSERT_NO_THROW(standardDeviationObject = stats.GetValueConverted<mitk::ImageStatisticsContainer::RealType>(mitk::ImageStatisticsConstants::STANDARDDEVIATION()));
 	CPPUNIT_ASSERT_NO_THROW(minObject = stats.GetValueConverted<mitk::ImageStatisticsContainer::RealType>(mitk::ImageStatisticsConstants::MINIMUM()));
 	CPPUNIT_ASSERT_NO_THROW(maxObject = stats.GetValueConverted<mitk::ImageStatisticsContainer::RealType>(mitk::ImageStatisticsConstants::MAXIMUM()));
 	CPPUNIT_ASSERT_NO_THROW(rmsObject = stats.GetValueConverted<mitk::ImageStatisticsContainer::RealType>(mitk::ImageStatisticsConstants::RMS()));
 	CPPUNIT_ASSERT_NO_THROW(minIndexObject = stats.GetValueConverted<mitk::ImageStatisticsContainer::IndexType>(mitk::ImageStatisticsConstants::MINIMUMPOSITION()));
 	CPPUNIT_ASSERT_NO_THROW(maxIndexObject = stats.GetValueConverted<mitk::ImageStatisticsContainer::IndexType>(mitk::ImageStatisticsConstants::MAXIMUMPOSITION()));
 
 	CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", numberOfVoxelsObject - N == 0);
 	CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(meanObject - mean) < mitk::eps);
 	// in three test cases MPP is None because the roi has no positive pixels
 	if (!std::isnan(mppObject))
 	{
 		CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(mppObject - MPP) < mitk::eps);
 	}
 	CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(skewnessObject - skewness) < mitk::eps);
 	CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(kurtosisObject - kurtosis) < mitk::eps);
 	CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(varianceObject - variance) < mitk::eps);
 	CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(standardDeviationObject - stdev) < mitk::eps);
 	CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(minObject - min) < mitk::eps);
 	CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(maxObject - max) < mitk::eps);
 	CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(rmsObject - RMS) < mitk::eps);
 	for (unsigned int i = 0; i < minIndex.size(); ++i)
 	{
 		CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(minIndexObject[i] - minIndex[i]) < mitk::eps);
 	}
 	for (unsigned int i = 0; i < maxIndex.size(); ++i)
 	{
 		CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(maxIndexObject[i] - maxIndex[i]) < mitk::eps);
 	}
 }
 MITK_TEST_SUITE_REGISTRATION(mitkImageStatisticsCalculator)
diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp
index 879348f407..7171121930 100644
--- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp
+++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp
@@ -1,426 +1,461 @@
 /*===================================================================
 
 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.
 
 ===================================================================*/
 
 
 // Blueberry
 #include <berryISelectionService.h>
 #include <berryIWorkbenchWindow.h>
 
 // Qmitk
 #include "QmitkFiberQuantificationView.h"
 
 // Qt
 #include <QMessageBox>
 
 // MITK
 #include <mitkNodePredicateDataType.h>
 #include <mitkNodePredicateDimension.h>
 #include <mitkNodePredicateAnd.h>
 #include <mitkImageCast.h>
 #include <mitkPeakImage.h>
 #include <mitkLabelSetImage.h>
 #include <mitkDICOMSegmentationConstants.h>
 #include <mitkDICOMSegmentationPropertyHelper.cpp>
 #include <mitkDICOMQIPropertyHelper.h>
 
 // ITK
 #include <itkTractDensityImageFilter.h>
 #include <itkTractsToRgbaImageFilter.h>
 #include <itkTractsToVectorImageFilter.h>
 #include <itkTractsToFiberEndingsImageFilter.h>
 
 #include <mitkLexicalCast.h>
 
 const std::string QmitkFiberQuantificationView::VIEW_ID = "org.mitk.views.fiberquantification";
 using namespace mitk;
 
 QmitkFiberQuantificationView::QmitkFiberQuantificationView()
   : QmitkAbstractView()
   , m_Controls( 0 )
   , m_UpsamplingFactor(5)
   , m_Visible(false)
 {
 
 }
 
 // Destructor
 QmitkFiberQuantificationView::~QmitkFiberQuantificationView()
 {
 
 }
 
 void QmitkFiberQuantificationView::CreateQtPartControl( QWidget *parent )
 {
   // build up qt view, unless already done
   if ( !m_Controls )
   {
     // create GUI widgets from the Qt Designer's .ui file
     m_Controls = new Ui::QmitkFiberQuantificationViewControls;
     m_Controls->setupUi( parent );
 
     connect( m_Controls->m_ProcessFiberBundleButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedBundles()) );
     connect( m_Controls->m_ExtractFiberPeaks, SIGNAL(clicked()), this, SLOT(CalculateFiberDirections()) );
 
     m_Controls->m_TractBox->SetDataStorage(this->GetDataStorage());
     mitk::TNodePredicateDataType<mitk::FiberBundle>::Pointer isFib = mitk::TNodePredicateDataType<mitk::FiberBundle>::New();
     m_Controls->m_TractBox->SetPredicate( isFib );
 
     m_Controls->m_ImageBox->SetDataStorage(this->GetDataStorage());
     m_Controls->m_ImageBox->SetZeroEntryText("--");
     mitk::TNodePredicateDataType<mitk::Image>::Pointer isImagePredicate = mitk::TNodePredicateDataType<mitk::Image>::New();
     mitk::NodePredicateDimension::Pointer is3D = mitk::NodePredicateDimension::New(3);
     m_Controls->m_ImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, is3D) );
 
     connect( (QObject*)(m_Controls->m_TractBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()));
     connect( (QObject*)(m_Controls->m_ImageBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()));
   }
 }
 
 void QmitkFiberQuantificationView::Activated()
 {
 
 }
 
 void QmitkFiberQuantificationView::Deactivated()
 {
 
 }
 
 void QmitkFiberQuantificationView::Visible()
 {
   m_Visible = true;
 }
 
 void QmitkFiberQuantificationView::Hidden()
 {
   m_Visible = false;
 }
 
 void QmitkFiberQuantificationView::SetFocus()
 {
   m_Controls->m_ProcessFiberBundleButton->setFocus();
 }
 
 void QmitkFiberQuantificationView::CalculateFiberDirections()
 {
   typedef itk::Image<unsigned char, 3>                                            ItkUcharImgType;
 
   // load fiber bundle
   mitk::FiberBundle::Pointer inputTractogram = dynamic_cast<mitk::FiberBundle*>(m_SelectedFB.back()->GetData());
 
   itk::TractsToVectorImageFilter<float>::Pointer fOdfFilter = itk::TractsToVectorImageFilter<float>::New();
   if (m_SelectedImage.IsNotNull())
   {
     ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New();
     mitk::CastToItkImage<ItkUcharImgType>(m_SelectedImage, itkMaskImage);
     fOdfFilter->SetMaskImage(itkMaskImage);
   }
 
   // extract directions from fiber bundle
   fOdfFilter->SetFiberBundle(inputTractogram);
   fOdfFilter->SetAngularThreshold(cos(m_Controls->m_AngularThreshold->value()*itk::Math::pi/180));
   switch (m_Controls->m_FiberDirNormBox->currentIndex())
   {
   case 0:
     fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter<float>::NormalizationMethods::GLOBAL_MAX);
     break;
   case 1:
     fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter<float>::NormalizationMethods::SINGLE_VEC_NORM);
     break;
   case 2:
     fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter<float>::NormalizationMethods::MAX_VEC_NORM);
     break;
   }
   fOdfFilter->SetSizeThreshold(m_Controls->m_PeakThreshold->value());
   fOdfFilter->SetMaxNumDirections(m_Controls->m_MaxNumDirections->value());
   fOdfFilter->Update();
 
   QString name = m_SelectedFB.back()->GetName().c_str();
 
   if (m_Controls->m_NumDirectionsBox->isChecked())
   {
     mitk::Image::Pointer mitkImage = mitk::Image::New();
     mitkImage->InitializeByItk( fOdfFilter->GetNumDirectionsImage().GetPointer() );
     mitkImage->SetVolume( fOdfFilter->GetNumDirectionsImage()->GetBufferPointer() );
 
     mitk::DataNode::Pointer node = mitk::DataNode::New();
     node->SetData(mitkImage);
     node->SetName((name+"_NUM_DIRECTIONS").toStdString().c_str());
     GetDataStorage()->Add(node, m_SelectedFB.back());
   }
 
   Image::Pointer mitkImage = dynamic_cast<Image*>(PeakImage::New().GetPointer());
   mitk::CastToMitkImage(fOdfFilter->GetDirectionImage(), mitkImage);
   mitkImage->SetVolume(fOdfFilter->GetDirectionImage()->GetBufferPointer());
 
   mitk::DataNode::Pointer node = mitk::DataNode::New();
   node->SetData(mitkImage);
   node->SetName( (name+"_DIRECTIONS").toStdString().c_str());
   GetDataStorage()->Add(node, m_SelectedFB.back());
 }
 
 void QmitkFiberQuantificationView::UpdateGui()
 {
   m_SelectedFB.clear();
   if (m_Controls->m_TractBox->GetSelectedNode().IsNotNull())
     m_SelectedFB.push_back(m_Controls->m_TractBox->GetSelectedNode());
 
   m_SelectedImage = nullptr;
   if (m_Controls->m_ImageBox->GetSelectedNode().IsNotNull())
     m_SelectedImage = dynamic_cast<mitk::Image*>(m_Controls->m_ImageBox->GetSelectedNode()->GetData());
 
   m_Controls->m_ProcessFiberBundleButton->setEnabled(!m_SelectedFB.empty());
   m_Controls->m_ExtractFiberPeaks->setEnabled(!m_SelectedFB.empty());
 }
 
 void QmitkFiberQuantificationView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList<mitk::DataNode::Pointer>& )
 {
   UpdateGui();
 }
 
 void QmitkFiberQuantificationView::ProcessSelectedBundles()
 {
   if ( m_SelectedFB.empty() ){
     QMessageBox::information( nullptr, "Warning", "No fibe bundle selected!");
     MITK_WARN("QmitkFiberQuantificationView") << "no fibe bundle selected";
     return;
   }
 
   int generationMethod = m_Controls->m_GenerationBox->currentIndex();
 
   for( unsigned int i=0; i<m_SelectedFB.size(); i++ )
   {
     mitk::DataNode::Pointer node = m_SelectedFB[i];
     if (node.IsNotNull() && dynamic_cast<mitk::FiberBundle*>(node->GetData()))
     {
       mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
       QString name(node->GetName().c_str());
       DataNode::Pointer newNode = nullptr;
       switch(generationMethod){
       case 0:
-        newNode = GenerateTractDensityImage(fib, false, true);
+        newNode = GenerateTractDensityImage(fib, false, true, node->GetName());
         name += "_TDI";
         break;
       case 1:
-        newNode = GenerateTractDensityImage(fib, false, false);
+        newNode = GenerateTractDensityImage(fib, false, false, node->GetName());
         name += "_TDI";
         break;
       case 2:
-        newNode = GenerateTractDensityImage(fib, true, false);
+        newNode = GenerateTractDensityImage(fib, true, false, node->GetName());
         name += "_envelope";
         break;
       case 3:
         newNode = GenerateColorHeatmap(fib);
         break;
       case 4:
         newNode = GenerateFiberEndingsImage(fib);
         name += "_fiber_endings";
         break;
       case 5:
         newNode = GenerateFiberEndingsPointSet(fib);
         name += "_fiber_endings";
         break;
       }
       if (newNode.IsNotNull())
       {
         newNode->SetName(name.toStdString());
         GetDataStorage()->Add(newNode);
       }
     }
   }
 }
 
 // generate pointset displaying the fiber endings
 mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateFiberEndingsPointSet(mitk::FiberBundle::Pointer fib)
 {
   mitk::PointSet::Pointer pointSet = mitk::PointSet::New();
   vtkSmartPointer<vtkPolyData> fiberPolyData = fib->GetFiberPolyData();
 
   int count = 0;
   int numFibers = fib->GetNumFibers();
   for( int i=0; i<numFibers; i++ )
   {
     vtkCell* cell = fiberPolyData->GetCell(i);
     int numPoints = cell->GetNumberOfPoints();
     vtkPoints* points = cell->GetPoints();
 
     if (numPoints>0)
     {
       double* point = points->GetPoint(0);
       itk::Point<float,3> itkPoint = mitk::imv::GetItkPoint(point);
       pointSet->InsertPoint(count, itkPoint);
       count++;
     }
     if (numPoints>2)
     {
       double* point = points->GetPoint(numPoints-1);
       itk::Point<float,3> itkPoint = mitk::imv::GetItkPoint(point);
       pointSet->InsertPoint(count, itkPoint);
       count++;
     }
   }
 
   mitk::DataNode::Pointer node = mitk::DataNode::New();
   node->SetData( pointSet );
   return node;
 }
 
 // generate image displaying the fiber endings
 mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateFiberEndingsImage(mitk::FiberBundle::Pointer fib)
 {
   typedef unsigned int OutPixType;
   typedef itk::Image<OutPixType,3> OutImageType;
 
   typedef itk::TractsToFiberEndingsImageFilter< OutImageType > ImageGeneratorType;
   ImageGeneratorType::Pointer generator = ImageGeneratorType::New();
   generator->SetFiberBundle(fib);
   generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value());
   if (m_SelectedImage.IsNotNull())
   {
     OutImageType::Pointer itkImage = OutImageType::New();
     CastToItkImage(m_SelectedImage, itkImage);
     generator->SetInputImage(itkImage);
     generator->SetUseImageGeometry(true);
   }
   generator->Update();
 
   // get output image
   OutImageType::Pointer outImg = generator->GetOutput();
   mitk::Image::Pointer img = mitk::Image::New();
   img->InitializeByItk(outImg.GetPointer());
   img->SetVolume(outImg->GetBufferPointer());
 
   // init data node
   mitk::DataNode::Pointer node = mitk::DataNode::New();
   node->SetData(img);
   return node;
 }
 
 // generate rgba heatmap from fiber bundle
 mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateColorHeatmap(mitk::FiberBundle::Pointer fib)
 {
   typedef itk::RGBAPixel<unsigned char> OutPixType;
   typedef itk::Image<OutPixType, 3> OutImageType;
   typedef itk::TractsToRgbaImageFilter< OutImageType > ImageGeneratorType;
   ImageGeneratorType::Pointer generator = ImageGeneratorType::New();
   generator->SetFiberBundle(fib);
   generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value());
   if (m_SelectedImage.IsNotNull())
   {
     itk::Image<unsigned char, 3>::Pointer itkImage = itk::Image<unsigned char, 3>::New();
     CastToItkImage(m_SelectedImage, itkImage);
     generator->SetInputImage(itkImage);
     generator->SetUseImageGeometry(true);
   }
   generator->Update();
 
   // get output image
   typedef itk::Image<OutPixType,3> OutType;
   OutType::Pointer outImg = generator->GetOutput();
   mitk::Image::Pointer img = mitk::Image::New();
   img->InitializeByItk(outImg.GetPointer());
   img->SetVolume(outImg->GetBufferPointer());
 
   // init data node
   mitk::DataNode::Pointer node = mitk::DataNode::New();
   node->SetData(img);
   return node;
 }
 
 // generate tract density image from fiber bundle
-mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, bool binary, bool absolute)
+mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, bool binary, bool absolute, std::string name)
 {
   mitk::DataNode::Pointer node = mitk::DataNode::New();
   if (binary)
   {
     typedef unsigned char OutPixType;
     typedef itk::Image<OutPixType, 3> OutImageType;
 
     itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New();
     generator->SetFiberBundle(fib);
     generator->SetBinaryOutput(binary);
     generator->SetOutputAbsoluteValues(absolute);
     generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value());
     if (m_SelectedImage.IsNotNull())
     {
       OutImageType::Pointer itkImage = OutImageType::New();
       CastToItkImage(m_SelectedImage, itkImage);
       generator->SetInputImage(itkImage);
       generator->SetUseImageGeometry(true);
 
     }
     generator->Update();
 
     // get output image
     typedef itk::Image<OutPixType,3> OutType;
     OutType::Pointer outImg = generator->GetOutput();
     mitk::Image::Pointer img = mitk::Image::New();
     img->InitializeByItk(outImg.GetPointer());
     img->SetVolume(outImg->GetBufferPointer());
 
 
     if (m_SelectedImage.IsNotNull())
     {
       mitk::LabelSetImage::Pointer multilabelImage = mitk::LabelSetImage::New();
       multilabelImage->InitializeByLabeledImage(img);
       multilabelImage->GetActiveLabelSet()->SetActiveLabel(1);
       mitk::Label::Pointer label = multilabelImage->GetActiveLabel();
       label->SetName("Tractogram");
+
+      // Add Segmented Property Category Code Sequence tags (0062, 0003): Sequence defining the general category of this
+      // segment.
+      // (0008,0100) Code Value
+      label->SetProperty(
+            DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_CATEGORY_CODE_VALUE_PATH()).c_str(),
+            TemporoSpatialStringProperty::New("T-D000A"));
+
+      // (0008,0102) Coding Scheme Designator
+      label->SetProperty(
+            DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_CATEGORY_CODE_SCHEME_PATH()).c_str(),
+            TemporoSpatialStringProperty::New("SRT"));
+
+      // (0008,0104) Code Meaning
+      label->SetProperty(
+            DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_CATEGORY_CODE_MEANING_PATH()).c_str(),
+            TemporoSpatialStringProperty::New("Anatomical Structure"));
+      //------------------------------------------------------------
+      // Add Segmented Property Type Code Sequence (0062, 000F): Sequence defining the specific property type of this
+      // segment.
+      // (0008,0100) Code Value
+      label->SetProperty(
+            DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_TYPE_CODE_VALUE_PATH()).c_str(),
+            TemporoSpatialStringProperty::New("DUMMY"));
+
+      // (0008,0102) Coding Scheme Designator
+      label->SetProperty(
+            DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_TYPE_CODE_SCHEME_PATH()).c_str(),
+            TemporoSpatialStringProperty::New("SRT"));
+
+      // (0008,0104) Code Meaning
+      label->SetProperty(
+            DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_TYPE_CODE_MEANING_PATH()).c_str(),
+            TemporoSpatialStringProperty::New(name));
+
       mitk::DICOMQIPropertyHandler::DeriveDICOMSourceProperties(m_SelectedImage, multilabelImage);
 
       // init data node
       node->SetData(multilabelImage);
     }
     else
     {
       // init data node
       node->SetData(img);
     }
 
   }
   else
   {
     typedef float OutPixType;
     typedef itk::Image<OutPixType, 3> OutImageType;
 
     itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New();
     generator->SetFiberBundle(fib);
     generator->SetBinaryOutput(binary);
     generator->SetOutputAbsoluteValues(absolute);
     generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value());
     if (m_SelectedImage.IsNotNull())
     {
       OutImageType::Pointer itkImage = OutImageType::New();
       CastToItkImage(m_SelectedImage, itkImage);
       generator->SetInputImage(itkImage);
       generator->SetUseImageGeometry(true);
 
     }
     //generator->SetDoFiberResampling(false);
     generator->Update();
 
     // get output image
     typedef itk::Image<OutPixType,3> OutType;
     OutType::Pointer outImg = generator->GetOutput();
     mitk::Image::Pointer img = mitk::Image::New();
     img->InitializeByItk(outImg.GetPointer());
     img->SetVolume(outImg->GetBufferPointer());
 
     // init data node
     node->SetData(img);
   }
   return node;
 }
diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.h
index 87968f7086..cfdd9a446a 100644
--- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.h
+++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.h
@@ -1,86 +1,86 @@
 /*===================================================================
 
 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.
 
 ===================================================================*/
 
 #ifndef QmitkFiberQuantificationView_h
 #define QmitkFiberQuantificationView_h
 
 #include <QmitkAbstractView.h>
 #include "ui_QmitkFiberQuantificationViewControls.h"
 
 #include <mitkFiberBundle.h>
 #include <mitkPointSet.h>
 #include <itkCastImageFilter.h>
 #include <mitkILifecycleAwarePart.h>
 
 /*!
 \brief Generation of images from fiber bundles (TDI, envelopes, endpoint distribution) and extraction of principal fiber directions from tractograms.
 
 */
 class QmitkFiberQuantificationView : public QmitkAbstractView, public mitk::ILifecycleAwarePart
 {
   // this is needed for all Qt objects that should have a Qt meta-object
   // (everything that derives from QObject and wants to have signal/slots)
   Q_OBJECT
 
 public:
 
   typedef itk::Image< unsigned char, 3 >    itkUCharImageType;
 
   static const std::string VIEW_ID;
 
   QmitkFiberQuantificationView();
   virtual ~QmitkFiberQuantificationView();
 
   virtual void CreateQtPartControl(QWidget *parent) override;
 
   ///
   /// Sets the focus to an internal widget.
   ///
   virtual void SetFocus() override;
 
   virtual void Activated() override;
   virtual void Deactivated() override;
   virtual void Visible() override;
   virtual void Hidden() override;
 
 protected slots:
 
   void ProcessSelectedBundles();    ///< start selected operation on fiber bundle (e.g. tract density image generation)
   void CalculateFiberDirections();  ///< Calculate main fiber directions from tractogram
   void UpdateGui();     ///< update button activity etc. dpending on current datamanager selection
 
 protected:
 
   /// \brief called by QmitkAbstractView when DataManager's selection has changed
   virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList<mitk::DataNode::Pointer>& nodes) override;
 
   Ui::QmitkFiberQuantificationViewControls* m_Controls;
 
   std::vector<mitk::DataNode::Pointer>  m_SelectedFB;       ///< selected fiber bundle nodes
   mitk::Image::Pointer                  m_SelectedImage;
   float                                 m_UpsamplingFactor; ///< upsampling factor for all image generations
 
-  mitk::DataNode::Pointer GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, bool binary, bool absolute);
+  mitk::DataNode::Pointer GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, bool binary, bool absolute, std::string name);
   mitk::DataNode::Pointer GenerateColorHeatmap(mitk::FiberBundle::Pointer fib);
   mitk::DataNode::Pointer GenerateFiberEndingsImage(mitk::FiberBundle::Pointer fib);
   mitk::DataNode::Pointer GenerateFiberEndingsPointSet(mitk::FiberBundle::Pointer fib);
   bool  m_Visible;
 };
 
 
 
 #endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED
 
diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.python/resources/tractseg.py b/Plugins/org.mitk.gui.qt.diffusionimaging.python/resources/tractseg.py
index e598a1cc1c..8acde177b5 100644
--- a/Plugins/org.mitk.gui.qt.diffusionimaging.python/resources/tractseg.py
+++ b/Plugins/org.mitk.gui.qt.diffusionimaging.python/resources/tractseg.py
@@ -1,112 +1,115 @@
 error_string = None
 del error_string
 try:
     import nibabel as nib
     import numpy as np
     from tractseg.python_api import run_tractseg
     from tractseg.python_api import img_utils
 
     data = sitk.GetArrayFromImage(in_image)
     data = np.nan_to_num(data)
     swapaxes = False
     if data.shape != (sx, sy, sz):
         print("Swapping axes!")
         data = np.swapaxes(data, 0, 2)
         swapaxes = True
 
     affine = np.zeros((4, 4))
     affine[0, 0] = in_image.GetDirection()[0] * in_image.GetSpacing()[0]
     affine[1, 0] = in_image.GetDirection()[1] * in_image.GetSpacing()[0]
     affine[2, 0] = in_image.GetDirection()[2] * in_image.GetSpacing()[0]
     affine[0, 1] = in_image.GetDirection()[3] * in_image.GetSpacing()[1]
     affine[1, 1] = in_image.GetDirection()[4] * in_image.GetSpacing()[1]
     affine[2, 1] = in_image.GetDirection()[5] * in_image.GetSpacing()[1]
     affine[0, 2] = in_image.GetDirection()[6] * in_image.GetSpacing()[2]
     affine[1, 2] = in_image.GetDirection()[7] * in_image.GetSpacing()[2]
     affine[2, 2] = in_image.GetDirection()[8] * in_image.GetSpacing()[2]
     affine[0, 3] = in_image.GetOrigin()[0]
     affine[1, 3] = in_image.GetOrigin()[1]
     affine[2, 3] = in_image.GetOrigin()[2]
     affine[3, 3] = 1
 
-    data, flip_axis = img_utils.flip_peaks_to_correct_orientation_if_needed(nib.Nifti1Image(data, affine=affine), do_flip=True)
+    data, flip_axis = img_utils.flip_peaks_to_correct_orientation_if_needed(nib.Nifti1Image(data, affine=affine),
+                                                                            do_flip=True)
     print('flip_axis', flip_axis)
 
     print('output_type', output_type)
     print('get_probs', get_probs)
     print('dropout_sampling', dropout_sampling)
     print('threshold', threshold)
-    seg = run_tractseg(data=data, output_type=output_type, input_type="peaks", verbose=verbose, get_probs=get_probs, dropout_sampling=dropout_sampling, threshold=threshold, postprocess=True)
 
-    bla = nib.Nifti1Image(seg, affine)
-    nib.save(bla, '/home/neher/test.nii.gz')
+    seg = run_tractseg(data=data, output_type=output_type, input_type="peaks", verbose=verbose, get_probs=get_probs,
+                       dropout_sampling=dropout_sampling, threshold=threshold, postprocess=False)
+
+    # bla = nib.Nifti1Image(seg, affine)
+    # nib.save(bla, '/home/neher/test.nii.gz')
     if swapaxes:
+        print("Swapping axes back!")
         seg = np.swapaxes(seg, 0, 2)
 
+    print('Output shape: ' + str(seg.shape))
     if output_type == "tract_segmentation":
 
-        print("tract_segmentation", get_probs, dropout_sampling)
         if not get_probs and not dropout_sampling:
             if collapse:
                 temp = np.zeros((seg.shape[0], seg.shape[1], seg.shape[2]))
                 for l in range(seg.shape[3]):
                     temp[np.where(seg[:, :, :, l] > 0)] = l + 1
                 seg = temp
             segmentation = sitk.GetImageFromArray(seg.astype(np.uint8))
         else:
             segmentation = sitk.GetImageFromArray(seg)
 
         segmentation.SetOrigin(in_image.GetOrigin())
         segmentation.SetSpacing(in_image.GetSpacing())
         segmentation.SetDirection(in_image.GetDirection())
 
     elif output_type == "endings_segmentation":
 
         print("endings_segmentation", get_probs, dropout_sampling)
         if not get_probs and not dropout_sampling:
 
             # merge start and end into labelmap
             temp = np.zeros((seg.shape[0], seg.shape[1], seg.shape[2], seg.shape[3] // 2))
             for l in range(temp.shape[3]):
                 temp[:, :, :, l][np.where(seg[:, :, :, l * 2] > 0)] = 1
                 temp[:, :, :, l][np.where(seg[:, :, :, l * 2 + 1] > 0)] = 2
             seg = temp
 
             if collapse:
                 temp = np.zeros((seg.shape[0], seg.shape[1], seg.shape[2]))
                 for l in range(seg.shape[3]):
                     temp[np.where(seg[:, :, :, l] == 1)] = 2*l + 1
                     temp[np.where(seg[:, :, :, l] == 2)] = 2*l + 2
                 seg = temp
             segmentation = sitk.GetImageFromArray(seg.astype(np.uint8))
         else:
 
             temp = np.zeros((seg.shape[0], seg.shape[1], seg.shape[2], seg.shape[3] // 2))
             for l in range(temp.shape[3]):
                 temp[:, :, :, l] += seg[:, :, :, l * 2]
                 temp[:, :, :, l] += seg[:, :, :, l * 2 + 1]
             seg = temp
             segmentation = sitk.GetImageFromArray(seg)
 
         segmentation.SetOrigin(in_image.GetOrigin())
         segmentation.SetSpacing(in_image.GetSpacing())
         segmentation.SetDirection(in_image.GetDirection())
 
     elif output_type == "TOM":
-
         if not collapse:
             for x in range(0, 72):
                 tom = sitk.GetImageFromArray(seg[:, :, :, x * 3:x * 3 + 3])
                 tom.SetOrigin(in_image.GetOrigin())
                 tom.SetSpacing(in_image.GetSpacing())
                 tom.SetDirection(in_image.GetDirection())
                 globals()['tom%s' % x] = tom
         else:
             tom0 = sitk.GetImageFromArray(seg)
             tom0.SetOrigin(in_image.GetOrigin())
             tom0.SetSpacing(in_image.GetSpacing())
             tom0.SetDirection(in_image.GetDirection())
 
 except Exception as e:
     error_string = str(e)
     print(error_string)