/**************************************************************************
 *                                                                        *
 * 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.                        *
 *                                                                        *
 **************************************************************************/


/*	convolves data files with Gaussian, boxcar, or user specified kernels

    SMOOTH file [box] [fw=] [fwr=] [fwc=] [outfile=] [boxcar=] [verbose] [nodc=] [help]

	stolen from the FORTRAN code by Tod Lauer and Kate Ebneter

	11/24/92 repaired a bug which moved the image 1 row by removing the
			complex code needed to prevent pages swaps (which don't seem
			to occur on the Sun) -rrt
	8/19/93 add smooth: to error message -rrt
	7/17/95 force prototype
	1/8/96  - removed 'file=' option, since the new 'ring' program
	                allows users to supply their own kernals for smoothing 
	                --- MWR
	1/25/96 - oops; 'ring' program does something completely different
	                It really would be useful to have a 'file=' option,
	                after all.  Well, the old version wasn't actually
	                implemented, so I don't feel so bad.  MWR
*/
#include "pcvista.h"
#include "fits.h"
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>	/* for exit */

/* 
 * I don't know why we need to include this under SunOs 5.5
 */
#ifndef LINUX
extern double hypot(double x, double y);
#endif



/*	possible methods of specifying the convolution kernel
*/
#define NO_KERNEL               0
#define BOXCAR_KERNEL           1
#define GAUSSIAN_KERNEL         3
#define NODC_KERNEL             4

#define MAXHOLD 		50

typedef short IMAGE_DATA;		/* type of data used to hold image */
typedef float BUFFER_DATA;		/* type of data used to hold intermediaries */
typedef int KERNEL;

static int verbose;

int main(int, char **);
static int check_kernel(int, int);
static int gauss_kernel(double , double *);
static int boxcar_kernel(double , double *);
static int nodc_kernel(double , double ,KERNEL ***);
static void normalize_kernel( double *);
static void row_convolve(IMAGE_DATA **, int, int ,double *, int);
static void col_convolve( IMAGE_DATA **, int, int, double *, int);
static void both_convolve( IMAGE_DATA **, int, int, KERNEL **, int);
static void convolve(IMAGE_DATA *, BUFFER_DATA *, int, double *, int);

static char *usage = 
   "usage: smooth file [box] [fw=] [fwr=] [fwc=] [outfile=] [boxcar=] \n                   [verbose] [nodc=] [help]";

int
main(argc,argv)
int argc;
char *argv[];
{
	FITS_HANDLE fin, fout;
	static int i, boxno=(-1);	/* ask for entire image */
	int nrow, ncol, nr, nc;
	int npr, npc, npb;		/* number of points in kernels */
	IMAGE_DATA **data;
	double fw, fwr, fwc, fwb;
	double filling=0.5;		/* used in nodc */
	static int col_flag, row_flag, both_flag; 
	static int kernel_mode=NO_KERNEL;
	static int normalize_flag=1;
	static int got_in_name, got_out_name;
	KERNEL **funcb;
	char *gotit, in_name[NBUF], out_name[NBUF];
	static double funcr[MAXHOLD], funcc[MAXHOLD];

	if (argc == 1)
		error(-1, usage);

	/* placate compiler */
	npr = 0;
	npc = 0;
	npb = 0;
	fw = 0.0;
	fwr = 0.0;
	fwc = 0.0;
	fwb = 0.0;

	for(i = 1; i < argc; i++) {
		if (keyword("help", argv[i])) {
			printf("%s\n", usage);
			exit(0);
		}
		if ((gotit = find("box", argv[i])) != NULL) {
			boxno = (int)evaluate(gotit);
			continue;
		}
		if ((gotit = find("fw", argv[i])) != NULL) {
			fwr= fwc = evaluate(gotit);
			kernel_mode=check_kernel(kernel_mode, GAUSSIAN_KERNEL);
			row_flag=col_flag=1;
			continue;
		}
		if ((gotit = find("fwc", argv[i])) != NULL) {
			fwc = evaluate(gotit);
			kernel_mode=check_kernel(kernel_mode, GAUSSIAN_KERNEL);
			col_flag=1;
			continue;
		}
		if ((gotit = find("fwr", argv[i])) != NULL) {
			fwr = evaluate(gotit);
			kernel_mode=check_kernel(kernel_mode, GAUSSIAN_KERNEL);
			row_flag=1;
			continue;
		}
		if ((gotit = find("boxcar", argv[i])) != NULL) {
			fwb = evaluate(gotit);
			kernel_mode=check_kernel(kernel_mode, BOXCAR_KERNEL);
			row_flag=col_flag=1;
			continue;
		}
		if ((gotit = find("nodc", argv[i])) != NULL) {
			fw = evaluate(gotit);
			kernel_mode=check_kernel(kernel_mode, NODC_KERNEL);
			both_flag=1;
			normalize_flag=0;
			continue;
		}
		if ((gotit = find("filling", argv[i])) != NULL) {
			filling = evaluate(gotit);
			if(filling <= 0. || filling >=1.0)
				error(1,"smooth: filling must 0.< f < 1.0");
			continue;
		}
		if ((gotit = find("outfile", argv[i])) != NULL) {
			strcpy(out_name, gotit);
			got_out_name=1;
			continue;
		}
		if(keyword("verbose",argv[i])){
			verbose=1;
			continue;
		}
		if(!got_in_name++)
			strcpy(in_name,argv[i]);
		else{
			fprintf(stderr, "unknown argument '%s'\n",argv[i]);
			exit(1);
		}
	}
	if(!got_in_name)
		error(1,"smooth: no input file specified");
	if(!got_out_name)
		fout=fin = fits_open(in_name, "r+", &nrow, &ncol);
	else{
		fin = fits_open(in_name, "r", &nrow, &ncol);
		fout= fits_open(out_name, "w", &nrow, &ncol);
	}

	if( ! (row_flag || col_flag || both_flag ) )
		error(1,"smooth: neither rows nor columns were specified"); 
			

	switch(kernel_mode){
		case GAUSSIAN_KERNEL:
			npr=gauss_kernel(fwr, funcr);
			npc=gauss_kernel(fwc, funcc);
			break;
		case BOXCAR_KERNEL:
			npr=boxcar_kernel(fwb, funcr);
			npc=boxcar_kernel(fwb, funcc);
			break;
		case NODC_KERNEL:
			npb=nodc_kernel(fw, filling,&funcb);
			break;
		default:
			error(1,"smooth: unknown kernel");
	}

	if(normalize_flag){
		normalize_kernel(funcr);
		normalize_kernel(funcc);
	}
		

	data = boxdata( fin , boxno , nrow, ncol, &nr, &nc);	/* read data */

	if(col_flag) 		/* perform column convolution if needed */
		col_convolve(data , nr, nc, funcc, npc);

	if(row_flag)		/* perform row convolution if needed*/
		row_convolve(data , nr, nc, funcr, npr);
			
	if(both_flag)
		both_convolve(data, nr, nc, funcb, npb);

	for(i=0; i<nr; i++)
		fits_put_data(fout,i,0,data[i],nc);

	fits_close(fin);
	fits_close(fout);

	exit(0);
}

/*	checks conflicts with present kernel selection
*/

static int check_kernel(mode, new)
int mode,new;
{
	if( mode != NO_KERNEL && mode != new)
		error(1,"smooth: conflicting kernel modes ");

	return new;
}

/*	generates a one-dimensional gaussian kernel
*/

static int gauss_kernel(fw, f)
double fw, *f;
{
	int i, np, mid;
	double x, fac;

	if(fw < 0.){
		fprintf(stderr,"fwhm cannot be less than zero ");
		exit(1);
	}
	if( fw < SMALL){
		f[0]=0;
		return 1.;
	}
	np = 2 * ((3 * (int)fw) /2) +1;
	if( np >= MAXHOLD)
		error(1,"smooth: gaussian width too large");

	mid = np/2;

	fac=(-4.0) * log( 2.0) /(fw * fw);
	for(i=0 ; i< np; i++){
		x=(double) (i-mid);
		f[i]= exp(x*x*fac);
	}
	return np;
}		


static int boxcar_kernel(w, f)
double w, *f;
{
	int i, np;

	if(w < 0.){
		fprintf(stderr,"fwhm cannot be less than zero ");
		exit(1);
	}

	np=2*(int)((w + 0.9999)/2.0) + 1;

	for(i=0; i<np; i++)
		f[i]=1;	

	f[0]=f[np-1]=(w - (double)(np-2))/2.;	/* end point values */

	return np;
}

	
/*	normalize so that the sum == 1.0
*/

static void normalize_kernel(f)
double *f;
{
	int i;
	double sum=0.;
	
	for(i=0; i< MAXHOLD; i++)
		sum += f[i];

	for(i=0; i< MAXHOLD; i++)
		f[i]/= sum;
}

/*	convolve each image row with the column convolving function
	pack the result back into the image
*/
static void col_convolve(d, nr, nc, f, npc)
IMAGE_DATA **d;
int nr,nc;
double *f;			/* kernel */
int npc;			/* number of points in convolving kernel */
{
	int i, j;
	BUFFER_DATA *hold;

	hold=(BUFFER_DATA *) malloc(nc * sizeof(BUFFER_DATA));
	if(hold == NULL)
		error(1,"smooth: can't allocate column buffer");

	for(i=0; i < nr ; i++){
		convolve(d[i] , hold, nc, f , npc);
		for( j=0; j < nc; j++)
			d[i][j]=(IMAGE_DATA)(hold[j]+0.5);
	}
	free(hold);
}

/*	Do the row convolution.  
*/

static void row_convolve(d, nr, nc, f, npc)
IMAGE_DATA **d;
int nr,nc;
double *f;			/* kernel */
int npc;			/* number of points in convolving kernel */
{
	int i, j;
	BUFFER_DATA *hold;
	IMAGE_DATA *im;

	hold=(BUFFER_DATA *) malloc(nr * sizeof(BUFFER_DATA));
	if(hold == NULL)
		error(1,"smooth: can't allocate row buffer");

	im=(IMAGE_DATA *) malloc(nr * sizeof(IMAGE_DATA));
	if(hold == NULL)
		error(1,"smooth: can't allocate row buffer2");

	for(i=0; i < nc ; i++){
		for(j=0; j<nr; j++)
			im[j]=d[j][i];

		convolve(im , hold, nr, f , npc);
		for( j=0; j < nr; j++)
			d[j][i]=(IMAGE_DATA)(hold[j]+0.5);
	}
	free(hold);
}

/*	performs the 1 dimensional convolution
		with normalization correction at edges
*/

static void convolve(d , b, n, psf, np)
IMAGE_DATA *d;			/* input data vector */
BUFFER_DATA *b;			/* output convolved data */
int n;					/* number of points in data array */
double *psf;			/* convolving kernel */
int np;					/* number of point in kernel array */
{
	int i, k, mid;
	int js, je;
	double sum;

	mid= np/2;

	for(i=0; i< n; i++ ){
		b[i]=(BUFFER_DATA) 0;			/* clear the output array */		
		js= i-mid;			/* find the limits of the convolution */
		je= js+np;
		if( js < 0) {		/* renormalize on edges */
			for(sum=1.,k=0 ; js < 0 ; js++,k++)
				sum -= psf[k];
			while(js< je)
				b[i] += psf[k++]*d[js++];
			if(sum > SMALL) 
				b[i] /= sum;
		}else{
			if(je >= n){
				for( k=0; js < n; js++,k++)
					b[i] += psf[k]*d[js];
				for(sum=1. ; js < je ; js++,k++)
					sum -= psf[k];
				if(sum > SMALL) 
					b[i] /= sum;
			}else{
				for( k=0; js < je; js++,k++)
					b[i] += psf[k]*d[js];
			}
		}
	}
}

/*	generates the kernel function with no net content
		inputs: full width of smoothing function
				pointer to array of KERNEL for kernel
		return: the size of the kernel

	the kernel function has an odd number of points so that it
	doesn't move the centroid.  It has a gaussian form which goes
	to zero at the half points, and the rest of the array is
	adjusted so that there is a zero mean.
*/

#define SCALE 	1000
#define FILL_FACTOR 1.0

static int nodc_kernel(fw, filling , p)
double fw,filling;
KERNEL ***p;
{
	int i,j,mid,np;
	static KERNEL *kernel[MAXHOLD],isum;
	double *temp[MAXHOLD];
	double fac,norm,dist,y;
	double sump=0., summ=0.;

	if(fw < 0.){
		fprintf(stderr,"fw cannot be less than zero\n");
		exit(1);
	}
	np = 2 * (int)(FILL_FACTOR*fw) +1;
	if(np >= MAXHOLD)
		error(1,"smooth: kernel too large ");

	for(i=0; i<np; i++){
		temp[i]=(double *) malloc(np * sizeof(double));
		kernel[i]=(KERNEL *) malloc(np * sizeof(KERNEL));
		if(temp[i] == NULL || kernel[i]== NULL)
			error(1,"smooth: can't allocate dc_kernel2");
	}

	fac=(-4.0) * log( 2.0) /(fw * fw);
	mid=np/2;

	for(i=0; i<np; i++){
		for(j=0; j< np; j++){
			dist=hypot((double)(i-mid), (double)(j-mid));
			y=SCALE*(exp(dist*dist*fac)-filling);
			if(y>0.){
				sump+=y;
			}else{
				summ-=y;
			}
			temp[i][j]=(KERNEL)(y);
		}
	}	
	/* now normalize so that it has 0.0 sum */
	norm=sump/summ;
	if(norm <=0.)
		error(1,"smooth: error normalizing .. bad filling factor");

	isum=0;
	for(i=0; i<np; i++){
		for(j=0; j < np; j++){
			if(temp[i][j]< 0.)
				temp[i][j]*=norm;
			kernel[i][j]=(KERNEL)(temp[i][j]+0.5);
			isum+=kernel[i][j];
		}
		free(temp[i]);
	}
	kernel[mid][mid]-=isum;		/* adjust central peak to get zero sum */

	if(verbose){
		for(i=0; i<np; i++){
			for(j=0; j < np; j++){
				printf("%7d ",kernel[i][j]);
			}
			printf("\n");
		}
		fflush(stdout);
	}

	*p=kernel;
	return np;
}

/*	does the two dimensional convolution
		edges set to zero
*/

static void both_convolve(d, nr, nc, f, n)
IMAGE_DATA **d;
int nr,nc;
KERNEL **f;
int n;
{
	int i,j,half,row,col;
	KERNEL **hold,*acc;
	IMAGE_DATA *pd;

	hold=(KERNEL **)malloc(nr * sizeof(KERNEL *));
	if( hold == NULL)
		error(1,"error allocating buffer");
	for(i=0; i<nr ; i++){
		hold[i]=(KERNEL *) malloc(nc * sizeof(KERNEL));
		if(hold[i] == NULL)
			error(1,"error allocating buffers2");
		for(j=0; j<nc; j++)
			hold[i][j]=(KERNEL) 0.;
	}

	half=(n/2);

	for( row=half ; row < nr-half ; row++){
		for(col=half ; col< nc-half ; col++){
			acc=(&hold[row][col]);
			for(i=0; i<n; i++){
				pd=(&d[row-half+i][col-half]);
				for(j=0; j<n; j++,pd++){
					(*acc)+=f[i][j]*(*pd);
				}
			}
		}
	}

	for(i=0; i<nr; i++){
		for(j=0; j<nc; j++){
			d[i][j]=(IMAGE_DATA) (hold[i][j]/SCALE);
		}
		free(hold[i]);
	}
}