#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <string.h>
#include <float.h>
#include <time.h>

/* Interface for the PROJ library of geographic projections */
#include <proj_api.h>

#include "bufrlib.h"
#include "bufr_io.h"

#include "cree_png.h"

/*===========================================================================*/

static char * create_filename(char * filename, char * suffix, radar_data_t *od)
{
  char *fmt = "-xx:xx:xx-yy.y";
  char * complete_filename = malloc(strlen(filename) + strlen(fmt) + strlen(suffix) + 1);
  strcpy(complete_filename, filename);
  if (od) sprintf(&complete_filename[strlen(filename)], "-%02d:%02d:%02d-%04.1f", od->meta.hour, od->meta.min, od->meta.sec, od->polar.elangle);
  strcat(complete_filename, suffix);
  return complete_filename;
}

static FILE* initPngFile(char *filename, png_structpp png_ptr_ptr, png_infopp info_ptr_ptr){
  FILE *fp = NULL;
  if (filename != NULL) {
    fp = fopen(filename , "wb");
    if (!fp) return fp;
  }
  else fp = stdout;  
  *png_ptr_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
  if (!png_ptr_ptr){
    fclose(fp);
    return NULL;
  }
  *info_ptr_ptr = png_create_info_struct(*png_ptr_ptr);
  if (!info_ptr_ptr){
    png_destroy_write_struct(png_ptr_ptr, (png_infopp)NULL);
    fclose(fp);
    return NULL;
  }
 
  if (setjmp(png_jmpbuf(*png_ptr_ptr))){  /*set longjump on error */
    png_destroy_write_struct(png_ptr_ptr, info_ptr_ptr);
    fclose(fp);
    return NULL;
  }
  png_init_io(*png_ptr_ptr, fp);
 
  return fp;
}

static void deinit(png_structpp png_ptr_ptr, png_infopp info_ptr_ptr){
  png_write_end(*png_ptr_ptr, *info_ptr_ptr);
  png_destroy_write_struct(png_ptr_ptr, info_ptr_ptr);
}

static unsigned short lookup(double Z)
{
  int i = 0;
  unsigned short res = 255;
  while (i < vraie_taille_palette) {
    if (Z < seuils_france54[i]) {
      res = i;
      return res;
    } else {
      i++;
    }
  }
  return res;
}

static void writeInfo(png_structp png_ptr, png_infop info_ptr, int width, int height){
  int count;
  int size_palette;
  png_byte trans[256];

  png_set_invert_alpha(png_ptr);
  png_set_IHDR(png_ptr, info_ptr, width, height, 8,
      PNG_COLOR_TYPE_PALETTE, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);

  size_palette = sizeof(palette)/sizeof(png_color);
    
  for (count = vraie_taille_palette; count < size_palette; count++) {
    palette[count] = no_value; 
  }
  for (count = 0; count < size_palette; count++) {
    /* KML files are half transparent */
    trans[count] = (kml ? 127 : 0);
  }
  png_set_tRNS(png_ptr, info_ptr, trans, size_palette, NULL);
  png_set_PLTE(png_ptr, info_ptr, palette, size_palette);
  png_write_info(png_ptr, info_ptr);
}

#define TEXT_SIZE 5

static void writeTable(png_structp png_ptr,png_infop info_ptr, char * pngfile)
{
  int height = (vraie_taille_palette+1)*TEXT_SIZE;
  int width = 50;
  unsigned char * row2 = calloc(height*width, sizeof(unsigned char));

  int i, j, k;
  char str[10000], str2[50];
  str[0] = '\0';

  for (i = 0; i < (vraie_taille_palette+1); i++) {
    if ((i < vraie_taille_palette) && ((i%4) == 0)) {
      sprintf(str2, " -fill blue -draw \"text 5,%d '< %g'\"", i*TEXT_SIZE, seuils_france54[i]);
    }
    if (i == vraie_taille_palette) {
      sprintf(str2, " -fill blue -draw \"text 3,%d 'Manque'\"", i*TEXT_SIZE);
    }
    strcat(str, str2);
    for (j = 0; j < TEXT_SIZE; j++) {
      for (k = 0; k < 50; k++) {
	row2[(i*TEXT_SIZE+j)*50+k] = i;
      }
    }
  }

  fprintf(stderr, "convert %s %s -\n", pngfile, str);
  writeInfo(png_ptr, info_ptr, width, height);

  png_bytep row_pointer = row2;
  int count = 0;

  for (count = 0; count < height; count++){
    png_write_row(png_ptr, row_pointer);
    row_pointer += width;
  }
  free(row2);
}

