// -------------------------------------------------------------------------
// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
// -------------------------------------------------------------------------

#ifndef __cpExtensions__Algorithms__DiscontinuityMapImageFilter__hxx__
#define __cpExtensions__Algorithms__DiscontinuityMapImageFilter__hxx__

#include <cmath>
#include <itkConstNeighborhoodIterator.h>
#include <itkImageRegionIterator.h>
#include <itkNeighborhoodAlgorithm.h>
#include <itkZeroFluxNeumannBoundaryCondition.h>

// -------------------------------------------------------------------------
template< class _TInputImage >
cpExtensions::Algorithms::DiscontinuityMapImageFilter< _TInputImage >::
DiscontinuityMapImageFilter( )
  : Superclass( ),
    m_UseSquareRoot( false )
{
}

// -------------------------------------------------------------------------
template< class _TInputImage >
cpExtensions::Algorithms::DiscontinuityMapImageFilter< _TInputImage >::
~DiscontinuityMapImageFilter( )
{
}

// -------------------------------------------------------------------------
template< class _TInputImage >
void cpExtensions::Algorithms::DiscontinuityMapImageFilter< _TInputImage >::
ThreadedGenerateData( const TRegion& region, itk::ThreadIdType threadId )
{
  unsigned int i;

  itk::ZeroFluxNeumannBoundaryCondition< TInputImage > nbc;
  itk::ConstNeighborhoodIterator< TInputImage > bit;
  itk::ImageRegionIterator< TOutputImage > it;

  // Allocate output
  typename TOutputImage::Pointer output = this->GetOutput( );
  typename TInputImage::ConstPointer input = this->GetInput( );

  // Find the data-set boundary "faces"
  typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< TInputImage >::FaceListType faceList;
  itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< TInputImage > bC;
  faceList = bC( input, region, this->GetRadius( ) );

  typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< TInputImage >::FaceListType::iterator fit;

  TOutputPixel sum1, sum2;
  TOutputPixel neighborhoodSize;

  // Process each of the boundary faces.  These are N-d regions which border
  // the edge of the buffer.
  for( fit = faceList.begin( ); fit != faceList.end( ); ++fit )
  {
    bit = itk::ConstNeighborhoodIterator< TInputImage >( this->GetRadius( ), input, *fit );
    neighborhoodSize = bit.Size( );
    it = itk::ImageRegionIterator< TOutputImage >( output, *fit );
    bit.OverrideBoundaryCondition( &nbc );
    bit.GoToBegin( );
    while( !bit.IsAtEnd( ) )
    {
      sum1 = itk::NumericTraits< TOutputPixel >::ZeroValue( );
      sum2 = itk::NumericTraits< TOutputPixel >::ZeroValue( );
      for ( i = 0; i < neighborhoodSize; ++i )
      {
        TOutputPixel v = TOutputPixel( bit.GetPixel( i ) );
        sum1 += v;
        sum2 += v * v;
      }
      TOutputPixel mean = ( sum1 / TOutputPixel( neighborhoodSize ) );
      TOutputPixel var =
        ( sum2 - ( ( sum1 * sum1 ) / TOutputPixel( neighborhoodSize ) ) ) /
        ( TOutputPixel( neighborhoodSize ) - TOutputPixel( 1 ) );
      if( this->m_UseSquareRoot )
        it.Set( std::sqrt( double( var ) ) );
      else
        it.Set( var );
      ++bit;
      ++it;

    } // elihw

  } // rof
}

#endif // __cpExtensions__Algorithms__DiscontinuityMapImageFilter__hxx__

// eof - $RCSfile$