From e69d65c084b579e4014cf35688eae78a9a397834 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Leonardo=20Fl=C3=B3rez-Valencia?= Date: Tue, 11 Jul 2017 17:05:10 -0500 Subject: [PATCH] ... --- appli/CTBronchi/CMakeLists.txt | 5 + appli/CTBronchi/CTBronchi_process.sh | 69 ++++++++++ appli/CTBronchi/MoriLabelling.h | 97 ++++++-------- appli/CTBronchi/MoriLabelling.hxx | 186 ++++++++++++++++++++++++++- appli/CTBronchi/RandomWalker.cxx | 40 +++++- 5 files changed, 325 insertions(+), 72 deletions(-) create mode 100755 appli/CTBronchi/CTBronchi_process.sh diff --git a/appli/CTBronchi/CMakeLists.txt b/appli/CTBronchi/CMakeLists.txt index 1d7e316..5b20064 100644 --- a/appli/CTBronchi/CMakeLists.txt +++ b/appli/CTBronchi/CMakeLists.txt @@ -1,5 +1,10 @@ option(fpa_BUILD_CTBronchi "Build bronchi analysis from CT images applications?" OFF) if(fpa_BUILD_CTBronchi) + configure_file( + CTBronchi_process.sh + ${PROJECT_BINARY_DIR}/CTBronchi_process.sh + COPYONLY + ) include_directories( ${PROJECT_SOURCE_DIR}/lib ${PROJECT_BINARY_DIR}/lib diff --git a/appli/CTBronchi/CTBronchi_process.sh b/appli/CTBronchi/CTBronchi_process.sh new file mode 100755 index 0000000..15bd591 --- /dev/null +++ b/appli/CTBronchi/CTBronchi_process.sh @@ -0,0 +1,69 @@ +#!/bin/bash + +## Command line arguments +if [ "$#" -lt 5 ]; then + echo "Usage: $0 input_image [index/point] seed_x seed_y seed_z" + exit 1 +fi + +exec_dir=`dirname $0` +mori_seg=`dirname $0`/fpa_CTBronchi_MoriSegmentation +mori_lab=`dirname $0`/fpa_CTBronchi_MoriLabelling +random_walker=`dirname $0`/fpa_CTBronchi_RandomWalker +input_image=$1 +seed_type=$2 +seed_x=$3 +seed_y=$4 +seed_z=$5 + +base_name=`dirname $input_image`/`basename $input_image .mhd` + +mori_output_image="$base_name"_mori.mhd +mori_output_signal="$base_name"_mori_signal.txt +mori_init_threshold=-1024 +mori_end_threshold=0 +mori_delta=1 +mori_minimum_threshold=-850 +mori_inside_value=255 +mori_outside_value=0 +mori_signal_kernel_size=20 +mori_signal_threshold=500 +mori_signal_influence=500 + +labels_output_image="$base_name"_labels.mhd +label_upper_threshold=-600 +label_inside=1 +label_outside=2 + +random_walker_output_image="$base_name"_rw.mhd +random_walker_alpha=0 +random_walker_beta=20 + +$mori_seg \ + $input_image $mori_output_image $mori_output_signal \ + $mori_init_threshold \ + $mori_end_threshold \ + $mori_delta \ + $mori_minimum_threshold \ + $mori_inside_value \ + $mori_outside_value \ + $mori_signal_kernel_size \ + $mori_signal_threshold \ + $mori_signal_influence \ + $seed_type \ + $seed_x $seed_y $seed_z + +$mori_lab \ + $input_image $mori_output_image $labels_output_image \ + $label_upper_threshold \ + $mori_inside_value \ + $label_inside \ + $label_outside + +$random_walker \ + $input_image $labels_output_image $random_walker_output_image \ + $label_inside \ + $random_walker_alpha \ + $random_walker_beta + +## eof - $RCSfile$ diff --git a/appli/CTBronchi/MoriLabelling.h b/appli/CTBronchi/MoriLabelling.h index f1ac2ea..b51cfeb 100644 --- a/appli/CTBronchi/MoriLabelling.h +++ b/appli/CTBronchi/MoriLabelling.h @@ -6,12 +6,11 @@ #ifndef __CTBronchi__MoriLabelling__h__ #define __CTBronchi__MoriLabelling__h__ -#include -#include -#include - -#include - +#include +#include +#include +#include +#include namespace CTBronchi { @@ -19,35 +18,40 @@ namespace CTBronchi */ template< class _TInputImage, class _TLabelImage > class MoriLabelling - : public itk::ImageToImageFilter< _TLabelImage, _TLabelImage > + : public fpa::Base::RegionGrow< fpa::Image::Algorithm< _TInputImage, _TLabelImage, fpa::Base::MarksInterface< typename _TInputImage::IndexType >, fpa::Image::LabelledSeedsInterface< typename _TInputImage::IndexType, typename _TInputImage::PointType, typename _TInputImage::PixelType, typename _TLabelImage::PixelType, typename _TLabelImage::PixelType, typename _TInputImage::IndexType::LexicographicCompare > > > { public: - typedef MoriLabelling Self; - typedef itk::ImageToImageFilter< _TLabelImage, _TLabelImage > Superclass; - typedef itk::SmartPointer< Self > Pointer; - typedef itk::SmartPointer< const Self > ConstPointer; - typedef _TInputImage TInputImage; typedef _TLabelImage TLabelImage; - typedef typename TInputImage::PixelType TPixel; - typedef typename TLabelImage::PixelType TLabel; - typedef typename TLabelImage::RegionType TRegion; + typedef typename TInputImage::PixelType TInputValue; + typedef typename TInputImage::PointType TPoint; + typedef typename TInputImage::IndexType TVertex; + typedef typename TLabelImage::PixelType TOutputValue; + typedef typename TVertex::LexicographicCompare TVertexCompare; - public: - itkNewMacro( Self ); - itkTypeMacro( MoriLabelling, itk::ImageToImageFilter ); + typedef fpa::Base::MarksInterface< TVertex > TMarksInterface; + typedef fpa::Image::LabelledSeedsInterface< TVertex, TPoint, TInputValue, TOutputValue, TOutputValue, TVertexCompare > TSeedsInterface; + typedef fpa::Image::Algorithm< TInputImage, TLabelImage, TMarksInterface, TSeedsInterface > TAlgorithm; - itkGetConstMacro( UpperThreshold, TPixel ); - itkSetMacro( UpperThreshold, TPixel ); + typedef MoriLabelling Self; + typedef fpa::Base::RegionGrow< TAlgorithm > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; - itkGetConstMacro( InsideValue, TLabel ); - itkSetMacro( InsideValue, TLabel ); + typedef typename TSeedsInterface::TNode TNode; + typedef typename TSeedsInterface::TNodes TNodes; - itkGetConstMacro( InsideLabel, TLabel ); - itkSetMacro( InsideLabel, TLabel ); + typedef fpa::Base::Functors::RegionGrow::BinaryThreshold< TInputValue > TThresholdFunction; - itkGetConstMacro( OutsideLabel, TLabel ); - itkSetMacro( OutsideLabel, TLabel ); + public: + itkNewMacro( Self ); + itkTypeMacro( MoriLabelling, fpa::Base::RegionGrow ); + + itkGetConstMacro( InsideLabel, TOutputValue ); + itkSetMacro( InsideLabel, TOutputValue ); + + itkGetConstMacro( OutsideLabel, TOutputValue ); + itkSetMacro( OutsideLabel, TOutputValue ); public: const TLabelImage* GetInputLabelImage( ) const; @@ -56,40 +60,15 @@ namespace CTBronchi const TInputImage* GetInputRawImage( ) const; void SetInputRawImage( TInputImage* image ); + TInputValue GetUpperThreshold( ) const; + void SetUpperThreshold( TInputValue t ); + protected: MoriLabelling( ); virtual ~MoriLabelling( ); - virtual void BeforeThreadedGenerateData( ) override - { - this->Superclass::BeforeThreadedGenerateData( ); - - const TInputImage* raw = this->GetInputRawImage( ); - typename TInputImage::SpacingType spac = raw->GetSpacing( ); - double sigma = spac[ 0 ]; - for( unsigned int d = 1; d < TInputImage::ImageDimension; ++d ) - sigma = ( spac[ d ] < sigma )? spac[ d ]: sigma; - sigma *= 1.5; - - typedef itk::HessianRecursiveGaussianImageFilter< TInputImage > _THessian; - typename _THessian::Pointer hessian = _THessian::New( ); - hessian->SetInput( raw ); - hessian->SetSigma( sigma ); - - typedef itk::Hessian3DToVesselnessMeasureImageFilter< double > _TVesselness; - typename _TVesselness::Pointer vesselness = _TVesselness::New( ); - vesselness->SetInput( hessian->GetOutput( ) ); - vesselness->SetAlpha1( 0.5 ); - vesselness->SetAlpha2( 2.0 ); - - typename itk::ImageFileWriter< typename _TVesselness::OutputImageType >::Pointer w = - itk::ImageFileWriter< typename _TVesselness::OutputImageType >::New( ); - w->SetInput( vesselness->GetOutput( ) ); - w->SetFileName( "vessel.mhd" ); - w->Update( ); - } - - virtual void ThreadedGenerateData( const TRegion& region, itk::ThreadIdType threadId ) override; + virtual TNodes _UnifySeeds( ) override; + virtual void _UpdateOutputValue( TNode& n ) override; private: // Purposely not implemented. @@ -97,10 +76,8 @@ namespace CTBronchi Self& operator=( const Self& other ); protected: - TPixel m_UpperThreshold; - TLabel m_InsideValue; - TLabel m_InsideLabel; - TLabel m_OutsideLabel; + TOutputValue m_InsideLabel; + TOutputValue m_OutsideLabel; }; } // ecapseman diff --git a/appli/CTBronchi/MoriLabelling.hxx b/appli/CTBronchi/MoriLabelling.hxx index 5ca4fc4..bc38389 100644 --- a/appli/CTBronchi/MoriLabelling.hxx +++ b/appli/CTBronchi/MoriLabelling.hxx @@ -6,6 +6,181 @@ #ifndef __CTBronchi__MoriLabelling__hxx__ #define __CTBronchi__MoriLabelling__hxx__ +#include + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TLabelImage > +const typename CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: +TLabelImage* CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: +GetInputLabelImage( ) const +{ + return( this->GetLabels( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TLabelImage > +void CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: +SetInputLabelImage( TLabelImage* image ) +{ + this->SetLabels( image ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TLabelImage > +const typename CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: +TInputImage* CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: +GetInputRawImage( ) const +{ + return( this->GetInput( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TLabelImage > +void CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: +SetInputRawImage( TInputImage* image ) +{ + this->SetInput( image ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TLabelImage > +typename CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: +TInputValue CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: +GetUpperThreshold( ) const +{ + const TThresholdFunction* func = + dynamic_cast< const TThresholdFunction* >( this->GetValuePredicate( ) ); + if( func != NULL ) + return( func->GetUpper( ) ); + else + return( TInputValue( 0 ) ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TLabelImage > +void CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: +SetUpperThreshold( TInputValue t ) +{ + TThresholdFunction* func = + dynamic_cast< TThresholdFunction* >( this->GetValuePredicate( ) ); + if( func != NULL ) + func->SetUpper( t ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TLabelImage > +CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: +MoriLabelling( ) + : Superclass( ), + m_InsideLabel( TOutputValue( 0 ) ), + m_OutsideLabel( TOutputValue( 0 ) ) +{ + typename TThresholdFunction::Pointer func = TThresholdFunction::New( ); + this->SetPredicate( func ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TLabelImage > +CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: +~MoriLabelling( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TLabelImage > +typename CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: +TNodes CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: +_UnifySeeds( ) +{ + this->m_Seeds.clear( ); + const TLabelImage* lbl = this->GetLabels( ); + if( lbl == NULL ) + { + std::ostringstream msg; + msg << "itk::ERROR: CTBronchi::MoriLabelling (" << this + << "): Labelled image not defined."; + ::itk::ExceptionObject e( + __FILE__, __LINE__, msg.str( ).c_str( ), ITK_LOCATION + ); + throw e; + + } // fi + + // Iterate over labels + typename TLabelImage::RegionType reg = lbl->GetRequestedRegion( ); + itk::ImageRegionConstIteratorWithIndex< TLabelImage > lIt( lbl, reg ); + for( lIt.GoToBegin( ); !lIt.IsAtEnd( ); ++lIt ) + { + if( lIt.Get( ) > 0 ) + { + bool is_seed = false; + for( unsigned int d = 0; d < TLabelImage::ImageDimension; ++d ) + { + for( int s = -1; s <= 1; s += 2 ) + { + TVertex neigh = lIt.GetIndex( ); + neigh[ d ] += s; + if( reg.IsInside( neigh ) ) + is_seed |= ( lbl->GetPixel( neigh ) == 0 ); + + } // rof + + } // rof + + if( !is_seed ) + { + typename TSeedsInterface::TNode node; + node.Vertex = lIt.GetIndex( ); + node.Parent = lIt.GetIndex( ); + node.FrontId = lIt.Get( ); + node.Value = this->m_InsideLabel; + this->_Mark( node.Vertex, node.FrontId ); + this->_UpdateOutputValue( node ); + } + else + { + typename TSeedsInterface::TSeed seed; + seed.Vertex = lIt.GetIndex( ); + seed.IsPoint = false; + seed.FrontId = lIt.Get( ); + this->m_Seeds.push_back( seed ); + + } // fi + + } // fi + + } // rof + + // Ok, finish initialization + return( this->Superclass::_UnifySeeds( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TLabelImage > +void CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: +_UpdateOutputValue( TNode& n ) +{ + const TInputImage* input = this->GetInputRawImage( ); + const TLabelImage* labels = this->GetInputLabelImage( ); +} + + +/* TODO + TOutputValue m_InsideValue; + TOutputValue m_InsideLabel; + TOutputValue m_OutsideLabel; +*/ + + + + + + + + + + +/* TODO #include #include @@ -52,10 +227,10 @@ template< class _TInputImage, class _TLabelImage > CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: MoriLabelling( ) : Superclass( ), - m_UpperThreshold( TLabel( 0 ) ), - m_InsideValue( TLabel( 0 ) ), - m_InsideLabel( TLabel( 0 ) ), - m_OutsideLabel( TLabel( 0 ) ) + m_UpperThreshold( TOutputValue( 0 ) ), + m_InsideValue( TOutputValue( 0 ) ), + m_InsideLabel( TOutputValue( 0 ) ), + m_OutsideLabel( TOutputValue( 0 ) ) { this->SetNumberOfRequiredInputs( 2 ); } @@ -89,7 +264,7 @@ ThreadedGenerateData( const TRegion& region, itk::ThreadIdType threadId ) if( this->m_UpperThreshold < iIt.Get( ) ) oIt.Set( this->m_OutsideLabel ); else - oIt.Set( TLabel( 0 ) ); + oIt.Set( TOutputValue( 0 ) ); } else oIt.Set( this->m_InsideLabel ); @@ -99,6 +274,7 @@ ThreadedGenerateData( const TRegion& region, itk::ThreadIdType threadId ) } // elihw } +*/ #endif // __CTBronchi__MoriLabelling__hxx__ diff --git a/appli/CTBronchi/RandomWalker.cxx b/appli/CTBronchi/RandomWalker.cxx index f6e81fb..5704558 100644 --- a/appli/CTBronchi/RandomWalker.cxx +++ b/appli/CTBronchi/RandomWalker.cxx @@ -5,29 +5,36 @@ #include #include +#include +#include // ------------------------------------------------------------------------- const unsigned int Dim = 3; typedef short TPixel; typedef unsigned short TLabel; +typedef unsigned char TBinary; typedef double TScalar; typedef itk::Image< TPixel, Dim > TInputImage; typedef itk::Image< TLabel, Dim > TLabelImage; +typedef itk::Image< TBinary, Dim > TBinaryImage; typedef fpa::Image::RandomWalker< TInputImage, TLabelImage, TScalar > TFilter; typedef fpa::Image::Functors::Dijkstra::Gaussian< TInputImage, TScalar > TWeight; +typedef ivq::ITK::ExtractLabelFunction< TLabel, TBinary > TLabelFunction; +typedef ivq::ITK::ImageUnaryFunctionFilter< TLabelImage, TBinaryImage > TLabelExtractor; + // ------------------------------------------------------------------------- int main( int argc, char* argv[] ) { // Get arguments - if( argc < 6 ) + if( argc < 7 ) { std::cerr << "Usage: " << argv[ 0 ] << std::endl << " input_image label_image output_image" << std::endl - << " alpha(0) beta(100)" + << " label alpha(0) beta(100)" << std::endl; return( 1 ); @@ -35,8 +42,9 @@ int main( int argc, char* argv[] ) std::string input_image_filename = argv[ 1 ]; std::string label_image_filename = argv[ 2 ]; std::string output_image_filename = argv[ 3 ]; - double alpha = std::atof( argv[ 4 ] ); - double beta = std::atof( argv[ 5 ] ); + TLabel label = TLabel( std::atoi( argv[ 4 ] ) ); + double alpha = std::atof( argv[ 5 ] ); + double beta = std::atof( argv[ 6 ] ); // Read images itk::ImageFileReader< TInputImage >::Pointer input_image_reader = @@ -58,6 +66,17 @@ int main( int argc, char* argv[] ) filter->SetLabels( label_image_reader->GetOutput( ) ); filter->SetWeightFunction( weight ); + // Label extraction + TLabelFunction::Pointer label_function = TLabelFunction::New( ); + label_function->SetLabel( label ); + label_function->SetInsideValue( 255 ); + label_function->SetOutsideValue( 0 ); + + TLabelExtractor::Pointer label_extractor = TLabelExtractor::New( ); + label_extractor->SetInput( filter->GetMarks( ) ); + label_extractor->SetFunction( label_function ); + + // Show some information std::cout << "----------------------------------------------" << std::endl; std::cout << "Image: " << input_image_filename << std::endl; @@ -85,6 +104,13 @@ int main( int argc, char* argv[] ) tr = te - ts; std::cout << "Labelling time: " << tr.count( ) << " s" << std::endl; + + ts = std::chrono::high_resolution_clock::now( ); + label_extractor->Update( ); + te = std::chrono::high_resolution_clock::now( ); + tr = te - ts; + std::cout + << "Extract label time: " << tr.count( ) << " s" << std::endl; } catch( std::exception& err ) { @@ -94,9 +120,9 @@ int main( int argc, char* argv[] ) } // yrt // Save output image - itk::ImageFileWriter< TLabelImage >::Pointer output_image_writer = - itk::ImageFileWriter< TLabelImage >::New( ); - output_image_writer->SetInput( filter->GetMarks( ) ); + itk::ImageFileWriter< TBinaryImage >::Pointer output_image_writer = + itk::ImageFileWriter< TBinaryImage >::New( ); + output_image_writer->SetInput( label_extractor->GetOutput( ) ); output_image_writer->SetFileName( output_image_filename ); try { -- 2.45.0