ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/tonemap.c
Revision: 3.2
Committed: Tue Apr 15 20:04:43 1997 UTC (27 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.1: +41 -35 lines
Log Message:
changed TM_F_NOERRS to TM_F_NOSTDERR
skipped histogram clipping for TM_F_LINEAR in tmComputeMapping()

File Contents

# Content
1 /* Copyright (c) 1997 Regents of the University of California */
2
3 #ifndef lint
4 static char SCCSid[] = "$SunId$ LBL";
5 #endif
6
7 /*
8 * Tone mapping functions.
9 * See tonemap.h for detailed function descriptions.
10 */
11
12 #include <stdio.h>
13 #include <math.h>
14 #include "tmprivat.h"
15 #include "tmerrmsg.h"
16
17 #define exp10(x) exp(M_LN10*(x))
18
19 struct tmStruct *tmTop = NULL; /* current tone mapping stack */
20
21
22 int
23 tmErrorReturn(func, err) /* error return (with message) */
24 char *func;
25 int err;
26 {
27 if (tmTop != NULL && tmTop->flags & TM_F_NOSTDERR)
28 return(err);
29 fputs(func, stderr);
30 fputs(": ", stderr);
31 fputs(tmErrorMessage[err], stderr);
32 fputs("!\n", stderr);
33 return(err);
34 }
35
36
37 struct tmStruct *
38 tmInit(flags, monpri, gamval) /* initialize new tone mapping */
39 int flags;
40 RGBPRIMP monpri;
41 double gamval;
42 {
43 static char funcName[] = "tmInit";
44 COLORMAT cmat;
45 register struct tmStruct *tmnew;
46 register int i;
47 /* allocate structure */
48 tmnew = (struct tmStruct *)malloc(sizeof(struct tmStruct));
49 if (tmnew == NULL)
50 return(NULL);
51
52 tmnew->flags = flags & ~TM_F_UNIMPL;
53 /* set monitor transform */
54 if (monpri == NULL || monpri == stdprims || tmnew->flags & TM_F_BW) {
55 tmnew->monpri = stdprims;
56 tmnew->clf[RED] = rgb2xyzmat[1][0];
57 tmnew->clf[GRN] = rgb2xyzmat[1][1];
58 tmnew->clf[BLU] = rgb2xyzmat[1][2];
59 } else {
60 comprgb2xyzmat(cmat, tmnew->monpri=monpri);
61 tmnew->clf[RED] = cmat[1][0];
62 tmnew->clf[GRN] = cmat[1][1];
63 tmnew->clf[BLU] = cmat[1][2];
64 }
65 tmnew->clfb[RED] = 256.*tmnew->clf[RED] + .5;
66 tmnew->clfb[GRN] = 256.*tmnew->clf[GRN] + .5;
67 tmnew->clfb[BLU] = 256.*tmnew->clf[BLU] + .5;
68 tmnew->clfb[EXP] = COLXS;
69 /* set gamma value */
70 if (gamval < MINGAM)
71 tmnew->mongam = DEFGAM;
72 else
73 tmnew->mongam = gamval;
74 for (i = TM_GAMTSZ; i--; )
75 tmnew->gamb[i] = 256.*pow((i+.5)/TM_GAMTSZ, 1./tmnew->mongam);
76 /* set input transform */
77 tmnew->inppri = tmnew->monpri;
78 tmnew->cmat[0][0] = tmnew->cmat[1][1] = tmnew->cmat[2][2] =
79 tmnew->inpsf = WHTEFFICACY;
80 tmnew->inpsfb = TM_BRTSCALE*log(tmnew->inpsf) + .5;
81 tmnew->cmat[0][1] = tmnew->cmat[0][2] = tmnew->cmat[1][0] =
82 tmnew->cmat[1][2] = tmnew->cmat[2][0] = tmnew->cmat[2][1] = 0.;
83 tmnew->flags &= ~TM_F_NEEDMAT;
84 tmnew->brmin = tmnew->brmax = 0;
85 tmnew->histo = NULL;
86 tmnew->lumap = NULL;
87 tmnew->tmprev = NULL;
88
89 tmnew->flags |= TM_F_INITED;
90 /* make it current */
91 tmnew->tmprev = tmTop;
92 return(tmTop = tmnew);
93 }
94
95
96 int
97 tmSetSpace(pri, sf) /* set input color space for conversions */
98 RGBPRIMP pri;
99 double sf;
100 {
101 static char funcName[] = "tmSetSpace";
102 register int i, j;
103 /* error check */
104 if (tmTop == NULL)
105 returnErr(TM_E_TMINVAL);
106 if (sf <= 1e-12)
107 returnErr(TM_E_ILLEGAL);
108 /* check if no change */
109 if (pri == tmTop->inppri && FEQ(sf, tmTop->inpsf))
110 returnOK;
111 tmTop->inppri = pri; /* let's set it */
112 tmTop->inpsf = sf;
113 tmTop->inpsfb = TM_BRTSCALE*log(sf) + (sf>=1. ? .5 : -.5);
114
115 if (tmTop->flags & TM_F_BW) { /* color doesn't matter */
116 tmTop->monpri = tmTop->inppri; /* eliminate xform */
117 if (tmTop->inppri == TM_XYZPRIM) {
118 tmTop->clf[CIEX] = tmTop->clf[CIEZ] = 0.;
119 tmTop->clf[CIEY] = 1.;
120 tmTop->clfb[CIEX] = tmTop->clfb[CIEZ] = 0;
121 tmTop->clfb[CIEY] = 255;
122 } else {
123 comprgb2xyzmat(tmTop->cmat, tmTop->monpri);
124 tmTop->clf[RED] = tmTop->cmat[1][0];
125 tmTop->clf[GRN] = tmTop->cmat[1][1];
126 tmTop->clf[BLU] = tmTop->cmat[1][2];
127 }
128 tmTop->cmat[0][0] = tmTop->cmat[1][1] = tmTop->cmat[2][2] =
129 tmTop->inpsf;
130 tmTop->cmat[0][1] = tmTop->cmat[0][2] = tmTop->cmat[1][0] =
131 tmTop->cmat[1][2] = tmTop->cmat[2][0] = tmTop->cmat[2][1] = 0.;
132
133 } else if (tmTop->inppri == TM_XYZPRIM) /* input is XYZ */
134 compxyz2rgbmat(tmTop->cmat, tmTop->monpri);
135
136 else { /* input is RGB */
137 if (tmTop->inppri != tmTop->monpri &&
138 PRIMEQ(tmTop->inppri, tmTop->monpri))
139 tmTop->inppri = tmTop->monpri; /* no xform */
140 comprgb2rgbmat(tmTop->cmat, tmTop->inppri, tmTop->monpri);
141 }
142 for (i = 0; i < 3; i++)
143 for (j = 0; j < 3; j++)
144 tmTop->cmat[i][j] *= tmTop->inpsf;
145 if (tmTop->inppri != tmTop->monpri)
146 tmTop->flags |= TM_F_NEEDMAT;
147 else
148 tmTop->flags &= ~TM_F_NEEDMAT;
149 returnOK;
150 }
151
152
153 void
154 tmClearHisto() /* clear current histogram */
155 {
156 if (tmTop == NULL || tmTop->histo == NULL)
157 return;
158 free((char *)tmTop->histo);
159 tmTop->histo = NULL;
160 }
161
162
163 int
164 tmCvColors(ls, cs, scan, len) /* convert float colors */
165 TMbright *ls;
166 BYTE *cs;
167 COLOR *scan;
168 int len;
169 {
170 static char funcName[] = "tmCvColors";
171 COLOR cmon;
172 double lum, slum;
173 register double d;
174 register int i;
175
176 if (tmTop == NULL)
177 returnErr(TM_E_TMINVAL);
178 if (ls == NULL | scan == NULL | len <= 0)
179 returnErr(TM_E_ILLEGAL);
180 for (i = len; i--; ) {
181 if (tmTop->flags & TM_F_NEEDMAT) /* get monitor RGB */
182 colortrans(cmon, tmTop->cmat, scan[i]);
183 else {
184 cmon[RED] = tmTop->inpsf*scan[i][RED];
185 cmon[GRN] = tmTop->inpsf*scan[i][GRN];
186 cmon[BLU] = tmTop->inpsf*scan[i][BLU];
187 }
188 /* world luminance */
189 lum = tmTop->clf[RED]*cmon[RED] +
190 tmTop->clf[GRN]*cmon[GRN] +
191 tmTop->clf[BLU]*cmon[BLU] ;
192 /* check range */
193 if (clipgamut(cmon, lum, CGAMUT_LOWER, cblack, cwhite))
194 lum = tmTop->clf[RED]*cmon[RED] +
195 tmTop->clf[GRN]*cmon[GRN] +
196 tmTop->clf[BLU]*cmon[BLU] ;
197 if (lum < MINLUM) {
198 ls[i] = MINBRT-1; /* bogus value */
199 lum = MINLUM;
200 } else {
201 d = TM_BRTSCALE*log(lum); /* encode it */
202 ls[i] = d>0. ? (int)(d+.5) : (int)(d-.5);
203 }
204 if (cs == TM_NOCHROM) /* no color? */
205 continue;
206 if (tmTop->flags & TM_F_MESOPIC && lum < LMESUPPER) {
207 slum = scotlum(cmon); /* mesopic adj. */
208 if (lum < LMESLOWER)
209 cmon[RED] = cmon[GRN] = cmon[BLU] = slum;
210 else {
211 d = (lum - LMESLOWER)/(LMESUPPER - LMESLOWER);
212 if (tmTop->flags & TM_F_BW)
213 cmon[RED] = cmon[GRN] =
214 cmon[BLU] = d*lum;
215 else
216 scalecolor(cmon, d);
217 d = (1.-d)*slum;
218 cmon[RED] += d;
219 cmon[GRN] += d;
220 cmon[BLU] += d;
221 }
222 } else if (tmTop->flags & TM_F_BW) {
223 cmon[RED] = cmon[GRN] = cmon[BLU] = lum;
224 }
225 d = tmTop->clf[RED]*cmon[RED]/lum;
226 /* cs[3*i ] = d>.999 ? 255 : 256.*pow(d, 1./tmTop->mongam); */
227 cs[3*i ] = d>.999 ? 255 : tmTop->gamb[(int)(d*TM_GAMTSZ)];
228 d = tmTop->clf[GRN]*cmon[GRN]/lum;
229 /* cs[3*i+1] = d>.999 ? 255 : 256.*pow(d, 1./tmTop->mongam); */
230 cs[3*i+1] = d>.999 ? 255 : tmTop->gamb[(int)(d*TM_GAMTSZ)];
231 d = tmTop->clf[BLU]*cmon[BLU]/lum;
232 /* cs[3*i+2] = d>.999 ? 255 : 256.*pow(d, 1./tmTop->mongam); */
233 cs[3*i+2] = d>.999 ? 255 : tmTop->gamb[(int)(d*TM_GAMTSZ)];
234 }
235 returnOK;
236 }
237
238
239 int
240 tmAddHisto(ls, len, wt) /* add values to histogram */
241 register TMbright *ls;
242 int len;
243 int wt;
244 {
245 static char funcName[] = "tmAddHisto";
246 int sum, oldorig, oldlen, horig, hlen;
247 register int i, j;
248
249 if (len <= 0)
250 returnErr(TM_E_ILLEGAL);
251 if (tmTop == NULL)
252 returnErr(TM_E_TMINVAL);
253 /* first, grow limits */
254 if (tmTop->histo == NULL) {
255 for (i = len; i-- && ls[i] < MINBRT; )
256 ;
257 if (i < 0)
258 returnOK;
259 tmTop->brmin = tmTop->brmax = ls[i];
260 oldlen = 0;
261 } else {
262 oldorig = (tmTop->brmin-MINBRT)/HISTEP;
263 oldlen = (tmTop->brmax-MINBRT)/HISTEP + 1 - oldorig;
264 }
265 for (i = len; i--; ) {
266 if ((j = ls[i]) < MINBRT)
267 continue;
268 if (j < tmTop->brmin)
269 tmTop->brmin = j;
270 else if (j > tmTop->brmax)
271 tmTop->brmax = j;
272 }
273 horig = (tmTop->brmin-MINBRT)/HISTEP;
274 hlen = (tmTop->brmax-MINBRT)/HISTEP + 1 - horig;
275 if (hlen > oldlen) { /* (re)allocate histogram */
276 register int *newhist = (int *)calloc(hlen, sizeof(int));
277 if (newhist == NULL)
278 returnErr(TM_E_NOMEM);
279 if (oldlen) { /* copy and free old */
280 for (i = oldlen, j = i+oldorig-horig; i; )
281 newhist[--j] = tmTop->histo[--i];
282 free((char *)tmTop->histo);
283 }
284 tmTop->histo = newhist;
285 if (tmTop->lumap != NULL) { /* invalid tone map */
286 free((char *)tmTop->lumap);
287 tmTop->lumap = NULL;
288 }
289 }
290 if (wt == 0)
291 returnOK;
292 for (i = len; i--; ) /* add in new counts */
293 if (ls[i] >= MINBRT)
294 tmTop->histo[ (ls[i]-MINBRT)/HISTEP - horig ] += wt;
295 returnOK;
296 }
297
298
299 static double
300 htcontrs(La) /* human threshold contrast sensitivity, dL(La) */
301 double La;
302 {
303 double l10La, l10dL;
304 /* formula taken from Ferwerda et al. [SG96] */
305 if (La < 1.148e-4)
306 return(1.38e-3);
307 l10La = log10(La);
308 if (l10La < -1.44) /* rod response regime */
309 l10dL = pow(.405*l10La + 1.6, 2.18) - 2.86;
310 else if (l10La < -.0184)
311 l10dL = l10La - .395;
312 else if (l10La < 1.9) /* cone response regime */
313 l10dL = pow(.249*l10La + .65, 2.7) - .72;
314 else
315 l10dL = l10La - 1.255;
316
317 return(exp10(l10dL));
318 }
319
320
321 int
322 tmComputeMapping(gamval, Lddyn, Ldmax)
323 double gamval;
324 double Lddyn;
325 double Ldmax;
326 {
327 static char funcName[] = "tmComputeMapping";
328 int *histo;
329 float *cumf;
330 int brt0, histlen, histot, threshold, ceiling, trimmings;
331 double logLddyn, Ldmin, Ldavg, Lwavg, Tr, Lw, Ld;
332 int4 sum;
333 register double d;
334 register int i, j;
335
336 if (tmTop == NULL || tmTop->histo == NULL)
337 returnErr(TM_E_TMINVAL);
338 /* check arguments */
339 if (Lddyn < MINLDDYN) Lddyn = DEFLDDYN;
340 if (Ldmax < MINLDMAX) Ldmax = DEFLDMAX;
341 if (gamval < MINGAM) gamval = tmTop->mongam;
342 /* compute handy values */
343 Ldmin = Ldmax/Lddyn;
344 logLddyn = log(Lddyn);
345 Ldavg = sqrt(Ldmax*Ldmin);
346 i = (tmTop->brmin-MINBRT)/HISTEP;
347 brt0 = MINBRT + HISTEP/2 + i*HISTEP;
348 histlen = (tmTop->brmax-MINBRT)/HISTEP + 1 - i;
349 /* histogram total and mean */
350 histot = 0; sum = 0;
351 j = brt0 + histlen*HISTEP;
352 for (i = histlen; i--; ) {
353 histot += tmTop->histo[i];
354 sum += (j -= HISTEP) * tmTop->histo[i];
355 }
356 threshold = histot*.025 + .5;
357 if (threshold < 4)
358 returnErr(TM_E_TMFAIL);
359 Lwavg = tmLuminance( (double)sum / histot );
360 if (!(tmTop->flags & TM_F_LINEAR)) { /* clamp histogram */
361 histo = (int *)malloc(histlen*sizeof(int));
362 cumf = (float *)malloc((histlen+1)*sizeof(float));
363 if (histo == NULL | cumf == NULL)
364 returnErr(TM_E_NOMEM);
365 for (i = histlen; i--; ) /* make malleable copy */
366 histo[i] = tmTop->histo[i];
367 do { /* iterate to solution */
368 sum = 0; /* cumulative probability */
369 for (i = 0; i < histlen; i++) {
370 cumf[i] = (double)sum/histot;
371 sum += histo[i];
372 }
373 cumf[i] = 1.;
374 Tr = histot * (double)(tmTop->brmax - tmTop->brmin) /
375 ((double)histlen*TM_BRTSCALE) / logLddyn;
376 ceiling = Tr + 1.;
377 trimmings = 0; /* clip to envelope */
378 for (i = histlen; i--; ) {
379 if (tmTop->flags & TM_F_HCONTR) {
380 Lw = tmLuminance(brt0 + i*HISTEP);
381 Ld = Ldmin * exp( logLddyn *
382 .5*(cumf[i]+cumf[i+1]) );
383 ceiling = Tr * (htcontrs(Ld) * Lw) /
384 (htcontrs(Lw) * Ld) + 1.;
385 }
386 if (histo[i] > ceiling) {
387 trimmings += histo[i] - ceiling;
388 histo[i] = ceiling;
389 }
390 }
391 } while ((histot -= trimmings) > threshold &&
392 trimmings > threshold);
393 }
394 /* allocate luminance map */
395 if (tmTop->lumap == NULL) {
396 tmTop->lumap = (unsigned short *)malloc(
397 (tmTop->brmax-tmTop->brmin+1)*sizeof(unsigned short) );
398 if (tmTop->lumap == NULL)
399 returnErr(TM_E_NOMEM);
400 }
401 if (tmTop->flags & TM_F_LINEAR || histot <= threshold) {
402 /* linear tone mapping */
403 if (tmTop->flags & TM_F_HCONTR)
404 d = htcontrs(Ldavg) / htcontrs(Lwavg);
405 else
406 d = Ldavg / Lwavg;
407 d = log(d/Ldmax);
408 for (i = tmTop->brmax-tmTop->brmin+1; i--; )
409 tmTop->lumap[i] = 256. * exp(
410 ( d + (tmTop->brmin+i)/(double)TM_BRTSCALE )
411 / gamval );
412 } else {
413 /* histogram adjustment */
414 for (i = tmTop->brmax-tmTop->brmin+1; i--; ) {
415 j = d = (double)i/(tmTop->brmax-tmTop->brmin)*histlen;
416 d -= (double)j;
417 Ld = Ldmin*exp(logLddyn*((1.-d)*cumf[j]+d*cumf[j+1]));
418 d = (Ld - Ldmin)/(Ldmax - Ldmin);
419 tmTop->lumap[i] = 256.*pow(d, 1./gamval);
420 }
421 }
422 if (!(tmTop->flags & TM_F_LINEAR)) {
423 free((char *)histo);
424 free((char *)cumf);
425 }
426 returnOK;
427 }
428
429
430 int
431 tmMapPixels(ps, ls, cs, len)
432 register BYTE *ps;
433 TMbright *ls;
434 register BYTE *cs;
435 int len;
436 {
437 static char funcName[] = "tmMapPixels";
438 int rdiv, gdiv, bdiv;
439 register int li, pv;
440
441 if (tmTop == NULL || tmTop->lumap == NULL)
442 returnErr(TM_E_TMINVAL);
443 if (ps == NULL | ls == NULL | len <= 0)
444 returnErr(TM_E_ILLEGAL);
445 rdiv = tmTop->gamb[((int4)TM_GAMTSZ*tmTop->clfb[RED])>>8];
446 gdiv = tmTop->gamb[((int4)TM_GAMTSZ*tmTop->clfb[GRN])>>8];
447 bdiv = tmTop->gamb[((int4)TM_GAMTSZ*tmTop->clfb[BLU])>>8];
448 while (len--) {
449 if ((li = *ls++) < tmTop->brmin)
450 li = tmTop->brmin;
451 else if (li > tmTop->brmax)
452 li = tmTop->brmax;
453 li = tmTop->lumap[li - tmTop->brmin];
454 if (cs == TM_NOCHROM)
455 *ps++ = li>255 ? 255 : li;
456 else {
457 pv = *cs++ * li / rdiv;
458 *ps++ = pv>255 ? 255 : pv;
459 pv = *cs++ * li / gdiv;
460 *ps++ = pv>255 ? 255 : pv;
461 pv = *cs++ * li / bdiv;
462 *ps++ = pv>255 ? 255 : pv;
463 }
464 }
465 returnOK;
466 }
467
468
469 struct tmStruct *
470 tmPop() /* pop top tone mapping off stack */
471 {
472 register struct tmStruct *tms;
473
474 if ((tms = tmTop) != NULL)
475 tmTop = tms->tmprev;
476 return(tms);
477 }
478
479
480 int
481 tmPull(tms) /* pull a tone mapping from stack */
482 register struct tmStruct *tms;
483 {
484 register struct tmStruct *tms2;
485 /* special cases first */
486 if (tms == NULL | tmTop == NULL)
487 return(0);
488 if (tms == tmTop) {
489 tmTop = tms->tmprev;
490 tms->tmprev = NULL;
491 return(1);
492 }
493 for (tms2 = tmTop; tms2->tmprev != NULL; tms2 = tms2->tmprev)
494 if (tms == tms2->tmprev) { /* remove it */
495 tms2->tmprev = tms->tmprev;
496 tms->tmprev = NULL;
497 return(1);
498 }
499 return(0); /* not found on stack */
500 }
501
502
503 int
504 tmPush(tms) /* push tone mapping on top of stack */
505 register struct tmStruct *tms;
506 {
507 static char funcName[] = "tmPush";
508 /* check validity */
509 if (tms == NULL || !(tms->flags & TM_F_INITED))
510 returnErr(TM_E_ILLEGAL);
511 if (tms == tmTop) /* check necessity */
512 returnOK;
513 /* pull if already in stack */
514 (void)tmPull(tms);
515 /* push it on top */
516 tms->tmprev = tmTop;
517 tmTop = tms;
518 returnOK;
519 }
520
521
522 void
523 tmDone(tms) /* done with tone mapping -- destroy it */
524 register struct tmStruct *tms;
525 {
526 /* NULL arg. is equiv. to tmTop */
527 if (tms == NULL && (tms = tmTop) == NULL)
528 return;
529 /* take out of stack if present */
530 (void)tmPull(tms);
531 /* free tables */
532 if (tms->histo != NULL)
533 free((char *)tms->histo);
534 if (tms->lumap != NULL)
535 free((char *)tms->lumap);
536 tms->flags = 0;
537 free((char *)tms); /* free basic structure */
538 }