...

[FrontAlgorithms.git] / lib / fpa / Image / ExtractEndPointsAndBifurcationsFromMinimumSpanningTree.hxx

diff --git a/lib/fpa/Image/ExtractEndPointsAndBifurcationsFromMinimumSpanningTree.hxx b/lib/fpa/Image/ExtractEndPointsAndBifurcationsFromMinimumSpanningTree.hxx
index 0a33ee98fdf27dfc4f80952dbcf45122b427cdc4..789d392f6a65634e186ca372921404559c31bc6b 100644 (file)
--- a/lib/fpa/Image/ExtractEndPointsAndBifurcationsFromMinimumSpanningTree.hxx
+++ b/lib/fpa/Image/ExtractEndPointsAndBifurcationsFromMinimumSpanningTree.hxx
@@ -2,7 +2,7 @@
 #define __FPA__IMAGE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__
 
 #include <cmath>
-#include <cpPlugins_ITKInstances/ImageIterators.h>
+#include <itkImageRegionIteratorWithIndex.h>
 
 // -------------------------------------------------------------------------
 template< class _TImage, class _TMST >
@@ -30,6 +30,23 @@ GetDistanceMap( )
     );
 }
 
+// -------------------------------------------------------------------------
+template< class _TImage, class _TMST >
+typename
+fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
+TSkeleton*
+fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
+GetSkeleton( )
+{
+  return(
+    dynamic_cast< TSkeleton* >(
+      this->itk::ProcessObject::GetOutput(
+        this->GetNumberOfRequiredOutputs( ) - 1
+        )
+      )
+    );
+}
+
 // -------------------------------------------------------------------------
 template< class _TImage, class _TMST >
 void
@@ -56,6 +73,10 @@ ExtractEndPointsAndBifurcationsFromMinimumSpanningTree( )
     m_SquaredDistanceMap( false )
 {
   this->SetNumberOfRequiredInputs( 3 );
+  unsigned int nOuts = this->GetNumberOfRequiredOutputs( );
+  this->SetNumberOfRequiredOutputs( nOuts + 1 );
+  typename TSkeleton::Pointer sk = TSkeleton::New( );
+  this->itk::ProcessObject::SetNthOutput( nOuts, sk.GetPointer( ) );
 }
 
 // -------------------------------------------------------------------------
@@ -71,8 +92,8 @@ void
 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
 GenerateData( )
 {
+  // Create auxiliary objects
   auto image = this->GetCostsImage( );
-
   this->m_MarkImage = TMarkImage::New( );
   this->m_MarkImage->SetLargestPossibleRegion( image->GetLargestPossibleRegion( ) );
   this->m_MarkImage->SetRequestedRegion( image->GetRequestedRegion( ) );
@@ -82,7 +103,6 @@ GenerateData( )
   this->m_MarkImage->SetSpacing( image->GetSpacing( ) );
   this->m_MarkImage->Allocate( );
   this->m_MarkImage->FillBuffer( 0 );
-
   this->m_SkeletonImage = TMarkImage::New( );
   this->m_SkeletonImage->SetLargestPossibleRegion( image->GetLargestPossibleRegion( ) );
   this->m_SkeletonImage->SetRequestedRegion( image->GetRequestedRegion( ) );
@@ -93,7 +113,15 @@ GenerateData( )
   this->m_SkeletonImage->Allocate( );
   this->m_SkeletonImage->FillBuffer( 0 );
 
+  // Real execution
   this->Superclass::GenerateData( );
+
+  // Build skeleton
+  auto sk = this->GetSkeleton( );
+  sk->SetMinimumSpanningTree( this->GetMinimumSpanningTree( ) );
+  auto bIt = this->m_Branches.begin( );
+  for( ; bIt != this->m_Branches.end( ); ++bIt )
+    sk->AddBranch( bIt->second, bIt->first );
 }
 
 // -------------------------------------------------------------------------
@@ -116,6 +144,8 @@ _MarkSphere(
   const unsigned long& l
   )
 {
+  double rr = r * double( 1.5 );
+
   // Get marking region
   auto rreg = this->m_MarkImage->GetRequestedRegion( );
   auto spac = this->m_MarkImage->GetSpacing( );
@@ -126,7 +156,7 @@ _MarkSphere(
   bool all = false;
   for( unsigned int d = 0; d < _TImage::ImageDimension; ++d )
   {
-    unsigned long s = std::ceil( r / double( spac[ d ] ) );
+    unsigned long s = std::ceil( rr / double( spac[ d ] ) );
     if( s < 3 )
     {
       s = 3;
@@ -158,9 +188,11 @@ _MarkSphere(
   {
     typename _TImage::PointType pnt;
     this->m_MarkImage->TransformIndexToPhysicalPoint( spIt.GetIndex( ), pnt );
-    // TODO: if( double( pnt.EuclideanDistanceTo( c ) ) <= r || all )
-    // TODO: if( spIt.Get( ) == 0 )
-      spIt.Set( l );
+    /*
+      if( double( pnt.EuclideanDistanceTo( c ) ) <= rr || all )
+      if( spIt.Get( ) == 0 )
+    */
+    spIt.Set( l );
 
   } // rof
 }