1 #include "vtkDijkstraImageData.h"
3 #include "vtkObjectFactory.h"
4 #include "vtkIntArray.h"
5 #include "vtkPointData.h"
6 #include "vtkPriorityQueue.h"
8 #include "vtkImageData.h"
9 #include "vtkFloatArray.h"
13 vtkCxxSetObjectMacro(vtkDijkstraImageData,BoundaryScalars,vtkDataArray);
15 //----------------------------------------------------------------------------
16 vtkDijkstraImageData* vtkDijkstraImageData::New()
22 // First try to create the object from the vtkObjectFactory
23 vtkObject* ret = vtkObjectFactory::CreateInstance("vtkDijkstraImageData");
26 return (vtkDijkstraImageData*)ret;
28 // If the factory was unable to create the object, then create it here.
29 return new vtkDijkstraImageData;
35 //----------------------------------------------------------------------------
36 // Constructor sets default values
37 vtkDijkstraImageData::vtkDijkstraImageData()
41 this->NumberOfInputPoints = 0;
42 this->PathPointer = -1;
43 this->StopWhenEndReached = 1;
45 this->ShortestPathIdList = NULL;
49 this->BoundaryScalars = NULL;
51 this->UseInverseDistance = 0;
52 this->UseInverseSquaredDistance = 0;
53 this->UseInverseExponentialDistance = 1;
54 this->UseSquaredDistance = 0;
62 //if ((this->logger = freopen("c:\\creatis\\maracas\\dijkstra.txt","w", stdout)) == NULL)
68 //------------------------------------------------------------------------------
69 vtkDijkstraImageData::~vtkDijkstraImageData()
73 if (this->ShortestPathIdList)
74 this->ShortestPathIdList->Delete();
76 this->Parent->Delete();
77 if(this->BoundaryScalars)
78 this->BoundaryScalars->Delete();
80 // this->Visited->Delete();
82 // this->PQ->Delete();
84 //this->Heap->Delete();
87 //if(this->BoundaryScalars)
88 //this->BoundaryScalars->Delete();
90 //if (this->CumulativeWeightFromSource)
91 //this->CumulativeWeightFromSource->Delete();
95 //------------------------------------------------------------------------------
96 unsigned long vtkDijkstraImageData::GetMTime()
98 unsigned long mTime=this->MTime.GetMTime();
104 void vtkDijkstraImageData::SetInput(vtkImageData* data){
105 this->vtkProcessObject::SetNthInput(0, data);
108 vtkImageData* vtkDijkstraImageData::GetInput(){
109 return (vtkImageData *)(this->Inputs[0]);
112 //------------------------------------------------------------------------------
113 void vtkDijkstraImageData::init(vtkImageData *inData)
116 if (this->ShortestPathIdList)
117 this->ShortestPathIdList->Delete();
120 this->Parent->Delete();
122 this->Visited->Delete();
126 this->ShortestPathIdList = vtkIdList::New();
127 this->Parent = vtkIntArray::New();
128 this->Visited = vtkIntArray::New();
129 this->PQ = vtkPriorityQueue::New();
135 int numPoints = inData->GetNumberOfPoints();
137 this->Parent->SetNumberOfComponents(1);
138 this->Parent->SetNumberOfTuples(numPoints);
139 this->Visited->SetNumberOfComponents(1);
140 this->Visited->SetNumberOfTuples(numPoints);
144 void vtkDijkstraImageData::DeleteGraph()
147 /*const int npoints = this->GetNumberOfInputPoints();
151 for (int i = 0; i < npoints; i++)
153 if(this->Neighbors[i])
154 this->Neighbors[i]->Delete();
156 delete [] this->Neighbors;
158 this->Neighbors = NULL;
162 //------------------------------------------------------------------------------
163 void vtkDijkstraImageData::CreateGraph(vtkImageData *inData) {
167 //delete old arrays in case we are re-executing this filter
168 int numPoints = inData->GetNumberOfPoints();
169 this->SetNumberOfInputPoints(numPoints);
172 int *dim = inData->GetDimensions();
173 vtkDataArray *scalars = inData->GetPointData()->GetScalars();
174 vtkIdList *graphNodes = vtkIdList::New();
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++) {
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
188 graphNodes->InsertNextId(id);
195 this->SetNumberOfGraphNodes(graphSize);
196 //printf("graph size %i \n ",graphSize);
199 PQ->Allocate(graphSize);
200 for(int i=0; i<graphSize;i++) {
201 PQ->Insert(VTK_LARGE_FLOAT,graphNodes->GetId(i));
204 graphNodes->Delete();
209 //------------------------------------------------------------------------------
210 void vtkDijkstraImageData::RunDijkstra(vtkDataArray *scalars,int startv, int endv)
215 printf("tipo de ponderacion de pesos : linear %i, squared %i, exponential %i\n",this->UseInverseDistance,this->UseInverseSquaredDistance,this->UseInverseExponentialDistance);
219 //-------------------------------------------------------------------------------
221 //-------------------------------------------------------------------------------
223 // 1. Se marcan todos los vertices del grafo como no visitados.
224 InitSingleSource(startv);
226 // 2. Se marca el punto de arranque como visitado
227 this->Visited->SetValue(startv, 1);
229 // 3. Se obtiene el tamanio inicial de los vertices no visitados
230 int initialSize = PQ->GetNumberOfItems();
231 int size = initialSize;
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.
237 //------------------------------------------------------------------------
238 // PARA TODOS LOS VERTICES EN PQ (SIN EVALUAR)
239 //------------------------------------------------------------------------
240 while ((PQ->GetNumberOfItems() > 0) && !stop)
243 this->UpdateProgress (0.5 + (float) (initialSize - size) / ((float) 2 * initialSize));
246 vtkFloatingPointType u_weight;
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);
256 u = PQ->Pop(0, u_weight);
258 //printf("Evaluando el nodo %i con prioridad %.2f\n",u,u_weight);
262 //-----------------------------------------
263 //2. Se marca "u" como ya visitado.
264 //-----------------------------------------
265 this->Visited->SetValue(u, 1);
269 //------------------------------------------------------------------------------------
270 //3. Si u es el nodo final el algoritmo se detiene. (despues de evaluar sus vecinos)
271 //------------------------------------------------------------------------------------
272 if (u == endv && StopWhenEndReached)
275 // Update all vertices v neighbors to u
276 // find the neighbors of u
280 //-----------------------------------------
281 //4. Encontrar los vecinos de u
282 //-----------------------------------------
283 vtkIdList *list = vtkIdList::New();
284 this->FindNeighbors(list,u,scalars);
286 //printf("%i tiene %i vecinos\n",u,list->GetNumberOfIds());
288 //-----------------------------------------
289 // 5. Para cada vecino ....
290 //-----------------------------------------
291 for (i = 0; i < list->GetNumberOfIds(); i++)
294 //---------------------------------------------------------------------
295 // 5.1 Se obtiene su ID (el indice es el orden en la lista de vecinos)
296 //---------------------------------------------------------------------
298 //printf("---> evaluando el vecino % i\n",v);
300 // s is the set of vertices with determined shortest path...do not
303 //-----------------------------------------------------------------
304 // 5.2 Si ya ha sido visitado se descarta...
305 //-----------------------------------------------------------------
306 if (this->Visited->GetValue(v) != 1)
308 //printf("---> ---> el vecino %i no ha sido visitado\n",v);
309 // Only relax edges where the end is not in s and edge
311 //---------------------------------------------------
312 // 5.2.1 Se obtiene el COSTO W
313 //---------------------------------------------------
314 float w = EdgeCost(scalars, u, v);
316 //---------------------------------------------------
318 //---------------------------------------------------
319 float v_weight = this->PQ->GetPriority(v);
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);
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);
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))
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);
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);
343 // 5.2.3.3 se coloca V conectado con U (esto sirve para reconstruir el camino minimo)
345 this->Parent->SetValue(v, u);
346 //printf("setting parent of %i to be %i",v,u);
350 //printf("---> el vecino %i ya se visito\n",v);
354 //-------------------------------------------------------
355 // 6. Liberar la memoria ocupada por la lista de vecinos
356 //-------------------------------------------------------
361 this->Visited->Delete();
364 //----------------------------------------------------------------------------
365 void vtkDijkstraImageData::InitSingleSource(int startv)
367 for (int v = 0; v < this->GetNumberOfInputPoints(); v++)
369 this->Parent->SetValue(v, -1);
370 this->Visited->SetValue(v, 0);
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));
379 //----------------------------------------------------------------------------
380 void vtkDijkstraImageData::FindNeighbors(vtkIdList *list,int id, vtkDataArray *scalars) {
382 // find i, j, k for that node
383 vtkImageData *input = this->GetInput();
385 int *dim = input->GetDimensions();
386 int numPts = dim[0] * dim[1] * dim[2];
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++) {
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);
412 float vtkDijkstraImageData::fuerzaAtraccion(int u, float w)
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];
419 this->GetInput()->GetPoint(this->GetSourceID(), coordsSource);
420 this->GetInput()->GetPoint(this->GetSinkID(), coordsSink);
423 //2. Identificar las coordenadas del punto al cual se le quiere sacar el potencial electrico
424 this->GetInput()->GetPoint(u, coordsPoint);
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]));
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]));
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]));
444 if (r2 < d && r1 > d){
446 } else if ( r2 > d && r1 < d){
448 } else if (r1 < d && r2 < d){
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)
464 float dist2 = scalars->GetTuple1(v); //+ fuerzaAtraccion(v, scalars->GetTuple1(v));
465 float dist = sqrt(dist2);
466 if(this->UseInverseDistance)
468 else if(this->UseInverseSquaredDistance)
470 else if(this->UseInverseExponentialDistance)
472 else if(this->UseSquaredDistance)
480 //----------------------------------------------------------------------------
481 void vtkDijkstraImageData::BuildShortestPath(int start,int end)
489 this->puntos = vtkPoints::New();
490 this->lineas = vtkCellArray::New();
494 //printf("======================camino minimo========================\n");
495 //printf("Punto Inicial = %i, Punto Final = %i\n",start, end);
497 while (p != start && p > 0)
501 this->GetInput()->GetPoint(p,punto);
502 this->puntos->InsertPoint(i,punto); //puntos y lineas usados como resultado (vtkPolyData)
504 next = this->Parent->GetValue(p);
506 this->GetInput()->GetPoint(next,punto);
507 this->puntos->InsertPoint(i+1,punto);
509 this->lineas->InsertNextCell(2);
510 this->lineas->InsertCellPoint(i);
511 this->lineas->InsertCellPoint(i+1);
513 this->ShortestPathIdList->InsertNextId(p);
518 this->ShortestPathIdList->InsertNextId(p);
520 //printf("===========================================================\n");
524 this->GetOutput()->SetPoints (this->puntos);
525 this->GetOutput()->SetLines(this->lineas);
526 this->GetOutput()->Modified();
531 vtkIdList* vtkDijkstraImageData::GetShortestPathIdList()
533 return this->ShortestPathIdList;
538 void vtkDijkstraImageData::printPointData(int pointID)
542 this->GetInput()->GetPoint(pointID, coords);
543 double scalar = this->GetInput()->GetPointData()->GetScalars()->GetTuple1(pointID);
545 printf("Punto ID: %i ",pointID);
546 printf("Coords[%.0f, %.0f, %.0f] ", coords[0], coords[1], coords[2]);
547 printf("Scalar: %.2f\n", scalar);
553 //----------------------------------------------------------------------------
556 void vtkDijkstraImageData::InitTraversePath(){
557 this->PathPointer = -1;
560 //----------------------------------------------------------------------------
561 int vtkDijkstraImageData::GetNumberOfPathNodes(){
562 return this->ShortestPathIdList->GetNumberOfIds();
565 //----------------------------------------------------------------------------
566 int vtkDijkstraImageData::GetNextPathNode(){
567 this->PathPointer = this->PathPointer + 1;
569 if(this->PathPointer < this->GetNumberOfPathNodes()) {
570 //printf("this->GetPathNode(this->PathPointer) %i",this->GetPathNode(this->PathPointer));
571 return this->ShortestPathIdList->GetId(this->PathPointer);
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) {
582 int kFactor = dim0 * dim1;
585 int numPoints = kFactor * dim2;
586 vtkIdList* Q = vtkIdList::New();
592 while(Q->GetNumberOfIds() != 0) {
594 int current = Q->GetId(pointer);
595 pointer = pointer + 1;
596 //printf("we are looking at id %i \n",current);
598 // check to see if we found something in the graph
599 if(scalars->GetTuple1(current) >0) {
600 //printf("before return");
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
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);
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);
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);
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);
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);
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);
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()
665 vtkImageData *inData = this->GetInput();
666 vtkPolyData *outData = this->GetOutput();
671 // Components turned into x, y and z
672 if (inData->GetNumberOfScalarComponents() > 3)
674 vtkErrorMacro("This filter can handle upto 3 components");
679 //printf("*************** vtkDijkstraImageData Execute **************\n\n\n");
681 //printf("*************** after init ****************\n\n\n");
683 int *dim = inData->GetDimensions();
684 vtkDataArray *scalars = inData->GetPointData()->GetScalars();
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) {
690 vtkFloatArray *copyScalars = vtkFloatArray::New();
691 copyScalars->DeepCopy(inData->GetPointData()->GetScalars());
693 this->SetSourceID(this->findClosestPointInGraph(copyScalars,this->GetSourceID(),dim[0],dim[1],dim[2]));
695 copyScalars->Delete();
696 printf("NEW source ID %i value is %f\n",this->GetSourceID(),scalars->GetTuple1(this->GetSourceID()));
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();
707 printf("NEW sink ID %i value is %f\n",this->GetSinkID(),scalars->GetTuple1(this->GetSinkID()));
711 this->RunDijkstra(scalars,this->GetSourceID(),this->GetSinkID());
713 this->BuildShortestPath(this->GetSourceID(),this->GetSinkID());
721 //----------------------------------------------------------------------------
722 void vtkDijkstraImageData::PrintSelf(ostream& os, vtkIndent indent)
724 Superclass::PrintSelf(os,indent);
726 os << indent << "Source ID: ( "
727 << this->GetSourceID() << " )\n";
729 os << indent << "Sink ID: ( "
730 << this->GetSinkID() << " )\n";