src/hd/rhpict2.c

#ifndef lint
static const char       RCSid[] = "$Id: rhpict2.c,v 3.17 2018/10/05 19:19:16 greg Exp $";
#endif
/*
 * Rendering routines for rhpict.
 */

#include <string.h>

#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          5       /* 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)    memset((char *)(f),'\0',FL4NELS(n)*sizeof(int32))
#endif

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

#define isqrt(i2)       (isqrttab[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. */

extern void pixFinish(double ransamp);
extern void pixBeam(BEAM *bp, HDBEAMI *hb);

typedef int (sampfunc)(int h, int v, short nl[NNEIGH][2], int nd[NNEIGH],
                int n, double *rf);
static sampfunc kill_occl;
static sampfunc smooth_samp;
static sampfunc random_samp;
static void meet_neighbors(int occ, sampfunc *nf, double *dp);
static void reset_flags(void);
static void init_wfunc(void);
static int findneigh(short nl[NNEIGH][2], int nd[NNEIGH], int h, int v,
                short (*rnl)[NNEIGH]);


void
pixBeam(                        /* render a particular beam */
        BEAM    *bp,
        HDBEAMI *hb
)
{
        GCOORD  gc[2];
        RAYVAL  *rv;
        FVECT   rorg, rdir, wp, ip;
        double  d, prox;
        COLOR   col;
        int     n;
        int32   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.;
                }
                if (viewloc(ip, &myview, wp) != VL_GOOD)
                        continue;               /* frustum clipping */
                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];
        }
}


static int
kill_occl(              /* check for occlusion errors */
        int     h,
        int     v,
        short   nl[NNEIGH][2],
        int     nd[NNEIGH],
        int     n,
        double *rf
)
{
        short   forequad[2][2];
        int     d;
        int     i;
        int32   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 = isqrt(nd[i]);
                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);
}


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

        if (n <= 0)
                return(1);
        p = v*hres + h;                         /* build kernel values */
        maxr2 = nd[n-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--; ) {
                r = isqrt(nd[i]);
                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);
}


static int
random_samp(    /* gather samples randomly */
        int     h,
        int     v,
        short   nl[NNEIGH][2],
        int     nd[NNEIGH],
        int     n,
        double  *rf
)
{
        float   rnt[NNEIGH];
        double  rvar;
        int32   p, pn;
        int     ni;

        if (n <= 0)
                return(1);
        p = v*hres + h;
        if (*rf <= FTINY)               /* straight Voronoi regions */
                ni = 0;
        else {                          /* weighted choice */
                DCHECK(nd[n-1]>=MAXRAD2, CONSISTENCY, "out of range neighbor");
                rnt[0] = pixWeight[nd[0]];
                for (ni = 1; ni < n; ni++)
                        rnt[ni] = rnt[ni-1] + pixWeight[nd[ni]];
                rvar = rnt[n-1]*pow(frandom(), 1. / *rf);
                for (ni = 0; rvar > rnt[ni]+FTINY; ni++)
                        ;
        }
        pn = nl[ni][1]*hres + nl[ni][0];
        addcolor(mypixel[p], mypixel[pn]);
        myweight[p] += myweight[pn];
        return(1);
}


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


static void
reset_flags(void)                       /* allocate/set/reset occupancy flags */
{
        int32   p;

        if (pixFlags == NULL) {
                pixFlags = (int32 *)calloc(FL4NELS(hres*vres), sizeof(int32));
                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);
}


static void
init_wfunc(void)                        /* initialize weighting function */
{
        int     r2;
        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;
}


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

        nd[NNEIGH-1] = MAXRAD2;
        for (hoff = 0; 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 == 0) | (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);
}


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