/*
this c-header-file provides some useful functions and operators for 3d-graphics.
Copyright (C) 2006 Sebastian Mach ([email protected], http://greenhybrid.net)
license: GPL version 3
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License, Version 3.
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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
*/
// NOTE:
// Some of the functions are out commented since they apply
// to root-engine types (more info about the root-engine is available
// on the author's website). Nearly all functions here stem
// actually from the root-engine and have simply be partially
// rewritten (the types are for example RTVECTOR or RTPLANE, and
// the functions have the form e.g. RTPlaneNormalize).
// The evolution went further and the functions were used in
// xploderz (see website above), that's were the x stems from.
// This header has been used extensively in the gladius ray tracer
// (website above).
// You might expand and update this file and send the author a
// note.
// ONE QUICK HINT:
// FOR FULL PERFORMANCE, IN-PARAMS ARE ASSUMED TO BE NOT THE OUT PARAMS
// e.g. cross( &a, &a, &b ) might result in crap
#ifndef __xmath_h__included__
#define __xmath_h__included__
#include <memory.h>
namespace xmath
{
#define xPI05 (0.5*3.141592653589793238462643f)
#define xPI2 (2*3.141592653589793238462643f)
#define x1BYPI (1.0/3.141592653589793238462643f)
#define xPI 3.141592653589793238462643f
#define xToRadian( f ) ((f) * (xPI / 180.0f))
#define xToDegree( f ) ((f) * (180.0f / xPI))
#define SCALAR_MAX FLT_MAX
typedef float scalar;
typedef int integer;
typedef unsigned int u_integer;
#ifdef SSE
typedef struct __attribute__ ((aligned (16))) t_basevector
{
union
{
struct { scalar x,y,z; };
scalar m[3];
};
scalar _dummy;
} __attribute__ ((aligned (16))) basevector, __attribute__ ((aligned (16))) basevec;
typedef struct __attribute__ ((aligned (16))) t_vector : public t_basevector
{
t_vector( scalar a, scalar b, scalar c )
{
x = a;
y = b;
z = c;
}
t_vector() {}
operator basevector ()
{
basevector b;
b.x = x;
b.y = y;
b.z = z;
return b;
}
} __attribute__ ((aligned (16))) vector;
#else
typedef struct t_basevector
{
union
{
struct { scalar x,y,z; };
scalar m[3];
};
} basevector, basevec;
typedef struct t_vector : public t_basevector
{
t_vector( scalar a, scalar b, scalar c )
{
x = a;
y = b;
z = c;
}
t_vector() {}
operator basevector ()
{
basevector b;
b.x = x;
b.y = y;
b.z = z;
return b;
}
} vector;
#endif
typedef struct t_ray
{
vector position, direction;
} ray;
typedef struct t_plane
{
union{
struct {scalar a,b,c,d;};
scalar m[4];
basevector N;
};
} plane;
typedef struct t_sphere
{
vector position;
scalar radius;
} sphere;
typedef struct t_matrix
{
union
{
struct {
scalar _11,_12,_13,_14;
scalar _21,_22,_23,_24;
scalar _31,_32,_33,_34;
scalar _41,_42,_43,_44;
};
scalar m[4][4];
};
t_matrix() {}
} matrix, mat;
/***********************************************************/
/* VECTOR MATH */
/***********************************************************/
//{ BASIC OPS
inline vector operator + ( const vector& a, const vector& b )
{
return vector( a.x+b.x, a.y+b.y, a.z+b.z );
}
inline vector operator - ( const vector& a, const vector& b )
{
return vector( a.x-b.x, a.y-b.y, a.z-b.z );
}
inline vector operator * ( const vector& a, const vector& b )
{
return vector( a.x*b.x, a.y*b.y, a.z*b.z );
}
inline vector operator + ( const vector& a, scalar b )
{
return vector( a.x+b, a.y+b, a.z+b );
}
inline vector operator - ( const vector& a, scalar b )
{
return vector( a.x-b, a.y-b, a.z-b );
}
inline vector operator * ( const vector& a, scalar b )
{
return vector( a.x*b, a.y*b, a.z*b );
}
inline vector operator + ( scalar b, const vector& a )
{
return vector( a.x+b, a.y+b, a.z+b );
}
inline vector operator - ( scalar b, const vector& a )
{
return vector( a.x-b, a.y-b, a.z-b );
}
inline vector operator * ( scalar b, const vector& a )
{
return vector( a.x*b, a.y*b, a.z*b );
}
inline vector add ( const vector& a, const vector& b )
{
return vector( a.x+b.x, a.y+b.y, a.z+b.z );
}
inline vector sub ( const vector& a, const vector& b )
{
return vector( a.x-b.x, a.y-b.y, a.z-b.z );
}
inline vector mul ( const vector& a, const vector& b )
{
return vector( a.x*b.x, a.y*b.y, a.z*b.z );
}
inline vector add ( const vector& a, scalar b )
{
return vector( a.x+b, a.y+b, a.z+b );
}
inline vector sub ( const vector& a, scalar b )
{
return vector( a.x-b, a.y-b, a.z-b );
}
inline vector mul ( const vector& a, scalar b )
{
return vector( a.x*b, a.y*b, a.z*b );
}
inline vector add ( scalar b, const vector& a )
{
return vector( a.x+b, a.y+b, a.z+b );
}
inline vector sub ( scalar b, const vector& a )
{
return vector( a.x-b, a.y-b, a.z-b );
}
inline vector mul ( scalar b, const vector& a )
{
return vector( a.x*b, a.y*b, a.z*b );
}
inline scalar length( const vector *in )
{
return sqrtf( in->x*in->x + in->y*in->y + in->z*in->z );
}
inline scalar dot( const vector *a, const vector *b )
{
return a->x*b->x + a->y*b->y + a->z*b->z;
}
inline vector * normalize( vector *out, const vector *in )
{
register scalar l = 1.0 / sqrtf( in->x*in->x + in->y*in->y + in->z*in->z );
out->x = l * in->x;
out->y = l * in->y;
out->z = l * in->z;
return out;
}
//}
//{ ADVANCED OPS
inline vector * cross( vector *out, const vector *u, const vector *v )
{
out->x = u->y * v->z - u->z * v->y;
out->y = u->z * v->x - u->x * v->z;
out->z = u->x * v->y - u->y * v->x;
return out;
}
inline vector * madd( vector *out, const vector *base, const vector *dir, scalar f )
{
out->x = base->x + f * dir->x;
out->y = base->y + f * dir->y;
out->z = base->z + f * dir->z;
return out;
}
inline vector * subtractAndNormalize( vector *out, const vector *a, const vector *b )
{
out->x = a->x - b->x;
out->y = a->y - b->y;
out->z = a->z - b->z;
register scalar l = 1.0 / sqrtf( out->x*out->x + out->y*out->y + out->z*out->z );
out->x *= l;
out->y *= l;
out->z *= l;
return out;
}
inline vector * crossAndNormalize( vector *out, const vector *u, const vector *v )
{
out->x = u->y * v->z - u->z * v->y;
out->y = u->z * v->x - u->x * v->z;
out->z = u->x * v->y - u->y * v->x;
register scalar l = 1.0 / sqrtf( out->x*out->x + out->y*out->y + out->z*out->z );
out->x *= l;
out->y *= l;
out->z *= l;
return out;
}
inline void reflect( vector *pOut, const vector *n, const vector *incident)
{
scalar d = 2 * dot(n,incident);
pOut->x = incident->x - d * n->x;
pOut->y = incident->y - d * n->y;
pOut->z = incident->z - d * n->z;
}
inline void refract( vector *pOut, const vector *normal, const vector *incident, scalar n1, scalar n2 )
{
/**pOut = *incident;
return;*/
scalar r = n1 / n2;
scalar w = -dot(incident,normal) * r;
scalar k = 1 + (w-r)*(w+r);
if( k<0 ){
reflect( pOut, normal, incident );
return;
}
scalar f = w-sqrt(k);
pOut->x = r*incident->x+f*normal->x;
pOut->y = r*incident->y+f*normal->y;
pOut->z = r*incident->z+f*normal->z;
}
/*inline vector refract( const vector &normal, const vector &incident, scalar n1, scalar n2 )
{
scalar r = n1 / n2;
scalar w = -dot(incident,normal) * r;
scalar k = 1 + (w-r)*(w+r);
if( k<0 ) return vector::invalid();
return r*incident+(w-sqrt(k))*normal;
}*/
inline scalar schlick2( const vector *normal, const vector *incident, scalar n1, scalar n2 )
{
scalar r0 = ( n1 - n2 ) / ( n1 + n2 );
r0 *= r0;
scalar cosX = -dot( normal, incident );
if( n1 > n2 ){
scalar n = n1 / n2;
scalar sinT2 = n * n * ( 1.0 - cosX * cosX );
if( sinT2 > 1.0 ) return 1.0;
cosX = sqrt( 1.0 - sinT2 );
}
scalar x = 1.0 - cosX;
return r0 + ( 1.0 - r0 ) * x * x * x * x * x;
}
inline integer intersectRayBox( scalar *pt_min, scalar *pt_max,
const ray *pRay,
const vector *bbMin, const vector *bbMax )
{
const vector n = pRay->direction;
const vector pos = pRay->position;
vector Min = *bbMin, Max = *bbMax;
scalar tmin, tmax, tymin, tymax, tzmin, tzmax;
scalar inx = 1.0 / n.x;
if( n.x >= 0 )
{
tmin = ( Min.x - pos.x ) * inx;
tmax = ( Max.x - pos.x ) * inx;
}
else
{
tmin = ( Max.x - pos.x ) * inx;
tmax = ( Min.x - pos.x ) * inx;
}
scalar iny = 1.0 / n.y;
if( n.y >= 0 )
{
tymin = ( Min.y - pos.y ) * iny;
tymax = ( Max.y - pos.y ) * iny;
}
else
{
tymin = ( Max.y - pos.y ) * iny;
tymax = ( Min.y - pos.y ) * iny;
}
if( ( tmin > tymax ) || ( tymin > tmax ) )
return 0;
if( tymin > tmin )
tmin = tymin;
if( tymax < tmax )
tmax = tymax;
scalar inz = 1.0 / n.z;
if( n.z >= 0 )
{
tzmin = ( Min.z - pos.z ) * inz;
tzmax = ( Max.z - pos.z ) * inz;
}
else
{
tzmin = ( Max.z - pos.z ) * inz;
tzmax = ( Min.z - pos.z ) * inz;
}
if( ( tmin > tzmax ) || ( tzmin > tmax ) )
return 0;
if( tzmin > tmin )
tmin = tzmin;
if( tzmax < tmax )
tmax = tzmax;
*pt_min = tmin;
*pt_max = tmax;
return 1;
}
//}
/***********************************************************/
/* PLANE MATH */
/***********************************************************/
//{ BASIC OPS
inline plane *planeNormalize( plane *out, const plane *in )
{
scalar _d = (1.0f) / sqrtf( in->a*in->a + in->b*in->b + in->c*in->c );
out->a = in->a * _d;
out->b = in->b * _d;
out->c = in->c * _d;
out->d = in->d * _d;
return out;
}
inline scalar planeDotCoord( const plane *Plane, const vector *Coord )
{
return Plane->a*Coord->x + Plane->b*Coord->y + Plane->c*Coord->z + Plane->d;
}
inline scalar planeDotNormal( const plane *Plane, const vector *Coord )
{
return Plane->a*Coord->x + Plane->b*Coord->y + Plane->c*Coord->z;
}
inline plane *planeFromPointNormal( plane *out, const vector *Pos, const vector *N )
{
out->a = N->x;
out->b = N->y;
out->c = N->z;
out->d = -Pos->x*N->x - Pos->y*N->y - Pos->z*N->z;// == -RTVectorDot( Pos, N );
return out;
}
inline plane *planeFromPoints( plane *Plane, const vector *pA, const vector *pB, const vector *pC )
{
vector U, V;
U = sub( *pA, *pB );
V = sub( *pA, *pC );
cross( (vector*)Plane->m, &U, &V );
Plane->d = -pA->x*Plane->a - pA->y*Plane->b - pA->z*Plane->c;// == -RTVectorDot( pA, Plane );
planeNormalize( Plane, Plane );
return Plane;
}
inline scalar planeIntersectLine( const plane *Plane, const vector *pA, const vector *pB )
{
scalar a = Plane->a*pA->x + Plane->b*pA->y + Plane->c*pA->z + Plane->d;
scalar b = Plane->a*pB->x + Plane->b*pB->y + Plane->c*pB->z + Plane->d;
return a / (a - b);
}
inline scalar planeIntersectRay( const plane *Plane, const t_ray *ray )
{
/*
scalar a = Plane->a*ray->position.x + Plane->b*ray->position.y + Plane->c*ray->position.z + Plane->d;
scalar b = Plane->a*(ray->position.x+ray->direction.x)
+ Plane->b*(ray->position.y+ray->direction.y)
+ Plane->c*(ray->position.z+ray->direction.z)
+ Plane->d;
return a / (a - b);
*/
return -( Plane->a*ray->position.x + Plane->b*ray->position.y + Plane->c*ray->position.z + Plane->d )
/ ( Plane->a*ray->direction.x + Plane->b*ray->direction.y + Plane->c*ray->direction.z );
//(AX0 + BY0 + CZ0 + D) / (AXd + BYd + CZd)
}
/*
RTPLANE *RTPlaneTransform( RTPLANE *pOut, const RTPLANE *pIn, const RTMATRIX *mat )
{
RTVECTOR Src;
Src.x = pIn->a; Src.y = pIn->b; Src.z = pIn->c;
pOut->a = Src.x*mat->_11 + Src.y*mat->_21 + Src.z*mat->_31 + pIn->d*mat->_41;
pOut->b = Src.x*mat->_12 + Src.y*mat->_22 + Src.z*mat->_32 + pIn->d*mat->_42;
pOut->c = Src.x*mat->_13 + Src.y*mat->_23 + Src.z*mat->_33 + pIn->d*mat->_43;
pOut->d = Src.x*mat->_14 + Src.y*mat->_24 + Src.z*mat->_34 + pIn->d*mat->_44;
return pOut;
}*/
//}
/***********************************************************/
/* MATRIX MATH */
/***********************************************************/
//{ BASIC OPS
inline matrix *set( matrix *m,
scalar _11, scalar _12, scalar _13, scalar _14,
scalar _21, scalar _22, scalar _23, scalar _24,
scalar _31, scalar _32, scalar _33, scalar _34,
scalar _41, scalar _42, scalar _43, scalar _44 )
{
m->_11 = _11; m->_12 = _12; m->_13 = _13; m->_14 = _14;
m->_21 = _21; m->_22 = _22; m->_23 = _23; m->_24 = _24;
m->_31 = _31; m->_32 = _32; m->_33 = _33; m->_34 = _34;
m->_41 = _41; m->_42 = _42; m->_43 = _43; m->_44 = _44;
return m;
}
inline matrix *identity( matrix *m )
{
m->_11 = 1; m->_12 = 0; m->_13 = 0; m->_14 = 0;
m->_21 = 0; m->_22 = 1; m->_23 = 0; m->_24 = 0;
m->_31 = 0; m->_32 = 0; m->_33 = 1; m->_34 = 0;
m->_41 = 0; m->_42 = 0; m->_43 = 0; m->_44 = 1;
return m;
}
inline matrix *scaling( matrix *m, scalar s )
{
m->_11 = s; m->_12 = 0; m->_13 = 0; m->_14 = 0;
m->_21 = 0; m->_22 = s; m->_23 = 0; m->_24 = 0;
m->_31 = 0; m->_32 = 0; m->_33 = s; m->_34 = 0;
m->_41 = 0; m->_42 = 0; m->_43 = 0; m->_44 = 1;
return m;
}
inline matrix *scaling( matrix *m, scalar x, scalar y, scalar z )
{
m->_11 = x; m->_12 = 0; m->_13 = 0; m->_14 = 0;
m->_21 = 0; m->_22 = y; m->_23 = 0; m->_24 = 0;
m->_31 = 0; m->_32 = 0; m->_33 = z; m->_34 = 0;
m->_41 = 0; m->_42 = 0; m->_43 = 0; m->_44 = 1;
return m;
}
inline matrix *translation( matrix *m, scalar x, scalar y, scalar z )
{
m->_11 = 1; m->_12 = 0; m->_13 = 0; m->_14 = 0;
m->_21 = 0; m->_22 = 1; m->_23 = 0; m->_24 = 0;
m->_31 = 0; m->_32 = 0; m->_33 = 1; m->_34 = 0;
m->_41 = x; m->_42 = y; m->_43 = z; m->_44 = 1;
return m;
}
inline matrix *rotationX( matrix *m, scalar angle )
{
register scalar sina = sinf(angle);
register scalar cosa = cosf(angle);
m->_11 = 1; m->_12 = 0; m->_13 = 0; m->_14 = 0;
m->_21 = 0; m->_22 = cosa; m->_23 = sina; m->_24 = 0;
m->_31 = 0; m->_32 =-sina; m->_33 = cosa; m->_34 = 0;
m->_41 = 0; m->_42 = 0; m->_43 = 0; m->_44 = 1;
return m;
}
inline matrix *rotationY( matrix *m, scalar angle )
{
register scalar sina = sinf(angle);
register scalar cosa = cosf(angle);
m->_11 = cosa; m->_12 = 0; m->_13 =-sina; m->_14 = 0;
m->_21 = 0; m->_22 = 1; m->_23 = 0; m->_24 = 0;
m->_31 = sina; m->_32 = 0; m->_33 = cosa; m->_34 = 0;
m->_41 = 0; m->_42 = 0; m->_43 = 0; m->_44 = 1;
return m;
}
inline void transform( vector *pOut, const matrix *mat, const vector *vec )
{
pOut->x = vec->x*mat->_11 + vec->y*mat->_21 + vec->z*mat->_31 + mat->_41;
pOut->y = vec->x*mat->_12 + vec->y*mat->_22 + vec->z*mat->_32 + mat->_42;
pOut->z = vec->x*mat->_13 + vec->y*mat->_23 + vec->z*mat->_33 + mat->_43;
}
inline matrix *rotationZ( matrix *m, scalar angle )
{
register scalar sina = sinf(angle);
register scalar cosa = cosf(angle);
m->_11 = cosa; m->_12 = sina; m->_13 = 0; m->_14 = 0;
m->_21 =-sina; m->_22 = cosa; m->_23 = 0; m->_24 = 0;
m->_31 = 0; m->_32 = 0; m->_33 = 1; m->_34 = 0;
m->_41 = 0; m->_42 = 0; m->_43 = 0; m->_44 = 1;
return m;
}
inline matrix *multiply( matrix *res, const matrix *a, const matrix *b )
{
register integer u,v;
for( v=0; v<4; v++ ) for( u=0; u<4; u++ )
res->m[u][v] = a->m[u][0]*b->m[0][v] + a->m[u][1]*b->m[1][v] + a->m[u][2]*b->m[2][v] + a->m[u][3]*b->m[3][v];
return res;
}
inline matrix operator * ( const matrix &a, const matrix &b )
{
matrix res;
register integer u,v;
for( v=0; v<4; v++ ) for( u=0; u<4; u++ )
res.m[u][v] = a.m[u][0]*b.m[0][v] + a.m[u][1]*b.m[1][v] + a.m[u][2]*b.m[2][v] + a.m[u][3]*b.m[3][v];
return res;
}
//}
//{ ADVANCED OPS
inline matrix *matrixPerspectiveFovLH( matrix *out, scalar fov, scalar aspectRatio, scalar zn, scalar zf )
{
scalar h = 1.0 / tan( fov*0.5 );
scalar w = h / aspectRatio;
out->_11=w; out->_12=0; out->_13=0; out->_14=0;
out->_21=0; out->_22=h; out->_23=0; out->_24=0;
out->_31=0; out->_32=0; out->_33=zf/(zf-zn); out->_34=1;
out->_41=0; out->_42=0; out->_43=-zn*zf/(zf-zn); out->_44=0;
return out;
}
inline matrix *matrixLookAtLH( matrix *m, const vector *eye, const vector *at, const vector *up )
{
vector X,Y,Z;
subtractAndNormalize( &Z, at, eye );
crossAndNormalize( &X, up, &Z );
cross( &Y, &Z, &X );
m->_11 = X.x; m->_12 = Y.x; m->_13 = Z.x; m->_14 = 0;
m->_21 = X.y; m->_22 = Y.y; m->_23 = Z.y; m->_24 = 0;
m->_31 = X.z; m->_32 = Y.z; m->_33 = Z.z; m->_34 = 0;
m->_44 = 1;
m->_41 = -dot(&X,eye);
m->_42 = -dot(&Y,eye);
m->_43 = -dot(&Z,eye);
return m;
}
inline scalar sarrus(
scalar _11, scalar _12, scalar _13,
scalar _21, scalar _22, scalar _23,
scalar _31, scalar _32, scalar _33
)
{
return _11*_22*_33 - _11*_23*_32
+ _12*_23*_31 - _12*_21*_33
+ _13*_21*_32 - _13*_22*_31;
}
inline scalar subdet33( const matrix *mat, int v, int u )
{
int u1, u2, u3, v1, v2, v3;
switch( u )
{
case 0: u1=1; u2=2; u3=3; break;
case 1: u1=0; u2=2; u3=3; break;
case 2: u1=0; u2=1; u3=3; break;
default: u1=0; u2=1; u3=2; break;
};
switch( v )
{
case 0: v1=1; v2=2; v3=3; break;
case 1: v1=0; v2=2; v3=3; break;
case 2: v1=0; v2=1; v3=3; break;
default: v1=0; v2=1; v3=2; break;
};
return sarrus( mat->m[v1][u1], mat->m[v1][u2], mat->m[v1][u3],
mat->m[v2][u1], mat->m[v2][u2], mat->m[v2][u3],
mat->m[v3][u1], mat->m[v3][u2], mat->m[v3][u3] );
}
inline scalar determinant( const matrix *in )
{
scalar det11 = in->_11 * sarrus( in->_22, in->_23, in->_24,
in->_32, in->_33, in->_34,
in->_42, in->_43, in->_44
);
scalar det12 = in->_12 * sarrus( in->_21, in->_23, in->_24,
in->_31, in->_33, in->_34,
in->_41, in->_43, in->_44
);
scalar det13 = in->_13 * sarrus( in->_21, in->_22, in->_24,
in->_31, in->_32, in->_34,
in->_41, in->_42, in->_44
);
scalar det14 = in->_14 * sarrus( in->_21, in->_22, in->_23,
in->_31, in->_33, in->_33,
in->_41, in->_42, in->_43
);
return ( +det11 -det12 +det13 -det14 );
}
inline matrix* inverse( matrix *out, scalar *pDet, const matrix *in )
{
matrix tmp;
memcpy( &tmp, in, sizeof( matrix ) );
scalar det = determinant( &tmp );
if( pDet )
*pDet = det;
scalar detinv = 1.0 / det;
int u,v;
for( v=0; v<4; v++ )
for( u=0; u<4; u++ )
out->m[v][u] = detinv * powf(-1,(scalar)(u+v)) * subdet33( &tmp, u, v );
return out;
}
inline matrix* transpose( matrix *out, const matrix *in )
{
matrix Tmp;
int y,x;
for( y=0; y<4; y++ )
for( x=0; x<4; x++ )
Tmp.m[x][y] = in->m[y][x];
memcpy( out, &Tmp, sizeof( matrix ) );
return out;
}
//}
/***********************************************************/
/* HELPER FUNCTIONS */
/***********************************************************/
//{
inline int intersectTriangleBox(
const t_vector *a, const t_vector *b, const t_vector *c,
const t_vector *bbMin, const t_vector *bbMax )
{
unsigned int u;
for( u=0; u<3; u++ )
{
if( ( a->m[u] < bbMin->m[u] && b->m[u] < bbMin->m[u] && c->m[u] < bbMin->m[u] )
|| ( a->m[u] > bbMin->m[u] && b->m[u] > bbMin->m[u] && c->m[u] < bbMin->m[u] )
)
return 0;
}
return 1;
}
inline int intersectSphereBox( const t_vector *sphereCenter, scalar sphereRadius, const t_vector *bbMin, const t_vector *bbMax )
{
// from Jim Arvo/graphics gems
scalar r2 = sphereRadius*sphereRadius;
scalar dmin = 0;
/*for( i = 0; i < n; i++ ) {
if( C[i] < Bmin[i] ) dmin += SQR(C[i] - Bmin[i] ); else
if( C[i] > Bmax[i] ) dmin += SQR( C[i] - Bmax[i] );
}*/
if( sphereCenter->x < bbMin->x )
dmin += ( sphereCenter->x - bbMin->x ) * ( sphereCenter->x - bbMin->x );
else if( sphereCenter->x > bbMax->x )
dmin += ( sphereCenter->x - bbMax->x ) * ( sphereCenter->x - bbMax->x );
if( sphereCenter->y < bbMin->y )
dmin += ( sphereCenter->y - bbMin->y ) * ( sphereCenter->y - bbMin->y );
else if( sphereCenter->y > bbMax->y )
dmin += ( sphereCenter->y - bbMax->y ) * ( sphereCenter->y - bbMax->y );
if( sphereCenter->z < bbMin->z )
dmin += ( sphereCenter->z - bbMin->z ) * ( sphereCenter->z - bbMin->z );
else if( sphereCenter->z > bbMax->z )
dmin += ( sphereCenter->z - bbMax->z ) * ( sphereCenter->z - bbMax->z );
if( dmin <= r2 )
return 1;
return 0;
}
//}
//}
}//namespace xmath
#endif