lib/maracasVisuLib/src/kernel/include/itkFM3D.cxx

   1 //--------------------------------------------------------------------------------------------
   2 // INCLUDES
   3 //--------------------------------------------------------------------------------------------
   4 #include "itkFM3D.h"
   5
   6
   7 #include "itkImage.h"
   8 #include "itkImageFileReader.h"
   9 #include "itkImageFileWriter.h"
  10 #include "itkSigmoidImageFilter.h"
  11 #include "itkFastMarchingImageFilter.h"
  12 #include "itkCurvatureAnisotropicDiffusionImageFilter.h"
  13 #include "itkBinaryThresholdImageFilter.h"
  14 #include "itkVTKImageToImageFilter.h"
  15 #include "itkImageToVTKImageFilter.h"
  16 #include "itkVTKImageIO.h"
  17
  18 #include "vtkOtsuImageData.h"
  19 #include "vtkImageCast.h"
  20 #include "vtkMetaImageWriter.h"
  21
  22 #include "vtkDataSetWriter.h"
  23
  24
  25 #include <wx/arrstr.h>
  26 #include <wx/process.h>
  27
  28
  29 itkFM3D::itkFM3D(){
  30         usefilter       = 1; // Por defecto hacer el filtrado de ruido con la curvatura anisotropica.
  31         alpha           = 0.5;
  32         beta            = 128.0;
  33         stopTime        = 10; // Valor por defecto;
  34         x = 0;
  35         y = 0;
  36         z = 0;
  37 }
  38
  39 itkFM3D::~itkFM3D(){
  40
  41 }
  42
  43 void itkFM3D::AddSeed(int x, int y, int z){
  44         this->x = x;
  45         this->y = y;
  46         this->z = z;
  47 }
  48
  49 void itkFM3D::SetAlpha(double alpha){
  50
  51         this->alpha = alpha;
  52 }
  53
  54 void itkFM3D::SetBeta(double beta){
  55         this->beta = beta;
  56 }
  57
  58
  59 void itkFM3D::SetStopTime(double stopTime){
  60         this->stopTime = stopTime;
  61 }
  62
  63 void itkFM3D::CurvatureAnisotropicFiltertOn(){
  64         usefilter = 1;
  65 }
  66
  67 void itkFM3D::CurvatureAnisotropicFilterOff(){
  68         usefilter = 0;
  69 }
  70
  71
  72 vtkImageData* itkFM3D::segment(vtkImageData *volume){
  73
  74         volume->SetSpacing(1,1,1);
  75
  76
  77         double spc[3];
  78         volume->GetSpacing(spc);
  79
  80
  81         double xx = x; // * spc[0];
  82         double yy = y; // * spc[1];
  83         double zz = z; // * spc[2];
  84
  85
  86         /*FILE *ff;
  87         ff=fopen("programa/Ups.bat","w");
  88         fprintf(ff,"vtkViewer imagen.vtk %d %d %d %f %f %f \n",(int)xx,(int)yy,(int)zz,(float)alpha,(float)beta,(float)stopTime);
  89         fclose(ff);
  90     int code = wxExecute("programa/Ups.bat", wxEXEC_SYNC);*/
  91
  92
  93
  94
  95
  96
  97   //DEFINICION IMAGEN ENTRADA
  98   //--------------------------------------------------------------------------------------------
  99   typedef   double           InternalPixelType;
 100   const     unsigned int    DIMENSION_IMAGEN = 3;
 101   typedef itk::Image< InternalPixelType, DIMENSION_IMAGEN >  InternalImageType;
 102
 103
 104
 105
 106   // CONVERTIR LA IMAGEN DE ENTRADA A ITK PARA SER PROCESADA
 107   //--------------------------------------------------------------------------------------------
 108
 109   vtkImageCast *ic = vtkImageCast::New();
 110   ic->SetInput(volume);
 111   ic->SetOutputScalarTypeToDouble();
 112   ic->Update();
 113   typedef itk::VTKImageToImageFilter<InternalImageType> VtkToItkFilterType;
 114   VtkToItkFilterType::Pointer toItk = VtkToItkFilterType::New();
 115   toItk->SetInput( ic->GetOutput() );
 116   toItk->Update();
 117    const InternalImageType* itkImageData = toItk->GetOutput();
 118   itkImageData->Register();
 119
 120   //typedef itk::ImageFileReader<InternalImageType> ReaderType;
 121   //ReaderType::Pointer reader = ReaderType::New();
 122   //reader->SetFileName("imagen.vtk");
 123
 124   //reader->Update();
 125
 126   //InternalImageType* itkImageData = null; //reader->GetOutput();
 127
 128
 129
 130
 131   //DEFINICION WRITER
 132   //--------------------------------------------------------------------------------------------
 133  // typedef itk::ImageFileWriter<InternalImageType> WriterType;
 134
 135
 136
 137   //DEFINICION DEL FILTRADO ANISOTROPICO
 138   //--------------------------------------------------------------------------------------------
 139   typedef itk::CurvatureAnisotropicDiffusionImageFilter<
 140                                 InternalImageType,
 141                                 InternalImageType > SmoothingFilterType;
 142
 143   SmoothingFilterType::Pointer filtroAnisotropico = SmoothingFilterType::New();
 144
 145
 146
 147   //DEFINICION DEL FILTRO SIGMOIDE
 148   //--------------------------------------------------------------------------------------------
 149   typedef   itk::SigmoidImageFilter<
 150                                InternalImageType,
 151                                InternalImageType >  SigmoidFilterType;
 152
 153   SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
 154
 155   sigmoid->SetOutputMinimum(  0.0  );
 156   sigmoid->SetOutputMaximum(  255  );
 157   sigmoid->SetAlpha(this->alpha);
 158   sigmoid->SetBeta(this->beta);
 159
 160
 161
 162   //DECLARACION DEL FILTRO DE FAST MARCHING
 163   //--------------------------------------------------------------------------------------------
 164   typedef  itk::FastMarchingImageFilter< InternalImageType,
 165                               InternalImageType >    FastMarchingFilterType;
 166   FastMarchingFilterType::Pointer  fastMarching = FastMarchingFilterType::New();
 167
 168   //DECLARACION DEL FILTRO THRESHOLDING BINARIO
 169   //--------------------------------------------------------------------------------------------
 170   typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType > ThresholdingFilterType;
 171   ThresholdingFilterType::Pointer thresholdFilter = ThresholdingFilterType::New();
 172
 173
 174   //--------------------------------------------------------------------------------------------
 175   // ESTABLECIMIENTO DE LA SEMILLA DEL FAST MARCHING
 176   //--------------------------------------------------------------------------------------------
 177   typedef FastMarchingFilterType::NodeContainer           NodeContainer;
 178   typedef FastMarchingFilterType::NodeType                NodeType;
 179   NodeContainer::Pointer seeds = NodeContainer::New();
 180
 181   InternalImageType::IndexType  seedPosition;
 182   seedPosition[0] = this->x; // POSICION X DE LA SEMILLA
 183   seedPosition[1] = this->y; // POSICION Y DE LA SEMILLA
 184   seedPosition[2] = this->z; // POSICION Z DE LA SEMILLA
 185
 186   NodeType node;
 187   const double seedValue = 0.0;
 188
 189   node.SetValue( seedValue );
 190   node.SetIndex( seedPosition );
 191
 192   seeds->Initialize();
 193   seeds->InsertElement( 0, node );
 194
 195
 196   //COLOCAR VARIAS SEMILLAS
 197   fastMarching->SetTrialPoints(  seeds  );
 198   fastMarching->SetOutputSize(itkImageData->GetBufferedRegion().GetSize() );
 199   fastMarching->SetStoppingValue(this->stopTime);
 200
 201
 202   //DEFINICION WRITER FILTRADO ANISOTROPICO
 203   //WriterType::Pointer writerFiltroRuido = WriterType::New();
 204   //writerFiltroRuido->SetFileName("filtrado.vtk");
 205
 206
 207   //DEFINICION WRITER SIGMOIDE
 208   //--------------------------------------------------------------------------------------------
 209   //WriterType::Pointer writerSigmoide = WriterType::New();
 210   //writerSigmoide->SetFileName("sigmoid.vtk");
 211
 212
 213
 214   //DEFINICION WRITER FAST MARCHING MAPA DE TIEMPOS
 215   //--------------------------------------------------------------------------------------------
 216   /*typedef itk::VTKImageIO ImageIOType;
 217   ImageIOType::Pointer vtkIO = ImageIOType::New();
 218   vtkIO->SetFileTypeToASCII();
 219
 220   WriterType::Pointer writerMapaTiempos = WriterType::New();
 221   writerMapaTiempos->SetFileName("fastmarch.vtk");
 222   writerMapaTiempos->SetImageIO(vtkIO);*/
 223
 224
 225   //DEFINICION WRITER THRESHOLD BINARIO
 226   //--------------------------------------------------------------------------------------------
 227 //  WriterType::Pointer writerThreshold = WriterType::New();
 228  // writerThreshold->SetFileName("binario.vtk");
 229   //writerThreshold->SetImageIO(vtkIO);
 230
 231
 232   //--------------------------------------------------------------------------------------------
 233   // CONSTRUCCION DEL PIPELINE DE PROCESAMIENTO
 234   //--------------------------------------------------------------------------------------------
 235
 236   // READER -> FILTRO ANISOTROPICO
 237   //FILTRO ANISOTROPICO -> SIGMOIDE
 238   if (usefilter == 1){
 239         //filtroAnisotropico->SetInput(itkImageData);
 240         //sigmoid->SetInput(filtroAnisotropico->GetOutput());
 241           sigmoid->SetInput(itkImageData);
 242   }
 243   else {
 244           sigmoid->SetInput(itkImageData);
 245   }
 246
 247   // SIGMOIDE -> FAST MARCHING
 248   fastMarching->SetInput(sigmoid->GetOutput());
 249
 250
 251   // FAST MARCHING -> WRITER
 252   //writerFiltroRuido->SetInput(filtroAnisotropico->GetOutput());
 253   //writerSigmoide->SetInput(sigmoid->GetOutput());
 254   //writerMapaTiempos->SetInput(fastMarching->GetOutput());
 255   //writerFiltroRuido->Update();
 256
 257
 258   //----------------------------------------------------------------------------------------
 259   //Ejecucion de OTSU
 260   //----------------------------------------------------------------------------------------
 261
 262   //DEFINICION DE IMAGEN VTK FILTRADA CON ANISOTROPICO
 263   //--------------------------------------------------------------------------------------------
 264   typedef itk::ImageToVTKImageFilter<InternalImageType> ItkToVtkFilterType;
 265
 266
 267
 268   /*vtkOtsuImageData* otsu = new vtkOtsuImageData();
 269
 270   if(this->usefilter){
 271         printf("Calculando Th Otsu sobre la imagen filtrada\n");
 272          ItkToVtkFilterType::Pointer toVtk = ItkToVtkFilterType::New();
 273         toVtk->SetInput(filtroAnisotropico->GetOutput());
 274         toVtk->Update();
 275         vtkImageData* imagenFiltrada = toVtk->GetOutput();
 276         estimatedOtsuThreshold = otsu->calculateOptimalThreshold(imagenFiltrada);
 277   }
 278   else{
 279     printf("Calculando Th Otsu sobre la imagen original\n");
 280         estimatedOtsuThreshold = otsu->calculateOptimalThreshold(volume);
 281   }
 282
 283   printf("BETA = %.0f, OTSU = %d\n",beta,estimatedOtsuThreshold);*/
 284
 285   //sigmoid->SetBeta((double)ThresholdOtsu);
 286
 287   //----------------------------------------------------------------------------------------
 288
 289
 290   //--------------------------------------------------------------------------------------------
 291   // EJECUCION DEL PIPELINE
 292   //--------------------------------------------------------------------------------------------
 293   try
 294   {
 295           fastMarching->Update();
 296           //writerFiltroRuido->Update();
 297           //writerSigmoide->Update();
 298           //writerMapaTiempos->Update();
 299       printf("Writers (ruido,sigmoide, mapatiempos) procesados \n");
 300   }
 301   catch( itk::ExceptionObject & excep )
 302   {
 303     std::cerr << "Ocurrio un problema" << std::endl;
 304     std::cerr << excep << std::endl;
 305   }
 306
 307
 308   //---------------------------------------------------------
 309   // LEYENDO MAPA DE TIEMPOS
 310   //---------------------------------------------------------
 311   // En esta parte escribimos el mapa de tiempos a un archivo plano para
 312   // poder estudiarlo en Excel.
 313
 314     printf("Procesando el mapa de tiempos....");
 315
 316         //ACTIVAR PARA DEBUG
 317         //-----------------------------
 318         //FILE *logger = NULL;
 319         //      if ((logger = freopen("mapatiempos.txt","w", stdout)) == NULL)
 320         //              exit(-1);
 321
 322
 323
 324         InternalImageType::Pointer mapa = fastMarching->GetOutput();
 325
 326         typedef itk::ImageRegionConstIterator< InternalImageType > ConstIteratorType;
 327
 328         ConstIteratorType constIterator( mapa, mapa->GetLargestPossibleRegion() );
 329
 330
 331         // inicializar histogramas
 332         double valor = 0.0;
 333         int indice = 0;
 334
 335         int histoTiempos[1000];
 336         int histoAcumulado[1000];
 337
 338         for (indice = 0; indice <= 1000; indice++){
 339                 histoTiempos[indice] = 0;
 340                 histoAcumulado[indice] = 0;
 341         }
 342
 343
 344         // calcular histograma
 345         for (constIterator.GoToBegin(); !constIterator.IsAtEnd(); ++ constIterator){
 346                 valor = constIterator.Get();
 347                 if (valor <= 10){
 348                         indice = (int)(valor * 100);
 349                         histoTiempos[indice]++;
 350
 351                 }
 352         }
 353
 354         //calcular histograma acumulado
 355         for (indice = 1 ; indice <= 1000; indice++){
 356                 histoAcumulado[indice]= histoAcumulado[indice-1]+ histoTiempos[indice];
 357         }
 358
 359         //buscar la frecuencia minima = derivada de la funcion de frecuencia acumulada
 360         int frecuenciaMinima = 1000000000;
 361         int indiceBuscado = -1;
 362         for (indice = 0; indice <= 1000; indice++){
 363                 if (histoTiempos[indice] < frecuenciaMinima){
 364                         frecuenciaMinima = histoTiempos[indice];
 365                         indiceBuscado = indice;
 366                 }
 367
 368         }
 369
 370
 371         double valorThreshold = (double)indiceBuscado / 100.0;
 372
 373         printf("VALOR DE SEGMENTACION: %.2f\n\n",valorThreshold);
 374
 375
 376         //for (indice = 0; indice <= 1000; indice++){
 377         //      printf("%.2f:%d:%d\n",(double)indice/100.0, histoTiempos[indice], histoAcumulado[indice]);
 378         //}
 379
 380
 381         //fclose(logger);
 382         //logger = NULL;
 383
 384         printf("Mapa de tiempos OK");
 385
 386
 387         thresholdFilter->SetOutsideValue( 0 );
 388         thresholdFilter->SetInsideValue( 255 );
 389         thresholdFilter->SetLowerThreshold(0);
 390         thresholdFilter->SetUpperThreshold(valorThreshold);
 391
 392         thresholdFilter->SetInput(fastMarching->GetOutput());
 393 //      writerThreshold->SetInput(thresholdFilter->GetOutput());
 394
 395         try{
 396                 thresholdFilter->Update();
 397                 //writerThreshold->Update();
 398         }
 399         catch( itk::ExceptionObject & excep )
 400         {
 401                 std::cerr << "Ocurrio un problema" << std::endl;
 402                 std::cerr << excep << std::endl;
 403         }
 404
 405         ItkToVtkFilterType::Pointer resultado = ItkToVtkFilterType::New();
 406         resultado->SetInput(thresholdFilter->GetOutput());
 407         vtkImageData* salida = resultado->GetOutput();
 408         resultado->Register();
 409         salida->Update();
 410         salida->ReleaseDataFlagOff();
 411
 412
 413
 414
 415         return salida;
 416 }