ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/pcond.c
Revision: 3.1
Committed: Thu Oct 3 16:52:46 1996 UTC (27 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Initial revision

File Contents

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