/****************************************************************************
* quadrics.c
*
* This module implements the code for the quadric shape primitive.
*
* from Persistence of Vision(tm) Ray Tracer
* Copyright 1996 Persistence of Vision Team
*---------------------------------------------------------------------------
* NOTICE: This source code file is provided so that users may experiment
* with enhancements to POV-Ray and to port the software to platforms other
* than those supported by the POV-Ray Team. There are strict rules under
* which you are permitted to use this file. The rules are in the file
* named POVLEGAL.DOC which should be distributed with this file. If
* POVLEGAL.DOC is not available or for more info please contact the POV-Ray
* Team Coordinator by leaving a message in CompuServe's Graphics Developer's
* Forum. The latest version of POV-Ray may be found there as well.
*
* This program is based on the popular DKB raytracer version 2.12.
* DKBTrace was originally written by David K. Buck.
* DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
*
*****************************************************************************/
#include "frame.h"
#include "povray.h"
#include "vector.h"
#include "povproto.h"
#include "bbox.h"
#include "objects.h"
#include "matrices.h"
#include "planes.h"
#include "quadrics.h"
/*****************************************************************************
* Local preprocessor defines
******************************************************************************/
/* The following defines make typing much easier! [DB 7/94] */
#define Xd Ray->Direction[X]
#define Yd Ray->Direction[Y]
#define Zd Ray->Direction[Z]
#define Xo Ray->Initial[X]
#define Yo Ray->Initial[Y]
#define Zo Ray->Initial[Z]
#define QA Quadric->Square_Terms[X]
#define QE Quadric->Square_Terms[Y]
#define QH Quadric->Square_Terms[Z]
#define QB Quadric->Mixed_Terms[X]
#define QC Quadric->Mixed_Terms[Y]
#define QF Quadric->Mixed_Terms[Z]
#define QD Quadric->Terms[X]
#define QG Quadric->Terms[Y]
#define QI Quadric->Terms[Z]
#define QJ Quadric->Constant
/*****************************************************************************
* Static functions
******************************************************************************/
static int Intersect_Quadric PARAMS((RAY *Ray, QUADRIC *Quadric, DBL *Depth1, DBL *Depth2));
static void Quadric_To_Matrix PARAMS((QUADRIC *Quadric, MATRIX Matrix));
static void Matrix_To_Quadric PARAMS((MATRIX Matrix, QUADRIC *Quadric));
static int All_Quadric_Intersections PARAMS((OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack));
static int Inside_Quadric PARAMS((VECTOR IPoint, OBJECT *Object));
static void Quadric_Normal PARAMS((VECTOR Result, OBJECT *Object, INTERSECTION *Inter));
static void *Copy_Quadric PARAMS((OBJECT *Object));
static void Translate_Quadric PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
static void Rotate_Quadric PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
static void Scale_Quadric PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
static void Transform_Quadric PARAMS((OBJECT *Object, TRANSFORM *Trans));
static void Invert_Quadric PARAMS((OBJECT *Object));
static void Destroy_Quadric PARAMS((OBJECT *Object));
/*****************************************************************************
* Local variables
******************************************************************************/
METHODS Quadric_Methods =
{
All_Quadric_Intersections,
Inside_Quadric, Quadric_Normal,
Copy_Quadric,
Translate_Quadric, Rotate_Quadric,
Scale_Quadric, Transform_Quadric, Invert_Quadric,
Destroy_Quadric
};
/*****************************************************************************
*
* FUNCTION
*
* All_Quadric_Intersections
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* POV-Ray Team
*
* DESCRIPTION
*
* -
*
* CHANGES
*
* -
*
******************************************************************************/
static int All_Quadric_Intersections(Object, Ray, Depth_Stack)
OBJECT *Object;
RAY *Ray;
ISTACK *Depth_Stack;
{
DBL Depth1, Depth2;
VECTOR IPoint;
register int Intersection_Found;
Intersection_Found = FALSE;
if (Intersect_Quadric(Ray, (QUADRIC *)Object, &Depth1, &Depth2))
{
if ((Depth1 > Small_Tolerance) && (Depth1 < Max_Distance))
{
VEvaluateRay(IPoint, Ray->Initial, Depth1, Ray->Direction);
if (Point_In_Clip(IPoint, Object->Clip))
{
push_entry(Depth1, IPoint, Object, Depth_Stack);
Intersection_Found = TRUE;
}
}
if ((Depth2 > Small_Tolerance) && (Depth2 < Max_Distance))
{
VEvaluateRay(IPoint, Ray->Initial, Depth2, Ray->Direction);
if (Point_In_Clip(IPoint, Object->Clip))
{
push_entry(Depth2, IPoint, Object, Depth_Stack);
Intersection_Found = TRUE;
}
}
}
return(Intersection_Found);
}
/*****************************************************************************
*
* FUNCTION
*
* Intersect_Quadric
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* POV-Ray Team
*
* DESCRIPTION
*
* -
*
* CHANGES
*
* -
*
******************************************************************************/
static int Intersect_Quadric(Ray, Quadric, Depth1, Depth2)
RAY *Ray;
QUADRIC *Quadric;
DBL *Depth1, *Depth2;
{
register DBL a, b, c, d;
Increase_Counter(stats[Ray_Quadric_Tests]);
a = Xd * (QA * Xd + QB * Yd + QC * Zd) +
Yd * (QE * Yd + QF * Zd) +
QH * Zd * Zd;
b = Xd * (QA * Xo + 0.5 * (QB * Yo + QC * Zo + QD)) +
Yd * (QE * Yo + 0.5 * (QB * Xo + QF * Zo + QG)) +
Zd * (QH * Zo + 0.5 * (QC * Xo + QF * Yo + QI));
c = Xo * (QA * Xo + QB * Yo + QC * Zo + QD) +
Yo * (QE * Yo + QF * Zo + QG) +
Zo * (QH * Zo + QI) +
QJ;
if (a != 0.0)
{
/* The equation is quadratic - find its roots */
d = Sqr(b) - a * c;
if (d <= 0.0)
{
return(FALSE);
}
d = sqrt (d);
*Depth1 = (-b + d) / (a);
*Depth2 = (-b - d) / (a);
}
else
{
/* There are no quadratic terms. Solve the linear equation instead. */
if (b == 0.0)
{
return(FALSE);
}
*Depth1 = - 0.5 * c / b;
*Depth2 = Max_Distance;
}
Increase_Counter(stats[Ray_Quadric_Tests_Succeeded]);
return(TRUE);
}
/*****************************************************************************
*
* FUNCTION
*
* Inside_Quadric
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* POV-Ray Team
*
* DESCRIPTION
*
* -
*
* CHANGES
*
* -
*
******************************************************************************/
static int Inside_Quadric(IPoint, Object)
VECTOR IPoint;
OBJECT *Object;
{
QUADRIC *Quadric = (QUADRIC *)Object;
/* This is faster and shorter. [DB 7/94] */
return((IPoint[X] * (QA * IPoint[X] + QB * IPoint[Y] + QD) +
IPoint[Y] * (QE * IPoint[Y] + QF * IPoint[Z] + QG) +
IPoint[Z] * (QH * IPoint[Z] + QC * IPoint[X] + QI) + QJ) <= 0.0);
}
/*****************************************************************************
*
* FUNCTION
*
* Quadric_Normal
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* POV-Ray Team
*
* DESCRIPTION
*
* -
*
* CHANGES
*
* -
*
******************************************************************************/
static void Quadric_Normal(Result, Object, Inter)
OBJECT *Object;
VECTOR Result;
INTERSECTION *Inter;
{
QUADRIC *Quadric = (QUADRIC *) Object;
DBL Len;
/* This is faster and shorter. [DB 7/94] */
Result[X] = 2.0 * QA * Inter->IPoint[X] +
QB * Inter->IPoint[Y] +
QC * Inter->IPoint[Z] +
QD;
Result[Y] = QB * Inter->IPoint[X] +
2.0 * QE * Inter->IPoint[Y] +
QF * Inter->IPoint[Z] +
QG;
Result[Z] = QC * Inter->IPoint[X] +
QF * Inter->IPoint[Y] +
2.0 * QH * Inter->IPoint[Z] +
QI;
VLength(Len, Result);
if (Len == 0.0)
{
/* The normal is not defined at this point of the surface. */
/* Set it to any arbitrary direction. */
Make_Vector(Result, 1.0, 0.0, 0.0);
}
else
{
VInverseScaleEq(Result, Len);
}
}
/*****************************************************************************
*
* FUNCTION
*
* Transform_Quadric
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* POV-Ray Team
*
* DESCRIPTION
*
* -
*
* CHANGES
*
* -
*
******************************************************************************/
static void Transform_Quadric(Object, Trans)
OBJECT *Object;
TRANSFORM *Trans;
{
QUADRIC *Quadric=(QUADRIC *)Object;
MATRIX Quadric_Matrix, Transform_Transposed;
Quadric_To_Matrix (Quadric, Quadric_Matrix);
MTimes (Quadric_Matrix, Trans->inverse, Quadric_Matrix);
MTranspose (Transform_Transposed, Trans->inverse);
MTimes (Quadric_Matrix, Quadric_Matrix, Transform_Transposed);
Matrix_To_Quadric (Quadric_Matrix, Quadric);
Recompute_BBox(&Object->BBox, Trans);
}
/*****************************************************************************
*
* FUNCTION
*
* Quadric_To_Matrix
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* POV-Ray Team
*
* DESCRIPTION
*
* -
*
* CHANGES
*
* -
*
******************************************************************************/
static void Quadric_To_Matrix(Quadric, Matrix)
QUADRIC *Quadric;
MATRIX Matrix;
{
MZero (Matrix);
Matrix[0][0] = Quadric->Square_Terms[X];
Matrix[1][1] = Quadric->Square_Terms[Y];
Matrix[2][2] = Quadric->Square_Terms[Z];
Matrix[0][1] = Quadric->Mixed_Terms[X];
Matrix[0][2] = Quadric->Mixed_Terms[Y];
Matrix[0][3] = Quadric->Terms[X];
Matrix[1][2] = Quadric->Mixed_Terms[Z];
Matrix[1][3] = Quadric->Terms[Y];
Matrix[2][3] = Quadric->Terms[Z];
Matrix[3][3] = Quadric->Constant;
}
/*****************************************************************************
*
* FUNCTION
*
* Matrix_To_Quadric
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* POV-Ray Team
*
* DESCRIPTION
*
* -
*
* CHANGES
*
* -
*
******************************************************************************/
static void Matrix_To_Quadric(Matrix, Quadric)
MATRIX Matrix;
QUADRIC *Quadric;
{
Quadric->Square_Terms[X] = Matrix[0][0];
Quadric->Square_Terms[Y] = Matrix[1][1];
Quadric->Square_Terms[Z] = Matrix[2][2];
Quadric->Mixed_Terms[X] = Matrix[0][1] + Matrix[1][0];
Quadric->Mixed_Terms[Y] = Matrix[0][2] + Matrix[2][0];
Quadric->Mixed_Terms[Z] = Matrix[1][2] + Matrix[2][1];
Quadric->Terms[X] = Matrix[0][3] + Matrix[3][0];
Quadric->Terms[Y] = Matrix[1][3] + Matrix[3][1];
Quadric->Terms[Z] = Matrix[2][3] + Matrix[3][2];
Quadric->Constant = Matrix[3][3];
}
/*****************************************************************************
*
* FUNCTION
*
* Translate_Quadric
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* POV-Ray Team
*
* DESCRIPTION
*
* -
*
* CHANGES
*
* -
*
******************************************************************************/
static void Translate_Quadric(Object, Vector, Trans)
OBJECT *Object;
VECTOR Vector;
TRANSFORM *Trans;
{
Transform_Quadric (Object, Trans);
}
/*****************************************************************************
*
* FUNCTION
*
* Rotate_Quadric
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* POV-Ray Team
*
* DESCRIPTION
*
* -
*
* CHANGES
*
* -
*
******************************************************************************/
static void Rotate_Quadric(Object, Vector, Trans)
OBJECT *Object;
VECTOR Vector;
TRANSFORM *Trans;
{
Transform_Quadric (Object, Trans);
}
/*****************************************************************************
*
* FUNCTION
*
* Scale_Quadric
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* POV-Ray Team
*
* DESCRIPTION
*
* -
*
* CHANGES
*
* -
*
******************************************************************************/
static void Scale_Quadric(Object, Vector, Trans)
OBJECT *Object;
VECTOR Vector;
TRANSFORM *Trans;
{
Transform_Quadric (Object, Trans);
}
/*****************************************************************************
*
* FUNCTION
*
* Invert_Quadric
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* POV-Ray Team
*
* DESCRIPTION
*
* -
*
* CHANGES
*
* -
*
******************************************************************************/
static void Invert_Quadric(Object)
OBJECT *Object;
{
QUADRIC *Quadric = (QUADRIC *) Object;
VScaleEq(Quadric->Square_Terms, -1.0);
VScaleEq(Quadric->Mixed_Terms, -1.0);
VScaleEq(Quadric->Terms, -1.0);
Quadric->Constant *= -1.0;
Invert_Flag(Object, INVERTED_FLAG);
}
/*****************************************************************************
*
* FUNCTION
*
* Create_Quadric
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* POV-Ray Team
*
* DESCRIPTION
*
* -
*
* CHANGES
*
* -
*
******************************************************************************/
QUADRIC *Create_Quadric()
{
QUADRIC *New;
New = (QUADRIC *)POV_MALLOC(sizeof (QUADRIC), "quadric");
INIT_OBJECT_FIELDS(New, QUADRIC_OBJECT, &Quadric_Methods)
Make_Vector (New->Square_Terms, 1.0, 1.0, 1.0);
Make_Vector (New->Mixed_Terms, 0.0, 0.0, 0.0);
Make_Vector (New->Terms, 0.0, 0.0, 0.0);
New->Constant = 1.0;
return(New);
}
/*****************************************************************************
*
* FUNCTION
*
* Copy_Quadric
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* POV-Ray Team
*
* DESCRIPTION
*
* -
*
* CHANGES
*
* -
*
******************************************************************************/
static void *Copy_Quadric(Object)
OBJECT *Object;
{
QUADRIC *New;
New = Create_Quadric();
*New = *((QUADRIC *)Object);
return(New);
}
/*****************************************************************************
*
* FUNCTION
*
* Destroy_Quadric
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
* POV-Ray Team
*
* DESCRIPTION
*
* -
*
* CHANGES
*
* -
*
******************************************************************************/
static void Destroy_Quadric(Object)
OBJECT *Object;
{
POV_FREE (Object);
}
/*****************************************************************************
*
* FUNCTION
*
* Compute_Quadric_BBox
*
* INPUT
*
* Quadric - Qaudric object
*
* OUTPUT
*
* Quadric
*
* RETURNS
*
* AUTHOR
*
* Dieter Bayer
*
* DESCRIPTION
*
* Compute a bounding box around a quadric.
*
* This function calculates the bounding box of a quadric given in
* its normal form, i.e. f(x,y,z) = A*x^2 + E*y^2 + H*z^2 + J = 0.
*
* NOTE: Translated quadrics can also be bounded by this function.
*
* Supported: cones, cylinders, ellipsoids, hyperboloids.
*
* CHANGES
*
* May 1994 : Creation.
*
* Sep 1994 : Added support of hyperpoloids. Improved bounding of
* quadrics used in CSG intersections. [DB]
*
******************************************************************************/
void Compute_Quadric_BBox(Quadric, ClipMin, ClipMax)
QUADRIC *Quadric;
VECTOR ClipMin, ClipMax;
{
DBL A, B, C, D, E, F, G, H, I, J;
DBL a, b, c, d;
DBL rx, ry, rz, rx1, rx2, ry1, ry2, rz1, rz2, x, y, z;
DBL New_Volume, Old_Volume;
VECTOR Min, Max, TmpMin, TmpMax, NewMin, NewMax, T1;
BBOX Old;
OBJECT *Sib;
/*
* Check for 'normal' form. If the quadric isn't in it's normal form
* we can't do anything (we could, but that would be to tedious!
* Diagonalising the quadric's 4x4 matrix, i.e. finding its eigenvalues
* and eigenvectors -> solving a 4th order polynom).
*/
/* Get quadrics coefficients. */
A = Quadric->Square_Terms[X];
E = Quadric->Square_Terms[Y];
H = Quadric->Square_Terms[Z];
B = Quadric->Mixed_Terms[X] / 2.0;
C = Quadric->Mixed_Terms[Y] / 2.0;
F = Quadric->Mixed_Terms[Z] / 2.0;
D = Quadric->Terms[X] / 2.0;
G = Quadric->Terms[Y] / 2.0;
I = Quadric->Terms[Z] / 2.0;
J = Quadric->Constant;
/* Set small values to 0. */
if (fabs(A) < EPSILON) A = 0.0;
if (fabs(B) < EPSILON) B = 0.0;
if (fabs(C) < EPSILON) C = 0.0;
if (fabs(D) < EPSILON) D = 0.0;
if (fabs(E) < EPSILON) E = 0.0;
if (fabs(F) < EPSILON) F = 0.0;
if (fabs(G) < EPSILON) G = 0.0;
if (fabs(H) < EPSILON) H = 0.0;
if (fabs(I) < EPSILON) I = 0.0;
if (fabs(J) < EPSILON) J = 0.0;
/* Non-zero mixed terms --> return */
if ((B != 0.0) || (C != 0.0) || (F != 0.0))
{
return;
}
/* Non-zero linear terms --> get translation vector */
if ((D != 0.0) || (G != 0.0) || (I != 0.0))
{
if (A != 0.0)
{
T1[X] = -D / A;
}
else
{
if (D != 0.0)
{
T1[X] = J / (2.0 * D);
}
else
{
T1[X] = 0.0;
}
}
if (E != 0.0)
{
T1[Y] = -G / E;
}
else
{
if (G != 0.0)
{
T1[Y] = J / (2.0 * G);
}
else
{
T1[Y] = 0.0;
}
}
if (H != 0.0)
{
T1[Z] = -I / H;
}
else
{
if (I != 0.0)
{
T1[Z] = J / (2.0 * I);
}
else
{
T1[Z] = 0.0;
}
}
/* Recalculate coefficients. */
D += A * T1[X];
G += E * T1[Y];
I += H * T1[Z];
J -= T1[X]*(A*T1[X] + 2.0*D) + T1[Y]*(E*T1[Y] + 2.0*G) + T1[Z]*(H*T1[Z] + 2.0*I);
}
else
{
Make_Vector(T1, 0.0, 0.0, 0.0);
}
/* Get old bounding box. */
Old = Quadric->BBox;
/* Init new bounding box. */
NewMin[X] = NewMin[Y] = NewMin[Z] = -BOUND_HUGE/2;
NewMax[X] = NewMax[Y] = NewMax[Z] = BOUND_HUGE/2;
/* Get the bounding box of the clipping object. */
if (Quadric->Clip != NULL)
{
Min[X] = Min[Y] = Min[Z] = -BOUND_HUGE;
Max[X] = Max[Y] = Max[Z] = BOUND_HUGE;
/* Intersect the members bounding boxes. */
for (Sib = Quadric->Clip; Sib != NULL; Sib = Sib->Sibling)
{
if (!Test_Flag(Sib, INVERTED_FLAG))
{
if (Sib->Methods == &Plane_Methods)
{
Compute_Plane_Min_Max((PLANE *)Sib, TmpMin, TmpMax);
}
else
{
Make_min_max_from_BBox(TmpMin, TmpMax, Sib->BBox);
}
Min[X] = max(Min[X], TmpMin[X]);
Min[Y] = max(Min[Y], TmpMin[Y]);
Min[Z] = max(Min[Z], TmpMin[Z]);
Max[X] = min(Max[X], TmpMax[X]);
Max[Y] = min(Max[Y], TmpMax[Y]);
Max[Z] = min(Max[Z], TmpMax[Z]);
}
}
Assign_Vector(ClipMin, Min);
Assign_Vector(ClipMax, Max);
}
/* Translate clipping box. */
VSubEq(ClipMin, T1);
VSubEq(ClipMax, T1);
/* We want A to be non-negative. */
if (A < 0.0)
{
A = -A;
D = -D;
E = -E;
G = -G;
H = -H;
I = -I;
J = -J;
}
/*
*
* Check for ellipsoid.
*
* x*x y*y z*z
* ----- + ----- + ----- - 1 = 0
* a*a b*b c*c
*
*/
if ((A > 0.0) && (E > 0.0) && (H > 0.0) && (J < 0.0))
{
a = sqrt(-J/A);
b = sqrt(-J/E);
c = sqrt(-J/H);
NewMin[X] = -a;
NewMin[Y] = -b;
NewMin[Z] = -c;
NewMax[X] = a;
NewMax[Y] = b;
NewMax[Z] = c;
}
/*
*
* Check for cylinder (x-axis).
*
* y*y z*z
* ----- + ----- - 1 = 0
* b*b c*c
*
*/
if ((A == 0.0) && (E > 0.0) && (H > 0.0) && (J < 0.0))
{
b = sqrt(-J/E);
c = sqrt(-J/H);
NewMin[Y] = -b;
NewMin[Z] = -c;
NewMax[Y] = b;
NewMax[Z] = c;
}
/*
*
* Check for cylinder (y-axis).
*
* x*x z*z
* ----- + ----- - 1 = 0
* a*a c*c
*
*/
if ((A > 0.0) && (E == 0.0) && (H > 0.0) && (J < 0.0))
{
a = sqrt(-J/A);
c = sqrt(-J/H);
NewMin[X] = -a;
NewMin[Z] = -c;
NewMax[X] = a;
NewMax[Z] = c;
}
/*
*
* Check for cylinder (z-axis).
*
* x*x y*y
* ----- + ----- - 1 = 0
* a*a b*b
*
*/
if ((A > 0.0) && (E > 0.0) && (H == 0.0) && (J < 0.0))
{
a = sqrt(-J/A);
b = sqrt(-J/E);
NewMin[X] = -a;
NewMin[Y] = -b;
NewMax[X] = a;
NewMax[Y] = b;
}
/*
*
* Check for cone (x-axis).
*
* x*x y*y z*z
* ----- - ----- - ----- = 0
* a*a b*b c*c
*
*/
if ((A > 0.0) && (E < 0.0) && (H < 0.0) && (J == 0.0))
{
a = sqrt(1.0/A);
b = sqrt(-1.0/E);
c = sqrt(-1.0/H);
/* Get radii for lower x value. */
x = ClipMin[X];
ry1 = fabs(x * b / a);
rz1 = fabs(x * c / a);
/* Get radii for upper x value. */
x = ClipMax[X];
ry2 = fabs(x * b / a);
rz2 = fabs(x * c / a);
ry = max(ry1, ry2);
rz = max(rz1, rz2);
NewMin[Y] = -ry;
NewMin[Z] = -rz;
NewMax[Y] = ry;
NewMax[Z] = rz;
}
/*
*
* Check for cone (y-axis).
*
* x*x y*y z*z
* ----- - ----- + ----- = 0
* a*a b*b c*c
*
*/
if ((A > 0.0) && (E < 0.0) && (H > 0.0) && (J == 0.0))
{
a = sqrt(1.0/A);
b = sqrt(-1.0/E);
c = sqrt(1.0/H);
/* Get radii for lower y value. */
y = ClipMin[Y];
rx1 = fabs(y * a / b);
rz1 = fabs(y * c / b);
/* Get radii for upper y value. */
y = ClipMax[Y];
rx2 = fabs(y * a / b);
rz2 = fabs(y * c / b);
rx = max(rx1, rx2);
rz = max(rz1, rz2);
NewMin[X] = -rx;
NewMin[Z] = -rz;
NewMax[X] = rx;
NewMax[Z] = rz;
}
/*
*
* Check for cone (z-axis).
*
* x*x y*y z*z
* ----- + ----- - ----- = 0
* a*a b*b c*c
*
*/
if ((A > 0.0) && (E > 0.0) && (H < 0.0) && (J == 0.0))
{
a = sqrt(1.0/A);
b = sqrt(1.0/E);
c = sqrt(-1.0/H);
/* Get radii for lower z value. */
z = ClipMin[Z];
rx1 = fabs(z * a / c);
ry1 = fabs(z * b / c);
/* Get radii for upper z value. */
z = ClipMax[Z];
rx2 = fabs(z * a / c);
ry2 = fabs(z * b / c);
rx = max(rx1, rx2);
ry = max(ry1, ry2);
NewMin[X] = -rx;
NewMin[Y] = -ry;
NewMax[X] = rx;
NewMax[Y] = ry;
}
/*
*
* Check for hyperboloid (x-axis).
*
* x*x y*y z*z
* ----- - ----- - ----- + 1 = 0
* a*a b*b c*c
*
*/
if ((A > 0.0) && (E < 0.0) && (H < 0.0) && (J > 0.0))
{
/* Get radii for lower x value. */
x = ClipMin[X];
d = 1.0 + A * Sqr(x);
ry1 = sqrt(-d / E);
rz1 = sqrt(-d / H);
/* Get radii for upper x value. */
x = ClipMax[X];
d = 1.0 + A * Sqr(x);
ry2 = sqrt(-d / E);
rz2 = sqrt(-d / H);
ry = max(ry1, ry2);
rz = max(rz1, rz2);
NewMin[Y] = -ry;
NewMin[Z] = -rz;
NewMax[Y] = ry;
NewMax[Z] = rz;
}
/*
*
* Check for hyperboloid (y-axis).
*
* x*x y*y z*z
* ----- - ----- + ----- - 1 = 0
* a*a b*b c*c
*
*/
if ((A > 0.0) && (E < 0.0) && (H > 0.0) && (J < 0.0))
{
/* Get radii for lower y value. */
y = ClipMin[Y];
d = 1.0 - E * Sqr(y);
rx1 = sqrt(d / A);
rz1 = sqrt(d / H);
/* Get radii for upper y value. */
y = ClipMax[Y];
d = 1.0 - E * Sqr(y);
rx2 = sqrt(d / A);
rz2 = sqrt(d / H);
rx = max(rx1, rx2);
rz = max(rz1, rz2);
NewMin[X] = -rx;
NewMin[Z] = -rz;
NewMax[X] = rx;
NewMax[Z] = rz;
}
/*
*
* Check for hyperboloid (z-axis).
*
* x*x y*y z*z
* ----- + ----- - ----- - 1 = 0
* a*a b*b c*c
*
*/
if ((A > 0.0) && (E > 0.0) && (H < 0.0) && (J < 0.0))
{
/* Get radii for lower z value. */
z = ClipMin[Z];
d = 1.0 - H * Sqr(z);
rx1 = sqrt(d / A);
ry1 = sqrt(d / E);
/* Get radii for upper z value. */
z = ClipMax[Z];
d = 1.0 - H * Sqr(z);
rx2 = sqrt(d / A);
ry2 = sqrt(d / E);
rx = max(rx1, rx2);
ry = max(ry1, ry2);
NewMin[X] = -rx;
NewMin[Y] = -ry;
NewMax[X] = rx;
NewMax[Y] = ry;
}
/*
*
* Check for paraboloid (x-axis).
*
* y*y z*z
* x - ----- - ----- = 0
* b*b c*c
*
*/
if ((A == 0.0) && (D != 0.0) && (E != 0.0) && (H != 0.0) && (J == 0.0))
{
/* Get radii for lower x value. */
x = ClipMin[X];
ry1 = sqrt(fabs(2.0 * D * x / E));
rz1 = sqrt(fabs(2.0 * D * x / H));
/* Get radii for upper x value. */
x = ClipMax[X];
ry2 = sqrt(fabs(2.0 * D * x / E));
rz2 = sqrt(fabs(2.0 * D * x / H));
ry = max(ry1, ry2);
rz = max(rz1, rz2);
NewMin[Y] = -ry;
NewMin[Z] = -rz;
NewMax[Y] = ry;
NewMax[Z] = rz;
}
/*
*
* Check for paraboloid (y-axis).
*
* x*x z*z
* y - ----- - ----- = 0
* a*a c*c
*
*/
if ((E == 0.0) && (G != 0.0) && (A != 0.0) && (H != 0.0) && (J == 0.0))
{
/* Get radii for lower y-value. */
y = ClipMin[Y];
rx1 = sqrt(fabs(2.0 * G * y / A));
rz1 = sqrt(fabs(2.0 * G * y / H));
/* Get radii for upper y value. */
y = ClipMax[Y];
rx2 = sqrt(fabs(2.0 * G * y / A));
rz2 = sqrt(fabs(2.0 * G * y / H));
rx = max(rx1, rx2);
rz = max(rz1, rz2);
NewMin[X] = -rx;
NewMin[Z] = -rz;
NewMax[X] = rx;
NewMax[Z] = rz;
}
/*
*
* Check for paraboloid (z-axis).
*
* x*x y*y
* z - ----- - ----- = 0
* a*a b*b
*
*/
if ((H == 0.0) && (I != 0.0) && (A != 0.0) && (E != 0.0) && (J == 0.0))
{
/* Get radii for lower z-value. */
z = ClipMin[Z];
rx1 = sqrt(fabs(2.0 * I * z / A));
ry1 = sqrt(fabs(2.0 * I * z / E));
/* Get radii for upper z value. */
z = ClipMax[Z];
rx2 = sqrt(fabs(2.0 * I * z / A));
ry2 = sqrt(fabs(2.0 * I * z / E));
rx = max(rx1, rx2);
ry = max(ry1, ry2);
NewMin[X] = -rx;
NewMin[Y] = -ry;
NewMax[X] = rx;
NewMax[Y] = ry;
}
/* Intersect clipping object's and quadric's bounding boxes */
NewMin[X] = max(NewMin[X], ClipMin[X]);
NewMin[Y] = max(NewMin[Y], ClipMin[Y]);
NewMin[Z] = max(NewMin[Z], ClipMin[Z]);
NewMax[X] = min(NewMax[X], ClipMax[X]);
NewMax[Y] = min(NewMax[Y], ClipMax[Y]);
NewMax[Z] = min(NewMax[Z], ClipMax[Z]);
/* Use old or new bounding box? */
New_Volume = (NewMax[X] - NewMin[X]) * (NewMax[Y] - NewMin[Y]) * (NewMax[Z] - NewMin[Z]);
BOUNDS_VOLUME(Old_Volume, Old);
if (New_Volume < Old_Volume)
{
/* Add translation. */
VAddEq(NewMin, T1);
VAddEq(NewMax, T1);
Make_BBox_from_min_max(Quadric->BBox, NewMin, NewMax);
/* Beware of bounding boxes to large. */
if ((Quadric->BBox.Lengths[X] > CRITICAL_LENGTH) ||
(Quadric->BBox.Lengths[Y] > CRITICAL_LENGTH) ||
(Quadric->BBox.Lengths[Z] > CRITICAL_LENGTH))
{
Make_BBox(Quadric->BBox, -BOUND_HUGE/2, -BOUND_HUGE/2, -BOUND_HUGE/2,
BOUND_HUGE, BOUND_HUGE, BOUND_HUGE);
}
}
}
/*****************************************************************************
*
* FUNCTION
*
* Compute_Plane_Min_Max
*
* INPUT
*
* Plane - Plane
* Min, Max - Vectors containing plane's dimensions
*
* OUTPUT
*
* Min, Max
*
* RETURNS
*
* AUTHOR
*
* Dieter Bayer
*
* DESCRIPTION
*
* Compute min/max vectors for planes perpendicular to an axis.
*
* CHANGES
*
* May 1994 : Creation.
*
******************************************************************************/
void Compute_Plane_Min_Max(Plane, Min, Max)
PLANE *Plane;
VECTOR Min, Max;
{
DBL d;
VECTOR N;
Assign_Vector(N, Plane->Normal_Vector);
d = -Plane->Distance;
Min[X] = Min[Y] = Min[Z] = -BOUND_HUGE/2;
Max[X] = Max[Y] = Max[Z] = BOUND_HUGE/2;
/* y-z-plane */
if (fabs(1.0 - fabs(N[X])) < EPSILON)
{
if (N[X] > 0.0)
{
Max[X] = d;
}
else
{
Min[X] = -d;
}
}
/* x-z-plane */
if (fabs(1.0 - fabs(N[Y])) < EPSILON)
{
if (N[Y] > 0.0)
{
Max[Y] = d;
}
else
{
Min[Y] = -d;
}
}
/* x-y-plane */
if (fabs(1.0 - fabs(N[Z])) < EPSILON)
{
if (N[Z] > 0.0)
{
Max[Z] = d;
}
else
{
Min[Z] = -d;
}
}
}
