Back to index

radiance  4R0+20100331
pcond4.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: pcond4.c,v 3.19 2004/11/14 16:57:18 greg Exp $";
00003 #endif
00004 /*
00005  * Routines for veiling glare and loss of acuity.
00006  */
00007 
00008 #include "pcond.h"
00009 
00010 /************** VEILING STUFF *****************/
00011 
00012 #define       VADAPT        0.08          /* fraction of adaptation from veil */
00013 
00014 static COLOR  *veilimg = NULL;     /* veiling image */
00015 
00016 #define veilscan(y)  (veilimg+(y)*fvxr)
00017 
00018 static float  (*raydir)[3] = NULL; /* ray direction for each pixel */
00019 
00020 #define rdirscan(y)  (raydir+(y)*fvxr)
00021 
00022 static void compraydir(void);
00023 #if ADJ_VEIL
00024 static void smoothveil(void);
00025 #endif
00026 
00027 
00028 static void
00029 compraydir(void)                          /* compute ray directions */
00030 {
00031        FVECT  rorg, rdir;
00032        double h, v;
00033        register int  x, y;
00034 
00035        if (raydir != NULL)         /* already done? */
00036               return;
00037        raydir = (float (*)[3])malloc(fvxr*fvyr*3*sizeof(float));
00038        if (raydir == NULL)
00039               syserror("malloc");
00040 
00041        for (y = 0; y < fvyr; y++) {
00042               switch (inpres.rt) {
00043               case YMAJOR: case YMAJOR|XDECR:
00044                      v = (y+.5)/fvyr; break;
00045               case YMAJOR|YDECR: case YMAJOR|YDECR|XDECR:
00046                      v = 1. - (y+.5)/fvyr; break;
00047               case 0: case YDECR:
00048                      h = (y+.5)/fvyr; break;
00049               case XDECR: case XDECR|YDECR:
00050                      h = 1. - (y+.5)/fvyr; break;
00051               }
00052               for (x = 0; x < fvxr; x++) {
00053                      switch (inpres.rt) {
00054                      case YMAJOR: case YMAJOR|YDECR:
00055                             h = (x+.5)/fvxr; break;
00056                      case YMAJOR|XDECR: case YMAJOR|XDECR|YDECR:
00057                             h = 1. - (x+.5)/fvxr; break;
00058                      case 0: case XDECR:
00059                             v = (x+.5)/fvxr; break;
00060                      case YDECR: case YDECR|XDECR:
00061                             v = 1. - (x+.5)/fvxr; break;
00062                      }
00063                      if (viewray(rorg, rdir, &ourview, h, v)
00064                                    >= -FTINY) {
00065                             rdirscan(y)[x][0] = rdir[0];
00066                             rdirscan(y)[x][1] = rdir[1];
00067                             rdirscan(y)[x][2] = rdir[2];
00068                      } else {
00069                             rdirscan(y)[x][0] =
00070                             rdirscan(y)[x][1] =
00071                             rdirscan(y)[x][2] = 0.0;
00072                      }
00073               }
00074        }
00075 }
00076 
00077 
00078 extern void
00079 compveil(void)                            /* compute veiling image */
00080 {
00081        double t2, t2sum;
00082        COLOR  ctmp, vsum;
00083        int    px, py;
00084        register int  x, y;
00085 
00086        if (veilimg != NULL)        /* already done? */
00087               return;
00088                                    /* compute ray directions */
00089        compraydir();
00090                                    /* compute veil image */
00091        veilimg = (COLOR *)malloc(fvxr*fvyr*sizeof(COLOR));
00092        if (veilimg == NULL)
00093               syserror("malloc");
00094        for (py = 0; py < fvyr; py++)
00095               for (px = 0; px < fvxr; px++) {
00096                      t2sum = 0.;
00097                      setcolor(vsum, 0., 0., 0.);
00098                      for (y = 0; y < fvyr; y++)
00099                             for (x = 0; x < fvxr; x++) {
00100                                    if (x == px && y == py) continue;
00101                                    t2 = DOT(rdirscan(py)[px],
00102                                                  rdirscan(y)[x]);
00103                                    if (t2 <= FTINY) continue;
00104                                    /*     use approximation instead
00105                                    t3 = acos(t2);
00106                                    t2 = t2/(t3*t3);
00107                                    */
00108                                    t2 *= .5 / (1. - t2);
00109                                    copycolor(ctmp, fovscan(y)[x]);
00110                                    scalecolor(ctmp, t2);
00111                                    addcolor(vsum, ctmp);
00112                                    t2sum += t2;
00113                             }
00114                      /* VADAPT of original is subtracted in addveil() */
00115                      if (t2sum > FTINY)
00116                             scalecolor(vsum, VADAPT/t2sum);
00117                      copycolor(veilscan(py)[px], vsum);
00118               }
00119                                    /* modify FOV sample image */
00120        for (y = 0; y < fvyr; y++)
00121               for (x = 0; x < fvxr; x++) {
00122                      scalecolor(fovscan(y)[x], 1.-VADAPT);
00123                      addcolor(fovscan(y)[x], veilscan(y)[x]);
00124               }
00125        comphist();                 /* recompute histogram */
00126 }
00127 
00128 
00129 #if ADJ_VEIL
00130 /*
00131  * The following veil adjustment was added to compensate for
00132  * the fact that contrast reduction gets confused with veil
00133  * in the human visual system.  Therefore, we reduce the
00134  * veil in portions of the image where our mapping has
00135  * already reduced contrast below the target value.
00136  * This gets called after the intial veil has been computed
00137  * and added to the foveal image, and the mapping has been
00138  * determined.
00139  */
00140 extern void
00141 adjveil(void)                      /* adjust veil image */
00142 {
00143        float  *crfptr = crfimg;
00144        COLOR  *fovptr = fovimg;
00145        COLOR  *veilptr = veilimg;
00146        double s2nits = 1./inpexp;
00147        double vl, vl2, fovl, vlsum;
00148        double deltavc[3];
00149        int    i, j;
00150 
00151        if (lumf == rgblum)
00152               s2nits *= WHTEFFICACY;
00153 
00154        for (i = fvxr*fvyr; i--; crfptr++, fovptr++, veilptr++) {
00155               if (crfptr[0] >= 0.95)
00156                      continue;
00157               vl = plum(veilptr[0]);
00158               fovl = (plum(fovptr[0]) - vl) * (1./(1.-VADAPT));
00159               if (vl <= 0.05*fovl)
00160                      continue;
00161               vlsum = vl;
00162               for (j = 2; j < 11; j++) {
00163                      vlsum += crfptr[0]*vl - (1.0 - crfptr[0])*fovl;
00164                      vl2 = vlsum / (double)j;
00165                      if (vl2 < 0.0)
00166                             vl2 = 0.0;
00167                      crfptr[0] = crfactor(fovl + vl2);
00168               }
00169               /* desaturation code causes color fringes at this level */
00170               for (j = 3; j--; ) {
00171                      double vc = colval(veilptr[0],j);
00172                      double fovc = (colval(fovptr[0],j) - vc) *
00173                                           (1./(1.-VADAPT));
00174                      deltavc[j] = (1.-crfptr[0])*(fovl/s2nits - fovc);
00175                      if (vc + deltavc[j] < 0.0)
00176                             break;
00177               }
00178               if (j < 0)
00179                      addcolor(veilptr[0], deltavc);
00180               else
00181                      scalecolor(veilptr[0], vl2/vl);
00182        }
00183        smoothveil();               /* smooth our result */
00184 }
00185 
00186 
00187 static void
00188 smoothveil(void)                          /* smooth veil image */
00189 {
00190        COLOR  *nveilimg;
00191        COLOR  *ovptr, *nvptr;
00192        int    x, y, i;
00193 
00194        nveilimg = (COLOR *)malloc(fvxr*fvyr*sizeof(COLOR));
00195        if (nveilimg == NULL)
00196               return;
00197        for (y = 1; y < fvyr-1; y++) {
00198               ovptr = veilimg + y*fvxr + 1;
00199               nvptr = nveilimg + y*fvxr + 1;
00200               for (x = 1; x < fvxr-1; x++, ovptr++, nvptr++)
00201                      for (i = 3; i--; )
00202                             nvptr[0][i] = 0.5 * ovptr[0][i]
00203                                    + (1./12.) *
00204                             (ovptr[-1][i] + ovptr[-fvxr][i] +
00205                                    ovptr[1][i] + ovptr[fvxr][i])
00206                                    + (1./24.) *
00207                             (ovptr[-fvxr-1][i] + ovptr[-fvxr+1][i] +
00208                                    ovptr[fvxr-1][i] + ovptr[fvxr+1][i]);
00209        }
00210        ovptr = veilimg + 1;
00211        nvptr = nveilimg + 1;
00212        for (x = 1; x < fvxr-1; x++, ovptr++, nvptr++)
00213               for (i = 3; i--; )
00214                      nvptr[0][i] = 0.5 * ovptr[0][i]
00215                             + (1./9.) *
00216                      (ovptr[-1][i] + ovptr[1][i] + ovptr[fvxr][i])
00217                             + (1./12.) *
00218                      (ovptr[fvxr-1][i] + ovptr[fvxr+1][i]);
00219        ovptr = veilimg + (fvyr-1)*fvxr + 1;
00220        nvptr = nveilimg + (fvyr-1)*fvxr + 1;
00221        for (x = 1; x < fvxr-1; x++, ovptr++, nvptr++)
00222               for (i = 3; i--; )
00223                      nvptr[0][i] = 0.5 * ovptr[0][i]
00224                             + (1./9.) *
00225                      (ovptr[-1][i] + ovptr[1][i] + ovptr[-fvxr][i])
00226                             + (1./12.) *
00227                      (ovptr[-fvxr-1][i] + ovptr[-fvxr+1][i]);
00228        ovptr = veilimg + fvxr;
00229        nvptr = nveilimg + fvxr;
00230        for (y = 1; y < fvyr-1; y++, ovptr += fvxr, nvptr += fvxr)
00231               for (i = 3; i--; )
00232                      nvptr[0][i] = 0.5 * ovptr[0][i]
00233                             + (1./9.) *
00234                      (ovptr[-fvxr][i] + ovptr[1][i] + ovptr[fvxr][i])
00235                             + (1./12.) *
00236                      (ovptr[-fvxr+1][i] + ovptr[fvxr+1][i]);
00237        ovptr = veilimg + fvxr - 1;
00238        nvptr = nveilimg + fvxr - 1;
00239        for (y = 1; y < fvyr-1; y++, ovptr += fvxr, nvptr += fvxr)
00240               for (i = 3; i--; )
00241                      nvptr[0][i] = 0.5 * ovptr[0][i]
00242                             + (1./9.) *
00243                      (ovptr[-fvxr][i] + ovptr[-1][i] + ovptr[fvxr][i])
00244                             + (1./12.) *
00245                      (ovptr[-fvxr-1][i] + ovptr[fvxr-1][i]);
00246        for (i = 3; i--; ) {
00247               nveilimg[0][i] = veilimg[0][i];
00248               nveilimg[fvxr-1][i] = veilimg[fvxr-1][i];
00249               nveilimg[(fvyr-1)*fvxr][i] = veilimg[(fvyr-1)*fvxr][i];
00250               nveilimg[fvyr*fvxr-1][i] = veilimg[fvyr*fvxr-1][i];
00251        }
00252        free((void *)veilimg);
00253        veilimg = nveilimg;
00254 }
00255 #endif
00256 
00257 extern void
00258 addveil(                           /* add veil to scanline */
00259        COLOR  *sl,
00260        int    y
00261 )
00262 {
00263        int    vx, vy;
00264        double dx, dy;
00265        double lv, uv;
00266        register int  x, i;
00267 
00268        vy = dy = (y+.5)/numscans(&inpres)*fvyr - .5;
00269        while (vy >= fvyr-1) vy--;
00270        dy -= (double)vy;
00271        for (x = 0; x < scanlen(&inpres); x++) {
00272               vx = dx = (x+.5)/scanlen(&inpres)*fvxr - .5;
00273               while (vx >= fvxr-1) vx--;
00274               dx -= (double)vx;
00275               for (i = 0; i < 3; i++) {
00276                      lv = (1.-dy)*colval(veilscan(vy)[vx],i) +
00277                                    dy*colval(veilscan(vy+1)[vx],i);
00278                      uv = (1.-dy)*colval(veilscan(vy)[vx+1],i) +
00279                                    dy*colval(veilscan(vy+1)[vx+1],i);
00280                      colval(sl[x],i) = (1.-VADAPT)*colval(sl[x],i) +
00281                                    (1.-dx)*lv + dx*uv;
00282               }
00283        }
00284 }
00285 
00286 
00287 /****************** ACUITY STUFF *******************/
00288 
00289 typedef struct {
00290        short  sampe;        /* sample area size (exponent of 2) */
00291        short  nscans;              /* number of scanlines in this bar */
00292        int    len;          /* individual scanline length */
00293        int    nread;        /* number of scanlines loaded */
00294        COLOR  *sdata;              /* scanbar data */
00295 } SCANBAR;
00296 
00297 #define bscan(sb,y)  ((COLOR *)(sb)->sdata+((y)%(sb)->nscans)*(sb)->len)
00298 
00299 SCANBAR       *rootbar;            /* root scan bar (lowest resolution) */
00300 
00301 float  *inpacuD = NULL;     /* input acuity data (cycles/degree) */
00302 
00303 #define tsampr(x,y)  inpacuD[(y)*fvxr+(x)]
00304 
00305 static COLOR * getascan(SCANBAR    *sb, int      y);
00306 static void acusample(COLOR col, int      x, int y, double     sr);
00307 static void ascanval(COLOR  col, int      x, int y, SCANBAR    *sb);
00308 static SCANBAR       *sballoc(int  se, int       ns, int       sl);
00309 
00310 extern double
00311 hacuity(                    /* return visual acuity in cycles/degree */
00312        double La
00313 )
00314 {
00315                                    /* functional fit */
00316        return(17.25*atan(1.4*log10(La) + 0.35) + 25.72);
00317 }
00318 
00319 
00320 static COLOR *
00321 getascan(                   /* find/read scanline y for scanbar sb */
00322        register SCANBAR     *sb,
00323        int    y
00324 )
00325 {
00326        register COLOR       *sl0, *sl1, *mysl;
00327        register int  i;
00328 
00329        if (y < sb->nread - sb->nscans)                  /* too far back? */
00330               return(NULL);
00331        for ( ; y >= sb->nread; sb->nread++) {           /* read as necessary */
00332               mysl = bscan(sb, sb->nread);
00333               if (sb->sampe == 0) {
00334                      if (freadscan(mysl, sb->len, infp) < 0) {
00335                             fprintf(stderr, "%s: %s: scanline read error\n",
00336                                           progname, infn);
00337                             exit(1);
00338                      }
00339               } else {
00340                      sl0 = getascan(sb+1, 2*y);
00341                      if (sl0 == NULL)
00342                             return(NULL);
00343                      sl1 = getascan(sb+1, 2*y+1);
00344                      for (i = 0; i < sb->len; i++) {
00345                             copycolor(mysl[i], sl0[2*i]);
00346                             addcolor(mysl[i], sl0[2*i+1]);
00347                             addcolor(mysl[i], sl1[2*i]);
00348                             addcolor(mysl[i], sl1[2*i+1]);
00349                             scalecolor(mysl[i], 0.25);
00350                      }
00351               }
00352        }
00353        return(bscan(sb, y));
00354 }
00355 
00356 
00357 extern void
00358 acuscan(             /* get acuity-sampled scanline */
00359        COLOR  *scln,
00360        int    y
00361 )
00362 {
00363        double sr;
00364        double dx, dy;
00365        int    ix, iy;
00366        register int  x;
00367        
00368        if (inpacuD == NULL)
00369               return;
00370                                    /* compute foveal y position */
00371        iy = dy = (y+.5)/numscans(&inpres)*fvyr - .5;
00372        while (iy >= fvyr-1) iy--;
00373        dy -= (double)iy;
00374        for (x = 0; x < scanlen(&inpres); x++) {
00375                                    /* compute foveal x position */
00376               ix = dx = (x+.5)/scanlen(&inpres)*fvxr - .5;
00377               while (ix >= fvxr-1) ix--;
00378               dx -= (double)ix;
00379                                    /* interpolate sample rate */
00380               sr = (1.-dy)*((1.-dx)*tsampr(ix,iy) + dx*tsampr(ix+1,iy)) +
00381                      dy*((1.-dx)*tsampr(ix,iy+1) + dx*tsampr(ix+1,iy+1));
00382 
00383               acusample(scln[x], x, y, sr);      /* compute sample */
00384        }
00385 }
00386 
00387 
00388 static void
00389 acusample(    /* interpolate sample at (x,y) using rate sr */
00390        COLOR  col,
00391        int    x,
00392        int    y,
00393        double sr
00394 )
00395 {
00396        COLOR  c1;
00397        double d;
00398        register SCANBAR     *sb0;
00399 
00400        for (sb0 = rootbar; sb0->sampe != 0 && 1<<sb0[1].sampe > sr; sb0++)
00401               ;
00402        ascanval(col, x, y, sb0);
00403        if (sb0->sampe == 0)        /* don't extrapolate highest */
00404               return;
00405        ascanval(c1, x, y, sb0+1);
00406        d = ((1<<sb0->sampe) - sr)/(1<<sb0[1].sampe);
00407        scalecolor(col, 1.-d);
00408        scalecolor(c1, d);
00409        addcolor(col, c1);
00410 }
00411 
00412 
00413 static void
00414 ascanval(            /* interpolate scanbar at orig. coords (x,y) */
00415        COLOR  col,
00416        int    x,
00417        int    y,
00418        SCANBAR       *sb
00419 )
00420 {
00421        COLOR  *sl0, *sl1, c1, c1y;
00422        double dx, dy;
00423        int    ix, iy;
00424 
00425        if (sb->sampe == 0) {              /* no need to interpolate */
00426               sl0 = getascan(sb, y);
00427               copycolor(col, sl0[x]);
00428               return;
00429        }
00430                                    /* compute coordinates for sb */
00431        ix = dx = (x+.5)/(1<<sb->sampe) - .5;
00432        while (ix >= sb->len-1) ix--;
00433        dx -= (double)ix;
00434        iy = dy = (y+.5)/(1<<sb->sampe) - .5;
00435        while (iy >= (numscans(&inpres)>>sb->sampe)-1) iy--;
00436        dy -= (double)iy;
00437                                    /* get scanlines */
00438        sl0 = getascan(sb, iy);
00439 #ifdef DEBUG
00440        if (sl0 == NULL)
00441               error(INTERNAL, "cannot backspace in ascanval");
00442 #endif
00443        sl1 = getascan(sb, iy+1);
00444                                    /* 2D linear interpolation */
00445        copycolor(col, sl0[ix]);
00446        scalecolor(col, 1.-dx);
00447        copycolor(c1, sl0[ix+1]);
00448        scalecolor(c1, dx);
00449        addcolor(col, c1);
00450        copycolor(c1y, sl1[ix]);
00451        scalecolor(c1y, 1.-dx);
00452        copycolor(c1, sl1[ix+1]);
00453        scalecolor(c1, dx);
00454        addcolor(c1y, c1);
00455        scalecolor(col, 1.-dy);
00456        scalecolor(c1y, dy);
00457        addcolor(col, c1y);
00458        for (ix = 0; ix < 3; ix++)  /* make sure no negative */
00459               if (colval(col,ix) < 0.)
00460                      colval(col,ix) = 0.;
00461 }
00462 
00463 
00464 static SCANBAR       *
00465 sballoc(             /* allocate scanbar */
00466        int    se,           /* sampling rate exponent */
00467        int    ns,           /* number of scanlines */
00468        int    sl            /* original scanline length */
00469 )
00470 {
00471        SCANBAR       *sbarr;
00472        register SCANBAR     *sb;
00473 
00474        sbarr = sb = (SCANBAR *)malloc((se+1)*sizeof(SCANBAR));
00475        if (sb == NULL)
00476               syserror("malloc");
00477        do {
00478               sb->len = sl>>se;
00479               if (sb->len <= 0)
00480                      continue;
00481               sb->sampe = se;
00482               sb->nscans = ns;
00483               sb->sdata = (COLOR *)malloc(sb->len*ns*sizeof(COLOR));
00484               if (sb->sdata == NULL)
00485                      syserror("malloc");
00486               sb->nread = 0;
00487               ns <<= 1;
00488               sb++;
00489        } while (--se >= 0);
00490        return(sbarr);
00491 }
00492 
00493 
00494 extern void
00495 initacuity(void)                   /* initialize variable acuity sampling */
00496 {
00497        FVECT  diffx, diffy, cp;
00498        double omega, maxsr;
00499        register int  x, y, i;
00500        
00501        compraydir();               /* compute ray directions */
00502 
00503        inpacuD = (float *)malloc(fvxr*fvyr*sizeof(float));
00504        if (inpacuD == NULL)
00505               syserror("malloc");
00506        maxsr = 1.;                 /* compute internal sample rates */
00507        for (y = 1; y < fvyr-1; y++)
00508               for (x = 1; x < fvxr-1; x++) {
00509                      for (i = 0; i < 3; i++) {
00510                             diffx[i] = 0.5*fvxr/scanlen(&inpres) *
00511                                           (rdirscan(y)[x+1][i] -
00512                                           rdirscan(y)[x-1][i]);
00513                             diffy[i] = 0.5*fvyr/numscans(&inpres) *
00514                                           (rdirscan(y+1)[x][i] -
00515                                           rdirscan(y-1)[x][i]);
00516                      }
00517                      fcross(cp, diffx, diffy);
00518                      omega = 0.5 * sqrt(DOT(cp,cp));
00519                      if (omega <= FTINY*FTINY)
00520                             tsampr(x,y) = 1.;
00521                      else if ((tsampr(x,y) = PI/180. / sqrt(omega) /
00522                                    hacuity(plum(fovscan(y)[x]))) > maxsr)
00523                             maxsr = tsampr(x,y);
00524               }
00525                                    /* copy perimeter (easier) */
00526        for (x = 1; x < fvxr-1; x++) {
00527               tsampr(x,0) = tsampr(x,1);
00528               tsampr(x,fvyr-1) = tsampr(x,fvyr-2);
00529        }
00530        for (y = 0; y < fvyr; y++) {
00531               tsampr(0,y) = tsampr(1,y);
00532               tsampr(fvxr-1,y) = tsampr(fvxr-2,y);
00533        }
00534                                    /* initialize with next power of two */
00535        rootbar = sballoc((int)(log(maxsr)/log(2.))+1, 2, scanlen(&inpres));
00536 }