src/px/pf3.c

#ifndef lint
static const char RCSid[] = "$Id: pf3.c,v 2.17 2003/07/27 22:12:03 schorsch Exp $";
#endif
/*
 *  pf3.c - routines for gaussian and box filtering
 *
 *     8/13/86
 */

#include  "standard.h"

#include  <string.h>

#include  "color.h"
#include  "pfilt.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 */

float  *gausstable;             /* gauss lookup table */

float  *ringsum;                /* sum of ring values */
short  *ringwt;                 /* weight (count) of ring values */
short  *ringndx;                /* ring index table */
float  *warr;                   /* array of pixel weights */

#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);


extern void
initmask(void)                  /* initialize gaussian lookup table */
{
        int  gtabsiz;
        double  gaussN;
        double  d;
        register int  x;

        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);
        for (x = 1; x < gtabsiz; x++)
                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*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++)
                gausstable[x] *= gaussN;
                                        /* create ring averages table */
        ringndx = (short *)malloc((2*orad*orad+1)*sizeof(short));
        if (ringndx == NULL)
                goto memerr;
        for (x = 2*orad*orad+1; --x > orad*orad; )
                ringndx[x] = -1;
        do
                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))
                goto memerr;
        return;
memerr:
        fprintf(stderr, "%s: out of memory in initmask\n", progname);
        quit(1);
}


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, offs;
        register COLOR  *scan;
        
        wsum = 0;
        setcolor(csum, 0.0, 0.0, 0.0);
        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+.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++) {
                        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+offs]);
                }
        }
        if (wsum > 1) {
                d = 1.0/wsum;
                scalecolor(csum, d);
        }
}


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, offs;
        register COLOR  *scan;

        wsum = FTINY;
        setcolor(csum, 0.0, 0.0, 0.0);
        for (y = ycent-yrad; y <= ycent+yrad; y++) {
                if (y < 0) continue;
                if (y >= yres) break;
                dy = (y_r*(y+.5) - (r+.5))/rad;
                scan = scanin[y%barsize];
                for (x = xcent-xrad; x <= xcent+xrad; x++) {
                        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+offs]);
                        scalecolor(ctmp, weight);
                        addcolor(csum, ctmp);
                }
        }
        weight = 1.0/wsum;
        scalecolor(csum, weight);
}


extern void
dothresh(       /* gaussian threshold filter */
        int  xcent,
        int  ycent,
        int  ccent,
        int  rcent
)
{
        double  d;
        int  r, y, offs;
        register int  c, x;
        register float  *gscan;
                                        /* compute ring sums */
        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 >= nrows) break;
                gscan = greybar[(rcent+r)%obarsize];
                for (c = -orad; c <= orad; c++) {
                        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+offs];
                        ringwt[x]++;
                }
        }
                                        /* filter each subpixel */
        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+.5) : (double)(y - ycent);
                if (d < -0.5) continue;
                if (d >= 0.5) break;
                for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
                        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;
                        sumans(x, y, rcent, ccent,
                        pickfilt((*ourbright)(scanin[y%barsize][x+offs])));
                }
        }
}


