Path: utzoo!utgpu!jarvis.csri.toronto.edu!mailrus!ncar!boulder!ccncsu!handel.colostate.edu!bullerj
From: bullerj@handel.colostate.edu (Jon Buller)
Newsgroups: comp.graphics
Subject: Re: Pixar's noise function
Summary: Heres my code...
Keywords: random numbers, Byte, elephants
Message-ID: <1540@ccncsu.ColoState.EDU>
Date: 31 Mar 89 07:02:06 GMT
References: <2553@ssc-vax.UUCP> <3599@pixar.UUCP>
Sender: news@ccncsu.ColoState.EDU
Reply-To: bullerj@handel.colostate.edu.UUCP (Jon Buller)
Organization: Colorado State University, Ft. Collins CO 80523
Lines: 305

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 belive I write
better Pascal than C...  One final note, this was Lightspeed Pascal 2.0 for
the Macintosh, but things have been reformatted slighly to get it on the net.
I hope this is what you all wanted.  Enjoy,

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.

(* ---------- 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 Noise (Loc: Vect;
                Pts: PtsTyp): Extended;
   const
      P1 = 173;
      P2 = 263;   {Reasonable numbers to (attempt to) avoid diagonal repetition}
      P3 = 337;
      phi = 0.618033988749894842;   {phi = (Sqrt(5)-1)/2}
   var
      xi, yi, zi: Integer;
      xa, xb, xc, ya, yb, yc, za, zb, zc: Integer;
      xf, yf, zf: Extended;
      x2, x1, x0, y2, y1, y0, z2, z1, z0: Extended;
      p000, p100, p200, p010, p110, p210, p020, p120, p220: Extended;
      p001, p101, p201, p011, p111, p211, p021, p121, p221: Extended;
      p002, p102, p202, p012, p112, p212, p022, p122, p222: Extended;
begin
   xi := Trunc(Loc[1]);
   xa := Trunc(P1 * (xi * phi - Trunc(xi * phi)));
   xb := Trunc(P1 * ((xi + 1) * phi - Trunc((xi + 1) * phi)));
   xc := Trunc(P1 * ((xi + 2) * phi - Trunc((xi + 2) * phi)));

   yi := Trunc(Loc[2]);
   ya := Trunc(P2 * (yi * phi - Trunc(yi * phi)));
   yb := Trunc(P2 * ((yi + 1) * phi - Trunc((yi + 1) * phi)));
   yc := Trunc(P2 * ((yi + 2) * phi - Trunc((yi + 2) * phi)));

   zi := Trunc(Loc[3]);
   za := Trunc(P3 * (zi * phi - Trunc(zi * phi)));
   zb := Trunc(P3 * ((zi + 1) * phi - Trunc((zi + 1) * phi)));
   zc := Trunc(P3 * ((zi + 2) * phi - Trunc((zi + 2) * phi)));

   p000 := Pts[BitAnd(xa + ya + za, MPM1)];{p000 := Pts[(xa+ya+za) mod MaxPts]}
   p100 := Pts[BitAnd(xb + ya + za, MPM1)];
   p200 := Pts[BitAnd(xc + ya + za, MPM1)];
   p010 := Pts[BitAnd(xa + yb + za, MPM1)];
   p110 := Pts[BitAnd(xb + yb + za, MPM1)];
   p210 := Pts[BitAnd(xc + yb + za, MPM1)];
   p020 := Pts[BitAnd(xa + yc + za, MPM1)];
   p120 := Pts[BitAnd(xb + yc + za, MPM1)];
   p220 := Pts[BitAnd(xc + yc + za, MPM1)];
   p001 := Pts[BitAnd(xa + ya + zb, MPM1)];
   p101 := Pts[BitAnd(xb + ya + zb, MPM1)];
   p201 := Pts[BitAnd(xc + ya + zb, MPM1)];
   p011 := Pts[BitAnd(xa + yb + zb, MPM1)];
   p111 := Pts[BitAnd(xb + yb + zb, MPM1)];
   p211 := Pts[BitAnd(xc + yb + zb, MPM1)];
   p021 := Pts[BitAnd(xa + yc + zb, MPM1)];
   p121 := Pts[BitAnd(xb + yc + zb, MPM1)];
   p221 := Pts[BitAnd(xc + yc + zb, MPM1)];
   p002 := Pts[BitAnd(xa + ya + zc, MPM1)];
   p102 := Pts[BitAnd(xb + ya + zc, MPM1)];
   p202 := Pts[BitAnd(xc + ya + zc, MPM1)];
   p012 := Pts[BitAnd(xa + yb + zc, MPM1)];
   p112 := Pts[BitAnd(xb + yb + zc, MPM1)];
   p212 := Pts[BitAnd(xc + yb + zc, MPM1)];
   p022 := Pts[BitAnd(xa + yc + zc, MPM1)];
   p122 := Pts[BitAnd(xb + yc + zc, MPM1)];
   p222 := Pts[BitAnd(xc + yc + zc, MPM1)];

   xf := Loc[1] - xi;
   x1 := xf * xf;
   x2 := 0.5 * x1;
   x1 := 0.5 + xf - x1;
   x0 := 0.5 - xf + x2;

   yf := Loc[2] - yi;
   y1 := yf * yf;
   y2 := 0.5 * y1;
   y1 := 0.5 + yf - y1;
   y0 := 0.5 - yf + y2;

   zf := Loc[3] - zi;
   z1 := zf * zf;
   z2 := 0.5 * z1;
   z1 := 0.5 + zf - z1;
   z0 := 0.5 - zf + z2;

   Noise := 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));
end;


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 DNoise (Loc: Vect;
                 Pts: PtsTyp): Vect;
   const
      P1 = 173;
      P2 = 263;
      P3 = 337;
      phi = 0.618033988749894842;
   var
      xi, yi, zi: Integer;
      xa, xb, xc, ya, yb, yc, za, zb, zc: Integer;
      xf, yf, zf: Extended;
      x2, x1, x0, y2, y1, y0, z2, z1, z0: Extended;
      xd2, xd1, xd0, yd2, yd1, yd0, zd2, zd1, zd0: Extended;
      p000, p100, p200, p010, p110, p210, p020, p120, p220: Extended;
      p001, p101, p201, p011, p111, p211, p021, p121, p221: Extended;
      p002, p102, p202, p012, p112, p212, p022, p122, p222: Extended;
begin
   xi := Trunc(Loc[1]);
   xa := Trunc(P1 * (xi * phi - Trunc(xi * phi)));
   xb := Trunc(P1 * ((xi + 1) * phi - Trunc((xi + 1) * phi)));
   xc := Trunc(P1 * ((xi + 2) * phi - Trunc((xi + 2) * phi)));

   yi := Trunc(Loc[2]);
   ya := Trunc(P2 * (yi * phi - Trunc(yi * phi)));
   yb := Trunc(P2 * ((yi + 1) * phi - Trunc((yi + 1) * phi)));
   yc := Trunc(P2 * ((yi + 2) * phi - Trunc((yi + 2) * phi)));

   zi := Trunc(Loc[3]);
   za := Trunc(P3 * (zi * phi - Trunc(zi * phi)));
   zb := Trunc(P3 * ((zi + 1) * phi - Trunc((zi + 1) * phi)));
   zc := Trunc(P3 * ((zi + 2) * phi - Trunc((zi + 2) * phi)));

   p000 := Pts[BitAnd(xa + ya + za, MPM1)];
   p100 := Pts[BitAnd(xb + ya + za, MPM1)];
   p200 := Pts[BitAnd(xc + ya + za, MPM1)];
   p010 := Pts[BitAnd(xa + yb + za, MPM1)];
   p110 := Pts[BitAnd(xb + yb + za, MPM1)];
   p210 := Pts[BitAnd(xc + yb + za, MPM1)];
   p020 := Pts[BitAnd(xa + yc + za, MPM1)];
   p120 := Pts[BitAnd(xb + yc + za, MPM1)];
   p220 := Pts[BitAnd(xc + yc + za, MPM1)];
   p001 := Pts[BitAnd(xa + ya + zb, MPM1)];
   p101 := Pts[BitAnd(xb + ya + zb, MPM1)];
   p201 := Pts[BitAnd(xc + ya + zb, MPM1)];
   p011 := Pts[BitAnd(xa + yb + zb, MPM1)];
   p111 := Pts[BitAnd(xb + yb + zb, MPM1)];
   p211 := Pts[BitAnd(xc + yb + zb, MPM1)];
   p021 := Pts[BitAnd(xa + yc + zb, MPM1)];
   p121 := Pts[BitAnd(xb + yc + zb, MPM1)];
   p221 := Pts[BitAnd(xc + yc + zb, MPM1)];
   p002 := Pts[BitAnd(xa + ya + zc, MPM1)];
   p102 := Pts[BitAnd(xb + ya + zc, MPM1)];
   p202 := Pts[BitAnd(xc + ya + zc, MPM1)];
   p012 := Pts[BitAnd(xa + yb + zc, MPM1)];
   p112 := Pts[BitAnd(xb + yb + zc, MPM1)];
   p212 := Pts[BitAnd(xc + yb + zc, MPM1)];
   p022 := Pts[BitAnd(xa + yc + zc, MPM1)];
   p122 := Pts[BitAnd(xb + yc + zc, MPM1)];
   p222 := Pts[BitAnd(xc + yc + zc, MPM1)];

   xf := Loc[1] - 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 := Loc[2] - 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 := Loc[3] - 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;

   DNoise[1] := 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));
                                  
   DNoise[2] := 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));
                                  
   DNoise[3] := 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));
   DNoise[4] := 0.0;
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;
-------------------------------------------------------------------------------
Jon Buller                                      FROM fortune IMPORT quote;
..!ccncsu!handel!bullerj                       FROM lawyers IMPORT disclaimer;