ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/noise3.c
(Generate patch)

Comparing ray/src/rt/noise3.c (file contents):
Revision 1.2 by greg, Fri Mar 3 12:25:34 1989 UTC vs.
Revision 2.8 by greg, Mon Aug 4 22:37:53 2003 UTC

# Line 1 | Line 1
1 /* Copyright (c) 1988 Regents of the University of California */
2
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ LBL";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6
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)
12 *
13 *     4/15/86
14 *     5/19/88  Added fractal noise function
9   */
10  
11 + #include "copyright.h"
12  
13 + #include  "calcomp.h"
14 +
15 + #include  <math.h>
16 +
17   #define  A              0
18   #define  B              1
19   #define  C              2
# Line 25 | Line 24 | static char SCCSid[] = "$SunId$ LBL";
24   #define  rand3c(x,y,z)  frand(89*(x)+97*(y)+101*(z))
25   #define  rand3d(x,y,z)  frand(103*(x)+107*(y)+109*(z))
26  
27 < #define  hermite(p0,p1,r0,r1,t)  (      p0*((2.0*t-3.0)*t*t+1.0) + \
28 <                                        p1*(-2.0*t+3.0)*t*t + \
29 <                                        r0*((t-2.0)*t+1.0)*t + \
30 <                                        r1*(t-1.0)*t*t )
27 > #define  hpoly1(t)      ((2.0*t-3.0)*t*t+1.0)
28 > #define  hpoly2(t)      (-2.0*t+3.0)*t*t
29 > #define  hpoly3(t)      ((t-2.0)*t+1.0)*t
30 > #define  hpoly4(t)      (t-1.0)*t*t
31  
32 < double  *noise3(), noise3coef(), argument(), frand();
32 > #define  hermite(p0,p1,r0,r1,t)  (      p0*hpoly1(t) + \
33 >                                        p1*hpoly2(t) + \
34 >                                        r0*hpoly3(t) + \
35 >                                        r1*hpoly4(t) )
36  
37 + static char  noise_name[4][8] = {"noise3x", "noise3y", "noise3z", "noise3"};
38 + static char  fnoise_name[] = "fnoise3";
39 + static char  hermite_name[] = "hermite";
40 +
41 + double  *noise3(), fnoise3(), frand();
42 + static  interpolate();
43 +
44   static long  xlim[3][2];
45   static double  xarg[3];
46  
47 < #define  EPSILON        .0001           /* error allowed in fractal */
47 > #define  EPSILON        .001            /* error allowed in fractal */
48  
49 < #define  frand3(x,y,z)  frand((long)((12.38*(x)-22.30*(y)-42.63*(z))/EPSILON))
49 > #define  frand3(x,y,z)  frand(17*(x)+23*(y)+29*(z))
50  
42 double  fnoise3();
51  
52 <
53 < double
54 < l_noise3()                      /* compute 3-dimensional noise function */
52 > static double
53 > l_noise3(nam)                   /* compute a noise function */
54 > register char  *nam;
55   {
56 <        return(noise3coef(D));
56 >        register int  i;
57 >        double  x[3];
58 >                                        /* get point */
59 >        x[0] = argument(1);
60 >        x[1] = argument(2);
61 >        x[2] = argument(3);
62 >                                        /* make appropriate call */
63 >        if (nam == fnoise_name)
64 >                return(fnoise3(x));
65 >        i = 4;
66 >        while (i--)
67 >                if (nam == noise_name[i])
68 >                        return(noise3(x)[i]);
69 >        eputs(nam);
70 >        eputs(": called l_noise3!\n");
71 >        quit(1);
72   }
73  
74  
75   double
76 < l_noise3a()                     /* compute x slope of noise function */
76 > l_hermite(char *nm)             /* library call for hermite interpolation */
77   {
78 <        return(noise3coef(A));
78 >        double  t;
79 >        
80 >        t = argument(5);
81 >        return( hermite(argument(1), argument(2),
82 >                        argument(3), argument(4), t) );
83   }
84  
85  
86 < double
60 < l_noise3b()                     /* compute y slope of noise function */
86 > setnoisefuncs()                 /* add noise functions to library */
87   {
88 <        return(noise3coef(B));
63 < }
88 >        register int  i;
89  
90 <
91 < double
92 < l_noise3c()                     /* compute z slope of noise function */
93 < {
94 <        return(noise3coef(C));
90 >        funset(hermite_name, 5, ':', l_hermite);
91 >        funset(fnoise_name, 3, ':', l_noise3);
92 >        i = 4;
93 >        while (i--)
94 >                funset(noise_name[i], 3, ':', l_noise3);
95   }
96  
97  
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
98   double *
99   noise3(xnew)                    /* compute the noise function */
100   register double  xnew[3];
101   {
104        extern double  floor();
102          static double  x[3] = {-100000.0, -100000.0, -100000.0};
103          static double  f[4];
104  
# Line 124 | Line 121 | interpolate(f, i, n)
121   double  f[4];
122   register int  i, n;
123   {
124 <        double  f0[4], f1[4];
124 >        double  f0[4], f1[4], hp1, hp2;
125  
126          if (n == 0) {
127                  f[A] = rand3a(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
# Line 135 | Line 132 | register int  i, n;
132                  n--;
133                  interpolate(f0, i, n);
134                  interpolate(f1, i | 1<<n, n);
135 <                f[A] = (1.0-xarg[n])*f0[A] + xarg[n]*f1[A];
136 <                f[B] = (1.0-xarg[n])*f0[B] + xarg[n]*f1[B];
137 <                f[C] = (1.0-xarg[n])*f0[C] + xarg[n]*f1[C];
138 <                f[D] = hermite(f0[D], f1[D], f0[n], f1[n], xarg[n]);
135 >                hp1 = hpoly1(xarg[n]); hp2 = hpoly2(xarg[n]);
136 >                f[A] = f0[A]*hp1 + f1[A]*hp2;
137 >                f[B] = f0[B]*hp1 + f1[B]*hp2;
138 >                f[C] = f0[C]*hp1 + f1[C]*hp2;
139 >                f[D] = f0[D]*hp1 + f1[D]*hp2 +
140 >                                f0[n]*hpoly3(xarg[n]) + f1[n]*hpoly4(xarg[n]);
141          }
142   }
143  
# Line 153 | Line 152 | register long  s;
152  
153  
154   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
155   fnoise3(p)                      /* compute fractal noise function */
156 < register double  p[3];
156 > double  p[3];
157   {
158 <        double  floor();
159 <        double  v[3], beg[3], fval[8], s, fc;
160 <        int  closing, branch;
158 >        long  t[3], v[3], beg[3];
159 >        double  fval[8], fc;
160 >        int  branch;
161 >        register long  s;
162          register int  i, j;
163                                                  /* get starting cube */
164 <        for (i = 0; i < 3; i++)
165 <                beg[i] = floor(p[i]);
164 >        s = (long)(1.0/EPSILON);
165 >        for (i = 0; i < 3; i++) {
166 >                t[i] = s*p[i];
167 >                beg[i] = s*floor(p[i]);
168 >        }
169          for (j = 0; j < 8; j++) {
170                  for (i = 0; i < 3; i++) {
171                          v[i] = beg[i];
172                          if (j & 1<<i)
173 <                                v[i] += 1.0;
173 >                                v[i] += s;
174                  }
175                  fval[j] = frand3(v[0],v[1],v[2]);
176          }
185        s = 1.0;
177                                                  /* compute fractal */
178          for ( ; ; ) {
179 <                s *= 0.5;
179 >                fc = 0.0;
180 >                for (j = 0; j < 8; j++)
181 >                        fc += fval[j];
182 >                fc *= 0.125;
183 >                if ((s >>= 1) == 0)
184 >                        return(fc);             /* close enough */
185                  branch = 0;
190                closing = 0;
186                  for (i = 0; i < 3; i++) {       /* do center */
187                          v[i] = beg[i] + s;
188 <                        if (p[i] > v[i]) {
188 >                        if (t[i] > v[i]) {
189                                  branch |= 1<<i;
190 <                                if (p[i] - v[i] > EPSILON)
196 <                                        closing++;
197 <                        } else if (v[i] - p[i] > EPSILON)
198 <                                closing++;
190 >                        }
191                  }
192 <                fc = 0.0;
201 <                for (j = 0; j < 8; j++)
202 <                        fc += fval[j];
203 <                fc = 0.125*fc + s*frand3(v[0],v[1],v[2]);
204 <                if (closing == 0)
205 <                        return(fc);             /* close enough */
192 >                fc += s*EPSILON*frand3(v[0],v[1],v[2]);
193                  fval[~branch & 7] = fc;
194                  for (i = 0; i < 3; i++) {       /* do faces */
195                          if (branch & 1<<i)
# Line 213 | Line 200 | register double  p[3];
200                          for (j = 0; j < 8; j++)
201                                  if (~(j^branch) & 1<<i)
202                                          fc += fval[j];
203 <                        fc = 0.25*fc + s*frand3(v[0],v[1],v[2]);
203 >                        fc = 0.25*fc + s*EPSILON*frand3(v[0],v[1],v[2]);
204                          fval[~(branch^1<<i) & 7] = fc;
205                          v[i] = beg[i] + s;
206                  }
# Line 230 | Line 217 | register double  p[3];
217                                  v[j] -= s;
218                          fc = fval[branch & ~(1<<i)];
219                          fc += fval[branch | 1<<i];
220 <                        fc = 0.5*fc + s*frand3(v[0],v[1],v[2]);
220 >                        fc = 0.5*fc + s*EPSILON*frand3(v[0],v[1],v[2]);
221                          fval[branch^1<<i] = fc;
222                          j = (i+1)%3;
223                          v[j] = beg[j] + s;

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines