src/hd/rhpict2.c

/* Copyright (c) 1999 Silicon Graphics, Inc. */

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

/*
 * Rendering routines for rhpict.
 */

#include "holo.h"
#include "view.h"
#include "random.h"

#ifndef DEPS
#define DEPS            0.02    /* depth epsilon */
#endif
#ifndef PEPS
#define PEPS            0.04    /* pixel value epsilon */
#endif
#ifndef MAXRAD
#define MAXRAD          64      /* maximum kernel radius */
#endif
#ifndef NNEIGH
#define NNEIGH          7       /* find this many neighbors */
#endif

#define NINF            16382

#define MAXRAD2         (MAXRAD*MAXRAD+1)

#define G0NORM          0.286   /* ground zero normalization (1/x integral) */

#ifndef FL4OP
#define FL4OP(f,i,op)   ((f)[(i)>>5] op (1L<<((i)&0x1f)))
#define CHK4(f,i)       FL4OP(f,i,&)
#define SET4(f,i)       FL4OP(f,i,|=)
#define CLR4(f,i)       FL4OP(f,i,&=~)
#define TGL4(f,i)       FL4OP(f,i,^=)
#define FL4NELS(n)      (((n)+0x1f)>>5)
#define CLR4ALL(f,n)    bzero((char *)(f),FL4NELS(n)*sizeof(int4))
#endif

static int4     *pixFlags;      /* pixel occupancy flags */
static float    pixWeight[MAXRAD2];     /* pixel weighting function */
static short    isqrttab[MAXRAD2];      /* integer square root table */

#define isqrt(i2)       ((int)isqrttab[(int)(i2)])

extern VIEW     myview;         /* current output view */
extern COLOR    *mypixel;       /* pixels being rendered */
extern float    *myweight;      /* weights (used to compute final pixels) */
extern float    *mydepth;       /* depth values (visibility culling) */
extern int      hres, vres;     /* current horizontal and vertical res. */


pixBeam(bp, hb)                 /* render a particular beam */
BEAM    *bp;
register HDBEAMI        *hb;
{
        GCOORD  gc[2];
        register RAYVAL *rv;
        FVECT   rorg, rdir, wp, ip;
        double  d, prox;
        COLOR   col;
        int     n;
        register int4   p;

        if (!hdbcoord(gc, hb->h, hb->b))
                error(CONSISTENCY, "bad beam in render_beam");
        for (n = bp->nrm, rv = hdbray(bp); n--; rv++) {
                                                /* reproject each sample */
                hdray(rorg, rdir, hb->h, gc, rv->r);
                if (rv->d < DCINF) {
                        d = hddepth(hb->h, rv->d);
                        VSUM(wp, rorg, rdir, d);
                        VSUB(ip, wp, myview.vp);
                        d = DOT(ip,rdir);
                        prox = d*d/DOT(ip,ip);  /* cos(diff_angle)^32 */
                        prox *= prox; prox *= prox; prox *= prox; prox *= prox;
                } else {
                        if (myview.type == VT_PAR || myview.vaft > FTINY)
                                continue;               /* inf. off view */
                        VSUM(wp, myview.vp, rdir, FHUGE);
                        prox = 1.;
                }
                viewloc(ip, &myview, wp);       /* frustum clipping */
                if (ip[2] < 0.)
                        continue;
                if (ip[0] < 0. || ip[0] >= 1.)
                        continue;
                if (ip[1] < 0. || ip[1] >= 1.)
                        continue;
                if (myview.vaft > FTINY && ip[2] > myview.vaft - myview.vfore)
                        continue;               /* not exact for VT_PER */
                p = (int)(ip[1]*vres)*hres + (int)(ip[0]*hres);
                if (mydepth[p] > FTINY) {       /* check depth */
                        if (ip[2] > mydepth[p]*(1.+DEPS))
                                continue;
                        if (ip[2] < mydepth[p]*(1.-DEPS)) {
                                setcolor(mypixel[p], 0., 0., 0.);
                                myweight[p] = 0.;
                        }
                }
                colr_color(col, rv->v);
                scalecolor(col, prox);
                addcolor(mypixel[p], col);
                myweight[p] += prox;
                mydepth[p] = ip[2];
        }
}


int
kill_occl(h, v, nl, n)          /* check for occlusion errors */
int     h, v;
register short  nl[NNEIGH][2];
int     n;
{
        short   forequad[2][2];
        int     d;
        register int4   i, p;

        if (n <= 0) {
#ifdef DEBUG
                error(WARNING, "neighborless sample in kill_occl");
#endif
                return(1);
        }
        p = v*hres + h;
        forequad[0][0] = forequad[0][1] = forequad[1][0] = forequad[1][1] = 0;
        for (i = n; i--; ) {
                d = (h-nl[i][0])*(h-nl[i][0]) + (v-nl[i][1])*(v-nl[i][1]);
                d = isqrt(d);
                if (mydepth[nl[i][1]*hres+nl[i][0]]*(1.+DEPS*d) < mydepth[p])
                        forequad[nl[i][0]<h][nl[i][1]<v] = 1;
        }
        if (forequad[0][0]+forequad[0][1]+forequad[1][0]+forequad[1][1] > 2) {
                setcolor(mypixel[p], 0., 0., 0.);
                myweight[p] = 0.;       /* occupancy reset afterwards */
        }
        return(1);
}


int
smooth_samp(h, v, nl, n)        /* grow sample point smoothly */
int     h, v;
register short  nl[NNEIGH][2];
int     n;
{
        int     dis[NNEIGH], ndis;
        COLOR   mykern[MAXRAD2];
        int4    maxr2;
        double  d;
        register int4   p, r2;
        int     i, r, maxr, h2, v2;

        if (n <= 0)
                return(1);
        p = v*hres + h;                         /* build kernel values */
        maxr2 = (h-nl[n-1][0])*(h-nl[n-1][0]) + (v-nl[n-1][1])*(v-nl[n-1][1]);
        DCHECK(maxr2>=MAXRAD2, CONSISTENCY, "out of range neighbor");
        maxr = isqrt(maxr2);
        for (v2 = 1; v2 <= maxr; v2++)
                for (h2 = 0; h2 <= v2; h2++) {
                        r2 = h2*h2 + v2*v2;
                        if (r2 > maxr2) break;
                        copycolor(mykern[r2], mypixel[p]);
                        scalecolor(mykern[r2], pixWeight[r2]);
                }
        ndis = 0;                               /* find discontinuities */
        for (i = n; i--; ) {
                r2 = (h-nl[i][0])*(h-nl[i][0]) + (v-nl[i][1])*(v-nl[i][1]);
                r = isqrt(r2);
                d = mydepth[nl[i][1]*hres+nl[i][0]] / mydepth[p];
                d = d>=1. ? d-1. : 1.-d;
                if (d > r*DEPS || bigdiff(mypixel[p],
                                mypixel[nl[i][1]*hres+nl[i][0]], r*PEPS))
                        dis[ndis++] = i;
        }
                                                /* stamp out that kernel */
        for (v2 = v-maxr; v2 <= v+maxr; v2++) {
                if (v2 < 0) v2 = 0;
                else if (v2 >= vres) break;
                for (h2 = h-maxr; h2 <= h+maxr; h2++) {
                        if (h2 < 0) h2 = 0;
                        else if (h2 >= hres) break;
                        r2 = (h2-h)*(h2-h) + (v2-v)*(v2-v);
                        if (r2 > maxr2) continue;
                        if (CHK4(pixFlags, v2*hres+h2))
                                continue;       /* occupied */
                        for (i = ndis; i--; ) {
                                r = (h2-nl[dis[i]][0])*(h2-nl[dis[i]][0]) +
                                        (v2-nl[dis[i]][1])*(v2-nl[dis[i]][1]);
                                if (r < r2) break;
                        }
                        if (i >= 0) continue;   /* outside edge */
                        addcolor(mypixel[v2*hres+h2], mykern[r2]);
                        myweight[v2*hres+h2] += pixWeight[r2] * myweight[p];
                }
        }
        return(1);
}


int
random_samp(h, v, nl, n, rf)    /* grow sample point noisily */
int     h, v;
register short  nl[NNEIGH][2];
int     n;
double  *rf;
{
        double  tv = *rf * (1.-1./NNEIGH);
        int4    maxr2;
        register int4   p, r2;
        register int    i;
        int     maxr, h2, v2;
        COLOR   ctmp;

        if (n <= 0)
                return(1);
        p = v*hres + h;                         /* compute kernel radius */
        maxr2 = (h-nl[n-1][0])*(h-nl[n-1][0]) + (v-nl[n-1][1])*(v-nl[n-1][1]);
        DCHECK(maxr2>=MAXRAD2, CONSISTENCY, "out of range neighbor");
        maxr = isqrt(maxr2);
                                                /* sample kernel */
        for (v2 = v-maxr; v2 <= v+maxr; v2++) {
                if (v2 < 0) v2 = 0;
                else if (v2 >= vres) break;
                for (h2 = h-maxr; h2 <= h+maxr; h2++) {
                        if (h2 < 0) h2 = 0;
                        else if (h2 >= hres) break;
                        r2 = (h2-h)*(h2-h) + (v2-v)*(v2-v);
                        if (r2 > maxr2) continue;
                        if (CHK4(pixFlags, v2*hres+h2))
                                continue;       /* occupied */
                        if (frandom() < tv) {   /* use neighbor instead? */
                                i = random() % n;
                                r2 = nl[i][1]*hres + nl[i][0];
                                copycolor(ctmp, mypixel[r2]);
                                r2 = (h2-nl[i][0])*(h2-nl[i][0]) +
                                                (v2-nl[i][1])*(v2-nl[i][1]);
                                if (r2 < MAXRAD2) {
                                        scalecolor(ctmp, pixWeight[r2]);
                                        addcolor(mypixel[v2*hres+h2], ctmp);
                                        myweight[v2*hres+h2] += pixWeight[r2] *
                                              myweight[nl[i][1]*hres+nl[i][0]];
                                        continue;
                                }
                        }
                                                /* use central sample */
                        copycolor(ctmp, mypixel[p]);
                        scalecolor(ctmp, pixWeight[r2]);
                        addcolor(mypixel[v2*hres+h2], ctmp);
                        myweight[v2*hres+h2] += pixWeight[r2] * myweight[p];
                }
        }
        return(1);
}


pixFinish(ransamp)              /* done with beams -- compute pixel values */
double  ransamp;
{
        if (pixWeight[0] <= FTINY)
                init_wfunc();           /* initialize weighting function */
        reset_flags();                  /* set occupancy flags */
        meet_neighbors(kill_occl,NULL); /* eliminate occlusion errors */
        reset_flags();                  /* reset occupancy flags */
        if (ransamp >= 0.)              /* spread samples over image */
                meet_neighbors(random_samp,&ransamp);
        else
                meet_neighbors(smooth_samp,NULL);
        free((char *)pixFlags);         /* free pixel flags */
        pixFlags = NULL;
}


reset_flags()                   /* allocate/set/reset occupancy flags */
{
        register int4   p;

        if (pixFlags == NULL) {
                pixFlags = (int4 *)calloc(FL4NELS(hres*vres), sizeof(int4));
                CHECK(pixFlags==NULL, SYSTEM, "out of memory in reset_flags");
        } else
                CLR4ALL(pixFlags, hres*vres);
        for (p = hres*vres; p--; )
                if (myweight[p] > FTINY)
                        SET4(pixFlags, p);
}


init_wfunc()                    /* initialize weighting function */
{
        register int4   r2;
        register double d;

        for (r2 = MAXRAD2; --r2; ) {
                d = sqrt((double)r2);
                pixWeight[r2] = G0NORM/d;
                isqrttab[r2] = d + 0.99;
        }
        pixWeight[0] = 1.;
        isqrttab[0] = 0;
}


int
findneigh(nl, h, v, rnl)        /* find NNEIGH neighbors for pixel */
short   nl[NNEIGH][2];
int     h, v;
register short  (*rnl)[NNEIGH];
{
        int     nn = 0;
        int4    d, nd[NNEIGH];
        int     n, hoff;
        register int    h2, n2;

        nd[NNEIGH-1] = MAXRAD2;
        for (hoff = 1; hoff < hres; hoff = (hoff<0) - hoff) {
                h2 = h + hoff;
                if (h2 < 0 | h2 >= hres)
                        continue;
                if ((h2-h)*(h2-h) >= nd[NNEIGH-1])
                        break;
                for (n = 0; n < NNEIGH && rnl[h2][n] < NINF; n++) {
                        d = (h2-h)*(h2-h) + (v-rnl[h2][n])*(v-rnl[h2][n]);
                        if (d >= nd[NNEIGH-1])
                                continue;
                        if (nn < NNEIGH)        /* insert neighbor */
                                nn++;
                        for (n2 = nn; n2--; )   {
                                if (!n2 || d >= nd[n2-1]) {
                                        nd[n2] = d;
                                        nl[n2][0] = h2;
                                        nl[n2][1] = rnl[h2][n];
                                        break;
                                }
                                nd[n2] = nd[n2-1];
                                nl[n2][0] = nl[n2-1][0];
                                nl[n2][1] = nl[n2-1][1];
                        }
                }
        }
        return(nn);
}


