[Insight-developers] Some Speedup in itkConstNeighboorhoodIterator
and AnisotropicDiffusion
Luis Ibanez
luis.ibanez at kitware.com
Wed Jul 13 13:08:42 EDT 2005
Hi Karl,
The Experimental build went green,
so your changes were committed this morning
with the appropriate credits in the CVS logs.
Thanks a lot !
Please let us know if you find other improvements to make.
Regards,
Luis
--------------------
Karl Krissian wrote:
>
> 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
>
>
> ------------------------------------------------------------------------
>
> /*=========================================================================
>
> 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
>
>
> ------------------------------------------------------------------------
>
> /*=========================================================================
>
> 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
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> Insight-developers mailing list
> Insight-developers at itk.org
> http://www.itk.org/mailman/listinfo/insight-developers
More information about the Insight-developers
mailing list