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