1 /*
2 Copyright 2015 Esri
3 
4 Licensed under the Apache License, Version 2.0 (the "License");
5 you may not use this file except in compliance with the License.
6 You may obtain a copy of the License at
7 
8 http://www.apache.org/licenses/LICENSE-2.0
9 
10 Unless required by applicable law or agreed to in writing, software
11 distributed under the License is distributed on an "AS IS" BASIS,
12 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 See the License for the specific language governing permissions and
14 limitations under the License.
15 
16 A local copy of the license and additional notices are located with the
17 source distribution at:
18 
19 http://github.com/Esri/lerc/
20 
21 Contributors:  Thomas Maurer
22 */
23 
24 #include "Defines.h"
25 #include "Lerc.h"
26 #include "Lerc2.h"
27 #include <typeinfo>
28 #include <limits>
29 
30 #ifdef HAVE_LERC1_DECODE
31 #include "Lerc1Decode/CntZImage.h"
32 #endif
33 
34 using namespace std;
35 USING_NAMESPACE_LERC
36 
37 // -------------------------------------------------------------------------- ;
38 
ComputeCompressedSize(const void * pData,int version,DataType dt,int nDim,int nCols,int nRows,int nBands,const BitMask * pBitMask,double maxZErr,unsigned int & numBytesNeeded)39 ErrCode Lerc::ComputeCompressedSize(const void* pData, int version, DataType dt, int nDim, int nCols, int nRows, int nBands,
40   const BitMask* pBitMask, double maxZErr, unsigned int& numBytesNeeded)
41 {
42   switch (dt)
43   {
44   case DT_Char:    return ComputeCompressedSizeTempl((const signed char*)pData, version, nDim, nCols, nRows, nBands, pBitMask, maxZErr, numBytesNeeded);
45   case DT_Byte:    return ComputeCompressedSizeTempl((const Byte*)pData, version, nDim, nCols, nRows, nBands, pBitMask, maxZErr, numBytesNeeded);
46   case DT_Short:   return ComputeCompressedSizeTempl((const short*)pData, version, nDim, nCols, nRows, nBands, pBitMask, maxZErr, numBytesNeeded);
47   case DT_UShort:  return ComputeCompressedSizeTempl((const unsigned short*)pData, version, nDim, nCols, nRows, nBands, pBitMask, maxZErr, numBytesNeeded);
48   case DT_Int:     return ComputeCompressedSizeTempl((const int*)pData, version, nDim, nCols, nRows, nBands, pBitMask, maxZErr, numBytesNeeded);
49   case DT_UInt:    return ComputeCompressedSizeTempl((const unsigned int*)pData, version, nDim, nCols, nRows, nBands, pBitMask, maxZErr, numBytesNeeded);
50   case DT_Float:   return ComputeCompressedSizeTempl((const float*)pData, version, nDim, nCols, nRows, nBands, pBitMask, maxZErr, numBytesNeeded);
51   case DT_Double:  return ComputeCompressedSizeTempl((const double*)pData, version, nDim, nCols, nRows, nBands, pBitMask, maxZErr, numBytesNeeded);
52 
53   default:
54     return ErrCode::WrongParam;
55   }
56 }
57 
58 // -------------------------------------------------------------------------- ;
59 
Encode(const void * pData,int version,DataType dt,int nDim,int nCols,int nRows,int nBands,const BitMask * pBitMask,double maxZErr,Byte * pBuffer,unsigned int numBytesBuffer,unsigned int & numBytesWritten)60 ErrCode Lerc::Encode(const void* pData, int version, DataType dt, int nDim, int nCols, int nRows, int nBands,
61   const BitMask* pBitMask, double maxZErr, Byte* pBuffer, unsigned int numBytesBuffer, unsigned int& numBytesWritten)
62 {
63   switch (dt)
64   {
65   case DT_Char:    return EncodeTempl((const signed char*)pData, version, nDim, nCols, nRows, nBands, pBitMask, maxZErr, pBuffer, numBytesBuffer, numBytesWritten);
66   case DT_Byte:    return EncodeTempl((const Byte*)pData, version, nDim, nCols, nRows, nBands, pBitMask, maxZErr, pBuffer, numBytesBuffer, numBytesWritten);
67   case DT_Short:   return EncodeTempl((const short*)pData, version, nDim, nCols, nRows, nBands, pBitMask, maxZErr, pBuffer, numBytesBuffer, numBytesWritten);
68   case DT_UShort:  return EncodeTempl((const unsigned short*)pData, version, nDim, nCols, nRows, nBands, pBitMask, maxZErr, pBuffer, numBytesBuffer, numBytesWritten);
69   case DT_Int:     return EncodeTempl((const int*)pData, version, nDim, nCols, nRows, nBands, pBitMask, maxZErr, pBuffer, numBytesBuffer, numBytesWritten);
70   case DT_UInt:    return EncodeTempl((const unsigned int*)pData, version, nDim, nCols, nRows, nBands, pBitMask, maxZErr, pBuffer, numBytesBuffer, numBytesWritten);
71   case DT_Float:   return EncodeTempl((const float*)pData, version, nDim, nCols, nRows, nBands, pBitMask, maxZErr, pBuffer, numBytesBuffer, numBytesWritten);
72   case DT_Double:  return EncodeTempl((const double*)pData, version, nDim, nCols, nRows, nBands, pBitMask, maxZErr, pBuffer, numBytesBuffer, numBytesWritten);
73 
74   default:
75     return ErrCode::WrongParam;
76   }
77 }
78 
79 // -------------------------------------------------------------------------- ;
80 
GetLercInfo(const Byte * pLercBlob,unsigned int numBytesBlob,struct LercInfo & lercInfo)81 ErrCode Lerc::GetLercInfo(const Byte* pLercBlob, unsigned int numBytesBlob, struct LercInfo& lercInfo)
82 {
83   lercInfo.RawInit();
84 
85   // first try Lerc2
86   struct Lerc2::HeaderInfo lerc2Info;
87   if (Lerc2::GetHeaderInfo(pLercBlob, numBytesBlob, lerc2Info))
88   {
89     lercInfo.version = lerc2Info.version;
90     lercInfo.nDim = lerc2Info.nDim;
91     lercInfo.nCols = lerc2Info.nCols;
92     lercInfo.nRows = lerc2Info.nRows;
93     lercInfo.numValidPixel = lerc2Info.numValidPixel;
94     lercInfo.nBands = 1;
95     lercInfo.blobSize = lerc2Info.blobSize;
96     lercInfo.dt = (DataType)lerc2Info.dt;
97     lercInfo.zMin = lerc2Info.zMin;
98     lercInfo.zMax = lerc2Info.zMax;
99     lercInfo.maxZError = lerc2Info.maxZError;
100 
101     if (lercInfo.blobSize > (int)numBytesBlob)    // truncated blob, we won't be able to read this band
102       return ErrCode::BufferTooSmall;
103 
104     struct Lerc2::HeaderInfo hdInfo;
105     while (Lerc2::GetHeaderInfo(pLercBlob + lercInfo.blobSize, numBytesBlob - lercInfo.blobSize, hdInfo))
106     {
107       if (hdInfo.nDim != lercInfo.nDim
108        || hdInfo.nCols != lercInfo.nCols
109        || hdInfo.nRows != lercInfo.nRows
110        || hdInfo.numValidPixel != lercInfo.numValidPixel
111        || (int)hdInfo.dt != (int)lercInfo.dt)
112        //|| hdInfo.maxZError != lercInfo.maxZError)  // with the new bitplane compression, maxZError can vary between bands
113       {
114         return ErrCode::Failed;
115       }
116 
117       if (lercInfo.blobSize > std::numeric_limits<int>::max() - hdInfo.blobSize)
118         return ErrCode::Failed;
119 
120       lercInfo.blobSize += hdInfo.blobSize;
121 
122       if (lercInfo.blobSize > (int)numBytesBlob)    // truncated blob, we won't be able to read this band
123         return ErrCode::BufferTooSmall;
124 
125       lercInfo.nBands++;
126       lercInfo.zMin = min(lercInfo.zMin, hdInfo.zMin);
127       lercInfo.zMax = max(lercInfo.zMax, hdInfo.zMax);
128       lercInfo.maxZError = max(lercInfo.maxZError, hdInfo.maxZError);  // with the new bitplane compression, maxZError can vary between bands
129     }
130 
131     return ErrCode::Ok;
132   }
133 
134 
135 #ifdef HAVE_LERC1_DECODE
136   // only if not Lerc2, try legacy Lerc1
137   unsigned int numBytesHeaderBand0 = CntZImage::computeNumBytesNeededToReadHeader(false);
138   unsigned int numBytesHeaderBand1 = CntZImage::computeNumBytesNeededToReadHeader(true);
139   Byte* pByte = const_cast<Byte*>(pLercBlob);
140 
141   lercInfo.zMin =  FLT_MAX;
142   lercInfo.zMax = -FLT_MAX;
143 
144   CntZImage cntZImg;
145   if (numBytesHeaderBand0 <= numBytesBlob && cntZImg.read(&pByte, 1e12, true))    // read just the header
146   {
147     size_t nBytesRead = pByte - pLercBlob;
148     size_t nBytesNeeded = 10 + 4 * sizeof(int) + 1 * sizeof(double);
149 
150     if (nBytesRead < nBytesNeeded)
151       return ErrCode::Failed;
152 
153     Byte* ptr = const_cast<Byte*>(pLercBlob);
154     ptr += 10 + 2 * sizeof(int);
155 
156     int height(0), width(0);
157     memcpy(&height, ptr, sizeof(int));  ptr += sizeof(int);
158     memcpy(&width,  ptr, sizeof(int));  ptr += sizeof(int);
159     double maxZErrorInFile(0);
160     memcpy(&maxZErrorInFile, ptr, sizeof(double));
161 
162     if (height > 20000 || width > 20000)    // guard against bogus numbers; size limitation for old Lerc1
163       return ErrCode::Failed;
164 
165     lercInfo.nDim = 1;
166     lercInfo.nCols = width;
167     lercInfo.nRows = height;
168     lercInfo.dt = Lerc::DT_Float;
169     lercInfo.maxZError = maxZErrorInFile;
170 
171     Byte* pByte = const_cast<Byte*>(pLercBlob);
172     bool onlyZPart = false;
173 
174     while (lercInfo.blobSize + numBytesHeaderBand1 < numBytesBlob)    // means there could be another band
175     {
176       if (!cntZImg.read(&pByte, 1e12, false, onlyZPart))
177         return (lercInfo.nBands > 0) ? ErrCode::Ok : ErrCode::Failed;    // no other band, we are done
178 
179       onlyZPart = true;
180 
181       lercInfo.nBands++;
182       lercInfo.blobSize = (int)(pByte - pLercBlob);
183 
184       // now that we have decoded it, we can go the extra mile and collect some extra info
185       int numValidPixels = 0;
186       float zMin =  FLT_MAX;
187       float zMax = -FLT_MAX;
188 
189       for (int i = 0; i < height; i++)
190       {
191         for (int j = 0; j < width; j++)
192           if (cntZImg(i, j).cnt > 0)
193           {
194             numValidPixels++;
195             float z = cntZImg(i, j).z;
196             zMax = max(zMax, z);
197             zMin = min(zMin, z);
198           }
199       }
200 
201       lercInfo.numValidPixel = numValidPixels;
202       lercInfo.zMin = std::min(lercInfo.zMin, (double)zMin);
203       lercInfo.zMax = std::max(lercInfo.zMax, (double)zMax);
204     }
205 
206     return ErrCode::Ok;
207   }
208 #endif
209 
210   return ErrCode::Failed;
211 }
212 
213 // -------------------------------------------------------------------------- ;
214 
Decode(const Byte * pLercBlob,unsigned int numBytesBlob,BitMask * pBitMask,int nDim,int nCols,int nRows,int nBands,DataType dt,void * pData)215 ErrCode Lerc::Decode(const Byte* pLercBlob, unsigned int numBytesBlob, BitMask* pBitMask,
216   int nDim, int nCols, int nRows, int nBands, DataType dt, void* pData)
217 {
218   switch (dt)
219   {
220   case DT_Char:    return DecodeTempl((signed char*)pData, pLercBlob, numBytesBlob, nDim, nCols, nRows, nBands, pBitMask);
221   case DT_Byte:    return DecodeTempl((Byte*)pData, pLercBlob, numBytesBlob, nDim, nCols, nRows, nBands, pBitMask);
222   case DT_Short:   return DecodeTempl((short*)pData, pLercBlob, numBytesBlob, nDim, nCols, nRows, nBands, pBitMask);
223   case DT_UShort:  return DecodeTempl((unsigned short*)pData, pLercBlob, numBytesBlob, nDim, nCols, nRows, nBands, pBitMask);
224   case DT_Int:     return DecodeTempl((int*)pData, pLercBlob, numBytesBlob, nDim, nCols, nRows, nBands, pBitMask);
225   case DT_UInt:    return DecodeTempl((unsigned int*)pData, pLercBlob, numBytesBlob, nDim, nCols, nRows, nBands, pBitMask);
226   case DT_Float:   return DecodeTempl((float*)pData, pLercBlob, numBytesBlob, nDim, nCols, nRows, nBands, pBitMask);
227   case DT_Double:  return DecodeTempl((double*)pData, pLercBlob, numBytesBlob, nDim, nCols, nRows, nBands, pBitMask);
228 
229   default:
230     return ErrCode::WrongParam;
231   }
232 }
233 
234 // -------------------------------------------------------------------------- ;
235 
ConvertToDouble(const void * pDataIn,DataType dt,size_t nDataValues,double * pDataOut)236 ErrCode Lerc::ConvertToDouble(const void* pDataIn, DataType dt, size_t nDataValues, double* pDataOut)
237 {
238   switch (dt)
239   {
240   case DT_Char:    return ConvertToDoubleTempl((const signed char*)pDataIn, nDataValues, pDataOut);
241   case DT_Byte:    return ConvertToDoubleTempl((const Byte*)pDataIn, nDataValues, pDataOut);
242   case DT_Short:   return ConvertToDoubleTempl((const short*)pDataIn, nDataValues, pDataOut);
243   case DT_UShort:  return ConvertToDoubleTempl((const unsigned short*)pDataIn, nDataValues, pDataOut);
244   case DT_Int:     return ConvertToDoubleTempl((const int*)pDataIn, nDataValues, pDataOut);
245   case DT_UInt:    return ConvertToDoubleTempl((const unsigned int*)pDataIn, nDataValues, pDataOut);
246   case DT_Float:   return ConvertToDoubleTempl((const float*)pDataIn, nDataValues, pDataOut);
247   //case DT_Double:  no convert double to double
248 
249   default:
250     return ErrCode::WrongParam;
251   }
252 }
253 
254 // -------------------------------------------------------------------------- ;
255 // -------------------------------------------------------------------------- ;
256 
257 template<class T>
ComputeCompressedSizeTempl(const T * pData,int version,int nDim,int nCols,int nRows,int nBands,const BitMask * pBitMask,double maxZErr,unsigned int & numBytesNeeded)258 ErrCode Lerc::ComputeCompressedSizeTempl(const T* pData, int version, int nDim, int nCols, int nRows, int nBands,
259   const BitMask* pBitMask, double maxZErr, unsigned int& numBytesNeeded)
260 {
261   numBytesNeeded = 0;
262 
263   if (!pData || nDim <= 0 || nCols <= 0 || nRows <= 0 || nBands <= 0 || maxZErr < 0)
264     return ErrCode::WrongParam;
265 
266   if (pBitMask && (pBitMask->GetHeight() != nRows || pBitMask->GetWidth() != nCols))
267     return ErrCode::WrongParam;
268 
269   Lerc2 lerc2;
270   if (version >= 0 && !lerc2.SetEncoderToOldVersion(version))
271     return ErrCode::WrongParam;
272 
273   bool rv = pBitMask ? lerc2.Set(nDim, nCols, nRows, pBitMask->Bits()) : lerc2.Set(nDim, nCols, nRows);
274   if (!rv)
275     return ErrCode::Failed;
276 
277   // loop over the bands
278   for (int iBand = 0; iBand < nBands; iBand++)
279   {
280     bool encMsk = (iBand == 0);    // store bit mask with first band only
281     const T* arr = pData + nDim * nCols * nRows * iBand;
282 
283     ErrCode errCode = CheckForNaN(arr, nDim, nCols, nRows, pBitMask);
284     if (errCode != ErrCode::Ok)
285       return errCode;
286 
287     unsigned int nBytes = lerc2.ComputeNumBytesNeededToWrite(arr, maxZErr, encMsk);
288     if (nBytes <= 0)
289       return ErrCode::Failed;
290 
291     numBytesNeeded += nBytes;
292   }
293 
294   return ErrCode::Ok;
295 }
296 
297 // -------------------------------------------------------------------------- ;
298 
299 template<class T>
EncodeTempl(const T * pData,int version,int nDim,int nCols,int nRows,int nBands,const BitMask * pBitMask,double maxZErr,Byte * pBuffer,unsigned int numBytesBuffer,unsigned int & numBytesWritten)300 ErrCode Lerc::EncodeTempl(const T* pData, int version, int nDim, int nCols, int nRows, int nBands,
301   const BitMask* pBitMask, double maxZErr, Byte* pBuffer, unsigned int numBytesBuffer, unsigned int& numBytesWritten)
302 {
303   numBytesWritten = 0;
304 
305   if (!pData || nDim <= 0 || nCols <= 0 || nRows <= 0 || nBands <= 0 || maxZErr < 0 || !pBuffer || !numBytesBuffer)
306     return ErrCode::WrongParam;
307 
308   if (pBitMask && (pBitMask->GetHeight() != nRows || pBitMask->GetWidth() != nCols))
309     return ErrCode::WrongParam;
310 
311   Lerc2 lerc2;
312   if (version >= 0 && !lerc2.SetEncoderToOldVersion(version))
313     return ErrCode::WrongParam;
314 
315   bool rv = pBitMask ? lerc2.Set(nDim, nCols, nRows, pBitMask->Bits()) : lerc2.Set(nDim, nCols, nRows);
316   if (!rv)
317     return ErrCode::Failed;
318 
319   Byte* pByte = pBuffer;
320 
321   // loop over the bands, encode into array of single band Lerc blobs
322   for (int iBand = 0; iBand < nBands; iBand++)
323   {
324     bool encMsk = (iBand == 0);    // store bit mask with first band only
325     const T* arr = pData + nDim * nCols * nRows * iBand;
326 
327     ErrCode errCode = CheckForNaN(arr, nDim, nCols, nRows, pBitMask);
328     if (errCode != ErrCode::Ok)
329       return errCode;
330 
331     unsigned int nBytes = lerc2.ComputeNumBytesNeededToWrite(arr, maxZErr, encMsk);
332     if (nBytes == 0)
333       return ErrCode::Failed;
334 
335     unsigned int nBytesAlloc = nBytes;
336 
337     if ((size_t)(pByte - pBuffer) + nBytesAlloc > numBytesBuffer)    // check we have enough space left
338       return ErrCode::BufferTooSmall;
339 
340     if (!lerc2.Encode(arr, &pByte))
341       return ErrCode::Failed;
342   }
343 
344   numBytesWritten = (unsigned int)(pByte - pBuffer);
345   return ErrCode::Ok;
346 }
347 
348 // -------------------------------------------------------------------------- ;
349 
350 template<class T>
DecodeTempl(T * pData,const Byte * pLercBlob,unsigned int numBytesBlob,int nDim,int nCols,int nRows,int nBands,BitMask * pBitMask)351 ErrCode Lerc::DecodeTempl(T* pData, const Byte* pLercBlob, unsigned int numBytesBlob,
352   int nDim, int nCols, int nRows, int nBands, BitMask* pBitMask)
353 {
354   if (!pData || nDim <= 0 || nCols <= 0 || nRows <= 0 || nBands <= 0 || !pLercBlob || !numBytesBlob)
355     return ErrCode::WrongParam;
356 
357   if (pBitMask && (pBitMask->GetHeight() != nRows || pBitMask->GetWidth() != nCols))
358     return ErrCode::WrongParam;
359 
360   const Byte* pByte = pLercBlob;
361 #ifdef HAVE_LERC1_DECODE
362   Byte* pByte1 = const_cast<Byte*>(pLercBlob);
363 #endif
364   Lerc2::HeaderInfo hdInfo;
365 
366   if (Lerc2::GetHeaderInfo(pByte, numBytesBlob, hdInfo) && hdInfo.version >= 1)    // is Lerc2
367   {
368     size_t nBytesRemaining = numBytesBlob;
369     Lerc2 lerc2;
370 
371     for (int iBand = 0; iBand < nBands; iBand++)
372     {
373       if (((size_t)(pByte - pLercBlob) < numBytesBlob) && Lerc2::GetHeaderInfo(pByte, nBytesRemaining, hdInfo))
374       {
375         if (hdInfo.nDim != nDim || hdInfo.nCols != nCols || hdInfo.nRows != nRows)
376           return ErrCode::Failed;
377 
378         if ((pByte - pLercBlob) + (size_t)hdInfo.blobSize > numBytesBlob)
379           return ErrCode::BufferTooSmall;
380 
381         T* arr = pData + nDim * nCols * nRows * iBand;
382 
383         if (!lerc2.Decode(&pByte, nBytesRemaining, arr, (pBitMask && iBand == 0) ? pBitMask->Bits() : nullptr))
384           return ErrCode::Failed;
385       }
386     }
387   }
388 
389   else    // might be old Lerc1
390   {
391 #ifdef HAVE_LERC1_DECODE
392     unsigned int numBytesHeaderBand0 = CntZImage::computeNumBytesNeededToReadHeader(false);
393     unsigned int numBytesHeaderBand1 = CntZImage::computeNumBytesNeededToReadHeader(true);
394     CntZImage zImg;
395 
396     for (int iBand = 0; iBand < nBands; iBand++)
397     {
398       unsigned int numBytesHeader = iBand == 0 ? numBytesHeaderBand0 : numBytesHeaderBand1;
399       if ((size_t)(pByte - pLercBlob) + numBytesHeader > numBytesBlob)
400         return ErrCode::BufferTooSmall;
401 
402       bool onlyZPart = iBand > 0;
403       if (!zImg.read(&pByte1, 1e12, false, onlyZPart))
404         return ErrCode::Failed;
405 
406       if (zImg.getWidth() != nCols || zImg.getHeight() != nRows)
407         return ErrCode::Failed;
408 
409       T* arr = pData + nCols * nRows * iBand;
410 
411       if (!Convert(zImg, arr, pBitMask))
412         return ErrCode::Failed;
413     }
414 #else
415     return ErrCode::Failed;
416 #endif
417   }
418 
419   return ErrCode::Ok;
420 }
421 
422 // -------------------------------------------------------------------------- ;
423 // -------------------------------------------------------------------------- ;
424 
425 #ifdef HAVE_LERC1_DECODE
426 template<class T>
Convert(const CntZImage & zImg,T * arr,BitMask * pBitMask)427 bool Lerc::Convert(const CntZImage& zImg, T* arr, BitMask* pBitMask)
428 {
429   if (!arr || !zImg.getSize())
430     return false;
431 
432   const bool fltPnt = (typeid(*arr) == typeid(double)) || (typeid(*arr) == typeid(float));
433 
434   int h = zImg.getHeight();
435   int w = zImg.getWidth();
436 
437   if (pBitMask && (pBitMask->GetHeight() != h || pBitMask->GetWidth() != w))
438     return false;
439 
440   if (pBitMask)
441     pBitMask->SetAllValid();
442 
443   const CntZ* srcPtr = zImg.getData();
444   T* dstPtr = arr;
445   int num = w * h;
446   for (int k = 0; k < num; k++)
447   {
448     if (srcPtr->cnt > 0)
449       *dstPtr = fltPnt ? (T)srcPtr->z : (T)floor(srcPtr->z + 0.5);
450     else if (pBitMask)
451       pBitMask->SetInvalid(k);
452 
453     srcPtr++;
454     dstPtr++;
455   }
456 
457   return true;
458 }
459 #endif
460 
461 // -------------------------------------------------------------------------- ;
462 
463 template<class T>
ConvertToDoubleTempl(const T * pDataIn,size_t nDataValues,double * pDataOut)464 ErrCode Lerc::ConvertToDoubleTempl(const T* pDataIn, size_t nDataValues, double* pDataOut)
465 {
466   if (!pDataIn || !nDataValues || !pDataOut)
467     return ErrCode::WrongParam;
468 
469   for (size_t k = 0; k < nDataValues; k++)
470     pDataOut[k] = pDataIn[k];
471 
472   return ErrCode::Ok;
473 }
474 
475 // -------------------------------------------------------------------------- ;
476 
CheckForNaN(const T * arr,int nDim,int nCols,int nRows,const BitMask * pBitMask)477 template<class T> ErrCode Lerc::CheckForNaN(const T* arr, int nDim, int nCols, int nRows, const BitMask* pBitMask)
478 {
479   if (!arr || nDim <= 0 || nCols <= 0 || nRows <= 0)
480     return ErrCode::WrongParam;
481 
482 #ifdef CHECK_FOR_NAN
483 
484   if (typeid(T) == typeid(double) || typeid(T) == typeid(float))
485   {
486     bool foundNaN = false;
487 
488     for (int k = 0, i = 0; i < nRows; i++)
489     {
490       const T* rowArr = &(arr[i * nCols * nDim]);
491 
492       if (!pBitMask)    // all valid
493       {
494         for (int n = 0, j = 0; j < nCols; j++, n += nDim)
495           for (int m = 0; m < nDim; m++)
496             if (std::isnan((double)rowArr[n + m]))
497               foundNaN = true;
498       }
499       else    // not all valid
500       {
501         for (int n = 0, j = 0; j < nCols; j++, k++, n += nDim)
502           if (pBitMask->IsValid(k))
503           {
504             for (int m = 0; m < nDim; m++)
505               if (std::isnan((double)rowArr[n + m]))
506                 foundNaN = true;
507           }
508       }
509 
510       if (foundNaN)
511         return ErrCode::NaN;
512     }
513   }
514 
515 #endif
516 
517   return ErrCode::Ok;
518 }
519 
520 // -------------------------------------------------------------------------- ;
521 
522