[Insight-developers] Some Speedup in itkConstNeighboorhoodIterator
and AnisotropicDiffusion
Karl Krissian
karl at bwh.harvard.edu
Tue Jul 12 19:12:56 EDT 2005
Hi,
After some profiling the itkGradientAnisotropicDiffusion on a
256x256x100, I found the following
speed improvements:
1. In Common/itkConstNeighborhoodIterator.h:
adding 3 lines in GetPixel allows speeding up the processing in the
central region that does not
need checking for boundaries (basically putting this part of the code
inline):
virtual PixelType GetPixel(const unsigned i) const
{
if (!m_NeedToUseBoundaryCondition)
return (*(this->operator[](i)));
else {
bool inbounds;
return this->GetPixel(i, inbounds);
}
}
This speeds up approx 25%
2. Avoiding using the NeighboorhoodInnerProduct to compute finite
differences
can speed up another 30%,
I attach the modified files
BasicFilters/itkGradientNDAnisotropicDiffusionFunction.txx
BasicFilters/itkScalarAnisotropicDiffusionFunction.txx
A test using the itkTimeProbe on the 256x256x100 with 5 iterations gives
the following timings:
initial GradientAD 46 sec
modif 1 35 sec
modif 1+2 25 sec
which means almost x2 speed up.
The methods that take most time in the fastest version are:
- 31.7% GradientNDAnisotropicDiffusionFunction ComputeUpdate
- 24.7% ConstNeighborhoodIterator
GetPixel(unsigned)
- 11.3% ScalarAnisotropicDiffusionFunction
CalculateAverageGradientMagnitudeSquared
- 8.3% ConstNeighborhoodIterator
ComputeInternalIndex
- 6.4% ConstNeighborhoodIterator operator++
- 3.5% ConstNeighborhoodIterator
GetPixel(unsigned, bool&)
- 3.5% ConstNeighborhoodIterator Inbounds()
Karl
--
Karl Krissian, PhD
Instructor in Radiology, Harvard Medical School
Laboratory of Mathematics in Imaging, Brigham and Women's Hospital
Tel:617-525-6232, Fax:617-525-6220
-------------- next part --------------
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkGradientNDAnisotropicDiffusionFunction.txx,v $
Language: C++
Date: $Date: 2003/09/10 14:28:48 $
Version: $Revision: 1.10 $
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 __itkGradientNDAnisotropicDiffusionFunction_txx_
#define __itkGradientNDAnisotropicDiffusionFunction_txx_
#include "itkNumericTraits.h"
namespace itk {
template<class TImage>
double GradientNDAnisotropicDiffusionFunction<TImage>
::m_MIN_NORM = 1.0e-10;
template<class TImage>
GradientNDAnisotropicDiffusionFunction<TImage>
::GradientNDAnisotropicDiffusionFunction()
{
unsigned int i, j;
RadiusType r;
for (i = 0; i < ImageDimension; ++i)
{
r[i] = 1;
}
this->SetRadius(r);
// Dummy neighborhood used to set up the slices.
Neighborhood<PixelType, ImageDimension> it;
it.SetRadius(r);
// Slice the neighborhood
m_Center = it.Size() / 2;
for (i = 0; i< ImageDimension; ++i)
{ m_Stride[i] = it.GetStride(i); }
for (i = 0; i< ImageDimension; ++i)
{ x_slice[i] = std::slice( m_Center - m_Stride[i], 3, m_Stride[i]); }
for (i = 0; i< ImageDimension; ++i)
{
for (j = 0; j < ImageDimension; ++j)
{
// For taking derivatives in the i direction that are offset one
// pixel in the j direction.
xa_slice[i][j]
= std::slice((m_Center + m_Stride[j])-m_Stride[i], 3, m_Stride[i]);
xd_slice[i][j]
= std::slice((m_Center - m_Stride[j])-m_Stride[i], 3, m_Stride[i]);
}
}
// Allocate the derivative operator.
dx_op.SetDirection(0); // Not relevant, will be applied in a slice-based
// fashion.
dx_op.SetOrder(1);
dx_op.CreateDirectional();
}
template<class TImage>
typename GradientNDAnisotropicDiffusionFunction<TImage>::PixelType
GradientNDAnisotropicDiffusionFunction<TImage>
::ComputeUpdate(const NeighborhoodType &it, void *,
const FloatOffsetType&)
{
unsigned int i, j;
double accum;
double accum_d;
double Cx;
double Cxd;
// PixelType is scalar in this context
typedef typename NumericTraits<PixelType>::RealType PixelRealType;
PixelRealType delta;
PixelRealType dx_forward;
PixelRealType dx_backward;
PixelRealType dx[ImageDimension];
PixelRealType dx_aug;
PixelRealType dx_dim;
delta = NumericTraits<PixelRealType>::Zero;
// Calculate the centralized derivatives for each dimension.
for (i = 0; i < ImageDimension; i++)
{ dx[i] = (it.GetPixel(m_Center + m_Stride[i])-it.GetPixel(m_Center - m_Stride[i]))/2.0f; }
for (i = 0; i < ImageDimension; i++)
{
// ``Half'' directional derivatives
dx_forward = it.GetPixel(m_Center + m_Stride[i])
- it.GetPixel(m_Center);
dx_backward = it.GetPixel(m_Center)
- it.GetPixel(m_Center - m_Stride[i]);
// Calculate the conductance terms. Conductance varies with each
// dimension because the gradient magnitude approximation is different
// along each dimension.
accum = 0.0;
accum_d = 0.0;
for (j = 0; j < ImageDimension; j++)
{
if (j != i)
{
//dx_aug = m_InnerProduct(xa_slice[j][i], it, dx_op);
//dx_dim = m_InnerProduct(xd_slice[j][i], it, dx_op);
dx_aug = (it.GetPixel(m_Center + m_Stride[i] + m_Stride[j]) -
it.GetPixel(m_Center + m_Stride[i] - m_Stride[j]))/2.0f;
dx_dim = (it.GetPixel(m_Center - m_Stride[i] + m_Stride[j]) -
it.GetPixel(m_Center - m_Stride[i] - m_Stride[j]))/2.0f;
accum += 0.25f * vnl_math_sqr( dx[j] + dx_aug );
accum_d += 0.25f * vnl_math_sqr( dx[j] + dx_dim );
}
}
if (m_K == 0.0)
{
Cx = 0.0;
Cxd = 0.0;
}
else
{
Cx = exp(( vnl_math_sqr( dx_forward ) + accum) / m_K );
Cxd= exp(( vnl_math_sqr( dx_backward) + accum_d)/ m_K );
}
// Conductance modified first order derivatives.
dx_forward = dx_forward * Cx;
dx_backward = dx_backward * Cxd;
// Conductance modified second order derivative.
delta += dx_forward - dx_backward;
}
return static_cast<PixelType>(delta);
}
} // end namespace itk
#endif
-------------- next part --------------
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkScalarAnisotropicDiffusionFunction.txx,v $
Language: C++
Date: $Date: 2004/12/21 22:47:31 $
Version: $Revision: 1.8 $
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 __itkScalarAnisotropicDiffusionFunction_txx_
#define __itkScalarAnisotropicDiffusionFunction_txx_
#include "itkZeroFluxNeumannBoundaryCondition.h"
#include "itkConstNeighborhoodIterator.h"
#include "itkNeighborhoodInnerProduct.h"
#include "itkNeighborhoodAlgorithm.h"
#include "itkDerivativeOperator.h"
//#include <iostream.h>
namespace itk {
template <class TImage>
void
ScalarAnisotropicDiffusionFunction<TImage>
::CalculateAverageGradientMagnitudeSquared(TImage *ip)
{
typedef ConstNeighborhoodIterator<TImage> RNI_type;
typedef ConstNeighborhoodIterator<TImage> SNI_type;
typedef NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<TImage> BFC_type;
unsigned int i;
ZeroFluxNeumannBoundaryCondition<TImage> bc;
PixelType accumulator;
PixelType val;
PixelType counter;
BFC_type bfc;
typename BFC_type::FaceListType faceList;
typename RNI_type::RadiusType radius;
typename BFC_type::FaceListType::iterator fit;
NeighborhoodInnerProduct<TImage> SIP;
NeighborhoodInnerProduct<TImage> IP;
RNI_type iterator_list[ImageDimension];
SNI_type face_iterator_list[ImageDimension];
DerivativeOperator<PixelType,
ImageDimension> operator_list[ImageDimension];
unsigned long Stride[ImageDimension];
unsigned long Center[ImageDimension];
// Set up the derivative operators, one for each dimension
for (i = 0; i < ImageDimension; ++i)
{
operator_list[i].SetOrder(1);
operator_list[i].SetDirection(i);
operator_list[i].CreateDirectional();
radius[i] = operator_list[i].GetRadius()[i];
}
// Get the various region "faces" that are on the data set boundary.
faceList = bfc(ip, ip->GetRequestedRegion(), radius);
fit = faceList.begin();
// Now do the actual processing
accumulator = NumericTraits<PixelType>::Zero;
counter = NumericTraits<PixelType>::Zero;
// First process the non-boundary region
// Instead of maintaining a single N-d neighborhood of pointers,
// we maintain a list of 1-d neighborhoods along each axial direction.
// This is more efficient for higher dimensions.
for (i = 0; i < ImageDimension; ++i)
{
iterator_list[i]=RNI_type(operator_list[i].GetRadius(), ip, *fit);
iterator_list[i].GoToBegin();
Center[i]=iterator_list[i].Size()/2;
Stride[i]=iterator_list[i].GetStride(i);
}
while ( !iterator_list[0].IsAtEnd() )
{
counter += NumericTraits<PixelType>::One;
for (i = 0; i < ImageDimension; ++i)
{
val = static_cast<PixelType> (iterator_list[i].GetPixel(Center[i]+Stride[i]))-
static_cast<PixelType> (iterator_list[i].GetPixel(Center[i]-Stride[i]));
val = val/-2.0f;
val = static_cast<PixelType>(static_cast<double>(val) * this->m_ScaleCoefficients[i]);
accumulator += val * val;
++iterator_list[i];
}
}
// Go on to the next region(s). These are on the boundary faces.
++fit;
while ( fit != faceList.end() )
{
for (i = 0; i < ImageDimension; ++i)
{
face_iterator_list[i]=SNI_type(operator_list[i].GetRadius(), ip,
*fit);
face_iterator_list[i].OverrideBoundaryCondition(&bc);
face_iterator_list[i].GoToBegin();
Center[i]=face_iterator_list[i].Size()/2;
Stride[i]=face_iterator_list[i].GetStride(i);
}
while ( ! face_iterator_list[0].IsAtEnd() )
{
counter += NumericTraits<PixelType>::One;
for (i = 0; i < ImageDimension; ++i)
{
val = static_cast<PixelType> (face_iterator_list[i].GetPixel(Center[i]+Stride[i]))-
static_cast<PixelType> (face_iterator_list[i].GetPixel(Center[i]-Stride[i]));
val = val/-2.0f;
val = static_cast<PixelType>(static_cast<double>(val) * this->m_ScaleCoefficients[i]);
accumulator += val * val;
++face_iterator_list[i];
}
}
++fit;
}
this->SetAverageGradientMagnitudeSquared( (double) (accumulator / counter) );
}
}// end namespace itk
#endif
More information about the Insight-developers
mailing list