1 //--------------------------------------------------------------------------------------------
3 //--------------------------------------------------------------------------------------------
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"
18 #include "vtkOtsuImageData.h"
19 #include "vtkImageCast.h"
20 #include "vtkMetaImageWriter.h"
22 #include "vtkDataSetWriter.h"
25 #include <wx/arrstr.h>
26 #include <wx/process.h>
30 usefilter = 1; // Por defecto hacer el filtrado de ruido con la curvatura anisotropica.
33 stopTime = 10; // Valor por defecto;
43 void itkFM3D::AddSeed(int x, int y, int z){
49 void itkFM3D::SetAlpha(double alpha){
54 void itkFM3D::SetBeta(double beta){
59 void itkFM3D::SetStopTime(double stopTime){
60 this->stopTime = stopTime;
63 void itkFM3D::CurvatureAnisotropicFiltertOn(){
67 void itkFM3D::CurvatureAnisotropicFilterOff(){
72 vtkImageData* itkFM3D::segment(vtkImageData *volume){
74 volume->SetSpacing(1,1,1);
78 volume->GetSpacing(spc);
81 double xx = x; // * spc[0];
82 double yy = y; // * spc[1];
83 double zz = z; // * spc[2];
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);
90 int code = wxExecute("programa/Ups.bat", wxEXEC_SYNC);*/
97 //DEFINICION IMAGEN ENTRADA
98 //--------------------------------------------------------------------------------------------
99 typedef double InternalPixelType;
100 const unsigned int DIMENSION_IMAGEN = 3;
101 typedef itk::Image< InternalPixelType, DIMENSION_IMAGEN > InternalImageType;
106 // CONVERTIR LA IMAGEN DE ENTRADA A ITK PARA SER PROCESADA
107 //--------------------------------------------------------------------------------------------
109 vtkImageCast *ic = vtkImageCast::New();
110 ic->SetInput(volume);
111 ic->SetOutputScalarTypeToDouble();
113 typedef itk::VTKImageToImageFilter<InternalImageType> VtkToItkFilterType;
114 VtkToItkFilterType::Pointer toItk = VtkToItkFilterType::New();
115 toItk->SetInput( ic->GetOutput() );
117 const InternalImageType* itkImageData = toItk->GetOutput();
118 itkImageData->Register();
120 //typedef itk::ImageFileReader<InternalImageType> ReaderType;
121 //ReaderType::Pointer reader = ReaderType::New();
122 //reader->SetFileName("imagen.vtk");
126 //InternalImageType* itkImageData = null; //reader->GetOutput();
132 //--------------------------------------------------------------------------------------------
133 // typedef itk::ImageFileWriter<InternalImageType> WriterType;
137 //DEFINICION DEL FILTRADO ANISOTROPICO
138 //--------------------------------------------------------------------------------------------
139 typedef itk::CurvatureAnisotropicDiffusionImageFilter<
141 InternalImageType > SmoothingFilterType;
143 SmoothingFilterType::Pointer filtroAnisotropico = SmoothingFilterType::New();
147 //DEFINICION DEL FILTRO SIGMOIDE
148 //--------------------------------------------------------------------------------------------
149 typedef itk::SigmoidImageFilter<
151 InternalImageType > SigmoidFilterType;
153 SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
155 sigmoid->SetOutputMinimum( 0.0 );
156 sigmoid->SetOutputMaximum( 255 );
157 sigmoid->SetAlpha(this->alpha);
158 sigmoid->SetBeta(this->beta);
162 //DECLARACION DEL FILTRO DE FAST MARCHING
163 //--------------------------------------------------------------------------------------------
164 typedef itk::FastMarchingImageFilter< InternalImageType,
165 InternalImageType > FastMarchingFilterType;
166 FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();
168 //DECLARACION DEL FILTRO THRESHOLDING BINARIO
169 //--------------------------------------------------------------------------------------------
170 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType > ThresholdingFilterType;
171 ThresholdingFilterType::Pointer thresholdFilter = ThresholdingFilterType::New();
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();
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
187 const double seedValue = 0.0;
189 node.SetValue( seedValue );
190 node.SetIndex( seedPosition );
193 seeds->InsertElement( 0, node );
196 //COLOCAR VARIAS SEMILLAS
197 fastMarching->SetTrialPoints( seeds );
198 fastMarching->SetOutputSize(itkImageData->GetBufferedRegion().GetSize() );
199 fastMarching->SetStoppingValue(this->stopTime);
202 //DEFINICION WRITER FILTRADO ANISOTROPICO
203 //WriterType::Pointer writerFiltroRuido = WriterType::New();
204 //writerFiltroRuido->SetFileName("filtrado.vtk");
207 //DEFINICION WRITER SIGMOIDE
208 //--------------------------------------------------------------------------------------------
209 //WriterType::Pointer writerSigmoide = WriterType::New();
210 //writerSigmoide->SetFileName("sigmoid.vtk");
214 //DEFINICION WRITER FAST MARCHING MAPA DE TIEMPOS
215 //--------------------------------------------------------------------------------------------
216 /*typedef itk::VTKImageIO ImageIOType;
217 ImageIOType::Pointer vtkIO = ImageIOType::New();
218 vtkIO->SetFileTypeToASCII();
220 WriterType::Pointer writerMapaTiempos = WriterType::New();
221 writerMapaTiempos->SetFileName("fastmarch.vtk");
222 writerMapaTiempos->SetImageIO(vtkIO);*/
225 //DEFINICION WRITER THRESHOLD BINARIO
226 //--------------------------------------------------------------------------------------------
227 // WriterType::Pointer writerThreshold = WriterType::New();
228 // writerThreshold->SetFileName("binario.vtk");
229 //writerThreshold->SetImageIO(vtkIO);
232 //--------------------------------------------------------------------------------------------
233 // CONSTRUCCION DEL PIPELINE DE PROCESAMIENTO
234 //--------------------------------------------------------------------------------------------
236 // READER -> FILTRO ANISOTROPICO
237 //FILTRO ANISOTROPICO -> SIGMOIDE
239 //filtroAnisotropico->SetInput(itkImageData);
240 //sigmoid->SetInput(filtroAnisotropico->GetOutput());
241 sigmoid->SetInput(itkImageData);
244 sigmoid->SetInput(itkImageData);
247 // SIGMOIDE -> FAST MARCHING
248 fastMarching->SetInput(sigmoid->GetOutput());
251 // FAST MARCHING -> WRITER
252 //writerFiltroRuido->SetInput(filtroAnisotropico->GetOutput());
253 //writerSigmoide->SetInput(sigmoid->GetOutput());
254 //writerMapaTiempos->SetInput(fastMarching->GetOutput());
255 //writerFiltroRuido->Update();
258 //----------------------------------------------------------------------------------------
260 //----------------------------------------------------------------------------------------
262 //DEFINICION DE IMAGEN VTK FILTRADA CON ANISOTROPICO
263 //--------------------------------------------------------------------------------------------
264 typedef itk::ImageToVTKImageFilter<InternalImageType> ItkToVtkFilterType;
268 /*vtkOtsuImageData* otsu = new vtkOtsuImageData();
271 printf("Calculando Th Otsu sobre la imagen filtrada\n");
272 ItkToVtkFilterType::Pointer toVtk = ItkToVtkFilterType::New();
273 toVtk->SetInput(filtroAnisotropico->GetOutput());
275 vtkImageData* imagenFiltrada = toVtk->GetOutput();
276 estimatedOtsuThreshold = otsu->calculateOptimalThreshold(imagenFiltrada);
279 printf("Calculando Th Otsu sobre la imagen original\n");
280 estimatedOtsuThreshold = otsu->calculateOptimalThreshold(volume);
283 printf("BETA = %.0f, OTSU = %d\n",beta,estimatedOtsuThreshold);*/
285 //sigmoid->SetBeta((double)ThresholdOtsu);
287 //----------------------------------------------------------------------------------------
290 //--------------------------------------------------------------------------------------------
291 // EJECUCION DEL PIPELINE
292 //--------------------------------------------------------------------------------------------
295 fastMarching->Update();
296 //writerFiltroRuido->Update();
297 //writerSigmoide->Update();
298 //writerMapaTiempos->Update();
299 printf("Writers (ruido,sigmoide, mapatiempos) procesados \n");
301 catch( itk::ExceptionObject & excep )
303 std::cerr << "Ocurrio un problema" << std::endl;
304 std::cerr << excep << std::endl;
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.
314 printf("Procesando el mapa de tiempos....");
317 //-----------------------------
318 //FILE *logger = NULL;
319 // if ((logger = freopen("mapatiempos.txt","w", stdout)) == NULL)
324 InternalImageType::Pointer mapa = fastMarching->GetOutput();
326 typedef itk::ImageRegionConstIterator< InternalImageType > ConstIteratorType;
328 ConstIteratorType constIterator( mapa, mapa->GetLargestPossibleRegion() );
331 // inicializar histogramas
335 int histoTiempos[1000];
336 int histoAcumulado[1000];
338 for (indice = 0; indice <= 1000; indice++){
339 histoTiempos[indice] = 0;
340 histoAcumulado[indice] = 0;
344 // calcular histograma
345 for (constIterator.GoToBegin(); !constIterator.IsAtEnd(); ++ constIterator){
346 valor = constIterator.Get();
348 indice = (int)(valor * 100);
349 histoTiempos[indice]++;
354 //calcular histograma acumulado
355 for (indice = 1 ; indice <= 1000; indice++){
356 histoAcumulado[indice]= histoAcumulado[indice-1]+ histoTiempos[indice];
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;
371 double valorThreshold = (double)indiceBuscado / 100.0;
373 printf("VALOR DE SEGMENTACION: %.2f\n\n",valorThreshold);
376 //for (indice = 0; indice <= 1000; indice++){
377 // printf("%.2f:%d:%d\n",(double)indice/100.0, histoTiempos[indice], histoAcumulado[indice]);
384 printf("Mapa de tiempos OK");
387 thresholdFilter->SetOutsideValue( 0 );
388 thresholdFilter->SetInsideValue( 255 );
389 thresholdFilter->SetLowerThreshold(0);
390 thresholdFilter->SetUpperThreshold(valorThreshold);
392 thresholdFilter->SetInput(fastMarching->GetOutput());
393 // writerThreshold->SetInput(thresholdFilter->GetOutput());
396 thresholdFilter->Update();
397 //writerThreshold->Update();
399 catch( itk::ExceptionObject & excep )
401 std::cerr << "Ocurrio un problema" << std::endl;
402 std::cerr << excep << std::endl;
405 ItkToVtkFilterType::Pointer resultado = ItkToVtkFilterType::New();
406 resultado->SetInput(thresholdFilter->GetOutput());
407 vtkImageData* salida = resultado->GetOutput();
408 resultado->Register();
410 salida->ReleaseDataFlagOff();