//--------------------------------------------------------------------------------------------
// INCLUDES
//--------------------------------------------------------------------------------------------
#include "itkFM3D.h"


#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkSigmoidImageFilter.h"
#include "itkFastMarchingImageFilter.h"
#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkVTKImageToImageFilter.h"
#include "itkImageToVTKImageFilter.h"
#include "itkVTKImageIO.h"

#include "vtkOtsuImageData.h"
#include "vtkImageCast.h"
#include "vtkMetaImageWriter.h"

#include "vtkDataSetWriter.h"


#include <wx/arrstr.h>
#include <wx/process.h>


itkFM3D::itkFM3D(){
	usefilter	= 1; // Por defecto hacer el filtrado de ruido con la curvatura anisotropica.
	alpha		= 0.5;
	beta		= 128.0;
	stopTime	= 10; // Valor por defecto;
	x = 0;
	y = 0;
	z = 0;
}
	
itkFM3D::~itkFM3D(){

}
	
void itkFM3D::AddSeed(int x, int y, int z){
	this->x = x;
	this->y = y;
	this->z = z;
}

void itkFM3D::SetAlpha(double alpha){

	this->alpha = alpha;
}

void itkFM3D::SetBeta(double beta){
	this->beta = beta;
}


void itkFM3D::SetStopTime(double stopTime){
	this->stopTime = stopTime;
}

void itkFM3D::CurvatureAnisotropicFiltertOn(){
	usefilter = 1;
}

void itkFM3D::CurvatureAnisotropicFilterOff(){
	usefilter = 0;
}


