1 // =========================================================================
2 // @author Leonardo Florez Valencia
3 // @email florez-l@javeriana.edu.co
4 // =========================================================================
5 #ifndef __fpa__Common__SliceBySliceRandomWalker__hxx__
6 #define __fpa__Common__SliceBySliceRandomWalker__hxx__
9 #include <itkBinaryThresholdImageFilter.h>
10 #include <itkExtractImageFilter.h>
11 #include <itkImageRegionConstIterator.h>
12 #include <itkImageRegionIterator.h>
13 #include <itkJoinSeriesImageFilter.h>
14 #include <itkMaximumImageFilter.h>
15 #include <itkMinimumMaximumImageCalculator.h>
16 #include <fpa/Common/RandomWalker.h>
17 #include <fpa/Functors/Dijkstra/Image/Gaussian.h>
19 // -------------------------------------------------------------------------
20 template< class _TImage, class _TLabels, class _TScalarImage >
21 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
22 SliceBySliceRandomWalker( )
24 m_Epsilon( TScalar( 1e-5 ) ),
25 m_Beta( TScalar( 1 ) ),
26 m_VesselnessThreshold( TScalar( 5 ) )
28 ivqITKInputConfigureMacro( InputLabels, TLabels );
29 ivqITKInputConfigureMacro( InputVesselness, TScalarImage );
32 // -------------------------------------------------------------------------
33 template< class _TImage, class _TLabels, class _TScalarImage >
34 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
35 ~SliceBySliceRandomWalker( )
39 // -------------------------------------------------------------------------
40 template< class _TImage, class _TLabels, class _TScalarImage >
42 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
45 typedef itk::Image< TPixel, TImage::ImageDimension - 1 > _TSliceImage;
46 typedef itk::Image< TLabel, TImage::ImageDimension - 1 > _TSliceLabels;
47 typedef itk::Image< TScalar, TImage::ImageDimension - 1 > _TSliceScalarImage;
50 const TImage* input = this->GetInput( );
51 const TLabels* labels = this->GetInputLabels( );
52 const TScalarImage* vesselness = this->GetInputVesselness( );
53 this->AllocateOutputs( );
55 // Get vesselness range
56 typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
57 typename _TMinMax::Pointer vMinMax = _TMinMax::New( );
58 vMinMax->SetImage( vesselness );
60 TScalar vMax = vMinMax->GetMaximum( );
63 typename TImage::RegionType region = input->GetRequestedRegion( );
64 typename TImage::IndexType regionIndex = region.GetIndex( );
65 typename TImage::SizeType sliceSize = region.GetSize( );
66 int numSlices = sliceSize[ 2 ];
68 this->m_VesselnessValue = TScalar( 0.01 );
69 this->m_VesselnessValue *= this->m_VesselnessThreshold * vMax;
72 typename TScalarImage::Pointer composite;
73 this->_Composite( composite, labels, vesselness, vMax );
76 std::vector< typename _TSliceImage::Pointer > data3D( numSlices );
77 std::vector< typename _TSliceLabels::Pointer > binaryTree( numSlices );
78 std::vector< typename _TSliceScalarImage::Pointer > fusion( numSlices );
79 for( int i = 0; i < numSlices; ++i )
81 // Prepare extraction region
82 typename TImage::IndexType sliceIndex = regionIndex;
84 typename TImage::RegionType sliceRegion( sliceIndex, sliceSize );
87 typename _TSliceImage::Pointer input_slice;
88 typename _TSliceLabels::Pointer labels_slice;
89 typename _TSliceScalarImage::Pointer composite_slice;
90 this->_Slice( input_slice, input, sliceRegion );
91 this->_Slice( labels_slice, labels, sliceRegion );
92 this->_Slice( composite_slice, composite.GetPointer( ), sliceRegion );
94 // Update vectors with each image
95 data3D[ i ] = input_slice;
96 binaryTree[ i ] = labels_slice;
97 fusion[ i ] = composite_slice;
101 // Random walker slice-by-slice
102 this->_RandomWalker( binaryTree, data3D, fusion, true );
103 this->_RandomWalker( binaryTree, data3D, fusion, false );
106 typedef itk::JoinSeriesImageFilter< _TSliceLabels, TLabels > _TJoin;
107 typename _TJoin::Pointer join = _TJoin::New( );
108 for( int i = 0; i < numSlices; ++i )
109 join->SetInput( i, binaryTree[ i ] );
110 join->SetOrigin( input->GetOrigin( )[ 2 ] );
111 join->SetSpacing( input->GetSpacing( )[ 2 ] );
113 this->GetOutput( )->Graft( join->GetOutput( ) );
116 // -------------------------------------------------------------------------
117 template< class _TImage, class _TLabels, class _TScalarImage >
119 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
121 typename TScalarImage::Pointer& composite,
122 const TLabels* labels,
123 const TScalarImage* vesselness,
124 const TScalar& maxVess
128 typename TScalarImage::RegionType region = labels->GetRequestedRegion( );
131 composite = TScalarImage::New( );
132 composite->CopyInformation( vesselness );
133 composite->SetBufferedRegion( vesselness->GetBufferedRegion( ) );
134 composite->Allocate( );
135 itk::ImageRegionConstIterator< TScalarImage > v( vesselness, region );
136 itk::ImageRegionConstIterator< TLabels > l( labels, region );
137 itk::ImageRegionIterator< TScalarImage > c( composite, region );
141 for( ; !v.IsAtEnd( ); ++v, ++l, ++c )
142 c.Set( ( l.Get( ) != 0 )? maxVess: v.Get( ) );
145 // -------------------------------------------------------------------------
146 template< class _TImage, class _TLabels, class _TScalarImage >
147 template< class _TSlicePtr, class _TInput >
149 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
152 const _TInput* input,
153 typename _TInput::RegionType region
156 typedef typename _TSlicePtr::ObjectType _TSlice;
157 typedef itk::ExtractImageFilter< _TInput, _TSlice > _TExtract;
159 typename _TExtract::Pointer extract = _TExtract::New( );
160 extract->SetInput( input );
161 extract->SetExtractionRegion( region );
162 extract->SetDirectionCollapseToIdentity( );
164 slice = extract->GetOutput( );
165 slice->DisconnectPipeline( );
168 // -------------------------------------------------------------------------
169 template< class _TImage, class _TLabels, class _TScalarImage >
170 template< class _TBinaryTree, class _TData3D, class _TFusion >
172 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
174 _TBinaryTree& binaryTree,
175 const _TData3D& data3D,
176 const _TFusion& fusion,
180 typedef typename _TBinaryTree::value_type _TSliceBinaryPtr;
181 typedef typename _TData3D::value_type _TSliceData3DPtr;
182 typedef typename _TFusion::value_type _TSliceFusionPtr;
183 typedef typename _TSliceBinaryPtr::ObjectType _TSliceBinary;
184 typedef typename _TSliceData3DPtr::ObjectType _TSliceData3D;
185 typedef typename _TSliceFusionPtr::ObjectType _TSliceFusion;
186 typedef typename _TSliceFusion::PixelType _TSliceScalar;
188 int numSlices = binaryTree.size( );
189 int z = ( down )? -1: numSlices;
190 int o = ( down )? 1: -1;
191 while( ( down && z < numSlices - 2 ) || ( !down && z > 1 ) )
194 _TSliceBinaryPtr tmp = binaryTree[ z ];
195 _TSliceBinaryPtr next = binaryTree[ z + o ];
196 _TSliceData3DPtr data = data3D[ z + o ];
197 _TSliceFusionPtr vess = fusion[ z + o ];
198 typename _TSliceBinary::RegionType region = tmp->GetRequestedRegion( );
200 // Create seeds and background
201 _TSliceBinaryPtr seeds = _TSliceBinary::New( );
202 seeds->CopyInformation( tmp );
203 seeds->SetBufferedRegion( region );
205 seeds->FillBuffer( 0 );
207 itk::ImageRegionConstIterator< _TSliceBinary > tIt( tmp, region );
208 itk::ImageRegionConstIterator< _TSliceFusion > vIt( vess, region );
209 itk::ImageRegionIterator< _TSliceBinary > sIt( seeds, region );
216 for( ; !tIt.IsAtEnd( ); ++tIt, ++vIt, ++sIt )
218 if( tIt.Get( ) != 0 )
221 if( vIt.Get( ) > this->m_VesselnessValue )
230 if( vIt.Get( ) < this->m_VesselnessValue )
242 // Random walker function
243 typedef fpa::Functors::Dijkstra::Image::Gaussian< _TSliceData3D, _TSliceScalar > _TFunction;
244 typename _TFunction::Pointer edge = _TFunction::New( );
245 edge->SetBeta( this->m_Beta );
246 edge->SetEpsilon( this->m_Epsilon );
247 edge->TreatAsWeightOff( );
249 // Real random walker
250 typedef fpa::Common::RandomWalker< _TSliceData3D, _TSliceBinary, _TSliceScalar > _TRandomWalker;
251 typename _TRandomWalker::Pointer rw = _TRandomWalker::New( );
252 rw->SetInput( data );
253 rw->SetInputLabels( seeds );
254 rw->SetEdgeFunction( edge );
257 typedef itk::BinaryThresholdImageFilter< _TSliceBinary, _TSliceBinary > _TExtract;
258 typename _TExtract::Pointer extract = _TExtract::New( );
259 extract->SetInput( rw->GetOutput( ) );
260 extract->SetInsideValue( 1 );
261 extract->SetOutsideValue( 0 );
262 extract->SetLowerThreshold( 1 );
263 extract->SetUpperThreshold( 1 );
266 typedef itk::MaximumImageFilter< _TSliceBinary > _TMax;
267 typename _TMax::Pointer max = _TMax::New();
268 max->SetInput( 0, extract->GetOutput( ) );
269 max->SetInput( 1, next );
271 binaryTree[ z + o ] = max->GetOutput( );
272 binaryTree[ z + o ]->DisconnectPipeline( );
277 #endif // __fpa__Common__SliceBySliceRandomWalker__hxx__