xref: /dports/graphics/openrm/rm160/rm/rmmatrix.c

/*
 * Copyright (C) 1997-2005, 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
 * Additional Contributor(s):
 *
 * The OpenRM project is located at http://openrm.sourceforge.net/.
 */
/*
 * $Id: rmmatrix.c,v 1.7 2005/07/04 16:49:38 wes Exp $
 * Version: $Name: OpenRM-1-6-0-2-RC2 $
 * $Revision: 1.7 $
 * $Log: rmmatrix.c,v $
 * Revision 1.7  2005/07/04 16:49:38  wes
 * Minor change to private_rmNewWCfromDCOffset.
 *
 * Revision 1.6  2005/06/26 18:52:35  wes
 * Streamlined computation of viewport matrix: created single routine that
 * computes the viewport matrix, invoke that routine where needed rather
 * than having duplicate copies of viewport matrix computation code.
 *
 * Revision 1.5  2005/02/19 16:22:50  wes
 * Distro sync and consolidation.
 *
 * Revision 1.4  2005/01/23 17:04:03  wes
 * Copyright update to 2005.
 *
 * Revision 1.3  2004/01/16 16:46:09  wes
 * Updated copyright line for 2004.
 *
 * Revision 1.2  2003/02/02 02:07:15  wes
 * Updated copyright to 2003.
 *
 * Revision 1.1.1.1  2003/01/28 02:15:23  wes
 * Manual rebuild of rm150 repository.
 *
 * Revision 1.10  2003/01/16 22:21:17  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.9  2002/04/30 19:32:22  wes
 * Updated copyright dates.
 *
 * Revision 1.8  2001/10/15 00:08:40  wes
 * Minor change in private_rmDC2fromWC() needed to display image
 * data. Problem was that the modelview and projection matrices
 * definitions changed slightly in the RMstate structure during
 * the multistage rendering work.
 *
 * Revision 1.7  2001/03/31 17:12:38  wes
 * v1.4.0-alpha-2 checkin.
 *
 * Revision 1.6  2000/12/03 22:35:38  wes
 * Mods for thread safety.
 *
 * Revision 1.5  2000/08/19 16:09:06  wes
 * Fixed logic bug in rmMatrixGetValue().
 *
 * Revision 1.4  2000/07/05 04:59:50  wes
 * Fixed index range check error in rmMatrixGetValue.
 *
 * Revision 1.3  2000/04/20 16:29:47  wes
 * Documentation additions/enhancements, some code rearragement.
 *
 * Revision 1.2  2000/02/29 23:43:53  wes
 * Compile warning cleanups.
 *
 * 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.
 *
 */

#include <rm/rm.h>
#include "rmprivat.h"

/*
 * ----------------------------------------------------
 * @Name rmMatrixNew
 @pstart
 RMmatrix * rmMatrixNew (void)
 @pend

 @astart
 No arguments.
 @aend

 @dstart

 Creates a new RMmatrix matrix, returning a handle to the new object
 to the caller upon success, or NULL upon failure.  Use rmMatrixDelete
 to free an RMmatrix object.

 The RMmatrix object is a 4x4 array of floating point values.
 Applications have two alternative approaches that may be used to set
 the individual values in the RMmatrix: either by using
 rmMatrixSetValue()/rmMatrixGetValue(), or through direct C-struct
 access.

 The RMmatrix object consists of a single 4x4 array of floating point
 values. The array is indexed with A[row][column], so that A[3][0]
 refers to an entry that is in the third row, and zero'th column. That
 particular entry controls translations along the X axis; in RM,
 matrix multiplies occur from the left.

 The general form of RMmatrices is
 <pre>
       |  1  0  0  0  |       |  Sx 0  0  0  |
       |  0  1  0  0  |       |  0  Sy 0  0  |
  T =  |  0  0  1  0  |   S = |  0  0  Sz 0  |
       |  Tx Ty Tz 1  |       |  0  0  0  1  |
 </pre>

 Where T and S are the translation and scale matrices, respectively.

 Access to a specific matrix entity using the C-struct interface would
 be performed as follows: to access the Tx, or X-axis translation in T
 above, Tx = matrix->m[3][0], where "matrix" is a handle to the
 RMmatrix object. The "m" entity is the single matrix struct inside
 the RMmatrix object.

 Access to the same element through the API would be:
 Tx = rmMatrixGetValue(matrix,3,0).

 @dend
 * ----------------------------------------------------
 */
RMmatrix *
rmMatrixNew (void)
{
    RMmatrix *m = (RMmatrix *)malloc(sizeof(RMmatrix));
    rmMatrixIdentity(m);
    return(m);
}


/*
 * ----------------------------------------------------
 * @Name rmMatrixCopy
 @pstart
 RMenum rmMatrixCopy (RMmatrix *dst,
	              const RMmatrix *src)
 @pend

 @astart
 RMmatrix *dst - a handle to the "copy to" RMmatrix (modified).

 const RMmatrix *src - a handle to the "copy from" RMmatrix (input).
 @aend

 @dstart

 Copies the contents of one RMmatrix object to another. Returns
 RM_CHILL upon success, or RM_WHACKED upon failure.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmMatrixCopy (RMmatrix *d,
	      const RMmatrix *s)
{
    if ((RM_ASSERT(d, "rmMatrixCopy() error: the dest RMmatrix pointer is NULL.") == RM_WHACKED) ||
	(RM_ASSERT(s, "rmMatrixCopy() error: the src RMmatrix pointer is NULL.") == RM_WHACKED))
	return (RM_WHACKED);

    memcpy((char *)d, (char *)s, sizeof(RMmatrix));
    return(RM_CHILL);
}


/*
 * ----------------------------------------------------
 * @Name rmMatrixDelete
 @pstart
 RMenum rmMatrixDelete (RMmatrix *toDelete)
 @pend

 @astart
 RMmatrix *toDelete - a handle to an RMmatrix object to delete
    (modified).
 @aend

 @dstart

 Releases resources associated with an RMmatrix object. Returns
 RM_CHILL upon success, or RM_WHACKED upon failure.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmMatrixDelete (RMmatrix *m)
{
    if (RM_ASSERT(m, "rmMatrixDelete() error: the input RMmatrix pointer is NULL.") == RM_WHACKED)
	return(RM_WHACKED);

    free(m);
    return(RM_CHILL);
}


/*
 * ----------------------------------------------------
 * @Name rmMatrixIdentity
 @pstart
 RMenum rmMatrixIdentity (RMmatrix *toModify)
 @pend

 @astart
 RMmatrix *toModify - a handle to an RMmatrix object (modified).
 @aend

 @dstart

 This routine will set the RMmatrix "toModify" to contain the Identity
 matrix. Returns RM_CHILL upon success, or RM_WHACKED upon failure.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmMatrixIdentity (RMmatrix *m)
{
    int i, j;

    if (RM_ASSERT(m, "rmMatrixIdentity() error: the input RMmatrix is NULL.") == RM_WHACKED)
	return(RM_WHACKED);

    for (j = 0; j < 4; j++)
	for (i = 0; i < 4; i++)
	    m->m[i][j] = (i == j) ? 1.0 : 0.0;

    return(RM_CHILL);
}


/*
 * ----------------------------------------------------
 * @Name rmMatrixSetValue
 @pstart
 RMenum rmMatrixSetValue (RMmatrix *toModify,
		          int row,
			  int col,
			  float newValue)
 @pend

 @astart
 RMmatrix *toModify - a handle to an RMmatrix object (modified).

 int row, col - integer values, treated as indices. Must both be in
    the range 0..3, otherwise an error condition results.

 float newValue - a floating point value to place into the RMmatrix
    object (input).
 @aend

 @dstart

 Assigns a floating point value (newValue) to a single entry in the
 4x4 matrix contained in the RMmatrix object, returning RM_CHILL upon
 success or RM_WHACKED upon failure.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmMatrixSetValue (RMmatrix *m,
		  int row,
		  int col,
		  float newValue)
{
    RMenum rstat = RM_CHILL;

    if (RM_ASSERT(m, "rmMatrixSetValue() error: the input RMmatrix is NULL.") == RM_WHACKED)
	return(RM_WHACKED);

    if (!(row < 0 || row > 3 || col < 0 || col > 3))
    {
	m->m[row][col] = newValue;
    }
    else
    {
	rmWarning("rmMatrixSetValue warning: either the row or column input parameters are out of range.");
        rstat = RM_WHACKED;
    }
    return(rstat);
}


/*
 * ----------------------------------------------------
 * @Name rmMatrixGetValue
 @pstart
 float rmMatrixGetValue (const RMmatrix *toQuery,
		         int row,
			 int col)
 @pend


 @astart const RMmatrix *toQuery - a handle to an RMmatrix object to
    query (input).

 int row, col - integer values assumed to be indices into the 4x4
    matrix contained within the RMmatrix object. Must be in the range
    0..3 else an error condition is detected.
 @aend

 @dstart

 Use this routine to obtain the value of a single entry from a 4x4
 matrix within an RMmatrix object. Returns the entry at [row][col]
 from the matrix upon success, otherwise an error message is issued
 and 0.0 is returned.

 @dend
 * ----------------------------------------------------
 */
float
rmMatrixGetValue (const RMmatrix *m,
		  int row,
		  int col)
{
    float rval;

    if (RM_ASSERT(m, "rmMatrixGetValue() error: the input RMmatrix is NULL.") == RM_WHACKED)
	return(0.0F);

    if (!((row < 0) || (row > 3) || (col < 0) || (col > 3)))
    {
	rval = m->m[row][col];
    }
    else
    {
	rval = 0.0F;
	rmWarning("rmMatrixGetValue warning: either the row or column input parameters are out of range.");
    }
    return(rval);
}


/*
 * ----------------------------------------------------
 * @Name rmMatrixTranspose
 @pstart
 RMenum rmMatrixTranspose (const RMmatrix *src,
		           RMmatrix *dst)
 @pend

 @astart
 const RMmatrix *src - a handle to an RMmatrix object (input).

 RMmatrix *dst - a handle to an RMmatrix object (modified).
 @aend

 @dstart

 Transposes the input matrix, copying the results back into the "dst"
 RMmatrix object. Returns RM_CHILL upon success, otherwise RM_WHACKED
 is returned.

 This routine is written such that "src" and "dst" may be the same
 RMmatrix object.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmMatrixTranspose (const RMmatrix *s,
		   RMmatrix *d)
{
    int      i, j;
    RMmatrix t;

    if ((RM_ASSERT(s, "rmMatrixTranspose() error: the input RMmatrix object is NULL") == RM_WHACKED) ||
	(RM_ASSERT(d, "rmMatrixTranspose() error: the dest RMmatrix object is NULL")) == RM_WHACKED)
	return(RM_WHACKED);

    for (j = 0; j < 4; j++)
      for (i = 0; i < 4; i++)
        t.m[j][i] = s->m[i][j];
    memcpy(d, &t, sizeof(RMmatrix));

    return(RM_CHILL);
}


/*
 * ----------------------------------------------------
 * @Name rmMatrixMultiply
 @pstart
 RMenum rmMatrixMultiply (const RMmatrix *srcA,
		          const RMmatrix *srcB,
			  RMmatrix *dst)
 @pend

 @astart
 const RMmatrix *srcA, *srcB - handles to RMmatrix objects (input).

 RMmatrix *dst - a handle to an RMmatrix object (modified).
 @aend

 @dstart

 Computes the matrix product srcA * srcB, placing the result into
 "dst".  Returns RM_CHILL upon success, or RM_WHACKED upon failure.

 Internally, rmMatrixMultiply is written such that "dst" may be the
 same RMmatrix object as either srcA or srcB.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmMatrixMultiply (const RMmatrix *a,
		  const RMmatrix *b,
		  RMmatrix *c)
 /* take product of a and b, put result in c */
{
    int      row, col;
    RMmatrix tmp;

    if ((RM_ASSERT(a, "rmMatrixMultiply() error: the input A RMmatrix is NULL") == RM_WHACKED) ||
	(RM_ASSERT(b, "rmMatrixMultiply() error: the input B RMmatrix is NULL") == RM_WHACKED) ||
	(RM_ASSERT(c, "rmMatrixMultiply() error: the destination RMmatrix is NULL") == RM_WHACKED))
	return(RM_WHACKED);

    for (row = 0; row < 4; row++)
    {
	for (col = 0; col < 4; col++)
	    tmp.m[row][col] = (a->m[row][0] * b->m[0][col]) + (a->m[row][1] * b->m[1][col]) + (a->m[row][2] * b->m[2][col]) + (a->m[row][3] * b->m[3][col]);
    }
    for (row = 0; row < 4; row++)
	for (col = 0; col < 4; col++)
	    c->m[row][col] = tmp.m[row][col];

    return(RM_CHILL);
}


/*
 * ----------------------------------------------------
 * @Name rmMatrixInverse
 @pstart
 RMenum rmMatrixInverse (const RMmatrix *src,
		         RMmatrix *dst)
 @pend

 @astart
 const RMmatrix *src - a handle to an RMmatrix object (input).

 RMmatrix *dst - a handle to an RMmatrix object (result).
 @aend

 @dstart

 Will compute the inverse of the 4x4 matrix "src," placing the result
 into "dst". Returns RM_CHILL upon success, or RM_WHACKED upon
 failure.

 Matrix inversion is performed using Gaussian elimination using a
 customized version of SGEFA from the Linpack distribution. If the
 input matrix is singular, an error is reported with rmWarning(), and
 RM_WHACKED is returned.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmMatrixInverse (const RMmatrix *in,
		 RMmatrix *out)
{
    /* the following is a software interface to some routines from LINPACK which are called to compute a matrix inverse */
    int    lda, n, ipvt[5], info, job;
    int   *ip;
    float  m[16];
    float  work[4*4], det[4][2];
    float *a;

    /* external LINPACK-hack */
    extern void sgefa (float *a, int *lda, int *n, int *ipvt, int *info);
    extern void sgedi (float *a, int *lda, int *n, int *ipvt, float *det, float *work, int *job);

    if ((RM_ASSERT(in, "rmMatrixInverse() error: the src RMmatrix pointer is NULL") == RM_WHACKED) ||
	(RM_ASSERT(out, "rmMatrixInverse() error: the dst RMmatrix pointer is NULL") == RM_WHACKED))
	return(RM_WHACKED);

    memcpy((char *)m, (char *)in, sizeof(RMmatrix));
    job = 1;
    lda = 4;
    n = 4;
    a = m;
    ip = ipvt;

    sgefa(a, &lda, &n, ip, &info);

    if (info == 0)
	sgedi(a, &lda, &n, ip, &(det[0][0]), work, &job);
    else
    {
	rmWarning("rmMatrixInverse - the input matrix is singular. ");
	return(RM_WHACKED);
    }
    memcpy((char *)out, (char *)m, sizeof(RMmatrix));

    return(RM_CHILL);
}


/*
 * ----------------------------------------------------
 * @Name rmPoint4MatrixTransform
 @pstart
 RMenum rmPoint4MatrixTransform (const float *src,
			         const RMmatrix *matrix,
				 float *dst)
 @pend

 @astart
 const float *src - a handle to an array of floats of length 4
    (input).

 const RMmatrix *matrix - a handle to an RMmatrix object.

 float *dst - a handle to an array of floats of length 4 (ouput).
 @aend

 @dstart

 rmPoint4MatrixTransform is one of a family of routines that perform
 vector-matrix multiplication. This routine multiplies the vector
 "src" by the matrix, placing the result into vector "dst." It is
 assumed that both "src" and "dst" are floating point arrays of length
 4.  Returns RM_CHILL upon success, or RM_WHACKED upon failure.

 This routine is written such that "src" and "dst" may be the same
 array from the caller's perspective.

 This routine implements the following operation:
 <pre>

                 |  a  b  c  d  |
                 |  e  f  g  h  |
 [s0 s1 s2 s3] * |  i  j  k  l  | = [d0 d1 d2 d3]
		 |  m  n  o  p  |

 </pre>

 The "src" vector multiplies the matrix from the left (pre
 multiplication) such that each element of "dst" is the dot product of
 "src" with a column from the matrix.

 See also rmPointMatrixTransform.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmPoint4MatrixTransform (const float *v,
			 const RMmatrix *m,
			 float *d)
{
    float t[4];

    if ((RM_ASSERT(v, "rmPoint4MatrixTransform() error: the input src vector is NULL") == RM_WHACKED) ||
	(RM_ASSERT(m, "rmPoint4MatrixTransform() error: the input RMmatrix object is NULL") == RM_WHACKED) ||
	(RM_ASSERT(d, "rmPoint4MatrixTransform() error: the result dst vector is NULL") == RM_WHACKED))
	return(RM_WHACKED);

    t[0] = (v[0] * m->m[0][0]) + (v[1] * m->m[1][0]) + (v[2] * m->m[2][0]) + (v[3] * m->m[3][0]);
    t[1] = (v[0] * m->m[0][1]) + (v[1] * m->m[1][1]) + (v[2] * m->m[2][1]) + (v[3] * m->m[3][1]);
    t[2] = (v[0] * m->m[0][2]) + (v[1] * m->m[1][2]) + (v[2] * m->m[2][2]) + (v[3] * m->m[3][2]);
    t[3] = (v[0] * m->m[0][3]) + (v[1] * m->m[1][3]) + (v[2] * m->m[2][3]) + (v[3] * m->m[3][3]);

    memcpy(d, t, sizeof(float) * 4);

    return(RM_CHILL);
}


/*
 * ----------------------------------------------------
 * @Name rmPointMatrixTransform
 @pstart
 RMenum rmPointMatrixTransform (const RMvertex3D *src,
			        const RMmatrix *matrix,
			        RMvertex3D *dst)
 @pend

 @astart
 const RMvertex3D *src - a handle to an RMvertex3D object (input).

 const RMmatrix *matrix - a handle to an RMmatrix object.

 RMvertex3D *dst - a handle to an RMvertex3D (ouput).
 @aend

 @dstart

 rmPointMatrixTransform is one of a family of routines that perform
 vector-matrix multiplication. This routine multiplies the RMvertex3D
 "src" by the matrix, placing the result into RMvertex3D "dst."

 This routine differs from rmPoint4MatrixTransform in that it is
 assumed that the w-coordinate of both "src" and "dst" is 1.0, and
 that the resulting w-coordinate is not available to the caller upon
 return.

 Returns RM_CHILL upon success, or RM_WHACKED upon failure.

 This routine is written such that "src" and "dst" may be the same
 array from the caller's perspective.

 This routine implements the following operation:
 <pre>

                 |  a  b  c  d  |
                 |  e  f  g  h  |
 [s0 s1 s2 1]  * |  i  j  k  l  | = [d0 d1 d2 n/a]
		 |  m  n  o  p  |

 </pre>

 The "src" vector multiplies the matrix from the left (pre
 multiplication) such that each element of "dst" is the dot product of
 "src" with a column from the matrix.

 See also rmPoint4MatrixTransform.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmPointMatrixTransform (const RMvertex3D *s,
		        const RMmatrix *m,
		        RMvertex3D *d)
{
    RMvertex3D t;

    if ((RM_ASSERT(s, "rmPointMatrixTransform error: the input S RMvertex3D pointer is NULL") == RM_WHACKED) ||
	(RM_ASSERT(m, "rmPointMatrixTransform error: the input RMmatrix pointer is NULL") == RM_WHACKED) ||
	(RM_ASSERT(d, "rmPointMatrixTransform error: the destination RMvertex3D pointer is NULL") == RM_WHACKED))
	return(RM_WHACKED);

    t.x = (s->x * m->m[0][0]) + (s->y * m->m[1][0]) + (s->z * m->m[2][0]) + m->m[3][0];
    t.y = (s->x * m->m[0][1]) + (s->y * m->m[1][1]) + (s->z * m->m[2][1]) + m->m[3][1];
    t.z = (s->x * m->m[0][2]) + (s->y * m->m[1][2]) + (s->z * m->m[2][2]) + m->m[3][2];

    VCOPY(&t, d);

    return(RM_CHILL);
}


/*
 * ----------------------------------------------------
 * @Name rmPrintMatrix
 @pstart
 RMenum rmPrintMatrix (const RMmatrix *toPrint)
 @pend

 @astart
 const RMmatrix *toPrint - a handle to an RMmatrix that will be
    printed (input).
 @aend

 @dstart

 Prints the contents of an RMmatrix object using C-printf's. Returns
 RM_CHILL upon success or RM_WHACKED upon failure.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmPrintMatrix (const RMmatrix *m)
{
    int i, j;

    if (RM_ASSERT(m, "rmPrintMatrix() error: the input RMmatrix object is NULL") == RM_WHACKED)
	return(RM_WHACKED);

    for (j = 0; j < 4; j++)
    {
	for (i = 0; i < 4; i++)
	    printf("\t%g", m->m[j][i]);
	printf("\n");
    }
    return(RM_CHILL);
}


/*
 * ----------------------------------------------------
 * @Name rmPointMin
 @pstart
 RMenum rmPointMin (const float *input,
	            int count,
		    int vdims,
		    int stride,
		    RMvertex3D *minReturn)
 @pend

 @astart
 const float *input - a handle to an array of floating point values
    (input).

 int count - an integer value indicating the number of "vertices" that
    will be processed (input).

 int vdims - an integer value indicating the cardinality, or number of
    dimensions, of each vertex. Should be either 1, 2 or 3 (input).

 int stride - an integer value specifying the stride length in floats
    from one vertex to the next (input).

 RMvertex3D *minReturn - a handle to an RMvertex3D object (modified).
 @aend

 @dstart

 Finds the minimum point, or vertex value, from a flat array of
 floating point values. Upon success, the minimum vertex value is
 copied into the caller-supplied "minReturn," and RM_CHILL is
 returned. Otherwise, RM_WHACKED is returned upon failure.

 The parameter "vdims" defines the number of coordinates in each
 vertex. Use a value of 2 to specify two-dimensional vertices (x and y
 values), or 3 for three-dimensional vertices (x, y and z values).  A
 value of 1 is also acceptable.

 The parameter "count" specifies how many such vertices are present in
 the input array. Note that this routine assumes that "input" is just
 a flat array of floats, and will be logically interpreted as groups
 of vertices, each of which consists of "vdims" floats per vertex. The
 coordinates of each vertex are assumed to be contiguously located.

 The parameter "stride" tells how many floating point values to skip
 from one vertex to the next.  If we assume that "input" consists of a
 flat array of RMvertex3D objects (cast to float * when calling this
 routine), the best way to compute "stride" is:

 <pre>
       stride = sizeof(RMvertex3D)/sizeof(float)
 </pre>

 The reason for "stride" is because some platforms force padding of
 C-structures to the nearest 8-byte boundary. Hence,
 sizeof(RMvertex3D) is not always equal to sizeof(float)*3.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmPointMin (const float *input,
	    int count,
	    int vdims,
	    int stride,
	    RMvertex3D *ret)
{
    register int          i;
    register float        mx, my, mz;
    register float        v1;
    register const float *v;

    mx = my = mz = RM_MAXFLOAT;

    if ((RM_ASSERT(input, "rmPointMin error: the input floating point array is NULL") == RM_WHACKED) ||
	(RM_ASSERT(ret, "rmPointMin error: the return RMvertex3D * is NULL.") == RM_WHACKED))
	return(RM_WHACKED);

    stride = stride / sizeof(float);

    v = input;

    for (i = 0; i < count; i++, v += stride)
    {
        /* x-component */
        v1 = v[0];
	if (v1 < mx)
	    mx = v1;

	if (vdims > 1)
	{
	    /* y-component */
	  v1 = v[1];
	  if (v1 < my)
	    my = v1;
	}
	else
	    my = 0.0;

	if (vdims > 2)
	{
	    /* z-component */
	  v1 = v[2];
	  if (v1 < mz)
	      mz = v1;
	}
	else
	    mz = 0.0;
    }
    ret->x = mx;
    ret->y = my;
    ret->z = mz;

    return(RM_CHILL);
}


/*
 * ----------------------------------------------------
 * @Name rmPointMax
 @pstart
 RMenum rmPointMax (const float *input,
	            int count,
		    int vdims,
		    int stride,
		    RMvertex3D *maxReturn)
 @pend

 @astart
 const float *input - a handle to an array of floating point values
    (input).

 int count - an integer value indicating the number of "vertices" that
    will be processed (input).

 int vdims - an integer value indicating the cardinality, or number of
    dimensions, of each vertex. Should be either 1, 2 or 3 (input).

 int stride - an integer value specifying the stride length in floats
    from one vertex to the next (input).

 RMvertex3D *maxReturn - a handle to an RMvertex3D object (modified).
 @aend

 @dstart

 Finds the maximum point, or vertex value, from a flat array of
 floating point values. Upon success, the maximum vertex value is
 copied into the caller-supplied "maxReturn," and RM_CHILL is
 returned. Otherwise, RM_WHACKED is returned upon failure.

 The parameter "vdims" defines the number of coordinates in each
 vertex. Use a value of 2 to specify two-dimensional vertices (x and y
 values), or 3 for three-dimensional vertices (x, y and z values).  A
 value of 1 is also acceptable.

 The parameter "count" specifies how many such vertices are present in
 the input array. Note that this routine assumes that "input" is just
 a flat array of floats, and will be logically interpreted as groups
 of vertices, each of which consists of "vdims" floats per vertex. The
 coordinates of each vertex are assumed to be contiguously located.

 The parameter "stride" tells how many floating point values to skip
 from one vertex to the next.  If we assume that "input" consists of a
 flat array of RMvertex3D objects (cast to float * when calling this
 routine), the best way to compute "stride" is:

 <pre>
       stride = sizeof(RMvertex3D)/sizeof(float)
 </pre>

 The reason for "stride" is because some platforms force padding of
 C-structures to the nearest 8-byte boundary. Hence,
 sizeof(RMvertex3D) is not always equal to sizeof(float)*3.

 @dend
 * ----------------------------------------------------
*/
RMenum
rmPointMax (const float *input,
	    int count,
	    int vdims,
	    int stride,
	    RMvertex3D *ret)
{
    register int          i;
    register float        mx, my, mz;
    register float        v1;
    register const float *v;

    mx = my = mz = -1. * RM_MAXFLOAT;

    if ((RM_ASSERT(input, "rmPointMax error: the input floating point array is NULL") == RM_WHACKED) ||
	(RM_ASSERT(ret, "rmPointMax error: the return RMvertex3D * is NULL.") == RM_WHACKED))
	return(RM_WHACKED);

    stride = stride/sizeof(float);

    /* assume the RMvertex3D structure has 3 floats, and the RMvertex2d structure has 2 floats */
    v = input;

    for (i = 0; i < count; i++, v += stride)
    {
        /* x-component */
        v1 = v[0];
	if (v1 > mx)
	    mx = v1;

	if (vdims > 1)
	{
	  /* y-component */
	  v1 = v[1];
	  if (v1 > my)
	    my = v1;
	}
	else
	    my = 0.0;

	if (vdims > 2)
	{
	  /* z-component */
	  v1 = v[2];
	  if (v1 > mz)
	    mz = v1;
	}
	else
	    mz = 0.0;
    }
    ret->x = mx;
    ret->y = my;
    ret->z = mz;

    return(RM_CHILL);
}


/*
 * ----------------------------------------------------
 * @Name rmPointMinMax
 @pstart
 RMenum rmPointMinMax (const float *input,
	               int count,
		       int vdims,
		       int stride,
		       RMvertex3D *minReturn,
		       RMvertex3D *maxReturn)
 @pend


 @astart
 const float *input - a handle to an array of floating point values
    (input).

 int count - an integer value indicating the number of "vertices" that
    will be processed (input).

 int vdims - an integer value indicating the cardinality, or number of
    dimensions, of each vertex. Should be either 1, 2 or 3 (input).

 int stride - an integer value specifying the stride length in floats
    from one vertex to the next (input).

 RMvertex3D *minReturn, *maxReturn - handles to RMvertex3D objects
    (modified).
 @aend

 @dstart

 Finds the minimum and maximum points, or vertex values, from a flat
 array of floating point values. Upon success, the minimum and maximum
 vertex values are copied into the caller-supplied "minReturn" and
 "maxReturn," respectively, and RM_CHILL is returned. Otherwise,
 RM_WHACKED is returned upon failure.

 The parameter "vdims" defines the number of coordinates in each
 vertex. Use a value of 2 to specify two-dimensional vertices (x and y
 values), or 3 for three-dimensional vertices (x, y and z values).  A
 value of 1 is also acceptable.

 The parameter "count" specifies how many such vertices are present in
 the input array. Note that this routine assumes that "input" is just
 a flat array of floats, and will be logically interpreted as groups
 of vertices, each of which consists of "vdims" floats per vertex. The
 coordinates of each vertex are assumed to be contiguously located.

 The parameter "stride" tells how many floating point values to skip
 from one vertex to the next.  If we assume that "input" consists of a
 flat array of RMvertex3D objects (cast to float * when calling this
 routine), the best way to compute "stride" is:

 <pre>
       stride = sizeof(RMvertex3D)/sizeof(float)
 </pre>

 The reason for "stride" is because some platforms force padding of
 C-structures to the nearest 8-byte boundary. Hence,
 sizeof(RMvertex3D) is not always equal to sizeof(float)*3.

 This routine internally calls rmPointMin, followed by rmPointMax.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmPointMinMax (const float *input,
	       int count,
	       int vdims,
	       int stride,
	        RMvertex3D *min,
	      RMvertex3D *max)
{
    if ((rmPointMin(input, count, vdims, stride, min) != RM_WHACKED) &&	(rmPointMax(input, count, vdims, stride, max) != RM_WHACKED))
	return(RM_CHILL);

    return(RM_WHACKED);
}


/*
 * ----------------------------------------------------
 * @Name rmVertex3DMag
 @pstart
 double rmVertex3DMag (const RMvertex3D *src)
 @pend

 @astart
 const RMvertex3D *src - a handle to an RMvertex3D object (input).
 @aend

 @dstart

 Computes the vector magnitude of the input RMvertex3D object (a
 3-component vector), returning the result to the caller. If the input
 is NULL, an error message is reported, and -1.0 is returned.

 @dend
 * ----------------------------------------------------
 */
