ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/tonemap.c
(Generate patch)

Comparing ray/src/common/tonemap.c (file contents):
Revision 3.8 by gwlarson, Mon Aug 17 17:58:47 1998 UTC vs.
Revision 3.12 by greg, Tue Jul 15 00:30:09 2003 UTC

# Line 1 | Line 1
1 /* Copyright (c) 1998 Silicon Graphics, Inc. */
2
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ SGI";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6
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"
# Line 26 | Line 28 | int    tmLastError;                    /* last error incurred by library *
28   char    *tmLastFunction;                /* error-generating function name */
29  
30  
29 int
30 tmErrorReturn(func, err)                /* error return (with message) */
31 char    *func;
32 int     err;
33 {
34        tmLastFunction = func;
35        tmLastError = err;
36        if (tmTop != NULL && tmTop->flags & TM_F_NOSTDERR)
37                return(err);
38        fputs(func, stderr);
39        fputs(": ", stderr);
40        fputs(tmErrorMessage[err], stderr);
41        fputs("!\n", stderr);
42        return(err);
43 }
44
45
31   struct tmStruct *
32   tmInit(flags, monpri, gamval)           /* initialize new tone mapping */
33   int     flags;
34   RGBPRIMP        monpri;
35   double  gamval;
36   {
52        static char     funcName[] = "tmInit";
37          COLORMAT        cmat;
38          register struct tmStruct        *tmnew;
39          register int    i;
# Line 66 | Line 50 | double gamval;
50                  tmnew->clf[GRN] = rgb2xyzmat[1][1];
51                  tmnew->clf[BLU] = rgb2xyzmat[1][2];
52          } else {
53 <                comprgb2xyzmat(cmat, tmnew->monpri=monpri);
53 >                comprgb2xyzWBmat(cmat, tmnew->monpri=monpri);
54                  tmnew->clf[RED] = cmat[1][0];
55                  tmnew->clf[GRN] = cmat[1][1];
56                  tmnew->clf[BLU] = cmat[1][2];
# Line 86 | Line 70 | double gamval;
70                          tmnew->inpsf = WHTEFFICACY;
71          tmnew->cmat[0][1] = tmnew->cmat[0][2] = tmnew->cmat[1][0] =
72          tmnew->cmat[1][2] = tmnew->cmat[2][0] = tmnew->cmat[2][1] = 0.;
73 <        tmnew->hbrmin = tmnew->hbrmax = 0;
73 >        tmnew->hbrmin = 10; tmnew->hbrmax = -10;
74          tmnew->histo = NULL;
75 <        tmnew->mbrmin = tmnew->mbrmax = 0;
75 >        tmnew->mbrmin = 10; tmnew->mbrmax = -10;
76          tmnew->lumap = NULL;
77                                                  /* zero private data */
78          for (i = TM_MAXPKG; i--; )
# Line 123 | Line 107 | double sf;
107                          tmTop->clf[CIEX] = tmTop->clf[CIEZ] = 0.;
108                          tmTop->clf[CIEY] = 1.;
109                  } else {
110 <                        comprgb2xyzmat(tmTop->cmat, tmTop->monpri);
110 >                        comprgb2xyzWBmat(tmTop->cmat, tmTop->monpri);
111                          tmTop->clf[RED] = tmTop->cmat[1][0];
112                          tmTop->clf[GRN] = tmTop->cmat[1][1];
113                          tmTop->clf[BLU] = tmTop->cmat[1][2];
# Line 134 | Line 118 | double sf;
118                  tmTop->cmat[1][2] = tmTop->cmat[2][0] = tmTop->cmat[2][1] = 0.;
119  
120          } else if (tmTop->inppri == TM_XYZPRIM) /* input is XYZ */
121 <                compxyz2rgbmat(tmTop->cmat, tmTop->monpri);
121 >                compxyz2rgbWBmat(tmTop->cmat, tmTop->monpri);
122  
123          else {                                  /* input is RGB */
124                  if (tmTop->inppri != tmTop->monpri &&
125                                  PRIMEQ(tmTop->inppri, tmTop->monpri))
126                          tmTop->inppri = tmTop->monpri;  /* no xform */
127 <                comprgb2rgbmat(tmTop->cmat, tmTop->inppri, tmTop->monpri);
127 >                comprgb2rgbWBmat(tmTop->cmat, tmTop->inppri, tmTop->monpri);
128          }
129          for (i = 0; i < 3; i++)
130                  for (j = 0; j < 3; j++)
# Line 178 | Line 162 | COLOR  *scan;
162   int     len;
163   {
164          static char     funcName[] = "tmCvColors";
165 <        static COLOR    csmall = {1e-6, 1e-6, 1e-6};
165 >        static COLOR    csmall = {.5*MINLUM, .5*MINLUM, .5*MINLUM};
166          COLOR   cmon;
167          double  lum, slum;
168          register double d;
# Line 248 | Line 232 | int    len;
232  
233  
234   int
235 + tmCvGrays(ls, scan, len)                /* convert float gray values */
236 + TMbright        *ls;
237 + float   *scan;
238 + int     len;
239 + {
240 +        static char     funcName[] = "tmCvGrays";
241 +        register double d;
242 +        register int    i;
243 +
244 +        if (tmTop == NULL)
245 +                returnErr(TM_E_TMINVAL);
246 +        if (ls == NULL | scan == NULL | len < 0)
247 +                returnErr(TM_E_ILLEGAL);
248 +        for (i = len; i--; )
249 +                if (scan[i] <= TM_NOLUM) {
250 +                        ls[i] = TM_NOBRT;               /* bogus value */
251 +                } else {
252 +                        d = TM_BRTSCALE*log(scan[i]);   /* encode it */
253 +                        ls[i] = d>0. ? (int)(d+.5) : (int)(d-.5);
254 +                }
255 +        returnOK;
256 + }
257 +
258 +
259 + int
260   tmAddHisto(ls, len, wt)                 /* add values to histogram */
261   register TMbright       *ls;
262   int     len;
263   int     wt;
264   {
265          static char     funcName[] = "tmAddHisto";
266 <        int     oldorig, oldlen, horig, hlen;
266 >        int     oldorig=0, oldlen, horig, hlen;
267          register int    i, j;
268  
269          if (tmTop == NULL)
# Line 327 | Line 336 | double La;
336   }
337  
338  
339 + static int
340 + tmNewMap()
341 + {
342 +        if (tmTop->lumap != NULL && (tmTop->mbrmax - tmTop->mbrmin) !=
343 +                                        (tmTop->hbrmax - tmTop->hbrmin)) {
344 +                free((MEM_PTR)tmTop->lumap);
345 +                tmTop->lumap = NULL;
346 +        }
347 +        tmTop->mbrmin = tmTop->hbrmin;
348 +        tmTop->mbrmax = tmTop->hbrmax;
349 +        if (tmTop->mbrmin > tmTop->mbrmax)
350 +                return 0;
351 +        if (tmTop->lumap == NULL)
352 +                tmTop->lumap = (unsigned short *)malloc(sizeof(unsigned short)*
353 +                                        (tmTop->mbrmax-tmTop->mbrmin+1));
354 +        return(tmTop->lumap != NULL);
355 + }
356 +
357 +
358   int
359 + tmFixedMapping(expmult, gamval)
360 + double  expmult;
361 + double  gamval;
362 + {
363 +        static char     funcName[] = "tmFixedMapping";
364 +        double          d;
365 +        register int    i;
366 +        
367 +        if (!tmNewMap())
368 +                returnErr(TM_E_NOMEM);
369 +        if (expmult <= .0)
370 +                expmult = 1.;
371 +        if (gamval < MINGAM)
372 +                gamval = tmTop->mongam;
373 +        d = log(expmult/tmTop->inpsf);
374 +        for (i = tmTop->mbrmax-tmTop->mbrmin+1; i--; )
375 +                tmTop->lumap[i] = 256. * exp(
376 +                        ( d + (tmTop->mbrmin+i)*(1./TM_BRTSCALE) )
377 +                        / gamval );
378 +        returnOK;
379 + }
380 +
381 +
382 + int
383   tmComputeMapping(gamval, Lddyn, Ldmax)
384   double  gamval;
385   double  Lddyn;
# Line 336 | Line 388 | double Ldmax;
388          static char     funcName[] = "tmComputeMapping";
389          int     *histo;
390          float   *cumf;
391 <        int     brt0, histlen, histot, threshold, ceiling, trimmings;
391 >        int     brt0, histlen, threshold, ceiling, trimmings;
392          double  logLddyn, Ldmin, Ldavg, Lwavg, Tr, Lw, Ld;
393 <        int4    sum;
393 >        int32   histot;
394 >        double  sum;
395          register double d;
396          register int    i, j;
397  
# Line 366 | Line 419 | double Ldmax;
419          if (threshold < 4)
420                  returnErr(TM_E_TMFAIL);
421          Lwavg = tmLuminance( (double)sum / histot );
422 <        if (!(tmTop->flags & TM_F_LINEAR)) {    /* clamp histogram */
423 <                histo = (int *)malloc(histlen*sizeof(int));
424 <                cumf = (float *)malloc((histlen+1)*sizeof(float));
425 <                if (histo == NULL | cumf == NULL)
426 <                        returnErr(TM_E_NOMEM);
427 <                for (i = histlen; i--; )        /* make malleable copy */
428 <                        histo[i] = tmTop->histo[i];
429 <                do {                            /* iterate to solution */
430 <                        sum = 0;                /* cumulative probability */
431 <                        for (i = 0; i < histlen; i++) {
432 <                                cumf[i] = (double)sum/histot;
433 <                                sum += histo[i];
422 >                                        /* allocate space for mapping */
423 >        if (!tmNewMap())
424 >                returnErr(TM_E_NOMEM);
425 >                                        /* use linear tone mapping? */
426 >        if (tmTop->flags & TM_F_LINEAR)
427 >                goto linearmap;
428 >                                        /* clamp histogram */
429 >        histo = (int *)malloc(histlen*sizeof(int));
430 >        cumf = (float *)malloc((histlen+2)*sizeof(float));
431 >        if (histo == NULL | cumf == NULL)
432 >                returnErr(TM_E_NOMEM);
433 >        cumf[histlen+1] = 1.;           /* guard for assignment code */
434 >        for (i = histlen; i--; )        /* make malleable copy */
435 >                histo[i] = tmTop->histo[i];
436 >        do {                            /* iterate to solution */
437 >                sum = 0;                /* cumulative probability */
438 >                for (i = 0; i < histlen; i++) {
439 >                        cumf[i] = (double)sum/histot;
440 >                        sum += histo[i];
441 >                }
442 >                cumf[histlen] = 1.;
443 >                Tr = histot * (double)(tmTop->hbrmax - tmTop->hbrmin) /
444 >                        ((double)histlen*TM_BRTSCALE) / logLddyn;
445 >                ceiling = Tr + 1.;
446 >                trimmings = 0;          /* clip to envelope */
447 >                for (i = histlen; i--; ) {
448 >                        if (tmTop->flags & TM_F_HCONTR) {
449 >                                Lw = tmLuminance(brt0 + i*HISTEP);
450 >                                Ld = Ldmin * exp( logLddyn *
451 >                                        .5*(cumf[i]+cumf[i+1]) );
452 >                                ceiling = Tr * (htcontrs(Ld) * Lw) /
453 >                                        (htcontrs(Lw) * Ld) + 1.;
454                          }
455 <                        cumf[i] = 1.;
456 <                        Tr = histot * (double)(tmTop->hbrmax - tmTop->hbrmin) /
457 <                                ((double)histlen*TM_BRTSCALE) / logLddyn;
385 <                        ceiling = Tr + 1.;
386 <                        trimmings = 0;                  /* clip to envelope */
387 <                        for (i = histlen; i--; ) {
388 <                                if (tmTop->flags & TM_F_HCONTR) {
389 <                                        Lw = tmLuminance(brt0 + i*HISTEP);
390 <                                        Ld = Ldmin * exp( logLddyn *
391 <                                                .5*(cumf[i]+cumf[i+1]) );
392 <                                        ceiling = Tr * (htcontrs(Ld) * Lw) /
393 <                                                (htcontrs(Lw) * Ld) + 1.;
394 <                                }
395 <                                if (histo[i] > ceiling) {
396 <                                        trimmings += histo[i] - ceiling;
397 <                                        histo[i] = ceiling;
398 <                                }
455 >                        if (histo[i] > ceiling) {
456 >                                trimmings += histo[i] - ceiling;
457 >                                histo[i] = ceiling;
458                          }
400                } while ((histot -= trimmings) > threshold &&
401                                                trimmings > threshold);
402        }
403                                                /* allocate luminance map */
404        if (tmTop->lumap != NULL)
405                free((MEM_PTR)tmTop->lumap);
406        tmTop->mbrmin = tmTop->hbrmin;
407        tmTop->mbrmax = tmTop->hbrmax;
408        tmTop->lumap = (unsigned short *)malloc(
409                (tmTop->mbrmax-tmTop->mbrmin+1)*sizeof(unsigned short) );
410        if (tmTop->lumap == NULL)
411                returnErr(TM_E_NOMEM);
412        if (tmTop->flags & TM_F_LINEAR || histot <= threshold) {
413                                                /* linear tone mapping */
414                if (tmTop->flags & TM_F_HCONTR)
415                        d = htcontrs(Ldavg) / htcontrs(Lwavg);
416                else
417                        d = Ldavg / Lwavg;
418                d = log(d/Ldmax);
419                for (i = tmTop->mbrmax-tmTop->mbrmin+1; i--; )
420                        tmTop->lumap[i] = 256. * exp(
421                                ( d + (tmTop->mbrmin+i)/(double)TM_BRTSCALE )
422                                / gamval );
423        } else {
424                                                /* histogram adjustment */
425                for (i = tmTop->mbrmax-tmTop->mbrmin+1; i--; ) {
426                        j = d = (double)i/(tmTop->mbrmax-tmTop->mbrmin)*histlen;
427                        d -= (double)j;
428                        Ld = Ldmin*exp(logLddyn*((1.-d)*cumf[j]+d*cumf[j+1]));
429                        d = (Ld - Ldmin)/(Ldmax - Ldmin);
430                        tmTop->lumap[i] = 256.*pow(d, 1./gamval);
459                  }
460 +                                        /* check if we're out of data */
461 +                if ((histot -= trimmings) <= threshold) {
462 +                        free((MEM_PTR)histo);
463 +                        free((MEM_PTR)cumf);
464 +                        goto linearmap;
465 +                }
466 +        } while (trimmings > threshold);
467 +                                        /* assign tone-mapping */
468 +        for (i = tmTop->mbrmax-tmTop->mbrmin+1; i--; ) {
469 +                j = d = (double)i/(tmTop->mbrmax-tmTop->mbrmin)*histlen;
470 +                d -= (double)j;
471 +                Ld = Ldmin*exp(logLddyn*((1.-d)*cumf[j]+d*cumf[j+1]));
472 +                d = (Ld - Ldmin)/(Ldmax - Ldmin);
473 +                tmTop->lumap[i] = 256.*pow(d, 1./gamval);
474          }
475 <        if (!(tmTop->flags & TM_F_LINEAR)) {
476 <                free((MEM_PTR)histo);
435 <                free((MEM_PTR)cumf);
436 <        }
475 >        free((MEM_PTR)histo);           /* clean up and return */
476 >        free((MEM_PTR)cumf);
477          returnOK;
478 + linearmap:                              /* linear tone-mapping */
479 +        if (tmTop->flags & TM_F_HCONTR)
480 +                d = htcontrs(Ldavg) / htcontrs(Lwavg);
481 +        else
482 +                d = Ldavg / Lwavg;
483 +        return(tmFixedMapping(tmTop->inpsf*d/Ldmax, gamval));
484   }
485  
486  
# Line 446 | Line 492 | register BYTE  *cs;
492   int     len;
493   {
494          static char     funcName[] = "tmMapPixels";
495 <        register int4   li, pv;
495 >        register int32  li, pv;
496  
497          if (tmTop == NULL || tmTop->lumap == NULL)
498                  returnErr(TM_E_TMINVAL);
499          if (ps == NULL | ls == NULL | len < 0)
500                  returnErr(TM_E_ILLEGAL);
501          while (len--) {
502 <                if ((li = *ls++) < tmTop->mbrmin)
503 <                        li = tmTop->mbrmin;
504 <                else if (li > tmTop->mbrmax)
505 <                        li = tmTop->mbrmax;
506 <                li = tmTop->lumap[li - tmTop->mbrmin];
502 >                if ((li = *ls++) < tmTop->mbrmin) {
503 >                        li = 0;
504 >                } else {
505 >                        if (li > tmTop->mbrmax)
506 >                                li = tmTop->mbrmax;
507 >                        li = tmTop->lumap[li - tmTop->mbrmin];
508 >                }
509                  if (cs == TM_NOCHROM)
510                          *ps++ = li>255 ? 255 : li;
511                  else {
# Line 583 | Line 631 | register struct tmStruct       *tms;
631                  if (tms->pd[i] != NULL)
632                          (*tmPkg[i]->Free)(tms->pd[i]);
633          free((MEM_PTR)tms);             /* free basic structure */
634 + }
635 +
636 + /******************** Shared but Private library routines *********************/
637 +
638 + BYTE    tmMesofact[BMESUPPER-BMESLOWER];
639 +
640 + void
641 + tmMkMesofact()                          /* build mesopic lookup factor table */
642 + {
643 +        register int    i;
644 +
645 +        if (tmMesofact[BMESUPPER-BMESLOWER-1])
646 +                return;
647 +
648 +        for (i = BMESLOWER; i < BMESUPPER; i++)
649 +                tmMesofact[i-BMESLOWER] = 256. *
650 +                                (tmLuminance(i) - LMESLOWER) /
651 +                                (LMESUPPER - LMESLOWER);
652 + }
653 +
654 +
655 + int
656 + tmErrorReturn(func, err)                /* error return (with message) */
657 + char    *func;
658 + int     err;
659 + {
660 +        tmLastFunction = func;
661 +        tmLastError = err;
662 +        if (tmTop != NULL && tmTop->flags & TM_F_NOSTDERR)
663 +                return(err);
664 +        fputs(func, stderr);
665 +        fputs(": ", stderr);
666 +        fputs(tmErrorMessage[err], stderr);
667 +        fputs("!\n", stderr);
668 +        return(err);
669   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines