/************/
/* mkmpng.c */
/****************************************************/
/*  to compile:  gcc mkmpng.c -o mkmpng -lm -lpng   */
/****************************************************/
/* Usage:   mkmpng  llx:lly:msz infile [outfile]    */
/*--------------------------------------------------*/
/*  drawing parameters (1st argument)               */
/*      llx :  Bottom X-coordinate (m)              */
/*      lly :  Left side Y-coordinate (m)           */
/*      msz :  Unit raster size (m)                 */
/*  infile  : Input data file path                  */
/*  outfile : PNG output filename                   */
/*             if omitted, output to stdout stream  */
/****************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <fcntl.h>
#include <math.h>
#include <png.h>
#define  LMX  4001
#define  LMY  4001

static int tcw[22][20] = {
          { 0x024059, 0x024a66, 0x025474, 0x025e82, 0x03678f,
            0x03719d, 0x037bab, 0x0485b8, 0x048fc6, 0x0499d4,
            0x0ca1dc, 0x1ca7de, 0x2cade0, 0x3bb3e3, 0x4bb9e5,
            0x5abfe7, 0x6ac5e9, 0x7acbeb, 0x89d1ee, 0x99d7f0 },
          { 0x004c4e, 0x00585a, 0x006466, 0x007072, 0x007c7e,
            0x00878a, 0x009397, 0x009fa3, 0x00abaf, 0x00b7bb,
            0x07bfc3, 0x17c3c7, 0x27c7cb, 0x37cbce, 0x47cfd2,
            0x57d3d6, 0x67d7da, 0x77dbde, 0x87e0e2, 0x97e4e5 },
          { 0x015742, 0x01644d, 0x017157, 0x017f61, 0x028c6c,
            0x029a76, 0x02a780, 0x02b58a, 0x02c295, 0x03cf9f,
            0x0bd7a7, 0x1adaad, 0x2adcb2, 0x3adfb8, 0x4ae1be,
            0x59e4c3, 0x69e6c9, 0x79e9cf, 0x89ecd4, 0x98eeda },
          { 0x055f36, 0x066d3f, 0x077c47, 0x078a50, 0x089958,
            0x09a861, 0x0ab669, 0x0bc572, 0x0cd37a, 0x0de282,
            0x15ea8a, 0x24eb92, 0x33ed99, 0x42eea1, 0x51efa8,
            0x60f1b0, 0x6ff2b7, 0x7ef3bf, 0x8df5c6, 0x9cf6ce },
          { 0x0b642b, 0x0d7331, 0x0f8338, 0x11923f, 0x13a245,
            0x14b14c, 0x16c053, 0x18d059, 0x1adf60, 0x1cef67,
            0x24f76f, 0x32f778, 0x40f881, 0x4ef88a, 0x5cf994,
            0x6af99d, 0x78faa6, 0x87fab0, 0x95fbb9, 0xa3fbc2 },
          { 0x136720, 0x167625, 0x19862a, 0x1d962f, 0x20a634,
            0x23b639, 0x26c63e, 0x29d643, 0x2ce548, 0x2ff54d,
            0x37fd55, 0x44fd60, 0x51fd6b, 0x5dfe76, 0x6afe81,
            0x77fe8c, 0x84fe97, 0x91fea1, 0x9efeac, 0xabfeb7 },
          { 0x1d6716, 0x21771a, 0x25871d, 0x2a9721, 0x2ea724,
            0x33b728, 0x37c62b, 0x3cd62f, 0x40e633, 0x45f636,
            0x4dfe3e, 0x58fe4a, 0x64fe57, 0x6ffe63, 0x7bfe70,
            0x86fe7c, 0x92fe89, 0x9dfe95, 0xa8fea1, 0xb4feae },
          { 0x26650e, 0x2c7511, 0x328413, 0x389415, 0x3ea417,
            0x44b31a, 0x4ac31c, 0x50d31e, 0x56e220, 0x5cf223,
            0x64fa2b, 0x6efa38, 0x78fb46, 0x82fb54, 0x8cfb61,
            0x96fb6f, 0xa0fc7d, 0xaafc8a, 0xb4fc98, 0xbefda6 },
          { 0x316108, 0x387109, 0x40800a, 0x478f0c, 0x4f9e0d,
            0x56ad0e, 0x5ebc10, 0x65cb11, 0x6dda12, 0x74e913,
            0x7cf11b, 0x85f22a, 0x8df339, 0x96f447, 0x9ef556,
            0xa6f565, 0xaff673, 0xb7f782, 0xc0f891, 0xc8f99f },
          { 0x3b5c03, 0x446a04, 0x4d7904, 0x568705, 0x5f9506,
            0x68a306, 0x71b207, 0x7ac007, 0x83ce08, 0x8cdc08,
            0x94e410, 0x9be620, 0xa2e82f, 0xa9e93e, 0xb0eb4e,
            0xb7ed5d, 0xbdee6d, 0xc4f07c, 0xcbf28b, 0xd2f39b },
          { 0x445500, 0x4f6201, 0x597001, 0x647d01, 0x6f8a01,
            0x799701, 0x84a401, 0x8eb201, 0x99bf02, 0xa3cc02,
            0xabd40a, 0xb1d719, 0xb6d929, 0xbbdc39, 0xc1df49,
            0xc6e259, 0xcce468, 0xd1e778, 0xd6ea88, 0xdced98 },
          { 0x4d4d00, 0x595900, 0x656500, 0x717100, 0x7d7d00,
            0x898900, 0x959500, 0xa1a100, 0xadad00, 0xb9b900,
            0xc1c107, 0xc5c517, 0xc9c927, 0xcdcd37, 0xd1d147,
            0xd5d557, 0xd9d967, 0xdddd77, 0xe1e187, 0xe5e597 },
          { 0x554400, 0x624f01, 0x705901, 0x7d6401, 0x8a6f01,
            0x977901, 0xa48401, 0xb28e01, 0xbf9902, 0xcca302,
            0xd4ab0a, 0xd7b119, 0xd9b629, 0xdcbb39, 0xdfc149,
            0xe2c659, 0xe4cc68, 0xe7d178, 0xead688, 0xeddc98 },
          { 0x5c3b03, 0x6a4404, 0x794d04, 0x875605, 0x955f06,
            0xa36806, 0xb27107, 0xc07a07, 0xce8308, 0xdc8c08,
            0xe49410, 0xe69b20, 0xe8a22f, 0xe9a93e, 0xebb04e,
            0xedb75d, 0xeebd6d, 0xf0c47c, 0xf2cb8b, 0xf3d29b },
          { 0x613108, 0x713809, 0x80400a, 0x8f470c, 0x9e4f0d,
            0xad560e, 0xbc5e10, 0xcb6511, 0xda6d12, 0xe97413,
            0xf17c1b, 0xf2852a, 0xf38d39, 0xf49647, 0xf59e56,
            0xf5a665, 0xf6af73, 0xf7b782, 0xf8c091, 0xf9c89f },
          { 0x65260e, 0x752c11, 0x843213, 0x943815, 0xa43e17,
            0xb3441a, 0xc34a1c, 0xd3501e, 0xe25620, 0xf25c23,
            0xfa642b, 0xfa6e38, 0xfb7846, 0xfb8254, 0xfb8c61,
            0xfb966f, 0xfca07d, 0xfcaa8a, 0xfcb498, 0xfdbea6 },
          { 0x671d16, 0x77211a, 0x87251d, 0x972a21, 0xa72e24,
            0xb73328, 0xc6372b, 0xd63c2f, 0xe64033, 0xf64536,
            0xfe4d3e, 0xfe584a, 0xfe6457, 0xfe6f63, 0xfe7b70,
            0xfe867c, 0xfe9289, 0xfe9d95, 0xfea8a1, 0xfeb4ae },
          { 0x671320, 0x761625, 0x86192a, 0x961d2f, 0xa62034,
            0xb62339, 0xc6263e, 0xd62943, 0xe52c48, 0xf52f4d,
            0xfd3755, 0xfd4460, 0xfd516b, 0xfe5d76, 0xfe6a81,
            0xfe778c, 0xfe8497, 0xfe91a1, 0xfe9eac, 0xfeabb7 },
          { 0x640b2b, 0x730d31, 0x830f38, 0x92113f, 0xa21345,
            0xb1144c, 0xc01653, 0xd01859, 0xdf1a60, 0xef1c67,
            0xf7246f, 0xf73278, 0xf84081, 0xf84e8a, 0xf95c94,
            0xf96a9d, 0xfa78a6, 0xfa87b0, 0xfb95b9, 0xfba3c2 },
          { 0x5f0536, 0x6d063f, 0x7c0747, 0x8a0750, 0x990858,
            0xa80961, 0xb60a69, 0xc50b72, 0xd30c7a, 0xe20d82,
            0xea158a, 0xeb2492, 0xed3399, 0xee42a1, 0xef51a8,
            0xf160b0, 0xf26fb7, 0xf37ebf, 0xf58dc6, 0xf69cce },
          { 0x570142, 0x64014d, 0x710157, 0x7f0161, 0x8c026c,
            0x9a0276, 0xa70280, 0xb5028a, 0xc20295, 0xcf039f,
            0xd70ba7, 0xda1aad, 0xdc2ab2, 0xdf3ab8, 0xe14abe,
            0xe459c3, 0xe669c9, 0xe979cf, 0xec89d4, 0xee98da },
          { 0x4c004e, 0x58005a, 0x640066, 0x700072, 0x7c007e,
            0x87008a, 0x930097, 0x9f00a3, 0xab00af, 0xb700bb,
            0xbf07c3, 0xc317c7, 0xc727cb, 0xcb37ce, 0xcf47d2,
            0xd357d6, 0xd767da, 0xdb77de, 0xe087e2, 0xe497e5 } } ;

main(int argc, char *argv[])
{ int i, j, n, ii, jj, mx, my, ix, iy, lv, ls, col;
  int llx=0, lly=0, msz=200, ixs, iys, mszx, mszy, nmx, nmy;
  float vlow, vstp, p, q, r, s;
  float fi, fj, fx, fy, cx=+0.7071, cy=-0.7071, cz=0., vs=10.;
  float vnul, vsum, vmin, vmax, vdev, vv, vx, vy;
  float v00, v01, v10, v11, va, vb, vc, *vm, *pv;
  char buf[81], fnam[81];
  FILE *fin, *fout;
  png_structp png_ptr;
  png_infop info_ptr;
  png_bytep row_pointers[441];
  png_color_8 sig_bit;
  unsigned char *bm, *pb, cr,cg,cb;

  if (argc < 3) { printf("mkmpng: parameters missing\n"); exit(1); }
  sscanf(argv[1], "%d:%d:%d", &llx,&lly,&msz);
  if ((fin=fopen(argv[2],"r")) == NULL)
    { printf("mkmpng: input-file open fail\n"); exit(1); }

  if ((vm=(float *)malloc(sizeof(float)*LMX*LMY)) == NULL)
    { fclose(fin); printf("mkmpng: malloc (indata) fail\n"); exit(1); }
  if ((bm=(unsigned char *)malloc(3*651*441)) == NULL)
    { fclose(fin); printf("mkmpng: malloc (bitmap) fail\n"); exit(1); }
  row_pointers[440] = bm;
  for (i=440; i>0; i--) { row_pointers[i-1] = row_pointers[i] + 3*651; }
  sig_bit.red = 8; sig_bit.green = 8; sig_bit.blue = 8; sig_bit.alpha = 0;

  while (1) {
    if (fgets(buf,81,fin) == NULL)
      { fclose(fin); printf("mkmpng: no data found\n"); exit(1); }
    if (buf[0] != '#') break;
  }
  fgets(buf,81,fin);
  sscanf(buf, "%d %d %d %d %d %d %f", &ixs,&iys,&mszx,&mszy,&nmx,&nmy,&vnul);
  if ((nmx>LMX) || (nmy>LMY))
    { fclose(fin); printf("mkmpng: too large array size\n"); exit(1); }
  if ((nmx<2) || (nmy<2))
    { fclose(fin); printf("mkmpng: illegal array size\n"); exit(1); }
  if ((mszx<1) || (mszy<1))
    { fclose(fin); printf("mkmpng: invalid mesh-size\n"); exit(1); }
  n = 0; vsum = 0.;
  for (pv=vm,j=0; j<nmy; j++) {
    for (i=0; i<nmx; pv++,i++) {
      fscanf(fin, "%f", &vv); *pv = vv;
      if (vv != vnul) {
        if (n++ == 0) { vmin = vv; vmax = vv; vsum = vv; }
        else { vsum += vv; if (vv<vmin) vmin = vv; if (vv>vmax) vmax = vv; }
      }
    }
  }
  vsum /= (float)n; vdev = 0.;
  for (pv=vm,j=0; j<nmy; j++) {
    for (i=0; i<nmx; pv++,i++) {
      vv = *pv; if (vv != vnul) { vdev += (vv-vsum)*(vv-vsum); }
    }
  }
  vdev = sqrtf(vdev / (float)n);
  vlow = vsum - vdev*2.; vstp = vdev / 5.;
  fclose(fin);

  fx = (vs/1000.) * (float)mszx; fy = (vs/1000.) * (float)mszx;
  pb = bm; mx = nmx-1; my = nmy-1;
  for (ii=-20; ii<0; ii++) for (jj=-20; jj<631; jj++)
    { *pb++ = 0xc0;  *pb++ = 0xc0;  *pb++ = 0xc0; }
  for (ix=(llx-ixs),ii=0; ii<401; ix+=msz,ii++) {
    if (ix<0) { i = -1; }
    else { i = ix/mszx; fi = (float)(ix%mszx)/(float)mszx; }
    if ((i==mx) && (fi== 0.)) { i--; fi = 1.; }
    for (jj=-20; jj<0; jj++) { *pb++ = 0xc0;  *pb++ = 0xc0;  *pb++ = 0xc0; }
    for (iy=(lly-iys),jj=0; jj<401; iy+=msz,jj++) {
      if (iy<0) { j = -1; }
      else { j = iy/mszy; fj = (float)(iy%mszy)/(float)mszy; }
      if ((j==my) && (fj== 0.)) { j--; fj = 1.; }
      if ((i<0) || (j<0) || ((i+1)>mx) || ((j+1)>my))
        { *pb++ = 0xff;  *pb++ = 0xff;  *pb++ = 0xff; }
      else {
        pv = vm + nmx*j + i;
        v00 = *pv; v01 = *(pv+1); v10 = *(pv+nmx); v11 = *(pv+nmx+1);
        if ((v00==vnul) || (v01==vnul) || (v10==vnul) || (v11==vnul)) 
          { *pb++ = 0xff;  *pb++ = 0xff;  *pb++ = 0xff; }
        else {
          va = v01-v00; vb = v10-v00; vc = v00-v01-v10+v11;
          vx = va + vc*fj; vy = vb + vc*fi; vv = v00 + vb*fj + vx*fi;
          lv = (int)((vv-vlow)/vstp + 1.);
          if (lv<0) lv = 0; if (lv>21) lv = 21;
          p = vx/fx; q = vy/fy; r = sqrtf(p*p + q*q + 1.);
          s = (p*cx + q*cy - cz) / r;
          if (s<0.) { ls = (int)(15. - s*5.); if (ls>19) ls = 19; }
          else { ls = (int)((1.-s) * 15.); }
          col = tcw[lv][ls];
          *pb++ = (col>>16 & 0xff); *pb++ = (col>>8 & 0xff);
          *pb++ = (col & 0xff);
        }
      }
    }
    for (jj=401; jj<631; jj++) { *pb++ = 0xc0;  *pb++ = 0xc0;  *pb++ = 0xc0; }
  }
  for (ii=401; ii<421; ii++) for (jj=-20; jj<631; jj++)
    { *pb++ = 0xc0;  *pb++ = 0xc0;  *pb++ = 0xc0; }

  if ((fin=fopen("setgp.par", "r")) == NULL)
    { printf("mkmpng: aux-file open fail\n"); exit(1); }
  while (1) {
    if (fgets(buf,81,fin) == NULL)
      { fclose(fin); printf("mkmpng: aux-file error1\n"); exit(1); }
    sscanf(buf, " %d %d", &ixs,&iys);
    if ((ixs == 0) && (iys == 0)) break;
    for (n=0; n<10; n++,iys++) {
      fgets(buf,81,fin); ls = strlen(buf);
      for (i=0; i<ls; i++) {
        if (buf[i] == 'x') {
          pb = bm + (i + ixs + (441-iys)*651)*3;
          *pb++ = 0x00; *pb++ = 0x00; *pb++ = 0x00; 
        }
      }
    }
  }
  while (1) {
    if (fgets(buf,81,fin) == NULL)
      { fclose(fin); printf("mkmpng: aux-file error2\n"); exit(1); }
    sscanf(buf, " %d %d %d", &ixs,&iys,&n);
    if ((ixs == 0) && (iys == 0)) break;
    if (n == 0) {
      for (i=0; i<5; i++) for (j=iys; j>(iys-2-i-i); j--) {
        pb = bm + (i + ixs + (441-j)*651)*3;
        *pb++ = 0x00; *pb++ = 0x00; *pb++ = 0xff; 
      }
      for (i=5; i<10; i++) for (j=iys; j>(iys-20+i+i); j--) {
        pb = bm + (i + ixs + (441-j)*651)*3;
        *pb++ = 0x00; *pb++ = 0x00; *pb++ = 0xff; 
      }
    } else if (n == 1) {
      for (j=0; j<5; j++) for (i=ixs; i<(ixs+2+j+j); i++) {
        pb = bm + (i + (441-j-iys)*651)*3;
        *pb++ = 0x00; *pb++ = 0x00; *pb++ = 0xff; 
      }
      for (j=5; j<10; j++) for (i=ixs; i<(ixs+20-j-j); i++) {
        pb = bm + (i + (441-j-iys)*651)*3;
        *pb++ = 0x00; *pb++ = 0x00; *pb++ = 0xff; 
      }
    } else if (n == 2) {
      for (i=0; i<5; i++) for (j=iys; j<(iys+2+i+i); j++) {
        pb = bm + (i + ixs + (441-j)*651)*3;
        *pb++ = 0x00; *pb++ = 0x00; *pb++ = 0xff; 
      }
      for (i=5; i<10; i++) for (j=iys; j<(iys+20-i-i); j++) {
        pb = bm + (i + ixs + (441-j)*651)*3;
        *pb++ = 0x00; *pb++ = 0x00; *pb++ = 0xff; 
      }
    } else if (n == 3) {
      for (j=0; j<5; j++) for (i=ixs; i>(ixs-2-j-j); i--) {
        pb = bm + (i + (441-j-iys)*651)*3;
        *pb++ = 0x00; *pb++ = 0x00; *pb++ = 0xff; 
      }
      for (j=5; j<10; j++) for (i=ixs; i>(ixs-20+j+j); i--) {
        pb = bm + (i + (441-j-iys)*651)*3;
        *pb++ = 0x00; *pb++ = 0x00; *pb++ = 0xff; 
      }
    }
  }
  while (1) {
    if (fgets(buf,81,fin) == NULL)
      { fclose(fin); printf("mkmpng: aux-file error3\n"); exit(1); }
    sscanf(buf, " %d %d %d %d", &ixs,&iys,&ii,&jj);
    if ((ixs == 0) && (iys == 0)) break;
    for (j=0; j<jj; j++) for (i=0; i<ii; i++) {
      pb = bm + (i + ixs + (441-j-iys)*651)*3;
      *pb++ = 0xe0; *pb++ = 0xe0; *pb++ = 0xe0; 
    }
  }
  while (1) {
    if (fgets(buf,81,fin) == NULL)
      { fclose(fin); printf("mkmpng: aux-file error4\n"); exit(1); }
    sscanf(buf, " %d %d %d %d", &ixs,&iys,&ii,&jj);
    if ((ixs == 0) && (iys == 0)) break;
    if (jj == 8) { cr = cg = 0x00; cb = 0xff; }
    else if (jj == 0) cr = cg = cb = 0x00;
    else if (jj == 1) cr = cg = cb = 0x70;
    else if (jj == 2) cr = cg = cb = 0xc0;
    else if (jj == 3) cr = cg = cb = 0xe0;
    else  cr = cg = cb = 0xff;
    if (ii < 0) {
      for (jj=-ii,j=0; j<jj; j++) {
        pb = bm + (ixs + (441-j-iys)*651)*3;
        *pb++ = cr; *pb++ = cg; *pb++ = cb;
      }
    } else {
      for (i=0; i<ii; i++) {
        pb = bm + (i + ixs + (441-iys)*651)*3;
        *pb++ = cr; *pb++ = cg; *pb++ = cb;
      }
    }
  }

  png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
  if (png_ptr == NULL) { printf("mkmpng: create png-struc fail\n"); exit(1); }
  info_ptr = png_create_info_struct(png_ptr);
  if (info_ptr == NULL) {
    png_destroy_write_struct(&png_ptr, (png_infopp)NULL);
    printf("mkmpng: create png-struc fail\n"); exit(1); 
  }
  if (argc > 3) {
    if ((fout=fopen(argv[3],"wb")) == NULL)
      { printf("mkmpng: output-file open fail\n"); exit(1); }
  } else {
    fout = stdout;
  }
  png_init_io(png_ptr, fout);
  png_set_IHDR(png_ptr, info_ptr, 651, 441, 8, PNG_COLOR_TYPE_RGB,
               PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT,
               PNG_FILTER_TYPE_DEFAULT);
  png_set_sBIT(png_ptr, info_ptr, &sig_bit);
  png_write_info(png_ptr, info_ptr);
  png_write_image(png_ptr, row_pointers);
  png_write_end(png_ptr, info_ptr);
  png_destroy_write_struct(&png_ptr, &info_ptr);
  if (argc > 3) { fclose(fout); }
  free(bm); free(vm);
  exit(0);
}