/*
 * Copyright (c) ICG. All rights reserved.
 * See copyright.txt for more information.
 *
 * Institute for Computer Graphics and Vision
 * Graz, University of Technology / Austria
 *
 *
 * This software is distributed WITHOUT ANY WARRANTY; without even
 * the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 * PURPOSE.  See the above copyright notices for more information.
 *
 *
 * Project     : MIPItkProjects
 * Module      : CommonItkUtilities
 * Class       : $RCSfile: SmallUtilityMethods.txx,v $
 * Language    : C++
 * Description : 
 *
 * Author     : Martin Urschler
 * EMail      : urschler@icg.tu-graz.ac.at
 * Date       : $Date: 2007-04-23 11:43:48 $
 * Version    : $Revision: 1.35 $
 * Full Id    : $Id: SmallUtilityMethods.txx,v 1.35 2007-04-23 11:43:48 urschler Exp $
 *
 */


#ifndef __SMALL_UTILITY_METHODS_TXX__
#define __SMALL_UTILITY_METHODS_TXX__

 
#include "itkShrinkImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkAbsoluteValueDifferenceImageFilter.h"
#include "itkSquaredDifferenceImageFilter.h"
#include "itkSubtractImageFilter.h"
#include "itkMinimumMaximumImageFilter.h"
#include "itkCurvatureFlowImageFilter.h"
#include "itkConnectedComponentImageFilter.h"
#include "itkRelabelComponentImageFilter.h"
#include "itkMultiplyImageFilter.h"
#include "itkShiftScaleImageFilter.h"
#include "itkConnectedThresholdImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkRegionOfInterestImageFilter.h"
#include "itkImageRegionConstIteratorWithIndex.h"
#include "itkRecursiveGaussianImageFilter.h"
#include "itkDiscreteGaussianImageFilter.h"
#include "itkIdentityTransform.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkWindowedSincInterpolateImageFunction.h"
#include "itkResampleImageFilter.h"





template< typename TImageType >
typename SmallUtilityMethods<TImageType>::TImageTypePointer
SmallUtilityMethods<TImageType>::shrinkVolume(
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image, 
   unsigned int shrink_factor, bool performSmoothing )
{
   return SmallUtilityMethods<TImageType>::shrinkVolume(
      image, shrink_factor, shrink_factor, shrink_factor, performSmoothing);
}


template< typename TImageType >
typename SmallUtilityMethods<TImageType>::TImageTypePointer 
SmallUtilityMethods<TImageType>::shrinkVolume( 
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image,
   unsigned int shrink_factor_x, 
   unsigned int shrink_factor_y, 
   unsigned int shrink_factor_z,
   bool performSmoothing )
{
   if (shrink_factor_x <= 1 && shrink_factor_y <= 1 && shrink_factor_z <= 1)
   {
      // do nothing, shrink factors are erroneous (0) or equal to 1, meaning
      // that we don't want to shrink
      
      return image;
   }
   
   typedef itk::ShrinkImageFilter< TImageType, TImageType > ShrinkFilterType;
   typename ShrinkFilterType::Pointer shrinker = ShrinkFilterType::New();

   std::string message;

   if (performSmoothing)
   {
      // sigmas are calculated in real world units!
      const double sigma_x = 
         0.5 * shrink_factor_x * image->GetSpacing()[0];
      const double sigma_y = 
         0.5 * shrink_factor_y * image->GetSpacing()[1];
      const double sigma_z = 
         0.5 * shrink_factor_z * image->GetSpacing()[2];

      TImageTypePointer smoothed_image = 
         gaussSmoothVolume( image, sigma_x, sigma_y, sigma_z );

      shrinker->SetInput( smoothed_image );

      message = std::string("Apply shrinker WITH gaussian smoothing.");
   }
   else
   {
      shrinker->SetInput( image );

      message = std::string("Apply shrinker WITHOUT smoothing.");
   }

   shrinker->SetShrinkFactor( 0, shrink_factor_x );
   shrinker->SetShrinkFactor( 1, shrink_factor_y );
   shrinker->SetShrinkFactor( 2, shrink_factor_z );
   
   ITK_EXCEPTION_CHECKED( message.c_str(), shrinker->Update(), 0 );

   return shrinker->GetOutput();
}


template< typename TImageType >
typename SmallUtilityMethods<TImageType>::TImageTypePointer 
SmallUtilityMethods<TImageType>::resampleVolumeWithoutSmoothing( 
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image,
   const typename SmallUtilityMethods<TImageType>::TImageSizeType& new_size )
{
   typename TImageType::SizeType original_size = 
      image->GetLargestPossibleRegion().GetSize();

   const unsigned int Dimensions = 
      itk::GetImageDimension<TImageType>::ImageDimension;

   double scale_factor[Dimensions];
   scale_factor[0] = 
      static_cast<double>(original_size[0]) / static_cast<double>(new_size[0]);
   scale_factor[1] = 
      static_cast<double>(original_size[1]) / static_cast<double>(new_size[1]);
   scale_factor[2] = 
      static_cast<double>(original_size[2]) / static_cast<double>(new_size[2]);

   typedef itk::ResampleImageFilter< TImageType, TImageType >  FilterType;
   typename FilterType::Pointer filter = FilterType::New();

   typedef itk::IdentityTransform< double, Dimensions > TransformType;
   typename TransformType::Pointer transform = TransformType::New();

   filter->SetDefaultPixelValue( 
      SmallUtilityMethods<TImageType>::getMinimumVoxelValue(image) );

   typename TImageType::SpacingType spacing = image->GetSpacing();
   spacing[0] *= scale_factor[0];
   spacing[1] *= scale_factor[1];
   spacing[2] *= scale_factor[2];

   std::cout << "Resample Scale Factors: ";
   std::cout << scale_factor[0] << " " << 
                scale_factor[1] << " " << 
                scale_factor[2] << std::endl;

   filter->SetOutputSpacing( spacing );
   filter->SetOutputOrigin( image->GetOrigin() );
   filter->SetSize( new_size );
   filter->SetOutputStartIndex( image->GetLargestPossibleRegion().GetIndex() );

   filter->SetInput( image );
   filter->SetTransform( transform );

   const unsigned int WindowedSincRadius = 4;
   typedef itk::WindowedSincInterpolateImageFunction< 
      TImageType,  
      WindowedSincRadius,
      itk::Function::LanczosWindowFunction<WindowedSincRadius>,
      itk::ConstantBoundaryCondition<TImageType>,
      double > 
      InterpolatorType;
   typename InterpolatorType::Pointer interpolator = InterpolatorType::New();

   filter->SetInterpolator( interpolator );

   try
   {
      std::cout << "Resampling input image!" << std::endl;
      filter->Update();
   }
   catch ( itk::MemoryAllocationError& err )
   {
      std::cout << "ITK Memory Allocation Error caught!" << std::endl; 
      std::cout << err << std::endl; 
      std::cout << "We don't have enough memory to perform windowed sinc " <<
                   "interpolation, let's try a linear interpolation!"
                << std::endl;
      
      typedef itk::LinearInterpolateImageFunction< TImageType, double >  
         LinearInterpolatorType;
      typename LinearInterpolatorType::Pointer interpolator = 
         LinearInterpolatorType::New();

      filter->SetInterpolator( interpolator );
      
      ITK_EXCEPTION_CHECKED( "Resampling input image!", filter->Update(), 0 );
   }

   return filter->GetOutput();
}





// standard recursive gaussian image filter

template< typename TImageType >
typename SmallUtilityMethods<TImageType>::TImageTypePointer 
SmallUtilityMethods<TImageType>::gaussSmoothVolume( 
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image, 
   double sigma_x, 
   double sigma_y, 
   double sigma_z )
{
   typename SmallUtilityMethods<TImageType>::TImageTypePointer 
      result = createNewImage( image );
   typename SmallUtilityMethods<TImageType>::TImageTypePointer 
      temporary_result = createNewImage( image );

   typedef itk::RecursiveGaussianImageFilter<TImageType, TImageType> 
      SmoothingFilterType;

   std::cout << "Recursive gauss-smoothing sigmas: " << sigma_x << " ; " << 
      sigma_y << " ; " << sigma_z << std::endl;

   // implement the x convolution given input 'image' and store the filter 
   // output in 'result'
   {
      typename SmoothingFilterType::Pointer filter = 
         SmoothingFilterType::New();

      filter->SetDirection( 0 ); // 0 --> X direction
      filter->SetSigma( sigma_x );
      filter->SetOrder( SmoothingFilterType::ZeroOrder );
      filter->SetNormalizeAcrossScale( false );
      filter->SetInput( image );
      filter->SetReleaseDataFlag(true);
      filter->GraftOutput( result );
      //filter->Print(std::cout);

      ITK_EXCEPTION_CHECKED( "Smoothing image in direction X!", 
         filter->Update(), 0 );
   }

   // implement the y convolution given input 'result' and store the filter 
   // output in 'temporary_result'
   {
      typename SmoothingFilterType::Pointer filter = 
         SmoothingFilterType::New();

      filter->SetDirection( 1 ); // 1 --> Y direction
      filter->SetSigma( sigma_y );
      filter->SetOrder( SmoothingFilterType::ZeroOrder );
      filter->SetNormalizeAcrossScale( false );
      filter->SetInput( result );
      filter->SetReleaseDataFlag(true);
      filter->GraftOutput( temporary_result );
      //filter->Print(std::cout);

      ITK_EXCEPTION_CHECKED( "Smoothing image in direction Y!", 
         filter->Update(), 0 );
   }

   // implement the z convolution given input 'temporary_result' and store the 
   // filter output in 'result'
   {
      typename SmoothingFilterType::Pointer filter = 
         SmoothingFilterType::New();

      filter->SetDirection( 2 ); // 2 --> Z direction
      filter->SetSigma( sigma_z );
      filter->SetOrder( SmoothingFilterType::ZeroOrder );
      filter->SetNormalizeAcrossScale( false );
      filter->SetReleaseDataFlag(true);
      filter->SetInput( temporary_result );
      filter->GraftOutput( result );
      //filter->Print(std::cout);

      ITK_EXCEPTION_CHECKED( "Smoothing image in direction Z!", 
         filter->Update(), 0 );
   }
      
   return result;
}

template< typename TImageType >
typename SmallUtilityMethods<TImageType>::TImageTypePointer 
SmallUtilityMethods<TImageType>::gaussSmoothDerivativeVolume( 
  const typename SmallUtilityMethods<TImageType>::TImageTypeConstPointer& image, 
  double sigma,
  unsigned int direction ) 
{
   if (direction > 2)
   {
      std::cout << "Direction larger than 2, doesn't make sense!" << std::endl;
      return 0;
   }
   
   typedef itk::RecursiveGaussianImageFilter<TImageType, TImageType> 
      DerivativeFilterType;

   typename DerivativeFilterType::Pointer filter = DerivativeFilterType::New();
   filter->SetDirection( direction );
   filter->SetOrder( DerivativeFilterType::FirstOrder );
   filter->SetSigma( sigma );
   filter->SetNormalizeAcrossScale( false );
   filter->SetInput( image );
   ITK_EXCEPTION_CHECKED( "Derivative of image!",
         filter->Update(), 0 );

   return filter->GetOutput();
}

template< typename TImageType >
typename SmallUtilityMethods<TImageType>::TImageTypePointer 
SmallUtilityMethods<TImageType>::gaussSmoothDerivativeVolume( 
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image, 
   double sigma,
   unsigned int direction ) 
{
   if (direction > 2)
   {
      std::cout << "Direction larger than 2, doesn't make sense!" << std::endl;
      return 0;
   }
   
   typedef itk::RecursiveGaussianImageFilter<TImageType, TImageType> 
      DerivativeFilterType;

   typename DerivativeFilterType::Pointer filter = DerivativeFilterType::New();
   filter->SetDirection( direction );
   filter->SetOrder( DerivativeFilterType::FirstOrder );
   filter->SetSigma( sigma );
   filter->SetNormalizeAcrossScale( false );
   filter->SetInput( image );
   ITK_EXCEPTION_CHECKED( "Derivative of image!",
         filter->Update(), 0 );

   return filter->GetOutput();
}

template< typename TImageType >
typename SmallUtilityMethods<TImageType>::TImageTypePointer 
SmallUtilityMethods<TImageType>::discreteGaussSmoothVolume( 
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image, 
   double sigma_x, 
   double sigma_y, 
   double sigma_z )
{
   typedef itk::DiscreteGaussianImageFilter<TImageType, TImageType> 
      SmoothingFilterType;

   typename SmoothingFilterType::Pointer smoother = 
      SmoothingFilterType::New();
   
   const unsigned int Dimensions = 
      itk::GetImageDimension<TImageType>::ImageDimension;
      
   smoother->SetInput( image );
   smoother->SetUseImageSpacingOn();
   
   std::cout << "Discrete gauss-smoothing sigmas in mm: " << sigma_x << " ; " << 
      sigma_y << " ; " << sigma_z << std::endl;

   double variances[Dimensions];
   variances[0] = vnl_math_sqr( sigma_x ); 
   variances[1] = vnl_math_sqr( sigma_y ); 
   variances[2] = vnl_math_sqr( sigma_z );
   smoother->SetVariance( variances );

   //std::cout << "variance x = " << variances[0] << " mm^2" << std::endl;
   //std::cout << "variance y = " << variances[1] << " mm^2" << std::endl;
   //std::cout << "variance z = " << variances[2] << " mm^2" << std::endl;

   //smoother->Print(std::cout); 

   ITK_EXCEPTION_CHECKED( "Smoothing image with discrete gauss filter!", 
      smoother->Update(), 0 );

   return smoother->GetOutput();
}





// some old code

//   double shrink_factor_xy = pow( 2, shrink_level );
//
//   double shrink_factor_z;
//   if (application_type == CTDemonsRegistration::LUNG_MODE) 
//      shrink_factor_z = shrink_factor_xy;
//   else if (application_type == CTDemonsRegistration::LIVER_MODE)
//      shrink_factor_z = 1;
//
//   std::cout << "variance: " << vnl_math_sqr( 0.5 * shrink_factor_xy ) 
//      << std::endl;
//
//   SmoothingFilterType::Pointer fixed_smoother = SmoothingFilterType::New();
//   fixed_smoother->SetUseImageSpacingOff();
//   fixed_smoother->SetVariance( vnl_math_sqr( 0.5 * shrink_factor_xy ) );
//
//   SmoothingFilterType::Pointer moving_smoother = SmoothingFilterType::New();
//   moving_smoother->SetUseImageSpacingOff();
//   moving_smoother->SetVariance( vnl_math_sqr( 0.5 * shrink_factor_xy ) );



template< typename TImageType >
typename SmallUtilityMethods<TImageType>::TImageTypePointer 
SmallUtilityMethods<TImageType>::curvFlowSmoothVolume( 
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image,
   unsigned int nr_iterations, 
   double time_step )
{
   // cast input type to float type

   typedef itk::CastImageFilter< TImageType, FloatImageType > CasterType;
   typename CasterType::Pointer caster = CasterType::New();
   caster->SetInput( image );
   
   // setup the Image Smoother

   typedef itk::CurvatureFlowImageFilter< FloatImageType, FloatImageType >
      CurvatureFlowImageFilterType;

   CurvatureFlowImageFilterType::Pointer smoother = 
      CurvatureFlowImageFilterType::New();

   smoother->SetNumberOfIterations( nr_iterations );
   smoother->SetTimeStep( time_step );

   smoother->SetInput( caster->GetOutput() );

   // cast back to input type
   
   typedef itk::CastImageFilter< FloatImageType, TImageType > BackCasterType;
   typename BackCasterType::Pointer back_caster = BackCasterType::New();
   back_caster->SetInput( smoother->GetOutput() );
   
   ITK_EXCEPTION_CHECKED( "Image smoothing using curvature flow filter (input image was first of all cast to a float image!).",
      back_caster->Update(), 0 );

   return back_caster->GetOutput();
}


