5 #include "vtkAppendFilter.h"
7 #include "vtkCallbackCommand.h"
8 #include "vtkCellArray.h"
9 #include "vtkCylinderSource.h"
10 #include "vtkClosePolyData.h"
11 #include "vtkCommand.h"
12 #include "vtkDoubleArray.h"
13 #include "vtkExtractVOI.h"
14 #include "vtkImageReader.h"
15 #include "vtkImageViewer.h"
16 #include "vtkImageViewer2.h"
17 #include "vtkImageToStructuredPoints.h"
18 #include "vtkImageThreshold.h"
19 #include "vtkImageWriter.h"
20 #include "vtkImageClip.h"
21 #include "vtkImageResample.h"
22 #include "vtkImageReslice.h"
23 #include "vtkImageResample.h"
24 #include "vtkImageThreshold.h"
25 #include "vtkImageCast.h"
26 #include "vtkImageData.h"
28 #include "vtkMarchingCubes.h"
29 #include "vtkObjectFactory.h"
30 #include "vtkPolyDataMapper.h"
31 #include "vtkRenderer.h"
32 #include "vtkRenderWindow.h"
33 #include "vtkRenderWindowInteractor.h"
34 #include "vtkPolyDataConnectivityFilter.h"
35 #include "vtkProperty.h"
36 #include "vtkPriorityQueue.h"
37 #include "vtkPoints.h"
38 #include "vtkPolyData.h"
39 #include "vtkPolyDataMapper.h"
40 #include "vtkPolyDataWriter.h"
41 #include "vtkPolyDataReader.h"
42 #include "vtkPointData.h"
43 #include "vtkSphereSource.h"
44 #include "vtkSTLWriter.h"
45 #include "vtkStripper.h"
46 #include "vtkTriangleFilter.h"
47 #include "vtkTransform.h"
49 #include "wxPathologyWidget_01.h"
50 #include "kernel/vtkOtsuSphereSource.h"
51 #include <wx/splitter.h>
56 //-------------------------------------------------------------------
57 //-------------------------------------------------------------------
58 //-------------------------------------------------------------------
59 wxPathologyWidget_01::wxPathologyWidget_01(wxWindow *parent, marInterface* mar)
60 : wxPanel( parent, -1)
65 wxBoxSizer *sizer = new wxBoxSizer(wxVERTICAL );
66 wxSplitterWindow *pnlSplitter = new wxSplitterWindow( this , -1);
67 wxPanel *viewPanel = CreateViewPanel(pnlSplitter);
68 wxPanel *controlPanel = CreateControlPanel(pnlSplitter);
70 sizer -> Add( pnlSplitter ,1,wxGROW ,0);
71 pnlSplitter -> SetMinimumPaneSize( 150 );
72 pnlSplitter -> SplitVertically( viewPanel, controlPanel, 600 );
73 this -> SetSizer(sizer);
84 p1SphereSource = NULL;
88 p2SphereSource = NULL;
92 patSphereSource = NULL;
97 pathologyFrame = NULL;
105 outlineFilter = NULL;
106 outlineMapper = NULL;
110 p3SphereSource = NULL;
116 //-------------------------------------------------------------------
117 wxPathologyWidget_01::~wxPathologyWidget_01(){
120 //-------------------------------------------------------------------
121 wxPanel* wxPathologyWidget_01::CreateViewPanel(wxWindow *parent)
123 wxPanel *panel = new wxPanel(parent,-1);
124 wxBoxSizer *sizer = new wxBoxSizer(wxVERTICAL);
125 _maracasSurfaceWidget = new wxSurfaceWidget(panel);
126 sizer->Add(_maracasSurfaceWidget , 1, wxEXPAND, 0);
127 panel->SetSizer(sizer);
128 panel->SetAutoLayout(true);
129 panel->SetSize(400,400);
133 //-------------------------------------------------------------------
134 wxPanel* wxPathologyWidget_01::CreateControlPanel(wxWindow *parent)
136 wxPanel *panel = new wxPanel(parent,-1);
137 wxFlexGridSizer *sizer = new wxFlexGridSizer(2);
140 //-- PATHOLOGY WIDGETS
142 sizer->Add(new wxStaticText(panel,-1,_T(" ")));
143 sizer->Add(new wxStaticText(panel,-1,_T(" ")));
144 sizer->Add(new wxStaticText(panel,-1,_T(" - - - PATHOLOGY EXTRACTION - - - ")));
145 sizer->Add(new wxStaticText(panel,-1,_T(" ")));
146 wxButton *btnP1 = new wxButton(panel,-1,_T("Set P1"));
147 wxButton *btnP2 = new wxButton(panel,-1,_T("Set P2"));
148 wxButton *btnPat = new wxButton(panel,-1,_T("Set Region"));
149 wxButton *btnExtract = new wxButton(panel,-1,_T("Extract Region"));
150 wxButton *btnAxis = new wxButton(panel,-1, _T("Axis"));
153 sizer->Add(btnPat);sizer->Add(new wxStaticText(panel,-1,_T(" ")));
154 sizer->Add(btnExtract);sizer->Add(new wxStaticText(panel,-1,_T(" ")));
156 patSliderOpacity = new wxSlider( panel, -1, 50 , 0, 100, wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS );
157 sizer->Add(new wxStaticText(panel, -1,_T("Opacity Pathological VOI")));
158 sizer->Add(patSliderOpacity);
159 patSliderMarchingCubes = new wxSlider( panel, -1, 254 , 0, 1000, wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS );
160 sizer->Add(new wxStaticText(panel, -1,_T(" Marching Cubes Level Pathological VOI")));
161 sizer->Add(patSliderMarchingCubes);
162 sizer->Add(btnAxis);sizer->Add(new wxStaticText(panel,-1,_T(" ")));
163 //-- PATHOLOGY WIDGETS
169 sizer->Add(new wxStaticText(panel,-1,_T(" ")));
170 sizer->Add(new wxStaticText(panel,-1,_T(" ")));
171 sizer->Add(new wxStaticText(panel,-1,_T(" - - - STL Diego - - - ")));
172 sizer->Add(new wxStaticText(panel,-1,_T(" ")));
173 stlSliderDeltaGauss = new wxSlider( panel, -1, 100 , 0, 300 , wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS );
174 stlSliderMarchingCubes= new wxSlider( panel, -1, 128 , 0, 256 , wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS );
176 sizer->Add(new wxStaticText(panel,-1,_T(" Delta Gauss")));
177 sizer->Add(stlSliderDeltaGauss);
179 sizer->Add(new wxStaticText(panel,-1,_T(" Marching Cubes Level")));
180 sizer->Add(stlSliderMarchingCubes);
182 stlSliderOpacityInternal = new wxSlider(panel, -1, 100,0,100, wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS);
183 stlSliderOpacityExternal = new wxSlider(panel, -1, 100,0,100, wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS);
185 sizer->Add(new wxStaticText(panel, -1,_T(" Opacity STL Internal")));
186 sizer->Add(stlSliderOpacityInternal);
188 sizer->Add(new wxStaticText(panel, -1,_T(" Opacity STL External")));
189 sizer->Add(stlSliderOpacityExternal);
191 wxButton *btnFileSTL = new wxButton(panel,-1,_T("Generate STL files"));
192 sizer->Add(btnFileSTL);
193 sizer->Add(new wxStaticText(panel,-1,_T(" ")));
204 panel->SetSizer(sizer);
205 panel->SetAutoLayout(true);
206 panel->SetSize(400,600);
210 // -- PATHOLOGY CONNECT WIDGETS
211 Connect(btnP1->GetId(),wxEVT_COMMAND_BUTTON_CLICKED, (wxObjectEventFunction) &wxPathologyWidget_01::OnBtnSetP1);
212 Connect(btnP2->GetId(),wxEVT_COMMAND_BUTTON_CLICKED, (wxObjectEventFunction) &wxPathologyWidget_01::OnBtnSetP2);
213 Connect(btnPat->GetId(),wxEVT_COMMAND_BUTTON_CLICKED,(wxObjectEventFunction) &wxPathologyWidget_01::OnBtnSetPat);
214 Connect(btnExtract->GetId(),wxEVT_COMMAND_BUTTON_CLICKED,(wxObjectEventFunction) &wxPathologyWidget_01::OnBtnExtractPat);
215 Connect(patSliderOpacity->GetId(), wxEVT_COMMAND_SLIDER_UPDATED , (wxObjectEventFunction) &wxPathologyWidget_01::OnChangePatOpacity);
216 Connect(patSliderMarchingCubes->GetId(), wxEVT_COMMAND_SLIDER_UPDATED , (wxObjectEventFunction) &wxPathologyWidget_01::OnChangePatMarchingCubes);
217 Connect(btnAxis->GetId(), wxEVT_COMMAND_BUTTON_CLICKED , (wxObjectEventFunction) &wxPathologyWidget_01::OnBtnExtractAxis);
220 // -- PATHOLOGY CONNECT WIDGETS
223 // -- STL CONNECT WIDGETS
224 Connect(btnFileSTL->GetId() , wxEVT_COMMAND_BUTTON_CLICKED , (wxObjectEventFunction) &wxPathologyWidget_01::OnBtnFileSTL );
225 Connect(stlSliderDeltaGauss->GetId() , wxEVT_SCROLL_THUMBRELEASE , (wxObjectEventFunction) &wxPathologyWidget_01::OnChangeSTLGaussLevel );
226 Connect(stlSliderMarchingCubes->GetId() , wxEVT_SCROLL_THUMBRELEASE , (wxObjectEventFunction) &wxPathologyWidget_01::OnChangeSTLMarchingCubesLevel);
227 Connect(stlSliderOpacityInternal->GetId(), wxEVT_COMMAND_SLIDER_UPDATED , (wxObjectEventFunction) &wxPathologyWidget_01::OnOpacitySTLInternal);
228 Connect(stlSliderOpacityExternal->GetId(), wxEVT_COMMAND_SLIDER_UPDATED , (wxObjectEventFunction) &wxPathologyWidget_01::OnOpacitySTLExternal);
229 // -- STL CONNECT WIDGETS
237 //------------------------------------------------------------------------
238 void wxPathologyWidget_01::Refresh()
240 _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->Render();
243 //------------------------------------------------------------------------
244 void wxPathologyWidget_01::ConfigureVTK()
248 vtkImageData *imagedata = _mar->_experiment->getDynData( )->getVolume( )->castVtk();
250 _maracasSurfaceWidget->ShowMARACASData( _mar );
252 //CONFIGURACION ADICIONAL
253 this->ConfigureSTL();
254 this->ConfigurePathologyExtraction();
259 void wxPathologyWidget_01::ConfigureVTK(vtkImageData *imagedata, int x, int y, int z, double param){
266 _maracasSurfaceWidget->ShowMARACASData( _mar );
269 //CONFIGURACION ADICIONAL
270 this->ConfigureMPR();
271 this->ConfigureSTL();
272 this->ConfigurePathologyExtraction();
279 // ==================================================================================================
280 // CODIGO NUEVO ADICIONADO POR DIEGO CANTOR
281 // ==================================================================================================
284 // ------------------------------------------------------------------------
285 // START MPR FUNCTIONS - DHC
286 // ------------------------------------------------------------------------
288 void wxPathologyWidget_01::ConfigureMPR(){
290 mprSphereSource = vtkSphereSource::New();
291 mprMapper = vtkPolyDataMapper::New();
292 mprActor = vtkActor::New();
294 mprMapper->SetInput(mprSphereSource->GetOutput());
295 mprActor->SetMapper(mprMapper);
296 mprActor->PickableOff( );
298 mprSphereSource->SetCenter(this->px,this->py,this->pz);
299 mprSphereSource->SetRadius(2.5);
301 mprActor->GetProperty()->SetColor( 1.0, 0.0, 0.0 );
302 vtkRenderer *ren = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();
303 ren->AddActor(mprActor);
307 // ------------------------------------------------------------------------
308 // END MPR FUNCTIONS - DHC
309 // ------------------------------------------------------------------------
313 // ------------------------------------------------------------------------
314 // START STL FUNCTIONS - DHC
315 // ------------------------------------------------------------------------
317 void wxPathologyWidget_01::ConfigureSTL()
319 stlExterna = vtkPolyData::New();
320 stlInterna = vtkPolyData::New();
322 dsm1 = vtkPolyDataMapper ::New();
323 dsm1->SetInput (stlInterna);
324 dsm1->ScalarVisibilityOff();
326 actorInternal = vtkActor::New();
327 actorInternal->SetMapper (dsm1);
328 actorInternal->GetProperty()->SetColor (0,1,0);
330 dsm2 = vtkPolyDataMapper ::New();
331 dsm2->SetInput (stlExterna);
332 dsm2->ScalarVisibilityOff();
334 actorExternal= vtkActor::New();
335 actorExternal->SetMapper (dsm2);
336 actorExternal->GetProperty()->SetRepresentationToWireframe();
337 vtkRenderer *ren = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();
338 ren->AddActor(actorInternal);
339 ren->AddActor(actorExternal);
341 stlExtractor = new vtkSTLExtractor();
344 void wxPathologyWidget_01::generateSTLSurfaces()
346 stlExtractor->setVolume(stlImageData);
347 stlExtractor->setSigmaLevel(stlDeltaGaussLevel);
348 stlExtractor->setMarchingCubesLevel(stlMarchingCubesLevel);
349 stlExtractor->calculate();
351 stlInterna->DeepCopy(stlExtractor->getInnerSurface());
352 stlExterna->DeepCopy(stlExtractor->getOuterSurface());
356 void wxPathologyWidget_01::OnOpacitySTLExternal(wxScrollEvent& event){
357 double value = ((double)stlSliderOpacityExternal->GetValue())/100;
358 actorExternal->GetProperty( )->SetOpacity( value );
363 void wxPathologyWidget_01::OnOpacitySTLInternal(wxScrollEvent& event){
364 double value = ((double)stlSliderOpacityInternal->GetValue())/100;
365 actorInternal->GetProperty( )->SetOpacity( value );
369 void wxPathologyWidget_01::OnBtnFileSTL(wxCommandEvent& event)
372 wxString dirSTL = _mar->_parameters->getStringParam(
373 marParameters::e_installation_directory );
374 dirSTL = ( dirSTL == _T("NO_DIRECTORY") ) ? wxGetHomeDir( ) : dirSTL;
375 wxDirDialog dialog( this, _T("Choose a directory..."), ( !dirSTL.IsEmpty( ) )?
376 dirSTL: wxGetHomeDir( ) );
378 if( dialog.ShowModal( ) == wxID_OK )
382 // ------------------------------------------------------------------------
383 // 1. GENERATE STL FILES
384 // ------------------------------------------------------------------------
385 const char* fileprefix = "c:\\Creatis\\";
386 std::string prefix = fileprefix;
389 // 1.1. Se hace un filtro triangular puesto que el stl writer solo recibe poligonos triangulares.
391 vtkTriangleFilter *filtro = vtkTriangleFilter::New();
392 filtro->SetInput(stlInterna);
393 vtkPolyDataConnectivityFilter *pdcf = vtkPolyDataConnectivityFilter::New();
394 pdcf->SetInput( filtro->GetOutput() );
395 vtkClosePolyData *cpd = vtkClosePolyData::New();
396 cpd->SetInput( pdcf->GetOutput() );
398 // 1.2 se escribe a disco el archivo stl de la superficie interna
400 vtkSTLWriter *writerI = vtkSTLWriter::New();
401 writerI->SetInput( cpd->GetOutput() );
403 prefix += "internal.stl";
404 writerI->SetFileName(prefix.c_str());
405 writerI->SetFileTypeToASCII();
409 // 1.3 se escribe a disco el archivo stl de la superficie externa
410 filtro->SetInput(stlExterna);
412 vtkSTLWriter *writerE = vtkSTLWriter::New();
413 writerE->SetInput( cpd->GetOutput() );
415 prefix += "external.stl";
416 writerE->SetFileName( prefix.c_str() );
417 writerE->SetFileTypeToASCII();
426 //By default *always* update e_installation_directory:
427 _mar->_parameters->setStringParam( marParameters::e_installation_directory, dialog.GetPath( ) );
428 _mar->saveParameters( );
432 void wxPathologyWidget_01::OnChangeSTLGaussLevel(wxScrollEvent& event)
434 stlDeltaGaussLevel = ((double)stlSliderDeltaGauss->GetValue())/100;
435 generateSTLSurfaces();
440 void wxPathologyWidget_01::OnChangeSTLMarchingCubesLevel(wxScrollEvent& event)
442 stlMarchingCubesLevel = ((double)stlSliderMarchingCubes->GetValue());
443 generateSTLSurfaces();
449 // ------------------------------------------------------------------------
450 // END STL FUNCTIONS - DHC
451 // ------------------------------------------------------------------------
454 // ------------------------------------------------------------------------
455 // START PATHOLOGY FUNCTIONS - DHC
456 // ------------------------------------------------------------------------
458 double GetDistanciaEuclideana(double* a, double* b){
459 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])));
463 double GetValorPromedio(vtkImageData *data){
471 data->GetExtent(ext);
474 for(i=ext[0];i<=ext[1];i++){
475 for(j=ext[2];j<=ext[3];j++){
476 for(k=ext[4];k<=ext[5];k++){
477 valor =data->GetScalarComponentAsDouble(i,j,k,0);
483 return (double)promedio/(double)numVoxels;
487 void wxPathologyWidget_01::ConfigurePathologyExtraction()
491 void wxPathologyWidget_01::generatePathologySurface(){
496 void wxPathologyWidget_01::OnBtnSetP1(wxCommandEvent &event){
497 vtk3DSurfaceWidget *sw = _maracasSurfaceWidget->GetVtk3DSurfaceWidget();
500 sw->GetRenderer()->RemoveActor(p1Actor);
504 sw->GetSphereCenter(p1Center);
506 p1SphereSource = vtkSphereSource::New();
507 p1Mapper = vtkPolyDataMapper::New();
508 p1Actor = vtkActor::New();
510 p1Mapper->SetInput(p1SphereSource->GetOutput());
511 p1Actor->SetMapper(p1Mapper);
512 p1Actor->PickableOff( );
514 p1SphereSource->SetCenter(p1Center[0],p1Center[1],p1Center[2]);
515 p1SphereSource->SetRadius(0.5);
517 p1Actor->GetProperty()->SetColor( 0.0, 1.0, 0.0);
519 sw->GetRenderer()->AddActor(p1Actor);
524 void wxPathologyWidget_01::OnBtnSetP2(wxCommandEvent &event){
526 vtk3DSurfaceWidget *sw = _maracasSurfaceWidget->GetVtk3DSurfaceWidget();
530 sw->GetRenderer()->RemoveActor(p1Actor);
533 sw->GetSphereCenter(p2Center);
535 p2SphereSource = vtkSphereSource::New();
536 p2Mapper = vtkPolyDataMapper::New();
537 p2Actor = vtkActor::New();
539 p2Mapper->SetInput(p2SphereSource->GetOutput());
540 p2Actor->SetMapper(p2Mapper);
541 p2Actor->PickableOff( );
543 p2SphereSource->SetCenter(p2Center[0],p2Center[1],p2Center[2]);
544 p2SphereSource->SetRadius(0.5);
546 p2Actor->GetProperty()->SetColor( 0.0, 1.0, 0.0);
548 sw->GetRenderer()->AddActor(p2Actor);
555 void wxPathologyWidget_01::OnBtnSetPat(wxCommandEvent &event){
557 if (p2SphereSource == NULL || p1SphereSource == NULL){
561 vtk3DSurfaceWidget *sw = _maracasSurfaceWidget->GetVtk3DSurfaceWidget();
564 sw->GetRenderer()->RemoveActor(patActor);
571 centro[0] = (p1Center[0] + p2Center[0]) / 2;
572 centro[1] = (p1Center[1] + p2Center[1]) / 2;
573 centro[2] = (p1Center[2] + p2Center[2]) / 2;
575 patSphereSource = vtkSphereSource::New();
576 patMapper = vtkPolyDataMapper::New();
577 patActor = vtkActor::New();
579 patMapper->SetInput(patSphereSource->GetOutput());
580 patActor->SetMapper(patMapper);
583 patSphereSource->SetCenter(centro[0],centro[1],centro[2]);
584 patSphereSource->SetRadius(GetDistanciaEuclideana(p1Center,p2Center));
586 patActor->GetProperty()->SetColor( 0.15, 0.75, 0.15 );
587 patActor->PickableOff( );
588 patActor->GetProperty( )->SetOpacity(0.15);
590 sw->GetRenderer()->AddActor(patActor);
597 void wxPathologyWidget_01::OnBtnExtractPat(wxCommandEvent &event){
600 //----------------------------------------------------------------------------------
601 vtkRenderer* renderer = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();
604 //----------------------------------------------------------------------------------
605 vtkImageData* originalData = vtkImageData::New();
606 originalData->DeepCopy(_mar->_experiment->getDynData()->getVolume()->castVtk());
609 originalData->GetSpacing(spacing);
612 originalData->GetExtent(extend);
616 //==================================================================================
617 // REGION DE INTERES EN LA IMAGEN ORIGINAL: patImageData
618 //==================================================================================
620 double p1crem[3], p2crem[3];
622 p1crem[0] = p1Center[0] / spacing[0];
623 p1crem[1] = p1Center[1] / spacing[1];
624 p1crem[2] = p1Center[2] / spacing[2];
626 p2crem[0] = p2Center[0] / spacing[0];
627 p2crem[1] = p2Center[1] / spacing[1];
628 p2crem[2] = p2Center[2] / spacing[2];
634 //COORDENADAS DE LA IMAGEN
636 double radio = GetDistanciaEuclideana(p1crem,p2crem) / 2;
639 centro[0] = (p1crem[0] + p2crem[0]) / 2;
640 centro[1] = (p1crem[1] + p2crem[1]) / 2;
641 centro[2] = (p1crem[2] + p2crem[2]) / 2;
645 //Obtener los ejes sobre la esfera para delimitar el area.
647 double lowX = centro[0] - (2*radio);
648 double highX = centro[0] + (2*radio);
649 double lowY = centro[1] - (2*radio);
650 double highY = centro[1] + (2*radio);
651 double lowZ = centro[2] - (2*radio);
652 double highZ = centro[2] + (2*radio);
656 //Extraer la region de interes:patImageData
657 //----------------------------------------------------------------------------------
659 voiFilter = vtkExtractVOI::New();
660 voiFilter->SetInput(originalData);
661 voiFilter->SetVOI((int)lowX, (int)highX, (int)lowY, (int)highY, (int)lowZ, (int)highZ);
663 patImageData = vtkImageData::New();
664 patImageData->DeepCopy(voiFilter->GetOutput());
665 // patImageData->SetSpacing(spacing);
667 //==================================================================================
668 // REGION DE INTERES ESPECIFICA
669 //==================================================================================
671 vtkOtsuSphereSource *otsu = new vtkOtsuSphereSource();
672 thresholdOTSU = otsu->calculateOptimalThreshold(patImageData);
675 //==================================================================================
677 //==================================================================================
680 cubesFilter = vtkMarchingCubes::New();
681 cubesFilter->SetInput(otsu->getImageForSegmentation());
682 cubesFilter->SetValue(0,patMCLevel);
683 cubesFilter->SetNumberOfContours(1);
685 //==================================================================================
687 //==================================================================================
689 polyUnico = vtkPolyDataConnectivityFilter::New();
690 polyUnico->SetInput(cubesFilter->GetOutput());
691 polyUnico->SetExtractionModeToLargestRegion();
693 stripper = vtkStripper::New();
694 stripper->SetInput(polyUnico->GetOutput());
697 //==================================================================================
699 //==================================================================================
701 dataMapper = vtkPolyDataMapper::New();
702 dataMapper->SetInput(stripper->GetOutput());
703 dataMapper->ScalarVisibilityOff();
705 //==================================================================================
707 //==================================================================================
710 dataActor = vtkActor::New();
711 dataActor->SetMapper(dataMapper);
712 dataActor->SetPosition(lowX*spacing[0],lowY*spacing[1],lowZ*spacing[2]);
713 dataActor->GetProperty()->SetOpacity(patOpacityLevel);
714 dataActor->GetProperty()->SetColor( 0.20, 0.20, 1.0 );
715 dataActor->GetProperty()->SetRepresentationToWireframe();
717 renderer->AddActor(dataActor);
721 outlineFilter = vtkOutlineFilter::New();
722 outlineFilter->SetInput(voiFilter->GetOutput());
724 vtkPolyDataMapper* outlineMapper = vtkPolyDataMapper::New();
725 outlineMapper->SetInput(outlineFilter->GetOutput());
728 renderer->RemoveActor(outlineActor);
732 renderer->RemoveActor(patActor);
736 renderer->RemoveActor(p3Actor);
740 outlineActor = vtkActor::New();
741 outlineActor->SetMapper(outlineMapper);
742 outlineActor->GetProperty()->SetOpacity(1);
745 /*p3SphereSource = vtkSphereSource::New();
746 p3Mapper = vtkPolyDataMapper::New();
747 p3Actor = vtkActor::New();
749 p3Mapper->SetInput(p3SphereSource->GetOutput());
750 p3Actor->SetMapper(p3Mapper);
751 p3Actor->PickableOff( );
753 p3SphereSource->SetCenter(p2Center[0]-(lowX*spacing[0]), p2Center[1]-(lowY*spacing[1]) , p2Center[2]-(lowZ*spacing[2]) );
754 p3SphereSource->SetRadius(0.5);
755 p3Actor->GetProperty()->SetColor( 0.9, 0.4, 0.2);*/
758 //renderer->AddActor(outlineActor);
759 //renderer->AddActor(p3Actor);
766 void wxPathologyWidget_01::OnChangePatOpacity(wxScrollEvent& event){
767 patOpacityLevel = ((double)patSliderOpacity->GetValue())/100.0;
768 dataActor->GetProperty( )->SetOpacity(patOpacityLevel);
772 void wxPathologyWidget_01::OnChangePatMarchingCubes(wxScrollEvent& event){
773 patMCLevel = ((double)patSliderMarchingCubes->GetValue());
774 cubesFilter->SetValue(0,patMCLevel);
775 cubesFilter->Update();
781 //Metodo auxiliar que devuelve la distancia minima de un punto a un polydata.
782 //Se evalua la distancia del punto a todos los centros de las celdas del polydata y se devuelve la menor
783 double GetDistanciaAPolydata(double* punto, vtkPoints* centers){
784 double distancia = 0.0;
785 double minima = 1e100;
787 int N = centers->GetNumberOfPoints();
789 for(int i = 0; i < N; i++){
790 centers->GetPoint(i,vertice);
791 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])));
793 if (distancia < minima){
801 double GetProductoPunto(double* vectorA, double* vectorB){
802 return (vectorA[0]*vectorB[0] + vectorA[1]*vectorB[1] + vectorA[2]*vectorB[2]);
806 double GetNormaVector(double* vector){
807 return sqrt(GetProductoPunto(vector, vector));
810 //construye el vector como la diferencia de dos puntos (vectores en el origen)
811 void GetVector(double* p0, double* p1, double* resultado){ //p1:destino - p0:origen
813 resultado[0] = p1[0] - p0[0];
814 resultado[1] = p1[1] - p0[1];
815 resultado[2] = p1[2] - p0[2];
820 void VectorUnitario(double *v){
821 double norma = GetNormaVector(v);
827 //calculo de las normales
828 void GetNormalCelda(double* p0, double* p1, double* p2, double* normal){
831 GetVector(p1,p2,U); //p2-p1
832 GetVector(p1,p0,V); //p0-p1
835 normal[0] = (U[1]*V[2])-(U[2]*V[1]);
836 normal[1] = -((U[0]*V[2])-(U[2]*V[0]));
837 normal[2] = (U[0]*V[1])-(U[1]*V[0]);
838 VectorUnitario(normal);
841 void GetNormal(vtkPoints* puntosCelda, double* normal){
842 double p0[3], p1[3], p2[3];
844 puntosCelda->GetPoint(0, p0);
845 puntosCelda->GetPoint(1, p1);
846 puntosCelda->GetPoint(2, p2);
848 GetNormalCelda(p0,p1,p2, normal);
852 // Obtiene el ID de la celda mas cercana
853 int GetClosestCellOnPolyData(double* punto, vtkPoints* centers){
854 double distancia = 0.0;
855 double minima = 1e100;
858 int N = centers->GetNumberOfPoints();
860 for(int i = 0; i < N; i++){
861 centers->GetPoint(i,vertice);
862 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])));
864 if (distancia < minima){
870 //printf("distancia del punto p: (%.2f,%.2f,%.2f) al polydata es %.2f (ID %i)\n", punto[0], punto[1], punto [2], minima, ID);
874 double GetProjectionOnNormal(double* vector, double* normal){
875 double pp = GetProductoPunto(vector, normal);
876 double vv = GetNormaVector(normal);
885 void wxPathologyWidget_01::OnBtnExtractAxis(wxCommandEvent& event){
889 //================================================================================
891 //================================================================================
892 vtkRenderer* renderer = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();
894 //===================================================================================
896 //===================================================================================
897 int extentOriginal[6];
898 this->patImageData->GetExtent(extentOriginal);
901 //===================================================================================
903 //===================================================================================
906 double resampling = 1.0/(double)factor;
910 // REMUESTREO DE LA IMAGEN ORIGINAL
911 //-----------------------------------------------------------------------------------
912 vtkImageResample* resamplingFilter = vtkImageResample::New();
913 resamplingFilter->SetAxisOutputSpacing(0,resampling);
914 resamplingFilter->SetAxisOutputSpacing(1,resampling);
915 resamplingFilter->SetAxisOutputSpacing(2,resampling);
916 resamplingFilter->ReleaseDataFlagOff();
917 resamplingFilter->SetInput(patImageData);
918 resamplingFilter->SetInterpolationModeToCubic();
919 resamplingFilter->Update();
922 //IMAGEN ORIGINAL REMUESTREADA
923 //-----------------------------------------------------------------------------------
924 vtkImageData* imagenRemuestreada = resamplingFilter->GetOutput();
925 imagenRemuestreada->ReleaseDataFlagOff();
929 //===================================================================================
930 // IMAGEN BINARIA = CRECIMIENTO DE REGIONES
931 //===================================================================================
933 this->isoBinaria = vtkImageData::New();
934 this->isoBinaria->CopyStructure(resamplingFilter->GetOutput());
935 this->isoBinaria->SetScalarTypeToDouble();
939 //CALCULO DE LAS DIMENSIONES Y EL EXTENT
940 //-----------------------------------------------------------------------------------
941 int extentBinaria[6];
942 this->isoBinaria->GetExtent(extentBinaria);
943 int *dim = isoBinaria->GetDimensions();
948 //-----------------------------------------------------------------------------------
949 int eX, eY, eZ, idpoint; // extent
952 FILE* logger = freopen("tests.txt","w",stdout);
956 for(eX = extentBinaria[0]; eX <= extentBinaria[1]; eX++){
957 for(eY = extentBinaria[2]; eY <= extentBinaria[3]; eY++){
958 for (eZ = extentBinaria[4]; eZ <= extentBinaria[5]; eZ++){
959 idpoint = isoBinaria->FindPoint(eX, eY, eZ);
960 ptr = (double *)isoBinaria->GetScalarPointer(eX,eY,eZ); // El apuntador se obtiene sobre el extent.
967 //===================================================================================
969 //===================================================================================
970 vtkPolyData *polyData= stripper->GetOutput(); //OJO SE TOMA SOLO LA SUPERFICIE MAS GRANDE (connectivity + stripper)
971 //===================================================================================
973 vtkPoints *puntosPolyData = polyData->GetPoints();
974 vtkPointData *pointData = polyData->GetPointData();
976 vtkPolyDataWriter* writePoly = vtkPolyDataWriter::New();
977 writePoly->SetFileName("poly.vtk");
978 writePoly->SetInput(polyData);
983 // PROCESAR LOS VERTICES
984 //----------------------------------------------------------------------------------
986 int numVertices = puntosPolyData->GetNumberOfPoints();
988 // mover el sistema de referencias del polydata para poder determinar lo que se encuentra
989 //por dentro segun las semillas.
992 /*for(j=0;j<numVertices;j++){
993 puntosPolyData->GetPoint(j,vertice);
994 vertice[0] = vertice[0] + extentOriginal[0];
995 vertice[1] = vertice[1] + extentOriginal[2];
996 vertice[2] = vertice[2] + extentOriginal[4];
997 puntosPolyData->SetPoint(j,vertice);
1002 // PROCESAR LAS CELDAS
1003 //----------------------------------------------------------------------------------
1006 //1. Obtener normales de la superficie (celdas)
1007 vtkPolyDataNormals* polyNormalsFilter = vtkPolyDataNormals::New(); //estas son las normales de los vertices NO de las celdas
1008 polyNormalsFilter->SetInput(polyData);
1009 polyNormalsFilter->SplittingOn();
1010 polyNormalsFilter->ConsistencyOn();
1011 polyNormalsFilter->NonManifoldTraversalOn();
1012 polyNormalsFilter->ComputeCellNormalsOn();
1013 polyNormalsFilter->Update();
1016 vtkPolyData* polyDataConNormales = polyNormalsFilter->GetOutput();
1019 //2. Obtener los centros de las celdas
1020 // Necesitamos el centro de las celdas para trazar el vector interno de un voxel
1021 //a una cara de la superficie.
1022 vtkCellCenters* cellCentersFilter = vtkCellCenters::New();
1023 cellCentersFilter->SetInput(polyDataConNormales);
1024 cellCentersFilter->VertexCellsOn();
1025 cellCentersFilter->Update();
1027 vtkPolyData* polydataCenters = cellCentersFilter->GetOutput();
1028 vtkPoints* polyCenterPoints = polydataCenters->GetPoints();
1030 vtkPoints* puntosCelda = NULL;
1031 vtkIdList* lista = NULL;
1032 vtkCell* celda = NULL;
1034 int numCeldas = polyCenterPoints->GetNumberOfPoints();
1037 vtkPoints *normales = vtkPoints::New();//almacena los componentes cartesianos de la normal de cada celda
1040 for(i = 0; i < numCeldas ; i++){
1042 polyCenterPoints->GetPoint(i,vertice);
1044 celda = polyDataConNormales->GetCell(i);
1045 puntosCelda = celda->GetPoints();
1046 GetNormal(puntosCelda, normal);
1047 normales->InsertPoint(i, normal);
1050 lista = celda->GetPointIds();
1051 for(j = 0; j < lista->GetNumberOfIds(); j++){
1052 int idPunto = lista->GetId(j);
1053 puntosCelda->GetPoint(j, vertice);
1061 //==========================================================================================
1062 //IMPLEMENTAR EL CRECIMIENTO DE REGIONES
1063 //==========================================================================================
1066 patImageData->GetSpacing(spacing);
1068 double p1crem[3], p2crem[3];
1070 p1crem[0] = p1Center[0] / spacing[0];
1071 p1crem[1] = p1Center[1] / spacing[1];
1072 p1crem[2] = p1Center[2] / spacing[2];
1074 p2crem[0] = p2Center[0] / spacing[0];
1075 p2crem[1] = p2Center[1] / spacing[1];
1076 p2crem[2] = p2Center[2] / spacing[2];
1082 //COORDENADAS DE LA IMAGEN
1084 double radio = GetDistanciaEuclideana(p1crem,p2crem) / 2;
1087 centro[0] = (p1crem[0] + p2crem[0]) / 2;
1088 centro[1] = (p1crem[1] + p2crem[1]) / 2;
1089 centro[2] = (p1crem[2] + p2crem[2]) / 2;
1093 //Obtener los ejes sobre la esfera para delimitar el area.
1095 double lowX = centro[0] - (2*radio);
1096 double highX = centro[0] + (2*radio);
1097 double lowY = centro[1] - (2*radio);
1098 double highY = centro[1] + (2*radio);
1099 double lowZ = centro[2] - (2*radio);
1100 double highZ = centro[2] + (2*radio);
1102 //================================================================================
1104 //================================================================================
1106 int numeroPuntos = dim[0] * dim[1] * dim[2];
1112 sem[0] = (int)(p2Center[0]-(lowX*spacing[0]));
1113 sem[1] = (int)(p2Center[1]-(lowY*spacing[1]));
1114 sem[2] = (int)(p2Center[2]-(lowZ*spacing[2]));
1117 int idSemilla = isoBinaria->FindPoint(sem[0],sem[1],sem[2]);
1119 if (idSemilla == -1) return;
1121 int graphSize = (extentBinaria[1]-extentBinaria[0]) * (extentBinaria[3]-extentBinaria[2])*(extentBinaria[5]-extentBinaria[4]);
1124 vtkPriorityQueue *colaEvaluacion = vtkPriorityQueue::New();
1126 colaEvaluacion->Allocate(graphSize);
1128 colaEvaluacion->Insert(0, idSemilla);
1130 vtkFloatingPointType prioridad;
1132 int longitudCola = colaEvaluacion->GetNumberOfItems();
1139 double* coordsVecino;
1141 double* coordsPunto;
1145 vtkPoints *resultado = vtkPoints::New();
1150 double centroCelda[3], vectorInterno[3], vectorNormal[3];
1152 while(longitudCola > 0){
1154 //================================================================================
1155 //1. Procesar el punto
1156 //================================================================================
1158 puntoID = colaEvaluacion->Pop(0,prioridad);
1160 coordsPunto = isoBinaria->GetPoint(puntoID);
1163 ptr = (double *)isoBinaria->GetScalarPointer((int)(coordsPunto[0]*factor), (int)(coordsPunto[1]*factor), (int)(coordsPunto[2]*factor));
1167 vtkImageWriter* writeImage = vtkImageWriter::New();
1168 writeImage->SetInput(isoBinaria);
1169 writeImage->SetFileName("sem.vtk");
1170 writeImage->Write();
1174 resultado->InsertPoint(cont, coordsPunto);
1177 //printf("iniciando en el punto (%.1f,%.1f,%.1f) \n", coordsPunto[0], coordsPunto[1], coordsPunto [2]);
1179 //================================================================================
1180 //2. Evaluar los vecinos del punto
1181 //================================================================================
1183 for(vk = -1; vk<2; vk++) {
1184 for(vj = -1; vj<2; vj++) {
1185 for(vi = -1; vi<2; vi++) {
1187 vecinoID = puntoID + (vk * dim[1]*dim[0]) + (vj * dim[0]) + vi;
1189 if( vecinoID >= 0 && vecinoID < numeroPuntos && vecinoID != 0 && vecinoID != puntoID ) {
1191 coordsVecino = isoBinaria->GetPoint(vecinoID);
1193 //ptr = (unsigned char *)isoBinaria->GetScalarPointer(coordsVecino[0],coordsVecino[1],coordsVecino[2]);
1194 ptr = (double *)isoBinaria->GetScalarPointer((int)(coordsVecino[0]*factor),(int)(coordsVecino[1]*factor),(int)(coordsVecino[2]*factor));
1196 //================================================================================
1197 // Para cada vecino no evaluado evaluar la distancia al polydata
1198 //================================================================================
1200 //printf("evaluando punto (%.1f,%.1f,%.1f) \n", coordsVecino[0], coordsVecino[1], coordsVecino [2]);
1202 int normalID = GetClosestCellOnPolyData(coordsVecino, polyCenterPoints);
1205 polyCenterPoints->GetPoint(normalID,centroCelda);
1206 normales->GetPoint(normalID, vectorNormal);
1208 GetVector(coordsVecino, centroCelda, vectorInterno);
1209 VectorUnitario(vectorInterno);
1211 double projection = GetProjectionOnNormal(vectorInterno, vectorNormal);
1212 ptr = (double *)isoBinaria->GetScalarPointer((int)(coordsVecino[0]*factor),(int)(coordsVecino[1]*factor),(int)(coordsVecino[2]*factor));
1214 if (projection > 0){
1216 //================================================================================
1217 //2.1 Si la proyeccion es positiva
1218 //colocar el vecino en la cola y poner el punto correspondiente en la imagen a 255.
1219 // (no ha pasado la frontera, el marcado a 255 se hace arriba en el while, ver paso 1)
1220 //================================================================================
1221 colaEvaluacion->Insert(prioridad,vecinoID);
1225 //================================================================================
1226 //2.2 Si la proyeccion es negativa con la normal mas cercana del polydata,
1227 //descartar el vecino == no poner en la cola == alcanzamos la superficie
1228 //================================================================================
1239 longitudCola = colaEvaluacion->GetNumberOfItems();
1241 colaEvaluacion->Delete();
1244 //================================================================================
1245 // MAPA DE DISTANCIAS
1246 //================================================================================
1248 this->isoDistance = vtkImageEuclideanDistance::New();
1249 this->isoDistance->SetInput(this->isoBinaria);
1250 this->isoDistance->InitializeOn();
1251 this->isoDistance->Update();
1252 vtkImageData *imDistancias = isoDistance->GetOutput();
1255 //================================================================================
1256 // CAMINO MINIMO SOBRE EL MAPA DE DISTANCIAS
1257 //================================================================================
1258 this->dijkstraFilter = vtkDijkstraImageData::New();
1259 this->dijkstraFilter->SetInput(isoDistance->GetOutput());
1261 int idA = (int)isoBinaria->FindPoint(p1Center[0],p1Center[1], p1Center[2]);
1262 int idB = (int)isoBinaria->FindPoint(p2Center[0],p2Center[1], p2Center[2]);
1264 dijkstraFilter->SourceID = idA;
1265 dijkstraFilter->SinkID = idB;
1266 dijkstraFilter->Update();
1268 caminoMapper = vtkPolyDataMapper::New();
1269 caminoMapper->SetInput(dijkstraFilter->GetOutput());
1270 caminoMapper->ScalarVisibilityOff();
1272 renderer->RemoveActor(dataActor);
1274 //=======================================================================================================
1275 //RECORRIDO DE LOS PUNTOS DEL EJE DETECTADO
1276 //=======================================================================================================
1277 vtkIdList *camino = dijkstraFilter->GetShortestPathIdList();
1279 int size = camino->GetNumberOfIds();
1283 for(i=0;i < size;i++){
1284 int pointId = camino->GetId(i);
1285 isoBinaria->GetPoint(pointId, coords);
1286 double scalar = imagenRemuestreada->GetPointData()->GetScalars()->GetTuple1(pointId);
1288 //----------------------------------------------------------------------------------------
1289 //Visualizacion del nodo del eje
1290 //----------------------------------------------------------------------------------------
1291 vtkSphereSource *_nodo = vtkSphereSource::New();
1292 _nodo->SetCenter(coords[0],coords[1],coords[2]);
1293 _nodo->SetRadius(0.2);
1295 vtkPolyDataMapper *_mapper = vtkPolyDataMapper::New();
1296 _mapper->SetInput(_nodo->GetOutput());
1299 vtkActor *_actor = vtkActor::New();
1300 _actor->GetProperty()->SetColor( 0.8, 0.2, 0.3);
1301 _actor->SetMapper(_mapper);
1302 _actor->PickableOff( );
1304 renderer->AddActor(_actor);
1311 // ------------------------------------------------------------------------
1312 // END PATHOLOGY FUNCTIONS - DHC
1313 // ------------------------------------------------------------------------