started a small OpenGL demo making use of Eigen's geometry features

11 files changed, 1410 insertions, 0 deletions

diff --git a/demos/CMakeLists.txt b/demos/CMakeLists.txt
index b1e206096..64298baf4 100644
--- a/demos/CMakeLists.txt
+++ b/demos/CMakeLists.txt
@@ -1,3 +1,4 @@
 IF(BUILD_DEMOS)
 ADD_SUBDIRECTORY(mandelbrot)
+ADD_SUBDIRECTORY(opengl)
 ENDIF(BUILD_DEMOS)
diff --git a/demos/opengl/CMakeLists.txt b/demos/opengl/CMakeLists.txt
new file mode 100644
index 000000000..8fd02551a
--- /dev/null
+++ b/demos/opengl/CMakeLists.txt
@@ -0,0 +1,17 @@
+
+FIND_PACKAGE(Qt4 REQUIRED)
+
+set(QT_USE_QTOPENGL TRUE)
+include(${QT_USE_FILE})
+
+SET(CMAKE_INCLUDE_CURRENT_DIR ON)
+
+INCLUDE_DIRECTORIES( ${QT_INCLUDE_DIR} )
+
+SET(quaternion_demo_SRCS  gpuhelper.cpp camera.cpp trackball.cpp quaternion_demo.cpp)
+
+QT4_AUTOMOC(${quaternion_demo_SRCS})
+
+ADD_EXECUTABLE(quaternion_demo ${quaternion_demo_SRCS})
+
+TARGET_LINK_LIBRARIES(quaternion_demo ${QT_QTCORE_LIBRARY} ${QT_QTGUI_LIBRARY}  ${QT_QTOPENGL_LIBRARY})
diff --git a/demos/opengl/README b/demos/opengl/README
new file mode 100644
index 000000000..8fb16496c
--- /dev/null
+++ b/demos/opengl/README
@@ -0,0 +1,13 @@
+
+Navigation:
+ left button:           rotate around the target
+ middle button:         zoom
+ left button + ctrl     quake rotate (rotate around camera position)
+ middle button + ctrl   walk (progress along camera's z direction)
+ left button:           pan (translate in the XY camera's plane)
+
+R : move the camera to initial position
+A : start/stop animation
+C : clear the animation
+G : add a key frame
+
diff --git a/demos/opengl/camera.cpp b/demos/opengl/camera.cpp
new file mode 100644
index 000000000..26118e6ba
--- /dev/null
+++ b/demos/opengl/camera.cpp
@@ -0,0 +1,284 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "camera.h"
+
+#include "gpuhelper.h"
+#include <GL/glu.h>
+
+#include "Eigen/LU"
+using namespace Eigen;
+
+Camera::Camera()
+    : mViewIsUptodate(false), mProjIsUptodate(false)
+{
+    mViewMatrix.setIdentity();
+    
+    mFovY = M_PI/3.;
+    mNearDist = 1.;
+    mFarDist = 50000.;
+    
+    mVpX = 0;
+    mVpY = 0;
+
+    setPosition(Vector3f::Constant(50.));
+
+    setTarget(Vector3f::Zero());
+}
+
+Camera& Camera::operator=(const Camera& other)
+{
+    mViewIsUptodate = false;
+    mProjIsUptodate = false;
+    
+    mVpX = other.mVpX;
+    mVpY = other.mVpY;
+    mVpWidth = other.mVpWidth;
+    mVpHeight = other.mVpHeight;
+
+    mTarget = other.mTarget;
+    mFovY = other.mFovY;
+    mNearDist = other.mNearDist;
+    mFarDist = other.mFarDist;
+    
+    mViewMatrix = other.mViewMatrix;
+    mProjectionMatrix = other.mProjectionMatrix;
+
+    return *this;
+}
+
+Camera::Camera(const Camera& other)
+{
+    *this = other;
+}
+
+Camera::~Camera()
+{
+}
+
+
+void Camera::setViewport(uint offsetx, uint offsety, uint width, uint height)
+{
+    mVpX = offsetx;
+    mVpY = offsety;
+    mVpWidth = width;
+    mVpHeight = height;
+    
+    mProjIsUptodate = false;
+}
+
+void Camera::setViewport(uint width, uint height)
+{
+    mVpWidth = width;
+    mVpHeight = height;
+    
+    mProjIsUptodate = false;
+}
+
+void Camera::setFovY(float value)
+{
+    mFovY = value;
+    mProjIsUptodate = false;
+}
+
+Vector3f Camera::direction(void) const
+{
+    return - (orientation() * Vector3f::UnitZ());
+}
+Vector3f Camera::up(void) const
+{
+    return orientation() * Vector3f::UnitY();
+}
+Vector3f Camera::right(void) const
+{
+    return orientation() * Vector3f::UnitX();
+}
+
+void Camera::setDirection(const Vector3f& newDirection)
+{
+    // TODO implement it computing the rotation between newDirection and current dir ?
+    Vector3f up = this->up();
+    
+    Matrix3f camAxes;
+
+    camAxes.col(2) = (-newDirection).normalized();
+    camAxes.col(0) = up.cross( camAxes.col(2) ).normalized();
+    camAxes.col(1) = camAxes.col(2).cross( camAxes.col(0) ).normalized();
+    setOrientation(Quaternionf(camAxes));
+    
+    mViewIsUptodate = false;
+}
+
+void Camera::setTarget(const Vector3f& target)
+{
+    mTarget = target;
+    if (!mTarget.isApprox(position()))
+    {
+        Vector3f newDirection = mTarget - position();
+        setDirection(newDirection.normalized());
+    }
+}
+
+void Camera::setPosition(const Vector3f& p)
+{
+    mFrame.position = p;
+    mViewIsUptodate = false;
+}
+
+void Camera::setOrientation(const Quaternionf& q)
+{
+    mFrame.orientation = q;
+    mViewIsUptodate = false;
+}
+
+void Camera::setFrame(const Frame& f)
+{
+  mFrame = f;
+  mViewIsUptodate = false;
+}
+
+void Camera::rotateAroundTarget(const Quaternionf& q)
+{
+    Matrix4f mrot, mt, mtm;
+    
+    // update the transform matrix
+    updateViewMatrix();
+    Vector3f t = mViewMatrix * mTarget;
+
+    mViewMatrix = Translation3f(t)
+                * q
+                * Translation3f(-t)
+                * mViewMatrix;
+    
+    Quaternionf qa(mViewMatrix.linear());
+    qa = qa.conjugate();
+    setOrientation(qa);
+    setPosition(- (qa * mViewMatrix.translation()) );
+
+    mViewIsUptodate = true;
+}
+
+void Camera::zoom(float d)
+{
+    float dist = (position() - mTarget).norm();
+    if(dist > d)
+    {
+        setPosition(position() + direction() * d);
+        mViewIsUptodate = false;
+    }
+}
+
+void Camera::localTranslate(const Vector3f& t)
+{
+  Vector3f trans = orientation() * t;
+  setPosition( position() + trans );
+  setTarget( mTarget + trans );
+
+  mViewIsUptodate = false;
+}
+
+void Camera::localRotate(float dTheta, float dPhi)
+{
+    float dist = (position() - mTarget).norm();
+
+    setOrientation(  AngleAxisf(dTheta, up())
+                   * AngleAxisf(dPhi, right())
+                   * orientation());
+    mTarget = position() + dist * direction();
+    
+    mViewIsUptodate = false;
+}
+
+void Camera::updateViewMatrix(void) const
+{
+    if(!mViewIsUptodate)
+    {
+        Quaternionf q = orientation().conjugate();
+        mViewMatrix.linear() = q.toRotationMatrix();
+        mViewMatrix.translation() = - (mViewMatrix.linear() * position());
+
+        mViewIsUptodate = true;
+    }
+}
+
+const Transform3f& Camera::viewMatrix(void) const
+{
+  updateViewMatrix();
+  return mViewMatrix;
+}
+
+void Camera::updateProjectionMatrix(void) const
+{
+  if(!mProjIsUptodate)
+  {
+    mProjectionMatrix.setIdentity();
+    float aspect = float(mVpWidth)/float(mVpHeight);
+    float theta = mFovY*0.5;
+    float range = mFarDist - mNearDist;
+    float invtan = 1./tan(theta);
+
+    mProjectionMatrix(0,0) = invtan / aspect;
+    mProjectionMatrix(1,1) = invtan;
+    mProjectionMatrix(2,2) = -(mNearDist + mFarDist) / range;
+    mProjectionMatrix(3,2) = -1;
+    mProjectionMatrix(2,3) = -2 * mNearDist * mFarDist / range;
+    mProjectionMatrix(3,3) = 0;
+    
+    mProjIsUptodate = true;
+  }
+}
+
+const Matrix4f& Camera::projectionMatrix(void) const
+{
+  updateProjectionMatrix();
+  return mProjectionMatrix;
+}
+
+void Camera::activateGL(void)
+{
+  glViewport(vpX(), vpY(), vpWidth(), vpHeight());
+  gpu.loadMatrix(projectionMatrix(),GL_PROJECTION);
+  gpu.loadMatrix(viewMatrix().matrix(),GL_MODELVIEW);
+}
+
+
+Vector3f Camera::unProject(const Vector2f& uv, float depth) const
+{
+    Matrix4f inv = mViewMatrix.inverse();
+    return unProject(uv, depth, inv);
+}
+
+Vector3f Camera::unProject(const Vector2f& uv, float depth, const Matrix4f& invModelview) const
+{
+    updateViewMatrix();
+    updateProjectionMatrix();
+    
+    Vector3f a(2.*uv.x()/float(mVpWidth)-1., 2.*uv.y()/float(mVpHeight)-1., 1.);
+    a.x() *= depth/mProjectionMatrix(0,0);
+    a.y() *= depth/mProjectionMatrix(1,1);
+    a.z() = -depth;
+    // FIXME /\/|
+    Vector4f b = invModelview * Vector4f(a.x(), a.y(), a.z(), 1.);
+    return Vector3f(b.x(), b.y(), b.z());
+}
diff --git a/demos/opengl/camera.h b/demos/opengl/camera.h
new file mode 100644
index 000000000..46c709a93
--- /dev/null
+++ b/demos/opengl/camera.h
@@ -0,0 +1,130 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_CAMERA_H
+#define EIGEN_CAMERA_H
+
+#include <Eigen/Geometry>
+#include <QObject>
+// #include <frame.h>
+
+class Frame
+{
+  public:
+    inline Frame(const Eigen::Vector3f& pos = Eigen::Vector3f::Zero(),
+                 const Eigen::Quaternionf& o = Eigen::Quaternionf())
+      : orientation(o), position(pos)
+    {}
+    Frame lerp(float alpha, const Frame& other) const
+    {
+      return Frame((1.f-alpha)*position + alpha * other.position,
+                   orientation.slerp(alpha,other.orientation));
+    }
+
+    Eigen::Quaternionf orientation;
+    Eigen::Vector3f position;
+};
+
+class Camera
+{
+  public:
+
+    Camera(void);
+    
+    Camera(const Camera& other);
+    
+    virtual ~Camera();
+    
+    Camera& operator=(const Camera& other);
+    
+    void setViewport(uint offsetx, uint offsety, uint width, uint height);
+    void setViewport(uint width, uint height);
+    
+    inline uint vpX(void) const { return mVpX; }
+    inline uint vpY(void) const { return mVpY; }
+    inline uint vpWidth(void) const { return mVpWidth; }
+    inline uint vpHeight(void) const { return mVpHeight; }
+
+    inline float fovY(void) const { return mFovY; }
+    void setFovY(float value);
+    
+    void setPosition(const Eigen::Vector3f& pos);
+    inline const Eigen::Vector3f& position(void) const { return mFrame.position; }
+
+    void setOrientation(const Eigen::Quaternionf& q);
+    inline const Eigen::Quaternionf& orientation(void) const { return mFrame.orientation; }
+
+    void setFrame(const Frame& f);
+    const Frame& frame(void) const { return mFrame; }
+    
+    void setDirection(const Eigen::Vector3f& newDirection);
+    Eigen::Vector3f direction(void) const;
+    void setUp(const Eigen::Vector3f& vectorUp);
+    Eigen::Vector3f up(void) const;
+    Eigen::Vector3f right(void) const;
+    
+    void setTarget(const Eigen::Vector3f& target);
+    inline const Eigen::Vector3f& target(void) { return mTarget; }
+    
+    const Eigen::Transform3f& viewMatrix(void) const;
+    const Eigen::Matrix4f& projectionMatrix(void) const;
+    
+    void rotateAroundTarget(const Eigen::Quaternionf& q);
+    void zoom(float d);
+    
+    void localTranslate(const Eigen::Vector3f& t);
+    void localRotate(float dTheta, float dPhi);
+    
+    /** Setup OpenGL matrices and viewport */
+    void activateGL(void);
+    
+    Eigen::Vector3f unProject(const Eigen::Vector2f& uv, float depth, const Eigen::Matrix4f& invModelview) const;
+    Eigen::Vector3f unProject(const Eigen::Vector2f& uv, float depth) const;
+    
+  protected:
+    void updateViewMatrix(void) const;
+    void updateProjectionMatrix(void) const;
+
+  protected:
+
+    uint mVpX, mVpY;
+    uint mVpWidth, mVpHeight;
+
+    Frame mFrame;
+    
+    mutable Eigen::Transform3f mViewMatrix;
+    mutable Eigen::Matrix4f mProjectionMatrix;
+
+    mutable bool mViewIsUptodate;
+    mutable bool mProjIsUptodate;
+
+    // used by rotateAroundTarget
+    Eigen::Vector3f mTarget;
+    
+    float mFovY;
+    float mNearDist;
+    float mFarDist;
+};
+
+#endif // EIGEN_CAMERA_H
diff --git a/demos/opengl/gpuhelper.cpp b/demos/opengl/gpuhelper.cpp
new file mode 100644
index 000000000..71348737c
--- /dev/null
+++ b/demos/opengl/gpuhelper.cpp
@@ -0,0 +1,211 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "gpuhelper.h"
+#include <GL/glu.h>
+// PLEASE don't look at this old code... ;)
+
+#include <fstream>
+#include <algorithm>
+
+GpuHelper gpu;
+
+//--------------------------------------------------------------------------------
+// icosahedron
+//--------------------------------------------------------------------------------
+#define X .525731112119133606
+#define Z .850650808352039932
+
+static GLfloat vdata[12][3] = {
+   {-X, 0.0, Z}, {X, 0.0, Z}, {-X, 0.0, -Z}, {X, 0.0, -Z},
+   {0.0, Z, X}, {0.0, Z, -X}, {0.0, -Z, X}, {0.0, -Z, -X},
+   {Z, X, 0.0}, {-Z, X, 0.0}, {Z, -X, 0.0}, {-Z, -X, 0.0}
+};
+
+static GLint tindices[20][3] = {
+   {0,4,1}, {0,9,4}, {9,5,4}, {4,5,8}, {4,8,1},
+   {8,10,1}, {8,3,10}, {5,3,8}, {5,2,3}, {2,7,3},
+   {7,10,3}, {7,6,10}, {7,11,6}, {11,0,6}, {0,1,6},
+   {6,1,10}, {9,0,11}, {9,11,2}, {9,2,5}, {7,2,11} };
+//--------------------------------------------------------------------------------
+
+
+GpuHelper::GpuHelper()
+{
+    mVpWidth = mVpHeight = 0;
+    mCurrentMatrixTarget = 0;
+    mInitialized = false;
+}
+
+GpuHelper::~GpuHelper()
+{
+}
+
+void GpuHelper::pushProjectionMode2D(ProjectionMode2D pm)
+{
+    // switch to 2D projection
+    pushMatrix(Matrix4f::Identity(),GL_PROJECTION);
+
+    if(pm==PM_Normalized)
+    {
+        //glOrtho(-1., 1., -1., 1., 0., 1.);
+    }
+    else if(pm==PM_Viewport)
+    {
+        GLint vp[4];
+        glGetIntegerv(GL_VIEWPORT, vp);
+        glOrtho(0., vp[2], 0., vp[3], -1., 1.);
+    }
+
+    pushMatrix(Matrix4f::Identity(),GL_MODELVIEW);
+}
+
+void GpuHelper::popProjectionMode2D(void)
+{
+    popMatrix(GL_PROJECTION);
+    popMatrix(GL_MODELVIEW);
+}
+
+void GpuHelper::drawVector(const Vector3f& position, const Vector3f& vec, const Color& color, float aspect /* = 50.*/)
+{
+    static GLUquadricObj *cylindre = gluNewQuadric();
+    glColor4fv(color.data());
+    float length = vec.norm();
+    pushMatrix(GL_MODELVIEW);
+    glTranslatef(position.x(), position.y(), position.z());
+    Vector3f ax = Matrix3f::Identity().col(2).cross(vec);
+    ax.normalize();
+    Vector3f tmp = vec;
+    tmp.normalize();
+    float angle = 180.f/M_PI * acos(tmp.z());
+    if (angle>1e-3)
+        glRotatef(angle, ax.x(), ax.y(), ax.z());
+    gluCylinder(cylindre, length/aspect, length/aspect, 0.8*length, 10, 10);
+    glTranslatef(0.0,0.0,0.8*length);
+    gluCylinder(cylindre, 2.0*length/aspect, 0.0, 0.2*length, 10, 10);
+
+    popMatrix(GL_MODELVIEW);
+}
+
+void GpuHelper::drawVectorBox(const Vector3f& position, const Vector3f& vec, const Color& color, float aspect)
+{
+    static GLUquadricObj *cylindre = gluNewQuadric();
+    glColor4fv(color.data());
+    float length = vec.norm();
+    pushMatrix(GL_MODELVIEW);
+    glTranslatef(position.x(), position.y(), position.z());
+    Vector3f ax = Matrix3f::Identity().col(2).cross(vec);
+    ax.normalize();
+    Vector3f tmp = vec;
+    tmp.normalize();
+    float angle = 180.f/M_PI * acos(tmp.z());
+    if (angle>1e-3)
+        glRotatef(angle, ax.x(), ax.y(), ax.z());
+    gluCylinder(cylindre, length/aspect, length/aspect, 0.8*length, 10, 10);
+    glTranslatef(0.0,0.0,0.8*length);
+    glScalef(4.0*length/aspect,4.0*length/aspect,4.0*length/aspect);
+    drawUnitCube();
+    popMatrix(GL_MODELVIEW);
+}
+
+void GpuHelper::drawUnitCube(void)
+{
+    static float vertices[][3] = {
+        {-0.5,-0.5,-0.5},
+        { 0.5,-0.5,-0.5},
+        {-0.5, 0.5,-0.5},
+        { 0.5, 0.5,-0.5},
+        {-0.5,-0.5, 0.5},
+        { 0.5,-0.5, 0.5},
+        {-0.5, 0.5, 0.5},
+        { 0.5, 0.5, 0.5}};
+
+    glBegin(GL_QUADS);
+    glNormal3f(0,0,-1); glVertex3fv(vertices[0]); glVertex3fv(vertices[2]); glVertex3fv(vertices[3]); glVertex3fv(vertices[1]);
+    glNormal3f(0,0, 1); glVertex3fv(vertices[4]); glVertex3fv(vertices[5]); glVertex3fv(vertices[7]); glVertex3fv(vertices[6]);
+    glNormal3f(0,-1,0); glVertex3fv(vertices[0]); glVertex3fv(vertices[1]); glVertex3fv(vertices[5]); glVertex3fv(vertices[4]);
+    glNormal3f(0, 1,0); glVertex3fv(vertices[2]); glVertex3fv(vertices[6]); glVertex3fv(vertices[7]); glVertex3fv(vertices[3]);
+    glNormal3f(-1,0,0); glVertex3fv(vertices[0]); glVertex3fv(vertices[4]); glVertex3fv(vertices[6]); glVertex3fv(vertices[2]);
+    glNormal3f( 1,0,0); glVertex3fv(vertices[1]); glVertex3fv(vertices[3]); glVertex3fv(vertices[7]); glVertex3fv(vertices[5]);
+    glEnd();
+}
+
+void _normalize(float* v)
+{
+    float s = 1.f/ei_sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
+    for (uint k=0; k<3; ++k)
+        v[k] *= s;
+}
+
+void _subdivide(float *v1, float *v2, float *v3, long depth)
+{
+    GLfloat v12[3], v23[3], v31[3];
+    GLint i;
+
+    if (depth == 0) {
+        //drawtriangle(v1, v2, v3);
+        glNormal3fv(v1);
+        glVertex3fv(v1);
+
+        glNormal3fv(v3);
+        glVertex3fv(v3);
+
+        glNormal3fv(v2);
+        glVertex3fv(v2);
+
+        return;
+    }
+    for (i = 0; i < 3; i++) {
+         v12[i] = v1[i]+v2[i];
+         v23[i] = v2[i]+v3[i];
+         v31[i] = v3[i]+v1[i];
+    }
+    _normalize(v12);
+    _normalize(v23);
+    _normalize(v31);
+    _subdivide(v1, v12, v31, depth-1);
+    _subdivide(v2, v23, v12, depth-1);
+    _subdivide(v3, v31, v23, depth-1);
+    _subdivide(v12, v23, v31, depth-1);
+}
+
+void GpuHelper::drawUnitLightSphere(int level)
+{
+  static int dlId = 0;
+  if (!dlId)
+  {
+    dlId = glGenLists(1);
+    glNewList(dlId, GL_COMPILE);
+    glBegin(GL_TRIANGLES);
+    for (int i = 0; i < 20; i++)
+    {
+      _subdivide(&vdata[tindices[i][0]][0], &vdata[tindices[i][1]][0], &vdata[tindices[i][2]][0], 1);
+    }
+    glEnd();
+    glEndList();
+  }
+  glCallList(dlId);
+}
+
+
diff --git a/demos/opengl/gpuhelper.h b/demos/opengl/gpuhelper.h
new file mode 100644
index 000000000..a3770aafc
--- /dev/null
+++ b/demos/opengl/gpuhelper.h
@@ -0,0 +1,222 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_GPUHELPER_H
+#define EIGEN_GPUHELPER_H
+
+#include <Eigen/Geometry>
+#include <GL/gl.h>
+#include <vector>
+
+using namespace Eigen;
+
+typedef Vector4f Color;
+
+class GpuHelper
+{
+  public:
+
+    GpuHelper();
+
+    ~GpuHelper();
+
+    enum ProjectionMode2D { PM_Normalized = 1, PM_Viewport = 2 };
+    void pushProjectionMode2D(ProjectionMode2D pm);
+    void popProjectionMode2D();
+
+    /** Multiply the OpenGL matrix \a matrixTarget by the matrix \a mat.
+        Essentially, this helper function automatically calls glMatrixMode(matrixTarget) if required
+        and does a proper call to the right glMultMatrix*() function according to the scalar type
+        and storage order.
+        \warning glMatrixMode() must never be called directly. If your're unsure, use forceMatrixMode().
+        \sa Matrix, loadMatrix(), forceMatrixMode()
+    */
+    template<typename Scalar, int _Flags>
+    void multMatrix(const Matrix<Scalar,4,4, _Flags, 4,4>& mat, GLenum matrixTarget);
+
+    /** Load the matrix \a mat to the OpenGL matrix \a matrixTarget.
+        Essentially, this helper function automatically calls glMatrixMode(matrixTarget) if required
+        and does a proper call to the right glLoadMatrix*() or glLoadIdentity() function according to the scalar type
+        and storage order.
+        \warning glMatrixMode() must never be called directly. If your're unsure, use forceMatrixMode().
+        \sa Matrix, multMatrix(), forceMatrixMode()
+    */
+    template<typename Scalar, int _Flags>
+    void loadMatrix(const Eigen::Matrix<Scalar,4,4, _Flags, 4,4>& mat, GLenum matrixTarget);
+
+    template<typename Scalar, typename Derived>
+    void loadMatrix(
+        const Eigen::CwiseNullaryOp<Eigen::ei_scalar_identity_op<Scalar>,Derived>&,
+        GLenum matrixTarget);
+
+    /** Make the matrix \a matrixTarget the current OpenGL matrix target.
+        Call this function before loadMatrix() or multMatrix() if you cannot guarrantee that glMatrixMode()
+        has never been called after the last loadMatrix() or multMatrix() calls.
+        \todo provides a debug mode checking the sanity of the cached matrix mode.
+    */
+    inline void forceMatrixTarget(GLenum matrixTarget) {glMatrixMode(mCurrentMatrixTarget=matrixTarget);}
+
+    inline void setMatrixTarget(GLenum matrixTarget);
+
+    /** Push the OpenGL matrix \a matrixTarget and load \a mat.
+    */
+    template<typename Scalar, int _Flags>
+    inline void pushMatrix(const Matrix<Scalar,4,4, _Flags, 4,4>& mat, GLenum matrixTarget);
+
+    template<typename Scalar, typename Derived>
+    void pushMatrix(
+        const Eigen::CwiseNullaryOp<Eigen::ei_scalar_identity_op<Scalar>,Derived>&,
+        GLenum matrixTarget);
+
+    /** Push and clone the OpenGL matrix \a matrixTarget
+    */
+    inline void pushMatrix(GLenum matrixTarget);
+
+    /** Pop the OpenGL matrix \a matrixTarget
+    */
+    inline void popMatrix(GLenum matrixTarget);
+
+    void drawVector(const Vector3f& position, const Vector3f& vec, const Color& color, float aspect = 50.);
+    void drawVectorBox(const Vector3f& position, const Vector3f& vec, const Color& color, float aspect = 50.);
+    void drawUnitCube(void);
+    void drawUnitLightSphere(int level=0);
+
+    /// draw the \a nofElement first elements
+    inline void draw(GLenum mode, uint nofElement);
+
+    /// draw a range of elements
+    inline void draw(GLenum mode, uint start, uint end);
+
+    /// draw an indexed subset
+    inline void draw(GLenum mode, const std::vector<uint>* pIndexes);
+
+protected:
+
+    void update(void);
+
+    GLuint mColorBufferId;
+    int mVpWidth, mVpHeight;
+    GLenum mCurrentMatrixTarget;
+    bool mInitialized;
+};
+
+/** Singleton shortcut
+*/
+extern GpuHelper gpu;
+
+
+/** \internal
+*/
+template<bool RowMajor, int _Flags> struct GlMatrixHelper;
+
+template<int _Flags> struct GlMatrixHelper<false,_Flags>
+{
+    static void loadMatrix(const Matrix<float, 4,4, _Flags, 4,4>&  mat) { glLoadMatrixf(mat.data()); }
+    static void loadMatrix(const Matrix<double,4,4, _Flags, 4,4>& mat) { glLoadMatrixd(mat.data()); }
+    static void multMatrix(const Matrix<float, 4,4, _Flags, 4,4>&  mat) { glMultMatrixf(mat.data()); }
+    static void multMatrix(const Matrix<double,4,4, _Flags, 4,4>& mat) { glMultMatrixd(mat.data()); }
+};
+
+template<int _Flags> struct GlMatrixHelper<true,_Flags>
+{
+    static void loadMatrix(const Matrix<float, 4,4, _Flags, 4,4>&  mat) { glLoadMatrixf(mat.transpose().eval().data()); }
+    static void loadMatrix(const Matrix<double,4,4, _Flags, 4,4>& mat) { glLoadMatrixd(mat.transpose().eval().data()); }
+    static void multMatrix(const Matrix<float, 4,4, _Flags, 4,4>&  mat) { glMultMatrixf(mat.transpose().eval().data()); }
+    static void multMatrix(const Matrix<double,4,4, _Flags, 4,4>& mat) { glMultMatrixd(mat.transpose().eval().data()); }
+};
+
+inline void GpuHelper::setMatrixTarget(GLenum matrixTarget)
+{
+    if (matrixTarget != mCurrentMatrixTarget)
+        glMatrixMode(mCurrentMatrixTarget=matrixTarget);
+}
+
+template<typename Scalar, int _Flags>
+void GpuHelper::multMatrix(const Matrix<Scalar,4,4, _Flags, 4,4>& mat, GLenum matrixTarget)
+{
+    setMatrixTarget(matrixTarget);
+    GlMatrixHelper<_Flags&Eigen::RowMajorBit, _Flags>::multMatrix(mat);
+}
+
+template<typename Scalar, typename Derived>
+void GpuHelper::loadMatrix(
+    const Eigen::CwiseNullaryOp<Eigen::ei_scalar_identity_op<Scalar>,Derived>&,
+    GLenum matrixTarget)
+{
+    setMatrixTarget(matrixTarget);
+    glLoadIdentity();
+}
+
+template<typename Scalar, int _Flags>
+void GpuHelper::loadMatrix(const Eigen::Matrix<Scalar,4,4, _Flags, 4,4>& mat, GLenum matrixTarget)
+{
+    setMatrixTarget(matrixTarget);
+    GlMatrixHelper<(_Flags&Eigen::RowMajorBit)!=0, _Flags>::loadMatrix(mat);
+}
+
+inline void GpuHelper::pushMatrix(GLenum matrixTarget)
+{
+    setMatrixTarget(matrixTarget);
+    glPushMatrix();
+}
+
+template<typename Scalar, int _Flags>
+inline void GpuHelper::pushMatrix(const Matrix<Scalar,4,4, _Flags, 4,4>& mat, GLenum matrixTarget)
+{
+    pushMatrix(matrixTarget);
+    GlMatrixHelper<_Flags&Eigen::RowMajorBit,_Flags>::loadMatrix(mat);
+}
+
+template<typename Scalar, typename Derived>
+void GpuHelper::pushMatrix(
+    const Eigen::CwiseNullaryOp<Eigen::ei_scalar_identity_op<Scalar>,Derived>&,
+    GLenum matrixTarget)
+{
+    pushMatrix(matrixTarget);
+    glLoadIdentity();
+}
+
+inline void GpuHelper::popMatrix(GLenum matrixTarget)
+{
+    setMatrixTarget(matrixTarget);
+    glPopMatrix();
+}
+
+inline void GpuHelper::draw(GLenum mode, uint nofElement)
+{
+    glDrawArrays(mode, 0, nofElement);
+}
+
+
+inline void GpuHelper::draw(GLenum mode, const std::vector<uint>* pIndexes)
+{
+    glDrawElements(mode, pIndexes->size(), GL_UNSIGNED_INT, &(pIndexes->front()));
+}
+
+inline void GpuHelper::draw(GLenum mode, uint start, uint end)
+{
+    glDrawArrays(mode, start, end-start);
+}
+
+#endif // EIGEN_GPUHELPER_H
diff --git a/demos/opengl/quaternion_demo.cpp b/demos/opengl/quaternion_demo.cpp
new file mode 100644
index 000000000..aef6083e2
--- /dev/null
+++ b/demos/opengl/quaternion_demo.cpp
@@ -0,0 +1,315 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "quaternion_demo.h"
+
+#include <Eigen/Array>
+#include <Eigen/QR>
+#include <Eigen/LU>
+
+#include <QEvent>
+#include <QMouseEvent>
+#include <QInputDialog>
+
+using namespace Eigen;
+
+
+
+template<typename T> T lerp(float t, const T& a, const T& b)
+{
+  return a*(1-t) + b*t;
+}
+
+template<> Quaternionf lerp(float t, const Quaternionf& a, const Quaternionf& b)
+{ return a.slerp(t,b); }
+
+template<> AngleAxisf lerp(float t, const AngleAxisf& a, const AngleAxisf& b)
+{
+  return AngleAxisf(lerp(t,a.angle(),b.angle()),
+                    lerp(t,a.axis(),b.axis()).normalized());
+}
+
+template<typename OrientationType>
+inline static Frame lerpFrame(float alpha, const Frame& a, const Frame& b)
+{
+  return Frame(::lerp(alpha,a.position,b.position),
+               Quaternionf(::lerp(alpha,OrientationType(a.orientation),OrientationType(b.orientation))));
+}
+
+QuaternionDemo::QuaternionDemo()
+{
+  mAnimate = false;
+  mTrackMode = TM_NO_TRACK;
+  mTrackball.setCamera(&mCamera);
+}
+
+void QuaternionDemo::grabFrame(void)
+{
+    // ask user for a time
+    bool ok = false;
+    double t = 0;
+    if (!m_timeline.empty())
+      t = (--m_timeline.end())->first + 1.;
+    t = QInputDialog::getDouble(this, "Eigen's QuaternionDemo", "time value: ",
+      t, 0, 1e3, 1, &ok);
+    if (ok)
+    {
+      Frame aux;
+      aux.orientation = mCamera.viewMatrix().linear();
+      aux.position = mCamera.viewMatrix().translation();
+      m_timeline[t] = aux;
+    }
+}
+
+void QuaternionDemo::drawScene()
+{
+  float length = 50;
+  gpu.drawVector(Vector3f::Zero(), length*Vector3f::UnitX(), Color(1,0,0,1));
+  gpu.drawVector(Vector3f::Zero(), length*Vector3f::UnitY(), Color(0,1,0,1));
+  gpu.drawVector(Vector3f::Zero(), length*Vector3f::UnitZ(), Color(0,0,1,1));
+}
+
+void QuaternionDemo::drawPath()
+{
+
+}
+
+void QuaternionDemo::animate()
+{
+  m_alpha += double(m_timer.interval()) * 1e-3;
+
+  TimeLine::const_iterator hi = m_timeline.upper_bound(m_alpha);
+  TimeLine::const_iterator lo = hi;
+  --lo;
+
+  Frame currentFrame;
+
+  if(hi==m_timeline.end())
+  {
+    // end
+    currentFrame = lo->second;
+    stopAnimation();
+  }
+  else if(hi==m_timeline.begin())
+  {
+    // start
+    currentFrame = hi->second;
+  }
+  else
+  {
+    float s = (m_alpha - lo->first)/(hi->first - lo->first);
+    currentFrame = ::lerpFrame<Eigen::Quaternionf>(s, lo->second, hi->second);
+    currentFrame.orientation.coeffs().normalize();
+  }
+
+  currentFrame.orientation = currentFrame.orientation.inverse();
+  currentFrame.position = - (currentFrame.orientation * currentFrame.position);
+  mCamera.setFrame(currentFrame);
+  
+  updateGL();
+}
+
+void QuaternionDemo::keyPressEvent(QKeyEvent * e)
+{
+    switch(e->key())
+    {
+        case Qt::Key_Up:
+            mCamera.zoom(2);
+            break;
+        case Qt::Key_Down:
+            mCamera.zoom(-2);
+            break;
+        // add a frame
+        case Qt::Key_G:
+            grabFrame();
+            break;
+        // clear the time line
+        case Qt::Key_C:
+            m_timeline.clear();
+            break;
+        // move the camera to initial pos
+        case Qt::Key_R:
+          {
+            if (mAnimate)
+              stopAnimation();
+            m_timeline.clear();
+            float duration = 3/*AngleAxisf(mCamera.orientation().inverse()
+                              * mInitFrame.orientation).angle()*/;
+            Frame aux = mCamera.frame();
+            aux.orientation = aux.orientation.inverse();
+            aux.position = mCamera.viewMatrix().translation();
+            m_timeline[0] = aux;
+            m_timeline[duration] = mInitFrame;
+          }
+        // start/stop the animation
+        case Qt::Key_A:
+            if (mAnimate)
+            {
+              stopAnimation();
+            }
+            else
+            {
+              m_alpha = 0;
+              connect(&m_timer, SIGNAL(timeout()), this, SLOT(animate()));
+              m_timer.start(1000/30);
+              mAnimate = true;
+            }
+            break;
+        default:
+            break;
+    }
+
+    updateGL();
+}
+
+void QuaternionDemo::stopAnimation()
+{
+  disconnect(&m_timer, SIGNAL(timeout()), this, SLOT(animate()));
+  m_timer.stop();
+  mAnimate = false;
+  m_alpha = 0;
+}
+
+void QuaternionDemo::mousePressEvent(QMouseEvent* e)
+{
+  mMouseCoords = Vector2i(e->pos().x(), e->pos().y());
+  switch(e->button())
+  {
+    case Qt::LeftButton:
+      if(e->modifiers()&Qt::ControlModifier)
+      {
+        mTrackMode = TM_QUAKE_ROTATE;
+      }
+      else
+      {
+        mTrackMode = TM_ROTATE_AROUND;
+        mTrackball.reset();
+        mTrackball.track(mMouseCoords);
+      }
+      break;
+    case Qt::MidButton:
+      if(e->modifiers()&Qt::ControlModifier)
+        mTrackMode = TM_QUAKE_WALK;
+      else
+        mTrackMode = TM_ZOOM;
+      break;
+    case Qt::RightButton:
+        mTrackMode = TM_QUAKE_PAN;
+      break;
+    default:
+      break;
+  }
+}
+void QuaternionDemo::mouseReleaseEvent(QMouseEvent*)
+{
+    mTrackMode = TM_NO_TRACK;
+    updateGL();
+}
+
+void QuaternionDemo::mouseMoveEvent(QMouseEvent* e)
+{
+    // tracking
+    if(mTrackMode != TM_NO_TRACK)
+    {
+        float dx =   float(e->x() - mMouseCoords.x()) / float(mCamera.vpWidth());
+        float dy = - float(e->y() - mMouseCoords.y()) / float(mCamera.vpHeight());
+
+        if(e->modifiers() & Qt::ShiftModifier)
+        {
+          dx *= 10.;
+          dy *= 10.;
+        }
+
+        switch(mTrackMode)
+        {
+          case TM_ROTATE_AROUND :
+            mTrackball.track(Vector2i(e->pos().x(), e->pos().y()));
+            break;
+          case TM_ZOOM :
+            mCamera.zoom(dy*50);
+            break;
+          case TM_QUAKE_WALK :
+            mCamera.localTranslate(Vector3f(0, 0, dy*100));
+            break;
+          case TM_QUAKE_PAN :
+            mCamera.localTranslate(Vector3f(dx*100, dy*100, 0));
+            break;
+          case TM_QUAKE_ROTATE :
+            mCamera.localRotate(-dx*M_PI, dy*M_PI);
+            break;
+          default:
+            break;
+        }
+
+        updateGL();
+    }
+
+    mMouseCoords = Vector2i(e->pos().x(), e->pos().y());
+}
+
+void QuaternionDemo::paintGL()
+{
+  glEnable(GL_DEPTH_TEST);
+  glDisable(GL_CULL_FACE);
+  glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
+  glDisable(GL_COLOR_MATERIAL);
+  glDisable(GL_BLEND);
+  glDisable(GL_ALPHA_TEST);
+  glDisable(GL_TEXTURE_1D);
+  glDisable(GL_TEXTURE_2D);
+  glDisable(GL_TEXTURE_3D);
+
+  // Clear buffers
+  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+
+  mCamera.activateGL();
+
+  drawScene();
+}
+
+void QuaternionDemo::initializeGL()
+{
+  glClearColor(1., 1., 1., 0.);
+  glLightModeli(GL_LIGHT_MODEL_LOCAL_VIEWER, 1);
+  glDepthMask(GL_TRUE);
+  glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE);
+
+  mInitFrame.orientation = mCamera.viewMatrix().linear();
+  mInitFrame.position = mCamera.viewMatrix().translation();
+}
+
+void QuaternionDemo::resizeGL(int width, int height)
+{
+    mCamera.setViewport(width,height);
+}
+
+int main(int argc, char *argv[])
+{
+  QApplication app(argc, argv);
+  QuaternionDemo demo;
+  demo.show();
+  return app.exec();
+}
+
+#include "quaternion_demo.moc"
diff --git a/demos/opengl/quaternion_demo.h b/demos/opengl/quaternion_demo.h
new file mode 100644
index 000000000..11b003097
--- /dev/null
+++ b/demos/opengl/quaternion_demo.h
@@ -0,0 +1,90 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_QUATERNION_DEMO_H
+#define EIGEN_QUATERNION_DEMO_H
+
+#include "gpuhelper.h"
+#include "camera.h"
+#include "trackball.h"
+#include <map>
+#include <QTimer>
+#include <QtGui/QApplication>
+#include <QtOpenGL/QGLWidget>
+
+class QuaternionDemo : public QGLWidget
+{
+  Q_OBJECT
+
+    typedef std::map<float,Frame> TimeLine;
+    TimeLine m_timeline;
+    Frame lerpFrame(float t);
+
+    Frame mInitFrame;
+    bool mAnimate;
+    float m_alpha;
+
+    
+    enum TrackMode {
+      TM_NO_TRACK=0, TM_ROTATE_AROUND, TM_ZOOM,
+      TM_QUAKE_ROTATE, TM_QUAKE_WALK, TM_QUAKE_PAN
+    };
+
+    Camera mCamera;
+    TrackMode mTrackMode;
+    Vector2i mMouseCoords;
+    Trackball mTrackball;
+
+    QTimer m_timer;
+
+    void setupCamera();
+
+  protected slots:
+
+    virtual void animate(void);
+    virtual void drawScene(void);
+    virtual void drawPath(void);
+
+    virtual void grabFrame(void);
+    virtual void stopAnimation();
+
+  protected:
+
+    virtual void initializeGL();
+    virtual void resizeGL(int width, int height);
+    virtual void paintGL();
+    
+    //--------------------------------------------------------------------------------
+    virtual void mousePressEvent(QMouseEvent * e);
+    virtual void mouseReleaseEvent(QMouseEvent * e);
+    virtual void mouseMoveEvent(QMouseEvent * e);
+    virtual void keyPressEvent(QKeyEvent * e);
+    //--------------------------------------------------------------------------------
+
+  public:
+    QuaternionDemo();
+    ~QuaternionDemo() { }
+};
+
+#endif // EIGEN_QUATERNION_DEMO_H
diff --git a/demos/opengl/trackball.cpp b/demos/opengl/trackball.cpp
new file mode 100644
index 000000000..f66243d0d
--- /dev/null
+++ b/demos/opengl/trackball.cpp
@@ -0,0 +1,72 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "trackball.h"
+#include "camera.h"
+
+using namespace Eigen;
+
+void Trackball::track(const Vector2i& newPoint2D)
+{
+  if (mpCamera==0)
+    return;
+  Vector3f newPoint3D;
+  bool newPointOk = mapToSphere(newPoint2D, newPoint3D);
+
+  if (mLastPointOk && newPointOk)
+  {
+    Vector3f axis = mLastPoint3D.cross(newPoint3D).normalized();
+    float cos_angle = mLastPoint3D.dot(newPoint3D);
+    if ( ei_abs(cos_angle) < 1.0 )
+    {
+      float angle = 2.0 * acos(cos_angle);
+      mpCamera->rotateAroundTarget(Quaternionf(AngleAxisf(angle, axis)));
+    }
+  }
+
+  mLastPoint2D = newPoint2D;
+  mLastPoint3D = newPoint3D;
+  mLastPointOk = newPointOk;
+}
+
+bool Trackball::mapToSphere(const Vector2i& p2, Vector3f& v3)
+{
+  if ((p2.x() >= 0) && (p2.x() <= int(mpCamera->vpWidth())) &&
+      (p2.y() >= 0) && (p2.y() <= int(mpCamera->vpHeight())) )
+  {
+    double x  = (double)(p2.x() - 0.5*mpCamera->vpWidth())  / (double)mpCamera->vpWidth();
+    double y  = (double)(0.5*mpCamera->vpHeight() - p2.y()) / (double)mpCamera->vpHeight();
+    double sinx         = sin(M_PI * x * 0.5);
+    double siny         = sin(M_PI * y * 0.5);
+    double sinx2siny2   = sinx * sinx + siny * siny;
+
+    v3.x() = sinx;
+    v3.y() = siny;
+    v3.z() = sinx2siny2 < 1.0 ? sqrt(1.0 - sinx2siny2) : 0.0;
+
+    return true;
+  }
+  else
+    return false;
+}
diff --git a/demos/opengl/trackball.h b/demos/opengl/trackball.h
new file mode 100644
index 000000000..29413becd
--- /dev/null
+++ b/demos/opengl/trackball.h
@@ -0,0 +1,55 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_TRACKBALL_H
+#define EIGEN_TRACKBALL_H
+
+#include <Eigen/Geometry>
+
+class Camera;
+
+class Trackball
+{
+  public:
+
+    Trackball() : mpCamera(0) {}
+
+    void reset() { mLastPointOk = false; }
+
+    void setCamera(Camera* pCam) { mpCamera = pCam; }
+
+    void track(const Eigen::Vector2i& newPoint2D);
+
+  protected:
+
+    bool mapToSphere( const Eigen::Vector2i& p2, Eigen::Vector3f& v3);
+
+    Camera* mpCamera;
+    Eigen::Vector2i mLastPoint2D;
+    Eigen::Vector3f mLastPoint3D;
+    bool mLastPointOk;
+
+};
+
+#endif // EIGEN_TRACKBALL_H