ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/tonemap.c
Revision: 3.4
Committed: Fri Apr 18 13:59:43 1997 UTC (27 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.3: +53 -45 lines
Log Message:
added notion of conversion packages with private data

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