/* FIXME : think about clipping.  currently, we clip to the 4 sides of
the viewing cone, throw away z==0 as a special case, and then clip
again to the 4 sides of the 2D screen.  Can we do this more efficiently
by just clipping to z>+EPSILON in 3D before doing the 2D clipping? */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include "render.h"
#define EPSILON 1.0E-20

void printbasis(PG_BASIS *b)
{
	printf("(%g,%g,%g)\n",
		b->x.x,
		b->x.y,
		b->x.z);
	printf("(%g,%g,%g)\n",
		b->y.x,
		b->y.y,
		b->y.z);
	printf("(%g,%g,%g)\n",
		b->z.x,
		b->z.y,
		b->z.z);
}

/*-------------------------------------------*/
#if UCHAR_MAX!=0xFF
	#error This code assumes unsigned char is 8 bits
#endif

static union tagoutbuf {
	PIX_RGB32 p;
#if ULONG_MAX==0xFFFFFFFFL
	unsigned long l;
#elif UINT_MAX==0xFFFFFFFFL
	unsigned int l;
#elif USHRT_MAX==0xFFFFFFFFL
	unsigned short l;
#else
	#error This code needs a 32 bit integral type here.
#endif
} *outbuf;

#define MAXSPANS 1024
static struct minmax {
	int minx,maxx;
} spans[MAXSPANS];

static int osx,osy,centerx,centery;
static float *zbuf;
static float zoom;
static float detail;
static long thismark=0;

#define MAXDEPTH 10

static struct {
    PG_VECTOR center;
    PG_BASIS basis;
} coordstack[MAXDEPTH];

static int ncoord;
static PG_PLANE clips[5]={{0.0,1.0,1.0,0.0},
													{0.0,-1.0,1.0,0.0},
													{1.0,0.0,1.0,0.0},
													{-1.0,0.0,1.0,0.0},
													{0.0,0.0,1.0,EPSILON}};
static PG_FOG myfog;
#define MAXLIGHTS 5
static PG_LIGHT mylights[MAXLIGHTS];
static PG_LIGHT theirlights[MAXLIGHTS];
static int nlights=0;
/*-------------------------------------------*/

/*5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{ChangeCamera}

This function copies a bunch of data from the camera structure into
static variables for later use by the library.  It also calculates the
clipping planes for the current camera position.  Currently the
polygon renderer uses an array of maximum size $1024$ to hold scan
line data, so ChangeCamera will fail if the bitmap height is larger than
$1024$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */

short ChangeCamera(PG_CAMERA *c)
{
float l;
int i;

	if (c->bmap->height>MAXSPANS) {
		return 0;
	}
	osx=c->bmap->dx;
	osy=c->bmap->height;
	outbuf=(void *)(c->bmap->buffer+c->bmap->top*osx+c->bmap->left);
	centerx=c->bmap->width>>1;
	centery=osy>>1;
	zbuf=c->zbuf;
	detail=c->detail;

/* set up coordinates based on the camera */
	zoom=c->zoom;
	ncoord=0;
	transpose(coordstack[ncoord].basis,c->basis);
	vecmult(coordstack[ncoord].center,c->center,coordstack[ncoord].basis);
	coordstack[ncoord].center.x=-coordstack[ncoord].center.x;
	coordstack[ncoord].center.y=-coordstack[ncoord].center.y;
	coordstack[ncoord].center.z=-coordstack[ncoord].center.z;

/* transform lights into new coordinates */
	for (i=0;i<nlights;i++) {
		vecmult(mylights[i].pos,theirlights[i].pos,coordstack[ncoord].basis);
		if (mylights[i].type==2) {
			mylights[i].pos.x+=coordstack[ncoord].center.x;
			mylights[i].pos.y+=coordstack[ncoord].center.y;
			mylights[i].pos.z+=coordstack[ncoord].center.z;
		}
	}

/* set up clipping */
	clips[0].y=zoom;
	clips[0].z=centery+0.5;
	clips[1].y=-zoom;
	clips[1].z=centery+0.5;
	clips[2].x=zoom;
	clips[2].z=centerx+0.5;
	clips[3].x=-zoom;
	clips[3].z=centerx+0.5;

	l=mysqrt(clips[0].y*clips[0].y+clips[0].z*clips[0].z);
	clips[0].y=clips[0].y/l;
	clips[0].z=clips[0].z/l;
	clips[1].y=clips[1].y/l;
	clips[1].z=clips[1].z/l;
	l=mysqrt(clips[3].x*clips[3].x+clips[3].z*clips[3].z);
	clips[2].x=clips[2].x/l;
	clips[2].z=clips[2].z/l;
	clips[3].x=clips[3].x/l;
	clips[3].z=clips[3].z/l;

	return 1;
}

/*5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{ChangeEnvironment}
This function copies lighting and fog data to static variables.  It
also translates the lighting vectors into local camera coordinates.
The lights are kept in both transformed versions for use by the
rendering routines, and untransformed versions which are needed by
the ChangeCamera function when it switches to a new coordinate system.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */

short ChangeEnvironment(PG_ENVIRONMENT *env)
{
int i;
	if (env) {
		nlights=env->nlights;
		if (nlights>MAXLIGHTS)
			nlights=MAXLIGHTS;
		memcpy(theirlights,env->lights,nlights*sizeof(PG_LIGHT));
		for (i=0;i<nlights;i++) {
			mylights[i].type=theirlights[i].type;
			mylights[i].red=theirlights[i].red;
			mylights[i].green=theirlights[i].green;
			mylights[i].blue=theirlights[i].blue;
			vecmult(mylights[i].pos,theirlights[i].pos,coordstack[ncoord].basis);
			if (mylights[i].type==2) {
				mylights[i].pos.x+=coordstack[ncoord].center.x;
				mylights[i].pos.y+=coordstack[ncoord].center.y;
				mylights[i].pos.z+=coordstack[ncoord].center.z;
			}
		}
		if (env->fog)
			memcpy(&myfog,env->fog,sizeof(PG_FOG));
		else
			myfog.type=0;
	}
	else {
		nlights=0;
		myfog.type=0;
	}
	return 1;
}

/*5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{BeginRender3D and EndRender3D}

BeginRender3D sets up the camera and environment bu calling
ChangeCamera and ChangeEnvironment, and sets all objects to the
unrendered state.  This is done even though the library doesn't have a
complete list of objects by changing the current value of thismark,
which is compared against the objects' marks to see if they have been
drawn already.  There is a slight problem with this method, because
an object might have a really old mark, so that a false match occurs
when thismark has rolled all the way around and happens to match the
old value.  This will occur only if an object remains undrawn for
exactly $2^{32}$ frames.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */

short BeginRender3D(PG_CAMERA *c,PG_ENVIRONMENT *env)
{

	nlights=0;
	if (!ChangeCamera(c))
		return 0;

	if (!ChangeEnvironment(env))
		return 0;

	/* adjust frame counter to tell window renderer that we're on a new frame */
	thismark++; 
	return 1;
}

void EndRender3D(PG_CAMERA *c)
{
	//This doesn't do anything right now.
}

/*5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{IsRendered, MarkRendered, and MarkUnrendered}
These manipulate the object's mark in various ways.  The library
assumes that an object has been rendered if and only if thismark
matches the object's mark.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */

int IsRendered(PG_OBJECT *obj)
{
	if (obj->mark==thismark)
		return 1;
	return 0;
}

void MarkRendered(PG_OBJECT *obj)
{
	obj->mark=thismark;
}

void MarkUnrendered(PG_OBJECT *obj)
{
	obj->mark=thismark-1;
}

/*5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{DrawBack}

This function renders a type $0$ background.  It doesn't actually compute
the color of each pixel correctly.  Instead it computes the color of the
center pixel, and the derivatives of the color components:

\begin{eqnarray*}
center\ color 
& = & 
(axis \cdot basis_z)\cdot (pole\ color)+(equator\ color)\\
\frac{\partial}{\partial x}color  
& = & 
\frac{(axis \cdot basis_x)\cdot (pole\ color)}{zoom}\\
\frac{\partial}{\partial y}color  
& = & 
\frac{(axis \cdot basis_y)\cdot (pole\ color)}{zoom}\\
\end{eqnarray*} 

These derivatives are used to linearly approximate the color of each pixel.
The approximation is better for larger values of \textc{zoom}. If 
$\sqrt{\textc{centerx}^2 + \textc{centery}^2} > \textc{zoom}$ then the
approximation can cause colors to be computed which are outside the range
$\left[\textc{equ} - \textc{pole} , \textc{equ} - \textc{pole}\right]$.

This function is a good candidate for optimization, because for simple scenes
it can represent the majority of the pixels drawn.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */

#define CENTER 0x10000L
//FIXME DrawIt is a stupid name.
/* DrawIt is an experimental function which uses random noise as a 
texture for the background.  It does the same thing as DrawTable,
except that it takes two tables, and chooses randomly for each pixel
which one it selects the color from.  This effectively gives each color
band a two value dither instead of a pure color. */
void DrawIt(PG_VECTOR *axis,PG_BASIS *basis,float zoom,
							PIX_RGB32 *col1,PIX_RGB32 *col2)
{
int x,y;
long v,dvx,dvy,t;
float dot,f;
long p;
long r1; //For random textures

	dot= axis->x*basis->z.x
		+axis->y*basis->z.y
		+axis->z*basis->z.z;
		
	v=CENTER+(long)(dot*CENTER);
	
	dot= axis->x*basis->x.x
		+axis->y*basis->x.y
		+axis->z*basis->x.z;
	dot/=zoom;
	dvx=dot*CENTER;

	dot= axis->x*basis->y.x
		+axis->y*basis->y.y
		+axis->z*basis->y.z;
	dot/=zoom;
	dvy=dot*CENTER;

/* because this is a linear approximation, the max and min of v
	may go slightly out of bounds.  Here we check for this and compensate */
	t=labs(dvx)*centerx+labs(dvy)*centery;
	if (t+labs(v-CENTER)>CENTER) {
		f=(CENTER-labs(v-CENTER))/(float)t;
		dvx=f*dvx; /*this relies on rounding towards 0 */
		dvy=f*dvy; /*this relies on rounding towards 0 */
	}
	
/* adjust the value of v so it represents the upper left corner of the image */
	v-=dvx*centerx;
	v+=dvy*centery;
	
	dvy+=dvx*(centerx<<1);
	
	p=0;
	r1=v;
	for (y=0;y<(centery<<1);y++) {
		for (x=0;x<(centerx<<1);x++) {
			r1=(r1*0x3877)+0x12345673; //for random textures

			/*
			outbuf[p+x].p.red=col1[v>>8].red+(col2[v>>8].red&(r1>>24));
			outbuf[p+x].p.green=col1[v>>8].green+(col2[v>>8].green&(r1>>24));
			outbuf[p+x].p.blue=col1[v>>8].blue+(col2[v>>8].blue&(r1>>24));
			*/
			if (r1&0x40000000L)
				outbuf[p+x].p=col1[v>>8];
			else
				outbuf[p+x].p=col2[v>>8];
			v+=dvx;
		}
		p+=osx;
		v-=dvy;
	}
}
void DrawTable(PG_VECTOR *axis,PG_BASIS *basis,float zoom,
							PIX_RGB32 *cols)
{
int x,y;
long v,dvx,dvy,t;
float dot,f;
long p;

	dot= axis->x*basis->z.x
		+axis->y*basis->z.y
		+axis->z*basis->z.z;
		
	v=CENTER+(long)(dot*CENTER);
	
	dot= axis->x*basis->x.x
		+axis->y*basis->x.y
		+axis->z*basis->x.z;
	dot/=zoom;
	dvx=dot*CENTER;

	dot= axis->x*basis->y.x
		+axis->y*basis->y.y
		+axis->z*basis->y.z;
	dot/=zoom;
	dvy=dot*CENTER;

/* because this is a linear approximation, the max and min of v
	may go slightly out of bounds.  Here we check for this and compensate */
	t=labs(dvx)*centerx+labs(dvy)*centery;
	if (t+labs(v-CENTER)>CENTER) {
		f=(CENTER-labs(v-CENTER))/(float)t;
		dvx=f*dvx; /*this relies on rounding towards 0 */
		dvy=f*dvy; /*this relies on rounding towards 0 */
	}
	
/* adjust the value of v so it represents the upper left corner of the image */
	v-=dvx*centerx;
	v+=dvy*centery;
	
	dvy+=dvx*(centerx<<1);
	
	p=0;
	for (y=0;y<(centery<<1);y++) {
		for (x=0;x<(centerx<<1);x++) {
			if (zbuf[p+x]==0) 
				outbuf[p+x].p=cols[v>>8];
			v+=dvx;
		}
		p+=osx;
		v-=dvy;
	}
}
void DrawBack(PG_VECTOR *axis,PG_BASIS *basis,float zoom,
							PIX_RGB32 *equ,PIX_RGB32 *pole)
{
static PIX_RGB32 mytable[513];
int i;
	for (i=0;i<513;i++) {
		mytable[i].red=equ->red+(((i-256)*(signed char)pole->red)>>8);
		mytable[i].green=equ->green+(((i-256)*(signed char)pole->green)>>8);
		mytable[i].blue=equ->blue+(((i-256)*(signed char)pole->blue)>>8);
	}
	DrawTable(axis,basis,zoom,mytable);
}

/*5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{dopoly}

The dopoly function scan converts polygons and draws them into the 
offscreen frame buffer.  It processes polygons in two steps.  First,
each edge of the polygon is converted into a set of points, one for each
scan line.  These points are accumulated in the \textc{spans} array.
Because the polygons are assumed to be convex, there are only two edges
which touch each scan line (except at verteces where two edges may share
a point) so only a left point and right point need to be collected for
each scan line.  The \textc{eline} function accumulates these points for
each edge.

\label{zbuf}
After each edge has been added to the spans array, the polygon is written
into the frame buffer line by line.  As each pixel of the line is written,
it is compared with the Z buffer to be sure that it is not occluded by
a previously drawn polygon.  The values in the Z buffer are actually 
$\frac{1}{z}$ values because $\frac{1}{z}$ changes linearly with changes in
screen coordinates, while $z$ itself does not.  Using $\frac{1}{z}$ allows
the polyon display routines to calculate the Z buffer value for each point
in the polygon by adding a constant.  There are two polygon writing
routines: \textc{scanpoly} which does opaque polygons, and
\textc{scantranspoly} which does partially transparent polygons.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */

static int miny,maxy;
static void eline(int,int,int,int);
static void initpoly(void);
static void scanpoly(PIX_RGB32,float,float,float);
static void scantranspoly(PIX_RGB32,float,float,float);
static short scanwin(float,float,float);

static void dopoly
				(PIX_RGB32 c,int n,PG_POINT2D *p,float z0,float dzx,float dzy)
{
int i;
	initpoly();
	eline(p[0].x+centerx,centery-p[0].y,p[n-1].x+centerx,centery-p[n-1].y);
	for (i=1;i<n;i++) {
		eline(p[i-1].x+centerx,centery-p[i-1].y,p[i].x+centerx,centery-p[i].y);
	}
	dzy=-dzy;
	z0-=dzx*centerx;
	z0-=dzy*centery;

/* FIXME kluge to deal with precision problem... */
/* what's going on here is that the integer calculation in the rasterizer
is sometimes almost a pixel off the correct position for an edge.  This makes
the Z buffer values almost a pixel's worth different from ideal in the worst
case.  If you've got a plane that's nearly edge-on, one pixel of change is A
LOT of delta Z, so the plane can appear to stick out in front of things which
are quite a ways in front of it.  We deal with this by subtracting one pixel's
worth of Z in both the X and Y directions.  Now we're never in front of where
we're supposed to be, but sometimes behind.  This looks better, because an
almost edge-on plane can disappear behind things when it's not supposed to
without looking bad, but if it appears in front of where it should be this
is BAD. */

/* FIXME  : this kluge causes another problem.  Any polygons which are only
a pixel or so away from the horizon will get clipped because we move their
farthest edge's z coordinate to less than 0.  This maks some strange things
go on with polygons that the eye point is close to, or that are really
far away. */

	if (dzy<0)
		z0+=dzy;
	else
		z0-=dzy;
	if (dzx<0)
		z0+=dzx;
	else
		z0-=dzx;
		
	if ((c.alpha&0x0F)==0)
		scanpoly(c,z0,dzx,dzy);
	else 
		scantranspoly(c,z0,dzx,dzy);
}

/*5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{dowin}

This function scans a window's polygon the same way that dopoly does for
a real polygon, and returns TRUE if the polygon is visible.  Most of the
work is actually done in the scanwin function.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
static short dowin(int n,PG_POINT2D *p,float z0,float dzx,float dzy)
{
int i;

	initpoly();
	eline(p[0].x+centerx,centery-p[0].y,p[n-1].x+centerx,centery-p[n-1].y);
	for (i=1;i<n;i++) {
		eline(p[i-1].x+centerx,centery-p[i-1].y,p[i].x+centerx,centery-p[i].y);
	}
	dzy=-dzy;
	z0-=dzx*centerx;
	z0-=dzy*centery;

/* FIXME kluge to deal with precision problem... */

	if (dzy<0)
		z0+=dzy;
	else
		z0-=dzy;
	if (dzx<0)
		z0+=dzx;
	else
		z0-=dzx;
		
	return scanwin(z0,dzx,dzy);
}

/*5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{eline}

This routine is used by dopoly to compute the points where the specified
edge intersects each scan line of the image buffer.  Edges are always drawn
from top to bottom.  The main loop logic is in the code in two slightly
different forms, one for lines with positive slope, and one for negative
slope.  The macro \textc{spancmp} adds each point to the \textc{spans}
array.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */

#define spancmp(x,y) if (spans[y].minx>(x)) spans[y].minx=(x);\
										 if (spans[y].maxx<(x)) spans[y].maxx=(x);

static void eline(int x0,int y0,int x1,int y1)
{
int dx,dy,delta,x,y,inc;

	if (y0>y1) {
		y=y0;
		y0=y1;
		y1=y;
		x=x0;
		x0=x1;
		x1=x;
	}

	dx=x1-x0;
	
	if (y0<0) {
		x0-=dx*y0/(y1-y0);
		y0=0;
	}
	//FIXME : this should be osy-1 if we render the bottom lines of any
	// polygons
	if (y1>osy) {
		y1=osy;
	}
	dy=y1-y0;

	if (dy==0)
		return;
	if (y0<miny)
		miny=y0;
	if (y1>maxy)
		maxy=y1;
		
	if (dx<0) {
		inc=(-dx)/dy;
		dx+=inc*dy;
		delta=dy>>1;
		
		x=x0;
		for (y=y0;y<y1;y++) {
			spancmp(x,y);
			x-=inc;
			delta+=dx;
			if (delta<0) {
				delta+=dy;
				x--;
			}
		}
	}
	else {
		inc=dx/dy;
		dx-=inc*dy;
		delta=-(dy>>1);
	
		x=x0;
		for (y=y0;y<y1;y++) {
			spancmp(x,y);
			x+=inc;
			delta+=dx;
			if (delta>0) {
				delta-=dy;
				x++;
			}
		}		
	}
	spancmp(x,y);
}

/*5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{initpoly}

This routine clears out the spans array, and the miny and maxy pointers
so that a new polygon can be started.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */

static void initpoly(void)
{
int i;
	for (i=0;i<osy;i++) {
		spans[i].minx=osx;
		spans[i].maxx=0;
	}
	miny=osy;
	maxy=0;
}

/*5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{scanpoly}

This function writes the polygon to the offscreen buffer, using the scan
line coordinates calculated by repeated calls to eline.  It also does 
comparisons with the Z buffer for hidden surface removal.  The \textc{z0},
\textc{dx}, and \textc{dy} parameters are the value of the polygon's
Z buffer key at the upper left corner of the screen, and the change in the
key for each one pixel move in the x and y directions.  These numbers
allow the inner loop to keep an accumulator which represents the Z key of
the current pixel using only an addition for each new pixel.  This 
implementation does both floating point and integer math in the inner loop,
so it will work faster on machines which have parallel floating point and
integer execution.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
/* there are some pieces of code commented out of scanpoly which
implement random textures.  This isn't a particularly good implementation,
though, because the amount of randomness is set to a constant value.
If you are doing this for real, make the bit mask which r1 is anded with
a parameter somehow.  Maybe encode it in the top bits of the color's
alpha channel? */
static void scanpoly(PIX_RGB32 c,float z0,float dx,float dy)
{
int y;
int x;
int tx;
long lp;
float lz,tz;

//long r1=0; //For random textures
	lp=miny*osx;
	lz=z0+dy*miny;

	for (y=miny;y<maxy;y++) {
		if (spans[y].minx<0)
			spans[y].minx=0;
		if (spans[y].maxx>(centerx<<1))
			spans[y].maxx=centerx<<1;

		tz=spans[y].minx*dx+lz;
		tx=spans[y].maxx;
		for (x=spans[y].minx;x<tx;x++) {
			if (tz>zbuf[lp+x]) {
				outbuf[lp+x].p=c;
				//outbuf[lp+x].l+=(r1&0x1F0F1F0F); //for random textures
				zbuf[lp+x]=tz;
			}
			tz+=dx;
			//r1=(r1*0x38697183)+0x12345673; //for random textures
		}
		lp+=osx;
		lz+=dy;
	}
}

/*5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{scantranspoly}

The logic of this function is identical to that of scanpoly, except that
any time a pixel is written it is combined with the previous value of the
frame buffer, using the transparency value to choose proportions.  The
combination is done using bit shifting and multiplication, and only keeps
the most significant 4 bits of both the new color and the previous color.
This sometimes produces noticable color quantization, particularly when an
almost totally transparent object is placed in front of the smoothly shaded
background, but it allows all three color components to be processed in a
single step.

A significant deficiency of the transparent rendering routines as a whole is
that transparent polygons only allow more distant polygons to show through if
the distant polygons were drawn first.  This means that it is best to put all
transparent polygons near the end of the polygon list for each object.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */

static void scantranspoly(PIX_RGB32 col,float z0,float dx,float dy)
{
int y;
int x;
int tx;
int trans;
long lp;
float lz,tz;
union tagoutbuf c;
	c.p=col;
	trans=c.p.alpha&0x0F;
	c.l=((c.l&0xF0F0F0F0)>>4)*(0x10-trans);
	lp=miny*osx;
	lz=z0+dy*miny;
	for (y=miny;y<maxy;y++) {

		if (spans[y].minx<0)
			spans[y].minx=0;
		if (spans[y].maxx>centerx<<1)
			spans[y].maxx=centerx<<1;

		tz=spans[y].minx*dx+lz;
		tx=spans[y].maxx;
		for (x=spans[y].minx;x<tx;x++) {
			if (tz>zbuf[lp+x]) {
				outbuf[lp+x].l=c.l+((outbuf[lp+x].l&0xF0F0F0F0)>>4)*trans;
				zbuf[lp+x]=tz;
			}
			tz+=dx;
		}
		lp+=osx;
		lz+=dy;
	}
}

/*5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{scanwin}

This traverses pixels in the same way as scanpoly, but doesn't write any
points.  It just compares the z~buffer, and returns TRUE if any points
are visible.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
   
static short scanwin(float z0,float dx,float dy)
{
int y;
int x;
int tx;
long lp;
float lz,tz;

	lp=miny*osx;
	lz=z0+dy*miny;
	for (y=miny;y<maxy;y++) {

		if (spans[y].minx<0)
			spans[y].minx=0;
		if (spans[y].maxx>centerx<<1)
			spans[y].maxx=centerx<<1;

		tz=spans[y].minx*dx+lz;
		tx=spans[y].maxx;
		for (x=spans[y].minx;x<tx;x++) {
			if (tz>zbuf[lp+x]) {
				return 1;
			}
			tz+=dx;
		}
		lp+=osx;
		lz+=dy;
	}
	return 0;
}


/*5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{DrawPoly}
	This function draws a polygon by doing the perspective geometry
calculations to project down to 2~Dimensions, then using \textc{dopoly}
to render it.

The dopoly function calculates Z~buffer values using data supplied
by DrawPoly.  The Z~buffer is filled with $\frac{zoom}{z}$, since
this varies linearly with pixel coordinates.  It can be calcualted
as follows.

One way to express the equation for the plane in which the polygon
lies is
$$\vec{P} \cdot \vec{N} = c$$
where $\vec{P}$ are points in the plane, $\vec{N}$ is the plane's
normal vector, and $c$ is a constant.  Because this holds for any point
on the plane, we can find $c$ by choosing any point $\vec{P}_0$ and
substituting.

If $\vec{P}=(x,y,z)$, we can solve for $\frac{1}{z}$ as follows.

$$x\cdot N_x + y \cdot N_y + z \cdot N_z = c$$

$$\frac{x}{z} N_x + \frac{y}{z} N_y + N_z = \frac{c}{z} $$

$$ \frac{1}{z} = \frac{x}{z} \cdot \frac{N_x}{c} + 
	 \frac{y}{z} \cdot \frac{N_y}{c} + \frac{N_z}{c} $$

Now, substitute in the perspective projections of $x$ and $y$,
$x'=\frac{zoom\cdot x}{z}$ and $y'=\frac{zoom \cdot y}{z}$ to get

$$ \frac{zoom}{z} = x' \cdot \frac{N_x}{c} + y' \cdot \frac{N_y}{c} +
			\frac{zoom \cdot N_z}{c}$$

which is linear in $x'$ and $y'$.  DrawPoly passes the constant term and 
the coefficients of $x'$ and $y'$ to dopoly.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
#define MAXPTS 100
static PG_POINT2D points2d[MAXPTS];

void DrawPoly(PIX_RGB32 c,
						int n,PG_VECTOR *pts,
						PG_VECTOR *norm,float zoom)
{
int i;
float pdist,zz;
		pdist=norm->x*pts[0].x+norm->y*pts[0].y+norm->z*pts[0].z;
		if (pdist>EPSILON)
			return;
/* FIXME: kluge.
	If we have a plane that's really close to edge on, we might get divide by
	zero in the slope calculations below.  We don't want to just refuse to draw
	these polygons, becuase this can lead to thin holes in some types of mesh
	surfaces, so we draw them with a fudged Z-Slope.  it's unclear whether this
	really works.
	*/

		if (n>MAXPTS)
			n=MAXPTS;
		for (i=0;i<n;i++) {
			/* this stops us from drawing any polygons which touch the eye point.
			If we weren't doing clipping earlier, this could also throw away
			polygons which we want to keep. */
			if (pts[i].z==0.0)
				return;
			zz=zoom/pts[i].z;
			points2d[i].x=pts[i].x*zz;
			points2d[i].y=pts[i].y*zz;
		}
		if (pdist<-EPSILON) {
			dopoly(c,n,points2d,zoom*norm->z/pdist,norm->x/pdist,norm->y/pdist);  
		}
		else {
			dopoly(c,n,points2d,zoom/pts[0].z,0.0,0.0);  
		}
}

/*5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{CheckWin}

This does the perspective projection for a windows's polygon the same way
that DrawPoly does a normal polygon.  It then uses dowin to determine
whether any pixels of the polygon are visible.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
short CheckWin(int n,PG_VECTOR *pts,
						PG_VECTOR *norm,float zoom)
{
int i;
float pdist,zz;
	pdist=norm->x*pts[0].x+norm->y*pts[0].y+norm->z*pts[0].z;
	if (pdist>-EPSILON)
		return 0;

	if (n>MAXPTS)
		n=MAXPTS;
	for (i=0;i<n;i++) {
		/* this stops us from drawing any polygons which touch the eye point.
		If we weren't doing clipping earlier, this could also throw away
		polygons which we want to keep. */
		if (pts[i].z==0.0)
			return 0;
		zz=zoom/pts[i].z;
		points2d[i].x=pts[i].x*zz;
		points2d[i].y=pts[i].y*zz;
	}
	return dowin(n,points2d,zoom*norm->z/pdist,norm->x/pdist,norm->y/pdist);  
}

static PG_VECTOR polypts[MAXPTS];
static PG_VECTOR clippts[MAXPTS];

#define MAXOBJPTS 1100
static PG_VECTOR mypts[MAXOBJPTS];

#define MAXWINS 31

void DoColor(PIX_RGB32 *tcol,PG_VECTOR *tnorm,PG_VECTOR *pos);

