[vtkusers] GenerateValues numerical error bug

Serge Lalonde serge at infolytica.com
Wed Oct 29 16:45:52 EDT 2008


Hello VTK users,

I thought that I'd post this little fix in case others might find it 
useful. The class vtkContourValues (in the Common library) is used to 
generate a set of N linearly spaced values between an upper and lower 
bound. I discovered that its implementation (in 5.2.0 and still in 
5.3.0) caused a numerical error to be introduced in all but the lower 
bound value. The effect being that values that should be exact would 
have a very small error introduced. For example, instead of 0, 0.1, 0.2, 
0.3, ..., 0.8, 0.9, 1 you would get 0, 0.1, 0.2, 0.3000000000013, 
0.40000000000026, ..., 0.8000000000045, 0.99999999999993 (or something 
very similar).

The problem is that the loop to calculate those values was 
pre-calculating the delta and then incrementing the value before each 
call to SetValue(). The more accurate method is to calculate the value 
from scratch each time so that the FPU can properly manage and round-off 
the errors in the calculation caused by the (finite) double-precision 
arithmetic. With this fix, the values are exact.

Here is the corrected version of GenerateValues() in vtkContourValues.cxx.

// Generate numContours equally spaced contour values between specified
// range. Contour values will include min/max range values.
// The value is calculated using the lower range values and the loop 
variable
// at each iteration to avoid accumulating the floating point error that
// happens when using using an intermediate calculation and/or when 
incrementing.
void vtkContourValues::GenerateValues(int numContours, double range[2])
{
  int i;

  this->SetNumberOfContours(numContours);
  if (numContours == 1)
    {
    this->SetValue(0, range[0]);
    }
  else
    {
    for (i= 0; i < numContours; ++i)
      {
      this->SetValue(i, range[0] + i*(range[1] - range[0])/(numContours 
- 1));
      }
    }
}

I will submit this change to the bug tracker as well so that it can 
(hopefully) be merged into the development branch.



More information about the vtkusers mailing list