#include <stdio.h>
#include <math.h>
#include <fcntl.h>
#include <sys/types.h>
#include <unistd.h>
#include <getopt.h>
#include <string.h>

#define ETOPO_XDIM 4320
#define ETOPO_YDIM 2160
#define ETOPO_ZMAX 10376    /* Max of absolute value of elevation */

#define SWAPINT(a,b)	{ int tmp; tmp=a; a=b; b=tmp; }

void usage(char *av0)
{
 fprintf(stderr,"Usage: %s [options] etopo5.bin lat0 lat1 lon0 lon1 out.wrl\n"
	"    Options: -z z-exaggeration\n"
	"	      -r resolution-scale\n"
	"	      -t texmap-name\n"
	"             -d dlat dlon dz\n"
	"	      -s scalefactor\n",av0);
 exit(-1);
}

void check_lat_lon(float lat0,float lat1,float lon0,float lon1)
{
 if ((lat0<-90.0) || (lat0>90.0) || (lat1<-90.0) || (lat1>90.0))
	{
	fprintf(stderr,"Latitude must be in the range [-90,90]\n");
	exit(1);
	}
#if 0
 if ((lon0<-180.0) || (lon0>180.0) || (lon1<-180.0) || (lon1>180.0))
	{
	fprintf(stderr,"Longitude must be in the range [-180,180]\n");
	exit(1);
	}
#endif
}

short indata[ETOPO_XDIM];

void read_row(int infd,int row,short *rowbuf,int left,int width,int resscale)
{
 int i,j;
 lseek(infd,row*ETOPO_XDIM*sizeof(short),SEEK_SET);
 read(infd,indata,sizeof(indata));
 for (i=0, j=left; (i<width) && (j<ETOPO_XDIM); i++, j+=resscale)
	rowbuf[i] = indata[j];
 for (j=0; i<width; i++, j+=resscale)
	rowbuf[i] = indata[j];
}

short **read_topo(char *fname,float lat0,float lat1,float lon0,float lon1,
		int *rwidth,int *rheight,int resscale)
{
 int infd,row0,row1,col0,col1,width,height,j;
 short **data;
 if ((infd = open(fname,O_RDONLY)) == -1)
	{
	perror(fname);
	exit(1);
	}
 row0 = (90.0-lat0)*(ETOPO_YDIM-1)/180.0 + 0.5;
 row1 = (90.0-lat1)*(ETOPO_YDIM-1)/180.0 + 0.5;
 if (row0 > row1) SWAPINT(row0,row1);
 height = row1-row0+1;
 height /= resscale;
 if (lon0<0) col0 = (lon0+360.0)*(ETOPO_XDIM-1)/360.0 + 0.5;
 else col0 = lon0*(ETOPO_XDIM-1)/360.0 + 0.5;
 if (lon1<0) col1 = (lon1+360.0)*(ETOPO_XDIM-1)/360.0 + 0.5;
 else col1 = lon1*(ETOPO_XDIM-1)/360.0 + 0.5;
 if (col1 > col0) width = col1-col0+1;
 else width = ETOPO_XDIM-col0 + col1+1;
 width /= resscale;
 data = (short **) malloc(height*sizeof(short *));
 for (j=0; j<height; j++)
	{
	data[j] = (short *) malloc(width*sizeof(short));
	read_row(infd,j*resscale+row0,data[j],col0,width,resscale);
	}
 close(infd);
 *rwidth = width;
 *rheight = height;
 return(data);
}

#include "dave.h"    /* See note below */
int read_colormap(char *fname,float **r,float **g,float **b,int *numcolors)
{
 UCHAR *cmap;
 int xdim,ydim,zdim,i;
/* NB: This function (dpReadSgiImageByte) is defined in dave.h & libdave.a,
   which can be found at http://www.evl.uic.edu/pape/sw/libdave.html.
   It also requires SGI's libimage.  If you're not using an SGI, replace
   it with your own colormap function. */
 if (!dpReadSgiImageByte(fname,&cmap,&xdim,&ydim)) return(0);
 *r = (float *) malloc(ydim * sizeof(float));
 *g = (float *) malloc(ydim * sizeof(float));
 *b = (float *) malloc(ydim * sizeof(float));
 for (i=0; i<ydim; i++)
	{
	(*r)[i] = cmap[i*xdim+1]/255.0;
	(*g)[i] = cmap[i*xdim+2]/255.0;
	(*b)[i] = cmap[i*xdim+3]/255.0;
	}
 free(cmap);
 *numcolors = ydim;
 return 1;
}

