diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp
index 7ede13e809..845c27ddfc 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp
@@ -1,781 +1,781 @@
 /*===================================================================
 
 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 __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_cpp
 #define __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_cpp
 
 #include "itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h"
 #include "itkImageRegionConstIterator.h"
 #include "itkImageRegionConstIteratorWithIndex.h"
 #include "itkImageRegionIterator.h"
 
 #include "vnl/vnl_matrix.h"
 #include "vnl/algo/vnl_symmetric_eigensystem.h"
 
 #include "itkRegularizedIVIMReconstructionFilter.h"
 
 #include <mitkLogMacros.h>
 
 #define IVIM_FOO -100000
 
 namespace itk {
 
 
 template< class TIn, class TOut>
 DiffusionIntravoxelIncoherentMotionReconstructionImageFilter<TIn, TOut>
 ::DiffusionIntravoxelIncoherentMotionReconstructionImageFilter() :
     m_GradientDirectionContainer(NULL),
     m_Method(IVIM_DSTAR_FIX),
     m_FitDStar(true),
     m_Verbose(false)
 {
     this->SetNumberOfRequiredInputs( 1 );
 
     this->SetNumberOfRequiredOutputs( 3 );
     typename OutputImageType::Pointer outputPtr1 = OutputImageType::New();
     this->SetNthOutput(0, outputPtr1.GetPointer());
     typename OutputImageType::Pointer outputPtr2 = OutputImageType::New();
     this->SetNthOutput(1, outputPtr2.GetPointer());
     typename OutputImageType::Pointer outputPtr3 = OutputImageType::New();
     this->SetNthOutput(2, outputPtr3.GetPointer());
 }
 
 
 
 template< class TIn, class TOut>
 void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter<TIn, TOut>
 ::BeforeThreadedGenerateData()
 {
 
     // Input must be an itk::VectorImage.
     std::string gradientImageClassName(
                 this->ProcessObject::GetInput(0)->GetNameOfClass());
     if ( strcmp(gradientImageClassName.c_str(),"VectorImage") != 0 )
     {
         itkExceptionMacro( <<
                            "There is only one Gradient image. I expect that to be a VectorImage. "
                            << "But its of type: " << gradientImageClassName );
     }
 
     // Compute the indicies of the baseline images and gradient images
     // If no b=0 mm/s² gradients ar found, the next lowest b-value is used.
     GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
     double minNorm = itk::NumericTraits<double>::max();
     while( gdcit != this->m_GradientDirectionContainer->End() )
     {
         double norm = gdcit.Value().one_norm();
         if (norm<minNorm)
             minNorm = norm;
         ++gdcit;
     }
     minNorm += 0.001;
     gdcit = this->m_GradientDirectionContainer->Begin();
     while( gdcit != this->m_GradientDirectionContainer->End() )
     {
         if(gdcit.Value().one_norm() <= minNorm)
         {
             m_Snap.baselineind.push_back(gdcit.Index());
         }
         else
         {
             m_Snap.gradientind.push_back(gdcit.Index());
             double twonorm = gdcit.Value().two_norm();
             m_Snap.bvals.push_back( m_BValue*twonorm*twonorm );
         }
         ++gdcit;
     }
     if (m_Snap.gradientind.size()==0)
         itkExceptionMacro("Only one b-value supplied. At least two needed for IVIM fit.");
 
     // check sind die grad und base gleichlang? alle grad gerade und base ungerade? dann iterierende aufnahme!!
     m_Snap.iterated_sequence = false;
     if(m_Snap.baselineind.size() == m_Snap.gradientind.size())
     {
         int size = m_Snap.baselineind.size();
         int sum_b = 0, sum_g = 0;
         for(int i=0; i<size; i++)
         {
             sum_b += m_Snap.baselineind.at(i) % 2;
             sum_g += m_Snap.gradientind.at(i) % 2;
         }
         if(  (sum_b == size || sum_b == 0)
              && (sum_g == size || sum_g == 0) )
         {
             m_Snap.iterated_sequence = true;
             if(m_Verbose)
             {
                 MITK_INFO << "Iterating b0 and diffusion weighted aquisition detected. Weighting each weighted measurement with its own b0.";
             }
         }
     }
 
     // number of measurements
     m_Snap.N = m_Snap.gradientind.size();
 
     // bvalue array
     m_Snap.bvalues.set_size(m_Snap.N);
     for(int i=0; i<m_Snap.N; i++)
     {
         m_Snap.bvalues[i] = m_Snap.bvals.at(i);
     }
 
     if(m_Verbose)
     {
         std::cout << "ref bval: " << m_BValue << "; b-values: ";
         for(int i=0; i<m_Snap.N; i++)
         {
             std::cout << m_Snap.bvalues[i] << "; ";
         }
         std::cout << std::endl;
     }
 
     // extract bvals higher than threshold
     if(m_Method == IVIM_D_THEN_DSTAR || m_Method == IVIM_LINEAR_D_THEN_F || m_Method == IVIM_REGULARIZED)
     {
         for(int i=0; i<m_Snap.N; i++)
         {
             if(m_Snap.bvalues[i]>m_BThres)
             {
                 m_Snap.high_indices.push_back(i);
             }
         }
     }
     m_Snap.Nhigh = m_Snap.high_indices.size();
     m_Snap.high_bvalues.set_size(m_Snap.Nhigh);
     m_Snap.high_meas.set_size(m_Snap.Nhigh);
     for(int i=0; i<m_Snap.Nhigh; i++)
     {
         m_Snap.high_bvalues[i] = m_Snap.bvalues[m_Snap.high_indices.at(i)];
     }
 
 }
 
 template< class TIn, class TOut>
 MeasAndBvals DiffusionIntravoxelIncoherentMotionReconstructionImageFilter<TIn, TOut>
 ::ApplyS0Threshold(vnl_vector<double> &meas, vnl_vector<double> &bvals)
 {
     std::vector<double> newmeas;
     std::vector<double> newbvals;
 
     int N = meas.size();
     for(int i=0; i<N; i++)
     {
         if( meas[i] != IVIM_FOO )
         {
             newmeas.push_back(meas[i]);
             newbvals.push_back(bvals[i]);
         }
     }
 
     MeasAndBvals retval;
 
     retval.N = newmeas.size();
 
     retval.meas.set_size(retval.N);
-    for(int i=0; i<newmeas.size(); i++)
+    for(size_t i=0; i<newmeas.size(); i++)
     {
         retval.meas[i] = newmeas[i];
     }
 
     retval.bvals.set_size(retval.N);
-    for(int i=0; i<newbvals.size(); i++)
+    for(size_t i=0; i<newbvals.size(); i++)
     {
         retval.bvals[i] = newbvals[i];
     }
 
     return retval;
 }
 
 template< class TIn, class TOut>
 void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter<TIn, TOut>
 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
                        ThreadIdType )
 {
 
     typename OutputImageType::Pointer outputImage =
             static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput());
     ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread);
     oit.GoToBegin();
 
     typename OutputImageType::Pointer dImage =
             static_cast< OutputImageType * >(this->ProcessObject::GetOutput(1));
     ImageRegionIterator< OutputImageType > oit1(dImage, outputRegionForThread);
     oit1.GoToBegin();
 
     typename OutputImageType::Pointer dstarImage =
             static_cast< OutputImageType * >(this->ProcessObject::GetOutput(2));
     ImageRegionIterator< OutputImageType > oit2(dstarImage, outputRegionForThread);
     oit2.GoToBegin();
 
     typedef ImageRegionConstIterator< InputImageType > InputIteratorType;
     typedef typename InputImageType::PixelType         InputVectorType;
     typename InputImageType::Pointer inputImagePointer = NULL;
 
     // Would have liked a dynamic_cast here, but seems SGI doesn't like it
     // The enum will DiffusionIntravoxelIncoherentMotionReconstructionImageFilterensure that an inappropriate cast is not done
     inputImagePointer = static_cast< InputImageType * >(
                 this->ProcessObject::GetInput(0) );
 
     InputIteratorType iit(inputImagePointer, outputRegionForThread );
     iit.GoToBegin();
 
     // init internal vector image for regularized fit
     m_InternalVectorImage = VectorImageType::New();
     m_InternalVectorImage->SetSpacing( inputImagePointer->GetSpacing() );   // Set the image spacing
     m_InternalVectorImage->SetOrigin( inputImagePointer->GetOrigin() );     // Set the image origin
     m_InternalVectorImage->SetDirection( inputImagePointer->GetDirection() );  // Set the image direction
     m_InternalVectorImage->SetRegions( inputImagePointer->GetLargestPossibleRegion() );
 
     m_InitialFitImage = InitialFitImageType::New();
     m_InitialFitImage->SetSpacing( inputImagePointer->GetSpacing() );   // Set the image spacing
     m_InitialFitImage->SetOrigin( inputImagePointer->GetOrigin() );     // Set the image origin
     m_InitialFitImage->SetDirection( inputImagePointer->GetDirection() );  // Set the image direction
     m_InitialFitImage->SetRegions( inputImagePointer->GetLargestPossibleRegion() );
 
     if(m_Method == IVIM_REGULARIZED)
     {
         m_InternalVectorImage->SetVectorLength(m_Snap.Nhigh);
         m_InternalVectorImage->Allocate();
         VectorImageType::PixelType varvec(m_Snap.Nhigh);
         for(int i=0; i<m_Snap.Nhigh; i++) varvec[i] = IVIM_FOO;
         m_InternalVectorImage->FillBuffer(varvec);
 
         m_InitialFitImage->Allocate();
         InitialFitImageType::PixelType vec;
         vec[0] = 0.5; vec[1] = 0.01; vec[2]=0.001;
         m_InitialFitImage->FillBuffer(vec);
     }
 
     typedef itk::ImageRegionIterator<VectorImageType> VectorIteratorType;
     VectorIteratorType vecit(m_InternalVectorImage, outputRegionForThread );
     vecit.GoToBegin();
 
     typedef itk::ImageRegionIterator<InitialFitImageType> InitIteratorType;
     InitIteratorType initit(m_InitialFitImage, outputRegionForThread );
     initit.GoToBegin();
 
     while( !iit.IsAtEnd() )
     {
         InputVectorType measvec = iit.Get();
 
         typename NumericTraits<InputPixelType>::AccumulateType b0 = NumericTraits<InputPixelType>::Zero;
 
         m_Snap.meas.set_size(m_Snap.N);
         m_Snap.allmeas.set_size(m_Snap.N);
         if(!m_Snap.iterated_sequence)
         {
             // Average the baseline image pixels
             for(unsigned int i = 0; i < m_Snap.baselineind.size(); ++i)
             {
                 b0 += measvec[m_Snap.baselineind[i]];
             }
             if(m_Snap.baselineind.size())
                 b0 /= m_Snap.baselineind.size();
 
             // measurement vector
             for(int i = 0; i < m_Snap.N; ++i)
             {
                 m_Snap.allmeas[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001);
 
                 if(measvec[m_Snap.gradientind[i]] > m_S0Thres)
                 {
                     m_Snap.meas[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001);
                 }
                 else
                 {
                     m_Snap.meas[i] = IVIM_FOO;
                 }
             }
         }
         else
         {
             // measurement vector
             for(int i = 0; i < m_Snap.N; ++i)
             {
                 b0 = measvec[m_Snap.baselineind[i]];
 
                 m_Snap.allmeas[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001);
 
                 if(measvec[m_Snap.gradientind[i]] > m_S0Thres)
                 {
                     m_Snap.meas[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001);
                 }
                 else
                 {
                     m_Snap.meas[i] = IVIM_FOO;
                 }
             }
         }
 
         m_Snap.currentF = 0;
         m_Snap.currentD = 0;
         m_Snap.currentDStar = 0;
 
         switch(m_Method)
         {
 
         case IVIM_D_THEN_DSTAR:
         {
 
             for(int i=0; i<m_Snap.Nhigh; i++)
             {
                 m_Snap.high_meas[i] = m_Snap.meas[m_Snap.high_indices.at(i)];
             }
 
             MeasAndBvals input = ApplyS0Threshold(m_Snap.high_meas, m_Snap.high_bvalues);
             m_Snap.bvals1 = input.bvals;
             m_Snap.meas1 = input.meas;
 
             if (input.N < 2)
             {
                 if (input.N == 1)
                 {
                     m_Snap.currentD = - log(input.meas[0]) / input.bvals[0];
                     m_Snap.currentF = 0;
                     m_Snap.currentDStar = 0;
                 }
                 break;
             }
 
             IVIM_d_and_f f_donly(input.N);
             f_donly.set_bvalues(input.bvals);
             f_donly.set_measurements(input.meas);
 
             vnl_vector< double > x_donly(2);
             x_donly[0] = 0.001;
             x_donly[1] = 0.1;
             // f 0.1 Dstar 0.01 D 0.001
 
             vnl_levenberg_marquardt lm_donly(f_donly);
             lm_donly.set_f_tolerance(0.0001);
             lm_donly.minimize(x_donly);
             m_Snap.currentD = x_donly[0];
             m_Snap.currentF = x_donly[1];
 
 
             if(m_FitDStar)
             {
                 MeasAndBvals input2 = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues);
                 m_Snap.bvals2 = input2.bvals;
                 m_Snap.meas2 = input2.meas;
                 if (input2.N < 2) break;
 
                 IVIM_dstar_only f_dstar_only(input2.N,m_Snap.currentD,m_Snap.currentF);
                 f_dstar_only.set_bvalues(input2.bvals);
                 f_dstar_only.set_measurements(input2.meas);
 
                 vnl_vector< double > x_dstar_only(1);
                 vnl_vector< double > fx_dstar_only(input2.N);
 
                 double opt = 1111111111111111.0;
                 int opt_idx = -1;
                 int num_its = 100;
                 double min_val = .001;
                 double max_val = .15;
                 for(int i=0; i<num_its; i++)
                 {
                     x_dstar_only[0] = min_val + i * ((max_val-min_val) / num_its);
                     f_dstar_only.f(x_dstar_only, fx_dstar_only);
                     double err = fx_dstar_only.two_norm();
                     if(err<opt)
                     {
                         opt = err;
                         opt_idx = i;
                     }
                 }
 
                 m_Snap.currentDStar = min_val + opt_idx * ((max_val-min_val) / num_its);
                 //          IVIM_fixd f_fixd(input2.N,m_Snap.currentD);
                 //          f_fixd.set_bvalues(input2.bvals);
                 //          f_fixd.set_measurements(input2.meas);
 
                 //          vnl_vector< double > x_fixd(2);
                 //          x_fixd[0] = 0.1;
                 //          x_fixd[1] = 0.01;
                 //          // f 0.1 Dstar 0.01 D 0.001
 
                 //          vnl_levenberg_marquardt lm_fixd(f_fixd);
                 //          lm_fixd.set_f_tolerance(0.0001);
                 //          lm_fixd.minimize(x_fixd);
 
                 //          m_Snap.currentF = x_fixd[0];
                 //          m_Snap.currentDStar = x_fixd[1];
             }
 
             break;
         }
 
         case IVIM_DSTAR_FIX:
         {
             MeasAndBvals input = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues);
             m_Snap.bvals1 = input.bvals;
             m_Snap.meas1 = input.meas;
             if (input.N < 2) break;
 
             IVIM_fixdstar f_fixdstar(input.N,m_DStar);
             f_fixdstar.set_bvalues(input.bvals);
             f_fixdstar.set_measurements(input.meas);
 
             vnl_vector< double > x(2);
             x[0] = 0.1;
             x[1] = 0.001;
             // f 0.1 Dstar 0.01 D 0.001
 
             vnl_levenberg_marquardt lm(f_fixdstar);
             lm.set_f_tolerance(0.0001);
             lm.minimize(x);
 
             m_Snap.currentF = x[0];
             m_Snap.currentD = x[1];
             m_Snap.currentDStar = m_DStar;
 
             break;
         }
 
         case IVIM_FIT_ALL:
         {
 
             MeasAndBvals input = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues);
             m_Snap.bvals1 = input.bvals;
             m_Snap.meas1 = input.meas;
             if (input.N < 3) break;
 
             IVIM_3param f_3param(input.N);
             f_3param.set_bvalues(input.bvals);
             f_3param.set_measurements(input.meas);
 
             vnl_vector< double > x(3);
             x[0] = 0.1;
             x[1] = 0.001;
             x[2] = 0.01;
             // f 0.1 Dstar 0.01 D 0.001
 
             vnl_levenberg_marquardt lm(f_3param);
             lm.set_f_tolerance(0.0001);
             lm.minimize(x);
 
             m_Snap.currentF = x[0];
             m_Snap.currentD = x[1];
             m_Snap.currentDStar = x[2];
 
             break;
         }
 
         case IVIM_LINEAR_D_THEN_F:
         {
 
             //          // neglect zero-measurements
             //          bool zero = false;
             //          for(int i=0; i<Nhigh; i++)
             //          {
             //            if( meas[high_indices.at(i)] == 0 )
             //            {
             //              f=0;
             //              zero = true;
             //              break;
             //            }
             //          }
             //          if(zero) break;
 
             for(int i=0; i<m_Snap.Nhigh; i++)
             {
                 m_Snap.high_meas[i] = m_Snap.meas[m_Snap.high_indices.at(i)];
             }
 
             MeasAndBvals input = ApplyS0Threshold(m_Snap.high_meas, m_Snap.high_bvalues);
             m_Snap.bvals1 = input.bvals;
             m_Snap.meas1 = input.meas;
 
             if (input.N < 2)
             {
                 if (input.N == 1)
                 {
                     m_Snap.currentD = - log(input.meas[0]) / input.bvals[0];
                     m_Snap.currentF = 0;
                     m_Snap.currentDStar = 0;
                 }
                 break;
             }
 
             for(int i=0; i<input.N; i++)
             {
                 input.meas[i] = log(input.meas[i]);
             }
 
             double bval_m = 0;
             double meas_m = 0;
             for(int i=0; i<input.N; i++)
             {
                 bval_m += input.bvals[i];
                 meas_m += input.meas[i];
             }
             bval_m /= input.N;
             meas_m /= input.N;
 
             vnl_matrix<double> X(input.N,2);
             for(int i=0; i<input.N; i++)
             {
                 X(i,0) = input.bvals[i] - bval_m;
                 X(i,1) = input.meas[i] - meas_m;
             }
 
             vnl_matrix<double> XX = X.transpose() * X;
             vnl_symmetric_eigensystem<double> eigs(XX);
 
             vnl_vector<double> eig;
             if(eigs.get_eigenvalue(0) > eigs.get_eigenvalue(1))
                 eig = eigs.get_eigenvector(0);
             else
                 eig = eigs.get_eigenvector(1);
 
             m_Snap.currentF = 1 - exp( meas_m - bval_m*(eig(1)/eig(0)) );
             m_Snap.currentD = -eig(1)/eig(0);
 
             if(m_FitDStar)
             {
                 MeasAndBvals input2 = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues);
                 m_Snap.bvals2 = input2.bvals;
                 m_Snap.meas2 = input2.meas;
                 if (input2.N < 2) break;
 
                 IVIM_dstar_only f_dstar_only(input2.N,m_Snap.currentD,m_Snap.currentF);
                 f_dstar_only.set_bvalues(input2.bvals);
                 f_dstar_only.set_measurements(input2.meas);
 
                 vnl_vector< double > x_dstar_only(1);
                 vnl_vector< double > fx_dstar_only(input2.N);
 
                 double opt = 1111111111111111.0;
                 int opt_idx = -1;
                 int num_its = 100;
                 double min_val = .001;
                 double max_val = .15;
                 for(int i=0; i<num_its; i++)
                 {
                     x_dstar_only[0] = min_val + i * ((max_val-min_val) / num_its);
                     f_dstar_only.f(x_dstar_only, fx_dstar_only);
                     double err = fx_dstar_only.two_norm();
                     if(err<opt)
                     {
                         opt = err;
                         opt_idx = i;
                     }
                 }
 
                 m_Snap.currentDStar = min_val + opt_idx * ((max_val-min_val) / num_its);
             }
             // MITK_INFO << "choosing " << opt_idx << " => " << DStar;
             //          x_dstar_only[0] = 0.01;
             //          // f 0.1 Dstar 0.01 D 0.001
 
             //          vnl_levenberg_marquardt lm_dstar_only(f_dstar_only);
             //          lm_dstar_only.set_f_tolerance(0.0001);
             //          lm_dstar_only.minimize(x_dstar_only);
 
             //          DStar = x_dstar_only[0];
 
             break;
         }
 
         case IVIM_REGULARIZED:
         {
             //m_Snap.high_meas, m_Snap.high_bvalues;
             for(int i=0; i<m_Snap.Nhigh; i++)
             {
                 m_Snap.high_meas[i] = m_Snap.meas[m_Snap.high_indices.at(i)];
             }
 
             MeasAndBvals input = ApplyS0Threshold(m_Snap.high_meas, m_Snap.high_bvalues);
 
             vnl_vector< double > x_donly(2);
             x_donly[0] = 0.001;
             x_donly[1] = 0.1;
 
             if(input.N >= 2)
             {
                 IVIM_d_and_f f_donly(input.N);
                 f_donly.set_bvalues(input.bvals);
                 f_donly.set_measurements(input.meas);
                 //MITK_INFO << "initial fit N=" << input.N << ", min-b = " << input.bvals[0] << ", max-b = " << input.bvals[input.N-1];
                 vnl_levenberg_marquardt lm_donly(f_donly);
                 lm_donly.set_f_tolerance(0.0001);
                 lm_donly.minimize(x_donly);
             }
 
             typename InitialFitImageType::PixelType initvec;
             initvec[0] = x_donly[1];
             initvec[1] = x_donly[0];
             initit.Set(initvec);
             //MITK_INFO << "Init vox " << initit.GetIndex() << " with " << initvec[0] << "; " << initvec[1];
             ++initit;
 
             int N = m_Snap.high_meas.size();
             typename VectorImageType::PixelType vec(N);
             for(int i=0; i<N; i++)
             {
                 vec[i] = m_Snap.high_meas[i];
                 //MITK_INFO << "vec" << i << " = " << m_Snap.high_meas[i];
             }
 
             vecit.Set(vec);
             ++vecit;
 
             if(!m_Verbose)
             {
                 // report the middle voxel
                 if(  vecit.GetIndex()[0] == m_CrossPosition[0]
                      && vecit.GetIndex()[0] == m_CrossPosition[1]
                      && vecit.GetIndex()[0] == m_CrossPosition[2] )
                 {
                     MeasAndBvals input = ApplyS0Threshold(m_Snap.high_meas, m_Snap.high_bvalues);
                     m_Snap.bvals1 = input.bvals;
                     m_Snap.meas1 = input.meas;
 
                     MeasAndBvals input2 = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues);
                     m_Snap.bvals2 = input2.bvals;
                     m_Snap.meas2 = input2.meas;
 
                     m_tmp_allmeas = m_Snap.allmeas;
                 }
             }
 
             break;
         }
         }
         m_Snap.currentFunceiled = m_Snap.currentF;
         IVIM_CEIL( m_Snap.currentF, 0.0, 1.0 );
 
         oit.Set( m_Snap.currentF );
         oit1.Set( m_Snap.currentD );
         oit2.Set( m_Snap.currentDStar );
 
         // std::cout << "\tf=" << x[0] << "\tD=" << x[1] << " ; "<<std::endl;
 
         ++oit;
         ++oit1;
         ++oit2;
         ++iit;
     }
 
     if(m_Verbose)
     {
         std::cout << "One Thread finished reconstruction" << std::endl;
     }
 }
 
 template< class TIn, class TOut>
 void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter<TIn, TOut>
 ::AfterThreadedGenerateData()
 {
     if(m_Method == IVIM_REGULARIZED)
     {
         typename OutputImageType::Pointer outputImage =
                 static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput());
         ImageRegionIterator< OutputImageType > oit0(outputImage, outputImage->GetLargestPossibleRegion());
         oit0.GoToBegin();
 
         typename OutputImageType::Pointer dImage =
                 static_cast< OutputImageType * >(this->ProcessObject::GetOutput(1));
         ImageRegionIterator< OutputImageType > oit1(dImage, dImage->GetLargestPossibleRegion());
         oit1.GoToBegin();
 
         typename OutputImageType::Pointer dstarImage =
                 static_cast< OutputImageType * >(this->ProcessObject::GetOutput(2));
         ImageRegionIterator< OutputImageType > oit2(dstarImage, dstarImage->GetLargestPossibleRegion());
         oit2.GoToBegin();
 
         typedef itk::RegularizedIVIMReconstructionFilter
                 <double,double,float> RegFitType;
         RegFitType::Pointer filter = RegFitType::New();
         filter->SetInput(m_InitialFitImage);
         filter->SetReferenceImage(m_InternalVectorImage);
         filter->SetBValues(m_Snap.high_bvalues);
         filter->SetNumberIterations(m_NumberIterations);
         filter->SetNumberOfThreads(1);
         filter->SetLambda(m_Lambda);
         filter->Update();
         typename RegFitType::OutputImageType::Pointer outimg = filter->GetOutput();
 
         ImageRegionConstIterator< RegFitType::OutputImageType > iit(outimg, outimg->GetLargestPossibleRegion());
         iit.GoToBegin();
 
         while( !iit.IsAtEnd() )
         {
             double f = iit.Get()[0];
             IVIM_CEIL( f, 0.0, 1.0 );
 
             oit0.Set( myround(f * 100.0) );
             oit1.Set( myround(iit.Get()[1] * 10000.0) );
             oit2.Set( myround(iit.Get()[2] * 1000.0) );
 
             if(!m_Verbose)
             {
                 // report the middle voxel
                 if(  iit.GetIndex()[0] == m_CrossPosition[0]
                      && iit.GetIndex()[1] == m_CrossPosition[1]
                      && iit.GetIndex()[2] == m_CrossPosition[2] )
                 {
                     m_Snap.currentF = f;
                     m_Snap.currentD = iit.Get()[1];
                     m_Snap.currentDStar = iit.Get()[2];
                     m_Snap.allmeas = m_tmp_allmeas;
                     MITK_INFO << "setting " << f << ";" << iit.Get()[1] << ";" << iit.Get()[2];
                 }
             }
 
             ++oit0;
             ++oit1;
             ++oit2;
             ++iit;
         }
     }
 }
 
 template< class TIn, class TOut>
 double DiffusionIntravoxelIncoherentMotionReconstructionImageFilter<TIn, TOut>
 ::myround(double number)
 {
     return number < 0.0 ? ceil(number - 0.5) : floor(number + 0.5);
 }
 
 template< class TIn, class TOut>
 void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter<TIn, TOut>
 ::SetGradientDirections( GradientDirectionContainerType *gradientDirection )
 {
     this->m_GradientDirectionContainer = gradientDirection;
     this->m_NumberOfGradientDirections = gradientDirection->Size();
 }
 
 template< class TIn, class TOut>
 void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter<TIn, TOut>
 ::PrintSelf(std::ostream& os, Indent indent) const
 {
     Superclass::PrintSelf(os,indent);
 
     if ( m_GradientDirectionContainer )
     {
         os << indent << "GradientDirectionContainer: "
            << m_GradientDirectionContainer << std::endl;
     }
     else
     {
         os << indent <<
               "GradientDirectionContainer: (Gradient directions not set)" << std::endl;
     }
 }
 
 }
 
 #endif // __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_cpp
diff --git a/Modules/DiffusionImaging/Quantification/Algorithms/itkRegularizedIVIMLocalVariationImageFilter.h b/Modules/DiffusionImaging/Quantification/Algorithms/itkRegularizedIVIMLocalVariationImageFilter.h
index ded47c327a..af2d0217c5 100644
--- a/Modules/DiffusionImaging/Quantification/Algorithms/itkRegularizedIVIMLocalVariationImageFilter.h
+++ b/Modules/DiffusionImaging/Quantification/Algorithms/itkRegularizedIVIMLocalVariationImageFilter.h
@@ -1,152 +1,152 @@
 /*===================================================================
 
 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 __itkRegularizedIVIMLocalVariationImageFilter_h
 #define __itkRegularizedIVIMLocalVariationImageFilter_h
 
 #include "itkImageToImageFilter.h"
 #include "itkImage.h"
 
 namespace itk
 {
 
   template<class TPixelType>
   class IVIMSquaredEuclideanMetric
   {
   public:
     static double Calc(TPixelType p)
     {
       return p*p;
     }
   };
 
   template<>
   class IVIMSquaredEuclideanMetric<itk::Vector<double,3> >
   {
   public:
     static double Calc(itk::Vector<double,3> p)
     {
       return p[1]*p[1];
     }
   };
 
   template<>
   class IVIMSquaredEuclideanMetric<itk::VariableLengthVector<float> >
   {
   public:
     static double Calc(itk::VariableLengthVector<float> p)
     {
       return p.GetSquaredNorm();
     }
   };
 
   template<>
   class IVIMSquaredEuclideanMetric<itk::VariableLengthVector<double> >
   {
   public:
     static double Calc(itk::VariableLengthVector<double> p)
     {
       return p.GetSquaredNorm();
     }
   };
 
 /** \class RegularizedIVIMLocalVariationImageFilter
  * \brief Calculates the local variation in each pixel
  *
  * Reference: Tony F. Chan et al., The digital TV filter and nonlinear denoising
  *
  * \sa Image
  * \sa Neighborhood
  * \sa NeighborhoodOperator
  * \sa NeighborhoodIterator
  *
  * \ingroup IntensityImageFilters
  */
 template <class TInputImage, class TOutputImage>
 class RegularizedIVIMLocalVariationImageFilter :
     public ImageToImageFilter< TInputImage, TOutputImage >
 {
 public:
   /** Extract dimension from input and output image. */
   itkStaticConstMacro(InputImageDimension, unsigned int,
                       TInputImage::ImageDimension);
   itkStaticConstMacro(OutputImageDimension, unsigned int,
                       TOutputImage::ImageDimension);
 
   /** Convenient typedefs for simplifying declarations. */
   typedef TInputImage InputImageType;
   typedef TOutputImage OutputImageType;
 
   /** Standard class typedefs. */
   typedef RegularizedIVIMLocalVariationImageFilter Self;
   typedef ImageToImageFilter< InputImageType, OutputImageType> Superclass;
   typedef SmartPointer<Self> Pointer;
   typedef SmartPointer<const Self>  ConstPointer;
 
   /** Method for creation through the object factory. */
   itkNewMacro(Self);
 
   /** Run-time type information (and related methods). */
   itkTypeMacro(RegularizedIVIMLocalVariationImageFilter, ImageToImageFilter);
 
   /** Image typedef support. */
   typedef typename InputImageType::PixelType InputPixelType;
   typedef typename OutputImageType::PixelType OutputPixelType;
 
   typedef typename InputImageType::RegionType InputImageRegionType;
   typedef typename OutputImageType::RegionType OutputImageRegionType;
 
   typedef typename InputImageType::SizeType InputSizeType;
 
   /** MedianImageFilter needs a larger input requested region than
    * the output requested region.  As such, MedianImageFilter needs
    * to provide an implementation for GenerateInputRequestedRegion()
    * in order to inform the pipeline execution model.
    *
    * \sa ImageToImageFilter::GenerateInputRequestedRegion() */
   virtual void GenerateInputRequestedRegion()
     throw(InvalidRequestedRegionError);
 
 protected:
   RegularizedIVIMLocalVariationImageFilter();
   virtual ~RegularizedIVIMLocalVariationImageFilter() {}
   void PrintSelf(std::ostream& os, Indent indent) const;
 
   /** MedianImageFilter can be implemented as a multithreaded filter.
    * Therefore, this implementation provides a ThreadedGenerateData()
    * routine which is called for each processing thread. The output
    * image data is allocated automatically by the superclass prior to
    * calling ThreadedGenerateData().  ThreadedGenerateData can only
    * write to the portion of the output image specified by the
    * parameter "outputRegionForThread"
    *
    * \sa ImageToImageFilter::ThreadedGenerateData(),
    *     ImageToImageFilter::GenerateData() */
   void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
-                            int threadId );
+                            ThreadIdType threadId );
 
 private:
   RegularizedIVIMLocalVariationImageFilter(const Self&); //purposely not implemented
   void operator=(const Self&); //purposely not implemented
 };
 
 } // end namespace itk
 
 #ifndef ITK_MANUAL_INSTANTIATION
 #include "itkRegularizedIVIMLocalVariationImageFilter.txx"
 #endif
 
 #endif //RegularizedIVIMLocalVariationImageFilter
diff --git a/Modules/DiffusionImaging/Quantification/Algorithms/itkRegularizedIVIMLocalVariationImageFilter.txx b/Modules/DiffusionImaging/Quantification/Algorithms/itkRegularizedIVIMLocalVariationImageFilter.txx
index 54acf9f02f..3acd73b82a 100644
--- a/Modules/DiffusionImaging/Quantification/Algorithms/itkRegularizedIVIMLocalVariationImageFilter.txx
+++ b/Modules/DiffusionImaging/Quantification/Algorithms/itkRegularizedIVIMLocalVariationImageFilter.txx
@@ -1,192 +1,192 @@
 /*===================================================================
 
 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 _itkRegularizedIVIMLocalVariationImageFilter_txx
 #define _itkRegularizedIVIMLocalVariationImageFilter_txx
 
 #include "itkConstShapedNeighborhoodIterator.h"
 #include "itkNeighborhoodInnerProduct.h"
 #include "itkImageRegionIterator.h"
 #include "itkNeighborhoodAlgorithm.h"
 #include "itkZeroFluxNeumannBoundaryCondition.h"
 #include "itkOffset.h"
 #include "itkProgressReporter.h"
 #include "itkVectorImage.h"
 
 #include <vector>
 #include <algorithm>
 
 namespace itk
 {
 
   template <class TInputImage, class TOutputImage>
   RegularizedIVIMLocalVariationImageFilter<TInputImage, TOutputImage>
     ::RegularizedIVIMLocalVariationImageFilter()
   {}
 
   template <class TInputImage, class TOutputImage>
   void
     RegularizedIVIMLocalVariationImageFilter<TInputImage, TOutputImage>
     ::GenerateInputRequestedRegion() throw (InvalidRequestedRegionError)
   {
     // call the superclass' implementation of this method
     Superclass::GenerateInputRequestedRegion();
 
     // get pointers to the input and output
     typename Superclass::InputImagePointer inputPtr =
       const_cast< TInputImage * >( this->GetInput() );
     typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
 
     if ( !inputPtr || !outputPtr )
     {
       return;
     }
 
     // get a copy of the input requested region (should equal the output
     // requested region)
     typename TInputImage::RegionType inputRequestedRegion;
     inputRequestedRegion = inputPtr->GetRequestedRegion();
 
     // pad the input requested region by 1
     inputRequestedRegion.PadByRadius( 1 );
 
     // crop the input requested region at the input's largest possible region
     if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
     {
       inputPtr->SetRequestedRegion( inputRequestedRegion );
       return;
     }
     else
     {
       // Couldn't crop the region (requested region is outside the largest
       // possible region).  Throw an exception.
 
       // store what we tried to request (prior to trying to crop)
       inputPtr->SetRequestedRegion( inputRequestedRegion );
 
       // build an exception
       InvalidRequestedRegionError e(__FILE__, __LINE__);
       e.SetLocation(ITK_LOCATION);
       e.SetDescription("Requested region outside possible region.");
       e.SetDataObject(inputPtr);
       throw e;
     }
   }
 
   template< class TInputImage, class TOutputImage>
   void
     RegularizedIVIMLocalVariationImageFilter< TInputImage, TOutputImage>
     ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
-    int threadId)
+    ThreadIdType threadId)
   {
 
     // Allocate output
     typename OutputImageType::Pointer output = this->GetOutput();
     typename  InputImageType::ConstPointer input  = this->GetInput();
 
     itk::Size<InputImageDimension> size;
     for( int i=0; i<InputImageDimension; i++)
       size[i] = 1;
 
     // Find the data-set boundary "faces"
     NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType> bC;
     typename NeighborhoodAlgorithm::
       ImageBoundaryFacesCalculator<InputImageType>::FaceListType
       faceList = bC(input, outputRegionForThread, size);
 
     // support progress methods/callbacks
     ProgressReporter progress(
       this, threadId, outputRegionForThread.GetNumberOfPixels());
 
     ZeroFluxNeumannBoundaryCondition<InputImageType> nbc;
     std::vector<InputPixelType> pixels;
 
     // Process each of the boundary faces.  These are N-d regions which border
     // the edge of the buffer.
     for ( typename NeighborhoodAlgorithm::
       ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator
       fit=faceList.begin(); fit != faceList.end(); ++fit)
     {
 
       // iterators over output and input
       ImageRegionIterator<OutputImageType>
         output_image_it(output, *fit);
       ImageRegionConstIterator<InputImageType>
         input_image_it(input.GetPointer(), *fit);
 
       // neighborhood iterator for input image
       ConstShapedNeighborhoodIterator<InputImageType>
         input_image_neighbors_it(size, input, *fit);
       typename ConstShapedNeighborhoodIterator<InputImageType>::
         OffsetType offset;
       input_image_neighbors_it.OverrideBoundaryCondition(&nbc);
       input_image_neighbors_it.ClearActiveList();
       for(int i=0; i<InputImageDimension; i++)
       {
         offset.Fill(0);
         offset[i] = -1;
         input_image_neighbors_it.ActivateOffset(offset);
         offset[i] = 1;
         input_image_neighbors_it.ActivateOffset(offset);
       }
       input_image_neighbors_it.GoToBegin();
       //const unsigned int neighborhoodSize = InputImageDimension*2;
 
       while ( ! input_image_neighbors_it.IsAtEnd() )
       {
         // collect all the pixels in the neighborhood, note that we use
         // GetPixel on the NeighborhoodIterator to honor the boundary conditions
         typename OutputImageType::PixelType locVariation = 0;
         typename ConstShapedNeighborhoodIterator<InputImageType>::
           ConstIterator input_neighbors_it;
         for (input_neighbors_it = input_image_neighbors_it.Begin();
           ! input_neighbors_it.IsAtEnd();
           input_neighbors_it++)
         {
             typename TInputImage::PixelType diffVec =
               input_neighbors_it.Get()-input_image_it.Get();
             locVariation += IVIMSquaredEuclideanMetric
               <typename TInputImage::PixelType>::Calc(diffVec);
         }
         locVariation = sqrt(locVariation + 0.0001);
         output_image_it.Set(locVariation);
 
         // update iterators
         ++input_image_neighbors_it;
         ++output_image_it;
         ++input_image_it;
 
         // report progress
         progress.CompletedPixel();
       }
     }
   }
 
   /**
   * Standard "PrintSelf" method
   */
   template <class TInputImage, class TOutput>
   void
     RegularizedIVIMLocalVariationImageFilter<TInputImage, TOutput>
     ::PrintSelf(
     std::ostream& os,
     Indent indent) const
   {
     Superclass::PrintSelf( os, indent );
   }
 
 } // end namespace itk
 
 #endif //_itkRegularizedIVIMLocalVariationImageFilter_txx
