/* img_erf.c <highpass_grav.img> -O<adjusted_grav.img> -S<sigma> [-K<kill>] [-P<pow>] [-T<scale>] [-V]
 *
 * Read an imgfile into memory.  Adjust values x[i] by a factor
 * erf(fabs(x[i]) / (sigma * M_SQRT2)).  Optionally [-K], if
 * fabs(x[i]) < kill, set x[i] = 0.0.  Optionally [-P<pow>]
 * adjust the factor by pow(factor, pow).
 * 
 * Idea is to adjust residual gravity so as to suppress noise while
 * passing real signals.
 *
 * Option -T<scale> will rescale resulting gravity by a factor (e.g.
 * to convert it to a topography perturbation for topo inversion).
 * Use the right scale for img gravity (mGal * 10) to img topo (m),
 * e.g. 1.454 for density 2.67 in seawater.
 
23 Feb 2009, WHFS:  Modified to use check_imgfile_size
 
 *
 */
 
#include "img_predict.h"


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

{
	double	sigma = 0.0, xkill = 0.0, p = 1.0, t = 1.0;
	double	erf_factor, x, y, z;
	short int	*s;
	int	verbose = 0, error = 0;
	int	i, nximg, nyimg, n;
	size_t	kbyte;
	char	*infile, *outfile;
	FILE	*fp;
	
	infile = (char *)NULL;
	outfile = (char *)NULL;
	

	/* Parse arg list:  */
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
				case 'K':
					xkill = atof(&argv[i][2]);
					break;
				case 'O':
					outfile = &argv[i][2];
					break;
				case 'P':
					p = atof(&argv[i][2]);
					break;
				case 'S':
					sigma = atof(&argv[i][2]);
					break;
				case 'T':
					t = atof(&argv[i][2]);
					break;
				case 'V':
					verbose = 1;
					break;
				default:
					fprintf (stderr, "img_erf:  Unsupported option:  %s\n", argv[i]);
					error ++;
					break;
			}
		}
		else {
			infile = argv[i];
		}
	}
	if (p <= 0.0) {
		fprintf (stderr, "img_erf:  Option -P<p> must be a positive value.\n");
		error ++;
	}
	if (sigma <= 0.0) {
		fprintf (stderr, "img_erf:  You must supply -S<sigma> in input data units (mGal * 10 ?)\n");
		error ++;
	}
	if (xkill < 0.0) {
		fprintf (stderr, "img_erf:  Option -K<xkill> must be positive, in input data units (mGal * 10 ?)\n");
		error ++;
	}
	if (infile == (char *)NULL) {
		fprintf (stderr, "img_erf:  You must supply an input file.\n");
		error ++;
	}
	if (outfile == (char *)NULL) {
		fprintf (stderr, "img_erf:  You must supply -O<output file>.\n");
		error ++;
	}
	
	if (error) {
		fprintf (stderr, "usage:  img_erf <infile> -O<outfile> -S<sigma> [-K<xkill>] [-P<p>] [-T<scale>] [-V]\n\n");
		fprintf (stderr, "\t<infile> and <outfile> are img files.\n");
		fprintf (stderr, "\t<sigma> is in input data units (mGal * 10, maybe).\n");
		fprintf (stderr, "\tOption -K<xkill> sets values within xkill of 0 to 0 (use data units).\n");
		fprintf (stderr, "\tOption -P<p> raises erf factor by pow(factor, p).\n");
		fprintf (stderr, "\tOption -T<scale> scales output.  E.g. use 1.454 for g to t.\n");
		fprintf (stderr, "\tOption -V reports verbose operations to stderr.\n");
		exit (-1);
	}

	if (check_imgfile_size (infile, &nximg, &nyimg) ) {
		fprintf (stderr, "img_erf:  FATAL ERROR.  Cannot make sense of size of %s\n", infile);
		exit (EXIT_FAILURE);
	}
	
	n = nyimg * nximg;	/* n of pixels in file  */


	/* Malloc stuff:  */
	kbyte = 2 * n;
	if ( (s = (short int *)malloc(kbyte) ) == NULL) {
		fprintf (stderr, "img_erf:  malloc failure for s array.\n");
		exit (-1);
	}
	if (verbose) printf("img_erf.c  Files are %d by %d; %6.1lf Mbytes used.\n",
		nximg, nyimg, kbyte/1048576.0);
		
	
	/* Read input  */
	if ( (fp = fopen(infile, "r")) == NULL) {
		fprintf (stderr, "img_erf:  Error.  Cannot open r %s\n", infile);
		exit (-1);
	}
	if ( (fread ( (void *)s, (size_t)2, (size_t)n, fp) ) != n) {
		fprintf (stderr, "img_erf:  Error.  Cannot read %s\n", infile);
		exit (-1);
	}
	fclose (fp);
	
	
	
	/* Here's the guts of it:  */

	erf_factor = 1.0 / (sigma * M_SQRT2);

	for (i = 0; i < n; i++) {
		x = (double)s[i];
		y = fabs(x);
		if (y <= xkill) {
			s[i] = 0;
			continue;
		}
		if (xkill > 0.0) y -= xkill;
		z = erf(y * erf_factor);
		if (p != 1.0) z = pow(z, p);
		y *= z;
		y *= t;
		y = copysign(y, x);
		s[i] = (short int) rint(y);
	}


	/* Write output  */
	if ( (fp = fopen(outfile, "w")) == NULL) {
		fprintf (stderr, "img_erf:  Error.  Cannot open w %s\n", outfile);
		exit (-1);
	}
	if ( (fwrite ( (void *)s, (size_t)2, (size_t)n, fp) ) != n) {
		fprintf (stderr, "img_erf:  Error.  Cannot write %s\n", outfile);
		exit (-1);
	}
	fclose (fp);

	
	free ( (void *)s);

	exit (0);
}