src/px/pf3.c

#ifndef lint
static const char RCSid[] = "$Id: pf3.c,v 2.21 2023/12/08 17:56:26 greg Exp $";
#endif
/*
 *  pf3.c - routines for gaussian and box filtering
 *
 *     8/13/86
 */

#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)(20.*(x)+.5)]

static double pickfilt(double  p0);
static void sumans(int  px, int  py, int  rcent, int  ccent, double  m);


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

        gtabsiz = 444*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.05 <= d)
                        gausstable[x] = gausstable[0];
                else
                        gausstable[x] = exp(-x*0.05);
        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);
}


void
dobox(                  /* simple box filter */
        SCOLOR  csum,
        int  xcent,
        int  ycent,
        int  c,
        int  r
)
{
        int  wsum;
        double  d;
        int  y;
        int  x, offs;
        COLORV  *scan, *scp;
        
        wsum = 0;
        scolorblack(csum);
        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 != 0) & !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++;
                        scp = scan + (x+offs)*NCSAMP;
                        saddscolor(csum, scp);
                }
        }
        if (wsum > 1) {
                d = 1.0/wsum;
                scalescolor(csum, d);
        }
}


void
dogauss(                /* gaussian filter */
        SCOLOR  csum,
        int  xcent,
        int  ycent,
        int  c,
        int  r
)
{
        double  dy, dx, weight, wsum;
        SCOLOR  ctmp;
        int  y;
        int  x, offs;
        COLORV  *scan;

        wsum = FTINY;
        scolorblack(csum);
        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 != 0) & !wrapfilt)
                                continue;
                        dx = (x_c*(x+.5) - (c+.5))/rad;
                        weight = lookgauss(dx*dx + dy*dy);
                        wsum += weight;
                        copyscolor(ctmp, scan+(x+offs)*NCSAMP);
                        scalescolor(ctmp, weight);
                        saddscolor(csum, ctmp);
                }
        }
        weight = 1.0/wsum;
        scalescolor(csum, weight);
}


void
dothresh(       /* gaussian threshold filter */
        int  xcent,
        int  ycent,
        int  ccent,
        int  rcent
)
{
        double  d;
        int  r, y, offs;
        int  c, x;
        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 != 0) & !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 != 0) & !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)*NCSAMP)));
                }
        }
}


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