glwindow.c

任意给定三维空间的点集

inline Mat4 Mat4::operator/(double s) const
{
    return Mat4(row[0]/s, row[1]/s, row[2]/s, row[3]/s);
}

inline GPoint4 Mat4::operator*(const GPoint4& v) const
{
    return GPoint4(row[0]*v, row[1]*v, row[2]*v, row[3]*v);
}
//
// Transform a homogeneous 3-vector and reproject into normal 3-space
//
inline GPoint3 Mat4::operator*(const GPoint3  & v) const
{
    GPoint4 u=GPoint4(v,1);
    double w=row[3].pnt*u;

    if(w==0.0)
	return GPoint3(row[0].pnt*u, row[1].pnt*u, row[2].pnt*u);
    else
	return GPoint3(row[0].pnt*u/w, row[1].pnt*u/w, row[2].pnt*u/w);
}

extern bool jacobi(const Mat4& m, GPoint4& vals, GPoint4 vecs[4]);

#ifdef MXGL_INCLUDED
inline void glGetMatrix(Mat4& m, GLenum which=GL_MODELVIEW_MATRIX)
{ Mat4 tmp;  glGetDoublev(which, &tmp(0,0));  m=tmp.transpose(); }

inline void glLoadMatrix(const Mat4& m)
{ Mat4 tmp = m.transpose();  glLoadMatrixd(&tmp(0,0)); }

inline void glMultMatrix(const Mat4& m)
{ Mat4 tmp = m.transpose();  glMultMatrixd(&tmp(0,0)); }
#endif



//#include "stdmix.h"
//#include "GL.h"
//#include "Arcball.h"

enum AxisSet {NoAxes, CameraAxes, BodyAxes, OtherAxes, NSets};
enum BallType { ArcBall, TrackBall };
typedef Mat4 HMatrix;



/* ***** BallAux.c ***** */

typedef  GPoint4  Quat;

Quat qOne(0,0,0,1);

/* Return quaternion product qL * qR.  Note: order is important!
 * To combine rotations, use the product Mul(qSecond, qFirst),
 * which gives the effect of rotating by qFirst then qSecond. */
Quat Qt_Mul(Quat& qL, Quat& qR)
{
    Quat qq;
    qq.pnt[W] = qL[W]*qR[W] - qL[X]*qR[X] - qL[Y]*qR[Y] - qL[Z]*qR[Z];
    qq.pnt[X] = qL[W]*qR[X] + qL[X]*qR[W] + qL[Y]*qR[Z] - qL[Z]*qR[Y];
    qq.pnt[Y] = qL[W]*qR[Y] + qL[Y]*qR[W] + qL[Z]*qR[X] - qL[X]*qR[Z];
    qq.pnt[Z] = qL[W]*qR[Z] + qL[Z]*qR[W] + qL[X]*qR[Y] - qL[Y]*qR[X];
    return (qq);
}


/* Construct rotation matrix from (possibly non-unit) quaternion.
 * Assumes matrix is used to multiply column vector on the left:
 * vnew = mat vold.  Works correctly for right-handed coordinate system
 * and right-handed rotations. */
HMatrix *Qt_ToMatrix(const Quat& q, HMatrix& out)
{
    double Nq = norm2(q);
    double s = (Nq > 0.0) ? (2.0 / Nq) : 0.0;

    double xs = q[X]*s,  ys = q[Y]*s,  zs = q[Z]*s;
    double wx = q[W]*xs, wy = q[W]*ys, wz = q[W]*zs;
    double xx = q[X]*xs, xy = q[X]*ys, xz = q[X]*zs;
    double yy = q[Y]*ys, yz = q[Y]*zs, zz = q[Z]*zs;

    out(X,X) = 1.0 - (yy + zz); out(Y,X) = xy + wz; out(Z,X) = xz - wy;
    out(X,Y) = xy - wz; out(Y,Y) = 1.0 - (xx + zz); out(Z,Y) = yz + wx;
    out(X,Z) = xz + wy; out(Y,Z) = yz - wx; out(Z,Z) = 1.0 - (xx + yy);
    out(X,W) = out(Y,W) = out(Z,W) = out(W,X) = out(W,Y) = out(W,Z) = 0.0;
    out(W,W) = 1.0;
    return ((HMatrix *)&out);
}

/* Return conjugate of quaternion. */
Quat Qt_Conj(Quat& q)
{
    Quat qq;
    qq.pnt[X] = -q[X];
    qq.pnt[Y] = -q[Y];
    qq.pnt[Z] = -q[Z];
    qq.pnt[W] = q[W];
    return (qq);
}

/* Return vector formed from components */
inline GPoint4 V3_(float x, float y, float z)
{
    return GPoint4(x, y, z, 0);
}

/* Return norm of v, defined as sum of squares of components */
inline double V3_Norm(const GPoint4& v)
{
    return ( v[X]*v[X] + v[Y]*v[Y] + v[Z]*v[Z] );
}

/* Return unit magnitude vector in direction of v */
GPoint4 V3_Unit(const GPoint4& v)
{
    static GPoint4 u(0, 0, 0, 0);
    float vlen = sqrt(V3_Norm(v));
    if (vlen != 0.0) {
	u.pnt[X] = v[X]/vlen;
        u.pnt[Y] = v[Y]/vlen; 
        u.pnt[Z] = v[Z]/vlen;
    }
    return (u);
}

/* Return version of v scaled by s */
inline GPoint4 V3_Scale(const GPoint4& v, float s)
{
    return GPoint4(s*v[X], s*v[Y], s*v[Z], v[W]);
}

/* Return negative of v */
inline GPoint4 V3_Negate(const GPoint4& v)
{
    return GPoint4(-v[X], -v[Y], -v[Z], 0);
}

/* Return sum of v1 and v2 */
inline GPoint4 V3_Add(const GPoint4& v1, const GPoint4& v2)
{
    return GPoint4(v1[X]+v2[X], v1[Y]+v2[Y], v1[Z]+v2[Z], 0);
}

/* Return difference of v1 minus v2 */
inline GPoint4 V3_Sub(const GPoint4& v1, const GPoint4& v2)
{
    return GPoint4(v1[X]-v2[X], v1[Y]-v2[Y], v1[Z]-v2[Z], 0);
}

/* Halve arc between unit vectors v0 and v1. */
GPoint4 V3_Bisect(const GPoint4& v0, const GPoint4& v1)
{
    GPoint4 v(0, 0, 0, 0);
    v = V3_Add(v0, v1);

    float Nv = V3_Norm(v);
    if (Nv < 1.0e-5) {
	v = V3_(0, 0, 1);
    } else {
	v = V3_Scale(v, 1/sqrt(Nv));
    }
    return (v);
}

/* Return dot product of v1 and v2 */
double V3_Dot(const GPoint4& v1, const GPoint4& v2)
{
    return (v1[X]*v2[X] + v1[Y]*v2[Y] + v1[Z]*v2[Z]);
}


/* Return cross product, v1 x v2 */
inline GPoint4 V3_Cross(const GPoint4& v1, const GPoint4& v2)
{
    return GPoint4(v1[Y]*v2[Z]-v1[Z]*v2[Y],
		v1[Z]*v2[X]-v1[X]*v2[Z],
		v1[X]*v2[Y]-v1[Y]*v2[X],
		0);
}



/* ***** Ball.h ***** */

typedef GPoint4 *ConstraintSet;
struct BallData {
    GPoint4 center;
    double radius, aspect;
    Quat qNow, qDown, qDrag;
    GPoint4 vNow, vDown, vFrom, vTo, vrFrom, vrTo;
    HMatrix mNow, mDown;
    int dragging;
    ConstraintSet sets[NSets];
    int setSizes[NSets];
    AxisSet axisSet;
    int axisIndex;
    BallType ballType;
};

/* Public routines */
extern BallData *Ball_Alloc();
extern void Ball_Free(BallData *);
extern void Ball_Init(BallData *ball, double aspect);
extern void Ball_Place(BallData *ball, GPoint4 center, double radius
		       /*, Quat orient*/);
extern void Ball_Mouse(BallData *ball, GPoint4 vNow);
extern void Ball_UseSet(BallData *ball, AxisSet axisSet);
extern void Ball_Update(BallData *ball, int al);
extern void Ball_Value(BallData *ball, HMatrix& mNow);
extern void Ball_BeginDrag(BallData *ball);
extern void Ball_EndDrag(BallData *ball);

extern int Ball_MouseOutsideSphere(BallData *ball, GPoint4 vMouse);
extern void Ball_SetBallType(BallData *ball, BallType bt);
/* Private routines */



/* ***** BallMath.c ***** */

/* Convert window coordinates to sphere coordinates. */
GPoint4 MouseOnSphere(GPoint4 mouse, GPoint4 ballCenter, double ballRadius)
{
    GPoint4 ballMouse;
    register double mag;
    ballMouse.pnt[X] = (mouse[X] - ballCenter[X]) / ballRadius;
    ballMouse.pnt[Y] = (mouse[Y] - ballCenter[Y]) / ballRadius;
    mag = ballMouse[X]*ballMouse[X] + ballMouse[Y]*ballMouse[Y];
    if (mag > 1.0) {
	register double scale = 1.0/sqrt(mag);
	ballMouse.pnt[X] *= scale; ballMouse.pnt[Y] *= scale;
	ballMouse.pnt[Z] = 0.0;
    } else {
	ballMouse.pnt[Z] = sqrt(1 - mag);
    }
    ballMouse.pnt[W] = 0.0;
    return (ballMouse);
}

/* Construct a unit quaternion from two points on unit sphere */
Quat Qt_FromBallPoints(GPoint4 from, GPoint4 to, BallType bt)
{
    Quat qu;
    qu.pnt[X] = from[Y]*to[Z] - from[Z]*to[Y];
    qu.pnt[Y] = from[Z]*to[X] - from[X]*to[Z];
    qu.pnt[Z] = from[X]*to[Y] - from[Y]*to[X];

    if (bt == ArcBall){
	qu.pnt[W] = from[X]*to[X] + from[Y]*to[Y] + from[Z]*to[Z];
    }
    else {
	float  a, s, mag;

	a = .5*acos(from[X]*to[X] + from[Y]*to[Y] + from[Z]*to[Z]);
	mag = sqrt(qu[X]*qu[X] + qu[Y]*qu[Y] + qu[Z]*qu[Z]);
	qu.pnt[W] = cos(a); s = sin(a);
	if (s != 0)
	  s /= mag;
	qu.pnt[X] *= s;  qu.pnt[Y] *= s; qu.pnt[Z] *= s;
    }

    return (qu);
}

/* Convert a unit quaternion to two points on unit sphere */
void Qt_ToBallPoints(Quat q, GPoint4 *arcFrom, GPoint4 *arcTo)
{
    double s;
    s = sqrt(q[X]*q[X] + q[Y]*q[Y]);
    if (s == 0.0) {
	*arcFrom = V3_(0.0, 1.0, 0.0);
    } else {
	*arcFrom = V3_(-q[Y]/s, q[X]/s, 0.0);
    }
    (*arcTo).pnt[X] = q[W]*(*arcFrom)[X] - q[Z]*(*arcFrom)[Y];
    (*arcTo).pnt[Y] = q[W]*(*arcFrom)[Y] + q[Z]*(*arcFrom)[X];
    (*arcTo).pnt[Z] = q[X]*(*arcFrom)[Y] - q[Y]*(*arcFrom)[X];
    if (q[W] < 0.0) *arcFrom = V3_(-(*arcFrom)[X], -(*arcFrom)[Y], 0.0);
}

/* Force sphere point to be perpendicular to axis. */
GPoint4 ConstrainToAxis(GPoint4& loose, GPoint4& axis)
{
    GPoint4 onPlane = V3_Sub(loose, V3_Scale(axis, V3_Dot(axis, loose)));
    register float norm = V3_Norm(onPlane);

    if (norm > 0.0) {
	if (onPlane[Z] < 0.0) onPlane = V3_Negate(onPlane);
	return ( V3_Scale(onPlane, 1/sqrt(norm)) );
    } /* else drop through */
    if (axis[Z] == 1) {
	onPlane = V3_(1.0, 0.0, 0.0);
    } else {
	onPlane = V3_Unit(V3_(-axis[Y], axis[X], 0.0));
    }
    return (onPlane);
}

/* Find the index of nearest arc of axis set. */
int NearestConstraintAxis(GPoint4& loose, GPoint4 *axes, int nAxes)
{
    GPoint4 onPlane;
    register float max, dot;
    register int i, nearest;
    max = -1; nearest = 0;
    for (i=0; i<nAxes; i++) {
	onPlane = ConstrainToAxis(loose, axes[i]);
	dot = V3_Dot(onPlane, loose);
	if (dot>max) {
	    max = dot; nearest = i;
	}
    }
    return (nearest);
}



