graph3d.c

Microwindows genesis was with the NanoGUI project, and is now the primary distribution for both th

/*
 * Copyright (c) 1999 Greg Haerr <[e-mail]>
 *
 * 3D Graphics Library for Micro-Windows
 */
#define MWINCLUDECOLORS
#include "windows.h"
#include "device.h"
#include "graph3d.h"
#define USEBLIT		1	/* =1 to use memDC's*/

static int	nxpix;
static int	nypix;
static vec1	xscale;
static vec1	yscale;
static vec3	eye;
static vec3	direct;
static double 	Q[5][5];
static HDC	hdc;
static HDC	hdcMem;
static HBITMAP	hbmp, hbmpOrg;

/* setup eye, direction, calc observation matrix Q*/
void
look3(vec1 x, vec1 y, vec1 z)
{
	eye.x = x;
	eye.y = y;
	eye.z = z;
	direct.x = -eye.x;
	direct.y = -eye.y;
	direct.z = -eye.z;
	findQ();
}

void
init3(HDC hDC, HWND memhwnd)
{
	HBRUSH	hbr;

	hdc = hDC;
	if(hdc) {
		nxpix = hdc->hwnd->clirect.right - hdc->hwnd->clirect.left;
		nypix = hdc->hwnd->clirect.bottom - hdc->hwnd->clirect.top;
		xscale = (vec1)(nxpix-1) / nxpix * nxpix/2;
		yscale = (vec1)(nypix-1) / nypix * nypix/2;

		if(memhwnd) {
			hdcMem = CreateCompatibleDC(NULL);
			if(hdcMem) {
				hbmp = CreateCompatibleBitmap(hdcMem,
					nxpix, nypix);
				hbmpOrg = SelectObject(hdcMem, hbmp);
				hdc = hdcMem;
			}
			hbr = (HBRUSH)GetClassLong(memhwnd, GCL_HBRBACKGROUND);
			FillRect(hdc, NULL, hbr);
		}
		/* create pen for setcolor3() color override*/
		SelectObject(hdc, CreatePen(PS_SOLID, 1, BLACK));
	}
}

void
paint3(HDC hDC)
{
	if(hdcMem) {
		BitBlt(hDC, 0, 0, nxpix, nypix, hdcMem, 0, 0, SRCCOPY);
		DeleteObject(SelectObject(hdcMem, hbmpOrg));
		DeleteDC(hdcMem);
	}
	hdcMem = NULL;
}

int
fx(vec1 x)
{
	return (int)(x * xscale + nxpix*0.5 - 0.5);
}

int
fy(vec1 y)
{
	return (int)(y * yscale + nypix*0.5 - 0.5);
}

void
moveto3(vec2 pt)
{
	MoveToEx(hdc, fx(pt.x), fy(pt.y), NULL);
}

void
setcolor3(MWCOLORVAL c)
{
	if(hdc)
		hdc->pen->color = c;
}

void
lineto3(vec2 pt)
{
	LineTo(hdc, fx(pt.x), fy(pt.y));
}

void
polyfill(int n, vec2 points[])
{
	int	i;
	int	xoff, yoff;
	MWPOINT	pv[MAXPOLY];

	if(!hdc)
		return;

	/* calc window offset*/
	xoff = hdc->hwnd->clirect.left;
	yoff = hdc->hwnd->clirect.top;

	/* only plot non-trivial polygons*/
	if(n > 2) {
		for(i=0; i<n; ++i) {
			pv[i].x = fx(points[i].x) + xoff;
			pv[i].y = fy(points[i].y) + yoff;
			/* fix: floating round error, y intercept difference
			 * with GdLine
			 */
			/*pv[i].x = fx(points[i].x + xoff);*/
			/*pv[i].y = fy(points[i].y + yoff);*/
		}
		GdSetForeground(GdFindColor(hdc->pen->color));
		GdFillPoly(hdc->psd, n, pv);
	}
}

void
square(void)
{
	vec2	pt0, pt1, pt2, pt3;

	pt0.x = -1; pt0.y = 1;
	pt1.x = -1; pt1.y = -1;
	pt2.x = 1; pt2.y = -1;
	pt3.x = 1; pt3.y = 1;
	moveto3(pt0);
	lineto3(pt1);
	lineto3(pt2);
	lineto3(pt3);
	lineto3(pt0);
}

void
circle3(vec1 r)
{
	vec1 	theta = 0;
	vec1	thinc = 2*pi/100;
	int 	i;
	vec2	pt;

	pt.x = r;
	pt.y = 0.0;
	moveto3(pt);

	for(i=0; i<100; ++i) {
		theta = theta + thinc;
		pt.x = r*cos(theta);
		pt.y = r*sin(theta);
		lineto3(pt);
	}
}

void
daisy(vec1 r,int points)
{
	int	i, j;
	vec1	theta = 0;
	vec1	thinc;
	vec2	pt[100];

	/* calculate n points on a circle*/
	thinc = 2*pi/points;
	for(i=0; i<points; ++i) {
		pt[i].x = r*cos(theta);
		pt[i].y = r*sin(theta);
		theta += thinc;
	}

	/* join point i to point j for all 0 <= i < j < n */
	for(i=0; i<points-1; ++i) {
		for(j=i+1; j<points; ++j) {
			moveto3(pt[i]);
			lineto3(pt[j]);
		}
	}
}

void
rose(vec1 r,int levels,int points)
{
	int	i, j, m, n;
	vec1	r1, theta, thinc;
	vec2	inner[100];
	vec2	outer[100];
	vec2	triangle[3];

	m = levels;
	n = points;
	thinc = 2*pi/n;

	/* initial inner circle*/
	for(i=0; i<n; ++i) {
		inner[i].x = 0.0;
		inner[i].y = 0.0;
	}

	/* loop thru m levels*/
	for(j=1; j<=m; ++j) {
		theta = -j*pi/n;
		r1 = r * (vec1)j/m;

		/* calc n points on outer circle*/
		for(i=0; i<n; ++i) {
			theta += thinc;
			outer[i].x = r1*cos(theta);
			outer[i].y = r1*sin(theta);
		}

		/* construct/draw triangles with vertices on
		 * inner and outer circles
		 */
		for(i=0; i<n; ++i) {
			triangle[0] = outer[i];
			triangle[1] = outer[(i+1) % n];
			triangle[2] = inner[i];

			/* fill triangle in red*/
			setcolor3(RED);
			polyfill(3, triangle);

#if 1
			/* outline triangle in white*/
			setcolor3(WHITE);
			moveto3(triangle[0]);
			lineto3(triangle[1]);
			lineto3(triangle[2]);
			lineto3(triangle[0]);
#endif
		}

		/* copy points on outer circle to inner arrays*/
		for(i=0; i<n; ++i)
			inner[i] = outer[i];
	}
}

/* draw a triangle with cordners v0, v1, v2*/
void
triangle(vec2 v0, vec2 v1, vec2 v2)
{
	vec2	poly[3];

	poly[0] = v0;
	poly[1] = v1;
	poly[2] = v2;

	setcolor3(GREEN);
	polyfill(3, poly);
	setcolor3(BLACK);
	moveto3(poly[2]);
	lineto3(poly[0]);
	lineto3(poly[1]);
	lineto3(poly[2]);
}

/* draw a quadrilateral with corners v0, v1, v2, v3*/
void
quadrilateral(vec2 v0, vec2 v1, vec2 v2, vec2 v3)
{
	vec2	poly[4];

	poly[0] = v0;
	poly[1] = v1;
	poly[2] = v2;
	poly[3] = v3;
	setcolor3(GREEN);
	polyfill(4, poly);
	setcolor3(BLACK);
	moveto3(poly[3]);
	lineto3(poly[0]);
	lineto3(poly[1]);
	lineto3(poly[2]);
	lineto3(poly[3]);
}