/*5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{DrawObject}

This function does all of the coordinate transformations, lighting
calculations and recursion through subobjects and windows necessary
to render an object.  It uses a static stack of basis matrices and
center coordinates to do coordinate transformations for subobjects.

The coordinate stack is initialized to the inverse of the camera's
coordinates by BeginRender3D or ChangeCamera.  The inverse of the
camera's coordinates relative to the world coordinate system are
the coordinates of the world in the camera's coordinate system, so
anything transformed using this initial corrdinate stack will be
converted from the world's coordinates to the camera's coordinates.

Before each object is drawn, the current top of the coordinate stack is
transformed using the object's basis and center to produce a new top
element on the stack.  This new element is the correct transformation
for converting points of that object into camera coordinates.  As 
subobjects are drawn, the stack gets deeper, and as objects are finished
elements are popped off of the stack.  This way, the entire chain of
transformations from camera coordinates to subobject coordinates does not
need to be recalculated for each subobject; there is only one transformation
step per object or subobject.

After setting up the coordinate stack, each object is tested for visibility
and a part structure is selected based on camera detail, object detail and
distance.  Each point in the part is then transformed to camera coordinates.

The polygons of the object are then clipped using ClipPoly, the colors
are adjusted according to the lighting model using the DoColor function,
and the polygon is rendered.

Subobjects are all drawn before the main objectg is drawn, which can be
the wrong thing to do when the subobjects have transparent parts.
Windows are all recursed into after the object is drawn.  The visibility
of each window is checked before any of the windows' objects are drawn,
because the window visibility check makes use of the transformed point
array which will be overwritten by processing a new object.  The
coordinate stack is popped before recursing into windows, which causes
the objects that are connected through windows to be drawn as if they are
siblings of the current object in the coordinate tree; that is, they act
as if the current object and the object seen through the window are
subobjects of the same parent.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
void DrawObject(PG_OBJECT *obj)
{
int i;
long winflags;
PG_OBJECTPART *part;
static int j,n;
static PG_POLYGON *poly;
static PG_POLYGON *win;
static float maxdetail;
static PG_VECTOR tnorm;
static PIX_RGB32 tcol;

	/* setting this mark prevents the window recursion code from redrawing
	this object. */
	obj->mark=thismark;

/* transform to camera coords: */
	if (ncoord>=MAXDEPTH-1) {
		puts("FAILED: too many levels of subobjects");
		return;  /* out of stack space */
	}
	n=ncoord;
	ncoord++;

	/* set up the new center coordinates */
	vecmult(coordstack[ncoord].center,obj->center,coordstack[n].basis);

	coordstack[ncoord].center.x+=coordstack[n].center.x;
	coordstack[ncoord].center.y+=coordstack[n].center.y;
	coordstack[ncoord].center.z+=coordstack[n].center.z;

	if (coordstack[ncoord].center.z>EPSILON) {
		maxdetail=zoom*detail/coordstack[ncoord].center.z;
	}
	else if (coordstack[ncoord].center.z<-EPSILON) {
		maxdetail=-zoom*detail/coordstack[ncoord].center.z;
	}
	else {
		maxdetail=zoom*detail*1000.0;
	}

	if (obj->detail>maxdetail) {
		ncoord--;
		return;
	}

	/* set up the new basis */
	matmult(coordstack[ncoord].basis,obj->basis,coordstack[n].basis);

	/* do a coarse visibility check */
	if ((dot(clips[0],coordstack[ncoord].center)<-obj->size)||
			(dot(clips[1],coordstack[ncoord].center)<-obj->size)||
			(dot(clips[2],coordstack[ncoord].center)<-obj->size)||
			(dot(clips[3],coordstack[ncoord].center)<-obj->size)) {
		ncoord--;
		return;
	}

/* do each subobject */
	for (i=0;i<obj->nobjects;i++) {
		DrawObject(obj->objects[i]);
	}

	if (obj->nparts==0) {
		ncoord--;
		return;
	}

	for (i=0;i<obj->nparts;i++) {
		part=&(obj->parts[i]);
		if (part->detail>maxdetail) {
			if (i>0)
				part=&(obj->parts[i-1]);
			break;
		}
	}

/* transform all points to new basis */
	if (part->npoints>MAXOBJPTS) {
		puts("FAILED: too many points per object");
		ncoord--;
		return;
	}
	for (i=0;i<part->npoints;i++) {
		vecmult(mypts[i],part->points[i],coordstack[ncoord].basis);
		mypts[i].x+=coordstack[ncoord].center.x;
		mypts[i].y+=coordstack[ncoord].center.y;
		mypts[i].z+=coordstack[ncoord].center.z;
	}

	for (i=0;i<part->npolys;i++) {
		poly=&(part->polys[i]);
	/* transform normal to new basis, check normal vector */
		vecmult(tnorm,part->norms[i],coordstack[ncoord].basis);

		if (poly->n>MAXPTS) {
			puts("FAILED: too many points per polygon");
			ncoord--;
			return;
		}
	/* copy polygon points from object point array */
		for (j=0;j<poly->n;j++) {
			polypts[j]=mypts[poly->points[j]];
		}
	/* clip to viewable range */
		n=ClipPoly(poly->n,polypts,clippts,4,clips);
		if (n<2)
			continue;

		tcol=part->colors[i];
		DoColor(&tcol,&tnorm,polypts);
		DrawPoly(tcol,n,clippts,&tnorm,zoom);
	}

/* do windows */
	if (obj->nwindows>MAXWINS) {
			puts("FAILED: too many windows per object");
			ncoord--;
			return;
	}
	/* find out for all windows whether they're visible or not first,
	then call DisplayObject on each of them.  This is because the vertices
	of the window are in a temp buffer from the transform step, and
	doing a new object will trash them. */
	
	winflags=0;
	for (i=0;i<obj->nwindows;i++) {
		if (obj->winobjects[i]==NULL) {
			continue;
		}
		
		win=&(part->windows[i]);
	/* transform normal to new basis */
		vecmult(tnorm,part->wnorms[i],coordstack[ncoord].basis);

		if (win->n>MAXPTS) {
			puts("FAILED: too many points per window");
			ncoord--;
			return;
		}
	/* copy polygon points from object point array */
		for (j=0;j<win->n;j++) {
			polypts[j]=mypts[win->points[j]];
		}
	/* clip to viewable range */
		n=ClipPoly(win->n,polypts,clippts,4,clips);
		if ((n>=2) && CheckWin(n,clippts,&tnorm,zoom)) {
			winflags|=1<<i;
		}
	}
	ncoord--;
	for (i=0;i<obj->nwindows;i++) {
		if (winflags&(1<<i)) {
			if (obj->winobjects[i]->mark!=thismark) {
				DrawObject(obj->winobjects[i]);
			}
		}
	}
	/* FIXME : there's a subtle bug here -- thismark rolls over once in a
			while, and this may cause windows to be ignored incorrectly, but
			on a 32 bit system, this will only happen if a window is out of
			view for exactly 2^32 frames.  I'm willing to accept a one frame
			graphics glitch every 4 billion frames rendered. */
}

/*5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{DoColor}

This function computes the change to the color of a polygon caused
by the current environment's lights and fog given it's normal vector
and position.  The position is given in camera coordinates.

The effect of ambient light is to multiply the three color components
of the polygon by the three color intensity values of the light.  For 
directed lights, the multiplier is the color intensity of the light
multiplied by the dot product of the normal vector and the light's
direction vector.  This produces bright illumination when the polygon
is facing the light, fading to $0$ as the polygon rotates to meet the
light edge-on.  If the dot product is less than $0$, the light makes
no contribution to the polygon's color.

Fog effects are computed using only the z coordinate of the polygon's
position, which is approximately the distance from the camera to the
polygon.  Fog effects are only computed at one point for each polygon,
which means that large polygons will be a single color rather than having
a smooth shading from normal to fogged as one might like.  A shaded effect
can be achieved by tiling the surface with many small polygons rather than
using one large one.  The particular case of an infinite flat plane with
fog can be simulated with proper use of the DrawTable function.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */

void DoColor(PIX_RGB32 *tcol,PG_VECTOR *tnorm,PG_VECTOR *pos)
{
int i;
float red,green,blue,tf;
PG_LIGHT *light;

	/* calculate illumination unless there are no lights, or
	unless the object is self-luminous */
	if ((nlights)&&((tcol->alpha&0x10)==0)) {
		red=0.0;
		green=0.0;
		blue=0.0;
		for (i=0;i<nlights;i++) {
			light=&(mylights[i]);
			switch (light->type) {
				case 0:
					red+=light->red;
					green+=light->green;
					blue+=light->blue;
					break;
				case 1:
					tf=light->pos.x*tnorm->x+light->pos.y*tnorm->y+light->pos.z*tnorm->z;
					if (tf<0)
					break;
					red+=tf*light->red;
					green+=tf*light->green;
					blue+=tf*light->blue;
					break;
// FIXME support case 3: local light source
			}
		}
	}
	else {
		if (myfog.type==0)
			return;
		red=1.0;
		green=1.0;
		blue=1.0;
	}

	red*=tcol->red;
	blue*=tcol->blue;
	green*=tcol->green;

/* do fog */
	switch (myfog.type) {
		case 1:
			if (pos->z>myfog.stop_dist) {
				*tcol=myfog.col;
				return;
			}
			if (pos->z>myfog.start_dist) {
				tf=1.0/(myfog.stop_dist-myfog.start_dist);
				red=(myfog.col.red*(pos->z-myfog.start_dist)+
						red*(myfog.stop_dist-pos->z))*tf;
				green=(myfog.col.green*(pos->z-myfog.start_dist)+
						green*(myfog.stop_dist-pos->z))*tf;
				blue=(myfog.col.blue*(pos->z-myfog.start_dist)+
						blue*(myfog.stop_dist-pos->z))*tf;
			}
	}
/* write values back */
	tcol->red=255;
	tcol->green=255;
	tcol->blue=255;
	
	i=(int)red;
	if (i<256)
		tcol->red=i;
	i=(int)blue;
	if (i<256)
		tcol->blue=i;
	i=(int)green;
	if (i<256)
		tcol->green=i;
}