]> Creatis software - FrontAlgorithms.git/commitdiff
Refactoring: gaussian model estimator
authorLeonardo Florez-Valencia <florez-l@javeriana.edu.co>
Mon, 6 Apr 2015 23:31:51 +0000 (18:31 -0500)
committerLeonardo Florez-Valencia <florez-l@javeriana.edu.co>
Mon, 6 Apr 2015 23:31:51 +0000 (18:31 -0500)
12 files changed:
appli/examples/CMakeLists.txt
appli/examples/example_Image_Dijkstra_DanielssonCost_TwoSeedsPath.cxx [new file with mode: 0644]
appli/examples/example_Image_RegionGrow_GaussianModelEstimation.cxx [moved from appli/examples/example_ImageAlgorithmRegionGrow_GaussianModelEstimator.cxx with 50% similarity]
lib/fpa/Base/Algorithm.hxx
lib/fpa/Base/MinimumSpanningTree.h [new file with mode: 0644]
lib/fpa/Base/MinimumSpanningTree.hxx [new file with mode: 0644]
lib/fpa/Base/TreeExtractor.h [deleted file]
lib/fpa/Base/TreeExtractor.hxx [deleted file]
lib/fpa/Image/Algorithm.h
lib/fpa/Image/Algorithm.hxx
lib/fpa/Image/MinimumSpanningTree.h [new file with mode: 0644]
lib/fpa/Image/MinimumSpanningTree.hxx [new file with mode: 0644]

index 985071e16c2f9032570a91adb05acd6070ad9227..5efc07716857475f7c6836e27a581e6145bd2ba8 100644 (file)
@@ -3,9 +3,11 @@ IF(USE_VTK)
     SIMPLE_VTK_EXAMPLES
     example_Image_RegionGrow_AllPixels
     example_Image_RegionGrow_AllRGBPixels
+    example_Image_RegionGrow_GaussianModelEstimation
     example_Image_Dijkstra_CostFromInput
     example_Image_Dijkstra_CostFromRGBInput
     example_Image_Dijkstra_DanielssonCost
+    example_Image_Dijkstra_DanielssonCost_TwoSeedsPath
     )
   FOREACH(EX ${SIMPLE_VTK_EXAMPLES})
     ADD_EXECUTABLE(${EX} ${EX}.cxx)
diff --git a/appli/examples/example_Image_Dijkstra_DanielssonCost_TwoSeedsPath.cxx b/appli/examples/example_Image_Dijkstra_DanielssonCost_TwoSeedsPath.cxx
new file mode 100644 (file)
index 0000000..8a5eab0
--- /dev/null
@@ -0,0 +1,253 @@
+#include <cstdlib>
+#include <iostream>
+#include <limits>
+#include <string>
+
+#include <itkImage.h>
+#include <itkRGBAPixel.h>
+#include <itkImageToVTKImageFilter.h>
+
+#include <itkImageFileReader.h>
+#include <itkMinimumMaximumImageCalculator.h>
+#include <itkInvertIntensityImageFilter.h>
+#include <itkDanielssonDistanceMapImageFilter.h>
+#include <itkImageFileWriter.h>
+
+#include <vtkCamera.h>
+#include <vtkImageActor.h>
+#include <vtkInteractorStyleImage.h>
+#include <vtkPointHandleRepresentation3D.h>
+#include <vtkProperty.h>
+#include <vtkRenderer.h>
+#include <vtkRenderWindow.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkSeedRepresentation.h>
+#include <vtkSeedWidget.h>
+#include <vtkSmartPointer.h>
+
+#include <fpa/Image/Dijkstra.h>
+#include <fpa/Base/Functors/InvertCostFunction.h>
+#include <fpa/VTK/Image2DObserver.h>
+
+// -------------------------------------------------------------------------
+const unsigned int Dim = 2;
+typedef unsigned char TInputPixel;
+typedef float         TOutputPixel;
+
+typedef itk::Image< TInputPixel, Dim >            TInputImage;
+typedef itk::Image< TOutputPixel, Dim >           TOutputImage;
+typedef itk::ImageToVTKImageFilter< TInputImage > TVTKInputImage;
+
+// -------------------------------------------------------------------------
+int main( int argc, char* argv[] )
+{
+  if( argc < 5 )
+  {
+    std::cerr
+      << "Usage: " << argv[ 0 ]
+      << " input_image output_image neighborhood_order stop_at_one_front"
+      << std::endl;
+    return( 1 );
+
+  } // fi
+  std::string input_image_fn = argv[ 1 ];
+  std::string output_image_fn = argv[ 2 ];
+  unsigned int neighborhood_order = std::atoi( argv[ 3 ] );
+  bool stop_at_one_front = ( std::atoi( argv[ 4 ] ) != 0 );
+
+  // Read image
+  itk::ImageFileReader< TInputImage >::Pointer input_image_reader =
+    itk::ImageFileReader< TInputImage >::New( );
+  input_image_reader->SetFileName( input_image_fn );
+  try
+  {
+    input_image_reader->Update( );
+  }
+  catch( itk::ExceptionObject& err )
+  {
+    std::cerr
+      << "Error while reading image from " << input_image_fn << ": "
+      << err << std::endl;
+    return( 1 );
+
+  } // yrt
+  TInputImage::ConstPointer input_image = input_image_reader->GetOutput( );
+  TVTKInputImage::Pointer vtk_input_image = TVTKInputImage::New( );
+  vtk_input_image->SetInput( input_image );
+  vtk_input_image->Update( );
+
+  // VTK visualization
+  vtkSmartPointer< vtkImageActor > actor =
+    vtkSmartPointer< vtkImageActor >::New( );
+  actor->SetInputData( vtk_input_image->GetOutput( ) );
+
+  vtkSmartPointer< vtkRenderer > renderer =
+    vtkSmartPointer< vtkRenderer >::New( );
+  renderer->SetBackground( 0.1, 0.2, 0.7 );
+  renderer->AddActor( actor );
+  vtkSmartPointer< vtkRenderWindow > window =
+    vtkSmartPointer< vtkRenderWindow >::New( );
+  window->SetSize( 800, 800 );
+  window->AddRenderer( renderer );
+
+  // Correct camera due to the loaded image
+  vtkCamera* camera = renderer->GetActiveCamera( );
+  camera->SetViewUp( 0, -1, 0 );
+  camera->SetPosition( 0, 0, -1 );
+  camera->SetFocalPoint( 0, 0, 0 );
+
+  // VTK interaction
+  vtkSmartPointer< vtkInteractorStyleImage > imageStyle =
+    vtkSmartPointer< vtkInteractorStyleImage >::New( );
+  vtkSmartPointer< vtkRenderWindowInteractor > interactor =
+    vtkSmartPointer< vtkRenderWindowInteractor >::New( );
+  interactor->SetInteractorStyle( imageStyle );
+  window->SetInteractor( interactor );
+
+  // Create the widget and its representation
+  vtkSmartPointer< vtkPointHandleRepresentation3D > handle =
+    vtkSmartPointer< vtkPointHandleRepresentation3D >::New( );
+  handle->GetProperty( )->SetColor( 1, 0, 0 );
+  vtkSmartPointer< vtkSeedRepresentation > rep =
+    vtkSmartPointer< vtkSeedRepresentation >::New( );
+  rep->SetHandleRepresentation( handle );
+
+  vtkSmartPointer< vtkSeedWidget > widget =
+    vtkSmartPointer< vtkSeedWidget >::New( );
+  widget->SetInteractor( interactor );
+  widget->SetRepresentation( rep );
+
+  // Let some interaction
+  interactor->Initialize( );
+  renderer->ResetCamera( );
+  window->Render( );
+  widget->On( );
+  interactor->Start( );
+
+  // Invert input image
+  itk::MinimumMaximumImageCalculator< TInputImage >::Pointer minmax =
+    itk::MinimumMaximumImageCalculator< TInputImage >::New( );
+  minmax->SetImage( input_image );
+  minmax->Compute( );
+
+  itk::InvertIntensityImageFilter< TInputImage >::Pointer invert =
+    itk::InvertIntensityImageFilter< TInputImage >::New( );
+  invert->SetInput( input_image );
+  invert->SetMaximum( minmax->GetMaximum( ) );
+
+  itk::DanielssonDistanceMapImageFilter< TInputImage, TOutputImage >::Pointer dmap =
+    itk::DanielssonDistanceMapImageFilter< TInputImage, TOutputImage >::New( );
+  dmap->SetInput( invert->GetOutput( ) );
+  dmap->InputIsBinaryOn( );
+  dmap->SquaredDistanceOn( );
+  dmap->UseImageSpacingOn( );
+  dmap->Update( );
+
+  typedef fpa::Base::Functors::InvertCostFunction< TOutputPixel > TFunction;
+  TFunction::Pointer function = TFunction::New( );
+
+  // Prepare region grow filter
+  typedef fpa::Image::Dijkstra< TOutputImage, TOutputImage > TFilter;
+  TFilter::Pointer filter = TFilter::New( );
+  filter->SetInput( dmap->GetOutput( ) );
+  filter->SetConversionFunction( function );
+  filter->SetNeighborhoodOrder( neighborhood_order );
+  filter->SetStopAtOneFront( stop_at_one_front );
+
+  // Get user-given seeds
+  for( unsigned int s = 0; s < rep->GetNumberOfSeeds( ); s++ )
+  {
+    double pos[ 3 ];
+    rep->GetSeedWorldPosition( s, pos );
+
+    TInputImage::PointType pnt;
+    pnt[ 0 ] = TInputImage::PointType::ValueType( pos[ 0 ] );
+    pnt[ 1 ] = TInputImage::PointType::ValueType( pos[ 1 ] );
+
+    TInputImage::IndexType idx;
+    if( input_image->TransformPhysicalPointToIndex( pnt, idx ) )
+      filter->AddSeed( idx, 0 );
+
+  } // rof
+
+  // Prepare graphical debugger
+  typedef fpa::VTK::Image2DObserver< TFilter, vtkRenderWindow > TDebugger;
+  TDebugger::Pointer debugger = TDebugger::New( );
+  debugger->SetRenderWindow( window );
+  debugger->SetRenderPercentage( 0.001 );
+  filter->AddObserver( itk::AnyEvent( ), debugger );
+  filter->ThrowEventsOn( );
+
+  // Go!
+  filter->Update( );
+
+  // Save final total cost map
+  itk::ImageFileWriter< TOutputImage >::Pointer output_image_writer =
+    itk::ImageFileWriter< TOutputImage >::New( );
+  output_image_writer->SetFileName( output_image_fn );
+  output_image_writer->SetInput( filter->GetOutput( ) );
+  try
+  {
+    output_image_writer->Update( );
+  }
+  catch( itk::ExceptionObject& err )
+  {
+    std::cerr
+      << "Error while writing image to " << output_image_fn << ": "
+      << err << std::endl;
+    return( 1 );
+
+  } // yrt
+
+  // Draw the path between the two seeds
+  if( filter->GetNumberOfSeeds( ) > 1 )
+  {
+    TInputImage::IndexType a = filter->GetSeed( 0 );
+    TInputImage::IndexType b =
+      filter->GetSeed( filter->GetNumberOfSeeds( ) - 1 );
+
+    std::vector< TInputImage::IndexType > path;
+    filter->GetMinimumSpanningTree( )->GetPath( path, a, b );
+
+    typedef itk::Image< itk::RGBAPixel< unsigned char >, Dim > TColorImage;
+    TColorImage::Pointer path_image = TColorImage::New( );
+    path_image->SetLargestPossibleRegion(
+      input_image->GetLargestPossibleRegion( )
+      );
+    path_image->SetRequestedRegion( input_image->GetRequestedRegion( ) );
+    path_image->SetBufferedRegion( input_image->GetBufferedRegion( ) );
+    path_image->SetSpacing( input_image->GetSpacing( ) );
+    path_image->SetOrigin( input_image->GetOrigin( ) );
+    path_image->SetDirection( input_image->GetDirection( ) );
+    path_image->Allocate( );
+
+    unsigned char background[ 4 ] = { 0, 0, 0, 0 };
+    unsigned char color[ 4 ] = { 255, 0, 0, 255 };
+    path_image->FillBuffer( itk::RGBAPixel< unsigned char >( background ) );
+
+    for( unsigned int pId = 0; pId < path.size( ); ++pId )
+      path_image->SetPixel(
+        path[ pId ], itk::RGBAPixel< unsigned char >( color )
+        );
+
+    itk::ImageToVTKImageFilter< TColorImage >::Pointer vtk_path_image =
+      itk::ImageToVTKImageFilter< TColorImage >::New( );
+    vtk_path_image->SetInput( path_image );
+    vtk_path_image->Update( );
+
+    vtkSmartPointer< vtkImageActor > path_image_actor =
+      vtkSmartPointer< vtkImageActor >::New( );
+    path_image_actor->SetInputData( vtk_path_image->GetOutput( ) );
+    renderer->AddActor( path_image_actor );
+
+    // More interaction
+    window->Render( );
+    widget->Off( );
+    interactor->Start( );
+
+  } // fi
+
+  return( 0 );
+}
+
+// eof - $RCSfile$
similarity index 50%
rename from appli/examples/example_ImageAlgorithmRegionGrow_GaussianModelEstimator.cxx
rename to appli/examples/example_Image_RegionGrow_GaussianModelEstimation.cxx
index 04fa936b43f9a1a0ca6bd8b95ac6836402ee9094..919409aa948c2838cc253fb39eb7d2f20c392383 100644 (file)
@@ -1,12 +1,15 @@
+#include <cstdlib>
 #include <iostream>
 #include <limits>
 #include <string>
 
 #include <itkImage.h>
-#include <itkImageFileReader.h>
-#include <itkSignedDanielssonDistanceMapImageFilter.h>
 #include <itkImageToVTKImageFilter.h>
+#include <itkRGBPixel.h>
+
+#include <itkImageFileReader.h>
 
+#include <vtkCamera.h>
 #include <vtkImageActor.h>
 #include <vtkInteractorStyleImage.h>
 #include <vtkPointHandleRepresentation3D.h>
 #include <vtkSeedWidget.h>
 #include <vtkSmartPointer.h>
 
