ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pcond4.c
Revision: 3.12
Committed: Fri Jan 10 13:06:54 1997 UTC (27 years, 3 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.11: +3 -3 lines
Log Message:
added cosine correction to veiling glare computation

File Contents

# User Rev Content
1 greg 3.1 /* Copyright (c) 1996 Regents of the University of California */
2    
3     #ifndef lint
4     static char SCCSid[] = "$SunId$ LBL";
5     #endif
6    
7     /*
8     * Routines for veiling glare and loss of acuity.
9     */
10    
11     #include "pcond.h"
12    
13 greg 3.3 /************** VEILING STUFF *****************/
14 greg 3.1
15     #define VADAPT 0.08 /* fraction of adaptation from veil */
16    
17 greg 3.11 static COLOR *veilimg = NULL; /* veiling image */
18 greg 3.1
19     #define veilscan(y) (veilimg+(y)*fvxr)
20    
21 greg 3.2 static float (*raydir)[3] = NULL; /* ray direction for each pixel */
22 greg 3.1
23     #define rdirscan(y) (raydir+(y)*fvxr)
24    
25    
26     compraydir() /* compute ray directions */
27     {
28 greg 3.2 FVECT rorg, rdir;
29 greg 3.1 double h, v;
30     register int x, y;
31    
32     if (raydir != NULL) /* already done? */
33     return;
34 greg 3.2 raydir = (float (*)[3])malloc(fvxr*fvyr*3*sizeof(float));
35 greg 3.1 if (raydir == NULL)
36     syserror("malloc");
37    
38     for (y = 0; y < fvyr; y++) {
39     switch (inpres.or) {
40     case YMAJOR: case YMAJOR|XDECR:
41     v = (y+.5)/fvyr; break;
42     case YMAJOR|YDECR: case YMAJOR|YDECR|XDECR:
43     v = 1. - (y+.5)/fvyr; break;
44     case 0: case YDECR:
45     h = (y+.5)/fvyr; break;
46     case XDECR: case XDECR|YDECR:
47     h = 1. - (y+.5)/fvyr; break;
48     }
49     for (x = 0; x < fvxr; x++) {
50     switch (inpres.or) {
51     case YMAJOR: case YMAJOR|YDECR:
52     h = (x+.5)/fvxr; break;
53     case YMAJOR|XDECR: case YMAJOR|XDECR|YDECR:
54     h = 1. - (x+.5)/fvxr; break;
55     case 0: case XDECR:
56     v = (x+.5)/fvxr; break;
57     case YDECR: case YDECR|XDECR:
58     v = 1. - (x+.5)/fvxr; break;
59     }
60 greg 3.2 if (viewray(rorg, rdir, &ourview, h, v)
61     >= -FTINY) {
62     rdirscan(y)[x][0] = rdir[0];
63     rdirscan(y)[x][1] = rdir[1];
64     rdirscan(y)[x][2] = rdir[2];
65     } else {
66 greg 3.1 rdirscan(y)[x][0] =
67     rdirscan(y)[x][1] =
68     rdirscan(y)[x][2] = 0.0;
69     }
70     }
71     }
72     }
73    
74    
75     compveil() /* compute veiling image */
76     {
77     double t2, t2sum;
78     COLOR ctmp, vsum;
79     int px, py;
80     register int x, y;
81 greg 3.11
82     if (veilimg != NULL) /* already done? */
83     return;
84 greg 3.1 /* compute ray directions */
85     compraydir();
86     /* compute veil image */
87     veilimg = (COLOR *)malloc(fvxr*fvyr*sizeof(COLOR));
88     if (veilimg == NULL)
89     syserror("malloc");
90     for (py = 0; py < fvyr; py++)
91     for (px = 0; px < fvxr; px++) {
92     t2sum = 0.;
93     setcolor(vsum, 0., 0., 0.);
94     for (y = 0; y < fvyr; y++)
95     for (x = 0; x < fvxr; x++) {
96     if (x == px && y == py) continue;
97     t2 = DOT(rdirscan(py)[px],
98     rdirscan(y)[x]);
99     if (t2 <= FTINY) continue;
100 greg 3.4 /* use approximation instead
101 greg 3.12 t3 = acos(t2);
102     t2 = t2/(t3*t3);
103 greg 3.4 */
104 greg 3.12 t2 *= .5 / (1. - t2);
105 greg 3.1 copycolor(ctmp, fovscan(y)[x]);
106     scalecolor(ctmp, t2);
107     addcolor(vsum, ctmp);
108     t2sum += t2;
109     }
110     /* VADAPT of original is subtracted in addveil() */
111     scalecolor(vsum, VADAPT/t2sum);
112     copycolor(veilscan(py)[px], vsum);
113     }
114 greg 3.11 /* modify FOV sample image */
115     for (y = 0; y < fvyr; y++)
116     for (x = 0; x < fvxr; x++) {
117     scalecolor(fovscan(y)[x], 1.-VADAPT);
118     addcolor(fovscan(y)[x], veilscan(y)[x]);
119     }
120     comphist(); /* recompute histogram */
121 greg 3.1 }
122    
123    
124     addveil(sl, y) /* add veil to scanline */
125     COLOR *sl;
126     int y;
127     {
128     int vx, vy;
129     double dx, dy;
130     double lv, uv;
131     register int x, i;
132    
133     vy = dy = (y+.5)/numscans(&inpres)*fvyr - .5;
134 greg 3.7 while (vy >= fvyr-1) vy--;
135 greg 3.1 dy -= (double)vy;
136     for (x = 0; x < scanlen(&inpres); x++) {
137     vx = dx = (x+.5)/scanlen(&inpres)*fvxr - .5;
138 greg 3.7 while (vx >= fvxr-1) vx--;
139 greg 3.1 dx -= (double)vx;
140     for (i = 0; i < 3; i++) {
141     lv = (1.-dy)*colval(veilscan(vy)[vx],i) +
142     dy*colval(veilscan(vy+1)[vx],i);
143     uv = (1.-dy)*colval(veilscan(vy)[vx+1],i) +
144     dy*colval(veilscan(vy+1)[vx+1],i);
145     colval(sl[x],i) = (1.-VADAPT)*colval(sl[x],i) +
146     (1.-dx)*lv + dx*uv;
147     }
148     }
149 greg 3.3 }
150    
151    
152     /****************** ACUITY STUFF *******************/
153    
154 greg 3.8 typedef struct {
155     short sampe; /* sample area size (exponent of 2) */
156 greg 3.3 short nscans; /* number of scanlines in this bar */
157     int len; /* individual scanline length */
158     int nread; /* number of scanlines loaded */
159 greg 3.8 COLOR *sdata; /* scanbar data */
160 greg 3.3 } SCANBAR;
161    
162 greg 3.8 #define bscan(sb,y) ((COLOR *)(sb)->sdata+((y)%(sb)->nscans)*(sb)->len)
163 greg 3.3
164     SCANBAR *rootbar; /* root scan bar (lowest resolution) */
165    
166     float *inpacuD; /* input acuity data (cycles/degree) */
167    
168     #define tsampr(x,y) inpacuD[(y)*fvxr+(x)]
169    
170    
171     double
172     hacuity(La) /* return visual acuity in cycles/degree */
173     double La;
174     { /* data due to S. Shaler (we should fit it!) */
175     #define NPOINTS 20
176     static float l10lum[NPOINTS] = {
177     -3.10503,-2.66403,-2.37703,-2.09303,-1.64403,-1.35803,
178     -1.07403,-0.67203,-0.38503,-0.10103,0.29397,0.58097,0.86497,
179     1.25697,1.54397,1.82797,2.27597,2.56297,2.84697,3.24897
180     };
181     static float resfreq[NPOINTS] = {
182     2.09,3.28,3.79,4.39,6.11,8.83,10.94,18.66,23.88,31.05,37.42,
183     37.68,41.60,43.16,45.30,47.00,48.43,48.32,51.06,51.09
184     };
185     double l10La;
186     register int i;
187 greg 3.5 /* check limits */
188     if (La <= 7.85e-4)
189     return(resfreq[0]);
190     if (La >= 1.78e3)
191     return(resfreq[NPOINTS-1]);
192     /* interpolate data */
193 greg 3.3 l10La = log10(La);
194     for (i = 0; i < NPOINTS-2 && l10lum[i+1] <= l10La; i++)
195     ;
196     return( ( (l10lum[i+1] - l10La)*resfreq[i] +
197     (l10La - l10lum[i])*resfreq[i+1] ) /
198     (l10lum[i+1] - l10lum[i]) );
199     #undef NPOINTS
200     }
201    
202    
203     COLOR *
204     getascan(sb, y) /* find/read scanline y for scanbar sb */
205     register SCANBAR *sb;
206     int y;
207     {
208     register COLOR *sl0, *sl1, *mysl;
209     register int i;
210    
211 greg 3.6 if (y < sb->nread - sb->nscans) /* too far back? */
212     return(NULL);
213 greg 3.3 for ( ; y >= sb->nread; sb->nread++) { /* read as necessary */
214     mysl = bscan(sb, sb->nread);
215 greg 3.8 if (sb->sampe == 0) {
216 greg 3.3 if (freadscan(mysl, sb->len, infp) < 0) {
217     fprintf(stderr, "%s: %s: scanline read error\n",
218     progname, infn);
219     exit(1);
220     }
221     } else {
222 greg 3.8 sl0 = getascan(sb+1, 2*y);
223 greg 3.6 if (sl0 == NULL)
224     return(NULL);
225 greg 3.8 sl1 = getascan(sb+1, 2*y+1);
226 greg 3.3 for (i = 0; i < sb->len; i++) {
227     copycolor(mysl[i], sl0[2*i]);
228     addcolor(mysl[i], sl0[2*i+1]);
229     addcolor(mysl[i], sl1[2*i]);
230     addcolor(mysl[i], sl1[2*i+1]);
231     scalecolor(mysl[i], 0.25);
232     }
233     }
234     }
235     return(bscan(sb, y));
236     }
237    
238    
239     acuscan(scln, y) /* get acuity-sampled scanline */
240     COLOR *scln;
241     int y;
242     {
243     double sr;
244     double dx, dy;
245     int ix, iy;
246     register int x;
247     /* compute foveal y position */
248     iy = dy = (y+.5)/numscans(&inpres)*fvyr - .5;
249 greg 3.9 while (iy >= fvyr-1) iy--;
250 greg 3.3 dy -= (double)iy;
251     for (x = 0; x < scanlen(&inpres); x++) {
252     /* compute foveal x position */
253     ix = dx = (x+.5)/scanlen(&inpres)*fvxr - .5;
254 greg 3.9 while (ix >= fvxr-1) ix--;
255 greg 3.3 dx -= (double)ix;
256     /* interpolate sample rate */
257     sr = (1.-dy)*((1.-dx)*tsampr(ix,iy) + dx*tsampr(ix+1,iy)) +
258     dy*((1.-dx)*tsampr(ix,iy+1) + dx*tsampr(ix+1,iy+1));
259    
260     acusample(scln[x], x, y, sr); /* compute sample */
261     }
262     }
263    
264    
265     acusample(col, x, y, sr) /* interpolate sample at (x,y) using rate sr */
266     COLOR col;
267     int x, y;
268     double sr;
269     {
270     COLOR c1;
271     double d;
272     register SCANBAR *sb0;
273    
274 greg 3.8 for (sb0 = rootbar; sb0->sampe != 0 && 1<<sb0[1].sampe > sr; sb0++)
275 greg 3.3 ;
276     ascanval(col, x, y, sb0);
277 greg 3.8 if (sb0->sampe == 0) /* don't extrapolate highest */
278 greg 3.3 return;
279 greg 3.8 ascanval(c1, x, y, sb0+1);
280     d = ((1<<sb0->sampe) - sr)/(1<<sb0[1].sampe);
281 greg 3.3 scalecolor(col, 1.-d);
282     scalecolor(c1, d);
283     addcolor(col, c1);
284     }
285    
286    
287     ascanval(col, x, y, sb) /* interpolate scanbar at orig. coords (x,y) */
288     COLOR col;
289     int x, y;
290     SCANBAR *sb;
291     {
292     COLOR *sl0, *sl1, c1, c1y;
293     double dx, dy;
294     int ix, iy;
295    
296 greg 3.8 if (sb->sampe == 0) { /* no need to interpolate */
297 greg 3.6 sl0 = getascan(sb, y);
298     copycolor(col, sl0[x]);
299     return;
300     }
301     /* compute coordinates for sb */
302 greg 3.8 ix = dx = (x+.5)/(1<<sb->sampe) - .5;
303 greg 3.7 while (ix >= sb->len-1) ix--;
304 greg 3.3 dx -= (double)ix;
305 greg 3.8 iy = dy = (y+.5)/(1<<sb->sampe) - .5;
306     while (iy >= (numscans(&inpres)>>sb->sampe)-1) iy--;
307 greg 3.3 dy -= (double)iy;
308     /* get scanlines */
309     sl0 = getascan(sb, iy);
310 greg 3.9 #ifdef DEBUG
311 greg 3.6 if (sl0 == NULL) {
312     fprintf(stderr, "%s: internal - cannot backspace in ascanval\n",
313     progname);
314 greg 3.9 abort();
315 greg 3.6 }
316 greg 3.9 #endif
317 greg 3.3 sl1 = getascan(sb, iy+1);
318     /* 2D linear interpolation */
319     copycolor(col, sl0[ix]);
320     scalecolor(col, 1.-dx);
321     copycolor(c1, sl0[ix+1]);
322     scalecolor(c1, dx);
323     addcolor(col, c1);
324     copycolor(c1y, sl1[ix]);
325     scalecolor(c1y, 1.-dx);
326     copycolor(c1, sl1[ix+1]);
327     scalecolor(c1, dx);
328     addcolor(c1y, c1);
329     scalecolor(col, 1.-dy);
330     scalecolor(c1y, dy);
331     addcolor(col, c1y);
332 greg 3.9 for (ix = 0; ix < 3; ix++) /* make sure no negative */
333     if (colval(col,ix) < 0.)
334     colval(col,ix) = 0.;
335 greg 3.3 }
336    
337    
338     SCANBAR *
339 greg 3.9 sballoc(se, ns, sl) /* allocate scanbar */
340     int se; /* sampling rate exponent */
341 greg 3.3 int ns; /* number of scanlines */
342     int sl; /* original scanline length */
343     {
344 greg 3.8 SCANBAR *sbarr;
345 greg 3.3 register SCANBAR *sb;
346    
347 greg 3.9 sbarr = sb = (SCANBAR *)malloc((se+1)*sizeof(SCANBAR));
348 greg 3.3 if (sb == NULL)
349     syserror("malloc");
350 greg 3.8 do {
351 greg 3.9 sb->sampe = se;
352     sb->len = sl>>se;
353     sb->nscans = ns;
354     sb->sdata = (COLOR *)malloc(sb->len*ns*sizeof(COLOR));
355 greg 3.8 if (sb->sdata == NULL)
356     syserror("malloc");
357     sb->nread = 0;
358     ns <<= 1;
359     sb++;
360 greg 3.9 } while (--se >= 0);
361 greg 3.8 return(sbarr);
362 greg 3.3 }
363    
364    
365     initacuity() /* initialize variable acuity sampling */
366     {
367     FVECT diffx, diffy, cp;
368     double omega, maxsr;
369     register int x, y, i;
370    
371     compraydir(); /* compute ray directions */
372    
373     inpacuD = (float *)malloc(fvxr*fvyr*sizeof(float));
374     if (inpacuD == NULL)
375     syserror("malloc");
376     maxsr = 1.; /* compute internal sample rates */
377     for (y = 1; y < fvyr-1; y++)
378     for (x = 1; x < fvxr-1; x++) {
379     for (i = 0; i < 3; i++) {
380     diffx[i] = 0.5*fvxr/scanlen(&inpres) *
381     (rdirscan(y)[x+1][i] -
382     rdirscan(y)[x-1][i]);
383     diffy[i] = 0.5*fvyr/numscans(&inpres) *
384     (rdirscan(y+1)[x][i] -
385     rdirscan(y-1)[x][i]);
386     }
387     fcross(cp, diffx, diffy);
388     omega = 0.5 * sqrt(DOT(cp,cp));
389 greg 3.10 if (omega <= FTINY*FTINY)
390 greg 3.4 tsampr(x,y) = 1.;
391     else if ((tsampr(x,y) = PI/180. / sqrt(omega) /
392     hacuity(plum(fovscan(y)[x]))) > maxsr)
393 greg 3.3 maxsr = tsampr(x,y);
394     }
395     /* copy perimeter (easier) */
396     for (x = 1; x < fvxr-1; x++) {
397     tsampr(x,0) = tsampr(x,1);
398     tsampr(x,fvyr-1) = tsampr(x,fvyr-2);
399     }
400     for (y = 0; y < fvyr; y++) {
401 greg 3.9 tsampr(0,y) = tsampr(1,y);
402     tsampr(fvxr-1,y) = tsampr(fvxr-2,y);
403 greg 3.3 }
404     /* initialize with next power of two */
405 greg 3.8 rootbar = sballoc((int)(log(maxsr)/log(2.))+1, 2, scanlen(&inpres));
406 greg 3.1 }