Back to index

radiance  4R0+20100331
pf3.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char RCSid[] = "$Id: pf3.c,v 2.18 2004/03/28 20:33:14 schorsch Exp $";
00003 #endif
00004 /*
00005  *  pf3.c - routines for gaussian and box filtering
00006  *
00007  *     8/13/86
00008  */
00009 
00010 #include  "standard.h"
00011 
00012 #include  <string.h>
00013 
00014 #include  "color.h"
00015 #include  "pfilt.h"
00016 
00017 #define  RSCA        1.13   /* square-radius multiplier: sqrt(4/PI) */
00018 #define  TEPS        0.2    /* threshold proximity goal */
00019 #define  REPS        0.1    /* radius proximity goal */
00020 
00021 float  *gausstable;         /* gauss lookup table */
00022 
00023 float  *ringsum;            /* sum of ring values */
00024 short  *ringwt;                    /* weight (count) of ring values */
00025 short  *ringndx;            /* ring index table */
00026 float  *warr;               /* array of pixel weights */
00027 
00028 #define        lookgauss(x)        gausstable[(int)(10.*(x)+.5)]
00029 
00030 static double pickfilt(double  p0);
00031 static void sumans(int  px, int  py, int  rcent, int  ccent, double  m);
00032 
00033 
00034 extern void
00035 initmask(void)                     /* initialize gaussian lookup table */
00036 {
00037        int  gtabsiz;
00038        double  gaussN;
00039        double  d;
00040        register int  x;
00041 
00042        gtabsiz = 111*CHECKRAD*CHECKRAD;
00043        gausstable = (float *)malloc(gtabsiz*sizeof(float));
00044        if (gausstable == NULL)
00045               goto memerr;
00046        d = x_c*y_r*0.25/(rad*rad);
00047        gausstable[0] = exp(-d);
00048        for (x = 1; x < gtabsiz; x++)
00049               if (x*0.1 <= d)
00050                      gausstable[x] = gausstable[0];
00051               else
00052                      gausstable[x] = exp(-x*0.1);
00053        if (obarsize == 0)
00054               return;
00055                                    /* compute integral of filter */
00056        gaussN = PI*d*exp(-d);                    /* plateau */
00057        for (d = sqrt(d)+0.05; d <= RSCA*CHECKRAD; d += 0.1)
00058               gaussN += 0.1*2.0*PI*d*exp(-d*d);
00059                                    /* normalize filter */
00060        gaussN = x_c*y_r/(rad*rad*gaussN);
00061        for (x = 0; x < gtabsiz; x++)
00062               gausstable[x] *= gaussN;
00063                                    /* create ring averages table */
00064        ringndx = (short *)malloc((2*orad*orad+1)*sizeof(short));
00065        if (ringndx == NULL)
00066               goto memerr;
00067        for (x = 2*orad*orad+1; --x > orad*orad; )
00068               ringndx[x] = -1;
00069        do
00070               ringndx[x] = sqrt((double)x);
00071        while (x--);
00072        ringsum = (float *)malloc((orad+1)*sizeof(float));
00073        ringwt = (short *)malloc((orad+1)*sizeof(short));
00074        warr = (float *)malloc(obarsize*obarsize*sizeof(float));
00075        if ((ringsum == NULL) | (ringwt == 0) | (warr == NULL))
00076               goto memerr;
00077        return;
00078 memerr:
00079        fprintf(stderr, "%s: out of memory in initmask\n", progname);
00080        quit(1);
00081 }
00082 
00083 
00084 extern void
00085 dobox(               /* simple box filter */
00086        COLOR  csum,
00087        int  xcent,
00088        int  ycent,
00089        int  c,
00090        int  r
00091 )
00092 {
00093        int  wsum;
00094        double  d;
00095        int  y;
00096        register int  x, offs;
00097        register COLOR       *scan;
00098        
00099        wsum = 0;
00100        setcolor(csum, 0.0, 0.0, 0.0);
00101        for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
00102               if (y < 0) continue;
00103               if (y >= yres) break;
00104               d = y_r < 1.0 ? y_r*y - (r+.5) : (double)(y - ycent);
00105               if (d < -0.5) continue;
00106               if (d >= 0.5) break;
00107               scan = scanin[y%barsize];
00108               for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
00109                      offs = x < 0 ? xres : x >= xres ? -xres : 0;
00110                      if (offs && !wrapfilt)
00111                             continue;
00112                      d = x_c < 1.0 ? x_c*x - (c+.5) : (double)(x - xcent);
00113                      if (d < -0.5) continue;
00114                      if (d >= 0.5) break;
00115                      wsum++;
00116                      addcolor(csum, scan[x+offs]);
00117               }
00118        }
00119        if (wsum > 1) {
00120               d = 1.0/wsum;
00121               scalecolor(csum, d);
00122        }
00123 }
00124 
00125 
00126 extern void
00127 dogauss(             /* gaussian filter */
00128        COLOR  csum,
00129        int  xcent,
00130        int  ycent,
00131        int  c,
00132        int  r
00133 )
00134 {
00135        double  dy, dx, weight, wsum;
00136        COLOR  ctmp;
00137        int  y;
00138        register int  x, offs;
00139        register COLOR       *scan;
00140 
00141        wsum = FTINY;
00142        setcolor(csum, 0.0, 0.0, 0.0);
00143        for (y = ycent-yrad; y <= ycent+yrad; y++) {
00144               if (y < 0) continue;
00145               if (y >= yres) break;
00146               dy = (y_r*(y+.5) - (r+.5))/rad;
00147               scan = scanin[y%barsize];
00148               for (x = xcent-xrad; x <= xcent+xrad; x++) {
00149                      offs = x < 0 ? xres : x >= xres ? -xres : 0;
00150                      if (offs && !wrapfilt)
00151                             continue;
00152                      dx = (x_c*(x+.5) - (c+.5))/rad;
00153                      weight = lookgauss(dx*dx + dy*dy);
00154                      wsum += weight;
00155                      copycolor(ctmp, scan[x+offs]);
00156                      scalecolor(ctmp, weight);
00157                      addcolor(csum, ctmp);
00158               }
00159        }
00160        weight = 1.0/wsum;
00161        scalecolor(csum, weight);
00162 }
00163 
00164 
00165 extern void
00166 dothresh(     /* gaussian threshold filter */
00167        int  xcent,
00168        int  ycent,
00169        int  ccent,
00170        int  rcent
00171 )
00172 {
00173        double  d;
00174        int  r, y, offs;
00175        register int  c, x;
00176        register float  *gscan;
00177                                    /* compute ring sums */
00178        memset((char *)ringsum, '\0', (orad+1)*sizeof(float));
00179        memset((char *)ringwt, '\0', (orad+1)*sizeof(short));
00180        for (r = -orad; r <= orad; r++) {
00181               if (rcent+r < 0) continue;
00182               if (rcent+r >= nrows) break;
00183               gscan = greybar[(rcent+r)%obarsize];
00184               for (c = -orad; c <= orad; c++) {
00185                      offs = ccent+c < 0 ? ncols :
00186                                    ccent+c >= ncols ? -ncols : 0;
00187                      if (offs && !wrapfilt)
00188                             continue;
00189                      x = ringndx[c*c + r*r];
00190                      if (x < 0) continue;
00191                      ringsum[x] += gscan[ccent+c+offs];
00192                      ringwt[x]++;
00193               }
00194        }
00195                                    /* filter each subpixel */
00196        for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
00197               if (y < 0) continue;
00198               if (y >= yres) break;
00199               d = y_r < 1.0 ? y_r*y - (rcent+.5) : (double)(y - ycent);
00200               if (d < -0.5) continue;
00201               if (d >= 0.5) break;
00202               for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
00203                      offs = x < 0 ? xres : x >= xres ? -xres : 0;
00204                      if (offs && !wrapfilt)
00205                             continue;
00206                      d = x_c < 1.0 ? x_c*x - (ccent+.5) : (double)(x - xcent);
00207                      if (d < -0.5) continue;
00208                      if (d >= 0.5) break;
00209                      sumans(x, y, rcent, ccent,
00210                      pickfilt((*ourbright)(scanin[y%barsize][x+offs])));
00211               }
00212        }
00213 }
00214 
00215 
00216 static double
00217 pickfilt(                   /* find filter multiplier for p0 */
00218        double  p0
00219 )
00220 {
00221        double  m = 1.0;
00222        double  t, num, denom, avg, wsum;
00223        double  mlimit[2];
00224        int  ilimit = 4.0/TEPS;
00225        register int  i;
00226                             /* iterative search for m */
00227        mlimit[0] = 1.0; mlimit[1] = orad/rad/CHECKRAD;
00228        do {
00229                                    /* compute grey weighted average */
00230               i = RSCA*CHECKRAD*rad*m + .5;
00231               if (i > orad) i = orad;
00232               avg = wsum = 0.0;
00233               while (i--) {
00234                      t = (i+.5)/(m*rad);
00235                      t = lookgauss(t*t);
00236                      avg += t*ringsum[i];
00237                      wsum += t*ringwt[i];
00238               }
00239               if (avg < 1e-20)            /* zero inclusive average */
00240                      return(1.0);
00241               avg /= wsum;
00242                                    /* check threshold */
00243               denom = m*m/gausstable[0] - p0/avg;
00244               if (denom <= FTINY) {              /* zero exclusive average */
00245                      if (m >= mlimit[1]-REPS)
00246                             break;
00247                      m = mlimit[1];
00248                      continue;
00249               }
00250               num = p0/avg - 1.0;
00251               if (num < 0.0) num = -num;
00252               t = num/denom;
00253               if (t <= thresh) {
00254                      if (m <= mlimit[0]+REPS || (thresh-t)/thresh <= TEPS)
00255                             break;
00256               } else {
00257                      if (m >= mlimit[1]-REPS || (t-thresh)/thresh <= TEPS)
00258                             break;
00259               }
00260               t = m;               /* remember current m */
00261                                    /* next guesstimate */
00262               m = sqrt(gausstable[0]*(num/thresh + p0/avg));
00263               if (m < t) {         /* bound it */
00264                      if (m <= mlimit[0]+FTINY)
00265                             m = 0.5*(mlimit[0] + t);
00266                      mlimit[1] = t;
00267               } else {
00268                      if (m >= mlimit[1]-FTINY)
00269                             m = 0.5*(mlimit[1] + t);
00270                      mlimit[0] = t;
00271               }
00272        } while (--ilimit > 0);
00273        return(m);
00274 }
00275 
00276 
00277 static void
00278 sumans(              /* sum input pixel to output */
00279        int  px,
00280        int  py,
00281        int  rcent,
00282        int  ccent,
00283        double  m
00284 )
00285 {
00286        double  dy2, dx;
00287        COLOR  pval, ctmp;
00288        int  ksiz, r, offs;
00289        double  pc, pr, norm;
00290        register int  i, c;
00291        register COLOR       *scan;
00292        /*
00293         * This normalization method fails at the picture borders because
00294         * a different number of input pixels contribute there.
00295         */
00296        scan = scanin[py%barsize] + (px < 0 ? xres : px >= xres ? -xres : 0);
00297        copycolor(pval, scan[px]);
00298        pc = x_c*(px+.5);
00299        pr = y_r*(py+.5);
00300        ksiz = CHECKRAD*m*rad + 1;
00301        if (ksiz > orad) ksiz = orad;
00302                                           /* compute normalization */
00303        norm = 0.0;
00304        i = 0;
00305        for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
00306               if (r < 0) continue;
00307               if (r >= nrows) break;
00308               dy2 = (pr - (r+.5))/(m*rad);
00309               dy2 *= dy2;
00310               for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
00311                      if (!wrapfilt) {
00312                             if (c < 0) continue;
00313                             if (c >= ncols) break;
00314                      }
00315                      dx = (pc - (c+.5))/(m*rad);
00316                      norm += warr[i++] = lookgauss(dx*dx + dy2);
00317               }
00318        }
00319        norm = 1.0/norm;
00320        if (x_c < 1.0) norm *= x_c;
00321        if (y_r < 1.0) norm *= y_r;
00322                                           /* sum pixels */
00323        i = 0;
00324        for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
00325               if (r < 0) continue;
00326               if (r >= nrows) break;
00327               scan = scoutbar[r%obarsize];
00328               for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
00329                      offs = c < 0 ? ncols : c >= ncols ? -ncols : 0;
00330                      if (offs && !wrapfilt)
00331                             continue;
00332                      copycolor(ctmp, pval);
00333                      dx = norm*warr[i++];
00334                      scalecolor(ctmp, dx);
00335                      addcolor(scan[c+offs], ctmp);
00336               }
00337        }
00338 }