/* try to correlate two lists of stellar positions, finding the offset
	   which produces the most matches between stars in the two lists */

#include <stdio.h>
#include <math.h>
#include "multipht.h"


	/* we use this struct to sort stars by mag in each image */
struct s_sort {
	int index;
	double mag;
};
static struct s_sort *sort_index[MAX_IMAGE];

static struct s_star temp[MAX_STAR];

static int find_offset( /* int nim1, nim2, double *dx, *dy */ );
static int quick_align( /* int nim1, nim2; double x, y, *dr, *dc */ );
static int try_align( /* struct s_image *im1, *im2; int max; 
                                  double x, y, *dr, *dc */ );
static int is_a_match( /* struct s_star *s1, *s2; double *dr, *dc */ );
static void shift_image( /* struct s_image *im; double dx, dy */ );
static void create_sort_arrays();
static int comp_ssort( /* struct s_sort *s1, *s2 */ );





	/* match up all images, finding the offset which brings them into 
	   register with all the others.  use the given file as the
	   template, if it is not ""; otherwise, use image with most stars */

do_match(template)
char *template;
{
	int i, max, max_image;
	double dx, dy;

	max_image = -1;
	/* find the exposure with the largest number of stars, and try to 
	   match all other exposures to it. */
	if (strcmp(template, "") == 0) {
		max = 0;
		for (i = 0; i < num_images; i++) {
			if (image_arr[i].num_star > max) {
				max = image_arr[i].num_star;
				max_image = i;
			}
		}
	}
	else {
		for (i = 0; i < num_images; i++) {
			if (strcmp(image_arr[i].name, template) == 0) {
				max_image = i;
				break;
			}
		}
		if (max_image == -1)
			max_image = 0;
	
	}

	create_sort_arrays();

	if (debug_flag)
		printf("about to match all to image %d %s: %d stars\n", 
				max_image, image_arr[max_image].name, 
				image_arr[max_image].num_star);	
	max = 0;
	for (i = 0; i < num_images; i++) {
		if (i != max_image) {
			if (debug_flag)
				printf("about to find offset for image %d of %d\n", 
								i, num_images);
			image_arr[i].num_match = find_offset(max_image, i,
						&dx, &dy);
			if (image_arr[i].num_match > max)
				max = image_arr[i].num_match;
			shift_image(&(image_arr[i]), dx, dy);
		}
	}
	image_arr[max_image].num_match = max;

}



	/* find the best offset between the two given images.
	   when found, print it out. 
	   return the number of stars matched, and place the best offset
	   factors in 'dx' and 'dy' */

static int
find_offset(nim1, nim2, dx, dy)
int nim1, nim2;
double *dx, *dy;
{
	int i, j, max, n, possible;
	double dr, dc, bestdr, bestdc, x, y, bestx, besty;
	struct s_star *star1, *star2;
	struct s_image *im1, *im2;

	im1 = &(image_arr[nim1]);
	im2 = &(image_arr[nim2]);

	max = 0;
	possible = (im1->num_star < im2->num_star ? im1->num_star : im2->num_star);

	x = 0.0;
	y = 0.0;
	n = try_align(im1, im2, x, y, max, &dr, &dc);
	max = 0;
	bestdr = 0.0;
	bestdc = 0.0;

	for (i = 0; i < im1->num_star; i++) {
		if (i > MAXMATCH) {
			break;
		}
		star1 = &(im1->stars[sort_index[nim1][i].index]);
		for (j = 0; j < im2->num_star; j++) {
			if (j > MAXMATCH) {
				break;
			}
			star2 = &(im2->stars[sort_index[nim2][j].index]);
			x = star1->xc - star2->xc;
			y = star1->yc - star2->yc;
			n = quick_align(nim1, nim2, x, y, &dr, &dc);
			if (debug_flag) {
				printf("aligning %4d (%5.0lf, %5.0lf), %4d (%5.0lf, %5.0lf) ",
						sort_index[nim1][i].index, star1->xc, star1->yc,
						sort_index[nim2][j].index, star2->xc, star2->yc);
				printf(" -> %4d ", n);
			}
			if (n > max) {
				max = n;
				bestdr = dr;
				bestdc = dc;
				bestx = x;
				besty = y;
				if (debug_flag) {
					printf("best!");
				}
			}
			if (debug_flag) {
				printf("\n");
			}
		}
	}
	if (max == 0)
		return(0);
	else {
		n = try_align(im1, im2, bestx, besty, max, &dr, &dc);
		printf("lists %s %s dr = %6.2lf dc = %6.2lf  %4d of %4d\n",
					im1->name, im2->name, bestdr, bestdc, n, possible);
		*dx = bestdr;
		*dy = bestdc;
		return(max);
	}
}

	/* using the given shifts in row and col,
	   compare the brightest MAXMATCH stars only in both lists to find
	   the number of other matches.  return the number of 'good' matches
	   (stars which fall within one critical distance of each other),
	   as well as the average 'dr' and 'dc' for all the good matches.

	   returning 0 means that there were no matches, which shouldn't
	   occur, since we explicitly match the two given.  */

