home *** CD-ROM | disk | FTP | other *** search
/ Borland Programmer's Resource / Borland_Programmers_Resource_CD_1995.iso / code / graphics / rtnews / rtnews7 < prev    next >
Encoding:
Text File  |  1995-05-19  |  23.2 KB  |  622 lines
  1.  _ __                 ______                         _ __
  2. ' )  )                  /                           ' )  )
  3.  /--' __.  __  ,     --/ __  __.  _. o ____  _,      /  / _  , , , _
  4. /  \_(_/|_/ (_/_    (_/ / (_(_/|_(__<_/ / <_(_)_    /  (_</_(_(_/_/_)_
  5.              /                               /|
  6.             '                               |/
  7.  
  8.             "Light Makes Right"
  9.  
  10.                May 12, 1989
  11.                 Volume 2, Number 3
  12.  
  13. Compiled by Eric Haines, 3D/Eye Inc, 2359 Triphammer Rd, Ithaca, NY 14850
  14.     607-257-1381, hpfcla!hpfcrs!eye!erich@hplabs.hp.com
  15. All contents are US copyright (c) 1989 by the individual authors
  16.  
  17. Contents:
  18.     Introduction
  19.     New People (Carl Bass, Paul Wanuga)
  20.     QRT Ray Tracer (and five other Amiga Ray Tracers) (Steve Koren)
  21.     New Version of MTV Ray Tracer (Mark VandeWettering)
  22.     Minimal Sphere Containing Three Points (Earl Culham)
  23.     Noise and Turbulence Function Code, Pascal and C (Jon Buller,
  24.     William Dirks)
  25.  
  26. -----------------------------------------------------------------------------
  27.  
  28. Introduction
  29.  
  30. Well, we're now at the point where there are six ray tracers for the Amiga.
  31. Interestingly, none of them have implicit efficiency schemes (i.e. where the
  32. user does not have to intervene and create the efficiency structure himself).
  33. Admittedly, efficiency schemes are more code, but I've found that I was getting
  34. a factor of three speed up for a simple scene (a ten ball sphereflake) by using
  35. an efficiency scheme vs. not using one.  When your computer is the speed of an
  36. Amiga, efficiency schemes become vital.
  37.  
  38. Next time I'll include "Tracing Tricks", an article I "edited" for the latest
  39. (and last) "Introduction to Ray Tracing" course notes.  The article is a "best
  40. of the RT News" compendium of efficiency tricks.  By the way, the course notes
  41. should be quite a bargain: they'll consist of the book of our notes by Academic
  42. Press, plus some new tidbits and reprints of "classic" articles.
  43.  
  44. I would like to put the "Ray Tracing News" back issues somewhere that people
  45. can FTP them.  Personally I don't have a computer that has an FTP site, so if
  46. there are any volunteers that would like to store the back issues, please
  47. contact me.  The entire archive is about 448K at this point (not including this
  48. issue), broken into 5 parts.  Can you volunteer?
  49.  
  50. -----------------------------------------------------------------------------
  51.  
  52. # Carl Bass - hybrid shading models and quick(er) hidden surface methods
  53. # Ithaca Software
  54. # 902 West Seneca Street
  55. # Ithaca, NY 14850
  56. # 607-273-3690
  57. alias    carl_bass    carl@mssun1.msi.cornell.edu
  58.  
  59. Carl is the co-founder of Ithaca Software Inc (once upon a time called "Flying
  60. Moose Inc"), which sells the HOOPS package for all kinds of computers.  This is
  61. an object-oriented system which I don't know much about beyond that their
  62. debugger is called WHOOPS.
  63.  
  64. --------
  65.  
  66. #
  67. # Paul Wanuga
  68. # Cornell Program of Computer Graphics
  69. # 120 Rand Hall
  70. # Ithaca, NY 14853
  71. # (607)-255-4880
  72. alias    paul_wanuga    phw@love.tn.cornell.edu
  73. Erich:
  74.  
  75.      Could you please include me in your list of wiz-bango ray tracers?  It
  76. appears Don has me slated for research in ray-tracing complex, realistic,
  77. non-procedural environments.
  78.  
  79. -----------------------------------------------------------------------------
  80.  
  81. QRT Ray Tracer (and five other Amiga Ray Tracers), Steve Koren
  82.  
  83.     This package appeared on comp.graphics a few months back.  I believe
  84. the latest version of the package is available from Mark VandeWettering's FTP
  85. site (see next article).  Write to Steve for more info:
  86.  
  87.     hpfela!koren@hplabs.hp.com
  88.  
  89. The software is written in C and worked just fine on my system.  Below is an
  90. excerpt of the UserMan.doc file of the QRT system (which Steve extensively
  91. documented).
  92.  
  93.  
  94. QRT is a ray tracing  image  rendering  system  that runs under a
  95. variety  of  operating  systems.   It  has  a  free  format input
  96. language   with   extensive   error   detection   and   reporting
  97. capabilities.
  98.     
  99. QRT was developed on the Amiga  personal  computer, so it will
  100. be compared to  other  Amiga  ray  tracers.  There  are, to my
  101. knowledge, five other  Amiga  ray  tracers,  each with its own
  102. strengths  and  weaknesses.    I  will  describe  each  system
  103. briefly, and compare it to QRT.  All the Amiga ray tracers can
  104. operate in HAM (4096 color) mode.
  105.  
  106.   RT: RT was the first ray  tracer  written for the Amiga, by
  107.       Eric    Graham.  It will model  a universe made of only
  108.       spheres, a   sky, and a checkered  or solid ground.  It
  109.       is  relatively  fast,  but  not  generally  useful  for
  110.       realistic modeling   because  of the sphere limitation.
  111.       The input language  is  cryptic,  although  some  error
  112.       checking is done.  The  system  will  only generate low
  113.       resolution images.
  114.  
  115.  SILVER: I have never seen SILVER, so I cannot say much about
  116.       this system.
  117.  
  118.  SCULPT-4D: This package  incorporates  an interactive editor
  119.       for  creating  objects,   and  is  capable  of  quickly
  120.       generating a preliminary  image  of  the scene by using
  121.       hidden surface techniques.  However, every primitive is
  122.       made of polygons, and some  primitives  such as spheres
  123.       require hundreds of  polygons  for a smooth texture, so
  124.       the ray tracing is very slow.   Also, the package takes
  125.       a large amount of memory to run, and is prone to system
  126.       crashes.  Its chief  feature  is  the ability to create
  127.       arbitrary shaped objects  using  a series of triangles.
  128.       Mirrored, dull, or shiny objects are supported.
  129.  
  130.  CLIGHT: This ray tracer also has  an interactive editor, but
  131.       produces very poor quality  images.   It is not capable
  132.       of any patterning or reflection characteristics.
  133.  
  134.  DBW: This is possibly the most  complete  ray tracer for the
  135.       Amiga.  It will support objects  with arbitrary degrees
  136.       of reflection and gloss,  depth  of field effects, some
  137.       texturing,   wavy   surfaces,   fractals,   transparent
  138.       surfaces, diffuse  propagation  of light from object to
  139.       object,  and  5  primitive   types  (sphere,  triangle,
  140.       parallelogram, fractal, and ring).  The input language,
  141.       however,   is    so    cryptic    as   to   be   nearly
  142.       incomprehensible, and  if  there  is  any  error in the
  143.       input file,  it  will  crash  the  system.   It is also
  144.       painfully slow;  some  images  take  16  to 24 hours to
  145.       complete.
  146.  
  147. QRT is meant to be a  compromise  between the fast, simple ray
  148. tracers and the slow powerful systems.   It compares favorably
  149. in speed to RT, and in power  to  Sculpt-3d  or DBW.  It has a
  150. very friendly input language  with  extensive  error checking.
  151. Here are some features of QRT:
  152.  
  153.    o  Multiple primitive types, including user defined
  154.       quadratic surfaces
  155.  
  156.    o  Arbitrary levels of diffuse reflection, spectral
  157.       reflection, transmission, ambient lighting, and gloss
  158.  
  159.    o  User defined pattern information for objects
  160.  
  161.    o  Bounding boxes for groups of objects
  162.  
  163.    o  Shadows
  164.  
  165.    o  Multiple light sources with different characteristics
  166.  
  167.    o  Arbitrary Phong spectral reflection coefficients
  168.  
  169.    o  Color dithering to increase the apparent number of
  170.       colors
  171.  
  172.    o  Easy to use, free format input language with error
  173.       checking.  Parameters are by keyword and may appear in
  174.       any order.
  175.  
  176.    o  Supports medium resolution (128k dots/screen)
  177.  
  178. -----------------------------------------------------------------------------
  179.  
  180. New Version of MTV Ray Tracer, Mark VandeWettering
  181.  
  182. Mark VandeWettering's nice little ray tracer (polygons, spheres, cones and
  183. cylinders, Kay/Kajiya efficiency scheme, yacc/lex parser for NFF format,
  184. otherwise written in C) was released on USENET in three parts on March 27.
  185. Others have interesting features, but the selection of primitives and the speed
  186. of the code of the MTV ray tracer is a big plus.  It's currently my favorite
  187. public domain ray tracer (the amazing BRL package I consider private).
  188.  
  189. The package is available at the usual comp.sources.unix archive sites or from
  190. Mark via anonymous ftp at drizzle.cs.uoregon.edu.  Mark's address is:
  191.  
  192.     markv@drizzle.cs.uoregon.edu
  193.  
  194. -----------------------------------------------------------------------------
  195.  
  196. Subject: Re: Ray Traced Bounding Spheres
  197. From: ECULHAM@UALTAVM.BITNET (Earl Culham)
  198. Organization: University of Alberta VM/CMS
  199.  
  200. In article <17241@versatc.UUCP>, ritter@versatc.UUCP (Jack Ritter) writes:
  201.  
  202. >Given a cluster of points in 3 space, is there
  203. >a good method for finding the minimum radius
  204. >sphere which encloses all the points?  If not
  205. >minimum, at least "small"?  Certainly it should
  206. >be tighter than the sphere which encloses the
  207. >minimum bounding box.
  208. >
  209. >I have a feeling the solution is iterative. If
  210. >so, I could provide a good initial guess for the
  211. >center & radius.
  212.  
  213. The solution is not iterative. A simple solution is available in a two step
  214. process, and is characterized by whether three, or only two of the points are
  215. on the surface of the optimal sphere.
  216.  
  217. Given the points, A, B, and C.
  218. Searching for the center of the optimal sphere, P, with the smallest
  219. radius, R.
  220.  
  221. Checking the midpoints.
  222.  
  223. set P = point halfway between A and B
  224. set R = 1/2 of length from A to B
  225. If length from C to P is less than or equal to R ---> done
  226. Repeat with line A-C, and B-C if needed.
  227.  
  228. Extending the midpoints at right angles.
  229.  
  230. build the line which intersects A-B at a right angle through the midpoint
  231.   of A-B, call the line D.
  232. build the line which intersects A-C at a right angle through the midpoint
  233.   of A-C, call it E.
  234. P is the point where D and E intersect. (Solve the simultaneous equations;
  235.   R is the length A-P)
  236.  
  237. -----------------------------------------------------------------------------
  238.  
  239. From: bullerj@handel.colostate.edu (Jon Buller)
  240. Subject: Re: Pixar's noise function
  241. Organization: Colorado State University, Ft. Collins CO 80523
  242.  
  243. In article <3599@pixar.UUCP> aaa@pixar.UUCP (Tony Apodaca) writes:
  244. >In article <2553@ssc-vax.UUCP> coy@ssc-vax.UUCP (Stephen B Coy) writes:
  245. >>          ...My question:  Does anyone out there know what this
  246. >>noise function really is?
  247. >
  248. >    Three-dimensional simulated natural textures using pseudorandom
  249. >functions were simultaneously and independently developed by Darwyn
  250. >Peachey and Ken Perlin in 1984-5.  Both have papers in the 1985
  251. >Siggraph Proceedings describing their systems.  Perlin, in particular,
  252. >describes in detail how noise() is implemented and can be used creatively
  253.  
  254.    [A description of the properties of the noise function goes here...]
  255.  
  256. >    If you have ever been interested in realistic computer graphics, do
  257. >whatever it takes to get a look at Perlin's paper.  In 1985, his pictures
  258. >were absolutely astounding.  In 1989, they are STILL astounding.
  259.  
  260.     No kidding, some of those pictures are INCREDIBLE.
  261.  
  262.      Here is my code for a look-alike to the Pixar Noise function, and while I
  263. can't say anything about exactly what Pixar's looks like, I think this is
  264. probably close.  After reading the 1985 SIGGRAPH papers on 3d texturing, (and
  265. seeing my prof's code to do a similar thing) I wrote this.  It uses a quadratic
  266. B-Spline instead of the cubic Hermite interpolant implied in the paper.  Also
  267. note that DNoise is just the x, y, and z derivatives of Noise (which are also
  268. B-Splines).  The hashing functions are from Knuth's Art of Computer Programming
  269. (I don't remember which volume though).
  270.  
  271. I know the code is Pascal, and all of you will hate it, but I believe I write
  272. better Pascal than C...  One final note, this was Lightspeed Pascal 2.0 for the
  273. Macintosh, but things have been reformatted slightly to get it on the net.  I
  274. hope this is what you all wanted.
  275.  
  276. --------
  277.  
  278. More:
  279.  
  280. One other thing you might notice, Noise is C1 continuous, DNoise is only C0.
  281. This means that DNoise will have creases in it (along the planes of the random
  282. grid.  To see this, crank out a square: 0<X<5, 0<Y<5, Z=0.  You will see smooth
  283. regions within each unit square, and creases between squares.  To avoid this,
  284. use a cubic B-Spline, or cubic Hermite (as hinted to in the SIGGRAPH
  285. proceedings) the problem there, is that you either need more data points (64
  286. instead of 27) for the B-Spline, or derivative info at each point of the grid
  287. (a normal plane, 4 floats instead of 1).  This would take too muh time for me
  288. to code up to be worth it, and would probably run too much slower (10min for a
  289. 200x200 pixel picture now, ug.)  If somebody wants to give me a cray-3 to play
  290. with, I'll write more accurate (and slower) code, until then... 8-)
  291.  
  292. Jon
  293. bullerj@handel.cs.colostate.edu          ...!ncar!ccncsu!handel!bullerj
  294. (These are my ideas (and code), nobody else SHOULD want these bugs)
  295. I'm just trying to graduate.  Apple, Pixar, HP, etc. take note, I would like
  296. your job offers, I have tired of the university life.
  297.  
  298.  
  299. [NOTE: I have attached the Pascal code for the Turbulence functions.  The Noise
  300. functions are in the next message in "C".  Sorry about the mixed languages, but
  301. I haven't nicely translated these. -- EAH]
  302.  
  303. (* ---------- cut here ---------- cut here ---------- cut here ---------- *)
  304.  
  305. const
  306.    MaxPts = 512;  { Must be 2^n}
  307.    MPM1 = MaxPts - 1;
  308. type
  309.    PtsTyp = array[0..MaxPts] of Extended;
  310. var
  311.    Points: PtsTyp;
  312.  
  313.  
  314. function Turb (Size: Integer;
  315.                ScaleFactor: Extended;
  316.                Loc: Vect;
  317.                Pts: PtsTyp): Extended;
  318.    var
  319.       Scale, Result: Extended;
  320.       Cur: Integer;
  321. begin
  322.    Result := Noise(Loc, Pts);
  323.    Scale := 1.0;
  324.  
  325.    Cur := 1;
  326.    while Cur < Size do begin
  327.       Cur := BSL(Cur, 1);          {Cur := Cur * 2}
  328.       Scale := Scale * ScaleFactor;
  329.       Loc := Scale_Vect(2.0, Loc);
  330.       Result := Result + Noise(Loc, Pts) * Scale;
  331.    end;
  332.    Turb := Result;
  333. end;
  334.  
  335.  
  336. function DTurb (Size: Integer;
  337.                 ScaleFactor: Extended;
  338.                 Loc: Vect;
  339.                 Pts: PtsTyp): Vect;
  340.    var
  341.       Result: Vect;
  342.       Scale: Extended;
  343.       Cur: Integer;
  344. begin
  345.    Result := DNoise(Loc, Pts);
  346.    Scale := 1.0;
  347.  
  348.    Cur := 1;
  349.    while Cur < Size do begin
  350.       Cur := BSL(Cur, 1);
  351.       Scale := Scale * ScaleFactor;
  352.       Loc := Scale_Vect(2.0, Loc);
  353.       Result := Add_Vect(Result, Scale_Vect(Scale, DNoise(Loc, Pts)));
  354.    end;
  355.    DTurb := Result;
  356. end;
  357.  
  358. -----------------------------------------------------------------------------
  359.  
  360. And the C Version...
  361.  
  362. From: dirks@oak.cis.ohio-state.edu (william r dirks)
  363. Subject: C Versions of Noise and DNoise Routines
  364. Organization: Ohio State University Computer and Information Science
  365.  
  366.  
  367. It was suggested to me that some of you would be interested in my 
  368. translated-into-C-and-corrected versions of noise() and Dnoise().
  369.  
  370. So, here they are.  
  371.  
  372. Note that this is only noise and Dnoise.  The turbulence routines are
  373. not included 'cause I haven't translated them yet.
  374.  
  375. Oh yeah, initnoise() fills the pts[] array with random numbers between
  376. 0 and 1.  Don't forget to initialize, or noise will always return 0.
  377. (That's been experimentally verified! :-))
  378.  
  379. --Bill--
  380.  
  381. [Note that you should look over the rand() function if you use this stuff.
  382. Some rand()'s need initialization (srand()), and some give numbers from 0
  383. to 32767, and so should use this as a divisor. -- EAH]
  384.  
  385.  
  386. ---------cut-here------------------------------------cut-here--------
  387. /*
  388. **    Noise and Dnoise routines
  389. *
  390. *    Many thanks to Jon Buller of Colorado State University for these
  391. *    routines.
  392. */
  393.  
  394.  
  395. typedef struct vector {
  396.    double x, y, z;
  397. } Vector;
  398.  
  399.  
  400. #define NUMPTS    512
  401. #define P1    173
  402. #define P2    263
  403. #define P3    337
  404. #define phi    0.6180339
  405.  
  406.  
  407. static double pts[NUMPTS];
  408.  
  409.  
  410. void initnoise()
  411. {
  412.    int i;
  413.    
  414.    for (i = 0; i < NUMPTS; ++i)
  415.       pts[i] = rand() / 2.147483e9;
  416. }
  417.  
  418.  
  419. double noise(p)
  420. Vector *p;
  421. {
  422.    int xi, yi, zi;
  423.    int xa, xb, xc, ya, yb, yc, za, zb, zc;
  424.    double xf, yf, zf;
  425.    double x2, x1, x0, y2, y1, y0, z2, z1, z0;
  426.    double p000, p100, p200, p010, p110, p210, p020, p120, p220;
  427.    double p001, p101, p201, p011, p111, p211, p021, p121, p221;
  428.    double p002, p102, p202, p012, p112, p212, p022, p122, p222;
  429.  
  430.    xi = floor(p->x);
  431.    xa = floor(P1 * (xi * phi - floor(xi * phi)));
  432.    xb = floor(P1 * ((xi + 1) * phi - floor((xi + 1) * phi)));
  433.    xc = floor(P1 * ((xi + 2) * phi - floor((xi + 2) * phi)));
  434.  
  435.    yi = floor(p->y);
  436.    ya = floor(P2 * (yi * phi - floor(yi * phi)));
  437.    yb = floor(P2 * ((yi + 1) * phi - floor((yi + 1) * phi)));
  438.    yc = floor(P2 * ((yi + 2) * phi - floor((yi + 2) * phi)));
  439.  
  440.    zi = floor(p->z);
  441.    za = floor(P3 * (zi * phi - floor(zi * phi)));
  442.    zb = floor(P3 * ((zi + 1) * phi - floor((zi + 1) * phi)));
  443.    zc = floor(P3 * ((zi + 2) * phi - floor((zi + 2) * phi)));
  444.  
  445.    p000 = pts[xa + ya + za & NUMPTS - 1];
  446.    p100 = pts[xb + ya + za & NUMPTS - 1];
  447.    p200 = pts[xc + ya + za & NUMPTS - 1];
  448.    p010 = pts[xa + yb + za & NUMPTS - 1];
  449.    p110 = pts[xb + yb + za & NUMPTS - 1];
  450.    p210 = pts[xc + yb + za & NUMPTS - 1];
  451.    p020 = pts[xa + yc + za & NUMPTS - 1];
  452.    p120 = pts[xb + yc + za & NUMPTS - 1];
  453.    p220 = pts[xc + yc + za & NUMPTS - 1];
  454.    p001 = pts[xa + ya + zb & NUMPTS - 1];
  455.    p101 = pts[xb + ya + zb & NUMPTS - 1];
  456.    p201 = pts[xc + ya + zb & NUMPTS - 1];
  457.    p011 = pts[xa + yb + zb & NUMPTS - 1];
  458.    p111 = pts[xb + yb + zb & NUMPTS - 1];
  459.    p211 = pts[xc + yb + zb & NUMPTS - 1];
  460.    p021 = pts[xa + yc + zb & NUMPTS - 1];
  461.    p121 = pts[xb + yc + zb & NUMPTS - 1];
  462.    p221 = pts[xc + yc + zb & NUMPTS - 1];
  463.    p002 = pts[xa + ya + zc & NUMPTS - 1];
  464.    p102 = pts[xb + ya + zc & NUMPTS - 1];
  465.    p202 = pts[xc + ya + zc & NUMPTS - 1];
  466.    p012 = pts[xa + yb + zc & NUMPTS - 1];
  467.    p112 = pts[xb + yb + zc & NUMPTS - 1];
  468.    p212 = pts[xc + yb + zc & NUMPTS - 1];
  469.    p022 = pts[xa + yc + zc & NUMPTS - 1];
  470.    p122 = pts[xb + yc + zc & NUMPTS - 1];
  471.    p222 = pts[xc + yc + zc & NUMPTS - 1];
  472.  
  473.    xf = p->x - xi;
  474.    x1 = xf * xf;
  475.    x2 = 0.5 * x1;
  476.    x1 = 0.5 + xf - x1;
  477.    x0 = 0.5 - xf + x2;
  478.  
  479.    yf = p->y - yi;
  480.    y1 = yf * yf;
  481.    y2 = 0.5 * y1;
  482.    y1 = 0.5 + yf - y1;
  483.    y0 = 0.5 - yf + y2;
  484.  
  485.    zf = p->z - zi;
  486.    z1 = zf * zf;
  487.    z2 = 0.5 * z1;
  488.    z1 = 0.5 + zf - z1;
  489.    z0 = 0.5 - zf + z2;
  490.  
  491.    return   z0 * (y0 * (x0 * p000 + x1 * p100 + x2 * p200) +
  492.                   y1 * (x0 * p010 + x1 * p110 + x2 * p210) +
  493.                   y2 * (x0 * p020 + x1 * p120 + x2 * p220)) +
  494.             z1 * (y0 * (x0 * p001 + x1 * p101 + x2 * p201) +
  495.                   y1 * (x0 * p011 + x1 * p111 + x2 * p211) +
  496.                   y2 * (x0 * p021 + x1 * p121 + x2 * p221)) +
  497.             z2 * (y0 * (x0 * p002 + x1 * p102 + x2 * p202) +
  498.                   y1 * (x0 * p012 + x1 * p112 + x2 * p212) +
  499.                   y2 * (x0 * p022 + x1 * p122 + x2 * p222));
  500. }
  501.  
  502.  
  503.  
  504. Vector Dnoise(p)
  505. Vector *p;
  506. {
  507.    Vector v;
  508.    int xi, yi, zi;
  509.    int xa, xb, xc, ya, yb, yc, za, zb, zc;
  510.    double xf, yf, zf;
  511.    double x2, x1, x0, y2, y1, y0, z2, z1, z0;
  512.    double xd2, xd1, xd0, yd2, yd1, yd0, zd2, zd1, zd0;
  513.    double p000, p100, p200, p010, p110, p210, p020, p120, p220;
  514.    double p001, p101, p201, p011, p111, p211, p021, p121, p221;
  515.    double p002, p102, p202, p012, p112, p212, p022, p122, p222;
  516.  
  517.    xi = floor(p->x);
  518.    xa = floor(P1 * (xi * phi - floor(xi * phi)));
  519.    xb = floor(P1 * ((xi + 1) * phi - floor((xi + 1) * phi)));
  520.    xc = floor(P1 * ((xi + 2) * phi - floor((xi + 2) * phi)));
  521.  
  522.    yi = floor(p->y);
  523.    ya = floor(P2 * (yi * phi - floor(yi * phi)));
  524.    yb = floor(P2 * ((yi + 1) * phi - floor((yi + 1) * phi)));
  525.    yc = floor(P2 * ((yi + 2) * phi - floor((yi + 2) * phi)));
  526.  
  527.    zi = floor(p->z);
  528.    za = floor(P3 * (zi * phi - floor(zi * phi)));
  529.    zb = floor(P3 * ((zi + 1) * phi - floor((zi + 1) * phi)));
  530.    zc = floor(P3 * ((zi + 2) * phi - floor((zi + 2) * phi)));
  531.  
  532.    p000 = pts[xa + ya + za & NUMPTS - 1];
  533.    p100 = pts[xb + ya + za & NUMPTS - 1];
  534.    p200 = pts[xc + ya + za & NUMPTS - 1];
  535.    p010 = pts[xa + yb + za & NUMPTS - 1];
  536.    p110 = pts[xb + yb + za & NUMPTS - 1];
  537.    p210 = pts[xc + yb + za & NUMPTS - 1];
  538.    p020 = pts[xa + yc + za & NUMPTS - 1];
  539.    p120 = pts[xb + yc + za & NUMPTS - 1];
  540.    p220 = pts[xc + yc + za & NUMPTS - 1];
  541.    p001 = pts[xa + ya + zb & NUMPTS - 1];
  542.    p101 = pts[xb + ya + zb & NUMPTS - 1];
  543.    p201 = pts[xc + ya + zb & NUMPTS - 1];
  544.    p011 = pts[xa + yb + zb & NUMPTS - 1];
  545.    p111 = pts[xb + yb + zb & NUMPTS - 1];
  546.    p211 = pts[xc + yb + zb & NUMPTS - 1];
  547.    p021 = pts[xa + yc + zb & NUMPTS - 1];
  548.    p121 = pts[xb + yc + zb & NUMPTS - 1];
  549.    p221 = pts[xc + yc + zb & NUMPTS - 1];
  550.    p002 = pts[xa + ya + zc & NUMPTS - 1];
  551.    p102 = pts[xb + ya + zc & NUMPTS - 1];
  552.    p202 = pts[xc + ya + zc & NUMPTS - 1];
  553.    p012 = pts[xa + yb + zc & NUMPTS - 1];
  554.    p112 = pts[xb + yb + zc & NUMPTS - 1];
  555.    p212 = pts[xc + yb + zc & NUMPTS - 1];
  556.    p022 = pts[xa + yc + zc & NUMPTS - 1];
  557.    p122 = pts[xb + yc + zc & NUMPTS - 1];
  558.    p222 = pts[xc + yc + zc & NUMPTS - 1];
  559.  
  560.    xf = p->x - xi;
  561.    x1 = xf * xf;
  562.    x2 = 0.5 * x1;
  563.    x1 = 0.5 + xf - x1;
  564.    x0 = 0.5 - xf + x2;
  565.    xd2 = xf;
  566.    xd1 = 1.0 - xf - xf;
  567.    xd0 = xf - 1.0;
  568.  
  569.    yf = p->y - yi;
  570.    y1 = yf * yf;
  571.    y2 = 0.5 * y1;
  572.    y1 = 0.5 + yf - y1;
  573.    y0 = 0.5 - yf + y2;
  574.    yd2 = yf;
  575.    yd1 = 1.0 - yf - yf;
  576.    yd0 = yf - 1.0;
  577.  
  578.    zf = p->z - zi;
  579.    z1 = zf * zf;
  580.    z2 = 0.5 * z1;
  581.    z1 = 0.5 + zf - z1;
  582.    z0 = 0.5 - zf + z2;
  583.    zd2 = zf;
  584.    zd1 = 1.0 - zf - zf;
  585.    zd0 = zf - 1.0;
  586.  
  587.    v.x        = z0 * (y0 * (xd0 * p000 + xd1 * p100 + xd2 * p200) +
  588.                       y1 * (xd0 * p010 + xd1 * p110 + xd2 * p210) +
  589.                       y2 * (xd0 * p020 + xd1 * p120 + xd2 * p220)) +
  590.                 z1 * (y0 * (xd0 * p001 + xd1 * p101 + xd2 * p201) +
  591.                       y1 * (xd0 * p011 + xd1 * p111 + xd2 * p211) +
  592.                       y2 * (xd0 * p021 + xd1 * p121 + xd2 * p221)) +
  593.                 z2 * (y0 * (xd0 * p002 + xd1 * p102 + xd2 * p202) +
  594.                       y1 * (xd0 * p012 + xd1 * p112 + xd2 * p212) +
  595.                       y2 * (xd0 * p022 + xd1 * p122 + xd2 * p222));
  596.                                   
  597.    v.y        = z0 * (yd0 * (x0 * p000 + x1 * p100 + x2 * p200) +
  598.                       yd1 * (x0 * p010 + x1 * p110 + x2 * p210) +
  599.                       yd2 * (x0 * p020 + x1 * p120 + x2 * p220)) +
  600.                 z1 * (yd0 * (x0 * p001 + x1 * p101 + x2 * p201) +
  601.                       yd1 * (x0 * p011 + x1 * p111 + x2 * p211) +
  602.                       yd2 * (x0 * p021 + x1 * p121 + x2 * p221)) +
  603.                 z2 * (yd0 * (x0 * p002 + x1 * p102 + x2 * p202) +
  604.                       yd1 * (x0 * p012 + x1 * p112 + x2 * p212) +
  605.                       yd2 * (x0 * p022 + x1 * p122 + x2 * p222));
  606.                                   
  607.    v.z        = zd0 * (y0 * (x0 * p000 + x1 * p100 + x2 * p200) +
  608.                        y1 * (x0 * p010 + x1 * p110 + x2 * p210) +
  609.                        y2 * (x0 * p020 + x1 * p120 + x2 * p220)) +
  610.                 zd1 * (y0 * (x0 * p001 + x1 * p101 + x2 * p201) +
  611.                        y1 * (x0 * p011 + x1 * p111 + x2 * p211) +
  612.                        y2 * (x0 * p021 + x1 * p121 + x2 * p221)) +
  613.                 zd2 * (y0 * (x0 * p002 + x1 * p102 + x2 * p202) +
  614.                        y1 * (x0 * p012 + x1 * p112 + x2 * p212) +
  615.                        y2 * (x0 * p022 + x1 * p122 + x2 * p222));
  616.    return v;
  617. }
  618.  
  619. -----------------------------------------------------------------------------
  620.  
  621. END OF RTNEWS
  622.