Hi Luke,
Fernando Cacciola wrote:
The library proposed by Lucannus, OTOH, at least acknowledges the problem appropriately and approaches a solution in the right direction (there is still some important design-influencing considerations but are not showstopers)
Thank you, Fernando. I agree with you 100% that robustness needs be designed in up front and not tacked onto an existing geometry code base as an afterthought. In particular, the robustness strategy needs to be comprehensive and not patchwork. I acknowledge that my library needs a lot of work in this area to make it robust for both integer and floating point, particularly since it was only designed with integer robustness in mind. I will have to redesign the low level stuff quite significantly to get both integer and floating point fully robust and sharing the same algorithm. I don't really have the knowledge I need to make my library floating point robust. I need to read some of the papers you mention.
I'll try to sketch the general ideas below to get you started.
Making the coordinate data type a template parameter and "allowing" the user to specify infinite precision numerical data type if they want robustness is not a solution
Precisely :) life would be too easy it that very already fast enough.
, because there are very simple and effective ways to make integer and floating point operations robust. I don't know if we need to design in points of customization to expose the ability to do that to the user, but it needs to be at least possible for the library developer to specify different strategies for different coordinate data types and different contexts.
Exactly... and the library as a service to other higher level libraries should and can provide a general methodology and the corresponding tools for doing that. Ultimately users shouldn't be concerned at all with all this, except when there is reasonable speed/robustness tradeoff that users would like to threshold themselves.
I think Brandon's library is trying to handle both integer and floating point robustness, and I've noticed he is at least as interested in multi-precision numerical data types as I am, which is a good sign that he's doing the right things. I'm interested in learning more about his approach.
As I'll try to sketch below, bigfloat are part of the toolset, but there is more to it.
Performance is fine as a goal, but if the algorithms fail (crash, hang, drop polygons or "gracefully" throw an exception) even a small percentage of the time due to robustness issues then they are worthless, regardless of how fast.
Exactly.
If we tried to use a non-robust library and it failed in a mission critical situation it costs us over one million dollars a day to fix the problem. Failure is not an option. When I say 100% integer robust I'm serious. Its all or nothing and the stakes for me are high.
And so are for end users.. it is the developers that think: "this works in all cases I tried, then that's good enough". But it isn't, and there is no need to compromise since lots of techniques exists and while they don't cover all of the possible algorithms, there is clearly little justification for not using the techniques where they are applicable. === LONG EXPLANTION BELOW === OK, I'll try to sketch the general idea behid floating point robustness in geometric computing. There is a lot on the subject on the literature, is just difficult to find, so consider this just an intrduction. On the one hand, if an algorithm only needs to perform algebraic operations, a "Bigfloat", i.e. a software-based adaptive unlimited precision floating point type, delivers totally exact results. If division is needed then you can use a rational of bifgloats. If other operations are needed, such as sqrt() or transcendentals (say cos, sin) they you can use a number type that tracks the expression trees and avoids actually performing the non-exact operations as much as possible, specially when you query for discrete properties such as its sign. On the other hand, it is possible to design and implement any geometric algorithm in terms of "predicates" and "constructions" and it is possible to consider the output as having separate topological and geometric parts. Let's consider an example: the straight skeleton (and the voronoi diagram in this case) of a square contains one and only one vertex in the center and one and only one edge from that center node to each square corner. That result has a topology (how vertices and skeleton nodes are connected) and a geometry (the coordinates of the center node). The topology might be correct even if the geometry is not: i.e. the center node is not exactly in the center, yet they are connected correctly. For many many geometric applications, it is the topology what needs to be correct, not the geometry. If you do clipping for instance, it might very well not matter at all (for instance in drawing applications) if the vertices are slightly off. But if an entire region is included in the result when it should not, then that is a disaster for whatever application. Since the input to an algorithm is error free no matter the representation (or more precisely: there is nothing the algorithm can do about its implicit error), it is easy to see that, if an algoritm is designed in such a way that the topological construction is based exclusively on the results of predicates called exclusively on the input (and never ever on an intermediate result), then the topology would be correct no matter what capabilities the underlying input arithmetic has, provided the predicates can do whatever magic is neccesary to produce exact and correct results given input of any type. That is, if there is a predicate that can take machine integer or floating point coordinates, and, in whatever way, produce a guaranteed error free result, then you can write the code of an algorithm totally agnostic of the limitations in the type of input coordinates. That is, the agorithm itsef can just pretend the computation is performed in the RealRAM. Researchers spent a lot of time figuring out how to implement exact predicates based on finite precision integer or floating-point input. While the results are limited and somewhat complicated from an engineering POV they do allow for incredibly fast algorithms that take simple native type coordinates (int or double) and produce 100% robust code. But even then, the *concept* of a black box predicate is still useful because there are other programming-specific techniques to provide simple-to-program but fast predicates: the so called "floating point filters" which I'll explan below. IMO, any geometric library should provide a toolbox of exact predicates (which shoud be exact independetly of the input type) and state as a guideline that any algorithm should only make critical decisions based on the result of predicates, and, unless totally impossible, the predicates should never be called with intermediate results. Here's a simple example: T d1 = distance(a,b) T d2 = distance(c,d) if ( d1 < d2 ) blabla that is fundamentally flawed beyond repair because rebustness becomes the resposability of T, and at the end of the day the poor user finds himself wrapping his input into a fancy numeric type just because of that flaw. If OTOH that were written like so: if ( compare_distace(a,b,c,d) ) blabla then the input could still be of any type yet some basic library tool could go to whatever extent is needed to give the correct answer. And if you are like me and love to exploit C++ unsurpreased expressive power, you can even write that like so: distance_result d1 = distance(a,b) distance_result d2 = distance(c,d) if ( d1 < d2 ) blabla where the comparison between distance_result types does the magic of producing an exact answer no matter what the type of input coordinates are. And this would even prevent the algorithms in the library to mistakely do the wrong thing. Following this methodology at large needs support from the bottom layer of the library, i.e. it needs a exact precicates toolbox that can operate on arbitrary input types and always yield correct results, but it also needs a general methodology: i.e. predicates needs to be used istead of "direct" logic and predicates need to be given only the input and never intermediate results (or as I'll explain below, that must be carefully cosidered and instrumented with its own special technique). These guidelines can be enforced by the lirary given the right design (i.e. having "disntace" return "distance_type" which actually stores the input) All that above allows algorithms to be called on arbitraty input types that don't need to carry the weight of providing robustness. i.e. machine integer and floating point numbers can be used. Couldn't be any simpler for the end user. Constructions OTOH are similarly back-boxes like predicates are, but serve different purposes since it is presumed that the algorithm is well written and never uses the results from constructions to feed predicates (or rather they do so only very carefully and by design, and only in those cases when it was impossible to express the algorithm output as a direct function of its input). That is, while the library provides, say, Point p = centroid(a,b,c) that 'p' is never used to drive the decisions on the algorithm (is never passed to a predicate), yet still that centroid can be a higher level facility that *interally* uses whatever it needs to come up with a reasonable approximated result. Whenever you can express an algorithm in terms of predicates and constructions such that predicates only operate on the input and never on the result of constructions, then your algorithm is robust automagically without any requirement on the input type and without the need to to anything at all to workaround numerical issues within the algorithm itself. That's the job of the predicate oracles. If an algorithm just can't call predicates with only the input and needs to step over intermediate results, things get more complicated. However, if all those intermediate results that would be feed into predicates are ONLY produced with "constructions", sound alternatives can be explored. One alternative, the one used the most, is to have the constructions produce higher level abstactions instead of actual numerical results. For instance, they can return the expression trees that refer to the input, and have the predicates evaluate that lazily on demand (so the predicate is still only really operating on the input). Again, this needs clever support from the library at design type in tgat, for example, the result from "centroid" would't be a point but a centroid_result, storing the input points as well as an approximated result, which could be perfecttly feed into predicates that could do the right thing having the input at hand. This works reasonably well given the inherent complexity but is still the subject of active research. Most of the time the use of DAGs for results is not really needed, so the ideal framework would be able to mix both native input and lazy expressions in a coherent manner. One starting idea (that I had many years ago but never had the time to evolve) would be to have a hybrid geometric type that uses a dynamic dual representation: input type or DAG. Now a word on predicates. In C++ it is actually somewhat trivial to have a blackbox predicate than can produce correct results in spite of the input type using the following technique known as a floating-point filter (and would be equally trivial *within* a virtual machine for any dynamic language FWIW): template<class T> bool pred( T a, T b ) { interval<T> ia = a ; interval<T> ib = b ; interval<T> r = some operations on ia and ib if ( certainly( r == 0 ) ) return true ; else if ( certainly( r != 0 ) ) return false ; else return pred<some_exact_type>(a,b); } That exploits the capability of interval arithmetic to implement certified logic. I won't explain how or why here, but "r == 0" doesn't return bool but TriBool, so you can explicitely query for the "I don't know" answer, and *only* in that case fall back to exact solution. All these machinery must ultimately use some mapping from input type and the interval and exact types needed. Designing that mapping needs to consder different aspects. For instance, not all algoritm require the same exact_type for instance: some might need an exact sqrt() for example. So, the mapping could need to be something like this: coordinate_trats<T>::template exact<requiremets_tag>::type Also, predicates and constructions serve different purposes so they also have differet requirement on the exact types. Going back to the previous straight skeleton example, while the topology must be right, the coordinates of the center nodes need not, BUT, a user could reasonably expect a black box costruction to neverthess come up with a decent approximation, to avoid for instance points totaly off because of floating point overflow. In this case the algorithm should use a "constructor" object that can somehow be told what exact type to use, which might be different than the one used for predicates. Thus, the traits might need to be more like this: coordinate_trats<T> ::template exact<for_predicate_tag, requiremets_tag>::type Best -- Fernando Cacciola SciSoft Consulting, Founder http://www.scisoft-consulting.com