diff --git a/Modules/DiffusionImaging/Quantification/Algorithms/itkRegularizedIVIMReconstructionSingleIteration.h b/Modules/DiffusionImaging/Quantification/Algorithms/itkRegularizedIVIMReconstructionSingleIteration.h
index 42701e4c20..3e5dbc7faf 100644
--- a/Modules/DiffusionImaging/Quantification/Algorithms/itkRegularizedIVIMReconstructionSingleIteration.h
+++ b/Modules/DiffusionImaging/Quantification/Algorithms/itkRegularizedIVIMReconstructionSingleIteration.h
@@ -1,141 +1,141 @@
 /*===================================================================
 
 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 __itkRegularizedIVIMReconstructionSingleIteration_h
 #define __itkRegularizedIVIMReconstructionSingleIteration_h
 
 #include "itkImageToImageFilter.h"
 #include "itkImage.h"
 #include "itkVectorImage.h"
 
 namespace itk
 {
 /** \class RegularizedIVIMReconstructionSingleIteration
  * \brief Applies a total variation denoising filter to an image
  *
  * Reference: Tony F. Chan et al., The digital TV filter and nonlinear denoising
  *
  * \sa Image
  * \sa Neighborhood
  * \sa NeighborhoodOperator
  * \sa NeighborhoodIterator
  *
  * \ingroup IntensityImageFilters
  */
 template <class TInputPixel, class TOutputPixel, class TRefPixelType>
 class RegularizedIVIMReconstructionSingleIteration :
     public ImageToImageFilter< itk::Image<itk::Vector<TInputPixel,3>, 3>, itk::Image<itk::Vector<TOutputPixel,3>, 3> >
 {
 public:
   /** Convenient typedefs for simplifying declarations. */
   typedef itk::Image<itk::Vector<TInputPixel,3>, 3> InputImageType;
   typedef itk::Image<itk::Vector<TOutputPixel,3>, 3> OutputImageType;
   typedef itk::VectorImage<TRefPixelType,3> RefImageType;
 
   /** Extract dimension from input and output image. */
   itkStaticConstMacro(InputImageDimension, unsigned int,
                       InputImageType::ImageDimension);
   itkStaticConstMacro(OutputImageDimension, unsigned int,
                       OutputImageType::ImageDimension);
 
   typedef itk::Image<float,InputImageDimension> LocalVariationImageType;
 
   /** Standard class typedefs. */
   typedef RegularizedIVIMReconstructionSingleIteration Self;
   typedef ImageToImageFilter< InputImageType, OutputImageType> Superclass;
   typedef SmartPointer<Self> Pointer;
   typedef SmartPointer<const Self>  ConstPointer;
 
   /** Method for creation through the object factory. */
   itkNewMacro(Self);
 
   /** Run-time type information (and related methods). */
   itkTypeMacro(RegularizedIVIMReconstructionSingleIteration, ImageToImageFilter);
 
   /** Image typedef support. */
   typedef typename InputImageType::PixelType InputPixelType;
   typedef typename OutputImageType::PixelType OutputPixelType;
 
   typedef typename InputImageType::RegionType InputImageRegionType;
   typedef typename OutputImageType::RegionType OutputImageRegionType;
 
   typedef typename InputImageType::SizeType InputSizeType;
 
   /** A larger input requested region than
    * the output requested region is required.
    * Therefore, an implementation for GenerateInputRequestedRegion()
    * is provided.
    *
    * \sa ImageToImageFilter::GenerateInputRequestedRegion() */
   virtual void GenerateInputRequestedRegion()
     throw(InvalidRequestedRegionError);
 
   itkSetMacro(Lambda, double);
   itkGetMacro(Lambda, double);
 
   void SetBValues(vnl_vector<double> bvals)
   { this->m_BValues = bvals; }
   vnl_vector<double> GetBValues()
   { return this->m_BValues; }
 
   void SetOriginalImage(RefImageType* in)
   { this->m_OriginalImage = in; }
   typename RefImageType::Pointer GetOriginialImage()
   { return this->m_OriginalImage; }
 
 protected:
   RegularizedIVIMReconstructionSingleIteration();
   virtual ~RegularizedIVIMReconstructionSingleIteration() {}
   void PrintSelf(std::ostream& os, Indent indent) const;
 
   /** MedianImageFilter can be implemented as a multithreaded filter.
    * Therefore, this implementation provides a ThreadedGenerateData()
    * routine which is called for each processing thread. The output
    * image data is allocated automatically by the superclass prior to
    * calling ThreadedGenerateData().  ThreadedGenerateData can only
    * write to the portion of the output image specified by the
    * parameter "outputRegionForThread"
    *
    * \sa ImageToImageFilter::ThreadedGenerateData(),
    *     ImageToImageFilter::GenerateData() */
   void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
-                            int threadId );
+                            ThreadIdType threadId );
 
   void BeforeThreadedGenerateData();
 
   typename LocalVariationImageType::Pointer m_LocalVariation;
 
   typename RefImageType::Pointer m_OriginalImage;
 
   double m_Lambda;
 
   vnl_vector<double> m_BValues;
 
 private:
 
   RegularizedIVIMReconstructionSingleIteration(const Self&);
 
   void operator=(const Self&);
 
 };
 
 } // end namespace itk
 
 #ifndef ITK_MANUAL_INSTANTIATION
 #include "itkRegularizedIVIMReconstructionSingleIteration.txx"
 #endif
 
 #endif //__itkRegularizedIVIMReconstructionSingleIteration__
