/*

    compile with:

  	cc -o cvhf cvhf.c -lgl_s

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

  	Cvhf
 
	(Contour 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 in 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".

   Useage:

	cvhf <height field file> <spacing> <window size>


		<height field file> is the name of the height field file.

		<spacing> is an integer that represents the change in
			altitude in thousands of a unit (ie. a spacing of
			100 will put lines at ...-0.2, -0.1, 0, 0.1, 0.2...)

		<window size> is the size of the window in pixels

	The program will open a square window and displays the altitude of
	each pixel as a color.  Red is the lowest, blue the highest with a
	gradient in between.  Pixels that fall on the specified spacing are
	white.  When the left mouse button is pressed within the window,
	a box will be shown that displays the x and y coordinates and the
	altitude at that point.  A small cross-hair is placed in the window
	centered on the pixel whose coordinates are in the box.  The
	cross-hair and box will follow the mouse around as the mouse is moved
	around the window with the left button pressed.  When the left mouse
	button is released, the last displayed coordinates will remain on the
	screen.
 
	I have tried cvhf with heightfield files with up to 1050625
	(((2**10)+1)**2) points with no noticeable problems.  However, I
	have found that if the height field file is made up of numerous
	whole integer altitudes, the points for the contour lines 
	do not form continuous lines; the color scheme holds up, as does
	the data for pixels picked by the mouse.


   This program is for displaying height field files in the format as described
	in the Rayshade 4.0 users guide.  Cvhf draws a contour map of the
	surface defined by the heightfield file.  Cvhf also colors in the
	areas between the contour lines with a shade between red and blue
	that corresponds to the altitude of that point. Red is the lowest
	points, blue is the highest.  The method by which cvhf draws the contour
	lines is probably pretty crude; I thought it up in a few hours and got
	rather nice results.  Cvhf scans each pixel and sends its x and y
	coordinates to a separate routine along with the array of height field
	values and the heightfield size.  The routine calculates and returns
	the altitude of the point the the main program.  If the
	((altitude*1000) or (altitude*1000)-1)%<spacing> is zero the pixel is
	colored white, else the pixel is colored from red to blue according to
	its altitude.

	The subroutine that calculates the altitude of a point, (represented
	by p7 in the figure below) does so by taking the altitudes of the
	four closest points defined in the heightfield, calculating the points
	on the lines from p1->p2 (point p5) and p3->p4 (point p6) that have
	the same x value as p7.  It then calculates the height of the line from
	p5->p6 where it passes through p7.  Note that the same altitude is
	obtained using p5~ and p6~ instead of p5 and p6.

	p3	   p6		     p4
	*----------*-----------------*
	|	   |		     |
	|	   |		     |
	|	   |		     |
	|	   |		     |
	|	   |		     |
	|	   |		     |
	|	   |		     |
	|	   |p7		     |
    p5~	*----------*-----------------*p6~
	|	   |		     |
	|	   |		     |
	|	   |		     |
	|	   |		     |
	|	   |		     |
	|	   |		     |
	|	   |		     |
	*----------*-----------------*
	p1	   p5		     p2

	Please send all comments, suggestions, and bug reports to:

	lcoffin@clciris.chem.umr.edu

*/

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


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;
}

float alt(float **,int,int,int);
void free_matrix(float**,int,int,int,int);
float **matrix(int,int,int,int);
void nrerror(char);

main(int argc, char *argv[])
{
	int i, j, x, y, cscale;
	int size, change, oldstry, oldstrx;
	float **al;
	FILE *hf;
	Object grid;
	int *color, xpos, ypos, strx, stry, mouse;
	int inwinx, inwiny;
	char coords[100];
	int scale;
	int xsize, ysize, dev, value;
	float altn, minal, maxal;
	scale = atoi(argv[3]);

	color = (int *)malloc((unsigned)((scale+1)*(scale+1)*sizeof(int)));

	if(argc < 4){
		printf("\nUsage: %s <file> <spacing> <window size>\n\n",argv[0]);
		printf("\t<file> is the name of the heightfield file.\n");
		printf("\t<spacing> is the change in altitude, in thousands of a unit,\n"); 
		printf("\t\tbetween contour lines\n");
		printf("\n\tie.a spacing of 100 will space the lines 0.1 units apart.\n");
		exit(1);
	}
/*
 * Open hieght field file and read size
 */
	if((hf = fopen(argv[1], "r")) == NULL){
		fprintf(stderr,"Unable to open %s\n",argv[1]);
		exit(1);
	}

	fread(&size,sizeof(int),1,hf);

/*
 * Allocate memory the array for the heightfield data and read the data
 */

	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);
/*
 * Find min and max altitude
 */

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


