ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/tonemap.c
Revision: 3.3
Committed: Wed Apr 16 20:28:03 1997 UTC (27 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.2: +34 -0 lines
Log Message:
minor bug fixes and added tmDup()

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