diff --git a/Modules/DiffusionImaging/Quantification/Algorithms/itkRegularizedIVIMReconstructionSingleIteration.txx b/Modules/DiffusionImaging/Quantification/Algorithms/itkRegularizedIVIMReconstructionSingleIteration.txx
index 09d0e5d527..adc1d08ac1 100644
--- a/Modules/DiffusionImaging/Quantification/Algorithms/itkRegularizedIVIMReconstructionSingleIteration.txx
+++ b/Modules/DiffusionImaging/Quantification/Algorithms/itkRegularizedIVIMReconstructionSingleIteration.txx
@@ -1,310 +1,310 @@
 /*===================================================================
 
 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 _itkRegularizedIVIMReconstructionSingleIteration_txx
 #define _itkRegularizedIVIMReconstructionSingleIteration_txx
 
 // itk includes
 #include "itkConstShapedNeighborhoodIterator.h"
 #include "itkNeighborhoodInnerProduct.h"
 #include "itkImageRegionIterator.h"
 #include "itkNeighborhoodAlgorithm.h"
 #include "itkZeroFluxNeumannBoundaryCondition.h"
 #include "itkOffset.h"
 #include "itkProgressReporter.h"
 #include "itkRegularizedIVIMLocalVariationImageFilter.h"
 
 // other includes
 #include <vector>
 #include <algorithm>
 
 #define IVIM_FOO -100000
 
 namespace itk
 {
 
   /**
   * constructor
   */
   template <class TInputPixel, class TOutputPixel, class TRefPixelType>
   RegularizedIVIMReconstructionSingleIteration<TInputPixel, TOutputPixel, TRefPixelType>
     ::RegularizedIVIMReconstructionSingleIteration()
   {
     m_Lambda = 1.0;
     m_LocalVariation = LocalVariationImageType::New();
   }
 
   /**
   * generate requested region
   */
   template <class TInputPixel, class TOutputPixel, class TRefPixelType>
   void
   RegularizedIVIMReconstructionSingleIteration<TInputPixel, TOutputPixel, TRefPixelType>
     ::GenerateInputRequestedRegion() throw (InvalidRequestedRegionError)
   {
     // call the superclass' implementation of this method
     Superclass::GenerateInputRequestedRegion();
 
     // get pointers to the input and output
     typename Superclass::InputImagePointer inputPtr =
       const_cast< InputImageType * >( this->GetInput() );
     typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
 
     if ( !inputPtr || !outputPtr )
     {
       return;
     }
 
     // get a copy of the input requested region (should equal the output
     // requested region)
     typename InputImageType::RegionType inputRequestedRegion;
     inputRequestedRegion = inputPtr->GetRequestedRegion();
 
     // pad the input requested region by 1
     inputRequestedRegion.PadByRadius( 1 );
 
     // crop the input requested region at the input's largest possible region
     if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
     {
       inputPtr->SetRequestedRegion( inputRequestedRegion );
       return;
     }
     else
     {
       // Couldn't crop the region (requested region is outside the largest
       // possible region).  Throw an exception.
 
       // store what we tried to request (prior to trying to crop)
       inputPtr->SetRequestedRegion( inputRequestedRegion );
 
       // build an exception
       InvalidRequestedRegionError e(__FILE__, __LINE__);
       e.SetLocation(ITK_LOCATION);
       e.SetDescription("Requested region outside possible region.");
       e.SetDataObject(inputPtr);
       throw e;
     }
   }
 
 
    /**
   * generate output
   */
   template <class TInputPixel, class TOutputPixel, class TRefPixelType>
   void
   RegularizedIVIMReconstructionSingleIteration<TInputPixel, TOutputPixel, TRefPixelType>
     ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
-    int threadId)
+    ThreadIdType)
   {
 
     typename OutputImageType::Pointer output = this->GetOutput();
     typename  InputImageType::ConstPointer input  = this->GetInput();
 
     // Find the data-set boundary "faces"
     itk::Size<InputImageDimension> size;
     for( int i=0; i<InputImageDimension; i++)
       size[i] = 1;
 
     NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType> bC;
     typename NeighborhoodAlgorithm::
       ImageBoundaryFacesCalculator<InputImageType>::FaceListType
       faceList = bC(input, outputRegionForThread, size);
 
     NeighborhoodAlgorithm::
       ImageBoundaryFacesCalculator<LocalVariationImageType> lv_bC;
     typename NeighborhoodAlgorithm::
       ImageBoundaryFacesCalculator<LocalVariationImageType>::FaceListType
       lv_faceList = lv_bC(m_LocalVariation, outputRegionForThread, size);
 
     ZeroFluxNeumannBoundaryCondition<InputImageType> nbc;
     ZeroFluxNeumannBoundaryCondition<LocalVariationImageType> lv_nbc;
     std::vector<double> ws;
     std::vector<double> hs;
 
     typename NeighborhoodAlgorithm::
       ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator
       lv_fit=lv_faceList.begin();
 
     // Process each of the boundary faces.  These are N-d regions which border
     // the edge of the buffer.
     for ( typename NeighborhoodAlgorithm::
       ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator
       fit=faceList.begin(); fit != faceList.end(); ++fit)
     {
 
       // iterators over output, input, original and local variation image
       ImageRegionIterator<OutputImageType> output_image_it =
         ImageRegionIterator<OutputImageType>(output, *fit);
       ImageRegionConstIterator<InputImageType> input_image_it =
         ImageRegionConstIterator<InputImageType>(input, *fit);
       ImageRegionConstIterator<RefImageType> orig_image_it =
         ImageRegionConstIterator<RefImageType>(m_OriginalImage, *fit);
       ImageRegionConstIterator<LocalVariationImageType> loc_var_image_it =
         ImageRegionConstIterator<LocalVariationImageType>(
         m_LocalVariation, *fit);
 
       // neighborhood in input image
       ConstShapedNeighborhoodIterator<InputImageType>
         input_image_neighbors_it(size, input, *fit);
       typename ConstShapedNeighborhoodIterator<InputImageType>::
         OffsetType offset;
       input_image_neighbors_it.OverrideBoundaryCondition(&nbc);
       input_image_neighbors_it.ClearActiveList();
       for(int i=0; i<InputImageDimension; i++)
       {
         offset.Fill(0);
         offset[i] = -1;
         input_image_neighbors_it.ActivateOffset(offset);
         offset[i] = 1;
         input_image_neighbors_it.ActivateOffset(offset);
       }
       input_image_neighbors_it.GoToBegin();
 
       // neighborhood in local variation image
       ConstShapedNeighborhoodIterator<LocalVariationImageType>
         loc_var_image_neighbors_it(size, m_LocalVariation, *lv_fit);
       loc_var_image_neighbors_it.OverrideBoundaryCondition(&lv_nbc);
       loc_var_image_neighbors_it.ClearActiveList();
       for(int i=0; i<InputImageDimension; i++)
       {
         offset.Fill(0);
         offset[i] = -1;
         loc_var_image_neighbors_it.ActivateOffset(offset);
         offset[i] = 1;
         loc_var_image_neighbors_it.ActivateOffset(offset);
       }
       loc_var_image_neighbors_it.GoToBegin();
 
       const unsigned int neighborhoodSize = InputImageDimension*2;
       ws.resize(neighborhoodSize);
 
       while ( ! output_image_it.IsAtEnd() )
       {
 
         //MITK_INFO << "looking at voxel " << output_image_it.GetIndex();
 
         //   1 / ||nabla_alpha(u)||_a
         double locvar_alpha_inv = 1.0/loc_var_image_it.Get();
         //MITK_INFO << "locvar: " << loc_var_image_it.Get();
         // compute w_alphabetas
         int count = 0;
         double wsum = 0;
         typename ConstShapedNeighborhoodIterator<LocalVariationImageType>::
           ConstIterator loc_var_neighbors_it;
         for (loc_var_neighbors_it = loc_var_image_neighbors_it.Begin();
           ! loc_var_neighbors_it.IsAtEnd();
           loc_var_neighbors_it++)
         {
           // w_alphabeta(u) =
           //   1 / ||nabla_alpha(u)||_a + 1 / ||nabla_beta(u)||_a
           ws[count] =
             locvar_alpha_inv + (1.0/(double)loc_var_neighbors_it.Get());
           wsum += ws[count++];
 
           //MITK_INFO << "nb var: " << count << ": " << loc_var_neighbors_it.Get();
         }
         //MITK_INFO << "wsum: " << wsum;
 
         // h_alphaalpha * u_alpha^zero
         typename RefImageType::PixelType orig = orig_image_it.Get();
 //        vnl_vector<double> estim(orig.GetSize());
 //        vnl_matrix<double> diff(orig.GetSize(),1);
 //        vnl_matrix<double> estimdash(orig.GetSize(),2);
         vnl_vector_fixed<double,2> step;
         step[0] = 0; step[1]=0;
-        for(int ind=0; ind<m_BValues.size(); ind++)
+        for(size_t ind=0; ind<m_BValues.size(); ind++)
         {
           //MITK_INFO << "refval: " << orig[ind];
           double estim = (1-input_image_it.Get()[0])*exp(-m_BValues[ind]*input_image_it.Get()[1]);
           double estimdash1 = exp(-m_BValues[ind]*input_image_it.Get()[1]);
           double estimdash2 = (-1.0) * (1.0-input_image_it.Get()[0]) * m_BValues[ind] * exp(-m_BValues[ind]*input_image_it.Get()[1]);
           //MITK_INFO << "estimdash1: " << estimdash1 << "; estimdash2: " << estimdash2 << "; diff: " << (double) orig[ind] - estim;
           if(orig[ind] != IVIM_FOO)
           {
             step[0] += ((double) orig[ind] - estim) * estimdash2;
             step[1] += ((double) orig[ind] - estim) * estimdash1;
           }
         }
 
         step[1] *= m_Lambda / (m_Lambda+wsum);
 
         // add the different h_alphabeta * u_beta
         count = 0;
         typename ConstShapedNeighborhoodIterator<InputImageType>::
           ConstIterator input_neighbors_it;
         for (input_neighbors_it = input_image_neighbors_it.Begin();
           ! input_neighbors_it.IsAtEnd();
           input_neighbors_it++)
         {
           step[1] += (input_neighbors_it.Get()[1] - input_image_it.Get()[1]) * (ws[count++] / (m_Lambda+wsum));
         }
 
         //MITK_INFO << "stepfinal: " << step[0] << "; " << step[1];
 
         // set output result
         OutputPixelType out;
         out[0] = input_image_it.Get()[0] + .001*step[0];
         out[1] = input_image_it.Get()[1] + .00001*step[1];
         output_image_it.Set( out );
 
         //MITK_INFO << "(" << input_image_it.Get()[0] << " ; " << input_image_it.Get()[1] << ") => (" << out[0] << " ; " << out[1] << ")";
         // increment iterators
         ++output_image_it;
         ++input_image_it;
         ++orig_image_it;
         ++loc_var_image_it;
         ++input_image_neighbors_it;
         ++loc_var_image_neighbors_it;
 
       }
 
       ++lv_fit;
     }
   }
 
   /**
   * first calculate local variation in the image
   */
   template <class TInputPixel, class TOutputPixel, class TRefPixelType>
   void
   RegularizedIVIMReconstructionSingleIteration<TInputPixel, TOutputPixel, TRefPixelType>
     ::BeforeThreadedGenerateData()
   {
     typedef typename itk::RegularizedIVIMLocalVariationImageFilter
       <InputImageType,LocalVariationImageType> FilterType;
     typename FilterType::Pointer filter = FilterType::New();
     filter->SetInput(this->GetInput(0));
     filter->SetNumberOfThreads(this->GetNumberOfThreads());
     filter->Update();
     this->m_LocalVariation = filter->GetOutput();
   }
 
   /**
   * Standard "PrintSelf" method
   */
   template <class TInputPixel, class TOutputPixel, class TRefPixelType>
   void
   RegularizedIVIMReconstructionSingleIteration<TInputPixel, TOutputPixel, TRefPixelType>
     ::PrintSelf(
     std::ostream& os,
     Indent indent) const
   {
     Superclass::PrintSelf( os, indent );
   }
 
 } // end namespace itk
 
 #endif