lib/maracasVisuLib/src/kernel/include/vtkDijkstraImageData.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 "vtkDijkstraImageData.h"
  27
  28 #include "vtkObjectFactory.h"
  29 #include "vtkIntArray.h"
  30 #include "vtkPointData.h"
  31 #include "vtkPriorityQueue.h"
  32 #include "vtkIdList.h"
  33 #include "vtkImageData.h"
  34 #include "vtkFloatArray.h"
  35  #include <math.h>
  36 #include <stdlib.h>
  37
  38 vtkCxxSetObjectMacro(vtkDijkstraImageData,BoundaryScalars,vtkDataArray);
  39
  40  //----------------------------------------------------------------------------
  41  vtkDijkstraImageData* vtkDijkstraImageData::New()
  42  {
  43
  44
  45
  46
  47    // First try to create the object from the vtkObjectFactory
  48    vtkObject* ret = vtkObjectFactory::CreateInstance("vtkDijkstraImageData");
  49    if(ret)
  50     {
  51      return (vtkDijkstraImageData*)ret;
  52      }
  53    // If the factory was unable to create the object, then create it here.
  54    return new vtkDijkstraImageData;
  55  }
  56
  57
  58
  59
  60 //----------------------------------------------------------------------------
  61 // Constructor sets default values
  62 vtkDijkstraImageData::vtkDijkstraImageData()
  63 {
  64   this->SourceID = 0;
  65   this->SinkID = 0;
  66   this->NumberOfInputPoints = 0;
  67   this->PathPointer = -1;
  68   this->StopWhenEndReached = 1;
  69
  70   this->ShortestPathIdList = NULL;
  71   this->Parent        = NULL;
  72   this->Visited          = NULL;
  73   this->PQ = NULL;
  74   this->BoundaryScalars = NULL;
  75
  76   this->UseInverseDistance = 0;
  77   this->UseInverseSquaredDistance = 0;
  78   this->UseInverseExponentialDistance = 1;
  79   this->UseSquaredDistance = 0;
  80
  81   this->puntos = NULL;
  82   this->lineas = NULL;
  83
  84
  85   this->logger = NULL;
  86
  87         //if ((this->logger = freopen("c:\\creatis\\maracas\\dijkstra.txt","w", stdout)) == NULL)
  88         //      exit(-1);
  89
  90 }
  91
  92
  93 //------------------------------------------------------------------------------
  94 vtkDijkstraImageData::~vtkDijkstraImageData()
  95 {
  96
  97
  98   if (this->ShortestPathIdList)
  99     this->ShortestPathIdList->Delete();
 100   if (this->Parent)
 101     this->Parent->Delete();
 102   if(this->BoundaryScalars)
 103     this->BoundaryScalars->Delete();
 104   //if (this->Visited)
 105   // this->Visited->Delete();
 106   //if (this->PQ)
 107   // this->PQ->Delete();
 108   //if (this->Heap)
 109   //this->Heap->Delete();
 110   //if (this->p)
 111   //this->p->Delete();
 112   //if(this->BoundaryScalars)
 113   //this->BoundaryScalars->Delete();
 114   //DeleteGraph();
 115   //if (this->CumulativeWeightFromSource)
 116   //this->CumulativeWeightFromSource->Delete();
 117
 118 }
 119
 120 //------------------------------------------------------------------------------
 121 unsigned long vtkDijkstraImageData::GetMTime()
 122 {
 123   unsigned long mTime=this->MTime.GetMTime();
 124
 125   return mTime;
 126 }
 127
 128
 129 void vtkDijkstraImageData::SetInput(vtkImageData* data){
 130         this->vtkProcessObject::SetNthInput(0, data);
 131 }
 132
 133 vtkImageData* vtkDijkstraImageData::GetInput(){
 134         return (vtkImageData *)(this->Inputs[0]);
 135 }
 136
 137 //------------------------------------------------------------------------------
 138 void vtkDijkstraImageData::init(vtkImageData *inData)
 139 {
 140
 141   if (this->ShortestPathIdList)
 142     this->ShortestPathIdList->Delete();
 143
 144   if (this->Parent)
 145     this->Parent->Delete();
 146   if (this->Visited)
 147     this->Visited->Delete();
 148   if (this->PQ)
 149   this->PQ->Delete();
 150
 151   this->ShortestPathIdList = vtkIdList::New();
 152   this->Parent        = vtkIntArray::New();
 153   this->Visited          = vtkIntArray::New();
 154   this->PQ = vtkPriorityQueue::New();
 155
 156
 157   CreateGraph(inData);
 158
 159
 160   int numPoints = inData->GetNumberOfPoints();
 161
 162   this->Parent->SetNumberOfComponents(1);
 163   this->Parent->SetNumberOfTuples(numPoints);
 164   this->Visited->SetNumberOfComponents(1);
 165   this->Visited->SetNumberOfTuples(numPoints);
 166
 167 }
 168
 169 void vtkDijkstraImageData::DeleteGraph()
 170 {
 171
 172   /*const int npoints = this->GetNumberOfInputPoints();
 173
 174     if (this->Neighbors)
 175     {
 176     for (int i = 0; i < npoints; i++)
 177     {
 178     if(this->Neighbors[i])
 179     this->Neighbors[i]->Delete();
 180     }
 181     delete [] this->Neighbors;
 182     }
 183     this->Neighbors = NULL;
 184   */
 185 }
 186
 187 //------------------------------------------------------------------------------
 188 void vtkDijkstraImageData::CreateGraph(vtkImageData *inData) {
 189
 190
 191   //DeleteGraph();
 192   //delete old arrays in case we are re-executing this filter
 193   int numPoints = inData->GetNumberOfPoints();
 194   this->SetNumberOfInputPoints(numPoints);
 195
 196   // initialization
 197   int *dim = inData->GetDimensions();
 198   vtkDataArray *scalars = inData->GetPointData()->GetScalars();
 199   vtkIdList *graphNodes = vtkIdList::New();
 200   int graphSize = 0;
 201   // create the graph
 202   for(int k = 0; k <dim[2]; k++) {
 203     this->UpdateProgress ((float) k / (2 * ((float) dim[2] - 1)));
 204     for(int j = 0; j <dim[1]; j++) {
 205       for(int i = 0; i <dim[0]; i++) {
 206
 207     int id = k*(dim[1]*dim[0]) + j*dim[0] + i;
 208     float maskValue = scalars->GetTuple1(id);
 209     // only add neighbor if it is in the graph
 210
 211     if(maskValue > 0) {
 212       // add to graph
 213       graphNodes->InsertNextId(id);
 214       graphSize++;
 215     }
 216       }
 217     }
 218   }
 219
 220   this->SetNumberOfGraphNodes(graphSize);
 221   //printf("graph size %i \n ",graphSize);
 222
 223   // fill the PQ
 224   PQ->Allocate(graphSize);
 225   for(int i=0; i<graphSize;i++) {
 226     PQ->Insert(VTK_LARGE_FLOAT,graphNodes->GetId(i));
 227   }
 228   // free some memory
 229   graphNodes->Delete();
 230
 231 }
 232
 233
 234 //------------------------------------------------------------------------------
 235 void vtkDijkstraImageData::RunDijkstra(vtkDataArray *scalars,int startv, int endv)
 236 {
 237
 238         int i, u, v;
 239
 240         printf("tipo de ponderacion de pesos : linear %i, squared %i, exponential %i\n",this->UseInverseDistance,this->UseInverseSquaredDistance,this->UseInverseExponentialDistance);
 241
 242
 243
 244         //-------------------------------------------------------------------------------
 245         // INICIALIZACION.
 246         //-------------------------------------------------------------------------------
 247
 248         // 1. Se marcan todos los vertices del grafo como no visitados.
 249         InitSingleSource(startv);
 250
 251         // 2. Se marca el punto de arranque como visitado
 252         this->Visited->SetValue(startv, 1);
 253
 254         // 3. Se obtiene el tamanio inicial de los vertices no visitados
 255         int initialSize = PQ->GetNumberOfItems();
 256         int size = initialSize;
 257         int stop = 0;
 258
 259         //Nota: PQ contiene los nodos del grafo no evaluados.
 260         //Segun Dijkstra, el algoritmo termina cuando todos los nodos tienen un valor de distancia minima.
 261
 262         //------------------------------------------------------------------------
 263         // PARA TODOS LOS VERTICES EN PQ (SIN EVALUAR)
 264         //------------------------------------------------------------------------
 265         while ((PQ->GetNumberOfItems() > 0) && !stop)
 266         {
 267
 268                 this->UpdateProgress (0.5 + (float) (initialSize - size) / ((float) 2 * initialSize));
 269
 270
 271                 vtkFloatingPointType u_weight;
 272
 273                 //----------------------------------------------------------------------
 274                 //1. Aqui obtiene el siguiente punto a evaluar. Al retirarlo de la cola
 275                 // se asume que ya ha sido evaluado, se evaluaran sus vecinos.
 276                 // u is now in S since the shortest path to u is determined
 277                 //----------------------------------------------------------------------
 278                 #if (VTK_MAJOR_VERSION == 4 && VTK_MINOR_VERSION == 0)
 279                 u = PQ->Pop(u_weight,0);
 280                 #else
 281                 u = PQ->Pop(0, u_weight);
 282                 #endif
 283                 //printf("Evaluando el nodo %i con prioridad %.2f\n",u,u_weight);
 284                 //printPointData(u);
 285
 286
 287         //-----------------------------------------
 288                 //2. Se marca "u" como ya visitado.
 289                 //-----------------------------------------
 290                 this->Visited->SetValue(u, 1);
 291
 292
 293
 294                 //------------------------------------------------------------------------------------
 295                 //3. Si u es el nodo final el algoritmo se detiene. (despues de evaluar sus vecinos)
 296                 //------------------------------------------------------------------------------------
 297                 if (u == endv && StopWhenEndReached)
 298                         stop = 1;
 299
 300                 // Update all vertices v neighbors to u
 301                 // find the neighbors of u
 302
 303
 304
 305                 //-----------------------------------------
 306                 //4. Encontrar los vecinos de u
 307                 //-----------------------------------------
 308                 vtkIdList *list = vtkIdList::New();
 309                 this->FindNeighbors(list,u,scalars);
 310
 311                 //printf("%i tiene %i vecinos\n",u,list->GetNumberOfIds());
 312
 313                 //-----------------------------------------
 314                 // 5. Para cada vecino ....
 315                 //-----------------------------------------
 316                 for (i = 0; i < list->GetNumberOfIds(); i++)
 317                 {
 318
 319                         //---------------------------------------------------------------------
 320                         // 5.1 Se obtiene su ID (el indice es el orden en la lista de vecinos)
 321                         //---------------------------------------------------------------------
 322                         v = list->GetId(i);
 323                         //printf("---> evaluando el vecino % i\n",v);
 324
 325                         // s is the set of vertices with determined shortest path...do not
 326                         // use them again
 327
 328                         //-----------------------------------------------------------------
 329                         // 5.2 Si ya ha sido visitado se descarta...
 330                         //-----------------------------------------------------------------
 331                         if (this->Visited->GetValue(v) != 1)
 332                         {
 333                                 //printf("---> ---> el vecino %i no ha sido visitado\n",v);
 334                                 // Only relax edges where the end is not in s and edge
 335                                 // is in the front
 336                                 //---------------------------------------------------
 337                                 // 5.2.1 Se obtiene el COSTO W
 338                                 //---------------------------------------------------
 339                                 float w = EdgeCost(scalars, u, v);
 340
 341                                 //---------------------------------------------------
 342                                 // 5.2.2
 343                                 //---------------------------------------------------
 344                                 float v_weight = this->PQ->GetPriority(v);
 345                                 if(v == endv) {
 346                                         //printf("we have found endv %i, its weight is %f, neighbor is %i , weight is %f, edge weight is %f\n",v,v_weight,u,u_weight,w);
 347                                         //stop=1;
 348                                         //break;
 349                                 }
 350
 351                                 //printf("---> ---> el costo de %i es %.2f\n",v,(w+u_weight));
 352                                 //printf("---> ---> el costo previo de %i es %.2f\n",v,v_weight);
 353
 354
 355                                 //-------------------------------------------------------------------------------
 356                                 // 5.2.3 Si se encuentra un camino menor, se reajusta el costo (Node Relaxation)
 357                                 //-------------------------------------------------------------------------------
 358                                 if (v_weight > (u_weight + w))
 359                                 {
 360                                         // 5.2.3.1 Se elimina V de la cola de nodos por evaluar (PQ). Porque se encontro un costo mas bajo
 361                                         this->PQ->DeleteId(v);
 362
 363                                         // 5.2.3.2 Se inserta nuevamente V en la cola de nodos por evaluar con el nuevo costo.
 364                                         this->PQ->Insert((u_weight + w),v);
 365                                         //printf("---> ---> reajustando el peso de %i a %.2f\n",v,u_weight + w);
 366                                         //printf("u_weight=%.2f\n",u_weight);
 367
 368                                         // 5.2.3.3 se coloca V conectado con U (esto sirve para reconstruir el camino minimo)
 369
 370                                         this->Parent->SetValue(v, u);
 371                                         //printf("setting parent of %i to be %i",v,u);
 372                                 }
 373                         }
 374                         else{
 375                                 //printf("---> el vecino %i ya se visito\n",v);
 376                         }
 377                 }
 378
 379                 //-------------------------------------------------------
 380                 // 6. Liberar la memoria ocupada por la lista de vecinos
 381                 //-------------------------------------------------------
 382                 list->Delete();
 383                 size--;
 384         }
 385         this->PQ->Delete();
 386         this->Visited->Delete();
 387 }
 388
 389 //----------------------------------------------------------------------------
 390 void vtkDijkstraImageData::InitSingleSource(int startv)
 391 {
 392   for (int v = 0; v < this->GetNumberOfInputPoints(); v++)
 393     {
 394       this->Parent->SetValue(v, -1);
 395       this->Visited->SetValue(v, 0);
 396     }
 397   PQ->DeleteId(startv);
 398   // re-insert the source with priority 0
 399   PQ->Insert(0,startv);
 400   //printf("priority of startv %f",PQ->GetPriority(startv));
 401 }
 402
 403
 404 //----------------------------------------------------------------------------
 405 void vtkDijkstraImageData::FindNeighbors(vtkIdList *list,int id, vtkDataArray *scalars) {
 406
 407   // find i, j, k for that node
 408   vtkImageData *input = this->GetInput();
 409
 410   int *dim = input->GetDimensions();
 411   int numPts = dim[0] * dim[1] * dim[2];
 412
 413   //printf("vecinos de %i :(",id);
 414   for(int vk = -1; vk<2; vk++) {
 415     for(int vj = -1; vj<2; vj++) {
 416       for(int vi = -1; vi<2; vi++) {
 417
 418                 int tmpID = id + (vk * dim[1]*dim[0]) + (vj * dim[0]) + vi;
 419                 // check we are in bounds (for volume faces)
 420                 if( tmpID >= 0 && tmpID < numPts && tmpID != 0) {
 421                         float mask = scalars->GetTuple1(tmpID);
 422                         // only add neighbor if it is in the graph
 423                         if(mask > 0 && tmpID != id) {
 424                                 list->InsertUniqueId(tmpID);
 425                                 //printf("%i,",tmpID);
 426                                 //printPointData(tmpID);
 427                         }
 428                 }
 429       }
 430     }
 431   }
 432   //printf(")\n");
 433 }
 434
 435
 436
 437 float vtkDijkstraImageData::fuerzaAtraccion(int u, float w)
 438 {
 439         //1. Identificar las coordenadas del punto de inicio y del punto final (polos)
 440         double coordsSource[3];
 441         double coordsSink[3];
 442         double coordsPoint[3];
 443
 444         this->GetInput()->GetPoint(this->GetSourceID(), coordsSource);
 445         this->GetInput()->GetPoint(this->GetSinkID(), coordsSink);
 446
 447
 448         //2. Identificar las coordenadas del punto al cual se le quiere sacar el potencial electrico
 449         this->GetInput()->GetPoint(u, coordsPoint);
 450
 451
 452         //3 Calcular r1 y r2
 453
 454         float r1 = sqrt((coordsSource[0]-coordsPoint[0])*(coordsSource[0]-coordsPoint[0]) +
 455                 (coordsSource[1]-coordsPoint[1])*(coordsSource[1]-coordsPoint[1]) +
 456                 (coordsSource[2]-coordsPoint[2])*(coordsSource[2]-coordsPoint[2]));
 457
 458
 459         float r2 = sqrt((coordsSink[0]-coordsPoint[0])*(coordsSink[0]-coordsPoint[0]) +
 460                 (coordsSink[1]-coordsPoint[1])*(coordsSink[1]-coordsPoint[1]) +
 461                 (coordsSink[2]-coordsPoint[2])*(coordsSink[2]-coordsPoint[2]));
 462
 463         float d = sqrt((coordsSink[0]-coordsSource[0])*(coordsSink[0]-coordsSource[0]) +
 464                 (coordsSink[1]-coordsSource[1])*(coordsSink[1]-coordsSource[1]) +
 465                 (coordsSink[2]-coordsSource[2])*(coordsSink[2]-coordsSource[2]));
 466
 467
 468
 469         if (r2 < d && r1 > d){
 470                 return w/2;
 471         } else if ( r2 > d && r1 < d){
 472                 return w/2;
 473         } else if (r1 < d && r2 < d){
 474                 return w*w;
 475         } else return w/2;
 476
 477
 478
 479 }
 480
 481 //----------------------------------------------------------------------------
 482 // The edge cost function should be implemented as a callback function to
 483 // allow more advanced weighting
 484 float vtkDijkstraImageData::EdgeCost(vtkDataArray *scalars, int u, int v)
 485 {
 486
 487   float w;
 488
 489     float dist2 = scalars->GetTuple1(v);  //+ fuerzaAtraccion(v, scalars->GetTuple1(v));
 490     float dist = sqrt(dist2);
 491     if(this->UseInverseDistance)
 492       w = (1.0/dist);
 493     else if(this->UseInverseSquaredDistance)
 494       w = (1.0/(dist2));
 495     else if(this->UseInverseExponentialDistance)
 496       w = (1.0/exp(dist));
 497     else if(this->UseSquaredDistance)
 498       w = dist2;
 499
 500
 501
 502   return w;
 503 }
 504
 505 //----------------------------------------------------------------------------
 506 void vtkDijkstraImageData::BuildShortestPath(int start,int end)
 507 {
 508
 509         int p = end;
 510         int next = 0;
 511         int i = 0;
 512         double punto[3];
 513
 514         this->puntos = vtkPoints::New();
 515         this->lineas = vtkCellArray::New();
 516
 517
 518
 519         //printf("======================camino minimo========================\n");
 520         //printf("Punto Inicial = %i,  Punto Final = %i\n",start, end);
 521
 522         while (p != start && p > 0)
 523     {
 524                 printPointData(p);
 525
 526                 this->GetInput()->GetPoint(p,punto);
 527                 this->puntos->InsertPoint(i,punto); //puntos y lineas usados como resultado (vtkPolyData)
 528
 529                 next = this->Parent->GetValue(p);
 530
 531                 this->GetInput()->GetPoint(next,punto);
 532                 this->puntos->InsertPoint(i+1,punto);
 533
 534                 this->lineas->InsertNextCell(2);
 535                 this->lineas->InsertCellPoint(i);
 536                 this->lineas->InsertCellPoint(i+1);
 537
 538                 this->ShortestPathIdList->InsertNextId(p);
 539                 p = next;
 540                 i++;
 541     }
 542
 543         this->ShortestPathIdList->InsertNextId(p);
 544         printPointData(p);
 545         //printf("===========================================================\n");
 546
 547
 548
 549         this->GetOutput()->SetPoints (this->puntos);
 550         this->GetOutput()->SetLines(this->lineas);
 551         this->GetOutput()->Modified();
 552     //fclose(logger);
 553
 554 }
 555
 556 vtkIdList* vtkDijkstraImageData::GetShortestPathIdList()
 557 {
 558         return this->ShortestPathIdList;
 559 }
 560
 561
 562
 563 void vtkDijkstraImageData::printPointData(int pointID)
 564 {
 565         double coords[3];
 566
 567         this->GetInput()->GetPoint(pointID, coords);
 568         double scalar = this->GetInput()->GetPointData()->GetScalars()->GetTuple1(pointID);
 569
 570         printf("Punto ID: %i ",pointID);
 571         printf("Coords[%.0f, %.0f, %.0f] ", coords[0], coords[1], coords[2]);
 572         printf("Scalar: %.2f\n", scalar);
 573
 574 }
 575
 576
 577
 578 //----------------------------------------------------------------------------
 579 // ITERATOR PART
 580
 581 void vtkDijkstraImageData::InitTraversePath(){
 582   this->PathPointer = -1;
 583 }
 584
 585 //----------------------------------------------------------------------------
 586 int vtkDijkstraImageData::GetNumberOfPathNodes(){
 587   return this->ShortestPathIdList->GetNumberOfIds();
 588 }
 589
 590 //----------------------------------------------------------------------------
 591 int vtkDijkstraImageData::GetNextPathNode(){
 592   this->PathPointer = this->PathPointer + 1;
 593
 594   if(this->PathPointer < this->GetNumberOfPathNodes()) {
 595     //printf("this->GetPathNode(this->PathPointer) %i",this->GetPathNode(this->PathPointer));
 596     return this->ShortestPathIdList->GetId(this->PathPointer);
 597   } else {
 598     return -1;
 599   }
 600 }
 601
 602 //----------------------------------------------------------------------------
 603 // find closest scalar to id that is non-zero
 604 int vtkDijkstraImageData::findClosestPointInGraph(vtkDataArray *scalars,int id,int dim0,int dim1, int dim2) {
 605
 606
 607   int kFactor = dim0 * dim1;
 608   int jFactor = dim0;
 609
 610   int numPoints = kFactor * dim2;
 611   vtkIdList* Q = vtkIdList::New();
 612   Q->InsertNextId(id);
 613
 614   int pointer = 0;
 615   int foundID = -1;
 616
 617   while(Q->GetNumberOfIds() != 0) {
 618
 619     int current = Q->GetId(pointer);
 620     pointer = pointer + 1;
 621     //printf("we are looking at id %i \n",current);
 622
 623     // check to see if we found something in the graph
 624     if(scalars->GetTuple1(current) >0) {
 625       //printf("before return");
 626       return current;
 627     } else {
 628       // set it to -1 to show that we already looked at it
 629       scalars->SetTuple1(current,-1);
 630       // put the neighbors on the stack
 631       // top
 632       if (current + kFactor <numPoints) {
 633     if(scalars->GetTuple1(current + kFactor) != -1){
 634       //printf("expand k+1 %i",current + kFactor);
 635       Q->InsertNextId(current+kFactor);
 636     }
 637       }
 638       // bottom
 639        if (current - kFactor >= 0){
 640      if(scalars->GetTuple1(current - kFactor) != -1){
 641        //printf("expand k-1 %i", current - kFactor);
 642        Q->InsertNextId(current-kFactor);
 643      }
 644        }
 645        // front
 646        if (current + jFactor < numPoints) {
 647      if(scalars->GetTuple1(current + jFactor) != -1){
 648        //printf("expand j+1 %i",current + jFactor);
 649        Q->InsertNextId(current + jFactor);
 650      }
 651        }
 652        // back
 653        if (current - jFactor >= 0) {
 654      if(scalars->GetTuple1(current - jFactor) != -1){
 655        //printf("expand j-1 %i",current - jFactor);
 656        Q->InsertNextId(current - jFactor);
 657      }
 658        }
 659        // left
 660        if (current+1 <numPoints){
 661      if(scalars->GetTuple1(current + 1) != -1){
 662        //printf("expand i+1 %i",current+1);
 663        Q->InsertNextId(current + 1);
 664      }
 665        }
 666        // right
 667        if (current -1 >= 0) {
 668      if(scalars->GetTuple1(current - 1) != -1){
 669        //printf("expand i-1 %i",current - 1);
 670        Q->InsertNextId(current - 1);
 671      }
 672        }
 673      }
 674    }
 675    Q->Delete();
 676    return foundID;
 677  }
 678
 679
 680
 681  //----------------------------------------------------------------------------
 682  // This method is passed a input and output Data, and executes the filter
 683  // algorithm to fill the output from the input.
 684  // It just executes a switch statement to call the correct function for
 685  // the Datas data types.
 686  // -- sp 2002-09-05 updated for vtk4
 687  void vtkDijkstraImageData::Execute()
 688  {
 689
 690    vtkImageData *inData = this->GetInput();
 691    vtkPolyData  *outData = this->GetOutput();
 692
 693
 694
 695
 696    // Components turned into x, y and z
 697    if (inData->GetNumberOfScalarComponents() > 3)
 698      {
 699      vtkErrorMacro("This filter can handle upto 3 components");
 700      return;
 701      }
 702
 703
 704    //printf("*************** vtkDijkstraImageData Execute **************\n\n\n");
 705    this->init(inData);
 706    //printf("*************** after init ****************\n\n\n");
 707
 708    int *dim = inData->GetDimensions();
 709    vtkDataArray *scalars = inData->GetPointData()->GetScalars();
 710
 711     // find closest point in graph to source and sink if their value is not 0
 712    printf("source ID %i value is %f\n",this->GetSourceID(),scalars->GetTuple1(this->GetSourceID()));
 713    if(scalars->GetTuple1(this->GetSourceID()) == 0) {
 714
 715      vtkFloatArray *copyScalars = vtkFloatArray::New();
 716      copyScalars->DeepCopy(inData->GetPointData()->GetScalars());
 717
 718      this->SetSourceID(this->findClosestPointInGraph(copyScalars,this->GetSourceID(),dim[0],dim[1],dim[2]));
 719
 720      copyScalars->Delete();
 721      printf("NEW source ID %i value is %f\n",this->GetSourceID(),scalars->GetTuple1(this->GetSourceID()));
 722
 723    }
 724
 725    printf("sink ID %i value is %f\n",this->GetSinkID(),scalars->GetTuple1(this->GetSinkID()));
 726    if(scalars->GetTuple1(this->GetSinkID()) == 0) {
 727      vtkFloatArray *copyScalars2 = vtkFloatArray::New();
 728      copyScalars2->DeepCopy(inData->GetPointData()->GetScalars());
 729      this->SetSinkID(this->findClosestPointInGraph(copyScalars2,this->GetSinkID(),dim[0],dim[1],dim[2]));
 730      copyScalars2->Delete();
 731
 732      printf("NEW sink ID %i value is %f\n",this->GetSinkID(),scalars->GetTuple1(this->GetSinkID()));
 733    }
 734
 735
 736    this->RunDijkstra(scalars,this->GetSourceID(),this->GetSinkID());
 737
 738    this->BuildShortestPath(this->GetSourceID(),this->GetSinkID());
 739
 740
 741
 742
 743  }
 744
 745
 746  //----------------------------------------------------------------------------
 747  void vtkDijkstraImageData::PrintSelf(ostream& os, vtkIndent indent)
 748  {
 749    Superclass::PrintSelf(os,indent);
 750
 751    os << indent << "Source ID: ( "
 752       << this->GetSourceID() << " )\n";
 753
 754    os << indent << "Sink ID: ( "
 755       << this->GetSinkID() << " )\n";
 756  }
 757
 758
 759