vtkImageData* itkFM3D::segment(vtkImageData *volume){

	volume->SetSpacing(1,1,1);


	double spc[3];
	volume->GetSpacing(spc);
	

	double xx = x; // * spc[0];
	double yy = y; // * spc[1];
	double zz = z; // * spc[2];


	/*FILE *ff;
	ff=fopen("programa/Ups.bat","w");
	fprintf(ff,"vtkViewer imagen.vtk %d %d %d %f %f %f \n",(int)xx,(int)yy,(int)zz,(float)alpha,(float)beta,(float)stopTime);
	fclose(ff);
    int code = wxExecute("programa/Ups.bat", wxEXEC_SYNC);*/

	




  //DEFINICION IMAGEN ENTRADA
  //--------------------------------------------------------------------------------------------
  typedef   double           InternalPixelType;
  const     unsigned int    DIMENSION_IMAGEN = 3;
  typedef itk::Image< InternalPixelType, DIMENSION_IMAGEN >  InternalImageType;
  
 


  // CONVERTIR LA IMAGEN DE ENTRADA A ITK PARA SER PROCESADA
  //--------------------------------------------------------------------------------------------

  vtkImageCast *ic = vtkImageCast::New();
  ic->SetInput(volume);
  ic->SetOutputScalarTypeToDouble();
  ic->Update();
  typedef itk::VTKImageToImageFilter<InternalImageType> VtkToItkFilterType;
  VtkToItkFilterType::Pointer toItk = VtkToItkFilterType::New();
  toItk->SetInput( ic->GetOutput() );
  toItk->Update();
   const InternalImageType* itkImageData = toItk->GetOutput();
  itkImageData->Register();

  //typedef itk::ImageFileReader<InternalImageType> ReaderType;
  //ReaderType::Pointer reader = ReaderType::New();
  //reader->SetFileName("imagen.vtk");

  //reader->Update();

  //InternalImageType* itkImageData = null; //reader->GetOutput();
	 
  


  //DEFINICION WRITER 
  //--------------------------------------------------------------------------------------------
 // typedef itk::ImageFileWriter<InternalImageType> WriterType;


 
  //DEFINICION DEL FILTRADO ANISOTROPICO
  //--------------------------------------------------------------------------------------------
  typedef itk::CurvatureAnisotropicDiffusionImageFilter<
                                InternalImageType,
                                InternalImageType > SmoothingFilterType;

  SmoothingFilterType::Pointer filtroAnisotropico = SmoothingFilterType::New();

  
   
  //DEFINICION DEL FILTRO SIGMOIDE 
  //--------------------------------------------------------------------------------------------
  typedef   itk::SigmoidImageFilter<                               
                               InternalImageType, 
                               InternalImageType >  SigmoidFilterType;

  SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();

  sigmoid->SetOutputMinimum(  0.0  );
  sigmoid->SetOutputMaximum(  255  );
  sigmoid->SetAlpha(this->alpha);
  sigmoid->SetBeta(this->beta);

  

  //DECLARACION DEL FILTRO DE FAST MARCHING
  //--------------------------------------------------------------------------------------------
  typedef  itk::FastMarchingImageFilter< InternalImageType, 
                              InternalImageType >    FastMarchingFilterType;
  FastMarchingFilterType::Pointer  fastMarching = FastMarchingFilterType::New();

  //DECLARACION DEL FILTRO THRESHOLDING BINARIO
  //--------------------------------------------------------------------------------------------
  typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType > ThresholdingFilterType;
  ThresholdingFilterType::Pointer thresholdFilter = ThresholdingFilterType::New();

 
  //--------------------------------------------------------------------------------------------   
  // ESTABLECIMIENTO DE LA SEMILLA DEL FAST MARCHING
  //--------------------------------------------------------------------------------------------
  typedef FastMarchingFilterType::NodeContainer           NodeContainer;
  typedef FastMarchingFilterType::NodeType                NodeType;
  NodeContainer::Pointer seeds = NodeContainer::New();  

  InternalImageType::IndexType  seedPosition;
  seedPosition[0] = this->x; // POSICION X DE LA SEMILLA
  seedPosition[1] = this->y; // POSICION Y DE LA SEMILLA
  seedPosition[2] = this->z; // POSICION Z DE LA SEMILLA

  NodeType node;
  const double seedValue = 0.0;
  
  node.SetValue( seedValue );
  node.SetIndex( seedPosition );

  seeds->Initialize();
  seeds->InsertElement( 0, node );


  //COLOCAR VARIAS SEMILLAS
  fastMarching->SetTrialPoints(  seeds  );
  fastMarching->SetOutputSize(itkImageData->GetBufferedRegion().GetSize() );
  fastMarching->SetStoppingValue(this->stopTime);


  //DEFINICION WRITER FILTRADO ANISOTROPICO
  //WriterType::Pointer writerFiltroRuido = WriterType::New();
  //writerFiltroRuido->SetFileName("filtrado.vtk");

  
  //DEFINICION WRITER SIGMOIDE
  //--------------------------------------------------------------------------------------------
  //WriterType::Pointer writerSigmoide = WriterType::New();
  //writerSigmoide->SetFileName("sigmoid.vtk");



  //DEFINICION WRITER FAST MARCHING MAPA DE TIEMPOS
  //--------------------------------------------------------------------------------------------
  /*typedef itk::VTKImageIO ImageIOType;
  ImageIOType::Pointer vtkIO = ImageIOType::New();
  vtkIO->SetFileTypeToASCII();

  WriterType::Pointer writerMapaTiempos = WriterType::New();
  writerMapaTiempos->SetFileName("fastmarch.vtk");
  writerMapaTiempos->SetImageIO(vtkIO);*/


  //DEFINICION WRITER THRESHOLD BINARIO
  //--------------------------------------------------------------------------------------------
//  WriterType::Pointer writerThreshold = WriterType::New();
 // writerThreshold->SetFileName("binario.vtk");
  //writerThreshold->SetImageIO(vtkIO);
  

  //--------------------------------------------------------------------------------------------
  // CONSTRUCCION DEL PIPELINE DE PROCESAMIENTO
  //--------------------------------------------------------------------------------------------
 
  // READER -> FILTRO ANISOTROPICO
  //FILTRO ANISOTROPICO -> SIGMOIDE
  if (usefilter == 1){
	//filtroAnisotropico->SetInput(itkImageData);
	//sigmoid->SetInput(filtroAnisotropico->GetOutput());
	  sigmoid->SetInput(itkImageData);
  }
  else {
	  sigmoid->SetInput(itkImageData);
  }

  // SIGMOIDE -> FAST MARCHING
  fastMarching->SetInput(sigmoid->GetOutput());
  
    
  // FAST MARCHING -> WRITER
  //writerFiltroRuido->SetInput(filtroAnisotropico->GetOutput());
  //writerSigmoide->SetInput(sigmoid->GetOutput());
  //writerMapaTiempos->SetInput(fastMarching->GetOutput());
  //writerFiltroRuido->Update();

  
  //----------------------------------------------------------------------------------------
  //Ejecucion de OTSU
  //----------------------------------------------------------------------------------------

  //DEFINICION DE IMAGEN VTK FILTRADA CON ANISOTROPICO
  //--------------------------------------------------------------------------------------------
  typedef itk::ImageToVTKImageFilter<InternalImageType> ItkToVtkFilterType;
  
 
  
  /*vtkOtsuImageData* otsu = new vtkOtsuImageData();

  if(this->usefilter){
	printf("Calculando Th Otsu sobre la imagen filtrada\n");
	 ItkToVtkFilterType::Pointer toVtk = ItkToVtkFilterType::New();
	toVtk->SetInput(filtroAnisotropico->GetOutput());
	toVtk->Update();
	vtkImageData* imagenFiltrada = toVtk->GetOutput();
	estimatedOtsuThreshold = otsu->calculateOptimalThreshold(imagenFiltrada);
  }
  else{
    printf("Calculando Th Otsu sobre la imagen original\n");
	estimatedOtsuThreshold = otsu->calculateOptimalThreshold(volume);
  }
  
  printf("BETA = %.0f, OTSU = %d\n",beta,estimatedOtsuThreshold);*/

  //sigmoid->SetBeta((double)ThresholdOtsu);
  
  //----------------------------------------------------------------------------------------
  
  
  //--------------------------------------------------------------------------------------------  
  // EJECUCION DEL PIPELINE
  //--------------------------------------------------------------------------------------------
  try
  {
	  fastMarching->Update();
	  //writerFiltroRuido->Update();
	  //writerSigmoide->Update();
	  //writerMapaTiempos->Update();
      printf("Writers (ruido,sigmoide, mapatiempos) procesados \n");
  }
  catch( itk::ExceptionObject & excep )
  {
    std::cerr << "Ocurrio un problema" << std::endl;
    std::cerr << excep << std::endl;
  }

  
  //---------------------------------------------------------
  // LEYENDO MAPA DE TIEMPOS
  //---------------------------------------------------------
  // En esta parte escribimos el mapa de tiempos a un archivo plano para
  // poder estudiarlo en Excel.

    printf("Procesando el mapa de tiempos....");

  	//ACTIVAR PARA DEBUG
	//-----------------------------
	//FILE *logger = NULL;
	//	if ((logger = freopen("mapatiempos.txt","w", stdout)) == NULL)
	//		exit(-1);

	
		
	InternalImageType::Pointer mapa = fastMarching->GetOutput();  

	typedef itk::ImageRegionConstIterator< InternalImageType > ConstIteratorType;

	ConstIteratorType constIterator( mapa, mapa->GetLargestPossibleRegion() );

	
	// inicializar histogramas
	double valor = 0.0;
	int indice = 0;

	int histoTiempos[1000];
	int histoAcumulado[1000];

	for (indice = 0; indice <= 1000; indice++){
		histoTiempos[indice] = 0;
		histoAcumulado[indice] = 0;
	}

	
	// calcular histograma
	for (constIterator.GoToBegin(); !constIterator.IsAtEnd(); ++ constIterator){
		valor = constIterator.Get();
		if (valor <= 10){
			indice = (int)(valor * 100);
			histoTiempos[indice]++;

		}
	}

	//calcular histograma acumulado
	for (indice = 1 ; indice <= 1000; indice++){
		histoAcumulado[indice]= histoAcumulado[indice-1]+ histoTiempos[indice];
	}

	//buscar la frecuencia minima = derivada de la funcion de frecuencia acumulada
	int frecuenciaMinima = 1000000000;
	int indiceBuscado = -1;
	for (indice = 0; indice <= 1000; indice++){
		if (histoTiempos[indice] < frecuenciaMinima){
			frecuenciaMinima = histoTiempos[indice];
			indiceBuscado = indice;
		}
		
	}

    
	double valorThreshold = (double)indiceBuscado / 100.0;
	
	printf("VALOR DE SEGMENTACION: %.2f\n\n",valorThreshold);
	

	//for (indice = 0; indice <= 1000; indice++){
	//	printf("%.2f:%d:%d\n",(double)indice/100.0, histoTiempos[indice], histoAcumulado[indice]);
	//}

	
	//fclose(logger);
	//logger = NULL;

	printf("Mapa de tiempos OK");


	thresholdFilter->SetOutsideValue( 0 );
	thresholdFilter->SetInsideValue( 255 );
	thresholdFilter->SetLowerThreshold(0);
	thresholdFilter->SetUpperThreshold(valorThreshold);

	thresholdFilter->SetInput(fastMarching->GetOutput());
//	writerThreshold->SetInput(thresholdFilter->GetOutput());

	try{
		thresholdFilter->Update();
		//writerThreshold->Update();
	}
	catch( itk::ExceptionObject & excep )
	{
		std::cerr << "Ocurrio un problema" << std::endl;
		std::cerr << excep << std::endl;
	}

	ItkToVtkFilterType::Pointer resultado = ItkToVtkFilterType::New();
	resultado->SetInput(thresholdFilter->GetOutput());
	vtkImageData* salida = resultado->GetOutput();
	resultado->Register();
	salida->Update();
	salida->ReleaseDataFlagOff();




	return salida;
}