From 1aa0bddc396d4ecec808dacdd14ebb5e1ec19e0a Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Wed, 24 Jun 2015 18:31:03 -0500 Subject: [PATCH] Iso slicer added --- .../Algorithms/BezierCurveFunction.hxx | 111 +++------ .../Extensions/Algorithms/IsoImageSlicer.h | 189 ++++++++++++++ .../Extensions/Algorithms/IsoImageSlicer.hxx | 233 ++++++++++++++++++ 3 files changed, 458 insertions(+), 75 deletions(-) create mode 100644 lib/cpPlugins/Extensions/Algorithms/IsoImageSlicer.h create mode 100644 lib/cpPlugins/Extensions/Algorithms/IsoImageSlicer.hxx diff --git a/lib/cpPlugins/Extensions/Algorithms/BezierCurveFunction.hxx b/lib/cpPlugins/Extensions/Algorithms/BezierCurveFunction.hxx index dad1482..3aedf8f 100644 --- a/lib/cpPlugins/Extensions/Algorithms/BezierCurveFunction.hxx +++ b/lib/cpPlugins/Extensions/Algorithms/BezierCurveFunction.hxx @@ -35,8 +35,9 @@ Evaluate( const TScalar& u ) const 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 ); @@ -67,81 +68,41 @@ EvaluateFrenetFrame( const TScalar& u ) const 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 ); } diff --git a/lib/cpPlugins/Extensions/Algorithms/IsoImageSlicer.h b/lib/cpPlugins/Extensions/Algorithms/IsoImageSlicer.h new file mode 100644 index 0000000..cda4130 --- /dev/null +++ b/lib/cpPlugins/Extensions/Algorithms/IsoImageSlicer.h @@ -0,0 +1,189 @@ +// ------------------------------------------------------------------------- +// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co) +// ------------------------------------------------------------------------- + +#ifndef __CPPLUGINS__EXTENSIONS__ALGORITHMS__ISOIMAGESLICER__H__ +#define __CPPLUGINS__EXTENSIONS__ALGORITHMS__ISOIMAGESLICER__H__ + +#include +#include +#include +#include +#include +#include +#include +#include + +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 + +#endif // __CPPLUGINS__EXTENSIONS__ALGORITHMS__ISOIMAGESLICER__H__ + +// eof - $RCSfile$ diff --git a/lib/cpPlugins/Extensions/Algorithms/IsoImageSlicer.hxx b/lib/cpPlugins/Extensions/Algorithms/IsoImageSlicer.hxx new file mode 100644 index 0000000..fd90d02 --- /dev/null +++ b/lib/cpPlugins/Extensions/Algorithms/IsoImageSlicer.hxx @@ -0,0 +1,233 @@ +// ------------------------------------------------------------------------- +// @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$ -- 2.45.1