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

Comparing ray/src/gen/mkillum2.c (file contents):
Revision 1.1 by greg, Tue Jul 23 15:42:42 1991 UTC vs.
Revision 2.16 by greg, Fri Sep 17 21:43:50 2004 UTC

# Line 1 | Line 1
1 /* Copyright (c) 1991 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 < * Routines to do the actual calcultion and output for mkillum
5 > * Routines to do the actual calculation for mkillum
6   */
7  
8 < #include  "mkillum.h"
8 > #include <string.h>
9  
10 + #include  "mkillum.h"
11   #include  "face.h"
14
12   #include  "cone.h"
13 + #include  "random.h"
14 + #include  "selcall.h"
15  
16  
17 < printobj(mod, obj)              /* print out an object */
18 < char  *mod;
19 < register OBJREC  *obj;
17 > static void mkaxes(FVECT u, FVECT v, FVECT n);
18 > static void rounddir(FVECT dv, double alt, double azi);
19 > static void flatdir(FVECT dv, double alt, double azi);
20 >
21 >
22 > static void
23 > rayclean(                       /* finish all pending rays */
24 >        struct rtproc *rt0
25 > )
26   {
27 <        register int  i;
27 >        rayflush(rt0, 1);
28 >        while (raywait(rt0) != NULL)
29 >                ;
30 > }
31  
32 <        printf("\n%s %s %s", mod, ofun[obj->otype].funame, obj->oname);
33 <        printf("\n%d", obj->oargs.nsargs);
34 <        for (i = 0; i < obj->oargs.nsargs; i++)
35 <                printf(" %s", obj->oargs.sarg[i]);
36 < #ifdef  IARGS
37 <        printf("\n%d", obj->oargs.niargs);
38 <        for (i = 0; i < obj->oargs.niargs; i++)
39 <                printf(" %d", obj->oargs.iarg[i]);
40 < #else
41 <        printf("\n0");
42 < #endif
43 <        printf("\n%d", obj->oargs.nfargs);
44 <        for (i = 0; i < obj->oargs.nfargs; i++) {
45 <                if (i%3 == 0)
46 <                        putchar('\n');
47 <                printf(" %18.12g", obj->oargs.farg[i]);
32 >
33 > int /* XXX type conflict with otypes.h */
34 > o_default(      /* default illum action */
35 >        OBJREC  *ob,
36 >        struct illum_args  *il,
37 >        struct rtproc  *rt0,
38 >        char  *nm
39 > )
40 > {
41 >        sprintf(errmsg, "(%s): cannot make illum for %s \"%s\"",
42 >                        nm, ofun[ob->otype].funame, ob->oname);
43 >        error(WARNING, errmsg);
44 >        printobj(il->altmat, ob);
45 >        return(1);
46 > }
47 >
48 >
49 > int
50 > o_face(         /* make an illum face */
51 >        OBJREC  *ob,
52 >        struct illum_args  *il,
53 >        struct rtproc  *rt0,
54 >        char  *nm
55 > )
56 > {
57 > #define MAXMISS         (5*n*il->nsamps)
58 >        int  dim[3];
59 >        int  n, nalt, nazi, h;
60 >        float  *distarr;
61 >        double  sp[2], r1, r2;
62 >        FVECT  dn, org, dir;
63 >        FVECT  u, v;
64 >        double  ur[2], vr[2];
65 >        int  nmisses;
66 >        register FACE  *fa;
67 >        register int  i, j;
68 >                                /* get/check arguments */
69 >        fa = getface(ob);
70 >        if (fa->area == 0.0) {
71 >                freeface(ob);
72 >                return(o_default(ob, il, rt0, nm));
73          }
74 <        putchar('\n');
74 >                                /* set up sampling */
75 >        if (il->sampdens <= 0)
76 >                nalt = nazi = 1;
77 >        else {
78 >                n = PI * il->sampdens;
79 >                nalt = sqrt(n/PI) + .5;
80 >                nazi = PI*nalt + .5;
81 >        }
82 >        n = nalt*nazi;
83 >        distarr = (float *)calloc(n, 3*sizeof(float));
84 >        if (distarr == NULL)
85 >                error(SYSTEM, "out of memory in o_face");
86 >                                /* take first edge longer than sqrt(area) */
87 >        for (j = fa->nv-1, i = 0; i < fa->nv; j = i++) {
88 >                u[0] = VERTEX(fa,i)[0] - VERTEX(fa,j)[0];
89 >                u[1] = VERTEX(fa,i)[1] - VERTEX(fa,j)[1];
90 >                u[2] = VERTEX(fa,i)[2] - VERTEX(fa,j)[2];
91 >                if ((r1 = DOT(u,u)) >= fa->area-FTINY)
92 >                        break;
93 >        }
94 >        if (i < fa->nv) {       /* got one! -- let's align our axes */
95 >                r2 = 1.0/sqrt(r1);
96 >                u[0] *= r2; u[1] *= r2; u[2] *= r2;
97 >                fcross(v, fa->norm, u);
98 >        } else                  /* oh well, we'll just have to wing it */
99 >                mkaxes(u, v, fa->norm);
100 >                                /* now, find limits in (u,v) coordinates */
101 >        ur[0] = vr[0] = FHUGE;
102 >        ur[1] = vr[1] = -FHUGE;
103 >        for (i = 0; i < fa->nv; i++) {
104 >                r1 = DOT(VERTEX(fa,i),u);
105 >                if (r1 < ur[0]) ur[0] = r1;
106 >                if (r1 > ur[1]) ur[1] = r1;
107 >                r2 = DOT(VERTEX(fa,i),v);
108 >                if (r2 < vr[0]) vr[0] = r2;
109 >                if (r2 > vr[1]) vr[1] = r2;
110 >        }
111 >        dim[0] = random();
112 >                                /* sample polygon */
113 >        nmisses = 0;
114 >        for (dim[1] = 0; dim[1] < nalt; dim[1]++)
115 >            for (dim[2] = 0; dim[2] < nazi; dim[2]++)
116 >                for (i = 0; i < il->nsamps; i++) {
117 >                                        /* random direction */
118 >                    h = ilhash(dim, 3) + i;
119 >                    multisamp(sp, 2, urand(h));
120 >                    r1 = (dim[1] + sp[0])/nalt;
121 >                    r2 = (dim[2] + sp[1] - .5)/nazi;
122 >                    flatdir(dn, r1, r2);
123 >                    for (j = 0; j < 3; j++)
124 >                        dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*fa->norm[j];
125 >                                        /* random location */
126 >                    do {
127 >                        multisamp(sp, 2, urand(h+4862+nmisses));
128 >                        r1 = ur[0] + (ur[1]-ur[0]) * sp[0];
129 >                        r2 = vr[0] + (vr[1]-vr[0]) * sp[1];
130 >                        for (j = 0; j < 3; j++)
131 >                            org[j] = r1*u[j] + r2*v[j]
132 >                                        + fa->offset*fa->norm[j];
133 >                    } while (!inface(org, fa) && nmisses++ < MAXMISS);
134 >                    if (nmisses > MAXMISS) {
135 >                        objerror(ob, WARNING, "bad aspect");
136 >                        rayclean(rt0);
137 >                        freeface(ob);
138 >                        free((void *)distarr);
139 >                        return(o_default(ob, il, rt0, nm));
140 >                    }
141 >                    for (j = 0; j < 3; j++)
142 >                        org[j] += .001*fa->norm[j];
143 >                                        /* send sample */
144 >                    raysamp(distarr+3*(dim[1]*nazi+dim[2]), org, dir, rt0);
145 >                }
146 >        rayclean(rt0);
147 >                                /* write out the face and its distribution */
148 >        if (average(il, distarr, nalt*nazi)) {
149 >                if (il->sampdens > 0)
150 >                        flatout(il, distarr, nalt, nazi, u, v, fa->norm);
151 >                illumout(il, ob);
152 >        } else
153 >                printobj(il->altmat, ob);
154 >                                /* clean up */
155 >        freeface(ob);
156 >        free((void *)distarr);
157 >        return(0);
158 > #undef MAXMISS
159   }
160  
161  
162 < mkillum(ob, il, rt)             /* make an illum object */
163 < OBJREC  *ob;
164 < struct illum_args  *il;
165 < struct rtproc  *rt;
162 > int
163 > o_sphere(       /* make an illum sphere */
164 >        register OBJREC  *ob,
165 >        struct illum_args  *il,
166 >        struct rtproc  *rt0,
167 >        char  *nm
168 > )
169   {
170 +        int  dim[3];
171 +        int  n, nalt, nazi;
172 +        float  *distarr;
173 +        double  sp[4], r1, r2, r3;
174 +        FVECT  org, dir;
175 +        FVECT  u, v;
176 +        register int  i, j;
177 +                                /* check arguments */
178 +        if (ob->oargs.nfargs != 4)
179 +                objerror(ob, USER, "bad # of arguments");
180 +                                /* set up sampling */
181 +        if (il->sampdens <= 0)
182 +                nalt = nazi = 1;
183 +        else {
184 +                n = 4.*PI * il->sampdens;
185 +                nalt = sqrt(2./PI*n) + .5;
186 +                nazi = PI/2.*nalt + .5;
187 +        }
188 +        n = nalt*nazi;
189 +        distarr = (float *)calloc(n, 3*sizeof(float));
190 +        if (distarr == NULL)
191 +                error(SYSTEM, "out of memory in o_sphere");
192 +        dim[0] = random();
193 +                                /* sample sphere */
194 +        for (dim[1] = 0; dim[1] < nalt; dim[1]++)
195 +            for (dim[2] = 0; dim[2] < nazi; dim[2]++)
196 +                for (i = 0; i < il->nsamps; i++) {
197 +                                        /* next sample point */
198 +                    multisamp(sp, 4, urand(ilhash(dim,3)+i));
199 +                                        /* random direction */
200 +                    r1 = (dim[1] + sp[0])/nalt;
201 +                    r2 = (dim[2] + sp[1] - .5)/nazi;
202 +                    rounddir(dir, r1, r2);
203 +                                        /* random location */
204 +                    mkaxes(u, v, dir);          /* yuck! */
205 +                    r3 = sqrt(sp[2]);
206 +                    r2 = 2.*PI*sp[3];
207 +                    r1 = r3*ob->oargs.farg[3]*cos(r2);
208 +                    r2 = r3*ob->oargs.farg[3]*sin(r2);
209 +                    r3 = ob->oargs.farg[3]*sqrt(1.01-r3*r3);
210 +                    for (j = 0; j < 3; j++) {
211 +                        org[j] = ob->oargs.farg[j] + r1*u[j] + r2*v[j] +
212 +                                        r3*dir[j];
213 +                        dir[j] = -dir[j];
214 +                    }
215 +                                        /* send sample */
216 +                    raysamp(distarr+3*(dim[1]*nazi+dim[2]), org, dir, rt0);
217 +                }
218 +        rayclean(rt0);
219 +                                /* write out the sphere and its distribution */
220 +        if (average(il, distarr, nalt*nazi)) {
221 +                if (il->sampdens > 0)
222 +                        roundout(il, distarr, nalt, nazi);
223 +                else
224 +                        objerror(ob, WARNING, "diffuse distribution");
225 +                illumout(il, ob);
226 +        } else
227 +                printobj(il->altmat, ob);
228 +                                /* clean up */
229 +        free((void *)distarr);
230 +        return(1);
231 + }
232 +
233 +
234 + int
235 + o_ring(         /* make an illum ring */
236 +        OBJREC  *ob,
237 +        struct illum_args  *il,
238 +        struct rtproc  *rt0,
239 +        char  *nm
240 + )
241 + {
242 +        int  dim[3];
243 +        int  n, nalt, nazi;
244 +        float  *distarr;
245 +        double  sp[4], r1, r2, r3;
246 +        FVECT  dn, org, dir;
247 +        FVECT  u, v;
248 +        register CONE  *co;
249 +        register int  i, j;
250 +                                /* get/check arguments */
251 +        co = getcone(ob, 0);
252 +                                /* set up sampling */
253 +        if (il->sampdens <= 0)
254 +                nalt = nazi = 1;
255 +        else {
256 +                n = PI * il->sampdens;
257 +                nalt = sqrt(n/PI) + .5;
258 +                nazi = PI*nalt + .5;
259 +        }
260 +        n = nalt*nazi;
261 +        distarr = (float *)calloc(n, 3*sizeof(float));
262 +        if (distarr == NULL)
263 +                error(SYSTEM, "out of memory in o_ring");
264 +        mkaxes(u, v, co->ad);
265 +        dim[0] = random();
266 +                                /* sample disk */
267 +        for (dim[1] = 0; dim[1] < nalt; dim[1]++)
268 +            for (dim[2] = 0; dim[2] < nazi; dim[2]++)
269 +                for (i = 0; i < il->nsamps; i++) {
270 +                                        /* next sample point */
271 +                    multisamp(sp, 4, urand(ilhash(dim,3)+i));
272 +                                        /* random direction */
273 +                    r1 = (dim[1] + sp[0])/nalt;
274 +                    r2 = (dim[2] + sp[1] - .5)/nazi;
275 +                    flatdir(dn, r1, r2);
276 +                    for (j = 0; j < 3; j++)
277 +                        dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*co->ad[j];
278 +                                        /* random location */
279 +                    r3 = sqrt(CO_R0(co)*CO_R0(co) +
280 +                            sp[2]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
281 +                    r2 = 2.*PI*sp[3];
282 +                    r1 = r3*cos(r2);
283 +                    r2 = r3*sin(r2);
284 +                    for (j = 0; j < 3; j++)
285 +                        org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j] +
286 +                                        .001*co->ad[j];
287 +
288 +                                        /* send sample */
289 +                    raysamp(distarr+3*(dim[1]*nazi+dim[2]), org, dir, rt0);
290 +                }
291 +        rayclean(rt0);
292 +                                /* write out the ring and its distribution */
293 +        if (average(il, distarr, nalt*nazi)) {
294 +                if (il->sampdens > 0)
295 +                        flatout(il, distarr, nalt, nazi, u, v, co->ad);
296 +                illumout(il, ob);
297 +        } else
298 +                printobj(il->altmat, ob);
299 +                                /* clean up */
300 +        freecone(ob);
301 +        free((void *)distarr);
302 +        return(1);
303 + }
304 +
305 +
306 + void
307 + raysamp(        /* queue a ray sample */
308 +        float  res[3],
309 +        FVECT  org,
310 +        FVECT  dir,
311 +        struct rtproc *rt0
312 + )
313 + {
314 +        register struct rtproc  *rt;
315 +        register float  *fp;
316 +
317 +        for (rt = rt0; rt != NULL; rt = rt->next)
318 +                if (rt->nrays < rt->bsiz && rt->dest[rt->nrays] == NULL)
319 +                        break;
320 +        if (rt == NULL)         /* need to free up buffer? */
321 +                rt = raywait(rt0);
322 +        if (rt == NULL)
323 +                error(SYSTEM, "raywait() returned NULL");
324 +        fp = rt->buf + 6*rt->nrays;
325 +        *fp++ = org[0]; *fp++ = org[1]; *fp++ = org[2];
326 +        *fp++ = dir[0]; *fp++ = dir[1]; *fp = dir[2];
327 +        rt->dest[rt->nrays++] = res;
328 +        if (rt->nrays == rt->bsiz)
329 +                rayflush(rt, 0);
330 + }
331 +
332 +
333 + void
334 + rayflush(                       /* flush queued rays to rtrace */
335 +        register struct rtproc  *rt,
336 +        int  doall
337 + )
338 + {
339 +        int     nw;
340 +
341 +        do {
342 +                if (rt->nrays <= 0)
343 +                        continue;
344 +                memset(rt->buf+6*rt->nrays, 0, 6*sizeof(float));
345 +                nw = 6*sizeof(float)*(rt->nrays+1);
346 +                errno = 0;
347 +                if (writebuf(rt->pd.w, (char *)rt->buf, nw) < nw)
348 +                        error(SYSTEM, "error writing to rtrace process");
349 +                rt->nrays = 0;          /* flag buffer as flushed */
350 +        } while (doall && (rt = rt->next) != NULL);
351 + }
352 +
353 +
354 + struct rtproc *
355 + raywait(                        /* retrieve rtrace results */
356 +        struct rtproc *rt0
357 + )
358 + {
359 +        fd_set  readset, errset;
360 +        int     nr;
361 +        struct rtproc   *rt, *rtfree;
362 +        register int    n;
363 +                                /* prepare select call */
364 +        FD_ZERO(&readset); FD_ZERO(&errset); n = 0;
365 +        nr = 0;
366 +        for (rt = rt0; rt != NULL; rt = rt->next) {
367 +                if (rt->nrays == 0 && rt->dest[0] != NULL) {
368 +                        FD_SET(rt->pd.r, &readset);
369 +                        ++nr;
370 +                }
371 +                FD_SET(rt->pd.r, &errset);
372 +                if (rt->pd.r >= n)
373 +                        n = rt->pd.r + 1;
374 +        }
375 +        if (!nr)                /* no rays pending */
376 +                return(NULL);
377 +        if (nr > 1)             /* call select for multiple processes */
378 +                n = select(n, &readset, (fd_set *)NULL, &errset,
379 +                                        (struct timeval *)NULL);
380 +        else
381 +                FD_ZERO(&errset);
382 +        if (n <= 0)
383 +                return(NULL);
384 +        rtfree = NULL;          /* read from ready process(es) */
385 +        for (rt = rt0; rt != NULL; rt = rt->next) {
386 +                if (!FD_ISSET(rt->pd.r, &readset) &&
387 +                                !FD_ISSET(rt->pd.r, &errset))
388 +                        continue;
389 +                for (n = 0; n < rt->bsiz && rt->dest[n] != NULL; n++)
390 +                        ;
391 +                errno = 0;
392 +                nr = read(rt->pd.r, (char *)rt->buf, 3*sizeof(float)*(n+1));
393 +                if (nr < 0)
394 +                        error(SYSTEM, "read error in raywait()");
395 +                if (nr == 0)                            /* unexpected EOF */
396 +                        error(USER, "rtrace process died");
397 +                if (nr < 3*sizeof(float)*(n+1)) {       /* read the rest */
398 +                        nr = readbuf(rt->pd.r, (char *)rt->buf,
399 +                                        3*sizeof(float)*(n+1) - nr);
400 +                        if (nr < 0)
401 +                                error(USER, "readbuf error in raywait()");
402 +                }
403 +                while (n-- > 0) {
404 +                        rt->dest[n][0] += rt->buf[3*n];
405 +                        rt->dest[n][1] += rt->buf[3*n+1];
406 +                        rt->dest[n][2] += rt->buf[3*n+2];
407 +                        rt->dest[n] = NULL;
408 +                }
409 +                rtfree = rt;
410 +        }
411 +        return(rtfree);
412 + }
413 +
414 +
415 + static void
416 + mkaxes(                 /* compute u and v to go with n */
417 +        FVECT  u,
418 +        FVECT  v,
419 +        FVECT  n
420 + )
421 + {
422 +        register int  i;
423 +
424 +        v[0] = v[1] = v[2] = 0.0;
425 +        for (i = 0; i < 3; i++)
426 +                if (n[i] < 0.6 && n[i] > -0.6)
427 +                        break;
428 +        v[i] = 1.0;
429 +        fcross(u, v, n);
430 +        normalize(u);
431 +        fcross(v, n, u);
432 + }
433 +
434 +
435 + static void
436 + rounddir(               /* compute uniform spherical direction */
437 +        register FVECT  dv,
438 +        double  alt,
439 +        double  azi
440 + )
441 + {
442 +        double  d1, d2;
443 +
444 +        dv[2] = 1. - 2.*alt;
445 +        d1 = sqrt(1. - dv[2]*dv[2]);
446 +        d2 = 2.*PI * azi;
447 +        dv[0] = d1*cos(d2);
448 +        dv[1] = d1*sin(d2);
449 + }
450 +
451 +
452 + static void
453 + flatdir(                /* compute uniform hemispherical direction */
454 +        register FVECT  dv,
455 +        double  alt,
456 +        double  azi
457 + )
458 + {
459 +        double  d1, d2;
460 +
461 +        d1 = sqrt(alt);
462 +        d2 = 2.*PI * azi;
463 +        dv[0] = d1*cos(d2);
464 +        dv[1] = d1*sin(d2);
465 +        dv[2] = sqrt(1. - alt);
466   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines