Floating point precission issues.

1,512

Author

February 14, 2016 09:46 PM

I ran into problem when two vertices that are on diffrent side of a plane (n,d) do not produce intersection by float point precission issue:














template <class T>
bool  SegmentPlaneIntersection(t3dpoint<T> n, T d, t3dpoint<T> rA, t3dpoint<T> rB, t3dpoint<T> & res)
{
t3dpoint<T> PointOnPlane;
long double originDistance = (long double)d;
long double		distance1 = ((n.x * rA.x)  + (n.y * rA.y)  + (n.z * rA.z)) + originDistance;
long double		distance2 = ((n.x * rB.x)  + (n.y * rB.y)  + (n.z * rB.z)) + originDistance;

if(distance1 * distance2 > 0.0) return false; //here
t3dpoint<T> vLineDir = Normalize(rB - rA);

long double 	Numerator = -1.0 * (n.x * rA.x + n.y * rA.y + n.z * rA.z + originDistance);
long double 	Denominator = ( (n.x * vLineDir.x) + (n.y * vLineDir.y) + (n.z * vLineDir.z) );

	if( Denominator == 0.0) PointOnPlane = rA;
else 
{
long double		dist = Numerator / Denominator;
  PointOnPlane = (rA + (vLineDir * T(dist)));
}
res = PointOnPlane;
return true;

}

code that is responsible for this is here:












if(distance1 * distance2 > 0.0) return false;

