
/*
 * c3dtransform.c --
 *
 *     This file contains code to calculate and apply transform matrices to
 *     vertices. This includes projection transformations.
 */

#include "c3d.h"
#include <assert.h>
#include <math.h>
#include <string.h>
#ifndef M_PI
#define M_PI (3.14159265358979323846)
#endif
#ifdef CAN3D_WGL
#define strcasecmp  stricmp
#endif

#ifndef NDEBUG
static void printTransform(C3dTransform *pTransform);
static void printMatrix(GLdouble *m);
#endif

typedef struct C3dVector C3dVector;
struct C3dVector {
    float x;
    float y;
    float z;
    float w;
};

/*
 * A C3dTransform is a 4x4 matrix. The transformation is defined as:
 *
 *     v' = Mv
 *
 * where M is the 4x4 matrix, v is a point in 3d space, and v' is the
 * transformed point in 3d space.
 *
 * A point in space is represented by a C3dVertex structure. The 'w'
 * component is implicitly 1.0 in all cases. i.e.:
 *
 *          /  x  \
 *     v =  |  y  |
 *          |  z  |
 *          \ 1.0 /
 *
 * (where x, y and z are member variables of an C3dVertex struct).
 *
 * and
 *
 *          / v1.x  v1.y  v1.z  v1.w \
 *     M =  | v2.x  v2.y  v2.z  v2.w |
 *          | v3.x  v3.y  v3.z  v3.w |
 *          \ v4.x  v4.y  v4.z  v4.w /
 *
 * (where v1 through v4 are member variables of a C3dTransform struct).
 *
 *           / x*v1.x + y*v1.y + z*v1.z + 1.0*v1.w  \
 *     v' =  | x*v2.x + y*v2.y + z*v2.z + 1.0*v2.w  |
 *           | x*v3.x + y*v3.y + z*v3.z + 1.0*v3.w  |
 *           \ x*v4.x + y*v4.y + z*v4.z + 1.0*v4.w  /
 *
 * Pointers to a C3dTransform structure are used opaquely by code outside
 * of this file, but only this file has direct access to structure data.
 */
struct C3dTransform {
    C3dVector v1;
    C3dVector v2;
    C3dVector v3;
    C3dVector v4;
};

/*
 * The three transformation primitives that everything else is built in
 * terms of. The following three values may be passed as the 'clientData'
 * argument to transformPrimitive().
 */
#define TRANSFORM_SCALE  1
#define TRANSFORM_ROTATE 2
#define TRANSFORM_MOVE   4

/*
 *---------------------------------------------------------------------------
 *
 * readVertex --
 *
 *     Read a vertex argument from a transform specification and store the
 *     results in (*pX, *pY, *pZ). A transform specification is a flat list
 *     of tokens, stored in the apElem array passed as a function argument.
 *     A vertex may be specified in one of two ways:
 *
 *         (a) As three floating point numbers, i.e. "3.0 4.0 5.0", or
 *         (b) By one of the tokens in the following table:
 *
 *         -----------------------------------------------------------
 *         center     | The value of the -cameracenter config option.
 *         up         | Vector pointing in the "up" direction of the viewport.
 *         down       | Additive inverse of 'up'.
 *         left       | Vector pointing in the "left" direction of the viewport.
 *         right      | Additive inverse of 'left'.
 *         location   | The calue of the -cameralocation config option.
 *         lineofsight| The vector between 'location' and 'center'.
 *         los        | Alias for 'lineofsight'.
 *
 *     Preceding any of the above symbols with a '-' character returns the
 *     additive inverse. (i.e. "-right" is a synonym for "left".)
 *
 *     Following either form (a) or (b) of the vertex may be a single floating
 *     point number, the multiplier (or W dimension of the vector, if you
 *     prefer). If present, each vertex coordinate is multiplied by the
 *     multiplier before returning. This allows transforms such as:
 *
 *         "move left 1.2"
 *
 * Results:
 *     TCL result. If TCL_ERROR is returned, then an error message is
 *     written to pCanvas->interp.
 *
 * Side effects:
 *     Modifies *piCurrent if successful.
 *
 *---------------------------------------------------------------------------
 */
static int
readVertex(
    Tcl_Interp *interp,
    C3dWidget *pCanvas,
    int *piCurrent,
    Tcl_Obj **apElem,
    int nElem,
    double *pX,
    double *pY,
    double *pZ
) {
    static const char *zMagic[] = {
        "center", "-center",
        "up", "-up", "down", "-down",
        "left", "-left", "right", "-right",
        "location", "-location",
        "lineofsight", "-lineofsight",
        "los", "-los",
        0
    };
    enum MagicVectors {
        CENTER, MINUS_CENTER,
        UP, MINUS_UP, DOWN, MINUS_DOWN,
        LEFT, MINUS_LEFT, RIGHT, MINUS_RIGHT,
        LOCATION, MINUS_LOCATION,
        LINEOFSIGHT, MINUS_LINEOFSIGHT,
        LOS, MINUS_LOS,
    };
    int i = *piCurrent;
    int c;
    double mul;

    if (i >= nElem) {
        return TCL_ERROR;
    }

    if (pCanvas) {
        for (c = 0; zMagic[c]; c++) {
            if (0 == strcasecmp(zMagic[c], Tcl_GetString(apElem[i]))) {
                break;
            }
        }
    }

    if (pCanvas && zMagic[c]) {
        enum MagicVectors choice = (enum MagicVectors)c;
        int up = 0;
        switch (choice) {
            /* "center" is easy, just copy the values from the
             * -cameracenter configuration option.
             */
            case CENTER:
            case MINUS_CENTER:
                *pX = pCanvas->options.cameracenter->x;
                *pY = pCanvas->options.cameracenter->y;
                *pZ = pCanvas->options.cameracenter->z;
                break;

            /* "location" is easy, just copy the values from the
             * -cameralocation configuration option.
             */
            case LOCATION:
            case MINUS_LOCATION:
                *pX = pCanvas->options.cameralocation->x;
                *pY = pCanvas->options.cameralocation->y;
                *pZ = pCanvas->options.cameralocation->z;
                break;

            case LINEOFSIGHT:
            case MINUS_LINEOFSIGHT:
            case MINUS_LOS:
            case LOS:
                *pX = pCanvas->options.cameralocation->x -
                      pCanvas->options.cameracenter->x;
                *pY = pCanvas->options.cameralocation->y -
                      pCanvas->options.cameracenter->y;
                *pZ = pCanvas->options.cameralocation->z -
                      pCanvas->options.cameracenter->z;
                break;

            case UP:
            case MINUS_UP:
            case DOWN:
            case MINUS_DOWN:
                up = 1;

            /* The "left" vector points to the left of the viewport. It is
             * calculated as the cross product of the "up" vector and the
             * vector between the camera location and center.
             */
            case LEFT:
            case MINUS_LEFT:
            case RIGHT:
            case MINUS_RIGHT: {
                float x1 = (pCanvas->options.cameralocation->x -
                            pCanvas->options.cameracenter->x);
                float y1 = (pCanvas->options.cameralocation->y -
                            pCanvas->options.cameracenter->y);
                float z1 = (pCanvas->options.cameralocation->z -
                            pCanvas->options.cameracenter->z);
                float x2 = pCanvas->options.cameraup->x;
                float y2 = pCanvas->options.cameraup->y;
                float z2 = pCanvas->options.cameraup->z;
                float x, y, z;
                C3dCross(x1, y1, z1, x2, y2, z2, &x, &y, &z);
                if (up) {
                    C3dCross(x, y, z, x1, y1, z1, &x, &y, &z);
                }
                C3dNormalize(&x, &y, &z, 0);
                *pX = x;
                *pY = y;
                *pZ = z;
                break;
            }
        }
        if (choice == MINUS_CENTER ||
            choice == DOWN || choice == MINUS_UP ||
            choice == RIGHT || choice == MINUS_LEFT ||
            choice == MINUS_LOCATION ||
            choice == MINUS_LINEOFSIGHT || choice == MINUS_LOS
        ) {
            *pX = *pX * -1.0;
            *pY = *pY * -1.0;
            *pZ = *pZ * -1.0;
        }
        i++;
    } else {
        if (i + 3 > nElem) {
            return TCL_ERROR;
        }
        if (Tcl_GetDoubleFromObj(interp, apElem[i], pX) ||
            Tcl_GetDoubleFromObj(interp, apElem[i+1], pY) ||
            Tcl_GetDoubleFromObj(interp, apElem[i+2], pZ)
        ) {
            return TCL_ERROR;
        }
        i += 3;
    }

    /* Now check if there is a trailing multiplier (a <w> field). */
    if (i < nElem && TCL_OK == Tcl_GetDoubleFromObj(0, apElem[i], &mul)) {
        i++;
        *pX = *pX * mul;
        *pY = *pY * mul;
        *pZ = *pZ * mul;
    }

    *piCurrent = i;
    return TCL_OK;
}

/*
 *---------------------------------------------------------------------------
 *
 * transformCombine --
 *
 *     Combine two transformation matrices into a single matrix.  If
 *     transformation M2 is performed on point v immediately after M1:
 *
 *         v' = M1v             and          v'' = M2v'
 *
 *         v'' = M2(M1v)
 *         v'' = (M2 M1)v
 *
 *     (since matrix multiplication is associative).
 *
 *     Assuming arguments pM1 and pM2 point to matrices M1 and M2
 *     respectively, this function writes the product (M2 M1) into the
 *     space pointed to by argument pOut.
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     Sets all values in *pOut.
 *
 *---------------------------------------------------------------------------
 */
static void
transformCombine(
    C3dTransform *pM1,
    C3dTransform *pM2,
    C3dTransform *pOut
) {
    assert(pM1!=pOut && pM2!=pOut);

    pOut->v1.x =
            pM2->v1.x * pM1->v1.x + pM2->v1.y * pM1->v2.x +
            pM2->v1.z * pM1->v3.x + pM2->v1.w * pM1->v4.x;
    pOut->v2.x =
            pM2->v2.x * pM1->v1.x + pM2->v2.y * pM1->v2.x +
            pM2->v2.z * pM1->v3.x + pM2->v2.w * pM1->v4.x;
    pOut->v3.x =
            pM2->v3.x * pM1->v1.x + pM2->v3.y * pM1->v2.x +
            pM2->v3.z * pM1->v3.x + pM2->v3.w * pM1->v4.x;
    pOut->v4.x =
            pM2->v4.x * pM1->v1.x + pM2->v4.y * pM1->v2.x +
            pM2->v4.z * pM1->v3.x + pM2->v4.w * pM1->v4.x;

    pOut->v1.y =
            pM2->v1.x * pM1->v1.y + pM2->v1.y * pM1->v2.y +
            pM2->v1.z * pM1->v3.y + pM2->v1.w * pM1->v4.y;
    pOut->v2.y =
            pM2->v2.x * pM1->v1.y + pM2->v2.y * pM1->v2.y +
            pM2->v2.z * pM1->v3.y + pM2->v2.w * pM1->v4.y;
    pOut->v3.y =
            pM2->v3.x * pM1->v1.y + pM2->v3.y * pM1->v2.y +
            pM2->v3.z * pM1->v3.y + pM2->v3.w * pM1->v4.y;
    pOut->v4.y =
            pM2->v4.x * pM1->v1.y + pM2->v4.y * pM1->v2.y +
            pM2->v4.z * pM1->v3.y + pM2->v4.w * pM1->v4.y;

    pOut->v1.z =
            pM2->v1.x * pM1->v1.z + pM2->v1.y * pM1->v2.z +
            pM2->v1.z * pM1->v3.z + pM2->v1.w * pM1->v4.z;
    pOut->v2.z =
            pM2->v2.x * pM1->v1.z + pM2->v2.y * pM1->v2.z +
            pM2->v2.z * pM1->v3.z + pM2->v2.w * pM1->v4.z;
    pOut->v3.z =
            pM2->v3.x * pM1->v1.z + pM2->v3.y * pM1->v2.z +
            pM2->v3.z * pM1->v3.z + pM2->v3.w * pM1->v4.z;
    pOut->v4.z =
            pM2->v4.x * pM1->v1.z + pM2->v4.y * pM1->v2.z +
            pM2->v4.z * pM1->v3.z + pM2->v4.w * pM1->v4.z;

    pOut->v1.w =
            pM2->v1.x * pM1->v1.w + pM2->v1.y * pM1->v2.w +
            pM2->v1.z * pM1->v3.w + pM2->v1.w * pM1->v4.w;
    pOut->v2.w =
            pM2->v2.x * pM1->v1.w + pM2->v2.y * pM1->v2.w +
            pM2->v2.z * pM1->v3.w + pM2->v2.w * pM1->v4.w;
    pOut->v3.w =
            pM2->v3.x * pM1->v1.w + pM2->v3.y * pM1->v2.w +
            pM2->v3.z * pM1->v3.w + pM2->v3.w * pM1->v4.w;
    pOut->v4.w =
            pM2->v4.x * pM1->v1.w + pM2->v4.y * pM1->v2.w +
            pM2->v4.z * pM1->v3.w + pM2->v4.w * pM1->v4.w;
}

/*
 *---------------------------------------------------------------------------
 *
 * transformMove --
 *
 *     Write the matrix corresponding to a translation of (x, y, z) to
 *     *pTransform.
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     Modifies *pTransform.
 *
 *---------------------------------------------------------------------------
 */
static void
transformMove(
    double x,
    double y,
    double z,
    C3dTransform *pTransform
) {
    memset(pTransform, 0, sizeof(C3dTransform));
    pTransform->v1.x = 1.0;
    pTransform->v1.w = x;
    pTransform->v2.y = 1.0;
    pTransform->v2.w = y;
    pTransform->v3.z = 1.0;
    pTransform->v3.w = z;
    pTransform->v4.w = 1.0;
}

/*
 *---------------------------------------------------------------------------
 *
 * transformRotate --
 *
 *     Write the matrix corresponding to a rotation of <angle> degrees
 *     around the origin point in the plane perpendicular to the vector
 *     (x, y, z) to *pTransform. (x, y, z) is assumed to be a unit vector.
 *     the rotation is right-handed, so if the unit vector points directly
 *     into the camera objects appear to rotate counter-clockwise.
 *
 *     Assuming (x, y, z) is a unit vector in the direction of the
 *     axis of rotation and:
 *
 *        c = cos(<angle>)   s = sin(<angle>)   t = 1 - cos(<angle>)
 *
 *            / t*x*x+c   t*x*y+sz  t*x*z-sy  0.0 \
 *            | t*y*x-sz  t*y*y+c   t*y*z+sx  0.0 |
 *            | t*z*x+sy  t*z*y-sx  t*z*z+c   0.0 |
 *            \ 0.0       0.0       0.0       1.0 /
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     Modifies *pTransform.
 *
 *---------------------------------------------------------------------------
 */
static void
transformRotate(
    double angle,
    double x,
    double y,
    double z,
    C3dTransform *pTransform
) {
    double s, c, t;
    memset(pTransform, 0, sizeof(C3dTransform));

    angle = (angle / 180.0) * M_PI * -1.0;
    s = sin(angle);
    c = cos(angle);
    t = 1.0 - c;

    pTransform->v1.x = t*x*x+c;
    pTransform->v1.y = t*x*y+s*z;
    pTransform->v1.z = t*x*z-s*y;

    pTransform->v2.x = t*y*x-s*z;
    pTransform->v2.y = t*y*y+c;
    pTransform->v2.z = t*y*z+s*x;

    pTransform->v3.x = t*z*x+s*y;
    pTransform->v3.y = t*z*y-s*x;
    pTransform->v3.z = t*z*z+c;

    pTransform->v4.w = 1.0;
}

/*
 *---------------------------------------------------------------------------
 *
 * transformScale --
 *
 *      Write the matrix corresponding to scaling the X, Y and Z components
 *      of vertices by the coefficients x, y and z passed as arguments to
 *      *pTransform.
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     Modifies *pTransform.
 *
 *---------------------------------------------------------------------------
 */
static void
transformScale(
    double x,
    double y,
    double z,
    C3dTransform *pTransform
) {
    memset(pTransform, 0, sizeof(C3dTransform));
    pTransform->v1.x = x;
    pTransform->v2.y = y;
    pTransform->v3.z = z;
    pTransform->v4.w = 1.0;
}

/*
 *---------------------------------------------------------------------------
 *
 * transformPrimitive --
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     None.
 *
 *---------------------------------------------------------------------------
 */
static int
transformPrimitive(
    ClientData clientData,
    Tcl_Interp *interp,
    C3dWidget *pCanvas,
    Tcl_Obj ** apElem,
    int nElem,
    int *piCurrent,
    C3dTransform *pTransform
) {
    long primitive = (long)clientData;
    int i = *piCurrent;
    double x, y, z, angle;

    i++;
    switch (primitive) {
        case TRANSFORM_MOVE:
            if (readVertex(interp, pCanvas, &i, apElem, nElem, &x, &y, &z)) {
                return TCL_ERROR;
            }
            transformMove(x, y, z, pTransform);
            break;

        case TRANSFORM_SCALE:
            if (readVertex(interp, pCanvas, &i, apElem, nElem, &x, &y, &z)) {
                return TCL_ERROR;
            }
            transformScale(x, y, z, pTransform);
            break;

        case TRANSFORM_ROTATE: {
            double magnitude;

            if (i == nElem) {
                Tcl_ResetResult(interp);
                Tcl_AppendResult(
                        interp, "Expected double but got <end of list>", 0);
                return TCL_OK;
            }
            if (Tcl_GetDoubleFromObj(interp, apElem[i++], &angle) ||
                readVertex(interp, pCanvas, &i, apElem, nElem, &x, &y, &z)
            ) {
                return TCL_ERROR;
            }

            /* Normalize (x, y, z) */
            magnitude = sqrt(x*x + y*y + z*z);
            x = x / magnitude;
            y = y / magnitude;
            z = z / magnitude;

            transformRotate(angle, x, y, z, pTransform);
            break;
        }

        default: assert(!"Can't happen");
    }

    *piCurrent = i;
    return TCL_OK;
}

/*
 *---------------------------------------------------------------------------
 *
 * transformShortcut --
 *
 *     This function implements the shortcut transforms, for example
 *     "orbitup".
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     None.
 *
 *---------------------------------------------------------------------------
 */
static int
transformShortcut(
    ClientData clientData,
    Tcl_Interp *interp,
    C3dWidget *pCanvas,
    Tcl_Obj ** apElem,
    int nElem,
    int *piCurrent,
    C3dTransform *pTransform
) {
    const char *zPrintf = (const char *)clientData;
    int i = *piCurrent;

    C3dTransform *pTmp;
    char zBuf[256];
    Tcl_Obj *pBuf;
    double arg;

    assert(i < nElem);

    /* Each of the shortcut transforms takes a single floating point
     * argument.
     */
    if (i == (nElem - 1)) {
        const char *z = Tcl_GetString(apElem[i]);
        Tcl_ResetResult(interp);
        Tcl_AppendResult(interp,"Transform \"", z, "\" requires an argument",0);
        return TCL_ERROR;
    }
    if (Tcl_GetDoubleFromObj(interp, apElem[i+1], &arg)) {
        return TCL_ERROR;
    }
    i += 2;

    if (zPrintf) {
        sprintf(zBuf, zPrintf, arg);
        pBuf = Tcl_NewStringObj(zBuf, -1);
    } else {
        pBuf = Tcl_NewStringObj("move -center scale ", -1);
        Tcl_AppendObjToObj(pBuf, Tcl_NewDoubleObj(arg));
        Tcl_AppendObjToObj(pBuf, Tcl_NewStringObj(" ", -1));
        Tcl_AppendObjToObj(pBuf, Tcl_NewDoubleObj(arg));
        Tcl_AppendObjToObj(pBuf, Tcl_NewStringObj(" ", -1));
        Tcl_AppendObjToObj(pBuf, Tcl_NewDoubleObj(arg));
        Tcl_AppendObjToObj(pBuf, Tcl_NewStringObj(" move center", -1));
    }
    C3dIncrRefCount(pBuf);
    pTmp = C3dTransformCompile(interp, pCanvas, pBuf);
    C3dDecrRefCount(pBuf);

    if (pTmp) {
        memcpy(pTransform, pTmp, sizeof(C3dTransform));
        C3dTransformDelete(pTmp);
    } else {
        return TCL_ERROR;
    }

    *piCurrent = i;
    return TCL_OK;
}

/*
 *---------------------------------------------------------------------------
 *
 * transformIdentity --
 *
 *     Modify the transformation matrix pointed to by pTransform so that it
 *     contains the identity transformation.
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     None.
 *
 *---------------------------------------------------------------------------
 */
static void
transformIdentity(
    C3dTransform *pTransform
) {
    memset(pTransform, 0, sizeof(C3dTransform));
    pTransform->v1.x = 1.0;
    pTransform->v2.y = 1.0;
    pTransform->v3.z = 1.0;
    pTransform->v4.w = 1.0;
}

/*
 *---------------------------------------------------------------------------
 *
 * transformLookat --
 *
 *     The search expression is used to find a bounding sphere that
 *     encapsulates all the objects to "look at". Logically, the camera is
 *     rotated in place until it points at the center of the sphere, and
 *     then moved away from or toward the center of the sphere so that the
 *     sphere consumes only slightly less than the entire viewport.
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     None.
 *
 *---------------------------------------------------------------------------
 */
static int
transformLookat(
    ClientData clientData,
    Tcl_Interp *interp,
    C3dWidget *pCanvas,
    Tcl_Obj ** apElem,
    int nElem,
    int *piCurrent,
    C3dTransform *pTransform
) {
    float r, x, y, z;             /* Bounding sphere definition */
    int valid;                    /* True if sphere is valid */
    float d;                      /* Distance from sphere center to camera */
    float a;                      /* Visible-angle in radians */
    float cr, cx, cy, cz;         /* Camera location relative to center */
    float nr, nx, ny, nz;
    float rx, ry, rz;
    float rotation;
    float rmag;
    int n;
    long flags = (long)clientData;
    char zTransform[1000];
    Tcl_Obj *pSearch;
    Tcl_Obj *pObj;
    C3dTransform *pTmp;

    C3dVertex *pLoc = pCanvas->options.cameralocation;
    C3dVertex *pCenter = pCanvas->options.cameracenter;

    int i = *piCurrent;

    if (i == (nElem - 1)) {
        const char *z = Tcl_GetString(apElem[i]);
        Tcl_ResetResult(interp);
        Tcl_AppendResult(interp,"Transform \"", z, "\" requires an argument",0);
        return TCL_ERROR;
    }
    pSearch = apElem[i + 1];
    i += 2;
    *piCurrent = i;

    if (C3dBoundingSphere(pCanvas, pSearch, &valid, &r, &x, &y, &z)) {
        return TCL_ERROR;
    }

    /* If the search returned no items, then return the identity matrix * */
    if (!valid) {
        transformIdentity(pTransform);
        return TCL_OK;
    }

    /* Figure out the current visible angle. If the width of the widget
     * viewport is greater than the height, then this is the configured
     * -visibleangle option. Otherwise it is the calculated horizontal
     * visible angle.
     */
    a = pCanvas->options.visibleangle * M_PI / 180.0;
    if (pCanvas->options.width < pCanvas->options.height) {
        float w = pCanvas->options.width;
        float h = pCanvas->options.height;
        a = 2.0 * atan(tan(a/2.0) * w / h);
    }

    /* If A is the current visible-angle, R is the radius of the sphere
     * and D is the required distance between the center of the sphere and
     * the camera, then:
     *
     *     sin(A/2) = R / D
     *     D        = R / sin(A/2)
     */
    d = r / sin(a/2.0);

    /* Set (cx, cy, cz) to the vector from the current location to the
     * current center. (nx, ny, nz) is set to the vector from the current
     * location to the new center.
     *
     * We cross these two vectors to get (rx, ry, rz), the vector normal to
     * the plane defined by the two vectors. If (cx, cy, cz) is parallel to
     * (nx, ny, nz), then the result of the cross will have a magnitude of
     * zero. In this case we find (rx, ry, rz), a vector normal to
     * (cx, cy, cz) by crossing it with the vector (1, 0, 0) or (0, 0, 1).
     */
    cx = pCenter->x - pLoc->x;
    cy = pCenter->y - pLoc->y;
    cz = pCenter->z - pLoc->z;
    nx = x - pLoc->x;
    ny = y - pLoc->y;
    nz = z - pLoc->z;

    C3dNormalize(&cx, &cy, &cz, &cr);
    C3dNormalize(&nx, &ny, &nz, &nr);

    C3dCross(cx, cy, cz, nx, ny, nz, &rx, &ry, &rz);
    rmag = sqrt(rx*rx + ry*ry + rz*rz);
    if (rmag < 0.0001) {
        /* If this happens the the vectors (cx, cy, cz) and (nx, ny, nz)
	 * are very close to parallel. We can approximate the camera
	 * rotation required as either 0 degrees or 180 degrees. We rotate
	 * around the positive Z-axis in this case.
         */
        rx = 0.0;
        ry = 0.0;
        rz = 1.0;
        if ((nx*cx + ny*cy + nz*cz) < 0.0) {
            rotation = 180.0;
        } else {
            rotation = 0.0;
        }
    } else {
        rotation = acos(cx*nx + cy*ny + cz*nz) * 180.0 / M_PI;
    }

    n = 0;
    zTransform[0] = 0;
    if( flags & TRANSFORM_ROTATE ){
        sprintf(zTransform,
            "move %f %f %f "        /* Move camera location to origin. */
            "rotate %f %f %f %f "   /* Rotate camera without moving it. */
            "move %f %f %f ",       /* Move camera back to original location. */

            -1.0 * pLoc->x, -1.0 * pLoc->y, -1.0 * pLoc->z,
            rotation, rx, ry, rz,
            pLoc->x, pLoc->y, pLoc->z
        );
        n += strlen(zTransform);
    }
    if( flags & TRANSFORM_SCALE ){
        sprintf(&zTransform[n],
            "move %f %f %f "      /* Move camera center to origin. */
            "scale %f %f %f "     /* Scale camera location relative to center */
            "move %f %f %f ",     /* Move camera center to center. */

            nx * (nr - cr) - x, ny * (nr - cr) - y, nz * (nr - cr) - z,
            d / cr, d / cr, d / cr,
            x, y, z
        );
    }

    pObj = Tcl_NewStringObj(zTransform, -1);
    C3dIncrRefCount(pObj);
    pTmp = C3dTransformCompile(interp, pCanvas, pObj);
    memcpy(pTransform, pTmp, sizeof(C3dTransform));
    C3dTransformDelete(pTmp);
    C3dDecrRefCount(pObj);

    return TCL_OK;
}

/*
 *---------------------------------------------------------------------------
 *
 * C3dCross --
 *
 *     Compute the cross product of the vectors (x1, y1, z1) and
 *     (x2, y2, z2). Return the resulting vector in *pX, *pY, *pZ.
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     Modifies *pX, *pY and *pZ.
 *
 *---------------------------------------------------------------------------
 */
void
C3dCross(
    float x1,
    float y1,
    float z1,
    float x2,
    float y2,
    float z2,
    float *pX,
    float *pY,
    float *pZ
) {
    *pX = y1*z2 - z1*y2;
    *pY = z1*x2 - x1*z2;
    *pZ = x1*y2 - y1*x2;
}

/*
 *---------------------------------------------------------------------------
 *
 * C3dNormalize --
 *
 *     Normalize the vector (*pX, *pY, *pZ).
 *
 *     If pM is not NULL, write the magnitude of the vector before
 *     normalization to *pM.
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     Modifies *pX, *pY and *pZ, unless the vector is already normalized.
 *
 *---------------------------------------------------------------------------
 */
void
C3dNormalize(
    float *pX,
    float *pY,
    float *pZ,
    float *pM
) {
    float x = *pX;
    float y = *pY;
    float z = *pZ;
    float mag = sqrt(x*x + y*y + z*z);
    float scale = 1.0 / mag;

    *pX = x * scale;
    *pY = y * scale;
    *pZ = z * scale;
    if (pM) {
        *pM = mag;
    }
}

/*
 *---------------------------------------------------------------------------
 *
 * C3dDot --
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     None.
 *
 *---------------------------------------------------------------------------
 */
float
C3dDot(
    float x1,
    float y1,
    float z1,
    float x2,
    float y2,
    float z2
) {
    return (x1*x2 + y1*y2 + z1*z2);
}

/*
 *---------------------------------------------------------------------------
 *
 * C3dMagnitude --
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     None.
 *
 *---------------------------------------------------------------------------
 */
float C3dMagnitude(
    float x,
    float y,
    float z
) {
    return sqrt(x*x + y*y + z*z);
}

/*
 *---------------------------------------------------------------------------
 *
 * C3dTransformCompile --
 *
 *     Compile the transformation specification contained in argument pObj.
 *
 *     If pObj represents a valid specification, a C3dTransform structure
 *     is allocated, populated with the corresponding transformation matrix
 *     and returned. The caller is responsible for eventually freeing the
 *     allocated structure by passing it to C3dTransformDelete(). Vertices
 *     may be transformed by the matrix using C3dTransformApply().
 *
 *     If pObj does not represent a valid transformation, NULL is returned
 *     and an error written to pCanvas->interp.
 *
 * Results:
 *     If successful, pointer to a newly allocated and populated
 *     C3dTransform structure. Otherwise NULL.
 *
 * Side effects:
 *     None.
 *
 *---------------------------------------------------------------------------
 */
C3dTransform * C3dTransformCompile(
    Tcl_Interp *interp,
    C3dWidget *pCanvas,
    Tcl_Obj *pObj
) {
    Tcl_Obj **apElem;
    int nElem;
    int i = 0;
    C3dTransform *pTransform;

#define SHORTCUT(x, y) {#x, transformShortcut, y}
    static struct _Transforms {
        const char *zTransform;
        int (*xTransform)(
            ClientData,
            Tcl_Interp *,
            C3dWidget *,
            Tcl_Obj **,
            int,
            int *,
            C3dTransform *
        );
        ClientData clientData;
    } transforms[] = {
        /* The three basic transformations: move, rotate and scale */
        {"move",       transformPrimitive, (ClientData)TRANSFORM_MOVE   },
        {"rotate",     transformPrimitive, (ClientData)TRANSFORM_ROTATE },
        {"scale",      transformPrimitive, (ClientData)TRANSFORM_SCALE  },

        /* Shortcut transformations designed for camera manipulation */
        SHORTCUT(orbitup,    "move -center rotate %f left  move center"     ),
        SHORTCUT(orbitdown,  "move -center rotate %f right move center"     ),
        SHORTCUT(orbitright, "move -center rotate %f up    move center"     ),
        SHORTCUT(orbitleft,  "move -center rotate %f down  move center"     ),
        SHORTCUT(panup,      "move -location rotate %f right move location" ),
        SHORTCUT(pandown,    "move -location rotate %f left  move location" ),
        SHORTCUT(panright,   "move -location rotate %f down  move location" ),
        SHORTCUT(panleft,    "move -location rotate %f up    move location" ),
        SHORTCUT(movein,     0),
        SHORTCUT(twistleft,  "move -location rotate %f los move location"   ),
        SHORTCUT(twistright, "move -location rotate %f -los move location"  ),

        /* The magical lookat transformation. */
        {"lookat",     transformLookat, (ClientData)(TRANSFORM_ROTATE|TRANSFORM_SCALE) },
        {"looktoward", transformLookat, (ClientData)(TRANSFORM_ROTATE) },
        {"zoomto",     transformLookat, (ClientData)(TRANSFORM_SCALE) },
        {0, 0, 0}
    };
    if (Tcl_ListObjGetElements(interp, pObj, &nElem, &apElem)) {
        return 0;
    }

    /* Initialise *pTransform to the identity matrix. */
    pTransform = (C3dTransform *)C3dAlloc(sizeof(C3dTransform));
    memset(pTransform, 0, sizeof(C3dTransform));
    pTransform->v1.x = 1.0;
    pTransform->v2.y = 1.0;
    pTransform->v3.z = 1.0;
    pTransform->v4.w = 1.0;

    /* This loop runs until all the tokens in the transform specification
     * have been processed. If a syntax error is encountered, an error is
     * left in pCanvas->interp and 0 is returned.
     */
    while (i < nElem) {
        C3dTransform sTransform;
        C3dTransform sCurrent;
        int j;
        const char *zOp = Tcl_GetString(apElem[i]);

        memset(&sTransform, 0, sizeof(C3dTransform));

        for (j = 0; transforms[j].zTransform; j++) {
            if (0 == strcmp(transforms[j].zTransform, zOp)) {
                break;
            }
        }

        if (!transforms[j].zTransform) {
            goto no_such_transform;
        }

        if (transforms[j].xTransform(
                transforms[j].clientData, interp, pCanvas,
                apElem, nElem, &i, &sTransform)
        ) {
            return 0;
        }

        memcpy(&sCurrent, pTransform, sizeof(C3dTransform));
        transformCombine(&sCurrent, &sTransform, pTransform);
    }

    return pTransform;

no_such_transform: {
        const char *zOp = Tcl_GetString(apElem[i]);
        int j;

        Tcl_ResetResult(interp);
        Tcl_AppendResult(interp, "Bad transform \"", zOp, "\" : must be ", 0);
        for (j = 0; transforms[j].zTransform; j++) {
             if (transforms[j+1].zTransform) {
                 Tcl_AppendResult(interp, transforms[j].zTransform, ", ", 0);
             } else {
                 Tcl_AppendResult(interp, "or ", transforms[j].zTransform, 0);
             }
        }
    }
    return 0;
}

/*
 *---------------------------------------------------------------------------
 *
 * C3dTransformApply --
 *
 *     Apply the transformation pTransform to vertex pVertex. See comments
 *     above the C3dTransform structure in this file for a description of
 *     the arithmatic this entails.
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     Modifies the contents of *pVertex.
 *
 *---------------------------------------------------------------------------
 */
void
C3dTransformApply(
    C3dTransform *pTransform,
    C3dVertex *pVertex
) {
    float x;
    float y;
    float z;
    float w;

    x = pVertex->x;
    y = pVertex->y;
    z = pVertex->z;

    pVertex->x =
            x * pTransform->v1.x + y * pTransform->v1.y +
            z * pTransform->v1.z + pTransform->v1.w;
    pVertex->y =
            x * pTransform->v2.x + y * pTransform->v2.y +
            z * pTransform->v2.z + pTransform->v2.w;
    pVertex->z =
            x * pTransform->v3.x + y * pTransform->v3.y +
            z * pTransform->v3.z + pTransform->v3.w;
    w =
            x * pTransform->v4.x + y * pTransform->v4.y +
            z * pTransform->v4.z + pTransform->v4.w;

    pVertex->x = pVertex->x / w;
    pVertex->y = pVertex->y / w;
    pVertex->z = pVertex->z / w;
}

/*
 *---------------------------------------------------------------------------
 *
 * C3dTransformDelete --
 *
 *     Delete a transform previously allocated by C3dTransformCompile. The
 *     pointer pTransform is unusable after this call.
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     Frees memory associated with matrix pTransform.
 *
 *---------------------------------------------------------------------------
 */
void
C3dTransformDelete(
  C3dTransform *pTransform
) {
    if (pTransform) {
        C3dFree((char *)pTransform);
    }
}

/*
 *---------------------------------------------------------------------------
 *
 * projectFrustum --
 *
 *     Set *pTransform to the same matrix that would be produced by the
 *     OpenGL function glFrustum(), using the current viewing angle and
 *     clipping planes.
 *
 * Results:
 *     TCL result.
 *
 * Side effects:
 *     Modifies *pTransform.
 *
 *---------------------------------------------------------------------------
 */
