ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/Development/ray/src/rt/noise3.c

Comparing ray/src/rt/noise3.c (file contents):
Revision 2.8 by greg, Mon Aug 4 22:37:53 2003 UTC vs.
Revision 2.14 by greg, Wed Nov 21 20:23:15 2018 UTC

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

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)