/**************************************************************************
 *                                                                        *
 * Copyright (c) 1996 Michael Richmond and Richard Treffers               *
 *                                                                        *
 *                    This software may be copied and distributed for     *
 *                    educational, research and not for profit services   *
 *                    provided that this copyright and statement are      *
 *                    included in all such copies.                        *
 *                                                                        *
 **************************************************************************/


/*	
	shows CCD image on screen in contour-map form

		CONTOUR file [box=] [low=] [diff=] [ratio=] [levels=] [zoom=]
		             [xo=] [yo=] [erase]

	contour levels specified three different ways
			1. LEVELS=a,b,d ...
			2. LOW=a     DIFF=b    which makes level=low+n*diff
			3. LOW=a     RATIO=b   which makes level=low*(ratio^n)

	HARDCOPY not ready yet

	1/30/91 - added 'floor' to "low" and "diff" expressions  - MWR
	11/17/92 - removed call to exten -rr
*/

#include <stdio.h>
#include <ctype.h>
#include <limits.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "pcvista.h"
#include "pcvimage.h"
#include "fits.h"

#define NLEVELS_MAX	100
#define DEF_RATIO   2.0			/* use ratio=DEF_RATIO for default */

#ifdef PROTO
static int crack_levels(char *, int *);
static int gen_even_levels(int, int, int *);
static int gen_ratio_levels(int, double, int *);
static int change_level(int);
static void rebin_data(int, int, int, int *);
static void dotit(int, int);
void main(int, char **);
#else
static int crack_levels();
static int gen_even_levels();
static int gen_ratio_levels();
static int change_level();
static void rebin_data();
static void dotit();
#endif

static int got_low, got_levels, got_diff, got_ratio;
static int nlevels;
static level_index;
static int zoom = 1, bin = 1;
static int x_origin = 0, y_origin = 0;
static rebin_flag, hard_flag;
static char buf[NBUF];
static int16 data1[NMAX], data2[NMAX];

