ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/noise3.c
Revision: 2.9
Committed: Tue Mar 30 16:13:01 2004 UTC (20 years 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.

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: noise3.c,v 2.8 2003/08/04 22:37:53 greg Exp $";
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 #include "copyright.h"
12
13 #include <math.h>
14
15 #include "calcomp.h"
16 #include "func.h"
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 hpoly1(t) ((2.0*t-3.0)*t*t+1.0)
29 #define hpoly2(t) (-2.0*t+3.0)*t*t
30 #define hpoly3(t) ((t-2.0)*t+1.0)*t
31 #define hpoly4(t) (t-1.0)*t*t
32
33 #define hermite(p0,p1,r0,r1,t) ( p0*hpoly1(t) + \
34 p1*hpoly2(t) + \
35 r0*hpoly3(t) + \
36 r1*hpoly4(t) )
37
38 static char noise_name[4][8] = {"noise3x", "noise3y", "noise3z", "noise3"};
39 static char fnoise_name[] = "fnoise3";
40 static char hermite_name[] = "hermite";
41
42 //double *noise3(), fnoise3(), frand();
43 //static interpolate();
44
45 static long xlim[3][2];
46 static double xarg[3];
47
48 #define EPSILON .001 /* error allowed in fractal */
49
50 #define frand3(x,y,z) frand(17*(x)+23*(y)+29*(z))
51
52 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
60 static double
61 l_noise3( /* compute a noise function */
62 register char *nam
63 )
64 {
65 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 eputs(nam);
79 eputs(": called l_noise3!\n");
80 quit(1);
81 return 1; /* pro forma return */
82 }
83
84
85 static double
86 l_hermite(char *nm) /* library call for hermite interpolation */
87 {
88 double t;
89
90 t = argument(5);
91 return( hermite(argument(1), argument(2),
92 argument(3), argument(4), t) );
93 }
94
95
96 extern void
97 setnoisefuncs(void) /* add noise functions to library */
98 {
99 register int i;
100
101 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 }
107
108
109 static double *
110 noise3( /* compute the noise function */
111 register double xnew[3]
112 )
113 {
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 static void
132 interpolate(
133 double f[4],
134 register int i,
135 register int n
136 )
137 {
138 double f0[4], f1[4], hp1, hp2;
139
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 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 }
156 }
157
158
159 static double
160 frand( /* get random number from seed */
161 register long s
162 )
163 {
164 s = s<<13 ^ s;
165 return(1.0-((s*(s*s*15731+789221)+1376312589)&0x7fffffff)/1073741824.0);
166 }
167
168
169 static double
170 fnoise3( /* compute fractal noise function */
171 double p[3]
172 )
173 {
174 long t[3], v[3], beg[3];
175 double fval[8], fc;
176 int branch;
177 register long s;
178 register int i, j;
179 /* get starting cube */
180 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 for (j = 0; j < 8; j++) {
186 for (i = 0; i < 3; i++) {
187 v[i] = beg[i];
188 if (j & 1<<i)
189 v[i] += s;
190 }
191 fval[j] = frand3(v[0],v[1],v[2]);
192 }
193 /* compute fractal */
194 for ( ; ; ) {
195 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 branch = 0;
202 for (i = 0; i < 3; i++) { /* do center */
203 v[i] = beg[i] + s;
204 if (t[i] > v[i]) {
205 branch |= 1<<i;
206 }
207 }
208 fc += s*EPSILON*frand3(v[0],v[1],v[2]);
209 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 fc = 0.25*fc + s*EPSILON*frand3(v[0],v[1],v[2]);
220 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 fc = 0.5*fc + s*EPSILON*frand3(v[0],v[1],v[2]);
237 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 }