From e9083d9f5f381f258f994fa9bbbe39a897f97c5b Mon Sep 17 00:00:00 2001 From: =?utf8?q?Leonardo=20Fl=C3=B3rez-Valencia?= Date: Tue, 14 Feb 2017 17:38:56 -0500 Subject: [PATCH] ... --- data/test_skeleton_2D.png | Bin 0 -> 4886 bytes examples/CMakeLists.txt | 3 +- examples/Skeleton_00.cxx | 49 ++++ lib/fpa/Base/Algorithm.h | 2 +- lib/fpa/Base/Dijkstra.h | 14 +- lib/fpa/Base/Functors/GaussianModel.h | 2 +- lib/fpa/Base/Functors/Inverse.h | 2 +- lib/fpa/Base/RegionGrow.h | 10 +- lib/fpa/Config.h | 14 -- lib/fpa/Image/Algorithm.h | 14 +- lib/fpa/Image/Dijkstra.h | 2 +- lib/fpa/Image/Functors/SimpleDijkstraCost.h | 2 +- lib/fpa/Image/Functors/SimpleNeighborhood.h | 2 +- lib/fpa/Image/MinimumSpanningTree.h | 8 +- lib/fpa/Image/RegionGrow.h | 2 +- lib/fpa/Image/SkeletonFilter.h | 13 +- lib/fpa/Image/SkeletonFilter.hxx | 243 ++++++++++---------- 17 files changed, 214 insertions(+), 168 deletions(-) create mode 100644 data/test_skeleton_2D.png create mode 100644 examples/Skeleton_00.cxx diff --git a/data/test_skeleton_2D.png b/data/test_skeleton_2D.png new file mode 100644 index 0000000000000000000000000000000000000000..17feca40a907c02414f327462964bdbad74f4d08 GIT binary patch literal 4886 zcmZ`-c{o(<`yU)Slc6(Xh!khWs7bG~RkYyD!dOF#EkcyEkTvSHoW>YqX(*&ckv2=> zl{d;#q_LDnma>#=g*Kur$@i%1_xt|!o$ET+`ON3O@6Y`_=RD7SKQr-mw&rq6S1iTh zaB>zLmIDqa1utG8M&N#_I+tQMd;rHe2#1qZSiB@~7cMDd6p#&mPj=pwt_p+fyZ&`K6@55_E{SQ5@YY56 z%0pZ#W-x^WSL^NA8+s)$1Mm(68AQxL-)bT5I28?(CH_QUQhyAK5yG`a7%&J-W7^=qK@K1OYRqTmcTtV64isQ^X`6 z9Xt(bNkLqV=_4}Aza^&cE)cT$pA|6b?i7aV^?#2&fca8|LDW9@dlr4XP>j9s4C3{@ zbjQ$^-^dI@UCoQ!dJIMUMoJi3VNK-q4gAvMQo5QTLEpfz1D6tEwa9hoVb(?TtqSA1 zEHT7#G-439Uk90`=QIe9+*V=8;G70oY@4hD{0KR#K@eLf>##OpPCcC#7}79$4E>#p zxsN{~3OsfwFpp=)kv;iJftWIf$n&sTV8mwsP0r>IAIGwqV~qJ2T{39^WKJeg$>NOy zJ1)G^6C{--QeWr@bXe7N2C=svWMi*@2k*5MP>^&Fwl5ueeo35F=rbA?cDyJzGJ1zl;#0C~GTy@QgE&N&@6% zmYlH=Tkh86A8;w|I}-5WJZYZw7S#hrEg-@D(=`(B?(4!8T1a1FK*xFcq}}E`-$E@C zE9jkkzm_CEr`paG4PZ@yISqd1(9u0c=!eo0T~?7|65Lhj1&EIiQTbV?;XV+VOw(v& zLC+jDP~YmM#5#b!3~7$;*Qq1Z<&fr>$50(D8oml5rbJB&Kim#N&x1oSKa+F%ZG!aE z1R$o^4v#7m)+t@mm0#MsQB_a_z4Uy72Gy%%|0Bb~CTs|oPOZ7ORKNd;!r+``EM>e4 zKjC)y@wpQtz)YF7r|rW4^pxXu9lHCE6G+FvmnF)IMHR{Ln7!&VORZD#fsgtKnU_-3 z<2k~^mgIUi>TyZ{U}NDq6@R{Mr>ufU->oG9F#~CoJp~t+lIzV-kMCKEK+Lx^3P7Y1 zD14~ARtSS!(Ii9>9kr-|@=X^yv=!R?Bm0w{t`QYZd2+ssoq*R}QM}wM32b5Er)50a z6~R#IvB=LWAAIcqH0Wais7g!u2gf*X=n5T$iH>(_}H5Y_R% zLwT>(Ebw{eS9#!aVFi}^GL!Z^lm`<=4vVDK#SoWLzFO~fEi`D)j|b<2GSKd^mte4y1Gp~e^9kXD{l(5FIM)r+?q!NL#`Xp$tn+r!SCI1icJ@G zN45wa%=>m-3^b90lF`I*M7V7zn^V@k1PM%wqu{4t&W#B>+7)%6VV4@u`{il4yixy& za8=ao@wD?e)NU}?XP}|TT}Sg2Q%klnXmM&lgNqs^&;sd5fVWF$dW)9jt?aHa7nTlX=8AsACMP}9;Jleg)EoTw zdU9o<2B&|bS#Pj9Z&mkJb7AvP<^fT)m73A4nXqjr)0cCnOgnOL{WA9VP0uzQJ>Paz zQG*|upE0qc9DP+%7G{l#pzBnA**Y>};GPNQOr0p<0G=R-ZEY(5SN<}@8kbcdkKj;4 zT00NEfm^gm@*qGI6|k~#=sJSx0KX{+pmJdk5gR~ocu67y(#d_E%oc9X*QS9i&ip@# zFb-kyK?=xJWR@3Y0yDABZLTrxRP73})K}f)XI4Mu1xUeLxci{S>Iqbi@oCE`sNd1oA z$pw)YMl2{6NvjM21wSbK_3;TK0Ow)%eG+RsIb@qp6)=UIRo}S}o|)MD5!p zND>t8dOtx)A5S%P`RwfIwiiUTR##sW`7A%VE**Ef)OY^l)V}Ge<;aUqdL`>aTkTzM zvF!cQ!*;#vFVaE+x|NGm+Y{n)X?#^q`|k=2CNrshllX}}eG6M?sfYan#;t$)8-_+N2ArWy=$7H-GyEIMZl zNPy@v1uuUtfb_H`odI;*gkk3+3{NQ_MoprWLkfvo^QrL=e6 zP^G{NDA(`&dyhGE-{n3R^7ezeim`R!y)*A-52>T=7zAYONMDrzYu|~5!69!+#5$;? zvlv^CXCWA92zD>;fVk>~`}iV?2>z1ZPpn%Hpy$Rl@~TN&m5`rh>Jax7ESgr4Poi*& zHw}(VM0|al3je(8?HLW%F9%+anYTr5%nu6#ORA-eWUlMXwpr@QnGrBfJ-f(_!7u zvKIqTyO)}A@X(S5=*Nav*@H&avlafeB{|N;YWgfb#EQ<1kbG}~Sbi$i%MB?maC&l~ zD)r}P8ee>lLf7V$Mhb_%uHkDp1+s+ri%Ht~hqGy%uKJ=mNfsmr(S}vL?%=7+OLMoK zhINCq2W~7xC>!Ei4>!47d%6oX0EOLoD)~F-1vZy}29wa_ljtM}U5SH$?SvG|LTUUN z8I=4C5m?YAX#-J?Vv zrJ!Y$F0u?tEb31)2r~CK3Ef=Fu3JS5r|SWyf%vd_558*!FjMvEfgF(vGz6_IO5Guc z9@(p^#MEI3r$uK7%WPBf-m~)o=kQ8TKA`K(;KGVhw1!c|kkMEH%9@IlC89Mp*rvlC zQErPgEn_?3_U_Y6`&e*xjVK3PElVl{ys-VuheiTSS88o>BwiugS_o4wmmMSyW;<&| zkmTR|Gx~6%Xsxz$*rfO6Blq_aSixNRTFo|2gqZBA zd(Q(#PIUdukRu;jPMK1mb;|D`;vFUnoR@{qjO_NlKt+Tx68&XwK252Oo~X?D{uj1M z%Lj@NU~N6~pFH=ngVIc})!3BYS})vhi+5B!xn$AOoWluc`pK5i9}C4OD_q|s^n?IY zDwhbRm^@*|wQ`ozY_Jt#Ppmxg9Ne>(@Hy zvFd-e@(7Hl@#D^)wE}WT!ZtZs zsY6djH%j&Hq=65{1!~TtN){K=Tg!IouSLqNJ6t~1SB&_~-!h7sGEoe0EjljA5#|U$ zB3fE~T$i>mw$-xZX44P1OL(HXs-;k3O;6?*&t3iHXfr50E129&15y3z*NCk`?v(bt z-~k^}p^L7{jh&-LsyU*Ghc1yVRH)L@hqiLa7RuB(bk`p5#tj*=C|`YT>Nw4?QE~sq zFDh(@SiV)n2sx9s6RP8IQjYKn{6d`K8y}Xn)aDrmU9WE!>^3s1m$t#IKV6CW9mkmR znzvLn5-3tT2Jg=uADef3amy~@hW=QFw^el$v`%T(Wg`nx zRf^muBX$r@ME(6dxxF2ul^t9m)jj!r>ko)@@7ss|wK-&+(y$8Y8c^)%d5se1)m?0v zHvdja|1_L(OU?J4V=30GCtqmtS6Ow=f?l2X{nv%8^k>jUQ#9SheSZzofw5otpC5NKY=_HcmEqEc7tjYdZoOc1mZnGr0(43 zDTTTZjSoICsYs3R0$Olj-x9^SH3-&dJaQ=c_u%|K#|h&vS@n13cx=>Kv&;2g7uk{v zk=kve?{pghgH7JOdQMO8N=fV_m;0hs)8i49>D?h2+t;nnx&ol*ZNnxXnvd!B@0N^} zcAlGesxx~7IM)+^nQY2o&44yV;Pz*D!HHtrx20#w+U=VACw~05N+%qWW%RmE7xB14 zf75=7anQV*Q=P;hqJcmB{>`2B-Jm3}Q5 zq*MDjoPB7YKKp~`nHqn`@9m$YRW{)Ra|++t6?-8s`W+u8NZ%x2eH4dvTS>>o7K+^s zdhaeCqvW0kqS0x;3Fomz;@Aivp}BY%kmSNTcL52@ISyp{s(<<8L<`6ig->$Z)O7KZ z9Oy`buwP*{b}&TUhz@_(&_HYgPHLbeDMFeaEPF`DjZ`-yk48t!69-*NT$*043|6n56o%+%ciJeo^Czv6%Wcl3XJ*Mv!MZ;!bXm(jgQS+H$c JcNiY%e*o#7=kovn literal 0 HcmV?d00001 diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 4a8ae35..acb1bc5 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -5,6 +5,7 @@ IF(BUILD_Examples) _examples RegionGrow_00 MoriRegionGrow_00 + Skeleton_00 ) INCLUDE_DIRECTORIES( ${PROJECT_SOURCE_DIR}/lib @@ -12,7 +13,7 @@ IF(BUILD_Examples) ) FOREACH(_example ${_examples}) ADD_EXECUTABLE(fpa_example_${_example} ${_example}.cxx) - TARGET_LINK_LIBRARIES(fpa_example_${_example} ${ITK_LIBRARIES}) + TARGET_LINK_LIBRARIES(fpa_example_${_example} ${ITK_LIBRARIES} ${VTK_LIBRARIES}) ENDFOREACH(_example) ENDIF(BUILD_Examples) diff --git a/examples/Skeleton_00.cxx b/examples/Skeleton_00.cxx new file mode 100644 index 0000000..5ddb5e4 --- /dev/null +++ b/examples/Skeleton_00.cxx @@ -0,0 +1,49 @@ +#include +#include +#include +#include +#include + +// ------------------------------------------------------------------------- +typedef itk::Image< double, 2 > TImage; +typedef itk::ImageFileReader< TImage > TReader; +typedef itk::SignedMaurerDistanceMapImageFilter< TImage, TImage > TDMap; +typedef itk::MinimumMaximumImageCalculator< TImage > TMinMax; +typedef fpa::Image::SkeletonFilter< TImage > TFilter; + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + if( argc < 2 ) + { + std::cerr + << "Usage: " << argv[ 0 ] + << " input_filename" << std::endl; + return( 1 ); + + } // fi + std::string in_fname = argv[ 1 ]; + + TReader::Pointer reader = TReader::New( ); + reader->SetFileName( in_fname ); + + TDMap::Pointer dmap = TDMap::New( ); + dmap->SetInput( reader->GetOutput( ) ); + dmap->InsideIsPositiveOn( ); + dmap->SquaredDistanceOff( ); + dmap->UseImageSpacingOn( ); + dmap->Update( ); + + TMinMax::Pointer minmax = TMinMax::New( ); + minmax->SetImage( dmap->GetOutput( ) ); + minmax->Compute( ); + + TFilter::Pointer filter = TFilter::New( ); + filter->SetInput( dmap->GetOutput( ) ); + filter->AddSeed( minmax->GetIndexOfMaximum( ), 0 ); + filter->Update( ); + + return( 0 ); +} + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/Algorithm.h b/lib/fpa/Base/Algorithm.h index 297727a..da7bee8 100644 --- a/lib/fpa/Base/Algorithm.h +++ b/lib/fpa/Base/Algorithm.h @@ -75,7 +75,7 @@ namespace fpa Algorithm( ); virtual ~Algorithm( ); - virtual void GenerateData( ) fpa_OVERRIDE; + virtual void GenerateData( ) override; // Particular methods virtual bool _ContinueGenerateData( ); diff --git a/lib/fpa/Base/Dijkstra.h b/lib/fpa/Base/Dijkstra.h index 5633f22..aa62f90 100644 --- a/lib/fpa/Base/Dijkstra.h +++ b/lib/fpa/Base/Dijkstra.h @@ -56,16 +56,16 @@ namespace fpa Dijkstra( ); virtual ~Dijkstra( ); - virtual void _AfterGenerateData( ) fpa_OVERRIDE; + virtual void _AfterGenerateData( ) override; - virtual void _UpdateResult( const _TQueueNode& n ) fpa_OVERRIDE; + virtual void _UpdateResult( const _TQueueNode& n ) override; virtual bool _UpdateValue( _TQueueNode& v, const _TQueueNode& p - ) fpa_OVERRIDE; - virtual unsigned long _QueueSize( ) const fpa_OVERRIDE; - virtual void _QueueClear( ) fpa_OVERRIDE; - virtual void _QueuePush( const _TQueueNode& node ) fpa_OVERRIDE; - virtual _TQueueNode _QueuePop( ) fpa_OVERRIDE; + ) override; + virtual unsigned long _QueueSize( ) const override; + virtual void _QueueClear( ) override; + virtual void _QueuePush( const _TQueueNode& node ) override; + virtual _TQueueNode _QueuePop( ) override; private: // Purposely not defined diff --git a/lib/fpa/Base/Functors/GaussianModel.h b/lib/fpa/Base/Functors/GaussianModel.h index a74e517..8d32f51 100644 --- a/lib/fpa/Base/Functors/GaussianModel.h +++ b/lib/fpa/Base/Functors/GaussianModel.h @@ -45,7 +45,7 @@ namespace fpa */ public: - virtual TOutput Evaluate( const TInput& x ) const fpa_OVERRIDE; + virtual TOutput Evaluate( const TInput& x ) const override; protected: GaussianModel( ); diff --git a/lib/fpa/Base/Functors/Inverse.h b/lib/fpa/Base/Functors/Inverse.h index dc8c7a2..af68022 100644 --- a/lib/fpa/Base/Functors/Inverse.h +++ b/lib/fpa/Base/Functors/Inverse.h @@ -33,7 +33,7 @@ namespace fpa itkSetMacro( NegativeValue, _TOutput ); public: - virtual TOutput Evaluate( const TInput& x ) const fpa_OVERRIDE; + virtual TOutput Evaluate( const TInput& x ) const override; protected: Inverse( ); diff --git a/lib/fpa/Base/RegionGrow.h b/lib/fpa/Base/RegionGrow.h index 2ef0703..a8fcc93 100644 --- a/lib/fpa/Base/RegionGrow.h +++ b/lib/fpa/Base/RegionGrow.h @@ -47,11 +47,11 @@ namespace fpa virtual bool _UpdateValue( _TQueueNode& v, const _TQueueNode& p - ) fpa_OVERRIDE; - virtual unsigned long _QueueSize( ) const fpa_OVERRIDE; - virtual void _QueueClear( ) fpa_OVERRIDE; - virtual void _QueuePush( const _TQueueNode& node ) fpa_OVERRIDE; - virtual _TQueueNode _QueuePop( ) fpa_OVERRIDE; + ) override; + virtual unsigned long _QueueSize( ) const override; + virtual void _QueueClear( ) override; + virtual void _QueuePush( const _TQueueNode& node ) override; + virtual _TQueueNode _QueuePop( ) override; private: // Purposely not defined diff --git a/lib/fpa/Config.h b/lib/fpa/Config.h index c59b1be..9728f0b 100644 --- a/lib/fpa/Config.h +++ b/lib/fpa/Config.h @@ -7,20 +7,6 @@ * ========================================================================= */ -#if __cplusplus >= 201103L -# define fpa_OVERRIDE override -# define fpa_DELETE_FUNCTION =delete -# define fpa_NULLPTR nullptr -# define fpa_NOEXCEPT noexcept -# define fpa_HAS_CXX11_STATIC_ASSERT -# define fpa_HAS_CXX11_RVREF -#else -# define fpa_OVERRIDE -# define fpa_DELETE_FUNCTION -# define fpa_NULLPTR NULL -# define fpa_NOEXCEPT throw() -#endif - #endif // __fpa__Config__h__ // eof - $RCSfile$ diff --git a/lib/fpa/Image/Algorithm.h b/lib/fpa/Image/Algorithm.h index 70a87a6..d613ce3 100644 --- a/lib/fpa/Image/Algorithm.h +++ b/lib/fpa/Image/Algorithm.h @@ -40,13 +40,13 @@ namespace fpa Algorithm( ); virtual ~Algorithm( ); - virtual void _BeforeGenerateData( ) fpa_OVERRIDE; - virtual void _InitMarks( ) fpa_OVERRIDE; - virtual void _InitResults( const TOutput& init_value ) fpa_OVERRIDE; - virtual bool _IsMarked( const TVertex& v ) const fpa_OVERRIDE; - virtual void _Mark( const _TQueueNode& n ) fpa_OVERRIDE; - virtual TFrontId _GetMark( const TVertex& v ) const fpa_OVERRIDE; - virtual void _UpdateResult( const _TQueueNode& n ) fpa_OVERRIDE; + virtual void _BeforeGenerateData( ) override; + virtual void _InitMarks( ) override; + virtual void _InitResults( const TOutput& init_value ) override; + virtual bool _IsMarked( const TVertex& v ) const override; + virtual void _Mark( const _TQueueNode& n ) override; + virtual TFrontId _GetMark( const TVertex& v ) const override; + virtual void _UpdateResult( const _TQueueNode& n ) override; private: // Purposely not defined diff --git a/lib/fpa/Image/Dijkstra.h b/lib/fpa/Image/Dijkstra.h index 253f973..96c642b 100644 --- a/lib/fpa/Image/Dijkstra.h +++ b/lib/fpa/Image/Dijkstra.h @@ -37,7 +37,7 @@ namespace fpa Dijkstra( ); virtual ~Dijkstra( ); - virtual void _BeforeGenerateData( ) fpa_OVERRIDE; + virtual void _BeforeGenerateData( ) override; private: // Purposely not defined diff --git a/lib/fpa/Image/Functors/SimpleDijkstraCost.h b/lib/fpa/Image/Functors/SimpleDijkstraCost.h index d9a0405..bad9500 100644 --- a/lib/fpa/Image/Functors/SimpleDijkstraCost.h +++ b/lib/fpa/Image/Functors/SimpleDijkstraCost.h @@ -38,7 +38,7 @@ namespace fpa public: virtual TOutput Evaluate( const TIndex& a, const TIndex& b - ) const fpa_OVERRIDE; + ) const override; protected: SimpleDijkstraCost( ); diff --git a/lib/fpa/Image/Functors/SimpleNeighborhood.h b/lib/fpa/Image/Functors/SimpleNeighborhood.h index 1fd8547..9f777b0 100644 --- a/lib/fpa/Image/Functors/SimpleNeighborhood.h +++ b/lib/fpa/Image/Functors/SimpleNeighborhood.h @@ -36,7 +36,7 @@ namespace fpa itkSetMacro( Order, unsigned int ); public: - virtual TOutput Evaluate( const TIndex& center ) const fpa_OVERRIDE; + virtual TOutput Evaluate( const TIndex& center ) const override; protected: SimpleNeighborhood( ); diff --git a/lib/fpa/Image/MinimumSpanningTree.h b/lib/fpa/Image/MinimumSpanningTree.h index 4a4e0e6..7d3941d 100644 --- a/lib/fpa/Image/MinimumSpanningTree.h +++ b/lib/fpa/Image/MinimumSpanningTree.h @@ -32,15 +32,15 @@ namespace fpa ); public: - virtual TVertex GetParent( const TVertex& v ) const fpa_OVERRIDE; - virtual void SetParent( const TVertex& v, const TVertex& p ) fpa_OVERRIDE; + virtual TVertex GetParent( const TVertex& v ) const override; + virtual void SetParent( const TVertex& v, const TVertex& p ) override; virtual void GetPath( typename TPath::Pointer& path, const TVertex& a - ) const fpa_OVERRIDE; + ) const override; virtual void GetPath( typename TPath::Pointer& path, const TVertex& a, const TVertex& b - ) const fpa_OVERRIDE; + ) const override; protected: MinimumSpanningTree( ); diff --git a/lib/fpa/Image/RegionGrow.h b/lib/fpa/Image/RegionGrow.h index ad41e80..dbc57f0 100644 --- a/lib/fpa/Image/RegionGrow.h +++ b/lib/fpa/Image/RegionGrow.h @@ -34,7 +34,7 @@ namespace fpa RegionGrow( ); virtual ~RegionGrow( ); - virtual void _BeforeGenerateData( ) fpa_OVERRIDE; + virtual void _BeforeGenerateData( ) override; private: // Purposely not defined diff --git a/lib/fpa/Image/SkeletonFilter.h b/lib/fpa/Image/SkeletonFilter.h index 5a25634..0fbb3d5 100644 --- a/lib/fpa/Image/SkeletonFilter.h +++ b/lib/fpa/Image/SkeletonFilter.h @@ -1,6 +1,7 @@ #ifndef __fpa__Image__SkeletonFilter__h__ #define __fpa__Image__SkeletonFilter__h__ +#include #include #include #include @@ -35,6 +36,7 @@ namespace fpa protected: typedef typename Superclass::_TQueueNode _TQueueNode; typedef std::multimap< double, TIndex, std::greater< double > > _TSkeletonQueue; + typedef std::map< TIndex, TIndex, typename TIndex::LexicographicCompare > _TAdjacencies; public: itkNewMacro( Self ); @@ -48,12 +50,13 @@ namespace fpa SkeletonFilter( ); virtual ~SkeletonFilter( ); - virtual void _BeforeGenerateData( ) fpa_OVERRIDE; - virtual void _UpdateResult( const _TQueueNode& n ) fpa_OVERRIDE; - virtual void _AfterGenerateData( ) fpa_OVERRIDE; + virtual void _BeforeGenerateData( ) override; + virtual void _UpdateResult( const _TQueueNode& n ) override; + virtual void _AfterGenerateData( ) override; - void _EndPoints( std::vector< TIndex >& end_points ); - void _Skeleton( const std::vector< TIndex >& end_points ); + void _EndPoints( std::vector< TIndex >& end_points, _TAdjacencies& A ); + void _Skeleton( const std::vector< TIndex >& end_points, _TAdjacencies& A ); + void _MarkSphere( const TIndex& idx ); private: // Purposely not defined diff --git a/lib/fpa/Image/SkeletonFilter.hxx b/lib/fpa/Image/SkeletonFilter.hxx index 069eea0..d8dc66a 100644 --- a/lib/fpa/Image/SkeletonFilter.hxx +++ b/lib/fpa/Image/SkeletonFilter.hxx @@ -93,36 +93,32 @@ _AfterGenerateData( ) { this->Superclass::_AfterGenerateData( ); + _TAdjacencies A; std::vector< TIndex > end_points; - this->_EndPoints( end_points ); - this->_Skeleton( end_points ); + this->_EndPoints( end_points, A ); + this->_Skeleton( end_points, A ); } // ------------------------------------------------------------------------- template< class _TImage > void fpa::Image::SkeletonFilter< _TImage >:: -_EndPoints( std::vector< TIndex >& end_points ) +_EndPoints( std::vector< TIndex >& end_points, _TAdjacencies& A ) { - typedef itk::ImageRegionIteratorWithIndex< TMarks > _TMarksIt; - - static const double _eps = std::sqrt( double( Self::Dimension + 1 ) ); - auto input = this->GetInput( ); - auto mst = this->GetMinimumSpanningTree( ); auto marks = this->GetMarks( ); + auto mst = this->GetMinimumSpanningTree( ); + auto spac = marks->GetSpacing( ); // Some values - marks->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) ); - marks->SetRequestedRegion( input->GetRequestedRegion( ) ); - marks->SetBufferedRegion( input->GetBufferedRegion( ) ); - marks->SetSpacing( input->GetSpacing( ) ); - marks->SetOrigin( input->GetOrigin( ) ); - marks->SetDirection( input->GetDirection( ) ); + marks->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) ); + marks->SetRequestedRegion( mst->GetRequestedRegion( ) ); + marks->SetBufferedRegion( mst->GetBufferedRegion( ) ); + marks->SetSpacing( mst->GetSpacing( ) ); + marks->SetOrigin( mst->GetOrigin( ) ); + marks->SetDirection( mst->GetDirection( ) ); marks->Allocate( ); marks->FillBuffer( 0 ); - auto spac = marks->GetSpacing( ); // BFS from maximum queue - auto region = input->GetRequestedRegion( ); while( this->m_SkeletonQueue.size( ) > 0 ) { // Get node @@ -136,66 +132,20 @@ _EndPoints( std::vector< TIndex >& end_points ) continue; marks->SetPixel( n.second, 1 ); end_points.push_back( n.second ); + std::cout << this->m_SkeletonQueue.size( ) << std::endl; // Mark path TIndex it = n.second; TIndex p = mst->GetParent( it ); while( it != p ) { - // TODO: tags->SetPixel( it, tags->GetPixel( it ) + 1 ); + this->_MarkSphere( it ); it = p; p = mst->GetParent( it ); } // elihw - // TODO: tags->SetPixel( it, tags->GetPixel( it ) + 1 ); - - /* TODO: use mst directly rather than computing paths - typename TMST::TPath::Pointer path; - mst->GetPath( path, n.second ); - for( unsigned long i = 0; i < path->GetSize( ); ++i ) - { - TIndex idx = path->GetVertex( i ); - typename _TImage::PointType cnt; - input->TransformIndexToPhysicalPoint( idx, cnt ); - double r = double( input->GetPixel( idx ) ) * _eps; - - TIndex i0, i1; - for( unsigned int d = 0; d < Self::Dimension; ++d ) - { - long off = long( std::ceil( r / double( spac[ d ] ) ) ); - if( off < 3 ) - off = 3; - i0[ d ] = idx[ d ] - off; - i1[ d ] = idx[ d ] + off; - - if( i0[ d ] < region.GetIndex( )[ d ] ) - i0[ d ] = region.GetIndex( )[ d ]; - - if( i1[ d ] >= region.GetIndex( )[ d ] + region.GetSize( )[ d ] ) - i1[ d ] = region.GetIndex( )[ d ] + region.GetSize( )[ d ] - 1; - - } // rof - - typename _TImage::SizeType size; - for( unsigned int d = 0; d < Self::Dimension; ++d ) - size[ d ] = i1[ d ] - i0[ d ] + 1; - - typename _TImage::RegionType neighRegion; - neighRegion.SetIndex( i0 ); - neighRegion.SetSize( size ); - - _TMarksIt mIt( marks, neighRegion ); - for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt ) - { - TIndex w = mIt.GetIndex( ); - typename _TImage::PointType p; - marks->TransformIndexToPhysicalPoint( w, p ); - mIt.Set( 1 ); - - } // rof - - } // rof - */ + this->_MarkSphere( it ); + A[ n.second ] = it; } // elihw } @@ -203,67 +153,124 @@ _EndPoints( std::vector< TIndex >& end_points ) // ------------------------------------------------------------------------- template< class _TImage > void fpa::Image::SkeletonFilter< _TImage >:: -_Skeleton( const std::vector< TIndex >& end_points ) +_Skeleton( const std::vector< TIndex >& end_points, _TAdjacencies& A ) { - typedef itk::Image< unsigned long, Self::Dimension > _TTagsImage; - typedef typename TMST::TPath _TPath; + std::cout << end_points.size( ) << " " << A.size( ) << std::endl; + + /* TODO + typedef itk::Image< unsigned long, Self::Dimension > _TTagsImage; + typedef typename TMST::TPath _TPath; + + auto mst = this->GetMinimumSpanningTree( ); + auto skeleton = this->GetSkeleton( ); + + // Tag branches + typename _TTagsImage::Pointer tags = _TTagsImage::New( ); + tags->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) ); + tags->SetRequestedRegion( mst->GetRequestedRegion( ) ); + tags->SetBufferedRegion( mst->GetBufferedRegion( ) ); + tags->Allocate( ); + tags->FillBuffer( 0 ); + for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt ) + { + TIndex it = *eIt; + TIndex p = mst->GetParent( it ); + while( it != p ) + { + tags->SetPixel( it, tags->GetPixel( it ) + 1 ); + it = p; + p = mst->GetParent( it ); + + } // elihw + tags->SetPixel( it, tags->GetPixel( it ) + 1 ); + + } // rof + + // Build paths (branches) + for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt ) + { + TIndex it = *eIt; + TIndex p = mst->GetParent( it ); + TIndex sIdx = *eIt; + typename _TPath::Pointer path = _TPath::New( ); + path->SetReferenceImage( mst ); + while( it != p ) + { + if( tags->GetPixel( sIdx ) != tags->GetPixel( it ) ) + { + // Ok, a new branch has been added + path->AddVertex( it ); + skeleton->AddBranch( path ); + + // Update a new starting index + path = _TPath::New( ); + path->SetReferenceImage( mst ); + sIdx = it; + } + else + path->AddVertex( it ); + it = p; + p = mst->GetParent( it ); + + } // elihw + + // Finally, add last branch + path->AddVertex( it ); + skeleton->AddBranch( path ); + + } // rof + */ +} - auto mst = this->GetMinimumSpanningTree( ); - auto skeleton = this->GetSkeleton( ); - - // Tag branches - typename _TTagsImage::Pointer tags = _TTagsImage::New( ); - tags->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) ); - tags->SetRequestedRegion( mst->GetRequestedRegion( ) ); - tags->SetBufferedRegion( mst->GetBufferedRegion( ) ); - tags->Allocate( ); - tags->FillBuffer( 0 ); - for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt ) +// ------------------------------------------------------------------------- +template< class _TImage > +void fpa::Image::SkeletonFilter< _TImage >:: +_MarkSphere( const TIndex& idx ) +{ + typedef itk::ImageRegionIteratorWithIndex< TMarks > _TMarksIt; + + static const double _eps = std::sqrt( double( Self::Dimension + 1 ) ); + auto input = this->GetInput( ); + auto marks = this->GetMarks( ); + auto spac = input->GetSpacing( ); + auto region = input->GetRequestedRegion( ); + + typename _TImage::PointType cnt; + input->TransformIndexToPhysicalPoint( idx, cnt ); + double r = double( input->GetPixel( idx ) ) * _eps; + + TIndex i0, i1; + for( unsigned int d = 0; d < Self::Dimension; ++d ) { - TIndex it = *eIt; - TIndex p = mst->GetParent( it ); - while( it != p ) - { - tags->SetPixel( it, tags->GetPixel( it ) + 1 ); - it = p; - p = mst->GetParent( it ); + long off = long( std::ceil( r / double( spac[ d ] ) ) ); + if( off < 3 ) + off = 3; + i0[ d ] = idx[ d ] - off; + i1[ d ] = idx[ d ] + off; - } // elihw - tags->SetPixel( it, tags->GetPixel( it ) + 1 ); + if( i0[ d ] < region.GetIndex( )[ d ] ) + i0[ d ] = region.GetIndex( )[ d ]; + + if( i1[ d ] >= region.GetIndex( )[ d ] + region.GetSize( )[ d ] ) + i1[ d ] = region.GetIndex( )[ d ] + region.GetSize( )[ d ] - 1; } // rof - // Build paths (branches) - for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt ) - { - TIndex it = *eIt; - TIndex p = mst->GetParent( it ); - TIndex sIdx = *eIt; - typename _TPath::Pointer path = _TPath::New( ); - path->SetReferenceImage( mst ); - while( it != p ) - { - if( tags->GetPixel( sIdx ) != tags->GetPixel( it ) ) - { - // Ok, a new branch has been added - path->AddVertex( it ); - skeleton->AddBranch( path ); - - // Update a new starting index - path = _TPath::New( ); - path->SetReferenceImage( mst ); - sIdx = it; - } - else - path->AddVertex( it ); - it = p; - p = mst->GetParent( it ); + typename _TImage::SizeType size; + for( unsigned int d = 0; d < Self::Dimension; ++d ) + size[ d ] = i1[ d ] - i0[ d ] + 1; - } // elihw + typename _TImage::RegionType neighRegion; + neighRegion.SetIndex( i0 ); + neighRegion.SetSize( size ); - // Finally, add last branch - path->AddVertex( it ); - skeleton->AddBranch( path ); + _TMarksIt mIt( marks, neighRegion ); + for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt ) + { + TIndex w = mIt.GetIndex( ); + typename _TImage::PointType p; + marks->TransformIndexToPhysicalPoint( w, p ); + mIt.Set( 1 ); } // rof } -- 2.45.0