static void writeData(png_structp png_ptr,png_infop info_ptr,radar_data_t *od, char * filename, FILE * fp) {
  unsigned short *p = od->img.data;
  
  int height = od->img.nrows;
  int width = od->img.ncols;

  unsigned char * row = calloc(height*width, sizeof(unsigned char));
  int i;

  for (i = 0; i < height*width; i++) {
    double Z;
    if (p[i] == 255) {
      row[i] = vraie_taille_palette;
      continue;
    }
    assert(p[i] < od->img.scale.nvals);
    Z = od->img.scale.vals[p[i]];
    row[i] = lookup(Z);
  }
  
  /* -------------- */

  unsigned char * row2;
  if (kml) {
    varfl north = -DBL_MAX, south = DBL_MAX, east = -DBL_MAX, west = DBL_MAX;
    projPJ pj, pj4;
    projUV pin, pout, ppNW, ppSW, ppNE, ppSE; 

    char str[200];
    sprintf(str, "+proj=aeqd +lat_0=%g +lon_0=%g +ellps=WGS84", 
	    od->meta.radar.lat, od->meta.radar.lon);
    if (!(pj = pj_init_plus(str))) {
      fprintf(stderr, "FATAL: error init proj %s (%s)\n", str, pj_strerrno(pj_errno));
      exit(1);
    }

    pj4 = pj_latlong_from_proj(pj); 

    /* bounding box calculation */
    pout.v = od->meta.radar.lat*DEG_TO_RAD;
    pout.u = od->meta.radar.lon*DEG_TO_RAD;
    pin = pj_fwd(pout,pj);

    pin.u -= od->img.ncols*od->img.psizex/2;
    pin.v += od->img.nrows*od->img.psizey/2;

    ppNW = pj_inv(pin,pj);
    fprintf(stderr, "NW corner: %g %g\n", ppNW.u/DEG_TO_RAD, ppNW.v/DEG_TO_RAD);
    
    pin.u += od->img.ncols*od->img.psizex;

    ppNE = pj_inv(pin,pj);
    fprintf(stderr, "NE corner: %g %g\n", ppNE.u/DEG_TO_RAD, ppNE.v/DEG_TO_RAD);

    pin.v -= od->img.nrows*od->img.psizey;
  
    ppSE = pj_inv(pin,pj);
    fprintf(stderr, "SE corner: %g %g\n", ppSE.u/DEG_TO_RAD, ppSE.v/DEG_TO_RAD);
  
    pin.u -= od->img.ncols*od->img.psizex;
  
    ppSW = pj_inv(pin,pj);
    fprintf(stderr, "SW corner: %g %g\n", ppSW.u/DEG_TO_RAD, ppSW.v/DEG_TO_RAD);

    pin.v += od->img.nrows*od->img.psizey;

    if (north < ppNW.v) north = ppNW.v;
    if (north < ppNE.v) north = ppNE.v;
    if (north < ppSE.v) north = ppSE.v;
    if (north < ppSW.v) north = ppSW.v;

    if (south > ppNW.v) south = ppNW.v;
    if (south > ppNE.v) south = ppNE.v;
    if (south > ppSE.v) south = ppSE.v;
    if (south > ppSW.v) south = ppSW.v;
    ppNW.v = ppNE.v = north;
    ppSE.v = ppSW.v = south;

    if (east < ppNW.u) east = ppNW.u;
    if (east < ppNE.u) east = ppNE.u;
    if (east < ppSE.u) east = ppSE.u;
    if (east < ppSW.u) east = ppSW.u;

    if (west > ppNW.u) west = ppNW.u;
    if (west > ppNE.u) west = ppNE.u;
    if (west > ppSE.u) west = ppSE.u;
    if (west > ppSW.u) west = ppSW.u;
    ppSE.u = ppNE.u = east;
    ppNW.u = ppSW.u = west;
  
    fprintf(fp, "  <GroundOverlay>\n");
    fprintf(fp, "    <name>Composite reflectivity</name>\n");
    fprintf(fp, "    <description>Composite reflectivity</description>\n");
    fprintf(fp, "    <Icon>\n");
    fprintf(fp, "      <href>%s.png</href>\n", (strrchr(filename, '/') == NULL ? filename : strrchr(filename, '/')+1));
    fprintf(fp, "    </Icon>\n");
    fprintf(fp, "    <LatLonBox>\n");
    fprintf(fp, "      <north>%g</north>\n", north/DEG_TO_RAD);
    fprintf(fp, "      <south>%g</south>\n", south/DEG_TO_RAD);
    fprintf(fp, "      <east>%g</east>\n", east/DEG_TO_RAD);
    fprintf(fp, "      <west>%g</west>\n", west/DEG_TO_RAD);
    fprintf(fp, "    </LatLonBox>\n");
    fprintf(fp, "  </GroundOverlay>\n");
    fprintf(fp, "  <Placemark><name>NW</name>\n");
    fprintf(fp, "  <description>Corner of original domain</description>\n");
    fprintf(fp, "  <Point>\n");
    fprintf(fp, "    <coordinates>%g,%g,0</coordinates>\n", ppNW.u/DEG_TO_RAD, ppNW.v/DEG_TO_RAD);
    fprintf(fp, "  </Point>\n");
    fprintf(fp, "  </Placemark>\n");
    fprintf(fp, "  <Placemark><name>NE</name>\n");
    fprintf(fp, "  <description>Corner of original domain</description>\n");
    fprintf(fp, "  <Point>\n");
    fprintf(fp, "    <coordinates>%g,%g,0</coordinates>\n", ppNE.u/DEG_TO_RAD, ppNE.v/DEG_TO_RAD);
    fprintf(fp, "  </Point>\n");
    fprintf(fp, "  </Placemark>\n");
    fprintf(fp, "  <Placemark><name>SE</name>\n");
    fprintf(fp, "  <description>Corner of original domain</description>\n");
    fprintf(fp, "  <Point>\n");
    fprintf(fp, "    <coordinates>%g,%g,0</coordinates>\n", ppSE.u/DEG_TO_RAD, ppSE.v/DEG_TO_RAD);
    fprintf(fp, "  </Point>\n");
    fprintf(fp, "  </Placemark>\n");
    fprintf(fp, "  <Placemark><name>SW</name>\n");
    fprintf(fp, "  <description>Corner of original domain</description>\n");
    fprintf(fp, "  <Point>\n");
    fprintf(fp, "    <coordinates>%g,%g,0</coordinates>\n", ppSW.u/DEG_TO_RAD, ppSW.v/DEG_TO_RAD);
    fprintf(fp, "  </Point>\n");
    fprintf(fp, "  </Placemark>\n");

    pj_transform(pj4, pj, 1, 1, &ppNW.u, &ppNW.v, NULL );
    pj_transform(pj4, pj, 1, 1, &ppNE.u, &ppNE.v, NULL );
    pj_transform(pj4, pj, 1, 1, &ppSW.u, &ppSW.v, NULL );
    fprintf(stderr, "x resolution for %d pixels: %g\n", width, (ppNE.u - ppNW.u)/width);
    fprintf(stderr, "y resolution for %d pixels: %g\n", height, (ppNW.v - ppSW.v)/height);
    
    // {
    //   int i, j;

    //   /* oversampling */
    //   int width0 = width*1.3; /*(ppNE.u - ppNW.u)/od->where.psizex;*/
    //   int height0 = height*1.3; /*(ppNW.v - ppSW.v)/od->where.psizey;*/

    //   fprintf(stderr, "x resolution for %d pixels: %g\n", width0, (ppNE.u - ppNW.u)/width0);
    //   fprintf(stderr, "y resolution for %d pixels: %g\n", height0, (ppNW.v - ppSW.v)/height0);

    //   fprintf(stderr, "False easting/northing: +x_0=%g +y_0=%g\n", ppSW.u, ppSW.v);

    //   fprintf(stderr, "Image projection: %s\n", pj_get_def(pj, 0));
    //   fprintf(stderr, "Latlon PROJ equivalent definition: %s\n", pj_get_def(pj4, 0));


    //   row2 = calloc(height0*width0, sizeof(unsigned char));
    //   for (i = 0; i < width0; i++) {
    // 	for (j = 0; j < height0; j++) {
    // 	  int k, l, jj, ll;

    // 	  pin.u = (i+0.5)*(ppNE.u - ppNW.u)/width0 + ppSW.u;
    // 	  pin.v = (j+0.5)*(ppNW.v - ppSW.v)/height0 + ppSW.v;
    // 	  pj_transform(pj4, pj, 1, 1, &pin.u, &pin.v, NULL);

    // 	  k = (pin.u - ppSW.u)/od->img.psizex - 0.5;
    // 	  l = (pin.v - ppSW.v)/od->img.psizey - 0.5;
 
    // 	  jj = height0 - j - 1;
    // 	  ll = height - l - 1;
	
    // 	  row2[i + jj*width0] = row[k + ll*width];
    // 	}
    //   }

    //   width = width0;
    //   height = height0;
    // }
  }
  //free(row);
  row2 = row;

  /* -------------- */

  writeInfo(png_ptr, info_ptr, width, height);

  png_bytep row_pointer = row2;
  int count = 0;

  for (count = 0; count < height; count++){
    png_write_row(png_ptr, row_pointer);
    row_pointer += width;
  }
  free(row2);
}

