/*
  "fts2gif.c" is my veeery simple program that reads the FITS file,
  computes IRAF's zscale params and pipes raw gray data to "raw2gif"
  utility which creates GIF image. There is also a temporary
  file "%%fts2gif.PIDNO" created in /tmp, used by "raw2gif" and then
  deleted.
  
  Only 2-dimensional images are supported, no extensions, no binary
  tables, just simple images.
  
  Usage:   fts2gif  [-ngray nlevels]  [fts_file [gif_file]]
  
  If argument(s) is omitted, uses stdin/stdout.
  The only option, "-ngray" allows to change default number of levels of
  gray (32).

  "fts2gif" is pretty simple. Also, it can be a memory hog for big images
  as it reads the whole image into double array.

  Feel free to modify, enhance etc. the code.
  If you add/change something important, let me know.

  Michal Szymanski (msz@astrouw.edu.pl)

*/

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <signal.h>

#define MAX(A,B) ((A)>(B)?(A):(B))
#define MIN(A,B) ((A)<(B)?(A):(B))

#define MAX_IMAGE_ROWS      8192
#define ZSPOINTS 613

static double *pix[MAX_IMAGE_ROWS], *pixbuf=NULL;
static int pix_signed=1;
static int image_w, image_h;
static double bscale, b_zero;
static double contrast=0.5;
static FILE *tmpf=NULL;
static char tmpfn[80];

double pixel_value(int, int);

void usage()
{
  fprintf(stderr, "usage: fts2gif [-ngray nlevels] [fts_file [gif_file]]\n");
  exit(1);
}

void clean(int sig)
{
  if ( tmpf  ) fclose(tmpf);
  if ( tmpfn[0] ) unlink(tmpfn);
  exit(sig);
}

int needswab()   /* returns 1 on Intel architecture machines, 0 on Suns */
{
  union { short s; char c[2]; } u;
  u.s = 1;
  return (int)(u.c[0]);
}

int
bell()
{
  fprintf(stderr,"\007");
  return 0;
}

void
error2(char *s1, char *s2)
{
  fprintf(stderr,"%s%s\n",s1,s2);
  if ( tmpf ) clean(0);
  exit(1);
}

void swab2(char *adr, int n)
{
  int i;
  char c;
  for(i=0;i<n;i++) {
    c = *adr;
    *adr = *(adr+1);
    *(adr+1) = c;
    adr += 2;
  }
}
void swab4(char *adr, int n)
     /* when called, adr points to real/long array of */
     /*  n elements. The array CAN be misaligned !     */
{
  int i;
  char c, *cp;
  cp=adr;
  for(i=0;i<n;i++) {
    c=*cp;
    *cp=*(cp+3);
    *(cp+3)=c;
    c=*(++cp);
    *cp=*(cp+1);
    *(++cp)=c;
    cp+=2;
  }
}

void swab8(char *adr, int n)
{
  char c;
  int i,j,k=0;
  for(i=0;i<n;i++) {
    for(j=0;j<4;j++) {
      c=adr[k+j];
      adr[k+j]=adr[k+7-j];
      adr[k+7-j]=c;
    }
    adr+=8;
  }
}

void dsort(int n, double *ra)
{
  int l,ir,i,j;
  double rra;

  if ( n <= 1 ) return;
  l=n/2;
  ir=n-1;
  while ( 10 ) {
    if ( l > 0 ) {
      l--;
      rra=ra[l];
    } else {
      rra=ra[ir];
      ra[ir]=ra[0];
      ir--;
      if ( ir == 0 ) {
	ra[0]=rra;
	return;
      }
    }
    i=l;
    j=l+l+1;
    while ( j <= ir ) {
      if ( j < ir ) 
	if ( ra[j] < ra[j+1] ) j++;
      if ( rra < ra[j] ) {
	ra[i]=ra[j];
	i=j;
	j+=j+1;
      } else {
	j=ir+1;
      }
    }
    ra[i]=rra;
  }
}

