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