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 907465 : for (int i = 0; i < size; i++) // add all leaf nodes
46 903934 : if (histo[i] > 0)
47 124195 : 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 120663 : Node* child1 = new Node(pq.top());
58 120666 : numNodes++;
59 120666 : pq.pop();
60 120665 : 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 517 : bool Huffman::ReadCodeTable(const Byte** ppByte, size_t& nBytesRemainingInOut, int lerc2Version)
171 : {
172 517 : if (!ppByte || !(*ppByte))
173 0 : return false;
174 :
175 517 : const Byte* ptr = *ppByte;
176 517 : size_t nBytesRemaining = nBytesRemainingInOut;
177 :
178 1034 : vector<int> intVec(4, 0);
179 517 : size_t len = intVec.size() * sizeof(int);
180 :
181 517 : if (nBytesRemaining < len)
182 0 : return false;
183 :
184 517 : memcpy(&intVec[0], ptr, len);
185 517 : ptr += len;
186 517 : nBytesRemaining -= len;
187 :
188 517 : int version = intVec[0];
189 :
190 517 : if (version < 2) // allow forward compatibility; for updates that break old decoders increase Lerc2 version number;
191 0 : return false;
192 :
193 517 : const int size = intVec[1];
194 517 : const int i0 = intVec[2];
195 517 : const int i1 = intVec[3];
196 :
197 517 : if (i0 >= i1 || i0 < 0 || size < 0 || size > (int)m_maxHistoSize)
198 0 : return false;
199 :
200 517 : if (GetIndexWrapAround(i0, size) >= size || GetIndexWrapAround(i1 - 1, size) >= size)
201 0 : return false;
202 :
203 : try
204 : {
205 1034 : vector<unsigned int> dataVec(i1 - i0, 0);
206 1034 : BitStuffer2 bitStuffer2;
207 517 : if (!bitStuffer2.Decode(&ptr, nBytesRemaining, dataVec, dataVec.size(), lerc2Version)) // unstuff the code lengths
208 0 : return false;
209 :
210 517 : if (dataVec.size() != static_cast<size_t>(i1 - i0))
211 0 : return false;
212 :
213 517 : m_codeTable.resize(size);
214 517 : std::fill(m_codeTable.begin(), m_codeTable.end(),
215 517 : std::pair<unsigned short, unsigned int>((short)0, 0));
216 :
217 97019 : for (int i = i0; i < i1; i++)
218 : {
219 96502 : int k = GetIndexWrapAround(i, size);
220 96502 : m_codeTable[k].first = (unsigned short)dataVec[i - i0];
221 : }
222 :
223 517 : if (!BitUnStuffCodes(&ptr, nBytesRemaining, i0, i1)) // unstuff the codes
224 0 : return false;
225 :
226 517 : *ppByte = ptr;
227 517 : nBytesRemainingInOut = nBytesRemaining;
228 517 : return true;
229 : }
230 0 : catch (std::exception&)
231 : {
232 0 : return false;
233 : }
234 : }
235 :
236 : // -------------------------------------------------------------------------- ;
237 :
238 517 : bool Huffman::BuildTreeFromCodes(int& numBitsLUT)
239 : {
240 517 : int i0 = 0, i1 = 0, maxLen = 0;
241 517 : if (!GetRange(i0, i1, maxLen))
242 0 : return false;
243 :
244 : // build decode LUT using max of 12 bits
245 517 : int size = (int)m_codeTable.size();
246 517 : int minNumZeroBits = 32;
247 :
248 517 : bool bNeedTree = maxLen > m_maxNumBitsLUT;
249 517 : numBitsLUT = min(maxLen, m_maxNumBitsLUT);
250 :
251 517 : int sizeLUT = 1 << numBitsLUT;
252 :
253 517 : m_decodeLUT.clear();
254 517 : m_decodeLUT.assign((size_t)sizeLUT, pair<short, short>((short)-1, (short)-1));
255 :
256 96978 : for (int i = i0; i < i1; i++)
257 : {
258 96461 : int k = GetIndexWrapAround(i, size);
259 96448 : int len = m_codeTable[k].first;
260 :
261 96472 : if (len == 0)
262 64513 : continue;
263 :
264 31959 : unsigned int code = m_codeTable[k].second;
265 :
266 31962 : if (len <= numBitsLUT)
267 : {
268 30774 : code <<= (numBitsLUT - len);
269 30774 : unsigned int numEntries = 1 << (numBitsLUT - len);
270 :
271 689062 : for (unsigned int j = 0; j < numEntries; j++)
272 : {
273 658376 : auto& entry = m_decodeLUT[code | j];
274 658288 : entry.first = (short)len; // add the duplicates
275 658288 : 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 1188 : int shift = 1;
281 4952 : while (code >>= 1) shift++; // large canonical codes start with zero's
282 1188 : minNumZeroBits = min(minNumZeroBits, len - shift);
283 : }
284 : }
285 :
286 517 : m_numBitsToSkipInTree = bNeedTree? minNumZeroBits : 0;
287 :
288 517 : if (!bNeedTree) // decode LUT covers it all, no tree needed
289 475 : 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 7811 : void Huffman::Clear()
338 : {
339 7811 : m_codeTable.clear();
340 7811 : m_decodeLUT.clear();
341 7811 : ClearTree();
342 7811 : }
343 :
344 : // -------------------------------------------------------------------------- ;
345 :
346 7853 : void Huffman::ClearTree()
347 : {
348 7853 : 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 7853 : }
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 4688 : bool Huffman::GetRange(int& i0, int& i1, int& maxCodeLength) const
387 : {
388 4688 : 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 4688 : int size = (int)m_codeTable.size();
393 : {
394 4688 : int i = 0;
395 8297 : while (i < size && m_codeTable[i].first == 0) i++;
396 4688 : i0 = i;
397 4688 : i = size - 1;
398 30639 : while (i >= 0 && m_codeTable[i].first == 0) i--;
399 4688 : i1 = i + 1; // exclusive
400 : }
401 :
402 4688 : if (i1 <= i0)
403 0 : return false;
404 :
405 : // second, cover the common case that the peak is close to 0
406 4688 : pair<int, int> segm(0, 0);
407 4688 : int j = 0;
408 85385 : while (j < size) // find the largest stretch of 0's, if any
409 : {
410 264612 : while (j < size && m_codeTable[j].first > 0) j++;
411 80697 : int k0 = j;
412 1096910 : while (j < size && m_codeTable[j].first == 0) j++;
413 80697 : int k1 = j;
414 :
415 80697 : if (k1 - k0 > segm.second)
416 17604 : segm = pair<int, int>(k0, k1 - k0);
417 : }
418 :
419 4688 : if (size - segm.second < i1 - i0)
420 : {
421 4313 : i0 = segm.first + segm.second;
422 4313 : i1 = segm.first + size; // do wrap around
423 : }
424 :
425 4688 : if (i1 <= i0)
426 0 : return false;
427 :
428 4688 : int maxLen = 0;
429 658560 : for (int i = i0; i < i1; i++)
430 : {
431 653872 : int k = GetIndexWrapAround(i, size);
432 653872 : int len = m_codeTable[k].first;
433 653872 : maxLen = max(maxLen, len);
434 : }
435 :
436 4688 : if (maxLen <= 0 || maxLen > 32)
437 0 : return false;
438 :
439 4688 : maxCodeLength = maxLen;
440 4688 : 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 96762 : for (int i = i0; i < i1; i++)
454 : {
455 96122 : int k = GetIndexWrapAround(i, size);
456 96122 : int len = m_codeTable[k].first;
457 96122 : 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 517 : bool Huffman::BitUnStuffCodes(const Byte** ppByte, size_t& nBytesRemainingInOut, int i0, int i1)
475 : {
476 517 : if (!ppByte || !(*ppByte))
477 0 : return false;
478 :
479 517 : size_t nBytesRemaining = nBytesRemainingInOut;
480 :
481 517 : const Byte* ptr0 = *ppByte;
482 517 : const Byte* ptr = ptr0;
483 :
484 517 : const size_t s4 = sizeof(unsigned int);
485 :
486 517 : int size = (int)m_codeTable.size();
487 517 : int bitPos = 0;
488 :
489 97019 : for (int i = i0; i < i1; i++)
490 : {
491 96502 : int k = GetIndexWrapAround(i, size);
492 96502 : int len = m_codeTable[k].first;
493 96502 : if (len > 0)
494 : {
495 31975 : if (nBytesRemaining < s4 || len > 32)
496 0 : return false;
497 :
498 31975 : unsigned int temp(0);
499 31975 : memcpy(&temp, ptr, s4);
500 31975 : m_codeTable[k].second = (temp << bitPos) >> (32 - len);
501 :
502 31975 : if (32 - bitPos >= len)
503 : {
504 25177 : bitPos += len;
505 25177 : if (bitPos == 32)
506 : {
507 876 : bitPos = 0;
508 876 : ptr += s4;
509 876 : nBytesRemaining -= s4;
510 : }
511 : }
512 : else
513 : {
514 6798 : bitPos += len - 32;
515 6798 : ptr += s4;
516 6798 : nBytesRemaining -= s4;
517 :
518 6798 : if (nBytesRemaining < s4)
519 0 : return false;
520 :
521 6798 : memcpy(&temp, ptr, s4);
522 6798 : m_codeTable[k].second |= temp >> (32 - bitPos); // bitPos > 0
523 : }
524 : }
525 : }
526 :
527 517 : size_t len = (ptr - ptr0) + (bitPos > 0 ? s4 : 0);
528 :
529 517 : if (nBytesRemainingInOut < len)
530 0 : return false;
531 :
532 517 : *ppByte += len;
533 517 : nBytesRemainingInOut -= len;
534 :
535 517 : if (nBytesRemaining != nBytesRemainingInOut
536 504 : && nBytesRemaining != nBytesRemainingInOut + s4) // the real check
537 0 : return false;
538 :
539 517 : 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 3531 : 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 3531 : unsigned int tableSize = (unsigned int)m_codeTable.size();
558 3531 : if (tableSize == 0)
559 0 : return true;
560 3531 : 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 907467 : for (unsigned int i = 0; i < tableSize; i++)
564 903936 : if (m_codeTable[i].first > 0)
565 124197 : 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 5734520 : [](const pair<int, unsigned int>& p0,
572 5734520 : 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 : // -------------------------------------------------------------------------- ;
|
