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 )

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

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

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

#define  frand3(x,y,z)  frand((long)((12.38*(x)-22.30*(y)-42.63*(z))/EPSILON))

double  fnoise3();


double
l_noise3()                      /* compute 3-dimensional noise function */
{
        return(noise3coef(D));
}


double
l_noise3a()                     /* compute x slope of noise function */
{
        return(noise3coef(A));
}


double
l_noise3b()                     /* compute y slope of noise function */
{
        return(noise3coef(B));
}


double
l_noise3c()                     /* compute z slope of noise function */
{
        return(noise3coef(C));
}


double
l_fnoise3()                     /* compute fractal noise function */
{
        double  x[3];

        x[0] = argument(1);
        x[1] = argument(2);
        x[2] = argument(3);

        return(fnoise3(x));
}


static double
noise3coef(coef)                /* return coefficient of noise function */
int  coef;
{
        double  x[3];

        x[0] = argument(1);
        x[1] = argument(2);
        x[2] = argument(3);

        return(noise3(x)[coef]);
}


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
l_hermite()                     /* library call for hermite interpolation */
{
        double  t;
        
        t = argument(5);
        return( hermite(argument(1), argument(2),
                        argument(3), argument(4), t) );
}


double
fnoise3(p)                      /* compute fractal noise function */
register double  p[3];
{
        double  floor();
        double  v[3], beg[3], fval[8], s, fc;
        int  closing, branch;
        register int  i, j;
                                                /* get starting cube */
        for (i = 0; i < 3; i++)
                beg[i] = 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] += 1.0;
                }
                fval[j] = frand3(v[0],v[1],v[2]);
        }
        s = 1.0;
                                                /* compute fractal */
        for ( ; ; ) {
                s *= 0.5;
                branch = 0;
                closing = 0;
                for (i = 0; i < 3; i++) {       /* do center */
                        v[i] = beg[i] + s;
                        if (p[i] > v[i]) {
                                branch |= 1<<i;
                                if (p[i] - v[i] > EPSILON)
                                        closing++;
                        } else if (v[i] - p[i] > EPSILON)
                                closing++;
                }
                fc = 0.0;
                for (j = 0; j < 8; j++)
                        fc += fval[j];
                fc = 0.125*fc + s*frand3(v[0],v[1],v[2]);
                if (closing == 0)
                        return(fc);             /* close enough */
                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*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*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.2
Committed:	Fri Mar 3 12:25:34 1989 UTC (35 years, 2 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 1.1:	+1 -1 lines
Log Message:	smaller EPSILON for fnoise3
#	User	Rev	Content
1	greg	1.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			double *noise3(), noise3coef(), argument(), frand();
34
35			static long xlim[3][2];
36			static double xarg[3];
37
38	greg	1.2	#define EPSILON .0001 /* error allowed in fractal */
39	greg	1.1
40			#define frand3(x,y,z) frand((long)((12.38(x)-22.30(y)-42.63*(z))/EPSILON))
41
42			double fnoise3();
43
44
45			double
46			l_noise3() /* compute 3-dimensional noise function */
47			{
48			return(noise3coef(D));
49			}
50
51
52			double
53			l_noise3a() /* compute x slope of noise function */
54			{
55			return(noise3coef(A));
56			}
57
58
59			double
60			l_noise3b() /* compute y slope of noise function */
61			{
62			return(noise3coef(B));
63			}
64
65
66			double
67			l_noise3c() /* compute z slope of noise function */
68			{
69			return(noise3coef(C));
70			}
71
72
73			double
74			l_fnoise3() /* compute fractal noise function */
75			{
76			double x[3];
77
78			x[0] = argument(1);
79			x[1] = argument(2);
80			x[2] = argument(3);
81
82			return(fnoise3(x));
83			}
84
85
86			static double
87			noise3coef(coef) /* return coefficient of noise function */
88			int coef;
89			{
90			double x[3];
91
92			x[0] = argument(1);
93			x[1] = argument(2);
94			x[2] = argument(3);
95
96			return(noise3(x)[coef]);
97			}
98
99
100			double *
101			noise3(xnew) /* compute the noise function */
102			register double xnew[3];
103			{
104			extern double floor();
105			static double x[3] = {-100000.0, -100000.0, -100000.0};
106			static double f[4];
107
108			if (x[0]==xnew[0] && x[1]==xnew[1] && x[2]==xnew[2])
109			return(f);
110			x[0] = xnew[0]; x[1] = xnew[1]; x[2] = xnew[2];
111			xlim[0][0] = floor(x[0]); xlim[0][1] = xlim[0][0] + 1;
112			xlim[1][0] = floor(x[1]); xlim[1][1] = xlim[1][0] + 1;
113			xlim[2][0] = floor(x[2]); xlim[2][1] = xlim[2][0] + 1;
114			xarg[0] = x[0] - xlim[0][0];
115			xarg[1] = x[1] - xlim[1][0];
116			xarg[2] = x[2] - xlim[2][0];
117			interpolate(f, 0, 3);
118			return(f);
119			}
120
121
122			static
123			interpolate(f, i, n)
124			double f[4];
125			register int i, n;
126			{
127			double f0[4], f1[4];
128
129			if (n == 0) {
130			f[A] = rand3a(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
131			f[B] = rand3b(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
132			f[C] = rand3c(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
133			f[D] = rand3d(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
134			} else {
135			n--;
136			interpolate(f0, i, n);
137			interpolate(f1, i \| 1<<n, n);
138			f[A] = (1.0-xarg[n])f0[A] + xarg[n]f1[A];
139			f[B] = (1.0-xarg[n])f0[B] + xarg[n]f1[B];
140			f[C] = (1.0-xarg[n])f0[C] + xarg[n]f1[C];
141			f[D] = hermite(f0[D], f1[D], f0[n], f1[n], xarg[n]);
142			}
143			}
144
145
146			double
147			frand(s) /* get random number from seed */
148			register long s;
149			{
150			s = s<<13 ^ s;
151			return(1.0-((s(ss*15731+789221)+1376312589)&0x7fffffff)/1073741824.0);
152			}
153
154
155			double
156			l_hermite() /* library call for hermite interpolation */
157			{
158			double t;
159
160			t = argument(5);
161			return( hermite(argument(1), argument(2),
162			argument(3), argument(4), t) );
163			}
164
165
166			double
167			fnoise3(p) /* compute fractal noise function */
168			register double p[3];
169			{
170			double floor();
171			double v[3], beg[3], fval[8], s, fc;
172			int closing, branch;
173			register int i, j;
174			/* get starting cube */
175			for (i = 0; i < 3; i++)
176			beg[i] = floor(p[i]);
177			for (j = 0; j < 8; j++) {
178			for (i = 0; i < 3; i++) {
179			v[i] = beg[i];
180			if (j & 1<<i)
181			v[i] += 1.0;
182			}
183			fval[j] = frand3(v[0],v[1],v[2]);
184			}
185			s = 1.0;
186			/* compute fractal */
187			for ( ; ; ) {
188			s *= 0.5;
189			branch = 0;
190			closing = 0;
191			for (i = 0; i < 3; i++) { /* do center */
192			v[i] = beg[i] + s;
193			if (p[i] > v[i]) {
194			branch \|= 1<<i;
195			if (p[i] - v[i] > EPSILON)
196			closing++;
197			} else if (v[i] - p[i] > EPSILON)
198			closing++;
199			}
200			fc = 0.0;
201			for (j = 0; j < 8; j++)
202			fc += fval[j];
203			fc = 0.125fc + sfrand3(v[0],v[1],v[2]);
204			if (closing == 0)
205			return(fc); /* close enough */
206			fval[~branch & 7] = fc;
207			for (i = 0; i < 3; i++) { /* do faces */
208			if (branch & 1<<i)
209			v[i] += s;
210			else
211			v[i] -= s;
212			fc = 0.0;
213			for (j = 0; j < 8; j++)
214			if (~(j^branch) & 1<<i)
215			fc += fval[j];
216			fc = 0.25fc + sfrand3(v[0],v[1],v[2]);
217			fval[~(branch^1<<i) & 7] = fc;
218			v[i] = beg[i] + s;
219			}
220			for (i = 0; i < 3; i++) { /* do edges */
221			j = (i+1)%3;
222			if (branch & 1<<j)
223			v[j] += s;
224			else
225			v[j] -= s;
226			j = (i+2)%3;
227			if (branch & 1<<j)
228			v[j] += s;
229			else
230			v[j] -= s;
231			fc = fval[branch & ~(1<<i)];
232			fc += fval[branch \| 1<<i];
233			fc = 0.5fc + sfrand3(v[0],v[1],v[2]);
234			fval[branch^1<<i] = fc;
235			j = (i+1)%3;
236			v[j] = beg[j] + s;
237			j = (i+2)%3;
238			v[j] = beg[j] + s;
239			}
240			for (i = 0; i < 3; i++) /* new cube */
241			if (branch & 1<<i)
242			beg[i] += s;
243			}
244			}