#define __fpa__Common__OriginalRandomWalker__hxx__
#include <cmath>
-#include <map>
#include <itkImageRegionConstIteratorWithIndex.h>
+#include <itkImageRegionIteratorWithIndex.h>
#include <Eigen/Sparse>
// -------------------------------------------------------------------------
m_Epsilon( TScalar( 1e-5 ) ),
m_NormalizeWeights( true )
{
- fpaFilterOptionalInputConfigureMacro( InputLabels, TLabels );
+ fpaFilterInputConfigureMacro( InputLabels, TLabels );
fpaFilterOutputConfigureMacro( OutputProbabilities, TScalarImage );
}
void fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
GenerateData( )
{
- typedef Eigen::SparseMatrix< TScalar > _TSparseMatrix;
- typedef Eigen::SimplicialLDLT< _TSparseMatrix > _TSparseSolver;
- typedef Eigen::Triplet< TScalar > _TTriplet;
- typedef std::vector< _TTriplet > _TTriplets;
- typedef std::map< unsigned long, unsigned long > _TMap;
- typedef _TMap::value_type _TMapValue;
-
- // Some input values
- const TImage* in = this->GetInput( );
- TRegion region = in->GetRequestedRegion( );
+ // Useful typedefs
+ typedef Eigen::Triplet< TScalar > _TTriplet;
+ typedef std::vector< _TTriplet > _TTriplets;
+ typedef Eigen::SparseMatrix< TScalar > _TMatrix;
+ typedef Eigen::SimplicialLDLT< _TMatrix > _TSolver;
+
+ // Allocate outputs
+ this->AllocateOutputs( );
+
+ // Build boundary triplets and count labels
+ _TTriplets St, Bt;
+ std::map< TLabel, unsigned long > labels;
+ this->_Boundary( St, labels );
+ struct _TTripletsOrd
+ {
+ bool operator()( const _TTriplet& a, const _TTriplet& b )
+ {
+ return( a.row( ) < b.row( ) );
+ }
+ };
+ std::sort( St.begin( ), St.end( ), _TTripletsOrd( ) );
+ for( unsigned long i = 0; i < St.size( ); ++i )
+ Bt.push_back( _TTriplet( i, labels[ St[ i ].col( ) ], St[ i ].value( ) ) );
+
+ // Laplacian triplets
+ _TTriplets At, Rt;
+ this->_Laplacian( At, Rt, St );
+
+ // Matrices
+ TRegion region = this->GetInput( )->GetRequestedRegion( );
+ unsigned long nSeeds = St.size( );
+ unsigned long nLabels = labels.size( );
unsigned long N = region.GetNumberOfPixels( );
- // Prepare seeds for linear algebra
- _TMap seeds;
- for( unsigned int i = 0; i < this->m_Seeds.size( ); ++i )
- seeds.insert(
- _TMapValue(
- Self::_1D( this->m_Seeds[ i ], region ),
- this->m_Labels[ i ]
- )
- );
+ std::vector< TLabel > invLabels( nLabels );
+ for( typename std::map< TLabel, unsigned long >::value_type s: labels )
+ invLabels[ s.second ] = s.first;
- // Use input labels
- /* TODO
- const TLabels* in_labels = this->GetInputLabels( );
- if( in_labels != NULL )
- {
- itk::ImageRegionConstIteratorWithIndex< TLabels > lIt( in_labels, region );
- for( lIt.GoToBegin( ); !lIt.IsAtEnd( ); ++lIt )
- if( lIt.Get( ) != 0 )
- seeds.insert(
- _TMapValue( Self::_1D( lIt.GetIndex( ), region ), lIt.Get( ) )
- );
-
- } // fi
- */
-
- // Prepare label tables
- _TMap labels, inv_labels, seeds_indexes;
- _TMap::const_iterator hIt = seeds.begin( );
- for( ; hIt != seeds.end( ); ++hIt )
- {
- seeds_indexes[ hIt->first ] = seeds_indexes.size( ) - 1;
- if( labels.find( hIt->second ) == labels.end( ) )
- {
- labels[ hIt->second ] = labels.size( ) - 1;
- inv_labels[ labels[ hIt->second ] ] = hIt->second;
+ _TMatrix B( nSeeds, nLabels );
+ B.setFromTriplets( Bt.begin( ), Bt.end( ) );
+ B.makeCompressed( );
- } // fi
+ _TMatrix R( N - nSeeds, nSeeds );
+ R.setFromTriplets( Rt.begin( ), Rt.end( ) );
+ R.makeCompressed( );
- } // rof
+ _TMatrix A( N - nSeeds, N - nSeeds );
+ A.setFromTriplets( At.begin( ), At.end( ) );
+ A.makeCompressed( );
- // Prepare matrix/image index conversion
- _TMap indexes;
- unsigned long o = 0;
- for( unsigned long n = 0; n < N; ++n )
+ // Solve dirichlet problem
+ _TSolver solver;
+ solver.compute( A );
+ if( solver.info( ) != Eigen::Success )
{
- if( seeds.find( n ) == seeds.end( ) )
- indexes.insert( _TMap::value_type( n, n - o ) );
- else
- o++;
+ std::cerr << "Error computing." << std::endl;
+ } // fi
+ _TMatrix x = solver.solve( R * B );
+ if( solver.info( ) != Eigen::Success )
+ {
+ std::cerr << "Error solving." << std::endl;
+ } // fi
- } // rof
- unsigned long nLabels = labels.size( );
+ // Fill outputs
+ this->_Output( x, St, invLabels );
+}
- // Boundary matrix
- _TTriplets Bt;
- _TMap::const_iterator mIt = seeds.begin( );
- for( unsigned long i = 0; mIt != seeds.end( ); ++mIt, ++i )
- Bt.push_back( _TTriplet( i, labels[ mIt->second ], TScalar( 1 ) ) );
- _TSparseMatrix B( seeds.size( ), nLabels );
- B.setFromTriplets( Bt.begin( ), Bt.end( ) );
- B.makeCompressed( );
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalar >
+_TScalar fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
+_W( const TIndex& i, const TIndex& j )
+{
+ const TImage* in = this->GetInput( );
+ TScalar a = TScalar( in->GetPixel( i ) );
+ TScalar b = TScalar( in->GetPixel( j ) );
+ TScalar v = std::exp( -this->m_Beta * std::fabs( a - b ) );
+ if( v < this->m_Epsilon )
+ return( this->m_Epsilon );
+ else
+ return( v );
+}
- // Laplacian matrix
- _TTriplets Lt;
- itk::ImageRegionConstIteratorWithIndex< TImage > it( in, region );
- TScalar maxV = TScalar( -1 );
- for( it.GoToBegin( ); !it.IsAtEnd( ); ++it )
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalar >
+_TScalar fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
+_L( const TIndex& i, const TIndex& j )
+{
+ if( i == j )
{
- TIndex idx = it.GetIndex( );
- TScalar vidx = TScalar( it.Get( ) );
- unsigned long iidx = Self::_1D( idx, region );
-
- // Neighbors
+ TRegion r = this->GetInput( )->GetRequestedRegion( );
TScalar s = TScalar( 0 );
for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
{
- TIndex jdx;
- for( int l = -1; l <= 1; l += 2 )
+ for( int n = -1; n <= 1; n += 2 )
{
- jdx = idx;
- jdx[ d ] += l;
- if( region.IsInside( jdx ) )
- {
- TScalar vjdx = TScalar( in->GetPixel( jdx ) );
- unsigned long ijdx = Self::_1D( jdx, region );
- TScalar v = std::fabs( vidx - vjdx );
- Lt.push_back( _TTriplet( iidx, ijdx, v ) );
- if( maxV < v )
- maxV = v;
-
- } // fi
+ TIndex k = i;
+ k[ d ] += n;
+ if( r.IsInside( k ) )
+ s += this->_W( i, k );
} // rof
} // rof
+ return( s );
+ }
+ else
+ return( this->_W( i, j ) * TScalar( -1 ) );
+}
- } // rof
- std::vector< TScalar > diag( N, TScalar( 0 ) );
- typename _TTriplets::iterator tIt;
- TScalar betaV = -this->m_Beta;
- if( this->m_NormalizeWeights )
- betaV /= maxV;
- for( tIt = Lt.begin( ); tIt != Lt.end( ); ++tIt )
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalar >
+template< class _TTriplets >
+void fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
+_Boundary( _TTriplets& B, std::map< TLabel, unsigned long >& labels )
+{
+ B.clear( );
+
+ // Set up the multithreaded processing
+ _TBoundaryThreadStruct thrStr;
+ thrStr.Filter = this;
+ thrStr.Triplets = reinterpret_cast< void* >( &B );
+ thrStr.Labels = &labels;
+
+ // Configure threader
+ const TLabels* in_labels = this->GetInputLabels( );
+ const itk::ImageRegionSplitterBase* split = this->GetImageRegionSplitter( );
+ const unsigned int nThreads =
+ split->GetNumberOfSplits(
+ in_labels->GetRequestedRegion( ), this->GetNumberOfThreads( )
+ );
+
+ itk::MultiThreader::Pointer threads = itk::MultiThreader::New( );
+ threads->SetNumberOfThreads( nThreads );
+ threads->SetSingleMethod( this->_BoundaryCbk< _TTriplets >, &thrStr );
+
+ // Execute threader
+ threads->SingleMethodExecute( );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalar >
+template< class _TTriplets >
+ITK_THREAD_RETURN_TYPE
+fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
+_BoundaryCbk( void* arg )
+{
+ _TBoundaryThreadStruct* thrStr;
+ itk::ThreadIdType total, thrId, thrCount;
+ itk::MultiThreader::ThreadInfoStruct* thrInfo =
+ reinterpret_cast< itk::MultiThreader::ThreadInfoStruct* >( arg );
+ thrId = thrInfo->ThreadID;
+ thrCount = thrInfo->NumberOfThreads;
+ thrStr = reinterpret_cast< _TBoundaryThreadStruct* >( thrInfo->UserData );
+
+ TRegion region;
+ total = thrStr->Filter->SplitRequestedRegion( thrId, thrCount, region );
+ if( thrId < total )
+ thrStr->Filter->_ThreadedBoundary(
+ region, thrId,
+ reinterpret_cast< _TTriplets* >( thrStr->Triplets ),
+ thrStr->Labels
+ );
+ return( ITK_THREAD_RETURN_VALUE );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalar >
+template< class _TTriplets >
+void fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
+_ThreadedBoundary(
+ const TRegion& region, const itk::ThreadIdType& id,
+ _TTriplets* B,
+ std::map< TLabel, unsigned long >* labels
+ )
+{
+ typedef itk::ImageRegionConstIteratorWithIndex< TLabels > _TIt;
+ typedef typename std::map< TLabel, unsigned long >::value_type _TMapValue;
+ typedef typename std::map< unsigned long, TLabel >::value_type _TInvValue;
+ typedef typename _TTriplets::value_type _TTriplet;
+
+ const TLabels* in_labels = this->GetInputLabels( );
+ TRegion reqRegion = in_labels->GetRequestedRegion( );
+ _TIt it( in_labels, region );
+ for( it.GoToBegin( ); !it.IsAtEnd( ); ++it )
{
- TScalar v = std::exp( betaV * tIt->value( ) );
- if( v < this->m_Epsilon )
- v = this->m_Epsilon;
- *tIt = _TTriplet( tIt->row( ), tIt->col( ), -v );
- diag[ tIt->col( ) ] += v;
+ if( it.Get( ) != 0 )
+ {
+ unsigned long i = Self::_1D( it.GetIndex( ), reqRegion );
+ this->m_Mutex.Lock( );
+ B->push_back( _TTriplet( i, it.Get( ), TScalar( 1 ) ) );
+ if( labels->find( it.Get( ) ) == labels->end( ) )
+ labels->insert( _TMapValue( it.Get( ), labels->size( ) ) );
+ this->m_Mutex.Unlock( );
+
+ } // fi
} // rof
- for( unsigned long i = 0; i < diag.size( ); ++i )
- Lt.push_back( _TTriplet( i, i, diag[ i ] ) );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalar >
+template< class _TTriplets >
+void fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
+_Laplacian( _TTriplets& A, _TTriplets& R, const _TTriplets& B )
+{
+ A.clear( );
+ R.clear( );
+
+ // Set up the multithreaded processing
+ _TLaplacianThreadStruct thrStr;
+ thrStr.Filter = this;
+ thrStr.A = reinterpret_cast< void* >( &A );
+ thrStr.R = reinterpret_cast< void* >( &R );
+ thrStr.B = reinterpret_cast< const void* >( &B );
+
+ // Configure threader
+ const TImage* in = this->GetInputLabels( );
+ const itk::ImageRegionSplitterBase* split = this->GetImageRegionSplitter( );
+ const unsigned int nThreads =
+ split->GetNumberOfSplits(
+ in->GetRequestedRegion( ), this->GetNumberOfThreads( )
+ );
- // Compute R and A
- _TTriplets Rt, At;
- for( tIt = Lt.begin( ); tIt != Lt.end( ); ++tIt )
+ itk::MultiThreader::Pointer threads = itk::MultiThreader::New( );
+ threads->SetNumberOfThreads( nThreads );
+ threads->SetSingleMethod( this->_LaplacianCbk< _TTriplets >, &thrStr );
+
+ // Execute threader
+ threads->SingleMethodExecute( );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalar >
+template< class _TTriplets >
+ITK_THREAD_RETURN_TYPE
+fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
+_LaplacianCbk( void* arg )
+{
+ _TLaplacianThreadStruct* thrStr;
+ itk::ThreadIdType total, thrId, thrCount;
+ itk::MultiThreader::ThreadInfoStruct* thrInfo =
+ reinterpret_cast< itk::MultiThreader::ThreadInfoStruct* >( arg );
+ thrId = thrInfo->ThreadID;
+ thrCount = thrInfo->NumberOfThreads;
+ thrStr = reinterpret_cast< _TLaplacianThreadStruct* >( thrInfo->UserData );
+
+ TRegion region;
+ total = thrStr->Filter->SplitRequestedRegion( thrId, thrCount, region );
+ if( thrId < total )
+ thrStr->Filter->_ThreadedLaplacian(
+ region, thrId,
+ reinterpret_cast< _TTriplets* >( thrStr->A ),
+ reinterpret_cast< _TTriplets* >( thrStr->R ),
+ reinterpret_cast< const _TTriplets* >( thrStr->B )
+ );
+ return( ITK_THREAD_RETURN_VALUE );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalar >
+template< class _TTriplets >
+void fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
+_ThreadedLaplacian(
+ const TRegion& region, const itk::ThreadIdType& id,
+ _TTriplets* A, _TTriplets* R, const _TTriplets* B
+ )
+{
+ typedef itk::ImageRegionConstIteratorWithIndex< TImage > _TIt;
+ typedef typename _TTriplets::value_type _TTriplet;
+
+ const TImage* in = this->GetInput( );
+ const TLabels* in_labels = this->GetInputLabels( );
+ TRegion rqRegion = in->GetRequestedRegion( );
+ _TIt it( in, region );
+ for( it.GoToBegin( ); !it.IsAtEnd( ); ++it )
{
- _TMap::const_iterator cIt, rIt;
- cIt = seeds.find( tIt->col( ) );
- rIt = seeds.find( tIt->row( ) );
- if( cIt != seeds.end( ) )
- {
- if( rIt == seeds.end( ) )
- {
- _TMap::const_iterator iIt = indexes.find( tIt->row( ) );
- _TMap::const_iterator jIt = seeds_indexes.find( cIt->first );
- Rt.push_back( _TTriplet( iIt->second, jIt->second, -tIt->value( ) ) );
+ TIndex idx = it.GetIndex( );
+ bool iSeed = ( in_labels->GetPixel( idx ) != 0 );
+ unsigned long i = Self::_1D( idx, rqRegion );
+ unsigned long si;
- } // fi
+ // A's diagonal values
+ if( !iSeed )
+ {
+ si = Self::_NearSeedIndex( i, *B );
+ this->m_Mutex.Lock( );
+ A->push_back( _TTriplet( si, si, this->_L( idx, idx ) ) );
+ this->m_Mutex.Unlock( );
}
else
+ si = Self::_SeedIndex( i, *B );
+
+ // Neighbors (final matrix is symmetric)
+ for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
{
- if( rIt == seeds.end( ) )
+ for( int s = -1; s <= 1; s += 2 )
{
- _TMap::const_iterator iIt = indexes.find( tIt->row( ) );
- _TMap::const_iterator jIt = indexes.find( tIt->col( ) );
- At.push_back( _TTriplet( iIt->second, jIt->second, tIt->value( ) ) );
+ TIndex jdx = idx;
+ jdx[ d ] += s;
+ if( rqRegion.IsInside( jdx ) )
+ {
+ TScalar L = this->_L( idx, jdx );
+ unsigned long j = Self::_1D( jdx, rqRegion );
+ bool jSeed = ( in_labels->GetPixel( jdx ) != 0 );
+ if( !jSeed )
+ {
+ unsigned long sj = Self::_NearSeedIndex( j, *B );
+ if( !iSeed )
+ {
+ this->m_Mutex.Lock( );
+ A->push_back( _TTriplet( si, sj, L ) );
+ this->m_Mutex.Unlock( );
+ }
+ else
+ {
+ this->m_Mutex.Lock( );
+ R->push_back( _TTriplet( sj, si, -L ) );
+ this->m_Mutex.Unlock( );
+
+ } // fi
+
+ } // fi
+
+ } // fi
- } // fi
+ } // rof
- } // fi
+ } // rof
} // rof
+}
- _TSparseMatrix R( N - seeds.size( ), seeds.size( ) );
- R.setFromTriplets( Rt.begin( ), Rt.end( ) );
- R.makeCompressed( );
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalar >
+template< class _TMatrix, class _TTriplets >
+void fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
+_Output(
+ const _TMatrix& X, const _TTriplets& S, const std::vector< TLabel >& invLabels
+ )
+{
+ // Set up the multithreaded processing
+ _TOutputThreadStruct thrStr;
+ thrStr.Filter = this;
+ thrStr.X = reinterpret_cast< const void* >( &X );
+ thrStr.S = reinterpret_cast< const void* >( &S );
+ thrStr.InvLabels = &invLabels;
+
+ // Configure threader
+ const TLabels* out = this->GetOutput( );
+ const itk::ImageRegionSplitterBase* split = this->GetImageRegionSplitter( );
+ const unsigned int nThreads =
+ split->GetNumberOfSplits(
+ out->GetRequestedRegion( ), this->GetNumberOfThreads( )
+ );
- _TSparseMatrix A( N - seeds.size( ), N - seeds.size( ) );
- A.setFromTriplets( At.begin( ), At.end( ) );
- A.makeCompressed( );
+ itk::MultiThreader::Pointer threads = itk::MultiThreader::New( );
+ threads->SetNumberOfThreads( nThreads );
+ threads->SetSingleMethod(
+ this->_OutputCbk< _TMatrix, _TTriplets >, &thrStr
+ );
- // Solve dirichlet problem
- _TSparseSolver solver;
- solver.compute( A );
- if( solver.info( ) != Eigen::Success )
- {
- std::cerr << "Error computing." << std::endl;
- } // fi
- _TSparseMatrix x = solver.solve( R * B );
- if( solver.info( ) != Eigen::Success )
- {
- std::cerr << "Error solving." << std::endl;
- } // fi
+ // Execute threader
+ threads->SingleMethodExecute( );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalar >
+template< class _TMatrix, class _TTriplets >
+ITK_THREAD_RETURN_TYPE
+fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
+_OutputCbk( void* arg )
+{
+ _TOutputThreadStruct* thrStr;
+ itk::ThreadIdType total, thrId, thrCount;
+ itk::MultiThreader::ThreadInfoStruct* thrInfo =
+ reinterpret_cast< itk::MultiThreader::ThreadInfoStruct* >( arg );
+ thrId = thrInfo->ThreadID;
+ thrCount = thrInfo->NumberOfThreads;
+ thrStr = reinterpret_cast< _TOutputThreadStruct* >( thrInfo->UserData );
+
+ TRegion region;
+ total = thrStr->Filter->SplitRequestedRegion( thrId, thrCount, region );
+ if( thrId < total )
+ thrStr->Filter->_ThreadedOutput(
+ region, thrId,
+ reinterpret_cast< const _TMatrix* >( thrStr->X ),
+ reinterpret_cast< const _TTriplets* >( thrStr->S ),
+ thrStr->InvLabels
+ );
+ return( ITK_THREAD_RETURN_VALUE );
+}
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalar >
+template< class _TMatrix, class _TTriplets >
+void fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
+_ThreadedOutput(
+ const TRegion& region, const itk::ThreadIdType& id,
+ const _TMatrix* X, const _TTriplets* S,
+ const std::vector< TLabel >* invLabels
+ )
+{
// Fill outputs
+ const TLabels* in_labels = this->GetInputLabels( );
TLabels* out_labels = this->GetOutput( );
- out_labels->SetLargestPossibleRegion( in->GetLargestPossibleRegion( ) );
- out_labels->SetRequestedRegion( in->GetRequestedRegion( ) );
- out_labels->SetBufferedRegion( in->GetBufferedRegion( ) );
- out_labels->SetSpacing( in->GetSpacing( ) );
- out_labels->SetOrigin( in->GetOrigin( ) );
- out_labels->SetDirection( in->GetDirection( ) );
- out_labels->Allocate( );
-
TScalarImage* out_probs = this->GetOutputProbabilities( );
- out_probs->SetLargestPossibleRegion( in->GetLargestPossibleRegion( ) );
- out_probs->SetRequestedRegion( in->GetRequestedRegion( ) );
- out_probs->SetBufferedRegion( in->GetBufferedRegion( ) );
- out_probs->SetSpacing( in->GetSpacing( ) );
- out_probs->SetOrigin( in->GetOrigin( ) );
- out_probs->SetDirection( in->GetDirection( ) );
- out_probs->Allocate( );
-
- mIt = seeds.begin( );
- o = 0;
- unsigned long j = 0;
- for( unsigned long j = 0; j < N; ++j )
+ itk::ImageRegionConstIteratorWithIndex< TLabels > iIt( in_labels, region );
+ itk::ImageRegionIteratorWithIndex< TLabels > oIt( out_labels, region );
+ itk::ImageRegionIteratorWithIndex< TScalarImage > pIt( out_probs, region );
+ iIt.GoToBegin( );
+ oIt.GoToBegin( );
+ pIt.GoToBegin( );
+ for( ; !iIt.IsAtEnd( ); ++iIt, ++oIt, ++pIt )
{
- TIndex idx;
- TLabel lbl;
- TScalar p;
- if( mIt->first != j )
+ if( iIt.Get( ) == 0 )
{
- idx = Self::_ND( j, region );
- TScalar maxP = x.coeffRef( j - o, 0 );
+ unsigned long i = Self::_1D( iIt.GetIndex( ), region );
+ unsigned long j = Self::_NearSeedIndex( i, *S );
+ TScalar maxP = X->coeff( j, 0 );
unsigned long maxL = 0;
- for( unsigned long l = 1; l < nLabels; ++l )
+ for( unsigned int s = 1; s < invLabels->size( ); ++s )
{
- TScalar vp = x.coeffRef( j - o, l );
- if( maxP < vp )
+ TScalar p = X->coeff( j, s );
+ if( maxP <= p )
{
- maxP = vp;
- maxL = l;
+ maxP = p;
+ maxL = s;
} // fi
} // rof
- lbl = inv_labels[ maxL ];
- p = maxP;
+ oIt.Set( ( *invLabels )[ maxL ] );
+ pIt.Set( maxP );
}
else
{
- idx = Self::_ND( mIt->first, region );
- lbl = mIt->second;
- p = TScalar( 1 );
- mIt++;
- o++;
+ oIt.Set( iIt.Get( ) );
+ pIt.Set( TScalar( 1 ) );
} // fi
- out_labels->SetPixel( idx, lbl );
- out_probs->SetPixel( idx, p );
} // rof
}
// -------------------------------------------------------------------------
template< class _TImage, class _TLabels, class _TScalar >
-typename fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
-TIndex fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
-_ND( const unsigned long& i, const TRegion& region )
+template< class _TTriplets >
+unsigned long
+fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
+_SeedIndex( const unsigned long& i, const _TTriplets& t )
{
- typename TRegion::SizeType size = region.GetSize( );
-
- unsigned long j = i;
- TIndex idx;
- if( TIndex::Dimension == 3 )
+ unsigned long s = 0;
+ unsigned long f = t.size( );
+ unsigned long e = f - 1;
+ while( e > s && f == t.size( ) )
{
- unsigned long z = size[ 0 ] * size[ 1 ];
- idx[ 2 ] = j / z;
- j -= idx[ 2 ] * z;
+ if( e > s + 1 )
+ {
+ unsigned long h = ( e + s ) >> 1;
+ if ( i < t[ h ].row( ) ) e = h;
+ else if( t[ h ].row( ) < i ) s = h;
+ else f = h;
+ }
+ else
+ f = ( t[ s ].row( ) == i )? s: e;
- } // fi
- idx[ 1 ] = j / size[ 0 ];
- idx[ 0 ] = j % size[ 0 ];
- return( idx );
+ } // elihw
+ return( f );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage, class _TLabels, class _TScalar >
+template< class _TTriplets >
+unsigned long
+fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar >::
+_NearSeedIndex( const unsigned long& i, const _TTriplets& t )
+{
+ long s = 0;
+ long e = t.size( ) - 1;
+ while( e > s + 1 )
+ {
+ long h = ( e + s ) >> 1;
+ if ( i < t[ h ].row( ) ) e = h;
+ else if( t[ h ].row( ) < i ) s = h;
+
+ } // elihw
+ long d;
+ if( i < t[ s ].row( ) )
+ d = -1;
+ else if( t[ s ].row( ) < i && i < t[ e ].row( ) )
+ d = s + 1;
+ else
+ d = e + 1;
+ if( d < 0 ) d = 0;
+ return( i - d );
}
#endif // __fpa__Common__OriginalRandomWalker__hxx__
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / commitdiff