--- ray/src/px/pf3.c 1993/07/16 11:32:06 2.11 +++ ray/src/px/pf3.c 2003/02/25 00:26:05 2.15 @@ -1,9 +1,6 @@ -/* Copyright (c) 1993 Regents of the University of California */ - #ifndef lint -static char SCCSid[] = "$SunId$ LBL"; +static const char RCSid[] = "$Id: pf3.c,v 2.15 2003/02/25 00:26:05 greg Exp $"; #endif - /* * pf3.c - routines for gaussian and box filtering * @@ -14,6 +11,7 @@ static char SCCSid[] = "$SunId$ LBL"; #include "color.h" +#define RSCA 1.13 /* square-radius multiplier: sqrt(4/PI) */ #define TEPS 0.2 /* threshold proximity goal */ #define REPS 0.1 /* radius proximity goal */ @@ -43,6 +41,8 @@ extern float **greybar; /* grey-averaged input values extern int obarsize; /* size of output scan bar */ extern int orad; /* output window radius */ +extern int wrapfilt; /* wrap filter horizontally? */ + extern char *progname; float *gausstable; /* gauss lookup table */ @@ -52,6 +52,8 @@ short *ringwt; /* weight (count) of ring values */ short *ringndx; /* ring index table */ float *warr; /* array of pixel weights */ +extern double (*ourbright)(); /* brightness computation function */ + double pickfilt(); #define lookgauss(x) gausstable[(int)(10.*(x)+.5)] @@ -64,21 +66,23 @@ initmask() /* initialize gaussian lookup table */ double d; register int x; - gtabsiz = 150*CHECKRAD; + gtabsiz = 111*CHECKRAD*CHECKRAD; gausstable = (float *)malloc(gtabsiz*sizeof(float)); if (gausstable == NULL) goto memerr; d = x_c*y_r*0.25/(rad*rad); - gausstable[0] = exp(-d)/sqrt(d); + gausstable[0] = exp(-d); for (x = 1; x < gtabsiz; x++) - if ((gausstable[x] = exp(-x*0.1)/sqrt(x*0.1)) > gausstable[0]) + if (x*0.1 <= d) gausstable[x] = gausstable[0]; + else + gausstable[x] = exp(-x*0.1); if (obarsize == 0) return; /* compute integral of filter */ - gaussN = PI*sqrt(d)*exp(-d); /* plateau */ - for (d = sqrt(d)+0.05; d <= 1.25*CHECKRAD; d += 0.1) - gaussN += 0.1*2.0*PI*exp(-d*d); + gaussN = PI*d*exp(-d); /* plateau */ + for (d = sqrt(d)+0.05; d <= RSCA*CHECKRAD; d += 0.1) + gaussN += 0.1*2.0*PI*d*exp(-d*d); /* normalize filter */ gaussN = x_c*y_r/(rad*rad*gaussN); for (x = 0; x < gtabsiz; x++) @@ -112,7 +116,7 @@ int c, r; int wsum; double d; int y; - register int x; + register int x, offs; register COLOR *scan; wsum = 0; @@ -120,22 +124,25 @@ int c, r; for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) { if (y < 0) continue; if (y >= yres) break; - d = y_r < 1.0 ? y_r*y - r : (double)(y - ycent); + d = y_r < 1.0 ? y_r*y - (r+.5) : (double)(y - ycent); if (d < -0.5) continue; if (d >= 0.5) break; scan = scanin[y%barsize]; for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) { - if (x < 0) continue; - if (x >= xres) break; - d = x_c < 1.0 ? x_c*x - c : (double)(x - xcent); + offs = x < 0 ? xres : x >= xres ? -xres : 0; + if (offs && !wrapfilt) + continue; + d = x_c < 1.0 ? x_c*x - (c+.5) : (double)(x - xcent); if (d < -0.5) continue; if (d >= 0.5) break; wsum++; - addcolor(csum, scan[x]); + addcolor(csum, scan[x+offs]); } } - if (wsum > 1) - scalecolor(csum, 1.0/wsum); + if (wsum > 1) { + d = 1.0/wsum; + scalecolor(csum, d); + } } @@ -147,7 +154,7 @@ int c, r; double dy, dx, weight, wsum; COLOR ctmp; int y; - register int x; + register int x, offs; register COLOR *scan; wsum = FTINY; @@ -158,17 +165,19 @@ int c, r; dy = (y_r*(y+.5) - (r+.5))/rad; scan = scanin[y%barsize]; for (x = xcent-xrad; x <= xcent+xrad; x++) { - if (x < 0) continue; - if (x >= xres) break; + offs = x < 0 ? xres : x >= xres ? -xres : 0; + if (offs && !wrapfilt) + continue; dx = (x_c*(x+.5) - (c+.5))/rad; weight = lookgauss(dx*dx + dy*dy); wsum += weight; - copycolor(ctmp, scan[x]); + copycolor(ctmp, scan[x+offs]); scalecolor(ctmp, weight); addcolor(csum, ctmp); } } - scalecolor(csum, 1.0/wsum); + weight = 1.0/wsum; + scalecolor(csum, weight); } @@ -177,10 +186,9 @@ int xcent, ycent; int ccent, rcent; { double d; - int r, y; + int r, y, offs; register int c, x; register float *gscan; -#define pval gscan /* compute ring sums */ bzero((char *)ringsum, (orad+1)*sizeof(float)); bzero((char *)ringwt, (orad+1)*sizeof(short)); @@ -189,11 +197,13 @@ int ccent, rcent; if (rcent+r >= nrows) break; gscan = greybar[(rcent+r)%obarsize]; for (c = -orad; c <= orad; c++) { - if (ccent+c < 0) continue; - if (ccent+c >= ncols) break; + offs = ccent+c < 0 ? ncols : + ccent+c >= ncols ? -ncols : 0; + if (offs && !wrapfilt) + continue; x = ringndx[c*c + r*r]; if (x < 0) continue; - ringsum[x] += gscan[ccent+c]; + ringsum[x] += gscan[ccent+c+offs]; ringwt[x]++; } } @@ -201,20 +211,20 @@ int ccent, rcent; for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) { if (y < 0) continue; if (y >= yres) break; - d = y_r < 1.0 ? y_r*y - rcent : (double)(y - ycent); + d = y_r < 1.0 ? y_r*y - (rcent+.5) : (double)(y - ycent); if (d < -0.5) continue; if (d >= 0.5) break; for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) { - if (x < 0) continue; - if (x >= xres) break; - d = x_c < 1.0 ? x_c*x - ccent : (double)(x - xcent); + offs = x < 0 ? xres : x >= xres ? -xres : 0; + if (offs && !wrapfilt) + continue; + d = x_c < 1.0 ? x_c*x - (ccent+.5) : (double)(x - xcent); if (d < -0.5) continue; if (d >= 0.5) break; - pval = scanin[y%barsize][x]; - sumans(x, y, rcent, ccent, pickfilt(bright(pval))); + sumans(x, y, rcent, ccent, + pickfilt((*ourbright)(scanin[y%barsize][x+offs]))); } } -#undef pval } @@ -225,13 +235,13 @@ double p0; double m = 1.0; double t, num, denom, avg, wsum; double mlimit[2]; - int ilimit = 4/TEPS; + int ilimit = 4.0/TEPS; register int i; /* iterative search for m */ mlimit[0] = 1.0; mlimit[1] = orad/rad/CHECKRAD; do { /* compute grey weighted average */ - i = 1.25*CHECKRAD*rad*m + .5; + i = RSCA*CHECKRAD*rad*m + .5; if (i > orad) i = orad; avg = wsum = 0.0; while (i--) { @@ -283,14 +293,18 @@ int px, py; int rcent, ccent; double m; { - double dy, dx; + double dy2, dx; COLOR pval, ctmp; - int ksiz, r; + int ksiz, r, offs; double pc, pr, norm; register int i, c; register COLOR *scan; - - copycolor(pval, scanin[py%barsize][px]); + /* + * This normalization method fails at the picture borders because + * a different number of input pixels contribute there. + */ + scan = scanin[py%barsize] + (px < 0 ? xres : px >= xres ? -xres : 0); + copycolor(pval, scan[px]); pc = x_c*(px+.5); pr = y_r*(py+.5); ksiz = CHECKRAD*m*rad + 1; @@ -301,12 +315,15 @@ double m; for (r = rcent-ksiz; r <= rcent+ksiz; r++) { if (r < 0) continue; if (r >= nrows) break; - dy = (pr - (r+.5))/(m*rad); + dy2 = (pr - (r+.5))/(m*rad); + dy2 *= dy2; for (c = ccent-ksiz; c <= ccent+ksiz; c++) { - if (c < 0) continue; - if (c >= ncols) break; + if (!wrapfilt) { + if (c < 0) continue; + if (c >= ncols) break; + } dx = (pc - (c+.5))/(m*rad); - norm += warr[i++] = lookgauss(dx*dx + dy*dy); + norm += warr[i++] = lookgauss(dx*dx + dy2); } } norm = 1.0/norm; @@ -319,12 +336,13 @@ double m; if (r >= nrows) break; scan = scoutbar[r%obarsize]; for (c = ccent-ksiz; c <= ccent+ksiz; c++) { - if (c < 0) continue; - if (c >= ncols) break; + offs = c < 0 ? ncols : c >= ncols ? -ncols : 0; + if (offs && !wrapfilt) + continue; copycolor(ctmp, pval); dx = norm*warr[i++]; scalecolor(ctmp, dx); - addcolor(scan[c], ctmp); + addcolor(scan[c+offs], ctmp); } } }