double
rmVertex3DMag (const RMvertex3D *v)
{
    double d;

    if (RM_ASSERT(v, "rmVertex3DMag() error: the input RMvertex3D object is NULL") == RM_WHACKED)
	return(-1.0);

    d = (v->x * v->x) + (v->y * v->y) + (v->z * v->z);

    return(sqrt(d));
}


/*
 * ----------------------------------------------------
 * @Name rmVertex3DSum
 @pstart
 RMenum rmVertex3DSum (const RMvertex3D *a,
	               const RMvertex3D *b,
		       RMvertex3D *dst)
 @pend

 @astart
 const RMvertex3D *a, *b - handles to RMvertex3D objects (input).

 RMvertex3D *dst - a handle to an RMvertex3D object (result).
 @aend

 @dstart

 Computes the vector sum (a + b), placing the result into "dst".  This
 is computed such that either a or b may be used as the result c.
 Returns RM_CHILL upon success, or RM_WHACKED upon failure.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmVertex3DSum (const RMvertex3D *a,
	       const RMvertex3D *b,
	       RMvertex3D *c)
{
    if ((RM_ASSERT(a, "rmVertex3DSum() error: the input RMvertex3D object A is NULL") == RM_WHACKED) ||
	(RM_ASSERT(b, "rmVertex3DSum() error: the input RMvertex3D object B is NULL") == RM_WHACKED) ||
	(RM_ASSERT(c, "rmVertex3DSum() error: the dest RMvertex3D object is NULL") == RM_WHACKED))
	return(RM_WHACKED);

    /* compute c = a + b */
    c->x = a->x + b->x;
    c->y = a->y + b->y;
    c->z = a->z + b->z;

    return(RM_CHILL);
}


/*
 * ----------------------------------------------------
 * @Name rmVertex3DDiff
 @pstart
 RMenum rmVertex3DDiff (const RMvertex3D *a,
	                const RMvertex3D *b,
		        RMvertex3D *dst)
 @pend

 @astart
 const RMvertex3D *a, *b - handles to RMvertex3D objects (input).

 RMvertex3D *dst - a handle to an RMvertex3D object (result).
 @aend

 @dstart

 Computes the vector difference (a - b), placing the result into
 "dst".  This is computed such that either a or b may be used as the
 result c.  Returns RM_CHILL upon success, or RM_WHACKED upon failure.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmVertex3DDiff (const RMvertex3D *a,
	        const RMvertex3D *b,
	        RMvertex3D *c)
{
    if ((RM_ASSERT(a, "rmVertex3DDiff() error: the input RMvertex3D object A is NULL") == RM_WHACKED) ||
	(RM_ASSERT(b, "rmVertex3DDiff() error: the input RMvertex3D object B is NULL") == RM_WHACKED) ||
	(RM_ASSERT(c, "rmVertex3DDiff() error: the dest RMvertex3D object is NULL") == RM_WHACKED))
	return(RM_WHACKED);

    /* compute c = a - b */
    c->x = a->x - b->x;
    c->y = a->y - b->y;
    c->z = a->z - b->z;

    return(RM_CHILL);
}


/*
 * ----------------------------------------------------
 * @Name rmVertex3DDot
 @pstart
 double rmVertex3DDot (const RMvertex3D *a,
	               const RMvertex3D *b)
 @pend

 @astart
 const RMvertex3D *a, *b - handles to RMvertex3D objects (input).
 @aend

 @dstart

 Computes the vector dot product of the input RMvertex3D objects "a"
 and "b," returning the result. Upon failure, an error message is
 issued through RMwarning() and 0.0 is returned.

 @dend
 * ----------------------------------------------------
 */
double
rmVertex3DDot(const RMvertex3D *a,
	      const RMvertex3D *b)
{
    double t;

    if ((RM_ASSERT(a, "rmVertex3DDot() error: the input RMvertex3D object A is NULL") == RM_WHACKED) ||
	(RM_ASSERT(b, "rmVertex3DDot() error: the input RMvertex3D object B is NULL") == RM_WHACKED))
	return(0.0);

    /* compute t = <a, b> */
    t = (a->x * b->x) + (a->y * b->y) + (a->z * b->z);
    if (fabs(t) < 0.0001)
	t = 0.0;
    return(t);
}


/*
 * ----------------------------------------------------
 * @Name rmVertex3DCross
 @pstart
 RMenum rmVertex3DCross (RMvertex3D *p,
		         RMvertex3D *r,
			 RMvertex3D *result)
 @pend

 @astart
 RMvertex3D *p, *r - handles to RMvertex3D objects (input).

 RMvertex3D *result - a handle to an RMvertex3D object (result).
 @aend

 @dstart

 Computes the vector cross product of the two input vectors "P" and
 "R", placing the resulting vector into "result". Returns RM_CHILL
 upon success or RM_WHACKED upon failure.

 This routine is written such that "result" may be the same object as
 either P or R.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmVertex3DCross (RMvertex3D *p,
		 RMvertex3D *r,
		 RMvertex3D *result)
{
    RMvertex3D t;

    if ((RM_ASSERT(p, "rmVertex3DCross() error: the input P vector is NULL") == RM_WHACKED) ||
	(RM_ASSERT(r, "rmVertex3DCross() error: the input R vector is NULL") == RM_WHACKED) ||
	(RM_ASSERT(result, "rmVertex3DCross() error: the result vector is NULL") == RM_WHACKED))
	return(RM_WHACKED);

    /* compute c = p x r */
    t.x = (p->y * r->z) - (p->z * r->y);
    t.y = (p->z * r->x) - (p->x * r->z);
    t.z = (p->x * r->y) - (p->y * r->x);

    memcpy((void *)result, (void *)&t, sizeof(RMvertex3D));

    return(RM_CHILL);
}


