Blender uses (single precision) float for coordinates and pretty much everything else that can have fractional values. The BLI_math and BLI_math_… routines are mostly routines that take floating point arguments, do floating point calculations, and return floating point results. There are a few integer versions of routines and a very few double versions of routines.
A question for discussion, in two parts:
A) Should everything change from float to double?
B) If not, are there other places in the code where calculations could or should be done in double precision, or even exact multiprecision rational arithmetic?
Some thoughts I have on these questions.
Question A: everything double?
Some advantages of this:
a) Occasionally users complain that for very big scenes, the difference in scale between the smallest size of features and the size of coordinates to position things in that space lead to problems when we only have floating point precision. With single precision we can only really deal with about 6 decimal orders of magnitude in scale difference, whereas with double precision we could deal with about 15 orders of magnitude.
b) Similar problems that come into play depending on the view plane vs the 3d space scales. As an example, numerous bugs have been reported with the knife tool that eventually turned out to be related to limited precision problems, and got fixed with very careful epsilon tweaking (that is, making comparisons against some small “epsilon” instead of zero, and adjusting the size of epsilon until things “work”).
c) Presumably other 3d content creation packages use doubles (anyone know which?) and there will be a precision loss when exporting from them and importing into Blender - possibly causing geometric anomalies in Blender when none existed in the original.
d) We could stop having to use “f” suffixes on all our floating point constants, and stop needing (float) casts when using some standard library routine that only deals with doubles
e) Some external packages we’ve incorporated into Blender use doubles internally. This incurs costs like: conversion back and forth; extra storage to store double versions of our data in order to pass to those packages; losing the full advantage of the precision of those calculations by forcing the answers back into floats.
Some disadvantages I see (I’m sure there are more: discuss)
a) Large effort to switch. First, mechanical changes like float → double, sinf → sin, glVertex3f → glVertex3d, 0.0f → 0.0, etc. Probably a script can do that for the most part. Then harder things like making sure that nobody did arithmetic or union/punning or alignment assumptions that assumed float sizes and alignment. Then the even harder problem of finding all the “epsilon tweaking” parts of the code and trying to reverse engineer why the existing epsilons were chosen and what should change if precision is double. And finally, the enormous burden on external users of Blender’s API: all the extant plugins, external renderers, etc. likely have dependencies on the existing float interface. One could perhaps shim that for a while with converters to the float-based API, but long term that is an awful solution.
b) Non-backward-compatible change to .blend in a fairly big way. Maybe OK for 2.8, but either need massive do_version support or else need to provide a way to convert all existing pre-2.8 .blend files into new format.
c) Increased memory usage: probably almost doubles the size of needed main memory for large models. Similarly for disk memory.
d) Performance penalty? I’m less sure here; would be good to see some benchmark comparisons. Maybe the biggest effect is just that more bandwidth between cpu-memory-gpu combinations could hurt performance. Less sure about actual computation cost.
I’m guessing disadvantage (a) is the biggest show-stopper and is what probably makes this whole idea a non-starter. But interested in other’s thoughts.
Question B: places in the code where calculations could or should be done in double precision, or even exact multiprecision rational arithmetic?
I came upon this question in thinking about Boolean operations on meshes. Any time one has to calculate the position of new vertices based on geometric and matrix computations over other vertex positions, we get the problem that limited precision arithmetic means that (a) we have to be careful and probably special case things when we get things near but not quite zero - certainly whenever division takes place; and (b) things that are true based on pure mathematical and geometric theorems about algorithms when all calculations are exact may become untrue when arithmetic is done and results are stored only approximately. For instance, one calculation might say a point is inside a given triangle and another mathematically equivalent one may come to the opposite conclusion. It is hard to write algorithms that always work (and don’t crash!) when building on the quicksand of not being able to trust your math and geometry theorems. The Boolean operations case is a notorious example of this. I can find dozens of papers in the literature about how “other implementations fail on these cases because they don’t use robust arithmetic, so here is our way to solve the problem”. Blender’s current BMesh-based boolean code has some careful sorting-vertex and epsilon-tweaked code to try to solve these problems, but it seems likely that we could easily make it fail by putting in some of the failure examples from the literature (usually cases where intersections lead to very small features - like long triangle slivers).
In general, the code that involves intersection (intersect, intersect (boolean), knife, project knife) will all have these problems lurking in the background, waiting to bite us when someone tries some new example.
What approaches do people take to the problems caused by finite precision arithmetic in geometric algorithms? Here is a brief survey:
-
Use exact arithmetic.
This means using a multiprecision package (either multiprecision floats, or multiprecision integers in a rational number package). There are tricks to speed this up: using double calculations first and some technique to decide whether or not one has enough precision or whether one needs to go to multiprecision; lazy evaluation and caching. This approach is used by the most robust Boolean packages reported in the literature. Maybe problem of inconsistent decisions so hard to solve in general that it is necessary to do this. But there are some drawbacks. Speed, obviously, is a big one. It seems like things might go 20x slower and a non-exact floating point solution, even using some of the tricks mentioned early. This is an important issue for Blender, where people want interactive speeds even for operations involving huge models. Another drawback is that it isn’t by itself the complete solution to problem A: when calculated exactly, many cases that are “close enough” will not actually generate coincidences, but instead will create vertices, edges, and faces that are very close to other geometry or very small in length or area. I suppose one would solve this by a post-processing step that dissolves tiny-extent geometry and snaps other very-close geometry to that geometry. Another drawback is just the amount of extra code that would have to be written or pulled into Blender in order to do the exact arithmetic. This is a non-negligeable consideration too. -
Perturb the input so that everything is in general position and such that we aren’t going to run into cases where inconsistent decisions are made.
It may not be easy to figure out how much perturbation is enough for this to be true. Anyway, this breaks the user requirement that we recognize and incorporate special cases into the output. Maybe post-processing can fix, but I have my doubts about that. -
Restrict input precision such and then do arithmetic in fixed floating point precision with enough bits to make every decision properly (e.g., look at how many bits are needed to calculate a determinant given a fixed precision input). I think for this to work, you have to not only fix the number of mantissa bits, but you have to restrict the actual scale so that we are essentially worked with scaled integers. If it is possible to make this work, it shares the problems of exact arithmetic except that it could be a lot faster.
-
Snap all vertex calculations to a grid.
This is kind of a disguised version of the previous, I think. It makes it easy to find coincident vertices (or does it? if snap is by rounding, then you could have very near things snapped to two different sides), but not sure it solves much else. -
Do coincidence finding as part of intersection routines, using epsilon tests. This is what current BMesh intersect does, and indeed, what most code in Blender does. It is easy enough to program, and will find most coincidences (others may have to be found by preprocessing and/or postprocessing). But it is kind of a cross-your-fingers-and-hope approach to the problems of finite precision arithmetic. So could lead to crashes or nonsense results in unforeseen inputs. Maybe we combine this with use double precision arithmetic in code like this where we are more likely to run into trouble.
What are people’s opinions on this? Is this a problem others have run into, and if so, do you have a favorite solution?
I tried making double versions of all the most basic BLI_math… routines and stashed them in a patch - see Blender Archive - developer.blender.org - but it seems annoyingly invasive to our code, and so I am not recommending that we incorporate that patch, at least not until and if we reach consensus that there are a number of places in Blender where we’d like to use double math, while maintaining the core of Blender using floats.
If we wanted to try multiprecision arithmetic, there are some open source packages like GMP (https://gmplib.org/) and GMPFR (http://www.mpfr.org/) layered on top of GMP that we could consider pulling into Blender. I haven’t investigated how much code that would add.