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