src/rt/noise3.c

#ifndef lint
static const char       RCSid[] = "$Id: noise3.c,v 2.8 2003/08/04 22:37:53 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)
 */

#include "copyright.h"

#include  <math.h>

#include  "calcomp.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";

//double  *noise3(), fnoise3(), frand();
//static  interpolate();

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

#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(                       /* 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(nam);
        eputs(": called l_noise3!\n");
        quit(1);
        return 1; /* pro forma return */
}


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
setnoisefuncs(void)                     /* 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);
}


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
)
{
        s = s<<13 ^ s;
        return(1.0-((s*(s*s*15731+789221)+1376312589)&0x7fffffff)/1073741824.0);
}


static double
fnoise3(                        /* compute fractal noise function */
        double  p[3]
)
{
        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:	2.9
Committed:	Tue Mar 30 16:13:01 2004 UTC (20 years, 1 month ago) by schorsch
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad3R7P2, rad3R7P1, rad3R6, rad3R6P1
Changes since 2.8:	+38 -22 lines
Log Message:	Continued ANSIfication. There are only bits and pieces left now.
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	schorsch	2.9	static const char RCSid[] = "$Id: noise3.c,v 2.8 2003/08/04 22:37:53 greg Exp $";
3	greg	1.1	#endif
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)
9	greg	2.6	*/
10
11	greg	2.7	#include "copyright.h"
12	greg	1.1
13	schorsch	2.9	#include <math.h>
14
15	greg	2.8	#include "calcomp.h"
16	schorsch	2.9	#include "func.h"
17	greg	1.1
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	greg	1.7	#define hpoly1(t) ((2.0t-3.0)t*t+1.0)
29			#define hpoly2(t) (-2.0t+3.0)t*t
30			#define hpoly3(t) ((t-2.0)t+1.0)t
31			#define hpoly4(t) (t-1.0)tt
32	greg	1.1
33	greg	1.7	#define hermite(p0,p1,r0,r1,t) ( p0*hpoly1(t) + \
34			p1*hpoly2(t) + \
35			r0*hpoly3(t) + \
36			r1*hpoly4(t) )
37
38	greg	2.5	static char noise_name[4][8] = {"noise3x", "noise3y", "noise3z", "noise3"};
39	greg	1.5	static char fnoise_name[] = "fnoise3";
40			static char hermite_name[] = "hermite";
41	greg	1.1
42	schorsch	2.9	//double *noise3(), fnoise3(), frand();
43			//static interpolate();
44	greg	1.5
45	greg	1.1	static long xlim[3][2];
46			static double xarg[3];
47
48	greg	2.4	#define EPSILON .001 /* error allowed in fractal */
49	greg	1.1
50	greg	1.3	#define frand3(x,y,z) frand(17(x)+23(y)+29*(z))
51	greg	1.1
52	schorsch	2.9	static double l_noise3(char *nam);
53			static double l_hermite(char *nm);
54			static double * noise3(double xnew[3]);
55			static void interpolate(double f[4], int i, int n);
56			static double frand(long s);
57			static double fnoise3(double p[3]);
58
59	greg	1.1
60	greg	1.5	static double
61	schorsch	2.9	l_noise3( /* compute a noise function */
62			register char *nam
63			)
64	greg	1.1	{
65	greg	1.5	register int i;
66			double x[3];
67			/* get point */
68			x[0] = argument(1);
69			x[1] = argument(2);
70			x[2] = argument(3);
71			/* make appropriate call */
72			if (nam == fnoise_name)
73			return(fnoise3(x));
74			i = 4;
75			while (i--)
76			if (nam == noise_name[i])
77			return(noise3(x)[i]);
78	greg	1.6	eputs(nam);
79			eputs(": called l_noise3!\n");
80	greg	1.5	quit(1);
81	schorsch	2.9	return 1; /* pro forma return */
82	greg	1.1	}
83
84
85	schorsch	2.9	static double
86	greg	2.8	l_hermite(char nm) / library call for hermite interpolation */
87	greg	1.1	{
88	greg	1.5	double t;
89
90			t = argument(5);
91			return( hermite(argument(1), argument(2),
92			argument(3), argument(4), t) );
93	greg	1.1	}
94
95
96	schorsch	2.9	extern void
97			setnoisefuncs(void) /* add noise functions to library */
98	greg	1.1	{
99	greg	1.5	register int i;
100	greg	1.1
101	greg	1.5	funset(hermite_name, 5, ':', l_hermite);
102			funset(fnoise_name, 3, ':', l_noise3);
103			i = 4;
104			while (i--)
105			funset(noise_name[i], 3, ':', l_noise3);
106	greg	1.1	}
107
108
109	schorsch	2.9	static double *
110			noise3( /* compute the noise function */
111			register double xnew[3]
112			)
113	greg	1.1	{
114			static double x[3] = {-100000.0, -100000.0, -100000.0};
115			static double f[4];
116
117			if (x[0]==xnew[0] && x[1]==xnew[1] && x[2]==xnew[2])
118			return(f);
119			x[0] = xnew[0]; x[1] = xnew[1]; x[2] = xnew[2];
120			xlim[0][0] = floor(x[0]); xlim[0][1] = xlim[0][0] + 1;
121			xlim[1][0] = floor(x[1]); xlim[1][1] = xlim[1][0] + 1;
122			xlim[2][0] = floor(x[2]); xlim[2][1] = xlim[2][0] + 1;
123			xarg[0] = x[0] - xlim[0][0];
124			xarg[1] = x[1] - xlim[1][0];
125			xarg[2] = x[2] - xlim[2][0];
126			interpolate(f, 0, 3);
127			return(f);
128			}
129
130
131	schorsch	2.9	static void
132			interpolate(
133			double f[4],
134			register int i,
135			register int n
136			)
137	greg	1.1	{
138	greg	1.7	double f0[4], f1[4], hp1, hp2;
139	greg	1.1
140			if (n == 0) {
141			f[A] = rand3a(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
142			f[B] = rand3b(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
143			f[C] = rand3c(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
144			f[D] = rand3d(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
145			} else {
146			n--;
147			interpolate(f0, i, n);
148			interpolate(f1, i \| 1<<n, n);
149	greg	1.7	hp1 = hpoly1(xarg[n]); hp2 = hpoly2(xarg[n]);
150			f[A] = f0[A]hp1 + f1[A]hp2;
151			f[B] = f0[B]hp1 + f1[B]hp2;
152			f[C] = f0[C]hp1 + f1[C]hp2;
153			f[D] = f0[D]hp1 + f1[D]hp2 +
154			f0[n]hpoly3(xarg[n]) + f1[n]hpoly4(xarg[n]);
155	greg	1.1	}
156			}
157
158
159	schorsch	2.9	static double
160			frand( /* get random number from seed */
161			register long s
162			)
163	greg	1.1	{
164			s = s<<13 ^ s;
165			return(1.0-((s(ss*15731+789221)+1376312589)&0x7fffffff)/1073741824.0);
166			}
167
168
169	schorsch	2.9	static double
170			fnoise3( /* compute fractal noise function */
171			double p[3]
172			)
173	greg	1.1	{
174	greg	1.4	long t[3], v[3], beg[3];
175	greg	1.3	double fval[8], fc;
176			int branch;
177	greg	1.4	register long s;
178	greg	1.1	register int i, j;
179			/* get starting cube */
180	greg	1.3	s = (long)(1.0/EPSILON);
181			for (i = 0; i < 3; i++) {
182			t[i] = s*p[i];
183			beg[i] = s*floor(p[i]);
184			}
185	greg	1.1	for (j = 0; j < 8; j++) {
186			for (i = 0; i < 3; i++) {
187			v[i] = beg[i];
188			if (j & 1<<i)
189	greg	1.3	v[i] += s;
190	greg	1.1	}
191			fval[j] = frand3(v[0],v[1],v[2]);
192			}
193			/* compute fractal */
194			for ( ; ; ) {
195	greg	1.4	fc = 0.0;
196			for (j = 0; j < 8; j++)
197			fc += fval[j];
198			fc *= 0.125;
199			if ((s >>= 1) == 0)
200			return(fc); /* close enough */
201	greg	1.1	branch = 0;
202			for (i = 0; i < 3; i++) { /* do center */
203			v[i] = beg[i] + s;
204	greg	1.3	if (t[i] > v[i]) {
205	greg	1.1	branch \|= 1<<i;
206	greg	1.3	}
207	greg	1.1	}
208	greg	1.3	fc += sEPSILONfrand3(v[0],v[1],v[2]);
209	greg	1.1	fval[~branch & 7] = fc;
210			for (i = 0; i < 3; i++) { /* do faces */
211			if (branch & 1<<i)
212			v[i] += s;
213			else
214			v[i] -= s;
215			fc = 0.0;
216			for (j = 0; j < 8; j++)
217			if (~(j^branch) & 1<<i)
218			fc += fval[j];
219	greg	1.3	fc = 0.25fc + sEPSILON*frand3(v[0],v[1],v[2]);
220	greg	1.1	fval[~(branch^1<<i) & 7] = fc;
221			v[i] = beg[i] + s;
222			}
223			for (i = 0; i < 3; i++) { /* do edges */
224			j = (i+1)%3;
225			if (branch & 1<<j)
226			v[j] += s;
227			else
228			v[j] -= s;
229			j = (i+2)%3;
230			if (branch & 1<<j)
231			v[j] += s;
232			else
233			v[j] -= s;
234			fc = fval[branch & ~(1<<i)];
235			fc += fval[branch \| 1<<i];
236	greg	1.3	fc = 0.5fc + sEPSILON*frand3(v[0],v[1],v[2]);
237	greg	1.1	fval[branch^1<<i] = fc;
238			j = (i+1)%3;
239			v[j] = beg[j] + s;
240			j = (i+2)%3;
241			v[j] = beg[j] + s;
242			}
243			for (i = 0; i < 3; i++) /* new cube */
244			if (branch & 1<<i)
245			beg[i] += s;
246			}
247			}