From 9817c556a0b8b5e3b332d45f07faa84d91afb2d0 Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Mon, 16 Oct 2017 21:25:33 -0500 Subject: [PATCH] ... --- appli/CTBronchi/CTBronchi_process.sh | 26 --- appli/CTBronchi/MoriLabelling.cxx | 2 +- appli/CTBronchi/MoriLabelling.h | 4 +- appli/CTBronchi/Process.cxx | 101 ----------- appli/CTBronchi/Process.sh | 168 ++++++++++++++++++ appli/CTBronchi/SliceBySliceRandomWalker.cxx | 30 ++-- .../RandomWalker/CompareRandomWalkers.cxx | 4 +- examples/image/RandomWalker/Original.cxx | 4 +- lib/fpa/Common/SliceBySliceRandomWalker.h | 19 +- lib/fpa/Common/SliceBySliceRandomWalker.hxx | 43 +++-- 10 files changed, 229 insertions(+), 172 deletions(-) delete mode 100755 appli/CTBronchi/CTBronchi_process.sh delete mode 100644 appli/CTBronchi/Process.cxx create mode 100755 appli/CTBronchi/Process.sh diff --git a/appli/CTBronchi/CTBronchi_process.sh b/appli/CTBronchi/CTBronchi_process.sh deleted file mode 100755 index 1f144db..0000000 --- a/appli/CTBronchi/CTBronchi_process.sh +++ /dev/null @@ -1,26 +0,0 @@ -#!/bin/bash - -input=$1 -dir=`dirname $input` -base=`basename $1 .mhd`_florez-l -output="$dir/$base"_segmentation.mhd -mori="$dir/$base"_mori.mhd -signal="$dir/$base"_mori_signal.txt -labels="$dir/$base"_labels.mhd -diff="$dir/$base"_diff.mhd -seed="$dir"/seed_florez-l.txt - -(>&2 echo -n "Processing $input... ") -echo "************************************************" -./fpa_CTBronchi_Process \ - -in $input \ - -out $output \ - -seed `cat $seed` \ - -out_mori $mori \ - -out_signal $signal \ - -out_labels $labels \ - -out_diff $diff -echo "************************************************" -(>&2 echo "done.") - -## eof - $RCSfile$ diff --git a/appli/CTBronchi/MoriLabelling.cxx b/appli/CTBronchi/MoriLabelling.cxx index b2c901b..70ad18f 100644 --- a/appli/CTBronchi/MoriLabelling.cxx +++ b/appli/CTBronchi/MoriLabelling.cxx @@ -30,7 +30,7 @@ int main( int argc, char* argv[] ) _TStringArg labels( "l", "labels", "Input labels", true, "", "file" ); _TStringArg vesselness( "v", "vesselness", "Input vesselness", true, "", "file" ); _TStringArg out( "o", "output", "Output image", true, "", "file" ); - _TRealArg vThr( "a", "vesselness_thr", "Vesselness threshold", false, 0.05, "value (0.05)" ); + _TRealArg vThr( "a", "vesselness_thr", "Vesselness threshold (%)", false, 5, "value (5)" ); _TPixelArg uThr( "u", "upper_thr", "Upper threshold", false, -400, "value (-400)" ); try { diff --git a/appli/CTBronchi/MoriLabelling.h b/appli/CTBronchi/MoriLabelling.h index ff6553e..f9adb5c 100644 --- a/appli/CTBronchi/MoriLabelling.h +++ b/appli/CTBronchi/MoriLabelling.h @@ -63,7 +63,7 @@ namespace CTBronchi protected: MoriLabelling( ) : Superclass( ), - m_VesselnessThreshold( 0.05 ), + m_VesselnessThreshold( 5 ), m_UpperThreshold( -650 ) { ivqITKInputConfigureMacro( InputLabels, _TLabels ); @@ -91,7 +91,7 @@ namespace CTBronchi minMax->SetImage( this->GetInputVesselness( ) ); minMax->Compute( ); this->m_MinVesselness = - ( 1.0 - this->m_VesselnessThreshold ) * + ( 1.0 - ( this->m_VesselnessThreshold / double( 100 ) ) ) * double( minMax->GetMaximum( ) ); } diff --git a/appli/CTBronchi/Process.cxx b/appli/CTBronchi/Process.cxx deleted file mode 100644 index 19d16d5..0000000 --- a/appli/CTBronchi/Process.cxx +++ /dev/null @@ -1,101 +0,0 @@ -// ========================================================================= -// @author Leonardo Florez Valencia -// @email florez-l@javeriana.edu.co -// ========================================================================= - -#include -#include -#include -#include -#include -#include -#include - -// ------------------------------------------------------------------------- -const unsigned int Dim = 3; -typedef short TInputPixel; -typedef unsigned char TLabelPixel; -typedef itk::Image< TInputPixel, Dim > TImage; -typedef itk::Image< TLabelPixel, Dim > TLabelImage; - -// ------------------------------------------------------------------------- -int main( int argc, char* argv[] ) -{ - std::map< std::string, std::string > args; - try - { - if( CTBronchi::ParseArgs( args, argc, argv ) ) - { - // Parse seed - TImage::PointType seed; - char* str = new char[ args[ "seed" ].size( ) + 1 ]; - strcpy( str, args[ "seed" ].c_str( ) ); - seed[ 0 ] = std::atof( strtok( str, ";" ) ); - seed[ 1 ] = std::atof( strtok( NULL, ";" ) ); - seed[ 2 ] = std::atof( strtok( NULL, ";" ) ); - - // Read input image - TImage::Pointer input_image; - CTBronchi::ReadImage( input_image, args[ "in" ] ); - - // Mori segmentation - TLabelImage::Pointer mori; - TInputPixel opt_thr = CTBronchi::Mori( mori, input_image, seed, args ); - - // Label image - TLabelImage::Pointer labels; - CTBronchi::Label( labels, input_image, mori, opt_thr, args ); - - // Final labels - TLabelImage::Pointer final_labels; - CTBronchi::RandomWalker( final_labels, input_image, labels, args ); - - // Extract label 1 - typedef itk::BinaryThresholdImageFilter< TLabelImage, TLabelImage > TBinThr; - TBinThr::Pointer bin_thr = TBinThr::New( ); - bin_thr->SetInput( final_labels ); - bin_thr->SetLowerThreshold( 1 ); - bin_thr->SetUpperThreshold( 1 ); - bin_thr->SetInsideValue( 1 ); - bin_thr->SetOutsideValue( 0 ); - double t = CTBronchi::MeasureTime( bin_thr ); - std::cout << "Label extracted in " << t << " s" << std::endl; - TLabelImage::Pointer segmentation = bin_thr->GetOutput( ); - - // Save - CTBronchi::WriteImage( segmentation, args[ "out" ] ); - - // Save difference - std::map< std::string, std::string >::const_iterator i = - args.find( "out_diff" ); - if( i != args.end( ) ) - { - typedef itk::AbsoluteValueDifferenceImageFilter< TLabelImage, TLabelImage, TLabelImage > TDiff; - TDiff::Pointer diff = TDiff::New( ); - diff->SetInput1( mori ); - diff->SetInput2( segmentation ); - t = CTBronchi::MeasureTime( diff ); - std::cout << "Difference computed in " << t << " s" << std::endl; - TLabelImage::Pointer difference = diff->GetOutput( ); - CTBronchi::WriteImage( difference, i->second ); - - } // fi - } - else - return( 1 ); - } - catch( std::exception& err ) - { - std::cerr - << "===============================" << std::endl - << "Error caught: " << std::endl - << err.what( ) - << "===============================" << std::endl - << std::endl; - return( 1 ); - - } // yrt - return( 0 ); -} - -// eof - $RCSfile$ diff --git a/appli/CTBronchi/Process.sh b/appli/CTBronchi/Process.sh new file mode 100755 index 0000000..d590b3d --- /dev/null +++ b/appli/CTBronchi/Process.sh @@ -0,0 +1,168 @@ +#!/bin/bash + +## -- Command line options +vesselness_sigma="0.5" +vesselness_alpha1="0.5" +vesselness_alpha2="2" +mori_min_thr="-850" +mori_kernel="20" +mori_signal_thr="100" +mori_signal_influence="0.5" +mori_lower="-1024" +mori_upper="0" +mori_delta="1" +labels_vesselness_thr="5" +labels_upper_thr="-400" +beta="2.5" +epsilon="1e-5" +while [[ "$#" -gt 0 ]]; do + key="$1" + case $key in + -input) + input=`realpath $2` + shift + ;; + -force) + force="1" + ;; + -vesselness_sigma) + vesselness_sigma="$2" + shift + ;; + -vesselness_alpha1) + vesselness_alpha1="$2" + shift + ;; + -vesselness_alpha2) + vesselness_alpha2="$2" + shift + ;; + -mori_min_thr) + mori_min_thr="$2" + shift + ;; + -mori_kernel) + mori_kernel="$2" + shift + ;; + -mori_signal_thr) + mori_signal_thr="$2" + shift + ;; + -mori_signal_influence) + mori_signal_influence="$2" + shift + ;; + -mori_lower) + mori_lower="$2" + shift + ;; + -mori_upper) + mori_upper="$2" + shift + ;; + -mori_delta) + mori_delta="$2" + shift + ;; + -labels_vesselness_thr) + labels_vesselness_thr="$2" + shift + ;; + -labels_upper_thr) + labels_upper_thr="$2" + shift + ;; + -beta) + beta="$2" + shift + ;; + -epsilon) + epsilon="$2" + shift + ;; + -seed) + seed="$2 $3 $4" + shift + shift + shift + ;; + *) + # Do nothing + ;; + esac + shift +done + +## -- Check command line options +if [ -z "$input" ] || [ -z "$seed" ] ; then + (>&2 echo "Usage: $0 -input -seed [-force]") + exit 1 +fi + +base_dir=`dirname $input | xargs realpath` +base_name=$base_dir/`basename $input .mhd` + +vesselness=$base_name"_vesselness.mhd" +mori=$base_name"_mori.mhd" +mori_signal=$base_name"_mori_signal.txt" +labels=$base_name"_labels.mhd" +fastrw=$base_name"_fastrw.mhd" +slicerw=$base_name"_slicerw.mhd" + +(>&2 echo "Processing $input... ") +echo "************************************************" +if [ ! -f $vesselness ] || [ -n "$force" ] ; then + ./fpa_CTBronchi_Vesselness \ + -i $input -o $vesselness \ + -s $vesselness_sigma -a $vesselness_alpha1 -b $vesselness_alpha2 +fi + +if [ ! -f $mori ] || [ -n "$force" ] ; then + ./fpa_CTBronchi_MoriSegmentation \ + -i $input \ + -o $mori \ + -s $mori_signal \ + -p "$seed" \ + -t $mori_min_thr \ + -k $mori_kernel \ + -r $mori_signal_thr \ + -f $mori_signal_influence \ + -l $mori_lower \ + -u $mori_upper \ + -d $mori_delta +fi + +if [ ! -f $labels ] || [ -n "$force" ] ; then + ./fpa_CTBronchi_MoriLabelling \ + -i $input \ + -l $mori \ + -v $vesselness \ + -o $labels \ + -a $labels_vesselness_thr \ + -u $labels_upper_thr +fi + +if [ ! -f $fastrw ] || [ -n "$force" ] ; then + ./fpa_CTBronchi_FastRandomWalker \ + -i $input \ + -l $labels \ + -o $fastrw \ + -b $beta \ + -e $epsilon +fi + +if [ ! -f $slicerw ] || [ -n "$force" ] ; then + ./fpa_CTBronchi_SliceBySliceRandomWalker \ + -i $input \ + -l $mori \ + -v $vesselness \ + -o $slicerw \ + -t $labels_vesselness_thr \ + -b $beta \ + -e $epsilon +fi +echo "************************************************" +(>&2 echo "done.") + +## eof - $RCSfile$ diff --git a/appli/CTBronchi/SliceBySliceRandomWalker.cxx b/appli/CTBronchi/SliceBySliceRandomWalker.cxx index b67e258..69895bf 100644 --- a/appli/CTBronchi/SliceBySliceRandomWalker.cxx +++ b/appli/CTBronchi/SliceBySliceRandomWalker.cxx @@ -31,10 +31,12 @@ int main( int argc, char* argv[] ) _TStringArg out( "o", "output", "Output image", true, "", "file" ); _TRealArg beta( "b", "beta", "Beta", false, 2.5, "value (2.5)" ); _TRealArg eps( "e", "epsilon", "Epsilon", false, 1e-5, "value (1e-5)" ); + _TRealArg vThr( "t", "vthreshold", "Vesselness thresnold (%)", false, 5, "value (5)" ); try { TCLAP::CmdLine cmd( "FastRandomWalker", ' ', "1.0.0" ); cmd.add( eps ); + cmd.add( vThr ); cmd.add( beta ); cmd.add( out ); cmd.add( vesselness ); @@ -74,26 +76,16 @@ int main( int argc, char* argv[] ) filter->SetInput( input_image ); filter->SetInputLabels( input_labels ); filter->SetInputVesselness( input_vesselness ); + filter->SetBeta( beta.getValue( ) ); + filter->SetVesselnessThreshold( vThr.getValue( ) ); + filter->SetEpsilon( eps.getValue( ) ); + double t = CTBronchi::MeasureTime( filter ); + std::cout + << "SliceBySliceRandomWalker executed in " + << t << " s" << std::endl; - // Prepare weight functor - /* TODO - typedef fpa::Functors::Dijkstra::Image::Gaussian< TImage, TScalar > TWeight; - TWeight::Pointer weight = TWeight::New( ); - weight->SetBeta( beta.getValue( ) ); - weight->SetEpsilon( eps.getValue( ) ); - - // Random walk - typedef fpa::Filters::Image::RandomWalker< TImage, TLabels, TScalar > TFilter; - TFilter::Pointer filter = TFilter::New( ); - filter->SetInputImage( input_image ); - filter->SetInputLabels( input_labels ); - filter->SetWeightFunction( weight ); - double t = CTBronchi::MeasureTime( filter ); - std::cout << "FastRandomWalk executed in " << t << " s" << std::endl; - - // Write result - CTBronchi::WriteImage( filter->GetOutputLabels( ), out.getValue( ) ); - */ + // Write result + CTBronchi::WriteImage( filter->GetOutput( ), out.getValue( ) ); } catch( std::exception& err ) { diff --git a/examples/image/RandomWalker/CompareRandomWalkers.cxx b/examples/image/RandomWalker/CompareRandomWalkers.cxx index b3177a9..34f19a0 100644 --- a/examples/image/RandomWalker/CompareRandomWalkers.cxx +++ b/examples/image/RandomWalker/CompareRandomWalkers.cxx @@ -11,7 +11,7 @@ #include #include -#include +#include #include #include @@ -81,7 +81,7 @@ void Original( fpa::Functors::Dijkstra::Image::Gaussian< _TImage, _TScalar > _TFunction; typedef - fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar > + fpa::Common::RandomWalker< _TImage, _TLabels, _TScalar > _TFilter; // Random walker function diff --git a/examples/image/RandomWalker/Original.cxx b/examples/image/RandomWalker/Original.cxx index 4bc9fbe..e68d0f8 100644 --- a/examples/image/RandomWalker/Original.cxx +++ b/examples/image/RandomWalker/Original.cxx @@ -6,7 +6,7 @@ #include #include #include -#include +#include #include // ------------------------------------------------------------------------- @@ -60,7 +60,7 @@ int main( int argc, char* argv[] ) edge->TreatAsWeightOff( ); // Random walker - typedef fpa::Common::OriginalRandomWalker< TImage, TLabels, TScalar > TFilter; + typedef fpa::Common::RandomWalker< TImage, TLabels, TScalar > TFilter; TFilter::Pointer rw = TFilter::New( ); rw->SetInput( input_image_reader->GetOutput( ) ); rw->SetInputLabels( input_labels_reader->GetOutput( ) ); diff --git a/lib/fpa/Common/SliceBySliceRandomWalker.h b/lib/fpa/Common/SliceBySliceRandomWalker.h index f16fe64..ea88fe3 100644 --- a/lib/fpa/Common/SliceBySliceRandomWalker.h +++ b/lib/fpa/Common/SliceBySliceRandomWalker.h @@ -38,6 +38,15 @@ namespace fpa fpa::Common::SliceBySliceRandomWalker, itk::ImageToImageFilter ); + itkGetConstMacro( Epsilon, TScalar ); + itkSetMacro( Epsilon, TScalar ); + + itkGetConstMacro( Beta, TScalar ); + itkSetMacro( Beta, TScalar ); + + itkGetConstMacro( VesselnessThreshold, TScalar ); + itkSetMacro( VesselnessThreshold, TScalar ); + ivqITKInputMacro( InputLabels, TLabels ); ivqITKInputMacro( InputVesselness, TScalarImage ); @@ -50,7 +59,8 @@ namespace fpa void _Composite( typename TScalarImage::Pointer& composite, const TLabels* labels, - const TScalarImage* vesselness + const TScalarImage* vesselness, + const TScalar& maxVess ); template< class _TSlicePtr, class _TInput > @@ -72,6 +82,13 @@ namespace fpa // Purposely not implemented SliceBySliceRandomWalker( const Self& other ); Self& operator=( const Self& other ); + + protected: + TScalar m_Epsilon; + TScalar m_Beta; + TScalar m_VesselnessThreshold; + + TScalar m_VesselnessValue; }; } // ecapseman diff --git a/lib/fpa/Common/SliceBySliceRandomWalker.hxx b/lib/fpa/Common/SliceBySliceRandomWalker.hxx index 1cad2a1..c0bc3ff 100644 --- a/lib/fpa/Common/SliceBySliceRandomWalker.hxx +++ b/lib/fpa/Common/SliceBySliceRandomWalker.hxx @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -18,7 +19,10 @@ template< class _TImage, class _TLabels, class _TScalarImage > fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >:: SliceBySliceRandomWalker( ) - : Superclass( ) + : Superclass( ), + m_Epsilon( TScalar( 1e-5 ) ), + m_Beta( TScalar( 1 ) ), + m_VesselnessThreshold( TScalar( 5 ) ) { ivqITKInputConfigureMacro( InputLabels, TLabels ); ivqITKInputConfigureMacro( InputVesselness, TScalarImage ); @@ -47,18 +51,25 @@ GenerateData( ) const TScalarImage* vesselness = this->GetInputVesselness( ); this->AllocateOutputs( ); + // Get vesselness range + typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax; + typename _TMinMax::Pointer vMinMax = _TMinMax::New( ); + vMinMax->SetImage( vesselness ); + vMinMax->Compute( ); + TScalar vMax = vMinMax->GetMaximum( ); + // Some values typename TImage::RegionType region = input->GetRequestedRegion( ); - typename TImage::SizeType regionIndex = region.GetIndex( ); + typename TImage::IndexType regionIndex = region.GetIndex( ); typename TImage::SizeType sliceSize = region.GetSize( ); int numSlices = sliceSize[ 2 ]; sliceSize[ 2 ] = 0; - this->m_VesselnessValue = - this->m_VesselnessThreshold * this->m_VesselnessMax / TScalar( 100 ); + this->m_VesselnessValue = TScalar( 0.01 ); + this->m_VesselnessValue *= this->m_VesselnessThreshold * vMax; // Composite image typename TScalarImage::Pointer composite; - this->_Composite( composite, labels, vesselness ); + this->_Composite( composite, labels, vesselness, vMax ); // Extract slices std::vector< typename _TSliceImage::Pointer > data3D( numSlices ); @@ -77,7 +88,7 @@ GenerateData( ) typename _TSliceScalarImage::Pointer composite_slice; this->_Slice( input_slice, input, sliceRegion ); this->_Slice( labels_slice, labels, sliceRegion ); - this->_Slice( composite_slice, composite, sliceRegion ); + this->_Slice( composite_slice, composite.GetPointer( ), sliceRegion ); // Update vectors with each image data3D[ i ] = input_slice; @@ -91,7 +102,7 @@ GenerateData( ) this->_RandomWalker( binaryTree, data3D, fusion, false ); // Join results - typedef itk::JoinSeriesImageFilter< _TSliceImage, TImage > _TJoin; + typedef itk::JoinSeriesImageFilter< _TSliceLabels, TLabels > _TJoin; typename _TJoin::Pointer join = _TJoin::New( ); for( int i = 0; i < numSlices; ++i ) join->SetInput( i, binaryTree[ i ] ); @@ -106,15 +117,11 @@ fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >:: _Composite( typename TScalarImage::Pointer& composite, const TLabels* labels, - const TScalarImage* vesselness + const TScalarImage* vesselness, + const TScalar& maxVess ) { - // Min-Max vesselness values - typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax; - typename _TMinMax::Pointer minMax = _TMinMax::New( ); - minMax->SetImage( vesselness ); - minMax->Compute( ); - TScalar maxVess = minMax->GetMaximum( ); + // Some values typename TScalarImage::RegionType region = labels->GetRequestedRegion( ); // Composite image @@ -143,7 +150,7 @@ _Slice( typename _TInput::RegionType region ) { - typedef typename _TSlicePTr::ObjectType _TSlice; + typedef typename _TSlicePtr::ObjectType _TSlice; typedef itk::ExtractImageFilter< _TInput, _TSlice > _TExtract; typename _TExtract::Pointer extract = _TExtract::New( ); @@ -231,9 +238,9 @@ _RandomWalker( // Random walker function typedef fpa::Functors::Dijkstra::Image::Gaussian< _TSliceData3D, _TSliceScalar > _TFunction; - typename TFunction::Pointer edge = TFunction::New( ); + typename _TFunction::Pointer edge = _TFunction::New( ); edge->SetBeta( this->m_Beta ); - edge->SetEpsilon( this->m_Eps ); + edge->SetEpsilon( this->m_Epsilon ); edge->TreatAsWeightOff( ); // Real random walker @@ -245,7 +252,7 @@ _RandomWalker( // Keep maximum typedef itk::MaximumImageFilter< _TSliceBinary > _TMax; - _TMax::Pointer max = _TMax::New(); + typename _TMax::Pointer max = _TMax::New(); max->SetInput( 0, rw->GetOutput( ) ); max->SetInput( 1, next ); max->Update( ); -- 2.47.1