int write_png(radar_data_t * od, char * file)
{
  png_structp png_ptr = NULL;
  png_infop info_ptr = NULL;

  FILE * fp;

  if (od && kml) { 
    char * kmlfile = create_filename(file, ".kml", od);
    fp = fopen(kmlfile, "wb");
    if (!fp) return;
    free(kmlfile);
    fprintf(fp, "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n");
    fprintf(fp, "<kml xmlns=\"http://www.opengis.net/kml/2.2\">\n");
    fprintf(fp, "<Folder>\n");
    fprintf(fp, "  <name>Radar data at station %02d%03d</name>\n", 
	    od->wmoblock, od->wmostat);
    {
      char str[100];
      fprintf(fp, "  <description>Radar reflectivity at date: %04d/%02d/%02d %02d:%02d:%02d and elevation %g</description>\n", 
	      od->meta.year, od->meta.month, od->meta.day, od->meta.hour, od->meta.min, od->meta.sec, od->polar.elangle);
    }
  }

  {
    char * pngfile = create_filename(file, ".png", od);
    FILE* fp_png=initPngFile(pngfile, &png_ptr, &info_ptr); 
    if (!fp_png) return 1;

    if (od) {
      writeData(png_ptr, info_ptr, od, file, fp);
    } else {
      writeTable(png_ptr, info_ptr, pngfile);
    }
    deinit(&png_ptr, &info_ptr);
    free(pngfile);
    fclose(fp_png);
  }

  if (od && kml) { 
    fprintf(fp, "</Folder>\n");
    fprintf(fp, "</kml>\n");
    fclose(fp);
  }

  return 1;
}

