text: #ifndef lint
  static const char	RCSid[] = "$Id: rhpict2.c,v 3.18 2020/05/14 20:58:03 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 */
  }
  