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.7 by greg, Sat Jan 11 10:39:25 1997 UTC vs.
Revision 3.18 by greg, Sat Jun 27 03:37:24 2020 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. */
15  
16 < #define exp10(x)        exp(2.302585093*(x))
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 < mkcumf()                        /* make cumulative distribution function */
32 >
33 >
34 > void
35 > getfixations(           /* load fixation history list */
36 > FILE    *fp
37 > )
38   {
39 <        register int    i;
40 <        register double sum;
39 > #define FIXHUNK         128
40 >        RESOLU  fvres;
41 >        int     pos[2];
42 >        int     px, py, i;
43 >                                /* initialize our resolution struct */
44 >        if ((fvres.rt=inpres.rt)&YMAJOR) {
45 >                fvres.xr = fvxr;
46 >                fvres.yr = fvyr;
47 >        } else {
48 >                fvres.xr = fvyr;
49 >                fvres.yr = fvxr;
50 >        }
51 >                                /* read each picture position */
52 >        while (fscanf(fp, "%d %d", &pos[0], &pos[1]) == 2) {
53 >                                /* convert to closest index in foveal image */
54 >                loc2pix(pos, &fvres,
55 >                                (pos[0]+.5)/inpres.xr, (pos[1]+.5)/inpres.yr);
56 >                                /* include nine neighborhood samples */
57 >                for (px = pos[0]-1; px <= pos[0]+1; px++) {
58 >                        if (px < 0 || px >= fvxr)
59 >                                continue;
60 >                        for (py = pos[1]-1; py <= pos[1]+1; py++) {
61 >                                if (py < 0 || py >= fvyr)
62 >                                        continue;
63 >                                for (i = nfixations; i-- > 0; )
64 >                                        if (fixlst[i][0] == px &&
65 >                                                        fixlst[i][1] == py)
66 >                                                break;
67 >                                if (i >= 0)
68 >                                        continue;       /* already there */
69 >                                if (nfixations % FIXHUNK == 0) {
70 >                                        if (nfixations)
71 >                                                fixlst = (short (*)[2])
72 >                                                        realloc((void *)fixlst,
73 >                                                        (nfixations+FIXHUNK)*
74 >                                                        2*sizeof(short));
75 >                                        else
76 >                                                fixlst = (short (*)[2])malloc(
77 >                                                        FIXHUNK*2*sizeof(short)
78 >                                                        );
79 >                                        if (fixlst == NULL)
80 >                                                syserror("malloc");
81 >                                }
82 >                                fixlst[nfixations][0] = px;
83 >                                fixlst[nfixations][1] = py;
84 >                                nfixations++;
85 >                        }
86 >                }
87 >        }
88 >        if (!feof(fp)) {
89 >                fprintf(stderr, "%s: format error reading fixation data\n",
90 >                                progname);
91 >                exit(1);
92 >        }
93 > #undef  FIXHUNK
94 > }
95  
96 +
97 + void
98 + gethisto(                       /* load precomputed luminance histogram */
99 +        FILE    *fp
100 + )
101 + {
102 +        double  histo[MAXPREHIST];
103 +        double  histart, histep;
104 +        double  b, lastb, w;
105 +        int     n;
106 +        int     i;
107 +                                        /* load data */
108 +        for (i = 0; i < MAXPREHIST &&
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 +                        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)) {
126 +                fprintf(stderr,
127 +                "%s: format/length error loading histogram (log10L %f at %d)\n",
128 +                                progname, b, i);
129 +                exit(1);
130 +        }
131 +        n = i;
132 +        histart *= LN_10;
133 +        histep *= LN_10;
134 +                                        /* find extrema */
135 +        for (i = 0; i < n && histo[i] <= FTINY; i++)
136 +                ;
137 +        bwmin = histart + (i-.001)*histep;
138 +        for (i = n; i-- && histo[i] <= FTINY; )
139 +                ;
140 +        bwmax = histart + (i+1.001)*histep;
141 +        if (bwmax > Bl(LMAX))
142 +                bwmax = Bl(LMAX);
143 +        if (bwmin < Bl(LMIN))
144 +                bwmin = Bl(LMIN);
145 +        else                            /* duplicate bottom bin */
146 +                bwmin = bwmax - (bwmax-bwmin)*HISTRES/(HISTRES-1);
147 +                                        /* convert histogram */
148 +        bwavg = 0.; histot = 0.;
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+FTINY) continue;
153 +                if (b >= bwmax-FTINY) break;
154 +                w = histo[i];
155 +                bwavg += w*b;
156 +                bwhist[bwhi(b)] += w;
157 +                histot += w;
158 +        }
159 +        bwavg /= histot;
160 +        if (bwmin > Bl(LMIN)+FTINY) {   /* add false samples at bottom */
161 +                bwhist[1] *= 0.5;
162 +                bwhist[0] += bwhist[1];
163 +        }
164 + }
165 +
166 +
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 */
175 +        xr = (x - .5*(fvxr-1))/90.;     /* 180 degree fisheye has fv?r == 90 */
176 +        yr = (y - .5*(fvyr-1))/90.;
177 +        p = 1. - xr*xr - yr*yr;
178 +        return(p < 0. ? 0. : p);
179 + }
180 +
181 +
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 && freadcolrs(cscan, x=scanlen(&inpres), infp) < 0) {
200 +                        fprintf(stderr, "%s: %s: scanline read error\n",
201 +                                        progname, infn);
202 +                        exit(1);
203 +                }
204 +                rval = frandom();       /* random distance to next sample */
205 +                while ((cumprob += (1.-cumprob)*n/(npleft-sv)) < rval)
206 +                        sv++;
207 +                if (x < ++sv) {         /* out of pixels in this scanline */
208 +                        npleft -= x;
209 +                        x = 0;
210 +                        continue;
211 +                }
212 +                x -= sv;
213 +                colr_color(col, cscan[x]);
214 +                ls[--n] = plum(col);
215 +                npleft -= sv;
216 +        }
217 +        free(cscan);                    /* clean up and reset file pointer */
218 +        if (fseek(infp, startpos, SEEK_SET) < 0)
219 +                syserror("fseek");
220 +        return(ls);
221 + }
222 +
223 +
224 + void
225 + comphist(void)                  /* create foveal sampling histogram */
226 + {
227 +        double  l, b, w, lwmin, lwmax;
228 +        float   *lumsamp;
229 +        int     nlumsamp;
230 +        int     x, y;
231 +                                        /* check for precalculated histogram */
232 +        if (what2do&DO_PREHIST)
233 +                return;
234 +        lwmin = 1e10;                   /* find extrema */
235 +        lwmax = 0.;
236 +        for (y = 0; y < fvyr; y++)
237 +                for (x = 0; x < fvxr; x++) {
238 +                        l = plum(fovscan(y)[x]);
239 +                        if (l < lwmin) lwmin = l;
240 +                        if (l > lwmax) lwmax = l;
241 +                }
242 +                                        /* sample luminance pixels */
243 +        nlumsamp = fvxr*fvyr*16;
244 +        if (nlumsamp > inpres.xr*inpres.yr)
245 +                nlumsamp = inpres.xr*inpres.yr;
246 +        lumsamp = getlumsamp(nlumsamp);
247 +        for (x = nlumsamp; x--; ) {
248 +                l = lumsamp[x];
249 +                if (l < lwmin) lwmin = l;
250 +                if (l > lwmax) lwmax = l;
251 +        }
252 +        lwmax *= 1.01;
253 +        if (lwmax > LMAX)
254 +                lwmax = LMAX;
255 +        bwmax = Bl(lwmax);
256 +        if (lwmin < LMIN) {
257 +                lwmin = LMIN;
258 +                bwmin = Bl(LMIN);
259 +        } else {                        /* duplicate bottom bin */
260 +                bwmin = bwmax - (bwmax-Bl(lwmin))*HISTRES/(HISTRES-1);
261 +                lwmin = Lb(bwmin);
262 +        }
263 +                                        /* (re)compute histogram */
264 +        bwavg = 0.;
265 +        histot = 0.;
266 +        memset(bwhist, 0, sizeof(bwhist));
267 +                                        /* global average */
268 +        if (!(what2do&DO_FIXHIST) || fixfrac < 1.-FTINY) {
269 +                for (y = 0; y < fvyr; y++)
270 +                        for (x = 0; x < fvxr; x++) {
271 +                                l = plum(fovscan(y)[x]);
272 +                                if (l < lwmin+FTINY) continue;
273 +                                if (l >= lwmax-FTINY) continue;
274 +                                b = Bl(l);
275 +                                w = what2do&DO_CWEIGHT ? centprob(x,y) : 1.;
276 +                                bwavg += w*b;
277 +                                bwhist[bwhi(b)] += w;
278 +                                histot += w;
279 +                        }
280 +                                        /* weight for point luminances */
281 +                w = 1. * histot / nlumsamp;
282 +                for (x = nlumsamp; x--; ) {
283 +                        l = lumsamp[x];
284 +                        if (l < lwmin+FTINY) continue;
285 +                        if (l >= lwmax-FTINY) continue;
286 +                        b = Bl(l);
287 +                        bwavg += w*b;
288 +                        bwhist[bwhi(b)] += w;
289 +                        histot += w;
290 +                }
291 +        }
292 +                                        /* average fixation points */
293 +        if (what2do&DO_FIXHIST && nfixations > 0) {
294 +                if (histot > FTINY)
295 +                        w = fixfrac/(1.-fixfrac)*histot/nfixations;
296 +                else
297 +                        w = 1.;
298 +                for (x = 0; x < nfixations; x++) {
299 +                        l = plum(fovscan(fixlst[x][1])[fixlst[x][0]]);
300 +                        if (l < lwmin+FTINY) continue;
301 +                        if (l >= lwmax-FTINY) continue;
302 +                        b = Bl(l);
303 +                        bwavg += w*b;
304 +                        bwhist[bwhi(b)] += w;
305 +                        histot += w;
306 +                }
307 +        }
308 +        bwavg /= histot;
309 +        if (lwmin > LMIN+FTINY) {       /* add false samples at bottom */
310 +                bwhist[1] *= 0.5;
311 +                bwhist[0] += bwhist[1];
312 +        }
313 +        free(lumsamp);
314 + }
315 +
316 +
317 + static void
318 + mkcumf(void)                    /* make cumulative distribution function */
319 + {
320 +        int     i;
321 +        double  sum;
322 +
323          mhistot = 0.;           /* compute modified total */
324          for (i = 0; i < HISTRES; i++)
325                  mhistot += modhist[i];
# Line 38 | Line 333 | mkcumf()                       /* make cumulative distribution function */
333   }
334  
335  
336 < double
337 < cf(b)                           /* return cumulative function at b */
338 < double  b;
336 > static double
337 > cf(                             /* return cumulative function at b */
338 >        double  b
339 > )
340   {
341          double  x;
342 <        register int    i;
342 >        int     i;
343  
344          i = x = HISTRES*(b - bwmin)/(bwmax - bwmin);
345          x -= (double)i;
# Line 51 | Line 347 | double b;
347   }
348  
349  
350 < double
351 < BLw(Lw)                         /* map world luminance to display brightness */
352 < double  Lw;
350 > static double
351 > BLw(                            /* map world luminance to display brightness */
352 >        double  Lw
353 > )
354   {
355          double  b;
356  
# Line 66 | Line 363 | double Lw;
363  
364  
365   double
366 < htcontrs(La)            /* human threshold contrast sensitivity, dL(La) */
367 < double  La;
366 > htcontrs(               /* human threshold contrast sensitivity, dL(La) */
367 >        double  La
368 > )
369   {
370          double  l10La, l10dL;
371                                  /* formula taken from Ferwerda et al. [SG96] */
# Line 88 | Line 386 | double La;
386  
387  
388   double
389 < clampf(Lw)              /* derivative clamping function */
390 < double  Lw;
389 > clampf(                 /* histogram clamping function */
390 >        double  Lw
391 > )
392   {
393          double  bLw, ratio;
394  
# Line 98 | Line 397 | double Lw;
397          return(ratio/(Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)));
398   }
399  
400 + double
401 + crfactor(                       /* contrast reduction factor */
402 +        double  Lw
403 + )
404 + {
405 +        int     i = HISTRES*(Bl(Lw) - bwmin)/(bwmax - bwmin);
406 +        double  bLw, ratio, Tdb;
407  
408 +        if (i <= 0)
409 +                return(1.0);
410 +        if (i >= HISTRES)
411 +                return(1.0);
412 +        bLw = BLw(Lw);
413 +        ratio = what2do&DO_HSENS ? htcontrs(Lb(bLw))/htcontrs(Lw) : Lb(bLw)/Lw;
414 +        Tdb = mhistot * (bwmax - bwmin) / HISTRES;
415 +        return(modhist[i]*Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)/(Tdb*ratio));
416 + }
417 +
418 +
419 + #if ADJ_VEIL
420 + static void
421 + mkcrfimage(void)                        /* compute contrast reduction factor image */
422 + {
423 +        int     i;
424 +        float   *crfptr;
425 +        COLOR   *fovptr;
426 +
427 +        if (crfimg == NULL)
428 +                crfimg = (float *)malloc(fvxr*fvyr*sizeof(float));
429 +        if (crfimg == NULL)
430 +                syserror("malloc");
431 +        crfptr = crfimg;
432 +        fovptr = fovimg;
433 +        for (i = fvxr*fvyr; i--; crfptr++, fovptr++)
434 +                crfptr[0] = crfactor(plum(fovptr[0]));
435 + }
436 + #endif
437 +
438 +
439   int
440 < mkbrmap()                       /* make dynamic range map */
440 > mkbrmap(void)                   /* make dynamic range map */
441   {
442 <        double  T, b, s;
442 >        double  Tdb, b, s;
443          double  ceiling, trimmings;
444 <        register int    i;
444 >        int     i;
445                                          /* copy initial histogram */
446 <        bcopy((char *)bwhist, (char *)modhist, sizeof(modhist));
447 <        s = (bwmax - bwmin)/HISTRES;
446 >        memcpy((void *)modhist, (void *)bwhist, sizeof(modhist));
447 >        s = (bwmax - bwmin)/HISTRES;    /* s is delta b */
448                                          /* loop until satisfactory */
449          do {
450                  mkcumf();                       /* sync brightness mapping */
451                  if (mhistot <= histot*CVRATIO)
452                          return(-1);             /* no compression needed! */
453 <                T = mhistot * (bwmax - bwmin) / HISTRES;
453 >                Tdb = mhistot * s;
454                  trimmings = 0.;                 /* clip to envelope */
455                  for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
456 <                        ceiling = T*clampf(Lb(b));
456 >                        ceiling = Tdb*clampf(Lb(b));
457                          if (modhist[i] > ceiling) {
458                                  trimmings += modhist[i] - ceiling;
459                                  modhist[i] = ceiling;
# Line 124 | Line 461 | mkbrmap()                      /* make dynamic range map */
461                  }
462          } while (trimmings > histot*CVRATIO);
463  
464 + #if ADJ_VEIL
465 +        mkcrfimage();                   /* contrast reduction image */
466 + #endif
467 +
468          return(0);                      /* we got it */
469   }
470  
471  
472 < scotscan(scan, xres)            /* apply scotopic color sensitivity loss */
473 < COLOR   *scan;
474 < int     xres;
472 > void
473 > scotscan(               /* apply scotopic color sensitivity loss */
474 >        COLOR   *scan,
475 >        int     xres
476 > )
477   {
478          COLOR   ctmp;
479          double  incolor, b, Lw;
480 <        register int    i;
480 >        int     i;
481  
482          for (i = 0; i < xres; i++) {
483                  Lw = plum(scan[i]);
# Line 159 | Line 502 | int    xres;
502   }
503  
504  
505 < mapscan(scan, xres)             /* apply tone mapping operator to scanline */
506 < COLOR   *scan;
507 < int     xres;
505 > void
506 > mapscan(                /* apply tone mapping operator to scanline */
507 >        COLOR   *scan,
508 >        int     xres
509 > )
510   {
511          double  mult, Lw, b;
512 <        register int    i;
512 >        int     x;
513  
514 <        for (i = 0; i < xres; i++) {
515 <                Lw = plum(scan[i]);
514 >        for (x = 0; x < xres; x++) {
515 >                Lw = plum(scan[x]);
516                  if (Lw < LMIN) {
517 <                        setcolor(scan[i], 0., 0., 0.);
517 >                        setcolor(scan[x], 0., 0., 0.);
518                          continue;
519                  }
520 <                b = BLw(Lw);
521 <                mult = (Lb(b) - ldmin)/(ldmax - ldmin) / (Lw*inpexp);
520 >                b = BLw(Lw);            /* apply brightness mapping */
521 >                mult = (Lb(b) - ldmin)/(ldmax - ldmin)/(Lw*inpexp);
522                  if (lumf == rgblum) mult *= WHTEFFICACY;
523 <                scalecolor(scan[i], mult);
523 >                scalecolor(scan[x], mult);
524          }
525   }
526  
527  
528 < putmapping(fp)                  /* put out mapping function */
529 < FILE    *fp;
528 > void
529 > putmapping(                     /* put out mapping function */
530 >        FILE    *fp
531 > )
532   {
533          double  b, s;
534 <        register int    i;
534 >        int     i;
535          double  wlum, sf, dlum;
536  
537          sf = scalef*inpexp;

Diff Legend

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