template< typename TImageType >
typename SmallUtilityMethods<TImageType>::TImageTypePointer 
SmallUtilityMethods<TImageType>::removeBackground( 
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image,
   MaskImageType::Pointer& mask_image,
   const double replaceValue,
   const double lowerThreshold,
   const double upperThreshold )
{
   // setup the Region Grower for Background Removal

   typename TImageType::SizeType image_size = 
      image->GetLargestPossibleRegion().GetSize();

   typedef itk::ConnectedThresholdImageFilter< TImageType, TImageType >
        RegionGrowFilterType;

   typename RegionGrowFilterType::Pointer regiongrower = 
      RegionGrowFilterType::New();

   typename TImageType::IndexType seed_index1;
   seed_index1[0] = 0;
   seed_index1[1] = 0;
   seed_index1[2] = 0;
   regiongrower->SetSeed( seed_index1 );

   regiongrower->SetReplaceValue( 
      static_cast<typename TImageType::PixelType>(replaceValue) );
   regiongrower->SetLower( 
      static_cast<typename TImageType::PixelType>(lowerThreshold) );
   regiongrower->SetUpper( 
      static_cast<typename TImageType::PixelType>(upperThreshold) );

   //std::cout << image_size[0] << std::endl;
   //std::cout << image_size[1] << std::endl;
   //std::cout << image_size[2] << std::endl;
   
   typename TImageType::IndexType seed_index2;
   seed_index2[0] = image_size[0]-1;
   seed_index2[1] = image_size[1]-1;
   seed_index2[2] = 0;

   regiongrower->AddSeed( seed_index2 );

   typename TImageType::IndexType seed_index3;
   seed_index3[0] = image_size[0]/2;
   seed_index3[1] = image_size[1]-1;
   seed_index3[2] = 0;

   regiongrower->AddSeed( seed_index3 );

   typename TImageType::IndexType seed_index4;
   seed_index4[0] = 0;
   seed_index4[1] = image_size[1]-1;
   seed_index4[2] = 0;

   regiongrower->AddSeed( seed_index4 );

   typename TImageType::IndexType seed_index5;
   seed_index5[0] = image_size[0]-1;
   seed_index5[1] = 0;
   seed_index5[2] = 0;

   regiongrower->AddSeed( seed_index5 );   


   // get the minimum from the smoother output image
   typename TImageType::PixelType input_image_minimum =
      SmallUtilityMethods<TImageType>::getMinimumVoxelValue( image );
   


   // setup the shift scaler

   typedef itk::ShiftScaleImageFilter< TImageType, TImageType > 
      ShiftScaleImageFilterType;

   typename ShiftScaleImageFilterType::Pointer positive_shifter = 
      ShiftScaleImageFilterType::New();
   typename ShiftScaleImageFilterType::Pointer negative_shifter = 
      ShiftScaleImageFilterType::New();
     
   positive_shifter->SetScale(1.0);
   negative_shifter->SetScale(1.0);
   
   if (input_image_minimum < 0)
   {
      positive_shifter->SetShift(-input_image_minimum);
      negative_shifter->SetShift( input_image_minimum);
   }
   else
   {
      positive_shifter->SetShift(0.0);
      negative_shifter->SetShift(0.0);
   }


   //typedef unsigned short InternalPixelMaskType;
   //typedef itk::Image< InternalPixelMaskType, _DIMENSION_ > 
   //   InternalPixelMaskImageType;

   //typedef unsigned char InternalPixelMask8BitType;
   //typedef itk::Image< InternalPixelMask8BitType, _DIMENSION_ > 
   //   InternalPixelMask8BitImageType;

   typedef itk::RescaleIntensityImageFilter< TImageType, UInt16ImageType >
      PixelMaskCastImageFilterType;

   typename PixelMaskCastImageFilterType::Pointer mask_caster = 
      PixelMaskCastImageFilterType::New();
   mask_caster->SetOutputMinimum( 0 );
   mask_caster->SetOutputMaximum( 1 );


   typedef itk::ConnectedComponentImageFilter<UInt16ImageType, UInt16ImageType>
      PixelMaskComponentLabeler;
   PixelMaskComponentLabeler::Pointer mask_labeler = 
      PixelMaskComponentLabeler::New();

   typedef itk::RelabelComponentImageFilter< UInt16ImageType, UInt16ImageType >
      PixelMaskRelabeler;
   PixelMaskRelabeler::Pointer mask_relabeler = PixelMaskRelabeler::New();
   mask_relabeler->InPlaceOn();

   typedef itk::BinaryThresholdImageFilter< UInt16ImageType, MaskImageType >
      RelabeledDataThresholderType;

   RelabeledDataThresholderType::Pointer mask_relabeled_thresholder = 
      RelabeledDataThresholderType::New();
   
   // only take the largest object!
   mask_relabeled_thresholder->SetLowerThreshold( 1 ); 
   mask_relabeled_thresholder->SetUpperThreshold( 1 );
   mask_relabeled_thresholder->SetInsideValue( 1 );
   mask_relabeled_thresholder->SetOutsideValue( 0 );

   typedef itk::MultiplyImageFilter< TImageType, MaskImageType, TImageType >
      MultiplyImageFilterType;

   typename MultiplyImageFilterType::Pointer multiplier = 
      MultiplyImageFilterType::New();



   // construct the pipeline
      
   regiongrower->SetInput( image );
   positive_shifter->SetInput( image );

   mask_caster->SetInput( regiongrower->GetOutput() );
   mask_labeler->SetInput( mask_caster->GetOutput() );
   mask_relabeler->SetInput( mask_labeler->GetOutput() );
   mask_relabeled_thresholder->SetInput( mask_relabeler->GetOutput() );

   multiplier->SetInput1( positive_shifter->GetOutput() );
   multiplier->SetInput2( mask_relabeled_thresholder->GetOutput() );

   negative_shifter->SetInput( multiplier->GetOutput() );

   ITK_EXCEPTION_CHECKED( "Apply whole bunch of preprocessing steps to remove background from CT image.",
      negative_shifter->Update(), 0 );

   // assign mask image to the output parameter
   mask_image = mask_relabeled_thresholder->GetOutput();
   
   return negative_shifter->GetOutput();
}





template< typename TImageType >
unsigned long
SmallUtilityMethods<TImageType>::countNonZeroVoxels( 
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image )
{
   typedef itk::ImageRegionConstIterator<TImageType> IteratorType;

   unsigned long count = 0;

   IteratorType it = IteratorType ( image, image->GetRequestedRegion() );
   it.GoToBegin();
   while( !it.IsAtEnd())
   {
     if ( it.Value() != 0 )
        ++count;

     ++it;
   }

   return count;
}



template< typename TImageType >
typename SmallUtilityMethods<TImageType>::TImagePixelType 
SmallUtilityMethods<TImageType>::getMinimumVoxelValue( 
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image )
{
   typename itk::MinimumMaximumImageCalculator< TImageType >::Pointer 
      min_calculator = itk::MinimumMaximumImageCalculator< TImageType >::New();

   min_calculator->SetImage( image );
   min_calculator->ComputeMinimum();
   return min_calculator->GetMinimum();
}

template< typename TImageType >
typename SmallUtilityMethods<TImageType>::TImagePixelType 
SmallUtilityMethods<TImageType>::getMaximumVoxelValue( 
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image )
{
   typename itk::MinimumMaximumImageCalculator< TImageType >::Pointer 
      min_calculator = itk::MinimumMaximumImageCalculator< TImageType >::New();

   min_calculator->SetImage( image );
   min_calculator->ComputeMaximum();
   return min_calculator->GetMaximum();
}


template< typename TImageType >
double
SmallUtilityMethods<TImageType>::getAverageGreyValueAtSamplePoints(
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image, 
	const std::vector< typename SmallUtilityMethods<TImageType>::
      TImageIndexType >& seeds
   )
{
   double avg_val = 0.0;
   
   for (int i=0; i<seeds.size(); ++i)
   {
      avg_val += image->GetPixel( seeds[i] );
   }

   avg_val /= static_cast<double>( seeds.size() );

   return avg_val;
}



template< typename TImageType >
typename SmallUtilityMethods<TImageType>::TImageTypePointer 
SmallUtilityMethods<TImageType>::calculateDifferenceImage( 
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image1,
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image2 )
{
   typedef itk::SubtractImageFilter< 
                                  TImageType, 
                                  TImageType, 
                                  TImageType > DifferenceFilterType;

   typename DifferenceFilterType::Pointer difference = 
      DifferenceFilterType::New();
   difference->SetInput1( image1 );
   difference->SetInput2( image2 );

   ITK_EXCEPTION_CHECKED( "Apply difference image filter.", 
      difference->Update(), 0 );

   return difference->GetOutput();
}



template< typename TImageType >
typename SmallUtilityMethods<TImageType>::TImageTypePointer
SmallUtilityMethods<TImageType>::createNewImage( 
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& 
      prototypeImage )
{
   // create result image
   typename SmallUtilityMethods<TImageType>::TImageTypePointer result = 
      TImageType::New();

   result->SetLargestPossibleRegion( 
      prototypeImage->GetLargestPossibleRegion() );
   result->SetRequestedRegion( prototypeImage->GetRequestedRegion() );
   result->SetBufferedRegion( prototypeImage->GetBufferedRegion() );
   
   result->SetSpacing( prototypeImage->GetSpacing() );
   result->SetOrigin( prototypeImage->GetOrigin() ); 
   result->SetDirection( prototypeImage->GetDirection() );
   
   result->Allocate();

   return result;
}

template< typename TImageType >
DeformationFieldType::Pointer
SmallUtilityMethods<TImageType>::createNewZeroedDeformationField( 
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& 
      prototypeImage )
{ 
   // create result image
   DeformationFieldType::Pointer result = DeformationFieldType::New();

   result->SetLargestPossibleRegion( 
      prototypeImage->GetLargestPossibleRegion() );
   result->SetRequestedRegion( prototypeImage->GetRequestedRegion() );
   result->SetBufferedRegion( prototypeImage->GetBufferedRegion() );
   
   result->SetSpacing( prototypeImage->GetSpacing() );
   result->SetOrigin( prototypeImage->GetOrigin() ); 
   result->SetDirection( prototypeImage->GetDirection() );
   
   result->Allocate();
   
   VectorPixelType zeroDeformation = 
      VectorPixelConstants<VectorPixelType::ValueType, 3>::zero();
   result->FillBuffer( zeroDeformation );

   return result;
}

template< typename TImageType >
typename SmallUtilityMethods<TImageType>::TImageTypePointer 
SmallUtilityMethods<TImageType>::calculateCheckerboardImage( 
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image1,
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image2 )
{
   typename TImageType::SizeType image1_size  = 
      image1->GetLargestPossibleRegion().GetSize();
   typename TImageType::SizeType image2_size  = 
      image2->GetLargestPossibleRegion().GetSize();

   if (image1_size != image2_size)
   {
      std::cout << "Error image sizes not equal!" << std::endl;
      return 0;
   }

   typename SmallUtilityMethods<TImageType>::TImageTypePointer result =
      createNewImage( image1 );

   typedef itk::ImageRegionConstIterator<TImageType> ConstIteratorType;
   typedef itk::ImageRegionIterator<TImageType> IteratorType;


   ConstIteratorType in_it1 = ConstIteratorType ( image1, 
      image1->GetLargestPossibleRegion() );
   ConstIteratorType in_it2 = ConstIteratorType ( image2, 
      image2->GetLargestPossibleRegion() );
   IteratorType out_it = IteratorType ( result, 
      result->GetLargestPossibleRegion() );

   in_it1.GoToBegin();
   in_it2.GoToBegin();
   out_it.GoToBegin();

   unsigned long count = 0;
   unsigned long slice = 0;
   unsigned long row = 0;
   unsigned long col = 0;

   unsigned long slice_size = image1_size[0] * image1_size[1];

   while( !out_it.IsAtEnd())
   {
	  slice = (count / slice_size) / 4;
	  unsigned long in_slice_count = count % slice_size;
	  row = (in_slice_count / image1_size[0]) / 16;
	  col = (in_slice_count % image1_size[0]) / 16;

     //std::cout << "(slice,row,col) = (" << slice << "," << row << "," << 
     //   col << ")" << std::endl;

     if (slice % 2)
     {
        if ( (row+col) % 2 )
           out_it.Set( in_it1.Get() );
        else
           out_it.Set( in_it2.Get() );
     }
     else
     {
        if ( (row+col) % 2 )
           out_it.Set( in_it2.Get() );
        else
           out_it.Set( in_it1.Get() );
     }

     ++in_it1;
     ++in_it2;
     ++out_it;

     ++count;
   }

   return result;
}



template< typename TImageType >
void
SmallUtilityMethods<TImageType>::calculateBoundingBoxHavingHomogenousBackground(
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image, 
	typename SmallUtilityMethods<TImageType>::TImageIndexType& 
      bound_box_min_index,
   typename SmallUtilityMethods<TImageType>::TImageIndexType& 
      bound_box_max_index )
{
   typedef itk::ImageRegionConstIterator<TImageType> ConstIteratorType;

   const unsigned int Dimensions = 
      itk::GetImageDimension<TImageType>::ImageDimension;

   // initialize
   
   for (int i=0; i<Dimensions; ++i)
   {
      bound_box_min_index[i] = itk::NumericTraits<TImageIndexValueType>::max();
      bound_box_max_index[i] = itk::NumericTraits<TImageIndexValueType>::min();
   }

   typename TImageType::PixelType minimum_grey_val = 
      getMinimumVoxelValue( image );
      
   typename TImageType::IndexType index;

   ConstIteratorType in_it = ConstIteratorType ( image, 
      image->GetLargestPossibleRegion() );
      
   in_it.GoToBegin();
   while (!in_it.IsAtEnd())
   {
      if (in_it.Value() > minimum_grey_val)
      {
         index = in_it.GetIndex();

         for (int i=0; i<Dimensions; ++i)
         {
            if (index[i] < bound_box_min_index[i])
               bound_box_min_index[i] = index[i];

            if (index[i] > bound_box_max_index[i])
               bound_box_max_index[i] = index[i];
         }
      }

      ++in_it;
   }
}


template< typename TImageType >
typename SmallUtilityMethods<TImageType>::TImageTypePointer
SmallUtilityMethods<TImageType>::cropVolumeHavingHomogenousBackground(
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& 
      image_with_homogenous_background )
{
   typename TImageType::SizeType original_size =
      image_with_homogenous_background->GetLargestPossibleRegion().GetSize();
   
   std::cout << "Find a bounding box for image." << std::endl;

   typename TImageType::IndexType bound_box_min;
   typename TImageType::IndexType bound_box_max;

   const unsigned int Dimensions = 
      itk::GetImageDimension<TImageType>::ImageDimension;

   calculateBoundingBoxHavingHomogenousBackground( 
      image_with_homogenous_background, bound_box_min, bound_box_max );

   std::cout << "bound_box_min: " <<  bound_box_min << std::endl;
   std::cout << "bound_box_max: " <<  bound_box_max << std::endl;

   std::cout << "Create new image and copy the image according to bounding box."
      << std::endl;

   typedef itk::RegionOfInterestImageFilter< TImageType,
                                             TImageType > 
                                           ROIExtractionFilterType;
   typename ROIExtractionFilterType::Pointer roi_extracter = 
      ROIExtractionFilterType::New();

   typename TImageType::IndexType start;
   typename TImageType::SizeType size;
   for (int i=0; i<Dimensions; ++i)
   {
      start[i] = bound_box_min[i];
      size[i] = bound_box_max[i] - bound_box_min[i] + 1;
   }

   std::cout << "start: " << start << std::endl;
   std::cout << "size: " << size  << std::endl;

   // correct size to be able to downsample images up to a factor of 8 (2^3)
   for (int j=0; j<Dimensions; ++j)
   {
      typename TImageType::SizeValueType temp = original_size[j];
      
      if ( (size[j] % 8) != 0)
         size[j] += ( 8 - (size[j] % 8) );
         
      if (size[j] > temp)
         size[j] = temp;
   }

   std::cout << "corrected size: " << size  << std::endl;

   typename TImageType::RegionType desiredRegion;
   desiredRegion.SetSize(  size  );
   desiredRegion.SetIndex( start );

   roi_extracter->SetRegionOfInterest( desiredRegion );
   roi_extracter->SetInput( image_with_homogenous_background );

   ITK_EXCEPTION_CHECKED( "Extract ROI according to bounding box.", 
      roi_extracter->Update(), 0 );

   return roi_extracter->GetOutput();
}


template< typename TImageType >
typename SmallUtilityMethods<TImageType>::TImageIndexType
SmallUtilityMethods<TImageType>::calculateCenterInBinaryImage(
   const typename SmallUtilityMethods<TImageType>::TImageTypePointer& image )
{
   double center_x = 0.0;
   double center_y = 0.0;
   double center_z = 0.0;
      
   itk::ImageRegionConstIteratorWithIndex<TImageType> it( image, 
      image->GetLargestPossibleRegion() );
      
   it.Begin();

   unsigned long counter = 0L;
   
   while ( !it.IsAtEnd() )
   {
      if (it.Get() > 0)
      {
         typename TImageType::IndexType index = it.GetIndex();
         
         center_x += index[0];
         center_y += index[1];
         center_z += index[2];

         ++counter;
      }
      ++it;
   }

   center_x /= static_cast<double>( counter );
   center_y /= static_cast<double>( counter );
   center_z /= static_cast<double>( counter );

   typename TImageType::IndexType center;
   center[0] = ROUND_DBL< TImageIndexValueType >( center_x );
   center[1] = ROUND_DBL< TImageIndexValueType >( center_y );
   center[2] = ROUND_DBL< TImageIndexValueType >( center_z );
   
   return center;
}



#endif //  __SMALL_UTILITY_METHODS_TXX__