#include <stdio.h>
#include <math.h>
#include <string.h>
#include <strings.h>
#include <fcntl.h>

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

main (ac, av)
int ac;
char **av;
	{
	int nx, ny, i, index, fd, namel;
	float val, max, rms;
	char root[120], name[120];
	if (ac < 2)
		{
		fprintf(stderr, "3compimage rootname >outfile.argb\n");
		exit(0);
		}
	/* 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 */
	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 */
	nx = readintelint(fd);
	ny = readintelint(fd);
	if (nx < 1 || ny < 1 || nx > 1000000 || ny > 1000000)
		{
		fprintf(stderr, "3compimage: nx=%d ny=%d improper, abort\n",
			nx, ny);
		exit(-1);
		}
	/* read five more integers from front (20 extra bytes); the last is a one */
	for (i=0; i<5; i++)
		index = readintelint(fd);
	if (index != 1)
		fprintf(stderr,
		"convimage: warning, unexpected seventh integer=%d\n", index);
	fprintf(stderr, "convimage: nx=%d ny=%d\n", nx, ny);
	xbuf = (float *)malloc(nx*ny*sizeof(float));
	for (i=0; i < nx*ny; i++)
		{
		val = readintelfloat(fd);
		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 */
	index = readintelint(fd);
	if (index != nx)
		{
		fprintf(stderr, "3compimage: nx != %d in %s, found %d, abort\n",
			nx, name, index);
		exit(-1);
		}
	index = readintelint(fd);
	if (index != ny)
		{
		fprintf(stderr, "3compimage: ny != %d in %s, found %d, abort\n",
			ny, name, index);
		exit(-1);
		}
	/* read five more integers from front (20 extra bytes); the last is a one */
	for (i=0; i<5; i++)
		index = readintelint(fd);
	if (index != 1)
		fprintf(stderr,
		"convimage: warning, unexpected seventh integer=%d\n", index);
	ybuf = (float *)malloc(nx*ny*sizeof(float));
	for (i=0; i < nx*ny; i++)
		{
		val = readintelfloat(fd);
		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 */
	index = readintelint(fd);
	if (index != nx)
		{
		fprintf(stderr, "3compimage: nx != %d in %s, found %d, abort\n",
			nx, name, index);
		exit(-1);
		}
	index = readintelint(fd);
	if (index != ny)
		{
		fprintf(stderr, "3compimage: ny != %d in %s, found %d, abort\n",
			ny, name, index);
		exit(-1);
		}
	/* read five more integers from front (20 extra bytes); the last is a one */
	for (i=0; i<5; i++)
		index = readintelint(fd);
	if (index != 1)
		fprintf(stderr,
		"convimage: warning, unexpected seventh integer=%d\n", index);
	zbuf = (float *)malloc(nx*ny*sizeof(float));
	for (i=0; i < nx*ny; i++)
		{
		val = readintelfloat(fd);
		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 */

	rms = getrms(nx*ny);
	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);
		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(f)
int f;
	{
	read(f, &c4, 1);
	read(f, &c3, 1);
	read(f, &c2, 1);
	read(f, &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(f)
int f;
	{
	unsigned ibuf;
	float fval, *fv;
	read(f, &c4, 1);
	read(f, &c3, 1);
	read(f, &c2, 1);
	read(f, &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))));
	}