Some tests

diff --git a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx
index 46bfe03442d7dfe835c2e64b34367396ca51f8bf..11fa961e078b81c35c49dca1e94fbd2be91eb12e 100644 (file)
--- a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx
+++ b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx
@@ -1,11 +1,13 @@
+#include <cmath>
 #include <ctime>
+#include <deque>
 #include <iostream>
 #include <limits>
 #include <map>
 #include <string>
-#include <deque>
 
 #include <itkConstNeighborhoodIterator.h>
+#include <itkNeighborhoodIterator.h>
 #include <itkDanielssonDistanceMapImageFilter.h>
 #include <itkImage.h>
 #include <itkImageFileReader.h>
@@ -110,9 +112,24 @@ protected:
 
   virtual void _AfterMainLoop( )
     {
+      typedef itk::Image< unsigned char, TImage::ImageDimension > _TMarkImage;
+
       this->Superclass::_AfterMainLoop( );
       if( this->m_Candidates.size( ) == 0 )
         return;
+
+      const TImage* input = this->GetInput( );
+      TImage::SpacingType spacing = input->GetSpacing( );
+
+      _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( 0 );
       
       this->m_Axes = vtkSmartPointer< vtkPolyData >::New( );
       vtkSmartPointer< vtkPoints > points =
@@ -123,41 +140,91 @@ protected:
       _TCandidates::const_reverse_iterator cIt =
         this->m_Candidates.rbegin( );
 
-      TVertices pS;
-      TVertex vIt = cIt->second;
-      TVertex parent = this->_Parent( vIt );
-      while( parent != vIt )
+      for( unsigned long leo = 0; leo < 50 && cIt != this->m_Candidates.rend( ); ++cIt )
       {
-        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 )
+        TVertex vIt = cIt->second;
+        if( marks->GetPixel( vIt ) != 0 )
+          continue;
+
+        leo++;
+        std::cout << "Leaf: " << leo << " " << cIt->first << " " << vIt << std::endl;
+        bool start = true;
+        do
         {
-          lines->InsertNextCell( 2 );
-          lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
-          lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
+          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 );
+        */
 
-        } // 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 );
+      } // rof
 
-        } // fi
+      itk::ImageFileWriter< _TMarkImage >::Pointer w =
+        itk::ImageFileWriter< _TMarkImage >::New( );
+      w->SetInput( marks );
+      w->SetFileName( "marks.mhd" );
+      w->Update( );
 
-        this->m_Axes->SetPoints( points );
-        this->m_Axes->SetLines( lines );
     }
 
   virtual bool _UpdateNeigh( _TNode& nn, const _TNode& n )
@@ -179,9 +246,10 @@ protected:
       if( ret )
       {
         TCost nc = this->_Cost( n.Vertex, n.Parent );
-        if( TCost( 0 ) < nc )
+        TCost nc2 = nc * nc * nc * nc;
+        if( TCost( 0 ) < nc2 )
         {
-          n.Result = n.Cost / ( nc * nc );
+          n.Result = n.Cost / nc2;
           this->GetOutput( )->SetPixel( n.Vertex, n.Result );
           this->m_Candidates.insert( _TCandidate( n.Result, n.Vertex ) );