...

diff --git a/appli/CTBronchi/MoriLabelling.h b/appli/CTBronchi/MoriLabelling.h
index d6f1790d7e21499fbf8c1a8785cfb1768422da6b..aa65e97daef5a9e8914357e5b48887ff88d7dab6 100644 (file)
--- a/appli/CTBronchi/MoriLabelling.h
+++ b/appli/CTBronchi/MoriLabelling.h
@@ -5,7 +5,7 @@
 #ifndef __CTBronchi__MoriLabelling__h__
 #define __CTBronchi__MoriLabelling__h__
 
-#include <itkStatisticsImageFilter.h>
+#include <itkLabelStatisticsImageFilter.h>
 #include <fpa/Filters/BaseMarksInterface.h>
 #include <fpa/Filters/Image/SeedsFromLabelsInterface.h>
 #include <fpa/Filters/Image/DefaultTraits.h>
@@ -86,35 +86,34 @@ namespace CTBronchi
 
         this->m_Functor->SetUpperThreshold( this->m_UpperThreshold );
 
-        typedef itk::StatisticsImageFilter< _TScalarImage > _TStats;
+        typedef itk::LabelStatisticsImageFilter< _TScalarImage, _TLabels > _TStats;
         typename _TStats::Pointer stats = _TStats::New( );
         stats->SetInput( this->GetInputVesselness( ) );
+        stats->SetLabelInput( this->GetInputLabels( ) );
         stats->Update( );
-        double vAvg = double( stats->GetMean( ) );
-        double vStd = double( stats->GetSigma( ) );
-        double vMin = double( stats->GetMinimum( ) );
-        double vMax = double( stats->GetMaximum( ) );
-
-        this->m_MinVesselness = vAvg + ( vStd * double( 5 ) );
-        /* TODO
-           std::cout
-           << vAvg << std::endl
-           << vStd << std::endl
-           << vMin << std::endl
-           << vMax << std::endl;
-
-           std::cout
-           << ( this->m_VesselnessThreshold / double( 100 ) ) * double( vMax )
-           << std::endl
-           << vAvg + ( vStd * double( 3 ) ) << std::endl;
-
-           std::exit( 1 );
-        */
-        /* TODO
-           this->m_MinVesselness =
-           ( this->m_VesselnessThreshold / double( 100 ) ) *
-           double( stats->GetMaximum( ) );
-        */
+
+        typedef typename _TStats::ValidLabelValuesContainerType _TValidLbl;
+        typedef typename _TStats::LabelPixelType                _TLabelPx;
+        unsigned long minCount = std::numeric_limits< unsigned long >::max( );
+        double vAvg = 0, vStd = 0;
+        typename _TValidLbl::const_iterator vIt = stats->GetValidLabelValues( ).begin( );
+        for( ; vIt != stats->GetValidLabelValues( ).end( ); ++vIt )
+        {
+          if( stats->HasLabel( *vIt ) )
+          {
+            unsigned long c = stats->GetCount( *vIt );
+            if( c < minCount )
+            {
+              minCount = c;
+              vAvg = stats->GetMean( *vIt );
+              vStd = stats->GetSigma( *vIt );
+
+            } // fi
+ 
+          } // fi
+
+        } // rof
+        this->m_MinVesselness = vAvg + ( vStd * this->m_VesselnessThreshold );
       }
 
     virtual void _PostComputeOutputValue( TNode& n ) override

diff --git a/appli/CTBronchi/Process.cxx b/appli/CTBronchi/Process.cxx
index d1498ea2f7ca62fe146cfca09cbc571590fe40a5..18a150732f3b159279bd03bf0251d3d93bfd63c4 100644 (file)
--- a/appli/CTBronchi/Process.cxx
+++ b/appli/CTBronchi/Process.cxx
@@ -59,7 +59,7 @@ Process( )
   this->m_DblArgs[ "mori_upper" ] = TDblArg( 0, "l", false );
   this->m_DblArgs[ "mori_delta" ] = TDblArg( 1, "m", false );
   this->m_DblArgs[ "slicerw_thr" ] = TDblArg( 5, "n", false );
-  this->m_DblArgs[ "fastrw_thr" ] = TDblArg( 65, "o", false );
+  this->m_DblArgs[ "fastrw_thr" ] = TDblArg( 5, "o", false );
   this->m_DblArgs[ "labels_upper_thr" ] = TDblArg( -400, "p", false );
 }
 
