lib/maracasVisuLib/src/interface/wxWindows/widgets/wxPathologyWidget_01.cxx

   1
   2
   3
   4 #include "vtkActor.h"
   5 #include "vtkAppendFilter.h"
   6 #include "vtkCamera.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"
  27 #include "vtkMath.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"
  48
  49 #include "wxPathologyWidget_01.h"
  50 #include "kernel/vtkOtsuSphereSource.h"
  51 #include <wx/splitter.h>
  52
  53
  54
  55
  56 //-------------------------------------------------------------------
  57 //-------------------------------------------------------------------
  58 //-------------------------------------------------------------------
  59 wxPathologyWidget_01::wxPathologyWidget_01(wxWindow *parent, marInterface* mar)
  60 : wxPanel( parent, -1)
  61 {
  62
  63
  64         _mar                                                            = mar;
  65         wxBoxSizer                      *sizer                  = new wxBoxSizer(wxVERTICAL  );
  66     wxSplitterWindow    *pnlSplitter    = new wxSplitterWindow( this , -1);
  67         wxPanel                         *viewPanel              = CreateViewPanel(pnlSplitter);
  68         wxPanel                         *controlPanel   = CreateControlPanel(pnlSplitter);
  69
  70         sizer           -> Add( pnlSplitter ,1,wxGROW  ,0);
  71         pnlSplitter     -> SetMinimumPaneSize( 150 );
  72     pnlSplitter -> SplitVertically( viewPanel, controlPanel, 600 );
  73         this            -> SetSizer(sizer);
  74
  75
  76
  77         //DHC STL SURFACES
  78         stlInterna = NULL;
  79         stlExterna = NULL;
  80         stlImageData = NULL;
  81         workingDir = NULL;
  82
  83         //PATHOLOGY
  84         p1SphereSource = NULL;
  85         p1Mapper = NULL;
  86         p1Actor = NULL;
  87
  88         p2SphereSource = NULL;
  89         p2Mapper = NULL;
  90         p2Actor = NULL;
  91
  92         patSphereSource = NULL;
  93         patMapper = NULL;
  94         patActor = NULL;
  95
  96
  97         pathologyFrame  = NULL;
  98
  99         patImageData    = NULL;
 100         voiFilter               = NULL;
 101         cubesFilter             = NULL;
 102         dataMapper              = NULL;
 103         dataActor               = NULL;
 104
 105         outlineFilter   = NULL;
 106         outlineMapper   = NULL;
 107         outlineActor    = NULL;
 108
 109
 110         p3SphereSource  = NULL;
 111         p3Mapper                = NULL;
 112         p3Actor                 = NULL;
 113
 114
 115 }
 116 //-------------------------------------------------------------------
 117 wxPathologyWidget_01::~wxPathologyWidget_01(){
 118
 119 }
 120 //-------------------------------------------------------------------
 121 wxPanel* wxPathologyWidget_01::CreateViewPanel(wxWindow *parent)
 122 {
 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);
 130         panel->Layout();
 131         return panel;
 132 }
 133 //-------------------------------------------------------------------
 134 wxPanel* wxPathologyWidget_01::CreateControlPanel(wxWindow *parent)
 135 {
 136         wxPanel *panel          = new wxPanel(parent,-1);
 137         wxFlexGridSizer *sizer = new wxFlexGridSizer(2);
 138
 139
 140         //-- PATHOLOGY WIDGETS
 141
 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"));
 151         sizer->Add(btnP1);
 152         sizer->Add(btnP2);
 153         sizer->Add(btnPat);sizer->Add(new wxStaticText(panel,-1,_T("  ")));
 154         sizer->Add(btnExtract);sizer->Add(new wxStaticText(panel,-1,_T("  ")));
 155
 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
 164
 165
 166
 167         //-- STL WIDGETS
 168
 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 );
 175
 176         sizer->Add(new wxStaticText(panel,-1,_T(" Delta Gauss")));
 177         sizer->Add(stlSliderDeltaGauss);
 178
 179         sizer->Add(new wxStaticText(panel,-1,_T(" Marching Cubes Level")));
 180         sizer->Add(stlSliderMarchingCubes);
 181
 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);
 184
 185         sizer->Add(new wxStaticText(panel, -1,_T(" Opacity STL Internal")));
 186         sizer->Add(stlSliderOpacityInternal);
 187
 188         sizer->Add(new wxStaticText(panel, -1,_T(" Opacity STL External")));
 189         sizer->Add(stlSliderOpacityExternal);
 190
 191         wxButton *btnFileSTL = new wxButton(panel,-1,_T("Generate STL files"));
 192         sizer->Add(btnFileSTL);
 193         sizer->Add(new wxStaticText(panel,-1,_T("  ")));
 194          //-- STL WIDGETS
 195
 196
 197
 198
 199
 200
 201
 202
 203
 204         panel->SetSizer(sizer);
 205         panel->SetAutoLayout(true);
 206         panel->SetSize(400,600);
 207         panel->Layout();
 208
 209
 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);
 218
 219
 220         // -- PATHOLOGY CONNECT WIDGETS
 221
 222
 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
 230
 231
 232
 233
 234
 235         return panel;
 236 }
 237 //------------------------------------------------------------------------
 238 void wxPathologyWidget_01::Refresh()
 239 {
 240         _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->Render();
 241 }
 242
 243 //------------------------------------------------------------------------
 244 void wxPathologyWidget_01::ConfigureVTK()
 245 {
 246     wxBusyCursor wait;
 247
 248         vtkImageData    *imagedata              = _mar->_experiment->getDynData( )->getVolume( )->castVtk();
 249         //-Maracas-
 250         _maracasSurfaceWidget->ShowMARACASData( _mar );
 251
 252         //CONFIGURACION ADICIONAL
 253     this->ConfigureSTL();
 254         this->ConfigurePathologyExtraction();
 255
 256 }
 257
 258
 259 void wxPathologyWidget_01::ConfigureVTK(vtkImageData *imagedata, int x, int y, int z, double param){
 260
 261         this->px = x;
 262         this->py = y;
 263         this->pz = z;
 264
 265         wxBusyCursor wait;
 266         _maracasSurfaceWidget->ShowMARACASData( _mar );
 267
 268
 269         //CONFIGURACION ADICIONAL
 270         this->ConfigureMPR();
 271         this->ConfigureSTL();
 272         this->ConfigurePathologyExtraction();
 273 }
 274
 275
 276
 277
 278
 279 // ==================================================================================================
 280 // CODIGO NUEVO ADICIONADO POR DIEGO CANTOR
 281 // ==================================================================================================
 282
 283
 284 // ------------------------------------------------------------------------
 285 // START MPR FUNCTIONS - DHC
 286 // ------------------------------------------------------------------------
 287
 288 void wxPathologyWidget_01::ConfigureMPR(){
 289
 290         mprSphereSource = vtkSphereSource::New();
 291     mprMapper = vtkPolyDataMapper::New();
 292     mprActor = vtkActor::New();
 293
 294         mprMapper->SetInput(mprSphereSource->GetOutput());
 295         mprActor->SetMapper(mprMapper);
 296         mprActor->PickableOff( );
 297
 298         mprSphereSource->SetCenter(this->px,this->py,this->pz);
 299         mprSphereSource->SetRadius(2.5);
 300
 301         mprActor->GetProperty()->SetColor( 1.0, 0.0, 0.0 );
 302     vtkRenderer *ren = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();
 303         ren->AddActor(mprActor);
 304
 305 }
 306
 307 // ------------------------------------------------------------------------
 308 // END MPR FUNCTIONS - DHC
 309 // ------------------------------------------------------------------------
 310
 311
 312
 313 // ------------------------------------------------------------------------
 314 // START STL FUNCTIONS - DHC
 315 // ------------------------------------------------------------------------
 316
 317 void wxPathologyWidget_01::ConfigureSTL()
 318 {
 319         stlExterna = vtkPolyData::New();
 320         stlInterna = vtkPolyData::New();
 321
 322         dsm1 = vtkPolyDataMapper ::New();
 323     dsm1->SetInput (stlInterna);
 324     dsm1->ScalarVisibilityOff();
 325
 326     actorInternal = vtkActor::New();
 327     actorInternal->SetMapper (dsm1);
 328     actorInternal->GetProperty()->SetColor (0,1,0);
 329
 330     dsm2 = vtkPolyDataMapper ::New();
 331     dsm2->SetInput (stlExterna);
 332     dsm2->ScalarVisibilityOff();
 333
 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);
 340
 341         stlExtractor = new vtkSTLExtractor();
 342 }
 343
 344 void wxPathologyWidget_01::generateSTLSurfaces()
 345 {
 346         stlExtractor->setVolume(stlImageData);
 347         stlExtractor->setSigmaLevel(stlDeltaGaussLevel);
 348         stlExtractor->setMarchingCubesLevel(stlMarchingCubesLevel);
 349         stlExtractor->calculate();
 350
 351         stlInterna->DeepCopy(stlExtractor->getInnerSurface());
 352         stlExterna->DeepCopy(stlExtractor->getOuterSurface());
 353 }
 354
 355
 356 void wxPathologyWidget_01::OnOpacitySTLExternal(wxScrollEvent& event){
 357         double value = ((double)stlSliderOpacityExternal->GetValue())/100;
 358     actorExternal->GetProperty( )->SetOpacity( value );
 359         Refresh();
 360 }
 361
 362
 363 void wxPathologyWidget_01::OnOpacitySTLInternal(wxScrollEvent& event){
 364         double value = ((double)stlSliderOpacityInternal->GetValue())/100;
 365     actorInternal->GetProperty( )->SetOpacity( value );
 366         Refresh();
 367 }
 368
 369 void wxPathologyWidget_01::OnBtnFileSTL(wxCommandEvent& event)
 370 {
 371         wxBusyCursor wait;
 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( ) );
 377
 378         if( dialog.ShowModal( ) == wxID_OK )
 379         {
 380
 381
 382                 // ------------------------------------------------------------------------
 383                 //  1.  GENERATE STL FILES
 384                 // ------------------------------------------------------------------------
 385                 const char* fileprefix = "c:\\Creatis\\";
 386                 std::string prefix = fileprefix;
 387
 388
 389                 // 1.1. Se hace un filtro triangular puesto que el stl writer solo recibe poligonos triangulares.
 390
 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() );
 397
 398                 // 1.2 se escribe a disco el archivo stl de la superficie interna
 399         cpd->Update();
 400         vtkSTLWriter *writerI = vtkSTLWriter::New();
 401         writerI->SetInput( cpd->GetOutput() );
 402         prefix = fileprefix;
 403         prefix += "internal.stl";
 404         writerI->SetFileName(prefix.c_str());
 405         writerI->SetFileTypeToASCII();
 406         writerI->Write();
 407         writerI->Delete();
 408
 409                 // 1.3 se escribe a disco el archivo stl de la superficie externa
 410                 filtro->SetInput(stlExterna);
 411         cpd->Update();
 412         vtkSTLWriter *writerE = vtkSTLWriter::New();
 413         writerE->SetInput( cpd->GetOutput() );
 414         prefix = fileprefix;
 415         prefix += "external.stl";
 416         writerE->SetFileName( prefix.c_str() );
 417         writerE->SetFileTypeToASCII();
 418         writerE->Write();
 419         writerE->Delete();
 420
 421         filtro->Delete();
 422         cpd->Delete();
 423         pdcf->Delete();
 424         }
 425
 426         //By default *always* update e_installation_directory:
 427         _mar->_parameters->setStringParam( marParameters::e_installation_directory, dialog.GetPath( ) );
 428         _mar->saveParameters( );
 429 }
 430
 431
 432 void wxPathologyWidget_01::OnChangeSTLGaussLevel(wxScrollEvent& event)
 433 {
 434         stlDeltaGaussLevel  = ((double)stlSliderDeltaGauss->GetValue())/100;
 435         generateSTLSurfaces();
 436         Refresh();
 437 }
 438
 439
 440 void wxPathologyWidget_01::OnChangeSTLMarchingCubesLevel(wxScrollEvent& event)
 441 {
 442         stlMarchingCubesLevel = ((double)stlSliderMarchingCubes->GetValue());
 443         generateSTLSurfaces();
 444         Refresh();
 445
 446 }
 447
 448
 449 // ------------------------------------------------------------------------
 450 // END STL FUNCTIONS - DHC
 451 // ------------------------------------------------------------------------
 452
 453
 454 // ------------------------------------------------------------------------
 455 // START PATHOLOGY FUNCTIONS - DHC
 456 // ------------------------------------------------------------------------
 457
 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])));
 460 }
 461
 462
 463 double GetValorPromedio(vtkImageData *data){
 464         int ext[6];
 465         int i, j, k;
 466         int numVoxels = 0;
 467         double promedio = 0;
 468         double valor;
 469
 470
 471         data->GetExtent(ext);
 472
 473
 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);
 478                                 promedio += valor;
 479                                 numVoxels++;
 480                         }
 481                 }
 482         }
 483         return (double)promedio/(double)numVoxels;
 484 }
 485
 486
 487 void wxPathologyWidget_01::ConfigurePathologyExtraction()
 488 {
 489 }
 490
 491 void wxPathologyWidget_01::generatePathologySurface(){
 492         // 1. Inicializacion
 493 }
 494
 495
 496 void wxPathologyWidget_01::OnBtnSetP1(wxCommandEvent &event){
 497         vtk3DSurfaceWidget *sw = _maracasSurfaceWidget->GetVtk3DSurfaceWidget();
 498
 499         if(p1Actor){
 500                 sw->GetRenderer()->RemoveActor(p1Actor);
 501         }
 502
 503
 504         sw->GetSphereCenter(p1Center);
 505
 506         p1SphereSource = vtkSphereSource::New();
 507     p1Mapper = vtkPolyDataMapper::New();
 508     p1Actor = vtkActor::New();
 509
 510         p1Mapper->SetInput(p1SphereSource->GetOutput());
 511         p1Actor->SetMapper(p1Mapper);
 512         p1Actor->PickableOff( );
 513
 514         p1SphereSource->SetCenter(p1Center[0],p1Center[1],p1Center[2]);
 515         p1SphereSource->SetRadius(0.5);
 516
 517         p1Actor->GetProperty()->SetColor( 0.0, 1.0, 0.0);
 518
 519         sw->GetRenderer()->AddActor(p1Actor);
 520         Refresh();
 521
 522 }
 523
 524 void wxPathologyWidget_01::OnBtnSetP2(wxCommandEvent &event){
 525
 526         vtk3DSurfaceWidget *sw = _maracasSurfaceWidget->GetVtk3DSurfaceWidget();
 527
 528
 529         if(p2Actor){
 530                 sw->GetRenderer()->RemoveActor(p1Actor);
 531         }
 532
 533         sw->GetSphereCenter(p2Center);
 534
 535         p2SphereSource = vtkSphereSource::New();
 536     p2Mapper = vtkPolyDataMapper::New();
 537     p2Actor = vtkActor::New();
 538
 539         p2Mapper->SetInput(p2SphereSource->GetOutput());
 540         p2Actor->SetMapper(p2Mapper);
 541         p2Actor->PickableOff( );
 542
 543         p2SphereSource->SetCenter(p2Center[0],p2Center[1],p2Center[2]);
 544         p2SphereSource->SetRadius(0.5);
 545
 546         p2Actor->GetProperty()->SetColor( 0.0, 1.0, 0.0);
 547
 548         sw->GetRenderer()->AddActor(p2Actor);
 549         Refresh();
 550 }
 551
 552
 553
 554
 555 void wxPathologyWidget_01::OnBtnSetPat(wxCommandEvent &event){
 556
 557         if (p2SphereSource == NULL || p1SphereSource == NULL){
 558                 return;
 559         }
 560
 561         vtk3DSurfaceWidget *sw = _maracasSurfaceWidget->GetVtk3DSurfaceWidget();
 562
 563         if(patActor){
 564                 sw->GetRenderer()->RemoveActor(patActor);
 565         }
 566
 567
 568         double centro[3];
 569
 570
 571         centro[0] = (p1Center[0] + p2Center[0]) / 2;
 572         centro[1] = (p1Center[1] + p2Center[1]) / 2;
 573         centro[2] = (p1Center[2] + p2Center[2]) / 2;
 574
 575         patSphereSource = vtkSphereSource::New();
 576     patMapper = vtkPolyDataMapper::New();
 577     patActor = vtkActor::New();
 578
 579         patMapper->SetInput(patSphereSource->GetOutput());
 580         patActor->SetMapper(patMapper);
 581
 582
 583         patSphereSource->SetCenter(centro[0],centro[1],centro[2]);
 584         patSphereSource->SetRadius(GetDistanciaEuclideana(p1Center,p2Center));
 585
 586         patActor->GetProperty()->SetColor( 0.15, 0.75, 0.15 );
 587         patActor->PickableOff( );
 588         patActor->GetProperty( )->SetOpacity(0.15);
 589
 590         sw->GetRenderer()->AddActor(patActor);
 591         Refresh();
 592 }
 593
 594
 595
 596
 597 void wxPathologyWidget_01::OnBtnExtractPat(wxCommandEvent &event){
 598
 599         //RENDERER
 600         //----------------------------------------------------------------------------------
 601         vtkRenderer* renderer = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();
 602
 603         //IMAGEN ORIGINAL
 604         //----------------------------------------------------------------------------------
 605         vtkImageData* originalData = vtkImageData::New();
 606         originalData->DeepCopy(_mar->_experiment->getDynData()->getVolume()->castVtk());
 607
 608         double spacing[3];
 609         originalData->GetSpacing(spacing);
 610
 611         int extend[6];
 612         originalData->GetExtent(extend);
 613
 614
 615
 616         //==================================================================================
 617         // REGION DE INTERES EN LA IMAGEN ORIGINAL: patImageData
 618         //==================================================================================
 619
 620         double p1crem[3], p2crem[3];
 621
 622         p1crem[0] = p1Center[0] / spacing[0];
 623         p1crem[1] = p1Center[1] / spacing[1];
 624         p1crem[2] = p1Center[2] / spacing[2];
 625
 626         p2crem[0] = p2Center[0] / spacing[0];
 627         p2crem[1] = p2Center[1] / spacing[1];
 628         p2crem[2] = p2Center[2] / spacing[2];
 629
 630
 631
 632
 633
 634     //COORDENADAS DE LA IMAGEN
 635
 636         double radio = GetDistanciaEuclideana(p1crem,p2crem) / 2;
 637         double centro[3];
 638
 639         centro[0] = (p1crem[0] + p2crem[0]) / 2;
 640         centro[1] = (p1crem[1] + p2crem[1]) / 2;
 641         centro[2] = (p1crem[2] + p2crem[2]) / 2;
 642
 643
 644
 645         //Obtener los ejes sobre la esfera para delimitar el area.
 646
 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);
 653
 654
 655
 656         //Extraer la region de interes:patImageData
 657         //----------------------------------------------------------------------------------
 658
 659         voiFilter = vtkExtractVOI::New();
 660         voiFilter->SetInput(originalData);
 661         voiFilter->SetVOI((int)lowX, (int)highX, (int)lowY, (int)highY, (int)lowZ, (int)highZ);
 662         voiFilter->Update();
 663         patImageData = vtkImageData::New();
 664         patImageData->DeepCopy(voiFilter->GetOutput());
 665 //      patImageData->SetSpacing(spacing);
 666
 667         //==================================================================================
 668         // REGION DE INTERES ESPECIFICA
 669         //==================================================================================
 670
 671         vtkOtsuSphereSource *otsu = new vtkOtsuSphereSource();
 672         thresholdOTSU = otsu->calculateOptimalThreshold(patImageData);
 673
 674
 675         //==================================================================================
 676         // MARCHING CUBES
 677         //==================================================================================
 678
 679         patMCLevel = 128.0;
 680         cubesFilter = vtkMarchingCubes::New();
 681         cubesFilter->SetInput(otsu->getImageForSegmentation());
 682         cubesFilter->SetValue(0,patMCLevel);
 683         cubesFilter->SetNumberOfContours(1);
 684
 685         //==================================================================================
 686         // POLYDATA
 687         //==================================================================================
 688
 689     polyUnico = vtkPolyDataConnectivityFilter::New();
 690         polyUnico->SetInput(cubesFilter->GetOutput());
 691         polyUnico->SetExtractionModeToLargestRegion();
 692
 693         stripper = vtkStripper::New();
 694         stripper->SetInput(polyUnico->GetOutput());
 695
 696
 697         //==================================================================================
 698         // DATA MAPPING
 699         //==================================================================================
 700
 701         dataMapper = vtkPolyDataMapper::New();
 702         dataMapper->SetInput(stripper->GetOutput());
 703         dataMapper->ScalarVisibilityOff();
 704
 705         //==================================================================================
 706         // ACTORES
 707         //==================================================================================
 708
 709
 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();
 716
 717         renderer->AddActor(dataActor);
 718
 719
 720
 721         outlineFilter = vtkOutlineFilter::New();
 722     outlineFilter->SetInput(voiFilter->GetOutput());
 723
 724         vtkPolyDataMapper* outlineMapper = vtkPolyDataMapper::New();
 725         outlineMapper->SetInput(outlineFilter->GetOutput());
 726
 727         if(outlineActor){
 728                 renderer->RemoveActor(outlineActor);
 729         }
 730
 731         if(patActor){
 732                 renderer->RemoveActor(patActor);
 733         }
 734
 735         if(p3Actor){
 736                 renderer->RemoveActor(p3Actor);
 737         }
 738
 739
 740         outlineActor = vtkActor::New();
 741         outlineActor->SetMapper(outlineMapper);
 742         outlineActor->GetProperty()->SetOpacity(1);
 743
 744
 745         /*p3SphereSource = vtkSphereSource::New();
 746     p3Mapper = vtkPolyDataMapper::New();
 747     p3Actor = vtkActor::New();
 748
 749         p3Mapper->SetInput(p3SphereSource->GetOutput());
 750         p3Actor->SetMapper(p3Mapper);
 751         p3Actor->PickableOff( );
 752
 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);*/
 756
 757
 758         //renderer->AddActor(outlineActor);
 759         //renderer->AddActor(p3Actor);
 760
 761
 762         Refresh();
 763 }
 764
 765
 766 void wxPathologyWidget_01::OnChangePatOpacity(wxScrollEvent& event){
 767         patOpacityLevel = ((double)patSliderOpacity->GetValue())/100.0;
 768     dataActor->GetProperty( )->SetOpacity(patOpacityLevel);
 769         Refresh();
 770 }
 771
 772 void wxPathologyWidget_01::OnChangePatMarchingCubes(wxScrollEvent& event){
 773         patMCLevel = ((double)patSliderMarchingCubes->GetValue());
 774         cubesFilter->SetValue(0,patMCLevel);
 775         cubesFilter->Update();
 776         Refresh();
 777 }
 778
 779
 780
 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;
 786         double vertice[3];
 787         int N = centers->GetNumberOfPoints();
 788
 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])));
 792
 793                  if (distancia < minima){
 794                          minima = distancia;
 795                  }
 796         }
 797         return minima;
 798 }
 799
 800
 801 double GetProductoPunto(double* vectorA, double* vectorB){
 802         return (vectorA[0]*vectorB[0] +  vectorA[1]*vectorB[1] + vectorA[2]*vectorB[2]);
 803 }
 804
 805
 806 double GetNormaVector(double* vector){
 807         return sqrt(GetProductoPunto(vector, vector));
 808 }
 809
 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
 812
 813         resultado[0] = p1[0] - p0[0];
 814         resultado[1] = p1[1] - p0[1];
 815         resultado[2] = p1[2] - p0[2];
 816
 817 }
 818
 819
 820 void VectorUnitario(double *v){
 821         double norma = GetNormaVector(v);
 822         v[0] = v[0] / norma;
 823         v[1] = v[1] / norma;
 824         v[2] = v[2] / norma;
 825 }
 826
 827 //calculo de las normales
 828 void GetNormalCelda(double* p0, double* p1, double* p2, double* normal){
 829
 830         double U[3],V[3];
 831         GetVector(p1,p2,U); //p2-p1
 832         GetVector(p1,p0,V); //p0-p1
 833
 834
 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);
 839 }
 840
 841 void GetNormal(vtkPoints* puntosCelda, double* normal){
 842         double p0[3], p1[3], p2[3];
 843
 844         puntosCelda->GetPoint(0, p0);
 845         puntosCelda->GetPoint(1, p1);
 846         puntosCelda->GetPoint(2, p2);
 847
 848         GetNormalCelda(p0,p1,p2, normal);
 849 }
 850
 851
 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;
 856         int ID = -1;
 857         double vertice[3];
 858         int N = centers->GetNumberOfPoints();
 859
 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])));
 863
 864                  if (distancia < minima){
 865                          minima = distancia;
 866                          ID = i;
 867                  }
 868         }
 869
 870         //printf("distancia del punto p: (%.2f,%.2f,%.2f) al polydata es %.2f (ID %i)\n", punto[0], punto[1], punto [2], minima, ID);
 871         return ID;
 872 }
 873
 874 double GetProjectionOnNormal(double* vector, double* normal){
 875         double pp = GetProductoPunto(vector, normal);
 876         double vv = GetNormaVector(normal);
 877         double pr = pp / vv;
 878         return pr;
 879 }
 880
 881
 882
 883
 884
 885 void wxPathologyWidget_01::OnBtnExtractAxis(wxCommandEvent& event){
 886
 887
 888
 889         //================================================================================
 890         // RENDERER
 891         //================================================================================
 892         vtkRenderer* renderer = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();
 893
 894         //===================================================================================
 895         // IMAGEN ORIGINAL.
 896         //===================================================================================
 897         int extentOriginal[6];
 898         this->patImageData->GetExtent(extentOriginal);
 899
 900
 901         //===================================================================================
 902         // REMUESTREO
 903         //===================================================================================
 904
 905         int factor = 1;
 906         double resampling = 1.0/(double)factor;
 907
 908
 909
 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();
 920
 921
 922         //IMAGEN ORIGINAL REMUESTREADA
 923         //-----------------------------------------------------------------------------------
 924         vtkImageData* imagenRemuestreada = resamplingFilter->GetOutput();
 925         imagenRemuestreada->ReleaseDataFlagOff();
 926
 927
 928
 929         //===================================================================================
 930         // IMAGEN BINARIA = CRECIMIENTO DE REGIONES
 931         //===================================================================================
 932
 933         this->isoBinaria = vtkImageData::New();
 934         this->isoBinaria->CopyStructure(resamplingFilter->GetOutput());
 935         this->isoBinaria->SetScalarTypeToDouble();
 936
 937
 938
 939         //CALCULO DE LAS DIMENSIONES Y EL EXTENT
 940         //-----------------------------------------------------------------------------------
 941         int extentBinaria[6];
 942         this->isoBinaria->GetExtent(extentBinaria);
 943         int *dim = isoBinaria->GetDimensions();
 944
 945
 946
 947         // IMAGEN EN NEGRO.
 948         //-----------------------------------------------------------------------------------
 949         int eX, eY, eZ, idpoint; // extent
 950         double *ptr = NULL;
 951
 952         FILE* logger = freopen("tests.txt","w",stdout);
 953
 954
 955
 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.
 961                                 *ptr = 0;
 962                         }
 963                 }
 964         }
 965
 966
 967     //===================================================================================
 968         // POLYDATA
 969         //===================================================================================
 970      vtkPolyData *polyData= stripper->GetOutput();  //OJO SE TOMA SOLO LA SUPERFICIE MAS GRANDE (connectivity + stripper)
 971          //===================================================================================
 972
 973          vtkPoints *puntosPolyData = polyData->GetPoints();
 974          vtkPointData *pointData = polyData->GetPointData();
 975
 976          vtkPolyDataWriter* writePoly = vtkPolyDataWriter::New();
 977          writePoly->SetFileName("poly.vtk");
 978          writePoly->SetInput(polyData);
 979          writePoly->Write();
 980
 981
 982
 983         // PROCESAR LOS VERTICES
 984         //----------------------------------------------------------------------------------
 985
 986          int numVertices = puntosPolyData->GetNumberOfPoints();
 987
 988          // mover el sistema de referencias del polydata para poder determinar lo que se encuentra
 989          //por dentro segun las semillas.
 990          int j;
 991          double vertice[3];
 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);
 998          }*/
 999
1000
1001
1002         // PROCESAR LAS CELDAS
1003         //----------------------------------------------------------------------------------
1004
1005
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();
1014
1015
1016          vtkPolyData* polyDataConNormales = polyNormalsFilter->GetOutput();
1017
1018
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();
1026
1027          vtkPolyData* polydataCenters = cellCentersFilter->GetOutput();
1028          vtkPoints*   polyCenterPoints  = polydataCenters->GetPoints();
1029
1030          vtkPoints*   puntosCelda = NULL;
1031          vtkIdList*   lista = NULL;
1032          vtkCell*     celda = NULL;
1033      double normal[3];
1034          int numCeldas = polyCenterPoints->GetNumberOfPoints();
1035
1036
1037          vtkPoints *normales = vtkPoints::New();//almacena los componentes cartesianos de la normal de cada celda
1038
1039      int i;
1040      for(i = 0; i < numCeldas ; i++){
1041
1042                 polyCenterPoints->GetPoint(i,vertice);
1043
1044             celda = polyDataConNormales->GetCell(i);
1045                 puntosCelda = celda->GetPoints();
1046         GetNormal(puntosCelda, normal);
1047                 normales->InsertPoint(i, normal);
1048
1049
1050                 lista = celda->GetPointIds();
1051                 for(j = 0; j < lista->GetNumberOfIds(); j++){
1052                         int idPunto = lista->GetId(j);
1053                         puntosCelda->GetPoint(j, vertice);
1054
1055                 }
1056
1057          }
1058
1059
1060
1061          //==========================================================================================
1062          //IMPLEMENTAR EL CRECIMIENTO DE REGIONES
1063          //==========================================================================================
1064
1065         double spacing[3];
1066         patImageData->GetSpacing(spacing);
1067
1068         double p1crem[3], p2crem[3];
1069
1070         p1crem[0] = p1Center[0] / spacing[0];
1071         p1crem[1] = p1Center[1] / spacing[1];
1072         p1crem[2] = p1Center[2] / spacing[2];
1073
1074         p2crem[0] = p2Center[0] / spacing[0];
1075         p2crem[1] = p2Center[1] / spacing[1];
1076         p2crem[2] = p2Center[2] / spacing[2];
1077
1078
1079
1080
1081
1082     //COORDENADAS DE LA IMAGEN
1083
1084         double radio = GetDistanciaEuclideana(p1crem,p2crem) / 2;
1085         double centro[3];
1086
1087         centro[0] = (p1crem[0] + p2crem[0]) / 2;
1088         centro[1] = (p1crem[1] + p2crem[1]) / 2;
1089         centro[2] = (p1crem[2] + p2crem[2]) / 2;
1090
1091
1092
1093         //Obtener los ejes sobre la esfera para delimitar el area.
1094
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);
1101
1102         //================================================================================
1103         //0. Inicializacion
1104         //================================================================================
1105
1106         int numeroPuntos = dim[0] * dim[1] * dim[2];
1107
1108
1109         int sem[3];
1110
1111
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]));
1115
1116
1117         int idSemilla = isoBinaria->FindPoint(sem[0],sem[1],sem[2]);
1118
1119         if (idSemilla == -1) return;
1120
1121         int graphSize = (extentBinaria[1]-extentBinaria[0]) * (extentBinaria[3]-extentBinaria[2])*(extentBinaria[5]-extentBinaria[4]);
1122
1123
1124         vtkPriorityQueue *colaEvaluacion = vtkPriorityQueue::New();
1125
1126         colaEvaluacion->Allocate(graphSize);
1127
1128         colaEvaluacion->Insert(0, idSemilla);
1129
1130         vtkFloatingPointType prioridad;
1131
1132         int longitudCola = colaEvaluacion->GetNumberOfItems();
1133
1134
1135         int vecinoID = -1;
1136
1137         int puntoID = -1;
1138
1139         double* coordsVecino;
1140
1141         double* coordsPunto;
1142
1143         int cont=0;
1144
1145         vtkPoints *resultado = vtkPoints::New();
1146
1147         int premiere = 1;
1148
1149
1150         double centroCelda[3], vectorInterno[3], vectorNormal[3];
1151
1152         while(longitudCola > 0){
1153
1154                         //================================================================================
1155                         //1. Procesar el punto
1156                         //================================================================================
1157
1158                     puntoID = colaEvaluacion->Pop(0,prioridad);
1159
1160                         coordsPunto = isoBinaria->GetPoint(puntoID);
1161
1162
1163                     ptr = (double *)isoBinaria->GetScalarPointer((int)(coordsPunto[0]*factor), (int)(coordsPunto[1]*factor), (int)(coordsPunto[2]*factor));
1164                         *ptr = 255.0;
1165
1166                         if(premiere == 1){
1167                                 vtkImageWriter* writeImage = vtkImageWriter::New();
1168                                 writeImage->SetInput(isoBinaria);
1169                                 writeImage->SetFileName("sem.vtk");
1170                                 writeImage->Write();
1171                                 premiere = 0;
1172                         }
1173
1174                         resultado->InsertPoint(cont, coordsPunto);
1175                         cont++;
1176
1177                         //printf("iniciando en el punto (%.1f,%.1f,%.1f) \n", coordsPunto[0], coordsPunto[1], coordsPunto [2]);
1178
1179                         //================================================================================
1180                         //2. Evaluar los vecinos del punto
1181                         //================================================================================
1182                         int vk,vj,vi;
1183                         for(vk = -1; vk<2; vk++) {
1184                                 for(vj = -1; vj<2; vj++) {
1185                                         for(vi = -1; vi<2; vi++) {
1186
1187                                                 vecinoID = puntoID + (vk * dim[1]*dim[0]) + (vj * dim[0]) + vi;
1188
1189                                                 if( vecinoID >= 0 && vecinoID < numeroPuntos && vecinoID != 0 && vecinoID != puntoID ) {
1190
1191                                                                 coordsVecino = isoBinaria->GetPoint(vecinoID);
1192
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));
1195
1196                                                                 //================================================================================
1197                                                                 // Para cada vecino no evaluado evaluar la distancia al polydata
1198                                                                 //================================================================================
1199                                                                 if (*ptr == 0){
1200                                                                         //printf("evaluando punto (%.1f,%.1f,%.1f) \n", coordsVecino[0], coordsVecino[1], coordsVecino [2]);
1201
1202                                                                         int normalID = GetClosestCellOnPolyData(coordsVecino, polyCenterPoints);
1203
1204
1205                                                                         polyCenterPoints->GetPoint(normalID,centroCelda);
1206                                                                         normales->GetPoint(normalID, vectorNormal);
1207
1208                                                                         GetVector(coordsVecino, centroCelda, vectorInterno);
1209                                                                         VectorUnitario(vectorInterno);
1210
1211                                                                         double projection = GetProjectionOnNormal(vectorInterno, vectorNormal);
1212                                                                         ptr = (double *)isoBinaria->GetScalarPointer((int)(coordsVecino[0]*factor),(int)(coordsVecino[1]*factor),(int)(coordsVecino[2]*factor));
1213
1214                                                                         if (projection > 0){
1215
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);
1222
1223                                                                         } else {
1224
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                                                                         //================================================================================
1229
1230                                                                                 *ptr = 0;
1231                                                                         }
1232                                                                 }
1233
1234                                                 }
1235                                         }
1236                                 }
1237                         }
1238
1239                         longitudCola = colaEvaluacion->GetNumberOfItems();
1240         }
1241         colaEvaluacion->Delete();
1242
1243
1244         //================================================================================
1245         //  MAPA DE DISTANCIAS
1246         //================================================================================
1247
1248         this->isoDistance = vtkImageEuclideanDistance::New();
1249         this->isoDistance->SetInput(this->isoBinaria);
1250         this->isoDistance->InitializeOn();
1251         this->isoDistance->Update();
1252         vtkImageData *imDistancias = isoDistance->GetOutput();
1253
1254
1255         //================================================================================
1256         //  CAMINO MINIMO SOBRE EL MAPA DE DISTANCIAS
1257         //================================================================================
1258         this->dijkstraFilter = vtkDijkstraImageData::New();
1259         this->dijkstraFilter->SetInput(isoDistance->GetOutput());
1260
1261         int idA = (int)isoBinaria->FindPoint(p1Center[0],p1Center[1], p1Center[2]);
1262         int idB = (int)isoBinaria->FindPoint(p2Center[0],p2Center[1], p2Center[2]);
1263
1264         dijkstraFilter->SourceID = idA;
1265         dijkstraFilter->SinkID = idB;
1266         dijkstraFilter->Update();
1267
1268         caminoMapper = vtkPolyDataMapper::New();
1269         caminoMapper->SetInput(dijkstraFilter->GetOutput());
1270         caminoMapper->ScalarVisibilityOff();
1271
1272         renderer->RemoveActor(dataActor);
1273
1274         //=======================================================================================================
1275         //RECORRIDO DE LOS PUNTOS DEL EJE DETECTADO
1276         //=======================================================================================================
1277         vtkIdList *camino = dijkstraFilter->GetShortestPathIdList();
1278
1279         int size = camino->GetNumberOfIds();
1280         double coords[3];
1281
1282
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);
1287
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);
1294
1295                 vtkPolyDataMapper *_mapper = vtkPolyDataMapper::New();
1296                 _mapper->SetInput(_nodo->GetOutput());
1297
1298
1299         vtkActor          *_actor  = vtkActor::New();
1300                 _actor->GetProperty()->SetColor( 0.8, 0.2, 0.3);
1301             _actor->SetMapper(_mapper);
1302                 _actor->PickableOff( );
1303
1304                 renderer->AddActor(_actor);
1305
1306
1307         }
1308
1309
1310 }
1311 // ------------------------------------------------------------------------
1312 // END PATHOLOGY FUNCTIONS - DHC
1313 // ------------------------------------------------------------------------
1314
1315
1316