...

[FrontAlgorithms.git] / appli / CTBronchi / Process.cxx

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$