From: Leonardo Florez-Valencia Date: Mon, 23 Feb 2015 02:38:22 +0000 (-0500) Subject: Complete skeletonization application added X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=73e3158dab6ffe30c215b2899fd80ee868f330ce;p=FrontAlgorithms.git Complete skeletonization application added --- diff --git a/appli/examples/CMakeLists.txt b/appli/examples/CMakeLists.txt index dff8fe0..520a4a3 100644 --- a/appli/examples/CMakeLists.txt +++ b/appli/examples/CMakeLists.txt @@ -19,6 +19,7 @@ IF(BUILD_EXAMPLES) example_ImageAlgorithmDijkstra_01 example_ImageAlgorithmDijkstra_02 example_ImageAlgorithmFastMarching_01 + example_ImageAlgorithm_Skeletonization ) FOREACH(APP ${vtk_APPLIS}) diff --git a/appli/examples/example_ImageAlgorithmRegionGrow_MultipleThresholds.cxx b/appli/examples/example_ImageAlgorithmRegionGrow_MultipleThresholds.cxx index e9b06cc..0216bed 100644 --- a/appli/examples/example_ImageAlgorithmRegionGrow_MultipleThresholds.cxx +++ b/appli/examples/example_ImageAlgorithmRegionGrow_MultipleThresholds.cxx @@ -110,6 +110,8 @@ int main( int argc, char* argv[] ) algorithm->SetInput( input_image ); algorithm->SetNeighborhoodOrder( 1 ); algorithm->SetDifferenceThreshold( double( 3 ) ); + algorithm->SetInsideValue( 0 ); + algorithm->SetOutsideValue( 1 ); if( visual_debug ) { diff --git a/appli/examples/example_ImageAlgorithm_Skeletonization.cxx b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx new file mode 100644 index 0000000..156469e --- /dev/null +++ b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx @@ -0,0 +1,222 @@ +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +/* + #include + #include + + #include + + + #include + #include + #include + #include + #include + #include + #include + #include + +*/ + +// ------------------------------------------------------------------------- +const unsigned int Dim = 3; +typedef short TPixel; +typedef double TScalar; +typedef itk::Image< TPixel, Dim > TImage; +typedef itk::Image< TScalar, Dim > TDistanceMap; + +/* + typedef itk::ImageToVTKImageFilter< TImage > TVTKImage; + typedef itk::ImageFileReader< TImage > TImageReader; + typedef itk::ImageFileWriter< TImage > TImageWriter; + typedef + fpa::Image::RegionGrowWithMultipleThresholds< TImage > + TSegmentor; + typedef + + TObserver; +*/ + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + if( argc < 6 ) + { + std::cerr + << "Usage: " << argv[ 0 ] + << " input_image thr_0 thr_1 step output_image" + << " visual_debug" + << std::endl; + return( 1 ); + + } // fi + std::string input_image_fn = argv[ 1 ]; + TPixel thr_0 = TPixel( std::atof( argv[ 2 ] ) ); + TPixel thr_1 = TPixel( std::atof( argv[ 3 ] ) ); + unsigned int step = std::atoi( argv[ 4 ] ); + std::string output_image_fn = argv[ 5 ]; + bool visual_debug = false; + if( argc > 6 ) + visual_debug = ( std::atoi( argv[ 6 ] ) == 1 ); + + // Read image + itk::ImageFileReader< TImage >::Pointer input_image_reader = + itk::ImageFileReader< TImage >::New( ); + input_image_reader->SetFileName( input_image_fn ); + try + { + input_image_reader->Update( ); + } + catch( itk::ExceptionObject& err ) + { + std::cerr << "Error caught: " << err << std::endl; + return( 1 ); + + } // yrt + TImage::ConstPointer input_image = input_image_reader->GetOutput( ); + + // Show image + itk::ImageToVTKImageFilter< TImage >::Pointer vtk_image = + itk::ImageToVTKImageFilter< TImage >::New( ); + vtk_image->SetInput( input_image ); + vtk_image->Update( ); + + fpa::VTK::ImageMPR view; + view.SetBackground( 0.3, 0.2, 0.8 ); + view.SetSize( 800, 800 ); + view.SetImage( vtk_image->GetOutput( ) ); + + // Wait for a seed to be given + while( view.GetNumberOfSeeds( ) == 0 ) + view.Start( ); + + // Get seed from user + double seed[ 3 ]; + view.GetSeed( 0, seed ); + TImage::PointType seed_pnt; + seed_pnt[ 0 ] = seed[ 0 ]; + seed_pnt[ 1 ] = seed[ 1 ]; + seed_pnt[ 2 ] = seed[ 2 ]; + TImage::IndexType seed_idx; + input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx ); + + // Segment input image + fpa::Image::RegionGrowWithMultipleThresholds< TImage >::Pointer segmentor = + fpa::Image::RegionGrowWithMultipleThresholds< TImage >::New( ); + segmentor->AddThresholds( thr_0, thr_1, step ); + segmentor->AddSeed( seed_idx, 0 ); + segmentor->SetInput( input_image ); + segmentor->SetNeighborhoodOrder( 1 ); + segmentor->SetDifferenceThreshold( double( 3 ) ); + segmentor->SetInsideValue( TPixel( 0 ) ); // WARNING: This is to optimize + segmentor->SetOutsideValue( TPixel( 1 ) ); // distance map calculation + if( visual_debug ) + { + // Configure observer + fpa::VTK::Image3DObserver< fpa::Image::RegionGrowWithMultipleThresholds< TImage >, vtkRenderWindow >::Pointer obs = + fpa::VTK::Image3DObserver< fpa::Image::RegionGrowWithMultipleThresholds< TImage >, vtkRenderWindow >::New( ); + obs->SetImage( input_image, view.GetWindow( ) ); + segmentor->AddObserver( itk::AnyEvent( ), obs ); + segmentor->ThrowEventsOn( ); + } + else + segmentor->ThrowEventsOff( ); + std::clock_t start = std::clock( ); + segmentor->Update( ); + std::clock_t end = std::clock( ); + double seconds = double( end - start ) / double( CLOCKS_PER_SEC ); + std::cout << "Segmentation time = " << seconds << std::endl; + + // Extract distance map + itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap >::Pointer + distanceMap = + itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap >::New( ); + distanceMap->SetInput( segmentor->GetOutput( ) ); + distanceMap->InputIsBinaryOn( ); + start = std::clock( ); + distanceMap->Update( ); + end = std::clock( ); + seconds = double( end - start ) / double( CLOCKS_PER_SEC ); + std::cout << "Distance map time = " << seconds << std::endl; + + // Extract paths + fpa::Image::Dijkstra< TDistanceMap, TScalar >::Pointer paths = + fpa::Image::Dijkstra< TDistanceMap, TScalar >::New( ); + paths->AddSeed( segmentor->GetSeed( 0 ), TScalar( 0 ) ); + paths->SetInput( distanceMap->GetOutput( ) ); + paths->SetNeighborhoodOrder( 1 ); + + // Associate an inversion cost function + fpa::Base::Functors::InvertCostFunction< TScalar >::Pointer + cost_function = + fpa::Base::Functors::InvertCostFunction< TScalar >::New( ); + paths->SetCostConversion( cost_function ); + + if( visual_debug ) + { + // Configure observer + fpa::VTK::Image3DObserver< fpa::Image::Dijkstra< TDistanceMap, TScalar >, vtkRenderWindow >::Pointer obs = + fpa::VTK::Image3DObserver< fpa::Image::Dijkstra< TDistanceMap, TScalar >, vtkRenderWindow >::New( ); + obs->SetImage( distanceMap->GetOutput( ), view.GetWindow( ) ); + paths->AddObserver( itk::AnyEvent( ), obs ); + paths->ThrowEventsOn( ); + } + else + paths->ThrowEventsOff( ); + start = std::clock( ); + paths->Update( ); + end = std::clock( ); + seconds = double( end - start ) / double( CLOCKS_PER_SEC ); + std::cout << "Paths extraction time = " << seconds << std::endl; + + // Show result + /* + TVTKImage::Pointer output_image_vtk = TVTKImage::New( ); + output_image_vtk->SetInput( segmentor->GetOutput( ) ); + output_image_vtk->Update( ); + + vtkSmartPointer< vtkImageMarchingCubes > mc = + vtkSmartPointer< vtkImageMarchingCubes >::New( ); + mc->SetInputData( output_image_vtk->GetOutput( ) ); + mc->SetValue( + 0, + double( segmentor->GetInsideValue( ) ) * double( 0.95 ) + ); + mc->Update( ); + + // Let some interaction and close program + view.AddPolyData( mc->GetOutput( ), 0.1, 0.6, 0.8, 0.5 ); + view.Start( ); + + // Write resulting image + TImageWriter::Pointer output_image_writer = TImageWriter::New( ); + output_image_writer->SetInput( segmentor->GetOutput( ) ); + output_image_writer->SetFileName( output_image_fn ); + try + { + output_image_writer->Update( ); + } + catch( itk::ExceptionObject& err ) + { + std::cerr << "Error caught: " << err << std::endl; + return( 1 ); + + } // yrt + */ + return( 0 ); +} + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/Functors/InvertCostFunction.h b/lib/fpa/Base/Functors/InvertCostFunction.h index faf4fb2..278ae55 100644 --- a/lib/fpa/Base/Functors/InvertCostFunction.h +++ b/lib/fpa/Base/Functors/InvertCostFunction.h @@ -29,7 +29,7 @@ namespace fpa public: virtual C Evaluate( const C& input ) const { - return( C( 1 ) / ( C( 1 ) + C( input ) ) ); + return( C( 1 ) / ( C( 0 ) + C( input ) ) ); } protected: