From 3e5846a571cd387c59348a7748a21bc3bf620d95 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Leonardo=20Fl=C3=B3rez-Valencia?= Date: Mon, 10 Oct 2016 14:48:23 -0500 Subject: [PATCH] ... --- lib/Instances/fpa_DataObjects.i | 4 + lib/fpa/Base/MinimumSpanningTree.hxx | 131 ++++++++++---------- lib/fpa/Image/Algorithm.hxx | 11 +- lib/fpa/Image/Dijkstra.hxx | 8 ++ lib/fpa/Image/EndPointsFilter.hxx | 121 ++++++++---------- plugins/Plugins/ImageDijkstra.cxx | 9 +- plugins/Plugins/SimpleImageDijkstraCost.cxx | 3 +- plugins/Plugins/SimpleImageNeighborhood.cxx | 3 +- 8 files changed, 150 insertions(+), 140 deletions(-) diff --git a/lib/Instances/fpa_DataObjects.i b/lib/Instances/fpa_DataObjects.i index 231d59f..62280ab 100644 --- a/lib/Instances/fpa_DataObjects.i +++ b/lib/Instances/fpa_DataObjects.i @@ -1,8 +1,12 @@ i cpPlugins_Paths.h +t fpa/Image/Functors/SimpleNeighborhood +t fpa/Image/Functors/SimpleDijkstraCost t fpa/Image/MinimumSpanningTree t fpa/Base/MinimumSpanningTree +c fpa::Image::Functors::SimpleNeighborhood< itk::ImageBase< #process_dims# > > +c fpa::Image::Functors::SimpleDijkstraCost< itk::Image< #pixels#, #process_dims# >, #reals# > c fpa::Base::MinimumSpanningTree< itk::Index< #process_dims# >, cpExtensions::DataStructures::PolyLineParametricPath< #process_dims# >, itk::Image< itk::Offset< #process_dims# >, #process_dims# > > c fpa::Image::MinimumSpanningTree< #process_dims# > diff --git a/lib/fpa/Base/MinimumSpanningTree.hxx b/lib/fpa/Base/MinimumSpanningTree.hxx index 4bf0135..c85c283 100644 --- a/lib/fpa/Base/MinimumSpanningTree.hxx +++ b/lib/fpa/Base/MinimumSpanningTree.hxx @@ -114,7 +114,8 @@ GetPath( typename _TPath::Pointer& path, const _TVertex& a ) const } // elihw vertices.push_back( it ); - path = _TPath::New( ); + if( path.IsNull( ) ) + path = _TPath::New( ); for( auto v = vertices.begin( ); v != vertices.end( ); ++v ) path->AddVertex( *v ); } @@ -126,97 +127,101 @@ GetPath( typename _TPath::Pointer& path, const _TVertex& a, const _TVertex& b ) const { -/* static const unsigned long _inf = std::numeric_limits< unsigned long >::max( ); - _TPath path; - _TPath pa = this->GetPath( a ); - _TPath pb = this->GetPath( b ); - if( pa.size( ) > 0 && pb.size( ) > 0 ) + if( path.IsNull( ) ) + path = _TPath::New( ); + typename _TPath::Pointer pa, pb; + this->GetPath( pa, a ); + this->GetPath( pb, b ); + if( pa->GetSize( ) > 0 && pb->GetSize( ) > 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.front( ) ) + if( this->m_Seeds[ i ] == pa->GetVertex( 0 ) ) ia = i; - if( this->m_Seeds[ i ] == pb.front( ) ) + if( this->m_Seeds[ i ] == pb->GetVertex( 0 ) ) ib = i; } // rof + // Check if there is a front-jump between given seeds if( ia != ib ) { - // Use this->m_FrontPaths from Floyd-Warshall - if( this->m_FrontPaths[ ia ][ ib ] < _inf ) - { - // 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 - path = 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 - ); - path.insert( path.end( ), ipath.begin( ), ipath.end( ) ); - - } // rof - - // Final path: from last collision to end point - TVertices lpath = - this->GetPath( - this->m_Collisions[ fpath[ N - 1 ] ][ fpath[ N - 2 ] ].first, b - ); - path.insert( path.end( ), lpath.begin( ), lpath.end( ) ); - - } // fi - - } // fi + // Compute front path + /* TODO + 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 + this->GetPath( + path, a, + this->m_Collisions[ fpath[ 0 ] ][ fpath[ 1 ] ].first + ); + + // Intermediary paths + for( unsigned int i = 1; i < N - 1; ++i ) + { + TPathTVertices ipath = + this->GetPath( + this->m_Collisions[ fpath[ i ] ][ fpath[ i - 1 ] ].first, + this->m_Collisions[ fpath[ i ] ][ fpath[ i + 1 ] ].first + ); + path.insert( path.end( ), ipath.begin( ), ipath.end( ) ); + + } // rof + + // Final path: from last collision to end point + TVertices lpath = + this->GetPath( + this->m_Collisions[ fpath[ N - 1 ] ][ fpath[ N - 2 ] ].first, b + ); + path.insert( path.end( ), lpath.begin( ), lpath.end( ) ); + + } // fi + + } // fi + */ } else { // Ignore common part: find common ancestor - auto aIt = pa.begin( ); - auto bIt = pb.begin( ); - while( *aIt == *bIt && aIt != pa.end( ) && bIt != pb.end( ) ) + long aIt = pa->GetSize( ) - 1; + long bIt = pb->GetSize( ) - 1; + bool cont = true; + while( aIt >= 0 && bIt >= 0 && cont ) { - ++aIt; - ++bIt; + cont = ( pa->GetVertex( aIt ) == pb->GetVertex( bIt ) ); + aIt--; + bIt--; } // elihw + aIt++; + bIt++; // Glue both parts - for( --aIt; aIt != pa.end( ); ++aIt ) - path.push_front( *aIt ); - for( ; bIt != pb.end( ); ++bIt ) - path.push_back( *bIt ); + for( long cIt = 0; cIt <= aIt; ++cIt ) + path->AddVertex( pa->GetVertex( cIt ) ); + for( ; bIt >= 0; --bIt ) + path->AddVertex( pb->GetVertex( bIt ) ); } // fi } // fi - return( path ); -*/ } // ------------------------------------------------------------------------- diff --git a/lib/fpa/Image/Algorithm.hxx b/lib/fpa/Image/Algorithm.hxx index 9eb2640..ed73421 100644 --- a/lib/fpa/Image/Algorithm.hxx +++ b/lib/fpa/Image/Algorithm.hxx @@ -4,6 +4,7 @@ // Send Piotr's code to Anna #include +#include // ------------------------------------------------------------------------- template< class _TInputImage, class _TOutputImage > @@ -11,11 +12,17 @@ fpa::Image::Algorithm< _TInputImage, _TOutputImage >:: Algorithm( ) : Superclass( ) { + typedef itk::ImageBase< _TInputImage::ImageDimension > _TImageBase; + typedef fpa::Image::Functors::SimpleNeighborhood< _TImageBase > _TNeigh; typedef itk::Image< TFrontId, _TInputImage::ImageDimension > _TMarks; + this->m_MarksIdx = this->GetNumberOfRequiredOutputs( ); this->itk::ProcessObject::SetNumberOfRequiredOutputs( this->m_MarksIdx + 1 ); typename _TMarks::Pointer marks = _TMarks::New( ); this->SetNthOutput( this->m_MarksIdx, marks ); + + typename _TNeigh::Pointer neigh = _TNeigh::New( ); + this->SetNeighborhoodFunction( neigh ); } // ------------------------------------------------------------------------- @@ -34,7 +41,9 @@ _BeforeGenerateData( ) this->AllocateOutputs( ); TNeighborhoodFunction* neighFunc = - dynamic_cast< TNeighborhoodFunction* >( this->GetNeighborhoodFunction( ) ); + dynamic_cast< TNeighborhoodFunction* >( + this->GetNeighborhoodFunction( ) + ); if( neighFunc == NULL ) itkExceptionMacro( << "NeighborhoodFunction not well defined." ); neighFunc->SetImage( this->GetInput( ) ); diff --git a/lib/fpa/Image/Dijkstra.hxx b/lib/fpa/Image/Dijkstra.hxx index 2831159..bf6b5f0 100644 --- a/lib/fpa/Image/Dijkstra.hxx +++ b/lib/fpa/Image/Dijkstra.hxx @@ -1,12 +1,20 @@ #ifndef __fpa__Image__Dijkstra__hxx__ #define __fpa__Image__Dijkstra__hxx__ +#include + // ------------------------------------------------------------------------- template< class _TInputImage, class _TOutputImage > fpa::Image::Dijkstra< _TInputImage, _TOutputImage >:: Dijkstra( ) : Superclass( ) { + typedef + fpa::Image::Functors::SimpleDijkstraCost< _TInputImage, TOutput > + _TCost; + + typename _TCost::Pointer cost = _TCost::New( ); + this->SetCostFunction( cost ); } // ------------------------------------------------------------------------- diff --git a/lib/fpa/Image/EndPointsFilter.hxx b/lib/fpa/Image/EndPointsFilter.hxx index a4c170e..5600e86 100644 --- a/lib/fpa/Image/EndPointsFilter.hxx +++ b/lib/fpa/Image/EndPointsFilter.hxx @@ -33,10 +33,20 @@ Compute( ) typedef itk::ImageRegionConstIteratorWithIndex< _TCostMap > _TCostMapIt; typedef std::multimap< double, TIndex, std::greater< double > > _TQueue; typedef typename _TQueue::value_type _TQueueValue; + typedef itk::Image< unsigned char, _TCostMap::ImageDimension > _TMarks; // Some values typename _TDistanceMap::RegionType region = this->m_DistanceMap->GetRequestedRegion( ); + typename _TMarks::Pointer marks = _TMarks::New( ); + marks->SetLargestPossibleRegion( this->m_DistanceMap->GetLargestPossibleRegion( ) ); + marks->SetRequestedRegion( this->m_DistanceMap->GetRequestedRegion( ) ); + marks->SetBufferedRegion( this->m_DistanceMap->GetBufferedRegion( ) ); + marks->SetSpacing( this->m_DistanceMap->GetSpacing( ) ); + marks->SetOrigin( this->m_DistanceMap->GetOrigin( ) ); + marks->SetDirection( this->m_DistanceMap->GetDirection( ) ); + marks->Allocate( ); + marks->FillBuffer( 0 ); // Create queue _TQueue queue; @@ -57,18 +67,20 @@ Compute( ) } // rof // BFS from maximum queue - TIndices marks; while( queue.size( ) > 0 ) { // Get node auto nIt = queue.begin( ); auto n = *nIt; queue.erase( nIt ); - if( marks.find( n.second ) != marks.end( ) ) + + unsigned char m = marks->GetPixel( n.second ); + if( ( m & 0x01 ) == 0x01 ) continue; // Mark it - marks.insert( n.second ); + m |= 0x01; + marks->SetPixel( n.second, m ); this->m_EndPoints.insert( n.second ); // Get path @@ -76,78 +88,49 @@ Compute( ) this->m_MST->GetPath( path, n.second ); for( unsigned long i = 0; i < path->GetSize( ); ++i ) { - typename TMST::TPath::TContinuousIndex cidx = path->GetVertex( i ); + TIndex idx = path->GetVertex( path->GetSize( ) - 1 - i ); typename _TCostMap::PointType cnt; - this->m_CostMap->TransformContinuousIndexToPhysicalPoint( cidx, cnt ); - TIndex idx; - this->m_CostMap->TransformPhysicalPointToIndex( cnt, idx ); + this->m_CostMap->TransformIndexToPhysicalPoint( idx, cnt ); double d = double( this->m_DistanceMap->GetPixel( idx ) ); - - /* TODO - TIndex idx; - for( unsigned int d = 0; d < _TCostMap::ImageDimension; ++d ) - idx[ d ] = cidx[ d ]; - */ + d *= double( 2 ); + + // Mark sphere + std::queue< TIndex > q; + q.push( idx ); + while( q.size( ) > 0 ) + { + TIndex v = q.front( ); + q.pop( ); + unsigned char m = marks->GetPixel( v ); + if( ( m & 0x02 ) == 0x02 ) + continue; + m |= 0x03; + marks->SetPixel( v, m ); + + for( unsigned int x = 0; x < _TCostMap::ImageDimension; ++x ) + { + for( int y = -1; y <= 1; y += 2 ) + { + TIndex w = v; + w[ x ] += y; + if( region.IsInside( w ) ) + { + typename _TCostMap::PointType p; + this->m_CostMap->TransformIndexToPhysicalPoint( w, p ); + if( cnt.EuclideanDistanceTo( p ) <= d ) + q.push( w ); + + } // fi + + } // rof + + } // rof + + } // elihw } // rof - // TODO: temporary - queue.clear( ); - } // elihw - - /* TODO - while( queue.size( ) > 0 ) - { - marks.insert( n.second ); - this->m_EndPoints.insert( n.second ); - auto path = this->m_MST->GetPath( n.second ); - std::cout << path.size( ) << std::endl; - for( auto pIt = path.begin( ); pIt != path.end( ); ++pIt ) - { - double d = double( this->m_DistanceMap->GetPixel( *pIt ) ); - d = std::sqrt( std::fabs( d ) ); - typename _TCostMap::PointType center; - - - std::queue< TIndex > q; - TIndices m; - q.push( *pIt ); - while( q.size( ) > 0 ) - { - TIndex idx = q.front( ); - q.pop( ); - if( m.find( idx ) != m.end( ) ) - continue; - m.insert( idx ); - marks.insert( idx ); - for( unsigned int x = 0; x < _TCostMap::ImageDimension; ++x ) - { - for( int y = -1; y <= 1; y += 2 ) - { - TIndex idx2 = idx; - idx2[ x ] += y; - typename _TCostMap::PointType c; - this->m_CostMap->TransformIndexToPhysicalPoint( idx2, c ); - if( this->m_CostMap->GetRequestedRegion( ).IsInside( idx2 ) ) - { - if( center.EuclideanDistanceTo( c ) <= ( d * 1.5 ) ) - q.push( idx2 ); - - } // fi - - } // rof - - } // rof - - } // elihw - - } // rof - - } // fi - - } // elihw - */ } // ------------------------------------------------------------------------- diff --git a/plugins/Plugins/ImageDijkstra.cxx b/plugins/Plugins/ImageDijkstra.cxx index c5d34a4..d00636e 100644 --- a/plugins/Plugins/ImageDijkstra.cxx +++ b/plugins/Plugins/ImageDijkstra.cxx @@ -71,9 +71,12 @@ _GD1( _TInputImage* image ) // Configure filter auto filter = this->_CreateITK< _TFilter >( ); filter->SetInput( image ); - filter->SetNeighborhoodFunction( neig ); - filter->SetCostFunction( cost ); - filter->SetCostConversionFunction( conv ); + if( neig != NULL ) + filter->SetNeighborhoodFunction( neig ); + if( cost != NULL ) + filter->SetCostFunction( cost ); + if( conv != NULL ) + filter->SetCostConversionFunction( conv ); filter->SetStopAtOneFront( this->m_Parameters.GetBool( "StopAtOneFront" ) ); // Assign seeds diff --git a/plugins/Plugins/SimpleImageDijkstraCost.cxx b/plugins/Plugins/SimpleImageDijkstraCost.cxx index 6557a7f..c87b912 100644 --- a/plugins/Plugins/SimpleImageDijkstraCost.cxx +++ b/plugins/Plugins/SimpleImageDijkstraCost.cxx @@ -1,8 +1,7 @@ #include #include -#include -#include +#include // ------------------------------------------------------------------------- fpaPlugins::SimpleImageDijkstraCost:: diff --git a/plugins/Plugins/SimpleImageNeighborhood.cxx b/plugins/Plugins/SimpleImageNeighborhood.cxx index 0fc9837..11f2658 100644 --- a/plugins/Plugins/SimpleImageNeighborhood.cxx +++ b/plugins/Plugins/SimpleImageNeighborhood.cxx @@ -1,8 +1,7 @@ #include #include -#include -#include +#include // ------------------------------------------------------------------------- fpaPlugins::SimpleImageNeighborhood:: -- 2.47.1