itkFM3D.cxx

Go to the documentation of this file.

//--------------------------------------------------------------------------------------------
// 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;
}