ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/tonemap.c
Revision: 3.1
Committed: Tue Apr 15 16:53:01 1997 UTC (27 years ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Initial revision

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     if (tmTop != NULL && tmTop->flags & TM_F_NOERRS)
28     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     for (i = oldlen, j = i+oldorig-horig; i--; )
281     newhist[--j] = tmTop->histo[i];
282     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     /* allocate temporary tables */
350     histo = (int *)malloc(histlen*sizeof(int));
351     cumf = (float *)malloc((histlen+1)*sizeof(float));
352     if (histo == NULL | cumf == NULL)
353     returnErr(TM_E_NOMEM);
354     /* histogram total and mean */
355     histot = 0; sum = 0;
356     j = brt0 + histlen*HISTEP;
357     for (i = histlen; i--; ) {
358     histot += (histo[i] = tmTop->histo[i]);
359     sum += (j -= HISTEP) * histo[i];
360     }
361     threshold = histot*.025 + .5;
362     if (threshold < 4)
363     returnErr(TM_E_TMFAIL);
364     Lwavg = tmLuminance( (double)sum / histot );
365     do { /* iterate to solution */
366     sum = 0; /* compute cumulative probability */
367     for (i = 0; i < histlen; i++) {
368     cumf[i] = (double)sum/histot;
369     sum += histo[i];
370     }
371     cumf[i] = 1.;
372     Tr = histot * (double)(tmTop->brmax - tmTop->brmin) /
373     ((double)histlen*TM_BRTSCALE) / logLddyn;
374     ceiling = Tr + 1.;
375     trimmings = 0; /* clip to envelope */
376     for (i = histlen; i--; ) {
377     if (tmTop->flags & TM_F_HCONTR) {
378     Lw = tmLuminance(brt0 + i*HISTEP);
379     Ld = Ldmin * exp( logLddyn *
380     .5*(cumf[i]+cumf[i+1]) );
381     ceiling = Tr * (htcontrs(Ld) * Lw) /
382     (htcontrs(Lw) * Ld) + 1.;
383     }
384     if (histo[i] > ceiling) {
385     trimmings += histo[i] - ceiling;
386     histo[i] = ceiling;
387     }
388     }
389     } while ((histot -= trimmings) > threshold && trimmings > threshold);
390    
391     if (tmTop->lumap == NULL) { /* allocate luminance map */
392     tmTop->lumap = (unsigned short *)malloc(
393     (tmTop->brmax-tmTop->brmin+1)*sizeof(unsigned short) );
394     if (tmTop->lumap == NULL)
395     returnErr(TM_E_NOMEM);
396     }
397     if (tmTop->flags & TM_F_LINEAR || histot <= threshold) {
398     /* linear tone mapping */
399     if (tmTop->flags & TM_F_HCONTR)
400     d = htcontrs(Ldavg) / htcontrs(Lwavg);
401     else
402     d = Ldavg / Lwavg;
403     d = log(d/Ldmax);
404     for (i = tmTop->brmax-tmTop->brmin+1; i--; )
405     tmTop->lumap[i] = 256. * exp(
406     ( d + (tmTop->brmin+i)/(double)TM_BRTSCALE )
407     / gamval );
408     } else {
409     /* histogram adjustment */
410     for (i = tmTop->brmax-tmTop->brmin+1; i--; ) {
411     j = d = (double)i/(tmTop->brmax-tmTop->brmin)*histlen;
412     d -= (double)j;
413     Ld = Ldmin*exp(logLddyn*((1.-d)*cumf[j]+d*cumf[j+1]));
414     d = (Ld - Ldmin)/(Ldmax - Ldmin);
415     tmTop->lumap[i] = 256.*pow(d, 1./gamval);
416     }
417     }
418     free((char *)histo);
419     free((char *)cumf);
420     returnOK;
421     }
422    
423    
424     int
425     tmMapPixels(ps, ls, cs, len)
426     register BYTE *ps;
427     TMbright *ls;
428     register BYTE *cs;
429     int len;
430     {
431     static char funcName[] = "tmMapPixels";
432     int rdiv, gdiv, bdiv;
433     register int li, pv;
434    
435     if (tmTop == NULL || tmTop->lumap == NULL)
436     returnErr(TM_E_TMINVAL);
437     if (ps == NULL | ls == NULL | len <= 0)
438     returnErr(TM_E_ILLEGAL);
439     rdiv = tmTop->gamb[((int4)TM_GAMTSZ*tmTop->clfb[RED])>>8];
440     gdiv = tmTop->gamb[((int4)TM_GAMTSZ*tmTop->clfb[GRN])>>8];
441     bdiv = tmTop->gamb[((int4)TM_GAMTSZ*tmTop->clfb[BLU])>>8];
442     while (len--) {
443     if ((li = *ls++) < tmTop->brmin)
444     li = tmTop->brmin;
445     else if (li > tmTop->brmax)
446     li = tmTop->brmax;
447     li = tmTop->lumap[li - tmTop->brmin];
448     if (cs == TM_NOCHROM)
449     *ps++ = li>255 ? 255 : li;
450     else {
451     pv = *cs++ * li / rdiv;
452     *ps++ = pv>255 ? 255 : pv;
453     pv = *cs++ * li / gdiv;
454     *ps++ = pv>255 ? 255 : pv;
455     pv = *cs++ * li / bdiv;
456     *ps++ = pv>255 ? 255 : pv;
457     }
458     }
459     returnOK;
460     }
461    
462    
463     struct tmStruct *
464     tmPop() /* pop top tone mapping off stack */
465     {
466     register struct tmStruct *tms;
467    
468     if ((tms = tmTop) != NULL)
469     tmTop = tms->tmprev;
470     return(tms);
471     }
472    
473    
474     int
475     tmPull(tms) /* pull a tone mapping from stack */
476     register struct tmStruct *tms;
477     {
478     register struct tmStruct *tms2;
479     /* special cases first */
480     if (tms == NULL | tmTop == NULL)
481     return(0);
482     if (tms == tmTop) {
483     tmTop = tms->tmprev;
484     tms->tmprev = NULL;
485     return(1);
486     }
487     for (tms2 = tmTop; tms2->tmprev != NULL; tms2 = tms2->tmprev)
488     if (tms == tms2->tmprev) { /* remove it */
489     tms2->tmprev = tms->tmprev;
490     tms->tmprev = NULL;
491     return(1);
492     }
493     return(0); /* not found on stack */
494     }
495    
496    
497     int
498     tmPush(tms) /* push tone mapping on top of stack */
499     register struct tmStruct *tms;
500     {
501     static char funcName[] = "tmPush";
502     /* check validity */
503     if (tms == NULL || !(tms->flags & TM_F_INITED))
504     returnErr(TM_E_ILLEGAL);
505     if (tms == tmTop) /* check necessity */
506     returnOK;
507     /* pull if already in stack */
508     (void)tmPull(tms);
509     /* push it on top */
510     tms->tmprev = tmTop;
511     tmTop = tms;
512     returnOK;
513     }
514    
515    
516     void
517     tmDone(tms) /* done with tone mapping -- destroy it */
518     register struct tmStruct *tms;
519     {
520     /* NULL arg. is equiv. to tmTop */
521     if (tms == NULL && (tms = tmTop) == NULL)
522     return;
523     /* take out of stack if present */
524     (void)tmPull(tms);
525     /* free tables */
526     if (tms->histo != NULL)
527     free((char *)tms->histo);
528     if (tms->lumap != NULL)
529     free((char *)tms->lumap);
530     tms->flags = 0;
531     free((char *)tms); /* free basic structure */
532     }