ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/mkillum2.c
Revision: 1.10
Committed: Tue Aug 13 13:45:18 1991 UTC (32 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.9: +26 -34 lines
Log Message:
changed from urind() to peano() for n-dimensional sampling

File Contents

# Content
1 /* Copyright (c) 1991 Regents of the University of California */
2
3 #ifndef lint
4 static char SCCSid[] = "$SunId$ LBL";
5 #endif
6
7 /*
8 * Routines to do the actual calculation for mkillum
9 */
10
11 #include "mkillum.h"
12
13 #include "face.h"
14
15 #include "cone.h"
16
17 #include "random.h"
18
19
20 o_default(ob, il, rt, nm) /* default illum action */
21 OBJREC *ob;
22 struct illum_args *il;
23 struct rtproc *rt;
24 char *nm;
25 {
26 sprintf(errmsg, "(%s): cannot make illum for %s \"%s\"",
27 nm, ofun[ob->otype].funame, ob->oname);
28 error(WARNING, errmsg);
29 if (!(il->flags & IL_LIGHT))
30 printobj(il->altmat, ob);
31 }
32
33
34 o_face(ob, il, rt, nm) /* make an illum face */
35 OBJREC *ob;
36 struct illum_args *il;
37 struct rtproc *rt;
38 char *nm;
39 {
40 #define MAXMISS (5*n*il->nsamps)
41 int dim[3];
42 int n, nalt, nazi, h;
43 float *distarr;
44 double sp[2], r1, r2;
45 FVECT dn, org, dir;
46 FVECT u, v;
47 double ur[2], vr[2];
48 int nmisses;
49 register FACE *fa;
50 register int i, j;
51 /* get/check arguments */
52 fa = getface(ob);
53 if (fa->area == 0.0) {
54 freeface(ob);
55 o_default(ob, il, rt, nm);
56 return;
57 }
58 /* set up sampling */
59 n = PI * il->sampdens;
60 nalt = sqrt(n/PI) + .5;
61 nazi = PI*nalt + .5;
62 n = nalt*nazi;
63 distarr = (float *)calloc(n, 3*sizeof(float));
64 if (distarr == NULL)
65 error(SYSTEM, "out of memory in o_face");
66 mkaxes(u, v, fa->norm);
67 ur[0] = vr[0] = FHUGE;
68 ur[1] = vr[1] = -FHUGE;
69 for (i = 0; i < fa->nv; i++) {
70 r1 = DOT(VERTEX(fa,i),u);
71 if (r1 < ur[0]) ur[0] = r1;
72 if (r1 > ur[1]) ur[1] = r1;
73 r2 = DOT(VERTEX(fa,i),v);
74 if (r2 < vr[0]) vr[0] = r2;
75 if (r2 > vr[1]) vr[1] = r2;
76 }
77 dim[0] = random();
78 /* sample polygon */
79 nmisses = 0;
80 for (dim[1] = 0; dim[1] < nalt; dim[1]++)
81 for (dim[2] = 0; dim[2] < nazi; dim[2]++)
82 for (i = 0; i < il->nsamps; i++) {
83 /* random direction */
84 h = ilhash(dim, 3) + i;
85 peano(sp, 2, urand(h), .02);
86 r1 = (dim[1] + sp[0])/nalt;
87 r2 = (dim[2] + sp[1])/nazi;
88 flatdir(dn, r1, r2);
89 for (j = 0; j < 3; j++)
90 dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*fa->norm[j];
91 /* random location */
92 do {
93 peano(sp, 2, urand(h+nmisses), .01);
94 r1 = ur[0] + (ur[1]-ur[0]) * sp[0];
95 r2 = vr[0] + (vr[1]-vr[0]) * sp[1];
96 for (j = 0; j < 3; j++)
97 org[j] = r1*u[j] + r2*v[j]
98 + fa->offset*fa->norm[j];
99 } while (!inface(org, fa) && nmisses++ < MAXMISS);
100 if (nmisses > MAXMISS) {
101 objerror(ob, WARNING, "bad aspect");
102 rt->nrays = 0;
103 freeface(ob);
104 free((char *)distarr);
105 o_default(ob, il, rt, nm);
106 return;
107 }
108 for (j = 0; j < 3; j++)
109 org[j] += .001*fa->norm[j];
110 /* send sample */
111 raysamp(distarr+3*(dim[1]*nazi+dim[2]), org, dir, rt);
112 }
113 rayflush(rt);
114 /* write out the face w/ distribution */
115 flatout(il, distarr, nalt, nazi, u, v, fa->norm);
116 illumout(il, ob);
117 /* clean up */
118 freeface(ob);
119 free((char *)distarr);
120 #undef MAXMISS
121 }
122
123
124 o_sphere(ob, il, rt, nm) /* make an illum sphere */
125 register OBJREC *ob;
126 struct illum_args *il;
127 struct rtproc *rt;
128 char *nm;
129 {
130 int dim[3];
131 int n, nalt, nazi;
132 float *distarr;
133 double sp[4], r1, r2, r3;
134 FVECT org, dir;
135 FVECT u, v;
136 register int i, j;
137 /* check arguments */
138 if (ob->oargs.nfargs != 4)
139 objerror(ob, USER, "bad # of arguments");
140 /* set up sampling */
141 n = 4.*PI * il->sampdens;
142 nalt = sqrt(n/PI) + .5;
143 nazi = PI*nalt + .5;
144 n = nalt*nazi;
145 distarr = (float *)calloc(n, 3*sizeof(float));
146 if (distarr == NULL)
147 error(SYSTEM, "out of memory in o_sphere");
148 dim[0] = random();
149 /* sample sphere */
150 for (dim[1] = 0; dim[1] < nalt; dim[1]++)
151 for (dim[2] = 0; dim[2] < nazi; dim[2]++)
152 for (i = 0; i < il->nsamps; i++) {
153 /* next sample point */
154 peano(sp, 4, urand(ilhash(dim,3)+i), .02);
155 /* random direction */
156 r1 = (dim[1] + sp[0])/nalt;
157 r2 = (dim[2] + sp[1])/nazi;
158 rounddir(dir, r1, r2);
159 /* random location */
160 mkaxes(u, v, dir); /* yuck! */
161 r3 = sqrt(sp[2]);
162 r2 = 2.*PI*sp[3];
163 r1 = r3*ob->oargs.farg[3]*cos(r2);
164 r2 = r3*ob->oargs.farg[3]*sin(r2);
165 r3 = ob->oargs.farg[3]*sqrt(1.01-r3*r3);
166 for (j = 0; j < 3; j++) {
167 org[j] = ob->oargs.farg[j] + r1*u[j] + r2*v[j] +
168 r3*dir[j];
169 dir[j] = -dir[j];
170 }
171 /* send sample */
172 raysamp(distarr+3*(dim[1]*nazi+dim[2]), org, dir, rt);
173 }
174 rayflush(rt);
175 /* write out the sphere w/ distribution */
176 roundout(il, distarr, nalt, nazi);
177 illumout(il, ob);
178 /* clean up */
179 free((char *)distarr);
180 }
181
182
183 o_ring(ob, il, rt, nm) /* make an illum ring */
184 OBJREC *ob;
185 struct illum_args *il;
186 struct rtproc *rt;
187 char *nm;
188 {
189 int dim[3];
190 int n, nalt, nazi;
191 float *distarr;
192 double sp[4], r1, r2, r3;
193 FVECT dn, org, dir;
194 FVECT u, v;
195 register CONE *co;
196 register int i, j;
197 /* get/check arguments */
198 co = getcone(ob, 0);
199 /* set up sampling */
200 n = PI * il->sampdens;
201 nalt = sqrt(n/PI) + .5;
202 nazi = PI*nalt + .5;
203 n = nalt*nazi;
204 distarr = (float *)calloc(n, 3*sizeof(float));
205 if (distarr == NULL)
206 error(SYSTEM, "out of memory in o_ring");
207 mkaxes(u, v, co->ad);
208 dim[0] = random();
209 /* sample disk */
210 for (dim[1] = 0; dim[1] < nalt; dim[1]++)
211 for (dim[2] = 0; dim[2] < nazi; dim[2]++)
212 for (i = 0; i < il->nsamps; i++) {
213 /* next sample point */
214 peano(sp, 4, urand(ilhash(dim,3)+i), .02);
215 /* random direction */
216 r1 = (dim[1] + sp[0])/nalt;
217 r2 = (dim[2] + sp[1])/nalt;
218 flatdir(dn, r1, r2);
219 for (j = 0; j < 3; j++)
220 dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*co->ad[j];
221 /* random location */
222 r3 = sqrt(CO_R0(co)*CO_R0(co) +
223 sp[2]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
224 r2 = 2.*PI*sp[3];
225 r1 = r3*cos(r2);
226 r2 = r3*sin(r2);
227 for (j = 0; j < 3; j++)
228 org[j] = CO_P0(co)[j] + r1*u[j] + r1*v[j] +
229 .001*co->ad[j];
230
231 /* send sample */
232 raysamp(distarr+3*(dim[1]*nazi+dim[2]), org, dir, rt);
233 }
234 rayflush(rt);
235 /* write out the ring w/ distribution */
236 flatout(il, distarr, nalt, nazi, u, v, co->ad);
237 illumout(il, ob);
238 /* clean up */
239 freecone(ob);
240 free((char *)distarr);
241 }
242
243
244 raysamp(res, org, dir, rt) /* compute a ray sample */
245 float res[3];
246 FVECT org, dir;
247 register struct rtproc *rt;
248 {
249 register float *fp;
250
251 if (rt->nrays == rt->bsiz)
252 rayflush(rt);
253 rt->dest[rt->nrays] = res;
254 fp = rt->buf + 6*rt->nrays++;
255 *fp++ = org[0]; *fp++ = org[1]; *fp++ = org[2];
256 *fp++ = dir[0]; *fp++ = dir[1]; *fp = dir[2];
257 }
258
259
260 rayflush(rt) /* flush buffered rays */
261 register struct rtproc *rt;
262 {
263 register int i;
264
265 if (rt->nrays <= 0)
266 return;
267 bzero(rt->buf+6*rt->nrays, 6*sizeof(float));
268 if ( process(rt->pd, (char *)rt->buf, (char *)rt->buf,
269 3*sizeof(float)*rt->nrays,
270 6*sizeof(float)*(rt->nrays+1)) <
271 3*sizeof(float)*rt->nrays )
272 error(SYSTEM, "error reading from rtrace process");
273 i = rt->nrays;
274 while (i--) {
275 rt->dest[i][0] += rt->buf[3*i];
276 rt->dest[i][1] += rt->buf[3*i+1];
277 rt->dest[i][2] += rt->buf[3*i+2];
278 }
279 rt->nrays = 0;
280 }
281
282
283 mkaxes(u, v, n) /* compute u and v to go with n */
284 FVECT u, v, n;
285 {
286 register int i;
287
288 v[0] = v[1] = v[2] = 0.0;
289 for (i = 0; i < 3; i++)
290 if (n[i] < 0.6 && n[i] > -0.6)
291 break;
292 v[i] = 1.0;
293 fcross(u, v, n);
294 normalize(u);
295 fcross(v, n, u);
296 }
297
298
299 rounddir(dv, alt, azi) /* compute uniform spherical direction */
300 register FVECT dv;
301 double alt, azi;
302 {
303 double d1, d2;
304
305 dv[2] = 1. - 2.*alt;
306 d1 = sqrt(1. - dv[2]*dv[2]);
307 d2 = 2.*PI * azi;
308 dv[0] = d1*cos(d2);
309 dv[1] = d1*sin(d2);
310 }
311
312
313 flatdir(dv, alt, azi) /* compute uniform hemispherical direction */
314 register FVECT dv;
315 double alt, azi;
316 {
317 double d1, d2;
318
319 d1 = sqrt(alt);
320 d2 = 2.*PI * azi;
321 dv[0] = d1*cos(d2);
322 dv[1] = d1*sin(d2);
323 dv[2] = sqrt(1. - alt);
324 }