ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/tonemap.c
Revision: 3.4
Committed: Fri Apr 18 13:59:43 1997 UTC (27 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.3: +53 -45 lines
Log Message:
added notion of conversion packages with private data

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