-#include <fpa/Image/RegionGrow.h>
-#include <fpa/Image/Functors/RegionGrowAllBelongsFunction.h>
-#include <fpa/VTK/Image2DObserver.h>
 #include <cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.h>
 #include <cpPlugins/Extensions/Algorithms/RGBToYPbPrFunction.h>
 
+#include <fpa/Image/RegionGrow.h>
+#include <fpa/Base/Functors/TautologyFunction.h>
+#include <fpa/VTK/Image2DObserver.h>
+
 // -------------------------------------------------------------------------
 const unsigned int Dim = 2;
-typedef unsigned char TChannel;
-typedef itk::RGBPixel< TChannel > TPixel;
-typedef double TScalar;
-typedef itk::Image< TChannel, Dim > TImage;
-typedef itk::Image< TPixel, Dim >   TColorImage;
-typedef itk::ImageToVTKImageFilter< TColorImage > TVTKImage;
-
-typedef itk::ImageFileReader< TColorImage >   TColorImageReader;
-typedef fpa::Image::RegionGrow< TColorImage, TImage, fpa::Image::Functors::CastVertexValueToConstantCost< TPixel, bool > > TFrontAlgorithm;
+typedef unsigned char TPixel;
+typedef float         TScalar;
 
-typedef
-fpa::VTK::Image2DObserver< TFrontAlgorithm, vtkRenderWindow >
-TObserver;
+typedef itk::Image< itk::RGBPixel< TPixel >, Dim > TColorImage;
+typedef itk::Image< TPixel, Dim >                  TImage;
+typedef itk::ImageToVTKImageFilter< TColorImage >  TVTKImage;
 
 // -------------------------------------------------------------------------
 template< class I, class S >
 class GaussianFunction
-  : public fpa::Image::Functors::RegionGrowAllBelongsFunction< I >
+  : public itk::FunctionBase< typename I::PixelType, bool >
 {
 public:
   // Type-related and pointers
-  typedef GaussianFunction                                        Self;
-  typedef fpa::Image::Functors::RegionGrowAllBelongsFunction< I > Superclass;
-  typedef itk::SmartPointer< Self >                               Pointer;
-  typedef itk::SmartPointer< const Self >                         ConstPointer;
+  typedef GaussianFunction                                 Self;
+  typedef itk::FunctionBase< typename I::PixelType, bool > Superclass;
+  typedef itk::SmartPointer< Self >                        Pointer;
+  typedef itk::SmartPointer< const Self >                  ConstPointer;
 
   // Superclass' types
-  typedef typename Superclass::TInputImage  TInputImage;
-  typedef typename Superclass::TOutputValue TOutputValue;
-  typedef typename Superclass::TIndex       TIndex;
-
   typedef cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, 3 > TEstimator;
   typedef cpPlugins::Extensions::Algorithms::RGBToYPbPrFunction< S > TYPbPrFunction;
 
 public:
   itkNewMacro( Self );
-  itkTypeMacro(
-    GaussianFunction,
-    fpa_Image_Functors_RegionGrowAllBelongsFunction
-    );
+  itkTypeMacro( GaussianFunction, itkFunctionBase );
 
   itkGetConstMacro( ModelSupport, unsigned long );
   itkSetMacro( ModelSupport, unsigned long );
 
-  typedef itk::Image< bool, I::ImageDimension > TMarks;
-
 public:
-  virtual void SetInputImage( TInputImage* img )
+  virtual bool Evaluate( const typename I::PixelType& rgb ) const
     {
-      std::cout << "Set!!!" << std::endl;
-      this->Superclass::SetInputImage( img );
-      this->m_Marks = TMarks::New( );
-      this->m_Marks->SetLargestPossibleRegion( img->GetLargestPossibleRegion( ) );
-      this->m_Marks->SetRequestedRegion( img->GetRequestedRegion( ) );
-      this->m_Marks->SetBufferedRegion( img->GetBufferedRegion( ) );
-      this->m_Marks->SetOrigin( img->GetOrigin( ) );
-      this->m_Marks->SetSpacing( img->GetSpacing( ) );
-      this->m_Marks->SetDirection( img->GetDirection( ) );
-      this->m_Marks->Allocate( );
-      this->m_Marks->FillBuffer( false );
-    }
-
-
-  virtual TOutputValue Evaluate( const TIndex& idx ) const
-    {
-      const TInputImage* img = this->GetInputImage( );
-      typename TInputImage::PixelType rgb = img->GetPixel( idx );
-
+      /* TODO
       if( !this->m_Estimating )
       {
-        if( !( this->m_Marks->GetPixel( idx ) ) )
+        this->m_Estimator->AddSample( this->m_YPbPrFunction( rgb ) );
+        if( this->m_Estimator->GetNumberOfSamples( ) == this->m_ModelSupport )
         {
-          this->m_Estimator->AddSample( this->m_YPbPrFunction( rgb ) );
-          this->m_Marks->SetPixel( idx, true );
-          if( this->m_Estimator->GetNumberOfSamples( ) == this->m_ModelSupport )
-          {
-            this->m_Estimating = true;
-            this->m_Estimator->UpdateModel( );
-            std::cout << this->m_Estimator->GetMinimumProbability( ) << std::endl;
-            std::cout << this->m_Estimator->GetMaximumProbability( ) << std::endl;
-
-          } // fi
+          this->m_Estimating = true;
+          this->m_Estimator->UpdateModel( );
+          std::cout << this->m_Estimator->GetMinimumProbability( ) << std::endl;
+          std::cout << this->m_Estimator->GetMaximumProbability( ) << std::endl;
 
         } // fi
         return( true );
@@ -118,6 +83,8 @@ public:
         return( p > this->m_Estimator->GetMinimumProbability( ) );
 
       } // fi
+      */
+      return( true );
     }
 
 protected:
@@ -143,27 +110,27 @@ protected:
 
   unsigned long m_ModelSupport;
   mutable bool m_Estimating;
-  typename TMarks::Pointer m_Marks;
 };
 
