Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Virtual GDAL Datasets
4 : * Purpose: Implementation of some filter types.
5 : * Author: Frank Warmerdam <warmerdam@pobox.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2003, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_port.h"
15 : #include "vrtdataset.h"
16 :
17 : #include <algorithm>
18 : #include <cmath>
19 : #include <cstddef>
20 : #include <cstdlib>
21 : #include <cstring>
22 : #include <limits>
23 : #include <map>
24 :
25 : #include "cpl_conv.h"
26 : #include "cpl_error.h"
27 : #include "cpl_minixml.h"
28 : #include "cpl_string.h"
29 : #include "cpl_vsi.h"
30 : #include "gdal.h"
31 : #include "gdal_priv.h"
32 :
33 : /*! @cond Doxygen_Suppress */
34 :
35 : /************************************************************************/
36 : /* ==================================================================== */
37 : /* VRTFilteredSource */
38 : /* ==================================================================== */
39 : /************************************************************************/
40 :
41 : /************************************************************************/
42 : /* VRTFilteredSource() */
43 : /************************************************************************/
44 :
45 42 : VRTFilteredSource::VRTFilteredSource()
46 42 : : m_nSupportedTypesCount(1), m_nExtraEdgePixels(0)
47 : {
48 882 : for (size_t i = 0; i < CPL_ARRAYSIZE(m_aeSupportedTypes); ++i)
49 840 : m_aeSupportedTypes[i] = GDT_Unknown;
50 :
51 42 : m_aeSupportedTypes[0] = GDT_Float32;
52 42 : }
53 :
54 : /************************************************************************/
55 : /* ~VRTFilteredSource() */
56 : /************************************************************************/
57 :
58 42 : VRTFilteredSource::~VRTFilteredSource()
59 : {
60 42 : }
61 :
62 : /************************************************************************/
63 : /* SetExtraEdgePixels() */
64 : /************************************************************************/
65 :
66 40 : void VRTFilteredSource::SetExtraEdgePixels(int nEdgePixels)
67 :
68 : {
69 40 : m_nExtraEdgePixels = nEdgePixels;
70 40 : }
71 :
72 : /************************************************************************/
73 : /* SetFilteringDataTypesSupported() */
74 : /************************************************************************/
75 :
76 42 : void VRTFilteredSource::SetFilteringDataTypesSupported(int nTypeCount,
77 : GDALDataType *paeTypes)
78 :
79 : {
80 42 : if (nTypeCount >
81 : static_cast<int>(sizeof(m_aeSupportedTypes) / sizeof(GDALDataType)))
82 : {
83 0 : CPLAssert(false);
84 : nTypeCount =
85 : static_cast<int>(sizeof(m_aeSupportedTypes) / sizeof(GDALDataType));
86 : }
87 :
88 42 : m_nSupportedTypesCount = nTypeCount;
89 42 : memcpy(m_aeSupportedTypes, paeTypes, sizeof(GDALDataType) * nTypeCount);
90 42 : }
91 :
92 : /************************************************************************/
93 : /* IsTypeSupported() */
94 : /************************************************************************/
95 :
96 88 : int VRTFilteredSource::IsTypeSupported(GDALDataType eTestType) const
97 :
98 : {
99 176 : for (int i = 0; i < m_nSupportedTypesCount; i++)
100 : {
101 88 : if (eTestType == m_aeSupportedTypes[i])
102 0 : return TRUE;
103 : }
104 :
105 88 : return FALSE;
106 : }
107 :
108 : /************************************************************************/
109 : /* RasterIO() */
110 : /************************************************************************/
111 :
112 44 : CPLErr VRTFilteredSource::RasterIO(GDALDataType eVRTBandDataType, int nXOff,
113 : int nYOff, int nXSize, int nYSize,
114 : void *pData, int nBufXSize, int nBufYSize,
115 : GDALDataType eBufType, GSpacing nPixelSpace,
116 : GSpacing nLineSpace,
117 : GDALRasterIOExtraArg *psExtraArg,
118 : WorkingState &oWorkingState)
119 :
120 : {
121 : /* -------------------------------------------------------------------- */
122 : /* For now we don't support filtered access to non-full */
123 : /* resolution requests. Just collect the data directly without */
124 : /* any operator. */
125 : /* -------------------------------------------------------------------- */
126 44 : if (nBufXSize != nXSize || nBufYSize != nYSize)
127 : {
128 0 : return VRTComplexSource::RasterIO(
129 : eVRTBandDataType, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize,
130 : nBufYSize, eBufType, nPixelSpace, nLineSpace, psExtraArg,
131 0 : oWorkingState);
132 : }
133 :
134 44 : double dfXOff = nXOff;
135 44 : double dfYOff = nYOff;
136 44 : double dfXSize = nXSize;
137 44 : double dfYSize = nYSize;
138 44 : if (psExtraArg != nullptr && psExtraArg->bFloatingPointWindowValidity)
139 : {
140 4 : dfXOff = psExtraArg->dfXOff;
141 4 : dfYOff = psExtraArg->dfYOff;
142 4 : dfXSize = psExtraArg->dfXSize;
143 4 : dfYSize = psExtraArg->dfYSize;
144 : }
145 :
146 : // The window we will actually request from the source raster band.
147 44 : double dfReqXOff = 0.0;
148 44 : double dfReqYOff = 0.0;
149 44 : double dfReqXSize = 0.0;
150 44 : double dfReqYSize = 0.0;
151 44 : int nReqXOff = 0;
152 44 : int nReqYOff = 0;
153 44 : int nReqXSize = 0;
154 44 : int nReqYSize = 0;
155 :
156 : // The window we will actual set _within_ the pData buffer.
157 44 : int nOutXOff = 0;
158 44 : int nOutYOff = 0;
159 44 : int nOutXSize = 0;
160 44 : int nOutYSize = 0;
161 :
162 44 : bool bError = false;
163 44 : if (!GetSrcDstWindow(dfXOff, dfYOff, dfXSize, dfYSize, nBufXSize, nBufYSize,
164 : &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize,
165 : &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
166 : &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize, bError))
167 : {
168 0 : return bError ? CE_Failure : CE_None;
169 : }
170 :
171 44 : pData = reinterpret_cast<GByte *>(pData) + nPixelSpace * nOutXOff +
172 44 : nLineSpace * nOutYOff;
173 :
174 : /* -------------------------------------------------------------------- */
175 : /* Determine the data type we want to request. We try to match */
176 : /* the source or destination request, and if both those fail we */
177 : /* fallback to the first supported type at least as expressive */
178 : /* as the request. */
179 : /* -------------------------------------------------------------------- */
180 44 : GDALDataType eOperDataType = GDT_Unknown;
181 :
182 44 : if (IsTypeSupported(eBufType))
183 0 : eOperDataType = eBufType;
184 :
185 44 : auto l_band = GetRasterBand();
186 44 : if (!l_band)
187 0 : return CE_Failure;
188 :
189 88 : if (eOperDataType == GDT_Unknown &&
190 44 : IsTypeSupported(l_band->GetRasterDataType()))
191 0 : eOperDataType = l_band->GetRasterDataType();
192 :
193 44 : if (eOperDataType == GDT_Unknown)
194 : {
195 88 : for (int i = 0; i < m_nSupportedTypesCount; i++)
196 : {
197 44 : if (GDALDataTypeUnion(m_aeSupportedTypes[i], eBufType) ==
198 44 : m_aeSupportedTypes[i])
199 : {
200 7 : eOperDataType = m_aeSupportedTypes[i];
201 : }
202 : }
203 : }
204 :
205 44 : if (eOperDataType == GDT_Unknown)
206 : {
207 37 : eOperDataType = m_aeSupportedTypes[0];
208 :
209 37 : for (int i = 1; i < m_nSupportedTypesCount; i++)
210 : {
211 0 : if (GDALGetDataTypeSizeBytes(m_aeSupportedTypes[i]) >
212 0 : GDALGetDataTypeSizeBytes(eOperDataType))
213 : {
214 0 : eOperDataType = m_aeSupportedTypes[i];
215 : }
216 : }
217 : }
218 :
219 : /* -------------------------------------------------------------------- */
220 : /* Allocate the buffer of data into which our imagery will be */
221 : /* read, with the extra edge pixels as well. This will be the */
222 : /* source data fed into the filter. */
223 : /* -------------------------------------------------------------------- */
224 44 : if (nOutXSize > INT_MAX - 2 * m_nExtraEdgePixels ||
225 44 : nOutYSize > INT_MAX - 2 * m_nExtraEdgePixels)
226 : {
227 0 : return CE_Failure;
228 : }
229 44 : const int nExtraXSize = nOutXSize + 2 * m_nExtraEdgePixels;
230 44 : const int nExtraYSize = nOutYSize + 2 * m_nExtraEdgePixels;
231 :
232 44 : GByte *pabyWorkData = static_cast<GByte *>(VSI_MALLOC3_VERBOSE(
233 : nExtraXSize, nExtraYSize, GDALGetDataTypeSizeBytes(eOperDataType)));
234 :
235 44 : if (pabyWorkData == nullptr)
236 : {
237 0 : return CE_Failure;
238 : }
239 :
240 44 : const GPtrDiff_t nPixelOffset = GDALGetDataTypeSizeBytes(eOperDataType);
241 44 : const GPtrDiff_t nLineOffset = nPixelOffset * nExtraXSize;
242 :
243 44 : int bHasNoData = false;
244 44 : const double dfSrcNoDataValue = l_band->GetNoDataValue(&bHasNoData);
245 44 : if (bHasNoData)
246 : {
247 5 : GDALCopyWords64(&dfSrcNoDataValue, GDT_Float64, 0, pabyWorkData,
248 : eOperDataType, static_cast<int>(nPixelOffset),
249 5 : static_cast<size_t>(nExtraXSize) * nExtraYSize);
250 : }
251 : else
252 : {
253 39 : memset(pabyWorkData, 0, nLineOffset * nExtraYSize);
254 : }
255 :
256 : /* -------------------------------------------------------------------- */
257 : /* Allocate the output buffer in the same dimensions as the work */
258 : /* buffer. This allows the filter process to write edge pixels */
259 : /* if needed for two-pass (separable) filtering. */
260 : /* -------------------------------------------------------------------- */
261 : GByte *pabyOutData = static_cast<GByte *>(
262 44 : VSI_MALLOC3_VERBOSE(nExtraXSize, nExtraYSize, nPixelOffset));
263 44 : if (pabyOutData == nullptr)
264 : {
265 0 : CPLFree(pabyWorkData);
266 0 : return CE_Failure;
267 : }
268 :
269 : /* -------------------------------------------------------------------- */
270 : /* Figure out the extended window that we want to load. Note */
271 : /* that we keep track of the file window as well as the amount */
272 : /* we will need to edge fill past the edge of the source dataset. */
273 : /* -------------------------------------------------------------------- */
274 44 : int nFileXOff = nReqXOff - m_nExtraEdgePixels;
275 44 : int nFileYOff = nReqYOff - m_nExtraEdgePixels;
276 44 : int nFileXSize = nExtraXSize;
277 44 : int nFileYSize = nExtraYSize;
278 :
279 44 : int nTopFill = 0;
280 44 : int nLeftFill = 0;
281 44 : int nRightFill = 0;
282 44 : int nBottomFill = 0;
283 :
284 44 : if (nFileXOff < 0)
285 : {
286 44 : nLeftFill = -nFileXOff;
287 44 : nFileXOff = 0;
288 44 : nFileXSize -= nLeftFill;
289 : }
290 :
291 44 : if (nFileYOff < 0)
292 : {
293 35 : nTopFill = -nFileYOff;
294 35 : nFileYOff = 0;
295 35 : nFileYSize -= nTopFill;
296 : }
297 :
298 44 : if (nFileXOff + nFileXSize > l_band->GetXSize())
299 : {
300 44 : nRightFill = nFileXOff + nFileXSize - l_band->GetXSize();
301 44 : nFileXSize -= nRightFill;
302 : }
303 :
304 44 : if (nFileYOff + nFileYSize > l_band->GetYSize())
305 : {
306 35 : nBottomFill = nFileYOff + nFileYSize - l_band->GetYSize();
307 35 : nFileYSize -= nBottomFill;
308 : }
309 :
310 : /* -------------------------------------------------------------------- */
311 : /* Load the data. */
312 : /* -------------------------------------------------------------------- */
313 : {
314 : GDALRasterIOExtraArg sExtraArgs;
315 44 : INIT_RASTERIO_EXTRA_ARG(sExtraArgs);
316 : const bool bIsComplex =
317 44 : CPL_TO_BOOL(GDALDataTypeIsComplex(eOperDataType));
318 44 : const CPLErr eErr = VRTComplexSource::RasterIOInternal<float>(
319 : l_band, eVRTBandDataType, nFileXOff, nFileYOff, nFileXSize,
320 : nFileYSize,
321 44 : pabyWorkData + nLineOffset * nTopFill + nPixelOffset * nLeftFill,
322 : nFileXSize, nFileYSize, eOperDataType, nPixelOffset, nLineOffset,
323 : &sExtraArgs, bIsComplex ? GDT_CFloat32 : GDT_Float32,
324 : oWorkingState);
325 :
326 44 : if (eErr != CE_None)
327 : {
328 0 : VSIFree(pabyWorkData);
329 0 : VSIFree(pabyOutData);
330 :
331 0 : return eErr;
332 : }
333 : }
334 :
335 : /* -------------------------------------------------------------------- */
336 : /* Fill in missing areas. Note that we replicate the edge */
337 : /* valid values out. We don't using "mirroring" which might be */
338 : /* more suitable for some times of filters. We also don't mark */
339 : /* these pixels as "nodata" though perhaps we should. */
340 : /* -------------------------------------------------------------------- */
341 44 : if (nLeftFill != 0 || nRightFill != 0)
342 : {
343 2191 : for (int i = nTopFill; i < nExtraYSize - nBottomFill; i++)
344 : {
345 2147 : if (nLeftFill != 0)
346 2147 : GDALCopyWords(
347 2147 : pabyWorkData + nPixelOffset * nLeftFill + i * nLineOffset,
348 2147 : eOperDataType, 0, pabyWorkData + i * nLineOffset,
349 : eOperDataType, static_cast<int>(nPixelOffset), nLeftFill);
350 :
351 2147 : if (nRightFill != 0)
352 2147 : GDALCopyWords(pabyWorkData + i * nLineOffset +
353 2147 : nPixelOffset * (nExtraXSize - nRightFill - 1),
354 : eOperDataType, 0,
355 2147 : pabyWorkData + i * nLineOffset +
356 2147 : nPixelOffset * (nExtraXSize - nRightFill),
357 : eOperDataType, static_cast<int>(nPixelOffset),
358 : nRightFill);
359 : }
360 : }
361 :
362 86 : for (int i = 0; i < nTopFill; i++)
363 : {
364 42 : memcpy(pabyWorkData + i * nLineOffset,
365 42 : pabyWorkData + nTopFill * nLineOffset, nLineOffset);
366 : }
367 :
368 86 : for (int i = nExtraYSize - nBottomFill; i < nExtraYSize; i++)
369 : {
370 42 : memcpy(pabyWorkData + i * nLineOffset,
371 42 : pabyWorkData + (nExtraYSize - nBottomFill - 1) * nLineOffset,
372 : nLineOffset);
373 : }
374 :
375 : /* -------------------------------------------------------------------- */
376 : /* Filter the data. */
377 : /* -------------------------------------------------------------------- */
378 88 : const CPLErr eErr = FilterData(nExtraXSize, nExtraYSize, eOperDataType,
379 44 : pabyWorkData, pabyOutData);
380 :
381 44 : VSIFree(pabyWorkData);
382 44 : if (eErr != CE_None)
383 : {
384 0 : VSIFree(pabyOutData);
385 :
386 0 : return eErr;
387 : }
388 :
389 : /* -------------------------------------------------------------------- */
390 : /* Copy from work buffer to target buffer. */
391 : /* -------------------------------------------------------------------- */
392 44 : GByte *pabySrcRow =
393 44 : pabyOutData + (nLineOffset + nPixelOffset) * m_nExtraEdgePixels;
394 44 : GByte *pabyDstRow = reinterpret_cast<GByte *>(pData);
395 :
396 2173 : for (int i = 0; i < nOutYSize;
397 2129 : i++, pabySrcRow += nLineOffset, pabyDstRow += nLineSpace)
398 : {
399 2129 : GDALCopyWords(pabySrcRow, eOperDataType, static_cast<int>(nPixelOffset),
400 : pabyDstRow, eBufType, static_cast<int>(nPixelSpace),
401 : nOutXSize);
402 : }
403 :
404 44 : VSIFree(pabyOutData);
405 :
406 44 : return CE_None;
407 : }
408 :
409 : /************************************************************************/
410 : /* ==================================================================== */
411 : /* VRTKernelFilteredSource */
412 : /* ==================================================================== */
413 : /************************************************************************/
414 :
415 : /************************************************************************/
416 : /* VRTKernelFilteredSource() */
417 : /************************************************************************/
418 :
419 42 : VRTKernelFilteredSource::VRTKernelFilteredSource()
420 : {
421 42 : GDALDataType aeSupTypes[] = {GDT_Float32};
422 42 : SetFilteringDataTypesSupported(1, aeSupTypes);
423 42 : }
424 :
425 : /************************************************************************/
426 : /* GetType() */
427 : /************************************************************************/
428 :
429 83 : const char *VRTKernelFilteredSource::GetType() const
430 : {
431 : static const char *TYPE = "KernelFilteredSource";
432 83 : return TYPE;
433 : }
434 :
435 : /************************************************************************/
436 : /* SetNormalized() */
437 : /************************************************************************/
438 :
439 40 : void VRTKernelFilteredSource::SetNormalized(bool bNormalizedIn)
440 :
441 : {
442 40 : m_bNormalized = bNormalizedIn;
443 40 : }
444 :
445 : /************************************************************************/
446 : /* SetKernel() */
447 : /************************************************************************/
448 :
449 : CPLErr
450 40 : VRTKernelFilteredSource::SetKernel(int nNewKernelSize, bool bSeparable,
451 : const std::vector<double> &adfNewCoefs)
452 :
453 : {
454 40 : if (nNewKernelSize < 1 || (nNewKernelSize % 2) != 1)
455 : {
456 0 : CPLError(CE_Failure, CPLE_AppDefined,
457 : "Illegal filtering kernel size %d, "
458 : "must be odd positive number.",
459 : nNewKernelSize);
460 0 : return CE_Failure;
461 : }
462 40 : if (adfNewCoefs.size() !=
463 40 : static_cast<size_t>(nNewKernelSize) * (bSeparable ? 1 : nNewKernelSize))
464 : {
465 0 : CPLError(CE_Failure, CPLE_AppDefined,
466 : "adfNewCoefs[] is not of expected size");
467 0 : return CE_Failure;
468 : }
469 :
470 40 : m_nKernelSize = nNewKernelSize;
471 40 : m_bSeparable = bSeparable;
472 40 : m_adfKernelCoefs = adfNewCoefs;
473 :
474 40 : SetExtraEdgePixels((nNewKernelSize - 1) / 2);
475 :
476 40 : return CE_None;
477 : }
478 :
479 : /************************************************************************/
480 : /* FilterData() */
481 : /************************************************************************/
482 :
483 44 : CPLErr VRTKernelFilteredSource::FilterData(int nXSize, int nYSize,
484 : GDALDataType eType,
485 : GByte *pabySrcData,
486 : GByte *pabyDstData)
487 :
488 : {
489 : /* -------------------------------------------------------------------- */
490 : /* Validate data type. */
491 : /* -------------------------------------------------------------------- */
492 44 : if (eType != GDT_Float32)
493 : {
494 0 : CPLError(CE_Failure, CPLE_AppDefined,
495 : "Unsupported data type (%s) in "
496 : "VRTKernelFilteredSource::FilterData()",
497 : GDALGetDataTypeName(eType));
498 0 : return CE_Failure;
499 : }
500 :
501 44 : CPLAssert(m_nExtraEdgePixels * 2 + 1 == m_nKernelSize ||
502 : (m_nKernelSize == 0 && m_nExtraEdgePixels == 0));
503 :
504 44 : const bool bMin = m_function == "min";
505 44 : const bool bMax = m_function == "max";
506 44 : const bool bStdDev = m_function == "stddev";
507 44 : const bool bMedian = m_function == "median";
508 44 : const bool bMode = m_function == "mode";
509 :
510 : /* -------------------------------------------------------------------- */
511 : /* Float32 case. */
512 : /* -------------------------------------------------------------------- */
513 : // if( eType == GDT_Float32 )
514 : {
515 44 : float *pafSrcData = reinterpret_cast<float *>(pabySrcData);
516 44 : float *pafDstData = reinterpret_cast<float *>(pabyDstData);
517 :
518 44 : int bHasNoData = FALSE;
519 44 : auto l_poBand = GetRasterBand();
520 44 : if (!l_poBand)
521 0 : return CE_Failure;
522 : const float fNoData =
523 44 : static_cast<float>(l_poBand->GetNoDataValue(&bHasNoData));
524 :
525 44 : const int nAxisCount = m_bSeparable ? 2 : 1;
526 :
527 89 : for (int nAxis = 0; nAxis < nAxisCount; ++nAxis)
528 : {
529 45 : const int nISize = nAxis == 0 ? nYSize : nXSize;
530 45 : const int nJSize = nAxis == 0 ? nXSize : nYSize;
531 45 : const int nIStride = nAxis == 0 ? nXSize : 1;
532 45 : const int nJStride = nAxis == 0 ? 1 : nXSize;
533 :
534 45 : const int nIMin = m_nExtraEdgePixels;
535 45 : const int nIMax = nISize - m_nExtraEdgePixels;
536 45 : const int nJMin = (m_bSeparable ? 0 : m_nExtraEdgePixels);
537 45 : const int nJMax = nJSize - (m_bSeparable ? 0 : m_nExtraEdgePixels);
538 :
539 7048 : for (GPtrDiff_t iJ = nJMin; iJ < nJMax; ++iJ)
540 : {
541 7003 : if (nAxis == 1)
542 62 : memcpy(pafSrcData + iJ * nJStride,
543 62 : pafDstData + iJ * nJStride, sizeof(float) * nXSize);
544 :
545 768120 : for (int iI = nIMin; iI < nIMax; ++iI)
546 : {
547 761117 : const GPtrDiff_t iIndex =
548 761117 : static_cast<GPtrDiff_t>(iI) * nIStride + iJ * nJStride;
549 :
550 761117 : if (bHasNoData && pafSrcData[iIndex] == fNoData)
551 : {
552 204 : pafDstData[iIndex] = fNoData;
553 204 : continue;
554 : }
555 :
556 760913 : double dfSum = 0.0, dfKernSum = 0.0;
557 760913 : size_t nValidCount = 0;
558 : double dfRes =
559 1521820 : bMin ? std::numeric_limits<double>::infinity()
560 760904 : : bMax ? -std::numeric_limits<double>::infinity()
561 760913 : : 0.0;
562 760913 : double dfMean = 0.0;
563 760913 : double dfM2 = 0.0;
564 1521830 : std::vector<double> adfVals;
565 760913 : std::map<double, size_t> mapValToCount;
566 760913 : size_t maxCount = 0;
567 :
568 3107250 : for (GPtrDiff_t iII = -m_nExtraEdgePixels, iK = 0;
569 3107250 : iII <= m_nExtraEdgePixels; ++iII)
570 : {
571 9232160 : for (GPtrDiff_t iJJ =
572 2346340 : (m_bSeparable ? 0 : -m_nExtraEdgePixels);
573 9232160 : iJJ <= (m_bSeparable ? 0 : m_nExtraEdgePixels);
574 : ++iJJ, ++iK)
575 : {
576 6885820 : const float *pfData = pafSrcData + iIndex +
577 6885820 : iII * nIStride +
578 6885820 : iJJ * nJStride;
579 6885820 : if (bHasNoData &&
580 6885820 : (*pfData == fNoData || std::isnan(*pfData)))
581 : {
582 976182 : continue;
583 : }
584 6884960 : if (m_adfKernelCoefs[iK] == 0.0)
585 : {
586 975326 : continue;
587 : }
588 11819300 : const double dfVal = static_cast<double>(*pfData) *
589 5909640 : m_adfKernelCoefs[iK];
590 5909640 : ++nValidCount;
591 :
592 5909640 : if (bMax)
593 : {
594 282 : if (dfVal > dfRes)
595 175 : dfRes = dfVal;
596 : }
597 5909350 : else if (bMin)
598 : {
599 81 : if (dfVal < dfRes)
600 9 : dfRes = dfVal;
601 : }
602 5909270 : else if (bStdDev)
603 : {
604 81 : const double dfDelta = dfVal - dfMean;
605 81 : dfMean += dfDelta / nValidCount;
606 81 : dfM2 += dfDelta * (dfVal - dfMean);
607 : }
608 5909190 : else if (bMedian)
609 : {
610 148 : adfVals.push_back(dfVal);
611 : }
612 5909040 : else if (bMode)
613 : {
614 81 : const size_t nCountVal = ++mapValToCount[dfVal];
615 81 : if (nCountVal > maxCount)
616 : {
617 26 : maxCount = nCountVal;
618 26 : dfRes = dfVal;
619 : }
620 : }
621 : else
622 : {
623 5908960 : dfSum += dfVal;
624 5908960 : dfKernSum += m_adfKernelCoefs[iK];
625 : }
626 : }
627 : }
628 :
629 : float fResult;
630 760913 : if (bMax || bMin || bMode)
631 : {
632 51 : fResult = nValidCount ? static_cast<float>(dfRes)
633 0 : : bHasNoData ? fNoData
634 0 : : 0.0f;
635 : }
636 760862 : else if (bStdDev)
637 : {
638 9 : fResult =
639 : nValidCount
640 9 : ? static_cast<float>(sqrt(
641 9 : dfM2 / static_cast<double>(nValidCount)))
642 0 : : bHasNoData ? fNoData
643 0 : : 0.0f;
644 : }
645 760853 : else if (bMedian)
646 : {
647 17 : if (!adfVals.empty())
648 : {
649 17 : if ((adfVals.size() % 2) == 1)
650 : {
651 32 : std::nth_element(adfVals.begin(),
652 0 : adfVals.begin() +
653 16 : adfVals.size() / 2,
654 : adfVals.end());
655 16 : dfRes = adfVals[adfVals.size() / 2];
656 : }
657 : else
658 : {
659 2 : std::nth_element(adfVals.begin(),
660 0 : adfVals.begin() +
661 1 : adfVals.size() / 2 - 1,
662 : adfVals.end());
663 1 : dfRes = adfVals[adfVals.size() / 2 - 1];
664 2 : std::nth_element(adfVals.begin(),
665 0 : adfVals.begin() +
666 1 : adfVals.size() / 2,
667 : adfVals.end());
668 1 : dfRes =
669 1 : (dfRes + adfVals[adfVals.size() / 2]) / 2;
670 : }
671 17 : fResult = static_cast<float>(dfRes);
672 : }
673 : else
674 0 : fResult = bHasNoData ? fNoData : 0.0f;
675 : }
676 760836 : else if (!m_bNormalized)
677 511627 : fResult = static_cast<float>(dfSum);
678 249209 : else if (nValidCount == 0 || dfKernSum == 0.0)
679 0 : fResult = bHasNoData ? fNoData : 0.0f;
680 : else
681 249209 : fResult = static_cast<float>(dfSum / dfKernSum);
682 760913 : pafDstData[iIndex] = fResult;
683 : }
684 : }
685 : }
686 : }
687 :
688 44 : return CE_None;
689 : }
690 :
691 : /************************************************************************/
692 : /* XMLInit() */
693 : /************************************************************************/
694 :
695 : CPLErr
696 14 : VRTKernelFilteredSource::XMLInit(const CPLXMLNode *psTree,
697 : const char *pszVRTPath,
698 : VRTMapSharedResources &oMapSharedSources)
699 :
700 : {
701 : {
702 : const CPLErr eErr =
703 14 : VRTFilteredSource::XMLInit(psTree, pszVRTPath, oMapSharedSources);
704 14 : if (eErr != CE_None)
705 0 : return eErr;
706 : }
707 :
708 14 : const int nNewKernelSize = atoi(CPLGetXMLValue(psTree, "Kernel.Size", "0"));
709 :
710 14 : if (nNewKernelSize == 0)
711 0 : return CE_None;
712 : // To prevent a integer overflow when computing nNewKernelSize *
713 : // nNewKernelSize
714 27 : if (nNewKernelSize < 0 ||
715 13 : nNewKernelSize > static_cast<int>(std::sqrt(static_cast<double>(
716 13 : std::numeric_limits<int>::max()))))
717 : {
718 2 : CPLError(CE_Failure, CPLE_IllegalArg,
719 : "Invalid value for kernel size: %d", nNewKernelSize);
720 2 : return CE_Failure;
721 : }
722 :
723 : const CPLStringList aosCoefItems(
724 24 : CSLTokenizeString(CPLGetXMLValue(psTree, "Kernel.Coefs", "")));
725 :
726 12 : const int nCoefs = aosCoefItems.size();
727 :
728 12 : const bool bSquare = nCoefs == nNewKernelSize * nNewKernelSize;
729 12 : const bool bSeparable = nCoefs == nNewKernelSize && nCoefs != 1;
730 :
731 12 : if (!bSquare && !bSeparable)
732 : {
733 0 : CPLError(CE_Failure, CPLE_AppDefined,
734 : "Got wrong number of filter kernel coefficients (%s). "
735 : "Expected %d or %d, got %d.",
736 : CPLGetXMLValue(psTree, "Kernel.Coefs", ""),
737 : nNewKernelSize * nNewKernelSize, nNewKernelSize, nCoefs);
738 0 : return CE_Failure;
739 : }
740 :
741 24 : std::vector<double> adfNewCoefs;
742 12 : adfNewCoefs.reserve(nCoefs);
743 122 : for (int i = 0; i < nCoefs; i++)
744 110 : adfNewCoefs.push_back(CPLAtof(aosCoefItems[i]));
745 :
746 12 : const CPLErr eErr = SetKernel(nNewKernelSize, bSeparable, adfNewCoefs);
747 12 : if (eErr == CE_None)
748 : {
749 12 : SetNormalized(atoi(CPLGetXMLValue(psTree, "Kernel.normalized", "0")) !=
750 : 0);
751 : }
752 :
753 12 : const char *pszFunction = CPLGetXMLValue(psTree, "Function", nullptr);
754 12 : if (pszFunction)
755 : {
756 0 : if (!EQUAL(pszFunction, "max") && !EQUAL(pszFunction, "min") &&
757 0 : !EQUAL(pszFunction, "median") && !EQUAL(pszFunction, "mode") &&
758 0 : !EQUAL(pszFunction, "stddev"))
759 : {
760 0 : CPLError(CE_Failure, CPLE_IllegalArg, "Unsupported function: %s",
761 : pszFunction);
762 0 : return CE_Failure;
763 : }
764 0 : SetFunction(pszFunction);
765 : }
766 :
767 12 : return eErr;
768 : }
769 :
770 : /************************************************************************/
771 : /* SerializeToXML() */
772 : /************************************************************************/
773 :
774 2 : CPLXMLNode *VRTKernelFilteredSource::SerializeToXML(const char *pszVRTPath)
775 :
776 : {
777 2 : CPLXMLNode *psSrc = VRTFilteredSource::SerializeToXML(pszVRTPath);
778 :
779 2 : if (psSrc == nullptr)
780 0 : return nullptr;
781 :
782 2 : CPLFree(psSrc->pszValue);
783 2 : psSrc->pszValue = CPLStrdup("KernelFilteredSource");
784 :
785 2 : if (m_nKernelSize == 0)
786 0 : return psSrc;
787 :
788 2 : CPLXMLNode *psKernel = CPLCreateXMLNode(psSrc, CXT_Element, "Kernel");
789 :
790 2 : CPLCreateXMLNode(CPLCreateXMLNode(psKernel, CXT_Attribute, "normalized"),
791 2 : CXT_Text, m_bNormalized ? "1" : "0");
792 :
793 2 : std::string osCoefs;
794 14 : for (auto dfVal : m_adfKernelCoefs)
795 : {
796 12 : if (!osCoefs.empty())
797 10 : osCoefs += ' ';
798 12 : osCoefs += CPLSPrintf("%.8g", dfVal);
799 : }
800 :
801 2 : CPLSetXMLValue(psKernel, "Size", CPLSPrintf("%d", m_nKernelSize));
802 2 : CPLSetXMLValue(psKernel, "Coefs", osCoefs.c_str());
803 :
804 2 : if (!m_function.empty())
805 0 : CPLCreateXMLElementAndValue(psSrc, "Function", m_function.c_str());
806 :
807 2 : return psSrc;
808 : }
809 :
810 : /************************************************************************/
811 : /* VRTParseFilterSources() */
812 : /************************************************************************/
813 :
814 14 : VRTSource *VRTParseFilterSources(const CPLXMLNode *psChild,
815 : const char *pszVRTPath,
816 : VRTMapSharedResources &oMapSharedSources)
817 :
818 : {
819 14 : if (EQUAL(psChild->pszValue, "KernelFilteredSource"))
820 : {
821 14 : VRTSource *poSrc = new VRTKernelFilteredSource();
822 14 : if (poSrc->XMLInit(psChild, pszVRTPath, oMapSharedSources) == CE_None)
823 12 : return poSrc;
824 :
825 2 : delete poSrc;
826 : }
827 :
828 2 : return nullptr;
829 : }
830 :
831 : /*! @endcond */
|
