#include "i3_des_cat.h"
#include "i3_logging.h"
#include "fitsio.h"
#include "i3_math.h"
#include "stdbool.h"


static const int DES_CATALOG_HDU=3;

i3_des_catalog * i3_des_catalog_stars(char * filename, float class_star_threshold){

	i3_print(i3_verb_debug,"Loading stars (threshold=%f) from file %s",class_star_threshold,filename);

	//filename is a FITS file.
	//open the fits file.
	fitsfile * input_file=NULL;
	int status=0;
	fits_open_file(&input_file, filename, READONLY, &status);
	if (status){
		i3_print(i3_verb_quiet,"Unable to open FITS file %s",filename);
		fits_report_error(stdout,status);
		return NULL;
	}
	
	//go to the right extension
	long total_rows;
	int hdu=DES_CATALOG_HDU;
	fits_movabs_hdu(input_file, hdu, NULL, &status);
	
	fits_get_num_rows(input_file, &total_rows, &status);
	int class_star_column=0;
	fits_get_colnum(input_file, CASEINSEN, "CLASS_STAR", &class_star_column, &status);
	if (status){
		fprintf(stderr,"Unable to read correct DES cat structure in FITS file %s",filename);
		fits_report_error(stderr,status);
		fits_close_file(input_file,&status);
		return NULL;
	}
	
	LONGLONG firstrow=1;
	LONGLONG firstelem=1;
	int anynul=0;
	float * class_star = malloc(sizeof(float)*total_rows);
	fits_read_col(input_file, TFLOAT, class_star_column, firstrow, firstelem, 
		total_rows, NULL, class_star, &anynul, &status);
		
	if (status){
		fprintf(stderr,"Unable to read correct DES cat columns in FITS file %s",filename);
		fits_report_error(stderr,status);
		fits_close_file(input_file,&status);
		return NULL;
	}
		
	
	//count the number of elements where class_star>threshold
	int number_stars=0;
	for(int i=0;i<total_rows;i++) if (class_star[i]>class_star_threshold) number_stars++;
	//create the catalog of the right size
	i3_print(i3_verb_extreme,"Found %d stars out of %d objects (threshold=%f) from file %s",number_stars,total_rows,class_star_threshold,filename);
	i3_des_catalog * cat = i3_des_catalog_create(number_stars);

	if (number_stars==0){
		i3_print(i3_verb_standard,"No stars beating threshold %f found in file %s",class_star_threshold,filename);
		fits_close_file(input_file,&status);
		return cat;
	}

	//pull out the right columns into arrays
	int column;
	float * x = malloc(sizeof(float)*total_rows);
	float * y = malloc(sizeof(float)*total_rows);
	float * xmin = malloc(sizeof(float)*total_rows);
	float * xmax = malloc(sizeof(float)*total_rows);
	float * ymin = malloc(sizeof(float)*total_rows);
	float * ymax = malloc(sizeof(float)*total_rows);
	
	fits_get_colnum(input_file, CASEINSEN, "X_IMAGE", &column, &status);
	fits_read_col(input_file, TFLOAT, column, firstrow, firstelem, 
		total_rows, NULL, x, &anynul, &status);
		
	fits_get_colnum(input_file, CASEINSEN, "Y_IMAGE", &column, &status);
	fits_read_col(input_file, TFLOAT, column, firstrow, firstelem, 
		total_rows, NULL, y, &anynul, &status);
		
	fits_get_colnum(input_file, CASEINSEN, "XMIN_IMAGE", &column, &status);
	fits_read_col(input_file, TFLOAT, column, firstrow, firstelem, 
		total_rows, NULL, xmin, &anynul, &status);
		
	fits_get_colnum(input_file, CASEINSEN, "XMAX_IMAGE", &column, &status);
	fits_read_col(input_file, TFLOAT, column, firstrow, firstelem, 
		total_rows, NULL, xmax, &anynul, &status);
		
	fits_get_colnum(input_file, CASEINSEN, "YMIN_IMAGE", &column, &status);
	fits_read_col(input_file, TFLOAT, column, firstrow, firstelem, 
		total_rows, NULL, ymin, &anynul, &status);
		
	fits_get_colnum(input_file, CASEINSEN, "YMAX_IMAGE", &column, &status);
	fits_read_col(input_file, TFLOAT, column, firstrow, firstelem, 
		total_rows, NULL, ymax, &anynul, &status);
		
	if (status){
		fprintf(stderr,"Unable to read correct DES col names in FITS file %s",filename);
		fits_report_error(stderr,status);
		fits_close_file(input_file,&status);
		return NULL;
	}
		
		
	int row=0;
	for (int i=0;i<total_rows;i++){
		if (class_star[i]>class_star_threshold){
			cat->x[row] = x[i];
			cat->y[row] = y[i];
			i3_flt xrange = xmax[i]-xmin[i];
			i3_flt yrange = ymax[i]-ymin[i];
			cat->size[row] = xrange>yrange ? xrange : yrange; //max of the two sizes
			row++;
		}
	}
	
	fits_close_file(input_file,&status);
	return cat;
}

