ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/tonemap.c
Revision: 3.18
Committed: Fri Jan 7 22:05:30 2005 UTC (19 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R7P2, rad3R7P1
Changes since 3.17: +6 -3 lines
Log Message:
Added client data pointer to tmSetSpace() call

File Contents

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