Sign in to follow this  
Finalspace

Proper 2d contact generation for speculative contacts

Recommended Posts

Hi there,
 
i am having trouble implementing 2d rigid body physics for box-vs-box or box vs plane using speculative contacts -> http://www.wildbunny.co.uk/blog/2011/03/25/speculative-contacts-an-continuous-collision-engine-approach-part-1/

There is always an impulse going on which results in my box falling on the plane instantly starting spinning and flying around (Will post video showing the behavior when i get home)


I assume that the contact points are not correct, because paul stated we should just use the closest points - but unfortunatly he never explained how to get more than one contact point :-(

So i implemented a hybrid approach, use a single closest point in case of a separation and use clipping to get the contact points for penetration.
 
What i dont know if this is the correct way to this. Do i need two contact points in a separation case as well???

If someone want to have a look how i implemented it, here it is:

Plane vs edge/box:
 
internal Face GetFace(U32 vertexCount, Vec2f *verts, Vec2f normal) {
	S32 firstIndex = 0;
	F32 firstDistance = Dot(verts[0], normal);
	for (U32 vertexIndex = 1; vertexIndex < vertexCount; vertexIndex++) {
		Vec2f v = verts[vertexIndex];
		F32 p = Dot(v, normal);
		if (p > firstDistance) {
			firstDistance = p;
			firstIndex = vertexIndex;
		}
	}

	S32 secondIndex = -1;
	F32 secondDistance = 0;
	for (U32 vertexIndex = 0; vertexIndex < vertexCount; vertexIndex++) {
		if (vertexIndex != firstIndex) {
			Vec2f v = verts[vertexIndex];
			F32 p = Dot(v, normal);
			if (secondIndex == -1 || p > secondDistance) {
				secondDistance = p;
				secondIndex = vertexIndex;
			}
		}
	}
	Assert(secondIndex > -1);
	Face result = {};
	result.index = firstIndex;
	result.points[0] = verts[firstIndex];
	result.points[1] = verts[secondIndex];
	return (result);
}

internal U32 PlaneEdgeContactGenerator(Physics *world, Transform transformA, Transform transformB, Shape *shapeA, Shape *shapeB, U32 offset, Contact *contacts) {
	U32 result = 0;

	Vec2f normal = transformA.rot.col1;
	Vec2f posA = transformA.pos;
	EdgeShape *edgeB = GetEdgeShape(shapeB);
	Vec2f* localVertsB = edgeB->localVerts;
	U32 numVerticesB = edgeB->numVertices;

	Mat2f transposeB = Transpose2(transformB.rot);
	Vec2f normalB = -(normal * transposeB);

	Vec2f pointsOnA[2];
	F32 distanceToA[2];
	U32 clipTypeOnB[2];
	U32 numPoints = 0;

	Face faceB = GetFace(numVerticesB, localVertsB, normalB);
	Vec2f supportPointB = faceB.points[0] * transformB;
	Vec2f planePoint = posA;
	Vec2f distanceToPlane = supportPointB - planePoint;
	F32 d = Dot(distanceToPlane, normal);
	if (d > 0) {
		pointsOnA[0] = supportPointB + normal * -d;
		distanceToA[0] = d;
		numPoints = 1;
		clipTypeOnB[0] = 0;
	} else {
		Vec2f p0 = faceB.points[0] * transformB;
		Vec2f p1 = faceB.points[1] * transformB;
		Vec2f distance0 = p0 - planePoint;
		F32 d0 = Dot(distance0, normal);
		if (d0 <= 0) {
			U32 idx = numPoints++;
			pointsOnA[idx] = p0 + normal * -d0;
			distanceToA[idx] = d0;
			clipTypeOnB[idx] = 1;
		}

		Vec2f distance1 = p1 - planePoint;
		F32 d1 = Dot(distance1, normal);
		if (d0 <= 0) {
			U32 idx = numPoints++;
			pointsOnA[idx] = p1 + normal * -d1;
			distanceToA[idx] = d1;
			clipTypeOnB[idx] = 2;
		}
	}

	for (U32 pointIndex = 0; pointIndex < numPoints; pointIndex++) {
		// NOTE(final): Plane has just one face - therefore only one feature
		U32 feature = 0 * CONTACT_FEATURE_FACEA_PRIME + (faceB.index + 1) * CONTACT_FEATURE_FACEB_PRIME + clipTypeOnB[pointIndex] * CONTACT_FEATURE_CLIP_PRIME;
		Contact *contact = contacts + (offset + result++);
		SetContact(contact, distanceToA[pointIndex], normal, pointsOnA[pointIndex] - posA, feature);
	}

	return(result);
}
Initializing the contact with mass ratio and build contact arms:
 
// Get center of masses
Vec2f comA = bodyA->position + bodyA->centroid;
Vec2f comB = bodyB->position + bodyB->centroid;

// Get contact points in world space (Point is already rotated, B is calculated on the fly)
Vec2f contactWorldA = bodyA->position + shapeA->localTransform.pos + contact->point;
Vec2f contactWorldB = contactWorldA + contact->normal * contact->distance;

// Get contact arms
contact->rA = contactWorldA - comA;
contact->rB = contactWorldB - comB;

// Get normal rotation on the z-plane
F32 rnA = Cross(contact->rA, contact->normal);
F32 rnB = Cross(contact->rB, contact->normal);

// Build normal matrix
F32 normalTensor = bodyA->invMass + bodyB->invMass + bodyA->invInertia * rnA * rnA + bodyB->invInertia * rnB * rnB;
contact->normalMass = normalTensor > 0 ? 1.0f / normalTensor : 0;

Solving the Impulse for each contact:
 
// Relative velocity on contact = vB + cross(wB, rB) - vA - cross(wA, rA);
Vec2f vAB = *vB + Cross(bodyB->angularVelocity, contact->rB) - *vA - Cross(bodyA->angularVelocity, contact->rA);

// Project relative velocity on normal
F32 projNormalVel = Dot(vAB, normal);

// Normal impulse
F32 remove = projNormalVel + contact->distance / dt;
F32 normalImpulse = remove * contact->normalMass;
F32 newImpulse = Min(contact->normalImpulse + normalImpulse, 0.0f);
F32 impulseChange = newImpulse - contact->normalImpulse;
contact->normalImpulse = newImpulse;

// Apply Impulses
Vec2f impulseA = normal * impulseChange;
bodyA->velocity += impulseA * bodyA->invMass;
bodyA->angularVelocity += Cross(contact->rA, impulseA) * bodyA->invInertia;
Vec2f impulseB = -normal * impulseChange;
bodyB->velocity += impulseB * bodyB->invMass;
bodyB->angularVelocity += Cross(contact->rB, impulseB) * bodyB->invInertia;
Edited by Finalspace

Share this post


Link to post
Share on other sites

Screw it, i rewrite the entire physics system. Seems everything is wrong, but wrong in a way that it never explodes but adds tiny little movement/rotation to the system so that nothing gets to a rest state. So sad. It would be so great when i can do that kind of programming at full time, but always pausing and getting back after weeks/months just kills it...

Edited by Finalspace

Share this post


Link to post
Share on other sites

I really recommend taking Box2D Lite as either a reference or by porting it to 3D. In particular if you are on and off on the project as it provides some kind of red line to stick to, Once you have the basic architecture you can add more functionality. There is also a lot of information from the GDC math and physics tutorial you can access here:

 

http://essentialmath.com/tutorial.htm

http://box2d.org/downloads/

 

Speculative contacts really create more problems than they solve from my experience. Continuous Collision Detection (CCD) and Continuous Physics (CP) are among the most difficult problems I can think of and there doesn't exist a simple *general* solution. Box2D has a pretty good CCD and CP solution. I would try to understand sub-stepping and its problems first before using Speculative Contacts.

Edited by Dirk Gregorius

Share this post


Link to post
Share on other sites

I really recommend taking Box2D Lite as either a reference or by porting it to 3D. In particular if you are on and off on the project as it provides some kind of red line to stick to, Once you have the basic architecture you can add more functionality. There is also a lot of information from the GDC math and physics tutorial you can access here:
 
http://essentialmath.com/tutorial.htm
http://box2d.org/downloads/
 
Speculative contacts really create more problems than they solve from my experience. Continuous Collision Detection (CCD) and Continuous Physics (CP) are among the most difficult problems I can think of and there doesn't exist a simple *general* solution. Box2D has a pretty good CCD and CP solution. I would try to understand sub-stepping and its problems first before using Speculative Contacts.

 
Oh i have done this a lot of times, ported box2d lite to other languages like javascript, java, c, pascal, c# mostly with success: The basic idea of box2d lite is adding the position correction directly to the impulses and using transpose() to get a reverse rotation for bringing the normal of B into the space A - so that the support function works with non-transformed vertices. Also it simplifies a lot by just using box vs box and uses simple sat for getting the normal/distances, but it has problems: Some weird stability issues in some cases - which i always have :(, especially when different masses are used, no restitution support (But can be easily added but with the cost of stability), Non moving penetrated / teleported bodies will not be corrected.
 
Also most of that math and physics papers / presentations i have already read more than once. I know all the theory behind, what a reference / incident face is, finding support points, separating axis theorem + support sat, planes, how clipping works, how to bring a vector into the space of another, how rotations work, impulse equations, interpolation, vector math, closest points, contact generation, collision detection, you name it. But it seems that i am always implementing systems which comes with more downsides than others, so i mostly ending up fighting with them.
 

Speculative contacts really create more problems than they solve from my experience. Continuous Collision Detection (CCD) and Continuous Physics (CP) are among the most difficult problems I can think of and there doesn't exist a simple *general* solution. Box2D has a pretty good CCD and CP solution. I would try to understand sub-stepping and its problems first before using Speculative Contacts.


I disagree, speculative contacts in general are very easy to understand. I have implemented it without rotation dynamics a ton of times and always got a very stable simulation out-of-it. So its not a bad technique - the stability compared to other impulse systems is incredible, but it has its downsides: Hard to get rotation dynamics working, internal edges, ghost collisions, no restitution. Edited by Finalspace

Share this post


Link to post
Share on other sites

I disagree, speculative contacts in general are very easy to understand...


It's probably a problem if you are disagreeing with Dirk. Not to be rude, but if someone with a lot of professional experience comes into a forum and grants you free advice it's best to take it with a large dose of humility. The advice to look at Box2D Lite is where you probably need to return to. Since you're having trouble implementing something based on the wildbunny link, you definitely did not understand much behind Box2D Lite, regardless of how many names of algorithms you ported or can recite.
 

The basic idea of box2d lite is adding the position correction directly to the impulses and using transpose() to get a reverse rotation for bringing the normal of B into the space A - so that the support function works with non-transformed vertices. Also it simplifies a lot by just using box vs box and uses simple sat for getting the normal/distances, but it has problems: Some weird stability issues in some cases - which i always have :(, especially when different masses are used, no restitution support (But can be easily added but with the cost of stability), Non moving penetrated / teleported bodies will not be corrected.


You described Baumgarte and change of basis transformation. That doesn't mean you know why Baumgarte works the way it does, or that you know why the impulse solver converges at all. If anything it shows you only carry a basic surface level understanding unsuited for creating a working simulation. Box2D Lite doesn't have "stability issues", at least not like you are describing. These stability problems are just bugs in your code because you didn't fully understand what was happening, and it sounds like you gave up and never figured it out. You can't just blame the systems and say they are flawed, or they have weaknesses, when every state of the art physics engine employs, in essence, the same techniques Box2D Lite presents. It's incredibly arrogant to blame tried and vetted algorithms for your own bugs and lack of understanding, and worse to come in and post in an adversarial manner towards very helpful and experienced devs like Dirk.

This kind of attitude is exactly why posters like Dirk are incredibly rare, and they eventually stop posting forever at some point.

Anyway, I have some professional experience with speculative contacts in a major physics engine, and I can say the experience sucked. There were tons of bugs and any performance gains were completely swamped by false positives, collision detection bugs, and required expensive content work-arounds or very hacky solutions. All in all, I would never use another speculative contact scheme for anything after such an experience. Speculative contacts seems extremely complicated and don't really offer much a benefit to justify the complexity. Edited by Randy Gaul

Share this post


Link to post
Share on other sites
Before diving into CCD and CP 
 
Read the maths and physics sections in the book Essential Math for Games and Interactive Applications
Study all GDC presentations by Erin and Dirk
Implement a decent debug drawing interface
Port Box2D Lite to 3D
Learn and implement the GJK and SAT algorithms
Study sub-stepping
Read Brian Mirtich CCD and CP thesis
Study root-finding methods, it helps to understand the general problem of TOI computation 
 
Make sure you understand every line of code of your physics engine and the limits of the engine. Comment the important parts, specially why and not how you wrote something. If you have a feature that you don't know how or why it worked (maybe you was too anxious to see the results on the screen or didn't have time to read about it properly) then remove it from your project and add it only after understanding. Usually you can see a well architected physics engine when the author understands the problems.  
 
You can learn a lot when creating something from scratch using only the literature and skipping existing implementation, but currently this is impractical and error prone since there is a great deal of intuitive resources and discussions out there (e.g. gamedev/Bullet forums).
 
Box2D Lite doesn't have a state of the art dynamics for games implemented because of its position corrector. Once you're familiar with it you can try porting *Box2D* position correctors to 3D. 
Edited by Irlan Robson

Share this post


Link to post
Share on other sites

First of all, i am really sorry! I never would dream about being rude to any professionals here.
I am really grateful to all the resources and explainations you guys provide. I learned so much from that: Thanks for that.
It just seems i am really frustrated right now and was going the wrong direction since i started doing physics programming.

So i get the point: drop the idea about speculative contacts entirely and start implementing sequential impulse based on box2d lite and build up from there.
Now if i think about that, i always had issues with speculative contacts in the past years and switching to sequential impulse had always got me a working solution.
Why was i so focused on that technique for so long... maybe because its not bad for platformers which is my target genre after all and i liked the physics in little big planet a lot.

And yes i dont and will never understand how the math behind the impulses really works, i never learned any advanced math at all.
So differential equations, linear systems, linear algebra, constraints will ever be a blackbox for me.
But i somehow get results when you look at things i have made (xenorate media player, fluid sandbox, tons of tech demos and physics stuff) - so i am not totally a lost sheep the dark.

Share this post


Link to post
Share on other sites
Well i fixed it - now it works perfectly with plane vs polygon.
The contact points generated was rotated in the wrong direction and i was storing them in world space.

This is what i usually do:

vec2 normalA = transformA.rotation.col1; // normal from plane A
mat2 rotAToB = transpose(transformB.rotation);
vec2 normalB = normalA * rotAToB;
vec2 supportPoint1 = GetSupportPoint(-normalB, localVertsB, numVertsB) * transformB.rotation + transformB.translation;

Using opengl mathematics library i only get proper results this way:

vec2 normalA = transformA.rotation.col1; // normal from plane A
vec2 normalB = normalA * transformB.rotation; // No transpose here
vec2 supportPoint1 = GetSupportPoint(-normalB, localVertsB, numVertsB) * glm::transpose(transformB.rotation) + transformB.translation; // But here will switch from column to row order

For rendering i also use glm::transpose(transform.rotation) + transform.translation.
But it works, the box falls to the plane and does not jump/penetrate at all. Edited by Finalspace

Share this post


Link to post
Share on other sites

Seems like i fixed the "rewrite" - polygon vs polygon now works and is kind of stable. But i havent adressed my initial issue yet i have asked here for help in the first place.

But like i have sayed in my apology i will drop speculative contacts entirely and implement the system used in box2d lite.

 

Anyway if someone want to play with what i made, here are the full source (one cpp file -> win32, glm as its only dependency):

http://pastebin.com/n4y5uZw1

Edited by Finalspace

Share this post


Link to post
Share on other sites
So after hours of debugging i found the culprits. It was just 3 issues boiled down to classic porting bugs.
I ported it over from my own java physics system i have written a year ago - but these do not had either rotation dynamics nor operator overloading in the math.

Long Story short: Now everything works perfectly with a sequential impulse system.
I also extended it to support multiple shapes per body and improved the broadphase.
The only thing i dont not like is this:
 
internal AABB ComputePlaneAABB(PlaneShape *plane, const Transform &transform, F32 tolerance) {
	AABB result = {};

	Vec2f e = V2(tolerance, tolerance);
	Vec2f center = transform.pos;
	Vec2f normal = transform.rot.col1;
	Vec2f tangent = Cross(normal, 1.0f);

	Vec2f startPoint = center + tangent * plane->len * 0.5f;
	Vec2f endPoint = center + tangent * -plane->len * 0.5f;

	result.min = V2Min(startPoint, endPoint);
	result.max = V2Max(startPoint, endPoint);
	result.min -= e;
	result.max += e;

	// TODO(final): What have i done here??? (ported from old java source - need to revist that soon)
	Vec2f depth = V2(tolerance * 10.0f, tolerance * 10.0f);
	if (normal.x > 0.99f || normal.y > 0.99f) {
		result.min += normal * -Dot(normal, depth);
	} else if (normal.x < 0.99f || normal.y < 0.99f) {
		result.max += normal * Dot(normal, depth);
	}

	return(result);
}

internal void ComputeBodyAABB(Body *body, const Transform &prevWorldTransform, const Transform &nextWorldTransform, F32 tolerance) {
	// TODO(final): This must be optimized!
	body->aabb.min = V2(F32_MAX, F32_MAX);
	body->aabb.max = V2(-F32_MAX, -F32_MAX);
	for (U32 shapeIndex = 0; shapeIndex < body->numShapes; ++shapeIndex) {
		Shape *shape = body->shapes + shapeIndex;
		Transform prevTransform = shape->localTransform * prevWorldTransform;
		Transform nextTransform = shape->localTransform * nextWorldTransform;
		AABB prevShapeAABB;
		AABB nextShapeAABB;
		switch (shape->type) {
			case ShapeType::ShapeType_Box:
				prevShapeAABB = ComputeEdgeAABB((EdgeShape *)&shape->box, prevTransform, tolerance);
				nextShapeAABB = ComputeEdgeAABB((EdgeShape *)&shape->box, nextTransform, tolerance);
				break;
			case ShapeType::ShapeType_Polygon:
				prevShapeAABB = ComputeEdgeAABB((EdgeShape *)&shape->polygon, prevTransform, tolerance);
				nextShapeAABB = ComputeEdgeAABB((EdgeShape *)&shape->polygon, nextTransform, tolerance);
				break;
			case ShapeType::ShapeType_LineSegment:
				prevShapeAABB = ComputeEdgeAABB((EdgeShape *)&shape->lineSegment, prevTransform, tolerance);
				nextShapeAABB = ComputeEdgeAABB((EdgeShape *)&shape->lineSegment, nextTransform, tolerance);
				break;
			case ShapeType::ShapeType_Circle:
				prevShapeAABB = ComputeCircleAABB(&shape->circle, prevTransform, tolerance);
				nextShapeAABB = ComputeCircleAABB(&shape->circle, nextTransform, tolerance);
				break;
			case ShapeType::ShapeType_Plane:
				prevShapeAABB = ComputePlaneAABB(&shape->plane, prevTransform, tolerance);
				nextShapeAABB = ComputePlaneAABB(&shape->plane, nextTransform, tolerance);
				break;
			default:
				Assert(!"Invalid shape type for aabb computation!");
				break;
		}
		AABB combinedShapeAABB = CombineAABB(prevShapeAABB, nextShapeAABB);
		body->aabb.min = V2Min(body->aabb.min, combinedShapeAABB.min);
		body->aabb.max = V2Max(body->aabb.max, combinedShapeAABB.max);
	}
}

If someone is interested i can modify the pastebin source for using sequential impulse.

Anyway, i lastly want to give my thanks to dirk and randy for the help - even when i dont deserve it. I know my place now. Edited by Finalspace

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this