#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "fitsio.h"

#define LINEPIX 2100
#define NOPIX  -1e10

void printerror( int status );

int main(int argc, char *argv[])
{
  fitsfile *fp;
  int i,j,status,nfound,anynull,bitpix=-32;
  long fpixel, nelements,naxis=2;
  long naxes[2]={LINEPIX, LINEPIX},nbuffer;
  float nullval;
  static float buffer[LINEPIX],z[LINEPIX][LINEPIX];
  static int m[LINEPIX][LINEPIX],m2[LINEPIX][LINEPIX];
  int mode,is,ie,js,je,n1,n2,keyc=1,ncomb,fmode;
  double ave1,ave2,std1,std2,zis,zie,zjs,zje,zdif;

  if(argc!=3 && argc!=4){
    fprintf(stderr,"Usage : %s input_file output_file [x/y]\n",argv[0]);
    exit(0);
  }
  if(argc==4){
    if(strcmp(argv[3],"x")==0) fmode=1;
    else if(strcmp(argv[3],"y")==0) fmode=2;
    else fmode=0;
  } else fmode=0;

  for(i=0;i<LINEPIX;i++){
    for(j=0;j<LINEPIX;j++) {m[i][j]=1;m2[i][j]=1;}
  }

  fprintf(stderr,"Read mask.fits  ... ");
  status=0;
  if(fits_open_file(&fp,"mask.fits", READONLY, &status))
    printerror(status);
  if(fits_read_keys_lng(fp, "NAXIS", 1, 2, naxes, &nfound, &status))
    printerror(status);
  if(fits_read_key(fp, TINT, "NCOMBINE", &ncomb, NULL, &status)) keyc=0;
  status=0;
  fprintf(stderr,"%ld x %ld format ",naxes[0],naxes[1]);
  if((naxes[0]>LINEPIX)||(naxes[1]>LINEPIX)){
    fprintf(stderr,"\nIllegal image format\n");
    exit(0);
  }
  fpixel = 1;
  nullval = 0;
  nbuffer = naxes[0];
  for(i=0;i<naxes[1];i++){
    if (fits_read_img(fp, TFLOAT, fpixel, nbuffer, &nullval,
                      buffer, &anynull, &status))
      printerror(status);
    for(j=0;j<naxes[0];j++){
      if(buffer[j]<0.5) m[i][j]=0;
    }
    fpixel+=naxes[0];
  }
  if(fits_close_file(fp, &status))
    printerror(status);
  fprintf(stderr,"OK\n");

  if(fits_open_file(&fp,"mask2.fits", READONLY, &status)) mode=0;
  else {
    mode=1;
    if(fits_close_file(fp, &status)) printerror(status);
  }
  if(mode){
    fprintf(stderr,"Read mask2.fits  ... ");
    status=0;
    if(fits_open_file(&fp,"mask2.fits", READONLY, &status))
      printerror(status);
    if(fits_read_keys_lng(fp, "NAXIS", 1, 2, naxes, &nfound, &status))
      printerror(status);
    fprintf(stderr,"%ld x %ld format ",naxes[0],naxes[1]);
    if((naxes[0]>LINEPIX)||(naxes[1]>LINEPIX)){
      fprintf(stderr,"\nIllegal image format\n");
      exit(0);
    }
    fpixel = 1;
    nullval = 0;
    nbuffer = naxes[0];
    for(i=0;i<naxes[1];i++){
      if (fits_read_img(fp, TFLOAT, fpixel, nbuffer, &nullval,
			buffer, &anynull, &status))
	printerror(status);
      for(j=0;j<naxes[0];j++){
	if(buffer[j]<0.5) m2[i][j]=0;
      }
      fpixel+=naxes[0];
    }
    if(fits_close_file(fp, &status))
      printerror(status);
    fprintf(stderr,"OK\n");
  }

  if(mode){n1=0;n2=0;ave1=0.;ave2=0.;std1=0.;std2=0.;}
  fprintf(stderr,"Read %s  ... ",argv[1]);

  status=0;
  if(fits_open_file(&fp,argv[1], READONLY, &status))
    printerror(status);
  if(fits_read_keys_lng(fp, "NAXIS", 1, 2, naxes, &nfound, &status))
    printerror(status);
  fprintf(stderr,"%ld x %ld format ",naxes[0],naxes[1]);
  if((naxes[0]>LINEPIX)||(naxes[1]>LINEPIX)){
    fprintf(stderr,"\nIllegal image format\n");
    exit(0);
  }
  fpixel = 1;
  nullval = 0;
  nbuffer = naxes[0];
  for(i=0;i<naxes[1];i++){
    if (fits_read_img(fp, TFLOAT, fpixel, nbuffer, &nullval,
                      buffer, &anynull, &status))
      printerror(status);
    for(j=0;j<naxes[0];j++) z[i][j]=buffer[j];
    fpixel+=naxes[0];
  }
  if(fits_close_file(fp, &status))
    printerror(status);
  fprintf(stderr,"OK\n");

  if(mode){
    for(i=0;i<naxes[1];i++){
      for(j=0;j<naxes[0];j++){
	if(m[i][j]){
	  n1++;ave1+=z[i][j];std1+=z[i][j]*z[i][j];
	  if(m2[i][j]==0){
	    is=i;
	    while(m2[is][j]==0){
	      is--;
	      if(is==-1) break;
	    }
	    ie=i;
	    while(m2[ie][j]==0){
	      ie++;
	      if(ie==naxes[1]) break;
	    }
	    js=j;
	    while(m2[i][js]==0){
	      js--;
	      if(js==-1) break;
	    }
	    je=j;
	    while(m2[i][je]==0){
	      je++;
	      if(je==naxes[0]) break;
	    }
	    if(is==-1 && ie==naxes[1]){
	      fprintf(stderr,"All the pixel in column %d is bad pixels!\n",j+1);
	      exit(1);
	    }
	    else if(is==-1){zis=z[ie][j];zie=zis;}
	    else if(ie==naxes[1]){zis=z[is][j];zie=zis;}
	    else {zis=z[is][j];zie=z[ie][j];}
	    if(js==-1 && je==naxes[0]){
	      fprintf(stderr,"All the pixel in line %d is bad pixels!\n",i+1);
	      exit(1);
	    }
	    else if(js==-1){zjs=z[i][je];zje=zjs;}
	    else if(je==naxes[0]){zjs=z[i][js];zje=zjs;}
	    else {zjs=z[i][js];zje=z[i][je];}
	    zdif=z[i][j]-((zis*(ie-i)+zie*(i-is))/(ie-is)*(je-js)+(zjs*(je-j)+zje*(j-js))/(je-js)*(ie-is))/(ie-is+je-js);
	    n2++;
	    ave2+=zdif;
	    std2+=zdif*zdif;
	    if(0){
	      fprintf(stderr,"\t%d\t %d\t %d\t %d\t %d\t %d\n",is,i,ie,js,j,je);
	      fprintf(stderr,"\t%g\t %g\t %g\t %g\t %g\t %g\n",zis,z[i][j]-zdif,zie,zjs,z[i][j]-zdif,zje);
	      fprintf(stderr,"\t%g\t %g\t %d\t %g\n",z[i][j],zdif,n2,ave2);getchar();
	    }
	  }
	}
      }
    }

    ave1/=n1;ave2/=n2;std1/=n1;std2/=n2;
    std1=sqrt(std1-ave1*ave1);std2=sqrt(std2-ave2*ave2);
    fprintf(stderr,"Average/Stddev/Npix(all, mask2)=%g/%g/%d, %g/%g/%d\n",ave1,std1,n1,ave2,std2,n2);
    if(fabs(ave2) > std1/sqrt(n2)*20){
      fprintf(stderr,"#### mask2.fits is applied. ###\n");
      //fprintf(stderr,"#### mask2.fits will be applied. OK? (Y/N) ###\n");
      //gc=0;
      //while(gc!=89 && gc!=121 && gc!=78 && gc!=110){
      //gc=getchar();getchar();//fprintf(stderr,"%d\n",gc);
      //}
      //if(gc==89 || gc==121){
	for(i=0;i<LINEPIX;i++){
	  for(j=0;j<LINEPIX;j++) m[i][j]*=m2[i][j];
	}
	//}
    }
  }

  for(i=0;i<LINEPIX;i++){
    for(j=0;j<LINEPIX;j++){
      if(m[i][j]==0){
	is=i;
	while(m[is][j]==0){
	  is--;
	  if(is==-1) break;
	}
	ie=i;
	while(m[ie][j]==0){
	  ie++;
	  if(ie==naxes[1]) break;
	}
	js=j;
	while(m[i][js]==0){
	  js--;
	  if(js==-1) break;
	}
	je=j;
	while(m[i][je]==0){
	  je++;
	  if(je==naxes[0]) break;
	}
	if(is==-1 && ie==naxes[1]){
	  if(fmode!=1){
	    fprintf(stderr,"All the pixel in column %d is bad pixels!\n",j+1);
	    exit(1);
	  }
	}
	else if(is==-1){zis=z[ie][j];zie=zis;}
	else if(ie==naxes[1]){zis=z[is][j];zie=zis;}
	else {zis=z[is][j];zie=z[ie][j];}
	if(zis<NOPIX*0.99) zis=zie;
	else if(zie<NOPIX*0.99) zie=zis;
	if(js==-1 && je==naxes[0]){
	  if(fmode!=2){
	    fprintf(stderr,"All the pixel in line %d is bad pixels!\n",i+1);
	    exit(1);
	  }
	}
	else if(js==-1){zjs=z[i][je];zje=zjs;}
	else if(je==naxes[0]){zjs=z[i][js];zje=zjs;}
	else {zjs=z[i][js];zje=z[i][je];}
	if(zjs<NOPIX*0.99) zjs=zje;
	else if(zje<NOPIX*0.99) zje=zjs;
	if(fmode==1 || zis<NOPIX*0.99) z[i][j]=(zjs*(je-j)+zje*(j-js))/(je-js);
	else if(fmode==2 || zjs<NOPIX*0.99) z[i][j]=(zis*(ie-i)+zie*(i-is))/(ie-is);
	else z[i][j]=((zis*(ie-i)+zie*(i-is))/(ie-is)*(je-js)+(zjs*(je-j)+zje*(j-js))/(je-js)*(ie-is))/(ie-is+je-js);
	if(0){
	  fprintf(stderr,"\t%d\t%d\t %d\t %d\t %d\t %d\t %d\n",fmode,is,i,ie,js,j,je);
	  fprintf(stderr,"\t%g\t %g\t %g\t %g\t %g\t %g\n",zis,z[i][j],zie,zjs,z[i][j],zje);getchar();
	}
      }
    }
  }

  fprintf(stderr,"Write %s  ... ",argv[2]);
  status=0;
  if (fits_create_file(&fp, argv[2], &status))
    printerror(status);
  if (fits_create_img(fp, bitpix, naxis, naxes, &status))
    printerror(status);
  if(keyc){
    if (fits_write_key(fp, TINT, "NCOMBINE", &ncomb, "", &status))
      printerror(status);
  }
  fpixel=1;
  nelements=naxes[0];
  for(i=0;i<naxes[1];i++){
    if (fits_write_img(fp, TFLOAT, fpixel, nelements, z[i], &status))
      printerror(status);
    fpixel+=naxes[0];
  }
  if (fits_close_file(fp, &status))
    printerror(status);
  fprintf(stderr,"OK\n");
  return(0);
}

void printerror( int status)
{
    /*****************************************************/
    /* Print out cfitsio error messages and exit program */
    /*****************************************************/

    char status_str[FLEN_STATUS], errmsg[FLEN_ERRMSG];

    if (status)
      fprintf(stderr, "\n*** Error occurred during program execution ***\n");

    fits_get_errstatus(status, status_str);   /* get the error description */
    fprintf(stderr, "\nstatus = %d: %s\n", status, status_str);

    /* get first message; null if stack is empty */
    if ( fits_read_errmsg(errmsg) )
    {
         fprintf(stderr, "\nError message stack:\n");
         fprintf(stderr, " %s\n", errmsg);

         while ( fits_read_errmsg(errmsg) )  /* get remaining messages */
             fprintf(stderr, " %s\n", errmsg);
    }

    exit( status );       /* terminate the program, returning error status */
}