/* find intersection of lines v0 to v1 and v2 to v3*/
static int
patch(vec2 v0, vec2 v1, vec2 v2, vec2 v3)
{
	vec1	denom;
	vec1	mu;
	vec2	v4;

	denom = (v1.x-v0.x)*(v3.y-v2.y) - (v1.y-v0.y)*(v3.x-v2.x);
	if(fabs(denom) > epsilon) {
		mu = ((v2.x-v0.x)*(v3.y-v2.y) - (v2.y-v0.y)*(v3.x-v2.x))/denom;

		/* if intersection between lines v0 to v1 and v2 to v3,
		 * call it v4 and form triangles v0,v2,v4 and v1,v3,v4
		 */
		if(mu >= 0 && mu <= 1) {
			v4.x = (1-mu)*v0.x + mu*v1.x;
			v4.y = (1-mu)*v0.y + mu*v1.y;
			triangle(v0, v2, v4);
			triangle(v1, v3, v4);
			return 0;
		}
	}

	/* else find intersection of lines v0 to v2 and v1 to v3*/
	denom = (v2.x-v0.x)*(v3.y-v1.y) - (v2.y-v0.y)*(v3.x-v1.x);
	if(fabs(denom) > epsilon) {
		mu = ((v1.x-v0.x)*(v3.y-v1.y) - (v1.y-v0.y)*(v3.x-v1.x))/denom;

		/* if intersection between v0 and v1, call it v4
		 * and form triangles v0,v1,v4 and v2,v3,v4
		 */
		if(mu >= 0 && mu <= 1) {
			v4.x = (1-mu)*v0.x + mu*v2.x;
			v4.y = (1-mu)*v0.y + mu*v2.y;
			triangle(v0, v1, v4);
			triangle(v2, v3, v4);
			return 0;
		}
	}

	/* there are no proper intersections so form quadrilateral v0,v1,v3,v2*/
	quadrilateral(v0, v1, v3, v2);
	return 1;
}

/* plotted function*/
static vec1
plotfn(vec1 x, vec1 z)
{
	vec1	t;

	/* y = 4sin(sqrt(x*x+z*z))/sqrt(x*x+z*z) */
	t = sqrt(x*x + z*z);
	if(fabs(t) < epsilon)
		return 4.0;
	return 4.0 * sin(t) / t;
}

/* draw mathematical function plotfn*/
void
drawgrid(vec1 xmin, vec1 xmax, int nx, vec1 zmin, vec1 zmax, int nz)
{
	int	i, j;
	vec1	xi, xstep, yij;
	vec1	zj, zstep;
	vec2	v[2][100];
	double	S[5][5];

	/* scale it down*/
	scale3(1.0/(xmax-xmin)*2, 1.0/(xmax-xmin)*2, 1.0/(zmax-zmin), S);
	mult3(Q, S, Q);

	/* grid from xmin to xmax in nx steps and zmin to xmax in nz steps*/
	xstep = (xmax-xmin)/nx;
	zstep = (zmax-zmin)/nz;
	xi = xmin;
	zj = zmin;

	/* calc grid points on first fixed-z line, fine the y-height
	 * and transfrorm the points (xi,yij,zj) into observed
	 * position.  Observed first set stored in v[0,1..nx]
	 */
	for(i=0; i<=nx; ++i) {
		yij = plotfn(xi, zj);
		v[0][i].x = Q[1][1]*xi + Q[1][2]*yij + Q[1][3]*zj;
		v[0][i].y = Q[2][1]*xi + Q[2][2]*yij + Q[2][3]*zj;
		xi += xstep;
	}

	/* run thru consecutive fixed-z lines (the second set)*/
	for(j=0; j<nz; ++j) {
		xi = xmin;
		zj += zstep;

		/* calc grid points on this second set, find the
		 * y-height and transform the points (xi,yij,zj)
		 * into observed position.  Observed second set
		 * stored in v[1,0..nx]
		 */
		for(i=0; i<=nx; ++i) {
			yij = plotfn(xi, zj);
			v[1][i].x = Q[1][1]*xi + Q[1][2]*yij + Q[1][3]*zj;
			v[1][i].y = Q[2][1]*xi + Q[2][2]*yij + Q[2][3]*zj;
			xi += xstep;
		}

		/* run thru the nx patches formed by these two sets*/
		for(i=0; i<nx; ++i)
			patch(v[0][i], v[0][i+1], v[1][i], v[1][i+1]);

		/* copy second set into first set*/
		for(i=0; i<=nx; ++i)
			v[0][i] = v[1][i];
	}
}

/* returns the angle whose tangent is y/x.
 * all anomalies such as x=0 are also checked
 */
vec1
angle(vec1 x, vec1 y)
{
	if(fabs(x) < epsilon)
		if(fabs(y) < epsilon)
			return 0.0;
		else 
			if(y > 0.0)
				return pi*0.5;
			else return pi*1.5;
	else
		if(x < 0.0)
			return atan(y/x) + pi;
		else return atan(y/x);
}

/* calc 3d scaling matrix A giving scaling vector sx,sy,sz.
 * one unit on the x axis becomes sx units, one unit on y, sy,
 * and one unit on the z axis becomes sz units
 */
void
scale3(vec1 sx, vec1 sy, vec1 sz, double A[][5])
{
	int	i, j;

	for(i=1; i<5; ++i)
		for(j=1; j<5; ++j)
			A[i][j] = 0.0;
	A[1][1] = sx;
	A[2][2] = sy;
	A[3][3] = sz;
	A[4][4] = 1.0;
}

/* calc 3d axes translation matrix A
 * origin translated by vectdor tx,ty,tz
 */
void
tran3(vec1 tx, vec1 ty, vec1 tz, double A[][5])
{
	int	i, j;

	for(i=1; i<5; ++i) {
		for(j=1; j<5; ++j)
			A[i][j] = 0.0;
		A[i][i] = 1.0;
	}
	A[1][4] = -tx;
	A[2][4] = -ty;
	A[3][4] = -tz;
}

/* calc 3d axes rotation matrix A.  The axes are
 * rotated anti-clockwise through an angle theta radians
 * about an axis specified by m: m=1 means x, m=2 y, m=3 z axis
 */
void
rot3(int m, vec1 theta, double A[][5])
{
	int	i, j, m1, m2;
	vec1	c, s;

	for(i=1; i<5; ++i) 
		for(j=1; j<5; ++j)
			A[i][j] = 0.0;
	A[m][m] = 1.0;
	A[4][4] = 1.0;
	m1 = (m % 3) + 1;
	m2 = (m1 % 3) + 1;
	c = cos(theta);
	s = sin(theta);
	A[m1][m1] = c;
	A[m2][m2] = c;
	A[m1][m2] = s;
	A[m2][m1] = s;
}

/* calc the matrix product C of two matrices A and B*/
void
mult3(double A[][5], double B[][5], double C[][5])
{
	int	i, j, k;
	vec1	ab;

	for(i=1; i<5; ++i) 
		for(j=1; j<5; ++j) {
			ab = 0;
			for(k=1; k<5; ++k)
				ab += A[i][k] * B[k][j];
			C[i][j] = ab;
		}
}

/* calc observation matrix Q for given observer*/
void
findQ(void)
{
	vec1	alpha, beta, gamma, v, w;
	double	E[5][5];
	double	F[5][5];
	double	G[5][5];
	double	H[5][5];
	double	U[5][5];

	/* calc translation matrix F*/
	tran3(eye.x, eye.y, eye.z, F);

	/* calc rotation matrix G*/
	alpha = angle(-direct.x, -direct.y);
	rot3(3, alpha, G);

	/* calc rotation matrix H*/
	v = sqrt(direct.x*direct.x + direct.y*direct.y);
	beta = angle(-direct.z, v);
	rot3(2, beta, H);

	/* calc rotation matrix U*/
	w = sqrt(v*v + direct.z*direct.z);
	gamma = angle(-direct.x*w, direct.y*direct.z);
	rot3(3, -gamma, U);

	/* combine the transformations to find Q*/
	mult3(G, F, Q);
	mult3(H, Q, E);
	mult3(U, E, Q);
}