even when one point is lets say 50, 1.1, 50 and second is 50, -0.8, 50 it produces positive dot product (that was odd really

exact points of triangle:












------------------
ADDED WHOLE FACE AS BELOW WATERLINE:
TESTING POINT A: 62.9247131347656  -0.0228385925292969  47.2107238769531  dot: -0.0228385925292969
TESTING POINT B: 40.1480102539063  1.10841751098633  44.3472938537598  dot: 1.10841751098633
TESTING POINT C: 40.1480102539063  -0.918651580810547  49.4782066345215  dot: -0.918651580810547
----------------------

note AB, BC they are clearly between that plane but dot in this function failed gratefully (please do not tell me that its the resaon of c casting i made a function like this: and the result was the same no interesection

[spoiler]














bool  SegmentPlaneIntersectionf(vec3 n, float d, vec3 rA, vec3 rB, vec3 & res)
{

vec3 PointOnPlane;


  float		distance1 = ((n.x * rA.x)  +
				 (n.y * rA.y)  +
				 (n.z * rA.z)) + d;

  float		distance2 = ((n.x * rB.x)  +
				 (n.y * rB.y)  +
				 (n.z * rB.z)) + d;

	if (distance1 * distance2 > 0.0) return false;


 vec3 vLineDir = Normalize(rB - rA);


	float 	Numerator = -( dot(n, rA) + d);




 float 	Denominator = dot(n, vLineDir);

if( Denominator == 0.0)
PointOnPlane = rA;

else
{
float		dist = Numerator / Denominator;
  PointOnPlane = (rA + (vLineDir * dist));
}
res = PointOnPlane;
return true;

}

[/spoiler]

Now to see maybe why is this happening i went to function that calls it: (its a long function so i put parts of it)

full

[spoiler]














void SliceHullByWater(TWaterWaves * ocean)
{



int arri = -1;

		TPolygon<float> sliced_area;
		float surface_area = 0.0;
		t3dpoint<float> A;
		t3dpoint<float> B;
		t3dpoint<float> output;
for (int f=0; f < hull->VBO_BE.Length; f++)
{
float dst;
vec3 pop;
vec3 normal = vec3(0.0, 1.0, 0.0);//FindOceanTriangleNormal(pos + ROTATION_MAT * hull->FACE_CENTER_POINT[f], ocean, 0.0, dst, pop);
dst = getplaneD(normal, vec3(0.0,0.0,0.0));
 int ahue = hull->VBO_BE[f].INDEX_START;
 int flen = hull->VBO_BE[f].length;


for (int i=ahue; i < ahue+flen; i++) //through all vertices in that poly..., 3 only
{
sliced_area.Count = 0;
		A = pos + ROTATION_MAT * hull->VBO_V[i];

		int next = i + 1;
		if (next >= ahue+flen) next = ahue;

		B = pos + ROTATION_MAT * hull->VBO_V[next];

		if ( SegmentPlaneIntersectionf(normal, dst, A, B, output) )
			sliced_area.AddVertex( output );


}


	//~~~~~~~~~~ Test if Area is fully submerged or not at all
		if (sliced_area.Count == 0)
		{
		if (dot( pos + ROTATION_MAT * hull->VBO_V[ahue], normal) + dst < 0.0) //fully submerged triangle
		{
		arri = arri + 1;
		ACTUAL_FRAME[arri].V[0] = pos + ROTATION_MAT * hull->VBO_V[ahue];
		ACTUAL_FRAME[arri].V[1] = pos + ROTATION_MAT * hull->VBO_V[ahue+1];
		ACTUAL_FRAME[arri].V[2] = pos + ROTATION_MAT * hull->VBO_V[ahue+2];
		ACTUAL_FRAME[arri].normal 	= normal;
		ACTUAL_FRAME[arri].dst 		= dst;
		ACTUAL_FRAME[arri].pop 		= pop;

		ALOG(""); ALOG("");ALOG("------------------");
		ALOG("ADDED WHOLE FACE AS BELOW WATERLINE:");
		ALOG("TESTING POINT: "+POINT_TO_TEXT(ACTUAL_FRAME[arri].V[0]) + "  dot: " +FloatToStr(dot( ACTUAL_FRAME[arri].V[0], normal) + dst));
		ALOG("TESTING POINT: "+POINT_TO_TEXT(ACTUAL_FRAME[arri].V[1]) + "  dot: " +FloatToStr(dot( ACTUAL_FRAME[arri].V[1], normal) + dst));
		ALOG("TESTING POINT: "+POINT_TO_TEXT(ACTUAL_FRAME[arri].V[2]) + "  dot: " +FloatToStr(dot( ACTUAL_FRAME[arri].V[2], normal) + dst));
		ALOG("----------------------");
		ALOG("");
		}
		continue;   //go to next face
		}
	//**************************if partially submerged

	//see how many vertices are behind to know what type of cut we have
	int backfaces = 0;
 ALOG("");ALOG("");
for (int i=ahue; i < ahue + flen; i++)
{
ALOG("TESTING POINT: "+POINT_TO_TEXT(pos + ROTATION_MAT * hull->VBO_V[i]) + "  dot: " +FloatToStr(dot( pos + ROTATION_MAT * hull->VBO_V[i], normal) + dst));

if (dot( pos + ROTATION_MAT * hull->VBO_V[i], normal) + dst < 0.0)
{
//	if (classify_point_plane( pos + ROTATION_MAT * hull->VBO_V[i], normal, dst) == isBack)
ALOG("ADDING AS BACKFACE");
		backfaces = backfaces + 1;
}
}
//see how many vertices are behind to know what type of cut we have

	sliced_area.Count = 0;

	for (int i=ahue; i < ahue + flen; i++)
	{

		A = pos + ROTATION_MAT * hull->VBO_V[i];

		if (dot( A, normal ) + dst < 0.0)
		{
		sliced_area.AddVertex( A );
		ALOG("ADDING POINT: " +POINT_TO_TEXT(A));
		}


		int next = i + 1;
		if (next >= ahue+flen) next = ahue;

		B = pos + ROTATION_MAT * hull->VBO_V[next];

		if ( SegmentPlaneIntersection(normal, dst, A, B, output) )
		{
		sliced_area.AddVertex( output );
		ALOG("ADDING SLICER POINT: " +POINT_TO_TEXT(output));
		}
		else
		{
		 float a = dot( A, normal ) + dst;
		 float b = dot( B, normal ) + dst;
		 if (!( ( (a <= 0) && (b <= 0) ) || ( (a > 0) && (b > 0) ) ))
		 ShowMessage("Definetly a hit but i cant see it");

		}

	}

	vec3 cp;

	if (backfaces == 1)
	{


			arri = arri + 1;
		ACTUAL_FRAME[arri].V[0] = sliced_area.V[0];
		ACTUAL_FRAME[arri].V[1] = sliced_area.V[1];
		ACTUAL_FRAME[arri].V[2] = sliced_area.V[2];
		ACTUAL_FRAME[arri].normal 	= normal;
		ACTUAL_FRAME[arri].dst 		= dst;
		ACTUAL_FRAME[arri].pop 		= pop;
   ALOG("FIRST TRI: "+IntToStr(sliced_area.Count));
if (dot(ACTUAL_FRAME[arri].V[0],normal) + dst > 0) ALOG("Vert 0: above waterline");
if (dot(ACTUAL_FRAME[arri].V[1],normal) + dst > 0) ALOG("Vert 1: above waterline");
if (dot(ACTUAL_FRAME[arri].V[2],normal) + dst > 0) ALOG("Vert 2: above waterline");


	} else //two triangle vertices under the water
	{
	if (backfaces < 2)
	{
	ShowMessage("DUD KURWA WTF!");
	}
				arri = arri + 1;
		ACTUAL_FRAME[arri].V[0] = sliced_area.V[0];
		ACTUAL_FRAME[arri].V[1] = sliced_area.V[1];
		ACTUAL_FRAME[arri].V[2] = sliced_area.V[2];
		ACTUAL_FRAME[arri].normal 	= normal;
		ACTUAL_FRAME[arri].dst 		= dst;
		ACTUAL_FRAME[arri].pop 		= pop;
		ALOG("FIRST TRI");
if (dot(ACTUAL_FRAME[arri].V[0],normal) + dst > 0) ALOG("Vert 0: above waterline");
if (dot(ACTUAL_FRAME[arri].V[1],normal) + dst > 0) ALOG("Vert 1: above waterline");
if (dot(ACTUAL_FRAME[arri].V[2],normal) + dst > 0) ALOG("Vert 2: above waterline");


        			arri = arri + 1;
		ACTUAL_FRAME[arri].V[0] = sliced_area.V[2];
		ACTUAL_FRAME[arri].V[1] = sliced_area.V[3];
		ACTUAL_FRAME[arri].V[2] = sliced_area.V[0];
		ACTUAL_FRAME[arri].normal 	= normal;
		ACTUAL_FRAME[arri].dst 		= dst;
		ACTUAL_FRAME[arri].pop 		= pop;
		ALOG("SECOND TRI");
if (dot(ACTUAL_FRAME[arri].V[0],normal) + dst > 0) ALOG("Vert 0: above waterline");
if (dot(ACTUAL_FRAME[arri].V[1],normal) + dst > 0) ALOG("Vert 1: above waterline");
if (dot(ACTUAL_FRAME[arri].V[2],normal) + dst > 0) ALOG("Vert 2: above waterline");
	}


}




  submerged_tri_count = arri + 1;
}

[/spoiler]

now to the code that was failing:

click for description:

[spoiler]

Function name is: void SliceHullByWater(TWaterWaves * ocean) so it slices hull model with water planes (but for simplification i used one plane for all triangles

- Ok so first of all I loop through all faces in the model

- then I go through all vertices in that model and check whenever they are on different sides of plane: if they are all on one side of plane we check if first vertex of the face is below that plane (water plane) if it is then we add this triangle to the stack. (ACTUAL_FRAME), then skip this face iteration. If not then we do nothing and skip this face iteration (because triangle is above water and we return only submerged part)

[/spoiler]










Ahue is index in vertex array (index of first vertex of that face, flen is face length)

for (int i=ahue; i < ahue+flen; i++) //through all vertices in that poly..., 3 only
{
sliced_area.Count = 0;
		A = pos + ROTATION_MAT * hull->VBO_V[i];

		int next = i + 1;
		if (next >= ahue+flen) next = ahue;

		B = pos + ROTATION_MAT * hull->VBO_V[next];
//here we get vector AB (for each side of triangle)

//Then we test if theres intersection if there is add it to some temp polygon
//Like I am trying to explain it fails here like crap
		if ( SegmentPlaneIntersectionf(normal, dst, A, B, output) )
			sliced_area.AddVertex( output );


}

Now If there was no intersections at all 
	//~~~~~~~~~~ Test if Area is fully submerged or not at all
		if (sliced_area.Count == 0)
		{
		if (dot( pos + ROTATION_MAT * hull->VBO_V[ahue], normal) + dst < 0.0) //fully submerged triangle it test first vertex of face (if its below waterline)
		{
		arri = arri + 1;
		ACTUAL_FRAME[arri].V[0] = pos + ROTATION_MAT * hull->VBO_V[ahue];
		ACTUAL_FRAME[arri].V[1] = pos + ROTATION_MAT * hull->VBO_V[ahue+1];
		ACTUAL_FRAME[arri].V[2] = pos + ROTATION_MAT * hull->VBO_V[ahue+2];
		ACTUAL_FRAME[arri].normal 	= normal;
		ACTUAL_FRAME[arri].dst 		= dst;
		ACTUAL_FRAME[arri].pop 		= pop;


Log which you saw above tells that first vertex in this trinagle
hull->VBO_V[ahue] is below waterline so that’s why i posted the log

		ALOG(""); ALOG("");ALOG("------------------");
		ALOG("ADDED WHOLE FACE AS BELOW WATERLINE:");
		ALOG("TESTING POINT: "+POINT_TO_TEXT(ACTUAL_FRAME[arri].V[0]) + "  dot: " +FloatToStr(dot( ACTUAL_FRAME[arri].V[0], normal) + dst));
		ALOG("TESTING POINT: "+POINT_TO_TEXT(ACTUAL_FRAME[arri].V[1]) + "  dot: " +FloatToStr(dot( ACTUAL_FRAME[arri].V[1], normal) + dst));
		ALOG("TESTING POINT: "+POINT_TO_TEXT(ACTUAL_FRAME[arri].V[2]) + "  dot: " +FloatToStr(dot( ACTUAL_FRAME[arri].V[2], normal) + dst));
		ALOG("----------------------");
		ALOG("");
		}
		continue;   //go to next face
		}

and i got this note that triangle 0,1,2 is actually this triangle:

and the result:

SegmentPlaneIntersection doesnt see anything for 0,1,2 vertices so sliced_area.Count == 0

then it tests dot(vertex 0, waterline_normal) + waterline_dst < 0 since it is below waterline it adds whole triangle to submerged triangle list (idea was good but SegmentPlaneIntersection fails in this part of code or maybe the problem is somewhere else) but then SegmentPlaneIntersection should work...

the hting is in the second part of full function o do the same loop with same code but second one works but not this first. D:!

after thinking why the hell (dot(a,n)+d ) * (dot(b,n)+d) fails i come with different solution to the problem (note that this dot * dot thing is very important here because i willl show you that it may work in another part of the function.

here is fixed part of code:












bool fully_submerged = true;
for (int i=ahue; i < ahue+flen; i++) //through all vertices in that poly..., 3 only
{

		A = pos + ROTATION_MAT * hull->VBO_V[i];
   
		int next = i + 1;
		if (next >= ahue+flen) next = ahue;
		
		B = pos + ROTATION_MAT * hull->VBO_V[next];

	  int aside = int(PLUSORMINUS(dot(A, normal) + dst));
	  int bside = int(PLUSORMINUS(dot(B, normal) + dst));
if ((aside*bside) < 0) fully_submerged = false;


}




	//~~~~~~~~~~ Test if Area is fully submerged or not at all

		if (fully_submerged) //fully submerged triangle
	   if ((dot(pos + ROTATION_MAT * hull->VBO_V[ahue], normal) + dst) < 0)
		{
		arri = arri + 1;
		ACTUAL_FRAME[arri].V[0] = pos + ROTATION_MAT * hull->VBO_V[ahue];
		ACTUAL_FRAME[arri].V[1] = pos + ROTATION_MAT * hull->VBO_V[ahue+1];
		ACTUAL_FRAME[arri].V[2] = pos + ROTATION_MAT * hull->VBO_V[ahue+2];
		ACTUAL_FRAME[arri].normal 	= normal;
		ACTUAL_FRAME[arri].dst 		= dst;
		ACTUAL_FRAME[arri].pop 		= pop;

		ALOG(""); ALOG("");ALOG("------------------");
		ALOG("ADDED WHOLE FACE AS BELOW WATERLINE:");
		ALOG("TESTING POINT: "+POINT_TO_TEXT(ACTUAL_FRAME[arri].V[0]) + "  dot: " +FloatToStr(dot( ACTUAL_FRAME[arri].V[0], normal) + dst));
		ALOG("TESTING POINT: "+POINT_TO_TEXT(ACTUAL_FRAME[arri].V[1]) + "  dot: " +FloatToStr(dot( ACTUAL_FRAME[arri].V[1], normal) + dst));
		ALOG("TESTING POINT: "+POINT_TO_TEXT(ACTUAL_FRAME[arri].V[2]) + "  dot: " +FloatToStr(dot( ACTUAL_FRAME[arri].V[2], normal) + dst));
		ALOG("----------------------");
		ALOG("");
				continue;   //go to next face
		} else continue;

where












template <class type> type  PLUSORMINUS(type k)
{
	if (k >= 0) return type(1.5); else return -type(1.5);
}

i basically checked dots and whenever dot was 0 or greater i returned integer = 1

if it was below i returned integer -1

so multiplying these two will always give me proper positive / negative values

it gave me this:

now theres important part of code i didint mention earlier the one that uses SegmentPlaneIntersection with crappy dot*dot thing.

but before i post part of the code and whole function code i want to explain what whole function does and what does that important second part of code:

Function name is: void SliceHullByWater(TWaterWaves * ocean) so it slices hull model with water planes (but for simplification i used one plane for all triangles

- Ok so first of all I loop through all faces in the model

- then I go through all vertices in that model and check whenever they are on different sides of plane: if they are all on one side of plane we check if first vertex of the face is below that plane (water plane) if it is then we add this triangle to the stack. (ACTUAL_FRAME), then skip this face iteration. If not then we do nothing and skip this face iteration (because triangle is above water and we return only submerged part)

But whenever we have points on both sides of water plane then we need to perform cut.

And here is the second part I use that same crappy dot*dot function that determines whenever theres intersection or not.












     sliced_area.Count = 0;

            for (int i=ahue; i < ahue + flen; i++) //throug all face vertices
            {

                        A = pos + ROTATION_MAT * hull->VBO_V[i];

                        if (dot( A, normal ) + dst < 0.0) // add to temp polygon as new underwater triangle vertex
                        sliced_area.AddVertex( A );


                        int next = i + 1;
                        if (next >= ahue+flen) next = ahue;

                        B = pos + ROTATION_MAT * hull->VBO_V[next];

//here check for intersection if theres intersection  then put that to output var and add that output to as another vertex to temp polygon
                        if ( SegmentPlaneIntersection(normal, dst, A, B, output) )
                        sliced_area.AddVertex( output );
            }

Like you see I am again using SegmentPlaneIntersection but now it works even when I test same sides, it just works....

now i am really confused coz i spent last 3 days thinking why the hell in first loop SegmentPlaneIntersection function didin't work and in this part of code is working. (my first thoughts are because of some floating point precision issues - since i dont see any bugs in original function (that whole code i posted in spoiler) i am not quite sure what is going on and this is my question is it really a FP issue or i just simply wrote (first part) it wrong? and i dont see some statement that causes magic to happen..

Im confused because it may work on windows machine but when i put that into android phone - depending on phone it can slice properly or wrongly liek in the first picture.

Master & Mentor
https://sites.google.com/site/customprog/

JoeJ

4,183

February 15, 2016 12:08 PM

Did not read through this all

The line segement plane intersection seems more complex than necessary (e.g. normalization).
Below the code i'm using without sqrt.

You wanna do Buoyancy simulation and you clip the whole detailed ship by the waterplane?
Why don't you use a low res approximization, like a box? Should be good enough?


inline float IntersectRayPlane (qVec3 &rO, qVec3 &rD, qVec3 &pO, qVec3 &pN)
	{
		// rD does not need to be unit length
		float d = pN.Dot(rD);
		float n = pN.Dot(pO - rO);

		if (fabs(d) < FP_EPSILON) // ray parallel to plane
		{           
			if (fabs(n) < FP_EPSILON) return 0; // ray in plane
			else return FLT_MAX; // no intersection
		}
		float t = n / d;
		//qVec3 intersectionPoint = rP + rD * t;
		return t;
	}

_WeirdCat_

1,512

Author

February 15, 2016 12:24 PM

when using boxes i cant calculate drag acting on a boat properly without using some imaginary coefficients that act on specified angle at specified power, additionally it was flipping the boat during turns(but thats not related to box thing) any way slicing box and hull uses the same algorithm

Master & Mentor
https://sites.google.com/site/customprog/

JoeJ

4,183

February 15, 2016 12:44 PM

if( Denominator == 0.0) PointOnPlane = rA;
else

That's very very bad - you compare a floating point number against zero.
You should test against a small epsilon instead, like i do here: if (fabs(d) < FP_EPSILON)
(FP_EPSILON = 1.0e-5f for single precision float)

The problem is, a tiny number larger than zero but smaller than a reasonable small number like 0.00001
may result in garbage after the division - testing for exactly zero is not enough.

There is almost never a reason to compere floating point against a exact value.
You need to use epsilons everywhere.

_WeirdCat_

1,512

Author

February 15, 2016 02:48 PM

main reason: if(distance1 * distance2 > 0.0) return false;

Master & Mentor
https://sites.google.com/site/customprog/

JoeJ

4,183

February 15, 2016 03:43 PM

main reason: if(distance1 * distance2 > 0.0) return false;

That's the same bad practice. First you should classify if the points are ON the plane (fabs(distance) < epsilon).
If one point is on the plane you can return, and afterwards your test would be ok, because the signs are guaranteed to make sense.
For robust polygon clipping you need to handle special cases like point on plane, using doubles won't fix it.

I can't tell if this causes your actual problem because you forgot to give numbers for the plane.
However - seems you're not willing to listen anyways.

Randy Gaul

2,835

February 15, 2016 08:08 PM

You should also look up Sutherland Hodgman clipping inside of Christer Ericson's "Real-Time Collision Detection" book. He has some text that explains planes, thick planes, clipping and how to avoid common bugs/problems like you are running into. It seems you are missing some fundamentals of using floats for geometry stuff on the CPU, and these links and this book will help clear that up. Like Joe was saying if we deal with planes with floats we must think in terms of "close enough" with error tolerances. These resources are great reads and directly help with what you are working on

Bouyancy demo (exactly what you are doing): http://box2d.org/files/Misc/Buoyancy.zip

Float tolerances: http://realtimecollisiondetection.net/blog/?p=89

_WeirdCat_

1,512

Author

February 15, 2016 11:16 PM

i saw that buoyancy.zip some time ago, it actually isn't the same thing because i clip whole model then make tetrahedrons where all of them share same apex in center of mass and then i clip these tetrahedrons then i compute center of buoyancy for each tetrahedron (or a clipped solid) (by adding vertices of each shape and dividing by vert amount) not quite sure if thats the way of doing it but, from what i understand i need a mass of every tetrahedron (even when sliced to another shape - basically had to compute a mass at wich i apply a force thats a bit stupid from now since i began to understand that after computing inertia matrix i would only add moments as: Moments = moments + vector_buo_center_2_center_of_mass x force, so i wont have to know the mass of each sliced tetrahedron/solid. but i figoured that after making that ;x

now i am a bit confused if i should put apex of all tetrahedrons in a water plane and calc each tetrahedroon buoyancy force or to calc whole buoyant force for them and somehow apply torque without specyfying point of contact (where force is applied like in this demo but it uses the mass of the submerged part)

or do it like i do: every tetrahedron for each submerged triangle share on apex (at hulls cetner of mass) then i slice again all tetrahedrons to properly calculate volume of submerged part of a hull. after having list of solids that are inside ship hull and they are below waterline i compute buoyant force for each one (and applying torque (after specifying buoyancy center as each solid geometric center point) i will also haveto add form drag to each submerged triangle of a hull that i slice but ill post that after i finish calculation of inertia matrix and fixing whole rotational part of the 'rigid body of mine' - basically i would love to know how to calculate form drag (as far as i remember normalized(-(linear vel + ang vel)) was the initial dragforce direction, since i have a shape i could do that even when it would be a really retarded approximation but at least in that case i won't have to use some hardcoded drag coefficients.

Master & Mentor
https://sites.google.com/site/customprog/

Randy Gaul

2,835

February 16, 2016 02:26 AM

I think you're overlooking the Catto bouyancy demo. Did you read his article in game programming gems 6? He's doing the same thing you are, and speaks about where to place the apex. He also covers drag.

And, if you need mass you can calculate the volume of your mesh and multiply by an appropriate density value, or just assign a custom mass yourself.