/* test program for catmull-rom bicubic splines. */

#include "exec/types.h"
#include "intuition/intuition.h"
#include <macros.h>

struct Library *IntuitionBase=NULL, *GfxBase=NULL;
void *OpenLibrary();

struct Window *Window=NULL, *OpenWindow();
struct Screen *Screen=NULL, *OpenScreen();
struct RastPort *rp;							/* rast port */

struct TmpRas TmpRas;
struct AreaInfo AreaInfo;
PLANEPTR TRBuf;
UWORD AIBuf[20];

USHORT colors[32] = {
	0x000,0xFFF,0x000,0x111,0x222,0x333,0x444,0x555,
	0x666,0x777,0x888,0x999,0xAAA,0xBBB,0xCCC,0xDDD,
	0xEEE,0xFFF,0xB86,0x963,0x740,0x400,0xD80,0x666,
	0xDB8,0xABF,0xDB8,0xC96,0xB86,0xBA8,0xA00,0xC08
};

struct NewScreen NewScreen = { 0,0,320,200, 5, 0,1, NULL, CUSTOMSCREEN, };

struct NewWindow back_ground =
{	0,0,320,200, 0,1,
	MENUPICK | MOUSEMOVE,
	NOCAREREFRESH | ACTIVATE | SMART_REFRESH | BORDERLESS | REPORTMOUSE,
	NULL,NULL,NULL,NULL,NULL,320,200,320,200,CUSTOMSCREEN
};

#define DF3	ITEMTEXT | ITEMENABLED | HIGHCOMP | COMMSEQ

struct IntuiText quit_tx = { 5,1,JAM2,4,1,NULL,(UBYTE *)"Exit" };
struct MenuItem quit_mi = { NULL,0,20,190,10,DF3,0,(APTR)&quit_tx,0,'X' };
struct Menu menu1 = { NULL, 1,0,56,10, MENUENABLED | MIDRAWN, "File", &quit_mi };

#define SCAT	goto exit_pgm;

struct Point {
	WORD	x,y;
};

typedef struct {
	float		x,y,z;
} Point3;

typedef struct {							/* varying point					*/
	Point3		p,du,dv;
} Point3V;

Point3 control_hull[10][10];
float m[4][4];							/* spline basis weighting factors */

struct Point coord2[4];

extern UWORD misc_grid_data[4][8];

long color_x;
short color = 30;

transform(Point3 *p, Point *s)
{	s->x = 160 + (WORD)((p->x - p->y) * 60.0);
	s->y = 260 - (WORD)((p->x + p->y + p->z*2.5) * 25.0);
}

struct BitMap *bm_draw;

WORD	ustep = 1, vstep = 1;

calc_cardinal(float tension)
{
	m[0][0] = -tension;
	m[0][1] = 2.0 - tension;
	m[0][2] = tension - 2.0;
	m[0][3] = tension;

	m[1][0] = 2*tension;
	m[1][1] = tension - 3.0;
	m[1][2] = 3.0 - 2.0*tension;
	m[1][3] = -tension;

	m[2][0] = -tension;
	m[2][1] = 0.0;
	m[2][2] = tension;
	m[2][3] = 0.0;

	m[3][0] = 0.0;
	m[3][1] = 1.0;
	m[3][2] = 0.0;
	m[3][3] = 0.0;
}