int i3_des_catalog_save(i3_des_catalog * cat, char * filename){
	i3_print(i3_verb_debug,"Saving catalog to file %s",filename);
	FILE * f = fopen(filename,"w");
	if (!f){
		i3_print(i3_verb_quiet,"Failed to save catalog to file %s",filename);
		return I3_PATH_ERROR;
	}
	i3_des_catalog_print(cat,f);
	fclose(f);
	return I3_OK;
}

void i3_des_catalog_print(i3_des_catalog * cat, FILE * output_file)
{
	fprintf(output_file,"#index        x            y            size\n");
	for (int i=0;i<cat->n;i++){
		fprintf(output_file,"%-12d %-12f %-12f %-12f\n",i,cat->x[i],cat->y[i],cat->size[i]);
	}
}

i3_des_catalog * i3_des_catalog_create(int n)
{
	i3_des_catalog * cat = malloc(sizeof(i3_des_catalog));
	cat->n=n;
	if (n>0){
		cat->x = malloc(sizeof(i3_flt)*n);
		cat->y = malloc(sizeof(i3_flt)*n);
		cat->size = malloc(sizeof(i3_flt)*n);
	}
	else{
		cat->x=NULL;
		cat->y=NULL;
		cat->size=NULL;
	}
	return cat;
}


i3_image * i3_des_catalog_get_raw_stamp(i3_des_catalog * cat, i3_image * image, int i)
{
	return i3_des_catalog_get_stamp(cat, image, i, 0);
}

i3_image * i3_des_catalog_get_stamp(i3_des_catalog * cat, i3_image * image, int i, int min_size)
{
	
	// Some checks just in case something is wrong.
	if (cat==NULL){
		i3_print(i3_verb_standard,"Passed a NULL catalog into i3_des_catalog_get_image");
		return NULL;
	}
	
	if (i<0 || i>=cat->n){
		i3_print(i3_verb_standard,"Tried to get non-existent image numbered %d from catalog with n=%d",i,cat->n);
		return NULL;
	}
	
	i3_print(i3_verb_extreme,"Cutting catalog %d of %d image out of image. x,y,size=%f,%f,%f",i,cat->n,cat->x[i],cat->y[i],cat->size[i]);
	
	
	//Now cut the object out, making sure not to cross over the edges.
	int x_size =  (int) cat->size[i];
	int y_size = (int) cat->size[i];
	
	if (x_size<min_size) x_size=min_size;
	if (y_size<min_size) y_size=min_size;
	
	int x_start = (int) floor(cat->x[i] - x_size/2.);
	int y_start = (int) floor(cat->y[i] - y_size/2.);
	
	

	
	if (x_start<0) x_start=0;
	if (y_start<0) y_start=0;

	if (x_start+x_size>=image->nx) x_size = image->nx-x_start-1;
	if (y_start+y_size>=image->ny) y_size = image->ny-y_start-1;
	i3_print(i3_verb_crazy,"Final cut-out size %d %d %d %d",x_start,y_start,x_size,y_size);
	
	i3_image * result = i3_image_copy_part(image, x_start, y_start, x_size, y_size);
	if (!result){
		i3_print(i3_verb_debug,"Image catalog cut out failed - edge cross or similar.  (%d  %d) + (%d  %d).  Image dim (%d  %d)",x_start,y_start,x_size,y_size,image->nx,image->ny);
	}
	return result;
}