aee.c 17.6 KB
Newer Older
1 2 3 4 5
/* Adaptive Entropy Encoder            */
/* CCSDS 121.0-B-1 and CCSDS 120.0-G-2 */

#include <stdio.h>
#include <stdlib.h>
6
#include <unistd.h>
7 8 9 10 11
#include <inttypes.h>
#include <string.h>

#include "libae.h"

12
#define ROS -1
13 14

#define MIN(a, b) (((a) < (b))? (a): (b))
15
#define MAX(a, b) (((a) > (b))? (a): (b))
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

enum
{
    M_NEW_BLOCK,
    M_GET_BLOCK,
    M_CHECK_ZERO_BLOCK,
    M_SELECT_CODE_OPTION,
    M_ENCODE_SPLIT,
    M_FLUSH_BLOCK,
    M_FLUSH_BLOCK_LOOP,
    M_ENCODE_UNCOMP,
    M_ENCODE_SE,
    M_ENCODE_ZERO,
};

typedef struct internal_state {
32 33 34
    int id_len;             /* bit length of code option identification key */
    int64_t last_in;        /* previous input for preprocessing */
    int64_t (*get_sample)(ae_streamp);
35 36 37
    int64_t xmin;           /* minimum integer for preprocessing */
    int64_t xmax;           /* maximum integer for preprocessing */
    int mode;               /* current mode of FSM */
38 39
    int i;                  /* counter for samples */
    int64_t *block_in;      /* input block buffer */
40 41
    uint8_t *block_out;     /* output block buffer */
    uint8_t *bp_out;        /* pointer to current output */
42
    int64_t total_blocks;
43 44
    int bitp;               /* bit pointer to the next unused bit in accumulator */
    int block_deferred;     /* there is a block in the input buffer
45
                               but we first have to emit a zero block */
46 47 48 49 50 51
    int ref;                /* length of reference sample in current block
                               i.e. 0 or 1 depending on whether the block has
                               a reference sample or not */
    int zero_ref;           /* current zero block has a reference sample */
    int64_t zero_ref_sample;/* reference sample of zero block */
    int zero_blocks;        /* number of contiguous zero blocks */
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
52 53 54
#ifdef PROFILE
    int *prof;
#endif
55 56
} encode_state;

Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
57 58 59
static int64_t get_lsb_32(ae_streamp strm)
{
    int64_t data;
60

61 62 63 64
    data = (strm->next_in[3] << 24)
        | (strm->next_in[2] << 16)
        | (strm->next_in[1] << 8)
        | strm->next_in[0];
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
65

66
    strm->next_in += 4;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
67
    strm->total_in += 4;
68
    strm->avail_in -= 4;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
69 70 71 72 73 74 75
    return data;
}

static int64_t get_lsb_16(ae_streamp strm)
{
    int64_t data;

76
    data = (strm->next_in[1] << 8) | strm->next_in[0];
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
77

78
    strm->next_in += 2;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
79
    strm->total_in += 2;
80
    strm->avail_in -= 2;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
81 82 83 84 85 86
    return data;
}

static int64_t get_msb_32(ae_streamp strm)
{
    int64_t data;
87

88 89 90 91
    data = (strm->next_in[0] << 24)
        | (strm->next_in[1] << 16)
        | (strm->next_in[2] << 8)
        | strm->next_in[3];
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
92

93
    strm->next_in += 4;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
94
    strm->total_in += 4;
95
    strm->avail_in -= 4;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
96 97 98 99 100 101 102
    return data;
}

static int64_t get_msb_16(ae_streamp strm)
{
    int64_t data;

103
    data = (strm->next_in[0] << 8) | strm->next_in[1];
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
104

105
    strm->next_in += 2;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
106
    strm->total_in += 2;
107
    strm->avail_in -= 2;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
108 109 110 111 112 113 114 115 116
    return data;
}

static int64_t get_8(ae_streamp strm)
{
    strm->avail_in--;
    strm->total_in++;
    return *strm->next_in++;
}
117

118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
int ae_encode_init(ae_streamp strm)
{
    int blklen;
    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;
    }
    strm->state = state;

137 138
    if (strm->bit_per_sample > 16)
    {
139
        state->id_len = 5;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
140 141
        if (strm->flags & AE_DATA_MSB)
            state->get_sample = get_msb_32;
142
        else
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
143
            state->get_sample = get_lsb_32;
144 145 146
    }
    else if (strm->bit_per_sample > 8)
    {
147
        state->id_len = 4;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
148 149
        if (strm->flags & AE_DATA_MSB)
            state->get_sample = get_msb_16;
150
        else
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
151
            state->get_sample = get_lsb_16;
152
    }
153
    else
154
    {
155
        state->id_len = 3;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
156
        state->get_sample = get_8;
157 158 159 160 161 162 163 164 165 166 167 168
    }

    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;
    }
169

Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
170 171 172 173 174 175 176 177 178
#ifdef PROFILE
    state->prof = (int *)malloc((strm->bit_per_sample + 2) * sizeof(int));
    if (state->prof == NULL)
    {
        return AE_MEM_ERROR;
    }
    memset(state->prof, 0, (strm->bit_per_sample + 2) * sizeof(int));
#endif

179
    state->block_in = (int64_t *)malloc(strm->block_size * sizeof(int64_t));
180 181 182 183 184
    if (state->block_in == NULL)
    {
        return AE_MEM_ERROR;
    }

185 186 187 188 189
    /* Zero blocks can span a segment and thus need up to segment_size
       bits in encoded block */
    blklen = MAX(strm->block_size * strm->bit_per_sample,
                 strm->segment_size + 10);
    blklen = (blklen + state->id_len) / 8 + 3;
190 191 192 193 194 195 196 197 198 199 200 201 202
    state->block_out = (uint8_t *)malloc(blklen);
    if (state->block_out == NULL)
    {
        return AE_MEM_ERROR;
    }
    state->bp_out = state->block_out;
    state->bitp = 8;

    strm->total_in = 0;
    strm->total_out = 0;

    state->mode = M_NEW_BLOCK;

Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
203
    state->total_blocks = 0;
204
    state->block_deferred = 0;
205
    state->zero_blocks = 0;
206 207 208 209 210 211
    state->zero_ref = 0;
    state->ref = 0;

    return AE_OK;
}

212 213 214 215 216 217 218 219 220 221 222 223 224 225
int ae_encode_end(ae_streamp strm)
{
    encode_state *state;

    state = strm->state;
#ifdef PROFILE
    free(state->prof);
#endif
    free(state->block_in);
    free(state->block_out);
    free(state);
    return AE_OK;
}

226
static inline void emit(encode_state *state, int64_t data, int bits)
227 228 229
{
    while(bits)
    {
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
230
        data &= ((1ULL << bits) - 1);
231 232 233
        if (bits <= state->bitp)
        {
            state->bitp -= bits;
234
            *state->bp_out += data << state->bitp;
235 236 237 238 239
            bits = 0;
        }
        else
        {
            bits -= state->bitp;
240
            *state->bp_out += data >> bits;
241 242 243 244 245 246
            *++state->bp_out = 0;
            state->bitp = 8;
        }
    }
}

247
static inline void emitfs(encode_state *state, int fs)
248
{
249
    emit(state, 1, fs + 1);
250
}
251

Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276
#ifdef PROFILE
static inline void profile_print(ae_streamp strm)
{
    int i, total;
    encode_state *state;

    state = strm->state;
    fprintf(stderr, "Blocks encoded by each coding option\n");
    fprintf(stderr, "Zero blocks:  %i\n", state->prof[0]);
    total = state->prof[0];
    fprintf(stderr, "Second Ext.:  %i\n", state->prof[strm->bit_per_sample+1]);
    total += state->prof[strm->bit_per_sample+1];
    fprintf(stderr, "FS:           %i\n", state->prof[1]);
    total += state->prof[1];
    for (i = 2; i < strm->bit_per_sample - 1; i++)
    {
        fprintf(stderr, "k = %02i:       %i\n", i-1, state->prof[i]);
        total += state->prof[i];
    }
    fprintf(stderr, "Uncompressed: %i\n", state->prof[strm->bit_per_sample]);
    total += state->prof[strm->bit_per_sample];
    fprintf(stderr, "Total blocks: %i\n", total);
}
#endif

277 278 279 280 281 282 283
int ae_encode(ae_streamp strm, int flush)
{
    /**
       Finite-state machine implementation of the adaptive entropy
       encoder.
    */

284 285 286
    int i, j, k, zb, this_bs;
    int64_t split_len;
    int64_t split_len_min, se_len, fs_len;
287
    int64_t d;
288
    int64_t theta, Delta;
289
    size_t avail_out, total_out;
290 291 292 293

    encode_state *state;

    state = strm->state;
294 295
    total_out = strm->total_out;
    avail_out = strm->avail_out;
296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327

    for (;;)
    {
        switch(state->mode)
        {
        case M_NEW_BLOCK:
            if (state->zero_blocks == 0)
            {
                /* copy leftover from last block */
                *state->block_out = *state->bp_out;
                state->bp_out = state->block_out;
            }

            if(state->block_deferred)
            {
                state->block_deferred = 0;
                state->mode = M_SELECT_CODE_OPTION;
                break;
            }

            state->i = 0;
            state->mode = M_GET_BLOCK;

        case M_GET_BLOCK:
            do
            {
                if (strm->avail_in == 0)
                {
                    if (flush == AE_FLUSH)
                    {
                        if (state->i > 0)
                        {
328
                            /* pad block with last sample if we have
329
                               a partial block */
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
330
                            state->block_in[state->i] = state->block_in[state->i - 1];
331 332 333 334 335 336 337
                        }
                        else
                        {
                            /* Pad last output byte with 1 bits
                               if user wants to flush, i.e. we got
                               all input there is.
                            */
338
                            emit(state, 0xff, state->bitp);
339
                            *strm->next_out++ = *state->bp_out;
340 341
                            avail_out--;
                            total_out++;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
342 343 344 345
#ifdef PROFILE
                            profile_print(strm);
#endif
                            goto req_buffer;
346 347
                        }
                    }
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
348 349 350 351
                    else
                    {
                        goto req_buffer;
                    }
352 353 354
                }
                else
                {
355
                    state->block_in[state->i] = state->get_sample(strm);
356 357 358 359
                }
            }
            while (++state->i < strm->block_size);

360 361
            state->total_blocks++;

362
            /* preprocess block if needed */
363
            if (strm->flags & AE_DATA_PREPROCESS)
364 365 366 367
            {
                /* If this is the first block in a segment
                   then we need to insert a reference sample.
                */
368
                if(state->total_blocks % strm->segment_size == 1)
369 370 371 372 373 374 375 376 377 378 379 380 381
                {
                    state->ref = 1;
                    state->last_in = state->block_in[0];
                }
                else
                {
                    state->ref = 0;
                }

                for (i = state->ref; i < strm->block_size; i++)
                {
                    theta = MIN(state->last_in - state->xmin,
                                state->xmax - state->last_in);
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
382
                    Delta = state->block_in[i] - state->last_in;
383
                    state->last_in = state->block_in[i];
384
                    if (0 <= Delta && Delta <= theta)
385 386 387
                    {
                        state->block_in[i] = 2 * Delta;
                    }
388
                    else if (-theta <= Delta && Delta < 0)
389 390 391 392
                    {
                        d = Delta < 0 ? -(uint64_t)Delta : Delta;
                        state->block_in[i] = 2 * d - 1;
                    }
393
                    else
394 395 396 397
                    {
                        state->block_in[i] = theta +
                            (Delta < 0 ? -(uint64_t)Delta : Delta);
                    }
398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417
                }
            }
            state->mode = M_CHECK_ZERO_BLOCK;

        case M_CHECK_ZERO_BLOCK:
            zb = 1;
            for (i = state->ref; i < strm->block_size && zb; i++)
                if (state->block_in[i] != 0) zb = 0;

            if (zb)
            {
                /* remember ref on first zero block */
                if (state->zero_blocks == 0)
                {
                    state->zero_ref = state->ref;
                    state->zero_ref_sample = state->block_in[0];
                }

                state->zero_blocks++;

418
                if (state->total_blocks % strm->segment_size == 0)
419
                {
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
420 421 422
#ifdef PROFILE
                    state->prof[0] += state->zero_blocks;
#endif
423 424
                    if (state->zero_blocks > 4)
                        state->zero_blocks = ROS;
425 426 427 428 429 430 431 432
                    state->mode = M_ENCODE_ZERO;
                    break;
                }
                state->mode = M_NEW_BLOCK;
                break;
            }
            else if (state->zero_blocks)
            {
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
433 434 435
#ifdef PROFILE
                state->prof[0] += state->zero_blocks;
#endif
436 437 438 439 440 441 442 443 444 445 446 447
                state->mode = M_ENCODE_ZERO;
                /* The current block isn't zero but we have to
                   emit a previous zero block first. The
                   current block has to be handled later.
                */
                state->block_deferred = 1;
                break;
            }

            state->mode = M_SELECT_CODE_OPTION;

        case M_SELECT_CODE_OPTION:
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
448 449 450
            /* If zero block isn't an option then count length of
               sample splitting options */

451
            /* Baseline is the size of an uncompressed block */
452
            split_len_min = (strm->block_size - state->ref) * strm->bit_per_sample;
453 454
            k = strm->bit_per_sample;

455 456 457 458
            /* Length of this block minus reference sample if present */
            this_bs = strm->block_size - state->ref;

            /* Add FS encoded to unencoded parts */
459 460
            for (j = 0; j < strm->bit_per_sample - 2; j++)
            {
461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477
#ifdef UNROLL_BLOCK_8
                fs_len = (state->block_in[1] >> j)
                    + (state->block_in[2] >> j)
                    + (state->block_in[3] >> j)
                    + (state->block_in[4] >> j)
                    + (state->block_in[5] >> j)
                    + (state->block_in[6] >> j)
                    + (state->block_in[7] >> j);
                if (state->ref == 0)
                    fs_len += (state->block_in[0] >> j);
#else
                fs_len = 0;
                for (i = state->ref; i < strm->block_size; i++)
                    fs_len += state->block_in[i] >> j;
#endif
                split_len = fs_len + this_bs * (j + 1);
                if (split_len < split_len_min)
478
                {
479
                    split_len_min = split_len;
480
                    k = j;
481

482
#if 0
483 484 485 486 487 488 489
                    if (fs_len < this_bs)
                    {
                        /* Next can't get better because what we lose
                           by additional uncompressed bits isn't
                           compensated by a smaller FS part. */
                        break;
                    }
490
                }
491 492
                else
                    break;
493 494 495
#else
            }
#endif
496 497 498
            }

            /* Count bits for 2nd extension */
499 500
            se_len = 1;
            for (i = 0; i < strm->block_size && split_len_min > se_len; i+= 2)
501 502
            {
                d = state->block_in[i] + state->block_in[i + 1];
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
503
                /* we have to worry about overflow here */
504
                if (d > split_len_min)
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
505 506 507
                    se_len = d;
                else
                    se_len += d * (d + 1) / 2 + state->block_in[i + 1];
508 509 510
            }

            /* Decide which option to use */
511
            if (split_len_min <= se_len)
512 513 514
            {
                if (k == strm->bit_per_sample)
                {
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
515 516 517
#ifdef PROFILE
                    state->prof[k]++;
#endif
518 519 520 521 522
                    state->mode = M_ENCODE_UNCOMP;
                    break;
                }
                else
                {
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
523 524 525
#ifdef PROFILE
                    state->prof[k + 1]++;
#endif
526 527 528 529 530
                    state->mode = M_ENCODE_SPLIT;
                }
            }
            else
            {
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
531 532 533
#ifdef PROFILE
                state->prof[strm->bit_per_sample + 1]++;
#endif
534 535 536 537
                state->mode = M_ENCODE_SE;
                break;
            }

538
            emit(state, k + 1, state->id_len);
539
            if (state->ref)
540
                emit(state, state->block_in[0], strm->bit_per_sample);
541 542

            for (i = state->ref; i < strm->block_size; i++)
543
                emitfs(state, state->block_in[i] >> k);
544 545

            for (i = state->ref; i < strm->block_size; i++)
546
                emit(state, state->block_in[i], k);
547 548 549 550 551 552 553

            state->mode = M_FLUSH_BLOCK;

        case M_FLUSH_BLOCK:
            if (strm->avail_in == 0 && flush == AE_FLUSH)
            {
                /* pad last byte with 1 bits */
554
                emit(state, 0xff, state->bitp);
555 556 557 558 559 560 561
            }
            state->i = 0;
            state->mode = M_FLUSH_BLOCK_LOOP;

        case M_FLUSH_BLOCK_LOOP:
            while(state->block_out + state->i < state->bp_out)
            {
562
                if (avail_out == 0)
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
563 564 565 566
                {
#ifdef PROFILE
                    profile_print(strm);
#endif
567
                    goto req_buffer;
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
568
                }
569

570
                *strm->next_out++ = state->block_out[state->i];
571 572
                avail_out--;
                total_out++;
573 574 575 576 577 578
                state->i++;
            }
            state->mode = M_NEW_BLOCK;
            break;

        case M_ENCODE_UNCOMP:
579
            emit(state, 0x1f, state->id_len);
580
            for (i = 0; i < strm->block_size; i++)
581
                emit(state, state->block_in[i], strm->bit_per_sample);
582 583 584 585 586

            state->mode = M_FLUSH_BLOCK;
            break;

        case M_ENCODE_SE:
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
587
            emit(state, 1, state->id_len + 1);
588
            if (state->ref)
589
                emit(state, state->block_in[0], strm->bit_per_sample);
590 591 592 593

            for (i = 0; i < strm->block_size; i+= 2)
            {
                d = state->block_in[i] + state->block_in[i + 1];
594
                emitfs(state, d * (d + 1) / 2 + state->block_in[i + 1]);
595 596 597 598 599 600
            }

            state->mode = M_FLUSH_BLOCK;
            break;

        case M_ENCODE_ZERO:
601
            emit(state, 0, state->id_len + 1);
602 603
            if (state->zero_ref)
            {
604
                emit(state, state->zero_ref_sample, strm->bit_per_sample);
605
            }
606 607 608 609 610 611 612 613
            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);
614 615 616 617 618 619 620 621 622 623
            state->zero_blocks = 0;
            state->mode = M_FLUSH_BLOCK;
            break;

        default:
            return AE_STREAM_ERROR;
        }
    }

req_buffer:
624 625
    strm->total_out = total_out;
    strm->avail_out = avail_out;
626 627
    return AE_OK;
}