text: /*=========================================================================
   *
   *  Copyright Insight Software Consortium
   *
   *  Licensed under the Apache License, Version 2.0 (the "License");
   *  you may not use this file except in compliance with the License.
   *  You may obtain a copy of the License at
   *
   *         http://www.apache.org/licenses/LICENSE-2.0.txt
   *
   *  Unless required by applicable law or agreed to in writing, software
   *  distributed under the License is distributed on an "AS IS" BASIS,
   *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   *  See the License for the specific language governing permissions and
   *  limitations under the License.
   *
   *=========================================================================*/
  
  #include "itkImageFileReader.h"
  #include "itkImageFileWriter.h"
  #include "itkDiscreteGaussianDerivativeImageFunction.h"
  #include "itkRescaleIntensityImageFilter.h"
  
  template < int VDimension >
  int itkDiscreteGaussianDerivativeImageFunctionTestND( int argc, char* argv[] )
  {
  
    // Verify the number of parameters in the command line
    if( argc < 5 )
      {
      std::cerr << "Usage: " << std::endl;
      std::cerr << argv[0] << "inputFileName outputFileName order sigma (maximum_error) (maximum_kernel_width)" << std::endl;
      return EXIT_FAILURE;
      }
  
    // Define the dimension of the images
    const unsigned int Dimension = VDimension;
    typedef float                            PixelType;
    typedef itk::Image<PixelType, Dimension> ImageType;
  
    // Read input
    typedef itk::ImageFileReader< ImageType > ReaderType;
    typename ReaderType::Pointer reader = ReaderType::New();
    reader->SetFileName( argv[1] );
    try
      {
      reader->Update();
      }
    catch ( itk::ExceptionObject &err)
      {
      std::cout << "ExceptionObject caught !" << std::endl;
      std::cout << err << std::endl;
      return EXIT_FAILURE;
      }
    ImageType *inputImage = reader->GetOutput();
  
    // Create image for storing result
    typename ImageType::Pointer output = ImageType::New();
    output->SetSpacing( inputImage->GetSpacing() );
    output->SetOrigin( inputImage->GetOrigin() );
    output->SetDirection( inputImage->GetDirection() );
    output->SetLargestPossibleRegion( inputImage->GetLargestPossibleRegion() );
    output->SetRequestedRegion( inputImage->GetRequestedRegion() );
    output->SetBufferedRegion( inputImage->GetBufferedRegion() );
    output->Allocate();
    output->FillBuffer( itk::NumericTraits<PixelType>::ZeroValue() );
  
    // Setup operator parameters
    unsigned int order[Dimension];
    for (unsigned int i = 0; i < Dimension; i++)
      {
      order[i] = atoi( argv[3] );
      }
  
    double variance = atof( argv[4] );
    variance *= variance;
  
    double maxError = 0.001;
    unsigned int maxKernelWidth = 100;
    if( argc == 6 )
      {
      maxError = atof( argv[5] );
      }
    else if( argc > 6 )
      {
      maxError = atof( argv[5] );
      maxKernelWidth = atoi( argv[6] );
      }
  
    // Create function
    typedef itk::DiscreteGaussianDerivativeImageFunction< ImageType, PixelType >
      GaussianDerivativeImageFunctionType;
    typename GaussianDerivativeImageFunctionType::Pointer function =
      GaussianDerivativeImageFunctionType::New();
    function->SetInputImage( inputImage );
    function->SetMaximumError( maxError );
    function->SetMaximumKernelWidth( maxKernelWidth );
    function->SetVariance( variance );
    function->SetOrder( order );
    function->SetNormalizeAcrossScale( true );
    function->SetUseImageSpacing( true );
    function->SetInterpolationMode( GaussianDerivativeImageFunctionType::NearestNeighbourInterpolation );
    function->Initialize( );
  
    // Step over input and output images
    typedef itk::ImageRegionConstIterator< ImageType > ConstIteratorType;
    typedef itk::ImageRegionIterator< ImageType >      IteratorType;
  
    ConstIteratorType it ( inputImage, inputImage->GetRequestedRegion() );
    it.GoToBegin();
    IteratorType out( output, output->GetRequestedRegion() );
    out.GoToBegin();
  
    typedef typename GaussianDerivativeImageFunctionType::PointType  PointType;
    PointType point;
    typedef typename GaussianDerivativeImageFunctionType::ContinuousIndexType ContinuousIndexType;
    ContinuousIndexType cindex;
    const unsigned long nop = inputImage->GetRequestedRegion().GetNumberOfPixels();
    unsigned long pixelNumber = 0;
    while( !it.IsAtEnd() )
      {
      // To test all available Evaluate functions, we split it in three parts.
      if ( pixelNumber < nop / 3 )
        {
        out.Set( function->EvaluateAtIndex( it.GetIndex() ) );
        }
      else if ( pixelNumber < nop * 2 / 3 )
        {
        inputImage->TransformIndexToPhysicalPoint( it.GetIndex(), point );
        out.Set( function->Evaluate( point ) );
        }
      else
        {
        inputImage->TransformIndexToPhysicalPoint( it.GetIndex(), point );
        inputImage->TransformPhysicalPointToContinuousIndex( point, cindex );
        out.Set( function->EvaluateAtContinuousIndex( cindex ) );
        }
      ++it;
      ++out;
      ++pixelNumber;
      }
  
    // Rescale output
    typedef unsigned char                                                   OutputPixelType;
    typedef itk::Image< OutputPixelType, Dimension >                        OutputImageType;
    typedef itk::RescaleIntensityImageFilter< ImageType, OutputImageType >  RescaleType;
  
    typename RescaleType::Pointer rescaler = RescaleType::New();
    rescaler->SetInput( output );
    rescaler->SetOutputMinimum( itk::NumericTraits<OutputPixelType>::min() );
    rescaler->SetOutputMaximum( itk::NumericTraits<OutputPixelType>::max() );
  
    // Write output
    typedef itk::ImageFileWriter< OutputImageType > WriterType;
    typename WriterType::Pointer writer = WriterType::New();
    writer->SetFileName( argv[2] );
    writer->SetInput( rescaler->GetOutput() );
    try
      {
      writer->Update();
      }
    catch ( itk::ExceptionObject &err)
      {
      std::cout << "ExceptionObject caught !" << std::endl;
      std::cout << err << std::endl;
      return EXIT_FAILURE;
      }
  
    // Test some functions
    typedef typename GaussianDerivativeImageFunctionType::VarianceArrayType VarianceArrayType;
    VarianceArrayType varReturned = function->GetVariance();
    for ( unsigned int i = 0; i < Dimension; ++i )
    {
      if ( varReturned[ i ] != variance )
      {
        std::cout << "GetVariance()[" << i << "] failed. Expected: "
          << variance
          << " but got: "
          << varReturned[ i ] << std::endl;
        return EXIT_FAILURE;
      }
    }
    typedef typename GaussianDerivativeImageFunctionType::OrderArrayType  OrderArrayType;
    OrderArrayType orderReturned = function->GetOrder();
    for ( unsigned int i = 0; i < Dimension; ++i )
    {
      if ( orderReturned[ i ] != order[ i ] )
      {
        std::cout << "GetOrder()[" << i << "] failed. Expected: "
          << order[ i ]
          << " but got: "
          << orderReturned[ i ] << std::endl;
        return EXIT_FAILURE;
      }
    }
    if ( function->GetMaximumError() != maxError )
    {
      std::cout << "GetMaximumError failed. Expected: "
        << maxError
        << " but got: "
        << function->GetMaximumError() << std::endl;
      return EXIT_FAILURE;
    }
    if ( function->GetNormalizeAcrossScale() != true )
    {
      std::cout << "GetNormalizeAcrossScale failed. Expected: "
        << true
        << " but got: "
        << function->GetNormalizeAcrossScale() << std::endl;
      return EXIT_FAILURE;
    }
    if ( function->GetUseImageSpacing() != true )
    {
      std::cout << "GetUseImageSpacing failed. Expected: "
        << true
        << " but got: "
        << function->GetUseImageSpacing() << std::endl;
      return EXIT_FAILURE;
    }
    if ( function->GetMaximumKernelWidth() != maxKernelWidth )
    {
      std::cout << "GetMaximumKernelWidth failed. Expected: "
        << maxKernelWidth
        << " but got: "
        << function->GetMaximumKernelWidth() << std::endl;
      return EXIT_FAILURE;
    }
    if ( function->GetInterpolationMode() != GaussianDerivativeImageFunctionType::NearestNeighbourInterpolation )
    {
      std::cout << "GetInterpolationMode failed. Expected: "
        << GaussianDerivativeImageFunctionType::NearestNeighbourInterpolation
        << " but got: "
        << function->GetInterpolationMode() << std::endl;
      return EXIT_FAILURE;
    }
  
    // Call PrintSelf.
    function->Print( std::cout );
  
    return EXIT_SUCCESS;
  }
  
  int itkDiscreteGaussianDerivativeImageFunctionTest(int argc, char* argv[] )
  {
    return itkDiscreteGaussianDerivativeImageFunctionTestND< 2 >( argc, argv );
  }
  