--- ray/src/px/pf3.c 1993/06/25 15:01:20 2.7 +++ ray/src/px/pf3.c 2004/03/28 20:33:14 2.18 @@ -1,9 +1,6 @@ -/* Copyright (c) 1992 Regents of the University of California */ - #ifndef lint -static char SCCSid[] = "$SunId$ LBL"; +static const char RCSid[] = "$Id: pf3.c,v 2.18 2004/03/28 20:33:14 schorsch Exp $"; #endif - /* * pf3.c - routines for gaussian and box filtering * @@ -12,38 +9,15 @@ static char SCCSid[] = "$SunId$ LBL"; #include "standard.h" +#include + #include "color.h" +#include "pfilt.h" -#define TEPS 0.25 /* threshold proximity goal */ +#define RSCA 1.13 /* square-radius multiplier: sqrt(4/PI) */ +#define TEPS 0.2 /* threshold proximity goal */ +#define REPS 0.1 /* radius proximity goal */ -extern double CHECKRAD; /* radius over which gaussian is summed */ - -extern double rad; /* output pixel radius for filtering */ - -extern double thresh; /* maximum contribution for subpixel */ - -extern int nrows; /* number of rows for output */ -extern int ncols; /* number of columns for output */ - -extern int xres, yres; /* resolution of input */ - -extern double x_c, y_r; /* conversion factors */ - -extern int xrad; /* x search radius */ -extern int yrad; /* y search radius */ -extern int xbrad; /* x box size */ -extern int ybrad; /* y box size */ - -extern int barsize; /* size of input scan bar */ -extern COLOR **scanin; /* input scan bar */ -extern COLOR *scanout; /* output scan line */ -extern COLOR **scoutbar; /* output scan bar (if thresh > 0) */ -extern float **greybar; /* grey-averaged input values */ -extern int obarsize; /* size of output scan bar */ -extern int orad; /* output window radius */ - -extern char *progname; - float *gausstable; /* gauss lookup table */ float *ringsum; /* sum of ring values */ @@ -51,33 +25,37 @@ short *ringwt; /* weight (count) of ring values */ short *ringndx; /* ring index table */ float *warr; /* array of pixel weights */ -double pickfilt(); - #define lookgauss(x) gausstable[(int)(10.*(x)+.5)] +static double pickfilt(double p0); +static void sumans(int px, int py, int rcent, int ccent, double m); -initmask() /* initialize gaussian lookup table */ + +extern void +initmask(void) /* initialize gaussian lookup table */ { int gtabsiz; double gaussN; 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++) @@ -89,12 +67,12 @@ initmask() /* initialize gaussian lookup table */ for (x = 2*orad*orad+1; --x > orad*orad; ) ringndx[x] = -1; do - ringndx[x] = sqrt((double)x) + 0.5; + ringndx[x] = sqrt((double)x); while (x--); ringsum = (float *)malloc((orad+1)*sizeof(float)); ringwt = (short *)malloc((orad+1)*sizeof(short)); warr = (float *)malloc(obarsize*obarsize*sizeof(float)); - if (ringsum == NULL | ringwt == 0 | warr == NULL) + if ((ringsum == NULL) | (ringwt == 0) | (warr == NULL)) goto memerr; return; memerr: @@ -103,15 +81,19 @@ memerr: } -dobox(csum, xcent, ycent, c, r) /* simple box filter */ -COLOR csum; -int xcent, ycent; -int c, r; +extern void +dobox( /* simple box filter */ + COLOR csum, + int xcent, + int ycent, + int c, + int r +) { int wsum; double d; int y; - register int x; + register int x, offs; register COLOR *scan; wsum = 0; @@ -119,34 +101,41 @@ 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); + } } -dogauss(csum, xcent, ycent, c, r) /* gaussian filter */ -COLOR csum; -int xcent, ycent; -int c, r; +extern void +dogauss( /* gaussian filter */ + COLOR csum, + int xcent, + int ycent, + int c, + int r +) { double dy, dx, weight, wsum; COLOR ctmp; int y; - register int x; + register int x, offs; register COLOR *scan; wsum = FTINY; @@ -157,42 +146,49 @@ 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); } -dothresh(xcent, ycent, ccent, rcent) /* gaussian threshold filter */ -int xcent, ycent; -int ccent, rcent; +extern void +dothresh( /* gaussian threshold filter */ + int xcent, + int ycent, + int ccent, + int 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)); + memset((char *)ringsum, '\0', (orad+1)*sizeof(float)); + memset((char *)ringwt, '\0', (orad+1)*sizeof(short)); for (r = -orad; r <= orad; r++) { if (rcent+r < 0) continue; - if (rcent+r >= ncols) break; + 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]++; } } @@ -200,76 +196,105 @@ 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 } -double -pickfilt(p0) /* find filter multiplier for p0 */ -double p0; +static double +pickfilt( /* find filter multiplier for p0 */ + double p0 +) { double m = 1.0; - double t, avg, wsum; - int ilimit = 5/TEPS; + double t, num, denom, avg, wsum; + double mlimit[2]; + 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; - if (i > orad) i = orad+1; + i = RSCA*CHECKRAD*rad*m + .5; + if (i > orad) i = orad; avg = wsum = 0.0; while (i--) { - t = i/(m*rad); + t = (i+.5)/(m*rad); t = lookgauss(t*t); avg += t*ringsum[i]; wsum += t*ringwt[i]; } + if (avg < 1e-20) /* zero inclusive average */ + return(1.0); avg /= wsum; /* check threshold */ - t = p0 - avg; - if (t < 0.0) t = -t; - t *= gausstable[0]/(m*m*avg); - if (t <= thresh && (m <= 1.0+FTINY || - (thresh-t)/thresh <= TEPS)) - break; - if (t > thresh && (m*CHECKRAD*rad >= orad-FTINY || - (t-thresh)/thresh <= TEPS)) - break; + denom = m*m/gausstable[0] - p0/avg; + if (denom <= FTINY) { /* zero exclusive average */ + if (m >= mlimit[1]-REPS) + break; + m = mlimit[1]; + continue; + } + num = p0/avg - 1.0; + if (num < 0.0) num = -num; + t = num/denom; + if (t <= thresh) { + if (m <= mlimit[0]+REPS || (thresh-t)/thresh <= TEPS) + break; + } else { + if (m >= mlimit[1]-REPS || (t-thresh)/thresh <= TEPS) + break; + } + t = m; /* remember current m */ /* next guesstimate */ - m *= sqrt(t/thresh); - if (m < 1.0) m = 1.0; - else if (m*CHECKRAD*rad > orad) m = orad/rad/CHECKRAD; + m = sqrt(gausstable[0]*(num/thresh + p0/avg)); + if (m < t) { /* bound it */ + if (m <= mlimit[0]+FTINY) + m = 0.5*(mlimit[0] + t); + mlimit[1] = t; + } else { + if (m >= mlimit[1]-FTINY) + m = 0.5*(mlimit[1] + t); + mlimit[0] = t; + } } while (--ilimit > 0); return(m); } -sumans(px, py, rcent, ccent, m) /* sum input pixel to output */ -int px, py; -int rcent, ccent; -double m; +static void +sumans( /* sum input pixel to output */ + int px, + int py, + int rcent, + int 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; @@ -280,12 +305,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; @@ -298,12 +326,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); } } }