Havok Physics Engine Tutorial Series: Chain

In this tutorial, I will show you how to create a rigid body chain using the simple distance constraint which is wrapped into the hkpStiffSpringConstraint object in the Havok Physics SDK. This tutorial will be building on top of to the Simple Distance Constraint Tutorial and Picking Tutorial that we did earlier.

Creating a Chain using a Simple Distance Constraint
In the previous tutorial, we saw how we could use a distance constraint to limit a pair of rigid bodies. In this tutorial, I will show you how to use a pair of distance constraints to create a chain of rigid bodies.You can pick any rigid body and move it around using the mouse.  So lets get started.

The AddRigidBodies function
This function is defined as follows

	void AddRigidBodies() {
	//add the falling box
	{
	hkVector4 halfExtents(0.5f, 0.5f, 0.125f);
	hkpBoxShape* boxShape = new hkpBoxShape(halfExtents);
	boxShape->setRadius(0.001f);
	hkpRigidBodyCinfo ci;
	ci.m_shape = boxShape;
	ci.m_motionType = hkpMotion::MOTION_DYNAMIC;
	const hkReal boxMass(10.0f);
	hkMassProperties massProps;
	hkpInertiaTensorComputer::computeShapeVolumeMassProperties(boxShape, boxMass, massProps);
	ci.setMassProperties(massProps);
	for(int i=0;i<MAX_BOXES;++i) {
	if(i==(MAX_BOXES-1))
	ci.m_motionType = hkpMotion::MOTION_FIXED;
	ci.m_position = hkVector4(0.0f,1.0f+i*1.125f,0.0f);
	hkpRigidBody* rigidBody = new hkpRigidBody(ci);
	boxes.push_back(static_cast<hkpRigidBody*>(g_pWorld->addEntity(rigidBody)));
	}
	boxShape->removeReference();
	//now add springs between adjacent boxes
	hkpStiffSpringConstraintData* spring = new hkpStiffSpringConstraintData();
	hkpStiffSpringConstraintData* spring2 = new hkpStiffSpringConstraintData();
	for(size_t i=0;i<boxes.size()-1;i++) {
	hkTransform b1 = boxes[i]->getTransform();
	hkTransform b2 = boxes[i+1]->getTransform();
	hkVector4f t1 = b1.getTranslation();
	hkVector4f t2 = b2.getTranslation();
	t1.add(hkVector4f(-0.5,0.5,0));
	t2.add(hkVector4f(-0.5,-0.5,0));
	spring->setInWorldSpace(b1,b2, t1, t2);
	{
	hkpConstraintInstance* constraint = new hkpConstraintInstance(boxes[i], boxes[i+1], spring );
	g_pWorld->addConstraint(constraint);
	constraint->removeReference();
	}
	t1.add(hkVector4f(1,0,0));
	t2.add(hkVector4f(1,0,0));
	spring2->setInWorldSpace(b1,b2, t1, t2);
	{
	hkpConstraintInstance* constraint = new hkpConstraintInstance(boxes[i], boxes[i+1], spring2 );
	g_pWorld->addConstraint(constraint);
	constraint->removeReference();
	}
	}
	spring->removeReference();
	spring2->removeReference();
	}
	//create the static box where the smaller box will fall
	{
	hkVector4 halfExtents(20.0f, 2.0f, 20.f);
	hkpBoxShape* boxShape = new hkpBoxShape(halfExtents);
	hkpRigidBodyCinfo ci;
	ci.m_shape = boxShape;
	ci.m_position = hkVector4(0, -2, 0);
	ci.m_motionType = hkpMotion::MOTION_FIXED;
	boxShape->setRadius(0.001f);
	hkpRigidBody* rigidBody = new hkpRigidBody(ci);
	boxShape->removeReference();
	g_pWorld->addEntity(rigidBody)->removeReference();
	}
	}

view raw Chain_01 hosted with ❤ by GitHub

Lets have a look at this function piece by piece. The AddRigidBodies function first creates a number of blocks. This is similar to how we created the boxes in the Multiple Bouncing Box Tutorial.

	hkVector4 halfExtents(0.5f, 0.5f, 0.125f);
	hkpBoxShape* boxShape = new hkpBoxShape(halfExtents);
	boxShape->setRadius(0.001f);
	hkpRigidBodyCinfo ci;
	ci.m_shape = boxShape;
	ci.m_motionType = hkpMotion::MOTION_DYNAMIC;
	const hkReal boxMass(10.0f);
	hkMassProperties massProps;
	hkpInertiaTensorComputer::computeShapeVolumeMassProperties(boxShape, boxMass, massProps);
	ci.setMassProperties(massProps);
	for(int i=0;i<MAX_BOXES;++i) {
	if(i==(MAX_BOXES-1))
	ci.m_motionType = hkpMotion::MOTION_FIXED;
	ci.m_position = hkVector4(0.0f,1.0f+i*1.125f,0.0f);
	hkpRigidBody* rigidBody = new hkpRigidBody(ci);
	boxes.push_back(static_cast<hkpRigidBody*>(g_pWorld->addEntity(rigidBody)));
	}
	boxShape->removeReference();

view raw Chain_02 hosted with ❤ by GitHub

Next, we create another loop but this time we create a pair of stiff spring constraints. The constraints are placed in parallel between each pair of boxes. The offsets are calculated based on the size of the box.

	//now add springs between adjacent boxes
	hkpStiffSpringConstraintData* spring = new hkpStiffSpringConstraintData();
	hkpStiffSpringConstraintData* spring2 = new hkpStiffSpringConstraintData();
	for(size_t i=0;i<boxes.size()-1;i++) {
	hkTransform b1 = boxes[i]->getTransform();
	hkTransform b2 = boxes[i+1]->getTransform();
	hkVector4f t1 = b1.getTranslation();
	hkVector4f t2 = b2.getTranslation();
	t1.add(hkVector4f(-0.5,0.5,0));
	t2.add(hkVector4f(-0.5,-0.5,0));
	spring->setInWorldSpace(b1,b2, t1, t2);
	{
	hkpConstraintInstance* constraint = new hkpConstraintInstance(boxes[i], boxes[i+1], spring );
	g_pWorld->addConstraint(constraint);
	constraint->removeReference();
	}
	t1.add(hkVector4f(1,0,0));
	t2.add(hkVector4f(1,0,0));
	spring2->setInWorldSpace(b1,b2, t1, t2);
	{
	hkpConstraintInstance* constraint = new hkpConstraintInstance(boxes[i], boxes[i+1], spring2 );
	g_pWorld->addConstraint(constraint);
	constraint->removeReference();
	}
	}
	spring->removeReference();
	spring2->removeReference();
	}

view raw Chain_03 hosted with ❤ by GitHub

Finally, we create the ground rigid body as was seen in earlier tutorials.

	//create the static box where the smaller box will fall
	{
	hkVector4 halfExtents(20.0f, 2.0f, 20.f);
	hkpBoxShape* boxShape = new hkpBoxShape(halfExtents);
	hkpRigidBodyCinfo ci;
	ci.m_shape = boxShape;
	ci.m_position = hkVector4(0, -2, 0);
	ci.m_motionType = hkpMotion::MOTION_FIXED;
	boxShape->setRadius(0.001f);
	hkpRigidBody* rigidBody = new hkpRigidBody(ci);
	boxShape->removeReference();
	g_pWorld->addEntity(rigidBody)->removeReference();
	}

view raw Chain_04 hosted with ❤ by GitHub

The picking and rigid body manipulation functions are virtually unchanged. I simply copied the code from the picking tutorial and it worked without a problem.

Output
That's all, after compiling and building the given code, you will see a chain of boxes. The top most box is fixed to ensure that the chain does not fall down. The rest of the boxes can be picked. You can pick the box by left clicking and dragging the mouse to reposition it as shown in the following figure.

You can get the full source code from my github repo https://github.com/mmmovania/HavokPhysicsTutorials

Controls:
Left click to rotate, left click on box to pick and reposition
Middle click to zoom
Right click to pan

What's next:
In the next tutorial, I will see if I can do a simple cloth using the simple distance constraint object.

Havok Physics Engine Tutorial Series: Simple Distance Joint

In this tutorial, I will show you how to render a simple distance joint (a basic stiff spring). We will add to the picking code we covered in the last tutorial so that we can displace a given box using the mouse.

The scene in this tutorial contains two boxes: a static box and a dynamic box. The dynamic box and the static box are linked with a distance constraint. We can move the dynamic box using mouse and due to the spring constraint, it maintains a certain distance from the static box. So lets get started. For this demo, the AddRigidBodies function is as follow.

	void AddRigidBodies() {
	//add the falling movingBox
	{
	hkVector4 halfExtents(0.5f, 0.5f, 0.5f);
	hkpBoxShape* boxShape = new hkpBoxShape(halfExtents);
	hkpRigidBodyCinfo ci;
	ci.m_shape = boxShape;
	ci.m_position = pos2;
	ci.m_motionType = hkpMotion::MOTION_DYNAMIC;
	boxShape->setRadius(0.001f);
	const hkReal boxMass(10.0f);
	hkMassProperties massProps;
	hkpInertiaTensorComputer::computeShapeVolumeMassProperties(boxShape, boxMass, massProps);
	ci.setMassProperties(massProps);
	hkpRigidBody* rigidBody = new hkpRigidBody(ci);
	movingBox = static_cast<hkpRigidBody*>(g_pWorld->addEntity(rigidBody));
	rigidBody->removeReference();
	//create the fixed box
	ci.m_position = pos1;
	ci.m_motionType = hkpMotion::MOTION_FIXED;
	rigidBody = new hkpRigidBody(ci);
	fixedBox = static_cast<hkpRigidBody*>(g_pWorld->addEntity(rigidBody));
	rigidBody->removeReference();
	boxShape->removeReference();
	}
	//create the static box where the smaller movingBox will fall
	{
	hkVector4 halfExtents(20.0f, 2.0f, 20.f);
	hkpBoxShape* boxShape = new hkpBoxShape(halfExtents);
	hkpRigidBodyCinfo ci;
	ci.m_shape = boxShape;
	ci.m_position = hkVector4(0, -2, 0);
	ci.m_motionType = hkpMotion::MOTION_FIXED;
	boxShape->setRadius(0.001f);
	hkpRigidBody* rigidBody = new hkpRigidBody(ci);
	boxShape->removeReference();
	g_pWorld->addEntity(rigidBody)->removeReference();
	}
	//
	// CREATE STIFF SPRING CONSTRAINT
	//
	{
	hkpStiffSpringConstraintData* spring = new hkpStiffSpringConstraintData();
	// Create constraint
	spring->setInWorldSpace(movingBox->getTransform(), fixedBox->getTransform(), pos2, pos1);
	//
	// Create and add the constraint
	//
	{
	hkpConstraintInstance* constraint = new hkpConstraintInstance(movingBox, fixedBox, spring );
	g_pWorld->addConstraint(constraint);
	constraint->removeReference();
	}
	spring->removeReference();
	}
	}

view raw Distance_01 hosted with ❤ by GitHub

Adding a spring constraint
The only difference here is the addition of the stiff spring constraint. To create it, we first create the hkpStiffSpringConstraintData object. We then pass it the world space position of where the spring is located along with the two rigid bodies between which this constraint is created. Next,
the hkpStiffSpringConstraintData object and the two rigid bodies are passed to the hkpConstraintInstance object and then hkpWorld::addConstraint function is called passing it the hkpConstraintInstance object as shown in the following code snippet.

	hkpStiffSpringConstraintData* spring = new hkpStiffSpringConstraintData();
	// Create constraint
	spring->setInWorldSpace(movingBox->getTransform(), fixedBox->getTransform(), pos2, pos1);
	// Create and add the constraint
	{
	hkpConstraintInstance* constraint = new hkpConstraintInstance(movingBox, fixedBox, spring );
	g_pWorld->addConstraint(constraint);
	constraint->removeReference();
	}
	spring->removeReference();

view raw Distance_02 hosted with ❤ by GitHub

Rendering of the spring constraint
Another change we did is in the render function. It is now changed to render the spring constraint as well along with the two rigid bodies. The picked dynamic rigid body is rendered green when picked.

	//draw spring
	glBegin(GL_LINES);
	glColor3f(0,1,0);
	glVertex3f(pos1.getComponent(0), pos1.getComponent(1), pos1.getComponent(2));
	glVertex3f(pos2.getComponent(0), pos2.getComponent(1), pos2.getComponent(2));
	glEnd();
	//draw boxes
	glEnable(GL_LIGHTING);
	glColor3f(1,1,1);
	DrawBox(fixedBox);
	if(wasHit)
	glColor3f(0,1,0);
	DrawBox(movingBox, true);
	glDisable(GL_LIGHTING);

view raw Distance_03 hosted with ❤ by GitHub

Output
That's all, after compiling and running, you will see two boxes linked with a spring constraint. You can pick the lower (dynamic) box by left clicking and dragging the mouse to reposition it as shown in the following figure.

After displacement, the spring constraint acts on the dynamic rigid body to ensure that it remains at the fixed distance that was given at the time of instantiation of the constraint.

You can get the full source code from my github repo https://github.com/mmmovania/HavokPhysicsTutorials

Controls:
Left click to rotate, left click on box to pick and reposition
Middle click to zoom
Right click to pan

What's next:
In the next tutorial, I will show you how to create a simple chain.

Havok Physics Engine Tutorial Series: Picking

In this tutorial, I will show you how to pick a rigid body using picking features available in the Havok Physics SDK. This tutorial will be adding to the SimpleBox tutorial that we did earlier. We will augment it with picking using the mouse.

Picking can be implemented in numerous ways. We will cast a ray from the clicked position on the screen to the 3D world using the hkpWorld::castRay function provided by the Havok Physics SDK. The user clicks on the screen at a position(x,y). We implement this in the PickActor function. The first thing we need is to determine the ray that has its start position at the clicked position(x,y)on screen and an end position that is at the far clip plane.

Determining the ray from the clicked position
We first determine the z position of the clicked screen space point. We know from the projection transformation that a given point (x,y), has z=0 at the near clip plane and z=1 at the far clip plane. We can use this knowledge to create two 3D points from the given screen space point, one at near clip plane (x,y,0) and another at the far clip plane (x,y,1). We then unproject these two 3D points to obtain the world space points for the two 3D points. These are the ray start and the ray end points that we will pass to the hkpWorld::castRay function.

	void PickActor(int x, int y) {
	hkVector4f rayStart(0,0,0,0);
	hkVector4f rayEnd(0,0,0,0);
	ViewUnProject(x, y, 0.0f, rayStart);
	ViewUnProject(x, y, 1.0f, rayEnd);

view raw Picking_01 hosted with ❤ by GitHub

The ViewUnProject is a utility function that we define as follows. Note that this function definition will be modified if you use DirectX or any other graphics library.

	void ViewUnProject(int xi, int yi, float depth, hkVector4 &v)
	{
	yi = viewport[3] - yi - 1;
	GLdouble wx=0, wy=0, wz=0;
	gluUnProject((GLdouble) xi, (GLdouble) yi, (GLdouble) depth,
	MV, P, viewport, &wx, &wy, &wz);
	v.set( wx, wy, wz, 1);
	}

view raw Picking_02 hosted with ❤ by GitHub

This function uses the current viewport to determine the screen size. It then passes the given x,y and z position, the current modelview and projection matrices and the current viewport to the gluUnProject function from the OpenGL utility library. Of course if you are using any other library, you will have to find the appropriate function or create your own function. After the unproject call, we get the world space point for the given point.

At the end of these calls, we have the ray start and ray end points.

Casting picking ray using hkpWorld::castRay function
OK now we have our ray, we can call g_pWorld->castRay function. If we look at the Havok SDK reference, we can see that the castRay function takes two parameters. hkpWorldRaycastInput and hkpClosestRayHitCollector. You will have to include the following new headers for these

	#include <Physics2012\Collide\Query\CastUtil\hkpWorldRayCastInput.h>
	#include <Physics2012\Collide\Query\CastUtil\hkpWorldRayCastOutput.h>
	#include <Physics2012\Collide/Query/Multithreaded/RayCastQuery/hkpRayCastQueryJobs.h>
	//for raycasthitcollector
	#include <Physics2012/Collide/Query/Collector/RayCollector/hkpClosestRayHitCollector.h>

view raw Picking_03 hosted with ❤ by GitHub

The hkpRaycastInput object has three fields, m_from which is the ray start point, m_to which is the ray end point, m_filterInfo which we pass 0. We then lock the hkpWorld, castRay and then unlock hkpWorld as shown below.

	hkpWorldRayCastInput input;
	hkpClosestRayHitCollector result;
	input.m_from = rayStart;
	input.m_to = rayEnd;
	input.m_filterInfo = 0;
	g_pWorld->lock();
	g_pWorld->castRay(input, result);
	g_pWorld->unlock();

view raw Picking_04 hosted with ❤ by GitHub

If any rigidbody in the physics world intersects with the given ray, the hkpClosestRayHitCollector::hasHit() function returns true. We call this function to determine the possible hit condition. Calling hkpClosestRayHitCollector::getHit() function returns the hkpWorldRaycastOutput object. This object contains the m_rootCollidable field which is the intersected rigid body. To get a rigidbody pointer from the m_rootCollidable object, we call hkpGetRigidBody function. Now since our world contains both static and dynamic rigid bodies, we are only interested in the dynamic rigid bodies so we determine the motion type field of the returned rigid body.

After we are sure that we have a dynamic rigid body, we then find the intersected point on the rigid body and then project the point to get the 3D hit point. This is just to render the intersected point. Next, we set the motion state of the picked rigid body as MOTION_FIXED to make the picked rigid body a static rigid body. This is to ensure that when we move the rigidbody around, it does not undergo motion. 

	if(result.hasHit()) {
	hkpRigidBody* pRB = hkpGetRigidBody(result.getHit().m_rootCollidable);
	wasHit = pRB->getMotionType() != hkpMotion::MOTION_FIXED;
	if(wasHit) {
	gSelectedActor = pRB;
	hkVector4f pos2 = gSelectedActor->getTransform().getTranslation();
	intersectionPointWorld.setInterpolate4( input.m_from, input.m_to, result.getHit().m_hitFraction );
	int hitx=0, hity=0;
	ViewProject(intersectionPointWorld, hitx, hity, hitz);
	g_pWorld->lock();
	gSelectedActor->setMotionType(hkpMotion::MOTION_FIXED);
	g_pWorld->unlock();
	}
	}

view raw Picking_05 hosted with ❤ by GitHub

In the mouse move event handler, we determine if the user has picked a rigid body. If there is a picked rigid body, we set the translation field of the rigid body's transform to the unprojected mouse point as was shown earlier.

	void MoveActor(int x, int y) {
	if (!wasHit)
	return;
	hkVector4 pos;
	ViewUnProject(x,y, hitz, pos);
	g_pWorld->lock();
	hkTransform trans = gSelectedActor->getTransform();
	trans.setTranslation(pos);
	gSelectedActor->setTransform(trans);
	g_pWorld->unlock();
	}

view raw Picking_06 hosted with ❤ by GitHub

Finally, when the mouse button up event is raised, we set the motion state of the picked rigid body to MOTION_DYNAMIC.to ensure that the picked rigid body is physically simulated again.

	if (s == GLUT_UP ) {
	if(wasHit) {
	g_pWorld->lock();
	gSelectedActor->setMotionType(hkpMotion::MOTION_DYNAMIC);
	g_pWorld->unlock();
	}
	wasHit = false;
	state=-1;
	}

view raw Picking_07 hosted with ❤ by GitHub

Output
That's all, you will see a box falling due to gravity. You can pick the box by left clicking and dragging the mouse to reposition it as shown in the following figure.

You can get the full source code from my github repo https://github.com/mmmovania/HavokPhysicsTutorials

Controls:
Left click to rotate, left click on box to pick and reposition
Middle click to zoom
Right click to pan

What's next:
In the next tutorial, I will show you how to create a simple distance joint.

Havok Physics Engine Tutorial Series: Multiple Bouncing Boxes

OK Now that we know hot to create one box. We will now add multiple boxes and let them fall under gravity and collide with each other.  The InitializeHavok functions remains the same. The changes here are in the AddRigidBodies function, the ShutdownHavok function and the OnRender function. Lets go through each of these one by one.

We create a global vector to store all of our rigid bodies (boxes).

	const int MAX_BOXES = 100;
	std::vector<hkpRigidBody*> boxes;

view raw MultipleBoxes_01 hosted with ❤ by GitHub

The AddRigidBodies function
The AddRigidBodies function is changed to the following.

	void AddRigidBodies() {
	//add the falling box
	{
	hkVector4 halfExtents(0.5f, 0.5f, 0.5f);
	hkpBoxShape* boxShape = new hkpBoxShape(halfExtents);
	boxShape->setRadius(0.001f);
	hkpRigidBodyCinfo ci;
	ci.m_shape = boxShape;
	ci.m_motionType = hkpMotion::MOTION_DYNAMIC;
	const hkReal boxMass(10.0f);
	hkMassProperties massProps;
	hkpInertiaTensorComputer::computeShapeVolumeMassProperties(boxShape, boxMass, massProps);
	ci.setMassProperties(massProps);
	for(int i=0;i<MAX_BOXES;++i) {
	ci.m_position = hkVector4(0.0f,5.0f+5*i,0.0f);
	hkpRigidBody* rigidBody = new hkpRigidBody(ci);
	boxes.push_back(static_cast<hkpRigidBody*>(g_pWorld->addEntity(rigidBody)));
	}
	boxShape->removeReference();
	}
	//create the static box where the smaller box will fall
	{
	hkVector4 halfExtents(20.0f, 2.0f, 20.f);
	hkpBoxShape* boxShape = new hkpBoxShape(halfExtents);
	hkpRigidBodyCinfo ci;
	ci.m_shape = boxShape;
	ci.m_position = hkVector4(0, -2, 0);
	ci.m_motionType = hkpMotion::MOTION_FIXED;
	boxShape->setRadius(0.001f);
	hkpRigidBody* rigidBody = new hkpRigidBody(ci);
	boxShape->removeReference();
	g_pWorld->addEntity(rigidBody)->removeReference();
	}
	}

view raw MultipleBoxes_02 hosted with ❤ by GitHub

The only difference here from the previous SimpleBox demo is that now we reuse the box shape to create several rigid bodies. We alter their positions by using the loop variable. We add the rigid body to the global vector as shown in the code snippet below.

	for(int i=0;i<MAX_BOXES;++i) {
	ci.m_position = hkVector4(0.0f,5.0f+5*i,0.0f);
	hkpRigidBody* rigidBody = new hkpRigidBody(ci);
	boxes.push_back(static_cast<hkpRigidBody*>(g_pWorld->addEntity(rigidBody)));
	}

view raw MultipleBoxes_03 hosted with ❤ by GitHub

The ShutdownHavok function
This function is changed to the following.

	void ShutdownHavok() {
	g_pWorld->markForWrite();
	for(int i=0;i<MAX_BOXES;++i) {
	boxes[i]->removeReference();
	}
	g_pWorld->unmarkForWrite();
	if(g_pWorld)
	{
	g_pWorld->markForWrite();
	g_pWorld->removeReference();
	g_pWorld = HK_NULL;
	}
	delete g_pJobQueue;
	// Clean up the thread pool
	g_pThreadPool->removeReference();
	if (g_bVdbEnabled)
	ShutdownVDB();
	boxes.clear();
	hkBaseSystem::quit();
	hkMemoryInitUtil::quit();
	}

view raw MultipleBoxes_04 hosted with ❤ by GitHub

We simply run a loop to remove reference of each rigidbody one by one.

The OnRender function
The only difference in this function is that instead of one box we run a loop to access the ith box rigid body and pass it to DrawBox function as follows.

	glEnable(GL_LIGHTING);
	for(int i=0;i<MAX_BOXES;++i) {
	DrawBox(boxes[i]);
	}
	glDisable(GL_LIGHTING);

view raw MultipleBoxes_05 hosted with ❤ by GitHub

Output
That's all, the accompanying tutorial code should give you a number of boxes falling under gravity on the grid plane colliding with each other as shown in the following figure.  

You can get the full source code from my github repo https://github.com/mmmovania/HavokPhysicsTutorials

Controls:
Left click to rotate
Middle click to zoom
Right click to pan

What's next:
In the next tutorial, I will show you how to pick a rigid body.


Enjoy!!!

Havok Physics Engine Tutorial Series: A Simple Bouncing Box

OK Now that we know how to setup Havok SDK on VisualStudio 2012. If not, you may want to go through my first tutorial which tells you how to get started with the new Havok SDK. In this tutorial, similar to my PhysX Tutorials, I start with a very simple box falling on the floor due to gravity.

In this tutorial, the InitializeHavok function is changed to the following.

	void InitializeHavok() {
	InitMemory();
	InitPhysicalWorld();
	if (g_bVdbEnabled)
	InitVDB();
	//new additions in SimpleBox
	g_pWorld->lock();
	AddRigidBodies();
	g_pWorld->unlock();
	}

view raw SimpleBox_01 hosted with ❤ by GitHub

I have highlighted the new additions in bold. Similar to the first tutorial, we initialize the memory, physics world settings and the visual debugger. After that, we lock the Havok physics world, then add our rigid bodies to the world and then we unlock the world. This is done to ensure that no two threads may access the world at the same time. All functions that may effect the physics world should be called by sandwiching them inside the world lock/unlock call.

Adding rigid bodies:
Similar to other physics engines like PhysX and Bullet, Havok also adds rigid bodies to the physics world by first specifying shape and its properties like moment of inertia tensor. In the tutorial, I add these in using the AddRigidBodies function which is defined as follows.

	void AddRigidBodies() {
	//add the falling box
	{
	hkVector4 halfExtents(0.5f, 0.5f, 0.5f);
	hkpBoxShape* boxShape = new hkpBoxShape(halfExtents);
	hkpRigidBodyCinfo ci;
	ci.m_shape = boxShape;
	ci.m_position = hkVector4(0, 10, 0);
	ci.m_motionType = hkpMotion::MOTION_DYNAMIC;
	boxShape->setRadius(0.001f);
	const hkReal boxMass(10.0f);
	hkMassProperties massProps;
	hkpInertiaTensorComputer::computeShapeVolumeMassProperties(boxShape, boxMass, massProps);
	ci.setMassProperties(massProps);
	hkpRigidBody* rigidBody = new hkpRigidBody(ci);
	boxShape->removeReference();
	box = static_cast<hkpRigidBody*>(g_pWorld->addEntity(rigidBody)));
	rigidBody->removeReference();
	}
	//create the static box where the smaller box will fall
	{
	hkVector4 halfExtents(20.0f, 2.0f, 20.f);
	hkpBoxShape* boxShape = new hkpBoxShape(halfExtents);
	hkpRigidBodyCinfo ci;
	ci.m_shape = boxShape;
	ci.m_position = hkVector4(0, -2, 0);
	ci.m_motionType = hkpMotion::MOTION_FIXED;
	boxShape->setRadius(0.001f);
	hkpRigidBody* rigidBody = new hkpRigidBody(ci);
	boxShape->removeReference();
	g_pWorld->addEntity(rigidBody)->removeReference();
	}
	}

view raw SimpleBox_02 hosted with ❤ by GitHub

Lets take a closer look at this function piece by piece. In this tutorial, we will allow a box to fall under gravity on a flat floor. We have purposely created two scopes to highlight the two added rigid bodies. In the first scope, we create a dynamic rigid body, a simple box that falls under gravity.

	hkVector4 halfExtents(0.5f, 0.5f, 0.5f);
	hkpBoxShape* boxShape = new hkpBoxShape(halfExtents);

view raw SimpleBox_03 hosted with ❤ by GitHub

In the above lines, we first specify the box's shape(hkpBoxShape) by passing it the half extents of the box.

	hkpRigidBodyCinfo ci;
	ci.m_shape = boxShape;
	ci.m_position = hkVector4(0, 10, 0);
	ci.m_motionType = hkpMotion::MOTION_DYNAMIC;
	boxShape->setRadius(0.001f);

view raw SimpleBox_04 hosted with ❤ by GitHub

In the above lines, after we specify the box's shape, we have to fill in the rigid body info structure which stores the shape as well as other rigid body dynamics properties like mass, moment of inertia, position, motion type etc. We then set its motion state to be dynamic (MOTION_DYNAMIC). If we want to make the rigid body static, we set the motion state to MOTION_FIXED. For all physics object, a collision margin is specified which controls the offset at which the collision is detected. We set the value to 0.001 which creates a very thin offset. If we left this line out, the box will collide at a much higher position than what is visible.

	const hkReal boxMass(10.0f);
	hkMassProperties massProps;
	hkpInertiaTensorComputer::computeShapeVolumeMassProperties(boxShape, boxMass, massProps);
	ci.setMassProperties(massProps);

view raw SimpleBox_05 hosted with ❤ by GitHub

In the above lines, we first calculate the box's moment of inertia and then specify the mass properties.
Once the rigidbody cinfo structure is filled, we can create a rigid body from it using the following code.

	hkpRigidBody* rigidBody = new hkpRigidBody(ci);
	boxShape->removeReference();
	box = static_cast<hkpRigidBody*>(g_pWorld->addEntity(rigidBody));
	rigidBody->removeReference();

view raw SimpleBox_06 hosted with ❤ by GitHub

Note that once the rigid body has been created, we can delete the box shape. We do so by making a call to removeReference function. This is the recommended approach to delete objects in Havok rather than calling the delete function.  After this call, we add the rigid body to the physics world by calling g_pWorld->addEntity function passing it the rigid body. On success, the returned reference contains the created rigid body entity. After this call, we can safely remove the reference to the rigid body to ensure that we always have a single reference to the rigid body.

In case of the static box (our floor), the only difference is in the motion state which is set at MOTION_FIXED. The rest of the code is same.

Obtaining the box's transform:
If all goes well, once we add the rigid body into the world, it starts to simulate. The step function calculates the next transform for all rigid bodies in the world taking into account their collisions and collision responses. For rendering or any other purposes, we need to know the local transform of the dynamic box. In order to prevent two threads from accessing the same object, Havok mandates that all world objects and their properties should be read inside a hkpWorld->lockForRead/unlockForRead fuction pair as follows.

	g_pWorld->markForRead();
	{
	//code for reading of rigid bodies and their properties
	}
	g_pWorld->unmarkForRead();

view raw SimpleBox_07 hosted with ❤ by GitHub

Now to access the matrix of the given box, we first call the box rigid body's approxCurrentTransform function that calculates the box's transformation matrix. Next, we store the matrix into a local variable. Using the local matrix varible, we fill our float array which we then pass to the rendering API. 

	hkMatrix4 m;
	float mat[16]={};
	hkTransform transform;
	box->approxCurrentTransform( transform );
	m.set( transform );
	memcpy(mat, &m, sizeof(float)*16);

view raw SimpleBox_08 hosted with ❤ by GitHub

Drawing the box:
If we are successful, the mat array contains our 4x4 matrix. In OpenGL, we can then multiply this matrix to the current modelview matrix and then the object is placed and oriented in the 3D graphics world using the transform calculated by the Havok physics engine.

	glPushMatrix();
	glMultMatrixf(mat);
	glutSolidCube(1);
	glPopMatrix();

view raw SimpleBox_09 hosted with ❤ by GitHub

For all those modern OpenGL lovers, I know this is legacy OpenGL code but my point is to preset the concept. You can convert this to modern OpenGL if required without a problem by using math libraries like glm.

Output
That's all, the accompanying tutorial code should give you a simple box falling under gravity on the grid plane as shown in the following figure.  

You can get the full source code from my github repo https://github.com/mmmovania/HavokPhysicsTutorials

Controls:
Left click to rotate
Middle click to zoom
Right click to pan

What's next:
In the next tutorial, I will show you how to add multiple boxes.


Enjoy!!!

Havok Physics Engine Tutorial Series: Getting Started

Hi all,
I am starting a new tutorial series on Havok which is an industry standard physics engine. A couple of its components namely the Havok Physics Engine and the Havok Animation SDKs were recently released free (binary only) with sponsorship of Intel under the following terms and conditions.

Havok's Intel® sponsored free binary-only PC download can be used during development to evaluate, prototype, and commercially release any PC game. There are some basic licensing rules to follow:

• PC titles sold for a retail value of less than $10.00 USD do not require a Havok distribution license to be executed.
• PC titles sold for a retail value of more than $10.00 USD or more do require a Havok license to be executed but at no additional cost.

Details here: http://software.intel.com/sites/havok/en/

As always, here is the disclaimer, I am not a Havok employee nor do I represent Havok. I am a hobbyist programmer who is trying to fill the online void in the OpenGL world. I have noticed this in the Havok physics engine case as well. So I am trying to make it easier for other OpenGL programmers to get up and running with Havok Physics SDK. All information contained in these tutorials is based on concepts gained from, the two tutorials cited below, the excellent Havok physics user guide and sample demos. These tutorials are written with clarity in mind showing clearly what is required to get started with Havok Physics SDK in VisualStudio 2012 on Windows 7. Note that there might be better and more optimized paths for these tutorials and I hope users will spot those in the comments below the tutorials.

When I started out with Havok Physics SDK, I was really surprised with the detailed documentation given with the Havok Physics SDK. This includes a lot of sample demos which show a lot of varied real-time physics concepts. Unfortunately though, the Havok Physics SDK uses its own DirectX based framework. There are no OpenGL demos in there. Ofcourse hiding the details behind a framework is good but it makes understanding of minute details difficult and you have to dive in code to know what is really required. Therefore, I like other programmers started out but then was frustrated to get up. I went online to find some information and luckily got these two links

http://piotrpluta.opol.pl/programming/havok-physics/
http://dalyup.wordpress.com/2012/02/07/havok-tutorial-01-getting-started/

Both of these cover the basics really well including how to get started from scratch with Havok Physics. I will hope that the readers of this blogs will follow these two tutorials first before proceeding forward. The issue with these is that the Havok sdk has changed a bit and there are some more changes that are required to get up and running with the latest free Havok physics sdk, So here are the missing links.

For all of the tutorial series, I will assume that VisualStudio2012 is used and that you have downloaded the Havok sdk and freeglut libraries to some place on your harddisk. To make it  smoother to follow, I suggest you create two environment variables

1. HAVOK_ROOT (pointing to the root folder of Havok sdk typically named by date for e.g. hk2013_1_0_r1)
2. LIBRARIES_ROOT  (generic folder for e.g. E:\Libraries containing the freeglut root folder e.g. E:\Libraries\freeglut-2.8.1) 

Compiler Settings
OK once this is done, you need to add the following paths to the includes directory (
C/C++ General->Additional Include Directories)
 $(LIBRARIES_ROOT)\freeglut-2.8.1\include;$(HAVOK_ROOT)/Source;%(AdditionalIncludeDirectories)

In addition, add HK_CONFIG_SIMD=1 in preprocessor definitions (C/C++->Preprocessor->Preprocessor Definitions)

Also change the C/C++->Code Generation page so that it appears as shown in the figure below

Linker Settings
Change Additional Libraries Directories (Linker->General->Additional Libraries Directories) to $(LIBRARIES_ROOT)\freeglut-2.8.1\lib\x86\Debug\;$(HAVOK_ROOT)\Lib\win32_vs2012_win7\debug_dll;%(AdditionalLibraryDirectories)

In Linker->Input->Additional Dependencies, add
hkBase.lib
hkCompat.lib
hkGeometryUtilities.lib
hkInternal.lib
hkSerialize.lib
hkVisualize.lib
hkaInternal.lib
hkaAnimation.lib
hkaPhysics2012Bridge.lib
hkpConstraint.lib
hkpConstraintSolver.lib
hkpCollide.lib
hkpDynamics.lib
hkpInternal.lib
hkpUtilities.lib
hkpVehicle.lib
hkcdCollide.lib
hkcdInternal.lib

That's it for the compiler and linker settings now to the real code. First, the include files

Include Files
To get our very first tutorial up, we need to include the following headers

	#include <GL/freeglut.h>
	#include <iostream>
	//Havok headers
	// Keycode
	#include <Common/Base/keycode.cxx>
	#include <Common/Base/Config/hkProductFeatures.cxx>
	// Math and base includes
	#include <Common/Base/hkBase.h>
	#include <Common/Base/System/hkBaseSystem.h>
	#include <Common/Base/System/Error/hkDefaultError.h>
	#include <Common/Base/Memory/System/hkMemorySystem.h>
	#include <Common/Base/Memory/System/Util/hkMemoryInitUtil.h>
	#include <Common/Base/Memory/Allocator/Malloc/hkMallocAllocator.h>
	#include <Common/Base/Thread/Job/ThreadPool/Cpu/hkCpuJobThreadPool.h>
	#include <Physics2012/Dynamics/World/hkpWorld.h>
	#include <Physics2012/Collide/Dispatch/hkpAgentRegisterUtil.h>
	// Visual Debugger includes
	#include <Common/Visualize/hkVisualDebugger.h>
	#include <Physics2012/Utilities/VisualDebugger/hkpPhysicsContext.h>

view raw GettingStarted_13 hosted with ❤ by GitHub

Havok Initialization
A lot of code is required to initialize Havok. I will go through these step by step. Basically all Havok codes need to initialize at least two main components (Havok Memory Settings and Havok Physical World Settings) and an optional third component (the Visual Debugger). I create an InitialzieHavok function which is defined as follows

	void InitializeHavok() {
	InitMemory();
	InitPhysicalWorld();
	if (g_bVdbEnabled)
	InitVDB();
	}

view raw GettingStarted_01 hosted with ❤ by GitHub

Basically, this code is just calling the three initialization functions. Here we will look at each of these one by one.

(a) Initializing Havok Memory routines       
I create a simple function InitMemory which does the memory initialization.  Here is how the InitMemory function is defined.

	void InitMemory() {
	#if defined(HK_COMPILER_HAS_INTRINSICS_IA32) &&
	HK_CONFIG_SIMD == HK_CONFIG_SIMD_ENABLED
	// Flush all denormal/subnormal numbers to zero.
	// Operations on denormals are very slow,
	// up to 100 times slower than normal numbers.
	_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
	#endif
	// Initialize the base system including our memory system
	// Allocate 0.5MB of physics solver buffer.
	hkMemoryRouter* memoryRouter = hkMemoryInitUtil::initDefault( hkMallocAllocator::m_defaultMallocAllocator, hkMemorySystem::FrameInfo( 500* 1024 ) );
	hkBaseSystem::init( memoryRouter, OnError );
	// We can cap the number of threads used - here we use the maximum for whatever multithreaded platform we are running on. This variable is
	// set in the following code sections.
	int totalNumThreadsUsed;
	// Get the number of physical threads available on the system
	hkHardwareInfo hwInfo;
	hkGetHardwareInfo(hwInfo);
	totalNumThreadsUsed = hwInfo.m_numThreads;
	// We use one less than this for our thread pool, because we must also use this thread for our simulation
	hkCpuJobThreadPoolCinfo threadPoolCinfo;
	threadPoolCinfo.m_numThreads = totalNumThreadsUsed - 1;
	// This line enables timers collection, by allocating 200 Kb per thread. If you leave this at its default (0),
	// timer collection will not be enabled.
	threadPoolCinfo.m_timerBufferPerThreadAllocation = 200000;
	g_pThreadPool = new hkCpuJobThreadPool( threadPoolCinfo );
	// We also need to create a Job queue. This job queue will be used by all Havok modules to run multithreaded work.
	// Here we only use it for physics.
	hkJobQueueCinfo info;
	info.m_jobQueueHwSetup.m_numCpuThreads = totalNumThreadsUsed;
	g_pJobQueue = new hkJobQueue(info);
	// Enable monitors for this thread.
	// Monitors have been enabled for thread pool threads already (see above comment).
	hkMonitorStream::getInstance().resize(200000);
	}

view raw GettingStarted_02 hosted with ❤ by GitHub

Lets see the details of the function piece by piece. First is a macro call _MM_SET_FLUSH_ZERO_MODE 
which is used to ensure that there are no subnormal numbers (numbers that are very small close to zero) because they can lead to slower performance. If you want to know more about this macro, have a look at the wikipedia entry   
http://en.wikipedia.org/wiki/Denormal_number

	hkMemoryRouter* memoryRouter = hkMemoryInitUtil::initDefault( hkMallocAllocator::m_defaultMallocAllocator, hkMemorySystem::FrameInfo( 500* 1024 ) );
	hkBaseSystem::init( memoryRouter, OnError );

view raw GettingStarted_12 hosted with ❤ by GitHub

The above calls initialize the memory subsystem by allocating 0.5MB for physics calculation with a memory router using the default allocator. We also call the Havok base system initialization function passing it our memory router and an error callback function. I just name this error callback OnError and I define it globally as follows which just dumps the passed msg to standard error stream (std::err)

	void OnError(const char* msg, void* userArgGivenToInit) {
	std::cerr < < "Report: " < < msg < < std::endl;
	}

view raw GettingStarted_03 hosted with ❤ by GitHub

Next, we create a thread pool with the given number of threads. We get the current device's capabilities to obtain the maximum number of parallel threads available on the multithreaded platform we are running on. 

	int totalNumThreadsUsed;
	// Get the number of physical threads available on the system
	hkHardwareInfo hwInfo;
	hkGetHardwareInfo(hwInfo);
	totalNumThreadsUsed = hwInfo.m_numThreads;
	// We use one less than this for our thread pool, because we must also use this thread for our simulation
	hkCpuJobThreadPoolCinfo threadPoolCinfo;
	threadPoolCinfo.m_numThreads = totalNumThreadsUsed - 1;
	// This line enables timers collection, by allocating 200 Kb per thread. If you leave this at its default (0),
	// timer collection will not be enabled.
	threadPoolCinfo.m_timerBufferPerThreadAllocation = 200000;
	g_pThreadPool = new hkCpuJobThreadPool( threadPoolCinfo );

view raw GettingStarted_04 hosted with ❤ by GitHub

We then create a job queue for Havok modules to run multithreaded work. Finally, the function ends with the creation of a stream of monitors for the thread pool.

	hkJobQueueCinfo info;
	info.m_jobQueueHwSetup.m_numCpuThreads = totalNumThreadsUsed;
	g_pJobQueue = new hkJobQueue(info);
	hkMonitorStream::getInstance().resize(200000);

view raw GettingStarted_05 hosted with ❤ by GitHub

  
(b) Initialize Physical World
The InitPhysicalWorld is defined as follows

	void InitPhysicalWorld() {
	g_pWorldInfo.m_simulationType = hkpWorldCinfo::SIMULATION_TYPE_MULTITHREADED;
	g_pWorldInfo.m_broadPhaseBorderBehaviour = hkpWorldCinfo::BROADPHASE_BORDER_REMOVE_ENTITY;
	g_pWorld = new hkpWorld(g_pWorldInfo);
	g_pWorld->m_wantDeactivation = false;
	g_pWorld->markForWrite();
	hkpAgentRegisterUtil::registerAllAgents( g_pWorld->getCollisionDispatcher() );
	g_pWorld->registerWithJobQueue( g_pJobQueue );
	}

view raw GettingStarted_06 hosted with ❤ by GitHub

In Havok, the physics simulation world is represented by a hkpWorld* object. This object is initialized by calling the hkpWorld constructor passing it the hkpWorldCInfo structure. This structure stores the global physics world settings like gravity etc. We first set the simulation type to a multi-threaded simulation. We then set the broadphase border behavior (which tells to the Havok physics engine to remove an entity if it goes out of the border). We pass the modified hkpWorldCInfo structure to the hkpWorld constructor to create our Havok physics world object. After this call, we set the deactivation flag of the hkpWorld to false to ensure that there is no deactivation of rigid bodies. At this point we have our physics world object created.

The next few calls modify the hkpWorld. To ensure that no two threads modify the shared hkpWorld instance at the same time, we issue a call to hkpWorld::markForWrite function. After this call we can issue all calls that modify the state of the physics world. We register collision dispatchers and the created job queue. Note that hkpWorld::markForWrite function call is paired with hkpWorld::unmarkForWrite call which is issued in the InitVDB detailed below. 

(c) Initialize Visual Debugger
Usually, you will need some mechanism to ensure that your physics world is behaving as expected. For debugging purposes or for checking physics simulation states, Havok provides a very useful application called the Visual Debugger in the SDK. We need to establish a connection to the running instance of the Havok Visual Debugger. This connection is established by creating an instance of the hkVisualDebugger object. This is done in the InitVDB function which is defined as follows.

	void InitVDB() {
	hkArray contexts;
	{
	g_pContext = new hkpPhysicsContext();
	hkpPhysicsContext::registerAllPhysicsProcesses();
	g_pContext->addWorld(g_pWorld);
	contexts.pushBack(g_pContext);
	g_pWorld->unmarkForWrite();
	}
	g_pVdb = new hkVisualDebugger(contexts);
	g_pVdb->serve();
	}

view raw GettingStarted_07 hosted with ❤ by GitHub

We first create a Havok Physics context object (hkpPhysicsContext). We then call the static function registerAllPhysicsProcesses function. We then add the Havok physics world to the created hkpPhysicsContext by calling hkpPhysicsContext::addWorld function. We then store the hkpPhysicsContext pointer in an hkArray object. We then call the hkpWorld::unmarkForWrite function that was paired with the hkpWorld::markForWrite in the InitPhysicalWorld function. The hkArray containing the hkpPhysicsContext object is passed to the hkVisualDebugger constructor and then the hkpVisualDebugger::serve function is called to initialize the connection with the hkVisualDebugger. The connection will be established with the running instance of the visual debugger.

Stepping the Havok Physics Engine and the Havok Visual Debugger 
In order to mode the Havok physics engine and the visual debugger forward in time, we need to make a call to the step function. This call is made in each frame before calling the render function. I name this function StepHavok and it is defined as follows.

	void StepHavok() {
	g_pWorld->stepMultithreaded( g_pJobQueue, g_pThreadPool, timeStep);
	if(g_bVdbEnabled)
	StepVDB();
	hkMonitorStream::getInstance().reset();
	g_pThreadPool->clearTimerData();
	}

view raw GettingStarted_08 hosted with ❤ by GitHub

We first call hkpWorld::stepMultithreaded function passing it the job queue and the timestep value which is given a constant step size of 1/60. Next, if the visual debugger is enabled, we step visual debugger using StepVDB function which is implemented as follows.

	void StepVDB() {
	g_pContext->syncTimers(g_pThreadPool);
	g_pVdb->step(timeStep);
	}

view raw GettingStarted_09 hosted with ❤ by GitHub

The StepVDB function first syncs timers in the thread pool and then calls the hkVisualDebugger::step function passing it the time step value which is also a constant step size of 1/60. Finally, the hkMonitorStream is reset and the time data values in the thread pool are cleared.

Havok Shutdown
The Havok Physics engine shutdown is carried out in the ShutdownHavok function. This function is defined as follows.

	void ShutdownHavok() {
	g_pWorld->markForWrite();
	g_pWorld->removeReference();
	delete g_pJobQueue;
	// Clean up the thread pool
	g_pThreadPool->removeReference();
	if (g_bVdbEnabled)
	ShutdownVDB();
	hkBaseSystem::quit();
	hkMemoryInitUtil::quit();
	}

view raw GettingStarted_10 hosted with ❤ by GitHub

We first ensure that the is thread safe deletion by calling hkpWorld::markForWrite function. Then, we call hkpWorld::removeReference which deletes the hkpWorld object. Since Havok internally keeps reference counts, the recommended approach to delete all Havok objects is to call removeReference on the object pointer instead of the delete function. Next, we delete the job queue and then call removeReference on the thread pool object. If the hkVisualDebugger is enabled, we call the ShutdownVDB function. Finally, we call hkBaseSystem and hkMemoryInitUtil interface's quit functions. The ShutdownVDB function is defined as follows.

	void ShutdownVDB() {
	g_pVdb->removeReference();
	g_pContext->removeReference();
	}

view raw GettingStarted_11 hosted with ❤ by GitHub

We first call, hkVisualDebugger::removeReference function to delete the connection to the hkVisualDebugger instance. We then delete the context pointer again by calling hkpPhysicsContext::removeReference function. 

Output
Running this tutorial does not show anything interesting. We just get a simple 3D grid rendered on screen as shown below. 

The console output shows the Havok initialization messages as shown below. It should only display the messages shown in the following figure. 

If you get any other message like some errors or stacktrace information, you are probably doing something incorrectly.

That's it for the first getting started tutorial. You can get the full source code from my github repository here: https://github.com/mmmovania/HavokPhysicsTutorials

Thanks,
Mobeen




 

GPU ray march renderer for PVR (Production Volume Rendering)

I have managed to integrate my GPU ray marcher with the awesome Production Volume Renderer (PVR) by Magnus Wrenninge on Windows 7 using Visual Studio 2012. The results are amazing as can be seen in the image below.

GPU based Ray Marcher for Production Volume Renderer

Here is the video

 

GPU Ray marching for PVR (Production Volume Rendering)

I had to put the video on vimeo as youtube is occassionally blocked in my country. The PVR code to model this volume is follows. This is based of the example python code from Chapter 1.

#include < pvr/Modeler.h >
#include < pvr/Primitives/Rasterization/PyroclasticPoint.h >
#include < pvr/Renderer.h >
#include < pvr/RaymarchSamplers/PhysicalSampler.h >
#include < pvr/Raymarchers/UniformRaymarcher.h >
#include < pvr/Camera.h >
#include < pvr/Occluders/OtfTransmittanceMapOccluder.h >
#include < pvr/Lights/PointLight.h >
#include < pvr/Volumes/VoxelVolume.h >
#include < Imath/ImathVec.h >
#include < pvr/VoxelBuffer.h >

void GenerateVolumeData(std::vector < GLubyte > & buffer, int& xdim, int& ydim, int& zdim) {
    pvr::Model::Modeler::Ptr modeler = pvr::Model::Modeler::create();
    pvr::Model::ModelerInput::Ptr input = pvr::Model::ModelerInput::create();
    pvr::Geo::Particles::Ptr parts = pvr::Geo::Particles::create();
    pvr::Geo::Geometry::Ptr geo = pvr::Geo::Geometry::create();
    pvr::Model::Prim::Rast::PyroclasticPoint::Ptr prim =    pvr::Model::Prim::Rast::PyroclasticPoint::create();
   
    parts - > add(1);
    geo - > setParticles(parts);       
     
    pvr::Util::ParamMap map;
    map.floatMap["amplitude"] = 0.5f;
    prim - > setParams(map);

    input - > setGeometry(geo);
    input - > setVolumePrimitive(prim);

    modeler - > addInput(input);   
    modeler - > updateBounds();
    modeler - > setResolution(200);
    modeler - > execute();

    pvr::VoxelBuffer::Ptr buf =    modeler - > buffer();
    Imath::V3i res = buf - > dataResolution(); 
   
    xdim = res.x;
    ydim = res.y;
    zdim = res.z;
   
    int index = 0; 
    for(pvr::VoxelBuffer::iterator i = buf - > begin(); i!= buf - > end(); ++i)
    {
        Imath::V3f value = *i; 
        buffer.push_back((GLubyte)(value.x*255));
        buffer.push_back((GLubyte)(value.y*255));
        buffer.push_back((GLubyte)(value.z*255)); 
    }
}

Making an OpenGL object look at another object in three different ways: quaternions, matrices and gluLookAt

Recently, I was trying to see how I could make one object look at another object in OpenGL. As always, I started out with google which give me the first link which listed to use quaternions. The second link showed me to use vector and matrices. I wanted to use gluLookAt but to my surprise none of the google links seem to provide me the answer. After some time, I got it working so here are the three methods that I have found and how they worked for me. I will be using glm library for math utilities. 

Method 1: Using Quaternions
Main reference: http://www.opengl-tutorial.org/intermediate-tutorials/tutorial-17-quaternions/
You first need to provide a function RotationBetweenVectors that returns the rotation as a quaternion. We define this function as given in the link above.

glm::quat RotationBetweenVectors(glm::vec3 start, glm::vec3 dest){
   start = glm::normalize(start);
   dest = glm::normalize(dest);

   float cosTheta = glm::dot(start, dest);
   glm::vec3 rotationAxis;

   if (cosTheta < -1 + 0.001f){
     rotationAxis = glm::cross(glm::vec3(0.0f, 0.0f, 1.0f), start);
     if (glm::length(rotationAxis) < 0.01 ) 
      rotationAxis = glm::cross(glm::vec3(1.0f, 0.0f, 0.0f), start);
      rotationAxis = glm::normalize(rotationAxis);
      return glm::angleAxis(180.0f, rotationAxis);
    }

    rotationAxis = glm::cross(start, dest);

    float s = sqrt( (1+cosTheta)*2 );
    float invs = 1 / s;

    return glm::quat(
        s * 0.5f,
        rotationAxis.x * invs,
        rotationAxis.y * invs,
        rotationAxis.z * invs
    );
}

After the RotationBetweenVectors function is defined, we can use it as follows. Target is the object you want to look at and object position is the position of the object that is going to look at the target.

