aee.c 18 KB
Newer Older
1 2 3 4 5 6 7 8 9
/**
 * @file aee.c
 * @author Mathis Rosenhauer, Deutsches Klimarechenzentrum
 * @section DESCRIPTION
 *
 * Adaptive Entropy Encoder
 * Based on CCSDS documents 121.0-B-1 and 120.0-G-2
 *
 */
10 11 12

#include <stdio.h>
#include <stdlib.h>
13
#include <unistd.h>
14 15 16 17
#include <inttypes.h>
#include <string.h>

#include "libae.h"
18 19
#include "aee.h"
#include "aee_mutators.h"
20

21
/* Marker for Remainder Of Segment condition in zero block encoding */
22
#define ROS -1
23 24 25

#define MIN(a, b) (((a) < (b))? (a): (b))

26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
static int m_get_block(ae_streamp strm);
static int m_get_block_cautious(ae_streamp strm);
static int m_check_zero_block(ae_streamp strm);
static int m_select_code_option(ae_streamp strm);
static int m_flush_block(ae_streamp strm);
static int m_flush_block_cautious(ae_streamp strm);
static int m_encode_splitting(ae_streamp strm);
static int m_encode_uncomp(ae_streamp strm);
static int m_encode_se(ae_streamp strm);
static int m_encode_zero(ae_streamp strm);

/*
 *
 * Bit emitters
 *
 */
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
42

43
static inline void emit(encode_state *state, uint64_t data, int bits)
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
44
{
45
    if (bits <= state->bitp)
46
    {
47 48 49 50 51 52 53 54 55
        state->bitp -= bits;
        *state->out_bp += data << state->bitp;
    }
    else
    {
        bits -= state->bitp;
        *state->out_bp++ += data >> bits;

        while (bits & ~7)
56
        {
57 58
            bits -= 8;
            *state->out_bp++ = data >> bits;
59
        }
60 61 62

        state->bitp = 8 - bits;
        *state->out_bp = data << state->bitp;
63
    }
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
64 65
}

66
static inline void emitfs(encode_state *state, int fs)
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
67
{
68 69
    /**
       Emits a fundamental sequence.
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
70

71 72
       fs zero bits followed by one 1 bit.
     */
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
73

74 75
    for(;;)
    {
76
        if (fs < state->bitp)
77
        {
78
            state->bitp -= fs + 1;
79 80 81 82 83 84 85 86 87 88
            *state->out_bp += 1 << state->bitp;
            break;
        }
        else
        {
            fs -= state->bitp;
            *++state->out_bp = 0;
            state->bitp = 8;
        }
    }
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
89
}
90

Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123
static inline void emitblock(ae_streamp strm, int k, int skip)
{
    int i;
    uint64_t acc;
    encode_state *state = strm->state;
    int64_t *in = state->in_block + skip;
    int64_t *in_end = state->in_block + strm->block_size;
    int64_t mask = (1ULL << k) - 1;
    uint8_t *out = state->out_bp;

    acc = *out;

    while(in < in_end)
    {
        acc <<= 56;
        state->bitp = (state->bitp % 8) + 56;

        while (state->bitp > k && in < in_end)
        {
            state->bitp -= k;
            acc += (*in++ & mask) << state->bitp;
        }

        for (i = 56; i > (state->bitp & ~7); i -= 8)
            *out++ = acc >> i;
        acc >>= i;
    }

    *out = acc;
    state->out_bp = out;
    state->bitp %= 8;
}

124
static inline void preprocess(ae_streamp strm)
125
{
126 127 128
    int i;
    int64_t theta, d, Delta;
    encode_state *state = strm->state;
129

130 131 132
    /* Insert reference samples into first block of Reference Sample
     * Interval.
     */
133
    if(state->in_total_blocks % strm->rsi == 0)
134
    {
135 136
        state->ref = 1;
        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 146 147 148 149 150 151 152 153 154 155 156 157
        theta = MIN(state->last_in - state->xmin,
                    state->xmax - state->last_in);
        Delta = state->in_block[i] - state->last_in;
        state->last_in = state->in_block[i];
        if (0 <= Delta && Delta <= theta)
        {
            state->in_block[i] = 2 * Delta;
        }
        else if (-theta <= Delta && Delta < 0)
        {
            d = Delta < 0 ? -(uint64_t)Delta : Delta;
            state->in_block[i] = 2 * d - 1;
        }
158
        else
159 160 161 162
        {
            state->in_block[i] = theta +
                (Delta < 0 ? -(uint64_t)Delta : Delta);
        }
163
    }
164
}
165

166 167 168 169 170 171 172 173 174 175 176
/*
 *
 * FSM functions
 *
 */

static int m_get_block(ae_streamp strm)
{
    encode_state *state = strm->state;

    if (strm->avail_out > state->out_blklen)
177
    {
178 179 180 181 182 183
        if (!state->out_direct)
        {
            state->out_direct = 1;
            *strm->next_out = *state->out_bp;
            state->out_bp = strm->next_out;
        }
184 185 186
    }
    else
    {
187 188 189 190 191 192 193
        if (state->zero_blocks == 0 || state->out_direct)
        {
            /* copy leftover from last block */
            *state->out_block = *state->out_bp;
            state->out_bp = state->out_block;
        }
        state->out_direct = 0;
194
    }
195

196
    if(state->block_deferred)
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
197
    {
198 199 200
        state->block_deferred = 0;
        state->mode = m_select_code_option;
        return M_CONTINUE;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
201 202
    }

203
    if (strm->avail_in >= state->in_blklen)
204
    {
205
        state->get_block(strm);
206

207 208 209 210 211 212 213
        if (strm->flags & AE_DATA_PREPROCESS)
            preprocess(strm);

        state->in_total_blocks++;
        return m_check_zero_block(strm);
    }
    else
214
    {
215 216
        state->i = 0;
        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 231 232 233
            if (state->flush == AE_FLUSH)
            {
                if (state->i > 0)
                {
234 235 236
                    /* Pad block with last sample if we have a partial
                     * block.
                     */
237 238 239 240 241 242 243 244 245 246
                    state->in_block[state->i] = state->in_block[state->i - 1];
                }
                else
                {
                    if (state->zero_blocks)
                    {
                        /* Output any remaining zero blocks */
                        state->mode = m_encode_zero;
                        return M_CONTINUE;
                    }
247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264

                    if ((strm->flags & AE_DATA_SZ_COMPAT)
                        && (state->in_total_blocks % strm->rsi != 0))
                    {
                        /* If user wants szip copatibility then we
                         * have to pad until but not including the
                         * next reference sample.
                         */
                        pad = 64 - (state->in_total_blocks % strm->rsi % 64);
                        state->in_total_blocks += pad;
                        state->zero_blocks = (pad > 4)? ROS: pad;
                        state->mode = m_encode_zero;
                        return M_CONTINUE;
                    }

                    /* Pad last output byte with 0 bits if user wants
                     * to flush, i.e. we got all input there is.
                     */
265 266 267 268 269
                    emit(state, 0, state->bitp);
                    if (state->out_direct == 0)
                        *strm->next_out++ = *state->out_bp;
                    strm->avail_out--;
                    strm->total_out++;
270

271 272 273 274 275 276 277
                    return M_EXIT;
                }
            }
            else
            {
                return M_EXIT;
            }
278 279 280
        }
        else
        {
281
            state->in_block[state->i] = state->get_sample(strm);
282 283
        }
    }
284 285 286 287 288 289 290
    while (++state->i < strm->block_size);

    if (strm->flags & AE_DATA_PREPROCESS)
        preprocess(strm);

    state->in_total_blocks++;
    return m_check_zero_block(strm);
291 292
}

293
static inline int m_check_zero_block(ae_streamp strm)
294
{
295 296 297 298 299 300 301 302
    int i;
    encode_state *state = strm->state;

    i = state->ref;
    while(i < strm->block_size && state->in_block[i] == 0)
        i++;

    if (i == strm->block_size)
303
    {
304 305 306 307 308 309
        /* remember ref on first zero block */
        if (state->zero_blocks == 0)
        {
            state->zero_ref = state->ref;
            state->zero_ref_sample = state->in_block[0];
        }
310

311
        state->zero_blocks++;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
312

313 314 315 316 317 318 319 320 321 322 323
        if (state->in_total_blocks % strm->rsi % 64 == 0)
        {
            if (state->zero_blocks > 4)
                state->zero_blocks = ROS;
            state->mode = m_encode_zero;
            return M_CONTINUE;
        }
        state->mode = m_get_block;
        return M_CONTINUE;
    }
    else if (state->zero_blocks)
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
324
    {
325 326 327 328
        /* The current block isn't zero but we have to emit a previous
         * zero block first. The current block will be handled
         * later.
         */
329 330 331
        state->block_deferred = 1;
        state->mode = m_encode_zero;
        return M_CONTINUE;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
332
    }
333 334
    state->mode = m_select_code_option;
    return M_CONTINUE;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
335 336
}

337
static inline int m_select_code_option(ae_streamp strm)
338
{
339 340
    int i, k, this_bs, looked_bothways, direction;
    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 345
    /* Length of this block minus reference sample (if present) */
    this_bs = strm->block_size - state->ref;
346

347 348 349 350
    split_len_min = INT64_MAX;
    i = state->k;
    direction = 1;
    looked_bothways = 0;
351

352 353 354
    /* Starting with splitting position of last block look left and
     * possibly right to find new minimum.
     */
355 356
    for (;;)
    {
357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380
        fs_len = (state->in_block[1] >> i)
            + (state->in_block[2] >> i)
            + (state->in_block[3] >> i)
            + (state->in_block[4] >> i)
            + (state->in_block[5] >> i)
            + (state->in_block[6] >> i)
            + (state->in_block[7] >> i);

        if (state->ref == 0)
            fs_len += (state->in_block[0] >> i);

        if (strm->block_size == 16)
            fs_len += (state->in_block[8] >> i)
                + (state->in_block[9] >> i)
                + (state->in_block[10] >> i)
                + (state->in_block[11] >> i)
                + (state->in_block[12] >> i)
                + (state->in_block[13] >> i)
                + (state->in_block[14] >> i)
                + (state->in_block[15] >> i);

        split_len = fs_len + this_bs * (i + 1);

        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 386
                 * the other direction.
                 */
387
                looked_bothways = 1;
388
            }
389 390
            split_len_min = split_len;
            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 398 399 400
                     * additional uncompressed bits isn't compensated
                     * by a smaller FS part. Vice versa if we are
                     * coming from the other direction.
                     */
401
                    if (looked_bothways)
402
                    {
403
                        break;
404
                    }
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
405 406
                    else
                    {
407 408 409
                        direction = -direction;
                        looked_bothways = 1;
                        i = state->k;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
410
                    }
411 412 413
                }
                else
                {
414 415 416 417 418
                    while (fs_len > 5 * this_bs)
                    {
                        i++;
                        fs_len /= 5;
                    }
419 420
                }
            }
421
            else if (fs_len > this_bs)
422
            {
423
                /* Since we started looking the other way there is no
424 425
                 * need to turn back.
                 */
426 427 428 429 430
                break;
            }
        }
        else
        {
431 432 433
            /* Stop looking for better option if we don't see any
             * improvement.
             */
434
                if (looked_bothways)
435
                {
436
                    break;
437 438 439
                }
                else
                {
440 441 442
                    direction = -direction;
                    looked_bothways = 1;
                    i = state->k;
443
                }
444 445 446 447 448 449
        }
        if (i + direction < 0
            || i + direction >= strm->bit_per_sample - 2)
        {
            if (looked_bothways)
                break;
450

451 452 453 454
            direction = -direction;
            looked_bothways = 1;
            i = state->k;
        }
455

456 457 458
        i += direction;
    }
    state->k = k;
459

460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475
    /* Count bits for 2nd extension */
    se_len = 1;
    for (i = 0; i < strm->block_size; i+= 2)
    {
        d = state->in_block[i] + state->in_block[i + 1];
        /* we have to worry about overflow here */
        if (d > split_len_min)
        {
            se_len = d;
            break;
        }
        else
        {
            se_len += d * (d + 1) / 2 + state->in_block[i + 1];
        }
    }
476

477 478
    /* Length of uncompressed block */
    uncomp_len = this_bs * strm->bit_per_sample;
479

480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504
    /* Decide which option to use */
    if (split_len_min < uncomp_len)
    {
        if (split_len_min <= se_len)
        {
            /* Splitting won - the most common case. */
            return m_encode_splitting(strm);
        }
        else
        {
            return m_encode_se(strm);
        }
    }
    else
    {
        if (uncomp_len <= se_len)
        {
            return m_encode_uncomp(strm);
        }
        else
        {
            return m_encode_se(strm);
        }
    }
}
505

506 507 508 509 510
static inline int m_encode_splitting(ae_streamp strm)
{
    int i;
    encode_state *state = strm->state;
    int k = state->k;
511

512 513 514
    emit(state, k + 1, state->id_len);
    if (state->ref)
        emit(state, state->in_block[0], strm->bit_per_sample);
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
515

516 517
    for (i = state->ref; i < strm->block_size; i++)
        emitfs(state, state->in_block[i] >> k);
518

Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
519 520
    if (k)
        emitblock(strm, k, state->ref);
521

522 523
    return m_flush_block(strm);
}
524

525 526 527
static inline int m_encode_uncomp(ae_streamp strm)
{
    encode_state *state = strm->state;
528

529
    emit(state, (1 << state->id_len) - 1, state->id_len);
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
530
    emitblock(strm, strm->bit_per_sample, 0);
531

532 533
    return m_flush_block(strm);
}
534

535 536 537 538 539
static inline int m_encode_se(ae_streamp strm)
{
    int i;
    int64_t d;
    encode_state *state = strm->state;
540

541 542 543
    emit(state, 1, state->id_len + 1);
    if (state->ref)
        emit(state, state->in_block[0], strm->bit_per_sample);
544

545 546 547 548 549
    for (i = 0; i < strm->block_size; i+= 2)
    {
        d = state->in_block[i] + state->in_block[i + 1];
        emitfs(state, d * (d + 1) / 2 + state->in_block[i + 1]);
    }
550

551 552
    return m_flush_block(strm);
}
553

554 555 556
static inline int m_encode_zero(ae_streamp strm)
{
    encode_state *state = strm->state;
557

558
    emit(state, 0, state->id_len + 1);
559

560 561
    if (state->zero_ref)
        emit(state, state->zero_ref_sample, strm->bit_per_sample);
562

563 564 565 566 567 568
    if (state->zero_blocks == ROS)
        emitfs(state, 4);
    else if (state->zero_blocks >= 5)
        emitfs(state, state->zero_blocks);
    else
        emitfs(state, state->zero_blocks - 1);
569

570 571 572
    state->zero_blocks = 0;
    return m_flush_block(strm);
}
573

574 575 576 577
static inline int m_flush_block(ae_streamp strm)
{
    int n;
    encode_state *state = strm->state;
578

579 580 581 582 583 584 585 586 587
    if (state->out_direct)
    {
        n = state->out_bp - strm->next_out;
        strm->next_out += n;
        strm->avail_out -= n;
        strm->total_out += n;
        state->mode = m_get_block;
        return M_CONTINUE;
    }
588

589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617
    state->i = 0;
    state->mode = m_flush_block_cautious;
    return M_CONTINUE;
}

static inline int m_flush_block_cautious(ae_streamp strm)
{
    encode_state *state = strm->state;

    /* Slow restartable flushing */
    while(state->out_block + state->i < state->out_bp)
    {
        if (strm->avail_out == 0)
            return M_EXIT;

        *strm->next_out++ = state->out_block[state->i];
        strm->avail_out--;
        strm->total_out++;
        state->i++;
    }
    state->mode = m_get_block;
    return M_CONTINUE;
}

/*
 *
 * API functions
 *
 */
618

619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662
int ae_encode_init(ae_streamp strm)
{
    encode_state *state;

    /* Some sanity checks */
    if (strm->bit_per_sample > 32 || strm->bit_per_sample == 0)
    {
        return AE_ERRNO;
    }

    /* Internal state for encoder */
    state = (encode_state *) malloc(sizeof(encode_state));
    if (state == NULL)
    {
        return AE_MEM_ERROR;
    }
    memset(state, 0, sizeof(encode_state));
    strm->state = state;

    if (strm->bit_per_sample > 16)
    {
        /* 32 bit settings */
        state->id_len = 5;
        state->in_blklen = 4 * strm->block_size;

        if (strm->flags & AE_DATA_MSB)
            state->get_sample = get_msb_32;
        else
            state->get_sample = get_lsb_32;
    }
    else if (strm->bit_per_sample > 8)
    {
        /* 16 bit settings */
        state->id_len = 4;
        state->in_blklen = 2 * strm->block_size;

        if (strm->flags & AE_DATA_MSB)
        {
            state->get_sample = get_msb_16;

            if (strm->block_size == 8)
                state->get_block = get_block_msb_16_bs_8;
            else
                state->get_block = get_block_msb_16_bs_16;
663
        }
664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703
        else
            state->get_sample = get_lsb_16;
    }
    else
    {
        /* 8 bit settings */
        state->in_blklen = strm->block_size;
        state->id_len = 3;

        state->get_sample = get_8;

        if (strm->block_size == 8)
            state->get_block = get_block_8_bs_8;
        else
            state->get_block = get_block_8_bs_16;
    }

    if (strm->flags & AE_DATA_SIGNED)
    {
        state->xmin = -(1ULL << (strm->bit_per_sample - 1));
        state->xmax = (1ULL << (strm->bit_per_sample - 1)) - 1;
    }
    else
    {
        state->xmin = 0;
        state->xmax = (1ULL << strm->bit_per_sample) - 1;
    }

    state->in_block = (int64_t *)malloc(strm->block_size * sizeof(int64_t));
    if (state->in_block == NULL)
    {
        return AE_MEM_ERROR;
    }

    /* Largest possible block according to specs */
    state->out_blklen = (5 + 16 * 32) / 8 + 3;
    state->out_block = (uint8_t *)malloc(state->out_blklen);
    if (state->out_block == NULL)
    {
        return AE_MEM_ERROR;
704 705
    }

706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746
    strm->total_in = 0;
    strm->total_out = 0;

    state->out_bp = state->out_block;
    *state->out_bp = 0;
    state->bitp = 8;
    state->mode = m_get_block;

    return AE_OK;
}

int ae_encode(ae_streamp strm, int flush)
{
    /**
       Finite-state machine implementation of the adaptive entropy
       encoder.
    */

    encode_state *state;
    state = strm->state;
    state->flush = flush;

    while (state->mode(strm) == M_CONTINUE);

    if (state->out_direct)
    {
        m_flush_block(strm);
        *state->out_block = *state->out_bp;
        state->out_bp = state->out_block;
        state->out_direct = 0;
    }
    return AE_OK;
}

int ae_encode_end(ae_streamp strm)
{
    encode_state *state = strm->state;

    free(state->in_block);
    free(state->out_block);
    free(state);
747 748
    return AE_OK;
}