...

diff --git a/appli/CTBronchi/Process.cxx b/appli/CTBronchi/Process.cxx
index 18a150732f3b159279bd03bf0251d3d93bfd63c4..10565399b2d82c83ab8a0b0f8db9937c1563e954 100644 (file)
--- a/appli/CTBronchi/Process.cxx
+++ b/appli/CTBronchi/Process.cxx
@@ -3,6 +3,7 @@
 // =========================================================================
 
 #include <fstream>
+#include <random>
 #include <streambuf>
 #include <tclap/CmdLine.h>
 
@@ -42,6 +43,7 @@ Process( )
   this->m_StrArgs[ "fastrw_points" ] = TStrArg( "", "B", false );
   this->m_StrArgs[ "slicerw_points" ] = TStrArg( "", "C", false );
   this->m_StrArgs[ "andrw_points" ] = TStrArg( "", "D", false );
+  this->m_StrArgs[ "points" ] = TStrArg( "", "E", false );
   this->m_StrArgs[ "seed_file" ] = TStrArg( "", "P", false );
   this->m_StrArgs[ "seed" ] = TStrArg( "", "S", false );
   this->m_StrArgs[ "seed_type" ] = TStrArg( "point", "T", false );
@@ -144,7 +146,7 @@ ParseArguments( int argc, char* argv[] )
   TString bname = std::get< 0 >( this->m_StrArgs[ "input" ] );
   bname = bname.substr( 0, bname.find_last_of( "." ) );
   const unsigned int N = 6;
-  const unsigned int M = 6;
+  const unsigned int M = 7;
   TString names[ N + M ] =
     {
       "vesselness",
@@ -158,7 +160,8 @@ ParseArguments( int argc, char* argv[] )
       "andrw_skeleton",
       "fastrw_points",
       "slicerw_points",
-      "andrw_points"
+      "andrw_points",
+      "points"
     };
   for( unsigned int i = 0; i < N; ++i )
     if( std::get< 0 >( this->m_StrArgs[ names[ i ] ] ) == "" )
@@ -186,20 +189,46 @@ Update( )
     );
   this->_AndImages( this->m_FastRW, this->m_SliceRW, this->m_AndRW );
   this->_Skeleton(
-    this->m_SliceRW, this->m_SliceRWSkeleton,
-    std::get< 0 >( this->m_StrArgs[ "slicerw_skeleton" ] ),
-    std::get< 0 >( this->m_StrArgs[ "slicerw_points" ] )
+    this->m_FastRW, this->m_FastRWSkeleton, this->m_FastRWPoints, 2,
+    std::get< 0 >( this->m_StrArgs[ "fastrw_skeleton" ] )
     );
   this->_Skeleton(
-    this->m_AndRW, this->m_AndRWSkeleton,
-    std::get< 0 >( this->m_StrArgs[ "andrw_skeleton" ] ),
-    std::get< 0 >( this->m_StrArgs[ "andrw_points" ] )
+    this->m_SliceRW, this->m_SliceRWSkeleton, this->m_SliceRWPoints, 1,
+    std::get< 0 >( this->m_StrArgs[ "slicerw_skeleton" ] )
     );
   this->_Skeleton(
-    this->m_FastRW, this->m_FastRWSkeleton,
-    std::get< 0 >( this->m_StrArgs[ "fastrw_skeleton" ] ),
-    std::get< 0 >( this->m_StrArgs[ "fastrw_points" ] )
+    this->m_AndRW, this->m_AndRWSkeleton, this->m_AndRWPoints, 0,
+    std::get< 0 >( this->m_StrArgs[ "andrw_skeleton" ] )
     );
+
+  TEndPoints eval_points;
+  this->_Points(
+    this->m_FastRWPoints, 3,
+    this->m_SliceRWPoints, 3,
+    this->m_AndRWPoints, 4,
+    eval_points
+    );
+
+  if( eval_points.size( ) > 0 )
+  {
+    std::stringstream ePointsStr;
+    for( TEndPoint ep: eval_points )
+    {
+      for( unsigned int d = 0; d < Self::Dim; ++d )
+        ePointsStr << ep.first[ d ] << " ";
+      ePointsStr << ep.second << std::endl;
+
+    } // rof
+
+    std::string pname = std::get< 0 >( this->m_StrArgs[ "points" ] );
+    std::ofstream ePointsF( pname.c_str( ), std::ofstream::binary );
+    if( !ePointsF )
+      throw std::runtime_error(
+        TString( "Unable to write skeleton to \"" ) + pname + "\""
+        );
+    ePointsF.write( ePointsStr.str( ).c_str( ), ePointsStr.str( ).size( ) );
+
+  } // fi
 }
 
 // -------------------------------------------------------------------------
@@ -508,12 +537,12 @@ _AndImages( _TInput& a, _TInput& b, _TInput& c )
 }
 
 // -------------------------------------------------------------------------
