veclib4d.c

[Game.Programming].Academic - Graphics Gems (6 books source code)

#include <stdio.h>
#include <stdlib.h>	/* for malloc() definition */
#include <math.h>
#include "GGems.h"
#include "GGems4d.h"

/* 4D VECTOR LIBRARY.
 *
 * Steve Hill, Computing Laboratory, University of Kent, UK
 *
 * Email: [e-mail]
 *
 * Included with Graphics Gems V, Academic Press ed. Alan Paeth.
 *
 * Extends the standard Graphics Gems library to four dimensions
 * and provides functions to create and manipulate homogeneous
 * transformations for 4-vectors (including projections).
 */

/* V4SquaredLength()
 *
 * Compute the square of the magnitude of a 4-vector
 */

double
V4SquaredLength(v)
Vector4 *v;
{
	return v->x * v->x +
	       v->y * v->y +
	       v->z * v->z +
	       v->w * v->w;
}

/* V4Length()
 *
 * Compute the magnitude of a 4-vector
 * If square of magnitude is required use V4SquaredLength().
 */

double
V4Length(v)
Vector4	*v;
{
	return sqrt(V4SquaredLength(v));
}

/* V4Negate()
 *
 * Negate vector elementwise.
 */

Vector4 *
V4Negate(v)
Vector4	*v;
{
	v->x = -v->x;
	v->y = -v->y;
	v->z = -v->z;
	v->w = -v->w;

	return v;
}

/* V4Normalize()
 *
 * Convert vector to unit magnitude.
 * A zero-length vector is unchanged.
 */

Vector4 *
V4Normalize(v)
Vector4	*v;
{
	double	len = V4Length(v);

	if (len != 0.0)
	{
		v->x /= len;
		v->y /= len;
		v->z /= len;
		v->w /= len;
	}

	return v;
}

/* V4Scale()
 *
 * Scale vector to requested magnitude.
 * Zero length vectors unchanged.
 */

Vector4 *
V4Scale(v, newlen)
Vector4	*v;
double	newlen;
{
	double	len = V4Length(v);

	if (len != 0.0)
	{
		v->x *= newlen/len;
		v->y *= newlen/len;
		v->z *= newlen/len;
		v->w *= newlen/len;
	}

	return v;
}

/* V4Add()
 *
 * Add (elementwise) two 4-vectors
 */

Vector4 *
V4Add(a, b, result)
Vector4	*a, *b, *result;
{
	result->x = a->x + b-> x;
	result->y = a->y + b-> y;
	result->z = a->z + b-> z;
	result->w = a->w + b-> w;

	return result;
}

/* V4Sub()
 *
 * Subtract (elementwise) vector b from vector a.
 */

Vector4 *
V4Sub(a, b, result)
Vector4	*a, *b, *result;
{
	result->x = a->x - b-> x;
	result->y = a->y - b-> y;
	result->z = a->z - b-> z;
	result->w = a->w - b-> w;

	return result;
}

/* V4Dot()
 *
 * Returns the 4-space dot-product of two vectors.
 */

double
V4Dot(a, b)
Vector4	*a, *b;
{
	return	a->x * b->x +
		a->y * b->y +
		a->z * b->z +
		a->w * b->w;
}

/* V4Lerp()
 *
 * Calculates vector that is a linear interpolation between
 * lo and hi, accoreding to alpha.
 * See also GraphicsGems.h (GGems.h)
 */

Vector4 *
V4Lerp(lo, hi, alpha, result)
Vector4	*lo, *hi, *result;
double	alpha;
{
	result->x = LERP(alpha, lo->x, hi->x);
	result->y = LERP(alpha, lo->y, hi->y);
	result->z = LERP(alpha, lo->z, hi->z);
	result->w = LERP(alpha, lo->w, hi->w);

	return result;
}

/* V4Combine()
 *
 * Calculate a linear combination of two vectors.
 */

Vector4 *
V4Combine(a, b, result, ascl, bscl)
Vector4	*a, *b, *result;
double	ascl, bscl;
{
	result->x = ascl * a->x + bscl * b->x;
	result->y = ascl * a->y + bscl * b->y;
	result->z = ascl * a->z + bscl * b->z;
	result->w = ascl * a->w + bscl * b->w;

	return result;
}