/*===========================================================================*/

/* check for missing values and flag tables */
static char* fmt (int ind, varfl val, char *sval)
{
            if (val == MISSVAL) {
                strcpy (sval, "      missing");
            }
#if BUFR_OUT_BIN
            else if (desc_is_flagtable (ind)) {
                place_bin_values (val, ind, sval);
            }
#endif
            else {
                sprintf (sval, "%15.7f", val);
            }
	    return sval;
}

static bufrval_t * locvals = NULL;

static void bufr_close_locvals () {

    if (locvals == (bufrval_t*) NULL) return;
        
    if (locvals->vals != (varfl*) NULL) {
        free (locvals->vals);
        locvals->vals = (varfl*) NULL;
    }
    free (locvals);
    locvals = (bufrval_t*) NULL;
}

static bufrval_t* bufr_open_locvals () {


    if (locvals != (bufrval_t*)  NULL) {
        fprintf (stderr, "Value array not empty!\n");
        return (bufrval_t*) NULL;
    }
    locvals = malloc (sizeof (bufrval_t));

    if (locvals == (bufrval_t*)  NULL) {
        fprintf (stderr, "Error allocating memory for Value array!\n");
        return (bufrval_t*) NULL;
    }
    memset (locvals, 0, sizeof (bufrval_t));
    return locvals;
}

