From: Leonardo Flórez-Valencia Date: Wed, 20 Sep 2017 15:32:14 +0000 (+0200) Subject: ... X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=34f4ff5d31c70f1127d43865c61e9f57a7071190;p=FrontAlgorithms.git ... --- diff --git a/lib/fpa/Common/OriginalRandomWalker.h b/lib/fpa/Common/OriginalRandomWalker.h index 915cb19..d4af732 100644 --- a/lib/fpa/Common/OriginalRandomWalker.h +++ b/lib/fpa/Common/OriginalRandomWalker.h @@ -6,9 +6,11 @@ #define __fpa__Common__OriginalRandomWalker__h__ #include +#include #include #include #include +#include namespace fpa { @@ -36,6 +38,30 @@ namespace fpa typedef typename TImage::RegionType TRegion; typedef typename TLabels::PixelType TLabel; + protected: + struct _TBoundaryThreadStruct + { + Pointer Filter; + void* Triplets; + std::map< TLabel, unsigned long >* Labels; + }; + + struct _TLaplacianThreadStruct + { + Pointer Filter; + void* A; + void* R; + const void* B; + }; + + struct _TOutputThreadStruct + { + Pointer Filter; + const void* X; + const void* S; + const std::vector< TLabel >* InvLabels; + }; + public: itkNewMacro( Self ); itkTypeMacro( @@ -64,8 +90,70 @@ namespace fpa virtual void GenerateData( ) override; + virtual TScalar _W( const TIndex& i, const TIndex& j ); + virtual TScalar _L( const TIndex& i, const TIndex& j ); + + // Boundary construction + template< class _TTriplets > + void _Boundary( + _TTriplets& B, + std::map< TLabel, unsigned long >& labels + ); + + template< class _TTriplets > + static ITK_THREAD_RETURN_TYPE _BoundaryCbk( void* arg ); + + template< class _TTriplets > + void _ThreadedBoundary( + const TRegion& region, const itk::ThreadIdType& id, + _TTriplets* B, + std::map< TLabel, unsigned long >* labels + ); + + // Laplacian construction + template< class _TTriplets > + void _Laplacian( + _TTriplets& A, _TTriplets& R, const _TTriplets& B + ); + + template< class _TTriplets > + static ITK_THREAD_RETURN_TYPE _LaplacianCbk( void* arg ); + + template< class _TTriplets > + void _ThreadedLaplacian( + const TRegion& region, const itk::ThreadIdType& id, + _TTriplets* A, _TTriplets* R, const _TTriplets* B + ); + + // Output filling + template< class _TMatrix, class _TTriplets > + void _Output( + const _TMatrix& X, const _TTriplets& S, + const std::vector< TLabel >& invLabels + ); + + template< class _TMatrix, class _TTriplets > + static ITK_THREAD_RETURN_TYPE _OutputCbk( void* arg ); + + template< class _TMatrix, class _TTriplets > + void _ThreadedOutput( + const TRegion& region, const itk::ThreadIdType& id, + const _TMatrix* X, const _TTriplets* S, + const std::vector< TLabel >* invLabels + ); + + // Array-Image conversion static unsigned long _1D( const TIndex& idx, const TRegion& region ); - static TIndex _ND( const unsigned long& i, const TRegion& region ); + + template< class _TTriplets > + static unsigned long _SeedIndex( + const unsigned long& i, const _TTriplets& t + ); + + template< class _TTriplets > + static unsigned long _NearSeedIndex( + const unsigned long& i, const _TTriplets& t + ); private: // Purposely not implemented @@ -79,6 +167,8 @@ namespace fpa TScalar m_Beta; TScalar m_Epsilon; bool m_NormalizeWeights; + + itk::SimpleFastMutexLock m_Mutex; }; } // ecapseman diff --git a/lib/fpa/Common/OriginalRandomWalker.hxx b/lib/fpa/Common/OriginalRandomWalker.hxx index 7c29595..9bc08d4 100644 --- a/lib/fpa/Common/OriginalRandomWalker.hxx +++ b/lib/fpa/Common/OriginalRandomWalker.hxx @@ -6,8 +6,8 @@ #define __fpa__Common__OriginalRandomWalker__hxx__ #include -#include #include +#include #include // ------------------------------------------------------------------------- @@ -29,7 +29,7 @@ OriginalRandomWalker( ) m_Epsilon( TScalar( 1e-5 ) ), m_NormalizeWeights( true ) { - fpaFilterOptionalInputConfigureMacro( InputLabels, TLabels ); + fpaFilterInputConfigureMacro( InputLabels, TLabels ); fpaFilterOutputConfigureMacro( OutputProbabilities, TScalarImage ); } @@ -45,241 +45,446 @@ template< class _TImage, class _TLabels, class _TScalar > 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 } @@ -304,24 +509,55 @@ _1D( const TIndex& idx, const TRegion& region ) // ------------------------------------------------------------------------- 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__ diff --git a/lib/fpa/Config.h.in b/lib/fpa/Config.h.in index e850355..831fb72 100644 --- a/lib/fpa/Config.h.in +++ b/lib/fpa/Config.h.in @@ -53,13 +53,6 @@ this->m_##__n__##Idx + 1 \ ) -// ------------------------------------------------------------------------- -#define fpaFilterOptionalInputConfigureMacro( __n__, __t__ ) \ - this->m_##__n__##Idx = this->GetNumberOfIndexedInputs( ); \ - this->itk::ProcessObject::SetNumberOfIndexedInputs( \ - this->m_##__n__##Idx + 1 \ - ) - // ------------------------------------------------------------------------- #define fpaFilterOutputMacro( __n__, __t__ ) \ private: \ diff --git a/tests/image/RandomWalker/Original.cxx b/tests/image/RandomWalker/Original.cxx index 5a6c6df..46ba198 100644 --- a/tests/image/RandomWalker/Original.cxx +++ b/tests/image/RandomWalker/Original.cxx @@ -9,7 +9,7 @@ #include // ------------------------------------------------------------------------- -const unsigned int Dim = 3; +const unsigned int Dim = 2; typedef unsigned char TPixel; typedef unsigned char TLabel; @@ -19,45 +19,47 @@ typedef itk::Image< TLabel, Dim > TLabels; // ------------------------------------------------------------------------- int main( int argc, char* argv[] ) { - // Get arguments - if( argc < 4 ) + // Get input arguments + if( argc < 8 ) { std::cerr << "Usage: " << argv[ 0 ] - << " input_image output_labels output_probabilities" + << " input_image input_labels output_labels output_probabilities" + << " beta eps normalize" << std::endl; return( 1 ); } // fi - std::string input_image_filename = argv[ 1 ]; - std::string output_labels_filename = argv[ 2 ]; - std::string output_probabilities_filename = argv[ 3 ]; - double beta = 90; - double eps = 1e-5; - - TImage::IndexType s1, s2; - s1.Fill( 0 ); - s2.Fill( 0 ); - s1[ 0 ] = 131; - s1[ 1 ] = 150; - s2[ 0 ] = 200; - s2[ 1 ] = 200; + std::string input_image_filename, input_labels_filename; + std::string output_labels_filename, output_probabilities_filename; + double beta, eps; + bool normalize; + input_image_filename = argv[ 1 ]; + input_labels_filename = argv[ 2 ]; + output_labels_filename = argv[ 3 ]; + output_probabilities_filename = argv[ 4 ]; + beta = std::atof( argv[ 5 ] ); + eps = std::atof( argv[ 6 ] ); + normalize = ( argv[ 7 ][ 0 ] == '1' ); // Read image - TImage::Pointer input; typedef itk::ImageFileReader< TImage > TImageReader; TImageReader::Pointer input_image_reader = TImageReader::New( ); input_image_reader->SetFileName( input_image_filename ); + // Read labels + typedef itk::ImageFileReader< TLabels > TLabelsReader; + TLabelsReader::Pointer input_labels_reader = TLabelsReader::New( ); + input_labels_reader->SetFileName( input_labels_filename ); + // Random walker typedef fpa::Common::OriginalRandomWalker< TImage, TLabels > TFilter; TFilter::Pointer rw = TFilter::New( ); rw->SetInput( input_image_reader->GetOutput( ) ); + rw->SetInputLabels( input_labels_reader->GetOutput( ) ); rw->SetBeta( beta ); rw->SetEpsilon( eps ); - rw->AddSeed( s1, 100 ); - rw->AddSeed( s2, 150 ); - rw->NormalizeWeightsOn( ); + rw->SetNormalizeWeights( normalize ); // Save labels typedef itk::ImageFileWriter< TFilter::TLabels > TLabelsWriter;