4 * @author Mathis Rosenhauer, Deutsches Klimarechenzentrum
5 * @author Moritz Hanke, Deutsches Klimarechenzentrum
6 * @author Joerg Behrens, Deutsches Klimarechenzentrum
7 * @author Luis Kornblueh, Max-Planck-Institut fuer Meteorologie
12 * Mathis Rosenhauer, Luis Kornblueh
16 * Deutsches Klimarechenzentrum GmbH Max-Planck-Institut fuer Meteorologie
17 * Bundesstr. 45a Bundesstr. 53
18 * 20146 Hamburg 20146 Hamburg
21 * All rights reserved.
23 * Redistribution and use in source and binary forms, with or without
24 * modification, are permitted provided that the following conditions
27 * 1. Redistributions of source code must retain the above copyright
28 * notice, this list of conditions and the following disclaimer.
29 * 2. Redistributions in binary form must reproduce the above
30 * copyright notice, this list of conditions and the following
31 * disclaimer in the documentation and/or other materials provided
32 * with the distribution.
34 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
35 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
36 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
37 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
38 * COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
39 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
40 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
41 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
42 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
43 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
44 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
45 * OF THE POSSIBILITY OF SUCH DAMAGE.
47 * @section DESCRIPTION
49 * Adaptive Entropy Encoder
50 * Based on CCSDS documents 121.0-B-2 and 120.0-G-2
67 #include "encode_accessors.h"
69 /* Marker for Remainder Of Segment condition in zero block encoding */
72 static int m_get_block(struct aec_stream *strm);
74 static inline void emit(struct internal_state *state,
75 uint32_t data, int bits)
78 Emit sequence of bits.
81 if (bits <= state->bits) {
83 *state->cds += data << state->bits;
86 *state->cds++ += (uint64_t)data >> bits;
90 *state->cds++ = data >> bits;
93 state->bits = 8 - bits;
94 *state->cds = data << state->bits;
98 static inline void emitfs(struct internal_state *state, int fs)
101 Emits a fundamental sequence.
103 fs zero bits followed by one 1 bit.
107 if (fs < state->bits) {
108 state->bits -= fs + 1;
109 *state->cds += 1U << state->bits;
119 static inline void copy64(uint8_t *dst, uint64_t src)
131 static inline void emitblock_fs(struct aec_stream *strm, int k, int ref)
134 int used; /* used bits in 64 bit accumulator */
135 uint64_t acc; /* accumulator */
136 struct internal_state *state = strm->state;
138 acc = (uint64_t)*state->cds << 56;
139 used = 7 - state->bits;
141 for (i = ref; i < strm->block_size; i++) {
142 used += (state->block[i] >> k) + 1;
144 copy64(state->cds, acc);
149 acc |= 1ULL << (63 - used);
152 copy64(state->cds, acc);
153 state->cds += used >> 3;
154 state->bits = 7 - (used & 7);
157 static inline void emitblock(struct aec_stream *strm, int k, int ref)
160 Emit the k LSB of a whole block of input data.
164 struct internal_state *state = strm->state;
165 uint32_t *in = state->block + ref;
166 uint32_t *in_end = state->block + strm->block_size;
167 uint64_t mask = (1ULL << k) - 1;
168 uint8_t *o = state->cds;
177 while (p > k && in < in_end) {
179 a += ((uint64_t)(*in++) & mask) << p;
248 static void preprocess_unsigned(struct aec_stream *strm)
251 Preprocess RSI of unsigned samples.
253 Combining preprocessing and converting to uint32_t in one loop
254 is slower due to the data dependance on x_i-1.
258 struct internal_state *state = strm->state;
259 const uint32_t *x = state->data_raw;
260 uint32_t *d = state->data_pp;
261 uint32_t xmax = state->xmax;
262 uint32_t rsi = strm->rsi * strm->block_size - 1;
274 if (D <= xmax - x[0])
283 state->uncomp_len = (strm->block_size - 1) * strm->bits_per_sample;
286 static void preprocess_signed(struct aec_stream *strm)
289 Preprocess RSI of signed samples.
293 struct internal_state *state = strm->state;
294 uint32_t *d = state->data_pp;
295 int32_t *x = (int32_t *)state->data_raw;
296 uint64_t m = 1ULL << (strm->bits_per_sample - 1);
297 int64_t xmax = state->xmax;
298 int64_t xmin = state->xmin;
299 uint32_t rsi = strm->rsi * strm->block_size - 1;
301 *d++ = (uint32_t)x[0];
302 x[0] = (x[0] ^ m) - m;
305 x[1] = (x[1] ^ m) - m;
307 D = (int64_t)x[0] - x[1];
308 if (D <= xmax - x[0])
313 D = (int64_t)x[1] - x[0];
314 if (D <= x[0] - xmin)
323 state->uncomp_len = (strm->block_size - 1) * strm->bits_per_sample;
326 static uint64_t block_fs(struct aec_stream *strm, int k)
329 Sum FS of all samples in block for given splitting position.
334 struct internal_state *state = strm->state;
336 for (i = 0; i < strm->block_size; i += 8)
338 (uint64_t)(state->block[i + 0] >> k)
339 + (uint64_t)(state->block[i + 1] >> k)
340 + (uint64_t)(state->block[i + 2] >> k)
341 + (uint64_t)(state->block[i + 3] >> k)
342 + (uint64_t)(state->block[i + 4] >> k)
343 + (uint64_t)(state->block[i + 5] >> k)
344 + (uint64_t)(state->block[i + 6] >> k)
345 + (uint64_t)(state->block[i + 7] >> k);
348 fs -= (uint64_t)(state->block[0] >> k);
353 static uint32_t assess_splitting_option(struct aec_stream *strm)
356 Length of CDS encoded with splitting option and optimal k.
358 In Rice coding each sample in a block of samples is split at
359 the same position into k LSB and bits_per_sample - k MSB. The
360 LSB part is left binary and the MSB part is coded as a
361 fundamental sequence a.k.a. unary (see CCSDS 121.0-B-2). The
362 function of the length of the Coded Data Set (CDS) depending on
363 k has exactly one minimum (see A. Kiely, IPN Progress Report
366 To find that minimum with only a few costly evaluations of the
367 CDS length, we start with the k of the previous CDS. K is
368 increased and the CDS length evaluated. If the CDS length gets
369 smaller, then we are moving towards the minimum. If the length
370 increases, then the minimum will be found with smaller k.
372 For increasing k we know that we will gain block_size bits in
373 length through the larger binary part. If the FS lenth is less
374 than the block size then a reduced FS part can't compensate the
375 larger binary part. So we know that the CDS for k+1 will be
376 larger than for k without actually computing the length. An
377 analogue check can be done for decreasing k.
382 int this_bs; /* Block size of current block */
383 int no_turn; /* 1 if we shouldn't reverse */
384 int dir; /* Direction, 1 means increasing k, 0 decreasing k */
385 uint64_t len; /* CDS length for current k */
386 uint64_t len_min; /* CDS length minimum so far */
387 uint64_t fs_len; /* Length of FS part (not including 1s) */
389 struct internal_state *state = strm->state;
391 this_bs = strm->block_size - state->ref;
392 len_min = UINT64_MAX;
393 k = k_min = state->k;
398 fs_len = block_fs(strm, k);
399 len = fs_len + this_bs * (k + 1);
402 if (len_min < UINT64_MAX)
409 if (fs_len < this_bs || k >= state->kmax) {
419 if (fs_len >= this_bs || k == 0)
436 static uint32_t assess_se_option(struct aec_stream *strm)
439 Length of CDS encoded with Second Extension option.
441 If length is above limit just return UINT32_MAX.
447 struct internal_state *state = strm->state;
451 for (i = 0; i < strm->block_size; i+= 2) {
452 d = (uint64_t)state->block[i]
453 + (uint64_t)state->block[i + 1];
454 /* we have to worry about overflow here */
455 if (d > state->uncomp_len) {
459 len += d * (d + 1) / 2 + state->block[i + 1];
465 static void init_output(struct aec_stream *strm)
468 Direct output to next_out if next_out can hold a Coded Data
469 Set, use internal buffer otherwise.
472 struct internal_state *state = strm->state;
474 if (strm->avail_out > state->cds_len) {
475 if (!state->direct_out) {
476 state->direct_out = 1;
477 *strm->next_out = *state->cds;
478 state->cds = strm->next_out;
481 if (state->zero_blocks == 0 || state->direct_out) {
482 /* copy leftover from last block */
483 *state->cds_buf = *state->cds;
484 state->cds = state->cds_buf;
486 state->direct_out = 0;
496 static int m_flush_block_resumable(struct aec_stream *strm)
499 Slow and restartable flushing
501 struct internal_state *state = strm->state;
503 int n = MIN(state->cds - state->cds_buf - state->i, strm->avail_out);
504 memcpy(strm->next_out, state->cds_buf + state->i, n);
506 strm->avail_out -= n;
509 if (strm->avail_out == 0) {
512 state->mode = m_get_block;
517 static int m_flush_block(struct aec_stream *strm)
520 Flush block in direct_out mode by updating counters.
522 Fall back to slow flushing if in buffered mode.
525 struct internal_state *state = strm->state;
527 if (state->direct_out) {
528 n = state->cds - strm->next_out;
530 strm->avail_out -= n;
531 state->mode = m_get_block;
536 state->mode = m_flush_block_resumable;
540 static int m_encode_splitting(struct aec_stream *strm)
542 struct internal_state *state = strm->state;
545 emit(state, k + 1, state->id_len);
548 emit(state, state->block[0], strm->bits_per_sample);
550 emitblock_fs(strm, k, state->ref);
552 emitblock(strm, k, state->ref);
554 return m_flush_block(strm);
557 static int m_encode_uncomp(struct aec_stream *strm)
559 struct internal_state *state = strm->state;
561 emit(state, (1U << state->id_len) - 1, state->id_len);
562 emitblock(strm, strm->bits_per_sample, 0);
564 return m_flush_block(strm);
567 static int m_encode_se(struct aec_stream *strm)
571 struct internal_state *state = strm->state;
573 emit(state, 1, state->id_len + 1);
575 emit(state, state->block[0], strm->bits_per_sample);
577 for (i = 0; i < strm->block_size; i+= 2) {
578 d = state->block[i] + state->block[i + 1];
579 emitfs(state, d * (d + 1) / 2 + state->block[i + 1]);
582 return m_flush_block(strm);
585 static int m_encode_zero(struct aec_stream *strm)
587 struct internal_state *state = strm->state;
589 emit(state, 0, state->id_len + 1);
592 emit(state, state->zero_ref_sample, strm->bits_per_sample);
594 if (state->zero_blocks == ROS)
596 else if (state->zero_blocks >= 5)
597 emitfs(state, state->zero_blocks);
599 emitfs(state, state->zero_blocks - 1);
601 state->zero_blocks = 0;
602 return m_flush_block(strm);
605 static int m_select_code_option(struct aec_stream *strm)
608 Decide which code option to use.
613 struct internal_state *state = strm->state;
615 split_len = assess_splitting_option(strm);
616 se_len = assess_se_option(strm);
618 if (split_len < state->uncomp_len) {
619 if (split_len < se_len)
620 return m_encode_splitting(strm);
622 return m_encode_se(strm);
624 if (state->uncomp_len <= se_len)
625 return m_encode_uncomp(strm);
627 return m_encode_se(strm);
631 static int m_check_zero_block(struct aec_stream *strm)
634 Check if input block is all zero.
636 Aggregate consecutive zero blocks until we find !0 or reach the
637 end of a segment or RSI.
640 struct internal_state *state = strm->state;
641 uint32_t *p = state->block + state->ref;
642 uint32_t *end = state->block + strm->block_size;
644 while(p < end && *p == 0)
648 if (state->zero_blocks) {
649 /* The current block isn't zero but we have to emit a
650 * previous zero block first. The current block will be
651 * flagged and handled later.
653 state->block_nonzero = 1;
654 state->mode = m_encode_zero;
657 state->mode = m_select_code_option;
660 state->zero_blocks++;
661 if (state->zero_blocks == 1) {
662 state->zero_ref = state->ref;
663 state->zero_ref_sample = state->block[0];
665 if (state->blocks_avail == 0
666 || (strm->rsi - state->blocks_avail) % 64 == 0) {
667 if (state->zero_blocks > 4)
668 state->zero_blocks = ROS;
669 state->mode = m_encode_zero;
672 state->mode = m_get_block;
677 static int m_get_rsi_resumable(struct aec_stream *strm)
680 Get RSI while input buffer is short.
682 Let user provide more input. Once we got all input pad buffer
686 struct internal_state *state = strm->state;
689 if (strm->avail_in >= state->bytes_per_sample) {
690 state->data_raw[state->i] = state->get_sample(strm);
692 if (state->flush == AEC_FLUSH) {
695 state->data_raw[state->i] =
696 state->data_raw[state->i - 1];
697 while(++state->i < strm->rsi * strm->block_size);
699 emit(state, 0, state->bits);
700 if (strm->avail_out > 0) {
701 if (!state->direct_out)
702 *strm->next_out++ = *state->cds;
711 } while (++state->i < strm->rsi * strm->block_size);
713 if (strm->flags & AEC_DATA_PREPROCESS)
714 state->preprocess(strm);
716 return m_check_zero_block(strm);
719 static int m_get_block(struct aec_stream *strm)
722 Provide the next block of preprocessed input data.
724 Pull in a whole Reference Sample Interval (RSI) of data if
725 block buffer is empty.
728 struct internal_state *state = strm->state;
732 if (state->block_nonzero) {
733 state->block_nonzero = 0;
734 state->mode = m_select_code_option;
738 if (state->blocks_avail == 0) {
739 state->blocks_avail = strm->rsi - 1;
740 state->block = state->data_pp;
742 if (strm->avail_in >= state->rsi_len) {
743 state->get_rsi(strm);
744 if (strm->flags & AEC_DATA_PREPROCESS)
745 state->preprocess(strm);
747 return m_check_zero_block(strm);
750 state->mode = m_get_rsi_resumable;
755 state->uncomp_len = strm->block_size * strm->bits_per_sample;
757 state->block += strm->block_size;
758 state->blocks_avail--;
759 return m_check_zero_block(strm);
770 int aec_encode_init(struct aec_stream *strm)
772 struct internal_state *state;
774 if (strm->bits_per_sample > 32 || strm->bits_per_sample == 0)
775 return AEC_CONF_ERROR;
777 if (strm->block_size != 8
778 && strm->block_size != 16
779 && strm->block_size != 32
780 && strm->block_size != 64)
781 return AEC_CONF_ERROR;
783 if (strm->rsi > 4096)
784 return AEC_CONF_ERROR;
786 state = malloc(sizeof(struct internal_state));
788 return AEC_MEM_ERROR;
790 memset(state, 0, sizeof(struct internal_state));
793 if (strm->bits_per_sample > 16) {
794 /* 24/32 input bit settings */
797 if (strm->bits_per_sample <= 24
798 && strm->flags & AEC_DATA_3BYTE) {
799 state->bytes_per_sample = 3;
800 if (strm->flags & AEC_DATA_MSB) {
801 state->get_sample = aec_get_msb_24;
802 state->get_rsi = aec_get_rsi_msb_24;
804 state->get_sample = aec_get_lsb_24;
805 state->get_rsi = aec_get_rsi_lsb_24;
808 state->bytes_per_sample = 4;
809 if (strm->flags & AEC_DATA_MSB) {
810 state->get_sample = aec_get_msb_32;
811 state->get_rsi = aec_get_rsi_msb_32;
813 state->get_sample = aec_get_lsb_32;
814 state->get_rsi = aec_get_rsi_lsb_32;
818 else if (strm->bits_per_sample > 8) {
819 /* 16 bit settings */
821 state->bytes_per_sample = 2;
823 if (strm->flags & AEC_DATA_MSB) {
824 state->get_sample = aec_get_msb_16;
825 state->get_rsi = aec_get_rsi_msb_16;
827 state->get_sample = aec_get_lsb_16;
828 state->get_rsi = aec_get_rsi_lsb_16;
833 state->bytes_per_sample = 1;
835 state->get_sample = aec_get_8;
836 state->get_rsi = aec_get_rsi_8;
838 state->rsi_len = strm->rsi * strm->block_size * state->bytes_per_sample;
840 if (strm->flags & AEC_DATA_SIGNED) {
841 state->xmin = -(1ULL << (strm->bits_per_sample - 1));
842 state->xmax = (1ULL << (strm->bits_per_sample - 1)) - 1;
843 state->preprocess = preprocess_signed;
846 state->xmax = (1ULL << strm->bits_per_sample) - 1;
847 state->preprocess = preprocess_unsigned;
850 state->kmax = (1U << state->id_len) - 3;
852 state->data_pp = malloc(strm->rsi
855 if (state->data_pp == NULL)
856 return AEC_MEM_ERROR;
858 if (strm->flags & AEC_DATA_PREPROCESS) {
859 state->data_raw = malloc(strm->rsi
862 if (state->data_raw == NULL)
863 return AEC_MEM_ERROR;
865 state->data_raw = state->data_pp;
868 state->block = state->data_pp;
870 /* Largest possible CDS according to specs */
871 state->cds_len = (5 + 64 * 32) / 8 + 3;
872 state->cds_buf = malloc(state->cds_len);
873 if (state->cds_buf == NULL)
874 return AEC_MEM_ERROR;
879 state->cds = state->cds_buf;
882 state->mode = m_get_block;
887 int aec_encode(struct aec_stream *strm, int flush)
890 Finite-state machine implementation of the adaptive entropy
894 struct internal_state *state = strm->state;
896 state->flush = flush;
897 strm->total_in += strm->avail_in;
898 strm->total_out += strm->avail_out;
900 while (state->mode(strm) == M_CONTINUE);
902 if (state->direct_out) {
903 n = state->cds - strm->next_out;
905 strm->avail_out -= n;
907 *state->cds_buf = *state->cds;
908 state->cds = state->cds_buf;
909 state->direct_out = 0;
911 strm->total_in -= strm->avail_in;
912 strm->total_out -= strm->avail_out;
916 int aec_encode_end(struct aec_stream *strm)
918 struct internal_state *state = strm->state;
920 if (strm->flags & AEC_DATA_PREPROCESS)
921 free(state->data_raw);
922 free(state->data_pp);
923 free(state->cds_buf);
928 int aec_buffer_encode(struct aec_stream *strm)
932 status = aec_encode_init(strm);
933 if (status != AEC_OK)
935 status = aec_encode(strm, AEC_FLUSH);
936 if (strm->avail_in > 0)
937 status = AEC_DATA_ERROR;
939 aec_encode_end(strm);