/*
 * Usage:
 * domap [file1] [file2] [file3] [...] [options] | pen xcenter=0 ycenter=0
 *						   ^^^
 *						   A Vplot filter
 *
 * If project=round, this program projects the Earth onto a plane tangent
 * to its surface, as seen from your "eye", where the position of your eye
 * is given by the following parameters:
 * 
 * If project=round
 *    lat=37.45       <- Latitude of eye
 *    long=-122.21    <- Longitude of eye
 *    dist=100.       <- Distance of eye above Earth, in miles
 *    inch=110.	      <- Miles per inch at center of projection
 *
 * If project=flat, a crude latitude-longitude plot is made instead.
 *
 * If project=flat
 *    long=-122.21    <- Longitude centered in the plot
 *
 * file1, file2, file3, etc, are map information files
 * such as are in ./Namer.Map/*.cbd, etc.
 *
 * A slew of other options are used to control what sorts of objects
 * get plotted, and in what colors. These are:
 * 
 * pby1=6   <- Vplot color to make this type of object; 0 means omit
 * pby2=6
 * pby3=5
 * 
 * bdy1=7
 * bdy2=7
 * bdy3=6
 * 
 * riv1=5
 * riv2=1
 * riv3=1
 * riv4=1
 * riv5=5
 * riv6=1
 * riv7=1
 * riv8=1
 * riv9=1
 * riv10=1
 * riv11=1
 * riv12=1
 * 
 * cil1=7
 * cil2=4
 * cil3=1
 * cil4=1
 * cil5=1
 * cil6=1
 * cil7=1
 * cil8=1
 * cil9=1
 * cil10=1
 * cil11=1
 * cil12=1
 * cil13=1
 * cil14=1
 * cil15=1
 *
 * The meanings of these different objects are defined in the
 * "Database_info" file.
 *
 * To read the data files properly, you may need to rewrite the
 * "twiddle32" and "twiddle16" subroutines. The data files are
 * integers, and so proper reading of them depends on the byte
 * ordering of your machine.
 *
 * If you don't have vplot, you will have to convert the "vp_*" routines
 * into appropriate equivalents. Here are examples of the various routines
 * used by domap:
 *
 *  vp_filep (stdout);
 *		Where the output file is to go. Here "standard out".
 *
 *  vp_erase (stdout);
 *		Clear the page.
 *
 *  vp_scale (1. / viewscale, 1. / viewscale);
 *		Sets a global scale factor for the plotting coordinates.
 *		(Here just divides all coordinates by "viewscale" before
 *		using them.)
 *
 *  vp_color (WHITE);
 *		Sets the output plotting color.
 *
 *  vp_penup ();
 *		Lift the pen.... next line-drawing command will be a move.
 *
 *  vp_upendn (horizon * cos (j * CONV / 10.), horizon * sin (j * CONV / 10.));
 *		Draw a line to this coordinate, if pen currently down.
 *		If pen currently up, move to this coordinate.
 *
 *  vp_fat (1);
 *		Set the current line-drawing fatness.
 *
 *  vp_move (-5., -2.5);
 *		Move current location to this coordinate.
 *  vp_draw (5., -2.5);
 *		Draw line from current location to this coordinate.
 */

#include <windows.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <fcntl.h>
#include <io.h>
#include <math.h>
#include <winlib/winlib.h>

/*
 * Much of this program was written by Brian Reid at the Dec
 * Western Research Lab. The rest was written by Joe Dellinger
 * at Stanford University. Neither of us much cares what becomes
 * of it. Sat Apr  1 01:26:51 PST 1989
 *
 * Slightly revised by Joe Dellinger
 * Fri Feb 12 21:57:56 HST 1993
 */

/* ----------------------------------------------------------------------
 * Read and process one binary map file.
 */


#define SHORTFLAG 0x4000
#define CBD_MAGIC 0x20770002
#define BIT32 int
#define BIT16 short

#define CONV (3.141592654 / 180.)

/*
 * Input parameters:
 * viewlat, viewlong = lat, long in center of picture;
 * viewdist = distance of eye from earth, in miles;
 * viewscale = miles per inch at center of picture
 *
 * For project=round, all of these are ignored except
 * viewlong.
 */
double viewlat, viewlong, viewdist, viewscale;
char projection[80];

/*
 * for getpar
 */
int xargc;
char **xargv;

struct cbdhead
{
    BIT32           magic;	/* Magic number */
    BIT32           dictaddr;	/* Offset of segment dictionary in file */
    BIT32           segcount;	/* Number of segments in file */
    BIT32           segsize;	/* Size of segment dictionary */
    BIT32           segmax;	/* Size of largest segment's strokes, /2 */
    BIT32           fill[(40 / sizeof (BIT32)) - 5];	/* Filler */
};


double my_acos(double x)
{
	if (fabs (x) >= 1.)
		return 0.;
	else
		return acos (x);
}


HWND hwndMain, hwndStatus;
HDC hdc;
int penUp = TRUE;


/*
 *  vp_filep (stdout);
 *		Where the output file is to go. Here "standard out".
 */

void vp_filep(FILE *fd)
{

}


/*
 *  vp_erase ();
 *		Clear the page.
 */

void vp_erase(void)
{

}


/*
 *  vp_scale (1. / viewscale, 1. / viewscale);
 *		Sets a global scale factor for the plotting coordinates.
 *		(Here just divides all coordinates by "viewscale" before
 *		using them.)
 */

void vp_scale(double x, double y)
{

}


/*
 *  vp_color (WHITE);
 *		Sets the output plotting color.
 */
 
void vp_color(COLORREF color)
{
	static HBRUSH hBrush = NULL;
	static HPEN hPen = NULL;

	if (hBrush) {
		DeleteObject(SelectObject(hdc, hBrush));
		DeleteObject(SelectObject(hdc, hPen));
	}

	hBrush = SelectObject(hdc, CreateSolidBrush(color));
	hPen = SelectObject(hdc, CreatePen(PS_SOLID, 0, color));
}


/*
 *  vp_penup ();
 *		Lift the pen.... next line-drawing command will be a move.
 */

void vp_penup (void)
{
	penUp = TRUE;
}


/* 
 * vp_upendn (horizon * cos (j * CONV / 10.), horizon * sin (j * CONV / 10.));
 *		Draw a line to this coordinate, if pen currently down.
 *		If pen currently up, move to this coordinate.
 */

void vp_upendn(double x, double y)
{
	POINT old;

	if (penUp) {
		MoveToEx(hdc, (int)x, (int)y, &old);
		penUp = FALSE;
	} else
		LineTo(hdc, (int)x, (int)y);
}


/*
 *  vp_fat (1);
 *		Set the current line-drawing fatness.
 */
 
void vp_fat(int n)
{

}


/*
 *  vp_move (-5., -2.5);
 *		Move current location to this coordinate.
 */
 
void vp_move(double x, double y)
{
	POINT old;

	penUp = TRUE;
	MoveToEx(hdc, (int)x, (int)y, &old);
}


/*
 *  vp_draw (5., -2.5);
 *		Draw line from current location to this coordinate.
 */
 
void vp_draw(double x, double y)
{
	penUp = FALSE;
	LineTo(hdc, (int)x, (int)y);
}


void project(double lati, double longi)
{
	double xx, yy, zz;
	double x, y, z;
	double dot, factor;
	static double xlast = 0.;
	char tmp[80];
	static int n = 1;

	if (projection[0] == 'r') {
		x = sin (longi - viewlong) * cos (lati);
		y = cos (longi - viewlong) * cos (lati);
		z = sin (lati);

		yy = y * cos (viewlat) + z * sin (viewlat);
		zz = -y * sin (viewlat) + z * cos (viewlat);
		xx = x;

		dot = xx * xx + zz * zz - yy * (viewdist + (1. - yy));

		if (dot > 0.) {
			vp_penup ();
			return;
		}

		factor = viewdist / (viewdist + (1. - yy));
		xx *= factor;
		zz *= factor;

		vp_upendn (xx, zz);
	} else {
		x = longi - viewlong;
		while (x >= 180 * CONV)
			x -= 360 * CONV;
		while (x < -180 * CONV)
			x += 360 * CONV;
		y = lati;
		
		if (fabs (xlast - x) > 180. * CONV)
			vp_penup ();
		xlast = x;
/*
		sprintf(tmp, "%f %f", x, y);
		ShowStatus(hwndStatus, 0, tmp);
*/
		vp_upendn (x, y);
	}
}


LRESULT CALLBACK MainWndProc(HWND hwnd, UINT msg, WPARAM wParam,
		LPARAM lParam)
{
	switch (msg) {
		case WM_DESTROY:
			PostQuitMessage(0);
			return 0;
	}
	return DefWindowProc(hwnd,msg,wParam,lParam);
}


BOOL MakeWindow(HINSTANCE hInstance, int nCmdShow)
{
	WNDCLASS wc;
	
	wc.style = CS_HREDRAW | CS_VREDRAW;
	wc.lpfnWndProc = (WNDPROC)MainWndProc;
	wc.cbClsExtra = 0;
	wc.cbWndExtra = 0;
	wc.hInstance = hInstance;
	wc.hIcon = LoadIcon(hInstance, IDI_WINLOGO);
	wc.hCursor = LoadCursor(NULL, IDC_ARROW);
	wc.hbrBackground = (HBRUSH)(COLOR_WINDOW + 1);
	wc.lpszMenuName = NULL;
	wc.lpszClassName = "MainWndClass";
	if (! RegisterClass(&wc))
		return 0;

	hwndMain = CreateWindow(
		"MainWndClass",
		"DoMap",
		WS_OVERLAPPEDWINDOW,
		CW_USEDEFAULT, 0, CW_USEDEFAULT, 0,
		NULL, NULL, hInstance, NULL);
		
	hwndStatus = CreateStatusWnd(hInstance, hwndMain, 2);

	ShowWindow(hwndMain, nCmdShow);
	UpdateWindow(hwndMain);
	InvalidateRect(hwndStatus, NULL, TRUE);
	UpdateWindow(hwndStatus);
	return TRUE;
}


int WINAPI WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance,
		LPSTR lpCmdLine, INT nCmdShow)
{
	int Ifile, segcount, idx, idy, segbufsize, olt, oln, j, k, jstroke, iseg;
	COLORREF colfile[20];
	char string[80], string2[80], *sptr;
	char filename[20][80];
	int filenum;
	BIT32 i32;
	BIT16 *segbuf;
	double x, y, z, x1, x2, xc, y1, y2, yc, z1, z2, zc;
	double dist, cent, clipdist, horizon, my_acos ();
	double longi, lati;
	RECT rect;
	MSG msg;

	struct segdict {
		BIT32 segid;
		BIT32 maxlat, minlat, maxlong, minlong;
		BIT32 absaddr;
		BIT16 nbytes;
		BIT16 rank;
	} *sd, *sdbuf;

	struct segment {
		BIT32 orgx, orgy;
		BIT32 id;
		BIT16 nstrokes;
		BIT16 dummy;
	} sb;

	struct cbdhead  header;

#ifdef DEBUG
	int deg, min, sec, count;
#endif

	if (! MakeWindow(hInstance, nCmdShow))
		return 0;

    viewlat = 58.37;
    viewlong = 16.1;
    viewdist = 100000.;
    viewscale = 100.;

	viewlat *= CONV;
	viewlong *= CONV;
	viewdist /= 3950.;
    viewscale /= 3950.;

	hdc = GetDC(hwndMain);

	GetClientRect(hwndMain, &rect);
	SetMapMode (hdc, MM_ISOTROPIC);
//	SetWindowExtEx (hdc, 5, 5, NULL);
	SetWindowExtEx (hdc, 180, 90, NULL);
	SetViewportExtEx (hdc, rect.right / 2, -rect.bottom / 2, NULL);
	SetViewportOrgEx (hdc, rect.right / 2,  rect.bottom / 2, NULL);

	vp_filep (stdout);
	vp_erase ();

	strcpy (projection, "flat");

	if (projection[0] == 'r')
		vp_scale (1. / viewscale, 1. / viewscale);
    else
		vp_scale (5. / (180. * CONV), 5. / (180. * CONV));

	clipdist = my_acos (1. / (1. + viewdist));

	vp_color(0x00000000);
	vp_penup ();
	vp_fat (1);

	if (projection[0] == 'r') {
		horizon = tan (clipdist / 2.);
		for (j = 0; j <= 360 * 10; j++) {
			vp_upendn (horizon * cos (j * CONV / 10.),
					horizon * sin (j * CONV / 10.));
		}
	} else {
/*
		vp_move (-5., -2.5);
		vp_draw (5., -2.5);
		vp_draw (5., 2.5);
		vp_draw (-5., 2.5);
		vp_draw (-5., -2.5);
*/
		vp_move (-100., -50);
		vp_draw (100., -50);
		vp_draw (100., 50);
		vp_draw (-100., 50);
		vp_draw (-100., -50);
	}

	strcpy(filename[0], "d:/graphics/gis/wdb-ii/europe/cil.cbd");
	strcpy(filename[1], "");

	for (filenum = 0; strlen(filename[filenum]); filenum++) {
		Ifile = open (filename[filenum], O_RDONLY | O_BINARY);
		if (Ifile == -1)
			continue;

		for (j = 0; j < 20; j++)
			colfile[j] = 0x00000000;
			
		sptr = strchr (filename[filenum], '.');
		if (sptr == NULL)
			continue;
		if (strcmp (sptr, ".cbd") != 0)
			continue;

		strcpy (string2, sptr - 3);
		string2[3] = '\0';

		colfile[1] = 0x00ff0000;
		colfile[2] = 0x00ff0000;

/*
 * Read in the file header, check it for the correct magic number,
 * and learn the address of the segment dictionary
 */
		read (Ifile, &header, sizeof header);

		if (header.magic != CBD_MAGIC) {
			MBPrintf ("", "File has bad magic number %X != %X\n",
					header.magic, CBD_MAGIC);
			goto error;
		}

/* allocate space for the segment buffer */
		segbufsize = 2 * header.segmax;
		segbuf = (BIT16 *) malloc (50 + segbufsize);

/* allocate space for the segment dictionary */
		sdbuf = (struct segdict *) malloc (100 + header.segsize);
		sd = sdbuf;
		sd++;

/* Go read in the segment dictionary (it's at the end of the file) */
		lseek (Ifile, header.dictaddr, SEEK_SET);
		j = read (Ifile, sd, header.segsize);

		if (j < header.segsize) {
			MBPrintf("", "File segment dictionary expecting %d got %d\n",
					header.segsize, j);
			close (Ifile);
			free (sdbuf);
			free (segbuf);
			goto error;
		}

/*
 * Now look at each segment and decide if we're
 * going to keep it or not
 */
		segcount = header.segcount;

		sd = sdbuf;
		for (iseg = 1; iseg <= segcount; iseg++) {
			sd++;		/* does this really work? wow! */

			if (colfile[sd->rank] == 0)
				continue;

			if (projection[0] == 'r') {
				longi = sd->maxlong * CONV / 3600.;
				lati = sd->maxlat * CONV / 3600;

				x = sin (longi - viewlong) * cos (lati);
				y = cos (longi - viewlong) * cos (lati);
				z = sin (lati);

				y1 = y * cos (viewlat) + z * sin (viewlat);
				z1 = -y * sin (viewlat) + z * cos (viewlat);
				x1 = x;

				longi = sd->minlong * CONV / 3600;
				lati = sd->minlat * CONV / 3600;
				x = sin (longi - viewlong) * cos (lati);
				y = cos (longi - viewlong) * cos (lati);
				z = sin (lati);

				y2 = y * cos (viewlat) + z * sin (viewlat);
				z2 = -y * sin (viewlat) + z * cos (viewlat);
				x2 = x;

				xc = (x1 + x2) / 2.;
				yc = (y1 + y2) / 2.;
				zc = (z1 + z2) / 2.;
				x = sqrt (xc * xc + yc * yc + zc * zc);
				xc /= x;
				yc /= x;
				zc /= x;

				dist = my_acos (xc * x1 + yc * y1 + zc * z1);
				cent = my_acos (yc);

				if (cent - dist > clipdist)
					continue;
			}

			vp_color (colfile[sd->rank]);
			if (sd->rank == 1)
				vp_fat (1);
			else
				vp_fat (0);

			lseek (Ifile, sd->absaddr, SEEK_SET);
			read (Ifile, &sb, sizeof sb);

			if (sd->nbytes > segbufsize) {
				MBPrintf("", "Segment %d needs %d bytes; buffer limit is %d.\n", iseg, sd->nbytes, segbufsize);
				close (Ifile);
				free (sdbuf);
				free (segbuf);
				goto error;
			}
			read (Ifile, segbuf, sd->nbytes);

			k = 0;

			oln = sb.orgx;
			olt = sb.orgy;

#ifdef DEBUG
			count = 1;
			deg = abs (oln) / 3600;
			min = (abs (oln) - deg * 3600) / 60;
			sec = abs (oln) - deg * 3600 - min * 60;
			if (oln > 0)
				fprintf (stderr, "+ %d %d %d   %d\n", deg, min, sec, count);
			else
				fprintf (stderr, "- %d %d %d   %d\n", deg, min, sec, count);
#endif

			vp_penup ();
//			project (CONV * olt / 3600., CONV * oln / 3600.);
			project (olt / 3600., oln / 3600.);

			for (jstroke = 1; jstroke <= sb.nstrokes; jstroke++) {
				if (segbuf[k] & SHORTFLAG) {
#ifdef DEBUG
					fprintf (stderr, "s ");
#endif
/* Flag bit on: unpack a 16-bit field into dx and dy */
					i32 = segbuf[k++];
					if (i32 > 0)
					i32 &= ~SHORTFLAG;
					idy = i32 & 0xFF;
					if (idy & 0x80)
						idy |= ~0xFF;	/* extend sign */
					idx = i32 >> 8;
					if (idx & 0x80)
						idx |= ~0xBF;	/* extend sign */
				} else {
#ifdef DEBUG
					fprintf (stderr, "L ");
#endif
/* Flag bit off: take dx and dy from 32-bit fields. */
					idx = segbuf[k++];
					if (idx < 0)
						idx |= SHORTFLAG;
/*
 * You said you wanted the fix for the problems of Wrangell island, et. al.  I
 * fixed the problem, and I'm including the fix...
 *  
 * The problem is in domap.c, not the database.  The following change should be
 * made:
 *  
 *                         idx |= SHORTFLAG;
 * >                   idx = (idx << 16) | segbuf[k];
 *                     k++;
 *  
 *                         idx |= SHORTFLAG;
 *                     twiddle16 (&segbuf[k]);
 * >                   idx = (idx << 16) | (unsigned short)segbuf[k];
 *                     k++;
 *  
 * The sign bit of the lower 16 bits was extended, overwriting the upper 16
 * bits from the previous word.  This caused a very large positive number (to
 * go from -180 West hemisphere to 180 East hemisphere) to become a smaller
 * negative number, causing the longitude to go below -180 degrees.
 *  
 * Hope you are still interested!
 *  
 * Bill Marsh
 * bmarsh@nosc.mil
 * 
 */
					idx = (idx << 16) | (unsigned short) segbuf[k];
					k++;

					idy = segbuf[k];
					k++;
					if (idy < 0)
						idy |= SHORTFLAG;
					idy = (idy << 16) | segbuf[k];
					k++;
				}
				olt = olt + idy;
				oln = oln + idx;

#ifdef DEBUG
				count++;
				deg = abs (oln) / 3600;
				min = (abs (oln) - deg * 3600) / 60;
				sec = abs (oln) - deg * 3600 - min * 60;
				if (oln > 0)
					fprintf (stderr, "  + %d %d %d   %d\n", deg, min, sec, count);
				else
					fprintf (stderr, "  - %d %d %d   %d\n", deg, min, sec, count);
#endif

//				project (CONV * olt / 3600., CONV * oln / 3600.);
				project (olt / 3600., oln / 3600.);
			}
		}
		close (Ifile);
		free (sdbuf);
		free (segbuf);
    }

	ReleaseDC(hwndMain, hdc);

	while (GetMessage (&msg, NULL, 0, 0)) {
		TranslateMessage(&msg); 
		DispatchMessage(&msg); 
	}

	return msg.wParam;

error:
	ReleaseDC(hwndMain, hdc);
	return 0;
}

.