src/aee.c

   1 /**
   2  * @file aee.c
   3  * @author Mathis Rosenhauer, Deutsches Klimarechenzentrum
   4  * @section DESCRIPTION
   5  *
   6  * Adaptive Entropy Encoder
   7  * Based on CCSDS documents 121.0-B-1 and 120.0-G-2
   8  *
   9  */
  10
  11 #include <stdio.h>
  12 #include <stdlib.h>
  13 #include <unistd.h>
  14 #include <inttypes.h>
  15 #include <string.h>
  16
  17 #include "libae.h"
  18 #include "aee.h"
  19 #include "aee_mutators.h"
  20
  21 /* Marker for Remainder Of Segment condition in zero block encoding */
  22 #define ROS -1
  23
  24 #define MIN(a, b) (((a) < (b))? (a): (b))
  25
  26 static int m_get_block(ae_streamp strm);
  27 static int m_get_block_cautious(ae_streamp strm);
  28 static int m_check_zero_block(ae_streamp strm);
  29 static int m_select_code_option(ae_streamp strm);
  30 static int m_flush_block(ae_streamp strm);
  31 static int m_flush_block_cautious(ae_streamp strm);
  32 static int m_encode_splitting(ae_streamp strm);
  33 static int m_encode_uncomp(ae_streamp strm);
  34 static int m_encode_se(ae_streamp strm);
  35 static int m_encode_zero(ae_streamp strm);
  36
  37 /*
  38  *
  39  * Bit emitters
  40  *
  41  */
  42
  43 static inline void emit(encode_state *state, uint64_t data, int bits)
  44 {
  45     if (bits <= state->bitp)
  46     {
  47         state->bitp -= bits;
  48         *state->out_bp += data << state->bitp;
  49     }
  50     else
  51     {
  52         bits -= state->bitp;
  53         *state->out_bp++ += data >> bits;
  54
  55         while (bits & ~7)
  56         {
  57             bits -= 8;
  58             *state->out_bp++ = data >> bits;
  59         }
  60
  61         state->bitp = 8 - bits;
  62         *state->out_bp = data << state->bitp;
  63     }
  64 }
  65
  66 static inline void emitfs(encode_state *state, int fs)
  67 {
  68     /**
  69        Emits a fundamental sequence.
  70
  71        fs zero bits followed by one 1 bit.
  72      */
  73
  74     for(;;)
  75     {
  76         if (fs < state->bitp)
  77         {
  78             state->bitp -= fs + 1;
  79             *state->out_bp += 1 << state->bitp;
  80             break;
  81         }
  82         else
  83         {
  84             fs -= state->bitp;
  85             *++state->out_bp = 0;
  86             state->bitp = 8;
  87         }
  88     }
  89 }
  90
  91 static inline void emitblock(ae_streamp strm, int k, int skip)
  92 {
  93     int i;
  94     uint64_t acc;
  95     encode_state *state = strm->state;
  96     int64_t *in = state->in_block + skip;
  97     int64_t *in_end = state->in_block + strm->block_size;
  98     int64_t mask = (1ULL << k) - 1;
  99     uint8_t *out = state->out_bp;
 100
 101     acc = *out;
 102
 103     while(in < in_end)
 104     {
 105         acc <<= 56;
 106         state->bitp = (state->bitp % 8) + 56;
 107
 108         while (state->bitp > k && in < in_end)
 109         {
 110             state->bitp -= k;
 111             acc += (*in++ & mask) << state->bitp;
 112         }
 113
 114         for (i = 56; i > (state->bitp & ~7); i -= 8)
 115             *out++ = acc >> i;
 116         acc >>= i;
 117     }
 118
 119     *out = acc;
 120     state->out_bp = out;
 121     state->bitp %= 8;
 122 }
 123
 124 static inline void preprocess(ae_streamp strm)
 125 {
 126     int i;
 127     int64_t theta, d, Delta;
 128     encode_state *state = strm->state;
 129
 130     /* Insert reference samples into first block of Reference Sample
 131      * Interval.
 132      */
 133     if(state->in_total_blocks % strm->rsi == 0)
 134     {
 135         state->ref = 1;
 136         state->last_in = state->in_block[0];
 137     }
 138     else
 139     {
 140         state->ref = 0;
 141     }
 142
 143     for (i = state->ref; i < strm->block_size; i++)
 144     {
 145         theta = MIN(state->last_in - state->xmin,
 146                     state->xmax - state->last_in);
 147         Delta = state->in_block[i] - state->last_in;
 148         state->last_in = state->in_block[i];
 149         if (0 <= Delta && Delta <= theta)
 150         {
 151             state->in_block[i] = 2 * Delta;
 152         }
 153         else if (-theta <= Delta && Delta < 0)
 154         {
 155             d = Delta < 0 ? -(uint64_t)Delta : Delta;
 156             state->in_block[i] = 2 * d - 1;
 157         }
 158         else
 159         {
 160             state->in_block[i] = theta +
 161                 (Delta < 0 ? -(uint64_t)Delta : Delta);
 162         }
 163     }
 164 }
 165
 166 /*
 167  *
 168  * FSM functions
 169  *
 170  */
 171
 172 static int m_get_block(ae_streamp strm)
 173 {
 174     encode_state *state = strm->state;
 175
 176     if (strm->avail_out > state->out_blklen)
 177     {
 178         if (!state->out_direct)
 179         {
 180             state->out_direct = 1;
 181             *strm->next_out = *state->out_bp;
 182             state->out_bp = strm->next_out;
 183         }
 184     }
 185     else
 186     {
 187         if (state->zero_blocks == 0 || state->out_direct)
 188         {
 189             /* copy leftover from last block */
 190             *state->out_block = *state->out_bp;
 191             state->out_bp = state->out_block;
 192         }
 193         state->out_direct = 0;
 194     }
 195
 196     if(state->block_deferred)
 197     {
 198         state->block_deferred = 0;
 199         state->mode = m_select_code_option;
 200         return M_CONTINUE;
 201     }
 202
 203     if (strm->avail_in >= state->in_blklen)
 204     {
 205         state->get_block(strm);
 206
 207         if (strm->flags & AE_DATA_PREPROCESS)
 208             preprocess(strm);
 209
 210         state->in_total_blocks++;
 211         return m_check_zero_block(strm);
 212     }
 213     else
 214     {
 215         state->i = 0;
 216         state->mode = m_get_block_cautious;
 217     }
 218     return M_CONTINUE;
 219 }
 220
 221 static int m_get_block_cautious(ae_streamp strm)
 222 {
 223     int pad;
 224     encode_state *state = strm->state;
 225
 226     do
 227     {
 228         if (strm->avail_in == 0)
 229         {
 230             if (state->flush == AE_FLUSH)
 231             {
 232                 if (state->i > 0)
 233                 {
 234                     /* Pad block with last sample if we have a partial
 235                      * block.
 236                      */
 237                     state->in_block[state->i] = state->in_block[state->i - 1];
 238                 }
 239                 else
 240                 {
 241                     if (state->zero_blocks)
 242                     {
 243                         /* Output any remaining zero blocks */
 244                         state->mode = m_encode_zero;
 245                         return M_CONTINUE;
 246                     }
 247
 248                     if ((strm->flags & AE_DATA_SZ_COMPAT)
 249                         && (state->in_total_blocks % strm->rsi != 0))
 250                     {
 251                         /* If user wants szip copatibility then we
 252                          * have to pad until but not including the
 253                          * next reference sample.
 254                          */
 255                         pad = 64 - (state->in_total_blocks % strm->rsi % 64);
 256                         state->in_total_blocks += pad;
 257                         state->zero_blocks = (pad > 4)? ROS: pad;
 258                         state->mode = m_encode_zero;
 259                         return M_CONTINUE;
 260                     }
 261
 262                     /* Pad last output byte with 0 bits if user wants
 263                      * to flush, i.e. we got all input there is.
 264                      */
 265                     emit(state, 0, state->bitp);
 266                     if (state->out_direct == 0)
 267                         *strm->next_out++ = *state->out_bp;
 268                     strm->avail_out--;
 269                     strm->total_out++;
 270
 271                     return M_EXIT;
 272                 }
 273             }
 274             else
 275             {
 276                 return M_EXIT;
 277             }
 278         }
 279         else
 280         {
 281             state->in_block[state->i] = state->get_sample(strm);
 282         }
 283     }
 284     while (++state->i < strm->block_size);
 285
 286     if (strm->flags & AE_DATA_PREPROCESS)
 287         preprocess(strm);
 288
 289     state->in_total_blocks++;
 290     return m_check_zero_block(strm);
 291 }
 292
 293 static inline int m_check_zero_block(ae_streamp strm)
 294 {
 295     int i;
 296     encode_state *state = strm->state;
 297
 298     i = state->ref;
 299     while(i < strm->block_size && state->in_block[i] == 0)
 300         i++;
 301
 302     if (i == strm->block_size)
 303     {
 304         /* remember ref on first zero block */
 305         if (state->zero_blocks == 0)
 306         {
 307             state->zero_ref = state->ref;
 308             state->zero_ref_sample = state->in_block[0];
 309         }
 310
 311         state->zero_blocks++;
 312
 313         if (state->in_total_blocks % strm->rsi % 64 == 0)
 314         {
 315             if (state->zero_blocks > 4)
 316                 state->zero_blocks = ROS;
 317             state->mode = m_encode_zero;
 318             return M_CONTINUE;
 319         }
 320         state->mode = m_get_block;
 321         return M_CONTINUE;
 322     }
 323     else if (state->zero_blocks)
 324     {
 325         /* The current block isn't zero but we have to emit a previous
 326          * zero block first. The current block will be handled
 327          * later.
 328          */
 329         state->block_deferred = 1;
 330         state->mode = m_encode_zero;
 331         return M_CONTINUE;
 332     }
 333     state->mode = m_select_code_option;
 334     return M_CONTINUE;
 335 }
 336
 337 static inline int m_select_code_option(ae_streamp strm)
 338 {
 339     int i, k, this_bs, looked_bothways, direction;
 340     int64_t d, split_len, uncomp_len;
 341     int64_t split_len_min, se_len, fs_len;
 342     encode_state *state = strm->state;
 343
 344     /* Length of this block minus reference sample (if present) */
 345     this_bs = strm->block_size - state->ref;
 346
 347     split_len_min = INT64_MAX;
 348     i = state->k;
 349     direction = 1;
 350     looked_bothways = 0;
 351
 352     /* Starting with splitting position of last block look left and
 353      * possibly right to find new minimum.
 354      */
 355     for (;;)
 356     {
 357         fs_len = (state->in_block[1] >> i)
 358             + (state->in_block[2] >> i)
 359             + (state->in_block[3] >> i)
 360             + (state->in_block[4] >> i)
 361             + (state->in_block[5] >> i)
 362             + (state->in_block[6] >> i)
 363             + (state->in_block[7] >> i);
 364
 365         if (state->ref == 0)
 366             fs_len += (state->in_block[0] >> i);
 367
 368         if (strm->block_size == 16)
 369             fs_len += (state->in_block[8] >> i)
 370                 + (state->in_block[9] >> i)
 371                 + (state->in_block[10] >> i)
 372                 + (state->in_block[11] >> i)
 373                 + (state->in_block[12] >> i)
 374                 + (state->in_block[13] >> i)
 375                 + (state->in_block[14] >> i)
 376                 + (state->in_block[15] >> i);
 377
 378         split_len = fs_len + this_bs * (i + 1);
 379
 380         if (split_len < split_len_min)
 381         {
 382             if (split_len_min < INT64_MAX)
 383             {
 384                 /* We are moving towards the minimum so it cant be in
 385                  * the other direction.
 386                  */
 387                 looked_bothways = 1;
 388             }
 389             split_len_min = split_len;
 390             k = i;
 391
 392             if (direction == 1)
 393             {
 394                 if (fs_len < this_bs)
 395                 {
 396                     /* Next can't get better because what we lose by
 397                      * additional uncompressed bits isn't compensated
 398                      * by a smaller FS part. Vice versa if we are
 399                      * coming from the other direction.
 400                      */
 401                     if (looked_bothways)
 402                     {
 403                         break;
 404                     }
 405                     else
 406                     {
 407                         direction = -direction;
 408                         looked_bothways = 1;
 409                         i = state->k;
 410                     }
 411                 }
 412                 else
 413                 {
 414                     while (fs_len > 5 * this_bs)
 415                     {
 416                         i++;
 417                         fs_len /= 5;
 418                     }
 419                 }
 420             }
 421             else if (fs_len > this_bs)
 422             {
 423                 /* Since we started looking the other way there is no
 424                  * need to turn back.
 425                  */
 426                 break;
 427             }
 428         }
 429         else
 430         {
 431             /* Stop looking for better option if we don't see any
 432              * improvement.
 433              */
 434                 if (looked_bothways)
 435                 {
 436                     break;
 437                 }
 438                 else
 439                 {
 440                     direction = -direction;
 441                     looked_bothways = 1;
 442                     i = state->k;
 443                 }
 444         }
 445         if (i + direction < 0
 446             || i + direction >= strm->bit_per_sample - 2)
 447         {
 448             if (looked_bothways)
 449                 break;
 450
 451             direction = -direction;
 452             looked_bothways = 1;
 453             i = state->k;
 454         }
 455
 456         i += direction;
 457     }
 458     state->k = k;
 459
 460     /* Count bits for 2nd extension */
 461     se_len = 1;
 462     for (i = 0; i < strm->block_size; i+= 2)
 463     {
 464         d = state->in_block[i] + state->in_block[i + 1];
 465         /* we have to worry about overflow here */
 466         if (d > split_len_min)
 467         {
 468             se_len = d;
 469             break;
 470         }
 471         else
 472         {
 473             se_len += d * (d + 1) / 2 + state->in_block[i + 1];
 474         }
 475     }
 476
 477     /* Length of uncompressed block */
 478     uncomp_len = this_bs * strm->bit_per_sample;
 479
 480     /* Decide which option to use */
 481     if (split_len_min < uncomp_len)
 482     {
 483         if (split_len_min <= se_len)
 484         {
 485             /* Splitting won - the most common case. */
 486             return m_encode_splitting(strm);
 487         }
 488         else
 489         {
 490             return m_encode_se(strm);
 491         }
 492     }
 493     else
 494     {
 495         if (uncomp_len <= se_len)
 496         {
 497             return m_encode_uncomp(strm);
 498         }
 499         else
 500         {
 501             return m_encode_se(strm);
 502         }
 503     }
 504 }
 505
 506 static inline int m_encode_splitting(ae_streamp strm)
 507 {
 508     int i;
 509     encode_state *state = strm->state;
 510     int k = state->k;
 511
 512     emit(state, k + 1, state->id_len);
 513     if (state->ref)
 514         emit(state, state->in_block[0], strm->bit_per_sample);
 515
 516     for (i = state->ref; i < strm->block_size; i++)
 517         emitfs(state, state->in_block[i] >> k);
 518
 519     if (k)
 520         emitblock(strm, k, state->ref);
 521
 522     return m_flush_block(strm);
 523 }
 524
 525 static inline int m_encode_uncomp(ae_streamp strm)
 526 {
 527     encode_state *state = strm->state;
 528
 529     emit(state, (1 << state->id_len) - 1, state->id_len);
 530     emitblock(strm, strm->bit_per_sample, 0);
 531
 532     return m_flush_block(strm);
 533 }
 534
 535 static inline int m_encode_se(ae_streamp strm)
 536 {
 537     int i;
 538     int64_t d;
 539     encode_state *state = strm->state;
 540
 541     emit(state, 1, state->id_len + 1);
 542     if (state->ref)
 543         emit(state, state->in_block[0], strm->bit_per_sample);
 544
 545     for (i = 0; i < strm->block_size; i+= 2)
 546     {
 547         d = state->in_block[i] + state->in_block[i + 1];
 548         emitfs(state, d * (d + 1) / 2 + state->in_block[i + 1]);
 549     }
 550
 551     return m_flush_block(strm);
 552 }
 553
 554 static inline int m_encode_zero(ae_streamp strm)
 555 {
 556     encode_state *state = strm->state;
 557
 558     emit(state, 0, state->id_len + 1);
 559
 560     if (state->zero_ref)
 561         emit(state, state->zero_ref_sample, strm->bit_per_sample);
 562
 563     if (state->zero_blocks == ROS)
 564         emitfs(state, 4);
 565     else if (state->zero_blocks >= 5)
 566         emitfs(state, state->zero_blocks);
 567     else
 568         emitfs(state, state->zero_blocks - 1);
 569
 570     state->zero_blocks = 0;
 571     return m_flush_block(strm);
 572 }
 573
 574 static inline int m_flush_block(ae_streamp strm)
 575 {
 576     int n;
 577     encode_state *state = strm->state;
 578
 579     if (state->out_direct)
 580     {
 581         n = state->out_bp - strm->next_out;
 582         strm->next_out += n;
 583         strm->avail_out -= n;
 584         strm->total_out += n;
 585         state->mode = m_get_block;
 586         return M_CONTINUE;
 587     }
 588
 589     state->i = 0;
 590     state->mode = m_flush_block_cautious;
 591     return M_CONTINUE;
 592 }
 593
 594 static inline int m_flush_block_cautious(ae_streamp strm)
 595 {
 596     encode_state *state = strm->state;
 597
 598     /* Slow restartable flushing */
 599     while(state->out_block + state->i < state->out_bp)
 600     {
 601         if (strm->avail_out == 0)
 602             return M_EXIT;
 603
 604         *strm->next_out++ = state->out_block[state->i];
 605         strm->avail_out--;
 606         strm->total_out++;
 607         state->i++;
 608     }
 609     state->mode = m_get_block;
 610     return M_CONTINUE;
 611 }
 612
 613 /*
 614  *
 615  * API functions
 616  *
 617  */
 618
 619 int ae_encode_init(ae_streamp strm)
 620 {
 621     encode_state *state;
 622
 623     /* Some sanity checks */
 624     if (strm->bit_per_sample > 32 || strm->bit_per_sample == 0)
 625     {
 626         return AE_ERRNO;
 627     }
 628
 629     /* Internal state for encoder */
 630     state = (encode_state *) malloc(sizeof(encode_state));
 631     if (state == NULL)
 632     {
 633         return AE_MEM_ERROR;
 634     }
 635     memset(state, 0, sizeof(encode_state));
 636     strm->state = state;
 637
 638     if (strm->bit_per_sample > 16)
 639     {
 640         /* 32 bit settings */
 641         state->id_len = 5;
 642         state->in_blklen = 4 * strm->block_size;
 643
 644         if (strm->flags & AE_DATA_MSB)
 645             state->get_sample = get_msb_32;
 646         else
 647             state->get_sample = get_lsb_32;
 648     }
 649     else if (strm->bit_per_sample > 8)
 650     {
 651         /* 16 bit settings */
 652         state->id_len = 4;
 653         state->in_blklen = 2 * strm->block_size;
 654
 655         if (strm->flags & AE_DATA_MSB)
 656         {
 657             state->get_sample = get_msb_16;
 658
 659             if (strm->block_size == 8)
 660                 state->get_block = get_block_msb_16_bs_8;
 661             else
 662                 state->get_block = get_block_msb_16_bs_16;
 663         }
 664         else
 665             state->get_sample = get_lsb_16;
 666     }
 667     else
 668     {
 669         /* 8 bit settings */
 670         state->in_blklen = strm->block_size;
 671         state->id_len = 3;
 672
 673         state->get_sample = get_8;
 674
 675         if (strm->block_size == 8)
 676             state->get_block = get_block_8_bs_8;
 677         else
 678             state->get_block = get_block_8_bs_16;
 679     }
 680
 681     if (strm->flags & AE_DATA_SIGNED)
 682     {
 683         state->xmin = -(1ULL << (strm->bit_per_sample - 1));
 684         state->xmax = (1ULL << (strm->bit_per_sample - 1)) - 1;
 685     }
 686     else
 687     {
 688         state->xmin = 0;
 689         state->xmax = (1ULL << strm->bit_per_sample) - 1;
 690     }
 691
 692     state->in_block = (int64_t *)malloc(strm->block_size * sizeof(int64_t));
 693     if (state->in_block == NULL)
 694     {
 695         return AE_MEM_ERROR;
 696     }
 697
 698     /* Largest possible block according to specs */
 699     state->out_blklen = (5 + 16 * 32) / 8 + 3;
 700     state->out_block = (uint8_t *)malloc(state->out_blklen);
 701     if (state->out_block == NULL)
 702     {
 703         return AE_MEM_ERROR;
 704     }
 705
 706     strm->total_in = 0;
 707     strm->total_out = 0;
 708
 709     state->out_bp = state->out_block;
 710     *state->out_bp = 0;
 711     state->bitp = 8;
 712     state->mode = m_get_block;
 713
 714     return AE_OK;
 715 }
 716
 717 int ae_encode(ae_streamp strm, int flush)
 718 {
 719     /**
 720        Finite-state machine implementation of the adaptive entropy
 721        encoder.
 722     */
 723
 724     encode_state *state;
 725     state = strm->state;
 726     state->flush = flush;
 727
 728     while (state->mode(strm) == M_CONTINUE);
 729
 730     if (state->out_direct)
 731     {
 732         m_flush_block(strm);
 733         *state->out_block = *state->out_bp;
 734         state->out_bp = state->out_block;
 735         state->out_direct = 0;
 736     }
 737     return AE_OK;
 738 }
 739
 740 int ae_encode_end(ae_streamp strm)
 741 {
 742     encode_state *state = strm->state;
 743
 744     free(state->in_block);
 745     free(state->out_block);
 746     free(state);
 747     return AE_OK;
 748 }