/** Rotates the matrix according to a fictitious trackball, placed in
    the middle of the given window.
    The trackball is approximated by a Gaussian curve.
    The trackball coordinate system is: x=right, y=up, z=to viewer<BR>
    The origin of the mouse coordinates zero (0,0) is considered to be top left.
  @param width, height  window size in pixels
  @param fromX, fromY   mouse starting position in pixels
  @param toX, toY       mouse end position in pixels
*/
Matrix4 Matrix4::trackballRotation(int width, int height, int fromX, int fromY, int toX, int toY)
{
  const float TRACKBALL_SIZE = 1.3f;              // virtual trackball size (empirical value)
  Matrix4 mInv;                                   // inverse of ObjectView matrix
  Vector3 v1, v2;                                 // mouse drag positions in normalized 3D space
  float smallSize;                                // smaller window size between width and height
  float halfWidth, halfHeight;                    // half window sizes
  float angle;                                    // rotational angle
  float d;                                        // distance

  // Compute mouse coordinates in window and normalized to -1..1
  // ((0,0)=window center, (-1,-1) = bottom left, (1,1) = top right)
  halfWidth   = (float)width  / 2.0f;
  halfHeight  = (float)height / 2.0f;
  smallSize   = (halfWidth < halfHeight) ? halfWidth : halfHeight;
  v1[0]       = ((float)fromX - halfWidth)  / smallSize;
  v1[1]       = ((float)(height-fromY) - halfHeight) / smallSize;
  v2[0]       = ((float)toX   - halfWidth)  / smallSize;
  v2[1]       = ((float)(height-toY)   - halfHeight) / smallSize;

  // Compute z-coordinates on Gaussian trackball:
  d       = sqrtf(v1[0] * v1[0] + v1[1] * v1[1]);
  v1[2]   = expf(-TRACKBALL_SIZE * d * d);
  d       = sqrtf(v2[0] * v2[0] + v2[1] * v2[1]);
  v2[2]   = expf(-TRACKBALL_SIZE * d * d);

  // Compute rotational angle:
  angle = v1.angle(&v2);                          // angle = angle between v1 and v2

  // Compute rotational axis:
  v2.cross(&v1);                                  // v2 = v2 x v1 (cross product)

  // Convert axis coordinates (v2) from WCS to OCS:
  mInv.identity();
  mInv.copyRot(this);                             // copy rotational part of mv to mInv
  mInv.invertOrtho();                             // invert orthogonal matrix mInv
  v2.multiply(&mInv);                             // v2 = v2 x mInv (matrix multiplication)
  v2.normalize();                                 // normalize v2 before rotation

  // Perform acutal model view matrix modification:
  return rotate(-angle, v2[0], v2[1], v2[2]);      // rotate model view matrix
}