src/px/pf3.c

/* Copyright (c) 1992 Regents of the University of California */

#ifndef lint
static char SCCSid[] = "$SunId$ LBL";
#endif

/*
 *  pf3.c - routines for gaussian and box filtering
 *
 *     8/13/86
 */

#include  <stdio.h>

#include  <math.h>

#include  "color.h"

#define  FTINY          1e-7

extern double  rad;             /* output pixel radius for filtering */

extern int  nrows;              /* number of rows for output */
extern int  ncols;              /* number of columns for output */

extern int  boxfilt;            /* do box filtering? */

extern int  xres, yres;         /* resolution of input */

extern double  x_c, y_r;        /* conversion factors */

extern int  xrad;               /* x window size */
extern int  yrad;               /* y window size */

extern int  barsize;            /* size of input scan bar */
extern COLOR  **scanin;         /* input scan bar */
extern COLOR  *scanout;         /* output scan line */

extern char  *progname;

float  *gausstable;             /* gauss lookup table */

#define  lookgauss(x)           gausstable[(int)(-10.*(x)+.5)]


initmask()                      /* initialize gaussian lookup table */
{
        extern char  *malloc();
        double  d;
        register int  x;

        gausstable = (float *)malloc(100*sizeof(float));
        if (gausstable == NULL) {
                fprintf(stderr, "%s: out of memory in initmask\n", progname);
                quit(1);
        }
        d = x_c*y_r*0.25/rad/rad;
        gausstable[0] = exp(-d)/sqrt(d);
        for (x = 1; x < 100; x++)
                if ((gausstable[x] = exp(-x*0.1)/sqrt(x*0.1)) > gausstable[0])
                        gausstable[x] = gausstable[0];
}


dobox(csum, xcent, ycent, c, r)                 /* simple box filter */
COLOR  csum;
int  xcent, ycent;
int  c, r;
{
        static int  wsum;
        static double  d;
        static int  y;
        register int  x;
        register COLOR  *scan;

        wsum = 0;
        setcolor(csum, 0.0, 0.0, 0.0);
        for (y = ycent+1-yrad; y <= ycent+yrad; y++) {
                if (y < 0) continue;
                if (y >= yres) break;
                d = y_r < 1.0 ? y_r*y - r : (double)(y - ycent);
                if (d > 0.5+FTINY || d < -0.5-FTINY)
                        continue;
                scan = scanin[y%barsize];
                for (x = xcent+1-xrad; x <= xcent+xrad; x++) {
                        if (x < 0) continue;
                        if (x >= xres) break;
                        d = x_c < 1.0 ? x_c*x - c : (double)(x - xcent);
                        if (d > 0.5+FTINY || d < -0.5-FTINY)
                                continue;
                        wsum++;
                        addcolor(csum, scan[x]);
                }
        }
        if (wsum > 1)
                scalecolor(csum, 1.0/wsum);
}


dogauss(csum, xcent, ycent, c, r)               /* gaussian filter */
COLOR  csum;
int  xcent, ycent;
int  c, r;
{
        static double  dy, dx, weight, wsum;
        static COLOR  ctmp;
        static int  y;
        register int  x;
        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++) {
                        if (x < 0) continue;
                        if (x >= xres) break;
                        dx = (x_c*(x+.5) - (c+.5))/rad;
                        weight = lookgauss(-(dx*dx + dy*dy));
                        wsum += weight;
                        copycolor(ctmp, scan[x]);
                        scalecolor(ctmp, weight);
                        addcolor(csum, ctmp);
                }
        }
        scalecolor(csum, 1.0/wsum);
}
Revision:	2.5
Committed:	Fri Jun 18 13:36:20 1993 UTC (32 years, 4 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.4:	+11 -7 lines
Log Message:	fixed long-standing bug in normalization of Gaussian filter
#	User	Rev	Content
1	greg	2.2	/* Copyright (c) 1992 Regents of the University of California */
2	greg	1.1
3			#ifndef lint
4			static char SCCSid[] = "$SunId$ LBL";
5			#endif
6
7			/*
8			* pf3.c - routines for gaussian and box filtering
9			*
10			* 8/13/86
11			*/
12
13			#include <stdio.h>
14
15	greg	2.3	#include <math.h>
16
17	greg	1.1	#include "color.h"
18
19	greg	2.2	#define FTINY 1e-7
20	greg	1.1
21			extern double rad; /* output pixel radius for filtering */
22
23			extern int nrows; /* number of rows for output */
24			extern int ncols; /* number of columns for output */
25
26			extern int boxfilt; /* do box filtering? */
27
28			extern int xres, yres; /* resolution of input */
29
30			extern double x_c, y_r; /* conversion factors */
31
32			extern int xrad; /* x window size */
33			extern int yrad; /* y window size */
34
35			extern int barsize; /* size of input scan bar */
36			extern COLOR *scanin; / input scan bar */
37			extern COLOR scanout; / output scan line */
38
39			extern char *progname;
40
41	greg	2.5	float gausstable; / gauss lookup table */
42	greg	1.1
43	greg	2.5	#define lookgauss(x) gausstable[(int)(-10.*(x)+.5)]
44	greg	1.1
45
46			initmask() /* initialize gaussian lookup table */
47			{
48			extern char *malloc();
49	greg	2.5	double d;
50	greg	1.1	register int x;
51
52	greg	2.5	gausstable = (float )malloc(100sizeof(float));
53			if (gausstable == NULL) {
54	greg	1.1	fprintf(stderr, "%s: out of memory in initmask\n", progname);
55			quit(1);
56			}
57	greg	2.5	d = x_cy_r0.25/rad/rad;
58			gausstable[0] = exp(-d)/sqrt(d);
59			for (x = 1; x < 100; x++)
60			if ((gausstable[x] = exp(-x0.1)/sqrt(x0.1)) > gausstable[0])
61			gausstable[x] = gausstable[0];
62	greg	1.1	}
63
64
65			dobox(csum, xcent, ycent, c, r) /* simple box filter */
66			COLOR csum;
67			int xcent, ycent;
68			int c, r;
69			{
70			static int wsum;
71			static double d;
72			static int y;
73			register int x;
74	greg	2.2	register COLOR *scan;
75	greg	1.1
76			wsum = 0;
77			setcolor(csum, 0.0, 0.0, 0.0);
78			for (y = ycent+1-yrad; y <= ycent+yrad; y++) {
79	greg	2.4	if (y < 0) continue;
80			if (y >= yres) break;
81			d = y_r < 1.0 ? y_r*y - r : (double)(y - ycent);
82	greg	1.1	if (d > 0.5+FTINY \|\| d < -0.5-FTINY)
83			continue;
84			scan = scanin[y%barsize];
85			for (x = xcent+1-xrad; x <= xcent+xrad; x++) {
86	greg	2.4	if (x < 0) continue;
87			if (x >= xres) break;
88			d = x_c < 1.0 ? x_c*x - c : (double)(x - xcent);
89	greg	1.1	if (d > 0.5+FTINY \|\| d < -0.5-FTINY)
90			continue;
91			wsum++;
92			addcolor(csum, scan[x]);
93			}
94			}
95			if (wsum > 1)
96			scalecolor(csum, 1.0/wsum);
97			}
98
99
100			dogauss(csum, xcent, ycent, c, r) /* gaussian filter */
101			COLOR csum;
102			int xcent, ycent;
103			int c, r;
104			{
105			static double dy, dx, weight, wsum;
106			static COLOR ctmp;
107			static int y;
108			register int x;
109	greg	2.2	register COLOR *scan;
110	greg	1.1
111			wsum = FTINY;
112			setcolor(csum, 0.0, 0.0, 0.0);
113	greg	1.2	for (y = ycent-yrad; y <= ycent+yrad; y++) {
114	greg	2.4	if (y < 0) continue;
115			if (y >= yres) break;
116	greg	1.2	dy = (y_r*(y+.5) - (r+.5))/rad;
117	greg	1.1	scan = scanin[y%barsize];
118	greg	1.2	for (x = xcent-xrad; x <= xcent+xrad; x++) {
119	greg	2.4	if (x < 0) continue;
120			if (x >= xres) break;
121	greg	1.2	dx = (x_c*(x+.5) - (c+.5))/rad;
122	greg	2.5	weight = lookgauss(-(dxdx + dydy));
123	greg	1.1	wsum += weight;
124			copycolor(ctmp, scan[x]);
125			scalecolor(ctmp, weight);
126			addcolor(csum, ctmp);
127			}
128			}
129			scalecolor(csum, 1.0/wsum);
130			}