Custom border extrapolation of ShapedImageNeighborhoodRange by ImageNeighborhoodPixelAccessPolicy


(Niels Dekker) #1

Hereby I would like explain the idea behind the patch I submitted last Saturday, ENH: Added ImageNeighborhoodPixelAccessPolicy for custom border extrapolation

I guess you have seen, last week, the ShapedImageNeighborhoodRange class for iteration over a neighborhood of pixels was merged onto the master branch of ITK. ShapedImageNeighborhoodRange has one specific way to deals with neighborhood pixel locations outside the image boundaries: it replicates the image border, just like itk::ZeroFluxNeumannBoundaryCondition. This border extrapolation method was implemented by its operator*() overload.

The new patch allows ITK users and developers to specify a different border extrapolation approach. It adds an extra template parameter, TImageNeighborhoodPixelAccessPolicy, to ShapedImageNeighborhoodRange, which specifies the “pixel access policy”: it specifies both the access to the pixel data and the border extrapolation method.

By default, BorderReplicatingImageNeighborhoodPixelAccessPolicy is used. But anyone can provide a custom policy class, as long as it has a GetPixelValue, a SetPixelValue, and a constructor, as follows:

template <typename TImage> class MyPolicy
{
  using NeighborhoodAccessorFunctorType = typename TImage::NeighborhoodAccessorFunctorType;
  using PixelType = typename TImage::PixelType;
  using InternalPixelType = typename TImage::InternalPixelType;
  using ImageDimensionType = typename TImage::ImageDimensionType;
  static constexpr ImageDimensionType ImageDimension = TImage::ImageDimension;

public:
  MyPolicy(
    const Size<ImageDimension> & imageSize,
    const OffsetType& offsetTable,
    const NeighborhoodAccessorFunctorType& neighborhoodAccessor,
    const IndexType& pixelIndex) ITK_NOEXCEPT

  PixelType GetPixelValue(
    const InternalPixelType* const imageBufferPointer) const ITK_NOEXCEPT;

  void SetPixelValue(
    InternalPixelType* const imageBufferPointer,
    const PixelType& pixelValue) const ITK_NOEXCEPT;
};

Please have a look at how BorderReplicatingImageNeighborhoodPixelAccessPolicy was implemented: http://review.source.kitware.com/#/c/23366/4/Modules/Core/Common/include/itkBorderReplicatingImageNeighborhoodPixelAccessPolicy.h

In a similar way, other boundary conditions could also be implemented. With zero-cost runtime overhead, because the selection of the policy class is done at compile-time :slight_smile:

Of course, I think it would be nice if this patch could still be included with ITK 5.0.0. But I can understand it it’s now too late in the game. It also depends on whether users would want to have this feature quickly. What do you think?

[Update, May 4, 2018: adapted text to Patch Set 4.]


(Niels Dekker) #2

Regarding the reply from @matt.mccormick at at http://review.source.kitware.com/#/c/23366/3

I chose the first part of the class name BorderReplicatingImageNeighborhoodPixelAccessPolicy, based on the term “border replication”, that is used by OpenCV (which defines an enum BORDER_REPLICATE, described at their tutorial) and MATLAB (whose imfilter function has a border replication option). On the other hand, I never really understood the term “zero-flux Neumann boundary condition”. Do you still think this border extrapolation method is best described by “zero-flux Neumann”?

The other part of the name that I chose, PixelAccessPolicy was based on the fact that it is a policy class, very much like the policy classes Andrei Alexandrescu presented in his book Modern C++ Design, https://en.wikipedia.org/wiki/Policy-based_design And it specifies the pixel access, by providing a GetPixelValue and a SetPixelValue member function. (Internally BorderReplicatingImageNeighborhoodPixelAccessPolicy uses the NeighborhoodAccessor of the image, but that’s an implementation detail!)

If you still don’t like the name BorderReplicatingImageNeighborhoodPixelAccessPolicy after this explanation, please let me know if you see other, more technical flaws in the design as well! Because in that case, things might need to be changed anyway… :thinking:


(Matt McCormick) #3

There will be other alpha’s and release candidates for ITK 5, and other opportunities past ITK 5.0.0 – we should not rush it in.

Yes, this is zero-flux Neumann behavior. We need to consistently use ITK naming conventions and semantics within ITK to avoid confusion.

