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