ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/tonemap.c
Revision: 3.16
Committed: Fri Jan 7 20:33:02 2005 UTC (19 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.15: +211 -243 lines
Log Message:
Modernized tone-mapping routines with structure pointer r.t. stack

File Contents

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