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

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