ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/mkillum2.c
Revision: 2.7
Committed: Tue Nov 7 17:50:44 1995 UTC (28 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.6: +3 -3 lines
Log Message:
fixed aspect on spherical divisions

File Contents

# Content
1 /* Copyright (c) 1995 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 printobj(il->altmat, ob);
30 }
31
32
33 o_face(ob, il, rt, nm) /* make an illum face */
34 OBJREC *ob;
35 struct illum_args *il;
36 struct rtproc *rt;
37 char *nm;
38 {
39 #define MAXMISS (5*n*il->nsamps)
40 int dim[3];
41 int n, nalt, nazi, h;
42 float *distarr;
43 double sp[2], r1, r2;
44 FVECT dn, org, dir;
45 FVECT u, v;
46 double ur[2], vr[2];
47 int nmisses;
48 register FACE *fa;
49 register int i, j;
50 /* get/check arguments */
51 fa = getface(ob);
52 if (fa->area == 0.0) {
53 freeface(ob);
54 o_default(ob, il, rt, nm);
55 return;
56 }
57 /* set up sampling */
58 if (il->sampdens <= 0)
59 nalt = nazi = 1;
60 else {
61 n = PI * il->sampdens;
62 nalt = sqrt(n/PI) + .5;
63 nazi = PI*nalt + .5;
64 }
65 n = nalt*nazi;
66 distarr = (float *)calloc(n, 3*sizeof(float));
67 if (distarr == NULL)
68 error(SYSTEM, "out of memory in o_face");
69 /* take first edge longer than sqrt(area) */
70 for (j = fa->nv-1, i = 0; i < fa->nv; j = i++) {
71 u[0] = VERTEX(fa,i)[0] - VERTEX(fa,j)[0];
72 u[1] = VERTEX(fa,i)[1] - VERTEX(fa,j)[1];
73 u[2] = VERTEX(fa,i)[2] - VERTEX(fa,j)[2];
74 if ((r1 = DOT(u,u)) >= fa->area-FTINY)
75 break;
76 }
77 if (i < fa->nv) { /* got one! -- let's align our axes */
78 r2 = 1.0/sqrt(r1);
79 u[0] *= r2; u[1] *= r2; u[2] *= r2;
80 fcross(v, fa->norm, u);
81 } else /* oh well, we'll just have to wing it */
82 mkaxes(u, v, fa->norm);
83 /* now, find limits in (u,v) coordinates */
84 ur[0] = vr[0] = FHUGE;
85 ur[1] = vr[1] = -FHUGE;
86 for (i = 0; i < fa->nv; i++) {
87 r1 = DOT(VERTEX(fa,i),u);
88 if (r1 < ur[0]) ur[0] = r1;
89 if (r1 > ur[1]) ur[1] = r1;
90 r2 = DOT(VERTEX(fa,i),v);
91 if (r2 < vr[0]) vr[0] = r2;
92 if (r2 > vr[1]) vr[1] = r2;
93 }
94 dim[0] = random();
95 /* sample polygon */
96 nmisses = 0;
97 for (dim[1] = 0; dim[1] < nalt; dim[1]++)
98 for (dim[2] = 0; dim[2] < nazi; dim[2]++)
99 for (i = 0; i < il->nsamps; i++) {
100 /* random direction */
101 h = ilhash(dim, 3) + i;
102 multisamp(sp, 2, urand(h));
103 r1 = (dim[1] + sp[0])/nalt;
104 r2 = (dim[2] + sp[1] - .5)/nazi;
105 flatdir(dn, r1, r2);
106 for (j = 0; j < 3; j++)
107 dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*fa->norm[j];
108 /* random location */
109 do {
110 multisamp(sp, 2, urand(h+4862+nmisses));
111 r1 = ur[0] + (ur[1]-ur[0]) * sp[0];
112 r2 = vr[0] + (vr[1]-vr[0]) * sp[1];
113 for (j = 0; j < 3; j++)
114 org[j] = r1*u[j] + r2*v[j]
115 + fa->offset*fa->norm[j];
116 } while (!inface(org, fa) && nmisses++ < MAXMISS);
117 if (nmisses > MAXMISS) {
118 objerror(ob, WARNING, "bad aspect");
119 rt->nrays = 0;
120 freeface(ob);
121 free((char *)distarr);
122 o_default(ob, il, rt, nm);
123 return;
124 }
125 for (j = 0; j < 3; j++)
126 org[j] += .001*fa->norm[j];
127 /* send sample */
128 raysamp(distarr+3*(dim[1]*nazi+dim[2]), org, dir, rt);
129 }
130 rayflush(rt);
131 /* write out the face and its distribution */
132 if (average(il, distarr, nalt*nazi)) {
133 if (il->sampdens > 0)
134 flatout(il, distarr, nalt, nazi, u, v, fa->norm);
135 illumout(il, ob);
136 } else
137 printobj(il->altmat, ob);
138 /* clean up */
139 freeface(ob);
140 free((char *)distarr);
141 #undef MAXMISS
142 }
143
144
145 o_sphere(ob, il, rt, nm) /* make an illum sphere */
146 register OBJREC *ob;
147 struct illum_args *il;
148 struct rtproc *rt;
149 char *nm;
150 {
151 int dim[3];
152 int n, nalt, nazi;
153 float *distarr;
154 double sp[4], r1, r2, r3;
155 FVECT org, dir;
156 FVECT u, v;
157 register int i, j;
158 /* check arguments */
159 if (ob->oargs.nfargs != 4)
160 objerror(ob, USER, "bad # of arguments");
161 /* set up sampling */
162 if (il->sampdens <= 0)
163 nalt = nazi = 1;
164 else {
165 n = 4.*PI * il->sampdens;
166 nalt = sqrt(2./PI*n) + .5;
167 nazi = PI/2.*nalt + .5;
168 }
169 n = nalt*nazi;
170 distarr = (float *)calloc(n, 3*sizeof(float));
171 if (distarr == NULL)
172 error(SYSTEM, "out of memory in o_sphere");
173 dim[0] = random();
174 /* sample sphere */
175 for (dim[1] = 0; dim[1] < nalt; dim[1]++)
176 for (dim[2] = 0; dim[2] < nazi; dim[2]++)
177 for (i = 0; i < il->nsamps; i++) {
178 /* next sample point */
179 multisamp(sp, 4, urand(ilhash(dim,3)+i));
180 /* random direction */
181 r1 = (dim[1] + sp[0])/nalt;
182 r2 = (dim[2] + sp[1] - .5)/nazi;
183 rounddir(dir, r1, r2);
184 /* random location */
185 mkaxes(u, v, dir); /* yuck! */
186 r3 = sqrt(sp[2]);
187 r2 = 2.*PI*sp[3];
188 r1 = r3*ob->oargs.farg[3]*cos(r2);
189 r2 = r3*ob->oargs.farg[3]*sin(r2);
190 r3 = ob->oargs.farg[3]*sqrt(1.01-r3*r3);
191 for (j = 0; j < 3; j++) {
192 org[j] = ob->oargs.farg[j] + r1*u[j] + r2*v[j] +
193 r3*dir[j];
194 dir[j] = -dir[j];
195 }
196 /* send sample */
197 raysamp(distarr+3*(dim[1]*nazi+dim[2]), org, dir, rt);
198 }
199 rayflush(rt);
200 /* write out the sphere and its distribution */
201 if (average(il, distarr, nalt*nazi)) {
202 if (il->sampdens > 0)
203 roundout(il, distarr, nalt, nazi);
204 else
205 objerror(ob, WARNING, "diffuse distribution");
206 illumout(il, ob);
207 } else
208 printobj(il->altmat, ob);
209 /* clean up */
210 free((char *)distarr);
211 }
212
213
214 o_ring(ob, il, rt, nm) /* make an illum ring */
215 OBJREC *ob;
216 struct illum_args *il;
217 struct rtproc *rt;
218 char *nm;
219 {
220 int dim[3];
221 int n, nalt, nazi;
222 float *distarr;
223 double sp[4], r1, r2, r3;
224 FVECT dn, org, dir;
225 FVECT u, v;
226 register CONE *co;
227 register int i, j;
228 /* get/check arguments */
229 co = getcone(ob, 0);
230 /* set up sampling */
231 if (il->sampdens <= 0)
232 nalt = nazi = 1;
233 else {
234 n = PI * il->sampdens;
235 nalt = sqrt(n/PI) + .5;
236 nazi = PI*nalt + .5;
237 }
238 n = nalt*nazi;
239 distarr = (float *)calloc(n, 3*sizeof(float));
240 if (distarr == NULL)
241 error(SYSTEM, "out of memory in o_ring");
242 mkaxes(u, v, co->ad);
243 dim[0] = random();
244 /* sample disk */
245 for (dim[1] = 0; dim[1] < nalt; dim[1]++)
246 for (dim[2] = 0; dim[2] < nazi; dim[2]++)
247 for (i = 0; i < il->nsamps; i++) {
248 /* next sample point */
249 multisamp(sp, 4, urand(ilhash(dim,3)+i));
250 /* random direction */
251 r1 = (dim[1] + sp[0])/nalt;
252 r2 = (dim[2] + sp[1] - .5)/nazi;
253 flatdir(dn, r1, r2);
254 for (j = 0; j < 3; j++)
255 dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*co->ad[j];
256 /* random location */
257 r3 = sqrt(CO_R0(co)*CO_R0(co) +
258 sp[2]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
259 r2 = 2.*PI*sp[3];
260 r1 = r3*cos(r2);
261 r2 = r3*sin(r2);
262 for (j = 0; j < 3; j++)
263 org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j] +
264 .001*co->ad[j];
265
266 /* send sample */
267 raysamp(distarr+3*(dim[1]*nazi+dim[2]), org, dir, rt);
268 }
269 rayflush(rt);
270 /* write out the ring and its distribution */
271 if (average(il, distarr, nalt*nazi)) {
272 if (il->sampdens > 0)
273 flatout(il, distarr, nalt, nazi, u, v, co->ad);
274 illumout(il, ob);
275 } else
276 printobj(il->altmat, ob);
277 /* clean up */
278 freecone(ob);
279 free((char *)distarr);
280 }
281
282
283 raysamp(res, org, dir, rt) /* compute a ray sample */
284 float res[3];
285 FVECT org, dir;
286 register struct rtproc *rt;
287 {
288 register float *fp;
289
290 if (rt->nrays == rt->bsiz)
291 rayflush(rt);
292 rt->dest[rt->nrays] = res;
293 fp = rt->buf + 6*rt->nrays++;
294 *fp++ = org[0]; *fp++ = org[1]; *fp++ = org[2];
295 *fp++ = dir[0]; *fp++ = dir[1]; *fp = dir[2];
296 }
297
298
299 rayflush(rt) /* flush buffered rays */
300 register struct rtproc *rt;
301 {
302 register int i;
303
304 if (rt->nrays <= 0)
305 return;
306 bzero(rt->buf+6*rt->nrays, 6*sizeof(float));
307 errno = 0;
308 if ( process(rt->pd, (char *)rt->buf, (char *)rt->buf,
309 3*sizeof(float)*rt->nrays,
310 6*sizeof(float)*(rt->nrays+1)) <
311 3*sizeof(float)*rt->nrays )
312 error(SYSTEM, "error reading from rtrace process");
313 i = rt->nrays;
314 while (i--) {
315 rt->dest[i][0] += rt->buf[3*i];
316 rt->dest[i][1] += rt->buf[3*i+1];
317 rt->dest[i][2] += rt->buf[3*i+2];
318 }
319 rt->nrays = 0;
320 }
321
322
323 mkaxes(u, v, n) /* compute u and v to go with n */
324 FVECT u, v, n;
325 {
326 register int i;
327
328 v[0] = v[1] = v[2] = 0.0;
329 for (i = 0; i < 3; i++)
330 if (n[i] < 0.6 && n[i] > -0.6)
331 break;
332 v[i] = 1.0;
333 fcross(u, v, n);
334 normalize(u);
335 fcross(v, n, u);
336 }
337
338
339 rounddir(dv, alt, azi) /* compute uniform spherical direction */
340 register FVECT dv;
341 double alt, azi;
342 {
343 double d1, d2;
344
345 dv[2] = 1. - 2.*alt;
346 d1 = sqrt(1. - dv[2]*dv[2]);
347 d2 = 2.*PI * azi;
348 dv[0] = d1*cos(d2);
349 dv[1] = d1*sin(d2);
350 }
351
352
353 flatdir(dv, alt, azi) /* compute uniform hemispherical direction */
354 register FVECT dv;
355 double alt, azi;
356 {
357 double d1, d2;
358
359 d1 = sqrt(alt);
360 d2 = 2.*PI * azi;
361 dv[0] = d1*cos(d2);
362 dv[1] = d1*sin(d2);
363 dv[2] = sqrt(1. - alt);
364 }