double *readfits(FILE *inf, char *infn, int *ncols, int *nrows, 
		 char *object, double *bscale, double *bzero,
		 double **pixbuf)
{
  int nx, ny, bitpix, i, j, npix;
  int eoh=0, size, SWAB=needswab();
  long ofs=0L;
  char *cp, buf[2880];
  unsigned char *crow;
  short *hrow;
  float *frow;
  double *drow;
  int *lrow;
  double *shbuf;
  
  if ( object ) *object = 0;
  *bscale = 1;
  *bzero = 0;
  while ( ! eoh ) {
    if ( fread(buf,2880,1,inf) != 1 ) {
      fprintf(stderr,"Error reading header from %s\nExiting...\n",infn);
      fclose(inf);
      return NULL;
    }
    for(i=0;i<2880;i+=80) {
      if ( strncmp(buf+i,"BITPIX  ",8) == 0 ) {
	sscanf(buf+i+9," %d",&bitpix);
      }
      else if ( strncmp(buf+i,"NAXIS1  ",8) == 0 ) {
	sscanf(buf+i+9," %d",&nx);
      }
      else if ( strncmp(buf+i,"NAXIS2  ",8) == 0 ) {
	sscanf(buf+i+9," %d",&ny);
      }
      else if ( strncmp(buf+i,"BSCALE  ",8) == 0 ) {
	sscanf(buf+i+9," %lf", bscale);
      }
      else if ( strncmp(buf+i,"BZERO   ",8) == 0 ) {
	sscanf(buf+i+9," %lf", bzero);
      }
      else if ( object && strncmp(buf+i,"OBJECT  ",8) == 0 ) {
	sscanf(buf+i+9," ' %[^']",object);
	cp = object + strlen(object) - 1;
	while ( cp >= object && *cp == ' ' ) *cp-- = 0;
      }
      else if ( strncmp(buf+i,"END     ",8) == 0 ) {
	eoh=1;
	break;
      }
    }
    ofs+=2880L;
  }
  npix = nx*ny;
  if (*pixbuf == NULL) {
    if ( *pixbuf ) free(*pixbuf);
    *pixbuf = (double *) malloc(npix*sizeof(double));
  }
  if ( *pixbuf == NULL ) {
    fprintf(stderr,"Error allocating memory for image: %s\n",infn);
    fclose(inf);
    return NULL;
  }
  fseek(inf,ofs,0);
  /*
  **  now read the pixel data
  */
  shbuf = (*pixbuf);
  switch ( bitpix ) {
  case 8:  /* byte data */
    crow = (unsigned char *) malloc(nx);
    for (i=0; i<ny; i++) {
      if ( fread(crow,1,nx,inf) != nx ) {
	fprintf(stderr,"Error reading pixel data from %s\n",infn);
	fclose(inf);
	free(*pixbuf);
	free(crow);
	*pixbuf = NULL;
	return NULL;
      }
      for (j=0; j<nx; j++) {
	*shbuf = (*bzero)+(*bscale)*crow[j];
	++shbuf;
      }
    }
    pix_signed = 1;
    free(crow);
    break;
  case -16:    /* -16 is probably illegal but some crazy guys use it ;-) */
  case 16:
    hrow = (short *) malloc(2*nx);
    for (i=0; i<ny; i++) {
      if ( fread(hrow,2,nx,inf) != nx ) {
	fprintf(stderr,"Error reading pixel data from %s\n",infn);
	fclose(inf);
	free(*pixbuf);
	free(hrow);
	*pixbuf = NULL;
	return NULL;
      }
      if ( SWAB ) swab2((char *)hrow, nx);
      for (j=0; j<nx; j++) {
	*shbuf = (*bzero)+(*bscale)*hrow[j];
	++shbuf;
      }
    }
    free(hrow);
    break;
  case 32:  /* 32-bit integer data */
    lrow = (int *) malloc(4*nx);
    for (i=0; i<ny; i++) {
      if ( fread(lrow,4,nx,inf) != nx ) {
	fprintf(stderr,"Error reading pixel data from %s\n",infn);
	fclose(inf);
	free(*pixbuf);
	free(lrow);
	*pixbuf = NULL;
	return NULL;
      }
      if ( SWAB ) swab4((char *)lrow, nx);
      for (j=0; j<nx; j++) {
	*shbuf = (*bzero)+(*bscale)*lrow[j];
	++shbuf;
      }
    }
    free(lrow);
    break;
  case -32:  /* (single) real data */
    frow = (float *) malloc(4*nx);
    for (i=0; i<ny; i++) {
      if ( fread(frow,4,nx,inf) != nx ) {
	fprintf(stderr,"Error reading pixel data from %s\n",infn);
	fclose(inf);
	free(*pixbuf);
	free(frow);
	*pixbuf = NULL;
	return NULL;
      }
      if ( SWAB ) swab4((char *)frow, nx);
      for (j=0; j<nx; j++) {
	*shbuf = frow[j];
	++shbuf;
      }
    }
    free(frow);
    break;
  case -64:  /* (double) real data */
    drow = (double *) malloc(8*nx);
    for (i=0; i<ny; i++) {
      if ( fread(drow,8,nx,inf) != nx ) {
	fprintf(stderr,"Error reading pixel data from %s\n",infn);
	fclose(inf);
	free(*pixbuf);
	free(drow);
	*pixbuf = NULL;
	return NULL;
      }
      if ( SWAB ) swab8((char *)drow, nx);
      for (j=0; j<nx; j++) {
	*shbuf = drow[j];
	++shbuf;
      }
    }
    free(drow);
    break;
  default:
    fprintf(stderr, "Unsupported BITPIX parameter: %d\n", bitpix);
    fclose(inf);
    free(*pixbuf);
    *pixbuf = NULL;
    return NULL;
  }
  *ncols = nx;
  *nrows = ny;
  fclose(inf);
  return *pixbuf;
}

double fit(double *points, int zspoints, int n1, int n2, double *a, double *b)
{
  double x, y;
  int i, sh=zspoints/2;
  double N, Sx=0, Sy=0, Sxx=0, Sxy=0;
  double Syy=0;
  double D, xm, ym, r;

  for(i=n1; i<n2; i++) {
    x = i - sh;
    y = points[i];
    Sx += x;
    Sy += y;
    Sxy += x*y;
    Sxx += x*x;
    Syy += y*y;
  }
  N = n2 - n1;
  D = (N*Sxx-Sx*Sx);
  *a = (Sxx*Sy-Sx*Sxy)/D;
  *b = (N*Sxy-Sx*Sy)/D;
  xm = Sx / N;
  ym = Sy / N;
  r = Sxy - xm*Sy - ym*Sx + N*xm*ym;
  r /= sqrt(Sxx - 2*xm*Sx + N*xm*xm);
  r /= sqrt(Syy - 2*ym*Sy + N*ym*ym);
  return r;
}

int zscale(double *buf, int npix, double *z1, double *z2)
{
  int i, k, j, step, midpoint, zspoints;
  double a, slope, r, pix, points[ZSPOINTS];

  zspoints = MIN(npix, ZSPOINTS);
  step = npix/zspoints;
  for (i=j=k=0; i<zspoints; i++,j+=step) {
    points[k++] = buf[j];
  }
  if ( k < zspoints ) {
    for (i=k,j=0; i<zspoints; i++,j+=(npix/(zspoints-k))) {
      points[i] = buf[j];
    }
  }
  dsort(zspoints, points);
  for (i=20; i<zspoints/4; i+=10) {
    r = fit(points, zspoints, i, zspoints-i, &a, &slope);
    if ( r >= 0.9 ) break;
  }
  midpoint = zspoints/2;
  *z1 = points[midpoint] + (slope / contrast) * (1 - midpoint);
  *z2 = points[midpoint] + (slope / contrast) * (zspoints - midpoint);
/*
  fprintf(stderr,"z1=%d z2=%d i=%d r=%5.2f slope=%5.2f medpix=%d\n",
          *z1,*z2,i,r,slope,points[midpoint]);
*/
  return 0;
}

int main(int argc, char **argv)
{
  FILE *f, *inf;
  static char object[20], label[48], fname[128], cmd[80];
  char *gifname=NULL;
  int i, j;
  unsigned char pixel;
  double value, *ret;
  int ng1, ngray=32;
  double image_z1, image_z2, z_diff;
  static unsigned char g, gray[256];

  while (argc > 1 && argv[1][0] == '-') {
    if (argc > 2 && strcmp(argv[1],"-ngray")==0) {
      ngray = atoi(argv[2]);
      if (ngray<2) ngray=2;
      if (ngray>256) ngray=256;
      --argc;
      ++argv;
    }
    else {
      usage();
    }
    --argc;
    ++argv;
  }
  if (argc < 2) {
    strcpy(fname, "stdin");
    inf=stdin;
    setbuf(inf, NULL);
  }
  else {
    strcpy(fname, argv[1]);
    inf = fopen(fname, "r");
    if (inf == NULL) error2("error opening FITS file: ", fname);
  }
  if (argc > 2) gifname=argv[2];
  ret = readfits(inf, fname, &image_w, &image_h, object, &bscale, &b_zero,
		 &pixbuf);
  if (ret == NULL) error2("error reading FITS file: ", fname);
  signal(SIGINT, clean);
  signal(SIGTERM, clean);
  sprintf(tmpfn, "/tmp/%%fts2gif.%d", getpid());
  tmpf = fopen(tmpfn, "w");
  if ( tmpf == NULL ) error2("error opening tmp file: ", tmpfn);
  for (i=0; i<ngray; i++) {
    g = gray[i] = (i*255)/(ngray-1);
    fprintf(tmpf, "%d %d %d %d\n", i+1, g, g, g);
  }
  fclose(tmpf);
  tmpf = NULL;
  for (i=0; i<image_h; i++)  pix[i] = pixbuf+i*image_w;
  zscale(pixbuf, image_w*image_h, &image_z1, &image_z2);
  ng1 = ngray-1;
  z_diff = image_z2-image_z1;
  if ( gifname )
    sprintf(cmd, "/usr/local/cfitsio/raw2gif -p %s -s %d %d > %s", tmpfn, image_w, image_h,
	    gifname);
  else
    sprintf(cmd, "/usr/local/cfitsio/raw2gif -p %s -s %d %d", tmpfn, image_w, image_h);
  f = popen(cmd, "w");
  if ( f == NULL ) error2("error opening pipe to raw2gif program", "");
  for (i=image_h-1; i>=0; i--) for(j=0; j<image_w; j++) {
    value = pix[i][j];
    if ( value <= image_z1 ) pixel = 0;
    else if ( value >= image_z2 ) pixel = ngray-1;
    else pixel = (value-image_z1)*ng1/z_diff+0.5;
    fwrite(&pixel, 1, 1, f);
  }
  pclose(f);
  clean(0);
  return 0;
}