--- ray/src/rt/noise3.c	2011/02/22 16:45:12	2.12
+++ ray/src/rt/noise3.c	2013/10/08 18:59:44	2.13
@@ -1,11 +1,12 @@
 #ifndef lint
-static const char	RCSid[] = "$Id: noise3.c,v 2.12 2011/02/22 16:45:12 greg Exp $";
+static const char	RCSid[] = "$Id: noise3.c,v 2.13 2013/10/08 18:59:44 greg Exp $";
 #endif
 /*
  *  noise3.c - noise functions for random textures.
  *
- *     Credit for the smooth algorithm goes to Ken Perlin.
- *     (ref. SIGGRAPH Vol 19, No 3, pp 287-96)
+ *     Credit for the smooth algorithm goes to Ken Perlin, and code
+ *	translation/implementation to Rahul Narain.
+ *	(ref. Improving Noise, Computer Graphics; Vol. 35 No. 3., 2002)
  */
 
 #include "copyright.h"
@@ -13,51 +14,27 @@ static const char	RCSid[] = "$Id: noise3.c,v 2.12 2011
 #include  "ray.h"
 #include  "func.h"
 
-#define  A		0
-#define  B		1
-#define  C		2
-#define  D		3
-
-#define  rand3a(x,y,z)	frand(67*(x)+59*(y)+71*(z))
-#define  rand3b(x,y,z)	frand(73*(x)+79*(y)+83*(z))
-#define  rand3c(x,y,z)	frand(89*(x)+97*(y)+101*(z))
-#define  rand3d(x,y,z)	frand(103*(x)+107*(y)+109*(z))
-
-#define  hpoly1(t)	((2.0*t-3.0)*t*t+1.0)
-#define  hpoly2(t)	(-2.0*t+3.0)*t*t
-#define  hpoly3(t)	((t-2.0)*t+1.0)*t
-#define  hpoly4(t)	(t-1.0)*t*t
-
-#define  hermite(p0,p1,r0,r1,t)  (	p0*hpoly1(t) + \
-					p1*hpoly2(t) + \
-					r0*hpoly3(t) + \
-					r1*hpoly4(t) )
-
 static char  noise_name[4][8] = {"noise3x", "noise3y", "noise3z", "noise3"};
 static char  fnoise_name[] = "fnoise3";
-static char  hermite_name[] = "hermite";
 
-static long  xlim[3][2];
-static double  xarg[3];
+#define  EPSILON	.0005		/* error allowed in fractal */
 
-#define  EPSILON	.001		/* error allowed in fractal */
-
 #define  frand3(x,y,z)	frand(17*(x)+23*(y)+29*(z))
 
-static double l_noise3(char  *nam);
-static double l_hermite(char *nm);
-static double * noise3(double  xnew[3]);
-static void interpolate(double  f[4], int  i, int  n);
-static double frand(long  s);
-static double fnoise3(double  p[3]);
+static double l_noise3(char *nam);
+static double noise3(double xnew[3], int i);
+static double noise3partial(double f3, double x[3], int i);
+static double perlin_noise (double x, double y, double z);
+static double frand(long s);
+static double fnoise3(double x[3]);
 
 
 static double
 l_noise3(			/* compute a noise function */
-	register char  *nam
+	char  *nam
 )
 {
-	register int  i;
+	int  i;
 	double  x[3];
 					/* get point */
 	x[0] = argument(1);
@@ -69,7 +46,7 @@ l_noise3(			/* compute a noise function */
 	i = 4;
 	while (i--)
 		if (nam == noise_name[i])
-			return(noise3(x)[i]);
+			return(noise3(x,i));
 	eputs(nam);
 	eputs(": called l_noise3!\n");
 	quit(1);
@@ -77,23 +54,11 @@ l_noise3(			/* compute a noise function */
 }
 
 
-static double
-l_hermite(char *nm)		/* library call for hermite interpolation */
-{
-	double  t;
-	
-	t = argument(5);
-	return( hermite(argument(1), argument(2),
-			argument(3), argument(4), t) );
-}
-
-
-extern void
+void
 setnoisefuncs(void)			/* add noise functions to library */
 {
-	register int  i;
+	int  i;
 
-	funset(hermite_name, 5, ':', l_hermite);
 	funset(fnoise_name, 3, ':', l_noise3);
 	i = 4;
 	while (i--)
@@ -101,59 +66,9 @@ setnoisefuncs(void)			/* add noise functions to librar
 }
 
 
-static double *
-noise3(			/* compute the noise function */
-	register double  xnew[3]
-)
-{
-	static double  x[3] = {-100000.0, -100000.0, -100000.0};
-	static double  f[4];
-
-	if (x[0]==xnew[0] && x[1]==xnew[1] && x[2]==xnew[2])
-		return(f);
-	x[0] = xnew[0]; x[1] = xnew[1]; x[2] = xnew[2];
-	xlim[0][0] = floor(x[0]); xlim[0][1] = xlim[0][0] + 1;
-	xlim[1][0] = floor(x[1]); xlim[1][1] = xlim[1][0] + 1;
-	xlim[2][0] = floor(x[2]); xlim[2][1] = xlim[2][0] + 1;
-	xarg[0] = x[0] - xlim[0][0];
-	xarg[1] = x[1] - xlim[1][0];
-	xarg[2] = x[2] - xlim[2][0];
-	interpolate(f, 0, 3);
-	return(f);
-}
-
-
-static void
-interpolate(
-	double  f[4],
-	register int  i,
-	register int  n
-)
-{
-	double  f0[4], f1[4], hp1, hp2;
-
-	if (n == 0) {
-		f[A] = rand3a(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
-		f[B] = rand3b(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
-		f[C] = rand3c(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
-		f[D] = rand3d(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
-	} else {
-		n--;
-		interpolate(f0, i, n);
-		interpolate(f1, i | 1<<n, n);
-		hp1 = hpoly1(xarg[n]); hp2 = hpoly2(xarg[n]);
-		f[A] = f0[A]*hp1 + f1[A]*hp2;
-		f[B] = f0[B]*hp1 + f1[B]*hp2;
-		f[C] = f0[C]*hp1 + f1[C]*hp2;
-		f[D] = f0[D]*hp1 + f1[D]*hp2 +
-				f0[n]*hpoly3(xarg[n]) + f1[n]*hpoly4(xarg[n]);
-	}
-}
-
-
 static double
 frand(			/* get random number from seed */
-	register long  s
+	long  s
 )
 {
 	s = s<<13 ^ s;
@@ -163,19 +78,19 @@ frand(			/* get random number from seed */
 
 static double
 fnoise3(			/* compute fractal noise function */
-	double  p[3]
+	double  x[3]
 )
 {
 	long  t[3], v[3], beg[3];
 	double  fval[8], fc;
 	int  branch;
-	register long  s;
-	register int  i, j;
+	long  s;
+	int  i, j;
 						/* get starting cube */
 	s = (long)(1.0/EPSILON);
 	for (i = 0; i < 3; i++) {
-		t[i] = s*p[i];
-		beg[i] = s*floor(p[i]);
+		t[i] = s*x[i];
+		beg[i] = s*floor(x[i]);
 	}
 	for (j = 0; j < 8; j++) {
 		for (i = 0; i < 3; i++) {
@@ -239,4 +154,110 @@ fnoise3(			/* compute fractal noise function */
 			if (branch & 1<<i)
 				beg[i] += s;
 	}
+}
+
+
+static double
+noise3(			/* compute the revised Perlin noise function */
+	double  xnew[3], int i
+)
+{
+	static int     gotV;
+	static double  x[3];
+	static double  f[4];
+
+	if (gotV && x[0]==xnew[0] && (x[1]==xnew[1]) & (x[2]==xnew[2])) {
+		if (!(gotV>>i & 1)) {
+			f[i] = noise3partial(f[3], x, i);
+			gotV |= 1<<i;
+		}
+		return(f[i]);
+	}
+	gotV = 0x8;
+	return(f[3] = perlin_noise(x[0]=xnew[0], x[1]=xnew[1], x[2]=xnew[2]));
+}
+
+static double
+noise3partial(		/* compute partial derivative for ith coordinate */
+	double f3, double x[3], int i
+)
+{
+	double	fc;
+
+	switch (i) {
+	case 0:
+		fc = perlin_noise(x[0]-EPSILON, x[1], x[2]);
+		break;
+	case 1:
+		fc = perlin_noise(x[0], x[1]-EPSILON, x[2]);
+		break;
+	case 2:
+		fc = perlin_noise(x[0], x[1], x[2]-EPSILON);
+		break;
+	default:
+		return(.0);
+	}
+	return((f3 - fc)/EPSILON);
+}
+
+#define fade(t) ((t)*(t)*(t)*((t)*((t)*6. - 15.) + 10.))
+
+static double lerp(double t, double a, double b) {return a + t*(b - a);}
+
+static double
+grad(int hash, double x, double y, double z)
+{
+    int h = hash & 15;                      // CONVERT LO 4 BITS OF HASH CODE
+    double u = h<8 ? x : y,                 // INTO 12 GRADIENT DIRECTIONS.
+                 v = h<4 ? y : h==12|h==14 ? x : z;
+    return (!(h&1) ? u : -u) + (!(h&2) ? v : -v);
+}
+
+static const int permutation[256] = {151,160,137,91,90,15,
+    131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
+    190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
+    88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
+    77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
+    102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
+    135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
+    5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
+    223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
+    129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
+    251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
+    49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
+    138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180
+};
+
+#define	p(i)	permutation[(i)&0xff]
+
+static double
+perlin_noise(double x, double y, double z) 
+{
+    int		X, Y, Z;
+    double	u, v, w;
+    int		A, AA, AB, B, BA, BB;
+
+    X = (int)x-(x<0);                    // FIND UNIT CUBE THAT
+    Y = (int)y-(y<0);                    // CONTAINS POINT.
+    Z = (int)z-(z<0);
+    x -= (double)X;                          // FIND RELATIVE X,Y,Z
+    y -= (double)Y;                          // OF POINT IN CUBE.
+    z -= (double)Z;
+    X &= 0xff; Y &= 0xff; Z &= 0xff;
+
+    u = fade(x);                                  // COMPUTE FADE CURVES
+    v = fade(y);                                  // FOR EACH OF X,Y,Z.
+    w = fade(z);
+
+    A = p(X  )+Y; AA = p(A)+Z; AB = p(A+1)+Z;          // HASH COORDINATES OF
+    B = p(X+1)+Y; BA = p(B)+Z; BB = p(B+1)+Z;          // THE 8 CUBE CORNERS,
+
+    return lerp(w, lerp(v, lerp(u, grad(p(AA  ), x  , y  , z  ),  // AND ADD
+                                   grad(p(BA  ), x-1, y  , z  )), // BLENDED
+                           lerp(u, grad(p(AB  ), x  , y-1, z  ),  // RESULTS
+                                   grad(p(BB  ), x-1, y-1, z  ))),// FROM  8
+                   lerp(v, lerp(u, grad(p(AA+1), x  , y  , z-1),  // CORNERS
+                                   grad(p(BA+1), x-1, y  , z-1)), // OF CUBE
+                           lerp(u, grad(p(AB+1), x  , y-1, z-1),
+                                   grad(p(BB+1), x-1, y-1, z-1))));
 }