void write_header(FILE *outfp)
{
 fprintf(outfp,"#VRML V2.0 utf8\n\n");
 fprintf(outfp,"Group {\n");
 fprintf(outfp,"children [");
 fprintf(outfp,"DirectionalLight { direction 0 -1 0 }\n");
 fprintf(outfp,"Shape {\n");
}

void write_appearance(FILE *outfp,char *texmap)
{
 fprintf(outfp,"appearance Appearance { material Material { diffuseColor .7 .7 .7 }");
 if (texmap)
	{
	fprintf(outfp," texture ImageTexture { url \"%s\" }",texmap);
	fprintf(outfp," textureTransform TextureTransform { scale 1 -1 }");
	}
 fprintf(outfp,"}\n");
}

void write_grid(FILE *outfp,short **topo,int width,int height,float zscale,
		int resscale,float scalefactor,int writeColors)
{
 float dx,dy,heightScale = zscale*scalefactor;
 int i,j;
 dx = dy = 40000000.0/ETOPO_XDIM * resscale;
 fprintf(outfp,"geometry ElevationGrid {\n"
		"	xDimension %d\n"
		"	xSpacing %f\n"
		"	zDimension %d\n"
		"	zSpacing %f\n",
	width,dx*scalefactor,height,dy*scalefactor);
 fprintf(outfp,"height [\n");
 for (j=0; j<height; j++)
    {
    for (i=0; i<width; i++)
	fprintf(outfp,"%.1f\n",topo[j][i]*heightScale);
    }
 fprintf(outfp,"]\n");
 if (writeColors)
	{
	float *r,*g,*b;
	int numcolors,index;
	if (!read_colormap("etopo.M",&r,&g,&b,&numcolors))
		;	/* Put dummy cmap here */
	fprintf(outfp,"color Color { color [\n");
	for (j=0; j<height; j++)
		{
		for (i=0; i<width; i++)
			{
			if (topo[j][i]>0)
				index = topo[j][i]*(numcolors/2-1)/ETOPO_ZMAX+numcolors/2;
			else
				index = (topo[j][i]+ETOPO_ZMAX)*(numcolors/2-1)/ETOPO_ZMAX;
			fprintf(outfp,"%.3f %.3f %.3f,\n",r[index],g[index],b[index]);
			}
		}
	fprintf(outfp,"] }\n");
	}
 fprintf(outfp,"}\n");
}


main(int argc,char **argv)
{
 float lat0,lat1,lon0,lon1,zscale=1.0,scalefactor=1.0,step;
 FILE *outfp;
 short **topo;
 int width,height,opt,tx,ty,resscale=1;
 char *texmap=NULL;
 while ((opt = getopt(argc,argv,"z:r:t:s:")) != -1)
	switch (opt)
		{
		case 'z': zscale = atof(optarg); break;
		case 'r': resscale = atoi(optarg); break;
		case 't': texmap = strdup(optarg); break;
		case 's': scalefactor = atof(optarg); break;
		}
 if ((argc-optind<6) || (resscale < 1))
	usage(argv[0]);
 lat0 = atof(argv[optind+1]);
 lat1 = atof(argv[optind+2]);
 lon0 = atof(argv[optind+3]);
 lon1 = atof(argv[optind+4]);
 check_lat_lon(lat0,lat1,lon0,lon1);
 topo = read_topo(argv[optind],lat0,lat1,lon0,lon1,&width,&height,resscale);
 if ((outfp = fopen(argv[optind+5],"w")) == NULL)
	{
	perror(argv[optind+5]);
	exit(1);
	}
 write_header(outfp);
 write_appearance(outfp,texmap);
 write_grid(outfp,topo,width,height,zscale,resscale,scalefactor,(texmap==NULL));
 fprintf(outfp,"}\n");  /* Close the Shape node */
 fprintf(outfp,"]\n}\n");  /* Close the Group node */
 step = 40000000.0/ETOPO_XDIM * resscale * scalefactor;
 fprintf(outfp,"Viewpoint { position %f %f %f  orientation 1 0 0 -1.57"
	" description \"Above\" }\n",
	width*step/2.0,width*step,height*step/2.0);
 fprintf(outfp,"Viewpoint { position %f %f %f description \"Center\" }\n",
	width*step/2.0,1000.0*scalefactor,height*step/2.0);
 fprintf(outfp,"NavigationInfo { speed %f type \"EXAMINE\" }\n",100000.0*scalefactor);
 fclose(outfp);
}