MantisBT - VTK
View Issue Details
IDProjectCategoryView StatusDate SubmittedLast Update
0007887VTK(No Category)public2008-10-29 16:502008-10-31 15:06
ReporterSerge Lalonde 
Assigned ToFrançois Bertel 
PrioritynormalSeverityminorReproducibilityalways
StatusclosedResolutionfixed 
PlatformOSOS Version
Product Version 
Target VersionFixed in Version 
Project
Type
Summary0007887: Calculated values from vtkContourValues::GenerateValues() contain round-off errors
DescriptionThe 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));
     }
   }
}
Steps To Reproduce
Additional Information
TagsNo tags attached.
Relationships
Attached Files
Issue History
Date ModifiedUsernameFieldChange
2008-10-29 16:50Serge LalondeNew Issue
2008-10-29 16:56François BertelStatusbacklog => tabled
2008-10-29 16:56François BertelAssigned To => François Bertel
2008-10-29 16:57François BertelNote Added: 0013978
2008-10-29 17:06François BertelNote Added: 0013979
2008-10-29 17:06François BertelStatustabled => @80@
2008-10-29 17:06François BertelResolutionopen => fixed
2008-10-30 08:42Serge LalondeNote Added: 0013985
2008-10-30 08:42Serge LalondeStatus@80@ => @20@
2008-10-30 08:42Serge LalondeResolutionfixed => reopened
2008-10-30 10:07François BertelNote Added: 0013987
2008-10-30 11:34Serge LalondeNote Added: 0013988
2008-10-30 12:42François BertelNote Added: 0013989
2008-10-30 12:42François BertelNote Added: 0013990
2008-10-30 12:42François BertelStatus@20@ => @80@
2008-10-30 12:42François BertelResolutionreopened => fixed
2008-10-31 15:06François BertelStatus@80@ => closed
2011-06-16 13:11Zack GalbreathCategory => (No Category)

Notes
(0013978)
François Bertel   
2008-10-29 16:57   
Make sense. Thanks.
(0013979)
François Bertel   
2008-10-29 17:06   
Fixed in revision 1.21 of VTK/Common/vtkContourValues.cxx

There is one difference with the actual patch and proposed patch:

(range[1] - range[0])/(numContours - 1) has been saved in a local variable out of
the loop.
(0013985)
Serge Lalonde   
2008-10-30 08:42   
Hi Francois,

It's nice to see such a fast turnaround. One problem though. Unfortunately, you can't store
   (range[1] - range[0])/(numContours - 1)
in a local variable because that will prevent the calculation from being exact. I know because that was the first change that I tried. It is also why I put that comment at the start of the GenerateValues() code in the bug description. The entire expression must be evaluated from scratch for the round off errors to be evaluated properly by the FPU. I know that it goes against normal programming methodology, but in numerical programming such things are sometimes needed.
(0013987)
François Bertel   
2008-10-30 10:07   
OK, I'll change that.

Just curious, on which platform did you notice this effect?

Is it on Linux, on x86 (not x86_64)? If so, here is the explanation
(default mode is extended precision on Linux x86 instead of double-precision rounding mode)

http://www.network-theory.co.uk/docs/gccintro/gccintro_70.html [^]
(0013988)
Serge Lalonde   
2008-10-30 11:34   
I tested on Windows using Visual Studio 2005 SP1. I built the DLLs using the Win32 32-bit platform configuration, but my OS is actually Windows XP x64 Professional (but I don't think that the 64-bit OS is a problem).
The extended precision issues have been an issue for us in the past so we try to remain vigilant about them. The VTK builds use the same floating point model (fp:precise) as we do in our own code.
Thanks.
(0013989)
François Bertel   
2008-10-30 12:42   
You are right, if you are using /fp:precise should get rid of the extended vs double precision mismatch. I just adding this quote for reference:

"With /fp:precise on x86 processors, the compiler will perform rounding on variables of type float to the proper precision for assignments and casts and when passing parameters to a function. This rounding guarantees that the data does not retain any significance greater than the capacity of its type."

ref: http://msdn.microsoft.com/en-us/library/e7s85ffb(VS.80).aspx [^]
(0013990)
François Bertel   
2008-10-30 12:42   
Really fixed in revision 1.22.