diff --git a/Modules/ImageDenoising/Testing/itkTotalVariationDenoisingImageFilterTest.cpp b/Modules/ImageDenoising/Testing/itkTotalVariationDenoisingImageFilterTest.cpp
index 2c8888cc59..2b0622a437 100644
--- a/Modules/ImageDenoising/Testing/itkTotalVariationDenoisingImageFilterTest.cpp
+++ b/Modules/ImageDenoising/Testing/itkTotalVariationDenoisingImageFilterTest.cpp
@@ -1,271 +1,277 @@
 /*============================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center (DKFZ)
 All rights reserved.
 
 Use of this source code is governed by a 3-clause BSD license that can be
 found in the LICENSE file.
 
 ============================================================================*/
 
 #include "itkImageRegionIterator.h"
 #include "itkLocalVariationImageFilter.h"
 #include "itkTotalVariationDenoisingImageFilter.h"
 #include "itkTotalVariationSingleIterationImageFilter.h"
 
 // image typedefs
 typedef itk::Image<float, 3> ImageType;
 typedef itk::ImageRegionIterator<ImageType> IteratorType;
 
 // vector image typedefs
 typedef itk::Vector<float, 2> VectorPixelType;
 typedef itk::Image<VectorPixelType, 3> VectorImageType;
 typedef itk::ImageRegionIterator<VectorImageType> VectorIteratorType;
 
 /**
 * 3x3x3 test image
 */
 ImageType::Pointer GenerateTestImage()
 {
   // init
   ImageType::Pointer image = ImageType::New();
   ;
 
   // spacing
   ImageType::SpacingType spacing;
   spacing[0] = 1;
   spacing[1] = 1;
   spacing[2] = 1;
   image->SetSpacing(spacing);
 
   // extent
   ImageType::RegionType largestPossibleRegion;
   ImageType::SizeType size = {{3, 3, 1}};
   largestPossibleRegion.SetSize(size);
   ImageType::IndexType index = {{0, 0, 0}};
   largestPossibleRegion.SetIndex(index);
   image->SetLargestPossibleRegion(largestPossibleRegion);
   image->SetBufferedRegion(largestPossibleRegion);
 
   // allocate memory
   image->Allocate();
 
   int i = 0;
   IteratorType it(image, largestPossibleRegion);
   it.GoToBegin();
   while (!it.IsAtEnd())
   {
     it.Set((float)i++);
     ++it;
   }
 
   return image;
 }
 
 VectorImageType::Pointer GenerateVectorTestImage()
 {
   // init
   VectorImageType::Pointer image = VectorImageType::New();
   ;
 
   // spacing
   VectorImageType::SpacingType spacing;
   spacing[0] = 1;
   spacing[1] = 1;
   spacing[2] = 1;
   image->SetSpacing(spacing);
 
   // extent
   VectorImageType::RegionType largestPossibleRegion;
   VectorImageType::SizeType size = {{3, 3, 1}};
   largestPossibleRegion.SetSize(size);
   VectorImageType::IndexType index = {{0, 0, 0}};
   largestPossibleRegion.SetIndex(index);
   image->SetLargestPossibleRegion(largestPossibleRegion);
   image->SetBufferedRegion(largestPossibleRegion);
 
   // allocate memory
   image->Allocate();
 
   int i = 0;
   VectorIteratorType it(image, largestPossibleRegion);
   it.GoToBegin();
   while (!it.IsAtEnd())
   {
     VectorPixelType vec;
     vec[0] = (float)i;
     vec[1] = (float)i++;
     it.Set(vec);
     ++it;
   }
 
   return image;
 }
 
 void PrintImage(ImageType::Pointer image)
 {
   IteratorType it(image, image->GetLargestPossibleRegion());
   for (it.GoToBegin(); !it.IsAtEnd(); ++it)
   {
     std::cout << it.Get() << " ";
   }
   std::cout << std::endl;
 }
 
 void PrintVectorImage(VectorImageType::Pointer image)
 {
   VectorIteratorType it(image, image->GetLargestPossibleRegion());
   for (it.GoToBegin(); !it.IsAtEnd(); ++it)
   {
     std::cout << it.Get() << " ";
   }
   std::cout << std::endl;
 }
 
 /**
 * todo
 */
 int itkTotalVariationDenoisingImageFilterTest(int /*argc*/, char * /*argv*/ [])
 {
   ImageType::Pointer image = GenerateTestImage();
   PrintImage(image);
 
   double precision = 0.01;
   ImageType::IndexType index = {{1, 1, 0}};
   VectorImageType::IndexType vecIndex = {{1, 1, 0}};
 
   try
   {
     typedef itk::LocalVariationImageFilter<ImageType, ImageType> LocalFilterType;
     LocalFilterType::Pointer filter = LocalFilterType::New();
     filter->SetInput(image);
     filter->SetNumberOfWorkUnits(1);
     filter->Update();
     ImageType::Pointer outImage = filter->GetOutput();
 
     PrintImage(outImage);
     if (fabs(outImage->GetPixel(index) - 4.472) > precision)
     {
       return EXIT_FAILURE;
     }
   }
-  catch (...)
+  catch (const itk::ExceptionObject& e)
   {
+    e.Print(std::cerr);
     return EXIT_FAILURE;
   }
 
   try
   {
     typedef itk::TotalVariationSingleIterationImageFilter<ImageType, ImageType> SingleFilterType;
     SingleFilterType::Pointer sFilter = SingleFilterType::New();
     sFilter->SetInput(image);
     sFilter->SetOriginalImage(GenerateTestImage());
     sFilter->SetLambda(0.5);
     sFilter->SetNumberOfWorkUnits(1);
     sFilter->Update();
     ImageType::Pointer outImageS = sFilter->GetOutput();
 
     PrintImage(outImageS);
     if (fabs(outImageS->GetPixel(index) - 4.0) > precision)
     {
       return EXIT_FAILURE;
     }
   }
-  catch (...)
+  catch (const itk::ExceptionObject& e)
   {
+    e.Print(std::cerr);
     return EXIT_FAILURE;
   }
 
   try
   {
     typedef itk::TotalVariationDenoisingImageFilter<ImageType, ImageType> TVFilterType;
     TVFilterType::Pointer tvFilter = TVFilterType::New();
     tvFilter->SetInput(image);
     tvFilter->SetNumberIterations(30);
     tvFilter->SetNumberOfWorkUnits(1);
     tvFilter->SetLambda(0.1);
     tvFilter->Update();
     ImageType::Pointer outImageTV = tvFilter->GetOutput();
 
     PrintImage(outImageTV);
     if (fabs(outImageTV->GetPixel(index) - 4.0) > precision)
     {
       return EXIT_FAILURE;
     }
   }
-  catch (...)
+  catch (const itk::ExceptionObject& e)
   {
+    e.Print(std::cerr);
     return EXIT_FAILURE;
   }
 
   VectorImageType::Pointer vecImage = GenerateVectorTestImage();
   PrintVectorImage(vecImage);
 
   try
   {
     typedef itk::LocalVariationImageFilter<VectorImageType, ImageType> LocalVecFilterType;
     LocalVecFilterType::Pointer vecFilter = LocalVecFilterType::New();
     vecFilter->SetInput(vecImage);
     vecFilter->SetNumberOfWorkUnits(1);
     vecFilter->Update();
     ImageType::Pointer outVecImage = vecFilter->GetOutput();
 
     PrintImage(outVecImage);
     if (fabs(outVecImage->GetPixel(index) - 6.324) > precision)
     {
       return EXIT_FAILURE;
     }
   }
-  catch (...)
+  catch (const itk::ExceptionObject& e)
   {
+    e.Print(std::cerr);
     return EXIT_FAILURE;
   }
 
   try
   {
     typedef itk::TotalVariationSingleIterationImageFilter<VectorImageType, VectorImageType> SingleVecFilterType;
     SingleVecFilterType::Pointer sVecFilter = SingleVecFilterType::New();
     sVecFilter->SetInput(vecImage);
     sVecFilter->SetOriginalImage(vecImage);
     sVecFilter->SetLambda(0.5);
     sVecFilter->SetNumberOfWorkUnits(1);
     sVecFilter->UpdateLargestPossibleRegion();
     VectorImageType::Pointer outVecImageS = sVecFilter->GetOutput();
 
     PrintVectorImage(outVecImageS);
     if (fabs(outVecImageS->GetPixel(vecIndex)[1] - 4.0) > precision)
     {
       return EXIT_FAILURE;
     }
   }
-  catch (...)
+  catch (const itk::ExceptionObject& e)
   {
+    e.Print(std::cerr);
     return EXIT_FAILURE;
   }
 
   try
   {
     typedef itk::TotalVariationDenoisingImageFilter<VectorImageType, VectorImageType> TVVectorFilterType;
     TVVectorFilterType::Pointer tvVecFilter = TVVectorFilterType::New();
     tvVecFilter->SetInput(vecImage);
     tvVecFilter->SetNumberIterations(30);
     tvVecFilter->SetNumberOfWorkUnits(1);
     tvVecFilter->SetLambda(0.1);
     tvVecFilter->Update();
     VectorImageType::Pointer outVecImageTV = tvVecFilter->GetOutput();
 
     PrintVectorImage(outVecImageTV);
     if (fabs(outVecImageTV->GetPixel(vecIndex)[1] - 4.0) > precision)
     {
       return EXIT_FAILURE;
     }
   }
-  catch (...)
+  catch (const itk::ExceptionObject& e)
   {
+    e.Print(std::cerr);
     return EXIT_FAILURE;
   }
 
   return EXIT_SUCCESS;
 }
diff --git a/Modules/ImageDenoising/itkLocalVariationImageFilter.txx b/Modules/ImageDenoising/itkLocalVariationImageFilter.txx
index 679247d56b..be84c917bb 100644
--- a/Modules/ImageDenoising/itkLocalVariationImageFilter.txx
+++ b/Modules/ImageDenoising/itkLocalVariationImageFilter.txx
@@ -1,188 +1,189 @@
 /*============================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center (DKFZ)
 All rights reserved.
 
 Use of this source code is governed by a 3-clause BSD license that can be
 found in the LICENSE file.
 
 ============================================================================*/
 
 /*===================================================================
 
 This file is based heavily on a corresponding ITK filter.
 
 ===================================================================*/
 #ifndef _itkLocalVariationImageFilter_txx
 #define _itkLocalVariationImageFilter_txx
 #include "itkLocalVariationImageFilter.h"
 
 #include "itkConstShapedNeighborhoodIterator.h"
 #include "itkImageRegionIterator.h"
 #include "itkNeighborhoodAlgorithm.h"
 #include "itkNeighborhoodInnerProduct.h"
 #include "itkOffset.h"
 #include "itkProgressReporter.h"
 #include "itkVectorImage.h"
 #include "itkZeroFluxNeumannBoundaryCondition.h"
 
 #include <algorithm>
 #include <vector>
 
 namespace itk
 {
   template <class TInputImage, class TOutputImage>
   LocalVariationImageFilter<TInputImage, TOutputImage>::LocalVariationImageFilter()
   {
+    this->DynamicMultiThreadingOff();
   }
 
   template <class TInputImage, class TOutputImage>
   void LocalVariationImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRegion()
   {
     // 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 <>
   double SquaredEuclideanMetric<itk::VariableLengthVector<float>>::Calc(itk::VariableLengthVector<float> p)
   {
     return p.GetSquaredNorm();
   }
 
   template <>
   double SquaredEuclideanMetric<itk::VariableLengthVector<double>>::Calc(itk::VariableLengthVector<double> p)
   {
     return p.GetSquaredNorm();
   }
 
   template <class TPixelType>
   double SquaredEuclideanMetric<TPixelType>::Calc(TPixelType p)
   {
     return p * p;
   }
 
   template <class TInputImage, class TOutputImage>
   void LocalVariationImageFilter<TInputImage, TOutputImage>::ThreadedGenerateData(
     const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId)
   {
     // Allocate output
     typename OutputImageType::Pointer output = this->GetOutput();
     typename InputImageType::ConstPointer input = this->GetInput();
 
     itk::Size<InputImageDimension> size;
     for (unsigned 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 (auto 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 (unsigned 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 += SquaredEuclideanMetric<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 LocalVariationImageFilter<TInputImage, TOutput>::PrintSelf(std::ostream &os, Indent indent) const
   {
     Superclass::PrintSelf(os, indent);
   }
 
 } // end namespace itk
 
 #endif //_itkLocalVariationImageFilter_txx
diff --git a/Modules/ImageDenoising/itkTotalVariationSingleIterationImageFilter.txx b/Modules/ImageDenoising/itkTotalVariationSingleIterationImageFilter.txx
index b409b2500f..a069c1d957 100644
--- a/Modules/ImageDenoising/itkTotalVariationSingleIterationImageFilter.txx
+++ b/Modules/ImageDenoising/itkTotalVariationSingleIterationImageFilter.txx
@@ -1,251 +1,252 @@
 /*============================================================================
 
 The Medical Imaging Interaction Toolkit (MITK)
 
 Copyright (c) German Cancer Research Center (DKFZ)
 All rights reserved.
 
 Use of this source code is governed by a 3-clause BSD license that can be
 found in the LICENSE file.
 
 ============================================================================*/
 
 /*===================================================================
 
 This file is based heavily on a corresponding ITK filter.
 
 ===================================================================*/
 
 #ifndef _itkTotalVariationSingleIterationImageFilter_txx
 #define _itkTotalVariationSingleIterationImageFilter_txx
 
 #include "itkTotalVariationSingleIterationImageFilter.h"
 
 // itk includes
 #include "itkConstShapedNeighborhoodIterator.h"
 #include "itkImageRegionIterator.h"
 #include "itkLocalVariationImageFilter.h"
 #include "itkNeighborhoodAlgorithm.h"
 #include "itkNeighborhoodInnerProduct.h"
 #include "itkOffset.h"
 #include "itkProgressReporter.h"
 #include "itkZeroFluxNeumannBoundaryCondition.h"
 
 // other includes
 #include <algorithm>
 #include <vector>
 
 namespace itk
 {
   /**
   * constructor
   */
   template <class TInputImage, class TOutputImage>
   TotalVariationSingleIterationImageFilter<TInputImage, TOutputImage>::TotalVariationSingleIterationImageFilter()
   {
+    this->DynamicMultiThreadingOff();
     m_Lambda = 1.0;
     m_LocalVariation = LocalVariationImageType::New();
   }
 
   /**
  * generate requested region
  */
   template <class TInputImage, class TOutputImage>
   void TotalVariationSingleIterationImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRegion()
   {
     // 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;
     }
   }
 
   /**
  * generate output
  */
   template <class TInputImage, class TOutputImage>
   void TotalVariationSingleIterationImageFilter<TInputImage, TOutputImage>::ThreadedGenerateData(
     const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId)
   {
     typename OutputImageType::Pointer output = this->GetOutput();
     typename InputImageType::ConstPointer input = this->GetInput();
 
     // Find the data-set boundary "faces"
     itk::Size<InputImageDimension> size;
     for (unsigned 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);
 
     // support progress methods/callbacks
     ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
 
     ZeroFluxNeumannBoundaryCondition<InputImageType> nbc;
     ZeroFluxNeumannBoundaryCondition<LocalVariationImageType> lv_nbc;
     std::vector<double> ws;
     std::vector<double> hs;
 
     auto lv_fit = lv_faceList.begin();
 
     // Process each of the boundary faces.  These are N-d regions which border
     // the edge of the buffer.
     for (auto 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<InputImageType> orig_image_it =
         ImageRegionConstIterator<InputImageType>(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 (unsigned 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 (unsigned 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())
       {
         //   1 / ||nabla_alpha(u)||_a
         double locvar_alpha_inv = 1.0 / 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++];
         }
 
         // h_alphaalpha * u_alpha^zero
         typename OutputImageType::PixelType res = static_cast<typename OutputImageType::PixelType>(
           ((typename OutputImageType::PixelType)orig_image_it.Get()) * (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++)
         {
           res += input_neighbors_it.Get() * (ws[count++] / (m_Lambda + wsum));
         }
 
         // set output result
         output_image_it.Set(res);
 
         // increment iterators
         ++output_image_it;
         ++input_image_it;
         ++orig_image_it;
         ++loc_var_image_it;
         ++input_image_neighbors_it;
         ++loc_var_image_neighbors_it;
 
         // report progress
         progress.CompletedPixel();
       }
 
       ++lv_fit;
     }
   }
 
   /**
   * first calculate local variation in the image
   */
   template <class TInputImage, class TOutputImage>
   void TotalVariationSingleIterationImageFilter<TInputImage, TOutputImage>::BeforeThreadedGenerateData()
   {
     typedef typename itk::LocalVariationImageFilter<TInputImage, LocalVariationImageType> FilterType;
     typename FilterType::Pointer filter = FilterType::New();
     filter->SetInput(this->GetInput(0));
     filter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
     filter->Update();
     this->m_LocalVariation = filter->GetOutput();
   }
 
   /**
   * Standard "PrintSelf" method
   */
   template <class TInputImage, class TOutput>
   void TotalVariationSingleIterationImageFilter<TInputImage, TOutput>::PrintSelf(std::ostream &os, Indent indent) const
   {
     Superclass::PrintSelf(os, indent);
   }
 
 } // end namespace itk
 
 #endif