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 #define EMITBLOCK_FS(ref) \
132 static inline void emitblock_fs_##ref(struct aec_stream *strm, \
136 int used; /* used bits in 64 bit accumulator */ \
137 uint64_t acc; /* accumulator */ \
138 struct internal_state *state = strm->state; \
140 acc = (uint64_t)*state->cds << 56; \
141 used = 7 - state->bits; \
143 for (i = ref; i < strm->block_size; i++) { \
144 used += (state->block[i] >> k) + 1; \
146 copy64(state->cds, acc); \
151 acc |= 1ULL << (63 - used); \
154 copy64(state->cds, acc); \
155 state->cds += used >> 3; \
156 state->bits = 7 - (used & 7); \
162 #define EMITBLOCK(ref) \
163 static inline void emitblock_##ref(struct aec_stream *strm, \
167 Emit the k LSB of a whole block of input data. \
172 struct internal_state *state = strm->state; \
173 uint32_t *in = state->block + ref; \
174 uint32_t *in_end = state->block + strm->block_size; \
175 uint64_t mask = (1ULL << k) - 1; \
176 uint8_t *o = state->cds; \
177 int p = state->bits; \
181 while(in < in_end) { \
185 while (p > k && in < in_end) { \
187 a += ((uint64_t)(*in++) & mask) << p; \
190 for (b = 56; b > (p & ~7); b -= 8) \
197 state->bits = p % 8; \
203 static void preprocess_unsigned(struct aec_stream *strm)
206 Preprocess RSI of unsigned samples.
208 Combining preprocessing and converting to uint32_t in one loop
209 is slower due to the data dependance on x_i-1.
213 struct internal_state *state = strm->state;
214 const uint32_t *x = state->data_raw;
215 uint32_t *d = state->data_pp;
216 uint32_t xmax = state->xmax;
217 uint32_t rsi = strm->rsi * strm->block_size - 1;
229 if (D <= xmax - x[0])
240 static void preprocess_signed(struct aec_stream *strm)
243 Preprocess RSI of signed samples.
247 struct internal_state *state = strm->state;
248 uint32_t *d = state->data_pp;
249 int32_t *x = (int32_t *)state->data_raw;
250 uint64_t m = 1ULL << (strm->bits_per_sample - 1);
251 int64_t xmax = state->xmax;
252 int64_t xmin = state->xmin;
253 uint32_t rsi = strm->rsi * strm->block_size - 1;
255 *d++ = (uint32_t)x[0];
256 x[0] = (x[0] ^ m) - m;
259 x[1] = (x[1] ^ m) - m;
261 D = (int64_t)x[0] - x[1];
262 if (D <= xmax - x[0])
267 D = (int64_t)x[1] - x[0];
268 if (D <= x[0] - xmin)
279 static uint64_t block_fs(struct aec_stream *strm, int k)
282 Sum FS of all samples in block for given splitting position.
287 struct internal_state *state = strm->state;
289 fs = (uint64_t)(state->block[1] >> k)
290 + (uint64_t)(state->block[2] >> k)
291 + (uint64_t)(state->block[3] >> k)
292 + (uint64_t)(state->block[4] >> k)
293 + (uint64_t)(state->block[5] >> k)
294 + (uint64_t)(state->block[6] >> k)
295 + (uint64_t)(state->block[7] >> k);
297 if (strm->block_size > 8)
298 for (j = 8; j < strm->block_size; j += 8)
300 (uint64_t)(state->block[j + 0] >> k)
301 + (uint64_t)(state->block[j + 1] >> k)
302 + (uint64_t)(state->block[j + 2] >> k)
303 + (uint64_t)(state->block[j + 3] >> k)
304 + (uint64_t)(state->block[j + 4] >> k)
305 + (uint64_t)(state->block[j + 5] >> k)
306 + (uint64_t)(state->block[j + 6] >> k)
307 + (uint64_t)(state->block[j + 7] >> k);
310 fs += (uint64_t)(state->block[0] >> k);
315 static int assess_splitting_option(struct aec_stream *strm)
318 Length of CDS encoded with splitting option and optimal k.
320 In Rice coding each sample in a block of samples is split at
321 the same position into k LSB and bits_per_sample - k MSB. The
322 LSB part is left binary and the MSB part is coded as a
323 fundamental sequence a.k.a. unary (see CCSDS 121.0-B-2). The
324 function of the length of the Coded Data Set (CDS) depending on
325 k has exactly one minimum (see A. Kiely, IPN Progress Report
328 To find that minimum with only a few costly evaluations of the
329 CDS length, we start with the k of the previous CDS. K is
330 increased and the CDS length evaluated. If the CDS length gets
331 smaller, then we are moving towards the minimum. If the length
332 increases, then the minimum will be found with smaller k.
334 For increasing k we know that we will gain block_size bits in
335 length through the larger binary part. If the FS lenth is less
336 than the block size then a reduced FS part can't compensate the
337 larger binary part. So we know that the CDS for k+1 will be
338 larger than for k without actually computing the length. An
339 analogue check can be done for decreasing k.
344 int this_bs; /* Block size of current block */
345 int no_turn; /* 1 if we shouldn't reverse */
346 int dir; /* Direction, 1 means increasing k, 0 decreasing k */
347 uint64_t len; /* CDS length for current k */
348 uint64_t len_min; /* CDS length minimum so far */
349 uint64_t fs_len; /* Length of FS part (not including 1s) */
351 struct internal_state *state = strm->state;
353 this_bs = strm->block_size - state->ref;
354 len_min = UINT64_MAX;
355 k = k_min = state->k;
356 no_turn = (k == 0) ? 1 : 0;
360 fs_len = block_fs(strm, k);
361 len = fs_len + this_bs * (k + 1);
364 if (len_min < UINT64_MAX)
371 if (fs_len < this_bs || k >= state->kmax) {
381 if (fs_len >= this_bs || k == 0)
398 static int assess_se_option(uint64_t limit, struct aec_stream *strm)
401 Length of CDS encoded with Second Extension option.
403 If length is above limit just return UINT64_MAX.
409 struct internal_state *state = strm->state;
413 for (i = 0; i < strm->block_size; i+= 2) {
414 d = (uint64_t)state->block[i]
415 + (uint64_t)state->block[i + 1];
416 /* we have to worry about overflow here */
421 len += d * (d + 1) / 2
422 + (uint64_t)state->block[i + 1];
434 static int m_flush_block_resumable(struct aec_stream *strm)
437 Slow and restartable flushing
439 struct internal_state *state = strm->state;
441 while(state->cds_buf + state->i < state->cds) {
442 if (strm->avail_out == 0)
445 *strm->next_out++ = state->cds_buf[state->i];
450 state->mode = m_get_block;
454 static int m_flush_block(struct aec_stream *strm)
457 Flush block in direct_out mode by updating counters.
459 Fall back to slow flushing if in buffered mode.
462 struct internal_state *state = strm->state;
464 if (state->direct_out) {
465 n = state->cds - strm->next_out;
467 strm->avail_out -= n;
468 strm->total_out += n;
469 state->mode = m_get_block;
474 state->mode = m_flush_block_resumable;
478 static int m_encode_splitting(struct aec_stream *strm)
480 struct internal_state *state = strm->state;
483 emit(state, k + 1, state->id_len);
487 emit(state, state->block[0], strm->bits_per_sample);
488 emitblock_fs_1(strm, k);
490 emitblock_1(strm, k);
494 emitblock_fs_0(strm, k);
496 emitblock_0(strm, k);
499 return m_flush_block(strm);
502 static int m_encode_uncomp(struct aec_stream *strm)
504 struct internal_state *state = strm->state;
506 emit(state, (1U << state->id_len) - 1, state->id_len);
507 emitblock_0(strm, strm->bits_per_sample);
509 return m_flush_block(strm);
512 static int m_encode_se(struct aec_stream *strm)
516 struct internal_state *state = strm->state;
518 emit(state, 1, state->id_len + 1);
520 emit(state, state->block[0], strm->bits_per_sample);
522 for (i = 0; i < strm->block_size; i+= 2) {
523 d = state->block[i] + state->block[i + 1];
524 emitfs(state, d * (d + 1) / 2 + state->block[i + 1]);
527 return m_flush_block(strm);
530 static int m_encode_zero(struct aec_stream *strm)
532 struct internal_state *state = strm->state;
534 emit(state, 0, state->id_len + 1);
537 emit(state, state->zero_ref_sample, strm->bits_per_sample);
539 if (state->zero_blocks == ROS)
541 else if (state->zero_blocks >= 5)
542 emitfs(state, state->zero_blocks);
544 emitfs(state, state->zero_blocks - 1);
546 state->zero_blocks = 0;
547 return m_flush_block(strm);
550 static int m_select_code_option(struct aec_stream *strm)
553 Decide which code option to use.
559 struct internal_state *state = strm->state;
561 uncomp_len = (strm->block_size - state->ref)
562 * strm->bits_per_sample;
563 split_len = assess_splitting_option(strm);
564 se_len = assess_se_option(split_len, strm);
566 if (split_len < uncomp_len) {
567 if (split_len < se_len)
568 return m_encode_splitting(strm);
570 return m_encode_se(strm);
572 if (uncomp_len <= se_len)
573 return m_encode_uncomp(strm);
575 return m_encode_se(strm);
579 static int m_check_zero_block(struct aec_stream *strm)
582 Check if input block is all zero.
584 Aggregate consecutive zero blocks until we find !0 or reach the
585 end of a segment or RSI.
588 struct internal_state *state = strm->state;
589 uint32_t *p = state->block + state->ref;
590 uint32_t *end = state->block + strm->block_size;
592 while(p < end && *p == 0)
596 if (state->zero_blocks) {
597 /* The current block isn't zero but we have to emit a
598 * previous zero block first. The current block will be
601 state->block -= strm->block_size;
602 state->blocks_avail++;
603 state->mode = m_encode_zero;
606 state->mode = m_select_code_option;
609 state->zero_blocks++;
610 if (state->zero_blocks == 1) {
611 state->zero_ref = state->ref;
612 state->zero_ref_sample = state->block[0];
614 if (state->blocks_avail == 0
615 || (strm->rsi - state->blocks_avail) % 64 == 0) {
616 if (state->zero_blocks > 4)
617 state->zero_blocks = ROS;
618 state->mode = m_encode_zero;
621 state->mode = m_get_block;
626 static int m_get_rsi_resumable(struct aec_stream *strm)
629 Get RSI while input buffer is short.
631 Let user provide more input. Once we got all input pad buffer
636 struct internal_state *state = strm->state;
639 if (strm->avail_in > 0) {
640 state->data_raw[state->i] = state->get_sample(strm);
642 if (state->flush == AEC_FLUSH) {
644 for (j = state->i; j < strm->rsi * strm->block_size; j++)
645 state->data_raw[j] = state->data_raw[state->i - 1];
646 state->i = strm->rsi * strm->block_size;
648 if (state->zero_blocks) {
649 state->mode = m_encode_zero;
653 emit(state, 0, state->bits);
654 if (strm->avail_out > 0) {
655 if (!state->direct_out)
656 *strm->next_out++ = *state->cds;
666 } while (++state->i < strm->rsi * strm->block_size);
668 if (strm->flags & AEC_DATA_PREPROCESS)
669 state->preprocess(strm);
671 return m_check_zero_block(strm);
674 static int m_get_block(struct aec_stream *strm)
677 Provide the next block of preprocessed input data.
679 Pull in a whole Reference Sample Interval (RSI) of data if
680 block buffer is empty.
683 struct internal_state *state = strm->state;
685 if (strm->avail_out > state->cds_len) {
686 if (!state->direct_out) {
687 state->direct_out = 1;
688 *strm->next_out = *state->cds;
689 state->cds = strm->next_out;
692 if (state->zero_blocks == 0 || state->direct_out) {
693 /* copy leftover from last block */
694 *state->cds_buf = *state->cds;
695 state->cds = state->cds_buf;
697 state->direct_out = 0;
700 if (state->blocks_avail == 0) {
701 state->blocks_avail = strm->rsi - 1;
702 state->block = state->data_pp;
704 if (strm->avail_in >= state->rsi_len) {
705 state->get_rsi(strm);
706 if (strm->flags & AEC_DATA_PREPROCESS)
707 state->preprocess(strm);
709 return m_check_zero_block(strm);
712 state->mode = m_get_rsi_resumable;
716 state->block += strm->block_size;
717 state->blocks_avail--;
718 return m_check_zero_block(strm);
729 int aec_encode_init(struct aec_stream *strm)
731 struct internal_state *state;
733 if (strm->bits_per_sample > 32 || strm->bits_per_sample == 0)
734 return AEC_CONF_ERROR;
736 if (strm->block_size != 8
737 && strm->block_size != 16
738 && strm->block_size != 32
739 && strm->block_size != 64)
740 return AEC_CONF_ERROR;
742 if (strm->rsi > 4096)
743 return AEC_CONF_ERROR;
745 state = malloc(sizeof(struct internal_state));
747 return AEC_MEM_ERROR;
749 memset(state, 0, sizeof(struct internal_state));
752 if (strm->bits_per_sample > 16) {
753 /* 24/32 input bit settings */
756 if (strm->bits_per_sample <= 24
757 && strm->flags & AEC_DATA_3BYTE) {
759 if (strm->flags & AEC_DATA_MSB) {
760 state->get_sample = aec_get_msb_24;
761 state->get_rsi = aec_get_rsi_msb_24;
763 state->get_sample = aec_get_lsb_24;
764 state->get_rsi = aec_get_rsi_lsb_24;
768 if (strm->flags & AEC_DATA_MSB) {
769 state->get_sample = aec_get_msb_32;
770 state->get_rsi = aec_get_rsi_msb_32;
772 state->get_sample = aec_get_lsb_32;
773 state->get_rsi = aec_get_rsi_lsb_32;
777 else if (strm->bits_per_sample > 8) {
778 /* 16 bit settings */
782 if (strm->flags & AEC_DATA_MSB) {
783 state->get_sample = aec_get_msb_16;
784 state->get_rsi = aec_get_rsi_msb_16;
786 state->get_sample = aec_get_lsb_16;
787 state->get_rsi = aec_get_rsi_lsb_16;
794 state->get_sample = aec_get_8;
795 state->get_rsi = aec_get_rsi_8;
797 state->rsi_len *= strm->rsi * strm->block_size;
799 if (strm->flags & AEC_DATA_SIGNED) {
800 state->xmin = -(1ULL << (strm->bits_per_sample - 1));
801 state->xmax = (1ULL << (strm->bits_per_sample - 1)) - 1;
802 state->preprocess = preprocess_signed;
805 state->xmax = (1ULL << strm->bits_per_sample) - 1;
806 state->preprocess = preprocess_unsigned;
809 state->kmax = (1U << state->id_len) - 3;
811 state->data_pp = malloc(strm->rsi
814 if (state->data_pp == NULL)
815 return AEC_MEM_ERROR;
817 if (strm->flags & AEC_DATA_PREPROCESS) {
818 state->data_raw = malloc(strm->rsi
821 if (state->data_raw == NULL)
822 return AEC_MEM_ERROR;
824 state->data_raw = state->data_pp;
827 state->block = state->data_pp;
829 /* Largest possible CDS according to specs */
830 state->cds_len = (5 + 64 * 32) / 8 + 3;
831 state->cds_buf = malloc(state->cds_len);
832 if (state->cds_buf == NULL)
833 return AEC_MEM_ERROR;
838 state->cds = state->cds_buf;
841 state->mode = m_get_block;
846 int aec_encode(struct aec_stream *strm, int flush)
849 Finite-state machine implementation of the adaptive entropy
853 struct internal_state *state = strm->state;
855 state->flush = flush;
857 while (state->mode(strm) == M_CONTINUE);
859 if (state->direct_out) {
860 n = state->cds - strm->next_out;
862 strm->avail_out -= n;
863 strm->total_out += n;
865 *state->cds_buf = *state->cds;
866 state->cds = state->cds_buf;
867 state->direct_out = 0;
872 int aec_encode_end(struct aec_stream *strm)
874 struct internal_state *state = strm->state;
876 if (strm->flags & AEC_DATA_PREPROCESS)
877 free(state->data_raw);
878 free(state->data_pp);
879 free(state->cds_buf);
884 int aec_buffer_encode(struct aec_stream *strm)
888 status = aec_encode_init(strm);
889 if (status != AEC_OK)
891 status = aec_encode(strm, AEC_FLUSH);
892 if (strm->avail_in > 0 || strm->avail_out == 0)
893 status = AEC_DATA_ERROR;
895 aec_encode_end(strm);