main()
{	int i,j,k;

	if ((IntuitionBase = OpenLibrary("intuition.library",0)) == NULL) SCAT;
	if ((GfxBase = OpenLibrary("graphics.library",0)) == NULL) SCAT;
	if ((Screen = OpenScreen(&NewScreen)) == NULL) SCAT;
	LoadRGB4(&(Screen->ViewPort),colors,32);

	back_ground.Screen = Screen;
	if ((Window = OpenWindow(&back_ground)) == NULL) SCAT;
	rp = Window->RPort;

	TRBuf = (void *)AllocRaster(320,200);
	if (!TRBuf) SCAT;
	InitTmpRas(&TmpRas,TRBuf,RASSIZE(320,200));
	rp->TmpRas = &TmpRas;
	InitArea(&AreaInfo,AIBuf,5);
	rp->AreaInfo = &AreaInfo;

	SetMenuStrip(Window,&menu1);

	SetAPen(rp,1); SetRast(rp,0);
	bm_draw = rp->BitMap;

	WaitBlit();

	calc_cardinal(0.5);

		/* initialize control points */

	for (i=0; i<8; i++)
	{	for (j=0; j<8; j++)
		{	control_hull[i][j].x = (float)i;
			control_hull[i][j].y = (float)j;
			control_hull[i][j].z = 0.0;
		}
	}

	control_hull[2][2].z = 0.9;
	control_hull[2][3].z = 1.5;
	control_hull[3][3].z = 0.5;
	control_hull[3][4].z = -0.5;
	control_hull[4][2].z = 0.7;

	for (i=3; i>=0; i--)
	{	for (j=3; j>=0; j--) do_patch(i,j);
	}

	for (;;)
	{	struct IntuiMessage *message, *GetMsg();

		Wait(1 << Window->UserPort->mp_SigBit);

		while (message = GetMsg(Window->UserPort))
		{	ULONG class; USHORT code;
			class = message->Class;
			code = message->Code;
			ReplyMsg(message);

			if (class == MENUPICK && code != MENUNULL) SCAT;
		}
	}
exit_pgm:
	if (TRBuf) { FreeRaster(TRBuf,320,200); rp->TmpRas = NULL; rp->AreaInfo = NULL; }
	if (Window){ ClearMenuStrip(Window); CloseWindow(Window); }
	if (Screen) CloseScreen(Screen); 
	if (GfxBase) CloseLibrary(GfxBase);
	if (IntuitionBase) CloseLibrary(IntuitionBase);
}

typedef Point3		SplineBasis[4][4];

interpolate_points(i,j,u,v,p) WORD i,j; float u,v; Point3 *p;
{	float uv0, uv1, uv2, uv3;

	uv0 = u * v;
	uv1 = (1.0 - u) * v;
	uv2 = (1.0 - u) * (1.0 - v);
	uv3 = u * (1.0 - v);

		/* first point */

	p->x =  control_hull[i  ][j  ].x * uv0 + 
			control_hull[i+1][j  ].x * uv1 + 
			control_hull[i+1][j+1].x * uv2 + 
			control_hull[i  ][j+1].x * uv3;

	p->y =  control_hull[i  ][j  ].y * uv0 + 
			control_hull[i+1][j  ].y * uv1 + 
			control_hull[i+1][j+1].y * uv2 + 
			control_hull[i  ][j+1].y * uv3;

	p->z =  control_hull[i  ][j  ].z * uv0 + 
			control_hull[i+1][j  ].z * uv1 + 
			control_hull[i+1][j+1].z * uv2 + 
			control_hull[i  ][j+1].z * uv3;
}

	/* evaluate a spline at a particular (u,v) values. This function expects

		A SplineBasis array, which is a pre-computed combination of the
		16 control points combined with the 4x4 Basis Function.

		This function returns the point, and the derivatives of that point
			(one with respect to u, the other with respect to v), for the
			given parameters.
	*/

spline_eval(SplineBasis *b,float u,float v,Point3 *p,WORD dflag)
{	WORD	j;
	float	u2 = u*u,						/* u squared						*/
			u3 = u2*u,						/* u cubed							*/
			v2 = v*v,						/* v squared						*/
			v3 = v2*v,						/* v cubed							*/
			du2 = u * 2.0,					/* derivative of u squared			*/
			du3 = u*u*3.0,					/* derivative of u cubed			*/
			dv2 = v * 2.0,					/* derivative of v squared			*/
			dv3 = v*v*3.0;					/* derivative of v cubed			*/
	Point3	*du = p + 1,					/* derivatives						*/
			*dv = p + 2;

	p->x = p->y = p->z = 0.0;				/* set point to zero				*/
	*du = *p;								/* set derivatives to zero			*/
	*dv = *p;

	for (j=0; j<4; j++)						/* iterate through Y control points */
	{	Point3			sum,
						dsum;
		Point3			*cpoint;
		float			basis;

			/* calculate each of the four u-splines */
			/* also calculate the derivative with respect to u */

		cpoint = &(*b)[j][3];				/* first row * u^0					*/

		sum.x = cpoint->x;
		sum.y = cpoint->y;
		sum.z = cpoint->z;

		cpoint = &(*b)[j][2];				/* second row * u^1					*/

		sum.x += u * cpoint->x;
		sum.y += u * cpoint->y;
		sum.z += u * cpoint->z;

		dsum.x = cpoint->x;
		dsum.y = cpoint->y;
		dsum.z = cpoint->z;

		cpoint = &(*b)[j][1];				/* third row * u^2					*/

		sum.x += u2 * cpoint->x;
		sum.y += u2 * cpoint->y;
		sum.z += u2 * cpoint->z;

		dsum.x += du2 * cpoint->x;
		dsum.y += du2 * cpoint->y;
		dsum.z += du2 * cpoint->z;

		cpoint = &(*b)[j][0];				/* fourth row * u^3					*/

		sum.x += u3 * cpoint->x;
		sum.y += u3 * cpoint->y;
		sum.z += u3 * cpoint->z;

		dsum.x += du3 * cpoint->x;
		dsum.y += du3 * cpoint->y;
		dsum.z += du3 * cpoint->z;

#if 0
		cpoint = &(*b)[j][2];				/* second row * du/du (1.0)			*/
		cpoint = &(*b)[j][1];				/* third row * du^2/d2 (2u)			*/
		cpoint = &(*b)[j][0];				/* fourth row * du^3/du				*/
#endif
			/* calculate a single vspline that runs across the four u-splines */

		basis = (m[3][j])  +
				(m[2][j] * v) +
				(m[1][j] * v2) +
				(m[0][j] * v3);

		p->x += basis * sum.x;
		p->y += basis * sum.y;
		p->z += basis * sum.z;

		if (dflag)
		{	du->x += basis * dsum.x;
			du->y += basis * dsum.y;
			du->z += basis * dsum.z;

				/* calculate the derivative with respect to v, similarly */

			basis = (m[2][j]) +
					(m[1][j] * dv2) +
					(m[0][j] * dv3);

			dv->x += basis * sum.x;
			dv->y += basis * sum.y;
			dv->z += basis * sum.z;
		}
	}
}

void cross (Point3 *a, Point3 *b, Point3 *c)
{	c->x = a->y * b->z - b->y * a->z;
	c->y = a->z * b->x - b->z * a->x;
	c->z = a->x * b->y - b->x * a->y;
}

void vsubtract (Point3 *a, Point3 *b, Point3 *c)
{	c->x = a->x - b->x;
	c->y = a->y - b->y;
	c->z = a->z - b->z;
}

float length(Point3 *v)
{	float l;
	l = v->x*v->x + v->y*v->y + v->z*v->z;
	l = sqrt(l);
	return l;
}

void normalize(Point3 *v)
{	float l = length(v);
	v->x /= l;
	v->y /= l;
	v->z /= l;
}

#if 0
sub_patch(SplineBasis *b, float u1, float u2, float v1, float v2, WORD level)
{	float	um,vm;

	um = (u1 + u2) / 2.0;
	vm = (v1 + v2) / 2.0;

	if (level < 3)
	{	level++;
		sub_patch(b, u1,um,v1,vm,level);
		sub_patch(b, u1,um,vm,v2,level);
		sub_patch(b, um,u2,v1,vm,level);
		sub_patch(b, um,u2,vm,v2,level);
		level--;
	}
	else
	{	Point		s0,s1,s2,s3;
		WORD		icolor;
		Point3		snormal, du, dv;
		Point3V		p0,p1,p2,p3;

		spline_eval( b, u1, v1, &p0.p, TRUE );
		spline_eval( b, u2, v1, &p1.p, TRUE );
		spline_eval( b, u2, v2, &p2.p, TRUE );
		spline_eval( b, u1, v2, &p3.p, TRUE );

/*		interpolate_points( i,j, u1,v1, &p0 );
		interpolate_points( i,j, u2,v1, &p1 );
		interpolate_points( i,j, u2,v2, &p2 );
		interpolate_points( i,j, u1,v2, &p3 );
*/
		transform( &p0.p, &s0 );
		transform( &p1.p, &s1 );
		transform( &p2.p, &s2 );
		transform( &p3.p, &s3 );

		cross(&p3.du,&p3.dv,&snormal);
		normalize(&snormal);

		icolor = (long)(8.0 - snormal.x * 8.0);

		icolor = clamp(2,icolor,15);
		SetAPen(rp,icolor);

		AreaMove(rp,s0.x,s0.y);
		AreaDraw(rp,s1.x,s1.y);
		AreaDraw(rp,s2.x,s2.y);
		AreaDraw(rp,s3.x,s3.y);
		AreaEnd(rp);
	}
}
#endif

