src/rt/noise3.c

/* Copyright (c) 1988 Regents of the University of California */

#ifndef lint
static char SCCSid[] = "$SunId$ LBL";
#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)
 *
 *     4/15/86
 *     5/19/88  Added fractal noise function
 */


#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  hermite(p0,p1,r0,r1,t)  (      p0*((2.0*t-3.0)*t*t+1.0) + \
                                        p1*(-2.0*t+3.0)*t*t + \
                                        r0*((t-2.0)*t+1.0)*t + \
                                        r1*(t-1.0)*t*t )

static char  *noise_name[4] = {"noise3a", "noise3b", "noise3c", "noise3"};
static char  fnoise_name[] = "fnoise3";
static char  hermite_name[] = "hermite";

double  *noise3(), fnoise3(), argument(), frand();

static long  xlim[3][2];
static double  xarg[3];

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

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


static double
l_noise3(nam)                   /* compute a noise function */
register char  *nam;
{
        register 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("Bad call to l_noise3()!\n");
        quit(1);
}


double
l_hermite()                     /* library call for hermite interpolation */
{
        double  t;
        
        t = argument(5);
        return( hermite(argument(1), argument(2),
                        argument(3), argument(4), t) );
}


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

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


double *
noise3(xnew)                    /* compute the noise function */
register double  xnew[3];
{
        extern double  floor();
        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
interpolate(f, i, n)
double  f[4];
register int  i, n;
{
        double  f0[4], f1[4];

        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);
                f[A] = (1.0-xarg[n])*f0[A] + xarg[n]*f1[A];
                f[B] = (1.0-xarg[n])*f0[B] + xarg[n]*f1[B];
                f[C] = (1.0-xarg[n])*f0[C] + xarg[n]*f1[C];
                f[D] = hermite(f0[D], f1[D], f0[n], f1[n], xarg[n]);
        }
}


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


