ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pcond.c
Revision: 3.8
Committed: Thu Jan 9 13:56:23 1997 UTC (27 years, 3 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.7: +33 -26 lines
Log Message:
changed to histogram truncation
made veiling reflection modify FOV sample image

File Contents

# User Rev Content
1 greg 3.7 /* Copyright (c) 1997 Regents of the University of California */
2 greg 3.1
3     #ifndef lint
4     static char SCCSid[] = "$SunId$ LBL";
5     #endif
6    
7     /*
8     * Condition Radiance picture for display/output
9     */
10    
11     #include "pcond.h"
12    
13    
14     #define LDMAX 100 /* default max. display luminance */
15     #define LDMINF 0.01 /* default min. display lum. factor */
16    
17     int what2do = 0; /* desired adjustments */
18    
19     double ldmax = LDMAX; /* maximum output luminance */
20     double ldmin = 0.; /* minimum output luminance */
21     double Bldmin, Bldmax; /* Bl(ldmin) and Bl(ldmax) */
22    
23     char *progname; /* global argv[0] */
24    
25     char *infn; /* input file name */
26     FILE *infp; /* input stream */
27 greg 3.6 FILE *mapfp = NULL; /* tone-mapping function stream */
28 greg 3.1 VIEW ourview = STDVIEW; /* picture view */
29     int gotview = 0; /* picture has view */
30     double pixaspect = 1.0; /* pixel aspect ratio */
31     RESOLU inpres; /* input picture resolution */
32    
33     COLOR *fovimg; /* foveal (1 degree) averaged image */
34     short fvxr, fvyr; /* foveal image resolution */
35 greg 3.8 float bwhist[HISTRES]; /* luminance histogram */
36     double histot; /* total count of histogram */
37 greg 3.1 double bwmin, bwmax; /* histogram limits */
38     double bwavg; /* mean brightness */
39    
40     double scalef = 0.; /* linear scaling factor */
41    
42    
43     main(argc, argv)
44     int argc;
45     char *argv[];
46     {
47     static RGBPRIMS outprimS;
48     int i;
49     #define bool(flg) switch (argv[i][2]) { \
50     case '\0': what2do ^= flg; break; \
51     case 'y': case 'Y': case 't': case 'T': \
52     case '+': case '1': what2do |= flg; break; \
53     case 'n': case 'N': case 'f': case 'F': \
54     case '-': case '0': what2do &= ~(flg); break; \
55     default: goto userr; }
56    
57     progname = argv[0];
58    
59     for (i = 1; i < argc && argv[i][0] == '-'; i++)
60     switch (argv[i][1]) {
61     case 'h':
62     bool(DO_HUMAN);
63     break;
64     case 'a':
65     bool(DO_ACUITY);
66     break;
67     case 'v':
68     bool(DO_VEIL);
69     break;
70     case 's':
71     bool(DO_HSENS);
72     break;
73     case 'c':
74     bool(DO_COLOR);
75     break;
76     case 'w':
77     bool(DO_CWEIGHT);
78     break;
79     case 'l':
80     bool(DO_LINEAR);
81     break;
82     case 'p':
83     if (i+8 >= argc) goto userr;
84     outprimS[RED][CIEX] = atof(argv[++i]);
85     outprimS[RED][CIEY] = atof(argv[++i]);
86     outprimS[GRN][CIEX] = atof(argv[++i]);
87     outprimS[GRN][CIEY] = atof(argv[++i]);
88     outprimS[BLU][CIEX] = atof(argv[++i]);
89     outprimS[BLU][CIEY] = atof(argv[++i]);
90     outprimS[WHT][CIEX] = atof(argv[++i]);
91     outprimS[WHT][CIEY] = atof(argv[++i]);
92     outprims = outprimS;
93     break;
94     case 'e':
95     if (i+1 >= argc) goto userr;
96     scalef = atof(argv[++i]);
97     if (argv[i][0] == '+' | argv[i][0] == '-')
98     scalef = pow(2.0, scalef);
99     what2do |= DO_LINEAR;
100     break;
101     case 'f':
102     if (i+1 >= argc) goto userr;
103     mbcalfile = argv[++i];
104     break;
105     case 't':
106     if (i+1 >= argc) goto userr;
107     ldmax = atof(argv[++i]);
108     if (ldmax <= FTINY)
109     goto userr;
110     break;
111     case 'b':
112     if (i+1 >= argc) goto userr;
113     ldmin = atof(argv[++i]);
114     break;
115 greg 3.6 case 'm':
116     if (i+1 >= argc) goto userr;
117     if ((mapfp = fopen(argv[++i], "w")) == NULL) {
118     fprintf(stderr,
119     "%s: cannot open for writing\n",
120     argv[i]);
121     exit(1);
122     }
123     break;
124 greg 3.1 default:
125     goto userr;
126     }
127     if (mbcalfile != NULL & outprims != stdprims) {
128     fprintf(stderr, "%s: only one of -p or -f option supported\n",
129     progname);
130     exit(1);
131     }
132     if (outprims == stdprims & inprims != stdprims)
133     outprims = inprims;
134     if (ldmin <= FTINY)
135     ldmin = ldmax*LDMINF;
136     else if (ldmin >= ldmax) {
137     fprintf(stderr, "%s: Ldmin (%f) >= Ldmax (%f)!\n", progname,
138     ldmin, ldmax);
139     exit(1);
140     }
141     Bldmin = Bl(ldmin);
142     Bldmax = Bl(ldmax);
143     if (i >= argc || i+2 < argc)
144     goto userr;
145     if ((infp = fopen(infn=argv[i], "r")) == NULL)
146     syserror(infn);
147     if (i+2 == argc && freopen(argv[i+1], "w", stdout) == NULL)
148     syserror(argv[i+1]);
149     #ifdef MSDOS
150     setmode(fileno(infp), O_BINARY);
151     setmode(fileno(stdout), O_BINARY);
152     #endif
153     getahead(); /* load input header */
154     printargs(argc, argv, stdout); /* add to output header */
155     if (outprims != inprims)
156     fputprims(outprims, stdout);
157     mapimage(); /* map the picture */
158 greg 3.6 if (mapfp != NULL) /* write out basic mapping */
159     putmapping(mapfp);
160 greg 3.1 exit(0);
161     userr:
162 greg 3.6 fprintf(stderr, "Usage: %s [-{h|a|v|s|c|l|w}[+-]][-e ev][-p xr yr xg yg xb yb xw yw|-f mbf.cal][-t Ldmax][-b Ldmin][-m mapfile] inpic [outpic]\n",
163 greg 3.1 progname);
164     exit(1);
165     #undef bool
166     }
167    
168    
169     syserror(s) /* report system error and exit */
170     char *s;
171     {
172     fprintf(stderr, "%s: ", progname);
173     perror(s);
174     exit(2);
175     }
176    
177    
178     headline(s) /* process header line */
179     char *s;
180     {
181     static RGBPRIMS inprimS;
182     char fmt[32];
183    
184     if (formatval(fmt, s)) { /* check if format string */
185     if (!strcmp(fmt,COLRFMT)) lumf = rgblum;
186     else if (!strcmp(fmt,CIEFMT)) lumf = cielum;
187     else lumf = NULL;
188     return; /* don't echo */
189     }
190     if (isprims(s)) { /* get input primaries */
191     primsval(inprimS, s);
192     inprims= inprimS;
193     return; /* don't echo */
194     }
195     if (isexpos(s)) { /* picture exposure */
196     inpexp *= exposval(s);
197     return; /* don't echo */
198     }
199     if (isaspect(s)) /* pixel aspect ratio */
200     pixaspect *= aspectval(s);
201     if (isview(s)) /* image view */
202     gotview += sscanview(&ourview, s);
203     fputs(s, stdout);
204     }
205    
206    
207     getahead() /* load picture header */
208     {
209     char *err;
210    
211     getheader(infp, headline, NULL);
212     if (lumf == NULL || !fgetsresolu(&inpres, infp)) {
213     fprintf(stderr, "%s: %s: not a Radiance picture\n",
214     progname, infn);
215     exit(1);
216     }
217     if (lumf == rgblum)
218     comprgb2xyzmat(inrgb2xyz, inprims);
219     else if (mbcalfile != NULL) {
220     fprintf(stderr, "%s: macbethcal only works with RGB pictures\n",
221     progname);
222     exit(1);
223     }
224     if (!gotview || ourview.type == VT_PAR) {
225     copystruct(&ourview, &stdview);
226     ourview.type = VT_PER;
227     if (pixaspect*inpres.yr < inpres.xr) {
228     ourview.horiz = 40.0;
229     ourview.vert = 2.*180./PI *
230     atan(.3639702*pixaspect*inpres.yr/inpres.xr);
231     } else {
232     ourview.vert = 40.0;
233     ourview.horiz = 2.*180./PI *
234     atan(.3639702*inpres.xr/pixaspect/inpres.yr);
235     }
236     }
237     if ((err = setview(&ourview)) != NULL) {
238     fprintf(stderr, "%s: view error in picture \"%s\": %s\n",
239     progname, infn, err);
240     exit(1);
241     }
242     }
243    
244    
245     mapimage() /* map picture and send to stdout */
246     {
247     COLOR *scan;
248    
249     #ifdef DEBUG
250     fprintf(stderr, "%s: generating histogram...", progname);
251     #endif
252 greg 3.8 getfovimg(); /* get foveal sample image */
253     comphist(); /* generate adaptation histogram */
254 greg 3.1 #ifdef DEBUG
255     fputs("done\n", stderr);
256     #endif
257     check2do(); /* modify what2do flags */
258     if (what2do&DO_VEIL) {
259     #ifdef DEBUG
260     fprintf(stderr, "%s: computing veiling...", progname);
261     #endif
262     compveil();
263     #ifdef DEBUG
264     fputs("done\n", stderr);
265     #endif
266     }
267     #ifdef DEBUG
268     fprintf(stderr, "%s: computing brightness mapping...", progname);
269     #endif
270     if (!(what2do&DO_LINEAR) && mkbrmap() < 0) { /* make tone map */
271     what2do |= DO_LINEAR; /* use linear scaling */
272     #ifdef DEBUG
273     fputs("failed!\n", stderr);
274     } else
275     fputs("done\n", stderr);
276     #else
277     }
278     #endif
279     if (what2do&DO_LINEAR) {
280     if (scalef <= FTINY) {
281     if (what2do&DO_HSENS)
282     scalef = htcontrs(Lb(0.5*(Bldmax+Bldmin))) /
283     htcontrs(Lb(bwavg));
284     else
285     scalef = Lb(0.5*(Bldmax+Bldmin)) / Lb(bwavg);
286     scalef *= WHTEFFICACY/(inpexp*ldmax);
287     }
288     #ifdef DEBUG
289     fprintf(stderr, "%s: linear scaling factor = %f\n",
290     progname, scalef);
291     #endif
292 greg 3.4 fputexpos(inpexp*scalef, stdout); /* record exposure */
293 greg 3.1 if (lumf == cielum) scalef /= WHTEFFICACY;
294     }
295     putchar('\n'); /* complete header */
296     fputsresolu(&inpres, stdout); /* resolution doesn't change */
297    
298     for (scan = firstscan(); scan != NULL; scan = nextscan())
299     if (fwritescan(scan, scanlen(&inpres), stdout) < 0) {
300     fprintf(stderr, "%s: scanline write error\n",
301     progname);
302     exit(1);
303     }
304     }
305    
306    
307     double
308     centprob(x, y) /* center-weighting probability function */
309     int x, y;
310     {
311 greg 3.8 double xr, yr, p;
312     /* paraboloid, 0 at 90 degrees from center */
313     xr = (x - .5*(fvxr-1))/90.; /* 180 degree fisheye has fv?r == 90 */
314     yr = (y - .5*(fvyr-1))/90.;
315     p = 1. - xr*xr - yr*yr;
316     return(p < 0. ? 0. : p);
317 greg 3.1 }
318    
319    
320 greg 3.8 getfovimg() /* load foveal sampled image */
321 greg 3.1 {
322     extern FILE *popen();
323     char combuf[128];
324     FILE *fp;
325     int x, y;
326 greg 3.8 /* compute image size */
327 greg 3.1 fvxr = sqrt(ourview.hn2)/FOVDIA + 0.5;
328     if (fvxr < 2) fvxr = 2;
329     fvyr = sqrt(ourview.vn2)/FOVDIA + 0.5;
330     if (fvyr < 2) fvyr = 2;
331     if (!(inpres.or & YMAJOR)) { /* picture is rotated? */
332     y = fvyr;
333     fvyr = fvxr;
334     fvxr = y;
335     }
336     if ((fovimg = (COLOR *)malloc(fvxr*fvyr*sizeof(COLOR))) == NULL)
337     syserror("malloc");
338     sprintf(combuf, "pfilt -1 -b -x %d -y %d %s", fvxr, fvyr, infn);
339     if ((fp = popen(combuf, "r")) == NULL)
340     syserror("popen");
341     getheader(fp, NULL, NULL); /* skip header */
342     if (fgetresolu(&x, &y, fp) < 0 || x != fvxr | y != fvyr)
343     goto readerr;
344     for (y = 0; y < fvyr; y++)
345     if (freadscan(fovscan(y), fvxr, fp) < 0)
346     goto readerr;
347     pclose(fp);
348 greg 3.8 return;
349     readerr:
350     fprintf(stderr, "%s: error reading from pfilt process in fovimage\n",
351     progname);
352     exit(1);
353     }
354    
355    
356     comphist() /* create foveal sampled image and histogram */
357     {
358     double l, b, lwmin, lwmax;
359     register int x, y;
360    
361 greg 3.1 lwmin = 1e10; /* find extrema */
362     lwmax = 0.;
363     for (y = 0; y < fvyr; y++)
364     for (x = 0; x < fvxr; x++) {
365     l = plum(fovscan(y)[x]);
366     if (l < lwmin) lwmin = l;
367     if (l > lwmax) lwmax = l;
368     }
369 greg 3.8 lwmin -= FTINY;
370     lwmax += FTINY;
371 greg 3.1 if (lwmin < LMIN) lwmin = LMIN;
372     if (lwmax > LMAX) lwmax = LMAX;
373     /* compute histogram */
374 greg 3.7 bwmin = Bl(lwmin);
375     bwmax = Bl(lwmax);
376 greg 3.1 bwavg = 0.;
377     for (y = 0; y < fvyr; y++)
378     for (x = 0; x < fvxr; x++) {
379     l = plum(fovscan(y)[x]);
380 greg 3.7 if (l < lwmin) continue;
381     if (l > lwmax) continue;
382 greg 3.1 b = Bl(l);
383     bwavg += b;
384 greg 3.8 l = what2do&DO_CWEIGHT ? centprob(x,y) : 1.;
385     bwhist[bwhi(b)] += l;
386     histot += l;
387 greg 3.1 }
388 greg 3.8 bwavg /= histot;
389 greg 3.1 }
390    
391    
392     check2do() /* check histogram to see what isn't worth doing */
393     {
394 greg 3.8 double sum;
395 greg 3.1 double b, l;
396     register int i;
397    
398     /* check for within display range */
399 greg 3.3 l = Lb(bwmax)/Lb(bwmin);
400     if (l <= ldmax/ldmin)
401 greg 3.1 what2do |= DO_LINEAR;
402 greg 3.3 /* determine if veiling significant */
403     if (l < 100.) /* heuristic */
404     what2do &= ~DO_VEIL;
405 greg 3.1
406 greg 3.3 if (!(what2do & (DO_ACUITY|DO_COLOR)))
407 greg 3.1 return;
408     /* find 5th percentile */
409 greg 3.8 sum = histot*0.05;
410 greg 3.1 for (i = 0; i < HISTRES; i++)
411     if ((sum -= bwhist[i]) <= 0)
412     break;
413     b = (i+.5)*(bwmax-bwmin)/HISTRES + bwmin;
414     l = Lb(b);
415     /* determine if acuity adj. useful */
416     if (what2do&DO_ACUITY &&
417     hacuity(l) >= (inpres.xr/sqrt(ourview.hn2) +
418 greg 3.5 inpres.yr/sqrt(ourview.vn2))/(2.*180./PI))
419 greg 3.1 what2do &= ~DO_ACUITY;
420     /* color sensitivity loss? */
421 greg 3.5 if (l >= TopMesopic)
422 greg 3.1 what2do &= ~DO_COLOR;
423     }