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 Lucian Plesea (provided checksum code)
23 */
24
25 #include <typeinfo>
26 #include "Defines.h"
27 #include "Lerc2.h"
28 #include "Huffman.h"
29 #include "RLE.h"
30
31 USING_NAMESPACE_LERC
32 using namespace std;
33
34 // -------------------------------------------------------------------------- ;
35
Lerc2()36 Lerc2::Lerc2()
37 {
38 Init();
39 }
40
41 // -------------------------------------------------------------------------- ;
42
Lerc2(int nDim,int nCols,int nRows,const Byte * pMaskBits)43 Lerc2::Lerc2(int nDim, int nCols, int nRows, const Byte* pMaskBits)
44 {
45 Init();
46 Set(nDim, nCols, nRows, pMaskBits);
47 }
48
49 // -------------------------------------------------------------------------- ;
50
SetEncoderToOldVersion(int version)51 bool Lerc2::SetEncoderToOldVersion(int version)
52 {
53 if (version < 2 || version > CurrentVersion())
54 return false;
55
56 if (version < 4 && m_headerInfo.nDim > 1)
57 return false;
58
59 m_headerInfo.version = version;
60
61 return true;
62 }
63
64 // -------------------------------------------------------------------------- ;
65
Init()66 void Lerc2::Init()
67 {
68 m_microBlockSize = 8;
69 m_maxValToQuantize = 0;
70 m_encodeMask = true;
71 m_writeDataOneSweep = false;
72 m_imageEncodeMode = IEM_Tiling;
73
74 m_headerInfo.RawInit();
75 m_headerInfo.version = CurrentVersion();
76 m_headerInfo.microBlockSize = m_microBlockSize;
77 }
78
79 // -------------------------------------------------------------------------- ;
80
Set(int nDim,int nCols,int nRows,const Byte * pMaskBits)81 bool Lerc2::Set(int nDim, int nCols, int nRows, const Byte* pMaskBits)
82 {
83 if (nDim > 1 && m_headerInfo.version < 4)
84 return false;
85
86 if (!m_bitMask.SetSize(nCols, nRows))
87 return false;
88
89 if (pMaskBits)
90 {
91 memcpy(m_bitMask.Bits(), pMaskBits, m_bitMask.Size());
92 m_headerInfo.numValidPixel = m_bitMask.CountValidBits();
93 }
94 else
95 {
96 m_headerInfo.numValidPixel = nCols * nRows;
97 m_bitMask.SetAllValid();
98 }
99
100 m_headerInfo.nDim = nDim;
101 m_headerInfo.nCols = nCols;
102 m_headerInfo.nRows = nRows;
103
104 return true;
105 }
106
107 // -------------------------------------------------------------------------- ;
108
109 template<class T>
ComputeNumBytesNeededToWrite(const T * arr,double maxZError,bool encodeMask)110 unsigned int Lerc2::ComputeNumBytesNeededToWrite(const T* arr, double maxZError, bool encodeMask)
111 {
112 if (!arr || !IsLittleEndianSystem())
113 return 0;
114
115 // header
116 unsigned int nBytesHeaderMask = ComputeNumBytesHeaderToWrite(m_headerInfo);
117
118 // valid / invalid mask
119 int numValid = m_headerInfo.numValidPixel;
120 int numTotal = m_headerInfo.nCols * m_headerInfo.nRows;
121
122 bool needMask = numValid > 0 && numValid < numTotal;
123
124 m_encodeMask = encodeMask;
125
126 nBytesHeaderMask += 1 * sizeof(int); // the mask encode numBytes
127
128 if (needMask && encodeMask)
129 {
130 RLE rle;
131 size_t n = rle.computeNumBytesRLE((const Byte*)m_bitMask.Bits(), m_bitMask.Size());
132 nBytesHeaderMask += (unsigned int)n;
133 }
134
135 m_headerInfo.dt = GetDataType(arr[0]);
136
137 if (m_headerInfo.dt == DT_Undefined)
138 return 0;
139
140 if (maxZError == 777) // cheat code
141 maxZError = -0.01;
142
143 if (m_headerInfo.dt < DT_Float) // integer types
144 {
145 // interpret a negative maxZError as bit plane epsilon; dflt = 0.01;
146 if (maxZError < 0 && (!TryBitPlaneCompression(arr, -maxZError, maxZError)))
147 maxZError = 0;
148
149 maxZError = std::max(0.5, floor(maxZError));
150 }
151 else // float types
152 {
153 if (maxZError < 0) // don't allow bit plane compression for float or double
154 return 0;
155
156 double maxZErrorNew = maxZError;
157 if (TryRaiseMaxZError(arr, maxZErrorNew))
158 {
159 //printf("%f --> %f\n", maxZError, maxZErrorNew);
160 maxZError = maxZErrorNew;
161 }
162 }
163
164 m_headerInfo.maxZError = maxZError;
165 m_headerInfo.zMin = 0;
166 m_headerInfo.zMax = 0;
167 m_headerInfo.microBlockSize = m_microBlockSize;
168 m_headerInfo.blobSize = nBytesHeaderMask;
169
170 if (numValid == 0)
171 return nBytesHeaderMask;
172
173 m_maxValToQuantize = GetMaxValToQuantize(m_headerInfo.dt);
174
175 Byte* ptr = nullptr; // only emulate the writing and just count the bytes needed
176 int nBytesTiling = 0;
177
178 if (!ComputeMinMaxRanges(arr, m_zMinVec, m_zMaxVec)) // need this for diff encoding before WriteTiles()
179 return 0;
180
181 m_headerInfo.zMin = *std::min_element(m_zMinVec.begin(), m_zMinVec.end());
182 m_headerInfo.zMax = *std::max_element(m_zMaxVec.begin(), m_zMaxVec.end());
183
184 if (m_headerInfo.zMin == m_headerInfo.zMax) // image is const
185 return nBytesHeaderMask;
186
187 int nDim = m_headerInfo.nDim;
188
189 if (m_headerInfo.version >= 4)
190 {
191 // add the min max ranges behind the mask and before the main data;
192 // so we do not write it if no valid pixel or all same value const
193 m_headerInfo.blobSize += 2 * nDim * sizeof(T);
194
195 bool minMaxEqual = false;
196 if (!CheckMinMaxRanges(minMaxEqual))
197 return 0;
198
199 if (minMaxEqual)
200 return m_headerInfo.blobSize; // all nDim bands are const
201 }
202
203 // data
204 if (!WriteTiles(arr, &ptr, nBytesTiling))
205 return 0;
206
207 m_imageEncodeMode = IEM_Tiling;
208 int nBytesData = nBytesTiling;
209 int nBytesHuffman = 0;
210
211 if (m_headerInfo.TryHuffman())
212 {
213 ImageEncodeMode huffmanEncMode;
214 ComputeHuffmanCodes(arr, nBytesHuffman, huffmanEncMode, m_huffmanCodes); // save Huffman codes for later use
215
216 if (!m_huffmanCodes.empty() && nBytesHuffman < nBytesTiling)
217 {
218 m_imageEncodeMode = huffmanEncMode;
219 nBytesData = nBytesHuffman;
220 }
221 else
222 m_huffmanCodes.resize(0);
223 }
224
225 m_writeDataOneSweep = false;
226 int nBytesDataOneSweep = (int)(numValid * nDim * sizeof(T));
227
228 {
229 // try with double block size to reduce block header overhead, if
230 if ((nBytesTiling * 8 < numTotal * nDim * 1.5) // resulting bit rate < x (2 bpp)
231 && (nBytesTiling < 4 * nBytesDataOneSweep) // bit stuffing is effective
232 && (nBytesHuffman == 0 || nBytesTiling < 2 * nBytesHuffman) // not much worse than huffman (otherwise huffman wins anyway)
233 && (m_headerInfo.nRows > m_microBlockSize || m_headerInfo.nCols > m_microBlockSize))
234 {
235 m_headerInfo.microBlockSize = m_microBlockSize * 2;
236
237 int nBytes2 = 0;
238 if (!WriteTiles(arr, &ptr, nBytes2)) // no huffman in here anymore
239 return 0;
240
241 if (nBytes2 <= nBytesData)
242 {
243 nBytesData = nBytes2;
244 m_imageEncodeMode = IEM_Tiling;
245 m_huffmanCodes.resize(0);
246 }
247 else
248 {
249 m_headerInfo.microBlockSize = m_microBlockSize; // reset to orig
250 }
251 }
252 }
253
254 if (m_headerInfo.TryHuffman())
255 nBytesData += 1; // flag for image encode mode
256
257 if (nBytesDataOneSweep <= nBytesData)
258 {
259 m_writeDataOneSweep = true; // fallback: write data binary uncompressed in one sweep
260 m_headerInfo.blobSize += 1 + nBytesDataOneSweep; // header, mask, min max ranges, flag, data one sweep
261 }
262 else
263 {
264 m_writeDataOneSweep = false;
265 m_headerInfo.blobSize += 1 + nBytesData; // header, mask, min max ranges, flag(s), data
266 }
267
268 return m_headerInfo.blobSize;
269 }
270
271 // -------------------------------------------------------------------------- ;
272
273 template unsigned int Lerc2::ComputeNumBytesNeededToWrite<signed char>(const signed char* arr, double maxZError, bool encodeMask);
274 template unsigned int Lerc2::ComputeNumBytesNeededToWrite<Byte>(const Byte* arr, double maxZError, bool encodeMask);
275 template unsigned int Lerc2::ComputeNumBytesNeededToWrite<short>(const short* arr, double maxZError, bool encodeMask);
276 template unsigned int Lerc2::ComputeNumBytesNeededToWrite<unsigned short>(const unsigned short* arr, double maxZError, bool encodeMask);
277 template unsigned int Lerc2::ComputeNumBytesNeededToWrite<int>(const int* arr, double maxZError, bool encodeMask);
278 template unsigned int Lerc2::ComputeNumBytesNeededToWrite<unsigned int>(const unsigned int* arr, double maxZError, bool encodeMask);
279 template unsigned int Lerc2::ComputeNumBytesNeededToWrite<float>(const float* arr, double maxZError, bool encodeMask);
280 template unsigned int Lerc2::ComputeNumBytesNeededToWrite<double>(const double* arr, double maxZError, bool encodeMask);
281
282 // -------------------------------------------------------------------------- ;
283
284 template<class T>
Encode(const T * arr,Byte ** ppByte)285 bool Lerc2::Encode(const T* arr, Byte** ppByte)
286 {
287 if (!arr || !ppByte || !IsLittleEndianSystem())
288 return false;
289
290 Byte* ptrBlob = *ppByte; // keep a ptr to the start of the blob
291
292 if (!WriteHeader(ppByte, m_headerInfo))
293 return false;
294
295 if (!WriteMask(ppByte))
296 return false;
297
298 if (m_headerInfo.numValidPixel == 0 || m_headerInfo.zMin == m_headerInfo.zMax)
299 {
300 return DoChecksOnEncode(ptrBlob, *ppByte);
301 }
302
303 if (m_headerInfo.version >= 4)
304 {
305 if (!WriteMinMaxRanges(arr, ppByte))
306 return false;
307
308 bool minMaxEqual = false;
309 if (!CheckMinMaxRanges(minMaxEqual))
310 return false;
311
312 if (minMaxEqual)
313 return DoChecksOnEncode(ptrBlob, *ppByte);
314 }
315
316 **ppByte = m_writeDataOneSweep ? 1 : 0; // write flag
317 (*ppByte)++;
318
319 if (!m_writeDataOneSweep)
320 {
321 if (m_headerInfo.TryHuffman())
322 {
323 **ppByte = (Byte)m_imageEncodeMode; // Huffman or tiling encode mode
324 (*ppByte)++;
325
326 if (!m_huffmanCodes.empty()) // Huffman, no tiling
327 {
328 if (m_imageEncodeMode != IEM_DeltaHuffman && m_imageEncodeMode != IEM_Huffman)
329 return false;
330
331 if (!EncodeHuffman(arr, ppByte)) // data bit stuffed
332 return false;
333
334 return DoChecksOnEncode(ptrBlob, *ppByte);
335 }
336 }
337
338 int numBytes = 0;
339 if (!WriteTiles(arr, ppByte, numBytes))
340 return false;
341 }
342 else
343 {
344 if (!WriteDataOneSweep(arr, ppByte))
345 return false;
346 }
347
348 return DoChecksOnEncode(ptrBlob, *ppByte);
349 }
350
351 // -------------------------------------------------------------------------- ;
352
353 template bool Lerc2::Encode<signed char>(const signed char* arr, Byte** ppByte);
354 template bool Lerc2::Encode<Byte>(const Byte* arr, Byte** ppByte);
355 template bool Lerc2::Encode<short>(const short* arr, Byte** ppByte);
356 template bool Lerc2::Encode<unsigned short>(const unsigned short* arr, Byte** ppByte);
357 template bool Lerc2::Encode<int>(const int* arr, Byte** ppByte);
358 template bool Lerc2::Encode<unsigned int>(const unsigned int* arr, Byte** ppByte);
359 template bool Lerc2::Encode<float>(const float* arr, Byte** ppByte);
360 template bool Lerc2::Encode<double>(const double* arr, Byte** ppByte);
361
362 // -------------------------------------------------------------------------- ;
363
GetHeaderInfo(const Byte * pByte,size_t nBytesRemaining,struct HeaderInfo & hd,bool & bHasMask)364 bool Lerc2::GetHeaderInfo(const Byte* pByte, size_t nBytesRemaining, struct HeaderInfo& hd, bool& bHasMask)
365 {
366 if (!pByte || !IsLittleEndianSystem())
367 return false;
368
369 if (!ReadHeader(&pByte, nBytesRemaining, hd))
370 return false;
371
372 int numBytesMask(0);
373 if (nBytesRemaining < sizeof(int) || !memcpy(&numBytesMask, pByte, sizeof(int)))
374 return false;
375
376 bHasMask = numBytesMask > 0;
377
378 return true;
379 }
380
381 // -------------------------------------------------------------------------- ;
382
383 template<class T>
Decode(const Byte ** ppByte,size_t & nBytesRemaining,T * arr,Byte * pMaskBits)384 bool Lerc2::Decode(const Byte** ppByte, size_t& nBytesRemaining, T* arr, Byte* pMaskBits)
385 {
386 if (!arr || !ppByte || !IsLittleEndianSystem())
387 return false;
388
389 const Byte* ptrBlob = *ppByte; // keep a ptr to the start of the blob
390 size_t nBytesRemaining00 = nBytesRemaining;
391
392 if (!ReadHeader(ppByte, nBytesRemaining, m_headerInfo))
393 return false;
394
395 if (nBytesRemaining00 < (size_t)m_headerInfo.blobSize)
396 return false;
397
398 if (m_headerInfo.version >= 3)
399 {
400 int nBytes = (int)(FileKey().length() + sizeof(int) + sizeof(unsigned int)); // start right after the checksum entry
401 if (m_headerInfo.blobSize < nBytes)
402 return false;
403 unsigned int checksum = ComputeChecksumFletcher32(ptrBlob + nBytes, m_headerInfo.blobSize - nBytes);
404
405 if (checksum != m_headerInfo.checksum)
406 return false;
407 }
408
409 if (!ReadMask(ppByte, nBytesRemaining))
410 return false;
411
412 if (pMaskBits) // return proper mask bits even if they were not stored
413 memcpy(pMaskBits, m_bitMask.Bits(), m_bitMask.Size());
414
415 memset(arr, 0, m_headerInfo.nCols * m_headerInfo.nRows * m_headerInfo.nDim * sizeof(T));
416
417 if (m_headerInfo.numValidPixel == 0)
418 return true;
419
420 if (m_headerInfo.zMin == m_headerInfo.zMax) // image is const
421 {
422 if (!FillConstImage(arr))
423 return false;
424
425 return true;
426 }
427
428 if (m_headerInfo.version >= 4)
429 {
430 if (!ReadMinMaxRanges(ppByte, nBytesRemaining, arr))
431 return false;
432
433 bool minMaxEqual = false;
434 if (!CheckMinMaxRanges(minMaxEqual))
435 return false;
436
437 if (minMaxEqual) // if all bands are const, fill outgoing and done
438 {
439 if (!FillConstImage(arr))
440 return false;
441
442 return true; // done
443 }
444 }
445
446 if (nBytesRemaining < 1)
447 return false;
448
449 Byte readDataOneSweep = **ppByte; // read flag
450 (*ppByte)++;
451 nBytesRemaining--;
452
453 if (!readDataOneSweep)
454 {
455 if (m_headerInfo.TryHuffman())
456 {
457 if (nBytesRemaining < 1)
458 return false;
459
460 Byte flag = **ppByte; // read flag Huffman / Lerc2
461 (*ppByte)++;
462 nBytesRemaining--;
463
464 if (flag > 2 || (m_headerInfo.version < 4 && flag > 1))
465 return false;
466
467 m_imageEncodeMode = (ImageEncodeMode)flag;
468
469 if (m_imageEncodeMode == IEM_DeltaHuffman || m_imageEncodeMode == IEM_Huffman)
470 {
471 if (!DecodeHuffman(ppByte, nBytesRemaining, arr))
472 return false;
473
474 return true; // done.
475 }
476 }
477
478 if (!ReadTiles(ppByte, nBytesRemaining, arr))
479 return false;
480 }
481 else
482 {
483 if (!ReadDataOneSweep(ppByte, nBytesRemaining, arr))
484 return false;
485 }
486
487 return true;
488 }
489
490 // -------------------------------------------------------------------------- ;
491
492 template bool Lerc2::Decode<signed char>(const Byte** ppByte, size_t& nBytesRemaining, signed char* arr, Byte* pMaskBits);
493 template bool Lerc2::Decode<Byte>(const Byte** ppByte, size_t& nBytesRemaining, Byte* arr, Byte* pMaskBits);
494 template bool Lerc2::Decode<short>(const Byte** ppByte, size_t& nBytesRemaining, short* arr, Byte* pMaskBits);
495 template bool Lerc2::Decode<unsigned short>(const Byte** ppByte, size_t& nBytesRemaining, unsigned short* arr, Byte* pMaskBits);
496 template bool Lerc2::Decode<int>(const Byte** ppByte, size_t& nBytesRemaining, int* arr, Byte* pMaskBits);
497 template bool Lerc2::Decode<unsigned int>(const Byte** ppByte, size_t& nBytesRemaining, unsigned int* arr, Byte* pMaskBits);
498 template bool Lerc2::Decode<float>(const Byte** ppByte, size_t& nBytesRemaining, float* arr, Byte* pMaskBits);
499 template bool Lerc2::Decode<double>(const Byte** ppByte, size_t& nBytesRemaining, double* arr, Byte* pMaskBits);
500
501 // -------------------------------------------------------------------------- ;
502 // -------------------------------------------------------------------------- ;
503
ComputeNumBytesHeaderToWrite(const struct HeaderInfo & hd)504 unsigned int Lerc2::ComputeNumBytesHeaderToWrite(const struct HeaderInfo& hd)
505 {
506 unsigned int numBytes = (unsigned int)FileKey().length();
507 numBytes += 1 * sizeof(int);
508 numBytes += (hd.version >= 3 ? 1 : 0) * sizeof(unsigned int);
509 numBytes += (hd.version >= 4 ? 7 : 6) * sizeof(int);
510 numBytes += 3 * sizeof(double);
511 return numBytes;
512 }
513
514 // -------------------------------------------------------------------------- ;
515
WriteHeader(Byte ** ppByte,const struct HeaderInfo & hd)516 bool Lerc2::WriteHeader(Byte** ppByte, const struct HeaderInfo& hd)
517 {
518 if (!ppByte)
519 return false;
520
521 Byte* ptr = *ppByte;
522
523 string fileKey = FileKey();
524 size_t len = fileKey.length();
525 memcpy(ptr, fileKey.c_str(), len);
526 ptr += len;
527
528 memcpy(ptr, &hd.version, sizeof(int));
529 ptr += sizeof(int);
530
531 if (hd.version >= 3)
532 {
533 unsigned int checksum = 0;
534 memcpy(ptr, &checksum, sizeof(unsigned int)); // place holder to be filled by the real check sum later
535 ptr += sizeof(unsigned int);
536 }
537
538 vector<int> intVec;
539 intVec.push_back(hd.nRows);
540 intVec.push_back(hd.nCols);
541
542 if (hd.version >= 4)
543 {
544 intVec.push_back(hd.nDim);
545 }
546
547 intVec.push_back(hd.numValidPixel);
548 intVec.push_back(hd.microBlockSize);
549 intVec.push_back(hd.blobSize);
550 intVec.push_back((int)hd.dt);
551
552 len = intVec.size() * sizeof(int);
553 memcpy(ptr, &intVec[0], len);
554 ptr += len;
555
556 vector<double> dblVec;
557 dblVec.push_back(hd.maxZError);
558 dblVec.push_back(hd.zMin);
559 dblVec.push_back(hd.zMax);
560
561 len = dblVec.size() * sizeof(double);
562 memcpy(ptr, &dblVec[0], len);
563 ptr += len;
564
565 *ppByte = ptr;
566 return true;
567 }
568
569 // -------------------------------------------------------------------------- ;
570
ReadHeader(const Byte ** ppByte,size_t & nBytesRemainingInOut,struct HeaderInfo & hd)571 bool Lerc2::ReadHeader(const Byte** ppByte, size_t& nBytesRemainingInOut, struct HeaderInfo& hd)
572 {
573 if (!ppByte || !*ppByte)
574 return false;
575
576 const Byte* ptr = *ppByte;
577 size_t nBytesRemaining = nBytesRemainingInOut;
578
579 string fileKey = FileKey();
580 size_t keyLen = fileKey.length();
581
582 hd.RawInit();
583
584 if (nBytesRemaining < keyLen || memcmp(ptr, fileKey.c_str(), keyLen))
585 return false;
586
587 ptr += keyLen;
588 nBytesRemaining -= keyLen;
589
590 if (nBytesRemaining < sizeof(int) || !memcpy(&(hd.version), ptr, sizeof(int)))
591 return false;
592
593 ptr += sizeof(int);
594 nBytesRemaining -= sizeof(int);
595
596 if (hd.version < 0 || hd.version > CurrentVersion()) // this reader is outdated
597 return false;
598
599 if (hd.version >= 3)
600 {
601 if (nBytesRemaining < sizeof(unsigned int) || !memcpy(&(hd.checksum), ptr, sizeof(unsigned int)))
602 return false;
603
604 ptr += sizeof(unsigned int);
605 nBytesRemaining -= sizeof(unsigned int);
606 }
607
608 int nInts = (hd.version >= 4) ? 7 : 6;
609 vector<int> intVec(nInts, 0);
610 vector<double> dblVec(3, 0);
611
612 size_t len = sizeof(int) * intVec.size();
613
614 if (nBytesRemaining < len || !memcpy(&intVec[0], ptr, len))
615 return false;
616
617 ptr += len;
618 nBytesRemaining -= len;
619
620 len = sizeof(double) * dblVec.size();
621
622 if (nBytesRemaining < len || !memcpy(&dblVec[0], ptr, len))
623 return false;
624
625 ptr += len;
626 nBytesRemaining -= len;
627
628 int i = 0;
629 hd.nRows = intVec[i++];
630 hd.nCols = intVec[i++];
631 hd.nDim = (hd.version >= 4) ? intVec[i++] : 1;
632 hd.numValidPixel = intVec[i++];
633 hd.microBlockSize = intVec[i++];
634 hd.blobSize = intVec[i++];
635 const int dt = intVec[i++];
636 if (dt < DT_Char || dt > DT_Double)
637 return false;
638 hd.dt = static_cast<DataType>(dt);
639
640 hd.maxZError = dblVec[0];
641 hd.zMin = dblVec[1];
642 hd.zMax = dblVec[2];
643
644 if (hd.nRows <= 0 || hd.nCols <= 0 || hd.nDim <= 0 || hd.numValidPixel < 0 || hd.microBlockSize <= 0 || hd.blobSize <= 0
645 || hd.numValidPixel > hd.nRows * hd.nCols)
646 return false;
647
648 *ppByte = ptr;
649 nBytesRemainingInOut = nBytesRemaining;
650
651 return true;
652 }
653
654 // -------------------------------------------------------------------------- ;
655
WriteMask(Byte ** ppByte) const656 bool Lerc2::WriteMask(Byte** ppByte) const
657 {
658 if (!ppByte)
659 return false;
660
661 int numValid = m_headerInfo.numValidPixel;
662 int numTotal = m_headerInfo.nCols * m_headerInfo.nRows;
663
664 bool needMask = numValid > 0 && numValid < numTotal;
665
666 Byte* ptr = *ppByte;
667
668 if (needMask && m_encodeMask)
669 {
670 Byte* pArrRLE;
671 size_t numBytesRLE;
672 RLE rle;
673 if (!rle.compress((const Byte*)m_bitMask.Bits(), m_bitMask.Size(), &pArrRLE, numBytesRLE, false))
674 return false;
675
676 int numBytesMask = (int)numBytesRLE;
677 memcpy(ptr, &numBytesMask, sizeof(int)); // num bytes for compressed mask
678 ptr += sizeof(int);
679 memcpy(ptr, pArrRLE, numBytesRLE);
680 ptr += numBytesRLE;
681
682 delete[] pArrRLE;
683 }
684 else
685 {
686 memset(ptr, 0, sizeof(int)); // indicates no mask stored
687 ptr += sizeof(int);
688 }
689
690 *ppByte = ptr;
691 return true;
692 }
693
694 // -------------------------------------------------------------------------- ;
695
ReadMask(const Byte ** ppByte,size_t & nBytesRemainingInOut)696 bool Lerc2::ReadMask(const Byte** ppByte, size_t& nBytesRemainingInOut)
697 {
698 if (!ppByte)
699 return false;
700
701 int numValid = m_headerInfo.numValidPixel;
702 int w = m_headerInfo.nCols;
703 int h = m_headerInfo.nRows;
704
705 const Byte* ptr = *ppByte;
706 size_t nBytesRemaining = nBytesRemainingInOut;
707
708 int numBytesMask;
709 if (nBytesRemaining < sizeof(int) || !memcpy(&numBytesMask, ptr, sizeof(int)))
710 return false;
711
712 ptr += sizeof(int);
713 nBytesRemaining -= sizeof(int);
714
715 if (numValid == 0 || numValid == w * h) // there should be no mask
716 {
717 if (numBytesMask != 0)
718 return false;
719 }
720
721 if (!m_bitMask.SetSize(w, h))
722 return false;
723
724 if (numValid == 0)
725 m_bitMask.SetAllInvalid();
726 else if (numValid == w * h)
727 m_bitMask.SetAllValid();
728 else if (numBytesMask > 0) // read it in
729 {
730 if (nBytesRemaining < static_cast<size_t>(numBytesMask))
731 return false;
732
733 RLE rle;
734 if (!rle.decompress(ptr, nBytesRemaining, m_bitMask.Bits(), m_bitMask.Size()))
735 return false;
736
737 ptr += numBytesMask;
738 nBytesRemaining -= numBytesMask;
739 }
740 // else use previous mask
741
742 *ppByte = ptr;
743 nBytesRemainingInOut = nBytesRemaining;
744
745 return true;
746 }
747
748 // -------------------------------------------------------------------------- ;
749
DoChecksOnEncode(Byte * pBlobBegin,Byte * pBlobEnd) const750 bool Lerc2::DoChecksOnEncode(Byte* pBlobBegin, Byte* pBlobEnd) const
751 {
752 if ((size_t)(pBlobEnd - pBlobBegin) != (size_t)m_headerInfo.blobSize)
753 return false;
754
755 if (m_headerInfo.version >= 3)
756 {
757 int blobSize = (int)(pBlobEnd - pBlobBegin);
758 int nBytes = (int)(FileKey().length() + sizeof(int) + sizeof(unsigned int)); // start right after the checksum entry
759 if (blobSize < nBytes)
760 return false;
761 unsigned int checksum = ComputeChecksumFletcher32(pBlobBegin + nBytes, blobSize - nBytes);
762
763 nBytes -= sizeof(unsigned int);
764 memcpy(pBlobBegin + nBytes, &checksum, sizeof(unsigned int));
765 }
766
767 return true;
768 }
769
770 // -------------------------------------------------------------------------- ;
771
772 // from https://en.wikipedia.org/wiki/Fletcher's_checksum
773 // modified from ushorts to bytes (by Lucian Plesea)
774
ComputeChecksumFletcher32(const Byte * pByte,int len)775 unsigned int Lerc2::ComputeChecksumFletcher32(const Byte* pByte, int len)
776 {
777 unsigned int sum1 = 0xffff, sum2 = 0xffff;
778 unsigned int words = len / 2;
779
780 while (words)
781 {
782 unsigned int tlen = (words >= 359) ? 359 : words;
783 words -= tlen;
784 do {
785 sum1 += (*pByte++ << 8);
786 sum2 += sum1 += *pByte++;
787 } while (--tlen);
788
789 sum1 = (sum1 & 0xffff) + (sum1 >> 16);
790 sum2 = (sum2 & 0xffff) + (sum2 >> 16);
791 }
792
793 // add the straggler byte if it exists
794 if (len & 1)
795 sum2 += sum1 += (*pByte << 8);
796
797 // second reduction step to reduce sums to 16 bits
798 sum1 = (sum1 & 0xffff) + (sum1 >> 16);
799 sum2 = (sum2 & 0xffff) + (sum2 >> 16);
800
801 return sum2 << 16 | sum1;
802 }
803
804
805 // -------------------------------------------------------------------------- ;
806
807 // for the theory and math, see
808 // https://pdfs.semanticscholar.org/d064/2e2ad1a4c3b445b0d795770f604a5d9e269c.pdf
809
810 template<class T>
TryBitPlaneCompression(const T * data,double eps,double & newMaxZError) const811 bool Lerc2::TryBitPlaneCompression(const T* data, double eps, double& newMaxZError) const
812 {
813 newMaxZError = 0; // lossless is the obvious fallback
814
815 if (!data || eps <= 0)
816 return false;
817
818 const HeaderInfo& hd = m_headerInfo;
819 const int nDim = hd.nDim;
820 const int maxShift = 8 * GetDataTypeSize(hd.dt);
821 const int minCnt = 5000;
822
823 if (hd.numValidPixel < minCnt) // not enough data for good stats
824 return false;
825
826 std::vector<int> cntDiffVec(nDim * maxShift, 0);
827 int cnt = 0;
828
829 if (nDim == 1 && hd.numValidPixel == hd.nCols * hd.nRows) // special but common case
830 {
831 if (hd.dt == DT_Byte || hd.dt == DT_UShort || hd.dt == DT_UInt) // unsigned int
832 {
833 for (int i = 0; i < hd.nRows - 1; i++)
834 for (int k = i * hd.nCols, j = 0; j < hd.nCols - 1; j++, k++)
835 {
836 unsigned int c = ((unsigned int)data[k]) ^ ((unsigned int)data[k + 1]);
837 AddUIntToCounts(&cntDiffVec[0], c, maxShift);
838 cnt++;
839 c = ((unsigned int)data[k]) ^ ((unsigned int)data[k + hd.nCols]);
840 AddUIntToCounts(&cntDiffVec[0], c, maxShift);
841 cnt++;
842 }
843 }
844 else if (hd.dt == DT_Char || hd.dt == DT_Short || hd.dt == DT_Int) // signed int
845 {
846 for (int i = 0; i < hd.nRows - 1; i++)
847 for (int k = i * hd.nCols, j = 0; j < hd.nCols - 1; j++, k++)
848 {
849 int c = ((int)data[k]) ^ ((int)data[k + 1]);
850 AddIntToCounts(&cntDiffVec[0], c, maxShift);
851 cnt++;
852 c = ((int)data[k]) ^ ((int)data[k + hd.nCols]);
853 AddIntToCounts(&cntDiffVec[0], c, maxShift);
854 cnt++;
855 }
856 }
857 else
858 return false; // unsupported data type
859 }
860
861 else // general case: nDim > 1 or not all pixel valid
862 {
863 if (hd.dt == DT_Byte || hd.dt == DT_UShort || hd.dt == DT_UInt) // unsigned int
864 {
865 for (int k = 0, m0 = 0, i = 0; i < hd.nRows; i++)
866 for (int j = 0; j < hd.nCols; j++, k++, m0 += nDim)
867 if (m_bitMask.IsValid(k))
868 {
869 if (j < hd.nCols - 1 && m_bitMask.IsValid(k + 1)) // hori
870 {
871 for (int s0 = 0, iDim = 0; iDim < nDim; iDim++, s0 += maxShift)
872 {
873 unsigned int c = ((unsigned int)data[m0 + iDim]) ^ ((unsigned int)data[m0 + iDim + nDim]);
874 AddUIntToCounts(&cntDiffVec[s0], c, maxShift);
875 }
876 cnt++;
877 }
878 if (i < hd.nRows - 1 && m_bitMask.IsValid(k + hd.nCols)) // vert
879 {
880 for (int s0 = 0, iDim = 0; iDim < nDim; iDim++, s0 += maxShift)
881 {
882 unsigned int c = ((unsigned int)data[m0 + iDim]) ^ ((unsigned int)data[m0 + iDim + nDim * hd.nCols]);
883 AddUIntToCounts(&cntDiffVec[s0], c, maxShift);
884 }
885 cnt++;
886 }
887 }
888 }
889 else if (hd.dt == DT_Char || hd.dt == DT_Short || hd.dt == DT_Int) // signed int
890 {
891 for (int k = 0, m0 = 0, i = 0; i < hd.nRows; i++)
892 for (int j = 0; j < hd.nCols; j++, k++, m0 += nDim)
893 if (m_bitMask.IsValid(k))
894 {
895 if (j < hd.nCols - 1 && m_bitMask.IsValid(k + 1)) // hori
896 {
897 for (int s0 = 0, iDim = 0; iDim < nDim; iDim++, s0 += maxShift)
898 {
899 int c = ((int)data[m0 + iDim]) ^ ((int)data[m0 + iDim + nDim]);
900 AddIntToCounts(&cntDiffVec[s0], c, maxShift);
901 }
902 cnt++;
903 }
904 if (i < hd.nRows - 1 && m_bitMask.IsValid(k + hd.nCols)) // vert
905 {
906 for (int s0 = 0, iDim = 0; iDim < nDim; iDim++, s0 += maxShift)
907 {
908 int c = ((int)data[m0 + iDim]) ^ ((int)data[m0 + iDim + nDim * hd.nCols]);
909 AddIntToCounts(&cntDiffVec[s0], c, maxShift);
910 }
911 cnt++;
912 }
913 }
914 }
915 else
916 return false; // unsupported data type
917 }
918
919 if (cnt < minCnt) // not enough data for good stats
920 return false;
921
922 int nCutFound = 0, lastPlaneKept = 0;
923 //const bool printAll = false;
924
925 for (int s = maxShift - 1; s >= 0; s--)
926 {
927 //if (printAll) printf("bit plane %2d: ", s);
928 bool bCrit = true;
929
930 for (int iDim = 0; iDim < nDim; iDim++)
931 {
932 double x = cntDiffVec[iDim * maxShift + s];
933 double n = cnt;
934 double m = x / n;
935 //double stdDev = sqrt(x * x / n - m * m) / n;
936
937 //printf(" %.4f +- %.4f ", (float)(2 * m), (float)(2 * stdDev));
938 //if (printAll) printf(" %.4f ", (float)(2 * m));
939
940 if (fabs(1 - 2 * m) >= eps)
941 bCrit = false;
942 }
943 //if (printAll) printf("\n");
944
945 if (bCrit && nCutFound < 2)
946 {
947 if (nCutFound == 0)
948 lastPlaneKept = s;
949
950 if (nCutFound == 1 && s < lastPlaneKept - 1)
951 {
952 lastPlaneKept = s;
953 nCutFound = 0;
954 //if (printAll) printf(" reset ");
955 }
956
957 nCutFound++;
958 //if (printAll && nCutFound == 1) printf("\n");
959 }
960 }
961
962 lastPlaneKept = std::max(0, lastPlaneKept);
963 //if (printAll) printf("%d \n", lastPlaneKept);
964
965 newMaxZError = (1 << lastPlaneKept) >> 1; // turn lastPlaneKept into new maxZError
966
967 return true;
968 }
969
970 // -------------------------------------------------------------------------- ;
971
972 template<class T>
TryRaiseMaxZError(const T * data,double & maxZError) const973 bool Lerc2::TryRaiseMaxZError(const T* data, double& maxZError) const
974 {
975 if (!data || m_headerInfo.dt < DT_Float || m_headerInfo.numValidPixel == 0)
976 return false;
977
978 const HeaderInfo& hd = m_headerInfo;
979 const int nDim = hd.nDim;
980
981 std::vector<double> roundErr, zErr, zErrCand = { 1, 0.5, 0.1, 0.05, 0.01, 0.005, 0.001, 0.0005, 0.0001 };
982 std::vector<int> zFac, zFacCand = { 1, 2, 10, 20, 100, 200, 1000, 2000, 10000 };
983
984 for (size_t i = 0; i < zErrCand.size(); i++)
985 if (zErrCand[i] / 2 > maxZError)
986 {
987 zErr.push_back(zErrCand[i] / 2);
988 zFac.push_back(zFacCand[i]);
989 roundErr.push_back(0);
990 }
991
992 if (zErr.empty())
993 return false;
994
995 if (nDim == 1 && hd.numValidPixel == hd.nCols * hd.nRows) // special but common case
996 {
997 for (int i = 0; i < hd.nRows; i++)
998 {
999 size_t nCand = zErr.size();
1000
1001 for (int k = i * hd.nCols, j = 0; j < hd.nCols; j++, k++)
1002 {
1003 double x = data[k];
1004
1005 for (size_t n = 0; n < nCand; n++)
1006 {
1007 double z = x * zFac[n];
1008 if (z == (int)z)
1009 break;
1010
1011 double delta = fabs(floor(z + 0.5) - z);
1012 roundErr[n] = std::max(roundErr[n], delta);
1013 }
1014 }
1015
1016 if (!PruneCandidates(roundErr, zErr, zFac, maxZError))
1017 return false;
1018 }
1019 }
1020
1021 else // general case: nDim > 1 or not all pixel valid
1022 {
1023 for (int k = 0, m0 = 0, i = 0; i < hd.nRows; i++)
1024 {
1025 size_t nCand = zErr.size();
1026
1027 for (int j = 0; j < hd.nCols; j++, k++, m0 += nDim)
1028 if (m_bitMask.IsValid(k))
1029 for (int iDim = 0; iDim < nDim; iDim++)
1030 {
1031 double x = data[m0 + iDim];
1032
1033 for (size_t n = 0; n < nCand; n++)
1034 {
1035 double z = x * zFac[n];
1036 if (z == (int)z)
1037 break;
1038
1039 double delta = fabs(floor(z + 0.5) - z);
1040 roundErr[n] = std::max(roundErr[n], delta);
1041 }
1042 }
1043
1044 if (!PruneCandidates(roundErr, zErr, zFac, maxZError))
1045 return false;
1046 }
1047 }
1048
1049 for (size_t n = 0; n < zErr.size(); n++)
1050 if (roundErr[n] / zFac[n] <= maxZError)
1051 {
1052 maxZError = zErr[n];
1053 return true;
1054 }
1055
1056 return false;
1057 }
1058
1059 // -------------------------------------------------------------------------- ;
1060
PruneCandidates(std::vector<double> & roundErr,std::vector<double> & zErr,std::vector<int> & zFac,double maxZError)1061 bool Lerc2::PruneCandidates(std::vector<double>& roundErr, std::vector<double>& zErr,
1062 std::vector<int>& zFac, double maxZError)
1063 {
1064 size_t nCand = zErr.size();
1065
1066 if (nCand == 0 || roundErr.size() != nCand || zFac.size() != nCand || maxZError <= 0)
1067 return false;
1068
1069 for (int n = (int)(nCand - 1); n >= 0; n--)
1070 if (roundErr[n] / zFac[n] > maxZError)
1071 {
1072 roundErr.erase(roundErr.begin() + n);
1073 zErr.erase(zErr.begin() + n);
1074 zFac.erase(zFac.begin() + n);
1075 }
1076
1077 return (!zErr.empty());
1078 }
1079
1080 // -------------------------------------------------------------------------- ;
1081
1082 template<class T>
WriteDataOneSweep(const T * data,Byte ** ppByte) const1083 bool Lerc2::WriteDataOneSweep(const T* data, Byte** ppByte) const
1084 {
1085 if (!data || !ppByte)
1086 return false;
1087
1088 Byte* ptr = (*ppByte);
1089 const HeaderInfo& hd = m_headerInfo;
1090 int nDim = hd.nDim;
1091 int len = nDim * sizeof(T);
1092
1093 for (int k = 0, m0 = 0, i = 0; i < hd.nRows; i++)
1094 for (int j = 0; j < hd.nCols; j++, k++, m0 += nDim)
1095 if (m_bitMask.IsValid(k))
1096 {
1097 memcpy(ptr, &data[m0], len);
1098 ptr += len;
1099 }
1100
1101 (*ppByte) = ptr;
1102 return true;
1103 }
1104
1105 // -------------------------------------------------------------------------- ;
1106
1107 template<class T>
ReadDataOneSweep(const Byte ** ppByte,size_t & nBytesRemaining,T * data) const1108 bool Lerc2::ReadDataOneSweep(const Byte** ppByte, size_t& nBytesRemaining, T* data) const
1109 {
1110 if (!data || !ppByte || !(*ppByte))
1111 return false;
1112
1113 const Byte* ptr = (*ppByte);
1114 const HeaderInfo& hd = m_headerInfo;
1115 int nDim = hd.nDim;
1116 int len = nDim * sizeof(T);
1117
1118 size_t nValidPix = (size_t)m_bitMask.CountValidBits();
1119
1120 if (nBytesRemaining < nValidPix * len)
1121 return false;
1122
1123 for (int k = 0, m0 = 0, i = 0; i < hd.nRows; i++)
1124 for (int j = 0; j < hd.nCols; j++, k++, m0 += nDim)
1125 if (m_bitMask.IsValid(k))
1126 {
1127 memcpy(&data[m0], ptr, len);
1128 ptr += len;
1129 }
1130
1131 (*ppByte) = ptr;
1132 nBytesRemaining -= nValidPix * len;
1133
1134 return true;
1135 }
1136
1137 // -------------------------------------------------------------------------- ;
1138
1139 template<class T>
ComputeMinMaxRanges(const T * data,std::vector<double> & zMinVecA,std::vector<double> & zMaxVecA) const1140 bool Lerc2::ComputeMinMaxRanges(const T* data, std::vector<double>& zMinVecA, std::vector<double>& zMaxVecA) const
1141 {
1142 if (!data || m_headerInfo.numValidPixel == 0)
1143 return false;
1144
1145 const HeaderInfo& hd = m_headerInfo;
1146 const int nDim = hd.nDim;
1147 bool bInit = false;
1148
1149 zMinVecA.resize(nDim);
1150 zMaxVecA.resize(nDim);
1151
1152 std::vector<T> zMinVec(nDim, 0), zMaxVec(nDim, 0);
1153
1154 if (hd.numValidPixel == hd.nRows * hd.nCols) // all valid, no mask
1155 {
1156 bInit = true;
1157 for (int iDim = 0; iDim < nDim; iDim++)
1158 zMinVec[iDim] = zMaxVec[iDim] = data[iDim];
1159
1160 for (int m0 = 0, i = 0; i < hd.nRows; i++)
1161 for (int j = 0; j < hd.nCols; j++, m0 += nDim)
1162 for (int iDim = 0; iDim < nDim; iDim++)
1163 {
1164 T val = data[m0 + iDim];
1165
1166 if (val < zMinVec[iDim])
1167 zMinVec[iDim] = val;
1168 else if (val > zMaxVec[iDim])
1169 zMaxVec[iDim] = val;
1170 }
1171 }
1172 else
1173 {
1174 for (int k = 0, m0 = 0, i = 0; i < hd.nRows; i++)
1175 for (int j = 0; j < hd.nCols; j++, k++, m0 += nDim)
1176 if (m_bitMask.IsValid(k))
1177 {
1178 if (bInit)
1179 for (int iDim = 0; iDim < nDim; iDim++)
1180 {
1181 T val = data[m0 + iDim];
1182
1183 if (val < zMinVec[iDim])
1184 zMinVec[iDim] = val;
1185 else if (val > zMaxVec[iDim])
1186 zMaxVec[iDim] = val;
1187 }
1188 else
1189 {
1190 bInit = true;
1191 for (int iDim = 0; iDim < nDim; iDim++)
1192 zMinVec[iDim] = zMaxVec[iDim] = data[m0 + iDim];
1193 }
1194 }
1195 }
1196
1197 if (bInit)
1198 for (int iDim = 0; iDim < nDim; iDim++)
1199 {
1200 zMinVecA[iDim] = zMinVec[iDim];
1201 zMaxVecA[iDim] = zMaxVec[iDim];
1202 }
1203
1204 return bInit;
1205 }
1206
1207 // -------------------------------------------------------------------------- ;
1208
1209 template<class T>
WriteTiles(const T * data,Byte ** ppByte,int & numBytes) const1210 bool Lerc2::WriteTiles(const T* data, Byte** ppByte, int& numBytes) const
1211 {
1212 if (!data || !ppByte)
1213 return false;
1214
1215 numBytes = 0;
1216 int numBytesLerc = 0;
1217
1218 std::vector<unsigned int> quantVec, quantVecDiff;
1219 std::vector<std::pair<unsigned int, unsigned int> > sortedQuantVec, sortedQuantVecDiff;
1220
1221 const HeaderInfo& hd = m_headerInfo;
1222 int mbSize = hd.microBlockSize;
1223 int nDim = hd.nDim;
1224
1225 std::vector<T> dataVec(mbSize * mbSize, 0);
1226 T* dataBuf = &dataVec[0];
1227
1228 const bool bDtInt = (hd.dt < DT_Float);
1229 const bool bIntLossless = bDtInt && (hd.maxZError == 0.5);
1230 const bool bTryDiffEnc = (hd.version >= 5) && (nDim > 1);
1231
1232 const bool bCheckForIntOverflow = NeedToCheckForIntOverflow(hd);
1233 const bool bCheckForFltRndErr = NeedToCheckForFltRndErr(hd);
1234
1235 int mbDiff2 = bTryDiffEnc ? mbSize * mbSize : 0;
1236 std::vector<int> diffDataVecInt(mbDiff2, 0); // use fixed type (int) for difference of all int types
1237 std::vector<T> diffDataVecFlt(mbDiff2, 0), prevDataVec(mbDiff2, 0);
1238
1239 int numTilesVert = (hd.nRows + mbSize - 1) / mbSize;
1240 int numTilesHori = (hd.nCols + mbSize - 1) / mbSize;
1241
1242 for (int iTile = 0; iTile < numTilesVert; iTile++)
1243 {
1244 int tileH = mbSize;
1245 int i0 = iTile * tileH;
1246 if (iTile == numTilesVert - 1)
1247 tileH = hd.nRows - i0;
1248
1249 for (int jTile = 0; jTile < numTilesHori; jTile++)
1250 {
1251 int tileW = mbSize;
1252 int j0 = jTile * tileW;
1253 if (jTile == numTilesHori - 1)
1254 tileW = hd.nCols - j0;
1255
1256 for (int iDim = 0; iDim < nDim; iDim++)
1257 {
1258 T zMin = 0, zMax = 0;
1259 int numValidPixel = 0;
1260 bool bQuantizeDone = false;
1261 bool tryLut = false;
1262
1263 if (!GetValidDataAndStats(data, i0, i0 + tileH, j0, j0 + tileW, iDim, dataBuf, zMin, zMax, numValidPixel, tryLut))
1264 return false;
1265
1266 if (numValidPixel == 0 && !(*ppByte))
1267 {
1268 numBytesLerc += nDim; // 1 byte per empty block
1269 break; // iDim loop
1270 }
1271
1272 //tryLut = false; // always OFF
1273 //tryLut = NeedToQuantize(numValidPixel, zMin, zMax); // always ON
1274
1275 // if needed, quantize the data here once
1276 if (((*ppByte && iDim == 0) || tryLut) && NeedToQuantize(numValidPixel, zMin, zMax))
1277 {
1278 Quantize(dataBuf, numValidPixel, zMin, quantVec);
1279 bQuantizeDone = true;
1280 if (tryLut)
1281 SortQuantArray(quantVec, sortedQuantVec);
1282 }
1283
1284 BlockEncodeMode blockEncodeMode(BEM_RawBinary), blockEncodeModeDiff(BEM_RawBinary);
1285 int numBytesNeeded = NumBytesTile(numValidPixel, zMin, zMax, hd.dt, tryLut, blockEncodeMode, sortedQuantVec);
1286
1287 int numBytesNeededDiff = numBytesNeeded + 1;
1288 int zMinDiffInt(0), zMaxDiffInt(0);
1289 T zMinDiffFlt(0), zMaxDiffFlt(0);
1290 double zMinDiff(0), zMaxDiff(0);
1291 bool bQuantizeDoneDiff = false, tryLutDiff = false;
1292
1293 if (bTryDiffEnc && iDim > 0 && numValidPixel > 0)
1294 {
1295 bool rv = bDtInt ? ComputeDiffSliceInt(&dataVec[0], &prevDataVec[0], numValidPixel, bCheckForIntOverflow, hd.maxZError, diffDataVecInt, zMinDiffInt, zMaxDiffInt, tryLutDiff)
1296 : ComputeDiffSliceFlt(&dataVec[0], &prevDataVec[0], numValidPixel, bCheckForFltRndErr, hd.maxZError, diffDataVecFlt, zMinDiffFlt, zMaxDiffFlt, tryLutDiff);
1297
1298 zMinDiff = bDtInt ? zMinDiffInt : zMinDiffFlt;
1299 zMaxDiff = bDtInt ? zMaxDiffInt : zMaxDiffFlt;
1300
1301 if (rv) // can be false due to int overflow, then fall back to abs / regular encode
1302 {
1303 // tryLutDiff = false; // for faster encode?
1304
1305 if (tryLutDiff && NeedToQuantize(numValidPixel, zMinDiff, zMaxDiff))
1306 {
1307 bDtInt ? Quantize(&diffDataVecInt[0], numValidPixel, zMinDiffInt, quantVecDiff)
1308 : Quantize(&diffDataVecFlt[0], numValidPixel, zMinDiffFlt, quantVecDiff);
1309 bQuantizeDoneDiff = true;
1310 SortQuantArray(quantVecDiff, sortedQuantVecDiff);
1311 }
1312
1313 int nb = bDtInt ? NumBytesTile(numValidPixel, zMinDiffInt, zMaxDiffInt, DT_Int, tryLutDiff, blockEncodeModeDiff, sortedQuantVecDiff)
1314 : NumBytesTile(numValidPixel, zMinDiffFlt, zMaxDiffFlt, hd.dt, tryLutDiff, blockEncodeModeDiff, sortedQuantVecDiff);
1315 if (nb > 0)
1316 numBytesNeededDiff = nb;
1317 }
1318 }
1319
1320 numBytesLerc += std::min(numBytesNeeded, numBytesNeededDiff);
1321
1322 if (bTryDiffEnc && iDim < (nDim - 1) && numValidPixel > 0)
1323 {
1324 if (iDim == 0)
1325 prevDataVec.resize(numValidPixel);
1326
1327 if (!bIntLossless) // if not int lossless: do lossy quantize back and forth, then copy to buffer
1328 {
1329 double zMaxClamp = m_zMaxVec[iDim];
1330 bool bClampScaleBack = (zMax + 2 * hd.maxZError > zMaxClamp);
1331
1332 if (iDim == 0 || numBytesNeeded <= numBytesNeededDiff) // for regular or abs encode
1333 {
1334 if (bQuantizeDone || NeedToQuantize(numValidPixel, zMin, zMax))
1335 {
1336 if (!bQuantizeDone)
1337 Quantize(dataBuf, numValidPixel, zMin, quantVec);
1338
1339 bQuantizeDone = true;
1340 ScaleBack(&prevDataVec[0], quantVec, zMin, false, bClampScaleBack, zMaxClamp, hd.maxZError);
1341 }
1342 else if (zMin == zMax || ((hd.maxZError > 0) && (0 == (unsigned int)(ComputeMaxVal(zMin, zMax, hd.maxZError) + 0.5)))) // const block case
1343 prevDataVec.assign(numValidPixel, zMin);
1344 else
1345 copy(dataVec.begin(), dataVec.begin() + numValidPixel, prevDataVec.begin());
1346 }
1347 else // for diff or relative encode
1348 {
1349 if (bQuantizeDoneDiff || NeedToQuantize(numValidPixel, zMinDiff, zMaxDiff))
1350 {
1351 if (!bQuantizeDoneDiff)
1352 bDtInt ? Quantize(&diffDataVecInt[0], numValidPixel, zMinDiffInt, quantVecDiff)
1353 : Quantize(&diffDataVecFlt[0], numValidPixel, zMinDiffFlt, quantVecDiff);
1354
1355 bQuantizeDoneDiff = true;
1356 ScaleBack(&prevDataVec[0], quantVecDiff, zMinDiff, true, bClampScaleBack, zMaxClamp, hd.maxZError);
1357 }
1358 else if (zMinDiff == zMaxDiff || ((hd.maxZError > 0) && (0 == (unsigned int)(ComputeMaxVal(zMinDiff, zMaxDiff, hd.maxZError) + 0.5)))) // const block case
1359 {
1360 ScaleBackConstBlock(&prevDataVec[0], numValidPixel, zMinDiff, bClampScaleBack, zMaxClamp);
1361 }
1362 else
1363 copy(dataVec.begin(), dataVec.begin() + numValidPixel, prevDataVec.begin());
1364 }
1365 }
1366 else
1367 copy(dataVec.begin(), dataVec.begin() + numValidPixel, prevDataVec.begin());
1368 }
1369
1370 if (*ppByte)
1371 {
1372 int numBytesWritten = 0;
1373 bool rv = false;
1374
1375 if (iDim == 0 || numBytesNeeded <= numBytesNeededDiff)
1376 {
1377 if (!bQuantizeDone && NeedToQuantize(numValidPixel, zMin, zMax))
1378 Quantize(dataBuf, numValidPixel, zMin, quantVec);
1379
1380 rv = WriteTile(dataBuf, numValidPixel, ppByte, numBytesWritten, j0, zMin, zMax, hd.dt, false, quantVec, blockEncodeMode, sortedQuantVec);
1381 }
1382 else
1383 {
1384 if (!bQuantizeDoneDiff && NeedToQuantize(numValidPixel, zMinDiff, zMaxDiff))
1385 {
1386 bDtInt ? Quantize(&diffDataVecInt[0], numValidPixel, zMinDiffInt, quantVecDiff)
1387 : Quantize(&diffDataVecFlt[0], numValidPixel, zMinDiffFlt, quantVecDiff);
1388 }
1389 rv = bDtInt
1390 ? WriteTile(&diffDataVecInt[0], numValidPixel, ppByte, numBytesWritten, j0, zMinDiffInt, zMaxDiffInt, DT_Int, true, quantVecDiff, blockEncodeModeDiff, sortedQuantVecDiff)
1391 : WriteTile(&diffDataVecFlt[0], numValidPixel, ppByte, numBytesWritten, j0, zMinDiffFlt, zMaxDiffFlt, hd.dt, true, quantVecDiff, blockEncodeModeDiff, sortedQuantVecDiff);
1392 }
1393
1394 if (!rv || numBytesWritten != std::min(numBytesNeeded, numBytesNeededDiff))
1395 return false;
1396 }
1397 }
1398 }
1399 }
1400
1401 numBytes += numBytesLerc;
1402 return true;
1403 }
1404
1405 // -------------------------------------------------------------------------- ;
1406
1407 template<class T>
ReadTiles(const Byte ** ppByte,size_t & nBytesRemaining,T * data) const1408 bool Lerc2::ReadTiles(const Byte** ppByte, size_t& nBytesRemaining, T* data) const
1409 {
1410 if (!data || !ppByte || !(*ppByte))
1411 return false;
1412
1413 std::vector<unsigned int> bufferVec;
1414
1415 const HeaderInfo& hd = m_headerInfo;
1416 int mbSize = hd.microBlockSize;
1417 int nDim = hd.nDim;
1418
1419 if (mbSize > 32) // fail gracefully in case of corrupted blob for old version <= 2 which had no checksum
1420 return false;
1421
1422 int numTilesVert = (hd.nRows + mbSize - 1) / mbSize;
1423 int numTilesHori = (hd.nCols + mbSize - 1) / mbSize;
1424
1425 for (int iTile = 0; iTile < numTilesVert; iTile++)
1426 {
1427 int tileH = mbSize;
1428 int i0 = iTile * tileH;
1429 if (iTile == numTilesVert - 1)
1430 tileH = hd.nRows - i0;
1431
1432 for (int jTile = 0; jTile < numTilesHori; jTile++)
1433 {
1434 int tileW = mbSize;
1435 int j0 = jTile * tileW;
1436 if (jTile == numTilesHori - 1)
1437 tileW = hd.nCols - j0;
1438
1439 for (int iDim = 0; iDim < nDim; iDim++)
1440 {
1441 if (!ReadTile(ppByte, nBytesRemaining, data, i0, i0 + tileH, j0, j0 + tileW, iDim, bufferVec))
1442 return false;
1443 }
1444 }
1445 }
1446
1447 return true;
1448 }
1449
1450 // -------------------------------------------------------------------------- ;
1451
1452 template<class T>
GetValidDataAndStats(const T * data,int i0,int i1,int j0,int j1,int iDim,T * dataBuf,T & zMin,T & zMax,int & numValidPixel,bool & tryLut) const1453 bool Lerc2::GetValidDataAndStats(const T* data, int i0, int i1, int j0, int j1, int iDim,
1454 T* dataBuf, T& zMin, T& zMax, int& numValidPixel, bool& tryLut) const
1455 {
1456 const HeaderInfo& hd = m_headerInfo;
1457
1458 if (!data || i0 < 0 || j0 < 0 || i1 > hd.nRows || j1 > hd.nCols || i0 >= i1 || j0 >= j1 || iDim < 0 || iDim > hd.nDim || !dataBuf)
1459 return false;
1460
1461 zMin = zMax = 0;
1462 tryLut = false;
1463
1464 T prevVal = 0;
1465 int cnt = 0, cntSameVal = 0;
1466 int nDim = hd.nDim;
1467
1468 if (hd.numValidPixel == hd.nCols * hd.nRows) // all valid, no mask
1469 {
1470 int k0 = i0 * hd.nCols + j0;
1471 int m0 = k0 * nDim + iDim;
1472 zMin = zMax = data[m0]; // init
1473
1474 for (int i = i0; i < i1; i++)
1475 {
1476 int k = i * hd.nCols + j0;
1477 int m = k * nDim + iDim;
1478
1479 for (int j = j0; j < j1; j++, m += nDim)
1480 {
1481 T val = data[m];
1482 dataBuf[cnt] = val;
1483
1484 if (val < zMin)
1485 zMin = val;
1486 else if (val > zMax)
1487 zMax = val;
1488
1489 if (val == prevVal)
1490 cntSameVal++;
1491
1492 prevVal = val;
1493 cnt++;
1494 }
1495 }
1496 }
1497 else // not all valid, use mask
1498 {
1499 for (int i = i0; i < i1; i++)
1500 {
1501 int k = i * hd.nCols + j0;
1502 int m = k * nDim + iDim;
1503
1504 for (int j = j0; j < j1; j++, k++, m += nDim)
1505 if (m_bitMask.IsValid(k))
1506 {
1507 T val = data[m];
1508 dataBuf[cnt] = val;
1509
1510 if (cnt > 0)
1511 {
1512 if (val < zMin)
1513 zMin = val;
1514 else if (val > zMax)
1515 zMax = val;
1516
1517 if (val == prevVal)
1518 cntSameVal++;
1519 }
1520 else
1521 zMin = zMax = val; // init
1522
1523 prevVal = val;
1524 cnt++;
1525 }
1526 }
1527 }
1528
1529 if (cnt > 4)
1530 tryLut = (zMax > zMin + 3 * hd.maxZError) && (2 * cntSameVal > cnt);
1531
1532 numValidPixel = cnt;
1533 return true;
1534 }
1535
1536 // -------------------------------------------------------------------------- ;
1537
1538 template<class T>
ComputeDiffSliceInt(const T * data,const T * prevData,int numValidPixel,bool bCheckForIntOverflow,double maxZError,std::vector<int> & diffDataVec,int & zMin,int & zMax,bool & tryLut)1539 bool Lerc2::ComputeDiffSliceInt(const T* data, const T* prevData, int numValidPixel,
1540 bool bCheckForIntOverflow, double maxZError,
1541 std::vector<int>& diffDataVec, int& zMin, int& zMax, bool& tryLut)
1542 {
1543 if (numValidPixel <= 0)
1544 return false;
1545
1546 diffDataVec.resize(numValidPixel);
1547
1548 int prevVal(0), cnt(0), cntSameVal(0);
1549
1550 if (!bCheckForIntOverflow)
1551 {
1552 zMin = zMax = (int)data[0] - (int)prevData[0]; // init
1553
1554 for (int i = 0; i < numValidPixel; i++)
1555 {
1556 int val = (int)data[i] - (int)prevData[i];
1557
1558 diffDataVec[i] = val;
1559
1560 if (val < zMin)
1561 zMin = val;
1562 else if (val > zMax)
1563 zMax = val;
1564
1565 if (val == prevVal)
1566 cntSameVal++;
1567
1568 prevVal = val;
1569 cnt++;
1570 }
1571 }
1572 else
1573 {
1574 zMin = zMax = (int)((double)data[0] - (double)prevData[0]); // init
1575 const double zIntMax = 0x7FFFFFFF;
1576 const double zIntMin = -zIntMax - 1;
1577 bool bOverflow = false;
1578
1579 for (int i = 0; i < numValidPixel; i++)
1580 {
1581 double z = (double)data[i] - (double)prevData[i];
1582 int val = (int)z;
1583
1584 if (z < zIntMin || z > zIntMax) // check for int32 overflow
1585 bOverflow = true;
1586
1587 diffDataVec[i] = val;
1588
1589 if (val < zMin)
1590 zMin = val;
1591 else if (val > zMax)
1592 zMax = val;
1593
1594 if (val == prevVal)
1595 cntSameVal++;
1596
1597 prevVal = val;
1598 cnt++;
1599 }
1600
1601 if (bOverflow)
1602 return false;
1603 }
1604
1605 if (cnt > 4)
1606 tryLut = (zMax > zMin + 3 * maxZError) && (2 * cntSameVal > cnt);
1607
1608 return true;
1609 }
1610
1611 // -------------------------------------------------------------------------- ;
1612
1613 template<class T>
ComputeDiffSliceFlt(const T * data,const T * prevData,int numValidPixel,bool bCheckForFltRndErr,double maxZError,std::vector<T> & diffDataVec,T & zMin,T & zMax,bool & tryLut)1614 bool Lerc2::ComputeDiffSliceFlt(const T* data, const T* prevData, int numValidPixel,
1615 bool bCheckForFltRndErr, double maxZError,
1616 std::vector<T>& diffDataVec, T& zMin, T& zMax, bool& tryLut)
1617 {
1618 if (numValidPixel <= 0)
1619 return false;
1620
1621 diffDataVec.resize(numValidPixel);
1622
1623 zMin = zMax = (T)((double)data[0] - (double)prevData[0]); // init
1624 T prevVal(0);
1625 int cnt(0), cntSameVal(0);
1626
1627 if (!bCheckForFltRndErr)
1628 {
1629 for (int i = 0; i < numValidPixel; i++)
1630 {
1631 T val = (T)((double)data[i] - (double)prevData[i]);
1632
1633 diffDataVec[i] = val;
1634
1635 if (val < zMin)
1636 zMin = val;
1637 else if (val > zMax)
1638 zMax = val;
1639
1640 if (val == prevVal)
1641 cntSameVal++;
1642
1643 prevVal = val;
1644 cnt++;
1645 }
1646 }
1647 else
1648 {
1649 double maxRoundErr = 0;
1650
1651 for (int i = 0; i < numValidPixel; i++)
1652 {
1653 T val = (T)((double)data[i] - (double)prevData[i]);
1654
1655 double testVal = (double)val + (double)prevData[i]; // compute flt pnt rounding error
1656 maxRoundErr = (std::max)(fabs(testVal - (double)data[i]), maxRoundErr);
1657
1658 diffDataVec[i] = val;
1659
1660 if (val < zMin)
1661 zMin = val;
1662 else if (val > zMax)
1663 zMax = val;
1664
1665 if (val == prevVal)
1666 cntSameVal++;
1667
1668 prevVal = val;
1669 cnt++;
1670 }
1671
1672 if (maxRoundErr > maxZError / 8)
1673 return false;
1674 }
1675
1676 if (cnt > 4)
1677 tryLut = (zMax > zMin + 3 * maxZError) && (2 * cntSameVal > cnt);
1678
1679 return true;
1680 }
1681
1682 // -------------------------------------------------------------------------- ;
1683
1684 template<class T>
WriteTile(const T * dataBuf,int num,Byte ** ppByte,int & numBytesWritten,int j0,T zMin,T zMax,DataType dtZ,bool bDiffEnc,const std::vector<unsigned int> & quantVec,BlockEncodeMode blockEncodeMode,const std::vector<std::pair<unsigned int,unsigned int>> & sortedQuantVec) const1685 bool Lerc2::WriteTile(const T* dataBuf, int num, Byte** ppByte, int& numBytesWritten, int j0, T zMin, T zMax,
1686 DataType dtZ, bool bDiffEnc, const std::vector<unsigned int>& quantVec, BlockEncodeMode blockEncodeMode,
1687 const std::vector<std::pair<unsigned int, unsigned int> >& sortedQuantVec) const
1688 {
1689 Byte* ptr = *ppByte;
1690 Byte comprFlag = ((j0 >> 3) & 15) << 2; // use bits 2345 for integrity check
1691
1692 if (m_headerInfo.version >= 5)
1693 comprFlag = bDiffEnc ? (comprFlag | 4) : (comprFlag & (7 << 3)); // bit 2 now encodes diff encoding
1694
1695 if (num == 0 || (zMin == 0 && zMax == 0)) // special cases
1696 {
1697 *ptr++ = comprFlag | 2; // set compression flag to 2 to mark tile as constant 0
1698 numBytesWritten = 1;
1699 *ppByte = ptr;
1700 return true;
1701 }
1702
1703 if (blockEncodeMode == BEM_RawBinary)
1704 {
1705 if (bDiffEnc)
1706 return false; // doesn't make sense, should not happen
1707
1708 *ptr++ = comprFlag | 0; // write z's binary uncompressed
1709
1710 memcpy(ptr, dataBuf, num * sizeof(T));
1711 ptr += num * sizeof(T);
1712 }
1713 else
1714 {
1715 double maxVal = (m_headerInfo.maxZError > 0) ? ComputeMaxVal(zMin, zMax, m_headerInfo.maxZError) : 0;
1716
1717 // write z's as int arr bit stuffed
1718 unsigned int maxElem = (unsigned int)(maxVal + 0.5);
1719 if (maxElem == 0)
1720 comprFlag |= 3; // set compression flag to 3 to mark tile as constant zMin
1721 else
1722 comprFlag |= 1; // use bit stuffing
1723
1724 DataType dtReduced;
1725 int bits67 = ReduceDataType(zMin, dtZ, dtReduced);
1726 comprFlag |= bits67 << 6;
1727
1728 *ptr++ = comprFlag;
1729
1730 if (!WriteVariableDataType(&ptr, (double)zMin, dtReduced))
1731 return false;
1732
1733 if (maxElem > 0)
1734 {
1735 if ((int)quantVec.size() != num)
1736 return false;
1737
1738 if (blockEncodeMode == BEM_BitStuffSimple)
1739 {
1740 if (!m_bitStuffer2.EncodeSimple(&ptr, quantVec, m_headerInfo.version))
1741 return false;
1742 }
1743 else if (blockEncodeMode == BEM_BitStuffLUT)
1744 {
1745 if (!m_bitStuffer2.EncodeLut(&ptr, sortedQuantVec, m_headerInfo.version))
1746 return false;
1747 }
1748 else
1749 return false;
1750 }
1751 }
1752
1753 numBytesWritten = (int)(ptr - *ppByte);
1754 *ppByte = ptr;
1755 return true;
1756 }
1757
1758 // -------------------------------------------------------------------------- ;
1759
1760 template<class T>
ReadTile(const Byte ** ppByte,size_t & nBytesRemainingInOut,T * data,int i0,int i1,int j0,int j1,int iDim,std::vector<unsigned int> & bufferVec) const1761 bool Lerc2::ReadTile(const Byte** ppByte, size_t& nBytesRemainingInOut, T* data, int i0, int i1, int j0, int j1, int iDim,
1762 std::vector<unsigned int>& bufferVec) const
1763 {
1764 const Byte* ptr = *ppByte;
1765 size_t nBytesRemaining = nBytesRemainingInOut;
1766
1767 if (nBytesRemaining < 1)
1768 return false;
1769
1770 const HeaderInfo& hd = m_headerInfo;
1771 int nCols = hd.nCols;
1772 int nDim = hd.nDim;
1773
1774 Byte comprFlag = *ptr++;
1775 nBytesRemaining--;
1776
1777 const bool bDiffEnc = (hd.version >= 5) ? (comprFlag & 4) : false; // bit 2 now encodes diff or delta encoding
1778 const int pattern = (hd.version >= 5) ? 14 : 15;
1779
1780 if (((comprFlag >> 2) & pattern) != ((j0 >> 3) & pattern)) // use bits 2345 for integrity check
1781 return false;
1782
1783 if (bDiffEnc && iDim == 0)
1784 return false;
1785
1786 int bits67 = comprFlag >> 6;
1787 comprFlag &= 3;
1788
1789 if (comprFlag == 2) // entire tile is constant 0 (all the valid pixels)
1790 {
1791 for (int i = i0; i < i1; i++)
1792 {
1793 int k = i * nCols + j0;
1794 int m = k * nDim + iDim;
1795
1796 for (int j = j0; j < j1; j++, k++, m += nDim)
1797 if (m_bitMask.IsValid(k))
1798 data[m] = bDiffEnc ? data[m - 1] : 0;
1799 }
1800
1801 *ppByte = ptr;
1802 nBytesRemainingInOut = nBytesRemaining;
1803 return true;
1804 }
1805
1806 else if (comprFlag == 0) // read z's binary uncompressed
1807 {
1808 if (bDiffEnc)
1809 return false; // doesn't make sense, should not happen
1810
1811 const T* srcPtr = (const T*)ptr;
1812 int cnt = 0;
1813
1814 for (int i = i0; i < i1; i++)
1815 {
1816 int k = i * nCols + j0;
1817 int m = k * nDim + iDim;
1818
1819 for (int j = j0; j < j1; j++, k++, m += nDim)
1820 if (m_bitMask.IsValid(k))
1821 {
1822 if (nBytesRemaining < sizeof(T))
1823 return false;
1824
1825 data[m] = *srcPtr++;
1826 nBytesRemaining -= sizeof(T);
1827
1828 cnt++;
1829 }
1830 }
1831
1832 ptr += cnt * sizeof(T);
1833 }
1834 else
1835 {
1836 // read z's as int arr bit stuffed
1837 DataType dtUsed = GetDataTypeUsed((bDiffEnc && hd.dt < DT_Float) ? DT_Int : hd.dt, bits67);
1838 if (dtUsed == DT_Undefined)
1839 return false;
1840 size_t n = GetDataTypeSize(dtUsed);
1841 if (nBytesRemaining < n)
1842 return false;
1843
1844 double offset = ReadVariableDataType(&ptr, dtUsed);
1845 nBytesRemaining -= n;
1846
1847 double zMax = (hd.version >= 4 && nDim > 1) ? m_zMaxVec[iDim] : hd.zMax;
1848
1849 if (comprFlag == 3) // entire tile is constant zMin (all the valid pixels)
1850 {
1851 for (int i = i0; i < i1; i++)
1852 {
1853 int k = i * nCols + j0;
1854 int m = k * nDim + iDim;
1855
1856 if (!bDiffEnc)
1857 {
1858 T val = (T)offset;
1859 for (int j = j0; j < j1; j++, k++, m += nDim)
1860 if (m_bitMask.IsValid(k))
1861 data[m] = val;
1862 }
1863 else
1864 {
1865 for (int j = j0; j < j1; j++, k++, m += nDim)
1866 if (m_bitMask.IsValid(k))
1867 {
1868 double z = offset + data[m - 1];
1869 data[m] = (T)std::min(z, zMax);
1870 }
1871 }
1872 }
1873 }
1874 else
1875 {
1876 size_t maxElementCount = (i1 - i0) * (j1 - j0);
1877 if (!m_bitStuffer2.Decode(&ptr, nBytesRemaining, bufferVec, maxElementCount, hd.version))
1878 return false;
1879
1880 double invScale = 2 * hd.maxZError; // for int types this is int
1881 const unsigned int* srcPtr = bufferVec.data();
1882
1883 if (bufferVec.size() == maxElementCount) // all valid
1884 {
1885 for (int i = i0; i < i1; i++)
1886 {
1887 int k = i * nCols + j0;
1888 int m = k * nDim + iDim;
1889
1890 if (!bDiffEnc)
1891 {
1892 for (int j = j0; j < j1; j++, k++, m += nDim)
1893 {
1894 double z = offset + *srcPtr++ * invScale;
1895 data[m] = (T)std::min(z, zMax); // make sure we stay in the orig range
1896 }
1897 }
1898 else
1899 {
1900 for (int j = j0; j < j1; j++, k++, m += nDim)
1901 {
1902 double z = offset + *srcPtr++ * invScale + data[m - 1];
1903 data[m] = (T)std::min(z, zMax);
1904 }
1905 }
1906 }
1907 }
1908 else // not all valid
1909 {
1910 if (hd.version > 2)
1911 {
1912 for (int i = i0; i < i1; i++)
1913 {
1914 int k = i * nCols + j0;
1915 int m = k * nDim + iDim;
1916
1917 if (!bDiffEnc)
1918 {
1919 for (int j = j0; j < j1; j++, k++, m += nDim)
1920 if (m_bitMask.IsValid(k))
1921 {
1922 double z = offset + *srcPtr++ * invScale;
1923 data[m] = (T)std::min(z, zMax); // make sure we stay in the orig range
1924 }
1925 }
1926 else
1927 {
1928 for (int j = j0; j < j1; j++, k++, m += nDim)
1929 if (m_bitMask.IsValid(k))
1930 {
1931 double z = offset + *srcPtr++ * invScale + data[m - 1];
1932 data[m] = (T)std::min(z, zMax);
1933 }
1934 }
1935 }
1936 }
1937 else // fail gracefully in case of corrupted blob for old version <= 2 which had no checksum
1938 {
1939 size_t bufferVecIdx = 0;
1940
1941 for (int i = i0; i < i1; i++)
1942 {
1943 int k = i * nCols + j0;
1944 int m = k * nDim + iDim;
1945
1946 for (int j = j0; j < j1; j++, k++, m += nDim)
1947 if (m_bitMask.IsValid(k))
1948 {
1949 if (bufferVecIdx == bufferVec.size()) // fail gracefully in case of corrupted blob for old version <= 2 which had no checksum
1950 return false;
1951
1952 double z = offset + bufferVec[bufferVecIdx] * invScale;
1953 bufferVecIdx++;
1954 data[m] = (T)std::min(z, zMax); // make sure we stay in the orig range
1955 }
1956 }
1957 }
1958 }
1959 }
1960 }
1961
1962 *ppByte = ptr;
1963 nBytesRemainingInOut = nBytesRemaining;
1964 return true;
1965 }
1966
1967 // -------------------------------------------------------------------------- ;
1968
1969 template<class T>
GetDataType(T z)1970 Lerc2::DataType Lerc2::GetDataType(T z)
1971 {
1972 const std::type_info& ti = typeid(z);
1973
1974 if (ti == typeid(signed char)) return DT_Char;
1975 else if (ti == typeid(Byte)) return DT_Byte;
1976 else if (ti == typeid(short)) return DT_Short;
1977 else if (ti == typeid(unsigned short)) return DT_UShort;
1978 else if (ti == typeid(int) && sizeof(int) == 4) return DT_Int;
1979 else if (ti == typeid(long) && sizeof(long) == 4) return DT_Int;
1980 else if (ti == typeid(unsigned int) && sizeof(unsigned int) == 4) return DT_UInt;
1981 else if (ti == typeid(unsigned long) && sizeof(unsigned long) == 4) return DT_UInt;
1982 else if (ti == typeid(float)) return DT_Float;
1983 else if (ti == typeid(double)) return DT_Double;
1984 else
1985 return DT_Undefined;
1986 }
1987
1988 // -------------------------------------------------------------------------- ;
1989
SortQuantArray(const vector<unsigned int> & quantVec,vector<pair<unsigned int,unsigned int>> & sortedQuantVec)1990 void Lerc2::SortQuantArray(const vector<unsigned int>& quantVec, vector<pair<unsigned int, unsigned int> >& sortedQuantVec)
1991 {
1992 int numElem = (int)quantVec.size();
1993 sortedQuantVec.resize(numElem);
1994
1995 for (int i = 0; i < numElem; i++)
1996 sortedQuantVec[i] = pair<unsigned int, unsigned int>(quantVec[i], i);
1997
1998 std::sort(sortedQuantVec.begin(), sortedQuantVec.end(),
1999 [](const pair<unsigned int, unsigned int>& p0,
2000 const pair<unsigned int, unsigned int>& p1) { return p0.first < p1.first; });
2001 }
2002
2003 // -------------------------------------------------------------------------- ;
2004
2005 template<class T>
ComputeHuffmanCodes(const T * data,int & numBytes,ImageEncodeMode & imageEncodeMode,std::vector<std::pair<unsigned short,unsigned int>> & codes) const2006 void Lerc2::ComputeHuffmanCodes(const T* data, int& numBytes, ImageEncodeMode& imageEncodeMode, std::vector<std::pair<unsigned short, unsigned int> >& codes) const
2007 {
2008 std::vector<int> histo, deltaHisto;
2009 ComputeHistoForHuffman(data, histo, deltaHisto);
2010
2011 int nBytes0 = 0, nBytes1 = 0;
2012 double avgBpp0 = 0, avgBpp1 = 0;
2013 Huffman huffman0, huffman1;
2014
2015 if (m_headerInfo.version >= 4)
2016 {
2017 if (!huffman0.ComputeCodes(histo) || !huffman0.ComputeCompressedSize(histo, nBytes0, avgBpp0))
2018 nBytes0 = 0;
2019 }
2020
2021 if (!huffman1.ComputeCodes(deltaHisto) || !huffman1.ComputeCompressedSize(deltaHisto, nBytes1, avgBpp1))
2022 nBytes1 = 0;
2023
2024 if (nBytes0 > 0 && nBytes1 > 0) // regular case, pick the better of the two
2025 {
2026 imageEncodeMode = (nBytes0 <= nBytes1) ? IEM_Huffman : IEM_DeltaHuffman;
2027 codes = (nBytes0 <= nBytes1) ? huffman0.GetCodes() : huffman1.GetCodes();
2028 numBytes = (std::min)(nBytes0, nBytes1);
2029 }
2030 else if (nBytes0 == 0 && nBytes1 == 0) // rare case huffman cannot handle, fall back to tiling
2031 {
2032 imageEncodeMode = IEM_Tiling;
2033 codes.resize(0);
2034 numBytes = 0;
2035 }
2036 else // rare also, pick the valid one, the other is 0
2037 {
2038 imageEncodeMode = (nBytes0 > nBytes1) ? IEM_Huffman : IEM_DeltaHuffman;
2039 codes = (nBytes0 > nBytes1) ? huffman0.GetCodes() : huffman1.GetCodes();
2040 numBytes = (std::max)(nBytes0, nBytes1);
2041 }
2042 }
2043
2044 // -------------------------------------------------------------------------- ;
2045
2046 template<class T>
ComputeHistoForHuffman(const T * data,std::vector<int> & histo,std::vector<int> & deltaHisto) const2047 void Lerc2::ComputeHistoForHuffman(const T* data, std::vector<int>& histo, std::vector<int>& deltaHisto) const
2048 {
2049 histo.resize(256);
2050 deltaHisto.resize(256);
2051
2052 memset(&histo[0], 0, histo.size() * sizeof(int));
2053 memset(&deltaHisto[0], 0, deltaHisto.size() * sizeof(int));
2054
2055 int offset = (m_headerInfo.dt == DT_Char) ? 128 : 0;
2056 int height = m_headerInfo.nRows;
2057 int width = m_headerInfo.nCols;
2058 int nDim = m_headerInfo.nDim;
2059
2060 if (m_headerInfo.numValidPixel == width * height) // all valid
2061 {
2062 for (int iDim = 0; iDim < nDim; iDim++)
2063 {
2064 T prevVal = 0;
2065 for (int m = iDim, i = 0; i < height; i++)
2066 for (int j = 0; j < width; j++, m += nDim)
2067 {
2068 T val = data[m];
2069 T delta = val;
2070
2071 if (j > 0)
2072 delta -= prevVal; // use overflow
2073 else if (i > 0)
2074 delta -= data[m - width * nDim];
2075 else
2076 delta -= prevVal;
2077
2078 prevVal = val;
2079
2080 histo[offset + (int)val]++;
2081 deltaHisto[offset + (int)delta]++;
2082 }
2083 }
2084 }
2085 else // not all valid
2086 {
2087 for (int iDim = 0; iDim < nDim; iDim++)
2088 {
2089 T prevVal = 0;
2090 for (int k = 0, m = iDim, i = 0; i < height; i++)
2091 for (int j = 0; j < width; j++, k++, m += nDim)
2092 if (m_bitMask.IsValid(k))
2093 {
2094 T val = data[m];
2095 T delta = val;
2096
2097 if (j > 0 && m_bitMask.IsValid(k - 1))
2098 {
2099 delta -= prevVal; // use overflow
2100 }
2101 else if (i > 0 && m_bitMask.IsValid(k - width))
2102 {
2103 delta -= data[m - width * nDim];
2104 }
2105 else
2106 delta -= prevVal;
2107
2108 prevVal = val;
2109
2110 histo[offset + (int)val]++;
2111 deltaHisto[offset + (int)delta]++;
2112 }
2113 }
2114 }
2115 }
2116
2117 // -------------------------------------------------------------------------- ;
2118
2119 template<class T>
EncodeHuffman(const T * data,Byte ** ppByte) const2120 bool Lerc2::EncodeHuffman(const T* data, Byte** ppByte) const
2121 {
2122 if (!data || !ppByte)
2123 return false;
2124
2125 Huffman huffman;
2126 if (!huffman.SetCodes(m_huffmanCodes) || !huffman.WriteCodeTable(ppByte, m_headerInfo.version)) // header and code table
2127 return false;
2128
2129 int offset = (m_headerInfo.dt == DT_Char) ? 128 : 0;
2130 int height = m_headerInfo.nRows;
2131 int width = m_headerInfo.nCols;
2132 int nDim = m_headerInfo.nDim;
2133
2134 unsigned int* arr = (unsigned int*)(*ppByte);
2135 unsigned int* dstPtr = arr;
2136 int bitPos = 0;
2137
2138 if (m_imageEncodeMode == IEM_DeltaHuffman)
2139 {
2140 for (int iDim = 0; iDim < nDim; iDim++)
2141 {
2142 T prevVal = 0;
2143 for (int k = 0, m = iDim, i = 0; i < height; i++)
2144 for (int j = 0; j < width; j++, k++, m += nDim)
2145 if (m_bitMask.IsValid(k))
2146 {
2147 T val = data[m];
2148 T delta = val;
2149
2150 if (j > 0 && m_bitMask.IsValid(k - 1))
2151 {
2152 delta -= prevVal; // use overflow
2153 }
2154 else if (i > 0 && m_bitMask.IsValid(k - width))
2155 {
2156 delta -= data[m - width * nDim];
2157 }
2158 else
2159 delta -= prevVal;
2160
2161 prevVal = val;
2162
2163 // bit stuff the huffman code for this delta
2164 int kBin = offset + (int)delta;
2165 int len = m_huffmanCodes[kBin].first;
2166 if (len <= 0)
2167 return false;
2168
2169 unsigned int code = m_huffmanCodes[kBin].second;
2170
2171 if (32 - bitPos >= len)
2172 {
2173 if (bitPos == 0)
2174 *dstPtr = 0;
2175
2176 *dstPtr |= code << (32 - bitPos - len);
2177 bitPos += len;
2178 if (bitPos == 32)
2179 {
2180 bitPos = 0;
2181 dstPtr++;
2182 }
2183 }
2184 else
2185 {
2186 bitPos += len - 32;
2187 *dstPtr++ |= code >> bitPos;
2188 *dstPtr = code << (32 - bitPos);
2189 }
2190 }
2191 }
2192 }
2193
2194 else if (m_imageEncodeMode == IEM_Huffman)
2195 {
2196 for (int k = 0, m0 = 0, i = 0; i < height; i++)
2197 for (int j = 0; j < width; j++, k++, m0 += nDim)
2198 if (m_bitMask.IsValid(k))
2199 for (int m = 0; m < nDim; m++)
2200 {
2201 T val = data[m0 + m];
2202
2203 // bit stuff the huffman code for this val
2204 int kBin = offset + (int)val;
2205 int len = m_huffmanCodes[kBin].first;
2206 if (len <= 0)
2207 return false;
2208
2209 unsigned int code = m_huffmanCodes[kBin].second;
2210
2211 if (32 - bitPos >= len)
2212 {
2213 if (bitPos == 0)
2214 *dstPtr = 0;
2215
2216 *dstPtr |= code << (32 - bitPos - len);
2217 bitPos += len;
2218 if (bitPos == 32)
2219 {
2220 bitPos = 0;
2221 dstPtr++;
2222 }
2223 }
2224 else
2225 {
2226 bitPos += len - 32;
2227 *dstPtr++ |= code >> bitPos;
2228 *dstPtr = code << (32 - bitPos);
2229 }
2230 }
2231 }
2232
2233 else
2234 return false;
2235
2236 size_t numUInts = dstPtr - arr + (bitPos > 0 ? 1 : 0) + 1; // add one more as the decode LUT can read ahead
2237 *ppByte += numUInts * sizeof(unsigned int);
2238 return true;
2239 }
2240
2241 // -------------------------------------------------------------------------- ;
2242
2243 template<class T>
DecodeHuffman(const Byte ** ppByte,size_t & nBytesRemainingInOut,T * data) const2244 bool Lerc2::DecodeHuffman(const Byte** ppByte, size_t& nBytesRemainingInOut, T* data) const
2245 {
2246 if (!data || !ppByte || !(*ppByte))
2247 return false;
2248
2249 Huffman huffman;
2250 if (!huffman.ReadCodeTable(ppByte, nBytesRemainingInOut, m_headerInfo.version)) // header and code table
2251 return false;
2252
2253 int numBitsLUT = 0;
2254 if (!huffman.BuildTreeFromCodes(numBitsLUT))
2255 return false;
2256
2257 int offset = (m_headerInfo.dt == DT_Char) ? 128 : 0;
2258 int height = m_headerInfo.nRows;
2259 int width = m_headerInfo.nCols;
2260 int nDim = m_headerInfo.nDim;
2261
2262 const unsigned int* arr = (const unsigned int*)(*ppByte);
2263 const unsigned int* srcPtr = arr;
2264 int bitPos = 0;
2265 size_t nBytesRemaining = nBytesRemainingInOut;
2266
2267 if (m_headerInfo.numValidPixel == width * height) // all valid
2268 {
2269 if (m_imageEncodeMode == IEM_DeltaHuffman)
2270 {
2271 for (int iDim = 0; iDim < nDim; iDim++)
2272 {
2273 T prevVal = 0;
2274 for (int m = iDim, i = 0; i < height; i++)
2275 for (int j = 0; j < width; j++, m += nDim)
2276 {
2277 int val = 0;
2278 if (nBytesRemaining >= 4 * sizeof(unsigned int))
2279 {
2280 if (!huffman.DecodeOneValue_NoOverrunCheck(&srcPtr, nBytesRemaining, bitPos, numBitsLUT, val))
2281 return false;
2282 }
2283 else
2284 {
2285 if (!huffman.DecodeOneValue(&srcPtr, nBytesRemaining, bitPos, numBitsLUT, val))
2286 return false;
2287 }
2288
2289 T delta = (T)(val - offset);
2290
2291 if (j > 0)
2292 delta += prevVal; // use overflow
2293 else if (i > 0)
2294 delta += data[m - width * nDim];
2295 else
2296 delta += prevVal;
2297
2298 data[m] = delta;
2299 prevVal = delta;
2300 }
2301 }
2302 }
2303
2304 else if (m_imageEncodeMode == IEM_Huffman)
2305 {
2306 for (int k = 0, m0 = 0, i = 0; i < height; i++)
2307 for (int j = 0; j < width; j++, k++, m0 += nDim)
2308 for (int m = 0; m < nDim; m++)
2309 {
2310 int val = 0;
2311 if (nBytesRemaining >= 4 * sizeof(unsigned int))
2312 {
2313 if (!huffman.DecodeOneValue_NoOverrunCheck(&srcPtr, nBytesRemaining, bitPos, numBitsLUT, val))
2314 return false;
2315 }
2316 else
2317 {
2318 if (!huffman.DecodeOneValue(&srcPtr, nBytesRemaining, bitPos, numBitsLUT, val))
2319 return false;
2320 }
2321
2322 data[m0 + m] = (T)(val - offset);
2323 }
2324 }
2325
2326 else
2327 return false;
2328 }
2329
2330 else // not all valid
2331 {
2332 if (m_imageEncodeMode == IEM_DeltaHuffman)
2333 {
2334 for (int iDim = 0; iDim < nDim; iDim++)
2335 {
2336 T prevVal = 0;
2337 for (int k = 0, m = iDim, i = 0; i < height; i++)
2338 for (int j = 0; j < width; j++, k++, m += nDim)
2339 if (m_bitMask.IsValid(k))
2340 {
2341 int val = 0;
2342 if (nBytesRemaining >= 4 * sizeof(unsigned int))
2343 {
2344 if (!huffman.DecodeOneValue_NoOverrunCheck(&srcPtr, nBytesRemaining, bitPos, numBitsLUT, val))
2345 return false;
2346 }
2347 else
2348 {
2349 if (!huffman.DecodeOneValue(&srcPtr, nBytesRemaining, bitPos, numBitsLUT, val))
2350 return false;
2351 }
2352
2353 T delta = (T)(val - offset);
2354
2355 if (j > 0 && m_bitMask.IsValid(k - 1))
2356 {
2357 delta += prevVal; // use overflow
2358 }
2359 else if (i > 0 && m_bitMask.IsValid(k - width))
2360 {
2361 delta += data[m - width * nDim];
2362 }
2363 else
2364 delta += prevVal;
2365
2366 data[m] = delta;
2367 prevVal = delta;
2368 }
2369 }
2370 }
2371
2372 else if (m_imageEncodeMode == IEM_Huffman)
2373 {
2374 for (int k = 0, m0 = 0, i = 0; i < height; i++)
2375 for (int j = 0; j < width; j++, k++, m0 += nDim)
2376 if (m_bitMask.IsValid(k))
2377 for (int m = 0; m < nDim; m++)
2378 {
2379 int val = 0;
2380 if (nBytesRemaining >= 4 * sizeof(unsigned int))
2381 {
2382 if (!huffman.DecodeOneValue_NoOverrunCheck(&srcPtr, nBytesRemaining, bitPos, numBitsLUT, val))
2383 return false;
2384 }
2385 else
2386 {
2387 if (!huffman.DecodeOneValue(&srcPtr, nBytesRemaining, bitPos, numBitsLUT, val))
2388 return false;
2389 }
2390
2391 data[m0 + m] = (T)(val - offset);
2392 }
2393 }
2394
2395 else
2396 return false;
2397 }
2398
2399 size_t numUInts = srcPtr - arr + (bitPos > 0 ? 1 : 0) + 1; // add one more as the decode LUT can read ahead
2400 size_t len = numUInts * sizeof(unsigned int);
2401
2402 if (nBytesRemainingInOut < len)
2403 return false;
2404
2405 *ppByte += len;
2406 nBytesRemainingInOut -= len;
2407 return true;
2408 }
2409
2410 // -------------------------------------------------------------------------- ;
2411
2412 template<class T>
WriteMinMaxRanges(const T *,Byte ** ppByte) const2413 bool Lerc2::WriteMinMaxRanges(const T* /*data*/, Byte** ppByte) const
2414 {
2415 if (!ppByte || !(*ppByte))
2416 return false;
2417
2418 //printf("write min / max = %f %f\n", m_zMinVec[0], m_zMaxVec[0]);
2419
2420 int nDim = m_headerInfo.nDim;
2421 if (/* nDim < 2 || */ (int)m_zMinVec.size() != nDim || (int)m_zMaxVec.size() != nDim)
2422 return false;
2423
2424 std::vector<T> zVec(nDim);
2425 size_t len = nDim * sizeof(T);
2426
2427 for (int i = 0; i < nDim; i++)
2428 zVec[i] = (T)m_zMinVec[i];
2429
2430 memcpy(*ppByte, &zVec[0], len);
2431 (*ppByte) += len;
2432
2433 for (int i = 0; i < nDim; i++)
2434 zVec[i] = (T)m_zMaxVec[i];
2435
2436 memcpy(*ppByte, &zVec[0], len);
2437 (*ppByte) += len;
2438
2439 return true;
2440 }
2441
2442 // -------------------------------------------------------------------------- ;
2443
2444 template<class T>
ReadMinMaxRanges(const Byte ** ppByte,size_t & nBytesRemaining,const T *)2445 bool Lerc2::ReadMinMaxRanges(const Byte** ppByte, size_t& nBytesRemaining, const T* /*data*/)
2446 {
2447 if (!ppByte || !(*ppByte))
2448 return false;
2449
2450 int nDim = m_headerInfo.nDim;
2451
2452 m_zMinVec.resize(nDim);
2453 m_zMaxVec.resize(nDim);
2454
2455 std::vector<T> zVec(nDim);
2456 size_t len = nDim * sizeof(T);
2457
2458 if (nBytesRemaining < len || !memcpy(&zVec[0], *ppByte, len))
2459 return false;
2460
2461 (*ppByte) += len;
2462 nBytesRemaining -= len;
2463
2464 for (int i = 0; i < nDim; i++)
2465 m_zMinVec[i] = zVec[i];
2466
2467 if (nBytesRemaining < len || !memcpy(&zVec[0], *ppByte, len))
2468 return false;
2469
2470 (*ppByte) += len;
2471 nBytesRemaining -= len;
2472
2473 for (int i = 0; i < nDim; i++)
2474 m_zMaxVec[i] = zVec[i];
2475
2476 //printf("read min / max = %f %f\n", m_zMinVec[0], m_zMaxVec[0]);
2477
2478 return true;
2479 }
2480
2481 // -------------------------------------------------------------------------- ;
2482
2483 template<class T>
FillConstImage(T * data) const2484 bool Lerc2::FillConstImage(T* data) const
2485 {
2486 if (!data)
2487 return false;
2488
2489 const HeaderInfo& hd = m_headerInfo;
2490 int nCols = hd.nCols;
2491 int nRows = hd.nRows;
2492 int nDim = hd.nDim;
2493 T z0 = (T)hd.zMin;
2494
2495 if (nDim == 1)
2496 {
2497 for (int k = 0, i = 0; i < nRows; i++)
2498 for (int j = 0; j < nCols; j++, k++)
2499 if (m_bitMask.IsValid(k))
2500 data[k] = z0;
2501 }
2502 else
2503 {
2504 std::vector<T> zBufVec(nDim, z0);
2505
2506 if (hd.zMin != hd.zMax)
2507 {
2508 if ((int)m_zMinVec.size() != nDim)
2509 return false;
2510
2511 for (int m = 0; m < nDim; m++)
2512 zBufVec[m] = (T)m_zMinVec[m];
2513 }
2514
2515 int len = nDim * sizeof(T);
2516 for (int k = 0, m = 0, i = 0; i < nRows; i++)
2517 for (int j = 0; j < nCols; j++, k++, m += nDim)
2518 if (m_bitMask.IsValid(k))
2519 memcpy(&data[m], &zBufVec[0], len);
2520 }
2521
2522 return true;
2523 }
2524
2525 // -------------------------------------------------------------------------- ;
2526
2527