sub_patch(SplineBasis *b,
	float u1, float u2, float v1, float v2,
	Point3V *p0, Point3V *p1, Point3V *p2, Point3V *p3,
	WORD level)
{	float	um,vm;

	um = (u1 + u2) / 2.0;
	vm = (v1 + v2) / 2.0;

	if (level < 5)
	{	Point3V		p01,p12,p23,p30,pm;

		spline_eval( b, um, v1, &p01.p, TRUE );
		spline_eval( b, u2, vm, &p12.p, TRUE );
		spline_eval( b, um, v2, &p23.p, TRUE );
		spline_eval( b, u1, vm, &p30.p, TRUE );
		spline_eval( b, um, vm, &pm.p,  TRUE );

		level++;
		sub_patch(b, um,u2,vm,v2, &pm, &p12, p2, &p23, level);
		sub_patch(b, um,u2,v1,vm, &p01, p1, &p12,&pm,  level);
		sub_patch(b, u1,um,vm,v2, &p30,&pm, &p23, p3,  level);
		sub_patch(b, u1,um,v1,vm,  p0, &p01,&pm, &p30, level);
		level--;
	}
	else
	{	Point		s0,s1,s2,s3;
		WORD		icolor;
		Point3		snormal, du, dv;

/*		Point3V		p0,p1,p2,p3;

		spline_eval( b, u1, v1, &p0.p, TRUE );
		spline_eval( b, u2, v1, &p1.p, TRUE );
		spline_eval( b, u2, v2, &p2.p, TRUE );
		spline_eval( b, u1, v2, &p3.p, TRUE );
*/

/*		interpolate_points( i,j, u1,v1, &p0 );
		interpolate_points( i,j, u2,v1, &p1 );
		interpolate_points( i,j, u2,v2, &p2 );
		interpolate_points( i,j, u1,v2, &p3 );
*/
		transform( &p0->p, &s0 );
		transform( &p1->p, &s1 );
		transform( &p2->p, &s2 );
		transform( &p3->p, &s3 );

		cross(&p3->du,&p3->dv,&snormal);
		normalize(&snormal);

		icolor = (long)(8.0 - snormal.x * 8.0);

		icolor = clamp(2,icolor,15);
		SetAPen(rp,icolor);

		AreaMove(rp,s0.x,s0.y);
		AreaDraw(rp,s1.x,s1.y);
		AreaDraw(rp,s2.x,s2.y);
		AreaDraw(rp,s3.x,s3.y);
		AreaEnd(rp);
	}
}

#define ITERATIONS 8.0

do_patch(ibase,jbase) WORD ibase,jbase;
{	SplineBasis		PreCalc;
	WORD			i,j,k;

	for (j=0; j<4; j++)						/* iterate through Y control points */
	{	Point3		sum;

		for (k=0; k<4; k++)					/* iterate through X control points */
		{
			PreCalc[j][k].x = 0.0;
			PreCalc[j][k].y = 0.0;
			PreCalc[j][k].z = 0.0;

			for (i=0; i<4; i++)				/* iterate through powers of U		*/
			{	float basis;				/* weighting factor for averaging	*/
				Point3	*cpoint;

				basis = m[k][i];
				cpoint = &control_hull[i+ibase][j+jbase];

				PreCalc[j][k].x += basis * cpoint->x;
				PreCalc[j][k].y += basis * cpoint->y;
				PreCalc[j][k].z += basis * cpoint->z;
			}
		}
	}

#if 0
	for (i=ITERATIONS; i>0; i--)
	{	float u1,u2;

		u1 = (float)i / ITERATIONS;
		u2 = (float)(i-1) / ITERATIONS;

		for (j=ITERATIONS; j>0; j--)
		{	float		v1,v2;

			v1 = (float)j / ITERATIONS;
			v2 = (float)(j-1) / ITERATIONS;

			sub_patch(&PreCalc, u1, u2, v1, v2, 5);
		}
	}
#endif

/*	sub_patch(ibase,jbase, 0.0, 1.0, 0.0, 1.0, 0); */
/*	sub_patch(&PreCalc, 0.0, 1.0, 0.0, 1.0, 0); */

	{	Point3V	p0,p1,p2,p3;

		spline_eval( &PreCalc, 0.0, 0.0, &p0.p, TRUE );
		spline_eval( &PreCalc, 1.0, 0.0, &p1.p, TRUE );
		spline_eval( &PreCalc, 1.0, 1.0, &p2.p, TRUE );
		spline_eval( &PreCalc, 0.0, 1.0, &p3.p, TRUE );
		sub_patch(&PreCalc, 0.0, 1.0, 0.0, 1.0, &p0,&p1,&p2,&p3, 0);
	}
}