/*
 * Module: windowing
 *
 * Routines for reading/writing and convolving.
 *
 *
 * Author:
 *   Stuart Inglis (singlis@internz.co.nz)
 *   (c) 1998 
 *
 */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "windowing.h"
#include "utils.h"
#include <assert.h>



double window_rectangular_1D(int i, int N)
{
  i=i; /* to quieten compilers */
  N=N;
  return 1.0;
}

double window_circular_1D(int i, int N)
{
    if((i>=0) && (i<=N))
	return 1.0;
    else
	return 0.0;
}

double window_bartlett_1D(int i, int N)
{
  if((i>=0) && (i<=N)){
    if(N==0)
      return 1.0;
    else
      return 1.0-i/(double)N;
  }
  else 
    return 0.0;
}

double window_hanning_1D(int i, int N)
{
  if((i>=0) && (i<=N)){
    if(N==0)
      return 1.0;
    else
      return (1.0-cos(PI*(1.0+i/(double)(N))))/2.0;
  }
  else
    return 0.0;
}


double window_hamming_1D(int i, int N)
{
  if((i>=0) && (i<=N)){
    if(N==0)
      return 1.0;
    else 
      return (0.54-0.46*cos(PI*(1.0+i/(double)(N))));
  }
  else
    return 0.08;
}


double window_blackman_1D(int i, int N)
{
  if((i>=0) && (i<N)){
    if(N==0)
      return 1.0;
    else 
      return (0.42-0.5*cos(PI*(1.0+i/(double)(N)))+0.08*
	    cos(2.0*PI*(1.0+i/(double)(N))));
  }
  else
    return 0.0;
}


/* sd is a percentage of N/2, ie. 0.5 = 50% of N/2 */
double window_gaussian_1D(int i, int N, double sd)
{
    sd=N*sd;
    
    if(N==0)
      return 1.0;
    else 
      return exp(-( (double)i*i)/(2.0*sd*sd));
}







void create_window(
		   double *window,
		   int width,
		   int TYPE
		   )

{
    int r,i;

    assert(window);
    assert(width>0);

    if((width%2)==0)
	fprintf(stderr,"comment: usually recommended to have an 'odd' size window width\n");

    for(r=0;r<width;r++){
	i=abs(r-width/2);
	switch (TYPE) {
	  case RECTANGULAR: window[r]=window_rectangular_1D(i,width/2);
	    break;
	  case BARTLETT:    window[r]=window_bartlett_1D(i,width/2);
	    break;
	  case HANNING:     window[r]=window_hanning_1D(i,width/2);
	    break;
	  case HAMMING:     window[r]=window_hamming_1D(i,width/2);
	    break;
	  case BLACKMAN:    window[r]=window_blackman_1D(i,width/2);
	    break;
	  case GAUSSIAN:    window[r]=window_gaussian_1D(i,width/2,0.4);
	    break;
	default:{
	      fprintf(stderr,"create_window: don't know this window type!\n");
	      exit(1);
	    }
	}
    }

}

void create_window_2D(
		   double **window,
		   int width,
		   int TYPE
		   )

{
    int i,j;
    int r;
    int width_r;

    if((width%2)==0)
	fprintf(stderr,"comment: usually recommended to have an "
		"odd size window width\n");

    /* circumscribe */
    width_r=ROUND(sqrt((width/2*width/2)));

    /* at corners */
    /*    width_r=ROUND(sqrt((width/2*width/2) * 2));*/

    for(j=0;j<width;j++)
	for(i=0;i<width;i++){
	    r=ROUND(sqrt((double)(j-width/2)*(j-width/2)+
			 (double)(i-width/2)*(i-width/2)));
	    switch (TYPE) {
	      case RECTANGULAR: window[j][i]=window_rectangular_1D(r,width_r);
		break;
	      case BARTLETT:    window[j][i]=window_bartlett_1D(r,width_r);
		break;
	      case HANNING:     window[j][i]=window_hanning_1D(r,width_r);
		break;
	      case HAMMING:     window[j][i]=window_hamming_1D(r,width_r);
		break;
	      case BLACKMAN:    window[j][i]=window_blackman_1D(r,width_r);
		break;
	      case GAUSSIAN:    window[j][i]=window_gaussian_1D(r,width_r,0.4);
		break;
	    default:{
		  fprintf(stderr,"create_window: don't know this window type!\n");
		  exit(1);
		}
	    }
	}
}


