/*

   compile with:

  	cc -o viewhf viewhf.c -lgl_s


-----------------------------------------------------------------------------

  	Viewhf
 
	(View Height Field)
 	Copyright (C) 1992, Lawrence K. Coffin.
	Parts Copyright (C) 1989, 1991, Craig E. Kolb.
 	All rights reserved.
 
	Parts of this program are from "Numerical Recipes ic C" by William H.
	Press, Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling.

   This software may be freely copied, modified, and redistributed
	provided that this copyright notice is preserved on all copies.
 
   You may not distribute this software, in whole or in part, as part of
	any commercial product without the express consent of the authors.
  
   There is no warranty or other guarantee of fitness of this software
	for any purpose.  It is provided solely "as is".


   This program is for displaying height field files in the format as described
	in the Rayshade 4.0 users guide.  Viewhf displays the points
	in the height field file as a surface of dots in three dimensional
	space.  The points are colored from red to blue according to their
	altitude.  The surface is centered at the origin.  Viewhf allows you to
	zoom in and out from the origin, to rotate the surface around the
	z-axis, as well as to view it from the top (0 degrees), bottom (180
	degrees), or any angle inbetween.  A colored axis is included to
	help you keep track of your position, the arms of the axis extend from
	the origin out to 0.1 along each axis.  The surface extends from -0.5 to
	0.5 in both the x and y directions.

   Useage:

	viewhf [options] <height field file>

   Options:
	-h		help
	-x <xscale>	x scaling
	-y <yscale>	y scaling
	-z <zscale>	z scaling

	<xscale> is the scaling along the x axis in thousands of a unit.
		ie. 500 would be one half the original height.
	<yscale> is the scaling along the y axis in thousands of a unit.
	<zscale> is the scaling along the z axis in thousands of a unit.

	The program will open a 630 x 500 window centered on the origin.  It
	displays the surface and a colored axis.  The initial view is from along
	the negative y-axis.  Holding the first mouse button down over the
	window lets you rotate the scene.  Moving the mouse left or right from
	the vertical center line will cause the nearest size of the surface to
	rotate left or right (left is clockwize around the z-axis, right is
	counterclockwise.)  Moving the mouse up and down from the horizontal
	centerline will change the view from above the surface, at the bottom of
	the window, to below the surface, at the top of the window, keeping
	the eye at a fixed distance from the origin.  Thus, holding the left
	button down and moving the mouse along the horizontal centerline will
	rotate the surface with the eye piece in the x y plane.  Pressing the
	middle mouse down while the mouse is inside the window will zoom in and
	out from the origin, when the mouse is above the horizontal centerline
	the eye piece will zoom in, below the line will zoom out.

	The number of points should only be limited to available memory, and
	user patience.  I run it on a Personal Iris with 8M ram with 4225
	points at a nice (fast?) speed, when in increase it to 16641 points
	(the next step up on my fractal landscape generator), it slows down
	quite a bit, but it's surviveable.  At points > 16641 
	(((2**n)+1)**2 : n > 7) it CRAWWWWWWWLS on this machine.

	Points are colored in a liner scale from (255, 0, 0) to (0, 127, 255)
	with the red being the minimum altitude and the greenish-blue the
	maximum altitude. (greenish-blue was used because pure blue was too
	dark against the black background). The color scheme can easily be
	changed, just change the #defines below.
 

   Future options:
 

	Using a mesh instead of points to display the surface.
		(Perhaps polygons also, once I learn how to do lighting
		and zbuffering and depthcueing.  Its hard to distinguish
		contour features from pure white polygons. Hmm... maybe colored
		polygons would work.)

	The ability to move around the surface instead of just in and out from
	the origin, around the z-axis, and from top to bottom. The reason for
	the present set-up is that I am using polarview to set the eye piece
	location.
	
	Have the program generate its own fractal surface and write it to a
	height field file, possibly allowing the user to edit individual points
	before or after the fractalization.

	Display the exact coordinates of the eyepiece.

	Shadows, depthcueing, adjustable windows, more fun stuff.

   BUGS:

	For some reason, the "input focus" is constantly being sent to the
	program so that even when the mouse is not in the window, the program
	knows its location as well as whether the mouse buttons are pushed or
	not (maybe it's because I'm polling for the info, rather than queueing.)
	For now I tell it to ignore the mouse if it is not within the bounds
	of the window, but this still does not stop the program from changing
	the display if the mouse buttons are pressed in a window that is 
	on top of viewhf's window but within its boundaries.


*/

#include <stdio.h>
#include <gl.h>
#include <device.h>

#define CSCALE cscale = 255*(al[x][y]-minal)/(maxal-minal)
	/* cscale = f(altitude):0->255 alt:minal->maxal */

#define RGBCOLOR RGBcolor(255-cscale, cscale/2, cscale)

GraphicsGetMousePos(int *x, int *y) /*Copyright (C) 1989, 1991, Craig E. Kolb*/
{
	int xo, yo;

	getorigin(&xo, &yo);
	*x = getvaluator( MOUSEX ) - xo;
	*y = getvaluator( MOUSEY ) - yo;

}

GraphicsRedraw() /*Copyright (C) 1989, 1991, Craig E. Kolb*/
{
	Device dev;
	short val;

	while (qtest()) {
		dev = qread(&val);
		if (dev == REDRAW) {
			reshapeviewport();
			return TRUE;
		}
	}
	return FALSE;
}

Object makeaxis(int);
Object heightfield(float **,int,float,float,float);
float **matrix(int,int,int,int);
void free_matrix(float**,int,int,int,int);

main(int argc, char *argv[])
{
	int c, i, x, y, z, dev, value;
	float azim, inc, dist, dazim, dinc, ddist;
	int size, change, xmin, ymin, xpos, ypos;
	long *xorig, *yorig, *width, *height;
	float **al;
	float minal, maxal;
	FILE *hf;
	Object axis, field;
	extern char *optarg;
	extern int optind;
	int xsize, ysize, help;
	float xscale, yscale, zscale;
	char *infile;

	xsize = 630;
	ysize = 500;
	maxal = minal = 0;
	xscale = 1.0;
	yscale = 1.0;
	zscale = 1.0;
	help = 0;

	while((c = getopt(argc, argv, "hx:y:z:")) != EOF)
		switch(c) {
		case 'h': help = 1;
			break;
		case 'x': xscale = (float)atoi(optarg)/1000.0;
			break;
		case 'y': yscale = (float)atoi(optarg)/1000.0;
			break;
		case 'z': zscale = (float)atoi(optarg)/1000.0;
			break;
		case '?': help = 1;
	}
	if((infile = argv[optind]) == NULL && !(help)) {
		fprintf(stderr,"\nI need a file to open!\n");
		help = 1;
	}

	if(help){
		fprintf(stderr,"\nUsage: viewhf [options] <height field file>\n");
		fprintf(stderr,"\n");
		fprintf(stderr,"Options:\n");
		fprintf(stderr,"\t-h\t\thelp\n");
		fprintf(stderr,"\t-x <xscale>\tx scaling\n");
		fprintf(stderr,"\t-y <yscale>\ty scaling\n");
		fprintf(stderr,"\t-z <zscale>\tz scaling\n");
		fprintf(stderr,"\n");
		fprintf(stderr,"\t<xscale> is the scaling along the x axis in thousands of a unit.\n");
		fprintf(stderr,"\t\tie. 500 would be one half the original height.\n");
		fprintf(stderr,"\t<yscale> is the scaling along the y axis in thousands of a unit.\n");
		fprintf(stderr,"\t<zscale> is the scaling along the z axis in thousands of a unit.\n\n");
		exit(1);
	}
	
/*
 * Open and read the heightfield file
 */
	if((hf = fopen(infile, "r")) == NULL){
		fprintf(stderr,"Unable to open %s\n",infile);
		exit(1);
	}
	
	fread(&size,sizeof(int),1,hf);

/*
 * Allocate memory to hold the array
 */

	al = matrix(0, size, 0, size);

	for(y=0; y <= size - 1; y++){
		for(x=0; x <= size - 1; x++){
			fread(&al[x][y],sizeof(float),1,hf);
		}
	}
	fclose(hf);



	prefsize(xsize, ysize);
	winopen(argv[0]);
	perspective(900, 1.26, 0, 10000);
	viewport(0, xsize -1, 0, ysize -1);
	color(BLACK);
	clear();
	RGBmode();
	doublebuffer();
	gconfig();
	unqdevice(INPUTCHANGE);
	qdevice(LEFTMOUSE);  /* Pop a square on request */
	qdevice(MIDDLEMOUSE);
	qdevice(RIGHTMOUSE);
	qdevice(REDRAW);

/*
 * Create the heightfield surface
 */

	field = heightfield(al,size,xscale,yscale,zscale);

/*
 * Free the matrix memory
 */

	free_matrix(al,0, size, 0, size);
fprintf(stderr,"");	/* for some reason I get an Object ID error
				without this here/*

/*
 * Create axis
 */
	axis = makeaxis(1);

	cpack(0x0);
	clear();
	swapbuffers();
	cpack(0x0);
	clear();
	dist = 1;
	azim = 0;
	inc = 900;

/*
 * Draw the initial view
 */

	pushmatrix();
		polarview(dist, azim, inc, 0);
		pushmatrix();
			translate(-.5*xscale, -.5*yscale, 0);
			callobj(field);
		popmatrix();
		pushmatrix();
			scale(.1, .1, .1);
			callobj(axis);
		popmatrix();
	popmatrix();
	swapbuffers();

/*
 * Draw and manipulate the objects
 */
	
    for(;;){
		
	if(getbutton(LEFTMOUSE)){
		GraphicsGetMousePos(&xpos, &ypos);
		if(xpos > 0 && xpos < 630 && ypos > 0 && ypos < 500){
			change = 1;
			dazim = -((float)xpos - 315.0)/10.0;
			dinc = ((float)ypos - 20.0);
			azim = azim + dazim;
			if(azim > 3600.0) azim = (azim - 3600.0);
			if(azim < 0) azim = (azim + 3600.0);
			inc = 1800.0*dinc/460.0;
			if(inc > 1800.0) inc = 1800.0;
			if(inc < 0.0) inc = 0.0;
		}
	}
	else if(getbutton(MIDDLEMOUSE)){
		GraphicsGetMousePos(&xpos, &ypos);
		if(xpos > 0 && xpos < 630 && ypos > 0 && ypos < 500){
			change = 1;
			ddist = -((float)ypos - 250)/4000;
			dist = dist + ddist ;
			if(dist > 10) dist = 10;
			if(dist < 0) dist = 0;
		}
	}
	else if(GraphicsRedraw()) {
		change = 1;
	}
	else {
		xpos = 0;
		ypos = 0;
		change = 0;
	}

	if(change){

		cpack(0x0);
		clear();
		pushmatrix();
			polarview(dist, azim, inc, 0);
			pushmatrix();
				translate(-.5*xscale, -.5*yscale, 0);
				callobj(field);
			popmatrix();
			pushmatrix();
				scale(.1, .1, .1);
				callobj(axis);
			popmatrix();
			swapbuffers();
		popmatrix();
	}

    }

}

Object makeaxis(int width)
{
	Object axis;
	makeobj(axis);
		pushattributes();
		setlinestyle(0);
		linewidth(width);
		RGBcolor(255, 0, 0);
		movei(0, 0, 0);
		drawi(1, 0, 0);
		RGBcolor(0, 255, 0);
		movei(0, 0, 0);
		drawi(0, 1, 0);
		RGBcolor(0, 0, 255);
		movei(0, 0, 0);
		drawi(0, 0, 1);
		RGBcolor(255, 255, 255);
		cmov(1, 0, 0);
		charstr("X");
		RGBcolor(255, 255, 255);
		cmov(0, 1, 0);
		charstr("Y");
		RGBcolor(255, 255, 255);
		cmov(0, 0, 1);
		charstr("Z");
		popattributes();
	closeobj();
	return axis;
}
Object heightfield(float **al,int size,float xscale,float yscale,float zscale)
{
	int i, x, y;
	float xcoor, ycoor;
	Object field;
	float minal, maxal;
	int cscale;

	minal = maxal = 0;

	for(y=0; y <= size - 1; y++){
		for(x=0; x <= size - 1; x++){
			if(al[x][y] < minal) minal = al[x][y];
			else if(al[x][y] > maxal) maxal = al[x][y];
		}
	}
	i = 0;

	makeobj(field);
	for(y=0; y <= size - 1; y++){
		for(x = 0; x <= size - 1; x++){
			CSCALE;
			RGBCOLOR;

			xcoor = ((float)(i%size)/(float)(size-1));
			ycoor = ((float)(i/size)/(float)(size-1));
			pnt((xcoor*xscale),(ycoor*yscale),(al[x][y]*zscale));
			i++;
		}
	}
	closeobj();

	return field;

}

/*
 *	The following are from "Numerical Recipes in C"
 */

/*
 * Generate an error message if error in allocating matrix
 */

void nrerror(error_text)
char error_text[];
{

	fprintf(stderr,"Numerical Recipes run-time error...\n");
	fprintf(stderr,"%s\n",error_text);
	fprintf(stderr,"...now exiting to system...\n");
	exit(1);
}

/*
 * Allocate memory for a matrix
 */

float **matrix(nrl,nrh,ncl,nch)
int nrl,nrh,ncl,nch;
{
	int i;
	float **m;

	m=(float **) malloc((unsigned) (nrh-nrl+1)*sizeof(float*));
	if (!m) nrerror("allocation failure 1 in matrix()");
	m -= nrl;

	for(i=nrl;i<=nrh;i++) {
		m[i]=(float *) malloc((unsigned) (nch-ncl+1)*sizeof(float));
		if (!m[i]) nrerror("allocation failure 2 in matrix()");
		m[i] -= ncl;
	}
	return m;
}

/*
 * Free memory used by a matrix
 */

void free_matrix(m,nrl,nrh,ncl,nch)
float **m;
int nrl,nrh,ncl,nch;
{
	int i;

	for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
	free((char*) (m+nrl));
}