lib/cpExtensions/Algorithms/CPRFilter.hxx

   1 // -------------------------------------------------------------------------
   2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
   3 // -------------------------------------------------------------------------
   4
   5 #ifndef __CPEXTENSIONS__ALGORITHMS__CPRFILTER__HXX__
   6 #define __CPEXTENSIONS__ALGORITHMS__CPRFILTER__HXX__
   7
   8 #include <itkPoint.h>
   9 #include <itkMinimumMaximumImageCalculator.h>
  10 #include <itkNearestNeighborInterpolateImageFunction.h>
  11 #include <cpExtensions/Algorithms/BezierCurveFunction.h>
  12
  13 // -------------------------------------------------------------------------
  14 template< class I, class S >
  15 const typename cpExtensions::Algorithms::CPRFilter< I, S >::
  16 TAxis* cpExtensions::Algorithms::CPRFilter< I, S >::
  17 GetAxis( ) const
  18 {
  19   return( NULL );
  20 }
  21
  22 // -------------------------------------------------------------------------
  23 template< class I, class S >
  24 void cpExtensions::Algorithms::CPRFilter< I, S >::
  25 SetAxis( const TAxis* axis )
  26 {
  27 }
  28
  29 // -------------------------------------------------------------------------
  30 template< class I, class S >
  31 cpExtensions::Algorithms::CPRFilter< I, S >::
  32 CPRFilter( )
  33   : Superclass( )
  34 {
  35   this->m_Interpolator =
  36     itk::NearestNeighborInterpolateImageFunction< I, S >::New( );
  37 }
  38
  39 // -------------------------------------------------------------------------
  40 template< class I, class S >
  41 cpExtensions::Algorithms::CPRFilter< I, S >::
  42 ~CPRFilter( )
  43 {
  44 }
  45
  46 // -------------------------------------------------------------------------
  47 template< class I, class S >
  48 void cpExtensions::Algorithms::CPRFilter< I, S >::
  49 GenerateData( )
  50 {
  51   typedef itk::Point< S, VDimension >                         _P;
  52   typedef typename _P::VectorType                             _V;
  53   typedef cpExtensions::Algorithms::BezierCurveFunction< _V > _Bezier;
  54   typedef typename _Bezier::TFrame                            _M;
  55   typedef itk::MinimumMaximumImageCalculator< I >             _MinMax;
  56
  57   auto image = this->GetInput( );
  58   auto axis = this->GetAxis( );
  59
  60   // 0. Get image properties
  61   typename _MinMax::Pointer minmax = _MinMax::New( );
  62   minmax->SetImage( image );
  63   minmax->Compute( );
  64   auto background = minmax->GetMinimum( );
  65
  66   // 1. Construct Bezier curve
  67   typename _Bezier::Pointer bezier = _Bezier::New( );
  68   auto vertices = axis->GetVertexList( );
  69   for( unsigned int i = 0; i < vertices->Size( ); ++i )
  70   {
  71     auto index = vertices->GetElement( i );
  72     _P point;
  73     image->TransformContinuousIndexToPhysicalPoint( index, point );
  74     bezier->AddPoint( point.GetVectorFromOrigin( ) );
  75
  76   } // rof
  77
  78   // 2. Slice image and prepare join filter
  79   this->m_Join = TJoin::New( );
  80   unsigned int nSlices = this->m_NumberOfSlices;
  81   if( nSlices == 0 )
  82     nSlices = bezier->GetNumberOfPoints( );
  83   _V vpre;
  84   _M mpre;
  85   S len = S( 0 );
  86   for( unsigned int i = 0; i <= nSlices; ++i )
  87   {
  88     // Get geometrical data and correct Frenet frame
  89     S u = S( i ) / S( nSlices );
  90     _V vnex = bezier->Evaluate( u );
  91     _M mnex = bezier->EvaluateFrenetFrame( u );
  92     if( i > 0 )
  93     {
  94       len += ( vpre - vnex ).GetNorm( );
  95
  96       // TODO: Update Frenet frame orientation
  97
  98     } // fi
  99
 100     // Slice image
 101     typename TSlicer::Pointer slicer = TSlicer::New( );
 102     slicer->SetInput( image );
 103     slicer->SetRotation( mnex );
 104     slicer->SetTranslation( vnex );
 105     slicer->SetSize( S( this->m_SliceRadius ) );
 106     slicer->SpacingFromMinimumOn( );
 107     slicer->SetDefaultValue( background );
 108     slicer->SetInterpolator( this->m_Interpolator );
 109
 110     // Add output to join filter and keep track of the current slicer
 111     this->m_Join->SetInput( i, slicer->GetOutput( ) );
 112     this->m_Slices.push_back( slicer );
 113
 114     // Update values
 115     vpre = vnex;
 116     mpre = mnex;
 117
 118   } // rof
 119
 120   // 3. Finish join filter configuration
 121   this->m_Join->SetSpacing( double( len / S( nSlices + 1 ) ) );
 122   this->m_Join->SetOrigin( double( 0 ) );
 123
 124   // 4. Graft outputs
 125   this->m_Join->GraftOutput( this->GetOutput( ) );
 126   this->m_Join->UpdateLargestPossibleRegion( );
 127   this->GraftOutput( this->m_Join->GetOutput( ) );
 128 }
 129
 130 #endif // __CPEXTENSIONS__ALGORITHMS__CPRFILTER__HXX__
 131
 132 // eof - $RCSfile$