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