[Insight-users] CSwig Wrapping: itkRecursiveGaussianImageFilter
SetOrder
Luis Ibanez
luis . ibanez at kitware . com
Wed, 15 Oct 2003 10:42:01 -0400
Hua, Glen
Thanks for contributing this code.
The modifications look fine, they were just
commited to the repository.
http://www . itk . org/cgi-bin/cvsweb . cgi/Insight/Code/BasicFilters/itkRecursiveGaussianImageFilter . h?cvsroot=Insight&sortby=date
Please let us know if you find any
further problems.
Thanks
Luis
------------------
Hua Qian wrote:
> Hello All,
>
> We are having some difficulties with python wrapping
> of itkRecursiveGaussianImageFilter class, please see
> the attached conversions below.
>
> It seems that there is no way to use the SetOrder method
> of itkRecursiveGaussianImageFilter in Python/Tcl.
> To solve the problem, we ended up adding new methods
> SetOrderZero, SetOrderFirst and SetOrderSecond in the
> C++ source code (the modified .h and .txx files from Glen
> are attached). Anyone has better solution? Otherwise,
> if our solution is reasonable, could someone commit
> the change to cvs?
>
> Regards,
> Hua
>
>
>
> ------------------------------------------
>
> It works! Thanks, Hua.
>
> I've attached the modified .h and .cxx files. Can they be committed to
> the cvs?
>
> Thanks,
> Glen
>
> Hua Qian wrote:
>
>> Hi Glen,
>>
>> This is a tough one. I could not figure out any clean
>> solution so far.
>>
>> A quick and dirty solution is to modify the C++ code -;) :
>> add methods SetOrderToFirst, SetOrderToSecond, etc.
>> ... maybe, we could ask someone to commit the change to the
>> ITK cvs repository.
>>
>> Hua
>>
>>
>>
>> Glen Lehmann wrote:
>>
>>> Hi Hua,
>>>
>>> I'm starting to have some success with ITK in python. I have a
>>> question for you regarding itkRecursiveGaussianImageFilter. It
>>> defines an enumeration for the derivative order and I don't know how
>>> to access this in python.
>>>
>>> in c++;
>>>
>>> typedef itk::RecursiveGaussianImageFilter< ImageType ,ImageType>
>>> GaussianFilterType;
>>>
>>> recursiveFilter = GaussianFilterType::New();
>>> recursiveFilter->SetOrder( GaussianFilterType::FirstOrder );
>>>
>>> in python;
>>>
>>> recursiveFilter = itk.itkRecursiveGaussianImageFilterF2F2_New()
>>> recursiveFilter.SetOrder( ??? )
>>>
>>> What is the correct way to access these enumerations from python?
>>>
>>> Thanks,
>>> Glen
>>>
>>>
>>
>>
>>
>
>
>
> ------------------------------------------------------------------------
>
> /*=========================================================================
>
> Program: Insight Segmentation & Registration Toolkit
> Module: $RCSfile: itkRecursiveGaussianImageFilter.h,v $
> Language: C++
> Date: $Date: 2003/09/10 14:28:55 $
> Version: $Revision: 1.21 $
>
> Copyright (c) Insight Software Consortium. All rights reserved.
> See ITKCopyright.txt or http://www . itk . org/HTML/Copyright . htm for details.
>
> This software is distributed WITHOUT ANY WARRANTY; without even
> the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
> PURPOSE. See the above copyright notices for more information.
>
> =========================================================================*/
> #ifndef __itkRecursiveGaussianImageFilter_h
> #define __itkRecursiveGaussianImageFilter_h
>
> #include "itkRecursiveSeparableImageFilter.h"
>
> namespace itk
> {
>
> /** \class RecursiveGaussianImageFilter
> * \brief Base class for recursive convolution with Gaussian kernel.
> *
> * RecursiveGaussianImageFilter is the base class for recursive filters that
> * approximate convolution with the Gaussian kernel.
> * This class implements the recursive filtering
> * method proposed by R.Deriche in IEEE-PAMI
> * Vol.12, No.1, January 1990, pp 78-87.
> *
> * As compared to itk::DiscreteGaussianImageFilter, this filter tends
> * to be faster for large kernels, and it can take the derivative
> * of the blurred image in one step. Also, note that we have
> * itk::RecursiveGaussianImageFilter::SetSigma(), but
> * itk::DiscreteGaussianImageFilter::SetVariance().
> *
> * \ingroup ImageEnhancement Singlethreaded
> */
> template <typename TInputImage, typename TOutputImage=TInputImage>
> class ITK_EXPORT RecursiveGaussianImageFilter :
> public RecursiveSeparableImageFilter<TInputImage,TOutputImage>
> {
> public:
> /** Standard class typedefs. */
> typedef RecursiveGaussianImageFilter Self;
> typedef RecursiveSeparableImageFilter<TInputImage,TOutputImage> Superclass;
> typedef SmartPointer<Self> Pointer;
> typedef SmartPointer<const Self> ConstPointer;
>
> typedef typename Superclass::RealType RealType;
>
> /** Method for creation through the object factory. */
> itkNewMacro(Self);
>
> /** Type macro that defines a name for this class */
> itkTypeMacro( RecursiveGaussianImageFilter, RecursiveSeparableImageFilter );
>
> /** Set/Get the Sigma, measured in world coordinates, of the Gaussian
> * kernel. The default is 1.0. */
> itkGetMacro( Sigma, RealType );
> itkSetMacro( Sigma, RealType );
>
> /** Enum type that indicates if the filter applies the equivalent operation
> of convolving with a gaussian, first derivative of a gaussian or the
> second derivative of a gaussian. */
> typedef enum { ZeroOrder, FirstOrder, SecondOrder } OrderEnumType;
>
> /** Type of the output image */
> typedef TOutputImage OutputImageType;
>
>
> /** Set/Get the flag for normalizing the gaussian over scale Space
> When this flag is ON the filter will be normalized in such a way
> that larger sigmas will not result in the image fading away.
>
> \f[
> \frac{ 1 }{ \sigma \sqrt{ 2 \pi } };
> \f]
>
> When the flag is OFF the normalization will conserve contant the
> integral of the image intensity.
> \f[
> \frac{ 1 }{ \sigma^2 \sqrt{ 2 \pi } };
> \f]
> For analyzing an image across Scale Space you want to enable
> this flag. It is disabled by default. */
> itkSetMacro( NormalizeAcrossScale, bool );
> itkGetMacro( NormalizeAcrossScale, bool );
>
> /** Set/Get the Order of the Gaussian to convolve with.
> \li ZeroOrder is equivalent to convolving with a Gaussian. This
> is the default.
> \li FirstOrder is equivalet to convolving with the first derivative of a Gaussian
> \li SecondOrder is equivalet to convolving with the second derivative of a Gaussian
> */
> itkSetMacro( Order, OrderEnumType );
> itkGetMacro( Order, OrderEnumType );
>
> /** Explicitly set a zeroth order derivative */
> void SetZeroOrder();
>
> /** Explicitly set a first order derivative */
> void SetFirstOrder();
>
> /** Explicitly set a second order derivative */
> void SetSecondOrder();
>
>
> protected:
> RecursiveGaussianImageFilter();
> virtual ~RecursiveGaussianImageFilter() {};
> void PrintSelf(std::ostream& os, Indent indent) const;
>
> /** Set up the coefficients of the filter to approximate a specific kernel.
> * typically it can be used to approximate a gaussian or one of its
> * derivatives. */
> virtual void SetUp(void);
>
> /** Compute Recursive Filter Coefficients this method prepares the values of
> * the coefficients used for filtering the image. The symmetric flag is
> * used to enforce that the filter will be symmetric or antisymmetric. For
> * example, the Gaussian kernel is symmetric, while its first derivative is
> * antisymmetric. */
> void ComputeFilterCoefficients(bool symmetric);
>
> private:
> RecursiveGaussianImageFilter(const Self&); //purposely not implemented
> void operator=(const Self&); //purposely not implemented
>
> /** Sigma of the gaussian kernel. */
> RealType m_Sigma;
>
> /** Normalize the image across scale space */
> bool m_NormalizeAcrossScale;
>
> OrderEnumType m_Order;
>
> };
>
> } // end namespace itk
>
> #ifndef ITK_MANUAL_INSTANTIATION
> #include "itkRecursiveGaussianImageFilter.txx"
> #endif
>
> #endif
>
>
>
> ------------------------------------------------------------------------
>
> /*=========================================================================
>
> Program: Insight Segmentation & Registration Toolkit
> Module: $RCSfile: itkRecursiveGaussianImageFilter.txx,v $
> Language: C++
> Date: $Date: 2003/09/10 14:28:55 $
> Version: $Revision: 1.27 $
>
> Copyright (c) Insight Software Consortium. All rights reserved.
> See ITKCopyright.txt or http://www . itk . org/HTML/Copyright . htm for details.
>
> This software is distributed WITHOUT ANY WARRANTY; without even
> the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
> PURPOSE. See the above copyright notices for more information.
>
> =========================================================================*/
> #ifndef _itkRecursiveGaussianImageFilter_txx
> #define _itkRecursiveGaussianImageFilter_txx
>
> #include "itkRecursiveGaussianImageFilter.h"
> #include "itkObjectFactory.h"
> #include "itkImageLinearIteratorWithIndex.h"
> #include <new>
>
>
> namespace itk
> {
>
> template <typename TInputImage, typename TOutputImage>
> RecursiveGaussianImageFilter<TInputImage,TOutputImage>
> ::RecursiveGaussianImageFilter()
> {
> m_Sigma = 1.0;
> m_NormalizeAcrossScale = false;
> m_Order = ZeroOrder;
> }
>
> /**
> * Explicitly set a zeroth order derivative
> */
> template <typename TInputImage, typename TOutputImage>
> void
> RecursiveGaussianImageFilter<TInputImage,TOutputImage>
> ::SetZeroOrder()
> {
> m_Order = ZeroOrder;
> }
>
> /**
> * Explicitly set a first order derivative
> */
> template <typename TInputImage, typename TOutputImage>
> void
> RecursiveGaussianImageFilter<TInputImage,TOutputImage>
> ::SetFirstOrder()
> {
> m_Order = FirstOrder;
> }
>
> /**
> * Explicitly set a second order derivative
> */
> template <typename TInputImage, typename TOutputImage>
> void
> RecursiveGaussianImageFilter<TInputImage,TOutputImage>
> ::SetSecondOrder()
> {
> m_Order = SecondOrder;
> }
>
>
> /**
> * Compute filter for Gaussian kernel
> */
> template <typename TInputImage, typename TOutputImage>
> void
> RecursiveGaussianImageFilter<TInputImage,TOutputImage>
> ::SetUp(void)
> {
>
> const RealType spacingTolerance = 1e-8;
>
> RealType direction = 1.0;
> if( m_Spacing < 0.0 )
> {
> direction = -1.0;
> m_Spacing = -m_Spacing;
> }
>
> if( m_Spacing < spacingTolerance )
> {
> itkExceptionMacro(<<"The spacing " << m_Spacing << "is suspiciosly small in this image");
> }
>
> const RealType sigmad = m_Sigma/m_Spacing;
>
> if( this->GetNormalizeAcrossScale() )
> {
> m_K = 1.0 / ( sigmad * sqrt( 2.0 * ( 4.0 * atan( 1.0f ) ) ) );
> }
> else
> {
> m_K = 1.0 / ( sigmad * sigmad * sqrt( 2.0 * ( 4.0 * atan( 1.0f ) ) ) );
> }
>
> m_K *= direction; // take into account the sign of the spacing.
>
> switch( m_Order )
> {
> case ZeroOrder: // equivalent to convolution with a gaussian
> {
> m_A0 = static_cast<RealType>( 1.680 );
> m_A1 = static_cast<RealType>( 3.735 );
> m_B0 = static_cast<RealType>( 1.783 );
> m_B1 = static_cast<RealType>( 1.723 );
> m_C0 = static_cast<RealType>( -0.6803 );
> m_C1 = static_cast<RealType>( -0.2598 );
> m_W0 = static_cast<RealType>( 0.6318 );
> m_W1 = static_cast<RealType>( 1.9970 );
> const bool symmetric = true;
> this->ComputeFilterCoefficients(symmetric);
> break;
> }
> case FirstOrder: // equivalent to convolution with
> // the first derivative of a gaussian
> {
> m_A0 = static_cast<RealType>( -0.6472 );
> m_A1 = static_cast<RealType>( -4.5310 );
> m_B0 = static_cast<RealType>( 1.5270 );
> m_B1 = static_cast<RealType>( 1.5160 );
> m_C0 = static_cast<RealType>( 0.6494 );
> m_C1 = static_cast<RealType>( 0.9557 );
> m_W0 = static_cast<RealType>( 0.6719 );
> m_W1 = static_cast<RealType>( 2.0720 );
> const bool symmetric = false;
> this->ComputeFilterCoefficients(symmetric);
> break;
> }
> case SecondOrder: // equivalent to convolution with
> // the second derivative of a gaussian
> {
> m_A0 = static_cast<RealType>( -1.3310 );
> m_A1 = static_cast<RealType>( 3.6610 );
> m_B0 = static_cast<RealType>( 1.2400 );
> m_B1 = static_cast<RealType>( 1.3140 );
> m_C0 = static_cast<RealType>( 0.3225 );
> m_C1 = static_cast<RealType>( -1.7380 );
> m_W0 = static_cast<RealType>( 0.7480 );
> m_W1 = static_cast<RealType>( 2.1660 );
> const bool symmetric = true;
> this->ComputeFilterCoefficients(symmetric);
> break;
> }
> default:
> {
> itkExceptionMacro(<<"Unknown Order");
> return;
> }
> }
>
>
> }
>
>
>
> /**
> * Compute Recursive Filter Coefficients
> */
> template <typename TInputImage, typename TOutputImage>
> void
> RecursiveGaussianImageFilter<TInputImage,TOutputImage>
> ::ComputeFilterCoefficients(bool symmetric)
> {
>
> const RealType spacing = (m_Spacing>0.0) ? m_Spacing : -m_Spacing;
>
> const RealType sigmad = m_Sigma / spacing;
>
> m_N00 = m_A0 + m_C0;
> m_N11 = exp(-m_B1/sigmad)*(m_C1*sin(m_W1/sigmad)-(m_C0+2*m_A0)*cos(m_W1/sigmad));
> m_N11 += exp(-m_B0/sigmad)*(m_A1*sin(m_W0/sigmad)-(m_A0+2*m_C0)*cos(m_W0/sigmad));
> m_N22 = ((m_A0+m_C0)*cos(m_W1/sigmad)*cos(m_W0/sigmad));
> m_N22 -= (m_A1*cos(m_W1/sigmad)*sin(m_W0/sigmad)+m_C1*cos(m_W0/sigmad)*sin(m_W1/sigmad));
> m_N22 *= 2*exp(-(m_B0+m_B1)/sigmad);
> m_N22 += m_C0*exp(-2*m_B0/sigmad) + m_A0*exp(-2*m_B1/sigmad);
> m_N33 = exp(-(m_B1+2*m_B0)/sigmad)*(m_C1*sin(m_W1/sigmad)-m_C0*cos(m_W1/sigmad));
> m_N33 += exp(-(m_B0+2*m_B1)/sigmad)*(m_A1*sin(m_W0/sigmad)-m_A0*cos(m_W0/sigmad));
>
> m_D44 = exp(-2*(m_B0+m_B1)/sigmad);
> m_D33 = -2*cos(m_W0/sigmad)*exp(-(m_B0+2*m_B1)/sigmad);
> m_D33 += -2*cos(m_W1/sigmad)*exp(-(m_B1+2*m_B0)/sigmad);
> m_D22 = 4*cos(m_W1/sigmad)*cos(m_W0/sigmad)*exp(-(m_B0+m_B1)/sigmad);
> m_D22 += exp(-2*m_B1/sigmad)+exp(-2*m_B0/sigmad);
> m_D11 = -2*exp(-m_B1/sigmad)*cos(m_W1/sigmad)-2*exp(-m_B0/sigmad)*cos(m_W0/sigmad);
>
> if( symmetric )
> {
> m_M11 = m_N11 - m_D11 * m_N00;
> m_M22 = m_N22 - m_D22 * m_N00;
> m_M33 = m_N33 - m_D33 * m_N00;
> m_M44 = - m_D44 * m_N00;
> }
> else
> {
> m_M11 = -( m_N11 - m_D11 * m_N00 );
> m_M22 = -( m_N22 - m_D22 * m_N00 );
> m_M33 = -( m_N33 - m_D33 * m_N00 );
> m_M44 = m_D44 * m_N00;
> }
>
> // Compute Coefficients to be used at the boundaries
> // in order to prevent border effects
> const RealType SumOfNCoefficients = m_N00 + m_N11 + m_N22 + m_N33;
> const RealType SumOfMCoefficients = m_M11 + m_M22 + m_M33 + m_M44;
> const RealType SumOfDCoefficients = m_D11 + m_D22 + m_D33 + m_D44;
> const RealType CoefficientNormN = SumOfNCoefficients / ( 1.0 + SumOfDCoefficients );
> const RealType CoefficientNormM = SumOfMCoefficients / ( 1.0 + SumOfDCoefficients );
>
> m_BN1 = m_D11 * CoefficientNormN;
> m_BN2 = m_D22 * CoefficientNormN;
> m_BN3 = m_D33 * CoefficientNormN;
> m_BN4 = m_D44 * CoefficientNormN;
>
> m_BM1 = m_D11 * CoefficientNormM;
> m_BM2 = m_D22 * CoefficientNormM;
> m_BM3 = m_D33 * CoefficientNormM;
> m_BM4 = m_D44 * CoefficientNormM;
>
> }
>
>
> template <typename TInputImage, typename TOutputImage>
> void
> RecursiveGaussianImageFilter<TInputImage,TOutputImage>
> ::PrintSelf(std::ostream& os, Indent indent) const
> {
> Superclass::PrintSelf(os,indent);
>
> os << "Sigma: " << m_Sigma << std::endl;
> os << "Order: " << m_Order << std::endl;
> os << "NormalizeAcrossScale: " << m_NormalizeAcrossScale << std::endl;
> }
>
> } // end namespace itk
>
> #endif
>