/*
 * Copyright (C) 1997-2009, R3vis Corporation.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library 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
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA, or visit http://www.gnu.org/copyleft/lgpl.html.
 *
 * Original Contributor:
 *   Wes Bethel, R3vis Corporation, Marin County, California
 *   http://www.r3vis.com/
 * Additional Contributor(s):
 *
 * The OpenRM project is located at http://openrm.sourceforge.net/.
 */
/*
 * $Id: rmarcball.c,v 1.6 2005/02/26 16:10:24 wes Exp $
 * Version: $Name: v180-alpha-02 $
 * $Revision: 1.6 $
 * $Log: rmarcball.c,v $
 * Revision 1.6  2005/02/26 16:10:24  wes
 * Modified rmauxArcBall, rmauxDolly, rmauxTranslate to perform more thorough
 * error checking upon entry.
 *
 * Revision 1.5  2005/02/19 16:08:45  wes
 * Distro sync and consolidation.
 *
 * Revision 1.4  2005/01/23 17:11:49  wes
 * Copyright update to 2005.
 *
 * Revision 1.3  2004/01/17 04:07:53  wes
 * Updated copyright line for 2004.
 *
 * Revision 1.2  2003/02/02 02:07:18  wes
 * Updated copyright to 2003.
 *
 * Revision 1.1.1.1  2003/01/28 02:15:23  wes
 * Manual rebuild of rm150 repository.
 *
 * Revision 1.5  2003/01/16 22:21:18  wes
 * Updated all source files to reflect new organization of header files:
 * all header files formerly located in include/rmaux, include/rmi, include/rmv
 * are now located in include/rm.
 *
 * Revision 1.4  2002/04/30 19:37:30  wes
 * Updated copyright dates.
 *
 * Revision 1.3  2001/03/31 17:09:31  wes
 * v1.4.0-alpha2 checkin.
 *
 * Revision 1.2  2000/08/28 01:31:28  wes
 * Modified rmauxDolly, rmauxTranslate to take RMcamera3D parms as input -
 * thereby moving the responsiblity of scene parameter updates up a level
 * to the caller. Reinstated rmauxUI(). Added isometric scaling back in.
 *
 * Revision 1.1  2000/08/27 19:29:36  wes
 * Initial entry. Code contribution by jdb.
 *
 * Revision 1.2  2000/04/17 00:06:51  wes
 * Numerous documentation updates and code reorganization courtesy of jdb.
 *
 * Revision 1.1.1.1  2000/02/28 21:29:40  wes
 * OpenRM 1.2 Checkin
 *
 * Revision 1.1.1.1  2000/02/28 17:18:48  wes
 * Initial entry - pre-RM120 release, source base for OpenRM 1.2.
 *
 */

/*
 * the routines in this file were inspired by a well-known
 * quaternion/trackball tutorial. there are only two routines
 * in this file that are intended to be called from outside
 * (rmauxTrackball & rmauxQuatToRotationMatrix),
 * and both of these are called from rmauxEventLoop to construct
 * the virtual trackball interface. developers are certainly
 * free to use these two routines if useful.
 *
 * 1/12/2000 w.bethel
 */

/*
 * MODIFICATIONS:
 *
 * 1. added a view plane translation (rmauxTranslate) and dolly function (rmauxDolly)
 * 2. recoded function rmauxTrackball to rmauxArcBall (a true Arcball implementation)
 * 3. updated function names and reogranized file order
 * 4. removed `static' function types
 * 5. removed rmauxQuatToRotationMatrix from public view (this is specific to Arcball)
 * 6. removed normalize_quat(), add_quats(), axis_to_quat() [deprecated for this UI]
 *
 * summary: new interface based on an Arcball implementation, along with image plane
 *          translation and camera dolly.  3 public functions (rmauxTranslate, rmauxArcBall
 *	    and rmauxDolly).
 *
 * 08/16/2000 jdb
 */

#include <rm/rm.h>
#include <rm/rmaux.h>
#include <string.h>

/* radius of virtual sphere superimposed on viewport for rotation manipulation
 * within the superimposed sphere, quaternion Arcball rotation occurs
 * outside the superimposed sphere, view plane normal rotation occurs
 */
#define ARCBALLSIZE (0.8)

/* PRIVATE declarations (vector routines [10]) */
void  vzero (float *v);
void  vset (float *v, float x, float y, float z);
void  vcopy (const float *v1, float *v2);
void  vadd (const float *src1, const float *src2, float *dst);
void  vsub (const float *src1, const float *src2, float *dst);
float vlength (const float *v);
void  vnormal (float *v);
void  vscale (float *v, float s);
float vdot (const float *v1, const float *v2);
void  vcross (const float *v1, const float *v2, float *cross);

/* PRIVATE declarations (quaternion/rotation routines [2]) */
void  project2sphere (float x, float y, float *p);
void  quat2rotation (float *q, float m[4][4]);


/*
 * ----------------------------------------------------
 * @Name rmauxArcBall
 @pstart
 void rmauxArcBall (float *x1, 
		    float *y1, 
		    float *x2, 
		    float *y2,
		    RMmatrix *result)
 @pend

 @astart
 
 float *x1, *y1 - pointers to floating point values in the range
    (-1..1).  These values define a pixel coordinate position within
    the viewport (NDC).  This point is the starting point for when
    computing a rotation. (input)

 float *x2, *y2 - pointers to floating point values in the range
    (-1..1).  These values define a pixel coordinate position within
    the viewport (NDC).  This point is the ending point for computing
    a rotation. (input)

 RMmatrix *result - pointer to an RMmatrix (result). The rotation matrix
    computed by rmauxArcBall is copied into this return parameter.
    
 @aend

 @dstart

 The virtual arcball interface implements a rotation model such that
 changes about the horizontal axis map to X-axis rotations (roughly)
 and changes about the vertical axis map (roughly) to Y-axis
 rotations.

 The points (x1,y1) and (x2,y2) are projected onto the surface of a
 sphere, then a quaternion (the result of this routine) is computed
 that reflects a rotation from point (x1,y1) to point (x2,y2) along a
 great-arc path on the surface of a sphere.

 The rotation matrix computed from (x1,y1) and (x2,y2) is copied into
 the RMmatrix parameter "return".
 
 @dend
 * ----------------------------------------------------
 */
void
rmauxArcBall (float *x1, 
              float *y1, 
              float *x2, 
              float *y2,
	      RMmatrix *result)
{
    float       p1[3], p2[3], q[4];

    if ((x1 == NULL) || (x2 == NULL) || (y1 == NULL) || (y2 == NULL) || (result == NULL))
    {
	rmError("rmauxArcBall error: one of the input parameters is NULL.");
	return;
    }

    /* compute spherical projections for (x1,y1) & (x2,y2) */
    project2sphere(*x1, *y1, p1);
    project2sphere(*x2, *y2, p2);

    /* find the axis and angle of rotation */
    vcross(p2, p1, q);
    q[3] = vdot(p1, p2);

    /* rotate UI target node */
    quat2rotation(q, result->m);
}

/*
 * ----------------------------------------------------
 * @Name rmauxDolly
 @pstart
 void 
 rmauxDolly (RMcamera3D *toModify,
             float *x1,
	     float *y1, 
	     float *x2, 
	     float *y2)
 @pend

 @astart

 RMcamera3D *toModify - a handle to an RMcamera3D. the eye point of
    the RMcamera3D is modified by this routine (modified).
    
 float *x1, *y1 - pointers to floating point values in the range
    (-1..1).  These values define a pixel coordinate position within
    the viewport (NDC).  This point is the starting point for
    computing the dolly translation. (input)

 float *x2, *y2 - pointers to floating point values in the range
    (-1..1).  These values define a pixel coordinate position within
    the viewport (NDC).  This point is the ending point for computing
    the dolly translation. (input)
 @aend

 @dstart

 The virtual dolly interface implements a camera dollying motion such
 that changes in the vertical directions exponentially control the
 camera dolly translation.  Positive motions up dolly the camera
 towards the look-at and negative motions down dolly the camera away
 from the camera look-at.

 This routine modifies the eye point of the toModify RMcamera3D
 parameter.

 @dend
 * ----------------------------------------------------
 */
void 
rmauxDolly (RMcamera3D *c,
	    float *x1, 
	    float *y1, 
	    float *x2, 
	    float *y2)
{
    float       d;

    if ((x1 == NULL) || (x2 == NULL) || (y1 == NULL) || (y2 == NULL) || (c == NULL))
    {
	rmError("rmauxDolly error: one of the input parameters is NULL.");
	return;
    }
    
    /* compute signed exponential dolly scale */
    d = (float)exp((double)((*y2 - *y1) * 2.0));

    /* scale eye point along eye->at direction */
    vsub((const float *)&(c->eye), (const float *)&(c->at), (float *)&(c->eye));
    vscale((float *)&(c->eye), d);
    vadd((const float *)&(c->eye), (const float *)&(c->at), (float *)&(c->eye));
}


/*
 * ----------------------------------------------------
 * @Name rmauxTranslate 
 @pstart
 void 
 rmauxTranslate (RMcamera3D *toModify,
                 float *x1, 
	         float *y1, 
	         float *x2, 
	         float *y2)
 @pend

 @astart

 RMcamera3D *toModify - a handle to an RMcamera3D object (modified).
    The eye and look-at points of this parameter are modified by
    this routine.

 float *x1, *y1 - pointers to floating point values in the range
    (-1..1).  These values define a pixel coordinate position within
    the viewport (NDC).  This point is the starting point for
    computing the image plane translation. (input)

 float *x2, *y2 - pointers to floating point values in the range
    (-1..1).  These values define a pixel coordinate position within
    the viewport (NDC).  This point is the ending point for computing
    the image plane translation. (input)
 @aend

 @dstart

 The virtual image plane translation implements a translation motion 
 in the image plane of the camera.  Motions in the window correspond 
 directly to the image plane translation applied to the camera eye and
 look-at points.

 This routine will modify the eye and look-at points of the input
 RMcamera3D object.

 @dend
 * ----------------------------------------------------
 */
void 
rmauxTranslate (RMcamera3D *c,
		float *x1, 
	        float *y1, 
	        float *x2, 
	        float *y2)
{
    float       dx, dy, s, d[3], u[3], v[3], norm[3];

    if ((x1 == NULL) || (x2 == NULL) || (y1 == NULL) || (y2 == NULL) || (c == NULL))
    {
	rmError("rmauxTranslate error: one of the input parameters is NULL.");
	return;
    }
    
    /* compute eye->at direction and scaled mouse motion in image plane */
    vsub((const float *)&(c->eye), (const float *)&(c->at), norm);
    s = vlength(norm) / 2;
    dx = s * (*x2 - *x1);
    dy = s * (*y2 - *y1);
    
    /* compute basis for translation [u,v,norm] */
    vnormal(norm);
    vcross((const float *)&(c->up), norm, u);
    vnormal(u);
    vcross(norm, u, v);

    /* compute image plane translation vector */
    vscale(u, dx);
    vscale(v, dy);
    vadd(u, v, d);

    /* translate eye & at points */
    vsub((const float *)&(c->eye), d, (float *)&(c->eye));
    vsub((const float *)&(c->at), d, (float *)&(c->at));
}


/* PRIVATE
 *
 * zeroes the (x,y,z) components of vector v
 */
void
vzero (float *v)
{
    v[0] = v[1] = v[2] = 0.0;
}


/* PRIVATE
 *
 * sets the (x,y,z) components of vector v
 */
void
vset (float *v, 
      float x, 
      float y, 
      float z)
{
    v[0] = x;
    v[1] = y;
    v[2] = z;
}


/* PRIVATE 
 *
 * copies vector v1 to vector v2
 */
void
vcopy (const float *v1, 
       float *v2)
{
    /* vectors are 3 floats each */
    memcpy(v2, v1, (3 * sizeof(float)));
}


/* PRIVATE
 *
 * computes the sum of vectors src1 & src2:
 * dst = src1 + src2
 */
void
vadd (const float *src1, 
      const float *src2, 
      float *dst)
{
   int i;

   for (i = 0; i < 3; i++)
       dst[i] = src1[i] + src2[i];
}


/* PRIVATE
 *
 * computes the difference of vectors src1 & src2:
 * dst = src1 - src2
 */
void
vsub (const float *src1, 
      const float *src2, 
      float *dst)
{
    int i;

    for (i = 0; i < 3; i++)
        dst[i] = src1[i] - src2[i];
}


/* PRIVATE
 *
 * computes magnitude (length) of vector v, ||v||
 */
float
vlength (const float *v)
{
    double t;

    /* compute as double for numerical accuracy */
    t = (v[0] * v[0]) + (v[1] * v[1]) + (v[2] * v[2]);
    t = sqrt(t);

    return((float)t);
}


/* PRIVATE
 *
 * normalizes vector v: ||v|| = 1
 */
void
vnormal (float *v)
{
    vscale(v, (1.0 / vlength(v)));
}


/* PRIVATE
 *
 * scales vector v by scalar s
 */
void
vscale (float *v, 
	float s)
{
    v[0] *= s;
    v[1] *= s;
    v[2] *= s;
}


/* PRIVATE
 *
 * computes the dot product of vectors v1 & v2:
 * (return) = <v1, v2>
 */
float
vdot (const float *v1, 
      const float *v2)
{
    return((v1[0] * v2[0]) + (v1[1] * v2[1]) + (v1[2] * v2[2]));
}


/* PRIVATE
 *
 * computes the cross product of vectors v1 & v2:
 * cross = v1 X c2
 */
void
vcross (const float *v1, 
	const float *v2, 
	float *cross)
{
    cross[0] = (v1[1] * v2[2]) - (v1[2] * v2[1]);
    cross[1] = (v1[2] * v2[0]) - (v1[0] * v2[2]);
    cross[2] = (v1[0] * v2[1]) - (v1[1] * v2[0]);
}


/* PRIVATE
 *
 * project a viewport (x,y) pair (NDC) onto a sphere of radius ARCBALLSIZE 
 * or a hyperbolic sheet if we are away from the silhouette of the sphere
 */
void
project2sphere (float x, 
		float y,
		float *p)
{
    double t;

    /* scale (x,y) to virtual sphere size */
    p[0] = x / ARCBALLSIZE;
    p[1] = y / ARCBALLSIZE;

    /* compute projection onto virtual sphere */
    t = (p[0] * p[0]) + (p[1] * p[1]);
    if (t > 1.0) /* "outside" sphere: trace silhouette */
    {
	t = 1.0 / sqrt(t);
	p[0] *= t;
	p[1] *= t;
	p[2] = 0.0;
    }
    else /* "inside" sphere */
	p[2] = sqrt(1.0 - t);
}


/* PRIVATE
 *
 * computes a rotation matrix m from a quaternion q.  
 * a pure rotation is computed from the quaternion,
 * but a homogeneous transformation is returned, with
 * identity translation and scaling.
 */
void
quat2rotation (float *q, float m[4][4])
{
    /* rotation from quaternion */
    m[0][0] = 1.0 - (2.0 * ((q[1] * q[1]) + (q[2] * q[2])));
    m[0][1] = 2.0 * ((q[0] * q[1]) - (q[2] * q[3]));
    m[0][2] = 2.0 * ((q[2] * q[0]) + (q[1] * q[3]));

    m[1][0] = 2.0 *  ((q[0] * q[1]) + (q[2] * q[3]));
    m[1][1] = 1.0 - (2.0 * ((q[2] * q[2]) + (q[0] * q[0])));
    m[1][2] = 2.0 * ((q[1] * q[2]) - (q[0] * q[3]));

    m[2][0] = 2.0 * ((q[2] * q[0]) - (q[1] * q[3]));
    m[2][1] = 2.0 * ((q[1] * q[2]) + (q[0] * q[3]));
    m[2][2] = 1.0 - (2.0 * ((q[1] * q[1]) + (q[0] * q[0])));

    /* identity translation and scaling */
    m[0][3] = m[1][3] = m[2][3] = m[3][0] = m[3][1] = m[3][2] = 0.0;
    m[3][3] = 1.0;
}


/* EOF */
