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