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{r5a`yb;`&t&k-7$U`)F>2yfW2tDtGYey@p~V&2MIoFBjZ?sWbuxJ4e zhf}iUah!2Dc`)ZN4u`}2(sC=o;c$c?o?9>ur>HvTa=5eSHE=lG{1975?)-OnWp!Pk zUiG93hr?-xa9u;#fd>xw9}K~<5BvLs`1^sdy&-#mjjbc!B}RE64o8;Qa?mYDJ3qaC znXGH0Z0*C%c{7`4zHr|DO|_n$&MO}DuX8b(_t7q8{fJ4fyVuIOG5&{nTMmn7d9Nnk z&eK;Uw=m!~ReIJpgt3&jeu&1mlN@!mxg-HasO=>YD0*af2?Irr4xj=QIXQxRps0N# zxC4qj`JfOK?QjBlAf?tBTmdQEEg&1D*lYzCKuY5_kO5MfUBEMrwS|9-X|IOn6@^sd zPWrgIlo_AWvE6%ivGuO-m=Clx#nVM70@dT@CsDM0Se1VNKl0x=HChAp&& ze6?smd9eEgqA!mKbSL)D)A%hhaeJ1SE&8k~(4DwFPRth7oBW|>|;dO{mjcmKwJA8S(^xGoy~LH3KIb>;x|$^5zrPnkc8d6zl;Uc&L%|M-8=-?=cmWaJ@9?7c~qeJDGQ;v-T&ETApok%Zp%v#i+c zzbVmHa3vq(a) zBaeYk%7M8h@idBzF6MI*7&QoXpOGWSXx9;2g)vy9ifou0BBT*z8goh%nIdJM-G_N4 z0vbhzh?yMAIk`)ss~;4&A+g{}Ryyq|Ear2x7&VCJC73%FO)I~wK#}QSe9nod*M1RW zMmGy+;ub`BeFIk`XN)LY7Q0(yO%hPo#%lZj$fi;*P>D{#LMEp{vj$PyEfU^X!w$(` z7e_WDFyE~oeKy43FJy8G2rN?CE>Y!jQ?}_Xvo{Y9Yj`)KeJi6_n3V{Dh@@xM*3(^|v6qLCrxdGk-&%Ivo4iFUacb?*=QaE3g? zd5alCL$-+c?ulx-cXy3Ixjx#R7}Rk_Icckna8H3gnG^g@xkn#Q$f@F6VL_a6IHw`N z205_J9Q#mGY{V&4O9Gt*zL4xtKTVW%0(8UZSekAl2YKwGi|wgeK&pcnDrj?b53`QK zP(qtyYLGg5G`Iw#RwO;DD4dUAXW)Jy%H*ARo1pM00g5T&6EF>;f%Dis?|jO zZ>I9~2G45+Zt-WUwzoRA2&qZ5$e#p^ncMEz=V>w{J<#eDs%q=hZfTJNk!XLj`XVgp z(#olcee1dur%qfu1&nhW{AbiH>C%oVIOW2X?JF_Az|V<0a23wqk5#~}C3h&{42xIj zD~`C}zRT<2j<3Hs*&|BJ9`yEmti;ADy`B5Mcg3kSo!uJQBEC1Xr}NxlOC=;3OB_KZ z>-)2LrCsyU!)b9;!Z_ev9p%$6>Oc)!w1s}pPJpsTlZTQeQB#M~&fqY9UxNYw4nxfZR%lL zv?Bp*P{{O?s@z=MRc<3G>CeoS{)kOZdaTQPGn#1J_wSYDiUM6;&uFu8U)9YeUF&Ql z&Hb6Xq*eCX=2O;^w*JgLyj!J)k$o#vxZhVlUU}e5+W|FQQDk1m=%zC4MR94E113eT z&_tyM6x7r!6V4ewTFirdV7LZRqltdysh|$HtbAn@hndmah2Sb~E`}E(kTfc2aby1# z6k|faDFmVkfVY$jp*RAb#6p<4Pm;Nk4S9xin8Ta-ClTOKjtHj0Yz-DtEmk0<(phsl z905}RykQiuvP^#i4$Ow}+YiFO0p=uxP{7hAkg#YfB*%t*AJz2#=>((EO-fXeJ7o zK>hQ0Z{e)R_vt)6`16n6!YF5BYdW0ub?3fB2;Ghbh&Bg4qrnUo^6~BM0DTssl}@3~ zE%^d)z0#-MSul4rKrk1^i;=##BWofkRFP8N!xLAw3BaCkJ(!=nzt#m|dO1>h6{poS2yvJOO`Q?NrM zrAon<^d3^(3J5zrqI zbzK$*ktBd4Ey4C3lAI+&UWG4N(O}Z%&VJ$3ed)mHv*NiRjQ>n84)o7&Kz^)zk=g3X)b3M>aqf|#eAkC(D^q{2p^IdvsSHD2Nu;Fz>oSpH z(_xN;a4yL(Z(lZ@*I8dUjprar2z@})_ZET1K0kfKb-*atu=nb0goYWRbzhVFeM#MTD)UG9f)lRN0tv^F22Y zatp8U5kW?NECEoHr#B3#9T|!hW1R6wMG{u+C;)hl2+3_xpl56#-q?DA?Gy`7EtBTJ zOBKllkS~lvB4jYgYO(%W7t%$_Mq5ehh0@)mzHB%B2(t2bB3 zB1iH>v14{7W#E!X%^u8pHQ#ZsFk83noc&be=n=~;b?=4H;L*;X8A_DB3#sF(^ykW( zNCfAxd^eTw$-%9DXKAP;My{vy&8Kngp`#TU-~Td5$CdRK?OINTGXJ^hRk~Y&4YyiY zF-J2w(s)Q-)ci!|r(!r1|e#f%z>4C92GScV=I(@;s7++lsrRR$T_dps_$K-${> zgit(h{h)|KV29ej>!`P0KY8bQw9H{$T&!wtN6XC>%bkrm^*>vMMAoDD5w}lD4U-E_ zwyhdg1=n=Sq7org0CR6vRM%k;w|;rg)q7#JV#<+(^-7BJdmj$1l7GFK4&NIQYr74p z+n!BtE!|?W94&R|aQ{?aKIlJl-8^R8QZ2}%@DQFS$q~aOtfcCY5q)-eon6PZrXQZ? z2_zjYJBi$~SD9aYw)B)?&9LN@cx(+F#!RYSq8n9n#|>x1_e6*$LwZSR_%tn2+XYM9 zd!AyeNmEeRyN*Y()u6>;TekT%uFO!xL^{h;N9bmaYCBhb(d0VEitHl>DVg-mNF9%t zvR{}VaKtsf@qTGb?M<`bEA{Q-t>)JC3XVt}@23YTuj3$FS?{`*ZUR+)Q{Ua`L&Gzk zPp|V6u9^&IxGc;RNwlLpDwZnx)=mw)Acn3#Al28CHeU4}u{AulXoJg-onB4%P1*q4 z_a?Bq`PD#Ax~OaENBtR}?rO=~&4UNtyCxU+^%)H`M#-}OeUs4_ z-}9baC$rdkvT8PpZqj$%!!hJr?|8%y8z0|}whiC^DYYXGZ`$($vlQGY*4Zf{2n z)HR$v5Up)_-79uUjqW;y{zcnr7Vh(zp9{;*47N3nYY#VxK)Cifb_IqvY z`slX2+-6GI9;|X=B%&g{>qy21gB4j9A?(b?u(A6#!$v(@@v#bS(=)Dh)^8xUdLlHL zP2Hy#)TRdA_zY%UsTSR9UyVC>UiSJVivL#OibJ!^U)C9+-WTa2SPrgXB6ZQV{9M<^ z&n3^z_jS`@DMqSJm8 zZo_lL!9hq)Z|*UG7XapMNX~AWhgiMnSy+>3i&$acn4nGDh=AuI`x7KRs!LHE0rMm| z|6NUoCs4#y8^)6(6c{0jdzqf(x=Cd$I$D|3=U(jIoDLnTKLNl+T_`g>pBv1~o$)3K zr~_hfv_oRClZ$K*gS@^x%%r^c(~0w&Qj_Rz-ud zsWbKQ;I2bTaxDkh?$2v-qrtbFatG-o$c51wnhay##HLv^tUI4puv-MouJOJRXz)g} v)~Mq!rlOZ=v=(|WoqlG~|M)ukk9Ce(hP-}d<0@cv&2hHejhtI7Z|wg7>7D2E 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.47.1