ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/tonemap.c
Revision: 3.17
Committed: Fri Jan 7 21:41:06 2005 UTC (19 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.16: +17 -16 lines
Log Message:
Moved error codes inside passed TMstruct

File Contents

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