Floating Point Accuracy Problems

Math and Physics Programming

Started by NMO August 02, 2012 02:53 PM

8 comments, last by Koder4Fun 11 years, 8 months ago

220

Author

August 02, 2012 02:53 PM

Hi!

I have been working on a BSP-Tree implementation. It works in simple cases but in more complicated cases the BSP-Tree fails because of rounding errors. I use single floating point precision.
The problem lies in the classification of the polygon:



BSPNode::classifyFace(Polygon* face, Polygon* newFrontFace, Polygon* newBackFace) {

vector<Vertex3f> *vec = face->getVertices();

for(int i = 0; i < face->numOfVertices(); i++) {

  Vertex3f v = (*vec);

  float res = plane.a*v.X + plane.b*v.Y + plane.c*v.Z + plane.d;

  res = isZero(res) ? 0.0f : res;

.

.

.

Description:
v will be tested against the plane equation (ax + by + cz + d = 0) of the polygon. If res < 0 then on one side, if res > 0 then on other side, else on plane. I use a function isZero() to test if the result is near zero (epsilon = 10^-5) to avoid problems near zero. So far so good.

Now the problem is that res sometimes has wrong values. I've added a screenshot to show what I mean. On the screenshot one can clearly see that there is an intersection between the polygons.
But res value for v1, v2 and v3 are very small negative numbers and will be rounded to zero by isZero() function. That is clearly not correct. What could be a solution to deal with those floating point rounding problems?

Thanks!

NMO

220

Author

August 04, 2012 03:45 PM

Really noone an idea? Maybe I should describe it better:
On the picture are two planes with four vertices each. v1 and v2 are clearly in front of other plane, v0 and v3 are in the back. Now because of rounding problems, v1, v2 and v3 are all zero (means on the plane and hence coplanar). This is clearly wrong because it is not possible that that edges v1-v2 AND v2-v3 are coplanar and v3-v0 and v0-v1 are not. So how to deal with such errors?

LorenzoGatti

4,648

August 04, 2012 08:20 PM

v1 and v2 are clearly in front of other plane, v0 and v3 are in the back. Now because of rounding problems, v1, v2 and v3 are all zero

This is clearly wrong: there is no good reason to round them to 0, just like there is no good reason to round to 0 the scalar product in your first post.
Are you sure your numbers fall into a reasonable range, with exponents close to 0 or at least close to each other?

Omae Wa Mou Shindeiru

jefferytitan

2,525

August 04, 2012 10:30 PM

There may be a coding error, but it's probably worth trying that test using doubles and a slightly lower epsilon value to reduce the chance of error. Even if your co-ordinates are floats, using doubles for intermediate calculations can be worthwhile. Also as the previous poster stated, make sure the model is not too large and not off-centre to ensure the best use of the floating point range.

Koder4Fun

147

August 06, 2012 09:22 AM

The IsZero works as espected, you simply need to change the algorithm to classify a polygon to a plane.
You must classify all vertices with res < -epsilon to be back and all vertices with res >= epsilon to be front, the remaining to be on plane.

Now, to check the polygon side, simply count the number of vertices in front and in back of plane:



if (frontVerts != 0 && backVerts != 0) polygon cross plane

if (frontVerts != 0 && backVerts == 0) polygon in front of plane

if (frontVerts == 0 && backVerts != 0) polygon in back of plane

else polygon on plane

A similar aproach must be used for all other routines.
You must not change the dot product value directly becose you need the real value to find the correct on-plane point on a split-polygon rountine.

I don't know the type of BSP tree you use, but, as a general rule, if you use tollerance values you need to implement it on all functions coded to classify, split, clip points/polygons/segments. You need coherent results on al tests. For example, during the construction of the tree, the partition plane selection routine must use these tollerances.

Using tollerance can help you make a smaller tree on complex scenes and throw away roundoff errors; you can search for "thick plane bsp tree" on google to find some interesting resources.

Stay on single float precision, is good enough if well implemented. You only need to find a good tollerance value; this is linked to the scene units, scene detail, and smallest objects/polygons size.

Also, if you can, decompose your scene into structural and detail geometry. Only the structural geometry is processed as a BSP tree, the detail geometry can be handled as a simple mesh and, if you use your tree to do visibility, computation you only need to link these meshes to the containiner tree leaves.

If you have some question, please tell me.

Sorry for the english is not my native language.

Please vote usefull replies.
Marco Sacchi
Coding is a challenge ... but solving problems is the fun part
My Blog - XNA Italian portal

NMO

220

Author

August 06, 2012 11:53 AM

Thank you for you answers guys!
Especially for your help Koder4Fun!

The IsZero works as espected, you simply need to change the algorithm to classify a polygon to a plane.
You must classify all vertices with res < -epsilon to be back and all vertices with res >= epsilon to be front, the remaining to be on plane.[/quote]

This is exactly what my isZero() function looks like with epsilon = 1.0e-4.

Let me show you an example problem case. I use the example vertices from my picture above:
As described above I calculate res like this:
float res = plane.a*v.X + plane.b*v.Y + plane.c*v.Z + plane.d;

res for v0: -0.00025367874 thats okay because v0 is in back (will not be rounded to zero by isZero() which is okay)
res for v1: -1.3646961e-009 will be rounded to zero that should be okay
res for v2: -1.9651565e-009 will be rounded to zero which should also be okay
so edge v0-v1 in the back and v1-v2 is on the plane, so far okay
res for v3: -8.6204658e-005: will be rounded to zero, so v2-v3 also on the plane ERROR, now I have two edges on the plane v1-v2 and v2-v3 that can't be correct
Thus v3 must no be rounded to zero.

BUT

Now, to check the polygon side, simply count the number of vertices in front and in back of plane:
if (frontVerts != 0 && backVerts != 0) polygon cross plane
if (frontVerts != 0 && backVerts == 0) polygon in front of plane
if (frontVerts == 0 && backVerts != 0) polygon in back of plane
else polygon on plane
[/quote]

your code might work, because it allows more than one coincident edges. I can't wait to test it after work.

You must not change the dot product value directly[/quote]
Which dot product do you mean? I only use the dot product to calc the intersection vertex:
if(res1 < 0 && res2 > 0) { // if vertex1 in front and vertex2 behind plane, then split Vector3f v = v2 - v1; float dist = intersect3(&v1, &v2, &plane); Vector3f newPoint = v1 + (v * dist); }

And I don't use my isZero() function here. Do I have to round here too?

I don't know the type of BSP tree you use[/quote]
I use a very similar implementation like this from BSP Tree FAQ: ftp://ftp.sgi.com/other/bspfaq/faq/bspfaq.html#7.txt

At the moment I simply use the first polygon from the list as the selection plane.

Sorry for the english is not my native language.[/quote]
Your English is very good.

Thank you!

NMO

220

Author

August 07, 2012 08:07 PM

Okay now I've read something about thick plane BSP trees and I've implemented the classification function like this:



Now, to check the polygon side, simply count the number of vertices in front and in back of plane:

if (frontVerts != 0 && backVerts != 0) polygon cross plane

if (frontVerts != 0 && backVerts == 0) polygon in front of plane

if (frontVerts == 0 && backVerts != 0) polygon in back of plane

else polygon on plane

The calculation of the splitting point looks like this:



	Vector3f v = v2 - v1;

	float dist = intersect3(&v1, &v2, &plane);

	Vector3f newPoint = v1 + (v * dist);

	back.push_back(newPoint);

	front.push_back(newPoint);

	front.push_back(v2);

	// assert that vertex is on the plane

	float res = plane.a*newPoint.X + plane.b*newPoint.Y + plane.c*newPoint.Z + plane.d;

	assert(isZero(res));

Simple examples work but more complicated examples do not. The problem is now that I don't know how to debug big examples with many vertices. How should I find the bugs?

LorenzoGatti

4,648

August 08, 2012 08:50 AM

Simple examples work but more complicated examples do not. The problem is now that I don't know how to debug big examples with many vertices. How should I find the bugs?

You shouldn't debug "big examples", only small test cases containing the minimum number of polygons that reproduces the problem.
Big examples should be an inspiration for tests: if you display your polygons and BSP tree you should be able to see errors and inspect the polygons that go wrong well enough to reproduce the trouble in a manageable test.
Reflecting on test cases, even before running them, is likely to provide insight into your bugs because you have to figure out how your code behaves. What's special about the data? What's needed to produce the observed bad output? Where does your code treat that case differently from the good ones?

Omae Wa Mou Shindeiru

NMO

220

Author

August 10, 2012 09:58 AM

@LorenzoGatti: You are right, this is the only possibility to find bugs. I already started to create small test cases.

Another BSP tree specific question: Wouldn't it be better to work with triangles instead of polygons? I think triangles should be easier to debug. The only problem I can think of is splitting but this is acutally also not problematic because when i split a triangle by a plane I get a triangle and a 4-gon. The 4-gon can be triangulated in linear time very easily. I drew a picture to show this.

Koder4Fun

147

August 11, 2012 11:12 AM

Which dot product do you mean? I only use the dot product to calc the intersection vertex:

I'm sorry, I intended distance from plane (that is a dot product + (or - depending on the coordinate system) plane distance form origin.

I know the document you refer for bsp and is not so good for a pratical implementation.
You can find on internet the source code of quake tools (qbsp light vis) to build levels; check Quake1 tools that are more simple.

Here you have a working implementation of a "thick plane bsp tree"; to be more speciifc a leafy bsp tree. The polygons are stored only on leafs, the polygons on partition planes are classified by a dot product of the polygon normal against the plane normal (you can get the normal of your plane ax+by+cz+d=0 by buinding a vector [a,b,c]) and are sended to back or front child nodes.

At the moment I simply use the first polygon from the list as the selection plane.

You can select the partition plane with in mind two methods: least splitting or best balancing.
The method you select depend on the use you make of the tree: if tree is used for rendering is better a best balancing approach, for a collision detection tree you can use a least splitting approach.

Another BSP tree specific question: Wouldn't it be better to work with triangles instead of polygons? I think triangles should be easier to debug. The only problem I can think of is splitting but this is acutally also not problematic because when i split a triangle by a plane I get a triangle and a 4-gon. The 4-gon can be triangulated in linear time very easily. I drew a picture to show this.

Triangles are not so good, you add more splits to the bsp creation.
Using planar polygons you can manage better the "thick plane method" and use epsilons values to check the planarity of a polygon.
Keep in mind that floating point round-off error generate non perfect planar polygons after several splits.

Please vote usefull replies.
Marco Sacchi
Coding is a challenge ... but solving problems is the fun part
My Blog - XNA Italian portal

Floating Point Accuracy Problems

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Floating Point Accuracy Problems

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Reticulating splines