/**************************************************************************** * poly.c * * This module implements the code for general 3 variable polynomial shapes * * This file was written by Alexander Enzmann. He wrote the code for * 4th - 6th order shapes and generously provided us these enhancements. * * 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 "vector.h" #include "povproto.h" #include "bbox.h" #include "polysolv.h" #include "matrices.h" #include "objects.h" #include "poly.h" #include "povray.h" /* * Basic form of a quartic equation: * * a00*x^4 + a01*x^3*y + a02*x^3*z + a03*x^3 + a04*x^2*y^2 + * a05*x^2*y*z+ a06*x^2*y + a07*x^2*z^2 + a08*x^2*z + a09*x^2 + * a10*x*y^3 + a11*x*y^2*z + a12*x*y^2 + a13*x*y*z^2 + a14*x*y*z + * a15*x*y + a16*x*z^3 + a17*x*z^2 + a18*x*z + a19*x + a20*y^4 + * a21*y^3*z + a22*y^3 + a23*y^2*z^2 + a24*y^2*z + a25*y^2 + a26*y*z^3 + * a27*y*z^2 + a28*y*z + a29*y + a30*z^4 + a31*z^3 + a32*z^2 + a33*z + a34 * */ /***************************************************************************** * Local preprocessor defines ******************************************************************************/ #define DEPTH_TOLERANCE 1.0e-4 #define INSIDE_TOLERANCE 1.0e-4 #define ROOT_TOLERANCE 1.0e-4 #define COEFF_LIMIT 1.0e-20 #define BINOMSIZE 40 /***************************************************************************** * Local typedefs ******************************************************************************/ /***************************************************************************** * Static functions ******************************************************************************/ static int intersect PARAMS((RAY *Ray, int Order, DBL *Coeffs, int Sturm_Flag, DBL *Depths)); static void normal0 PARAMS((VECTOR Result, int Order, DBL *Coeffs, VECTOR IPoint)); static void normal1 PARAMS((VECTOR Result, int Order, DBL *Coeffs, VECTOR IPoint)); static DBL inside PARAMS((VECTOR IPoint, int Order, DBL *Coeffs)); static int intersect_linear PARAMS((RAY *ray, DBL *Coeffs, DBL *Depths)); static int intersect_quadratic PARAMS((RAY *ray, DBL *Coeffs, DBL *Depths)); static int factor_out PARAMS((int n, int i, int *c, int *s)); static long binomial PARAMS((int n, int r)); static void factor1 PARAMS((int n, int *c, int *s)); /* unused static DBL evaluate_linear PARAMS((VECTOR P, DBL *a)); static DBL evaluate_quadratic PARAMS((VECTOR P, DBL *a)); */ static int All_Poly_Intersections PARAMS((OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack)); static int Inside_Poly PARAMS((VECTOR IPoint, OBJECT *Object)); static void Poly_Normal PARAMS((VECTOR Result, OBJECT *Object, INTERSECTION *Inter)); static void *Copy_Poly PARAMS((OBJECT *Object)); static void Translate_Poly PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)); static void Rotate_Poly PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)); static void Scale_Poly PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)); static void Transform_Poly PARAMS((OBJECT *Object, TRANSFORM *Trans)); static void Invert_Poly PARAMS((OBJECT *Object)); static void Destroy_Poly PARAMS((OBJECT *Object)); /***************************************************************************** * Local variables ******************************************************************************/ METHODS Poly_Methods = { All_Poly_Intersections, Inside_Poly, Poly_Normal, Copy_Poly, Translate_Poly, Rotate_Poly, Scale_Poly, Transform_Poly, Invert_Poly, Destroy_Poly }; /* The following table contains the binomial coefficients up to 15 */ static int binomials[15][15] = { {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 3, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 4, 6, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 5, 10, 10, 5, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 6, 15, 20, 15, 6, 1, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 7, 21, 35, 35, 21, 7, 1, 0, 0, 0, 0, 0, 0, 0}, {1, 8, 28, 56, 70, 56, 28, 8, 1, 0, 0, 0, 0, 0, 0}, {1, 9, 36, 84, 126, 126, 84, 36, 9, 1, 0, 0, 0, 0, 0}, {1, 10, 45,120, 210, 252, 210, 120, 45, 10, 1, 0, 0, 0, 0}, {1, 11, 55,165, 330, 462, 462, 330, 165, 55, 11, 1, 0, 0, 0}, {1, 12, 66,220, 495, 792, 924, 792, 495, 220, 66, 12, 1, 0, 0}, {1, 13, 78,286, 715,1287,1716,1716,1287, 715, 286, 78, 13, 1, 0}, {1, 14, 91,364,1001,2002,3003,3432,3003,2002,1001,364, 91, 14, 1} }; static DBL eqn_v[3][MAX_ORDER+1], eqn_vt[3][MAX_ORDER+1]; /***************************************************************************** * * FUNCTION * * All_Poly_Intersections * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * - * * CHANGES * * - * ******************************************************************************/ static int All_Poly_Intersections(Object, Ray, Depth_Stack) OBJECT *Object; RAY *Ray; ISTACK *Depth_Stack; { POLY *Poly = (POLY *) Object; DBL Depths[MAX_ORDER], len; VECTOR IPoint; int cnt, i, j, Intersection_Found, same_root; RAY New_Ray; /* Transform the ray into the polynomial's space */ MInvTransPoint(New_Ray.Initial, Ray->Initial, Poly->Trans); MInvTransDirection(New_Ray.Direction, Ray->Direction, Poly->Trans); VLength(len, New_Ray.Direction); VInverseScaleEq(New_Ray.Direction, len); Intersection_Found = FALSE; Increase_Counter(stats[Ray_Poly_Tests]); switch (Poly->Order) { case 1: cnt = intersect_linear(&New_Ray, Poly->Coeffs, Depths); break; case 2: cnt = intersect_quadratic(&New_Ray, Poly->Coeffs, Depths); break; default: cnt = intersect(&New_Ray, Poly->Order, Poly->Coeffs, Test_Flag(Poly, STURM_FLAG), Depths); } if (cnt > 0) { Increase_Counter(stats[Ray_Poly_Tests_Succeeded]); } for (i = 0; i < cnt; i++) { if (Depths[i] > DEPTH_TOLERANCE) { same_root = FALSE; for (j = 0; j < i; j++) { if (Depths[i] == Depths[j]) { same_root = TRUE; break; } } if (!same_root) { VEvaluateRay(IPoint, New_Ray.Initial, Depths[i], New_Ray.Direction); /* Transform the point into world space */ MTransPoint(IPoint, IPoint, Poly->Trans); if (Point_In_Clip(IPoint, Object->Clip)) { push_entry(Depths[i] / len,IPoint,Object,Depth_Stack); Intersection_Found = TRUE; } } } } return (Intersection_Found); } /***************************************************************************** * * FUNCTION * * evaluate_linear * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * - * * CHANGES * * - * ******************************************************************************/ /* For speedup of low order polynomials, expand out the terms involved in evaluating the poly. */ /* unused static DBL evaluate_linear(P, a) VECTOR P; DBL *a; { return(a[0] * P[X]) + (a[1] * P[Y]) + (a[2] * P[Z]) + a[3]; } */ /***************************************************************************** * * FUNCTION * * evaluate_quadratic * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * - * * CHANGES * * - * ******************************************************************************/ /* static DBL evaluate_quadratic(P, a) VECTOR P; DBL *a; { DBL x, y, z; x = P[X]; y = P[Y]; z = P[Z]; return(a[0] * x * x + a[1] * x * y + a[2] * x * z + a[3] * x + a[4] * y * y + a[5] * y * z + a[6] * y + a[7] * z * z + a[8] * z + a[9]); } */ /***************************************************************************** * * FUNCTION * * factor_out * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Remove all factors of i from n. * * CHANGES * * - * ******************************************************************************/ static int factor_out(n, i, c, s) int n, i, *c, *s; { while (!(n % i)) { n /= i; s[(*c)++] = i; } return(n); } /***************************************************************************** * * FUNCTION * * factor1 * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Find all prime factors of n. (Note that n must be less than 2^15. * * CHANGES * * - * ******************************************************************************/ static void factor1(n, c, s) int n, *c, *s; { int i,k; /* First factor out any 2s. */ n = factor_out(n, 2, c, s); /* Now any odd factors. */ k = (int)sqrt((DBL)n) + 1; for (i = 3; (n > 1) && (i <= k); i += 2) { if (!(n % i)) { n = factor_out(n, i, c, s); k = (int)sqrt((DBL)n)+1; } } if (n > 1) { s[(*c)++] = n; } } /***************************************************************************** * * FUNCTION * * binomial * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Calculate the binomial coefficent of n, r. * * CHANGES * * - * ******************************************************************************/ static long binomial(n, r) int n, r; { int h,i,j,k,l; unsigned long result; static int stack1[BINOMSIZE], stack2[BINOMSIZE]; if ((n < 0) || (r < 0) || (r > n)) { result = 0L; } else { if (r == n) { result = 1L; } else { if ((r < 15) && (n < 15)) { result = (long)binomials[n][r]; } else { j = 0; for (i = r + 1; i <= n; i++) { stack1[j++] = i; } for (i = 2; i <= (n-r); i++) { h = 0; factor1(i, &h, stack2); for (k = 0; k < h; k++) { for (l = 0; l < j; l++) { if (!(stack1[l] % stack2[k])) { stack1[l] /= stack2[k]; goto l1; } } } /* Error if we get here */ /* Debug_Info("Failed to factor\n");*/ l1:; } result = 1; for (i = 0; i < j; i++) { result *= stack1[i]; } } } } return(result); } /***************************************************************************** * * FUNCTION * * inside * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * - * * CHANGES * * - * ******************************************************************************/ static DBL inside(IPoint, Order, Coeffs) VECTOR IPoint; int Order; DBL *Coeffs; { DBL x[MAX_ORDER+1], y[MAX_ORDER+1]; DBL z[MAX_ORDER+1], c, Result; int i, j, k, term; x[0] = 1.0; y[0] = 1.0; z[0] = 1.0; x[1] = IPoint[X]; y[1] = IPoint[Y]; z[1] = IPoint[Z]; for (i = 2; i <= Order; i++) { x[i] = x[1] * x[i-1]; y[i] = y[1] * y[i-1]; z[i] = z[1] * z[i-1]; } Result = 0.0; term = 0; for (i = Order; i >= 0; i--) { for (j=Order-i;j>=0;j--) { for (k=Order-(i+j);k>=0;k--) { if ((c = Coeffs[term]) != 0.0) { Result += c * x[i] * y[j] * z[k]; } term++; } } } return(Result); } /***************************************************************************** * * FUNCTION * * intersect * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Intersection of a ray and an arbitrary polynomial function. * * CHANGES * * - * ******************************************************************************/ static int intersect(ray, Order, Coeffs, Sturm_Flag, Depths) RAY *ray; int Order, Sturm_Flag; DBL *Coeffs, *Depths; { DBL eqn[MAX_ORDER+1]; DBL t[3][MAX_ORDER+1]; VECTOR P, D; DBL val; int h, i, j, k, i1, j1, k1, term; int offset; /* First we calculate the values of the individual powers of x, y, and z as they are represented by the ray */ Assign_Vector(P,ray->Initial); Assign_Vector(D,ray->Direction); for (i = 0; i < 3; i++) { eqn_v[i][0] = 1.0; eqn_vt[i][0] = 1.0; } eqn_v[0][1] = P[X]; eqn_v[1][1] = P[Y]; eqn_v[2][1] = P[Z]; eqn_vt[0][1] = D[X]; eqn_vt[1][1] = D[Y]; eqn_vt[2][1] = D[Z]; for (i = 2; i <= Order; i++) { for (j=0;j<3;j++) { eqn_v[j][i] = eqn_v[j][1] * eqn_v[j][i-1]; eqn_vt[j][i] = eqn_vt[j][1] * eqn_vt[j][i-1]; } } for (i = 0; i <= Order; i++) { eqn[i] = 0.0; } /* Now walk through the terms of the polynomial. As we go we substitute the ray equation for each of the variables. */ term = 0; for (i = Order; i >= 0; i--) { for (h = 0; h <= i; h++) { t[0][h] = binomial(i, h) * eqn_vt[0][i-h] * eqn_v[0][h]; } for (j = Order-i; j >= 0; j--) { for (h = 0; h <= j; h++) { t[1][h] = binomial(j, h) * eqn_vt[1][j-h] * eqn_v[1][h]; } for (k = Order-(i+j); k >= 0; k--) { if (Coeffs[term] != 0) { for (h = 0; h <= k; h++) { t[2][h] = binomial(k, h) * eqn_vt[2][k-h] * eqn_v[2][h]; } /* Multiply together the three polynomials. */ offset = Order - (i + j + k); for (i1 = 0; i1 <= i; i1++) { for (j1=0;j1<=j;j1++) { for (k1=0;k1<=k;k1++) { h = offset + i1 + j1 + k1; val = Coeffs[term]; val *= t[0][i1]; val *= t[1][j1]; val *= t[2][k1]; eqn[h] += val; } } } } term++; } } } for (i = 0, j = Order; i <= Order; i++) { if (eqn[i] != 0.0) { break; } else { j--; } } if (j > 1) { return(Solve_Polynomial(j, &eqn[i], Depths, Sturm_Flag, ROOT_TOLERANCE)); } else { return(0); } } /***************************************************************************** * * FUNCTION * * intersect_linear * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Intersection of a ray and a quadratic. * * CHANGES * * - * ******************************************************************************/ static int intersect_linear(ray, Coeffs, Depths) RAY *ray; DBL *Coeffs, *Depths; { DBL t0, t1, *a = Coeffs; t0 = a[0] * ray->Initial[X] + a[1] * ray->Initial[Y] + a[2] * ray->Initial[Z]; t1 = a[0] * ray->Direction[X] + a[1] * ray->Direction[Y] + a[2] * ray->Direction[Z]; if (fabs(t1) < EPSILON) { return(0); } Depths[0] = -(a[3] + t0) / t1; return(1); } /***************************************************************************** * * FUNCTION * * intersect_quadratic * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Intersection of a ray and a quadratic. * * CHANGES * * - * ******************************************************************************/ static int intersect_quadratic(ray, Coeffs, Depths) RAY *ray; DBL *Coeffs, *Depths; { DBL x, y, z, x2, y2, z2; DBL xx, yy, zz, xx2, yy2, zz2; DBL *a, ac, bc, cc, d, t; x = ray->Initial[X]; y = ray->Initial[Y]; z = ray->Initial[Z]; xx = ray->Direction[X]; yy = ray->Direction[Y]; zz = ray->Direction[Z]; x2 = x * x; y2 = y * y; z2 = z * z; xx2 = xx * xx; yy2 = yy * yy; zz2 = zz * zz; a = Coeffs; /* * Determine the coefficients of t^n, where the line is represented * as (x,y,z) + (xx,yy,zz)*t. */ ac = (a[0]*xx2 + a[1]*xx*yy + a[2]*xx*zz + a[4]*yy2 + a[5]*yy*zz + a[7]*zz2); bc = (2*a[0]*x*xx + a[1]*(x*yy + xx*y) + a[2]*(x*zz + xx*z) + a[3]*xx + 2*a[4]*y*yy + a[5]*(y*zz + yy*z) + a[6]*yy + 2*a[7]*z*zz + a[8]*zz); cc = a[0]*x2 + a[1]*x*y + a[2]*x*z + a[3]*x + a[4]*y2 + a[5]*y*z + a[6]*y + a[7]*z2 + a[8]*z + a[9]; if (fabs(ac) < COEFF_LIMIT) { if (fabs(bc) < COEFF_LIMIT) { return(0); } Depths[0] = -cc / bc; return(1); } /* * Solve the quadratic formula & return results that are * within the correct interval. */ d = bc * bc - 4.0 * ac * cc; if (d < 0.0) { return(0); } d = sqrt(d); bc = -bc; t = 2.0 * ac; Depths[0] = (bc + d) / t; Depths[1] = (bc - d) / t; return(2); } /***************************************************************************** * * FUNCTION * * normal0 * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Normal to a polynomial (used for polynomials with order > 4). * * CHANGES * * - * ******************************************************************************/ static void normal0(Result, Order, Coeffs, IPoint) VECTOR Result; int Order; DBL *Coeffs; VECTOR IPoint; { int i, j, k, term; DBL val, *a, x[MAX_ORDER+1], y[MAX_ORDER+1], z[MAX_ORDER+1]; x[0] = 1.0; y[0] = 1.0; z[0] = 1.0; x[1] = IPoint[X]; y[1] = IPoint[Y]; z[1] = IPoint[Z]; for (i = 2; i <= Order; i++) { x[i] = IPoint[X] * x[i-1]; y[i] = IPoint[Y] * y[i-1]; z[i] = IPoint[Z] * z[i-1]; } a = Coeffs; term = 0; Make_Vector(Result, 0.0, 0.0, 0.0); for (i = Order; i >= 0; i--) { for (j = Order-i; j >= 0; j--) { for (k = Order-(i+j); k >= 0; k--) { if (i >= 1) { val = x[i-1] * y[j] * z[k]; Result[X] += i * a[term] * val; } if (j >= 1) { val = x[i] * y[j-1] * z[k]; Result[Y] += j * a[term] * val; } if (k >= 1) { val = x[i] * y[j] * z[k-1]; Result[Z] += k * a[term] * val; } term++; } } } } /***************************************************************************** * * FUNCTION * * nromal1 * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Normal to a polynomial (for polynomials of order <= 4). * * CHANGES * * - * ******************************************************************************/ static void normal1(Result, Order, Coeffs, IPoint) VECTOR Result; int Order; DBL *Coeffs; VECTOR IPoint; { DBL *a, x, y, z, x2, y2, z2, x3, y3, z3; a = Coeffs; x = IPoint[X]; y = IPoint[Y]; z = IPoint[Z]; switch (Order) { case 1: /* Linear partial derivatives */ Make_Vector(Result, a[0], a[1], a[2]) break; case 2: /* Quadratic partial derivatives */ Result[X] = 2*a[0]*x+a[1]*y+a[2]*z+a[3]; Result[Y] = a[1]*x+2*a[4]*y+a[5]*z+a[6]; Result[Z] = a[2]*x+a[5]*y+2*a[7]*z+a[8]; break; case 3: x2 = x * x; y2 = y * y; z2 = z * z; /* Cubic partial derivatives */ Result[X] = 3*a[0]*x2 + 2*x*(a[1]*y + a[2]*z + a[3]) + a[4]*y2 + y*(a[5]*z + a[6]) + a[7]*z2 + a[8]*z + a[9]; Result[Y] = a[1]*x2 + x*(2*a[4]*y + a[5]*z + a[6]) + 3*a[10]*y2 + 2*y*(a[11]*z + a[12]) + a[13]*z2 + a[14]*z + a[15]; Result[Z] = a[2]*x2 + x*(a[5]*y + 2*a[7]*z + a[8]) + a[11]*y2 + y*(2*a[13]*z + a[14]) + 3*a[16]*z2 + 2*a[17]*z + a[18]; break; case 4: /* Quartic partial derivatives */ x2 = x * x; y2 = y * y; z2 = z * z; x3 = x * x2; y3 = y * y2; z3 = z * z2; Result[X] = 4*a[ 0]*x3+3*x2*(a[ 1]*y+a[ 2]*z+a[ 3])+ 2*x*(a[ 4]*y2+y*(a[ 5]*z+a[ 6])+a[ 7]*z2+a[ 8]*z+a[ 9])+ a[10]*y3+y2*(a[11]*z+a[12])+y*(a[13]*z2+a[14]*z+a[15])+ a[16]*z3+a[17]*z2+a[18]*z+a[19]; Result[Y] = a[ 1]*x3+x2*(2*a[ 4]*y+a[ 5]*z+a[ 6])+ x*(3*a[10]*y2+2*y*(a[11]*z+a[12])+a[13]*z2+a[14]*z+a[15])+ 4*a[20]*y3+3*y2*(a[21]*z+a[22])+2*y*(a[23]*z2+a[24]*z+a[25])+ a[26]*z3+a[27]*z2+a[28]*z+a[29]; Result[Z] = a[ 2]*x3+x2*(a[ 5]*y+2*a[ 7]*z+a[ 8])+ x*(a[11]*y2+y*(2*a[13]*z+a[14])+3*a[16]*z2+2*a[17]*z+a[18])+ a[21]*y3+y2*(2*a[23]*z+a[24])+y*(3*a[26]*z2+2*a[27]*z+a[28])+ 4*a[30]*z3+3*a[31]*z2+2*a[32]*z+a[33]; } } /***************************************************************************** * * FUNCTION * * Inside_Poly * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * - * * CHANGES * * - * ******************************************************************************/ static int Inside_Poly (IPoint, Object) VECTOR IPoint; OBJECT *Object; { VECTOR New_Point; DBL Result; /* Transform the point into polynomial's space */ MInvTransPoint(New_Point, IPoint, ((POLY *)Object)->Trans); Result = inside(New_Point, ((POLY *)Object)->Order, ((POLY *)Object)->Coeffs); if (Result < INSIDE_TOLERANCE) { return ((int)(!Test_Flag(Object, INVERTED_FLAG))); } else { return ((int)Test_Flag(Object, INVERTED_FLAG)); } } /***************************************************************************** * * FUNCTION * * Poly_Normal * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Normal to a polynomial. * * CHANGES * * - * ******************************************************************************/ static void Poly_Normal(Result, Object, Inter) OBJECT *Object; VECTOR Result; INTERSECTION *Inter; { DBL val; VECTOR New_Point; POLY *Shape = (POLY *)Object; /* Transform the point into the polynomials space. */ MInvTransPoint(New_Point, Inter->IPoint, Shape->Trans); if (((POLY *)Object)->Order > 4) { normal0(Result, Shape->Order, Shape->Coeffs, New_Point); } else { normal1(Result, Shape->Order, Shape->Coeffs, New_Point); } /* Transform back to world space. */ MTransNormal(Result, Result, Shape->Trans); /* Normalize (accounting for the possibility of a 0 length normal). */ VDot(val, Result, Result); if (val > 0.0) { val = 1.0 / sqrt(val); VScaleEq(Result, val); } else { Make_Vector(Result, 1.0, 0.0, 0.0) } } /***************************************************************************** * * FUNCTION * * Translate_Poly * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * - * * CHANGES * * - * ******************************************************************************/ static void Translate_Poly (Object, Vector, Trans) OBJECT *Object; VECTOR Vector; TRANSFORM *Trans; { Transform_Poly(Object, Trans); } /***************************************************************************** * * FUNCTION * * Rotate_Poly * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * - * * CHANGES * * - * ******************************************************************************/ static void Rotate_Poly (Object, Vector, Trans) OBJECT *Object; VECTOR Vector; TRANSFORM *Trans; { Transform_Poly(Object, Trans); } /***************************************************************************** * * FUNCTION * * Scale_Poly * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * - * * CHANGES * * - * ******************************************************************************/ static void Scale_Poly (Object, Vector, Trans) OBJECT *Object; VECTOR Vector; TRANSFORM *Trans; { Transform_Poly(Object, Trans); } /***************************************************************************** * * FUNCTION * * Transform_Poly * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * - * * CHANGES * * - * ******************************************************************************/ static void Transform_Poly(Object,Trans) OBJECT *Object; TRANSFORM *Trans; { Compose_Transforms(((POLY *)Object)->Trans, Trans); Compute_Poly_BBox((POLY *)Object); } /***************************************************************************** * * FUNCTION * * Invert_Poly * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * - * * CHANGES * * - * ******************************************************************************/ static void Invert_Poly(Object) OBJECT *Object; { Invert_Flag(Object, INVERTED_FLAG); } /***************************************************************************** * * FUNCTION * * Create_Poly * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * - * * CHANGES * * - * ******************************************************************************/ POLY *Create_Poly(Order) int Order; { POLY *New; int i; New = (POLY *)POV_MALLOC(sizeof (POLY), "poly"); INIT_OBJECT_FIELDS(New,POLY_OBJECT, &Poly_Methods); New->Order = Order; New->Trans = Create_Transform(); New->Coeffs = (DBL *)POV_MALLOC(term_counts(Order) * sizeof(DBL), "coefficients for POLY"); for (i = 0; i < term_counts(Order); i++) { New->Coeffs[i] = 0.0; } return(New); } /***************************************************************************** * * FUNCTION * * Copy_Poly * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * Make a copy of a polynomial object. * * CHANGES * * - * ******************************************************************************/ static void *Copy_Poly(Object) OBJECT *Object; { POLY *New; int i; New = Create_Poly (((POLY *)Object)->Order); /* Get rid of transform created in Create_Poly. */ Destroy_Transform(New->Trans); Copy_Flag(New, Object, STURM_FLAG); Copy_Flag(New, Object, INVERTED_FLAG); New->Trans = Copy_Transform(((POLY *)Object)->Trans); for (i = 0; i < term_counts(New->Order); i++) { New->Coeffs[i] = ((POLY *)Object)->Coeffs[i]; } return (New); } /***************************************************************************** * * FUNCTION * * Destroy_Poly * * INPUT * * OUTPUT * * RETURNS * * AUTHOR * * Alexander Enzmann * * DESCRIPTION * * - * * CHANGES * * - * ******************************************************************************/ static void Destroy_Poly(Object) OBJECT *Object; { Destroy_Transform (((POLY *)Object)->Trans); POV_FREE (((POLY *)Object)->Coeffs); POV_FREE (Object); } /***************************************************************************** * * FUNCTION * * Compute_Poly_BBox * * INPUT * * Poly - Poly * * OUTPUT * * Poly * * RETURNS * * AUTHOR * * Dieter Bayer * * DESCRIPTION * * Calculate the bounding box of a poly. * * CHANGES * * Aug 1994 : Creation. * ******************************************************************************/ void Compute_Poly_BBox(Poly) POLY *Poly; { Make_BBox(Poly->BBox, -BOUND_HUGE/2, -BOUND_HUGE/2, -BOUND_HUGE/2, BOUND_HUGE, BOUND_HUGE, BOUND_HUGE); if (Poly->Clip != NULL) { Poly->BBox = Poly->Clip->BBox; } }