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