external/io_grib1/WGRIB/BDSunpk.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <stddef.h>
   4 #include <limits.h>
   5 #include "grib.h"
   6 #include "pds4.h"
   7 #include "bms.h"
   8 #include "bds.h"
   9
  10 /* 1996                         wesley ebisuzaki
  11  *
  12  * Unpack BDS section
  13  *
  14  * input: *bits, pointer to packed integer data
  15  *        *bitmap, pointer to bitmap (undefined data), NULL if none
  16  *        n_bits, number of bits per packed integer
  17  *        n, number of data points (includes undefined data)
  18  *        ref, scale: flt[] = ref + scale*packed_int
  19  * output: *flt, pointer to output array
  20  *        undefined values filled with UNDEFINED
  21  *
  22  * note: code assumes an integer > 32 bits
  23  *
  24  * 7/98 v1.2.1 fix bug for bitmaps and nbit >= 25 found by Larry Brasfield
  25  * 2/01 v1.2.2 changed jj from long int to double
  26  * 3/02 v1.2.3 added unpacking extensions for spectral data
  27  *             Luis Kornblueh, MPIfM
  28  */
  29
  30 static unsigned int mask[] = {0,1,3,7,15,31,63,127,255};
  31 static unsigned int map_masks[8] = {128, 64, 32, 16, 8, 4, 2, 1};
  32 static double shift[9] = {1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0, 128.0, 256.0};
  33
  34 void BDS_unpack(float *flt, unsigned char *bds, unsigned char *bitmap,
  35         int n_bits, int n, double ref, double scale) {
  36
  37     unsigned char *bits;
  38
  39     int i, mask_idx, t_bits, c_bits, j_bits;
  40     unsigned int j, map_mask, tbits, jmask, bbits;
  41     double jj;
  42
  43     if (BDS_Harmonic(bds)) {
  44         bits = bds + 15;
  45         /* fill in global mean */
  46         *flt++ = BDS_Harmonic_RefValue(bds);
  47         n -= 1;
  48     }
  49     else {
  50         bits = bds + 11;
  51     }
  52
  53     tbits = bbits = 0;
  54
  55     /* assume integer has 32+ bits */
  56     if (n_bits <= 25) {
  57         jmask = (1 << n_bits) - 1;
  58         t_bits = 0;
  59
  60         if (bitmap) {
  61             for (i = 0; i < n; i++) {
  62                 /* check bitmap */
  63                 mask_idx = i & 7;
  64                 if (mask_idx == 0) bbits = *bitmap++;
  65                 if ((bbits & map_masks[mask_idx]) == 0) {
  66                     *flt++ = UNDEFINED;
  67                     continue;
  68                 }
  69
  70                 while (t_bits < n_bits) {
  71                     tbits = (tbits * 256) + *bits++;
  72                     t_bits += 8;
  73                 }
  74                 t_bits -= n_bits;
  75                 j = (tbits >> t_bits) & jmask;
  76                 *flt++ = ref + scale*j;
  77             }
  78         }
  79         else {
  80             for (i = 0; i < n; i++) {
  81                 while (t_bits < n_bits) {
  82                     tbits = (tbits * 256) + *bits++;
  83                     t_bits += 8;
  84                 }
  85                 t_bits -= n_bits;
  86                 flt[i] = (tbits >> t_bits) & jmask;
  87             }
  88             /* at least this vectorizes :) */
  89             for (i = 0; i < n; i++) {
  90                 flt[i] = ref + scale*flt[i];
  91             }
  92         }
  93     }
  94     else {
  95         /* older unoptimized code, not often used */
  96         c_bits = 8;
  97         map_mask = 128;
  98         while (n-- > 0) {
  99             if (bitmap) {
 100                 j = (*bitmap & map_mask);
 101                 if ((map_mask >>= 1) == 0) {
 102                     map_mask = 128;
 103                     bitmap++;
 104                 }
 105                 if (j == 0) {
 106                     *flt++ = UNDEFINED;
 107                     continue;
 108                 }
 109             }
 110
 111             jj = 0.0;
 112             j_bits = n_bits;
 113             while (c_bits <= j_bits) {
 114                 if (c_bits == 8) {
 115                     jj = jj * 256.0  + (double) (*bits++);
 116                     j_bits -= 8;
 117                 }
 118                 else {
 119                     jj = (jj * shift[c_bits]) + (double) (*bits & mask[c_bits]);
 120                     bits++;
 121                     j_bits -= c_bits;
 122                     c_bits = 8;
 123                 }
 124             }
 125             if (j_bits) {
 126                 c_bits -= j_bits;
 127                 jj = (jj * shift[j_bits]) + (double) ((*bits >> c_bits) & mask[j_bits]);
 128             }
 129             *flt++ = ref + scale*jj;
 130         }
 131     }
 132     return;
 133 }