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

  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.

=========================================================================*/
#ifndef __itkLesionSegmentationImageFilter7_txx
#define __itkLesionSegmentationImageFilter7_txx
#include "itkLesionSegmentationImageFilter7.h"

#include "itkNumericTraits.h"
#include "itkProgressReporter.h"
#include "itkGradientMagnitudeImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

namespace itk
{
  
template <class TInputImage, class TOutputImage>
LesionSegmentationImageFilter7<TInputImage, TOutputImage>::
LesionSegmentationImageFilter7()
{
  m_GradientMagnitudeSigmoidFeatureGenerator = GradientMagnitudeSigmoidGeneratorType::New();
  m_LesionSegmentationMethod = LesionSegmentationMethodType::New();
  m_LungWallFeatureGenerator = LungWallGeneratorType::New();
  m_VesselnessFeatureGenerator = VesselnessGeneratorType::New();
  m_SigmoidFeatureGenerator = SigmoidFeatureGeneratorType::New();
  m_FeatureAggregator = FeatureAggregatorType::New();
  m_SegmentationModule = SegmentationModuleType::New();
  m_CropFilter = CropFilterType::New();
  m_InputSpatialObject = InputImageSpatialObjectType::New();

  // Report progress.
  m_CommandObserver    = CommandType::New();
  m_CommandObserver->SetCallbackFunction(
    this, &Self::ProgressUpdate );
  m_LungWallFeatureGenerator->AddObserver( 
      itk::ProgressEvent(), m_CommandObserver );
  m_SigmoidFeatureGenerator->AddObserver( 
      itk::ProgressEvent(), m_CommandObserver );
  m_VesselnessFeatureGenerator->AddObserver( 
      itk::ProgressEvent(), m_CommandObserver );
  m_GradientMagnitudeSigmoidFeatureGenerator->AddObserver( 
      itk::ProgressEvent(), m_CommandObserver );
  m_SegmentationModule->AddObserver( 
      itk::ProgressEvent(), m_CommandObserver );
  m_CropFilter->AddObserver( 
      itk::ProgressEvent(), m_CommandObserver );

  // Connect pipeline
  m_LungWallFeatureGenerator->SetInput( m_InputSpatialObject );
  m_SigmoidFeatureGenerator->SetInput( m_InputSpatialObject );
  m_VesselnessFeatureGenerator->SetInput( m_InputSpatialObject );
  m_GradientMagnitudeSigmoidFeatureGenerator->SetInput( m_InputSpatialObject );
  m_FeatureAggregator->AddFeatureGenerator( m_LungWallFeatureGenerator );
  m_FeatureAggregator->AddFeatureGenerator( m_VesselnessFeatureGenerator );
  m_FeatureAggregator->AddFeatureGenerator( m_SigmoidFeatureGenerator );
  m_FeatureAggregator->AddFeatureGenerator( m_GradientMagnitudeSigmoidFeatureGenerator );
  m_LesionSegmentationMethod->AddFeatureGenerator( m_FeatureAggregator );
  m_LesionSegmentationMethod->SetSegmentationModule( m_SegmentationModule );

  // Populate some parameters
  m_LungWallFeatureGenerator->SetLungThreshold( -400 );
  m_VesselnessFeatureGenerator->SetSigma( 1.0 );
  m_VesselnessFeatureGenerator->SetAlpha1( 0.5 );
  m_VesselnessFeatureGenerator->SetAlpha2( 2.0 );
  m_SigmoidFeatureGenerator->SetAlpha( 1.0 );
  m_GradientMagnitudeSigmoidFeatureGenerator->SetSigma(1.0);
  m_GradientMagnitudeSigmoidFeatureGenerator->SetAlpha( -0.1 );
  m_GradientMagnitudeSigmoidFeatureGenerator->SetBeta( 150.0 );
  m_FastMarchingStoppingTime = 5.0;
  m_FastMarchingDistanceFromSeeds = 2.0;  
  m_SigmoidBeta = -200.0;
  m_StatusMessage = "";
  m_SegmentationModule->SetCurvatureScaling(1.0);
  m_SegmentationModule->SetAdvectionScaling(50.0);
  m_SegmentationModule->SetPropagationScaling(10.0);
  m_SegmentationModule->SetMaximumRMSError(0.0002);
  m_SegmentationModule->SetMaximumNumberOfIterations(300);
}
 
template <class TInputImage, class TOutputImage>
void 
LesionSegmentationImageFilter7<TInputImage,TOutputImage>
::GenerateInputRequestedRegion() throw(InvalidRequestedRegionError)
{
  // call the superclass' implementation of this method
  Superclass::GenerateInputRequestedRegion();
}

template <class TInputImage, class TOutputImage>
void 
LesionSegmentationImageFilter7<TInputImage,TOutputImage>
::GenerateOutputInformation() 
{
  // get pointers to the input and output
  typename Superclass::OutputImagePointer      outputPtr = this->GetOutput();
  typename Superclass::InputImageConstPointer  inputPtr  = this->GetInput();

  if ( !outputPtr || !inputPtr)
    {
    return;
    }

  // Set the output image size to the same value as the region of interest.
  RegionType region;
  IndexType  start;
  start.Fill(0);

  region.SetSize( m_RegionOfInterest.GetSize() );
  region.SetIndex( start );
 
  // Copy Information without modification.
  outputPtr->CopyInformation( inputPtr );

  // Adjust output region
  outputPtr->SetLargestPossibleRegion(region);

  // Correct origin of the extracted region.
  IndexType roiStart( m_RegionOfInterest.GetIndex() );
  typename Superclass::OutputImageType::PointType  outputOrigin;
  typedef Image< ITK_TYPENAME TInputImage::PixelType,
    Superclass::InputImageDimension > ImageType;
  typename ImageType::ConstPointer imagePtr =
    dynamic_cast< const ImageType * >( inputPtr.GetPointer() );
  if ( imagePtr )
    {
    // Input image supports TransformIndexToContinuousPoint
    inputPtr->TransformIndexToPhysicalPoint( roiStart, outputOrigin);
    }
  else
    {
    // Generic type of image
    const typename Superclass::InputImageType::PointType&
      inputOrigin = inputPtr->GetOrigin();

    const typename Superclass::InputImageType::SpacingType&
      spacing = inputPtr->GetSpacing() ;
 
    for( unsigned int i=0; i<ImageDimension; i++)
      {
      outputOrigin[i] = inputOrigin[i] + roiStart[i] * spacing[i];
      }
    }
  
  outputPtr->SetOrigin( outputOrigin );  
}

  
template< class TInputImage, class TOutputImage >
void
LesionSegmentationImageFilter7< TInputImage, TOutputImage >
::GenerateData()
{
  m_SigmoidFeatureGenerator->SetBeta( m_SigmoidBeta );
  m_SegmentationModule->SetDistanceFromSeeds(m_FastMarchingDistanceFromSeeds);
  m_SegmentationModule->SetStoppingValue(m_FastMarchingStoppingTime);

  // Allocate the output
  this->GetOutput()->SetBufferedRegion( this->GetOutput()->GetRequestedRegion() );
  this->GetOutput()->Allocate();
 
  /* DEBUG REMOVEME
  typedef itk::ImageFileReader< OutputImageType > OutputImageReaderType;
  OutputImageReaderType::Pointer outputReader = OutputImageReaderType::New();
  outputReader->SetFileName("f:/Kitware/VolView/bin/Dashboard-Static/bin/Release/fastmarchingoutput.mha");
  outputReader->Update();
  this->GraftOutput(outputReader->GetOutput());
  return;*/
  // END DEBUG REMOVEME
  
  typename  InputImageType::ConstPointer  input  = this->GetInput();

  // Crop
  m_CropFilter->SetInput(input);
  m_CropFilter->SetRegionOfInterest(m_RegionOfInterest);
  m_CropFilter->Update(); 
 
  //this->GraftOutput(m_CropFilter->GetOutput());
  //return;

  // Convert the output of cropping to a spatial object that can be fed into
  // the lesion segmentation method
  typename InputImageType::Pointer inputImage = m_CropFilter->GetOutput();
  inputImage->DisconnectPipeline();
  m_InputSpatialObject->SetImage(inputImage);

  // Seeds
  typename SeedSpatialObjectType::Pointer seedSpatialObject =
    SeedSpatialObjectType::New();
  seedSpatialObject->SetPoints(m_Seeds);
  m_LesionSegmentationMethod->SetInitialSegmentation(seedSpatialObject);
  
  // Do the actual segmentation.
  m_LesionSegmentationMethod->Update();
  
  // Graft the output.
  typename SpatialObjectType::Pointer segmentation = 
    const_cast< SpatialObjectType * >(m_SegmentationModule->GetOutput());
  typename OutputSpatialObjectType::Pointer outputObject = 
    dynamic_cast< OutputSpatialObjectType * >( segmentation.GetPointer() );
  typename OutputImageType::Pointer outputImage = 
    const_cast< OutputImageType * >(outputObject->GetImage());
  outputImage->DisconnectPipeline();
  this->GraftOutput(outputImage);

  /*typedef itk::ImageFileWriter< OutputImageType > WriterType;
  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName("output.mha");
  writer->SetInput(outputImage);
  writer->Write();*/
}


template <class TInputImage, class TOutputImage>
void LesionSegmentationImageFilter7< TInputImage,TOutputImage >
::ProgressUpdate( Object * caller, 
                  const EventObject & e )
{
  if( typeid( itk::ProgressEvent ) == typeid( e ) )
    {  
    if (dynamic_cast< CropFilterType * >(caller))
      {
      this->m_StatusMessage = "Cropping data..";
      this->UpdateProgress( m_CropFilter->GetProgress() );
      }

    else if (dynamic_cast< LungWallGeneratorType * >(caller))
      {
      // Given its iterative nature.. a cranky heuristic here.
      this->m_StatusMessage = "Generating lung wall feature..";
      this->UpdateProgress( ((double)(((int)(
        m_LungWallFeatureGenerator->GetProgress()*500))%100))/100.0 );
      }

    else if (dynamic_cast< SigmoidFeatureGeneratorType * >(caller))
      {
      this->m_StatusMessage = "Generating intensity feature..";
      this->UpdateProgress( m_SigmoidFeatureGenerator->GetProgress() );
      }

    else if (dynamic_cast< GradientMagnitudeSigmoidGeneratorType * >(caller))
      {
      m_StatusMessage = "Generating edge feature..";
      this->UpdateProgress( m_GradientMagnitudeSigmoidFeatureGenerator->GetProgress());
      }

    else if (dynamic_cast< VesselnessGeneratorType * >(caller))
      {
      m_StatusMessage = "Generating vesselness feature (Sato et al.)..";
      this->UpdateProgress( m_LungWallFeatureGenerator->GetProgress() );
      }

    else if (dynamic_cast< SegmentationModuleType * >(caller))
      {
      m_StatusMessage = "Segmenting using geodesic active contour level set..";
      this->UpdateProgress( m_SegmentationModule->GetProgress() );
      }
    }
}

template <class TInputImage, class TOutputImage>
void 
LesionSegmentationImageFilter7<TInputImage,TOutputImage>
::PrintSelf(std::ostream& os, Indent indent) const
{
  Superclass::PrintSelf(os,indent);
}

}//end of itk namespace
#endif