diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.cxx b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.cxx
index 556ca7b45f..6cbae05b64 100644
--- a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.cxx
+++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.cxx
@@ -1,454 +1,457 @@
 /*===================================================================
 
 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 DIFFUSIONKURTOSISRECONSTRUCTIONIMAGEFILTER_CXX
 #define DIFFUSIONKURTOSISRECONSTRUCTIONIMAGEFILTER_CXX
 
 #include "itkDiffusionKurtosisReconstructionImageFilter.h"
 
 #include <itkImageRegionConstIterator.h>
 #include <itkImageRegionConstIteratorWithIndex.h>
 #include <itkImageRegionIteratorWithIndex.h>
 #include <itkImageRegionIterator.h>
 
 #include <itkVectorIndexSelectionCastImageFilter.h>
 #include <itkComposeImageFilter.h>
 #include <itkDiscreteGaussianImageFilter.h>
 
 template< class TInputPixelType>
 static void FitSingleVoxel( const itk::VariableLengthVector< TInputPixelType > &input, const vnl_vector<double>& bvalues, vnl_vector<double>& result, bool omit_bzero = true)
 {
   // check for length
   assert( input.Size() == bvalues.size() );
 
   // assembly data vectors for fitting
   auto bvalueIter = bvalues.begin();
   unsigned int unused_values = 0;
   while( bvalueIter != bvalues.end() )
   {
     if( *bvalueIter < vnl_math::eps && omit_bzero )
     {
       ++unused_values;
     }
     ++bvalueIter;
   }
 
   // initialize data vectors with the estimated size (after filtering)
   vnl_vector<double> fit_measurements( input.Size() - unused_values, 0 );
   vnl_vector<double> fit_bvalues( input.Size() - unused_values, 0 );
 
   bvalueIter = bvalues.begin();
   unsigned int running_index = 0;
   unsigned int skip_count = 0;
   while( bvalueIter != bvalues.end() )
   {
     if( *bvalueIter < vnl_math::eps && omit_bzero )
     {
       ++skip_count;
     }
     else
     {
       fit_measurements[ running_index - skip_count ]  = input.GetElement(running_index);
       fit_bvalues[ running_index - skip_count] = *bvalueIter;
     }
 
     ++running_index;
     ++bvalueIter;
   }
 
 
   MITK_DEBUG("KurtosisFilter.FitSingleVoxel.Meas") << fit_measurements;
   MITK_DEBUG("KurtosisFilter.FitSingleVoxel.Bval") << fit_bvalues;
 
   // perform fit on data vectors
   if( omit_bzero )
   {
    itk::kurtosis_fit_omit_unweighted kurtosis_cost_fn( fit_measurements.size() );
    kurtosis_cost_fn.initialize( fit_measurements, fit_bvalues );
 
    vnl_levenberg_marquardt nonlinear_fit( kurtosis_cost_fn );
    nonlinear_fit.minimize( result );
   }
   else
   {
    itk::kurtosis_fit_lsq_function kurtosis_cost_fn( fit_measurements.size() );
    kurtosis_cost_fn.initialize( fit_measurements, fit_bvalues );
    kurtosis_cost_fn.use_bounds();
 
    vnl_levenberg_marquardt nonlinear_fit( kurtosis_cost_fn );
    nonlinear_fit.minimize(result);
   }
 
   MITK_DEBUG("KurtosisFilter.FitSingleVoxel.Rslt") << result;
 }
 
 
 
 template< class TInputPixelType, class TOutputPixelType>
 itk::DiffusionKurtosisReconstructionImageFilter<TInputPixelType, TOutputPixelType>
 ::DiffusionKurtosisReconstructionImageFilter()
   : m_ReferenceBValue(-1),
     m_OmitBZero(false),
     m_MaskImage(nullptr),
     m_ApplyPriorSmoothing(false),
-    m_SmoothingSigma(1.5)
+    m_SmoothingSigma(1.5),
+    m_MaxFitBValue( 3000 )
 {
   this->m_InitialPosition = vnl_vector<double>(3, 0);
   this->m_InitialPosition[2] = 1000.0; // S_0
   this->m_InitialPosition[0] = 0.001; // D
   this->m_InitialPosition[1] = 1; // K
 
+  this->m_KurtosisBounds.fill(0);
+
   // single input image
   this->SetNumberOfRequiredInputs(1);
 
   // two output images (D, K)
   this->SetNumberOfRequiredOutputs(2);
   typename OutputImageType::Pointer outputPtr1 = OutputImageType::New();
   typename OutputImageType::Pointer outputPtr2 = OutputImageType::New();
 
   this->SetNthOutput(0, outputPtr1.GetPointer() );
   this->SetNthOutput(1, outputPtr2.GetPointer() );
 
 }
 
 template< class TInputPixelType, class TOutputPixelType>
 void itk::DiffusionKurtosisReconstructionImageFilter<TInputPixelType, TOutputPixelType>
 ::GenerateOutputInformation()
 {
   // first call superclass
   Superclass::GenerateOutputInformation();
 }
 
 template< class TInputPixelType, class TOutputPixelType>
 void itk::DiffusionKurtosisReconstructionImageFilter<TInputPixelType, TOutputPixelType>
 ::SetImageMask(MaskImageType::Pointer mask)
 {
   this->m_MaskImage = mask;
 }
 
 template< class TInputPixelType, class TOutputPixelType>
 void itk::DiffusionKurtosisReconstructionImageFilter<TInputPixelType, TOutputPixelType>
 ::BeforeThreadedGenerateData()
 {
   // if we have a region set, convert it to a mask image, which is to be used as default for telling
   // we need to set the image anyway, so by default the mask is overall 1
   if( m_MaskImage.IsNull() )
   {
     m_MaskImage = MaskImageType::New();
     m_MaskImage->SetRegions( this->GetInput()->GetLargestPossibleRegion() );
     m_MaskImage->CopyInformation( this->GetInput() );
     m_MaskImage->Allocate();
 
     if( this->m_MapOutputRegion.GetNumberOfPixels() > 0 )
     {
 
       m_MaskImage->FillBuffer(0);
 
       typedef itk::ImageRegionIteratorWithIndex< MaskImageType > MaskIteratorType;
       MaskIteratorType maskIter( this->m_MaskImage, this->m_MapOutputRegion );
       maskIter.GoToBegin();
 
       while( !maskIter.IsAtEnd() )
       {
         maskIter.Set( 1 );
         ++maskIter;
       }
 
     }
     else
     {
       m_MaskImage->FillBuffer(1);
     }
 
 
   }
 
 
   // apply smoothing to the input image
   if( this->m_ApplyPriorSmoothing )
   {
     // filter typedefs
     typedef itk::DiscreteGaussianImageFilter< itk::Image<TInputPixelType, 3 >, itk::Image<TInputPixelType, 3 > > GaussianFilterType;
     typedef itk::VectorIndexSelectionCastImageFilter< InputImageType, itk::Image<TInputPixelType, 3 > > IndexSelectionType;
     typedef itk::ComposeImageFilter< itk::Image<TInputPixelType, 3>, InputImageType > ComposeFilterType;
 
 
     auto vectorImage = this->GetInput();
     typename IndexSelectionType::Pointer indexSelectionFilter = IndexSelectionType::New();
     indexSelectionFilter->SetInput( vectorImage );
 
     typename ComposeFilterType::Pointer vec_composer = ComposeFilterType::New();
 
     for( unsigned int i=0; i<vectorImage->GetVectorLength(); ++i)
     {
       typename GaussianFilterType::Pointer gaussian_filter = GaussianFilterType::New();
 
       indexSelectionFilter->SetIndex( i );
 
       gaussian_filter->SetInput( indexSelectionFilter->GetOutput() );
       gaussian_filter->SetVariance( m_SmoothingSigma );
 
       vec_composer->SetInput(i, gaussian_filter->GetOutput() );
 
       gaussian_filter->Update();
     }
 
     try
     {
       vec_composer->Update();
     }
     catch(const itk::ExceptionObject &e)
     {
       mitkThrow() << "[VectorImage.GaussianSmoothing] !! Failed with ITK Exception: " << e.what();
     }
 
     this->m_ProcessedInputImage = vec_composer->GetOutput();
   }
   else
   {
     this->m_ProcessedInputImage = const_cast<InputImageType*>( this->GetInput() );
   }
 }
 
 template< class TInputPixelType, class TOutputPixelType>
 void itk::DiffusionKurtosisReconstructionImageFilter<TInputPixelType, TOutputPixelType>
 ::AfterThreadedGenerateData()
 {
  /* // initialize buffer to zero overall, but don't forget the requested region pointer
   for( unsigned int i=0; i<this->GetNumberOfOutputs(); ++i)
   {
     typename OutputImageType::Pointer output = this->GetOutput(i);
 
     // after generating, set the buffered region to max, we have taken care about filling it with valid data in
     // UpdateOutputInformation, if not set, a writer (or any filter using the LargestPossibleRegion() during update may
     // complain about not getting the requested region
     output->SetBufferedRegion( output->GetLargestPossibleRegion() );
 
     std::cout << "[DiffusionKurtosisReconstructionImageFilter.After]" << output->GetLargestPossibleRegion() << std::endl;
   }*/
 
 }
 
 template< class TInputPixelType, class TOutputPixelType>
 typename itk::DiffusionKurtosisReconstructionImageFilter<TInputPixelType, TOutputPixelType>::KurtosisSnapshot
 itk::DiffusionKurtosisReconstructionImageFilter<TInputPixelType, TOutputPixelType>
 ::GetSnapshot(const itk::VariableLengthVector<TInputPixelType> &input, GradientDirectionContainerType::Pointer gradients, float bvalue, bool omit_bzero)
 {
   // initialize bvalues from reference value and the gradients provided on input
   this->SetReferenceBValue(bvalue);
   this->SetGradientDirections( gradients );
 
   // call the other method
   return this->GetSnapshot( input, this->m_BValues, omit_bzero );
 }
 
 
 template< class TInputPixelType, class TOutputPixelType>
 typename itk::DiffusionKurtosisReconstructionImageFilter<TInputPixelType, TOutputPixelType>::KurtosisSnapshot
 itk::DiffusionKurtosisReconstructionImageFilter<TInputPixelType, TOutputPixelType>
 ::GetSnapshot(const itk::VariableLengthVector<TInputPixelType> &input, vnl_vector<double> bvalues, bool omit_bzero)
 {
   // initialize
   vnl_vector<double> initial_position;
   if( !omit_bzero )
   {
     initial_position.set_size(2);
     initial_position[0] = this->m_InitialPosition[0];
     initial_position[1] = this->m_InitialPosition[1];
   }
   else
   {
     initial_position.set_size(3);
     initial_position = this->m_InitialPosition;
   }
 
   // fit
   FitSingleVoxel( input, bvalues, initial_position, omit_bzero );
 
   // assembly result snapshot
   KurtosisSnapshot result;
   result.m_D = initial_position[0];
   result.m_K = initial_position[1];
 
   if( omit_bzero )
   {
     result.m_f = 1 - initial_position[2] / std::fmax(0.01, input.GetElement(0));
     result.m_BzeroFit = initial_position[2];
   }
   else
     result.m_f = 1;
 
   // assembly data vectors for fitting
   auto bvalueIter = bvalues.begin();
   unsigned int unused_values = 0;
   while( bvalueIter != bvalues.end() )
   {
     if( *bvalueIter < vnl_math::eps && omit_bzero )
     {
       ++unused_values;
     }
     ++bvalueIter;
   }
 
   // initialize data vectors with the estimated size (after filtering)
   vnl_vector<double> fit_measurements( input.Size() - unused_values, 0 );
   vnl_vector<double> fit_bvalues( input.Size() - unused_values, 0 );
 
   // original data
   vnl_vector<double> orig_measurements( input.Size(), 0 );
 
   bvalueIter = bvalues.begin();
   unsigned int running_index = 0;
   unsigned int skip_count = 0;
   while( bvalueIter != bvalues.end() )
   {
     if( *bvalueIter < vnl_math::eps && omit_bzero )
     {
       ++skip_count;
     }
     else
     {
       fit_measurements[ running_index - skip_count ]  = input.GetElement(running_index);
       fit_bvalues[ running_index - skip_count ] = *bvalueIter;
     }
 
     orig_measurements[ running_index ] = input.GetElement(running_index);
 
     ++running_index;
     ++bvalueIter;
   }
 
   result.fit_bvalues = fit_bvalues;
   result.fit_measurements = fit_measurements;
 
   result.bvalues = bvalues;
   result.measurements = orig_measurements;
 
   result.m_fittedBZero = omit_bzero;
 
   return result;
 }
 
 template< class TInputPixelType, class TOutputPixelType>
 void itk::DiffusionKurtosisReconstructionImageFilter<TInputPixelType, TOutputPixelType>
 ::SetGradientDirections(GradientDirectionContainerType::Pointer gradients)
 {
   if( this->m_ReferenceBValue < 0)
   {
     itkExceptionMacro( << "Specify reference b-value prior to setting the gradient directions." );
   }
 
   if( gradients->Size() == 0 )
   {
     itkExceptionMacro( << "Empty gradient directions container retrieved" );
   }
 
   vnl_vector<double> vnl_bvalues( gradients->Size(), 0 );
 
   // initialize the b-values
   auto grIter = gradients->Begin();
   unsigned int index = 0;
   while( grIter != gradients->End() )
   {
     const double twonorm = grIter.Value().two_norm();
     vnl_bvalues( index++ ) = this->m_ReferenceBValue * twonorm * twonorm;
 
     ++grIter;
   }
 
   this->m_BValues = vnl_bvalues;
 
 
 }
 
 template< class TInputPixelType, class TOutputPixelType>
 void itk::DiffusionKurtosisReconstructionImageFilter<TInputPixelType, TOutputPixelType>
 ::SetInitialSolution(const vnl_vector<double>& x0 )
 {
   unsigned int param_size = 2 + static_cast<int>( this->m_OmitBZero );
   assert( x0.size() == param_size );
 
   this->m_InitialPosition = x0;
 }
 
 template< class TInputPixelType, class TOutputPixelType>
 void itk::DiffusionKurtosisReconstructionImageFilter<TInputPixelType, TOutputPixelType>
 ::ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, ThreadIdType /*threadId*/)
 {
   std::cout << "[ThreadedGenerateData]" << outputRegionForThread;
 
   typename OutputImageType::Pointer dImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0));
   itk::ImageRegionIteratorWithIndex< OutputImageType > dImageIt(dImage, outputRegionForThread);
   dImageIt.GoToBegin();
 
   typename OutputImageType::Pointer kImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(1));
   itk::ImageRegionIteratorWithIndex< OutputImageType > kImageIt(kImage, outputRegionForThread);
   kImageIt.GoToBegin();
 
   typedef itk::ImageRegionConstIteratorWithIndex< InputImageType > InputIteratorType;
   InputIteratorType inputIter( m_ProcessedInputImage, outputRegionForThread );
   inputIter.GoToBegin();
 
   typedef itk::ImageRegionConstIteratorWithIndex< MaskImageType > MaskIteratorType;
   MaskIteratorType maskIter( this->m_MaskImage, outputRegionForThread );
   maskIter.GoToBegin();
 
   vnl_vector<double> initial_position;
   if( !this->m_OmitBZero )
   {
     initial_position.set_size(2);
     initial_position[0] = this->m_InitialPosition[0];
     initial_position[1] = this->m_InitialPosition[1];
 
   }
   else
   {
     initial_position.set_size(3);
     initial_position = this->m_InitialPosition;
   }
 
   while( !inputIter.IsAtEnd() )
   {
     // set (reset) each iteration
     vnl_vector<double> result = initial_position;
 
     // fit single voxel (if inside mask )
     if( maskIter.Get() > 0 )
     {
       FitSingleVoxel( inputIter.Get(), this->m_BValues, result, this->m_OmitBZero );
     }
     else
     {
       result.fill(0);
     }
 
     // regardless the fit type, the parameters are always in the first two position
     // of the results vector
     dImageIt.Set( result[0] );
     kImageIt.Set( result[1] );
 
     //std::cout << "[Kurtosis.Fit]" << inputIter.GetIndex() << " --> " << dImageIt.GetIndex() << " result: " << result << std::endl;
 
     ++maskIter;
     ++inputIter;
     ++dImageIt;
     ++kImageIt;
   }
 
 }
 
 #endif // guards
diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.h
index ada8a7de82..07f9c83339 100644
--- a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.h
+++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.h
@@ -1,331 +1,364 @@
 /*===================================================================
 
 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 DIFFUSIONKURTOSISRECONSTRUCTIONIMAGEFILTER_H
 #define DIFFUSIONKURTOSISRECONSTRUCTIONIMAGEFILTER_H
 
 #include "itkImageToImageFilter.h"
 #include "itkVectorImage.h"
 
 #include "mitkDiffusionPropertyHelper.h"
 
 // vnl includes
 #include <vnl/algo/vnl_levenberg_marquardt.h>
 #include <vnl/vnl_least_squares_function.h>
 
 namespace itk
 {
 
   // Fitting routines
 
   /** @struct Kurtosis_fit_lsq_function
       @brief A base least-squares function for the diffusion kurtosis fit (non-IVIM)
 
       The basic function fits the signal to S_0 = S * exp [ -b * D + -b^2 * D^2 * K^2 ]
       */
   struct kurtosis_fit_lsq_function :
       public vnl_least_squares_function
   {
   public:
     /** full lsq_function constructor */
     kurtosis_fit_lsq_function( unsigned int num_params, unsigned int num_measurements, UseGradient g=no_gradient)
       : vnl_least_squares_function( num_params, num_measurements, g),
         m_use_bounds(false)
     {}
 
     /** simplified constructor for the 2-parameters fit */
     kurtosis_fit_lsq_function( unsigned int number_measurements)
       : kurtosis_fit_lsq_function( 2, number_measurements, no_gradient )
     {}
 
 
     /** Initialize the function by setting measurements and the corresponding b-values */
     void initialize( vnl_vector< double > const& _meas, vnl_vector< double> const& _bvals )
     {
       meas = _meas;
       bvalues = _bvals;
     }
 
 
     /** use penalty terms on fitting to force the parameters stay within the default bounds */
     void use_bounds()
     {
       m_use_bounds = true;
 
       // initialize bounds
       double upper_bounds[2] = {4e-3, 4 };
       kurtosis_upper_bounds = vnl_vector<double>(2, 2, upper_bounds);
       kurtosis_lower_bounds = vnl_vector<double>(2, 0);
     }
 
     virtual void f(const vnl_vector<double> &x, vnl_vector<double> &fx)
     {
       for ( unsigned int s=0; s < fx.size(); s++ )
       {
         const double factor = ( meas[s] - M(x, s) );
         fx[s] = factor * factor + penalty_term(x);
       }
     }
 
   protected:
 
     /** Formula for diffusion term, use for internal computations */
     double Diff( double x1, double x2, double b)
     {
       const double quotient = -1. * b * x1 + b*b * x1 * x1 * x2 / 6;
       return exp(quotient);
     }
 
     /** The fitting measurement function, has to be reimplemented in the classes */
     virtual double M( vnl_vector<double> const& x, unsigned int idx )
     {
       const double bvalue = bvalues[idx];
       return meas[0] *  Diff( x[0], x[1], bvalue);
     }
 
     /** Penalty term on D and K during fitting, make sure the vector that is passed in contains (D, K) in this ordering */
     double penalty_term( vnl_vector<double> const& x)
     {
       double penalty = 0;
 
       // skip when turned off
       if( !m_use_bounds )
         return penalty;
 
       for( unsigned int i=0; i< x.size(); i++)
       {
         // 5% penalty boundary
         // use exponential function to scale the penalty (max when x[i] == bounds )
         double penalty_boundary = 0.05 * (kurtosis_upper_bounds[i] - kurtosis_lower_bounds[i]);
 
         if( x[i] < kurtosis_lower_bounds[i] + penalty_boundary )
         {
           penalty += 1e6 * exp( -1 * ( x[i] - kurtosis_lower_bounds[i]) / penalty_boundary );
         }
         else if ( x[i] > kurtosis_upper_bounds[i] - penalty_boundary )
         {
           penalty += 1e6 * exp( -1 * ( kurtosis_upper_bounds[i] - x[i]) / penalty_boundary );
         }
       }
 
       return penalty;
     }
 
     bool m_use_bounds;
 
     vnl_vector<double> kurtosis_upper_bounds;
     vnl_vector<double> kurtosis_lower_bounds;
 
     vnl_vector<double> meas;
     vnl_vector<double> bvalues;
 
   };
 
   /** @struct kurtosis_fit_omit_unweighted
       @brief A fitting function handling the unweighted signal b_0 as a fitted parameter */
   struct kurtosis_fit_omit_unweighted
       : public kurtosis_fit_lsq_function
   {
   public:
     /** simplified constructor for the 3-parameters fit */
     kurtosis_fit_omit_unweighted( unsigned int number_measurements)
       : kurtosis_fit_lsq_function( 3, number_measurements, no_gradient )
     {}
 
   protected:
     virtual double M(const vnl_vector<double> &x, unsigned int idx) override
     {
       const double bvalue = bvalues[idx];
       return x[2] *  Diff( x[0], x[1], bvalue);
     }
 
     virtual double penalty_term(const vnl_vector<double> &x)
     {
       // grab last two elements (D, K) from the 3-param fit and get penalty from superclass
       vnl_vector<double> xx(x.data_block()[0], 2) ;
       return kurtosis_fit_lsq_function::penalty_term(xx);
     }
   };
 
 
 /**
   @class DiffusionKurtosisReconstructionImageFilter
   @brief This filter provides the fit of the kurtosis (non-IVIM) signal to the data
 
   It has two main modes of operation, either as an image filter to compute the D and K maps, i.e. fitting the values to each voxel or a computation on a single voxel
   or a voxel group (with mask) can be triggered by @sa GetSnapshot, GetCurrentSnapshot methods.
   */
 template< class TInputPixelType, class TOutputPixelType >
 class DiffusionKurtosisReconstructionImageFilter :
     public ImageToImageFilter< VectorImage< TInputPixelType, 3>, Image<TOutputPixelType, 3> >
 {
 public:
 
   /**
     @struct KurtosisSnapsho
 
     @brief Struct describing a result (and the data) of a Kurtosis model fit
     */
   struct KurtosisSnapshot
   {
     // input data structures
     //vnl_vector<double> filtered_measurements;
     vnl_vector<double> bvalues;
     vnl_vector<double> measurements;
 
     vnl_vector<double> fit_bvalues;
     vnl_vector<double> fit_measurements;
     vnl_vector<unsigned int> weighted_image_indices;
 
     bool m_fittedBZero;
 
     // variables holding the fitted values
     double m_f;
     double m_BzeroFit;
     double m_D;
     double m_K;
   };
 
+  enum FitScale
+  {
+    STRAIGHT = 0,
+    LOGARITHMIC
+  };
+
   //-- class typedefs
   typedef DiffusionKurtosisReconstructionImageFilter Self;
   typedef SmartPointer<Self>                      Pointer;
   typedef SmartPointer<const Self>                ConstPointer;
   typedef ImageToImageFilter< VectorImage< TInputPixelType, 3>,
                               Image< TOutputPixelType,3 > >   Superclass;
 
   /** Method for creation through the object factory. */
   itkFactorylessNewMacro(Self)
   itkCloneMacro(Self)
 
   /** Runtime information support. */
   itkTypeMacro(DiffusionKurtosisReconstructionImageFilter, ImageToImageFilter)
 
   typedef TOutputPixelType                         OutputPixelType;
   typedef TInputPixelType                          InputPixelType;
 
   typedef typename Superclass::InputImageType      InputImageType;
   typedef Image< OutputPixelType, 3 >              OutputImageType;
 
   typedef itk::Image< short , 3>            MaskImageType;
 
   typedef typename Superclass::OutputImageRegionType   OutputImageRegionType;
 
 
   /** Holds each magnetic field gradient used to acquire one DWImage */
   typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType
                                                    GradientDirectionContainerType;
 
   // vector image typedefs for regularized fit
   typedef itk::VectorImage<float,3> VectorImageType;
   typedef itk::Image<itk::Vector<double, 3>, 3> InitialFitImageType;
 
   //-- input (Set) methods
   /** Set the initial solution for fitting, make sure the length and the values correspond to the parameters
     *   x0 = ( S_0, ADC_0, AKC_0 ) when also the S_0 is estimated
     *   x0 = ( ADC_0, AKC_0 ) when the S_0 is used in fitting
     */
   void SetInitialSolution(const vnl_vector<double>& x0 );
 
   /** Set whether the S_0 value is fitted or used in fitting */
   void SetOmitUnweightedValue( bool flag)
   { this->m_OmitBZero = flag; }
 
   /**
     Trigger a single computation of the Kurtosis values from the given input vector and the bvalues and returns the result as a KurtosisSnapshot object */
   KurtosisSnapshot GetSnapshot( const itk::VariableLengthVector< TInputPixelType > &input, vnl_vector<double> bvalues, bool omit_bzero);
 
   /**
     Trigger a single computation of the kurtosis values, first the bvalues vector is computed internally but then also stored into the returend snapshot */
   KurtosisSnapshot GetSnapshot( const itk::VariableLengthVector< TInputPixelType > &input, GradientDirectionContainerType::Pointer, float bvalue, bool omit_bzero);
 
   /**
     * Returns the value of the current data presented to the filter.
 
       If a mask is set, the voxels are first averaged before passed to the fitting procedure
     */
   KurtosisSnapshot GetCurrentSnapshot(bool omit_bzero);
 
   /** Set the reference bvalue of the input DW image */
   void SetReferenceBValue( double bvalue )
   { this->m_ReferenceBValue = bvalue; }
 
   /** Set the gradient directions */
   void SetGradientDirections( GradientDirectionContainerType::Pointer gradients );
 
   /** Restrict map generation to an image region */
   void SetMapOutputRegion( OutputImageRegionType region )
   {
     m_MapOutputRegion = region;
     this->m_ApplyPriorSmoothing = true;
   }
 
   void SetImageMask( MaskImageType::Pointer mask );
 
   /** Set smoothing sigma (default = 1.5 ), automatically enables smoothing prior to fitting */
   void SetSmoothingSigma( double sigma )
   {
     this->m_SmoothingSigma = sigma;
     this->m_ApplyPriorSmoothing = true;
   }
 
   /** Activate/Deactivate the gaussian smoothing applied to the input prior to fitting ( default = off ) */
   void SetUseSmoothingPriorToFitting( bool flag)
   {
     this->m_ApplyPriorSmoothing = flag;
   }
 
+  /** Set boundaries enforced by penalty terms in the fitting procedure */
+  void SetBoundariesForKurtosis( double lower, double upper )
+  {
+    m_KurtosisBounds[0] = lower; m_KurtosisBounds[1] = upper;
+  }
+
+  /** Exclude measurements associated with b-values higher than max_bvalue from fitting */
+  void SetMaximalBValueUsedForFitting( double max_bvalue )
+  {
+    m_MaxFitBValue = max_bvalue;
+  }
+
+  /** Select the method used in fitting of the data
+
+    STRAIHT - fit the exponential signal equation S / S_0 = exp [ ... ]
+    LOGARITHMIC - fit the logarithmic signal equation ln( S / S_0 ) = []
+    */
+  void SetFittingScale( FitScale scale )
+  {
+    m_ScaleForFitting = scale;
+  }
+
 protected:
   DiffusionKurtosisReconstructionImageFilter();
   virtual ~DiffusionKurtosisReconstructionImageFilter() {}
 
   void GenerateOutputInformation() override;
 
   void AfterThreadedGenerateData() override;
 
   void BeforeThreadedGenerateData() override;
 
   void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId) override;
 
   double m_ReferenceBValue;
 
   vnl_vector<double> m_BValues;
 
   vnl_vector<double> m_InitialPosition;
 
   bool m_OmitBZero;
 
   OutputImageRegionType m_MapOutputRegion;
 
   MaskImageType::Pointer m_MaskImage;
 
   typename InputImageType::Pointer m_ProcessedInputImage;
 
   bool m_ApplyPriorSmoothing;
   double m_SmoothingSigma;
 
+  vnl_vector_fixed<double, 2> m_KurtosisBounds;
+  double m_MaxFitBValue;
+
+  FitScale m_ScaleForFitting;
+
 private:
 
 
 };
 
 } //end namespace itk
 
 #ifndef ITK_MANUAL_INSTANTIATION
 #include "itkDiffusionKurtosisReconstructionImageFilter.cxx"
 #endif
 
 
 #endif // DIFFUSIONKURTOSISRECONSTRUCTIONIMAGEFILTER_H