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
    int i, j, 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
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
        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);

368
369
370
        if (strm->block_size > 8)
            for (j = 8; j < strm->block_size; j++)
                fs_len += state->in_block[j] >> i;
371
372
373
374

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

        if (split_len < split_len_min)
375
        {
376
            if (split_len_min < INT64_MAX)
377
            {
378
                /* We are moving towards the minimum so it cant be in
379
380
                 * the other direction.
                 */
381
                looked_bothways = 1;
382
            }
383
384
            split_len_min = split_len;
            k = i;
385

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

445
446
447
448
            direction = -direction;
            looked_bothways = 1;
            i = state->k;
        }
449

450
451
452
        i += direction;
    }
    state->k = k;
453

454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
    /* 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];
        }
    }
470

471
472
    /* Length of uncompressed block */
    uncomp_len = this_bs * strm->bit_per_sample;
473

474
475
476
    /* Decide which option to use */
    if (split_len_min < uncomp_len)
    {
477
        if (split_len_min < se_len)
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
        {
            /* 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);
        }
    }
}
499

500
501
502
503
504
static inline int m_encode_splitting(ae_streamp strm)
{
    int i;
    encode_state *state = strm->state;
    int k = state->k;
505

506
507
508
    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
509

510
511
    for (i = state->ref; i < strm->block_size; i++)
        emitfs(state, state->in_block[i] >> k);
512

Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
513
514
    if (k)
        emitblock(strm, k, state->ref);
515

516
517
    return m_flush_block(strm);
}
518

519
520
521
static inline int m_encode_uncomp(ae_streamp strm)
{
    encode_state *state = strm->state;
522

523
    emit(state, (1 << state->id_len) - 1, state->id_len);
Mathis Rosenhauer's avatar
Mathis Rosenhauer committed
524
    emitblock(strm, strm->bit_per_sample, 0);
525

526
527
    return m_flush_block(strm);
}
528

529
530
531
532
533
static inline int m_encode_se(ae_streamp strm)
{
    int i;
    int64_t d;
    encode_state *state = strm->state;
534

535
536
537
    emit(state, 1, state->id_len + 1);
    if (state->ref)
        emit(state, state->in_block[0], strm->bit_per_sample);
538

539
540
541
542
543
    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]);
    }
544

545
546
    return m_flush_block(strm);
}
547

548
549
550
static inline int m_encode_zero(ae_streamp strm)
{
    encode_state *state = strm->state;
551

552
    emit(state, 0, state->id_len + 1);
553

554
555
    if (state->zero_ref)
        emit(state, state->zero_ref_sample, strm->bit_per_sample);
556

557
558
559
560
561
562
    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);
563

564
565
566
    state->zero_blocks = 0;
    return m_flush_block(strm);
}
567

568
569
570
571
static inline int m_flush_block(ae_streamp strm)
{
    int n;
    encode_state *state = strm->state;
572

573
574
575
576
577
578
579
580
581
    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;
    }
582

583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
    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
 *
 */
612

613
614
615
616
617
618
619
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;
620
621
622
623
624
625
626
627
628

    if (strm->block_size != 8
        && strm->block_size != 16
        && strm->block_size != 32
        && strm->block_size != 64)
        return AE_ERRNO;

    if (strm->rsi > 4096)
        return AE_ERRNO;
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645

    /* 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)
646
        {
647
            state->get_sample = get_msb_32;
648
649
            state->get_block = get_block_msb_32;
        }
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
        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
666
                state->get_block = get_block_msb_16;
667
        }
668
669
670
671
672
673
674
675
676
677
678
679
680
681
        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
682
            state->get_block = get_block_8;
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
    }

    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 */
703
    state->out_blklen = (5 + 64 * 32) / 8 + 3;
704
705
706
707
    state->out_block = (uint8_t *)malloc(state->out_blklen);
    if (state->out_block == NULL)
    {
        return AE_MEM_ERROR;
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
747
748
749
750
    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);
751
752
    return AE_OK;
}