Back to index

radiance  4R0+20100331
pcond3.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: pcond3.c,v 3.15 2004/03/28 20:33:14 schorsch Exp $";
00003 #endif
00004 /*
00005  * Routines for computing and applying brightness mapping.
00006  */
00007 
00008 #include <string.h>
00009 
00010 #include "pcond.h"
00011 
00012 
00013 #define CVRATIO             0.025         /* fraction of samples allowed > env. */
00014 
00015 #define LN_10        2.30258509299404568402
00016 #define exp10(x)     exp(LN_10*(x))
00017 
00018 double modhist[HISTRES];           /* modified histogram */
00019 double mhistot;                    /* modified histogram total */
00020 double cumf[HISTRES+1];            /* cumulative distribution function */
00021 
00022 static double centprob(int  x, int y);
00023 static void mkcumf(void);
00024 static double cf(double     b);
00025 static double BLw(double    Lw);
00026 #if ADJ_VEIL
00027 static void mkcrfimage(void);
00028 #endif
00029 
00030 
00031 
00032 extern void
00033 getfixations(        /* load fixation history list */
00034 FILE   *fp
00035 )
00036 {
00037 #define       FIXHUNK              128
00038        RESOLU fvres;
00039        int    pos[2];
00040        register int  px, py, i;
00041                             /* initialize our resolution struct */
00042        if ((fvres.rt=inpres.rt)&YMAJOR) {
00043               fvres.xr = fvxr;
00044               fvres.yr = fvyr;
00045        } else {
00046               fvres.xr = fvyr;
00047               fvres.yr = fvxr;
00048        }
00049                             /* read each picture position */
00050        while (fscanf(fp, "%d %d", &pos[0], &pos[1]) == 2) {
00051                             /* convert to closest index in foveal image */
00052               loc2pix(pos, &fvres,
00053                             (pos[0]+.5)/inpres.xr, (pos[1]+.5)/inpres.yr);
00054                             /* include nine neighborhood samples */
00055               for (px = pos[0]-1; px <= pos[0]+1; px++) {
00056                      if (px < 0 || px >= fvxr)
00057                             continue;
00058                      for (py = pos[1]-1; py <= pos[1]+1; py++) {
00059                             if (py < 0 || py >= fvyr)
00060                                    continue;
00061                             for (i = nfixations; i-- > 0; )
00062                                    if (fixlst[i][0] == px &&
00063                                                  fixlst[i][1] == py)
00064                                           break;
00065                             if (i >= 0)
00066                                    continue;     /* already there */
00067                             if (nfixations % FIXHUNK == 0) {
00068                                    if (nfixations)
00069                                           fixlst = (short (*)[2])
00070                                                  realloc((void *)fixlst,
00071                                                  (nfixations+FIXHUNK)*
00072                                                  2*sizeof(short));
00073                                    else
00074                                           fixlst = (short (*)[2])malloc(
00075                                                  FIXHUNK*2*sizeof(short)
00076                                                  );
00077                                    if (fixlst == NULL)
00078                                           syserror("malloc");
00079                             }
00080                             fixlst[nfixations][0] = px;
00081                             fixlst[nfixations][1] = py;
00082                             nfixations++;
00083                      }
00084               }
00085        }
00086        if (!feof(fp)) {
00087               fprintf(stderr, "%s: format error reading fixation data\n",
00088                             progname);
00089               exit(1);
00090        }
00091 #undef FIXHUNK
00092 }
00093 
00094 
00095 extern void
00096 gethisto(                   /* load precomputed luminance histogram */
00097        FILE   *fp
00098 )
00099 {
00100        double histo[MAXPREHIST];
00101        double histart, histep;
00102        double b, lastb, w;
00103        int    n;
00104        register int  i;
00105                                    /* load data */
00106        for (i = 0; i < MAXPREHIST &&
00107                      fscanf(fp, "%lf %lf", &b, &histo[i]) == 2; i++) {
00108               if (i > 1 && fabs(b - lastb - histep) > .001) {
00109                      fprintf(stderr,
00110                             "%s: uneven step size in histogram data\n",
00111                                    progname);
00112                      exit(1);
00113               }
00114               if (i == 1)
00115                      if ((histep = b - (histart = lastb)) <= FTINY) {
00116                             fprintf(stderr,
00117                                    "%s: illegal step in histogram data\n",
00118                                           progname);
00119                             exit(1);
00120                      }
00121               lastb = b;
00122        }
00123        if (i < 2 || !feof(fp)) {
00124               fprintf(stderr,
00125               "%s: format/length error loading histogram (log10L %f at %d)\n",
00126                             progname, b, i);
00127               exit(1);
00128        }
00129        n = i;
00130        histart *= LN_10;
00131        histep *= LN_10;
00132                                    /* find extrema */
00133        for (i = 0; i < n && histo[i] <= FTINY; i++)
00134               ;
00135        bwmin = histart + (i-.001)*histep;
00136        for (i = n; i-- && histo[i] <= FTINY; )
00137               ;
00138        bwmax = histart + (i+1.001)*histep;
00139        if (bwmax > Bl(LMAX))
00140               bwmax = Bl(LMAX);
00141        if (bwmin < Bl(LMIN))
00142               bwmin = Bl(LMIN);
00143        else                        /* duplicate bottom bin */
00144               bwmin = bwmax - (bwmax-bwmin)*HISTRES/(HISTRES-1);
00145                                    /* convert histogram */
00146        bwavg = 0.; histot = 0.;
00147        for (i = 0; i < HISTRES; i++)
00148               bwhist[i] = 0.;
00149        for (i = 0, b = histart; i < n; i++, b += histep) {
00150               if (b < bwmin+FTINY) continue;
00151               if (b >= bwmax-FTINY) break;
00152               w = histo[i];
00153               bwavg += w*b;
00154               bwhist[bwhi(b)] += w;
00155               histot += w;
00156        }
00157        bwavg /= histot;
00158        if (bwmin > Bl(LMIN)+FTINY) {      /* add false samples at bottom */
00159               bwhist[1] *= 0.5;
00160               bwhist[0] += bwhist[1];
00161        }
00162 }
00163 
00164 
00165 static double
00166 centprob(                   /* center-weighting probability function */
00167        int    x,
00168        int    y
00169 )
00170 {
00171        double xr, yr, p;
00172                             /* paraboloid, 0 at 90 degrees from center */
00173        xr = (x - .5*(fvxr-1))/90.; /* 180 degree fisheye has fv?r == 90 */
00174        yr = (y - .5*(fvyr-1))/90.;
00175        p = 1. - xr*xr - yr*yr;
00176        return(p < 0. ? 0. : p);
00177 }
00178 
00179 
00180 extern void
00181 comphist(void)                     /* create foveal sampling histogram */
00182 {
00183        double l, b, w, lwmin, lwmax;
00184        register int  x, y;
00185                                    /* check for precalculated histogram */
00186        if (what2do&DO_PREHIST)
00187               return;
00188        lwmin = 1e10;               /* find extrema */
00189        lwmax = 0.;
00190        for (y = 0; y < fvyr; y++)
00191               for (x = 0; x < fvxr; x++) {
00192                      l = plum(fovscan(y)[x]);
00193                      if (l < lwmin) lwmin = l;
00194                      if (l > lwmax) lwmax = l;
00195               }
00196        lwmax *= 1.01;
00197        if (lwmax > LMAX)
00198               lwmax = LMAX;
00199        bwmax = Bl(lwmax);
00200        if (lwmin < LMIN) {
00201               lwmin = LMIN;
00202               bwmin = Bl(LMIN);
00203        } else {                    /* duplicate bottom bin */
00204               bwmin = bwmax - (bwmax-Bl(lwmin))*HISTRES/(HISTRES-1);
00205               lwmin = Lb(bwmin);
00206        }
00207                                    /* (re)compute histogram */
00208        bwavg = 0.;
00209        histot = 0.;
00210        for (x = 0; x < HISTRES; x++)
00211               bwhist[x] = 0.;
00212                                    /* global average */
00213        if (!(what2do&DO_FIXHIST) || fixfrac < 1.-FTINY)
00214               for (y = 0; y < fvyr; y++)
00215                      for (x = 0; x < fvxr; x++) {
00216                             l = plum(fovscan(y)[x]);
00217                             if (l < lwmin+FTINY) continue;
00218                             if (l >= lwmax-FTINY) continue;
00219                             b = Bl(l);
00220                             w = what2do&DO_CWEIGHT ? centprob(x,y) : 1.;
00221                             bwavg += w*b;
00222                             bwhist[bwhi(b)] += w;
00223                             histot += w;
00224                      }
00225                                    /* average fixation points */
00226        if (what2do&DO_FIXHIST && nfixations > 0) {
00227               if (histot > FTINY)
00228                      w = fixfrac/(1.-fixfrac)*histot/nfixations;
00229               else
00230                      w = 1.;
00231               for (x = 0; x < nfixations; x++) {
00232                      l = plum(fovscan(fixlst[x][1])[fixlst[x][0]]);
00233                      if (l < lwmin+FTINY) continue;
00234                      if (l >= lwmax-FTINY) continue;
00235                      b = Bl(l);
00236                      bwavg += w*b;
00237                      bwhist[bwhi(b)] += w;
00238                      histot += w;
00239               }
00240        }
00241        bwavg /= histot;
00242        if (lwmin > LMIN+FTINY) {   /* add false samples at bottom */
00243               bwhist[1] *= 0.5;
00244               bwhist[0] += bwhist[1];
00245        }
00246 }
00247 
00248 
00249 static void
00250 mkcumf(void)                /* make cumulative distribution function */
00251 {
00252        register int  i;
00253        register double      sum;
00254 
00255        mhistot = 0.;        /* compute modified total */
00256        for (i = 0; i < HISTRES; i++)
00257               mhistot += modhist[i];
00258 
00259        sum = 0.;            /* compute cumulative function */
00260        for (i = 0; i < HISTRES; i++) {
00261               cumf[i] = sum/mhistot;
00262               sum += modhist[i];
00263        }
00264        cumf[HISTRES] = 1.;
00265 }
00266 
00267 
00268 static double
00269 cf(                         /* return cumulative function at b */
00270        double b
00271 )
00272 {
00273        double x;
00274        register int  i;
00275 
00276        i = x = HISTRES*(b - bwmin)/(bwmax - bwmin);
00277        x -= (double)i;
00278        return(cumf[i]*(1.-x) + cumf[i+1]*x);
00279 }
00280 
00281 
00282 static double
00283 BLw(                        /* map world luminance to display brightness */
00284        double Lw
00285 )
00286 {
00287        double b;
00288 
00289        if (Lw <= LMIN || (b = Bl(Lw)) <= bwmin+FTINY)
00290               return(Bldmin);
00291        if (b >= bwmax-FTINY)
00292               return(Bldmax);
00293        return(Bldmin + cf(b)*(Bldmax-Bldmin));
00294 }
00295 
00296 
00297 extern double
00298 htcontrs(            /* human threshold contrast sensitivity, dL(La) */
00299        double La
00300 )
00301 {
00302        double l10La, l10dL;
00303                             /* formula taken from Ferwerda et al. [SG96] */
00304        if (La < 1.148e-4)
00305               return(1.38e-3);
00306        l10La = log10(La);
00307        if (l10La < -1.44)          /* rod response regime */
00308               l10dL = pow(.405*l10La + 1.6, 2.18) - 2.86;
00309        else if (l10La < -.0184)
00310               l10dL = l10La - .395;
00311        else if (l10La < 1.9)              /* cone response regime */
00312               l10dL = pow(.249*l10La + .65, 2.7) - .72;
00313        else
00314               l10dL = l10La - 1.255;
00315 
00316        return(exp10(l10dL));
00317 }
00318 
00319 
00320 extern double
00321 clampf(                     /* histogram clamping function */
00322        double Lw
00323 )
00324 {
00325        double bLw, ratio;
00326 
00327        bLw = BLw(Lw);              /* apply current brightness map */
00328        ratio = what2do&DO_HSENS ? htcontrs(Lb(bLw))/htcontrs(Lw) : Lb(bLw)/Lw;
00329        return(ratio/(Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)));
00330 }
00331 
00332 extern double
00333 crfactor(                   /* contrast reduction factor */
00334        double Lw
00335 )
00336 {
00337        int    i = HISTRES*(Bl(Lw) - bwmin)/(bwmax - bwmin);
00338        double bLw, ratio, Tdb;
00339 
00340        if (i <= 0)
00341               return(1.0);
00342        if (i >= HISTRES)
00343               return(1.0);
00344        bLw = BLw(Lw);
00345        ratio = what2do&DO_HSENS ? htcontrs(Lb(bLw))/htcontrs(Lw) : Lb(bLw)/Lw;
00346        Tdb = mhistot * (bwmax - bwmin) / HISTRES;
00347        return(modhist[i]*Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)/(Tdb*ratio));
00348 }
00349 
00350 
00351 #if ADJ_VEIL
00352 static void
00353 mkcrfimage(void)                   /* compute contrast reduction factor image */
00354 {
00355        int    i;
00356        float  *crfptr;
00357        COLOR  *fovptr;
00358 
00359        if (crfimg == NULL)
00360               crfimg = (float *)malloc(fvxr*fvyr*sizeof(float));
00361        if (crfimg == NULL)
00362               syserror("malloc");
00363        crfptr = crfimg;
00364        fovptr = fovimg;
00365        for (i = fvxr*fvyr; i--; crfptr++, fovptr++)
00366               crfptr[0] = crfactor(plum(fovptr[0]));
00367 }
00368 #endif
00369 
00370 
00371 extern int
00372 mkbrmap(void)               /* make dynamic range map */
00373 {
00374        double Tdb, b, s;
00375        double ceiling, trimmings;
00376        register int  i;
00377                                    /* copy initial histogram */
00378        memcpy((void *)modhist, (void *)bwhist, sizeof(modhist));
00379        s = (bwmax - bwmin)/HISTRES;       /* s is delta b */
00380                                    /* loop until satisfactory */
00381        do {
00382               mkcumf();                   /* sync brightness mapping */
00383               if (mhistot <= histot*CVRATIO)
00384                      return(-1);          /* no compression needed! */
00385               Tdb = mhistot * s;
00386               trimmings = 0.;                    /* clip to envelope */
00387               for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
00388                      ceiling = Tdb*clampf(Lb(b));
00389                      if (modhist[i] > ceiling) {
00390                             trimmings += modhist[i] - ceiling;
00391                             modhist[i] = ceiling;
00392                      }
00393               }
00394        } while (trimmings > histot*CVRATIO);
00395 
00396 #if ADJ_VEIL
00397        mkcrfimage();               /* contrast reduction image */
00398 #endif
00399 
00400        return(0);                  /* we got it */
00401 }
00402 
00403 
00404 extern void
00405 scotscan(            /* apply scotopic color sensitivity loss */
00406        COLOR  *scan,
00407        int    xres
00408 )
00409 {
00410        COLOR  ctmp;
00411        double incolor, b, Lw;
00412        register int  i;
00413 
00414        for (i = 0; i < xres; i++) {
00415               Lw = plum(scan[i]);
00416               if (Lw >= TopMesopic)
00417                      incolor = 1.;
00418               else if (Lw <= BotMesopic)
00419                      incolor = 0.;
00420               else
00421                      incolor = (Lw - BotMesopic) /
00422                                    (TopMesopic - BotMesopic);
00423               if (incolor < 1.-FTINY) {
00424                      b = (1.-incolor)*slum(scan[i])*inpexp/SWNORM;
00425                      if (lumf == rgblum) b /= WHTEFFICACY;
00426                      setcolor(ctmp, b, b, b);
00427                      if (incolor <= FTINY)
00428                             setcolor(scan[i], 0., 0., 0.);
00429                      else
00430                             scalecolor(scan[i], incolor);
00431                      addcolor(scan[i], ctmp);
00432               }
00433        }
00434 }
00435 
00436 
00437 extern void
00438 mapscan(             /* apply tone mapping operator to scanline */
00439        COLOR  *scan,
00440        int    xres
00441 )
00442 {
00443        double mult, Lw, b;
00444        register int  x;
00445 
00446        for (x = 0; x < xres; x++) {
00447               Lw = plum(scan[x]);
00448               if (Lw < LMIN) {
00449                      setcolor(scan[x], 0., 0., 0.);
00450                      continue;
00451               }
00452               b = BLw(Lw);         /* apply brightness mapping */
00453               mult = (Lb(b) - ldmin)/(ldmax - ldmin)/(Lw*inpexp);
00454               if (lumf == rgblum) mult *= WHTEFFICACY;
00455               scalecolor(scan[x], mult);
00456        }
00457 }
00458 
00459 
00460 extern void
00461 putmapping(                 /* put out mapping function */
00462        FILE   *fp
00463 )
00464 {
00465        double b, s;
00466        register int  i;
00467        double wlum, sf, dlum;
00468 
00469        sf = scalef*inpexp;
00470        if (lumf == cielum) sf *= WHTEFFICACY;
00471        s = (bwmax - bwmin)/HISTRES;
00472        for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
00473               wlum = Lb(b);
00474               if (what2do&DO_LINEAR) {
00475                      dlum = sf*wlum;
00476                      if (dlum > ldmax) dlum = ldmax;
00477                      else if (dlum < ldmin) dlum = ldmin;
00478                      fprintf(fp, "%e %e\n", wlum, dlum);
00479               } else
00480                      fprintf(fp, "%e %e\n", wlum, Lb(BLw(wlum)));
00481        }
00482 }