glm::vec3 delta =  (targetPosition-objectPosition);
glm::vec3 desiredUp(0,1,0.00001);
glm::quat rot1 = RotationBetweenVectors(glm::vec3(0,0,1), delta);
glm::vec3 right = glm::cross(delta, desiredUp);
desiredUp = glm::cross(right, delta);
glm::vec3 newUp = rot1 * glm::vec3(0.0f, 1.0f, 0.0f);
glm::quat rot2 = RotationBetweenVectors(newUp, desiredUp);
glm::quat targetOrientation = rot2 * rot1;
glm::mat4 M=glm::toMat4(targetOrientation);
M[3][0]=objectPosition.x;
M[3][1]=objectPosition.y;
M[3][2]=objectPosition.z;

Now the matrix M is the desired matrix. To use this matrix, you multiply the matrix M with the current modelview matrix. I do it as follows.

glPushMatrix();
   glMultMatrix(glm::value_ptr(M));
glPopMatrix(); 

Method 2: Matrix based approach
Main Reference: http://stackoverflow.com/questions/6992541/opengl-rotation-in-given-direction.
This method uses basic vector calcualtion as follows.

glm::vec3 delta = targetPosition-objectPosition;
glm::vec3 up;
glm::vec3 direction(glm::normalize(delta));
if(abs(direction.x)< 0.00001 && abs(direction.z) < 0.00001){ 
   if(direction.y > 0)
     up = glm::vec3(0.0, 0.0, -1.0); //if direction points in +y
   else
        up = glm::vec3(0.0, 0.0, 1.0); //if direction points in -y 
   } else {
        up = glm::vec3(0.0, 1.0, 0.0); //y-axis is the general up
   }
up=glm::normalize(up);
glm::vec3 right = glm::normalize(glm::cross(up,direction));
up= glm::normalize(glm::cross(direction, right));
   
return glm::mat4(right.x, right.y, right.z, 0.0f,     
        up.x, up.y, up.z, 0.0f,                      
        direction.x, direction.y, direction.z, 0.0f, 
        objectPosition.x, objectPosition.y, objectPosition.z, 1.0f);   
Now the matrix M is the desired matrix. To use this matrix, you multiply the matrix M with the current modelview matrix. I do it as follows.

glPushMatrix();
   glMultMatrix(glm::value_ptr(M));
glPopMatrix();

Once again, rather than me explaining the theory, I would ask you to go to the original link given above for details.

Method 3: Using gluLookAt
Main Reference: None :( I found it myself

Now this is the method I was trying to find online but none of the references seem to provide the details. I wanted to use the gluLookAt function to find the look at matrix. So here is how I implemented it. In order to get the matrix from OpenGL, I store the current matrix (glPushMatrix) clear it to identitfy (glLoadIdentity), then call the gluLookAt function. Then I extract the current modelview matrix as follows

glPushMatrix(); //store the current MV matrix
   glLoadIdentity(); //clear current MV matrix
   gluLookAt(objectPosition.x,objectPosition.y,objectPosition.z,
             targetPosition.x,targetPosition.y,targetPosition.z,
             0,1,0);
   GLfloat MV[16];
   glGetFloatv(GL_MODELVIEW_MATRIX, MV);
glPopMatrix();

The gluLookAt function calculates the orientation matrix I want but it is to orient my object in eye space. In order to get the object space matrix, I need to find the inverse of this matrix. I use glm library to calculate the inverse by passing it my MV matrix as follows (please note how the MV matrix is passed to glm)

glm::mat4 M(MV[0], MV[1], MV[2], MV[3],
            MV[4], MV[5], MV[6], MV[7],
            MV[8], MV[9], MV[10], MV[11],
            MV[12], MV[13], MV[14], MV[15]);
M = glm::inverse(M);


Another thing we need to do is to invert the X and Z axes which points in the -ve X and -ve Z axis in object space after the inverse. Thus, I do the following calls.

M[0][0] = -M[0][0];
M[0][1] = -M[0][1];
M[0][2] = -M[0][2];

M[2][0] = -M[2][0];
M[2][1] = -M[2][1];
M[2][2] = -M[2][2];

Now the matrix M is the desired matrix. To use this matrix, you multiply the matrix M with the current modelview matrix. I do it as follows.

glPushMatrix();
   glMultMatrix(glm::value_ptr(M));
glPopMatrix(); 

This produces the same result as the previous two methods.  For your convenience, here are the three methods in their separate functions.

glm::mat4 GetMatrixMethod1(const glm::vec3& object, const glm::vec3& target) {
    glm::vec3 delta =  (target-object);
    glm::vec3 desiredUp(0,1,0.00001);
    glm::quat rot1 = RotationBetweenVectors(glm::vec3(0,0,1), delta);
    glm::vec3 right = glm::cross(delta, desiredUp);
    desiredUp = glm::cross(right, delta);
 
    glm::vec3 newUp = rot1 * glm::vec3(0.0f, 1.0f, 0.0f);
    glm::quat rot2 = RotationBetweenVectors(newUp, desiredUp);
    glm::quat targetOrientation = rot2 * rot1;
    glm::mat4 M=glm::toMat4(targetOrientation);
    M[3][0] = object.x;
    M[3][1] = object.y;
    M[3][2] = object.z;
    return M;
}

glm::mat4 GetMatrixMethod2(const glm::vec3& object, const glm::vec3& target) {
    //second method
    glm::vec3 up;
    glm::vec3 direction(glm::normalize(target-object));
    if(abs(direction.x)< 0.00001 && abs(direction.z) < 0.00001){  
    if(direction.y > 0)
        up = glm::vec3(0.0, 0.0, -1.0);
    else
        up = glm::vec3(0.0, 0.0, 1.0); 
} else {
        up = glm::vec3(0.0, 1.0, 0.0);      }
    up=glm::normalize(up);
   
   glm::vec3 right = glm::normalize(glm::cross(up,direction));
   up= glm::normalize(glm::cross(direction, right));
   
   return glm::mat4(right.x, right.y, right.z, 0.0f,     
        up.x, up.y, up.z, 0.0f,                      
        direction.x, direction.y, direction.z, 0.0f, 
        object.x, object.y, object.z, 1.0f);  
}

glm::mat4 GetMatrixMethod3(const glm::vec3& object, const glm::vec3& target) {
    //assuming that the current matrix mode is modelview matrix
    glPushMatrix();
        glLoadIdentity();
        gluLookAt(object.x, object.y, object.z,
                  target.x, target.y, target.z,
                  0,1,0);

        GLfloat MV[16];
        glGetFloatv(GL_MODELVIEW_MATRIX, MV);       
    glPopMatrix();

    glm::mat4 T(MV[0], MV[1], MV[2], MV[3],
                MV[4], MV[5], MV[6], MV[7],
                MV[8], MV[9], MV[10], MV[11],
                MV[12], MV[13], MV[14], MV[15]);
    T = glm::inverse(T);
    T[2][0] = -T[2][0];
    T[2][1] = -T[2][1];
    T[2][2] = -T[2][2];


    T[0][0] = -T[0][0];
    T[0][1] = -T[0][1];
    T[0][2] = -T[0][2];    
    return  T ;
}

And their usage is as follows.

glm::mat4 M = GetMatrixMethod1(boxPosition, spherePosition);
//glm::mat4 M = GetMatrixMethod2(boxPosition, spherePosition);
//glm::mat4 M = GetMatrixMethod3(boxPosition, spherePosition);

glPushMatrix();
  glTranslatef(spherePosition[0], spherePosition[1], spherePosition[2]);
    DrawAxes();
    glutWireSphere(1,10,10);
glPopMatrix();

glPushMatrix();
  glMultMatrixf(glm::value_ptr(M));
           
  DrawAxes();
  glutWireCube(1);
glPopMatrix();

Thanks,
Mobeen