4 * Created on: May 3, 2014
5 * Author: caceres@creatis.insa-lyon.fr
6 * Modifications by: Alfredo Morales Pinzón
9 #ifndef _AIRWAYS_TREE_TXX_
10 #define _AIRWAYS_TREE_TXX_
12 #include "airwaysTree.h"
17 AirwaysTree::AirwaysTree() : m_root(NULL)
23 AirwaysTree::AirwaysTree(TInputImage::Pointer img, TInputImage::Pointer img_skele, Node* root, bool shouldUpdate){
25 this->m_root->SetLevel(0);
27 this->m_skeleton = img_skele;
28 //TODO: Should this calculate levels and populate edges and nodes vectors?
30 if(shouldUpdate) this->UpdateEdges();
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();
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);
50 AirwaysTree::AirwaysTree(TInputImage::Pointer img, TInputImage::Pointer img_sk, Node* root)
53 this->m_root->SetLevel(0);
54 this->m_skeleton = img_sk;
59 AirwaysTree::AirwaysTree(TInputImage::Pointer image_airwaySegmentation, TInputImage::Pointer image_skeletonDM, EdgeVector& vector_edges, Edge* edge_trachea)
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;
64 if(edge_trachea == NULL)
65 std::cout << "The trachea edge is NULL" << std::endl;
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);
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
76 // Insert the trachea edge
77 this->Insert(edge_trachea);
80 unsigned int iniSize = vector_edges.size();
82 while (!vector_edges.empty())
84 bool oneEdgeAdded = false;
85 EdgeVector::iterator it = vector_edges.begin();
86 while (it != vector_edges.end())
88 if (this->Insert(*it))
90 std::cout << "Edges added: [" << (*it)->GetSource()->GetCoords() << "]->[" << (*it)->GetTarget()->GetCoords() << "], missing: " << vector_edges.size()-1 << std::endl;
92 it = vector_edges.erase(it);
98 // If there are edges to be added but they could not be added
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;
106 // AMP: I did not find any explanation for this "exception"
108 if (count == iniSize)
110 std::cout << "handled exception: arwaysTree.txx, AirwaysTree" << std::endl;
117 std::cout << "AirwaysTree::AirwaysTree(TInputImage::Pointer img, TInputImage::Pointer img_sk, EdgeVector& tEdgeVector, Edge* trachea) ... OK" << std::endl;
120 AirwaysTree::~AirwaysTree()
125 Node* AirwaysTree::GetRoot() const
130 std::string AirwaysTree::GetImageName() const
132 return (this->m_img_name);
135 std::string AirwaysTree::GetResultPath() const
137 return (this->m_result_path);
140 TInputImage::Pointer AirwaysTree::GetSegmentedImage() const
145 TInputImage::Pointer AirwaysTree::GetSkeletonImage() const
147 return this->m_skeleton;
150 const Node* AirwaysTree::GetNode(const Vec3 nodeCoords) const
152 return this->GetNode(nodeCoords);
155 Node* AirwaysTree::GetNodeById(unsigned int nId)
157 return this->m_root->GetNodeById(nId);
160 int AirwaysTree::GetDepthById(unsigned int nId)
162 return this->m_root->GetDepthById(nId);
165 unsigned int AirwaysTree::CountMarkedNodes(Node* current) const
168 current = this->m_root;
170 bool found = current->FindMarked(aux);
171 if (!found && !current->IsMarked())
173 if (aux != NULL && found)
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);
182 unsigned int AirwaysTree::GetWeight( ) const
184 if (this->m_root == NULL)
187 return this->m_root->GetWeight();
190 unsigned int AirwaysTree::GetHeight(Node* current) const
193 current = this->m_root;
194 if (current->IsLeaf())
196 unsigned int max = 0;
197 vec_nodes::const_iterator it = current->GetChildren().begin();
198 for (; it != current->GetChildren().end(); ++it)
200 unsigned int cMax = this->GetHeight(*it);
204 if (current == this->m_root)
209 unsigned int AirwaysTree::GetLevels(Node* current, unsigned int cLevel) const
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)
219 unsigned int level = this->GetLevels(*it, maxLevel);
220 if (level > maxLevel)
226 unsigned int AirwaysTree::GetNumberOfLeafs( ) const
231 return this->m_root->GetNumberOfLeafs( );
234 const EdgeVector& AirwaysTree::GetEdges() const
236 return this->m_edges;
239 const vec_nodes& AirwaysTree::GetNodes() const
241 return this->m_nodes;
244 const Node* AirwaysTree::GetClosestBrachNodeToIndex(float point_x, float point_y, float point_z) const
247 Node* closestNode = this->m_root->GetClosestBranchNodeToPoint(point_x, point_y, point_z, distance);
251 void AirwaysTree::GetBrancheNodesInfluencingPoint(float point_x, float point_y, float point_z, vec_nodes& vec_nodes_influence)
255 std::cout << "AirwaysTree::GetBrancheNodesInfluencingPoint" << std::endl;
256 this->m_root->GetBrancheNodesInfluencingPoint(point_x, point_y, point_z, vec_nodes_influence);
260 AirwaysTree& AirwaysTree::GetTreeFromMarkedNodes(Node* currentNode, AirwaysTree* tree)
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);
277 if (currentNode->IsMarked())
279 Edge* edge = new Edge();
280 *edge = *currentNode->GetEdge();
281 Node* target = new Node(currentNode->GetCoords());
282 edge->m_target = target;
284 vec_nodes::const_iterator it = currentNode->GetChildren().begin();
285 for (; it != currentNode->GetChildren().end(); ++it)
286 this->GetTreeFromMarkedNodes(*it, tree);
288 return *(new AirwaysTree());
291 AirwaysTree& AirwaysTree::GetSubTreeFromNode(Node* node, AirwaysTree* tree)
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);
305 Edge* edge = new Edge();
306 *edge = *node->GetEdge();
307 Node* target = new Node(node->GetCoords());
308 edge->m_target = target;
310 vec_nodes::const_iterator it = node->GetChildren().begin();
311 for (; it != node->GetChildren().end(); ++it)
312 this->GetSubTreeFromNode(*it, tree);
314 return *(new AirwaysTree());
317 AirwaysTree& AirwaysTree::GetSubTreeByLevels(const int& level, Node* node,
323 tree = new AirwaysTree();
324 Node* current = new Node(node->GetCoords());
325 tree->SetRoot(current);
326 if (node->GetLevel() <= level)
329 vec_nodes::const_iterator it = node->GetChildren().begin();
330 for (; it != node->GetChildren().end(); ++it)
331 this->GetSubTreeByLevels(level, *it, tree);
333 if (node->GetLevel() == level)
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();
345 this->ImageReconstruction(tree->m_img);
346 //this->ImageReconstruction(tree->m_skeleton, true);
351 if (node->GetLevel() <= level)
354 Edge* edge = new Edge();
355 *edge = *node->GetEdge();
356 Node* target = new Node(node->GetCoords());
357 edge->m_target = target;
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)
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)
372 this->m_edges.erase(it2);
375 vec_nodes::iterator it3 = this->m_nodes.begin();
376 for (; it3 != this->m_nodes.end(); ++it3)
379 this->m_nodes.erase(it3);
384 return *(new AirwaysTree());
387 AirwaysTree& AirwaysTree::GetSubTreeByLength(const double& length, Node* node,
388 AirwaysTree* tree, double cLenght)
393 tree = new AirwaysTree();
394 Node* current = new Node(node->GetCoords());
395 tree->SetRoot(current);
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);
405 cLenght += node->GetEdge()->m_length;
406 if (cLenght <= length)
409 Edge* edge = new Edge();
410 *edge = *node->GetEdge();
411 Node* target = new Node(node->GetCoords());
412 edge->m_target = target;
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)
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)
427 this->m_edges.erase(it2);
430 vec_nodes::iterator it3 = this->m_nodes.begin();
431 for (; it3 != this->m_nodes.end(); ++it3)
434 this->m_nodes.erase(it3);
439 return *(new AirwaysTree());
442 AirwaysTree& AirwaysTree::GetSubTreeByRadius(const double& radius, Node* node,
443 AirwaysTree* tree, double cRadius)
448 tree = new AirwaysTree();
449 Node* current = new Node(node->GetCoords());
450 tree->SetRoot(current);
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);
460 cRadius += node->GetEdge()->m_aRadius;
461 if (cRadius >= radius)
464 Edge* edge = new Edge();
465 *edge = *node->GetEdge();
466 Node* target = new Node(node->GetCoords());
467 edge->m_target = target;
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)
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)
482 this->m_edges.erase(it2);
485 vec_nodes::iterator it3 = this->m_nodes.begin();
486 for (; it3 != this->m_nodes.end(); ++it3)
489 this->m_nodes.erase(it3);
494 return *(new AirwaysTree());
497 void AirwaysTree::MarkLevels(const int& level, Node* currentNode)
499 if (currentNode == NULL)
500 currentNode = this->m_root;
501 if (currentNode->GetLevel() <= level)
504 vec_nodes::const_iterator it = currentNode->GetChildren().begin();
505 for (; it != currentNode->GetChildren().end(); it++)
506 this->MarkLevels(level, *it);
510 bool AirwaysTree::IsEmpty() const
512 return this->m_root == NULL;
515 bool AirwaysTree::IsNodeInsidePath(unsigned int idNode_starting, unsigned int idNode_ending, unsigned int idNode_searched)
517 // Check if the starting point exist
518 Node* node_starting = this->m_root->GetNodeById(idNode_starting);
522 // Check if the searched node exist in the subtree of the starting node
523 Node* node_searched = node_starting->GetNodeById(idNode_searched);
527 // Check if the ending node exist in the subtree of the searched node
528 if(!node_searched->HasNodeId(idNode_ending))
531 // If all conditions were true then the node is in the path
535 void AirwaysTree::SetRoot(Node* root)
538 this->m_root->SetLevel(0);
539 this->m_nodes.push_back(this->m_root);
542 void AirwaysTree::SetImageName(const std::string& img_name)
544 this->m_img_name = img_name;
547 void AirwaysTree::SetResultPath(const std::string& folderPath)
549 this->m_result_path = folderPath;
552 bool AirwaysTree::Insert(Edge* edge)
554 if (this->m_root == NULL)
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());
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) )
564 else if (source == NULL && target != NULL) // If target was found
566 // Change the direction of the edge
567 edge->m_target = edge->m_source;
568 edge->m_source = target;
571 edge->m_source = source;
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);
579 // Add the new child node
580 this->m_nodes.push_back(edge->m_target);
581 this->m_edges.push_back(edge);
584 edge->m_target->SetId(this->m_nodes.size());
588 void AirwaysTree::UnMarkAll(Node* currentNode)
590 if (currentNode == NULL)
591 currentNode = this->m_root;
592 currentNode->UnMark();
594 vec_nodes::const_iterator it = currentNode->GetChildren().begin();
595 for (; it != currentNode->GetChildren().end(); ++it)
596 this->UnMarkAll(*it);
599 void AirwaysTree::UpdateLevels(unsigned int levels)
601 this->m_root->UpdateLevels(levels);
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)
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.
619 if( !this->IsEmpty() && !treeB.IsEmpty() )
622 // ----------------------------
623 // ---- Pre-Matching steps ----
625 // Get the root nodes
626 Node* root_A = this->m_root;
627 Node* root_B = treeB.m_root;
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);
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);
643 // ----------------------------
644 // ----- Run the comparison----
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);
650 std::cout << "Trees are empties." << std::endl;
653 void AirwaysTree::MarkPathFromNodeToNode(Node* node_begin, Node* node_end)
656 this->m_root->MarkPathFromNodeToNode(node_begin, node_end);
659 void AirwaysTree::CompareTrees(AirwaysTree& treeB)
661 // Compare all the "branches" (paths) in treeA with all branches in treeB
662 Vector_Pair_PairNodes_Score vectorPairNodesScore = this->CompareAllBranches(treeB);
664 std::cout << "Vector pairs size:" << vectorPairNodesScore.size() << std::endl;
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.
676 //AirwaysTree* tree_copyA = this->GetCopy();
677 //AirwaysTree* tree_copyB = treeB.GetCopy();
679 //NodeVector nonCommonNodes_A;
680 //NodeVector commonNodes_A;
681 //NodeVector nonCommonNodes_B;
682 //NodeVector commonNodes_B;
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;
689 // Until one of the trees is not empty
690 while( this->GetWeight() > 1 && treeB.GetWeight() > 1 )
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;
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;
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())
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)
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);
724 // If the score is greater
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();
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;
736 it=commonNodes_A.find(node_actualA->GetId());
737 if(it==commonNodes_A.end())
739 std::cout << "CommonNode A:" << node_actualA->GetId() << std::endl;
740 commonNodes_A[node_actualA->GetId()] = node_actualA;
743 //commonNodes_B.push_back(node_actualB);
744 //commonNodes_A.push_back(node_actualA);
746 this->DeleteNodeById(node_actualA->GetId());
747 this->DeleteEntriesByIdNode(vectorPairNodesScore, 0, node_actualA->GetId());
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)
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);
764 // If the score is greater
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);
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;
778 it=commonNodes_A.find(node_actualA->GetId());
779 if(it==commonNodes_A.end())
780 commonNodes_A[node_actualA->GetId()] = node_actualA;
782 treeB.DeleteNodeById(node_actualB->GetId());
783 this->DeleteEntriesByIdNode(vectorPairNodesScore, 1, node_actualB->GetId());
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;
791 // Get the branches with the lowest difference
793 double lowestValue = 0;
794 Pair_Node_Node bestPairNodes;
796 for(Vector_Pair_PairNodes_Score::iterator it = vectorPairNodesScore.begin(); it != vectorPairNodesScore.end(); ++it )
800 lowestValue = (*it).second;
801 bestPairNodes = (*it).first;
804 else if( (*it).second < lowestValue )
806 lowestValue = (*it).second;
807 bestPairNodes = (*it).first;
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;
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.
824 if (this->IsEmpty() || treeB.IsEmpty())
826 std::cout << "One or both airways are empty." << std::endl;
830 Edge actualRootEdge = this->m_root->GetEdge();
831 Edge comparedRootEdge = treeB.m_root->GetEdge();
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);
839 Pair_PairNodes_Score AirwaysTree::GetBestLeaf(Vector_Pair_PairNodes_Score vectorPairNodesScore)
842 Pair_PairNodes_Score pairNodesScore_final;
845 for(Vector_Pair_PairNodes_Score::iterator it = vectorPairNodesScore.begin(); it != vectorPairNodesScore.end(); ++it)
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() )
856 pairNodesScore_final = (*it);
858 else if((*it).second < pairNodesScore_final.second)
860 pairNodesScore_final = (*it);
866 std::cout << "Best leaf not found" << std::endl;
868 return pairNodesScore_final;
871 double AirwaysTree::GetPairNodes_Score(Vector_Pair_PairNodes_Score vectorPairNodesScore, const unsigned int idA, const unsigned int idB)
876 for(Vector_Pair_PairNodes_Score::iterator it = vectorPairNodesScore.begin(); it != vectorPairNodesScore.end(); ++it)
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 )
885 score = (*it).second;
890 std::cout << "Pair of Ids not found: " << idA << ", " << idB << std::endl;
895 bool AirwaysTree::DeleteNodeById(unsigned int nId)
897 bool deleted = false;
898 if(m_root->GetId() == nId)
905 return this->m_root->DeleteNodeById(nId);
908 unsigned int AirwaysTree::DeleteEntriesByIdNode(Vector_Pair_PairNodes_Score& vectorPairNodesScore, unsigned int component, unsigned int nId)
910 unsigned int pairs_deleted = 0;
912 Vector_Pair_PairNodes_Score::iterator it = vectorPairNodesScore.begin();
913 while(it != vectorPairNodesScore.end())
915 // Check that one of the node is a leaf
916 pair_nodes pair_actual = (*it).first;
917 unsigned int id_actual = -1;
919 id_actual = pair_actual.first->GetId();
921 id_actual = pair_actual.second->GetId();
923 if( id_actual == nId)
925 it = vectorPairNodesScore.erase(it);
932 std::cout << "Pairs deleted: " << pairs_deleted << std::endl;
934 return pairs_deleted;
937 AirwaysTree* AirwaysTree::GetCopy()
942 Vector_Pair_PairNodes_Score AirwaysTree::CompareAllBranches(AirwaysTree& treeB)
944 std::cout << "Compare all branches ..." << std::endl;
946 int nodesA = this->m_nodes.size();
947 int nodesB = treeB.m_nodes.size();
949 Vector_Pair_PairNodes_Score vectorPairNodesScore;
951 // For all the actual branches ( nodes different from the root)
952 for(unsigned int idA = 2; idA <= nodesA; ++idA)
954 // Get the actual branch in A
955 Edge* branchA = this->GetComposedEdge(idA);
957 //Get the final actual node in A
958 Node* node_actualA = this->GetNodeById(idA);
960 std::cout << "IdA: " << idA << std::endl;
962 // For all the branches in B tree
963 for(unsigned int idB = 2; idB <= nodesB; ++idB)
965 // Get the actual branch in B
966 Edge* branchB = treeB.GetComposedEdge(idB);
968 //Get the final actual node in B
969 Node* node_actualB = treeB.GetNodeById(idB);
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);
975 // Get the final score
978 //double score = scoreB2A > scoreA2B ? scoreB2A : scoreA2B ;
981 double score = scoreB2A + scoreA2B;
983 // Build the pair of nodes
984 pair_nodes actualPairNodes(node_actualA, node_actualB);
986 Pair_PairNodes_Score actualPairNodesScore(actualPairNodes, score);
987 vectorPairNodesScore.push_back(actualPairNodesScore);
989 //std::cout << "IdA: " << idA << ", IdB: " << idB << ", score: " << score << ", ... OK" << std::endl;
995 std::cout << "Compare all branches ... OK" << std::endl;
996 return vectorPairNodesScore;
999 Edge* AirwaysTree::GetComposedEdge(unsigned int idNode)
1001 return this->m_root->GetComposedEdge(idNode);
1004 AirwaysTree* AirwaysTree::GetSingleDetachedTreeNodeById(unsigned int nId)
1006 AirwaysTree* tree_detached = NULL;
1007 Node* node_searched = this->GetNodeById(nId);
1010 Node* node_root_detached = node_searched->GetDetachedRootCopy();
1011 tree_detached = new AirwaysTree();
1012 tree_detached->m_root = node_root_detached;
1015 std::cout << "Node not found" << std::endl;
1016 return tree_detached;
1019 void AirwaysTree::SubIsomorphism(AirwaysTree& treeB)
1021 if (this->IsEmpty() || treeB.IsEmpty())
1023 //Marking the root node
1024 this->m_root->Mark();
1025 treeB.m_root->Mark();
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)
1031 unsigned int totalCost = 0;
1032 Node* currentB = NULL;
1033 Node* currentA = *it;
1034 if (currentA->IsMarked())
1036 vec_nodes::const_iterator it2 = treeB.m_root->GetChildren().begin();
1037 for (; it2 != treeB.m_root->GetChildren().end(); ++it2)
1039 if ((*it2)->IsMarked())
1041 unsigned int currentCost = (*it)->GetCost(*it2);
1042 if (totalCost < currentCost)
1045 totalCost = currentCost;
1048 if (currentB != NULL)
1049 if (*currentA == *currentB)
1050 this->GetSubTreeByDetachingNode(currentA).SubIsomorphism(
1051 treeB.GetSubTreeByDetachingNode(currentB));
1055 it = this->m_root->GetChildren().begin();
1056 for (; it != this->m_root->GetChildren().end(); ++it)
1058 if ((*it)->IsMarked())
1060 vec_nodes::const_iterator it2 = treeB.m_root->GetChildren().begin();
1061 for (; it2 != treeB.m_root->GetChildren().end(); ++it2)
1063 Node* currentB = *it2;
1064 Node* currentA = *it;
1065 if (*currentB > *currentA)
1067 if (currentA->FindSimilarEdgeByAccumulation(currentB))
1069 AirwaysTree subTreeA = this->GetSubTreeByDetachingNode(
1071 AirwaysTree subTreeB = treeB.GetSubTreeByDetachingNode(
1073 subTreeB.UpdateLevels((*it2)->GetLevel());
1074 subTreeA.SubIsomorphism(subTreeB);
1075 treeB.m_root->UpdateLevels(treeB.m_root->GetLevel());
1077 if (!(*it2)->MarkPath(currentB))
1078 std::cout << "Error finding path" << std::endl;
1083 vec_nodes::const_iterator it2 = treeB.m_root->GetChildren().begin();
1084 for (; it2 != treeB.m_root->GetChildren().end(); ++it2)
1086 if ((*it2)->IsMarked())
1088 it = this->m_root->GetChildren().begin();
1089 for (; it != this->m_root->GetChildren().end(); ++it)
1091 Node* currentB = *it2;
1092 Node* currentA = *it;
1093 if (*currentA > *currentB)
1095 if (currentB->FindSimilarEdgeByAccumulation(currentA))
1097 AirwaysTree subTreeA = this->GetSubTreeByDetachingNode(
1099 AirwaysTree subTreeB = treeB.GetSubTreeByDetachingNode(
1101 subTreeA.UpdateLevels((*it)->GetLevel());
1102 subTreeA.SubIsomorphism(subTreeB);
1103 this->m_root->UpdateLevels(this->m_root->GetLevel());
1105 if (!(*it)->MarkPath(currentA))
1106 std::cout << "Error finding path" << std::endl;
1111 //End of second step
1114 bool AirwaysTree::Isomorphism(AirwaysTree& tree)
1116 if ((this->GetWeight() != tree.GetWeight())
1117 || (this->GetHeight() != tree.GetHeight()))
1119 return this->m_root->CompareNodes(tree.m_root);
1122 void AirwaysTree::ImageReconstructionFromNode(TInputImage::Pointer& result, Node* node, bool skeleton)
1126 //Finding the center of the sphere
1127 const Edge* edge = node->GetEdge();
1130 //making a copy of the current image
1131 TInputImage::Pointer copy;
1132 DuplicatorType::Pointer duplicator = DuplicatorType::New();
1134 duplicator->SetInputImage(this->m_skeleton);
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)
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)
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;
1159 this->m_skeleton->TransformPhysicalPointToIndex(point, indCenter);
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);
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;
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;
1186 TInputImage::RegionType region;
1187 region.SetSize(regionSize);
1188 region.SetIndex(regionIndex);
1189 RegionIterator it2(copy, region);
1190 while (!it2.IsAtEnd())
1192 TInputImage::IndexType ind = it2.GetIndex();
1193 SphereType::InputType point;
1197 SphereType::OutputType output = sphere->Evaluate(point);
1204 * First region growing - Segmenting the wanted branch to be removed
1205 * Removing the branch from the image
1207 if (!ComputeSimpleRegionGrowing(copy, this->m_root->GetCoords()))
1209 SubtractImageFilterType::Pointer subtractFilter =
1210 SubtractImageFilterType::New();
1212 subtractFilter->SetInput1(this->m_skeleton);
1214 subtractFilter->SetInput1(this->m_img);
1215 subtractFilter->SetInput2(copy);
1216 subtractFilter->Update();
1217 TInputImage::Pointer diff = subtractFilter->GetOutput();
1219 * Second region growing - Cleaning the resulting image
1221 Node* leaf = this->GetLeaf(node);
1222 if (ComputeSimpleRegionGrowing(diff, leaf->GetCoords()))
1225 std::cout << "problemas" << std::endl;
1226 //cin.ignore().get(); //Pause Command for Linux Terminal
1230 void AirwaysTree::ImageReconstruction(TInputImage::Pointer& result, bool skeleton, Node* node)
1234 DuplicatorType::Pointer duplicator = DuplicatorType::New();
1236 duplicator->SetInputImage(this->m_skeleton);
1238 duplicator->SetInputImage(this->m_img);
1239 duplicator->Update();
1240 result = duplicator->GetOutput();
1241 node = this->m_root;
1243 if (!node->IsMarked())
1245 this->RemoveBranchFromImage(result, node);
1248 vec_nodes::const_iterator it = node->GetChildren().begin();
1249 for (; it != node->GetChildren().end(); ++it)
1250 this->ImageReconstruction(result, skeleton, *it);
1254 void AirwaysTree::ImageReconstructionBySpheres(TInputImage::Pointer& result, Node* node, bool useMarks)
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())
1272 if (node->IsMarked() || !useMarks)
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);
1281 Node* AirwaysTree::GetNode(const Vec3 nodeCoords, Node* current)
1283 if (current == NULL)
1284 current = this->m_root;
1285 if (current->GetCoords() == nodeCoords)
1287 vec_nodes::const_iterator it = current->GetChildren().begin();
1288 for (; it != current->GetChildren().end(); ++it)
1290 Node* child = GetNode(nodeCoords, *it);
1297 Node* AirwaysTree::GetLeaf(Node* current)
1299 if (current == NULL)
1300 current = this->m_root;
1301 if (current->IsLeaf())
1303 vec_nodes::const_iterator it = current->GetChildren().begin();
1305 unsigned int maxHeight = 0;
1306 for (; it != current->GetChildren().end(); ++it)
1308 unsigned int cHeight = this->GetHeight(*it);
1309 if (maxHeight <= cHeight)
1311 maxHeight = cHeight;
1316 return this->GetLeaf(node);
1320 void AirwaysTree::UpdateEdges(Node* node)
1323 node = this->m_root;
1325 if (node->GetEdge() != NULL)
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();
1332 vec_nodes::const_iterator it = node->GetChildren().begin();
1333 for (; it != node->GetChildren().end(); ++it)
1334 this->UpdateEdges(*it);
1337 AirwaysTree& AirwaysTree::GetSubTreeByDetachingNode(Node* node)
1339 AirwaysTree* tree = new AirwaysTree();
1340 tree->m_root = node;
1344 DiGraph_EdgePair AirwaysTree::AddBoostEdge(const pair_posVox_rad& source, const pair_posVox_rad& target, DiGraph& graph)
1346 bool vFound[2] = { false, false };
1347 DiGraph_VertexDescriptor vDescriptor[2];
1349 //Remember that DiGraph_VertexIteratorPair, first = begin, second = end
1350 DiGraph_VertexIndexMap index = get(boost::vertex_index, graph);
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)
1355 // Get the current information [vector, intensity]
1356 pair_posVox_rad current = graph[index[*vPair.first]];
1358 // Check if current is equal to the source or target vertex
1359 if (current.first == source.first)
1361 vDescriptor[0] = index[*vPair.first];
1364 else if (current.first == target.first)
1366 vDescriptor[1] = index[*vPair.first];
1371 // If one of the vertices is not in the graph, adds the new vertex
1373 // Add the vertices if they are not in the graph
1375 vDescriptor[0] = boost::add_vertex(source, graph);
1377 vDescriptor[1] = boost::add_vertex(target, graph);
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);
1384 // If the edge does not exist, then add it
1385 if (!edgeC1.second && !edgeC2.second)
1387 DiGraph_EdgePair nEdge = boost::add_edge(vDescriptor[0], vDescriptor[1], graph);
1391 // If an edge has not been added then -> second = false
1392 // Remove the added vertices because the edge already exists.
1394 boost::remove_vertex(vDescriptor[0], graph);
1397 boost::remove_vertex(vDescriptor[1], graph);
1399 edgeC1.second = false;
1400 edgeC1.second = false;
1405 bool AirwaysTree::FindNeighborhood(DiGraph& graph, Vec3Queue& queue, Vec3List& used, const Vec3& end, ConstNeighborhoodIterator& iterator)
1410 // Get the next pixel in the queue
1411 Vec3 current = queue.front();
1413 // Put the pixel in the used list
1414 used.push_back(current);
1416 // Delete next pixel from the queu
1419 //AMP//iterator.GoToBegin();
1420 //AMP//while (!iterator.IsAtEnd())
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)
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];
1444 // Place the iterator in the actual index
1445 iterator.SetLocation(indexActual);
1447 // Create the pair [vector, intensity] for the actual index
1448 pair_posVox_rad src(current, iterator.GetPixel(iterator.GetCenterNeighborhoodIndex()));
1452 // Iterate over the neighbors
1453 for (unsigned int i = 0; i < iterator.Size(); i++)
1455 //AMP//if (iterator.GetPixel(i) != 0 && i != centerIndex)
1457 // If the intensity is not zero and it's the center
1458 if (iterator.GetPixel(i) != 0 && i != iterator.GetCenterNeighborhoodIndex())
1460 // Get the index of the iterator
1461 TInputImage::IndexType ind = iterator.GetIndex(i);
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));
1466 // Add the edge to the graph
1467 bool added = AddBoostEdge(src, target, graph).second;
1469 // If it could not be added then it is because it already exist, then go to next neighboor
1473 // The edge was added. Stop the adding process if the target vertex was reached.
1474 if (target.first == end)
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);
1488 DiGraph AirwaysTree::GetGraphShortestPath(const DiGraph& graph, const Vec3& begin, const Vec3& end)
1491 DiGraph ret; // Returning graph with the shortest path
1492 DiGraph_EdgeIterator ei, eo; // Initial and final iterators
1494 // Iterate over the edges of the graph
1495 boost::tie(ei, eo) = boost::edges(graph);
1496 for (; ei != eo; ++ei)
1498 // Get the source of the actual edge
1499 pair_posVox_rad src = graph[boost::source(*ei, graph)];
1501 // If the source corresponds to the begin pixel
1502 if (src.first == begin)
1504 // Get the target of the actual edge
1505 pair_posVox_rad target = graph[boost::target(*ei, graph)];
1507 // If the actual target corresponds to the end pixel
1508 if (target.first == end)
1510 // Add the edge to the final path
1511 AddBoostEdge(src, target, ret);
1514 std::cout << "Add final target [src]-[ target]:[" << src.first << "]-[" << target.first << "], numVertices:" << boost::num_vertices(ret) << std::endl;
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);
1525 // If the path is empty then the iteration continues
1526 if (boost::num_vertices(ret) == 0)
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);
1533 //std::cout << "Add found edge [src]-[ target]:[" << src.first << "]-[" << target.first << "], numVertices:" << boost::num_vertices(ret) << std::endl;
1542 vec_pair_posVox_rad AirwaysTree::GetGraphShortestPath_AMP(const DiGraph& graph, const Vec3& begin, const Vec3& end)
1544 //std::cout << "SKPairInfoVector AirwaysTree::GetGraphShortestPath_AMP(const DiGraph& graph, const Vec3& begin, const Vec3& end)" << std::endl;
1546 vec_pair_posVox_rad path;
1547 DiGraph ret; // Returning graph with the shortest path
1548 DiGraph_EdgeIterator ei, eo; // Initial and final iterators
1550 // Iterate over the edges of the graph
1551 boost::tie(ei, eo) = boost::edges(graph);
1552 for (; ei != eo; ++ei)
1554 // Get the source of the actual edge
1555 pair_posVox_rad src = graph[boost::source(*ei, graph)];
1557 // If the source corresponds to the begin pixel
1558 if (src.first == begin)
1560 // Get the target of the actual edge
1561 pair_posVox_rad target = graph[boost::target(*ei, graph)];
1563 // If the actual target corresponds to the end pixel
1564 if (target.first == end)
1566 // Add the edge to the final path
1567 AddBoostEdge(src, target, ret);
1568 path.push_back(target);
1569 path.push_back(src);
1572 //std::cout << "Add final target [src]-[ target]:[" << src.first << "]-[" << target.first << "], numVertices:" << boost::num_vertices(ret) << ", path_size:" << path.size() << std::endl;
1576 //std::cout << "SKPairInfoVector AirwaysTree::GetGraphShortestPath_AMP(const DiGraph& graph, const Vec3& begin, const Vec3& end) ... FIRST ... OK" << std::endl;
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;
1586 // If the path is empty then the iteration continues
1587 if (path.size() == 0)
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);
1595 //std::cout << "Add found edge [src]-[ target]:[" << src.first << "]-[" << target.first << "], numVertices:" << boost::num_vertices(ret) << ", path_size:" << path.size() << std::endl;
1597 //std::cout << "SKPairInfoVector AirwaysTree::GetGraphShortestPath_AMP(const DiGraph& graph, const Vec3& begin, const Vec3& end) ... SECOND ... OK" << std::endl;
1603 //std::cout << "SKPairInfoVector AirwaysTree::GetGraphShortestPath_AMP(const DiGraph& graph, const Vec3& begin, const Vec3& end) ... END ... OK" << std::endl;
1607 vec_pair_posVox_rad AirwaysTree::GetSkeletonInfoFromEdge(const Vec3& source, const Vec3& target)
1609 //converting to pixel
1610 TInputImage::PointType pointSource;
1611 pointSource[0] = source[0];
1612 pointSource[1] = source[1];
1613 pointSource[2] = source[2];
1615 TInputImage::PointType pointTarget;
1616 pointTarget[0] = target[0];
1617 pointTarget[1] = target[1];
1618 pointTarget[2] = target[2];
1620 TInputImage::IndexType indexSource;
1621 TInputImage::IndexType indexTarget;
1623 this->m_skeleton->TransformPhysicalPointToIndex(pointSource, indexSource);
1624 this->m_skeleton->TransformPhysicalPointToIndex(pointTarget, indexTarget);
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;
1634 Vec3 sourcePxl(indexSource[0], indexSource[1], indexSource[2]);
1635 Vec3 targetPxl(indexTarget[0], indexTarget[1], indexTarget[2]);
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;
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];
1653 index[0] -= index[0] <= 10 ? 0 : 10;
1654 index[1] -= index[1] <= 10 ? 0 : 10;
1655 index[2] -= index[2] <= 10 ? 0 : 10;
1659 TInputImage::SizeType size;
1660 size[0] = ceil(dstX);
1661 size[1] = ceil(dstY);
1662 size[2] = ceil(dstZ);
1664 TInputImage::RegionType region;
1665 region.SetSize(size);
1666 region.SetIndex(index);
1668 //ConstNeighborhoodIterator iterator(radius, this->m_skeleton, region);
1672 TInputImage::SizeType radius;
1676 ConstNeighborhoodIterator iterator(radius, this->m_skeleton, this->m_skeleton->GetLargestPossibleRegion());
1682 // Set the first pixel as starting pixel
1683 queue.push(sourcePxl);
1685 //FindNeighborhood(graph, queue, used, endPxl, iterator);
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);
1693 std::cout << "***************** END NOT FOUND!!!!" << std::endl;
1695 std::cout << "***************** END FOUND!!!! - Graph Vertices:" << boost::num_vertices(graph) << ", Graph Edges:" << boost::num_edges(graph) << std::endl;
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;
1704 bool beginPixelFound = false;
1705 bool endPixelFound = false;
1706 vec_pair_posVox_rad vector;
1708 for(vec_pair_posVox_rad::reverse_iterator it_path = path.rbegin(); it_path != path.rend(); ++it_path)
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 << "], ";
1716 // To avoid precision problems
1717 if (indexSource == ind)
1719 beginPixelFound = true;
1720 (*it_path).first = source;
1722 else if (indexTarget == ind)
1724 endPixelFound = true;
1725 (*it_path).first = target;
1729 TInputImage::PointType pnt;
1730 this->m_skeleton->TransformIndexToPhysicalPoint(ind, pnt);
1731 (*it_path).first = Vec3(pnt[0], pnt[1], pnt[2]);
1733 vector.push_back((*it_path));
1738 // Variables to iterate over the shortest path and save the results
1740 DiGraph_VertexIterator i, e;
1741 SKPairInfoVector vector;
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;
1747 //for (; e != i; --e)
1749 // Get the actual vertex
1750 SKPairInfo vertex = result[*i];
1751 //SKPairInfo vertex = result[*e];
1752 std::cout << "["<< vertex.first << "], ";
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];
1760 // To avoid precision problems
1761 if (indexBegin == ind)
1763 beginPixelFound = true;
1764 vertex.first = begin;
1766 else if (indexEnd == ind)
1768 endPixelFound = true;
1773 TInputImage::PointType pnt;
1774 this->m_skeleton->TransformIndexToPhysicalPoint(ind, pnt);
1775 vertex.first = Vec3(pnt[0], pnt[1], pnt[2]);
1777 vector.push_back(vertex);
1779 std::cout << "Inserting ... OK" << std::endl;
1785 if(!beginPixelFound)
1786 std::cout << "Begin not found: "<< source << std::endl;
1788 std::cout << "End not found:" << target << std::endl;
1793 bool AirwaysTree::ComputeSimpleRegionGrowing(TInputImage::Pointer& img,
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)
1804 //std::cout << "seed is 0" << std::endl;
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();
1819 void AirwaysTree::RemoveBranchFromImage(TInputImage::Pointer& img, Node* node)
1821 //typedef itk::ImageFileWriter<TInputImage> WriterType;
1822 //Finding the center of the sphere
1823 const Edge* edge = node->GetEdge();
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)
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)
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();
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];
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;
1896 regionSize[0] = size[0] - regionIndex[0];
1898 regionSize[1] = size[1] - regionIndex[1];
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())
1910 TInputImage::IndexType ind = it2.GetIndex();
1911 SphereType::InputType point;
1915 SphereType::OutputType output = sphere->Evaluate(point);
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)
1930 std::cout << "It entered here!!!!" << std::endl;
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();
1944 /*WriterType::Pointer writer = WriterType::New();
1945 writer->SetInput(copy);
1946 writer->SetFileName("output_antes.mhd");
1949 //std::cout << "salio" << std::endl;
1951 * First region growing - Segmenting the wanted branch to be removed
1952 * Removing the branch from the image
1954 Node* leaf = this->GetLeaf(node);
1955 if (!ComputeSimpleRegionGrowing(copy, leaf->GetCoords()))
1957 SubtractImageFilterType::Pointer subtractFilter =
1958 SubtractImageFilterType::New();
1959 subtractFilter->SetInput1(img);
1960 subtractFilter->SetInput2(copy);
1961 subtractFilter->Update();
1962 TInputImage::Pointer diff = subtractFilter->GetOutput();
1964 * Second region growing - Cleaning the resulting image
1966 if (ComputeSimpleRegionGrowing(diff, this->m_root->GetCoords()))
1969 std::cout << "problemas" << std::endl;*/
1970 //cin.ignore().get(); //Pause Command for Linux Terminal
1973 //Metodo obsoleto, ya no nos interesa la reconstruccion por esferas
1974 void AirwaysTree::CreateSpheres(TInputImage::Pointer& img, Node* node)
1976 if (node->GetEdge() == NULL)
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)
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);
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())
2012 TInputImage::IndexType ind = it2.GetIndex();
2013 SphereType::InputType point;
2017 SphereType::OutputType output = sphere->Evaluate(point);
2025 void AirwaysTree::printUpToLevel(int level)
2027 cout << "Printing up to level: " << level <<endl;
2028 this->m_root->printUpToLevel(level);
2029 cout << "Printing done" <<endl;
2032 void AirwaysTree::printNodeAndChildrenIds()
2034 cout << "Printing nodes and children: " <<endl;
2035 this->m_root->printNodeAndChildrenIds( );
2036 cout << "Printing nodes and children DONE" <<endl;
2039 void AirwaysTree::createQuaternionsUpToLevel(int level)
2041 cout << "Creating quaternions up to level: " << level <<endl;
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;
2053 this->m_root->createQuaternionsUpToLevel(level, level);
2056 cout << "Quaternions must be created from level 2." <<endl;
2057 cout << "Creating quaternions done" <<endl;
2060 void AirwaysTree::getStatisticsBifurcations()
2062 cout << "Statistics bifurcations ... " <<endl;
2064 int p[6] = {0,0,0,0,0,0};
2066 this->m_root->getStasticsBifurcations(p);
2068 cout << "Bifurcations:" << endl;
2069 for(int i = 0; i < 6; i++)
2071 cout << i+1 << "->" << p[i] << endl;
2074 cout << "Statistics bifurcations done " <<endl;
2077 void AirwaysTree::saveAirwayToImage(std::string filename, int dims[3], double spc[3], double nOrigin[3])
2079 TInputImage::Pointer imagePointer;
2080 imagePointer = TInputImage::New();
2082 TInputImage::IndexType start;
2087 TInputImage::SizeType size;
2092 TInputImage::RegionType region;
2093 region.SetSize( size );
2094 region.SetIndex( start );
2096 TInputImage::SpacingType spacing;
2097 spacing[0] = spc[0];
2098 spacing[1] = spc[1];
2099 spacing[2] = spc[2];
2101 TInputImage::PointType origin;
2102 origin[0] = nOrigin[0];
2103 origin[1] = nOrigin[1];
2104 origin[2] = nOrigin[2];
2106 imagePointer->SetRegions( region );
2107 imagePointer->SetSpacing( spacing );
2108 imagePointer->SetOrigin( origin );
2109 imagePointer->Allocate();
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)
2116 Vec3 coords_father = (*it_children)->GetEdge()->GetSource()->GetCoords();
2117 Vec3 coords_child = (*it_children)->GetCoords();
2118 drawLineInImage(coords_father, coords_child, imagePointer);
2121 cout << "Writing ..." << endl;
2122 typedef itk::ImageFileWriter<TInputImage> WriterType;
2123 WriterType::Pointer writer = WriterType::New();
2124 writer->SetInput(imagePointer);
2125 writer->SetFileName(filename);
2127 cout << "Writing ... done" << endl;
2130 void AirwaysTree::drawLineInImage(Vec3 initVoxel, Vec3 endVoxel, TInputImage::Pointer imagePointer)
2132 TInputImage::RegionType region = imagePointer->GetLargestPossibleRegion();
2134 TInputImage::SpacingType spacing = imagePointer->GetSpacing();
2136 TInputImage::PointType origin = imagePointer->GetOrigin();
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];
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];
2151 itk::LineIterator<TInputImage> it(imagePointer, init, end);
2153 while (!it.IsAtEnd())
2160 //To store the graph in a vtk
2161 void AirwaysTree::saveAirwayAsVTK(const std::string filename, bool common)
2163 vtkSmartPointer<vtkUnsignedCharArray> colors = vtkSmartPointer<vtkUnsignedCharArray>::New();
2164 colors->SetNumberOfComponents(3);
2165 colors->SetName("Colors");
2167 //pointer of lines and points
2168 vtkSmartPointer<vtkCellArray> lines = vtkSmartPointer<vtkCellArray>::New();
2169 vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
2171 Vec3 root = this->m_root->GetCoords();
2176 unsigned int id = 1;
2178 CalculateVTKLinesFromEdges(this->m_root, 0, id, pts, lines, colors, common);
2180 vtkSmartPointer<vtkPolyData> linesPolyData = vtkSmartPointer<vtkPolyData>::New();
2181 linesPolyData->SetPoints(pts);
2182 linesPolyData->SetLines(lines);
2184 linesPolyData->GetCellData()->SetScalars(colors);
2187 vtkSmartPointer<vtkPolyDataWriter> writer = vtkSmartPointer<vtkPolyDataWriter>::New();
2188 writer->SetFileName(filename.c_str());
2189 #if VTK_MAJOR_VERSION <= 5
2190 writer->SetInput(linesPolyData);
2192 writer->SetInputData(linesPolyData);
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)
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 };
2212 const vec_nodes children = node->GetChildren();
2214 for(vec_nodes::const_iterator it = children.begin(); it != children.end(); ++it)
2216 if (!(*it)->IsMarked() && common)
2219 const Edge* edge = (*it)->GetEdge();
2220 //pts->InsertNextPoint(edge->GetSource()->GetCoords().GetVec3());
2221 //pts->InsertNextPoint(edge->GetTarget()->GetCoords().GetVec3());
2223 pts->InsertNextPoint(node->GetCoords().GetVec3());
2224 pts->InsertNextPoint((*it)->GetCoords().GetVec3());
2226 // Set color to be used
2227 int numColor = node->GetLevel() % 4;
2229 colors->InsertNextTupleValue(green);
2230 else if(numColor == 1)
2231 colors->InsertNextTupleValue(red);
2232 else if(numColor == 2)
2233 colors->InsertNextTupleValue(blue);
2235 colors->InsertNextTupleValue(yellow);
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);
2247 } /* namespace airways */