/* assign unique IDs to all stars, based on their aligned positions */

#include <stdio.h>
#include <math.h>
#include <search.h>
#include "bait.h"
#include "pf.h"
#include "multipht.h"

#define HIST_BINS    20			/* number of bins in histogram */

static struct s_hash *hashtab;		/* hash table array */
static int nhash = 0;				/* number of hash table entries */

static void sort_stars( /* struct s_image *im */ );
static int comp_by_xc( /* struct s_star *s1, *s2 */ );
struct s_hash * lookup( /* struct s_star *star */ );
static struct s_hash * looklist( /* struct s_star *star; int index */ );
static void hash_insert( /* struct s_star *star */ );
static void hash_add( /* struct s_star *star; struct s_hash *p */ );
static void add_id( /* struct s_hash *p */ );
static void init_hash( /* */ );
static int is_a_match( /* struct s_star *s1; struct s_hash *p */ );
static void print_hash( /* */ );
static void print_by_id( /* */ );

	/* assign unique IDs to each star in each image, based upon its
	   position in the template coordinate system. 
	   Note that the stars should already be almost - but not
	   quite - sorted by their rows, so we have to sort them first.  

	   as each star is examined, the 'id' field of its structure is set
	   the master ID number for all stars in all images.  Moreover,
	   the hash table entry for that ID has the new star added to 
	   its linked list of stars that share the ID, and has its
	   'num' field incremented.
	   thus, we can go through all stars BY ID efficiently after 
	   having built up the table, and we can quickly find out how
	   many stars share the same ID (i.e. in how many images was
	   each star found).

	   return the total number of unique stars found */

assign_ids()
{
	int i, j;
	struct s_hash *p; 
	struct s_star *s;

	init_hash();
	for (i = 0; i < num_images; i++) {
		sort_stars(&(image_arr[i]));
		for (j = 0; j < image_arr[i].num_star; j++) {
			s = &(image_arr[i].stars[j]);
			if ((p = lookup(s)) != NULL)
				hash_add(s, p);
			else
				hash_insert(s);
		}
	}

	if (debug_flag) {
/****
		print_by_id();
		print_hash();
		print_star_hist();
		print_image_hist();
****/
	}

	return(num_ids);
}

	/* sort all the stars in the given image by their 'xc' value */

static void
sort_stars(im, how)
struct s_image *im;
{
	qsort((char *) &(im->stars[0]), im->num_star, sizeof(struct s_star),
				comp_by_xc);
}

	/* return -1 if the first star has an 'xc' coord less than second
	           0                                     equal to 
	          +1                                     greater than */

static int
comp_by_xc(s1, s2)
struct s_star *s1, *s2;
{
	if (s1->xc < s2->xc)
		return(-1);
	else if (s1->xc > s2->xc)
		return(1);
	else
		return(0);
}
	

	/* lookup a star in the hash table, using its 'xc' coord as a key.
	   return a pointer to the hash table entry, if found, or NULL 
	   if not found. */

struct s_hash *
lookup(star)
struct s_star *star;
{
	int i, index;
	struct s_hash *found;

	index = (int) star->xc;

	/* look first in the linked list of stars with the same 'xc' value */
	if ((found = looklist(star, index)) != (struct s_hash *) NULL)
		return(found);
	
	/* now look in the linked lists of stars with progressively smaller xc
	   values, down to 'xc - 2*RAD'. */
	for (i = index - 1; (i >= 0) && (i >= index - (int)(2*RADIUS)); i--)
		if ((found = looklist(star, i)) != (struct s_hash *) NULL)
			return(found);

	/* now look in the linked lists of stars with progressively larger xc
	   values, up to 'xc + RAD'. */
	for (i = index + 1; (i < nhash) && (i <= index + (int)(2*RADIUS)); i++)
		if ((found = looklist(star, i)) != (struct s_hash *) NULL)
			return(found);

	return((struct s_hash *) NULL);
}

	/* look for a star in the given linked list of hash entries. if found,
	   return a pointer to the hash entry.  if not, return NULL. */

static struct s_hash *
looklist(star, index)
struct s_star *star;
int index;
{
	struct s_hash *p;

	if ((index < 0) || (index >= nhash))
		return((struct s_hash *) NULL);
	for (p = &(hashtab[index]); (p != NULL) && (p->id != -1); p = p->next) {
		if (is_a_match(star, p))
			return(p);
	}

	return((struct s_hash *) NULL);
}

	/* insert a star into the hash table */

static void
hash_insert(star)
struct s_star *star;
{
	int index, flag;
	struct s_hash *p, *q;

	index = (int) star->xc;
	if (hashtab[index].id == -1)
		flag = 1;

	if (flag == 1) {
		p = &(hashtab[index]);
	}
	else {
		if ((p = (struct s_hash *) malloc(sizeof(struct s_hash))) == NULL) {
			fprintf(stderr, "hash_insert: can't malloc for new entry, num_ids=%d\n",
					num_ids);
			exit(-1);
		}
		p->num = 0;
	}
	p->xc = star->xc;
	p->yc = star->yc;
	p->id = num_ids++;
	star->id = p->id;
	p->num++;
	if ((p->starp = (struct s_starp *) malloc(sizeof(struct s_starp))) == NULL) {
		fprintf(stderr, "hash_insert: can't malloc for starp, num_ids=%d\n", 
							num_ids);
		exit(-1);
	}
	p->starp->star = star;
	p->starp->next = (struct s_starp *) NULL;
	p->next = (struct s_hash *) NULL;
	p->include = 1;

	if (flag == 0) {
		/* walk down the linked list of hash structures to the end, then
		   put the new structure at the end of the list */	
		for (q = &(hashtab[index]); q->next != NULL; q = q->next)
			;
		q->next = p;
	}

	add_id(p);
}
	

	/* add a star to the hash entry that matches it.  modify the hash entry's
	   'xc' and 'yc' values so that they are the average of all matching
	   star values. */

static void
hash_add(star, p)
struct s_star *star;
struct s_hash *p;
{
	double dn;
	struct s_starp *sp, *q;

	if ((sp = (struct s_starp *) malloc(sizeof(struct s_starp))) == NULL) {
		fprintf(stderr, 
                  "hash_add: can't malloc for starp, num_ids=%d\n", num_ids);
		exit(-1);
	}
	sp->star = star;
	sp->next = (struct s_starp *) NULL;
	p->num++;
	star->id = p->id;
	for (q = p->starp; (q != NULL) && (q->next != NULL); q = q->next)
		;
	q->next = sp;
	
	/* now modify the hash entry's xc and yc values */
	dn = p->num;
	p->xc = ((p->xc*dn) + star->xc)/(dn + 1.0);
	p->yc = ((p->yc*dn) + star->yc)/(dn + 1.0);
}	

	/* add a new entry to the linked list of hash entries which
	   are in order by their ID numbers.  the head of this list
	   is 'star_list', and by walking down it, you can walk through
	   all hash entries (and, hence, stars) by their ID number */

static void
add_id(p)
struct s_hash *p;
{
	static int first = 0;
	static struct s_hashp *last;
	struct s_hashp *hp;

	if ((hp = (struct s_hashp *) malloc(sizeof(struct s_hashp))) == NULL) {
		fprintf(stderr, "add_id: can't malloc for hashp, num_ids=%d\n", num_ids);
		exit(-1);
	}
	hp->hash = p;
	hp->next = (struct s_hashp *) NULL;
	if (first == 0) {
		first = 1;
		star_list = hp;
		last = star_list;
	}
	else {	
		last->next = hp;
		last = last->next;
	}
}

	/* initialize the hash table */

static void
init_hash()
{
	int i, rowmax, rowmin, colmax, colmin;
	pf_handle pfp;
	char c1, c2, val[PF_LEN + 1], str[PF_LEN + 1];

	if ((pfp = pf_open(CONFIG_FILE, 1, PF_EXIST)) < 0) {
		fprintf(stderr, "init_hash: can't open config file %s\n", 
				CONFIG_FILE);
		exit(-1);
	}
	if ((pf_getval(pfp, "CCDSEC", val) == NULL) ||
	    (pf_valtostr(val, str) == PF_ERROR)) {
		fprintf(stderr, "init_hash: can't read CCDSEC from config file\n");
		pf_close(pfp);
		exit(-1);
	}
	pf_close(pfp);

	if (sscanf(str, "%c%d:%d,%d:%d%c", &c1, &rowmin, &rowmax, &colmin,
						&colmax, &c2) != 6) {
		fprintf(stderr, "init_hash: can't parse CCDSEC str ..%s..\n", str);
		exit(-1);
	}
	nhash = rowmax;

	/* allocate an array large enough to fit 'rowmax' rows */
	if ((hashtab = (struct s_hash *) calloc(nhash, sizeof(struct s_hash))) == NULL) {
		fprintf(stderr, "init_hash: can't calloc for %d hash entries \n", 
							nhash);
		exit(-1);
	}

	/* initialize the values in the table */
	for (i = 0; i < nhash; i++) {
		hashtab[i].xc = -1.0;
		hashtab[i].yc = -1.0;
		hashtab[i].id = -1;
		hashtab[i].num = 0;
		hashtab[i].starp = (struct s_starp *) NULL;
		hashtab[i].next = (struct s_hash *) NULL;
		hashtab[i].include = 0;
	}
}

	/* compare the distance between the star and the hash-entry-centroid 
	   to the critical distance for a match; if the distance is less,
	   we have a match: return 1.  otherwise, return 0. */

static int
is_a_match(s1, p)
struct s_star *s1;
struct s_hash *p;
{
	double dr, dc;

	dr = (s1->xc - p->xc);
	dc = (s1->yc - p->yc);
	if ((dr*dr + dc*dc) <= RAD2) 
		return(1);
	else
		return(0);
} 

	/* print out the entire hash table - for debugging purposes only */

static void
print_hash()
{
	int i, n, nn;
	struct s_hash *p;
	struct s_starp *s;

	for (i = 0; i < nhash; i++) {
		if (hashtab[i].id == -1) {
			printf("hashtab[%4d]   ", i);
			printf("empty \n", i);
		}
		else {
			n = 0;
			printf("hashtab[%4d]   ", i);
			p = &(hashtab[i]);
			while (p != NULL) {
				if (n == 0)
					printf("entry %3d: ID=%5d  num=%3d  xc=%6.2lf yc=%6.2lf inc=%d\n", 
							n, p->id, p->num, p->xc, p->yc, p->include);
				else 
					printf("%15s entry %3d: ID=%5d  num=%3d  xc=%6.2lf yc=%6.2lf inc=%d\n", 
							" ", n, p->id, p->num, p->xc, p->yc, p->include);
				nn = 0;
				s = p->starp;
				while (s != NULL) {
					printf("%20s starp=%d  xc=%6.2lf yc=%6.2lf reqid=%s\n",
								" ", nn, s->star->xc, s->star->yc, 
								s->star->image->name);
					nn++;
					s = s->next;
				}
				n++;
				p = p->next;
			}
		}
	}
}

	/* print out a nice list of stars by their ID numbers */

static void
print_by_id()
{
	int nn;
	struct s_hashp *hp;
	struct s_hash *p;
	struct s_starp *s;

	for (hp = star_list; hp != NULL; hp = hp->next) {
		p = hp->hash;
		printf("ID=%5d  num=%3d  xc=%6.2lf yc=%6.2lf inc=%d\n", 
						p->id, p->num, p->xc, p->yc, p->include);
		nn = 0;
		s = p->starp;
		while (s != NULL) {
			printf("%9s starp=%d  xc=%6.2lf yc=%6.2lf w3=%d w4=%-8.1lf  reqid=%s\n", 
						" ", nn, s->star->xc, s->star->yc, s->star->w3,
						s->star->w4, s->star->image->name);
			printf("%23s peak=%6.0lf fwhm=%5.2lf sharp=%4.2lf round=%5.2lf\n",
						" ", s->star->peak, s->star->fwhm, s->star->sharp,
						s->star->round);
			nn++;
			s = s->next;
		}
	}

}

	/* print out a histogram of the number of images in which each
	   star appears */

void
print_star_hist()
{
	int i, div, hist[HIST_BINS], max;
	struct s_hashp *hp;
	struct s_hash *p;

	for (i = 0; i < HIST_BINS; i++)
		hist[i] = 0;

	max = 0;
	for (hp = star_list; hp != NULL; hp = hp->next) {
		p = hp->hash;
		if (p->num > max)
			max = p->num;
	}
	div = ((max - 1)/HIST_BINS) + 1;
	for (hp = star_list; hp != NULL; hp = hp->next) {
		p = hp->hash;
		hist[(p->num/div)]++;
	}
	for (i = 0; i < HIST_BINS; i++)
		printf("%5d  %3d\n", (int) ((i+0.5)*div), hist[i]);

}

	/* print out a histogram of the number of images with at least N
	   stars, vs. N */

void
print_image_hist()
{
	int i, div, hist[HIST_BINS], max;

	for (i = 0; i < HIST_BINS; i++)
		hist[i] = 0;

	max = 0;
	for (i = 0; i < num_images; i++) {
		if (image_arr[i].num_star > max)
			max = image_arr[i].num_star;
	}
	div = ((max - 1)/HIST_BINS) + 1;
	for (i = 0; i < num_images; i++) {
		hist[(image_arr[i].num_star/div)]++;
	}
	for (i = 0; i < HIST_BINS; i++)
		printf("%5d  %3d\n", (int) ((i+0.5)*div), hist[i]);

}