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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkAnytimeImageToImageMetric.txx,v $
  Language:  C++
  Date:      $Date: 2007/09/11 18:12:58 $
  Version:   $Revision: 1.18 $

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

#include "itkAnytimeImageToImageMetric.h"


namespace itk
{


template < class TFixedImage, class TMovingImage >
  double AnytimeImageToImageMetric<TFixedImage,TMovingImage>
  ::GetComputationLevel(int channel) const{
    if(channel>=NUMCHANNELS || channel<0) {
      itkExceptionMacro( "illegal Channel "<<channel);
    }
    return m_ComputationLevel[channel];
  }


template < class TFixedImage, class TMovingImage >
  void AnytimeImageToImageMetric<TFixedImage,TMovingImage>
  ::SetComputationLevel(double percent,int channel) const{
    if(channel>=NUMCHANNELS || channel<0) {
      itkExceptionMacro( "illegal Channel "<<channel);
    }
    if(percent>1) {
      itkWarningMacro( "Cannot do more than 100% computation");
    }
    //std::cout<<"Setting computation level to: "<<percent<<" was "<<m_ComputationLevel[channel]<<std::endl;
    m_ComputationLevel[channel]=percent;
  }



template < class TFixedImage, class TMovingImage >
  double AnytimeImageToImageMetric<TFixedImage,TMovingImage>
  ::IncrementComputationLevel(double percent,int channel) {
    if(channel>=NUMCHANNELS || channel<0) {
      itkExceptionMacro( "illegal Channel "<<channel);
    }
    if(percent>1) {
      itkExceptionMacro( "Cannot do more than 100% computation");
    }
    m_ComputationLevel[channel]+=percent;
    if(m_ComputationLevel[channel]>1.0) m_ComputationLevel[channel]=1.0;
    //std::cout<<"Setting computation level to: "<<m_ComputationLevel[channel]<<std::endl;
    return m_ComputationLevel[channel];
  }


template < class TFixedImage, class TMovingImage >
  void AnytimeImageToImageMetric<TFixedImage,TMovingImage>
  ::ResetComputation (const int channel) const {
    if(channel>=NUMCHANNELS || channel<0) {
      itkExceptionMacro( "illegal Channel "<<channel);
    }
    m_nFixedImageSamples[channel]=0;
    m_nSamples[channel]=0;
    //m_CachedParameters[channel]=NULL;
    // std::cout<<"Metric resetting channel "<<channel<<std::endl;
  }


template < class TFixedImage, class TMovingImage >
  void AnytimeImageToImageMetric<TFixedImage,TMovingImage>
  ::SetChannel (const int channel) {
    if(channel>=NUMCHANNELS || channel<0) {
      itkExceptionMacro( "illegal Channel "<<channel);
    }
    m_CurrentChannel=channel;
  }


///////////////////////////////////////////
/**
 * Get the match Measure
 */
template < class TFixedImage, class TMovingImage  >
typename AnytimeImageToImageMetric<TFixedImage,TMovingImage>::MeasureType
AnytimeImageToImageMetric<TFixedImage,TMovingImage>
::GetValue( const ParametersType& parameters ) const
{
  //this->ResetComputation(m_CurrentChannel);
  return static_cast<MeasureType>(this->GetValue(parameters,m_CurrentChannel));
}


///////////////////////////////////////////
/**
 * Get the match Measure
 */
template < class TFixedImage, class TMovingImage  >
typename AnytimeImageToImageMetric<TFixedImage,TMovingImage>::MeasureType
AnytimeImageToImageMetric<TFixedImage,TMovingImage>
::GetValue( const ParametersType& parameters, const int channel, const double level) const
{
  this->SetComputationLevel(level,channel);
  //this->ResetComputation(m_CurrentChannel);
  return static_cast<MeasureType>(this->GetValue(parameters,channel));
}


///////////////////////////////////////////
/**
 * Get the match Measure
 */
template < class TFixedImage, class TMovingImage  >
typename AnytimeImageToImageMetric<TFixedImage,TMovingImage>::MeasureType
AnytimeImageToImageMetric<TFixedImage,TMovingImage>
::GetValue( const ParametersType& parameters, const double level) const
{
  this->SetComputationLevel(level,m_CurrentChannel);
  //this->ResetComputation(m_CurrentChannel);
  return static_cast<MeasureType>(this->GetValue(parameters,m_CurrentChannel));
}



/** Get the derivatives of the match measure. */
template < class TFixedImage, class TMovingImage  >
void AnytimeImageToImageMetric<TFixedImage,TMovingImage>
::GetDerivative( const TransformParametersType & parameters,
		    DerivativeType & derivative, int channel,
		  double level ) const {
  this->SetComputationLevel(level,channel);
  this->GetDerivative(parameters,derivative,channel);
}
/**  Get value and derivatives for multiple valued optimizers. */
template < class TFixedImage, class TMovingImage  >
void AnytimeImageToImageMetric<TFixedImage,TMovingImage>
::GetValueAndDerivative( const TransformParametersType & parameters,
			    MeasureType& Value, DerivativeType& Derivative,
			    int channel,
			 double level ) const {
  this->SetComputationLevel(level,channel);
  this->GetValueAndDerivative(parameters,Value,Derivative,channel);
}
 
/** Get the derivatives of the match measure. */
template < class TFixedImage, class TMovingImage  >
void AnytimeImageToImageMetric<TFixedImage,TMovingImage>
::GetDerivative( const TransformParametersType & parameters,
		 DerivativeType & derivative, double level ) const {
  this->SetComputationLevel(level,m_CurrentChannel);
  this->GetDerivative(parameters,derivative,m_CurrentChannel);
}
/**  Get value and derivatives for multiple valued optimizers. */
template < class TFixedImage, class TMovingImage  >
void AnytimeImageToImageMetric<TFixedImage,TMovingImage>
::GetValueAndDerivative( const TransformParametersType & parameters,
			    MeasureType& Value, DerivativeType& Derivative,
			     double level ) const{
  this->SetComputationLevel(level,m_CurrentChannel);
  this->GetValueAndDerivative(parameters,Value,Derivative,m_CurrentChannel);
}
 
template < class TFixedImage, class TMovingImage  >
void AnytimeImageToImageMetric<TFixedImage,TMovingImage>
::GetDerivative( const TransformParametersType & parameters,
		 DerivativeType & derivative) const {
  this->GetDerivative(parameters,derivative,m_CurrentChannel);
}
/**  Get value and derivatives for multiple valued optimizers. */
template < class TFixedImage, class TMovingImage  >
void AnytimeImageToImageMetric<TFixedImage,TMovingImage>
::GetValueAndDerivative( const TransformParametersType & parameters,
			 MeasureType& Value, 
			 DerivativeType& Derivative) const {
  this->GetValueAndDerivative(parameters,Value,Derivative,m_CurrentChannel);
}



/**
 * Initialize
 */
template <class TFixedImage, class TMovingImage> 
void
AnytimeImageToImageMetric<TFixedImage,TMovingImage>
::Initialize(void) throw ( ExceptionObject )
{

  this->Superclass::Initialize();
  // Cache the number of transformation parameters
  this->m_NumberOfParameters = this->m_Transform->GetNumberOfParameters();
  m_CurrentChannel=0;
  //m_TotalPixelsProcessed=0;
  for (int i=0;i<NUMCHANNELS;i++) {
    m_nFixedImageSamples[i]=0;
    m_nSamples[i]=0;
    m_CachedParameters[i]=new ParametersType(m_NumberOfParameters);
    m_DerivativeCache[i]=new DerivativeType(m_NumberOfParameters);
    m_ComputationLevel[i]=1.0;
  }
  /**
   * Allocate memory for the fixed image sample container.
   */
  std::cout<<"Initializing space for "<< m_NumberOfSpatialSamples<< " samples...."<<std::endl;
  std::cout<<"Samples are "<< sizeof(FixedImageSpatialSample) << " bytes. total is "<<sizeof(FixedImageSpatialSample)*m_NumberOfSpatialSamples<<std::endl;
  try {
  m_FixedImageSamples.resize( m_NumberOfSpatialSamples);
  } catch (std::exception e) {
    itkExceptionMacro( "Resize for " << m_NumberOfSpatialSamples<<" samples failed: ");
  }
  //std::cout<<"done"<<std::endl;
  /** 
   * Uniformly sample the fixed image (within the fixed image region)
   * to create the sample points list.
   */
  this->SampleFixedImageDomain( m_FixedImageSamples );
}



/**
 * Uniformly sample the fixed image domain using a random walk
 */
template < class TFixedImage, class TMovingImage >
void
AnytimeImageToImageMetric<TFixedImage,TMovingImage>
::SampleFixedImageDomain( FixedImageSpatialSampleContainer& samples )
{

  // Set up a random interator within the user specified fixed image region.
  RandomIterator randIter( this->m_FixedImage, this->GetFixedImageRegion() );

  randIter.SetNumberOfSamples( this->GetFixedImageRegion().GetNumberOfPixels() );
  randIter.ReinitializeSeed(m_RandomSeed);
  if(m_PriorityImage!=NULL)    randIter.SetPriorityImage(m_PriorityImage);
  randIter.GoToBegin();
  m_NumberOfAvailableSamples=0;
  typename FixedImageSpatialSampleContainer::iterator iter;
  typename FixedImageSpatialSampleContainer::const_iterator end=samples.end();
  //std::cout<<"Trying to get "<<m_NumberOfSpatialSamples<<" samples"<<std::endl;

  if( this->m_FixedImageMask )
    {

    typename Superclass::InputPointType inputPoint;

    iter=samples.begin();
    while( iter != end && !randIter.IsAtEnd())
      {
      // Get sampled index
      FixedImageIndexType index = randIter.GetIndex();
      // Check if the Index is inside the mask, translate index to point
      this->m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );

      // If not inside the mask, ignore the point
      if( !this->m_FixedImageMask->IsInside( inputPoint ) )
        {
        ++randIter; // jump to another random position
        continue;
        }

      // Save index
      (*iter).FixedImageIndex = index;
      // Get sampled fixed image value
      (*iter).FixedImageValue = randIter.Get();
      // Translate index to point
      (*iter).FixedImagePointValue = inputPoint;

      // Jump to random position
      ++randIter;
      ++iter;
      ++m_NumberOfAvailableSamples;
      }
    }
  else
    {
    for( iter=samples.begin(); iter != end && !randIter.IsAtEnd(); ++iter,      ++m_NumberOfSpatialSamples )
      {
      // Get sampled index
      FixedImageIndexType index = randIter.GetIndex();
      // Save index
      (*iter).FixedImageIndex = index;
      // Get sampled fixed image value
      (*iter).FixedImageValue = randIter.Get();
      // Translate index to point
      this->m_FixedImage->TransformIndexToPhysicalPoint( index,
                                                   (*iter).FixedImagePointValue );
      // Jump to random position
      ++randIter;
      ++m_NumberOfAvailableSamples;

      }
    }
  m_FixedImageSamples.resize( m_NumberOfAvailableSamples);
  // m_NumberOfSpatialSamples=samples.size();
  //std::cout<<"Size of sample set"<<m_FixedImageSamples.size()<<std::endl;
}

/**
 * write out an image of the pixels used
 */
template < class TFixedImage, class TMovingImage >
void
AnytimeImageToImageMetric<TFixedImage,TMovingImage>
::WritePixelsUsed (const char* filename ){
	if(!this->m_FixedImage) {
    itkExceptionMacro( << "Fixed image has not been assigned" );
	}
  typename OutputPixelImageType::Pointer pixelImage = OutputPixelImageType::New();
  typename TFixedImage::RegionType region=this->m_FixedImage->GetLargestPossibleRegion();
  pixelImage->SetLargestPossibleRegion( region );
  pixelImage->SetBufferedRegion( region );
  pixelImage->SetRequestedRegion( region );
  pixelImage->Allocate();
  pixelImage->FillBuffer(0);

  pixelImage->SetSpacing(this->m_FixedImage->GetSpacing());
  pixelImage->SetOrigin(this->m_FixedImage->GetOrigin());
  pixelImage->SetDirection(this->m_FixedImage->GetDirection());
  typename FixedImageSpatialSampleContainer::const_iterator fiter=this->m_FixedImageSamples.begin();
  typename FixedImageSpatialSampleContainer::const_iterator fend = 
    this->m_FixedImageSamples.end();
  while (fiter!=fend) {
	  pixelImage->SetPixel((*fiter).FixedImageIndex,255);
	  ++fiter;
  }
  typename PixelImageWriterType::Pointer writer=PixelImageWriterType::New();
  writer->SetInput(pixelImage);
  writer->SetFileName(filename);
  writer->Update();


}
} // end namespace itk

#endif