/* V4Mul()
 *
 * Multiply (elementwise) two vectors.
 */

Vector4 *
V4Mul(a, b, result)
Vector4	*a, *b, *result;
{
	result->x = a->x * b->x;
	result->y = a->y * b->y;
	result->z = a->z * b->z;
	result->w = a->w * b->w;

	return result;
}

/* V4DistanceBetween2Points()
 *
 * Calculates the distance between two 4-space points.
 */

double
V4DistanceBetween2Points(a, b)
Point4	*a, *b;
{
	double	dx = a->x - b->x,
		dy = a->y - b->y,
		dz = a->z - b->z,
		dw = a->w - b->w;

	return sqrt(dx*dx + dy*dy + dz*dz + dw*dw);
}

/* V4Cross()
 *
 * Calculates the cross product of three 4-space vectors.
 */

Vector4 *
V4Cross(a, b, c, result)
Vector4	*a, *b, *c, *result;
{
	double	d1, d2, d3, d4, d5, d6;

	d1 = (b->z * c->w) - (b->w * c->z);
	d2 = (b->y * c->w) - (b->w * c->y);
	d3 = (b->y * c->z) - (b->z * c->y);
	d4 = (b->x * c->w) - (b->w * c->x);
	d5 = (b->x * c->z) - (b->z * c->x);
	d6 = (b->x * c->y) - (b->y * c->x);

	result->x = - a->y * d1 + a->z * d2 - a->w * d3;
	result->y =   a->x * d1 - a->z * d4 + a->w * d5;
	result->z = - a->x * d2 + a->y * d4 - a->w * d6;
	result->w =   a->x * d3 - a->y * d5 + a->z * d6;

	return result;
}

/* V4New()
 *
 * Allocate and initialise a new 4-vector
 */

Vector4 *
V4New(x, y, z, w)
double	x, y, z, w;
{
	Vector4	*v = NEWTYPE(Vector4);

	v->x = x;
	v->y = y;
	v->z = z;
	v->w = w;

	return v;
}

/* V4Duplicate()
 *
 * Create a copy of a 4-vector.
 */

Vector4 *
V4Duplicate(a)
Vector4	*a;
{
	Vector4	*v = NEWTYPE(Vector4);

	v->x = a->x;
	v->y = a->y;
	v->z = a->z;
	v->w = a->w;

	return v;
}

/* V4MatPrint()
 *
 * Diagnostic function to print out a 5x5 matrix
 *
 */

void
V4MatPrint(file, mat)
FILE	*file;
Matrix5	*mat;
{
	int	i, j;

	for (i = 0; i < 5; i += 1)
	{
		for (j = 0; j < 5; j += 1)
			fprintf(file, "%lf ", mat->element[i][j]);
		putc('\n', file);
	}
}

/* V4MatCopy()
 *
 * Copy a 5x5 matrix.
 */

Matrix5 *
V4MatCopy(from, to)
Matrix5	*from, *to;
{
	int	i, j;

	for (i = 0; i < 5; i += 1)
	for (j = 0; j < 5; j += 1)
		to->element[i][j] = from->element[i][j];

	return to;
}

/* V4MulPointByMatrix()
 *
 * Apply 4x4 transformation matrix to 4-space point.
 *
 * In common with the standard Graphics Gems library, points/vectors
 * are considered to be row vectors, hence multiplication is
 * post-multiplication.  The transformation T o S o R represents
 * a translation followed by a scale followed by a rotation.
 * Where o is matrix multiplication (see V4MatMul() and V4MatMul2() below).
 */

