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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: main.cxx,v $
  Language:  C++
  Date:      $Date: 2006/01/15 04:28:36 $
  Version:   $Revision: 1.12 $

  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.

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

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkImageSliceIteratorWithIndex.h"
#include "itkOrthogonalBisectImageFilter.h"
#include "itkSphereSpatialFunction.h"
#include "itkFloodFilledSpatialFunctionConditionalIterator.h"
#include "itkImageRegionIterator.h"

int main()
{
  // This is meant to validate the the filter. A 3D image consisting of 2 spheres 
  // is generated. The filter is used to bisect the image such that each output image
  // contains one of the spheres. Each output image is iterated over to 
  // verify that the pixel values correspond to the input image.
 
  // The image dimension. Eventhough imageDimension is set to 3, a 2D image can still be used.
  const unsigned int imageDimension = 3;
  
  // Image typedefs.
  typedef unsigned char PixelType;
  typedef itk::Image<PixelType, imageDimension> ImageType;
  
  // Initialize the image reader.
  typedef itk::ImageFileReader<ImageType> ImageReaderType;
  ImageReaderType::Pointer imageReader = ImageReaderType::New();
  
  // Initialize the image writer.
  typedef itk::ImageFileWriter<ImageType> FileWriterType;
  FileWriterType::Pointer fileWriter = FileWriterType::New();
  
  // Initialize the OrthogonalBisectImageFilter.
  typedef itk::OrthogonalBisectImageFilter<ImageType> FilterType;
  FilterType::Pointer filter = FilterType::New();
  
  // Create a 100x100x100 image containing two spheres
  ImageType::Pointer synthImage = ImageType::New();
  ImageType::RegionType synthRegion;
  ImageType::SizeType synthSize = {{100, 100, 100}};
  ImageType::IndexType synthIndex;
  double synthSpacing[] = {1.0, 1.0, 1.0};
  double synthOrigin[] = {0, 0, 0};
  
  // Voxel intensities to be used
  unsigned char white = 255;
  unsigned char gray = 128;
  unsigned char black = 0;
  
  // Initialize the image
  synthIndex.Fill(0);
  synthRegion.SetSize(synthSize);
  synthRegion.SetIndex(synthIndex);
  synthImage->SetOrigin(synthOrigin);
  synthImage->SetSpacing(synthSpacing);
  synthImage->SetLargestPossibleRegion(synthRegion);
  synthImage->SetBufferedRegion(synthRegion);
  synthImage->SetRequestedRegion(synthRegion);
  synthImage->Allocate();
  synthImage->FillBuffer(black);
  
  // Put a white sphere in the "left half" of the image.
  typedef itk::SphereSpatialFunction<imageDimension, ImageType::PointType> SphereSpatialFuncType;
  SphereSpatialFuncType::Pointer sphereSpatialFunc = SphereSpatialFuncType::New();

  double sphere1Radius = 28;
  typedef SphereSpatialFuncType::InputType SphereCenterType;
  SphereCenterType sphereCenter;
  sphereCenter[0] = 30;
  sphereCenter[1] = 30;
  sphereCenter[2] = 45;
  sphereSpatialFunc->SetCenter(sphereCenter);
  sphereSpatialFunc->SetRadius(sphere1Radius);
  
  typedef ImageType::IndexType IndexType;
  typedef IndexType::IndexValueType IndexValueType;
  IndexType seedPos;
  seedPos[0] = static_cast<IndexValueType>(sphereCenter[0]); 
  seedPos[1] = static_cast<IndexValueType>(sphereCenter[1]); 
  seedPos[2] = static_cast<IndexValueType>(sphereCenter[2]); 

  itk::FloodFilledSpatialFunctionConditionalIterator<ImageType, SphereSpatialFuncType> 
	sfi = itk::FloodFilledSpatialFunctionConditionalIterator<ImageType,
	SphereSpatialFuncType>(synthImage, sphereSpatialFunc, seedPos);
   
  sfi.SetOriginInclusionStrategy();
  // Iterate through the entire image and set interior pixels to 255  
  for(; !sfi.IsAtEnd(); ++sfi)
    sfi.Set(white);
      
  // Put a small gray sphere in the "right half" of the image.
  SphereSpatialFuncType::Pointer sphere2SpatialFunc = SphereSpatialFuncType::New();

  double sphere2Radius = 15;
  SphereCenterType sphere2Center;
  sphere2Center[0] = 80;
  sphere2Center[1] = 80;
  sphere2Center[2] = 80;
  sphere2SpatialFunc->SetCenter(sphere2Center);
  sphere2SpatialFunc->SetRadius(sphere2Radius);
  
  IndexType seed2Pos;
  seed2Pos[0] = static_cast<IndexValueType>(sphere2Center[0]); 
  seed2Pos[1] = static_cast<IndexValueType>(sphere2Center[1]); 
  seed2Pos[2] = static_cast<IndexValueType>(sphere2Center[2]); 
  
  // Iterate through the entire image and set interior pixels to 128  
  itk::FloodFilledSpatialFunctionConditionalIterator<ImageType, SphereSpatialFuncType> 
	sfi2 = itk::FloodFilledSpatialFunctionConditionalIterator<ImageType,
	SphereSpatialFuncType>(synthImage, sphere2SpatialFunc, seed2Pos);
   
  sfi2.SetOriginInclusionStrategy();
  for(; !sfi2.IsAtEnd(); ++sfi2)
	sfi2.Set(gray);
    
  // Write the synthetic image to disk.
  fileWriter->SetInput(synthImage);
  fileWriter->SetFileName("synthTestImage.img");
  fileWriter->Update();
  
  // Bisect the image between the spheres along the x-axis.
  filter->SetInput(synthImage); // Set the input image
  filter->SetOrthogonalDirection(std::string("x")); // Set the bisect direction
  int bisectIndex = 60;
  filter->SetBisectIndex(bisectIndex); // Set the bisect index.
  filter->SetOutput1FillIntensity(gray);  // Set the first output image's fill intensity to gray
  filter->SetOutput2FillIntensity(white);  // Set the second output image's fill intensity to white

  // Try and update the filter. If that fails, throw an exception.
  try
	{
	filter->Update();
	}
  catch(itk::ExceptionObject & e)
	{
	std::cerr << "Exception detected: "  << e.GetDescription();
	return EXIT_FAILURE;
	}
	
  // Print the member variables of the filter.
  std::cerr << filter << std::endl;
  
  // Save the output images to disk.
  fileWriter->SetInput(filter->GetOutput(0));
  fileWriter->SetFileName("synthOutput1.img");
  fileWriter->Update();
  fileWriter->SetInput(filter->GetOutput(1));
  fileWriter->SetFileName("synthOutput2.img");
  fileWriter->Update();

  // Iterate the output images and verify that pixels were copied 
  // with the correct value to the correct location
  typedef itk::ImageRegionIterator<ImageType> ImageIteratorType;
  ImageIteratorType imageIt(filter->GetOutput(0), filter->GetOutput(0)->GetLargestPossibleRegion());

  ImageType::IndexType tempIndex;
  SphereSpatialFuncType::InputType sfPos;
  bool inSphere = 0;
  for(imageIt.GoToBegin(); !imageIt.IsAtEnd(); ++imageIt)
	{
	tempIndex = imageIt.GetIndex();
	for(unsigned int j = 0; j < imageDimension; j++)
	  sfPos[j] = tempIndex[j];
          
	inSphere = sphereSpatialFunc->Evaluate(sfPos);
	if(inSphere == 1 && imageIt.Get() != white)
	  {
	  // In the sphere, the value should be 255 (sfi.Set(white)).
	  std::cerr << "Failure 1: This pixel is inside the sphere and should have a value of " << (int)white 
				<< " but it is " << (int)imageIt.Get() << std::endl;
	  return EXIT_FAILURE;
	  }
	else if(inSphere == 0 && tempIndex[0] <= bisectIndex && imageIt.Get() != black) 
	  {
	  // In the first image division and outside the sphere, the value should be 0 (synthImage->FillBuffer(black)).
	  std::cerr << "Failure 2: This pixel is outside the sphere and from the first bisected region. "
				<< "Its value should be " << (int)black << ", but it is " << (int)imageIt.Get() << std::endl;
	  return EXIT_FAILURE;
	  }
	else if(tempIndex[0] > bisectIndex && imageIt.Get() != gray) 
	  {
	  // In the second image division, the value should be 128 (filter->SetOutput1FillIntensity(gray)).
	  std::cerr << "Failure 3: This pixel is from the second bisected region. "
				<< "Its intensity should be " << (int)gray << ", but it is " << (int)imageIt.Get() << std::endl;
	  return EXIT_FAILURE;
	  }
	}
  
  std::cerr << "The first output image is correct!" << std::endl;

  typedef itk::ImageRegionIterator<ImageType> ImageIteratorType;
  ImageIteratorType imageIt2(filter->GetOutput(1), filter->GetOutput(1)->GetLargestPossibleRegion());
  
  for(imageIt2.GoToBegin(); !imageIt2.IsAtEnd(); ++imageIt2)
	{
	tempIndex = imageIt2.GetIndex();
	for(unsigned int j = 0; j < imageDimension; j++)
	  sfPos[j] = tempIndex[j];
          
	inSphere = sphere2SpatialFunc->Evaluate(sfPos);
	if(inSphere == 1 && imageIt2.Get() != gray)
	  {
	  // In the sphere, the value should be 128 (sfi2.Set(gray)).
	  std::cerr << "Failure 4: This pixel is inside the sphere and should have a value of " << (int)gray 
				<< " but it is " << (int)imageIt2.Get() << std::endl;
	  return EXIT_FAILURE;
	  }
	else if(inSphere == 0 && tempIndex[0] > bisectIndex && imageIt2.Get() != black) 
	  {
	  // In the second image division and outside the sphere, the value should be 0 (synthImage->FillBuffer(black)).
	  std::cerr << "Failure 5: This pixel is outside the sphere and from the first bisected region. "
				<< "Its value should be " << (int)black << ", but it is " << (int)imageIt2.Get() << std::endl;
	  return EXIT_FAILURE;
	  }
	else if(tempIndex[0] <= bisectIndex && imageIt2.Get() != white) 
	  {
	  // In the first image division, the value should be 255 (filter->SetOutput2FillIntensity(white)).
	  std::cerr << "Failure 6: This pixel is from the second bisected region. "
				<< "Its intensity should be " << (int)white << ", but it is " << (int)imageIt2.Get() << std::endl;
	  return EXIT_FAILURE;
	  }
    }
    
  std::cerr << "The second output image is correct!" << std::endl;
	
  return EXIT_SUCCESS;
}