src/rt/noise3.c

#ifndef lint
static const char       RCSid[] = "$Id: noise3.c,v 2.12 2011/02/22 16:45:12 greg Exp $";
#endif
/*
 *  noise3.c - noise functions for random textures.
 *
 *     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"

#include  "ray.h"
#include  "func.h"

static char  noise_name[4][8] = {"noise3x", "noise3y", "noise3z", "noise3"};
static char  fnoise_name[] = "fnoise3";

#define  EPSILON        .0005           /* error allowed in fractal */

#define  frand3(x,y,z)  frand(17*(x)+23*(y)+29*(z))

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 */
        char  *nam
)
{
        int  i;
        double  x[3];
                                        /* get point */
        x[0] = argument(1);
        x[1] = argument(2);
        x[2] = argument(3);
                                        /* make appropriate call */
        if (nam == fnoise_name)
                return(fnoise3(x));
        i = 4;
        while (i--)
                if (nam == noise_name[i])
                        return(noise3(x,i));
        eputs(nam);
        eputs(": called l_noise3!\n");
        quit(1);
        return 1; /* pro forma return */
}


void
setnoisefuncs(void)                     /* add noise functions to library */
{
        int  i;

        funset(fnoise_name, 3, ':', l_noise3);
        i = 4;
        while (i--)
                funset(noise_name[i], 3, ':', l_noise3);
}


static double
frand(                  /* get random number from seed */
        long  s
)
{
        s = s<<13 ^ s;
        return(1.0-((s*(s*s*15731+789221)+1376312589)&0x7fffffff)/1073741824.0);
}


static double
fnoise3(                        /* compute fractal noise function */
        double  x[3]
)
{
        long  t[3], v[3], beg[3];
        double  fval[8], fc;
        int  branch;
        long  s;
        int  i, j;
                                                /* get starting cube */
        s = (long)(1.0/EPSILON);
        for (i = 0; i < 3; i++) {
                t[i] = s*x[i];
                beg[i] = s*floor(x[i]);
        }
        for (j = 0; j < 8; j++) {
                for (i = 0; i < 3; i++) {
                        v[i] = beg[i];
                        if (j & 1<<i)
                                v[i] += s;
                }
                fval[j] = frand3(v[0],v[1],v[2]);
        }
                                                /* compute fractal */
        for ( ; ; ) {
                fc = 0.0;
                for (j = 0; j < 8; j++)
                        fc += fval[j];
                fc *= 0.125;
                if ((s >>= 1) == 0)
                        return(fc);             /* close enough */
                branch = 0;
                for (i = 0; i < 3; i++) {       /* do center */
                        v[i] = beg[i] + s;
                        if (t[i] > v[i]) {
                                branch |= 1<<i;
                        }
                }
                fc += s*EPSILON*frand3(v[0],v[1],v[2]);
                fval[~branch & 7] = fc;
                for (i = 0; i < 3; i++) {       /* do faces */
                        if (branch & 1<<i)
                                v[i] += s;
                        else
                                v[i] -= s;
                        fc = 0.0;
                        for (j = 0; j < 8; j++)
                                if (~(j^branch) & 1<<i)
                                        fc += fval[j];
                        fc = 0.25*fc + s*EPSILON*frand3(v[0],v[1],v[2]);
                        fval[~(branch^1<<i) & 7] = fc;
                        v[i] = beg[i] + s;
                }
                for (i = 0; i < 3; i++) {       /* do edges */
                        if ((j = i+1) == 3) j = 0;
                        if (branch & 1<<j)
                                v[j] += s;
                        else
                                v[j] -= s;
                        if (++j == 3) j = 0;
                        if (branch & 1<<j)
                                v[j] += s;
                        else
                                v[j] -= s;
                        fc = fval[branch & ~(1<<i)];
                        fc += fval[branch | 1<<i];
                        fc = 0.5*fc + s*EPSILON*frand3(v[0],v[1],v[2]);
                        fval[branch^1<<i] = fc;
                        if ((j = i+1) == 3) j = 0;
                        v[j] = beg[j] + s;
                        if (++j == 3) j = 0;
                        v[j] = beg[j] + s;
                }
                for (i = 0; i < 3; i++)         /* new cube */
                        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))));
}
Revision:	2.13
Committed:	Tue Oct 8 18:59:44 2013 UTC (12 years ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R2, rad5R0, rad5R1, rad4R2, rad4R2P1, rad4R2P2
Changes since 2.12:	+128 -107 lines
Log Message:	Implemented Perlin's improved noise function
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.13	static const char RCSid[] = "$Id: noise3.c,v 2.12 2011/02/22 16:45:12 greg Exp $";
3	greg	1.1	#endif
4			/*
5			* noise3.c - noise functions for random textures.
6			*
7	greg	2.13	* 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	greg	2.6	*/
11
12	greg	2.7	#include "copyright.h"
13	greg	1.1
14	greg	2.12	#include "ray.h"
15	schorsch	2.9	#include "func.h"
16	greg	1.1
17	greg	2.5	static char noise_name[4][8] = {"noise3x", "noise3y", "noise3z", "noise3"};
18	greg	1.5	static char fnoise_name[] = "fnoise3";
19	greg	1.1
20	greg	2.13	#define EPSILON .0005 /* error allowed in fractal */
21	greg	1.1
22	greg	1.3	#define frand3(x,y,z) frand(17(x)+23(y)+29*(z))
23	greg	1.1
24	greg	2.13	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	schorsch	2.9
31	greg	1.1
32	greg	1.5	static double
33	schorsch	2.9	l_noise3( /* compute a noise function */
34	greg	2.13	char *nam
35	schorsch	2.9	)
36	greg	1.1	{
37	greg	2.13	int i;
38	greg	1.5	double x[3];
39			/* get point */
40			x[0] = argument(1);
41			x[1] = argument(2);
42			x[2] = argument(3);
43			/* make appropriate call */
44			if (nam == fnoise_name)
45			return(fnoise3(x));
46			i = 4;
47			while (i--)
48			if (nam == noise_name[i])
49	greg	2.13	return(noise3(x,i));
50	greg	1.6	eputs(nam);
51			eputs(": called l_noise3!\n");
52	greg	1.5	quit(1);
53	schorsch	2.9	return 1; /* pro forma return */
54	greg	1.1	}
55
56
57	greg	2.13	void
58	schorsch	2.9	setnoisefuncs(void) /* add noise functions to library */
59	greg	1.1	{
60	greg	2.13	int i;
61	greg	1.1
62	greg	1.5	funset(fnoise_name, 3, ':', l_noise3);
63			i = 4;
64			while (i--)
65			funset(noise_name[i], 3, ':', l_noise3);
66	greg	1.1	}
67
68
69	schorsch	2.9	static double
70			frand( /* get random number from seed */
71	greg	2.13	long s
72	schorsch	2.9	)
73	greg	1.1	{
74			s = s<<13 ^ s;
75			return(1.0-((s(ss*15731+789221)+1376312589)&0x7fffffff)/1073741824.0);
76			}
77
78
79	schorsch	2.9	static double
80			fnoise3( /* compute fractal noise function */
81	greg	2.13	double x[3]
82	schorsch	2.9	)
83	greg	1.1	{
84	greg	1.4	long t[3], v[3], beg[3];
85	greg	1.3	double fval[8], fc;
86			int branch;
87	greg	2.13	long s;
88			int i, j;
89	greg	1.1	/* get starting cube */
90	greg	1.3	s = (long)(1.0/EPSILON);
91			for (i = 0; i < 3; i++) {
92	greg	2.13	t[i] = s*x[i];
93			beg[i] = s*floor(x[i]);
94	greg	1.3	}
95	greg	1.1	for (j = 0; j < 8; j++) {
96			for (i = 0; i < 3; i++) {
97			v[i] = beg[i];
98			if (j & 1<<i)
99	greg	1.3	v[i] += s;
100	greg	1.1	}
101			fval[j] = frand3(v[0],v[1],v[2]);
102			}
103			/* compute fractal */
104			for ( ; ; ) {
105	greg	1.4	fc = 0.0;
106			for (j = 0; j < 8; j++)
107			fc += fval[j];
108			fc *= 0.125;
109			if ((s >>= 1) == 0)
110			return(fc); /* close enough */
111	greg	1.1	branch = 0;
112			for (i = 0; i < 3; i++) { /* do center */
113			v[i] = beg[i] + s;
114	greg	1.3	if (t[i] > v[i]) {
115	greg	1.1	branch \|= 1<<i;
116	greg	1.3	}
117	greg	1.1	}
118	greg	1.3	fc += sEPSILONfrand3(v[0],v[1],v[2]);
119	greg	1.1	fval[~branch & 7] = fc;
120			for (i = 0; i < 3; i++) { /* do faces */
121			if (branch & 1<<i)
122			v[i] += s;
123			else
124			v[i] -= s;
125			fc = 0.0;
126			for (j = 0; j < 8; j++)
127			if (~(j^branch) & 1<<i)
128			fc += fval[j];
129	greg	1.3	fc = 0.25fc + sEPSILON*frand3(v[0],v[1],v[2]);
130	greg	1.1	fval[~(branch^1<<i) & 7] = fc;
131			v[i] = beg[i] + s;
132			}
133			for (i = 0; i < 3; i++) { /* do edges */
134	greg	2.10	if ((j = i+1) == 3) j = 0;
135	greg	1.1	if (branch & 1<<j)
136			v[j] += s;
137			else
138			v[j] -= s;
139	greg	2.10	if (++j == 3) j = 0;
140	greg	1.1	if (branch & 1<<j)
141			v[j] += s;
142			else
143			v[j] -= s;
144			fc = fval[branch & ~(1<<i)];
145			fc += fval[branch \| 1<<i];
146	greg	1.3	fc = 0.5fc + sEPSILON*frand3(v[0],v[1],v[2]);
147	greg	1.1	fval[branch^1<<i] = fc;
148	greg	2.10	if ((j = i+1) == 3) j = 0;
149	greg	1.1	v[j] = beg[j] + s;
150	greg	2.10	if (++j == 3) j = 0;
151	greg	1.1	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	greg	2.13
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 && x[0]==xnew[0] && (x[1]==xnew[1]) & (x[2]==xnew[2])) {
170			if (!(gotV>>i & 1)) {
171			f[i] = noise3partial(f[3], x, i);
172			gotV \|= 1<<i;
173			}
174			return(f[i]);
175			}
176			gotV = 0x8;
177			return(f[3] = perlin_noise(x[0]=xnew[0], x[1]=xnew[1], x[2]=xnew[2]));
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			}