/* img2min_rgb4min.c
 *
 * read the 2 minute img file world_grav.img.7.2 and create
 * a 4 minute file of r,g,b values.  Set the color to 0,0,0
 * on land areas (defined as 3 or 4 of the four pixels averaged
 * are land).  Compute the average gravity and convert it to
 * hue using a Laplace cdf with zero mean and sigma = 26.5 mGal.
 * Compute the average gradients N to S and W to E, and choose
 * whichever is large in magnitude.  Convert this to a shade
 * factor using a Laplace cdf with zero mean and sigma = 23.3
 * (mGal * 10)/pixel.  Write a raw rgb stream out.

****MODIFIED 15 May 1996 to change rgb to 128,128,128 over land,
for application by Brian Kelly for  NOAA flyer.
 *
 *
 * 30 January 1996
 */

#include "/amos1/gmt/src/gmt.h"

#define NXIN 10800		/* 360 x 30  */
#define NTOTIN 21600		/* NXIN * 2  */
#define NYOUT 3168
#define NXOUT 5400		/* 360 x 15  */
#define NBOXES 300

#define MASKSIZE 8553600

static unsigned char maskset[8] = {128, 64, 32, 16, 8, 4, 2, 1};

unsigned char landmask[MASKSIZE];

short int data[NTOTIN];

main(argc, argv)
int     argc;
char    **argv;
{
	int	i,j, ii, ij, k, east, south, deriv, rr, gg, bb, nland;
	short	a, b, c, d;
	FILE	*fp;
	double	grav_factor = 0.025 * M_SQRT2 / 26.5;
	double	grad_factor = M_SQRT2 / 23.3;
	double	xgrav, intensity, h, s, v;
	unsigned char	rgb[3]; 

	if (argc != 3) {
		fprintf(stderr,"usage:  img2min_rgb4min.c <imgfile> <bitmaskfile>\n");
		exit(-1);
	}

	/* Read bitmask:  */

	if ( (fp = fopen(argv[2], "r") ) == NULL) {
		fprintf(stderr,"ERROR opening file %s.\n", argv[2]);
		exit(-1);
	}
        if ( (fread(landmask, MASKSIZE, 1, fp)) != 1) {
                fprintf(stderr,"img2min_rgb4min:  Error reading maskfile.\n");
                exit(-1);
        }
        fclose (fp);

	/* Now get started:  */

	if ( (fp = fopen(argv[1], "r") ) == NULL) {
		fprintf(stderr,"ERROR opening file %s.\n", argv[1]);
		exit(-1);
	}
	for (j = 0, ij = 0; j < NYOUT; j++, ij += NTOTIN) {
		if ( (fread((char *)data, 2, NTOTIN, fp)) != NTOTIN) {
			fprintf(stderr,"ERROR reading input.\n");
			exit(-1);
		}
		/* Check data and demux the tracks:  */
		for (i = 0; i < NTOTIN; i++) {
			if (data[i] == 32767) {
				fprintf(stderr,"ERROR values have 32767.\n");
				exit(-1);
			}
			if (data[i]%2) data[i]--;
		}
		for (i = 0; i < NXOUT; i++, k++) {
			nland = 0;
			ii = 2*i;

			a = data[ii + 1];
			k = ij + ii + 1;
			if (landmask[k/8] & maskset [k%8]) nland++;

			b = data[ii];
			k = ij + ii;
			if (landmask[k/8] & maskset [k%8]) nland++;

			c = data[ii+NXIN];
			k = ij + ii + NXIN;
			if (landmask[k/8] & maskset [k%8]) nland++;

			d = data[ii+NXIN+1];
			k = ij + ii + NXIN + 1;
			if (landmask[k/8] & maskset [k%8]) nland++;

			if (nland < 3) {

				east  = ((a+d)-(b+c))/2;
				south = ((c+d)-(a+b))/2;
				if (abs(east)>abs(south)) {
					deriv = east;
				}
				else {
					deriv = south;
				}
				/* scale derivative into Laplace CDF in range 0.6 * [-1, 1].  */
				if (deriv < 0) {
					intensity = 0.6 * (exp(deriv*grad_factor) - 1.0);
				}
				else {
					intensity = 0.6 * (1.0 - exp(-deriv*grad_factor));
				}

				/* get average gravity and scale into hue in range 260 * [0, 1]:  */
				xgrav = (a + b + c + d) * grav_factor;
				if (xgrav < 0.0) {
					h = 260.0 * (1.0 - 0.5 * exp(xgrav));
				}
				else {
					h = 130.0 * exp(-xgrav);
				}
				s = 0.6;
				v = 1.0;
				if (intensity > 0) {
					if (s != 0.0) s = (1.0 - intensity) * s + intensity * gmtdefs.hsv_max_saturation;
					v = (1.0 - intensity) * v + intensity * gmtdefs.hsv_max_value;
				}
				else {
					if (s != 0.0) s = (1.0 + intensity) * s - intensity * gmtdefs.hsv_min_saturation;
					v = (1.0 + intensity) * v - intensity * gmtdefs.hsv_min_value;
				}
				if (v < 0.0) v = 0.0;
				if (s < 0.0) s = 0.0;
				if (v > 1.0) v = 1.0;
				if (s > 1.0) s = 1.0;
				gmt_hsv_to_rgb (&rr, &gg, &bb, h, s, v);
				rgb[0] = rr;
				rgb[1] = gg;
				rgb[2] = bb;
			}
			else {
				rgb[0] = rgb[1] = rgb[2] = 128;
			}

			fwrite(rgb, 1, 3, stdout);
		}
	}

	fclose (fp);

	exit (0);
}