/*
   program to read BYU SIR image files and write out a byte file

   Written by DGL March 26, 1997

   Written in ANSI C.  Define or undefine SWAP based on machine byte order

*/

#define SWAP 1                /* swap bytes (little endian: dec, pc) */
#undef SWAP                    /* no byte swapping (big endian: HP,SUN) */

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

#define REL_BEGIN 0           /* fseek relative to beginning of file */

#define SIR_HEADER_SIZE 512   /* SIR file header size in bytes */
#define BUFFER 6000           /* space in buffer */ 

void pixtolatlon(int, int, float *, float *);
void latlon2pix(float, float, float *, float *);
void f2ipix(float, float, int *, int *, int, int);

int read_sir_header(FILE *imf, int *nhead, int *ndes, int *nhtype,
		    int *idatatype, int *nsx, int *nsy, 
		    float *xdeg, float *ydeg, 
		    float *ascale, float *bscale, float *a0, float *b0, 
		    int *ioff, int *iscale, int *iyear, 
		    int *isday, int *ismin, int *ieday, int *iemin, 
		    int *iregion, int *itype, int *iopt,
		    int *ispare1, int *ispare2, int *ispare3,
		    float *anodata, float *v_min, float *v_max,
		    char *sensor, char *title, char *type, char *tag,
		    char *crproc, char *crtime, int maxdes, 
		    char *descrip, int *ldes, int maxi, short *iaopt,
		    int *nia);

#ifdef SWAP
void swap(short *, int);
#endif


