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.1 by greg, Thu Oct 3 16:52:50 1996 UTC vs.
Revision 3.20 by greg, Wed Sep 11 18:56:12 2024 UTC

# Line 1 | Line 1
1 /* Copyright (c) 1996 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 pixels allowed > env. */
14 > #define CVRATIO         0.025           /* fraction of samples allowed > env. */
15  
16 < #define BotMesopic      5.62e-3         /* top of scotopic range */
17 < #define TopMesopic      5.62            /* bottom of photopic range */
16 > #define LN_10           2.30258509299404568402
17 > #define exp10(x)        exp(LN_10*(x))
18  
19 < #define exp10(x)        exp(2.302585093*(x))
19 > double  modhist[HISTRES];               /* modified histogram */
20 > double  mhistot;                        /* modified histogram total */
21 > double  cumf[HISTRES+1];                /* cumulative distribution function */
22  
23 < int     modhist[HISTRES];               /* modified histogram */
24 < float   cumf[HISTRES+1];                /* cumulative distribution function */
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  
33 < mkcumf()                        /* make cumulative distribution function */
33 >
34 > void
35 > getfixations(           /* load fixation history list */
36 > FILE    *fp
37 > )
38   {
39 <        register int    i;
40 <        register long   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 <        cumf[0] = 0.;
97 <        sum = modhist[0];
98 <        for (i = 1; i < HISTRES; i++) {
99 <                cumf[i] = (double)sum/histot;
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 && 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 >        float   *lumsamp;
230 >        int     nlumsamp;
231 >        int     x, y;
232 >                                        /* check for precalculated histogram */
233 >        if (what2do&DO_PREHIST)
234 >                return;
235 >        lwmin = 1e10;                   /* find extrema */
236 >        lwmax = 0.;
237 >        for (y = 0; y < fvyr; y++)
238 >                for (x = 0; x < fvxr; x++) {
239 >                        l = plum(fovscan(y)[x]);
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;
256 >        bwmax = Bl(lwmax);
257 >        if (lwmin < LMIN) {
258 >                lwmin = LMIN;
259 >                bwmin = Bl(LMIN);
260 >        } else {                        /* duplicate bottom bin */
261 >                bwmin = bwmax - (bwmax-Bl(lwmin))*HISTRES/(HISTRES-1);
262 >                lwmin = Lb(bwmin);
263 >        }
264 >                                        /* (re)compute histogram */
265 >        bwavg = 0.;
266 >        histot = 0.;
267 >        memset(bwhist, 0, sizeof(bwhist));
268 >                                        /* global average */
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+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)
296 >                        w = fixfrac/(1.-fixfrac)*histot/nfixations;
297 >                else
298 >                        w = 1.;
299 >                for (x = 0; x < nfixations; x++) {
300 >                        l = plum(fovscan(fixlst[x][1])[fixlst[x][0]]);
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;
306 >                        histot += w;
307 >                }
308 >        }
309 >        bwavg /= histot;
310 >        if (lwmin > LMIN+FTINY) {       /* add false samples at bottom */
311 >                bwhist[1] *= 0.5;
312 >                bwhist[0] += bwhist[1];
313 >        }
314 >        free(lumsamp);
315 > }
316 >
317 >
318 > static void
319 > mkcumf(void)                    /* make cumulative distribution function */
320 > {
321 >        int     i;
322 >        double  sum;
323 >
324 >        mhistot = 0.;           /* compute modified total */
325 >        for (i = 0; i < HISTRES; i++)
326 >                mhistot += modhist[i];
327 >
328 >        sum = 0.;               /* compute cumulative function */
329 >        for (i = 0; i < HISTRES; i++) {
330 >                cumf[i] = sum/mhistot;
331                  sum += modhist[i];
332          }
333          cumf[HISTRES] = 1.;
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 50 | 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 60 | Line 359 | double Lw;
359                  return(Bldmin);
360          if (b >= bwmax-FTINY)
361                  return(Bldmax);
362 <        return(Bldmin + cf(Bl(Lw))*(Bldmax-Bldmin));
362 >        return(Bldmin + cf(b)*(Bldmax-Bldmin));
363   }
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 87 | 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 97 | 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 < int
410 < shiftdir(bw)            /* compute shift direction for histogram */
411 < double  bw;
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 <        if (what2do&DO_HSENS && cf(bw) - (bw - bwmin)/(Bldmax - bwmin))
425 <                return(1);
426 <        return(-1);
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 <        int     hdiffs[HISTRES], above, below;
444 <        double  T, b, s;
445 <        int     maxd, maxi, sd;
117 <        register int    i;
443 >        double  Tdb, b, s;
444 >        double  ceiling, trimmings;
445 >        int     i;
446                                          /* copy initial histogram */
447 <        for (i = 0; i < HISTRES; i++)
448 <                modhist[i] = bwhist[i];
121 <        T = histot * (bwmax - bwmin) / HISTRES;
122 <        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 <        for ( ; ; ) {
451 <                mkcumf();               /* sync brightness mapping */
452 <                above = below = 0;      /* compute visibility overflow */
450 >        do {
451 >                mkcumf();                       /* sync brightness mapping */
452 >                if (mhistot <= histot*CVRATIO)
453 >                        return(-1);             /* no compression needed! */
454 >                Tdb = mhistot * s;
455 >                trimmings = 0.;                 /* clip to envelope */
456                  for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
457 <                        hdiffs[i] = modhist[i] - (int)(T*clampf(Lb(b)) + .5);
458 <                        if (hdiffs[i] > 0) above += hdiffs[i];
459 <                        else below -= hdiffs[i];
457 >                        ceiling = Tdb*clampf(Lb(b));
458 >                        if (modhist[i] > ceiling) {
459 >                                trimmings += modhist[i] - ceiling;
460 >                                modhist[i] = ceiling;
461 >                        }
462                  }
463 <                if (above <= histot*CVRATIO)
464 <                        break;          /* close enough */
465 <                if (above-below >= 0)
466 <                        return(-1);     /* Houston, we have a problem.... */
467 <                /* original looped here as well (BEGIN_L2) */
468 <                maxd = 0;               /* find largest overvis */
469 <                for (i = 0; i < HISTRES; i++)
139 <                        if (hdiffs[i] > maxd)
140 <                                maxd = hdiffs[maxi=i];
141 <                /* broke loop here when (maxd == 0) (BREAK_L2) */
142 <                for (sd = shiftdir((maxi+.5)/HISTRES*(bwmax-bwmin)+bwmin);
143 <                                hdiffs[maxi] == maxd; sd = -sd)
144 <                        for (i = maxi+sd; i >= 0 & i < HISTRES; i += sd)
145 <                                if (hdiffs[i] < 0) {
146 <                                        if (hdiffs[i] <= -maxd) {
147 <                                                modhist[i] += maxd;
148 <                                                modhist[maxi] -= maxd;
149 <                                                hdiffs[i] += maxd;
150 <                                                hdiffs[maxi] = 0;
151 <                                        } else {
152 <                                                modhist[maxi] += hdiffs[i];
153 <                                                modhist[i] -= hdiffs[i];
154 <                                                hdiffs[maxi] += hdiffs[i];
155 <                                                hdiffs[i] = 0;
156 <                                        }
157 <                                        break;
158 <                                }
159 <                /* (END_L2) */
160 <        }
161 <        return(0);
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 180 | Line 490 | int    xres;
490                          incolor = (Lw - BotMesopic) /
491                                          (TopMesopic - BotMesopic);
492                  if (incolor < 1.-FTINY) {
493 +                        b = (1.-incolor)*slum(scan[i])*inpexp/SWNORM;
494 +                        if (lumf == rgblum) b /= WHTEFFICACY;
495 +                        setcolor(ctmp, b, b, b);
496                          if (incolor <= FTINY)
497                                  setcolor(scan[i], 0., 0., 0.);
498                          else
499                                  scalecolor(scan[i], incolor);
187                        b = (1.-incolor)*slum(scan[i])*inpexp/SWNORM;
188                        if (lumf == rgblum) b /= WHTEFFICACY;
189                        setcolor(ctmp, b, b, b);
500                          addcolor(scan[i], ctmp);
501                  }
502          }
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 < #ifdef DEBUG
530 < doplots()                       /* generate debugging plots */
529 > void
530 > putmapping(                     /* put out mapping function */
531 >        FILE    *fp
532 > )
533   {
534 <        double  T, b, s;
535 <        FILE    *fp;
536 <        char    fname[128];
223 <        register int    i;
534 >        double  b, s;
535 >        int     i;
536 >        double  wlum, sf, dlum;
537  
538 <        T = histot * (bwmax - bwmin) / HISTRES;
538 >        sf = scalef*inpexp;
539 >        if (lumf == cielum) sf *= WHTEFFICACY;
540          s = (bwmax - bwmin)/HISTRES;
541 <
542 <        sprintf(fname, "%s_hist.plt", infn);
543 <        if ((fp = fopen(fname, "w")) == NULL)
544 <                syserror(fname);
545 <        fputs("include=curve.plt\n", fp);
546 <        fputs("title=\"Brightness Frequency Distribution\"\n", fp);
547 <        fprintf(fp, "subtitle=%s\n", infn);
548 <        fputs("ymin=0\n", fp);
549 <        fputs("xlabel=\"Perceptual Brightness B(Lw)\"\n", fp);
550 <        fputs("ylabel=\"Frequency Count\"\n", fp);
237 <        fputs("Alabel=\"Histogram\"\n", fp);
238 <        fputs("Alintype=0\n", fp);
239 <        fputs("Blabel=\"Envelope\"\n", fp);
240 <        fputs("Bsymsize=0\n", fp);
241 <        fputs("Adata=\n", fp);
242 <        for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s)
243 <                fprintf(fp, "\t%f %d\n", b, modhist[i]);
244 <        fputs(";\nBdata=\n", fp);
245 <        for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s)
246 <                fprintf(fp, "\t%f %f\n", b, T*clampf(Lb(b)));
247 <        fputs(";\n", fp);
248 <        fclose(fp);
541 >        for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
542 >                wlum = Lb(b);
543 >                if (what2do&DO_LINEAR) {
544 >                        dlum = sf*wlum;
545 >                        if (dlum > ldmax) dlum = ldmax;
546 >                        else if (dlum < ldmin) dlum = ldmin;
547 >                        fprintf(fp, "%e %e\n", wlum, dlum);
548 >                } else
549 >                        fprintf(fp, "%e %e\n", wlum, Lb(BLw(wlum)));
550 >        }
551   }
250 #endif

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines