lib/Airways/AirwaysLib/airwaysTree.cxx

   1 /*
   2  * airwaysTree.txx
   3  *
   4  *  Created on: May 3, 2014
   5  *      Author: caceres@creatis.insa-lyon.fr
   6  *      Modifications by: Alfredo Morales Pinzón
   7  */
   8
   9 #ifndef _AIRWAYS_TREE_TXX_
  10 #define _AIRWAYS_TREE_TXX_
  11
  12 #include "airwaysTree.h"
  13
  14 namespace airways
  15 {
  16
  17 AirwaysTree::AirwaysTree() :  m_root(NULL)
  18 {
  19
  20 }
  21
  22 //Ceron
  23 AirwaysTree::AirwaysTree(TInputImage::Pointer img, TInputImage::Pointer img_skele, Node* root, bool shouldUpdate){
  24     this->m_root = root;
  25     this->m_root->SetLevel(0);
  26     this->m_img = img;
  27     this->m_skeleton = img_skele;
  28     //TODO: Should this calculate levels and populate edges and nodes vectors?
  29     FillTree(root);
  30     if(shouldUpdate) this->UpdateEdges();
  31 }
  32
  33 void AirwaysTree::FillTree(Node* node){
  34     //std::cout << "Calls fill tree for " << node->GetCoords() << std::endl;
  35     this->m_nodes.push_back(node);
  36     node->SetId(this->m_nodes.size());
  37     if(node->GetEdge() != NULL){
  38         this->m_edges.push_back(node->GetEdge());
  39         node->GetEdge()->UpdateEdgeInfo();
  40     }
  41     int level = node->GetLevel();
  42     //std::cout << "Node is at level " << level << std::endl;
  43     vec_nodes::const_iterator it = node->GetChildren().begin();
  44         for (; it != node->GetChildren().end(); ++it){
  45         (*it)->SetLevel(level+1);
  46         this->FillTree(*it);
  47     }
  48 }
  49
  50 AirwaysTree::AirwaysTree(TInputImage::Pointer img, TInputImage::Pointer img_sk, Node* root)
  51 {
  52         this->m_root = root;
  53         this->m_root->SetLevel(0);
  54         this->m_skeleton = img_sk;
  55         this->m_img = img;
  56         this->UpdateEdges();
  57 }
  58
  59 AirwaysTree::AirwaysTree(TInputImage::Pointer image_airwaySegmentation, TInputImage::Pointer image_skeletonDM, EdgeVector& vector_edges, Edge* edge_trachea)
  60 {
  61         std::cout << "AirwaysTree::AirwaysTree(TInputImage::Pointer img, TInputImage::Pointer img_sk, EdgeVector& tEdgeVector, Edge* trachea)" << std::endl;
  62         std::cout << "Edges to add:" << vector_edges.size() << std::endl;
  63
  64         if(edge_trachea == NULL)
  65                 std::cout << "The trachea edge is NULL" << std::endl;
  66
  67         // Set the attribtes
  68         this->m_skeleton = image_skeletonDM;
  69         this->m_img = image_airwaySegmentation;
  70         this->m_root = edge_trachea->m_source;
  71         this->m_root->SetLevel(0);
  72
  73         this->m_nodes.push_back(edge_trachea->m_source);        // Add the first node, the root to the nodes of the tree
  74         edge_trachea->m_source->SetId(this->m_nodes.size());// Set the root id
  75
  76         // Insert the trachea edge
  77         this->Insert(edge_trachea);
  78
  79         // Insert all edges
  80         unsigned int iniSize = vector_edges.size();
  81
  82         while (!vector_edges.empty())
  83         {
  84                 bool oneEdgeAdded = false;
  85                 EdgeVector::iterator it = vector_edges.begin();
  86                 while (it != vector_edges.end())
  87                 {
  88                         if (this->Insert(*it))
  89                         {
  90                                 std::cout << "Edges added: [" << (*it)->GetSource()->GetCoords() << "]->[" << (*it)->GetTarget()->GetCoords() << "], missing: " << vector_edges.size()-1 << std::endl;
  91                                 oneEdgeAdded = true;
  92                                 it = vector_edges.erase(it);
  93                         }
  94                         else
  95                                 ++it;
  96                 } //while
  97
  98                 // If there are edges to be added but they could not be added
  99                 if(!oneEdgeAdded)
 100                 {
 101                         std::cout << "Edges not added:" << vector_edges.size() << std::endl;
 102                         for (EdgeVector::iterator it = vector_edges.begin(); it != vector_edges.end(); ++it)
 103                                 std::cout << "Edge missing: [" << (*it)->GetSource()->GetCoords() << "]->[" << (*it)->GetTarget()->GetCoords() << "]" << std::endl;
 104                         break;
 105                 }
 106                 // AMP: I did not find any explanation for this "exception"
 107                 /*count++;
 108                 if (count == iniSize)
 109                 {
 110                         std::cout << "handled exception: arwaysTree.txx, AirwaysTree" << std::endl;
 111                         break;
 112                 } //if
 113                  */
 114         } //while
 115         this->UpdateEdges();
 116
 117         std::cout << "AirwaysTree::AirwaysTree(TInputImage::Pointer img, TInputImage::Pointer img_sk, EdgeVector& tEdgeVector, Edge* trachea) ...  OK" << std::endl;
 118 }
 119
 120 AirwaysTree::~AirwaysTree()
 121 {
 122
 123 }
 124
 125 Node* AirwaysTree::GetRoot() const
 126 {
 127         return this->m_root;
 128 }
 129
 130 std::string AirwaysTree::GetImageName() const
 131 {
 132         return (this->m_img_name);
 133 }
 134
 135 std::string AirwaysTree::GetResultPath() const
 136 {
 137         return (this->m_result_path);
 138 }
 139
 140 TInputImage::Pointer AirwaysTree::GetSegmentedImage() const
 141 {
 142         return this->m_img;
 143 }
 144
 145 TInputImage::Pointer AirwaysTree::GetSkeletonImage() const
 146 {
 147         return this->m_skeleton;
 148 }
 149
 150 const Node* AirwaysTree::GetNode(const Vec3 nodeCoords) const
 151 {
 152         return this->GetNode(nodeCoords);
 153 }
 154
 155 Node* AirwaysTree::GetNodeById(unsigned int nId)
 156 {
 157         return this->m_root->GetNodeById(nId);
 158 }
 159
 160 int AirwaysTree::GetDepthById(unsigned int nId)
 161 {
 162         return this->m_root->GetDepthById(nId);
 163 }
 164
 165 unsigned int AirwaysTree::CountMarkedNodes(Node* current) const
 166 {
 167         if (current == NULL)
 168                 current = this->m_root;
 169         Node* aux = NULL;
 170         bool found = current->FindMarked(aux);
 171         if (!found && !current->IsMarked())
 172                 return 0;
 173         if (aux != NULL && found)
 174                 current = aux;
 175         unsigned int count = 1;
 176         vec_nodes::const_iterator it = current->GetChildren().begin();
 177         for (; it != current->GetChildren().end(); ++it)
 178                 count += this->CountMarkedNodes(*it);
 179         return count;
 180 }
 181
 182 unsigned int AirwaysTree::GetWeight( ) const
 183 {
 184         if (this->m_root == NULL)
 185                 return 0;
 186
 187         return this->m_root->GetWeight();
 188 }
 189
 190 unsigned int AirwaysTree::GetHeight(Node* current) const
 191 {
 192         if (current == NULL)
 193                 current = this->m_root;
 194         if (current->IsLeaf())
 195                 return 1;
 196         unsigned int max = 0;
 197         vec_nodes::const_iterator it = current->GetChildren().begin();
 198         for (; it != current->GetChildren().end(); ++it)
 199         {
 200                 unsigned int cMax = this->GetHeight(*it);
 201                 if (cMax > max)
 202                         max = cMax;
 203         } //rof
 204         if (current == this->m_root)
 205                 return max;
 206         return ++max;
 207 }
 208
 209 unsigned int AirwaysTree::GetLevels(Node* current, unsigned int cLevel) const
 210 {
 211         if (current == NULL)
 212                 current = this->m_root;
 213         if (current->IsLeaf())
 214                 return current->GetLevel();
 215         unsigned int maxLevel = cLevel;
 216         vec_nodes::const_iterator it = current->GetChildren().begin();
 217         for (; it != current->GetChildren().end(); ++it)
 218         {
 219                 unsigned int level = this->GetLevels(*it, maxLevel);
 220                 if (level > maxLevel)
 221                         maxLevel = level;
 222         } //rof
 223         return maxLevel;
 224 }
 225
 226 unsigned int AirwaysTree::GetNumberOfLeafs( ) const
 227 {
 228         if(!this->m_root)
 229                 return 0;
 230
 231         return this->m_root->GetNumberOfLeafs( );
 232 }
 233
 234 const EdgeVector& AirwaysTree::GetEdges() const
 235 {
 236         return this->m_edges;
 237 }
 238
 239 const vec_nodes& AirwaysTree::GetNodes() const
 240 {
 241         return this->m_nodes;
 242 }
 243
 244 const Node* AirwaysTree::GetClosestBrachNodeToIndex(float point_x, float point_y, float point_z) const
 245 {
 246         float distance = 0;
 247         Node* closestNode = this->m_root->GetClosestBranchNodeToPoint(point_x, point_y, point_z, distance);
 248         return closestNode;
 249 }
 250
 251 void AirwaysTree::GetBrancheNodesInfluencingPoint(float point_x, float point_y, float point_z, vec_nodes& vec_nodes_influence)
 252 {
 253         if(this->m_root)
 254         {
 255                 std::cout << "AirwaysTree::GetBrancheNodesInfluencingPoint" << std::endl;
 256                 this->m_root->GetBrancheNodesInfluencingPoint(point_x, point_y, point_z, vec_nodes_influence);
 257         }
 258 }
 259
 260 AirwaysTree& AirwaysTree::GetTreeFromMarkedNodes(Node* currentNode, AirwaysTree* tree)
 261 {
 262         if (tree == NULL)
 263         {
 264                 if (!this->m_root->IsMarked())
 265                         return *(new AirwaysTree());
 266                 currentNode = this->m_root;
 267                 tree = new AirwaysTree();
 268                 Node* current = new Node(currentNode->GetCoords());
 269                 tree->SetRoot(current);
 270                 vec_nodes::const_iterator it = currentNode->GetChildren().begin();
 271                 for (; it != currentNode->GetChildren().end(); ++it)
 272                         this->GetTreeFromMarkedNodes(*it, tree);
 273                 this->ImageReconstruction(tree->m_img);
 274                 //this->ImageReconstruction(tree->m_skeleton);
 275                 return *tree;
 276         } //fi
 277         if (currentNode->IsMarked())
 278         {
 279                 Edge* edge = new Edge();
 280                 *edge = *currentNode->GetEdge();
 281                 Node* target = new Node(currentNode->GetCoords());
 282                 edge->m_target = target;
 283                 tree->Insert(edge);
 284                 vec_nodes::const_iterator it = currentNode->GetChildren().begin();
 285                 for (; it != currentNode->GetChildren().end(); ++it)
 286                         this->GetTreeFromMarkedNodes(*it, tree);
 287         } //fi
 288         return *(new AirwaysTree());
 289 }
 290
 291 AirwaysTree& AirwaysTree::GetSubTreeFromNode(Node* node, AirwaysTree* tree)
 292 {
 293         if (tree == NULL)
 294         {
 295                 tree = new AirwaysTree();
 296                 Node* current = new Node(node->GetCoords());
 297                 tree->SetRoot(current);
 298                 vec_nodes::const_iterator it = node->GetChildren().begin();
 299                 for (; it != node->GetChildren().end(); ++it)
 300                         this->GetSubTreeFromNode(*it, tree);
 301                 this->ImageReconstructionFromNode(tree->m_img, node);
 302                 //this->ImageReconstructionFromNode(tree->m_skeleton, node, true);
 303                 return *tree;
 304         } //fi
 305         Edge* edge = new Edge();
 306         *edge = *node->GetEdge();
 307         Node* target = new Node(node->GetCoords());
 308         edge->m_target = target;
 309         tree->Insert(edge);
 310         vec_nodes::const_iterator it = node->GetChildren().begin();
 311         for (; it != node->GetChildren().end(); ++it)
 312                 this->GetSubTreeFromNode(*it, tree);
 313
 314         return *(new AirwaysTree());
 315 }
 316
 317 AirwaysTree& AirwaysTree::GetSubTreeByLevels(const int& level, Node* node,
 318                 AirwaysTree* tree)
 319 {
 320         if (tree == NULL)
 321         {
 322                 node = this->m_root;
 323                 tree = new AirwaysTree();
 324                 Node* current = new Node(node->GetCoords());
 325                 tree->SetRoot(current);
 326                 if (node->GetLevel() <= level)
 327                 {
 328                         node->Mark();
 329                         vec_nodes::const_iterator it = node->GetChildren().begin();
 330                         for (; it != node->GetChildren().end(); ++it)
 331                                 this->GetSubTreeByLevels(level, *it, tree);
 332                 } //fi
 333                 if (node->GetLevel() == level)
 334                 {
 335                         DuplicatorType::Pointer duplicator = DuplicatorType::New();
 336                         duplicator->SetInputImage(this->m_img);
 337                         duplicator->Update();
 338                         tree->m_img = duplicator->GetOutput();
 339                         duplicator->SetInputImage(this->m_skeleton);
 340                         duplicator->Update();
 341                         tree->m_skeleton = duplicator->GetOutput();
 342                 } //if
 343                 else
 344                 {
 345                         this->ImageReconstruction(tree->m_img);
 346                         //this->ImageReconstruction(tree->m_skeleton, true);
 347                 } //else
 348                 this->UnMarkAll();
 349                 return *tree;
 350         } //fi
 351         if (node->GetLevel() <= level)
 352         {
 353                 node->Mark();
 354                 Edge* edge = new Edge();
 355                 *edge = *node->GetEdge();
 356                 Node* target = new Node(node->GetCoords());
 357                 edge->m_target = target;
 358                 tree->Insert(edge);
 359                 vec_nodes::const_iterator it = node->GetChildren().begin();
 360                 for (; it != node->GetChildren().end(); ++it)
 361                         this->GetSubTreeByLevels(level, *it, tree);
 362                 if (node->GetChildren().size() == 1)
 363                 {
 364                         Node* child = *(node->GetChildren().begin());
 365                         node->MergeNodes(child);
 366                         const Edge* edge = (child)->GetEdge();
 367                         //removing from nodeVector and edgeVector
 368                         EdgeVector::iterator it2 = this->m_edges.begin();
 369                         for (; it2 != this->m_edges.end(); ++it2)
 370                                 if (edge == *it2)
 371                                 {
 372                                         this->m_edges.erase(it2);
 373                                         break;
 374                                 } //if
 375                         vec_nodes::iterator it3 = this->m_nodes.begin();
 376                         for (; it3 != this->m_nodes.end(); ++it3)
 377                                 if (child == *it3)
 378                                 {
 379                                         this->m_nodes.erase(it3);
 380                                         break;
 381                                 } //if
 382                 } //if
 383         } //if
 384         return *(new AirwaysTree());
 385 }
 386
 387 AirwaysTree& AirwaysTree::GetSubTreeByLength(const double& length, Node* node,
 388                 AirwaysTree* tree, double cLenght)
 389 {
 390         if (tree == NULL)
 391         {
 392                 node = this->m_root;
 393                 tree = new AirwaysTree();
 394                 Node* current = new Node(node->GetCoords());
 395                 tree->SetRoot(current);
 396                 node->Mark();
 397                 vec_nodes::const_iterator it = node->GetChildren().begin();
 398                 for (; it != node->GetChildren().end(); ++it)
 399                         this->GetSubTreeByLength(length, *it, tree);
 400                 this->ImageReconstruction(tree->m_img);
 401                 //this->ImageReconstruction(tree->m_skeleton, true);
 402                 this->UnMarkAll();
 403                 return *tree;
 404         } //fi
 405         cLenght += node->GetEdge()->m_length;
 406         if (cLenght <= length)
 407         {
 408                 node->Mark();
 409                 Edge* edge = new Edge();
 410                 *edge = *node->GetEdge();
 411                 Node* target = new Node(node->GetCoords());
 412                 edge->m_target = target;
 413                 tree->Insert(edge);
 414                 vec_nodes::const_iterator it = node->GetChildren().begin();
 415                 for (; it != node->GetChildren().end(); ++it)
 416                         this->GetSubTreeByLength(length, *it, tree, cLenght);
 417                 if (node->GetChildren().size() == 1)
 418                 {
 419                         Node* child = *(node->GetChildren().begin());
 420                         node->MergeNodes(child);
 421                         const Edge* edge = (child)->GetEdge();
 422                         //removing from nodeVector and edgeVector
 423                         EdgeVector::iterator it2 = this->m_edges.begin();
 424                         for (; it2 != this->m_edges.end(); ++it2)
 425                                 if (edge == *it2)
 426                                 {
 427                                         this->m_edges.erase(it2);
 428                                         break;
 429                                 } //if
 430                         vec_nodes::iterator it3 = this->m_nodes.begin();
 431                         for (; it3 != this->m_nodes.end(); ++it3)
 432                                 if (child == *it3)
 433                                 {
 434                                         this->m_nodes.erase(it3);
 435                                         break;
 436                                 } //if
 437                 } //if
 438         } //if
 439         return *(new AirwaysTree());
 440 }
 441
 442 AirwaysTree& AirwaysTree::GetSubTreeByRadius(const double& radius, Node* node,
 443                 AirwaysTree* tree, double cRadius)
 444 {
 445         if (tree == NULL)
 446         {
 447                 node = this->m_root;
 448                 tree = new AirwaysTree();
 449                 Node* current = new Node(node->GetCoords());
 450                 tree->SetRoot(current);
 451                 node->Mark();
 452                 vec_nodes::const_iterator it = node->GetChildren().begin();
 453                 for (; it != node->GetChildren().end(); ++it)
 454                         this->GetSubTreeByRadius(radius, *it, tree);
 455                 this->ImageReconstruction(tree->m_img);
 456                 //this->ImageReconstruction(tree->m_skeleton, true);
 457                 this->UnMarkAll();
 458                 return *tree;
 459         } //fi
 460         cRadius += node->GetEdge()->m_aRadius;
 461         if (cRadius >= radius)
 462         {
 463                 node->Mark();
 464                 Edge* edge = new Edge();
 465                 *edge = *node->GetEdge();
 466                 Node* target = new Node(node->GetCoords());
 467                 edge->m_target = target;
 468                 tree->Insert(edge);
 469                 vec_nodes::const_iterator it = node->GetChildren().begin();
 470                 for (; it != node->GetChildren().end(); ++it)
 471                         this->GetSubTreeByRadius(radius, *it, tree, cRadius);
 472                 if (node->GetChildren().size() == 1)
 473                 {
 474                         Node* child = *(node->GetChildren().begin());
 475                         node->MergeNodes(child);
 476                         const Edge* edge = (child)->GetEdge();
 477                         //removing from nodeVector and edgeVector
 478                         EdgeVector::iterator it2 = this->m_edges.begin();
 479                         for (; it2 != this->m_edges.end(); ++it2)
 480                                 if (edge == *it2)
 481                                 {
 482                                         this->m_edges.erase(it2);
 483                                         break;
 484                                 } //if
 485                         vec_nodes::iterator it3 = this->m_nodes.begin();
 486                         for (; it3 != this->m_nodes.end(); ++it3)
 487                                 if (child == *it3)
 488                                 {
 489                                         this->m_nodes.erase(it3);
 490                                         break;
 491                                 } //if
 492                 } //if
 493         } //if
 494         return *(new AirwaysTree());
 495 }
 496
 497 void AirwaysTree::MarkLevels(const int& level, Node* currentNode)
 498 {
 499         if (currentNode == NULL)
 500                 currentNode = this->m_root;
 501         if (currentNode->GetLevel() <= level)
 502         {
 503                 currentNode->Mark();
 504                 vec_nodes::const_iterator it = currentNode->GetChildren().begin();
 505                 for (; it != currentNode->GetChildren().end(); it++)
 506                         this->MarkLevels(level, *it);
 507         } //if
 508 }
 509
 510 bool AirwaysTree::IsEmpty() const
 511 {
 512         return this->m_root == NULL;
 513 }
 514
 515 bool AirwaysTree::IsNodeInsidePath(unsigned int idNode_starting, unsigned int  idNode_ending, unsigned int  idNode_searched)
 516 {
 517         // Check if the starting point exist
 518         Node* node_starting = this->m_root->GetNodeById(idNode_starting);
 519         if(!node_starting)
 520                 return false;
 521
 522         // Check if the searched node exist in the subtree of the starting node
 523         Node* node_searched = node_starting->GetNodeById(idNode_searched);
 524         if(!node_searched)
 525                 return false;
 526
 527         // Check if the ending node exist in the subtree of the searched node
 528         if(!node_searched->HasNodeId(idNode_ending))
 529                 return false;
 530
 531         // If all conditions were true then the node is in the path
 532         return true;
 533 }
 534
 535 void AirwaysTree::SetRoot(Node* root)
 536 {
 537         this->m_root = root;
 538         this->m_root->SetLevel(0);
 539         this->m_nodes.push_back(this->m_root);
 540 }
 541
 542 void AirwaysTree::SetImageName(const std::string& img_name)
 543 {
 544         this->m_img_name = img_name;
 545 }
 546
 547 void AirwaysTree::SetResultPath(const std::string& folderPath)
 548 {
 549         this->m_result_path = folderPath;
 550 }
 551
 552 bool AirwaysTree::Insert(Edge* edge)
 553 {
 554         if (this->m_root == NULL)
 555                 return false;
 556         // Search if the source or target coordinates exist in the tree
 557         Node* source = this->GetNode(edge->m_source->GetCoords());
 558         Node* target = this->GetNode(edge->m_target->GetCoords());
 559
 560         // If the source and the target coordinates do not exist then the
 561         // edge can not be added, also if both exist in the tree.
 562         if ( ( source == NULL && target == NULL ) || ( source != NULL && target != NULL) )
 563                 return false;
 564         else if (source == NULL && target != NULL) // If target was found
 565         {
 566                 // Change the direction of the edge
 567                 edge->m_target = edge->m_source;
 568                 edge->m_source = target;
 569         }
 570         else
 571                 edge->m_source = source;
 572
 573         // Set the target and source nodes for the Edge
 574         edge->m_target->SetEdge(edge);
 575         edge->m_target->SetLevel(edge->m_source->GetLevel() + 1);
 576         edge->m_source->AddChild(edge->m_target);
 577
 578
 579         // Add the new child node
 580         this->m_nodes.push_back(edge->m_target);
 581         this->m_edges.push_back(edge);
 582
 583         // Set the id node
 584         edge->m_target->SetId(this->m_nodes.size());
 585         return true;
 586 }
 587
 588 void AirwaysTree::UnMarkAll(Node* currentNode)
 589 {
 590         if (currentNode == NULL)
 591                 currentNode = this->m_root;
 592         currentNode->UnMark();
 593
 594         vec_nodes::const_iterator it = currentNode->GetChildren().begin();
 595         for (; it != currentNode->GetChildren().end(); ++it)
 596                 this->UnMarkAll(*it);
 597 }
 598
 599 void AirwaysTree::UpdateLevels(unsigned int levels)
 600 {
 601         this->m_root->UpdateLevels(levels);
 602 }
 603
 604 void AirwaysTree::CompareTreesOrkiszMorales(AirwaysTree& treeB, unsigned int Q, unsigned int F, vec_nodes& commonA, vec_nodes& commonB, vec_nodes& nonCommonA, vec_nodes& nonCommonB, std::map< unsigned int, std::vector<Node*> >& map_A_to_B, std::map< unsigned int, std::vector<Node*> >& map_B_to_A, std::vector< std::pair< std::pair<Node*, Node*>, std::pair<Node*, Node*> > >& vector_pair_edges_A_to_B)
 605 {
 606         // Algorithm
 607         // The algorithm is recursive it is implemented the node class. The idea is to overlap the initial nodes of the subtrees
 608         // and find the best branch correspondence for the first level of the trees.
 609         // Hypothesis: It is supposed that the roots of the trees are common, so the algorithm does not run in the root nodes
 610         // 1. Translate one of the subtree so that the initial nodes overlap.
 611         // 2. While one of the trees has non-marked child
 612         //   2.1 Compare each child to all the branches is the other tree.
 613         //   2.2 Select the pair, child and branch in the other tree, that has the lowest distance function.
 614         //   2.3 Add the pair and the distance to the final result.
 615         //   2.4 Mark the nodes
 616         //   2.5 Find and mark the omitted nodes in the similar branch. If the selected branch has more that 2 nodes it means that intemediary nodes do not have correpondece.
 617         //   2.6 Run the process (First step) using the found pair of nodes.
 618
 619         if( !this->IsEmpty() && !treeB.IsEmpty() )
 620         {
 621
 622                 // ----------------------------
 623                 // ---- Pre-Matching steps ----
 624
 625                 // Get the root nodes
 626                 Node* root_A = this->m_root;
 627                 Node* root_B = treeB.m_root;
 628
 629                 // Mark the roots
 630                 root_A->Mark();
 631                 root_B->Mark();
 632
 633                 // Add the root match
 634                 map_A_to_B[root_A->GetId()].push_back(root_B);
 635                 map_B_to_A[root_B->GetId()].push_back(root_A);
 636
 637                 // Add the root nodes match
 638                 pair_nodes pair_edgeA(root_A, root_A);
 639                 pair_nodes pair_edgeB(root_B, root_B);
 640                 std::pair< pair_nodes, pair_nodes > pairRoorNodes_edge(pair_edgeA, pair_edgeB);
 641                 vector_pair_edges_A_to_B.push_back(pairRoorNodes_edge);
 642
 643                 // ----------------------------
 644                 // ----- Run the comparison----
 645
 646                 root_A->CompareSubTreesOrkiszMorales(root_B, Q, F, commonA, commonB, nonCommonA, nonCommonB, map_A_to_B, map_B_to_A, vector_pair_edges_A_to_B);
 647
 648         }
 649         else
 650                 std::cout << "Trees are empties." << std::endl;
 651 }
 652
 653 void AirwaysTree::MarkPathFromNodeToNode(Node* node_begin, Node* node_end)
 654 {
 655         if(this->m_root)
 656                 this->m_root->MarkPathFromNodeToNode(node_begin, node_end);
 657 }
 658
 659 void AirwaysTree::CompareTrees(AirwaysTree& treeB)
 660 {
 661         // Compare all the "branches" (paths) in treeA with all branches in treeB
 662         Vector_Pair_PairNodes_Score vectorPairNodesScore = this->CompareAllBranches(treeB);
 663
 664         std::cout << "Vector pairs size:" <<  vectorPairNodesScore.size() << std::endl;
 665
 666         // Algorithm to find the comment branches (paths)
 667         // Until one of the trees is empty
 668         // 1. Find the pair, with at least one leaf in the pair, with the lowest score
 669         // 2. Compare the score of the father node of the leaf, or one of the leafs, with the correspondent pair leaf
 670         //   If the score is lower that the score of the actual pair
 671         //      - The compared leaf is not common. Mark the compared leaf, delete the compared leave and save the leaf as not common.
 672         //   If the score is greater
 673         //      - The leafs are common then mark both leafs as commons. Mark both leaves, delete them and save them as common leafs.
 674
 675         // Variables
 676         //AirwaysTree* tree_copyA = this->GetCopy();
 677         //AirwaysTree* tree_copyB = treeB.GetCopy();
 678
 679         //NodeVector nonCommonNodes_A;
 680         //NodeVector commonNodes_A;
 681         //NodeVector nonCommonNodes_B;
 682         //NodeVector commonNodes_B;
 683
 684         std::map<unsigned int, Node*> nonCommonNodes_A;
 685         std::map<unsigned int, Node*> commonNodes_A;
 686         std::map<unsigned int, Node*> nonCommonNodes_B;
 687         std::map<unsigned int, Node*> commonNodes_B;
 688
 689         // Until one of the trees is not empty
 690         while( this->GetWeight() > 1 && treeB.GetWeight() > 1 )
 691         {
 692                 // Print the weight of the trees
 693                 std::cout << "Weight [A,B]: [" << this->GetWeight() << "," << treeB.GetWeight() << "]" << std::endl;
 694                 std::cout << "Leafs [A,B]: [" << this->GetNumberOfLeafs() << "," << treeB.GetNumberOfLeafs() << "]" << std::endl;
 695                 std::cout << "Vector pairs size:" <<  vectorPairNodesScore.size() << std::endl;
 696
 697
 698                 // 1. Find the pair, with at least one leaf in the pair, with the lowest score
 699                 Pair_PairNodes_Score pair_nonMarked = this->GetBestLeaf(vectorPairNodesScore);
 700                 Node* node_actualA = pair_nonMarked.first.first;
 701                 Node* node_actualB = pair_nonMarked.first.second;
 702                 double score_actual = pair_nonMarked.second;
 703
 704                 // 2. Take the leaf (L) and his father (F), the other node is called (O). Compare the score of the father node (F) with the correspondent pair leaf (O).
 705                 if(node_actualA->IsLeaf())
 706                 {
 707                         //const Edge* edgeActual = node_actualA->GetEdge();
 708                         //Node* targetActual = edgeActual->GetTarget();
 709                         //const unsigned int id_father = targetActual->GetId();
 710                         //double score_fatherA_to_B = this->GetPairNodes_Score(vectorPairNodesScore, 0, node_actualB->GetId());
 711                         double score_fatherA_to_B = this->GetPairNodes_Score(vectorPairNodesScore, node_actualA->GetEdge()->GetTarget()->GetId(), node_actualB->GetId());
 712                         //double score_fatherA_to_B =  pair_fatherA_to_B.second;
 713                         // If the score is lower than the score of the actual pair
 714                         if(score_fatherA_to_B < score_actual)
 715                         {
 716                                 //- The compared leaf is not common. Mark the compared leaf, delete the compared leave and save the leaf as not common.
 717                                 node_actualA->Mark();
 718                                 std::map<unsigned int, Node*>::iterator it;
 719                                 it=nonCommonNodes_A.find(node_actualA->GetId());
 720                                 if(it==nonCommonNodes_A.end())
 721                                         nonCommonNodes_A[node_actualA->GetId()] = node_actualA;
 722                                 //nonCommonNodes_A.push_back(node_actualA);
 723                         }
 724                         // If the score is greater
 725                         else
 726                         {
 727                                 // - The leafs are common then mark both leafs as commons. Mark both leaves, delete them and save them as common leafs.
 728                                 node_actualA->Mark();
 729                                 node_actualB->Mark();
 730
 731                                 std::map<unsigned int, Node*>::iterator it;
 732                                 it=commonNodes_B.find(node_actualB->GetId());
 733                                 if(it==commonNodes_B.end())
 734                                         commonNodes_B[node_actualB->GetId()] = node_actualB;
 735
 736                                 it=commonNodes_A.find(node_actualA->GetId());
 737                                 if(it==commonNodes_A.end())
 738                                 {
 739                                         std::cout << "CommonNode A:" << node_actualA->GetId() << std::endl;
 740                                         commonNodes_A[node_actualA->GetId()] = node_actualA;
 741                                 }
 742
 743                                 //commonNodes_B.push_back(node_actualB);
 744                                 //commonNodes_A.push_back(node_actualA);
 745                         }
 746                         this->DeleteNodeById(node_actualA->GetId());
 747                         this->DeleteEntriesByIdNode(vectorPairNodesScore, 0, node_actualA->GetId());
 748                 }
 749                 else
 750                 {
 751                         double score_fatherB_to_A = this->GetPairNodes_Score(vectorPairNodesScore, node_actualA->GetId(), node_actualB->GetEdge()->GetTarget()->GetId());
 752                         //double score_fatherB_to_A =  pair_fatherB_to_A.second;
 753                         // If the score is lower that the score of the actual pair
 754                         if(score_fatherB_to_A < score_actual)
 755                         {
 756                                 //- The compared leaf is not common. Mark the compared leaf, delete the compared leave and save the leaf as not common.
 757                                 node_actualB->Mark();
 758                                 std::map<unsigned int, Node*>::iterator it;
 759                                 it=nonCommonNodes_B.find(node_actualB->GetId());
 760                                 if(it==nonCommonNodes_B.end())
 761                                         nonCommonNodes_B[node_actualB->GetId()] = node_actualB;
 762                                 //nonCommonNodes_B.push_back(node_actualB);
 763                         }
 764                         // If the score is greater
 765                         else
 766                         {
 767                                 // - The leafs are common then mark both leafs as commons. Mark both leaves, delete them and save them as common leafs.
 768                                 node_actualA->Mark();
 769                                 node_actualB->Mark();
 770                                 //commonNodes_B.push_back(node_actualB);
 771                                 //commonNodes_A.push_back(node_actualA);
 772
 773                                 std::map<unsigned int, Node*>::iterator it;
 774                                 it=commonNodes_B.find(node_actualB->GetId());
 775                                 if(it==commonNodes_B.end())
 776                                         commonNodes_B[node_actualB->GetId()] = node_actualB;
 777
 778                                 it=commonNodes_A.find(node_actualA->GetId());
 779                                 if(it==commonNodes_A.end())
 780                                         commonNodes_A[node_actualA->GetId()] = node_actualA;
 781                         }
 782                         treeB.DeleteNodeById(node_actualB->GetId());
 783                         this->DeleteEntriesByIdNode(vectorPairNodesScore, 1, node_actualB->GetId());
 784                 }
 785         }
 786
 787         // Print the results
 788         std::cout << "Non-CommonNodes [A,B]" <<  nonCommonNodes_A.size() << ", " << nonCommonNodes_B.size() << std::endl;
 789         std::cout << "CommonNodes [A,B]" <<  commonNodes_A.size() << ", " << commonNodes_B.size() << std::endl;
 790
 791         // Get the branches with the lowest difference
 792         /*bool first = true;
 793         double lowestValue = 0;
 794         Pair_Node_Node bestPairNodes;
 795         int iteration = 1;
 796         for(Vector_Pair_PairNodes_Score::iterator it = vectorPairNodesScore.begin(); it != vectorPairNodesScore.end(); ++it )
 797         {
 798                 if(first)
 799                 {
 800                         lowestValue = (*it).second;
 801                         bestPairNodes = (*it).first;
 802                         first = false;
 803                 }
 804                 else if( (*it).second < lowestValue )
 805                 {
 806                         lowestValue = (*it).second;
 807                         bestPairNodes = (*it).first;
 808                 }
 809         }
 810
 811         // Print the best pair of nodes
 812         if(vectorPairNodesScore.size() > 0)
 813                 std::cout << "Best pair: Node A:" << (bestPairNodes).first->GetId() << ", Node B:" << (bestPairNodes).second->GetId() << std::endl;
 814          */
 815
 816         /*// AMP - This code was commented after a meeting with Leonardo Florez
 817         // We make the hypothesis that both roots are equal and
 818         // that the root has only one child.
 819         // The comparison of two branches implies the alignment
 820         // of the mother branches and then the comparison of the
 821         // daughters branches using intrinsic characteristics of
 822         // each branch and a distance function between both.
 823
 824         if (this->IsEmpty() || treeB.IsEmpty())
 825         {
 826                 std::cout << "One or both airways are empty." << std::endl;
 827                 return;
 828         }
 829
 830         Edge actualRootEdge = this->m_root->GetEdge();
 831         Edge comparedRootEdge = treeB.m_root->GetEdge();
 832
 833         // By the hypothesis, both root edges are similar. So we have to compared
 834         // the edged of the root edges.
 835         actualRootEdge.compareChildEdges(comparedRootEdge);
 836          */
 837 }
 838
 839 Pair_PairNodes_Score AirwaysTree::GetBestLeaf(Vector_Pair_PairNodes_Score vectorPairNodesScore)
 840 {
 841         // Variables
 842         Pair_PairNodes_Score pairNodesScore_final;
 843         bool first = true;
 844
 845         for(Vector_Pair_PairNodes_Score::iterator it = vectorPairNodesScore.begin(); it != vectorPairNodesScore.end(); ++it)
 846         {
 847                 // Check that one of the node is a leaf
 848                 pair_nodes pair_actual = (*it).first;
 849                 Node* a = pair_actual.first;
 850                 Node* b = pair_actual.second;
 851                 if( a->IsLeaf() || b->IsLeaf() )
 852                 {
 853                         if(first)
 854                         {
 855                                 first = false;
 856                                 pairNodesScore_final = (*it);
 857                         }
 858                         else if((*it).second < pairNodesScore_final.second)
 859                         {
 860                                 pairNodesScore_final = (*it);
 861                         }
 862                 }
 863         }
 864
 865         if(first)
 866                 std::cout << "Best leaf not found" << std::endl;
 867
 868         return pairNodesScore_final;
 869 }
 870
 871 double AirwaysTree::GetPairNodes_Score(Vector_Pair_PairNodes_Score vectorPairNodesScore, const unsigned int idA, const unsigned int idB)
 872 {
 873         bool found = false;
 874         double score = -1;
 875
 876         for(Vector_Pair_PairNodes_Score::iterator it = vectorPairNodesScore.begin(); it != vectorPairNodesScore.end(); ++it)
 877         {
 878                 // Check that one of the node is a leaf
 879                 pair_nodes pair_actual = (*it).first;
 880                 Node* a = pair_actual.first;
 881                 Node* b = pair_actual.second;
 882                 if( a->GetId() == idA && b->GetId() == idB )
 883                 {
 884                         found = true;
 885                         score = (*it).second;
 886                 }
 887         }
 888
 889         if(!found)
 890                 std::cout << "Pair of Ids not found: " << idA << ", " << idB << std::endl;
 891
 892         return score;
 893 }
 894
 895 bool AirwaysTree::DeleteNodeById(unsigned int nId)
 896 {
 897         bool deleted = false;
 898         if(m_root->GetId() == nId)
 899         {
 900                 delete(m_root);
 901                 m_root = NULL;
 902                 return true;
 903         }
 904
 905         return this->m_root->DeleteNodeById(nId);
 906 }
 907
 908 unsigned int AirwaysTree::DeleteEntriesByIdNode(Vector_Pair_PairNodes_Score& vectorPairNodesScore, unsigned int component, unsigned int nId)
 909 {
 910         unsigned int pairs_deleted = 0;
 911
 912         Vector_Pair_PairNodes_Score::iterator it = vectorPairNodesScore.begin();
 913         while(it != vectorPairNodesScore.end())
 914         {
 915                 // Check that one of the node is a leaf
 916                 pair_nodes pair_actual = (*it).first;
 917                 unsigned int id_actual = -1;
 918                 if(component == 0)
 919                         id_actual = pair_actual.first->GetId();
 920                 else
 921                         id_actual = pair_actual.second->GetId();
 922
 923                 if( id_actual == nId)
 924                 {
 925                         it = vectorPairNodesScore.erase(it);
 926                         ++pairs_deleted;
 927                 }
 928                 else
 929                         ++it;
 930         }
 931
 932         std::cout << "Pairs deleted: " << pairs_deleted << std::endl;
 933
 934         return pairs_deleted;
 935 }
 936
 937 AirwaysTree* AirwaysTree::GetCopy()
 938 {
 939         return NULL;
 940 }
 941
 942 Vector_Pair_PairNodes_Score AirwaysTree::CompareAllBranches(AirwaysTree& treeB)
 943 {
 944         std::cout << "Compare all branches  ..." << std::endl;
 945         // Variables
 946         int nodesA = this->m_nodes.size();
 947         int nodesB = treeB.m_nodes.size();
 948
 949         Vector_Pair_PairNodes_Score vectorPairNodesScore;
 950
 951         // For all the actual branches ( nodes different from the root)
 952         for(unsigned int idA = 2; idA <= nodesA; ++idA)
 953         {
 954                 // Get the actual branch in A
 955                 Edge* branchA = this->GetComposedEdge(idA);
 956
 957                 //Get the final actual node in A
 958                 Node* node_actualA = this->GetNodeById(idA);
 959
 960                 std::cout << "IdA: " << idA << std::endl;
 961
 962                 // For all the branches in B tree
 963                 for(unsigned int idB = 2; idB <= nodesB; ++idB)
 964                 {
 965                         // Get the actual branch in B
 966                         Edge* branchB = treeB.GetComposedEdge(idB);
 967
 968                         //Get the final actual node in B
 969                         Node* node_actualB = treeB.GetNodeById(idB);
 970
 971                         // Get the score from A to B and from B to A
 972                         double scoreA2B = branchA->CompareWith(branchB);
 973                         double scoreB2A = branchB->CompareWith(branchA);
 974
 975                         // Get the final score
 976                         // Options:
 977                         // A. Maximum score
 978                         //double score = scoreB2A > scoreA2B ? scoreB2A : scoreA2B ;
 979
 980                         // B. Sum of scores
 981                         double score = scoreB2A + scoreA2B;
 982
 983                         // Build the pair of nodes
 984                         pair_nodes actualPairNodes(node_actualA, node_actualB);
 985
 986                         Pair_PairNodes_Score actualPairNodesScore(actualPairNodes, score);
 987                         vectorPairNodesScore.push_back(actualPairNodesScore);
 988
 989                         //std::cout << "IdA: " << idA << ", IdB: " << idB << ", score: " << score << ", ... OK" << std::endl;
 990                         delete(branchB);
 991                 }
 992                 delete(branchA);
 993         }
 994
 995         std::cout << "Compare all branches  ... OK" << std::endl;
 996         return vectorPairNodesScore;
 997 }
 998
 999 Edge* AirwaysTree::GetComposedEdge(unsigned int idNode)
1000 {
1001         return this->m_root->GetComposedEdge(idNode);
1002 }
1003
1004 AirwaysTree* AirwaysTree::GetSingleDetachedTreeNodeById(unsigned int nId)
1005 {
1006         AirwaysTree* tree_detached = NULL;
1007         Node* node_searched = this->GetNodeById(nId);
1008         if(node_searched)
1009         {
1010                 Node* node_root_detached = node_searched->GetDetachedRootCopy();
1011                 tree_detached = new AirwaysTree();
1012                 tree_detached->m_root = node_root_detached;
1013         }
1014         else
1015                 std::cout << "Node not found" << std::endl;
1016         return tree_detached;
1017 }
1018
1019 void AirwaysTree::SubIsomorphism(AirwaysTree& treeB)
1020 {
1021         if (this->IsEmpty() || treeB.IsEmpty())
1022                 return;
1023         //Marking the root node
1024         this->m_root->Mark();
1025         treeB.m_root->Mark();
1026
1027         //iterating children - first step
1028         vec_nodes::const_iterator it = this->m_root->GetChildren().begin();
1029         for (; it != this->m_root->GetChildren().end(); ++it)
1030         {
1031                 unsigned int totalCost = 0;
1032                 Node* currentB = NULL;
1033                 Node* currentA = *it;
1034                 if (currentA->IsMarked())
1035                         continue;
1036                 vec_nodes::const_iterator it2 = treeB.m_root->GetChildren().begin();
1037                 for (; it2 != treeB.m_root->GetChildren().end(); ++it2)
1038                 {
1039                         if ((*it2)->IsMarked())
1040                                 continue;
1041                         unsigned int currentCost = (*it)->GetCost(*it2);
1042                         if (totalCost < currentCost)
1043                         {
1044                                 currentB = *it2;
1045                                 totalCost = currentCost;
1046                         } //if
1047                 } //for
1048                 if (currentB != NULL)
1049                         if (*currentA == *currentB)
1050                                 this->GetSubTreeByDetachingNode(currentA).SubIsomorphism(
1051                                                 treeB.GetSubTreeByDetachingNode(currentB));
1052         } //for
1053         //end of first step
1054         //second step
1055         it = this->m_root->GetChildren().begin();
1056         for (; it != this->m_root->GetChildren().end(); ++it)
1057         {
1058                 if ((*it)->IsMarked())
1059                         continue;
1060                 vec_nodes::const_iterator it2 = treeB.m_root->GetChildren().begin();
1061                 for (; it2 != treeB.m_root->GetChildren().end(); ++it2)
1062                 {
1063                         Node* currentB = *it2;
1064                         Node* currentA = *it;
1065                         if (*currentB > *currentA)
1066                                 continue;
1067                         if (currentA->FindSimilarEdgeByAccumulation(currentB))
1068                         {
1069                                 AirwaysTree subTreeA = this->GetSubTreeByDetachingNode(
1070                                                 currentA);
1071                                 AirwaysTree subTreeB = treeB.GetSubTreeByDetachingNode(
1072                                                 currentB);
1073                                 subTreeB.UpdateLevels((*it2)->GetLevel());
1074                                 subTreeA.SubIsomorphism(subTreeB);
1075                                 treeB.m_root->UpdateLevels(treeB.m_root->GetLevel());
1076                                 currentA->Mark();
1077                                 if (!(*it2)->MarkPath(currentB))
1078                                         std::cout << "Error finding path" << std::endl;
1079                                 break;
1080                         } //if
1081                 } //for
1082         } //for
1083         vec_nodes::const_iterator it2 = treeB.m_root->GetChildren().begin();
1084         for (; it2 != treeB.m_root->GetChildren().end(); ++it2)
1085         {
1086                 if ((*it2)->IsMarked())
1087                         continue;
1088                 it = this->m_root->GetChildren().begin();
1089                 for (; it != this->m_root->GetChildren().end(); ++it)
1090                 {
1091                         Node* currentB = *it2;
1092                         Node* currentA = *it;
1093                         if (*currentA > *currentB)
1094                                 continue;
1095                         if (currentB->FindSimilarEdgeByAccumulation(currentA))
1096                         {
1097                                 AirwaysTree subTreeA = this->GetSubTreeByDetachingNode(
1098                                                 currentA);
1099                                 AirwaysTree subTreeB = treeB.GetSubTreeByDetachingNode(
1100                                                 currentB);
1101                                 subTreeA.UpdateLevels((*it)->GetLevel());
1102                                 subTreeA.SubIsomorphism(subTreeB);
1103                                 this->m_root->UpdateLevels(this->m_root->GetLevel());
1104                                 currentB->Mark();
1105                                 if (!(*it)->MarkPath(currentA))
1106                                         std::cout << "Error finding path" << std::endl;
1107                                 break;
1108                         } //if
1109                 } //for
1110         }
1111         //End of second step
1112 }
1113
1114 bool AirwaysTree::Isomorphism(AirwaysTree& tree)
1115 {
1116         if ((this->GetWeight() != tree.GetWeight())
1117                         || (this->GetHeight() != tree.GetHeight()))
1118                 return false;
1119         return this->m_root->CompareNodes(tree.m_root);
1120 }
1121
1122 void AirwaysTree::ImageReconstructionFromNode(TInputImage::Pointer& result,     Node* node, bool skeleton)
1123 {
1124         if (node == NULL)
1125                 return;
1126         //Finding the center of the sphere
1127         const Edge* edge = node->GetEdge();
1128         if (edge == NULL)
1129                 return;
1130         //making a copy of the current image
1131         TInputImage::Pointer copy;
1132         DuplicatorType::Pointer duplicator = DuplicatorType::New();
1133         if (skeleton)
1134                 duplicator->SetInputImage(this->m_skeleton);
1135         else
1136                 duplicator->SetInputImage(this->m_img);
1137         duplicator->Update();
1138         copy = duplicator->GetOutput();
1139         //First remove the node brothers -- the following code in a way does not
1140         // make sense but I'll just comment it if it is necessary
1141         /*Node* parent = edge->m_source;
1142      NodeVector::iterator pIt = parent->m_children.begin();
1143      for (; pIt != parent->m_children.end(); ++pIt)
1144      if ((*pIt) != node)
1145      this->RemoveBranchFromImage(copy, node);
1146      //end of brother removal*/
1147         pair_posVox_rad center;
1148         vec_pair_posVox_rad::const_iterator it = edge->m_vec_pair_posVox_rad.begin();
1149         for (; it != edge->m_vec_pair_posVox_rad.end(); ++it)
1150                 if (edge->m_source->GetCoords() == (*it).first)
1151                         center = *it;
1152         //converting to index
1153         TInputImage::PointType point;
1154         point[0] = center.first[0];
1155         point[1] = center.first[1];
1156         point[2] = center.first[2];
1157         TInputImage::IndexType indCenter;
1158         if (skeleton)
1159                 this->m_skeleton->TransformPhysicalPointToIndex(point, indCenter);
1160         else
1161                 this->m_img->TransformPhysicalPointToIndex(point, indCenter);
1162         //creating itk sphere
1163         SphereType::Pointer sphere = SphereType::New();
1164         sphere->SetRadius(center.second * 4);
1165         SphereType::InputType sphereCenter;
1166         sphereCenter[0] = indCenter[0];
1167         sphereCenter[1] = indCenter[1];
1168         sphereCenter[2] = indCenter[2];
1169         sphere->SetCenter(sphereCenter);
1170
1171         //testing code
1172         //      std::cout<<"center = "<< indCenter << std::endl;
1173         //std::cout << "Radio = " << center.second << std::endl;
1174         TInputImage::SizeType regionSize;
1175         regionSize[0] = 8 * center.second;
1176         regionSize[1] = 8 * center.second;
1177         regionSize[2] = 8 * center.second;
1178         //std::cout << "RegionSize = " << regionSize << std::endl;
1179
1180         TInputImage::IndexType regionIndex;
1181         regionIndex[0] = indCenter[0] - center.second * 4;
1182         regionIndex[1] = indCenter[1] - center.second * 4;
1183         regionIndex[2] = indCenter[2] - center.second * 4;
1184         //      std::cout << "regionIndex = " << regionIndex << std::endl;
1185
1186         TInputImage::RegionType region;
1187         region.SetSize(regionSize);
1188         region.SetIndex(regionIndex);
1189         RegionIterator it2(copy, region);
1190         while (!it2.IsAtEnd())
1191         {
1192                 TInputImage::IndexType ind = it2.GetIndex();
1193                 SphereType::InputType point;
1194                 point[0] = ind[0];
1195                 point[1] = ind[1];
1196                 point[2] = ind[2];
1197                 SphereType::OutputType output = sphere->Evaluate(point);
1198                 if (output)
1199                         it2.Set(0);
1200                 ++it2;
1201         }
1202
1203         /*
1204          * First region growing - Segmenting the wanted branch to be removed
1205          * Removing the branch from the image
1206          */
1207         if (!ComputeSimpleRegionGrowing(copy, this->m_root->GetCoords()))
1208                 return;
1209         SubtractImageFilterType::Pointer subtractFilter =
1210                         SubtractImageFilterType::New();
1211         if (skeleton)
1212                 subtractFilter->SetInput1(this->m_skeleton);
1213         else
1214                 subtractFilter->SetInput1(this->m_img);
1215         subtractFilter->SetInput2(copy);
1216         subtractFilter->Update();
1217         TInputImage::Pointer diff = subtractFilter->GetOutput();
1218         /*
1219          * Second region growing - Cleaning the resulting image
1220          */
1221         Node* leaf = this->GetLeaf(node);
1222         if (ComputeSimpleRegionGrowing(diff, leaf->GetCoords()))
1223                 result = diff;
1224         else
1225                 std::cout << "problemas" << std::endl;
1226         //cin.ignore().get(); //Pause Command for Linux Terminal
1227
1228 }
1229
1230 void AirwaysTree::ImageReconstruction(TInputImage::Pointer& result, bool skeleton, Node* node)
1231 {
1232         if (node == NULL)
1233         {
1234                 DuplicatorType::Pointer duplicator = DuplicatorType::New();
1235                 if (skeleton)
1236                         duplicator->SetInputImage(this->m_skeleton);
1237                 else
1238                         duplicator->SetInputImage(this->m_img);
1239                 duplicator->Update();
1240                 result = duplicator->GetOutput();
1241                 node = this->m_root;
1242         } //if
1243         if (!node->IsMarked())
1244         {
1245                 this->RemoveBranchFromImage(result, node);
1246                 return;
1247         }
1248         vec_nodes::const_iterator it = node->GetChildren().begin();
1249         for (; it != node->GetChildren().end(); ++it)
1250                 this->ImageReconstruction(result, skeleton, *it);
1251
1252 }
1253
1254 void AirwaysTree::ImageReconstructionBySpheres(TInputImage::Pointer& result, Node* node, bool useMarks)
1255 {
1256         if (node == NULL)
1257         {
1258                 node = this->m_root;
1259                 DuplicatorType::Pointer duplicator = DuplicatorType::New();
1260                 duplicator->SetInputImage(m_img);
1261                 duplicator->Update();
1262                 result = duplicator->GetOutput();
1263                 //Removing white pixels -- cleaning image
1264                 TInputImage::RegionType lRegion = result->GetLargestPossibleRegion();
1265                 RegionIterator it(result, lRegion);
1266                 while (!it.IsAtEnd())
1267                 {
1268                         it.Set(0);
1269                         ++it;
1270                 } //while
1271         } //if
1272         if (node->IsMarked() || !useMarks)
1273         {
1274                 this->CreateSpheres(result, node);
1275                 vec_nodes::const_iterator it = node->GetChildren().begin();
1276                 for (; it != node->GetChildren().end(); ++it)
1277                         this->ImageReconstructionBySpheres(result, *it);
1278         } //if
1279 }
1280
1281 Node* AirwaysTree::GetNode(const Vec3 nodeCoords, Node* current)
1282 {
1283         if (current == NULL)
1284                 current = this->m_root;
1285         if (current->GetCoords() == nodeCoords)
1286                 return current;
1287         vec_nodes::const_iterator it = current->GetChildren().begin();
1288         for (; it != current->GetChildren().end(); ++it)
1289         {
1290                 Node* child = GetNode(nodeCoords, *it);
1291                 if (child != NULL)
1292                         return child;
1293         } //rof
1294         return NULL;
1295 }
1296
1297 Node* AirwaysTree::GetLeaf(Node* current)
1298 {
1299         if (current == NULL)
1300                 current = this->m_root;
1301         if (current->IsLeaf())
1302                 return current;
1303         vec_nodes::const_iterator it = current->GetChildren().begin();
1304         Node* node = NULL;
1305         unsigned int maxHeight = 0;
1306         for (; it != current->GetChildren().end(); ++it)
1307         {
1308                 unsigned int cHeight = this->GetHeight(*it);
1309                 if (maxHeight <= cHeight)
1310                 {
1311                         maxHeight = cHeight;
1312                         node = *it;
1313                 } //if
1314         } //for
1315         if (node != NULL)
1316                 return this->GetLeaf(node);
1317         return NULL;
1318 }
1319
1320 void AirwaysTree::UpdateEdges(Node* node)
1321 {
1322         if (node == NULL)
1323                 node = this->m_root;
1324
1325         if (node->GetEdge() != NULL)
1326         {
1327                 Edge* edge = node->GetEdge();
1328                 //edge->m_skInfoPairVector = this->GetSkeletonInfoFromEdge(edge->m_source->GetCoords(), edge->m_target->GetCoords());
1329                 edge->SetSkeletonPairVector( this->GetSkeletonInfoFromEdge(edge->m_source->GetCoords(), edge->m_target->GetCoords()) );
1330                 node->GetEdge()->UpdateEdgeInfo();
1331         }
1332         vec_nodes::const_iterator it = node->GetChildren().begin();
1333         for (; it != node->GetChildren().end(); ++it)
1334                 this->UpdateEdges(*it);
1335 }
1336
1337 AirwaysTree& AirwaysTree::GetSubTreeByDetachingNode(Node* node)
1338 {
1339         AirwaysTree* tree = new AirwaysTree();
1340         tree->m_root = node;
1341         return *tree;
1342 }
1343
1344 DiGraph_EdgePair AirwaysTree::AddBoostEdge(const pair_posVox_rad& source, const pair_posVox_rad& target, DiGraph& graph)
1345 {
1346         bool vFound[2] = { false, false };
1347         DiGraph_VertexDescriptor vDescriptor[2];
1348
1349         //Remember that DiGraph_VertexIteratorPair, first = begin, second = end
1350         DiGraph_VertexIndexMap index = get(boost::vertex_index, graph);
1351
1352         //Search the vertex in the graph
1353         for (DiGraph_VertexIteratorPair vPair = boost::vertices(graph); vPair.first != vPair.second && (!vFound[0] || !vFound[1]); ++vPair.first)
1354         {
1355                 // Get the current information [vector, intensity]
1356                 pair_posVox_rad current = graph[index[*vPair.first]];
1357
1358                 // Check if current is equal to the source or target vertex
1359                 if (current.first == source.first)
1360                 {
1361                         vDescriptor[0] = index[*vPair.first];
1362                         vFound[0] = true;
1363                 } //if
1364                 else if (current.first == target.first)
1365                 {
1366                         vDescriptor[1] = index[*vPair.first];
1367                         vFound[1] = true;
1368                 } //if
1369         } //for
1370
1371         // If one of the vertices is not in the graph, adds the new vertex
1372
1373         // Add the vertices if they are not in the graph
1374         if (!vFound[0])
1375                 vDescriptor[0] = boost::add_vertex(source, graph);
1376         if (!vFound[1])
1377                 vDescriptor[1] = boost::add_vertex(target, graph);
1378
1379         // Check if the edge exist in the graph
1380         // boost::edge returns the pair<edge_descriptor, bool>, with bool=true is the edge exists
1381         DiGraph_EdgePair edgeC1 = boost::edge(vDescriptor[0], vDescriptor[1],graph);
1382         DiGraph_EdgePair edgeC2 = boost::edge(vDescriptor[1], vDescriptor[0],graph);
1383
1384         // If the edge does not exist, then add it
1385         if (!edgeC1.second && !edgeC2.second)
1386         {
1387                 DiGraph_EdgePair nEdge = boost::add_edge(vDescriptor[0], vDescriptor[1], graph);
1388                 return nEdge;
1389         } //if
1390
1391         // If an edge has not been added then -> second = false
1392         // Remove the added vertices because the edge already exists.
1393         if (!vFound[0])
1394                 boost::remove_vertex(vDescriptor[0], graph);
1395
1396         if (!vFound[1])
1397                 boost::remove_vertex(vDescriptor[1], graph);
1398
1399         edgeC1.second = false;
1400         edgeC1.second = false;
1401
1402         return edgeC1;
1403 }
1404
1405 bool AirwaysTree::FindNeighborhood(DiGraph& graph, Vec3Queue& queue, Vec3List& used, const Vec3& end, ConstNeighborhoodIterator& iterator)
1406 {
1407         if (queue.empty())
1408                 return false;
1409
1410         // Get the next pixel in the queue
1411         Vec3 current = queue.front();
1412
1413         // Put the pixel in the used list
1414         used.push_back(current);
1415
1416         // Delete next pixel from the queu
1417         queue.pop();
1418
1419         //AMP//iterator.GoToBegin();
1420         //AMP//while (!iterator.IsAtEnd())
1421         //AMP//{
1422
1423         // AMP commented
1424         /*unsigned int centerIndex = iterator.GetCenterNeighborhoodIndex();
1425                 TInputImage::IndexType index = iterator.GetIndex(centerIndex);
1426                 SKPairInfo src(Vec3(index[0], index[1], index[2]),
1427                                 iterator.GetPixel(iterator.GetCenterNeighborhoodIndex()));
1428                 if (src.first != current)
1429                 {
1430                         ++iterator;
1431                         continue;
1432                 }       //if
1433          */
1434         //PMA commented
1435
1436         // AMP
1437
1438         // Set the current index to the actual pixel
1439         TInputImage::IndexType indexActual;
1440         indexActual[0] = current[0];
1441         indexActual[1] = current[1];
1442         indexActual[2] = current[2];
1443
1444         // Place the iterator in the actual index
1445         iterator.SetLocation(indexActual);
1446
1447         // Create the pair [vector, intensity] for the actual index
1448         pair_posVox_rad src(current, iterator.GetPixel(iterator.GetCenterNeighborhoodIndex()));
1449
1450         // PMA
1451
1452         // Iterate over the neighbors
1453         for (unsigned int i = 0; i < iterator.Size(); i++)
1454         {
1455                 //AMP//if (iterator.GetPixel(i) != 0 && i != centerIndex)
1456
1457                 // If the intensity is not zero and it's the center
1458                 if (iterator.GetPixel(i) != 0 && i != iterator.GetCenterNeighborhoodIndex())
1459                 {
1460                         // Get the index of the iterator
1461                         TInputImage::IndexType ind = iterator.GetIndex(i);
1462
1463                         // Create the target pair [vector, intensity] for the target
1464                         pair_posVox_rad target(Vec3(ind[0], ind[1], ind[2]),    iterator.GetPixel(i));
1465
1466                         // Add the edge to the graph
1467                         bool added = AddBoostEdge(src, target, graph).second;
1468
1469                         // If it could not be added then it is because it already exist, then go to next neighboor
1470                         if (!added)
1471                                 continue;
1472
1473                         // The edge was added. Stop the adding process if the target vertex was reached.
1474                         if (target.first == end)
1475                                 return true;
1476
1477                         // If it the edge was added and it is not the target then add it to the queue, if and only if
1478                         // it was not already used.
1479                         if (std::find(used.begin(), used.end(), target.first)== used.end())
1480                                 queue.push(target.first);
1481                 }//if
1482         }//for
1483         //AMP//break;
1484         //AMP//}        //while
1485         return false;
1486 }
1487
1488 DiGraph AirwaysTree::GetGraphShortestPath(const DiGraph& graph, const Vec3& begin, const Vec3& end)
1489 {
1490         // Variables
1491         DiGraph ret; // Returning graph with the shortest path
1492         DiGraph_EdgeIterator ei, eo; // Initial and final iterators
1493
1494         // Iterate over the edges of the graph
1495         boost::tie(ei, eo) = boost::edges(graph);
1496         for (; ei != eo; ++ei)
1497         {
1498                 // Get the source of the actual edge
1499                 pair_posVox_rad src = graph[boost::source(*ei, graph)];
1500
1501                 // If the source corresponds to the begin pixel
1502                 if (src.first == begin)
1503                 {
1504                         // Get the target of the actual edge
1505                         pair_posVox_rad target = graph[boost::target(*ei, graph)];
1506
1507                         // If the actual target corresponds to the end pixel
1508                         if (target.first == end)
1509                         {
1510                                 // Add the edge to the final path
1511                                 AddBoostEdge(src, target, ret);
1512
1513                                 //AMP
1514                                 std::cout << "Add final target [src]-[ target]:[" << src.first << "]-[" << target.first << "], numVertices:" << boost::num_vertices(ret) << std::endl;
1515                                 //PMA
1516
1517                                 // Return the path
1518                                 return ret;
1519                         }
1520
1521                         // If the actual target does not correspond to the end pixel, then add it to the returning path
1522                         // and run the shortest path with the target of the actual edge.
1523                         ret = GetGraphShortestPath(graph, target.first, end);
1524
1525                         // If the path is empty then the iteration continues
1526                         if (boost::num_vertices(ret) == 0)
1527                                 continue;
1528
1529                         // if not return ret + new edge
1530                         // If the path was found then add the actual edge and return the path
1531                         AddBoostEdge(src, target, ret);
1532                         //AMP
1533                         //std::cout << "Add found edge [src]-[ target]:[" << src.first << "]-[" << target.first << "], numVertices:" << boost::num_vertices(ret) << std::endl;
1534                         //PMA
1535                         return ret;
1536
1537                 }//fi
1538         }
1539         return ret;
1540 }
1541
1542 vec_pair_posVox_rad AirwaysTree::GetGraphShortestPath_AMP(const DiGraph& graph, const Vec3& begin, const Vec3& end)
1543 {
1544         //std::cout << "SKPairInfoVector AirwaysTree::GetGraphShortestPath_AMP(const DiGraph& graph, const Vec3& begin, const Vec3& end)" << std::endl;
1545         // Variables
1546         vec_pair_posVox_rad path;
1547         DiGraph ret; // Returning graph with the shortest path
1548         DiGraph_EdgeIterator ei, eo; // Initial and final iterators
1549
1550         // Iterate over the edges of the graph
1551         boost::tie(ei, eo) = boost::edges(graph);
1552         for (; ei != eo; ++ei)
1553         {
1554                 // Get the source of the actual edge
1555                 pair_posVox_rad src = graph[boost::source(*ei, graph)];
1556
1557                 // If the source corresponds to the begin pixel
1558                 if (src.first == begin)
1559                 {
1560                         // Get the target of the actual edge
1561                         pair_posVox_rad target = graph[boost::target(*ei, graph)];
1562
1563                         // If the actual target corresponds to the end pixel
1564                         if (target.first == end)
1565                         {
1566                                 // Add the edge to the final path
1567                                 AddBoostEdge(src, target, ret);
1568                                 path.push_back(target);
1569                                 path.push_back(src);
1570
1571                                 //AMP
1572                                 //std::cout << "Add final target [src]-[ target]:[" << src.first << "]-[" << target.first << "], numVertices:" << boost::num_vertices(ret) << ", path_size:" << path.size() << std::endl;
1573                                 //PMA
1574
1575                                 // Return the path
1576                                 //std::cout << "SKPairInfoVector AirwaysTree::GetGraphShortestPath_AMP(const DiGraph& graph, const Vec3& begin, const Vec3& end) ... FIRST ... OK" << std::endl;
1577                                 return path;
1578                         }
1579
1580                         //std::cout << "Not yet in the target voxel" << std::endl;
1581                         // If the actual target does not correspond to the end pixel, then add it to the returning path
1582                         // and run the shortest path with the target of the actual edge.
1583                         path = GetGraphShortestPath_AMP(graph, target.first, end);
1584                         //std::cout << "Path found, path size:" << path.size() << std::endl;
1585
1586                         // If the path is empty then the iteration continues
1587                         if (path.size() == 0)
1588                                 continue;
1589
1590                         // if not return ret + new edge
1591                         // If the path was found then add the actual edge and return the path
1592                         AddBoostEdge(src, target, ret);
1593                         path.push_back(src);
1594                         //AMP
1595                         //std::cout << "Add found edge [src]-[ target]:[" << src.first << "]-[" << target.first << "], numVertices:" << boost::num_vertices(ret)  << ", path_size:" << path.size() << std::endl;
1596                         //PMA
1597                         //std::cout << "SKPairInfoVector AirwaysTree::GetGraphShortestPath_AMP(const DiGraph& graph, const Vec3& begin, const Vec3& end) ... SECOND ... OK" << std::endl;
1598                         return path;
1599
1600                 }//fi
1601         }
1602
1603         //std::cout << "SKPairInfoVector AirwaysTree::GetGraphShortestPath_AMP(const DiGraph& graph, const Vec3& begin, const Vec3& end) ... END ... OK" << std::endl;
1604         return path;
1605 }
1606
1607 vec_pair_posVox_rad AirwaysTree::GetSkeletonInfoFromEdge(const Vec3& source, const Vec3& target)
1608 {
1609         //converting to pixel
1610         TInputImage::PointType pointSource;
1611         pointSource[0] = source[0];
1612         pointSource[1] = source[1];
1613         pointSource[2] = source[2];
1614
1615         TInputImage::PointType pointTarget;
1616         pointTarget[0] = target[0];
1617         pointTarget[1] = target[1];
1618         pointTarget[2] = target[2];
1619
1620         TInputImage::IndexType indexSource;
1621         TInputImage::IndexType indexTarget;
1622
1623         this->m_skeleton->TransformPhysicalPointToIndex(pointSource, indexSource);
1624         this->m_skeleton->TransformPhysicalPointToIndex(pointTarget, indexTarget);
1625         //end of conversion
1626
1627         // AMP // Print the initial and final voxels to verify their positions
1628         //std::cout << "Begin: [" << indexSource[0] << "," << indexSource[1] << "," << indexSource[2] << "] -- " << this->m_skeleton->GetPixel(indexSource) << std::endl;
1629         //std::cout << "End: [" << indexTarget[0] << "," << indexTarget[1] << "," << indexTarget[2] << "] -- " << this->m_skeleton->GetPixel(indexTarget) << std::endl;
1630         //std::cout << "BeginVoxel: [" << pointSource[0] << "," << pointSource[1] << "," << pointSource[2] << "]"<< std::endl;
1631         //std::cout << "EndVoxel: [" << pointTarget[0] << "," << pointTarget[1] << "," << pointTarget[2] << "]" << std::endl;
1632         // PMA //
1633
1634         Vec3 sourcePxl(indexSource[0], indexSource[1], indexSource[2]);
1635         Vec3 targetPxl(indexTarget[0], indexTarget[1], indexTarget[2]);
1636
1637         // AMP // Commented
1638         /*
1639         // Compute the distance in each direction and add "21" in each one
1640         float dstX = beginPxl.EculideanDistance(
1641                         Vec3(endPxl[0], beginPxl[1], beginPxl[2])) + 21;
1642         float dstY = beginPxl.EculideanDistance(
1643                         Vec3(beginPxl[0], endPxl[1], beginPxl[2])) + 21;
1644         float dstZ = beginPxl.EculideanDistance(
1645                         Vec3(beginPxl[0], beginPxl[1], endPxl[2])) + 21;
1646
1647         //start point
1648         TInputImage::IndexType index;
1649         index[0] = beginPxl[0] < endPxl[0] ? beginPxl[0] : endPxl[0];
1650         index[1] = beginPxl[1] < endPxl[1] ? beginPxl[1] : endPxl[1];
1651         index[2] = beginPxl[2] < endPxl[2] ? beginPxl[2] : endPxl[2];
1652
1653         index[0] -= index[0] <= 10 ? 0 : 10;
1654         index[1] -= index[1] <= 10 ? 0 : 10;
1655         index[2] -= index[2] <= 10 ? 0 : 10;
1656
1657
1658         //size of region
1659         TInputImage::SizeType size;
1660         size[0] = ceil(dstX);
1661         size[1] = ceil(dstY);
1662         size[2] = ceil(dstZ);
1663
1664         TInputImage::RegionType region;
1665         region.SetSize(size);
1666         region.SetIndex(index);
1667
1668         //ConstNeighborhoodIterator iterator(radius, this->m_skeleton, region);
1669         // PMA // Commented
1670          */
1671
1672         TInputImage::SizeType radius;
1673         radius[0] = 1;
1674         radius[1] = 1;
1675         radius[2] = 1;
1676         ConstNeighborhoodIterator iterator(radius, this->m_skeleton, this->m_skeleton->GetLargestPossibleRegion());
1677
1678         DiGraph graph;
1679         Vec3Queue queue;
1680         Vec3List used;
1681
1682         // Set the first pixel as starting pixel
1683         queue.push(sourcePxl);
1684
1685         //FindNeighborhood(graph, queue, used, endPxl, iterator);
1686
1687         // Add to the graph all the connected voxel until the ending pixel is found
1688         bool endFound = false;
1689         while (!queue.empty() && !endFound)
1690                 endFound = FindNeighborhood(graph, queue, used, targetPxl, iterator);
1691
1692         if(!endFound)
1693                 std::cout << "***************** END NOT FOUND!!!!" << std::endl;
1694         else
1695                 std::cout << "***************** END FOUND!!!! - Graph Vertices:" << boost::num_vertices(graph) << ", Graph Edges:" << boost::num_edges(graph) << std::endl;
1696
1697         // Get the shortest path between the begin and end pixels
1698         //std::cout << "GetGraphShortestPath [from]-[to]:[" << sourcePxl << "]-[" << targetPxl << "]" << std::endl;
1699         //DiGraph result = GetGraphShortestPath(graph, sourcePxl, targetPxl);
1700         vec_pair_posVox_rad path = GetGraphShortestPath_AMP(graph, sourcePxl, targetPxl);
1701         //std::cout << "GetGraphShortestPath ... OK , numVertices:" << boost::num_vertices(result) << std::endl;
1702
1703         //AMP
1704         bool beginPixelFound = false;
1705         bool endPixelFound = false;
1706         vec_pair_posVox_rad vector;
1707
1708         for(vec_pair_posVox_rad::reverse_iterator it_path = path.rbegin(); it_path != path.rend(); ++it_path)
1709         {
1710                 TInputImage::IndexType ind;
1711                 ind[0] = (*it_path).first[0];
1712                 ind[1] = (*it_path).first[1];
1713                 ind[2] = (*it_path).first[2];
1714                 //std::cout << "["<< (*it_path).first << "], ";
1715
1716                 // To avoid precision problems
1717                 if (indexSource == ind)
1718                 {
1719                         beginPixelFound = true;
1720                         (*it_path).first = source;
1721                 }
1722                 else if (indexTarget == ind)
1723                 {
1724                         endPixelFound = true;
1725                         (*it_path).first = target;
1726                 }
1727                 else
1728                 {
1729                         TInputImage::PointType pnt;
1730                         this->m_skeleton->TransformIndexToPhysicalPoint(ind, pnt);
1731                         (*it_path).first = Vec3(pnt[0], pnt[1], pnt[2]);
1732                 }
1733                 vector.push_back((*it_path));
1734         }
1735         //vector = path;
1736         //PMA
1737
1738         // Variables to iterate over the shortest path and save the results
1739         /*// AMP
1740         DiGraph_VertexIterator i, e;
1741         SKPairInfoVector vector;
1742
1743         // Iterate over the shortest path and add each vertex in real coordinates
1744         boost::tie(i, e) = boost::vertices(result);
1745         std::cout << "Inserting ... [begin,end]: [" << begin << "], [" << end << "]" << std::endl;
1746         for (; i != e; ++i)
1747         //for (; e != i; --e)
1748         {
1749                 // Get the actual vertex
1750                 SKPairInfo vertex = result[*i];
1751                 //SKPairInfo vertex = result[*e];
1752                 std::cout << "["<< vertex.first << "], ";
1753
1754                 // Get the index (voxel position) of the actual vertex
1755                 TInputImage::IndexType ind;
1756                 ind[0] = vertex.first[0];
1757                 ind[1] = vertex.first[1];
1758                 ind[2] = vertex.first[2];
1759
1760                 // To avoid precision problems
1761                 if (indexBegin == ind)
1762                 {
1763                         beginPixelFound = true;
1764                         vertex.first = begin;
1765                 }
1766                 else if (indexEnd == ind)
1767                 {
1768                         endPixelFound = true;
1769                         vertex.first = end;
1770                 }
1771                 else
1772                 {
1773                         TInputImage::PointType pnt;
1774                         this->m_skeleton->TransformIndexToPhysicalPoint(ind, pnt);
1775                         vertex.first = Vec3(pnt[0], pnt[1], pnt[2]);
1776                 }
1777                 vector.push_back(vertex);
1778         }
1779         std::cout << "Inserting ... OK" << std::endl;
1780
1781
1782          *///AMP
1783
1784
1785         if(!beginPixelFound)
1786                 std::cout << "Begin not found: "<< source << std::endl;
1787         if(!endPixelFound)
1788                 std::cout << "End not found:" << target << std::endl;
1789
1790         return vector;
1791 }
1792
1793 bool AirwaysTree::ComputeSimpleRegionGrowing(TInputImage::Pointer& img,
1794                 const Vec3 seedPtn)
1795 {
1796         TInputImage::PointType targetPtn;
1797         targetPtn[0] = seedPtn[0];
1798         targetPtn[1] = seedPtn[1];
1799         targetPtn[2] = seedPtn[2];
1800         TInputImage::IndexType seed;
1801         img->TransformPhysicalPointToIndex(targetPtn, seed);
1802         if (img->GetPixel(seed) == 0)
1803         {
1804                 //std::cout << "seed is 0" << std::endl;
1805                 return false;
1806         }
1807         ConnectedFilterType::Pointer neighborhoodConnected =
1808                         ConnectedFilterType::New();
1809         neighborhoodConnected->SetLower(1);
1810         neighborhoodConnected->SetUpper(1);
1811         neighborhoodConnected->SetReplaceValue(1);
1812         neighborhoodConnected->SetSeed(seed);
1813         neighborhoodConnected->SetInput(img);
1814         neighborhoodConnected->Update();
1815         img = neighborhoodConnected->GetOutput();
1816         return true;
1817 }
1818
1819 void AirwaysTree::RemoveBranchFromImage(TInputImage::Pointer& img, Node* node)
1820 {
1821         //typedef itk::ImageFileWriter<TInputImage> WriterType;
1822         //Finding the center of the sphere
1823         const Edge* edge = node->GetEdge();
1824         if (edge == NULL)
1825                 return;
1826         pair_posVox_rad center;
1827         vec_pair_posVox_rad::const_iterator it = edge->m_vec_pair_posVox_rad.begin();
1828         for (; it != edge->m_vec_pair_posVox_rad.end(); ++it)
1829                 if (edge->m_source->GetCoords() == (*it).first)
1830                         center = *it;
1831         //converting to index
1832         TInputImage::PointType point;
1833         point[0] = center.first[0];
1834         point[1] = center.first[1];
1835         point[2] = center.first[2];
1836         TInputImage::IndexType indCenter;
1837         img->TransformPhysicalPointToIndex(point, indCenter);
1838         double radius = center.second;
1839         //creating itk sphere
1840         SphereType::Pointer sphere = SphereType::New();
1841         SphereType::InputType sphereCenter;
1842         sphereCenter[0] = indCenter[0];
1843         sphereCenter[1] = indCenter[1];
1844         sphereCenter[2] = indCenter[2];
1845         sphere->SetCenter(sphereCenter);
1846         unsigned int objCount = 0;
1847         unsigned int obj3Count = 0; // to assure 3 Connected obj
1848         //Number of connected components, children_size + parent
1849         unsigned int nCConnex = edge->m_source->GetNumberOfChildren() + 1;
1850         TInputImage::Pointer copy;
1851         //std::cout << "entró" << std::endl;
1852         while (objCount < nCConnex || obj3Count <= nCConnex)
1853         {
1854                 if (objCount >= 3)
1855                         obj3Count++;
1856                 radius += 1.0;
1857                 sphere->SetRadius(radius);
1858                 //making a copy of the current image
1859                 DuplicatorType::Pointer duplicator = DuplicatorType::New();
1860                 duplicator->SetInputImage(img);
1861                 duplicator->Update();
1862                 copy = duplicator->GetOutput();
1863
1864                 //testing code
1865                 //      std::cout<<"center = "<< indCenter << std::endl;
1866                 //std::cout << "Radio = " << center.second << std::endl;
1867                 TInputImage::SizeType regionSize;
1868                 regionSize[0] = 4 * radius;
1869                 regionSize[1] = 4 * radius;
1870                 regionSize[2] = 4 * radius;
1871                 //std::cout << "RegionSize = " << regionSize << std::endl;
1872                 TInputImage::SizeType size = img->GetLargestPossibleRegion().GetSize();
1873                 TInputImage::IndexType regionIndex;
1874                 regionIndex[0] = indCenter[0] - (radius * 2);
1875                 if (regionIndex[0] < 0)
1876                         regionIndex[0] = indCenter[0];
1877                 else if (regionIndex[0] > (unsigned int) size[0])
1878                         regionIndex[0] = size[0];
1879                 regionIndex[1] = indCenter[1] - (radius * 2);
1880                 if (regionIndex[1] < 0)
1881                         regionIndex[1] = indCenter[1];
1882                 else if (regionIndex[1] > (unsigned int) size[1])
1883                         regionIndex[1] = size[1];
1884                 regionIndex[2] = indCenter[2] - (radius * 2);
1885                 if (regionIndex[2] < 0)
1886                         regionIndex[2] = indCenter[2];
1887                 else if (regionIndex[2] > (unsigned int) size[2])
1888                         regionIndex[2] = size[2];
1889                 int diff[3];
1890                 diff[0] = size[0] - (regionSize[0] + regionIndex[0]);
1891                 diff[1] = size[1] - (regionSize[1] + regionIndex[1]);
1892                 diff[2] = size[2] - (regionSize[2] + regionIndex[2]);
1893                 //std::cout << "(1) Region Index = " << regionIndex << std::endl;
1894                 //std::cout << "(1) Region size = " << regionSize << std::endl;
1895                 if (diff[0] < 0)
1896                         regionSize[0] = size[0] - regionIndex[0];
1897                 if (diff[1] < 0)
1898                         regionSize[1] = size[1] - regionIndex[1];
1899                 if (diff[2] < 0)
1900                         regionSize[2] = size[2] - regionIndex[2];
1901                 //      std::cout << "regionIndex = " << regionIndex << std::endl;
1902                 //std::cout << "(2) Region Index = " << regionIndex << std::endl;
1903                 //std::cout << "(2) Region size = " << regionSize << std::endl;
1904                 TInputImage::RegionType region;
1905                 region.SetSize(regionSize);
1906                 region.SetIndex(regionIndex);
1907                 RegionIterator it2(copy, region);
1908                 while (!it2.IsAtEnd())
1909                 {
1910                         TInputImage::IndexType ind = it2.GetIndex();
1911                         SphereType::InputType point;
1912                         point[0] = ind[0];
1913                         point[1] = ind[1];
1914                         point[2] = ind[2];
1915                         SphereType::OutputType output = sphere->Evaluate(point);
1916                         if (output)
1917                                 it2.Set(0);
1918                         ++it2;
1919                 }                               //while
1920                 if (node->IsLeaf())
1921                 {
1922                         TInputImage::PointType ptn;
1923                         ptn[0] = node->GetCoords()[0];
1924                         ptn[1] = node->GetCoords()[1];
1925                         ptn[2] = node->GetCoords()[2];
1926                         TInputImage::IndexType tgt;
1927                         img->TransformPhysicalPointToIndex(ptn, tgt);
1928                         if (img->GetPixel(tgt) == 0)
1929                         {
1930                                 std::cout << "It entered here!!!!" << std::endl;
1931                                 return;
1932                         }
1933                 }                               //if
1934                 CastFilterType::Pointer castFilter = CastFilterType::New();
1935                 castFilter->SetInput(copy);
1936                 castFilter->Update();
1937                 ConnectedComponentImageFilterType::Pointer connFilter =
1938                                 ConnectedComponentImageFilterType::New();
1939                 connFilter->SetFullyConnected(true);
1940                 connFilter->SetInput(castFilter->GetOutput());
1941                 connFilter->Update();
1942                 objCount = connFilter->GetObjectCount();
1943                 //test
1944                 /*WriterType::Pointer writer = WriterType::New();
1945          writer->SetInput(copy);
1946          writer->SetFileName("output_antes.mhd");
1947          writer->Update();*/
1948         }                               //while
1949         //std::cout << "salio" << std::endl;
1950         /*
1951          * First region growing - Segmenting the wanted branch to be removed
1952          * Removing the branch from the image
1953          */
1954         Node* leaf = this->GetLeaf(node);
1955         if (!ComputeSimpleRegionGrowing(copy, leaf->GetCoords()))
1956                 return;
1957         SubtractImageFilterType::Pointer subtractFilter =
1958                         SubtractImageFilterType::New();
1959         subtractFilter->SetInput1(img);
1960         subtractFilter->SetInput2(copy);
1961         subtractFilter->Update();
1962         TInputImage::Pointer diff = subtractFilter->GetOutput();
1963         /*
1964          * Second region growing - Cleaning the resulting image
1965          */
1966         if (ComputeSimpleRegionGrowing(diff, this->m_root->GetCoords()))
1967                 img = diff;
1968         /*else
1969      std::cout << "problemas" << std::endl;*/
1970         //cin.ignore().get(); //Pause Command for Linux Terminal
1971 }
1972
1973 //Metodo obsoleto, ya no nos interesa la reconstruccion por esferas
1974 void AirwaysTree::CreateSpheres(TInputImage::Pointer& img, Node* node)
1975 {
1976         if (node->GetEdge() == NULL)
1977                 return;
1978         const Edge* edge = node->GetEdge();
1979         vec_pair_posVox_rad::const_iterator it = edge->m_vec_pair_posVox_rad.begin();
1980         for (; it != edge->m_vec_pair_posVox_rad.end(); ++it)
1981         {
1982                 //converting to index
1983                 TInputImage::PointType point;
1984                 point[0] = (*it).first[0];
1985                 point[1] = (*it).first[1];
1986                 point[2] = (*it).first[2];
1987                 TInputImage::IndexType indexPoint;
1988                 img->TransformPhysicalPointToIndex(point, indexPoint);
1989                 //creating itk sphere
1990                 SphereType::Pointer sphere = SphereType::New();
1991                 sphere->SetRadius((*it).second);
1992                 SphereType::InputType sphereCenter;
1993                 sphereCenter[0] = indexPoint[0];
1994                 sphereCenter[1] = indexPoint[1];
1995                 sphereCenter[2] = indexPoint[2];
1996                 sphere->SetCenter(sphereCenter);
1997                 //creating region
1998                 TInputImage::SizeType regionSize;
1999                 regionSize[0] = 4 * (*it).second;
2000                 regionSize[1] = 4 * (*it).second;
2001                 regionSize[2] = 4 * (*it).second;
2002                 TInputImage::IndexType regionIndex;
2003                 regionIndex[0] = indexPoint[0] - 2 * ceil((*it).second);
2004                 regionIndex[1] = indexPoint[1] - 2 * ceil((*it).second);
2005                 regionIndex[2] = indexPoint[2] - 2 * ceil((*it).second);
2006                 TInputImage::RegionType region;
2007                 region.SetSize(regionSize);
2008                 region.SetIndex(regionIndex);
2009                 RegionIterator it2(img, region);
2010                 while (!it2.IsAtEnd())
2011                 {
2012                         TInputImage::IndexType ind = it2.GetIndex();
2013                         SphereType::InputType point;
2014                         point[0] = ind[0];
2015                         point[1] = ind[1];
2016                         point[2] = ind[2];
2017                         SphereType::OutputType output = sphere->Evaluate(point);
2018                         if (output)
2019                                 it2.Set(1);
2020                         ++it2;
2021                 }               //while
2022         }               //for
2023 }
2024
2025 void AirwaysTree::printUpToLevel(int level)
2026 {
2027         cout << "Printing up to level: " << level <<endl;
2028         this->m_root->printUpToLevel(level);
2029         cout << "Printing done" <<endl;
2030 }
2031
2032 void AirwaysTree::printNodeAndChildrenIds()
2033 {
2034         cout << "Printing nodes and children: " <<endl;
2035         this->m_root->printNodeAndChildrenIds( );
2036         cout << "Printing nodes and children DONE" <<endl;
2037 }
2038
2039 void AirwaysTree::createQuaternionsUpToLevel(int level)
2040 {
2041         cout << "Creating quaternions up to level: " << level <<endl;
2042         if(level >= 2)
2043         {
2044                 // Print header
2045                 cout << "Level" << " " << "SonPositionX" << " " << "SonPositionY" << " " << "SonPositionZ" << " "
2046                                 "ParentVectorX" << " " << "ParentVectorY" << " " << "ParentVectorZ" << " " <<
2047                                 "ChildVectorX" << " " << "ChildVectorY" << " " << "ChildVectorZ"<< " " <<
2048                                 "Qr" << " " << "Qx" << " " << "Qy" " " << "Qz" " " << "Qtheta" << " " <<
2049                                 "Wx" << " " << "Wy" << " " << "Wz" << " " <<
2050                                 "EuclDistSonParent" << " " <<"EuclDistSonGrandson" << " " <<
2051                                 "RadMeanSonParent" << " " << "RadMeanSonGrandson" <<endl;
2052
2053                 this->m_root->createQuaternionsUpToLevel(level, level);
2054         }
2055         else
2056                 cout << "Quaternions must be created from level 2." <<endl;
2057         cout << "Creating quaternions done" <<endl;
2058 }
2059
2060 void AirwaysTree::getStatisticsBifurcations()
2061 {
2062         cout << "Statistics bifurcations ... " <<endl;
2063
2064         int p[6] = {0,0,0,0,0,0};
2065
2066         this->m_root->getStasticsBifurcations(p);
2067
2068         cout << "Bifurcations:" << endl;
2069         for(int i = 0; i < 6; i++)
2070         {
2071                 cout << i+1 << "->" << p[i] << endl;
2072         }
2073
2074         cout << "Statistics bifurcations done " <<endl;
2075 }
2076
2077 void AirwaysTree::saveAirwayToImage(std::string filename, int dims[3], double spc[3], double nOrigin[3])
2078 {
2079         TInputImage::Pointer imagePointer;
2080         imagePointer = TInputImage::New();
2081
2082         TInputImage::IndexType start;
2083         start[0] = 0;
2084         start[1] = 0;
2085         start[2] = 0;
2086
2087         TInputImage::SizeType size;
2088         size[0] = dims[0];
2089         size[1] = dims[1];
2090         size[2] = dims[2];
2091
2092         TInputImage::RegionType region;
2093         region.SetSize( size );
2094         region.SetIndex( start );
2095
2096         TInputImage::SpacingType spacing;
2097         spacing[0] = spc[0];
2098         spacing[1] = spc[1];
2099         spacing[2] = spc[2];
2100
2101         TInputImage::PointType origin;
2102         origin[0] = nOrigin[0];
2103         origin[1] = nOrigin[1];
2104         origin[2] = nOrigin[2];
2105
2106         imagePointer->SetRegions( region );
2107         imagePointer->SetSpacing( spacing );
2108         imagePointer->SetOrigin( origin );
2109         imagePointer->Allocate();
2110
2111         //this->m_root->saveToImage(_imagePointer);
2112         vec_nodes nodes_children_root = this->m_root->GetChildren();
2113         vec_nodes::const_iterator it_children = nodes_children_root.begin();
2114         for( ; it_children != nodes_children_root.end(); ++it_children)
2115         {
2116                 Vec3 coords_father = (*it_children)->GetEdge()->GetSource()->GetCoords();
2117                 Vec3 coords_child = (*it_children)->GetCoords();
2118                 drawLineInImage(coords_father, coords_child, imagePointer);
2119         }
2120
2121         cout << "Writing ..." << endl;
2122         typedef itk::ImageFileWriter<TInputImage> WriterType;
2123         WriterType::Pointer writer = WriterType::New();
2124         writer->SetInput(imagePointer);
2125         writer->SetFileName(filename);
2126         writer->Update();
2127         cout << "Writing ... done" << endl;
2128 }
2129
2130 void AirwaysTree::drawLineInImage(Vec3 initVoxel, Vec3 endVoxel, TInputImage::Pointer imagePointer)
2131 {
2132         TInputImage::RegionType region = imagePointer->GetLargestPossibleRegion();
2133
2134         TInputImage::SpacingType spacing = imagePointer->GetSpacing();
2135
2136         TInputImage::PointType origin = imagePointer->GetOrigin();
2137
2138         TInputImage::IndexType init;
2139         init[0]=(initVoxel[0]-origin[0])/spacing[0];
2140         init[1]=(initVoxel[1]-origin[1])/spacing[1];
2141         init[2]=(initVoxel[2]-origin[2])/spacing[2];
2142
2143         TInputImage::IndexType end;
2144         end[0]=(endVoxel[0]-origin[0])/spacing[0];
2145         end[1]=(endVoxel[1]-origin[1])/spacing[1];
2146         end[2]=(endVoxel[2]-origin[2])/spacing[2];
2147
2148         TPixel pixel;
2149         pixel = 1;
2150
2151         itk::LineIterator<TInputImage> it(imagePointer, init, end);
2152         it.GoToBegin();
2153         while (!it.IsAtEnd())
2154         {
2155                 it.Set(pixel);
2156                 ++it;
2157         }
2158 }
2159
2160 //To store the graph in a vtk
2161 void AirwaysTree::saveAirwayAsVTK(const std::string filename, bool common)
2162 {
2163         vtkSmartPointer<vtkUnsignedCharArray> colors = vtkSmartPointer<vtkUnsignedCharArray>::New();
2164         colors->SetNumberOfComponents(3);
2165         colors->SetName("Colors");
2166
2167         //pointer of lines and points
2168         vtkSmartPointer<vtkCellArray> lines = vtkSmartPointer<vtkCellArray>::New();
2169         vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
2170
2171         Vec3 root = this->m_root->GetCoords();
2172
2173         //UndirectedGraph
2174
2175         srand(time(NULL));
2176         unsigned int id = 1;
2177
2178         CalculateVTKLinesFromEdges(this->m_root, 0, id, pts, lines, colors, common);
2179
2180         vtkSmartPointer<vtkPolyData> linesPolyData = vtkSmartPointer<vtkPolyData>::New();
2181         linesPolyData->SetPoints(pts);
2182         linesPolyData->SetLines(lines);
2183
2184         linesPolyData->GetCellData()->SetScalars(colors);
2185
2186         // Write the file
2187         vtkSmartPointer<vtkPolyDataWriter> writer = vtkSmartPointer<vtkPolyDataWriter>::New();
2188         writer->SetFileName(filename.c_str());
2189 #if VTK_MAJOR_VERSION <= 5
2190         writer->SetInput(linesPolyData);
2191 #else
2192         writer->SetInputData(linesPolyData);
2193 #endif
2194         writer->Write();
2195 }
2196
2197 //To store the graph in a vtk
2198 void AirwaysTree::CalculateVTKLinesFromEdges(const Node* node, const unsigned int& parentId,
2199                 unsigned int& cId, vtkSmartPointer<vtkPoints>& pts,
2200                 vtkSmartPointer<vtkCellArray>& lines,
2201                 vtkSmartPointer<vtkUnsignedCharArray>& colors, bool common)
2202 {
2203         if (node == NULL)
2204                 return;
2205
2206         // Colors
2207         unsigned char red[3] = { 255, 0, 0 };
2208         unsigned char green[3] = { 0, 255, 0 };
2209         unsigned char blue[3] = { 0, 0, 255 };
2210         unsigned char yellow[3] = { 255, 255, 0 };
2211
2212         const vec_nodes children = node->GetChildren();
2213
2214         for(vec_nodes::const_iterator it = children.begin(); it != children.end(); ++it)
2215         {
2216                 if (!(*it)->IsMarked() && common)
2217                         continue;
2218
2219                 const Edge* edge = (*it)->GetEdge();
2220                 //pts->InsertNextPoint(edge->GetSource()->GetCoords().GetVec3());
2221                 //pts->InsertNextPoint(edge->GetTarget()->GetCoords().GetVec3());
2222
2223                 pts->InsertNextPoint(node->GetCoords().GetVec3());
2224                 pts->InsertNextPoint((*it)->GetCoords().GetVec3());
2225
2226                 // Set color to be used
2227                 int numColor = node->GetLevel() % 4;
2228                 if(numColor == 0)
2229                         colors->InsertNextTupleValue(green);
2230                 else if(numColor == 1)
2231                         colors->InsertNextTupleValue(red);
2232                 else if(numColor == 2)
2233                         colors->InsertNextTupleValue(blue);
2234                 else
2235                         colors->InsertNextTupleValue(yellow);
2236
2237                 vtkSmartPointer<vtkLine> line = vtkSmartPointer<vtkLine>::New();
2238                 line->GetPointIds()->SetId(0, parentId);
2239                 line->GetPointIds()->SetId(1, cId++);
2240                 lines->InsertNextCell(line);
2241                 unsigned int newParent = cId++;
2242                 CalculateVTKLinesFromEdges(*it, newParent, cId, pts, lines, colors,     common);
2243         } //for
2244 }
2245
2246
2247 } /* namespace airways */
2248
2249
2250 #endif
2251