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