From b4ed6ddfa7e90e892f07cad9a760961bd4e84e6b Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Thu, 9 Apr 2015 18:48:36 -0500 Subject: [PATCH] I/O classes added --- appli/examples/CMakeLists.txt | 9 + ...ample_Image_Dijkstra_EndPointDetection.cxx | 265 +++++---- ..._Dijkstra_EndPointDetection_WithoutVTK.cxx | 249 ++++++++ lib/CMakeLists.txt | 3 + lib/fpa/Base/Dijkstra.h | 9 +- lib/fpa/Base/MatrixValuesContainer.h | 80 +++ lib/fpa/Base/MinimumSpanningTree.h | 2 +- lib/fpa/Base/MinimumSpanningTree.hxx | 1 + lib/fpa/Base/UniqueValuesContainer.h | 72 +++ lib/fpa/IO/MatrixValuesContainerWriter.h | 61 ++ lib/fpa/IO/MatrixValuesContainerWriter.hxx | 112 ++++ lib/fpa/IO/MinimumSpanningTreeWriter.h | 61 ++ lib/fpa/IO/MinimumSpanningTreeWriter.hxx | 118 ++++ lib/fpa/IO/UniqueValuesContainerWriter.h | 61 ++ lib/fpa/IO/UniqueValuesContainerWriter.hxx | 105 ++++ lib/fpa/Image/Algorithm.h | 9 +- lib/fpa/Image/Dijkstra.h | 11 +- lib/fpa/Image/DijkstraWithEndPointDetection.h | 53 +- .../Image/DijkstraWithEndPointDetection.hxx | 559 +++++++++++------- lib/fpa/VTK/ImageMPR.cxx | 5 +- 20 files changed, 1486 insertions(+), 359 deletions(-) create mode 100644 appli/examples/example_Image_Dijkstra_EndPointDetection_WithoutVTK.cxx create mode 100644 lib/fpa/Base/MatrixValuesContainer.h create mode 100644 lib/fpa/Base/UniqueValuesContainer.h create mode 100644 lib/fpa/IO/MatrixValuesContainerWriter.h create mode 100644 lib/fpa/IO/MatrixValuesContainerWriter.hxx create mode 100644 lib/fpa/IO/MinimumSpanningTreeWriter.h create mode 100644 lib/fpa/IO/MinimumSpanningTreeWriter.hxx create mode 100644 lib/fpa/IO/UniqueValuesContainerWriter.h create mode 100644 lib/fpa/IO/UniqueValuesContainerWriter.hxx diff --git a/appli/examples/CMakeLists.txt b/appli/examples/CMakeLists.txt index 61768cb..fd255c1 100644 --- a/appli/examples/CMakeLists.txt +++ b/appli/examples/CMakeLists.txt @@ -1,3 +1,12 @@ +SET( + SIMPLE_EXAMPLES + example_Image_Dijkstra_EndPointDetection_WithoutVTK + ) +FOREACH(EX ${SIMPLE_EXAMPLES}) + ADD_EXECUTABLE(${EX} ${EX}.cxx) + TARGET_LINK_LIBRARIES(${EX} FrontAlgorithms) +ENDFOREACH(EX) + IF(USE_VTK) SET( SIMPLE_VTK_EXAMPLES diff --git a/appli/examples/example_Image_Dijkstra_EndPointDetection.cxx b/appli/examples/example_Image_Dijkstra_EndPointDetection.cxx index b2d3e7c..a334de3 100644 --- a/appli/examples/example_Image_Dijkstra_EndPointDetection.cxx +++ b/appli/examples/example_Image_Dijkstra_EndPointDetection.cxx @@ -20,7 +20,12 @@ #include #include +#include +#include +#include + #include +#include #include #include #include @@ -50,11 +55,18 @@ void DistanceMap( // ------------------------------------------------------------------------- int main( int argc, char* argv[] ) { - if( argc < 5 ) + if( argc < 9 ) { std::cerr - << "Usage: " << argv[ 0 ] - << " input_image distancemap output_costmap output_labels" + << "Usage: " << argv[ 0 ] << std::endl + << " input_image" << std::endl + << " output_distancemap" << std::endl + << " output_costmap" << std::endl + << " output_labels" << std::endl + << " output_minimum_spanning_tree" << std::endl + << " output_endpoints" << std::endl + << " output_bifurcations" << std::endl + << " output_branches" << std::endl << std::endl; return( 1 ); @@ -63,6 +75,10 @@ int main( int argc, char* argv[] ) std::string distancemap_fn = argv[ 2 ]; std::string output_costmap_fn = argv[ 3 ]; std::string output_labels_fn = argv[ 4 ]; + std::string mst_output_fn = argv[ 5 ]; + std::string endpoints_output_fn = argv[ 6 ]; + std::string bifurcations_output_fn = argv[ 7 ]; + std::string branches_output_fn = argv[ 8 ]; // Read image TImage::Pointer input_image; @@ -147,142 +163,183 @@ int main( int argc, char* argv[] ) << std::difftime( end, start ) << " s." << std::endl; - // Minimum spanning tree - const TFilter::TMinimumSpanningTree* mst = - filter->GetMinimumSpanningTree( ); + // Outputs + const TFilter::TOutputImage* accumulated_costs = filter->GetOutput( ); + const TFilter::TLabelImage* labeled_image = filter->GetLabelImage( ); + const TFilter::TMinimumSpanningTree* mst = filter->GetMinimumSpanningTree( ); + const TFilter::TUniqueVertices* endpoints = filter->GetEndPoints( ); + const TFilter::TUniqueVertices* bifurcations = filter->GetBifurcations( ); + const TFilter::TBranches* branches = filter->GetBranches( ); + unsigned long nBranches = filter->GetNumberOfBranches( ); - // Build branches from endpoints to seed + // Save outputs + SaveImage( dmap.GetPointer( ), distancemap_fn ); + SaveImage( accumulated_costs, output_costmap_fn ); + SaveImage( labeled_image, output_labels_fn ); + + fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::Pointer + mst_writer = + fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::New( ); + mst_writer->SetInput( mst ); + mst_writer->SetFileName( mst_output_fn ); + mst_writer->Update( ); + + fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer + endpoints_writer = + fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( ); + endpoints_writer->SetInput( endpoints ); + endpoints_writer->SetFileName( endpoints_output_fn ); + endpoints_writer->Update( ); + + fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer + bifurcations_writer = + fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( ); + bifurcations_writer->SetInput( bifurcations ); + bifurcations_writer->SetFileName( bifurcations_output_fn ); + bifurcations_writer->Update( ); + + fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::Pointer + branches_writer = + fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::New( ); + branches_writer->SetInput( branches ); + branches_writer->SetFileName( branches_output_fn ); + branches_writer->Update( ); + + // Show endpoints vtkSmartPointer< vtkPoints > endpoints_points = vtkSmartPointer< vtkPoints >::New( ); vtkSmartPointer< vtkCellArray > endpoints_cells = vtkSmartPointer< vtkCellArray >::New( ); - vtkSmartPointer< vtkCellArray > endpoints_vertices = - vtkSmartPointer< vtkCellArray >::New( ); - - const TFilter::TUniqueVertices& endpoints = filter->GetEndPoints( ); - TFilter::TUniqueVertices::const_iterator eIt = endpoints.begin( ); - for( ; eIt != endpoints.end( ); ++eIt ) + for( + TFilter::TUniqueVertices::ConstIterator epIt = endpoints->Begin( ); + epIt != endpoints->End( ); + ++epIt + ) { - TFilter::TVertices path; - mst->GetPath( path, *eIt, filter->GetSeed( 0 ) ); - - std::cout << *eIt << std::endl; - TFilter::TVertices::const_iterator pIt = path.begin( ); - for( ; pIt != path.end( ); ++pIt ) - { - TImage::PointType pnt; - input_image->TransformIndexToPhysicalPoint( *pIt, pnt ); - endpoints_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); - - if( pIt != path.begin( ) ) - { - endpoints_cells->InsertNextCell( 2 ); - endpoints_cells->InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 2 ); - endpoints_cells->InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 1 ); - } - else - { - std::cout << "ok" << std::endl; - endpoints_vertices->InsertNextCell( 1 ); - endpoints_vertices->InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 1 ); - - } // fi - - - - - } // rof + TImage::PointType pnt; + input_image->TransformIndexToPhysicalPoint( *epIt, pnt ); + endpoints_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); + endpoints_cells->InsertNextCell( 1 ); + endpoints_cells-> + InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 1 ); } // rof - std::cout << endpoints.size( ) << std::endl; - vtkSmartPointer< vtkPolyData > endpoints_polydata = vtkSmartPointer< vtkPolyData >::New( ); endpoints_polydata->SetPoints( endpoints_points ); - endpoints_polydata->SetLines( endpoints_cells ); - endpoints_polydata->SetVerts( endpoints_vertices ); - - - - - - - - + endpoints_polydata->SetVerts( endpoints_cells ); + view.AddPolyData( endpoints_polydata, 0, 0, 1, 1 ); - - - // Bifurcations + // Show bifurcations vtkSmartPointer< vtkPoints > bifurcations_points = vtkSmartPointer< vtkPoints >::New( ); - vtkSmartPointer< vtkCellArray > bifurcations_vertices = + vtkSmartPointer< vtkCellArray > bifurcations_cells = vtkSmartPointer< vtkCellArray >::New( ); - - const TFilter::TUniqueVertices& bifurcations = filter->GetBifurcationPoints( ); - TFilter::TUniqueVertices::const_iterator bIt = bifurcations.begin( ); - for( ; bIt != bifurcations.end( ); ++bIt ) + for( + TFilter::TUniqueVertices::ConstIterator bfIt = bifurcations->Begin( ); + bfIt != bifurcations->End( ); + ++bfIt + ) { - std::cout << *bIt << std::endl; - TImage::PointType pnt; - input_image->TransformIndexToPhysicalPoint( *bIt, pnt ); + input_image->TransformIndexToPhysicalPoint( *bfIt, pnt ); bifurcations_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); - bifurcations_vertices->InsertNextCell( 1 ); - bifurcations_vertices->InsertCellPoint( bifurcations_points->GetNumberOfPoints( ) - 1 ); + bifurcations_cells->InsertNextCell( 1 ); + bifurcations_cells-> + InsertCellPoint( bifurcations_points->GetNumberOfPoints( ) - 1 ); } // rof - vtkSmartPointer< vtkPolyData > bifurcations_polydata = vtkSmartPointer< vtkPolyData >::New( ); bifurcations_polydata->SetPoints( bifurcations_points ); - bifurcations_polydata->SetVerts( bifurcations_vertices ); - - - + bifurcations_polydata->SetVerts( bifurcations_cells ); + view.AddPolyData( bifurcations_polydata, 0, 1, 0, 1 ); + // Show branches (simple and detailed) + vtkSmartPointer< vtkPoints > simple_branches_points = + vtkSmartPointer< vtkPoints >::New( ); + vtkSmartPointer< vtkCellArray > simple_branches_cells = + vtkSmartPointer< vtkCellArray >::New( ); - // Build branches from endpoints to seed - vtkSmartPointer< vtkPoints > branches_points = + vtkSmartPointer< vtkPoints > detailed_branches_points = vtkSmartPointer< vtkPoints >::New( ); - vtkSmartPointer< vtkCellArray > branches_cells = + vtkSmartPointer< vtkCellArray > detailed_branches_cells = vtkSmartPointer< vtkCellArray >::New( ); + vtkSmartPointer< vtkFloatArray > detailed_branches_scalars = + vtkSmartPointer< vtkFloatArray >::New( ); - const TFilter::TBranches& branches = filter->GetBranches( ); - TFilter::TBranches::const_iterator brIt = branches.begin( ); - for( ; brIt != branches.end( ); ++brIt ) + TFilter::TBranches::ConstIterator brIt = branches->Begin( ); + for( ; brIt != branches->End( ); ++brIt ) { - TFilter::TBranch::const_iterator brsIt = brIt->second.begin( ); - for( ; brsIt != brIt->second.end( ); ++brsIt ) + // Branch's first point + TImage::PointType first_point; + input_image->TransformIndexToPhysicalPoint( brIt->first, first_point ); + unsigned long first_id = simple_branches_points->GetNumberOfPoints( ); + simple_branches_points->InsertNextPoint( + first_point[ 0 ], first_point[ 1 ], first_point[ 2 ] + ); + + TFilter::TBranches::ConstRowIterator brRowIt = branches->Begin( brIt ); + for( ; brRowIt != branches->End( brIt ); ++brRowIt ) { - TImage::PointType pnt0, pnt1; - input_image->TransformIndexToPhysicalPoint( brIt->first, pnt0 ); - input_image->TransformIndexToPhysicalPoint( brsIt->first, pnt1 ); - branches_points->InsertNextPoint( pnt0[ 0 ], pnt0[ 1 ], pnt0[ 2 ] ); - branches_points->InsertNextPoint( pnt1[ 0 ], pnt1[ 1 ], pnt1[ 2 ] ); - - branches_cells->InsertNextCell( 2 ); - branches_cells->InsertCellPoint( branches_points->GetNumberOfPoints( ) - 2 ); - branches_cells->InsertCellPoint( branches_points->GetNumberOfPoints( ) - 1 ); + // Branch's second point + TImage::PointType second_point; + input_image-> + TransformIndexToPhysicalPoint( brRowIt->first, second_point ); + unsigned long second_id = simple_branches_points->GetNumberOfPoints( ); + simple_branches_points->InsertNextPoint( + second_point[ 0 ], second_point[ 1 ], second_point[ 2 ] + ); + simple_branches_cells->InsertNextCell( 2 ); + simple_branches_cells->InsertCellPoint( first_id ); + simple_branches_cells->InsertCellPoint( second_id ); + + // Detailed path + double pathId = double( brRowIt->second - 1 ) / double( nBranches - 1 ); + TFilter::TVertices path; + mst->GetPath( path, brIt->first, brRowIt->first ); + TFilter::TVertices::const_iterator pIt = path.begin( ); + for( ; pIt != path.end( ); ++pIt ) + { + TImage::PointType path_point; + input_image->TransformIndexToPhysicalPoint( *pIt, path_point ); + detailed_branches_points->InsertNextPoint( + path_point[ 0 ], path_point[ 1 ], path_point[ 2 ] + ); + detailed_branches_scalars->InsertNextTuple1( pathId ); + if( pIt != path.begin( ) ) + { + unsigned long nPoints = + detailed_branches_points->GetNumberOfPoints( ); + detailed_branches_cells->InsertNextCell( 2 ); + detailed_branches_cells->InsertCellPoint( nPoints - 2 ); + detailed_branches_cells->InsertCellPoint( nPoints - 1 ); + + } // fi + + } // rof } // rof } // rof + vtkSmartPointer< vtkPolyData > simple_branches_polydata = + vtkSmartPointer< vtkPolyData >::New( ); + simple_branches_polydata->SetPoints( simple_branches_points ); + simple_branches_polydata->SetLines( simple_branches_cells ); + view.AddPolyData( simple_branches_polydata, 1, 0, 1, 1 ); - vtkSmartPointer< vtkPolyData > branches_polydata = + vtkSmartPointer< vtkPolyData > detailed_branches_polydata = vtkSmartPointer< vtkPolyData >::New( ); - branches_polydata->SetPoints( branches_points ); - branches_polydata->SetLines( branches_cells ); + detailed_branches_polydata->SetPoints( detailed_branches_points ); + detailed_branches_polydata->SetLines( detailed_branches_cells ); + detailed_branches_polydata-> + GetPointData( )->SetScalars( detailed_branches_scalars ); + view.AddPolyData( detailed_branches_polydata, 1 ); - view.AddPolyData( endpoints_polydata, 1, 0, 0, 1 ); - view.AddPolyData( bifurcations_polydata, 0, 0, 1, 1 ); - view.AddPolyData( branches_polydata, 0, 1, 0, 1 ); + // Let some more interaction view.Render( ); view.Start( ); - // Save final total cost map - SaveImage( filter->GetOutput( ), output_costmap_fn ); - SaveImage( dmap.GetPointer( ), distancemap_fn ); - SaveImage( filter->GetLabelImage( ), output_labels_fn ); return( 0 ); } @@ -362,13 +419,3 @@ void DistanceMap( } // eof - $RCSfile$ - - - - - - - - - - diff --git a/appli/examples/example_Image_Dijkstra_EndPointDetection_WithoutVTK.cxx b/appli/examples/example_Image_Dijkstra_EndPointDetection_WithoutVTK.cxx new file mode 100644 index 0000000..1a8b476 --- /dev/null +++ b/appli/examples/example_Image_Dijkstra_EndPointDetection_WithoutVTK.cxx @@ -0,0 +1,249 @@ +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include + +#include + +#include +#include + +#include +#include +#include + +// ------------------------------------------------------------------------- +const unsigned int Dim = 3; +typedef unsigned char TPixel; +typedef float TScalar; + +typedef itk::Image< TPixel, Dim > TImage; +typedef itk::Image< TScalar, Dim > TScalarImage; + +// ------------------------------------------------------------------------- +template< class I > +void ReadImage( typename I::Pointer& image, const std::string& filename ); + +template< class I > +void SaveImage( const I* image, const std::string& filename ); + +template< class I, class O > +void DistanceMap( + const typename I::Pointer& input, typename O::Pointer& output + ); + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + if( argc < 12 ) + { + std::cerr + << "Usage: " << argv[ 0 ] << std::endl + << " input_image" << std::endl + << " seed_x seed_y seed_z" << std::endl + << " output_distancemap" << std::endl + << " output_costmap" << std::endl + << " output_labels" << std::endl + << " output_minimum_spanning_tree" << std::endl + << " output_endpoints" << std::endl + << " output_bifurcations" << std::endl + << " output_branches" << std::endl + << std::endl; + return( 1 ); + + } // fi + std::string input_image_fn = argv[ 1 ]; + double seed_x = std::atof( argv[ 2 ] ); + double seed_y = std::atof( argv[ 3 ] ); + double seed_z = std::atof( argv[ 4 ] ); + std::string distancemap_fn = argv[ 5 ]; + std::string output_costmap_fn = argv[ 6 ]; + std::string output_labels_fn = argv[ 7 ]; + std::string mst_output_fn = argv[ 8 ]; + std::string endpoints_output_fn = argv[ 9 ]; + std::string bifurcations_output_fn = argv[ 10 ]; + std::string branches_output_fn = argv[ 11 ]; + + // Read image + TImage::Pointer input_image; + try + { + ReadImage< TImage >( input_image, input_image_fn ); + } + catch( itk::ExceptionObject& err ) + { + std::cerr + << "Error caught while reading \"" + << input_image_fn << "\": " << err + << std::endl; + return( 1 ); + + } // yrt + + // Get seed + double p[ 3 ]; + TImage::PointType pnt; + TImage::IndexType seed; + pnt[ 0 ] = seed_x; + pnt[ 1 ] = seed_y; + pnt[ 2 ] = seed_z; + input_image->TransformPhysicalPointToIndex( pnt, seed ); + + // Compute squared distance map + TScalarImage::Pointer dmap; + DistanceMap< TImage, TScalarImage >( input_image, dmap ); + + // Prepare cost conversion function + typedef fpa::Base::Functors::InvertCostFunction< TScalar > TFunction; + TFunction::Pointer function = TFunction::New( ); + + // Prepare Dijkstra filter + typedef fpa::Image::DijkstraWithEndPointDetection< TScalarImage, TScalarImage > TFilter; + TFilter::Pointer filter = TFilter::New( ); + filter->SetInput( dmap ); + filter->SetConversionFunction( function ); + filter->SetNeighborhoodOrder( 2 ); + filter->StopAtOneFrontOff( ); + filter->AddSeed( seed, TScalar( 0 ) ); + + // Go! + std::time_t start, end; + std::time( &start ); + filter->Update( ); + std::time( &end ); + std::cout + << "Extraction time = " + << std::difftime( end, start ) + << " s." << std::endl; + + // Outputs + const TFilter::TOutputImage* accumulated_costs = filter->GetOutput( ); + const TFilter::TLabelImage* labeled_image = filter->GetLabelImage( ); + const TFilter::TMinimumSpanningTree* mst = filter->GetMinimumSpanningTree( ); + const TFilter::TUniqueVertices* endpoints = filter->GetEndPoints( ); + const TFilter::TUniqueVertices* bifurcations = filter->GetBifurcations( ); + const TFilter::TBranches* branches = filter->GetBranches( ); + unsigned long nBranches = filter->GetNumberOfBranches( ); + + // Save outputs + SaveImage( dmap.GetPointer( ), distancemap_fn ); + SaveImage( accumulated_costs, output_costmap_fn ); + SaveImage( labeled_image, output_labels_fn ); + + fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::Pointer + mst_writer = + fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::New( ); + mst_writer->SetInput( mst ); + mst_writer->SetFileName( mst_output_fn ); + mst_writer->Update( ); + + fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer + endpoints_writer = + fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( ); + endpoints_writer->SetInput( endpoints ); + endpoints_writer->SetFileName( endpoints_output_fn ); + endpoints_writer->Update( ); + + fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer + bifurcations_writer = + fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( ); + bifurcations_writer->SetInput( bifurcations ); + bifurcations_writer->SetFileName( bifurcations_output_fn ); + bifurcations_writer->Update( ); + + fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::Pointer + branches_writer = + fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::New( ); + branches_writer->SetInput( branches ); + branches_writer->SetFileName( branches_output_fn ); + branches_writer->Update( ); + + return( 0 ); +} + +// ------------------------------------------------------------------------- +template< class I > +void ReadImage( typename I::Pointer& image, const std::string& filename ) +{ + typename itk::ImageFileReader< I >::Pointer reader = + itk::ImageFileReader< I >::New( ); + reader->SetFileName( filename ); + reader->Update( ); + + typename itk::MinimumMaximumImageCalculator< I >::Pointer minmax = + itk::MinimumMaximumImageCalculator< I >::New( ); + minmax->SetImage( reader->GetOutput( ) ); + minmax->Compute( ); + double vmin = double( minmax->GetMinimum( ) ); + double vmax = double( minmax->GetMaximum( ) ); + + typename itk::ShiftScaleImageFilter< I, I >::Pointer shift = + itk::ShiftScaleImageFilter< I, I >::New( ); + shift->SetInput( reader->GetOutput( ) ); + shift->SetScale( vmax - vmin ); + shift->SetShift( vmin / ( vmax - vmin ) ); + shift->Update( ); + + image = shift->GetOutput( ); + image->DisconnectPipeline( ); +} + +// ------------------------------------------------------------------------- +template< class I > +void SaveImage( const I* image, const std::string& filename ) +{ + typename itk::ImageFileWriter< I >::Pointer writer = + itk::ImageFileWriter< I >::New( ); + writer->SetInput( image ); + writer->SetFileName( filename ); + try + { + writer->Update( ); + } + catch( itk::ExceptionObject& err ) + { + std::cerr + << "Error saving \"" << filename << "\": " << err + << std::endl; + + } // yrt +} + +// ------------------------------------------------------------------------- +template< class I, class O > +void DistanceMap( + const typename I::Pointer& input, typename O::Pointer& output + ) +{ + typename itk::SignedMaurerDistanceMapImageFilter< I, O >::Pointer dmap = + itk::SignedMaurerDistanceMapImageFilter< I, O >::New( ); + dmap->SetInput( input ); + dmap->SetBackgroundValue( ( typename I::PixelType )( 0 ) ); + dmap->InsideIsPositiveOn( ); + dmap->SquaredDistanceOn( ); + dmap->UseImageSpacingOn( ); + + std::time_t start, end; + std::time( &start ); + dmap->Update( ); + std::time( &end ); + std::cout + << "Distance map time = " + << std::difftime( end, start ) + << " s." << std::endl; + + output = dmap->GetOutput( ); + output->DisconnectPipeline( ); +} + +// eof - $RCSfile$ + diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index 28c4df3..744f619 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -11,10 +11,12 @@ CONFIGURE_FILE( FILE(GLOB ${LIB_NAME}_HEADERS "fpa/*.h" "fpa/*.hxx") FILE(GLOB ${LIB_NAME}_BASE_HEADERS "fpa/Base/*.h" "fpa/Base/*.hxx") +FILE(GLOB ${LIB_NAME}_IO_HEADERS "fpa/IO/*.h" "fpa/IO/*.hxx") FILE(GLOB ${LIB_NAME}_IMAGE_HEADERS "fpa/Image/*.h" "fpa/Image/*.hxx") FILE(GLOB ${LIB_NAME}_SOURCES "fpa/*.cxx") FILE(GLOB ${LIB_NAME}_BASE_SOURCES "fpa/Base/*.cxx") +FILE(GLOB ${LIB_NAME}_IO_SOURCES "fpa/IO/*.cxx") FILE(GLOB ${LIB_NAME}_IMAGE_SOURCES "fpa/Image/*.cxx") IF(USE_VTK) @@ -27,6 +29,7 @@ SET( ${PROJECT_BINARY_DIR}/lib/fpa/Common.cxx ${${LIB_NAME}_SOURCES} ${${LIB_NAME}_BASE_SOURCES} + ${${LIB_NAME}_IO_SOURCES} ${${LIB_NAME}_IMAGE_SOURCES} ${${LIB_NAME}_VTK_SOURCES} ) diff --git a/lib/fpa/Base/Dijkstra.h b/lib/fpa/Base/Dijkstra.h index cfa06b3..602608e 100644 --- a/lib/fpa/Base/Dijkstra.h +++ b/lib/fpa/Base/Dijkstra.h @@ -29,10 +29,11 @@ namespace fpa typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; - typedef typename Superclass::TVertex TVertex; - typedef typename Superclass::TValue TValue; - typedef typename Superclass::TResult TResult; - typedef typename Superclass::TVertexCompare TVertexCompare; + typedef typename Superclass::TVertex TVertex; + typedef typename Superclass::TValue TValue; + typedef typename Superclass::TResult TResult; + typedef typename Superclass::TVertexCompare TVertexCompare; + typedef typename Superclass::TMinimumSpanningTree TMinimumSpanningTree; typedef typename Superclass::TStartEvent TStartEvent; typedef typename Superclass::TStartLoopEvent TStartLoopEvent; diff --git a/lib/fpa/Base/MatrixValuesContainer.h b/lib/fpa/Base/MatrixValuesContainer.h new file mode 100644 index 0000000..c5e4e41 --- /dev/null +++ b/lib/fpa/Base/MatrixValuesContainer.h @@ -0,0 +1,80 @@ +#ifndef __FPA__BASE__MATRIXVALUESCONTAINER__H__ +#define __FPA__BASE__MATRIXVALUESCONTAINER__H__ + +#include +#include +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< class I, class V, class VI > + class MatrixValuesContainer + : public itk::SimpleDataObjectDecorator< std::map< I, std::map< I, V, VI >, VI > > + { + public: + typedef std::map< I, V, VI > TDecoratedRow; + typedef std::map< I, TDecoratedRow, VI > TDecorated; + typedef MatrixValuesContainer Self; + typedef itk::SimpleDataObjectDecorator< TDecorated > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef I TIndex; + typedef V TValue; + typedef VI TIndexCompare; + + typedef typename TDecoratedRow::iterator RowIterator; + typedef typename TDecoratedRow::const_iterator ConstRowIterator; + typedef typename TDecorated::iterator Iterator; + typedef typename TDecorated::const_iterator ConstIterator; + + public: + itkNewMacro( Self ); + itkTypeMacro( MatrixValuesContainer, itkSimpleDataObjectDecorator ); + + public: + void SetValue( const I& a, const I& b, const V& v ) + { this->Get( )[ a ][ b ] = v; this->Modified( ); } + void Clear( ) + { this->Get( ).clear( ); this->Modified( ); } + Iterator Begin( ) + { return( this->Get( ).begin( ) ); } + Iterator End( ) + { return( this->Get( ).end( ) ); } + ConstIterator Begin( ) const + { return( this->Get( ).begin( ) ); } + ConstIterator End( ) const + { return( this->Get( ).end( ) ); } + RowIterator Begin( const Iterator& i ) + { return( i->second.begin( ) ); } + RowIterator End( const Iterator& i ) + { return( i->second.end( ) ); } + ConstRowIterator Begin( const ConstIterator& i ) const + { return( i->second.begin( ) ); } + ConstRowIterator End( const ConstIterator& i ) const + { return( i->second.end( ) ); } + + protected: + MatrixValuesContainer( ) + : Superclass( ) + { } + virtual ~MatrixValuesContainer( ) + { } + + private: + // Purposely not implemented + MatrixValuesContainer( const Self& other ); + Self& operator=( const Self& other ); + }; + + } // ecapseman + +} // ecapseman + +#endif // __FPA__BASE__MATRIXVALUESCONTAINER__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/MinimumSpanningTree.h b/lib/fpa/Base/MinimumSpanningTree.h index 9adcd6d..be6dae4 100644 --- a/lib/fpa/Base/MinimumSpanningTree.h +++ b/lib/fpa/Base/MinimumSpanningTree.h @@ -35,7 +35,7 @@ namespace fpa public: itkNewMacro( Self ); - itkTypeMacro( MinimumSpanningTree, B ); + itkTypeMacro( MinimumSpanningTree, itkSimpleDataObjectDecorator ); itkGetConstMacro( Collisions, TCollisions ); diff --git a/lib/fpa/Base/MinimumSpanningTree.hxx b/lib/fpa/Base/MinimumSpanningTree.hxx index 6c9c790..d932247 100644 --- a/lib/fpa/Base/MinimumSpanningTree.hxx +++ b/lib/fpa/Base/MinimumSpanningTree.hxx @@ -90,6 +90,7 @@ GetPath( std::vector< V >& path, const V& a, const V& b ) const } // elihw if( raIt != ap.rbegin( ) ) --raIt; + if( rbIt != bp.rbegin( ) ) --rbIt; // Add part from a typename std::vector< V >::const_iterator iaIt = ap.begin( ); diff --git a/lib/fpa/Base/UniqueValuesContainer.h b/lib/fpa/Base/UniqueValuesContainer.h new file mode 100644 index 0000000..bc2b17c --- /dev/null +++ b/lib/fpa/Base/UniqueValuesContainer.h @@ -0,0 +1,72 @@ +#ifndef __FPA__BASE__UNIQUEVALUESCONTAINER__H__ +#define __FPA__BASE__UNIQUEVALUESCONTAINER__H__ + +#include +#include +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< class V, class VC > + class UniqueValuesContainer + : public itk::SimpleDataObjectDecorator< std::set< V, VC > > + { + public: + typedef std::set< V, VC > TDecorated; + typedef UniqueValuesContainer Self; + typedef itk::SimpleDataObjectDecorator< TDecorated > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef V TValue; + typedef VC TValueCompare; + + typedef typename TDecorated::iterator Iterator; + typedef typename TDecorated::const_iterator ConstIterator; + + public: + itkNewMacro( Self ); + itkTypeMacro( UniqueValuesContainer, itkSimpleDataObjectDecorator ); + + public: + void Insert( const V& v ) + { this->Get( ).insert( v ); this->Modified( ); } + void Erase( const V& v ) + { this->Get( ).erase( v ); this->Modified( ); } + void Clear( ) + { this->Get( ).clear( ); this->Modified( ); } + bool Find( const V& v ) const + { return( this->Get( ).find( v ) != this->Get( ).end( ) ); } + Iterator Begin( ) + { return( this->Get( ).begin( ) ); } + Iterator End( ) + { return( this->Get( ).end( ) ); } + ConstIterator Begin( ) const + { return( this->Get( ).begin( ) ); } + ConstIterator End( ) const + { return( this->Get( ).end( ) ); } + + protected: + UniqueValuesContainer( ) + : Superclass( ) + { } + virtual ~UniqueValuesContainer( ) + { } + + private: + // Purposely not implemented + UniqueValuesContainer( const Self& other ); + Self& operator=( const Self& other ); + }; + + } // ecapseman + +} // ecapseman + +#endif // __FPA__BASE__UNIQUEVALUESCONTAINER__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/IO/MatrixValuesContainerWriter.h b/lib/fpa/IO/MatrixValuesContainerWriter.h new file mode 100644 index 0000000..456716e --- /dev/null +++ b/lib/fpa/IO/MatrixValuesContainerWriter.h @@ -0,0 +1,61 @@ +#ifndef __FPA__IO__MATRIXVALUESCONTAINERWRITER__H__ +#define __FPA__IO__MATRIXVALUESCONTAINERWRITER__H__ + +#include +#include + +namespace fpa +{ + namespace IO + { + /** + */ + template< class T > + class MatrixValuesContainerWriter + : public itk::ProcessObject + { + public: + typedef MatrixValuesContainerWriter Self; + typedef itk::ProcessObject Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef T TTree; + + public: + itkNewMacro( Self ); + itkTypeMacro( MatrixValuesContainerWriter, itkProcessObject ); + + itkGetConstMacro( FileName, std::string ); + itkSetMacro( FileName, std::string ); + + public: + void SetInput( const T* input ); + T* GetInput( ); + const T* GetInput( ) const; + virtual void Update( ); + + protected: + MatrixValuesContainerWriter( ); + virtual ~MatrixValuesContainerWriter( ); + + virtual void GenerateData( ); + + private: + // Purposely not implemented + MatrixValuesContainerWriter( const Self& other ); + Self& operator=( const Self& other ); + + protected: + std::string m_FileName; + }; + + } // ecapseman + +} // ecapseman + +#include + +#endif // __FPA__IO__MATRIXVALUESCONTAINERWRITER__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/IO/MatrixValuesContainerWriter.hxx b/lib/fpa/IO/MatrixValuesContainerWriter.hxx new file mode 100644 index 0000000..dedc587 --- /dev/null +++ b/lib/fpa/IO/MatrixValuesContainerWriter.hxx @@ -0,0 +1,112 @@ +#ifndef __FPA__IO__MATRIXVALUESCONTAINERWRITER__HXX__ +#define __FPA__IO__MATRIXVALUESCONTAINERWRITER__HXX__ + +#include +#include + +// ------------------------------------------------------------------------- +template< class T > +void fpa::IO::MatrixValuesContainerWriter< T >:: +SetInput( const T* input ) +{ + this->itk::ProcessObject::SetNthInput( 0, const_cast< T* >( input ) ); +} + +// ------------------------------------------------------------------------- +template< class T > +T* fpa::IO::MatrixValuesContainerWriter< T >:: +GetInput( ) +{ + return( dynamic_cast< T* >( this->itk::ProcessObject::GetInput( 0 ) ) ); +} + +// ------------------------------------------------------------------------- +template< class T > +const T* fpa::IO::MatrixValuesContainerWriter< T >:: +GetInput( ) const +{ + return( + dynamic_cast< const T* >( this->itk::ProcessObject::GetInput( 0 ) ) + ); +} + +// ------------------------------------------------------------------------- +template< class T > +void fpa::IO::MatrixValuesContainerWriter< T >:: +Update( ) +{ + T* input = this->GetInput( ); + if( input != NULL ) + { + input->UpdateOutputInformation( ); + input->UpdateOutputData( ); + this->GenerateData( ); + this->ReleaseInputs( ); + + } // fi +} + +// ------------------------------------------------------------------------- +template< class T > +fpa::IO::MatrixValuesContainerWriter< T >:: +MatrixValuesContainerWriter( ) + : Superclass( ), + m_FileName( "" ) +{ + this->SetNumberOfRequiredInputs( 1 ); +} + +// ------------------------------------------------------------------------- +template< class T > +fpa::IO::MatrixValuesContainerWriter< T >:: +~MatrixValuesContainerWriter( ) +{ +} + +// ------------------------------------------------------------------------- +template< class T > +void fpa::IO::MatrixValuesContainerWriter< T >:: +GenerateData( ) +{ + const T* input = this->GetInput( ); + + std::ofstream out( this->m_FileName.c_str( ) ); + if( !out ) + { + itkExceptionMacro( + << "Error opening file to write a minimum spanning tree: \"" + << this->m_FileName + << "\"" + ); + return; + + } // fi + + out << T::TIndex::Dimension << std::endl; + + const typename T::TDecorated& real_input = input->Get( ); + out << real_input.size( ) << std::endl; + typename T::TDecorated::const_iterator cIt = real_input.begin( ); + for( ; cIt != real_input.end( ); ++cIt ) + { + out << cIt->second.size( ) << " "; + for( unsigned int d = 0; d < T::TIndex::Dimension; ++d ) + out << ( cIt->first )[ d ] << " "; + + typename T::TDecorated::value_type::second_type::const_iterator rIt = + cIt->second.begin( ); + for( ; rIt != cIt->second.end( ); ++rIt ) + { + for( unsigned int d = 0; d < T::TIndex::Dimension; ++d ) + out << ( rIt->first )[ d ] << " "; + out << rIt->second << std::endl; + + } // rof + + } // rof + out.close( ); +} + +#endif // __FPA__IO__MATRIXVALUESCONTAINERWRITER__HXX__ + +// eof - $RCSfile$ diff --git a/lib/fpa/IO/MinimumSpanningTreeWriter.h b/lib/fpa/IO/MinimumSpanningTreeWriter.h new file mode 100644 index 0000000..3c3e9d1 --- /dev/null +++ b/lib/fpa/IO/MinimumSpanningTreeWriter.h @@ -0,0 +1,61 @@ +#ifndef __FPA__IO__MINIMUMSPANNINGTREEWRITER__H__ +#define __FPA__IO__MINIMUMSPANNINGTREEWRITER__H__ + +#include +#include + +namespace fpa +{ + namespace IO + { + /** + */ + template< class T > + class MinimumSpanningTreeWriter + : public itk::ProcessObject + { + public: + typedef MinimumSpanningTreeWriter Self; + typedef itk::ProcessObject Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef T TTree; + + public: + itkNewMacro( Self ); + itkTypeMacro( MinimumSpanningTreeWriter, itkProcessObject ); + + itkGetConstMacro( FileName, std::string ); + itkSetMacro( FileName, std::string ); + + public: + void SetInput( const T* input ); + T* GetInput( ); + const T* GetInput( ) const; + virtual void Update( ); + + protected: + MinimumSpanningTreeWriter( ); + virtual ~MinimumSpanningTreeWriter( ); + + virtual void GenerateData( ); + + private: + // Purposely not implemented + MinimumSpanningTreeWriter( const Self& other ); + Self& operator=( const Self& other ); + + protected: + std::string m_FileName; + }; + + } // ecapseman + +} // ecapseman + +#include + +#endif // __FPA__IO__MINIMUMSPANNINGTREEWRITER__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/IO/MinimumSpanningTreeWriter.hxx b/lib/fpa/IO/MinimumSpanningTreeWriter.hxx new file mode 100644 index 0000000..f9b8533 --- /dev/null +++ b/lib/fpa/IO/MinimumSpanningTreeWriter.hxx @@ -0,0 +1,118 @@ +#ifndef __FPA__IO__MINIMUMSPANNINGTREEWRITER__HXX__ +#define __FPA__IO__MINIMUMSPANNINGTREEWRITER__HXX__ + +#include +#include + +// ------------------------------------------------------------------------- +template< class T > +void fpa::IO::MinimumSpanningTreeWriter< T >:: +SetInput( const T* input ) +{ + this->itk::ProcessObject::SetNthInput( 0, const_cast< T* >( input ) ); +} + +// ------------------------------------------------------------------------- +template< class T > +T* fpa::IO::MinimumSpanningTreeWriter< T >:: +GetInput( ) +{ + return( dynamic_cast< T* >( this->itk::ProcessObject::GetInput( 0 ) ) ); +} + +// ------------------------------------------------------------------------- +template< class T > +const T* fpa::IO::MinimumSpanningTreeWriter< T >:: +GetInput( ) const +{ + return( + dynamic_cast< const T* >( this->itk::ProcessObject::GetInput( 0 ) ) + ); +} + +// ------------------------------------------------------------------------- +template< class T > +void fpa::IO::MinimumSpanningTreeWriter< T >:: +Update( ) +{ + T* input = this->GetInput( ); + if( input != NULL ) + { + input->UpdateOutputInformation( ); + input->UpdateOutputData( ); + this->GenerateData( ); + this->ReleaseInputs( ); + + } // fi +} + +// ------------------------------------------------------------------------- +template< class T > +fpa::IO::MinimumSpanningTreeWriter< T >:: +MinimumSpanningTreeWriter( ) + : Superclass( ), + m_FileName( "" ) +{ + this->SetNumberOfRequiredInputs( 1 ); +} + +// ------------------------------------------------------------------------- +template< class T > +fpa::IO::MinimumSpanningTreeWriter< T >:: +~MinimumSpanningTreeWriter( ) +{ +} + +// ------------------------------------------------------------------------- +template< class T > +void fpa::IO::MinimumSpanningTreeWriter< T >:: +GenerateData( ) +{ + const T* input = this->GetInput( ); + + std::ofstream out( this->m_FileName.c_str( ) ); + if( !out ) + { + itkExceptionMacro( + << "Error opening file to write a minimum spanning tree: \"" + << this->m_FileName + << "\"" + ); + return; + + } // fi + + out << T::TVertex::Dimension << std::endl; + + const typename T::TCollisions& collisions = input->GetCollisions( ); + unsigned long nSeeds = collisions.size( ); + out << nSeeds << std::endl; + for( unsigned long i = 0; i < nSeeds; ++i ) + { + for( unsigned long j = 0; j < nSeeds; ++j ) + out << collisions[ i ][ j ].second << " "; + out << std::endl; + + } // rof + + const typename T::TDecorated& real_input = input->Get( ); + out << real_input.size( ) << std::endl; + for( + typename T::TDecorated::const_iterator iIt = real_input.begin( ); + iIt != real_input.end( ); + ++iIt + ) + { + for( unsigned int d = 0; d < T::TVertex::Dimension; ++d ) + out << iIt->first[ d ] << " "; + for( unsigned int d = 0; d < T::TVertex::Dimension; ++d ) + out << iIt->second.first[ d ] << " "; + out << iIt->second.second << std::endl; + + } // rof + out.close( ); +} + +#endif // __FPA__IO__MINIMUMSPANNINGTREEWRITER__HXX__ + +// eof - $RCSfile$ diff --git a/lib/fpa/IO/UniqueValuesContainerWriter.h b/lib/fpa/IO/UniqueValuesContainerWriter.h new file mode 100644 index 0000000..fc27012 --- /dev/null +++ b/lib/fpa/IO/UniqueValuesContainerWriter.h @@ -0,0 +1,61 @@ +#ifndef __FPA__IO__UNIQUEVALUESCONTAINERWRITER__H__ +#define __FPA__IO__UNIQUEVALUESCONTAINERWRITER__H__ + +#include +#include + +namespace fpa +{ + namespace IO + { + /** + */ + template< class T > + class UniqueValuesContainerWriter + : public itk::ProcessObject + { + public: + typedef UniqueValuesContainerWriter Self; + typedef itk::ProcessObject Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef T TTree; + + public: + itkNewMacro( Self ); + itkTypeMacro( UniqueValuesContainerWriter, itkProcessObject ); + + itkGetConstMacro( FileName, std::string ); + itkSetMacro( FileName, std::string ); + + public: + void SetInput( const T* input ); + T* GetInput( ); + const T* GetInput( ) const; + virtual void Update( ); + + protected: + UniqueValuesContainerWriter( ); + virtual ~UniqueValuesContainerWriter( ); + + virtual void GenerateData( ); + + private: + // Purposely not implemented + UniqueValuesContainerWriter( const Self& other ); + Self& operator=( const Self& other ); + + protected: + std::string m_FileName; + }; + + } // ecapseman + +} // ecapseman + +#include + +#endif // __FPA__IO__UNIQUEVALUESCONTAINERWRITER__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/IO/UniqueValuesContainerWriter.hxx b/lib/fpa/IO/UniqueValuesContainerWriter.hxx new file mode 100644 index 0000000..17eb1d1 --- /dev/null +++ b/lib/fpa/IO/UniqueValuesContainerWriter.hxx @@ -0,0 +1,105 @@ +#ifndef __FPA__IO__UNIQUEVALUESCONTAINERWRITER__HXX__ +#define __FPA__IO__UNIQUEVALUESCONTAINERWRITER__HXX__ + +#include +#include + +// ------------------------------------------------------------------------- +template< class T > +void fpa::IO::UniqueValuesContainerWriter< T >:: +SetInput( const T* input ) +{ + this->itk::ProcessObject::SetNthInput( 0, const_cast< T* >( input ) ); +} + +// ------------------------------------------------------------------------- +template< class T > +T* fpa::IO::UniqueValuesContainerWriter< T >:: +GetInput( ) +{ + return( dynamic_cast< T* >( this->itk::ProcessObject::GetInput( 0 ) ) ); +} + +// ------------------------------------------------------------------------- +template< class T > +const T* fpa::IO::UniqueValuesContainerWriter< T >:: +GetInput( ) const +{ + return( + dynamic_cast< const T* >( this->itk::ProcessObject::GetInput( 0 ) ) + ); +} + +// ------------------------------------------------------------------------- +template< class T > +void fpa::IO::UniqueValuesContainerWriter< T >:: +Update( ) +{ + T* input = this->GetInput( ); + if( input != NULL ) + { + input->UpdateOutputInformation( ); + input->UpdateOutputData( ); + this->GenerateData( ); + this->ReleaseInputs( ); + + } // fi +} + +// ------------------------------------------------------------------------- +template< class T > +fpa::IO::UniqueValuesContainerWriter< T >:: +UniqueValuesContainerWriter( ) + : Superclass( ), + m_FileName( "" ) +{ + this->SetNumberOfRequiredInputs( 1 ); +} + +// ------------------------------------------------------------------------- +template< class T > +fpa::IO::UniqueValuesContainerWriter< T >:: +~UniqueValuesContainerWriter( ) +{ +} + +// ------------------------------------------------------------------------- +template< class T > +void fpa::IO::UniqueValuesContainerWriter< T >:: +GenerateData( ) +{ + const T* input = this->GetInput( ); + + std::ofstream out( this->m_FileName.c_str( ) ); + if( !out ) + { + itkExceptionMacro( + << "Error opening file to write a minimum spanning tree: \"" + << this->m_FileName + << "\"" + ); + return; + + } // fi + + out << T::TValue::Dimension << std::endl; + + const typename T::TDecorated& real_input = input->Get( ); + out << real_input.size( ) << std::endl; + for( + typename T::TDecorated::const_iterator iIt = real_input.begin( ); + iIt != real_input.end( ); + ++iIt + ) + { + for( unsigned int d = 0; d < T::TValue::Dimension; ++d ) + out << ( *iIt )[ d ] << " "; + out << std::endl; + + } // rof + out.close( ); +} + +#endif // __FPA__IO__UNIQUEVALUESCONTAINERWRITER__HXX__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Image/Algorithm.h b/lib/fpa/Image/Algorithm.h index f9ec7c4..efb376b 100644 --- a/lib/fpa/Image/Algorithm.h +++ b/lib/fpa/Image/Algorithm.h @@ -29,10 +29,11 @@ namespace fpa typedef I TInputImage; typedef O TOutputImage; - typedef typename Superclass::TVertex TVertex; - typedef typename Superclass::TValue TValue; - typedef typename Superclass::TResult TResult; - typedef typename Superclass::TVertexCompare TVertexCompare; + typedef typename Superclass::TVertex TVertex; + typedef typename Superclass::TValue TValue; + typedef typename Superclass::TResult TResult; + typedef typename Superclass::TVertexCompare TVertexCompare; + typedef typename Superclass::TMinimumSpanningTree TMinimumSpanningTree; protected: typedef typename Superclass::_TVertices _TVertices; diff --git a/lib/fpa/Image/Dijkstra.h b/lib/fpa/Image/Dijkstra.h index 48d7a79..e38fd0d 100644 --- a/lib/fpa/Image/Dijkstra.h +++ b/lib/fpa/Image/Dijkstra.h @@ -27,11 +27,12 @@ namespace fpa typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; - typedef typename Superclass::TInputImage TInputImage; - typedef typename Superclass::TOutputImage TOutputImage; - typedef typename Superclass::TVertex TVertex; - typedef typename Superclass::TValue TValue; - typedef typename Superclass::TResult TResult; + typedef typename Superclass::TInputImage TInputImage; + typedef typename Superclass::TOutputImage TOutputImage; + typedef typename Superclass::TVertex TVertex; + typedef typename Superclass::TValue TValue; + typedef typename Superclass::TResult TResult; + typedef typename Superclass::TMinimumSpanningTree TMinimumSpanningTree; typedef typename Superclass::TStartEvent TStartEvent; typedef typename Superclass::TStartLoopEvent TStartLoopEvent; diff --git a/lib/fpa/Image/DijkstraWithEndPointDetection.h b/lib/fpa/Image/DijkstraWithEndPointDetection.h index c63f4f1..4288da3 100644 --- a/lib/fpa/Image/DijkstraWithEndPointDetection.h +++ b/lib/fpa/Image/DijkstraWithEndPointDetection.h @@ -4,6 +4,8 @@ #include #include #include +#include +#include namespace fpa { @@ -23,15 +25,16 @@ namespace fpa typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; - typedef typename Superclass::TInputImage TInputImage; - typedef typename Superclass::TOutputImage TOutputImage; - typedef typename Superclass::TVertex TVertex; - typedef typename Superclass::TVertexCompare TVertexCompare; - typedef typename Superclass::TValue TValue; - typedef typename Superclass::TResult TResult; - typedef typename Superclass::TCostFunction TCostFunction; - typedef typename Superclass::TConversionFunction TConversionFunction; - typedef typename Superclass::_TVertices TVertices; + typedef typename Superclass::TInputImage TInputImage; + typedef typename Superclass::TOutputImage TOutputImage; + typedef typename Superclass::TVertex TVertex; + typedef typename Superclass::TVertexCompare TVertexCompare; + typedef typename Superclass::TValue TValue; + typedef typename Superclass::TResult TResult; + typedef typename Superclass::TMinimumSpanningTree TMinimumSpanningTree; + typedef typename Superclass::TCostFunction TCostFunction; + typedef typename Superclass::TConversionFunction TConversionFunction; + typedef typename Superclass::_TVertices TVertices; typedef typename Superclass::TStartEvent TStartEvent; typedef typename Superclass::TStartLoopEvent TStartLoopEvent; @@ -51,9 +54,8 @@ namespace fpa typedef unsigned short TLabel; typedef itk::Image< TLabel, I::ImageDimension > TLabelImage; - typedef std::set< TVertex, TVertexCompare > TUniqueVertices; - typedef std::map< TVertex, TLabel, TVertexCompare > TBranch; - typedef std::map< TVertex, TBranch, TVertexCompare > TBranches; + typedef fpa::Base::UniqueValuesContainer< TVertex, TVertexCompare > TUniqueVertices; + typedef fpa::Base::MatrixValuesContainer< TVertex, TLabel, TVertexCompare > TBranches; protected: typedef typename Superclass::_TVertices _TVertices; @@ -74,16 +76,25 @@ namespace fpa itkNewMacro( Self ); itkTypeMacro( DijkstraWithEndPointDetection, Dijkstra ); - itkGetConstMacro( EndPoints, TUniqueVertices ); - itkGetConstMacro( BifurcationPoints, TUniqueVertices ); itkGetConstMacro( NumberOfBranches, unsigned long ); - itkGetConstMacro( Branches, TBranches ); public: TLabelImage* GetLabelImage( ); const TLabelImage* GetLabelImage( ) const; void GraftLabelImage( itk::DataObject* obj ); + TUniqueVertices* GetEndPoints( ); + const TUniqueVertices* GetEndPoints( ) const; + void GraftEndPoints( itk::DataObject* obj ); + + TUniqueVertices* GetBifurcations( ); + const TUniqueVertices* GetBifurcations( ) const; + void GraftBifurcations( itk::DataObject* obj ); + + TBranches* GetBranches( ); + const TBranches* GetBranches( ) const; + void GraftBranches( itk::DataObject* obj ); + protected: DijkstraWithEndPointDetection( ); virtual ~DijkstraWithEndPointDetection( ); @@ -92,6 +103,10 @@ namespace fpa virtual void _AfterGenerateData( ); virtual void _SetResult( const TVertex& v, const _TNode& n ); + void _EndPointsAndBifurcations( ); + void _FindBranches( ); + void _LabelAll( ); + _TRegion _Region( const TVertex& c, const double& r ); template< class _T > @@ -99,8 +114,6 @@ namespace fpa const _T* image, const TVertex& v, const double& r ); - void _Label( const TVertex& v, const TLabel& l ); - private: // Purposely not implemented DijkstraWithEndPointDetection( const Self& other ); @@ -108,12 +121,12 @@ namespace fpa protected: unsigned int m_LabelImageIndex; + unsigned int m_BifurcationsIndex; + unsigned int m_EndPointsIndex; + unsigned int m_BranchesIndex; _TCandidates m_Candidates; - TUniqueVertices m_BifurcationPoints; - TUniqueVertices m_EndPoints; unsigned long m_NumberOfBranches; - TBranches m_Branches; }; } // ecapseman diff --git a/lib/fpa/Image/DijkstraWithEndPointDetection.hxx b/lib/fpa/Image/DijkstraWithEndPointDetection.hxx index 8ddd9ba..76695c2 100644 --- a/lib/fpa/Image/DijkstraWithEndPointDetection.hxx +++ b/lib/fpa/Image/DijkstraWithEndPointDetection.hxx @@ -4,6 +4,7 @@ #include #include #include +#include // ------------------------------------------------------------------------- template< class I, class O > @@ -44,6 +45,127 @@ GraftLabelImage( itk::DataObject* obj ) this->GraftNthOutput( this->m_LabelImageIndex, lbl ); } + + + +// ------------------------------------------------------------------------- +template< class I, class O > +typename fpa::Image::DijkstraWithEndPointDetection< I, O >:: +TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GetEndPoints( ) +{ + return( + dynamic_cast< TUniqueVertices* >( + this->itk::ProcessObject::GetOutput( this->m_EndPointsIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +const typename fpa::Image::DijkstraWithEndPointDetection< I, O >:: +TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GetEndPoints( ) const +{ + return( + dynamic_cast< const TUniqueVertices* >( + this->itk::ProcessObject::GetOutput( this->m_EndPointsIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GraftEndPoints( itk::DataObject* obj ) +{ + TUniqueVertices* lbl = + dynamic_cast< TUniqueVertices* >( + this->itk::ProcessObject::GetOutput( this->m_EndPointsIndex ) + ); + if( lbl != NULL ) + this->GraftNthOutput( this->m_EndPointsIndex, lbl ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +typename fpa::Image::DijkstraWithEndPointDetection< I, O >:: +TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GetBifurcations( ) +{ + return( + dynamic_cast< TUniqueVertices* >( + this->itk::ProcessObject::GetOutput( this->m_BifurcationsIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +const typename fpa::Image::DijkstraWithEndPointDetection< I, O >:: +TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GetBifurcations( ) const +{ + return( + dynamic_cast< const TUniqueVertices* >( + this->itk::ProcessObject::GetOutput( this->m_BifurcationsIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GraftBifurcations( itk::DataObject* obj ) +{ + TUniqueVertices* lbl = + dynamic_cast< TUniqueVertices* >( + this->itk::ProcessObject::GetOutput( this->m_BifurcationsIndex ) + ); + if( lbl != NULL ) + this->GraftNthOutput( this->m_BifurcationsIndex, lbl ); +} + + +// ------------------------------------------------------------------------- +template< class I, class O > +typename fpa::Image::DijkstraWithEndPointDetection< I, O >:: +TBranches* fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GetBranches( ) +{ + return( + dynamic_cast< TBranches* >( + this->itk::ProcessObject::GetOutput( this->m_BranchesIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +const typename fpa::Image::DijkstraWithEndPointDetection< I, O >:: +TBranches* fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GetBranches( ) const +{ + return( + dynamic_cast< const TBranches* >( + this->itk::ProcessObject::GetOutput( this->m_BranchesIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GraftBranches( itk::DataObject* obj ) +{ + TBranches* lbl = + dynamic_cast< TBranches* >( + this->itk::ProcessObject::GetOutput( this->m_BranchesIndex ) + ); + if( lbl != NULL ) + this->GraftNthOutput( this->m_BranchesIndex, lbl ); +} + // ------------------------------------------------------------------------- template< class I, class O > fpa::Image::DijkstraWithEndPointDetection< I, O >:: @@ -51,10 +173,22 @@ DijkstraWithEndPointDetection( ) : Superclass( ) { this->m_LabelImageIndex = this->GetNumberOfRequiredOutputs( ); - this->SetNumberOfRequiredOutputs( this->m_LabelImageIndex + 1 ); + this->m_EndPointsIndex = this->m_LabelImageIndex + 1; + this->m_BifurcationsIndex = this->m_LabelImageIndex + 2; + this->m_BranchesIndex = this->m_LabelImageIndex + 3; + this->SetNumberOfRequiredOutputs( this->m_LabelImageIndex + 4 ); this->itk::ProcessObject::SetNthOutput( this->m_LabelImageIndex, TLabelImage::New( ) ); + this->itk::ProcessObject::SetNthOutput( + this->m_EndPointsIndex, TUniqueVertices::New( ) + ); + this->itk::ProcessObject::SetNthOutput( + this->m_BifurcationsIndex, TUniqueVertices::New( ) + ); + this->itk::ProcessObject::SetNthOutput( + this->m_BranchesIndex, TBranches::New( ) + ); } // ------------------------------------------------------------------------- @@ -69,31 +203,39 @@ template< class I, class O > void fpa::Image::DijkstraWithEndPointDetection< I, O >:: _BeforeGenerateData( ) { + const I* input = this->GetInput( ); + typename I::SpacingType spac = input->GetSpacing( ); + double max_spac = double( spac[ 0 ] ); + for( unsigned int d = 1; d < I::ImageDimension; ++d ) + max_spac = + ( double( spac[ d ] ) > max_spac )? double( spac[ d ] ): max_spac; + max_spac *= double( 3 ); + // Correct seeds - /* TODO - for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s ) - { - _TNode seed = this->m_Seeds[ s ]; - _TRegion region = this->_Region( seed.Vertex, max_spac ); - itk::ImageRegionConstIteratorWithIndex< I > iIt( img, region ); - - iIt.GoToBegin( ); - _TPixel max_value = iIt.Get( ); - for( ++iIt; !iIt.IsAtEnd( ); ++iIt ) - { - if( iIt.Get( ) > max_value ) - { - seed.Vertex = iIt.GetIndex( ); - seed.Parent = seed.Vertex; - max_value = iIt.Get( ); + for( unsigned int s = 0; s < this->m_SeedVertices.size( ); ++s ) + { + TVertex seed = this->m_SeedVertices[ s ]; + _TNode n = this->m_Seeds[ seed ]; + _TRegion region = this->_Region( seed, max_spac ); + itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region ); + + iIt.GoToBegin( ); + _TPixel max_value = iIt.Get( ); + for( ++iIt; !iIt.IsAtEnd( ); ++iIt ) + { + if( iIt.Get( ) > max_value ) + { + this->m_SeedVertices[ s ] = iIt.GetIndex( ); + max_value = iIt.Get( ); - } // fi + } // fi - } // rof - this->m_Seeds[ s ] = seed; + } // rof + this->m_Seeds.erase( seed ); + n.Parent = this->m_SeedVertices[ s ]; + this->m_Seeds[ this->m_SeedVertices[ s ] ] = n; - } // rof - */ + } // rof // End initialization this->Superclass::_BeforeGenerateData( ); @@ -108,240 +250,245 @@ _AfterGenerateData( ) // Finish base algorithm this->Superclass::_AfterGenerateData( ); - // Prepare backtracking objects - this->m_EndPoints.clear( ); - this->m_BifurcationPoints.clear( ); + // Check if backtracking has any sense if( this->m_Candidates.size( ) == 0 ) return; this->InvokeEvent( TEndEvent( ) ); this->InvokeEvent( TStartBacktrackingEvent( ) ); - // Get some input values - const I* input = this->GetInput( ); - typename I::SpacingType spac = input->GetSpacing( ); - double ms = double( spac[ 0 ] ); - for( unsigned int d = 1; d < I::ImageDimension; ++d ) - ms =( ms < double( spac[ d ] ) )? double( spac[ d ] ): ms; + // Detect endpoints and bifurcations + this->_EndPointsAndBifurcations( ); - // Prepare labels - TLabelImage* label = this->GetLabelImage( ); - label->FillBuffer( 0 ); + // Find branches + this->_FindBranches( ); + + // Label pixels and branches + this->_LabelAll( ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +_SetResult( const TVertex& v, const _TNode& n ) +{ + this->Superclass::_SetResult( v, n ); + this->m_Candidates.insert( _TCandidate( n.Result, v ) ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +_EndPointsAndBifurcations( ) +{ + // Prepare backtracking objects + TUniqueVertices* endpoints = this->GetEndPoints( ); + TUniqueVertices* bifurcations = this->GetBifurcations( ); + endpoints->Clear( ); + bifurcations->Clear( ); - // Object to hold marks - std::set< TVertex, TVertexCompare > tree_marks; + // Get some input values + TVertex seed = this->GetSeed( 0 ); + const I* input = this->GetInput( ); + const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( ); - // Object to hold all branches - this->m_Branches.clear( ); + // Prepare a pixel marks structure + typedef std::map< TVertex, bool, TVertexCompare > _TPixelMarks; + _TPixelMarks pixel_marks; - // First label - this->m_NumberOfBranches = 1; + // Object to hold tree-nodes marks + typedef std::set< TVertex, TVertexCompare > _TTreeMarks; + _TTreeMarks tree_marks; // Iterate over the candidates, starting from the candidates that - // are near thin branches + // are near thin branches (high accumulated cost) typename _TCandidates::const_reverse_iterator cIt = this->m_Candidates.rbegin( ); for( ; cIt != this->m_Candidates.rend( ); ++cIt ) { - // If pixel has been already labelled, pass + // If pixel has already been visited, pass TVertex v = cIt->second; - if( label->GetPixel( v ) != 0 ) + if( pixel_marks[ v ] ) continue; - // Compute nearest start candidate - TVertex endpoint = - this->_MaxInRegion( - input, v, - double( std::sqrt( input->GetPixel( v ) ) ) * double( 1.5 ) - ); + // Correct it to nearest start candidate (high distance value) + double vr = std::sqrt( double( input->GetPixel( v ) ) ); + v = this->_MaxInRegion( input, v, vr * double( 1.5 ) ); - // Re-check labelling - if( this->_Node( endpoint ).Label == Self::FarLabel ) + // Now, check for real marking conditions + // 1. Has it been visited by dijkstra? + if( this->_Node( v ).Label == Self::FarLabel ) continue; - if( label->GetPixel( endpoint ) != 0 ) + // 2. Is it already marked? + if( pixel_marks[ v ] ) continue; - if( this->m_EndPoints.find( endpoint ) != this->m_EndPoints.end( ) ) + // 3. Is it already an endpoint? + if( endpoints->Find( v ) ) continue; - this->m_EndPoints.insert( endpoint ); - std::cout << "endpoint " << endpoint << " inserted " << this->m_EndPoints.size( ) << std::endl; - - // Get the path all the way to seed - std::vector< TVertex > path; - this->GetMinimumSpanningTree( )-> - GetPath( path, endpoint, this->GetSeed( 0 ) ); - - // Mark branches - bool start = true; - bool change = false; - TVertex last_start = endpoint; - typename std::vector< TVertex >::const_iterator pIt = path.begin( ); - for( ; pIt != path.end( ); ++pIt ) + // 4. Ok, it is completely new! + endpoints->Insert( v ); + + // Get the path all the way to global seed + TVertices path; + mst->GetPath( path, v, seed ); + + // Backtracking to find endpoints and bifurcations + bool adding_new_points = true; + typename TVertices::const_iterator pIt = path.begin( ); + for( ; pIt != path.end( ) && adding_new_points; ++pIt ) { - this->InvokeEvent( - TBacktrackingEvent( *pIt, this->m_NumberOfBranches ) - ); - if( start ) + this->InvokeEvent( TBacktrackingEvent( *pIt, 1 ) ); + if( tree_marks.find( *pIt ) == tree_marks.end( ) ) { - if( tree_marks.find( *pIt ) == tree_marks.end( ) ) + // Mark current point as a tree point + tree_marks.insert( *pIt ); + + // Mark a region around current point as visited + vr = std::sqrt( double( input->GetPixel( *pIt ) ) ); + _TRegion region = this->_Region( *pIt, vr ); + if( region.GetNumberOfPixels( ) > 0 ) { - // Mark a region around current point as visited - tree_marks.insert( *pIt ); - this->_Label( *pIt, TLabel( this->m_NumberOfBranches ) ); + itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region ); + for( iIt.GoToBegin( ); !iIt.IsAtEnd( ); ++iIt ) + pixel_marks[ iIt.GetIndex( ) ] = true; } else - { - // A bifurcation point has been reached! - if( *pIt != this->GetSeed( 0 ) ) - { - this->m_Branches[ *pIt ][ last_start ] = this->m_NumberOfBranches; - last_start = *pIt; - this->m_BifurcationPoints.insert( *pIt ); - std::cout << "bifurcation = " << *pIt << " " << this->GetSeed( 0 ) << std::endl; - this->m_NumberOfBranches++; - start = false; - - } // fi - - } // fi + pixel_marks[ *pIt ] = true; } else { - if( - std::find( - this->m_BifurcationPoints.begin( ), - this->m_BifurcationPoints.end( ), - *pIt - ) == this->m_BifurcationPoints.end( ) - ) + // A bifurcation point has been reached! + if( *pIt != seed ) { - // Mark a sphere around current point as visited - this->_Label( *pIt, TLabel( this->m_NumberOfBranches ) ); - change = true; - } - else - { - // A bifurcation point has been reached! - if( *pIt != this->GetSeed( 0 ) ) - { - this->m_Branches[ *pIt ][ last_start ] = this->m_NumberOfBranches; - last_start = *pIt; - this->m_NumberOfBranches++; - - } // fi + bifurcations->Insert( *pIt ); + adding_new_points = false; } // fi } // fi } // rof - if( start || change ) - this->m_NumberOfBranches++; this->InvokeEvent( TEndBacktrackingEvent( ) ); } // rof +} - typename TBranches::const_iterator bIt = this->m_Branches.begin( ); - unsigned int leo = 0; - for( ; bIt != this->m_Branches.end( ); ++bIt ) +// ------------------------------------------------------------------------- +template< class I, class O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +_FindBranches( ) +{ + // Configure and obtain inputs + TVertex seed = this->GetSeed( 0 ); + const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( ); + const TUniqueVertices* endpoints = this->GetEndPoints( ); + const TUniqueVertices* bifurcations = this->GetBifurcations( ); + TBranches* branches = this->GetBranches( ); + branches->Clear( ); + + // Reconstruct pixels + typename TUniqueVertices::ConstIterator eIt = endpoints->Begin( ); + for( ; eIt != endpoints->End( ); ++eIt ) { - typename TBranch::const_iterator brIt = bIt->second.begin( ); - for( ; brIt != bIt->second.end( ); ++brIt ) + // Get the path all the way to global seed + TVertices path; + mst->GetPath( path, *eIt, seed ); + + TVertex start_vertex = *eIt; + typename TVertices::const_iterator pIt = path.begin( ); + for( ; pIt != path.end( ); ++pIt ) { - std::cout << bIt->first << " " << brIt->first << std::endl; - leo++; - } - } - - // Re-enumerate labels - std::map< TLabel, unsigned long > histo; - for( - typename std::set< TVertex, TVertexCompare >::iterator treeIt = tree_marks.begin( ); - treeIt != tree_marks.end( ); - ++treeIt - ) - histo[ label->GetPixel( *treeIt ) ]++; - - std::map< TLabel, TLabel > changes; - TLabel last_change = 1; - for( TLabel i = 1; i <= this->m_NumberOfBranches; ++i ) - { - if( histo[ i ] != 0 ) - changes[ i ] = last_change++; + if( bifurcations->Find( *pIt ) ) + { + branches->SetValue( start_vertex, *pIt, 1 ); + start_vertex = *pIt; - } // rof - this->m_NumberOfBranches = changes.size( ); + } // fi - std::cout << leo << " - " << this->m_EndPoints.size( ) << " - " << this->m_BifurcationPoints.size( ) << " - " << this->m_NumberOfBranches << std::endl; + } // rof - /* - for( - typename TTree::iterator treeIt = this->m_FullTree.begin( ); - treeIt != this->m_FullTree.end( ); - ++treeIt - ) - { - TMark old = treeIt->second.second; - treeIt->second.second = changes[ old ]; + // Finish with branches to global seed + branches->SetValue( start_vertex, seed, 1 ); - } // fi + } // rof +} - // Construct reduced tree - for( - typename TVertices::const_iterator eIt = this->m_EndPoints.begin( ); - eIt != this->m_EndPoints.end( ); - ++eIt - ) - { - typename TTree::const_iterator tIt = this->m_FullTree.find( *eIt ); - if( tIt != this->m_FullTree.end( ) ) - { - TMark id = tIt->second.second; - do +// ------------------------------------------------------------------------- +template< class I, class O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +_LabelAll( ) +{ + // Configure and obtain inputs + const I* input = this->GetInput( ); + const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( ); + TBranches* branches = this->GetBranches( ); + TLabelImage* label = this->GetLabelImage( ); + label->FillBuffer( 0 ); + + // Label branches + typename TBranches::Iterator bIt = branches->Begin( ); + TLabel actual_label = 0; + for( ; bIt != branches->End( ); ++bIt ) + { + typename TBranches::RowIterator brIt = branches->Begin( bIt ); + for( ; brIt != branches->End( bIt ); ++brIt ) { - tIt = this->m_FullTree.find( tIt->second.first ); - if( tIt == this->m_FullTree.end( ) ) - break; + actual_label++; + brIt->second = actual_label; - } while( tIt->second.second == id ); - this->m_ReducedTree[ *eIt ] = TTreeNode( tIt->first, id ); + TVertices path; + mst->GetPath( path, bIt->first, brIt->first ); + typename TVertices::const_iterator pIt = path.begin( ); + for( ; pIt != path.end( ); ++pIt ) + { + double d = std::sqrt( double( input->GetPixel( *pIt ) ) ); + _TRegion region = this->_Region( *pIt, d ); + itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region ); + itk::ImageRegionIteratorWithIndex< TLabelImage > lIt( label, region ); + iIt.GoToBegin( ); + lIt.GoToBegin( ); + for( ; !iIt.IsAtEnd( ); ++iIt, ++lIt ) + { + // Mask real input image + if( + !( iIt.Get( ) < ( typename I::PixelType )( 0 ) ) && + lIt.Get( ) == TLabel( 0 ) + ) + lIt.Set( actual_label ); + + } // rof - } // fi + } // rof } // rof - for( - typename TVertices::const_iterator bIt = this->m_BifurcationPoints.begin( ); - bIt != this->m_BifurcationPoints.end( ); - ++bIt - ) - { - typename TTree::const_iterator tIt = this->m_FullTree.find( *bIt ); - if( tIt != this->m_FullTree.end( ) ) - { - TMark id = tIt->second.second; - do - { - tIt = this->m_FullTree.find( tIt->second.first ); - if( tIt == this->m_FullTree.end( ) ) - break; + } // rof + this->m_NumberOfBranches = actual_label; - } while( tIt->second.second == id ); - this->m_ReducedTree[ *bIt ] = TTreeNode( tIt->first, id ); + // Label bifurcations + /* TODO: Is this necessary? + typename TUniqueVertices::const_iterator bifIt = + this->m_Bifurcations.begin( ); + for( ; bifIt != this->m_Bifurcations.end( ); ++bifIt ) + { + actual_label++; + double d = std::sqrt( double( input->GetPixel( *bifIt ) ) ); + _TRegion region = this->_Region( *bifIt, d ); + itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region ); + itk::ImageRegionIteratorWithIndex< TLabelImage > lIt( label, region ); + iIt.GoToBegin( ); + lIt.GoToBegin( ); + for( ; !iIt.IsAtEnd( ); ++iIt, ++lIt ) + { + // Mask real input image + if( !( iIt.Get( ) < ( typename I::PixelType )( 0 ) ) ) + lIt.Set( actual_label ); - } // fi + } // rof - } // rof + } // rof */ } -// ------------------------------------------------------------------------- -template< class I, class O > -void fpa::Image::DijkstraWithEndPointDetection< I, O >:: -_SetResult( const TVertex& v, const _TNode& n ) -{ - this->Superclass::_SetResult( v, n ); - this->m_Candidates.insert( _TCandidate( n.Result, v ) ); -} - // ------------------------------------------------------------------------- template< class I, class O > typename fpa::Image::DijkstraWithEndPointDetection< I, O >:: @@ -359,6 +506,7 @@ _Region( const TVertex& c, const double& r ) _TSize size; for( unsigned int d = 0; d < I::ImageDimension; ++d ) { + // NOTE: 3 is a minimum neighborhood size long s = long( std::ceil( r / double( spac[ d ] ) ) ) + 3; i0[ d ] = c[ d ] - s; i1[ d ] = c[ d ] + s; @@ -405,25 +553,6 @@ _MaxInRegion( const _T* image, const TVertex& v, const double& r ) return( max_vertex ); } -// ------------------------------------------------------------------------- -template< class I, class O > -void fpa::Image::DijkstraWithEndPointDetection< I, O >:: -_Label( const TVertex& v, const TLabel& l ) -{ - typedef itk::ImageRegionIteratorWithIndex< TLabelImage > _TIt; - - double d = std::sqrt( double( this->GetInput( )->GetPixel( v ) ) ); - _TRegion region = this->_Region( v, d ); - if( region.GetNumberOfPixels( ) > 0 ) - { - _TIt lIt( this->GetLabelImage( ), region ); - for( lIt.GoToBegin( ); !lIt.IsAtEnd( ); ++lIt ) - lIt.Set( l ); - } - else - this->GetLabelImage( )->SetPixel( v, l ); -} - #endif // __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__ // eof - $RCSfile$ diff --git a/lib/fpa/VTK/ImageMPR.cxx b/lib/fpa/VTK/ImageMPR.cxx index 29f4dc4..1cc4b85 100644 --- a/lib/fpa/VTK/ImageMPR.cxx +++ b/lib/fpa/VTK/ImageMPR.cxx @@ -234,6 +234,8 @@ AddPolyData( vtkPolyData* pd, double opacity ) this->m_Mappers[ i ]->SetInputData( pd ); this->m_Actors[ i ]->SetMapper( this->m_Mappers[ i ] ); this->m_Actors[ i ]->GetProperty( )->SetOpacity( opacity ); + this->m_Actors[ i ]->GetProperty( )->SetLineWidth( 3 ); + this->m_Actors[ i ]->GetProperty( )->SetPointSize( 10 ); this->m_Renderer->AddActor( this->m_Actors[ i ] ); } @@ -252,7 +254,8 @@ AddPolyData( vtkPolyData* pd, double r, double g, double b, double opacity ) this->m_Actors[ i ]->SetMapper( this->m_Mappers[ i ] ); this->m_Actors[ i ]->GetProperty( )->SetColor( r, g, b ); this->m_Actors[ i ]->GetProperty( )->SetOpacity( opacity ); - this->m_Actors[ i ]->GetProperty( )->SetPointSize( 20 ); + this->m_Actors[ i ]->GetProperty( )->SetLineWidth( 3 ); + this->m_Actors[ i ]->GetProperty( )->SetPointSize( 10 ); this->m_Renderer->AddActor( this->m_Actors[ i ] ); } -- 2.45.1