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 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "cpl_port.h"
31 : #include "vrtdataset.h"
32 :
33 : #include <cmath>
34 : #include <cstddef>
35 : #include <cstdlib>
36 : #include <cstring>
37 : #include <limits>
38 :
39 : #include "cpl_conv.h"
40 : #include "cpl_error.h"
41 : #include "cpl_minixml.h"
42 : #include "cpl_string.h"
43 : #include "cpl_vsi.h"
44 : #include "gdal.h"
45 : #include "gdal_priv.h"
46 :
47 : /*! @cond Doxygen_Suppress */
48 :
49 : /************************************************************************/
50 : /* ==================================================================== */
51 : /* VRTFilteredSource */
52 : /* ==================================================================== */
53 : /************************************************************************/
54 :
55 : /************************************************************************/
56 : /* VRTFilteredSource() */
57 : /************************************************************************/
58 :
59 13 : VRTFilteredSource::VRTFilteredSource()
60 13 : : m_nSupportedTypesCount(1), m_nExtraEdgePixels(0)
61 : {
62 273 : for (size_t i = 0; i < CPL_ARRAYSIZE(m_aeSupportedTypes); ++i)
63 260 : m_aeSupportedTypes[i] = GDT_Unknown;
64 :
65 13 : m_aeSupportedTypes[0] = GDT_Float32;
66 13 : }
67 :
68 : /************************************************************************/
69 : /* ~VRTFilteredSource() */
70 : /************************************************************************/
71 :
72 13 : VRTFilteredSource::~VRTFilteredSource()
73 : {
74 13 : }
75 :
76 : /************************************************************************/
77 : /* SetExtraEdgePixels() */
78 : /************************************************************************/
79 :
80 11 : void VRTFilteredSource::SetExtraEdgePixels(int nEdgePixels)
81 :
82 : {
83 11 : m_nExtraEdgePixels = nEdgePixels;
84 11 : }
85 :
86 : /************************************************************************/
87 : /* SetFilteringDataTypesSupported() */
88 : /************************************************************************/
89 :
90 13 : void VRTFilteredSource::SetFilteringDataTypesSupported(int nTypeCount,
91 : GDALDataType *paeTypes)
92 :
93 : {
94 13 : if (nTypeCount >
95 : static_cast<int>(sizeof(m_aeSupportedTypes) / sizeof(GDALDataType)))
96 : {
97 0 : CPLAssert(false);
98 : nTypeCount =
99 : static_cast<int>(sizeof(m_aeSupportedTypes) / sizeof(GDALDataType));
100 : }
101 :
102 13 : m_nSupportedTypesCount = nTypeCount;
103 13 : memcpy(m_aeSupportedTypes, paeTypes, sizeof(GDALDataType) * nTypeCount);
104 13 : }
105 :
106 : /************************************************************************/
107 : /* IsTypeSupported() */
108 : /************************************************************************/
109 :
110 28 : int VRTFilteredSource::IsTypeSupported(GDALDataType eTestType) const
111 :
112 : {
113 56 : for (int i = 0; i < m_nSupportedTypesCount; i++)
114 : {
115 28 : if (eTestType == m_aeSupportedTypes[i])
116 0 : return TRUE;
117 : }
118 :
119 28 : return FALSE;
120 : }
121 :
122 : /************************************************************************/
123 : /* RasterIO() */
124 : /************************************************************************/
125 :
126 14 : CPLErr VRTFilteredSource::RasterIO(GDALDataType eVRTBandDataType, int nXOff,
127 : int nYOff, int nXSize, int nYSize,
128 : void *pData, int nBufXSize, int nBufYSize,
129 : GDALDataType eBufType, GSpacing nPixelSpace,
130 : GSpacing nLineSpace,
131 : GDALRasterIOExtraArg *psExtraArg,
132 : WorkingState &oWorkingState)
133 :
134 : {
135 : /* -------------------------------------------------------------------- */
136 : /* For now we don't support filtered access to non-full */
137 : /* resolution requests. Just collect the data directly without */
138 : /* any operator. */
139 : /* -------------------------------------------------------------------- */
140 14 : if (nBufXSize != nXSize || nBufYSize != nYSize)
141 : {
142 0 : return VRTComplexSource::RasterIO(
143 : eVRTBandDataType, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize,
144 : nBufYSize, eBufType, nPixelSpace, nLineSpace, psExtraArg,
145 0 : oWorkingState);
146 : }
147 :
148 14 : double dfXOff = nXOff;
149 14 : double dfYOff = nYOff;
150 14 : double dfXSize = nXSize;
151 14 : double dfYSize = nYSize;
152 14 : if (psExtraArg != nullptr && psExtraArg->bFloatingPointWindowValidity)
153 : {
154 4 : dfXOff = psExtraArg->dfXOff;
155 4 : dfYOff = psExtraArg->dfYOff;
156 4 : dfXSize = psExtraArg->dfXSize;
157 4 : dfYSize = psExtraArg->dfYSize;
158 : }
159 :
160 : // The window we will actually request from the source raster band.
161 14 : double dfReqXOff = 0.0;
162 14 : double dfReqYOff = 0.0;
163 14 : double dfReqXSize = 0.0;
164 14 : double dfReqYSize = 0.0;
165 14 : int nReqXOff = 0;
166 14 : int nReqYOff = 0;
167 14 : int nReqXSize = 0;
168 14 : int nReqYSize = 0;
169 :
170 : // The window we will actual set _within_ the pData buffer.
171 14 : int nOutXOff = 0;
172 14 : int nOutYOff = 0;
173 14 : int nOutXSize = 0;
174 14 : int nOutYSize = 0;
175 :
176 14 : bool bError = false;
177 14 : if (!GetSrcDstWindow(dfXOff, dfYOff, dfXSize, dfYSize, nBufXSize, nBufYSize,
178 : &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize,
179 : &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
180 : &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize, bError))
181 : {
182 0 : return bError ? CE_Failure : CE_None;
183 : }
184 :
185 14 : pData = reinterpret_cast<GByte *>(pData) + nPixelSpace * nOutXOff +
186 14 : nLineSpace * nOutYOff;
187 :
188 : /* -------------------------------------------------------------------- */
189 : /* Determine the data type we want to request. We try to match */
190 : /* the source or destination request, and if both those fail we */
191 : /* fallback to the first supported type at least as expressive */
192 : /* as the request. */
193 : /* -------------------------------------------------------------------- */
194 14 : GDALDataType eOperDataType = GDT_Unknown;
195 :
196 14 : if (IsTypeSupported(eBufType))
197 0 : eOperDataType = eBufType;
198 :
199 14 : auto l_band = GetRasterBand();
200 14 : if (!l_band)
201 0 : return CE_Failure;
202 :
203 28 : if (eOperDataType == GDT_Unknown &&
204 14 : IsTypeSupported(l_band->GetRasterDataType()))
205 0 : eOperDataType = l_band->GetRasterDataType();
206 :
207 14 : if (eOperDataType == GDT_Unknown)
208 : {
209 28 : for (int i = 0; i < m_nSupportedTypesCount; i++)
210 : {
211 14 : if (GDALDataTypeUnion(m_aeSupportedTypes[i], eBufType) ==
212 14 : m_aeSupportedTypes[i])
213 : {
214 0 : eOperDataType = m_aeSupportedTypes[i];
215 : }
216 : }
217 : }
218 :
219 14 : if (eOperDataType == GDT_Unknown)
220 : {
221 14 : eOperDataType = m_aeSupportedTypes[0];
222 :
223 14 : for (int i = 1; i < m_nSupportedTypesCount; i++)
224 : {
225 0 : if (GDALGetDataTypeSize(m_aeSupportedTypes[i]) >
226 0 : GDALGetDataTypeSize(eOperDataType))
227 : {
228 0 : eOperDataType = m_aeSupportedTypes[i];
229 : }
230 : }
231 : }
232 :
233 : /* -------------------------------------------------------------------- */
234 : /* Allocate the buffer of data into which our imagery will be */
235 : /* read, with the extra edge pixels as well. This will be the */
236 : /* source data fed into the filter. */
237 : /* -------------------------------------------------------------------- */
238 14 : if (nOutXSize > INT_MAX - 2 * m_nExtraEdgePixels ||
239 14 : nOutYSize > INT_MAX - 2 * m_nExtraEdgePixels)
240 : {
241 0 : return CE_Failure;
242 : }
243 14 : const int nExtraXSize = nOutXSize + 2 * m_nExtraEdgePixels;
244 14 : const int nExtraYSize = nOutYSize + 2 * m_nExtraEdgePixels;
245 :
246 14 : GByte *pabyWorkData = static_cast<GByte *>(VSI_MALLOC3_VERBOSE(
247 : nExtraXSize, nExtraYSize, GDALGetDataTypeSizeBytes(eOperDataType)));
248 :
249 14 : if (pabyWorkData == nullptr)
250 : {
251 0 : return CE_Failure;
252 : }
253 :
254 14 : const GPtrDiff_t nPixelOffset = GDALGetDataTypeSizeBytes(eOperDataType);
255 14 : const GPtrDiff_t nLineOffset = nPixelOffset * nExtraXSize;
256 :
257 14 : memset(pabyWorkData, 0, nLineOffset * nExtraYSize);
258 :
259 : /* -------------------------------------------------------------------- */
260 : /* Allocate the output buffer in the same dimensions as the work */
261 : /* buffer. This allows the filter process to write edge pixels */
262 : /* if needed for two-pass (separable) filtering. */
263 : /* -------------------------------------------------------------------- */
264 : GByte *pabyOutData = static_cast<GByte *>(
265 14 : VSI_MALLOC3_VERBOSE(nExtraXSize, nExtraYSize, nPixelOffset));
266 14 : if (pabyOutData == nullptr)
267 : {
268 0 : CPLFree(pabyWorkData);
269 0 : return CE_Failure;
270 : }
271 :
272 : /* -------------------------------------------------------------------- */
273 : /* Figure out the extended window that we want to load. Note */
274 : /* that we keep track of the file window as well as the amount */
275 : /* we will need to edge fill past the edge of the source dataset. */
276 : /* -------------------------------------------------------------------- */
277 14 : int nFileXOff = nReqXOff - m_nExtraEdgePixels;
278 14 : int nFileYOff = nReqYOff - m_nExtraEdgePixels;
279 14 : int nFileXSize = nExtraXSize;
280 14 : int nFileYSize = nExtraYSize;
281 :
282 14 : int nTopFill = 0;
283 14 : int nLeftFill = 0;
284 14 : int nRightFill = 0;
285 14 : int nBottomFill = 0;
286 :
287 14 : if (nFileXOff < 0)
288 : {
289 14 : nLeftFill = -nFileXOff;
290 14 : nFileXOff = 0;
291 14 : nFileXSize -= nLeftFill;
292 : }
293 :
294 14 : if (nFileYOff < 0)
295 : {
296 8 : nTopFill = -nFileYOff;
297 8 : nFileYOff = 0;
298 8 : nFileYSize -= nTopFill;
299 : }
300 :
301 14 : if (nFileXOff + nFileXSize > l_band->GetXSize())
302 : {
303 14 : nRightFill = nFileXOff + nFileXSize - l_band->GetXSize();
304 14 : nFileXSize -= nRightFill;
305 : }
306 :
307 14 : if (nFileYOff + nFileYSize > l_band->GetYSize())
308 : {
309 8 : nBottomFill = nFileYOff + nFileYSize - l_band->GetYSize();
310 8 : nFileYSize -= nBottomFill;
311 : }
312 :
313 : /* -------------------------------------------------------------------- */
314 : /* Load the data. */
315 : /* -------------------------------------------------------------------- */
316 : {
317 : GDALRasterIOExtraArg sExtraArgs;
318 14 : INIT_RASTERIO_EXTRA_ARG(sExtraArgs);
319 : const bool bIsComplex =
320 14 : CPL_TO_BOOL(GDALDataTypeIsComplex(eOperDataType));
321 14 : const CPLErr eErr = VRTComplexSource::RasterIOInternal<float>(
322 : l_band, eVRTBandDataType, nFileXOff, nFileYOff, nFileXSize,
323 : nFileYSize,
324 14 : pabyWorkData + nLineOffset * nTopFill + nPixelOffset * nLeftFill,
325 : nFileXSize, nFileYSize, eOperDataType, nPixelOffset, nLineOffset,
326 : &sExtraArgs, bIsComplex ? GDT_CFloat32 : GDT_Float32,
327 : oWorkingState);
328 :
329 14 : if (eErr != CE_None)
330 : {
331 0 : VSIFree(pabyWorkData);
332 0 : VSIFree(pabyOutData);
333 :
334 0 : return eErr;
335 : }
336 : }
337 :
338 : /* -------------------------------------------------------------------- */
339 : /* Fill in missing areas. Note that we replicate the edge */
340 : /* valid values out. We don't using "mirroring" which might be */
341 : /* more suitable for some times of filters. We also don't mark */
342 : /* these pixels as "nodata" though perhaps we should. */
343 : /* -------------------------------------------------------------------- */
344 14 : if (nLeftFill != 0 || nRightFill != 0)
345 : {
346 1296 : for (int i = nTopFill; i < nExtraYSize - nBottomFill; i++)
347 : {
348 1282 : if (nLeftFill != 0)
349 1282 : GDALCopyWords(
350 1282 : pabyWorkData + nPixelOffset * nLeftFill + i * nLineOffset,
351 1282 : eOperDataType, 0, pabyWorkData + i * nLineOffset,
352 : eOperDataType, static_cast<int>(nPixelOffset), nLeftFill);
353 :
354 1282 : if (nRightFill != 0)
355 1282 : GDALCopyWords(pabyWorkData + i * nLineOffset +
356 1282 : nPixelOffset * (nExtraXSize - nRightFill - 1),
357 : eOperDataType, 0,
358 1282 : pabyWorkData + i * nLineOffset +
359 1282 : nPixelOffset * (nExtraXSize - nRightFill),
360 : eOperDataType, static_cast<int>(nPixelOffset),
361 : nRightFill);
362 : }
363 : }
364 :
365 27 : for (int i = 0; i < nTopFill; i++)
366 : {
367 13 : memcpy(pabyWorkData + i * nLineOffset,
368 13 : pabyWorkData + nTopFill * nLineOffset, nLineOffset);
369 : }
370 :
371 27 : for (int i = nExtraYSize - nBottomFill; i < nExtraYSize; i++)
372 : {
373 13 : memcpy(pabyWorkData + i * nLineOffset,
374 13 : pabyWorkData + (nExtraYSize - nBottomFill - 1) * nLineOffset,
375 : nLineOffset);
376 : }
377 :
378 : /* -------------------------------------------------------------------- */
379 : /* Filter the data. */
380 : /* -------------------------------------------------------------------- */
381 28 : const CPLErr eErr = FilterData(nExtraXSize, nExtraYSize, eOperDataType,
382 14 : pabyWorkData, pabyOutData);
383 :
384 14 : VSIFree(pabyWorkData);
385 14 : if (eErr != CE_None)
386 : {
387 0 : VSIFree(pabyOutData);
388 :
389 0 : return eErr;
390 : }
391 :
392 : /* -------------------------------------------------------------------- */
393 : /* Copy from work buffer to target buffer. */
394 : /* -------------------------------------------------------------------- */
395 14 : GByte *pabySrcRow =
396 14 : pabyOutData + (nLineOffset + nPixelOffset) * m_nExtraEdgePixels;
397 14 : GByte *pabyDstRow = reinterpret_cast<GByte *>(pData);
398 :
399 1284 : for (int i = 0; i < nOutYSize;
400 1270 : i++, pabySrcRow += nLineOffset, pabyDstRow += nLineSpace)
401 : {
402 1270 : GDALCopyWords(pabySrcRow, eOperDataType, static_cast<int>(nPixelOffset),
403 : pabyDstRow, eBufType, static_cast<int>(nPixelSpace),
404 : nOutXSize);
405 : }
406 :
407 14 : VSIFree(pabyOutData);
408 :
409 14 : return CE_None;
410 : }
411 :
412 : /************************************************************************/
413 : /* ==================================================================== */
414 : /* VRTKernelFilteredSource */
415 : /* ==================================================================== */
416 : /************************************************************************/
417 :
418 : /************************************************************************/
419 : /* VRTKernelFilteredSource() */
420 : /************************************************************************/
421 :
422 13 : VRTKernelFilteredSource::VRTKernelFilteredSource()
423 : : m_nKernelSize(0), m_bSeparable(FALSE), m_padfKernelCoefs(nullptr),
424 13 : m_bNormalized(FALSE)
425 : {
426 13 : GDALDataType aeSupTypes[] = {GDT_Float32};
427 13 : SetFilteringDataTypesSupported(1, aeSupTypes);
428 13 : }
429 :
430 : /************************************************************************/
431 : /* ~VRTKernelFilteredSource() */
432 : /************************************************************************/
433 :
434 26 : VRTKernelFilteredSource::~VRTKernelFilteredSource()
435 :
436 : {
437 13 : CPLFree(m_padfKernelCoefs);
438 26 : }
439 :
440 : /************************************************************************/
441 : /* SetNormalized() */
442 : /************************************************************************/
443 :
444 11 : void VRTKernelFilteredSource::SetNormalized(int bNormalizedIn)
445 :
446 : {
447 11 : m_bNormalized = bNormalizedIn;
448 11 : }
449 :
450 : /************************************************************************/
451 : /* SetKernel() */
452 : /************************************************************************/
453 :
454 11 : CPLErr VRTKernelFilteredSource::SetKernel(int nNewKernelSize, bool bSeparable,
455 : double *padfNewCoefs)
456 :
457 : {
458 11 : if (nNewKernelSize < 1 || (nNewKernelSize % 2) != 1)
459 : {
460 0 : CPLError(CE_Failure, CPLE_AppDefined,
461 : "Illegal filtering kernel size %d, "
462 : "must be odd positive number.",
463 : nNewKernelSize);
464 0 : return CE_Failure;
465 : }
466 :
467 11 : CPLFree(m_padfKernelCoefs);
468 11 : m_nKernelSize = nNewKernelSize;
469 11 : m_bSeparable = bSeparable;
470 :
471 11 : int nKernelBufferSize = m_nKernelSize * (m_bSeparable ? 1 : m_nKernelSize);
472 :
473 11 : m_padfKernelCoefs =
474 11 : static_cast<double *>(CPLMalloc(sizeof(double) * nKernelBufferSize));
475 11 : memcpy(m_padfKernelCoefs, padfNewCoefs, sizeof(double) * nKernelBufferSize);
476 :
477 11 : SetExtraEdgePixels((nNewKernelSize - 1) / 2);
478 :
479 11 : return CE_None;
480 : }
481 :
482 : /************************************************************************/
483 : /* FilterData() */
484 : /************************************************************************/
485 :
486 14 : CPLErr VRTKernelFilteredSource::FilterData(int nXSize, int nYSize,
487 : GDALDataType eType,
488 : GByte *pabySrcData,
489 : GByte *pabyDstData)
490 :
491 : {
492 : /* -------------------------------------------------------------------- */
493 : /* Validate data type. */
494 : /* -------------------------------------------------------------------- */
495 14 : if (eType != GDT_Float32)
496 : {
497 0 : CPLError(CE_Failure, CPLE_AppDefined,
498 : "Unsupported data type (%s) in "
499 : "VRTKernelFilteredSource::FilterData()",
500 : GDALGetDataTypeName(eType));
501 0 : return CE_Failure;
502 : }
503 :
504 14 : CPLAssert(m_nExtraEdgePixels * 2 + 1 == m_nKernelSize ||
505 : (m_nKernelSize == 0 && m_nExtraEdgePixels == 0));
506 :
507 : /* -------------------------------------------------------------------- */
508 : /* Float32 case. */
509 : /* -------------------------------------------------------------------- */
510 : // if( eType == GDT_Float32 )
511 : {
512 14 : float *pafSrcData = reinterpret_cast<float *>(pabySrcData);
513 14 : float *pafDstData = reinterpret_cast<float *>(pabyDstData);
514 :
515 14 : int bHasNoData = FALSE;
516 14 : auto l_poBand = GetRasterBand();
517 14 : if (!l_poBand)
518 0 : return CE_Failure;
519 : const float fNoData =
520 14 : static_cast<float>(l_poBand->GetNoDataValue(&bHasNoData));
521 :
522 14 : const int nAxisCount = m_bSeparable ? 2 : 1;
523 :
524 29 : for (int nAxis = 0; nAxis < nAxisCount; ++nAxis)
525 : {
526 15 : const int nISize = nAxis == 0 ? nYSize : nXSize;
527 15 : const int nJSize = nAxis == 0 ? nXSize : nYSize;
528 15 : const int nIStride = nAxis == 0 ? nXSize : 1;
529 15 : const int nJStride = nAxis == 0 ? 1 : nXSize;
530 :
531 15 : const int nIMin = m_nExtraEdgePixels;
532 15 : const int nIMax = nISize - m_nExtraEdgePixels;
533 15 : const int nJMin = (m_bSeparable ? 0 : m_nExtraEdgePixels);
534 15 : const int nJMax = nJSize - (m_bSeparable ? 0 : m_nExtraEdgePixels);
535 :
536 4359 : for (GPtrDiff_t iJ = nJMin; iJ < nJMax; ++iJ)
537 : {
538 4344 : if (nAxis == 1)
539 62 : memcpy(pafSrcData + iJ * nJStride,
540 62 : pafDstData + iJ * nJStride, sizeof(float) * nXSize);
541 :
542 520944 : for (int iI = nIMin; iI < nIMax; ++iI)
543 : {
544 516600 : const GPtrDiff_t iIndex =
545 516600 : static_cast<GPtrDiff_t>(iI) * nIStride + iJ * nJStride;
546 :
547 516600 : if (bHasNoData && pafSrcData[iIndex] == fNoData)
548 : {
549 200 : pafDstData[iIndex] = fNoData;
550 200 : continue;
551 : }
552 :
553 516400 : double dfSum = 0.0, dfKernSum = 0.0;
554 :
555 2127600 : for (GPtrDiff_t iII = -m_nExtraEdgePixels, iK = 0;
556 2127600 : iII <= m_nExtraEdgePixels; ++iII)
557 : {
558 6283600 : for (GPtrDiff_t iJJ =
559 1611200 : (m_bSeparable ? 0 : -m_nExtraEdgePixels);
560 6283600 : iJJ <= (m_bSeparable ? 0 : m_nExtraEdgePixels);
561 : ++iJJ, ++iK)
562 : {
563 4672400 : const float *pfData = pafSrcData + iIndex +
564 4672400 : iII * nIStride +
565 4672400 : iJJ * nJStride;
566 4672400 : if (bHasNoData && *pfData == fNoData)
567 836 : continue;
568 4671560 : dfSum += *pfData * m_padfKernelCoefs[iK];
569 4671560 : dfKernSum += m_padfKernelCoefs[iK];
570 : }
571 : }
572 :
573 : double fResult;
574 :
575 516400 : if (!m_bNormalized)
576 510000 : fResult = dfSum;
577 6400 : else if (dfKernSum == 0.0)
578 0 : fResult = 0.0;
579 : else
580 6400 : fResult = dfSum / dfKernSum;
581 516400 : pafDstData[iIndex] = static_cast<float>(fResult);
582 : }
583 : }
584 : }
585 : }
586 :
587 14 : return CE_None;
588 : }
589 :
590 : /************************************************************************/
591 : /* XMLInit() */
592 : /************************************************************************/
593 :
594 13 : CPLErr VRTKernelFilteredSource::XMLInit(
595 : const CPLXMLNode *psTree, const char *pszVRTPath,
596 : std::map<CPLString, GDALDataset *> &oMapSharedSources)
597 :
598 : {
599 : {
600 : const CPLErr eErr =
601 13 : VRTFilteredSource::XMLInit(psTree, pszVRTPath, oMapSharedSources);
602 13 : if (eErr != CE_None)
603 0 : return eErr;
604 : }
605 :
606 13 : const int nNewKernelSize = atoi(CPLGetXMLValue(psTree, "Kernel.Size", "0"));
607 :
608 13 : if (nNewKernelSize == 0)
609 0 : return CE_None;
610 : // To prevent a integer overflow when computing nNewKernelSize *
611 : // nNewKernelSize
612 25 : if (nNewKernelSize < 0 ||
613 12 : nNewKernelSize > static_cast<int>(std::sqrt(static_cast<double>(
614 12 : std::numeric_limits<int>::max()))))
615 : {
616 2 : CPLError(CE_Failure, CPLE_IllegalArg,
617 : "Invalid value for kernel size: %d", nNewKernelSize);
618 2 : return CE_Failure;
619 : }
620 :
621 : char **papszCoefItems =
622 11 : CSLTokenizeString(CPLGetXMLValue(psTree, "Kernel.Coefs", ""));
623 :
624 11 : const int nCoefs = CSLCount(papszCoefItems);
625 :
626 11 : const bool bSquare = nCoefs == nNewKernelSize * nNewKernelSize;
627 11 : const bool bSeparable = nCoefs == nNewKernelSize && nCoefs != 1;
628 :
629 11 : if (!bSquare && !bSeparable)
630 : {
631 0 : CSLDestroy(papszCoefItems);
632 0 : CPLError(CE_Failure, CPLE_AppDefined,
633 : "Got wrong number of filter kernel coefficients (%s). "
634 : "Expected %d or %d, got %d.",
635 : CPLGetXMLValue(psTree, "Kernel.Coefs", ""),
636 : nNewKernelSize * nNewKernelSize, nNewKernelSize, nCoefs);
637 0 : return CE_Failure;
638 : }
639 :
640 : double *padfNewCoefs =
641 11 : static_cast<double *>(CPLMalloc(sizeof(double) * nCoefs));
642 :
643 118 : for (int i = 0; i < nCoefs; i++)
644 107 : padfNewCoefs[i] = CPLAtof(papszCoefItems[i]);
645 :
646 11 : const CPLErr eErr = SetKernel(nNewKernelSize, bSeparable, padfNewCoefs);
647 :
648 11 : CPLFree(padfNewCoefs);
649 11 : CSLDestroy(papszCoefItems);
650 :
651 11 : SetNormalized(atoi(CPLGetXMLValue(psTree, "Kernel.normalized", "0")));
652 :
653 11 : return eErr;
654 : }
655 :
656 : /************************************************************************/
657 : /* SerializeToXML() */
658 : /************************************************************************/
659 :
660 1 : CPLXMLNode *VRTKernelFilteredSource::SerializeToXML(const char *pszVRTPath)
661 :
662 : {
663 1 : CPLXMLNode *psSrc = VRTFilteredSource::SerializeToXML(pszVRTPath);
664 :
665 1 : if (psSrc == nullptr)
666 0 : return nullptr;
667 :
668 1 : CPLFree(psSrc->pszValue);
669 1 : psSrc->pszValue = CPLStrdup("KernelFilteredSource");
670 :
671 1 : if (m_nKernelSize == 0)
672 0 : return psSrc;
673 :
674 1 : CPLXMLNode *psKernel = CPLCreateXMLNode(psSrc, CXT_Element, "Kernel");
675 :
676 1 : if (m_bNormalized)
677 0 : CPLCreateXMLNode(
678 : CPLCreateXMLNode(psKernel, CXT_Attribute, "normalized"), CXT_Text,
679 : "1");
680 : else
681 1 : CPLCreateXMLNode(
682 : CPLCreateXMLNode(psKernel, CXT_Attribute, "normalized"), CXT_Text,
683 : "0");
684 :
685 1 : const int nCoefCount = m_nKernelSize * m_nKernelSize;
686 1 : const size_t nBufLen = nCoefCount * 32;
687 1 : char *pszKernelCoefs = static_cast<char *>(CPLMalloc(nBufLen));
688 :
689 1 : strcpy(pszKernelCoefs, "");
690 10 : for (int iCoef = 0; iCoef < nCoefCount; iCoef++)
691 9 : CPLsnprintf(pszKernelCoefs + strlen(pszKernelCoefs),
692 9 : nBufLen - strlen(pszKernelCoefs), "%.8g ",
693 9 : m_padfKernelCoefs[iCoef]);
694 :
695 1 : CPLSetXMLValue(psKernel, "Size", CPLSPrintf("%d", m_nKernelSize));
696 1 : CPLSetXMLValue(psKernel, "Coefs", pszKernelCoefs);
697 :
698 1 : CPLFree(pszKernelCoefs);
699 :
700 1 : return psSrc;
701 : }
702 :
703 : /************************************************************************/
704 : /* VRTParseFilterSources() */
705 : /************************************************************************/
706 :
707 : VRTSource *
708 13 : VRTParseFilterSources(const CPLXMLNode *psChild, const char *pszVRTPath,
709 : std::map<CPLString, GDALDataset *> &oMapSharedSources)
710 :
711 : {
712 13 : if (EQUAL(psChild->pszValue, "KernelFilteredSource"))
713 : {
714 13 : VRTSource *poSrc = new VRTKernelFilteredSource();
715 13 : if (poSrc->XMLInit(psChild, pszVRTPath, oMapSharedSources) == CE_None)
716 11 : return poSrc;
717 :
718 2 : delete poSrc;
719 : }
720 :
721 2 : return nullptr;
722 : }
723 :
724 : /*! @endcond */
|
