for( unsigned int k = 1; k < n; k++ )
{
- // CM Fixed a bug appearing under Windows : changed the stopping condition from <= to <.
- // Otherwise, on the last step, an element out of the range of vector Q is accessed (Q[ i + 1 ])...
+ // CM Fixed a bug appearing under Windows : changed the stopping
+ // condition from <= to <. Otherwise, on the last step, an element out
+ // of the range of vector Q is accessed (Q[ i + 1 ])...
for( unsigned int i = 0; i < n - k; i++ )
Q[ i ] = ( Q[ i ] * _1u ) + ( Q[ i + 1 ] * u );
fr[ 0 ][ 1 ] = -vT[ 1 ];
fr[ 1 ][ 1 ] = vT[ 0 ];
+ }
+ else if( TVector::Dimension == 3 )
+ {
+ this->_UpdateDerivative( );
+ this->m_Derivative->_UpdateDerivative( );
+ TVector vT = this->m_Derivative->Evaluate( u );
+ TScalar nvT = vT.GetNorm( );
+ if( nvT > TScalar( 0 ) )
+ {
+ vT /= nvT;
+ TVector vN = this->m_Derivative->m_Derivative->Evaluate( u );
+ TScalar nvN = vN.GetNorm( );
+ if( nvT > TScalar( 0 ) )
+ {
+ vN /= nvN;
+ TVector vB;
+ vB[ 0 ] = ( vT[ 1 ] * vN[ 2 ] ) - ( vT[ 2 ] * vN[ 1 ] );
+ vB[ 1 ] = ( vT[ 2 ] * vN[ 0 ] ) - ( vT[ 0 ] * vN[ 2 ] );
+ vB[ 2 ] = ( vT[ 0 ] * vN[ 1 ] ) - ( vT[ 1 ] * vN[ 0 ] );
+
+ for( unsigned int d = 0; d < 3; d++ )
+ {
+ fr[ d ][ 0 ] = vT[ d ];
+ fr[ d ][ 1 ] = vN[ d ];
+ fr[ d ][ 2 ] = vB[ d ];
+
+ } // rof
+ }
+ else
+ std::cerr << "ERROR normal" << std::endl;
+ }
+ else
+ std::cerr << "ERROR tangent" << std::endl;
} // fi
-
- /* TODO
- if( TVector::Dimension == 3 )
- {
- this->_UpdateDerivative( );
- this->m_Derivative->_UpdateDerivative( );
-
- TVector vT = this->m_Derivative->Evaluate( u );
- TVector vN = this->m_Derivative->m_Derivative->Evaluate( u );
- TScalar nvT = vT.GetNorm( );
- TScalar nvN = vN.GetNorm( );
-
- if( nvT > TScalar( 0 ) && nvN > TScalar( 0 ) )
- {
- vT /= nvT;
- vN /= nvN;
-
- TVector vB;
- vB[ 0 ] = ( vT[ 1 ] * vN[ 2 ] ) - ( vT[ 2 ] * vN[ 1 ] );
- vB[ 1 ] = ( vT[ 2 ] * vN[ 0 ] ) - ( vT[ 0 ] * vN[ 2 ] );
- vB[ 2 ] = ( vT[ 0 ] * vN[ 1 ] ) - ( vT[ 1 ] * vN[ 0 ] );
-
- TScalar nvB = vB.GetNorm( );
- if( nvB > TScalar( 0 ) )
- {
- vB /= nvB;
-
- // WARNING: Paranoiac test
- vN[ 0 ] = ( vB[ 1 ] * vT[ 2 ] ) - ( vB[ 2 ] * vT[ 1 ] );
- vN[ 1 ] = ( vB[ 2 ] * vT[ 0 ] ) - ( vB[ 0 ] * vT[ 2 ] );
- vN[ 2 ] = ( vB[ 0 ] * vT[ 1 ] ) - ( vB[ 1 ] * vT[ 0 ] );
-
- for( unsigned int d = 0; d < 3; d++ )
- {
- fr[ d ][ 0 ] = vT[ d ];
- fr[ d ][ 1 ] = vN[ d ];
- fr[ d ][ 2 ] = vB[ d ];
-
- } // rof
- }
- else
- {
- // WARNING: Trick to avoid numerical instabilities
- // in straight lines.
- typedef itk::Vector< TScalar, 3 > _TVector3;
- typedef ext::VectorToFrameFunction< _TVector3 > _TFunction;
-
- _TVector3 vT3;
- vT3[ 0 ] = vT[ 0 ];
- vT3[ 1 ] = vT[ 1 ];
- vT3[ 2 ] = vT[ 2 ];
-
- typename _TFunction::Pointer fun = _TFunction::New( );
- typename _TFunction::TFrame ffr = fun->Evaluate( vT3 );
-
- fr[ 0 ][ 0 ] = ffr[ 0 ][ 0 ];
- fr[ 0 ][ 1 ] = ffr[ 0 ][ 1 ];
- fr[ 0 ][ 2 ] = ffr[ 0 ][ 2 ];
-
- fr[ 1 ][ 0 ] = ffr[ 1 ][ 0 ];
- fr[ 1 ][ 1 ] = ffr[ 1 ][ 1 ];
- fr[ 1 ][ 2 ] = ffr[ 1 ][ 2 ];
-
- fr[ 2 ][ 0 ] = ffr[ 2 ][ 0 ];
- fr[ 2 ][ 1 ] = ffr[ 2 ][ 1 ];
- fr[ 2 ][ 2 ] = ffr[ 2 ][ 2 ];
-
- } // fi
-
- } // fi
-
- } // fi
- */
return( fr );
}
--- /dev/null
+// -------------------------------------------------------------------------
+// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
+// -------------------------------------------------------------------------
+
+#ifndef __CPPLUGINS__EXTENSIONS__ALGORITHMS__ISOIMAGESLICER__H__
+#define __CPPLUGINS__EXTENSIONS__ALGORITHMS__ISOIMAGESLICER__H__
+
+#include <itkAffineTransform.h>
+#include <itkExtractImageFilter.h>
+#include <itkImage.h>
+#include <itkImageToImageFilter.h>
+#include <itkInterpolateImageFunction.h>
+#include <itkResampleImageFilter.h>
+#include <itkVectorResampleImageFilter.h>
+#include <itkVectorInterpolateImageFunction.h>
+
+namespace cpPlugins
+{
+ namespace Extensions
+ {
+ namespace Algorithms
+ {
+ /**
+ */
+ template< class R, class I >
+ class BaseImageSlicer
+ : public itk::ImageToImageFilter< typename R::InputImageType, itk::Image< typename R::InputImageType::PixelType, R::ImageDimension - 1 > >
+ {
+ public:
+ // Basic types
+ typedef BaseImageSlicer Self;
+ typedef R TSlicer;
+ typedef I TInterpolateFunction;
+ typedef typename R::InputImageType TImage;
+ typedef typename I::CoordRepType TScalar;
+ typedef typename TImage::PixelType TPixel;
+ enum
+ {
+ Dim = TImage::ImageDimension,
+ SliceDim = TImage::ImageDimension - 1
+ };
+ typedef itk::Image< TPixel, Self::SliceDim > TSliceImage;
+
+ // itk types
+ typedef itk::ImageToImageFilter< TImage, TSliceImage > Superclass;
+ typedef itk::SmartPointer< Self > Pointer;
+ typedef itk::SmartPointer< const Self > ConstPointer;
+
+ // Internal filters
+ typedef itk::ExtractImageFilter< TImage, TSliceImage > TCollapsor;
+
+ // Various types
+ typedef typename TImage::IndexType TIndex;
+ typedef typename TImage::RegionType TRegion;
+ typedef typename TImage::SizeType TSize;
+ typedef typename TImage::SpacingType TSpacing;
+ typedef typename TSpacing::ValueType TSpacingValue;
+
+ typedef itk::AffineTransform< TScalar, Self::Dim > TTransform;
+ typedef typename TTransform::MatrixType TMatrix;
+ typedef typename TTransform::OffsetType TVector;
+
+ public:
+ itkNewMacro( Self );
+ itkTypeMacro( BaseImageSlicer, itkImageToImageFilter );
+
+ itkBooleanMacro( SizeFromMaximum );
+ itkBooleanMacro( SizeFromMinimum );
+ itkBooleanMacro( SpacingFromMaximum );
+ itkBooleanMacro( SpacingFromMinimum );
+
+ itkGetConstObjectMacro( Transform, TTransform );
+ itkGetConstMacro( DefaultValue, TPixel );
+ itkGetConstMacro( Size, TVector );
+ itkGetConstMacro( SizeFromMaximum, bool );
+ itkGetConstMacro( SizeFromMinimum, bool );
+ itkGetConstMacro( Spacing, TSpacingValue );
+ itkGetConstMacro( SpacingFromMaximum, bool );
+ itkGetConstMacro( SpacingFromMinimum, bool );
+
+ itkSetObjectMacro( Transform, TTransform );
+ itkSetMacro( Size, TVector );
+ itkSetMacro( DefaultValue, TPixel );
+ itkSetMacro( SizeFromMaximum, bool );
+ itkSetMacro( SizeFromMinimum, bool );
+ itkSetMacro( Spacing, TSpacingValue );
+ itkSetMacro( SpacingFromMaximum, bool );
+ itkSetMacro( SpacingFromMinimum, bool );
+
+ public:
+ virtual unsigned long GetMTime( ) const;
+
+ const TInterpolateFunction* GetInterpolator( ) const;
+ const TMatrix& GetRotation( ) const;
+ const TVector& GetTranslation( ) const;
+
+ void SetInterpolator( TInterpolateFunction* f );
+
+ template< class M >
+ void SetRotation( const M& r );
+
+ template< class V >
+ void SetTranslation( const V& t );
+ void SetSize( TScalar s );
+
+ protected:
+ BaseImageSlicer( );
+ virtual ~BaseImageSlicer( );
+
+ virtual void GenerateOutputInformation( ); // TODO { }
+ virtual void GenerateInputRequestedRegion( );
+ virtual void GenerateData( );
+
+ private:
+ // Purposely not implemented
+ BaseImageSlicer( const Self& );
+ void operator=( const Self& );
+
+ protected:
+ typename TSlicer::Pointer m_Slicer;
+ typename TCollapsor::Pointer m_Collapsor;
+ typename TTransform::Pointer m_Transform;
+
+ TPixel m_DefaultValue;
+
+ TVector m_Size;
+ bool m_SizeFromMaximum;
+ bool m_SizeFromMinimum;
+
+ TSpacingValue m_Spacing;
+ bool m_SpacingFromMaximum;
+ bool m_SpacingFromMinimum;
+ };
+
+ } // ecapseman
+
+ } // ecapseman
+
+} // ecapseman
+
+// -------------------------------------------------------------------------
+#define CPPLUGINS_DEFINE_ISOIMAGESLICER( name, R, F ) \
+ template< class I, class S = double > \
+ class name \
+ : public BaseImageSlicer< R< I, I, S >, F< I, S > > \
+ { \
+ public: \
+ typedef BaseImageSlicer< R< I, I, S >, F< I, S > > Superclass; \
+ typedef name Self; \
+ typedef itk::SmartPointer< Self > Pointer; \
+ typedef itk::SmartPointer< const Self > ConstPointer; \
+ public: \
+ itkNewMacro( Self ); \
+ itkTypeMacro( name, BaseSlicer ); \
+ protected: \
+ name( ) : Superclass( ) { } \
+ virtual ~name( ) { } \
+ private: \
+ name( const Self& ); \
+ void operator=( const Self& ); \
+ };
+
+namespace cpPlugins
+{
+ namespace Extensions
+ {
+ namespace Algorithms
+ {
+ CPPLUGINS_DEFINE_ISOIMAGESLICER(
+ IsoImageSlicer,
+ itk::ResampleImageFilter,
+ itk::InterpolateImageFunction
+ );
+ CPPLUGINS_DEFINE_ISOIMAGESLICER(
+ VectorIsoImageSlicer,
+ itk::VectorResampleImageFilter,
+ itk::VectorInterpolateImageFunction
+ );
+ } // ecapseman
+
+ } // ecapseman
+
+} // ecapseman
+
+#include <cpPlugins/Extensions/Algorithms/IsoImageSlicer.hxx>
+
+#endif // __CPPLUGINS__EXTENSIONS__ALGORITHMS__ISOIMAGESLICER__H__
+
+// eof - $RCSfile$
--- /dev/null
+// -------------------------------------------------------------------------
+// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
+// -------------------------------------------------------------------------
+
+#ifndef __CPPLUGINS__EXTENSIONS__ALGORITHMS__ISOIMAGESLICER__HXX__
+#define __CPPLUGINS__EXTENSIONS__ALGORITHMS__ISOIMAGESLICER__HXX__
+
+// -------------------------------------------------------------------------
+template< class R, class I >
+unsigned long cpPlugins::Extensions::Algorithms::BaseImageSlicer< R, I >::
+GetMTime( ) const
+{
+ unsigned long t = this->Superclass::GetMTime( );
+ unsigned long sT = this->m_Slicer->GetMTime( );
+ unsigned long cT = this->m_Collapsor->GetMTime( );
+ unsigned long tT = this->m_Transform->GetMTime( );
+ t = ( sT > t )? sT: t;
+ t = ( cT > t )? cT: t;
+ t = ( tT > t )? tT: t;
+ return( t );
+}
+
+// -------------------------------------------------------------------------
+template< class R, class I >
+const typename cpPlugins::Extensions::Algorithms::BaseImageSlicer< R, I >::
+TInterpolateFunction*
+cpPlugins::Extensions::Algorithms::BaseImageSlicer< R, I >::
+GetInterpolator( ) const
+{
+ return( this->m_Slicer->GetInterpolator( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class R, class I >
+const typename cpPlugins::Extensions::Algorithms::BaseImageSlicer< R, I >::
+TMatrix& cpPlugins::Extensions::Algorithms::BaseImageSlicer< R, I >::
+GetRotation( ) const
+{
+ return( this->m_Transform->GetMatrix( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class R, class I >
+const typename cpPlugins::Extensions::Algorithms::BaseImageSlicer< R, I >::
+TVector& cpPlugins::Extensions::Algorithms::BaseImageSlicer< R, I >::
+GetTranslation( ) const
+{
+ return( this->m_Transform->GetOffset( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class R, class I >
+void cpPlugins::Extensions::Algorithms::BaseImageSlicer< R, I >::
+SetInterpolator( TInterpolateFunction* f )
+{
+ this->m_Slicer->SetInterpolator( f );
+ this->Modified( );
+}
+
+// -------------------------------------------------------------------------
+template< class R, class I >
+template< class M >
+void cpPlugins::Extensions::Algorithms::BaseImageSlicer< R, I >::
+SetRotation( const M& r )
+{
+ TMatrix rotation;
+ for( unsigned int i = 0; i < Self::Dim; ++i )
+ for( unsigned int j = 0; j < Self::Dim; ++j )
+ rotation[ i ][ j ] = r[ i ][ j ];
+ this->m_Transform->SetMatrix( rotation );
+ this->Modified( );
+}
+
+// -------------------------------------------------------------------------
+template< class R, class I >
+template< class V >
+void cpPlugins::Extensions::Algorithms::BaseImageSlicer< R, I >::
+SetTranslation( const V& t )
+{
+ TVector off;
+ for( unsigned int i = 0; i < Self::Dim; ++i )
+ off[ i ] = t[ i ];
+ this->m_Transform->SetOffset( off );
+ this->Modified( );
+}
+
+// -------------------------------------------------------------------------
+template< class R, class I >
+void cpPlugins::Extensions::Algorithms::BaseImageSlicer< R, I >::
+SetSize( TScalar s )
+{
+ this->m_Size.Fill( s );
+ this->Modified( );
+}
+
+// -------------------------------------------------------------------------
+template< class R, class I >
+cpPlugins::Extensions::Algorithms::BaseImageSlicer< R, I >::
+BaseImageSlicer( )
+ : Superclass( ),
+ m_SizeFromMaximum( false ),
+ m_SizeFromMinimum( false ),
+ m_Spacing( TSpacingValue( 1 ) ),
+ m_SpacingFromMaximum( false ),
+ m_SpacingFromMinimum( false )
+{
+ this->m_Size.Fill( TScalar( 1 ) );
+
+ // Slicer
+ this->m_Slicer = TSlicer::New( );
+
+ TIndex idx;
+ idx.Fill( 0 );
+ this->m_Slicer->SetOutputStartIndex( idx );
+
+ // Dimension collapsor
+ this->m_Collapsor = TCollapsor::New( );
+ this->m_Collapsor->SetInput( this->m_Slicer->GetOutput( ) );
+
+ this->m_Transform = TTransform::New( );
+ this->m_Transform->SetIdentity( );
+}
+
+// -------------------------------------------------------------------------
+template< class R, class I >
+cpPlugins::Extensions::Algorithms::BaseImageSlicer< R, I >::
+~BaseImageSlicer( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class R, class I >
+void cpPlugins::Extensions::Algorithms::BaseImageSlicer< R, I >::
+GenerateOutputInformation( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class R, class I >
+void cpPlugins::Extensions::Algorithms::BaseImageSlicer< R, I >::
+GenerateInputRequestedRegion( )
+{
+ TImage* input = const_cast< TImage* >( this->GetInput( ) );
+ if( input != NULL )
+ input->SetRequestedRegionToLargestPossibleRegion( );
+}
+
+// -------------------------------------------------------------------------
+template< class R, class I >
+void cpPlugins::Extensions::Algorithms::BaseImageSlicer< R, I >::
+GenerateData( )
+{
+ const TImage* input = this->GetInput( );
+
+ // Spacing
+ TSpacing spac;
+ if( this->m_SpacingFromMaximum || this->m_SpacingFromMinimum )
+ {
+ spac = input->GetSpacing( );
+ TSpacingValue minIso = spac[ 0 ];
+ TSpacingValue maxIso = spac[ 0 ];
+ for( unsigned int i = 1; i < Self::Dim; i++ )
+ {
+ minIso = ( spac[ i ] < minIso )? spac[ i ]: minIso;
+ maxIso = ( spac[ i ] > maxIso )? spac[ i ]: maxIso;
+
+ } // rof
+ this->m_Spacing = ( this->m_SpacingFromMinimum )? minIso: maxIso;
+
+ } // fi
+ spac.Fill( this->m_Spacing );
+
+ // Size and origin
+ if( this->m_SizeFromMaximum || this->m_SizeFromMinimum )
+ {
+ TSize iSize = input->GetRequestedRegion( ).GetSize( );
+ TSpacing iSpac = input->GetSpacing( );
+ TScalar minSize = TScalar( iSize[ 0 ] ) * TScalar( iSpac[ 0 ] );
+ TScalar maxSize = minSize;
+ for( unsigned int i = 1; i < Self::Dim; i++ )
+ {
+ TScalar v = TScalar( iSize[ i ] ) * TScalar( iSpac[ i ] );
+ minSize = ( v < minSize )? v: minSize;
+ maxSize = ( v > maxSize )? v: maxSize;
+
+ } // rof
+ if( this->m_SizeFromMaximum )
+ this->m_Size.Fill( maxSize );
+ else
+ this->m_Size.Fill( minSize );
+
+ } // fi
+
+ TSize size;
+ typename TSlicer::OriginPointType origin;
+ size[ 0 ] = 1;
+ origin[ 0 ] = 0;
+ for( unsigned int i = 1; i < Self::Dim; i++ )
+ {
+ double s = double( this->m_Size[ i ] ) / double( spac[ i ] );
+ size[ i ] = ( unsigned int )( s );
+ origin[ i ] = -( 0.5 * this->m_Size[ i ] );
+
+ } // rof
+
+ // Prepare slicer
+ this->m_Slicer->SetInput( input );
+ this->m_Slicer->SetTransform( this->m_Transform );
+ this->m_Slicer->SetOutputSpacing( spac );
+ this->m_Slicer->SetOutputOrigin( origin );
+ this->m_Slicer->SetSize( size );
+ this->m_Slicer->SetDefaultPixelValue( this->m_DefaultValue );
+
+ // Slice!
+ // Note: UpdateLargestPossibleRegion( ) is used since we need the
+ // output regions to be updated at each filter call.
+ this->m_Slicer->UpdateLargestPossibleRegion( );
+
+ // Collapse result
+ TRegion region = this->m_Slicer->GetOutput( )->GetRequestedRegion( );
+ TSize regionSize = region.GetSize( );
+ regionSize[ 0 ] = 0;
+ region.SetSize( regionSize );
+ this->m_Collapsor->SetExtractionRegion( region );
+
+ this->m_Collapsor->GraftOutput( this->GetOutput( ) );
+ this->m_Collapsor->UpdateLargestPossibleRegion( );
+ this->GraftOutput( this->m_Collapsor->GetOutput( ) );
+}
+
+#endif // __CPPLUGINS__EXTENSIONS__ALGORITHMS__ISOIMAGESLICER__HXX__
+
+// eof - $RCSfile$
git://git.creatis.insa-lyon.fr / cpPlugins.git / commitdiff