static int bufr_val_to_local (varfl val, int ind) {

    assert (locvals != (bufrval_t*) NULL);

    /* No output for special descriptors and sequence descs*/

    if (ind == _desc_special || des[ind]->id == SEQDESC) return 1;

    return bufr_val_to_array (&(locvals->vals), val, &(locvals->nvals));
}

static int bufr_val_to_local2 (varfl val, int ind) {

    assert (locvals != (bufrval_t*) NULL);

    return bufr_val_to_array (&(locvals->vals), val, &(locvals->nvals));
}

int bufr_header_to_file (sect_1_t* s1, const char* file) {

    FILE* fp;

    if (file) {
      fp = fopen (file, "w");
    } else {
      fp = stderr;
    }
    if (fp == NULL) {
        fprintf (stderr, "unable to open output file for section 1 !\n");
        return 0;
    }

    fprintf (fp, "%5d    BUFR edition                       \n", _bufr_edition);
    fprintf (fp, "%5d    master table used                  \n", s1->mtab);
    fprintf (fp, "%5d    originating subcenter              \n", s1->subcent);
    fprintf (fp, "%5d    generating center                  \n", s1->gencent);
    fprintf (fp, "%5d    original BUFR message              \n", s1->updsequ);
    fprintf (fp, "%5d    no optional section                \n", s1->opsec);
    fprintf (fp, "%5d    message type                       \n", s1->dcat);
    fprintf (fp, "%5d    local message subtype              \n", s1->dcatst);
    fprintf (fp, "%5d    version number of master table used\n", s1->vmtab);
    fprintf (fp, "%5d    version number of local table used \n", s1->vltab);
    fprintf (fp, "%5d    year                               \n", s1->year);
    fprintf (fp, "%5d    month                              \n", s1->mon);
    fprintf (fp, "%5d    day                                \n", s1->day);
    fprintf (fp, "%5d    hour                               \n", s1->hour);
    fprintf (fp, "%5d    minute                             \n", s1->min);
    /* new fields for bufr edition 4 */
    if (_bufr_edition >= 4) {
        fprintf (fp, "%5d    second                             \n", s1->sec);
        fprintf (fp, "%5d    international message subtype      \n", 
                 s1->idcatst);
    }

    if (file) fclose (fp);

    return 1;
}
 
static int bufr_char_to_file (varfl val, int ind)
{
    assert (ind == ccitt_special);

    if (val == 0) val = 0x20;
    
    fprintf (stderr, "%c", (int) val);
    return 1;
}