Point4 *
V4MulPointByMatrix(pin, m, pout)
Point4	*pin, *pout;
Matrix4	*m;
{
	pout->x = pin->x * m->element[0][0] +
		  pin->y * m->element[1][0] +
		  pin->z * m->element[2][0] +
		  pin->w * m->element[3][0];
	pout->y = pin->x * m->element[0][1] +
		  pin->y * m->element[1][1] +
		  pin->z * m->element[2][1] +
		  pin->w * m->element[3][1];
	pout->z = pin->x * m->element[0][2] +
		  pin->y * m->element[1][2] +
		  pin->z * m->element[2][2] +
		  pin->w * m->element[3][2];
	pout->w = pin->x * m->element[0][3] +
		  pin->y * m->element[1][3] +
		  pin->z * m->element[2][3] +
		  pin->w * m->element[3][3];

	return pout;
}

/* V4MulPointByProjMatrix()
 *
 * Apply 5x5 (projection) matrix to 4-space point
 *
 * For example, the projection matrix:
 *
 *   [  1  0  0  0    0    ]
 *   [  0  1  0  0    0    ]
 *   [  0  0  0  0   1/dz  ]
 *   [  0  0  0  0   1/dw  ]
 *   [  0  0  0  0    1    ]
 *
 * represents a combined perspective/parallel projection along
 * the Z and W axes.  As dz or dw tend to infinity the projections
 * along those axes becomes parallel.
 *
 * 5x5 matrices may also be used to perform translations.
 *
 *   [  1  0  0  0  0 ]
 *   [  0  1  0  0  0 ]
 *   [  0  0  1  0  0 ]
 *   [  0  0  0  1  0 ]
 *   [  x  y  z  w  1 ]
 */

Point4 *
V4MulPointByProjMatrix(pin, m, pout)
Point4	*pin, *pout;
Matrix5	*m;
{
	double	v;

	pout->x = pin->x * m->element[0][0] +
		  pin->y * m->element[1][0] +
		  pin->z * m->element[2][0] +
		  pin->w * m->element[3][0] +
		  m->element[4][0];
	pout->y = pin->x * m->element[0][1] +
		  pin->y * m->element[1][1] +
		  pin->z * m->element[2][1] +
		  pin->w * m->element[3][1] +
		  m->element[4][1];
	pout->z = pin->x * m->element[0][2] +
		  pin->y * m->element[1][2] +
		  pin->z * m->element[2][2] +
		  pin->w * m->element[3][2] +
		  m->element[4][2];
	pout->w = pin->x * m->element[0][3] +
		  pin->y * m->element[1][3] +
		  pin->z * m->element[2][3] +
		  pin->w * m->element[3][3] +
		  m->element[4][3];

	v = pin->x * m->element[0][4] +
	    pin->y * m->element[1][4] +
	    pin->z * m->element[2][4] +
	    pin->w * m->element[3][4] +
	    m->element[4][4];

	if (v != 0.0)
	{
		pout->x /= v;
		pout->y /= v;
		pout->z /= v;
		pout->w /= v;
	}

	return pout;
}

/* V4MatMul()
 *
 * Multiply two 5x5 matrices.  Result matrix must be distinct
 * from argument matrices.
 */

Matrix5 *
V4MatMul(a, b, result)
Matrix5	*a, *b, *result;
{
	int	i, j, k;

	for (i = 0; i < 5; i++)
	for (j = 0; j < 5; j++)
	{
		double t = 0.0;

		for (k = 0; k < 5; k++)
			t += a->element[i][k] * b->element[k][j];

		result->element[i][j] = t;
	}

	return result;
}

/* V4MatMul2
 *
 * Variation of V4MatMul() where result matrix may be the same
 * as matrices a or b.
 */

Matrix5 *
V4MatMul2(a, b, result)
Matrix5 *a, *b, *result;
{
	Matrix5	tmp;

	V4MatMul(a, b, &tmp);
	V4MatCopy(&tmp, result);
	return result;
}

/* 4D TRANSFORMATIONS
 *
 * The following functions provide various methods
 * to apply and manipulate transformations.
 */

/* V4MatZero()
 *
 * Clear a 5x5 matrix.
 */

Matrix5 *
V4MatZero(mat)
Matrix5	*mat;
{
	int	i, j;

	for (i = 0; i < 5; i += 1)
	for (j = 0; j < 5; j += 1)
		mat->element[i][j] = 0.0;
}

/* STANDARD TRANSFORMATIONS
 *
 * The following functions initialise the standard affine transformations