-template< class _TInput, class _TSkeleton >
+template< class _TInput, class _TSkeleton, class _TIndices >
 void CTBronchi::Process::
 _Skeleton(
   _TInput& a, _TSkeleton& s,
-  const TString& fname,
-  const TString& pname
+  _TIndices& e, unsigned int id,
+  const TString& fname
   )
 {
   double t = s.Load( fname );
@@ -539,26 +568,55 @@ _Skeleton(
     std::cout << "\"" << fname << "\" saved in " << t1 << " s." << std::endl;
 
     // End points
-    std::stringstream ePointsStr;
-    std::vector< typename _TInput::TImage::IndexType > ePoints =
-      filter.Get( )->GetEndPoints( );
-    for( typename _TInput::TImage::IndexType i: ePoints )
-    {
-      for( unsigned int d = 0; d < _TInput::VDim; ++d )
-        ePointsStr << i[ d ] << " ";
-      ePointsStr << std::endl;
-
-    } // rof
-
-    std::ofstream ePointsF( pname.c_str( ), std::ofstream::binary );
-    if( !ePointsF )
-      throw std::runtime_error(
-        TString( "Unable to write skeleton to \"" ) + pname + "\""
-        );
-    ePointsF.write( ePointsStr.str( ).c_str( ), ePointsStr.str( ).size( ) );
+    std::vector< TIndex > ePoints = filter.Get( )->GetEndPoints( );
+    e.clear( );
+    for( TIndex p: ePoints )
+      e.insert( TEndPoint( p, id ) );
 
   } // fi
   std::cout << "Skeleton computed in " << t << " s." << std::endl;
 }
 
+// -------------------------------------------------------------------------
+template< class _TIndices >
+void CTBronchi::Process::
+_Points(
+  _TIndices& a, unsigned int ca,
+  _TIndices& b, unsigned int cb,
+  _TIndices& c, unsigned int cc,
+  _TIndices& d
+  )
+{
+  std::random_device rd;
+  std::mt19937 gen( rd( ) );
+  d.clear( );
+  while( d.size( ) < ca )
+  {
+    std::uniform_int_distribution< > dis( 0, a.size( ) - 1 );
+    unsigned int N = dis( gen );
+    typename _TIndices::const_iterator it = a.begin( );
+    for( unsigned int i = 0; i < N; ++i, ++it );
+    d.insert( *it );
+      
+  } // fi
+  while( d.size( ) < ca + cb )
+  {
+    std::uniform_int_distribution< > dis( 0, b.size( ) - 1 );
+    unsigned int N = dis( gen );
+    typename _TIndices::const_iterator it = b.begin( );
+    for( unsigned int i = 0; i < N; ++i, ++it );
+    d.insert( *it );
+      
+  } // fi
+  while( d.size( ) < ca + cb + cc )
+  {
+    std::uniform_int_distribution< > dis( 0, c.size( ) - 1 );
+    unsigned int N = dis( gen );
+    typename _TIndices::const_iterator it = c.begin( );
+    for( unsigned int i = 0; i < N; ++i, ++it );
+    d.insert( *it );
+      
+  } // fi
+}
+
 // eof - $RCSfile$

diff --git a/appli/CTBronchi/Process.h b/appli/CTBronchi/Process.h
index 7d5e2dd70121cb8cbaf90d294d2d036fd7472842..b14db2b39c3b2e4a3181d2da1ce5e0402b515f3d 100644 (file)
--- a/appli/CTBronchi/Process.h
+++ b/appli/CTBronchi/Process.h
@@ -5,6 +5,7 @@
 #define __CTBronchi__Process__h__
 
 #include <map>
+#include <set>
 #include <tuple>
 #include <fpa_ctbronchi_export.h>
 #include "Image.h"
@@ -15,6 +16,8 @@ namespace CTBronchi
   class FPA_CTBRONCHI_EXPORT Process
   {
   public:
+    typedef Process Self;
+
     // Some types and values
     static const unsigned int Dim = 3;
     typedef short         TPixel;
@@ -37,9 +40,19 @@ namespace CTBronchi
     typedef CTBronchi::Skeleton< Dim > TSkeleton;
 
     // Seed
-    typedef TPixelImage::TImage::PointType TPoint;
-    typedef TPixelImage::TImage::IndexType TIndex;
-    typedef std::pair< TIndex, TPoint >    TSeed;
+    typedef TPixelImage::TImage::PointType    TPoint;
+    typedef TPixelImage::TImage::IndexType    TIndex;
+    typedef std::pair< TIndex, TPoint >       TSeed;
+    typedef std::pair< TIndex, unsigned int > TEndPoint;
+    struct TEndPointCompare
+    {
+      bool operator()( const TEndPoint& a, const TEndPoint& b ) const
+        {
+          return( cmp( a.first, b.first ) );
+        }
+      TIndex::LexicographicCompare cmp;
+    };
+    typedef std::set< TEndPoint, TEndPointCompare > TEndPoints;
 
   public:
     Process( );
@@ -74,11 +87,19 @@ namespace CTBronchi
     template< class _TInput >
     void _AndImages( _TInput& a, _TInput& b, _TInput& c );
 
-    template< class _TInput, class _TSkeleton >
+    template< class _TInput, class _TSkeleton, class _TIndices >
     void _Skeleton(
       _TInput& a, _TSkeleton& s,
-      const TString& fname,
-      const TString& pname
+      _TIndices& e, unsigned int id,
+      const TString& fname
+      );
+
+    template< class _TIndices >
+    void _Points(
+      _TIndices& a, unsigned int ca,
+      _TIndices& b, unsigned int cb,
+      _TIndices& c, unsigned int cc,
+      _TIndices& d
       );
 
   protected:
@@ -101,6 +122,10 @@ namespace CTBronchi
     TSkeleton m_FastRWSkeleton;
     TSkeleton m_SliceRWSkeleton;
     TSkeleton m_AndRWSkeleton;
+
+    TEndPoints m_FastRWPoints;
+    TEndPoints m_SliceRWPoints;
+    TEndPoints m_AndRWPoints;
   };
 
 } // ecapseman