#include "itkNumericTraits.h"

#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkImageSeriesReader.h"
#include "itkImageFileWriter.h"

#include "itkMetaDataDictionary.h"

#include "itkNaryElevateImageFilter.h"
#include "itkUnaryRetractImageFilter.h"

#include "itkStreamingImageFilter.h"

#include "functors.cxx"

int main(int argc, char* argv[])
{
  // check arguments
  if (argc!=2)
  {
    std::cerr << "Usage: executable <dicomDirectory>" << std::endl;
    return EXIT_FAILURE;
  }

  // declare images
  typedef signed short PixelType;
  typedef itk::Image<PixelType, 2> Image2DType;
  typedef itk::Image<PixelType, 3> Image3DType;
  typedef itk::Image<PixelType, 4> Image4DType;

  // create elevation filter
  typedef itk::NaryElevateImageFilter<Image3DType, Image4DType> NaryElevateFilterType;
  NaryElevateFilterType::Pointer elevateFilter=NaryElevateFilterType::New();

  // create name generator and attach to reader
  typedef itk::GDCMSeriesFileNames NamesGeneratorType;
  NamesGeneratorType::Pointer nameGenerator=NamesGeneratorType::New();
  nameGenerator->SetUseSeriesDetails(true);
  nameGenerator->AddSeriesRestriction("0020|0012"); // acquisition number
  nameGenerator->SetDirectory(argv[1]);

  // get series IDs
  typedef std::vector<std::string> SeriesIdContainer;
  const SeriesIdContainer &seriesUID=nameGenerator->GetSeriesUIDs();

  // create reader array
  typedef itk::ImageSeriesReader<Image3DType> ReaderType;
  ReaderType::Pointer *reader=new ReaderType::Pointer[seriesUID.size()];

  // declare series iterators
  SeriesIdContainer::const_iterator seriesItr=seriesUID.begin();
  SeriesIdContainer::const_iterator seriesEnd=seriesUID.end();
  int seriesNum=0;

  // connect each series to elevation filter
  while (seriesItr!=seriesEnd)
  {
    reader[seriesNum]=ReaderType::New();

    typedef itk::GDCMImageIO ImageIOType;
    ImageIOType::Pointer dicomIO=ImageIOType::New();
    reader[seriesNum]->SetImageIO(dicomIO);

    typedef std::vector<std::string> FileNamesContainer;
    FileNamesContainer fileNames;

    fileNames=nameGenerator->GetFileNames(seriesItr->c_str());
    reader[seriesNum]->SetFileNames(fileNames);

    elevateFilter->SetInput(seriesNum,reader[seriesNum]->GetOutput());
    reader[seriesNum]->ReleaseDataFlagOn();

    reader[seriesNum]->Update();
    reader[seriesNum]->GetOutput()->SetMetaDataDictionary(dicomIO->GetMetaDataDictionary());

    seriesNum++;

    ++seriesItr;
  }

  // get statistics from elevate filter
  elevateFilter->Update();
  float dataMin=elevateFilter->GetDataMin();
  float dataMax=elevateFilter->GetDataMax();
  float dataMean=elevateFilter->GetDataMean();
  float noiseThres=elevateFilter->GetNoiseThres();

  // create retraction filter #1
  typedef itk::UnaryRetractImageFilter<Image4DType, Image3DType, itk::Functor::CalcMin<typename Image4DType::PixelType, typename Image3DType::PixelType>> UnaryRetractFilterType1;
  UnaryRetractFilterType1::Pointer retractFilter1=UnaryRetractFilterType1::New();

  // create retraction filter #2
  typedef itk::UnaryRetractImageFilter<Image4DType, Image3DType, itk::Functor::CalcMinDev<typename Image4DType::PixelType, typename Image3DType::PixelType>> UnaryRetractFilterType2;
  UnaryRetractFilterType2::Pointer retractFilter2=UnaryRetractFilterType2::New();

  // connect elevation to retract filters
  retractFilter1->SetInput(elevateFilter->GetOutput());
  retractFilter2->SetInput(elevateFilter->GetOutput());

  // write data to multiframe DICOM
  typedef itk::ImageFileWriter<Image3DType> WriterType;
  WriterType::Pointer writer=WriterType::New();
  writer->SetFileName("out-MIPt.dcm");
  writer->SetInput(retractFilter1->GetOutput());
  writer->Update();
  retractFilter1=NULL;
  writer->SetFileName("out-WI.dcm");
  writer->SetInput(retractFilter2->GetOutput());
  writer->Update();
  retractFilter2=NULL;

  return EXIT_SUCCESS;
}