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 # ------------------------------------------------------------------------ */
27 #include "vtkAppendFilter.h"
28 #include "vtkCamera.h"
29 #include "vtkCallbackCommand.h"
30 #include "vtkCellArray.h"
31 #include "vtkCylinderSource.h"
32 #include "vtkClosePolyData.h"
33 #include "vtkCommand.h"
34 #include "vtkDoubleArray.h"
35 #include "vtkExtractVOI.h"
36 #include "vtkImageReader.h"
37 #include "vtkImageViewer.h"
38 #include "vtkImageViewer2.h"
39 #include "vtkImageToStructuredPoints.h"
40 #include "vtkImageThreshold.h"
41 #include "vtkImageWriter.h"
42 #include "vtkImageClip.h"
43 #include "vtkImageResample.h"
44 #include "vtkImageReslice.h"
45 #include "vtkImageResample.h"
46 #include "vtkImageThreshold.h"
47 #include "vtkImageCast.h"
48 #include "vtkImageData.h"
50 #include "vtkMarchingCubes.h"
51 #include "vtkObjectFactory.h"
52 #include "vtkPolyDataMapper.h"
53 #include "vtkRenderer.h"
54 #include "vtkRenderWindow.h"
55 #include "vtkRenderWindowInteractor.h"
56 #include "vtkPolyDataConnectivityFilter.h"
57 #include "vtkProperty.h"
58 #include "vtkPriorityQueue.h"
59 #include "vtkPoints.h"
60 #include "vtkPolyData.h"
61 #include "vtkPolyDataMapper.h"
62 #include "vtkPolyDataWriter.h"
63 #include "vtkPolyDataReader.h"
64 #include "vtkPointData.h"
65 #include "vtkSphereSource.h"
66 #include "vtkSTLWriter.h"
67 #include "vtkStripper.h"
68 #include "vtkTriangleFilter.h"
69 #include "vtkTransform.h"
71 #include "wxPathologyWidget_01.h"
72 #include "kernel/vtkOtsuSphereSource.h"
73 #include <wx/splitter.h>
78 //-------------------------------------------------------------------
79 //-------------------------------------------------------------------
80 //-------------------------------------------------------------------
81 wxPathologyWidget_01::wxPathologyWidget_01(wxWindow *parent, marInterface* mar)
82 : wxPanel( parent, -1)
87 wxBoxSizer *sizer = new wxBoxSizer(wxVERTICAL );
88 wxSplitterWindow *pnlSplitter = new wxSplitterWindow( this , -1);
89 wxPanel *viewPanel = CreateViewPanel(pnlSplitter);
90 wxPanel *controlPanel = CreateControlPanel(pnlSplitter);
92 sizer -> Add( pnlSplitter ,1,wxGROW ,0);
93 pnlSplitter -> SetMinimumPaneSize( 150 );
94 pnlSplitter -> SplitVertically( viewPanel, controlPanel, 600 );
95 this -> SetSizer(sizer);
106 p1SphereSource = NULL;
110 p2SphereSource = NULL;
114 patSphereSource = NULL;
119 pathologyFrame = NULL;
127 outlineFilter = NULL;
128 outlineMapper = NULL;
132 p3SphereSource = NULL;
138 //-------------------------------------------------------------------
139 wxPathologyWidget_01::~wxPathologyWidget_01(){
142 //-------------------------------------------------------------------
143 wxPanel* wxPathologyWidget_01::CreateViewPanel(wxWindow *parent)
145 wxPanel *panel = new wxPanel(parent,-1);
146 wxBoxSizer *sizer = new wxBoxSizer(wxVERTICAL);
147 _maracasSurfaceWidget = new wxSurfaceWidget(panel);
148 sizer->Add(_maracasSurfaceWidget , 1, wxEXPAND, 0);
149 panel->SetSizer(sizer);
150 panel->SetAutoLayout(true);
151 panel->SetSize(400,400);
155 //-------------------------------------------------------------------
156 wxPanel* wxPathologyWidget_01::CreateControlPanel(wxWindow *parent)
158 wxPanel *panel = new wxPanel(parent,-1);
159 wxFlexGridSizer *sizer = new wxFlexGridSizer(2);
162 //-- PATHOLOGY WIDGETS
164 sizer->Add(new wxStaticText(panel,-1,_T(" ")));
165 sizer->Add(new wxStaticText(panel,-1,_T(" ")));
166 sizer->Add(new wxStaticText(panel,-1,_T(" - - - PATHOLOGY EXTRACTION - - - ")));
167 sizer->Add(new wxStaticText(panel,-1,_T(" ")));
168 wxButton *btnP1 = new wxButton(panel,-1,_T("Set P1"));
169 wxButton *btnP2 = new wxButton(panel,-1,_T("Set P2"));
170 wxButton *btnPat = new wxButton(panel,-1,_T("Set Region"));
171 wxButton *btnExtract = new wxButton(panel,-1,_T("Extract Region"));
172 wxButton *btnAxis = new wxButton(panel,-1, _T("Axis"));
175 sizer->Add(btnPat);sizer->Add(new wxStaticText(panel,-1,_T(" ")));
176 sizer->Add(btnExtract);sizer->Add(new wxStaticText(panel,-1,_T(" ")));
178 patSliderOpacity = new wxSlider( panel, -1, 50 , 0, 100, wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS );
179 sizer->Add(new wxStaticText(panel, -1,_T("Opacity Pathological VOI")));
180 sizer->Add(patSliderOpacity);
181 patSliderMarchingCubes = new wxSlider( panel, -1, 254 , 0, 1000, wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS );
182 sizer->Add(new wxStaticText(panel, -1,_T(" Marching Cubes Level Pathological VOI")));
183 sizer->Add(patSliderMarchingCubes);
184 sizer->Add(btnAxis);sizer->Add(new wxStaticText(panel,-1,_T(" ")));
185 //-- PATHOLOGY WIDGETS
191 sizer->Add(new wxStaticText(panel,-1,_T(" ")));
192 sizer->Add(new wxStaticText(panel,-1,_T(" ")));
193 sizer->Add(new wxStaticText(panel,-1,_T(" - - - STL Diego - - - ")));
194 sizer->Add(new wxStaticText(panel,-1,_T(" ")));
195 stlSliderDeltaGauss = new wxSlider( panel, -1, 100 , 0, 300 , wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS );
196 stlSliderMarchingCubes= new wxSlider( panel, -1, 128 , 0, 256 , wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS );
198 sizer->Add(new wxStaticText(panel,-1,_T(" Delta Gauss")));
199 sizer->Add(stlSliderDeltaGauss);
201 sizer->Add(new wxStaticText(panel,-1,_T(" Marching Cubes Level")));
202 sizer->Add(stlSliderMarchingCubes);
204 stlSliderOpacityInternal = new wxSlider(panel, -1, 100,0,100, wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS);
205 stlSliderOpacityExternal = new wxSlider(panel, -1, 100,0,100, wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS);
207 sizer->Add(new wxStaticText(panel, -1,_T(" Opacity STL Internal")));
208 sizer->Add(stlSliderOpacityInternal);
210 sizer->Add(new wxStaticText(panel, -1,_T(" Opacity STL External")));
211 sizer->Add(stlSliderOpacityExternal);
213 wxButton *btnFileSTL = new wxButton(panel,-1,_T("Generate STL files"));
214 sizer->Add(btnFileSTL);
215 sizer->Add(new wxStaticText(panel,-1,_T(" ")));
226 panel->SetSizer(sizer);
227 panel->SetAutoLayout(true);
228 panel->SetSize(400,600);
232 // -- PATHOLOGY CONNECT WIDGETS
233 Connect(btnP1->GetId(),wxEVT_COMMAND_BUTTON_CLICKED, (wxObjectEventFunction) &wxPathologyWidget_01::OnBtnSetP1);
234 Connect(btnP2->GetId(),wxEVT_COMMAND_BUTTON_CLICKED, (wxObjectEventFunction) &wxPathologyWidget_01::OnBtnSetP2);
235 Connect(btnPat->GetId(),wxEVT_COMMAND_BUTTON_CLICKED,(wxObjectEventFunction) &wxPathologyWidget_01::OnBtnSetPat);
236 Connect(btnExtract->GetId(),wxEVT_COMMAND_BUTTON_CLICKED,(wxObjectEventFunction) &wxPathologyWidget_01::OnBtnExtractPat);
237 Connect(patSliderOpacity->GetId(), wxEVT_COMMAND_SLIDER_UPDATED , (wxObjectEventFunction) &wxPathologyWidget_01::OnChangePatOpacity);
238 Connect(patSliderMarchingCubes->GetId(), wxEVT_COMMAND_SLIDER_UPDATED , (wxObjectEventFunction) &wxPathologyWidget_01::OnChangePatMarchingCubes);
239 Connect(btnAxis->GetId(), wxEVT_COMMAND_BUTTON_CLICKED , (wxObjectEventFunction) &wxPathologyWidget_01::OnBtnExtractAxis);
242 // -- PATHOLOGY CONNECT WIDGETS
245 // -- STL CONNECT WIDGETS
246 Connect(btnFileSTL->GetId() , wxEVT_COMMAND_BUTTON_CLICKED , (wxObjectEventFunction) &wxPathologyWidget_01::OnBtnFileSTL );
247 Connect(stlSliderDeltaGauss->GetId() , wxEVT_SCROLL_THUMBRELEASE , (wxObjectEventFunction) &wxPathologyWidget_01::OnChangeSTLGaussLevel );
248 Connect(stlSliderMarchingCubes->GetId() , wxEVT_SCROLL_THUMBRELEASE , (wxObjectEventFunction) &wxPathologyWidget_01::OnChangeSTLMarchingCubesLevel);
249 Connect(stlSliderOpacityInternal->GetId(), wxEVT_COMMAND_SLIDER_UPDATED , (wxObjectEventFunction) &wxPathologyWidget_01::OnOpacitySTLInternal);
250 Connect(stlSliderOpacityExternal->GetId(), wxEVT_COMMAND_SLIDER_UPDATED , (wxObjectEventFunction) &wxPathologyWidget_01::OnOpacitySTLExternal);
251 // -- STL CONNECT WIDGETS
259 //------------------------------------------------------------------------
260 void wxPathologyWidget_01::Refresh()
262 _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->Render();
265 //------------------------------------------------------------------------
266 void wxPathologyWidget_01::ConfigureVTK()
270 vtkImageData *imagedata = _mar->_experiment->getDynData( )->getVolume( )->castVtk();
272 _maracasSurfaceWidget->ShowMARACASData( _mar );
274 //CONFIGURACION ADICIONAL
275 this->ConfigureSTL();
276 this->ConfigurePathologyExtraction();
281 void wxPathologyWidget_01::ConfigureVTK(vtkImageData *imagedata, int x, int y, int z, double param){
288 _maracasSurfaceWidget->ShowMARACASData( _mar );
291 //CONFIGURACION ADICIONAL
292 this->ConfigureMPR();
293 this->ConfigureSTL();
294 this->ConfigurePathologyExtraction();
301 // ==================================================================================================
302 // CODIGO NUEVO ADICIONADO POR DIEGO CANTOR
303 // ==================================================================================================
306 // ------------------------------------------------------------------------
307 // START MPR FUNCTIONS - DHC
308 // ------------------------------------------------------------------------
310 void wxPathologyWidget_01::ConfigureMPR(){
312 mprSphereSource = vtkSphereSource::New();
313 mprMapper = vtkPolyDataMapper::New();
314 mprActor = vtkActor::New();
316 mprMapper->SetInput(mprSphereSource->GetOutput());
317 mprActor->SetMapper(mprMapper);
318 mprActor->PickableOff( );
320 mprSphereSource->SetCenter(this->px,this->py,this->pz);
321 mprSphereSource->SetRadius(2.5);
323 mprActor->GetProperty()->SetColor( 1.0, 0.0, 0.0 );
324 vtkRenderer *ren = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();
325 ren->AddActor(mprActor);
329 // ------------------------------------------------------------------------
330 // END MPR FUNCTIONS - DHC
331 // ------------------------------------------------------------------------
335 // ------------------------------------------------------------------------
336 // START STL FUNCTIONS - DHC
337 // ------------------------------------------------------------------------
339 void wxPathologyWidget_01::ConfigureSTL()
341 stlExterna = vtkPolyData::New();
342 stlInterna = vtkPolyData::New();
344 dsm1 = vtkPolyDataMapper ::New();
345 dsm1->SetInput (stlInterna);
346 dsm1->ScalarVisibilityOff();
348 actorInternal = vtkActor::New();
349 actorInternal->SetMapper (dsm1);
350 actorInternal->GetProperty()->SetColor (0,1,0);
352 dsm2 = vtkPolyDataMapper ::New();
353 dsm2->SetInput (stlExterna);
354 dsm2->ScalarVisibilityOff();
356 actorExternal= vtkActor::New();
357 actorExternal->SetMapper (dsm2);
358 actorExternal->GetProperty()->SetRepresentationToWireframe();
359 vtkRenderer *ren = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();
360 ren->AddActor(actorInternal);
361 ren->AddActor(actorExternal);
363 stlExtractor = new vtkSTLExtractor();
366 void wxPathologyWidget_01::generateSTLSurfaces()
368 stlExtractor->setVolume(stlImageData);
369 stlExtractor->setSigmaLevel(stlDeltaGaussLevel);
370 stlExtractor->setMarchingCubesLevel(stlMarchingCubesLevel);
371 stlExtractor->calculate();
373 stlInterna->DeepCopy(stlExtractor->getInnerSurface());
374 stlExterna->DeepCopy(stlExtractor->getOuterSurface());
378 void wxPathologyWidget_01::OnOpacitySTLExternal(wxScrollEvent& event){
379 double value = ((double)stlSliderOpacityExternal->GetValue())/100;
380 actorExternal->GetProperty( )->SetOpacity( value );
385 void wxPathologyWidget_01::OnOpacitySTLInternal(wxScrollEvent& event){
386 double value = ((double)stlSliderOpacityInternal->GetValue())/100;
387 actorInternal->GetProperty( )->SetOpacity( value );
391 void wxPathologyWidget_01::OnBtnFileSTL(wxCommandEvent& event)
394 wxString dirSTL = _mar->_parameters->getStringParam(
395 marParameters::e_installation_directory );
396 dirSTL = ( dirSTL == _T("NO_DIRECTORY") ) ? wxGetHomeDir( ) : dirSTL;
397 wxDirDialog dialog( this, _T("Choose a directory..."), ( !dirSTL.IsEmpty( ) )?
398 dirSTL: wxGetHomeDir( ) );
400 if( dialog.ShowModal( ) == wxID_OK )
404 // ------------------------------------------------------------------------
405 // 1. GENERATE STL FILES
406 // ------------------------------------------------------------------------
407 const char* fileprefix = "c:\\Creatis\\";
408 std::string prefix = fileprefix;
411 // 1.1. Se hace un filtro triangular puesto que el stl writer solo recibe poligonos triangulares.
413 vtkTriangleFilter *filtro = vtkTriangleFilter::New();
414 filtro->SetInput(stlInterna);
415 vtkPolyDataConnectivityFilter *pdcf = vtkPolyDataConnectivityFilter::New();
416 pdcf->SetInput( filtro->GetOutput() );
417 vtkClosePolyData *cpd = vtkClosePolyData::New();
418 cpd->SetInput( pdcf->GetOutput() );
420 // 1.2 se escribe a disco el archivo stl de la superficie interna
422 vtkSTLWriter *writerI = vtkSTLWriter::New();
423 writerI->SetInput( cpd->GetOutput() );
425 prefix += "internal.stl";
426 writerI->SetFileName(prefix.c_str());
427 writerI->SetFileTypeToASCII();
431 // 1.3 se escribe a disco el archivo stl de la superficie externa
432 filtro->SetInput(stlExterna);
434 vtkSTLWriter *writerE = vtkSTLWriter::New();
435 writerE->SetInput( cpd->GetOutput() );
437 prefix += "external.stl";
438 writerE->SetFileName( prefix.c_str() );
439 writerE->SetFileTypeToASCII();
448 //By default *always* update e_installation_directory:
449 _mar->_parameters->setStringParam( marParameters::e_installation_directory, dialog.GetPath( ) );
450 _mar->saveParameters( );
454 void wxPathologyWidget_01::OnChangeSTLGaussLevel(wxScrollEvent& event)
456 stlDeltaGaussLevel = ((double)stlSliderDeltaGauss->GetValue())/100;
457 generateSTLSurfaces();
462 void wxPathologyWidget_01::OnChangeSTLMarchingCubesLevel(wxScrollEvent& event)
464 stlMarchingCubesLevel = ((double)stlSliderMarchingCubes->GetValue());
465 generateSTLSurfaces();
471 // ------------------------------------------------------------------------
472 // END STL FUNCTIONS - DHC
473 // ------------------------------------------------------------------------
476 // ------------------------------------------------------------------------
477 // START PATHOLOGY FUNCTIONS - DHC
478 // ------------------------------------------------------------------------
480 double GetDistanciaEuclideana(double* a, double* b){
481 return sqrt(((a[0]-b[0])*(a[0]-b[0]))+((a[1]-b[1])*(a[1]-b[1]))+((a[2]-b[2])*(a[2]-b[2])));
485 double GetValorPromedio(vtkImageData *data){
493 data->GetExtent(ext);
496 for(i=ext[0];i<=ext[1];i++){
497 for(j=ext[2];j<=ext[3];j++){
498 for(k=ext[4];k<=ext[5];k++){
499 valor =data->GetScalarComponentAsDouble(i,j,k,0);
505 return (double)promedio/(double)numVoxels;
509 void wxPathologyWidget_01::ConfigurePathologyExtraction()
513 void wxPathologyWidget_01::generatePathologySurface(){
518 void wxPathologyWidget_01::OnBtnSetP1(wxCommandEvent &event){
519 vtk3DSurfaceWidget *sw = _maracasSurfaceWidget->GetVtk3DSurfaceWidget();
522 sw->GetRenderer()->RemoveActor(p1Actor);
526 sw->GetSphereCenter(p1Center);
528 p1SphereSource = vtkSphereSource::New();
529 p1Mapper = vtkPolyDataMapper::New();
530 p1Actor = vtkActor::New();
532 p1Mapper->SetInput(p1SphereSource->GetOutput());
533 p1Actor->SetMapper(p1Mapper);
534 p1Actor->PickableOff( );
536 p1SphereSource->SetCenter(p1Center[0],p1Center[1],p1Center[2]);
537 p1SphereSource->SetRadius(0.5);
539 p1Actor->GetProperty()->SetColor( 0.0, 1.0, 0.0);
541 sw->GetRenderer()->AddActor(p1Actor);
546 void wxPathologyWidget_01::OnBtnSetP2(wxCommandEvent &event){
548 vtk3DSurfaceWidget *sw = _maracasSurfaceWidget->GetVtk3DSurfaceWidget();
552 sw->GetRenderer()->RemoveActor(p1Actor);
555 sw->GetSphereCenter(p2Center);
557 p2SphereSource = vtkSphereSource::New();
558 p2Mapper = vtkPolyDataMapper::New();
559 p2Actor = vtkActor::New();
561 p2Mapper->SetInput(p2SphereSource->GetOutput());
562 p2Actor->SetMapper(p2Mapper);
563 p2Actor->PickableOff( );
565 p2SphereSource->SetCenter(p2Center[0],p2Center[1],p2Center[2]);
566 p2SphereSource->SetRadius(0.5);
568 p2Actor->GetProperty()->SetColor( 0.0, 1.0, 0.0);
570 sw->GetRenderer()->AddActor(p2Actor);
577 void wxPathologyWidget_01::OnBtnSetPat(wxCommandEvent &event){
579 if (p2SphereSource == NULL || p1SphereSource == NULL){
583 vtk3DSurfaceWidget *sw = _maracasSurfaceWidget->GetVtk3DSurfaceWidget();
586 sw->GetRenderer()->RemoveActor(patActor);
593 centro[0] = (p1Center[0] + p2Center[0]) / 2;
594 centro[1] = (p1Center[1] + p2Center[1]) / 2;
595 centro[2] = (p1Center[2] + p2Center[2]) / 2;
597 patSphereSource = vtkSphereSource::New();
598 patMapper = vtkPolyDataMapper::New();
599 patActor = vtkActor::New();
601 patMapper->SetInput(patSphereSource->GetOutput());
602 patActor->SetMapper(patMapper);
605 patSphereSource->SetCenter(centro[0],centro[1],centro[2]);
606 patSphereSource->SetRadius(GetDistanciaEuclideana(p1Center,p2Center));
608 patActor->GetProperty()->SetColor( 0.15, 0.75, 0.15 );
609 patActor->PickableOff( );
610 patActor->GetProperty( )->SetOpacity(0.15);
612 sw->GetRenderer()->AddActor(patActor);
619 void wxPathologyWidget_01::OnBtnExtractPat(wxCommandEvent &event){
622 //----------------------------------------------------------------------------------
623 vtkRenderer* renderer = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();
626 //----------------------------------------------------------------------------------
627 vtkImageData* originalData = vtkImageData::New();
628 originalData->DeepCopy(_mar->_experiment->getDynData()->getVolume()->castVtk());
631 originalData->GetSpacing(spacing);
634 originalData->GetExtent(extend);
638 //==================================================================================
639 // REGION DE INTERES EN LA IMAGEN ORIGINAL: patImageData
640 //==================================================================================
642 double p1crem[3], p2crem[3];
644 p1crem[0] = p1Center[0] / spacing[0];
645 p1crem[1] = p1Center[1] / spacing[1];
646 p1crem[2] = p1Center[2] / spacing[2];
648 p2crem[0] = p2Center[0] / spacing[0];
649 p2crem[1] = p2Center[1] / spacing[1];
650 p2crem[2] = p2Center[2] / spacing[2];
656 //COORDENADAS DE LA IMAGEN
658 double radio = GetDistanciaEuclideana(p1crem,p2crem) / 2;
661 centro[0] = (p1crem[0] + p2crem[0]) / 2;
662 centro[1] = (p1crem[1] + p2crem[1]) / 2;
663 centro[2] = (p1crem[2] + p2crem[2]) / 2;
667 //Obtener los ejes sobre la esfera para delimitar el area.
669 double lowX = centro[0] - (2*radio);
670 double highX = centro[0] + (2*radio);
671 double lowY = centro[1] - (2*radio);
672 double highY = centro[1] + (2*radio);
673 double lowZ = centro[2] - (2*radio);
674 double highZ = centro[2] + (2*radio);
678 //Extraer la region de interes:patImageData
679 //----------------------------------------------------------------------------------
681 voiFilter = vtkExtractVOI::New();
682 voiFilter->SetInput(originalData);
683 voiFilter->SetVOI((int)lowX, (int)highX, (int)lowY, (int)highY, (int)lowZ, (int)highZ);
685 patImageData = vtkImageData::New();
686 patImageData->DeepCopy(voiFilter->GetOutput());
687 // patImageData->SetSpacing(spacing);
689 //==================================================================================
690 // REGION DE INTERES ESPECIFICA
691 //==================================================================================
693 vtkOtsuSphereSource *otsu = new vtkOtsuSphereSource();
694 thresholdOTSU = otsu->calculateOptimalThreshold(patImageData);
697 //==================================================================================
699 //==================================================================================
702 cubesFilter = vtkMarchingCubes::New();
703 cubesFilter->SetInput(otsu->getImageForSegmentation());
704 cubesFilter->SetValue(0,patMCLevel);
705 cubesFilter->SetNumberOfContours(1);
707 //==================================================================================
709 //==================================================================================
711 polyUnico = vtkPolyDataConnectivityFilter::New();
712 polyUnico->SetInput(cubesFilter->GetOutput());
713 polyUnico->SetExtractionModeToLargestRegion();
715 stripper = vtkStripper::New();
716 stripper->SetInput(polyUnico->GetOutput());
719 //==================================================================================
721 //==================================================================================
723 dataMapper = vtkPolyDataMapper::New();
724 dataMapper->SetInput(stripper->GetOutput());
725 dataMapper->ScalarVisibilityOff();
727 //==================================================================================
729 //==================================================================================
732 dataActor = vtkActor::New();
733 dataActor->SetMapper(dataMapper);
734 dataActor->SetPosition(lowX*spacing[0],lowY*spacing[1],lowZ*spacing[2]);
735 dataActor->GetProperty()->SetOpacity(patOpacityLevel);
736 dataActor->GetProperty()->SetColor( 0.20, 0.20, 1.0 );
737 dataActor->GetProperty()->SetRepresentationToWireframe();
739 renderer->AddActor(dataActor);
743 outlineFilter = vtkOutlineFilter::New();
744 outlineFilter->SetInput(voiFilter->GetOutput());
746 vtkPolyDataMapper* outlineMapper = vtkPolyDataMapper::New();
747 outlineMapper->SetInput(outlineFilter->GetOutput());
750 renderer->RemoveActor(outlineActor);
754 renderer->RemoveActor(patActor);
758 renderer->RemoveActor(p3Actor);
762 outlineActor = vtkActor::New();
763 outlineActor->SetMapper(outlineMapper);
764 outlineActor->GetProperty()->SetOpacity(1);
767 /*p3SphereSource = vtkSphereSource::New();
768 p3Mapper = vtkPolyDataMapper::New();
769 p3Actor = vtkActor::New();
771 p3Mapper->SetInput(p3SphereSource->GetOutput());
772 p3Actor->SetMapper(p3Mapper);
773 p3Actor->PickableOff( );
775 p3SphereSource->SetCenter(p2Center[0]-(lowX*spacing[0]), p2Center[1]-(lowY*spacing[1]) , p2Center[2]-(lowZ*spacing[2]) );
776 p3SphereSource->SetRadius(0.5);
777 p3Actor->GetProperty()->SetColor( 0.9, 0.4, 0.2);*/
780 //renderer->AddActor(outlineActor);
781 //renderer->AddActor(p3Actor);
788 void wxPathologyWidget_01::OnChangePatOpacity(wxScrollEvent& event){
789 patOpacityLevel = ((double)patSliderOpacity->GetValue())/100.0;
790 dataActor->GetProperty( )->SetOpacity(patOpacityLevel);
794 void wxPathologyWidget_01::OnChangePatMarchingCubes(wxScrollEvent& event){
795 patMCLevel = ((double)patSliderMarchingCubes->GetValue());
796 cubesFilter->SetValue(0,patMCLevel);
797 cubesFilter->Update();
803 //Metodo auxiliar que devuelve la distancia minima de un punto a un polydata.
804 //Se evalua la distancia del punto a todos los centros de las celdas del polydata y se devuelve la menor
805 double GetDistanciaAPolydata(double* punto, vtkPoints* centers){
806 double distancia = 0.0;
807 double minima = 1e100;
809 int N = centers->GetNumberOfPoints();
811 for(int i = 0; i < N; i++){
812 centers->GetPoint(i,vertice);
813 distancia = sqrt(((punto[0]-vertice[0])*(punto[0]-vertice[0]))+((punto[1]-vertice[1])*(punto[1]-vertice[1]))+((punto[2]-vertice[2])*(punto[2]-vertice[2])));
815 if (distancia < minima){
823 double GetProductoPunto(double* vectorA, double* vectorB){
824 return (vectorA[0]*vectorB[0] + vectorA[1]*vectorB[1] + vectorA[2]*vectorB[2]);
828 double GetNormaVector(double* vector){
829 return sqrt(GetProductoPunto(vector, vector));
832 //construye el vector como la diferencia de dos puntos (vectores en el origen)
833 void GetVector(double* p0, double* p1, double* resultado){ //p1:destino - p0:origen
835 resultado[0] = p1[0] - p0[0];
836 resultado[1] = p1[1] - p0[1];
837 resultado[2] = p1[2] - p0[2];
842 void VectorUnitario(double *v){
843 double norma = GetNormaVector(v);
849 //calculo de las normales
850 void GetNormalCelda(double* p0, double* p1, double* p2, double* normal){
853 GetVector(p1,p2,U); //p2-p1
854 GetVector(p1,p0,V); //p0-p1
857 normal[0] = (U[1]*V[2])-(U[2]*V[1]);
858 normal[1] = -((U[0]*V[2])-(U[2]*V[0]));
859 normal[2] = (U[0]*V[1])-(U[1]*V[0]);
860 VectorUnitario(normal);
863 void GetNormal(vtkPoints* puntosCelda, double* normal){
864 double p0[3], p1[3], p2[3];
866 puntosCelda->GetPoint(0, p0);
867 puntosCelda->GetPoint(1, p1);
868 puntosCelda->GetPoint(2, p2);
870 GetNormalCelda(p0,p1,p2, normal);
874 // Obtiene el ID de la celda mas cercana
875 int GetClosestCellOnPolyData(double* punto, vtkPoints* centers){
876 double distancia = 0.0;
877 double minima = 1e100;
880 int N = centers->GetNumberOfPoints();
882 for(int i = 0; i < N; i++){
883 centers->GetPoint(i,vertice);
884 distancia = sqrt(((punto[0]-vertice[0])*(punto[0]-vertice[0]))+((punto[1]-vertice[1])*(punto[1]-vertice[1]))+((punto[2]-vertice[2])*(punto[2]-vertice[2])));
886 if (distancia < minima){
892 //printf("distancia del punto p: (%.2f,%.2f,%.2f) al polydata es %.2f (ID %i)\n", punto[0], punto[1], punto [2], minima, ID);
896 double GetProjectionOnNormal(double* vector, double* normal){
897 double pp = GetProductoPunto(vector, normal);
898 double vv = GetNormaVector(normal);
907 void wxPathologyWidget_01::OnBtnExtractAxis(wxCommandEvent& event){
911 //================================================================================
913 //================================================================================
914 vtkRenderer* renderer = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();
916 //===================================================================================
918 //===================================================================================
919 int extentOriginal[6];
920 this->patImageData->GetExtent(extentOriginal);
923 //===================================================================================
925 //===================================================================================
928 double resampling = 1.0/(double)factor;
932 // REMUESTREO DE LA IMAGEN ORIGINAL
933 //-----------------------------------------------------------------------------------
934 vtkImageResample* resamplingFilter = vtkImageResample::New();
935 resamplingFilter->SetAxisOutputSpacing(0,resampling);
936 resamplingFilter->SetAxisOutputSpacing(1,resampling);
937 resamplingFilter->SetAxisOutputSpacing(2,resampling);
938 resamplingFilter->ReleaseDataFlagOff();
939 resamplingFilter->SetInput(patImageData);
940 resamplingFilter->SetInterpolationModeToCubic();
941 resamplingFilter->Update();
944 //IMAGEN ORIGINAL REMUESTREADA
945 //-----------------------------------------------------------------------------------
946 vtkImageData* imagenRemuestreada = resamplingFilter->GetOutput();
947 imagenRemuestreada->ReleaseDataFlagOff();
951 //===================================================================================
952 // IMAGEN BINARIA = CRECIMIENTO DE REGIONES
953 //===================================================================================
955 this->isoBinaria = vtkImageData::New();
956 this->isoBinaria->CopyStructure(resamplingFilter->GetOutput());
957 this->isoBinaria->SetScalarTypeToDouble();
961 //CALCULO DE LAS DIMENSIONES Y EL EXTENT
962 //-----------------------------------------------------------------------------------
963 int extentBinaria[6];
964 this->isoBinaria->GetExtent(extentBinaria);
965 int *dim = isoBinaria->GetDimensions();
970 //-----------------------------------------------------------------------------------
971 int eX, eY, eZ, idpoint; // extent
974 FILE* logger = freopen("tests.txt","w",stdout);
978 for(eX = extentBinaria[0]; eX <= extentBinaria[1]; eX++){
979 for(eY = extentBinaria[2]; eY <= extentBinaria[3]; eY++){
980 for (eZ = extentBinaria[4]; eZ <= extentBinaria[5]; eZ++){
981 idpoint = isoBinaria->FindPoint(eX, eY, eZ);
982 ptr = (double *)isoBinaria->GetScalarPointer(eX,eY,eZ); // El apuntador se obtiene sobre el extent.
989 //===================================================================================
991 //===================================================================================
992 vtkPolyData *polyData= stripper->GetOutput(); //OJO SE TOMA SOLO LA SUPERFICIE MAS GRANDE (connectivity + stripper)
993 //===================================================================================
995 vtkPoints *puntosPolyData = polyData->GetPoints();
996 vtkPointData *pointData = polyData->GetPointData();
998 vtkPolyDataWriter* writePoly = vtkPolyDataWriter::New();
999 writePoly->SetFileName("poly.vtk");
1000 writePoly->SetInput(polyData);
1005 // PROCESAR LOS VERTICES
1006 //----------------------------------------------------------------------------------
1008 int numVertices = puntosPolyData->GetNumberOfPoints();
1010 // mover el sistema de referencias del polydata para poder determinar lo que se encuentra
1011 //por dentro segun las semillas.
1014 /*for(j=0;j<numVertices;j++){
1015 puntosPolyData->GetPoint(j,vertice);
1016 vertice[0] = vertice[0] + extentOriginal[0];
1017 vertice[1] = vertice[1] + extentOriginal[2];
1018 vertice[2] = vertice[2] + extentOriginal[4];
1019 puntosPolyData->SetPoint(j,vertice);
1024 // PROCESAR LAS CELDAS
1025 //----------------------------------------------------------------------------------
1028 //1. Obtener normales de la superficie (celdas)
1029 vtkPolyDataNormals* polyNormalsFilter = vtkPolyDataNormals::New(); //estas son las normales de los vertices NO de las celdas
1030 polyNormalsFilter->SetInput(polyData);
1031 polyNormalsFilter->SplittingOn();
1032 polyNormalsFilter->ConsistencyOn();
1033 polyNormalsFilter->NonManifoldTraversalOn();
1034 polyNormalsFilter->ComputeCellNormalsOn();
1035 polyNormalsFilter->Update();
1038 vtkPolyData* polyDataConNormales = polyNormalsFilter->GetOutput();
1041 //2. Obtener los centros de las celdas
1042 // Necesitamos el centro de las celdas para trazar el vector interno de un voxel
1043 //a una cara de la superficie.
1044 vtkCellCenters* cellCentersFilter = vtkCellCenters::New();
1045 cellCentersFilter->SetInput(polyDataConNormales);
1046 cellCentersFilter->VertexCellsOn();
1047 cellCentersFilter->Update();
1049 vtkPolyData* polydataCenters = cellCentersFilter->GetOutput();
1050 vtkPoints* polyCenterPoints = polydataCenters->GetPoints();
1052 vtkPoints* puntosCelda = NULL;
1053 vtkIdList* lista = NULL;
1054 vtkCell* celda = NULL;
1056 int numCeldas = polyCenterPoints->GetNumberOfPoints();
1059 vtkPoints *normales = vtkPoints::New();//almacena los componentes cartesianos de la normal de cada celda
1062 for(i = 0; i < numCeldas ; i++){
1064 polyCenterPoints->GetPoint(i,vertice);
1066 celda = polyDataConNormales->GetCell(i);
1067 puntosCelda = celda->GetPoints();
1068 GetNormal(puntosCelda, normal);
1069 normales->InsertPoint(i, normal);
1072 lista = celda->GetPointIds();
1073 for(j = 0; j < lista->GetNumberOfIds(); j++){
1074 int idPunto = lista->GetId(j);
1075 puntosCelda->GetPoint(j, vertice);
1083 //==========================================================================================
1084 //IMPLEMENTAR EL CRECIMIENTO DE REGIONES
1085 //==========================================================================================
1088 patImageData->GetSpacing(spacing);
1090 double p1crem[3], p2crem[3];
1092 p1crem[0] = p1Center[0] / spacing[0];
1093 p1crem[1] = p1Center[1] / spacing[1];
1094 p1crem[2] = p1Center[2] / spacing[2];
1096 p2crem[0] = p2Center[0] / spacing[0];
1097 p2crem[1] = p2Center[1] / spacing[1];
1098 p2crem[2] = p2Center[2] / spacing[2];
1104 //COORDENADAS DE LA IMAGEN
1106 double radio = GetDistanciaEuclideana(p1crem,p2crem) / 2;
1109 centro[0] = (p1crem[0] + p2crem[0]) / 2;
1110 centro[1] = (p1crem[1] + p2crem[1]) / 2;
1111 centro[2] = (p1crem[2] + p2crem[2]) / 2;
1115 //Obtener los ejes sobre la esfera para delimitar el area.
1117 double lowX = centro[0] - (2*radio);
1118 double highX = centro[0] + (2*radio);
1119 double lowY = centro[1] - (2*radio);
1120 double highY = centro[1] + (2*radio);
1121 double lowZ = centro[2] - (2*radio);
1122 double highZ = centro[2] + (2*radio);
1124 //================================================================================
1126 //================================================================================
1128 int numeroPuntos = dim[0] * dim[1] * dim[2];
1134 sem[0] = (int)(p2Center[0]-(lowX*spacing[0]));
1135 sem[1] = (int)(p2Center[1]-(lowY*spacing[1]));
1136 sem[2] = (int)(p2Center[2]-(lowZ*spacing[2]));
1139 int idSemilla = isoBinaria->FindPoint(sem[0],sem[1],sem[2]);
1141 if (idSemilla == -1) return;
1143 int graphSize = (extentBinaria[1]-extentBinaria[0]) * (extentBinaria[3]-extentBinaria[2])*(extentBinaria[5]-extentBinaria[4]);
1146 vtkPriorityQueue *colaEvaluacion = vtkPriorityQueue::New();
1148 colaEvaluacion->Allocate(graphSize);
1150 colaEvaluacion->Insert(0, idSemilla);
1152 vtkFloatingPointType prioridad;
1154 int longitudCola = colaEvaluacion->GetNumberOfItems();
1161 double* coordsVecino;
1163 double* coordsPunto;
1167 vtkPoints *resultado = vtkPoints::New();
1172 double centroCelda[3], vectorInterno[3], vectorNormal[3];
1174 while(longitudCola > 0){
1176 //================================================================================
1177 //1. Procesar el punto
1178 //================================================================================
1180 puntoID = colaEvaluacion->Pop(0,prioridad);
1182 coordsPunto = isoBinaria->GetPoint(puntoID);
1185 ptr = (double *)isoBinaria->GetScalarPointer((int)(coordsPunto[0]*factor), (int)(coordsPunto[1]*factor), (int)(coordsPunto[2]*factor));
1189 vtkImageWriter* writeImage = vtkImageWriter::New();
1190 writeImage->SetInput(isoBinaria);
1191 writeImage->SetFileName("sem.vtk");
1192 writeImage->Write();
1196 resultado->InsertPoint(cont, coordsPunto);
1199 //printf("iniciando en el punto (%.1f,%.1f,%.1f) \n", coordsPunto[0], coordsPunto[1], coordsPunto [2]);
1201 //================================================================================
1202 //2. Evaluar los vecinos del punto
1203 //================================================================================
1205 for(vk = -1; vk<2; vk++) {
1206 for(vj = -1; vj<2; vj++) {
1207 for(vi = -1; vi<2; vi++) {
1209 vecinoID = puntoID + (vk * dim[1]*dim[0]) + (vj * dim[0]) + vi;
1211 if( vecinoID >= 0 && vecinoID < numeroPuntos && vecinoID != 0 && vecinoID != puntoID ) {
1213 coordsVecino = isoBinaria->GetPoint(vecinoID);
1215 //ptr = (unsigned char *)isoBinaria->GetScalarPointer(coordsVecino[0],coordsVecino[1],coordsVecino[2]);
1216 ptr = (double *)isoBinaria->GetScalarPointer((int)(coordsVecino[0]*factor),(int)(coordsVecino[1]*factor),(int)(coordsVecino[2]*factor));
1218 //================================================================================
1219 // Para cada vecino no evaluado evaluar la distancia al polydata
1220 //================================================================================
1222 //printf("evaluando punto (%.1f,%.1f,%.1f) \n", coordsVecino[0], coordsVecino[1], coordsVecino [2]);
1224 int normalID = GetClosestCellOnPolyData(coordsVecino, polyCenterPoints);
1227 polyCenterPoints->GetPoint(normalID,centroCelda);
1228 normales->GetPoint(normalID, vectorNormal);
1230 GetVector(coordsVecino, centroCelda, vectorInterno);
1231 VectorUnitario(vectorInterno);
1233 double projection = GetProjectionOnNormal(vectorInterno, vectorNormal);
1234 ptr = (double *)isoBinaria->GetScalarPointer((int)(coordsVecino[0]*factor),(int)(coordsVecino[1]*factor),(int)(coordsVecino[2]*factor));
1236 if (projection > 0){
1238 //================================================================================
1239 //2.1 Si la proyeccion es positiva
1240 //colocar el vecino en la cola y poner el punto correspondiente en la imagen a 255.
1241 // (no ha pasado la frontera, el marcado a 255 se hace arriba en el while, ver paso 1)
1242 //================================================================================
1243 colaEvaluacion->Insert(prioridad,vecinoID);
1247 //================================================================================
1248 //2.2 Si la proyeccion es negativa con la normal mas cercana del polydata,
1249 //descartar el vecino == no poner en la cola == alcanzamos la superficie
1250 //================================================================================
1261 longitudCola = colaEvaluacion->GetNumberOfItems();
1263 colaEvaluacion->Delete();
1266 //================================================================================
1267 // MAPA DE DISTANCIAS
1268 //================================================================================
1270 this->isoDistance = vtkImageEuclideanDistance::New();
1271 this->isoDistance->SetInput(this->isoBinaria);
1272 this->isoDistance->InitializeOn();
1273 this->isoDistance->Update();
1274 vtkImageData *imDistancias = isoDistance->GetOutput();
1277 //================================================================================
1278 // CAMINO MINIMO SOBRE EL MAPA DE DISTANCIAS
1279 //================================================================================
1280 this->dijkstraFilter = vtkDijkstraImageData::New();
1281 this->dijkstraFilter->SetInput(isoDistance->GetOutput());
1283 int idA = (int)isoBinaria->FindPoint(p1Center[0],p1Center[1], p1Center[2]);
1284 int idB = (int)isoBinaria->FindPoint(p2Center[0],p2Center[1], p2Center[2]);
1286 dijkstraFilter->SourceID = idA;
1287 dijkstraFilter->SinkID = idB;
1288 dijkstraFilter->Update();
1290 caminoMapper = vtkPolyDataMapper::New();
1291 caminoMapper->SetInput(dijkstraFilter->GetOutput());
1292 caminoMapper->ScalarVisibilityOff();
1294 renderer->RemoveActor(dataActor);
1296 //=======================================================================================================
1297 //RECORRIDO DE LOS PUNTOS DEL EJE DETECTADO
1298 //=======================================================================================================
1299 vtkIdList *camino = dijkstraFilter->GetShortestPathIdList();
1301 int size = camino->GetNumberOfIds();
1305 for(i=0;i < size;i++){
1306 int pointId = camino->GetId(i);
1307 isoBinaria->GetPoint(pointId, coords);
1308 double scalar = imagenRemuestreada->GetPointData()->GetScalars()->GetTuple1(pointId);
1310 //----------------------------------------------------------------------------------------
1311 //Visualizacion del nodo del eje
1312 //----------------------------------------------------------------------------------------
1313 vtkSphereSource *_nodo = vtkSphereSource::New();
1314 _nodo->SetCenter(coords[0],coords[1],coords[2]);
1315 _nodo->SetRadius(0.2);
1317 vtkPolyDataMapper *_mapper = vtkPolyDataMapper::New();
1318 _mapper->SetInput(_nodo->GetOutput());
1321 vtkActor *_actor = vtkActor::New();
1322 _actor->GetProperty()->SetColor( 0.8, 0.2, 0.3);
1323 _actor->SetMapper(_mapper);
1324 _actor->PickableOff( );
1326 renderer->AddActor(_actor);
1333 // ------------------------------------------------------------------------
1334 // END PATHOLOGY FUNCTIONS - DHC
1335 // ------------------------------------------------------------------------