#include <stdlib.h>
#include <math.h>
#include "grib2.h"

void pngpack(g2float *fld,g2int width,g2int height,g2int *idrstmpl,
             unsigned char *cpack,g2int *lcpack)
//$$$  SUBPROGRAM DOCUMENTATION BLOCK
//                .      .    .                                       .
// SUBPROGRAM:    pngpack
//   PRGMMR: Gilbert          ORG: W/NP11    DATE: 2003-08-27
//
// ABSTRACT: This subroutine packs up a data field into PNG image format.
//   After the data field is scaled, and the reference value is subtracted out,
//   it is treated as a grayscale image and passed to a PNG encoder.
//   It also fills in GRIB2 Data Representation Template 5.41 or 5.40010 with
//   the appropriate values.
//
// PROGRAM HISTORY LOG:
// 2003-08-27  Gilbert
//
// USAGE:    pngpack(g2float *fld,g2int width,g2int height,g2int *idrstmpl,
//                   unsigned char *cpack,g2int *lcpack);
//   INPUT ARGUMENT LIST:
//     fld[]    - Contains the data values to pack
//     width    - number of points in the x direction
//     height   - number of points in the y direction
//     idrstmpl - Contains the array of values for Data Representation
//                Template 5.41 or 5.40010
//                [0] = Reference value - ignored on input
//                [1] = Binary Scale Factor
//                [2] = Decimal Scale Factor
//                [3] = number of bits for each data value - ignored on input
//                [4] = Original field type - currently ignored on input
//                      Data values assumed to be reals.
//
//   OUTPUT ARGUMENT LIST:
//     idrstmpl - Contains the array of values for Data Representation
//                Template 5.41 or 5.40010
//                [0] = Reference value - set by pngpack routine.
//                [1] = Binary Scale Factor - unchanged from input
//                [2] = Decimal Scale Factor - unchanged from input
//                [3] = Number of bits containing each grayscale pixel value
//                [4] = Original field type - currently set = 0 on output.
//                      Data values assumed to be reals.
//     cpack    - The packed data field
//     lcpack   - length of packed field cpack.
//
// REMARKS: None
//
// ATTRIBUTES:
//   LANGUAGE: C
//   MACHINE:  IBM SP
//
//$$$
{
      g2int  *ifld;
      const g2float alog2=0.69314718f;       //  ln(2.0)
      g2int  j,nbits,imin,imax,maxdif;
      g2int  ndpts,nbytes;
      g2float  bscale,dscale,rmax,rmin,temp;
      unsigned char *ctemp;

      ifld=0;
      ndpts=width*height;
      bscale=(g2float)int_power(2.0,-idrstmpl[1]);
      dscale=(g2float)int_power(10.0,idrstmpl[2]);
//
//  Find max and min values in the data
//
      rmax=fld[0];
      rmin=fld[0];
      for (j=1;j<ndpts;j++) {
        if (fld[j] > rmax) rmax=fld[j];
        if (fld[j] < rmin) rmin=fld[j];
      }
      maxdif = (g2int)RINT( (rmax-rmin)*dscale*bscale );
//
//  If max and min values are not equal, pack up field.
//  If they are equal, we have a constant field, and the reference
//  value (rmin) is the value for each point in the field and
//  set nbits to 0.
//
      if (rmin != rmax  &&  maxdif != 0 ) {
        ifld=(g2int *)malloc(ndpts*sizeof(g2int));
        //
        //  Determine which algorithm to use based on user-supplied
        //  binary scale factor and number of bits.
        //
        if (idrstmpl[1] == 0) {
           //
           //  No binary scaling and calculate minimum number of
           //  bits in which the data will fit.
           //
           imin=(g2int)RINT(rmin*dscale);
           imax=(g2int)RINT(rmax*dscale);
           maxdif=imax-imin;
           temp=(g2float)log((double)(maxdif+1))/alog2;
           nbits=(g2int)ceil(temp);
           rmin=(g2float)imin;
           //   scale data
           for(j=0;j<ndpts;j++)
             ifld[j]=(g2int)RINT(fld[j]*dscale)-imin;
        }
        else {
           //
           //  Use binary scaling factor and calculate minimum number of
           //  bits in which the data will fit.
           //
           rmin=rmin*dscale;
           rmax=rmax*dscale;
           maxdif=(g2int)RINT((rmax-rmin)*bscale);
           temp=(g2float)log((double)(maxdif+1))/alog2;
           nbits=(g2int)ceil(temp);
           //   scale data
           for (j=0;j<ndpts;j++)
             ifld[j]=(g2int)RINT(((fld[j]*dscale)-rmin)*bscale);
        }
        //
        //  Pack data into full octets, then do PNG encode.
        //  and calculate the length of the packed data in bytes
        //
        if (nbits <= 8) {
            nbits=8;
        }
        else if (nbits <= 16) {
            nbits=16;
        }
        else if (nbits <= 24) {
            nbits=24;
        }
        else {
            nbits=32;
        }
        nbytes=(nbits/8)*ndpts;
        ctemp=calloc(nbytes,1);
        sbits(ctemp,ifld,0,nbits,0,ndpts);
        //
        //  Encode data into PNG Format.
        //
        *lcpack=(g2int)enc_png((char *)ctemp,width,height,nbits,(char *)cpack);
        if (*lcpack <= 0) {
           printf("pngpack: ERROR Packing PNG = %d\n",(int)*lcpack);
        }
        free(ctemp);

      }
      else {
        /* Force E and D to 0 to avoid compatibility issues */
        idrstmpl[1]=0;
        idrstmpl[2]=0;
        if( dscale != 1.0 )
        {
          rmin = (float)floor(rmin * dscale) / dscale;
        }
        nbits=0;
        *lcpack=0;
      }

//
//  Fill in ref value and number of bits in Template 5.0
//
      mkieee(&rmin,idrstmpl+0,1);   // ensure reference value is IEEE format
      idrstmpl[3]=nbits;
      idrstmpl[4]=0;         // original data were reals
      if (ifld != 0) free(ifld);

}