lib/maracasVisuLib/src/kernel/vtkDijkstraImageData.cxx

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