Line data    Source code

       1             : #include <stdio.h>
       2             : #include <stdlib.h>
       3             : #include <limits.h>
       4             : #include "grib2.h"
       5             : 
       6             : #include "libaec.h"
       7             : 
       8             : // Cf https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-42.shtml
       9             : // and https://github.com/erget/wgrib2/commit/07b0f6fcb9669e0e3285318f50d516731b6956b2
      10             : 
      11           1 : g2int aecunpack(unsigned char *cpack,g2int len,g2int *idrstmpl,g2int ndpts,
      12             :                 g2float *fld)
      13             : {
      14             : 
      15           1 :       g2int  iret = 0;
      16             :       g2float  refV;
      17             : 
      18           1 :       rdieee(idrstmpl+0,&refV,1);
      19           1 :       g2float bscale = (g2float)int_power(2.0,idrstmpl[1]);
      20           1 :       g2float dscale = (g2float)int_power(10.0,-idrstmpl[2]);
      21           1 :       g2float bdscale = bscale * dscale;
      22           1 :       g2float refD = refV * dscale;
      23             : 
      24           1 :       g2int nbits = idrstmpl[3];
      25             : //
      26             : //  if nbits equals 0, we have a constant field where the reference value
      27             : //  is the data value at each gridpoint
      28             : //
      29           1 :       if (nbits != 0) {
      30           1 :          int nbytes_per_sample = (nbits + 7) / 8;
      31           1 :          if( ndpts != 0 && nbytes_per_sample > INT_MAX / ndpts )
      32             :          {
      33           0 :              return 1;
      34             :          }
      35           1 :          g2int* ifld=(g2int *)calloc(ndpts,sizeof(g2int));
      36             :          // Was checked just before
      37             :          // coverity[integer_overflow,overflow_sink]
      38           1 :          unsigned char* ctemp=(unsigned char *)calloc((size_t)(ndpts) * nbytes_per_sample,1);
      39           1 :          if ( ifld == NULL || ctemp == NULL) {
      40           0 :             fprintf(stderr, "Could not allocate space in aecunpack.\n"
      41             :                     "Data field NOT unpacked.\n");
      42           0 :             free(ifld);
      43           0 :             free(ctemp);
      44           0 :             return(1);
      45             :          }
      46             : 
      47           1 :          struct aec_stream strm = {0};
      48           1 :          strm.flags = idrstmpl[5]; // CCSDS compression options mask
      49           1 :          strm.bits_per_sample = nbits;
      50           1 :          strm.block_size = idrstmpl[6];
      51           1 :          strm.rsi = idrstmpl[7]; // Restart interval
      52             : 
      53           1 :          strm.next_in = cpack;
      54           1 :          strm.avail_in = len;
      55           1 :          strm.next_out = ctemp;
      56           1 :          strm.avail_out = (size_t)(ndpts) * nbytes_per_sample;
      57             : 
      58             :          // Note: libaec doesn't seem to be very robust to invalid inputs...
      59           1 :          int status = aec_buffer_decode(&strm);
      60           1 :          if (status != AEC_OK)
      61             :          {
      62           0 :              fprintf(stderr, "aec_buffer_decode() failed with return code %d", status);
      63           0 :              iret = 1;
      64             :          }
      65             :          else
      66             :          {
      67           1 :              gbits(ctemp,ndpts * nbytes_per_sample,ifld,0,nbytes_per_sample*8,0,ndpts);
      68             :              g2int j;
      69      405901 :              for (j=0;j<ndpts;j++) {
      70      405900 :                 fld[j] = refD + bdscale*(g2float)(ifld[j]);
      71             :              }
      72             :          }
      73           1 :          free(ctemp);
      74           1 :          free(ifld);
      75             :       }
      76             :       else {
      77             :          g2int j;
      78           0 :          for (j=0;j<ndpts;j++) fld[j]=refD;
      79             :       }
      80             : 
      81           1 :       return(iret);
      82             : }