From bd0ef351d497debf21da4af475aaec5930b5a77f Mon Sep 17 00:00:00 2001 From: =?utf8?q?Leonardo=20Fl=C3=B3rez-Valencia?= Date: Thu, 23 Mar 2017 16:38:44 -0500 Subject: [PATCH] ... --- data/binary_test_2D_00.png | Bin 0 -> 5192 bytes examples/CMakeLists.txt | 1 + examples/Dijkstra_Maurer.cxx | 120 +++++++++++ libs/fpa/Base/Dijkstra.h | 17 +- libs/fpa/Base/Dijkstra.hxx | 129 +++++++++-- libs/fpa/Base/Functors/InvertValue.h | 80 +++++++ libs/fpa/Base/Functors/VertexParentBase.h | 58 +++++ libs/fpa/Base/MinimumSpanningTree.h | 77 +++++++ libs/fpa/Base/MinimumSpanningTree.hxx | 239 +++++++++++++++++++++ libs/fpa/Image/Dijkstra.h | 91 ++++++++ libs/fpa/Image/Functors/VertexIdentity.h | 67 ++++++ libs/fpa/Image/Functors/VertexParentBase.h | 63 ++++++ libs/fpa/Image/MinimumSpanningTree.h | 71 ++++++ 13 files changed, 993 insertions(+), 20 deletions(-) create mode 100644 data/binary_test_2D_00.png create mode 100644 examples/Dijkstra_Maurer.cxx create mode 100644 libs/fpa/Base/Functors/InvertValue.h create mode 100644 libs/fpa/Base/Functors/VertexParentBase.h create mode 100644 libs/fpa/Base/MinimumSpanningTree.h create mode 100644 libs/fpa/Base/MinimumSpanningTree.hxx create mode 100644 libs/fpa/Image/Dijkstra.h create mode 100644 libs/fpa/Image/Functors/VertexIdentity.h create mode 100644 libs/fpa/Image/Functors/VertexParentBase.h create mode 100644 libs/fpa/Image/MinimumSpanningTree.h diff --git a/data/binary_test_2D_00.png b/data/binary_test_2D_00.png new file mode 100644 index 0000000000000000000000000000000000000000..fbe9d39bd59189ecf7f8ef87e48abde2f3d01953 GIT binary patch literal 5192 zcmaJ_c{G$?_rKE*W zWDUhPDnu%>ge;Y%l-=)4-}5`a|9Hs6i%_YQ*g`DvU@kESG4jyKyy9oe*c&&*6dMxzg_x?|5Us{9> zJ8(DG*RzJK61Hd^RY*Nx-xKxtX5p2@rq}6*O>X7M{WbM>OhJm-h93j4L9N~2cB8#6 z@A25fUyr5*4Mq2w!x0Y>YOJVXy^+{VP3U-flE)6!p#5 zTzM-1Hi$=gf~STO4jqdISDK}=xvt(ZG3Q)Z6kM+;5in+Fqa z9~O!l7rZr?P4F6VWjs%>a3@h0IKo)Vxfdh$KD6ls;_Ie4XdG##&2`M1_?#du&fuj_$gT@;h3GX zMQIvBP%f{{uooe%QvyZ;1VdbjpN#BRdbx&t;oFJB985fIw!VH-BQeaib?(EXLC z!Mn{Zx4O=X$imb=?{pK?5-7>Uf+LZ^Z9$|?6R!GvJNn9{#Rn@`I(>5Kia;v0Wr@(5 zW%V=UpggJ52qtR|UNhu&vqxzCC%X7|NrMCm$E+uAe&NoCX&HkwJv4hMj!j(~^xf$H!^-k%{L8~DeAmOq^Yh&G@s#C;5 zvBjN}HL-co9g}$;#n0?gCCm!VJ8PA>>+5yV1BXOc^WQG_^_3>smuf~YN7u?m+IU1) zIhK@KT;3BE+xyckv04AgaGfIm_hvo$uUxc?8uvSug{^9U%VQUv*6(KNYq$}` zDCMqueej*XUmkS-b<4O7k!9(SW&wA+SKYby9$lOeZy&jeRR| zfybIuJZ55SG)tgBMs7?ZNJvG-cZoN*eG~)J^D~4t6+mhG#*RWE$+L8pc+(tSuzThy4|$dhB^ZjXxBn2 z$VY`72b;yHLhELt$$ntqE8Q3v>|bGG01?C$%#RzvT1${zUS}4fab4UXd^QdleED7X zC`Lk)mQodExk$>SvXnTu|9-NBK&z+BqAKZ=#O|^{wM(#ZFsLT5xa#qDDMJVl~ zII#)n^u4>avByfS@K41(mmdEA-sM6@BGz`KJLrSLhB~iwvFjTE7T9Z z**iX2F0hmK{@tl$``?ag*|Z|OrXsCu&}x~uMNk=iG|5Wj$J9i*XAkO~!A)F3N-FVN zqv<_BQ7G6O)sIsfiVO2bqAVFwcF@3=&J91B`BkLG^!tDlqr))>u4&%; z^$&wnsFc`B27((|6TEe_o*++`83x6W*D*)lsG_zN7ZsW$&r{6YfAwzXY?j)ytb_vC zQ_XE8v`k|C6F%*uG}z&UaDeLzwKW=XF)Q~$yu+_Cr)ghjrHw}-we^c01ge$_8uyN? zN*x~0uX=Icb8SyGb1OQ?HiU;#=C07-7(|}+rRYW#a`&G1pF3g(nBGZWnkAZR-duwA z6IngiRq6a0>nNIQ{g@NhiPM)z;S^7ZPLJHWnU!DVyql2~MoGU2l(Nj?p(o})@BbCBGOSq))u|3T`%r1Vp)8tkXAHJK}*A@T3R1mTx(M>h1)!*6C`TKs-|T=hct#2l@SZN zcdhKvDfybDMW>_?!r3g9(`QHJd=hXJYuYQ>eiGz^aGWr# z{$P12Vgxr{2I<0P;cDB8a}NI6;MYP&w&LPFbY|5fqWZ z05EVc76Ew0<-(g08icDhg28d_zo6Hntj}#c3na_raj>B;@e%?flcq2q4ufDAbWKNa zV4ktsDgjFm?lmENuW|6)=RK$tV9`(6>b>V8>>JS;u35$BGe*)kAcccOM+w zbRk8_ct4iiA(?{!XbcrX(6;eAVia)a%{bn|r!Px6rLY3WlnKZ1u+*6|SUpUK5AtE( zS|VT&f!+?*o^7`g)KEby?AwzOkdOVjgy-hGJQc4Q2F|)BfP6W*B=D{j!bNCVjEO_2 z{k4r8uw!4ke_udz5Pu8;7lC}TqJEgfr)>%x)y&O%Ae;bz#ta=@(qV0>uLmzj3M{ft zy(VD**f1s|qa6=7Y5zO<9xSQz;a9U#J z0pNtj%{W9Oo)<0+e!9r3v>C7C&66Hm|7Q5ooz>vyq2q5vBF4<*_=upKY)Q$=WjqcB zx_g;zjijzRZee`0C|0@H+F^`ME1^4!}t(Z_IB#*V13o{W>AA4L#l~Nb2PKo64Y9 zZjoaUj#o#$g7Dv^;eT8Wxd?!ci166F)3Cvx(k>|4G526(wVL?%A~DExUx7zW{lus( z27s|ZgE$siP6Pud1!0T>2ZmT7msrYXK#m86MSD79F>>bQ_K==#jw$audm- zUX;NhgGfv9vkO8}XrDxkh*h!uDznr(5lMHPTwREzbBU}K<($gUS#ak$Piq+yE;*|> zq6WGU(mdHktvV&$ar>ETjNI&ylS`RX8WLEnGb*shz9ejO+DbrkQ3RM*<2B;YeEMY? z^7FK97~G^d^&vWJCwz$f)Rld}Wn z+v#nHI6rQO{SX)FXy8ViY<(VstO~c93a!o3KNgWkzBWC3p}OcIn2Oy0$}P^6>&--* zo`n7?+2OaO6uK2vVYLWF#UQ?)#w7J`^jwaq8#W&>e&cyhpbhc=h)yR;6>@bMB48=}7wv?u)1v7~3~Mt85o#e>3IMrP<#~PqR-}xf))oc-`^) z#8*O)(}}zpM!Yt%bS6|H<#byoBH4D_=8@N2ql2a6s4A~tymw3+YyUN&82P!je!%oLJh1U-+Q_h2x zC0D;LGnPpY{1sfsWl z3!BvJM~>`oD$YMt2n@)glZ!g}nsVLdUpZ{#@kQ30&Fa@;*~|x7>5%AL6?W6xQB9W3 z_aLk7Bsv$zzOW|jLRic=-?KS$%7^fz-1dh`Jb}_;!x&pXEzLfuYWF~eL5xRU5^QE$ zBiP0FY#*rL7*LQ=s!m=JtEKL%_m4Aba~p^wtHGzbm0sj-+-sI;m+4QSv|PFRc&28J z@IfaoqI|l9b@lM9`Rksm(_ICIbE6JBiGfFbP53m zRGMA(735<;I5rFtC=f}LYla7u5}}e=V$%@6;fwub*RP74FT$p3jFbD(dKinLxWjtG z8PfYTn|_$V`-1nUBN}WPrZ<{jp=q#z3DYoIJQMZ|K4+Lp=)e)dfCvt}AK|ro7zV#j zyMwV11tQ_F`Yq=XdGIeXfKb3#WK&<}uUIf}I2Lh+{&V!dasBNgC;sdAJhee_G=CJ2 z)AJ0-eYcSaa*5?g-qZj0&-VWf)1Puf^ofj&#1`p`Y*#7{?v7li9 literal 0 HcmV?d00001 diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index f0c8826..297746a 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -5,6 +5,7 @@ IF(BUILD_EXAMPLES) RegionGrow_Tautology RegionGrow_BinaryThreshold RegionGrow_Mori + Dijkstra_Maurer #CreateMoriInputImage #BronchiiInitialSegmentationWithMori #BronchiiInitialSegmentationWithBinaryThresholdRegionGrow diff --git a/examples/Dijkstra_Maurer.cxx b/examples/Dijkstra_Maurer.cxx new file mode 100644 index 0000000..8875959 --- /dev/null +++ b/examples/Dijkstra_Maurer.cxx @@ -0,0 +1,120 @@ +#include +#include +#include +#include +#include + +#include +#include +#include + +// ------------------------------------------------------------------------- +static const unsigned int VDim = 2; +typedef short TPixel; +typedef double TScalar; +typedef itk::Image< TPixel, VDim > TImage; +typedef itk::Image< TScalar, VDim > TScalarImage; +typedef itk::ImageFileReader< TImage > TReader; +typedef itk::ImageFileWriter< TScalarImage > TWriter; +typedef fpa::Image::Dijkstra< TScalarImage, TScalarImage > TFilter; +typedef itk::MinimumMaximumImageCalculator< TImage > TMinMax; +typedef itk::SignedMaurerDistanceMapImageFilter< TImage, TScalarImage > TDMap; + +typedef fpa::Image::Functors::VertexIdentity< TScalarImage, TScalar > TVertexFunc; +typedef fpa::Base::Functors::InvertValue< TScalar, TScalar > TValueFunc; + +typedef TFilter::TMST TMST; +typedef itk::ImageFileWriter< TMST > TMSTWriter; + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + // Get arguments + if( argc < 5 + VDim ) + { + std::cerr + << "Usage: " << argv[ 0 ] + << " input_image output_image output_mst stop_at_one_front"; + for( unsigned int i = 0; i < VDim; ++i ) + std::cerr << " s_" << i; + std::cerr << " ..." << std::endl; + return( 1 ); + + } // fi + std::string input_image_filename = argv[ 1 ]; + std::string output_image_filename = argv[ 2 ]; + std::string output_mst_filename = argv[ 3 ]; + bool stop_at_one_front = ( std::atoi( argv[ 4 ] ) == 1 ); + + TReader::Pointer reader = TReader::New( ); + reader->SetFileName( input_image_filename ); + try + { + reader->Update( ); + } + catch( std::exception& err ) + { + std::cerr << "Error caught: " << err.what( ) << std::endl; + return( 1 ); + + } // yrt + + TMinMax::Pointer minmax = TMinMax::New( ); + minmax->SetImage( reader->GetOutput( ) ); + minmax->Compute( ); + + TDMap::Pointer dmap = TDMap::New( ); + dmap->SetInput( reader->GetOutput( ) ); + dmap->SetBackgroundValue( minmax->GetMinimum( ) ); + dmap->InsideIsPositiveOn( ); + dmap->UseImageSpacingOn( ); + dmap->SquaredDistanceOff( ); + + TFilter::Pointer filter = TFilter::New( ); + filter->SetInput( dmap->GetOutput( ) ); + filter->SetStopAtOneFront( stop_at_one_front ); + + for( int i = 5; i < argc; i += VDim ) + { + if( i + VDim <= argc ) + { + TImage::IndexType seed; + for( int j = 0; j < VDim; ++j ) + seed[ j ] = std::atoi( argv[ i + j ] ); + filter->AddSeed( seed ); + + } // fi + + } // rof + + TVertexFunc::Pointer vertex_func = TVertexFunc::New( ); + filter->SetFunctor( vertex_func ); + + TValueFunc::Pointer value_func = TValueFunc::New( ); + value_func->SetAlpha( 1 ); + value_func->SetBeta( 1 ); + filter->SetFunctor( value_func ); + + TWriter::Pointer writer = TWriter::New( ); + writer->SetInput( filter->GetOutput( ) ); + writer->SetFileName( output_image_filename ); + + TMSTWriter::Pointer mst_writer = TMSTWriter::New( ); + mst_writer->SetInput( filter->GetMinimumSpanningTree( ) ); + mst_writer->SetFileName( output_mst_filename ); + + try + { + writer->Update( ); + mst_writer->Update( ); + } + catch( std::exception& err ) + { + std::cerr << "ERROR: " << err.what( ) << std::endl; + return( 1 ); + + } // yrt + return( 0 ); +} + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/Dijkstra.h b/libs/fpa/Base/Dijkstra.h index 2c49c06..d231c7c 100644 --- a/libs/fpa/Base/Dijkstra.h +++ b/libs/fpa/Base/Dijkstra.h @@ -7,6 +7,8 @@ #define __fpa__Base__Dijkstra__h__ #include +#include +#include namespace fpa { @@ -14,7 +16,7 @@ namespace fpa { /** */ - template< class _TFilter, class _TMarksInterface, class _TSeedsInterface > + template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST > class Dijkstra : public _TFilter, public _TMarksInterface, @@ -25,6 +27,7 @@ namespace fpa typedef _TFilter Superclass; typedef _TMarksInterface TMarksInterface; typedef _TSeedsInterface TSeedsInterface; + typedef _TMST TMST; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; @@ -34,7 +37,7 @@ namespace fpa typedef typename Superclass::TVertices TVertices; typedef itk::FunctionBase< TInputValue, TOutputValue > TIntensityFunctor; - typedef itk::FunctionBase< TVertex, TOutputValue > TVertexFunctor; + typedef fpa::Base::Functors::VertexParentBase< TVertex, TOutputValue > TVertexFunctor; protected: struct _TNode @@ -42,15 +45,17 @@ namespace fpa TVertex Vertex; TVertex Parent; TOutputValue Cost; - _TNode( const TVertex& v, const TVertex& p ) + unsigned long FrontId; + _TNode( const TVertex& v, const TVertex& p, const unsigned long& fId ) { this->Vertex = v; this->Parent = p; + this->FrontId = fId; this->Cost = TOutputValue( 0 ); } bool operator<( const _TNode& b ) const { - return( this->Cost < b.Cost ); + return( b.Cost < this->Cost ); } }; @@ -58,6 +63,9 @@ namespace fpa itkTypeMacro( Dijkstra, TFilter ); public: + TMST* GetMinimumSpanningTree( ); + const TMST* GetMinimumSpanningTree( ) const; + const TIntensityFunctor* GetIntensityFunctor( ) const; const TVertexFunctor* GetVertexFunctor( ) const; @@ -77,6 +85,7 @@ namespace fpa protected: typename TIntensityFunctor::Pointer m_IntensityFunctor; typename TVertexFunctor::Pointer m_VertexFunctor; + unsigned long m_MSTIndex; }; } // ecapseman diff --git a/libs/fpa/Base/Dijkstra.hxx b/libs/fpa/Base/Dijkstra.hxx index 9e8e529..c5563c9 100644 --- a/libs/fpa/Base/Dijkstra.hxx +++ b/libs/fpa/Base/Dijkstra.hxx @@ -7,30 +7,58 @@ #define __fpa__Base__Dijkstra__hxx__ // ------------------------------------------------------------------------- -template< class _TFilter, class _TMarksInterface, class _TSeedsInterface > +template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST > +typename +fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >:: +TMST* fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >:: +GetMinimumSpanningTree( ) +{ + return( + dynamic_cast< TMST* >( + this->itk::ProcessObject::GetOutput( this->m_MSTIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST > +const typename +fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >:: +TMST* fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >:: +GetMinimumSpanningTree( ) const +{ + return( + dynamic_cast< const TMST* >( + this->itk::ProcessObject::GetOutput( this->m_MSTIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST > const typename -fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >:: +fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >:: TIntensityFunctor* -fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >:: +fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >:: GetIntensityFunctor( ) const { return( this->m_IntensityFunctor ); } // ------------------------------------------------------------------------- -template< class _TFilter, class _TMarksInterface, class _TSeedsInterface > +template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST > const typename -fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >:: +fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >:: TVertexFunctor* -fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >:: +fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >:: GetVertexFunctor( ) const { return( this->m_VertexFunctor ); } // ------------------------------------------------------------------------- -template< class _TFilter, class _TMarksInterface, class _TSeedsInterface > -void fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >:: +template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST > +void fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >:: SetFunctor( TIntensityFunctor* functor ) { if( this->m_IntensityFunctor.GetPointer( ) != functor ) @@ -42,8 +70,8 @@ SetFunctor( TIntensityFunctor* functor ) } // ------------------------------------------------------------------------- -template< class _TFilter, class _TMarksInterface, class _TSeedsInterface > -void fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >:: +template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST > +void fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >:: SetFunctor( TVertexFunctor* functor ) { if( this->m_VertexFunctor.GetPointer( ) != functor ) @@ -55,27 +83,96 @@ SetFunctor( TVertexFunctor* functor ) } // ------------------------------------------------------------------------- -template< class _TFilter, class _TMarksInterface, class _TSeedsInterface > -fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >:: +template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST > +fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >:: Dijkstra( ) : Superclass( ), _TMarksInterface( this ), _TSeedsInterface( this ) { + this->m_MSTIndex = this->GetNumberOfRequiredOutputs( ); + this->SetNumberOfRequiredOutputs( this->m_MSTIndex + 1 ); + this->itk::ProcessObject::SetNthOutput( this->m_MSTIndex, TMST::New( ) ); } // ------------------------------------------------------------------------- -template< class _TFilter, class _TMarksInterface, class _TSeedsInterface > -fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >:: +template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST > +fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >:: ~Dijkstra( ) { } // ------------------------------------------------------------------------- -template< class _TFilter, class _TMarksInterface, class _TSeedsInterface > -void fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface >:: +template< class _TFilter, class _TMarksInterface, class _TSeedsInterface, class _TMST > +void fpa::Base::Dijkstra< _TFilter, _TMarksInterface, _TSeedsInterface, _TMST >:: GenerateData( ) { + // Init objects + this->_ConfigureOutputs( TOutputValue( 0 ) ); + this->_InitMarks( this->GetNumberOfSeeds( ) ); + TMST* mst = this->GetMinimumSpanningTree( ); + + // Init queue + std::vector< _TNode > q; + unsigned long frontId = 1; + for( TVertex seed: this->GetSeeds( ) ) + q.push_back( _TNode( seed, seed, frontId++ ) ); + + // Main loop + while( q.size( ) > 0 ) + { + // Get next candidate + std::pop_heap( q.begin( ), q.end( ) ); + _TNode node = q.back( ); + q.pop_back( ); + if( this->_IsMarked( node.Vertex ) ) + continue; + this->_Mark( node.Vertex, node.FrontId ); + + // Ok, pixel lays inside region + this->_SetOutputValue( node.Vertex, node.Cost ); + mst->SetParent( node.Vertex, node.Parent ); + + // Add neighborhood + TVertices neighbors = this->_GetNeighbors( node.Vertex ); + for( TVertex neigh: neighbors ) + { + if( this->_IsMarked( neigh ) ) + { + // Invoke stop at collisions + unsigned long nColl = this->_Collisions( node.Vertex, neigh ); + if( + this->StopAtOneFront( ) && + this->GetNumberOfSeeds( ) > 1 && + nColl == 1 + ) + q.clear( ); + } + else + { + // Compute new cost + TOutputValue ncost = + this->m_VertexFunctor->Evaluate( neigh, node.Vertex ); + if( this->m_IntensityFunctor.IsNotNull( ) ) + ncost = this->m_IntensityFunctor->Evaluate( ncost ); + + // This algorithm only supports positive values + if( ncost >= TOutputValue( 0 ) ) + { + // Insert new node + _TNode nn( neigh, node.Vertex, node.FrontId ); + nn.Cost = node.Cost + ncost; + q.push_back( nn ); + std::push_heap( q.begin( ), q.end( ) ); + + } // fi + + } // fi + + } // rof + + } // elihw + this->_FreeMarks( ); } #endif // __fpa__Base__Dijkstra__hxx__ diff --git a/libs/fpa/Base/Functors/InvertValue.h b/libs/fpa/Base/Functors/InvertValue.h new file mode 100644 index 0000000..9da3221 --- /dev/null +++ b/libs/fpa/Base/Functors/InvertValue.h @@ -0,0 +1,80 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Base__Functors__InvertValue__h__ +#define __fpa__Base__Functors__InvertValue__h__ + +#include +#include + +namespace fpa +{ + namespace Base + { + namespace Functors + { + /** + */ + template< class _TInputValue, class _TOutputValue > + class InvertValue + : public itk::FunctionBase< _TInputValue, _TOutputValue > + { + public: + typedef InvertValue Self; + typedef itk::FunctionBase< _TInputValue, _TOutputValue > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TInputValue TInputValue; + typedef _TOutputValue TOutputValue; + + public: + itkNewMacro( Self ); + itkTypeMacro( InvertValue, itk::FunctionBase ); + + itkGetConstMacro( Alpha, double ); + itkGetConstMacro( Beta, double ); + itkSetMacro( Alpha, double ); + itkSetMacro( Beta, double ); + + public: + virtual TOutputValue Evaluate( const TInputValue& a ) const override + { + double d = this->m_Alpha; + d += std::pow( double( a ), this->m_Beta ); + double x = -1; + if( std::fabs( d ) > double( 0 ) ) + x = + double( 1 ) / + ( this->m_Alpha + std::pow( double( a ), this->m_Beta ) ); + return( TOutputValue( x ) ); + } + + protected: + InvertValue( ) + : Superclass( ), + m_Alpha( double( 1 ) ), + m_Beta( double( 1 ) ) + { } + virtual ~InvertValue( ) { } + + private: + InvertValue( const Self& other ); + Self& operator=( const Self& other ); + + protected: + double m_Alpha; + double m_Beta; + }; + + } // ecapseman + + } // ecapseman + +} // ecapseman + +#endif // __fpa__Base__Functors__InvertValue__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/Functors/VertexParentBase.h b/libs/fpa/Base/Functors/VertexParentBase.h new file mode 100644 index 0000000..2725063 --- /dev/null +++ b/libs/fpa/Base/Functors/VertexParentBase.h @@ -0,0 +1,58 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Base__Functors__VertexParentBase__h__ +#define __fpa__Base__Functors__VertexParentBase__h__ + +#include +#include + +namespace fpa +{ + namespace Base + { + namespace Functors + { + /** + */ + template< class _TVertex, class _TOutputValue > + class VertexParentBase + : public itk::Object + { + public: + typedef VertexParentBase Self; + typedef itk::Object Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TVertex TVertex; + typedef _TOutputValue TOutputValue; + + public: + itkTypeMacro( VertexParentBase, TFilter ); + + public: + virtual TOutputValue Evaluate( const TVertex& v, const TVertex& p ) const = 0; + + protected: + VertexParentBase( ) + : Superclass( ) + { } + virtual ~VertexParentBase( ) { } + + private: + VertexParentBase( const Self& other ); + Self& operator=( const Self& other ); + }; + + } // ecapseman + + } // ecapseman + +} // ecapseman + +#endif // __fpa__Base__Functors__VertexParentBase__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/MinimumSpanningTree.h b/libs/fpa/Base/MinimumSpanningTree.h new file mode 100644 index 0000000..35cb6e4 --- /dev/null +++ b/libs/fpa/Base/MinimumSpanningTree.h @@ -0,0 +1,77 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Base__MinimumSpanningTree__h__ +#define __fpa__Base__MinimumSpanningTree__h__ + +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< class _TVertex, class _Superclass > + class MinimumSpanningTree + : public _Superclass + { + public: + typedef MinimumSpanningTree Self; + typedef _Superclass Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TVertex TVertex; + typedef std::pair< TVertex, bool > TCollision; + typedef std::vector< TCollision > TCollisionsRow; + typedef std::vector< TCollisionsRow > TCollisions; + typedef std::vector< TVertex > TVertices; + + protected: + typedef std::vector< unsigned long > _TRow; + typedef std::vector< _TRow > _TMatrix; + + public: + itkTypeMacro( fpa::Base::MinimumSpanningTree, _Superclass ); + + public: + const TCollisions& GetCollisions( ) const; + void SetCollisions( const TCollisions& collisions ); + + void ClearSeeds( ); + void AddSeed( const TVertex& seed ); + + virtual TVertex GetParent( const TVertex& v ) const = 0; + virtual void SetParent( const TVertex& v, const TVertex& p ) = 0; + + virtual TVertices GetPath( const TVertex& a ) const; + virtual TVertices GetPath( const TVertex& a, const TVertex& b ) const; + + protected: + MinimumSpanningTree( ); + virtual ~MinimumSpanningTree( ); + + private: + MinimumSpanningTree( const Self& other ); + Self& operator=( const Self& other ); + + protected: + TCollisions m_Collisions; + _TMatrix m_FrontPaths; + std::vector< TVertex > m_Seeds; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Base__MinimumSpanningTree__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/MinimumSpanningTree.hxx b/libs/fpa/Base/MinimumSpanningTree.hxx new file mode 100644 index 0000000..79ddebe --- /dev/null +++ b/libs/fpa/Base/MinimumSpanningTree.hxx @@ -0,0 +1,239 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Base__MinimumSpanningTree__hxx__ +#define __fpa__Base__MinimumSpanningTree__hxx__ + +// ------------------------------------------------------------------------- +template< class _TVertex, class _Superclass > +const typename fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >:: +TCollisions& fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >:: +GetCollisions( ) const +{ + return( this->m_Collisions ); +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _Superclass > +void fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >:: +SetCollisions( const TCollisions& collisions ) +{ + static const unsigned long _inf = + std::numeric_limits< unsigned long >::max( ); + if( this->m_Collisions == collisions ) + return; + + this->m_Collisions = collisions; + + // Prepare a front graph + unsigned long N = this->m_Collisions.size( ); + _TMatrix dist( N, _TRow( N, _inf ) ); + this->m_FrontPaths = dist; + for( unsigned long i = 0; i < N; ++i ) + { + for( unsigned long j = 0; j < N; ++j ) + { + if( this->m_Collisions[ i ][ j ].second ) + { + dist[ i ][ j ] = 1; + dist[ j ][ i ] = 1; + this->m_FrontPaths[ i ][ j ] = j; + this->m_FrontPaths[ j ][ i ] = i; + + } // fi + + } // rof + dist[ i ][ i ] = 0; + this->m_FrontPaths[ i ][ i ] = i; + + } // rof + + // Use Floyd-Warshall to compute all possible paths between fronts + for( unsigned long k = 0; k < N; ++k ) + { + for( unsigned long i = 0; i < N; ++i ) + { + for( unsigned long j = 0; j < N; ++j ) + { + // WARNING: you don't want a numeric overflow!!! + unsigned long dik = dist[ i ][ k ]; + unsigned long dkj = dist[ k ][ j ]; + unsigned long sum = _inf; + if( dik < _inf && dkj < _inf ) + sum = dik + dkj; + + // Ok, continue Floyd-Warshall + if( sum < dist[ i ][ j ] ) + { + dist[ i ][ j ] = sum; + this->m_FrontPaths[ i ][ j ] = this->m_FrontPaths[ i ][ k ]; + + } // fi + + } // rof + + } // rof + + } // rof + this->Modified( ); +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _Superclass > +void fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >:: +ClearSeeds( ) +{ + this->m_Seeds.clear( ); + if( this->m_DataObject != NULL ) + this->m_DataObject->Modified( ); +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _Superclass > +void fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >:: +AddSeed( const _TVertex& seed ) +{ + this->m_Seeds.push_back( seed ); + if( this->m_DataObject != NULL ) + this->m_DataObject->Modified( ); +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _Superclass > +typename fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >:: +TVertices fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >:: +GetPath( const _TVertex& a ) const +{ + TVertices vertices; + _TVertex it = a; + _TVertex p = this->GetParent( it ); + while( it != p ) + { + vertices.push_back( it ); + it = p; + p = this->GetParent( it ); + + } // elihw + vertices.push_back( it ); + return( vertices ); +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _Superclass > +typename fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >:: +TVertices fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >:: +GetPath( const _TVertex& a, const _TVertex& b ) const +{ + static const unsigned long _inf = + std::numeric_limits< unsigned long >::max( ); + + TVertices vertices; + TVertices pa = this->GetPath( a ); + TVertices pb = this->GetPath( b ); + if( pa.size( ) > 0 && pb.size( ) > 0 ) + { + // Find front identifiers + unsigned long ia = _inf, ib = _inf; + unsigned long N = this->m_Seeds.size( ); + for( unsigned long i = 0; i < N; ++i ) + { + if( this->m_Seeds[ i ] == pa[ pa.size( ) - 1 ] ) + ia = i; + if( this->m_Seeds[ i ] == pb[ pb.size( ) - 1 ] ) + ib = i; + + } // rof + + // Check if there is a front-jump between given seeds + if( ia != ib ) + { + // Compute front path + std::vector< long > fpath; + fpath.push_back( ia ); + while( ia != ib ) + { + ia = this->m_FrontPaths[ ia ][ ib ]; + fpath.push_back( ia ); + + } // elihw + + // Continue only if both fronts are connected + unsigned int N = fpath.size( ); + if( N > 0 ) + { + // First path: from start vertex to first collision + vertices = this->GetPath( + a, this->m_Collisions[ fpath[ 0 ] ][ fpath[ 1 ] ].first + ); + + // Intermediary paths + for( unsigned int i = 1; i < N - 1; ++i ) + { + TVertices ipath = + this->GetPath( + this->m_Collisions[ fpath[ i ] ][ fpath[ i - 1 ] ].first, + this->m_Collisions[ fpath[ i ] ][ fpath[ i + 1 ] ].first + ); + for( long id = 0; id < ipath.size( ); ++id ) + vertices.push_back( ipath[ id ] ); + + } // rof + + // Final path: from last collision to end point + TVertices lpath = + this->GetPath( + this->m_Collisions[ fpath[ N - 1 ] ][ fpath[ N - 2 ] ].first, b + ); + for( long id = 0; id < lpath.size( ); ++id ) + vertices.push_back( lpath[ id ] ); + + } // fi + } + else + { + // Ignore common part: find common ancestor + long aIt = pa.size( ) - 1; + long bIt = pb.size( ) - 1; + bool cont = true; + while( aIt >= 0 && bIt >= 0 && cont ) + { + cont = ( pa[ aIt ] == pb[ bIt ] ); + aIt--; + bIt--; + + } // elihw + aIt++; + bIt++; + + // Glue both parts + for( long cIt = 0; cIt <= aIt; ++cIt ) + vertices.push_back( pa[ cIt ] ); + for( ; bIt >= 0; --bIt ) + vertices.push_back( pb[ bIt ] ); + + } // fi + + } // fi + return( vertices ); +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _Superclass > +fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >:: +MinimumSpanningTree( ) + : Superclass( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _Superclass > +fpa::Base::MinimumSpanningTree< _TVertex, _Superclass >:: +~MinimumSpanningTree( ) +{ +} + +#endif // __fpa__Base__MinimumSpanningTree__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/Dijkstra.h b/libs/fpa/Image/Dijkstra.h new file mode 100644 index 0000000..8516d7b --- /dev/null +++ b/libs/fpa/Image/Dijkstra.h @@ -0,0 +1,91 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Image__Dijkstra__h__ +#define __fpa__Image__Dijkstra__h__ + +#include +#include +#include +#include +#include +#include + +namespace fpa +{ + namespace Image + { + /** + */ + template< class _TInputImage, class _TOutputImage > + class Dijkstra + : public fpa::Base::Dijkstra< fpa::Image::Filter< _TInputImage, _TOutputImage >, fpa::Image::MarksInterface< _TInputImage::ImageDimension >, fpa::Base::SeedsInterface< typename _TInputImage::IndexType, typename _TInputImage::IndexType::LexicographicCompare >, fpa::Image::MinimumSpanningTree< _TInputImage::ImageDimension > > + { + public: + // Interfaces + typedef fpa::Image::Filter< _TInputImage, _TOutputImage > TFilter; + typedef fpa::Image::MarksInterface< _TInputImage::ImageDimension > TMarksInterface; + typedef fpa::Base::SeedsInterface< typename _TInputImage::IndexType, typename _TInputImage::IndexType::LexicographicCompare > TSeedsInterface; + typedef fpa::Image::MinimumSpanningTree< _TInputImage::ImageDimension > TMST; + + // Smart pointers + typedef Dijkstra Self; + typedef fpa::Base::Dijkstra< TFilter, TMarksInterface, TSeedsInterface, TMST > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef typename TFilter::TInputImage TInputImage; + typedef typename TFilter::TOutputValue TOutputValue; + typedef typename TFilter::TVertex TVertex; + + public: + itkNewMacro( Self ); + itkTypeMacro( fpa::Image::Dijkstra, fpa::Base::Dijkstra ); + + protected: + Dijkstra( ) : Superclass( ) { } + virtual ~Dijkstra( ) { } + + virtual void _ConfigureOutputs( const TOutputValue& init_value ) override + { + this->Superclass::_ConfigureOutputs( init_value ); + + typename TVertex::OffsetType o; + o.Fill( 0 ); + const TInputImage* input = this->GetInput( ); + TMST* mst = this->GetMinimumSpanningTree( ); + mst->CopyInformation( input ); + mst->SetBufferedRegion( input->GetRequestedRegion( ) ); + mst->Allocate( ); + mst->FillBuffer( o ); + } + + virtual void GenerateData( ) override + { + // Configure functors with input image + typedef typename TFilter::TInputImage _TInputIage; + typedef typename TFilter::TOutputValue _TOutputValue; + typedef fpa::Image::Functors::VertexParentBase< _TInputImage, _TOutputValue > _TVFunc; + _TVFunc* vfunc = + dynamic_cast< _TVFunc* >( this->m_VertexFunctor.GetPointer( ) ); + if( vfunc != NULL ) + vfunc->SetImage( this->GetInput( ) ); + + // Ok, continue + this->Superclass::GenerateData( ); + } + + private: + Dijkstra( const Self& other ); + Self& operator=( const Self& other ); + }; + + } // ecapseman + +} // ecapseman + +#endif // __fpa__Image__Dijkstra__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/Functors/VertexIdentity.h b/libs/fpa/Image/Functors/VertexIdentity.h new file mode 100644 index 0000000..716c2e7 --- /dev/null +++ b/libs/fpa/Image/Functors/VertexIdentity.h @@ -0,0 +1,67 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Image__Functors__VertexIdentity__h__ +#define __fpa__Image__Functors__VertexIdentity__h__ + +#include + +namespace fpa +{ + namespace Image + { + namespace Functors + { + /** + */ + template< class _TInputImage, class _TOutputValue > + class VertexIdentity + : public fpa::Image::Functors::VertexParentBase< _TInputImage, _TOutputValue > + { + public: + typedef _TInputImage TInputImage; + typedef _TOutputValue TOutputValue; + typedef VertexIdentity Self; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + typedef fpa::Image::Functors::VertexParentBase< TInputImage, TOutputValue > Superclass; + + typedef typename Superclass::TVertex TVertex; + + public: + itkNewMacro( Self ); + itkTypeMacro( + fpa::Image::Functors::VertexIdentity, + fpa::Image::Functors::VertexParentBase + ); + + public: + virtual TOutputValue Evaluate( + const TVertex& a, const TVertex& p + ) const override + { + return( TOutputValue( this->m_Image->GetPixel( a ) ) ); + } + + protected: + VertexIdentity( ) + : Superclass( ) + { } + virtual ~VertexIdentity( ) { } + + private: + VertexIdentity( const Self& other ); + Self& operator=( const Self& other ); + }; + + } // ecapseman + + } // ecapseman + +} // ecapseman + +#endif // __fpa__Image__Functors__VertexIdentity__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/Functors/VertexParentBase.h b/libs/fpa/Image/Functors/VertexParentBase.h new file mode 100644 index 0000000..20b76e0 --- /dev/null +++ b/libs/fpa/Image/Functors/VertexParentBase.h @@ -0,0 +1,63 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Image__Functors__VertexParentBase__h__ +#define __fpa__Image__Functors__VertexParentBase__h__ + +#include + +namespace fpa +{ + namespace Image + { + namespace Functors + { + /** + */ + template< class _TInputImage, class _TOutputValue > + class VertexParentBase + : public fpa::Base::Functors::VertexParentBase< typename _TInputImage::IndexType, _TOutputValue > + { + public: + typedef _TInputImage TInputImage; + typedef _TOutputValue TOutputValue; + typedef typename TInputImage::IndexType TVertex; + typedef VertexParentBase Self; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + typedef fpa::Base::Functors::VertexParentBase< TVertex, TOutputValue > Superclass; + + public: + itkTypeMacro( + fpa::Image::Functors::VertexParentBase, + fpa::Base::Functors::VertexParentBase + ); + + itkGetConstObjectMacro( Image, TInputImage ); + itkSetConstObjectMacro( Image, TInputImage ); + + protected: + VertexParentBase( ) + : Superclass( ) + { } + virtual ~VertexParentBase( ) { } + + private: + VertexParentBase( const Self& other ); + Self& operator=( const Self& other ); + + protected: + typename TInputImage::ConstPointer m_Image; + }; + + } // ecapseman + + } // ecapseman + +} // ecapseman + +#endif // __fpa__Image__Functors__VertexParentBase__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/MinimumSpanningTree.h b/libs/fpa/Image/MinimumSpanningTree.h new file mode 100644 index 0000000..a93157c --- /dev/null +++ b/libs/fpa/Image/MinimumSpanningTree.h @@ -0,0 +1,71 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Image__MinimumSpanningTree__h__ +#define __fpa__Image__MinimumSpanningTree__h__ + +#include +#include + +namespace fpa +{ + namespace Image + { + /** + */ + template< unsigned int _VDim > + class MinimumSpanningTree + : public fpa::Base::MinimumSpanningTree< itk::Index< _VDim >, itk::Image< itk::Offset< _VDim >, _VDim > > + { + public: + typedef itk::Index< _VDim > TVertex; + typedef itk::Image< itk::Offset< _VDim >, _VDim > TBaseImage; + + typedef MinimumSpanningTree Self; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + typedef fpa::Base::MinimumSpanningTree< TVertex, TBaseImage > Superclass; + + typedef typename Superclass::TCollision TCollision; + typedef typename Superclass::TCollisionsRow TCollisionsRow; + typedef typename Superclass::TCollisions TCollisions; + typedef typename Superclass::TVertices TVertices; + + public: + itkNewMacro( Self ); + itkTypeMacro( + fpa::Image::MinimumSpanningTree, + fpa::Base::MinimumSpanningTree + ); + + public: + virtual TVertex GetParent( const TVertex& v ) const override + { + return( v + this->GetPixel( v ) ); + } + virtual void SetParent( const TVertex& v, const TVertex& p ) override + { + this->SetPixel( v, p - v ); + } + + protected: + MinimumSpanningTree( ) + : Superclass( ) + { } + virtual ~MinimumSpanningTree( ) + { } + + private: + MinimumSpanningTree( const Self& other ); + Self& operator=( const Self& other ); + }; + + } // ecapseman + +} // ecapseman + +#endif // __fpa__Image__MinimumSpanningTree__h__ + +// eof - $RCSfile$ -- 2.47.1