tenebrae/tenebrae2
0
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity
tenebrae2/mathlib.c
cholleme 92f88376bc Added some usefull math routines
2003-05-04 21:43:46 +00:00

735 lines
15 KiB
C
Raw Blame History

/*
Copyright (C) 1996-1997 Id Software, Inc.
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
// mathlib.c -- math primitives
#include <math.h>
#include "quakedef.h"
void Sys_Error (char *error, ...);
vec3_t vec3_origin = {0,0,0};
int nanmask = 255<<23;
/*-----------------------------------------------------------------*/
void ProjectPointOnPlane( vec3_t dst, const vec3_t p, const vec3_t normal )
{
float d;
vec3_t n;
float inv_denom;
inv_denom = 1.0F / DotProduct( normal, normal );
d = DotProduct( normal, p ) * inv_denom;
n[0] = normal[0] * inv_denom;
n[1] = normal[1] * inv_denom;
n[2] = normal[2] * inv_denom;
dst[0] = p[0] - d * n[0];
dst[1] = p[1] - d * n[1];
dst[2] = p[2] - d * n[2];
}
/*
** assumes "src" is normalized
*/
/*void PerpendicularVector( vec3_t dst, const vec3_t src )
{
int pos;
int i;
float minelem = 1.0F;
vec3_t tempvec;
for ( pos = 0, i = 0; i < 3; i++ )
{
if ( fabs( src[i] ) < minelem )
{
pos = i;
minelem = fabs( src[i] );
}
}
tempvec[0] = tempvec[1] = tempvec[2] = 0.0F;
tempvec[pos] = 1.0F;
ProjectPointOnPlane( dst, tempvec, src );
VectorNormalize( dst );
}*/
void PerpendicularVector( vec3_t dst, const vec3_t src ) //Optimized a bit :) - Eradicator
{
int pos;
float minelem;
if (src[0])
{
dst[0] = 0;
if (src[1])
{
dst[1] = 0;
if (src[2])
{
dst[2] = 0;
pos = 0;
minelem = fabs(src[0]);
if (fabs(src[1]) < minelem)
{
pos = 1;
minelem = fabs(src[1]);
}
if (fabs(src[2]) < minelem)
pos = 2;
dst[pos] = 1;
dst[0] -= src[pos] * src[0];
dst[1] -= src[pos] * src[1];
dst[2] -= src[pos] * src[2];
VectorNormalize(dst);
}
else
dst[2] = 1;
}
else
{
dst[1] = 1;
dst[2] = 0;
}
}
else
{
dst[0] = 1;
dst[1] = 0;
dst[2] = 0;
}
}
#ifdef _WIN32
#pragma optimize( "", off )
#endif
void RotatePointAroundVector( vec3_t dst, const vec3_t dir, const vec3_t point, float degrees )
{
float m[3][3];
float im[3][3];
float zrot[3][3];
float tmpmat[3][3];
float rot[3][3];
int i;
vec3_t vr, vup, vf;
vf[0] = dir[0];
vf[1] = dir[1];
vf[2] = dir[2];
PerpendicularVector( vr, dir );
CrossProduct( vr, vf, vup );
m[0][0] = vr[0];
m[1][0] = vr[1];
m[2][0] = vr[2];
m[0][1] = vup[0];
m[1][1] = vup[1];
m[2][1] = vup[2];
m[0][2] = vf[0];
m[1][2] = vf[1];
m[2][2] = vf[2];
memcpy( im, m, sizeof( im ) );
im[0][1] = m[1][0];
im[0][2] = m[2][0];
im[1][0] = m[0][1];
im[1][2] = m[2][1];
im[2][0] = m[0][2];
im[2][1] = m[1][2];
memset( zrot, 0, sizeof( zrot ) );
zrot[0][0] = zrot[1][1] = zrot[2][2] = 1.0F;
zrot[0][0] = cos( DEG2RAD( degrees ) );
zrot[0][1] = sin( DEG2RAD( degrees ) );
zrot[1][0] = -sin( DEG2RAD( degrees ) );
zrot[1][1] = cos( DEG2RAD( degrees ) );
R_ConcatRotations( m, zrot, tmpmat );
R_ConcatRotations( tmpmat, im, rot );
for ( i = 0; i < 3; i++ )
{
dst[i] = rot[i][0] * point[0] + rot[i][1] * point[1] + rot[i][2] * point[2];
}
}
#ifdef _WIN32
#pragma optimize( "", on )
#endif
/*-----------------------------------------------------------------*/
float anglemod(float a)
{
#if 0
if (a >= 0)
a -= 360*(int)(a/360);
else
a += 360*( 1 + (int)(-a/360) );
#endif
a = (360.0/65536) * ((int)(a*(65536/360.0)) & 65535);
return a;
}
/*
==================
BOPS_Error
Split out like this for ASM to call.
==================
*/
void BOPS_Error (void)
{
Sys_Error ("BoxOnPlaneSide: Bad signbits");
}
#if !id386
/*
==================
BoxOnPlaneSide
Returns 1, 2, or 1 + 2
==================
*/
int BoxOnPlaneSide (vec3_t emins, vec3_t emaxs, mplane_t *p)
{
float dist1, dist2;
int sides;
#if 0 // this is done by the BOX_ON_PLANE_SIDE macro before calling this
// function
// fast axial cases
if (p->type < 3)
{
if (p->dist <= emins[p->type])
return 1;
if (p->dist >= emaxs[p->type])
return 2;
return 3;
}
#endif
// general case
switch (p->signbits)
{
case 0:
dist1 = p->normal[0]*emaxs[0] + p->normal[1]*emaxs[1] + p->normal[2]*emaxs[2];
dist2 = p->normal[0]*emins[0] + p->normal[1]*emins[1] + p->normal[2]*emins[2];
break;
case 1:
dist1 = p->normal[0]*emins[0] + p->normal[1]*emaxs[1] + p->normal[2]*emaxs[2];
dist2 = p->normal[0]*emaxs[0] + p->normal[1]*emins[1] + p->normal[2]*emins[2];
break;
case 2:
dist1 = p->normal[0]*emaxs[0] + p->normal[1]*emins[1] + p->normal[2]*emaxs[2];
dist2 = p->normal[0]*emins[0] + p->normal[1]*emaxs[1] + p->normal[2]*emins[2];
break;
case 3:
dist1 = p->normal[0]*emins[0] + p->normal[1]*emins[1] + p->normal[2]*emaxs[2];
dist2 = p->normal[0]*emaxs[0] + p->normal[1]*emaxs[1] + p->normal[2]*emins[2];
break;
case 4:
dist1 = p->normal[0]*emaxs[0] + p->normal[1]*emaxs[1] + p->normal[2]*emins[2];
dist2 = p->normal[0]*emins[0] + p->normal[1]*emins[1] + p->normal[2]*emaxs[2];
break;
case 5:
dist1 = p->normal[0]*emins[0] + p->normal[1]*emaxs[1] + p->normal[2]*emins[2];
dist2 = p->normal[0]*emaxs[0] + p->normal[1]*emins[1] + p->normal[2]*emaxs[2];
break;
case 6:
dist1 = p->normal[0]*emaxs[0] + p->normal[1]*emins[1] + p->normal[2]*emins[2];
dist2 = p->normal[0]*emins[0] + p->normal[1]*emaxs[1] + p->normal[2]*emaxs[2];
break;
case 7:
dist1 = p->normal[0]*emins[0] + p->normal[1]*emins[1] + p->normal[2]*emins[2];
dist2 = p->normal[0]*emaxs[0] + p->normal[1]*emaxs[1] + p->normal[2]*emaxs[2];
break;
default:
dist1 = dist2 = 0; // shut up compiler
BOPS_Error ();
break;
}
#if 0
int i;
vec3_t corners[2];
for (i=0 ; i<3 ; i++)
{
if (plane->normal[i] < 0)
{
corners[0][i] = emins[i];
corners[1][i] = emaxs[i];
}
else
{
corners[1][i] = emins[i];
corners[0][i] = emaxs[i];
}
}
dist = DotProduct (plane->normal, corners[0]) - plane->dist;
dist2 = DotProduct (plane->normal, corners[1]) - plane->dist;
sides = 0;
if (dist1 >= 0)
sides = 1;
if (dist2 < 0)
sides |= 2;
#endif
sides = 0;
if (dist1 >= p->dist)
sides = 1;
if (dist2 < p->dist)
sides |= 2;
#ifdef PARANOID
if (sides == 0)
Sys_Error ("BoxOnPlaneSide: sides==0");
#endif
return sides;
}
#endif
void AngleVectors (vec3_t angles, vec3_t forward, vec3_t right, vec3_t up)
{
float angle;
float sr, sp, sy, cr, cp, cy;
angle = angles[YAW] * (M_PI*2 / 360);
sy = sin(angle);
cy = cos(angle);
angle = angles[PITCH] * (M_PI*2 / 360);
sp = sin(angle);
cp = cos(angle);
angle = angles[ROLL] * (M_PI*2 / 360);
sr = sin(angle);
cr = cos(angle);
forward[0] = cp*cy;
forward[1] = cp*sy;
forward[2] = -sp;
right[0] = (-1*sr*sp*cy+-1*cr*-sy);
right[1] = (-1*sr*sp*sy+-1*cr*cy);
right[2] = -1*sr*cp;
up[0] = (cr*sp*cy+-sr*-sy);
up[1] = (cr*sp*sy+-sr*cy);
up[2] = cr*cp;
}
int VectorCompare (vec3_t v1, vec3_t v2)
{
int i;
for (i=0 ; i<3 ; i++)
if (v1[i] != v2[i])
return 0;
return 1;
}
void VectorMA (vec3_t veca, float scale, vec3_t vecb, vec3_t vecc)
{
vecc[0] = veca[0] + scale*vecb[0];
vecc[1] = veca[1] + scale*vecb[1];
vecc[2] = veca[2] + scale*vecb[2];
}
vec_t _DotProduct (vec3_t v1, vec3_t v2)
{
return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
}
void _VectorSubtract (vec3_t veca, vec3_t vecb, vec3_t out)
{
out[0] = veca[0]-vecb[0];
out[1] = veca[1]-vecb[1];
out[2] = veca[2]-vecb[2];
}
void _VectorAdd (vec3_t veca, vec3_t vecb, vec3_t out)
{
out[0] = veca[0]+vecb[0];
out[1] = veca[1]+vecb[1];
out[2] = veca[2]+vecb[2];
}
void _VectorCopy (vec3_t in, vec3_t out)
{
out[0] = in[0];
out[1] = in[1];
out[2] = in[2];
}
void VectorConstruct(float v1, float v2, float v3, vec3_t v) {
v[0] = v1;
v[1] = v2;
v[2] = v3;
}
void CrossProduct (vec3_t v1, vec3_t v2, vec3_t cross)
{
cross[0] = v1[1]*v2[2] - v1[2]*v2[1];
cross[1] = v1[2]*v2[0] - v1[0]*v2[2];
cross[2] = v1[0]*v2[1] - v1[1]*v2[0];
}
double sqrt(double x);
vec_t Length(vec3_t v)
{
int i;
float length;
length = 0;
for (i=0 ; i< 3 ; i++)
length += v[i]*v[i];
length = sqrt (length); // FIXME
return length;
}
vec_t Distance(vec3_t v, vec3_t v2)
{
int i;
float length, d;
length = 0;
for (i=0 ; i< 3 ; i++)
d = v[i] - v2[i];
length += d*d;
length = sqrt (length); // FIXME
return length;
}
float VectorNormalize (vec3_t v)
{
float length, ilength;
length = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
length = sqrt (length); // FIXME
if (length)
{
ilength = 1/length;
v[0] *= ilength;
v[1] *= ilength;
v[2] *= ilength;
}
return length;
}
void VectorInverse (vec3_t v)
{
v[0] = -v[0];
v[1] = -v[1];
v[2] = -v[2];
}
void VectorScale (vec3_t in, vec_t scale, vec3_t out)
{
out[0] = in[0]*scale;
out[1] = in[1]*scale;
out[2] = in[2]*scale;
}
int Q_log2(int val)
{
int answer=0;
while (val>>=1)
answer++;
return answer;
}
/*
================
R_ConcatRotations
================
*/
void R_ConcatRotations (float in1[3][3], float in2[3][3], float out[3][3])
{
out[0][0] = in1[0][0] * in2[0][0] + in1[0][1] * in2[1][0] +
in1[0][2] * in2[2][0];
out[0][1] = in1[0][0] * in2[0][1] + in1[0][1] * in2[1][1] +
in1[0][2] * in2[2][1];
out[0][2] = in1[0][0] * in2[0][2] + in1[0][1] * in2[1][2] +
in1[0][2] * in2[2][2];
out[1][0] = in1[1][0] * in2[0][0] + in1[1][1] * in2[1][0] +
in1[1][2] * in2[2][0];
out[1][1] = in1[1][0] * in2[0][1] + in1[1][1] * in2[1][1] +
in1[1][2] * in2[2][1];
out[1][2] = in1[1][0] * in2[0][2] + in1[1][1] * in2[1][2] +
in1[1][2] * in2[2][2];
out[2][0] = in1[2][0] * in2[0][0] + in1[2][1] * in2[1][0] +
in1[2][2] * in2[2][0];
out[2][1] = in1[2][0] * in2[0][1] + in1[2][1] * in2[1][1] +
in1[2][2] * in2[2][1];
out[2][2] = in1[2][0] * in2[0][2] + in1[2][1] * in2[1][2] +
in1[2][2] * in2[2][2];
}
/*
================
R_ConcatTransforms
================
*/
void R_ConcatTransforms (float in1[3][4], float in2[3][4], float out[3][4])
{
out[0][0] = in1[0][0] * in2[0][0] + in1[0][1] * in2[1][0] +
in1[0][2] * in2[2][0];
out[0][1] = in1[0][0] * in2[0][1] + in1[0][1] * in2[1][1] +
in1[0][2] * in2[2][1];
out[0][2] = in1[0][0] * in2[0][2] + in1[0][1] * in2[1][2] +
in1[0][2] * in2[2][2];
out[0][3] = in1[0][0] * in2[0][3] + in1[0][1] * in2[1][3] +
in1[0][2] * in2[2][3] + in1[0][3];
out[1][0] = in1[1][0] * in2[0][0] + in1[1][1] * in2[1][0] +
in1[1][2] * in2[2][0];
out[1][1] = in1[1][0] * in2[0][1] + in1[1][1] * in2[1][1] +
in1[1][2] * in2[2][1];
out[1][2] = in1[1][0] * in2[0][2] + in1[1][1] * in2[1][2] +
in1[1][2] * in2[2][2];
out[1][3] = in1[1][0] * in2[0][3] + in1[1][1] * in2[1][3] +
in1[1][2] * in2[2][3] + in1[1][3];
out[2][0] = in1[2][0] * in2[0][0] + in1[2][1] * in2[1][0] +
in1[2][2] * in2[2][0];
out[2][1] = in1[2][0] * in2[0][1] + in1[2][1] * in2[1][1] +
in1[2][2] * in2[2][1];
out[2][2] = in1[2][0] * in2[0][2] + in1[2][1] * in2[1][2] +
in1[2][2] * in2[2][2];
out[2][3] = in1[2][0] * in2[0][3] + in1[2][1] * in2[1][3] +
in1[2][2] * in2[2][3] + in1[2][3];
}
/*
===================
FloorDivMod
Returns mathematically correct (floor-based) quotient and remainder for
numer and denom, both of which should contain no fractional part. The
quotient must fit in 32 bits.
====================
*/
void FloorDivMod (double numer, double denom, int *quotient,
int *rem)
{
int q, r;
double x;
#ifndef PARANOID
if (denom <= 0.0)
Sys_Error ("FloorDivMod: bad denominator %d\n", denom);
// if ((floor(numer) != numer) || (floor(denom) != denom))
// Sys_Error ("FloorDivMod: non-integer numer or denom %f %f\n",
// numer, denom);
#endif
if (numer >= 0.0)
{
x = floor(numer / denom);
q = (int)x;
r = (int)floor(numer - (x * denom));
}
else
{
//
// perform operations with positive values, and fix mod to make floor-based
//
x = floor(-numer / denom);
q = -(int)x;
r = (int)floor(-numer - (x * denom));
if (r != 0)
{
q--;
r = (int)denom - r;
}
}
*quotient = q;
*rem = r;
}
/*
===================
GreatestCommonDivisor
====================
*/
int GreatestCommonDivisor (int i1, int i2)
{
if (i1 > i2)
{
if (i2 == 0)
return (i1);
return GreatestCommonDivisor (i2, i1 % i2);
}
else
{
if (i1 == 0)
return (i2);
return GreatestCommonDivisor (i1, i2 % i1);
}
}
#if !id386
// TODO: move to nonintel.c
/*
===================
Invert24To16
Inverts an 8.24 value to a 16.16 value
====================
*/
fixed16_t Invert24To16(fixed16_t val)
{
if (val < 256)
return (0xFFFFFFFF);
return (fixed16_t)
(((double)0x10000 * (double)0x1000000 / (double)val) + 0.5);
}
#endif
void Mat_Mul_1x4_4x4(matrix_1x4 a,
matrix_4x4 b,
matrix_1x4 result)
{
// this function multiplies a 1x4 by a 4x4 and stores the result in a 1x4
int index_j, // column index
index_k; // row index
float sum; // temp used to hold sum of products
// loop thru columns of b
for (index_j=0; index_j<4; index_j++)
{
// multiply ith row of a by jth column of b and store the sum
// of products in the position i,j of result
sum=0;
for (index_k=0; index_k<4; index_k++)
sum+=a[index_k]*b[index_k][index_j];
// store result
result[index_j] = sum;
} // end for index_j
} // end Mat_Mul_1x4_4x4
/*
Penta: Quake doesn't have matrix muliplication!!
this function multiplies a 4x4 by a 4x4 and stores the result in a 4x4
*/
void Mat_Mul_4x4_4x4(matrix_4x4 a, matrix_4x4 b, matrix_4x4 result) {
int i, j, k;
float sum;
for (i=0; i<4; i++) {
for (j=0; j<4; j++) {
sum=0;
for (k=0; k<4; k++)
sum+=a[i][k]*b[k][j];
result[i][j] = sum;
}
}
}
/*
PENTA: Easy & fast matrix inversions with some quirks (just what Carmack likes ;) )
Thnx to http://www.cs.unc.edu/~gotz/code/affinverse.html
Find the inverse of a matrix that is made up of only scales, rotations,
and translations.
*/
void MatrixAffineInverse( matrix_4x4 m, matrix_4x4 result )
{
float Tx, Ty, Tz;
// The rotational part of the matrix is simply the transpose of the
// original matrix.
result[0][0] = m[0][0];
result[1][0] = m[0][1];
result[2][0] = m[0][2];
result[0][1] = m[1][0];
result[1][1] = m[1][1];
result[2][1] = m[1][2];
result[0][2] = m[2][0];
result[1][2] = m[2][1];
result[2][2] = m[2][2];
// The right column vector of the matrix should always be [ 0 0 0 1 ]
// In most cases. . . you don't need this column at all because it'll
// never be used in the program, but since this code is used with GL
// and it does consider this column, it is here.
result[0][3] = result[1][3] = result[2][3] = 0;
result[3][3] = 1;
// The translation components of the original matrix.
Tx = m[3][0];
Ty = m[3][1];
Tz = m[3][2];
// Rresult = -(Tm * Rm) to get the translation part of the inverse
result[3][0] = -( m[0][0] * Tx + m[0][1] * Ty + m[0][2] * Tz );
result[3][1] = -( m[1][0] * Tx + m[1][1] * Ty + m[1][2] * Tz );
result[3][2] = -( m[2][0] * Tx + m[2][1] * Ty + m[2][2] * Tz );
}
Reference in a new issue View git blame Copy permalink