#ifndef _itkUnaryRetractImageFilter_txx
#define _itkUnaryRetractImageFilter_txx

#include "itkUnaryRetractImageFilter.h"

#include "itkImageLinearConstIteratorWithIndex.h"

namespace itk
{

template <class TInputImage, class TOutputImage, class TFunction>
UnaryRetractImageFilter<TInputImage, TOutputImage, TFunction>::
UnaryRetractImageFilter() {}

// check input
template <class TInputImage, class TOutputImage, class TFunction>
void
UnaryRetractImageFilter<TInputImage, TOutputImage, TFunction>::
GenerateOutputInformation()
{
  const unsigned int dimI=InputImageDimension;
  const unsigned int dimO=OutputImageDimension;

  if (dimI<1 || dimO<dimI-1) itkExceptionMacro(<< "inadequate output image dimension");

  OutputImagePointer outputPtr = this->GetOutput(0);

  InputImagePointer inputPtr =
    dynamic_cast<InputImageType *>(ProcessObject::GetInput(0));

  InputImageRegionType regionI=inputPtr->GetRequestedRegion();

  InputImageSizeType    sizeI    = regionI.GetSize();
  InputImageIndexType   indexI   = regionI.GetIndex();
  InputImageOriginType  originI  = inputPtr->GetOrigin();
  InputImageSpacingType spacingI = inputPtr->GetSpacing();

  OutputImageSizeType    sizeO;
  OutputImageIndexType   indexO;
  OutputImageOriginType  originO;
  OutputImageSpacingType spacingO;

  for (int i=0; i<dimI-1; i++)
  {
    sizeO[i]    = sizeI[i];
    indexO[i]   = indexI[i];
    originO[i]  = originI[i];
    spacingO[i] = spacingI[i];
  }

  for (int i=dimI; i<dimO; i++)
  {
    sizeO[i]    = 1;
    indexO[i]   = 0;
    originO[i]  = 0;
    spacingO[i] = 1;
  }

  OutputImageRegionType regionO;
  regionO.SetSize(sizeO);
  regionO.SetIndex(indexO);

  outputPtr->SetOrigin(originO);
  outputPtr->SetSpacing(spacingO);

  outputPtr->SetRegions(regionO);

  // pass tags:

  InputImagePointer inputPtr0 =
    dynamic_cast<InputImageType *>(ProcessObject::GetInput(0));

  outputPtr->SetMetaDataDictionary(inputPtr0->GetMetaDataDictionary());
}

// retract input
template <class TInputImage, class TOutputImage, class TFunction>
void
UnaryRetractImageFilter<TInputImage, TOutputImage, TFunction>::
GenerateData()
{
  const unsigned int dimI=InputImageDimension;
  const unsigned int dimO=OutputImageDimension;

  OutputImagePointer outputPtr = this->GetOutput(0);

  InputImagePointer inputPtr =
    dynamic_cast<InputImageType *>(ProcessObject::GetInput(0));

  inputPtr->Update();

  outputPtr->Allocate();

  const unsigned int runLength=inputPtr->GetBufferedRegion().GetSize()[dimI-1];

  // get time point list from meta dictionary:

  typedef MetaDataDictionary DictionaryType;
  DictionaryType &metaDict=inputPtr->GetMetaDataDictionary();

  std::string internalId="___internal-4D-filter-time-point-list"; // internal tag
  MetaDataDictionary::ConstIterator tagItr=metaDict.Find(internalId);

  ExceptionObject ex("invalid time point list");

  Array<double> timePoints;
  timePoints.SetSize(runLength);

  if (tagItr!=metaDict.End())
  {
    typedef MetaDataObject< std::list<double> > MetaDataListType;

    MetaDataObjectBase::Pointer entry = tagItr->second;
    MetaDataListType::Pointer entryvalue = dynamic_cast<MetaDataListType *>( entry.GetPointer() );

    // check whether or not the type of the entry value is correct
    if (entryvalue)
    {
      std::list<double> list;
      list=entryvalue->GetMetaDataObjectValue();

      if (list.size()!=runLength) throw ex;

      std::list<double>::const_iterator listItr=list.begin();

      for (int k=0; k<runLength; k++)
      {
        timePoints[k]=*listItr;
        listItr++;
      }
    }
    else throw ex;
  }
  else throw ex;

  typedef ImageLinearConstIteratorWithIndex<InputImageType> IteratorType;

  // process input by retrieving one time series after another:

  IteratorType itr(inputPtr, inputPtr->GetBufferedRegion());
  itr.SetDirection(dimI-1); // walk along highest dimension
  itr.GoToBegin();

  InputImageIndexType indexI;
  OutputImageIndexType indexO;

  for (int i=dimI-1; i<dimO; i++) indexO[i]=0;

  NaryArrayType naryInputArray;
  naryInputArray.SetSize(runLength);

  while(!itr.IsAtEnd())
  {
    itr.GoToBeginOfLine();

    unsigned int actualRun=0;

    // get actual time series
    while(!itr.IsAtEndOfLine())
    {
      naryInputArray[actualRun++] = itr.Get();
      ++itr;
    }

    indexI=itr.GetIndex();

    for (int i=0; i<dimI-1; i++) indexO[i]=indexI[i];

    // push actual time series through retract functor and store result
    outputPtr->SetPixel(indexO, static_cast<OutputImagePixelType>(m_Functor(naryInputArray,timePoints)));

    itr.NextLine();
  }
}

}

#endif