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.9 by greg, Sat Feb 22 02:07:22 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 + /* ====================================================================
13 + * The Radiance Software License, Version 1.0
14 + *
15 + * Copyright (c) 1990 - 2002 The Regents of the University of California,
16 + * through Lawrence Berkeley National Laboratory.   All rights reserved.
17 + *
18 + * Redistribution and use in source and binary forms, with or without
19 + * modification, are permitted provided that the following conditions
20 + * are met:
21 + *
22 + * 1. Redistributions of source code must retain the above copyright
23 + *         notice, this list of conditions and the following disclaimer.
24 + *
25 + * 2. Redistributions in binary form must reproduce the above copyright
26 + *       notice, this list of conditions and the following disclaimer in
27 + *       the documentation and/or other materials provided with the
28 + *       distribution.
29 + *
30 + * 3. The end-user documentation included with the redistribution,
31 + *           if any, must include the following acknowledgment:
32 + *             "This product includes Radiance software
33 + *                 (http://radsite.lbl.gov/)
34 + *                 developed by the Lawrence Berkeley National Laboratory
35 + *               (http://www.lbl.gov/)."
36 + *       Alternately, this acknowledgment may appear in the software itself,
37 + *       if and wherever such third-party acknowledgments normally appear.
38 + *
39 + * 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
40 + *       and "The Regents of the University of California" must
41 + *       not be used to endorse or promote products derived from this
42 + *       software without prior written permission. For written
43 + *       permission, please contact [email protected].
44 + *
45 + * 5. Products derived from this software may not be called "Radiance",
46 + *       nor may "Radiance" appear in their name, without prior written
47 + *       permission of Lawrence Berkeley National Laboratory.
48 + *
49 + * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
50 + * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
51 + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
52 + * DISCLAIMED.   IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
53 + * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
54 + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
55 + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
56 + * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
57 + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
58 + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
59 + * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
60 + * SUCH DAMAGE.
61 + * ====================================================================
62 + *
63 + * This software consists of voluntary contributions made by many
64 + * individuals on behalf of Lawrence Berkeley National Laboratory.   For more
65 + * information on Lawrence Berkeley National Laboratory, please see
66 + * <http://www.lbl.gov/>.
67 + */
68 +
69   #include        <stdio.h>
70   #include        <math.h>
71   #include        "tmprivat.h"
# Line 26 | Line 83 | int    tmLastError;                    /* last error incurred by library *
83   char    *tmLastFunction;                /* error-generating function name */
84  
85  
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 –
86   struct tmStruct *
87   tmInit(flags, monpri, gamval)           /* initialize new tone mapping */
88   int     flags;
89   RGBPRIMP        monpri;
90   double  gamval;
91   {
52 –        static char     funcName[] = "tmInit";
92          COLORMAT        cmat;
93          register struct tmStruct        *tmnew;
94          register int    i;
# Line 66 | Line 105 | double gamval;
105                  tmnew->clf[GRN] = rgb2xyzmat[1][1];
106                  tmnew->clf[BLU] = rgb2xyzmat[1][2];
107          } else {
108 <                comprgb2xyzmat(cmat, tmnew->monpri=monpri);
108 >                comprgb2xyzWBmat(cmat, tmnew->monpri=monpri);
109                  tmnew->clf[RED] = cmat[1][0];
110                  tmnew->clf[GRN] = cmat[1][1];
111                  tmnew->clf[BLU] = cmat[1][2];
# Line 86 | Line 125 | double gamval;
125                          tmnew->inpsf = WHTEFFICACY;
126          tmnew->cmat[0][1] = tmnew->cmat[0][2] = tmnew->cmat[1][0] =
127          tmnew->cmat[1][2] = tmnew->cmat[2][0] = tmnew->cmat[2][1] = 0.;
128 <        tmnew->hbrmin = tmnew->hbrmax = 0;
128 >        tmnew->hbrmin = 10; tmnew->hbrmax = -10;
129          tmnew->histo = NULL;
130 <        tmnew->mbrmin = tmnew->mbrmax = 0;
130 >        tmnew->mbrmin = 10; tmnew->mbrmax = -10;
131          tmnew->lumap = NULL;
132                                                  /* zero private data */
133          for (i = TM_MAXPKG; i--; )
# Line 123 | Line 162 | double sf;
162                          tmTop->clf[CIEX] = tmTop->clf[CIEZ] = 0.;
163                          tmTop->clf[CIEY] = 1.;
164                  } else {
165 <                        comprgb2xyzmat(tmTop->cmat, tmTop->monpri);
165 >                        comprgb2xyzWBmat(tmTop->cmat, tmTop->monpri);
166                          tmTop->clf[RED] = tmTop->cmat[1][0];
167                          tmTop->clf[GRN] = tmTop->cmat[1][1];
168                          tmTop->clf[BLU] = tmTop->cmat[1][2];
# Line 134 | Line 173 | double sf;
173                  tmTop->cmat[1][2] = tmTop->cmat[2][0] = tmTop->cmat[2][1] = 0.;
174  
175          } else if (tmTop->inppri == TM_XYZPRIM) /* input is XYZ */
176 <                compxyz2rgbmat(tmTop->cmat, tmTop->monpri);
176 >                compxyz2rgbWBmat(tmTop->cmat, tmTop->monpri);
177  
178          else {                                  /* input is RGB */
179                  if (tmTop->inppri != tmTop->monpri &&
180                                  PRIMEQ(tmTop->inppri, tmTop->monpri))
181                          tmTop->inppri = tmTop->monpri;  /* no xform */
182 <                comprgb2rgbmat(tmTop->cmat, tmTop->inppri, tmTop->monpri);
182 >                comprgb2rgbWBmat(tmTop->cmat, tmTop->inppri, tmTop->monpri);
183          }
184          for (i = 0; i < 3; i++)
185                  for (j = 0; j < 3; j++)
# Line 248 | Line 287 | int    len;
287  
288  
289   int
290 + tmCvGrays(ls, scan, len)                /* convert float gray values */
291 + TMbright        *ls;
292 + float   *scan;
293 + int     len;
294 + {
295 +        static char     funcName[] = "tmCvGrays";
296 +        register double d;
297 +        register int    i;
298 +
299 +        if (tmTop == NULL)
300 +                returnErr(TM_E_TMINVAL);
301 +        if (ls == NULL | scan == NULL | len < 0)
302 +                returnErr(TM_E_ILLEGAL);
303 +        for (i = len; i--; )
304 +                if (scan[i] <= TM_NOLUM) {
305 +                        ls[i] = TM_NOBRT;               /* bogus value */
306 +                } else {
307 +                        d = TM_BRTSCALE*log(scan[i]);   /* encode it */
308 +                        ls[i] = d>0. ? (int)(d+.5) : (int)(d-.5);
309 +                }
310 +        returnOK;
311 + }
312 +
313 +
314 + int
315   tmAddHisto(ls, len, wt)                 /* add values to histogram */
316   register TMbright       *ls;
317   int     len;
318   int     wt;
319   {
320          static char     funcName[] = "tmAddHisto";
321 <        int     oldorig, oldlen, horig, hlen;
321 >        int     oldorig=0, oldlen, horig, hlen;
322          register int    i, j;
323  
324          if (tmTop == NULL)
# Line 327 | Line 391 | double La;
391   }
392  
393  
394 + static int
395 + tmNewMap()
396 + {
397 +        if (tmTop->lumap != NULL && (tmTop->mbrmax - tmTop->mbrmin) !=
398 +                                        (tmTop->hbrmax - tmTop->hbrmin)) {
399 +                free((MEM_PTR)tmTop->lumap);
400 +                tmTop->lumap = NULL;
401 +        }
402 +        tmTop->mbrmin = tmTop->hbrmin;
403 +        tmTop->mbrmax = tmTop->hbrmax;
404 +        if (tmTop->mbrmin > tmTop->mbrmax)
405 +                return 0;
406 +        if (tmTop->lumap == NULL)
407 +                tmTop->lumap = (unsigned short *)malloc(sizeof(unsigned short)*
408 +                                        (tmTop->mbrmax-tmTop->mbrmin+1));
409 +        return(tmTop->lumap != NULL);
410 + }
411 +
412 +
413   int
414 + tmFixedMapping(expmult, gamval)
415 + double  expmult;
416 + double  gamval;
417 + {
418 +        static char     funcName[] = "tmFixedMapping";
419 +        double          d;
420 +        register int    i;
421 +        
422 +        if (!tmNewMap())
423 +                returnErr(TM_E_NOMEM);
424 +        if (expmult <= .0)
425 +                expmult = 1.;
426 +        if (gamval < MINGAM)
427 +                gamval = tmTop->mongam;
428 +        d = log(expmult/tmTop->inpsf);
429 +        for (i = tmTop->mbrmax-tmTop->mbrmin+1; i--; )
430 +                tmTop->lumap[i] = 256. * exp(
431 +                        ( d + (tmTop->mbrmin+i)*(1./TM_BRTSCALE) )
432 +                        / gamval );
433 +        returnOK;
434 + }
435 +
436 +
437 + int
438   tmComputeMapping(gamval, Lddyn, Ldmax)
439   double  gamval;
440   double  Lddyn;
# Line 336 | Line 443 | double Ldmax;
443          static char     funcName[] = "tmComputeMapping";
444          int     *histo;
445          float   *cumf;
446 <        int     brt0, histlen, histot, threshold, ceiling, trimmings;
446 >        int     brt0, histlen, threshold, ceiling, trimmings;
447          double  logLddyn, Ldmin, Ldavg, Lwavg, Tr, Lw, Ld;
448 <        int4    sum;
448 >        int4    histot;
449 >        double  sum;
450          register double d;
451          register int    i, j;
452  
# Line 366 | Line 474 | double Ldmax;
474          if (threshold < 4)
475                  returnErr(TM_E_TMFAIL);
476          Lwavg = tmLuminance( (double)sum / histot );
477 <        if (!(tmTop->flags & TM_F_LINEAR)) {    /* clamp histogram */
478 <                histo = (int *)malloc(histlen*sizeof(int));
479 <                cumf = (float *)malloc((histlen+1)*sizeof(float));
480 <                if (histo == NULL | cumf == NULL)
481 <                        returnErr(TM_E_NOMEM);
482 <                for (i = histlen; i--; )        /* make malleable copy */
483 <                        histo[i] = tmTop->histo[i];
484 <                do {                            /* iterate to solution */
485 <                        sum = 0;                /* cumulative probability */
486 <                        for (i = 0; i < histlen; i++) {
487 <                                cumf[i] = (double)sum/histot;
488 <                                sum += histo[i];
477 >                                        /* allocate space for mapping */
478 >        if (!tmNewMap())
479 >                returnErr(TM_E_NOMEM);
480 >                                        /* use linear tone mapping? */
481 >        if (tmTop->flags & TM_F_LINEAR)
482 >                goto linearmap;
483 >                                        /* clamp histogram */
484 >        histo = (int *)malloc(histlen*sizeof(int));
485 >        cumf = (float *)malloc((histlen+2)*sizeof(float));
486 >        if (histo == NULL | cumf == NULL)
487 >                returnErr(TM_E_NOMEM);
488 >        cumf[histlen+1] = 1.;           /* guard for assignment code */
489 >        for (i = histlen; i--; )        /* make malleable copy */
490 >                histo[i] = tmTop->histo[i];
491 >        do {                            /* iterate to solution */
492 >                sum = 0;                /* cumulative probability */
493 >                for (i = 0; i < histlen; i++) {
494 >                        cumf[i] = (double)sum/histot;
495 >                        sum += histo[i];
496 >                }
497 >                cumf[histlen] = 1.;
498 >                Tr = histot * (double)(tmTop->hbrmax - tmTop->hbrmin) /
499 >                        ((double)histlen*TM_BRTSCALE) / logLddyn;
500 >                ceiling = Tr + 1.;
501 >                trimmings = 0;          /* clip to envelope */
502 >                for (i = histlen; i--; ) {
503 >                        if (tmTop->flags & TM_F_HCONTR) {
504 >                                Lw = tmLuminance(brt0 + i*HISTEP);
505 >                                Ld = Ldmin * exp( logLddyn *
506 >                                        .5*(cumf[i]+cumf[i+1]) );
507 >                                ceiling = Tr * (htcontrs(Ld) * Lw) /
508 >                                        (htcontrs(Lw) * Ld) + 1.;
509                          }
510 <                        cumf[i] = 1.;
511 <                        Tr = histot * (double)(tmTop->hbrmax - tmTop->hbrmin) /
512 <                                ((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 <                                }
510 >                        if (histo[i] > ceiling) {
511 >                                trimmings += histo[i] - ceiling;
512 >                                histo[i] = ceiling;
513                          }
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);
514                  }
515 +                                        /* check if we're out of data */
516 +                if ((histot -= trimmings) <= threshold) {
517 +                        free((MEM_PTR)histo);
518 +                        free((MEM_PTR)cumf);
519 +                        goto linearmap;
520 +                }
521 +        } while (trimmings > threshold);
522 +                                        /* assign tone-mapping */
523 +        for (i = tmTop->mbrmax-tmTop->mbrmin+1; i--; ) {
524 +                j = d = (double)i/(tmTop->mbrmax-tmTop->mbrmin)*histlen;
525 +                d -= (double)j;
526 +                Ld = Ldmin*exp(logLddyn*((1.-d)*cumf[j]+d*cumf[j+1]));
527 +                d = (Ld - Ldmin)/(Ldmax - Ldmin);
528 +                tmTop->lumap[i] = 256.*pow(d, 1./gamval);
529          }
530 <        if (!(tmTop->flags & TM_F_LINEAR)) {
531 <                free((MEM_PTR)histo);
435 <                free((MEM_PTR)cumf);
436 <        }
530 >        free((MEM_PTR)histo);           /* clean up and return */
531 >        free((MEM_PTR)cumf);
532          returnOK;
533 + linearmap:                              /* linear tone-mapping */
534 +        if (tmTop->flags & TM_F_HCONTR)
535 +                d = htcontrs(Ldavg) / htcontrs(Lwavg);
536 +        else
537 +                d = Ldavg / Lwavg;
538 +        return(tmFixedMapping(tmTop->inpsf*d/Ldmax, gamval));
539   }
540  
541  
# Line 453 | Line 554 | int    len;
554          if (ps == NULL | ls == NULL | len < 0)
555                  returnErr(TM_E_ILLEGAL);
556          while (len--) {
557 <                if ((li = *ls++) < tmTop->mbrmin)
558 <                        li = tmTop->mbrmin;
559 <                else if (li > tmTop->mbrmax)
560 <                        li = tmTop->mbrmax;
561 <                li = tmTop->lumap[li - tmTop->mbrmin];
557 >                if ((li = *ls++) < tmTop->mbrmin) {
558 >                        li = 0;
559 >                } else {
560 >                        if (li > tmTop->mbrmax)
561 >                                li = tmTop->mbrmax;
562 >                        li = tmTop->lumap[li - tmTop->mbrmin];
563 >                }
564                  if (cs == TM_NOCHROM)
565                          *ps++ = li>255 ? 255 : li;
566                  else {
# Line 583 | Line 686 | register struct tmStruct       *tms;
686                  if (tms->pd[i] != NULL)
687                          (*tmPkg[i]->Free)(tms->pd[i]);
688          free((MEM_PTR)tms);             /* free basic structure */
689 + }
690 +
691 + /******************** Shared but Private library routines *********************/
692 +
693 + BYTE    tmMesofact[BMESUPPER-BMESLOWER];
694 +
695 + void
696 + tmMkMesofact()                          /* build mesopic lookup factor table */
697 + {
698 +        register int    i;
699 +
700 +        if (tmMesofact[BMESUPPER-BMESLOWER-1])
701 +                return;
702 +
703 +        for (i = BMESLOWER; i < BMESUPPER; i++)
704 +                tmMesofact[i-BMESLOWER] = 256. *
705 +                                (tmLuminance(i) - LMESLOWER) /
706 +                                (LMESUPPER - LMESLOWER);
707 + }
708 +
709 +
710 + int
711 + tmErrorReturn(func, err)                /* error return (with message) */
712 + char    *func;
713 + int     err;
714 + {
715 +        tmLastFunction = func;
716 +        tmLastError = err;
717 +        if (tmTop != NULL && tmTop->flags & TM_F_NOSTDERR)
718 +                return(err);
719 +        fputs(func, stderr);
720 +        fputs(": ", stderr);
721 +        fputs(tmErrorMessage[err], stderr);
722 +        fputs("!\n", stderr);
723 +        return(err);
724   }

Diff Legend

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