/*
 * ----------------------------------------------------
 * @Name rmVertex3DNormalize
 @pstart
 RMenum rmVertex3DNormalize (RMvertex3D *toNormalize)
 @pend

 @astart
 RMvertex3D *toNormalize - a handle to the RMvertex3D object
    (modified).
 @aend

 @dstart

 Use this routine to normalize the 3-component vector contained in an
 RMvertex3D object. Inside this routine, the vector magnitude of the
 RMvertex3D object is computed. Next, if the magnitude is non-zero,
 each component of the vector is divided by the magnitude, and
 RM_CHILL is returned. If the magnitude is zero, RM_WHACKED is
 returned.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmVertex3DNormalize (RMvertex3D *a)
{
    double mag;
    RMenum rstat = RM_CHILL;

    if (RM_ASSERT(a, "rmVertex3DNormalize() error: the input RMvertex3D object is NULL") == RM_WHACKED)
	return(RM_WHACKED);

    mag = rmVertex3DMag(a);
    if (mag != 0.0)
    {
        mag = 1.0 / mag;
	a->x *= mag;
	a->y *= mag;
	a->z *= mag;
    }
    else
	rstat = RM_WHACKED;

    return(rstat);
}


/*
 * ----------------------------------------------------
 * @Name rmVertex3DMagNormalize
 @pstart
 RMenum rmVertex3DMagNormalize (RMvertex3D *toNormalize,
		                double *magReturn)
 @pend

 @astart
 RMvertex3D *toNormalize - a handle to an RMvertex3D object
    (modified).

 double *magReturn - a pointer to a double (result).
 @aend

 @dstart

 Use this routine when you want to normalize a vector, and get the
 vector magnitude of the unnormalized vector in one call.

 A return value of RM_CHILL means that all's well, while a return
 value of RM_WHACKED means that either one of the input parameters was
 NULL, or that the input RMvertex3D had a zero-valued vector
 magnitude.

 Upon success, the input RMvertex3D is normalized in-situ, and the
 vector magnitude (used to normalize the vector) is returned in the
 caller-supplied memory.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmVertex3DMagNormalize (RMvertex3D *a,
		        double *m)
{
    double mag;
    RMenum rstat = RM_CHILL;

    if ((RM_ASSERT(a, "rmVertex3DMagNormalize() error: the input RMvertex3D object is NULL") == RM_WHACKED) ||
	(RM_ASSERT(m, "rmVertex3DMagNormalize() error: the input magReturn pointer is NULL") == RM_WHACKED))
	return(RM_WHACKED);

    *m = mag = rmVertex3DMag(a);
    if (mag != 0.0)
    {
        mag = 1.0 / mag;
	a->x *= mag;
	a->y *= mag;
	a->z *= mag;
    }
    else
	rstat = RM_WHACKED;

    return(rstat);
}


/*
 * ----------------------------------------------------
 * @Name rmVertex3DMidpoint
 @pstart
 RMenum rmVertex3DMidpoint (const RMvertex3D *a,
		            const RMvertex3D *b,
			    RMvertex3D *dst)
 @pend

 @astart
 const RMvertex3D *a,*b - handles to RMvertex3D objects (input).

 RMvertex3D *dst - a handle to an RMvertex3D object (result).
 @aend

 @dstart

 Computes the midpoint between the two RMvertex3D objects "a" and "b",
 placing the result into "dst". Returns RM_CHILL upon success, or
 RM_WHACKED upon failure.

 @dend
 * ----------------------------------------------------
 */
RMenum
rmVertex3DMidpoint (const RMvertex3D *a,
		    const RMvertex3D *b,
		    RMvertex3D *r)
{
    RMvertex3D t;

    if ((RM_ASSERT(a, "rmVertex3DMidpoint() error: the input A RMvertex3D is NULL") == RM_WHACKED) ||
	(RM_ASSERT(b, "rmVertex3DMidpoint() error: the input B RMvertex3D is NULL") == RM_WHACKED))
	return(RM_WHACKED);

    t.x = a->x + 0.5 * (b->x - a->x);
    t.y = a->y + 0.5 * (b->y - a->y);
    t.z = a->z + 0.5 * (b->z - a->z);

    memcpy((void *)r, (void *)&t, sizeof(RMvertex3D));

    return(RM_CHILL);
}


/*
 * ----------------------------------------------------
 * @Name rmDCFromWC3
 @pstart
 RMenum rmDCFromWC3 (const RMvertex3D *src,
	             RMvertex3D *dst,
		     int nPoints,
		     const RMcamera3D *cam3d,
		     const RMmatrix *model,
		     float viewPort[4],
		     int windowWidth,
		     int windowHeight)
 @pend

 @astart
 const RMvertex3D *src - an array of RMvertex3D objects (input).

 RMvertex3D *dst - an array of RMvertex3D objects (result).

 int nPoints - an integer value indicating the number of points from
    "src" that will be transformed.

 RMcamera3D *cam3d - a handle to an RMcamera3D object that specifies
    view geometry (input).

 RMmatrix *model - an RMmatrix object specifying a model
    transformation in world coordinates (input). This parameter is
    optional, and NULL may be used in place of an RMmatrix object.

 float viewPort[4] - an array of floating point values (length 4)
    specifying the viewport geometry within the display window. Each
    array value should be in the range 0..1, specifying the viewport
    origin and width/height in NDC coordinates. The array is in the
    following order: xmin, ymin, xmax,ymax. In other words, the
    viewport is specified using two coordinates defining the min/max
    corners of the viewport in NDC space.

 int windowWidth, windowHeight - integer values specifying the width
    and height in pixels of the display window (input).
 @aend

 @dstart

 Use this routine to compute the device coordinates from world
 coordinates.  RM_CHILL is returned upon success, or RM_WHACKED upon
 failure.

 Each RMvertex3D from "src", of which there are nPoints, is
 transformed through a derived projection matrix producing window
 (pixel) coordinates that are placed into "dst." Each dst[i]
 represents the projected window coordinates of src[i]. The derived
 projection matrix consists of an optional modeling transformation
 matrix (model), and a viewing matrix computed from the camera model
 defined by the input RMcamera3D object.

 Note that no window is required, the projected coordinates are
 computed using the window and viewport dimensions specified by the
 input, no actual rendering occurs (OpenGL is not involved with this
 routine).

 @dend
 * ----------------------------------------------------
 */
RMenum
rmDCFromWC3 (const RMvertex3D *a,	/* list of world coord points */
	     RMvertex3D *ret,    	/* returned device coords */
	     int n,			/* how many points? */
	     const RMcamera3D *cam3d,	/* a 3d camera  */
	     const RMmatrix *model,	/* world model matrix */
	     float vp[4],		/* viewport  */
	     int width,			/* window width in pixels */
	     int height)		/* windoe height in pixels */
{
    int      i;
    float    vpWidth, vpHeight;
    float    s[4];
    RMmatrix projection, vpm;
    RMmatrix forward;
    RMmatrix myModel, myView;

    /* given a point in 3D world coords, return its pixel coords inwindow space (replicates the OpenGL pipeline) */
    if ((RM_ASSERT(a, "rmDCFromWC3() error: the input list of coordinates is NULL") == RM_WHACKED) ||
	(RM_ASSERT(ret, "rmDCFromWC3() error: the return RMvertex3D handle is NULL") == RM_WHACKED) ||
	(RM_ASSERT(cam3d, "rmDCFromWC3() error: the input RMcamera3D object is NULL") == RM_WHACKED))
	return(RM_WHACKED);

    rmCamera3DComputeViewMatrix(cam3d, &myView, &projection);

    if (model == NULL)
	rmMatrixIdentity(&myModel);
    else
	rmMatrixCopy(&myModel, model);

/*    rmMatrixInverse(&projection, &projection); */
    rmMatrixMultiply(&myModel, &myView, &forward);
    rmMatrixMultiply(&forward, &projection, &forward);

    /* now, composite has the entire wc to pixel mapping */
    vpWidth = vp[2] - vp[0];
    vpHeight = vp[3] - vp[1];

    rmMatrixIdentity(&vpm);
    private_rmComputeViewportMatrix(vp, vpWidth, vpHeight, &vpm);

    VCOPY(a, (RMvertex3D *)s);

    for (i = 0; i < n; i++)
    {
	int   j;
	float w;
	float s[4];

	VCOPY((a + i), s);
	s[3] = 1.0;

	/* run the point through the model and projection matrices */
	rmPoint4MatrixTransform(s, &forward, s);

	/*
	 * the point is now in "clip coords".
	 * next, convert to ndc space. this is done by dividing by
	 * the homogenous-space coordinate.
	 */
	w = s[3];
	w = 1.0 / w;

	for (j = 0; j < 4; j++)
	    s[j] *= w;

	/* now transform through vpm. this gives us the window coords */
	rmPoint4MatrixTransform(s, &vpm, s);

	VCOPY(s, (ret + i));
    }
    return(RM_CHILL);
}


