lib/cpPlugins/Extensions/Algorithms/BezierCurveFunction.hxx

   1 // -------------------------------------------------------------------------
   2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
   3 // -------------------------------------------------------------------------
   4
   5 #ifndef __CPPLUGINS__EXTENSIONS__ALGORITHMS__BEZIERCURVEFUNCTION__HXX__
   6 #define __CPPLUGINS__EXTENSIONS__ALGORITHMS__BEZIERCURVEFUNCTION__HXX__
   7
   8 // -------------------------------------------------------------------------
   9 template< class V >
  10 void cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
  11 AddPoint( const TVector& v )
  12 {
  13   this->m_Vectors.push_back( v );
  14   this->m_DerivativeUpdated = false;
  15   this->Modified( );
  16 }
  17
  18 // -------------------------------------------------------------------------
  19 template< class V >
  20 unsigned int cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
  21 GetNumberOfPoints( ) const
  22 {
  23   return( this->m_Vectors.size( ) );
  24 }
  25
  26 // -------------------------------------------------------------------------
  27 template< class V >
  28 typename cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
  29 TVector cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
  30 Evaluate( const TScalar& u ) const
  31 {
  32   TVectorsContainer Q = this->m_Vectors;
  33   unsigned int n = Q.size( );
  34   TScalar _1u = TScalar( 1 ) - u;
  35
  36   for( unsigned int k = 1; k < n; k++ )
  37   {
  38     // CM Fixed a bug appearing under Windows : changed the stopping condition from <= to <.
  39     // Otherwise, on the last step, an element out of the range of vector Q is accessed (Q[ i + 1 ])...
  40     for( unsigned int i = 0; i < n - k; i++ )
  41       Q[ i ] = ( Q[ i ] * _1u ) + ( Q[ i + 1 ] * u );
  42
  43   } // rof
  44   return( Q[ 0 ] );
  45 }
  46
  47 // -------------------------------------------------------------------------
  48 template< class V >
  49 typename cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
  50 TFrame cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
  51 EvaluateFrenetFrame( const TScalar& u ) const
  52 {
  53   TFrame fr;
  54   fr.Fill( TScalar( 0 ) );
  55   if( TVector::Dimension == 2 )
  56   {
  57     this->_UpdateDerivative( );
  58     this->m_Derivative->_UpdateDerivative( );
  59
  60     TVector vT = this->m_Derivative->Evaluate( u );
  61     TScalar nvT = vT.GetNorm( );
  62     if( TScalar( 0 ) < nvT )
  63       vT /= nvT;
  64
  65     fr[ 0 ][ 0 ] = vT[ 0 ];
  66     fr[ 1 ][ 0 ] = vT[ 1 ];
  67
  68     fr[ 0 ][ 1 ] = -vT[ 1 ];
  69     fr[ 1 ][ 1 ] =  vT[ 0 ];
  70
  71   } // fi
  72
  73   /* TODO
  74      if( TVector::Dimension == 3 )
  75      {
  76      this->_UpdateDerivative( );
  77      this->m_Derivative->_UpdateDerivative( );
  78
  79      TVector vT = this->m_Derivative->Evaluate( u );
  80      TVector vN = this->m_Derivative->m_Derivative->Evaluate( u );
  81      TScalar nvT = vT.GetNorm( );
  82      TScalar nvN = vN.GetNorm( );
  83
  84      if( nvT > TScalar( 0 ) && nvN > TScalar( 0 ) )
  85      {
  86      vT /= nvT;
  87      vN /= nvN;
  88
  89      TVector vB;
  90      vB[ 0 ] = ( vT[ 1 ] * vN[ 2 ] ) - ( vT[ 2 ] * vN[ 1 ] );
  91      vB[ 1 ] = ( vT[ 2 ] * vN[ 0 ] ) - ( vT[ 0 ] * vN[ 2 ] );
  92      vB[ 2 ] = ( vT[ 0 ] * vN[ 1 ] ) - ( vT[ 1 ] * vN[ 0 ] );
  93
  94      TScalar nvB = vB.GetNorm( );
  95      if( nvB > TScalar( 0 ) )
  96      {
  97      vB /= nvB;
  98
  99      // WARNING: Paranoiac test
 100      vN[ 0 ] = ( vB[ 1 ] * vT[ 2 ] ) - ( vB[ 2 ] * vT[ 1 ] );
 101      vN[ 1 ] = ( vB[ 2 ] * vT[ 0 ] ) - ( vB[ 0 ] * vT[ 2 ] );
 102      vN[ 2 ] = ( vB[ 0 ] * vT[ 1 ] ) - ( vB[ 1 ] * vT[ 0 ] );
 103
 104      for( unsigned int d = 0; d < 3; d++ )
 105      {
 106      fr[ d ][ 0 ] = vT[ d ];
 107      fr[ d ][ 1 ] = vN[ d ];
 108      fr[ d ][ 2 ] = vB[ d ];
 109
 110      } // rof
 111      }
 112      else
 113      {
 114      // WARNING: Trick to avoid numerical instabilities
 115      //          in straight lines.
 116      typedef itk::Vector< TScalar, 3 >               _TVector3;
 117      typedef ext::VectorToFrameFunction< _TVector3 > _TFunction;
 118
 119      _TVector3 vT3;
 120      vT3[ 0 ] = vT[ 0 ];
 121      vT3[ 1 ] = vT[ 1 ];
 122      vT3[ 2 ] = vT[ 2 ];
 123
 124      typename _TFunction::Pointer fun = _TFunction::New( );
 125      typename _TFunction::TFrame ffr = fun->Evaluate( vT3 );
 126
 127      fr[ 0 ][ 0 ] = ffr[ 0 ][ 0 ];
 128      fr[ 0 ][ 1 ] = ffr[ 0 ][ 1 ];
 129      fr[ 0 ][ 2 ] = ffr[ 0 ][ 2 ];
 130
 131      fr[ 1 ][ 0 ] = ffr[ 1 ][ 0 ];
 132      fr[ 1 ][ 1 ] = ffr[ 1 ][ 1 ];
 133      fr[ 1 ][ 2 ] = ffr[ 1 ][ 2 ];
 134
 135      fr[ 2 ][ 0 ] = ffr[ 2 ][ 0 ];
 136      fr[ 2 ][ 1 ] = ffr[ 2 ][ 1 ];
 137      fr[ 2 ][ 2 ] = ffr[ 2 ][ 2 ];
 138
 139      } // fi
 140
 141      } // fi
 142
 143      } // fi
 144   */
 145   return( fr );
 146 }
 147
 148 // -------------------------------------------------------------------------
 149 template< class V >
 150 typename cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
 151 TScalar cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
 152 EvaluateLength( ) const
 153 {
 154   unsigned int n = this->GetNumberOfPoints( ) << 1;
 155   TScalar d = TScalar( 0 );
 156   TVector v0 = this->Evaluate( 0 );
 157   for( unsigned int i = 1; i < n; i++ )
 158   {
 159     TVector v1 = this->Evaluate( TScalar( i ) / TScalar( n - 1 ) );
 160     d += ( v1 - v0 ).GetNorm( );
 161     v0 = v1;
 162
 163   } // rof
 164   return( d );
 165 }
 166
 167 // -------------------------------------------------------------------------
 168 template< class V >
 169 cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
 170 BezierCurveFunction( )
 171   : Superclass( ),
 172     m_DerivativeUpdated( false )
 173 {
 174 }
 175
 176 // -------------------------------------------------------------------------
 177 template< class V >
 178 void cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
 179 _UpdateDerivative( ) const
 180 {
 181   if( this->m_DerivativeUpdated )
 182     return;
 183
 184   this->m_Derivative = Self::New( );
 185   unsigned int n = this->m_Vectors.size( ) - 1;
 186   for( unsigned int i = 0; i < n; i++ )
 187     this->m_Derivative->AddPoint(
 188       TScalar( n ) * ( this->m_Vectors[ i + 1 ] - this->m_Vectors[ i ] )
 189       );
 190 }
 191
 192 #endif // __CPPLUGINS__EXTENSIONS__ALGORITHMS__BEZIERCURVEFUNCTION__HXX__
 193
 194 // eof - $RCSfile$