-typedef GaussianFunction< TColorImage, TScalar > TMembershipFunction;
-
 // -------------------------------------------------------------------------
 int main( int argc, char* argv[] )
 {
-  if( argc < 3 )
+  if( argc < 4 )
   {
     std::cerr
       << "Usage: " << argv[ 0 ]
-      << " input_image support" << std::endl;
+      << " input_image neighborhood_order support"
+      << std::endl;
     return( 1 );
 
   } // fi
   std::string input_image_fn = argv[ 1 ];
-  unsigned long support = std::atoi( argv[ 2 ] );
+  unsigned int neighborhood_order = std::atoi( argv[ 2 ] );
+  unsigned long support = std::atoi( argv[ 3 ] );
 
   // Read image
-  TColorImageReader::Pointer input_image_reader = TColorImageReader::New( );
+  itk::ImageFileReader< TColorImage >::Pointer input_image_reader =
+    itk::ImageFileReader< TColorImage >::New( );
   input_image_reader->SetFileName( input_image_fn );
   try
   {
@@ -171,20 +138,21 @@ int main( int argc, char* argv[] )
   }
   catch( itk::ExceptionObject& err )
   {
-    std::cerr << "Error caught: " << err << std::endl;
+    std::cerr
+      << "Error while reading image from " << input_image_fn << ": "
+      << err << std::endl;
     return( 1 );
 
   } // yrt
   TColorImage::ConstPointer input_image = input_image_reader->GetOutput( );
-
-  TVTKImage::Pointer vtk_image = TVTKImage::New( );
-  vtk_image->SetInput( input_image );
-  vtk_image->Update( );
+  TVTKImage::Pointer vtk_input_image = TVTKImage::New( );
+  vtk_input_image->SetInput( input_image );
+  vtk_input_image->Update( );
 
   // VTK visualization
   vtkSmartPointer< vtkImageActor > actor =
     vtkSmartPointer< vtkImageActor >::New( );
-  actor->SetInputData( vtk_image->GetOutput( ) );
+  actor->SetInputData( vtk_input_image->GetOutput( ) );
 
   vtkSmartPointer< vtkRenderer > renderer =
     vtkSmartPointer< vtkRenderer >::New( );
@@ -195,6 +163,12 @@ int main( int argc, char* argv[] )
   window->SetSize( 800, 800 );
   window->AddRenderer( renderer );
 
+  // Correct camera due to the loaded image
+  vtkCamera* camera = renderer->GetActiveCamera( );
+  camera->SetViewUp( 0, -1, 0 );
+  camera->SetPosition( 0, 0, -1 );
+  camera->SetFocalPoint( 0, 0, 0 );
+
   // VTK interaction
   vtkSmartPointer< vtkInteractorStyleImage > imageStyle =
     vtkSmartPointer< vtkInteractorStyleImage >::New( );
@@ -202,7 +176,6 @@ int main( int argc, char* argv[] )
     vtkSmartPointer< vtkRenderWindowInteractor >::New( );
   interactor->SetInteractorStyle( imageStyle );
   window->SetInteractor( interactor );
-  window->Render( );
 
   // Create the widget and its representation
   vtkSmartPointer< vtkPointHandleRepresentation3D > handle =
@@ -219,44 +192,52 @@ int main( int argc, char* argv[] )
 
   // Let some interaction
   interactor->Initialize( );
+  renderer->ResetCamera( );
   window->Render( );
   widget->On( );
   interactor->Start( );
 
-  // Configure observer
-  TObserver::Pointer obs = TObserver::New( );
-  obs->SetImage( input_image, window );
-
-  // Configure membership function
-  TMembershipFunction::Pointer membership = TMembershipFunction::New( );
-  membership->SetModelSupport( support );
-  membership->SetInputImage( const_cast< TColorImage* >( input_image.GetPointer( ) ) );
-
-  // Configure algorithm
-  TFrontAlgorithm::Pointer algorithm = TFrontAlgorithm::New( );
+  // Prepare region grow function
+  typedef GaussianFunction< TColorImage, TScalar > TFunction;
+  TFunction::Pointer function = TFunction::New( );
+  function->SetModelSupport( support );
+
+  // Prepare region grow filter
+  typedef fpa::Image::RegionGrow< TColorImage, TImage > TFilter;
+  TFilter::Pointer filter = TFilter::New( );
+  filter->SetInput( input_image );
+  filter->SetMembershipFunction( function );
+  filter->SetNeighborhoodOrder( neighborhood_order );
+  filter->SetOutsideValue( TPixel( 0 ) );
+  filter->SetInsideValue( std::numeric_limits< TPixel >::max( ) );
+  filter->StopAtOneFrontOff( );
+
+  // Get user-given seeds
   for( unsigned int s = 0; s < rep->GetNumberOfSeeds( ); s++ )
   {
     double pos[ 3 ];
     rep->GetSeedWorldPosition( s, pos );
 
-    TColorImage::PointType pnt;
-    pnt[ 0 ] = TScalar( pos[ 0 ] );
-    pnt[ 1 ] = TScalar( pos[ 1 ] );
+    TImage::PointType pnt;
+    pnt[ 0 ] = TImage::PointType::ValueType( pos[ 0 ] );
+    pnt[ 1 ] = TImage::PointType::ValueType( pos[ 1 ] );
 
-    TColorImage::IndexType idx;
+    TImage::IndexType idx;
     if( input_image->TransformPhysicalPointToIndex( pnt, idx ) )
-      algorithm->AddSeed( idx, 0 );
+      filter->AddSeed( idx, 0 );
 
   } // rof
-  algorithm->AddObserver( itk::AnyEvent( ), obs );
-  algorithm->ThrowEventsOn( );
-  algorithm->SetInput( input_image );
-  algorithm->SetNeighborhoodOrder( 1 );
-  algorithm->SetMembershipFunction( membership );
-  algorithm->Update( );
-
-  // One last interaction
-  interactor->Start( );
+
+  // Prepare graphical debugger
+  typedef fpa::VTK::Image2DObserver< TFilter, vtkRenderWindow > TDebugger;
+  TDebugger::Pointer debugger = TDebugger::New( );
+  debugger->SetRenderWindow( window );
+  debugger->SetRenderPercentage( 0.001 );
+  filter->AddObserver( itk::AnyEvent( ), debugger );
+  filter->ThrowEventsOn( );
+
+  // Go!
+  filter->Update( );
 
   return( 0 );
 }
index 21670ef7396522c40c2a918432182d7d96037b34..eb307580601ef9a85042165a46cd93b934afa94f 100644 (file)
@@ -148,7 +148,7 @@ _Loop( )
 
       // Iterate over neighbors
       typename _TVertices::iterator nIt = neighborhood.begin( );
-      for( ; nIt != neighborhood.end( ); ++nIt )
+      while( nIt != neighborhood.end( ) )
       {
         _TNode neighbor = this->_Node( *nIt );
         neighbor.Vertex = *nIt;
@@ -159,6 +159,7 @@ _Loop( )
           {
             this->_QueueClear( );
             nIt = neighborhood.end( );
+            continue;
 
           } // fi
         }
@@ -182,8 +183,9 @@ _Loop( )
             this->InvokeEvent( TFreezeEvent( *nIt, candidate.FrontId ) );
 
         } // fi
+        ++nIt;
 
-      } // rof
+      } // elihw
     }
     else
       this->_QueueClear( );
@@ -230,13 +232,14 @@ _UpdateCollisions( const TVertex& a, const TVertex& b )
   _TNode na = this->_Node( a );
   _TNode nb = this->_Node( b );
   long fa = na.FrontId;
-  long fb = na.FrontId;
+  long fb = nb.FrontId;
 
   if( fa != fb )
   {
     // Mark collision, if it is new
     bool exists = this->m_Collisions[ fa ][ fb ].second;
     exists     &= this->m_Collisions[ fb ][ fa ].second;
+        
     if( !exists )
     {
       this->m_Collisions[ fa ][ fb ].first = na.Vertex;
diff --git a/lib/fpa/Base/MinimumSpanningTree.h b/lib/fpa/Base/MinimumSpanningTree.h
new file mode 100644 (file)
index 0000000..ec744e6
--- /dev/null
@@ -0,0 +1,69 @@
+#ifndef __FPA__BASE__MINIMUMSPANNINGTREE__H__
+#define __FPA__BASE__MINIMUMSPANNINGTREE__H__
+
+#include <vector>
+#include <itkSmartPointer.h>
+
+namespace fpa
+{
+  namespace Base
+  {
+    /**
+     */
+    template< class V, class C, class B >
+    class MinimumSpanningTree
+      : public B
+    {
+    public:
+      typedef MinimumSpanningTree             Self;
+      typedef B                               Superclass;
+      typedef itk::SmartPointer< Self >       Pointer;
+      typedef itk::SmartPointer< const Self > ConstPointer;
+
+      typedef V TVertex;
+      typedef C TCollisions;
+
+    protected:
+      typedef std::vector< unsigned long > _TRow;
+      typedef std::vector< _TRow >         _TMatrix;
+
+    public:
+      itkTypeMacro( MinimumSpanningTree, B );
+
+      itkGetConstMacro( Collisions, TCollisions );
+
+    public:
+      void SetCollisions( const TCollisions& collisions );
+
+      virtual void GetPath(
+        std::vector< V >& path, const V& a, const V& b
+        ) const;
+
+    protected:
+      MinimumSpanningTree( );
+      virtual ~MinimumSpanningTree( );
+
+      virtual void _Path( std::vector< V >& path, const V& a ) const;
+      virtual long _FrontId( const V& v ) const = 0;
+      virtual V _Parent( const V& v ) const = 0;
+
+    private:
+      // Purposely not implemented
+      MinimumSpanningTree( const Self& other );
+      Self& operator=( const Self& other );
+
+    protected:
+      TCollisions m_Collisions;
+      _TMatrix    m_FrontPaths;
+      static const unsigned long INF_VALUE;
+    };
+
+  } // ecapseman
+
+} // ecapseman
+
+#include <fpa/Base/MinimumSpanningTree.hxx>
+
+#endif // __FPA__BASE__MINIMUMSPANNINGTREE__H__
+
+// eof - $RCSfile$
diff --git a/lib/fpa/Base/MinimumSpanningTree.hxx b/lib/fpa/Base/MinimumSpanningTree.hxx
new file mode 100644 (file)
index 0000000..1140e40
--- /dev/null
@@ -0,0 +1,179 @@
+#ifndef __FPA__BASE__MINIMUMSPANNINGTREE__HXX__
+#define __FPA__BASE__MINIMUMSPANNINGTREE__HXX__
+
+#include <limits>
+
+// -------------------------------------------------------------------------
+template< class V, class C, class B >
+const unsigned long fpa::Base::MinimumSpanningTree< V, C, B >::INF_VALUE =
+  std::numeric_limits< unsigned long >::max( ) >> 1;
+
+// -------------------------------------------------------------------------
+template< class V, class C, class B >
+void fpa::Base::MinimumSpanningTree< V, C, B >::
+SetCollisions( const TCollisions& collisions )
+{
+  this->m_Collisions = collisions;
+
+  // Prepare fronts graph using Floyd-Warshall
+  unsigned long nSeeds = this->m_Collisions.size( );
+  _TMatrix dist( nSeeds, _TRow( nSeeds, Self::INF_VALUE ) );
+  this->m_FrontPaths = dist;
+
+  for( unsigned long i = 0; i < nSeeds; ++i )
+  {
+    for( unsigned long j = 0; j < nSeeds; ++j )
+    {
+      if( this->m_Collisions[ i ][ j ].second )
+      {
+        dist[ i ][ j ] = 1;
+        dist[ j ][ i ] = 1;
+        this->m_FrontPaths[ i ][ j ] = j;
+        this->m_FrontPaths[ j ][ i ] = i;
+
+      } // fi
+
+    } // rof
+    dist[ i ][ i ] = 0;
+    this->m_FrontPaths[ i ][ i ] = i;
+
+  } // rof
+  for( unsigned long k = 0; k < nSeeds; ++k )
+  {
+    for( unsigned long i = 0; i < nSeeds; ++i )
+    {
+      for( unsigned long j = 0; j < nSeeds; ++j )
+      {
+        if( ( dist[ i ][ k ] + dist[ k ][ j ] ) < dist[ i ][ j ] )
+        {
+          dist[ i ][ j ] = dist[ i ][ k ] + dist[ k ][ j ];
+          this->m_FrontPaths[ i ][ j ] = this->m_FrontPaths[ i ][ k ];
+
+        } // fi
+
+      } // rof
+
+    } // rof
+
+  } // rof
+
+  this->Modified( );
+}
+
+// -------------------------------------------------------------------------
+template< class V, class C, class B >
+void fpa::Base::MinimumSpanningTree< V, C, B >::
+GetPath( std::vector< V >& path, const V& a, const V& b ) const
+{
+  long fa = this->_FrontId( a );
+  long fb = this->_FrontId( b );
+
+  if( fa == fb )
+  {
+    std::vector< V > ap, bp;
+    this->_Path( ap, a );
+    this->_Path( bp, b );
+
+    // Ignore common part
+    typename std::vector< V >::const_reverse_iterator raIt, rbIt;
+    raIt = ap.rbegin( );
+    rbIt = bp.rbegin( );
+    while( *raIt == *rbIt && raIt != ap.rend( ) && rbIt != bp.rend( ) )
+    {
+      ++raIt;
+      ++rbIt;
+
+    } // elihw
+    if( raIt != ap.rbegin( ) ) --raIt;
+    if( rbIt != bp.rbegin( ) ) --rbIt;
+
+    // Add part from a
+    typename std::vector< V >::const_iterator iaIt = ap.begin( );
+    for( ; iaIt != ap.end( ) && *iaIt != *raIt; ++iaIt )
+      path.push_back( *iaIt );
+
+    // Add part from b
+    for( ; rbIt != bp.rend( ); ++rbIt )
+      path.push_back( *rbIt );
+  }
+  else
+  {
+    // Use this->m_FrontPaths from Floyd-Warshall
+    if( this->m_FrontPaths[ fa ][ fb ] != Self::INF_VALUE )
+    {
+      // Compute front path
+      std::vector< long > fpath;
+      fpath.push_back( fa );
+      while( fa != fb )
+      {
+        fa = this->m_FrontPaths[ fa ][ fb ];
+        fpath.push_back( fa );
+
+      } // elihw
+
+      // Continue only if both fronts are connected
+      unsigned int N = fpath.size( );
+      if( 0 < N )
+      {
+        // First path: from start vertex to first collision
+        this->GetPath(
+          path, a, this->m_Collisions[ fpath[ 0 ] ][ fpath[ 1 ] ].first
+          );
+
+        // Intermediary paths
+        for( unsigned int i = 1; i < N - 1; ++i )
+        {
+          this->GetPath(
+            path,
+            this->m_Collisions[ fpath[ i ] ][ fpath[ i - 1 ] ].first,
+            this->m_Collisions[ fpath[ i ] ][ fpath[ i + 1 ] ].first
+            );
+
+        } // rof
+
+        // Final path: from last collision to end point
+        this->GetPath(
+          path,
+          this->m_Collisions[ fpath[ N - 1 ] ][ fpath[ N - 2 ] ].first, b
+          );
+
+      } // fi
+
+    } // fi
+
+  } // fi
+}
+
+// -------------------------------------------------------------------------
+template< class V, class C, class B >
+fpa::Base::MinimumSpanningTree< V, C, B >::
+MinimumSpanningTree( )
+  : Superclass( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class V, class C, class B >
+fpa::Base::MinimumSpanningTree< V, C, B >::
+~MinimumSpanningTree( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class V, class C, class B >
+void fpa::Base::MinimumSpanningTree< V, C, B >::
+_Path( std::vector< V >& path, const V& a ) const
+{
+  V it = a;
+  do
+  {
+    path.push_back( it );
+    it = this->_Parent( it );
+
+  } while( it != this->_Parent( it ) );
+  path.push_back( it );
+}
+
+#endif // __FPA__BASE__MINIMUMSPANNINGTREE__HXX__
+
+// eof - $RCSfile$
diff --git a/lib/fpa/Base/TreeExtractor.h b/lib/fpa/Base/TreeExtractor.h
deleted file mode 100644 (file)
index ee6499b..0000000
+++ /dev/null
@@ -1,84 +0,0 @@
-#ifndef __FPA__BASE__TREEEXTRACTOR__H__
-#define __FPA__BASE__TREEEXTRACTOR__H__
-
-#include <map>
-#include <vector>
-
-namespace fpa
-{
-  namespace Base
-  {
-    /**
-     */
-    template< class A >
-    class TreeExtractor
-      : public A
-    {
-    public:
-      /// Standard class typdedefs
-      typedef TreeExtractor                   Self;
-      typedef A                               Superclass;
-      typedef itk::SmartPointer< Self >       Pointer;
-      typedef itk::SmartPointer< const Self > ConstPointer;
-
-      typedef A TAlgorithm;
-
-      typedef typename A::TCost                        TCost;
-      typedef typename A::TVertex                      TVertex;
-      typedef std::vector< TVertex >                   TVertices;
-      typedef typename A::TTraits::TVertexCmp          TVertexCmp;
-      typedef std::map< TVertex, TVertex, TVertexCmp > TTree;
-
-    protected:
-      typedef typename A::_TFrontId           _TFrontId;
-      typedef typename A::_TNode              _TNode;
-      typedef typename A::_TCollisionSitesRow _TCollisionSitesRow;
-      typedef typename A::_TCollisionSites    _TCollisionSites;
-      typedef std::vector< unsigned long >    _TRow;
-      typedef std::vector< _TRow >            _TMatrix;
-
-    public:
-      itkNewMacro( Self );
-      itkTypeMacro( TreeExtractor, TAlgorithm );
-
-      itkGetConstMacro( Root, TVertex );
-      itkGetConstMacro( Tree, TTree );
-      itkSetMacro( Root, TVertex );
-
-    public:
-      unsigned int AddLeaf( const TVertex& leaf );
-      unsigned int GetNumberOfLeafs( ) const;
-      const TVertex& GetLeaf( const unsigned int& id ) const;
-
-    protected:
-      TreeExtractor( );
-      virtual ~TreeExtractor( );
-
-      virtual void _BeforeMainLoop( );
-      virtual void _AfterMainLoop( );
-
-      TVertices _Path( const TVertex& s, const TVertex& e );
-
-    private:
-      // Purporsely not implemented
-      TreeExtractor( const Self& );
-      void operator=( const Self& );
-
-    protected:
-      TVertex   m_Root;
-      TVertices m_Leafs;
-      TTree     m_Tree;
-
-      _TMatrix m_FrontPaths;
-      static const unsigned long INF_VALUE;
-    };
-
-  } // ecapseman
-
-} // ecapseman
-
-#include <fpa/Base/TreeExtractor.hxx>
-
-#endif // __FPA__BASE__TREEEXTRACTOR__H__
-
-// eof - $RCSfile$
diff --git a/lib/fpa/Base/TreeExtractor.hxx b/lib/fpa/Base/TreeExtractor.hxx
deleted file mode 100644 (file)
index 7bad1ef..0000000
+++ /dev/null
@@ -1,250 +0,0 @@
-#ifndef __FPA__BASE__TREEEXTRACTOR__HXX__
-#define __FPA__BASE__TREEEXTRACTOR__HXX__
-
-#include <limits>
-
-// -------------------------------------------------------------------------
-template< class A >
-const unsigned long fpa::Base::TreeExtractor< A >::INF_VALUE =
-  std::numeric_limits< unsigned long >::max( ) >> 1;
-
-// -------------------------------------------------------------------------
-template< class A >
-unsigned int fpa::Base::TreeExtractor< A >::
-AddLeaf( const TVertex& leaf )
-{
-  this->m_Leafs.push_back( leaf );
-  this->Modified( );
-  return( this->m_Leafs.size( ) - 1 );
-}
-
-// -------------------------------------------------------------------------
-template< class A >
-unsigned int fpa::Base::TreeExtractor< A >::
-GetNumberOfLeafs( ) const
-{
-  return( this->m_Leafs.size( ) );
-}
-
-// -------------------------------------------------------------------------
-template< class A >
-const typename fpa::Base::TreeExtractor< A >::
-TVertex& fpa::Base::TreeExtractor< A >::
-GetLeaf( const unsigned int& id ) const
-{
-  static const TVertex zero;
-  if( id < this->m_Leafs.size( ) )
-    return( this->m_Leafs[ id ] );
-  else
-    return( zero );
-}
-
-// -------------------------------------------------------------------------
-template< class A >
-fpa::Base::TreeExtractor< A >::
-TreeExtractor( )
-  : Superclass( )
-{
-}
-
-// -------------------------------------------------------------------------
-template< class A >
-fpa::Base::TreeExtractor< A >::
-~TreeExtractor( )
-{
-}
-
-// -------------------------------------------------------------------------
-template< class A >
-void fpa::Base::TreeExtractor< A >::
-_BeforeMainLoop( )
-{
-  this->Superclass::_BeforeMainLoop( );
-  this->m_Tree.clear( );
-}
-
-// -------------------------------------------------------------------------
-template< class A >
-void fpa::Base::TreeExtractor< A >::
-_AfterMainLoop( )
-{
-  this->Superclass::_AfterMainLoop( );
-  
-  // Prepare fronts graph using Floyd-Warshall
-  unsigned long nSeeds = this->GetNumberOfSeeds( );
-  _TMatrix dist( nSeeds, _TRow( nSeeds, Self::INF_VALUE ) );
-  this->m_FrontPaths = dist;
-
-  for( unsigned long i = 0; i < nSeeds; ++i )
-  {
-    for( unsigned long j = 0; j < nSeeds; ++j )
-    {
-      if( this->m_CollisionSites[ i ][ j ].second )
-      {
-        dist[ i ][ j ] = 1;
-        dist[ j ][ i ] = 1;
-        this->m_FrontPaths[ i ][ j ] = j;
-        this->m_FrontPaths[ j ][ i ] = i;
-
-      } // fi
-
-    } // rof
-    dist[ i ][ i ] = 0;
-    this->m_FrontPaths[ i ][ i ] = i;
-
-  } // rof
-  for( unsigned long k = 0; k < nSeeds; ++k )
-  {
-    for( unsigned long i = 0; i < nSeeds; ++i )
-    {
-      for( unsigned long j = 0; j < nSeeds; ++j )
-      {
-        if( ( dist[ i ][ k ] + dist[ k ][ j ] ) < dist[ i ][ j ] )
-        {
-          dist[ i ][ j ] = dist[ i ][ k ] + dist[ k ][ j ];
-          this->m_FrontPaths[ i ][ j ] = this->m_FrontPaths[ i ][ k ];
-
-        } // fi
-
-      } // rof
-
-    } // rof
-
-  } // rof
-
-  // Build tree
-  this->m_Tree.clear( );
-  typename TVertices::const_iterator vIt = this->m_Leafs.begin( );
-  for( ; vIt != this->m_Leafs.end( ); ++vIt )
-  {
-    TVertices path = this->_Path( this->m_Root, *vIt );
-    typename TVertices::const_iterator pIt = path.begin( );
-    TVertex parent = *pIt;
-    for( ; pIt != path.end( ); ++pIt )
-    {
-      this->m_Tree[ *pIt ] = parent;
-      parent = *pIt;
-
-    } // rof
-
-  } // rof
-}
-
-// -------------------------------------------------------------------------
-template< class A >
-typename fpa::Base::TreeExtractor< A >::
-TVertices fpa::Base::TreeExtractor< A >::
-_Path( const TVertex& s, const TVertex& e )
-{
-  TVertices pS, pE;
-
-  // Check if both vertices belong to the same front or if collisions
-  // should be analyzed.
-  _TFrontId fs = this->_FrontId( s );
-  _TFrontId fe = this->_FrontId( e );
-  if( fs == fe )
-  {
-    // Path from start vertex until seed
-    TVertex vIt = s;
-    TVertex parent = this->_Parent( vIt );
-    while( parent != vIt )
-    {
-      pS.push_back( vIt );
-      vIt = parent;
-      parent = this->_Parent( vIt );
-
-    } // elihw
-    pS.push_back( vIt );
-
-    // Path from end vertex until seed
-    vIt = e;
-    parent = this->_Parent( vIt );
-    while( parent != vIt )
-    {
-      pE.push_back( vIt );
-      vIt = parent;
-      parent = this->_Parent( vIt );
-
-    } // elihw
-    pE.push_back( vIt );
-
-    // Erase duplicate vertices
-    if( pS.size( ) > 0 && pE.size( ) > 0 )
-    {
-      bool cont = ( pS.back( ) == pE.back( ) );
-      bool enter = false;
-      TVertex lastVertex = pS.back( );
-      while( cont )
-      {
-        enter = true;
-
-        lastVertex = pS.back( );
-        pS.pop_back( );
-        pE.pop_back( );
-
-        cont = ( pS.size( ) > 0 && pE.size( ) > 0 );
-        if( cont )
-          cont = ( pS.back( ) == pE.back( ) );
-
-      } // elihw
-      if( enter )
-        pS.push_back( lastVertex );
-
-    } // fi
-
-    // Contatenate both paths
-    pS.insert( pS.end( ), pE.rbegin( ), pE.rend( ) );
-  }
-  else
-  {
-    // Use this->m_FrontPaths from Floyd-Warshall
-    if( this->m_FrontPaths[ fs ][ fe ] != Self::INF_VALUE )
-    {
-      // Compute front path
-      std::vector< unsigned long > fpath;
-      fpath.push_back( fs );
-      while( fs != fe )
-      {
-        fs = this->m_FrontPaths[ fs ][ fe ];
-        fpath.push_back( fs );
-
-      } // elihw
-
-      // Continue only if both fronts are connected
-      unsigned int N = fpath.size( );
-      if( 0 < N )
-      {
-        // First path: from start vertex to first collision
-        pS = this->_Path(
-          s, this->m_CollisionSites[ fpath[ 0 ] ][ fpath[ 1 ] ].first
-          );
-
-        // Intermediary paths
-        for( unsigned int i = 1; i < N - 1; ++i )
-        {
-          pE = this->_Path(
-            this->m_CollisionSites[ fpath[ i ] ][ fpath[ i - 1 ] ].first,
-            this->m_CollisionSites[ fpath[ i ] ][ fpath[ i + 1 ] ].first
-            );
-          pS.insert( pS.end( ), pE.begin( ), pE.end( ) );
-
-        } // rof
-
-        // Final path: from last collision to end point
-        pE = this->_Path(
-          this->m_CollisionSites[ fpath[ N - 1 ] ][ fpath[ N - 2 ] ].first,
-          e
-          );
-        pS.insert( pS.end( ), pE.begin( ), pE.end( ) );
-
-      } // fi
-
-    } // fi
-
-  } // fi
-  return( pS );
-}
-
-#endif // __FPA__BASE__TREEEXTRACTOR__HXX__
-
-// eof - $RCSfile$
index cd218a6ff0d9ef1eaae6694206a27a37d126237b..e59bed2311c19912e26c6d75df7aa6c91079f7d9 100644 (file)
@@ -2,6 +2,7 @@
 #define __FPA__IMAGE__ALGORITHM__H__
 
 #include <itkImage.h>
+#include <fpa/Image/MinimumSpanningTree.h>
 
 namespace fpa
 {
@@ -41,7 +42,10 @@ namespace fpa
       typedef typename Superclass::_TNode          _TNode;
       typedef typename Superclass::_TNodes         _TNodes;
 
-      typedef itk::Image< _TNode, I::ImageDimension > _TMarks;
+      typedef fpa::Image::MinimumSpanningTree< TVertex, _TNode, _TCollisions, I::ImageDimension, Self::AliveLabel > _TMarks;
+
+    public:
+      typedef _TMarks TMinimumSpanningTree;
 
     public:
       itkTypeMacro( Algorithm, TAlgorithm );
@@ -50,11 +54,17 @@ namespace fpa
       itkGetConstMacro( NeighborhoodOrder, unsigned int );
       itkSetMacro( NeighborhoodOrder, unsigned int );
 
+    public:
+      TMinimumSpanningTree* GetMinimumSpanningTree( );
+      const TMinimumSpanningTree* GetMinimumSpanningTree( ) const;
+      void GraftMinimumSpanningTree( itk::DataObject* obj );
+
     protected:
       Algorithm( );
       virtual ~Algorithm( );
 
       virtual void _BeforeGenerateData( );
+      virtual void _AfterGenerateData( );
 
       // Graph-related abstract methods
       virtual unsigned long _NumberOfVertices( ) const;
@@ -84,7 +94,6 @@ namespace fpa
 
     protected:
       unsigned int m_NeighborhoodOrder;
-      typename _TMarks::Pointer m_Marks;
     };
 
   } // ecapseman
index 0d83b035e6f821ccb52fe9d3bc0a843f0c72ebe7..a5b1f8548e5e4dd481a93eb1fbb3e3f6561d4ff4 100644 (file)
@@ -4,6 +4,42 @@
 #include <cmath>
 #include <itkConstNeighborhoodIterator.h>
 
+// -------------------------------------------------------------------------
+template< class I, class O, class A >
+typename fpa::Image::Algorithm< I, O, A >::
+TMinimumSpanningTree* fpa::Image::Algorithm< I, O, A >::
+GetMinimumSpanningTree( )
+{
+  return(
+    dynamic_cast< TMinimumSpanningTree* >(
+      this->itk::ProcessObject::GetOutput( 1 )
+      )
+    );
+}
+
+// -------------------------------------------------------------------------
+template< class I, class O, class A >
+const typename fpa::Image::Algorithm< I, O, A >::
+TMinimumSpanningTree* fpa::Image::Algorithm< I, O, A >::
+GetMinimumSpanningTree( ) const
+{
+  return(
+    dynamic_cast< const TMinimumSpanningTree* >(
+      this->itk::ProcessObject::GetOutput( 1 )
+      )
+    );
+}
+
+// -------------------------------------------------------------------------
+template< class I, class O, class A >
+void fpa::Image::Algorithm< I, O, A >::
+GraftMinimumSpanningTree( itk::DataObject* obj )
+{
+  TMinimumSpanningTree* mst = dynamic_cast< TMinimumSpanningTree* >( obj );
+  if( mst != NULL )
+    this->GraftNthOutput( 1, mst );
+}
+
 // -------------------------------------------------------------------------
 template< class I, class O, class A >
 fpa::Image::Algorithm< I, O, A >::
@@ -11,6 +47,9 @@ Algorithm( )
   : Superclass( ),
     m_NeighborhoodOrder( 1 )
 {
+  this->itk::ProcessObject::SetNumberOfRequiredOutputs( 2 );
+  this->itk::ProcessObject::SetNthOutput( 0, O::New( ) );
+  this->itk::ProcessObject::SetNthOutput( 1, TMinimumSpanningTree::New( ) );
 }
 
 // -------------------------------------------------------------------------
@@ -29,6 +68,15 @@ _BeforeGenerateData( )
   this->AllocateOutputs( );
 }
 
+// -------------------------------------------------------------------------
+template< class I, class O, class A >
+void fpa::Image::Algorithm< I, O, A >::
+_AfterGenerateData( )
+{
+  this->Superclass::_AfterGenerateData( );
+  this->GetMinimumSpanningTree( )->SetCollisions( this->m_Collisions );
+}
+
 // -------------------------------------------------------------------------
 template< class I, class O, class A >
 unsigned long fpa::Image::Algorithm< I, O, A >::
@@ -145,7 +193,7 @@ const typename fpa::Image::Algorithm< I, O, A >::
 _TNode& fpa::Image::Algorithm< I, O, A >::
 _Node( const TVertex& v ) const
 {
-  return( this->m_Marks->GetPixel( v ) );
+  return( this->GetMinimumSpanningTree( )->GetPixel( v ) );
 }
 
 // -------------------------------------------------------------------------
@@ -153,20 +201,9 @@ template< class I, class O, class A >
 void fpa::Image::Algorithm< I, O, A >::
 _InitMarks( )
 {
-  const I* in = this->GetInput( );
-
-  this->m_Marks = _TMarks::New( );
-  this->m_Marks->SetLargestPossibleRegion( in->GetLargestPossibleRegion( ) );
-  this->m_Marks->SetRequestedRegion( in->GetRequestedRegion( ) );
-  this->m_Marks->SetBufferedRegion( in->GetBufferedRegion( ) );
-  this->m_Marks->SetOrigin( in->GetOrigin( ) );
-  this->m_Marks->SetSpacing( in->GetSpacing( ) );
-  this->m_Marks->SetDirection( in->GetDirection( ) );
-  this->m_Marks->Allocate( );
-
   _TNode far_node;
   far_node.Label = Self::FarLabel;
-  this->m_Marks->FillBuffer( far_node );
+  this->GetMinimumSpanningTree( )->FillBuffer( far_node );
 }
 
 // -------------------------------------------------------------------------
@@ -174,7 +211,7 @@ template< class I, class O, class A >
 void fpa::Image::Algorithm< I, O, A >::
 _Mark( const _TNode& node )
 {
-  this->m_Marks->SetPixel( node.Vertex, node );
+  this->GetMinimumSpanningTree( )->SetPixel( node.Vertex, node );
 }
 
 #endif // __FPA__IMAGE__ALGORITHM__HXX__
diff --git a/lib/fpa/Image/MinimumSpanningTree.h b/lib/fpa/Image/MinimumSpanningTree.h
new file mode 100644 (file)
index 0000000..652547d
--- /dev/null
@@ -0,0 +1,54 @@
+#ifndef __FPA__IMAGE__MINIMUMSPANNINGTREE__H__
+#define __FPA__IMAGE__MINIMUMSPANNINGTREE__H__
+
+#include <itkImage.h>
+#include <fpa/Base/MinimumSpanningTree.h>
+
+namespace fpa
+{
+  namespace Image
+  {
+    /**
+     */
+    template< class V, class N, class C, unsigned int D, unsigned long L >
+    class MinimumSpanningTree
+      : public fpa::Base::MinimumSpanningTree< V, C, itk::Image< N, D > >
+    {
+    public:
+      typedef fpa::Base::MinimumSpanningTree< V, C, itk::Image< N, D > > Superclass;
+      typedef MinimumSpanningTree             Self;
+      typedef itk::SmartPointer< Self >       Pointer;
+      typedef itk::SmartPointer< const Self > ConstPointer;
+
+      typedef V TVertex;
+      typedef C TCollisions;
+
+    protected:
+      typedef N _TNode;
+
+    public:
+      itkNewMacro( Self );
+      itkTypeMacro( MinimumSpanningTree, fpa_Base_MinimumSpanningTree );
+
+    protected:
+      MinimumSpanningTree( );
+      virtual ~MinimumSpanningTree( );
+
+      virtual long _FrontId( const V& v ) const;
+      virtual V _Parent( const V& v ) const;
+
+    private:
+      // Purposely not implemented
+      MinimumSpanningTree( const Self& other );
+      Self& operator=( const Self& other );
+    };
+
+  } // ecapseman
+
+} // ecapseman
+
+#include <fpa/Image/MinimumSpanningTree.hxx>
+
+#endif // __FPA__IMAGE__MINIMUMSPANNINGTREE__H__
+
+// eof - $RCSfile$
diff --git a/lib/fpa/Image/MinimumSpanningTree.hxx b/lib/fpa/Image/MinimumSpanningTree.hxx
new file mode 100644 (file)
index 0000000..48a43cb
--- /dev/null
@@ -0,0 +1,41 @@
+#ifndef __FPA__IMAGE__MINIMUMSPANNINGTREE__HXX__
+#define __FPA__IMAGE__MINIMUMSPANNINGTREE__HXX__
+
+// -------------------------------------------------------------------------
+template< class V, class N, class C, unsigned int D, unsigned long L >
+fpa::Image::MinimumSpanningTree< V, N, C, D, L >::
+MinimumSpanningTree( )
+  : Superclass( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class V, class N, class C, unsigned int D, unsigned long L >
+fpa::Image::MinimumSpanningTree< V, N, C, D, L >::
+~MinimumSpanningTree( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class V, class N, class C, unsigned int D, unsigned long L >
+long fpa::Image::MinimumSpanningTree< V, N, C, D, L >::
+_FrontId( const V& v ) const
+{
+  return( this->GetPixel( v ).FrontId );
+}
+
+// -------------------------------------------------------------------------
+template< class V, class N, class C, unsigned int D, unsigned long L >
+V fpa::Image::MinimumSpanningTree< V, N, C, D, L >::
+_Parent( const V& v ) const
+{
+  _TNode n = this->GetPixel( v );
+  if( n.Label == L )
+    return( n.Parent );
+  else
+    return( v );
+}
+
+#endif // __FPA__IMAGE__MINIMUMSPANNINGTREE__HXX__
+
+// eof - $RCSfile$