1 /*# ---------------------------------------------------------------------
3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
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
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.
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
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 # ------------------------------------------------------------------------ */
26 #include "vtkDijkstraImageData.h"
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"
38 vtkCxxSetObjectMacro(vtkDijkstraImageData,BoundaryScalars,vtkDataArray);
40 //----------------------------------------------------------------------------
41 vtkDijkstraImageData* vtkDijkstraImageData::New()
47 // First try to create the object from the vtkObjectFactory
48 vtkObject* ret = vtkObjectFactory::CreateInstance("vtkDijkstraImageData");
51 return (vtkDijkstraImageData*)ret;
53 // If the factory was unable to create the object, then create it here.
54 return new vtkDijkstraImageData;
60 //----------------------------------------------------------------------------
61 // Constructor sets default values
62 vtkDijkstraImageData::vtkDijkstraImageData()
66 this->NumberOfInputPoints = 0;
67 this->PathPointer = -1;
68 this->StopWhenEndReached = 1;
70 this->ShortestPathIdList = NULL;
74 this->BoundaryScalars = NULL;
76 this->UseInverseDistance = 0;
77 this->UseInverseSquaredDistance = 0;
78 this->UseInverseExponentialDistance = 1;
79 this->UseSquaredDistance = 0;
87 //if ((this->logger = freopen("c:\\creatis\\maracas\\dijkstra.txt","w", stdout)) == NULL)
93 //------------------------------------------------------------------------------
94 vtkDijkstraImageData::~vtkDijkstraImageData()
98 if (this->ShortestPathIdList)
99 this->ShortestPathIdList->Delete();
101 this->Parent->Delete();
102 if(this->BoundaryScalars)
103 this->BoundaryScalars->Delete();
105 // this->Visited->Delete();
107 // this->PQ->Delete();
109 //this->Heap->Delete();
112 //if(this->BoundaryScalars)
113 //this->BoundaryScalars->Delete();
115 //if (this->CumulativeWeightFromSource)
116 //this->CumulativeWeightFromSource->Delete();
120 //------------------------------------------------------------------------------
121 unsigned long vtkDijkstraImageData::GetMTime()
123 unsigned long mTime=this->MTime.GetMTime();
129 void vtkDijkstraImageData::SetInput(vtkImageData* data){
130 this->vtkProcessObject::SetNthInput(0, data);
133 vtkImageData* vtkDijkstraImageData::GetInput(){
134 return (vtkImageData *)(this->Inputs[0]);
137 //------------------------------------------------------------------------------
138 void vtkDijkstraImageData::init(vtkImageData *inData)
141 if (this->ShortestPathIdList)
142 this->ShortestPathIdList->Delete();
145 this->Parent->Delete();
147 this->Visited->Delete();
151 this->ShortestPathIdList = vtkIdList::New();
152 this->Parent = vtkIntArray::New();
153 this->Visited = vtkIntArray::New();
154 this->PQ = vtkPriorityQueue::New();
160 int numPoints = inData->GetNumberOfPoints();
162 this->Parent->SetNumberOfComponents(1);
163 this->Parent->SetNumberOfTuples(numPoints);
164 this->Visited->SetNumberOfComponents(1);
165 this->Visited->SetNumberOfTuples(numPoints);
169 void vtkDijkstraImageData::DeleteGraph()
172 /*const int npoints = this->GetNumberOfInputPoints();
176 for (int i = 0; i < npoints; i++)
178 if(this->Neighbors[i])
179 this->Neighbors[i]->Delete();
181 delete [] this->Neighbors;
183 this->Neighbors = NULL;
187 //------------------------------------------------------------------------------
188 void vtkDijkstraImageData::CreateGraph(vtkImageData *inData) {
192 //delete old arrays in case we are re-executing this filter
193 int numPoints = inData->GetNumberOfPoints();
194 this->SetNumberOfInputPoints(numPoints);
197 int *dim = inData->GetDimensions();
198 vtkDataArray *scalars = inData->GetPointData()->GetScalars();
199 vtkIdList *graphNodes = vtkIdList::New();
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++) {
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
213 graphNodes->InsertNextId(id);
220 this->SetNumberOfGraphNodes(graphSize);
221 //printf("graph size %i \n ",graphSize);
224 PQ->Allocate(graphSize);
225 for(int i=0; i<graphSize;i++) {
226 PQ->Insert(VTK_LARGE_FLOAT,graphNodes->GetId(i));
229 graphNodes->Delete();
234 //------------------------------------------------------------------------------
235 void vtkDijkstraImageData::RunDijkstra(vtkDataArray *scalars,int startv, int endv)
240 printf("tipo de ponderacion de pesos : linear %i, squared %i, exponential %i\n",this->UseInverseDistance,this->UseInverseSquaredDistance,this->UseInverseExponentialDistance);
244 //-------------------------------------------------------------------------------
246 //-------------------------------------------------------------------------------
248 // 1. Se marcan todos los vertices del grafo como no visitados.
249 InitSingleSource(startv);
251 // 2. Se marca el punto de arranque como visitado
252 this->Visited->SetValue(startv, 1);
254 // 3. Se obtiene el tamanio inicial de los vertices no visitados
255 int initialSize = PQ->GetNumberOfItems();
256 int size = initialSize;
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.
262 //------------------------------------------------------------------------
263 // PARA TODOS LOS VERTICES EN PQ (SIN EVALUAR)
264 //------------------------------------------------------------------------
265 while ((PQ->GetNumberOfItems() > 0) && !stop)
268 this->UpdateProgress (0.5 + (float) (initialSize - size) / ((float) 2 * initialSize));
271 vtkFloatingPointType u_weight;
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);
281 u = PQ->Pop(0, u_weight);
283 //printf("Evaluando el nodo %i con prioridad %.2f\n",u,u_weight);
287 //-----------------------------------------
288 //2. Se marca "u" como ya visitado.
289 //-----------------------------------------
290 this->Visited->SetValue(u, 1);
294 //------------------------------------------------------------------------------------
295 //3. Si u es el nodo final el algoritmo se detiene. (despues de evaluar sus vecinos)
296 //------------------------------------------------------------------------------------
297 if (u == endv && StopWhenEndReached)
300 // Update all vertices v neighbors to u
301 // find the neighbors of u
305 //-----------------------------------------
306 //4. Encontrar los vecinos de u
307 //-----------------------------------------
308 vtkIdList *list = vtkIdList::New();
309 this->FindNeighbors(list,u,scalars);
311 //printf("%i tiene %i vecinos\n",u,list->GetNumberOfIds());
313 //-----------------------------------------
314 // 5. Para cada vecino ....
315 //-----------------------------------------
316 for (i = 0; i < list->GetNumberOfIds(); i++)
319 //---------------------------------------------------------------------
320 // 5.1 Se obtiene su ID (el indice es el orden en la lista de vecinos)
321 //---------------------------------------------------------------------
323 //printf("---> evaluando el vecino % i\n",v);
325 // s is the set of vertices with determined shortest path...do not
328 //-----------------------------------------------------------------
329 // 5.2 Si ya ha sido visitado se descarta...
330 //-----------------------------------------------------------------
331 if (this->Visited->GetValue(v) != 1)
333 //printf("---> ---> el vecino %i no ha sido visitado\n",v);
334 // Only relax edges where the end is not in s and edge
336 //---------------------------------------------------
337 // 5.2.1 Se obtiene el COSTO W
338 //---------------------------------------------------
339 float w = EdgeCost(scalars, u, v);
341 //---------------------------------------------------
343 //---------------------------------------------------
344 float v_weight = this->PQ->GetPriority(v);
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);
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);
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))
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);
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);
368 // 5.2.3.3 se coloca V conectado con U (esto sirve para reconstruir el camino minimo)
370 this->Parent->SetValue(v, u);
371 //printf("setting parent of %i to be %i",v,u);
375 //printf("---> el vecino %i ya se visito\n",v);
379 //-------------------------------------------------------
380 // 6. Liberar la memoria ocupada por la lista de vecinos
381 //-------------------------------------------------------
386 this->Visited->Delete();
389 //----------------------------------------------------------------------------
390 void vtkDijkstraImageData::InitSingleSource(int startv)
392 for (int v = 0; v < this->GetNumberOfInputPoints(); v++)
394 this->Parent->SetValue(v, -1);
395 this->Visited->SetValue(v, 0);
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));
404 //----------------------------------------------------------------------------
405 void vtkDijkstraImageData::FindNeighbors(vtkIdList *list,int id, vtkDataArray *scalars) {
407 // find i, j, k for that node
408 vtkImageData *input = this->GetInput();
410 int *dim = input->GetDimensions();
411 int numPts = dim[0] * dim[1] * dim[2];
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++) {
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);
437 float vtkDijkstraImageData::fuerzaAtraccion(int u, float w)
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];
444 this->GetInput()->GetPoint(this->GetSourceID(), coordsSource);
445 this->GetInput()->GetPoint(this->GetSinkID(), coordsSink);
448 //2. Identificar las coordenadas del punto al cual se le quiere sacar el potencial electrico
449 this->GetInput()->GetPoint(u, coordsPoint);
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]));
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]));
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]));
469 if (r2 < d && r1 > d){
471 } else if ( r2 > d && r1 < d){
473 } else if (r1 < d && r2 < d){
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)
489 float dist2 = scalars->GetTuple1(v); //+ fuerzaAtraccion(v, scalars->GetTuple1(v));
490 float dist = sqrt(dist2);
491 if(this->UseInverseDistance)
493 else if(this->UseInverseSquaredDistance)
495 else if(this->UseInverseExponentialDistance)
497 else if(this->UseSquaredDistance)
505 //----------------------------------------------------------------------------
506 void vtkDijkstraImageData::BuildShortestPath(int start,int end)
514 this->puntos = vtkPoints::New();
515 this->lineas = vtkCellArray::New();
519 //printf("======================camino minimo========================\n");
520 //printf("Punto Inicial = %i, Punto Final = %i\n",start, end);
522 while (p != start && p > 0)
526 this->GetInput()->GetPoint(p,punto);
527 this->puntos->InsertPoint(i,punto); //puntos y lineas usados como resultado (vtkPolyData)
529 next = this->Parent->GetValue(p);
531 this->GetInput()->GetPoint(next,punto);
532 this->puntos->InsertPoint(i+1,punto);
534 this->lineas->InsertNextCell(2);
535 this->lineas->InsertCellPoint(i);
536 this->lineas->InsertCellPoint(i+1);
538 this->ShortestPathIdList->InsertNextId(p);
543 this->ShortestPathIdList->InsertNextId(p);
545 //printf("===========================================================\n");
549 this->GetOutput()->SetPoints (this->puntos);
550 this->GetOutput()->SetLines(this->lineas);
551 this->GetOutput()->Modified();
556 vtkIdList* vtkDijkstraImageData::GetShortestPathIdList()
558 return this->ShortestPathIdList;
563 void vtkDijkstraImageData::printPointData(int pointID)
567 this->GetInput()->GetPoint(pointID, coords);
568 double scalar = this->GetInput()->GetPointData()->GetScalars()->GetTuple1(pointID);
570 printf("Punto ID: %i ",pointID);
571 printf("Coords[%.0f, %.0f, %.0f] ", coords[0], coords[1], coords[2]);
572 printf("Scalar: %.2f\n", scalar);
578 //----------------------------------------------------------------------------
581 void vtkDijkstraImageData::InitTraversePath(){
582 this->PathPointer = -1;
585 //----------------------------------------------------------------------------
586 int vtkDijkstraImageData::GetNumberOfPathNodes(){
587 return this->ShortestPathIdList->GetNumberOfIds();
590 //----------------------------------------------------------------------------
591 int vtkDijkstraImageData::GetNextPathNode(){
592 this->PathPointer = this->PathPointer + 1;
594 if(this->PathPointer < this->GetNumberOfPathNodes()) {
595 //printf("this->GetPathNode(this->PathPointer) %i",this->GetPathNode(this->PathPointer));
596 return this->ShortestPathIdList->GetId(this->PathPointer);
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) {
607 int kFactor = dim0 * dim1;
610 int numPoints = kFactor * dim2;
611 vtkIdList* Q = vtkIdList::New();
617 while(Q->GetNumberOfIds() != 0) {
619 int current = Q->GetId(pointer);
620 pointer = pointer + 1;
621 //printf("we are looking at id %i \n",current);
623 // check to see if we found something in the graph
624 if(scalars->GetTuple1(current) >0) {
625 //printf("before return");
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
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);
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);
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);
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);
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);
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);
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()
690 vtkImageData *inData = this->GetInput();
691 vtkPolyData *outData = this->GetOutput();
696 // Components turned into x, y and z
697 if (inData->GetNumberOfScalarComponents() > 3)
699 vtkErrorMacro("This filter can handle upto 3 components");
704 //printf("*************** vtkDijkstraImageData Execute **************\n\n\n");
706 //printf("*************** after init ****************\n\n\n");
708 int *dim = inData->GetDimensions();
709 vtkDataArray *scalars = inData->GetPointData()->GetScalars();
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) {
715 vtkFloatArray *copyScalars = vtkFloatArray::New();
716 copyScalars->DeepCopy(inData->GetPointData()->GetScalars());
718 this->SetSourceID(this->findClosestPointInGraph(copyScalars,this->GetSourceID(),dim[0],dim[1],dim[2]));
720 copyScalars->Delete();
721 printf("NEW source ID %i value is %f\n",this->GetSourceID(),scalars->GetTuple1(this->GetSourceID()));
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();
732 printf("NEW sink ID %i value is %f\n",this->GetSinkID(),scalars->GetTuple1(this->GetSinkID()));
736 this->RunDijkstra(scalars,this->GetSourceID(),this->GetSinkID());
738 this->BuildShortestPath(this->GetSourceID(),this->GetSinkID());
746 //----------------------------------------------------------------------------
747 void vtkDijkstraImageData::PrintSelf(ostream& os, vtkIndent indent)
749 Superclass::PrintSelf(os,indent);
751 os << indent << "Source ID: ( "
752 << this->GetSourceID() << " )\n";
754 os << indent << "Sink ID: ( "
755 << this->GetSinkID() << " )\n";