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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: OrthogonalBisectImageFilter.txx,v $
  Language:  C++
  Date:      $Date: 2006/06/27 19:08:06 $
  Version:   $Revision: 1.00 $

  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.

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

#ifndef _itkOrthogonalBisectImageFilter_txx
#define _itkOrthogonalBisectImageFilter_txx

#include "itkOrthogonalBisectImageFilter.h"

namespace itk
{

/**
 * Constructor
 */
template<class TImage>
OrthogonalBisectImageFilter<TImage>
::OrthogonalBisectImageFilter()
{
  this->SetNumberOfRequiredInputs(1);
  this->SetNumberOfRequiredOutputs(2);

  m_OrthogonalDirection = std::string("x");
  m_BisectIndex = -1;
  m_Output1FillIntensity = 0;
  m_Output2FillIntensity = 0;  
}

/**
 * Return output images
 */
template<class TImage>
typename OrthogonalBisectImageFilter<TImage>::ImageType * 
OrthogonalBisectImageFilter<TImage>
::GetOutput(unsigned int whichOutput)
{
  if(whichOutput < this->GetNumberOfRequiredOutputs())
    return dynamic_cast<ImageType *>(this->ProcessObject::GetOutput(whichOutput));
}

/**
 * Print member variables
 */
template <class TImage>
void
OrthogonalBisectImageFilter<TImage>
::PrintSelf(std::ostream& os, Indent indent) const
{
  Superclass::PrintSelf(os,indent);
  
  os << std::endl;
  os << indent << "OrthogonalDirection: " << m_OrthogonalDirection << std::endl;
  os << indent << "BisectIndex: " << m_BisectIndex << std::endl;
  os << indent << "Output1FillIntensity: " << m_Output1FillIntensity << std::endl;
  os << indent << "Output2FillIntensity: " << m_Output2FillIntensity << std::endl;
}

/**
 * Go
 */
template <class TImage>
void
OrthogonalBisectImageFilter<TImage>
::GenerateData()
{
  // Set pointer to input image
  typename Superclass::InputImageConstPointer inputImage = this->GetInput();
 
  // Instantiate first output image
  ImagePointer outputImage1 = ImageType::New();
  this->SetNthOutput(0, outputImage1.GetPointer());
  outputImage1->SetLargestPossibleRegion(inputImage->GetLargestPossibleRegion());
  outputImage1->SetBufferedRegion(inputImage->GetBufferedRegion());
  outputImage1->SetRequestedRegion(inputImage->GetRequestedRegion());
  outputImage1->Allocate();
  outputImage1->FillBuffer(m_Output1FillIntensity);
  
  // Instantiate second output image
  ImagePointer outputImage2 = ImageType::New();
  this->SetNthOutput(1, outputImage2.GetPointer());
  outputImage2->SetLargestPossibleRegion(inputImage->GetLargestPossibleRegion());
  outputImage2->SetBufferedRegion(inputImage->GetBufferedRegion());
  outputImage2->SetRequestedRegion(inputImage->GetRequestedRegion());
  outputImage2->Allocate();
  outputImage2->FillBuffer(m_Output2FillIntensity);
  
  // Establish the vectors (directions[0] and directions[1]) that define 
  // the bisecting line or plane from the direction orthogonal to it (directions[2])
  unsigned int directions[3];
  if(m_OrthogonalDirection == std::string("x"))
    {
    directions[0] = 1;
    directions[1] = 2;
    directions[2] = 0;
    }
  else if(m_OrthogonalDirection == std::string("y"))
    { 
    directions[0] = 0;
    directions[1] = 2;
    directions[2] = 1;
    }
  else if(m_OrthogonalDirection == std::string("z"))
    {
    directions[0] = 0;
    directions[1] = 1;
    directions[2] = 2;
    }
  else
    {
	itk::ExceptionObject exceptionObject;
	exceptionObject.SetDescription("Invalid orthogonal direction!\n");
	throw exceptionObject;
    }
	
  // Initialize a default value for m_BisectIndex
  if(m_BisectIndex == -1 && m_OrthogonalDirection == std::string("x"))
    m_BisectIndex = (int) (0.5 * this->GetInput()->GetLargestPossibleRegion().GetSize(0) - 1);
  else if(m_BisectIndex == -1 && m_OrthogonalDirection == std::string("y"))
    m_BisectIndex = (int) (0.5 * this->GetInput()->GetLargestPossibleRegion().GetSize(1) - 1);
  else if(m_BisectIndex == -1 && m_OrthogonalDirection == std::string("z"))
    m_BisectIndex = (int) (0.5 * this->GetInput()->GetLargestPossibleRegion().GetSize(2) - 1);
	
  // Initialize slice iterator for traversing the input image
  typedef itk::ImageSliceConstIteratorWithIndex<ImageType> SliceIteratorType;
  SliceIteratorType sliceIt(inputImage, inputImage->GetLargestPossibleRegion());
  sliceIt.GoToBegin();
  sliceIt.SetFirstDirection(directions[0]); // 0 = x, 1 = y, 2 = z
  sliceIt.SetSecondDirection(directions[1]); // 0 = x, 1 = y, 2 = z

  unsigned int sliceNum = 0;
  typename ImageType::IndexType index;
  PixelType value = 0;
  
  // Iterate input image and copy each voxel, into the appropriate output image
  while(!sliceIt.IsAtEnd())
    {
	while(!sliceIt.IsAtEndOfSlice())
	  {
	  while(!sliceIt.IsAtEndOfLine())
	    {
        index = sliceIt.GetIndex();
		value = sliceIt.Value();
		
		if(sliceNum <= m_BisectIndex)
		  outputImage1->SetPixel(index, value);
		else
		  outputImage2->SetPixel(index, value);
		  		
		++sliceIt;
		}
		sliceIt.NextLine();
	  }
	  
	sliceIt.NextSlice();
  	++sliceNum;
  	}

}

} // end namespace itk


#endif