#include <vtkPoints.h>
#include <vtkCellArray.h>
+#include <vtkFloatArray.h>
#include <vtkPolyData.h>
#include <vtkSmartPointer.h>
-#include <fpa/Base/Functors/InvertCostFunction.h>
-#include <fpa/Image/Dijkstra.h>
-#include <fpa/Image/RegionGrowWithMultipleThresholds.h>
+#include <fpa/Image/DijkstraWithSphereBacktracking.h>
#include <fpa/VTK/ImageMPR.h>
#include <fpa/VTK/Image3DObserver.h>
typedef itk::ImageFileReader< TImage > TImageReader;
typedef itk::ImageFileWriter< TImage > TImageWriter;
-typedef fpa::Image::Dijkstra< TImage, TScalar > TBaseDijkstra;
+typedef fpa::Image::DijkstraWithSphereBacktracking< TImage, TScalar > TDijkstra;
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;
-
- typedef Superclass::TEndEvent TEndEvent;
- typedef Superclass::TBacktrackingEvent TBacktrackingEvent;
-
-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( )
- {
- typedef unsigned char _TMark;
- typedef itk::Image< _TMark, TImage::ImageDimension > _TMarkImage;
-
- this->Superclass::_AfterMainLoop( );
- if( this->m_Candidates.size( ) == 0 )
- return;
-
- this->InvokeEvent( TEndEvent( ) );
-
- const TImage* input = this->GetInput( );
- TImage::SpacingType spacing = input->GetSpacing( );
-
- // Prepare an object to hold marks
- _TMarkImage::Pointer marks = _TMarkImage::New( );
- marks->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
- marks->SetRequestedRegion( input->GetRequestedRegion( ) );
- marks->SetBufferedRegion( input->GetBufferedRegion( ) );
- marks->SetOrigin( input->GetOrigin( ) );
- marks->SetSpacing( spacing );
- marks->SetDirection( input->GetDirection( ) );
- marks->Allocate( );
- marks->FillBuffer( _TMark( 0 ) );
-
- // Iterate over the candidates, starting from the candidates that
- // are near thin branches
- _TCandidates::const_reverse_iterator cIt =
- this->m_Candidates.rbegin( );
- for( int backId = 0; backId < 100 && cIt != this->m_Candidates.rend( ); ++cIt )
- {
- // If pixel hasn't been visited, continue
- TVertex v = cIt->second;
- if( marks->GetPixel( v ) != _TMark( 0 ) )
- continue;
-
- // Compute nearest start candidate
- TImage::SizeType radius;
- radius.Fill( 3 );
- itk::ConstNeighborhoodIterator< TImage > iIt(
- radius, input, input->GetRequestedRegion( )
- );
- iIt.SetLocation( v );
- unsigned int nN = 1;
- for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
- nN *= iIt.GetSize( )[ d ];
- TVertex max_vertex = iIt.GetIndex( 0 );
- TImage::PixelType max_value = iIt.GetPixel( 0 );
- for( unsigned int i = 1; i < nN; ++i )
- {
- TImage::PixelType value = iIt.GetPixel( i );
- if( value > max_value )
- {
- max_value = value;
- max_vertex = iIt.GetIndex( i );
-
- } // fi
-
- } // rof
-
- if( marks->GetPixel( max_vertex ) != _TMark( 0 ) )
- continue;
- backId++;
- std::cout << "Leaf: " << backId << " " << max_value << " " << max_vertex << std::endl;
-
- bool start = true;
- do
- {
- double dist = double( 1.5 ) * std::sqrt( double( input->GetPixel( max_vertex ) ) );
- for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
- radius[ d ] =
- ( unsigned long )( dist / spacing[ d ] ) + 1;
- itk::NeighborhoodIterator< _TMarkImage > mIt(
- radius, marks, marks->GetRequestedRegion( )
- );
- mIt.SetLocation( max_vertex );
- nN = 1;
- for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
- nN *= mIt.GetSize( )[ d ];
- for( unsigned int i = 0; i < nN; ++i )
- if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) )
- {
- mIt.SetPixel( i, ( start )? 255: 100 );
- start = false;
- }
-
- /*
- TImage::SizeType radius;
- mIt.GoToBegin( );
- mIt.SetLocation( vIt );
-
- TImage::SizeType size = mIt.GetSize( );
- unsigned int nN = 1;
- for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
- nN *= size[ d ];
- for( unsigned int i = 0; i < nN; ++i )
- if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) )
- mIt.SetPixel( i, ( start )? 255: 100 );
-
- start = false;
- */
- // Next vertex in current path
- this->InvokeEvent( TBacktrackingEvent( max_vertex, backId ) );
- max_vertex = this->_Parent( max_vertex );
-
- } while( max_vertex != this->_Parent( max_vertex ) );
-
- } // rof
- /*
- else
- marks->SetPixel( v, _TMark( 1 ) );
- } // rof
- */
-
- /*
- itk::ImageFileWriter< _TMarkImage >::Pointer w =
- itk::ImageFileWriter< _TMarkImage >::New( );
- w->SetInput( marks );
- w->SetFileName( "marks.mhd" );
- w->Update( );
-
-
- this->m_Axes = vtkSmartPointer< vtkPolyData >::New( );
- vtkSmartPointer< vtkPoints > points =
- vtkSmartPointer< vtkPoints >::New( );
- vtkSmartPointer< vtkCellArray > lines =
- vtkSmartPointer< vtkCellArray >::New( );
-
- {
-
- backId++;
- std::cout << "Leaf: " << backId << " " << cIt->first << " " << vIt << std::endl;
- bool start = true;
- do
- {
- double dist = double( input->GetPixel( vIt ) );
- TImage::SizeType radius;
- for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
- radius[ d ] =
- ( unsigned long )( dist / spacing[ d ] ) + 1;
- itk::NeighborhoodIterator< _TMarkImage > mIt(
- radius, marks, marks->GetRequestedRegion( )
- );
- mIt.GoToBegin( );
- mIt.SetLocation( vIt );
-
- TImage::SizeType size = mIt.GetSize( );
- unsigned int nN = 1;
- for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
- nN *= size[ d ];
- for( unsigned int i = 0; i < nN; ++i )
- if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) )
- mIt.SetPixel( i, ( start )? 255: 100 );
-
- start = false;
-
- // Next vertex
- vIt = this->_Parent( vIt );
-
- } while( vIt != this->_Parent( vIt ) );
-
- // Last vertex
- // std::cout << vIt << " " << this->GetSeed( 0 ) << std::endl;
- */
- /* TODO
- TVertices pS;
- 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 );
- */
-
- /*
- } // rof
- */
- }
-
- 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 )
- {
- 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 );
- }
-
-private:
- TDijkstra( const Self& other );
- Self& operator=( const Self& other );
-
-protected:
- _TCandidates m_Candidates;
-public:
- vtkSmartPointer< vtkPolyData > m_Axes;
-};
+typedef fpa::VTK::Image3DObserver< TDijkstra, vtkRenderWindow > TDijkstraObs;
// -------------------------------------------------------------------------
int main( int argc, char* argv[] )
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( );
+ // Create polydata
+ vtkSmartPointer< vtkPoints > points =
+ vtkSmartPointer< vtkPoints >::New( );
+ vtkSmartPointer< vtkCellArray > cells =
+ vtkSmartPointer< vtkCellArray >::New( );
+ vtkSmartPointer< vtkFloatArray > scalars =
+ vtkSmartPointer< vtkFloatArray >::New( );
+
+ const TDijkstra::TVertices& endpoints = paths->GetEndPoints( );
+ const TDijkstra::TTree& tree = paths->GetFinalTree( );
+ TDijkstra::TVertices::const_iterator epIt = endpoints.begin( );
+ for( unsigned int epId = 0; epIt != endpoints.end( ); ++epIt, ++epId )
+ {
+ double pd = double( epId ) / double( endpoints.size( ) - 1 );
- TImageWriter::Pointer dijkstra_writer =
- TImageWriter::New( );
- dijkstra_writer->SetInput( paths->GetOutput( ) );
- dijkstra_writer->SetFileName( "dijkstra.mhd" );
- dijkstra_writer->Update( );
+ TDijkstra::TVertex idx = *epIt;
+ do
+ {
+ TImage::PointType pnt;
+ input_image->TransformIndexToPhysicalPoint( idx, pnt );
- // Show result
- /*
- TVTKImage::Pointer output_image_vtk = TVTKImage::New( );
- output_image_vtk->SetInput( segmentor->GetOutput( ) );
- output_image_vtk->Update( );
+ points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
+ scalars->InsertNextTuple1( pd );
+ if( idx != *epIt )
+ {
+ cells->InsertNextCell( 2 );
+ cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
+ cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
- vtkSmartPointer< vtkImageMarchingCubes > mc =
- vtkSmartPointer< vtkImageMarchingCubes >::New( );
- mc->SetInputData( output_image_vtk->GetOutput( ) );
- mc->SetValue(
- 0,
- double( segmentor->GetInsideValue( ) ) * double( 0.95 )
- );
- mc->Update( );
+ } // fi
+ idx = tree.find( idx )->second;
- // Let some interaction and close program
- view.AddPolyData( mc->GetOutput( ), 0.1, 0.6, 0.8, 0.5 );
- view.Start( );
+ } while( idx != tree.find( idx )->second );
- // 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 );
+ } // rof
- } // yrt
+ vtkSmartPointer< vtkPolyData > vtk_tree =
+ vtkSmartPointer< vtkPolyData >::New( );
+ vtk_tree->SetPoints( points );
+ vtk_tree->SetLines( cells );
+ vtk_tree->GetPointData( )->SetScalars( scalars );
+
+ view.AddPolyData( vtk_tree );
+ view.Render( );
+ view.Start( );
+
+ /* TODO
+ TImageWriter::Pointer dijkstra_writer =
+ TImageWriter::New( );
+ dijkstra_writer->SetInput( paths->GetOutput( ) );
+ dijkstra_writer->SetFileName( "dijkstra.mhd" );
+ dijkstra_writer->Update( );
*/
+
return( 0 );
}
#include <itkImageFileWriter.h>
#include <itkImageToVTKImageFilter.h>
+#include <vtkImageMarchingCubes.h>
+
#include <fpa/Image/RegionGrowWithMultipleThresholds.h>
+#include <fpa/Image/DijkstraWithSphereBacktracking.h>
#include <fpa/VTK/ImageMPR.h>
#include <fpa/VTK/Image3DObserver.h>
typedef itk::Image< TPixel, Dim > TImage;
typedef itk::Image< TScalar, Dim > TDistanceMap;
typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
+typedef itk::ImageToVTKImageFilter< TDistanceMap > TVTKDistanceMap;
typedef itk::ImageFileReader< TImage > TImageReader;
typedef itk::ImageFileWriter< TImage > TImageWriter;
typedef
itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap >
TDistanceFilter;
+typedef
+fpa::Image::DijkstraWithSphereBacktracking< TDistanceMap, TScalar >
+TDijkstra;
typedef fpa::VTK::ImageMPR TMPR;
typedef fpa::VTK::Image3DObserver< TSegmentor, vtkRenderWindow > TSegmentorObs;
+typedef fpa::VTK::Image3DObserver< TDijkstra, vtkRenderWindow > TDijkstraObs;
// -------------------------------------------------------------------------
int main( int argc, char* argv[] )
seconds = double( end - start ) / double( CLOCKS_PER_SEC );
std::cout << "Distance map time = " << seconds << std::endl;
- std::cout << "Final seed: " << seed_idx << std::endl;
+ TVTKDistanceMap::Pointer vtk_distanceMap = TVTKDistanceMap::New( );
+ vtk_distanceMap->SetInput( distanceMap->GetOutput( ) );
+ vtk_distanceMap->Update( );
+
+ vtkSmartPointer< vtkImageMarchingCubes > mc =
+ vtkSmartPointer< vtkImageMarchingCubes >::New( );
+ mc->SetInputData( vtk_distanceMap->GetOutput( ) );
+ mc->SetValue( 0, 1e-2 );
+ mc->Update( );
+ view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 );
+ view.Render( );
+
+ // Extract paths
+ TDijkstra::Pointer paths = TDijkstra::New( );
+ paths->AddSeed( segmentor->GetSeed( 0 ), TScalar( 0 ) );
+ paths->SetInput( distanceMap->GetOutput( ) );
+ paths->SetNeighborhoodOrder( 1 );
+
+ if( visual_debug )
+ {
+ // Configure observer
+ TDijkstraObs::Pointer obs = TDijkstraObs::New( );
+ obs->SetRenderWindow( 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;
+
+ std::cout << "Final seed: " << paths->GetSeed( 0 ) << std::endl;
+
+ // Create polydata
+ vtkSmartPointer< vtkPoints > points =
+ vtkSmartPointer< vtkPoints >::New( );
+ vtkSmartPointer< vtkCellArray > cells =
+ vtkSmartPointer< vtkCellArray >::New( );
+ vtkSmartPointer< vtkFloatArray > scalars =
+ vtkSmartPointer< vtkFloatArray >::New( );
+
+ const TDijkstra::TVertices& endpoints = paths->GetEndPoints( );
+ const TDijkstra::TTree& tree = paths->GetFinalTree( );
+ TDijkstra::TVertices::const_iterator epIt = endpoints.begin( );
+ for( unsigned int epId = 0; epIt != endpoints.end( ); ++epIt, ++epId )
+ {
+ double pd = double( epId ) / double( endpoints.size( ) - 1 );
+
+ TDijkstra::TVertex idx = *epIt;
+ do
+ {
+ TImage::PointType pnt;
+ input_image->TransformIndexToPhysicalPoint( idx, pnt );
+
+ points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
+ scalars->InsertNextTuple1( pd );
+ if( idx != *epIt )
+ {
+ cells->InsertNextCell( 2 );
+ cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
+ cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
+
+ } // fi
+ idx = tree.find( idx )->second;
+
+ } while( idx != tree.find( idx )->second );
+
+ } // rof
- TDistanceMapWriter::Pointer distancemap_writer =
- TDistanceMapWriter::New( );
- distancemap_writer->SetInput( distanceMap->GetOutput( ) );
- distancemap_writer->SetFileName( output_distancemap_fn );
- distancemap_writer->Update( );
+ vtkSmartPointer< vtkPolyData > vtk_tree =
+ vtkSmartPointer< vtkPolyData >::New( );
+ vtk_tree->SetPoints( points );
+ vtk_tree->SetLines( cells );
+ vtk_tree->GetPointData( )->SetScalars( scalars );
- TImageWriter::Pointer segmentation_writer =
- TImageWriter::New( );
- segmentation_writer->SetInput( segmentor->GetOutput( ) );
- segmentation_writer->SetFileName( output_segmentation_fn );
- segmentation_writer->Update( );
+ view.AddPolyData( vtk_tree );
+ view.Render( );
+ view.Start( );
+
+ /* TODO
+ TDistanceMapWriter::Pointer distancemap_writer =
+ TDistanceMapWriter::New( );
+ distancemap_writer->SetInput( distanceMap->GetOutput( ) );
+ distancemap_writer->SetFileName( output_distancemap_fn );
+ distancemap_writer->Update( );
+
+ TImageWriter::Pointer segmentation_writer =
+ TImageWriter::New( );
+ segmentation_writer->SetInput( segmentor->GetOutput( ) );
+ segmentation_writer->SetFileName( output_segmentation_fn );
+ segmentation_writer->Update( );
+ */
// Show result
/*
typedef std::vector< _TCollisionSitesRow > _TCollisionSites;
public:
- typedef BaseEvent< _TNode > TEvent;
- typedef FrontEvent< _TNode > TFrontEvent;
- typedef MarkEvent< _TNode > TMarkEvent;
- typedef CollisionEvent< _TNode > TCollisionEvent;
- typedef EndEvent< _TNode > TEndEvent;
- typedef BacktrackingEvent< TVertex > TBacktrackingEvent;
+ typedef BaseEvent< _TNode > TEvent;
+ typedef FrontEvent< _TNode > TFrontEvent;
+ typedef MarkEvent< _TNode > TMarkEvent;
+ typedef CollisionEvent< _TNode > TCollisionEvent;
+ typedef EndEvent< _TNode > TEndEvent;
+ typedef BacktrackingEvent< TVertex > TBacktrackingEvent;
+ typedef EndBacktrackingEvent< TVertex > TEndBacktrackingEvent;
public:
itkTypeMacro( Algorithm, itkProcessObject );
{ }
BacktrackingEvent( const N& n, const unsigned long& id )
: BaseEvent< N >( n ),
- BackId( id )
+ BackId( id )
{ }
virtual ~BacktrackingEvent( )
{ }
unsigned long BackId;
};
+ /**
+ */
+ template< class N >
+ class EndBacktrackingEvent
+ : public BaseEvent< N >
+ {
+ public:
+ EndBacktrackingEvent( )
+ : BaseEvent< N >( )
+ { }
+ EndBacktrackingEvent( const unsigned long& id )
+ : BaseEvent< N >( ),
+ BackId( id )
+ { }
+ virtual ~EndBacktrackingEvent( )
+ { }
+ const char* GetEventName( ) const
+ { return( "fpa::Base::EndBacktrackingEvent" ); }
+ bool CheckEvent( const itk::EventObject* e ) const
+ {
+ return(
+ dynamic_cast< const EndBacktrackingEvent< N >* >( e ) != NULL
+ );
+ }
+ itk::EventObject* MakeObject( ) const
+ { return( new EndBacktrackingEvent< N >( ) ); }
+
+ unsigned long BackId;
+ };
+
} // ecapseman
} // ecapseman
--- /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 typename Superclass::TTraits::TVertexCmp TVertexCmp;
+ typedef std::map< TVertex, TVertex, 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( FinalTree, TTree );
+ itkGetConstMacro( EndPoints, TVertices );
+
+ 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_FinalTree;
+ TVertices m_EndPoints;
+ };
+
+ } // 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 >
+fpa::Image::DijkstraWithSphereBacktracking< I, C >::
+DijkstraWithSphereBacktracking( )
+ : Superclass( )
+{
+}
+
+// -------------------------------------------------------------------------
+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( 3 );
+
+ // 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( )
+{
+ typedef itk::Image< bool, I::ImageDimension > _TMarkImage;
+
+ // Finish base algorithm
+ this->Superclass::_AfterMainLoop( );
+ this->m_FinalTree.clear( );
+ this->m_EndPoints.clear( );
+ if( this->m_Candidates.size( ) == 0 )
+ return;
+ this->InvokeEvent( TEndEvent( ) );
+
+ // Get 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 = _TMarkImage::New( );
+ marks->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
+ marks->SetRequestedRegion( input->GetRequestedRegion( ) );
+ marks->SetBufferedRegion( input->GetBufferedRegion( ) );
+ marks->SetOrigin( input->GetOrigin( ) );
+ marks->SetSpacing( input->GetSpacing( ) );
+ marks->SetDirection( input->GetDirection( ) );
+ marks->Allocate( );
+ marks->FillBuffer( false );
+
+ // Iterate over the candidates, starting from the candidates that
+ // are near thin branches
+ typename _TCandidates::const_reverse_iterator cIt =
+ this->m_Candidates.rbegin( );
+ for( int backId = 0; cIt != this->m_Candidates.rend( ); ++cIt )
+ {
+ // If pixel hasn't been visited, continue
+ TVertex v = cIt->second;
+ if( marks->GetPixel( v ) )
+ 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 ) )
+ continue;
+ backId++;
+ this->m_EndPoints.push_back( max_vertex );
+
+ // Construct path (at least the part that hasn't been iterated)
+ do
+ {
+ if( this->m_FinalTree.find( max_vertex ) == this->m_FinalTree.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.25 ) );
+ itk::ImageRegionIteratorWithIndex< _TMarkImage >
+ mIt( marks, region );
+ for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
+ mIt.Set( true );
+
+ // Next vertex in current path
+ this->InvokeEvent( TBacktrackingEvent( max_vertex, backId ) );
+ this->m_FinalTree[ max_vertex ] = this->_Parent( max_vertex );
+
+ } // fi
+ max_vertex = this->_Parent( max_vertex );
+
+ } while( max_vertex != this->_Parent( max_vertex ) );
+ this->m_FinalTree[ max_vertex ] = max_vertex;
+ this->InvokeEvent( TEndBacktrackingEvent( backId ) );
+
+ } // 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 )
+ {
+ nn.Cost = n.Cost + ( TCost( 1 ) / nc );
+ 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$
void fpa::VTK::Image3DObserver< F, R >::
Execute( const itk::Object* c, const itk::EventObject& e )
{
- typedef itk::ImageBase< 3 > _TImage;
- typedef typename F::TEvent _TEvent;
- typedef typename F::TFrontEvent _TFrontEvent;
- typedef typename F::TMarkEvent _TMarkEvent;
- typedef typename F::TCollisionEvent _TCollisionEvent;
- typedef typename F::TEndEvent _TEndEvent;
- typedef typename F::TBacktrackingEvent _TBacktrackingEvent;
+ typedef itk::ImageBase< 3 > _TImage;
+ typedef typename F::TEvent _TEvent;
+ typedef typename F::TFrontEvent _TFrontEvent;
+ typedef typename F::TMarkEvent _TMarkEvent;
+ typedef typename F::TCollisionEvent _TCollisionEvent;
+ typedef typename F::TEndEvent _TEndEvent;
+ typedef typename F::TBacktrackingEvent _TBacktrackingEvent;
+ typedef typename F::TEndBacktrackingEvent _TEndBacktrackingEvent;
// Check inputs
if( this->m_RenderWindow == NULL )
{
const _TBacktrackingEvent* bevt =
dynamic_cast< const _TBacktrackingEvent* >( &e );
+ const _TEndBacktrackingEvent* ebevt =
+ dynamic_cast< const _TEndBacktrackingEvent* >( &e );
if( bevt != NULL )
{
static const unsigned long nColors = 10;
this->m_Data->GetPointData( )->
GetScalars( )->InsertNextTuple1( back_id );
this->m_Data->Modified( );
- this->m_RenderWindow->Render( );
return;
} // fi
+ if( ebevt != NULL )
+ {
+ this->m_RenderWindow->Render( );
+ return;
+
+ } // fi
+
} // fi
}
} // fi
}
break;
+ case 'z':
+ case 'Z':
+ {
+ this->SeedWidget->ProcessEventsOff( );
+ this->WidgetX->InteractionOff( );
+ this->WidgetY->InteractionOff( );
+ this->WidgetZ->InteractionOff( );
+ }
+ break;
default:
break;
this->m_Window->SetSize( w, h );
}
+// -------------------------------------------------------------------------
+void fpa::VTK::ImageMPR::
+AddPolyData( vtkPolyData* pd, double opacity )
+{
+ unsigned int i = this->m_PolyDatas.size( );
+
+ this->m_PolyDatas.push_back( pd );
+ this->m_Mappers.push_back( vtkSmartPointer< vtkPolyDataMapper >::New( ) );
+ this->m_Actors.push_back( vtkSmartPointer< vtkActor >::New( ) );
+
+ this->m_Mappers[ i ]->SetInputData( pd );
+ this->m_Actors[ i ]->SetMapper( this->m_Mappers[ i ] );
+ this->m_Actors[ i ]->GetProperty( )->SetOpacity( opacity );
+ this->m_Renderer->AddActor( this->m_Actors[ i ] );
+}
+
// -------------------------------------------------------------------------
void fpa::VTK::ImageMPR::
AddPolyData( vtkPolyData* pd, double r, double g, double b, double opacity )
this->m_Interactor->Start( );
}
+// -------------------------------------------------------------------------
+void fpa::VTK::ImageMPR::
+Render( )
+{
+ this->m_Window->Render( );
+}
+
// eof - $RCSfile$
void SetBackground( double r, double g, double b );
void SetSize( unsigned int w, unsigned int h );
+ void AddPolyData( vtkPolyData* pd, double opacity = double( 1 ) );
void AddPolyData(
vtkPolyData* pd,
double r, double g, double b, double opacity = double( 1 )
vtkRenderer* GetRenderer( ) const;
void Start( );
+ void Render( );
protected:
vtkSmartPointer< vtkImageData > m_Image;
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / commitdiff