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

Comparing ray/src/px/pcond3.c (file contents):
Revision 3.10 by gwlarson, Thu Mar 12 15:47:34 1998 UTC vs.
Revision 3.15 by schorsch, Sun Mar 28 20:33:14 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   * Routines for computing and applying brightness mapping.
6   */
7  
8 + #include <string.h>
9 +
10   #include "pcond.h"
11  
12  
# Line 16 | Line 15 | static char SCCSid[] = "$SunId$ LBL";
15   #define LN_10           2.30258509299404568402
16   #define exp10(x)        exp(LN_10*(x))
17  
18 < float   modhist[HISTRES];               /* modified histogram */
18 > double  modhist[HISTRES];               /* modified histogram */
19   double  mhistot;                        /* modified histogram total */
20 < float   cumf[HISTRES+1];                /* cumulative distribution function */
20 > double  cumf[HISTRES+1];                /* cumulative distribution function */
21  
22 + static double centprob(int      x, int  y);
23 + static void mkcumf(void);
24 + static double cf(double b);
25 + static double BLw(double        Lw);
26 + #if ADJ_VEIL
27 + static void mkcrfimage(void);
28 + #endif
29  
30 < getfixations(fp)                /* load fixation history list */
31 < FILE    *fp;
30 >
31 >
32 > extern void
33 > getfixations(           /* load fixation history list */
34 > FILE    *fp
35 > )
36   {
37   #define FIXHUNK         128
38          RESOLU  fvres;
39          int     pos[2];
40          register int    px, py, i;
41                                  /* initialize our resolution struct */
42 <        if ((fvres.or=inpres.or)&YMAJOR) {
42 >        if ((fvres.rt=inpres.rt)&YMAJOR) {
43                  fvres.xr = fvxr;
44                  fvres.yr = fvyr;
45          } else {
# Line 57 | Line 67 | FILE   *fp;
67                                  if (nfixations % FIXHUNK == 0) {
68                                          if (nfixations)
69                                                  fixlst = (short (*)[2])
70 <                                                        realloc((char *)fixlst,
70 >                                                        realloc((void *)fixlst,
71                                                          (nfixations+FIXHUNK)*
72                                                          2*sizeof(short));
73                                          else
# Line 82 | Line 92 | FILE   *fp;
92   }
93  
94  
95 < gethisto(fp)                    /* load precomputed luminance histogram */
96 < FILE    *fp;
95 > extern void
96 > gethisto(                       /* load precomputed luminance histogram */
97 >        FILE    *fp
98 > )
99   {
100 <        float   histo[MAXPREHIST];
100 >        double  histo[MAXPREHIST];
101          double  histart, histep;
102 <        double  l, b, lastb, w;
102 >        double  b, lastb, w;
103          int     n;
104          register int    i;
105                                          /* load data */
106          for (i = 0; i < MAXPREHIST &&
107 <                        fscanf(fp, "%lf %f", &b, &histo[i]) == 2; i++) {
108 <                if (i > 1 && fabs(b - lastb - histep) > 1e-3) {
107 >                        fscanf(fp, "%lf %lf", &b, &histo[i]) == 2; i++) {
108 >                if (i > 1 && fabs(b - lastb - histep) > .001) {
109                          fprintf(stderr,
110                                  "%s: uneven step size in histogram data\n",
111                                          progname);
112                          exit(1);
113                  }
114                  if (i == 1)
115 <                        histep = b - (histart = lastb);
115 >                        if ((histep = b - (histart = lastb)) <= FTINY) {
116 >                                fprintf(stderr,
117 >                                        "%s: illegal step in histogram data\n",
118 >                                                progname);
119 >                                exit(1);
120 >                        }
121                  lastb = b;
122          }
123          if (i < 2 || !feof(fp)) {
# Line 115 | Line 132 | FILE   *fp;
132                                          /* find extrema */
133          for (i = 0; i < n && histo[i] <= FTINY; i++)
134                  ;
135 <        bwmin = histart + i*histep;
135 >        bwmin = histart + (i-.001)*histep;
136          for (i = n; i-- && histo[i] <= FTINY; )
137                  ;
138 <        bwmax = histart + i*histep;
138 >        bwmax = histart + (i+1.001)*histep;
139          if (bwmax > Bl(LMAX))
140                  bwmax = Bl(LMAX);
141          if (bwmin < Bl(LMIN))
# Line 130 | Line 147 | FILE   *fp;
147          for (i = 0; i < HISTRES; i++)
148                  bwhist[i] = 0.;
149          for (i = 0, b = histart; i < n; i++, b += histep) {
150 <                if (b < bwmin) continue;
151 <                if (b > bwmax) break;
150 >                if (b < bwmin+FTINY) continue;
151 >                if (b >= bwmax-FTINY) break;
152                  w = histo[i];
153                  bwavg += w*b;
154                  bwhist[bwhi(b)] += w;
# Line 145 | Line 162 | FILE   *fp;
162   }
163  
164  
165 < double
166 < centprob(x, y)                  /* center-weighting probability function */
167 < int     x, y;
165 > static double
166 > centprob(                       /* center-weighting probability function */
167 >        int     x,
168 >        int     y
169 > )
170   {
171          double  xr, yr, p;
172                                  /* paraboloid, 0 at 90 degrees from center */
# Line 158 | Line 177 | int    x, y;
177   }
178  
179  
180 < comphist()                      /* create foveal sampling histogram */
180 > extern void
181 > comphist(void)                  /* create foveal sampling histogram */
182   {
183          double  l, b, w, lwmin, lwmax;
184          register int    x, y;
# Line 194 | Line 214 | comphist()                     /* create foveal sampling histogram */
214                  for (y = 0; y < fvyr; y++)
215                          for (x = 0; x < fvxr; x++) {
216                                  l = plum(fovscan(y)[x]);
217 <                                if (l < lwmin) continue;
218 <                                if (l > lwmax) continue;
217 >                                if (l < lwmin+FTINY) continue;
218 >                                if (l >= lwmax-FTINY) continue;
219                                  b = Bl(l);
220                                  w = what2do&DO_CWEIGHT ? centprob(x,y) : 1.;
221                                  bwavg += w*b;
# Line 210 | Line 230 | comphist()                     /* create foveal sampling histogram */
230                          w = 1.;
231                  for (x = 0; x < nfixations; x++) {
232                          l = plum(fovscan(fixlst[x][1])[fixlst[x][0]]);
233 <                        if (l < lwmin) continue;
234 <                        if (l > lwmax) continue;
233 >                        if (l < lwmin+FTINY) continue;
234 >                        if (l >= lwmax-FTINY) continue;
235                          b = Bl(l);
236                          bwavg += w*b;
237                          bwhist[bwhi(b)] += w;
# Line 226 | Line 246 | comphist()                     /* create foveal sampling histogram */
246   }
247  
248  
249 < mkcumf()                        /* make cumulative distribution function */
249 > static void
250 > mkcumf(void)                    /* make cumulative distribution function */
251   {
252          register int    i;
253          register double sum;
# Line 244 | Line 265 | mkcumf()                       /* make cumulative distribution function */
265   }
266  
267  
268 < double
269 < cf(b)                           /* return cumulative function at b */
270 < double  b;
268 > static double
269 > cf(                             /* return cumulative function at b */
270 >        double  b
271 > )
272   {
273          double  x;
274          register int    i;
# Line 257 | Line 279 | double b;
279   }
280  
281  
282 < double
283 < BLw(Lw)                         /* map world luminance to display brightness */
284 < double  Lw;
282 > static double
283 > BLw(                            /* map world luminance to display brightness */
284 >        double  Lw
285 > )
286   {
287          double  b;
288  
# Line 271 | Line 294 | double Lw;
294   }
295  
296  
297 < double
298 < htcontrs(La)            /* human threshold contrast sensitivity, dL(La) */
299 < double  La;
297 > extern double
298 > htcontrs(               /* human threshold contrast sensitivity, dL(La) */
299 >        double  La
300 > )
301   {
302          double  l10La, l10dL;
303                                  /* formula taken from Ferwerda et al. [SG96] */
# Line 293 | Line 317 | double La;
317   }
318  
319  
320 < double
321 < clampf(Lw)              /* derivative clamping function */
322 < double  Lw;
320 > extern double
321 > clampf(                 /* histogram clamping function */
322 >        double  Lw
323 > )
324   {
325          double  bLw, ratio;
326  
# Line 304 | Line 329 | double Lw;
329          return(ratio/(Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)));
330   }
331  
332 + extern double
333 + crfactor(                       /* contrast reduction factor */
334 +        double  Lw
335 + )
336 + {
337 +        int     i = HISTRES*(Bl(Lw) - bwmin)/(bwmax - bwmin);
338 +        double  bLw, ratio, Tdb;
339  
340 < int
341 < mkbrmap()                       /* make dynamic range map */
340 >        if (i <= 0)
341 >                return(1.0);
342 >        if (i >= HISTRES)
343 >                return(1.0);
344 >        bLw = BLw(Lw);
345 >        ratio = what2do&DO_HSENS ? htcontrs(Lb(bLw))/htcontrs(Lw) : Lb(bLw)/Lw;
346 >        Tdb = mhistot * (bwmax - bwmin) / HISTRES;
347 >        return(modhist[i]*Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)/(Tdb*ratio));
348 > }
349 >
350 >
351 > #if ADJ_VEIL
352 > static void
353 > mkcrfimage(void)                        /* compute contrast reduction factor image */
354   {
355 <        double  T, b, s;
355 >        int     i;
356 >        float   *crfptr;
357 >        COLOR   *fovptr;
358 >
359 >        if (crfimg == NULL)
360 >                crfimg = (float *)malloc(fvxr*fvyr*sizeof(float));
361 >        if (crfimg == NULL)
362 >                syserror("malloc");
363 >        crfptr = crfimg;
364 >        fovptr = fovimg;
365 >        for (i = fvxr*fvyr; i--; crfptr++, fovptr++)
366 >                crfptr[0] = crfactor(plum(fovptr[0]));
367 > }
368 > #endif
369 >
370 >
371 > extern int
372 > mkbrmap(void)                   /* make dynamic range map */
373 > {
374 >        double  Tdb, b, s;
375          double  ceiling, trimmings;
376          register int    i;
377                                          /* copy initial histogram */
378 <        bcopy((char *)bwhist, (char *)modhist, sizeof(modhist));
379 <        s = (bwmax - bwmin)/HISTRES;
378 >        memcpy((void *)modhist, (void *)bwhist, sizeof(modhist));
379 >        s = (bwmax - bwmin)/HISTRES;    /* s is delta b */
380                                          /* loop until satisfactory */
381          do {
382                  mkcumf();                       /* sync brightness mapping */
383                  if (mhistot <= histot*CVRATIO)
384                          return(-1);             /* no compression needed! */
385 <                T = mhistot * (bwmax - bwmin) / HISTRES;
385 >                Tdb = mhistot * s;
386                  trimmings = 0.;                 /* clip to envelope */
387                  for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
388 <                        ceiling = T*clampf(Lb(b));
388 >                        ceiling = Tdb*clampf(Lb(b));
389                          if (modhist[i] > ceiling) {
390                                  trimmings += modhist[i] - ceiling;
391                                  modhist[i] = ceiling;
# Line 330 | Line 393 | mkbrmap()                      /* make dynamic range map */
393                  }
394          } while (trimmings > histot*CVRATIO);
395  
396 + #if ADJ_VEIL
397 +        mkcrfimage();                   /* contrast reduction image */
398 + #endif
399 +
400          return(0);                      /* we got it */
401   }
402  
403  
404 < scotscan(scan, xres)            /* apply scotopic color sensitivity loss */
405 < COLOR   *scan;
406 < int     xres;
404 > extern void
405 > scotscan(               /* apply scotopic color sensitivity loss */
406 >        COLOR   *scan,
407 >        int     xres
408 > )
409   {
410          COLOR   ctmp;
411          double  incolor, b, Lw;
# Line 365 | Line 434 | int    xres;
434   }
435  
436  
437 < mapscan(scan, xres)             /* apply tone mapping operator to scanline */
438 < COLOR   *scan;
439 < int     xres;
437 > extern void
438 > mapscan(                /* apply tone mapping operator to scanline */
439 >        COLOR   *scan,
440 >        int     xres
441 > )
442   {
443          double  mult, Lw, b;
444 <        register int    i;
444 >        register int    x;
445  
446 <        for (i = 0; i < xres; i++) {
447 <                Lw = plum(scan[i]);
446 >        for (x = 0; x < xres; x++) {
447 >                Lw = plum(scan[x]);
448                  if (Lw < LMIN) {
449 <                        setcolor(scan[i], 0., 0., 0.);
449 >                        setcolor(scan[x], 0., 0., 0.);
450                          continue;
451                  }
452 <                b = BLw(Lw);
453 <                mult = (Lb(b) - ldmin)/(ldmax - ldmin) / (Lw*inpexp);
452 >                b = BLw(Lw);            /* apply brightness mapping */
453 >                mult = (Lb(b) - ldmin)/(ldmax - ldmin)/(Lw*inpexp);
454                  if (lumf == rgblum) mult *= WHTEFFICACY;
455 <                scalecolor(scan[i], mult);
455 >                scalecolor(scan[x], mult);
456          }
457   }
458  
459  
460 < putmapping(fp)                  /* put out mapping function */
461 < FILE    *fp;
460 > extern void
461 > putmapping(                     /* put out mapping function */
462 >        FILE    *fp
463 > )
464   {
465          double  b, s;
466          register int    i;

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines