5 #include <boost/filesystem.hpp>
6 #include <boost/program_options.hpp>
9 #include "airwaysTree.h"
11 #include <cpPlugins/Interface.h>
12 #include <cpPlugins/Workspace.h>
14 using namespace airways;
15 namespace po = boost::program_options;
19 const size_t ERROR_IN_COMMAND_LINE = 1;
20 const size_t SUCCESS = 0;
21 const size_t ERROR_UNHANDLED_EXCEPTION = 2;
24 typedef std::vector< AirwaysTree > AirwaysVector;
26 // Auxiliar struct to save info for execution.
27 typedef void* TImagePointer;
31 std::string folderpath_pigResults;
33 std::string image_name;
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 };
44 cpPlugins::Interface myPlugins;
45 cpPlugins::Workspace myWorkspace;
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);
57 // -------------------------------------------------------------------------
58 int main( int argc, char* argv[] )
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")
75 vector<TreeInfo> treeInfoVector;
76 AirwaysVector aVector;
77 TImagePointer segmentationImage;
79 po::store(po::parse_command_line(argc, argv, desc), vm); // can throw
80 if (vm.count("use_file"))
82 fileName = vm["use_file"].as<std::string>();
83 treeInfoVector = ReadInputFile(fileName.c_str());
86 if(vm.count("segment_images"))
91 if (vm.count("build_trees"))
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);
100 for (unsigned int i = 0; i < aVector.size(); ++i)
102 CreateResultDirectory(aVector[i]);
103 std::string fullPath = aVector[i].GetResultPath() + "/" + aVector[i].GetImageName() + ".vtk";
104 DrawVTKLinesFromTree(aVector[i], fullPath, false);
109 // Subtree levels option
110 if (vm.count("subtree_levels"))
113 else if (vm.count("subtree_length"))
116 else if (vm.count("subtree_diameter"))
120 if (vm.count("compare_trees"))
122 std::cout << "Option: Compare trees" << std::endl;
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
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;
136 std::cout << "Tree A ... " << std::endl;
137 //aVector[0].printNodeAndChildrenIds();
138 std::cout << "Tree A ... OK" << std::endl;
140 std::cout << "Tree B ... " << std::endl;
141 //aVector[1].printNodeAndChildrenIds();
142 std::cout << "Tree B ... OK" << std::endl;
144 std::cout << "Comparing trees Orkisz-Morales..." << std::endl;
145 // Vectors to save the common and uncommon nodes
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;
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;
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;
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);
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)
176 if((*it_map_A_to_B).second.size() > 1)
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)
181 std::cout << "Id: " << (*it_node)->GetId() << std::endl;
185 std::cout << "Printing multiple relation A to B ... OK" << std::endl;
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)
190 if((*it_map_B_to_A).second.size() > 1)
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)
195 std::cout << "Id: " << (*it_node)->GetId() << std::endl;
199 std::cout << "Printing multiple relation B to A ... OK" << std::endl;
201 // Mark only the common nodes
202 aVector[0].UnMarkAll();
203 aVector[1].UnMarkAll();
205 // --------------------------------------
206 // ------ Get common paths marked -------
207 // --------------------------------------
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";
214 for(vec_nodes::iterator it_A = commonA.begin(); it_A != commonA.end(); ++it_A)
216 aVector[0].MarkPathFromNodeToNode(aVector[0].GetRoot(), (*it_A));
218 std::stringstream filepath_actualIteration;
219 filepath_actualIteration << path_folder << name_tree << "__" << iteration << "_idA_" << (*it_A)->GetId() << suffix_vtk;
221 DrawVTKLinesFromTree(aVector[0], filepath_actualIteration.str(), true);
226 name_tree = "TreeB_dP2P";
227 for(vec_nodes::iterator it_B = commonB.begin(); it_B != commonB.end(); ++it_B)
229 aVector[1].MarkPathFromNodeToNode(aVector[1].GetRoot(), (*it_B));
231 std::stringstream filepath_actualIteration;
232 filepath_actualIteration << path_folder << name_tree << "__" << iteration << "_idB_" << (*it_B)->GetId() << suffix_vtk;
234 DrawVTKLinesFromTree(aVector[1], filepath_actualIteration.str(), true);
239 // --------------------------------------
240 // --------------------------------------
243 //XXXXXXXXXXXXXXXXXXXXXXXXXXXX
244 //1. PRINT EACH NODE OF EACH TREE AND SAVE IT USING AS NAME ITS ID
247 for(int i = 2; i < aVector[0].GetWeight(); ++i)
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;
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);
260 std::cout << "Tree NULL" << std::endl;
264 for(int i = 2; i < aVector[1].GetWeight(); ++i)
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;
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);
277 std::cout << "Tree NULL" << std::endl;
280 std::cout << "Comparing trees Orkisz-Morales... OK" << std::endl;
284 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
285 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
286 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
289 aVector[0].SubIsomorphism(aVector[1]);
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);
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;
302 aVector[1].ImageReconstruction(imgB);
303 std::cout << "Flag after Reconstruction B" << std::endl;
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);
312 po::notify(vm); // throws on error, so do after help in case
314 catch (std::exception& e)
316 std::cerr << "Unhandled Exception reached the top of main: " << e.what() << ", application will now exit" << std::endl;
317 return ERROR_UNHANDLED_EXCEPTION;
320 catch (std::exception& e)
322 std::cerr << "Unhandled Exception reached the top of main: " << e.what() << ", application will now exit" << std::endl;
323 return ERROR_UNHANDLED_EXCEPTION;
328 // -------------------------------------------------------------------------
329 void Load_cpPlugins( const std::string& plugins, const std::string& ws )
331 myPlugins.LoadConfiguration( plugins );
332 myWorkspace.SetInterface( &myPlugins );
333 std::string err = myWorkspace.LoadWorkspace( ws );
336 std::cerr << "Error: " << err << std::endl;
342 // -------------------------------------------------------------------------
343 void CreateResultDirectory(AirwaysTree tree)
345 boost::filesystem::path dir(tree.GetResultPath());
346 if (boost::filesystem::create_directories(dir))
348 std::cout << "FolderName = " << tree.GetResultPath() << " has been created" << std::endl;
351 std::cout << "FolderName = " << tree.GetResultPath() << " already exists" << std::endl;
354 // -------------------------------------------------------------------------
355 void DrawVTKLinesFromTree(AirwaysTree& tree, const std::string filepath, bool common)
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");
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);
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);
379 // Create the polydata
380 vtkSmartPointer<vtkPolyData> linesPolyData = vtkSmartPointer<vtkPolyData>::New();
382 // Set the points, lines, and colors
383 linesPolyData->SetPoints(pts);
384 linesPolyData->SetLines(lines);
385 linesPolyData->GetCellData()->SetScalars(colors);
388 vtkSmartPointer<vtkPolyDataWriter> writer = vtkSmartPointer<vtkPolyDataWriter>::New();
389 writer->SetFileName(filepath.c_str());
391 #if VTK_MAJOR_VERSION <= 5
392 writer->SetInput(linesPolyData);
394 writer->SetInputData(linesPolyData);
398 //The following code is to test the vtk polydata and visualize it
400 vtkSmartPointer<vtkPolyDataMapper> mapper =
401 vtkSmartPointer<vtkPolyDataMapper>::New();
402 #if VTK_MAJOR_VERSION <= 5
403 mapper->SetInput(linesPolyData);
405 mapper->SetInputData(linesPolyData);
407 vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
408 actor->SetMapper(mapper);
411 vtkSmartPointer<vtkPolyDataMapper> mapperSphere = vtkSmartPointer<
412 vtkPolyDataMapper>::New();
413 mapperSphere->SetInputConnection(sphereSource->GetOutputPort());
415 vtkSmartPointer<vtkActor> actorSphere = vtkSmartPointer<vtkActor>::New();
416 actorSphere->SetMapper(mapperSphere);
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);
429 renderWindow->Render();
430 renderWindowInteractor->Start();*/
433 // -------------------------------------------------------------------------
434 #include <fpa/Base/ImageSkeleton.h>
435 #include <fpa/Image/MinimumSpanningTree.h>
436 #include <cpExtensions/DataStructures/ImageIndexesContainer.h>
437 #include <itkImage.h>
439 template< class TImage, class TVertex >
440 Node* FAVertexToNode( TVertex vertex, TImage* image )
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);
449 AirwaysTree& ConvertFilterToAirwaysTree( )
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;
461 fpa::Base::ImageSkeleton< fpa::Image::MinimumSpanningTree< 3 > >
463 typedef TFASkeleton::TVertex TVertexFA;
464 typedef TFASkeleton::TVertexCmp TVertexCompareFA;
465 typedef cpExtensions::DataStructures::ImageIndexesContainer< 3 > TVerticesFA;
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;
473 std::time_t start,end;
475 std::cout << "Starting conversion " << std::endl;
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 > >( );
485 TBranchesFA* branches = filter->GetBranches( );
486 TMst* mst = filter->GetMinimumSpanningTree();
488 int leoTreeWeigth = endpoints.size()+bifurcations.size()+1;
489 std::cout<< "Creates FA Tree with : "<<leoTreeWeigth << " nodes "<<std::endl;
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.
494 vertexMap[ seed0 ] = FAVertexToNode( seed0, image );
496 auto eIt = endpoints.begin();
497 for( ; eIt != endpoints.end(); ++eIt )
499 vertexMap[*eIt]=FAVertexToNode( *eIt,image );
502 auto biIt = bifurcations.begin();
503 for(; biIt != bifurcations.end(); ++biIt)
505 vertexMap[*biIt]=FAVertexToNode(*biIt,image);
508 //Now navigate branches to make the edges of the tree.
509 auto bIt = branches->Get( ).begin();
510 for(; bIt!=branches->Get( ).end(); ++bIt)
512 auto dVertex = bIt->first;
513 TVertexAirways* destination = vertexMap[dVertex];
514 if( destination == NULL )
516 vertexMap[dVertex]=FAVertexToNode(dVertex,image);
517 destination = vertexMap[dVertex];
519 auto brIt = bIt->second.begin( );
520 for( ; brIt != bIt->second.end( ); ++brIt )
522 auto sVertex = brIt->first;
523 TVertexAirways* source = vertexMap[sVertex];
526 vertexMap[sVertex]=FAVertexToNode(sVertex,image);
527 source = vertexMap[sVertex];
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.
536 auto edgePath = brIt->second->GetVertexList( );
537 for( unsigned int pIt = 0; pIt < edgePath->Size( ); ++pIt )
539 auto cidx = edgePath->GetElement(pIt);
540 itk::ImageBase< 3 >::PointType pnt;
541 image->TransformContinuousIndexToPhysicalPoint(cidx,pnt);
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);
548 destination->SetEdge(edge);
551 AirwaysTree* tree = new AirwaysTree(dynamic_cast< TInputImage* >( image ), NULL, vertexMap[seed0], false);
553 std::cout << "Finished conversion. New AlfTree has weight: "<<tree->GetWeight()<<". Takes "<<(end-start)<<" s." << std::endl;
557 // -------------------------------------------------------------------------
558 AirwaysTree& CreateAirwaysTreeFromSegmentation(Vec3 seed, TImagePointer input_image)
560 std::string err = myWorkspace.Execute( "eb" );
563 std::cerr << "Error: " << err << std::endl;
567 return( ConvertFilterToAirwaysTree( ) );
570 // -------------------------------------------------------------------------
571 vector<TreeInfo> ReadInputFile(const char* filename)
574 std::ifstream infile(filename);
576 std::string folderpath_allResults;
577 vector<TreeInfo> vectorInfo;
579 bool firstLine = false;
581 while (std::getline(infile, line))
583 // First line contains the output folder
584 std::istringstream iss(line);
587 if (!(iss >> folderpath_allResults))
589 std::cout << "no file" << std::endl;
594 // Other lines, not the first one, contain the information for the airways to be created
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))
602 std::cout << "There is no tree information in the file." << std::endl;
608 std::cout << "Point trachea:[" << point_trachea[0] << "," << point_trachea[1] << "," << point_trachea[2] << "]" << std::endl;
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;
622 treeInfo.image = ReadImage<TInputImage>(filepath_airwayImage);
624 Load_cpPlugins( "./plugins.cfg", "./workspace_airwaysappli.wxml" );
626 // Execute first pipeline's part
627 std::stringstream seed_stream;
629 << point_trachea[0] << " "
630 << point_trachea[1] << " "
632 myWorkspace.SetParameter( "FileNames@reader", filepath_airwayImage );
633 myWorkspace.SetParameter( "Text@seed", seed_stream.str( ) );
634 vectorInfo.push_back(treeInfo);
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)
644 std::cout << "printCommonTreeBetweenTwoTrees, edges:" << vector_pair_edges_A_to_B.size() << std::endl;
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");
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);
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)
661 points_common->SetPoint((*it_nodesA)->GetId(),(*it_nodesA)->GetCoords().GetVec3());
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)
668 points_common->SetPoint((*it_nodesB)->GetId()+1000,(*it_nodesB)->GetCoords().GetVec3());
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)
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;
681 edge_A.first->GetPathToNode(edge_A.second, path_A);
684 edge_B.first->GetPathToNode(edge_B.second, path_B);
686 // Set color to be used
687 int numColor = number_pairs % 6;
688 unsigned char color[3];
698 else if(numColor == 2)
704 else if(numColor == 3)
706 color[0] = yellow[0];
707 color[1] = yellow[1];
708 color[2] = yellow[2];
710 else if(numColor == 4)
712 color[0] = purple[0];
713 color[1] = purple[1];
714 color[2] = purple[2];
716 else if(numColor == 5)
726 //if(path_A.size() > 0 && number_pairs < 50)
729 vtkSmartPointer<vtkLine> line = vtkSmartPointer<vtkLine>::New();
730 lines_common->InsertNextCell(path_A.size());
732 for(vec_nodes::iterator it_pathA = path_A.begin(); it_pathA != path_A.end(); ++it_pathA)
734 lines_common->InsertCellPoint((*it_pathA)->GetId());
735 //line->GetPointIds()->SetId( id_point, (*it_pathA)->GetId() );
738 //lines_common->InsertNextCell(line);
739 colors_common->InsertNextTupleValue(color);
743 std::cout << "No path A" << std::endl;
748 //if(path_B.size() > 0 && number_pairs < 50)
751 vtkSmartPointer<vtkLine> line = vtkSmartPointer<vtkLine>::New();
752 lines_common->InsertNextCell(path_B.size());
754 for(vec_nodes::iterator it_pathB = path_B.begin(); it_pathB != path_B.end(); ++it_pathB)
756 lines_common->InsertCellPoint((*it_pathB)->GetId()+1000);
757 //line->GetPointIds()->SetId( id_point, (*it_pathB)->GetId()+1000 );
760 //lines_common->InsertNextCell(line);
761 colors_common->InsertNextTupleValue(color);
765 std::cout << "No path B" << std::endl;
770 // Create the polydata
771 vtkSmartPointer<vtkPolyData> linesPolyData = vtkSmartPointer<vtkPolyData>::New();
773 // Set the points, lines, and colors
774 linesPolyData->SetPoints(points_common);
775 linesPolyData->SetLines(lines_common);
776 linesPolyData->GetCellData()->SetScalars(colors_common);
778 // ------------------------------------------------
779 // Write the vtk file
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;
787 vtkSmartPointer<vtkPolyDataWriter> writer = vtkSmartPointer<vtkPolyDataWriter>::New();
788 writer->SetFileName(filepath_actualIteration.str().c_str());
790 // Set input and write
791 #if VTK_MAJOR_VERSION <= 5
792 writer->SetInput(linesPolyData);
794 writer->SetInputData(linesPolyData);
800 // ***************************************
801 // Save the links in a file
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";
808 std::ofstream file(filepath_evaluation.str().c_str());
812 //file << "Node_" << tree_A.GetImageName() << " " << "Node_" << tree_B.GetImageName()<< "\n";
814 "idLocalN1 idLocalN2 "
815 "Node1x Node1y Node1z "
816 "Node2x Node2y Node2z "
818 "idMatch" << std::endl;
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)
825 std::pair<Node*, Node*> edge_A = (*it_edges).first;
826 std::pair<Node*, Node*> edge_B = (*it_edges).second;
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";
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";
849 // ***************************************
854 std::cout << "printCommonTreeBetweenTwoTrees ... OK" << std::endl;
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)
860 std::cout << "Printing matching result ... " << std::endl;
862 // Variables and types
863 typedef std::map< unsigned int, vec_nodes > map_id_node;
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";
869 std::ofstream file(filepath_evaluation.str().c_str());
873 file << "TreeName idLocalN1 idLocalN2 "
874 "Node1x Node1y Node1z "
875 "Node2x Node2y Node2z "
878 "typeMatch_match_1_nonmatch_0 "
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)
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();
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( ) )
898 // Get the correspoding matching nodes and print them
899 vec_nodes nodes_B = (*it_a2b).second;
901 vec_nodes::iterator it_nodes_b = nodes_B.begin( );
902 for( ; it_nodes_b != nodes_B.end( ); ++it_nodes_b)
904 Vec3 coords_b = (*it_nodes_b)->GetCoords();
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] << " "
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" << " "
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)
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();
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( ) )
945 // Get the correspoding matching nodes and print them
946 vec_nodes nodes_A = (*it_b2a).second;
948 vec_nodes::iterator it_nodes_a = nodes_A.begin( );
949 for( ; it_nodes_a != nodes_A.end( ); ++it_nodes_a)
951 Vec3 coords_a = (*it_nodes_a)->GetCoords();
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] << " "
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" << " "
977 file.close(); // Close the file
979 std::cout << "Printing matching result DONE" << std::endl;
982 // -------------------------------------------------------------------------
983 void createLinesAndPointsForVTK(const Node* node, vtkSmartPointer<vtkPoints>& pts,vtkSmartPointer<vtkCellArray>& lines, vtkSmartPointer<vtkUnsignedCharArray>& colors, bool common, bool isRoot)
985 // Insert the actual point/node
986 //vtkIdType id_father = idNonRoot;
988 //id_father = pts->InsertNextPoint(node->GetCoords().GetVec3());
989 vtkIdType id_father = node->GetId();
991 pts->SetPoint(id_father,node->GetCoords().GetVec3());
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)
997 if (!(*it_child)->IsMarked() && common)
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());
1004 // Set color to be used
1005 int numColor = (*it_child)->GetLevel() % 4;
1007 colors->InsertNextTupleValue(green);
1008 else if(numColor == 1)
1009 colors->InsertNextTupleValue(red);
1010 else if(numColor == 2)
1011 colors->InsertNextTupleValue(red);
1013 colors->InsertNextTupleValue(red);
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);
1021 createLinesAndPointsForVTK(*it_child, pts, lines, colors, common, false);