appli/TempAirwaysAppli/TempAirwaysAppli.cxx

   1 #include <iostream>
   2 #include <cstdlib>
   3 #include <sstream>
   4
   5 #include <boost/filesystem.hpp>
   6 #include <boost/program_options.hpp>
   7
   8 #include "vec3.h"
   9 #include "airwaysTree.h"
  10
  11 #include <cpPlugins/Interface.h>
  12 #include <cpPlugins/Workspace.h>
  13
  14 using namespace airways;
  15 namespace po = boost::program_options;
  16
  17 namespace
  18 {
  19   const size_t ERROR_IN_COMMAND_LINE = 1;
  20   const size_t SUCCESS = 0;
  21   const size_t ERROR_UNHANDLED_EXCEPTION = 2;
  22 } // namespace
  23
  24 typedef std::vector< AirwaysTree > AirwaysVector;
  25
  26 #include <itkCastImageFilter.h>
  27
  28 // Auxiliar struct to save info for execution.
  29 typedef void* TImagePointer;
  30 struct TreeInfo{
  31
  32   TreeInfo( )
  33     {
  34       this->myWorkspace = new cpPlugins::Workspace( );
  35     }
  36   ~TreeInfo( )
  37     {
  38       /*
  39         if( this->IsMyWorkspace )
  40         delete this->myWorkspace;
  41       */
  42     }
  43
  44   void CastImage( )
  45     {
  46       auto image = this->myWorkspace->GetFilter( "reader" )->GetOutputData( "Output" )->GetITK< itk::Image< unsigned char, 3 > >( );
  47       typedef itk::CastImageFilter< itk::Image< unsigned char, 3 >, TInputImage > _TCast;
  48       _TCast::Pointer cast = _TCast::New( );
  49       cast->SetInput( image );
  50       cast->Update( );
  51       this->Image = cast->GetOutput( );
  52       this->Image->DisconnectPipeline( );
  53     }
  54
  55   TInputImage::Pointer Image;
  56   cpPlugins::Workspace* myWorkspace;
  57   bool IsMyWorkspace;
  58   Vec3 seed;
  59   std::string folderpath_pigResults;
  60   std::string pig_name;
  61   std::string image_name;
  62
  63 };
  64 // Constants
  65 unsigned char red[3]   = { 255, 0, 0 };
  66 unsigned char green[3] = { 0, 255, 0 };
  67 unsigned char blue[3]  = { 0, 0, 255 };
  68 unsigned char yellow[3]= { 255, 255, 0 };
  69 unsigned char purple[3]= { 255, 0, 255 };
  70 unsigned char cyan[3]= { 0, 255, 255 };
  71
  72 cpPlugins::Interface myPlugins;
  73
  74 // -------------------------------------------------------------------------
  75 void Load_cpPlugins( const std::string& plugins );
  76 void CreateResultDirectory(AirwaysTree tree);
  77 void DrawVTKLinesFromTree(AirwaysTree& tree, const std::string filename, bool common);
  78 AirwaysTree& CreateAirwaysTreeFromSegmentation(Vec3 seed, TInputImage* input_image, cpPlugins::Workspace& ws );
  79 vector<TreeInfo> ReadInputFile(const char* filename);
  80 void printCommonTreeBetweenTwoTrees(AirwaysTree tree_A, AirwaysTree tree_B, std::vector< std::pair< std::pair<Node*, Node*>, std::pair<Node*, Node*> > > vector_pair_edges_A_to_B, unsigned int Q, unsigned int F);
  81 void printMatchingResultToFile(AirwaysTree tree_A, AirwaysTree tree_B, std::map< unsigned int, std::vector<Node*> > map_A_to_B, std::map< unsigned int, std::vector<Node*> > map_B_to_A, unsigned int Q, unsigned int F);
  82 void createLinesAndPointsForVTK(const Node* node, vtkSmartPointer<vtkPoints>& pts,vtkSmartPointer<vtkCellArray>& lines, vtkSmartPointer<vtkUnsignedCharArray>& colors, bool common, bool isRoot);
  83
  84 // -------------------------------------------------------------------------
  85 int main( int argc, char* argv[] )
  86 {
  87   Load_cpPlugins( "./plugins.cfg" );
  88
  89   try
  90   {
  91     // Define and parse the program options
  92     po::options_description desc("Options");
  93     desc.add_options()("use_file", po::value<std::string>(), "Adds the filepath containing the file to be used in creaAirwaysTree -- Mandatory")
  94       ("build_trees","Creates trees from a given filepath (arg) and stores it in a .vtk file")
  95       ("subtree_levels", po::value<int>(),"Get a subtree(s) by levels - result: a vtk and reconstructed img .mhd")
  96       ("subtree_length", po::value<float>(),"Get a subtree(s) by length - result: a vtk and reconstructed img .mhd")
  97       ("subtree_diameter", po::value<float>(),"Get a subtree(s) by diameter - result: a vtk and reconstructed img .mhd")
  98       ("compare_trees", "Compare the trees given in the input file or the subtrees using an option - result: a vtk and reconstructed img .mhd");
  99     //TODO: Fix image segmentation. ("segment_images",  "Creates a binary image corresponding to the segmentation of of a given original image (arg) and seed (arg), saves it to a .mhd file and uses it for subsequent operations")
 100     po::variables_map vm;
 101     try
 102     {
 103       std::string fileName;
 104       vector<TreeInfo> treeInfoVector;
 105       AirwaysVector aVector;
 106       TImagePointer segmentationImage;
 107       Vec3 seed;
 108       po::store(po::parse_command_line(argc, argv, desc), vm); // can throw
 109       if (vm.count("use_file"))
 110       {
 111         fileName = vm["use_file"].as<std::string>();
 112         treeInfoVector = ReadInputFile(fileName.c_str());
 113       }
 114
 115       if(vm.count("segment_images"))
 116       {
 117       }
 118
 119       // Build trees option
 120       if (vm.count("build_trees"))
 121       {
 122         for(unsigned int i = 0; i < treeInfoVector.size(); ++ i){
 123           TreeInfo info = treeInfoVector[i];
 124           std::string err = info.myWorkspace->Execute( "eb" );
 125           if( err != "" )
 126           {
 127             std::cerr << "Error: " << err << std::endl;
 128             std::exit( 1 );
 129
 130           } // fi
 131           info.CastImage( );
 132           AirwaysTree tree = CreateAirwaysTreeFromSegmentation( info.seed, info.Image, *(info.myWorkspace ));
 133           tree.SetResultPath(info.folderpath_pigResults);
 134           tree.SetImageName(info.image_name);
 135           aVector.push_back(tree);
 136         }
 137         for (unsigned int i = 0; i < aVector.size(); ++i)
 138         {
 139           CreateResultDirectory(aVector[i]);
 140           std::string fullPath = aVector[i].GetResultPath() + "/" + aVector[i].GetImageName() + ".vtk";
 141           DrawVTKLinesFromTree(aVector[i], fullPath, false);
 142
 143         } //for
 144       } //if
 145
 146       // Subtree levels option
 147       if (vm.count("subtree_levels"))
 148       {
 149       } //if
 150       else if (vm.count("subtree_length"))
 151       {
 152       } //if
 153       else if (vm.count("subtree_diameter"))
 154       {
 155       } //if
 156
 157       if (vm.count("compare_trees"))
 158       {
 159         std::cout << "Option: Compare trees" << std::endl;
 160         /**
 161          * this piece of code only test in the case of two trees, so it can be replaced
 162          * to a code which does a cascade
 163          */
 164
 165         //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 166         //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 167         //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 168         //CODIGO DE COMPARACIÓN DE ARBOLES
 169         //std::cout << "Comparing trees ..." << std::endl;
 170         //aVector[0].CompareTrees(aVector[1]);
 171         //std::cout << "Comparing trees ... OK" << std::endl;
 172
 173         std::cout << "Tree A ... " << std::endl;
 174         //aVector[0].printNodeAndChildrenIds();
 175         std::cout << "Tree A ... OK" << std::endl;
 176
 177         std::cout << "Tree B ... " << std::endl;
 178         //aVector[1].printNodeAndChildrenIds();
 179         std::cout << "Tree B ... OK" << std::endl;
 180
 181         std::cout << "Comparing trees Orkisz-Morales..." << std::endl;
 182         // Vectors to save the common and uncommon nodes
 183         vec_nodes commonA;
 184         vec_nodes commonB;
 185         vec_nodes nonCommonA;
 186         vec_nodes nonCommonB;
 187         std::map< unsigned int, std::vector<Node*> > map_A_to_B;
 188         std::map< unsigned int, std::vector<Node*> > map_B_to_A;
 189         std::vector< std::pair< std::pair<Node*, Node*>, std::pair<Node*, Node*> > > vector_pair_edges_A_to_B;
 190
 191         // Input parameter for comparison
 192         unsigned int Q = 1; // Corresponds to the depth to select "fathers" nodes
 193         unsigned int F = 1; // Correspond to the depth to select "family" nodes
 194         std::cout << "Q = " << Q << std::endl;
 195         std::cout << "F = " << F << std::endl;
 196
 197         clock_t before_compare = clock();
 198         aVector[0].CompareTreesOrkiszMorales(aVector[1], Q, F, commonA, commonB, nonCommonA, nonCommonB, map_A_to_B, map_B_to_A, vector_pair_edges_A_to_B);
 199         clock_t after_compare = clock();
 200         double compare_time = double(after_compare - before_compare) / CLOCKS_PER_SEC;
 201         std::cout << "Matching time: " << compare_time << std::endl;
 202
 203         // Print the common tree with common edges
 204         printCommonTreeBetweenTwoTrees(aVector[0], aVector[1], vector_pair_edges_A_to_B, Q, F);
 205         printMatchingResultToFile(aVector[0], aVector[1], map_A_to_B, map_B_to_A, Q, F);
 206         //printSeparatedCommonTrees(aVector[0], aVector[1], vector_pair_edges_A_to_B, Q, F);
 207         //reconstructAndSaveCommonTrees(aVector[0], aVector[1], vector_pair_edges_A_to_B, Q, F);
 208
 209         // Print all the maps that have more that two connections
 210         std::cout << "Printing multiple relation A to B. Total relations: " << map_A_to_B.size() << std::endl;
 211         for(std::map< unsigned int, std::vector<Node*> >::iterator it_map_A_to_B = map_A_to_B.begin(); it_map_A_to_B != map_A_to_B.end(); ++it_map_A_to_B)
 212         {
 213           if((*it_map_A_to_B).second.size() > 1)
 214           {
 215             std::cout << "Multiple relation A to B with size:" << (*it_map_A_to_B).second.size() << ", from Id:" << (*it_map_A_to_B).first << std::endl;
 216             for(std::vector<Node*>::iterator it_node = (*it_map_A_to_B).second.begin(); it_node != (*it_map_A_to_B).second.end(); ++it_node)
 217             {
 218               std::cout << "Id: " << (*it_node)->GetId() << std::endl;
 219             }
 220           }
 221         }
 222         std::cout << "Printing multiple relation A to B ... OK" << std::endl;
 223
 224         std::cout << "Printing multiple relation B to A. Total relations: " << map_B_to_A.size() << std::endl;
 225         for(std::map< unsigned int, std::vector<Node*> >::iterator it_map_B_to_A = map_B_to_A.begin(); it_map_B_to_A != map_B_to_A.end(); ++it_map_B_to_A)
 226         {
 227           if((*it_map_B_to_A).second.size() > 1)
 228           {
 229             std::cout << "Multiple relation B to A with size:" << (*it_map_B_to_A).second.size() << ", from Id:" << (*it_map_B_to_A).first << std::endl;
 230             for(std::vector<Node*>::iterator it_node = (*it_map_B_to_A).second.begin(); it_node != (*it_map_B_to_A).second.end(); ++it_node)
 231             {
 232               std::cout << "Id: " << (*it_node)->GetId() << std::endl;
 233             }
 234           }
 235         }
 236         std::cout << "Printing multiple relation B to A  ... OK" << std::endl;
 237
 238         // Mark only the common nodes
 239         aVector[0].UnMarkAll();
 240         aVector[1].UnMarkAll();
 241
 242         // --------------------------------------
 243         // ------ Get common paths marked -------
 244         // --------------------------------------
 245
 246         /*
 247           std::string path_folder = "/run/media/alfredo/Data/Pulmones/results/airwaysGraphs/comparisonResutls/probes/";
 248           std::string suffix_vtk = ".vtk";
 249           std::string name_tree = "TreeA_dP2P";
 250           int iteration = 1;
 251           for(vec_nodes::iterator it_A = commonA.begin(); it_A != commonA.end(); ++it_A)
 252           {
 253           aVector[0].MarkPathFromNodeToNode(aVector[0].GetRoot(), (*it_A));
 254
 255           std::stringstream filepath_actualIteration;
 256           filepath_actualIteration << path_folder << name_tree << "__" << iteration << "_idA_" << (*it_A)->GetId() << suffix_vtk;
 257
 258           DrawVTKLinesFromTree(aVector[0], filepath_actualIteration.str(), true);
 259           ++iteration;
 260           }
 261
 262           iteration = 1;
 263           name_tree = "TreeB_dP2P";
 264           for(vec_nodes::iterator it_B = commonB.begin(); it_B != commonB.end(); ++it_B)
 265           {
 266           aVector[1].MarkPathFromNodeToNode(aVector[1].GetRoot(), (*it_B));
 267
 268           std::stringstream filepath_actualIteration;
 269           filepath_actualIteration << path_folder << name_tree << "__" << iteration << "_idB_" << (*it_B)->GetId() << suffix_vtk;
 270
 271           DrawVTKLinesFromTree(aVector[1], filepath_actualIteration.str(), true);
 272           ++iteration;
 273           }
 274         */
 275
 276         // --------------------------------------
 277         // --------------------------------------
 278
 279         /*
 280         //XXXXXXXXXXXXXXXXXXXXXXXXXXXX
 281         //1. PRINT EACH NODE OF EACH TREE AND SAVE IT USING AS NAME ITS ID
 282
 283         // Tree A
 284         for(int i = 2; i < aVector[0].GetWeight(); ++i)
 285         {
 286         //std::cout << "Writing idA: " << i << std::endl;
 287         AirwaysTree* tree_id = aVector[0].GetSingleDetachedTreeNodeById(i);
 288         //std::cout << "Got id: " << i << ", numberOfNodes: " << tree_id->GetWeight() <<  std::endl;
 289         if(tree_id)
 290         {
 291         std::stringstream filepath_actualIteration;
 292         filepath_actualIteration << aVector[0].GetResultPath() << "/" << aVector[0].GetImageName() << "_id_" << i << suffix_vtk;
 293         DrawVTKLinesFromTree(*tree_id, filepath_actualIteration.str(), false);
 294         delete(tree_id);
 295         }
 296         else
 297         std::cout << "Tree NULL" << std::endl;
 298         }
 299
 300         // Tree B
 301         for(int i = 2; i < aVector[1].GetWeight(); ++i)
 302         {
 303         //std::cout << "Writing idA: " << i << std::endl;
 304         AirwaysTree* tree_id = aVector[1].GetSingleDetachedTreeNodeById(i);
 305         //std::cout << "Got id: " << i << ", numberOfNodes: " << tree_id->GetWeight() <<  std::endl;
 306         if(tree_id)
 307         {
 308         std::stringstream filepath_actualIteration;
 309         filepath_actualIteration << aVector[1].GetResultPath() << "/" << aVector[1].GetImageName() << "_id_" << i << suffix_vtk;
 310         DrawVTKLinesFromTree(*tree_id, filepath_actualIteration.str(), false);
 311         delete(tree_id);
 312         }
 313         else
 314         std::cout << "Tree NULL" << std::endl;
 315         }
 316
 317         std::cout << "Comparing trees Orkisz-Morales... OK" << std::endl;
 318
 319         */
 320
 321         //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 322         //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 323         //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 324
 325         /*// AMP
 326           aVector[0].SubIsomorphism(aVector[1]);
 327
 328           std::cout << "Flag after SubIsomorphism" << std::endl;
 329           std::string fullPathA = aVector[0].GetResultPath() + "/" + aVector[0].GetImageName() + "_CMP.vtk";
 330           std::string fullPathB = aVector[1].GetResultPath() + "/" + aVector[1].GetImageName() + "_CMP.vtk";
 331           DrawVTKLinesFromTree(aVector[0], fullPathA, true);
 332           DrawVTKLinesFromTree(aVector[1], fullPathB, true);
 333
 334           std::cout << "Flag after write vtk" << std::endl;
 335           TInputImage::Pointer imgA, imgB;
 336           aVector[0].ImageReconstruction(imgA);
 337           std::cout << "Flag after Reconstruction A" << std::endl;
 338
 339           aVector[1].ImageReconstruction(imgB);
 340           std::cout << "Flag after Reconstruction B" << std::endl;
 341
 342           std::string fullPathA_mhd = aVector[0].GetResultPath() + "/" + aVector[0].GetImageName() + "_CMP.vtk";
 343           std::string fullPathB_mhd = aVector[1].GetResultPath() + "/" + aVector[1].GetImageName() + "_CMP.vtk";
 344           WriteImage(imgA, fullPathA_mhd);
 345           WriteImage(imgB, fullPathB_mhd);
 346
 347           //AMP*/
 348       }
 349       po::notify(vm); // throws on error, so do after help in case
 350     }
 351     catch (std::exception& e)
 352     {
 353       std::cerr << "Unhandled Exception reached the top of main: " << e.what() << ", application will now exit" << std::endl;
 354       return ERROR_UNHANDLED_EXCEPTION;
 355     } // yrt
 356   }
 357   catch (std::exception& e)
 358   {
 359     std::cerr << "Unhandled Exception reached the top of main: " << e.what() << ", application will now exit" << std::endl;
 360     return ERROR_UNHANDLED_EXCEPTION;
 361   } // yrt
 362   return( 0 );
 363 }
 364
 365 // -------------------------------------------------------------------------
 366 void Load_cpPlugins( const std::string& plugins )
 367 {
 368   myPlugins.LoadConfiguration( plugins );
 369 }
 370
 371 // -------------------------------------------------------------------------
 372 void CreateResultDirectory(AirwaysTree tree)
 373 {
 374   boost::filesystem::path dir(tree.GetResultPath());
 375   if (boost::filesystem::create_directories(dir))
 376   {
 377     std::cout << "FolderName = " << tree.GetResultPath() << " has been created" << std::endl;
 378   } //if
 379   else
 380     std::cout << "FolderName = " << tree.GetResultPath() << " already exists" << std::endl;
 381 }
 382
 383 // -------------------------------------------------------------------------
 384 void DrawVTKLinesFromTree(AirwaysTree& tree, const std::string filepath, bool common)
 385 {
 386   // Create the array of points, lines, and colors
 387   vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
 388   vtkSmartPointer<vtkCellArray> lines = vtkSmartPointer<vtkCellArray>::New();
 389   vtkSmartPointer<vtkUnsignedCharArray> colors = vtkSmartPointer<vtkUnsignedCharArray>::New();
 390   colors->SetNumberOfComponents(3);
 391   colors->SetName("Colors");
 392
 393   //vtk sphere for sources
 394   //vtkSmartPointer<vtkSphereSource> sphereSource = vtkSmartPointer<vtkSphereSource>::New();
 395   //Vec3 root = tree.GetRoot()->GetCoords();
 396   //sphereSource->SetCenter(root[0], root[1], root[2]);
 397   //sphereSource->SetRadius(1);
 398   //UndirectedGraph
 399
 400   srand(time(NULL));
 401   unsigned int id = 1;
 402
 403   // Create and fill the points, lines, and color vectors
 404   //CalculateVTKLinesFromEdges(tree.GetRoot(), 0, id, pts, lines, colors, common);
 405   pts->SetNumberOfPoints(tree.GetWeight()+1);
 406   createLinesAndPointsForVTK(tree.GetRoot(), pts, lines, colors, common, true);
 407
 408   // Create the polydata
 409   vtkSmartPointer<vtkPolyData> linesPolyData = vtkSmartPointer<vtkPolyData>::New();
 410
 411   // Set the points, lines, and colors
 412   linesPolyData->SetPoints(pts);
 413   linesPolyData->SetLines(lines);
 414   linesPolyData->GetCellData()->SetScalars(colors);
 415
 416   // Write the file
 417   vtkSmartPointer<vtkPolyDataWriter> writer = vtkSmartPointer<vtkPolyDataWriter>::New();
 418   writer->SetFileName(filepath.c_str());
 419
 420 #if VTK_MAJOR_VERSION <= 5
 421   writer->SetInput(linesPolyData);
 422 #else
 423   writer->SetInputData(linesPolyData);
 424 #endif
 425   writer->Write();
 426
 427   //The following code is to test the vtk polydata and visualize it
 428   /*    // Visualize
 429         vtkSmartPointer<vtkPolyDataMapper> mapper =
 430         vtkSmartPointer<vtkPolyDataMapper>::New();
 431         #if VTK_MAJOR_VERSION <= 5
 432         mapper->SetInput(linesPolyData);
 433         #else
 434         mapper->SetInputData(linesPolyData);
 435         #endif
 436         vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
 437         actor->SetMapper(mapper);
 438
 439         //sphere
 440         vtkSmartPointer<vtkPolyDataMapper> mapperSphere = vtkSmartPointer<
 441         vtkPolyDataMapper>::New();
 442         mapperSphere->SetInputConnection(sphereSource->GetOutputPort());
 443
 444         vtkSmartPointer<vtkActor> actorSphere = vtkSmartPointer<vtkActor>::New();
 445         actorSphere->SetMapper(mapperSphere);
 446         //end sphere
 447
 448         vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
 449         vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<
 450         vtkRenderWindow>::New();
 451         renderWindow->AddRenderer(renderer);
 452         vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
 453         vtkSmartPointer<vtkRenderWindowInteractor>::New();
 454         renderWindowInteractor->SetRenderWindow(renderWindow);
 455         renderer->AddActor(actorSphere);
 456         renderer->AddActor(actor);
 457
 458         renderWindow->Render();
 459         renderWindowInteractor->Start();*/
 460 }
 461
 462 // -------------------------------------------------------------------------
 463 #include <fpa/Base/ImageSkeleton.h>
 464 #include <fpa/Image/MinimumSpanningTree.h>
 465 #include <cpExtensions/DataStructures/ImageIndexesContainer.h>
 466 #include <itkImage.h>
 467
 468 template< class TImage, class TVertex >
 469 Node* FAVertexToNode( TVertex vertex, TImage* image )
 470 {
 471   //The FrontAlgorithms Vertex is an TImageType::IndexType
 472   typename TImage::PointType point;
 473   image->TransformIndexToPhysicalPoint( vertex, point );
 474   Vec3 alfPoint(point[0],point[1],point[2]);
 475   return new Node(alfPoint);
 476 }
 477
 478 AirwaysTree& ConvertFilterToAirwaysTree( TInputImage* input_image, cpPlugins::Workspace& ws )
 479 {
 480   /*
 481     typedef FilterType TFilter;
 482     typedef typename TFilter::TVertex TVertexFA;
 483     typedef typename TFilter::TVertices TVerticesFA;
 484     typedef typename TFilter::TUniqueVertices TUniqueVerticesFA;
 485     typedef typename TFilter::TVertexCompare TVertexCompareFA;
 486     typedef typename TFilter::TBranches TBranchesFA;
 487     typedef typename TFilter::TMinimumSpanningTree TMst;
 488   */
 489   typedef
 490     fpa::Base::ImageSkeleton< fpa::Image::MinimumSpanningTree< 3 > >
 491     TFASkeleton;
 492   typedef TFASkeleton::TVertex    TVertexFA;
 493   typedef TFASkeleton::TVertexCmp TVertexCompareFA;
 494   typedef cpExtensions::DataStructures::ImageIndexesContainer< 3 > TVerticesFA;
 495
 496   typedef Node TVertexAirways;
 497   typedef Edge TEdgeAirways;
 498   typedef pair_posVox_rad TSkelePoint;
 499   typedef vec_pair_posVox_rad TSkelePoints;
 500   typedef std::map< TVertexFA, TVertexAirways*, TVertexCompareFA > VertexMap;
 501   VertexMap vertexMap;
 502   std::time_t start,end;
 503   std::time(&start);
 504   std::cout << "Starting conversion " << std::endl;
 505
 506   auto w_filter = ws.GetFilter( "eb" );
 507   auto branches = w_filter->GetOutputData( "Skeleton" )->GetITK< TFASkeleton >( );
 508   auto& endpoints = w_filter->GetOutputData( "EndPoints" )->GetITK< TVerticesFA >( )->Get( );
 509   auto& bifurcations = w_filter->GetOutputData( "Bifurcations" )->GetITK< TVerticesFA >( )->Get( );
 510   auto image = ws.GetFilter( "reader" )->GetOutputData( "Output" )->GetITK< itk::ImageBase< 3 > >( );
 511   auto distance_map = ws.GetFilter( "dmap" )->GetOutputData( "Output" )->GetITK< itk::Image< float, 3 > >( );
 512
 513   /*
 514     TBranchesFA* branches = filter->GetBranches( );
 515     TMst* mst = filter->GetMinimumSpanningTree();
 516   */
 517   int leoTreeWeigth = endpoints.size()+bifurcations.size()+1;
 518   std::cout<< "Creates FA Tree with : "<<leoTreeWeigth << " nodes "<<std::endl;
 519
 520   auto seed0 = ws.GetFilter( "seed" )->GetOutputData( "Output" )->GetITK< TVerticesFA >( )->Get( )[ 0 ];
 521   //Fill vertex map. Gotta do it with bifurcations, endpoints and the root.
 522   //Do the root
 523   vertexMap[ seed0 ] = FAVertexToNode( seed0, image );
 524   //Do Endpoints
 525   auto eIt = endpoints.begin();
 526   for( ; eIt != endpoints.end(); ++eIt )
 527   {
 528     vertexMap[*eIt]=FAVertexToNode( *eIt, image );
 529   }
 530
 531   auto biIt = bifurcations.begin();
 532   for(; biIt != bifurcations.end(); ++biIt)
 533   {
 534     vertexMap[*biIt]=FAVertexToNode(*biIt,image);
 535   }
 536
 537   //Now navigate branches to make the edges of the tree.
 538   auto bIt = branches->Get( ).begin();
 539   for(; bIt!=branches->Get( ).end(); ++bIt)
 540   {
 541     auto dVertex = bIt->first;
 542     TVertexAirways* destination = vertexMap[dVertex];
 543     if( destination == NULL )
 544     {
 545       vertexMap[dVertex]=FAVertexToNode(dVertex,image);
 546       destination = vertexMap[dVertex];
 547     }
 548     auto brIt = bIt->second.begin( );
 549     for( ; brIt != bIt->second.end( ); ++brIt )
 550     {
 551       auto sVertex = brIt->first;
 552       TVertexAirways* source = vertexMap[sVertex];
 553       if( source == NULL )
 554       {
 555         vertexMap[sVertex]=FAVertexToNode(sVertex,image);
 556         source = vertexMap[sVertex];
 557       }
 558       destination->SetFather(source);
 559       source->AddChild(destination);
 560       TEdgeAirways* edge = new Edge();
 561       //Update path info for the edge. This is why we don't need to call UpdateEdges on the constructor of the tree.
 562
 563       auto edgePath = brIt->second->GetVertexList( );
 564       for( unsigned int pIt = 0; pIt < edgePath->Size( ); ++pIt )
 565       {
 566         itk::ImageBase< 3 >::PointType pnt;
 567         auto cidx = edgePath->GetElement( pIt );
 568         image->TransformContinuousIndexToPhysicalPoint(cidx,pnt);
 569         TVertexFA idx;
 570         image->TransformPhysicalPointToIndex(pnt,idx);
 571         Vec3 coords = FAVertexToNode(idx,image)->GetCoords();
 572         pair_posVox_rad skPair(coords,distance_map->GetPixel(idx));
 573         edge->AddSkeletonPairInfo(skPair);
 574       }
 575     }
 576   }
 577   AirwaysTree* tree = new AirwaysTree( input_image, NULL, vertexMap[seed0], false);
 578   std::time(&end);
 579   std::cout << input_image << std::endl;
 580   std::cout << "Finished conversion. New AlfTree has weight: "<<tree->GetWeight()<<". Takes "<<(end-start)<<" s." << std::endl;
 581   return *tree;
 582 }
 583
 584 // -------------------------------------------------------------------------
 585 AirwaysTree& CreateAirwaysTreeFromSegmentation(Vec3 seed, TInputImage* input_image, cpPlugins::Workspace& ws )
 586 {
 587   return( ConvertFilterToAirwaysTree( input_image, ws ) );
 588 }
 589
 590 // -------------------------------------------------------------------------
 591 vector<TreeInfo> ReadInputFile(const char* filename)
 592 {
 593   // Variables
 594   std::ifstream infile(filename);
 595   std::string line;
 596   std::string folderpath_allResults;
 597   vector<TreeInfo> vectorInfo;
 598   // Parse the file
 599   bool firstLine = false;
 600
 601   while (std::getline(infile, line))
 602   {
 603     // First line contains the output folder
 604     std::istringstream iss(line);
 605     if (!firstLine)
 606     {
 607       if (!(iss >> folderpath_allResults))
 608       {
 609         std::cout << "no file" << std::endl;
 610         return vectorInfo;
 611       } //if
 612       firstLine = true;
 613     } //if
 614     // Other lines, not the first one, contain the information for the airways to be created
 615     else
 616     {
 617       TreeInfo treeInfo;
 618       float point_trachea[3];
 619       std::string filepath_airwayImage, name_pig, name_image;
 620       if (!(iss >> point_trachea[0] >> point_trachea[1] >> point_trachea[2] >> filepath_airwayImage >> name_pig >> name_image))
 621       {
 622         std::cout << "There is no tree information in the file." << std::endl;
 623         break;
 624       } //if
 625       else
 626       {
 627         // Print information
 628         std::cout << "Point trachea:[" << point_trachea[0] << "," << point_trachea[1] << "," << point_trachea[2] << "]" << std::endl;
 629       }
 630
 631       //Save info.
 632
 633       //Save seed info.
 634       Vec3 seed(point_trachea[0], point_trachea[1], point_trachea[2]); //real coords seed
 635       treeInfo.seed = seed;
 636       // Create the outputs
 637       treeInfo.folderpath_pigResults = folderpath_allResults + "/" + name_pig + "/";
 638       treeInfo.pig_name = name_pig;
 639       treeInfo.image_name = name_image;
 640       // Read the image
 641       /*
 642         treeInfo.image = ReadImage<TInputImage>(filepath_airwayImage);
 643       */
 644       // Execute first pipeline's part
 645       std::stringstream seed_stream;
 646       seed_stream
 647         << point_trachea[0] << " "
 648         << point_trachea[1] << " "
 649         << point_trachea[2];
 650
 651       treeInfo.myWorkspace->SetInterface( &myPlugins );
 652       std::string err = treeInfo.myWorkspace->LoadWorkspace( "./workspace_airwaysappli.wxml" );
 653       if( err != "" )
 654       {
 655         std::cerr << "Error: " << err << std::endl;
 656         std::exit( 1 );
 657
 658       } // fi
 659       treeInfo.myWorkspace->SetParameter( "FileNames@reader", filepath_airwayImage );
 660       treeInfo.myWorkspace->SetParameter( "Text@seed", seed_stream.str( ) );
 661       vectorInfo.push_back(treeInfo);
 662
 663     } //else
 664   } //while
 665   return vectorInfo;
 666 }
 667
 668 // -------------------------------------------------------------------------
 669 void printCommonTreeBetweenTwoTrees(AirwaysTree tree_A, AirwaysTree tree_B, std::vector< std::pair< std::pair<Node*, Node*>, std::pair<Node*, Node*> > > vector_pair_edges_A_to_B, unsigned int Q, unsigned int F)
 670 {
 671   std::cout << "printCommonTreeBetweenTwoTrees, edges:" << vector_pair_edges_A_to_B.size() << std::endl;
 672
 673   // Vtk points, cell array, and colors
 674   vtkSmartPointer<vtkPoints> points_common = vtkSmartPointer<vtkPoints>::New();
 675   vtkSmartPointer<vtkCellArray> lines_common = vtkSmartPointer<vtkCellArray>::New();
 676   vtkSmartPointer<vtkUnsignedCharArray> colors_common = vtkSmartPointer<vtkUnsignedCharArray>::New();
 677   colors_common->SetNumberOfComponents(3);
 678   colors_common->SetName("Colors");
 679
 680   // Add all the points
 681   // 0 - 1000 points for tree A and from 1001 for tree B
 682   points_common->SetNumberOfPoints(1000 + tree_B.GetWeight() + 100);
 683
 684   // Add points for tree A
 685   vec_nodes vector_nodesA = tree_A.GetNodes();
 686   for(vec_nodes::iterator it_nodesA = vector_nodesA.begin(); it_nodesA != vector_nodesA.end(); ++it_nodesA)
 687   {
 688     points_common->SetPoint((*it_nodesA)->GetId(),(*it_nodesA)->GetCoords().GetVec3());
 689   }
 690
 691   // Add points for tree B
 692   vec_nodes vector_nodesB = tree_B.GetNodes();
 693   for(vec_nodes::iterator it_nodesB = vector_nodesB.begin(); it_nodesB != vector_nodesB.end(); ++it_nodesB)
 694   {
 695     points_common->SetPoint((*it_nodesB)->GetId()+1000,(*it_nodesB)->GetCoords().GetVec3());
 696   }
 697
 698   // Add the edges for both trees
 699   std::vector< std::pair< std::pair<Node*, Node*>, std::pair<Node*, Node*> > >::iterator it_edges = vector_pair_edges_A_to_B.begin();
 700   int number_pairs = 0;
 701   for(; it_edges != vector_pair_edges_A_to_B.end(); ++it_edges)
 702   {
 703     std::cout << "Pair:" << number_pairs << std::endl;
 704     std::pair<Node*, Node*> edge_A = (*it_edges).first;
 705     std::pair<Node*, Node*> edge_B = (*it_edges).second;
 706
 707     vec_nodes path_A;
 708     edge_A.first->GetPathToNode(edge_A.second, path_A);
 709
 710     vec_nodes path_B;
 711     edge_B.first->GetPathToNode(edge_B.second, path_B);
 712
 713     // Set color to be used
 714     int numColor = number_pairs % 6;
 715     unsigned char color[3];
 716     color[0] = green[0];
 717     color[1] = green[1];
 718     color[2] = green[2];
 719     if(numColor == 1)
 720     {
 721       color[0] = red[0];
 722       color[1] = red[1];
 723       color[2] = red[2];
 724     }
 725     else if(numColor == 2)
 726     {
 727       color[0] = blue[0];
 728       color[1] = blue[1];
 729       color[2] = blue[2];
 730     }
 731     else if(numColor == 3)
 732     {
 733       color[0] = yellow[0];
 734       color[1] = yellow[1];
 735       color[2] = yellow[2];
 736     }
 737     else if(numColor == 4)
 738     {
 739       color[0] = purple[0];
 740       color[1] = purple[1];
 741       color[2] = purple[2];
 742     }
 743     else if(numColor == 5)
 744     {
 745       color[0] = cyan[0];
 746       color[1] = cyan[1];
 747       color[2] = cyan[2];
 748     }
 749
 750     number_pairs++;
 751
 752     // Trace line A
 753     //if(path_A.size() > 0 && number_pairs < 50)
 754     if( path_A.size() )
 755     {
 756       vtkSmartPointer<vtkLine> line = vtkSmartPointer<vtkLine>::New();
 757       lines_common->InsertNextCell(path_A.size());
 758       int id_point = 0;
 759       for(vec_nodes::iterator it_pathA = path_A.begin(); it_pathA != path_A.end(); ++it_pathA)
 760       {
 761         lines_common->InsertCellPoint((*it_pathA)->GetId());
 762         //line->GetPointIds()->SetId( id_point, (*it_pathA)->GetId() );
 763         //id_point++;
 764       }
 765       //lines_common->InsertNextCell(line);
 766       colors_common->InsertNextTupleValue(color);
 767     }
 768     else
 769     {
 770       std::cout << "No path A" << std::endl;
 771     }
 772
 773
 774     // Trace line B
 775     //if(path_B.size() > 0 && number_pairs < 50)
 776     if( path_B.size() )
 777     {
 778       vtkSmartPointer<vtkLine> line = vtkSmartPointer<vtkLine>::New();
 779       lines_common->InsertNextCell(path_B.size());
 780       int id_point = 0;
 781       for(vec_nodes::iterator it_pathB = path_B.begin(); it_pathB != path_B.end(); ++it_pathB)
 782       {
 783         lines_common->InsertCellPoint((*it_pathB)->GetId()+1000);
 784         //line->GetPointIds()->SetId( id_point, (*it_pathB)->GetId()+1000 );
 785         id_point++;
 786       }
 787       //lines_common->InsertNextCell(line);
 788       colors_common->InsertNextTupleValue(color);
 789     }
 790     else
 791     {
 792       std::cout << "No path B" << std::endl;
 793     }
 794   }
 795
 796   // Save the polydata
 797   // Create the polydata
 798   vtkSmartPointer<vtkPolyData> linesPolyData = vtkSmartPointer<vtkPolyData>::New();
 799
 800   // Set the points, lines, and colors
 801   linesPolyData->SetPoints(points_common);
 802   linesPolyData->SetLines(lines_common);
 803   linesPolyData->GetCellData()->SetScalars(colors_common);
 804
 805   // ------------------------------------------------
 806   // Write the vtk file
 807
 808   // Create the pathfile to save
 809   std::stringstream filepath_actualIteration;
 810   filepath_actualIteration << tree_A.GetResultPath() << "/" << tree_A.GetImageName() << "_" << tree_B.GetImageName() << "_Q_" << Q << "_F_" << F << ".vtk";
 811   std::cout << "File to save:" << filepath_actualIteration.str() << std::endl;
 812
 813   // Create the writer
 814   vtkSmartPointer<vtkPolyDataWriter> writer = vtkSmartPointer<vtkPolyDataWriter>::New();
 815   writer->SetFileName(filepath_actualIteration.str().c_str());
 816
 817   // Set input and write
 818 #if VTK_MAJOR_VERSION <= 5
 819   writer->SetInput(linesPolyData);
 820 #else
 821   writer->SetInputData(linesPolyData);
 822 #endif
 823   writer->Write();
 824
 825
 826
 827   // ***************************************
 828   // Save the links in a file
 829   // *******
 830
 831   // Create the outputfile
 832   std::stringstream filepath_evaluation;
 833   filepath_evaluation << tree_A.GetResultPath() << "/" << tree_A.GetImageName() << "_" << tree_B.GetImageName() << "_Q_" << Q << "_F_" << F << ".csv";
 834
 835   std::ofstream file(filepath_evaluation.str().c_str());
 836   if(file.is_open())
 837   {
 838     // Save the header
 839     //file << "Node_" << tree_A.GetImageName() << " " << "Node_" << tree_B.GetImageName()<< "\n";
 840     file << "TreeName "
 841       "idLocalN1 idLocalN2 "
 842       "Node1x Node1y Node1z "
 843       "Node2x Node2y Node2z "
 844       "idGlobalN1 "
 845       "idMatch" << std::endl;
 846
 847
 848     std::vector< std::pair< std::pair<Node*, Node*>, std::pair<Node*, Node*> > >::iterator it_edges = vector_pair_edges_A_to_B.begin();
 849     for(; it_edges != vector_pair_edges_A_to_B.end(); ++it_edges)
 850     {
 851
 852       std::pair<Node*, Node*> edge_A = (*it_edges).first;
 853       std::pair<Node*, Node*> edge_B = (*it_edges).second;
 854
 855       //file << edge_A.second->GetId() << " " << edge_B.second->GetId()+1000 << "\n";
 856       Vec3 coord_EndNodeA = edge_A.second->GetCoords();
 857       Vec3 coord_EndNodeB = edge_B.second->GetCoords();
 858       file << tree_A.GetImageName() << " " <<
 859         edge_A.second->GetId() << " " << edge_B.second->GetId()+1000 << " " <<
 860         coord_EndNodeA[0] << " " << coord_EndNodeA[1] << " " << coord_EndNodeA[2] << " " <<
 861         coord_EndNodeB[0] << " " << coord_EndNodeB[1] << " " << coord_EndNodeB[2] << " " <<
 862         tree_A.GetImageName() << coord_EndNodeA[0] << coord_EndNodeA[1] << coord_EndNodeA[2] << " " <<
 863         tree_B.GetImageName() << coord_EndNodeB[0] << coord_EndNodeB[1] << coord_EndNodeB[2] <<"\n";
 864
 865       file << tree_B.GetImageName() << " " <<
 866         edge_B.second->GetId()+1000 << " " << edge_A.second->GetId() << " " <<
 867         coord_EndNodeB[0] << " " << coord_EndNodeB[1] << " " << coord_EndNodeB[2] << " " <<
 868         coord_EndNodeA[0] << " " << coord_EndNodeA[1] << " " << coord_EndNodeA[2] << " " <<
 869         tree_B.GetImageName() << coord_EndNodeB[0] << coord_EndNodeB[1] << coord_EndNodeB[2] << " " <<
 870         tree_A.GetImageName() << coord_EndNodeA[0] << coord_EndNodeA[1] << coord_EndNodeA[2] << "\n";
 871     }
 872   }
 873
 874   file.close();
 875   // *******
 876   // ***************************************
 877
 878
 879
 880
 881   std::cout << "printCommonTreeBetweenTwoTrees ... OK" << std::endl;
 882 }
 883
 884 // -------------------------------------------------------------------------
 885 void printMatchingResultToFile(AirwaysTree tree_A, AirwaysTree tree_B, std::map< unsigned int, std::vector<Node*> > map_A_to_B, std::map< unsigned int, std::vector<Node*> > map_B_to_A, unsigned int Q, unsigned int F)
 886 {
 887   std::cout << "Printing matching result ... " << std::endl;
 888
 889   // Variables and types
 890   typedef std::map< unsigned int, vec_nodes > map_id_node;
 891
 892   // Create the outputfile
 893   std::stringstream filepath_evaluation;
 894   filepath_evaluation << tree_A.GetResultPath() << "/" << tree_A.GetImageName() << "_" << tree_B.GetImageName() << "_Q_" << Q << "_F_" << F << "_matching.csv";
 895
 896   std::ofstream file(filepath_evaluation.str().c_str());
 897   if(file.is_open())
 898   {
 899     // Save the header
 900     file << "TreeName idLocalN1 idLocalN2 "
 901       "Node1x Node1y Node1z "
 902       "Node2x Node2y Node2z "
 903       "idGlobalN1 "
 904       "idMatch "
 905       "typeMatch_match_1_nonmatch_0 "
 906       "depth "
 907       "leaf"<< std::endl;
 908
 909     // ---------------
 910     // From A to B
 911
 912     // Save the match or non-match for each node from A to B
 913     for(int id_a=1; id_a <= tree_A.GetWeight( ); ++id_a)
 914     {
 915       // Get the actual node in tree A
 916       Node* node_a = tree_A.GetNodeById( id_a );
 917       Vec3 coords_a = node_a->GetCoords();
 918       unsigned int depth_a = tree_A.GetDepthById(id_a);
 919       bool leaf_a = node_a->IsLeaf();
 920
 921       // Check if the node was matched
 922       map_id_node::iterator it_a2b = map_A_to_B.find( id_a );
 923       if( it_a2b != map_A_to_B.end( ) )
 924       {
 925         // Get the correspoding matching nodes and print them
 926         vec_nodes nodes_B = (*it_a2b).second;
 927
 928         vec_nodes::iterator it_nodes_b = nodes_B.begin( );
 929         for( ; it_nodes_b != nodes_B.end( ); ++it_nodes_b)
 930         {
 931           Vec3 coords_b = (*it_nodes_b)->GetCoords();
 932
 933           file << tree_A.GetImageName() << " " << id_a << " " << (*it_nodes_b)->GetId()+1000 << " "
 934                << coords_a[0] << " " << coords_a[1] << " " << coords_a[2] << " "
 935                << coords_b[0] << " " << coords_b[1] << " " << coords_b[2] << " "
 936                << tree_A.GetImageName() << coords_a[0] << coords_a[1] << coords_a[2] << " "
 937                << tree_A.GetImageName() << coords_a[0] << coords_a[1] << coords_a[2] << coords_b[0] << coords_b[1] << coords_b[2] << " "
 938                << "1" << " "
 939                << depth_a << " "
 940                << leaf_a << "\n";
 941         }
 942       }
 943       else
 944       {
 945         file << tree_A.GetImageName() << " " << id_a << " " << "0" << " "
 946              << coords_a[0] << " " << coords_a[1] << " " << coords_a[2] << " "
 947              << "0" << " " << "0" << " " << "0" << " "
 948              << tree_A.GetImageName() << coords_a[0] << coords_a[1] << coords_a[2] << " "
 949              << tree_A.GetImageName() << coords_a[0] << coords_a[1] << coords_a[2] << "0" << "0" << "0" << " "
 950              << "0" << " "
 951              << depth_a << " "
 952              << leaf_a  << "\n";
 953       }
 954     }
 955
 956     // ---------------
 957     // From B to A
 958
 959     // Save the match or non-match for each node from A to B
 960     for(int id_b=1; id_b <= tree_B.GetWeight( ); ++id_b)
 961     {
 962       // Get the actual node in tree B
 963       Node* node_b = tree_B.GetNodeById( id_b );
 964       Vec3 coords_b = node_b->GetCoords();
 965       unsigned int depth_b = tree_B.GetDepthById(id_b);
 966       bool leaf_b = node_b->IsLeaf();
 967
 968       // Check if the node was matched
 969       map_id_node::iterator it_b2a = map_B_to_A.find( id_b );
 970       if( it_b2a != map_B_to_A.end( ) )
 971       {
 972         // Get the correspoding matching nodes and print them
 973         vec_nodes nodes_A = (*it_b2a).second;
 974
 975         vec_nodes::iterator it_nodes_a = nodes_A.begin( );
 976         for( ; it_nodes_a != nodes_A.end( ); ++it_nodes_a)
 977         {
 978           Vec3 coords_a = (*it_nodes_a)->GetCoords();
 979
 980           file << tree_B.GetImageName() << " " << id_b+1000 << " " << (*it_nodes_a)->GetId() << " "
 981                << coords_b[0] << " " << coords_b[1] << " " << coords_b[2] << " "
 982                << coords_a[0] << " " << coords_a[1] << " " << coords_a[2] << " "
 983                << tree_B.GetImageName() << coords_b[0] << coords_b[1] << coords_b[2] << " "
 984                << tree_B.GetImageName() << coords_b[0] << coords_b[1] << coords_b[2] << coords_a[0] << coords_a[1] << coords_a[2] << " "
 985                << "1" << " "
 986                << depth_b << " "
 987                << leaf_b  << "\n";
 988         }
 989       }
 990       else
 991       {
 992         file << tree_B.GetImageName() << " " << id_b+1000 << " " << "0" << " "
 993              << coords_b[0] << " " << coords_b[1] << " " << coords_b[2] << " "
 994              << "0" << " " << "0" << " " << "0" << " "
 995              << tree_B.GetImageName() << coords_b[0] << coords_b[1] << coords_b[2] << " "
 996              << tree_B.GetImageName() << coords_b[0] << coords_b[1] << coords_b[2] << "0" << "0" << "0" << " "
 997              << "0" << " "
 998              << depth_b << " "
 999              << leaf_b  << "\n";
1000       }//esle
1001     }//rof
1002   }//fi
1003
1004   file.close(); // Close the file
1005
1006   std::cout << "Printing matching result DONE" << std::endl;
1007 }
1008
1009 // -------------------------------------------------------------------------
1010 void createLinesAndPointsForVTK(const Node* node, vtkSmartPointer<vtkPoints>& pts,vtkSmartPointer<vtkCellArray>& lines, vtkSmartPointer<vtkUnsignedCharArray>& colors, bool common, bool isRoot)
1011 {
1012   // Insert the actual point/node
1013   //vtkIdType id_father = idNonRoot;
1014   //if(isRoot)
1015   //id_father = pts->InsertNextPoint(node->GetCoords().GetVec3());
1016   vtkIdType id_father = node->GetId();
1017   if(isRoot)
1018     pts->SetPoint(id_father,node->GetCoords().GetVec3());
1019
1020   // Iterate over the children
1021   const vec_nodes children = node->GetChildren();
1022   for (vec_nodes::const_iterator it_child = children.begin(); it_child != children.end(); ++it_child)
1023   {
1024     if (!(*it_child)->IsMarked() && common)
1025       continue;
1026
1027     //vtkIdType id_son = pts->InsertNextPoint((*it_child)->GetCoords().GetVec3());
1028     vtkIdType id_son = (*it_child)->GetId();
1029     pts->SetPoint(id_son,(*it_child)->GetCoords().GetVec3());
1030
1031     // Set color to be used
1032     int numColor = (*it_child)->GetLevel() % 4;
1033     if(numColor == 0)
1034       colors->InsertNextTupleValue(green);
1035     else if(numColor == 1)
1036       colors->InsertNextTupleValue(red);
1037     else if(numColor == 2)
1038       colors->InsertNextTupleValue(red);
1039     else
1040       colors->InsertNextTupleValue(red);
1041
1042     // Add the line
1043     vtkSmartPointer<vtkLine> line = vtkSmartPointer<vtkLine>::New();
1044     line->GetPointIds()->SetId(0, id_father);
1045     line->GetPointIds()->SetId(1, id_son);
1046     lines->InsertNextCell(line);
1047
1048     createLinesAndPointsForVTK(*it_child, pts, lines, colors, common, false);
1049   } //for
1050 }
1051
1052 // eof - $RCSfile$