From c7563ffe89c76b52736c3bd6bb2d5e6fac009cdd Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Mon, 23 Feb 2015 17:01:35 -0500 Subject: [PATCH] A dijkstra example added --- appli/examples/CMakeLists.txt | 1 + .../example_ImageAlgorithmDijkstra_03.cxx | 323 ++++++++++++++++++ ...example_ImageAlgorithm_Skeletonization.cxx | 112 ++---- .../RegionGrowWithMultipleThresholds.hxx | 1 + 4 files changed, 362 insertions(+), 75 deletions(-) create mode 100644 appli/examples/example_ImageAlgorithmDijkstra_03.cxx diff --git a/appli/examples/CMakeLists.txt b/appli/examples/CMakeLists.txt index 520a4a3..62f894b 100644 --- a/appli/examples/CMakeLists.txt +++ b/appli/examples/CMakeLists.txt @@ -18,6 +18,7 @@ IF(BUILD_EXAMPLES) example_ImageAlgorithmRegionGrow_MultipleThresholds example_ImageAlgorithmDijkstra_01 example_ImageAlgorithmDijkstra_02 + example_ImageAlgorithmDijkstra_03 example_ImageAlgorithmFastMarching_01 example_ImageAlgorithm_Skeletonization ) diff --git a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx new file mode 100644 index 0000000..46bfe03 --- /dev/null +++ b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx @@ -0,0 +1,323 @@ +#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 double TPixel; +typedef double TScalar; +typedef itk::Image< TPixel, Dim > TImage; +typedef itk::ImageToVTKImageFilter< TImage > TVTKImage; + +typedef itk::ImageFileReader< TImage > TImageReader; +typedef itk::ImageFileWriter< TImage > TImageWriter; +typedef fpa::Image::Dijkstra< TImage, TScalar > TBaseDijkstra; + +typedef fpa::VTK::ImageMPR TMPR; +typedef fpa::VTK::Image3DObserver< TBaseDijkstra, vtkRenderWindow > TDijkstraObs; + +// ------------------------------------------------------------------------- +class TDijkstra + : public TBaseDijkstra +{ +public: + typedef TDijkstra Self; + typedef TBaseDijkstra Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef Superclass::TCost TCost; + typedef Superclass::TVertex TVertex; + typedef Superclass::InputImageType TImage; + typedef std::deque< TVertex > TVertices; + +protected: + typedef std::pair< TCost, TVertex > _TCandidate; + typedef std::multimap< TCost, TVertex > _TCandidates; + typedef Superclass::_TNode _TNode; + +public: + itkNewMacro( Self ); + itkTypeMacro( TDijkstra, TBaseDijkstra ); + +protected: + TDijkstra( ) + : Superclass( ) + { } + virtual ~TDijkstra( ) + { } + + virtual void _BeforeMainLoop( ) + { + const TImage* img = this->GetInput( ); + + // Correct seeds + TImage::SizeType radius; + radius.Fill( 3 ); + itk::ConstNeighborhoodIterator< TImage > it( + radius, img, img->GetRequestedRegion( ) + ); + for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s ) + { + _TNode seed = this->m_Seeds[ s ]; + it.SetLocation( seed.Vertex ); + + TImage::SizeType size = it.GetSize( ); + unsigned int nN = 1; + for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) + nN *= size[ d ]; + TImage::PixelType max_value = img->GetPixel( seed.Vertex ); + for( unsigned int i = 0; i < nN; ++i ) + { + if( it.GetPixel( i ) > max_value ) + { + seed.Vertex = it.GetIndex( i ); + seed.Parent = seed.Vertex; + max_value = it.GetPixel( i ); + + } // fi + + } // rof + this->m_Seeds[ s ] = seed; + + } // rof + + // End initialization + this->Superclass::_BeforeMainLoop( ); + this->m_Candidates.clear( ); + } + + virtual void _AfterMainLoop( ) + { + this->Superclass::_AfterMainLoop( ); + if( this->m_Candidates.size( ) == 0 ) + return; + + this->m_Axes = vtkSmartPointer< vtkPolyData >::New( ); + vtkSmartPointer< vtkPoints > points = + vtkSmartPointer< vtkPoints >::New( ); + vtkSmartPointer< vtkCellArray > lines = + vtkSmartPointer< vtkCellArray >::New( ); + + _TCandidates::const_reverse_iterator cIt = + this->m_Candidates.rbegin( ); + + TVertices pS; + TVertex vIt = cIt->second; + TVertex parent = this->_Parent( vIt ); + while( parent != vIt ) + { + pS.push_front( vIt ); + vIt = parent; + parent = this->_Parent( vIt ); + + TImage::PointType pnt; + this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt ); + points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); + if( points->GetNumberOfPoints( ) > 1 ) + { + lines->InsertNextCell( 2 ); + lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 ); + lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 ); + + } // fi + + } // elihw + pS.push_front( vIt ); + TImage::PointType pnt; + this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt ); + points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); + if( points->GetNumberOfPoints( ) > 1 ) + { + lines->InsertNextCell( 2 ); + lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 ); + lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 ); + + } // fi + + this->m_Axes->SetPoints( points ); + this->m_Axes->SetLines( lines ); + } + + virtual bool _UpdateNeigh( _TNode& nn, const _TNode& n ) + { + TCost nc = this->_Cost( nn.Vertex, n.Vertex ); + if( TCost( 0 ) < nc ) + { + nn.Cost = n.Cost + ( TCost( 1 ) / nc ); + nn.Result = nn.Cost; + return( true ); + } + else + return( false ); + } + + virtual bool _UpdateResult( _TNode& n ) + { + bool ret = this->Superclass::_UpdateResult( n ); + if( ret ) + { + TCost nc = this->_Cost( n.Vertex, n.Parent ); + if( TCost( 0 ) < nc ) + { + n.Result = n.Cost / ( nc * nc ); + this->GetOutput( )->SetPixel( n.Vertex, n.Result ); + this->m_Candidates.insert( _TCandidate( n.Result, n.Vertex ) ); + + } // fi + + } // fi + return( ret ); + } + +private: + TDijkstra( const Self& other ); + Self& operator=( const Self& other ); + +protected: + _TCandidates m_Candidates; +public: + vtkSmartPointer< vtkPolyData > m_Axes; +}; + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + if( argc < 6 ) + { + std::cerr + << "Usage: " << argv[ 0 ] + << " input_image x y z output_image" + << " visual_debug" + << std::endl; + return( 1 ); + + } // fi + std::string input_image_fn = argv[ 1 ]; + TImage::IndexType seed_idx; + seed_idx[ 0 ] = std::atoi( argv[ 2 ] ); + seed_idx[ 1 ] = std::atoi( argv[ 3 ] ); + seed_idx[ 2 ] = 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 + TImageReader::Pointer input_image_reader = TImageReader::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 + TVTKImage::Pointer vtk_image = TVTKImage::New( ); + vtk_image->SetInput( input_image ); + vtk_image->Update( ); + + TMPR view; + view.SetBackground( 0.3, 0.2, 0.8 ); + view.SetSize( 800, 800 ); + view.SetImage( vtk_image->GetOutput( ) ); + + // Allow some interaction + view.Start( ); + + // Extract paths + TDijkstra::Pointer paths = TDijkstra::New( ); + paths->AddSeed( seed_idx, TScalar( 0 ) ); + paths->SetInput( input_image ); + paths->SetNeighborhoodOrder( 1 ); + + if( visual_debug ) + { + // Configure observer + TDijkstraObs::Pointer obs = TDijkstraObs::New( ); + obs->SetImage( input_image, view.GetWindow( ) ); + paths->AddObserver( itk::AnyEvent( ), obs ); + paths->ThrowEventsOn( ); + } + else + paths->ThrowEventsOff( ); + std::clock_t start = std::clock( ); + paths->Update( ); + std::clock_t end = std::clock( ); + double seconds = double( end - start ) / double( CLOCKS_PER_SEC ); + std::cout << "Paths extraction time = " << seconds << std::endl; + + view.AddPolyData( paths->m_Axes, 1, 0, 0 ); + view.Start( ); + + TImageWriter::Pointer dijkstra_writer = + TImageWriter::New( ); + dijkstra_writer->SetInput( paths->GetOutput( ) ); + dijkstra_writer->SetFileName( "dijkstra.mhd" ); + dijkstra_writer->Update( ); + + // 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/appli/examples/example_ImageAlgorithm_Skeletonization.cxx b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx index 156469e..f2b8279 100644 --- a/appli/examples/example_ImageAlgorithm_Skeletonization.cxx +++ b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx @@ -2,62 +2,45 @@ #include #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 itk::ImageFileWriter< TDistanceMap > TDistanceMapWriter; +typedef fpa::Image::RegionGrowWithMultipleThresholds< TImage > TSegmentor; +typedef +itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap > +TDistanceFilter; -/* - typedef itk::ImageToVTKImageFilter< TImage > TVTKImage; - typedef itk::ImageFileReader< TImage > TImageReader; - typedef itk::ImageFileWriter< TImage > TImageWriter; - typedef - fpa::Image::RegionGrowWithMultipleThresholds< TImage > - TSegmentor; - typedef - - TObserver; -*/ +typedef fpa::VTK::ImageMPR TMPR; +typedef fpa::VTK::Image3DObserver< TSegmentor, vtkRenderWindow > TSegmentorObs; // ------------------------------------------------------------------------- int main( int argc, char* argv[] ) { - if( argc < 6 ) + if( argc < 7 ) { std::cerr << "Usage: " << argv[ 0 ] - << " input_image thr_0 thr_1 step output_image" + << " input_image thr_0 thr_1 step" + << " output_segmentation output_distancemap" << " visual_debug" << std::endl; return( 1 ); @@ -67,14 +50,14 @@ int main( int argc, char* argv[] ) 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 ]; + std::string output_segmentation_fn = argv[ 5 ]; + std::string output_distancemap_fn = argv[ 6 ]; bool visual_debug = false; - if( argc > 6 ) - visual_debug = ( std::atoi( argv[ 6 ] ) == 1 ); + if( argc > 7 ) + visual_debug = ( std::atoi( argv[ 7 ] ) == 1 ); // Read image - itk::ImageFileReader< TImage >::Pointer input_image_reader = - itk::ImageFileReader< TImage >::New( ); + TImageReader::Pointer input_image_reader = TImageReader::New( ); input_image_reader->SetFileName( input_image_fn ); try { @@ -89,12 +72,11 @@ int main( int argc, char* argv[] ) TImage::ConstPointer input_image = input_image_reader->GetOutput( ); // Show image - itk::ImageToVTKImageFilter< TImage >::Pointer vtk_image = - itk::ImageToVTKImageFilter< TImage >::New( ); + TVTKImage::Pointer vtk_image = TVTKImage::New( ); vtk_image->SetInput( input_image ); vtk_image->Update( ); - fpa::VTK::ImageMPR view; + TMPR view; view.SetBackground( 0.3, 0.2, 0.8 ); view.SetSize( 800, 800 ); view.SetImage( vtk_image->GetOutput( ) ); @@ -114,8 +96,7 @@ int main( int argc, char* argv[] ) input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx ); // Segment input image - fpa::Image::RegionGrowWithMultipleThresholds< TImage >::Pointer segmentor = - fpa::Image::RegionGrowWithMultipleThresholds< TImage >::New( ); + TSegmentor::Pointer segmentor = TSegmentor::New( ); segmentor->AddThresholds( thr_0, thr_1, step ); segmentor->AddSeed( seed_idx, 0 ); segmentor->SetInput( input_image ); @@ -126,8 +107,7 @@ int main( int argc, char* argv[] ) if( visual_debug ) { // Configure observer - fpa::VTK::Image3DObserver< fpa::Image::RegionGrowWithMultipleThresholds< TImage >, vtkRenderWindow >::Pointer obs = - fpa::VTK::Image3DObserver< fpa::Image::RegionGrowWithMultipleThresholds< TImage >, vtkRenderWindow >::New( ); + TSegmentorObs::Pointer obs = TSegmentorObs::New( ); obs->SetImage( input_image, view.GetWindow( ) ); segmentor->AddObserver( itk::AnyEvent( ), obs ); segmentor->ThrowEventsOn( ); @@ -141,9 +121,7 @@ int main( int argc, char* argv[] ) std::cout << "Segmentation time = " << seconds << std::endl; // Extract distance map - itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap >::Pointer - distanceMap = - itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap >::New( ); + TDistanceFilter::Pointer distanceMap = TDistanceFilter::New( ); distanceMap->SetInput( segmentor->GetOutput( ) ); distanceMap->InputIsBinaryOn( ); start = std::clock( ); @@ -152,35 +130,19 @@ int main( int argc, char* argv[] ) 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 ); + std::cout << "Final seed: " << seed_idx << std::endl; - // Associate an inversion cost function - fpa::Base::Functors::InvertCostFunction< TScalar >::Pointer - cost_function = - fpa::Base::Functors::InvertCostFunction< TScalar >::New( ); - paths->SetCostConversion( cost_function ); + TDistanceMapWriter::Pointer distancemap_writer = + TDistanceMapWriter::New( ); + distancemap_writer->SetInput( distanceMap->GetOutput( ) ); + distancemap_writer->SetFileName( output_distancemap_fn ); + distancemap_writer->Update( ); - 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; + TImageWriter::Pointer segmentation_writer = + TImageWriter::New( ); + segmentation_writer->SetInput( segmentor->GetOutput( ) ); + segmentation_writer->SetFileName( output_segmentation_fn ); + segmentation_writer->Update( ); // Show result /* diff --git a/lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx b/lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx index e9ea456..2e8c995 100644 --- a/lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx +++ b/lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx @@ -93,6 +93,7 @@ _BeforeMainLoop( ) if( it.GetPixel( i ) < min_value ) { seed.Vertex = it.GetIndex( i ); + seed.Parent = seed.Vertex; min_value = it.GetPixel( i ); } // fi -- 2.47.1