Path: utzoo!utgpu!jarvis.csri.toronto.edu!rutgers!mit-eddie!uw-beaver!cornell!vax5!cnsy
From: cnsy@vax5.CIT.CORNELL.EDU
Newsgroups: comp.graphics
Subject: Ray Tracing News archive 7 of 7 (whew!)
Summary: echoes of echoes...
Message-ID: <18718@vax5.CIT.CORNELL.EDU>
Date: 1 Jun 89 14:57:57 GMT
Sender: news@vax5.CIT.CORNELL.EDU
Reply-To: cnsy@vax5.cit.cornell.edu (Eric Haines)
Organization: 3D/Eye Inc, Ithaca, NY
Lines: 632
References:

Well, now you're up to date.  Please forgive all the repeats of USENET
postings--they should be of interest to newcomers to comp.graphics, at least.
If you're interested in contributing and subscribing to the RT News, drop
me a line (notes to this Cornell account risk being ignored--please write to
the email address listed below).  I'll continue to post issues here, so if
you are not interested in contributing, please don't ask to subscribe!

Enjoy,

Eric Haines

 _ __                 ______                         _ __
' )  )                  /                           ' )  )
 /--' __.  __  ,     --/ __  __.  _. o ____  _,      /  / _  , , , _
/  \_(_/|_/ (_/_    (_/ / (_(_/|_(__<_/ / <_(_)_    /  (_</_(_(_/_/_)_
             /                               /|
            '                               |/

			"Light Makes Right"

			   May 12, 1989
		        Volume 2, Number 3

Compiled by Eric Haines, 3D/Eye Inc, 2359 Triphammer Rd, Ithaca, NY 14850
    607-257-1381, hpfcla!hpfcrs!eye!erich@hplabs.hp.com
All contents are US copyright (c) 1989 by the individual authors

Contents:
    Introduction
    New People (Carl Bass, Paul Wanuga)
    QRT Ray Tracer (and five other Amiga Ray Tracers) (Steve Koren)
    New Version of MTV Ray Tracer (Mark VandeWettering)
    Minimal Sphere Containing Three Points (Earl Culham)
    Noise and Turbulence Function Code, Pascal and C (Jon Buller,
	William Dirks)

-----------------------------------------------------------------------------

Introduction

Well, we're now at the point where there are six ray tracers for the Amiga.
Interestingly, none of them have implicit efficiency schemes (i.e. where the
user does not have to intervene and create the efficiency structure himself).
Admittedly, efficiency schemes are more code, but I've found that I was getting
a factor of three speed up for a simple scene (a ten ball sphereflake) by using
an efficiency scheme vs. not using one.  When your computer is the speed of an
Amiga, efficiency schemes become vital.

Next time I'll include "Tracing Tricks", an article I "edited" for the latest
(and last) "Introduction to Ray Tracing" course notes.  The article is a "best
of the RT News" compendium of efficiency tricks.  By the way, the course notes
should be quite a bargain: they'll consist of the book of our notes by Academic
Press, plus some new tidbits and reprints of "classic" articles.

I would like to put the "Ray Tracing News" back issues somewhere that people
can FTP them.  Personally I don't have a computer that has an FTP site, so if
there are any volunteers that would like to store the back issues, please
contact me.  The entire archive is about 448K at this point (not including this
issue), broken into 5 parts.  Can you volunteer?

-----------------------------------------------------------------------------

# Carl Bass - hybrid shading models and quick(er) hidden surface methods
# Ithaca Software
# 902 West Seneca Street
# Ithaca, NY 14850
# 607-273-3690
alias	carl_bass	carl@mssun1.msi.cornell.edu

Carl is the co-founder of Ithaca Software Inc (once upon a time called "Flying
Moose Inc"), which sells the HOOPS package for all kinds of computers.  This is
an object-oriented system which I don't know much about beyond that their
debugger is called WHOOPS.

--------

#
# Paul Wanuga
# Cornell Program of Computer Graphics
# 120 Rand Hall
# Ithaca, NY 14853
# (607)-255-4880
alias	paul_wanuga	phw@love.tn.cornell.edu
Erich:
 
     Could you please include me in your list of wiz-bango ray tracers?  It
appears Don has me slated for research in ray-tracing complex, realistic,
non-procedural environments.

-----------------------------------------------------------------------------

QRT Ray Tracer (and five other Amiga Ray Tracers), Steve Koren

	This package appeared on comp.graphics a few months back.  I believe
the latest version of the package is available from Mark VandeWettering's FTP
site (see next article).  Write to Steve for more info:

	hpfela!koren@hplabs.hp.com

The software is written in C and worked just fine on my system.  Below is an
excerpt of the UserMan.doc file of the QRT system (which Steve extensively
documented).


QRT is a ray tracing  image  rendering  system  that runs under a
variety  of  operating  systems.   It  has  a  free  format input
language   with   extensive   error   detection   and   reporting
capabilities.
    
QRT was developed on the Amiga  personal  computer, so it will
be compared to  other  Amiga  ray  tracers.  There  are, to my
knowledge, five other  Amiga  ray  tracers,  each with its own
strengths  and  weaknesses.    I  will  describe  each  system
briefly, and compare it to QRT.  All the Amiga ray tracers can
operate in HAM (4096 color) mode.

  RT: RT was the first ray  tracer  written for the Amiga, by
      Eric    Graham.  It will model  a universe made of only
      spheres, a   sky, and a checkered  or solid ground.  It
      is  relatively  fast,  but  not  generally  useful  for
      realistic modeling   because  of the sphere limitation.
      The input language  is  cryptic,  although  some  error
      checking is done.  The  system  will  only generate low
      resolution images.

 SILVER: I have never seen SILVER, so I cannot say much about
      this system.

 SCULPT-4D: This package  incorporates  an interactive editor
      for  creating  objects,   and  is  capable  of  quickly
      generating a preliminary  image  of  the scene by using
      hidden surface techniques.  However, every primitive is
      made of polygons, and some  primitives  such as spheres
      require hundreds of  polygons  for a smooth texture, so
      the ray tracing is very slow.   Also, the package takes
      a large amount of memory to run, and is prone to system
      crashes.  Its chief  feature  is  the ability to create
      arbitrary shaped objects  using  a series of triangles.
      Mirrored, dull, or shiny objects are supported.

 CLIGHT: This ray tracer also has  an interactive editor, but
      produces very poor quality  images.   It is not capable
      of any patterning or reflection characteristics.

 DBW: This is possibly the most  complete  ray tracer for the
      Amiga.  It will support objects  with arbitrary degrees
      of reflection and gloss,  depth  of field effects, some
      texturing,   wavy   surfaces,   fractals,   transparent
      surfaces, diffuse  propagation  of light from object to
      object,  and  5  primitive   types  (sphere,  triangle,
      parallelogram, fractal, and ring).  The input language,
      however,   is    so    cryptic    as   to   be   nearly
      incomprehensible, and  if  there  is  any  error in the
      input file,  it  will  crash  the  system.   It is also
      painfully slow;  some  images  take  16  to 24 hours to
      complete.

QRT is meant to be a  compromise  between the fast, simple ray
tracers and the slow powerful systems.   It compares favorably
in speed to RT, and in power  to  Sculpt-3d  or DBW.  It has a
very friendly input language  with  extensive  error checking.
Here are some features of QRT:

   o  Multiple primitive types, including user defined
      quadratic surfaces

   o  Arbitrary levels of diffuse reflection, spectral
      reflection, transmission, ambient lighting, and gloss

   o  User defined pattern information for objects

   o  Bounding boxes for groups of objects

   o  Shadows

   o  Multiple light sources with different characteristics

   o  Arbitrary Phong spectral reflection coefficients

   o  Color dithering to increase the apparent number of
      colors

   o  Easy to use, free format input language with error
      checking.  Parameters are by keyword and may appear in
      any order.

   o  Supports medium resolution (128k dots/screen)

-----------------------------------------------------------------------------

New Version of MTV Ray Tracer, Mark VandeWettering

Mark VandeWettering's nice little ray tracer (polygons, spheres, cones and
cylinders, Kay/Kajiya efficiency scheme, yacc/lex parser for NFF format,
otherwise written in C) was released on USENET in three parts on March 27.
Others have interesting features, but the selection of primitives and the speed
of the code of the MTV ray tracer is a big plus.  It's currently my favorite
public domain ray tracer (the amazing BRL package I consider private).

The package is available at the usual comp.sources.unix archive sites or from
Mark via anonymous ftp at drizzle.cs.uoregon.edu.  Mark's address is:

    markv@drizzle.cs.uoregon.edu

-----------------------------------------------------------------------------

Subject: Re: Ray Traced Bounding Spheres
From: ECULHAM@UALTAVM.BITNET (Earl Culham)
Organization: University of Alberta VM/CMS

In article <17241@versatc.UUCP>, ritter@versatc.UUCP (Jack Ritter) writes:
 
>Given a cluster of points in 3 space, is there
>a good method for finding the minimum radius
>sphere which encloses all the points?  If not
>minimum, at least "small"?  Certainly it should
>be tighter than the sphere which encloses the
>minimum bounding box.
>
>I have a feeling the solution is iterative. If
>so, I could provide a good initial guess for the
>center & radius.
 
The solution is not iterative. A simple solution is available in a two step
process, and is characterized by whether three, or only two of the points are
on the surface of the optimal sphere.
 
Given the points, A, B, and C.
Searching for the center of the optimal sphere, P, with the smallest
radius, R.
 
Checking the midpoints.
 
set P = point halfway between A and B
set R = 1/2 of length from A to B
If length from C to P is less than or equal to R ---> done
Repeat with line A-C, and B-C if needed.
 
Extending the midpoints at right angles.
 
build the line which intersects A-B at a right angle through the midpoint
  of A-B, call the line D.
build the line which intersects A-C at a right angle through the midpoint
  of A-C, call it E.
P is the point where D and E intersect. (Solve the simultaneous equations;
  R is the length A-P)

-----------------------------------------------------------------------------

From: bullerj@handel.colostate.edu (Jon Buller)
Subject: Re: Pixar's noise function
Organization: Colorado State University, Ft. Collins CO 80523

In article <3599@pixar.UUCP> aaa@pixar.UUCP (Tony Apodaca) writes:
>In article <2553@ssc-vax.UUCP> coy@ssc-vax.UUCP (Stephen B Coy) writes:
>>          ...My question:  Does anyone out there know what this
>>noise function really is?
>
>	Three-dimensional simulated natural textures using pseudorandom
>functions were simultaneously and independently developed by Darwyn
>Peachey and Ken Perlin in 1984-5.  Both have papers in the 1985
>Siggraph Proceedings describing their systems.  Perlin, in particular,
>describes in detail how noise() is implemented and can be used creatively

   [A description of the properties of the noise function goes here...]

>	If you have ever been interested in realistic computer graphics, do
>whatever it takes to get a look at Perlin's paper.  In 1985, his pictures
>were absolutely astounding.  In 1989, they are STILL astounding.

	No kidding, some of those pictures are INCREDIBLE.

     Here is my code for a look-alike to the Pixar Noise function, and while I
can't say anything about exactly what Pixar's looks like, I think this is
probably close.  After reading the 1985 SIGGRAPH papers on 3d texturing, (and
seeing my prof's code to do a similar thing) I wrote this.  It uses a quadratic
B-Spline instead of the cubic Hermite interpolant implied in the paper.  Also
note that DNoise is just the x, y, and z derivatives of Noise (which are also
B-Splines).  The hashing functions are from Knuth's Art of Computer Programming
(I don't remember which volume though).

I know the code is Pascal, and all of you will hate it, but I believe I write
better Pascal than C...  One final note, this was Lightspeed Pascal 2.0 for the
Macintosh, but things have been reformatted slightly to get it on the net.  I
hope this is what you all wanted.

--------

More:

One other thing you might notice, Noise is C1 continuous, DNoise is only C0.
This means that DNoise will have creases in it (along the planes of the random
grid.  To see this, crank out a square: 0<X<5, 0<Y<5, Z=0.  You will see smooth
regions within each unit square, and creases between squares.  To avoid this,
use a cubic B-Spline, or cubic Hermite (as hinted to in the SIGGRAPH
proceedings) the problem there, is that you either need more data points (64
instead of 27) for the B-Spline, or derivative info at each point of the grid
(a normal plane, 4 floats instead of 1).  This would take too muh time for me
to code up to be worth it, and would probably run too much slower (10min for a
200x200 pixel picture now, ug.)  If somebody wants to give me a cray-3 to play
with, I'll write more accurate (and slower) code, until then... 8-)

Jon
bullerj@handel.cs.colostate.edu          ...!ncar!ccncsu!handel!bullerj
(These are my ideas (and code), nobody else SHOULD want these bugs)
I'm just trying to graduate.  Apple, Pixar, HP, etc. take note, I would like
your job offers, I have tired of the university life.


[NOTE: I have attached the Pascal code for the Turbulence functions.  The Noise
functions are in the next message in "C".  Sorry about the mixed languages, but
I haven't nicely translated these. -- EAH]

(* ---------- cut here ---------- cut here ---------- cut here ---------- *)

const
   MaxPts = 512;  { Must be 2^n}
   MPM1 = MaxPts - 1;
type
   PtsTyp = array[0..MaxPts] of Extended;
var
   Points: PtsTyp;


function Turb (Size: Integer;
               ScaleFactor: Extended;
               Loc: Vect;
               Pts: PtsTyp): Extended;
   var
      Scale, Result: Extended;
      Cur: Integer;
begin
   Result := Noise(Loc, Pts);
   Scale := 1.0;

   Cur := 1;
   while Cur < Size do begin
      Cur := BSL(Cur, 1);          {Cur := Cur * 2}
      Scale := Scale * ScaleFactor;
      Loc := Scale_Vect(2.0, Loc);
      Result := Result + Noise(Loc, Pts) * Scale;
   end;
   Turb := Result;
end;


function DTurb (Size: Integer;
                ScaleFactor: Extended;
                Loc: Vect;
                Pts: PtsTyp): Vect;
   var
      Result: Vect;
      Scale: Extended;
      Cur: Integer;
begin
   Result := DNoise(Loc, Pts);
   Scale := 1.0;

   Cur := 1;
   while Cur < Size do begin
      Cur := BSL(Cur, 1);
      Scale := Scale * ScaleFactor;
      Loc := Scale_Vect(2.0, Loc);
      Result := Add_Vect(Result, Scale_Vect(Scale, DNoise(Loc, Pts)));
   end;
   DTurb := Result;
end;

-----------------------------------------------------------------------------

And the C Version...

From: dirks@oak.cis.ohio-state.edu (william r dirks)
Subject: C Versions of Noise and DNoise Routines
Organization: Ohio State University Computer and Information Science


It was suggested to me that some of you would be interested in my 
translated-into-C-and-corrected versions of noise() and Dnoise().

So, here they are.  

Note that this is only noise and Dnoise.  The turbulence routines are
not included 'cause I haven't translated them yet.

Oh yeah, initnoise() fills the pts[] array with random numbers between
0 and 1.  Don't forget to initialize, or noise will always return 0.
(That's been experimentally verified! :-))

--Bill--

[Note that you should look over the rand() function if you use this stuff.
Some rand()'s need initialization (srand()), and some give numbers from 0
to 32767, and so should use this as a divisor. -- EAH]


---------cut-here------------------------------------cut-here--------
/*
**	Noise and Dnoise routines
*
*	Many thanks to Jon Buller of Colorado State University for these
*	routines.
*/


typedef struct vector {
   double x, y, z;
} Vector;


#define NUMPTS	512
#define P1	173
#define P2	263
#define P3	337
#define phi	0.6180339


static double pts[NUMPTS];


void initnoise()
{
   int i;
   
   for (i = 0; i < NUMPTS; ++i)
      pts[i] = rand() / 2.147483e9;
}


double noise(p)
Vector *p;
{
   int xi, yi, zi;
   int xa, xb, xc, ya, yb, yc, za, zb, zc;
   double xf, yf, zf;
   double x2, x1, x0, y2, y1, y0, z2, z1, z0;
   double p000, p100, p200, p010, p110, p210, p020, p120, p220;
   double p001, p101, p201, p011, p111, p211, p021, p121, p221;
   double p002, p102, p202, p012, p112, p212, p022, p122, p222;

   xi = floor(p->x);
   xa = floor(P1 * (xi * phi - floor(xi * phi)));
   xb = floor(P1 * ((xi + 1) * phi - floor((xi + 1) * phi)));
   xc = floor(P1 * ((xi + 2) * phi - floor((xi + 2) * phi)));

   yi = floor(p->y);
   ya = floor(P2 * (yi * phi - floor(yi * phi)));
   yb = floor(P2 * ((yi + 1) * phi - floor((yi + 1) * phi)));
   yc = floor(P2 * ((yi + 2) * phi - floor((yi + 2) * phi)));

   zi = floor(p->z);
   za = floor(P3 * (zi * phi - floor(zi * phi)));
   zb = floor(P3 * ((zi + 1) * phi - floor((zi + 1) * phi)));
   zc = floor(P3 * ((zi + 2) * phi - floor((zi + 2) * phi)));

   p000 = pts[xa + ya + za & NUMPTS - 1];
   p100 = pts[xb + ya + za & NUMPTS - 1];
   p200 = pts[xc + ya + za & NUMPTS - 1];
   p010 = pts[xa + yb + za & NUMPTS - 1];
   p110 = pts[xb + yb + za & NUMPTS - 1];
   p210 = pts[xc + yb + za & NUMPTS - 1];
   p020 = pts[xa + yc + za & NUMPTS - 1];
   p120 = pts[xb + yc + za & NUMPTS - 1];
   p220 = pts[xc + yc + za & NUMPTS - 1];
   p001 = pts[xa + ya + zb & NUMPTS - 1];
   p101 = pts[xb + ya + zb & NUMPTS - 1];
   p201 = pts[xc + ya + zb & NUMPTS - 1];
   p011 = pts[xa + yb + zb & NUMPTS - 1];
   p111 = pts[xb + yb + zb & NUMPTS - 1];
   p211 = pts[xc + yb + zb & NUMPTS - 1];
   p021 = pts[xa + yc + zb & NUMPTS - 1];
   p121 = pts[xb + yc + zb & NUMPTS - 1];
   p221 = pts[xc + yc + zb & NUMPTS - 1];
   p002 = pts[xa + ya + zc & NUMPTS - 1];
   p102 = pts[xb + ya + zc & NUMPTS - 1];
   p202 = pts[xc + ya + zc & NUMPTS - 1];
   p012 = pts[xa + yb + zc & NUMPTS - 1];
   p112 = pts[xb + yb + zc & NUMPTS - 1];
   p212 = pts[xc + yb + zc & NUMPTS - 1];
   p022 = pts[xa + yc + zc & NUMPTS - 1];
   p122 = pts[xb + yc + zc & NUMPTS - 1];
   p222 = pts[xc + yc + zc & NUMPTS - 1];

   xf = p->x - xi;
   x1 = xf * xf;
   x2 = 0.5 * x1;
   x1 = 0.5 + xf - x1;
   x0 = 0.5 - xf + x2;

   yf = p->y - yi;
   y1 = yf * yf;
   y2 = 0.5 * y1;
   y1 = 0.5 + yf - y1;
   y0 = 0.5 - yf + y2;

   zf = p->z - zi;
   z1 = zf * zf;
   z2 = 0.5 * z1;
   z1 = 0.5 + zf - z1;
   z0 = 0.5 - zf + z2;

   return   z0 * (y0 * (x0 * p000 + x1 * p100 + x2 * p200) +
                  y1 * (x0 * p010 + x1 * p110 + x2 * p210) +
                  y2 * (x0 * p020 + x1 * p120 + x2 * p220)) +
            z1 * (y0 * (x0 * p001 + x1 * p101 + x2 * p201) +
                  y1 * (x0 * p011 + x1 * p111 + x2 * p211) +
                  y2 * (x0 * p021 + x1 * p121 + x2 * p221)) +
            z2 * (y0 * (x0 * p002 + x1 * p102 + x2 * p202) +
                  y1 * (x0 * p012 + x1 * p112 + x2 * p212) +
                  y2 * (x0 * p022 + x1 * p122 + x2 * p222));
}



Vector Dnoise(p)
Vector *p;
{
   Vector v;
   int xi, yi, zi;
   int xa, xb, xc, ya, yb, yc, za, zb, zc;
   double xf, yf, zf;
   double x2, x1, x0, y2, y1, y0, z2, z1, z0;
   double xd2, xd1, xd0, yd2, yd1, yd0, zd2, zd1, zd0;
   double p000, p100, p200, p010, p110, p210, p020, p120, p220;
   double p001, p101, p201, p011, p111, p211, p021, p121, p221;
   double p002, p102, p202, p012, p112, p212, p022, p122, p222;

   xi = floor(p->x);
   xa = floor(P1 * (xi * phi - floor(xi * phi)));
   xb = floor(P1 * ((xi + 1) * phi - floor((xi + 1) * phi)));
   xc = floor(P1 * ((xi + 2) * phi - floor((xi + 2) * phi)));

   yi = floor(p->y);
   ya = floor(P2 * (yi * phi - floor(yi * phi)));
   yb = floor(P2 * ((yi + 1) * phi - floor((yi + 1) * phi)));
   yc = floor(P2 * ((yi + 2) * phi - floor((yi + 2) * phi)));

   zi = floor(p->z);
   za = floor(P3 * (zi * phi - floor(zi * phi)));
   zb = floor(P3 * ((zi + 1) * phi - floor((zi + 1) * phi)));
   zc = floor(P3 * ((zi + 2) * phi - floor((zi + 2) * phi)));

   p000 = pts[xa + ya + za & NUMPTS - 1];
   p100 = pts[xb + ya + za & NUMPTS - 1];
   p200 = pts[xc + ya + za & NUMPTS - 1];
   p010 = pts[xa + yb + za & NUMPTS - 1];
   p110 = pts[xb + yb + za & NUMPTS - 1];
   p210 = pts[xc + yb + za & NUMPTS - 1];
   p020 = pts[xa + yc + za & NUMPTS - 1];
   p120 = pts[xb + yc + za & NUMPTS - 1];
   p220 = pts[xc + yc + za & NUMPTS - 1];
   p001 = pts[xa + ya + zb & NUMPTS - 1];
   p101 = pts[xb + ya + zb & NUMPTS - 1];
   p201 = pts[xc + ya + zb & NUMPTS - 1];
   p011 = pts[xa + yb + zb & NUMPTS - 1];
   p111 = pts[xb + yb + zb & NUMPTS - 1];
   p211 = pts[xc + yb + zb & NUMPTS - 1];
   p021 = pts[xa + yc + zb & NUMPTS - 1];
   p121 = pts[xb + yc + zb & NUMPTS - 1];
   p221 = pts[xc + yc + zb & NUMPTS - 1];
   p002 = pts[xa + ya + zc & NUMPTS - 1];
   p102 = pts[xb + ya + zc & NUMPTS - 1];
   p202 = pts[xc + ya + zc & NUMPTS - 1];
   p012 = pts[xa + yb + zc & NUMPTS - 1];
   p112 = pts[xb + yb + zc & NUMPTS - 1];
   p212 = pts[xc + yb + zc & NUMPTS - 1];
   p022 = pts[xa + yc + zc & NUMPTS - 1];
   p122 = pts[xb + yc + zc & NUMPTS - 1];
   p222 = pts[xc + yc + zc & NUMPTS - 1];

   xf = p->x - xi;
   x1 = xf * xf;
   x2 = 0.5 * x1;
   x1 = 0.5 + xf - x1;
   x0 = 0.5 - xf + x2;
   xd2 = xf;
   xd1 = 1.0 - xf - xf;
   xd0 = xf - 1.0;

   yf = p->y - yi;
   y1 = yf * yf;
   y2 = 0.5 * y1;
   y1 = 0.5 + yf - y1;
   y0 = 0.5 - yf + y2;
   yd2 = yf;
   yd1 = 1.0 - yf - yf;
   yd0 = yf - 1.0;

   zf = p->z - zi;
   z1 = zf * zf;
   z2 = 0.5 * z1;
   z1 = 0.5 + zf - z1;
   z0 = 0.5 - zf + z2;
   zd2 = zf;
   zd1 = 1.0 - zf - zf;
   zd0 = zf - 1.0;

   v.x        = z0 * (y0 * (xd0 * p000 + xd1 * p100 + xd2 * p200) +
                      y1 * (xd0 * p010 + xd1 * p110 + xd2 * p210) +
                      y2 * (xd0 * p020 + xd1 * p120 + xd2 * p220)) +
                z1 * (y0 * (xd0 * p001 + xd1 * p101 + xd2 * p201) +
                      y1 * (xd0 * p011 + xd1 * p111 + xd2 * p211) +
                      y2 * (xd0 * p021 + xd1 * p121 + xd2 * p221)) +
                z2 * (y0 * (xd0 * p002 + xd1 * p102 + xd2 * p202) +
                      y1 * (xd0 * p012 + xd1 * p112 + xd2 * p212) +
                      y2 * (xd0 * p022 + xd1 * p122 + xd2 * p222));
                                  
   v.y        = z0 * (yd0 * (x0 * p000 + x1 * p100 + x2 * p200) +
                      yd1 * (x0 * p010 + x1 * p110 + x2 * p210) +
                      yd2 * (x0 * p020 + x1 * p120 + x2 * p220)) +
                z1 * (yd0 * (x0 * p001 + x1 * p101 + x2 * p201) +
                      yd1 * (x0 * p011 + x1 * p111 + x2 * p211) +
                      yd2 * (x0 * p021 + x1 * p121 + x2 * p221)) +
                z2 * (yd0 * (x0 * p002 + x1 * p102 + x2 * p202) +
                      yd1 * (x0 * p012 + x1 * p112 + x2 * p212) +
                      yd2 * (x0 * p022 + x1 * p122 + x2 * p222));
                                  
   v.z        = zd0 * (y0 * (x0 * p000 + x1 * p100 + x2 * p200) +
                       y1 * (x0 * p010 + x1 * p110 + x2 * p210) +
                       y2 * (x0 * p020 + x1 * p120 + x2 * p220)) +
                zd1 * (y0 * (x0 * p001 + x1 * p101 + x2 * p201) +
                       y1 * (x0 * p011 + x1 * p111 + x2 * p211) +
                       y2 * (x0 * p021 + x1 * p121 + x2 * p221)) +
                zd2 * (y0 * (x0 * p002 + x1 * p102 + x2 * p202) +
                       y1 * (x0 * p012 + x1 * p112 + x2 * p212) +
                       y2 * (x0 * p022 + x1 * p122 + x2 * p222));
   return v;
}

-----------------------------------------------------------------------------

END OF RTNEWS