/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkOrthogonalSwath2DPathFilterTest.cxx,v $
  Language:  C++
  Date:      $Date: 2003-09-10 14:30:04 $
  Version:   $Revision: 1.6 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     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.

=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#include "itkPolyLineParametricPath.h"
#include "itkChainCodePath.h"
#include "itkFourierSeriesPath.h"
#include "itkPathToChainCodePathFilter.h"
#include "itkChainCodeToFourierSeriesPathFilter.h"

#include "itkImage.h"
#include "itkImageRegionIterator.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkDiscreteGaussianImageFilter.h"
#include "itkExtractOrthogonalSwath2DImageFilter.h"
#include "itkDerivativeImageFilter.h"

#include "itkOrthogonallyCorrected2DParametricPath.h"
#include "itkOrthogonalSwath2DPathFilter.h"

#include "itkPathToImageFilter.h"
#include "itkImageFileWriter.h"


int itkOrthogonalSwath2DPathFilterTest(int, char*[])
{
  typedef itk::Image<unsigned char, 2>                ImageType;         
  typedef itk::Image<double, 2>                       FloatImageType;
  typedef itk::PolyLineParametricPath<2>              InPathType;        
  typedef itk::ChainCodePath<2>                       ChainPathType;     
  typedef itk::FourierSeriesPath<2>                   FSPathType;        

  typedef InPathType::VertexType                      VertexType;        
  typedef InPathType::OffsetType                      OffsetType;        
  typedef InPathType::InputType                       InPathInputType;   

  typedef ImageType::IndexType                        IndexType;         
  typedef ImageType::SizeType                         SizeType;          

  typedef itk::OrthogonallyCorrected2DParametricPath  OutputPathType;    

  // pre-process the image
  typedef itk::RescaleIntensityImageFilter<ImageType,FloatImageType>      CastFilterType;
  typedef itk::DiscreteGaussianImageFilter<FloatImageType,FloatImageType> SmoothFilterType;
  typedef itk::ExtractOrthogonalSwath2DImageFilter<FloatImageType>        SwathFilterType;
  typedef itk::DerivativeImageFilter<FloatImageType,FloatImageType>       MeritFilterType;
  
  // pre-process the path
  typedef itk::PathToChainCodePathFilter<InPathType,ChainPathType>          PathFilter1Type;
  typedef itk::ChainCodeToFourierSeriesPathFilter<ChainPathType,FSPathType> PathFilter2Type;
  
  // test the filter
  typedef itk::OrthogonalSwath2DPathFilter<FSPathType,FloatImageType>       TestFilterType;
  
  // save the results
  typedef itk::PathToImageFilter< FSPathType,     ImageType >         Output1FilterType;
  typedef itk::PathToImageFilter< OutputPathType, ImageType >         Output2FilterType;
  typedef itk::RescaleIntensityImageFilter<FloatImageType,ImageType>  Output3FilterType;
  
  
  
  // Setup the path
  std::cout << "Making a square Path with v0 at (24,24) -> (24,104) -> (104,104) -> (104,24)" << std::endl;
  InPathType::Pointer inPath;
  VertexType        v;
  inPath = InPathType::New();
  v.Fill(24);
  inPath->AddVertex(v);
  v[0]=24;
  v[1]=104;
  inPath->AddVertex(v);
  v.Fill(104);
  inPath->AddVertex(v);
  v[0]=104;
  v[1]=24;
  inPath->AddVertex(v);
  v.Fill(24);
  inPath->AddVertex(v);
  
  // Setup the first path filter
  PathFilter1Type::Pointer pathFilter1 = PathFilter1Type::New();
  pathFilter1->SetInput(inPath);
  
  // Setup the second path filter
  PathFilter2Type::Pointer pathFilter2 = PathFilter2Type::New();
  pathFilter2->SetInput(pathFilter1->GetOutput());
  pathFilter2->SetNumberOfHarmonics(7); // make a nice, round, path for the swath
  
  
  
  // Setup the image
  std::cout << "Making a 64x64 black square centered in a 128x128 white image"<<std::endl;
  ImageType::Pointer  inImage = ImageType::New();
  IndexType start;
  start[0]=0;
  start[1]=0;
  ImageType::SizeType size;
  size[0]=128;
  size[1]=128;
  ImageType::RegionType region;
  region.SetSize(size);
  region.SetIndex(start);
  inImage->SetRegions(region);
  double spacing[ ImageType::ImageDimension ];
  spacing[0]=1.0;
  spacing[1]=1.0;
  inImage->SetSpacing(spacing);
  inImage->Allocate();
  typedef itk::ImageRegionIterator<ImageType> ImageItType;
  ImageItType it( inImage, inImage->GetRequestedRegion() );
  it.GoToBegin();
  IndexType pixelIndex;
  while( !it.IsAtEnd() )
    {
    pixelIndex = it.GetIndex();
    if( pixelIndex[0] >= int(size[0]/4) && pixelIndex[0] < int(size[0]*3/4) &&
        pixelIndex[1] >= int(size[1]/4) && pixelIndex[1] < int(size[1]*3/4) )
      it.Set(0);
    else
      it.Set(255);
    ++it;
    }
  
  // Cast the input image into a double image
  CastFilterType::Pointer castFilter = CastFilterType::New();
  castFilter->SetInput( inImage );
  castFilter->SetOutputMinimum(0);
  castFilter->SetOutputMaximum(1.0);

  // Smooth the (double pixel type) input image
  SmoothFilterType::Pointer smoothFilter = SmoothFilterType::New();
  smoothFilter->SetInput( castFilter->GetOutput() );
  double gaussianVariance = 1.0;
  // We want a fast 3x3 kernel.  Gausian Operator will not truncate its kernel
  // width to any less than a 5x5 kernel (kernel width of 3 for 1 center pixel +
  // 2 edge pixels).  However, GausianOperator always uses at least a 3x3
  // kernel, and so setting the maximum error to 1.0 (no limit) will make it
  // stop growing the kernel at the desired 3x3 size.
  double maxError = 0.9;  
  smoothFilter->SetUseImageSpacingOff();
  smoothFilter->SetVariance( gaussianVariance );
  smoothFilter->SetMaximumError( maxError );
  
  // Extract the swath image
  SwathFilterType::Pointer swathFilter = SwathFilterType::New();
  swathFilter->SetImageInput( smoothFilter->GetOutput() );
//  swathFilter->SetImageInput( castFilter->GetOutput() );
  swathFilter->SetPathInput( pathFilter2->GetOutput() );
  size[0]=512;
  size[1]=16*2+1; // the top 1 and bottom 1 rows are dropped when smoothing
  swathFilter->SetSize(size);
  
  // Find the vertical gradient of the swath image
  MeritFilterType::Pointer meritFilter = MeritFilterType::New();
  meritFilter->SetInput( swathFilter->GetOutput() );
  meritFilter->SetOrder( 1 ); // first partial derivative
  meritFilter->SetDirection( 1 ); // d/dy
  
  
  
  // Setup the test filter
  std::cerr << "Creating the test filter" << std::endl;
  TestFilterType::Pointer testFilter = TestFilterType::New();
  std::cerr << "Setting up the test filter" << std::endl;
  testFilter->SetPathInput(pathFilter2->GetOutput());
  testFilter->SetImageInput(meritFilter->GetOutput());
  OutputPathType::Pointer outPath = testFilter->GetOutput();
  
  // Setup the input & output path images
  std::cerr << "Setting up the output path image" << std::endl;
  Output1FilterType::Pointer output1Filter = Output1FilterType::New();
  Output2FilterType::Pointer output2Filter = Output2FilterType::New();
  size[0]=128;
  size[1]=128;
  output1Filter->SetSize(size);  // same size as the input image
  output2Filter->SetSize(size);  // same size as the input image
  output1Filter->SetPathValue(255);
  output2Filter->SetPathValue(255);
  output1Filter->SetInput( testFilter->GetPathInput() );
  output2Filter->SetInput( testFilter->GetOutput() );
  ImageType::Pointer inPathImage = output1Filter->GetOutput();
  ImageType::Pointer outImage =    output2Filter->GetOutput();
  
  // Setup the swath merit output image
  Output3FilterType::Pointer output3Filter = Output3FilterType::New();
  output3Filter->SetInput( meritFilter->GetOutput() );
  output3Filter->SetOutputMinimum(0);
  output3Filter->SetOutputMaximum(255);
  ImageType::Pointer swathMeritImage = output3Filter->GetOutput();
  
  
  // Testing PrintSelf
  std::cerr << testFilter << std::endl;
  
  // Update the pipeline
  std::cerr << "Running the Pipeline: ";
  outImage->Update();
  std::cerr << "[DONE]" << std::endl;
/*  
  // Save the output images
  itk::ImageFileWriter<ImageType>::Pointer writer
    = itk::ImageFileWriter<ImageType>::New();
  writer->SetInput( output1Filter->GetOutput() );
  writer->SetFileName( "OrthSwathFilterPathIn.png" );
  writer->Write();
  writer->SetInput( output2Filter->GetOutput() );
  writer->SetFileName( "OrthSwathFilterPathOut.png" );
  writer->Write();
  writer->SetInput( output3Filter->GetOutput() );
  writer->SetFileName( "OrthSwathFilterSwathMerit.png" );
  writer->Write();
*/  
  return EXIT_SUCCESS;
}