src/rt/noise3.c

#ifndef lint
static const char       RCSid[] = "$Id$";
#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)
 */

/* ====================================================================
 * The Radiance Software License, Version 1.0
 *
 * Copyright (c) 1990 - 2002 The Regents of the University of California,
 * through Lawrence Berkeley National Laboratory.   All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *         notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *       notice, this list of conditions and the following disclaimer in
 *       the documentation and/or other materials provided with the
 *       distribution.
 *
 * 3. The end-user documentation included with the redistribution,
 *           if any, must include the following acknowledgment:
 *             "This product includes Radiance software
 *                 (http://radsite.lbl.gov/)
 *                 developed by the Lawrence Berkeley National Laboratory
 *               (http://www.lbl.gov/)."
 *       Alternately, this acknowledgment may appear in the software itself,
 *       if and wherever such third-party acknowledgments normally appear.
 *
 * 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
 *       and "The Regents of the University of California" must
 *       not be used to endorse or promote products derived from this
 *       software without prior written permission. For written
 *       permission, please contact [email protected].
 *
 * 5. Products derived from this software may not be called "Radiance",
 *       nor may "Radiance" appear in their name, without prior written
 *       permission of Lawrence Berkeley National Laboratory.
 *
 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED.   IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
 * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
 * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 * ====================================================================
 *
 * This software consists of voluntary contributions made by many
 * individuals on behalf of Lawrence Berkeley National Laboratory.   For more
 * information on Lawrence Berkeley National Laboratory, please see
 * <http://www.lbl.gov/>.
 */

