home *** CD-ROM | disk | FTP | other *** search
/ rtsi.com / 2014.01.www.rtsi.com.tar / www.rtsi.com / OS9 / OSK / GRAPHICS / rayshade.lzh / superq.c < prev    next >
Text File  |  1990-05-08  |  3KB  |  160 lines
  1. /*
  2.  * superq.c
  3.  *
  4.  * Slightly modified version of Kuchkuda's superquadric code.
  5.  *
  6.  * $Id: superq.c,v 3.0 89/10/27 02:06:05 craig Exp $
  7.  *
  8.  * $Log:    superq.c,v $
  9.  * Revision 3.0  89/10/27  02:06:05  craig
  10.  * Baseline for first official release.
  11.  * 
  12.  */
  13. #include <math.h>
  14. #include "constants.h"
  15. #include "typedefs.h"
  16. #include "funcdefs.h"
  17.  
  18. #define ERRCONST 1.006
  19.  
  20. /*
  21.  * Create and return reference to a superquadric.
  22.  */
  23. Object *
  24. maksup(surf, x, y, z, xs, ys, zs, power)
  25. char    *surf;
  26. double    x, y, z, xs, ys, zs, power;
  27. {
  28.     double          maxval;
  29.     Superq       *super;
  30.     Primitive    *prim;
  31.     Object *newobj;
  32.  
  33.     prim = mallocprim();
  34.     prim->surf = find_surface(surf);
  35.     newobj = new_object(NULL, SUPERQ, (char *)prim, (Trans *)NULL);
  36.     prim->type = SUPERQ;
  37.     super = (Superq *) Malloc((unsigned)sizeof(Superq));
  38.     prim->objpnt.p_superq = super;
  39.     super->x = x;
  40.     super->y = y;
  41.     super->z = z;
  42.     super->bounds[LOW][X] = x - xs;
  43.     super->bounds[HIGH][X] = x + xs;
  44.     super->bounds[LOW][Y] = y - ys;
  45.     super->bounds[HIGH][Y] = y + ys;
  46.     super->bounds[LOW][Z] = z - zs;
  47.     super->bounds[HIGH][Z] = z + zs;
  48.  
  49.     super->pow = power;
  50.     maxval = xs;
  51.     if (ys > maxval)
  52.         maxval = ys;
  53.     if (zs > maxval)
  54.         maxval = zs;
  55.     super->a = xs / maxval;
  56.     super->b = ys / maxval;
  57.     super->c = zs / maxval;
  58.     super->r = pow(maxval, power);
  59.     super->err = pow((ERRCONST * maxval), power) - super->r;
  60.  
  61.     return newobj;
  62. }
  63.  
  64. double
  65. intsup(pos, ray, obj)
  66. Vector           *pos, *ray;
  67. Primitive       *obj;
  68. {
  69.     double          s, si, t;
  70.     double          xadj, yadj, zadj;
  71.     double          old, result, p;
  72.     Ray        tmpray;
  73.     Superq        *super;
  74.     extern        unsigned long primtests[];
  75.  
  76.     primtests[SUPERQ]++;
  77.     super = obj->objpnt.p_superq;
  78.     /* find box intersection */
  79.     if (OutOfBounds(pos, super->bounds)) {
  80.         tmpray.pos = *pos;
  81.         tmpray.dir = *ray;
  82.         s = IntBounds(&tmpray, super->bounds);
  83.         if (s <= 0.)
  84.             return 0;
  85.     }
  86.     else
  87.         s = 0.;
  88.  
  89.     xadj = pos->x - super->x;
  90.     yadj = pos->y - super->y;
  91.     zadj = pos->z - super->z;
  92.  
  93.     /* initial solution */
  94.     p = super->pow;
  95.     result = pow(fabs((xadj + ray->x * s) / super->a), p)
  96.         + pow(fabs((yadj + ray->y * s) / super->b), p)
  97.         + pow(fabs((zadj + ray->z * s) / super->c), p)
  98.         - super->r;
  99.  
  100.     si = s;
  101.     s = s + 0.001;
  102.     /* iterative refinement */
  103.     while (result > super->err) {
  104.         old = result;
  105.         result = pow(fabs((xadj + ray->x * s) / super->a), p)
  106.             + pow(fabs((yadj + ray->y * s) / super->b), p)
  107.             + pow(fabs((zadj + ray->z * s) / super->c), p)
  108.             - super->r;
  109.  
  110.         if (result >= old)
  111.             return (0.0);
  112.         t = (result * (s - si)) / (result - old);
  113.         si = s;
  114.         s -= t;
  115.     }
  116.  
  117.     return (s);
  118. }
  119.  
  120. nrmsup(pos, obj, nrm)
  121. Vector           *pos, *nrm;
  122. Primitive       *obj;
  123. {
  124.     double k;
  125.     Superq *super;
  126.  
  127.     super = obj->objpnt.p_superq;
  128.     nrm->x = (pos->x - super->x) / super->a;
  129.     nrm->y = (pos->y - super->y) / super->b;
  130.     nrm->z = (pos->z - super->z) / super->c;
  131.     k = super->pow - 1;
  132.     if (nrm->x > 0)
  133.         nrm->x = pow(nrm->x, k);
  134.     else
  135.         nrm->x = -pow(-nrm->x, k);
  136.     if (nrm->y > 0)
  137.         nrm->y = pow(nrm->y, k);
  138.     else
  139.         nrm->y = -pow(-nrm->y, k);
  140.     if (nrm->z > 0)
  141.         nrm->z = pow(nrm->z, k);
  142.     else
  143.         nrm->z = -pow(-nrm->z, k);
  144. }
  145.  
  146.  
  147. supextent(o, bounds)
  148. Primitive *o;
  149. double bounds[2][3];
  150. {
  151.     Superq *s = o->objpnt.p_superq;
  152.  
  153.     bounds[LOW][X] = s->bounds[LOW][X];
  154.     bounds[HIGH][X] = s->bounds[HIGH][X];
  155.     bounds[LOW][Y] = s->bounds[LOW][Y];
  156.     bounds[HIGH][Y] = s->bounds[HIGH][Y];
  157.     bounds[LOW][Z] = s->bounds[LOW][Z];
  158.     bounds[HIGH][Z] = s->bounds[HIGH][Z];
  159. }
  160.