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