ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/tonemap.c
Revision: 3.1
Committed: Tue Apr 15 16:53:01 1997 UTC (27 years ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Initial revision

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_NOERRS)
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 /* allocate temporary tables */
350 histo = (int *)malloc(histlen*sizeof(int));
351 cumf = (float *)malloc((histlen+1)*sizeof(float));
352 if (histo == NULL | cumf == NULL)
353 returnErr(TM_E_NOMEM);
354 /* histogram total and mean */
355 histot = 0; sum = 0;
356 j = brt0 + histlen*HISTEP;
357 for (i = histlen; i--; ) {
358 histot += (histo[i] = tmTop->histo[i]);
359 sum += (j -= HISTEP) * histo[i];
360 }
361 threshold = histot*.025 + .5;
362 if (threshold < 4)
363 returnErr(TM_E_TMFAIL);
364 Lwavg = tmLuminance( (double)sum / histot );
365 do { /* iterate to solution */
366 sum = 0; /* compute cumulative probability */
367 for (i = 0; i < histlen; i++) {
368 cumf[i] = (double)sum/histot;
369 sum += histo[i];
370 }
371 cumf[i] = 1.;
372 Tr = histot * (double)(tmTop->brmax - tmTop->brmin) /
373 ((double)histlen*TM_BRTSCALE) / logLddyn;
374 ceiling = Tr + 1.;
375 trimmings = 0; /* clip to envelope */
376 for (i = histlen; i--; ) {
377 if (tmTop->flags & TM_F_HCONTR) {
378 Lw = tmLuminance(brt0 + i*HISTEP);
379 Ld = Ldmin * exp( logLddyn *
380 .5*(cumf[i]+cumf[i+1]) );
381 ceiling = Tr * (htcontrs(Ld) * Lw) /
382 (htcontrs(Lw) * Ld) + 1.;
383 }
384 if (histo[i] > ceiling) {
385 trimmings += histo[i] - ceiling;
386 histo[i] = ceiling;
387 }
388 }
389 } while ((histot -= trimmings) > threshold && trimmings > threshold);
390
391 if (tmTop->lumap == NULL) { /* allocate luminance map */
392 tmTop->lumap = (unsigned short *)malloc(
393 (tmTop->brmax-tmTop->brmin+1)*sizeof(unsigned short) );
394 if (tmTop->lumap == NULL)
395 returnErr(TM_E_NOMEM);
396 }
397 if (tmTop->flags & TM_F_LINEAR || histot <= threshold) {
398 /* linear tone mapping */
399 if (tmTop->flags & TM_F_HCONTR)
400 d = htcontrs(Ldavg) / htcontrs(Lwavg);
401 else
402 d = Ldavg / Lwavg;
403 d = log(d/Ldmax);
404 for (i = tmTop->brmax-tmTop->brmin+1; i--; )
405 tmTop->lumap[i] = 256. * exp(
406 ( d + (tmTop->brmin+i)/(double)TM_BRTSCALE )
407 / gamval );
408 } else {
409 /* histogram adjustment */
410 for (i = tmTop->brmax-tmTop->brmin+1; i--; ) {
411 j = d = (double)i/(tmTop->brmax-tmTop->brmin)*histlen;
412 d -= (double)j;
413 Ld = Ldmin*exp(logLddyn*((1.-d)*cumf[j]+d*cumf[j+1]));
414 d = (Ld - Ldmin)/(Ldmax - Ldmin);
415 tmTop->lumap[i] = 256.*pow(d, 1./gamval);
416 }
417 }
418 free((char *)histo);
419 free((char *)cumf);
420 returnOK;
421 }
422
423
424 int
425 tmMapPixels(ps, ls, cs, len)
426 register BYTE *ps;
427 TMbright *ls;
428 register BYTE *cs;
429 int len;
430 {
431 static char funcName[] = "tmMapPixels";
432 int rdiv, gdiv, bdiv;
433 register int li, pv;
434
435 if (tmTop == NULL || tmTop->lumap == NULL)
436 returnErr(TM_E_TMINVAL);
437 if (ps == NULL | ls == NULL | len <= 0)
438 returnErr(TM_E_ILLEGAL);
439 rdiv = tmTop->gamb[((int4)TM_GAMTSZ*tmTop->clfb[RED])>>8];
440 gdiv = tmTop->gamb[((int4)TM_GAMTSZ*tmTop->clfb[GRN])>>8];
441 bdiv = tmTop->gamb[((int4)TM_GAMTSZ*tmTop->clfb[BLU])>>8];
442 while (len--) {
443 if ((li = *ls++) < tmTop->brmin)
444 li = tmTop->brmin;
445 else if (li > tmTop->brmax)
446 li = tmTop->brmax;
447 li = tmTop->lumap[li - tmTop->brmin];
448 if (cs == TM_NOCHROM)
449 *ps++ = li>255 ? 255 : li;
450 else {
451 pv = *cs++ * li / rdiv;
452 *ps++ = pv>255 ? 255 : pv;
453 pv = *cs++ * li / gdiv;
454 *ps++ = pv>255 ? 255 : pv;
455 pv = *cs++ * li / bdiv;
456 *ps++ = pv>255 ? 255 : pv;
457 }
458 }
459 returnOK;
460 }
461
462
463 struct tmStruct *
464 tmPop() /* pop top tone mapping off stack */
465 {
466 register struct tmStruct *tms;
467
468 if ((tms = tmTop) != NULL)
469 tmTop = tms->tmprev;
470 return(tms);
471 }
472
473
474 int
475 tmPull(tms) /* pull a tone mapping from stack */
476 register struct tmStruct *tms;
477 {
478 register struct tmStruct *tms2;
479 /* special cases first */
480 if (tms == NULL | tmTop == NULL)
481 return(0);
482 if (tms == tmTop) {
483 tmTop = tms->tmprev;
484 tms->tmprev = NULL;
485 return(1);
486 }
487 for (tms2 = tmTop; tms2->tmprev != NULL; tms2 = tms2->tmprev)
488 if (tms == tms2->tmprev) { /* remove it */
489 tms2->tmprev = tms->tmprev;
490 tms->tmprev = NULL;
491 return(1);
492 }
493 return(0); /* not found on stack */
494 }
495
496
497 int
498 tmPush(tms) /* push tone mapping on top of stack */
499 register struct tmStruct *tms;
500 {
501 static char funcName[] = "tmPush";
502 /* check validity */
503 if (tms == NULL || !(tms->flags & TM_F_INITED))
504 returnErr(TM_E_ILLEGAL);
505 if (tms == tmTop) /* check necessity */
506 returnOK;
507 /* pull if already in stack */
508 (void)tmPull(tms);
509 /* push it on top */
510 tms->tmprev = tmTop;
511 tmTop = tms;
512 returnOK;
513 }
514
515
516 void
517 tmDone(tms) /* done with tone mapping -- destroy it */
518 register struct tmStruct *tms;
519 {
520 /* NULL arg. is equiv. to tmTop */
521 if (tms == NULL && (tms = tmTop) == NULL)
522 return;
523 /* take out of stack if present */
524 (void)tmPull(tms);
525 /* free tables */
526 if (tms->histo != NULL)
527 free((char *)tms->histo);
528 if (tms->lumap != NULL)
529 free((char *)tms->lumap);
530 tms->flags = 0;
531 free((char *)tms); /* free basic structure */
532 }