/*
 * Generate altitudes and draw points
 */

	for(y = 0; y < (scale+1); y++){
		for(x = 0; x < (scale+1); x++){
			altn = alt(al,size,scale,x,y);
			if(((int)(1000*altn)%atoi(argv[2]))==0 ||
				(((int)(1000*altn)-1)%atoi(argv[2]))==0 ||
				(((int)(1000*altn)+1)%atoi(argv[2]))==0){
				color[y*(scale+1)+x] = 0xffffff;
			}
			else {
				cscale = 255*(altn - minal)/(maxal - minal);
				color[y*(scale+1)+x] = (255-cscale)+cscale*0x10000;
			}
		}
	}

/*
 * Open window and clear background
 */
	xsize = (scale+1);
	ysize = (scale+1);

	prefsize(xsize, ysize);
	winopen(argv[0]);
	ortho2(0,scale,0,scale);
	viewport(0, xsize -1, 0, ysize -1);
	RGBmode();
	doublebuffer();
	gconfig();

	setcross();

	qdevice(LEFTMOUSE);

	RGBcolor(0,0,0);
	clear();


	lrectwrite(0,0,scale,scale,color);
	swapbuffers();
	lrectwrite(0,0,scale,scale,color);


	stry = 20;	/* Initial values so that nothing bad happens */
	strx = 20;
	xpos = scale/2;
	ypos = scale/2;

	while(1){
		if(getbutton(LEFTMOUSE)){
			GraphicsGetMousePos(&xpos, &ypos);
			if(xpos > 0 && xpos < (scale+1) && ypos > 0 && ypos < (scale+1)){
				change = 1;
				mouse = 1;
				altn = alt(al,size,scale,xpos,ypos);
				sprintf(coords,"(%2.3f, %2.3f, %2.3f)",(float)xpos/scale,(float)ypos/scale,altn);
				oldstrx = strx;
				oldstry = stry;
				if(xpos < scale/2) strx = xpos + 15;
				if(xpos >= scale/2) strx = xpos-(strwidth(coords)+10);
				if(ypos < scale/2) stry = ypos+10;
				if(ypos >= scale/2) stry = ypos-15;

			}
		}
		else if(GraphicsRedraw()){
			change = 1;
		}
		else {
			change = 0;
			GraphicsGetMousePos(&inwinx, &inwiny);
			if(inwinx > 0 && inwinx < (scale+1) && inwiny > 0 && inwiny < (scale+1)){
				setcross();
			}
		}

		if(change){
			lrectwrite(0,0,scale,scale,color);
fprintf(stderr,"");/* Won't display the following if this is not here! */
			if(mouse){

				pushattributes();
					cpack(0x0);
					rectf(strx-2,stry-5,strx+strwidth(coords)+1,stry+14);

					linewidth(1);
					cpack(0xffffff);
					rect(strx-1,stry-4,strx+strwidth(coords),stry+13);
fprintf(stderr,"");/* Won't display the following if this is not here! */
					cpack(0x0000);
					move(xpos-5,ypos);
					draw(xpos-2,ypos);
					move(xpos+2,ypos);
					draw(xpos+5,ypos);
					move(xpos,ypos-5);
					draw(xpos,ypos-2);
					move(xpos,ypos+2);
					draw(xpos,ypos+5);
				popattributes();
				
				cpack(0xffff00);
				cmov(strx,stry);
				charstr(coords);

			}
			swapbuffers();
		}

	}
}

float alt(float **al,int size,int scale, int x, int y)
{
	float alt1, alt2, alt3, alt4, alt5, alt6, altn;
	float pntspace;
	int x1, x2, y1, y2, xn, yn;
	int i1, i2, j1, j2;

	pntspace = (scale*1.0/(float)(size-1));

	i1 = (int)(x/pntspace);
	i2 = i1 + 1;
	j1 = (int)(y/pntspace);
	j2 = j1 + 1;



	x1 = (int)(i1*pntspace);
	x2 = (int)(i2*pntspace);
	y1 = (int)(j1*pntspace);
	y2 = (int)(j2*pntspace);
	xn = x;
	yn = y;

	

	alt1 = al[i1][j1];
	alt2 = al[i2][j1];
	alt3 = al[i1][j2];
	alt4 = al[i2][j2];
/*
 * Here's the magic formula!!
 */

	altn = ((((alt4-alt3-alt2+alt1)*(xn-x1)/pntspace) + alt3 -
		alt1)*(yn - y1)/pntspace) + ((alt2 - alt1)*(xn-
		x1)/pntspace) + alt1;

	return altn;
}


setcross()
{
	static short cross[16] = {
		0x0100,	0x0100,	0x0100,	0x0100,
		0x0100,	0x0100,	0,	0xfc7e,
		0,	0x0100,	0x0100,	0x0100,
		0x0100,	0x0100,	0x0100,	0
	};


	curstype(C16X1);
	drawmode(CURSORDRAW);
	mapcolor(1,255,255,0);
	defcursor(1,cross);
	curorigin(1,7,7);
	setcursor(1,0,0);
	drawmode(NORMALDRAW);
}
		

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

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);
}

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;
}

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));
}