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

# User Rev Content
1 greg 3.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 greg 3.4 /* our list of conversion packages */
22     struct tmPackage *tmPkg[TM_MAXPKG];
23     int tmNumPkgs = 0; /* number of registered packages */
24 greg 3.1
25 gregl 3.5 int tmLastError; /* last error incurred by library */
26     char *tmLastFunction; /* error-generating function name */
27 greg 3.4
28 gregl 3.5
29 greg 3.1 int
30     tmErrorReturn(func, err) /* error return (with message) */
31     char *func;
32     int err;
33     {
34 gregl 3.5 tmLastFunction = func;
35     tmLastError = err;
36 greg 3.2 if (tmTop != NULL && tmTop->flags & TM_F_NOSTDERR)
37 greg 3.1 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 greg 3.4 /* set color divisors */
80     for (i = 0; i < 3; i++)
81     tmnew->cdiv[i] = 256.*pow(tmnew->clf[i], 1./tmnew->mongam);
82    
83 greg 3.1 /* 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 greg 3.4 /* zero private data */
93     for (i = TM_MAXPKG; i--; )
94     tmnew->pd[i] = NULL;
95     /* make tmnew current */
96 greg 3.1 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 greg 3.4 /* 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 greg 3.1 returnOK;
159     }
160    
161    
162     void
163     tmClearHisto() /* clear current histogram */
164     {
165     if (tmTop == NULL || tmTop->histo == NULL)
166     return;
167 greg 3.4 free((MEM_PTR)tmTop->histo);
168 greg 3.1 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 greg 3.4 static COLOR csmall = {1e-6, 1e-6, 1e-6};
181 greg 3.1 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 greg 3.4 if (tmNeedMatrix(tmTop)) /* get monitor RGB */
192 greg 3.1 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 greg 3.4 if (clipgamut(cmon, lum, CGAMUT_LOWER, csmall, cwhite))
204 greg 3.1 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 greg 3.4 cs[3*i ] = d>=.999 ? 255 :
237     (int)(256.*pow(d, 1./tmTop->mongam));
238 greg 3.1 d = tmTop->clf[GRN]*cmon[GRN]/lum;
239 greg 3.4 cs[3*i+1] = d>=.999 ? 255 :
240     (int)(256.*pow(d, 1./tmTop->mongam));
241 greg 3.1 d = tmTop->clf[BLU]*cmon[BLU]/lum;
242 greg 3.4 cs[3*i+2] = d>=.999 ? 255 :
243     (int)(256.*pow(d, 1./tmTop->mongam));
244 greg 3.1 }
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 greg 3.2 for (i = oldlen, j = i+oldorig-horig; i; )
291     newhist[--j] = tmTop->histo[--i];
292 greg 3.4 free((MEM_PTR)tmTop->histo);
293 greg 3.1 }
294     tmTop->histo = newhist;
295     if (tmTop->lumap != NULL) { /* invalid tone map */
296 greg 3.4 free((MEM_PTR)tmTop->lumap);
297 greg 3.1 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 greg 3.2 histot += tmTop->histo[i];
364     sum += (j -= HISTEP) * tmTop->histo[i];
365 greg 3.1 }
366     threshold = histot*.025 + .5;
367     if (threshold < 4)
368     returnErr(TM_E_TMFAIL);
369     Lwavg = tmLuminance( (double)sum / histot );
370 greg 3.2 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 greg 3.1 ((double)histlen*TM_BRTSCALE) / logLddyn;
386 greg 3.2 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 greg 3.1 .5*(cumf[i]+cumf[i+1]) );
393 greg 3.2 ceiling = Tr * (htcontrs(Ld) * Lw) /
394 greg 3.1 (htcontrs(Lw) * Ld) + 1.;
395 greg 3.2 }
396     if (histo[i] > ceiling) {
397     trimmings += histo[i] - ceiling;
398     histo[i] = ceiling;
399     }
400 greg 3.1 }
401 greg 3.2 } while ((histot -= trimmings) > threshold &&
402     trimmings > threshold);
403     }
404     /* allocate luminance map */
405     if (tmTop->lumap == NULL) {
406 greg 3.1 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 greg 3.2 if (!(tmTop->flags & TM_F_LINEAR)) {
433 greg 3.4 free((MEM_PTR)histo);
434     free((MEM_PTR)cumf);
435 greg 3.2 }
436 greg 3.1 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 greg 3.4 register int4 li, pv;
449 greg 3.1
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 greg 3.4 pv = *cs++ * li / tmTop->cdiv[RED];
464 greg 3.1 *ps++ = pv>255 ? 255 : pv;
465 greg 3.4 pv = *cs++ * li / tmTop->cdiv[GRN];
466 greg 3.1 *ps++ = pv>255 ? 255 : pv;
467 greg 3.4 pv = *cs++ * li / tmTop->cdiv[BLU];
468 greg 3.1 *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 greg 3.3 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 greg 3.4 /* clear package data */
539     for (i = tmNumPkgs; i--; )
540     tmnew->pd[i] = NULL;
541 greg 3.3 tmnew->tmprev = tmTop; /* make copy current */
542     return(tmTop = tmnew);
543     }
544    
545    
546 greg 3.1 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 greg 3.4 if (tms == NULL)
553 greg 3.1 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 greg 3.4 register int i;
570 greg 3.1 /* 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 greg 3.4 free((MEM_PTR)tms->histo);
578 greg 3.1 if (tms->lumap != NULL)
579 greg 3.4 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 greg 3.1 }