void write_window(
		  char *filename,
		  double *window,
		  int L,
		  int H
		  )

{
    int i;
    FILE *fp;
    int open=0;

    assert(window);
    assert(H>=L);

    if((!filename) || (strcmp(filename,"-")==0)){
	fp=(stdout);
    } else{
	fp=fopen(filename,"wb");
	if(!fp) {
	    fprintf(stderr,"write_window: can't create file '%s'\n",filename);
	    return ;
	}
	open=1;
    }

    for(i=L;i<=H;i++)
	fprintf(fp,"%d %g\n",i,window[i]);
    
    fflush(fp);
    if(open)
	fclose(fp);
}


void write_window_labelled(
		  char *filename,
		  double *window,
		  int L,
		  int H,
		  double lowrange,
		  double highrange,
                  double bin
		  )

{
    int i;
    double r;
    FILE *fp;
    int open=0;

    assert(window);
    assert(H>=L);

    highrange=highrange; /* to quieten compilers */

    if(bin<=0){
	fprintf(stderr,"write_window_labelled: warning, binsize<=0\n");
    }

    if((!filename) || (strcmp(filename,"-")==0))
	fp=(stdout);
    else{
	fp=fopen(filename,"wb");
	if(!fp) {
	    fprintf(stderr,"write_window: can't create file '%s'\n",filename);
	    return;
	}
	open=1;
    }

    r=lowrange;
    for(i=L ;i<=H; i++,	r+=bin){
/*	
	r=lowrange+(i-L)*(highrange-lowrange)/(H-L);
*/
	if(bin<0.01)
	    fprintf(fp,"%.3f %g\n",r,window[i]);
	else if(bin<0.1)
	    fprintf(fp,"%.2f %g\n",r,window[i]);
	else if(bin<1)
	    fprintf(fp,"%.1f %g\n",r,window[i]);
	else 
	    fprintf(fp,"%g %g\n",r,window[i]);
    }
    
    fflush(fp);
    if(open)
	fclose(fp);
}


void write_window_2D(
		  char *filename,
		  double **window,
		  int Ly,
		  int Hy,
		  int Lx,
		  int Hx
		  )

{
    int i,j;
    FILE *fp;
    int open=0;
    
    if((!filename) || (strcmp(filename,"-")==0))
	fp=(stdout);
    else{
	fp=fopen(filename,"wb");
	if(!fp) {
	    fprintf(stderr,"write_window: can't create file '%s'\n",filename);
	    return;
	}
	open=1;
    }
    
    for(j=Ly;j<=Hy;j++)
	for(i=Lx;i<=Hx;i++)
	    fprintf(fp,"%d\t%d\t%f\n",i,j,window[j][i]);
    
    fflush(fp);
    if(open)
	fclose(fp);
}


void write_window_2D_pgm_P2(
		  char *filename,
		  double **window,
		  int Ly,
		  int Hy,
		  int Lx,
		  int Hx,
		  int maxval,
		  char *commentstring
		  )

{
    double **copy;
    int j,i;
    int open=0;
    FILE *fp;

    CALLOC_2D(copy,Hy-Ly+1,Hx-Lx+1,double);

    for(j=Ly;j<=Hy;j++)
	for(i=Lx;i<=Hx;i++)
	    copy[j-Ly][i-Lx]=window[j][i];

    normalise_window_2D(copy,Ly,Hy,Lx,Hx);


    if((!filename) || (strcmp(filename,"-")==0)){
	fp=(stdout);
    } else{
	fp=fopen(filename,"wb");
	if(!fp) {
	    fprintf(stderr,"write_window_2D_pgm: can't create file '%s'\n",filename);
	    return;
	}
	open=1;
    }

    if(commentstring){
      if(commentstring[strlen(commentstring)-1]=='\n')
	commentstring[strlen(commentstring)-1]=0;
    }
    if(maxval<=255){
      fprintf(fp,"P5\n");
      if(commentstring){
	fprintf(fp,"# %s\n",commentstring);
      }
      fprintf(fp,"%d %d\n",Hx-Lx+1,Hy-Ly+1);
      fprintf(fp,"%d\n",maxval);
    } else {
      fprintf(fp,"P2\n");
      if(commentstring){
	fprintf(fp,"# %s\n",commentstring);
      }
      fprintf(fp,"%d %d\n",Hx-Lx+1,Hy-Ly+1);
      fprintf(fp,"%d\n",maxval);
    }
    
    for(j=Ly;j<=Hy;j++){
	for(i=Lx;i<=Hx;i++){
	    if((int)(copy[j-Ly][i-Lx]*maxval)<0)
		fprintf(stderr,"warning - file '%s' has negative values!!!\n",filename);
	    if(maxval<=255){
	      fputc((int)(copy[j-Ly][i-Lx]*maxval),fp);
	    } else {
	      fprintf(fp,"%d ",(int)(copy[j-Ly][i-Lx]*maxval));
	    }
	}
	if(maxval>255){
	  fprintf(fp,"\n");
	}
    }
    
    fflush(fp);
    if(open)
	fclose(fp);

    FREE_2D(copy,Hy-Ly+1);
    
}




void write_window_2D_pgm_P2_float(
		  char *filename,
		  float **window,
		  int Ly,
		  int Hy,
		  int Lx,
		  int Hx,
		  int maxval
		  )

{
    double **copy;
    int j,i;
    int open=0;
    FILE *fp;

    CALLOC_2D(copy,Hy-Ly+1,Hx-Lx+1,double);

    for(j=Ly;j<=Hy;j++)
	for(i=Lx;i<=Hx;i++)
	    copy[j-Ly][i-Lx]=window[j][i];

    normalise_window_2D(copy,Ly,Hy,Lx,Hx);


    if((!filename) || (strcmp(filename,"-")==0))
	fp=(stdout);
    else{
	fp=fopen(filename,"wb");
	if(!fp) {
	    fprintf(stderr,"write_window_2D_pgm: can't create file '%s'\n",filename);
	    return;
	}
	open=1;
    }

    fprintf(fp,"P2\n");
    fprintf(fp,"%d %d\n",Hx-Lx+1,Hy-Ly+1);
    fprintf(fp,"%d\n",maxval);
    
    for(j=Ly;j<=Hy;j++){
	for(i=Lx;i<=Hx;i++){
	    if((int)(copy[j-Ly][i-Lx]*maxval)<0)
		fprintf(stderr,"warning - file '%s' has negative values!!!\n",filename);
	    fprintf(fp,"%d ",(int)(copy[j-Ly][i-Lx]*maxval));
	}
	fprintf(fp,"\n");
    }
    
    fflush(fp);
    if(open)
	fclose(fp);

    FREE_2D(copy,Hy-Ly+1);
}


void write_window_2D_pgm_P5_float(
		  char *filename,
		  float **window,
		  int Ly,
		  int Hy,
		  int Lx,
		  int Hx,
		  int maxval
		  )

{
    double **copy;
    int j,i;
    int open=0;
    FILE *fp;

    CALLOC_2D(copy,Hy-Ly+1,Hx-Lx+1,double);

    for(j=Ly;j<=Hy;j++)
	for(i=Lx;i<=Hx;i++)
	    copy[j-Ly][i-Lx]=window[j][i];

    normalise_window_2D(copy,Ly,Hy,Lx,Hx);


    if((!filename) || (strcmp(filename,"-")==0))
	fp=(stdout);
    else{
	fp=fopen(filename,"wb");
	if(!fp) {
	    fprintf(stderr,"write_window_2D_pgm: can't create file '%s'\n",filename);
	    return;
	}
	open=1;
    }

    fprintf(fp,"P5\n");
    fprintf(fp,"%d %d\n",Hx-Lx+1,Hy-Ly+1);
    fprintf(fp,"%d\n",maxval);
    
    for(j=Ly;j<=Hy;j++){
	for(i=Lx;i<=Hx;i++){
	  if((int)(copy[j-Ly][i-Lx]*maxval)<0)
	    fprintf(stderr,"warning - file '%s' has negative values!!!\n",filename);
	  fputc((char)(copy[j-Ly][i-Lx]*maxval),fp);
	}
    }
    
    fflush(fp);
    if(open)
      fclose(fp);
    
    FREE_2D(copy,Hy-Ly+1);
}



	
void normalise_window(
		      double *window,
		      int L,
		      int H
		      )

{
    int i;
    double maxval,minval;
    double gap;


    minval=maxval=window[L];

    for(i=L; i<=H; i++){
	if(window[i]>maxval) 
	    maxval=window[i];
	if(window[i]<minval) 
	    minval=window[i];
    }

    gap=maxval-minval;
    if(gap!=0){
	for(i=L; i<=H; i++)
	    window[i]=(window[i]-minval)/gap;
    }
}

void normalise_window_2D(
		      double **window,
		      int Ly,
		      int Hy,
		      int Lx,
		      int Hx
		      )

{
    int i,j;
    double maxval,minval;
    double gap;
    
    minval=maxval=window[Ly][Lx];

    for(j=Ly; j<=Hy;j++)
	for(i=Lx; i<=Hx; i++){
	    if(window[j][i]>maxval) 
		maxval=window[j][i];
	    if(window[j][i]<minval) 
		minval=window[j][i];
	}
    
    gap=maxval-minval;
    
    if(gap!=0){
	for(j=Ly; j<=Hy;j++)
	    for(i=Lx; i<=Hx; i++){
		window[j][i]=(window[j][i]-minval)/gap;
	    }
    }

    for(j=Ly; j<=Hy;j++)
	for(i=Lx; i<=Hx; i++)
	    if(window[j][i]<0){
		fprintf(stderr,"warning - there are negative values in this window\n");return;}

}

void convolve(
	      double *data,
	      int L,
	      int H,
	      double *window,
	      int width,
	      int origin,
	      double *result,
	      int OPTIONS
	      )

{
    int n,k;
    double x,h;
    double sum;
    int b;
    
    for(n=L; n<=H ;n++){
	sum=0;
	for(k=-origin; k<width-origin; k++){
	    if(OPTIONS&CYCLIC){
		b=n+k;
		while(b<L) b+=(H-L+1);
		while(b>H) b-=(H-L+1);
		x=data[b];
	    }
	    else {
		if((n+k<L) || (n+k>H))
		    x=0.0;
		else
		    x=data[n+k];
	    }
	    
	    h=window[k+origin];
	    
	    sum+=x*h;
	}
	result[n]=sum;
    }
}


void convolve_2D(
	      double **data,
	      int Ly,
	      int Hy,
	      int Lx,
	      int Hx,
	      double **window,
	      int width,
	      int originy,
	      int originx,
	      double **result,
	      int OPTIONS
	      )

{
    int i,j;
    int lx,ly;
    double x;
    double sum;
    int bx,by;


    for(j=Ly;j<=Hy;j++)
	for(i=Lx; i<=Hx ;i++){
	    sum=0;

	    for(ly=-originy;ly<width-originy;ly++){
		for(lx=-originx;lx<width-originx;lx++){
		    if(OPTIONS&CYCLIC){
			bx=i+lx;
			while(bx<Lx) bx+=(Hx-Lx+1);
			while(bx>Hx) bx-=(Hx-Lx+1);
			by=j+ly;
			while(by<Ly) by+=(Hy-Ly+1);
			while(by>Hy) by-=(Hy-Ly+1);
			x=data[by][bx]*
			    window[ly+originy][lx+originx];
		    }
		    else {
			if((j+ly<Ly) || (j+ly>Hy) || (i+lx<Lx) || (i+lx>Hx))
			    x=0;
			else
			    x=data[j+ly][i+lx] *
				window[ly+originy][lx+originx];
		    }
		sum+=x;
		}
	    }
	    result[j][i]=sum;
	}
}