#ifndef __fpa__Image__EndPointsFilter__hxx__
#define __fpa__Image__EndPointsFilter__hxx__
+#include <functional>
#include <map>
#include <queue>
#include <itkImageRegionConstIteratorWithIndex.h>
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;
- _TDistMapIt dIt(
- this->m_DistanceMap, this->m_DistanceMap->GetRequestedRegion( )
- );
- _TCostMapIt cIt(
- this->m_CostMap, this->m_CostMap->GetRequestedRegion( )
- );
+ _TDistMapIt dIt( this->m_DistanceMap, region );
+ _TCostMapIt cIt( this->m_CostMap, region );
dIt.GoToBegin( );
cIt.GoToBegin( );
for( ; !dIt.IsAtEnd( ) && !cIt.IsAtEnd( ); ++dIt, ++cIt )
double d = double( dIt.Get( ) );
if( d > 0 )
{
- double v = double( cIt.Get( ) ) * d;
+ double v = double( cIt.Get( ) ) / ( d * d );
queue.insert( _TQueueValue( v, dIt.GetIndex( ) ) );
} // fi
} // rof
- TIndices marks;
+ // BFS from maximum queue
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
+ m |= 0x01;
+ marks->SetPixel( n.second, m );
+ this->m_EndPoints.insert( n.second );
+
+ // Get path
+ typename TMST::TPath::Pointer path;
+ this->m_MST->GetPath( path, n.second );
+ for( unsigned long i = 0; i < path->GetSize( ); ++i )
{
- std::cout << queue.size( ) << " " << n.first << std::endl;
- 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 )
+ TIndex idx = path->GetVertex( path->GetSize( ) - 1 - i );
+ typename _TCostMap::PointType cnt;
+ this->m_CostMap->TransformIndexToPhysicalPoint( idx, cnt );
+ double d = double( this->m_DistanceMap->GetPixel( idx ) );
+ d *= double( 2 );
+
+ // Mark sphere
+ std::queue< TIndex > q;
+ q.push( idx );
+ while( q.size( ) > 0 )
{
- double d = double( this->m_DistanceMap->GetPixel( *pIt ) );
- d = std::sqrt( std::fabs( d ) );
- typename _TCostMap::PointType center;
- this->m_CostMap->TransformIndexToPhysicalPoint( *pIt, center );
-
- std::queue< TIndex > q;
- TIndices m;
- q.push( *pIt );
- 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 )
{
- 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 )
{
- for( int y = -1; y <= 1; y += 2 )
+ TIndex w = v;
+ w[ x ] += y;
+ if( region.IsInside( w ) )
{
- 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 );
+ typename _TCostMap::PointType p;
+ this->m_CostMap->TransformIndexToPhysicalPoint( w, p );
+ if( cnt.EuclideanDistanceTo( p ) <= d )
+ q.push( w );
- } // fi
-
- } // rof
+ } // fi
} // rof
- } // elihw
+ } // rof
- } // rof
+ } // elihw
- } // fi
+ } // rof
} // elihw
}