Policy could be included in the name or in the comments. A lesson learned from itk::ImageBoundaryCondition is that object-oriented programming makes policy design potentially problematic after we start creating something complex. In particular, it can be easy to create incompatibilities because classes or their parents or their members were templated with different policies – but then an attempt is made to use these complex object systems together, and it causes issues. But, this could be mitigated by:

  1. Do not provide a default template parameter for the proxy. This will encourage being explicit with the policy in a consistent way, and it should reduce the number of times when the policy is the default and when the policy was explicitly set to something other than the default.
  2. State in the documentation that this class should only be used as a template parameter and not as an instance. Enforce this if possible, with whatever mechanisms C++ provides (remove default constructors, etc.) This will help avoid issues when a templated class is passed a pointer to a type with a different policy.

(Niels Dekker) #4

Thanks for your reply, Matt!

Personally I think the patch ENH: Added ImageNeighborhoodPixelAccessPolicy for custom border extrapolation is almost ready now, but it’s OK to me to postpone the feature until somebody actually needs it!

Honestly I always found both terms “boundary condition” and “zero-flux Neumann boundary condition” quite confusing. They seem to be from physics, rather than image processing. That’s why I searched for different terms, and found “border extrapolation” and “border replication”, which seem much clearer to me. But I guess you feel differently… no problem!

Anyway, even if the patch doesn’t make it into the next version of ITK, it can still serve as a “proof of principle”, showing that this feature can be implemented any time later, with zero run-time overhead :slight_smile:


(Matt McCormick) #5

I like “border”, too, but for the codebase to be readable and understandable, consistency throughout the toolkit has to take precedent over our personal preferences.


(Niels Dekker) #6

Update: Here’s an alternative patch, hoping to bypass our unresolved naming discussion (“border replication” versus “zero-flux Neumann”), by just leaving the default behavior unnamed:

ENH: Allowed custom PixelAccessPolicy for ShapedImageNeighborhoodRange

Although I still think the name BorderReplicatingImageNeighborhoodPixelAccessPolicy is fine (as I proposed at http://review.source.kitware.com/#/c/23366), I find it more important to allow user customization of this policy, one way or another.


(Pablo Hernandez-Cerdan) #7

Hi Niels! Awesome work on the iterators, I am looking forward to using them in the wild!

It doesn’t seem unresolved. Your pull-request has received reasoned feedback about why zero flux is preferred.

Why unnamed? Why is ZeroFlux the unnamed default? If any of the boundary conditions make sense (to me) as default would be the constant boundary condition with a background pixel value, but this has the problem of specifying the background value.
So I agree with @matt.mccormick, relying on default parameters has caused and will cause problems in future filters when they use these iterators but do not provide an interface to change the BoundaryCondition, just because there is a default the creators haven’t thought about it.

Maybe, to help with the usage, and also with the name, could we create type alias?

template<typename TImageType>
using ZeroFluxNeumannNeighborhoodRange<TImageType> =  ShapedImageNeighborhoodRange<TImageType, ZeroFluxNeumannBoundaryCondition>

template<typename TImageType>
using BorderReplicatingNeighborhoodRange<TImageType> =  ShapedImageNeighborhoodRange<TImageType, ZeroFluxNeumannBoundaryCondition>

template<typename TImageType>
using ConstantBoundaryNeighborhoodRange<TImageType> =  ShapedImageNeighborhoodRange<TImageType, ConstantBoundaryCondition>

However, what if the boundary condition or boundary condition policy use run-time arguments? Like in the case of the itk::ConstantBoundaryCondition?


(Niels Dekker) #8

Thanks for your reply, Pablo!

Sorry, but didn’t I give a reasonable explanation on why I prefer the term “border replication”? Of course, anyone of you can rename the policy class to something you prefer (if possible, in a separate commit).

Good question. I thought about this before, allowing to user to add run-time arguments tweak the pixel access. But I felt that it could make the design of the range class too complicated. Also I wasn’t sure if it would be an important use case. Is ConstantBoundaryCondition really used like that, in practice? I mean, isn’t the constant value used for ConstantBoundaryCondition usually known at compile-time, in practice?

By the way, I used the term policy class to indicate that it is a template argument that tweaks the pixel access entirely at compile-time. I think that’s according to the definition of a policy class in C++: https://en.wikipedia.org/wiki/Policy-based_design If the class should also support run-time arguments, the term policy class might not be entirely appropriate anymore :thinking:


(Pablo Hernandez-Cerdan) #9

In my view, not really Niels, you prefer one name over the other, I can understand that and even agree that the name is easier to grasp. But you don’t seem to put any weight or importance in the reasons given by the reviewers: consistency over the toolkit. There are currently these classes, named after that boundary condition:

./Modules/Filtering/ImageGrid/include/itkZeroFluxNeumannPadImageFilter.h
./Modules/Core/Common/include/itkZeroFluxNeumannBoundaryCondition.h

Should we deprecate them, and rename them to BorderReplicating? Or should both terms live at the same time, adding extra complexity? I don’t think it is worth the pain either, but we should
add BorderReplicating to the docs for people coming from others frameworks for sure.

These are the classes using ZeroFlux: ack --cpp -m1 "ZeroFlux" | grep -v test

Modules/Segmentation/LabelVoting/include/itkVotingBinaryHoleFillingImageFilter.hxx
Modules/Segmentation/LabelVoting/include/itkVotingBinaryImageFilter.hxx
Modules/Segmentation/LabelVoting/include/itkBinaryMedianImageFilter.hxx
Modules/Segmentation/LevelSetsv4/include/itkUpdateWhitakerSparseLevelSet.hxx
Modules/Segmentation/LevelSetsv4/include/itkLevelSetEquationCurvatureTerm.h
Modules/Segmentation/LevelSetsv4/include/itkUpdateShiSparseLevelSet.hxx
Modules/Segmentation/LevelSetsv4/include/itkUpdateMalcolmSparseLevelSet.hxx
Modules/Segmentation/LevelSetsv4/include/itkLevelSetEquationAdvectionTerm.h
Modules/Segmentation/LevelSetsv4/include/itkLevelSetEquationPropagationTerm.h
Modules/Segmentation/LevelSetsv4/include/itkLevelSetEquationLaplacianTerm.h
Modules/Filtering/ImageGrid/include/itkZeroFluxNeumannPadImageFilter.h
Modules/Segmentation/LevelSetsv4/include/itkBinaryImageToLevelSetImageAdaptor.hxx
Modules/Filtering/ImageGrid/include/itkZeroFluxNeumannPadImageFilter.hxx
Modules/Filtering/MathematicalMorphology/include/itkGrayscaleGeodesicDilateImageFilter.hxx
Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.hxx
Modules/Filtering/Denoising/include/itkPatchBasedDenoisingBaseImageFilter.h
Modules/Filtering/GPUImageFilterBase/include/itkGPUNeighborhoodOperatorImageFilter.h
Modules/Filtering/MathematicalMorphology/include/itkGrayscaleGeodesicErodeImageFilter.hxx
Modules/Filtering/Smoothing/include/itkBoxUtilities.h
Modules/Filtering/Smoothing/include/itkMeanImageFilter.hxx
Modules/Filtering/CurvatureFlow/include/itkCurvatureFlowFunction.h
Modules/Filtering/Smoothing/include/itkMedianImageFilter.hxx
Modules/Filtering/CurvatureFlow/include/itkBinaryMinMaxCurvatureFlowFunction.h
Modules/Filtering/AnisotropicSmoothing/include/itkVectorAnisotropicDiffusionFunction.hxx
Modules/Filtering/AnisotropicSmoothing/include/itkScalarAnisotropicDiffusionFunction.hxx
Modules/Filtering/Convolution/include/itkFFTConvolutionImageFilter.h
Modules/Filtering/Convolution/include/itkConvolutionImageFilter.h
Modules/Filtering/Convolution/include/itkConvolutionImageFilterBase.h
Modules/Filtering/FFT/include/itkFFTPadImageFilter.h
Modules/Filtering/DistanceMap/include/itkContourDirectedMeanDistanceImageFilter.hxx
Modules/Core/FiniteDifference/include/itkFiniteDifferenceFunction.h
Modules/Filtering/ImageFilterBase/include/itkNoiseImageFilter.hxx
Modules/Filtering/ImageGradient/include/itkGradientImageFilter.hxx
Modules/Filtering/ImageFilterBase/include/itkNeighborhoodOperatorImageFilter.h
Modules/Filtering/ImageGradient/include/itkVectorGradientMagnitudeImageFilter.hxx
Modules/Filtering/ImageGradient/include/itkGradientMagnitudeImageFilter.hxx
Modules/Filtering/CurvatureFlow/include/itkMinMaxCurvatureFlowFunction.h
Modules/Filtering/ImageFeature/include/itkZeroCrossingImageFilter.hxx
Modules/Core/TestKernel/include/itkTestingComparisonImageFilter.hxx
Modules/Filtering/ImageFeature/include/itkCannyEdgeDetectionImageFilter.hxx
Modules/Filtering/ImageFeature/include/itkSobelEdgeDetectionImageFilter.hxx
Modules/Filtering/ImageFeature/include/itkBilateralImageFilter.hxx
Modules/Filtering/ImageFeature/include/itkLaplacianSharpeningImageFilter.hxx
Modules/Filtering/ImageFeature/include/itkSimpleContourExtractorImageFilter.hxx
Modules/Filtering/ImageFeature/include/itkLaplacianImageFilter.hxx
Modules/Core/ImageFunction/include/itkMeanImageFunction.h
Modules/Core/ImageFunction/include/itkVarianceImageFunction.h
Modules/Core/ImageFunction/include/itkCovarianceImageFunction.h
Modules/Core/ImageFunction/include/itkVectorMeanImageFunction.h
Modules/Core/ImageFunction/include/itkSumOfSquaresImageFunction.h
Modules/Core/ImageFunction/include/itkScatterMatrixImageFunction.h
Modules/Core/ImageFunction/include/itkMedianImageFunction.h
Modules/Core/ImageFunction/include/itkWindowedSincInterpolateImageFunction.h
Modules/Core/Common/include/itkZeroFluxNeumannBoundaryCondition.h
Modules/Core/Common/include/itkConstNeighborhoodIterator.h
Modules/Core/Common/include/itkShapedNeighborhoodIterator.h
Modules/Core/Common/include/itkConstShapedNeighborhoodIterator.hxx
Modules/Core/Common/include/itkZeroFluxNeumannBoundaryCondition.hxx
Modules/Core/Common/include/itkNeighborhoodIterator.h
Modules/Core/Common/include/itkConstShapedNeighborhoodIterator.h
Modules/Nonunit/Review/include/itkWarpHarmonicEnergyCalculator.hxx
Modules/Filtering/ImageFeature/include/itkCannyEdgeDetectionImageFilter.h
Modules/Filtering/ImageFeature/include/itkDerivativeImageFilter.hxx
Modules/Registration/Common/include/itkGradientDifferenceImageToImageMetric.h

It would be cool to have the best of both approaches, consistency and a clear name. Maybe it is worth to explore a type alias, but again, not sure if it is worth, the user can just make a 1 to 1 map: BorderReplicating <-> ZeroFluxNeumann.


(Niels Dekker) #10

So what name would you propose for the class that I named BorderReplicatingImageNeighborhoodPixelAccessPolicy at http://review.source.kitware.com/#/c/23366 ? ZeroFluxNeumannBoundaryConditionPixelAccessPolicy?

I’m interested to hear if it’s important to you to allow run-time arguments to tweak the pixel access, as the constant value of itk::ConstantBoundaryCondition. My feeling is that the constant used with itk::ConstantBoundaryCondition is usually zero, or at least a compile-time constant. Did you see other use cases?


(Pablo Hernandez-Cerdan) #11

That sounds good to me.

You are right, maybe it is not important to generalize to accept run-time argument, for the case of itk::ConstantBoundaryCondition, I think providing compile-time values: Zero, One, Min, and Max (from itk:NumericalLimits<TImage::PixelType>) should cover 99% of the use cases of the ConstantBoundaryCondition.

So, Max/Min/Zero/One/CostantBoundaryConditionPixelAccessPolicy or similar?


(Niels Dekker) #12

Sure! The constant value could also be passed as an extra template argument to the policy class. That’s probably fine for most use cases, but it has a limitation, of course: then the constant value must be an integer. (C++ does not support float or double as template argument.)

Another possibility could be to allow the user to specify a functor as an extra template argument to the policy class. This functor could then produce the constant value for the policy, and would be called by PixelAccessPolicy::GetPixelValue whenever the index is outside the image borders. Something like this:

template <typename TImage, typename TFunctor>
class ConstantBoundaryConditionPixelAccessPolicy
{
public:
  PixelType GetPixelValue(const InternalPixelType* const imageBufferPointer) const ITK_NOEXCEPT
  {
    return
      m_IsOutsideImageBorders ?
      TFunctor{}() :
      m_NeighborhoodAccessor.Get(imageBufferPointer + m_PixelIndexValue);
  }
  ...
};

With a user-supplied functor like this:

struct MyFunctor
{
    double operator()() const ITK_NOEXCEPT
    {
        return 0.01;
    }
};

Then again, it would be handy to have some pre-defined Max/Min/Zero/One policies, as you’re suggesting.


(Pablo Hernandez-Cerdan) #13

Nice! The functor would solve the case when the user wants a specific value. Following the example you posted, I was originally thinking in a simpler case with no extra template parameter.

template <typename TImage>
class MaxConstantBoundaryConditionPixelAccessPolicy
{
public:
  PixelType GetPixelValue(const InternalPixelType* const imageBufferPointer) const ITK_NOEXCEPT
  {
    return
      m_IsOutsideImageBorders ?
      NumericTraits<typename TImage::PixelType>::max() :
      m_NeighborhoodAccessor.Get(imageBufferPointer + m_PixelIndexValue);
  }
  ...
};
template <typename TImage>
class ZeroConstantBoundaryConditionPixelAccessPolicy
{
public:
  PixelType GetPixelValue(const InternalPixelType* const imageBufferPointer) const ITK_NOEXCEPT
  {
    return
      m_IsOutsideImageBorders ?
      NumericTraits<typename TImage::PixelType>::ZeroValue() :
      m_NeighborhoodAccessor.Get(imageBufferPointer + m_PixelIndexValue);
  }
  ...
};

One question, why a functor, and not a static function, or a member?


(Niels Dekker) #14

First of all, I just hope we can avoid adding an extra run-time parameter (the “constant value”) to the constructor of ShapedImageNeighborhoodRange<TImage, TPixelAccessPolicy>. That should be possible as long as we may assume that the value of the argument is known at compile-time, in practice. But then this “constant value” should be encoded somehow in the type of the PixelAccessPolicy.

When Functor1 and Functor2 are two different functor classes, producing different “constant values”, the following two template instantiations would then specify the corresponding pixel access policies, without having to add extra runtime parameters:

ConstantBoundaryConditionPixelAccessPolicy<ImageType, Functor1>
ConstantBoundaryConditionPixelAccessPolicy<ImageType, Functor2>

Ideally, people should be able to specify such a functor by a lambda expression, but unfortunately in C++, the type of a lambda does have a default-constructor. (The GetPixelValue function I posted here does TFunctor{} so TFunctor needs to be default-constructible.)

How would you implement this feature (the ability to tweak the pixel access by such a “constant value”), using a static function or a member?


(Pablo Hernandez-Cerdan) #15

Thanks for the explanation, I know that the goal is to avoid run-time parameters. it was just a minor question out of curiosity on why do you prefer a functor over something like:

template<TPixelType>
struct MyConstant255Value
{
    static TPixelType Value()
    {
      return static_cast<TPixelType>(255)
    }
}
ConstantBoundaryConditionPixelAccessPolicy<ImageType, MyConstant255Value>

Is there any pros or cons or one approach over the other, or they are just syntactically different, but with the same result?


(Niels Dekker) #16

I think it’s a good question, actually :thinking: I think your approach with the static Value() member function should also work fine. I just thought that the TFunctor approach would be more generic, more generally usable. Ideally, I would like it to accept any type of function or function object that can be called with zero arguments, f(), and returns a value that can be converted to a pixel value. My initial idea was that it should also accept lambda’s, like this:

ConstantBoundaryConditionPixelAccessPolicy<ImageType, decltype([]{ return 255; })>

But then I realized that lambda’s aren’t DefaultConstructible, at least not according to the current C++ standard. While TFunctor has to be DefaultConstructible! So that won’t work for now.

To be continued…


(Niels Dekker) #17

So far we discussed two different ways to support specifying the constant value for a policy class like ConstantBoundaryConditionPixelAccessPolicy by means of a type template argument:

TFunctor: Allowing the policy to retrieve the constant value by TFunctor{}().

TStaticValueClass: Allowing the policy to retrieve the constant value by TStaticValueClass::Value().

Actually there’s a third possibility which I find quite reasonable, and quite appealing:

TValueConvertible: Allowing the policy to retrieve the constant value, just by TValueConvertible{}. (Without the extra parentheses needed for TFunctor{}()). TValueConvertible{} should then be convertible to the pixel value.

For a simple zero-constant policy and a simple built-in pixel type, TValueConvertible could just be int, as the expression int{} denotes the integer value zero. But in general, such a convertible type could be written, for instance, as follows:

struct My255ValueConvertible
{
  template <typename TPixelType> 
  operator TPixelType() const
  {
    return static_cast<TPixelType>(255);
  }
};

What do you think?


(Pablo Hernandez-Cerdan) #18

I don’t see much difference, Niels! These structs would probably only be used as template parameters of the Policies, so… whatever interface sounds good, the third option with the template in operator() is slightly strange for me, how are you calling it?

My255ValueConvertible my255;
int myValue = my255.operator()<int>

I am more familiar with:

int myValue = My255ValueStatic<int>::Value();

or, if you want to put the template in the operator:


struct MyConstant255Value
{
    template<TPixelType>
    static TPixelType Value()
    {
      return static_cast<TPixelType>(255)
    }
}
int myValue = My255ValueStatic::Value<int>();

But as I said, I think the all the options look good for providing the compile time value.


(Niels Dekker) #19

Thanks for your reply, Pablo.

No, simply as follows:

My255ValueConvertible my255;
int myValue = my255;

You may try it out at http://rextester.com/CWQJ19410 :smiley:


(Pablo Hernandez-Cerdan) #20

That’s super nice! :grin: