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.7 by gregl, Mon Dec 1 09:52:00 1997 UTC vs.
Revision 3.15 by greg, Thu Jul 22 17:48:52 2004 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1997 Regents of the University of California */
2 –
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ LBL";
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 59 | Line 43 | double gamval;
43                  return(NULL);
44  
45          tmnew->flags = flags & ~TM_F_UNIMPL;
46 +        if (tmnew->flags & TM_F_BW)
47 +                tmnew->flags &= ~TM_F_MESOPIC;
48                                                  /* set monitor transform */
49          if (monpri == NULL || monpri == stdprims || tmnew->flags & TM_F_BW) {
50                  tmnew->monpri = stdprims;
# Line 66 | Line 52 | double gamval;
52                  tmnew->clf[GRN] = rgb2xyzmat[1][1];
53                  tmnew->clf[BLU] = rgb2xyzmat[1][2];
54          } else {
55 <                comprgb2xyzmat(cmat, tmnew->monpri=monpri);
55 >                comprgb2xyzWBmat(cmat, tmnew->monpri=monpri);
56                  tmnew->clf[RED] = cmat[1][0];
57                  tmnew->clf[GRN] = cmat[1][1];
58                  tmnew->clf[BLU] = cmat[1][2];
# Line 86 | Line 72 | double gamval;
72                          tmnew->inpsf = WHTEFFICACY;
73          tmnew->cmat[0][1] = tmnew->cmat[0][2] = tmnew->cmat[1][0] =
74          tmnew->cmat[1][2] = tmnew->cmat[2][0] = tmnew->cmat[2][1] = 0.;
75 <        tmnew->hbrmin = tmnew->hbrmax = 0;
75 >        tmnew->hbrmin = 10; tmnew->hbrmax = -10;
76          tmnew->histo = NULL;
77 <        tmnew->mbrmin = tmnew->mbrmax = 0;
77 >        tmnew->mbrmin = 10; tmnew->mbrmax = -10;
78          tmnew->lumap = NULL;
79                                                  /* zero private data */
80          for (i = TM_MAXPKG; i--; )
# Line 123 | Line 109 | double sf;
109                          tmTop->clf[CIEX] = tmTop->clf[CIEZ] = 0.;
110                          tmTop->clf[CIEY] = 1.;
111                  } else {
112 <                        comprgb2xyzmat(tmTop->cmat, tmTop->monpri);
112 >                        comprgb2xyzWBmat(tmTop->cmat, tmTop->monpri);
113                          tmTop->clf[RED] = tmTop->cmat[1][0];
114                          tmTop->clf[GRN] = tmTop->cmat[1][1];
115                          tmTop->clf[BLU] = tmTop->cmat[1][2];
# Line 134 | Line 120 | double sf;
120                  tmTop->cmat[1][2] = tmTop->cmat[2][0] = tmTop->cmat[2][1] = 0.;
121  
122          } else if (tmTop->inppri == TM_XYZPRIM) /* input is XYZ */
123 <                compxyz2rgbmat(tmTop->cmat, tmTop->monpri);
123 >                compxyz2rgbWBmat(tmTop->cmat, tmTop->monpri);
124  
125          else {                                  /* input is RGB */
126                  if (tmTop->inppri != tmTop->monpri &&
127                                  PRIMEQ(tmTop->inppri, tmTop->monpri))
128                          tmTop->inppri = tmTop->monpri;  /* no xform */
129 <                comprgb2rgbmat(tmTop->cmat, tmTop->inppri, tmTop->monpri);
129 >                comprgb2rgbWBmat(tmTop->cmat, tmTop->inppri, tmTop->monpri);
130          }
131          for (i = 0; i < 3; i++)
132                  for (j = 0; j < 3; j++)
# Line 178 | Line 164 | COLOR  *scan;
164   int     len;
165   {
166          static char     funcName[] = "tmCvColors";
167 <        static COLOR    csmall = {1e-6, 1e-6, 1e-6};
167 >        static COLOR    csmall = {.5*MINLUM, .5*MINLUM, .5*MINLUM};
168          COLOR   cmon;
169          double  lum, slum;
170          register double d;
# Line 186 | Line 172 | int    len;
172  
173          if (tmTop == NULL)
174                  returnErr(TM_E_TMINVAL);
175 <        if (ls == NULL | scan == NULL | len <= 0)
175 >        if ((ls == NULL) | (scan == NULL) | (len < 0))
176                  returnErr(TM_E_ILLEGAL);
177          for (i = len; i--; ) {
178 <                if (tmNeedMatrix(tmTop))                /* get monitor RGB */
178 >                if (tmNeedMatrix(tmTop)) {              /* get monitor RGB */
179                          colortrans(cmon, tmTop->cmat, scan[i]);
180 <                else {
180 >                } else {
181                          cmon[RED] = tmTop->inpsf*scan[i][RED];
182                          cmon[GRN] = tmTop->inpsf*scan[i][GRN];
183                          cmon[BLU] = tmTop->inpsf*scan[i][BLU];
# Line 248 | Line 234 | int    len;
234  
235  
236   int
237 + tmCvGrays(ls, scan, len)                /* convert float gray values */
238 + TMbright        *ls;
239 + float   *scan;
240 + int     len;
241 + {
242 +        static char     funcName[] = "tmCvGrays";
243 +        register double d;
244 +        register int    i;
245 +
246 +        if (tmTop == NULL)
247 +                returnErr(TM_E_TMINVAL);
248 +        if ((ls == NULL) | (scan == NULL) | (len < 0))
249 +                returnErr(TM_E_ILLEGAL);
250 +        for (i = len; i--; )
251 +                if (scan[i] <= TM_NOLUM) {
252 +                        ls[i] = TM_NOBRT;               /* bogus value */
253 +                } else {
254 +                        d = TM_BRTSCALE*log(scan[i]);   /* encode it */
255 +                        ls[i] = d>0. ? (int)(d+.5) : (int)(d-.5);
256 +                }
257 +        returnOK;
258 + }
259 +
260 +
261 + int
262   tmAddHisto(ls, len, wt)                 /* add values to histogram */
263   register TMbright       *ls;
264   int     len;
265   int     wt;
266   {
267          static char     funcName[] = "tmAddHisto";
268 <        int     oldorig, oldlen, horig, hlen;
268 >        int     oldorig=0, oldlen, horig, hlen;
269          register int    i, j;
270  
260 –        if (len <= 0)
261 –                returnErr(TM_E_ILLEGAL);
271          if (tmTop == NULL)
272                  returnErr(TM_E_TMINVAL);
273 +        if (len < 0)
274 +                returnErr(TM_E_ILLEGAL);
275 +        if (len == 0)
276 +                returnOK;
277                                                  /* first, grow limits */
278          if (tmTop->histo == NULL) {
279                  for (i = len; i-- && ls[i] < MINBRT; )
# Line 325 | Line 338 | double La;
338   }
339  
340  
341 + static int
342 + tmNewMap()
343 + {
344 +        if (tmTop->lumap != NULL && (tmTop->mbrmax - tmTop->mbrmin) !=
345 +                                        (tmTop->hbrmax - tmTop->hbrmin)) {
346 +                free((MEM_PTR)tmTop->lumap);
347 +                tmTop->lumap = NULL;
348 +        }
349 +        tmTop->mbrmin = tmTop->hbrmin;
350 +        tmTop->mbrmax = tmTop->hbrmax;
351 +        if (tmTop->mbrmin > tmTop->mbrmax)
352 +                return 0;
353 +        if (tmTop->lumap == NULL)
354 +                tmTop->lumap = (unsigned short *)malloc(sizeof(unsigned short)*
355 +                                        (tmTop->mbrmax-tmTop->mbrmin+1));
356 +        return(tmTop->lumap != NULL);
357 + }
358 +
359 +
360   int
361 + tmFixedMapping(expmult, gamval)
362 + double  expmult;
363 + double  gamval;
364 + {
365 +        static char     funcName[] = "tmFixedMapping";
366 +        double          d;
367 +        register int    i;
368 +        
369 +        if (!tmNewMap())
370 +                returnErr(TM_E_NOMEM);
371 +        if (expmult <= .0)
372 +                expmult = 1.;
373 +        if (gamval < MINGAM)
374 +                gamval = tmTop->mongam;
375 +        d = log(expmult/tmTop->inpsf);
376 +        for (i = tmTop->mbrmax-tmTop->mbrmin+1; i--; )
377 +                tmTop->lumap[i] = 256. * exp(
378 +                        ( d + (tmTop->mbrmin+i)*(1./TM_BRTSCALE) )
379 +                        / gamval );
380 +        returnOK;
381 + }
382 +
383 +
384 + int
385   tmComputeMapping(gamval, Lddyn, Ldmax)
386   double  gamval;
387   double  Lddyn;
# Line 334 | Line 390 | double Ldmax;
390          static char     funcName[] = "tmComputeMapping";
391          int     *histo;
392          float   *cumf;
393 <        int     brt0, histlen, histot, threshold, ceiling, trimmings;
393 >        int     brt0, histlen, threshold, ceiling, trimmings;
394          double  logLddyn, Ldmin, Ldavg, Lwavg, Tr, Lw, Ld;
395 <        int4    sum;
395 >        int32   histot;
396 >        double  sum;
397          register double d;
398          register int    i, j;
399  
# Line 360 | Line 417 | double Ldmax;
417                  histot += tmTop->histo[i];
418                  sum += (j -= HISTEP) * tmTop->histo[i];
419          }
420 <        threshold = histot*.025 + .5;
420 >        threshold = histot*0.005 + .5;
421          if (threshold < 4)
422                  returnErr(TM_E_TMFAIL);
423          Lwavg = tmLuminance( (double)sum / histot );
424 <        if (!(tmTop->flags & TM_F_LINEAR)) {    /* clamp histogram */
425 <                histo = (int *)malloc(histlen*sizeof(int));
426 <                cumf = (float *)malloc((histlen+1)*sizeof(float));
427 <                if (histo == NULL | cumf == NULL)
428 <                        returnErr(TM_E_NOMEM);
429 <                for (i = histlen; i--; )        /* make malleable copy */
430 <                        histo[i] = tmTop->histo[i];
431 <                do {                            /* iterate to solution */
432 <                        sum = 0;                /* cumulative probability */
433 <                        for (i = 0; i < histlen; i++) {
434 <                                cumf[i] = (double)sum/histot;
435 <                                sum += histo[i];
424 >                                        /* allocate space for mapping */
425 >        if (!tmNewMap())
426 >                returnErr(TM_E_NOMEM);
427 >                                        /* use linear tone mapping? */
428 >        if (tmTop->flags & TM_F_LINEAR)
429 >                goto linearmap;
430 >                                        /* clamp histogram */
431 >        histo = (int *)malloc(histlen*sizeof(int));
432 >        cumf = (float *)malloc((histlen+2)*sizeof(float));
433 >        if ((histo == NULL) | (cumf == NULL))
434 >                returnErr(TM_E_NOMEM);
435 >        cumf[histlen+1] = 1.;           /* guard for assignment code */
436 >        for (i = histlen; i--; )        /* make malleable copy */
437 >                histo[i] = tmTop->histo[i];
438 >        do {                            /* iterate to solution */
439 >                sum = 0;                /* cumulative probability */
440 >                for (i = 0; i < histlen; i++) {
441 >                        cumf[i] = (double)sum/histot;
442 >                        sum += histo[i];
443 >                }
444 >                cumf[histlen] = 1.;
445 >                Tr = histot * (double)(tmTop->hbrmax - tmTop->hbrmin) /
446 >                        ((double)histlen*TM_BRTSCALE) / logLddyn;
447 >                ceiling = Tr + 1.;
448 >                trimmings = 0;          /* clip to envelope */
449 >                for (i = histlen; i--; ) {
450 >                        if (tmTop->flags & TM_F_HCONTR) {
451 >                                Lw = tmLuminance(brt0 + i*HISTEP);
452 >                                Ld = Ldmin * exp( logLddyn *
453 >                                        .5*(cumf[i]+cumf[i+1]) );
454 >                                ceiling = Tr * (htcontrs(Ld) * Lw) /
455 >                                        (htcontrs(Lw) * Ld) + 1.;
456                          }
457 <                        cumf[i] = 1.;
458 <                        Tr = histot * (double)(tmTop->hbrmax - tmTop->hbrmin) /
459 <                                ((double)histlen*TM_BRTSCALE) / logLddyn;
383 <                        ceiling = Tr + 1.;
384 <                        trimmings = 0;                  /* clip to envelope */
385 <                        for (i = histlen; i--; ) {
386 <                                if (tmTop->flags & TM_F_HCONTR) {
387 <                                        Lw = tmLuminance(brt0 + i*HISTEP);
388 <                                        Ld = Ldmin * exp( logLddyn *
389 <                                                .5*(cumf[i]+cumf[i+1]) );
390 <                                        ceiling = Tr * (htcontrs(Ld) * Lw) /
391 <                                                (htcontrs(Lw) * Ld) + 1.;
392 <                                }
393 <                                if (histo[i] > ceiling) {
394 <                                        trimmings += histo[i] - ceiling;
395 <                                        histo[i] = ceiling;
396 <                                }
457 >                        if (histo[i] > ceiling) {
458 >                                trimmings += histo[i] - ceiling;
459 >                                histo[i] = ceiling;
460                          }
398 –                } while ((histot -= trimmings) > threshold &&
399 –                                                trimmings > threshold);
400 –        }
401 –                                                /* allocate luminance map */
402 –        if (tmTop->lumap != NULL)
403 –                free((MEM_PTR)tmTop->lumap);
404 –        tmTop->mbrmin = tmTop->hbrmin;
405 –        tmTop->mbrmax = tmTop->hbrmax;
406 –        tmTop->lumap = (unsigned short *)malloc(
407 –                (tmTop->mbrmax-tmTop->mbrmin+1)*sizeof(unsigned short) );
408 –        if (tmTop->lumap == NULL)
409 –                returnErr(TM_E_NOMEM);
410 –        if (tmTop->flags & TM_F_LINEAR || histot <= threshold) {
411 –                                                /* linear tone mapping */
412 –                if (tmTop->flags & TM_F_HCONTR)
413 –                        d = htcontrs(Ldavg) / htcontrs(Lwavg);
414 –                else
415 –                        d = Ldavg / Lwavg;
416 –                d = log(d/Ldmax);
417 –                for (i = tmTop->mbrmax-tmTop->mbrmin+1; i--; )
418 –                        tmTop->lumap[i] = 256. * exp(
419 –                                ( d + (tmTop->mbrmin+i)/(double)TM_BRTSCALE )
420 –                                / gamval );
421 –        } else {
422 –                                                /* histogram adjustment */
423 –                for (i = tmTop->mbrmax-tmTop->mbrmin+1; i--; ) {
424 –                        j = d = (double)i/(tmTop->mbrmax-tmTop->mbrmin)*histlen;
425 –                        d -= (double)j;
426 –                        Ld = Ldmin*exp(logLddyn*((1.-d)*cumf[j]+d*cumf[j+1]));
427 –                        d = (Ld - Ldmin)/(Ldmax - Ldmin);
428 –                        tmTop->lumap[i] = 256.*pow(d, 1./gamval);
461                  }
462 +                                        /* check if we're out of data */
463 +                if ((histot -= trimmings) <= threshold) {
464 +                        free((MEM_PTR)histo);
465 +                        free((MEM_PTR)cumf);
466 +                        goto linearmap;
467 +                }
468 +        } while (trimmings > threshold);
469 +                                        /* assign tone-mapping */
470 +        for (i = tmTop->mbrmax-tmTop->mbrmin+1; i--; ) {
471 +                j = d = (double)i/(tmTop->mbrmax-tmTop->mbrmin)*histlen;
472 +                d -= (double)j;
473 +                Ld = Ldmin*exp(logLddyn*((1.-d)*cumf[j]+d*cumf[j+1]));
474 +                d = (Ld - Ldmin)/(Ldmax - Ldmin);
475 +                tmTop->lumap[i] = 256.*pow(d, 1./gamval);
476          }
477 <        if (!(tmTop->flags & TM_F_LINEAR)) {
478 <                free((MEM_PTR)histo);
433 <                free((MEM_PTR)cumf);
434 <        }
477 >        free((MEM_PTR)histo);           /* clean up and return */
478 >        free((MEM_PTR)cumf);
479          returnOK;
480 + linearmap:                              /* linear tone-mapping */
481 +        if (tmTop->flags & TM_F_HCONTR)
482 +                d = htcontrs(Ldavg) / htcontrs(Lwavg);
483 +        else
484 +                d = Ldavg / Lwavg;
485 +        return(tmFixedMapping(tmTop->inpsf*d/Ldmax, gamval));
486   }
487  
488  
# Line 444 | Line 494 | register BYTE  *cs;
494   int     len;
495   {
496          static char     funcName[] = "tmMapPixels";
497 <        register int4   li, pv;
497 >        register int32  li, pv;
498  
499          if (tmTop == NULL || tmTop->lumap == NULL)
500                  returnErr(TM_E_TMINVAL);
501 <        if (ps == NULL | ls == NULL | len <= 0)
501 >        if ((ps == NULL) | (ls == NULL) | (len < 0))
502                  returnErr(TM_E_ILLEGAL);
503          while (len--) {
504 <                if ((li = *ls++) < tmTop->mbrmin)
505 <                        li = tmTop->mbrmin;
506 <                else if (li > tmTop->mbrmax)
507 <                        li = tmTop->mbrmax;
508 <                li = tmTop->lumap[li - tmTop->mbrmin];
504 >                if ((li = *ls++) < tmTop->mbrmin) {
505 >                        li = 0;
506 >                } else {
507 >                        if (li > tmTop->mbrmax)
508 >                                li = tmTop->mbrmax;
509 >                        li = tmTop->lumap[li - tmTop->mbrmin];
510 >                }
511                  if (cs == TM_NOCHROM)
512                          *ps++ = li>255 ? 255 : li;
513                  else {
# Line 488 | Line 540 | register struct tmStruct       *tms;
540   {
541          register struct tmStruct        *tms2;
542                                          /* special cases first */
543 <        if (tms == NULL | tmTop == NULL)
543 >        if ((tms == NULL) | (tmTop == NULL))
544                  return(0);
545          if (tms == tmTop) {
546                  tmTop = tms->tmprev;
# Line 581 | Line 633 | register struct tmStruct       *tms;
633                  if (tms->pd[i] != NULL)
634                          (*tmPkg[i]->Free)(tms->pd[i]);
635          free((MEM_PTR)tms);             /* free basic structure */
636 + }
637 +
638 + /******************** Shared but Private library routines *********************/
639 +
640 + BYTE    tmMesofact[BMESUPPER-BMESLOWER];
641 +
642 + void
643 + tmMkMesofact()                          /* build mesopic lookup factor table */
644 + {
645 +        register int    i;
646 +
647 +        if (tmMesofact[BMESUPPER-BMESLOWER-1])
648 +                return;
649 +
650 +        for (i = BMESLOWER; i < BMESUPPER; i++)
651 +                tmMesofact[i-BMESLOWER] = 256. *
652 +                                (tmLuminance(i) - LMESLOWER) /
653 +                                (LMESUPPER - LMESLOWER);
654 + }
655 +
656 +
657 + int
658 + tmErrorReturn(func, err)                /* error return (with message) */
659 + char    *func;
660 + int     err;
661 + {
662 +        tmLastFunction = func;
663 +        tmLastError = err;
664 +        if (tmTop != NULL && tmTop->flags & TM_F_NOSTDERR)
665 +                return(err);
666 +        fputs(func, stderr);
667 +        fputs(": ", stderr);
668 +        fputs(tmErrorMessage[err], stderr);
669 +        fputs("!\n", stderr);
670 +        return(err);
671   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines