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 : #ifndef LERC2_H
25 : #define LERC2_H
26 :
27 : #include <algorithm>
28 : #include <cfloat>
29 : #include <cmath>
30 : #include <limits>
31 : #include <string>
32 : #include <typeinfo>
33 : #include "Defines.h"
34 : #include "BitMask.h"
35 : #include "BitStuffer2.h"
36 : #include "Huffman.h"
37 : #include "RLE.h"
38 :
39 : NAMESPACE_LERC_START
40 :
41 : /** Lerc2 v1
42 : *
43 : * -- allow for lossless compression of all common data types
44 : * -- avoid data type conversions and copies
45 : * -- optimized compression for segmented rasters (10-15x lossless)
46 : * -- micro block is 8x8 fixed, only gets doubled to 16x16 if bit rate < 1 bpp
47 : * -- cnt is replaced by bit mask
48 : * -- Lerc blob header has data range [min, max]
49 : * -- harden consistency checks to detect if the byte blob has been tampered with
50 : * -- drop support for big endian, this is legacy now
51 : *
52 : * Lerc2 v2
53 : *
54 : * -- add Huffman coding for better lossless compression of 8 bit data types Char, Byte
55 : *
56 : * Lerc2 v3
57 : *
58 : * -- add checksum for the entire byte blob, for more rigorous detection of compressed data corruption
59 : * -- for the main bit stuffing routine, use an extra uint buffer for guaranteed memory alignment
60 : * -- this also allows dropping the NumExtraBytesToAllocate functions
61 : *
62 : * Lerc2 v4
63 : *
64 : * -- allow array per pixel, nDim values per pixel. Such as RGB, complex number, or larger arrays per pixel
65 : * -- extend Huffman coding for 8 bit data types from delta only to trying both delta and orig
66 : * -- for integer data types, allow to drop bit planes containing only random noise
67 : *
68 : */
69 :
70 : class Lerc2
71 : {
72 : public:
73 : Lerc2();
74 : Lerc2(int nDim, int nCols, int nRows, const Byte* pMaskBits = nullptr); // valid / invalid bits as byte array
75 7323 : virtual ~Lerc2() {}
76 :
77 : bool SetEncoderToOldVersion(int version); // call this to encode compatible to an old decoder
78 :
79 : bool Set(int nDim, int nCols, int nRows, const Byte* pMaskBits = nullptr);
80 :
81 : template<class T>
82 : unsigned int ComputeNumBytesNeededToWrite(const T* arr, double maxZError, bool encodeMask);
83 :
84 : //static unsigned int MinNumBytesNeededToReadHeader();
85 :
86 : /// dst buffer already allocated; byte ptr is moved like a file pointer
87 : template<class T>
88 : bool Encode(const T* arr, Byte** ppByte);
89 :
90 : // data types supported by Lerc2
91 : enum DataType {DT_Char = 0, DT_Byte, DT_Short, DT_UShort, DT_Int, DT_UInt, DT_Float, DT_Double, DT_Undefined};
92 :
93 : struct HeaderInfo
94 : {
95 : int version;
96 : unsigned int checksum;
97 : int nRows,
98 : nCols,
99 : nDim,
100 : numValidPixel,
101 : microBlockSize,
102 : blobSize;
103 :
104 : DataType dt;
105 :
106 : double maxZError,
107 : zMin, // if nDim > 1, this is the overall range
108 : zMax;
109 :
110 14811 : void RawInit() { memset(this, 0, sizeof(struct HeaderInfo)); }
111 :
112 11928 : bool TryHuffman() const { return version > 1 && (dt == DT_Byte || dt == DT_Char) && maxZError == 0.5; }
113 : };
114 :
115 : static bool GetHeaderInfo(const Byte* pByte, size_t nBytesRemaining, struct HeaderInfo& headerInfo);
116 :
117 : /// dst buffer already allocated; byte ptr is moved like a file pointer
118 : template<class T>
119 : bool Decode(const Byte** ppByte, size_t& nBytesRemaining, T* arr, Byte* pMaskBits = nullptr); // if mask ptr is not 0, mask bits are returned (even if all valid or same as previous)
120 :
121 : private:
122 : static const int kCurrVersion = 4; // 2: added Huffman coding to 8 bit types DT_Char, DT_Byte;
123 : // 3: changed the bit stuffing to using a uint aligned buffer,
124 : // added Fletcher32 checksum
125 : // 4: allow nDim values per pixel
126 :
127 : enum ImageEncodeMode { IEM_Tiling = 0, IEM_DeltaHuffman, IEM_Huffman };
128 : enum BlockEncodeMode { BEM_RawBinary = 0, BEM_BitStuffSimple, BEM_BitStuffLUT };
129 :
130 : int m_microBlockSize,
131 : m_maxValToQuantize;
132 : BitMask m_bitMask;
133 : HeaderInfo m_headerInfo;
134 : BitStuffer2 m_bitStuffer2;
135 : bool m_encodeMask,
136 : m_writeDataOneSweep;
137 : ImageEncodeMode m_imageEncodeMode;
138 :
139 : std::vector<double> m_zMinVec, m_zMaxVec;
140 : std::vector<std::pair<unsigned short, unsigned int> > m_huffmanCodes; // <= 256 codes, 1.5 kB
141 :
142 : private:
143 18750 : static std::string FileKey() { return "Lerc2 "; }
144 17850 : static bool IsLittleEndianSystem() { int n = 1; return (1 == *((Byte*)&n)) && (4 == sizeof(int)); }
145 : void Init();
146 :
147 : static unsigned int ComputeNumBytesHeaderToWrite(const struct HeaderInfo& hd);
148 : static bool WriteHeader(Byte** ppByte, const struct HeaderInfo& hd);
149 : static bool ReadHeader(const Byte** ppByte, size_t& nBytesRemaining, struct HeaderInfo& hd);
150 :
151 : bool WriteMask(Byte** ppByte) const;
152 : bool ReadMask(const Byte** ppByte, size_t& nBytesRemaining);
153 :
154 : bool DoChecksOnEncode(Byte* pBlobBegin, Byte* pBlobEnd) const;
155 : static unsigned int ComputeChecksumFletcher32(const Byte* pByte, int len);
156 :
157 : static void AddUIntToCounts(int* pCounts, unsigned int val, int nBits);
158 : static void AddIntToCounts(int* pCounts, int val, int nBits);
159 :
160 : template<class T>
161 : bool TryBitPlaneCompression(const T* data, double eps, double& newMaxZError) const;
162 :
163 : template<class T>
164 : bool WriteDataOneSweep(const T* data, Byte** ppByte) const;
165 :
166 : template<class T>
167 : bool ReadDataOneSweep(const Byte** ppByte, size_t& nBytesRemaining, T* data) const;
168 :
169 : template<class T>
170 : bool WriteTiles(const T* data, Byte** ppByte, int& numBytes, std::vector<double>& zMinVec, std::vector<double>& zMaxVec) const;
171 :
172 : template<class T>
173 : bool ReadTiles(const Byte** ppByte, size_t& nBytesRemaining, T* data) const;
174 :
175 : template<class T>
176 : bool GetValidDataAndStats(const T* data, int i0, int i1, int j0, int j1, int iDim,
177 : T* dataBuf, T& zMinA, T& zMaxA, int& numValidPixel, bool& tryLutA) const;
178 :
179 : static double ComputeMaxVal(double zMin, double zMax, double maxZError);
180 :
181 : template<class T>
182 : bool NeedToQuantize(int numValidPixel, T zMin, T zMax) const;
183 :
184 : template<class T>
185 : bool Quantize(const T* dataBuf, int num, T zMin, std::vector<unsigned int>& quantVec) const;
186 :
187 : template<class T>
188 : int NumBytesTile(int numValidPixel, T zMin, T zMax, bool tryLut, BlockEncodeMode& blockEncodeMode,
189 : const std::vector<std::pair<unsigned int, unsigned int> >& sortedQuantVec) const;
190 :
191 : template<class T>
192 : bool WriteTile(const T* dataBuf, int num, Byte** ppByte, int& numBytesWritten, int j0, T zMin, T zMax,
193 : const std::vector<unsigned int>& quantVec, BlockEncodeMode blockEncodeMode,
194 : const std::vector<std::pair<unsigned int, unsigned int> >& sortedQuantVec) const;
195 :
196 : template<class T>
197 : bool ReadTile(const Byte** ppByte, size_t& nBytesRemaining, T* data, int i0, int i1, int j0, int j1, int iDim,
198 : std::vector<unsigned int>& bufferVec) const;
199 :
200 : template<class T>
201 : int TypeCode(T z, DataType& dtUsed) const;
202 :
203 : DataType GetDataTypeUsed(int typeCode) const;
204 :
205 : static DataType ValidateDataType(int dt);
206 :
207 : static bool WriteVariableDataType(Byte** ppByte, double z, DataType dtUsed);
208 :
209 : static double ReadVariableDataType(const Byte** ppByte, DataType dtUsed);
210 :
211 : template<class T> DataType GetDataType(T z) const;
212 :
213 : static unsigned int GetMaxValToQuantize(DataType dt);
214 :
215 : static unsigned int GetDataTypeSize(DataType dt);
216 :
217 : static void SortQuantArray(const std::vector<unsigned int>& quantVec,
218 : std::vector<std::pair<unsigned int, unsigned int> >& sortedQuantVec);
219 :
220 : template<class T>
221 : void ComputeHuffmanCodes(const T* data, int& numBytes, ImageEncodeMode& imageEncodeMode,
222 : std::vector<std::pair<unsigned short, unsigned int> >& codes) const;
223 :
224 : template<class T>
225 : void ComputeHistoForHuffman(const T* data, std::vector<int>& histo, std::vector<int>& deltaHisto) const;
226 :
227 : template<class T>
228 : bool EncodeHuffman(const T* data, Byte** ppByte) const;
229 :
230 : template<class T>
231 : bool DecodeHuffman(const Byte** ppByte, size_t& nBytesRemaining, T* data) const;
232 :
233 : template<class T>
234 : bool WriteMinMaxRanges(const T* data, Byte** ppByte) const;
235 :
236 : template<class T>
237 : bool ReadMinMaxRanges(const Byte** ppByte, size_t& nBytesRemaining, const T* data);
238 :
239 : bool CheckMinMaxRanges(bool& minMaxEqual) const;
240 :
241 : template<class T>
242 : bool FillConstImage(T* data) const;
243 : };
244 :
245 : // -------------------------------------------------------------------------- ;
246 : // -------------------------------------------------------------------------- ;
247 :
248 : template<class T>
249 5181 : unsigned int Lerc2::ComputeNumBytesNeededToWrite(const T* arr, double maxZError, bool encodeMask)
250 : {
251 5181 : if (!arr || !IsLittleEndianSystem())
252 0 : return 0;
253 :
254 : // header
255 5181 : unsigned int nBytesHeaderMask = ComputeNumBytesHeaderToWrite(m_headerInfo);
256 :
257 : // valid / invalid mask
258 5181 : int numValid = m_headerInfo.numValidPixel;
259 5181 : int numTotal = m_headerInfo.nCols * m_headerInfo.nRows;
260 :
261 5181 : bool needMask = numValid > 0 && numValid < numTotal;
262 :
263 5181 : m_encodeMask = encodeMask;
264 :
265 5181 : nBytesHeaderMask += 1 * sizeof(int); // the mask encode numBytes
266 :
267 5181 : if (needMask && encodeMask)
268 : {
269 9 : RLE rle;
270 9 : size_t n = rle.computeNumBytesRLE((const Byte*)m_bitMask.Bits(), m_bitMask.Size());
271 9 : nBytesHeaderMask += (unsigned int)n;
272 : }
273 :
274 5181 : m_headerInfo.dt = GetDataType(arr[0]);
275 :
276 5181 : if (m_headerInfo.dt == DT_Undefined)
277 0 : return 0;
278 :
279 5181 : if (maxZError == 777) // cheat code
280 0 : maxZError = -0.01;
281 :
282 5181 : if (m_headerInfo.dt < DT_Float) // integer types
283 : {
284 : // interpret a negative maxZError as bit plane epsilon; dflt = 0.01;
285 5128 : if (maxZError < 0 && (!TryBitPlaneCompression(arr, -maxZError, maxZError)))
286 0 : maxZError = 0;
287 :
288 5128 : maxZError = std::max(0.5, floor(maxZError));
289 : }
290 53 : else if (maxZError < 0) // don't allow bit plane compression for float or double yet
291 0 : return 0;
292 :
293 5181 : m_headerInfo.maxZError = maxZError;
294 5181 : m_headerInfo.zMin = 0;
295 5181 : m_headerInfo.zMax = 0;
296 5181 : m_headerInfo.microBlockSize = m_microBlockSize;
297 5181 : m_headerInfo.blobSize = nBytesHeaderMask;
298 :
299 5181 : if (numValid == 0)
300 9 : return nBytesHeaderMask;
301 :
302 5172 : m_maxValToQuantize = GetMaxValToQuantize(m_headerInfo.dt);
303 :
304 5172 : Byte* ptr = nullptr; // only emulate the writing and just count the bytes needed
305 5172 : int nBytesTiling = 0;
306 :
307 5172 : if (!WriteTiles(arr, &ptr, nBytesTiling, m_zMinVec, m_zMaxVec)) // also fills the min max ranges
308 0 : return 0;
309 :
310 5172 : m_headerInfo.zMin = *std::min_element(m_zMinVec.begin(), m_zMinVec.end());
311 5172 : m_headerInfo.zMax = *std::max_element(m_zMaxVec.begin(), m_zMaxVec.end());
312 :
313 5172 : if (m_headerInfo.zMin == m_headerInfo.zMax) // image is const
314 1699 : return nBytesHeaderMask;
315 :
316 3473 : int nDim = m_headerInfo.nDim;
317 :
318 3473 : if (m_headerInfo.version >= 4)
319 : {
320 : // add the min max ranges behind the mask and before the main data;
321 : // so we do not write it if no valid pixel or all same value const
322 342 : m_headerInfo.blobSize += 2 * nDim * sizeof(T);
323 :
324 342 : bool minMaxEqual = false;
325 342 : if (!CheckMinMaxRanges(minMaxEqual))
326 0 : return 0;
327 :
328 342 : if (minMaxEqual)
329 0 : return m_headerInfo.blobSize; // all nDim bands are const
330 : }
331 :
332 : // data
333 3473 : m_imageEncodeMode = IEM_Tiling;
334 3473 : int nBytesData = nBytesTiling;
335 3473 : int nBytesHuffman = 0;
336 :
337 3473 : if (m_headerInfo.TryHuffman())
338 : {
339 : ImageEncodeMode huffmanEncMode;
340 3327 : ComputeHuffmanCodes(arr, nBytesHuffman, huffmanEncMode, m_huffmanCodes); // save Huffman codes for later use
341 :
342 3327 : if (!m_huffmanCodes.empty() && nBytesHuffman < nBytesTiling)
343 : {
344 639 : m_imageEncodeMode = huffmanEncMode;
345 639 : nBytesData = nBytesHuffman;
346 : }
347 : else
348 2687 : m_huffmanCodes.resize(0);
349 : }
350 :
351 3472 : m_writeDataOneSweep = false;
352 3472 : int nBytesDataOneSweep = (int)(numValid * nDim * sizeof(T));
353 :
354 : {
355 : // try with double block size to reduce block header overhead, if
356 3472 : if ( (nBytesTiling * 8 < numTotal * nDim * 2) // resulting bit rate < x (2 bpp)
357 3154 : && (nBytesTiling < 4 * nBytesDataOneSweep) // bit stuffing is effective
358 3154 : && (nBytesHuffman == 0 || nBytesTiling < 2 * nBytesHuffman) ) // not much worse than huffman (otherwise huffman wins anyway)
359 : {
360 3154 : m_headerInfo.microBlockSize = m_microBlockSize * 2;
361 :
362 3154 : std::vector<double> zMinVec, zMaxVec;
363 3154 : int nBytes2 = 0;
364 3154 : if (!WriteTiles(arr, &ptr, nBytes2, zMinVec, zMaxVec)) // no huffman in here anymore
365 0 : return 0;
366 :
367 3154 : if (nBytes2 <= nBytesData)
368 : {
369 627 : nBytesData = nBytes2;
370 627 : m_imageEncodeMode = IEM_Tiling;
371 627 : m_huffmanCodes.resize(0);
372 : }
373 : else
374 : {
375 2527 : m_headerInfo.microBlockSize = m_microBlockSize; // reset to orig
376 : }
377 : }
378 : }
379 :
380 3472 : if (m_headerInfo.TryHuffman())
381 3327 : nBytesData += 1; // flag for image encode mode
382 :
383 3473 : if (nBytesDataOneSweep <= nBytesData)
384 : {
385 24 : m_writeDataOneSweep = true; // fallback: write data binary uncompressed in one sweep
386 24 : m_headerInfo.blobSize += 1 + nBytesDataOneSweep; // header, mask, min max ranges, flag, data one sweep
387 : }
388 : else
389 : {
390 3449 : m_writeDataOneSweep = false;
391 3449 : m_headerInfo.blobSize += 1 + nBytesData; // header, mask, min max ranges, flag(s), data
392 : }
393 :
394 3473 : return m_headerInfo.blobSize;
395 : }
396 :
397 : // -------------------------------------------------------------------------- ;
398 :
399 : template<class T>
400 5180 : bool Lerc2::Encode(const T* arr, Byte** ppByte)
401 : {
402 5180 : if (!arr || !ppByte || !IsLittleEndianSystem())
403 0 : return false;
404 :
405 5180 : Byte* ptrBlob = *ppByte; // keep a ptr to the start of the blob
406 :
407 5180 : if (!WriteHeader(ppByte, m_headerInfo))
408 0 : return false;
409 :
410 5180 : if (!WriteMask(ppByte))
411 0 : return false;
412 :
413 5180 : if (m_headerInfo.numValidPixel == 0 || m_headerInfo.zMin == m_headerInfo.zMax)
414 : {
415 1708 : return DoChecksOnEncode(ptrBlob, *ppByte);
416 : }
417 :
418 3472 : if (m_headerInfo.version >= 4)
419 : {
420 342 : if (!WriteMinMaxRanges(arr, ppByte))
421 0 : return false;
422 :
423 342 : bool minMaxEqual = false;
424 342 : if (!CheckMinMaxRanges(minMaxEqual))
425 0 : return false;
426 :
427 341 : if (minMaxEqual)
428 0 : return DoChecksOnEncode(ptrBlob, *ppByte);
429 : }
430 :
431 3471 : **ppByte = m_writeDataOneSweep ? 1 : 0; // write flag
432 3471 : (*ppByte)++;
433 :
434 3471 : if (!m_writeDataOneSweep)
435 : {
436 3448 : if (m_headerInfo.TryHuffman())
437 : {
438 3326 : **ppByte = (Byte)m_imageEncodeMode; // Huffman or tiling encode mode
439 3326 : (*ppByte)++;
440 :
441 3326 : if (!m_huffmanCodes.empty()) // Huffman, no tiling
442 : {
443 640 : if (m_imageEncodeMode != IEM_DeltaHuffman && m_imageEncodeMode != IEM_Huffman)
444 640 : return false;
445 :
446 640 : if (!EncodeHuffman(arr, ppByte)) // data bit stuffed
447 0 : return false;
448 :
449 640 : return DoChecksOnEncode(ptrBlob, *ppByte);
450 : }
451 : }
452 :
453 2809 : int numBytes = 0;
454 2809 : std::vector<double> zMinVec, zMaxVec;
455 2809 : if (!WriteTiles(arr, ppByte, numBytes, zMinVec, zMaxVec))
456 0 : return false;
457 : }
458 : else
459 : {
460 24 : if (!WriteDataOneSweep(arr, ppByte))
461 0 : return false;
462 : }
463 :
464 2833 : return DoChecksOnEncode(ptrBlob, *ppByte);
465 : }
466 :
467 : // -------------------------------------------------------------------------- ;
468 :
469 : template<class T>
470 2142 : bool Lerc2::Decode(const Byte** ppByte, size_t& nBytesRemaining, T* arr, Byte* pMaskBits)
471 : {
472 2142 : if (!arr || !ppByte || !IsLittleEndianSystem())
473 0 : return false;
474 :
475 2142 : const Byte* ptrBlob = *ppByte; // keep a ptr to the start of the blob
476 2142 : size_t nBytesRemaining00 = nBytesRemaining;
477 :
478 2142 : if (!ReadHeader(ppByte, nBytesRemaining, m_headerInfo))
479 0 : return false;
480 :
481 2142 : if (nBytesRemaining00 < (size_t)m_headerInfo.blobSize)
482 0 : return false;
483 :
484 2142 : if (m_headerInfo.version >= 3)
485 : {
486 531 : int nBytes = (int)(FileKey().length() + sizeof(int) + sizeof(unsigned int)); // start right after the checksum entry
487 531 : if (m_headerInfo.blobSize < nBytes)
488 0 : return false;
489 531 : unsigned int checksum = ComputeChecksumFletcher32(ptrBlob + nBytes, m_headerInfo.blobSize - nBytes);
490 : #ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
491 : // For fuzzing, ignore checksum verification
492 : (void)checksum;
493 : #else
494 531 : if (checksum != m_headerInfo.checksum)
495 0 : return false;
496 : #endif
497 : }
498 :
499 2142 : if (!ReadMask(ppByte, nBytesRemaining))
500 0 : return false;
501 :
502 2142 : if (pMaskBits) // return proper mask bits even if they were not stored
503 51 : memcpy(pMaskBits, m_bitMask.Bits(), m_bitMask.Size());
504 :
505 2142 : memset(arr, 0, m_headerInfo.nCols * m_headerInfo.nRows * m_headerInfo.nDim * sizeof(T));
506 :
507 2142 : if (m_headerInfo.numValidPixel == 0)
508 9 : return true;
509 :
510 2133 : if (m_headerInfo.zMin == m_headerInfo.zMax) // image is const
511 : {
512 576 : if (!FillConstImage(arr))
513 0 : return false;
514 :
515 576 : return true;
516 : }
517 :
518 1557 : if (m_headerInfo.version >= 4)
519 : {
520 505 : if (!ReadMinMaxRanges(ppByte, nBytesRemaining, arr))
521 0 : return false;
522 :
523 505 : bool minMaxEqual = false;
524 505 : if (!CheckMinMaxRanges(minMaxEqual))
525 0 : return false;
526 :
527 505 : if (minMaxEqual) // if all bands are const, fill outgoing and done
528 : {
529 0 : if (!FillConstImage(arr))
530 0 : return false;
531 :
532 0 : return true; // done
533 : }
534 : }
535 :
536 1557 : if (nBytesRemaining < 1)
537 0 : return false;
538 :
539 1557 : Byte readDataOneSweep = **ppByte; // read flag
540 1557 : (*ppByte)++;
541 1557 : nBytesRemaining--;
542 :
543 1557 : if (!readDataOneSweep)
544 : {
545 1533 : if (m_headerInfo.TryHuffman())
546 : {
547 1328 : if (nBytesRemaining < 1)
548 0 : return false;
549 :
550 1328 : Byte flag = **ppByte; // read flag Huffman / Lerc2
551 1328 : (*ppByte)++;
552 1328 : nBytesRemaining--;
553 :
554 1328 : if (flag > 2 || (m_headerInfo.version < 4 && flag > 1))
555 0 : return false;
556 :
557 1328 : m_imageEncodeMode = (ImageEncodeMode)flag;
558 :
559 1328 : if (m_imageEncodeMode == IEM_DeltaHuffman || m_imageEncodeMode == IEM_Huffman)
560 : {
561 409 : if (!DecodeHuffman(ppByte, nBytesRemaining, arr))
562 0 : return false;
563 :
564 409 : return true; // done.
565 : }
566 : }
567 :
568 1124 : if (!ReadTiles(ppByte, nBytesRemaining, arr))
569 0 : return false;
570 : }
571 : else
572 : {
573 24 : if (!ReadDataOneSweep(ppByte, nBytesRemaining, arr))
574 0 : return false;
575 : }
576 :
577 1148 : return true;
578 : }
579 :
580 : // -------------------------------------------------------------------------- ;
581 : // -------------------------------------------------------------------------- ;
582 :
583 : inline
584 0 : void Lerc2::AddUIntToCounts(int* pCounts, unsigned int val, int nBits)
585 : {
586 0 : pCounts[0] += val & 1;
587 0 : for (int i = 1; i < nBits; i++)
588 0 : pCounts[i] += (val >>= 1) & 1;
589 0 : }
590 :
591 : // -------------------------------------------------------------------------- ;
592 :
593 : inline
594 0 : void Lerc2::AddIntToCounts(int* pCounts, int val, int nBits)
595 : {
596 0 : pCounts[0] += val & 1;
597 0 : for (int i = 1; i < nBits; i++)
598 0 : pCounts[i] += (val >>= 1) & 1;
599 0 : }
600 :
601 : // -------------------------------------------------------------------------- ;
602 :
603 : // for the theory and math, see
604 : // https://pdfs.semanticscholar.org/d064/2e2ad1a4c3b445b0d795770f604a5d9e269c.pdf
605 :
606 : template<class T>
607 0 : bool Lerc2::TryBitPlaneCompression(const T* data, double eps, double& newMaxZError) const
608 : {
609 0 : newMaxZError = 0; // lossless is the obvious fallback
610 :
611 0 : if (!data || eps <= 0)
612 0 : return false;
613 :
614 0 : const HeaderInfo& hd = m_headerInfo;
615 0 : const int nDim = hd.nDim;
616 0 : const int maxShift = 8 * GetDataTypeSize(hd.dt);
617 0 : const int minCnt = 5000;
618 :
619 0 : if (hd.numValidPixel < minCnt) // not enough data for good stats
620 0 : return false;
621 :
622 0 : std::vector<int> cntDiffVec(nDim * maxShift, 0);
623 0 : int cnt = 0;
624 :
625 0 : if (nDim == 1 && hd.numValidPixel == hd.nCols * hd.nRows) // special but common case
626 : {
627 0 : if (hd.dt == DT_Byte || hd.dt == DT_UShort || hd.dt == DT_UInt) // unsigned int
628 : {
629 0 : for (int i = 0; i < hd.nRows - 1; i++)
630 0 : for (int k = i * hd.nCols, j = 0; j < hd.nCols - 1; j++, k++)
631 : {
632 0 : unsigned int c = ((unsigned int)data[k]) ^ ((unsigned int)data[k + 1]);
633 0 : AddUIntToCounts(&cntDiffVec[0], c, maxShift);
634 0 : cnt++;
635 0 : c = ((unsigned int)data[k]) ^ ((unsigned int)data[k + hd.nCols]);
636 0 : AddUIntToCounts(&cntDiffVec[0], c, maxShift);
637 0 : cnt++;
638 0 : }
639 : }
640 0 : else if (hd.dt == DT_Char || hd.dt == DT_Short || hd.dt == DT_Int) // signed int
641 : {
642 0 : for (int i = 0; i < hd.nRows - 1; i++)
643 0 : for (int k = i * hd.nCols, j = 0; j < hd.nCols - 1; j++, k++)
644 : {
645 0 : int c = ((int)data[k]) ^ ((int)data[k + 1]);
646 0 : AddIntToCounts(&cntDiffVec[0], c, maxShift);
647 0 : cnt++;
648 0 : c = ((int)data[k]) ^ ((int)data[k + hd.nCols]);
649 0 : AddIntToCounts(&cntDiffVec[0], c, maxShift);
650 0 : cnt++;
651 0 : }
652 : }
653 : else
654 0 : return false; // unsupported data type
655 : }
656 :
657 : else // general case: nDim > 1 or not all pixel valid
658 : {
659 0 : if (hd.dt == DT_Byte || hd.dt == DT_UShort || hd.dt == DT_UInt) // unsigned int
660 : {
661 0 : for (int k = 0, m0 = 0, i = 0; i < hd.nRows; i++)
662 0 : for (int j = 0; j < hd.nCols; j++, k++, m0 += nDim)
663 0 : if (m_bitMask.IsValid(k))
664 : {
665 0 : if (j < hd.nCols - 1 && m_bitMask.IsValid(k + 1)) // hori
666 : {
667 0 : for (int s0 = 0, iDim = 0; iDim < nDim; iDim++, s0 += maxShift)
668 : {
669 0 : unsigned int c = ((unsigned int)data[m0 + iDim]) ^ ((unsigned int)data[m0 + iDim + nDim]);
670 0 : AddUIntToCounts(&cntDiffVec[s0], c, maxShift);
671 : }
672 0 : cnt++;
673 : }
674 0 : if (i < hd.nRows - 1 && m_bitMask.IsValid(k + hd.nCols)) // vert
675 : {
676 0 : for (int s0 = 0, iDim = 0; iDim < nDim; iDim++, s0 += maxShift)
677 : {
678 0 : unsigned int c = ((unsigned int)data[m0 + iDim]) ^ ((unsigned int)data[m0 + iDim + nDim * hd.nCols]);
679 0 : AddUIntToCounts(&cntDiffVec[s0], c, maxShift);
680 : }
681 0 : cnt++;
682 : }
683 0 : }
684 : }
685 0 : else if (hd.dt == DT_Char || hd.dt == DT_Short || hd.dt == DT_Int) // signed int
686 : {
687 0 : for (int k = 0, m0 = 0, i = 0; i < hd.nRows; i++)
688 0 : for (int j = 0; j < hd.nCols; j++, k++, m0 += nDim)
689 0 : if (m_bitMask.IsValid(k))
690 : {
691 0 : if (j < hd.nCols - 1 && m_bitMask.IsValid(k + 1)) // hori
692 : {
693 0 : for (int s0 = 0, iDim = 0; iDim < nDim; iDim++, s0 += maxShift)
694 : {
695 0 : int c = ((int)data[m0 + iDim]) ^ ((int)data[m0 + iDim + nDim]);
696 0 : AddIntToCounts(&cntDiffVec[s0], c, maxShift);
697 : }
698 0 : cnt++;
699 : }
700 0 : if (i < hd.nRows - 1 && m_bitMask.IsValid(k + hd.nCols)) // vert
701 : {
702 0 : for (int s0 = 0, iDim = 0; iDim < nDim; iDim++, s0 += maxShift)
703 : {
704 0 : int c = ((int)data[m0 + iDim]) ^ ((int)data[m0 + iDim + nDim * hd.nCols]);
705 0 : AddIntToCounts(&cntDiffVec[s0], c, maxShift);
706 : }
707 0 : cnt++;
708 : }
709 0 : }
710 : }
711 : else
712 0 : return false; // unsupported data type
713 : }
714 :
715 0 : if (cnt < minCnt) // not enough data for good stats
716 0 : return false;
717 :
718 0 : int nCutFound = 0, lastPlaneKept = 0;
719 0 : const bool printAll = false;
720 :
721 0 : for (int s = maxShift - 1; s >= 0; s--)
722 : {
723 : if (printAll) printf("bit plane %2d: ", s);
724 0 : bool bCrit = true;
725 :
726 0 : for (int iDim = 0; iDim < nDim; iDim++)
727 : {
728 0 : double x = cntDiffVec[iDim * maxShift + s];
729 0 : double n = cnt;
730 0 : double m = x / n;
731 : //double stdDev = sqrt(x * x / n - m * m) / n;
732 :
733 : //printf(" %.4f +- %.4f ", (float)(2 * m), (float)(2 * stdDev));
734 : if (printAll) printf(" %.4f ", (float)(2 * m));
735 :
736 0 : if (fabs(1 - 2 * m) >= eps)
737 0 : bCrit = false;
738 : }
739 : if (printAll) printf("\n");
740 :
741 0 : if (bCrit && nCutFound < 2)
742 : {
743 0 : if (nCutFound == 0)
744 0 : lastPlaneKept = s;
745 :
746 0 : if (nCutFound == 1 && s < lastPlaneKept - 1)
747 : {
748 0 : lastPlaneKept = s;
749 0 : nCutFound = 0;
750 : if (printAll) printf(" reset ");
751 : }
752 :
753 0 : nCutFound++;
754 : if (printAll && nCutFound == 1) printf("\n");
755 : }
756 : }
757 :
758 0 : lastPlaneKept = std::max(0, lastPlaneKept);
759 : if (printAll) printf("%d \n", lastPlaneKept);
760 :
761 0 : newMaxZError = (1 << lastPlaneKept) >> 1; // turn lastPlaneKept into new maxZError
762 :
763 0 : return true;
764 : }
765 :
766 : // -------------------------------------------------------------------------- ;
767 :
768 : template<class T>
769 24 : bool Lerc2::WriteDataOneSweep(const T* data, Byte** ppByte) const
770 : {
771 24 : if (!data || !ppByte)
772 0 : return false;
773 :
774 24 : Byte* ptr = (*ppByte);
775 24 : const HeaderInfo& hd = m_headerInfo;
776 24 : int nDim = hd.nDim;
777 24 : int len = nDim * sizeof(T);
778 :
779 86 : for (int k = 0, m0 = 0, i = 0; i < hd.nRows; i++)
780 3492 : for (int j = 0; j < hd.nCols; j++, k++, m0 += nDim)
781 3430 : if (m_bitMask.IsValid(k))
782 : {
783 3430 : memcpy(ptr, &data[m0], len);
784 3430 : ptr += len;
785 : }
786 :
787 24 : (*ppByte) = ptr;
788 24 : return true;
789 : }
790 :
791 : // -------------------------------------------------------------------------- ;
792 :
793 : template<class T>
794 24 : bool Lerc2::ReadDataOneSweep(const Byte** ppByte, size_t& nBytesRemaining, T* data) const
795 : {
796 24 : if (!data || !ppByte || !(*ppByte))
797 0 : return false;
798 :
799 24 : const Byte* ptr = (*ppByte);
800 24 : const HeaderInfo& hd = m_headerInfo;
801 24 : int nDim = hd.nDim;
802 24 : int len = nDim * sizeof(T);
803 :
804 24 : size_t nValidPix = (size_t)m_bitMask.CountValidBits();
805 :
806 24 : if (nBytesRemaining < nValidPix * len)
807 0 : return false;
808 :
809 86 : for (int k = 0, m0 = 0, i = 0; i < hd.nRows; i++)
810 3492 : for (int j = 0; j < hd.nCols; j++, k++, m0 += nDim)
811 3430 : if (m_bitMask.IsValid(k))
812 : {
813 3430 : memcpy(&data[m0], ptr, len);
814 3430 : ptr += len;
815 : }
816 :
817 24 : (*ppByte) = ptr;
818 24 : nBytesRemaining -= nValidPix * len;
819 :
820 24 : return true;
821 : }
822 :
823 : // -------------------------------------------------------------------------- ;
824 :
825 : template<class T>
826 11135 : bool Lerc2::WriteTiles(const T* data, Byte** ppByte, int& numBytes, std::vector<double>& zMinVec, std::vector<double>& zMaxVec) const
827 : {
828 11135 : if (!data || !ppByte)
829 0 : return false;
830 :
831 11135 : numBytes = 0;
832 11135 : int numBytesLerc = 0;
833 :
834 22270 : std::vector<unsigned int> quantVec;
835 22270 : std::vector<std::pair<unsigned int, unsigned int> > sortedQuantVec;
836 :
837 11135 : const HeaderInfo& hd = m_headerInfo;
838 11135 : int mbSize = hd.microBlockSize;
839 11135 : int nDim = hd.nDim;
840 :
841 22271 : std::vector<T> dataVec(mbSize * mbSize, 0);
842 11135 : T* dataBuf = &dataVec[0];
843 :
844 11135 : zMinVec.assign(nDim, DBL_MAX);
845 11135 : zMaxVec.assign(nDim, -DBL_MAX);
846 :
847 11135 : int numTilesVert = (hd.nRows + mbSize - 1) / mbSize;
848 11135 : int numTilesHori = (hd.nCols + mbSize - 1) / mbSize;
849 :
850 155882 : for (int iTile = 0; iTile < numTilesVert; iTile++)
851 : {
852 144747 : int tileH = mbSize;
853 144747 : int i0 = iTile * tileH;
854 144747 : if (iTile == numTilesVert - 1)
855 11135 : tileH = hd.nRows - i0;
856 :
857 2319464 : for (int jTile = 0; jTile < numTilesHori; jTile++)
858 : {
859 2174719 : int tileW = mbSize;
860 2174719 : int j0 = jTile * tileW;
861 2174719 : if (jTile == numTilesHori - 1)
862 144747 : tileW = hd.nCols - j0;
863 :
864 4489554 : for (int iDim = 0; iDim < nDim; iDim++)
865 : {
866 2314845 : T zMin = 0, zMax = 0;
867 2314845 : int numValidPixel = 0;
868 2314845 : bool tryLut = false;
869 :
870 2314845 : if (!GetValidDataAndStats(data, i0, i0 + tileH, j0, j0 + tileW, iDim, dataBuf, zMin, zMax, numValidPixel, tryLut))
871 1 : return false;
872 :
873 2314835 : if (numValidPixel > 0)
874 : {
875 2314835 : zMinVec[iDim] = (std::min)(zMinVec[iDim], (double)zMin);
876 2314835 : zMaxVec[iDim] = (std::max)(zMaxVec[iDim], (double)zMax);
877 : }
878 :
879 : //tryLut = false;
880 :
881 : // if needed, quantize the data here once
882 2314835 : if ((*ppByte || tryLut) && NeedToQuantize(numValidPixel, zMin, zMax))
883 : {
884 728526 : if (!Quantize(dataBuf, numValidPixel, zMin, quantVec))
885 0 : return false;
886 :
887 728526 : if (tryLut)
888 709569 : SortQuantArray(quantVec, sortedQuantVec);
889 : }
890 :
891 : BlockEncodeMode blockEncodeMode;
892 2314835 : int numBytesNeeded = NumBytesTile(numValidPixel, zMin, zMax, tryLut, blockEncodeMode, sortedQuantVec);
893 2314845 : numBytesLerc += numBytesNeeded;
894 :
895 2314845 : if (*ppByte)
896 : {
897 629449 : int numBytesWritten = 0;
898 :
899 629449 : if (!WriteTile(dataBuf, numValidPixel, ppByte, numBytesWritten, j0, zMin, zMax, quantVec, blockEncodeMode, sortedQuantVec))
900 1 : return false;
901 :
902 629449 : if (numBytesWritten != numBytesNeeded)
903 1 : return false;
904 : }
905 : }
906 : }
907 : }
908 :
909 11135 : numBytes += numBytesLerc;
910 11135 : return true;
911 : }
912 :
913 : // -------------------------------------------------------------------------- ;
914 :
915 : template<class T>
916 1124 : bool Lerc2::ReadTiles(const Byte** ppByte, size_t& nBytesRemaining, T* data) const
917 : {
918 1124 : if (!data || !ppByte || !(*ppByte))
919 0 : return false;
920 :
921 2248 : std::vector<unsigned int> bufferVec;
922 :
923 1124 : const HeaderInfo& hd = m_headerInfo;
924 1124 : int mbSize = hd.microBlockSize;
925 1124 : int nDim = hd.nDim;
926 :
927 1124 : if (mbSize > 32) // fail gracefully in case of corrupted blob for old version <= 2 which had no checksum
928 0 : return false;
929 :
930 1124 : if( mbSize <= 0 || hd.nRows < 0 || hd.nCols < 0 ||
931 3372 : hd.nRows > std::numeric_limits<int>::max() - (mbSize - 1) ||
932 1124 : hd.nCols > std::numeric_limits<int>::max() - (mbSize - 1) )
933 : {
934 0 : return false;
935 : }
936 1124 : int numTilesVert = (hd.nRows + mbSize - 1) / mbSize;
937 1124 : int numTilesHori = (hd.nCols + mbSize - 1) / mbSize;
938 :
939 15523 : for (int iTile = 0; iTile < numTilesVert; iTile++)
940 : {
941 14399 : int tileH = mbSize;
942 14399 : int i0 = iTile * tileH;
943 14399 : if (iTile == numTilesVert - 1)
944 1124 : tileH = hd.nRows - i0;
945 :
946 251414 : for (int jTile = 0; jTile < numTilesHori; jTile++)
947 : {
948 237015 : int tileW = mbSize;
949 237015 : int j0 = jTile * tileW;
950 237015 : if (jTile == numTilesHori - 1)
951 14399 : tileW = hd.nCols - j0;
952 :
953 485152 : for (int iDim = 0; iDim < nDim; iDim++)
954 : {
955 248137 : if (!ReadTile(ppByte, nBytesRemaining, data, i0, i0 + tileH, j0, j0 + tileW, iDim, bufferVec))
956 0 : return false;
957 : }
958 : }
959 : }
960 :
961 1124 : return true;
962 : }
963 :
964 : // -------------------------------------------------------------------------- ;
965 :
966 : template<class T>
967 2314845 : bool Lerc2::GetValidDataAndStats(const T* data, int i0, int i1, int j0, int j1, int iDim,
968 : T* dataBuf, T& zMin, T& zMax, int& numValidPixel, bool& tryLut) const
969 : {
970 2314845 : const HeaderInfo& hd = m_headerInfo;
971 :
972 2314845 : if (!data || i0 < 0 || j0 < 0 || i1 > hd.nRows || j1 > hd.nCols || iDim < 0 || iDim > hd.nDim || !dataBuf)
973 1 : return false;
974 :
975 2314835 : zMin = 0;
976 2314835 : zMax = 0;
977 2314835 : tryLut = false;
978 :
979 2314835 : T prevVal = 0;
980 2314835 : int cnt = 0, cntSameVal = 0;
981 2314835 : int nDim = hd.nDim;
982 :
983 2314835 : if (hd.numValidPixel == hd.nCols * hd.nRows) // all valid, no mask
984 : {
985 23229467 : for (int i = i0; i < i1; i++)
986 : {
987 20914686 : int k = i * hd.nCols + j0;
988 20914686 : int m = k * nDim + iDim;
989 :
990 226878342 : for (int j = j0; j < j1; j++, k++, m += nDim)
991 : {
992 205962858 : T val = data[m];
993 205962858 : dataBuf[cnt] = val;
994 :
995 205962858 : if (cnt > 0)
996 : {
997 203649081 : if (val < zMin)
998 696827 : zMin = val;
999 202952350 : else if (val > zMax)
1000 714172 : zMax = val;
1001 :
1002 203649081 : if (val == prevVal)
1003 188847021 : cntSameVal++;
1004 : }
1005 : else
1006 2314651 : zMin = zMax = val; // init
1007 :
1008 205962858 : prevVal = val;
1009 205962858 : cnt++;
1010 : }
1011 : }
1012 : }
1013 : else // not all valid, use mask
1014 : {
1015 219 : for (int i = i0; i < i1; i++)
1016 : {
1017 110 : int k = i * hd.nCols + j0;
1018 110 : int m = k * nDim + iDim;
1019 :
1020 930 : for (int j = j0; j < j1; j++, k++, m += nDim)
1021 820 : if (m_bitMask.IsValid(k))
1022 : {
1023 413 : T val = data[m];
1024 413 : dataBuf[cnt] = val;
1025 :
1026 413 : if (cnt > 0)
1027 : {
1028 306 : if (val < zMin)
1029 0 : zMin = val;
1030 306 : else if (val > zMax)
1031 0 : zMax = val;
1032 :
1033 306 : if (val == prevVal)
1034 306 : cntSameVal++;
1035 : }
1036 : else
1037 107 : zMin = zMax = val; // init
1038 :
1039 413 : prevVal = val;
1040 413 : cnt++;
1041 : }
1042 : }
1043 : }
1044 :
1045 2314835 : if (cnt > 4)
1046 2314701 : tryLut = (zMax > zMin + hd.maxZError) && (2 * cntSameVal > cnt);
1047 :
1048 2314835 : numValidPixel = cnt;
1049 2314835 : return true;
1050 : }
1051 :
1052 : // -------------------------------------------------------------------------- ;
1053 :
1054 3559960 : inline double Lerc2::ComputeMaxVal(double zMin, double zMax, double maxZError)
1055 : {
1056 3559960 : double fac = 1 / (2 * maxZError); // must match the code in Decode(), don't touch it
1057 3559960 : return (zMax - zMin) * fac;
1058 : }
1059 :
1060 : // -------------------------------------------------------------------------- ;
1061 :
1062 : template<class T>
1063 1093786 : bool Lerc2::NeedToQuantize(int numValidPixel, T zMin, T zMax) const
1064 : {
1065 1093786 : if (numValidPixel == 0 || m_headerInfo.maxZError == 0)
1066 0 : return false;
1067 :
1068 1093786 : double maxVal = ComputeMaxVal(zMin, zMax, m_headerInfo.maxZError);
1069 1093786 : return !(maxVal > m_maxValToQuantize || (unsigned int)(maxVal + 0.5) == 0);
1070 : }
1071 :
1072 : // -------------------------------------------------------------------------- ;
1073 :
1074 : template<class T>
1075 728526 : bool Lerc2::Quantize(const T* dataBuf, int num, T zMin, std::vector<unsigned int>& quantVec) const
1076 : {
1077 728526 : quantVec.resize(num);
1078 :
1079 728526 : if (m_headerInfo.dt < DT_Float && m_headerInfo.maxZError == 0.5) // int lossless
1080 : {
1081 72097312 : for (int i = 0; i < num; i++)
1082 72097312 : quantVec[i] = (unsigned int)(dataBuf[i] - zMin); // ok, as char, short get promoted to int by C++ integral promotion rule
1083 : }
1084 : else // float and/or lossy
1085 : {
1086 15661 : double scale = 1 / (2 * m_headerInfo.maxZError);
1087 15661 : double zMinDbl = (double)zMin;
1088 :
1089 1062142 : for (int i = 0; i < num; i++)
1090 1046482 : quantVec[i] = (unsigned int)(((double)dataBuf[i] - zMinDbl) * scale + 0.5); // ok, consistent with ComputeMaxVal(...)
1091 : //quantVec[i] = (unsigned int)((dataBuf[i] - zMin) * scale + 0.5); // bad, not consistent with ComputeMaxVal(...)
1092 : }
1093 :
1094 728526 : return true;
1095 : }
1096 :
1097 : // -------------------------------------------------------------------------- ;
1098 :
1099 : template<class T>
1100 2314835 : int Lerc2::NumBytesTile(int numValidPixel, T zMin, T zMax, bool tryLut, BlockEncodeMode& blockEncodeMode,
1101 : const std::vector<std::pair<unsigned int, unsigned int> >& sortedQuantVec) const
1102 : {
1103 2314835 : blockEncodeMode = BEM_RawBinary;
1104 :
1105 2314835 : if (numValidPixel == 0 || (zMin == 0 && zMax == 0))
1106 376137 : return 1;
1107 :
1108 1938696 : double maxVal = 0, maxZError = m_headerInfo.maxZError;
1109 1938696 : int nBytesRaw = (int)(1 + numValidPixel * sizeof(T));
1110 :
1111 596 : if ((maxZError == 0 && zMax > zMin)
1112 1939292 : || (maxZError > 0 && (maxVal = ComputeMaxVal(zMin, zMax, maxZError)) > m_maxValToQuantize))
1113 : {
1114 436 : return nBytesRaw;
1115 : }
1116 : else
1117 : {
1118 : DataType dtUsed;
1119 1938260 : TypeCode(zMin, dtUsed);
1120 1938270 : int nBytes = 1 + GetDataTypeSize(dtUsed);
1121 :
1122 1938270 : unsigned int maxElem = (unsigned int)(maxVal + 0.5);
1123 1938270 : if (maxElem > 0)
1124 : {
1125 1480163 : nBytes += (!tryLut) ? m_bitStuffer2.ComputeNumBytesNeededSimple(numValidPixel, maxElem)
1126 709569 : : m_bitStuffer2.ComputeNumBytesNeededLut(sortedQuantVec, tryLut);
1127 : }
1128 :
1129 1938270 : if (nBytes < nBytesRaw)
1130 1918287 : blockEncodeMode = (!tryLut || maxElem == 0) ? BEM_BitStuffSimple : BEM_BitStuffLUT;
1131 : else
1132 19980 : nBytes = nBytesRaw;
1133 :
1134 1938270 : return nBytes;
1135 : }
1136 : }
1137 :
1138 : // -------------------------------------------------------------------------- ;
1139 :
1140 : template<class T>
1141 629449 : bool Lerc2::WriteTile(const T* dataBuf, int num, Byte** ppByte, int& numBytesWritten, int j0, T zMin, T zMax,
1142 : const std::vector<unsigned int>& quantVec, BlockEncodeMode blockEncodeMode,
1143 : const std::vector<std::pair<unsigned int, unsigned int> >& sortedQuantVec) const
1144 : {
1145 629449 : Byte* ptr = *ppByte;
1146 629449 : Byte comprFlag = ((j0 >> 3) & 15) << 2; // use bits 2345 for integrity check
1147 :
1148 629449 : if (num == 0 || (zMin == 0 && zMax == 0)) // special cases
1149 : {
1150 100003 : *ptr++ = comprFlag | 2; // set compression flag to 2 to mark tile as constant 0
1151 100003 : numBytesWritten = 1;
1152 100003 : *ppByte = ptr;
1153 100003 : return true;
1154 : }
1155 :
1156 529446 : if (blockEncodeMode == BEM_RawBinary)
1157 : {
1158 1369 : *ptr++ = comprFlag | 0; // write z's binary uncompressed
1159 :
1160 1369 : memcpy(ptr, dataBuf, num * sizeof(T));
1161 1369 : ptr += num * sizeof(T);
1162 : }
1163 : else
1164 : {
1165 528077 : double maxVal = (m_headerInfo.maxZError > 0) ? ComputeMaxVal(zMin, zMax, m_headerInfo.maxZError) : 0;
1166 :
1167 : // write z's as int arr bit stuffed
1168 528077 : unsigned int maxElem = (unsigned int)(maxVal + 0.5);
1169 528077 : if (maxElem == 0)
1170 265252 : comprFlag |= 3; // set compression flag to 3 to mark tile as constant zMin
1171 : else
1172 262825 : comprFlag |= 1; // use bit stuffing
1173 :
1174 : DataType dtUsed;
1175 528077 : int bits67 = TypeCode(zMin, dtUsed);
1176 528077 : comprFlag |= bits67 << 6;
1177 :
1178 528077 : *ptr++ = comprFlag;
1179 :
1180 528077 : if (!WriteVariableDataType(&ptr, (double)zMin, dtUsed))
1181 0 : return false;
1182 :
1183 528077 : if (maxElem > 0)
1184 : {
1185 262825 : if ((int)quantVec.size() != num)
1186 0 : return false;
1187 :
1188 262825 : if (blockEncodeMode == BEM_BitStuffSimple)
1189 : {
1190 41909 : if (!m_bitStuffer2.EncodeSimple(&ptr, quantVec, m_headerInfo.version))
1191 0 : return false;
1192 : }
1193 220916 : else if (blockEncodeMode == BEM_BitStuffLUT)
1194 : {
1195 220916 : if (!m_bitStuffer2.EncodeLut(&ptr, sortedQuantVec, m_headerInfo.version))
1196 0 : return false;
1197 : }
1198 : else
1199 0 : return false;
1200 : }
1201 : }
1202 :
1203 529446 : numBytesWritten = (int)(ptr - *ppByte);
1204 529446 : *ppByte = ptr;
1205 529446 : return true;
1206 : }
1207 :
1208 : // -------------------------------------------------------------------------- ;
1209 :
1210 : template<class T>
1211 248137 : bool Lerc2::ReadTile(const Byte** ppByte, size_t& nBytesRemainingInOut, T* data, int i0, int i1, int j0, int j1, int iDim,
1212 : std::vector<unsigned int>& bufferVec) const
1213 : {
1214 248137 : const Byte* ptr = *ppByte;
1215 248137 : size_t nBytesRemaining = nBytesRemainingInOut;
1216 :
1217 248137 : if (nBytesRemaining < 1)
1218 0 : return false;
1219 :
1220 248137 : Byte comprFlag = *ptr++;
1221 248137 : nBytesRemaining--;
1222 :
1223 248137 : int bits67 = comprFlag >> 6;
1224 248137 : int testCode = (comprFlag >> 2) & 15; // use bits 2345 for integrity check
1225 248137 : if (testCode != ((j0 >> 3) & 15))
1226 0 : return false;
1227 :
1228 248137 : const HeaderInfo& hd = m_headerInfo;
1229 248137 : int nCols = hd.nCols;
1230 248137 : int nDim = hd.nDim;
1231 :
1232 248137 : comprFlag &= 3;
1233 :
1234 248137 : if (comprFlag == 2) // entire tile is constant 0 (all the valid pixels)
1235 : {
1236 479349 : for (int i = i0; i < i1; i++)
1237 : {
1238 436220 : int k = i * nCols + j0;
1239 436220 : int m = k * nDim + iDim;
1240 :
1241 5408890 : for (int j = j0; j < j1; j++, k++, m += nDim)
1242 4972670 : if (m_bitMask.IsValid(k))
1243 4972670 : data[m] = 0;
1244 : }
1245 :
1246 43129 : *ppByte = ptr;
1247 43129 : nBytesRemainingInOut = nBytesRemaining;
1248 43129 : return true;
1249 : }
1250 :
1251 205008 : else if (comprFlag == 0) // read z's binary uncompressed
1252 : {
1253 1209 : const T* srcPtr = (const T*)ptr;
1254 1209 : int cnt = 0;
1255 :
1256 10713 : for (int i = i0; i < i1; i++)
1257 : {
1258 9504 : int k = i * nCols + j0;
1259 9504 : int m = k * nDim + iDim;
1260 :
1261 85536 : for (int j = j0; j < j1; j++, k++, m += nDim)
1262 76032 : if (m_bitMask.IsValid(k))
1263 : {
1264 76032 : if (nBytesRemaining < sizeof(T))
1265 0 : return false;
1266 :
1267 76032 : data[m] = *srcPtr++;
1268 76032 : nBytesRemaining -= sizeof(T);
1269 :
1270 76032 : cnt++;
1271 : }
1272 : }
1273 :
1274 1209 : ptr += cnt * sizeof(T);
1275 : }
1276 : else
1277 : {
1278 : // read z's as int arr bit stuffed
1279 203799 : DataType dtUsed = GetDataTypeUsed(bits67);
1280 203799 : if( dtUsed == DT_Undefined )
1281 0 : return false;
1282 203799 : size_t n = GetDataTypeSize(dtUsed);
1283 203799 : if (nBytesRemaining < n)
1284 0 : return false;
1285 :
1286 203799 : double offset = ReadVariableDataType(&ptr, dtUsed);
1287 203799 : nBytesRemaining -= n;
1288 :
1289 203799 : if (comprFlag == 3) // entire tile is constant zMin (all the valid pixels)
1290 : {
1291 905018 : for (int i = i0; i < i1; i++)
1292 : {
1293 814328 : int k = i * nCols + j0;
1294 814328 : int m = k * nDim + iDim;
1295 :
1296 8764730 : for (int j = j0; j < j1; j++, k++, m += nDim)
1297 7950400 : if (m_bitMask.IsValid(k))
1298 7950400 : data[m] = (T)offset;
1299 : }
1300 : }
1301 : else
1302 : {
1303 113109 : size_t maxElementCount = (i1 - i0) * (j1 - j0);
1304 113109 : if (!m_bitStuffer2.Decode(&ptr, nBytesRemaining, bufferVec, maxElementCount, hd.version))
1305 0 : return false;
1306 :
1307 113109 : double invScale = 2 * hd.maxZError; // for int types this is int
1308 113109 : double zMax = (hd.version >= 4 && nDim > 1) ? m_zMaxVec[iDim] : hd.zMax;
1309 113109 : const unsigned int* srcPtr = bufferVec.data();
1310 :
1311 113109 : if (bufferVec.size() == maxElementCount) // all valid
1312 : {
1313 1031152 : for (int i = i0; i < i1; i++)
1314 : {
1315 918040 : int k = i * nCols + j0;
1316 918040 : int m = k * nDim + iDim;
1317 :
1318 8486986 : for (int j = j0; j < j1; j++, k++, m += nDim)
1319 : {
1320 7568948 : double z = offset + *srcPtr++ * invScale;
1321 7568948 : data[m] = (T)std::min(z, zMax); // make sure we stay in the orig range
1322 : }
1323 : }
1324 : }
1325 : else // not all valid
1326 : {
1327 : #ifndef GDAL_COMPILATION
1328 : if (hd.version > 2)
1329 : {
1330 : for (int i = i0; i < i1; i++)
1331 : {
1332 : int k = i * nCols + j0;
1333 : int m = k * nDim + iDim;
1334 :
1335 : for (int j = j0; j < j1; j++, k++, m += nDim)
1336 : if (m_bitMask.IsValid(k))
1337 : {
1338 : double z = offset + *srcPtr++ * invScale;
1339 : data[m] = (T)std::min(z, zMax); // make sure we stay in the orig range
1340 : }
1341 : }
1342 : }
1343 : else // fail gracefully in case of corrupted blob for old version <= 2 which had no checksum
1344 : #endif
1345 : {
1346 0 : size_t bufferVecIdx = 0;
1347 :
1348 0 : for (int i = i0; i < i1; i++)
1349 : {
1350 0 : int k = i * nCols + j0;
1351 0 : int m = k * nDim + iDim;
1352 :
1353 0 : for (int j = j0; j < j1; j++, k++, m += nDim)
1354 0 : if (m_bitMask.IsValid(k))
1355 : {
1356 0 : if (bufferVecIdx == bufferVec.size()) // fail gracefully in case of corrupted blob for old version <= 2 which had no checksum
1357 0 : return false;
1358 :
1359 0 : double z = offset + bufferVec[bufferVecIdx] * invScale;
1360 0 : bufferVecIdx++;
1361 0 : data[m] = (T)std::min(z, zMax); // make sure we stay in the orig range
1362 : }
1363 : }
1364 : }
1365 : }
1366 : }
1367 : }
1368 :
1369 205008 : *ppByte = ptr;
1370 205008 : nBytesRemainingInOut = nBytesRemaining;
1371 205008 : return true;
1372 : }
1373 :
1374 : // -------------------------------------------------------------------------- ;
1375 :
1376 : template<class T>
1377 2466344 : int Lerc2::TypeCode(T z, DataType& dtUsed) const
1378 : {
1379 2466344 : Byte b = (Byte)z;
1380 2466344 : DataType dt = m_headerInfo.dt;
1381 2466344 : switch (dt)
1382 : {
1383 48 : case DT_Short:
1384 : {
1385 48 : signed char c = (signed char)z;
1386 48 : int tc = (T)c == z ? 2 : (T)b == z ? 1 : 0;
1387 48 : dtUsed = (DataType)(dt - tc);
1388 48 : return tc;
1389 : }
1390 48 : case DT_UShort:
1391 : {
1392 48 : int tc = (T)b == z ? 1 : 0;
1393 48 : dtUsed = (DataType)(dt - 2 * tc);
1394 48 : return tc;
1395 : }
1396 48 : case DT_Int:
1397 : {
1398 48 : short s = (short)z;
1399 48 : unsigned short us = (unsigned short)z;
1400 48 : int tc = (T)b == z ? 3 : (T)s == z ? 2 : (T)us == z ? 1 : 0;
1401 48 : dtUsed = (DataType)(dt - tc);
1402 48 : return tc;
1403 : }
1404 48 : case DT_UInt:
1405 : {
1406 48 : unsigned short us = (unsigned short)z;
1407 48 : int tc = (T)b == z ? 2 : (T)us == z ? 1 : 0;
1408 48 : dtUsed = (DataType)(dt - 2 * tc);
1409 48 : return tc;
1410 : }
1411 101 : case DT_Float:
1412 : {
1413 101 : short s = (short)z;
1414 101 : int tc = (T)b == z ? 2 : (T)s == z ? 1 : 0;
1415 101 : dtUsed = tc == 0 ? dt : (tc == 1 ? DT_Short : DT_Byte);
1416 101 : return tc;
1417 : }
1418 122 : case DT_Double:
1419 : {
1420 122 : short s = (short)z;
1421 122 : int l = (int)z;
1422 122 : float f = (float)z;
1423 122 : int tc = (T)s == z ? 3 : (T)l == z ? 2 : (T)f == z ? 1 : 0;
1424 122 : dtUsed = tc == 0 ? dt : (DataType)(dt - 2 * tc + 1);
1425 122 : return tc;
1426 : }
1427 2465929 : default:
1428 : {
1429 2465929 : dtUsed = dt;
1430 2465929 : return 0;
1431 : }
1432 : }
1433 : }
1434 :
1435 : // -------------------------------------------------------------------------- ;
1436 :
1437 60 : inline Lerc2::DataType Lerc2::ValidateDataType(int dt)
1438 : {
1439 60 : if( dt >= DT_Char && dt <= DT_Double )
1440 60 : return static_cast<DataType>(dt);
1441 0 : return DT_Undefined;
1442 : }
1443 :
1444 : // -------------------------------------------------------------------------- ;
1445 :
1446 : inline
1447 203799 : Lerc2::DataType Lerc2::GetDataTypeUsed(int tc) const
1448 : {
1449 203799 : DataType dt = m_headerInfo.dt;
1450 203799 : switch (dt)
1451 : {
1452 26 : case DT_Short:
1453 26 : case DT_Int: return ValidateDataType(dt - tc);
1454 26 : case DT_UShort:
1455 26 : case DT_UInt: return ValidateDataType(dt - 2 * tc);
1456 4 : case DT_Float: return tc == 0 ? dt : (tc == 1 ? DT_Short : DT_Byte);
1457 8 : case DT_Double: return tc == 0 ? dt : ValidateDataType(dt - 2 * tc + 1);
1458 203735 : default:
1459 203735 : return dt;
1460 : }
1461 : }
1462 :
1463 : // -------------------------------------------------------------------------- ;
1464 :
1465 : inline
1466 528077 : bool Lerc2::WriteVariableDataType(Byte** ppByte, double z, DataType dtUsed)
1467 : {
1468 528077 : Byte* ptr = *ppByte;
1469 :
1470 528077 : switch (dtUsed)
1471 : {
1472 13 : case DT_Char:
1473 : {
1474 13 : *((signed char*)ptr) = (signed char)z;
1475 13 : ptr++;
1476 13 : break;
1477 : }
1478 528056 : case DT_Byte:
1479 : {
1480 528056 : *((Byte*)ptr) = (Byte)z;
1481 528056 : ptr++;
1482 528056 : break;
1483 : }
1484 8 : case DT_Short:
1485 : {
1486 8 : short s = (short)z;
1487 8 : memcpy(ptr, &s, sizeof(short));
1488 8 : ptr += 2;
1489 8 : break;
1490 : }
1491 0 : case DT_UShort:
1492 : {
1493 0 : unsigned short us = (unsigned short)z;
1494 0 : memcpy(ptr, &us, sizeof(unsigned short));
1495 0 : ptr += 2;
1496 0 : break;
1497 : }
1498 0 : case DT_Int:
1499 : {
1500 0 : int i = (int)z;
1501 0 : memcpy(ptr, &i, sizeof(int));
1502 0 : ptr += 4;
1503 0 : break;
1504 : }
1505 0 : case DT_UInt:
1506 : {
1507 0 : unsigned int n = (unsigned int)z;
1508 0 : memcpy(ptr, &n, sizeof(unsigned int));
1509 0 : ptr += 4;
1510 0 : break;
1511 : }
1512 0 : case DT_Float:
1513 : {
1514 0 : float f = (float)z;
1515 0 : memcpy(ptr, &f, sizeof(float));
1516 0 : ptr += 4;
1517 0 : break;
1518 : }
1519 0 : case DT_Double:
1520 : {
1521 0 : memcpy(ptr, &z, sizeof(double));
1522 0 : ptr += 8;
1523 0 : break;
1524 : }
1525 :
1526 0 : default:
1527 0 : return false;
1528 : }
1529 :
1530 528077 : *ppByte = ptr;
1531 528077 : return true;
1532 : }
1533 :
1534 : // -------------------------------------------------------------------------- ;
1535 :
1536 : inline
1537 203799 : double Lerc2::ReadVariableDataType(const Byte** ppByte, DataType dtUsed)
1538 : {
1539 203799 : const Byte* ptr = *ppByte;
1540 :
1541 203799 : switch (dtUsed)
1542 : {
1543 13 : case DT_Char:
1544 : {
1545 13 : signed char c = *((signed char*)ptr);
1546 13 : *ppByte = ptr + 1;
1547 13 : return c;
1548 : }
1549 203778 : case DT_Byte:
1550 : {
1551 203778 : Byte b = *((Byte*)ptr);
1552 203778 : *ppByte = ptr + 1;
1553 203778 : return b;
1554 : }
1555 8 : case DT_Short:
1556 : {
1557 : short s;
1558 8 : memcpy(&s, ptr, sizeof(short));
1559 8 : *ppByte = ptr + 2;
1560 8 : return s;
1561 : }
1562 0 : case DT_UShort:
1563 : {
1564 : unsigned short us;
1565 0 : memcpy(&us, ptr, sizeof(unsigned short));
1566 0 : *ppByte = ptr + 2;
1567 0 : return us;
1568 : }
1569 0 : case DT_Int:
1570 : {
1571 : int i;
1572 0 : memcpy(&i, ptr, sizeof(int));
1573 0 : *ppByte = ptr + 4;
1574 0 : return i;
1575 : }
1576 0 : case DT_UInt:
1577 : {
1578 : unsigned int n;
1579 0 : memcpy(&n, ptr, sizeof(unsigned int));
1580 0 : *ppByte = ptr + 4;
1581 0 : return n;
1582 : }
1583 0 : case DT_Float:
1584 : {
1585 : float f;
1586 0 : memcpy(&f, ptr, sizeof(float));
1587 0 : *ppByte = ptr + 4;
1588 0 : return f;
1589 : }
1590 0 : case DT_Double:
1591 : {
1592 : double d;
1593 0 : memcpy(&d, ptr, sizeof(double));
1594 0 : *ppByte = ptr + 8;
1595 0 : return d;
1596 : }
1597 0 : default:
1598 0 : return 0;
1599 : }
1600 : }
1601 :
1602 : // -------------------------------------------------------------------------- ;
1603 :
1604 : template<class T>
1605 5181 : Lerc2::DataType Lerc2::GetDataType(T z) const
1606 : {
1607 5181 : const std::type_info& ti = typeid(z);
1608 :
1609 5181 : if (ti == typeid(signed char)) return DT_Char;
1610 5180 : else if (ti == typeid(Byte)) return DT_Byte;
1611 61 : else if (ti == typeid(short)) return DT_Short;
1612 59 : else if (ti == typeid(unsigned short)) return DT_UShort;
1613 57 : else if (ti == typeid(int ) && sizeof(int ) == 4) return DT_Int;
1614 55 : else if (ti == typeid(long) && sizeof(long) == 4) return DT_Int;
1615 55 : else if (ti == typeid(unsigned int ) && sizeof(unsigned int ) == 4) return DT_UInt;
1616 53 : else if (ti == typeid(unsigned long) && sizeof(unsigned long) == 4) return DT_UInt;
1617 53 : else if (ti == typeid(float)) return DT_Float;
1618 26 : else if (ti == typeid(double)) return DT_Double;
1619 : else
1620 0 : return DT_Undefined;
1621 : }
1622 :
1623 : // -------------------------------------------------------------------------- ;
1624 :
1625 : inline
1626 5172 : unsigned int Lerc2::GetMaxValToQuantize(DataType dt)
1627 : {
1628 5172 : switch (dt)
1629 : {
1630 5123 : case DT_Char:
1631 : case DT_Byte: //return (1 << 7) - 1; // disabled: allow LUT mode for 8 bit segmented
1632 : case DT_Short:
1633 5123 : case DT_UShort: return (1 << 15) - 1;
1634 :
1635 49 : case DT_Int:
1636 : case DT_UInt:
1637 : case DT_Float:
1638 49 : case DT_Double: return (1 << 30) - 1;
1639 :
1640 0 : default:
1641 0 : return 0;
1642 : }
1643 : }
1644 :
1645 : // -------------------------------------------------------------------------- ;
1646 :
1647 : inline
1648 2142070 : unsigned int Lerc2::GetDataTypeSize(DataType dt)
1649 : {
1650 2142070 : switch (dt)
1651 : {
1652 2141860 : case DT_Char:
1653 2141860 : case DT_Byte: return 1;
1654 42 : case DT_Short:
1655 42 : case DT_UShort: return 2;
1656 160 : case DT_Int:
1657 : case DT_UInt:
1658 160 : case DT_Float: return 4;
1659 0 : case DT_Double: return 8;
1660 :
1661 0 : default:
1662 0 : return 0;
1663 : }
1664 : }
1665 :
1666 : // -------------------------------------------------------------------------- ;
1667 :
1668 : template<class T>
1669 3327 : void Lerc2::ComputeHuffmanCodes(const T* data, int& numBytes, ImageEncodeMode& imageEncodeMode, std::vector<std::pair<unsigned short, unsigned int> >& codes) const
1670 : {
1671 6654 : std::vector<int> histo, deltaHisto;
1672 3327 : ComputeHistoForHuffman(data, histo, deltaHisto);
1673 :
1674 3327 : int nBytes0 = 0, nBytes1 = 0;
1675 3327 : double avgBpp0 = 0, avgBpp1 = 0;
1676 6654 : Huffman huffman0, huffman1;
1677 :
1678 3327 : if (m_headerInfo.version >= 4)
1679 : {
1680 204 : if (!huffman0.ComputeCodes(histo) || !huffman0.ComputeCompressedSize(histo, nBytes0, avgBpp0))
1681 0 : nBytes0 = 0;
1682 : }
1683 :
1684 3327 : if (!huffman1.ComputeCodes(deltaHisto) || !huffman1.ComputeCompressedSize(deltaHisto, nBytes1, avgBpp1))
1685 0 : nBytes1 = 0;
1686 :
1687 3327 : if (nBytes0 > 0 && nBytes1 > 0) // regular case, pick the better of the two
1688 : {
1689 204 : imageEncodeMode = (nBytes0 <= nBytes1) ? IEM_Huffman : IEM_DeltaHuffman;
1690 204 : codes = (nBytes0 <= nBytes1) ? huffman0.GetCodes() : huffman1.GetCodes();
1691 204 : numBytes = (std::min)(nBytes0, nBytes1);
1692 : }
1693 3123 : else if (nBytes0 == 0 && nBytes1 == 0) // rare case huffman cannot handle, fall back to tiling
1694 : {
1695 0 : imageEncodeMode = IEM_Tiling;
1696 0 : codes.resize(0);
1697 0 : numBytes = 0;
1698 : }
1699 : else // rare also, pick the valid one, the other is 0
1700 : {
1701 3123 : imageEncodeMode = (nBytes0 > nBytes1) ? IEM_Huffman : IEM_DeltaHuffman;
1702 3123 : codes = (nBytes0 > nBytes1) ? huffman0.GetCodes() : huffman1.GetCodes();
1703 3123 : numBytes = (std::max)(nBytes0, nBytes1);
1704 : }
1705 3327 : }
1706 :
1707 : // -------------------------------------------------------------------------- ;
1708 :
1709 : template<class T>
1710 3327 : void Lerc2::ComputeHistoForHuffman(const T* data, std::vector<int>& histo, std::vector<int>& deltaHisto) const
1711 : {
1712 3327 : histo.resize(256);
1713 3327 : deltaHisto.resize(256);
1714 :
1715 3327 : memset(&histo[0], 0, histo.size() * sizeof(int));
1716 3327 : memset(&deltaHisto[0], 0, deltaHisto.size() * sizeof(int));
1717 :
1718 3327 : int offset = (m_headerInfo.dt == DT_Char) ? 128 : 0;
1719 3327 : int height = m_headerInfo.nRows;
1720 3327 : int width = m_headerInfo.nCols;
1721 3327 : int nDim = m_headerInfo.nDim;
1722 :
1723 3327 : if (m_headerInfo.numValidPixel == width * height) // all valid
1724 : {
1725 6718 : for (int iDim = 0; iDim < nDim; iDim++)
1726 : {
1727 3391 : T prevVal = 0;
1728 423843 : for (int m = iDim, i = 0; i < height; i++)
1729 59177320 : for (int j = 0; j < width; j++, m += nDim)
1730 : {
1731 58756900 : T val = data[m];
1732 58756900 : T delta = val;
1733 :
1734 58756900 : if (j > 0)
1735 58336680 : delta -= prevVal; // use overflow
1736 420261 : else if (i > 0)
1737 417061 : delta -= data[m - width * nDim];
1738 : else
1739 3200 : delta -= prevVal;
1740 :
1741 58756900 : prevVal = val;
1742 :
1743 58756900 : histo[offset + (int)val]++;
1744 58755700 : deltaHisto[offset + (int)delta]++;
1745 : }
1746 : }
1747 : }
1748 : else // not all valid
1749 : {
1750 0 : for (int iDim = 0; iDim < nDim; iDim++)
1751 : {
1752 0 : T prevVal = 0;
1753 0 : for (int k = 0, m = iDim, i = 0; i < height; i++)
1754 0 : for (int j = 0; j < width; j++, k++, m += nDim)
1755 0 : if (m_bitMask.IsValid(k))
1756 : {
1757 0 : T val = data[m];
1758 0 : T delta = val;
1759 :
1760 0 : if (j > 0 && m_bitMask.IsValid(k - 1))
1761 : {
1762 0 : delta -= prevVal; // use overflow
1763 : }
1764 0 : else if (i > 0 && m_bitMask.IsValid(k - width))
1765 : {
1766 0 : delta -= data[m - width * nDim];
1767 : }
1768 : else
1769 0 : delta -= prevVal;
1770 :
1771 0 : prevVal = val;
1772 :
1773 0 : histo[offset + (int)val]++;
1774 0 : deltaHisto[offset + (int)delta]++;
1775 : }
1776 : }
1777 : }
1778 3327 : }
1779 :
1780 : // -------------------------------------------------------------------------- ;
1781 :
1782 : template<class T>
1783 640 : bool Lerc2::EncodeHuffman(const T* data, Byte** ppByte) const
1784 : {
1785 640 : if (!data || !ppByte)
1786 0 : return false;
1787 :
1788 1280 : Huffman huffman;
1789 640 : if (!huffman.SetCodes(m_huffmanCodes) || !huffman.WriteCodeTable(ppByte, m_headerInfo.version)) // header and code table
1790 0 : return false;
1791 :
1792 640 : int offset = (m_headerInfo.dt == DT_Char) ? 128 : 0;
1793 640 : int height = m_headerInfo.nRows;
1794 640 : int width = m_headerInfo.nCols;
1795 640 : int nDim = m_headerInfo.nDim;
1796 :
1797 640 : unsigned int* arr = (unsigned int*)(*ppByte);
1798 640 : unsigned int* dstPtr = arr;
1799 640 : int bitPos = 0;
1800 :
1801 640 : if (m_imageEncodeMode == IEM_DeltaHuffman)
1802 : {
1803 1039 : for (int iDim = 0; iDim < nDim; iDim++)
1804 : {
1805 534 : T prevVal = 0;
1806 62570 : for (int k = 0, m = iDim, i = 0; i < height; i++)
1807 8220130 : for (int j = 0; j < width; j++, k++, m += nDim)
1808 8158090 : if (m_bitMask.IsValid(k))
1809 : {
1810 8157910 : T val = data[m];
1811 8157910 : T delta = val;
1812 :
1813 8157910 : if (j > 0 && m_bitMask.IsValid(k - 1))
1814 : {
1815 8095680 : delta -= prevVal; // use overflow
1816 : }
1817 62052 : else if (i > 0 && m_bitMask.IsValid(k - width))
1818 : {
1819 61501 : delta -= data[m - width * nDim];
1820 : }
1821 : else
1822 551 : delta -= prevVal;
1823 :
1824 8157730 : prevVal = val;
1825 :
1826 : // bit stuff the huffman code for this delta
1827 8157730 : int kBin = offset + (int)delta;
1828 8157730 : int len = m_huffmanCodes[kBin].first;
1829 8157550 : if (len <= 0)
1830 0 : return false;
1831 :
1832 8157550 : unsigned int code = m_huffmanCodes[kBin].second;
1833 :
1834 8158100 : if (32 - bitPos >= len)
1835 : {
1836 8040170 : if (bitPos == 0)
1837 254850 : *dstPtr = 0;
1838 :
1839 8040170 : *dstPtr |= code << (32 - bitPos - len);
1840 8040170 : bitPos += len;
1841 8040170 : if (bitPos == 32)
1842 : {
1843 254357 : bitPos = 0;
1844 254357 : dstPtr++;
1845 : }
1846 : }
1847 : else
1848 : {
1849 117934 : bitPos += len - 32;
1850 117934 : *dstPtr++ |= code >> bitPos;
1851 117934 : *dstPtr = code << (32 - bitPos);
1852 : }
1853 : }
1854 : }
1855 : }
1856 :
1857 135 : else if (m_imageEncodeMode == IEM_Huffman)
1858 : {
1859 5583 : for (int k = 0, m0 = 0, i = 0; i < height; i++)
1860 1155140 : for (int j = 0; j < width; j++, k++, m0 += nDim)
1861 1149700 : if (m_bitMask.IsValid(k))
1862 2299390 : for (int m = 0; m < nDim; m++)
1863 : {
1864 1149700 : T val = data[m0 + m];
1865 :
1866 : // bit stuff the huffman code for this val
1867 1149700 : int kBin = offset + (int)val;
1868 1149700 : int len = m_huffmanCodes[kBin].first;
1869 1149700 : if (len <= 0)
1870 0 : return false;
1871 :
1872 1149700 : unsigned int code = m_huffmanCodes[kBin].second;
1873 :
1874 1149700 : if (32 - bitPos >= len)
1875 : {
1876 1013540 : if (bitPos == 0)
1877 36755 : *dstPtr = 0;
1878 :
1879 1013540 : *dstPtr |= code << (32 - bitPos - len);
1880 1013540 : bitPos += len;
1881 1013540 : if (bitPos == 32)
1882 : {
1883 36620 : bitPos = 0;
1884 36620 : dstPtr++;
1885 : }
1886 : }
1887 : else
1888 : {
1889 136154 : bitPos += len - 32;
1890 136154 : *dstPtr++ |= code >> bitPos;
1891 136154 : *dstPtr = code << (32 - bitPos);
1892 : }
1893 : }
1894 : }
1895 :
1896 : else
1897 0 : return false;
1898 :
1899 640 : size_t numUInts = dstPtr - arr + (bitPos > 0 ? 1 : 0) + 1; // add one more as the decode LUT can read ahead
1900 640 : *ppByte += numUInts * sizeof(unsigned int);
1901 640 : return true;
1902 : }
1903 :
1904 : // -------------------------------------------------------------------------- ;
1905 :
1906 : template<class T>
1907 409 : bool Lerc2::DecodeHuffman(const Byte** ppByte, size_t& nBytesRemainingInOut, T* data) const
1908 : {
1909 409 : if (!data || !ppByte || !(*ppByte))
1910 0 : return false;
1911 :
1912 1765 : Huffman huffman;
1913 409 : if (!huffman.ReadCodeTable(ppByte, nBytesRemainingInOut, m_headerInfo.version)) // header and code table
1914 0 : return false;
1915 :
1916 409 : int numBitsLUT = 0;
1917 409 : if (!huffman.BuildTreeFromCodes(numBitsLUT))
1918 0 : return false;
1919 :
1920 409 : int offset = (m_headerInfo.dt == DT_Char) ? 128 : 0;
1921 409 : int height = m_headerInfo.nRows;
1922 409 : int width = m_headerInfo.nCols;
1923 409 : int nDim = m_headerInfo.nDim;
1924 :
1925 409 : const unsigned int* arr = (const unsigned int*)(*ppByte);
1926 409 : const unsigned int* srcPtr = arr;
1927 409 : int bitPos = 0;
1928 409 : size_t nBytesRemaining = nBytesRemainingInOut;
1929 :
1930 409 : if (m_headerInfo.numValidPixel == width * height) // all valid
1931 : {
1932 409 : if (m_imageEncodeMode == IEM_DeltaHuffman)
1933 : {
1934 519 : for (int iDim = 0; iDim < nDim; iDim++)
1935 : {
1936 304 : T prevVal = 0;
1937 22880 : for (int m = iDim, i = 0; i < height; i++)
1938 3082030 : for (int j = 0; j < width; j++, m += nDim)
1939 : {
1940 3059450 : int val = 0;
1941 3059450 : if (nBytesRemaining >= 4 * sizeof(unsigned int))
1942 : {
1943 3051130 : if (!huffman.DecodeOneValue_NoOverrunCheck(&srcPtr, nBytesRemaining, bitPos, numBitsLUT, val))
1944 947 : return false;
1945 : }
1946 : else
1947 : {
1948 8324 : if (!huffman.DecodeOneValue(&srcPtr, nBytesRemaining, bitPos, numBitsLUT, val))
1949 947 : return false;
1950 : }
1951 :
1952 3059450 : T delta = (T)(val - offset);
1953 :
1954 3059450 : if (j > 0)
1955 3037860 : delta += prevVal; // use overflow
1956 21590 : else if (i > 0)
1957 22272 : delta += data[m - width * nDim];
1958 : else
1959 0 : delta += prevVal;
1960 :
1961 3059450 : data[m] = delta;
1962 3059450 : prevVal = delta;
1963 : }
1964 : }
1965 : }
1966 :
1967 194 : else if (m_imageEncodeMode == IEM_Huffman)
1968 : {
1969 7830 : for (int k = 0, m0 = 0, i = 0; i < height; i++)
1970 1618880 : for (int j = 0; j < width; j++, k++, m0 += nDim)
1971 3222500 : for (int m = 0; m < nDim; m++)
1972 : {
1973 1611250 : int val = 0;
1974 1611250 : if (nBytesRemaining >= 4 * sizeof(unsigned int))
1975 : {
1976 1609205 : if (!huffman.DecodeOneValue_NoOverrunCheck(&srcPtr, nBytesRemaining, bitPos, numBitsLUT, val))
1977 0 : return false;
1978 : }
1979 : else
1980 : {
1981 2042 : if (!huffman.DecodeOneValue(&srcPtr, nBytesRemaining, bitPos, numBitsLUT, val))
1982 0 : return false;
1983 : }
1984 :
1985 1611250 : data[m0 + m] = (T)(val - offset);
1986 : }
1987 : }
1988 :
1989 : else
1990 0 : return false;
1991 : }
1992 :
1993 : else // not all valid
1994 : {
1995 0 : if (m_imageEncodeMode == IEM_DeltaHuffman)
1996 : {
1997 0 : for (int iDim = 0; iDim < nDim; iDim++)
1998 : {
1999 0 : T prevVal = 0;
2000 0 : for (int k = 0, m = iDim, i = 0; i < height; i++)
2001 0 : for (int j = 0; j < width; j++, k++, m += nDim)
2002 0 : if (m_bitMask.IsValid(k))
2003 : {
2004 0 : int val = 0;
2005 0 : if (nBytesRemaining >= 4 * sizeof(unsigned int))
2006 : {
2007 0 : if (!huffman.DecodeOneValue_NoOverrunCheck(&srcPtr, nBytesRemaining, bitPos, numBitsLUT, val))
2008 0 : return false;
2009 : }
2010 : else
2011 : {
2012 0 : if (!huffman.DecodeOneValue(&srcPtr, nBytesRemaining, bitPos, numBitsLUT, val))
2013 0 : return false;
2014 : }
2015 :
2016 0 : T delta = (T)(val - offset);
2017 :
2018 0 : if (j > 0 && m_bitMask.IsValid(k - 1))
2019 : {
2020 0 : delta += prevVal; // use overflow
2021 : }
2022 0 : else if (i > 0 && m_bitMask.IsValid(k - width))
2023 : {
2024 0 : delta += data[m - width * nDim];
2025 : }
2026 : else
2027 0 : delta += prevVal;
2028 :
2029 0 : data[m] = delta;
2030 0 : prevVal = delta;
2031 : }
2032 : }
2033 : }
2034 :
2035 0 : else if (m_imageEncodeMode == IEM_Huffman)
2036 : {
2037 0 : for (int k = 0, m0 = 0, i = 0; i < height; i++)
2038 0 : for (int j = 0; j < width; j++, k++, m0 += nDim)
2039 0 : if (m_bitMask.IsValid(k))
2040 0 : for (int m = 0; m < nDim; m++)
2041 : {
2042 0 : int val = 0;
2043 0 : if (nBytesRemaining >= 4 * sizeof(unsigned int))
2044 : {
2045 0 : if (!huffman.DecodeOneValue_NoOverrunCheck(&srcPtr, nBytesRemaining, bitPos, numBitsLUT, val))
2046 0 : return false;
2047 : }
2048 : else
2049 : {
2050 0 : if (!huffman.DecodeOneValue(&srcPtr, nBytesRemaining, bitPos, numBitsLUT, val))
2051 0 : return false;
2052 : }
2053 :
2054 0 : data[m0 + m] = (T)(val - offset);
2055 : }
2056 : }
2057 :
2058 : else
2059 0 : return false;
2060 : }
2061 :
2062 409 : size_t numUInts = srcPtr - arr + (bitPos > 0 ? 1 : 0) + 1; // add one more as the decode LUT can read ahead
2063 409 : size_t len = numUInts * sizeof(unsigned int);
2064 :
2065 409 : if (nBytesRemainingInOut < len)
2066 0 : return false;
2067 :
2068 409 : *ppByte += len;
2069 409 : nBytesRemainingInOut -= len;
2070 409 : return true;
2071 : }
2072 :
2073 : // -------------------------------------------------------------------------- ;
2074 :
2075 : template<class T>
2076 341 : bool Lerc2::WriteMinMaxRanges(const T* /*data*/, Byte** ppByte) const
2077 : {
2078 341 : if (!ppByte || !(*ppByte))
2079 0 : return false;
2080 :
2081 : //printf("write min / max = %f %f\n", m_zMinVec[0], m_zMaxVec[0]);
2082 :
2083 341 : int nDim = m_headerInfo.nDim;
2084 341 : if (/* nDim < 2 || */ (int)m_zMinVec.size() != nDim || (int)m_zMaxVec.size() != nDim)
2085 0 : return false;
2086 :
2087 341 : std::vector<T> zVec(nDim);
2088 342 : size_t len = nDim * sizeof(T);
2089 :
2090 752 : for (int i = 0; i < nDim; i++)
2091 410 : zVec[i] = (T)m_zMinVec[i];
2092 :
2093 342 : memcpy(*ppByte, &zVec[0], len);
2094 341 : (*ppByte) += len;
2095 :
2096 752 : for (int i = 0; i < nDim; i++)
2097 410 : zVec[i] = (T)m_zMaxVec[i];
2098 :
2099 342 : memcpy(*ppByte, &zVec[0], len);
2100 341 : (*ppByte) += len;
2101 :
2102 341 : return true;
2103 : }
2104 :
2105 : // -------------------------------------------------------------------------- ;
2106 :
2107 : template<class T>
2108 505 : bool Lerc2::ReadMinMaxRanges(const Byte** ppByte, size_t& nBytesRemaining, const T* /*data*/)
2109 : {
2110 505 : if (!ppByte || !(*ppByte))
2111 0 : return false;
2112 :
2113 505 : int nDim = m_headerInfo.nDim;
2114 :
2115 505 : m_zMinVec.resize(nDim);
2116 505 : m_zMaxVec.resize(nDim);
2117 :
2118 1010 : std::vector<T> zVec(nDim);
2119 505 : size_t len = nDim * sizeof(T);
2120 :
2121 505 : if (nBytesRemaining < len || !memcpy(&zVec[0], *ppByte, len))
2122 0 : return false;
2123 :
2124 505 : (*ppByte) += len;
2125 505 : nBytesRemaining -= len;
2126 :
2127 1128 : for (int i = 0; i < nDim; i++)
2128 623 : m_zMinVec[i] = zVec[i];
2129 :
2130 505 : if (nBytesRemaining < len || !memcpy(&zVec[0], *ppByte, len))
2131 0 : return false;
2132 :
2133 505 : (*ppByte) += len;
2134 505 : nBytesRemaining -= len;
2135 :
2136 1128 : for (int i = 0; i < nDim; i++)
2137 623 : m_zMaxVec[i] = zVec[i];
2138 :
2139 : //printf("read min / max = %f %f\n", m_zMinVec[0], m_zMaxVec[0]);
2140 :
2141 505 : return true;
2142 : }
2143 :
2144 : // -------------------------------------------------------------------------- ;
2145 :
2146 : inline
2147 1189 : bool Lerc2::CheckMinMaxRanges(bool& minMaxEqual) const
2148 : {
2149 1189 : int nDim = m_headerInfo.nDim;
2150 1189 : if ((int)m_zMinVec.size() != nDim || (int)m_zMaxVec.size() != nDim)
2151 0 : return false;
2152 :
2153 1188 : minMaxEqual = (0 == memcmp(&m_zMinVec[0], &m_zMaxVec[0], nDim * sizeof(m_zMinVec[0])));
2154 1188 : return true;
2155 : }
2156 :
2157 : // -------------------------------------------------------------------------- ;
2158 :
2159 : template<class T>
2160 576 : bool Lerc2::FillConstImage(T* data) const
2161 : {
2162 576 : if (!data)
2163 0 : return false;
2164 :
2165 576 : const HeaderInfo& hd = m_headerInfo;
2166 576 : int nCols = hd.nCols;
2167 576 : int nRows = hd.nRows;
2168 576 : int nDim = hd.nDim;
2169 576 : T z0 = (T)hd.zMin;
2170 :
2171 576 : if (nDim == 1)
2172 : {
2173 72143 : for (int k = 0, i = 0; i < nRows; i++)
2174 9231438 : for (int j = 0; j < nCols; j++, k++)
2175 9159872 : if (m_bitMask.IsValid(k))
2176 9159468 : data[k] = z0;
2177 : }
2178 : else
2179 : {
2180 1 : std::vector<T> zBufVec(nDim, z0);
2181 :
2182 1 : if (hd.zMin != hd.zMax)
2183 : {
2184 0 : if ((int)m_zMinVec.size() != nDim)
2185 0 : return false;
2186 :
2187 0 : for (int m = 0; m < nDim; m++)
2188 0 : zBufVec[m] = (T)m_zMinVec[m];
2189 : }
2190 :
2191 1 : int len = nDim * sizeof(T);
2192 3 : for (int k = 0, m = 0, i = 0; i < nRows; i++)
2193 6 : for (int j = 0; j < nCols; j++, k++, m += nDim)
2194 4 : if (m_bitMask.IsValid(k))
2195 3 : memcpy(&data[m], &zBufVec[0], len);
2196 : }
2197 :
2198 576 : return true;
2199 : }
2200 :
2201 : // -------------------------------------------------------------------------- ;
2202 :
2203 : NAMESPACE_LERC_END
2204 : #endif
|
