home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 3 / Amiga Tools 3.iso / grafik / raytracing / rayshade-4.0.6.3 / libray / libobj / rotspline.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-09  |  6.2 KB  |  275 lines
  1. /*
  2.  * rotspline.c
  3.  *
  4.  * Copyright (C) 1991, Gerald Iles.
  5.  * All rights reserved.
  6.  *
  7.  * This software may be freely copied, modified, and redistributed
  8.  * provided that this copyright notice is preserved on all copies.
  9.  *
  10.  * You may not distribute this software, in whole or in part, as part of
  11.  * any commercial product without the express consent of the authors.
  12.  *
  13.  * There is no warranty or other guarantee of fitness of this software
  14.  * for any purpose.  It is provided solely "as is".
  15.  *
  16.  * rotspline.c,v 4.1 1994/08/09 07:59:58 explorer Exp
  17.  *
  18.  * rotspline.c,v
  19.  * Revision 4.1  1994/08/09  07:59:58  explorer
  20.  * Bump version to 4.1
  21.  *
  22.  * Revision 1.1.1.1  1994/08/08  04:52:10  explorer
  23.  * Initial import.  This is a prerelease of 4.0.6enh3, or 4.1 possibly.
  24.  *
  25.  */
  26. #include "geom.h"
  27. #include "rotspline.h"
  28.  
  29. static Methods *iRotsplineMethods = NULL;
  30. static char rotsplineName[] = "rotspline";
  31.  
  32. unsigned long RotsplineTests, RotsplineHits;
  33.  
  34. Rotspline *
  35. RotsplineCreate(bot, x0, d0, top, xw, dw, coeffs)
  36. Vector *bot, *top;
  37. Float x0, xw, d0, dw, *coeffs;
  38. {
  39.     Rotspline *rotspline;
  40.     Float w;
  41.     Vector axis, tmp;
  42.  
  43.     VecSub(*top, *bot, &axis);
  44.     w = VecNormalize(&axis);
  45.     if (w < EPSILON) {
  46.         RLerror(RL_WARN, "Degenerate rotspline.\n");
  47.         return (Rotspline *)NULL;
  48.     }
  49.  
  50.     rotspline = (Rotspline *)Malloc(sizeof(Rotspline));
  51.  
  52.     rotspline->w = w;
  53.  
  54.     if (coeffs == NULL)
  55.     {
  56.         rotspline->k = (Float *)Calloc(4, sizeof(Float));
  57.         rotspline->k[3] = 2.0 * (x0 * x0 - xw * xw + w * x0 * d0 + w * xw * dw) / (w * w * w);
  58.         rotspline->k[2] = (3.0 * (xw * xw - x0 * x0) - 4.0*w*x0*d0 - 2.0*w*xw*dw) / (w * w);
  59.         rotspline->k[1] = 2.0 * x0 * d0;
  60.         rotspline->k[0] = x0 * x0;
  61.     }
  62.     else
  63.         rotspline->k = coeffs;
  64.     CoordSysTransform(bot, &axis, 1.0, 1.0, &rotspline->trans);
  65.  
  66.     return rotspline;
  67. }
  68.  
  69. Methods *
  70. RotsplineMethods()
  71. {
  72.     if (iRotsplineMethods == (Methods *)NULL) {
  73.         iRotsplineMethods = MethodsCreate();
  74.         iRotsplineMethods->name = RotsplineName;
  75.         iRotsplineMethods->create = (GeomCreateFunc *)RotsplineCreate;
  76.         iRotsplineMethods->methods = RotsplineMethods;
  77.         iRotsplineMethods->intersect = RotsplineIntersect;
  78.         iRotsplineMethods->normal = RotsplineNormal;
  79.         iRotsplineMethods->uv = RotsplineUV;
  80.         iRotsplineMethods->bounds = RotsplineBounds;
  81.         iRotsplineMethods->stats = RotsplineStats;
  82.         iRotsplineMethods->checkbounds = TRUE;
  83.         iRotsplineMethods->closed = FALSE;
  84.     }
  85.     return iRotsplineMethods;
  86. }
  87.  
  88. int
  89. RotsplineIntersect(rotspline, ray, mindist, maxdist)
  90. Rotspline *rotspline;
  91. Ray *ray;
  92. Float mindist, *maxdist;
  93. {
  94.     double c[4], s[3], pp, pd, pz, dz, tmp, dist;
  95.     Ray newray;
  96.     Vector ndir, npos;
  97.     int num, i;
  98.  
  99.     RotsplineTests++;
  100.  
  101.     newray = *ray;
  102.     (void) RayTransform(&newray, &rotspline->trans.itrans);
  103.     ndir = newray.dir;
  104.     npos = newray.pos;
  105.  
  106.     pp = npos.x*npos.x + npos.y*npos.y +npos.z*npos.z;
  107.     pd = npos.x*ndir.x + npos.y*ndir.y +npos.z*ndir.z;
  108.     pz = npos.z;
  109.     dz = ndir.z;
  110.  
  111.     c[3] = rotspline->k[3] * dz * dz * dz;
  112.     c[2] = dz * dz *(3.0 * rotspline->k[3] * pz + rotspline->k[2] + 1.0) - 1.0;
  113.     c[1] = dz *(pz *(3.0 * rotspline->k[3] * pz + 2.0 * rotspline->k[2] + 2.0) + rotspline->k[1]) - 2.0 * pd;
  114.     c[0] = pz *(pz *(rotspline->k[3] * pz + rotspline->k[2] +1.0) + rotspline->k[1]) + rotspline->k[0] - pp;
  115.  
  116.     if (fabs(c[3]) > EPSILON)
  117.         num = SolveCubic(c, s);
  118.     else if (fabs(c[2]) > EPSILON)
  119.         num = SolveQuadric(c, s);
  120.     else if (fabs(c[1]) > EPSILON)
  121.     {
  122.         s[0] = -c[0]/c[1];
  123.         num = 1;
  124.     }
  125.     else
  126.         num = 0;
  127.         
  128.     if (num==0) return FALSE;
  129.     dist = 0.0;
  130.     for(i=0;i<num;i++)
  131.     {
  132.         if (s[i] > mindist)
  133.         {
  134.             tmp = pz + dz * s[i];
  135.             if (tmp > 0.0 && tmp < rotspline->w)
  136.             {
  137.                 if (dist == 0.0)
  138.                     dist = s[i];
  139.                 else if (s[i] < dist)
  140.                     dist = s[i];
  141.             }
  142.         }
  143.     }
  144.     if (dist > mindist && dist < *maxdist) {
  145.         *maxdist = dist;
  146.         RotsplineHits++;
  147.         return TRUE;
  148.     }
  149.     return FALSE;
  150. }
  151.  
  152. int
  153. RotsplineNormal(rotspline, pos, nrm, gnrm)
  154. Rotspline *rotspline;
  155. Vector *pos, *nrm, *gnrm;
  156. {
  157.     Vector npos;
  158.     Float du, r;
  159.  
  160.     npos = *pos;
  161.     PointTransform(&npos, &rotspline->trans.itrans);
  162.     
  163.     du = rotspline->k[3] * 3.0 * npos.z * npos.z
  164.        + rotspline->k[2] * 2.0 * npos.z
  165.        + rotspline->k[1];
  166.  
  167.     r = sqrt(npos.x * npos.x + npos.y * npos.y);
  168.  
  169.     nrm->x = r * npos.x;
  170.     nrm->y = r * npos.y;
  171.     nrm->z = -du / 2.0;
  172.  
  173.     NormalTransform(nrm, &rotspline->trans.itrans);
  174.     *gnrm = *nrm;
  175.     return FALSE;
  176. }
  177.  
  178. void
  179. RotsplineUV(rotspline, pos, norm, uv, dpdu, dpdv)
  180. Rotspline *rotspline;
  181. Vector *pos, *norm, *dpdu, *dpdv;
  182. Vec2d *uv;
  183. {
  184.     Vector npos;
  185.     Float val, d;
  186.  
  187.     npos = *pos;
  188.     PointTransform(&npos, &rotspline->trans.itrans);
  189.  
  190.     uv->v = npos.z / rotspline->w;
  191.     d = sqrt(npos.x * npos.x + npos.y * npos.y);
  192.     if (d < EPSILON)
  193.         uv->u = 0.;
  194.     else
  195.     {
  196.         val = npos.x / d;
  197.         if (val > 1.)
  198.             uv->u = 0.;
  199.         else if (val < -1.)
  200.             uv->u = 0.5;
  201.         else
  202.         {
  203.             uv->u = acos(val) / TWOPI;
  204.             if (npos.y < 0.)
  205.                 uv->u = 1. - uv->u;
  206.         }
  207.     }
  208.  
  209.     if (dpdu)
  210.     {
  211.         dpdu->x = -npos.y;
  212.         dpdu->y = npos.x;
  213.         dpdu->z = 0.;
  214.         dpdv->x = rotspline->k[3] * 3.0 * npos.z * npos.z
  215.                 + rotspline->k[2] * 2.0 * npos.z
  216.                 + rotspline->k[1];
  217.         dpdv->y = 0.;
  218.         dpdv->z = (d < EPSILON)?0.:d;
  219.         VecTransform(dpdu, &rotspline->trans.trans);
  220.         VecTransform(dpdv, &rotspline->trans.trans);
  221.         (void)VecNormalize(dpdu);
  222.         (void)VecNormalize(dpdv);
  223.     }
  224. }
  225.  
  226. void
  227. RotsplineBounds(rotspline, bounds)
  228. Rotspline *rotspline;
  229. Float bounds[2][3];
  230. {
  231.     Float r0, rw, rmax, d, z, rz;
  232.  
  233.     r0 = (rotspline->k[0] < EPSILON ? 0.0 : sqrt(rotspline->k[0]));
  234.     rw = EvalSpline(rotspline->k, rotspline->w);
  235.     rw = (rw < EPSILON ? 0.0 : sqrt(rw));
  236.     rmax = max(r0, rw);
  237.     d = 4.*rotspline->k[2]*rotspline->k[2] - 12.*rotspline->k[3]*rotspline->k[1];
  238.     if (d > 0.)
  239.     {
  240.         z = (-2. * rotspline->k[2] + sqrt(d))/(6. * rotspline->k[3]);
  241.         if (z > 0. && z < rotspline->w)
  242.         {
  243.             rz = EvalSpline(rotspline->k, z);
  244.             rz = (rz < EPSILON ? 0.0 : sqrt(rz));
  245.             rmax = max(rmax, rz);
  246.         }
  247.         z = (-2. * rotspline->k[2] - sqrt(d))/(6. * rotspline->k[3]);
  248.         if (z > 0. && z < rotspline->w)
  249.         {
  250.             rz = EvalSpline(rotspline->k, z);
  251.             rz = (rz < EPSILON ? 0.0 : sqrt(rz));
  252.             rmax = max(rmax, rz);
  253.         }
  254.     }
  255.     bounds[LOW][X] = bounds[LOW][Y] = -rmax;
  256.     bounds[HIGH][X] = bounds[HIGH][Y] = rmax;
  257.     bounds[LOW][Z] = 0.;
  258.     bounds[HIGH][Z] = rotspline->w;
  259.     BoundsTransform(&rotspline->trans.trans, bounds);
  260. }
  261.  
  262. char *
  263. RotsplineName()
  264. {
  265.     return rotsplineName;
  266. }
  267.  
  268. void
  269. RotsplineStats(tests, hits)
  270. unsigned long *tests, *hits;
  271. {
  272.     *tests = RotsplineTests;
  273.     *hits = RotsplineHits;
  274. }
  275.