/* ***** Ball.c ***** */
/* Ken Shoemake, 1993 */
/* Modified by Victor Ng, Jan. 1996 for OpenGL */

#define FALSE 0
#define TRUE 1

#define LG_NSEGS 4
#define NSEGS (1<<LG_NSEGS)
#define RIMCOLOR()    glColor3f(0.0f, 1.0f, 1.0f);
#define FARCOLOR()    glColor3f(0.8f, 0.5f, 0.0f);
#define NEARCOLOR()   glColor3f(1.0f, 0.8f, 0.0f);
#define DRAGCOLOR()   glColor3f(0.0f, 1.0f, 1.0f);
#define RESCOLOR()    glColor3f(0.8f, 0.0f, 0.0f);



GPoint4 otherAxis(-0.48, 0.80, 0.36, 1);

BallData *Ball_Alloc()
{
    return (BallData *)malloc( sizeof(BallData) );
}

void Ball_Free(BallData *ball)
{
    free(ball);
}

/* Establish reasonable initial values for controller. */
void Ball_Init(BallData *ball, double aspect)
{
    ball->center = qOne;
    ball->radius = 1.0;
    ball->aspect = aspect;
    ball->ballType = ArcBall;
    ball->vDown = ball->vNow = qOne;
    ball->qDown = ball->qNow = qOne;

    ball->mNow = Mat4::I;
    ball->mDown = Mat4::I;

    ball->dragging = FALSE;
    ball->axisSet = NoAxes;

    ball->sets[CameraAxes]     = (GPoint4 *)&Mat4::I(0,0);  // Mat4::I.row(0)
    ball->setSizes[CameraAxes] = 3;

    ball->sets[BodyAxes] = (GPoint4 *)&ball->mDown(0,0);  // mDown.row(0)
    ball->setSizes[BodyAxes] = 3;

    ball->sets[OtherAxes] = &otherAxis;
    ball->setSizes[OtherAxes] = 1;
}

/* Set the center and size of the controller. */
void Ball_Place(BallData *ball, GPoint4 center, double radius)
{
    ball->center = center;
    ball->radius = radius;
}

/* Incorporate new mouse position. */
void Ball_Mouse(BallData *ball, GPoint4 vNow)
{
    ball->vNow = vNow;
}

/* Choose a constraint set, or none. */
void Ball_UseSet(BallData *ball, AxisSet axisSet)
{
    /* if (!ball->dragging) */
	ball->axisSet = axisSet;
}

/* Using vDown, vNow, dragging, and axisSet, compute rotation etc. */
/* To use closest constraint axis, set aI to -1 */
void Ball_Update(BallData *ball, int aI)
{
    int setSize = ball->setSizes[ball->axisSet];
    GPoint4 *set = (GPoint4 *)(ball->sets[ball->axisSet]);
    ball->vFrom = MouseOnSphere(ball->vDown, ball->center, ball->radius);
    ball->vTo = MouseOnSphere(ball->vNow, ball->center, ball->radius);
    if (ball->dragging) {
        if (ball->axisSet!=NoAxes) {
	    aI = (aI==-1) ? ball->axisIndex : aI;
            ball->vFrom = ConstrainToAxis(ball->vFrom, set[aI]);
            ball->vTo = ConstrainToAxis(ball->vTo, set[aI]);
        }
        ball->qDrag =Qt_FromBallPoints(ball->vFrom, ball->vTo, ball->ballType);
        ball->qNow = Qt_Mul(ball->qDrag, ball->qDown);
    } else {
        if (ball->axisSet!=NoAxes) {
            ball->axisIndex = (aI==-1) ?
	                      NearestConstraintAxis(ball->vTo, set, setSize) :
	                      aI;
        }
    }
    Qt_ToBallPoints(ball->qDown, &ball->vrFrom, &ball->vrTo);
    Qt_ToMatrix(Qt_Conj(ball->qNow), ball->mNow); /* Gives transpose for GL. */
}

/* Returns TRUE iff vMouse is a point off the arcball sphere */
int Ball_MouseOutsideSphere(BallData *ball, GPoint4 vMouse)
{
    GPoint4  vTest = MouseOnSphere(vMouse, ball->center, ball->radius);

    return vTest[Z] == 0.0;
}