/*
 * ----------------------------------------------------
 * @Name rmDCFromWC2
 @pstart
 RMenum rmDCFromWC2 (const RMvertex2D *src,
	             RMvertex2D *dst,
		     int nPoints,
		     const RMcamera2D *cam3d,
		     const RMmatrix *model,
		     float viewPort[4],
		     int windowWidth,
		     int windowHeight)
 @pend

 @astart
 const RMvertex2D *src - an array of RMvertex2D objects (input).

 RMvertex2D *dst - an array of RMvertex2D objects (result).

 int nPoints - an integer value indicating the number of points from
    "src" that will be transformed.

 RMcamera2D *cam2d - a handle to an RMcamera2D object that specifies
    view geometry (input).

 RMmatrix *model - an RMmatrix object specifying a model
    transformation in world coordinates (input). This parameter is
    optional, and NULL may be used in place of an RMmatrix object.

 float viewPort[4] - an array of floating point values (length 4)
    specifying the viewport geometry within the display window. Each
    array value should be in the range 0..1. The array is in the
    following order: xmin, ymin, xmax,ymax. In other words, the
    viewport is specified using two coordinates defining the min/max
    corners of the viewport in NDC space.

 int windowWidth, windowHeight - integer values specifying the width
    and height in pixels of the display window (input).
 @aend

 @dstart

 Use this routine to compute the device coordinates from world
 coordinates.  RM_CHILL is returned upon success, or RM_WHACKED upon
 failure.

 Each RMvertex2D from "src", of which there are nPoints, is
 transformed through a derived projection matrix producing window
 (pixel) coordinates that are placed into "dst." Each dst[i]
 represents the projected window coordinates of src[i]. The derived
 projection matrix consists of an optional modeling transformation
 matrix (model), and a viewing matrix computed from the camera model
 defined by the input RMcamera2D object.

 Note that no window is required, the projected coordinates are
 computed using the window and viewport dimensions specified by the
 input, no actual rendering occurs (OpenGL is not involved with this
 routine).

 @dend
 * ----------------------------------------------------
 */
RMenum
rmDCFromWC2 (const RMvertex2D *a,	/* list of world coord points */
	     RMvertex2D *ret,      	/* returned device coords */
	     int n,			/* how many points? */
	     const RMcamera2D *cam2d,	/* a 3d camera  */
	     const RMmatrix *model,	/* world model matrix */
	     float vp[4],		/* viewport  */
	     int width,			/* window width in pixels */
	     int height)		/* windoe height in pixels */
{
    int      i;
    float    s[4];
    float    vpWidth, vpHeight;
    double   w;
    RMmatrix projection, vpm;
    RMmatrix forward;
    RMmatrix myModel;

    /* given a point in 2D world coords, return its pixel coords inwindow space (replicates the OpenGL pipeline) */
    if ((RM_ASSERT(a, "rmDCFromWC2() error: the input list of coordinates is NULL") == RM_WHACKED) ||
	(RM_ASSERT(ret, "rmDCFromWC2() error: the return RMvertex2D handle is NULL") == RM_WHACKED) ||
	(RM_ASSERT(cam2d, "rmDCFromWC3() error: the input RMcamera2D object is NULL") == RM_WHACKED))
	return(RM_WHACKED);


    rmCamera2DComputeViewMatrix(cam2d, &projection);

    if (model == NULL)
	rmMatrixIdentity(&myModel);
    else
	rmMatrixCopy(&myModel, model);

    rmMatrixMultiply(&myModel, &projection, &forward);

    /* now, composite has the entire wc to pixel mapping */
    vpWidth = vp[2] - vp[0];
    vpHeight = vp[3] - vp[1];

    rmMatrixIdentity(&vpm);
    private_rmComputeViewportMatrix(vp, vpWidth, vpHeight, &vpm);

    for (i = 0; i < n; i++)
    {
	int j;

	s[0] = a[i].x;
	s[1] = a[i].y;
	s[2] = 0.0;
	s[3] = 1.0;

	/* run the point through the model and projection matrices */
	rmPoint4MatrixTransform(s, &forward, s);

	/*
	 * the point is now in "clip coords".
	 * next, convert to ndc space. this is done by dividing by
	 * the homogenous-space coordinate.
	 */
	w = s[3];
	w = 1.0 / w;

	for (j = 0; j < 4; j++)
	    s[j] *= w;

	/* now transform through vpm. this gives us the window coords */
	rmPoint4MatrixTransform(s, &vpm, s);

	ret[i].x = s[0];
	ret[i].y = s[1];

    }
    return(RM_CHILL);
}


/* PRIVATE
 *
 * this routine is like the public routine, but the transformation
 * information is contained within the RMstate object. this routine is
 * useful from within inside a draw routine to obtain the projected
 * pixel coordinates of a vertex, perhaps to perform some view dependent
 * operations.
 */
void
private_rmDCFromWC3 (RMvertex3D *a,	/* list of world coord points */
		     RMvertex3D *ret,	/* returned device coords */
		     int n,		/* how many points? */
		     RMstate *state)
{
    int      i;
    float    s[4];
    RMmatrix vpm;
    RMmatrix forward;

    /* given a point in 3D world coords, return its pixel coords inwindow space (replicates the OpenGL pipeline) */
    /*    rmMatrixMultiply(&(state->model), &(state->view), &forward); */
    rmMatrixMultiply(&(state->modelView), &(state->projection), &forward);

    /* now, composite has the entire wc to pixel mapping */
    rmMatrixIdentity(&vpm);

    private_rmComputeViewportMatrix(state->vp, (float)state->w,
				    (float)state->h, &vpm);

    VCOPY(a,(RMvertex3D *)s);

    for (i = 0; i < n; i++)
    {
	int   j;
	float w;

	VCOPY((a + i), s);
	s[3] = 1.0;

	/* run the point through the model and projection matrices */
	rmPoint4MatrixTransform(s, &forward, s);

	/*
	 * the point is now in "clip coords".
	 * next, convert to ndc space. this is done by dividing by
	 * the homogenous-space coordinate.
	 */
	w = s[3];
	w = 1.0 / w;

	for (j = 0; j < 4; j++)
	    s[j] *= w;

	/* now transform through vpm. this gives us the window coords */
	rmPoint4MatrixTransform(s, &vpm, s);
	VCOPY(s, (ret + i));
    }
}


/* PRIVATE
 *
 * this routine is like the public routine, but the transformation
 * information is contained within the RMstate object. this routine is
 * useful from within inside a draw routine to obtain the projected
 * pixel coordinates of a vertex, perhaps to perform some view dependent
 * operations.
 */