static double
pickfilt(                       /* find filter multiplier for p0 */
        double  p0
)
{
        double  m = 1.0;
        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 = RSCA*CHECKRAD*rad*m + .5;
                if (i > orad) i = orad;
                avg = wsum = 0.0;
                while (i--) {
                        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 */
                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(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);
}


static void
sumans(         /* sum input pixel to output */
        int  px,
        int  py,
        int  rcent,
        int  ccent,
        double  m
)
{
        double  dy2, dx;
        COLOR  pval, ctmp;
        int  ksiz, r, offs;
        double  pc, pr, norm;
        register int  i, c;
        register COLOR  *scan;
        /*
         * 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;
        if (ksiz > orad) ksiz = orad;
                                                /* compute normalization */
        norm = 0.0;
        i = 0;
        for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
                if (r < 0) continue;
                if (r >= nrows) break;
                dy2 = (pr - (r+.5))/(m*rad);
                dy2 *= dy2;
                for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
                        if (!wrapfilt) {
                                if (c < 0) continue;
                                if (c >= ncols) break;
                        }
                        dx = (pc - (c+.5))/(m*rad);
                        norm += warr[i++] = lookgauss(dx*dx + dy2);
                }
        }
        norm = 1.0/norm;
        if (x_c < 1.0) norm *= x_c;
        if (y_r < 1.0) norm *= y_r;
                                                /* sum pixels */
        i = 0;
        for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
                if (r < 0) continue;
                if (r >= nrows) break;
                scan = scoutbar[r%obarsize];
                for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
                        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+offs], ctmp);
                }
        }
}
Revision:	2.18
Committed:	Sun Mar 28 20:33:14 2004 UTC (21 years, 7 months ago) by schorsch
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R0, rad5R1, rad4R2, rad3R7P2, rad3R7P1, rad4R1, rad4R0, rad3R6, rad3R6P1, rad3R8, rad3R9, rad4R2P1, rad4R2P2
Changes since 2.17:	+42 -54 lines
Log Message:	Continued ANSIfication, and other fixes and clarifications.
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	schorsch	2.18	static const char RCSid[] = "$Id: pf3.c,v 2.17 2003/07/27 22:12:03 schorsch Exp $";
3	greg	1.1	#endif
4			/*
5			* pf3.c - routines for gaussian and box filtering
6			*
7			* 8/13/86
8			*/
9
10	greg	2.7	#include "standard.h"
11	greg	1.1
12	schorsch	2.16	#include <string.h>
13
14	greg	1.1	#include "color.h"
15	schorsch	2.18	#include "pfilt.h"
16	greg	1.1
17	greg	2.14	#define RSCA 1.13 /* square-radius multiplier: sqrt(4/PI) */
18	greg	2.9	#define TEPS 0.2 /* threshold proximity goal */
19	greg	2.10	#define REPS 0.1 /* radius proximity goal */
20	greg	1.1
21	greg	2.5	float gausstable; / gauss lookup table */
22	greg	1.1
23	greg	2.7	float ringsum; / sum of ring values */
24			short ringwt; / weight (count) of ring values */
25			short ringndx; / ring index table */
26			float warr; / array of pixel weights */
27
28	schorsch	2.18	#define lookgauss(x) gausstable[(int)(10.*(x)+.5)]
29	greg	2.7
30	schorsch	2.18	static double pickfilt(double p0);
31			static void sumans(int px, int py, int rcent, int ccent, double m);
32	greg	1.1
33
34	schorsch	2.18	extern void
35			initmask(void) /* initialize gaussian lookup table */
36	greg	1.1	{
37	greg	2.7	int gtabsiz;
38			double gaussN;
39	greg	2.5	double d;
40	greg	1.1	register int x;
41
42	greg	2.15	gtabsiz = 111CHECKRADCHECKRAD;
43	greg	2.7	gausstable = (float )malloc(gtabsizsizeof(float));
44			if (gausstable == NULL)
45			goto memerr;
46			d = x_cy_r0.25/(rad*rad);
47	greg	2.14	gausstable[0] = exp(-d);
48	greg	2.7	for (x = 1; x < gtabsiz; x++)
49	greg	2.14	if (x*0.1 <= d)
50	greg	2.5	gausstable[x] = gausstable[0];
51	greg	2.14	else
52			gausstable[x] = exp(-x*0.1);
53	greg	2.7	if (obarsize == 0)
54			return;
55			/* compute integral of filter */
56	greg	2.14	gaussN = PIdexp(-d); /* plateau */
57			for (d = sqrt(d)+0.05; d <= RSCA*CHECKRAD; d += 0.1)
58			gaussN += 0.12.0PIdexp(-d*d);
59	greg	2.7	/* normalize filter */
60			gaussN = x_cy_r/(radrad*gaussN);
61			for (x = 0; x < gtabsiz; x++)
62			gausstable[x] *= gaussN;
63			/* create ring averages table */
64			ringndx = (short )malloc((2oradorad+1)sizeof(short));
65			if (ringndx == NULL)
66			goto memerr;
67			for (x = 2oradorad+1; --x > orad*orad; )
68			ringndx[x] = -1;
69			do
70	greg	2.8	ringndx[x] = sqrt((double)x);
71	greg	2.7	while (x--);
72			ringsum = (float )malloc((orad+1)sizeof(float));
73			ringwt = (short )malloc((orad+1)sizeof(short));
74			warr = (float )malloc(obarsizeobarsize*sizeof(float));
75	schorsch	2.17	if ((ringsum == NULL) \| (ringwt == 0) \| (warr == NULL))
76	greg	2.7	goto memerr;
77			return;
78			memerr:
79			fprintf(stderr, "%s: out of memory in initmask\n", progname);
80			quit(1);
81	greg	1.1	}
82
83
84	schorsch	2.18	extern void
85			dobox( /* simple box filter */
86			COLOR csum,
87			int xcent,
88			int ycent,
89			int c,
90			int r
91			)
92	greg	1.1	{
93	greg	2.6	int wsum;
94			double d;
95			int y;
96	greg	2.12	register int x, offs;
97	greg	2.2	register COLOR *scan;
98	greg	2.7
99	greg	1.1	wsum = 0;
100			setcolor(csum, 0.0, 0.0, 0.0);
101	greg	2.7	for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
102	greg	2.4	if (y < 0) continue;
103			if (y >= yres) break;
104	greg	2.13	d = y_r < 1.0 ? y_r*y - (r+.5) : (double)(y - ycent);
105	greg	2.7	if (d < -0.5) continue;
106			if (d >= 0.5) break;
107	greg	1.1	scan = scanin[y%barsize];
108	greg	2.7	for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
109	greg	2.12	offs = x < 0 ? xres : x >= xres ? -xres : 0;
110			if (offs && !wrapfilt)
111			continue;
112	greg	2.13	d = x_c < 1.0 ? x_c*x - (c+.5) : (double)(x - xcent);
113	greg	2.7	if (d < -0.5) continue;
114			if (d >= 0.5) break;
115	greg	1.1	wsum++;
116	greg	2.12	addcolor(csum, scan[x+offs]);
117	greg	1.1	}
118			}
119	greg	2.14	if (wsum > 1) {
120			d = 1.0/wsum;
121			scalecolor(csum, d);
122			}
123	greg	1.1	}
124
125
126	schorsch	2.18	extern void
127			dogauss( /* gaussian filter */
128			COLOR csum,
129			int xcent,
130			int ycent,
131			int c,
132			int r
133			)
134	greg	1.1	{
135	greg	2.6	double dy, dx, weight, wsum;
136			COLOR ctmp;
137			int y;
138	greg	2.12	register int x, offs;
139	greg	2.2	register COLOR *scan;
140	greg	1.1
141			wsum = FTINY;
142			setcolor(csum, 0.0, 0.0, 0.0);
143	greg	1.2	for (y = ycent-yrad; y <= ycent+yrad; y++) {
144	greg	2.4	if (y < 0) continue;
145			if (y >= yres) break;
146	greg	1.2	dy = (y_r*(y+.5) - (r+.5))/rad;
147	greg	1.1	scan = scanin[y%barsize];
148	greg	1.2	for (x = xcent-xrad; x <= xcent+xrad; x++) {
149	greg	2.12	offs = x < 0 ? xres : x >= xres ? -xres : 0;
150			if (offs && !wrapfilt)
151			continue;
152	greg	1.2	dx = (x_c*(x+.5) - (c+.5))/rad;
153	greg	2.6	weight = lookgauss(dxdx + dydy);
154	greg	1.1	wsum += weight;
155	greg	2.12	copycolor(ctmp, scan[x+offs]);
156	greg	1.1	scalecolor(ctmp, weight);
157			addcolor(csum, ctmp);
158			}
159			}
160	greg	2.14	weight = 1.0/wsum;
161			scalecolor(csum, weight);
162	greg	2.7	}
163
164
165	schorsch	2.18	extern void
166			dothresh( /* gaussian threshold filter */
167			int xcent,
168			int ycent,
169			int ccent,
170			int rcent
171			)
172	greg	2.7	{
173			double d;
174	greg	2.12	int r, y, offs;
175	greg	2.7	register int c, x;
176			register float *gscan;
177			/* compute ring sums */
178	schorsch	2.16	memset((char )ringsum, '\0', (orad+1)sizeof(float));
179			memset((char )ringwt, '\0', (orad+1)sizeof(short));
180	greg	2.7	for (r = -orad; r <= orad; r++) {
181			if (rcent+r < 0) continue;
182	greg	2.8	if (rcent+r >= nrows) break;
183	greg	2.7	gscan = greybar[(rcent+r)%obarsize];
184			for (c = -orad; c <= orad; c++) {
185	greg	2.12	offs = ccent+c < 0 ? ncols :
186			ccent+c >= ncols ? -ncols : 0;
187			if (offs && !wrapfilt)
188			continue;
189	greg	2.7	x = ringndx[cc + rr];
190			if (x < 0) continue;
191	greg	2.12	ringsum[x] += gscan[ccent+c+offs];
192	greg	2.7	ringwt[x]++;
193			}
194			}
195			/* filter each subpixel */
196			for (y = ycent+1-ybrad; y <= ycent+ybrad; y++) {
197			if (y < 0) continue;
198			if (y >= yres) break;
199	greg	2.13	d = y_r < 1.0 ? y_r*y - (rcent+.5) : (double)(y - ycent);
200	greg	2.7	if (d < -0.5) continue;
201			if (d >= 0.5) break;
202			for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
203	greg	2.12	offs = x < 0 ? xres : x >= xres ? -xres : 0;
204			if (offs && !wrapfilt)
205			continue;
206	greg	2.13	d = x_c < 1.0 ? x_c*x - (ccent+.5) : (double)(x - xcent);
207	greg	2.7	if (d < -0.5) continue;
208			if (d >= 0.5) break;
209	greg	2.12	sumans(x, y, rcent, ccent,
210			pickfilt((*ourbright)(scanin[y%barsize][x+offs])));
211	greg	2.7	}
212			}
213			}
214
215
216	schorsch	2.18	static double
217			pickfilt( /* find filter multiplier for p0 */
218			double p0
219			)
220	greg	2.7	{
221			double m = 1.0;
222	greg	2.10	double t, num, denom, avg, wsum;
223			double mlimit[2];
224	greg	2.14	int ilimit = 4.0/TEPS;
225	greg	2.7	register int i;
226			/* iterative search for m */
227	greg	2.11	mlimit[0] = 1.0; mlimit[1] = orad/rad/CHECKRAD;
228	greg	2.7	do {
229			/* compute grey weighted average */
230	greg	2.14	i = RSCACHECKRADrad*m + .5;
231	greg	2.8	if (i > orad) i = orad;
232	greg	2.7	avg = wsum = 0.0;
233			while (i--) {
234	greg	2.8	t = (i+.5)/(m*rad);
235	greg	2.7	t = lookgauss(t*t);
236			avg += t*ringsum[i];
237			wsum += t*ringwt[i];
238			}
239	greg	2.10	if (avg < 1e-20) /* zero inclusive average */
240			return(1.0);
241	greg	2.7	avg /= wsum;
242			/* check threshold */
243	greg	2.10	denom = m*m/gausstable[0] - p0/avg;
244	greg	2.11	if (denom <= FTINY) { /* zero exclusive average */
245			if (m >= mlimit[1]-REPS)
246			break;
247			m = mlimit[1];
248			continue;
249			}
250	greg	2.10	num = p0/avg - 1.0;
251			if (num < 0.0) num = -num;
252			t = num/denom;
253			if (t <= thresh) {
254			if (m <= mlimit[0]+REPS \|\| (thresh-t)/thresh <= TEPS)
255			break;
256			} else {
257			if (m >= mlimit[1]-REPS \|\| (t-thresh)/thresh <= TEPS)
258			break;
259			}
260			t = m; /* remember current m */
261	greg	2.7	/* next guesstimate */
262	greg	2.10	m = sqrt(gausstable[0]*(num/thresh + p0/avg));
263			if (m < t) { /* bound it */
264			if (m <= mlimit[0]+FTINY)
265			m = 0.5*(mlimit[0] + t);
266			mlimit[1] = t;
267			} else {
268			if (m >= mlimit[1]-FTINY)
269			m = 0.5*(mlimit[1] + t);
270			mlimit[0] = t;
271			}
272	greg	2.7	} while (--ilimit > 0);
273			return(m);
274			}
275
276
277	schorsch	2.18	static void
278			sumans( /* sum input pixel to output */
279			int px,
280			int py,
281			int rcent,
282			int ccent,
283			double m
284			)
285	greg	2.7	{
286	greg	2.14	double dy2, dx;
287	greg	2.7	COLOR pval, ctmp;
288	greg	2.12	int ksiz, r, offs;
289	greg	2.7	double pc, pr, norm;
290			register int i, c;
291			register COLOR *scan;
292	greg	2.12	/*
293			* This normalization method fails at the picture borders because
294			* a different number of input pixels contribute there.
295			*/
296			scan = scanin[py%barsize] + (px < 0 ? xres : px >= xres ? -xres : 0);
297			copycolor(pval, scan[px]);
298	greg	2.7	pc = x_c*(px+.5);
299			pr = y_r*(py+.5);
300			ksiz = CHECKRADmrad + 1;
301			if (ksiz > orad) ksiz = orad;
302			/* compute normalization */
303			norm = 0.0;
304			i = 0;
305			for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
306			if (r < 0) continue;
307			if (r >= nrows) break;
308	greg	2.14	dy2 = (pr - (r+.5))/(m*rad);
309			dy2 *= dy2;
310	greg	2.7	for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
311	greg	2.12	if (!wrapfilt) {
312			if (c < 0) continue;
313			if (c >= ncols) break;
314			}
315	greg	2.7	dx = (pc - (c+.5))/(m*rad);
316	greg	2.14	norm += warr[i++] = lookgauss(dx*dx + dy2);
317	greg	2.7	}
318			}
319			norm = 1.0/norm;
320			if (x_c < 1.0) norm *= x_c;
321			if (y_r < 1.0) norm *= y_r;
322			/* sum pixels */
323			i = 0;
324			for (r = rcent-ksiz; r <= rcent+ksiz; r++) {
325			if (r < 0) continue;
326			if (r >= nrows) break;
327			scan = scoutbar[r%obarsize];
328			for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
329	greg	2.12	offs = c < 0 ? ncols : c >= ncols ? -ncols : 0;
330			if (offs && !wrapfilt)
331			continue;
332	greg	2.7	copycolor(ctmp, pval);
333			dx = norm*warr[i++];
334			scalecolor(ctmp, dx);
335	greg	2.12	addcolor(scan[c+offs], ctmp);
336	greg	2.7	}
337			}
338	greg	1.1	}