static int
projectFrustum(
    C3dWidget *pCanvas,
    C3dTransform *pTransform
) {
    /* The matrix produced by this function is identical to that produced
     * by the OpenGL function:
     *
     *     glFrustum(l, r, b, t, n, f)
     */
    float n;                     /* Near   */
    float f;                     /* Far    */
    float l;                     /* Left   */
    float r;                     /* Right  */
    float b;                     /* Bottom */
    float t;                     /* Top    */

    float focus;

    Tk_Window win = pCanvas->tkwin;
    float aspect = (double)Tk_Width(win)/(double)Tk_Height(win);

    if (Tk_Width(win) == 0) {
        return TCL_ERROR;
    }

    focus = C3dMagnitude(
        pCanvas->options.cameracenter->x - pCanvas->options.cameralocation->x,
        pCanvas->options.cameracenter->y - pCanvas->options.cameralocation->y,
        pCanvas->options.cameracenter->z - pCanvas->options.cameralocation->z
    );
    n = focus / 30.0;
    f = focus * 30.0;

    t = n * tan(pCanvas->options.visibleangle * M_PI / 360.0);
    b = -1.0 * t;
    r = t * aspect;
    l = b * aspect;

    memset(pTransform, 0, sizeof(C3dTransform));
    pTransform->v1.x = 2.0 * n / (r - l);
    pTransform->v1.z = (r + l) / (r - l);
    pTransform->v2.y = -2.0 * n / (t - b);
    pTransform->v2.z = (t + b) / (t - b);
    pTransform->v3.z = -1.0 * (f + n) / (f - n);
    pTransform->v3.w = -2.0 * f * n / (f - n);
    pTransform->v4.z = -1.0;

    return TCL_OK;
}

/*
 *---------------------------------------------------------------------------
 *
 * projectLookat --
 *
 *     Set *pTransform to the same matrix that would be produced by the
 *     OpenGL function gluLookAt(), using the current camera configuration.
 *
 * Results:
 *     TCL result.
 *
 * Side effects:
 *     Modifies *pTransform.
 *
 *---------------------------------------------------------------------------
 */
static int
projectLookat(
    C3dWidget *pCanvas,
    C3dTransform *pTransform
) {
    C3dTransform sTransform;
    C3dTransform sTranslate;
    float leftx, lefty, leftz;
    float fx, fy, fz;
    float upx, upy, upz;

    memset(&sTransform, 0, sizeof(C3dTransform));
    memset(&sTranslate, 0, sizeof(C3dTransform));

    /* Variables (fx, fy, fz) store the unit vector pointing from the
     * center point to the camera location.
     */
    fx = pCanvas->options.cameracenter->x - pCanvas->options.cameralocation->x;
    fy = pCanvas->options.cameracenter->y - pCanvas->options.cameralocation->y;
    fz = pCanvas->options.cameracenter->z - pCanvas->options.cameralocation->z;
    C3dNormalize(&fx, &fy, &fz, 0);

    /* (upx, upy, upz) is the "up" vector. */
    upx = pCanvas->options.cameraup->x;
    upy = pCanvas->options.cameraup->y;
    upz = pCanvas->options.cameraup->z;

    /* (leftx, lefty, leftz) is a unit vector pointing left across the
     * viewport. Obtained by (F x U), where F is the "forward" vector and U
     * is the "up" vector.
     */
    C3dCross(fx, fy, fz, upx, upy, upz, &leftx, &lefty, &leftz);
    C3dNormalize(&leftx, &lefty, &leftz, 0);

    /* The "up" vector we have at the moment may be supplied by the user,
     * and is only required to be an approximation. Compute an exact "up"
     * vector as (L x F). Because both L and F are normalized, we don't
     * need to normalize the result of this cross.
     */
    C3dCross(leftx, lefty, leftz, fx, fy, fz, &upx, &upy, &upz);

    sTransform.v1.x = leftx;
    sTransform.v1.y = lefty;
    sTransform.v1.z = leftz;
    sTransform.v2.x = upx;
    sTransform.v2.y = upy;
    sTransform.v2.z = upz;
    sTransform.v3.x = -1.0 * fx;
    sTransform.v3.y = -1.0 * fy;
    sTransform.v3.z = -1.0 * fz;
    sTransform.v4.w = 1.0;

    transformMove(
        -1.0 * pCanvas->options.cameralocation->x,
        -1.0 * pCanvas->options.cameralocation->y,
        -1.0 * pCanvas->options.cameralocation->z,
        &sTranslate
    );
    transformCombine(&sTranslate, &sTransform, pTransform);

    return TCL_OK;
}

/*
 *---------------------------------------------------------------------------
 *
 * transformSelect --
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     None.
 *
 *---------------------------------------------------------------------------
 */
static void
transformSelect(
    C3dWidget *pCanvas,
    C3dTransform *pTransform
) {

    if (pCanvas->eSelectMode) {
        /* If we're in select mode, then we zoom in on a small rectangle
         * around the point selected. Anything that is drawn in this small
         * area will be considered selected.
         */
        C3dTransform sMove;
        C3dTransform sScale;
        C3dTransform sTmp;

        float x = (pCanvas->fSelectX1 + pCanvas->fSelectX2) / 2;
        float y = (pCanvas->fSelectY1 + pCanvas->fSelectY2) / 2;

        float xtol = (pCanvas->fSelectX2 - pCanvas->fSelectX1) / 2;
        float ytol = (pCanvas->fSelectY2 - pCanvas->fSelectY1) / 2;

        transformMove(x * -1.0,y *-1.0, 0.0, &sMove);
        transformScale(1.0 / xtol, 1.0 / ytol, 1.0, &sScale);
        transformCombine(pTransform, &sMove, &sTmp);
        transformCombine(&sTmp, &sScale, pTransform);
    }
}

/*
 *---------------------------------------------------------------------------
 *
 * C3dTransformProjection --
 *
 *     Return a transformation that mimics the projection of a vertex onto
 *     the viewport, according to the current camera configuration. The
 *     returned transformation should be eventually freed by the caller,
 *     using C3dTransformDelete().
 *
 *     Applying the returned transformation to a vertex in 3-D space
 *     results in a vertex with all three values in the range [-1.0, +1.0].
 *     The X and Y elements of the result vertex are the location of the
 *     projected vertex on the viewport, in a cartesian coordinate system
 *     with (-1.0, -1.0) at the top-left of the window and (+1.0, +1.0) at
 *     the bottom-right.
 *
 *     The Z element of the result vertex is the distance of the projected
 *     vertex from the camera location. A value of +1.0 represents an
 *     element right on the near clipping plane, a value of -1.0 is an
 *     object on the far clipping plane.
 *
 * Results:
 *     A compiled transform object.
 *
 * Side effects:
 *     None.
 *
 *---------------------------------------------------------------------------
 */
