Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: gdal "raster compare" subcommand
5 : * Author: Even Rouault <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdalalg_raster_compare.h"
14 :
15 : #include "cpl_conv.h"
16 : #include "cpl_vsi_virtual.h"
17 : #include "gdal_alg.h"
18 : #include "gdal_priv.h"
19 :
20 : #include <algorithm>
21 : #include <array>
22 : #include <cmath>
23 : #include <limits>
24 : #include <type_traits>
25 :
26 : #if defined(__x86_64__) || defined(_M_X64)
27 : #define USE_SSE2
28 : #include <emmintrin.h>
29 : #elif defined(USE_NEON_OPTIMIZATIONS)
30 : #define USE_SSE2
31 : #include "include_sse2neon.h"
32 : #endif
33 :
34 : //! @cond Doxygen_Suppress
35 :
36 : #ifndef _
37 : #define _(x) (x)
38 : #endif
39 :
40 : /************************************************************************/
41 : /* GDALRasterCompareAlgorithm::GDALRasterCompareAlgorithm() */
42 : /************************************************************************/
43 :
44 286 : GDALRasterCompareAlgorithm::GDALRasterCompareAlgorithm(bool standaloneStep)
45 : : GDALRasterPipelineStepAlgorithm(
46 : NAME, DESCRIPTION, HELP_URL,
47 0 : ConstructorOptions()
48 286 : .SetStandaloneStep(standaloneStep)
49 286 : .SetInputDatasetMaxCount(1)
50 572 : .SetInputDatasetHelpMsg(_("Input raster dataset"))
51 858 : .SetAddDefaultArguments(false))
52 : {
53 286 : if (standaloneStep)
54 : {
55 245 : AddProgressArg();
56 : }
57 : else
58 : {
59 41 : AddRasterHiddenInputDatasetArg();
60 : }
61 :
62 : auto &referenceDatasetArg = AddArg("reference", 0, _("Reference dataset"),
63 572 : &m_referenceDataset, GDAL_OF_RASTER)
64 286 : .SetPositional()
65 286 : .SetRequired();
66 :
67 286 : SetAutoCompleteFunctionForFilename(referenceDatasetArg, GDAL_OF_RASTER);
68 :
69 286 : if (standaloneStep)
70 : {
71 245 : AddRasterInputArgs(/* openForMixedRasterVector = */ false,
72 : /* hiddenForCLI = */ false);
73 : }
74 :
75 : AddArg("skip-all-optional", 0, _("Skip all optional comparisons"),
76 286 : &m_skipAllOptional);
77 286 : AddArg("skip-binary", 0, _("Skip binary file comparison"), &m_skipBinary);
78 286 : AddArg("skip-crs", 0, _("Skip CRS comparison"), &m_skipCRS);
79 : AddArg("skip-geotransform", 0, _("Skip geotransform comparison"),
80 286 : &m_skipGeotransform);
81 286 : AddArg("skip-overview", 0, _("Skip overview comparison"), &m_skipOverview);
82 286 : AddArg("skip-metadata", 0, _("Skip metadata comparison"), &m_skipMetadata);
83 286 : AddArg("skip-rpc", 0, _("Skip RPC metadata comparison"), &m_skipRPC);
84 : AddArg("skip-geolocation", 0, _("Skip Geolocation metadata comparison"),
85 286 : &m_skipGeolocation);
86 : AddArg("skip-subdataset", 0, _("Skip subdataset comparison"),
87 286 : &m_skipSubdataset);
88 :
89 286 : AddOutputStringArg(&m_output);
90 :
91 572 : AddArg("return-code", 0, _("Return code"), &m_retCode)
92 286 : .SetHiddenForCLI()
93 286 : .SetIsInput(false)
94 286 : .SetIsOutput(true);
95 286 : }
96 :
97 : /************************************************************************/
98 : /* GDALRasterCompareAlgorithm::BinaryComparison() */
99 : /************************************************************************/
100 :
101 12 : bool GDALRasterCompareAlgorithm::BinaryComparison(
102 : std::vector<std::string> &aosReport, GDALDataset *poRefDS,
103 : GDALDataset *poInputDS)
104 : {
105 12 : if (poRefDS->GetDescription()[0] == 0)
106 : {
107 1 : ReportError(
108 : CE_Warning, CPLE_AppDefined,
109 : "Reference dataset has no name. Skipping binary file comparison");
110 1 : return false;
111 : }
112 :
113 11 : auto poRefDrv = poRefDS->GetDriver();
114 11 : if (poRefDrv && EQUAL(poRefDrv->GetDescription(), "MEM"))
115 : {
116 1 : ReportError(CE_Warning, CPLE_AppDefined,
117 : "Reference dataset is a in-memory dataset. Skipping binary "
118 : "file comparison");
119 1 : return false;
120 : }
121 :
122 10 : if (poInputDS->GetDescription()[0] == 0)
123 : {
124 1 : ReportError(
125 : CE_Warning, CPLE_AppDefined,
126 : "Input dataset has no name. Skipping binary file comparison");
127 1 : return false;
128 : }
129 :
130 9 : auto poInputDrv = poInputDS->GetDriver();
131 9 : if (poInputDrv && EQUAL(poInputDrv->GetDescription(), "MEM"))
132 : {
133 1 : ReportError(CE_Warning, CPLE_AppDefined,
134 : "Input dataset is a in-memory dataset. Skipping binary "
135 : "file comparison");
136 1 : return false;
137 : }
138 :
139 16 : VSIVirtualHandleUniquePtr fpRef(VSIFOpenL(poRefDS->GetDescription(), "rb"));
140 : VSIVirtualHandleUniquePtr fpInput(
141 16 : VSIFOpenL(poInputDS->GetDescription(), "rb"));
142 8 : if (!fpRef)
143 : {
144 1 : ReportError(CE_Warning, CPLE_AppDefined,
145 : "Reference dataset '%s' is not a file. Skipping binary "
146 : "file comparison",
147 1 : poRefDS->GetDescription());
148 1 : return false;
149 : }
150 :
151 7 : if (!fpInput)
152 : {
153 1 : ReportError(
154 : CE_Warning, CPLE_AppDefined,
155 : "Input dataset '%s' is not a file. Skipping binary file comparison",
156 1 : poInputDS->GetDescription());
157 1 : return false;
158 : }
159 :
160 6 : fpRef->Seek(0, SEEK_END);
161 6 : fpInput->Seek(0, SEEK_END);
162 6 : const auto nRefSize = fpRef->Tell();
163 6 : const auto nInputSize = fpInput->Tell();
164 6 : if (nRefSize != nInputSize)
165 : {
166 4 : aosReport.push_back("Reference file has size " +
167 6 : std::to_string(nRefSize) +
168 4 : " bytes, whereas input file has size " +
169 8 : std::to_string(nInputSize) + " bytes.");
170 :
171 2 : return false;
172 : }
173 :
174 4 : constexpr size_t BUF_SIZE = 1024 * 1024;
175 8 : std::vector<GByte> abyRef(BUF_SIZE);
176 8 : std::vector<GByte> abyInput(BUF_SIZE);
177 :
178 4 : fpRef->Seek(0, SEEK_SET);
179 4 : fpInput->Seek(0, SEEK_SET);
180 :
181 0 : do
182 : {
183 4 : const size_t nRefRead = fpRef->Read(abyRef.data(), 1, BUF_SIZE);
184 4 : const size_t nInputRead = fpInput->Read(abyInput.data(), 1, BUF_SIZE);
185 :
186 4 : if (nRefRead != BUF_SIZE && fpRef->Tell() != nRefSize)
187 : {
188 0 : aosReport.push_back("Failed to fully read reference file");
189 0 : return false;
190 : }
191 :
192 4 : if (nInputRead != BUF_SIZE && fpInput->Tell() != nRefSize)
193 : {
194 0 : aosReport.push_back("Failed to fully read input file");
195 0 : return false;
196 : }
197 :
198 4 : if (abyRef != abyInput)
199 : {
200 1 : aosReport.push_back(
201 : "Reference file and input file differ at the binary level.");
202 1 : return false;
203 : }
204 3 : } while (fpRef->Tell() < nRefSize);
205 :
206 3 : return true;
207 : }
208 :
209 : /************************************************************************/
210 : /* GDALRasterCompareAlgorithm::CRSComparison() */
211 : /************************************************************************/
212 :
213 196 : void GDALRasterCompareAlgorithm::CRSComparison(
214 : std::vector<std::string> &aosReport, GDALDataset *poRefDS,
215 : GDALDataset *poInputDS)
216 : {
217 196 : const auto poRefCRS = poRefDS->GetSpatialRef();
218 196 : const auto poInputCRS = poInputDS->GetSpatialRef();
219 :
220 196 : if (poRefCRS == nullptr)
221 : {
222 177 : if (poInputCRS)
223 : {
224 1 : aosReport.push_back(
225 : "Reference dataset has no CRS, but input dataset has one.");
226 : }
227 195 : return;
228 : }
229 :
230 19 : if (poInputCRS == nullptr)
231 : {
232 1 : aosReport.push_back(
233 : "Reference dataset has a CRS, but input dataset has none.");
234 1 : return;
235 : }
236 :
237 18 : if (poRefCRS->IsSame(poInputCRS))
238 17 : return;
239 :
240 1 : const char *apszOptions[] = {"FORMAT=WKT2_2019", nullptr};
241 2 : const auto poRefWKT = poRefCRS->exportToWkt(apszOptions);
242 1 : const auto poInputWKT = poInputCRS->exportToWkt(apszOptions);
243 1 : aosReport.push_back(
244 2 : "Reference and input CRS are not equivalent. Reference one is '" +
245 2 : poRefWKT + "'. Input one is '" + poInputWKT + "'");
246 : }
247 :
248 : /************************************************************************/
249 : /* GDALRasterCompareAlgorithm::GeotransformComparison() */
250 : /************************************************************************/
251 :
252 196 : void GDALRasterCompareAlgorithm::GeoTransformComparison(
253 : std::vector<std::string> &aosReport, GDALDataset *poRefDS,
254 : GDALDataset *poInputDS)
255 : {
256 196 : GDALGeoTransform refGT;
257 196 : CPLErr eErr1 = poRefDS->GetGeoTransform(refGT);
258 196 : GDALGeoTransform inputGT;
259 196 : CPLErr eErr2 = poInputDS->GetGeoTransform(inputGT);
260 196 : if (eErr1 == CE_Failure && eErr2 == CE_Failure)
261 174 : return;
262 :
263 25 : if (eErr1 == CE_Failure && eErr2 == CE_None)
264 : {
265 1 : aosReport.push_back(
266 : "Reference dataset has no geotransform, but input one has one.");
267 1 : return;
268 : }
269 :
270 24 : if (eErr1 == CE_None && eErr2 == CE_Failure)
271 : {
272 1 : aosReport.push_back(
273 : "Reference dataset has a geotransform, but input one has none.");
274 1 : return;
275 : }
276 :
277 160 : for (int i = 0; i < 6; ++i)
278 : {
279 138 : if ((refGT[i] != 0 &&
280 275 : std::fabs(refGT[i] - inputGT[i]) > 1e-10 * std::fabs(refGT[i])) ||
281 137 : (refGT[i] == 0 && std::fabs(refGT[i] - inputGT[i]) > 1e-10))
282 : {
283 : std::string s = "Geotransform of reference and input dataset are "
284 1 : "not equivalent. Reference geotransform is (";
285 7 : for (int j = 0; j < 6; ++j)
286 : {
287 6 : if (j > 0)
288 5 : s += ',';
289 6 : s += std::to_string(refGT[j]);
290 : }
291 1 : s += "). Input geotransform is (";
292 7 : for (int j = 0; j < 6; ++j)
293 : {
294 6 : if (j > 0)
295 5 : s += ',';
296 6 : s += std::to_string(inputGT[j]);
297 : }
298 1 : s += ')';
299 1 : aosReport.push_back(std::move(s));
300 1 : return;
301 : }
302 : }
303 : }
304 :
305 : #if defined(__GNUC__) && !defined(__clang__)
306 : #pragma GCC push_options
307 : #pragma GCC optimize("O3")
308 : #endif
309 :
310 : /************************************************************************/
311 : /* Diff() */
312 : /************************************************************************/
313 :
314 2376184 : template <class T> CPL_NOSANITIZE_UNSIGNED_INT_OVERFLOW static T Diff(T a, T b)
315 : {
316 2376184 : return a - b;
317 : }
318 :
319 : /************************************************************************/
320 : /* CompareVectors() */
321 : /************************************************************************/
322 :
323 : template <class T, class Tdiff, bool bIsComplex>
324 1116 : static void CompareVectors(size_t nValCount, const T *refValues,
325 : const T *inputValues, uint64_t &countDiffPixels,
326 : Tdiff &maxDiffValue)
327 : {
328 : constexpr bool bIsFloatingPoint = std::is_floating_point_v<T>;
329 : if constexpr (bIsComplex)
330 : {
331 374 : for (size_t i = 0; i < nValCount; ++i)
332 : {
333 : if constexpr (bIsFloatingPoint)
334 : {
335 : static_assert(std::is_same_v<T, Tdiff>);
336 121 : if (std::isnan(refValues[2 * i]) &&
337 4 : std::isnan(inputValues[2 * i]) &&
338 123 : std::isnan(refValues[2 * i + 1]) &&
339 2 : std::isnan(inputValues[2 * i + 1]))
340 : {
341 2 : continue;
342 : }
343 : }
344 :
345 185 : if (refValues[2 * i] != inputValues[2 * i] ||
346 175 : refValues[2 * i + 1] != inputValues[2 * i + 1])
347 : {
348 : const Tdiff diff =
349 10 : std::hypot(static_cast<Tdiff>(refValues[2 * i]) -
350 : static_cast<Tdiff>(inputValues[2 * i]),
351 10 : static_cast<Tdiff>(refValues[2 * i + 1]) -
352 10 : static_cast<Tdiff>(inputValues[2 * i + 1]));
353 10 : ++countDiffPixels;
354 10 : if (diff > maxDiffValue)
355 10 : maxDiffValue = diff;
356 : }
357 : }
358 : }
359 : else
360 : {
361 : static_assert(sizeof(Tdiff) == sizeof(T));
362 : size_t i = 0;
363 : #ifdef USE_SSE2
364 : if constexpr (std::is_same_v<T, float>)
365 : {
366 : static_assert(std::is_same_v<T, Tdiff>);
367 :
368 : auto vMaxDiff = _mm_setzero_ps();
369 :
370 : // Mask for absolute value (clears the sign bit)
371 260 : const auto absMask = _mm_castsi128_ps(
372 : _mm_set1_epi32(std::numeric_limits<int32_t>::max()));
373 :
374 : constexpr size_t VALS_PER_REG = sizeof(vMaxDiff) / sizeof(T);
375 280 : while (i + VALS_PER_REG <= nValCount)
376 : {
377 : auto vCountDiff = _mm_setzero_si128();
378 :
379 : // We can do a maximum of std::numeric_limits<uint32_t>::max()
380 : // accumulations into vCountDiff
381 20 : const size_t nInnerLimit = [i, nValCount](size_t valsPerReg)
382 : {
383 : if constexpr (sizeof(size_t) > sizeof(uint32_t))
384 : {
385 : return std::min(
386 40 : nValCount - valsPerReg,
387 20 : i + std::numeric_limits<uint32_t>::max() *
388 20 : valsPerReg);
389 : }
390 : else
391 : {
392 : return nValCount - valsPerReg;
393 : }
394 20 : }(VALS_PER_REG);
395 :
396 70 : for (; i <= nInnerLimit; i += VALS_PER_REG)
397 : {
398 50 : const auto a = _mm_loadu_ps(refValues + i);
399 50 : const auto b = _mm_loadu_ps(inputValues + i);
400 :
401 : // Compute absolute value of difference
402 : const auto absDiff = _mm_and_ps(_mm_sub_ps(a, b), absMask);
403 :
404 : // Update vMaxDiff
405 : const auto aIsNan = _mm_cmpunord_ps(a, a);
406 : const auto bIsNan = _mm_cmpunord_ps(b, b);
407 : const auto valNotEqual = _mm_andnot_ps(
408 : _mm_or_ps(aIsNan, bIsNan), _mm_cmpneq_ps(a, b));
409 : vMaxDiff =
410 : _mm_max_ps(vMaxDiff, _mm_and_ps(absDiff, valNotEqual));
411 :
412 : // Update vCountDiff
413 : const auto nanMisMatch = _mm_xor_ps(aIsNan, bIsNan);
414 : // if nanMisMatch OR (both values not NaN and a != b)
415 : const auto maskIsDiff = _mm_or_ps(nanMisMatch, valNotEqual);
416 : const auto shiftedMaskDiff =
417 : _mm_srli_epi32(_mm_castps_si128(maskIsDiff), 31);
418 : vCountDiff = _mm_add_epi32(vCountDiff, shiftedMaskDiff);
419 : }
420 :
421 : // Horizontal add into countDiffPixels
422 : uint32_t anCountDiff[VALS_PER_REG];
423 : _mm_storeu_si128(reinterpret_cast<__m128i *>(anCountDiff),
424 : vCountDiff);
425 100 : for (size_t j = 0; j < VALS_PER_REG; ++j)
426 : {
427 80 : countDiffPixels += anCountDiff[j];
428 : }
429 : }
430 :
431 : // Horizontal max into maxDiffValue
432 : float afMaxDiffValue[VALS_PER_REG];
433 : _mm_storeu_ps(afMaxDiffValue, vMaxDiff);
434 1300 : for (size_t j = 0; j < VALS_PER_REG; ++j)
435 : {
436 1040 : CPLAssert(!std::isnan(afMaxDiffValue[j]));
437 1040 : maxDiffValue = std::max(maxDiffValue, afMaxDiffValue[j]);
438 : }
439 : }
440 : #endif
441 : if constexpr (bIsFloatingPoint)
442 : {
443 : static_assert(std::is_same_v<T, Tdiff>);
444 724 : for (; i < nValCount; ++i)
445 : {
446 362 : if (std::isnan(refValues[i]))
447 : {
448 14 : if (!std::isnan(inputValues[i]))
449 : {
450 7 : ++countDiffPixels;
451 : }
452 14 : continue;
453 : }
454 348 : else if (std::isnan(inputValues[i]))
455 : {
456 7 : ++countDiffPixels;
457 7 : continue;
458 : }
459 341 : else if (refValues[i] == inputValues[i])
460 : {
461 325 : continue;
462 : }
463 :
464 : const Tdiff diff =
465 : refValues[i] >= inputValues[i]
466 16 : ? Diff(static_cast<Tdiff>(refValues[i]),
467 : static_cast<Tdiff>(inputValues[i]))
468 8 : : Diff(static_cast<Tdiff>(inputValues[i]),
469 : static_cast<Tdiff>(refValues[i]));
470 16 : if (diff > 0)
471 : {
472 16 : ++countDiffPixels;
473 16 : if (diff > maxDiffValue)
474 12 : maxDiffValue = diff;
475 : }
476 : }
477 : }
478 : else
479 : {
480 : static_assert(std::is_unsigned_v<Tdiff>);
481 10397 : while (i < nValCount)
482 : {
483 : // Autovectorizer friendly inner loop (GCC, clang, ICX),
484 : // by making sure it increases countDiffLocal on the same size
485 : // as Tdiff.
486 :
487 : Tdiff countDiffLocal = 0;
488 9830 : const size_t innerLimit = [i, nValCount]()
489 : {
490 : if constexpr (sizeof(Tdiff) < sizeof(size_t))
491 : {
492 19422 : return std::min(nValCount - 1,
493 9711 : i + std::numeric_limits<Tdiff>::max());
494 : }
495 : else
496 : {
497 : (void)i;
498 119 : return nValCount - 1;
499 : }
500 9830 : }();
501 2385993 : for (; i <= innerLimit; ++i)
502 : {
503 2376166 : const Tdiff diff =
504 2376166 : refValues[i] >= inputValues[i]
505 2376166 : ? Diff(static_cast<Tdiff>(refValues[i]),
506 133 : static_cast<Tdiff>(inputValues[i]))
507 16 : : Diff(static_cast<Tdiff>(inputValues[i]),
508 4 : static_cast<Tdiff>(refValues[i]));
509 2376166 : countDiffLocal += (diff > 0);
510 2376166 : maxDiffValue = std::max(maxDiffValue, diff);
511 : }
512 9830 : countDiffPixels += countDiffLocal;
513 : }
514 : }
515 : }
516 1116 : }
517 :
518 : /************************************************************************/
519 : /* DatasetPixelComparison() */
520 : /************************************************************************/
521 :
522 : template <class T, class Tdiff, bool bIsComplex>
523 61 : static void DatasetPixelComparison(std::vector<std::string> &aosReport,
524 : GDALDataset *poRefDS, GDALDataset *poInputDS,
525 : GDALDataType eReqDT,
526 : GDALProgressFunc pfnProgress,
527 : void *pProgressData)
528 : {
529 122 : std::vector<T> refValues;
530 122 : std::vector<T> inputValues;
531 :
532 61 : CPLAssert(GDALDataTypeIsComplex(eReqDT) == bIsComplex);
533 :
534 61 : const uint64_t nTotalPixels =
535 61 : static_cast<uint64_t>(poRefDS->GetRasterXSize()) *
536 61 : poRefDS->GetRasterYSize();
537 : uint64_t nIterPixels = 0;
538 :
539 : constexpr int nValPerPixel = bIsComplex ? 2 : 1;
540 61 : const int nBands = poRefDS->GetRasterCount();
541 :
542 122 : std::vector<Tdiff> maxDiffValue(nBands, 0);
543 122 : std::vector<uint64_t> countDiffPixels(nBands, 0);
544 :
545 : size_t nMaxSize = 0;
546 61 : const GIntBig nUsableRAM = CPLGetUsablePhysicalRAM() / 10;
547 61 : if (nUsableRAM > 0)
548 : nMaxSize = static_cast<size_t>(nUsableRAM);
549 :
550 121 : for (const auto &window : GDALRasterBand::WindowIteratorWrapper(
551 61 : *(poRefDS->GetRasterBand(1)), *(poInputDS->GetRasterBand(1)),
552 : nMaxSize))
553 : {
554 61 : const size_t nValCount =
555 61 : static_cast<size_t>(window.nXSize) * window.nYSize;
556 61 : const size_t nArraySize = nValCount * nValPerPixel * nBands;
557 : try
558 : {
559 61 : if (refValues.size() < nArraySize)
560 : {
561 61 : refValues.resize(nArraySize);
562 61 : inputValues.resize(nArraySize);
563 : }
564 : }
565 0 : catch (const std::exception &)
566 : {
567 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
568 : "Out of memory allocating temporary arrays");
569 0 : aosReport.push_back("Out of memory allocating temporary arrays");
570 : return;
571 : }
572 :
573 : if (poRefDS->RasterIO(GF_Read, window.nXOff, window.nYOff,
574 : window.nXSize, window.nYSize, refValues.data(),
575 : window.nXSize, window.nYSize, eReqDT, nBands,
576 122 : nullptr, 0, 0, 0, nullptr) == CE_None &&
577 : poInputDS->RasterIO(
578 : GF_Read, window.nXOff, window.nYOff, window.nXSize,
579 : window.nYSize, inputValues.data(), window.nXSize, window.nYSize,
580 61 : eReqDT, nBands, nullptr, 0, 0, 0, nullptr) == CE_None)
581 : {
582 1037 : for (int i = 0; i < nBands; ++i)
583 : {
584 976 : CompareVectors<T, Tdiff, bIsComplex>(
585 976 : nValCount, refValues.data() + i * nValCount * nValPerPixel,
586 976 : inputValues.data() + i * nValCount * nValPerPixel,
587 976 : countDiffPixels[i], maxDiffValue[i]);
588 : }
589 : }
590 : else
591 : {
592 0 : aosReport.push_back("I/O error when comparing pixel values");
593 : }
594 :
595 61 : if (pfnProgress)
596 : {
597 2 : nIterPixels += nValCount;
598 2 : if (!pfnProgress(static_cast<double>(nIterPixels) /
599 : static_cast<double>(nTotalPixels),
600 : "", pProgressData))
601 : {
602 1 : CPLError(CE_Failure, CPLE_UserInterrupt, "Interrupted by user");
603 : break;
604 : }
605 : }
606 : }
607 1037 : for (int i = 0; i < nBands; ++i)
608 : {
609 976 : if (countDiffPixels[i])
610 : {
611 29 : aosReport.push_back(
612 : "Band " + std::to_string(i + 1) +
613 29 : ": pixels differing: " + std::to_string(countDiffPixels[i]));
614 29 : aosReport.push_back("Band " + std::to_string(i + 1) +
615 : ": maximum pixel value difference: " +
616 29 : std::to_string(maxDiffValue[i]));
617 : }
618 : }
619 : }
620 :
621 : /************************************************************************/
622 : /* GDALRasterCompareAlgorithm::DatasetComparison() */
623 : /************************************************************************/
624 :
625 205 : void GDALRasterCompareAlgorithm::DatasetComparison(
626 : std::vector<std::string> &aosReport, GDALDataset *poRefDS,
627 : GDALDataset *poInputDS, GDALProgressFunc pfnProgress, void *pProgressData)
628 : {
629 205 : if (!m_skipCRS)
630 : {
631 196 : CRSComparison(aosReport, poRefDS, poInputDS);
632 : }
633 :
634 205 : if (!m_skipGeotransform)
635 : {
636 196 : GeoTransformComparison(aosReport, poRefDS, poInputDS);
637 : }
638 :
639 : bool ret = true;
640 205 : if (poRefDS->GetRasterCount() != poInputDS->GetRasterCount())
641 : {
642 2 : aosReport.push_back("Reference dataset has " +
643 3 : std::to_string(poRefDS->GetRasterCount()) +
644 2 : " band(s), but input dataset has " +
645 2 : std::to_string(poInputDS->GetRasterCount()));
646 : ret = false;
647 : }
648 :
649 205 : if (poRefDS->GetRasterXSize() != poInputDS->GetRasterXSize())
650 : {
651 2 : aosReport.push_back("Reference dataset width is " +
652 3 : std::to_string(poRefDS->GetRasterXSize()) +
653 2 : ", but input dataset width is " +
654 2 : std::to_string(poInputDS->GetRasterXSize()));
655 : ret = false;
656 : }
657 :
658 205 : if (poRefDS->GetRasterYSize() != poInputDS->GetRasterYSize())
659 : {
660 2 : aosReport.push_back("Reference dataset height is " +
661 3 : std::to_string(poRefDS->GetRasterYSize()) +
662 2 : ", but input dataset height is " +
663 2 : std::to_string(poInputDS->GetRasterYSize()));
664 : ret = false;
665 : }
666 :
667 205 : if (!m_skipMetadata)
668 : {
669 196 : MetadataComparison(aosReport, "(dataset default metadata domain)",
670 196 : poRefDS->GetMetadata(), poInputDS->GetMetadata());
671 : }
672 :
673 205 : if (!m_skipRPC)
674 : {
675 196 : MetadataComparison(aosReport, "RPC", poRefDS->GetMetadata("RPC"),
676 196 : poInputDS->GetMetadata("RPC"));
677 : }
678 :
679 205 : if (!m_skipGeolocation)
680 : {
681 196 : MetadataComparison(aosReport, "GEOLOCATION",
682 196 : poRefDS->GetMetadata("GEOLOCATION"),
683 196 : poInputDS->GetMetadata("GEOLOCATION"));
684 : }
685 :
686 205 : if (!ret)
687 : return;
688 :
689 202 : const int nBands = poRefDS->GetRasterCount();
690 :
691 : bool doBandBasedPixelComparison = true;
692 : // Do not do band-by-band pixel difference if there are too many interleaved
693 : // bands as this could be extremely slow
694 202 : if (nBands > 10)
695 : {
696 : const char *pszRefInterleave =
697 61 : poRefDS->GetMetadataItem("INTERLEAVE", "IMAGE_STRUCTURE");
698 : const char *pszInputInterleave =
699 61 : poInputDS->GetMetadataItem("INTERLEAVE", "IMAGE_STRUCTURE");
700 61 : if ((pszRefInterleave && EQUAL(pszRefInterleave, "PIXEL")) ||
701 0 : (pszInputInterleave && EQUAL(pszInputInterleave, "PIXEL")))
702 : {
703 : doBandBasedPixelComparison = false;
704 : }
705 : }
706 :
707 1315 : for (int i = 0; i < nBands; ++i)
708 : {
709 1113 : void *pScaledProgress = GDALCreateScaledProgress(
710 1113 : static_cast<double>(i) / nBands,
711 1113 : static_cast<double>(i + 1) / nBands, pfnProgress, pProgressData);
712 3301 : BandComparison(
713 2226 : aosReport, std::to_string(i + 1), doBandBasedPixelComparison,
714 : poRefDS->GetRasterBand(i + 1), poInputDS->GetRasterBand(i + 1),
715 : pScaledProgress ? GDALScaledProgress : nullptr, pScaledProgress);
716 1113 : GDALDestroyScaledProgress(pScaledProgress);
717 : }
718 :
719 202 : if (!doBandBasedPixelComparison)
720 : {
721 : const auto eReqDT =
722 61 : GDALDataTypeUnion(poRefDS->GetRasterBand(1)->GetRasterDataType(),
723 : poInputDS->GetRasterBand(1)->GetRasterDataType());
724 61 : switch (eReqDT)
725 : {
726 5 : case GDT_UInt8:
727 5 : DatasetPixelComparison<uint8_t, uint8_t, false>(
728 : aosReport, poRefDS, poInputDS, eReqDT, pfnProgress,
729 : pProgressData);
730 5 : break;
731 4 : case GDT_Int8:
732 4 : DatasetPixelComparison<int8_t, uint8_t, false>(
733 : aosReport, poRefDS, poInputDS, eReqDT, pfnProgress,
734 : pProgressData);
735 4 : break;
736 3 : case GDT_UInt16:
737 3 : DatasetPixelComparison<uint16_t, uint16_t, false>(
738 : aosReport, poRefDS, poInputDS, eReqDT, pfnProgress,
739 : pProgressData);
740 3 : break;
741 4 : case GDT_Int16:
742 4 : DatasetPixelComparison<int16_t, uint16_t, false>(
743 : aosReport, poRefDS, poInputDS, eReqDT, pfnProgress,
744 : pProgressData);
745 4 : break;
746 3 : case GDT_UInt32:
747 3 : DatasetPixelComparison<uint32_t, uint32_t, false>(
748 : aosReport, poRefDS, poInputDS, eReqDT, pfnProgress,
749 : pProgressData);
750 3 : break;
751 4 : case GDT_Int32:
752 4 : DatasetPixelComparison<int32_t, uint32_t, false>(
753 : aosReport, poRefDS, poInputDS, eReqDT, pfnProgress,
754 : pProgressData);
755 4 : break;
756 3 : case GDT_UInt64:
757 3 : DatasetPixelComparison<uint64_t, uint64_t, false>(
758 : aosReport, poRefDS, poInputDS, eReqDT, pfnProgress,
759 : pProgressData);
760 3 : break;
761 4 : case GDT_Int64:
762 4 : DatasetPixelComparison<int64_t, uint64_t, false>(
763 : aosReport, poRefDS, poInputDS, eReqDT, pfnProgress,
764 : pProgressData);
765 4 : break;
766 14 : case GDT_Float16:
767 : case GDT_Float32:
768 14 : DatasetPixelComparison<float, float, false>(
769 : aosReport, poRefDS, poInputDS, GDT_Float32, pfnProgress,
770 : pProgressData);
771 14 : break;
772 6 : case GDT_Float64:
773 6 : DatasetPixelComparison<double, double, false>(
774 : aosReport, poRefDS, poInputDS, eReqDT, pfnProgress,
775 : pProgressData);
776 6 : break;
777 2 : case GDT_CInt16:
778 2 : DatasetPixelComparison<int16_t, float, true>(
779 : aosReport, poRefDS, poInputDS, eReqDT, pfnProgress,
780 : pProgressData);
781 2 : break;
782 2 : case GDT_CInt32:
783 2 : DatasetPixelComparison<int32_t, double, true>(
784 : aosReport, poRefDS, poInputDS, eReqDT, pfnProgress,
785 : pProgressData);
786 2 : break;
787 5 : case GDT_CFloat16:
788 : case GDT_CFloat32:
789 5 : DatasetPixelComparison<float, float, true>(
790 : aosReport, poRefDS, poInputDS, GDT_CFloat32, pfnProgress,
791 : pProgressData);
792 5 : break;
793 2 : case GDT_CFloat64:
794 2 : DatasetPixelComparison<double, double, true>(
795 : aosReport, poRefDS, poInputDS, eReqDT, pfnProgress,
796 : pProgressData);
797 2 : break;
798 : case GDT_Unknown:
799 : case GDT_TypeCount:
800 : break;
801 : }
802 : }
803 : }
804 :
805 : /************************************************************************/
806 : /* ComparePixels() */
807 : /************************************************************************/
808 :
809 : template <class T, class Tdiff, bool bIsComplex>
810 140 : static void ComparePixels(std::vector<std::string> &aosReport,
811 : const std::string &bandId, GDALRasterBand *poRefBand,
812 : GDALRasterBand *poInputBand, GDALDataType eReqDT,
813 : GDALProgressFunc pfnProgress, void *pProgressData)
814 : {
815 280 : std::vector<T> refValues;
816 280 : std::vector<T> inputValues;
817 140 : Tdiff maxDiffValue = 0;
818 140 : uint64_t countDiffPixels = 0;
819 :
820 140 : CPLAssert(GDALDataTypeIsComplex(eReqDT) == bIsComplex);
821 140 : const uint64_t nTotalPixels =
822 140 : static_cast<uint64_t>(poRefBand->GetXSize()) * poRefBand->GetYSize();
823 : uint64_t nIterPixels = 0;
824 :
825 : constexpr int nValPerPixel = bIsComplex ? 2 : 1;
826 :
827 : size_t nMaxSize = 0;
828 140 : const GIntBig nUsableRAM = CPLGetUsablePhysicalRAM() / 10;
829 140 : if (nUsableRAM > 0)
830 : nMaxSize = static_cast<size_t>(nUsableRAM);
831 :
832 279 : for (const auto &window : GDALRasterBand::WindowIteratorWrapper(
833 : *poRefBand, *poInputBand, nMaxSize))
834 : {
835 140 : const size_t nValCount =
836 140 : static_cast<size_t>(window.nXSize) * window.nYSize;
837 11 : const size_t nArraySize = nValCount * nValPerPixel;
838 : try
839 : {
840 140 : if (refValues.size() < nArraySize)
841 : {
842 140 : refValues.resize(nArraySize);
843 140 : inputValues.resize(nArraySize);
844 : }
845 : }
846 0 : catch (const std::exception &)
847 : {
848 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
849 : "Out of memory allocating temporary arrays");
850 0 : aosReport.push_back("Out of memory allocating temporary arrays");
851 : return;
852 : }
853 :
854 : if (poRefBand->RasterIO(GF_Read, window.nXOff, window.nYOff,
855 : window.nXSize, window.nYSize, refValues.data(),
856 : window.nXSize, window.nYSize, eReqDT, 0, 0,
857 280 : nullptr) == CE_None &&
858 : poInputBand->RasterIO(
859 : GF_Read, window.nXOff, window.nYOff, window.nXSize,
860 : window.nYSize, inputValues.data(), window.nXSize, window.nYSize,
861 140 : eReqDT, 0, 0, nullptr) == CE_None)
862 : {
863 140 : CompareVectors<T, Tdiff, bIsComplex>(nValCount, refValues.data(),
864 : inputValues.data(),
865 : countDiffPixels, maxDiffValue);
866 : }
867 : else
868 : {
869 0 : aosReport.push_back("I/O error when comparing pixel values");
870 : }
871 :
872 140 : if (pfnProgress)
873 : {
874 6 : nIterPixels += nValCount;
875 6 : if (!pfnProgress(static_cast<double>(nIterPixels) /
876 : static_cast<double>(nTotalPixels),
877 : "", pProgressData))
878 : {
879 1 : CPLError(CE_Failure, CPLE_UserInterrupt, "Interrupted by user");
880 : break;
881 : }
882 : }
883 : }
884 140 : if (countDiffPixels)
885 : {
886 46 : aosReport.push_back(
887 : bandId + ": pixels differing: " + std::to_string(countDiffPixels));
888 :
889 92 : std::string reportMessage(bandId);
890 46 : reportMessage += ": maximum pixel value difference: ";
891 : if constexpr (std::is_floating_point_v<T>)
892 : {
893 27 : if (std::isinf(maxDiffValue))
894 2 : reportMessage += "inf";
895 25 : else if (std::isnan(maxDiffValue))
896 0 : reportMessage += "nan";
897 : else
898 25 : reportMessage += std::to_string(maxDiffValue);
899 : }
900 : else
901 : {
902 19 : reportMessage += std::to_string(maxDiffValue);
903 : }
904 46 : aosReport.push_back(std::move(reportMessage));
905 : }
906 : }
907 :
908 : /************************************************************************/
909 : /* ComparePixels() */
910 : /************************************************************************/
911 :
912 140 : static void ComparePixels(std::vector<std::string> &aosReport,
913 : const std::string &bandId, GDALRasterBand *poRefBand,
914 : GDALRasterBand *poInputBand,
915 : GDALProgressFunc pfnProgress, void *pProgressData)
916 : {
917 140 : const auto eReqDT = GDALDataTypeUnion(poRefBand->GetRasterDataType(),
918 : poInputBand->GetRasterDataType());
919 140 : switch (eReqDT)
920 : {
921 60 : case GDT_UInt8:
922 60 : ComparePixels<uint8_t, uint8_t, false>(aosReport, bandId, poRefBand,
923 : poInputBand, eReqDT,
924 : pfnProgress, pProgressData);
925 60 : break;
926 4 : case GDT_Int8:
927 4 : ComparePixels<int8_t, uint8_t, false>(aosReport, bandId, poRefBand,
928 : poInputBand, eReqDT,
929 : pfnProgress, pProgressData);
930 4 : break;
931 4 : case GDT_UInt16:
932 4 : ComparePixels<uint16_t, uint16_t, false>(
933 : aosReport, bandId, poRefBand, poInputBand, eReqDT, pfnProgress,
934 : pProgressData);
935 4 : break;
936 5 : case GDT_Int16:
937 5 : ComparePixels<int16_t, uint16_t, false>(
938 : aosReport, bandId, poRefBand, poInputBand, eReqDT, pfnProgress,
939 : pProgressData);
940 5 : break;
941 3 : case GDT_UInt32:
942 3 : ComparePixels<uint32_t, uint32_t, false>(
943 : aosReport, bandId, poRefBand, poInputBand, eReqDT, pfnProgress,
944 : pProgressData);
945 3 : break;
946 4 : case GDT_Int32:
947 4 : ComparePixels<int32_t, uint32_t, false>(
948 : aosReport, bandId, poRefBand, poInputBand, eReqDT, pfnProgress,
949 : pProgressData);
950 4 : break;
951 3 : case GDT_UInt64:
952 3 : ComparePixels<uint64_t, uint64_t, false>(
953 : aosReport, bandId, poRefBand, poInputBand, eReqDT, pfnProgress,
954 : pProgressData);
955 3 : break;
956 4 : case GDT_Int64:
957 4 : ComparePixels<int64_t, uint64_t, false>(
958 : aosReport, bandId, poRefBand, poInputBand, eReqDT, pfnProgress,
959 : pProgressData);
960 4 : break;
961 36 : case GDT_Float16:
962 : case GDT_Float32:
963 36 : ComparePixels<float, float, false>(aosReport, bandId, poRefBand,
964 : poInputBand, GDT_Float32,
965 : pfnProgress, pProgressData);
966 36 : break;
967 6 : case GDT_Float64:
968 6 : ComparePixels<double, double, false>(aosReport, bandId, poRefBand,
969 : poInputBand, eReqDT,
970 : pfnProgress, pProgressData);
971 6 : break;
972 2 : case GDT_CInt16:
973 2 : ComparePixels<int16_t, float, true>(aosReport, bandId, poRefBand,
974 : poInputBand, eReqDT,
975 : pfnProgress, pProgressData);
976 2 : break;
977 2 : case GDT_CInt32:
978 2 : ComparePixels<int32_t, double, true>(aosReport, bandId, poRefBand,
979 : poInputBand, eReqDT,
980 : pfnProgress, pProgressData);
981 2 : break;
982 5 : case GDT_CFloat16:
983 : case GDT_CFloat32:
984 5 : ComparePixels<float, float, true>(aosReport, bandId, poRefBand,
985 : poInputBand, GDT_CFloat32,
986 : pfnProgress, pProgressData);
987 5 : break;
988 2 : case GDT_CFloat64:
989 2 : ComparePixels<double, double, true>(aosReport, bandId, poRefBand,
990 : poInputBand, eReqDT,
991 : pfnProgress, pProgressData);
992 2 : break;
993 : case GDT_Unknown:
994 : case GDT_TypeCount:
995 : break;
996 : }
997 140 : }
998 :
999 : #if defined(__GNUC__) && !defined(__clang__)
1000 : #pragma GCC pop_options
1001 : #endif
1002 :
1003 : /************************************************************************/
1004 : /* GDALRasterCompareAlgorithm::BandComparison() */
1005 : /************************************************************************/
1006 :
1007 1117 : void GDALRasterCompareAlgorithm::BandComparison(
1008 : std::vector<std::string> &aosReport, const std::string &bandId,
1009 : bool doBandBasedPixelComparison, GDALRasterBand *poRefBand,
1010 : GDALRasterBand *poInputBand, GDALProgressFunc pfnProgress,
1011 : void *pProgressData)
1012 : {
1013 1117 : bool ret = true;
1014 :
1015 1117 : if (poRefBand->GetXSize() != poInputBand->GetXSize())
1016 : {
1017 2 : aosReport.push_back("Reference band width is " +
1018 3 : std::to_string(poRefBand->GetXSize()) +
1019 2 : ", but input band width is " +
1020 2 : std::to_string(poInputBand->GetXSize()));
1021 1 : ret = false;
1022 : }
1023 :
1024 1117 : if (poRefBand->GetYSize() != poInputBand->GetYSize())
1025 : {
1026 2 : aosReport.push_back("Reference band height is " +
1027 3 : std::to_string(poRefBand->GetYSize()) +
1028 2 : ", but input band height is " +
1029 2 : std::to_string(poInputBand->GetYSize()));
1030 1 : ret = false;
1031 : }
1032 :
1033 1117 : if (strcmp(poRefBand->GetDescription(), poInputBand->GetDescription()) != 0)
1034 : {
1035 3 : aosReport.push_back("Reference band " + bandId + " has description " +
1036 4 : std::string(poRefBand->GetDescription()) +
1037 2 : ", but input band has description " +
1038 2 : std::string(poInputBand->GetDescription()));
1039 : }
1040 :
1041 1117 : if (poRefBand->GetRasterDataType() != poInputBand->GetRasterDataType())
1042 : {
1043 2 : aosReport.push_back(
1044 4 : "Reference band " + bandId + " has data type " +
1045 8 : std::string(GDALGetDataTypeName(poRefBand->GetRasterDataType())) +
1046 4 : ", but input band has data type " +
1047 4 : std::string(GDALGetDataTypeName(poInputBand->GetRasterDataType())));
1048 : }
1049 :
1050 1117 : int bRefHasNoData = false;
1051 1117 : const double dfRefNoData = poRefBand->GetNoDataValue(&bRefHasNoData);
1052 1117 : int bInputHasNoData = false;
1053 1117 : const double dfInputNoData = poInputBand->GetNoDataValue(&bInputHasNoData);
1054 1117 : if (!bRefHasNoData && !bInputHasNoData)
1055 : {
1056 : // ok
1057 : }
1058 6 : else if (bRefHasNoData && !bInputHasNoData)
1059 : {
1060 3 : aosReport.push_back("Reference band " + bandId + " has nodata value " +
1061 4 : std::to_string(dfRefNoData) +
1062 : ", but input band has none.");
1063 : }
1064 5 : else if (!bRefHasNoData && bInputHasNoData)
1065 : {
1066 3 : aosReport.push_back("Reference band " + bandId +
1067 2 : " has no nodata value, " +
1068 2 : "but input band has no data value " +
1069 4 : std::to_string(dfInputNoData) + ".");
1070 : }
1071 4 : else if ((std::isnan(dfRefNoData) && std::isnan(dfInputNoData)) ||
1072 : dfRefNoData == dfInputNoData)
1073 : {
1074 : // ok
1075 : }
1076 : else
1077 : {
1078 6 : aosReport.push_back("Reference band " + bandId + " has nodata value " +
1079 8 : std::to_string(dfRefNoData) +
1080 4 : ", but input band has no data value " +
1081 8 : std::to_string(dfInputNoData) + ".");
1082 : }
1083 :
1084 1117 : if (poRefBand->GetColorInterpretation() !=
1085 1117 : poInputBand->GetColorInterpretation())
1086 : {
1087 3 : aosReport.push_back("Reference band " + bandId +
1088 2 : " has color interpretation " +
1089 3 : std::string(GDALGetColorInterpretationName(
1090 3 : poRefBand->GetColorInterpretation())) +
1091 2 : ", but input band has color interpretation " +
1092 3 : std::string(GDALGetColorInterpretationName(
1093 1 : poInputBand->GetColorInterpretation())));
1094 : }
1095 :
1096 1117 : if (!ret)
1097 1 : return;
1098 :
1099 : const uint64_t nBasePixels =
1100 1116 : static_cast<uint64_t>(poRefBand->GetXSize()) * poRefBand->GetYSize();
1101 1116 : uint64_t nTotalPixels = nBasePixels;
1102 1116 : const int nOvrCount = poRefBand->GetOverviewCount();
1103 1118 : for (int i = 0; i < nOvrCount; ++i)
1104 : {
1105 2 : auto poOvrBand = poRefBand->GetOverview(i);
1106 : const uint64_t nOvrPixels =
1107 2 : static_cast<uint64_t>(poOvrBand->GetXSize()) *
1108 2 : poOvrBand->GetYSize();
1109 2 : nTotalPixels += nOvrPixels;
1110 : }
1111 :
1112 1116 : if (doBandBasedPixelComparison)
1113 : {
1114 : void *pScaledProgress =
1115 280 : GDALCreateScaledProgress(0.0,
1116 140 : static_cast<double>(nBasePixels) /
1117 140 : static_cast<double>(nTotalPixels),
1118 : pfnProgress, pProgressData);
1119 140 : ComparePixels(aosReport, bandId, poRefBand, poInputBand,
1120 : pScaledProgress ? GDALScaledProgress : nullptr,
1121 : pScaledProgress);
1122 140 : GDALDestroyScaledProgress(pScaledProgress);
1123 : }
1124 :
1125 1116 : if (!m_skipOverview)
1126 : {
1127 1108 : if (nOvrCount != poInputBand->GetOverviewCount())
1128 : {
1129 1 : aosReport.push_back(
1130 2 : "Reference band " + bandId + " has " +
1131 4 : std::to_string(nOvrCount) +
1132 2 : " overview band(s), but input band has " +
1133 2 : std::to_string(poInputBand->GetOverviewCount()));
1134 : }
1135 : else
1136 : {
1137 1107 : uint64_t nIterPixels = nBasePixels;
1138 :
1139 1109 : for (int i = 0; i < nOvrCount; ++i)
1140 : {
1141 2 : GDALRasterBand *poOvrBand = poRefBand->GetOverview(i);
1142 : const uint64_t nOvrPixels =
1143 2 : static_cast<uint64_t>(poOvrBand->GetXSize()) *
1144 2 : poOvrBand->GetYSize();
1145 4 : void *pScaledProgress = GDALCreateScaledProgress(
1146 2 : static_cast<double>(nIterPixels) /
1147 2 : static_cast<double>(nTotalPixels),
1148 2 : static_cast<double>(nIterPixels + nOvrPixels) /
1149 2 : static_cast<double>(nTotalPixels),
1150 : pfnProgress, pProgressData);
1151 4 : BandComparison(aosReport, "overview of band " + bandId,
1152 : doBandBasedPixelComparison, poOvrBand,
1153 2 : poInputBand->GetOverview(i),
1154 : pScaledProgress ? GDALScaledProgress : nullptr,
1155 : pScaledProgress);
1156 2 : GDALDestroyScaledProgress(pScaledProgress);
1157 2 : nIterPixels += nOvrPixels;
1158 : }
1159 : }
1160 : }
1161 :
1162 1116 : if (poRefBand->GetMaskFlags() != poInputBand->GetMaskFlags())
1163 : {
1164 6 : aosReport.push_back("Reference band " + bandId + " has mask flags = " +
1165 8 : std::to_string(poRefBand->GetMaskFlags()) +
1166 4 : " , but input band has mask flags = " +
1167 4 : std::to_string(poInputBand->GetMaskFlags()));
1168 : }
1169 1114 : else if (poRefBand->GetMaskFlags() == GMF_PER_DATASET)
1170 : {
1171 2 : BandComparison(aosReport, "mask of band " + bandId, true,
1172 2 : poRefBand->GetMaskBand(), poInputBand->GetMaskBand(),
1173 : nullptr, nullptr);
1174 : }
1175 :
1176 1116 : if (!m_skipMetadata)
1177 : {
1178 1108 : MetadataComparison(aosReport, "(band default metadata domain)",
1179 1108 : poRefBand->GetMetadata(),
1180 1108 : poInputBand->GetMetadata());
1181 : }
1182 : }
1183 :
1184 : /************************************************************************/
1185 : /* GDALRasterCompareAlgorithm::MetadataComparison() */
1186 : /************************************************************************/
1187 :
1188 1696 : void GDALRasterCompareAlgorithm::MetadataComparison(
1189 : std::vector<std::string> &aosReport, const std::string &metadataDomain,
1190 : CSLConstList aosRef, CSLConstList aosInput)
1191 : {
1192 3392 : std::map<std::string, std::string> oMapRef;
1193 3392 : std::map<std::string, std::string> oMapInput;
1194 :
1195 1696 : std::array<const char *, 3> ignoredKeys = {
1196 : "backend", // from gdalcompare.py. Not sure why
1197 : "ERR_BIAS", // RPC optional key
1198 : "ERR_RAND", // RPC optional key
1199 : };
1200 :
1201 1774 : for (const auto &[key, value] : cpl::IterateNameValue(aosRef))
1202 : {
1203 78 : const char *pszKey = key;
1204 234 : const auto eq = [pszKey](const char *s)
1205 234 : { return strcmp(pszKey, s) == 0; };
1206 78 : auto it = std::find_if(ignoredKeys.begin(), ignoredKeys.end(), eq);
1207 78 : if (it == ignoredKeys.end())
1208 : {
1209 77 : oMapRef[key] = value;
1210 : }
1211 : }
1212 :
1213 1774 : for (const auto &[key, value] : cpl::IterateNameValue(aosInput))
1214 : {
1215 78 : const char *pszKey = key;
1216 233 : const auto eq = [pszKey](const char *s)
1217 233 : { return strcmp(pszKey, s) == 0; };
1218 78 : auto it = std::find_if(ignoredKeys.begin(), ignoredKeys.end(), eq);
1219 78 : if (it == ignoredKeys.end())
1220 : {
1221 77 : oMapInput[key] = value;
1222 : }
1223 : }
1224 :
1225 4 : const auto strip = [](const std::string &s)
1226 : {
1227 4 : const auto posBegin = s.find_first_not_of(' ');
1228 4 : if (posBegin == std::string::npos)
1229 0 : return std::string();
1230 4 : const auto posEnd = s.find_last_not_of(' ');
1231 4 : return s.substr(posBegin, posEnd - posBegin + 1);
1232 : };
1233 :
1234 1773 : for (const auto &sKeyValuePair : oMapRef)
1235 : {
1236 77 : const auto oIter = oMapInput.find(sKeyValuePair.first);
1237 77 : if (oIter == oMapInput.end())
1238 : {
1239 6 : aosReport.push_back("Reference metadata " + metadataDomain +
1240 6 : " contains key '" + sKeyValuePair.first +
1241 : "' but input metadata does not.");
1242 : }
1243 : else
1244 : {
1245 : // this will always have the current date set
1246 75 : if (sKeyValuePair.first == "NITF_FDT")
1247 2 : continue;
1248 :
1249 146 : std::string ref = oIter->second;
1250 146 : std::string input = sKeyValuePair.second;
1251 73 : if (metadataDomain == "RPC")
1252 : {
1253 : // _RPC.TXT files and in-file have a difference
1254 : // in white space that is not otherwise meaningful.
1255 2 : ref = strip(ref);
1256 2 : input = strip(input);
1257 : }
1258 73 : if (ref != input)
1259 : {
1260 2 : aosReport.push_back(
1261 4 : "Reference metadata " + metadataDomain + " has value '" +
1262 6 : ref + "' for key '" + sKeyValuePair.first +
1263 4 : "' but input metadata has value '" + input + "'.");
1264 : }
1265 : }
1266 : }
1267 :
1268 1773 : for (const auto &sKeyValuePair : oMapInput)
1269 : {
1270 77 : if (!cpl::contains(oMapRef, sKeyValuePair.first))
1271 : {
1272 6 : aosReport.push_back("Input metadata " + metadataDomain +
1273 6 : " contains key '" + sKeyValuePair.first +
1274 : "' but reference metadata does not.");
1275 : }
1276 : }
1277 1696 : }
1278 :
1279 : /************************************************************************/
1280 : /* GDALRasterCompareAlgorithm::RunStep() */
1281 : /************************************************************************/
1282 :
1283 205 : bool GDALRasterCompareAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
1284 : {
1285 205 : auto poRefDS = m_referenceDataset.GetDatasetRef();
1286 205 : CPLAssert(poRefDS);
1287 :
1288 205 : CPLAssert(m_inputDataset.size() == 1);
1289 205 : auto poInputDS = m_inputDataset[0].GetDatasetRef();
1290 205 : CPLAssert(poInputDS);
1291 :
1292 205 : if (m_skipAllOptional)
1293 : {
1294 8 : m_skipBinary = true;
1295 8 : m_skipCRS = true;
1296 8 : m_skipGeotransform = true;
1297 8 : m_skipOverview = true;
1298 8 : m_skipMetadata = true;
1299 8 : m_skipRPC = true;
1300 8 : m_skipGeolocation = true;
1301 8 : m_skipSubdataset = true;
1302 : }
1303 :
1304 410 : std::vector<std::string> aosReport;
1305 :
1306 205 : if (!m_skipBinary)
1307 : {
1308 12 : if (BinaryComparison(aosReport, poRefDS, poInputDS))
1309 : {
1310 3 : return true;
1311 : }
1312 : }
1313 :
1314 : CSLConstList papszSubDSRef =
1315 202 : m_skipSubdataset ? nullptr : poRefDS->GetMetadata("SUBDATASETS");
1316 202 : const int nCountRef = CSLCount(papszSubDSRef) / 2;
1317 : CSLConstList papszSubDSInput =
1318 202 : m_skipSubdataset ? nullptr : poInputDS->GetMetadata("SUBDATASETS");
1319 202 : const int nCountInput = CSLCount(papszSubDSInput) / 2;
1320 :
1321 202 : if (!m_skipSubdataset)
1322 : {
1323 193 : if (nCountRef != nCountInput)
1324 : {
1325 2 : aosReport.push_back("Reference dataset has " +
1326 3 : std::to_string(nCountRef) +
1327 2 : " subdataset(s) whereas input dataset has " +
1328 4 : std::to_string(nCountInput) + " one(s).");
1329 1 : m_skipSubdataset = true;
1330 : }
1331 : }
1332 :
1333 : // Compute total number of pixels, including in subdatasets
1334 : const uint64_t nBasePixels =
1335 202 : static_cast<uint64_t>(poRefDS->GetRasterXSize()) *
1336 202 : poRefDS->GetRasterYSize() * poRefDS->GetRasterCount();
1337 202 : uint64_t nTotalPixels = nBasePixels;
1338 202 : if (ctxt.m_pfnProgress && !m_skipSubdataset)
1339 : {
1340 7 : for (int i = 0; i < nCountRef; ++i)
1341 : {
1342 1 : const char *pszRef = CSLFetchNameValue(
1343 : papszSubDSRef, CPLSPrintf("SUBDATASET_%d_NAME", i + 1));
1344 1 : const char *pszInput = CSLFetchNameValue(
1345 : papszSubDSInput, CPLSPrintf("SUBDATASET_%d_NAME", i + 1));
1346 1 : if (pszRef && pszInput)
1347 : {
1348 : auto poSubRef = std::unique_ptr<GDALDataset>(
1349 2 : GDALDataset::Open(pszRef, GDAL_OF_RASTER));
1350 : auto poSubInput = std::unique_ptr<GDALDataset>(
1351 2 : GDALDataset::Open(pszInput, GDAL_OF_RASTER));
1352 1 : if (poSubRef && poSubInput)
1353 : {
1354 : const uint64_t nSubDSPixels =
1355 1 : static_cast<uint64_t>(poSubRef->GetRasterXSize()) *
1356 1 : poSubRef->GetRasterYSize() * poSubRef->GetRasterCount();
1357 1 : nTotalPixels += nSubDSPixels;
1358 : }
1359 : }
1360 : }
1361 : }
1362 :
1363 : {
1364 : void *pScaledProgress =
1365 404 : GDALCreateScaledProgress(0.0,
1366 202 : static_cast<double>(nBasePixels) /
1367 202 : static_cast<double>(nTotalPixels),
1368 : ctxt.m_pfnProgress, ctxt.m_pProgressData);
1369 202 : DatasetComparison(aosReport, poRefDS, poInputDS,
1370 : pScaledProgress ? GDALScaledProgress : nullptr,
1371 : pScaledProgress);
1372 202 : GDALDestroyScaledProgress(pScaledProgress);
1373 : }
1374 :
1375 202 : if (!m_skipSubdataset)
1376 : {
1377 192 : uint64_t nIterPixels = nBasePixels;
1378 195 : for (int i = 0; i < nCountRef; ++i)
1379 : {
1380 3 : const char *pszRef = CSLFetchNameValue(
1381 : papszSubDSRef, CPLSPrintf("SUBDATASET_%d_NAME", i + 1));
1382 3 : const char *pszInput = CSLFetchNameValue(
1383 : papszSubDSInput, CPLSPrintf("SUBDATASET_%d_NAME", i + 1));
1384 3 : if (pszRef && pszInput)
1385 : {
1386 : auto poSubRef = std::unique_ptr<GDALDataset>(GDALDataset::Open(
1387 6 : pszRef, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR));
1388 : auto poSubInput =
1389 : std::unique_ptr<GDALDataset>(GDALDataset::Open(
1390 6 : pszInput, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR));
1391 3 : if (poSubRef && poSubInput)
1392 : {
1393 : const uint64_t nSubDSPixels =
1394 3 : static_cast<uint64_t>(poSubRef->GetRasterXSize()) *
1395 3 : poSubRef->GetRasterYSize() * poSubRef->GetRasterCount();
1396 6 : void *pScaledProgress = GDALCreateScaledProgress(
1397 3 : static_cast<double>(nIterPixels) /
1398 3 : static_cast<double>(nTotalPixels),
1399 3 : static_cast<double>(nIterPixels + nSubDSPixels) /
1400 3 : static_cast<double>(nTotalPixels),
1401 : ctxt.m_pfnProgress, ctxt.m_pProgressData);
1402 3 : DatasetComparison(
1403 : aosReport, poSubRef.get(), poSubInput.get(),
1404 : pScaledProgress ? GDALScaledProgress : nullptr,
1405 : pScaledProgress);
1406 3 : GDALDestroyScaledProgress(pScaledProgress);
1407 3 : nIterPixels += nSubDSPixels;
1408 : }
1409 : }
1410 : }
1411 : }
1412 :
1413 384 : for (const auto &s : aosReport)
1414 : {
1415 182 : m_output += s;
1416 182 : m_output += '\n';
1417 : }
1418 :
1419 202 : m_retCode = static_cast<int>(aosReport.size());
1420 :
1421 202 : return true;
1422 : }
1423 :
1424 : /************************************************************************/
1425 : /* ~GDALRasterCompareAlgorithmStandalone() */
1426 : /************************************************************************/
1427 :
1428 : GDALRasterCompareAlgorithmStandalone::~GDALRasterCompareAlgorithmStandalone() =
1429 : default;
1430 :
1431 : //! @endcond
|
