Line data Source code
1 : /*
2 : Copyright 2015 Esri
3 :
4 : Licensed under the Apache License, Version 2.0 (the "License");
5 : you may not use this file except in compliance with the License.
6 : You may obtain a copy of the License at
7 :
8 : http://www.apache.org/licenses/LICENSE-2.0
9 :
10 : Unless required by applicable law or agreed to in writing, software
11 : distributed under the License is distributed on an "AS IS" BASIS,
12 : WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 : See the License for the specific language governing permissions and
14 : limitations under the License.
15 :
16 : A local copy of the license and additional notices are located with the
17 : source distribution at:
18 :
19 : http://github.com/Esri/lerc/
20 :
21 : Contributors: Thomas Maurer
22 : */
23 :
24 : #include <algorithm>
25 : #include <queue>
26 : #include "Defines.h"
27 : #include "Huffman.h"
28 : #include "BitStuffer2.h"
29 :
30 : using namespace std;
31 : USING_NAMESPACE_LERC
32 :
33 : // -------------------------------------------------------------------------- ;
34 :
35 3531 : bool Huffman::ComputeCodes(const vector<int>& histo)
36 : {
37 3531 : if (histo.empty() || histo.size() >= m_maxHistoSize)
38 0 : return false;
39 :
40 7062 : priority_queue<Node, vector<Node>, less<Node> > pq;
41 :
42 3531 : int numNodes = 0;
43 :
44 3531 : int size = (int)histo.size();
45 907467 : for (int i = 0; i < size; i++) // add all leaf nodes
46 903936 : if (histo[i] > 0)
47 124197 : pq.push(Node((short)i, histo[i]));
48 :
49 3531 : if (pq.size() < 2) // histo has only 0 or 1 bin that is not empty; quit Huffman and give it to Lerc
50 0 : return false;
51 :
52 124197 : while (pq.size() > 1) // build the Huffman tree
53 : {
54 120666 : Node* child0 = new Node(pq.top());
55 120666 : numNodes++;
56 120666 : pq.pop();
57 120666 : Node* child1 = new Node(pq.top());
58 120666 : numNodes++;
59 120666 : pq.pop();
60 120666 : pq.push(Node(child0, child1));
61 : }
62 :
63 3531 : m_codeTable.resize(size);
64 3531 : std::fill(m_codeTable.begin(), m_codeTable.end(),
65 3531 : std::pair<unsigned short, unsigned int>((short)0, 0));
66 :
67 3531 : if (!pq.top().TreeToLUT(0, 0, m_codeTable)) // fill the LUT
68 0 : return false;
69 :
70 : //pq.top().FreeTree(numNodes); // Linux compiler complains
71 3531 : Node nodeNonConst = pq.top();
72 3531 : nodeNonConst.FreeTree(numNodes); // free all the nodes
73 :
74 3531 : if (numNodes != 0) // check the ref count
75 0 : return false;
76 :
77 3531 : if (!ConvertCodesToCanonical())
78 0 : return false;
79 :
80 3531 : return true;
81 : }
82 :
83 : // -------------------------------------------------------------------------- ;
84 :
85 3531 : bool Huffman::ComputeCompressedSize(const std::vector<int>& histo, int& numBytes, double& avgBpp) const
86 : {
87 3531 : if (histo.empty() || histo.size() >= m_maxHistoSize)
88 0 : return false;
89 :
90 3531 : numBytes = 0;
91 3531 : if (!ComputeNumBytesCodeTable(numBytes)) // header and code table
92 0 : return false;
93 :
94 3531 : int numBits = 0, numElem = 0;
95 3531 : int size = (int)histo.size();
96 907467 : for (int i = 0; i < size; i++)
97 903936 : if (histo[i] > 0)
98 : {
99 124197 : numBits += histo[i] * m_codeTable[i].first;
100 124197 : numElem += histo[i];
101 : }
102 :
103 3531 : if (numElem == 0)
104 0 : return false;
105 :
106 3531 : int numUInts = ((((numBits + 7) >> 3) + 3) >> 2) + 1; // add one more as the decode LUT can read ahead
107 3531 : numBytes += 4 * numUInts; // data huffman coded
108 3531 : avgBpp = 8 * numBytes / (double)numElem;
109 :
110 3531 : return true;
111 : }
112 :
113 : // -------------------------------------------------------------------------- ;
114 :
115 640 : bool Huffman::SetCodes(const vector<pair<unsigned short, unsigned int> >& codeTable)
116 : {
117 640 : if (codeTable.empty() || codeTable.size() >= m_maxHistoSize)
118 0 : return false;
119 :
120 640 : m_codeTable = codeTable;
121 640 : return true;
122 : }
123 :
124 : // -------------------------------------------------------------------------- ;
125 :
126 640 : bool Huffman::WriteCodeTable(Byte** ppByte, int lerc2Version) const
127 : {
128 640 : if (!ppByte)
129 0 : return false;
130 :
131 : int i0, i1, maxLen;
132 640 : if (!GetRange(i0, i1, maxLen))
133 0 : return false;
134 :
135 640 : int size = (int)m_codeTable.size();
136 1280 : vector<unsigned int> dataVec(i1 - i0, 0);
137 :
138 96762 : for (int i = i0; i < i1; i++)
139 : {
140 96122 : int k = GetIndexWrapAround(i, size);
141 96122 : dataVec[i - i0] = m_codeTable[k].first;
142 : }
143 :
144 : // header
145 1280 : vector<int> intVec;
146 640 : intVec.push_back(4); // huffman version; 4 guarantees canonical codes
147 640 : intVec.push_back(size);
148 640 : intVec.push_back(i0); // code range
149 640 : intVec.push_back(i1);
150 :
151 640 : Byte* ptr = *ppByte;
152 :
153 640 : size_t len = intVec.size() * sizeof(int);
154 640 : memcpy(ptr, &intVec[0], len);
155 640 : ptr += len;
156 :
157 1280 : BitStuffer2 bitStuffer2;
158 640 : if (!bitStuffer2.EncodeSimple(&ptr, dataVec, lerc2Version)) // code lengths, bit stuffed
159 0 : return false;
160 :
161 640 : if (!BitStuffCodes(&ptr, i0, i1)) // variable length codes, bit stuffed
162 0 : return false;
163 :
164 640 : *ppByte = ptr;
165 640 : return true;
166 : }
167 :
168 : // -------------------------------------------------------------------------- ;
169 :
170 409 : bool Huffman::ReadCodeTable(const Byte** ppByte, size_t& nBytesRemainingInOut, int lerc2Version)
171 : {
172 409 : if (!ppByte || !(*ppByte))
173 0 : return false;
174 :
175 409 : const Byte* ptr = *ppByte;
176 409 : size_t nBytesRemaining = nBytesRemainingInOut;
177 :
178 818 : vector<int> intVec(4, 0);
179 409 : size_t len = intVec.size() * sizeof(int);
180 :
181 409 : if (nBytesRemaining < len)
182 0 : return false;
183 :
184 409 : memcpy(&intVec[0], ptr, len);
185 409 : ptr += len;
186 409 : nBytesRemaining -= len;
187 :
188 409 : int version = intVec[0];
189 :
190 409 : if (version < 2) // allow forward compatibility; for updates that break old decoders increase Lerc2 version number;
191 0 : return false;
192 :
193 409 : const int size = intVec[1];
194 409 : const int i0 = intVec[2];
195 409 : const int i1 = intVec[3];
196 :
197 409 : if (i0 >= i1 || i0 < 0 || size < 0 || size > (int)m_maxHistoSize)
198 0 : return false;
199 :
200 409 : if (GetIndexWrapAround(i0, size) >= size || GetIndexWrapAround(i1 - 1, size) >= size)
201 0 : return false;
202 :
203 : try
204 : {
205 818 : vector<unsigned int> dataVec(i1 - i0, 0);
206 818 : BitStuffer2 bitStuffer2;
207 409 : if (!bitStuffer2.Decode(&ptr, nBytesRemaining, dataVec, dataVec.size(), lerc2Version)) // unstuff the code lengths
208 0 : return false;
209 :
210 409 : if (dataVec.size() != static_cast<size_t>(i1 - i0))
211 0 : return false;
212 :
213 409 : m_codeTable.resize(size);
214 409 : std::fill(m_codeTable.begin(), m_codeTable.end(),
215 409 : std::pair<unsigned short, unsigned int>((short)0, 0));
216 :
217 75770 : for (int i = i0; i < i1; i++)
218 : {
219 75361 : int k = GetIndexWrapAround(i, size);
220 75361 : m_codeTable[k].first = (unsigned short)dataVec[i - i0];
221 : }
222 :
223 409 : if (!BitUnStuffCodes(&ptr, nBytesRemaining, i0, i1)) // unstuff the codes
224 0 : return false;
225 :
226 409 : *ppByte = ptr;
227 409 : nBytesRemainingInOut = nBytesRemaining;
228 409 : return true;
229 : }
230 0 : catch (std::exception&)
231 : {
232 0 : return false;
233 : }
234 : }
235 :
236 : // -------------------------------------------------------------------------- ;
237 :
238 409 : bool Huffman::BuildTreeFromCodes(int& numBitsLUT)
239 : {
240 409 : int i0 = 0, i1 = 0, maxLen = 0;
241 409 : if (!GetRange(i0, i1, maxLen))
242 0 : return false;
243 :
244 : // build decode LUT using max of 12 bits
245 409 : int size = (int)m_codeTable.size();
246 409 : int minNumZeroBits = 32;
247 :
248 409 : bool bNeedTree = maxLen > m_maxNumBitsLUT;
249 409 : numBitsLUT = min(maxLen, m_maxNumBitsLUT);
250 :
251 409 : int sizeLUT = 1 << numBitsLUT;
252 :
253 409 : m_decodeLUT.clear();
254 409 : m_decodeLUT.assign((size_t)sizeLUT, pair<short, short>((short)-1, (short)-1));
255 :
256 75770 : for (int i = i0; i < i1; i++)
257 : {
258 75361 : int k = GetIndexWrapAround(i, size);
259 75361 : int len = m_codeTable[k].first;
260 :
261 75361 : if (len == 0)
262 56096 : continue;
263 :
264 19265 : unsigned int code = m_codeTable[k].second;
265 :
266 19265 : if (len <= numBitsLUT)
267 : {
268 18376 : code <<= (numBitsLUT - len);
269 18376 : unsigned int numEntries = 1 << (numBitsLUT - len);
270 :
271 472641 : for (unsigned int j = 0; j < numEntries; j++)
272 : {
273 454265 : auto& entry = m_decodeLUT[code | j];
274 454265 : entry.first = (short)len; // add the duplicates
275 454265 : entry.second = (short)k; // add the duplicates
276 : }
277 : }
278 : else // for the codes too long for the LUT, count how many leading bits are 0
279 : {
280 889 : int shift = 1;
281 3753 : while (code >>= 1) shift++; // large canonical codes start with zero's
282 889 : minNumZeroBits = min(minNumZeroBits, len - shift);
283 : }
284 : }
285 :
286 409 : m_numBitsToSkipInTree = bNeedTree? minNumZeroBits : 0;
287 :
288 409 : if (!bNeedTree) // decode LUT covers it all, no tree needed
289 377 : return true;
290 :
291 : //m_numBitsToSkipInTree = 0; // to disable skipping the 0 bits
292 :
293 32 : ClearTree(); // if there
294 :
295 32 : Node emptyNode((short)-1, 0);
296 32 : m_root = new Node(emptyNode);
297 :
298 7978 : for (int i = i0; i < i1; i++)
299 : {
300 7946 : int k = GetIndexWrapAround(i, size);
301 7946 : int len = m_codeTable[k].first;
302 :
303 7946 : if (len > 0 && len > numBitsLUT) // add only codes not in the decode LUT
304 : {
305 889 : unsigned int code = m_codeTable[k].second;
306 889 : Node* node = m_root;
307 889 : int j = len - m_numBitsToSkipInTree; // reduce len by number of leading 0 bits from above
308 :
309 5907 : while (--j >= 0) // go over the bits
310 : {
311 5018 : if (code & (1 << j))
312 : {
313 2156 : if (!node->child1)
314 857 : node->child1 = new Node(emptyNode);
315 :
316 2156 : node = node->child1;
317 : }
318 : else
319 : {
320 2862 : if (!node->child0)
321 898 : node->child0 = new Node(emptyNode);
322 :
323 2862 : node = node->child0;
324 : }
325 :
326 5018 : if (j == 0) // last bit, leaf node
327 889 : node->value = (short)k; // set the value
328 : }
329 : }
330 : }
331 :
332 32 : return true;
333 : }
334 :
335 : // -------------------------------------------------------------------------- ;
336 :
337 7703 : void Huffman::Clear()
338 : {
339 7703 : m_codeTable.clear();
340 7703 : m_decodeLUT.clear();
341 7703 : ClearTree();
342 7703 : }
343 :
344 : // -------------------------------------------------------------------------- ;
345 :
346 7735 : void Huffman::ClearTree()
347 : {
348 7735 : if (m_root)
349 : {
350 32 : int n = 0;
351 32 : m_root->FreeTree(n);
352 32 : delete m_root;
353 32 : m_root = nullptr;
354 : }
355 7735 : }
356 :
357 : // -------------------------------------------------------------------------- ;
358 : // -------------------------------------------------------------------------- ;
359 :
360 3531 : bool Huffman::ComputeNumBytesCodeTable(int& numBytes) const
361 : {
362 : int i0, i1, maxLen;
363 3531 : if (!GetRange(i0, i1, maxLen))
364 0 : return false;
365 :
366 3531 : int size = (int)m_codeTable.size();
367 3531 : int sum = 0;
368 464779 : for (int i = i0; i < i1; i++)
369 : {
370 461248 : int k = GetIndexWrapAround(i, size);
371 461248 : sum += m_codeTable[k].first;
372 : }
373 :
374 3531 : numBytes = 4 * sizeof(int); // version, size, first bin, (last + 1) bin
375 :
376 3531 : BitStuffer2 bitStuffer2;
377 3531 : numBytes += bitStuffer2.ComputeNumBytesNeededSimple((unsigned int)(i1 - i0), (unsigned int)maxLen); // code lengths
378 3531 : int numUInts = (((sum + 7) >> 3) + 3) >> 2;
379 3531 : numBytes += 4 * numUInts; // byte array with the codes bit stuffed
380 :
381 3531 : return true;
382 : }
383 :
384 : // -------------------------------------------------------------------------- ;
385 :
386 4580 : bool Huffman::GetRange(int& i0, int& i1, int& maxCodeLength) const
387 : {
388 4580 : if (m_codeTable.empty() || m_codeTable.size() >= m_maxHistoSize)
389 0 : return false;
390 :
391 : // first, check for peak somewhere in the middle with 0 stretches left and right
392 4580 : int size = (int)m_codeTable.size();
393 : {
394 4580 : int i = 0;
395 7421 : while (i < size && m_codeTable[i].first == 0) i++;
396 4580 : i0 = i;
397 4580 : i = size - 1;
398 28503 : while (i >= 0 && m_codeTable[i].first == 0) i--;
399 4580 : i1 = i + 1; // exclusive
400 : }
401 :
402 4580 : if (i1 <= i0)
403 0 : return false;
404 :
405 : // second, cover the common case that the peak is close to 0
406 4580 : pair<int, int> segm(0, 0);
407 4580 : int j = 0;
408 82745 : while (j < size) // find the largest stretch of 0's, if any
409 : {
410 249370 : while (j < size && m_codeTable[j].first > 0) j++;
411 78165 : int k0 = j;
412 1079440 : while (j < size && m_codeTable[j].first == 0) j++;
413 78165 : int k1 = j;
414 :
415 78165 : if (k1 - k0 > segm.second)
416 17106 : segm = pair<int, int>(k0, k1 - k0);
417 : }
418 :
419 4580 : if (size - segm.second < i1 - i0)
420 : {
421 4249 : i0 = segm.first + segm.second;
422 4249 : i1 = segm.first + size; // do wrap around
423 : }
424 :
425 4580 : if (i1 <= i0)
426 0 : return false;
427 :
428 4580 : int maxLen = 0;
429 637311 : for (int i = i0; i < i1; i++)
430 : {
431 632731 : int k = GetIndexWrapAround(i, size);
432 632731 : int len = m_codeTable[k].first;
433 632731 : maxLen = max(maxLen, len);
434 : }
435 :
436 4580 : if (maxLen <= 0 || maxLen > 32)
437 0 : return false;
438 :
439 4580 : maxCodeLength = maxLen;
440 4580 : return true;
441 : }
442 :
443 : // -------------------------------------------------------------------------- ;
444 :
445 640 : bool Huffman::BitStuffCodes(Byte** ppByte, int i0, int i1) const
446 : {
447 640 : if (!ppByte)
448 0 : return false;
449 :
450 640 : unsigned int* arr = (unsigned int*)(*ppByte);
451 640 : unsigned int* dstPtr = arr;
452 640 : int size = (int)m_codeTable.size();
453 640 : int bitPos = 0;
454 :
455 96762 : for (int i = i0; i < i1; i++)
456 : {
457 96122 : int k = GetIndexWrapAround(i, size);
458 96122 : int len = m_codeTable[k].first;
459 96122 : if (len > 0)
460 : {
461 27743 : unsigned int val = m_codeTable[k].second;
462 27743 : if (32 - bitPos >= len)
463 : {
464 22449 : if (bitPos == 0)
465 1382 : *dstPtr = 0;
466 :
467 22449 : *dstPtr |= val << (32 - bitPos - len);
468 22449 : bitPos += len;
469 22449 : if (bitPos == 32)
470 : {
471 760 : bitPos = 0;
472 760 : dstPtr++;
473 : }
474 : }
475 : else
476 : {
477 5294 : bitPos += len - 32;
478 5294 : *dstPtr++ |= val >> bitPos; // bitPos > 0
479 5294 : *dstPtr = val << (32 - bitPos);
480 : }
481 : }
482 : }
483 :
484 640 : size_t numUInts = dstPtr - arr + (bitPos > 0 ? 1 : 0);
485 640 : *ppByte += numUInts * sizeof(unsigned int);
486 640 : return true;
487 : }
488 :
489 : // -------------------------------------------------------------------------- ;
490 :
491 409 : bool Huffman::BitUnStuffCodes(const Byte** ppByte, size_t& nBytesRemainingInOut, int i0, int i1)
492 : {
493 409 : if (!ppByte || !(*ppByte))
494 0 : return false;
495 :
496 409 : size_t nBytesRemaining = nBytesRemainingInOut;
497 :
498 409 : const unsigned int* arr = (const unsigned int*)(*ppByte);
499 409 : const unsigned int* srcPtr = arr;
500 409 : const size_t sizeUInt = sizeof(*srcPtr);
501 :
502 409 : int size = (int)m_codeTable.size();
503 409 : int bitPos = 0;
504 :
505 75770 : for (int i = i0; i < i1; i++)
506 : {
507 75361 : int k = GetIndexWrapAround(i, size);
508 75361 : int len = m_codeTable[k].first;
509 75361 : if (len > 0)
510 : {
511 19265 : if (nBytesRemaining < sizeUInt || len > 32)
512 0 : return false;
513 :
514 19265 : m_codeTable[k].second = ((*srcPtr) << bitPos) >> (32 - len);
515 :
516 19265 : if (32 - bitPos >= len)
517 : {
518 15344 : bitPos += len;
519 15344 : if (bitPos == 32)
520 : {
521 480 : bitPos = 0;
522 480 : srcPtr++;
523 480 : nBytesRemaining -= sizeUInt;
524 : }
525 : }
526 : else
527 : {
528 3921 : bitPos += len - 32;
529 3921 : srcPtr++;
530 3921 : nBytesRemaining -= sizeUInt;
531 :
532 3921 : if (nBytesRemaining < sizeUInt)
533 0 : return false;
534 :
535 3921 : m_codeTable[k].second |= (*srcPtr) >> (32 - bitPos); // bitPos > 0
536 : }
537 : }
538 : }
539 :
540 409 : size_t numUInts = srcPtr - arr + (bitPos > 0 ? 1 : 0);
541 409 : size_t len = numUInts * sizeUInt;
542 :
543 409 : if (nBytesRemainingInOut < len)
544 0 : return false;
545 :
546 409 : *ppByte += len;
547 409 : nBytesRemainingInOut -= len;
548 :
549 409 : if (nBytesRemaining != nBytesRemainingInOut && nBytesRemaining != nBytesRemainingInOut + sizeUInt) // the real check
550 0 : return false;
551 :
552 409 : return true;
553 : }
554 :
555 : // -------------------------------------------------------------------------- ;
556 :
557 : //struct MyLargerThanOp
558 : //{
559 : // inline bool operator() (const pair<int, unsigned int>& p0,
560 : // const pair<int, unsigned int>& p1) { return p0.first > p1.first; }
561 : //};
562 :
563 : // -------------------------------------------------------------------------- ;
564 :
565 3531 : bool Huffman::ConvertCodesToCanonical()
566 : {
567 : // from the non canonical code book, create an array to be sorted in descending order:
568 : // codeLength * tableSize - index
569 :
570 3531 : unsigned int tableSize = (unsigned int)m_codeTable.size();
571 3531 : if (tableSize == 0)
572 0 : return true;
573 3531 : vector<pair<int, unsigned int> > sortVec(tableSize, pair<int, unsigned int>(0, 0));
574 : //memset(&sortVec[0], 0, tableSize * sizeof(pair<int, unsigned int>));
575 :
576 907467 : for (unsigned int i = 0; i < tableSize; i++)
577 903936 : if (m_codeTable[i].first > 0)
578 124197 : sortVec[i] = pair<int, unsigned int>(m_codeTable[i].first * tableSize - i, i);
579 :
580 : // sort descending
581 : //std::sort(sortVec.begin(), sortVec.end(), MyLargerThanOp());
582 :
583 3531 : std::sort(sortVec.begin(), sortVec.end(),
584 5734690 : [](const pair<int, unsigned int>& p0,
585 5734690 : const pair<int, unsigned int>& p1) { return p0.first > p1.first; });
586 :
587 : // create canonical codes and assign to orig code table
588 3531 : unsigned int index = sortVec[0].second;
589 3531 : unsigned short codeLen = m_codeTable[index].first; // max code length for this table
590 3531 : unsigned int i = 0, codeCanonical = 0;
591 :
592 127728 : while (i < tableSize && sortVec[i].first > 0)
593 : {
594 124197 : index = sortVec[i++].second;
595 124197 : short delta = codeLen - m_codeTable[index].first; // difference of 2 consecutive code lengths, >= 0 as sorted
596 124197 : codeCanonical >>= delta;
597 124197 : codeLen -= delta;
598 124197 : m_codeTable[index].second = codeCanonical++;
599 : }
600 :
601 3531 : return true;
602 : }
603 :
604 : // -------------------------------------------------------------------------- ;
|
