jpcpack.c

Go to the documentation of this file.

 
#include <stdlib.h>
#include <math.h>
#include "grib2_int.h"
 
static void
jpcpack_int(void *fld, int fld_is_double, g2int width, g2int height, g2int *idrstmpl,
        unsigned char *cpack, g2int *lcpack)
{
    g2int  *ifld = NULL;
    static float alog2 = ALOG2;       /*  ln(2.0) */
    g2int  j, nbits, imin, imax, maxdif;
    g2int  ndpts, nbytes, nsize, retry;
    float  bscale, dscale, rmax, rmin, temp;
    double  rmaxd, rmind;
    unsigned char *ctemp;
    float *ffld = fld;
    double *dfld = fld;
 
    LOG((2, "jpcpack_int() fld_is_double %d width %ld height %ld idrstmpl[1] %ld *lcpack %ld",
     fld_is_double, width, height, idrstmpl[1], *lcpack));
    
    ndpts = width * height;
    bscale = int_power(2.0, -idrstmpl[1]);
    dscale = int_power(10.0, idrstmpl[2]);
    LOG((3, "ndpts %ld bscale %g dscale %g", ndpts, bscale, dscale));
    
    /* Find max and min values in the data. */
    rmaxd = dfld[0];
    rmind = dfld[0];
    rmax = ffld[0];
    rmin = ffld[0];
    if (fld_is_double)
    {
    for (j = 1; j < ndpts; j++)
    {
        if (dfld[j] > rmaxd)
        rmaxd = dfld[j];
        if (dfld[j] < rmind)
        rmind = dfld[j];
    }
    if (idrstmpl[1] == 0)
        maxdif = (g2int)(rint(rmaxd * dscale) - rint(rmind * dscale));
    else
        maxdif = (g2int)rint((rmaxd - rmind) * dscale * bscale);
    }
    else
    {
    for (j = 1; j < ndpts; j++)
    {
        if (ffld[j] > rmax)
        rmax = ffld[j];
        if (ffld[j] < rmin)
        rmin = ffld[j];
    }
    if (idrstmpl[1] == 0)
        maxdif = (g2int)(rint(rmax * dscale) - rint(rmin * dscale));
    else
        maxdif = (g2int)rint((rmax - rmin) * dscale * bscale);
    }
    LOG((3, "rmax %g rmaxd %g rmin %g rmind %g", rmax, rmaxd, rmin, rmind));
 
    /* 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 (((fld_is_double && rmind != rmaxd) || (!fld_is_double && rmin != rmax)) && maxdif != 0)
    {
        ifld = 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 minumum number of bits
             *  in which the data will fit. */
            imin = (g2int)rint((fld_is_double ? rmind : rmin) * dscale);
            imax = (g2int)rint((fld_is_double ? rmaxd : rmax) * dscale);
            maxdif = imax - imin;
            temp = log((double)(maxdif + 1)) / alog2;
            nbits = (g2int)ceil(temp);
            /*   scale data */
        if (fld_is_double)
        {
        rmind = (float)imin;
        for(j = 0; j < ndpts; j++)
            ifld[j] = (g2int)rint(dfld[j] * dscale) - imin;
        }
        else
        {
        rmin = (float)imin;
        for(j = 0; j < ndpts; j++)
            ifld[j] = (g2int)rint(ffld[j] * dscale) - imin;
        }
        }
        else
        {
            /* Use binary scaling factor and calculate minumum number
             * of bits in which the data will fit. */
        if (fld_is_double)
        {
        rmind = rmind * dscale;
        rmaxd = rmaxd * dscale;
        maxdif = (g2int)rint((rmaxd - rmind) * bscale);
        }
        else
        {
        rmin = rmin * dscale;
        rmax = rmax * dscale;
        maxdif = (g2int)rint((rmax - rmin) * bscale);
        }
        
            temp = log((double)(maxdif + 1)) / alog2;
            nbits = (g2int)ceil(temp);
            /*   scale data */
        if (fld_is_double)
        {
        for (j = 0; j < ndpts; j++)
            ifld[j] = (g2int)rint(((dfld[j] * dscale) - rmind) * bscale);
        }
        else
        {
        for (j = 0; j < ndpts; j++)
            ifld[j] = (g2int)rint(((ffld[j] * dscale) - rmin) * bscale);
        }       
        }
 
        /* Pack data into full octets, then do JPEG 2000 encode and
         * calculate the length of the packed data in bytes. */
        retry = 0;
        nbytes = (nbits + 7) / 8;
        nsize = *lcpack;          /* needed for input to enc_jpeg2000 */
        ctemp = calloc(ndpts, nbytes);
        sbits(ctemp, ifld, 0, nbytes * 8, 0, ndpts);
        if ((*lcpack = (g2int)enc_jpeg2000(ctemp, width, height, nbits, idrstmpl[5],
                                           idrstmpl[6], retry, (char *)cpack, nsize)) <= 0)
        {
            printf("jpcpack: ERROR Packing JPC = %d\n", (int)*lcpack);
            if (*lcpack == -3)
            {
                retry = 1;
                if ((*lcpack = (g2int)enc_jpeg2000(ctemp, width, height, nbits, idrstmpl[5],
                                                   idrstmpl[6], retry, (char *)cpack, nsize)) <= 0)
                    printf("jpcpack: Retry Failed.\n");
                else
                    printf("jpcpack: Retry Successful.\n");
            }
        }
        free(ctemp);
    }
    else
    {
        nbits = 0;
        *lcpack = 0;
    }
 
    /* Fill in ref value and number of bits in Template 5.0. */
    if (fld_is_double)
    rmin = (float)rmind;
    mkieee(&rmin, idrstmpl, 1);   /* ensure reference value is IEEE format. */
    idrstmpl[3] = nbits;
    idrstmpl[4] = 0;         /* original data were reals */
    if (idrstmpl[5] == 0)
        idrstmpl[6] = 255;       /* lossy not used */
    if (ifld)
        free(ifld);
}
 
void
jpcpack(float *fld, g2int width, g2int height, g2int *idrstmpl,
        unsigned char *cpack, g2int *lcpack)
{
    jpcpack_int(fld, 0, width, height, idrstmpl, cpack, lcpack);
}
 
void
jpcpackd(double *fld, g2int width, g2int height, g2int *idrstmpl,
        unsigned char *cpack, g2int *lcpack)
{
    jpcpack_int(fld, 1, width, height, idrstmpl, cpack, lcpack);
}