void
private_rmDCFromWC2 (RMvertex3D *a,	/* list of world coord points */
		     RMvertex3D *ret,	/* returned device coords */
		     int n,		/* how many points? */
		     RMstate *state)
{
    RMmatrix vpm;
    RMmatrix forward;
    int      i;
    float    s[4];

    /* given a point in 2D world coords, return its pixel coords inwindow space (replicates the OpenGL pipeline) */
#if 0
    rmMatrixMultiply(&(state->model), &(state->view), &forward);
    rmMatrixMultiply(&forward, &(state->projection), &forward);
#endif

    /* 7/28/01 - use the modelView matrix, rather than the separate
     model & view matrices. The separate model & view matrices are no
     longer guaranteed to be valid in the multistage rendering, but
     modelView will be valid. There are no explicit commands to update
     the model & view matrices separately; these are subsumed by a
     single modelView matrix in order to more closely model the OpenGL
     matrix stack. */

    rmMatrixMultiply(&(state->modelView), &(state->projection), &forward);

    /* now, composite has the entire wc to pixel mapping */
    rmMatrixIdentity(&vpm);

    private_rmComputeViewportMatrix(state->vp, (float)state->w,
				    (float)state->h, &vpm);

#if 0
    vpm.m[0][0] = state->vp[2] * 0.5; /*  * state->w; */
    vpm.m[1][1] = state->vp[3] * 0.5; /*  * state->h; */
    vpm.m[3][0] = vpm.m[0][0] + ((state->w  * state->vp[0]) / state->w);
    vpm.m[3][1] = vpm.m[1][1] + ((state->h  * state->vp[1]) / state->h);
#endif

    for (i = 0; i < n; i++)
    {
	int   j;
	float w;

	s[0] = a[i].x;
	s[1] = a[i].y;
	s[2] = 0.0;
	s[3] = 1.0;

	/* run the point through the model and projection matrices */
	rmPoint4MatrixTransform(s,&forward,s);

	/*
	 * the point is now in "clip coords".
	 * next, convert to ndc space. this is done by dividing by
	 * the homogenous-space coordinate.
	 */
	w = s[3];
	w = 1.0 / w;

	for (j = 0; j < 4; j++)
	    s[j] *= w;

	/* now transform through vpm. this gives us the window coords */
	rmPoint4MatrixTransform(s, &vpm, s);
	ret[i].x = s[0];
	ret[i].y = s[1];
	ret[i].z = s[2];
    }
}


/*
 * the following routines do conversions to/from floating point
 * and fixed point.
 *
 * these conversion routines are designed to represent numbers
 * in the range -2048..2048. we use 11 bits for representing the
 * integer part, and 8 bits for representing the fractional part..
 * that leads to a fractional resolution (potential error) of 0.00390625.
 */

/* PRIVATE */
int
private_rmFixedFromFloat (float a)
{
    /* we only do positive numbers */
    float        whole, frac;
    unsigned int iwhole, ifrac, isign;

    if (a < 0.0)
    {
	isign = (1 << 19);
	a = -1.0 * a;
    }
    else
	isign = 0;

    whole = (float)((int)a);
    frac = a - whole;

    frac = frac * 256.0;
    ifrac = (int)frac;
    ifrac &= 0x0ff;		/* keep 8 bits of fraction */

    iwhole = (int)whole;
    iwhole = iwhole & 0x07FF;	/* keep 11 bits */
    iwhole = iwhole << 8;	/* make room for the fraction */

    return(iwhole | ifrac | isign);
}


/* PRIVATE */
float
private_rmFloatFromFixed (int a)
{
    float frac, whole, result;
    int   ifrac, iwhole, isign;

    isign = a & (1 << 19);
    ifrac = a & 0x0ff;		/* grab the fraction */
    iwhole = (a >> 8) & 0x07ff;	/* grab the integer part */

    frac = (float)ifrac * 0.003921568; /* 1/255 */
    whole = (float)iwhole;

    result = whole + frac;

    if (isign)
	result *= -1.0;

    return(result);
}

/* PRIVATE */
void
private_rmNewWCfromDCOffset(RMvertex3D *source,	/* input 3d coordinate */
			    float xmove, 	/* pixel offset */
			    float ymove, 	/* pixel offset */
			    RMvertex3D *dest, 	/* result */
			    RMmatrix *vpm, 	/* viewport matrix */
			    RMmatrix *vpm_inv, 	/* inverse viewport matrix */
			    RMmatrix *forward, 	/* model+view+projection */
			    RMmatrix *inv) 	/* model/view/projection inverse */
{
    /*
     * given a world coordinate "source", what is the new
     * world coordinate that is the sum of "source" and
     * (xmove, ymove) where xmove, ymove are pixel offsets from
     * the projected point "source"?
     *
     * in some sense, this problem is underspecified because we have
     * no Z information in the delta. we just assume the resulting
     * coordinate will have the same z-value as the input coordinate.
     */
    float      f[4], f2[4];
    float      save_z, w;
    RMvertex3D t2;

    VCOPY(source, (RMvertex3D *)f);
    f[3] = 1.0;

    rmPoint4MatrixTransform(f, forward, f);

    save_z = f[2];
    w = f[3];

    VCOPY(f,&t2);

    t2.x = t2.x / w;
    t2.y = t2.y / w;
    t2.z = t2.z / w;

    rmPointMatrixTransform(&t2, vpm, &t2);

    t2.x += xmove;
    t2.y += ymove;

    rmPointMatrixTransform(&t2, vpm_inv, (RMvertex3D *)f);
    f[0] *= w;
    f[1] *= w;
    f[2] *= w;
    f[3] = w;
    rmPoint4MatrixTransform(f, inv, f2);
    VCOPY(f2, dest);
}

/* PRIVATE */
void
private_rmNewWCfromDC(RMvertex3D *source,	/* input 3d coordinate */
		      RMvertex3D *dest, 	/* result */
		      RMmatrix *vpm, 	/* viewport matrix */
		      RMmatrix *vpm_inv, 	/* inverse viewport matrix */
		      RMmatrix *forward, 	/* model+view+projection */
		      RMmatrix *inv)  /* model/view/projection inverse */
{
    /*
     * given an input device coordinate in "source", what is the
     * corresponding world coordinate?
     */
    float      f[4], f2[4];
    float      save_z, w;
    RMvertex3D t2;

    VCOPY(source, (RMvertex3D *)f);
    f[3] = 1.0;

    rmPoint4MatrixTransform(f, forward, f);

    save_z = f[2];
    w = f[3];

    VCOPY(f,&t2);

    t2.x = t2.x / w;
    t2.y = t2.y / w;
    t2.z = t2.z / w;

    rmPointMatrixTransform(&t2, vpm, &t2);

    rmPointMatrixTransform(&t2, vpm_inv, (RMvertex3D *)f);
    f[0] *= w;
    f[1] *= w;
    f[2] = save_z;
    f[3] = w;
    rmPoint4MatrixTransform(f, inv, f2);
    VCOPY(f2, dest);
}

/* PRIVATE */
void
private_rmComputeViewportMatrix(float vp[4],
				float fw,
				float fh,
				RMmatrix *returnMatrix)
{
    RMmatrix *m = returnMatrix;

    /*
     * the following formulation should work in both cases, where:
     * (a) the vp AND fw/fh are NDC coordinates/values, and
     * (b) the vp AND fw/fh are pixel values.
     *
     * vp should be xmin, ymin, xwidth, yheight
     */
    m->m[0][0] = vp[2] * 0.5F;
    m->m[1][1] = vp[3] * 0.5F;
    m->m[3][0] = m->m[0][0] + fw * vp[0];
    m->m[3][1] = m->m[1][1] + fh * vp[1];
}
/* EOF */