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 3530 : int numNodes = 0;
43 :
44 3530 : int size = (int)histo.size();
45 907219 : for (int i = 0; i < size; i++) // add all leaf nodes
46 903688 : if (histo[i] > 0)
47 124171 : 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 124177 : while (pq.size() > 1) // build the Huffman tree
53 : {
54 120646 : Node* child0 = new Node(pq.top());
55 120657 : numNodes++;
56 120657 : pq.pop();
57 120655 : Node* child1 = new Node(pq.top());
58 120663 : numNodes++;
59 120663 : pq.pop();
60 120663 : pq.push(Node(child0, child1));
61 : }
62 :
63 3532 : m_codeTable.resize(size);
64 3531 : std::fill(m_codeTable.begin(), m_codeTable.end(),
65 3530 : 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 3530 : if (numNodes != 0) // check the ref count
75 0 : return false;
76 :
77 3530 : 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 907352 : for (int i = 0; i < size; i++)
97 903821 : if (histo[i] > 0)
98 : {
99 124143 : numBits += histo[i] * m_codeTable[i].first;
100 124139 : 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 639 : bool Huffman::SetCodes(const vector<pair<unsigned short, unsigned int> >& codeTable)
116 : {
117 639 : if (codeTable.empty() || codeTable.size() >= m_maxHistoSize)
118 0 : return false;
119 :
120 639 : 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 96741 : for (int i = i0; i < i1; i++)
139 : {
140 96101 : int k = GetIndexWrapAround(i, size);
141 96094 : 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 639 : intVec.push_back(i0); // code range
149 640 : intVec.push_back(i1);
150 :
151 639 : Byte* ptr = *ppByte;
152 :
153 639 : size_t len = intVec.size() * sizeof(int);
154 639 : memcpy(ptr, &intVec[0], len);
155 639 : ptr += len;
156 :
157 1279 : BitStuffer2 bitStuffer2;
158 639 : 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 515 : bool Huffman::ReadCodeTable(const Byte** ppByte, size_t& nBytesRemainingInOut, int lerc2Version)
171 : {
172 515 : if (!ppByte || !(*ppByte))
173 2 : return false;
174 :
175 513 : const Byte* ptr = *ppByte;
176 513 : size_t nBytesRemaining = nBytesRemainingInOut;
177 :
178 1029 : vector<int> intVec(4, 0);
179 516 : size_t len = intVec.size() * sizeof(int);
180 :
181 514 : if (nBytesRemaining < len)
182 0 : return false;
183 :
184 514 : memcpy(&intVec[0], ptr, len);
185 515 : ptr += len;
186 515 : nBytesRemaining -= len;
187 :
188 515 : int version = intVec[0];
189 :
190 513 : if (version < 2) // allow forward compatibility; for updates that break old decoders increase Lerc2 version number;
191 0 : return false;
192 :
193 513 : const int size = intVec[1];
194 515 : const int i0 = intVec[2];
195 513 : const int i1 = intVec[3];
196 :
197 515 : if (i0 >= i1 || i0 < 0 || size < 0 || size > (int)m_maxHistoSize)
198 0 : return false;
199 :
200 515 : if (GetIndexWrapAround(i0, size) >= size || GetIndexWrapAround(i1 - 1, size) >= size)
201 0 : return false;
202 :
203 : try
204 : {
205 1030 : vector<unsigned int> dataVec(i1 - i0, 0);
206 1030 : BitStuffer2 bitStuffer2;
207 516 : if (!bitStuffer2.Decode(&ptr, nBytesRemaining, dataVec, dataVec.size(), lerc2Version)) // unstuff the code lengths
208 0 : return false;
209 :
210 516 : if (dataVec.size() != static_cast<size_t>(i1 - i0))
211 0 : return false;
212 :
213 516 : m_codeTable.resize(size);
214 513 : std::fill(m_codeTable.begin(), m_codeTable.end(),
215 513 : std::pair<unsigned short, unsigned int>((short)0, 0));
216 :
217 96541 : for (int i = i0; i < i1; i++)
218 : {
219 96025 : int k = GetIndexWrapAround(i, size);
220 96011 : m_codeTable[k].first = (unsigned short)dataVec[i - i0];
221 : }
222 :
223 516 : if (!BitUnStuffCodes(&ptr, nBytesRemaining, i0, i1)) // unstuff the codes
224 0 : return false;
225 :
226 516 : *ppByte = ptr;
227 516 : nBytesRemainingInOut = nBytesRemaining;
228 516 : return true;
229 : }
230 0 : catch (std::exception&)
231 : {
232 0 : return false;
233 : }
234 : }
235 :
236 : // -------------------------------------------------------------------------- ;
237 :
238 516 : bool Huffman::BuildTreeFromCodes(int& numBitsLUT)
239 : {
240 516 : int i0 = 0, i1 = 0, maxLen = 0;
241 516 : if (!GetRange(i0, i1, maxLen))
242 0 : return false;
243 :
244 : // build decode LUT using max of 12 bits
245 516 : int size = (int)m_codeTable.size();
246 516 : int minNumZeroBits = 32;
247 :
248 516 : bool bNeedTree = maxLen > m_maxNumBitsLUT;
249 516 : numBitsLUT = min(maxLen, m_maxNumBitsLUT);
250 :
251 516 : int sizeLUT = 1 << numBitsLUT;
252 :
253 516 : m_decodeLUT.clear();
254 516 : m_decodeLUT.assign((size_t)sizeLUT, pair<short, short>((short)-1, (short)-1));
255 :
256 96801 : for (int i = i0; i < i1; i++)
257 : {
258 96285 : int k = GetIndexWrapAround(i, size);
259 96270 : int len = m_codeTable[k].first;
260 :
261 96296 : if (len == 0)
262 64462 : continue;
263 :
264 31834 : unsigned int code = m_codeTable[k].second;
265 :
266 31831 : if (len <= numBitsLUT)
267 : {
268 30642 : code <<= (numBitsLUT - len);
269 30642 : unsigned int numEntries = 1 << (numBitsLUT - len);
270 :
271 685212 : for (unsigned int j = 0; j < numEntries; j++)
272 : {
273 654625 : auto& entry = m_decodeLUT[code | j];
274 654570 : entry.first = (short)len; // add the duplicates
275 654570 : 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 1189 : int shift = 1;
281 4953 : while (code >>= 1) shift++; // large canonical codes start with zero's
282 1189 : minNumZeroBits = min(minNumZeroBits, len - shift);
283 : }
284 : }
285 :
286 516 : m_numBitsToSkipInTree = bNeedTree? minNumZeroBits : 0;
287 :
288 516 : if (!bNeedTree) // decode LUT covers it all, no tree needed
289 474 : return true;
290 :
291 : //m_numBitsToSkipInTree = 0; // to disable skipping the 0 bits
292 :
293 42 : ClearTree(); // if there
294 :
295 42 : Node emptyNode((short)-1, 0);
296 42 : m_root = new Node(emptyNode);
297 :
298 10428 : for (int i = i0; i < i1; i++)
299 : {
300 10386 : int k = GetIndexWrapAround(i, size);
301 10386 : int len = m_codeTable[k].first;
302 :
303 10386 : if (len > 0 && len > numBitsLUT) // add only codes not in the decode LUT
304 : {
305 1189 : unsigned int code = m_codeTable[k].second;
306 1189 : Node* node = m_root;
307 1189 : int j = len - m_numBitsToSkipInTree; // reduce len by number of leading 0 bits from above
308 :
309 7707 : while (--j >= 0) // go over the bits
310 : {
311 6518 : if (code & (1 << j))
312 : {
313 2866 : if (!node->child1)
314 1147 : node->child1 = new Node(emptyNode);
315 :
316 2866 : node = node->child1;
317 : }
318 : else
319 : {
320 3652 : if (!node->child0)
321 1198 : node->child0 = new Node(emptyNode);
322 :
323 3652 : node = node->child0;
324 : }
325 :
326 6518 : if (j == 0) // last bit, leaf node
327 1189 : node->value = (short)k; // set the value
328 : }
329 : }
330 : }
331 :
332 42 : return true;
333 : }
334 :
335 : // -------------------------------------------------------------------------- ;
336 :
337 7808 : void Huffman::Clear()
338 : {
339 7808 : m_codeTable.clear();
340 7809 : m_decodeLUT.clear();
341 7807 : ClearTree();
342 7807 : }
343 :
344 : // -------------------------------------------------------------------------- ;
345 :
346 7849 : void Huffman::ClearTree()
347 : {
348 7849 : if (m_root)
349 : {
350 42 : int n = 0;
351 42 : m_root->FreeTree(n);
352 42 : delete m_root;
353 42 : m_root = nullptr;
354 : }
355 7849 : }
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 464748 : for (int i = i0; i < i1; i++)
369 : {
370 461217 : int k = GetIndexWrapAround(i, size);
371 461202 : 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 4687 : bool Huffman::GetRange(int& i0, int& i1, int& maxCodeLength) const
387 : {
388 4687 : 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 4687 : int size = (int)m_codeTable.size();
393 : {
394 4687 : int i = 0;
395 8296 : while (i < size && m_codeTable[i].first == 0) i++;
396 4687 : i0 = i;
397 4687 : i = size - 1;
398 30532 : while (i >= 0 && m_codeTable[i].first == 0) i--;
399 4681 : i1 = i + 1; // exclusive
400 : }
401 :
402 4681 : if (i1 <= i0)
403 0 : return false;
404 :
405 : // second, cover the common case that the peak is close to 0
406 4681 : pair<int, int> segm(0, 0);
407 4681 : int j = 0;
408 85352 : while (j < size) // find the largest stretch of 0's, if any
409 : {
410 264399 : while (j < size && m_codeTable[j].first > 0) j++;
411 80655 : int k0 = j;
412 1096490 : while (j < size && m_codeTable[j].first == 0) j++;
413 80641 : int k1 = j;
414 :
415 80641 : if (k1 - k0 > segm.second)
416 17597 : segm = pair<int, int>(k0, k1 - k0);
417 : }
418 :
419 4687 : if (size - segm.second < i1 - i0)
420 : {
421 4312 : i0 = segm.first + segm.second;
422 4312 : i1 = segm.first + size; // do wrap around
423 : }
424 :
425 4687 : if (i1 <= i0)
426 0 : return false;
427 :
428 4687 : int maxLen = 0;
429 658107 : for (int i = i0; i < i1; i++)
430 : {
431 653420 : int k = GetIndexWrapAround(i, size);
432 653418 : int len = m_codeTable[k].first;
433 653397 : maxLen = max(maxLen, len);
434 : }
435 :
436 4687 : if (maxLen <= 0 || maxLen > 32)
437 0 : return false;
438 :
439 4687 : maxCodeLength = maxLen;
440 4687 : 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 : int size = (int)m_codeTable.size();
451 640 : int bitPos = 0;
452 :
453 96748 : for (int i = i0; i < i1; i++)
454 : {
455 96108 : int k = GetIndexWrapAround(i, size);
456 96104 : int len = m_codeTable[k].first;
457 96109 : if (len > 0)
458 : {
459 27743 : unsigned int val = m_codeTable[k].second;
460 :
461 27743 : if (!Huffman::PushValue(ppByte, bitPos, val, len))
462 0 : return false;
463 : }
464 : }
465 :
466 640 : size_t numUInts = (bitPos > 0 ? 1 : 0);
467 640 : *ppByte += numUInts * sizeof(unsigned int);
468 :
469 640 : return true;
470 : }
471 :
472 : // -------------------------------------------------------------------------- ;
473 :
474 516 : bool Huffman::BitUnStuffCodes(const Byte** ppByte, size_t& nBytesRemainingInOut, int i0, int i1)
475 : {
476 516 : if (!ppByte || !(*ppByte))
477 0 : return false;
478 :
479 516 : size_t nBytesRemaining = nBytesRemainingInOut;
480 :
481 516 : const Byte* ptr0 = *ppByte;
482 516 : const Byte* ptr = ptr0;
483 :
484 516 : const size_t s4 = sizeof(unsigned int);
485 :
486 516 : int size = (int)m_codeTable.size();
487 516 : int bitPos = 0;
488 :
489 96719 : for (int i = i0; i < i1; i++)
490 : {
491 96203 : int k = GetIndexWrapAround(i, size);
492 96160 : int len = m_codeTable[k].first;
493 96231 : if (len > 0)
494 : {
495 31770 : if (nBytesRemaining < s4 || len > 32)
496 0 : return false;
497 :
498 31770 : unsigned int temp(0);
499 31770 : memcpy(&temp, ptr, s4);
500 31770 : m_codeTable[k].second = (temp << bitPos) >> (32 - len);
501 :
502 31744 : if (32 - bitPos >= len)
503 : {
504 24975 : bitPos += len;
505 24975 : if (bitPos == 32)
506 : {
507 872 : bitPos = 0;
508 872 : ptr += s4;
509 872 : nBytesRemaining -= s4;
510 : }
511 : }
512 : else
513 : {
514 6769 : bitPos += len - 32;
515 6769 : ptr += s4;
516 6769 : nBytesRemaining -= s4;
517 :
518 6769 : if (nBytesRemaining < s4)
519 0 : return false;
520 :
521 6769 : memcpy(&temp, ptr, s4);
522 6769 : m_codeTable[k].second |= temp >> (32 - bitPos); // bitPos > 0
523 : }
524 : }
525 : }
526 :
527 516 : size_t len = (ptr - ptr0) + (bitPos > 0 ? s4 : 0);
528 :
529 516 : if (nBytesRemainingInOut < len)
530 0 : return false;
531 :
532 516 : *ppByte += len;
533 516 : nBytesRemainingInOut -= len;
534 :
535 516 : if (nBytesRemaining != nBytesRemainingInOut
536 503 : && nBytesRemaining != nBytesRemainingInOut + s4) // the real check
537 0 : return false;
538 :
539 516 : return true;
540 : }
541 :
542 : // -------------------------------------------------------------------------- ;
543 :
544 : //struct MyLargerThanOp
545 : //{
546 : // inline bool operator() (const pair<int, unsigned int>& p0,
547 : // const pair<int, unsigned int>& p1) { return p0.first > p1.first; }
548 : //};
549 :
550 : // -------------------------------------------------------------------------- ;
551 :
552 3530 : bool Huffman::ConvertCodesToCanonical()
553 : {
554 : // from the non canonical code book, create an array to be sorted in descending order:
555 : // codeLength * tableSize - index
556 :
557 3530 : unsigned int tableSize = (unsigned int)m_codeTable.size();
558 3530 : if (tableSize == 0)
559 0 : return true;
560 3530 : vector<pair<int, unsigned int> > sortVec(tableSize, pair<int, unsigned int>(0, 0));
561 : //memset(&sortVec[0], 0, tableSize * sizeof(pair<int, unsigned int>));
562 :
563 907250 : for (unsigned int i = 0; i < tableSize; i++)
564 903719 : if (m_codeTable[i].first > 0)
565 124178 : sortVec[i] = pair<int, unsigned int>(m_codeTable[i].first * tableSize - i, i);
566 :
567 : // sort descending
568 : //std::sort(sortVec.begin(), sortVec.end(), MyLargerThanOp());
569 :
570 3531 : std::sort(sortVec.begin(), sortVec.end(),
571 5733400 : [](const pair<int, unsigned int>& p0,
572 5733400 : const pair<int, unsigned int>& p1) { return p0.first > p1.first; });
573 :
574 : // create canonical codes and assign to orig code table
575 3531 : unsigned int index = sortVec[0].second;
576 3531 : unsigned short codeLen = m_codeTable[index].first; // max code length for this table
577 3531 : unsigned int i = 0, codeCanonical = 0;
578 :
579 127728 : while (i < tableSize && sortVec[i].first > 0)
580 : {
581 124197 : index = sortVec[i++].second;
582 124197 : short delta = codeLen - m_codeTable[index].first; // difference of 2 consecutive code lengths, >= 0 as sorted
583 124197 : codeCanonical >>= delta;
584 124197 : codeLen -= delta;
585 124197 : m_codeTable[index].second = codeCanonical++;
586 : }
587 :
588 3531 : return true;
589 : }
590 :
591 : // -------------------------------------------------------------------------- ;
|
