SEHexception when using vtkContour with vtkCutter

Development
1

I’ve discovered a bug that is probably 2 bugs in one when slicing a 3D scalar field using vtkCutter using the plane implicit function point=(0,0,0) normal=(0,0,1) with vtkContour that causes it to throw an SEHexception.
I’m using VTK 8.1.1.

Summary of issues:

  1. vtkBoundingBox::ComputeDivisions() is not robust with very small bounding box lengths
  2. Is there a numerical error (rouding, precision, etc) problem with vtkCutter?

I traced the code and found the problem to be an extremely large memory allocation done in
vtkPointLocator::InitPointInsertion()
at line 1191 of vtkPointLocator.cxx
this->HashTable = new vtkIdListPtr[this->NumberOfBuckets];
because NumberofBuckets=4,342,286,870,571

Digging a little deeper, I found that
vtkBoundingBox::ComputeDivisions()
handles 0 bounding box dimension length (X, Y or Z), but not very small values close to 0, such as
4.4408920985006262e-16
in Z. The other dimension lengths are LengthX=40 and LengthY=30.
This causes a division in ComputeDivisions() to produce a very large number (I won’t get into the details), which eventually produces the large value for NumberOfBuckets.

This can be fixed easily by snapping values of the length to 0 after calling this->GetLengths(lengths) (line 436) and before the look checking for 0 length (line 438) like this:
const relative_tolerance= this->GetMaxLength()*1e-10; (is there an absolute VTK tolerance defined?)
for (int i=0; i<3; ++i)
{
if (lengths[i] <= relative_tolerance)
lengths[i]= 0;
}

This would certainly make the code in ComputeDivisions() more robust.
Has anyone seen this before and has a similar fix already been submitted?

But the second (and stranger) bug is why am I getting points off of the Z=0 plane at all?
According to the bounding box of the points, I’m getting a point with a MinZ at -2.2204460492503131e-16 and a MaxZ at 2.2204460492503131e-16.
Is there a precision problem with vtkCutter?

Thanks.

2

Interestingly, I just discovered that 2.2204460492503131e-016 is DBL_EPSILON.

3

It seems that this was discovered. A very similar fix to what I described is in VTK 8.2.

// Use a finite tolerance when detecting zero width sides to ensure that
// numerical noise doesn’t cause an explosion later on. We’ll consider any
// length that’s less than 0.1% of the average length to be zero:
double totLen = lengths[0] + lengths[1] + lengths[2];
const double zeroDetectionTolerance = totLen * (0.001 / 3.);

I would argue that 1e-3 is too big. We typically use 1e-6 in similar checks in our FE code. But I’m fairly certain that this will work.

4

The question remains though of why vtkCutter generates points with Z=DBL_EPSILON and Z=-DBL_EPSILON when the implicit function’s plane definition is exact.

Does anyone have any ideas?

5

For posterity, I traced the generation of the DBL_EPSILON value in Z and it happens in
vtkTetra::Contour() as part of the interpolation.
So it looks like numerical noise that has to be handled by ComputeDivisions().
It certainly won’t make a difference to the rendering code that it isn’t exactly 0.