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

Program:   Insight Segmentation & Registration Toolkit
Module:    $RCSfile: itkEllipseBoundaryToImageFilter.txx,v $
Language:  C++
Date:      $Date: 2005/08/29 20:06:17 $
Version:   $Revision: 1.20 $

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 _itkEllipseBoundaryToImageFilter_txx
#define _itkEllipseBoundaryToImageFilter_txx

#include <itkEllipseBoundaryToImageFilter.h>

namespace itk
{

/** Constructor */
template < class TInputSpatialObject , class TOutputImage >
EllipseBoundaryToImageFilter<TInputSpatialObject, TOutputImage>
::EllipseBoundaryToImageFilter():Superclass ( ) 
  {
  m_PointSetFilter = PointSetFilterType::New();
  }

/** Destructor */
template < class TInputSpatialObject , class TOutputImage >
EllipseBoundaryToImageFilter<TInputSpatialObject , TOutputImage>
::~EllipseBoundaryToImageFilter()
{
}

//----------------------------------------------------------------------------

/** Update */
template < class TInputSpatialObject , class TOutputImage >
void
EllipseBoundaryToImageFilter<TInputSpatialObject , TOutputImage>
::GenerateData(void)
{
  unsigned int i;
  itkDebugMacro(<< "EllipseBounaryToImageFilter::Update() called");

  // Get the input and output pointers 
  const InputSpatialObjectType * InputObject  = this->GetInput();

  // Generate the image
  SizeType size;

  InputObject->ComputeBoundingBox();
  for(i=0;i<ObjectDimension;i++)
    {
    size[i] = (long unsigned int)(InputObject->GetBoundingBox()->GetMaximum()[i]
                                  - InputObject->GetBoundingBox()->GetMinimum()[i]);
    }
  
  // If the size of the output has been explicitly specified, the filter
  // will set the output size to the explicit size, otherwise the size from the spatial
  // object's bounding box will be used as default.

  bool specified = false;
  for (i = 0; i < OutputImageDimension; i++)
    {
    if (m_Size[i] != 0)
      {
      specified = true;
      break;
      }
    }

  if (!specified)
    {
    m_Size = size ;
    }

  // If the spacing has been explicitly specified, the filter
  // will set the output spacing to that explicit spacing, otherwise the spacing from
  // the spatial object is used as default.
  
  specified = false;
  for (i = 0; i < OutputImageDimension; i++)
    {
    if (m_Spacing[i] != 0)
      {
      specified = true;
      break;
      }
    }

  if (!specified)
    {
    for (i = 0; i < OutputImageDimension; i++)
      {
      m_Spacing [ i ] = 
        InputObject->GetIndexToObjectTransform()->GetScaleComponent() [ i ] ;   // set spacing
      }
    }

  const TransformType * constTransform = InputObject -> GetObjectToWorldTransform ( )  ; 
  const double PI = 3.1415926 ; 
  InputSpatialObjectType:: ArrayType radiuses = InputObject -> GetRadius ( ) ;//in pixel

  PointSetType::Pointer pointsSet = PointSetType::New();

  typedef PointSetType::PointType PointType;
   
  int counter = 0; 
  for ( double phi = 0 ; phi <=  PI ; phi += PI / 100 )
    {
    for( double theta = 0 ; theta < 2 * PI ; theta += PI / 100)
      {
      TransformType::InputPointType inputPoint ;
      inputPoint [ 0 ] = static_cast< long int > (
        radiuses [ 0 ] * cos ( theta ) * sin ( phi )  * m_Spacing [ 0 ] ) ;
      inputPoint [ 1 ] = static_cast< long int > (
        radiuses [ 1 ] * sin ( theta ) * sin ( phi ) * m_Spacing [ 1 ] ) ;
      inputPoint [ 2 ] = static_cast< long int > (
        radiuses [ 2 ] * cos ( phi )  * m_Spacing [ 2 ] ) ;

      TransformType::OutputPointType transformedPoint = constTransform
        -> TransformPoint ( inputPoint ) ;
      PointType point;
      for (i = 0; i < OutputImageDimension; i++)
        {
        point [ i ] = transformedPoint [ i ] ;
        }
      pointsSet->SetPoint( counter , point );
      counter ++ ;
      }
    }

  PointSetType::Pointer tmpPointSet = pointsSet;

  m_PointSetFilter->SetOrigin( m_Origin );
  m_PointSetFilter->SetSize( m_Size );
  m_PointSetFilter->SetSpacing( m_Spacing );
  m_PointSetFilter->SetInsideValue ( m_InsideValue );
  m_PointSetFilter->SetOutsideValue ( m_OutsideValue );
  m_PointSetFilter->SetInput( tmpPointSet );
  try
    {
    m_PointSetFilter->Update();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Exception Caught !" << std::endl;
    std::cerr << excep << std::endl;
    }


  this -> GraftOutput ( m_PointSetFilter -> GetOutput ( ) ) ;
  itkDebugMacro(<< "EllispeBoundaryToImageFilter::Update() finished");

} // end update function  


} // end namespace itk

#endif