#include <stdio.h>
#include <math.h>
#include <string.h>
#include <strings.h>
#include <fcntl.h>
#define PRO 64

unsigned char c1, c2, c3, c4, *argbbuf;
unsigned char r, g, b, a;
unsigned int i1, i2, i3, i4;
float readintelfloat(), *xbuf, *ybuf, *zbuf, *backbuf, getrms();

main (ac, av)
int ac;
char **av;
	{
	int nx, ny, i, index, fd, bgfd, namel;
	float val, max, rms;
	char root[120], name[120], backfile[120];
	if (ac < 2)
		{
		fprintf(stderr, "3compimage rootname [background.flt] >outfile.argb\n");
		exit(0);
		}
	/* Open and binary float (for gray) background file, if named */
	bgfd = -1;
	if (ac > 2)
		{
		sprintf(name, "%s", av[2]);
		bgfd = open(name, O_RDONLY, 0664);
		if (bgfd < 0)
			{
			fprintf(stderr,
			"3compimage: can't open background image %s, abort\n",
				name);
			exit(-1);
			}
		}
		
	/* strip .x, .y, or .z from root name, if present */
	namel = strlen(av[1]);
	strcpy(root, av[1]);
	if (namel > 2 && root[namel-2] == '.')
		if (root[namel-1] == 'x' || root[namel-1] == 'y' || root[namel-1] == 'z')
			root[namel-2] = '\0';
	max = 0;
	/* Read in x-component file */
	sprintf(name, "%s.x", root);
	fd = open(name, O_RDONLY, 0664);
	if (fd < 0)
		{
		fprintf(stderr, "3compimage: can't open %s, abort\n", name);
		exit(-1);
		}
	/* read two integers from front of file, nx and ny */
	read(fd, &nx, sizeof(int));
	read(fd, &ny, sizeof(int));
	if (nx < 1 || ny < 1 || nx > 1000000 || ny > 1000000)
		{
		fprintf(stderr, "3compimage: nx=%d ny=%d improper, abort\n",
			nx, ny);
		exit(-1);
		}
	xbuf = (float *)malloc(nx*ny*sizeof(float));
	read(fd, xbuf, nx*ny*sizeof(float));

	/* Read background float image (for grayscale), if named */
	if (bgfd >= 0)
		{
		backbuf = (float *)malloc(nx*ny*sizeof(float));
		read(bgfd, backbuf, nx*ny*sizeof(float));
		close(bgfd);
		}

	/* Scale x-component to red intensity */
	for (i=0; i < nx*ny; i++)
		{
		val = xbuf[i];
		if (val < -1e30 || val > 1e30)
			{
			fprintf(stderr, "3compimage: i=%d val=%g improper, abort\n",
				i, val);
			exit(-1);
			}
		val = (val<0 ? -1 : 1)*
			(float)(pow((double)(val<0 ? -val : val), (double)(0.8)));
		if ((val<0 ? -val : val) > max)
			max = val<0 ? -val : val;
		xbuf[i] = val;
		}
	close(fd);
	/* Done reading x-component */

	/* Read in y-component */
	sprintf(name, "%s.y", root);
	fd = open(name, O_RDONLY, 0664);
	if (fd < 0)
		{
		fprintf(stderr, "3compimage: can't open %s, abort\n", name);
		exit(-1);
		}
	/* read two integers from front of file, nx and ny */
	read(fd, &index, sizeof(int));
	if (index != nx)
		{
		fprintf(stderr, "3compimage: nx != %d in %s, found %d, abort\n",
			nx, name, index);
		exit(-1);
		}
	read(fd, &index, sizeof(int));
	if (index != ny)
		{
		fprintf(stderr, "3compimage: ny != %d in %s, found %d, abort\n",
			ny, name, index);
		exit(-1);
		}
	ybuf = (float *)malloc(nx*ny*sizeof(float));
	read(fd, ybuf, nx*ny*sizeof(float));
	/* Scale y-component to red intensity */
	for (i=0; i < nx*ny; i++)
		{
		val = ybuf[i];
		if (val < -1e30 || val > 1e30)
			{
			fprintf(stderr, "3compimage: i=%d val=%g improper, abort\n",
				i, val);
			exit(-1);
			}
		val = (val<0 ? -1 : 1)*
			(float)(pow((double)(val<0 ? -val : val), (double)(0.8)));
		if ((val<0 ? -val : val) > max)
			max = val<0 ? -val : val;
		ybuf[i] = val;
		}
	close(fd);
	/* Done reading y-component */
	/* Read in y-component */
	sprintf(name, "%s.z", root);
	fd = open(name, O_RDONLY, 0664);
	if (fd < 0)
		{
		fprintf(stderr, "3compimage: can't open %s, abort\n", name);
		exit(-1);
		}
	/* read two integers from front of file, nx and ny */
	read(fd, &index, sizeof(int));
	if (index != nx)
		{
		fprintf(stderr, "3compimage: nx != %d in %s, found %d, abort\n",
			nx, name, index);
		exit(-1);
		}
	read(fd, &index, sizeof(int));
	if (index != ny)
		{
		fprintf(stderr, "3compimage: ny != %d in %s, found %d, abort\n",
			ny, name, index);
		exit(-1);
		}
	zbuf = (float *)malloc(nx*ny*sizeof(float));
	read(fd, zbuf, nx*ny*sizeof(float));
	/* Scale x-component to red intensity */
	for (i=0; i < nx*ny; i++)
		{
		val = zbuf[i];
		if (val < -1e30 || val > 1e30)
			{
			fprintf(stderr, "3compimage: i=%d val=%g improper, abort\n",
				i, val);
			exit(-1);
			}
		val = (val<0 ? -1 : 1)*
			(float)(pow((double)(val<0 ? -val : val), (double)(0.8)));
		if ((val<0 ? -val : val) > max)
			max = val<0 ? -val : val;
		zbuf[i] = val;
		}
	close(fd);
	/* Done reading z-component */

	/* Create scaled 3-color intensity from x,y,z components */
	rms = getrms(nx*ny);
	if (bgfd >= 0)
		scalebg(nx*ny, av[2]);
	argbbuf = (unsigned char *)malloc(nx*ny*4);
	for (i=0; i < nx*ny; i++)
		{
		a = (unsigned char)(255);
		val = xbuf[i];
		val = val<0 ? -val : val;
		index = (int)(val/3/rms*255);
		if (index < 0)
			index = 0;
		else if (index > 255)
			index = 255;
		r = (unsigned char)(index);
		val = ybuf[i];
		val = val<0 ? -val : val;
		index = (int)(val/3/rms*255);
		if (index < 0)
			index = 0;
		else if (index > 255)
			index = 255;
		g = (unsigned char)(index);
		val = zbuf[i];
		val = val<0 ? -val : val;
		index = (int)(val/3/rms*255);
		if (index < 0)
			index = 0;
		else if (index > 255)
			index = 255;
		b = (unsigned char)(index);
		/* Combine background image */
		if (bgfd >= 0)
			{
			val = (r + g + b)/3;
			if (val < PRO)
				{
				index = (PRO - val)/PRO*backbuf[i] + r;
				if (index < 0)
					index = 0;
				else if (index > 255)
					index = 255;
				r = (unsigned char)(index);
				index = (PRO - val)/PRO*backbuf[i] + g;
				if (index < 0)
					index = 0;
				else if (index > 255)
					index = 255;
				g = (unsigned char)(index);
				index = (PRO - val)/PRO*backbuf[i] + b;
				if (index < 0)
					index = 0;
				else if (index > 255)
					index = 255;
				b = (unsigned char)(index);
				}
			}
		argbbuf[4*i] = a;
		argbbuf[4*i + 1] = r;
		argbbuf[4*i + 2] = g;
		argbbuf[4*i + 3] = b;
		}
	write(1, argbbuf, nx*ny*4);
	}

int readintelint()
	{
	read(0, &c4, 1);
	read(0, &c3, 1);
	read(0, &c2, 1);
	read(0, &c1, 1);
	i1 = c1;
	i2 = c2;
	i3 = c3;
	i4 = c4;
	/*fprintf(stderr, "%d %d %d %d\n", i1<<24, i2<<16, i3<<8, i4);*/
	return((i1<<24) + (i2<<16) + (i3<<8) + i4);
	/*return((i1*16777216) + (i2*65536) + (i3*256) + i4);*/
	}

float readintelfloat()
	{
	unsigned ibuf;
	float fval, *fv;
	read(0, &c4, 1);
	read(0, &c3, 1);
	read(0, &c2, 1);
	read(0, &c1, 1);
	i1 = c1;
	i2 = c2;
	i3 = c3;
	i4 = c4;
	ibuf = i1<<24 | i2<<16 | i3<<8 | i4;
	fv = (float *)(&ibuf);
	fval = *fv;
	return(fval);
	}

writeargb()
	{
	write(1, &a, 1);
	write(1, &r, 1);
	write(1, &g, 1);
	write(1, &b, 1);
	}

float getrms(n)
int n;
	{
	int j;
	float sum;
	sum = 0;
	for (j=0; j<n; j++)
		sum += xbuf[j] * xbuf[j];
	for (j=0; j<n; j++)
		sum += ybuf[j] * ybuf[j];
	for (j=0; j<n; j++)
		sum += zbuf[j] * zbuf[j];
	sum /= 3*n;
	return((float)(sqrt((double)(sum))));
	}
	
scalebg(n, name)
int n;
char *name;
	{
	int j;
	float min, max, val;
	min = backbuf[0];
	max = backbuf[0];
	for (j=0; j<n; j++)
		{
		val = backbuf[j];
		if (val < -1e30 || val > 1e30)
			{
			fprintf(stderr, "3compimage: i=%d val=%g improper in background %s, abort\n",
				j, val, name);
			exit(-1);
			}
		if (val < min)
			min = val;
		if (val > max)
			max = val;
		}
	if (min == max)
		{
		fprintf(stderr,
		"3compimage: min=max=%g improper in background %s, abort\n",
			min, name);
		exit(-1);
		}
	for (j=0; j<n; j++)
		{
		backbuf[j] = (backbuf[j] - min)/(max - min)*255;
		}
	}