From ae58efdf42a219a951e60338f30c40e995befea7 Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Sun, 26 Nov 2017 22:02:11 -0500 Subject: [PATCH] ... --- appli/CTBronchi/CMakeLists.txt | 1 + appli/CTBronchi/Process.sh | 26 +++-- appli/CTBronchi/Skeleton.cxx | 119 ++++++++++++++++++++++ lib/fpa/Filters/Image/ExtractSkeleton.h | 3 + lib/fpa/Filters/Image/ExtractSkeleton.hxx | 1 + 5 files changed, 143 insertions(+), 7 deletions(-) create mode 100644 appli/CTBronchi/Skeleton.cxx diff --git a/appli/CTBronchi/CMakeLists.txt b/appli/CTBronchi/CMakeLists.txt index 5f10a8e..b5f19be 100644 --- a/appli/CTBronchi/CMakeLists.txt +++ b/appli/CTBronchi/CMakeLists.txt @@ -13,6 +13,7 @@ if(fpa_BUILD_CTBronchi) FastRandomWalker SliceBySliceRandomWalker AndSegmentations + Skeleton ) foreach(_e ${_examples}) BuildApplication( diff --git a/appli/CTBronchi/Process.sh b/appli/CTBronchi/Process.sh index 618a3bc..5f42202 100755 --- a/appli/CTBronchi/Process.sh +++ b/appli/CTBronchi/Process.sh @@ -2,6 +2,7 @@ ## -- Command line options curr_dir=`dirname $0` +ext="mhd" vesselness_sigma="0.5" vesselness_alpha1="0.5" vesselness_alpha2="2" @@ -13,7 +14,7 @@ mori_lower="-1024" mori_upper="0" mori_delta="1" slice_by_slice_vesselness_thr="5" -labels_vesselness_thr="65" +fast_vesselness_thr="65" labels_upper_thr="-400" beta="2.5" epsilon="1e-5" @@ -111,14 +112,17 @@ if [ -z "$input" ] || [ -z "$seed" ] ; then fi base_dir=`dirname $input | xargs realpath` -base_name=$base_dir/`basename $input .mha` +base_name=$base_dir/`basename $input .$ext` -vesselness=$base_name"_vesselness.mha" -mori=$base_name"_mori.mha" +vesselness=$base_name"_vesselness.$ext" +mori=$base_name"_mori.$ext" mori_signal=$base_name"_mori_signal.txt" -labels=$base_name"_labels.mha" -fastrw=$base_name"_fastrw.mha" -slicerw=$base_name"_slicerw.mha" +labels=$base_name"_labels.$ext" +fastrw=$base_name"_fastrw.$ext" +slicerw=$base_name"_slicerw.$ext" +skeleton_fastrw=$base_name"_fastrw_skeleton.txt" +skeleton_slicerw=$base_name"_slicerw_skeleton.txt" +points=$base_name"_points.txt" echo "************************************************" (>&2 echo "Processing $input... ") @@ -173,6 +177,14 @@ if [ ! -f $slicerw ] || [ -n "$force" ] ; then -b $beta \ -e $epsilon fi + +if [ ! -f $skeleton_slicerw ] || [ -n "$force" ] ; then + $curr_dir/fpa_CTBronchi_Skeleton \ + -i $slicerw \ + -o $skeleton_slicerw \ + -e $points \ + -p "$seed" +fi (>&2 echo "done.") echo "done." echo "************************************************" diff --git a/appli/CTBronchi/Skeleton.cxx b/appli/CTBronchi/Skeleton.cxx new file mode 100644 index 0000000..4110825 --- /dev/null +++ b/appli/CTBronchi/Skeleton.cxx @@ -0,0 +1,119 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#include +#include +#include +#include +#include +#include +#include "Functions.h" + +// ------------------------------------------------------------------------- +const unsigned int Dim = 3; +typedef unsigned char TPixel; +typedef itk::Image< TPixel, Dim > TImage; + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + typedef TCLAP::ValueArg< std::string > _TStringArg; + typedef TCLAP::ValueArg< double > _TRealArg; + typedef TCLAP::ValueArg< unsigned long > _TUIntArg; + typedef TCLAP::ValueArg< TPixel > _TPixelArg; + + // Parse input line + _TStringArg in( "i", "input", "Input image", true, "", "file" ); + _TStringArg out( "o", "output", "Output skeleton", true, "", "file" ); + _TStringArg oute( "e", "end_points", "Output end points", false, "", "file" ); + _TStringArg aSeed( "p", "seed", "Input seed", true, "", "\"x y z\"" ); + try + { + TCLAP::CmdLine cmd( "Skeleton", ' ', "1.0.0" ); + cmd.add( aSeed ); + cmd.add( oute ); + cmd.add( out ); + cmd.add( in ); + cmd.parse( argc, argv ); + } + catch( TCLAP::ArgException& err ) + { + std::cerr + << "===============================" << std::endl + << "Error caught: " << std::endl + << err.error( ) << " " << err.argId( ) + << "===============================" << std::endl + << std::endl; + return( 1 ); + + } // yrt + + // Get seed + std::istringstream sSeed( aSeed.getValue( ) ); + TImage::PointType pnt; + sSeed >> pnt[ 0 ] >> pnt[ 1 ] >> pnt[ 2 ]; + + try + { + // Read image + TImage::Pointer input_image; + CTBronchi::ReadImage( input_image, in.getValue( ) ); + + TImage::IndexType seed; + input_image->TransformPhysicalPointToIndex( pnt, seed ); + + // Skeleton extraction + typedef fpa::Filters::Image::ExtractSkeleton< TImage > TFilter; + TFilter::Pointer filter = TFilter::New( ); + filter->SetInput( input_image ); + filter->SeedFromMaximumDistanceOff( ); + filter->SetSeed( seed ); + filter->GetDistanceMap( )->InsideIsPositiveOn( ); + filter->GetDistanceMap( )->SquaredDistanceOff( ); + filter->GetDistanceMap( )->UseImageSpacingOn( ); + double t = CTBronchi::MeasureTime( filter ); + std::cout << "Skeleton executed in " << t << " s" << std::endl; + + // Save results + typedef ivq::ITK::ImageSkeletonWriter< TFilter::TSkeleton > TWriter; + TWriter::Pointer skeleton_writer = TWriter::New( ); + skeleton_writer->SetInput( filter->GetOutput( ) ); + skeleton_writer->SetFileName( out.getValue( ) ); + skeleton_writer->Update( ); + + // Write end points + if( oute.getValue( ) != "" ) + { + std::stringstream eStr; + auto pnts = filter->GetEndPoints( ); + for( unsigned int i = 0; i < 10; ++i ) + { + eStr + << pnts[ i ][ 0 ] << " " + << pnts[ i ][ 1 ] << " " + << pnts[ i ][ 2 ] << " " << i % 3 << std::endl; + + } // rof + std::ofstream oStrE( oute.getValue( ).c_str( ) ); + oStrE << eStr.str( ); + oStrE.close( ); + + } // fi + } + 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/lib/fpa/Filters/Image/ExtractSkeleton.h b/lib/fpa/Filters/Image/ExtractSkeleton.h index 9e7f948..2312403 100644 --- a/lib/fpa/Filters/Image/ExtractSkeleton.h +++ b/lib/fpa/Filters/Image/ExtractSkeleton.h @@ -97,6 +97,8 @@ namespace fpa itkGetConstMacro( Seed, TIndex ); itkSetMacro( Seed, TIndex ); + itkGetConstMacro( EndPoints, std::vector< TIndex > ); + public: virtual itk::ModifiedTimeType GetMTime( ) const override; @@ -136,6 +138,7 @@ namespace fpa typename TDistanceMap::Pointer m_DistanceMap; bool m_SeedFromMaximumDistance; TIndex m_Seed; + std::vector< TIndex > m_EndPoints; }; } // ecapseman diff --git a/lib/fpa/Filters/Image/ExtractSkeleton.hxx b/lib/fpa/Filters/Image/ExtractSkeleton.hxx index 38b3df4..5904cd4 100644 --- a/lib/fpa/Filters/Image/ExtractSkeleton.hxx +++ b/lib/fpa/Filters/Image/ExtractSkeleton.hxx @@ -172,6 +172,7 @@ GenerateData( ) end_points, this->m_DistanceMap->GetOutput( ), mst, dijkstra->GetSkeletonQueue( ) ); + this->m_EndPoints = end_points; // Compute symbolic branches typedef std::map< TIndex, TIndex, typename TIndex::LexicographicCompare > _TTags; -- 2.45.1