From 6c5c1fc0b8eefdd093932cc11136e3e1c1c5e4e8 Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Wed, 6 Dec 2017 21:52:50 -0500 Subject: [PATCH] ... --- appli/CTBronchi/MoriLabelling.h | 53 ++++++++++++++++----------------- appli/CTBronchi/Process.cxx | 42 +++++++++++++++++++++----- appli/CTBronchi/Process.h | 6 +++- 3 files changed, 65 insertions(+), 36 deletions(-) diff --git a/appli/CTBronchi/MoriLabelling.h b/appli/CTBronchi/MoriLabelling.h index d6f1790..aa65e97 100644 --- a/appli/CTBronchi/MoriLabelling.h +++ b/appli/CTBronchi/MoriLabelling.h @@ -5,7 +5,7 @@ #ifndef __CTBronchi__MoriLabelling__h__ #define __CTBronchi__MoriLabelling__h__ -#include +#include #include #include #include @@ -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 d1498ea..18a1507 100644 --- 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 d9b939b..7d5e2dd 100644 --- 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; -- 2.45.1