[Insight-developers] BUG in BinaryBallStructuringElement

[Padfield, Dirk R (GE Global Research)]

Hi All,

I am writing to ask your advice about a bug I found in BinaryBallStructuringElement.

For a while, I have been bothered by the fact that the BinaryBallStructuringElement return balls that are larger than the specified radius.  For example, when given a radius of 1, it returns the structuring element:
1    1    1
1    1    1
1    1    1

But this structuring element has a radius that is more than 1!  If it truly had a radius of 1, it would be a cross shape in this case.

When choosing a larger radius, the problem is more obvious.  Setting radius = 5 (leading to a structuring element size of 11x11) results in:
0    0    0    1    1    1    1    1    0    0    0
0    0    1    1    1    1    1    1    1    0    0
0    1    1    1    1    1    1    1    1    1    0
1    1    1    1    1    1    1    1    1    1    1
1    1    1    1    1    1    1    1    1    1    1
1    1    1    1    1    1    1    1    1    1    1
1    1    1    1    1    1    1    1    1    1    1
1    1    1    1    1    1    1    1    1    1    1
0    1    1    1    1    1    1    1    1    1    0
0    0    1    1    1    1    1    1    1    0    0
0    0    0    1    1    1    1    1    0    0    0

This is clearly not an ellipse/circle with radius 5 because the interior ellipse/circle is touching each image border at five points rather than just one.  As it turns out, the code is actually defining a radius that is "X + 0.5", where X is the radius that is requested!

The problem is in the specification of the ellipse axes on lines 70-76 of itkBinaryBallStructuringElement.hxx:
  // Define and set the axes lengths for the ellipsoid
  typename EllipsoidType::InputType axes;
  for ( i = 0; i < VDimension; i++ )
    {
    axes[i] = this->GetSize(i);
    }
  spatialFunction->SetAxes(axes);

In this case, "this->GetSize()" is equal to radius*2+1.  But, an ellipse/circle with radius X should have axes length 2X, not 2X+1!  In the implementation, the center of the ellipse is properly accounted for by setting it to "this->GetRadius+1", but the size of the ellipse is not correct!

To correct this, we can make a simple change either
    axes[i] = this->GetSize(i) - 1;
or
    axes[i] = this->GetRadius(i) * 2;

The second is probably more intuitive.

With this fix, we get for radius=1:
0    1    0
1    1    1
0    1    0

and for radius=5:
0    0    0    0    0    1    0    0    0    0    0
0    0    1    1    1    1    1    1    1    0    0
0    1    1    1    1    1    1    1    1    1    0
0    1    1    1    1    1    1    1    1    1    0
0    1    1    1    1    1    1    1    1    1    0
1    1    1    1    1    1    1    1    1    1    1
0    1    1    1    1    1    1    1    1    1    0
0    1    1    1    1    1    1    1    1    1    0
0    1    1    1    1    1    1    1    1    1    0
0    0    1    1    1    1    1    1    1    0    0
0    0    0    0    0    1    0    0    0    0    0

This is a true circle with radius 5!

My questions are:
1) Is anyone else bothered by this bug?  I imagine that many users expect the corrected version and don't realize they are getting the incorrect one.
2) Do others agree with this fix?
3) Can we make this fix given the number of filters/applications that will change slightly as a result of this fix?

Many thanks,
Dirk