[Insight-users] CSwig Wrapping: itkRecursiveGaussianImageFilter SetOrder
Hua Qian
hqian at imaging . robarts . ca
Tue, 14 Oct 2003 16:16:14 -0400
This is a multi-part message in MIME format.
--------------010303020207090502070506
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit
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
>>
>>
>
>
>
--------------010303020207090502070506
Content-Type: text/plain;
name="itkRecursiveGaussianImageFilter.h"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="itkRecursiveGaussianImageFilter.h"
/*=========================================================================
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
--------------010303020207090502070506
Content-Type: text/plain;
name="itkRecursiveGaussianImageFilter.txx"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="itkRecursiveGaussianImageFilter.txx"
/*=========================================================================
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
--------------010303020207090502070506--