C3dTransform *
C3dTransformProjection(
    C3dWidget *pCanvas
) {
    C3dTransform *pTransform;
    C3dTransform sFrustum;
    C3dTransform sLookat;

    pTransform = (C3dTransform *)C3dAlloc(sizeof(C3dTransform));

    if (projectFrustum(pCanvas, &sFrustum) ||
        projectLookat(pCanvas, &sLookat)
    ) {
        return 0;
    }
    transformCombine(&sLookat, &sFrustum, pTransform);
    transformSelect(pCanvas, pTransform);

    return pTransform;
}

/*
 *---------------------------------------------------------------------------
 *
 * C3dTransformOverview --
 *
 *     Return the transformation matrix used to project overlay items. This
 *     function does for overlay items what C3dTransformProjection does for
 *     three dimensional items.
 *
 *     The value returned is usually used as a projection matrix.
 *
 *     In normal mode (when C3dWidget.eSelectMode is false), this function
 *     returns a matrix that makes OpenGL coordinates match screen
 *     coordinates. (i.e. the OpenGL vertex (0.0, 0.0, 0.0) corresponds to
 *     the top left of the viewport, and (<width>-1, <height>-1, 0.0) is
 *     the bottom right. Everything is orthogonally projected. The near and
 *     far clipping planes are defined by (z = -1) and (z = 1),
 *     respectively.
 *
 *     In Select mode, the matrix returned projects only the small
 *     rectangle around the selection point into the viewport region.
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     None.
 *
 *---------------------------------------------------------------------------
 */
C3dTransform *
C3dTransformOverview(
    C3dWidget *pCanvas
) {
    /* The matrix produced by this function is identical to that produced
     * by the OpenGL function:
     *
     *     glOrtho(l, r, b, t, n, f)
     */
    float n;                     /* Near   */
    float f;                     /* Far    */
    float l;                     /* Left   */
    float r;                     /* Right  */
    float b;                     /* Bottom */
    float t;                     /* Top    */
    C3dTransform *pTransform;

    pTransform = (C3dTransform *)C3dAlloc(sizeof(C3dTransform));
    memset(pTransform, 0, sizeof(C3dTransform));
    n = -1.0;
    f = 1.0;
    l = 0.0;
    r = (float)pCanvas->options.width;
    b = 0.0;
    t = (float)pCanvas->options.height;

    pTransform->v1.x = 2.0 / (r - l);
    pTransform->v1.w = -(r + l) / (r - l);
    pTransform->v2.y = 2.0 / (t - b);
    pTransform->v2.w = -(t + b) / (t - b);
    pTransform->v3.z = -2.0 / (f - n);
    pTransform->v3.w = -(f + n) / (f - n);
    pTransform->v4.w = 1.0;

    transformSelect(pCanvas, pTransform);

    return pTransform;
}

/*
 *---------------------------------------------------------------------------
 *
 * C3dTransformToMatrix --
 *
 *     Copy the transformation pTransform to pMatrix, which must be a
 *     pointer to an array of type GLdouble, size 16. The array may then be
 *     used as an OpenGL matrix with commands like glMultMatrix().
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     None.
 *
 *---------------------------------------------------------------------------
 */
void
C3dTransformToMatrix(
    C3dTransform *pTransform,
    GLdouble *pMatrix
) {
    #define M(row,col)  pMatrix[col*4+row]
    M(0, 0) = pTransform->v1.x;
    M(0, 1) = pTransform->v1.y;
    M(0, 2) = pTransform->v1.z;
    M(0, 3) = pTransform->v1.w;

    M(1, 0) = -1.0 * pTransform->v2.x;
    M(1, 1) = -1.0 * pTransform->v2.y;
    M(1, 2) = -1.0 * pTransform->v2.z;
    M(1, 3) = -1.0 * pTransform->v2.w;

    M(2, 0) = pTransform->v3.x;
    M(2, 1) = pTransform->v3.y;
    M(2, 2) = pTransform->v3.z;
    M(2, 3) = pTransform->v3.w;

    M(3, 0) = pTransform->v4.x;
    M(3, 1) = pTransform->v4.y;
    M(3, 2) = pTransform->v4.z;
    M(3, 3) = pTransform->v4.w;
    #undef M
}

/*
 *---------------------------------------------------------------------------
 *
 * printTransform --
 *
 *     This procedure prints a transformation matrix to standard output. It
 *     is not used by the widget code, but was useful at one point to call
 *     from within a debugger.
 *
 * Results:
 *     None.
 *
 * Side effects:
 *     Prints using printf().
 *
 *---------------------------------------------------------------------------
 */
#ifndef NDEBUG
static void
printTransform(
    C3dTransform *pTransform
) {
    printf("%.4f %.4f %.4f %.4f\n",
        pTransform->v1.x, pTransform->v1.y, pTransform->v1.z, pTransform->v1.w
    );
    printf("%.4f %.4f %.4f %.4f\n",
        pTransform->v2.x, pTransform->v2.y, pTransform->v2.z, pTransform->v2.w
    );
    printf("%.4f %.4f %.4f %.4f\n",
        pTransform->v3.x, pTransform->v3.y, pTransform->v3.z, pTransform->v3.w
    );
    printf("%.4f %.4f %.4f %.4f\n",
        pTransform->v4.x, pTransform->v4.y, pTransform->v4.z, pTransform->v4.w
    );
    printf("\n");

    /* Get rid of the warning that printTransform is defined but not used. */
    if (0) { printTransform(pTransform); }
}
static void
printMatrix(m)
    GLdouble *m;
{
    printf("%.4f %.4f %.4f %.4f\n", m[0], m[4], m[8], m[12]);
    printf("%.4f %.4f %.4f %.4f\n", m[1], m[5], m[9], m[13]);
    printf("%.4f %.4f %.4f %.4f\n", m[2], m[6], m[10], m[14]);
    printf("%.4f %.4f %.4f %.4f\n", m[3], m[7], m[11], m[15]);
    printf("\n");

    /* Get rid of the warning that printMatrix is defined but not used. */
    if (0) { printMatrix(m); }
}
#endif
