+ // Create polydata
+ vtkSmartPointer< vtkPoints > points =
+ vtkSmartPointer< vtkPoints >::New( );
+ vtkSmartPointer< vtkCellArray > cells =
+ vtkSmartPointer< vtkCellArray >::New( );
+ vtkSmartPointer< vtkFloatArray > scalars =
+ vtkSmartPointer< vtkFloatArray >::New( );
+
+ TDijkstra::TMarkImage::Pointer new_marks = TDijkstra::TMarkImage::New( );
+ new_marks->SetLargestPossibleRegion( input_image->GetLargestPossibleRegion( ) );
+ new_marks->SetRequestedRegion( input_image->GetRequestedRegion( ) );
+ new_marks->SetBufferedRegion( input_image->GetBufferedRegion( ) );
+ new_marks->SetOrigin( input_image->GetOrigin( ) );
+ new_marks->SetSpacing( input_image->GetSpacing( ) );
+ new_marks->SetDirection( input_image->GetDirection( ) );
+ new_marks->Allocate( );
+ new_marks->FillBuffer( 0 );
+
+ const TDijkstra::TMarkImage* marks = paths->GetOutputMarkImage( );
+ TDijkstra::TMark max_mark = paths->GetNumberOfBranches( );
+ const TDijkstra::TVertices& endpoints = paths->GetEndPoints( );
+ const TDijkstra::TVertices& bifurcations = paths->GetBifurcationPoints( );
+ const TDijkstra::TTree& tree = paths->GetFullTree( );
+ const TDijkstra::TTree& reduced_tree = paths->GetReducedTree( );
+ TDijkstra::TVertices::const_iterator epIt = endpoints.begin( );
+
+ TBranches branches;
+ TImage::PointType ori = input_image->GetOrigin( );
+
+ TDijkstra::TTree::const_iterator rIt = reduced_tree.begin( );
+ for( ; rIt != reduced_tree.end( ); ++rIt )
+ {
+ TDijkstra::TTree::const_iterator tIt = tree.find( rIt->first );
+
+ // Prepare branch "a la ImageJ"
+ TBranch branch;
+ TImage::PointType p1, p2;
+ input_image->TransformIndexToPhysicalPoint( rIt->second.first, p1 );
+ input_image->TransformIndexToPhysicalPoint( rIt->first, p2 );
+ branch.V1 = p1 - ori;
+ branch.V2 = p2 - ori;
+ branch.Length = double( 0 );
+
+ double pd =
+ ( double( tIt->second.second ) - double( 1 ) ) /
+ ( double( max_mark ) - double( 1 ) );
+ bool start = true;
+ TImage::PointType prev;
+ do
+ {
+ TDijkstra::TVertex idx = tIt->first;
+ new_marks->SetPixel( idx, 255 );
+ TImage::PointType pnt;
+ input_image->TransformIndexToPhysicalPoint( idx, pnt );
+ points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
+ scalars->InsertNextTuple1( pd );
+ if( !start )
+ {
+ cells->InsertNextCell( 2 );
+ cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
+ cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
+ branch.Length += prev.EuclideanDistanceTo( pnt );
+
+ } // fi
+ start = false;
+ prev = pnt;
+ tIt = tree.find( tIt->second.first );
+
+ } while( tIt->first != rIt->second.first );
+
+ // Last point
+ TDijkstra::TVertex idx = tIt->first;
+ new_marks->SetPixel( idx, 255 );
+ TImage::PointType pnt;
+ input_image->TransformIndexToPhysicalPoint( idx, pnt );
+ points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
+ scalars->InsertNextTuple1( pd );
+ if( !start )
+ {
+ cells->InsertNextCell( 2 );
+ cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
+ cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
+ branch.Length += prev.EuclideanDistanceTo( pnt );
+
+ } // fi
+ branches.insert( branch );
+
+ } // rof
+
+ std::ofstream bout( output_imagej_fn.c_str( ) );
+ if( !bout )
+ {
+ std::cerr << "Error: could not open " << output_imagej_fn << std::endl;
+ return( 1 );
+
+ } // fi
+ int i = 1;
+ for( TBranches::const_iterator bIt = branches.begin( ); bIt != branches.end( ); ++bIt, ++i )
+ {
+ bout
+ << std::setiosflags( std::ios::fixed) << std::setprecision( 3 )
+ << i << "\t1.000\t"
+ << bIt->Length << "\t"
+ << bIt->V1[ 0 ] << "\t"
+ << bIt->V1[ 1 ] << "\t"
+ << bIt->V1[ 2 ] << "\t"
+ << bIt->V2[ 0 ] << "\t"
+ << bIt->V2[ 1 ] << "\t"
+ << bIt->V2[ 2 ] << "\t"
+ << ( bIt->V2 - bIt->V1 ).GetNorm( )
+ << std::endl;
+
+ } // rof
+ bout.close( );
+
+ vtkSmartPointer< vtkPolyData > vtk_tree =
+ vtkSmartPointer< vtkPolyData >::New( );
+ vtk_tree->SetPoints( points );
+ vtk_tree->SetLines( cells );
+ vtk_tree->GetPointData( )->SetScalars( scalars );
+
+ view.AddPolyData( vtk_tree );
+ view.Render( );
+ if( visual_debug )
+ view.Start( );
+
+ vtkSmartPointer< vtkPolyDataWriter > writer =
+ vtkSmartPointer< vtkPolyDataWriter >::New( );
+ writer->SetInputData( vtk_tree );
+ writer->SetFileName( output_vtk_fn.c_str( ) );
+ writer->Update( );
+
+ itk::ImageFileWriter< TDijkstra::TMarkImage >::Pointer marks_writer =
+ itk::ImageFileWriter< TDijkstra::TMarkImage >::New( );
+ marks_writer->SetInput( new_marks );
+ marks_writer->SetFileName( output_image_fn );
+ marks_writer->Update( );