#include  <math.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(), argument(), 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(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(nam);
        eputs(": called 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];
{
        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], 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]);
        }
}


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];
{
        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.6
Committed:	Sat Feb 22 02:07:29 2003 UTC (21 years, 2 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.5:	+57 -6 lines
Log Message:	Changes and check-in for 3.5 release Includes new source files and modifications not recorded for many years See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release
#	Content
1	#ifndef lint
2	static const char RCSid[] = "$Id$";
3	#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	*/
10
11	/* ====================================================================
12	* The Radiance Software License, Version 1.0
13	*
14	* Copyright (c) 1990 - 2002 The Regents of the University of California,
15	* through Lawrence Berkeley National Laboratory. All rights reserved.
16	*
17	* Redistribution and use in source and binary forms, with or without
18	* modification, are permitted provided that the following conditions
19	* are met:
20	*
21	* 1. Redistributions of source code must retain the above copyright
22	* notice, this list of conditions and the following disclaimer.
23	*
24	* 2. Redistributions in binary form must reproduce the above copyright
25	* notice, this list of conditions and the following disclaimer in
26	* the documentation and/or other materials provided with the
27	* distribution.
28	*
29	* 3. The end-user documentation included with the redistribution,
30	* if any, must include the following acknowledgment:
31	* "This product includes Radiance software
32	* (http://radsite.lbl.gov/)
33	* developed by the Lawrence Berkeley National Laboratory
34	* (http://www.lbl.gov/)."
35	* Alternately, this acknowledgment may appear in the software itself,
36	* if and wherever such third-party acknowledgments normally appear.
37	*
38	* 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
39	* and "The Regents of the University of California" must
40	* not be used to endorse or promote products derived from this
41	* software without prior written permission. For written
42	* permission, please contact [email protected].
43	*
44	* 5. Products derived from this software may not be called "Radiance",
45	* nor may "Radiance" appear in their name, without prior written
46	* permission of Lawrence Berkeley National Laboratory.
47	*
48	* THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
49	* WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
50	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
51	* DISCLAIMED. IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
52	* ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
53	* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
54	* LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
55	* USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
56	* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
57	* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
58	* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
59	* SUCH DAMAGE.
60	* ====================================================================
61	*
62	* This software consists of voluntary contributions made by many
63	* individuals on behalf of Lawrence Berkeley National Laboratory. For more
64	* information on Lawrence Berkeley National Laboratory, please see
65	* <http://www.lbl.gov/>.
66	*/
67
68	#include <math.h>
69
70	#define A 0
71	#define B 1
72	#define C 2
73	#define D 3
74
75	#define rand3a(x,y,z) frand(67(x)+59(y)+71*(z))
76	#define rand3b(x,y,z) frand(73(x)+79(y)+83*(z))
77	#define rand3c(x,y,z) frand(89(x)+97(y)+101*(z))
78	#define rand3d(x,y,z) frand(103(x)+107(y)+109*(z))
79
80	#define hpoly1(t) ((2.0t-3.0)t*t+1.0)
81	#define hpoly2(t) (-2.0t+3.0)t*t
82	#define hpoly3(t) ((t-2.0)t+1.0)t
83	#define hpoly4(t) (t-1.0)tt
84
85	#define hermite(p0,p1,r0,r1,t) ( p0*hpoly1(t) + \
86	p1*hpoly2(t) + \
87	r0*hpoly3(t) + \
88	r1*hpoly4(t) )
89
90	static char noise_name[4][8] = {"noise3x", "noise3y", "noise3z", "noise3"};
91	static char fnoise_name[] = "fnoise3";
92	static char hermite_name[] = "hermite";
93
94	double *noise3(), fnoise3(), argument(), frand();
95	static interpolate();
96
97	static long xlim[3][2];
98	static double xarg[3];
99
100	#define EPSILON .001 /* error allowed in fractal */
101
102	#define frand3(x,y,z) frand(17(x)+23(y)+29*(z))
103
104
105	static double
106	l_noise3(nam) /* compute a noise function */
107	register char *nam;
108	{
109	register int i;
110	double x[3];
111	/* get point */
112	x[0] = argument(1);
113	x[1] = argument(2);
114	x[2] = argument(3);
115	/* make appropriate call */
116	if (nam == fnoise_name)
117	return(fnoise3(x));
118	i = 4;
119	while (i--)
120	if (nam == noise_name[i])
121	return(noise3(x)[i]);
122	eputs(nam);
123	eputs(": called l_noise3!\n");
124	quit(1);
125	}
126
127
128	double
129	l_hermite() /* library call for hermite interpolation */
130	{
131	double t;
132
133	t = argument(5);
134	return( hermite(argument(1), argument(2),
135	argument(3), argument(4), t) );
136	}
137
138
139	setnoisefuncs() /* add noise functions to library */
140	{
141	register int i;
142
143	funset(hermite_name, 5, ':', l_hermite);
144	funset(fnoise_name, 3, ':', l_noise3);
145	i = 4;
146	while (i--)
147	funset(noise_name[i], 3, ':', l_noise3);
148	}
149
150
151	double *
152	noise3(xnew) /* compute the noise function */
153	register double xnew[3];
154	{
155	static double x[3] = {-100000.0, -100000.0, -100000.0};
156	static double f[4];
157
158	if (x[0]==xnew[0] && x[1]==xnew[1] && x[2]==xnew[2])
159	return(f);
160	x[0] = xnew[0]; x[1] = xnew[1]; x[2] = xnew[2];
161	xlim[0][0] = floor(x[0]); xlim[0][1] = xlim[0][0] + 1;
162	xlim[1][0] = floor(x[1]); xlim[1][1] = xlim[1][0] + 1;
163	xlim[2][0] = floor(x[2]); xlim[2][1] = xlim[2][0] + 1;
164	xarg[0] = x[0] - xlim[0][0];
165	xarg[1] = x[1] - xlim[1][0];
166	xarg[2] = x[2] - xlim[2][0];
167	interpolate(f, 0, 3);
168	return(f);
169	}
170
171
172	static
173	interpolate(f, i, n)
174	double f[4];
175	register int i, n;
176	{
177	double f0[4], f1[4], hp1, hp2;
178
179	if (n == 0) {
180	f[A] = rand3a(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
181	f[B] = rand3b(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
182	f[C] = rand3c(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
183	f[D] = rand3d(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
184	} else {
185	n--;
186	interpolate(f0, i, n);
187	interpolate(f1, i \| 1<<n, n);
188	hp1 = hpoly1(xarg[n]); hp2 = hpoly2(xarg[n]);
189	f[A] = f0[A]hp1 + f1[A]hp2;
190	f[B] = f0[B]hp1 + f1[B]hp2;
191	f[C] = f0[C]hp1 + f1[C]hp2;
192	f[D] = f0[D]hp1 + f1[D]hp2 +
193	f0[n]hpoly3(xarg[n]) + f1[n]hpoly4(xarg[n]);
194	}
195	}
196
197
198	double
199	frand(s) /* get random number from seed */
200	register long s;
201	{
202	s = s<<13 ^ s;
203	return(1.0-((s(ss*15731+789221)+1376312589)&0x7fffffff)/1073741824.0);
204	}
205
206
207	double
208	fnoise3(p) /* compute fractal noise function */
209	double p[3];
210	{
211	long t[3], v[3], beg[3];
212	double fval[8], fc;
213	int branch;
214	register long s;
215	register int i, j;
216	/* get starting cube */
217	s = (long)(1.0/EPSILON);
218	for (i = 0; i < 3; i++) {
219	t[i] = s*p[i];
220	beg[i] = s*floor(p[i]);
221	}
222	for (j = 0; j < 8; j++) {
223	for (i = 0; i < 3; i++) {
224	v[i] = beg[i];
225	if (j & 1<<i)
226	v[i] += s;
227	}
228	fval[j] = frand3(v[0],v[1],v[2]);
229	}
230	/* compute fractal */
231	for ( ; ; ) {
232	fc = 0.0;
233	for (j = 0; j < 8; j++)
234	fc += fval[j];
235	fc *= 0.125;
236	if ((s >>= 1) == 0)
237	return(fc); /* close enough */
238	branch = 0;
239	for (i = 0; i < 3; i++) { /* do center */
240	v[i] = beg[i] + s;
241	if (t[i] > v[i]) {
242	branch \|= 1<<i;
243	}
244	}
245	fc += sEPSILONfrand3(v[0],v[1],v[2]);
246	fval[~branch & 7] = fc;
247	for (i = 0; i < 3; i++) { /* do faces */
248	if (branch & 1<<i)
249	v[i] += s;
250	else
251	v[i] -= s;
252	fc = 0.0;
253	for (j = 0; j < 8; j++)
254	if (~(j^branch) & 1<<i)
255	fc += fval[j];
256	fc = 0.25fc + sEPSILON*frand3(v[0],v[1],v[2]);
257	fval[~(branch^1<<i) & 7] = fc;
258	v[i] = beg[i] + s;
259	}
260	for (i = 0; i < 3; i++) { /* do edges */
261	j = (i+1)%3;
262	if (branch & 1<<j)
263	v[j] += s;
264	else
265	v[j] -= s;
266	j = (i+2)%3;
267	if (branch & 1<<j)
268	v[j] += s;
269	else
270	v[j] -= s;
271	fc = fval[branch & ~(1<<i)];
272	fc += fval[branch \| 1<<i];
273	fc = 0.5fc + sEPSILON*frand3(v[0],v[1],v[2]);
274	fval[branch^1<<i] = fc;
275	j = (i+1)%3;
276	v[j] = beg[j] + s;
277	j = (i+2)%3;
278	v[j] = beg[j] + s;
279	}
280	for (i = 0; i < 3; i++) /* new cube */
281	if (branch & 1<<i)
282	beg[i] += s;
283	}
284	}