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.4 by greg, Fri Apr 18 13:59:43 1997 UTC vs.
Revision 3.14 by greg, Thu Sep 18 20:18:19 2003 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 22 | Line 24 | struct tmStruct        *tmTop = NULL;          /* current tone mappin
24   struct tmPackage        *tmPkg[TM_MAXPKG];
25   int     tmNumPkgs = 0;                  /* number of registered packages */
26  
27 + int     tmLastError;                    /* last error incurred by library */
28 + char    *tmLastFunction;                /* error-generating function name */
29  
26 – int
27 – tmErrorReturn(func, err)                /* error return (with message) */
28 – char    *func;
29 – int     err;
30 – {
31 –        if (tmTop != NULL && tmTop->flags & TM_F_NOSTDERR)
32 –                return(err);
33 –        fputs(func, stderr);
34 –        fputs(": ", stderr);
35 –        fputs(tmErrorMessage[err], stderr);
36 –        fputs("!\n", stderr);
37 –        return(err);
38 – }
30  
40 –
31   struct tmStruct *
32   tmInit(flags, monpri, gamval)           /* initialize new tone mapping */
33   int     flags;
34   RGBPRIMP        monpri;
35   double  gamval;
36   {
47 –        static char     funcName[] = "tmInit";
37          COLORMAT        cmat;
38          register struct tmStruct        *tmnew;
39          register int    i;
# Line 61 | 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 81 | 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->brmin = tmnew->brmax = 0;
73 >        tmnew->hbrmin = 10; tmnew->hbrmax = -10;
74          tmnew->histo = NULL;
75 +        tmnew->mbrmin = 10; tmnew->mbrmax = -10;
76          tmnew->lumap = NULL;
77                                                  /* zero private data */
78          for (i = TM_MAXPKG; i--; )
# Line 117 | 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 128 | 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 172 | 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 180 | Line 170 | int    len;
170  
171          if (tmTop == NULL)
172                  returnErr(TM_E_TMINVAL);
173 <        if (ls == NULL | scan == NULL | len <= 0)
173 >        if ((ls == NULL) | (scan == NULL) | (len < 0))
174                  returnErr(TM_E_ILLEGAL);
175          for (i = len; i--; ) {
176                  if (tmNeedMatrix(tmTop))                /* get monitor RGB */
# Line 242 | 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     sum, oldorig, oldlen, horig, hlen;
266 >        int     oldorig=0, oldlen, horig, hlen;
267          register int    i, j;
268  
254 –        if (len <= 0)
255 –                returnErr(TM_E_ILLEGAL);
269          if (tmTop == NULL)
270                  returnErr(TM_E_TMINVAL);
271 +        if (len < 0)
272 +                returnErr(TM_E_ILLEGAL);
273 +        if (len == 0)
274 +                returnOK;
275                                                  /* first, grow limits */
276          if (tmTop->histo == NULL) {
277                  for (i = len; i-- && ls[i] < MINBRT; )
278                          ;
279                  if (i < 0)
280                          returnOK;
281 <                tmTop->brmin = tmTop->brmax = ls[i];
281 >                tmTop->hbrmin = tmTop->hbrmax = ls[i];
282                  oldlen = 0;
283          } else {
284 <                oldorig = (tmTop->brmin-MINBRT)/HISTEP;
285 <                oldlen = (tmTop->brmax-MINBRT)/HISTEP + 1 - oldorig;
284 >                oldorig = (tmTop->hbrmin-MINBRT)/HISTEP;
285 >                oldlen = (tmTop->hbrmax-MINBRT)/HISTEP + 1 - oldorig;
286          }
287          for (i = len; i--; ) {
288                  if ((j = ls[i]) < MINBRT)
289                          continue;
290 <                if (j < tmTop->brmin)
291 <                        tmTop->brmin = j;
292 <                else if (j > tmTop->brmax)
293 <                        tmTop->brmax = j;
290 >                if (j < tmTop->hbrmin)
291 >                        tmTop->hbrmin = j;
292 >                else if (j > tmTop->hbrmax)
293 >                        tmTop->hbrmax = j;
294          }
295 <        horig = (tmTop->brmin-MINBRT)/HISTEP;
296 <        hlen = (tmTop->brmax-MINBRT)/HISTEP + 1 - horig;
295 >        horig = (tmTop->hbrmin-MINBRT)/HISTEP;
296 >        hlen = (tmTop->hbrmax-MINBRT)/HISTEP + 1 - horig;
297          if (hlen > oldlen) {                    /* (re)allocate histogram */
298                  register int    *newhist = (int *)calloc(hlen, sizeof(int));
299                  if (newhist == NULL)
# Line 287 | Line 304 | int    wt;
304                          free((MEM_PTR)tmTop->histo);
305                  }
306                  tmTop->histo = newhist;
290 –                if (tmTop->lumap != NULL) {     /* invalid tone map */
291 –                        free((MEM_PTR)tmTop->lumap);
292 –                        tmTop->lumap = NULL;
293 –                }
307          }
308          if (wt == 0)
309                  returnOK;
# Line 323 | 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 332 | 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 348 | Line 405 | double Ldmax;
405          Ldmin = Ldmax/Lddyn;
406          logLddyn = log(Lddyn);
407          Ldavg = sqrt(Ldmax*Ldmin);
408 <        i = (tmTop->brmin-MINBRT)/HISTEP;
408 >        i = (tmTop->hbrmin-MINBRT)/HISTEP;
409          brt0 = MINBRT + HISTEP/2 + i*HISTEP;
410 <        histlen = (tmTop->brmax-MINBRT)/HISTEP + 1 - i;
410 >        histlen = (tmTop->hbrmax-MINBRT)/HISTEP + 1 - i;
411                                          /* histogram total and mean */
412          histot = 0; sum = 0;
413          j = brt0 + histlen*HISTEP;
# Line 358 | Line 415 | double Ldmax;
415                  histot += tmTop->histo[i];
416                  sum += (j -= HISTEP) * tmTop->histo[i];
417          }
418 <        threshold = histot*.025 + .5;
418 >        threshold = histot*0.005 + .5;
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->brmax - tmTop->brmin) /
457 <                                ((double)histlen*TM_BRTSCALE) / logLddyn;
381 <                        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 <                                                .5*(cumf[i]+cumf[i+1]) );
388 <                                        ceiling = Tr * (htcontrs(Ld) * Lw) /
389 <                                                (htcontrs(Lw) * Ld) + 1.;
390 <                                }
391 <                                if (histo[i] > ceiling) {
392 <                                        trimmings += histo[i] - ceiling;
393 <                                        histo[i] = ceiling;
394 <                                }
455 >                        if (histo[i] > ceiling) {
456 >                                trimmings += histo[i] - ceiling;
457 >                                histo[i] = ceiling;
458                          }
396 –                } while ((histot -= trimmings) > threshold &&
397 –                                                trimmings > threshold);
398 –        }
399 –                                                /* allocate luminance map */
400 –        if (tmTop->lumap == NULL) {
401 –                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);
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);
429 <                free((MEM_PTR)cumf);
430 <        }
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 440 | 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)
499 >        if ((ps == NULL) | (ls == NULL) | (len < 0))
500                  returnErr(TM_E_ILLEGAL);
501          while (len--) {
502 <                if ((li = *ls++) < tmTop->brmin)
503 <                        li = tmTop->brmin;
504 <                else if (li > tmTop->brmax)
505 <                        li = tmTop->brmax;
506 <                li = tmTop->lumap[li - tmTop->brmin];
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 484 | Line 538 | register struct tmStruct       *tms;
538   {
539          register struct tmStruct        *tms2;
540                                          /* special cases first */
541 <        if (tms == NULL | tmTop == NULL)
541 >        if ((tms == NULL) | (tmTop == NULL))
542                  return(0);
543          if (tms == tmTop) {
544                  tmTop = tms->tmprev;
# Line 515 | Line 569 | tmDup()                                /* duplicate top tone mapping */
569                  return(NULL);
570          *tmnew = *tmTop;                /* copy everything */
571          if (tmnew->histo != NULL) {     /* duplicate histogram */
572 <                len = (tmnew->brmax-MINBRT)/HISTEP + 1 -
573 <                                (tmnew->brmin-MINBRT)/HISTEP;
572 >                len = (tmnew->hbrmax-MINBRT)/HISTEP + 1 -
573 >                                (tmnew->hbrmin-MINBRT)/HISTEP;
574                  tmnew->histo = (int *)malloc(len*sizeof(int));
575                  if (tmnew->histo != NULL)
576                          for (i = len; i--; )
577                                  tmnew->histo[i] = tmTop->histo[i];
578          }
579          if (tmnew->lumap != NULL) {     /* duplicate luminance mapping */
580 <                len = tmnew->brmax-tmnew->brmin+1;
580 >                len = tmnew->mbrmax-tmnew->mbrmin+1;
581                  tmnew->lumap = (unsigned short *)malloc(
582                                                  len*sizeof(unsigned short) );
583                  if (tmnew->lumap != NULL)
# Line 577 | 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