int main(int argc, char **argv)
{
  FILE *imf, *omf;
  int i, j, k;
  int ix, iy, ix1, iy1;
  float x, y, x1, y1, thelon, thelat, alon, alat, am, amag;
  char outfname[80];
	
  float smin, smax;
  float s, soff;
  float scale, scaleoffset;

/* SIR file header information */

  float xdeg, ydeg, ascale, bscale, a0, b0;
  int iopt;
  float v_min, v_max, anodata;
  int nsx, nsy, ioff, iscale, iyear, isday, ismin, ieday, iemin;
  int iregion, itype, nhead, ndes, nhtype, idatatype, ldes, nia;
  int ipol, ifreqhm, ispare1;
  char title[101], sensor[41], crproc[101], type[139], tag[101], crtime[29];

#define MAXDES 1024
  char descrip[MAXDES+1];
#define MAXI 128
  short iaopt[MAXI];

  float *stval;

  int ierr;
  short *a;
  char *data;
  int nlines, width;

  union 
  {
    short i2[2];
    float f2;
  } un;

  fprintf(stdout,"BYU SIR sir2byte c program\n") ;
  if(argc < 3) {
    fprintf(stdout,"\nusage: %s file out1 <dmin> <dmax>\n\n",argv[0]);
    fprintf(stdout," input parameters:\n");
    fprintf(stdout,"   file  =  input SIR file\n");
    fprintf(stdout,"   out   =  output byte file\n");
    fprintf(stdout,"   dmin  =  min saturation value (optional)\n");
    fprintf(stdout,"   dmax  =  max saturation value (optional)\n"); 
    return(0);
  }

  imf = fopen(argv[1],"r"); 
  if (imf == NULL) {
     fprintf(stdout,"ERROR: cannot open input file: %s\n",argv[1]); 
     exit(-1);
  }

/* get SIR image header information */

  ierr = read_sir_header(imf, &nhead, &ndes, &nhtype,
		    &idatatype, &nsx, &nsy, 
		    &xdeg, &ydeg, &ascale, &bscale, &a0, &b0, 
		    &ioff, &iscale, &iyear, &isday, &ismin, &ieday, &iemin, 
		    &iregion, &itype, &iopt, &ipol, &ifreqhm, &ispare1,
		    &anodata, &v_min, &v_max,
		    sensor, title, type, tag, crproc, crtime, 
		    MAXDES, descrip, &ldes, MAXI, iaopt, &nia);  
  if (ierr < 0) {
     fprintf(stdout,"*** Error reading SIR header from file '%s'\n",argv[1]);
     exit(-1);
  }

/* write out SIR file header information */

  fprintf(stdout,"\nSIR file header: '%s'\n",argv[1]);
  if (nhtype < 20) fprintf(stdout," (Old style header) %d\n",nhtype);
  fprintf(stdout,"  Title:   '%s'\n",title);
  fprintf(stdout,"  Sensor:  '%s'\n",sensor);
  if (nhtype > 16) {
     fprintf(stdout,"  Type:    '%s'\n",type);
     fprintf(stdout,"  Tag:     '%s'\n",tag);
     fprintf(stdout,"  Creator: '%s'\n",crproc);
     fprintf(stdout,"  Created: '%s'\n",crtime);
  }
  fprintf(stdout,"  Size: %d x %d    Total:%d\n",nsx,nsy, nsx*nsy);
  fprintf(stdout,"  Offset: %d  Scale: %d\n",ioff,iscale);
  fprintf(stdout,"  Year: %d  JD range: %d-%d\n",iyear,isday,ieday);
  fprintf(stdout,"  Region Number: %d  Type: %d  Form: %d\n",iregion,itype,iopt);
  if (nhtype > 16) {
     fprintf(stdout,"  Polarization: %d  Frequency: %f MHz\n",ipol,ifreqhm*0.1);
     fprintf(stdout,"  Datatype: %d  Headers: %d  Ver:%d\n",idatatype,nhead,nhtype);
     fprintf(stdout,"  Nodata: %f\n   Vmin: %f  Vmax: %f\n",anodata,v_min,v_max);
     if (ldes > 0)
	fprintf(stdout,"  Description: (%d) '%s'\n",ldes,descrip);
     if (nia > 0) {
	fprintf(stdout,"  Extra Ints: %d\n",nia);
	for (i=0; i < nia; i++)
	   fprintf(stdout,"     %d %d",i,iaopt[i]);
     }	
  }

  switch(iopt) {
   case -1:
     fprintf(stdout,"  Rectangular image-only: \n");
     fprintf(stdout,"  Xspan,  Yspan:  %f , %f\n",xdeg,ydeg);
     fprintf(stdout,"  Xscale, Yscale: %f , %f\n",ascale,bscale);
     fprintf(stdout,"  Xorg,   Yorg:   %f , %f\n",a0,b0);
     break;

   case 0:
     fprintf(stdout,"  Rectangular Lat/Long form: \n");
     fprintf(stdout,"  Size (deg):     %f , %f\n",xdeg,ydeg);
     fprintf(stdout,"  Lon, Lat scale: %f , %f (pix/deg)\n",ascale,bscale);
     fprintf(stdout,"  Offsets:        %f , %f\n",a0,b0);
     break;

   case 2:
     fprintf(stdout,"  Lambert form: (local radius used)\n");
   case 1:
     if (iopt==1) fprintf(stdout,"  Lambert form:\n");
     fprintf(stdout,"  Center point:      %f , %f\n",xdeg,ydeg);
     fprintf(stdout,"  Lon, Lat scale:    %f , %f (km/pix)\n",1./ascale,1./bscale);
     fprintf(stdout,"  Lower-Left Corner: %f , %f\n",a0,b0);
     break;

   case 5:
     fprintf(stdout,"  Polar sterographic form: \n");
     fprintf(stdout,"  Center Lon,Lat:    %f , %f\n",xdeg,ydeg);
     fprintf(stdout,"  X,Y scales:        %f , %f (km/pix)\n",ascale,bscale);
     fprintf(stdout,"  Lower-Left Corner: %f , %f\n",a0,b0);
     break;

   case 11:
   case 12:
     fprintf(stdout,"  EASE polar azimuthal form: \n");
     fprintf(stdout,"  Map center (col,row): %f , %f\n",xdeg,ydeg);
     fprintf(stdout,"  A,B scales:           %f , %f\n",ascale,bscale);
     fprintf(stdout,"  Map origin (col,row): %f , %f\n",a0,b0);
     break;

   case 13:
     fprintf(stdout,"  EASE cylindrical form: \n");
     fprintf(stdout,"  Map center (col,row): %f , %f\n",xdeg,ydeg);
     fprintf(stdout,"  A,B scales:           %f , %f\n",ascale,bscale);
     fprintf(stdout,"  Map origin (col,row): %f , %f\n",a0,b0);
     break;

   default:
     fprintf(stdout,"  Unrecognized SIR file option: \n");
     fprintf(stdout,"  Xspan,  Yspan:  %f , %f\n",xdeg,ydeg);
     fprintf(stdout,"  Xscale, Yscale: %f , %f\n",ascale,bscale);
     fprintf(stdout,"  Xorg,   Yorg:   %f , %f\n",a0,b0);
     break;
  }

  fprintf(stdout,"\n  Image Min, Max: %f , %f\n\n",v_min,v_max);
  rewind(imf);

  smin = v_min;
  smax = v_max;

  if (argc > 2) sscanf(argv[3],"%f",&smin); 
  if (argc > 3) sscanf(argv[4],"%f",&smax); 
  if (smax - smin == 0.0) smax = smin + 1.0;
  fprintf(stdout,"\nByte Min, Max: %f , %f\n\n",smin,smax);

  width = nsx;
  nlines = nsy;

  omf = fopen(argv[2],"w"); 
  if (omf == NULL) {
     fprintf(stderr,"ERROR: cannot open output byte file: '%s'\n",argv[2]); 
     exit(-1);
  }
  fprintf(stdout,"Writing output byte file '%s'\n", argv[2]);

  a =(short *) malloc(sizeof(short)*width*2);
  data = (char *) malloc(sizeof(char)*width*nlines);
  stval = (float *) malloc(sizeof(float)*width*nlines);

  if (a == NULL || data == NULL || stval == NULL) {
     fprintf(stderr,"*** ERROR: memory allocation failure...\n");
     exit(-1);
  }

  s = 1.0 / ((float) iscale);
  soff = 32767.0/ ((float) iscale) + ((float) ioff);
  if (idatatype == 1) soff = 128.0/ ((float) iscale) + ((float) ioff);

  scale = (smax-smin);
  if (scale > 0.) scale = 255./ scale;
  scaleoffset = smin;

/* read file, converting two byte integers to floats then to bytes and
   writing bytes to a separate file.  For later use, copy floating
   point valuses to an array */

  fseek(imf,SIR_HEADER_SIZE*nhead,REL_BEGIN);

  k = 0;
  for (i = 0; i < nlines; i++) {

    switch(idatatype) {
    case 1:     /* bytes */
      fread((char *)a, sizeof(char), width, imf) ;
      for (j = 0; j < width; j++){
	amag = ((float) a[j]) * s + soff; /* pixel value */
	*(stval+k+j) = amag;              /* floating point array */
      }
      break;
      
    case 4:     /* floats */
#ifdef SWAP
      fread((char *)a, sizeof(short), 2*width, imf) ;
      swap(a,width);
      for (j = 0; j < width; j++){
	un.i2[1]=a[j*2];
	un.i2[0]=a[j*2+1];
	amag = un.f2;
	*(stval+k+j) = amag;              /* floating point array */
      }
#else
      fread((char *) (stval+k), sizeof(float), width, imf);
#endif
      break;
      
    default:    /* 0 or 2: read two byte integers */
      fread((char *)a, sizeof(short), width, imf) ;
#ifdef SWAP	
      swap(a,width);
#endif
      for (j = 0; j < width; j++){
	amag = ((float) a[j]) * s + soff; /* pixel value */
	*(stval+k+j) = amag;              /* floating point array */
      }
    }
    
    for (j = 0; j < width; j++){      
      amag = *(stval+k);

      /* scale floating point to byte values */
      am = scale * (amag - scaleoffset);
      if(am > 255.) am = 255.;		/* check overflow */
      if(am < 0.) am = 0.;		/* check underflow */
      *(data+k) = (char)((int)(am));    /* byte array */

      k++;
    }
  }
  fprintf(stdout,"data file read complete\n");
  
  fwrite(data, sizeof(char), width*nlines, omf);
  fclose(imf);
  fclose(omf);
  fprintf(stdout,"SIR input file successfully read and\n");
  fprintf(stdout,"BYTE output file successfully written\n");

  return(0);
}

	    
int read_sir_header(FILE *imf, int *nhead, int *ndes, int *nhtype,
		    int *idatatype, int *nsx, int *nsy, 
		    float *xdeg, float *ydeg, 
		    float *ascale, float *bscale, float *a0, float *b0, 
		    int *ioff, int *iscale, int *iyear, 
		    int *isday, int *ismin, int *ieday, int *iemin, 
		    int *iregion, int *itype, int *iopt,
		    int *ipol, int *ifreqhm, int *ispare1,
		    float *anodata, float *v_min, float *v_max,
		    char *sensor, char *title, char *type, char *tag,
		    char *crproc, char *crtime, int maxdes, 
		    char *descrip, int *ldes, int maxi, short *iaopt,
		    int *nia)
{
   /* v2 read header of BYU SIR file --- written by D.Long March 1997 */

/*  Header consists of a variable number of 512 byte blocks.  The
    block contains scaling information and strings:
    note that character strings may not be null terminated in file

    The first block consists of 256 short integers.  Indexing them
    from 1..256 in an array temp, they contain:

	temp(1)	= nsx			! pixels in x direction
	temp(2) = nsy			! pixels in y direction
        temp(3) <= xdeg			! span of deg in x
	temp(4) <= ydeg			! span of deg in y
	temp(5) = nhtype                ! header type (old<15)
	temp(6) <= ascale               ! x scaling
	temp(7) <= bscale               ! y scaling
	temp(8) <= a0                   ! x (lon) origin
         note: longitudes should be in the range -180 to 180
	temp(9) <= b0                   ! y (lat) origin
	temp(10) = ioff			! offset to be added to scaled val
	temp(11) = iscale		! scale factor ival=(val-ioff)/iscale
	temp(12) = iyear		! year for data used
	temp(13) = isday		! starting JD
	temp(14) = ismin		! time of day for first data (in min)
	temp(15) = ieday		! ending JD
	temp(16) = iemin		! time of day for last data (in min)
	temp(17) = iopt			! projection type
					!  -1 = no projection, image only
					!   0 = rectalinear lat/lon
					!   1 = lambert equal area
					!   2 = lambert equal area (local rad)
					!   5 = polar stereographic
					!  11 = EASE north equal area grid
					!  12 = EASE south equal area grid
					!  13 = EASE cylindrical grid
	temp(18) = iregion		! region id code
	temp(19) = itype		! image type code
                                        ! standard values: 0=unknown or n/a
                                        ! 1 = scatterometer A (dB)
                                        ! 2 = scatterometer B (dB/deg)
                                        ! 3 = radiometer Tb (K)
                                        ! 9 = topography (m)
	temp(20)-temp(39) 40 chars of sensor
        temp(40) = 0 
        temp(41) = nhead                ! number of 512 byte header blocks
        temp(42) = ndes                 ! number of 512 byte blocks description
        temp(43) = ldes                 ! number of bytes of description
        temp(44) = nia                  ! number of optional integers
        temp(45) = ipol                 ! polarization (0=n/a,1=H,2=V)
        temp(46) = ifreqhm              ! frequency in 100's of MHz
        temp(47) = ispare1              ! spare
        temp(48) = idatatype            ! data type code 0,2=i*2,1=i*1,4=f

       the value of idata type determines how data is stored and how
       anodata, vmin, and vmax are stored.

        if idatatype = 1 data is stored as bytes and minv=128
        if idatatype = 2 data is stored as 2 byte integers and minv=32766
        if idatatype = 4 data is stored as IEEE floating point

        if idatatype = 1,2 anodata,vmin,vmax are stored as 2 byte integers
         in temp(49)..temp(51)  minv, ioff and iscal used to convert
         integers or bytes into floating point values
         nodata, vmin, and vmax must be representable with ioff and iscale
            temp(*) = (value-ioff)*iscale-minv
            value = float(temp(*)+minv)/float(iscale)+ioff
        idatatype=2 is considered the SIR standard format

        if idatatype = f anodata,vmin,vmax are stored as floating points
         in temp(52)..temp(57) and minv, ioff and iscale are ignored here
         and when reading the file.
         floating point numbers are NOT standard across platforms and
         is not recommended

        temp(49) <= anodata            ! value representing no data
        temp(50) <= vmin               ! minimum useful value from creator prg
        temp(51) <= vmax               ! maximum useful value from creator prg
        temp(52,53) = anodata          ! IEEE floating value of no data
        temp(54,55) = vmin             ! IEEE floating minimum useful value
        temp(56,57) = vmax             ! IEEE floating maximum useful value

	temp(58)-temp(126) 138 chars of type
        temp(127) = 0 
        temp(128) = 0 
	temp(129)-temp(168) 80 chars of title
        temp(169) = 0 
	temp(170)-temp(189) 40 chars of tag
        temp(190) = 0 
	temp(191)-temp(240) 100 chars of crproc
        temp(241) = 0 
	temp(242)-temp(255) 28 chars of crtime
        temp(256) = 0 

     optional header blocks:

	ndes header blocks of 512 bytes: chars of description
	nhead-ndes-1 header blocks of 512 bytes: values of iaopt
        by convention, first value iaopt is a code telling how to interpret
        the rest of the array if nia>0

        remainder of file is image data in a multiple of 512 byte blocks
        two byte integer and byte scaling is
          intval = (fvalue-ioff)*iscale-minv
          fvalue = float(intval+minv)/float(iscale)+ioff

*/

   short in, temp[256];
   int nch,i;
   float f,soff;

   union 
     {
       short i2[2];
       float f2;
   } un;

   rewind(imf);

/* read first header block */

   if (fread(temp, sizeof(short), 256, imf) == 0) return(-1);
#ifdef SWAP
   swap(temp, 256);
#endif

/* decode nsx */

   *nsx = temp[0];
   *nsy = temp[1];
   *xdeg = 0.01 * temp[2];
   *ydeg = 0.01 * temp[3];
   *nhtype = temp[4];
   if (*nhtype < 20) *nhtype=1;
   *ascale = 0.001 * temp[5];
   *bscale = 0.001 * temp[6];
   *a0 = 0.01 * temp[7];
   *b0 = 0.01 * temp[8];
   *ioff = temp[9];
   *iscale = temp[10];
   if (*iscale == 0) *iscale=1;
   *iyear = temp[11];
   *isday = temp[12];
   *ismin = temp[13];
   *ieday = temp[14];
   *iemin = temp[15];
   *iopt  = temp[16];
   *iregion = temp[17];
   *itype = temp[18];

   switch (*iopt){
    case -1:
      *xdeg = *xdeg * 10.0;
      *ydeg = *ydeg * 10.0;
      break;
    case 0:
      *xdeg = *xdeg + 100.0;
      break;
    case 1:
    case 2:
      *ascale = 1./ *ascale;
      *bscale = 1./ *bscale;
      *a0 = *a0 * 100.0;
      *b0 = *b0 * 100.0;
      break;
    case 5:
      *ascale = *ascale * 10.0;
      *bscale = *bscale * 10.0;
      *a0 = *a0 *100.0;
      *b0 = *b0 *100.0;
      *xdeg = *xdeg + 100.0;
      break;
    case 11:
    case 12:
    case 13:
      *ascale = (float) (((double) *ascale) * (double) 2.0 * 6371.228/25.067525);
      *bscale = (float) (((double) *bscale) * (double) 2.0 * 25.067525);
      *xdeg = *xdeg * 10.0;
      *ydeg = *ydeg * 10.0;
      *a0 = *a0 * 10.0;
      *b0 = *b0 * 10.0;
      break;
    default:
      fprintf(stderr,"\n *** Unrecognized SIR option in read_sir_header ***\n");
   }

   *nhead = temp[40];
   if (*nhead == 0) *nhead=1;
   *ndes = temp[41];
   *ldes = temp[42];
   *nia = temp[43];
   *ipol = temp[44];
   *ifreqhm = temp[45];
   *ispare1 = temp[46];
   *idatatype = temp[47];
   if (*idatatype == 0) *idatatype = 2;

   if (*iscale == 0) *iscale=1;
   soff = 32767.0/(float) *iscale;
   if (*idatatype == 1) soff = 128.0/(float) *iscale;

   *anodata = (float) temp[48]/(float) *iscale + *ioff + soff;
   *v_min = (float) temp[49]/(float) *iscale + *ioff + soff;
   *v_max = (float) temp[50]/(float) *iscale + *ioff + soff;

   if (*idatatype == 4) {
#ifdef SWAP
     un.i2[1]=temp[51];
     un.i2[0]=temp[52];
     *anodata = un.f2;
     un.i2[1]=temp[53];
     un.i2[0]=temp[54];
     *v_min = un.f2;
     un.i2[1]=temp[55];
     un.i2[0]=temp[56];
     *v_max = un.f2;
#else
     un.i2[0]=temp[51];
     un.i2[1]=temp[52];
     *anodata = un.f2;
     un.i2[0]=temp[53];
     un.i2[1]=temp[54];
     *v_min = un.f2;
     un.i2[0]=temp[55];
     un.i2[1]=temp[56];
     *v_max = un.f2;
#endif
   }
   
   for (in = 0; in < 20; in++) {
      sensor[2*in] = temp[19+in] % 256;
      sensor[2*in+1] = temp[19+in]/256; 
   }
   sensor[40] = '\0';

   for (in = 0; in < 69; in++) {
      type[2*in] = temp[57+in] % 256;
      type[2*in+1] = temp[57+in]/256;
   }
   type[138] = '\0';

   for (in = 0; in < 40; in++) {
      title[2*in] = temp[128+in] % 256;
      title[2*in+1] = temp[128+in]/256;
   }
   title[80] = '\0';

   for (in = 0; in < 20; in++) {
      tag[2*in] = temp[169+in] % 256;
      tag[2*in+1] = temp[169+in]/256;
   }
   tag[40] = '\0';

   for (in = 0; in < 50; in++) {
      crproc[2*in] = temp[190+in] % 256;
      crproc[2*in+1] = temp[190+in]/256;
   }
   crproc[100] = '\0';

   for (in = 0; in < 14; in++) {
      crtime[2*in] = temp[241+in] % 256;
      crtime[2*in+1] = temp[241+in]/256;
   }
   crtime[28] = '\0';

   if (*nhtype == 1) {  /* old style header */
     *nhead=1;
     *ndes=0;
     *ldes=0;
     *nia=0;
     type[0]='\0';
     tag[0]='\0';
     crproc[0]='\0';
     crtime[0]='\0';
   }

   if (*nhead > 1) {   /* read additional header blocks */
      if (*ndes > 0) {
	 fseek(imf, 512, REL_BEGIN);
	 nch=(maxdes < *ldes ? maxdes : *ldes);
	 nch=0;
	 for (i = 0; i < *ndes; i++) {
	   if (fread(temp, sizeof(short), 256, imf) == 0) return(-1);
#ifdef SWAP
	   swap(temp, 256);
#endif
	   for (in = 0; in < 256; in++) {
	     if (nch < maxdes) descrip[nch] = temp[in] % 256;
	     nch++;
	     if (nch < maxdes) descrip[nch] = temp[in]/256;
	     nch++;
	   }
	 }
      }
      if (maxdes > 0) descrip[maxdes-1] = '\0';

      if (*nhead-*ndes-1 > 0) {
	 fseek(imf, 512*(*ndes+1), REL_BEGIN);
	 nch=(maxi < *nia ? maxi : *nia);
	 if (fread(iaopt, sizeof(short), nch, imf) == 0) return(-1);
#ifdef SWAP
	 swap(iaopt, nch);
#endif
      }
   }

   return(0);
}

#ifdef SWAP

void swap(short *i, int n)
{
   char *s, t;
   int j;
   
   for (j = 0; j < n; j++) {
     s = (char *) (i+j);
     t = *s;
     *s = *(s+1);
     *(s+1) = t;
   }
   return;
}

#endif