Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: High Performance Image Reprojector
4 : * Purpose: Implementation of the GDALWarpOperation class.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2003, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2007-2012, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_port.h"
15 : #include "gdalwarper.h"
16 :
17 : #include <cctype>
18 : #include <climits>
19 : #include <cmath>
20 : #include <cstddef>
21 : #include <cstdlib>
22 : #include <cstring>
23 :
24 : #include <algorithm>
25 : #include <limits>
26 : #include <map>
27 : #include <memory>
28 : #include <mutex>
29 :
30 : #include "cpl_config.h"
31 : #include "cpl_conv.h"
32 : #include "cpl_error.h"
33 : #include "cpl_error_internal.h"
34 : #include "cpl_mask.h"
35 : #include "cpl_multiproc.h"
36 : #include "cpl_string.h"
37 : #include "cpl_vsi.h"
38 : #include "gdal.h"
39 : #include "gdal_priv.h"
40 : #include "gdal_alg_priv.h"
41 : #include "ogr_api.h"
42 : #include "ogr_core.h"
43 :
44 : struct _GDALWarpChunk
45 : {
46 : int dx, dy, dsx, dsy;
47 : int sx, sy, ssx, ssy;
48 : double sExtraSx, sExtraSy;
49 : };
50 :
51 : struct GDALWarpPrivateData
52 : {
53 : int nStepCount = 0;
54 : std::vector<int> abSuccess{};
55 : std::vector<double> adfDstX{};
56 : std::vector<double> adfDstY{};
57 : };
58 :
59 : static std::mutex gMutex{};
60 : static std::map<GDALWarpOperation *, std::unique_ptr<GDALWarpPrivateData>>
61 : gMapPrivate{};
62 :
63 : static GDALWarpPrivateData *
64 1260 : GetWarpPrivateData(GDALWarpOperation *poWarpOperation)
65 : {
66 2520 : std::lock_guard<std::mutex> oLock(gMutex);
67 1260 : auto oItem = gMapPrivate.find(poWarpOperation);
68 1260 : if (oItem != gMapPrivate.end())
69 : {
70 948 : return oItem->second.get();
71 : }
72 : else
73 : {
74 312 : gMapPrivate[poWarpOperation] =
75 624 : std::unique_ptr<GDALWarpPrivateData>(new GDALWarpPrivateData());
76 312 : return gMapPrivate[poWarpOperation].get();
77 : }
78 : }
79 :
80 : /************************************************************************/
81 : /* ==================================================================== */
82 : /* GDALWarpOperation */
83 : /* ==================================================================== */
84 : /************************************************************************/
85 :
86 : /**
87 : * \class GDALWarpOperation "gdalwarper.h"
88 : *
89 : * High level image warping class.
90 :
91 : <h2>Warper Design</h2>
92 :
93 : The overall GDAL high performance image warper is split into a few components.
94 :
95 : - The transformation between input and output file coordinates is handled
96 : via GDALTransformerFunc() implementations such as the one returned by
97 : GDALCreateGenImgProjTransformer(). The transformers are ultimately responsible
98 : for translating pixel/line locations on the destination image to pixel/line
99 : locations on the source image.
100 :
101 : - In order to handle images too large to hold in RAM, the warper needs to
102 : segment large images. This is the responsibility of the GDALWarpOperation
103 : class. The GDALWarpOperation::ChunkAndWarpImage() invokes
104 : GDALWarpOperation::WarpRegion() on chunks of output and input image that
105 : are small enough to hold in the amount of memory allowed by the application.
106 : This process is described in greater detail in the <b>Image Chunking</b>
107 : section.
108 :
109 : - The GDALWarpOperation::WarpRegion() function creates and loads an output
110 : image buffer, and then calls WarpRegionToBuffer().
111 :
112 : - GDALWarpOperation::WarpRegionToBuffer() is responsible for loading the
113 : source imagery corresponding to a particular output region, and generating
114 : masks and density masks from the source and destination imagery using
115 : the generator functions found in the GDALWarpOptions structure. Binds this
116 : all into an instance of GDALWarpKernel on which the
117 : GDALWarpKernel::PerformWarp() method is called.
118 :
119 : - GDALWarpKernel does the actual image warping, but is given an input image
120 : and an output image to operate on. The GDALWarpKernel does no IO, and in
121 : fact knows nothing about GDAL. It invokes the transformation function to
122 : get sample locations, builds output values based on the resampling algorithm
123 : in use. It also takes any validity and density masks into account during
124 : this operation.
125 :
126 : <h3>Chunk Size Selection</h3>
127 :
128 : The GDALWarpOptions ChunkAndWarpImage() method is responsible for invoking
129 : the WarpRegion() method on appropriate sized output chunks such that the
130 : memory required for the output image buffer, input image buffer and any
131 : required density and validity buffers is less than or equal to the application
132 : defined maximum memory available for use.
133 :
134 : It checks the memory required by walking the edges of the output region,
135 : transforming the locations back into source pixel/line coordinates and
136 : establishing a bounding rectangle of source imagery that would be required
137 : for the output area. This is actually accomplished by the private
138 : GDALWarpOperation::ComputeSourceWindow() method.
139 :
140 : Then memory requirements are used by totaling the memory required for all
141 : output bands, input bands, validity masks and density masks. If this is
142 : greater than the GDALWarpOptions::dfWarpMemoryLimit then the destination
143 : region is divided in two (splitting the longest dimension), and
144 : ChunkAndWarpImage() recursively invoked on each destination subregion.
145 :
146 : <h3>Validity and Density Masks Generation</h3>
147 :
148 : Fill in ways in which the validity and density masks may be generated here.
149 : Note that detailed semantics of the masks should be found in
150 : GDALWarpKernel.
151 : */
152 :
153 : /************************************************************************/
154 : /* GDALWarpOperation() */
155 : /************************************************************************/
156 :
157 : GDALWarpOperation::GDALWarpOperation() = default;
158 :
159 : /************************************************************************/
160 : /* ~GDALWarpOperation() */
161 : /************************************************************************/
162 :
163 1831 : GDALWarpOperation::~GDALWarpOperation()
164 :
165 : {
166 : {
167 3662 : std::lock_guard<std::mutex> oLock(gMutex);
168 1831 : auto oItem = gMapPrivate.find(this);
169 1831 : if (oItem != gMapPrivate.end())
170 : {
171 312 : gMapPrivate.erase(oItem);
172 : }
173 : }
174 :
175 1831 : WipeOptions();
176 :
177 1831 : if (hIOMutex != nullptr)
178 : {
179 6 : CPLDestroyMutex(hIOMutex);
180 6 : CPLDestroyMutex(hWarpMutex);
181 : }
182 :
183 1831 : WipeChunkList();
184 1831 : if (psThreadData)
185 1827 : GWKThreadsEnd(psThreadData);
186 1831 : }
187 :
188 : /************************************************************************/
189 : /* GetOptions() */
190 : /************************************************************************/
191 :
192 : /** Return warp options */
193 4372 : const GDALWarpOptions *GDALWarpOperation::GetOptions()
194 :
195 : {
196 4372 : return psOptions;
197 : }
198 :
199 : /************************************************************************/
200 : /* WipeOptions() */
201 : /************************************************************************/
202 :
203 1835 : void GDALWarpOperation::WipeOptions()
204 :
205 : {
206 1835 : if (psOptions != nullptr)
207 : {
208 1831 : GDALDestroyWarpOptions(psOptions);
209 1831 : psOptions = nullptr;
210 : }
211 1835 : }
212 :
213 : /************************************************************************/
214 : /* ValidateOptions() */
215 : /************************************************************************/
216 :
217 1831 : int GDALWarpOperation::ValidateOptions()
218 :
219 : {
220 1831 : if (psOptions == nullptr)
221 : {
222 0 : CPLError(CE_Failure, CPLE_IllegalArg,
223 : "GDALWarpOptions.Validate(): "
224 : "no options currently initialized.");
225 0 : return FALSE;
226 : }
227 :
228 1831 : if (psOptions->dfWarpMemoryLimit < 100000.0)
229 : {
230 0 : CPLError(CE_Failure, CPLE_IllegalArg,
231 : "GDALWarpOptions.Validate(): "
232 : "dfWarpMemoryLimit=%g is unreasonably small.",
233 0 : psOptions->dfWarpMemoryLimit);
234 0 : return FALSE;
235 : }
236 :
237 1831 : if (psOptions->eResampleAlg != GRA_NearestNeighbour &&
238 1048 : psOptions->eResampleAlg != GRA_Bilinear &&
239 624 : psOptions->eResampleAlg != GRA_Cubic &&
240 329 : psOptions->eResampleAlg != GRA_CubicSpline &&
241 265 : psOptions->eResampleAlg != GRA_Lanczos &&
242 207 : psOptions->eResampleAlg != GRA_Average &&
243 109 : psOptions->eResampleAlg != GRA_RMS &&
244 99 : psOptions->eResampleAlg != GRA_Mode &&
245 53 : psOptions->eResampleAlg != GRA_Max &&
246 47 : psOptions->eResampleAlg != GRA_Min &&
247 42 : psOptions->eResampleAlg != GRA_Med &&
248 36 : psOptions->eResampleAlg != GRA_Q1 &&
249 26 : psOptions->eResampleAlg != GRA_Q3 && psOptions->eResampleAlg != GRA_Sum)
250 : {
251 0 : CPLError(CE_Failure, CPLE_IllegalArg,
252 : "GDALWarpOptions.Validate(): "
253 : "eResampleArg=%d is not a supported value.",
254 0 : psOptions->eResampleAlg);
255 0 : return FALSE;
256 : }
257 :
258 1831 : if (static_cast<int>(psOptions->eWorkingDataType) < 1 ||
259 1831 : static_cast<int>(psOptions->eWorkingDataType) >= GDT_TypeCount)
260 : {
261 0 : CPLError(CE_Failure, CPLE_IllegalArg,
262 : "GDALWarpOptions.Validate(): "
263 : "eWorkingDataType=%d is not a supported value.",
264 0 : psOptions->eWorkingDataType);
265 0 : return FALSE;
266 : }
267 :
268 2084 : if (GDALDataTypeIsComplex(psOptions->eWorkingDataType) != 0 &&
269 253 : (psOptions->eResampleAlg == GRA_Max ||
270 253 : psOptions->eResampleAlg == GRA_Min ||
271 253 : psOptions->eResampleAlg == GRA_Med ||
272 253 : psOptions->eResampleAlg == GRA_Q1 ||
273 253 : psOptions->eResampleAlg == GRA_Q3))
274 : {
275 :
276 0 : CPLError(CE_Failure, CPLE_NotSupported,
277 : "GDALWarpOptions.Validate(): "
278 : "min/max/qnt not supported for complex valued data.");
279 0 : return FALSE;
280 : }
281 :
282 1831 : if (psOptions->hSrcDS == nullptr)
283 : {
284 0 : CPLError(CE_Failure, CPLE_IllegalArg,
285 : "GDALWarpOptions.Validate(): "
286 : "hSrcDS is not set.");
287 0 : return FALSE;
288 : }
289 :
290 1831 : if (psOptions->nBandCount == 0)
291 : {
292 0 : CPLError(CE_Failure, CPLE_IllegalArg,
293 : "GDALWarpOptions.Validate(): "
294 : "nBandCount=0, no bands configured!");
295 0 : return FALSE;
296 : }
297 :
298 1831 : if (psOptions->panSrcBands == nullptr)
299 : {
300 0 : CPLError(CE_Failure, CPLE_IllegalArg,
301 : "GDALWarpOptions.Validate(): "
302 : "panSrcBands is NULL.");
303 0 : return FALSE;
304 : }
305 :
306 1831 : if (psOptions->hDstDS != nullptr && psOptions->panDstBands == nullptr)
307 : {
308 0 : CPLError(CE_Failure, CPLE_IllegalArg,
309 : "GDALWarpOptions.Validate(): "
310 : "panDstBands is NULL.");
311 0 : return FALSE;
312 : }
313 :
314 4509 : for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
315 : {
316 5356 : if (psOptions->panSrcBands[iBand] < 1 ||
317 2678 : psOptions->panSrcBands[iBand] >
318 2678 : GDALGetRasterCount(psOptions->hSrcDS))
319 : {
320 0 : CPLError(CE_Failure, CPLE_IllegalArg,
321 : "panSrcBands[%d] = %d ... out of range for dataset.",
322 0 : iBand, psOptions->panSrcBands[iBand]);
323 0 : return FALSE;
324 : }
325 5343 : if (psOptions->hDstDS != nullptr &&
326 2665 : (psOptions->panDstBands[iBand] < 1 ||
327 2665 : psOptions->panDstBands[iBand] >
328 2665 : GDALGetRasterCount(psOptions->hDstDS)))
329 : {
330 0 : CPLError(CE_Failure, CPLE_IllegalArg,
331 : "panDstBands[%d] = %d ... out of range for dataset.",
332 0 : iBand, psOptions->panDstBands[iBand]);
333 0 : return FALSE;
334 : }
335 :
336 5343 : if (psOptions->hDstDS != nullptr &&
337 2665 : GDALGetRasterAccess(GDALGetRasterBand(
338 2665 : psOptions->hDstDS, psOptions->panDstBands[iBand])) ==
339 : GA_ReadOnly)
340 : {
341 0 : CPLError(CE_Failure, CPLE_IllegalArg,
342 : "Destination band %d appears to be read-only.",
343 0 : psOptions->panDstBands[iBand]);
344 0 : return FALSE;
345 : }
346 : }
347 :
348 1831 : if (psOptions->nBandCount == 0)
349 : {
350 0 : CPLError(CE_Failure, CPLE_IllegalArg,
351 : "GDALWarpOptions.Validate(): "
352 : "nBandCount=0, no bands configured!");
353 0 : return FALSE;
354 : }
355 :
356 1831 : if (psOptions->pfnProgress == nullptr)
357 : {
358 0 : CPLError(CE_Failure, CPLE_IllegalArg,
359 : "GDALWarpOptions.Validate(): "
360 : "pfnProgress is NULL.");
361 0 : return FALSE;
362 : }
363 :
364 1831 : if (psOptions->pfnTransformer == nullptr)
365 : {
366 0 : CPLError(CE_Failure, CPLE_IllegalArg,
367 : "GDALWarpOptions.Validate(): "
368 : "pfnTransformer is NULL.");
369 0 : return FALSE;
370 : }
371 :
372 : {
373 3662 : CPLStringList aosWO(CSLDuplicate(psOptions->papszWarpOptions));
374 : // A few internal/undocumented options
375 1831 : aosWO.SetNameValue("EXTRA_ELTS", nullptr);
376 1831 : aosWO.SetNameValue("USE_GENERAL_CASE", nullptr);
377 1831 : aosWO.SetNameValue("ERROR_THRESHOLD", nullptr);
378 1831 : aosWO.SetNameValue("ERROR_OUT_IF_EMPTY_SOURCE_WINDOW", nullptr);
379 1831 : aosWO.SetNameValue("MULT_FACTOR_VERTICAL_SHIFT_PIPELINE", nullptr);
380 1831 : aosWO.SetNameValue("SRC_FILL_RATIO_HEURISTICS", nullptr);
381 1831 : GDALValidateOptions(nullptr, GDALWarpGetOptionList(), aosWO.List(),
382 : "option", "warp options");
383 : }
384 :
385 : const char *pszSampleSteps =
386 1831 : CSLFetchNameValue(psOptions->papszWarpOptions, "SAMPLE_STEPS");
387 1831 : if (pszSampleSteps)
388 : {
389 8 : if (!EQUAL(pszSampleSteps, "ALL") && atoi(pszSampleSteps) < 2)
390 : {
391 0 : CPLError(CE_Failure, CPLE_IllegalArg,
392 : "GDALWarpOptions.Validate(): "
393 : "SAMPLE_STEPS warp option has illegal value.");
394 0 : return FALSE;
395 : }
396 : }
397 :
398 1831 : if (psOptions->nSrcAlphaBand > 0)
399 : {
400 214 : if (psOptions->hSrcDS == nullptr ||
401 107 : psOptions->nSrcAlphaBand > GDALGetRasterCount(psOptions->hSrcDS))
402 : {
403 0 : CPLError(CE_Failure, CPLE_IllegalArg,
404 : "nSrcAlphaBand = %d ... out of range for dataset.",
405 0 : psOptions->nSrcAlphaBand);
406 0 : return FALSE;
407 : }
408 : }
409 :
410 1831 : if (psOptions->nDstAlphaBand > 0)
411 : {
412 808 : if (psOptions->hDstDS == nullptr ||
413 404 : psOptions->nDstAlphaBand > GDALGetRasterCount(psOptions->hDstDS))
414 : {
415 0 : CPLError(CE_Failure, CPLE_IllegalArg,
416 : "nDstAlphaBand = %d ... out of range for dataset.",
417 0 : psOptions->nDstAlphaBand);
418 0 : return FALSE;
419 : }
420 : }
421 :
422 1831 : if (psOptions->nSrcAlphaBand > 0 &&
423 107 : psOptions->pfnSrcDensityMaskFunc != nullptr)
424 : {
425 0 : CPLError(CE_Failure, CPLE_IllegalArg,
426 : "GDALWarpOptions.Validate(): "
427 : "pfnSrcDensityMaskFunc provided as well as a SrcAlphaBand.");
428 0 : return FALSE;
429 : }
430 :
431 1831 : if (psOptions->nDstAlphaBand > 0 &&
432 404 : psOptions->pfnDstDensityMaskFunc != nullptr)
433 : {
434 0 : CPLError(CE_Failure, CPLE_IllegalArg,
435 : "GDALWarpOptions.Validate(): "
436 : "pfnDstDensityMaskFunc provided as well as a DstAlphaBand.");
437 0 : return FALSE;
438 : }
439 :
440 : GDALRasterBandH hSrcBand =
441 1831 : GDALGetRasterBand(psOptions->hSrcDS, psOptions->panSrcBands[0]);
442 1833 : if (GDALGetMaskFlags(hSrcBand) == GMF_PER_DATASET &&
443 2 : psOptions->padfSrcNoDataReal != nullptr)
444 : {
445 1 : CPLError(
446 : CE_Warning, CPLE_AppDefined,
447 : "Source dataset has both a per-dataset mask band and the warper "
448 : "has been also configured with a source nodata value. Only taking "
449 : "into account the latter (i.e. ignoring the per-dataset mask "
450 : "band)");
451 : }
452 :
453 3662 : const bool bErrorOutIfEmptySourceWindow = CPLFetchBool(
454 1831 : psOptions->papszWarpOptions, "ERROR_OUT_IF_EMPTY_SOURCE_WINDOW", true);
455 2219 : if (!bErrorOutIfEmptySourceWindow &&
456 388 : CSLFetchNameValue(psOptions->papszWarpOptions, "INIT_DEST") == nullptr)
457 : {
458 0 : CPLError(CE_Failure, CPLE_IllegalArg,
459 : "GDALWarpOptions.Validate(): "
460 : "ERROR_OUT_IF_EMPTY_SOURCE_WINDOW=FALSE can only be used "
461 : "if INIT_DEST is set");
462 0 : return FALSE;
463 : }
464 :
465 1831 : return TRUE;
466 : }
467 :
468 : /************************************************************************/
469 : /* SetAlphaMax() */
470 : /************************************************************************/
471 :
472 510 : static void SetAlphaMax(GDALWarpOptions *psOptions, GDALRasterBandH hBand,
473 : const char *pszKey)
474 : {
475 : const char *pszNBits =
476 510 : GDALGetMetadataItem(hBand, GDALMD_NBITS, GDAL_MDD_IMAGE_STRUCTURE);
477 510 : const char *pszAlphaMax = nullptr;
478 510 : if (pszNBits)
479 : {
480 4 : pszAlphaMax = CPLSPrintf("%u", (1U << atoi(pszNBits)) - 1U);
481 : }
482 506 : else if (GDALGetRasterDataType(hBand) == GDT_Int16)
483 : {
484 20 : pszAlphaMax = "32767";
485 : }
486 486 : else if (GDALGetRasterDataType(hBand) == GDT_UInt16)
487 : {
488 22 : pszAlphaMax = "65535";
489 : }
490 :
491 510 : if (pszAlphaMax != nullptr)
492 46 : psOptions->papszWarpOptions =
493 46 : CSLSetNameValue(psOptions->papszWarpOptions, pszKey, pszAlphaMax);
494 : else
495 464 : CPLDebug("WARP", "SetAlphaMax: AlphaMax not set.");
496 510 : }
497 :
498 : /************************************************************************/
499 : /* SetTieStrategy() */
500 : /************************************************************************/
501 :
502 1831 : static void SetTieStrategy(GDALWarpOptions *psOptions, CPLErr *peErr)
503 : {
504 1831 : if (const char *pszTieStrategy =
505 1831 : CSLFetchNameValue(psOptions->papszWarpOptions, "MODE_TIES"))
506 : {
507 14 : if (EQUAL(pszTieStrategy, "FIRST"))
508 : {
509 4 : psOptions->eTieStrategy = GWKTS_First;
510 : }
511 10 : else if (EQUAL(pszTieStrategy, "MIN"))
512 : {
513 4 : psOptions->eTieStrategy = GWKTS_Min;
514 : }
515 6 : else if (EQUAL(pszTieStrategy, "MAX"))
516 : {
517 4 : psOptions->eTieStrategy = GWKTS_Max;
518 : }
519 : else
520 : {
521 2 : CPLError(CE_Failure, CPLE_IllegalArg,
522 : "Unknown value of MODE_TIES: %s", pszTieStrategy);
523 2 : *peErr = CE_Failure;
524 : }
525 : }
526 1831 : }
527 :
528 : /************************************************************************/
529 : /* Initialize() */
530 : /************************************************************************/
531 :
532 : /**
533 : * \fn CPLErr GDALWarpOperation::Initialize( const GDALWarpOptions * );
534 : *
535 : * This method initializes the GDALWarpOperation's concept of the warp
536 : * options in effect. It creates an internal copy of the GDALWarpOptions
537 : * structure and defaults a variety of additional fields in the internal
538 : * copy if not set in the provided warp options.
539 : *
540 : * Defaulting operations include:
541 : * - If the nBandCount is 0, it will be set to the number of bands in the
542 : * source image (which must match the output image) and the panSrcBands
543 : * and panDstBands will be populated.
544 : *
545 : * @param psNewOptions input set of warp options. These are copied and may
546 : * be destroyed after this call by the application.
547 : * @param pfnTransformer Transformer function that this GDALWarpOperation must use
548 : * and own, or NULL. When pfnTransformer is not NULL, this implies that
549 : * psNewOptions->pfnTransformer is NULL
550 : * @param psOwnedTransformerArg Transformer argument that this GDALWarpOperation
551 : * must use, and own, or NULL. When psOwnedTransformerArg is set, this implies that
552 : * psNewOptions->pTransformerArg is NULL
553 : *
554 : * @return CE_None on success or CE_Failure if an error occurs.
555 : */
556 :
557 : CPLErr
558 1831 : GDALWarpOperation::Initialize(const GDALWarpOptions *psNewOptions,
559 : GDALTransformerFunc pfnTransformer,
560 : GDALTransformerArgUniquePtr psOwnedTransformerArg)
561 :
562 : {
563 : /* -------------------------------------------------------------------- */
564 : /* Copy the passed in options. */
565 : /* -------------------------------------------------------------------- */
566 1831 : if (psOptions != nullptr)
567 0 : WipeOptions();
568 :
569 1831 : CPLErr eErr = CE_None;
570 :
571 1831 : psOptions = GDALCloneWarpOptions(psNewOptions);
572 :
573 1831 : if (psOptions->pfnTransformer)
574 : {
575 854 : CPLAssert(pfnTransformer == nullptr);
576 854 : CPLAssert(psOwnedTransformerArg.get() == nullptr);
577 : }
578 : else
579 : {
580 977 : m_psOwnedTransformerArg = std::move(psOwnedTransformerArg);
581 977 : psOptions->pfnTransformer = pfnTransformer;
582 977 : psOptions->pTransformerArg = m_psOwnedTransformerArg.get();
583 : }
584 :
585 3662 : psOptions->papszWarpOptions =
586 1831 : CSLSetNameValue(psOptions->papszWarpOptions, "EXTRA_ELTS",
587 : CPLSPrintf("%d", WARP_EXTRA_ELTS));
588 :
589 : /* -------------------------------------------------------------------- */
590 : /* Default band mapping if missing. */
591 : /* -------------------------------------------------------------------- */
592 0 : if (psOptions->nBandCount == 0 && psOptions->hSrcDS != nullptr &&
593 1831 : psOptions->hDstDS != nullptr &&
594 0 : GDALGetRasterCount(psOptions->hSrcDS) ==
595 0 : GDALGetRasterCount(psOptions->hDstDS))
596 : {
597 0 : GDALWarpInitDefaultBandMapping(psOptions,
598 0 : GDALGetRasterCount(psOptions->hSrcDS));
599 : }
600 :
601 1831 : GDALWarpResolveWorkingDataType(psOptions);
602 1831 : SetTieStrategy(psOptions, &eErr);
603 :
604 : /* -------------------------------------------------------------------- */
605 : /* Default memory available. */
606 : /* */
607 : /* For now we default to 64MB of RAM, but eventually we should */
608 : /* try various schemes to query physical RAM. This can */
609 : /* certainly be done on Win32 and Linux. */
610 : /* -------------------------------------------------------------------- */
611 1831 : if (psOptions->dfWarpMemoryLimit == 0.0)
612 : {
613 1622 : psOptions->dfWarpMemoryLimit = 64.0 * 1024 * 1024;
614 : }
615 :
616 : /* -------------------------------------------------------------------- */
617 : /* Are we doing timings? */
618 : /* -------------------------------------------------------------------- */
619 1831 : bReportTimings =
620 1831 : CPLFetchBool(psOptions->papszWarpOptions, "REPORT_TIMINGS", false);
621 :
622 : /* -------------------------------------------------------------------- */
623 : /* Support creating cutline from text warpoption. */
624 : /* -------------------------------------------------------------------- */
625 : const char *pszCutlineWKT =
626 1831 : CSLFetchNameValue(psOptions->papszWarpOptions, "CUTLINE");
627 :
628 1831 : if (pszCutlineWKT && psOptions->hCutline == nullptr)
629 : {
630 44 : char *pszWKTTmp = const_cast<char *>(pszCutlineWKT);
631 88 : if (OGR_G_CreateFromWkt(&pszWKTTmp, nullptr,
632 : reinterpret_cast<OGRGeometryH *>(
633 44 : &(psOptions->hCutline))) != OGRERR_NONE)
634 : {
635 2 : eErr = CE_Failure;
636 2 : CPLError(CE_Failure, CPLE_AppDefined,
637 : "Failed to parse CUTLINE geometry wkt.");
638 : }
639 : }
640 : const char *pszBD =
641 1831 : CSLFetchNameValue(psOptions->papszWarpOptions, "CUTLINE_BLEND_DIST");
642 1831 : if (pszBD)
643 0 : psOptions->dfCutlineBlendDist = CPLAtof(pszBD);
644 :
645 : /* -------------------------------------------------------------------- */
646 : /* Set SRC_ALPHA_MAX if not provided. */
647 : /* -------------------------------------------------------------------- */
648 1831 : if (psOptions->hSrcDS != nullptr && psOptions->nSrcAlphaBand > 0 &&
649 3769 : psOptions->nSrcAlphaBand <= GDALGetRasterCount(psOptions->hSrcDS) &&
650 107 : CSLFetchNameValue(psOptions->papszWarpOptions, "SRC_ALPHA_MAX") ==
651 : nullptr)
652 : {
653 : GDALRasterBandH hSrcAlphaBand =
654 107 : GDALGetRasterBand(psOptions->hSrcDS, psOptions->nSrcAlphaBand);
655 107 : SetAlphaMax(psOptions, hSrcAlphaBand, "SRC_ALPHA_MAX");
656 : }
657 :
658 : /* -------------------------------------------------------------------- */
659 : /* Set DST_ALPHA_MAX if not provided. */
660 : /* -------------------------------------------------------------------- */
661 1818 : if (psOptions->hDstDS != nullptr && psOptions->nDstAlphaBand > 0 &&
662 4053 : psOptions->nDstAlphaBand <= GDALGetRasterCount(psOptions->hDstDS) &&
663 404 : CSLFetchNameValue(psOptions->papszWarpOptions, "DST_ALPHA_MAX") ==
664 : nullptr)
665 : {
666 : GDALRasterBandH hDstAlphaBand =
667 403 : GDALGetRasterBand(psOptions->hDstDS, psOptions->nDstAlphaBand);
668 403 : SetAlphaMax(psOptions, hDstAlphaBand, "DST_ALPHA_MAX");
669 : }
670 :
671 : /* -------------------------------------------------------------------- */
672 : /* If the options don't validate, then wipe them. */
673 : /* -------------------------------------------------------------------- */
674 1831 : if (!ValidateOptions())
675 0 : eErr = CE_Failure;
676 :
677 1831 : if (eErr != CE_None)
678 : {
679 4 : WipeOptions();
680 : }
681 : else
682 : {
683 3654 : psThreadData = GWKThreadsCreate(psOptions->papszWarpOptions,
684 1827 : psOptions->pfnTransformer,
685 1827 : psOptions->pTransformerArg);
686 1827 : if (psThreadData == nullptr)
687 0 : eErr = CE_Failure;
688 :
689 : /* --------------------------------------------------------------------
690 : */
691 : /* Compute dstcoordinates of a few special points. */
692 : /* --------------------------------------------------------------------
693 : */
694 :
695 : // South and north poles. Do not exactly take +/-90 as the
696 : // round-tripping of the longitude value fails with some projections.
697 5481 : for (double dfY : {-89.9999, 89.9999})
698 : {
699 3654 : double dfX = 0;
700 3654 : if ((GDALIsTransformer(psOptions->pTransformerArg,
701 3034 : GDAL_APPROX_TRANSFORMER_CLASS_NAME) &&
702 3034 : GDALTransformLonLatToDestApproxTransformer(
703 7308 : psOptions->pTransformerArg, &dfX, &dfY)) ||
704 2090 : (GDALIsTransformer(psOptions->pTransformerArg,
705 314 : GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME) &&
706 314 : GDALTransformLonLatToDestGenImgProjTransformer(
707 314 : psOptions->pTransformerArg, &dfX, &dfY)))
708 : {
709 : aDstXYSpecialPoints.emplace_back(
710 1620 : std::pair<double, double>(dfX, dfY));
711 : }
712 : }
713 :
714 1827 : m_bIsTranslationOnPixelBoundaries =
715 3654 : GDALTransformIsTranslationOnPixelBoundaries(
716 2090 : psOptions->pfnTransformer, psOptions->pTransformerArg) &&
717 263 : CPLTestBool(
718 : CPLGetConfigOption("GDAL_WARP_USE_TRANSLATION_OPTIM", "YES"));
719 1827 : if (m_bIsTranslationOnPixelBoundaries)
720 : {
721 257 : CPLDebug("WARP",
722 : "Using translation-on-pixel-boundaries optimization");
723 : }
724 : }
725 :
726 1831 : if (eErr == CE_None && psOptions->hDstDS)
727 : {
728 : const auto oResetDestPixels =
729 : cpl::strict_parse<bool>(CSLFetchNameValueDef(
730 1814 : psOptions->papszWarpOptions, "RESET_DEST_PIXELS", "NO"));
731 :
732 1814 : if (!oResetDestPixels.has_value())
733 : {
734 1 : CPLError(CE_Failure, CPLE_IllegalArg,
735 : "Invalid value of RESET_DEST_PIXELS");
736 1 : return CE_Failure;
737 : }
738 :
739 1813 : bool bResetDestPixels = false;
740 : try
741 : {
742 1813 : bResetDestPixels = oResetDestPixels.value();
743 : }
744 0 : catch (const std::exception &)
745 : {
746 : // to make Coverity Scan happy
747 : }
748 :
749 1813 : if (bResetDestPixels)
750 : {
751 4 : for (int i = 0; eErr == CE_None && i < psOptions->nBandCount; ++i)
752 : {
753 2 : eErr =
754 5 : GDALFillRaster(GDALGetRasterBand(psOptions->hDstDS,
755 2 : psOptions->panDstBands[i]),
756 2 : psOptions->padfDstNoDataReal
757 1 : ? psOptions->padfDstNoDataReal[i]
758 : : 0.0,
759 2 : psOptions->padfDstNoDataImag
760 0 : ? psOptions->padfDstNoDataImag[i]
761 : : 0.0);
762 : }
763 : }
764 : }
765 :
766 1830 : return eErr;
767 : }
768 :
769 : /**
770 : * \fn void* GDALWarpOperation::CreateDestinationBuffer(
771 : int nDstXSize, int nDstYSize, int *pbInitialized);
772 : *
773 : * This method creates a destination buffer for use with WarpRegionToBuffer.
774 : * The output is initialized based on the INIT_DEST settings.
775 : *
776 : * @param nDstXSize Width of output window on destination buffer to be produced.
777 : * @param nDstYSize Height of output window on destination buffer to be
778 : produced.
779 : * @param pbInitialized Filled with boolean indicating if the buffer was
780 : initialized.
781 : *
782 : * @return Buffer capable for use as a warp operation output destination
783 : */
784 3386 : void *GDALWarpOperation::CreateDestinationBuffer(int nDstXSize, int nDstYSize,
785 : int *pbInitialized)
786 : {
787 :
788 : /* -------------------------------------------------------------------- */
789 : /* Allocate block of memory large enough to hold all the bands */
790 : /* for this block. */
791 : /* -------------------------------------------------------------------- */
792 3386 : const int nWordSize = GDALGetDataTypeSizeBytes(psOptions->eWorkingDataType);
793 :
794 3386 : void *pDstBuffer = VSI_MALLOC3_VERBOSE(
795 : cpl::fits_on<int>(nWordSize * psOptions->nBandCount), nDstXSize,
796 : nDstYSize);
797 3386 : if (pDstBuffer)
798 : {
799 3386 : auto eErr = InitializeDestinationBuffer(pDstBuffer, nDstXSize,
800 : nDstYSize, pbInitialized);
801 3386 : if (eErr != CE_None)
802 : {
803 2 : CPLFree(pDstBuffer);
804 2 : return nullptr;
805 : }
806 : }
807 3384 : return pDstBuffer;
808 : }
809 :
810 : /**
811 : * This method initializes a destination buffer for use with WarpRegionToBuffer.
812 : *
813 : * It is initialized based on the INIT_DEST settings.
814 : *
815 : * This method is called by CreateDestinationBuffer().
816 : * It is meant at being used by callers that have already allocated the
817 : * destination buffer without using CreateDestinationBuffer().
818 : *
819 : * @param pDstBuffer Buffer of size
820 : * GDALGetDataTypeSizeBytes(psOptions->eWorkingDataType) *
821 : * nDstXSize * nDstYSize * psOptions->nBandCount bytes.
822 : * @param nDstXSize Width of output window on destination buffer to be produced.
823 : * @param nDstYSize Height of output window on destination buffer to be
824 : * produced.
825 : * @param pbInitialized Filled with boolean indicating if the buffer was
826 : * initialized.
827 : * @since 3.10
828 : */
829 3485 : CPLErr GDALWarpOperation::InitializeDestinationBuffer(void *pDstBuffer,
830 : int nDstXSize,
831 : int nDstYSize,
832 : int *pbInitialized) const
833 : {
834 3485 : const int nWordSize = GDALGetDataTypeSizeBytes(psOptions->eWorkingDataType);
835 :
836 3485 : const GPtrDiff_t nBandSize =
837 3485 : static_cast<GPtrDiff_t>(nWordSize) * nDstXSize * nDstYSize;
838 :
839 : /* -------------------------------------------------------------------- */
840 : /* Initialize if requested in the options */
841 : /* -------------------------------------------------------------------- */
842 : const char *pszInitDest =
843 3485 : CSLFetchNameValue(psOptions->papszWarpOptions, "INIT_DEST");
844 :
845 3485 : if (pszInitDest == nullptr || EQUAL(pszInitDest, ""))
846 : {
847 384 : if (pbInitialized != nullptr)
848 : {
849 384 : *pbInitialized = FALSE;
850 : }
851 384 : return CE_None;
852 : }
853 :
854 3101 : if (pbInitialized != nullptr)
855 : {
856 1938 : *pbInitialized = TRUE;
857 : }
858 :
859 : CPLStringList aosInitValues(
860 6202 : CSLTokenizeStringComplex(pszInitDest, ",", FALSE, FALSE));
861 3101 : const int nInitCount = aosInitValues.Count();
862 :
863 8306 : for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
864 : {
865 5207 : double adfInitRealImag[2] = {0.0, 0.0};
866 : const char *pszBandInit =
867 5207 : aosInitValues[std::min(iBand, nInitCount - 1)];
868 :
869 5207 : if (EQUAL(pszBandInit, "NO_DATA"))
870 : {
871 709 : if (psOptions->padfDstNoDataReal == nullptr)
872 : {
873 1 : CPLError(CE_Failure, CPLE_AppDefined,
874 : "INIT_DEST was set to NO_DATA, but a NoData value was "
875 : "not defined.");
876 : }
877 : else
878 : {
879 708 : adfInitRealImag[0] = psOptions->padfDstNoDataReal[iBand];
880 708 : if (psOptions->padfDstNoDataImag != nullptr)
881 : {
882 601 : adfInitRealImag[1] = psOptions->padfDstNoDataImag[iBand];
883 : }
884 : }
885 : }
886 : else
887 : {
888 4498 : if (CPLStringToComplex(pszBandInit, &adfInitRealImag[0],
889 4498 : &adfInitRealImag[1]) != CE_None)
890 : {
891 2 : CPLError(CE_Failure, CPLE_AppDefined,
892 : "Error parsing INIT_DEST");
893 2 : return CE_Failure;
894 : }
895 : }
896 :
897 5205 : GByte *pBandData = static_cast<GByte *>(pDstBuffer) + iBand * nBandSize;
898 :
899 5205 : if (psOptions->eWorkingDataType == GDT_UInt8)
900 : {
901 9112 : memset(pBandData,
902 : std::max(
903 4556 : 0, std::min(255, static_cast<int>(adfInitRealImag[0]))),
904 : nBandSize);
905 : }
906 1293 : else if (!std::isnan(adfInitRealImag[0]) && adfInitRealImag[0] == 0.0 &&
907 1293 : !std::isnan(adfInitRealImag[1]) && adfInitRealImag[1] == 0.0)
908 : {
909 559 : memset(pBandData, 0, nBandSize);
910 : }
911 90 : else if (!std::isnan(adfInitRealImag[1]) && adfInitRealImag[1] == 0.0)
912 : {
913 90 : GDALCopyWords64(&adfInitRealImag, GDT_Float64, 0, pBandData,
914 90 : psOptions->eWorkingDataType, nWordSize,
915 90 : static_cast<GPtrDiff_t>(nDstXSize) * nDstYSize);
916 : }
917 : else
918 : {
919 0 : GDALCopyWords64(&adfInitRealImag, GDT_CFloat64, 0, pBandData,
920 0 : psOptions->eWorkingDataType, nWordSize,
921 0 : static_cast<GPtrDiff_t>(nDstXSize) * nDstYSize);
922 : }
923 : }
924 :
925 3099 : return CE_None;
926 : }
927 :
928 : /**
929 : * \fn void GDALWarpOperation::DestroyDestinationBuffer( void *pDstBuffer )
930 : *
931 : * This method destroys a buffer previously retrieved from
932 : * CreateDestinationBuffer
933 : *
934 : * @param pDstBuffer destination buffer to be destroyed
935 : *
936 : */
937 3384 : void GDALWarpOperation::DestroyDestinationBuffer(void *pDstBuffer)
938 : {
939 3384 : VSIFree(pDstBuffer);
940 3384 : }
941 :
942 : /************************************************************************/
943 : /* GDALCreateWarpOperation() */
944 : /************************************************************************/
945 :
946 : /**
947 : * @see GDALWarpOperation::Initialize()
948 : */
949 :
950 149 : GDALWarpOperationH GDALCreateWarpOperation(const GDALWarpOptions *psNewOptions)
951 : {
952 149 : GDALWarpOperation *poOperation = new GDALWarpOperation;
953 149 : if (poOperation->Initialize(psNewOptions) != CE_None)
954 : {
955 0 : delete poOperation;
956 0 : return nullptr;
957 : }
958 :
959 149 : return reinterpret_cast<GDALWarpOperationH>(poOperation);
960 : }
961 :
962 : /************************************************************************/
963 : /* GDALDestroyWarpOperation() */
964 : /************************************************************************/
965 :
966 : /**
967 : * @see GDALWarpOperation::~GDALWarpOperation()
968 : */
969 :
970 149 : void GDALDestroyWarpOperation(GDALWarpOperationH hOperation)
971 : {
972 149 : if (hOperation)
973 149 : delete static_cast<GDALWarpOperation *>(hOperation);
974 149 : }
975 :
976 : /************************************************************************/
977 : /* CollectChunkList() */
978 : /************************************************************************/
979 :
980 1226 : void GDALWarpOperation::CollectChunkList(int nDstXOff, int nDstYOff,
981 : int nDstXSize, int nDstYSize)
982 :
983 : {
984 : /* -------------------------------------------------------------------- */
985 : /* Collect the list of chunks to operate on. */
986 : /* -------------------------------------------------------------------- */
987 1226 : WipeChunkList();
988 1226 : CollectChunkListInternal(nDstXOff, nDstYOff, nDstXSize, nDstYSize);
989 :
990 : // Sort chunks from top to bottom, and for equal y, from left to right.
991 1226 : if (nChunkListCount > 1)
992 : {
993 56 : std::sort(pasChunkList, pasChunkList + nChunkListCount,
994 7201 : [](const GDALWarpChunk &a, const GDALWarpChunk &b)
995 : {
996 7201 : if (a.dy < b.dy)
997 3220 : return true;
998 3981 : if (a.dy > b.dy)
999 1415 : return false;
1000 2566 : return a.dx < b.dx;
1001 : });
1002 : }
1003 :
1004 : /* -------------------------------------------------------------------- */
1005 : /* Find the global source window. */
1006 : /* -------------------------------------------------------------------- */
1007 :
1008 1226 : const int knIntMax = std::numeric_limits<int>::max();
1009 1226 : const int knIntMin = std::numeric_limits<int>::min();
1010 1226 : int nSrcXOff = knIntMax;
1011 1226 : int nSrcYOff = knIntMax;
1012 1226 : int nSrcX2Off = knIntMin;
1013 1226 : int nSrcY2Off = knIntMin;
1014 1226 : double dfApproxAccArea = 0;
1015 3548 : for (int iChunk = 0; pasChunkList != nullptr && iChunk < nChunkListCount;
1016 : iChunk++)
1017 : {
1018 2322 : GDALWarpChunk *pasThisChunk = pasChunkList + iChunk;
1019 2322 : nSrcXOff = std::min(nSrcXOff, pasThisChunk->sx);
1020 2322 : nSrcYOff = std::min(nSrcYOff, pasThisChunk->sy);
1021 2322 : nSrcX2Off = std::max(nSrcX2Off, pasThisChunk->sx + pasThisChunk->ssx);
1022 2322 : nSrcY2Off = std::max(nSrcY2Off, pasThisChunk->sy + pasThisChunk->ssy);
1023 2322 : dfApproxAccArea +=
1024 2322 : static_cast<double>(pasThisChunk->ssx) * pasThisChunk->ssy;
1025 : }
1026 1226 : if (nSrcXOff < nSrcX2Off)
1027 : {
1028 1218 : const double dfTotalArea =
1029 1218 : static_cast<double>(nSrcX2Off - nSrcXOff) * (nSrcY2Off - nSrcYOff);
1030 : // This is really a gross heuristics, but should work in most cases
1031 1218 : if (dfApproxAccArea >= dfTotalArea * 0.80)
1032 : {
1033 1218 : GDALDataset::FromHandle(psOptions->hSrcDS)
1034 1218 : ->AdviseRead(nSrcXOff, nSrcYOff, nSrcX2Off - nSrcXOff,
1035 : nSrcY2Off - nSrcYOff, nDstXSize, nDstYSize,
1036 1218 : psOptions->eWorkingDataType, psOptions->nBandCount,
1037 1218 : psOptions->panSrcBands, nullptr);
1038 : }
1039 : }
1040 1226 : }
1041 :
1042 : /************************************************************************/
1043 : /* ChunkAndWarpImage() */
1044 : /************************************************************************/
1045 :
1046 : /**
1047 : * \fn CPLErr GDALWarpOperation::ChunkAndWarpImage(
1048 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize );
1049 : *
1050 : * This method does a complete warp of the source image to the destination
1051 : * image for the indicated region with the current warp options in effect.
1052 : * Progress is reported to the installed progress monitor, if any.
1053 : *
1054 : * This function will subdivide the region and recursively call itself
1055 : * until the total memory required to process a region chunk will all fit
1056 : * in the memory pool defined by GDALWarpOptions::dfWarpMemoryLimit.
1057 : *
1058 : * Once an appropriate region is selected GDALWarpOperation::WarpRegion()
1059 : * is invoked to do the actual work.
1060 : *
1061 : * @param nDstXOff X offset to window of destination data to be produced.
1062 : * @param nDstYOff Y offset to window of destination data to be produced.
1063 : * @param nDstXSize Width of output window on destination file to be produced.
1064 : * @param nDstYSize Height of output window on destination file to be produced.
1065 : *
1066 : * @return CE_None on success or CE_Failure if an error occurs.
1067 : */
1068 :
1069 1220 : CPLErr GDALWarpOperation::ChunkAndWarpImage(int nDstXOff, int nDstYOff,
1070 : int nDstXSize, int nDstYSize)
1071 :
1072 : {
1073 : /* -------------------------------------------------------------------- */
1074 : /* Collect the list of chunks to operate on. */
1075 : /* -------------------------------------------------------------------- */
1076 1220 : CollectChunkList(nDstXOff, nDstYOff, nDstXSize, nDstYSize);
1077 :
1078 : /* -------------------------------------------------------------------- */
1079 : /* Total up output pixels to process. */
1080 : /* -------------------------------------------------------------------- */
1081 1220 : double dfTotalPixels = 0.0;
1082 :
1083 3531 : for (int iChunk = 0; pasChunkList != nullptr && iChunk < nChunkListCount;
1084 : iChunk++)
1085 : {
1086 2311 : GDALWarpChunk *pasThisChunk = pasChunkList + iChunk;
1087 2311 : const double dfChunkPixels =
1088 2311 : pasThisChunk->dsx * static_cast<double>(pasThisChunk->dsy);
1089 :
1090 2311 : dfTotalPixels += dfChunkPixels;
1091 : }
1092 :
1093 : /* -------------------------------------------------------------------- */
1094 : /* Process them one at a time, updating the progress */
1095 : /* information for each region. */
1096 : /* -------------------------------------------------------------------- */
1097 1220 : double dfPixelsProcessed = 0.0;
1098 :
1099 3525 : for (int iChunk = 0; pasChunkList != nullptr && iChunk < nChunkListCount;
1100 : iChunk++)
1101 : {
1102 2311 : GDALWarpChunk *pasThisChunk = pasChunkList + iChunk;
1103 2311 : const double dfChunkPixels =
1104 2311 : pasThisChunk->dsx * static_cast<double>(pasThisChunk->dsy);
1105 :
1106 2311 : const double dfProgressBase = dfPixelsProcessed / dfTotalPixels;
1107 2311 : const double dfProgressScale = dfChunkPixels / dfTotalPixels;
1108 :
1109 2311 : CPLErr eErr = WarpRegion(
1110 : pasThisChunk->dx, pasThisChunk->dy, pasThisChunk->dsx,
1111 : pasThisChunk->dsy, pasThisChunk->sx, pasThisChunk->sy,
1112 : pasThisChunk->ssx, pasThisChunk->ssy, pasThisChunk->sExtraSx,
1113 : pasThisChunk->sExtraSy, dfProgressBase, dfProgressScale);
1114 :
1115 2311 : if (eErr != CE_None)
1116 6 : return eErr;
1117 :
1118 2305 : dfPixelsProcessed += dfChunkPixels;
1119 : }
1120 :
1121 1214 : WipeChunkList();
1122 :
1123 1214 : psOptions->pfnProgress(1.0, "", psOptions->pProgressArg);
1124 :
1125 1214 : return CE_None;
1126 : }
1127 :
1128 : /************************************************************************/
1129 : /* GDALChunkAndWarpImage() */
1130 : /************************************************************************/
1131 :
1132 : /**
1133 : * @see GDALWarpOperation::ChunkAndWarpImage()
1134 : */
1135 :
1136 149 : CPLErr GDALChunkAndWarpImage(GDALWarpOperationH hOperation, int nDstXOff,
1137 : int nDstYOff, int nDstXSize, int nDstYSize)
1138 : {
1139 149 : VALIDATE_POINTER1(hOperation, "GDALChunkAndWarpImage", CE_Failure);
1140 :
1141 : return reinterpret_cast<GDALWarpOperation *>(hOperation)
1142 149 : ->ChunkAndWarpImage(nDstXOff, nDstYOff, nDstXSize, nDstYSize);
1143 : }
1144 :
1145 : /************************************************************************/
1146 : /* ChunkThreadMain() */
1147 : /************************************************************************/
1148 :
1149 : struct ChunkThreadData
1150 : {
1151 : GDALWarpOperation *poOperation = nullptr;
1152 : GDALWarpChunk *pasChunkInfo = nullptr;
1153 : CPLJoinableThread *hThreadHandle = nullptr;
1154 : CPLErr eErr = CE_None;
1155 : double dfProgressBase = 0;
1156 : double dfProgressScale = 0;
1157 : CPLMutex *hIOMutex = nullptr;
1158 :
1159 : CPLMutex *hCondMutex = nullptr;
1160 : volatile int bIOMutexTaken = 0;
1161 : CPLCond *hCond = nullptr;
1162 :
1163 : CPLErrorAccumulator *poErrorAccumulator = nullptr;
1164 : };
1165 :
1166 11 : static void ChunkThreadMain(void *pThreadData)
1167 :
1168 : {
1169 11 : volatile ChunkThreadData *psData =
1170 : static_cast<volatile ChunkThreadData *>(pThreadData);
1171 :
1172 11 : GDALWarpChunk *pasChunkInfo = psData->pasChunkInfo;
1173 :
1174 : /* -------------------------------------------------------------------- */
1175 : /* Acquire IO mutex. */
1176 : /* -------------------------------------------------------------------- */
1177 11 : if (!CPLAcquireMutex(psData->hIOMutex, 600.0))
1178 : {
1179 0 : CPLError(CE_Failure, CPLE_AppDefined,
1180 : "Failed to acquire IOMutex in WarpRegion().");
1181 0 : psData->eErr = CE_Failure;
1182 : }
1183 : else
1184 : {
1185 11 : if (psData->hCond != nullptr)
1186 : {
1187 6 : CPLAcquireMutex(psData->hCondMutex, 1.0);
1188 6 : psData->bIOMutexTaken = TRUE;
1189 6 : CPLCondSignal(psData->hCond);
1190 6 : CPLReleaseMutex(psData->hCondMutex);
1191 : }
1192 :
1193 : auto oAccumulator =
1194 22 : psData->poErrorAccumulator->InstallForCurrentScope();
1195 11 : CPL_IGNORE_RET_VAL(oAccumulator);
1196 :
1197 22 : psData->eErr = psData->poOperation->WarpRegion(
1198 : pasChunkInfo->dx, pasChunkInfo->dy, pasChunkInfo->dsx,
1199 : pasChunkInfo->dsy, pasChunkInfo->sx, pasChunkInfo->sy,
1200 : pasChunkInfo->ssx, pasChunkInfo->ssy, pasChunkInfo->sExtraSx,
1201 11 : pasChunkInfo->sExtraSy, psData->dfProgressBase,
1202 11 : psData->dfProgressScale);
1203 :
1204 : /* --------------------------------------------------------------------
1205 : */
1206 : /* Release the IO mutex. */
1207 : /* --------------------------------------------------------------------
1208 : */
1209 11 : CPLReleaseMutex(psData->hIOMutex);
1210 : }
1211 11 : }
1212 :
1213 : /************************************************************************/
1214 : /* ChunkAndWarpMulti() */
1215 : /************************************************************************/
1216 :
1217 : /**
1218 : * \fn CPLErr GDALWarpOperation::ChunkAndWarpMulti(
1219 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize );
1220 : *
1221 : * This method does a complete warp of the source image to the destination
1222 : * image for the indicated region with the current warp options in effect.
1223 : * Progress is reported to the installed progress monitor, if any.
1224 : *
1225 : * Externally this method operates the same as ChunkAndWarpImage(), but
1226 : * internally this method uses multiple threads to interleave input/output
1227 : * for one region while the processing is being done for another.
1228 : *
1229 : * @param nDstXOff X offset to window of destination data to be produced.
1230 : * @param nDstYOff Y offset to window of destination data to be produced.
1231 : * @param nDstXSize Width of output window on destination file to be produced.
1232 : * @param nDstYSize Height of output window on destination file to be produced.
1233 : *
1234 : * @return CE_None on success or CE_Failure if an error occurs.
1235 : */
1236 :
1237 6 : CPLErr GDALWarpOperation::ChunkAndWarpMulti(int nDstXOff, int nDstYOff,
1238 : int nDstXSize, int nDstYSize)
1239 :
1240 : {
1241 6 : hIOMutex = CPLCreateMutex();
1242 6 : hWarpMutex = CPLCreateMutex();
1243 :
1244 6 : CPLReleaseMutex(hIOMutex);
1245 6 : CPLReleaseMutex(hWarpMutex);
1246 :
1247 6 : CPLCond *hCond = CPLCreateCond();
1248 6 : CPLMutex *hCondMutex = CPLCreateMutex();
1249 6 : CPLReleaseMutex(hCondMutex);
1250 :
1251 : /* -------------------------------------------------------------------- */
1252 : /* Collect the list of chunks to operate on. */
1253 : /* -------------------------------------------------------------------- */
1254 6 : CollectChunkList(nDstXOff, nDstYOff, nDstXSize, nDstYSize);
1255 :
1256 : /* -------------------------------------------------------------------- */
1257 : /* Process them one at a time, updating the progress */
1258 : /* information for each region. */
1259 : /* -------------------------------------------------------------------- */
1260 6 : ChunkThreadData volatile asThreadData[2] = {};
1261 6 : CPLErrorAccumulator oErrorAccumulator;
1262 18 : for (int i = 0; i < 2; ++i)
1263 : {
1264 12 : asThreadData[i].poOperation = this;
1265 12 : asThreadData[i].hIOMutex = hIOMutex;
1266 12 : asThreadData[i].poErrorAccumulator = &oErrorAccumulator;
1267 : }
1268 :
1269 6 : double dfPixelsProcessed = 0.0;
1270 6 : double dfTotalPixels = static_cast<double>(nDstXSize) * nDstYSize;
1271 :
1272 6 : CPLErr eErr = CE_None;
1273 22 : for (int iChunk = 0; iChunk < nChunkListCount + 1; iChunk++)
1274 : {
1275 17 : int iThread = iChunk % 2;
1276 :
1277 : /* --------------------------------------------------------------------
1278 : */
1279 : /* Launch thread for this chunk. */
1280 : /* --------------------------------------------------------------------
1281 : */
1282 17 : if (pasChunkList != nullptr && iChunk < nChunkListCount)
1283 : {
1284 11 : GDALWarpChunk *pasThisChunk = pasChunkList + iChunk;
1285 11 : const double dfChunkPixels =
1286 11 : pasThisChunk->dsx * static_cast<double>(pasThisChunk->dsy);
1287 :
1288 11 : asThreadData[iThread].dfProgressBase =
1289 11 : dfPixelsProcessed / dfTotalPixels;
1290 11 : asThreadData[iThread].dfProgressScale =
1291 11 : dfChunkPixels / dfTotalPixels;
1292 :
1293 11 : dfPixelsProcessed += dfChunkPixels;
1294 :
1295 11 : asThreadData[iThread].pasChunkInfo = pasThisChunk;
1296 :
1297 11 : if (iChunk == 0)
1298 : {
1299 6 : asThreadData[iThread].hCond = hCond;
1300 6 : asThreadData[iThread].hCondMutex = hCondMutex;
1301 : }
1302 : else
1303 : {
1304 5 : asThreadData[iThread].hCond = nullptr;
1305 5 : asThreadData[iThread].hCondMutex = nullptr;
1306 : }
1307 11 : asThreadData[iThread].bIOMutexTaken = FALSE;
1308 :
1309 11 : CPLDebug("GDAL", "Start chunk %d / %d.", iChunk, nChunkListCount);
1310 22 : asThreadData[iThread].hThreadHandle = CPLCreateJoinableThread(
1311 : ChunkThreadMain,
1312 11 : const_cast<ChunkThreadData *>(&asThreadData[iThread]));
1313 11 : if (asThreadData[iThread].hThreadHandle == nullptr)
1314 : {
1315 0 : CPLError(
1316 : CE_Failure, CPLE_AppDefined,
1317 : "CPLCreateJoinableThread() failed in ChunkAndWarpMulti()");
1318 0 : eErr = CE_Failure;
1319 0 : break;
1320 : }
1321 :
1322 : // Wait that the first thread has acquired the IO mutex before
1323 : // proceeding. This will ensure that the first thread will run
1324 : // before the second one.
1325 11 : if (iChunk == 0)
1326 : {
1327 6 : CPLAcquireMutex(hCondMutex, 1.0);
1328 11 : while (asThreadData[iThread].bIOMutexTaken == FALSE)
1329 5 : CPLCondWait(hCond, hCondMutex);
1330 6 : CPLReleaseMutex(hCondMutex);
1331 : }
1332 : }
1333 :
1334 : /* --------------------------------------------------------------------
1335 : */
1336 : /* Wait for previous chunks thread to complete. */
1337 : /* --------------------------------------------------------------------
1338 : */
1339 17 : if (iChunk > 0)
1340 : {
1341 11 : iThread = (iChunk - 1) % 2;
1342 :
1343 : // Wait for thread to finish.
1344 11 : CPLJoinThread(asThreadData[iThread].hThreadHandle);
1345 11 : asThreadData[iThread].hThreadHandle = nullptr;
1346 :
1347 11 : CPLDebug("GDAL", "Finished chunk %d / %d.", iChunk - 1,
1348 : nChunkListCount);
1349 :
1350 11 : eErr = asThreadData[iThread].eErr;
1351 :
1352 11 : if (eErr != CE_None)
1353 1 : break;
1354 : }
1355 : }
1356 :
1357 : /* -------------------------------------------------------------------- */
1358 : /* Wait for all threads to complete. */
1359 : /* -------------------------------------------------------------------- */
1360 18 : for (int iThread = 0; iThread < 2; iThread++)
1361 : {
1362 12 : if (asThreadData[iThread].hThreadHandle)
1363 0 : CPLJoinThread(asThreadData[iThread].hThreadHandle);
1364 : }
1365 :
1366 6 : CPLDestroyCond(hCond);
1367 6 : CPLDestroyMutex(hCondMutex);
1368 :
1369 6 : WipeChunkList();
1370 :
1371 6 : oErrorAccumulator.ReplayErrors();
1372 :
1373 6 : psOptions->pfnProgress(1.0, "", psOptions->pProgressArg);
1374 :
1375 12 : return eErr;
1376 : }
1377 :
1378 : /************************************************************************/
1379 : /* GDALChunkAndWarpMulti() */
1380 : /************************************************************************/
1381 :
1382 : /**
1383 : * @see GDALWarpOperation::ChunkAndWarpMulti()
1384 : */
1385 :
1386 0 : CPLErr GDALChunkAndWarpMulti(GDALWarpOperationH hOperation, int nDstXOff,
1387 : int nDstYOff, int nDstXSize, int nDstYSize)
1388 : {
1389 0 : VALIDATE_POINTER1(hOperation, "GDALChunkAndWarpMulti", CE_Failure);
1390 :
1391 : return reinterpret_cast<GDALWarpOperation *>(hOperation)
1392 0 : ->ChunkAndWarpMulti(nDstXOff, nDstYOff, nDstXSize, nDstYSize);
1393 : }
1394 :
1395 : /************************************************************************/
1396 : /* WipeChunkList() */
1397 : /************************************************************************/
1398 :
1399 4277 : void GDALWarpOperation::WipeChunkList()
1400 :
1401 : {
1402 4277 : CPLFree(pasChunkList);
1403 4277 : pasChunkList = nullptr;
1404 4277 : nChunkListCount = 0;
1405 4277 : nChunkListMax = 0;
1406 4277 : }
1407 :
1408 : /************************************************************************/
1409 : /* GetWorkingMemoryForWindow() */
1410 : /************************************************************************/
1411 :
1412 : /** Returns the amount of working memory, in bytes, required to process
1413 : * a warped window of source dimensions nSrcXSize x nSrcYSize and target
1414 : * dimensions nDstXSize x nDstYSize.
1415 : */
1416 4113 : double GDALWarpOperation::GetWorkingMemoryForWindow(int nSrcXSize,
1417 : int nSrcYSize,
1418 : int nDstXSize,
1419 : int nDstYSize) const
1420 : {
1421 : /* -------------------------------------------------------------------- */
1422 : /* Based on the types of masks in use, how many bits will each */
1423 : /* source pixel cost us? */
1424 : /* -------------------------------------------------------------------- */
1425 : int nSrcPixelCostInBits =
1426 4113 : GDALGetDataTypeSizeBits(psOptions->eWorkingDataType) *
1427 4113 : psOptions->nBandCount;
1428 :
1429 4113 : if (psOptions->pfnSrcDensityMaskFunc != nullptr)
1430 0 : nSrcPixelCostInBits += 32; // Float mask?
1431 :
1432 4113 : GDALRasterBandH hSrcBand = nullptr;
1433 4113 : if (psOptions->nBandCount > 0)
1434 : hSrcBand =
1435 4113 : GDALGetRasterBand(psOptions->hSrcDS, psOptions->panSrcBands[0]);
1436 :
1437 4113 : if (psOptions->nSrcAlphaBand > 0 || psOptions->hCutline != nullptr)
1438 118 : nSrcPixelCostInBits += 32; // UnifiedSrcDensity float mask.
1439 7990 : else if (hSrcBand != nullptr &&
1440 3995 : (GDALGetMaskFlags(hSrcBand) & GMF_PER_DATASET))
1441 6 : nSrcPixelCostInBits += 1; // UnifiedSrcValid bit mask.
1442 :
1443 4113 : if (psOptions->papfnSrcPerBandValidityMaskFunc != nullptr ||
1444 4113 : psOptions->padfSrcNoDataReal != nullptr)
1445 211 : nSrcPixelCostInBits += psOptions->nBandCount; // Bit/band mask.
1446 :
1447 4113 : if (psOptions->pfnSrcValidityMaskFunc != nullptr)
1448 0 : nSrcPixelCostInBits += 1; // Bit mask.
1449 :
1450 : /* -------------------------------------------------------------------- */
1451 : /* What about the cost for the destination. */
1452 : /* -------------------------------------------------------------------- */
1453 : int nDstPixelCostInBits =
1454 4113 : GDALGetDataTypeSizeBits(psOptions->eWorkingDataType) *
1455 4113 : psOptions->nBandCount;
1456 :
1457 4113 : if (psOptions->pfnDstDensityMaskFunc != nullptr)
1458 0 : nDstPixelCostInBits += 32;
1459 :
1460 4113 : if (psOptions->padfDstNoDataReal != nullptr ||
1461 2435 : psOptions->pfnDstValidityMaskFunc != nullptr)
1462 1678 : nDstPixelCostInBits += psOptions->nBandCount;
1463 :
1464 4113 : if (psOptions->nDstAlphaBand > 0)
1465 256 : nDstPixelCostInBits += 32; // DstDensity float mask.
1466 :
1467 4113 : const double dfTotalMemoryUse =
1468 4113 : (static_cast<double>(nSrcPixelCostInBits) * nSrcXSize * nSrcYSize +
1469 4113 : static_cast<double>(nDstPixelCostInBits) * nDstXSize * nDstYSize) /
1470 : 8.0;
1471 4113 : return dfTotalMemoryUse;
1472 : }
1473 :
1474 : /************************************************************************/
1475 : /* CollectChunkListInternal() */
1476 : /************************************************************************/
1477 :
1478 4402 : CPLErr GDALWarpOperation::CollectChunkListInternal(int nDstXOff, int nDstYOff,
1479 : int nDstXSize, int nDstYSize)
1480 :
1481 : {
1482 : /* -------------------------------------------------------------------- */
1483 : /* Compute the bounds of the input area corresponding to the */
1484 : /* output area. */
1485 : /* -------------------------------------------------------------------- */
1486 4402 : int nSrcXOff = 0;
1487 4402 : int nSrcYOff = 0;
1488 4402 : int nSrcXSize = 0;
1489 4402 : int nSrcYSize = 0;
1490 4402 : double dfSrcXExtraSize = 0.0;
1491 4402 : double dfSrcYExtraSize = 0.0;
1492 4402 : double dfSrcFillRatio = 0.0;
1493 : CPLErr eErr;
1494 : {
1495 4402 : CPLTurnFailureIntoWarningBackuper oBackuper;
1496 4402 : eErr = ComputeSourceWindow(nDstXOff, nDstYOff, nDstXSize, nDstYSize,
1497 : &nSrcXOff, &nSrcYOff, &nSrcXSize, &nSrcYSize,
1498 : &dfSrcXExtraSize, &dfSrcYExtraSize,
1499 : &dfSrcFillRatio);
1500 : }
1501 :
1502 4402 : if (eErr != CE_None)
1503 : {
1504 : const bool bErrorOutIfEmptySourceWindow =
1505 3 : CPLFetchBool(psOptions->papszWarpOptions,
1506 : "ERROR_OUT_IF_EMPTY_SOURCE_WINDOW", true);
1507 3 : if (bErrorOutIfEmptySourceWindow)
1508 : {
1509 3 : CPLError(CE_Warning, CPLE_AppDefined,
1510 : "Unable to compute source region for "
1511 : "output window %d,%d,%d,%d, skipping.",
1512 : nDstXOff, nDstYOff, nDstXSize, nDstYSize);
1513 : }
1514 : else
1515 : {
1516 0 : CPLDebug("WARP",
1517 : "Unable to compute source region for "
1518 : "output window %d,%d,%d,%d, skipping.",
1519 : nDstXOff, nDstYOff, nDstXSize, nDstYSize);
1520 : }
1521 : }
1522 :
1523 : /* -------------------------------------------------------------------- */
1524 : /* If we are allowed to drop no-source regions, do so now if */
1525 : /* appropriate. */
1526 : /* -------------------------------------------------------------------- */
1527 5527 : if ((nSrcXSize == 0 || nSrcYSize == 0) &&
1528 1125 : CPLFetchBool(psOptions->papszWarpOptions, "SKIP_NOSOURCE", false))
1529 492 : return CE_None;
1530 :
1531 : /* -------------------------------------------------------------------- */
1532 : /* Does the cost of the current rectangle exceed our memory */
1533 : /* limit? If so, split the destination along the longest */
1534 : /* dimension and recurse. */
1535 : /* -------------------------------------------------------------------- */
1536 : const double dfTotalMemoryUse =
1537 3910 : GetWorkingMemoryForWindow(nSrcXSize, nSrcYSize, nDstXSize, nDstYSize);
1538 :
1539 : // If size of working buffers need exceed the allow limit, then divide
1540 : // the target area
1541 : // Do it also if the "fill ratio" of the source is too low (#3120), but
1542 : // only if there's at least some source pixel intersecting. The
1543 : // SRC_FILL_RATIO_HEURISTICS warping option is undocumented and only here
1544 : // in case the heuristics would cause issues.
1545 : #if DEBUG_VERBOSE
1546 : CPLDebug("WARP",
1547 : "dst=(%d,%d,%d,%d) src=(%d,%d,%d,%d) srcfillratio=%.17g, "
1548 : "dfTotalMemoryUse=%.1f MB",
1549 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, nSrcXOff, nSrcYOff,
1550 : nSrcXSize, nSrcYSize, dfSrcFillRatio,
1551 : dfTotalMemoryUse / (1024 * 1024));
1552 : #endif
1553 870 : if ((dfTotalMemoryUse > psOptions->dfWarpMemoryLimit &&
1554 7820 : (nDstXSize > 2 || nDstYSize > 2)) ||
1555 3040 : (dfSrcFillRatio > 0 && dfSrcFillRatio < 0.5 &&
1556 350 : (nDstXSize > 100 || nDstYSize > 100) &&
1557 720 : CPLFetchBool(psOptions->papszWarpOptions, "SRC_FILL_RATIO_HEURISTICS",
1558 : true)))
1559 : {
1560 1589 : int nBlockXSize = 1;
1561 1589 : int nBlockYSize = 1;
1562 1589 : if (psOptions->hDstDS)
1563 : {
1564 1589 : GDALGetBlockSize(GDALGetRasterBand(psOptions->hDstDS, 1),
1565 : &nBlockXSize, &nBlockYSize);
1566 : }
1567 :
1568 1589 : int bStreamableOutput = CPLFetchBool(psOptions->papszWarpOptions,
1569 1589 : "STREAMABLE_OUTPUT", false);
1570 : const char *pszOptimizeSize =
1571 1589 : CSLFetchNameValue(psOptions->papszWarpOptions, "OPTIMIZE_SIZE");
1572 1589 : const bool bOptimizeSizeAuto =
1573 1589 : !pszOptimizeSize || EQUAL(pszOptimizeSize, "AUTO");
1574 : const bool bOptimizeSize =
1575 4436 : !bStreamableOutput &&
1576 97 : ((pszOptimizeSize && !bOptimizeSizeAuto &&
1577 1588 : CPLTestBool(pszOptimizeSize)) ||
1578 : // Auto-enable optimize-size mode if output region is at least
1579 : // 2x2 blocks large and the shapes of the source and target regions
1580 : // are not excessively different. All those thresholds are a bit
1581 : // arbitrary
1582 1491 : (bOptimizeSizeAuto && nSrcXSize > 0 && nDstYSize > 0 &&
1583 1259 : (nDstXSize > nDstYSize ? fabs(double(nDstXSize) / nDstYSize -
1584 506 : double(nSrcXSize) / nSrcYSize) <
1585 506 : 5 * double(nDstXSize) / nDstYSize
1586 753 : : fabs(double(nDstYSize) / nDstXSize -
1587 753 : double(nSrcYSize) / nSrcXSize) <
1588 753 : 5 * double(nDstYSize) / nDstXSize) &&
1589 1256 : nDstXSize / 2 >= nBlockXSize && nDstYSize / 2 >= nBlockYSize));
1590 :
1591 : // If the region width is greater than the region height,
1592 : // cut in half in the width. When we want to optimize the size
1593 : // of a compressed output dataset, do this only if each half part
1594 : // is at least as wide as the block width.
1595 1589 : bool bHasDivided = false;
1596 1589 : CPLErr eErr2 = CE_None;
1597 1589 : if (nDstXSize > nDstYSize &&
1598 658 : ((!bOptimizeSize && !bStreamableOutput) ||
1599 88 : (bOptimizeSize &&
1600 89 : (nDstXSize / 2 >= nBlockXSize || nDstYSize == 1)) ||
1601 1 : (bStreamableOutput && nDstXSize / 2 >= nBlockXSize &&
1602 1 : nDstYSize == nBlockYSize)))
1603 : {
1604 611 : bHasDivided = true;
1605 611 : int nChunk1 = nDstXSize / 2;
1606 :
1607 : // In the optimize size case, try to stick on target block
1608 : // boundaries.
1609 611 : if ((bOptimizeSize || bStreamableOutput) && nChunk1 > nBlockXSize)
1610 42 : nChunk1 = (nChunk1 / nBlockXSize) * nBlockXSize;
1611 :
1612 611 : int nChunk2 = nDstXSize - nChunk1;
1613 :
1614 611 : eErr = CollectChunkListInternal(nDstXOff, nDstYOff, nChunk1,
1615 : nDstYSize);
1616 :
1617 611 : eErr2 = CollectChunkListInternal(nDstXOff + nChunk1, nDstYOff,
1618 611 : nChunk2, nDstYSize);
1619 : }
1620 978 : else if (!(bStreamableOutput && nDstYSize / 2 < nBlockYSize))
1621 : {
1622 977 : bHasDivided = true;
1623 977 : int nChunk1 = nDstYSize / 2;
1624 :
1625 : // In the optimize size case, try to stick on target block
1626 : // boundaries.
1627 977 : if ((bOptimizeSize || bStreamableOutput) && nChunk1 > nBlockYSize)
1628 77 : nChunk1 = (nChunk1 / nBlockYSize) * nBlockYSize;
1629 :
1630 977 : const int nChunk2 = nDstYSize - nChunk1;
1631 :
1632 977 : eErr = CollectChunkListInternal(nDstXOff, nDstYOff, nDstXSize,
1633 : nChunk1);
1634 :
1635 977 : eErr2 = CollectChunkListInternal(nDstXOff, nDstYOff + nChunk1,
1636 : nDstXSize, nChunk2);
1637 : }
1638 :
1639 1589 : if (bHasDivided)
1640 : {
1641 1588 : if (eErr == CE_None)
1642 1588 : return eErr2;
1643 : else
1644 0 : return eErr;
1645 : }
1646 : }
1647 :
1648 : /* -------------------------------------------------------------------- */
1649 : /* OK, everything fits, so add to the chunk list. */
1650 : /* -------------------------------------------------------------------- */
1651 2322 : if (nChunkListCount == nChunkListMax)
1652 : {
1653 1387 : nChunkListMax = nChunkListMax * 2 + 1;
1654 1387 : pasChunkList = static_cast<GDALWarpChunk *>(
1655 1387 : CPLRealloc(pasChunkList, sizeof(GDALWarpChunk) * nChunkListMax));
1656 : }
1657 :
1658 2322 : pasChunkList[nChunkListCount].dx = nDstXOff;
1659 2322 : pasChunkList[nChunkListCount].dy = nDstYOff;
1660 2322 : pasChunkList[nChunkListCount].dsx = nDstXSize;
1661 2322 : pasChunkList[nChunkListCount].dsy = nDstYSize;
1662 2322 : pasChunkList[nChunkListCount].sx = nSrcXOff;
1663 2322 : pasChunkList[nChunkListCount].sy = nSrcYOff;
1664 2322 : pasChunkList[nChunkListCount].ssx = nSrcXSize;
1665 2322 : pasChunkList[nChunkListCount].ssy = nSrcYSize;
1666 2322 : pasChunkList[nChunkListCount].sExtraSx = dfSrcXExtraSize;
1667 2322 : pasChunkList[nChunkListCount].sExtraSy = dfSrcYExtraSize;
1668 :
1669 2322 : nChunkListCount++;
1670 :
1671 2322 : return CE_None;
1672 : }
1673 :
1674 : /************************************************************************/
1675 : /* WarpRegion() */
1676 : /************************************************************************/
1677 :
1678 : /**
1679 : * This method requests the indicated region of the output file be generated.
1680 : *
1681 : * Note that WarpRegion() will produce the requested area in one low level warp
1682 : * operation without verifying that this does not exceed the stated memory
1683 : * limits for the warp operation. Applications should take care not to call
1684 : * WarpRegion() on too large a region! This function
1685 : * is normally called by ChunkAndWarpImage(), the normal entry point for
1686 : * applications. Use it instead if staying within memory constraints is
1687 : * desired.
1688 : *
1689 : * Progress is reported from dfProgressBase to dfProgressBase + dfProgressScale
1690 : * for the indicated region.
1691 : *
1692 : * @param nDstXOff X offset to window of destination data to be produced.
1693 : * @param nDstYOff Y offset to window of destination data to be produced.
1694 : * @param nDstXSize Width of output window on destination file to be produced.
1695 : * @param nDstYSize Height of output window on destination file to be produced.
1696 : * @param nSrcXOff source window X offset (computed if window all zero)
1697 : * @param nSrcYOff source window Y offset (computed if window all zero)
1698 : * @param nSrcXSize source window X size (computed if window all zero)
1699 : * @param nSrcYSize source window Y size (computed if window all zero)
1700 : * @param dfProgressBase minimum progress value reported
1701 : * @param dfProgressScale value such as dfProgressBase + dfProgressScale is the
1702 : * maximum progress value reported
1703 : *
1704 : * @return CE_None on success or CE_Failure if an error occurs.
1705 : */
1706 :
1707 0 : CPLErr GDALWarpOperation::WarpRegion(int nDstXOff, int nDstYOff, int nDstXSize,
1708 : int nDstYSize, int nSrcXOff, int nSrcYOff,
1709 : int nSrcXSize, int nSrcYSize,
1710 : double dfProgressBase,
1711 : double dfProgressScale)
1712 : {
1713 0 : return WarpRegion(nDstXOff, nDstYOff, nDstXSize, nDstYSize, nSrcXOff,
1714 : nSrcYOff, nSrcXSize, nSrcYSize, 0, 0, dfProgressBase,
1715 0 : dfProgressScale);
1716 : }
1717 :
1718 : /**
1719 : * This method requests the indicated region of the output file be generated.
1720 : *
1721 : * Note that WarpRegion() will produce the requested area in one low level warp
1722 : * operation without verifying that this does not exceed the stated memory
1723 : * limits for the warp operation. Applications should take care not to call
1724 : * WarpRegion() on too large a region! This function
1725 : * is normally called by ChunkAndWarpImage(), the normal entry point for
1726 : * applications. Use it instead if staying within memory constraints is
1727 : * desired.
1728 : *
1729 : * Progress is reported from dfProgressBase to dfProgressBase + dfProgressScale
1730 : * for the indicated region.
1731 : *
1732 : * @param nDstXOff X offset to window of destination data to be produced.
1733 : * @param nDstYOff Y offset to window of destination data to be produced.
1734 : * @param nDstXSize Width of output window on destination file to be produced.
1735 : * @param nDstYSize Height of output window on destination file to be produced.
1736 : * @param nSrcXOff source window X offset (computed if window all zero)
1737 : * @param nSrcYOff source window Y offset (computed if window all zero)
1738 : * @param nSrcXSize source window X size (computed if window all zero)
1739 : * @param nSrcYSize source window Y size (computed if window all zero)
1740 : * @param dfSrcXExtraSize Extra pixels (included in nSrcXSize) reserved
1741 : * for filter window. Should be ignored in scale computation
1742 : * @param dfSrcYExtraSize Extra pixels (included in nSrcYSize) reserved
1743 : * for filter window. Should be ignored in scale computation
1744 : * @param dfProgressBase minimum progress value reported
1745 : * @param dfProgressScale value such as dfProgressBase + dfProgressScale is the
1746 : * maximum progress value reported
1747 : *
1748 : * @return CE_None on success or CE_Failure if an error occurs.
1749 : */
1750 :
1751 2322 : CPLErr GDALWarpOperation::WarpRegion(
1752 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize, int nSrcXOff,
1753 : int nSrcYOff, int nSrcXSize, int nSrcYSize, double dfSrcXExtraSize,
1754 : double dfSrcYExtraSize, double dfProgressBase, double dfProgressScale)
1755 :
1756 : {
1757 2322 : ReportTiming(nullptr);
1758 :
1759 : /* -------------------------------------------------------------------- */
1760 : /* Allocate the output buffer. */
1761 : /* -------------------------------------------------------------------- */
1762 2322 : int bDstBufferInitialized = FALSE;
1763 : void *pDstBuffer =
1764 2322 : CreateDestinationBuffer(nDstXSize, nDstYSize, &bDstBufferInitialized);
1765 2322 : if (pDstBuffer == nullptr)
1766 : {
1767 2 : return CE_Failure;
1768 : }
1769 :
1770 : /* -------------------------------------------------------------------- */
1771 : /* If we aren't doing fixed initialization of the output buffer */
1772 : /* then read it from disk so we can overlay on existing imagery. */
1773 : /* -------------------------------------------------------------------- */
1774 2320 : GDALDataset *poDstDS = GDALDataset::FromHandle(psOptions->hDstDS);
1775 2320 : if (!bDstBufferInitialized)
1776 : {
1777 384 : CPLErr eErr = CE_None;
1778 384 : if (psOptions->nBandCount == 1)
1779 : {
1780 : // Particular case to simplify the stack a bit.
1781 : // TODO(rouault): Need an explanation of what and why r34502 helps.
1782 360 : eErr = poDstDS->GetRasterBand(psOptions->panDstBands[0])
1783 720 : ->RasterIO(GF_Read, nDstXOff, nDstYOff, nDstXSize,
1784 : nDstYSize, pDstBuffer, nDstXSize, nDstYSize,
1785 360 : psOptions->eWorkingDataType, 0, 0, nullptr);
1786 : }
1787 : else
1788 : {
1789 24 : eErr = poDstDS->RasterIO(GF_Read, nDstXOff, nDstYOff, nDstXSize,
1790 : nDstYSize, pDstBuffer, nDstXSize,
1791 24 : nDstYSize, psOptions->eWorkingDataType,
1792 24 : psOptions->nBandCount,
1793 24 : psOptions->panDstBands, 0, 0, 0, nullptr);
1794 : }
1795 :
1796 384 : if (eErr != CE_None)
1797 : {
1798 0 : DestroyDestinationBuffer(pDstBuffer);
1799 0 : return eErr;
1800 : }
1801 :
1802 384 : ReportTiming("Output buffer read");
1803 : }
1804 :
1805 : /* -------------------------------------------------------------------- */
1806 : /* Perform the warp. */
1807 : /* -------------------------------------------------------------------- */
1808 : CPLErr eErr = nSrcXSize == 0
1809 2320 : ? CE_None
1810 1919 : : WarpRegionToBuffer(
1811 : nDstXOff, nDstYOff, nDstXSize, nDstYSize,
1812 1919 : pDstBuffer, psOptions->eWorkingDataType, nSrcXOff,
1813 : nSrcYOff, nSrcXSize, nSrcYSize, dfSrcXExtraSize,
1814 2320 : dfSrcYExtraSize, dfProgressBase, dfProgressScale);
1815 :
1816 : /* -------------------------------------------------------------------- */
1817 : /* Write the output data back to disk if all went well. */
1818 : /* -------------------------------------------------------------------- */
1819 2320 : if (eErr == CE_None)
1820 : {
1821 2315 : if (psOptions->nBandCount == 1)
1822 : {
1823 : // Particular case to simplify the stack a bit.
1824 2101 : eErr = poDstDS->GetRasterBand(psOptions->panDstBands[0])
1825 4202 : ->RasterIO(GF_Write, nDstXOff, nDstYOff, nDstXSize,
1826 : nDstYSize, pDstBuffer, nDstXSize, nDstYSize,
1827 2101 : psOptions->eWorkingDataType, 0, 0, nullptr);
1828 : }
1829 : else
1830 : {
1831 214 : eErr = poDstDS->RasterIO(GF_Write, nDstXOff, nDstYOff, nDstXSize,
1832 : nDstYSize, pDstBuffer, nDstXSize,
1833 214 : nDstYSize, psOptions->eWorkingDataType,
1834 214 : psOptions->nBandCount,
1835 214 : psOptions->panDstBands, 0, 0, 0, nullptr);
1836 : }
1837 :
1838 4630 : if (eErr == CE_None &&
1839 2315 : CPLFetchBool(psOptions->papszWarpOptions, "WRITE_FLUSH", false))
1840 : {
1841 0 : const CPLErr eOldErr = CPLGetLastErrorType();
1842 0 : const CPLString osLastErrMsg = CPLGetLastErrorMsg();
1843 0 : GDALFlushCache(psOptions->hDstDS);
1844 0 : const CPLErr eNewErr = CPLGetLastErrorType();
1845 0 : if (eNewErr != eOldErr ||
1846 0 : osLastErrMsg.compare(CPLGetLastErrorMsg()) != 0)
1847 0 : eErr = CE_Failure;
1848 : }
1849 2315 : ReportTiming("Output buffer write");
1850 : }
1851 :
1852 : /* -------------------------------------------------------------------- */
1853 : /* Cleanup and return. */
1854 : /* -------------------------------------------------------------------- */
1855 2320 : DestroyDestinationBuffer(pDstBuffer);
1856 :
1857 2320 : return eErr;
1858 : }
1859 :
1860 : /************************************************************************/
1861 : /* GDALWarpRegion() */
1862 : /************************************************************************/
1863 :
1864 : /**
1865 : * @see GDALWarpOperation::WarpRegion()
1866 : */
1867 :
1868 0 : CPLErr GDALWarpRegion(GDALWarpOperationH hOperation, int nDstXOff, int nDstYOff,
1869 : int nDstXSize, int nDstYSize, int nSrcXOff, int nSrcYOff,
1870 : int nSrcXSize, int nSrcYSize)
1871 :
1872 : {
1873 0 : VALIDATE_POINTER1(hOperation, "GDALWarpRegion", CE_Failure);
1874 :
1875 : return reinterpret_cast<GDALWarpOperation *>(hOperation)
1876 0 : ->WarpRegion(nDstXOff, nDstYOff, nDstXSize, nDstYSize, nSrcXOff,
1877 0 : nSrcYOff, nSrcXSize, nSrcYSize);
1878 : }
1879 :
1880 : /************************************************************************/
1881 : /* WarpRegionToBuffer() */
1882 : /************************************************************************/
1883 :
1884 : /**
1885 : * This method requests that a particular window of the output dataset
1886 : * be warped and the result put into the provided data buffer. The output
1887 : * dataset doesn't even really have to exist to use this method as long as
1888 : * the transformation function in the GDALWarpOptions is setup to map to
1889 : * a virtual pixel/line space.
1890 : *
1891 : * This method will do the whole region in one chunk, so be wary of the
1892 : * amount of memory that might be used.
1893 : *
1894 : * @param nDstXOff X offset to window of destination data to be produced.
1895 : * @param nDstYOff Y offset to window of destination data to be produced.
1896 : * @param nDstXSize Width of output window on destination file to be produced.
1897 : * @param nDstYSize Height of output window on destination file to be produced.
1898 : * @param pDataBuf the data buffer to place result in, of type eBufDataType.
1899 : * @param eBufDataType the type of the output data buffer. For now this
1900 : * must match GDALWarpOptions::eWorkingDataType.
1901 : * @param nSrcXOff source window X offset (computed if window all zero)
1902 : * @param nSrcYOff source window Y offset (computed if window all zero)
1903 : * @param nSrcXSize source window X size (computed if window all zero)
1904 : * @param nSrcYSize source window Y size (computed if window all zero)
1905 : * @param dfProgressBase minimum progress value reported
1906 : * @param dfProgressScale value such as dfProgressBase + dfProgressScale is the
1907 : * maximum progress value reported
1908 : *
1909 : * @return CE_None on success or CE_Failure if an error occurs.
1910 : */
1911 :
1912 1976 : CPLErr GDALWarpOperation::WarpRegionToBuffer(
1913 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize, void *pDataBuf,
1914 : GDALDataType eBufDataType, int nSrcXOff, int nSrcYOff, int nSrcXSize,
1915 : int nSrcYSize, double dfProgressBase, double dfProgressScale)
1916 : {
1917 1976 : return WarpRegionToBuffer(nDstXOff, nDstYOff, nDstXSize, nDstYSize,
1918 : pDataBuf, eBufDataType, nSrcXOff, nSrcYOff,
1919 : nSrcXSize, nSrcYSize, 0, 0, dfProgressBase,
1920 1976 : dfProgressScale);
1921 : }
1922 :
1923 : /**
1924 : * This method requests that a particular window of the output dataset
1925 : * be warped and the result put into the provided data buffer. The output
1926 : * dataset doesn't even really have to exist to use this method as long as
1927 : * the transformation function in the GDALWarpOptions is setup to map to
1928 : * a virtual pixel/line space.
1929 : *
1930 : * This method will do the whole region in one chunk, so be wary of the
1931 : * amount of memory that might be used.
1932 : *
1933 : * @param nDstXOff X offset to window of destination data to be produced.
1934 : * @param nDstYOff Y offset to window of destination data to be produced.
1935 : * @param nDstXSize Width of output window on destination file to be produced.
1936 : * @param nDstYSize Height of output window on destination file to be produced.
1937 : * @param pDataBuf the data buffer to place result in, of type eBufDataType.
1938 : * @param eBufDataType the type of the output data buffer. For now this
1939 : * must match GDALWarpOptions::eWorkingDataType.
1940 : * @param nSrcXOff source window X offset (computed if window all zero)
1941 : * @param nSrcYOff source window Y offset (computed if window all zero)
1942 : * @param nSrcXSize source window X size (computed if window all zero)
1943 : * @param nSrcYSize source window Y size (computed if window all zero)
1944 : * @param dfSrcXExtraSize Extra pixels (included in nSrcXSize) reserved
1945 : * for filter window. Should be ignored in scale computation
1946 : * @param dfSrcYExtraSize Extra pixels (included in nSrcYSize) reserved
1947 : * for filter window. Should be ignored in scale computation
1948 : * @param dfProgressBase minimum progress value reported
1949 : * @param dfProgressScale value such as dfProgressBase + dfProgressScale is the
1950 : * maximum progress value reported
1951 : *
1952 : * @return CE_None on success or CE_Failure if an error occurs.
1953 : */
1954 :
1955 3895 : CPLErr GDALWarpOperation::WarpRegionToBuffer(
1956 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize, void *pDataBuf,
1957 : // Only in a CPLAssert.
1958 : CPL_UNUSED GDALDataType eBufDataType, int nSrcXOff, int nSrcYOff,
1959 : int nSrcXSize, int nSrcYSize, double dfSrcXExtraSize,
1960 : double dfSrcYExtraSize, double dfProgressBase, double dfProgressScale)
1961 :
1962 : {
1963 3895 : const int nWordSize = GDALGetDataTypeSizeBytes(psOptions->eWorkingDataType);
1964 :
1965 3895 : CPLAssert(eBufDataType == psOptions->eWorkingDataType);
1966 :
1967 : /* -------------------------------------------------------------------- */
1968 : /* If not given a corresponding source window compute one now. */
1969 : /* -------------------------------------------------------------------- */
1970 3895 : if (nSrcXSize == 0 && nSrcYSize == 0)
1971 : {
1972 : // TODO: This taking of the warp mutex is suboptimal. We could get rid
1973 : // of it, but that would require making sure ComputeSourceWindow()
1974 : // uses a different pTransformerArg than the warp kernel.
1975 1777 : if (hWarpMutex != nullptr && !CPLAcquireMutex(hWarpMutex, 600.0))
1976 : {
1977 0 : CPLError(CE_Failure, CPLE_AppDefined,
1978 : "Failed to acquire WarpMutex in WarpRegion().");
1979 0 : return CE_Failure;
1980 : }
1981 : const CPLErr eErr =
1982 1777 : ComputeSourceWindow(nDstXOff, nDstYOff, nDstXSize, nDstYSize,
1983 : &nSrcXOff, &nSrcYOff, &nSrcXSize, &nSrcYSize,
1984 : &dfSrcXExtraSize, &dfSrcYExtraSize, nullptr);
1985 1777 : if (hWarpMutex != nullptr)
1986 0 : CPLReleaseMutex(hWarpMutex);
1987 1777 : if (eErr != CE_None)
1988 : {
1989 : const bool bErrorOutIfEmptySourceWindow =
1990 36 : CPLFetchBool(psOptions->papszWarpOptions,
1991 : "ERROR_OUT_IF_EMPTY_SOURCE_WINDOW", true);
1992 36 : if (!bErrorOutIfEmptySourceWindow)
1993 36 : return CE_None;
1994 0 : return eErr;
1995 : }
1996 : }
1997 :
1998 : /* -------------------------------------------------------------------- */
1999 : /* Prepare a WarpKernel object to match this operation. */
2000 : /* -------------------------------------------------------------------- */
2001 7718 : GDALWarpKernel oWK;
2002 :
2003 3859 : oWK.eResample = m_bIsTranslationOnPixelBoundaries ? GRA_NearestNeighbour
2004 3544 : : psOptions->eResampleAlg;
2005 3859 : oWK.eTieStrategy = psOptions->eTieStrategy;
2006 3859 : oWK.nBands = psOptions->nBandCount;
2007 3859 : oWK.eWorkingDataType = psOptions->eWorkingDataType;
2008 :
2009 3859 : oWK.pfnTransformer = psOptions->pfnTransformer;
2010 3859 : oWK.pTransformerArg = psOptions->pTransformerArg;
2011 :
2012 3859 : oWK.pfnProgress = psOptions->pfnProgress;
2013 3859 : oWK.pProgress = psOptions->pProgressArg;
2014 3859 : oWK.dfProgressBase = dfProgressBase;
2015 3859 : oWK.dfProgressScale = dfProgressScale;
2016 :
2017 3859 : oWK.papszWarpOptions = psOptions->papszWarpOptions;
2018 3859 : oWK.psThreadData = psThreadData;
2019 :
2020 3859 : oWK.padfDstNoDataReal = psOptions->padfDstNoDataReal;
2021 :
2022 : /* -------------------------------------------------------------------- */
2023 : /* Setup the source buffer. */
2024 : /* */
2025 : /* Eventually we may need to take advantage of pixel */
2026 : /* interleaved reading here. */
2027 : /* -------------------------------------------------------------------- */
2028 3859 : oWK.nSrcXOff = nSrcXOff;
2029 3859 : oWK.nSrcYOff = nSrcYOff;
2030 3859 : oWK.nSrcXSize = nSrcXSize;
2031 3859 : oWK.nSrcYSize = nSrcYSize;
2032 3859 : oWK.dfSrcXExtraSize = dfSrcXExtraSize;
2033 3859 : oWK.dfSrcYExtraSize = dfSrcYExtraSize;
2034 :
2035 : // Check for overflows in computation of nAlloc
2036 7115 : if (nSrcYSize > 0 &&
2037 3256 : ((static_cast<size_t>(nSrcXSize) >
2038 3256 : (std::numeric_limits<size_t>::max() - WARP_EXTRA_ELTS) / nSrcYSize) ||
2039 3256 : (static_cast<size_t>(nSrcXSize) * nSrcYSize + WARP_EXTRA_ELTS >
2040 3256 : std::numeric_limits<size_t>::max() /
2041 3256 : (nWordSize * psOptions->nBandCount))))
2042 : {
2043 0 : CPLError(CE_Failure, CPLE_AppDefined,
2044 : "WarpRegionToBuffer(): Integer overflow : nWordSize(=%d) * "
2045 : "(nSrcXSize(=%d) * nSrcYSize(=%d) + WARP_EXTRA_ELTS(=%d)) * "
2046 : "nBandCount(=%d)",
2047 : nWordSize, nSrcXSize, nSrcYSize, WARP_EXTRA_ELTS,
2048 0 : psOptions->nBandCount);
2049 0 : return CE_Failure;
2050 : }
2051 :
2052 3859 : const size_t nAlloc =
2053 3859 : nWordSize *
2054 3859 : (static_cast<size_t>(nSrcXSize) * nSrcYSize + WARP_EXTRA_ELTS) *
2055 3859 : psOptions->nBandCount;
2056 :
2057 3859 : oWK.papabySrcImage = static_cast<GByte **>(
2058 3859 : CPLCalloc(sizeof(GByte *), psOptions->nBandCount));
2059 3859 : oWK.papabySrcImage[0] = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nAlloc));
2060 :
2061 3859 : CPLErr eErr =
2062 3257 : nSrcXSize != 0 && nSrcYSize != 0 && oWK.papabySrcImage[0] == nullptr
2063 7116 : ? CE_Failure
2064 : : CE_None;
2065 :
2066 11309 : for (int i = 0; i < psOptions->nBandCount && eErr == CE_None; i++)
2067 7450 : oWK.papabySrcImage[i] =
2068 7450 : reinterpret_cast<GByte *>(oWK.papabySrcImage[0]) +
2069 7450 : nWordSize *
2070 7450 : (static_cast<GPtrDiff_t>(nSrcXSize) * nSrcYSize +
2071 7450 : WARP_EXTRA_ELTS) *
2072 7450 : i;
2073 :
2074 3859 : if (eErr == CE_None && nSrcXSize > 0 && nSrcYSize > 0)
2075 : {
2076 3246 : GDALDataset *poSrcDS = GDALDataset::FromHandle(psOptions->hSrcDS);
2077 3246 : if (psOptions->nBandCount == 1)
2078 : {
2079 : // Particular case to simplify the stack a bit.
2080 2023 : eErr = poSrcDS->GetRasterBand(psOptions->panSrcBands[0])
2081 4046 : ->RasterIO(GF_Read, nSrcXOff, nSrcYOff, nSrcXSize,
2082 2023 : nSrcYSize, oWK.papabySrcImage[0], nSrcXSize,
2083 2023 : nSrcYSize, psOptions->eWorkingDataType, 0, 0,
2084 : nullptr);
2085 : }
2086 : else
2087 : {
2088 1223 : eErr = poSrcDS->RasterIO(
2089 : GF_Read, nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
2090 1223 : oWK.papabySrcImage[0], nSrcXSize, nSrcYSize,
2091 1223 : psOptions->eWorkingDataType, psOptions->nBandCount,
2092 1223 : psOptions->panSrcBands, 0, 0,
2093 1223 : nWordSize * (static_cast<GPtrDiff_t>(nSrcXSize) * nSrcYSize +
2094 : WARP_EXTRA_ELTS),
2095 : nullptr);
2096 : }
2097 : }
2098 :
2099 3859 : ReportTiming("Input buffer read");
2100 :
2101 : /* -------------------------------------------------------------------- */
2102 : /* Initialize destination buffer. */
2103 : /* -------------------------------------------------------------------- */
2104 3859 : oWK.nDstXOff = nDstXOff;
2105 3859 : oWK.nDstYOff = nDstYOff;
2106 3859 : oWK.nDstXSize = nDstXSize;
2107 3859 : oWK.nDstYSize = nDstYSize;
2108 :
2109 3859 : oWK.papabyDstImage = reinterpret_cast<GByte **>(
2110 3859 : CPLCalloc(sizeof(GByte *), psOptions->nBandCount));
2111 :
2112 11301 : for (int i = 0; i < psOptions->nBandCount && eErr == CE_None; i++)
2113 : {
2114 7442 : oWK.papabyDstImage[i] =
2115 7442 : static_cast<GByte *>(pDataBuf) +
2116 7442 : i * static_cast<GPtrDiff_t>(nDstXSize) * nDstYSize * nWordSize;
2117 : }
2118 :
2119 : /* -------------------------------------------------------------------- */
2120 : /* Eventually we need handling for a whole bunch of the */
2121 : /* validity and density masks here. */
2122 : /* -------------------------------------------------------------------- */
2123 :
2124 : // TODO
2125 :
2126 : /* -------------------------------------------------------------------- */
2127 : /* Generate a source density mask if we have a source alpha band */
2128 : /* -------------------------------------------------------------------- */
2129 3859 : if (eErr == CE_None && psOptions->nSrcAlphaBand > 0 && nSrcXSize > 0 &&
2130 142 : nSrcYSize > 0)
2131 : {
2132 142 : CPLAssert(oWK.pafUnifiedSrcDensity == nullptr);
2133 :
2134 142 : eErr = CreateKernelMask(&oWK, 0 /* not used */, "UnifiedSrcDensity");
2135 :
2136 142 : if (eErr == CE_None)
2137 : {
2138 142 : int bOutAllOpaque = FALSE;
2139 284 : eErr = GDALWarpSrcAlphaMasker(
2140 142 : psOptions, psOptions->nBandCount, psOptions->eWorkingDataType,
2141 : oWK.nSrcXOff, oWK.nSrcYOff, oWK.nSrcXSize, oWK.nSrcYSize,
2142 142 : oWK.papabySrcImage, TRUE, oWK.pafUnifiedSrcDensity,
2143 : &bOutAllOpaque);
2144 142 : if (bOutAllOpaque)
2145 : {
2146 : #if DEBUG_VERBOSE
2147 : CPLDebug("WARP",
2148 : "No need for a source density mask as all values "
2149 : "are opaque");
2150 : #endif
2151 35 : CPLFree(oWK.pafUnifiedSrcDensity);
2152 35 : oWK.pafUnifiedSrcDensity = nullptr;
2153 : }
2154 : }
2155 : }
2156 :
2157 : /* -------------------------------------------------------------------- */
2158 : /* Generate a source density mask if we have a source cutline. */
2159 : /* -------------------------------------------------------------------- */
2160 3859 : if (eErr == CE_None && psOptions->hCutline != nullptr && nSrcXSize > 0 &&
2161 44 : nSrcYSize > 0)
2162 : {
2163 44 : const bool bUnifiedSrcDensityJustCreated =
2164 44 : (oWK.pafUnifiedSrcDensity == nullptr);
2165 44 : if (bUnifiedSrcDensityJustCreated)
2166 : {
2167 : eErr =
2168 44 : CreateKernelMask(&oWK, 0 /* not used */, "UnifiedSrcDensity");
2169 :
2170 44 : if (eErr == CE_None)
2171 : {
2172 44 : for (GPtrDiff_t j = 0;
2173 2929360 : j < static_cast<GPtrDiff_t>(oWK.nSrcXSize) * oWK.nSrcYSize;
2174 : j++)
2175 2929320 : oWK.pafUnifiedSrcDensity[j] = 1.0;
2176 : }
2177 : }
2178 :
2179 44 : int nValidityFlag = 0;
2180 44 : if (eErr == CE_None)
2181 44 : eErr = GDALWarpCutlineMaskerEx(
2182 44 : psOptions, psOptions->nBandCount, psOptions->eWorkingDataType,
2183 : oWK.nSrcXOff, oWK.nSrcYOff, oWK.nSrcXSize, oWK.nSrcYSize,
2184 44 : oWK.papabySrcImage, TRUE, oWK.pafUnifiedSrcDensity,
2185 : &nValidityFlag);
2186 44 : if (nValidityFlag == GCMVF_CHUNK_FULLY_WITHIN_CUTLINE &&
2187 : bUnifiedSrcDensityJustCreated)
2188 : {
2189 8 : VSIFree(oWK.pafUnifiedSrcDensity);
2190 8 : oWK.pafUnifiedSrcDensity = nullptr;
2191 : }
2192 : }
2193 :
2194 : /* -------------------------------------------------------------------- */
2195 : /* Generate a destination density mask if we have a destination */
2196 : /* alpha band. */
2197 : /* -------------------------------------------------------------------- */
2198 3859 : if (eErr == CE_None && psOptions->nDstAlphaBand > 0)
2199 : {
2200 1748 : CPLAssert(oWK.pafDstDensity == nullptr);
2201 :
2202 1748 : eErr = CreateKernelMask(&oWK, 0 /* not used */, "DstDensity");
2203 :
2204 1748 : if (eErr == CE_None)
2205 1748 : eErr = GDALWarpDstAlphaMasker(
2206 1748 : psOptions, psOptions->nBandCount, psOptions->eWorkingDataType,
2207 : oWK.nDstXOff, oWK.nDstYOff, oWK.nDstXSize, oWK.nDstYSize,
2208 1748 : oWK.papabyDstImage, TRUE, oWK.pafDstDensity);
2209 : }
2210 :
2211 : /* -------------------------------------------------------------------- */
2212 : /* If we have source nodata values create the validity mask. */
2213 : /* -------------------------------------------------------------------- */
2214 3859 : if (eErr == CE_None && psOptions->padfSrcNoDataReal != nullptr &&
2215 232 : nSrcXSize > 0 && nSrcYSize > 0)
2216 : {
2217 214 : CPLAssert(oWK.papanBandSrcValid == nullptr);
2218 :
2219 214 : bool bAllBandsAllValid = true;
2220 539 : for (int i = 0; i < psOptions->nBandCount && eErr == CE_None; i++)
2221 : {
2222 325 : eErr = CreateKernelMask(&oWK, i, "BandSrcValid");
2223 325 : if (eErr == CE_None)
2224 : {
2225 325 : double adfNoData[2] = {psOptions->padfSrcNoDataReal[i],
2226 325 : psOptions->padfSrcNoDataImag != nullptr
2227 325 : ? psOptions->padfSrcNoDataImag[i]
2228 325 : : 0.0};
2229 :
2230 325 : int bAllValid = FALSE;
2231 650 : eErr = GDALWarpNoDataMasker(
2232 325 : adfNoData, 1, psOptions->eWorkingDataType, oWK.nSrcXOff,
2233 : oWK.nSrcYOff, oWK.nSrcXSize, oWK.nSrcYSize,
2234 325 : &(oWK.papabySrcImage[i]), FALSE, oWK.papanBandSrcValid[i],
2235 : &bAllValid);
2236 325 : if (!bAllValid)
2237 188 : bAllBandsAllValid = false;
2238 : }
2239 : }
2240 :
2241 : // Optimization: if all pixels in all bands are valid,
2242 : // we don't need a mask.
2243 214 : if (bAllBandsAllValid)
2244 : {
2245 : #if DEBUG_VERBOSE
2246 : CPLDebug(
2247 : "WARP",
2248 : "No need for a source nodata mask as all values are valid");
2249 : #endif
2250 231 : for (int k = 0; k < psOptions->nBandCount; k++)
2251 136 : CPLFree(oWK.papanBandSrcValid[k]);
2252 95 : CPLFree(oWK.papanBandSrcValid);
2253 95 : oWK.papanBandSrcValid = nullptr;
2254 : }
2255 :
2256 : /* --------------------------------------------------------------------
2257 : */
2258 : /* If there's just a single band, then transfer */
2259 : /* papanBandSrcValid[0] as panUnifiedSrcValid. */
2260 : /* --------------------------------------------------------------------
2261 : */
2262 214 : if (oWK.papanBandSrcValid != nullptr && psOptions->nBandCount == 1)
2263 : {
2264 73 : oWK.panUnifiedSrcValid = oWK.papanBandSrcValid[0];
2265 73 : CPLFree(oWK.papanBandSrcValid);
2266 73 : oWK.papanBandSrcValid = nullptr;
2267 : }
2268 :
2269 : /* --------------------------------------------------------------------
2270 : */
2271 : /* Compute a unified input pixel mask if and only if all bands */
2272 : /* nodata is true. That is, we only treat a pixel as nodata if */
2273 : /* all bands match their respective nodata values. */
2274 : /* --------------------------------------------------------------------
2275 : */
2276 141 : else if (oWK.papanBandSrcValid != nullptr && eErr == CE_None)
2277 : {
2278 46 : bool bAtLeastOneBandAllValid = false;
2279 162 : for (int k = 0; k < psOptions->nBandCount; k++)
2280 : {
2281 116 : if (oWK.papanBandSrcValid[k] == nullptr)
2282 : {
2283 0 : bAtLeastOneBandAllValid = true;
2284 0 : break;
2285 : }
2286 : }
2287 :
2288 92 : const char *pszUnifiedSrcNoData = CSLFetchNameValue(
2289 46 : psOptions->papszWarpOptions, "UNIFIED_SRC_NODATA");
2290 84 : if (!bAtLeastOneBandAllValid && (pszUnifiedSrcNoData == nullptr ||
2291 38 : CPLTestBool(pszUnifiedSrcNoData)))
2292 : {
2293 45 : auto nMaskBits =
2294 45 : static_cast<GPtrDiff_t>(oWK.nSrcXSize) * oWK.nSrcYSize;
2295 :
2296 : eErr =
2297 45 : CreateKernelMask(&oWK, 0 /* not used */, "UnifiedSrcValid");
2298 :
2299 45 : if (eErr == CE_None)
2300 : {
2301 45 : CPLMaskClearAll(oWK.panUnifiedSrcValid, nMaskBits);
2302 :
2303 158 : for (int k = 0; k < psOptions->nBandCount; k++)
2304 : {
2305 113 : CPLMaskMerge(oWK.panUnifiedSrcValid,
2306 113 : oWK.papanBandSrcValid[k], nMaskBits);
2307 : }
2308 :
2309 : // If UNIFIED_SRC_NODATA is set, then we will ignore the
2310 : // individual nodata status of each band. If it is not set,
2311 : // both mechanism apply:
2312 : // - if panUnifiedSrcValid[] indicates a pixel is invalid
2313 : // (that is all its bands are at nodata), then the output
2314 : // pixel will be invalid
2315 : // - otherwise, the status band per band will be check with
2316 : // papanBandSrcValid[iBand][], and the output pixel will
2317 : // be valid
2318 45 : if (pszUnifiedSrcNoData != nullptr &&
2319 37 : !EQUAL(pszUnifiedSrcNoData, "PARTIAL"))
2320 : {
2321 123 : for (int k = 0; k < psOptions->nBandCount; k++)
2322 87 : CPLFree(oWK.papanBandSrcValid[k]);
2323 36 : CPLFree(oWK.papanBandSrcValid);
2324 36 : oWK.papanBandSrcValid = nullptr;
2325 : }
2326 : }
2327 : }
2328 : }
2329 : }
2330 :
2331 : /* -------------------------------------------------------------------- */
2332 : /* Generate a source validity mask if we have a source mask for */
2333 : /* the whole input dataset (and didn't already treat it as */
2334 : /* alpha band). */
2335 : /* -------------------------------------------------------------------- */
2336 : GDALRasterBandH hSrcBand =
2337 3859 : psOptions->nBandCount < 1
2338 3859 : ? nullptr
2339 3859 : : GDALGetRasterBand(psOptions->hSrcDS, psOptions->panSrcBands[0]);
2340 :
2341 3855 : if (eErr == CE_None && oWK.pafUnifiedSrcDensity == nullptr &&
2342 3712 : oWK.panUnifiedSrcValid == nullptr && psOptions->nSrcAlphaBand <= 0 &&
2343 3531 : (GDALGetMaskFlags(hSrcBand) & GMF_PER_DATASET)
2344 : // Need to double check for -nosrcalpha case.
2345 5 : && !(GDALGetMaskFlags(hSrcBand) & GMF_ALPHA) &&
2346 7716 : psOptions->padfSrcNoDataReal == nullptr && nSrcXSize > 0 &&
2347 2 : nSrcYSize > 0)
2348 :
2349 : {
2350 2 : eErr = CreateKernelMask(&oWK, 0 /* not used */, "UnifiedSrcValid");
2351 :
2352 2 : if (eErr == CE_None)
2353 2 : eErr = GDALWarpSrcMaskMasker(
2354 2 : psOptions, psOptions->nBandCount, psOptions->eWorkingDataType,
2355 : oWK.nSrcXOff, oWK.nSrcYOff, oWK.nSrcXSize, oWK.nSrcYSize,
2356 2 : oWK.papabySrcImage, FALSE, oWK.panUnifiedSrcValid);
2357 : }
2358 :
2359 : /* -------------------------------------------------------------------- */
2360 : /* If we have destination nodata values create the */
2361 : /* validity mask. We set the DstValid for any pixel that we */
2362 : /* do no have valid data in *any* of the source bands. */
2363 : /* */
2364 : /* Note that we don't support any concept of unified nodata on */
2365 : /* the destination image. At some point that should be added */
2366 : /* and then this logic will be significantly different. */
2367 : /* -------------------------------------------------------------------- */
2368 3859 : if (eErr == CE_None && psOptions->padfDstNoDataReal != nullptr)
2369 : {
2370 587 : CPLAssert(oWK.panDstValid == nullptr);
2371 :
2372 587 : const GPtrDiff_t nMaskBits =
2373 587 : static_cast<GPtrDiff_t>(oWK.nDstXSize) * oWK.nDstYSize;
2374 :
2375 587 : eErr = CreateKernelMask(&oWK, 0 /* not used */, "DstValid");
2376 : GUInt32 *panBandMask =
2377 587 : eErr == CE_None ? CPLMaskCreate(nMaskBits, true) : nullptr;
2378 :
2379 587 : if (eErr == CE_None && panBandMask != nullptr)
2380 : {
2381 1241 : for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
2382 : {
2383 712 : CPLMaskSetAll(panBandMask, nMaskBits);
2384 :
2385 712 : double adfNoData[2] = {psOptions->padfDstNoDataReal[iBand],
2386 712 : psOptions->padfDstNoDataImag != nullptr
2387 712 : ? psOptions->padfDstNoDataImag[iBand]
2388 712 : : 0.0};
2389 :
2390 712 : int bAllValid = FALSE;
2391 1424 : eErr = GDALWarpNoDataMasker(
2392 712 : adfNoData, 1, psOptions->eWorkingDataType, oWK.nDstXOff,
2393 : oWK.nDstYOff, oWK.nDstXSize, oWK.nDstYSize,
2394 712 : oWK.papabyDstImage + iBand, FALSE, panBandMask, &bAllValid);
2395 :
2396 : // Optimization: if there's a single band and all pixels are
2397 : // valid then we don't need a mask.
2398 712 : if (bAllValid && psOptions->nBandCount == 1)
2399 : {
2400 : #if DEBUG_VERBOSE
2401 : CPLDebug("WARP", "No need for a destination nodata mask as "
2402 : "all values are valid");
2403 : #endif
2404 58 : CPLFree(oWK.panDstValid);
2405 58 : oWK.panDstValid = nullptr;
2406 58 : break;
2407 : }
2408 :
2409 654 : CPLMaskMerge(oWK.panDstValid, panBandMask, nMaskBits);
2410 : }
2411 587 : CPLFree(panBandMask);
2412 : }
2413 : }
2414 :
2415 : /* -------------------------------------------------------------------- */
2416 : /* Release IO Mutex, and acquire warper mutex. */
2417 : /* -------------------------------------------------------------------- */
2418 3859 : if (hIOMutex != nullptr)
2419 : {
2420 11 : CPLReleaseMutex(hIOMutex);
2421 11 : if (!CPLAcquireMutex(hWarpMutex, 600.0))
2422 : {
2423 0 : CPLError(CE_Failure, CPLE_AppDefined,
2424 : "Failed to acquire WarpMutex in WarpRegion().");
2425 0 : return CE_Failure;
2426 : }
2427 : }
2428 :
2429 : /* -------------------------------------------------------------------- */
2430 : /* Optional application provided prewarp chunk processor. */
2431 : /* -------------------------------------------------------------------- */
2432 3859 : if (eErr == CE_None && psOptions->pfnPreWarpChunkProcessor != nullptr)
2433 0 : eErr = psOptions->pfnPreWarpChunkProcessor(
2434 0 : &oWK, psOptions->pPreWarpProcessorArg);
2435 :
2436 : /* -------------------------------------------------------------------- */
2437 : /* Perform the warp. */
2438 : /* -------------------------------------------------------------------- */
2439 3859 : if (eErr == CE_None)
2440 : {
2441 3855 : eErr = oWK.PerformWarp();
2442 3855 : ReportTiming("In memory warp operation");
2443 : }
2444 :
2445 : /* -------------------------------------------------------------------- */
2446 : /* Optional application provided postwarp chunk processor. */
2447 : /* -------------------------------------------------------------------- */
2448 3859 : if (eErr == CE_None && psOptions->pfnPostWarpChunkProcessor != nullptr)
2449 0 : eErr = psOptions->pfnPostWarpChunkProcessor(
2450 0 : &oWK, psOptions->pPostWarpProcessorArg);
2451 :
2452 : /* -------------------------------------------------------------------- */
2453 : /* Release Warp Mutex, and acquire io mutex. */
2454 : /* -------------------------------------------------------------------- */
2455 3859 : if (hIOMutex != nullptr)
2456 : {
2457 11 : CPLReleaseMutex(hWarpMutex);
2458 11 : if (!CPLAcquireMutex(hIOMutex, 600.0))
2459 : {
2460 0 : CPLError(CE_Failure, CPLE_AppDefined,
2461 : "Failed to acquire IOMutex in WarpRegion().");
2462 0 : return CE_Failure;
2463 : }
2464 : }
2465 :
2466 : /* -------------------------------------------------------------------- */
2467 : /* Write destination alpha if available. */
2468 : /* -------------------------------------------------------------------- */
2469 3859 : if (eErr == CE_None && psOptions->nDstAlphaBand > 0)
2470 : {
2471 1748 : eErr = GDALWarpDstAlphaMasker(
2472 1748 : psOptions, -psOptions->nBandCount, psOptions->eWorkingDataType,
2473 : oWK.nDstXOff, oWK.nDstYOff, oWK.nDstXSize, oWK.nDstYSize,
2474 1748 : oWK.papabyDstImage, TRUE, oWK.pafDstDensity);
2475 : }
2476 :
2477 : /* -------------------------------------------------------------------- */
2478 : /* Cleanup. */
2479 : /* -------------------------------------------------------------------- */
2480 3859 : CPLFree(oWK.papabySrcImage[0]);
2481 3859 : CPLFree(oWK.papabySrcImage);
2482 3859 : CPLFree(oWK.papabyDstImage);
2483 :
2484 3859 : if (oWK.papanBandSrcValid != nullptr)
2485 : {
2486 39 : for (int i = 0; i < oWK.nBands; i++)
2487 29 : CPLFree(oWK.papanBandSrcValid[i]);
2488 10 : CPLFree(oWK.papanBandSrcValid);
2489 : }
2490 3859 : CPLFree(oWK.panUnifiedSrcValid);
2491 3859 : CPLFree(oWK.pafUnifiedSrcDensity);
2492 3859 : CPLFree(oWK.panDstValid);
2493 3859 : CPLFree(oWK.pafDstDensity);
2494 :
2495 3859 : return eErr;
2496 : }
2497 :
2498 : /************************************************************************/
2499 : /* GDALWarpRegionToBuffer() */
2500 : /************************************************************************/
2501 :
2502 : /**
2503 : * @see GDALWarpOperation::WarpRegionToBuffer()
2504 : */
2505 :
2506 0 : CPLErr GDALWarpRegionToBuffer(GDALWarpOperationH hOperation, int nDstXOff,
2507 : int nDstYOff, int nDstXSize, int nDstYSize,
2508 : void *pDataBuf, GDALDataType eBufDataType,
2509 : int nSrcXOff, int nSrcYOff, int nSrcXSize,
2510 : int nSrcYSize)
2511 :
2512 : {
2513 0 : VALIDATE_POINTER1(hOperation, "GDALWarpRegionToBuffer", CE_Failure);
2514 :
2515 : return reinterpret_cast<GDALWarpOperation *>(hOperation)
2516 0 : ->WarpRegionToBuffer(nDstXOff, nDstYOff, nDstXSize, nDstYSize, pDataBuf,
2517 : eBufDataType, nSrcXOff, nSrcYOff, nSrcXSize,
2518 0 : nSrcYSize);
2519 : }
2520 :
2521 : /************************************************************************/
2522 : /* CreateKernelMask() */
2523 : /* */
2524 : /* If mask does not yet exist, create it. Supported types are */
2525 : /* the name of the variable in question. That is */
2526 : /* "BandSrcValid", "UnifiedSrcValid", "UnifiedSrcDensity", */
2527 : /* "DstValid", and "DstDensity". */
2528 : /************************************************************************/
2529 :
2530 2893 : CPLErr GDALWarpOperation::CreateKernelMask(GDALWarpKernel *poKernel, int iBand,
2531 : const char *pszType)
2532 :
2533 : {
2534 2893 : void **ppMask = nullptr;
2535 2893 : int nXSize = 0;
2536 2893 : int nYSize = 0;
2537 2893 : int nBitsPerPixel = 0;
2538 2893 : int nDefault = 0;
2539 2893 : int nExtraElts = 0;
2540 2893 : bool bDoMemset = true;
2541 :
2542 : /* -------------------------------------------------------------------- */
2543 : /* Get particulars of mask to be updated. */
2544 : /* -------------------------------------------------------------------- */
2545 2893 : if (EQUAL(pszType, "BandSrcValid"))
2546 : {
2547 325 : if (poKernel->papanBandSrcValid == nullptr)
2548 214 : poKernel->papanBandSrcValid = static_cast<GUInt32 **>(
2549 214 : CPLCalloc(sizeof(void *), poKernel->nBands));
2550 :
2551 325 : ppMask =
2552 325 : reinterpret_cast<void **>(&(poKernel->papanBandSrcValid[iBand]));
2553 325 : nExtraElts = WARP_EXTRA_ELTS;
2554 325 : nXSize = poKernel->nSrcXSize;
2555 325 : nYSize = poKernel->nSrcYSize;
2556 325 : nBitsPerPixel = 1;
2557 325 : nDefault = 0xff;
2558 : }
2559 2568 : else if (EQUAL(pszType, "UnifiedSrcValid"))
2560 : {
2561 47 : ppMask = reinterpret_cast<void **>(&(poKernel->panUnifiedSrcValid));
2562 47 : nExtraElts = WARP_EXTRA_ELTS;
2563 47 : nXSize = poKernel->nSrcXSize;
2564 47 : nYSize = poKernel->nSrcYSize;
2565 47 : nBitsPerPixel = 1;
2566 47 : nDefault = 0xff;
2567 : }
2568 2521 : else if (EQUAL(pszType, "UnifiedSrcDensity"))
2569 : {
2570 186 : ppMask = reinterpret_cast<void **>(&(poKernel->pafUnifiedSrcDensity));
2571 186 : nExtraElts = WARP_EXTRA_ELTS;
2572 186 : nXSize = poKernel->nSrcXSize;
2573 186 : nYSize = poKernel->nSrcYSize;
2574 186 : nBitsPerPixel = 32;
2575 186 : nDefault = 0;
2576 186 : bDoMemset = false;
2577 : }
2578 2335 : else if (EQUAL(pszType, "DstValid"))
2579 : {
2580 587 : ppMask = reinterpret_cast<void **>(&(poKernel->panDstValid));
2581 587 : nXSize = poKernel->nDstXSize;
2582 587 : nYSize = poKernel->nDstYSize;
2583 587 : nBitsPerPixel = 1;
2584 587 : nDefault = 0;
2585 : }
2586 1748 : else if (EQUAL(pszType, "DstDensity"))
2587 : {
2588 1748 : ppMask = reinterpret_cast<void **>(&(poKernel->pafDstDensity));
2589 1748 : nXSize = poKernel->nDstXSize;
2590 1748 : nYSize = poKernel->nDstYSize;
2591 1748 : nBitsPerPixel = 32;
2592 1748 : nDefault = 0;
2593 1748 : bDoMemset = false;
2594 : }
2595 : else
2596 : {
2597 0 : CPLError(CE_Failure, CPLE_AppDefined,
2598 : "Internal error in CreateKernelMask(%s).", pszType);
2599 0 : return CE_Failure;
2600 : }
2601 :
2602 : /* -------------------------------------------------------------------- */
2603 : /* Allocate if needed. */
2604 : /* -------------------------------------------------------------------- */
2605 2893 : if (*ppMask == nullptr)
2606 : {
2607 2893 : const GIntBig nBytes =
2608 : nBitsPerPixel == 32
2609 2893 : ? (static_cast<GIntBig>(nXSize) * nYSize + nExtraElts) * 4
2610 959 : : (static_cast<GIntBig>(nXSize) * nYSize + nExtraElts + 31) / 8;
2611 :
2612 2893 : const size_t nByteSize_t = static_cast<size_t>(nBytes);
2613 : #if SIZEOF_VOIDP == 4
2614 : if (static_cast<GIntBig>(nByteSize_t) != nBytes)
2615 : {
2616 : CPLError(CE_Failure, CPLE_OutOfMemory,
2617 : "Cannot allocate " CPL_FRMT_GIB " bytes", nBytes);
2618 : return CE_Failure;
2619 : }
2620 : #endif
2621 :
2622 2893 : *ppMask = VSI_MALLOC_VERBOSE(nByteSize_t);
2623 :
2624 2893 : if (*ppMask == nullptr)
2625 : {
2626 0 : return CE_Failure;
2627 : }
2628 :
2629 2893 : if (bDoMemset)
2630 959 : memset(*ppMask, nDefault, nByteSize_t);
2631 : }
2632 :
2633 2893 : return CE_None;
2634 : }
2635 :
2636 : /************************************************************************/
2637 : /* ComputeSourceWindowStartingFromSource() */
2638 : /************************************************************************/
2639 :
2640 : constexpr int DEFAULT_STEP_COUNT = 21;
2641 :
2642 1260 : void GDALWarpOperation::ComputeSourceWindowStartingFromSource(
2643 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize,
2644 : double *padfSrcMinX, double *padfSrcMinY, double *padfSrcMaxX,
2645 : double *padfSrcMaxY)
2646 : {
2647 1260 : const int nSrcRasterXSize = GDALGetRasterXSize(psOptions->hSrcDS);
2648 1260 : const int nSrcRasterYSize = GDALGetRasterYSize(psOptions->hSrcDS);
2649 1260 : if (nSrcRasterXSize == 0 || nSrcRasterYSize == 0)
2650 0 : return;
2651 :
2652 1260 : GDALWarpPrivateData *privateData = GetWarpPrivateData(this);
2653 1260 : if (privateData->nStepCount == 0)
2654 : {
2655 312 : int nStepCount = DEFAULT_STEP_COUNT;
2656 312 : std::vector<double> adfDstZ{};
2657 :
2658 : const char *pszSampleSteps =
2659 312 : CSLFetchNameValue(psOptions->papszWarpOptions, "SAMPLE_STEPS");
2660 312 : constexpr int knIntMax = std::numeric_limits<int>::max();
2661 312 : if (pszSampleSteps && !EQUAL(pszSampleSteps, "ALL"))
2662 : {
2663 0 : nStepCount = atoi(
2664 0 : CSLFetchNameValue(psOptions->papszWarpOptions, "SAMPLE_STEPS"));
2665 0 : nStepCount = std::max(2, nStepCount);
2666 : }
2667 :
2668 312 : const double dfStepSize = 1.0 / (nStepCount - 1);
2669 312 : if (nStepCount > knIntMax - 2 ||
2670 312 : (nStepCount + 2) > knIntMax / (nStepCount + 2))
2671 : {
2672 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps : %d",
2673 : nStepCount);
2674 0 : return;
2675 : }
2676 312 : const int nSampleMax = (nStepCount + 2) * (nStepCount + 2);
2677 :
2678 : try
2679 : {
2680 312 : privateData->abSuccess.resize(nSampleMax);
2681 312 : privateData->adfDstX.resize(nSampleMax);
2682 312 : privateData->adfDstY.resize(nSampleMax);
2683 312 : adfDstZ.resize(nSampleMax);
2684 : }
2685 0 : catch (const std::exception &)
2686 : {
2687 0 : return;
2688 : }
2689 :
2690 : /* --------------------------------------------------------------------
2691 : */
2692 : /* Setup sample points on a grid pattern throughout the source */
2693 : /* raster. */
2694 : /* --------------------------------------------------------------------
2695 : */
2696 312 : int iPoint = 0;
2697 7488 : for (int iY = 0; iY < nStepCount + 2; iY++)
2698 : {
2699 14040 : const double dfRatioY = (iY == 0) ? 0.5 / nSrcRasterYSize
2700 6864 : : (iY <= nStepCount)
2701 6864 : ? (iY - 1) * dfStepSize
2702 312 : : 1 - 0.5 / nSrcRasterYSize;
2703 172224 : for (int iX = 0; iX < nStepCount + 2; iX++)
2704 : {
2705 322920 : const double dfRatioX = (iX == 0) ? 0.5 / nSrcRasterXSize
2706 157872 : : (iX <= nStepCount)
2707 157872 : ? (iX - 1) * dfStepSize
2708 7176 : : 1 - 0.5 / nSrcRasterXSize;
2709 165048 : privateData->adfDstX[iPoint] = dfRatioX * nSrcRasterXSize;
2710 165048 : privateData->adfDstY[iPoint] = dfRatioY * nSrcRasterYSize;
2711 165048 : iPoint++;
2712 : }
2713 : }
2714 312 : CPLAssert(iPoint == nSampleMax);
2715 :
2716 : /* --------------------------------------------------------------------
2717 : */
2718 : /* Transform them to the output pixel coordinate space */
2719 : /* --------------------------------------------------------------------
2720 : */
2721 312 : psOptions->pfnTransformer(psOptions->pTransformerArg, FALSE, nSampleMax,
2722 : privateData->adfDstX.data(),
2723 : privateData->adfDstY.data(), adfDstZ.data(),
2724 : privateData->abSuccess.data());
2725 312 : privateData->nStepCount = nStepCount;
2726 : }
2727 :
2728 : /* -------------------------------------------------------------------- */
2729 : /* Collect the bounds, ignoring any failed points. */
2730 : /* -------------------------------------------------------------------- */
2731 1260 : const int nStepCount = privateData->nStepCount;
2732 1260 : const double dfStepSize = 1.0 / (nStepCount - 1);
2733 1260 : int iPoint = 0;
2734 : #ifdef DEBUG
2735 1260 : const size_t nSampleMax =
2736 1260 : static_cast<size_t>(nStepCount + 2) * (nStepCount + 2);
2737 1260 : CPL_IGNORE_RET_VAL(nSampleMax);
2738 1260 : CPLAssert(privateData->adfDstX.size() == nSampleMax);
2739 1260 : CPLAssert(privateData->adfDstY.size() == nSampleMax);
2740 1260 : CPLAssert(privateData->abSuccess.size() == nSampleMax);
2741 : #endif
2742 30240 : for (int iY = 0; iY < nStepCount + 2; iY++)
2743 : {
2744 56700 : const double dfRatioY = (iY == 0) ? 0.5 / nSrcRasterYSize
2745 : : (iY <= nStepCount)
2746 27720 : ? (iY - 1) * dfStepSize
2747 1260 : : 1 - 0.5 / nSrcRasterYSize;
2748 695520 : for (int iX = 0; iX < nStepCount + 2; iX++)
2749 : {
2750 666540 : if (privateData->abSuccess[iPoint] &&
2751 631068 : privateData->adfDstX[iPoint] >= nDstXOff &&
2752 458825 : privateData->adfDstX[iPoint] <= nDstXOff + nDstXSize &&
2753 1528080 : privateData->adfDstY[iPoint] >= nDstYOff &&
2754 230469 : privateData->adfDstY[iPoint] <= nDstYOff + nDstYSize)
2755 : {
2756 343563 : const double dfRatioX = (iX == 0) ? 0.5 / nSrcRasterXSize
2757 : : (iX <= nStepCount)
2758 168009 : ? (iX - 1) * dfStepSize
2759 7543 : : 1 - 0.5 / nSrcRasterXSize;
2760 175554 : double dfSrcX = dfRatioX * nSrcRasterXSize;
2761 175554 : double dfSrcY = dfRatioY * nSrcRasterYSize;
2762 175554 : *padfSrcMinX = std::min(*padfSrcMinX, dfSrcX);
2763 175554 : *padfSrcMinY = std::min(*padfSrcMinY, dfSrcY);
2764 175554 : *padfSrcMaxX = std::max(*padfSrcMaxX, dfSrcX);
2765 175554 : *padfSrcMaxY = std::max(*padfSrcMaxY, dfSrcY);
2766 : }
2767 666540 : iPoint++;
2768 : }
2769 : }
2770 : }
2771 :
2772 : /************************************************************************/
2773 : /* ComputeSourceWindowTransformPoints() */
2774 : /************************************************************************/
2775 :
2776 6703 : bool GDALWarpOperation::ComputeSourceWindowTransformPoints(
2777 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize, bool bUseGrid,
2778 : bool bAll, int nStepCount, bool bTryWithCheckWithInvertProj,
2779 : double &dfMinXOut, double &dfMinYOut, double &dfMaxXOut, double &dfMaxYOut,
2780 : int &nSamplePoints, int &nFailedCount)
2781 : {
2782 6703 : nSamplePoints = 0;
2783 6703 : nFailedCount = 0;
2784 :
2785 6703 : const double dfStepSize = bAll ? 0 : 1.0 / (nStepCount - 1);
2786 6703 : constexpr int knIntMax = std::numeric_limits<int>::max();
2787 6703 : int nSampleMax = 0;
2788 6703 : if (bUseGrid)
2789 : {
2790 1322 : if (bAll)
2791 : {
2792 0 : if (nDstXSize > knIntMax - 1 ||
2793 0 : nDstYSize > knIntMax / (nDstXSize + 1) - 1)
2794 : {
2795 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps");
2796 0 : return false;
2797 : }
2798 0 : nSampleMax = (nDstXSize + 1) * (nDstYSize + 1);
2799 : }
2800 : else
2801 : {
2802 1322 : if (nStepCount > knIntMax - 2 ||
2803 1322 : (nStepCount + 2) > knIntMax / (nStepCount + 2))
2804 : {
2805 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps : %d",
2806 : nStepCount);
2807 0 : return false;
2808 : }
2809 1322 : nSampleMax = (nStepCount + 2) * (nStepCount + 2);
2810 : }
2811 : }
2812 : else
2813 : {
2814 5381 : if (bAll)
2815 : {
2816 160 : if (nDstXSize > knIntMax / 2 - nDstYSize)
2817 : {
2818 : // Extremely unlikely !
2819 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps");
2820 0 : return false;
2821 : }
2822 160 : nSampleMax = 2 * (nDstXSize + nDstYSize);
2823 : }
2824 : else
2825 : {
2826 5221 : if (nStepCount > knIntMax / 4)
2827 : {
2828 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps : %d * 4",
2829 : nStepCount);
2830 0 : return false;
2831 : }
2832 5221 : nSampleMax = nStepCount * 4;
2833 : }
2834 : }
2835 :
2836 : int *pabSuccess =
2837 6703 : static_cast<int *>(VSI_MALLOC2_VERBOSE(sizeof(int), nSampleMax));
2838 : double *padfX = static_cast<double *>(
2839 6703 : VSI_MALLOC2_VERBOSE(sizeof(double) * 3, nSampleMax));
2840 6703 : if (pabSuccess == nullptr || padfX == nullptr)
2841 : {
2842 0 : CPLFree(padfX);
2843 0 : CPLFree(pabSuccess);
2844 0 : return false;
2845 : }
2846 6703 : double *padfY = padfX + nSampleMax;
2847 6703 : double *padfZ = padfX + static_cast<size_t>(nSampleMax) * 2;
2848 :
2849 : /* -------------------------------------------------------------------- */
2850 : /* Setup sample points on a grid pattern throughout the area. */
2851 : /* -------------------------------------------------------------------- */
2852 6703 : if (bUseGrid)
2853 : {
2854 1322 : if (bAll)
2855 : {
2856 0 : for (int iY = 0; iY <= nDstYSize; ++iY)
2857 : {
2858 0 : for (int iX = 0; iX <= nDstXSize; ++iX)
2859 : {
2860 0 : padfX[nSamplePoints] = nDstXOff + iX;
2861 0 : padfY[nSamplePoints] = nDstYOff + iY;
2862 0 : padfZ[nSamplePoints++] = 0.0;
2863 : }
2864 : }
2865 : }
2866 : else
2867 : {
2868 31728 : for (int iY = 0; iY < nStepCount + 2; iY++)
2869 : {
2870 59490 : const double dfRatioY = (iY == 0) ? 0.5 / nDstXSize
2871 : : (iY <= nStepCount)
2872 29084 : ? (iY - 1) * dfStepSize
2873 1322 : : 1 - 0.5 / nDstXSize;
2874 729744 : for (int iX = 0; iX < nStepCount + 2; iX++)
2875 : {
2876 1368270 : const double dfRatioX = (iX == 0) ? 0.5 / nDstXSize
2877 : : (iX <= nStepCount)
2878 668932 : ? (iX - 1) * dfStepSize
2879 30406 : : 1 - 0.5 / nDstXSize;
2880 699338 : padfX[nSamplePoints] = dfRatioX * nDstXSize + nDstXOff;
2881 699338 : padfY[nSamplePoints] = dfRatioY * nDstYSize + nDstYOff;
2882 699338 : padfZ[nSamplePoints++] = 0.0;
2883 : }
2884 : }
2885 : }
2886 : }
2887 : /* -------------------------------------------------------------------- */
2888 : /* Setup sample points all around the edge of the output raster. */
2889 : /* -------------------------------------------------------------------- */
2890 : else
2891 : {
2892 5381 : if (bAll)
2893 : {
2894 68927 : for (int iX = 0; iX <= nDstXSize; ++iX)
2895 : {
2896 : // Along top
2897 68767 : padfX[nSamplePoints] = nDstXOff + iX;
2898 68767 : padfY[nSamplePoints] = nDstYOff;
2899 68767 : padfZ[nSamplePoints++] = 0.0;
2900 :
2901 : // Along bottom
2902 68767 : padfX[nSamplePoints] = nDstXOff + iX;
2903 68767 : padfY[nSamplePoints] = nDstYOff + nDstYSize;
2904 68767 : padfZ[nSamplePoints++] = 0.0;
2905 : }
2906 :
2907 44154 : for (int iY = 1; iY < nDstYSize; ++iY)
2908 : {
2909 : // Along left
2910 43994 : padfX[nSamplePoints] = nDstXOff;
2911 43994 : padfY[nSamplePoints] = nDstYOff + iY;
2912 43994 : padfZ[nSamplePoints++] = 0.0;
2913 :
2914 : // Along right
2915 43994 : padfX[nSamplePoints] = nDstXOff + nDstXSize;
2916 43994 : padfY[nSamplePoints] = nDstYOff + iY;
2917 43994 : padfZ[nSamplePoints++] = 0.0;
2918 : }
2919 : }
2920 : else
2921 : {
2922 114862 : for (double dfRatio = 0.0; dfRatio <= 1.0 + dfStepSize * 0.5;
2923 109641 : dfRatio += dfStepSize)
2924 : {
2925 : // Along top
2926 109641 : padfX[nSamplePoints] = dfRatio * nDstXSize + nDstXOff;
2927 109641 : padfY[nSamplePoints] = nDstYOff;
2928 109641 : padfZ[nSamplePoints++] = 0.0;
2929 :
2930 : // Along bottom
2931 109641 : padfX[nSamplePoints] = dfRatio * nDstXSize + nDstXOff;
2932 109641 : padfY[nSamplePoints] = nDstYOff + nDstYSize;
2933 109641 : padfZ[nSamplePoints++] = 0.0;
2934 :
2935 : // Along left
2936 109641 : padfX[nSamplePoints] = nDstXOff;
2937 109641 : padfY[nSamplePoints] = dfRatio * nDstYSize + nDstYOff;
2938 109641 : padfZ[nSamplePoints++] = 0.0;
2939 :
2940 : // Along right
2941 109641 : padfX[nSamplePoints] = nDstXSize + nDstXOff;
2942 109641 : padfY[nSamplePoints] = dfRatio * nDstYSize + nDstYOff;
2943 109641 : padfZ[nSamplePoints++] = 0.0;
2944 : }
2945 : }
2946 : }
2947 :
2948 6703 : CPLAssert(nSamplePoints == nSampleMax);
2949 :
2950 : /* -------------------------------------------------------------------- */
2951 : /* Transform them to the input pixel coordinate space */
2952 : /* -------------------------------------------------------------------- */
2953 :
2954 138 : const auto RefreshTransformer = [this]()
2955 : {
2956 46 : if (GDALIsTransformer(psOptions->pTransformerArg,
2957 : GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
2958 : {
2959 0 : GDALRefreshGenImgProjTransformer(psOptions->pTransformerArg);
2960 : }
2961 46 : else if (GDALIsTransformer(psOptions->pTransformerArg,
2962 : GDAL_APPROX_TRANSFORMER_CLASS_NAME))
2963 : {
2964 46 : GDALRefreshApproxTransformer(psOptions->pTransformerArg);
2965 : }
2966 6749 : };
2967 :
2968 6703 : if (bTryWithCheckWithInvertProj)
2969 : {
2970 23 : CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ", "YES");
2971 23 : RefreshTransformer();
2972 : }
2973 6703 : psOptions->pfnTransformer(psOptions->pTransformerArg, TRUE, nSamplePoints,
2974 : padfX, padfY, padfZ, pabSuccess);
2975 6703 : if (bTryWithCheckWithInvertProj)
2976 : {
2977 23 : CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ", nullptr);
2978 23 : RefreshTransformer();
2979 : }
2980 :
2981 : /* -------------------------------------------------------------------- */
2982 : /* Collect the bounds, ignoring any failed points. */
2983 : /* -------------------------------------------------------------------- */
2984 1370130 : for (int i = 0; i < nSamplePoints; i++)
2985 : {
2986 1363420 : if (!pabSuccess[i])
2987 : {
2988 112179 : nFailedCount++;
2989 112179 : continue;
2990 : }
2991 :
2992 : // If this happens this is likely the symptom of a bug somewhere.
2993 1251240 : if (std::isnan(padfX[i]) || std::isnan(padfY[i]))
2994 : {
2995 : static bool bNanCoordFound = false;
2996 0 : if (!bNanCoordFound)
2997 : {
2998 0 : CPLDebug("WARP",
2999 : "ComputeSourceWindow(): "
3000 : "NaN coordinate found on point %d.",
3001 : i);
3002 0 : bNanCoordFound = true;
3003 : }
3004 0 : nFailedCount++;
3005 0 : continue;
3006 : }
3007 :
3008 1251240 : dfMinXOut = std::min(dfMinXOut, padfX[i]);
3009 1251240 : dfMinYOut = std::min(dfMinYOut, padfY[i]);
3010 1251240 : dfMaxXOut = std::max(dfMaxXOut, padfX[i]);
3011 1251240 : dfMaxYOut = std::max(dfMaxYOut, padfY[i]);
3012 : }
3013 :
3014 6703 : CPLFree(padfX);
3015 6703 : CPLFree(pabSuccess);
3016 6703 : return true;
3017 : }
3018 :
3019 : /************************************************************************/
3020 : /* ComputeSourceWindow() */
3021 : /************************************************************************/
3022 :
3023 : /** Given a target window starting at pixel (nDstOff, nDstYOff) and of
3024 : * dimension (nDstXSize, nDstYSize), compute the corresponding window in
3025 : * the source raster, and return the source position in (*pnSrcXOff, *pnSrcYOff),
3026 : * the source dimension in (*pnSrcXSize, *pnSrcYSize).
3027 : * If pdfSrcXExtraSize is not null, its pointed value will be filled with the
3028 : * number of extra source pixels in X dimension to acquire to take into account
3029 : * the size of the resampling kernel. Similarly for pdfSrcYExtraSize for the
3030 : * Y dimension.
3031 : * If pdfSrcFillRatio is not null, its pointed value will be filled with the
3032 : * the ratio of the clamped source raster window size over the unclamped source
3033 : * raster window size. When this ratio is too low, this might be an indication
3034 : * that it might be beneficial to split the target window to avoid requesting
3035 : * too many source pixels.
3036 : */
3037 6382 : CPLErr GDALWarpOperation::ComputeSourceWindow(
3038 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize, int *pnSrcXOff,
3039 : int *pnSrcYOff, int *pnSrcXSize, int *pnSrcYSize, double *pdfSrcXExtraSize,
3040 : double *pdfSrcYExtraSize, double *pdfSrcFillRatio)
3041 :
3042 : {
3043 : /* -------------------------------------------------------------------- */
3044 : /* Figure out whether we just want to do the usual "along the */
3045 : /* edge" sampling, or using a grid. The grid usage is */
3046 : /* important in some weird "inside out" cases like WGS84 to */
3047 : /* polar stereographic around the pole. Also figure out the */
3048 : /* sampling rate. */
3049 : /* -------------------------------------------------------------------- */
3050 6382 : int nStepCount = DEFAULT_STEP_COUNT;
3051 6382 : bool bAll = false;
3052 :
3053 : bool bUseGrid =
3054 6382 : CPLFetchBool(psOptions->papszWarpOptions, "SAMPLE_GRID", false);
3055 :
3056 : const char *pszSampleSteps =
3057 6382 : CSLFetchNameValue(psOptions->papszWarpOptions, "SAMPLE_STEPS");
3058 6382 : if (pszSampleSteps)
3059 : {
3060 94 : if (EQUAL(pszSampleSteps, "ALL"))
3061 : {
3062 94 : bAll = true;
3063 : }
3064 : else
3065 : {
3066 0 : nStepCount = atoi(pszSampleSteps);
3067 0 : nStepCount = std::max(2, nStepCount);
3068 : }
3069 : }
3070 6288 : else if (!bUseGrid)
3071 : {
3072 : // Detect if at least one of the 4 corner in destination raster fails
3073 : // to project back to source.
3074 : // Helps for long-lat to orthographic on areas that are partly in
3075 : // space / partly on Earth. Cf https://github.com/OSGeo/gdal/issues/9056
3076 : double adfCornerX[4];
3077 : double adfCornerY[4];
3078 5287 : double adfCornerZ[4] = {0, 0, 0, 0};
3079 5287 : int anCornerSuccess[4] = {FALSE, FALSE, FALSE, FALSE};
3080 5287 : adfCornerX[0] = nDstXOff;
3081 5287 : adfCornerY[0] = nDstYOff;
3082 5287 : adfCornerX[1] = nDstXOff + nDstXSize;
3083 5287 : adfCornerY[1] = nDstYOff;
3084 5287 : adfCornerX[2] = nDstXOff;
3085 5287 : adfCornerY[2] = nDstYOff + nDstYSize;
3086 5287 : adfCornerX[3] = nDstXOff + nDstXSize;
3087 5287 : adfCornerY[3] = nDstYOff + nDstYSize;
3088 5287 : if (!psOptions->pfnTransformer(psOptions->pTransformerArg, TRUE, 4,
3089 : adfCornerX, adfCornerY, adfCornerZ,
3090 5221 : anCornerSuccess) ||
3091 10508 : !anCornerSuccess[0] || !anCornerSuccess[1] || !anCornerSuccess[2] ||
3092 5221 : !anCornerSuccess[3])
3093 : {
3094 66 : bAll = true;
3095 : }
3096 : }
3097 :
3098 6382 : bool bTryWithCheckWithInvertProj = false;
3099 6382 : double dfMinXOut = std::numeric_limits<double>::infinity();
3100 6382 : double dfMinYOut = std::numeric_limits<double>::infinity();
3101 6382 : double dfMaxXOut = -std::numeric_limits<double>::infinity();
3102 6382 : double dfMaxYOut = -std::numeric_limits<double>::infinity();
3103 :
3104 6382 : int nSamplePoints = 0;
3105 6382 : int nFailedCount = 0;
3106 6382 : if (!ComputeSourceWindowTransformPoints(
3107 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, bUseGrid, bAll,
3108 : nStepCount, bTryWithCheckWithInvertProj, dfMinXOut, dfMinYOut,
3109 : dfMaxXOut, dfMaxYOut, nSamplePoints, nFailedCount))
3110 : {
3111 0 : return CE_Failure;
3112 : }
3113 :
3114 : // Use grid sampling as soon as a special point falls into the extent of
3115 : // the target raster.
3116 6382 : if (!bUseGrid && psOptions->hDstDS)
3117 : {
3118 13029 : for (const auto &xy : aDstXYSpecialPoints)
3119 : {
3120 17525 : if (0 <= xy.first &&
3121 1737 : GDALGetRasterXSize(psOptions->hDstDS) >= xy.first &&
3122 10427 : 0 <= xy.second &&
3123 796 : GDALGetRasterYSize(psOptions->hDstDS) >= xy.second)
3124 : {
3125 233 : bUseGrid = true;
3126 233 : bAll = false;
3127 233 : if (!ComputeSourceWindowTransformPoints(
3128 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, bUseGrid,
3129 : bAll, nStepCount, bTryWithCheckWithInvertProj,
3130 : dfMinXOut, dfMinYOut, dfMaxXOut, dfMaxYOut,
3131 : nSamplePoints, nFailedCount))
3132 : {
3133 0 : return CE_Failure;
3134 : }
3135 233 : break;
3136 : }
3137 : }
3138 : }
3139 :
3140 6382 : const int nRasterXSize = GDALGetRasterXSize(psOptions->hSrcDS);
3141 6382 : const int nRasterYSize = GDALGetRasterYSize(psOptions->hSrcDS);
3142 :
3143 : // Try to detect crazy values coming from reprojection that would not
3144 : // have resulted in a PROJ error. Could happen for example with PROJ
3145 : // <= 4.9.2 with inverse UTM/tmerc (Snyder approximation without sanity
3146 : // check) when being far away from the central meridian. But might be worth
3147 : // keeping that even for later versions in case some exotic projection isn't
3148 : // properly sanitized.
3149 6308 : if (nFailedCount == 0 && !bTryWithCheckWithInvertProj &&
3150 6308 : (dfMinXOut < -1e6 || dfMinYOut < -1e6 ||
3151 12694 : dfMaxXOut > nRasterXSize + 1e6 || dfMaxYOut > nRasterYSize + 1e6) &&
3152 23 : !CPLTestBool(CPLGetConfigOption("CHECK_WITH_INVERT_PROJ", "NO")))
3153 : {
3154 23 : CPLDebug("WARP",
3155 : "ComputeSourceWindow(): bogus source dataset window "
3156 : "returned. Trying again with CHECK_WITH_INVERT_PROJ=YES");
3157 23 : bTryWithCheckWithInvertProj = true;
3158 :
3159 : // We should probably perform the coordinate transformation in the
3160 : // warp kernel under CHECK_WITH_INVERT_PROJ too...
3161 23 : if (!ComputeSourceWindowTransformPoints(
3162 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, bUseGrid, bAll,
3163 : nStepCount, bTryWithCheckWithInvertProj, dfMinXOut, dfMinYOut,
3164 : dfMaxXOut, dfMaxYOut, nSamplePoints, nFailedCount))
3165 : {
3166 0 : return CE_Failure;
3167 : }
3168 : }
3169 :
3170 : /* -------------------------------------------------------------------- */
3171 : /* If we got any failures when not using a grid, we should */
3172 : /* really go back and try again with the grid. Sorry for the */
3173 : /* goto. */
3174 : /* -------------------------------------------------------------------- */
3175 6382 : if (!bUseGrid && nFailedCount > 0)
3176 : {
3177 65 : bUseGrid = true;
3178 65 : bAll = false;
3179 65 : if (!ComputeSourceWindowTransformPoints(
3180 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, bUseGrid, bAll,
3181 : nStepCount, bTryWithCheckWithInvertProj, dfMinXOut, dfMinYOut,
3182 : dfMaxXOut, dfMaxYOut, nSamplePoints, nFailedCount))
3183 : {
3184 0 : return CE_Failure;
3185 : }
3186 : }
3187 :
3188 : /* -------------------------------------------------------------------- */
3189 : /* If we get hardly any points (or none) transforming, we give */
3190 : /* up. */
3191 : /* -------------------------------------------------------------------- */
3192 6382 : if (nFailedCount > nSamplePoints - 5)
3193 : {
3194 : const bool bErrorOutIfEmptySourceWindow =
3195 39 : CPLFetchBool(psOptions->papszWarpOptions,
3196 : "ERROR_OUT_IF_EMPTY_SOURCE_WINDOW", true);
3197 39 : if (bErrorOutIfEmptySourceWindow)
3198 : {
3199 3 : CPLError(CE_Failure, CPLE_AppDefined,
3200 : "Too many points (%d out of %d) failed to transform, "
3201 : "unable to compute output bounds.",
3202 : nFailedCount, nSamplePoints);
3203 : }
3204 : else
3205 : {
3206 36 : CPLDebug("WARP", "Cannot determine source window for %d,%d,%d,%d",
3207 : nDstXOff, nDstYOff, nDstXSize, nDstYSize);
3208 : }
3209 39 : return CE_Failure;
3210 : }
3211 :
3212 6343 : if (nFailedCount > 0)
3213 39 : CPLDebug("GDAL",
3214 : "GDALWarpOperation::ComputeSourceWindow() %d out of %d "
3215 : "points failed to transform.",
3216 : nFailedCount, nSamplePoints);
3217 :
3218 : /* -------------------------------------------------------------------- */
3219 : /* In some cases (see https://github.com/OSGeo/gdal/issues/862) */
3220 : /* the reverse transform does not work at some points, so try by */
3221 : /* transforming from source raster space to target raster space and */
3222 : /* see which source coordinates end up being in the AOI in the target */
3223 : /* raster space. */
3224 : /* -------------------------------------------------------------------- */
3225 6343 : if (bUseGrid)
3226 : {
3227 1260 : ComputeSourceWindowStartingFromSource(nDstXOff, nDstYOff, nDstXSize,
3228 : nDstYSize, &dfMinXOut, &dfMinYOut,
3229 : &dfMaxXOut, &dfMaxYOut);
3230 : }
3231 :
3232 : /* -------------------------------------------------------------------- */
3233 : /* Early exit to avoid crazy values to cause a huge nResWinSize that */
3234 : /* would result in a result window wrongly covering the whole raster. */
3235 : /* -------------------------------------------------------------------- */
3236 6343 : if (dfMinXOut > nRasterXSize || dfMaxXOut < 0 || dfMinYOut > nRasterYSize ||
3237 5273 : dfMaxYOut < 0)
3238 : {
3239 1531 : *pnSrcXOff = 0;
3240 1531 : *pnSrcYOff = 0;
3241 1531 : *pnSrcXSize = 0;
3242 1531 : *pnSrcYSize = 0;
3243 1531 : if (pdfSrcXExtraSize)
3244 1531 : *pdfSrcXExtraSize = 0.0;
3245 1531 : if (pdfSrcYExtraSize)
3246 1531 : *pdfSrcYExtraSize = 0.0;
3247 1531 : if (pdfSrcFillRatio)
3248 949 : *pdfSrcFillRatio = 0.0;
3249 1531 : return CE_None;
3250 : }
3251 :
3252 : // For scenarios where warping is used as a "decoration", try to clamp
3253 : // source pixel coordinates to integer when very close.
3254 19248 : const auto roundIfCloseEnough = [](double dfVal)
3255 : {
3256 19248 : const double dfRounded = std::round(dfVal);
3257 19248 : if (std::fabs(dfRounded - dfVal) < 1e-6)
3258 15213 : return dfRounded;
3259 4035 : return dfVal;
3260 : };
3261 :
3262 4812 : dfMinXOut = roundIfCloseEnough(dfMinXOut);
3263 4812 : dfMinYOut = roundIfCloseEnough(dfMinYOut);
3264 4812 : dfMaxXOut = roundIfCloseEnough(dfMaxXOut);
3265 4812 : dfMaxYOut = roundIfCloseEnough(dfMaxYOut);
3266 :
3267 4812 : if (m_bIsTranslationOnPixelBoundaries)
3268 : {
3269 376 : CPLAssert(dfMinXOut == std::round(dfMinXOut));
3270 376 : CPLAssert(dfMinYOut == std::round(dfMinYOut));
3271 376 : CPLAssert(dfMaxXOut == std::round(dfMaxXOut));
3272 376 : CPLAssert(dfMaxYOut == std::round(dfMaxYOut));
3273 376 : CPLAssert(std::round(dfMaxXOut - dfMinXOut) == nDstXSize);
3274 376 : CPLAssert(std::round(dfMaxYOut - dfMinYOut) == nDstYSize);
3275 : }
3276 :
3277 : /* -------------------------------------------------------------------- */
3278 : /* How much of a window around our source pixel might we need */
3279 : /* to collect data from based on the resampling kernel? Even */
3280 : /* if the requested central pixel falls off the source image, */
3281 : /* we may need to collect data if some portion of the */
3282 : /* resampling kernel could be on-image. */
3283 : /* -------------------------------------------------------------------- */
3284 4812 : const int nResWinSize = m_bIsTranslationOnPixelBoundaries
3285 4812 : ? 0
3286 4436 : : GWKGetFilterRadius(psOptions->eResampleAlg);
3287 :
3288 : // Take scaling into account.
3289 : // Avoid ridiculous small scaling factors to avoid potential further integer
3290 : // overflows
3291 9624 : const double dfXScale = std::max(1e-3, static_cast<double>(nDstXSize) /
3292 4812 : (dfMaxXOut - dfMinXOut));
3293 9624 : const double dfYScale = std::max(1e-3, static_cast<double>(nDstYSize) /
3294 4812 : (dfMaxYOut - dfMinYOut));
3295 4812 : int nXRadius = dfXScale < 0.95
3296 4812 : ? static_cast<int>(ceil(nResWinSize / dfXScale))
3297 : : nResWinSize;
3298 4812 : int nYRadius = dfYScale < 0.95
3299 4812 : ? static_cast<int>(ceil(nResWinSize / dfYScale))
3300 : : nResWinSize;
3301 :
3302 : /* -------------------------------------------------------------------- */
3303 : /* Allow addition of extra sample pixels to source window to */
3304 : /* avoid missing pixels due to sampling error. In fact, */
3305 : /* fallback to adding a bit to the window if any points failed */
3306 : /* to transform. */
3307 : /* -------------------------------------------------------------------- */
3308 4812 : if (const char *pszSourceExtra =
3309 4812 : CSLFetchNameValue(psOptions->papszWarpOptions, "SOURCE_EXTRA"))
3310 : {
3311 81 : int nSrcExtra = cpl::strict_parse<int>(pszSourceExtra).value_or(-1);
3312 :
3313 81 : if (nSrcExtra < 0)
3314 : {
3315 : // no point raising CE_Failure because it will get converted into
3316 : // a warning at an outer scope
3317 1 : CPLError(CE_Warning, CPLE_IllegalArg,
3318 : "SOURCE_EXTRA must be a positive integer or zero.");
3319 1 : nSrcExtra = 0;
3320 : }
3321 :
3322 81 : nXRadius += nSrcExtra;
3323 81 : nYRadius += nSrcExtra;
3324 : }
3325 4731 : else if (nFailedCount > 0)
3326 : {
3327 33 : nXRadius += 10;
3328 33 : nYRadius += 10;
3329 : }
3330 :
3331 : /* -------------------------------------------------------------------- */
3332 : /* return bounds. */
3333 : /* -------------------------------------------------------------------- */
3334 : #if DEBUG_VERBOSE
3335 : CPLDebug("WARP",
3336 : "dst=(%d,%d,%d,%d) raw "
3337 : "src=(minx=%.17g,miny=%.17g,maxx=%.17g,maxy=%.17g)",
3338 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, dfMinXOut, dfMinYOut,
3339 : dfMaxXOut, dfMaxYOut);
3340 : #endif
3341 4812 : const int nMinXOutClamped = static_cast<int>(std::max(0.0, dfMinXOut));
3342 4812 : const int nMinYOutClamped = static_cast<int>(std::max(0.0, dfMinYOut));
3343 : const int nMaxXOutClamped = static_cast<int>(
3344 4812 : std::min(ceil(dfMaxXOut), static_cast<double>(nRasterXSize)));
3345 : const int nMaxYOutClamped = static_cast<int>(
3346 4812 : std::min(ceil(dfMaxYOut), static_cast<double>(nRasterYSize)));
3347 :
3348 : const double dfSrcXSizeRaw = std::max(
3349 14436 : 0.0, std::min(static_cast<double>(nRasterXSize - nMinXOutClamped),
3350 4812 : dfMaxXOut - dfMinXOut));
3351 : const double dfSrcYSizeRaw = std::max(
3352 14436 : 0.0, std::min(static_cast<double>(nRasterYSize - nMinYOutClamped),
3353 4812 : dfMaxYOut - dfMinYOut));
3354 :
3355 : // If we cover more than 90% of the width, then use it fully (helps for
3356 : // anti-meridian discontinuities)
3357 4812 : if (nMaxXOutClamped - nMinXOutClamped > 0.9 * nRasterXSize)
3358 : {
3359 1595 : *pnSrcXOff = 0;
3360 1595 : *pnSrcXSize = nRasterXSize;
3361 : }
3362 : else
3363 : {
3364 3217 : *pnSrcXOff =
3365 3217 : std::max(0, std::min(nMinXOutClamped - nXRadius, nRasterXSize));
3366 3217 : *pnSrcXSize =
3367 9651 : std::max(0, std::min(nRasterXSize - *pnSrcXOff,
3368 3217 : nMaxXOutClamped - *pnSrcXOff + nXRadius));
3369 : }
3370 :
3371 4812 : if (nMaxYOutClamped - nMinYOutClamped > 0.9 * nRasterYSize)
3372 : {
3373 1493 : *pnSrcYOff = 0;
3374 1493 : *pnSrcYSize = nRasterYSize;
3375 : }
3376 : else
3377 : {
3378 3319 : *pnSrcYOff =
3379 3319 : std::max(0, std::min(nMinYOutClamped - nYRadius, nRasterYSize));
3380 3319 : *pnSrcYSize =
3381 9957 : std::max(0, std::min(nRasterYSize - *pnSrcYOff,
3382 3319 : nMaxYOutClamped - *pnSrcYOff + nYRadius));
3383 : }
3384 :
3385 4812 : if (pdfSrcXExtraSize)
3386 4812 : *pdfSrcXExtraSize = *pnSrcXSize - dfSrcXSizeRaw;
3387 4812 : if (pdfSrcYExtraSize)
3388 4812 : *pdfSrcYExtraSize = *pnSrcYSize - dfSrcYSizeRaw;
3389 :
3390 : // Computed the ratio of the clamped source raster window size over
3391 : // the unclamped source raster window size.
3392 4812 : if (pdfSrcFillRatio)
3393 3653 : *pdfSrcFillRatio =
3394 7306 : static_cast<double>(*pnSrcXSize) * (*pnSrcYSize) /
3395 3653 : std::max(1.0, (dfMaxXOut - dfMinXOut + 2 * nXRadius) *
3396 3653 : (dfMaxYOut - dfMinYOut + 2 * nYRadius));
3397 :
3398 4812 : return CE_None;
3399 : }
3400 :
3401 : /************************************************************************/
3402 : /* ReportTiming() */
3403 : /************************************************************************/
3404 :
3405 12735 : void GDALWarpOperation::ReportTiming(const char *pszMessage)
3406 :
3407 : {
3408 12735 : if (!bReportTimings)
3409 12735 : return;
3410 :
3411 0 : const unsigned long nNewTime = VSITime(nullptr);
3412 :
3413 0 : if (pszMessage != nullptr)
3414 : {
3415 0 : CPLDebug("WARP_TIMING", "%s: %lds", pszMessage,
3416 0 : static_cast<long>(nNewTime - nLastTimeReported));
3417 : }
3418 :
3419 0 : nLastTimeReported = nNewTime;
3420 : }
|
