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