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 <itkExtractImageFilter.h>
10 #include <itkImageRegionConstIterator.h>
11 #include <itkImageRegionIterator.h>
12 #include <itkJoinSeriesImageFilter.h>
13 #include <itkMinimumMaximumImageCalculator.h>
14 #include <fpa/Common/RandomWalker.h>
15 #include <fpa/Functors/Dijkstra/Image/Gaussian.h>
17 // -------------------------------------------------------------------------
18 template< class _TImage, class _TLabels, class _TScalarImage >
19 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
20 SliceBySliceRandomWalker( )
23 ivqITKInputConfigureMacro( InputLabels, TLabels );
24 ivqITKInputConfigureMacro( InputVesselness, TScalarImage );
27 // -------------------------------------------------------------------------
28 template< class _TImage, class _TLabels, class _TScalarImage >
29 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
30 ~SliceBySliceRandomWalker( )
34 // -------------------------------------------------------------------------
35 template< class _TImage, class _TLabels, class _TScalarImage >
37 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
40 typedef itk::Image< TPixel, TImage::ImageDimension - 1 > _TSliceImage;
41 typedef itk::Image< TLabel, TImage::ImageDimension - 1 > _TSliceLabels;
42 typedef itk::Image< TScalar, TImage::ImageDimension - 1 > _TSliceScalarImage;
45 const TImage* input = this->GetInput( );
46 const TLabels* labels = this->GetInputLabels( );
47 const TScalarImage* vesselness = this->GetInputVesselness( );
48 this->AllocateOutputs( );
51 typename TImage::RegionType region = input->GetRequestedRegion( );
52 typename TImage::SizeType regionIndex = region.GetIndex( );
53 typename TImage::SizeType sliceSize = region.GetSize( );
54 int numSlices = sliceSize[ 2 ];
56 this->m_VesselnessValue =
57 this->m_VesselnessThreshold * this->m_VesselnessMax / TScalar( 100 );
60 typename TScalarImage::Pointer composite = TScalarImage::New( );
61 composite->CopyInformation( vesselness );
62 composite->SetBufferedRegion( vesselness->GetBufferedRegion( ) );
63 composite->Allocate( );
64 this->_Composite( composite, labels, vesselness );
67 std::vector< typename _TSliceImage::Pointer > data3D( numSlices );
68 std::vector< typename _TSliceLabels::Pointer > binaryTree( numSlices );
69 std::vector< typename _TSliceScalarImage::Pointer > fusion( numSlices );
70 for( int i = 0; i < numSlices; ++i )
72 // Prepare extraction region
73 typename TImage::IndexType sliceIndex = regionIndex;
75 typename TImage::RegionType sliceRegion( sliceIndex, sliceSize );
78 typename _TSliceImage::Pointer input_slice;
79 typename _TSliceLabels::Pointer labels_slice;
80 typename _TSliceScalarImage::Pointer composite_slice;
81 this->_Slice( input_slice, input, sliceRegion );
82 this->_Slice( labels_slice, labels, sliceRegion );
83 this->_Slice( composite_slice, composite, sliceRegion );
85 // Update vectors with each image
86 data3D[ i ] = input_slice;
87 binaryTree[ i ] = labels_slice;
88 fusion[ i ] = composite_slice;
92 // Random walker slice-by-slice
93 this->_GoDown( binaryTree, data3D, fusion, numSlices );
94 this->_GoUp( binaryTree, data3D, fusion, numSlices );
97 typedef itk::JoinSeriesImageFilter< _TSliceImage, TImage > _TJoin;
98 typename _TJoin::Pointer join = _TJoin::New( );
99 for( int i = 0; i < numSlices; ++i )
100 join->SetInput( i, binaryTree[ i ] );
102 this->GetOutput( )->Graft( join->GetOutput( ) );
105 // -------------------------------------------------------------------------
106 template< class _TImage, class _TLabels, class _TScalarImage >
108 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
110 typename TScalarImage::Pointer& composite,
111 const TLabels* labels,
112 const TScalarImage* vesselness
115 // Min-Max vesselness values
116 typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
117 typename _TMinMax::Pointer minMax = _TMinMax::New( );
118 minMax->SetImage( vesselness );
120 TScalar maxVess = minMax->GetMaximum( );
121 typename TScalarImage::RegionType region = labels->GetRequestedRegion( );
124 typename TScalarImage::Pointer composite = TScalarImage::New( );
125 composite->CopyInformation( vesselness );
126 composite->SetBufferedRegion( vesselness->GetBufferedRegion( ) );
127 composite->Allocate( );
128 itk::ImageRegionConstIterator< TScalarImage > v( vesselness, region );
129 itk::ImageRegionConstIterator< TLabels > l( labels, region );
130 itk::ImageRegionIterator< TScalarImage > c( composite, region );
134 for( ; !v.IsAtEnd( ); ++v, ++l, ++c )
135 c.Set( ( l.Get( ) != 0 )? maxVess: v.Get( ) );
138 // -------------------------------------------------------------------------
139 template< class _TImage, class _TLabels, class _TScalarImage >
140 template< class _TSlicePtr, class _TInput >
142 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
145 const _TInput* input,
146 typename _TInput::RegionType region
149 typedef typename _TSlicePTr::ObjectType _TSlice;
150 typedef itk::ExtractImageFilter< _TInput, _TSlice > _TExtract;
152 typename _TExtract::Pointer extract = _TExtract::New( );
153 extract->SetInput( input );
154 extract->SetExtractionRegion( region );
155 extract->SetDirectionCollapseToIdentity( );
157 slice = extract->GetOutput( );
158 slice->DisconnectPipeline( );
161 // -------------------------------------------------------------------------
162 template< class _TImage, class _TLabels, class _TScalarImage >
163 template< class _TBinaryTree, class _TData3D, class _TFusion >
165 fpa::Common::SliceBySliceRandomWalker< _TImage, _TLabels, _TScalarImage >::
167 _TBinaryTree& binaryTree,
168 const _TData3D& data3D,
169 const _TFusion& fusion,
173 typedef typename _TBinaryTree::value_type _TSliceBinaryPtr;
174 typedef typename _TData3D::value_type _TSliceData3DPtr;
175 typedef typename _TFusion::value_type _TSliceFusionPtr;
176 typedef typename _TSliceBinaryPtr::ObjectType _TSliceBinary;
177 typedef typename _TSliceData3DPtr::ObjectType _TSliceData3D;
178 typedef typename _TSliceFusionPtr::ObjectType _TSliceFusion;
179 typedef typename _TSliceFusion::PixelType _TSliceScalar;
182 while( z < numSlices - 2 )
185 _TSliceBinaryPtr tmp = binaryTree[ z ];
186 _TSliceBinaryPtr next = binaryTree[ z + 1 ];
187 _TSliceData3DPtr data = data3D[ z + 1 ];
188 _TSliceFusionPtr vess = fusion[ z + 1 ];
189 typename _TSliceBinary::RegionType region = tmp->GetRequestedRegion( );
191 // Create seeds and background
192 _TSliceBinaryPtr seeds = _TSliceBinary::New( );
193 seeds->CopyInformation( tmp );
194 seeds->SetRequestedRegion( region );
196 seeds->FillBuffer( 0 );
198 itk::ImageRegionConstIterator< _TSliceBinary > tIt( tmp, region );
199 itk::ImageRegionConstIterator< _TSliceFusion > vIt( vess, region );
200 itk::ImageRegionIterator< _TSliceBinary > sIt( seeds, region );
207 for( ; !tIt.IsAtEnd( ); ++tIt, ++vIt, ++sIt )
209 if( tIt.Get( ) != 0 )
212 if( vIt.Get( ) > this->m_VesselnessValue )
221 if( vIt.Get( ) < this->m_VesselnessValue )
233 // Random walker function
234 typedef fpa::Functors::Dijkstra::Image::Gaussian< _TSliceData3D, _TSliceScalar > _TFunction;
235 typename TFunction::Pointer edge = TFunction::New( );
236 edge->SetBeta( this->m_Beta );
237 edge->SetEpsilon( this->m_Eps );
238 edge->TreatAsWeightOff( );
240 // Real random walker
241 typedef fpa::Common::RandomWalker< _TSliceData3D, _TSliceBinary, _TSliceScalar > _TRandomWalker;
242 typename _TRandomWalker::Pointer rw = _TRandomWalker::New( );
243 rw->SetInput( data );
244 rw->SetInputLabels( seeds );
245 rw->SetEdgeFunction( edge );
249 typedef itk::MaximumImageFilter<SliceType> MaxFilterType;
250 MaxFilterType::Pointer maxFilter = MaxFilterType::New();
251 maxFilter->SetInput(0, img2);
252 maxFilter->SetInput(1, next);
255 SliceType::Pointer sth = maxFilter->GetOutput();
256 ImageAlgorithm::Copy(sth.GetPointer(), binaryTree[z + 1].GetPointer(), sth->GetRequestedRegion(), binaryTree[z + 1]->GetRequestedRegion());
262 // -------------------------------------------------------------------------
263 template< class _TImage, class _TLabels, class _TScalarImage >
264 template< class _TBinaryTree, class _TData3D, class _TFusion >
267 _TBinaryTree& binaryTree,
268 const _TData3D& data3D,
269 const _TFusion& fusion,
275 #endif // __fpa__Common__SliceBySliceRandomWalker__hxx__