Difference between ITK and elastix BSplineTransform dealing with input points at the border

[N-Dekker]

It appears that the elastix version of BSplineTransform excludes points at the border of the TransformDomainPhysicalDimensions, while ITK does include them, with the transformation.

For example when the OptimizerParameters are {0.75, 0.75, ... , 0.75}, and the input point is like {1.0, y} or {x, 1.0}, then elastix TransformPoint(inputPoint) just returns the input point, while ITK TransformPoint(inputPoint) yields {x + 0.75, y + 0.75}.

ITK TransformPoint(inputPoint) even transforms an input point to {x + 0.75, y + 0.75} when the point is an epsilon beyond the border.

[N-Dekker]

@mstaring @stefanklein These look like minor rounding errors to me, because the BSplineTransform implementations of ITK and elastix are slightly different. However, in order to smoothly exchange ITK Transform files with elastix, would it be worth to get elastix BSpline "in sync" with ITK?

[mstaring]

If the difference is only at the border, could it be due to the decision when we are at the border and when not?

This could be an region->IsInside( x ) difference for example.

Let's discuss at the sprint.

[N-Dekker]

It appears to be caused by a difference between ITK BSplineTransform::InsideValidRegion(ContinuousIndexType &index ) const and elastix AdvancedBSplineDeformableTransformBase::InsideValidRegion(const ContinuousIndexType & index) const.

Elastix: https://github.com/SuperElastix/elastix/blob/5.0.1/Common/Transforms/itkAdvancedBSplineDeformableTransformBase.hxx#L588-L605

ITK: https://github.com/InsightSoftwareConsortium/ITK/blob/v5.2.0/Modules/Core/Transform/include/itkBSplineTransform.hxx#L479-L509

Which does:

index[j] = Math::FloatAddULP(maxLimit, -6);

ITK started modifying the parameter, index, apparently, with this commit: https://github.com/InsightSoftwareConsortium/ITK/commit/054e6b70430bb478a33f31241cf186fb5dd1a93c Nick Tustison,10 June 2011, "BUG: Fixing transform domain initializer and transform domain definition."

Specifically on this line of code:

  index[j] -= 1e-6;

https://github.com/InsightSoftwareConsortium/ITK/blob/054e6b70430bb478a33f31241cf186fb5dd1a93c/Modules/Core/Transform/include/itkBSplineDeformableTransform.hxx#L685

[mstaring]

Then I propose to follow the ITK. Can we do it in such a way that we automatically pick up new modifications as well, e.g. by using inheritance?

[N-Dekker]

Then I propose to follow the ITK.

Honestly I don't really like the ITK approach. In my opinion, InsideValidRegion(index) should not modify the index. The adjustment of setting index[j] = Math::FloatAddULP(maxLimit, -6), as ITK now does, looks like some kind of workaround to me.

Can we do it in such a way that we automatically pick up new modifications as well, e.g. by using inheritance?

I think that would be a major redesign. But it may be worth it 🤔 having to keep our elastix transforms manually "in sync" with the corresponding ITK transforms (by copy-paste) does not seem workable in the long term.

[mstaring]

ok, let's discuss during our sprint. The major re-design may be put in an issue to work on later.

[N-Dekker]

I just made a few pull request to significantly improve the performance of ITK's BSpline implementation, especially https://github.com/InsightSoftwareConsortium/ITK/pull/2716 (BSplineKernelFunction::FastEvaluate(u): a faster way to evaluate BSplineKernelFunction) and https://github.com/InsightSoftwareConsortium/ITK/pull/2712 (Use FixedArray for BSplineInterpolationWeightFunction OutputType)

It's really a pity that those improvements do not get into elastix automatically.

[N-Dekker]

First attempt to implement elastix AdvancedTranslationTransform in terms of ITK TranslationTransform: https://github.com/SuperElastix/elastix/tree/AdvancedTranslationTransform-ITK-TranslationTransform Certainly a code reduction 😃

[N-Dekker]

The two InsideValidRegion(const ContinuousIndexType & index) member functions in Common/Transforms look identical, but not entirely!

https://github.com/SuperElastix/elastix/blob/5.0.1/Common/Transforms/itkAdvancedBSplineDeformableTransformBase.hxx#L589-L605

does for( unsigned int j = 0; j < SpaceDimension; j++ )

Whereas https://github.com/SuperElastix/elastix/blob/5.0.1/Common/Transforms/itkCyclicBSplineDeformableTransform.hxx#L68-L84

does for( unsigned int j = 0; j < SpaceDimension - 1; j++ )

Would that difference be intended, or is it a bug in CyclicBSplineDeformableTransform?

[mstaring]

that is intended due to the nature of the cyclic transform