static int
quick_align(nim1, nim2, x, y, dr, dc)
int nim1, nim2;
double x, y, *dr, *dc;
{
	int i, j, k, num_match, bestj;
	double z, sumdr, sumdc, tdr, tdc, bestdr, bestdc;
	struct s_image *im1, *im2;

	im1 = &(image_arr[nim1]);
	im2 = &(image_arr[nim2]);

	for (j = 0; j < im2->num_star; j++) {
		if (j > MAXMATCH) {
			break;
		}
		temp[j].xc = im2->stars[sort_index[nim2][j].index].xc + x;
		temp[j].yc = im2->stars[sort_index[nim2][j].index].yc + y;
	}
		
	sumdr = 0.0;
	sumdc = 0.0;
	num_match = 0;
	for (i = 0; i < im1->num_star; i++) {
		if (i > MAXMATCH) {
			break;
		}
		bestdr = BIG;
		bestdc = BIG;
		bestj = 0;
		for (j = 0; j < im2->num_star; j++) {
			if (j > MAXMATCH) {
				break;
			}
			if (is_a_match(&(im1->stars[sort_index[nim1][i].index]), 
			                &(temp[j]), &tdr, &tdc)) {
				if (tdr*tdr + tdc*tdc < bestdr*bestdr + bestdc*bestdc) {
					bestdr = tdr;
					bestdc = tdc;
					bestj = j;
				}
			}
		}
		if (bestdr < BIG) {
			sumdr += bestdr;
			sumdc += bestdc;
			num_match++;
		}
	}

	if (num_match > 0) {
		*dr = x + sumdr/(double)num_match;
		*dc = y + sumdc/(double)num_match;
	}
	return(num_match);
}




	/* using the given shifts in row and col,
	   compare all the other stars in both lists to find
	   the number of other matches.  return the number of 'good' matches
	   (stars which fall within one critical distance of each other),
	   as well as the average 'dr' and 'dc' for all the good matches.

	   the maximum number 
	   matched so far is 'max'.  Use this number in an intelligent
	   fashion to realize when going further is not going to beat the
	   best match. 

	   returning 0 means that there were no matches, which shouldn't
	   occur, since we explicitly match the two given.  */

static int
try_align(im1, im2, x, y, max, dr, dc)
int max;
double x, y, *dr, *dc;
struct s_image *im1, *im2;
{
	int i, j, k, num_match, bestj, giveupi;
	double z, sumdr, sumdc, tdr, tdc, bestdr, bestdc;

	giveupi = im1->num_star - max;

	for (i = 0; i < im2->num_star; i++) {
		temp[i].xc = im2->stars[i].xc + x;
		temp[i].yc = im2->stars[i].yc + y;
	}
		
	sumdr = 0.0;
	sumdc = 0.0;
	num_match = 0;
	for (i = 0; i < im1->num_star; i++) {
		bestdr = BIG;
		bestdc = BIG;
		bestj = 0;
		for (j = 0; j < im2->num_star; j++) {
			if (is_a_match(&(im1->stars[i]), &(temp[j]), &tdr, &tdc)) {
				if (tdr*tdr + tdc*tdc < bestdr*bestdr + bestdc*bestdc) {
					bestdr = tdr;
					bestdc = tdc;
					bestj = j;
				}
			}
		}
		if (bestdr < BIG) {
			sumdr += bestdr;
			sumdc += bestdc;
			num_match++;
		}
		if (i - num_match > giveupi)
			break;
	}

	if (num_match > 0) {
		*dr = x + sumdr/(double)num_match;
		*dc = y + sumdc/(double)num_match;
	}
	return(num_match);
}

	/* given two stars, find the difference in their row and col positions,
	   placing the results in the passed parameters 'dr' and 'dc'.  Then
	   compare the distance between the two stars to the critical distance for
	   a match; if the distance is less (and we have a match), return 1.
	   otherwise, return 0. */

static int
is_a_match(s1, s2, dr, dc)
struct s_star *s1, *s2;
double *dr, *dc;
{
	*dr = (s1->xc - s2->xc);
	*dc = (s1->yc - s2->yc);
	if (((*dr)*(*dr) + (*dc)*(*dc)) <= RAD2) 
		return(1);
	else
		return(0);
} 

	/* shift all the stars in the given image by the given offsets */

static void
shift_image(im, dx, dy)
struct s_image *im;
double dx, dy;
{
	int i;

	for (i = 0; i < im->num_star; i++) {
		im->stars[i].xc += dx;
		im->stars[i].yc += dy;
	}
}

	/* 
	 * create an array of indices for each image; sort the array
	 * of indices so that sort_index[i][j] is the index of the j'th
	 * brightest star in image [i].
	 *
	 * Use the zero'th aperture mag for sorting.
	 */

static void
create_sort_arrays()
{
	int i, j;
	struct s_star *sp;

	for (i = 0; i < num_images; i++) {


		/* allocate space for the sort_index for this image */
		if ((sort_index[i] = (struct s_sort *) 
		       malloc(image_arr[i].num_star*sizeof(struct s_sort))) == NULL) {
			fprintf(stderr, 
                    "create_sort_array: can't alloc for index array %d\n", i);
			exit(1);
		}

		/* fill the sort_index with data for this image */
		for (j = 0; j < image_arr[i].num_star; j++) {
			sort_index[i][j].index = j;
			sort_index[i][j].mag = image_arr[i].stars[j].mag[0];
		}

#if 1
		printf("before sorting, first 5 are \n");
		for (j = 0; j < 5; j++) {
			printf("%5d %7.2lf\n", j, sort_index[i][j].mag);
		}
#endif

		/* sort the sort_index, placing bright stars in the first elements */
		qsort((char *) sort_index[i], image_arr[i].num_star,
		               sizeof(struct s_sort), comp_ssort);

#if 1
		printf("after  sorting, first 5 are \n");
		for (j = 0; j < 5; j++) {
			printf("%5d %7.2lf\n", j, sort_index[i][j].mag);
		}
#endif

	}

}
		

	/* 
	 * return -1 if the first star has a brighter mag than the second
	 *         0                        the same 
	 *         1                         fainter mag  
	 */

static int
comp_ssort(s1, s2)
struct s_sort *s1, *s2;
{
        if (s1->mag < s2->mag)
                return(-1);
        else if (s1->mag > s2->mag)
                return(1);
        else
                return(0);
}