/* img_drape.c <topo_lowpass.img> -O<grav_drape.img> -M<min> -N<swath> [-V]
 *
 * Expect to find grav0.img, grav1.img, ..., grav6.img in current
 * working directory.  Read these and topo_lowpass.img into memory,
 * using swaths of -N<swath> rows.  All files are -M<min> files.
 * Then create gravity draped over lowpass topo by interpolating
 * among down?.img files.  Write result to <grav_drape.img> file.
 *
 
23 February 2009 adjusted to use check_imgfile_size and eliminate
switches for pixel size in minutes and swath size.  -- WHFS

 */
 
#include "img_predict.h"

main (argc, argv)
int	argc;
char	**argv;

{
	double	z, z1, t;
	short int	*s;
	int	verbose = 0, error = 0;
	int	i, j, k, nximg, nyimg, nystrip, iswath, nswaths, n, n6;
	size_t	kbyte;
	char	*infile, *outfile, downfile[16];
	FILE	*fpin[8], *fpout;
	
	infile = (char *)NULL;
	outfile = (char *)NULL;
	

	/* Parse arg list:  */
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
				case 'O':
					outfile = &argv[i][2];
					break;
				case 'V':
					verbose = 1;
					break;
				default:
					fprintf (stderr, "img_drape:  Unsupported option:  %s\n", argv[i]);
					error ++;
					break;
			}
		}
		else {
			infile = argv[i];
		}
	}
	if (infile == (char *)NULL) {
		fprintf (stderr, "img_drape:  You must supply an input file.\n");
		error ++;
	}
	if (outfile == (char *)NULL) {
		fprintf (stderr, "img_drape:  You must supply -O<output file>.\n");
		error ++;
	}
	
	if (error) {
		fprintf (stderr, "usage:  img_drape <infile> -O<outfile> [-V]\n\n");
		fprintf (stderr, "\t<infile> and <outfile> are img files.\n");
		fprintf (stderr, "\tOption -V reports verbose operations to stderr.\n");
		exit (-1);
	}


	if (check_imgfile_size (infile, &nximg, &nyimg) ) {
		fprintf (stderr, "img_drape:  FATAL ERROR.  Cannot make sense of size of %s\n", infile);
		exit (EXIT_FAILURE);
	}

	/* Set up some stuff  */	
	nystrip = 576;
	
	if ( nyimg%nystrip != 0) {
		fprintf (stderr, "img_drape:  nystrip = %d  does not fit nyimg = %d\n",
			nystrip, nyimg);
		exit (-1);
	}
	nswaths = nyimg / nystrip;
	n = nystrip * nximg;	/* n of pixels in a swath  */
	kbyte = n * 16;		/* 8 groups of short ints  */
	n6 = n * 6;


	/* Malloc stuff:  */
	if ( (s = (short int *)malloc(kbyte) ) == NULL) {
		fprintf (stderr, "img_drape:  malloc failure for s array.\n");
		exit (-1);
	}
	if (verbose) printf("img_drape.c  files are %d by %d; nystrip is %d; %6.1lf Mbytes used.\n",
		nximg, nyimg, nystrip, kbyte/1048576.0);


	/* Open files:  */
	for (i = 0; i <=6; i++) {
		sprintf(downfile, "down%1.1d.img", i);
		if ( (fpin[i] = fopen(downfile, "r")) == NULL) {
			fprintf (stderr, "img_drape:  Error.  Cannot open r %s\n", infile);
			exit (-1);
		}
	}
	if ( (fpin[7] = fopen(infile, "r")) == NULL) {
		fprintf (stderr, "img_drape:  Error.  Cannot open r %s\n", infile);
		exit (-1);
	}
	if ( (fpout = fopen(outfile, "w")) == NULL) {
		fprintf (stderr, "img_drape:  Error.  Cannot open w %s\n", outfile);
		exit (-1);
	}


	/* Do swaths  */
	for (iswath = 0; iswath < nswaths; iswath++) {
		if (verbose) printf ("Strip %d of %d.\n", iswath, nswaths);
		for (i = 0; i < 8; i++) {
			k = i * n;
			if ( (fread( (char *)&s[k], (size_t)2, (size_t)n, fpin[i]) ) != n) {
				fprintf (stderr, "img_drape:  Error reading input %d\n", i);
				exit (-1);
			}
		}
		/* Overwrite first group with new value:  */
		for (i = 0; i < n; i++) {
			z = 0.001 * s[k + i];
			if (z >= 0.0) continue;
			if (z <= -6.0) {
				s[i] = s[n6 + i];
				continue;
			}
			z = -z;
			z1 = floor(z);
			t = z - z1;
			j = n * (int)z1;
			s[i] = (short int) rint(t*s[j+n+i] + (1.0 - t)*s[j+i]);
		}
		
		/* Write out swath  */
		if ( (fwrite( (char *)s, (size_t)2, (size_t)n, fpout) ) != n) {
			fprintf (stderr, "img_drape:  Error writing output\n");
			exit (-1);
		}
	}
	/* close files and free mallocs  */
	exit (0);
}