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

   1 /*# ---------------------------------------------------------------------
   2 #
   3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
   4 #                        pour la Sant�)
   5 # Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
   6 # Previous Authors : Laurent Guigues, Jean-Pierre Roux
   7 # CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil
   8 #
   9 #  This software is governed by the CeCILL-B license under French law and
  10 #  abiding by the rules of distribution of free software. You can  use,
  11 #  modify and/ or redistribute the software under the terms of the CeCILL-B
  12 #  license as circulated by CEA, CNRS and INRIA at the following URL
  13 #  http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  14 #  or in the file LICENSE.txt.
  15 #
  16 #  As a counterpart to the access to the source code and  rights to copy,
  17 #  modify and redistribute granted by the license, users are provided only
  18 #  with a limited warranty  and the software's author,  the holder of the
  19 #  economic rights,  and the successive licensors  have only  limited
  20 #  liability.
  21 #
  22 #  The fact that you are presently reading this means that you have had
  23 #  knowledge of the CeCILL-B license and that you accept its terms.
  24 # ------------------------------------------------------------------------ */
  25
  26 #include "vtkActor.h"
  27 #include "vtkAppendFilter.h"
  28 #include "vtkCamera.h"
  29 #include "vtkCallbackCommand.h"
  30 #include "vtkCellArray.h"
  31 #include "vtkCylinderSource.h"
  32 #include "vtkClosePolyData.h"
  33 #include "vtkCommand.h"
  34 #include "vtkDoubleArray.h"
  35 #include "vtkExtractVOI.h"
  36 #include "vtkImageReader.h"
  37 #include "vtkImageViewer.h"
  38 #include "vtkImageViewer2.h"
  39 #include "vtkImageToStructuredPoints.h"
  40 #include "vtkImageThreshold.h"
  41 #include "vtkImageWriter.h"
  42 #include "vtkImageClip.h"
  43 #include "vtkImageResample.h"
  44 #include "vtkImageReslice.h"
  45 #include "vtkImageResample.h"
  46 #include "vtkImageThreshold.h"
  47 #include "vtkImageCast.h"
  48 #include "vtkImageData.h"
  49 #include "vtkMath.h"
  50 #include "vtkMarchingCubes.h"
  51 #include "vtkObjectFactory.h"
  52 #include "vtkPolyDataMapper.h"
  53 #include "vtkRenderer.h"
  54 #include "vtkRenderWindow.h"
  55 #include "vtkRenderWindowInteractor.h"
  56 #include "vtkPolyDataConnectivityFilter.h"
  57 #include "vtkProperty.h"
  58 #include "vtkPriorityQueue.h"
  59 #include "vtkPoints.h"
  60 #include "vtkPolyData.h"
  61 #include "vtkPolyDataMapper.h"
  62 #include "vtkPolyDataWriter.h"
  63 #include "vtkPolyDataReader.h"
  64 #include "vtkPointData.h"
  65 #include "vtkSphereSource.h"
  66 #include "vtkSTLWriter.h"
  67 #include "vtkStripper.h"
  68 #include "vtkTriangleFilter.h"
  69 #include "vtkTransform.h"
  70
  71 #include "wxPathologyWidget_01.h"
  72 #include "kernel/vtkOtsuSphereSource.h"
  73 #include <wx/splitter.h>
  74
  75
  76
  77
  78 //-------------------------------------------------------------------
  79 //-------------------------------------------------------------------
  80 //-------------------------------------------------------------------
  81 wxPathologyWidget_01::wxPathologyWidget_01(wxWindow *parent, marInterface* mar)
  82 : wxPanel( parent, -1)
  83 {
  84
  85
  86         _mar                                                            = mar;
  87         wxBoxSizer                      *sizer                  = new wxBoxSizer(wxVERTICAL  );
  88     wxSplitterWindow    *pnlSplitter    = new wxSplitterWindow( this , -1);
  89         wxPanel                         *viewPanel              = CreateViewPanel(pnlSplitter);
  90         wxPanel                         *controlPanel   = CreateControlPanel(pnlSplitter);
  91
  92         sizer           -> Add( pnlSplitter ,1,wxGROW  ,0);
  93         pnlSplitter     -> SetMinimumPaneSize( 150 );
  94     pnlSplitter -> SplitVertically( viewPanel, controlPanel, 600 );
  95         this            -> SetSizer(sizer);
  96
  97
  98
  99         //DHC STL SURFACES
 100         stlInterna = NULL;
 101         stlExterna = NULL;
 102         stlImageData = NULL;
 103         workingDir = NULL;
 104
 105         //PATHOLOGY
 106         p1SphereSource = NULL;
 107         p1Mapper = NULL;
 108         p1Actor = NULL;
 109
 110         p2SphereSource = NULL;
 111         p2Mapper = NULL;
 112         p2Actor = NULL;
 113
 114         patSphereSource = NULL;
 115         patMapper = NULL;
 116         patActor = NULL;
 117
 118
 119         pathologyFrame  = NULL;
 120
 121         patImageData    = NULL;
 122         voiFilter               = NULL;
 123         cubesFilter             = NULL;
 124         dataMapper              = NULL;
 125         dataActor               = NULL;
 126
 127         outlineFilter   = NULL;
 128         outlineMapper   = NULL;
 129         outlineActor    = NULL;
 130
 131
 132         p3SphereSource  = NULL;
 133         p3Mapper                = NULL;
 134         p3Actor                 = NULL;
 135
 136
 137 }
 138 //-------------------------------------------------------------------
 139 wxPathologyWidget_01::~wxPathologyWidget_01(){
 140
 141 }
 142 //-------------------------------------------------------------------
 143 wxPanel* wxPathologyWidget_01::CreateViewPanel(wxWindow *parent)
 144 {
 145         wxPanel *panel          =       new wxPanel(parent,-1);
 146     wxBoxSizer *sizer   =       new wxBoxSizer(wxVERTICAL);
 147         _maracasSurfaceWidget = new wxSurfaceWidget(panel);
 148         sizer->Add(_maracasSurfaceWidget        , 1, wxEXPAND, 0);
 149         panel->SetSizer(sizer);
 150         panel->SetAutoLayout(true);
 151         panel->SetSize(400,400);
 152         panel->Layout();
 153         return panel;
 154 }
 155 //-------------------------------------------------------------------
 156 wxPanel* wxPathologyWidget_01::CreateControlPanel(wxWindow *parent)
 157 {
 158         wxPanel *panel          = new wxPanel(parent,-1);
 159         wxFlexGridSizer *sizer = new wxFlexGridSizer(2);
 160
 161
 162         //-- PATHOLOGY WIDGETS
 163
 164         sizer->Add(new wxStaticText(panel,-1,_T("  ")));
 165         sizer->Add(new wxStaticText(panel,-1,_T("  ")));
 166     sizer->Add(new wxStaticText(panel,-1,_T("  - - -  PATHOLOGY EXTRACTION - - - ")));
 167         sizer->Add(new wxStaticText(panel,-1,_T("  ")));
 168         wxButton *btnP1 = new wxButton(panel,-1,_T("Set P1"));
 169         wxButton *btnP2 = new wxButton(panel,-1,_T("Set P2"));
 170         wxButton *btnPat = new wxButton(panel,-1,_T("Set Region"));
 171         wxButton *btnExtract = new wxButton(panel,-1,_T("Extract Region"));
 172         wxButton *btnAxis = new wxButton(panel,-1, _T("Axis"));
 173         sizer->Add(btnP1);
 174         sizer->Add(btnP2);
 175         sizer->Add(btnPat);sizer->Add(new wxStaticText(panel,-1,_T("  ")));
 176         sizer->Add(btnExtract);sizer->Add(new wxStaticText(panel,-1,_T("  ")));
 177
 178         patSliderOpacity        = new wxSlider( panel, -1, 50 , 0, 100, wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS );
 179         sizer->Add(new wxStaticText(panel, -1,_T("Opacity Pathological VOI")));
 180         sizer->Add(patSliderOpacity);
 181         patSliderMarchingCubes  = new wxSlider( panel, -1, 254 , 0, 1000, wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS );
 182         sizer->Add(new wxStaticText(panel, -1,_T(" Marching Cubes Level Pathological VOI")));
 183         sizer->Add(patSliderMarchingCubes);
 184         sizer->Add(btnAxis);sizer->Add(new wxStaticText(panel,-1,_T("  ")));
 185         //-- PATHOLOGY WIDGETS
 186
 187
 188
 189         //-- STL WIDGETS
 190
 191         sizer->Add(new wxStaticText(panel,-1,_T("  ")));
 192         sizer->Add(new wxStaticText(panel,-1,_T("  ")));
 193     sizer->Add(new wxStaticText(panel,-1,_T("  - - -  STL Diego - - - ")));
 194     sizer->Add(new wxStaticText(panel,-1,_T("  ")));
 195         stlSliderDeltaGauss     = new wxSlider( panel, -1, 100 , 0, 300 , wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS );
 196         stlSliderMarchingCubes= new wxSlider( panel, -1, 128 , 0, 256 , wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS );
 197
 198         sizer->Add(new wxStaticText(panel,-1,_T(" Delta Gauss")));
 199         sizer->Add(stlSliderDeltaGauss);
 200
 201         sizer->Add(new wxStaticText(panel,-1,_T(" Marching Cubes Level")));
 202         sizer->Add(stlSliderMarchingCubes);
 203
 204         stlSliderOpacityInternal = new wxSlider(panel, -1, 100,0,100, wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS);
 205         stlSliderOpacityExternal = new wxSlider(panel, -1, 100,0,100, wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS);
 206
 207         sizer->Add(new wxStaticText(panel, -1,_T(" Opacity STL Internal")));
 208         sizer->Add(stlSliderOpacityInternal);
 209
 210         sizer->Add(new wxStaticText(panel, -1,_T(" Opacity STL External")));
 211         sizer->Add(stlSliderOpacityExternal);
 212
 213         wxButton *btnFileSTL = new wxButton(panel,-1,_T("Generate STL files"));
 214         sizer->Add(btnFileSTL);
 215         sizer->Add(new wxStaticText(panel,-1,_T("  ")));
 216          //-- STL WIDGETS
 217
 218
 219
 220
 221
 222
 223
 224
 225
 226         panel->SetSizer(sizer);
 227         panel->SetAutoLayout(true);
 228         panel->SetSize(400,600);
 229         panel->Layout();
 230
 231
 232         // -- PATHOLOGY CONNECT WIDGETS
 233         Connect(btnP1->GetId(),wxEVT_COMMAND_BUTTON_CLICKED, (wxObjectEventFunction) &wxPathologyWidget_01::OnBtnSetP1);
 234         Connect(btnP2->GetId(),wxEVT_COMMAND_BUTTON_CLICKED, (wxObjectEventFunction) &wxPathologyWidget_01::OnBtnSetP2);
 235         Connect(btnPat->GetId(),wxEVT_COMMAND_BUTTON_CLICKED,(wxObjectEventFunction) &wxPathologyWidget_01::OnBtnSetPat);
 236         Connect(btnExtract->GetId(),wxEVT_COMMAND_BUTTON_CLICKED,(wxObjectEventFunction) &wxPathologyWidget_01::OnBtnExtractPat);
 237         Connect(patSliderOpacity->GetId(), wxEVT_COMMAND_SLIDER_UPDATED , (wxObjectEventFunction) &wxPathologyWidget_01::OnChangePatOpacity);
 238         Connect(patSliderMarchingCubes->GetId(), wxEVT_COMMAND_SLIDER_UPDATED   , (wxObjectEventFunction) &wxPathologyWidget_01::OnChangePatMarchingCubes);
 239         Connect(btnAxis->GetId(), wxEVT_COMMAND_BUTTON_CLICKED  , (wxObjectEventFunction) &wxPathologyWidget_01::OnBtnExtractAxis);
 240
 241
 242         // -- PATHOLOGY CONNECT WIDGETS
 243
 244
 245         // -- STL CONNECT WIDGETS
 246         Connect(btnFileSTL->GetId()                             , wxEVT_COMMAND_BUTTON_CLICKED  , (wxObjectEventFunction) &wxPathologyWidget_01::OnBtnFileSTL           );
 247         Connect(stlSliderDeltaGauss->GetId()    , wxEVT_SCROLL_THUMBRELEASE             , (wxObjectEventFunction) &wxPathologyWidget_01::OnChangeSTLGaussLevel  );
 248         Connect(stlSliderMarchingCubes->GetId() , wxEVT_SCROLL_THUMBRELEASE             , (wxObjectEventFunction) &wxPathologyWidget_01::OnChangeSTLMarchingCubesLevel);
 249         Connect(stlSliderOpacityInternal->GetId(), wxEVT_COMMAND_SLIDER_UPDATED , (wxObjectEventFunction) &wxPathologyWidget_01::OnOpacitySTLInternal);
 250         Connect(stlSliderOpacityExternal->GetId(), wxEVT_COMMAND_SLIDER_UPDATED , (wxObjectEventFunction) &wxPathologyWidget_01::OnOpacitySTLExternal);
 251         // -- STL CONNECT WIDGETS
 252
 253
 254
 255
 256
 257         return panel;
 258 }
 259 //------------------------------------------------------------------------
 260 void wxPathologyWidget_01::Refresh()
 261 {
 262         _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->Render();
 263 }
 264
 265 //------------------------------------------------------------------------
 266 void wxPathologyWidget_01::ConfigureVTK()
 267 {
 268     wxBusyCursor wait;
 269
 270         vtkImageData    *imagedata              = _mar->_experiment->getDynData( )->getVolume( )->castVtk();
 271         //-Maracas-
 272         _maracasSurfaceWidget->ShowMARACASData( _mar );
 273
 274         //CONFIGURACION ADICIONAL
 275     this->ConfigureSTL();
 276         this->ConfigurePathologyExtraction();
 277
 278 }
 279
 280
 281 void wxPathologyWidget_01::ConfigureVTK(vtkImageData *imagedata, int x, int y, int z, double param){
 282
 283         this->px = x;
 284         this->py = y;
 285         this->pz = z;
 286
 287         wxBusyCursor wait;
 288         _maracasSurfaceWidget->ShowMARACASData( _mar );
 289
 290
 291         //CONFIGURACION ADICIONAL
 292         this->ConfigureMPR();
 293         this->ConfigureSTL();
 294         this->ConfigurePathologyExtraction();
 295 }
 296
 297
 298
 299
 300
 301 // ==================================================================================================
 302 // CODIGO NUEVO ADICIONADO POR DIEGO CANTOR
 303 // ==================================================================================================
 304
 305
 306 // ------------------------------------------------------------------------
 307 // START MPR FUNCTIONS - DHC
 308 // ------------------------------------------------------------------------
 309
 310 void wxPathologyWidget_01::ConfigureMPR(){
 311
 312         mprSphereSource = vtkSphereSource::New();
 313     mprMapper = vtkPolyDataMapper::New();
 314     mprActor = vtkActor::New();
 315
 316         mprMapper->SetInput(mprSphereSource->GetOutput());
 317         mprActor->SetMapper(mprMapper);
 318         mprActor->PickableOff( );
 319
 320         mprSphereSource->SetCenter(this->px,this->py,this->pz);
 321         mprSphereSource->SetRadius(2.5);
 322
 323         mprActor->GetProperty()->SetColor( 1.0, 0.0, 0.0 );
 324     vtkRenderer *ren = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();
 325         ren->AddActor(mprActor);
 326
 327 }
 328
 329 // ------------------------------------------------------------------------
 330 // END MPR FUNCTIONS - DHC
 331 // ------------------------------------------------------------------------
 332
 333
 334
 335 // ------------------------------------------------------------------------
 336 // START STL FUNCTIONS - DHC
 337 // ------------------------------------------------------------------------
 338
 339 void wxPathologyWidget_01::ConfigureSTL()
 340 {
 341         stlExterna = vtkPolyData::New();
 342         stlInterna = vtkPolyData::New();
 343
 344         dsm1 = vtkPolyDataMapper ::New();
 345     dsm1->SetInput (stlInterna);
 346     dsm1->ScalarVisibilityOff();
 347
 348     actorInternal = vtkActor::New();
 349     actorInternal->SetMapper (dsm1);
 350     actorInternal->GetProperty()->SetColor (0,1,0);
 351
 352     dsm2 = vtkPolyDataMapper ::New();
 353     dsm2->SetInput (stlExterna);
 354     dsm2->ScalarVisibilityOff();
 355
 356     actorExternal= vtkActor::New();
 357     actorExternal->SetMapper (dsm2);
 358     actorExternal->GetProperty()->SetRepresentationToWireframe();
 359     vtkRenderer *ren = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();
 360     ren->AddActor(actorInternal);
 361     ren->AddActor(actorExternal);
 362
 363         stlExtractor = new vtkSTLExtractor();
 364 }
 365
 366 void wxPathologyWidget_01::generateSTLSurfaces()
 367 {
 368         stlExtractor->setVolume(stlImageData);
 369         stlExtractor->setSigmaLevel(stlDeltaGaussLevel);
 370         stlExtractor->setMarchingCubesLevel(stlMarchingCubesLevel);
 371         stlExtractor->calculate();
 372
 373         stlInterna->DeepCopy(stlExtractor->getInnerSurface());
 374         stlExterna->DeepCopy(stlExtractor->getOuterSurface());
 375 }
 376
 377
 378 void wxPathologyWidget_01::OnOpacitySTLExternal(wxScrollEvent& event){
 379         double value = ((double)stlSliderOpacityExternal->GetValue())/100;
 380     actorExternal->GetProperty( )->SetOpacity( value );
 381         Refresh();
 382 }
 383
 384
 385 void wxPathologyWidget_01::OnOpacitySTLInternal(wxScrollEvent& event){
 386         double value = ((double)stlSliderOpacityInternal->GetValue())/100;
 387     actorInternal->GetProperty( )->SetOpacity( value );
 388         Refresh();
 389 }
 390
 391 void wxPathologyWidget_01::OnBtnFileSTL(wxCommandEvent& event)
 392 {
 393         wxBusyCursor wait;
 394         wxString dirSTL = _mar->_parameters->getStringParam(
 395     marParameters::e_installation_directory );
 396         dirSTL = ( dirSTL == _T("NO_DIRECTORY") ) ? wxGetHomeDir( ) : dirSTL;
 397         wxDirDialog dialog( this, _T("Choose a directory..."), ( !dirSTL.IsEmpty( ) )?
 398         dirSTL: wxGetHomeDir( ) );
 399
 400         if( dialog.ShowModal( ) == wxID_OK )
 401         {
 402
 403
 404                 // ------------------------------------------------------------------------
 405                 //  1.  GENERATE STL FILES
 406                 // ------------------------------------------------------------------------
 407                 const char* fileprefix = "c:\\Creatis\\";
 408                 std::string prefix = fileprefix;
 409
 410
 411                 // 1.1. Se hace un filtro triangular puesto que el stl writer solo recibe poligonos triangulares.
 412
 413         vtkTriangleFilter *filtro = vtkTriangleFilter::New();
 414                 filtro->SetInput(stlInterna);
 415                 vtkPolyDataConnectivityFilter *pdcf = vtkPolyDataConnectivityFilter::New();
 416         pdcf->SetInput( filtro->GetOutput() );
 417         vtkClosePolyData *cpd = vtkClosePolyData::New();
 418         cpd->SetInput( pdcf->GetOutput() );
 419
 420                 // 1.2 se escribe a disco el archivo stl de la superficie interna
 421         cpd->Update();
 422         vtkSTLWriter *writerI = vtkSTLWriter::New();
 423         writerI->SetInput( cpd->GetOutput() );
 424         prefix = fileprefix;
 425         prefix += "internal.stl";
 426         writerI->SetFileName(prefix.c_str());
 427         writerI->SetFileTypeToASCII();
 428         writerI->Write();
 429         writerI->Delete();
 430
 431                 // 1.3 se escribe a disco el archivo stl de la superficie externa
 432                 filtro->SetInput(stlExterna);
 433         cpd->Update();
 434         vtkSTLWriter *writerE = vtkSTLWriter::New();
 435         writerE->SetInput( cpd->GetOutput() );
 436         prefix = fileprefix;
 437         prefix += "external.stl";
 438         writerE->SetFileName( prefix.c_str() );
 439         writerE->SetFileTypeToASCII();
 440         writerE->Write();
 441         writerE->Delete();
 442
 443         filtro->Delete();
 444         cpd->Delete();
 445         pdcf->Delete();
 446         }
 447
 448         //By default *always* update e_installation_directory:
 449         _mar->_parameters->setStringParam( marParameters::e_installation_directory, dialog.GetPath( ) );
 450         _mar->saveParameters( );
 451 }
 452
 453
 454 void wxPathologyWidget_01::OnChangeSTLGaussLevel(wxScrollEvent& event)
 455 {
 456         stlDeltaGaussLevel  = ((double)stlSliderDeltaGauss->GetValue())/100;
 457         generateSTLSurfaces();
 458         Refresh();
 459 }
 460
 461
 462 void wxPathologyWidget_01::OnChangeSTLMarchingCubesLevel(wxScrollEvent& event)
 463 {
 464         stlMarchingCubesLevel = ((double)stlSliderMarchingCubes->GetValue());
 465         generateSTLSurfaces();
 466         Refresh();
 467
 468 }
 469
 470
 471 // ------------------------------------------------------------------------
 472 // END STL FUNCTIONS - DHC
 473 // ------------------------------------------------------------------------
 474
 475
 476 // ------------------------------------------------------------------------
 477 // START PATHOLOGY FUNCTIONS - DHC
 478 // ------------------------------------------------------------------------
 479
 480 double GetDistanciaEuclideana(double* a, double* b){
 481         return sqrt(((a[0]-b[0])*(a[0]-b[0]))+((a[1]-b[1])*(a[1]-b[1]))+((a[2]-b[2])*(a[2]-b[2])));
 482 }
 483
 484
 485 double GetValorPromedio(vtkImageData *data){
 486         int ext[6];
 487         int i, j, k;
 488         int numVoxels = 0;
 489         double promedio = 0;
 490         double valor;
 491
 492
 493         data->GetExtent(ext);
 494
 495
 496         for(i=ext[0];i<=ext[1];i++){
 497                 for(j=ext[2];j<=ext[3];j++){
 498                         for(k=ext[4];k<=ext[5];k++){
 499                                 valor =data->GetScalarComponentAsDouble(i,j,k,0);
 500                                 promedio += valor;
 501                                 numVoxels++;
 502                         }
 503                 }
 504         }
 505         return (double)promedio/(double)numVoxels;
 506 }
 507
 508
 509 void wxPathologyWidget_01::ConfigurePathologyExtraction()
 510 {
 511 }
 512
 513 void wxPathologyWidget_01::generatePathologySurface(){
 514         // 1. Inicializacion
 515 }
 516
 517
 518 void wxPathologyWidget_01::OnBtnSetP1(wxCommandEvent &event){
 519         vtk3DSurfaceWidget *sw = _maracasSurfaceWidget->GetVtk3DSurfaceWidget();
 520
 521         if(p1Actor){
 522                 sw->GetRenderer()->RemoveActor(p1Actor);
 523         }
 524
 525
 526         sw->GetSphereCenter(p1Center);
 527
 528         p1SphereSource = vtkSphereSource::New();
 529     p1Mapper = vtkPolyDataMapper::New();
 530     p1Actor = vtkActor::New();
 531
 532         p1Mapper->SetInput(p1SphereSource->GetOutput());
 533         p1Actor->SetMapper(p1Mapper);
 534         p1Actor->PickableOff( );
 535
 536         p1SphereSource->SetCenter(p1Center[0],p1Center[1],p1Center[2]);
 537         p1SphereSource->SetRadius(0.5);
 538
 539         p1Actor->GetProperty()->SetColor( 0.0, 1.0, 0.0);
 540
 541         sw->GetRenderer()->AddActor(p1Actor);
 542         Refresh();
 543
 544 }
 545
 546 void wxPathologyWidget_01::OnBtnSetP2(wxCommandEvent &event){
 547
 548         vtk3DSurfaceWidget *sw = _maracasSurfaceWidget->GetVtk3DSurfaceWidget();
 549
 550
 551         if(p2Actor){
 552                 sw->GetRenderer()->RemoveActor(p1Actor);
 553         }
 554
 555         sw->GetSphereCenter(p2Center);
 556
 557         p2SphereSource = vtkSphereSource::New();
 558     p2Mapper = vtkPolyDataMapper::New();
 559     p2Actor = vtkActor::New();
 560
 561         p2Mapper->SetInput(p2SphereSource->GetOutput());
 562         p2Actor->SetMapper(p2Mapper);
 563         p2Actor->PickableOff( );
 564
 565         p2SphereSource->SetCenter(p2Center[0],p2Center[1],p2Center[2]);
 566         p2SphereSource->SetRadius(0.5);
 567
 568         p2Actor->GetProperty()->SetColor( 0.0, 1.0, 0.0);
 569
 570         sw->GetRenderer()->AddActor(p2Actor);
 571         Refresh();
 572 }
 573
 574
 575
 576
 577 void wxPathologyWidget_01::OnBtnSetPat(wxCommandEvent &event){
 578
 579         if (p2SphereSource == NULL || p1SphereSource == NULL){
 580                 return;
 581         }
 582
 583         vtk3DSurfaceWidget *sw = _maracasSurfaceWidget->GetVtk3DSurfaceWidget();
 584
 585         if(patActor){
 586                 sw->GetRenderer()->RemoveActor(patActor);
 587         }
 588
 589
 590         double centro[3];
 591
 592
 593         centro[0] = (p1Center[0] + p2Center[0]) / 2;
 594         centro[1] = (p1Center[1] + p2Center[1]) / 2;
 595         centro[2] = (p1Center[2] + p2Center[2]) / 2;
 596
 597         patSphereSource = vtkSphereSource::New();
 598     patMapper = vtkPolyDataMapper::New();
 599     patActor = vtkActor::New();
 600
 601         patMapper->SetInput(patSphereSource->GetOutput());
 602         patActor->SetMapper(patMapper);
 603
 604
 605         patSphereSource->SetCenter(centro[0],centro[1],centro[2]);
 606         patSphereSource->SetRadius(GetDistanciaEuclideana(p1Center,p2Center));
 607
 608         patActor->GetProperty()->SetColor( 0.15, 0.75, 0.15 );
 609         patActor->PickableOff( );
 610         patActor->GetProperty( )->SetOpacity(0.15);
 611
 612         sw->GetRenderer()->AddActor(patActor);
 613         Refresh();
 614 }
 615
 616
 617
 618
 619 void wxPathologyWidget_01::OnBtnExtractPat(wxCommandEvent &event){
 620
 621         //RENDERER
 622         //----------------------------------------------------------------------------------
 623         vtkRenderer* renderer = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();
 624
 625         //IMAGEN ORIGINAL
 626         //----------------------------------------------------------------------------------
 627         vtkImageData* originalData = vtkImageData::New();
 628         originalData->DeepCopy(_mar->_experiment->getDynData()->getVolume()->castVtk());
 629
 630         double spacing[3];
 631         originalData->GetSpacing(spacing);
 632
 633         int extend[6];
 634         originalData->GetExtent(extend);
 635
 636
 637
 638         //==================================================================================
 639         // REGION DE INTERES EN LA IMAGEN ORIGINAL: patImageData
 640         //==================================================================================
 641
 642         double p1crem[3], p2crem[3];
 643
 644         p1crem[0] = p1Center[0] / spacing[0];
 645         p1crem[1] = p1Center[1] / spacing[1];
 646         p1crem[2] = p1Center[2] / spacing[2];
 647
 648         p2crem[0] = p2Center[0] / spacing[0];
 649         p2crem[1] = p2Center[1] / spacing[1];
 650         p2crem[2] = p2Center[2] / spacing[2];
 651
 652
 653
 654
 655
 656     //COORDENADAS DE LA IMAGEN
 657
 658         double radio = GetDistanciaEuclideana(p1crem,p2crem) / 2;
 659         double centro[3];
 660
 661         centro[0] = (p1crem[0] + p2crem[0]) / 2;
 662         centro[1] = (p1crem[1] + p2crem[1]) / 2;
 663         centro[2] = (p1crem[2] + p2crem[2]) / 2;
 664
 665
 666
 667         //Obtener los ejes sobre la esfera para delimitar el area.
 668
 669         double lowX     = centro[0] - (2*radio);
 670         double highX    = centro[0] + (2*radio);
 671         double lowY     = centro[1] - (2*radio);
 672         double highY    = centro[1] + (2*radio);
 673         double lowZ     = centro[2] - (2*radio);
 674         double highZ    = centro[2] + (2*radio);
 675
 676
 677
 678         //Extraer la region de interes:patImageData
 679         //----------------------------------------------------------------------------------
 680
 681         voiFilter = vtkExtractVOI::New();
 682         voiFilter->SetInput(originalData);
 683         voiFilter->SetVOI((int)lowX, (int)highX, (int)lowY, (int)highY, (int)lowZ, (int)highZ);
 684         voiFilter->Update();
 685         patImageData = vtkImageData::New();
 686         patImageData->DeepCopy(voiFilter->GetOutput());
 687 //      patImageData->SetSpacing(spacing);
 688
 689         //==================================================================================
 690         // REGION DE INTERES ESPECIFICA
 691         //==================================================================================
 692
 693         vtkOtsuSphereSource *otsu = new vtkOtsuSphereSource();
 694         thresholdOTSU = otsu->calculateOptimalThreshold(patImageData);
 695
 696
 697         //==================================================================================
 698         // MARCHING CUBES
 699         //==================================================================================
 700
 701         patMCLevel = 128.0;
 702         cubesFilter = vtkMarchingCubes::New();
 703         cubesFilter->SetInput(otsu->getImageForSegmentation());
 704         cubesFilter->SetValue(0,patMCLevel);
 705         cubesFilter->SetNumberOfContours(1);
 706
 707         //==================================================================================
 708         // POLYDATA
 709         //==================================================================================
 710
 711     polyUnico = vtkPolyDataConnectivityFilter::New();
 712         polyUnico->SetInput(cubesFilter->GetOutput());
 713         polyUnico->SetExtractionModeToLargestRegion();
 714
 715         stripper = vtkStripper::New();
 716         stripper->SetInput(polyUnico->GetOutput());
 717
 718
 719         //==================================================================================
 720         // DATA MAPPING
 721         //==================================================================================
 722
 723         dataMapper = vtkPolyDataMapper::New();
 724         dataMapper->SetInput(stripper->GetOutput());
 725         dataMapper->ScalarVisibilityOff();
 726
 727         //==================================================================================
 728         // ACTORES
 729         //==================================================================================
 730
 731
 732         dataActor = vtkActor::New();
 733         dataActor->SetMapper(dataMapper);
 734         dataActor->SetPosition(lowX*spacing[0],lowY*spacing[1],lowZ*spacing[2]);
 735         dataActor->GetProperty()->SetOpacity(patOpacityLevel);
 736         dataActor->GetProperty()->SetColor( 0.20, 0.20, 1.0 );
 737         dataActor->GetProperty()->SetRepresentationToWireframe();
 738
 739         renderer->AddActor(dataActor);
 740
 741
 742
 743         outlineFilter = vtkOutlineFilter::New();
 744     outlineFilter->SetInput(voiFilter->GetOutput());
 745
 746         vtkPolyDataMapper* outlineMapper = vtkPolyDataMapper::New();
 747         outlineMapper->SetInput(outlineFilter->GetOutput());
 748
 749         if(outlineActor){
 750                 renderer->RemoveActor(outlineActor);
 751         }
 752
 753         if(patActor){
 754                 renderer->RemoveActor(patActor);
 755         }
 756
 757         if(p3Actor){
 758                 renderer->RemoveActor(p3Actor);
 759         }
 760
 761
 762         outlineActor = vtkActor::New();
 763         outlineActor->SetMapper(outlineMapper);
 764         outlineActor->GetProperty()->SetOpacity(1);
 765
 766
 767         /*p3SphereSource = vtkSphereSource::New();
 768     p3Mapper = vtkPolyDataMapper::New();
 769     p3Actor = vtkActor::New();
 770
 771         p3Mapper->SetInput(p3SphereSource->GetOutput());
 772         p3Actor->SetMapper(p3Mapper);
 773         p3Actor->PickableOff( );
 774
 775         p3SphereSource->SetCenter(p2Center[0]-(lowX*spacing[0]), p2Center[1]-(lowY*spacing[1]) , p2Center[2]-(lowZ*spacing[2]) );
 776         p3SphereSource->SetRadius(0.5);
 777         p3Actor->GetProperty()->SetColor( 0.9, 0.4, 0.2);*/
 778
 779
 780         //renderer->AddActor(outlineActor);
 781         //renderer->AddActor(p3Actor);
 782
 783
 784         Refresh();
 785 }
 786
 787
 788 void wxPathologyWidget_01::OnChangePatOpacity(wxScrollEvent& event){
 789         patOpacityLevel = ((double)patSliderOpacity->GetValue())/100.0;
 790     dataActor->GetProperty( )->SetOpacity(patOpacityLevel);
 791         Refresh();
 792 }
 793
 794 void wxPathologyWidget_01::OnChangePatMarchingCubes(wxScrollEvent& event){
 795         patMCLevel = ((double)patSliderMarchingCubes->GetValue());
 796         cubesFilter->SetValue(0,patMCLevel);
 797         cubesFilter->Update();
 798         Refresh();
 799 }
 800
 801
 802
 803 //Metodo auxiliar que devuelve la distancia minima de un punto a un polydata.
 804 //Se evalua la distancia del punto a todos los centros de las celdas del polydata y se devuelve la menor
 805 double GetDistanciaAPolydata(double* punto, vtkPoints* centers){
 806         double distancia = 0.0;
 807         double minima = 1e100;
 808         double vertice[3];
 809         int N = centers->GetNumberOfPoints();
 810
 811         for(int i = 0; i < N; i++){
 812                  centers->GetPoint(i,vertice);
 813                  distancia =  sqrt(((punto[0]-vertice[0])*(punto[0]-vertice[0]))+((punto[1]-vertice[1])*(punto[1]-vertice[1]))+((punto[2]-vertice[2])*(punto[2]-vertice[2])));
 814
 815                  if (distancia < minima){
 816                          minima = distancia;
 817                  }
 818         }
 819         return minima;
 820 }
 821
 822
 823 double GetProductoPunto(double* vectorA, double* vectorB){
 824         return (vectorA[0]*vectorB[0] +  vectorA[1]*vectorB[1] + vectorA[2]*vectorB[2]);
 825 }
 826
 827
 828 double GetNormaVector(double* vector){
 829         return sqrt(GetProductoPunto(vector, vector));
 830 }
 831
 832 //construye el vector como la diferencia de dos puntos (vectores en el origen)
 833 void GetVector(double* p0, double* p1, double* resultado){ //p1:destino - p0:origen
 834
 835         resultado[0] = p1[0] - p0[0];
 836         resultado[1] = p1[1] - p0[1];
 837         resultado[2] = p1[2] - p0[2];
 838
 839 }
 840
 841
 842 void VectorUnitario(double *v){
 843         double norma = GetNormaVector(v);
 844         v[0] = v[0] / norma;
 845         v[1] = v[1] / norma;
 846         v[2] = v[2] / norma;
 847 }
 848
 849 //calculo de las normales
 850 void GetNormalCelda(double* p0, double* p1, double* p2, double* normal){
 851
 852         double U[3],V[3];
 853         GetVector(p1,p2,U); //p2-p1
 854         GetVector(p1,p0,V); //p0-p1
 855
 856
 857         normal[0] = (U[1]*V[2])-(U[2]*V[1]);
 858         normal[1] = -((U[0]*V[2])-(U[2]*V[0]));
 859         normal[2] = (U[0]*V[1])-(U[1]*V[0]);
 860         VectorUnitario(normal);
 861 }
 862
 863 void GetNormal(vtkPoints* puntosCelda, double* normal){
 864         double p0[3], p1[3], p2[3];
 865
 866         puntosCelda->GetPoint(0, p0);
 867         puntosCelda->GetPoint(1, p1);
 868         puntosCelda->GetPoint(2, p2);
 869
 870         GetNormalCelda(p0,p1,p2, normal);
 871 }
 872
 873
 874 // Obtiene el ID de la celda mas cercana
 875 int GetClosestCellOnPolyData(double* punto, vtkPoints* centers){
 876         double distancia = 0.0;
 877         double minima = 1e100;
 878         int ID = -1;
 879         double vertice[3];
 880         int N = centers->GetNumberOfPoints();
 881
 882         for(int i = 0; i < N; i++){
 883                  centers->GetPoint(i,vertice);
 884                  distancia =  sqrt(((punto[0]-vertice[0])*(punto[0]-vertice[0]))+((punto[1]-vertice[1])*(punto[1]-vertice[1]))+((punto[2]-vertice[2])*(punto[2]-vertice[2])));
 885
 886                  if (distancia < minima){
 887                          minima = distancia;
 888                          ID = i;
 889                  }
 890         }
 891
 892         //printf("distancia del punto p: (%.2f,%.2f,%.2f) al polydata es %.2f (ID %i)\n", punto[0], punto[1], punto [2], minima, ID);
 893         return ID;
 894 }
 895
 896 double GetProjectionOnNormal(double* vector, double* normal){
 897         double pp = GetProductoPunto(vector, normal);
 898         double vv = GetNormaVector(normal);
 899         double pr = pp / vv;
 900         return pr;
 901 }
 902
 903
 904
 905
 906
 907 void wxPathologyWidget_01::OnBtnExtractAxis(wxCommandEvent& event){
 908
 909
 910
 911         //================================================================================
 912         // RENDERER
 913         //================================================================================
 914         vtkRenderer* renderer = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();
 915
 916         //===================================================================================
 917         // IMAGEN ORIGINAL.
 918         //===================================================================================
 919         int extentOriginal[6];
 920         this->patImageData->GetExtent(extentOriginal);
 921
 922
 923         //===================================================================================
 924         // REMUESTREO
 925         //===================================================================================
 926
 927         int factor = 1;
 928         double resampling = 1.0/(double)factor;
 929
 930
 931
 932         // REMUESTREO DE LA IMAGEN ORIGINAL
 933         //-----------------------------------------------------------------------------------
 934         vtkImageResample* resamplingFilter = vtkImageResample::New();
 935         resamplingFilter->SetAxisOutputSpacing(0,resampling);
 936         resamplingFilter->SetAxisOutputSpacing(1,resampling);
 937         resamplingFilter->SetAxisOutputSpacing(2,resampling);
 938         resamplingFilter->ReleaseDataFlagOff();
 939         resamplingFilter->SetInput(patImageData);
 940         resamplingFilter->SetInterpolationModeToCubic();
 941         resamplingFilter->Update();
 942
 943
 944         //IMAGEN ORIGINAL REMUESTREADA
 945         //-----------------------------------------------------------------------------------
 946         vtkImageData* imagenRemuestreada = resamplingFilter->GetOutput();
 947         imagenRemuestreada->ReleaseDataFlagOff();
 948
 949
 950
 951         //===================================================================================
 952         // IMAGEN BINARIA = CRECIMIENTO DE REGIONES
 953         //===================================================================================
 954
 955         this->isoBinaria = vtkImageData::New();
 956         this->isoBinaria->CopyStructure(resamplingFilter->GetOutput());
 957         this->isoBinaria->SetScalarTypeToDouble();
 958
 959
 960
 961         //CALCULO DE LAS DIMENSIONES Y EL EXTENT
 962         //-----------------------------------------------------------------------------------
 963         int extentBinaria[6];
 964         this->isoBinaria->GetExtent(extentBinaria);
 965         int *dim = isoBinaria->GetDimensions();
 966
 967
 968
 969         // IMAGEN EN NEGRO.
 970         //-----------------------------------------------------------------------------------
 971         int eX, eY, eZ, idpoint; // extent
 972         double *ptr = NULL;
 973
 974         FILE* logger = freopen("tests.txt","w",stdout);
 975
 976
 977
 978         for(eX = extentBinaria[0]; eX <= extentBinaria[1]; eX++){
 979                 for(eY = extentBinaria[2]; eY <= extentBinaria[3]; eY++){
 980                         for (eZ = extentBinaria[4]; eZ <= extentBinaria[5]; eZ++){
 981                                 idpoint = isoBinaria->FindPoint(eX, eY, eZ);
 982                                 ptr = (double *)isoBinaria->GetScalarPointer(eX,eY,eZ); // El apuntador se obtiene sobre el extent.
 983                                 *ptr = 0;
 984                         }
 985                 }
 986         }
 987
 988
 989     //===================================================================================
 990         // POLYDATA
 991         //===================================================================================
 992      vtkPolyData *polyData= stripper->GetOutput();  //OJO SE TOMA SOLO LA SUPERFICIE MAS GRANDE (connectivity + stripper)
 993          //===================================================================================
 994
 995          vtkPoints *puntosPolyData = polyData->GetPoints();
 996          vtkPointData *pointData = polyData->GetPointData();
 997
 998          vtkPolyDataWriter* writePoly = vtkPolyDataWriter::New();
 999          writePoly->SetFileName("poly.vtk");
1000          writePoly->SetInput(polyData);
1001          writePoly->Write();
1002
1003
1004
1005         // PROCESAR LOS VERTICES
1006         //----------------------------------------------------------------------------------
1007
1008          int numVertices = puntosPolyData->GetNumberOfPoints();
1009
1010          // mover el sistema de referencias del polydata para poder determinar lo que se encuentra
1011          //por dentro segun las semillas.
1012          int j;
1013          double vertice[3];
1014          /*for(j=0;j<numVertices;j++){
1015                  puntosPolyData->GetPoint(j,vertice);
1016                  vertice[0] = vertice[0] + extentOriginal[0];
1017                  vertice[1] = vertice[1] + extentOriginal[2];
1018                  vertice[2] = vertice[2] + extentOriginal[4];
1019                  puntosPolyData->SetPoint(j,vertice);
1020          }*/
1021
1022
1023
1024         // PROCESAR LAS CELDAS
1025         //----------------------------------------------------------------------------------
1026
1027
1028          //1. Obtener normales de la superficie (celdas)
1029          vtkPolyDataNormals* polyNormalsFilter = vtkPolyDataNormals::New(); //estas son las normales de los vertices NO de las celdas
1030          polyNormalsFilter->SetInput(polyData);
1031          polyNormalsFilter->SplittingOn();
1032          polyNormalsFilter->ConsistencyOn();
1033          polyNormalsFilter->NonManifoldTraversalOn();
1034          polyNormalsFilter->ComputeCellNormalsOn();
1035          polyNormalsFilter->Update();
1036
1037
1038          vtkPolyData* polyDataConNormales = polyNormalsFilter->GetOutput();
1039
1040
1041          //2. Obtener los centros de las celdas
1042          // Necesitamos el centro de las celdas para trazar el vector interno de un voxel
1043          //a una cara de la superficie.
1044          vtkCellCenters* cellCentersFilter = vtkCellCenters::New();
1045          cellCentersFilter->SetInput(polyDataConNormales);
1046          cellCentersFilter->VertexCellsOn();
1047          cellCentersFilter->Update();
1048
1049          vtkPolyData* polydataCenters = cellCentersFilter->GetOutput();
1050          vtkPoints*   polyCenterPoints  = polydataCenters->GetPoints();
1051
1052          vtkPoints*   puntosCelda = NULL;
1053          vtkIdList*   lista = NULL;
1054          vtkCell*     celda = NULL;
1055      double normal[3];
1056          int numCeldas = polyCenterPoints->GetNumberOfPoints();
1057
1058
1059          vtkPoints *normales = vtkPoints::New();//almacena los componentes cartesianos de la normal de cada celda
1060
1061      int i;
1062      for(i = 0; i < numCeldas ; i++){
1063
1064                 polyCenterPoints->GetPoint(i,vertice);
1065
1066             celda = polyDataConNormales->GetCell(i);
1067                 puntosCelda = celda->GetPoints();
1068         GetNormal(puntosCelda, normal);
1069                 normales->InsertPoint(i, normal);
1070
1071
1072                 lista = celda->GetPointIds();
1073                 for(j = 0; j < lista->GetNumberOfIds(); j++){
1074                         int idPunto = lista->GetId(j);
1075                         puntosCelda->GetPoint(j, vertice);
1076
1077                 }
1078
1079          }
1080
1081
1082
1083          //==========================================================================================
1084          //IMPLEMENTAR EL CRECIMIENTO DE REGIONES
1085          //==========================================================================================
1086
1087         double spacing[3];
1088         patImageData->GetSpacing(spacing);
1089
1090         double p1crem[3], p2crem[3];
1091
1092         p1crem[0] = p1Center[0] / spacing[0];
1093         p1crem[1] = p1Center[1] / spacing[1];
1094         p1crem[2] = p1Center[2] / spacing[2];
1095
1096         p2crem[0] = p2Center[0] / spacing[0];
1097         p2crem[1] = p2Center[1] / spacing[1];
1098         p2crem[2] = p2Center[2] / spacing[2];
1099
1100
1101
1102
1103
1104     //COORDENADAS DE LA IMAGEN
1105
1106         double radio = GetDistanciaEuclideana(p1crem,p2crem) / 2;
1107         double centro[3];
1108
1109         centro[0] = (p1crem[0] + p2crem[0]) / 2;
1110         centro[1] = (p1crem[1] + p2crem[1]) / 2;
1111         centro[2] = (p1crem[2] + p2crem[2]) / 2;
1112
1113
1114
1115         //Obtener los ejes sobre la esfera para delimitar el area.
1116
1117         double lowX             = centro[0] - (2*radio);
1118         double highX    = centro[0] + (2*radio);
1119         double lowY             = centro[1] - (2*radio);
1120         double highY    = centro[1] + (2*radio);
1121         double lowZ             = centro[2] - (2*radio);
1122         double highZ    = centro[2] + (2*radio);
1123
1124         //================================================================================
1125         //0. Inicializacion
1126         //================================================================================
1127
1128         int numeroPuntos = dim[0] * dim[1] * dim[2];
1129
1130
1131         int sem[3];
1132
1133
1134         sem[0] = (int)(p2Center[0]-(lowX*spacing[0]));
1135         sem[1] = (int)(p2Center[1]-(lowY*spacing[1]));
1136         sem[2] = (int)(p2Center[2]-(lowZ*spacing[2]));
1137
1138
1139         int idSemilla = isoBinaria->FindPoint(sem[0],sem[1],sem[2]);
1140
1141         if (idSemilla == -1) return;
1142
1143         int graphSize = (extentBinaria[1]-extentBinaria[0]) * (extentBinaria[3]-extentBinaria[2])*(extentBinaria[5]-extentBinaria[4]);
1144
1145
1146         vtkPriorityQueue *colaEvaluacion = vtkPriorityQueue::New();
1147
1148         colaEvaluacion->Allocate(graphSize);
1149
1150         colaEvaluacion->Insert(0, idSemilla);
1151
1152         vtkFloatingPointType prioridad;
1153
1154         int longitudCola = colaEvaluacion->GetNumberOfItems();
1155
1156
1157         int vecinoID = -1;
1158
1159         int puntoID = -1;
1160
1161         double* coordsVecino;
1162
1163         double* coordsPunto;
1164
1165         int cont=0;
1166
1167         vtkPoints *resultado = vtkPoints::New();
1168
1169         int premiere = 1;
1170
1171
1172         double centroCelda[3], vectorInterno[3], vectorNormal[3];
1173
1174         while(longitudCola > 0){
1175
1176                         //================================================================================
1177                         //1. Procesar el punto
1178                         //================================================================================
1179
1180                     puntoID = colaEvaluacion->Pop(0,prioridad);
1181
1182                         coordsPunto = isoBinaria->GetPoint(puntoID);
1183
1184
1185                     ptr = (double *)isoBinaria->GetScalarPointer((int)(coordsPunto[0]*factor), (int)(coordsPunto[1]*factor), (int)(coordsPunto[2]*factor));
1186                         *ptr = 255.0;
1187
1188                         if(premiere == 1){
1189                                 vtkImageWriter* writeImage = vtkImageWriter::New();
1190                                 writeImage->SetInput(isoBinaria);
1191                                 writeImage->SetFileName("sem.vtk");
1192                                 writeImage->Write();
1193                                 premiere = 0;
1194                         }
1195
1196                         resultado->InsertPoint(cont, coordsPunto);
1197                         cont++;
1198
1199                         //printf("iniciando en el punto (%.1f,%.1f,%.1f) \n", coordsPunto[0], coordsPunto[1], coordsPunto [2]);
1200
1201                         //================================================================================
1202                         //2. Evaluar los vecinos del punto
1203                         //================================================================================
1204                         int vk,vj,vi;
1205                         for(vk = -1; vk<2; vk++) {
1206                                 for(vj = -1; vj<2; vj++) {
1207                                         for(vi = -1; vi<2; vi++) {
1208
1209                                                 vecinoID = puntoID + (vk * dim[1]*dim[0]) + (vj * dim[0]) + vi;
1210
1211                                                 if( vecinoID >= 0 && vecinoID < numeroPuntos && vecinoID != 0 && vecinoID != puntoID ) {
1212
1213                                                                 coordsVecino = isoBinaria->GetPoint(vecinoID);
1214
1215                                                                 //ptr = (unsigned char *)isoBinaria->GetScalarPointer(coordsVecino[0],coordsVecino[1],coordsVecino[2]);
1216                                                                 ptr = (double *)isoBinaria->GetScalarPointer((int)(coordsVecino[0]*factor),(int)(coordsVecino[1]*factor),(int)(coordsVecino[2]*factor));
1217
1218                                                                 //================================================================================
1219                                                                 // Para cada vecino no evaluado evaluar la distancia al polydata
1220                                                                 //================================================================================
1221                                                                 if (*ptr == 0){
1222                                                                         //printf("evaluando punto (%.1f,%.1f,%.1f) \n", coordsVecino[0], coordsVecino[1], coordsVecino [2]);
1223
1224                                                                         int normalID = GetClosestCellOnPolyData(coordsVecino, polyCenterPoints);
1225
1226
1227                                                                         polyCenterPoints->GetPoint(normalID,centroCelda);
1228                                                                         normales->GetPoint(normalID, vectorNormal);
1229
1230                                                                         GetVector(coordsVecino, centroCelda, vectorInterno);
1231                                                                         VectorUnitario(vectorInterno);
1232
1233                                                                         double projection = GetProjectionOnNormal(vectorInterno, vectorNormal);
1234                                                                         ptr = (double *)isoBinaria->GetScalarPointer((int)(coordsVecino[0]*factor),(int)(coordsVecino[1]*factor),(int)(coordsVecino[2]*factor));
1235
1236                                                                         if (projection > 0){
1237
1238                                                                         //================================================================================
1239                                                                         //2.1 Si la proyeccion es positiva
1240                                                                         //colocar el vecino en la cola y poner el punto correspondiente en la imagen a 255.
1241                                                                         // (no ha pasado la frontera, el marcado a 255 se hace arriba en el while, ver paso 1)
1242                                                                         //================================================================================
1243                                                                                 colaEvaluacion->Insert(prioridad,vecinoID);
1244
1245                                                                         } else {
1246
1247                                                                         //================================================================================
1248                                                                         //2.2 Si la proyeccion es negativa con la normal mas cercana del polydata,
1249                                                                         //descartar el vecino == no poner en la cola == alcanzamos la superficie
1250                                                                         //================================================================================
1251
1252                                                                                 *ptr = 0;
1253                                                                         }
1254                                                                 }
1255
1256                                                 }
1257                                         }
1258                                 }
1259                         }
1260
1261                         longitudCola = colaEvaluacion->GetNumberOfItems();
1262         }
1263         colaEvaluacion->Delete();
1264
1265
1266         //================================================================================
1267         //  MAPA DE DISTANCIAS
1268         //================================================================================
1269
1270         this->isoDistance = vtkImageEuclideanDistance::New();
1271         this->isoDistance->SetInput(this->isoBinaria);
1272         this->isoDistance->InitializeOn();
1273         this->isoDistance->Update();
1274         vtkImageData *imDistancias = isoDistance->GetOutput();
1275
1276
1277         //================================================================================
1278         //  CAMINO MINIMO SOBRE EL MAPA DE DISTANCIAS
1279         //================================================================================
1280         this->dijkstraFilter = vtkDijkstraImageData::New();
1281         this->dijkstraFilter->SetInput(isoDistance->GetOutput());
1282
1283         int idA = (int)isoBinaria->FindPoint(p1Center[0],p1Center[1], p1Center[2]);
1284         int idB = (int)isoBinaria->FindPoint(p2Center[0],p2Center[1], p2Center[2]);
1285
1286         dijkstraFilter->SourceID = idA;
1287         dijkstraFilter->SinkID = idB;
1288         dijkstraFilter->Update();
1289
1290         caminoMapper = vtkPolyDataMapper::New();
1291         caminoMapper->SetInput(dijkstraFilter->GetOutput());
1292         caminoMapper->ScalarVisibilityOff();
1293
1294         renderer->RemoveActor(dataActor);
1295
1296         //=======================================================================================================
1297         //RECORRIDO DE LOS PUNTOS DEL EJE DETECTADO
1298         //=======================================================================================================
1299         vtkIdList *camino = dijkstraFilter->GetShortestPathIdList();
1300
1301         int size = camino->GetNumberOfIds();
1302         double coords[3];
1303
1304
1305         for(i=0;i < size;i++){
1306                 int pointId = camino->GetId(i);
1307                 isoBinaria->GetPoint(pointId, coords);
1308                 double scalar = imagenRemuestreada->GetPointData()->GetScalars()->GetTuple1(pointId);
1309
1310                 //----------------------------------------------------------------------------------------
1311                 //Visualizacion del nodo del eje
1312                 //----------------------------------------------------------------------------------------
1313                 vtkSphereSource   *_nodo   = vtkSphereSource::New();
1314                 _nodo->SetCenter(coords[0],coords[1],coords[2]);
1315                 _nodo->SetRadius(0.2);
1316
1317                 vtkPolyDataMapper *_mapper = vtkPolyDataMapper::New();
1318                 _mapper->SetInput(_nodo->GetOutput());
1319
1320
1321         vtkActor          *_actor  = vtkActor::New();
1322                 _actor->GetProperty()->SetColor( 0.8, 0.2, 0.3);
1323             _actor->SetMapper(_mapper);
1324                 _actor->PickableOff( );
1325
1326                 renderer->AddActor(_actor);
1327
1328
1329         }
1330
1331
1332 }
1333 // ------------------------------------------------------------------------
1334 // END PATHOLOGY FUNCTIONS - DHC
1335 // ------------------------------------------------------------------------
1336
1337
1338