meet_neighbors(nf, dp)          /* run through samples and their neighbors */
int     (*nf)();
char    *dp;
{
        short   ln[NNEIGH][2];
        int     h, v, n, v2;
        register short  (*rnl)[NNEIGH];
                                        /* initialize bottom row list */
        rnl = (short (*)[NNEIGH])malloc(NNEIGH*sizeof(short)*hres);
        CHECK(rnl==NULL, SYSTEM, "out of memory in meet_neighbors");
        for (h = 0; h < hres; h++) {
                for (n = v = 0; v < vres; v++)
                        if (CHK4(pixFlags, v*hres+h)) {
                                rnl[h][n++] = v;
                                if (n >= NNEIGH)
                                        break;
                        }
                while (n < NNEIGH)
                        rnl[h][n++] = NINF;
        }
        v = 0;                          /* do each row */
        for ( ; ; ) {
                for (h = 0; h < hres; h++) {
                        if (!CHK4(pixFlags, v*hres+h))
                                continue;       /* no one home */
                        n = findneigh(ln, h, v, rnl);
                        (*nf)(h, v, ln, n, dp); /* call on neighbors */
                }
                if (++v >= vres)                /* reinitialize row list */
                        break;
                for (h = 0; h < hres; h++)
                        for (v2 = rnl[h][NNEIGH-1]+1; v2 < vres; v2++) {
                                if (v2 - v > v - rnl[h][0])
                                        break;          /* not close enough */
                                if (CHK4(pixFlags, v2*hres+h)) {
                                        for (n = 0; n < NNEIGH-1; n++)
                                                rnl[h][n] = rnl[h][n+1];
                                        rnl[h][NNEIGH-1] = v2;
                                }
                        }
        }
        free((char *)rnl);              /* free row list */
}
Revision:	3.8
Committed:	Tue Mar 9 15:08:22 1999 UTC (25 years, 1 month ago) by gwlarson
Content type:	text/plain
Branch:	MAIN
Changes since 3.7:	+12 -11 lines
Log Message:	final cut on random resampling algorithm -- fairer choice algorithm
#	User	Rev	Content
1	gwlarson	3.1	/* Copyright (c) 1999 Silicon Graphics, Inc. */
2
3			#ifndef lint
4			static char SCCSid[] = "$SunId$ SGI";
5			#endif
6
7			/*
8			* Rendering routines for rhpict.
9			*/
10
11			#include "holo.h"
12			#include "view.h"
13	gwlarson	3.6	#include "random.h"
14	gwlarson	3.1
15	gwlarson	3.2	#ifndef DEPS
16	gwlarson	3.3	#define DEPS 0.02 /* depth epsilon */
17	gwlarson	3.2	#endif
18	gwlarson	3.3	#ifndef PEPS
19			#define PEPS 0.04 /* pixel value epsilon */
20			#endif
21	gwlarson	3.2	#ifndef MAXRAD
22			#define MAXRAD 64 /* maximum kernel radius */
23			#endif
24			#ifndef NNEIGH
25			#define NNEIGH 7 /* find this many neighbors */
26			#endif
27
28			#define NINF 16382
29
30			#define MAXRAD2 (MAXRAD*MAXRAD+1)
31
32			#define G0NORM 0.286 /* ground zero normalization (1/x integral) */
33
34			#ifndef FL4OP
35			#define FL4OP(f,i,op) ((f)[(i)>>5] op (1L<<((i)&0x1f)))
36			#define CHK4(f,i) FL4OP(f,i,&)
37			#define SET4(f,i) FL4OP(f,i,\|=)
38			#define CLR4(f,i) FL4OP(f,i,&=~)
39			#define TGL4(f,i) FL4OP(f,i,^=)
40			#define FL4NELS(n) (((n)+0x1f)>>5)
41			#define CLR4ALL(f,n) bzero((char )(f),FL4NELS(n)sizeof(int4))
42			#endif
43
44			static int4 pixFlags; / pixel occupancy flags */
45			static float pixWeight[MAXRAD2]; /* pixel weighting function */
46	gwlarson	3.3	static short isqrttab[MAXRAD2]; /* integer square root table */
47	gwlarson	3.2
48	gwlarson	3.3	#define isqrt(i2) ((int)isqrttab[(int)(i2)])
49
50	gwlarson	3.1	extern VIEW myview; /* current output view */
51			extern COLOR mypixel; / pixels being rendered */
52			extern float myweight; / weights (used to compute final pixels) */
53	gwlarson	3.2	extern float mydepth; / depth values (visibility culling) */
54	gwlarson	3.1	extern int hres, vres; /* current horizontal and vertical res. */
55
56
57	gwlarson	3.2	pixBeam(bp, hb) /* render a particular beam */
58	gwlarson	3.1	BEAM *bp;
59			register HDBEAMI *hb;
60			{
61			GCOORD gc[2];
62			register RAYVAL *rv;
63			FVECT rorg, rdir, wp, ip;
64			double d, prox;
65			COLOR col;
66	gwlarson	3.2	int n;
67			register int4 p;
68	gwlarson	3.1
69			if (!hdbcoord(gc, hb->h, hb->b))
70			error(CONSISTENCY, "bad beam in render_beam");
71			for (n = bp->nrm, rv = hdbray(bp); n--; rv++) {
72			/* reproject each sample */
73			hdray(rorg, rdir, hb->h, gc, rv->r);
74			if (rv->d < DCINF) {
75			d = hddepth(hb->h, rv->d);
76			VSUM(wp, rorg, rdir, d);
77			VSUB(ip, wp, myview.vp);
78			d = DOT(ip,rdir);
79	gwlarson	3.2	prox = dd/DOT(ip,ip); / cos(diff_angle)^32 */
80			prox = prox; prox = prox; prox = prox; prox = prox;
81	gwlarson	3.1	} else {
82			if (myview.type == VT_PAR \|\| myview.vaft > FTINY)
83			continue; /* inf. off view */
84			VSUM(wp, myview.vp, rdir, FHUGE);
85	gwlarson	3.2	prox = 1.;
86	gwlarson	3.1	}
87			viewloc(ip, &myview, wp); /* frustum clipping */
88			if (ip[2] < 0.)
89			continue;
90			if (ip[0] < 0. \|\| ip[0] >= 1.)
91			continue;
92			if (ip[1] < 0. \|\| ip[1] >= 1.)
93			continue;
94			if (myview.vaft > FTINY && ip[2] > myview.vaft - myview.vfore)
95	gwlarson	3.2	continue; /* not exact for VT_PER */
96	gwlarson	3.1	p = (int)(ip[1]vres)hres + (int)(ip[0]*hres);
97	gwlarson	3.2	if (mydepth[p] > FTINY) { /* check depth */
98			if (ip[2] > mydepth[p]*(1.+DEPS))
99			continue;
100			if (ip[2] < mydepth[p]*(1.-DEPS)) {
101			setcolor(mypixel[p], 0., 0., 0.);
102			myweight[p] = 0.;
103			}
104			}
105	gwlarson	3.1	colr_color(col, rv->v);
106	gwlarson	3.2	scalecolor(col, prox);
107	gwlarson	3.1	addcolor(mypixel[p], col);
108	gwlarson	3.2	myweight[p] += prox;
109			mydepth[p] = ip[2];
110	gwlarson	3.1	}
111	gwlarson	3.2	}
112
113
114			int
115			kill_occl(h, v, nl, n) /* check for occlusion errors */
116			int h, v;
117	gwlarson	3.3	register short nl[NNEIGH][2];
118	gwlarson	3.2	int n;
119			{
120			short forequad[2][2];
121			int d;
122	gwlarson	3.3	register int4 i, p;
123	gwlarson	3.2
124	gwlarson	3.5	if (n <= 0) {
125			#ifdef DEBUG
126			error(WARNING, "neighborless sample in kill_occl");
127			#endif
128	gwlarson	3.2	return(1);
129	gwlarson	3.5	}
130	gwlarson	3.3	p = v*hres + h;
131	gwlarson	3.2	forequad[0][0] = forequad[0][1] = forequad[1][0] = forequad[1][1] = 0;
132			for (i = n; i--; ) {
133			d = (h-nl[i][0])(h-nl[i][0]) + (v-nl[i][1])(v-nl[i][1]);
134	gwlarson	3.3	d = isqrt(d);
135			if (mydepth[nl[i][1]hres+nl[i][0]](1.+DEPS*d) < mydepth[p])
136	gwlarson	3.2	forequad[nl[i][0]<h][nl[i][1]<v] = 1;
137			}
138	gwlarson	3.3	if (forequad[0][0]+forequad[0][1]+forequad[1][0]+forequad[1][1] > 2) {
139			setcolor(mypixel[p], 0., 0., 0.);
140			myweight[p] = 0.; /* occupancy reset afterwards */
141	gwlarson	3.2	}
142			return(1);
143			}
144
145
146			int
147	gwlarson	3.6	smooth_samp(h, v, nl, n) /* grow sample point smoothly */
148	gwlarson	3.2	int h, v;
149			register short nl[NNEIGH][2];
150			int n;
151			{
152	gwlarson	3.3	int dis[NNEIGH], ndis;
153	gwlarson	3.2	COLOR mykern[MAXRAD2];
154			int4 maxr2;
155	gwlarson	3.5	double d;
156	gwlarson	3.2	register int4 p, r2;
157	gwlarson	3.3	int i, r, maxr, h2, v2;
158	gwlarson	3.2
159			if (n <= 0)
160			return(1);
161			p = vhres + h; / build kernel values */
162			maxr2 = (h-nl[n-1][0])(h-nl[n-1][0]) + (v-nl[n-1][1])(v-nl[n-1][1]);
163			DCHECK(maxr2>=MAXRAD2, CONSISTENCY, "out of range neighbor");
164	gwlarson	3.3	maxr = isqrt(maxr2);
165			for (v2 = 1; v2 <= maxr; v2++)
166			for (h2 = 0; h2 <= v2; h2++) {
167			r2 = h2h2 + v2v2;
168			if (r2 > maxr2) break;
169			copycolor(mykern[r2], mypixel[p]);
170			scalecolor(mykern[r2], pixWeight[r2]);
171			}
172			ndis = 0; /* find discontinuities */
173			for (i = n; i--; ) {
174			r2 = (h-nl[i][0])(h-nl[i][0]) + (v-nl[i][1])(v-nl[i][1]);
175			r = isqrt(r2);
176			d = mydepth[nl[i][1]*hres+nl[i][0]] / mydepth[p];
177			d = d>=1. ? d-1. : 1.-d;
178			if (d > r*DEPS \|\| bigdiff(mypixel[p],
179			mypixel[nl[i][1]hres+nl[i][0]], rPEPS))
180			dis[ndis++] = i;
181	gwlarson	3.2	}
182	gwlarson	3.3	/* stamp out that kernel */
183	gwlarson	3.2	for (v2 = v-maxr; v2 <= v+maxr; v2++) {
184	gwlarson	3.3	if (v2 < 0) v2 = 0;
185			else if (v2 >= vres) break;
186	gwlarson	3.2	for (h2 = h-maxr; h2 <= h+maxr; h2++) {
187	gwlarson	3.3	if (h2 < 0) h2 = 0;
188			else if (h2 >= hres) break;
189			r2 = (h2-h)(h2-h) + (v2-v)(v2-v);
190	gwlarson	3.2	if (r2 > maxr2) continue;
191			if (CHK4(pixFlags, v2*hres+h2))
192			continue; /* occupied */
193	gwlarson	3.3	for (i = ndis; i--; ) {
194			r = (h2-nl[dis[i]][0])*(h2-nl[dis[i]][0]) +
195			(v2-nl[dis[i]][1])*(v2-nl[dis[i]][1]);
196			if (r < r2) break;
197			}
198			if (i >= 0) continue; /* outside edge */
199	gwlarson	3.2	addcolor(mypixel[v2*hres+h2], mykern[r2]);
200	gwlarson	3.5	myweight[v2hres+h2] += pixWeight[r2] myweight[p];
201	gwlarson	3.2	}
202			}
203			return(1);
204			}
205
206
207	gwlarson	3.6	int
208	gwlarson	3.7	random_samp(h, v, nl, n, rf) /* grow sample point noisily */
209	gwlarson	3.6	int h, v;
210			register short nl[NNEIGH][2];
211			int n;
212	gwlarson	3.7	double *rf;
213	gwlarson	3.2	{
214	gwlarson	3.8	double tv = rf (1.-1./NNEIGH);
215	gwlarson	3.6	int4 maxr2;
216			register int4 p, r2;
217			register int i;
218			int maxr, h2, v2;
219			COLOR ctmp;
220
221			if (n <= 0)
222			return(1);
223			p = vhres + h; / compute kernel radius */
224			maxr2 = (h-nl[n-1][0])(h-nl[n-1][0]) + (v-nl[n-1][1])(v-nl[n-1][1]);
225			DCHECK(maxr2>=MAXRAD2, CONSISTENCY, "out of range neighbor");
226			maxr = isqrt(maxr2);
227			/* sample kernel */
228			for (v2 = v-maxr; v2 <= v+maxr; v2++) {
229			if (v2 < 0) v2 = 0;
230			else if (v2 >= vres) break;
231			for (h2 = h-maxr; h2 <= h+maxr; h2++) {
232			if (h2 < 0) h2 = 0;
233			else if (h2 >= hres) break;
234			r2 = (h2-h)(h2-h) + (v2-v)(v2-v);
235			if (r2 > maxr2) continue;
236			if (CHK4(pixFlags, v2*hres+h2))
237			continue; /* occupied */
238	gwlarson	3.8	if (frandom() < tv) { /* use neighbor instead? */
239	gwlarson	3.6	i = random() % n;
240			r2 = nl[i][1]*hres + nl[i][0];
241			copycolor(ctmp, mypixel[r2]);
242	gwlarson	3.8	r2 = (h2-nl[i][0])*(h2-nl[i][0]) +
243			(v2-nl[i][1])*(v2-nl[i][1]);
244			if (r2 < MAXRAD2) {
245			scalecolor(ctmp, pixWeight[r2]);
246			addcolor(mypixel[v2*hres+h2], ctmp);
247			myweight[v2hres+h2] += pixWeight[r2]
248			myweight[nl[i][1]*hres+nl[i][0]];
249			continue;
250			}
251	gwlarson	3.6	}
252	gwlarson	3.8	/* use central sample */
253	gwlarson	3.7	copycolor(ctmp, mypixel[p]);
254			scalecolor(ctmp, pixWeight[r2]);
255			addcolor(mypixel[v2*hres+h2], ctmp);
256			myweight[v2hres+h2] += pixWeight[r2] myweight[p];
257	gwlarson	3.6	}
258			}
259			return(1);
260			}
261
262
263			pixFinish(ransamp) /* done with beams -- compute pixel values */
264	gwlarson	3.7	double ransamp;
265	gwlarson	3.6	{
266	gwlarson	3.5	if (pixWeight[0] <= FTINY)
267			init_wfunc(); /* initialize weighting function */
268	gwlarson	3.2	reset_flags(); /* set occupancy flags */
269	gwlarson	3.7	meet_neighbors(kill_occl,NULL); /* eliminate occlusion errors */
270	gwlarson	3.2	reset_flags(); /* reset occupancy flags */
271	gwlarson	3.7	if (ransamp >= 0.) /* spread samples over image */
272			meet_neighbors(random_samp,&ransamp);
273	gwlarson	3.6	else
274	gwlarson	3.7	meet_neighbors(smooth_samp,NULL);
275	gwlarson	3.2	free((char )pixFlags); / free pixel flags */
276			pixFlags = NULL;
277			}
278
279
280			reset_flags() /* allocate/set/reset occupancy flags */
281			{
282	gwlarson	3.4	register int4 p;
283	gwlarson	3.2
284			if (pixFlags == NULL) {
285			pixFlags = (int4 )calloc(FL4NELS(hresvres), sizeof(int4));
286			CHECK(pixFlags==NULL, SYSTEM, "out of memory in reset_flags");
287			} else
288			CLR4ALL(pixFlags, hres*vres);
289			for (p = hres*vres; p--; )
290			if (myweight[p] > FTINY)
291			SET4(pixFlags, p);
292	gwlarson	3.4	}
293
294
295			init_wfunc() /* initialize weighting function */
296			{
297			register int4 r2;
298			register double d;
299
300	gwlarson	3.6	for (r2 = MAXRAD2; --r2; ) {
301			d = sqrt((double)r2);
302			pixWeight[r2] = G0NORM/d;
303			isqrttab[r2] = d + 0.99;
304			}
305	gwlarson	3.4	pixWeight[0] = 1.;
306			isqrttab[0] = 0;
307	gwlarson	3.2	}
308
309
310			int
311			findneigh(nl, h, v, rnl) /* find NNEIGH neighbors for pixel */
312			short nl[NNEIGH][2];
313			int h, v;
314			register short (*rnl)[NNEIGH];
315			{
316			int nn = 0;
317	gwlarson	3.3	int4 d, nd[NNEIGH];
318	gwlarson	3.2	int n, hoff;
319			register int h2, n2;
320
321	gwlarson	3.3	nd[NNEIGH-1] = MAXRAD2;
322	gwlarson	3.2	for (hoff = 1; hoff < hres; hoff = (hoff<0) - hoff) {
323			h2 = h + hoff;
324			if (h2 < 0 \| h2 >= hres)
325			continue;
326	gwlarson	3.3	if ((h2-h)*(h2-h) >= nd[NNEIGH-1])
327	gwlarson	3.2	break;
328			for (n = 0; n < NNEIGH && rnl[h2][n] < NINF; n++) {
329			d = (h2-h)(h2-h) + (v-rnl[h2][n])(v-rnl[h2][n]);
330	gwlarson	3.3	if (d >= nd[NNEIGH-1])
331	gwlarson	3.2	continue;
332	gwlarson	3.3	if (nn < NNEIGH) /* insert neighbor */
333			nn++;
334			for (n2 = nn; n2--; ) {
335	gwlarson	3.2	if (!n2 \|\| d >= nd[n2-1]) {
336			nd[n2] = d;
337			nl[n2][0] = h2;
338			nl[n2][1] = rnl[h2][n];
339			break;
340			}
341			nd[n2] = nd[n2-1];
342			nl[n2][0] = nl[n2-1][0];
343			nl[n2][1] = nl[n2-1][1];
344			}
345			}
346			}
347			return(nn);
348			}
349
350
351	gwlarson	3.7	meet_neighbors(nf, dp) /* run through samples and their neighbors */
352	gwlarson	3.2	int (*nf)();
353	gwlarson	3.7	char *dp;
354	gwlarson	3.2	{
355			short ln[NNEIGH][2];
356			int h, v, n, v2;
357			register short (*rnl)[NNEIGH];
358			/* initialize bottom row list */
359			rnl = (short ()[NNEIGH])malloc(NNEIGHsizeof(short)*hres);
360			CHECK(rnl==NULL, SYSTEM, "out of memory in meet_neighbors");
361			for (h = 0; h < hres; h++) {
362			for (n = v = 0; v < vres; v++)
363			if (CHK4(pixFlags, v*hres+h)) {
364			rnl[h][n++] = v;
365			if (n >= NNEIGH)
366			break;
367			}
368			while (n < NNEIGH)
369			rnl[h][n++] = NINF;
370			}
371			v = 0; /* do each row */
372			for ( ; ; ) {
373			for (h = 0; h < hres; h++) {
374			if (!CHK4(pixFlags, v*hres+h))
375			continue; /* no one home */
376			n = findneigh(ln, h, v, rnl);
377	gwlarson	3.7	(nf)(h, v, ln, n, dp); / call on neighbors */
378	gwlarson	3.2	}
379			if (++v >= vres) /* reinitialize row list */
380			break;
381			for (h = 0; h < hres; h++)
382			for (v2 = rnl[h][NNEIGH-1]+1; v2 < vres; v2++) {
383			if (v2 - v > v - rnl[h][0])
384			break; /* not close enough */
385			if (CHK4(pixFlags, v2*hres+h)) {
386			for (n = 0; n < NNEIGH-1; n++)
387			rnl[h][n] = rnl[h][n+1];
388			rnl[h][NNEIGH-1] = v2;
389			}
390			}
391			}
392			free((char )rnl); / free row list */
393	gwlarson	3.1	}