ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/tonemap.c
Revision: 3.5
Committed: Wed Nov 12 13:25:48 1997 UTC (26 years, 5 months ago) by gregl
Content type: text/plain
Branch: MAIN
Changes since 3.4: +5 -0 lines
Log Message:
added tmLastError and tmLastFunction globals (set by tmErrorReturn)

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