int bufr_file_out (varfl val, int ind)
{
    int depth = 1, nv, ok;
    char sval[80];
    char fname[512], tmp[80];
    char* unit;
    dd* d;
    static int nchars = 0;    /* number of characters for ccitt output */
    static int in_seq = 0;    /* flag to indicate sequences */
    static int count = 0;     /* counter for image files */
    bufrval_t* vals;
    int got_elt = 0;

    /* sanity checks */

    if (des[ind] == (desc*) NULL) {
        fprintf (stderr, "Data not available for bufr_file_out!\n");
        return 0;
    }

    /* element descriptor */

    if (des[ind]->id == ELDESC) {

        d = &(des[ind]->el->d);

	  got_elt = get_elt (d, val);

        /* output descriptor if not inside a sequence */

	  if (got_elt && !in_seq && ind != ccitt_special && !_replicating 
            && ind != add_f_special) {
	  fprintf (stderr, "%2d %2d %3d ", d->f, d->x, d->y);
	}

        /* descriptor without data (1.x.y, 2.x.y) or ascii) */

	{ 
            /* do we have a descriptor before the data element? */

	  if (got_elt) {
	    if (!in_seq && !_replicating && ind != add_f_special) {
	      fprintf (stderr, "%s            %s\n", 
		       fmt(ind, val, sval), des[ind]->el->elname);
            }
	  }
        }
	got_elt = 0;
    } /* end if ("Element descriptor") */

    /* sequence descriptor */

    else if (des[ind]->id == SEQDESC) {
      int got_seq;
        d = &(des[ind]->seq->d);

            in_seq ++;

	    locvals = bufr_open_locvals ();
            if (locvals == (bufrval_t*) NULL) return 0;

            if (!bufr_parse_out (des[ind]->seq->del, 0, des[ind]->seq->nel - 1,
                                 bufr_val_to_local, 0)) {
                bufr_close_locvals ();
                return 0;
            }

	    got_seq = get_seq (d, locvals->vals);
	    in_seq --;

        /* output descriptor if not inside another sequence descriptor */

	    if (got_seq && !in_seq && !_replicating) {
	      fprintf (stderr, "%2d %2d %3d ", d->f, d->x, d->y);

	      for (nv = 0; nv <  des[ind]->seq->nel; nv++) {
		int k = get_index(ELDESC, &des[ind]->seq->del[nv]);
		if (nv) fprintf(stderr, "          ");
		fprintf (stderr, "%2d %2d %3d %s  %s\n", des[ind]->seq->del[nv].f, des[ind]->seq->del[nv].x, 
			 des[ind]->seq->del[nv].y, fmt (ind, locvals->vals[nv], sval), des[k]->el->elname);
	      }
	    }
	    bufr_close_locvals ();

    } /* if ("seqdesc") */
    return 1;
}

int bufr_to_data (char* file, bufr_t* msg) {

    dd *dds = NULL;
    int ok;
    int ndescs, desch;

    /* open bitstreams for section 3 and 4 */

    desch = bufr_open_descsec_r(msg, &_subsets);

    ok = (desch >= 0);

    if (ok)
        ok = (bufr_open_datasect_r(msg) >= 0);

    /* calculate number of data descriptors  */
    
    ndescs = bufr_get_ndescs (msg);

    /* allocate memory and read data descriptors from bitstream */

    if (ok) 
        ok = bufr_in_descsec (&dds, ndescs, desch);


    /* output data and descriptors */

    if (ok)
      while (_subsets--) { 
        ok = bufr_parse_out (dds, 0, ndescs - 1, bufr_file_out, 1);
	subset_courant++;
      }
        
    /* close bitstreams and free descriptor array */

    if (dds != (dd*) NULL)
        free (dds);
    bufr_close_descsec_r (desch);
    bufr_close_datasect_r ();

    return ok;
}

int bufr_read_file2 (bufr_t* msg, const char* file) {

    FILE* fp;           /* file pointer to bufr file */
    char* bm;           /* pointer to memory holding bufr file */
    int len;

    /* open file */

    if ((strlen(file) == 1) && (file[0] == '-')) {
      fp = stdin;
    } else {  
      fp = fopen (file, "rb");
    }
    if (fp == NULL) {
        fprintf (stderr, "unable to open file '%s'\n", file);
        return 0;
    }

    /* get length of message */

    fseek (fp, 0L, SEEK_END);
    len = ftell (fp);
    fseek (fp, 0L, SEEK_SET);

    /* allocate memory and read message */

    bm = (char *) malloc ((size_t) len);
    if (bm == NULL) {
        fprintf (stderr, 
                 "unable to allocate %d bytes to hold BUFR-message !\n", len);
        fclose (fp);
        return 0;
    }
    if (fread (bm, 1, (size_t) len, fp) != (size_t) len) {
        fprintf (stderr, "Error reading BUFR message from file!\n");
        fclose (fp);
        free (bm);
        return 0;
    }

    fclose (fp);

    /* get raw bufr data */

    if (!bufr_get_sections (bm, len, msg)) {
        free (bm);
        return 0;
    }

    free (bm);
    return 1;
}
