Test

diff --git a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx
index 11fa961e078b81c35c49dca1e94fbd2be91eb12e..f2d140be01904e105c78c80299241cd3cbf1ab25 100644 (file)
--- a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx
+++ b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx
@@ -112,7 +112,8 @@ protected:
 
   virtual void _AfterMainLoop( )
     {
-      typedef itk::Image< unsigned char, TImage::ImageDimension > _TMarkImage;
+      typedef unsigned char                                _TMark;
+      typedef itk::Image< _TMark, TImage::ImageDimension > _TMarkImage;
 
       this->Superclass::_AfterMainLoop( );
       if( this->m_Candidates.size( ) == 0 )
@@ -121,6 +122,7 @@ protected:
       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( ) );
@@ -129,93 +131,97 @@ protected:
       marks->SetSpacing( spacing );
       marks->SetDirection( input->GetDirection( ) );
       marks->Allocate( );
-      marks->FillBuffer( 0 );
-      
-      this->m_Axes = vtkSmartPointer< vtkPolyData >::New( );
-      vtkSmartPointer< vtkPoints > points =
-        vtkSmartPointer< vtkPoints >::New( );
-      vtkSmartPointer< vtkCellArray > lines =
-        vtkSmartPointer< vtkCellArray >::New( );
+      marks->FillBuffer( _TMark( 0 ) );
 
+      // Iterate over the candidates, starting fromt the candidates that
+      // are near thin branches
       _TCandidates::const_reverse_iterator cIt =
         this->m_Candidates.rbegin( );
-
-      for( unsigned long leo = 0; leo < 50 && cIt != this->m_Candidates.rend( ); ++cIt )
+      for( int leo = 0; leo < 100 && cIt != this->m_Candidates.rend( ); ++cIt )
       {
-        TVertex vIt = cIt->second;
-        if( marks->GetPixel( vIt ) != 0 )
+        // If pixel hasn't been visited, continue
+        TVertex v = cIt->second;
+        if( marks->GetPixel( v ) != _TMark( 0 ) )
           continue;
 
-        leo++;
-        std::cout << "Leaf: " << leo << " " << cIt->first << " " << vIt << std::endl;
-        bool start = true;
-        do
-        {
-          double dist = double( input->GetPixel( vIt ) );
+        // Compute nearest start candidate
         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 );
+        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 ];
+        TImage::PixelType min_value = iIt.GetPixel( 0 );
+        TVertex min_vertex = iIt.GetIndex( 0 );
+        if( !( min_value > TImage::PixelType( 0 ) ) )
+          min_value = std::numeric_limits< TImage::PixelType >::max( );
+        for( unsigned int i = 1; i < nN; ++i )
+        {
+          TImage::PixelType value = iIt.GetPixel( i );
+          if( !( value > TImage::PixelType( 0 ) ) )
+            value = std::numeric_limits< TImage::PixelType >::max( );
+          if( value < min_value )
+          {
+            min_value = value;
+            min_vertex = iIt.GetIndex( i );
 
-          start = false;
+          } // fi
 
-          // Next vertex
-          vIt = this->_Parent( vIt );
+        } // rof
 
-        } while( vIt != this->_Parent( vIt ) );
+        if( min_value < std::numeric_limits< TImage::PixelType >::max( ) )
+        {
+          if( marks->GetPixel( min_vertex ) != _TMark( 0 ) )
+            continue;
+          leo++;
+          std::cout << "Leaf: " << leo << " " << min_value << " " << min_vertex << std::endl;
 
-        // Last vertex
-        // std::cout << vIt << " " << this->GetSeed( 0 ) << std::endl;
+          bool start = true;
+          do
+          {
+            double dist = double( 1.5 ) * std::sqrt( double( input->GetPixel( min_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( min_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 );
 
-        /* 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 );
-        */
+              start = false;
+            */
+            // Next vertex in current path
+            min_vertex = this->_Parent( min_vertex );
+
+          } while( min_vertex != this->_Parent( min_vertex ) );
+        }
+        else
+          marks->SetPixel( v, _TMark( 1 ) );
 
       } // rof
 
@@ -224,7 +230,92 @@ protected:
       w->SetInput( marks );
       w->SetFileName( "marks.mhd" );
       w->Update( );
+      std::exit( 1 );
 
+      /*
+
+        this->m_Axes = vtkSmartPointer< vtkPolyData >::New( );
+        vtkSmartPointer< vtkPoints > points =
+        vtkSmartPointer< vtkPoints >::New( );
+        vtkSmartPointer< vtkCellArray > lines =
+        vtkSmartPointer< vtkCellArray >::New( );
+
+        {
+
+        leo++;
+        std::cout << "Leaf: " << leo << " " << 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 )
@@ -246,13 +337,13 @@ protected:
       if( ret )
       {
         TCost nc = this->_Cost( n.Vertex, n.Parent );
-        TCost nc2 = nc * nc * nc * nc;
-        if( TCost( 0 ) < nc2 )
+        if( TCost( 0 ) < nc )
         {
-          n.Result = n.Cost / nc2;
-          this->GetOutput( )->SetPixel( n.Vertex, n.Result );
-          this->m_Candidates.insert( _TCandidate( n.Result, n.Vertex ) );
-
+          TCost cc = n.Cost / nc;
+          this->m_Candidates.insert( _TCandidate( cc, n.Vertex ) );
+          /* TODO: debug code
+             this->GetOutput( )->SetPixel( n.Vertex, cc );
+          */
         } // fi
 
       } // fi

diff --git a/appli/examples/example_ImageAlgorithm_Skeletonization.cxx b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx
index f2b82794894656dd0aa8a3158e47257cf5efb49e..98023c5353fec0bae3708ffd7045a55df4cbd29c 100644 (file)
--- a/appli/examples/example_ImageAlgorithm_Skeletonization.cxx
+++ b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx
@@ -124,6 +124,7 @@ int main( int argc, char* argv[] )
   TDistanceFilter::Pointer distanceMap = TDistanceFilter::New( );
   distanceMap->SetInput( segmentor->GetOutput( ) );
   distanceMap->InputIsBinaryOn( );
+  distanceMap->SquaredDistanceOn( );
   start = std::clock( );
   distanceMap->Update( );
   end = std::clock( );