| 1 | /* license:BSD-3-Clause |
| 2 | * copyright-holders:Aaron Giles |
| 3 | **************************************************************************** |
| 4 | |
| 5 | huffman.c |
| 6 | |
| 7 | Static Huffman compression and decompression helpers. |
| 8 | |
| 9 | **************************************************************************** |
| 10 | |
| 11 | Maximum codelength is officially (alphabetsize - 1). This would be 255 bits |
| 12 | (since we use 1 byte values). However, it is also dependent upon the number |
| 13 | of samples used, as follows: |
| 14 | |
| 15 | 2 bits -> 3..4 samples |
| 16 | 3 bits -> 5..7 samples |
| 17 | 4 bits -> 8..12 samples |
| 18 | 5 bits -> 13..20 samples |
| 19 | 6 bits -> 21..33 samples |
| 20 | 7 bits -> 34..54 samples |
| 21 | 8 bits -> 55..88 samples |
| 22 | 9 bits -> 89..143 samples |
| 23 | 10 bits -> 144..232 samples |
| 24 | 11 bits -> 233..376 samples |
| 25 | 12 bits -> 377..609 samples |
| 26 | 13 bits -> 610..986 samples |
| 27 | 14 bits -> 987..1596 samples |
| 28 | 15 bits -> 1597..2583 samples |
| 29 | 16 bits -> 2584..4180 samples -> note that a 4k data size guarantees codelength <= 16 bits |
| 30 | 17 bits -> 4181..6764 samples |
| 31 | 18 bits -> 6765..10945 samples |
| 32 | 19 bits -> 10946..17710 samples |
| 33 | 20 bits -> 17711..28656 samples |
| 34 | 21 bits -> 28657..46367 samples |
| 35 | 22 bits -> 46368..75024 samples |
| 36 | 23 bits -> 75025..121392 samples |
| 37 | 24 bits -> 121393..196417 samples |
| 38 | 25 bits -> 196418..317810 samples |
| 39 | 26 bits -> 317811..514228 samples |
| 40 | 27 bits -> 514229..832039 samples |
| 41 | 28 bits -> 832040..1346268 samples |
| 42 | 29 bits -> 1346269..2178308 samples |
| 43 | 30 bits -> 2178309..3524577 samples |
| 44 | 31 bits -> 3524578..5702886 samples |
| 45 | 32 bits -> 5702887..9227464 samples |
| 46 | |
| 47 | Looking at it differently, here is where powers of 2 fall into these buckets: |
| 48 | |
| 49 | 256 samples -> 11 bits max |
| 50 | 512 samples -> 12 bits max |
| 51 | 1k samples -> 14 bits max |
| 52 | 2k samples -> 15 bits max |
| 53 | 4k samples -> 16 bits max |
| 54 | 8k samples -> 18 bits max |
| 55 | 16k samples -> 19 bits max |
| 56 | 32k samples -> 21 bits max |
| 57 | 64k samples -> 22 bits max |
| 58 | 128k samples -> 24 bits max |
| 59 | 256k samples -> 25 bits max |
| 60 | 512k samples -> 27 bits max |
| 61 | 1M samples -> 28 bits max |
| 62 | 2M samples -> 29 bits max |
| 63 | 4M samples -> 31 bits max |
| 64 | 8M samples -> 32 bits max |
| 65 | |
| 66 | **************************************************************************** |
| 67 | |
| 68 | Delta-RLE encoding works as follows: |
| 69 | |
| 70 | Starting value is assumed to be 0. All data is encoded as a delta |
| 71 | from the previous value, such that final[i] = final[i - 1] + delta. |
| 72 | Long runs of 0s are RLE-encoded as follows: |
| 73 | |
| 74 | 0x100 = repeat count of 8 |
| 75 | 0x101 = repeat count of 9 |
| 76 | 0x102 = repeat count of 10 |
| 77 | 0x103 = repeat count of 11 |
| 78 | 0x104 = repeat count of 12 |
| 79 | 0x105 = repeat count of 13 |
| 80 | 0x106 = repeat count of 14 |
| 81 | 0x107 = repeat count of 15 |
| 82 | 0x108 = repeat count of 16 |
| 83 | 0x109 = repeat count of 32 |
| 84 | 0x10a = repeat count of 64 |
| 85 | 0x10b = repeat count of 128 |
| 86 | 0x10c = repeat count of 256 |
| 87 | 0x10d = repeat count of 512 |
| 88 | 0x10e = repeat count of 1024 |
| 89 | 0x10f = repeat count of 2048 |
| 90 | |
| 91 | Note that repeat counts are reset at the end of a row, so if a 0 run |
| 92 | extends to the end of a row, a large repeat count may be used. |
| 93 | |
| 94 | The reason for starting the run counts at 8 is that 0 is expected to |
| 95 | be the most common symbol, and is typically encoded in 1 or 2 bits. |
| 96 | |
| 97 | ***************************************************************************/ |
| 98 | |
| 99 | #include <stdlib.h> |
| 100 | #include <assert.h> |
| 101 | #include <stdio.h> |
| 102 | #include <string.h> |
| 103 | |
| 104 | #include "huffman.h" |
| 105 | |
| 106 | #define MAX(x,y) ((x) > (y) ? (x) : (y)) |
| 107 | |
| 108 | /*************************************************************************** |
| 109 | * MACROS |
| 110 | *************************************************************************** |
| 111 | */ |
| 112 | |
| 113 | #define MAKE_LOOKUP(code,bits) (((code) << 5) | ((bits) & 0x1f)) |
| 114 | |
| 115 | /*************************************************************************** |
| 116 | * IMPLEMENTATION |
| 117 | *************************************************************************** |
| 118 | */ |
| 119 | |
| 120 | /*------------------------------------------------- |
| 121 | * huffman_context_base - create an encoding/ |
| 122 | * decoding context |
| 123 | *------------------------------------------------- |
| 124 | */ |
| 125 | |
| 126 | struct huffman_decoder* create_huffman_decoder(int numcodes, int maxbits) |
| 127 | { |
| 128 | /* limit to 24 bits */ |
| 129 | if (maxbits > 24) |
| 130 | return NULL; |
| 131 | |
| 132 | struct huffman_decoder* decoder = (struct huffman_decoder*)malloc(sizeof(struct huffman_decoder)); |
| 133 | decoder->numcodes = numcodes; |
| 134 | decoder->maxbits = maxbits; |
| 135 | decoder->lookup = (lookup_value*)malloc(sizeof(lookup_value) * (1 << maxbits)); |
| 136 | decoder->huffnode = (struct node_t*)malloc(sizeof(struct node_t) * numcodes); |
| 137 | decoder->datahisto = NULL; |
| 138 | decoder->prevdata = 0; |
| 139 | decoder->rleremaining = 0; |
| 140 | return decoder; |
| 141 | } |
| 142 | |
| 143 | /*------------------------------------------------- |
| 144 | * decode_one - decode a single code from the |
| 145 | * huffman stream |
| 146 | *------------------------------------------------- |
| 147 | */ |
| 148 | |
| 149 | uint32_t huffman_decode_one(struct huffman_decoder* decoder, struct bitstream* bitbuf) |
| 150 | { |
| 151 | /* peek ahead to get maxbits worth of data */ |
| 152 | uint32_t bits = bitstream_peek(bitbuf, decoder->maxbits); |
| 153 | |
| 154 | /* look it up, then remove the actual number of bits for this code */ |
| 155 | lookup_value lookup = decoder->lookup[bits]; |
| 156 | bitstream_remove(bitbuf, lookup & 0x1f); |
| 157 | |
| 158 | /* return the value */ |
| 159 | return lookup >> 5; |
| 160 | } |
| 161 | |
| 162 | /*------------------------------------------------- |
| 163 | * import_tree_rle - import an RLE-encoded |
| 164 | * huffman tree from a source data stream |
| 165 | *------------------------------------------------- |
| 166 | */ |
| 167 | |
| 168 | enum huffman_error huffman_import_tree_rle(struct huffman_decoder* decoder, struct bitstream* bitbuf) |
| 169 | { |
| 170 | /* bits per entry depends on the maxbits */ |
| 171 | int numbits; |
| 172 | if (decoder->maxbits >= 16) |
| 173 | numbits = 5; |
| 174 | else if (decoder->maxbits >= 8) |
| 175 | numbits = 4; |
| 176 | else |
| 177 | numbits = 3; |
| 178 | |
| 179 | /* loop until we read all the nodes */ |
| 180 | int curnode; |
| 181 | for (curnode = 0; curnode < decoder->numcodes; ) |
| 182 | { |
| 183 | /* a non-one value is just raw */ |
| 184 | int nodebits = bitstream_read(bitbuf, numbits); |
| 185 | if (nodebits != 1) |
| 186 | decoder->huffnode[curnode++].numbits = nodebits; |
| 187 | |
| 188 | /* a one value is an escape code */ |
| 189 | else |
| 190 | { |
| 191 | /* a double 1 is just a single 1 */ |
| 192 | nodebits = bitstream_read(bitbuf, numbits); |
| 193 | if (nodebits == 1) |
| 194 | decoder->huffnode[curnode++].numbits = nodebits; |
| 195 | |
| 196 | /* otherwise, we need one for value for the repeat count */ |
| 197 | else |
| 198 | { |
| 199 | int repcount = bitstream_read(bitbuf, numbits) + 3; |
| 200 | while (repcount--) |
| 201 | decoder->huffnode[curnode++].numbits = nodebits; |
| 202 | } |
| 203 | } |
| 204 | } |
| 205 | |
| 206 | /* make sure we ended up with the right number */ |
| 207 | if (curnode != decoder->numcodes) |
| 208 | return HUFFERR_INVALID_DATA; |
| 209 | |
| 210 | /* assign canonical codes for all nodes based on their code lengths */ |
| 211 | enum huffman_error error = huffman_assign_canonical_codes(decoder); |
| 212 | if (error != HUFFERR_NONE) |
| 213 | return error; |
| 214 | |
| 215 | /* build the lookup table */ |
| 216 | huffman_build_lookup_table(decoder); |
| 217 | |
| 218 | /* determine final input length and report errors */ |
| 219 | return bitstream_overflow(bitbuf) ? HUFFERR_INPUT_BUFFER_TOO_SMALL : HUFFERR_NONE; |
| 220 | } |
| 221 | |
| 222 | |
| 223 | /*------------------------------------------------- |
| 224 | * import_tree_huffman - import a huffman-encoded |
| 225 | * huffman tree from a source data stream |
| 226 | *------------------------------------------------- |
| 227 | */ |
| 228 | |
| 229 | enum huffman_error huffman_import_tree_huffman(struct huffman_decoder* decoder, struct bitstream* bitbuf) |
| 230 | { |
| 231 | /* start by parsing the lengths for the small tree */ |
| 232 | struct huffman_decoder* smallhuff = create_huffman_decoder(24, 6); |
| 233 | smallhuff->huffnode[0].numbits = bitstream_read(bitbuf, 3); |
| 234 | int start = bitstream_read(bitbuf, 3) + 1; |
| 235 | int count = 0; |
| 236 | for (int index = 1; index < 24; index++) |
| 237 | { |
| 238 | if (index < start || count == 7) |
| 239 | smallhuff->huffnode[index].numbits = 0; |
| 240 | else |
| 241 | { |
| 242 | count = bitstream_read(bitbuf, 3); |
| 243 | smallhuff->huffnode[index].numbits = (count == 7) ? 0 : count; |
| 244 | } |
| 245 | } |
| 246 | |
| 247 | /* then regenerate the tree */ |
| 248 | enum huffman_error error = huffman_assign_canonical_codes(smallhuff); |
| 249 | if (error != HUFFERR_NONE) |
| 250 | return error; |
| 251 | huffman_build_lookup_table(smallhuff); |
| 252 | |
| 253 | /* determine the maximum length of an RLE count */ |
| 254 | uint32_t temp = decoder->numcodes - 9; |
| 255 | uint8_t rlefullbits = 0; |
| 256 | while (temp != 0) |
| 257 | temp >>= 1, rlefullbits++; |
| 258 | |
| 259 | /* now process the rest of the data */ |
| 260 | int last = 0; |
| 261 | int curcode; |
| 262 | for (curcode = 0; curcode < decoder->numcodes; ) |
| 263 | { |
| 264 | int value = huffman_decode_one(smallhuff, bitbuf); |
| 265 | if (value != 0) |
| 266 | decoder->huffnode[curcode++].numbits = last = value - 1; |
| 267 | else |
| 268 | { |
| 269 | int count = bitstream_read(bitbuf, 3) + 2; |
| 270 | if (count == 7+2) |
| 271 | count += bitstream_read(bitbuf, rlefullbits); |
| 272 | for ( ; count != 0 && curcode < decoder->numcodes; count--) |
| 273 | decoder->huffnode[curcode++].numbits = last; |
| 274 | } |
| 275 | } |
| 276 | |
| 277 | /* make sure we ended up with the right number */ |
| 278 | if (curcode != decoder->numcodes) |
| 279 | return HUFFERR_INVALID_DATA; |
| 280 | |
| 281 | /* assign canonical codes for all nodes based on their code lengths */ |
| 282 | error = huffman_assign_canonical_codes(decoder); |
| 283 | if (error != HUFFERR_NONE) |
| 284 | return error; |
| 285 | |
| 286 | /* build the lookup table */ |
| 287 | huffman_build_lookup_table(decoder); |
| 288 | |
| 289 | /* determine final input length and report errors */ |
| 290 | return bitstream_overflow(bitbuf) ? HUFFERR_INPUT_BUFFER_TOO_SMALL : HUFFERR_NONE; |
| 291 | } |
| 292 | |
| 293 | /*------------------------------------------------- |
| 294 | * compute_tree_from_histo - common backend for |
| 295 | * computing a tree based on the data histogram |
| 296 | *------------------------------------------------- |
| 297 | */ |
| 298 | |
| 299 | enum huffman_error huffman_compute_tree_from_histo(struct huffman_decoder* decoder) |
| 300 | { |
| 301 | /* compute the number of data items in the histogram */ |
| 302 | uint32_t sdatacount = 0; |
| 303 | for (int i = 0; i < decoder->numcodes; i++) |
| 304 | sdatacount += decoder->datahisto[i]; |
| 305 | |
| 306 | /* binary search to achieve the optimum encoding */ |
| 307 | uint32_t lowerweight = 0; |
| 308 | uint32_t upperweight = sdatacount * 2; |
| 309 | while (1) |
| 310 | { |
| 311 | /* build a tree using the current weight */ |
| 312 | uint32_t curweight = (upperweight + lowerweight) / 2; |
| 313 | int curmaxbits = huffman_build_tree(decoder, sdatacount, curweight); |
| 314 | |
| 315 | /* apply binary search here */ |
| 316 | if (curmaxbits <= decoder->maxbits) |
| 317 | { |
| 318 | lowerweight = curweight; |
| 319 | |
| 320 | /* early out if it worked with the raw weights, or if we're done searching */ |
| 321 | if (curweight == sdatacount || (upperweight - lowerweight) <= 1) |
| 322 | break; |
| 323 | } |
| 324 | else |
| 325 | upperweight = curweight; |
| 326 | } |
| 327 | |
| 328 | /* assign canonical codes for all nodes based on their code lengths */ |
| 329 | return huffman_assign_canonical_codes(decoder); |
| 330 | } |
| 331 | |
| 332 | /*************************************************************************** |
| 333 | * INTERNAL FUNCTIONS |
| 334 | *************************************************************************** |
| 335 | */ |
| 336 | |
| 337 | /*------------------------------------------------- |
| 338 | * tree_node_compare - compare two tree nodes |
| 339 | * by weight |
| 340 | *------------------------------------------------- |
| 341 | */ |
| 342 | |
| 343 | static int huffman_tree_node_compare(const void *item1, const void *item2) |
| 344 | { |
| 345 | const struct node_t *node1 = *(const struct node_t **)item1; |
| 346 | const struct node_t *node2 = *(const struct node_t **)item2; |
| 347 | if (node2->weight != node1->weight) |
| 348 | return node2->weight - node1->weight; |
| 349 | if (node2->bits - node1->bits == 0) |
| 350 | fprintf(stderr, "identical node sort keys, should not happen!\n"); |
| 351 | return (int)node1->bits - (int)node2->bits; |
| 352 | } |
| 353 | |
| 354 | /*------------------------------------------------- |
| 355 | * build_tree - build a huffman tree based on the |
| 356 | * data distribution |
| 357 | *------------------------------------------------- |
| 358 | */ |
| 359 | |
| 360 | int huffman_build_tree(struct huffman_decoder* decoder, uint32_t totaldata, uint32_t totalweight) |
| 361 | { |
| 362 | /* make a list of all non-zero nodes */ |
| 363 | struct node_t** list = (struct node_t**)malloc(sizeof(struct node_t*) * decoder->numcodes * 2); |
| 364 | int listitems = 0; |
| 365 | memset(decoder->huffnode, 0, decoder->numcodes * sizeof(decoder->huffnode[0])); |
| 366 | for (int curcode = 0; curcode < decoder->numcodes; curcode++) |
| 367 | if (decoder->datahisto[curcode] != 0) |
| 368 | { |
| 369 | list[listitems++] = &decoder->huffnode[curcode]; |
| 370 | decoder->huffnode[curcode].count = decoder->datahisto[curcode]; |
| 371 | decoder->huffnode[curcode].bits = curcode; |
| 372 | |
| 373 | /* scale the weight by the current effective length, ensuring we don't go to 0 */ |
| 374 | decoder->huffnode[curcode].weight = ((uint64_t)decoder->datahisto[curcode]) * ((uint64_t)totalweight) / ((uint64_t)totaldata); |
| 375 | if (decoder->huffnode[curcode].weight == 0) |
| 376 | decoder->huffnode[curcode].weight = 1; |
| 377 | } |
| 378 | |
| 379 | #if 0 |
| 380 | fprintf(stderr, "Pre-sort:\n"); |
| 381 | for (int i = 0; i < listitems; i++) { |
| 382 | fprintf(stderr, "weight: %d code: %d\n", list[i]->m_weight, list[i]->m_bits); |
| 383 | } |
| 384 | #endif |
| 385 | |
| 386 | /* sort the list by weight, largest weight first */ |
| 387 | qsort(&list[0], listitems, sizeof(list[0]), huffman_tree_node_compare); |
| 388 | |
| 389 | #if 0 |
| 390 | fprintf(stderr, "Post-sort:\n"); |
| 391 | for (int i = 0; i < listitems; i++) { |
| 392 | fprintf(stderr, "weight: %d code: %d\n", list[i]->m_weight, list[i]->m_bits); |
| 393 | } |
| 394 | fprintf(stderr, "===================\n"); |
| 395 | #endif |
| 396 | |
| 397 | /* now build the tree */ |
| 398 | int nextalloc = decoder->numcodes; |
| 399 | while (listitems > 1) |
| 400 | { |
| 401 | /* remove lowest two items */ |
| 402 | struct node_t* node1 = &(*list[--listitems]); |
| 403 | struct node_t* node0 = &(*list[--listitems]); |
| 404 | |
| 405 | /* create new node */ |
| 406 | struct node_t* newnode = &decoder->huffnode[nextalloc++]; |
| 407 | newnode->parent = NULL; |
| 408 | node0->parent = node1->parent = newnode; |
| 409 | newnode->weight = node0->weight + node1->weight; |
| 410 | |
| 411 | /* insert into list at appropriate location */ |
| 412 | int curitem; |
| 413 | for (curitem = 0; curitem < listitems; curitem++) |
| 414 | if (newnode->weight > list[curitem]->weight) |
| 415 | { |
| 416 | memmove(&list[curitem+1], &list[curitem], (listitems - curitem) * sizeof(list[0])); |
| 417 | break; |
| 418 | } |
| 419 | list[curitem] = newnode; |
| 420 | listitems++; |
| 421 | } |
| 422 | |
| 423 | /* compute the number of bits in each code, and fill in another histogram */ |
| 424 | int maxbits = 0; |
| 425 | for (int curcode = 0; curcode < decoder->numcodes; curcode++) |
| 426 | { |
| 427 | struct node_t* node = &decoder->huffnode[curcode]; |
| 428 | node->numbits = 0; |
| 429 | node->bits = 0; |
| 430 | |
| 431 | /* if we have a non-zero weight, compute the number of bits */ |
| 432 | if (node->weight > 0) |
| 433 | { |
| 434 | /* determine the number of bits for this node */ |
| 435 | for (struct node_t *curnode = node; curnode->parent != NULL; curnode = curnode->parent) |
| 436 | node->numbits++; |
| 437 | if (node->numbits == 0) |
| 438 | node->numbits = 1; |
| 439 | |
| 440 | /* keep track of the max */ |
| 441 | maxbits = MAX(maxbits, ((int)node->numbits)); |
| 442 | } |
| 443 | } |
| 444 | return maxbits; |
| 445 | } |
| 446 | |
| 447 | /*------------------------------------------------- |
| 448 | * assign_canonical_codes - assign canonical codes |
| 449 | * to all the nodes based on the number of bits |
| 450 | * in each |
| 451 | *------------------------------------------------- |
| 452 | */ |
| 453 | |
| 454 | enum huffman_error huffman_assign_canonical_codes(struct huffman_decoder* decoder) |
| 455 | { |
| 456 | /* build up a histogram of bit lengths */ |
| 457 | uint32_t bithisto[33] = { 0 }; |
| 458 | for (int curcode = 0; curcode < decoder->numcodes; curcode++) |
| 459 | { |
| 460 | struct node_t* node = &decoder->huffnode[curcode]; |
| 461 | if (node->numbits > decoder->maxbits) |
| 462 | return HUFFERR_INTERNAL_INCONSISTENCY; |
| 463 | if (node->numbits <= 32) |
| 464 | bithisto[node->numbits]++; |
| 465 | } |
| 466 | |
| 467 | /* for each code length, determine the starting code number */ |
| 468 | uint32_t curstart = 0; |
| 469 | for (int codelen = 32; codelen > 0; codelen--) |
| 470 | { |
| 471 | uint32_t nextstart = (curstart + bithisto[codelen]) >> 1; |
| 472 | if (codelen != 1 && nextstart * 2 != (curstart + bithisto[codelen])) |
| 473 | return HUFFERR_INTERNAL_INCONSISTENCY; |
| 474 | bithisto[codelen] = curstart; |
| 475 | curstart = nextstart; |
| 476 | } |
| 477 | |
| 478 | /* now assign canonical codes */ |
| 479 | for (int curcode = 0; curcode < decoder->numcodes; curcode++) |
| 480 | { |
| 481 | struct node_t* node = &decoder->huffnode[curcode]; |
| 482 | if (node->numbits > 0) |
| 483 | node->bits = bithisto[node->numbits]++; |
| 484 | } |
| 485 | return HUFFERR_NONE; |
| 486 | } |
| 487 | |
| 488 | /*------------------------------------------------- |
| 489 | * build_lookup_table - build a lookup table for |
| 490 | * fast decoding |
| 491 | *------------------------------------------------- |
| 492 | */ |
| 493 | |
| 494 | void huffman_build_lookup_table(struct huffman_decoder* decoder) |
| 495 | { |
| 496 | /* iterate over all codes */ |
| 497 | for (int curcode = 0; curcode < decoder->numcodes; curcode++) |
| 498 | { |
| 499 | /* process all nodes which have non-zero bits */ |
| 500 | struct node_t* node = &decoder->huffnode[curcode]; |
| 501 | if (node->numbits > 0) |
| 502 | { |
| 503 | /* set up the entry */ |
| 504 | lookup_value value = MAKE_LOOKUP(curcode, node->numbits); |
| 505 | |
| 506 | /* fill all matching entries */ |
| 507 | int shift = decoder->maxbits - node->numbits; |
| 508 | lookup_value *dest = &decoder->lookup[node->bits << shift]; |
| 509 | lookup_value *destend = &decoder->lookup[((node->bits + 1) << shift) - 1]; |
| 510 | while (dest <= destend) |
| 511 | *dest++ = value; |
| 512 | } |
| 513 | } |
| 514 | } |