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.20 by greg, Wed Sep 11 18:56:12 2024 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 + #include "random.h"
12  
13  
14   #define CVRATIO         0.025           /* fraction of samples allowed > env. */
# Line 16 | Line 16 | static char SCCSid[] = "$SunId$ LBL";
16   #define LN_10           2.30258509299404568402
17   #define exp10(x)        exp(LN_10*(x))
18  
19 < float   modhist[HISTRES];               /* modified histogram */
19 > double  modhist[HISTRES];               /* modified histogram */
20   double  mhistot;                        /* modified histogram total */
21 < float   cumf[HISTRES+1];                /* cumulative distribution function */
21 > double  cumf[HISTRES+1];                /* cumulative distribution function */
22  
23 + static double centprob(int x, int y);
24 + static void mkcumf(void);
25 + static double cf(double b);
26 + static double BLw(double Lw);
27 + static float *getlumsamp(int n);
28 + #if ADJ_VEIL
29 + static void mkcrfimage(void);
30 + #endif
31  
32 < getfixations(fp)                /* load fixation history list */
33 < FILE    *fp;
32 >
33 >
34 > void
35 > getfixations(           /* load fixation history list */
36 > FILE    *fp
37 > )
38   {
39   #define FIXHUNK         128
40          RESOLU  fvres;
41          int     pos[2];
42 <        register int    px, py, i;
42 >        int     px, py, i;
43                                  /* initialize our resolution struct */
44 <        if ((fvres.or=inpres.or)&YMAJOR) {
44 >        if ((fvres.rt=inpres.rt)&YMAJOR) {
45                  fvres.xr = fvxr;
46                  fvres.yr = fvyr;
47          } else {
# Line 57 | Line 69 | FILE   *fp;
69                                  if (nfixations % FIXHUNK == 0) {
70                                          if (nfixations)
71                                                  fixlst = (short (*)[2])
72 <                                                        realloc((char *)fixlst,
72 >                                                        realloc((void *)fixlst,
73                                                          (nfixations+FIXHUNK)*
74                                                          2*sizeof(short));
75                                          else
# Line 82 | Line 94 | FILE   *fp;
94   }
95  
96  
97 < gethisto(fp)                    /* load precomputed luminance histogram */
98 < FILE    *fp;
97 > void
98 > gethisto(                       /* load precomputed luminance histogram */
99 >        FILE    *fp
100 > )
101   {
102 <        float   histo[MAXPREHIST];
102 >        double  histo[MAXPREHIST];
103          double  histart, histep;
104 <        double  l, b, lastb, w;
104 >        double  b, lastb, w;
105          int     n;
106 <        register int    i;
106 >        int     i;
107                                          /* load data */
108          for (i = 0; i < MAXPREHIST &&
109 <                        fscanf(fp, "%lf %f", &b, &histo[i]) == 2; i++) {
110 <                if (i > 1 && fabs(b - lastb - histep) > 1e-3) {
109 >                        fscanf(fp, "%lf %lf", &b, &histo[i]) == 2; i++) {
110 >                if (i > 1 && fabs(b - lastb - histep) > .001) {
111                          fprintf(stderr,
112                                  "%s: uneven step size in histogram data\n",
113                                          progname);
114                          exit(1);
115                  }
116                  if (i == 1)
117 <                        histep = b - (histart = lastb);
117 >                        if ((histep = b - (histart = lastb)) <= FTINY) {
118 >                                fprintf(stderr,
119 >                                        "%s: illegal step in histogram data\n",
120 >                                                progname);
121 >                                exit(1);
122 >                        }
123                  lastb = b;
124          }
125          if (i < 2 || !feof(fp)) {
# Line 115 | Line 134 | FILE   *fp;
134                                          /* find extrema */
135          for (i = 0; i < n && histo[i] <= FTINY; i++)
136                  ;
137 <        bwmin = histart + i*histep;
137 >        bwmin = histart + (i-.001)*histep;
138          for (i = n; i-- && histo[i] <= FTINY; )
139                  ;
140 <        bwmax = histart + i*histep;
140 >        bwmax = histart + (i+1.001)*histep;
141          if (bwmax > Bl(LMAX))
142                  bwmax = Bl(LMAX);
143          if (bwmin < Bl(LMIN))
# Line 130 | Line 149 | FILE   *fp;
149          for (i = 0; i < HISTRES; i++)
150                  bwhist[i] = 0.;
151          for (i = 0, b = histart; i < n; i++, b += histep) {
152 <                if (b < bwmin) continue;
153 <                if (b > bwmax) break;
152 >                if (b < bwmin+FTINY) continue;
153 >                if (b >= bwmax-FTINY) break;
154                  w = histo[i];
155                  bwavg += w*b;
156                  bwhist[bwhi(b)] += w;
# Line 145 | Line 164 | FILE   *fp;
164   }
165  
166  
167 < double
168 < centprob(x, y)                  /* center-weighting probability function */
169 < int     x, y;
167 > static double
168 > centprob(                       /* center-weighting probability function */
169 >        int     x,
170 >        int     y
171 > )
172   {
173          double  xr, yr, p;
174                                  /* paraboloid, 0 at 90 degrees from center */
# Line 158 | Line 179 | int    x, y;
179   }
180  
181  
182 < comphist()                      /* create foveal sampling histogram */
182 > static float *
183 > getlumsamp(int n)               /* get set of random point sample luminances */
184   {
185 +        float   *ls = (float *)malloc(n*sizeof(float));
186 +        COLR    *cscan = (COLR *)malloc(scanlen(&inpres)*sizeof(COLR));
187 +        long    startpos = ftell(infp);
188 +        long    npleft = (long)inpres.xr*inpres.yr;
189 +        int     x;
190 +        
191 +        if ((ls == NULL) | (cscan == NULL))
192 +                syserror("malloc");
193 +        x = 0;                          /* read/convert samples */
194 +        while (n > 0) {
195 +                COLOR   col;
196 +                int     sv = 0;
197 +                double  rval, cumprob = 0;
198 +
199 +                if (x <= 0 && fread2colrs(cscan, x=scanlen(&inpres), infp,
200 +                                                NCSAMP, WLPART) < 0) {
201 +                        fprintf(stderr, "%s: %s: scanline read error\n",
202 +                                        progname, infn);
203 +                        exit(1);
204 +                }
205 +                rval = frandom();       /* random distance to next sample */
206 +                while ((cumprob += (1.-cumprob)*n/(npleft-sv)) < rval)
207 +                        sv++;
208 +                if (x < ++sv) {         /* out of pixels in this scanline */
209 +                        npleft -= x;
210 +                        x = 0;
211 +                        continue;
212 +                }
213 +                x -= sv;
214 +                colr_color(col, cscan[x]);
215 +                ls[--n] = plum(col);
216 +                npleft -= sv;
217 +        }
218 +        free(cscan);                    /* clean up and reset file pointer */
219 +        if (fseek(infp, startpos, SEEK_SET) < 0)
220 +                syserror("fseek");
221 +        return(ls);
222 + }
223 +
224 +
225 + void
226 + comphist(void)                  /* create foveal sampling histogram */
227 + {
228          double  l, b, w, lwmin, lwmax;
229 <        register int    x, y;
229 >        float   *lumsamp;
230 >        int     nlumsamp;
231 >        int     x, y;
232                                          /* check for precalculated histogram */
233          if (what2do&DO_PREHIST)
234                  return;
# Line 173 | Line 240 | comphist()                     /* create foveal sampling histogram */
240                          if (l < lwmin) lwmin = l;
241                          if (l > lwmax) lwmax = l;
242                  }
243 +                                        /* sample luminance pixels */
244 +        nlumsamp = fvxr*fvyr*16;
245 +        if (nlumsamp > inpres.xr*inpres.yr)
246 +                nlumsamp = inpres.xr*inpres.yr;
247 +        lumsamp = getlumsamp(nlumsamp);
248 +        for (x = nlumsamp; x--; ) {
249 +                l = lumsamp[x];
250 +                if (l < lwmin) lwmin = l;
251 +                if (l > lwmax) lwmax = l;
252 +        }
253          lwmax *= 1.01;
254          if (lwmax > LMAX)
255                  lwmax = LMAX;
# Line 187 | Line 264 | comphist()                     /* create foveal sampling histogram */
264                                          /* (re)compute histogram */
265          bwavg = 0.;
266          histot = 0.;
267 <        for (x = 0; x < HISTRES; x++)
191 <                bwhist[x] = 0.;
267 >        memset(bwhist, 0, sizeof(bwhist));
268                                          /* global average */
269 <        if (!(what2do&DO_FIXHIST) || fixfrac < 1.-FTINY)
269 >        if (!(what2do&DO_FIXHIST) || fixfrac < 1.-FTINY) {
270                  for (y = 0; y < fvyr; y++)
271                          for (x = 0; x < fvxr; x++) {
272                                  l = plum(fovscan(y)[x]);
273 <                                if (l < lwmin) continue;
274 <                                if (l > lwmax) continue;
273 >                                if (l < lwmin+FTINY) continue;
274 >                                if (l >= lwmax-FTINY) continue;
275                                  b = Bl(l);
276                                  w = what2do&DO_CWEIGHT ? centprob(x,y) : 1.;
277                                  bwavg += w*b;
278                                  bwhist[bwhi(b)] += w;
279                                  histot += w;
280                          }
281 +                                        /* weight for point luminances */
282 +                w = 1. * histot / nlumsamp;
283 +                for (x = nlumsamp; x--; ) {
284 +                        l = lumsamp[x];
285 +                        if (l < lwmin+FTINY) continue;
286 +                        if (l >= lwmax-FTINY) continue;
287 +                        b = Bl(l);
288 +                        bwavg += w*b;
289 +                        bwhist[bwhi(b)] += w;
290 +                        histot += w;
291 +                }
292 +        }
293                                          /* average fixation points */
294          if (what2do&DO_FIXHIST && nfixations > 0) {
295                  if (histot > FTINY)
# Line 210 | Line 298 | comphist()                     /* create foveal sampling histogram */
298                          w = 1.;
299                  for (x = 0; x < nfixations; x++) {
300                          l = plum(fovscan(fixlst[x][1])[fixlst[x][0]]);
301 <                        if (l < lwmin) continue;
302 <                        if (l > lwmax) continue;
301 >                        if (l < lwmin+FTINY) continue;
302 >                        if (l >= lwmax-FTINY) continue;
303                          b = Bl(l);
304                          bwavg += w*b;
305                          bwhist[bwhi(b)] += w;
# Line 223 | Line 311 | comphist()                     /* create foveal sampling histogram */
311                  bwhist[1] *= 0.5;
312                  bwhist[0] += bwhist[1];
313          }
314 +        free(lumsamp);
315   }
316  
317  
318 < mkcumf()                        /* make cumulative distribution function */
318 > static void
319 > mkcumf(void)                    /* make cumulative distribution function */
320   {
321 <        register int    i;
322 <        register double sum;
321 >        int     i;
322 >        double  sum;
323  
324          mhistot = 0.;           /* compute modified total */
325          for (i = 0; i < HISTRES; i++)
# Line 244 | Line 334 | mkcumf()                       /* make cumulative distribution function */
334   }
335  
336  
337 < double
338 < cf(b)                           /* return cumulative function at b */
339 < double  b;
337 > static double
338 > cf(                             /* return cumulative function at b */
339 >        double  b
340 > )
341   {
342          double  x;
343 <        register int    i;
343 >        int     i;
344  
345          i = x = HISTRES*(b - bwmin)/(bwmax - bwmin);
346          x -= (double)i;
# Line 257 | Line 348 | double b;
348   }
349  
350  
351 < double
352 < BLw(Lw)                         /* map world luminance to display brightness */
353 < double  Lw;
351 > static double
352 > BLw(                            /* map world luminance to display brightness */
353 >        double  Lw
354 > )
355   {
356          double  b;
357  
# Line 272 | Line 364 | double Lw;
364  
365  
366   double
367 < htcontrs(La)            /* human threshold contrast sensitivity, dL(La) */
368 < double  La;
367 > htcontrs(               /* human threshold contrast sensitivity, dL(La) */
368 >        double  La
369 > )
370   {
371          double  l10La, l10dL;
372                                  /* formula taken from Ferwerda et al. [SG96] */
# Line 294 | Line 387 | double La;
387  
388  
389   double
390 < clampf(Lw)              /* derivative clamping function */
391 < double  Lw;
390 > clampf(                 /* histogram clamping function */
391 >        double  Lw
392 > )
393   {
394          double  bLw, ratio;
395  
# Line 304 | Line 398 | double Lw;
398          return(ratio/(Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)));
399   }
400  
401 + double
402 + crfactor(                       /* contrast reduction factor */
403 +        double  Lw
404 + )
405 + {
406 +        int     i = HISTRES*(Bl(Lw) - bwmin)/(bwmax - bwmin);
407 +        double  bLw, ratio, Tdb;
408  
409 +        if (i <= 0)
410 +                return(1.0);
411 +        if (i >= HISTRES)
412 +                return(1.0);
413 +        bLw = BLw(Lw);
414 +        ratio = what2do&DO_HSENS ? htcontrs(Lb(bLw))/htcontrs(Lw) : Lb(bLw)/Lw;
415 +        Tdb = mhistot * (bwmax - bwmin) / HISTRES;
416 +        return(modhist[i]*Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)/(Tdb*ratio));
417 + }
418 +
419 +
420 + #if ADJ_VEIL
421 + static void
422 + mkcrfimage(void)                        /* compute contrast reduction factor image */
423 + {
424 +        int     i;
425 +        float   *crfptr;
426 +        COLOR   *fovptr;
427 +
428 +        if (crfimg == NULL)
429 +                crfimg = (float *)malloc(fvxr*fvyr*sizeof(float));
430 +        if (crfimg == NULL)
431 +                syserror("malloc");
432 +        crfptr = crfimg;
433 +        fovptr = fovimg;
434 +        for (i = fvxr*fvyr; i--; crfptr++, fovptr++)
435 +                crfptr[0] = crfactor(plum(fovptr[0]));
436 + }
437 + #endif
438 +
439 +
440   int
441 < mkbrmap()                       /* make dynamic range map */
441 > mkbrmap(void)                   /* make dynamic range map */
442   {
443 <        double  T, b, s;
443 >        double  Tdb, b, s;
444          double  ceiling, trimmings;
445 <        register int    i;
445 >        int     i;
446                                          /* copy initial histogram */
447 <        bcopy((char *)bwhist, (char *)modhist, sizeof(modhist));
448 <        s = (bwmax - bwmin)/HISTRES;
447 >        memcpy((void *)modhist, (void *)bwhist, sizeof(modhist));
448 >        s = (bwmax - bwmin)/HISTRES;    /* s is delta b */
449                                          /* loop until satisfactory */
450          do {
451                  mkcumf();                       /* sync brightness mapping */
452                  if (mhistot <= histot*CVRATIO)
453                          return(-1);             /* no compression needed! */
454 <                T = mhistot * (bwmax - bwmin) / HISTRES;
454 >                Tdb = mhistot * s;
455                  trimmings = 0.;                 /* clip to envelope */
456                  for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
457 <                        ceiling = T*clampf(Lb(b));
457 >                        ceiling = Tdb*clampf(Lb(b));
458                          if (modhist[i] > ceiling) {
459                                  trimmings += modhist[i] - ceiling;
460                                  modhist[i] = ceiling;
461                          }
462                  }
463 <        } while (trimmings > histot*CVRATIO);
463 >        } while (trimmings > mhistot*CVRATIO);
464  
465 + #if ADJ_VEIL
466 +        mkcrfimage();                   /* contrast reduction image */
467 + #endif
468 +
469          return(0);                      /* we got it */
470   }
471  
472  
473 < scotscan(scan, xres)            /* apply scotopic color sensitivity loss */
474 < COLOR   *scan;
475 < int     xres;
473 > void
474 > scotscan(               /* apply scotopic color sensitivity loss */
475 >        COLOR   *scan,
476 >        int     xres
477 > )
478   {
479          COLOR   ctmp;
480          double  incolor, b, Lw;
481 <        register int    i;
481 >        int     i;
482  
483          for (i = 0; i < xres; i++) {
484                  Lw = plum(scan[i]);
# Line 365 | Line 503 | int    xres;
503   }
504  
505  
506 < mapscan(scan, xres)             /* apply tone mapping operator to scanline */
507 < COLOR   *scan;
508 < int     xres;
506 > void
507 > mapscan(                /* apply tone mapping operator to scanline */
508 >        COLOR   *scan,
509 >        int     xres
510 > )
511   {
512          double  mult, Lw, b;
513 <        register int    i;
513 >        int     x;
514  
515 <        for (i = 0; i < xres; i++) {
516 <                Lw = plum(scan[i]);
515 >        for (x = 0; x < xres; x++) {
516 >                Lw = plum(scan[x]);
517                  if (Lw < LMIN) {
518 <                        setcolor(scan[i], 0., 0., 0.);
518 >                        setcolor(scan[x], 0., 0., 0.);
519                          continue;
520                  }
521 <                b = BLw(Lw);
522 <                mult = (Lb(b) - ldmin)/(ldmax - ldmin) / (Lw*inpexp);
521 >                b = BLw(Lw);            /* apply brightness mapping */
522 >                mult = (Lb(b) - ldmin)/(ldmax - ldmin)/(Lw*inpexp);
523                  if (lumf == rgblum) mult *= WHTEFFICACY;
524 <                scalecolor(scan[i], mult);
524 >                scalecolor(scan[x], mult);
525          }
526   }
527  
528  
529 < putmapping(fp)                  /* put out mapping function */
530 < FILE    *fp;
529 > void
530 > putmapping(                     /* put out mapping function */
531 >        FILE    *fp
532 > )
533   {
534          double  b, s;
535 <        register int    i;
535 >        int     i;
536          double  wlum, sf, dlum;
537  
538          sf = scalef*inpexp;

Diff Legend

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