double
fnoise3(p)                      /* compute fractal noise function */
double  p[3];
{
        double  floor();
        long  t[3], v[3], beg[3];
        double  fval[8], fc;
        int  branch;
        register long  s;
        register 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]);
        }
        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 */
                        j = (i+1)%3;
                        if (branch & 1<<j)
                                v[j] += s;
                        else
                                v[j] -= s;
                        j = (i+2)%3;
                        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;
                        j = (i+1)%3;
                        v[j] = beg[j] + s;
                        j = (i+2)%3;
                        v[j] = beg[j] + s;
                }
                for (i = 0; i < 3; i++)         /* new cube */
                        if (branch & 1<<i)
                                beg[i] += s;
        }
}
Revision:	1.5
Committed:	Fri May 24 14:23:33 1991 UTC (33 years, 5 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 1.4:	+36 -55 lines
Log Message:	changed noise3 initialization
#	Content
1	/* Copyright (c) 1988 Regents of the University of California */
2
3	#ifndef lint
4	static char SCCSid[] = "$SunId$ LBL";
5	#endif
6
7	/*
8	* noise3.c - noise functions for random textures.
9	*
10	* Credit for the smooth algorithm goes to Ken Perlin.
11	* (ref. SIGGRAPH Vol 19, No 3, pp 287-96)
12	*
13	* 4/15/86
14	* 5/19/88 Added fractal noise function
15	*/
16
17
18	#define A 0
19	#define B 1
20	#define C 2
21	#define D 3
22
23	#define rand3a(x,y,z) frand(67(x)+59(y)+71*(z))
24	#define rand3b(x,y,z) frand(73(x)+79(y)+83*(z))
25	#define rand3c(x,y,z) frand(89(x)+97(y)+101*(z))
26	#define rand3d(x,y,z) frand(103(x)+107(y)+109*(z))
27
28	#define hermite(p0,p1,r0,r1,t) ( p0((2.0t-3.0)tt+1.0) + \
29	p1(-2.0t+3.0)tt + \
30	r0((t-2.0)t+1.0)*t + \
31	r1(t-1.0)t*t )
32
33	static char *noise_name[4] = {"noise3a", "noise3b", "noise3c", "noise3"};
34	static char fnoise_name[] = "fnoise3";
35	static char hermite_name[] = "hermite";
36
37	double *noise3(), fnoise3(), argument(), frand();
38
39	static long xlim[3][2];
40	static double xarg[3];
41
42	#define EPSILON .0001 /* error allowed in fractal */
43
44	#define frand3(x,y,z) frand(17(x)+23(y)+29*(z))
45
46
47	static double
48	l_noise3(nam) /* compute a noise function */
49	register char *nam;
50	{
51	register int i;
52	double x[3];
53	/* get point */
54	x[0] = argument(1);
55	x[1] = argument(2);
56	x[2] = argument(3);
57	/* make appropriate call */
58	if (nam == fnoise_name)
59	return(fnoise3(x));
60	i = 4;
61	while (i--)
62	if (nam == noise_name[i])
63	return(noise3(x)[i]);
64	eputs("Bad call to l_noise3()!\n");
65	quit(1);
66	}
67
68
69	double
70	l_hermite() /* library call for hermite interpolation */
71	{
72	double t;
73
74	t = argument(5);
75	return( hermite(argument(1), argument(2),
76	argument(3), argument(4), t) );
77	}
78
79
80	setnoisefuncs() /* add noise functions to library */
81	{
82	register int i;
83
84	funset(hermite_name, 5, ':', l_hermite);
85	funset(fnoise_name, 3, ':', l_noise3);
86	i = 4;
87	while (i--)
88	funset(noise_name[i], 3, ':', l_noise3);
89	}
90
91
92	double *
93	noise3(xnew) /* compute the noise function */
94	register double xnew[3];
95	{
96	extern double floor();
97	static double x[3] = {-100000.0, -100000.0, -100000.0};
98	static double f[4];
99
100	if (x[0]==xnew[0] && x[1]==xnew[1] && x[2]==xnew[2])
101	return(f);
102	x[0] = xnew[0]; x[1] = xnew[1]; x[2] = xnew[2];
103	xlim[0][0] = floor(x[0]); xlim[0][1] = xlim[0][0] + 1;
104	xlim[1][0] = floor(x[1]); xlim[1][1] = xlim[1][0] + 1;
105	xlim[2][0] = floor(x[2]); xlim[2][1] = xlim[2][0] + 1;
106	xarg[0] = x[0] - xlim[0][0];
107	xarg[1] = x[1] - xlim[1][0];
108	xarg[2] = x[2] - xlim[2][0];
109	interpolate(f, 0, 3);
110	return(f);
111	}
112
113
114	static
115	interpolate(f, i, n)
116	double f[4];
117	register int i, n;
118	{
119	double f0[4], f1[4];
120
121	if (n == 0) {
122	f[A] = rand3a(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
123	f[B] = rand3b(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
124	f[C] = rand3c(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
125	f[D] = rand3d(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
126	} else {
127	n--;
128	interpolate(f0, i, n);
129	interpolate(f1, i \| 1<<n, n);
130	f[A] = (1.0-xarg[n])f0[A] + xarg[n]f1[A];
131	f[B] = (1.0-xarg[n])f0[B] + xarg[n]f1[B];
132	f[C] = (1.0-xarg[n])f0[C] + xarg[n]f1[C];
133	f[D] = hermite(f0[D], f1[D], f0[n], f1[n], xarg[n]);
134	}
135	}
136
137
138	double
139	frand(s) /* get random number from seed */
140	register long s;
141	{
142	s = s<<13 ^ s;
143	return(1.0-((s(ss*15731+789221)+1376312589)&0x7fffffff)/1073741824.0);
144	}
145
146
147	double
148	fnoise3(p) /* compute fractal noise function */
149	double p[3];
150	{
151	double floor();
152	long t[3], v[3], beg[3];
153	double fval[8], fc;
154	int branch;
155	register long s;
156	register int i, j;
157	/* get starting cube */
158	s = (long)(1.0/EPSILON);
159	for (i = 0; i < 3; i++) {
160	t[i] = s*p[i];
161	beg[i] = s*floor(p[i]);
162	}
163	for (j = 0; j < 8; j++) {
164	for (i = 0; i < 3; i++) {
165	v[i] = beg[i];
166	if (j & 1<<i)
167	v[i] += s;
168	}
169	fval[j] = frand3(v[0],v[1],v[2]);
170	}
171	/* compute fractal */
172	for ( ; ; ) {
173	fc = 0.0;
174	for (j = 0; j < 8; j++)
175	fc += fval[j];
176	fc *= 0.125;
177	if ((s >>= 1) == 0)
178	return(fc); /* close enough */
179	branch = 0;
180	for (i = 0; i < 3; i++) { /* do center */
181	v[i] = beg[i] + s;
182	if (t[i] > v[i]) {
183	branch \|= 1<<i;
184	}
185	}
186	fc += sEPSILONfrand3(v[0],v[1],v[2]);
187	fval[~branch & 7] = fc;
188	for (i = 0; i < 3; i++) { /* do faces */
189	if (branch & 1<<i)
190	v[i] += s;
191	else
192	v[i] -= s;
193	fc = 0.0;
194	for (j = 0; j < 8; j++)
195	if (~(j^branch) & 1<<i)
196	fc += fval[j];
197	fc = 0.25fc + sEPSILON*frand3(v[0],v[1],v[2]);
198	fval[~(branch^1<<i) & 7] = fc;
199	v[i] = beg[i] + s;
200	}
201	for (i = 0; i < 3; i++) { /* do edges */
202	j = (i+1)%3;
203	if (branch & 1<<j)
204	v[j] += s;
205	else
206	v[j] -= s;
207	j = (i+2)%3;
208	if (branch & 1<<j)
209	v[j] += s;
210	else
211	v[j] -= s;
212	fc = fval[branch & ~(1<<i)];
213	fc += fval[branch \| 1<<i];
214	fc = 0.5fc + sEPSILON*frand3(v[0],v[1],v[2]);
215	fval[branch^1<<i] = fc;
216	j = (i+1)%3;
217	v[j] = beg[j] + s;
218	j = (i+2)%3;
219	v[j] = beg[j] + s;
220	}
221	for (i = 0; i < 3; i++) /* new cube */
222	if (branch & 1<<i)
223	beg[i] += s;
224	}
225	}