#define __fpa__Common__SliceBySliceRandomWalker__hxx__
#include <vector>
+#include <itkBinaryThresholdImageFilter.h>
#include <itkExtractImageFilter.h>
#include <itkImageRegionConstIterator.h>
#include <itkImageRegionIterator.h>
#include <itkJoinSeriesImageFilter.h>
+#include <itkMaximumImageFilter.h>
#include <itkMinimumMaximumImageCalculator.h>
#include <fpa/Common/RandomWalker.h>
#include <fpa/Functors/Dijkstra/Image/Gaussian.h>
template< class _TImage, class _TLabels, class _TScalarImage >
fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
SliceBySliceRandomWalker( )
- : Superclass( )
+ : Superclass( ),
+ m_Epsilon( TScalar( 1e-5 ) ),
+ m_Beta( TScalar( 1 ) ),
+ m_VesselnessThreshold( TScalar( 5 ) )
{
ivqITKInputConfigureMacro( InputLabels, TLabels );
ivqITKInputConfigureMacro( InputVesselness, TScalarImage );
const TScalarImage* vesselness = this->GetInputVesselness( );
this->AllocateOutputs( );
+ // Get vesselness range
+ typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
+ typename _TMinMax::Pointer vMinMax = _TMinMax::New( );
+ vMinMax->SetImage( vesselness );
+ vMinMax->Compute( );
+ TScalar vMax = vMinMax->GetMaximum( );
+
// Some values
typename TImage::RegionType region = input->GetRequestedRegion( );
- typename TImage::SizeType regionIndex = region.GetIndex( );
+ typename TImage::IndexType regionIndex = region.GetIndex( );
typename TImage::SizeType sliceSize = region.GetSize( );
int numSlices = sliceSize[ 2 ];
sliceSize[ 2 ] = 0;
- this->m_VesselnessValue =
- this->m_VesselnessThreshold * this->m_VesselnessMax / TScalar( 100 );
+ this->m_VesselnessValue = TScalar( 0.01 );
+ this->m_VesselnessValue *= this->m_VesselnessThreshold * vMax;
// Composite image
- typename TScalarImage::Pointer composite = TScalarImage::New( );
- composite->CopyInformation( vesselness );
- composite->SetBufferedRegion( vesselness->GetBufferedRegion( ) );
- composite->Allocate( );
- this->_Composite( composite, labels, vesselness );
+ typename TScalarImage::Pointer composite;
+ this->_Composite( composite, labels, vesselness, vMax );
// Extract slices
std::vector< typename _TSliceImage::Pointer > data3D( numSlices );
typename _TSliceScalarImage::Pointer composite_slice;
this->_Slice( input_slice, input, sliceRegion );
this->_Slice( labels_slice, labels, sliceRegion );
- this->_Slice( composite_slice, composite, sliceRegion );
+ this->_Slice( composite_slice, composite.GetPointer( ), sliceRegion );
// Update vectors with each image
data3D[ i ] = input_slice;
} // rof
// Random walker slice-by-slice
- this->_GoDown( binaryTree, data3D, fusion, numSlices );
- this->_GoUp( binaryTree, data3D, fusion, numSlices );
+ this->_RandomWalker( binaryTree, data3D, fusion, true );
+ this->_RandomWalker( binaryTree, data3D, fusion, false );
// Join results
- typedef itk::JoinSeriesImageFilter< _TSliceImage, TImage > _TJoin;
+ typedef itk::JoinSeriesImageFilter< _TSliceLabels, TLabels > _TJoin;
typename _TJoin::Pointer join = _TJoin::New( );
for( int i = 0; i < numSlices; ++i )
join->SetInput( i, binaryTree[ i ] );
+ join->SetOrigin( input->GetOrigin( )[ 2 ] );
+ join->SetSpacing( input->GetSpacing( )[ 2 ] );
join->Update( );
this->GetOutput( )->Graft( join->GetOutput( ) );
}
_Composite(
typename TScalarImage::Pointer& composite,
const TLabels* labels,
- const TScalarImage* vesselness
+ const TScalarImage* vesselness,
+ const TScalar& maxVess
)
{
- // Min-Max vesselness values
- typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
- typename _TMinMax::Pointer minMax = _TMinMax::New( );
- minMax->SetImage( vesselness );
- minMax->Compute( );
- TScalar maxVess = minMax->GetMaximum( );
+ // Some values
typename TScalarImage::RegionType region = labels->GetRequestedRegion( );
// Composite image
- typename TScalarImage::Pointer composite = TScalarImage::New( );
+ composite = TScalarImage::New( );
composite->CopyInformation( vesselness );
composite->SetBufferedRegion( vesselness->GetBufferedRegion( ) );
composite->Allocate( );
typename _TInput::RegionType region
)
{
- typedef typename _TSlicePTr::ObjectType _TSlice;
+ typedef typename _TSlicePtr::ObjectType _TSlice;
typedef itk::ExtractImageFilter< _TInput, _TSlice > _TExtract;
typename _TExtract::Pointer extract = _TExtract::New( );
template< class _TBinaryTree, class _TData3D, class _TFusion >
void
fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
-_GoDown(
+_RandomWalker(
_TBinaryTree& binaryTree,
const _TData3D& data3D,
const _TFusion& fusion,
- int numSlices
+ bool down
)
{
typedef typename _TBinaryTree::value_type _TSliceBinaryPtr;
typedef typename _TSliceFusionPtr::ObjectType _TSliceFusion;
typedef typename _TSliceFusion::PixelType _TSliceScalar;
- int z = -1;
- while( z < numSlices - 2 )
+ int numSlices = binaryTree.size( );
+ int z = ( down )? -1: numSlices;
+ int o = ( down )? 1: -1;
+ while( ( down && z < numSlices - 2 ) || ( !down && z > 1 ) )
{
- z++;
+ z += o;
_TSliceBinaryPtr tmp = binaryTree[ z ];
- _TSliceBinaryPtr next = binaryTree[ z + 1 ];
- _TSliceData3DPtr data = data3D[ z + 1 ];
- _TSliceFusionPtr vess = fusion[ z + 1 ];
+ _TSliceBinaryPtr next = binaryTree[ z + o ];
+ _TSliceData3DPtr data = data3D[ z + o ];
+ _TSliceFusionPtr vess = fusion[ z + o ];
typename _TSliceBinary::RegionType region = tmp->GetRequestedRegion( );
// Create seeds and background
_TSliceBinaryPtr seeds = _TSliceBinary::New( );
seeds->CopyInformation( tmp );
- seeds->SetRequestedRegion( region );
+ seeds->SetBufferedRegion( region );
seeds->Allocate( );
seeds->FillBuffer( 0 );
// Random walker function
typedef fpa::Functors::Dijkstra::Image::Gaussian< _TSliceData3D, _TSliceScalar > _TFunction;
- typename TFunction::Pointer edge = TFunction::New( );
+ typename _TFunction::Pointer edge = _TFunction::New( );
edge->SetBeta( this->m_Beta );
- edge->SetEpsilon( this->m_Eps );
+ edge->SetEpsilon( this->m_Epsilon );
edge->TreatAsWeightOff( );
// Real random walker
rw->SetInputLabels( seeds );
rw->SetEdgeFunction( edge );
- // Keep maximum
- /* TODO
- typedef itk::MaximumImageFilter<SliceType> MaxFilterType;
- MaxFilterType::Pointer maxFilter = MaxFilterType::New();
- maxFilter->SetInput(0, img2);
- maxFilter->SetInput(1, next);
- maxFilter->Update();
-
- SliceType::Pointer sth = maxFilter->GetOutput();
- ImageAlgorithm::Copy(sth.GetPointer(), binaryTree[z + 1].GetPointer(), sth->GetRequestedRegion(), binaryTree[z + 1]->GetRequestedRegion());
- delete[] probs;
- */
- }
-}
+ // Extract label
+ typedef itk::BinaryThresholdImageFilter< _TSliceBinary, _TSliceBinary > _TExtract;
+ typename _TExtract::Pointer extract = _TExtract::New( );
+ extract->SetInput( rw->GetOutput( ) );
+ extract->SetInsideValue( 1 );
+ extract->SetOutsideValue( 0 );
+ extract->SetLowerThreshold( 1 );
+ extract->SetUpperThreshold( 1 );
-// -------------------------------------------------------------------------
-template< class _TImage, class _TLabels, class _TScalarImage >
-template< class _TBinaryTree, class _TData3D, class _TFusion >
-void
-_GoUp(
- _TBinaryTree& binaryTree,
- const _TData3D& data3D,
- const _TFusion& fusion,
- int numSlices
- )
-{
+ // Keep maximum
+ typedef itk::MaximumImageFilter< _TSliceBinary > _TMax;
+ typename _TMax::Pointer max = _TMax::New();
+ max->SetInput( 0, extract->GetOutput( ) );
+ max->SetInput( 1, next );
+ max->Update( );
+ binaryTree[ z + o ] = max->GetOutput( );
+ binaryTree[ z + o ]->DisconnectPipeline( );
+
+ } // elihw
}
#endif // __fpa__Common__SliceBySliceRandomWalker__hxx__