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