Big bug stamped on wall

[FrontAlgorithms.git] / appli / examples / example_ImageAlgorithmDijkstra_03.cxx

diff --git a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx
index 613378b2c1e8db3a6aa130b99334e24a9ce9fce9..5fe6b24c27ea109d329f71d895dc62682b3b48a8 100644 (file)
--- a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx
+++ b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx
@@ -2,6 +2,7 @@
 #include <ctime>
 #include <deque>
 #include <iostream>
+#include <iomanip>
 #include <limits>
 #include <map>
 #include <string>
@@ -44,6 +45,19 @@ typedef fpa::Image::DijkstraWithSphereBacktracking< TImage, TScalar > TDijkstra;
 typedef fpa::VTK::ImageMPR TMPR;
 typedef fpa::VTK::Image3DObserver< TDijkstra, vtkRenderWindow > TDijkstraObs;
 
+struct TBranch
+{
+  double Length;
+  TImage::PointType V1;
+  TImage::PointType V2;
+
+  bool operator<( const TBranch& other ) const
+    {
+      return( other.Length < this->Length );
+    }
+};
+typedef std::set< TBranch > TBranches;
+
 // -------------------------------------------------------------------------
 int main( int argc, char* argv[] )
 {
@@ -92,12 +106,14 @@ int main( int argc, char* argv[] )
   view.SetSize( 800, 800 );
   view.SetImage( vtk_image->GetOutput( ) );
 
+  /*
   vtkSmartPointer< vtkImageMarchingCubes > mc =
     vtkSmartPointer< vtkImageMarchingCubes >::New( );
   mc->SetInputData( vtk_image->GetOutput( ) );
   mc->SetValue( 0, 1e-2 );
   mc->Update( );
   view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 );
+  */
 
   // Allow some interaction
   view.Render( );
@@ -143,47 +159,97 @@ int main( int argc, char* argv[] )
   const TDijkstra::TMarkImage* marks = paths->GetOutputMarkImage( );
   TDijkstra::TMark max_mark = paths->GetNumberOfBranches( );
   const TDijkstra::TVertices& endpoints = paths->GetEndPoints( );
+  const TDijkstra::TVertices& bifurcations = paths->GetBifurcationPoints( );
   const TDijkstra::TTree& tree = paths->GetFinalTree( );
   TDijkstra::TVertices::const_iterator epIt = endpoints.begin( );
-  for( unsigned int epId = 0; epIt != endpoints.end( ); ++epIt, ++epId )
+
+  TDijkstra::TMarkImage::Pointer new_marks = TDijkstra::TMarkImage::New( );
+  new_marks->SetLargestPossibleRegion( marks->GetLargestPossibleRegion( ) );
+  new_marks->SetRequestedRegion( marks->GetRequestedRegion( ) );
+  new_marks->SetBufferedRegion( marks->GetBufferedRegion( ) );
+  new_marks->SetOrigin( marks->GetOrigin( ) );
+  new_marks->SetSpacing( marks->GetSpacing( ) );
+  new_marks->SetDirection( marks->GetDirection( ) );
+  new_marks->Allocate( );
+  new_marks->FillBuffer( 0 );
+
+  std::cout << max_mark << std::endl;
+  std::cout << endpoints.size( ) << std::endl;
+
+  TBranches branches;
+  TBranch branch;
+  for( ; epIt != endpoints.end( ); ++epIt )
   {
     TDijkstra::TVertex idx = *epIt;
+    TDijkstra::TTree::const_iterator tIt = tree.find( idx );
+    TImage::PointType pnt, prev;
+    input_image->TransformIndexToPhysicalPoint( idx, pnt );
     do
     {
-      TImage::PointType pnt;
       input_image->TransformIndexToPhysicalPoint( idx, pnt );
-
-      TDijkstra::TMark mark = marks->GetPixel( idx );
-      double pd = double( mark );
+      branch.V2 = pnt;
 
       points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
-      scalars->InsertNextTuple1( pd );
+      scalars->InsertNextTuple1( double( tIt->second.second ) );
+      new_marks->SetPixel( idx, tIt->second.second );
+
       if( idx != *epIt )
       {
         cells->InsertNextCell( 2 );
         cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 );
         cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 );
+        branch.Length += prev.EuclideanDistanceTo( pnt );
 
       } // fi
-      idx = tree.find( idx )->second;
+      prev = pnt;
+      idx = tIt->second.first;
+
+      if( std::find( bifurcations.begin( ), bifurcations.end( ), idx ) != bifurcations.end( ) )
+      {
+        branches.insert( branch );
+        prev = pnt;
+        branch.V1 = pnt;
+        branch.Length = double( 0 );
+      }
+
+      tIt = tree.find( idx );
+
+    } while( idx != tIt->second.first );
 
-    } while( idx != tree.find( idx )->second );
 
   } // rof
 
+  /* TODO
+     int i = 1;
+     TImage::PointType ori = input_image->GetOrigin( );
+     std::cout << ori << std::endl;
+     for( TBranches::const_iterator bIt = branches.begin( ); bIt != branches.end( ); ++bIt, ++i )
+     {
+     TImage::PointType::VectorType v1 = bIt->V1 - ori;
+     TImage::PointType::VectorType v2 = bIt->V2 - ori;
+     std::cout
+     << std::setiosflags( std::ios::fixed) << std::setprecision( 3 )
+     << i << "\t1.000\t"
+     << bIt->Length << "\t"
+     << v1[ 0 ] << "\t"
+     << v1[ 1 ] << "\t"
+     << v1[ 2 ] << "\t"
+     << v2[ 0 ] << "\t"
+     << v2[ 1 ] << "\t"
+     << v2[ 2 ] << "\t"
+     << ( v2 - v1 ).GetNorm( )
+     << std::endl;
+
+     } // rof
+  */
+
   vtkSmartPointer< vtkPolyData > vtk_tree =
     vtkSmartPointer< vtkPolyData >::New( );
   vtk_tree->SetPoints( points );
   vtk_tree->SetLines( cells );
   vtk_tree->GetPointData( )->SetScalars( scalars );
 
-  vtkSmartPointer<vtkLookupTable> lut =
-    vtkSmartPointer<vtkLookupTable>::New( );
-  lut->SetNumberOfTableValues( max_mark + 1 );
-  lut->SetTableRange( 0, max_mark );
-  lut->Build( );
-
-  view.AddPolyData( vtk_tree, lut );
+  view.AddPolyData( vtk_tree );
   view.Render( );
   view.Start( );
 
@@ -195,7 +261,7 @@ int main( int argc, char* argv[] )
 
   itk::ImageFileWriter< TDijkstra::TMarkImage >::Pointer marks_writer =
     itk::ImageFileWriter< TDijkstra::TMarkImage >::New( );
-  marks_writer->SetInput( marks );
+  marks_writer->SetInput( new_marks );
   marks_writer->SetFileName( "marks.mhd" );
   marks_writer->Update( );