ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/mkillum2.c
Revision: 2.12
Committed: Sun Nov 16 10:29:38 2003 UTC (20 years, 4 months ago) by schorsch
Content type: text/plain
Branch: MAIN
Changes since 2.11: +69 -36 lines
Log Message:
Continued ANSIfication and reduced other compile warnings.

File Contents

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