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

  Copyright (c) Kitware, Inc.
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/VolViewCopyright.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 notice for more information.

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

#include "itkPhantomLocationCalculator.h"
#include "itkImageFileReader.h"
#include <iostream>

int main( int argc, char *argv[] )
{
  if (argc < 10)
    {
    std::cerr << "Usage args : PhantomCTData.mha NTeflonPoints NDelrinPoints TeflonPoint_x TeflonPoint_y TeflonPoint_z .....  DelrinPoint1_x DelrinPoint1_y DelrinPoint1_z ......" << std::endl;
    return EXIT_FAILURE;
    }

  const unsigned int ImageDimension = 3;
  typedef itk::Image< short, ImageDimension > ImageType;
  typedef itk::PhantomLocationCalculator< ImageType > CalculatorType;
  typedef CalculatorType::PointType PointType;
  typedef CalculatorType::PointContainerType PointContainerType;

  itk::ImageFileReader< ImageType >::Pointer reader = 
    itk::ImageFileReader< ImageType >::New();
  reader->SetFileName( argv[1] );
  reader->Update();

  CalculatorType::Pointer calculator = CalculatorType::New();
  calculator->SetImage( reader->GetOutput() );
  calculator->SetLowerThreshold( 20 );
  calculator->SetUpperThreshold( 1200 );
  calculator->SetSearchRadius(80);

  PointType point;
  const unsigned int nTeflonPoints = atoi(argv[2]);
  const unsigned int nDelrinPoints = atoi(argv[3]);
  if (argc < (nTeflonPoints + nDelrinPoints + 4))
    {
    std::cerr << "Expected " << (nTeflonPoints + nDelrinPoints + 4) << " arguments, but got " << argc << " arguments." << std::endl;
    return EXIT_FAILURE;
    }

  PointContainerType teflonPoints(nTeflonPoints), delrinPoints(nDelrinPoints);
  
  for (unsigned int i = 0; i < nTeflonPoints; i++)
    {
    for (unsigned int j = 0; j < ImageDimension; j++)
      {
      teflonPoints[i][j] = atof(argv[4 + ImageDimension*i + j]);
      }
    std::cout << "Adding teflon point: " << teflonPoints[i] << std::endl;
    }

  for (unsigned int i = 0; i < nDelrinPoints; i++)
    {
    for (unsigned int j = 0; j < ImageDimension; j++)
      {
      delrinPoints[i][j] = atof(argv[4 + ImageDimension*i + j + nTeflonPoints*ImageDimension]);
      }
    std::cout << "Adding delrin point: " << delrinPoints[i] << std::endl;
    }
  
  calculator->SetUnOrderedTeflonPositions( teflonPoints );
  calculator->SetUnOrderedDelrinPositions( delrinPoints );

  calculator->Update();

  PointContainerType ordererdTePoints = calculator->GetOrderedTeflonPositions();
  PointContainerType ordererdDePoints = calculator->GetOrderedDelrinPositions();

  if (ordererdDePoints.size() != ordererdTePoints.size())
    {
    std::cerr << "Must be equal in pairs." << std::endl;
    return EXIT_FAILURE;
    }

  for (unsigned int i = 0; i < ordererdTePoints.size(); i++)
    {
    std::cout << "Teflon: " << ordererdTePoints[i] << " Delrin: " << ordererdDePoints[i] << std::endl;
    }
  
  return EXIT_SUCCESS;
}