Back to index

radiance  4R0+20100331
rhpict2.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: rhpict2.c,v 3.15 2004/01/01 11:21:55 schorsch Exp $";
00003 #endif
00004 /*
00005  * Rendering routines for rhpict.
00006  */
00007 
00008 #include <string.h>
00009 
00010 #include "holo.h"
00011 #include "view.h"
00012 #include "random.h"
00013 
00014 #ifndef DEPS
00015 #define DEPS         0.02   /* depth epsilon */
00016 #endif
00017 #ifndef PEPS
00018 #define PEPS         0.04   /* pixel value epsilon */
00019 #endif
00020 #ifndef MAXRAD
00021 #define MAXRAD              64     /* maximum kernel radius */
00022 #endif
00023 #ifndef NNEIGH
00024 #define NNEIGH              5      /* find this many neighbors */
00025 #endif
00026 
00027 #define NINF         16382
00028 
00029 #define MAXRAD2             (MAXRAD*MAXRAD+1)
00030 
00031 #define G0NORM              0.286  /* ground zero normalization (1/x integral) */
00032 
00033 #ifndef FL4OP
00034 #define FL4OP(f,i,op)       ((f)[(i)>>5] op (1L<<((i)&0x1f)))
00035 #define CHK4(f,i)    FL4OP(f,i,&)
00036 #define SET4(f,i)    FL4OP(f,i,|=)
00037 #define CLR4(f,i)    FL4OP(f,i,&=~)
00038 #define TGL4(f,i)    FL4OP(f,i,^=)
00039 #define FL4NELS(n)   (((n)+0x1f)>>5)
00040 #define CLR4ALL(f,n) memset((char *)(f),'\0',FL4NELS(n)*sizeof(int32))
00041 #endif
00042 
00043 static int32  *pixFlags;    /* pixel occupancy flags */
00044 static float  pixWeight[MAXRAD2];  /* pixel weighting function */
00045 static short  isqrttab[MAXRAD2];   /* integer square root table */
00046 
00047 #define isqrt(i2)    (isqrttab[i2])
00048 
00049 extern VIEW   myview;              /* current output view */
00050 extern COLOR  *mypixel;     /* pixels being rendered */
00051 extern float  *myweight;    /* weights (used to compute final pixels) */
00052 extern float  *mydepth;     /* depth values (visibility culling) */
00053 extern int    hres, vres;   /* current horizontal and vertical res. */
00054 
00055 extern void pixFinish(double ransamp);
00056 extern void pixBeam(BEAM *bp, HDBEAMI *hb);
00057 
00058 typedef int (sampfunc)(int h, int v, short nl[NNEIGH][2], int nd[NNEIGH],
00059               int n, double *rf);
00060 static sampfunc kill_occl;
00061 static sampfunc smooth_samp;
00062 static sampfunc random_samp;
00063 static void meet_neighbors(int occ, sampfunc *nf, double *dp);
00064 static void reset_flags(void);
00065 static void init_wfunc(void);
00066 static int findneigh(short nl[NNEIGH][2], int nd[NNEIGH], int h, int v,
00067               register short (*rnl)[NNEIGH]);
00068 
00069 
00070 extern void
00071 pixBeam(                    /* render a particular beam */
00072        BEAM   *bp,
00073        register HDBEAMI     *hb
00074 )
00075 {
00076        GCOORD gc[2];
00077        register RAYVAL      *rv;
00078        FVECT  rorg, rdir, wp, ip;
00079        double d, prox;
00080        COLOR  col;
00081        int    n;
00082        register int32       p;
00083 
00084        if (!hdbcoord(gc, hb->h, hb->b))
00085               error(CONSISTENCY, "bad beam in render_beam");
00086        for (n = bp->nrm, rv = hdbray(bp); n--; rv++) {
00087                                           /* reproject each sample */
00088               hdray(rorg, rdir, hb->h, gc, rv->r);
00089               if (rv->d < DCINF) {
00090                      d = hddepth(hb->h, rv->d);
00091                      VSUM(wp, rorg, rdir, d);
00092                      VSUB(ip, wp, myview.vp);
00093                      d = DOT(ip,rdir);
00094                      prox = d*d/DOT(ip,ip);      /* cos(diff_angle)^32 */
00095                      prox *= prox; prox *= prox; prox *= prox; prox *= prox;
00096               } else {
00097                      if (myview.type == VT_PAR || myview.vaft > FTINY)
00098                             continue;            /* inf. off view */
00099                      VSUM(wp, myview.vp, rdir, FHUGE);
00100                      prox = 1.;
00101               }
00102               viewloc(ip, &myview, wp);   /* frustum clipping */
00103               if (ip[2] < 0.)
00104                      continue;
00105               if (ip[0] < 0. || ip[0] >= 1.)
00106                      continue;
00107               if (ip[1] < 0. || ip[1] >= 1.)
00108                      continue;
00109               if (myview.vaft > FTINY && ip[2] > myview.vaft - myview.vfore)
00110                      continue;            /* not exact for VT_PER */
00111               p = (int)(ip[1]*vres)*hres + (int)(ip[0]*hres);
00112               if (mydepth[p] > FTINY) {   /* check depth */
00113                      if (ip[2] > mydepth[p]*(1.+DEPS))
00114                             continue;
00115                      if (ip[2] < mydepth[p]*(1.-DEPS)) {
00116                             setcolor(mypixel[p], 0., 0., 0.);
00117                             myweight[p] = 0.;
00118                      }
00119               }
00120               colr_color(col, rv->v);
00121               scalecolor(col, prox);
00122               addcolor(mypixel[p], col);
00123               myweight[p] += prox;
00124               mydepth[p] = ip[2];
00125        }
00126 }
00127 
00128 
00129 static int
00130 kill_occl(           /* check for occlusion errors */
00131        int    h,
00132        int    v,
00133        register short       nl[NNEIGH][2],
00134        int    nd[NNEIGH],
00135        int    n,
00136        double *rf
00137 )
00138 {
00139        short  forequad[2][2];
00140        int    d;
00141        register int  i;
00142        register int32       p;
00143 
00144        if (n <= 0) {
00145 #ifdef DEBUG
00146               error(WARNING, "neighborless sample in kill_occl");
00147 #endif
00148               return(1);
00149        }
00150        p = v*hres + h;
00151        forequad[0][0] = forequad[0][1] = forequad[1][0] = forequad[1][1] = 0;
00152        for (i = n; i--; ) {
00153               d = isqrt(nd[i]);
00154               if (mydepth[nl[i][1]*hres+nl[i][0]]*(1.+DEPS*d) < mydepth[p])
00155                      forequad[nl[i][0]<h][nl[i][1]<v] = 1;
00156        }
00157        if (forequad[0][0]+forequad[0][1]+forequad[1][0]+forequad[1][1] > 2) {
00158               setcolor(mypixel[p], 0., 0., 0.);
00159               myweight[p] = 0.;    /* occupancy reset afterwards */
00160        }
00161        return(1);
00162 }
00163 
00164 
00165 static int
00166 smooth_samp(  /* grow sample point smoothly */
00167        int    h,
00168        int    v,
00169        register short       nl[NNEIGH][2],
00170        int    nd[NNEIGH],
00171        int    n,
00172        double *rf
00173 )
00174 {
00175        int    dis[NNEIGH], ndis;
00176        COLOR  mykern[MAXRAD2];
00177        int    maxr2;
00178        double d;
00179        register int32       p;
00180        register int  r2;
00181        int    i, r, maxr, h2, v2;
00182 
00183        if (n <= 0)
00184               return(1);
00185        p = v*hres + h;                           /* build kernel values */
00186        maxr2 = nd[n-1];
00187        DCHECK(maxr2>=MAXRAD2, CONSISTENCY, "out of range neighbor");
00188        maxr = isqrt(maxr2);
00189        for (v2 = 1; v2 <= maxr; v2++)
00190               for (h2 = 0; h2 <= v2; h2++) {
00191                      r2 = h2*h2 + v2*v2;
00192                      if (r2 > maxr2) break;
00193                      copycolor(mykern[r2], mypixel[p]);
00194                      scalecolor(mykern[r2], pixWeight[r2]);
00195               }
00196        ndis = 0;                          /* find discontinuities */
00197        for (i = n; i--; ) {
00198               r = isqrt(nd[i]);
00199               d = mydepth[nl[i][1]*hres+nl[i][0]] / mydepth[p];
00200               d = d>=1. ? d-1. : 1.-d;
00201               if (d > r*DEPS || bigdiff(mypixel[p],
00202                             mypixel[nl[i][1]*hres+nl[i][0]], r*PEPS))
00203                      dis[ndis++] = i;
00204        }
00205                                           /* stamp out that kernel */
00206        for (v2 = v-maxr; v2 <= v+maxr; v2++) {
00207               if (v2 < 0) v2 = 0;
00208               else if (v2 >= vres) break;
00209               for (h2 = h-maxr; h2 <= h+maxr; h2++) {
00210                      if (h2 < 0) h2 = 0;
00211                      else if (h2 >= hres) break;
00212                      r2 = (h2-h)*(h2-h) + (v2-v)*(v2-v);
00213                      if (r2 > maxr2) continue;
00214                      if (CHK4(pixFlags, v2*hres+h2))
00215                             continue;     /* occupied */
00216                      for (i = ndis; i--; ) {
00217                             r = (h2-nl[dis[i]][0])*(h2-nl[dis[i]][0]) +
00218                                    (v2-nl[dis[i]][1])*(v2-nl[dis[i]][1]);
00219                             if (r < r2) break;
00220                      }
00221                      if (i >= 0) continue;       /* outside edge */
00222                      addcolor(mypixel[v2*hres+h2], mykern[r2]);
00223                      myweight[v2*hres+h2] += pixWeight[r2] * myweight[p];
00224               }
00225        }
00226        return(1);
00227 }
00228 
00229 
00230 static int
00231 random_samp(  /* gather samples randomly */
00232        int    h,
00233        int    v,
00234        register short       nl[NNEIGH][2],
00235        int    nd[NNEIGH],
00236        int    n,
00237        double *rf
00238 )
00239 {
00240        float  rnt[NNEIGH];
00241        double rvar;
00242        register int32       p, pn;
00243        register int  ni;
00244 
00245        if (n <= 0)
00246               return(1);
00247        p = v*hres + h;
00248        if (*rf <= FTINY)           /* straight Voronoi regions */
00249               ni = 0;
00250        else {                      /* weighted choice */
00251               DCHECK(nd[n-1]>=MAXRAD2, CONSISTENCY, "out of range neighbor");
00252               rnt[0] = pixWeight[nd[0]];
00253               for (ni = 1; ni < n; ni++)
00254                      rnt[ni] = rnt[ni-1] + pixWeight[nd[ni]];
00255               rvar = rnt[n-1]*pow(frandom(), 1. / *rf);
00256               for (ni = 0; rvar > rnt[ni]+FTINY; ni++)
00257                      ;
00258        }
00259        pn = nl[ni][1]*hres + nl[ni][0];
00260        addcolor(mypixel[p], mypixel[pn]);
00261        myweight[p] += myweight[pn];
00262        return(1);
00263 }
00264 
00265 
00266 extern void
00267 pixFinish(           /* done with beams -- compute pixel values */
00268        double ransamp
00269 )
00270 {
00271        if (pixWeight[0] <= FTINY)
00272               init_wfunc();        /* initialize weighting function */
00273        reset_flags();                     /* set occupancy flags */
00274        meet_neighbors(1,kill_occl,NULL); /* identify occlusion errors */
00275        reset_flags();                     /* reset occupancy flags */
00276        if (ransamp >= 0.)          /* spread samples over image */
00277               meet_neighbors(0,random_samp,&ransamp);
00278        else
00279               meet_neighbors(1,smooth_samp,NULL);
00280        free((void *)pixFlags);            /* free pixel flags */
00281        pixFlags = NULL;
00282 }
00283 
00284 
00285 static void
00286 reset_flags(void)                  /* allocate/set/reset occupancy flags */
00287 {
00288        register int32       p;
00289 
00290        if (pixFlags == NULL) {
00291               pixFlags = (int32 *)calloc(FL4NELS(hres*vres), sizeof(int32));
00292               CHECK(pixFlags==NULL, SYSTEM, "out of memory in reset_flags");
00293        } else
00294               CLR4ALL(pixFlags, hres*vres);
00295        for (p = hres*vres; p--; )
00296               if (myweight[p] > FTINY)
00297                      SET4(pixFlags, p);
00298 }
00299 
00300 
00301 static void
00302 init_wfunc(void)                   /* initialize weighting function */
00303 {
00304        register int  r2;
00305        register double      d;
00306 
00307        for (r2 = MAXRAD2; --r2; ) {
00308               d = sqrt((double)r2);
00309               pixWeight[r2] = G0NORM/d;
00310               isqrttab[r2] = d + 0.99;
00311        }
00312        pixWeight[0] = 1.;
00313        isqrttab[0] = 0;
00314 }
00315 
00316 
00317 static int
00318 findneigh(    /* find NNEIGH neighbors for pixel */
00319        short  nl[NNEIGH][2],
00320        int    nd[NNEIGH],
00321        int    h,
00322        int    v,
00323        register short       (*rnl)[NNEIGH]
00324 )
00325 {
00326        int    nn = 0;
00327        int    d, n, hoff;
00328        register int  h2, n2;
00329 
00330        nd[NNEIGH-1] = MAXRAD2;
00331        for (hoff = 0; hoff < hres; hoff = (hoff<=0) - hoff) {
00332               h2 = h + hoff;
00333               if ((h2 < 0) | (h2 >= hres))
00334                      continue;
00335               if ((h2-h)*(h2-h) >= nd[NNEIGH-1])
00336                      break;
00337               for (n = 0; n < NNEIGH && rnl[h2][n] < NINF; n++) {
00338                      d = (h2-h)*(h2-h) + (v-rnl[h2][n])*(v-rnl[h2][n]);
00339                      if ((d == 0) | (d >= nd[NNEIGH-1]))
00340                             continue;
00341                      if (nn < NNEIGH)     /* insert neighbor */
00342                             nn++;
00343                      for (n2 = nn; n2--; )       {
00344                             if (!n2 || d >= nd[n2-1]) {
00345                                    nd[n2] = d;
00346                                    nl[n2][0] = h2;
00347                                    nl[n2][1] = rnl[h2][n];
00348                                    break;
00349                             }
00350                             nd[n2] = nd[n2-1];
00351                             nl[n2][0] = nl[n2-1][0];
00352                             nl[n2][1] = nl[n2-1][1];
00353                      }
00354               }
00355        }
00356        return(nn);
00357 }
00358 
00359 
00360 static void
00361 meet_neighbors(      /* run through samples and their neighbors */
00362        int    occ,
00363        sampfunc *nf,
00364        double *dp
00365 )
00366 {
00367        short  ln[NNEIGH][2];
00368        int    nd[NNEIGH];
00369        int    h, v, n, v2;
00370        register short       (*rnl)[NNEIGH];
00371                                    /* initialize bottom row list */
00372        rnl = (short (*)[NNEIGH])malloc(NNEIGH*sizeof(short)*hres);
00373        CHECK(rnl==NULL, SYSTEM, "out of memory in meet_neighbors");
00374        for (h = 0; h < hres; h++) {
00375               for (n = v = 0; v < vres; v++)
00376                      if (CHK4(pixFlags, v*hres+h)) {
00377                             rnl[h][n++] = v;
00378                             if (n >= NNEIGH)
00379                                    break;
00380                      }
00381               while (n < NNEIGH)
00382                      rnl[h][n++] = NINF;
00383        }
00384        v = 0;                      /* do each row */
00385        for ( ; ; ) {
00386               for (h = 0; h < hres; h++) {
00387                      if (!CHK4(pixFlags, v*hres+h) != !occ)
00388                             continue;     /* occupancy mismatch */
00389                                           /* find neighbors */
00390                      n = findneigh(ln, nd, h, v, rnl);
00391                                           /* call on neighbors */
00392                      (*nf)(h, v, ln, nd, n, dp);
00393               }
00394               if (++v >= vres)            /* reinitialize row list */
00395                      break;
00396               for (h = 0; h < hres; h++)
00397                      for (v2 = rnl[h][NNEIGH-1]+1; v2 < vres; v2++) {
00398                             if (v2 - v > v - rnl[h][0])
00399                                    break;        /* not close enough */
00400                             if (CHK4(pixFlags, v2*hres+h)) {
00401                                    for (n = 0; n < NNEIGH-1; n++)
00402                                           rnl[h][n] = rnl[h][n+1];
00403                                    rnl[h][NNEIGH-1] = v2;
00404                             }
00405                      }
00406        }
00407        free((void *)rnl);          /* free row list */
00408 }