example_Image_Dijkstra_DanielssonCost
example_Image_Dijkstra_DanielssonCost_TwoSeedsPath
example_Image_Dijkstra_EndPointDetection
+ example_Image_Dijkstra_LabelSkeleton
example_ShowSkeleton
)
FOREACH(EX ${SIMPLE_VTK_EXAMPLES})
--- /dev/null
+#include <cstdlib>
+#include <ctime>
+#include <iostream>
+#include <limits>
+#include <string>
+
+#include <itkImage.h>
+#include <itkImageToVTKImageFilter.h>
+
+#include <itkImageFileReader.h>
+#include <itkImageFileWriter.h>
+
+#include <fpa/Image/Functors/ImageCostFunction.h>
+#include <fpa/Image/DijkstraWithEndPointDetection.h>
+#include <fpa/VTK/ImageMPR.h>
+#include <fpa/VTK/Image3DObserver.h>
+
+#include <fpa/IO/MinimumSpanningTreeWriter.h>
+#include <fpa/IO/UniqueValuesContainerWriter.h>
+#include <fpa/IO/MatrixValuesContainerWriter.h>
+
+#include <vtkCellArray.h>
+#include <vtkFloatArray.h>
+#include <vtkImageMarchingCubes.h>
+#include <vtkPoints.h>
+#include <vtkPolyData.h>
+#include <vtkSmartPointer.h>
+
+// -------------------------------------------------------------------------
+const unsigned int Dim = 3;
+typedef unsigned char TPixel;
+typedef float TScalar;
+
+typedef itk::Image< TPixel, Dim > TImage;
+typedef itk::Image< TScalar, Dim > TScalarImage;
+typedef itk::ImageToVTKImageFilter< TImage > TVTKInputImage;
+
+// -------------------------------------------------------------------------
+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 F >
+class SkeletonCostFunction
+ : public fpa::Image::Functors::ImageCostFunction< typename F::TInputImage, typename F::TResult >
+{
+public:
+ typedef fpa::Image::Functors::ImageCostFunction< typename F::TInputImage, typename F::TResult > Superclass;
+ typedef SkeletonCostFunction Self;
+ typedef itk::SmartPointer< Self > Pointer;
+ typedef itk::SmartPointer< const Self > ConstPointer;
+
+ typedef typename Superclass::TInputImage TInputImage;
+ typedef typename Superclass::TResult TResult;
+ typedef typename Superclass::TIndex TIndex;
+
+public:
+ itkNewMacro( Self );
+ itkTypeMacro( SkeletonCostFunction, fpa_Image_Functors_ImageCostFunction );
+
+public:
+ virtual void SetInputImage( const TInputImage* img )
+ {
+ this->Superclass::SetInputImage( img );
+
+ typename TInputImage::SpacingType spac = img->GetSpacing( );
+ typename TInputImage::SpacingType::ValueType min_spac = spac[ 0 ];
+ for( unsigned int d = 1; d < TInputImage::ImageDimension; ++d )
+ min_spac = ( spac[ d ] < min_spac )? spac[ d ]: min_spac;
+ this->m_Distance = TResult( min_spac );
+ }
+
+ virtual TResult Evaluate( const TIndex& v, const TIndex& p ) const
+ {
+ TResult res = this->Superclass::Evaluate( v, p );
+ if( res > TResult( 0 ) )
+ return( this->m_Distance );
+ else
+ return( TResult( -1 ) );
+ }
+
+protected:
+ SkeletonCostFunction( )
+ : Superclass( )
+ { }
+ virtual ~SkeletonCostFunction( )
+ { }
+
+private:
+ // Purposely not implemented
+ SkeletonCostFunction( const Self& );
+ void operator=( const Self& );
+
+protected:
+ TResult m_Distance;
+};
+
+// -------------------------------------------------------------------------
+int main( int argc, char* argv[] )
+{
+ if( argc < 11 )
+ {
+ std::cerr
+ << "Usage: " << argv[ 0 ] << std::endl
+ << " input_image" << std::endl
+ << " seed_x seed_y seed_z" << 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 ];
+ std::string output_costmap_fn = argv[ 5 ];
+ std::string output_labels_fn = argv[ 6 ];
+ std::string mst_output_fn = argv[ 7 ];
+ std::string endpoints_output_fn = argv[ 8 ];
+ std::string bifurcations_output_fn = argv[ 9 ];
+ std::string branches_output_fn = argv[ 10 ];
+
+ // Get seed
+ TImage::PointType pnt;
+ pnt[ 0 ] = TImage::PointType::ValueType( std::atof( argv[ 2 ] ) );
+ pnt[ 1 ] = TImage::PointType::ValueType( std::atof( argv[ 3 ] ) );
+ pnt[ 2 ] = TImage::PointType::ValueType( std::atof( argv[ 4 ] ) );
+
+ // 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
+
+ TVTKInputImage::Pointer vtk_input_image = TVTKInputImage::New( );
+ vtk_input_image->SetInput( input_image );
+ vtk_input_image->Update( );
+
+ // Show input image and let some interaction
+ fpa::VTK::ImageMPR view;
+ view.SetBackground( 0.3, 0.2, 0.8 );
+ view.SetSize( 500, 500 );
+ view.SetImage( vtk_input_image->GetOutput( ) );
+
+ vtkSmartPointer< vtkImageMarchingCubes > mc =
+ vtkSmartPointer< vtkImageMarchingCubes >::New( );
+ mc->SetInputData( vtk_input_image->GetOutput( ) );
+ mc->SetValue( 0, 1e-1 );
+ mc->Update( );
+ view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 );
+
+ // Allow some interaction and wait for, at least, one seed
+ view.Render( );
+ view.Start( );
+
+ // Transform seed
+ TImage::IndexType seed;
+ input_image->TransformPhysicalPointToIndex( pnt, seed );
+
+ // Prepare Dijkstra filter
+ typedef fpa::Image::DijkstraWithEndPointDetection< TImage, TScalarImage > TFilter;
+ typedef SkeletonCostFunction< TFilter > TFunction;
+
+ TFunction::Pointer function = TFunction::New( );
+
+ TFilter::Pointer filter = TFilter::New( );
+ filter->SetInput( input_image );
+ filter->SetCostFunction( function );
+ filter->SetNeighborhoodOrder( 2 );
+ filter->SetSafetyNeighborhoodSize( 0 );
+ filter->StopAtOneFrontOff( );
+ filter->CorrectSeedsOff( );
+ filter->CorrectEndPointsOff( );
+ filter->AddSeed( seed, TScalar( 0 ) );
+
+ // Prepare graphical debugger
+ typedef fpa::VTK::Image3DObserver< TFilter, vtkRenderWindow > TDebugger;
+ TDebugger::Pointer debugger = TDebugger::New( );
+ debugger->SetRenderWindow( view.GetWindow( ) );
+ debugger->SetRenderPercentage( 0.0000001 );
+ filter->AddObserver( itk::AnyEvent( ), debugger );
+ filter->ThrowEventsOn( );
+
+ // 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( 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( );
+ for(
+ TFilter::TUniqueVertices::ConstIterator epIt = endpoints->Begin( );
+ epIt != endpoints->End( );
+ ++epIt
+ )
+ {
+ 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
+ vtkSmartPointer< vtkPolyData > endpoints_polydata =
+ vtkSmartPointer< vtkPolyData >::New( );
+ endpoints_polydata->SetPoints( endpoints_points );
+ endpoints_polydata->SetVerts( endpoints_cells );
+ view.AddPolyData( endpoints_polydata, 0, 0, 1, 1 );
+
+ // Show bifurcations
+ vtkSmartPointer< vtkPoints > bifurcations_points =
+ vtkSmartPointer< vtkPoints >::New( );
+ vtkSmartPointer< vtkCellArray > bifurcations_cells =
+ vtkSmartPointer< vtkCellArray >::New( );
+ for(
+ TFilter::TUniqueVertices::ConstIterator bfIt = bifurcations->Begin( );
+ bfIt != bifurcations->End( );
+ ++bfIt
+ )
+ {
+ TImage::PointType pnt;
+ input_image->TransformIndexToPhysicalPoint( *bfIt, pnt );
+ bifurcations_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
+ 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_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( );
+
+ vtkSmartPointer< vtkPoints > detailed_branches_points =
+ vtkSmartPointer< vtkPoints >::New( );
+ vtkSmartPointer< vtkCellArray > detailed_branches_cells =
+ vtkSmartPointer< vtkCellArray >::New( );
+ vtkSmartPointer< vtkFloatArray > detailed_branches_scalars =
+ vtkSmartPointer< vtkFloatArray >::New( );
+
+ TFilter::TBranches::ConstIterator brIt = branches->Begin( );
+ for( ; brIt != branches->End( ); ++brIt )
+ {
+ // 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 )
+ {
+ // 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 > detailed_branches_polydata =
+ vtkSmartPointer< vtkPolyData >::New( );
+ 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 );
+
+ // Let some more interaction
+ view.Render( );
+ view.Start( );
+
+ 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( );
+ image = reader->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
+}
+
+// eof - $RCSfile$
itkNewMacro( Self );
itkTypeMacro( DijkstraWithEndPointDetection, Dijkstra );
+ itkBooleanMacro( CorrectSeeds );
+ itkBooleanMacro( CorrectEndPoints );
+
+ itkGetConstMacro( CorrectSeeds, bool );
+ itkGetConstMacro( CorrectEndPoints, bool );
+ itkGetConstMacro( SafetyNeighborhoodSize, unsigned int );
+
itkGetConstMacro( NumberOfBranches, unsigned long );
+ itkSetMacro( CorrectSeeds, bool );
+ itkSetMacro( CorrectEndPoints, bool );
+ itkSetMacro( SafetyNeighborhoodSize, unsigned int );
+
public:
TLabelImage* GetLabelImage( );
const TLabelImage* GetLabelImage( ) const;
unsigned int m_EndPointsIndex;
unsigned int m_BranchesIndex;
- _TCandidates m_Candidates;
- unsigned long m_NumberOfBranches;
+ bool m_CorrectSeeds;
+ bool m_CorrectEndPoints;
+ unsigned int m_SafetyNeighborhoodSize;
+
+ _TCandidates m_Candidates;
+ unsigned long m_NumberOfBranches;
};
} // ecapseman
template< class I, class O >
fpa::Image::DijkstraWithEndPointDetection< I, O >::
DijkstraWithEndPointDetection( )
- : Superclass( )
+ : Superclass( ),
+ m_CorrectSeeds( true ),
+ m_CorrectEndPoints( true ),
+ m_SafetyNeighborhoodSize( 3 )
{
this->m_LabelImageIndex = this->GetNumberOfRequiredOutputs( );
this->m_EndPointsIndex = this->m_LabelImageIndex + 1;
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
- for( unsigned int s = 0; s < this->m_SeedVertices.size( ); ++s )
+ if( this->m_CorrectSeeds )
{
- 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 )
+ 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
+ for( unsigned int s = 0; s < this->m_SeedVertices.size( ); ++s )
{
- if( iIt.Get( ) > max_value )
+ 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 )
{
- this->m_SeedVertices[ s ] = iIt.GetIndex( );
- max_value = iIt.Get( );
+ if( iIt.Get( ) > max_value )
+ {
+ this->m_SeedVertices[ s ] = iIt.GetIndex( );
+ max_value = iIt.Get( );
- } // fi
+ } // fi
+
+ } // rof
+ this->m_Seeds.erase( seed );
+ n.Parent = this->m_SeedVertices[ s ];
+ this->m_Seeds[ this->m_SeedVertices[ s ] ] = n;
} // rof
- this->m_Seeds.erase( seed );
- n.Parent = this->m_SeedVertices[ s ];
- this->m_Seeds[ this->m_SeedVertices[ s ] ] = n;
- } // rof
+ } // fi
// End initialization
this->Superclass::_BeforeGenerateData( );
// 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 ) );
+ if( this->m_CorrectEndPoints )
+ v = this->_MaxInRegion( input, v, vr * double( 1.5 ) );
// Now, check for real marking conditions
// 1. Has it been visited by dijkstra?
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;
+ long s =
+ long( std::ceil( r / double( spac[ d ] ) ) ) +
+ long( this->m_SafetyNeighborhoodSize );
i0[ d ] = c[ d ] - s;
i1[ d ] = c[ d ] + s;
+++ /dev/null
-#ifndef __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__H__
-#define __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__H__
-
-#include <map>
-#include <deque>
-#include <utility>
-#include <fpa/Image/Dijkstra.h>
-
-namespace fpa
-{
- namespace Image
- {
- /**
- * @param I Input image type
- */
- template< class I, class C >
- class DijkstraWithSphereBacktracking
- : public fpa::Image::Dijkstra< I, C >
- {
- public:
- typedef DijkstraWithSphereBacktracking Self;
- typedef fpa::Image::Dijkstra< I, C > Superclass;
- typedef itk::SmartPointer< Self > Pointer;
- typedef itk::SmartPointer< const Self > ConstPointer;
-
- typedef typename Superclass::TCost TCost;
- typedef typename Superclass::TVertex TVertex;
- typedef typename Superclass::InputImageType TImage;
- typedef std::deque< TVertex > TVertices;
-
- typedef unsigned short TMark;
- typedef itk::Image< TMark, I::ImageDimension > TMarkImage;
-
- typedef typename Superclass::TTraits::TVertexCmp TVertexCmp;
- typedef std::pair< TVertex, TMark > TTreeNode;
- typedef std::map< TVertex, TTreeNode, TVertexCmp > TTree;
-
- typedef typename Superclass::TEndEvent TEndEvent;
- typedef typename Superclass::TBacktrackingEvent TBacktrackingEvent;
- typedef typename Superclass::TEndBacktrackingEvent TEndBacktrackingEvent;
-
- protected:
- typedef std::pair< TCost, TVertex > _TCandidate;
- typedef std::multimap< TCost, TVertex > _TCandidates;
- typedef typename Superclass::_TNode _TNode;
-
- typedef typename I::PixelType _TPixel;
- typedef typename I::RegionType _TRegion;
- typedef typename I::SizeType _TSize;
-
- public:
- itkNewMacro( Self );
- itkTypeMacro( DijkstraWithSphereBacktracking, Dijkstra );
-
- itkGetConstMacro( FullTree, TTree );
- itkGetConstMacro( ReducedTree, TTree );
- itkGetConstMacro( EndPoints, TVertices );
- itkGetConstMacro( BifurcationPoints, TVertices );
- itkGetConstMacro( NumberOfBranches, TMark );
-
- public:
- TMarkImage* GetOutputMarkImage( );
- const TMarkImage* GetOutputMarkImage( ) const;
-
- protected:
- DijkstraWithSphereBacktracking( );
- virtual ~DijkstraWithSphereBacktracking( );
-
- virtual void _BeforeMainLoop( );
- virtual void _AfterMainLoop( );
- virtual bool _UpdateNeigh( _TNode& nn, const _TNode& n );
- virtual bool _UpdateResult( _TNode& n );
-
- _TRegion _Region( const TVertex& c, const double& r );
-
- private:
- DijkstraWithSphereBacktracking( const Self& other );
- Self& operator=( const Self& other );
-
- protected:
- _TCandidates m_Candidates;
- TTree m_FullTree;
- TTree m_ReducedTree;
- TVertices m_BifurcationPoints;
- TVertices m_EndPoints;
- TMark m_NumberOfBranches;
- };
-
- } // ecapseman
-
-} // ecapseman
-
-#include <fpa/Image/DijkstraWithSphereBacktracking.hxx>
-
-#endif // __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__H__
-
-// eof - $RCSfile$
+++ /dev/null
-#ifndef __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__HXX__
-#define __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__HXX__
-
-#include <cmath>
-#include <itkImageRegionConstIteratorWithIndex.h>
-#include <itkImageRegionIteratorWithIndex.h>
-
-// -------------------------------------------------------------------------
-template< class I, class C >
-typename fpa::Image::DijkstraWithSphereBacktracking< I, C >::
-TMarkImage* fpa::Image::DijkstraWithSphereBacktracking< I, C >::
-GetOutputMarkImage( )
-{
- return(
- dynamic_cast< TMarkImage* >(
- this->itk::ProcessObject::GetOutput( 1 )
- )
- );
-}
-
-// -------------------------------------------------------------------------
-template< class I, class C >
-const typename fpa::Image::DijkstraWithSphereBacktracking< I, C >::
-TMarkImage* fpa::Image::DijkstraWithSphereBacktracking< I, C >::
-GetOutputMarkImage( ) const
-{
- return(
- dynamic_cast< const TMarkImage* >(
- this->itk::ProcessObject::GetOutput( 1 )
- )
- );
-}
-
-// -------------------------------------------------------------------------
-template< class I, class C >
-fpa::Image::DijkstraWithSphereBacktracking< I, C >::
-DijkstraWithSphereBacktracking( )
- : Superclass( ),
- m_NumberOfBranches( 0 )
-{
- this->SetNumberOfRequiredOutputs( 2 );
- this->SetNthOutput( 0, I::New( ) );
- this->SetNthOutput( 1, TMarkImage::New( ) );
-}
-
-// -------------------------------------------------------------------------
-template< class I, class C >
-fpa::Image::DijkstraWithSphereBacktracking< I, C >::
-~DijkstraWithSphereBacktracking( )
-{
-}
-
-// -------------------------------------------------------------------------
-template< class I, class C >
-void fpa::Image::DijkstraWithSphereBacktracking< I, C >::
-_BeforeMainLoop( )
-{
- const I* img = this->GetInput( );
- typename I::SpacingType spac = img->GetSpacing( );
- double max_spac = spac[ 0 ];
- for( unsigned int d = 1; d < I::ImageDimension; ++d )
- max_spac =
- ( max_spac < double( spac[ d ] ) )? double( spac[ d ] ): max_spac;
- max_spac *= double( 30 );
-
- // Correct seeds
- 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( );
-
- } // fi
-
- } // rof
- this->m_Seeds[ s ] = seed;
-
- } // rof
-
- // End initialization
- this->Superclass::_BeforeMainLoop( );
- this->m_Candidates.clear( );
-}
-
-// -------------------------------------------------------------------------
-template< class I, class C >
-void fpa::Image::DijkstraWithSphereBacktracking< I, C >::
-_AfterMainLoop( )
-{
- // Finish base algorithm
- this->Superclass::_AfterMainLoop( );
- this->m_FullTree.clear( );
- this->m_ReducedTree.clear( );
- this->m_EndPoints.clear( );
- this->m_BifurcationPoints.clear( );
- if( this->m_Candidates.size( ) == 0 )
- return;
- this->InvokeEvent( TEndEvent( ) );
-
- // Get some input values
- const I* input = this->GetInput( );
- typename I::SpacingType spac = input->GetSpacing( );
- double max_spac = spac[ 0 ];
- for( unsigned int d = 1; d < I::ImageDimension; ++d )
- max_spac =
- ( max_spac < double( spac[ d ] ) )? double( spac[ d ] ): max_spac;
- max_spac *= double( 3 );
-
- // Prepare an object to hold marks
- typename TMarkImage::Pointer marks = this->GetOutputMarkImage( );
- marks->FillBuffer( 0 );
-
- // Iterate over the candidates, starting from the candidates that
- // are near thin branches
- typename _TCandidates::const_reverse_iterator cIt =
- this->m_Candidates.rbegin( );
- this->m_NumberOfBranches = 1;
- for( ; cIt != this->m_Candidates.rend( ); ++cIt )
- {
- // If pixel hasn't been visited, continue
- TVertex v = cIt->second;
- if( marks->GetPixel( v ) != 0 )
- continue;
-
- // Compute nearest start candidate
- _TRegion region = this->_Region( v, max_spac );
- itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
- iIt.GoToBegin( );
- TVertex max_vertex = iIt.GetIndex( );
- _TPixel max_value = iIt.Get( );
- for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
- {
- _TPixel value = iIt.Get( );
- if( value > max_value )
- {
- max_value = value;
- max_vertex = iIt.GetIndex( );
-
- } // fi
-
- } // rof
-
- // Keep endpoint
- if( marks->GetPixel( max_vertex ) != 0 )
- continue;
-
- this->m_EndPoints.push_back( max_vertex );
-
- // Construct path (at least the part that hasn't been iterated)
- bool start = true;
- bool change = false;
- do
- {
- if( start )
- {
- if( this->m_FullTree.find( max_vertex ) == this->m_FullTree.end( ) )
- {
- // Mark a sphere around current point as visited
- double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) );
- region = this->_Region( max_vertex, dist * double( 1.5 ) );
- itk::ImageRegionIteratorWithIndex< TMarkImage >
- mIt( marks, region );
- for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
- mIt.Set( this->m_NumberOfBranches );
-
- // Next vertex in current path
- this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) );
- this->m_FullTree[ max_vertex ] =
- TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
- // std::cout << "New: " << this->m_NumberOfBranches << std::endl;
- }
- else
- {
- // A bifurcation point has been reached!
- this->m_BifurcationPoints.push_back( max_vertex );
- this->m_NumberOfBranches++;
-
- this->m_FullTree[ max_vertex ] =
- TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
- // std::cout << "Bifurcation: " << this->m_NumberOfBranches << std::endl;
-
- start = false;
-
- } // fi
- }
- else
- {
- if(
- std::find(
- this->m_BifurcationPoints.begin( ),
- this->m_BifurcationPoints.end( ),
- max_vertex
- ) == this->m_BifurcationPoints.end( )
- )
- {
- // Mark a sphere around current point as visited
- double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) );
- region = this->_Region( max_vertex, dist * double( 1.5 ) );
- itk::ImageRegionIteratorWithIndex< TMarkImage >
- mIt( marks, region );
- for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
- mIt.Set( this->m_NumberOfBranches );
-
- // Next vertex in current path
- this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) );
- this->m_FullTree[ max_vertex ] =
- TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
- change = true;
- // std::cout << "Change: " << this->m_NumberOfBranches << std::endl;
- }
- else
- {
- // A bifurcation point has been reached!
- // TODO: this->m_BifurcationPoints.push_back( max_vertex );
- this->m_NumberOfBranches++;
- this->m_FullTree[ max_vertex ] =
- TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
- // std::cout << "Change bifurcation: " << this->m_NumberOfBranches << std::endl;
-
- } // fi
-
- } // fi
- max_vertex = this->_Parent( max_vertex );
-
- } while( max_vertex != this->_Parent( max_vertex ) );
- if( start || change )
- this->m_NumberOfBranches++;
- this->m_FullTree[ max_vertex ] = TTreeNode( max_vertex, this->m_NumberOfBranches );
-
- this->InvokeEvent( TEndBacktrackingEvent( ) );
-
- } // rof
-
- std::map< TMark, unsigned long > histo;
- for(
- typename TTree::iterator treeIt = this->m_FullTree.begin( );
- treeIt != this->m_FullTree.end( );
- ++treeIt
- )
- histo[ treeIt->second.second ]++;
-
- std::map< TMark, TMark > changes;
- TMark last_change = 1;
- for( TMark i = 1; i <= this->m_NumberOfBranches; ++i )
- {
- if( histo[ i ] != 0 )
- changes[ i ] = last_change++;
-
- } // rof
- this->m_NumberOfBranches = changes.size( );
-
- 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 ];
-
- } // fi
-
- // 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
- {
- tIt = this->m_FullTree.find( tIt->second.first );
- if( tIt == this->m_FullTree.end( ) )
- break;
-
- } while( tIt->second.second == id );
- this->m_ReducedTree[ *eIt ] = TTreeNode( tIt->first, id );
-
- } // fi
-
- } // 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;
-
- } while( tIt->second.second == id );
- this->m_ReducedTree[ *bIt ] = TTreeNode( tIt->first, id );
-
- } // fi
-
- } // rof
-}
-
-// -------------------------------------------------------------------------
-template< class I, class C >
-bool fpa::Image::DijkstraWithSphereBacktracking< I, C >::
-_UpdateNeigh( _TNode& nn, const _TNode& n )
-{
- C nc = this->_Cost( nn.Vertex, n.Vertex );
- if( TCost( 0 ) < nc )
- {
- typename I::PointType pnn, pn;
- this->GetInput( )->TransformIndexToPhysicalPoint( nn.Vertex, pnn );
- this->GetInput( )->TransformIndexToPhysicalPoint( n.Vertex, pn );
-
-
- nc += TCost( 1 );
- nn.Cost = n.Cost + ( TCost( pnn.EuclideanDistanceTo( pn ) ) / std::pow( nc, 4 ) );
- nn.Result = nn.Cost;
- return( true );
- }
- else
- return( false );
-}
-
-// -------------------------------------------------------------------------
-template< class I, class C >
-bool fpa::Image::DijkstraWithSphereBacktracking< I, C >::
-_UpdateResult( _TNode& n )
-{
- bool ret = this->Superclass::_UpdateResult( n );
- if( ret )
- {
- TCost nc = this->_Cost( n.Vertex, n.Parent );
- if( TCost( 0 ) < nc )
- {
- TCost cc = n.Cost / nc;
- this->m_Candidates.insert( _TCandidate( cc, n.Vertex ) );
- /* TODO: debug code
- this->GetOutput( )->SetPixel( n.Vertex, cc );
- */
- } // fi
-
- } // fi
- return( ret );
-}
-
-// -------------------------------------------------------------------------
-template< class I, class C >
-typename fpa::Image::DijkstraWithSphereBacktracking< I, C >::
-_TRegion fpa::Image::DijkstraWithSphereBacktracking< I, C >::
-_Region( const TVertex& c, const double& r )
-{
- typename I::ConstPointer input = this->GetInput( );
- typename I::SpacingType spac = input->GetSpacing( );
- _TRegion region = input->GetLargestPossibleRegion( );
- typename I::IndexType idx0 = region.GetIndex( );
- typename I::IndexType idx1 = idx0 + region.GetSize( );
-
- // Compute region size and index
- typename I::IndexType i0, i1;
- _TSize size;
- for( unsigned int d = 0; d < I::ImageDimension; ++d )
- {
- long s = long( std::ceil( r / double( spac[ d ] ) ) );
- i0[ d ] = c[ d ] - s;
- i1[ d ] = c[ d ] + s;
-
- if( i0[ d ] < idx0[ d ] ) i0[ d ] = idx0[ d ];
- if( i1[ d ] < idx0[ d ] ) i1[ d ] = idx0[ d ];
- if( i0[ d ] > idx1[ d ] ) i0[ d ] = idx1[ d ];
- if( i1[ d ] > idx1[ d ] ) i1[ d ] = idx1[ d ];
- size[ d ] = i1[ d ] - i0[ d ];
-
- } // rof
-
- // Prepare region and return it
- region.SetIndex( i0 );
- region.SetSize( size );
- return( region );
-}
-
-#endif // __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__HXX__
-
-// eof - $RCSfile$
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / commitdiff