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

  Program:   ShapeDetectionLevelSetImageFilter + GeodesicActiveContourImageFilter
  Language:  C++
  Author:	 Amit Mukherjee
  Date:      $Date: 2007/10/10
  Copyright (c) Amit Mukherjee 2007,
=========================================================================
  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.

=========================================================================*/
//ShapeDetectionLevelSetImageFilter + GeodesicActiveContourImageFilter
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
//  Software Guide : BeginCommandLineArgs
#ifdef __BORLANDC__
#define ITK_LEAN_AND_MEAN
#endif


#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"

#include "itkImage.h"
#include "itkFastMarchingImageFilter.h"
#include "itkShapeDetectionLevelSetImageFilter.h"
#include "itkGeodesicActiveContourLevelSetImageFilter.h"

#include "itkBinaryThresholdImageFilter.h"

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkRescaleIntensityImageFilter.h"
#include "itkGeodesicActiveContourLevelSetImageFilter.h"

#define N_SEEDS 21

int main( int argc, char *argv[] )	{
	
	
	if( argc < 6 )    {
			std::cerr << "Missing Parameters " << std::endl;
	    	std::cerr << "Usage: " << argv[0];
	    	std::cerr << " inputImage  outputImage";
	    	std::cerr << " switch [0-ShapeDetectioLevelsetFilter 1-GeodesicActiveContourLevelsetFilter]";
	    	std::cerr << " curvatureScaling propagationScaling" << std::endl;
	    	return 1;
    }
    
	long int seed_points[N_SEEDS][3] = {
										{	123,  166, 45	},
										{	66,  180,  45	},
										{	62,  165,   45	},
										{	62,  166,  45	},
										{	98,  112,  45	},
										{	92,  109,  45	},
										{	68,  183,  50	},
										{	98,  115,  50	},
										{	90,  114,  50	},
										{	109, 96,   50	},
										{	79,  95,   50	},
										{	99,  109,  56	},
										{	99,  109,  56	},
										{	90,  110,  56	},		
										{	70,  156,  57	},
										{	113,  158,  57	},
										{	115,  160,  60	},
										{	115,  161,  60	},
										{	70,  156,  60	},
										{	109, 108,  63	},
										{	112,  159,  63	}
									};
	
	  	
    
	typedef   float           InternalPixelType;
	const     unsigned int    Dimension = 3;
	typedef itk::Image< InternalPixelType, Dimension >  InternalImageType;
	typedef unsigned char OutputPixelType;
	typedef itk::Image< OutputPixelType, Dimension > OutputImageType;


	typedef  itk::ImageFileReader< InternalImageType > ReaderType;
	typedef  itk::ImageFileWriter<  OutputImageType  > WriterType;
	typedef  itk::ImageFileWriter<  InternalImageType  > WriterInternalImageType;

	ReaderType::Pointer reader = ReaderType::New();
	WriterType::Pointer writer = WriterType::New();
	WriterInternalImageType:: Pointer wr0 = WriterInternalImageType::New();
	WriterInternalImageType:: Pointer wr1 = WriterInternalImageType::New();
	WriterInternalImageType:: Pointer wr2 = WriterInternalImageType::New();
	WriterInternalImageType:: Pointer wr3 = WriterInternalImageType::New();


	reader->SetFileName( argv[1] );
	writer->SetFileName( argv[2] );

    
    typedef itk::BinaryThresholdImageFilter< InternalImageType, OutputImageType >    ThresholdingFilterType;
    ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
    thresholder->SetLowerThreshold( -1000.0 );
	thresholder->SetUpperThreshold(     0.0 );
	thresholder->SetOutsideValue(  0  );
	thresholder->SetInsideValue(  255 );

    //Smoothing filter settings    
    typedef   itk::CurvatureAnisotropicDiffusionImageFilter< InternalImageType, InternalImageType >  SmoothingFilterType;
    SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
    smoothing->SetTimeStep( 0.06 );
	smoothing->SetNumberOfIterations(  5 );
  	smoothing->SetConductanceParameter( 9.0 );
    
    //Gradient magnitude filter settings
    typedef   itk::GradientMagnitudeRecursiveGaussianImageFilter< InternalImageType, InternalImageType >  GradientFilterType;
	GradientFilterType::Pointer  gradientMagnitude = GradientFilterType::New();
	const double sigma = 1.0;
	gradientMagnitude->SetSigma(  sigma  );
	
	//Sigmoid filter settings
	typedef   itk::SigmoidImageFilter< InternalImageType, InternalImageType >  SigmoidFilterType;
	SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
	sigmoid->SetOutputMinimum(  0.0  );
  	sigmoid->SetOutputMaximum(  1.0  );
    const double alpha =  4;
  	const double beta  =  400;	
  	sigmoid->SetAlpha( alpha );
  	sigmoid->SetBeta(  beta  );
  	
  	// Fast marching filter settings
  	typedef  itk::FastMarchingImageFilter< InternalImageType, InternalImageType >   FastMarchingFilterType;
	FastMarchingFilterType::Pointer  fastMarching = FastMarchingFilterType::New();
	
	typedef FastMarchingFilterType::NodeContainer           NodeContainer;
  	typedef FastMarchingFilterType::NodeType                NodeType;
  	NodeContainer::Pointer seeds = NodeContainer::New();
 	
 	
 	InternalImageType::IndexType  seedPosition, seedPosition2;	
 	const double initialDistance = 3.0;
	NodeType node, node1, node2;
	const double seedValue = -1*initialDistance;
 	
 	for (int i = 0; i<N_SEEDS; i++)	{
		seedPosition[0] = seed_points[i][0];
		seedPosition[1] = seed_points[i][1];
		seedPosition[2] = seed_points[i][2];
		node.SetValue( seedValue );
		node.SetIndex( seedPosition );
		seeds->InsertElement( i, node );
		std::cout << "Seed Position " << seedPosition << std::endl;
	}
 	
 	
  	
 	fastMarching->SetTrialPoints(  seeds  );
 	fastMarching->SetSpeedConstant( 1.0 );
 	reader->Update();
 	fastMarching->SetOutputSize( reader->GetOutput()->GetBufferedRegion().GetSize() );
 	std::cout << "Image Size  :" << reader->GetOutput()->GetBufferedRegion().GetSize() << std::endl;
 	
// 	const double curvatureScaling   = 1;
// 	const double propagationScaling = 50;
 	const double curvatureScaling   = atof(argv[4]); 
 	const double propagationScaling = atof(argv[5]);
 	
 	typedef  itk::ShapeDetectionLevelSetImageFilter< InternalImageType, InternalImageType >    ShapeDetectionFilterType;
 	ShapeDetectionFilterType::Pointer   shapeDetection = ShapeDetectionFilterType::New();
 	
 	typedef  itk::GeodesicActiveContourLevelSetImageFilter< InternalImageType, InternalImageType >    GeodesicActiveContourFilterType;
	GeodesicActiveContourFilterType::Pointer geodesicActiveContour = GeodesicActiveContourFilterType::New();
		
 	if (atoi(argv[3])== 0)	{
		//Shape detection Levelset filter parameters
		shapeDetection->SetPropagationScaling(  propagationScaling );
  		shapeDetection->SetCurvatureScaling( curvatureScaling );
		shapeDetection->SetMaximumRMSError( 0.002 );
  		shapeDetection->SetNumberOfIterations( 800 );
	}
    else if (atoi(argv[3])== 1)	{	
		//Geodesic active contour level set filter setting
		geodesicActiveContour->SetPropagationScaling( propagationScaling );
		geodesicActiveContour->SetCurvatureScaling( curvatureScaling );
		geodesicActiveContour->SetAdvectionScaling( 1.0 );
		geodesicActiveContour->SetMaximumRMSError( 0.002 );
		geodesicActiveContour->SetNumberOfIterations( 800 );
	}
 	
	//Pipeline
	smoothing->SetInput( reader->GetOutput() );
	smoothing->Update(); 
	
	gradientMagnitude->SetInput( smoothing->GetOutput() );
	gradientMagnitude->Update();
	std::cout << "Gradient Computation done .." << std::endl;
	
	sigmoid->SetInput( gradientMagnitude->GetOutput() );
	sigmoid->Update();
	std::cout << "Sigmoid Computation done .." << std::endl;
	
	if (atoi(argv[3])== 0)	{
		shapeDetection->SetInput( fastMarching->GetOutput() );
		shapeDetection->SetFeatureImage( sigmoid->GetOutput() );
		shapeDetection->Update();
		std::cout << "Fast Marching done .." << std::endl;
	
		thresholder->SetInput( shapeDetection->GetOutput() );
		thresholder->Update();
		std::cout << "Shape Detection Levelset Computation done .." << std::endl;
	}
	else if (atoi(argv[3])== 1)	{
		geodesicActiveContour->SetInput(  fastMarching->GetOutput() );
		geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );
		geodesicActiveContour->Update();
		std::cout << "Fast Marching done .." << std::endl;
	
		thresholder->SetInput( geodesicActiveContour->GetOutput() );
		geodesicActiveContour->Update();
		std::cout << "Geodesic Active Contour Computation done .." << std::endl;
	}

	writer->SetInput( thresholder->GetOutput() );
	
	try		{
		writer->Update();
	}
		
	catch( itk::ExceptionObject & excep )    {
	    std::cerr << "Exception caught !" << std::endl;
	    std::cerr << excep << std::endl;
    }
	
	std::cout << "ALL Computation done .." << std::endl;
	
	
	wr0->SetInput( gradientMagnitude->GetOutput());
	wr0->SetFileName("GradientMagnitude.mha");
	wr0->Update();
	
	
	wr1->SetInput( sigmoid->GetOutput());
	wr1->SetFileName("SigmoidFunction.mha");
    wr1->Update();
	
	
	wr2->SetInput( fastMarching->GetOutput());
	wr2->SetFileName("FastMarchignOpt.mha");
	wr2->Update();
  	
	
  	if (atoi(argv[3])== 0)	{
		wr3->SetInput( shapeDetection->GetOutput());
		wr3->SetFileName("ShapeDetection.mha");
		wr3->Update();
		std::cout << std::endl;
		std::cout << "Max. no. iterations: " << shapeDetection->GetNumberOfIterations() << std::endl;
		std::cout << "Max. RMS error: " << shapeDetection->GetMaximumRMSError() << std::endl;
		std::cout << std::endl;
		std::cout << "No. elpased iterations: " << shapeDetection->GetElapsedIterations() << std::endl;
  		std::cout << "RMS change: " << shapeDetection->GetRMSChange() << std::endl;
	}

	else if (atoi(argv[3])== 1)	{
		wr3->SetInput( geodesicActiveContour->GetOutput());
		wr3->SetFileName("GeodesicActiveContour.mha");
		wr3->Update();
		std::cout << std::endl;
		std::cout << "Max. no. iterations: " << geodesicActiveContour->GetNumberOfIterations() << std::endl;
		std::cout << "Max. RMS error: " << geodesicActiveContour->GetMaximumRMSError() << std::endl;
		std::cout << std::endl;
		std::cout << "No. elpased iterations: " << geodesicActiveContour->GetElapsedIterations() << std::endl;
		std::cout << "RMS change: " << geodesicActiveContour->GetRMSChange() << std::endl;
	}
	return 0;
}