1 // =========================================================================
2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
3 // =========================================================================
4 #ifndef __RandomWalkSegmentation__hxx__
5 #define __RandomWalkSegmentation__hxx__
7 #include <itkBinaryThresholdImageFilter.h>
8 #include <itkSmoothingRecursiveGaussianImageFilter.h>
9 #include <ivq/ITK/RegionOfInterestWithPaddingImageFilter.h>
10 #include <fpa/Filters/Image/ExtractAxis.h>
11 #include <fpa/Filters/Image/RandomWalker.h>
12 #include "DijkstraWithMeanAndVariance.h"
13 #include "RandomWalkLabelling.h"
16 #include <itkImageRegionConstIterator.h>
17 #include <itkImageRegionIterator.h>
20 #include <fpa/Functors/Dijkstra/Image/Gaussian.h>
22 #include <itkImageFileWriter.h>
25 // -------------------------------------------------------------------------
26 template< class _TInputImage, class _TOutputImage >
27 void RandomWalkSegmentation< _TInputImage, _TOutputImage >::
28 AddSeed( const TIndex& s )
33 this->m_Seeds.push_back( seed );
37 // -------------------------------------------------------------------------
38 template< class _TInputImage, class _TOutputImage >
39 void RandomWalkSegmentation< _TInputImage, _TOutputImage >::
40 AddSeed( const TPoint& s )
45 this->m_Seeds.push_back( seed );
49 // -------------------------------------------------------------------------
50 template< class _TInputImage, class _TOutputImage >
51 void RandomWalkSegmentation< _TInputImage, _TOutputImage >::
54 this->m_Seeds.clear( );
58 // -------------------------------------------------------------------------
59 template< class _TInputImage, class _TOutputImage >
60 unsigned long RandomWalkSegmentation< _TInputImage, _TOutputImage >::
61 GetNumberOfSeeds( ) const
63 return( this->m_Seeds.size( ) );
66 // -------------------------------------------------------------------------
67 template< class _TInputImage, class _TOutputImage >
68 void RandomWalkSegmentation< _TInputImage, _TOutputImage >::
69 SetStartSeed( const TIndex& s )
71 this->m_StartSeed.Index = s;
72 this->m_StartSeed.IsPoint = false;
76 // -------------------------------------------------------------------------
77 template< class _TInputImage, class _TOutputImage >
78 void RandomWalkSegmentation< _TInputImage, _TOutputImage >::
79 SetStartSeed( const TPoint& s )
81 this->m_StartSeed.Point = s;
82 this->m_StartSeed.IsPoint = true;
86 // -------------------------------------------------------------------------
87 template< class _TInputImage, class _TOutputImage >
88 void RandomWalkSegmentation< _TInputImage, _TOutputImage >::
89 SetEndSeed( const TIndex& s )
91 this->m_EndSeed.Index = s;
92 this->m_EndSeed.IsPoint = false;
96 // -------------------------------------------------------------------------
97 template< class _TInputImage, class _TOutputImage >
98 void RandomWalkSegmentation< _TInputImage, _TOutputImage >::
99 SetEndSeed( const TPoint& s )
101 this->m_EndSeed.Point = s;
102 this->m_EndSeed.IsPoint = true;
106 // -------------------------------------------------------------------------
107 template< class _TInputImage, class _TOutputImage >
108 RandomWalkSegmentation< _TInputImage, _TOutputImage >::
109 RandomWalkSegmentation( )
111 m_Beta( double( 20 ) ),
112 m_Sigma( double( 10 ) ),
113 m_Radius( double( 5 ) )
115 fpaFilterOutputConfigureMacro( OutputAxis, TPath );
118 // -------------------------------------------------------------------------
119 template< class _TInputImage, class _TOutputImage >
120 RandomWalkSegmentation< _TInputImage, _TOutputImage >::
121 ~RandomWalkSegmentation( )
125 // -------------------------------------------------------------------------
126 template< class _TInputImage, class _TOutputImage >
127 void RandomWalkSegmentation< _TInputImage, _TOutputImage >::
130 typedef unsigned char _TLabel;
131 typedef typename TOutputImage::PixelType _TScalar;
132 typedef itk::Image< _TLabel, TInputImage::ImageDimension > _TLabels;
134 // Prepare initial seeds
135 const TInputImage* input = this->GetInput( );
136 this->_SynchSeed( input, this->m_StartSeed );
137 this->_SynchSeed( input, this->m_EndSeed );
138 for( TSeed& seed: this->m_Seeds )
139 this->_SynchSeed( input, seed );
142 typename TOutputImage::Pointer smooth_in;
143 this->_Smooth( this->GetInput( ), smooth_in );
145 // Initial segmentation
146 typename TOutputImage::Pointer init_seg;
147 typename TPath::Pointer init_axis;
148 _TScalar init_mean, init_std;
149 typename TInputImage::RegionType roi =
150 this->_RawSegmentation(
151 smooth_in.GetPointer( ), this->m_Seeds,
152 this->m_StartSeed.Index, this->m_EndSeed.Index,
154 init_seg, init_axis, init_mean, init_std
156 init_std *= _TScalar( this->m_Sigma );
158 // Extract input ROIs
159 typename TOutputImage::Pointer smooth_in_roi, init_seg_roi;
160 roi = this->_ROI( smooth_in.GetPointer( ), roi, 10, smooth_in_roi );
161 this->_ROI( init_seg.GetPointer( ), roi, 0, init_seg_roi );
162 typename TPath::Pointer init_axis_roi = TPath::New( );
163 init_axis_roi->SetReferenceImage( smooth_in_roi.GetPointer( ) );
164 this->_AxisROI( init_axis.GetPointer( ), roi, init_axis_roi );
167 typename _TLabels::Pointer init_labels;
168 _TScalar radius = _TScalar( this->m_Radius );
170 smooth_in_roi.GetPointer( ),
171 init_seg_roi.GetPointer( ),
172 init_axis_roi.GetPointer( ),
173 init_mean, init_std, radius,
178 typename _TLabels::Pointer rw_seg;
180 smooth_in_roi.GetPointer( ),
181 init_labels.GetPointer( ),
182 init_std / _TScalar( 2 ),
187 typename TOutputImage::Pointer out_dist;
188 typename TPath::Pointer out_axis;
189 this->_DistanceAndAxis(
190 rw_seg.GetPointer( ),
191 this->m_Seeds, this->m_StartSeed, this->m_EndSeed,
195 // Put everything back to requested region
197 std::cout << "6" << std::endl;
199 TOutputImage* output = this->GetOutput( );
200 output->SetBufferedRegion( output->GetRequestedRegion( ) );
202 output->FillBuffer( -std::numeric_limits< _TScalar >::max( ) );
204 itk::ImageRegionConstIterator< TOutputImage > rIt(
205 output_roi, output_roi->GetRequestedRegion( )
207 itk::ImageRegionIterator< TOutputImage > oIt( output, roi );
210 for( ; !rIt.IsAtEnd( ); ++rIt, ++oIt )
211 oIt.Set( rIt.Get( ) );
213 TPath* output_axis = this->GetOutputAxis( );
214 output_axis->SetReferenceImage( output );
215 for( unsigned long i = 0; i < output_axis_roi->GetSize( ); ++i )
217 TIndex v = output_axis_roi->GetVertex( i );
218 for( unsigned int d = 0; d < TInputImage::ImageDimension; ++d )
219 v[ d ] += roi.GetIndex( )[ d ];
220 output_axis->AddVertex( v );
225 std::cout << "7" << std::endl;
229 // -------------------------------------------------------------------------
230 template< class _TInputImage, class _TOutputImage >
231 template< class _TIn, class _TSeed >
232 void RandomWalkSegmentation< _TInputImage, _TOutputImage >::
233 _SynchSeed( const _TIn* in, _TSeed& seed )
236 in->TransformPhysicalPointToIndex( seed.Point, seed.Index );
238 in->TransformIndexToPhysicalPoint( seed.Index, seed.Point );
242 // -------------------------------------------------------------------------
243 template< class _TInputImage, class _TOutputImage >
244 template< class _TIn, class _TOutPtr >
245 void RandomWalkSegmentation< _TInputImage, _TOutputImage >::
246 _Smooth( const _TIn* in, _TOutPtr& out )
248 typedef typename _TOutPtr::ObjectType _TOut;
249 typedef itk::SmoothingRecursiveGaussianImageFilter< _TIn, _TOut > _TSmooth;
251 typename _TSmooth::Pointer smooth = _TSmooth::New( );
252 smooth->SetInput( in );
253 smooth->SetNormalizeAcrossScale( true );
254 smooth->SetSigmaArray( in->GetSpacing( ) * double( 2 ) );
256 out = smooth->GetOutput( );
257 out->DisconnectPipeline( );
260 // -------------------------------------------------------------------------
261 template< class _TInputImage, class _TOutputImage >
262 template< class _TIn, class _TOutPtr, class _TAxisPtr, class _TSeeds >
263 typename _TIn::RegionType
264 RandomWalkSegmentation< _TInputImage, _TOutputImage >::
266 const _TIn* in, const _TSeeds& seeds,
267 const typename _TIn::IndexType& s0,
268 const typename _TIn::IndexType& s1,
269 const typename _TOutPtr::ObjectType::PixelType& beta,
270 _TOutPtr& out, _TAxisPtr& out_axis,
271 typename _TOutPtr::ObjectType::PixelType& oMean,
272 typename _TOutPtr::ObjectType::PixelType& oSTD
275 typedef typename _TOutPtr::ObjectType _TOut;
276 typedef typename _TOut::PixelType _TScalar;
277 typedef DijkstraWithMeanAndVariance< _TIn, _TOut > _TInit;
278 typedef fpa::Functors::Dijkstra::Image::Gaussian< _TIn, _TScalar > _TFun;
280 typename _TFun::Pointer fun = _TFun::New( );
281 fun->SetBeta( beta );
283 typename _TInit::Pointer init = _TInit::New( );
284 init->SetInput( in );
285 init->SetWeightFunction( fun );
286 init->StopAtOneFrontOn( );
287 for( typename _TSeeds::value_type seed: seeds )
288 init->AddSeed( seed.Point );
291 // Get initial values
292 oMean = _TScalar( init->GetMean( ) );
293 oSTD = _TScalar( init->GetDeviation( ) );
295 // Get initial objects
296 init->GetMinimumSpanningTree( )->GetPath( out_axis, s0, s1 );
297 out = init->GetOutput( );
298 out->DisconnectPipeline( );
301 typename _TIn::IndexType min_idx = init->GetMinVertex( );
302 typename _TIn::IndexType max_idx = init->GetMaxVertex( );
303 typename _TIn::SizeType roi_size;
304 for( unsigned int i = 0; i < TInputImage::ImageDimension; ++i )
305 roi_size[ i ] = max_idx[ i ] - min_idx[ i ] + 1;
306 typename _TIn::RegionType roi;
307 roi.SetIndex( min_idx );
308 roi.SetSize( roi_size );
312 // -------------------------------------------------------------------------
313 template< class _TInputImage, class _TOutputImage >
314 template< class _TIn, class _TOutPtr >
315 typename _TIn::RegionType
316 RandomWalkSegmentation< _TInputImage, _TOutputImage >::
318 const _TIn* in, const typename _TIn::RegionType& roi, unsigned int pad,
322 typedef typename _TOutPtr::ObjectType _TOut;
323 typedef ivq::ITK::RegionOfInterestWithPaddingImageFilter< _TIn, _TOut > _TROI;
325 typename _TROI::Pointer filter = _TROI::New( );
326 filter->SetInput( in );
327 filter->SetRegionOfInterest( roi );
328 filter->SetPadding( pad );
330 out = filter->GetOutput( );
331 out->DisconnectPipeline( );
332 return( filter->GetRegionOfInterest( ) );
335 // -------------------------------------------------------------------------
336 template< class _TInputImage, class _TOutputImage >
337 template< class _TAxisPtr, class _TRegion >
338 void RandomWalkSegmentation< _TInputImage, _TOutputImage >::
340 const typename _TAxisPtr::ObjectType* in, const _TRegion& roi,
344 typedef typename _TAxisPtr::ObjectType _TAxis;
346 for( unsigned long i = 0; i < in->GetSize( ); ++i )
348 typename _TAxis::TIndex v = in->GetVertex( i );
349 for( unsigned int d = 0; d < _TAxis::Dimension; ++d )
350 v[ d ] -= roi.GetIndex( )[ d ];
356 // -------------------------------------------------------------------------
357 template< class _TInputImage, class _TOutputImage >
358 template< class _TInRaw, class _TInCosts, class _TAxis, class _TScalar, class _TOutPtr >
359 void RandomWalkSegmentation< _TInputImage, _TOutputImage >::
361 const _TInRaw* raw, const _TInCosts* costs, const _TAxis* axis,
362 const _TScalar& mean, const _TScalar& dev, const _TScalar& radius,
366 typedef typename _TOutPtr::ObjectType _TOut;
367 typedef RandomWalkLabelling< _TInRaw, _TInCosts, _TOut > _TLabel;
369 typename _TLabel::Pointer label = _TLabel::New( );
370 label->SetInputImage( raw );
371 label->SetInputCosts( costs );
372 label->SetInputPath( axis );
373 label->SetInsideLabel( 1 );
374 label->SetOutsideLabel( 2 );
375 label->SetLowerLabel( 3 );
376 label->SetUpperLabel( 4 );
377 label->SetRadius( radius );
378 label->SetLowerThreshold( mean - dev );
379 label->SetUpperThreshold( mean + dev );
381 out = label->GetOutputLabels( );
382 out->DisconnectPipeline( );
385 // -------------------------------------------------------------------------
386 template< class _TInputImage, class _TOutputImage >
387 template< class _TIn, class _TOutPtr >
388 void RandomWalkSegmentation< _TInputImage, _TOutputImage >::
390 const _TIn* in, const typename _TOutPtr::ObjectType* labels,
391 const typename _TIn::PixelType& beta,
395 typedef typename _TIn::PixelType _TScalar;
396 typedef typename _TOutPtr::ObjectType _TOut;
397 typedef fpa::Functors::Dijkstra::Image::Gaussian< _TIn, _TScalar > _TFun;
398 typedef fpa::Filters::Image::RandomWalker< _TIn, _TOut, _TScalar > _TRandomWalker;
399 typedef itk::BinaryThresholdImageFilter< _TOut, _TOut > _TExtract;
401 typename _TFun::Pointer fun = _TFun::New( );
402 fun->SetBeta( beta );
404 typename _TRandomWalker::Pointer rw = _TRandomWalker::New( );
405 rw->SetInputImage( in );
406 rw->SetInputLabels( labels );
407 rw->SetWeightFunction( fun );
409 typename _TExtract::Pointer extract = _TExtract::New( );
410 extract->SetInput( rw->GetOutputLabels( ) );
411 extract->SetInsideValue( 1 );
412 extract->SetOutsideValue( 0 );
413 extract->SetLowerThreshold( 1 );
414 extract->SetUpperThreshold( 1 );
416 out = extract->GetOutput( );
417 out->DisconnectPipeline( );
420 // -------------------------------------------------------------------------
421 template< class _TInputImage, class _TOutputImage >
422 template< class _TIn, class _TSeeds, class _TSeed, class _TOutPtr, class _TAxisPtr >
423 void RandomWalkSegmentation< _TInputImage, _TOutputImage >::
425 const _TIn* in, const _TSeeds& seeds,
426 const _TSeed& p0, const _TSeed& p1,
427 _TOutPtr& out_dist, _TAxisPtr& out_axis
430 typedef typename _TOutPtr::ObjectType _TOut;
431 typedef typename _TOut::PixelType _TScalar;
432 typedef fpa::Filters::Image::ExtractAxis< _TIn, _TScalar > _TExtract;
434 // Prepare output values
435 typename _TIn::IndexType s0, s1;
436 in->TransformPhysicalPointToIndex( p0.Point, s0 );
437 in->TransformPhysicalPointToIndex( p1.Point, s1 );
439 typename _TExtract::Pointer extract = _TExtract::New( );
440 extract->SetInput( const_cast< _TIn* >( in ) );
441 extract->AddSeed( s0 );
442 extract->AddSeed( s1 );
443 for( typename _TSeeds::value_type seed: seeds )
444 extract->AddSeed( seed.Point );
445 extract->SetStartIndex( s0 );
446 extract->SetEndIndex( s1 );
448 out_axis = TPath::New( );
449 out_axis->Graft( extract->GetOutput( ) );
450 this->_Smooth( extract->GetCenterness( )->GetOutput( ), out_dist );
453 #endif // __RandomWalkSegmentation__hxx__