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