void
main(argc,argv)
int argc;
char *argv[];
{
	FITS_HANDLE fh;
	int *a, *b, *temp;
	int sr = 0, sc = 0, endrow, endcol;
	int i, k, boxno, low, diff;
	int x, row, test;
	int nr, nc, nrow, ncol, dir, temp_span, binrow;
	int levels[NLEVELS_MAX];
	double ratio;
	char *gotit, fname[NBUF];

	if (argc == 1)
		error(-1, "CONTOUR file [xo= ] [yo=] [erase] [box=] [low=] \n             [diff=] [ratio] [levels=a,b,] [hard]");	

	fh = fits_open(argv[1], "r", &nrow, &ncol);
	nr=nrow;	
	nc=ncol;

	window_init();
	image_init(TV_GREYSCALE);
	
	for (i = 2; i < argc; i++) {
		if (gotit = find("box", argv[i])) {
			boxno = (int) evaluate(gotit);
			if (getbox(boxno, &sr, &sc, &nr, &nc))
				error(-1, "box not found");
			check_box(boxno, nrow, ncol);
			continue;
		}
		if (gotit = find("levels", argv[i])) {	/* generate levels explicitly */
			nlevels = crack_levels(gotit, levels);
			got_levels = 1;
			continue;
		}
		if (gotit = find("low", argv[i])) {
			low = (int) floor(evaluate(gotit) + 0.5);
			got_low = 1;
			continue;
		}
		if (gotit = find("diff", argv[i])) {
			diff = (int) floor(evaluate(gotit) + 0.5);
			got_diff = 1;
			if (diff == 0)
				error(-1, "can't have zero difference");
			continue;
		}
		if (gotit = find("ratio", argv[i])) {
			ratio = evaluate(gotit);
			got_ratio = 1;
			if (ratio <= 1.0)
				error(-1, "ratio must be greater than 1");
			continue;
		}
		if (gotit = find("xo", argv[i])) { 
			x_origin = (int) evaluate(gotit);
			continue;
		}
		if (gotit = find("yo", argv[i])) {
			y_origin = (int) evaluate(gotit);
			continue;
		}
		if (gotit = find("zoom", argv[i])) {
			double val;
			val = evaluate(gotit);
			if (val <= 0.0)
				error(-1, "zoom must be greater than zero");
			if (val > 0.5)
				zoom = (int) (val + SMALL);
			else
				bin = (int) ((1.0/val) + SMALL);
			continue;
		}
		if (keyword("erase", argv[i])) {
			(*image_clr)();
			flush_windows();
			continue;
		}
		if (keyword("hard", argv[i])) {
			hard_flag = 1;
			error(-1, "hard copy not ready yet");
			continue;
		}
	}
	if ((bin > 1) || (zoom > 1))
		rebin_flag = 1;

	if (got_ratio + got_diff + got_levels > 1)
		error(-1, "specify only one of {diff, ratio, levels} methods");
	if ((got_ratio || got_diff) && (!got_low)) {
		autospan(fh, &temp_span, &low, sr, sc, nr, nc);
		sprintf(buf, "using low=%5d\n", low);
		(*screen_puts)(SCREEN_ADD, buf);
		got_low = 1;
	}

	/* if nothing is specfied, use DEF_RATIO as the ratio,
	   starting at the the first level above zero. */

	if (!got_diff && !got_levels && !got_ratio) {
		autospan(fh, &temp_span, &low, sr, sc, nr, nc);
		ratio = DEF_RATIO;
		low *= DEF_RATIO; 
		sprintf(buf, "using low=%5d ratio=%5lf\n", low, ratio);
		(*screen_puts)(SCREEN_ADD, buf);
		got_ratio = 1;
		got_low = 1;
	}

	if (got_diff)
		nlevels = gen_even_levels(low, diff, levels);
	if (got_ratio)
		nlevels = gen_ratio_levels(low, ratio, levels);
		
	trim_window(zoom, bin, x_origin, y_origin, &nr, &nc);
	endcol = sc + nc;
	endrow = sr + nr;

	a = data1;			/* two lines are stored in alternating data buffers */
	b = data2;

	fits_get_data(fh, sr, sc, a, nc);
	if (rebin_flag)
		rebin_data(bin, zoom, nc, a);

	for (row = sr, x = x_origin; row < (endrow - 1); row += bin) {
		binrow = (row/bin)*zoom;
		fits_get_data(fh, row + 1, sc, b, nc);
		if (rebin_flag)
			rebin_data(bin, zoom, nc, b);

		for (i = (sc/bin)*zoom; i < ((endcol - 1)/bin)*zoom; i++) {
			do {
				test = levels[level_index];
				if (a[i] == test) {
					for (k = 0; k < zoom; k++)
						dotit(binrow + k, i);
					dir = 0;
					continue;
				}
				if (a[i] < test) {
					if ((a[i+1] > test) || (b[i] > test)) {
						for (k = 0; k < zoom; k++)
							dotit(binrow + k, i);
						dir = 0;
						continue;
					} 
					else
						dir = -1;
				}
				else {
					if ((a[i+1] < test) || (b[i] < test)) {
						for (k = 0; k < zoom; k++)
							dotit(binrow + k, i);
						dir = 0;
						continue;
					}
					else
						dir = 1;
				}
			} while (change_level(dir));
		}	
		temp = b;		/* swap the buffers */
		b = a;
		a = temp;
	}

	strcpy(fname, argv[1]);			/* record the image to the manager */
	store_window(fname, x_origin, y_origin, zoom, sr, endrow, sc, endcol, bin);
	(*screen_puts)(SCREEN_PAUSE,"");
}


	/* take a string of numbers separated by commas, like this,
	         100,101,300,500,...
	   and get the numbers out of it. Put them into the array
	   called 'array' */

static int
crack_levels(string, array)
char *string;
int *array;
{
	int i = 0, j = 0, start = 0;

	if (string[0] == 0)
		error(-1, "no levels specified after LEVEL keyword");

	do {
		i++;
		if ((string[i] == ',') || (string[i] == 0)) {
			array[j] = atoi(&string[start]);
			start = i + 1;
			j++;
			if (j == NLEVELS_MAX)
				error(-1, "too many levels");
		}
	} while(string[i]);
	return(j);
}

	/* generate evenly-spaced levels, starting at 'low' and separated
	   by a constant difference of 'diff', to go into 'array'. return the 
	   number of levels generated. */

static int
gen_even_levels(low, diff, array)
int low, diff, *array;
{
	int i, n;
	long test;

	if (got_levels && (got_diff || got_low))
		error(-1, "you can't have levels specified by level= and low=, diff=");
	if (!got_diff)
		error(-1, "no difference specified");
	if (!got_low)
		error(-1, "no low value specified");
			
	test = (long) INT_MAX/diff;
	if (test > NLEVELS_MAX) {
		(*screen_puts)(SCREEN_ADD, "warning: not all contour levels will be shown\n");
		n = NLEVELS_MAX;
	}
	else
		n = (int) test;
	
	for (i = 0; i < n; i++)
		array[i] = low + i*diff;
	
	return(n);
}

	/* generate evenly-ratioed levels, starting at 'low' and separated
	   by a constant multiplicative factor of 'ratio', to go into 'array'. 
	   return the number of levels generated. */

static int
gen_ratio_levels(low, ratio, array)
int low, *array;
double ratio;
{
	int i;
	long lev;

	if (got_levels && (got_ratio || got_low))
		error(-1, "you can't have levels specified by level= and low=, ratio=");
	if (!got_ratio)
		error(-1, "no ratio specified");
	if (!got_low)
		error(-1, "no low value specified");

	lev = (long) low;
	for (i = 0; ((i < NLEVELS_MAX) && (lev < (long) INT_MAX)); i++) {
		array[i] = (int) lev;
		lev *= ratio;
	}
	if (i == NLEVELS_MAX)
		(*screen_puts)(SCREEN_ADD, "warning: not all contour levels will be shown\n");

	return(i);
}

	/* keep track of whether pixel values are increasing or decreasing,
	   so that we know when to draw a contour. change the value of
	   'level_index' as needed, and return the new value of 'dir'. */

static int 
change_level(dir)
int dir;
{
	static old_dir;

	if (dir == 1) {
		if (old_dir == -1) {			/* if direction flips cancel */
			dir = 0;
		}
		else {
			if (level_index < (nlevels - 1)) {	/* increment level if poss */
				level_index++;
			}
			else
				dir=0;
		}
	}
	else {
		if (dir == -1) {
			if (old_dir == 1) {			/* prohibit direction flipping */
				dir = 0;
			}
			else {			
				if (level_index) {		/* decrement if possible */
					level_index--;
				}
				else
					dir = 0;
			}
		}
	}
	old_dir = dir;
	return(dir);
}


static void
dotit(x, y)
int x, y;
{
	extern int image_matrix[];
	int cx, cy;

	cx = x_origin + image_matrix[0]*x + image_matrix[1]*y;
	cy = y_origin + image_matrix[2]*x + image_matrix[3]*y;
	if (!hard_flag)
		(*image_dot)(cx, cy);
}
	
		
		/* rebin the data in the array 'buf' to either
		   crunch it down (if bin > 1) or expand it out
		   (if zoom > 1). or, if neither is > 1, do 
		   nothing. */

		/* assumes that 'buf' is at least NMAX in size */

static void
rebin_data(bin, zoom, nc, buf)
int bin, zoom, nc, *buf;
{
	int i, j, l, top;

	if (bin > 1) {
		for (i = 1, j = bin; i < nc; i++, j += bin)
			buf[i] = buf[j];
	}
	else if (zoom > 1) {
		if ((top = nc*zoom - 1) >= NMAX)
			top = ((NMAX/zoom)*zoom) - 1;	

		for (i = top, j = top/zoom; j >= 0; j--)
			for (l = 0; l < zoom; l++)
				buf[i--] = buf[j];

		/* may be a few rows left between 'top' and NMAX if 
		   nc*zoom > NMAX, so take care of them. */
		if (top + zoom >= NMAX) {
			j = (top + 1)/zoom;
			for (i = top + 1; i < NMAX; i++)
				buf[i] = buf[j];
		}
	}
}