1 /*# ---------------------------------------------------------------------
3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
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
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.
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
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 # ------------------------------------------------------------------------ */
26 //--------------------------------------------------------------------------------------------
28 //--------------------------------------------------------------------------------------------
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"
43 #include "vtkOtsuImageData.h"
44 #include "vtkImageCast.h"
45 #include "vtkMetaImageWriter.h"
47 #include "vtkDataSetWriter.h"
50 #include <wx/arrstr.h>
51 #include <wx/process.h>
55 usefilter = 1; // Por defecto hacer el filtrado de ruido con la curvatura anisotropica.
58 stopTime = 10; // Valor por defecto;
68 void itkFM3D::AddSeed(int x, int y, int z){
74 void itkFM3D::SetAlpha(double alpha){
79 void itkFM3D::SetBeta(double beta){
84 void itkFM3D::SetStopTime(double stopTime){
85 this->stopTime = stopTime;
88 void itkFM3D::CurvatureAnisotropicFiltertOn(){
92 void itkFM3D::CurvatureAnisotropicFilterOff(){
97 vtkImageData* itkFM3D::segment(vtkImageData *volume){
99 volume->SetSpacing(1,1,1);
103 volume->GetSpacing(spc);
106 double xx = x; // * spc[0];
107 double yy = y; // * spc[1];
108 double zz = z; // * spc[2];
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);
115 int code = wxExecute("programa/Ups.bat", wxEXEC_SYNC);*/
122 //DEFINICION IMAGEN ENTRADA
123 //--------------------------------------------------------------------------------------------
124 typedef double InternalPixelType;
125 const unsigned int DIMENSION_IMAGEN = 3;
126 typedef itk::Image< InternalPixelType, DIMENSION_IMAGEN > InternalImageType;
131 // CONVERTIR LA IMAGEN DE ENTRADA A ITK PARA SER PROCESADA
132 //--------------------------------------------------------------------------------------------
134 vtkImageCast *ic = vtkImageCast::New();
135 ic->SetInput(volume);
136 ic->SetOutputScalarTypeToDouble();
138 typedef itk::VTKImageToImageFilter<InternalImageType> VtkToItkFilterType;
139 VtkToItkFilterType::Pointer toItk = VtkToItkFilterType::New();
140 toItk->SetInput( ic->GetOutput() );
142 const InternalImageType* itkImageData = toItk->GetOutput();
143 itkImageData->Register();
145 //typedef itk::ImageFileReader<InternalImageType> ReaderType;
146 //ReaderType::Pointer reader = ReaderType::New();
147 //reader->SetFileName("imagen.vtk");
151 //InternalImageType* itkImageData = null; //reader->GetOutput();
157 //--------------------------------------------------------------------------------------------
158 // typedef itk::ImageFileWriter<InternalImageType> WriterType;
162 //DEFINICION DEL FILTRADO ANISOTROPICO
163 //--------------------------------------------------------------------------------------------
164 typedef itk::CurvatureAnisotropicDiffusionImageFilter<
166 InternalImageType > SmoothingFilterType;
168 SmoothingFilterType::Pointer filtroAnisotropico = SmoothingFilterType::New();
172 //DEFINICION DEL FILTRO SIGMOIDE
173 //--------------------------------------------------------------------------------------------
174 typedef itk::SigmoidImageFilter<
176 InternalImageType > SigmoidFilterType;
178 SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
180 sigmoid->SetOutputMinimum( 0.0 );
181 sigmoid->SetOutputMaximum( 255 );
182 sigmoid->SetAlpha(this->alpha);
183 sigmoid->SetBeta(this->beta);
187 //DECLARACION DEL FILTRO DE FAST MARCHING
188 //--------------------------------------------------------------------------------------------
189 typedef itk::FastMarchingImageFilter< InternalImageType,
190 InternalImageType > FastMarchingFilterType;
191 FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();
193 //DECLARACION DEL FILTRO THRESHOLDING BINARIO
194 //--------------------------------------------------------------------------------------------
195 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType > ThresholdingFilterType;
196 ThresholdingFilterType::Pointer thresholdFilter = ThresholdingFilterType::New();
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();
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
212 const double seedValue = 0.0;
214 node.SetValue( seedValue );
215 node.SetIndex( seedPosition );
218 seeds->InsertElement( 0, node );
221 //COLOCAR VARIAS SEMILLAS
222 fastMarching->SetTrialPoints( seeds );
223 fastMarching->SetOutputSize(itkImageData->GetBufferedRegion().GetSize() );
224 fastMarching->SetStoppingValue(this->stopTime);
227 //DEFINICION WRITER FILTRADO ANISOTROPICO
228 //WriterType::Pointer writerFiltroRuido = WriterType::New();
229 //writerFiltroRuido->SetFileName("filtrado.vtk");
232 //DEFINICION WRITER SIGMOIDE
233 //--------------------------------------------------------------------------------------------
234 //WriterType::Pointer writerSigmoide = WriterType::New();
235 //writerSigmoide->SetFileName("sigmoid.vtk");
239 //DEFINICION WRITER FAST MARCHING MAPA DE TIEMPOS
240 //--------------------------------------------------------------------------------------------
241 /*typedef itk::VTKImageIO ImageIOType;
242 ImageIOType::Pointer vtkIO = ImageIOType::New();
243 vtkIO->SetFileTypeToASCII();
245 WriterType::Pointer writerMapaTiempos = WriterType::New();
246 writerMapaTiempos->SetFileName("fastmarch.vtk");
247 writerMapaTiempos->SetImageIO(vtkIO);*/
250 //DEFINICION WRITER THRESHOLD BINARIO
251 //--------------------------------------------------------------------------------------------
252 // WriterType::Pointer writerThreshold = WriterType::New();
253 // writerThreshold->SetFileName("binario.vtk");
254 //writerThreshold->SetImageIO(vtkIO);
257 //--------------------------------------------------------------------------------------------
258 // CONSTRUCCION DEL PIPELINE DE PROCESAMIENTO
259 //--------------------------------------------------------------------------------------------
261 // READER -> FILTRO ANISOTROPICO
262 //FILTRO ANISOTROPICO -> SIGMOIDE
264 //filtroAnisotropico->SetInput(itkImageData);
265 //sigmoid->SetInput(filtroAnisotropico->GetOutput());
266 sigmoid->SetInput(itkImageData);
269 sigmoid->SetInput(itkImageData);
272 // SIGMOIDE -> FAST MARCHING
273 fastMarching->SetInput(sigmoid->GetOutput());
276 // FAST MARCHING -> WRITER
277 //writerFiltroRuido->SetInput(filtroAnisotropico->GetOutput());
278 //writerSigmoide->SetInput(sigmoid->GetOutput());
279 //writerMapaTiempos->SetInput(fastMarching->GetOutput());
280 //writerFiltroRuido->Update();
283 //----------------------------------------------------------------------------------------
285 //----------------------------------------------------------------------------------------
287 //DEFINICION DE IMAGEN VTK FILTRADA CON ANISOTROPICO
288 //--------------------------------------------------------------------------------------------
289 typedef itk::ImageToVTKImageFilter<InternalImageType> ItkToVtkFilterType;
293 /*vtkOtsuImageData* otsu = new vtkOtsuImageData();
296 printf("Calculando Th Otsu sobre la imagen filtrada\n");
297 ItkToVtkFilterType::Pointer toVtk = ItkToVtkFilterType::New();
298 toVtk->SetInput(filtroAnisotropico->GetOutput());
300 vtkImageData* imagenFiltrada = toVtk->GetOutput();
301 estimatedOtsuThreshold = otsu->calculateOptimalThreshold(imagenFiltrada);
304 printf("Calculando Th Otsu sobre la imagen original\n");
305 estimatedOtsuThreshold = otsu->calculateOptimalThreshold(volume);
308 printf("BETA = %.0f, OTSU = %d\n",beta,estimatedOtsuThreshold);*/
310 //sigmoid->SetBeta((double)ThresholdOtsu);
312 //----------------------------------------------------------------------------------------
315 //--------------------------------------------------------------------------------------------
316 // EJECUCION DEL PIPELINE
317 //--------------------------------------------------------------------------------------------
320 fastMarching->Update();
321 //writerFiltroRuido->Update();
322 //writerSigmoide->Update();
323 //writerMapaTiempos->Update();
324 printf("Writers (ruido,sigmoide, mapatiempos) procesados \n");
326 catch( itk::ExceptionObject & excep )
328 std::cerr << "Ocurrio un problema" << std::endl;
329 std::cerr << excep << std::endl;
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.
339 printf("Procesando el mapa de tiempos....");
342 //-----------------------------
343 //FILE *logger = NULL;
344 // if ((logger = freopen("mapatiempos.txt","w", stdout)) == NULL)
349 InternalImageType::Pointer mapa = fastMarching->GetOutput();
351 typedef itk::ImageRegionConstIterator< InternalImageType > ConstIteratorType;
353 ConstIteratorType constIterator( mapa, mapa->GetLargestPossibleRegion() );
356 // inicializar histogramas
360 int histoTiempos[1000];
361 int histoAcumulado[1000];
363 for (indice = 0; indice <= 1000; indice++){
364 histoTiempos[indice] = 0;
365 histoAcumulado[indice] = 0;
369 // calcular histograma
370 for (constIterator.GoToBegin(); !constIterator.IsAtEnd(); ++ constIterator){
371 valor = constIterator.Get();
373 indice = (int)(valor * 100);
374 histoTiempos[indice]++;
379 //calcular histograma acumulado
380 for (indice = 1 ; indice <= 1000; indice++){
381 histoAcumulado[indice]= histoAcumulado[indice-1]+ histoTiempos[indice];
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;
396 double valorThreshold = (double)indiceBuscado / 100.0;
398 printf("VALOR DE SEGMENTACION: %.2f\n\n",valorThreshold);
401 //for (indice = 0; indice <= 1000; indice++){
402 // printf("%.2f:%d:%d\n",(double)indice/100.0, histoTiempos[indice], histoAcumulado[indice]);
409 printf("Mapa de tiempos OK");
412 thresholdFilter->SetOutsideValue( 0 );
413 thresholdFilter->SetInsideValue( 255 );
414 thresholdFilter->SetLowerThreshold(0);
415 thresholdFilter->SetUpperThreshold(valorThreshold);
417 thresholdFilter->SetInput(fastMarching->GetOutput());
418 // writerThreshold->SetInput(thresholdFilter->GetOutput());
421 thresholdFilter->Update();
422 //writerThreshold->Update();
424 catch( itk::ExceptionObject & excep )
426 std::cerr << "Ocurrio un problema" << std::endl;
427 std::cerr << excep << std::endl;
430 ItkToVtkFilterType::Pointer resultado = ItkToVtkFilterType::New();
431 resultado->SetInput(thresholdFilter->GetOutput());
432 vtkImageData* salida = resultado->GetOutput();
433 resultado->Register();
435 salida->ReleaseDataFlagOff();