ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pcond4.c
Revision: 3.19
Committed: Sun Nov 14 16:57:18 2004 UTC (19 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad3R7P2, rad3R7P1, rad4R2, rad4R1, rad4R0, rad3R8, rad3R9, rad4R2P1
Changes since 3.18: +2 -6 lines
Log Message:
Changed logic so small view angles are prevented along with parallel views

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: pcond4.c,v 3.18 2004/11/08 15:50:59 greg Exp $";
3 #endif
4 /*
5 * Routines for veiling glare and loss of acuity.
6 */
7
8 #include "pcond.h"
9
10 /************** VEILING STUFF *****************/
11
12 #define VADAPT 0.08 /* fraction of adaptation from veil */
13
14 static COLOR *veilimg = NULL; /* veiling image */
15
16 #define veilscan(y) (veilimg+(y)*fvxr)
17
18 static float (*raydir)[3] = NULL; /* ray direction for each pixel */
19
20 #define rdirscan(y) (raydir+(y)*fvxr)
21
22 static void compraydir(void);
23 #if ADJ_VEIL
24 static void smoothveil(void);
25 #endif
26
27
28 static void
29 compraydir(void) /* compute ray directions */
30 {
31 FVECT rorg, rdir;
32 double h, v;
33 register int x, y;
34
35 if (raydir != NULL) /* already done? */
36 return;
37 raydir = (float (*)[3])malloc(fvxr*fvyr*3*sizeof(float));
38 if (raydir == NULL)
39 syserror("malloc");
40
41 for (y = 0; y < fvyr; y++) {
42 switch (inpres.rt) {
43 case YMAJOR: case YMAJOR|XDECR:
44 v = (y+.5)/fvyr; break;
45 case YMAJOR|YDECR: case YMAJOR|YDECR|XDECR:
46 v = 1. - (y+.5)/fvyr; break;
47 case 0: case YDECR:
48 h = (y+.5)/fvyr; break;
49 case XDECR: case XDECR|YDECR:
50 h = 1. - (y+.5)/fvyr; break;
51 }
52 for (x = 0; x < fvxr; x++) {
53 switch (inpres.rt) {
54 case YMAJOR: case YMAJOR|YDECR:
55 h = (x+.5)/fvxr; break;
56 case YMAJOR|XDECR: case YMAJOR|XDECR|YDECR:
57 h = 1. - (x+.5)/fvxr; break;
58 case 0: case XDECR:
59 v = (x+.5)/fvxr; break;
60 case YDECR: case YDECR|XDECR:
61 v = 1. - (x+.5)/fvxr; break;
62 }
63 if (viewray(rorg, rdir, &ourview, h, v)
64 >= -FTINY) {
65 rdirscan(y)[x][0] = rdir[0];
66 rdirscan(y)[x][1] = rdir[1];
67 rdirscan(y)[x][2] = rdir[2];
68 } else {
69 rdirscan(y)[x][0] =
70 rdirscan(y)[x][1] =
71 rdirscan(y)[x][2] = 0.0;
72 }
73 }
74 }
75 }
76
77
78 extern void
79 compveil(void) /* compute veiling image */
80 {
81 double t2, t2sum;
82 COLOR ctmp, vsum;
83 int px, py;
84 register int x, y;
85
86 if (veilimg != NULL) /* already done? */
87 return;
88 /* compute ray directions */
89 compraydir();
90 /* compute veil image */
91 veilimg = (COLOR *)malloc(fvxr*fvyr*sizeof(COLOR));
92 if (veilimg == NULL)
93 syserror("malloc");
94 for (py = 0; py < fvyr; py++)
95 for (px = 0; px < fvxr; px++) {
96 t2sum = 0.;
97 setcolor(vsum, 0., 0., 0.);
98 for (y = 0; y < fvyr; y++)
99 for (x = 0; x < fvxr; x++) {
100 if (x == px && y == py) continue;
101 t2 = DOT(rdirscan(py)[px],
102 rdirscan(y)[x]);
103 if (t2 <= FTINY) continue;
104 /* use approximation instead
105 t3 = acos(t2);
106 t2 = t2/(t3*t3);
107 */
108 t2 *= .5 / (1. - t2);
109 copycolor(ctmp, fovscan(y)[x]);
110 scalecolor(ctmp, t2);
111 addcolor(vsum, ctmp);
112 t2sum += t2;
113 }
114 /* VADAPT of original is subtracted in addveil() */
115 if (t2sum > FTINY)
116 scalecolor(vsum, VADAPT/t2sum);
117 copycolor(veilscan(py)[px], vsum);
118 }
119 /* modify FOV sample image */
120 for (y = 0; y < fvyr; y++)
121 for (x = 0; x < fvxr; x++) {
122 scalecolor(fovscan(y)[x], 1.-VADAPT);
123 addcolor(fovscan(y)[x], veilscan(y)[x]);
124 }
125 comphist(); /* recompute histogram */
126 }
127
128
129 #if ADJ_VEIL
130 /*
131 * The following veil adjustment was added to compensate for
132 * the fact that contrast reduction gets confused with veil
133 * in the human visual system. Therefore, we reduce the
134 * veil in portions of the image where our mapping has
135 * already reduced contrast below the target value.
136 * This gets called after the intial veil has been computed
137 * and added to the foveal image, and the mapping has been
138 * determined.
139 */
140 extern void
141 adjveil(void) /* adjust veil image */
142 {
143 float *crfptr = crfimg;
144 COLOR *fovptr = fovimg;
145 COLOR *veilptr = veilimg;
146 double s2nits = 1./inpexp;
147 double vl, vl2, fovl, vlsum;
148 double deltavc[3];
149 int i, j;
150
151 if (lumf == rgblum)
152 s2nits *= WHTEFFICACY;
153
154 for (i = fvxr*fvyr; i--; crfptr++, fovptr++, veilptr++) {
155 if (crfptr[0] >= 0.95)
156 continue;
157 vl = plum(veilptr[0]);
158 fovl = (plum(fovptr[0]) - vl) * (1./(1.-VADAPT));
159 if (vl <= 0.05*fovl)
160 continue;
161 vlsum = vl;
162 for (j = 2; j < 11; j++) {
163 vlsum += crfptr[0]*vl - (1.0 - crfptr[0])*fovl;
164 vl2 = vlsum / (double)j;
165 if (vl2 < 0.0)
166 vl2 = 0.0;
167 crfptr[0] = crfactor(fovl + vl2);
168 }
169 /* desaturation code causes color fringes at this level */
170 for (j = 3; j--; ) {
171 double vc = colval(veilptr[0],j);
172 double fovc = (colval(fovptr[0],j) - vc) *
173 (1./(1.-VADAPT));
174 deltavc[j] = (1.-crfptr[0])*(fovl/s2nits - fovc);
175 if (vc + deltavc[j] < 0.0)
176 break;
177 }
178 if (j < 0)
179 addcolor(veilptr[0], deltavc);
180 else
181 scalecolor(veilptr[0], vl2/vl);
182 }
183 smoothveil(); /* smooth our result */
184 }
185
186
187 static void
188 smoothveil(void) /* smooth veil image */
189 {
190 COLOR *nveilimg;
191 COLOR *ovptr, *nvptr;
192 int x, y, i;
193
194 nveilimg = (COLOR *)malloc(fvxr*fvyr*sizeof(COLOR));
195 if (nveilimg == NULL)
196 return;
197 for (y = 1; y < fvyr-1; y++) {
198 ovptr = veilimg + y*fvxr + 1;
199 nvptr = nveilimg + y*fvxr + 1;
200 for (x = 1; x < fvxr-1; x++, ovptr++, nvptr++)
201 for (i = 3; i--; )
202 nvptr[0][i] = 0.5 * ovptr[0][i]
203 + (1./12.) *
204 (ovptr[-1][i] + ovptr[-fvxr][i] +
205 ovptr[1][i] + ovptr[fvxr][i])
206 + (1./24.) *
207 (ovptr[-fvxr-1][i] + ovptr[-fvxr+1][i] +
208 ovptr[fvxr-1][i] + ovptr[fvxr+1][i]);
209 }
210 ovptr = veilimg + 1;
211 nvptr = nveilimg + 1;
212 for (x = 1; x < fvxr-1; x++, ovptr++, nvptr++)
213 for (i = 3; i--; )
214 nvptr[0][i] = 0.5 * ovptr[0][i]
215 + (1./9.) *
216 (ovptr[-1][i] + ovptr[1][i] + ovptr[fvxr][i])
217 + (1./12.) *
218 (ovptr[fvxr-1][i] + ovptr[fvxr+1][i]);
219 ovptr = veilimg + (fvyr-1)*fvxr + 1;
220 nvptr = nveilimg + (fvyr-1)*fvxr + 1;
221 for (x = 1; x < fvxr-1; x++, ovptr++, nvptr++)
222 for (i = 3; i--; )
223 nvptr[0][i] = 0.5 * ovptr[0][i]
224 + (1./9.) *
225 (ovptr[-1][i] + ovptr[1][i] + ovptr[-fvxr][i])
226 + (1./12.) *
227 (ovptr[-fvxr-1][i] + ovptr[-fvxr+1][i]);
228 ovptr = veilimg + fvxr;
229 nvptr = nveilimg + fvxr;
230 for (y = 1; y < fvyr-1; y++, ovptr += fvxr, nvptr += fvxr)
231 for (i = 3; i--; )
232 nvptr[0][i] = 0.5 * ovptr[0][i]
233 + (1./9.) *
234 (ovptr[-fvxr][i] + ovptr[1][i] + ovptr[fvxr][i])
235 + (1./12.) *
236 (ovptr[-fvxr+1][i] + ovptr[fvxr+1][i]);
237 ovptr = veilimg + fvxr - 1;
238 nvptr = nveilimg + fvxr - 1;
239 for (y = 1; y < fvyr-1; y++, ovptr += fvxr, nvptr += fvxr)
240 for (i = 3; i--; )
241 nvptr[0][i] = 0.5 * ovptr[0][i]
242 + (1./9.) *
243 (ovptr[-fvxr][i] + ovptr[-1][i] + ovptr[fvxr][i])
244 + (1./12.) *
245 (ovptr[-fvxr-1][i] + ovptr[fvxr-1][i]);
246 for (i = 3; i--; ) {
247 nveilimg[0][i] = veilimg[0][i];
248 nveilimg[fvxr-1][i] = veilimg[fvxr-1][i];
249 nveilimg[(fvyr-1)*fvxr][i] = veilimg[(fvyr-1)*fvxr][i];
250 nveilimg[fvyr*fvxr-1][i] = veilimg[fvyr*fvxr-1][i];
251 }
252 free((void *)veilimg);
253 veilimg = nveilimg;
254 }
255 #endif
256
257 extern void
258 addveil( /* add veil to scanline */
259 COLOR *sl,
260 int y
261 )
262 {
263 int vx, vy;
264 double dx, dy;
265 double lv, uv;
266 register int x, i;
267
268 vy = dy = (y+.5)/numscans(&inpres)*fvyr - .5;
269 while (vy >= fvyr-1) vy--;
270 dy -= (double)vy;
271 for (x = 0; x < scanlen(&inpres); x++) {
272 vx = dx = (x+.5)/scanlen(&inpres)*fvxr - .5;
273 while (vx >= fvxr-1) vx--;
274 dx -= (double)vx;
275 for (i = 0; i < 3; i++) {
276 lv = (1.-dy)*colval(veilscan(vy)[vx],i) +
277 dy*colval(veilscan(vy+1)[vx],i);
278 uv = (1.-dy)*colval(veilscan(vy)[vx+1],i) +
279 dy*colval(veilscan(vy+1)[vx+1],i);
280 colval(sl[x],i) = (1.-VADAPT)*colval(sl[x],i) +
281 (1.-dx)*lv + dx*uv;
282 }
283 }
284 }
285
286
287 /****************** ACUITY STUFF *******************/
288
289 typedef struct {
290 short sampe; /* sample area size (exponent of 2) */
291 short nscans; /* number of scanlines in this bar */
292 int len; /* individual scanline length */
293 int nread; /* number of scanlines loaded */
294 COLOR *sdata; /* scanbar data */
295 } SCANBAR;
296
297 #define bscan(sb,y) ((COLOR *)(sb)->sdata+((y)%(sb)->nscans)*(sb)->len)
298
299 SCANBAR *rootbar; /* root scan bar (lowest resolution) */
300
301 float *inpacuD = NULL; /* input acuity data (cycles/degree) */
302
303 #define tsampr(x,y) inpacuD[(y)*fvxr+(x)]
304
305 static COLOR * getascan(SCANBAR *sb, int y);
306 static void acusample(COLOR col, int x, int y, double sr);
307 static void ascanval(COLOR col, int x, int y, SCANBAR *sb);
308 static SCANBAR *sballoc(int se, int ns, int sl);
309
310 extern double
311 hacuity( /* return visual acuity in cycles/degree */
312 double La
313 )
314 {
315 /* functional fit */
316 return(17.25*atan(1.4*log10(La) + 0.35) + 25.72);
317 }
318
319
320 static COLOR *
321 getascan( /* find/read scanline y for scanbar sb */
322 register SCANBAR *sb,
323 int y
324 )
325 {
326 register COLOR *sl0, *sl1, *mysl;
327 register int i;
328
329 if (y < sb->nread - sb->nscans) /* too far back? */
330 return(NULL);
331 for ( ; y >= sb->nread; sb->nread++) { /* read as necessary */
332 mysl = bscan(sb, sb->nread);
333 if (sb->sampe == 0) {
334 if (freadscan(mysl, sb->len, infp) < 0) {
335 fprintf(stderr, "%s: %s: scanline read error\n",
336 progname, infn);
337 exit(1);
338 }
339 } else {
340 sl0 = getascan(sb+1, 2*y);
341 if (sl0 == NULL)
342 return(NULL);
343 sl1 = getascan(sb+1, 2*y+1);
344 for (i = 0; i < sb->len; i++) {
345 copycolor(mysl[i], sl0[2*i]);
346 addcolor(mysl[i], sl0[2*i+1]);
347 addcolor(mysl[i], sl1[2*i]);
348 addcolor(mysl[i], sl1[2*i+1]);
349 scalecolor(mysl[i], 0.25);
350 }
351 }
352 }
353 return(bscan(sb, y));
354 }
355
356
357 extern void
358 acuscan( /* get acuity-sampled scanline */
359 COLOR *scln,
360 int y
361 )
362 {
363 double sr;
364 double dx, dy;
365 int ix, iy;
366 register int x;
367
368 if (inpacuD == NULL)
369 return;
370 /* compute foveal y position */
371 iy = dy = (y+.5)/numscans(&inpres)*fvyr - .5;
372 while (iy >= fvyr-1) iy--;
373 dy -= (double)iy;
374 for (x = 0; x < scanlen(&inpres); x++) {
375 /* compute foveal x position */
376 ix = dx = (x+.5)/scanlen(&inpres)*fvxr - .5;
377 while (ix >= fvxr-1) ix--;
378 dx -= (double)ix;
379 /* interpolate sample rate */
380 sr = (1.-dy)*((1.-dx)*tsampr(ix,iy) + dx*tsampr(ix+1,iy)) +
381 dy*((1.-dx)*tsampr(ix,iy+1) + dx*tsampr(ix+1,iy+1));
382
383 acusample(scln[x], x, y, sr); /* compute sample */
384 }
385 }
386
387
388 static void
389 acusample( /* interpolate sample at (x,y) using rate sr */
390 COLOR col,
391 int x,
392 int y,
393 double sr
394 )
395 {
396 COLOR c1;
397 double d;
398 register SCANBAR *sb0;
399
400 for (sb0 = rootbar; sb0->sampe != 0 && 1<<sb0[1].sampe > sr; sb0++)
401 ;
402 ascanval(col, x, y, sb0);
403 if (sb0->sampe == 0) /* don't extrapolate highest */
404 return;
405 ascanval(c1, x, y, sb0+1);
406 d = ((1<<sb0->sampe) - sr)/(1<<sb0[1].sampe);
407 scalecolor(col, 1.-d);
408 scalecolor(c1, d);
409 addcolor(col, c1);
410 }
411
412
413 static void
414 ascanval( /* interpolate scanbar at orig. coords (x,y) */
415 COLOR col,
416 int x,
417 int y,
418 SCANBAR *sb
419 )
420 {
421 COLOR *sl0, *sl1, c1, c1y;
422 double dx, dy;
423 int ix, iy;
424
425 if (sb->sampe == 0) { /* no need to interpolate */
426 sl0 = getascan(sb, y);
427 copycolor(col, sl0[x]);
428 return;
429 }
430 /* compute coordinates for sb */
431 ix = dx = (x+.5)/(1<<sb->sampe) - .5;
432 while (ix >= sb->len-1) ix--;
433 dx -= (double)ix;
434 iy = dy = (y+.5)/(1<<sb->sampe) - .5;
435 while (iy >= (numscans(&inpres)>>sb->sampe)-1) iy--;
436 dy -= (double)iy;
437 /* get scanlines */
438 sl0 = getascan(sb, iy);
439 #ifdef DEBUG
440 if (sl0 == NULL)
441 error(INTERNAL, "cannot backspace in ascanval");
442 #endif
443 sl1 = getascan(sb, iy+1);
444 /* 2D linear interpolation */
445 copycolor(col, sl0[ix]);
446 scalecolor(col, 1.-dx);
447 copycolor(c1, sl0[ix+1]);
448 scalecolor(c1, dx);
449 addcolor(col, c1);
450 copycolor(c1y, sl1[ix]);
451 scalecolor(c1y, 1.-dx);
452 copycolor(c1, sl1[ix+1]);
453 scalecolor(c1, dx);
454 addcolor(c1y, c1);
455 scalecolor(col, 1.-dy);
456 scalecolor(c1y, dy);
457 addcolor(col, c1y);
458 for (ix = 0; ix < 3; ix++) /* make sure no negative */
459 if (colval(col,ix) < 0.)
460 colval(col,ix) = 0.;
461 }
462
463
464 static SCANBAR *
465 sballoc( /* allocate scanbar */
466 int se, /* sampling rate exponent */
467 int ns, /* number of scanlines */
468 int sl /* original scanline length */
469 )
470 {
471 SCANBAR *sbarr;
472 register SCANBAR *sb;
473
474 sbarr = sb = (SCANBAR *)malloc((se+1)*sizeof(SCANBAR));
475 if (sb == NULL)
476 syserror("malloc");
477 do {
478 sb->len = sl>>se;
479 if (sb->len <= 0)
480 continue;
481 sb->sampe = se;
482 sb->nscans = ns;
483 sb->sdata = (COLOR *)malloc(sb->len*ns*sizeof(COLOR));
484 if (sb->sdata == NULL)
485 syserror("malloc");
486 sb->nread = 0;
487 ns <<= 1;
488 sb++;
489 } while (--se >= 0);
490 return(sbarr);
491 }
492
493
494 extern void
495 initacuity(void) /* initialize variable acuity sampling */
496 {
497 FVECT diffx, diffy, cp;
498 double omega, maxsr;
499 register int x, y, i;
500
501 compraydir(); /* compute ray directions */
502
503 inpacuD = (float *)malloc(fvxr*fvyr*sizeof(float));
504 if (inpacuD == NULL)
505 syserror("malloc");
506 maxsr = 1.; /* compute internal sample rates */
507 for (y = 1; y < fvyr-1; y++)
508 for (x = 1; x < fvxr-1; x++) {
509 for (i = 0; i < 3; i++) {
510 diffx[i] = 0.5*fvxr/scanlen(&inpres) *
511 (rdirscan(y)[x+1][i] -
512 rdirscan(y)[x-1][i]);
513 diffy[i] = 0.5*fvyr/numscans(&inpres) *
514 (rdirscan(y+1)[x][i] -
515 rdirscan(y-1)[x][i]);
516 }
517 fcross(cp, diffx, diffy);
518 omega = 0.5 * sqrt(DOT(cp,cp));
519 if (omega <= FTINY*FTINY)
520 tsampr(x,y) = 1.;
521 else if ((tsampr(x,y) = PI/180. / sqrt(omega) /
522 hacuity(plum(fovscan(y)[x]))) > maxsr)
523 maxsr = tsampr(x,y);
524 }
525 /* copy perimeter (easier) */
526 for (x = 1; x < fvxr-1; x++) {
527 tsampr(x,0) = tsampr(x,1);
528 tsampr(x,fvyr-1) = tsampr(x,fvyr-2);
529 }
530 for (y = 0; y < fvyr; y++) {
531 tsampr(0,y) = tsampr(1,y);
532 tsampr(fvxr-1,y) = tsampr(fvxr-2,y);
533 }
534 /* initialize with next power of two */
535 rootbar = sballoc((int)(log(maxsr)/log(2.))+1, 2, scanlen(&inpres));
536 }