/*--------------------------------------------------------------------
 *	$Id: gmt2rgb.c,v 1.15 2004/09/28 19:24:13 pwessel Exp $
 *
 *	Copyright (c) 1991-2004 by P. Wessel and W. H. F. Smith
 *	See COPYING file for copying and redistribution conditions.
 *
 *	This program is free software; you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
 *	the Free Software Foundation; version 2 of the License.
 *
 *	This program is distributed in the hope that it will be useful,
 *	but WITHOUT ANY WARRANTY; without even the implied warranty of
 *	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *	GNU General Public License for more details.
 *
 *	Contact info: gmt.soest.hawaii.edu
 *--------------------------------------------------------------------*/
/*
 * gmt2rgb reads either (1) an 8, 24, or 32 bit Sun rasterfile and writes out the
 * red, green, and blue components in three separate grd files, or (2) a z grd
 * file and a cpt file and compute r, g, b and write these out instead.
 *
 * Author:	Paul Wessel
 * Date:	17-SEP-2001
 * Version:	4
 *
 */

#include "gmt.h"

unsigned char *loadraw (char *file, struct rasterfile *header, int byte_per_pixel, int nx, int ny);
int     get_prime_factors(int n, int *f);
void guess_width (char *file, int byte_per_pixel);

int	raw_nx = 0, raw_ny = 0;

main (int argc, char **argv)
{
	int i, k, k3, nm, error = 0, index, irgb[3], byte_per_pixel = 3;
        int entry, one_or_zero;
	BOOLEAN pixel = FALSE, do_ztorgb = FALSE, given_I = FALSE, raw = FALSE, guess = FALSE;
	float *z;
	double w, e, s, n, dx = 0.0, dy = 0.0;
	struct GRD_HEADER grd;
	struct rasterfile header;
	char rgb[3] = {'r', 'g', 'b'}, *comp[3] = {"red", "green", "blue"};
	char *file = CNULL, *format, grdfile[BUFSIZ], *cpt_file = CNULL, layer = 0, *ptr;
	unsigned char *picture;

	argc = GMT_begin (argc, argv);
	
	/* Check and interpret the command line arguments */
	
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch(argv[i][1]) {
		
				/* Common parameters */
			
				case 'R':
				case 'V':
				case '\0':
					error += GMT_get_common_args (argv[i], &w, &e, &s, &n);
					break;
				case 'C':
					cpt_file = &argv[i][2];
					do_ztorgb = TRUE;
					break;
				case 'F':
					pixel = TRUE;
					break;
				case 'G':
					format = &argv[i][2];
					break;
				case 'I':
					GMT_getinc (&argv[i][2], &dx, &dy);
					given_I = TRUE;
					break;
				case 'L':
					layer = argv[i][2];
					break;
				case 'W':
					raw = TRUE;
					guess = TRUE;
					ptr = strtok (&argv[i][2], "/");
					entry = 0;
					while (ptr) {
						if (ptr[0] != '=') {
							switch (entry) {
								case 0:
									raw_nx = atoi(ptr);
									guess = FALSE;
									break;
								case 1:
									raw_ny = atoi(ptr);
									break;
								case 2:
									byte_per_pixel = atoi(ptr);
									break;
								default:
									break;
							}
						}
						ptr = strtok (CNULL, "/");
						entry++;
					}
					break;
					
				/* Options not recognized */
						
				default:
					error = TRUE;
					fprintf (stderr, "GMT SYNTAX ERROR:  Unrecognized option -%c\n", argv[i][1]);
					break;
			}
		}
		else
			file = argv[i];
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr,"%s %s - Write r/g/b gridfiles from a gridfile, a raw, or SUN rasterfile\n\n", GMT_program, GMT_VERSION);
		fprintf (stderr,"usage: %s <rasterfile|grdfile> -G<nametemplate> [-C<cptfile>] [-Idx[/<dy>]]\n", GMT_program);
		fprintf (stderr,"\t[-L<layer> [-R<w/e/s/n>] [-W<width/height>[/<n_bytes>]] [-V]\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr,"\t<infile> can be one of three different intput files:\n");
		fprintf (stderr,"\t  (1) An 8, 24, or 32-bit Sun rasterfile.  Use -I, -R, and -F to change the\n");
		fprintf (stderr,"\t      the default values of dx = dy = 1 and region 1/ncols/1/nrows.\n");
		fprintf (stderr,"\t  (2) A regular z grdfile.  Use -C to provide a cpt file with which\n");
		fprintf (stderr,"\t      to convert z to r/g/b triplets. -R, -I, and -F are ignored.\n");
		fprintf (stderr,"\t  (3) A RGB or RGBA raw rasterfile. Since raw rastefiles have no header, you have to\n");
		fprintf (stderr,"\t      give the image dimensions via -W\n");
		fprintf (stderr,"\t      However, you may take the chance of letting the program try to\n");
		fprintf (stderr,"\t      guess the image dimensions.\n");
		fprintf (stderr,"\t-G Give outputfile name template for the three red, green, blue gridfiles.\n");
		fprintf (stderr,"\t   The template MUST contain the format code %%c which will be replaced with r, g, and b.\n");
		fprintf (stderr,"\n\tOPTIONS:\n");
		fprintf (stderr,"\t-C color palette file to convert z to rgb.  If given, we assume a z grdfile is provided,\n");
		fprintf (stderr,"\t   else we will try to read a Sun rasterfile.\n");
		fprintf (stderr, "\t-F will force pixel registration [Default is grid registration]\n");
		fprintf (stderr, "\t-I specifies grid size(s).  Append m (or c) to <dx> and/or <dy> for minutes (or seconds)\n");
		fprintf (stderr, "\t-L Only output the given layer (r, g, or b) [Default output all three]\n");
		GMT_explain_option ('R');
		GMT_explain_option ('V');
		fprintf (stderr, "\t-W sets the size of the raw raster file. By default an RGB file (which has 3 bytes/pixel)\n");
		fprintf (stderr, "\t   is assumed. For RGBA files use n_bytes = 4\n");
		fprintf (stderr, "\t   Use -W for guessing the image size of a RGB raw file, and -W=/=/4\n");
		fprintf (stderr, "\t   if the raw image is of the RGBA type. Notice that this might be a\n");
		fprintf (stderr, "\t   bit slow because the guessing algorithm makes uses of FFTs.\n");
		exit (EXIT_FAILURE);
	}

	/* Check that the options selected are mutually consistent */
	
	if (!do_ztorgb) {
		if (!file) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input raster file\n", GMT_program);
			error++;
		}
		if ((dx == 0.0 || dy == 0.0) && !raw) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -Idx/dy\n", GMT_program);
			error++;
		}
		if (byte_per_pixel != 3 && byte_per_pixel != 4) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR: byte_per_pixel must be either 3 or 4\n", GMT_program);
			error++;
		}
		if (guess) 
			guess_width (file, byte_per_pixel);

		if (raw && raw_nx <= 0) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR:  Witdth of raw raster file must be a positive integer. Not %d\n", GMT_program, raw_nx);
			error++;
		}
		if (raw && raw_ny <= 0) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR:  Height of raw raster file must be a positive integer. Not %d\n", GMT_program, raw_ny);
			error++;
		}
	}
	else {
		if (!file) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input z grd file\n", GMT_program);
			error++;
		}
		if (!cpt_file) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify cpt file\n", GMT_program);
			error++;
		}
	}
	if (!strstr (format, "%c")) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  output template must contain %%c\n", GMT_program);
		error++;
	}
	if (!strchr ("rgb", layer)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  -L layer must be one of r, g, or b\n", GMT_program);
		error++;
	}
	if (error) exit (EXIT_FAILURE);

	GMT_put_history (argc, argv);	/* Update .gmtcommands4 */
	
	GMT_grd_init (&grd, argc, argv, FALSE);

	if (do_ztorgb) {
		GMT_read_cpt (cpt_file);
		if (GMT_read_grd_info (file, &grd)) {
			fprintf (stderr, "%s: Error opening file %s\n", GMT_program, file);
			exit (EXIT_FAILURE);
		}
		nm = grd.nx * grd.ny;
		z = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
		for (i = 0; i < 3; i++) {	/* Do the r, g, and b channels */
			if (layer && layer != rgb[i]) continue;
			if (gmtdefs.verbose) fprintf (stderr, "%s: Processing the %s components\n", GMT_program, comp[i]);
			if (GMT_read_grd (file, &grd, z, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
				fprintf (stderr, "%s: Error reading file %s\n", GMT_program, file);
				exit (EXIT_FAILURE);
			}
			sprintf (grdfile, format, rgb[i]);
			sprintf (grd.remark, "Grid of %s components in the 0-255 range", comp[i]);
			for (k = 0; k < nm; k++) {
				index = GMT_get_rgb24 (z[k], irgb);
				z[k] = (float)irgb[i];
			}
			if (GMT_write_grd (grdfile, &grd, z, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
				fprintf (stderr, "%s: Error writing file %s\n", GMT_program, grdfile);
				exit (EXIT_FAILURE);
			}
		}
	}
	else {
		if (access (file, R_OK)) {
			fprintf (stderr, "%s: Cannot find/open/read file %s\n", GMT_program, file);
			exit (EXIT_FAILURE);
		}

		if (!raw)
			picture = ps_loadraster (file, &header, FALSE, FALSE, FALSE, NULL, NULL);
		else
			picture = loadraw (file, &header, byte_per_pixel, raw_nx, raw_ny);

		if (!picture) {
			fprintf (stderr, "%s: Trouble loading/converting Sun rasterfile!\n", GMT_program);
			exit (EXIT_FAILURE);
		}
		if (header.ras_depth < 8) {
			fprintf (stderr, "%s: Sun rasterfile must be at least 8 bits deep\n", GMT_program);
			exit (EXIT_FAILURE);
		}
	
		if (pixel) {
			grd.node_offset = 1;
			one_or_zero = 0;
		}
		else {
			grd.node_offset = 0;
			one_or_zero = 1;
		}
		if (!given_I) {
			if (gmtdefs.verbose) fprintf (stderr, "%s: Assign default dx = dy = 1\n", GMT_program);
			dx = dy = 1.0;
		}
		if (w == e && s == n) {	/* R not given, provide default */
			if (gmtdefs.verbose) fprintf (stderr, "%s: Assign default -R1/%d/1/%d\n", GMT_program, header.ras_width, header.ras_height);
			w = s = 1.0;
			e = header.ras_width;
			n = header.ras_height;
		}

		grd.nx = irint ((e-w)/dx) + one_or_zero;
		grd.ny = irint ((n-s)/dy) + one_or_zero;
		if (raw && !given_I) {		/* This isn't correct because it doesn't deal with -F */
			grd.nx = raw_nx;
			grd.ny = raw_ny;
			dx = (e-w)/(raw_nx - 1);
			dy = (n-s)/(raw_ny - 1);
		}
		if (header.ras_width != grd.nx) {
			fprintf (stderr, "%s: Sun rasterfile width and -R -I do not match (%d versus %d)  Need -F?\n", GMT_program, header.ras_width, grd.nx);
			exit (EXIT_FAILURE);
		}
		if (header.ras_height != grd.ny) {
			fprintf (stderr, "%s: Sun rasterfile height and -R -I do not match (%d versus %d)  Need -F?\n", GMT_program, header.ras_height, grd.ny);
			exit (EXIT_FAILURE);
		}
		grd.x_min = w;	grd.x_max = e;
		grd.y_min = s;	grd.y_max = n;
		grd.x_inc = dx;	grd.y_inc = dy;
		nm = grd.nx * grd.ny;
	
		GMT_grd_RI_verify (&grd, 1);
	
		if (gmtdefs.verbose) fprintf (stderr, "%s: nx = %d  ny = %d\n", GMT_program, grd.nx, grd.ny);
	
		z = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);

		for (i = 0; i < 3; i++) {	/* Do the r, g, and b channels */
			if (layer && layer != rgb[i]) continue;
			if (gmtdefs.verbose) fprintf (stderr, "%s: Processing the %s components\n", GMT_program, comp[i]);
			sprintf (grdfile, format, rgb[i]);
			sprintf (grd.remark, "Grid of %s components in the 0-255 range", comp[i]);
			k3 = i;
			for (k = 0; k < nm; k++) {
				if (header.ras_depth == 8)	/* Gray ramp */
					z[k] = (float)picture[k];
				else {				/* 24-bit image */
					z[k] = (float)picture[k3];
					k3 += 3;
				}
			}
			if (GMT_write_grd (grdfile, &grd, z, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
				fprintf (stderr, "%s: Error writing file %s\n", GMT_program, grdfile);
				exit (EXIT_FAILURE);
			}
		}
		GMT_free ((void *)picture);
	}
	GMT_free ((void *)z);
	
	GMT_end (argc, argv);
}


unsigned char *loadraw (char *file, struct rasterfile *header, int byte_per_pixel, int nx, int ny) {
	/* loadraw reads a raw binary grb or rgba rasterfile of depth 24, or 32 into memory */

	int j, i, nm;
	unsigned char *buffer;
	FILE *fp;

	if ((fp = GMT_fopen (file, "rb")) == NULL) {
		fprintf (stderr, "%s: Cannot open rasterfile %s!\n", GMT_program, file);
		exit (EXIT_FAILURE);
	}

	/* Lets pretend that the raw file is a sunraster file. This way the gmt2rgb code 
	   can be used with very little changes */
	header->ras_depth = 24; 
	header->ras_width = nx; 
	header->ras_height = ny; 
	nm = nx * ny * byte_per_pixel;
	header->ras_length = nm; 

	buffer = (unsigned char *) GMT_memory (VNULL, (size_t)nm, sizeof (unsigned char), GMT_program);
	if (fread ((void *)buffer, (size_t)1, (size_t)nm, fp) != (size_t)nm) {
		if (byte_per_pixel == 3)
			fprintf (stderr, "%s: Trouble reading raw 24-bit rasterfile!\n", GMT_program);
		if (byte_per_pixel == 4)
			fprintf (stderr, "%s: Trouble reading raw 32-bit rasterfile!\n", GMT_program);
		exit (EXIT_FAILURE);
	}

	if (byte_per_pixel == 4) {		/* RGBA */
		for (i = 3, j = 4; j < nm; i += 3, j += 4) {
			buffer[i] = buffer[j];
			buffer[i+1] = buffer[j+1];
			buffer[i+2] = buffer[j+2];
		}
	}

	GMT_fclose (fp);
	return (buffer);	
}

void guess_width (char *file, int byte_per_pixel) {
	unsigned char *buffer;
	float	*work, *datac, *img_pow, pow_max = -FLT_MAX, pm;
	int	narray, img_size, k = 0, j, rgb[3], inc, n_pix, i, l, even;
	FILE *fp;

	if ((fp = GMT_fopen (file, "rb")) == NULL) {
		fprintf (stderr, "%s: Cannot open rasterfile %s!\n", GMT_program, file);
		exit (EXIT_FAILURE);
	}

	fseek (fp, (long)0, SEEK_END);
	img_size = ftell (fp);
	rewind(fp);

	n_pix = img_size / byte_per_pixel;

	buffer = (unsigned char *) GMT_memory (VNULL, (size_t)img_size, sizeof (unsigned char), GMT_program);
	datac = (float *) GMT_memory (VNULL, (size_t)2*n_pix, sizeof (float), GMT_program);
	work = (float *) GMT_memory (VNULL, (size_t)2*n_pix, sizeof(float), GMT_program);
	img_pow = (float *) GMT_memory (VNULL, (size_t)n_pix/2, sizeof (float), GMT_program);
	memset ((char *)work, 0, (size_t)(2*n_pix * sizeof(float)));

	if (fread ((void *)buffer, (size_t)1, (size_t)img_size, fp) != (size_t)img_size) {
		if (byte_per_pixel == 3)
			fprintf (stderr, "%s: Trouble_ reading raw 24-bit rasterfile!\n", GMT_program);
		if (byte_per_pixel == 4)
			fprintf (stderr, "%s: Trouble_ reading raw 32-bit rasterfile!\n", GMT_program);
		exit (EXIT_FAILURE);
	}

	inc = (byte_per_pixel == 3) ? 3: 4;
	for (j = 0; j < img_size; j += inc) {
		rgb[0] = buffer[j];
		rgb[1] = buffer[j+1];
		rgb[2] = buffer[j+2];
		/* Convert rgb to gray using the YIQ transformation */
		datac[k] = (float) YIQ(rgb);
		k += 2;
	}

	narray = n_pix;
	GMT_fourt (datac, &narray, 1, -1, 1, work);

	/* Now compute the image's power spectrum */
	for (k = 0, j = 0; k < n_pix; k+= 2, j++) {
		img_pow[j] = (datac[k]*datac[k] + datac[k+1]*datac[k+1]) / n_pix; /* I*I-conj = power */
	}

	/* I'll assume that searching on one fifth of the spectrum is enough to find the
	   line frequency. */
	for (k = 5; k < n_pix/10; k++) {
		if (img_pow[k] > pow_max) {
			pow_max = img_pow[k];	 j = k+1;
		}
	}

	/* That's the way it should be but, I don't know why, the result is transposed. Instead
	   of the number of lines I get number of columns. This is very weard and smels	BUG */
	/*raw_nx = j;		raw_ny = irint((float)n_pix / raw_nx);*/

	/* So be it */
	raw_ny = j;		raw_nx = irint((float)n_pix / raw_ny);

	if (raw_nx * raw_ny != n_pix) {
		/* Let's make another attempt to find the right nx * ny combination. The idea is that we
	   	failed by a little, so we'll look arround the approximate solution adding 1 to nx and
	   	subtracting 1 to ny. Then we revert (subtract 1 to nx and add 1 to ny). Next apply the
	   	same test with an offset of 2, and so on until the offset is 10. */
		fprintf (stderr, "%s WARNING: first test based on FFT failed to guess image dimensions.\n\tI'll do now a second try\t", GMT_program);
		k = 1;		pm = 1;		l = 1;
		while (k < 41) {
			i = raw_ny + (int)copysign((float)l,pm);
			pm *= -1.;
			j = raw_nx + (int)copysign((float)l,pm);
			if (i*j == n_pix) {	/* Got a good candidate */
				raw_ny = i;	raw_nx = j;
				fprintf (stderr, "... SUCESS (W = %d, H = %d)\n", raw_nx, raw_ny);
				break;
			}
			even = (k%2 == 0) ? 1: 0;
			if (even) l++;
			k++;
		}
	}
	else
		if (gmtdefs.verbose) fprintf (stderr, "File %s has %d Lines and %d Cols\n", file, raw_ny, raw_nx);

	/* If both attempts failed */
	if (raw_nx * raw_ny != n_pix) {
		fprintf (stderr, "FAILURE while guessing image dimensions (W = %d, H = %d)\n", raw_nx, raw_ny);
		exit (EXIT_FAILURE);
	}

	GMT_fclose (fp);
	GMT_free ((void *)buffer);
	GMT_free ((void *)datac);
	GMT_free ((void *)work);
	GMT_free ((void *)img_pow);
}