@@ -187,15 +187,18 @@ 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_skeleton" ] ),
+    std::get< 0 >( this->m_StrArgs[ "slicerw_points" ] )
     );
   this->_Skeleton(
     this->m_AndRW, this->m_AndRWSkeleton,
-    std::get< 0 >( this->m_StrArgs[ "andrw_skeleton" ] )
+    std::get< 0 >( this->m_StrArgs[ "andrw_skeleton" ] ),
+    std::get< 0 >( this->m_StrArgs[ "andrw_points" ] )
     );
   this->_Skeleton(
     this->m_FastRW, this->m_FastRWSkeleton,
-    std::get< 0 >( this->m_StrArgs[ "fastrw_skeleton" ] )
+    std::get< 0 >( this->m_StrArgs[ "fastrw_skeleton" ] ),
+    std::get< 0 >( this->m_StrArgs[ "fastrw_points" ] )
     );
 }
 
@@ -285,7 +288,7 @@ _Vesselness( _TInput& input, _TVesselness& vesselness )
     std::cout << "Hessian executed in " << t1 << " s" << std::endl;
 
     // Vesselness
-    typedef CTBronchi::Filter< itk::Hessian3DToVesselnessMeasureImageFilter< typename _TVesselness::TImage::PixelType > > _TVesselnessFilter;
+    typedef CTBronchi::Filter< itk::Hessian3DToVesselnessMeasureImageFilter< typename _TVesselness::TPixel > > _TVesselnessFilter;
     _TVesselnessFilter vFilter;
     vFilter.Get( )->SetInput( hessian.Get( )->GetOutput( ) );
     vFilter.Get( )->SetAlpha1(
@@ -404,13 +407,13 @@ _FastRW( _TInput& input, _TLabels& labels, _TFastRW& fastrw )
     t = 0;
 
     // Prepare weight functor
-    typedef fpa::Functors::Dijkstra::Image::Gaussian< typename _TInput::TImage, typename _TFastRW::TImage::PixelType > TWeight;
+    typedef fpa::Functors::Dijkstra::Image::Gaussian< typename _TInput::TImage, typename _TFastRW::TPixel > TWeight;
     typename TWeight::Pointer weight = TWeight::New( );
     weight->SetBeta( std::get< 0 >( this->m_DblArgs[ "beta" ] ) );
     weight->SetEpsilon( std::get< 0 >( this->m_DblArgs[ "epsilon" ] ) );
 
     // Random walk
-    typedef CTBronchi::Filter< fpa::Filters::Image::RandomWalker< typename _TInput::TImage, typename _TLabels::TImage, typename _TFastRW::TImage::PixelType > > TFilter;
+    typedef CTBronchi::Filter< fpa::Filters::Image::RandomWalker< typename _TInput::TImage, typename _TLabels::TImage, typename _TFastRW::TPixel > > TFilter;
     TFilter filter;
     filter.Get( )->SetInputImage( input.Get( ) );
     filter.Get( )->SetInputLabels( labels.Get( ) );
@@ -507,7 +510,11 @@ _AndImages( _TInput& a, _TInput& b, _TInput& c )
 // -------------------------------------------------------------------------
 template< class _TInput, class _TSkeleton >
 void CTBronchi::Process::
-_Skeleton( _TInput& a, _TSkeleton& s, const TString& fname )
+_Skeleton(
+  _TInput& a, _TSkeleton& s,
+  const TString& fname,
+  const TString& pname
+  )
 {
   double t = s.Load( fname );
   if( t < 0 )
@@ -531,6 +538,25 @@ _Skeleton( _TInput& a, _TSkeleton& s, const TString& fname )
     t += t1;
     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( ) );
+
   } // fi
   std::cout << "Skeleton computed in " << t << " s." << std::endl;
 }

diff --git a/appli/CTBronchi/Process.h b/appli/CTBronchi/Process.h
index d9b939b2c0da79b8902da8256f6619ab29a27019..7d5e2dd70121cb8cbaf90d294d2d036fd7472842 100644 (file)
--- a/appli/CTBronchi/Process.h
+++ b/appli/CTBronchi/Process.h
@@ -75,7 +75,11 @@ namespace CTBronchi
     void _AndImages( _TInput& a, _TInput& b, _TInput& c );
 
     template< class _TInput, class _TSkeleton >
-    void _Skeleton( _TInput& a, _TSkeleton& s, const TString& fname );
+    void _Skeleton(
+      _TInput& a, _TSkeleton& s,
+      const TString& fname,
+      const TString& pname
+      );
 
   protected:
     TStrArgs m_StrArgs;