/* Set the BallType: 'ArcBall' rotates through _twice_ the angle
 * between the two points on the sphere, and exhibits no
 * hysteresis. 'TrackBall' rotates through exactly the angle between
 * the two points on the sphere, like a physical trackball, and
 * consequently exhibits hysteresis. */
void Ball_SetBallType(BallData *ball, BallType bt)
{
    ball->ballType = bt;
}

/* Return rotation matrix defined by controller use. */
void Ball_Value(BallData *ball, HMatrix& mNow)
{
    mNow = ball->mNow;
}


/* Begin drag sequence. */
void Ball_BeginDrag(BallData *ball)
{
    ball->dragging = TRUE;
    ball->vDown = ball->vNow;
}

/* Stop drag sequence. */
void Ball_EndDrag(BallData *ball)
{
    ball->dragging = FALSE;
    ball->qDown = ball->qNow;

    ball->mDown = ball->mNow;
}



Arcball::~Arcball()
{
    if( ball ) Ball_Free(ball);
}

void Arcball::init(const GBBox   & bounds, float aspect)
{
    if( ball ) Ball_Free(ball);  // Don't try to reuse old Balls
    ball = Ball_Alloc();
    Ball_Init(ball, aspect);

    Ball_UseSet(ball, NoAxes);
    Ball_SetBallType(ball, ArcBall);

    GPoint3  tmp_center;
    bounds.center( tmp_center );

    center = GPoint4( tmp_center, 1 );

    radius = bounds.get_radius();
    cam[X] = cam[Y] = cam[Z] = 0.0f;
    dragging = NoDrag;

    Ball_Place(ball, qOne, 0.9 /* 0.75 */);
}

bool Arcball::mouse_down(int button, int x, int y, 
                         int  width, int  height )
{
    mousex = startx = x;
    mousey = starty = y;

    //printf( "mouse_down( button: %d, x: %d, y: %d, width: %d, height: %d\n",
    //        button, x, y, width, height );
    if( button == 1 ) { 
        //int width, height;
        //glGetViewport(NULL, NULL, &width, &height);

        vNow.pnt[X] = 2.0*startx/(float)width-1.0;
        vNow.pnt[Y] = -2.0*starty/(float)height+1.0;

        Ball_Mouse(ball, vNow);

        Ball_BeginDrag(ball);

        dragging = Spin;
        //printf/( "spinning!\n" );
    }
    else if( button == 2 )
    {
        dragging = Dolly;
    }
    else if( button == 3 )
    {
        dragging = Zoom;
    }

    return true;
}

bool Arcball::mouse_up(int /*button*/, int /*x*/, int /*y*/)
{
    startx = mousex = 0;
    starty = mousey = 0;

    if( dragging == Spin )
        Ball_EndDrag(ball);

    dragging = NoDrag;

    return false;
}

bool Arcball::motion(int x, int y, int  width, int  height )
{
    if( !dragging )
        return false;

    if( dragging == Spin )
    {
        mousex = x;
        mousey = y;

        //int width, height;
        //glGetViewport(NULL, NULL, &width, &height);

        vNow.pnt[X] = 2.0*mousex/(float)width-1.0;
        vNow.pnt[Y] = -2.0*mousey/(float)height+1.0;

        Ball_Mouse(ball, vNow);
    }
    else if( dragging == Dolly )
    {
	//int width, height;
        //	glGetViewport(NULL, NULL, &width, &height);

	cam[X] += radius * (float)(x - mousex) / (float)width;
	cam[Y] += radius * (float)(mousey - y) / (float)height;

        mousex = x;
        mousey = y;
    }
    else if( dragging == Zoom )
    {
        cam[Z] += 0.02*radius*(float)(y - mousey);
        mousey = y;
    }

    return true;
}

void Arcball::apply()
{
    Ball_Update(ball, -1);

    glPushMatrix();
    glTranslatef(cam[X], cam[Y], cam[Z]);

    HMatrix mNow;
    Ball_Value(ball, mNow);

    glTranslatef(center[X], center[Y], center[Z]);
    glMultMatrixd((double *)&mNow);
    glTranslatef(-center[X], -center[Y], -center[Z]);
}

void Arcball::unapply()
{
    glPopMatrix();
}


/* glbox.C - End of File ------------------------------------------*/