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 1261 : GetWarpPrivateData(GDALWarpOperation *poWarpOperation)
65 : {
66 2522 : std::lock_guard<std::mutex> oLock(gMutex);
67 1261 : auto oItem = gMapPrivate.find(poWarpOperation);
68 1261 : if (oItem != gMapPrivate.end())
69 : {
70 949 : 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 : if (oResetDestPixels.value())
740 : {
741 4 : for (int i = 0; eErr == CE_None && i < psOptions->nBandCount; ++i)
742 : {
743 2 : eErr =
744 5 : GDALFillRaster(GDALGetRasterBand(psOptions->hDstDS,
745 2 : psOptions->panDstBands[i]),
746 2 : psOptions->padfDstNoDataReal
747 1 : ? psOptions->padfDstNoDataReal[i]
748 : : 0.0,
749 2 : psOptions->padfDstNoDataImag
750 0 : ? psOptions->padfDstNoDataImag[i]
751 : : 0.0);
752 : }
753 : }
754 : }
755 :
756 1830 : return eErr;
757 : }
758 :
759 : /**
760 : * \fn void* GDALWarpOperation::CreateDestinationBuffer(
761 : int nDstXSize, int nDstYSize, int *pbInitialized);
762 : *
763 : * This method creates a destination buffer for use with WarpRegionToBuffer.
764 : * The output is initialized based on the INIT_DEST settings.
765 : *
766 : * @param nDstXSize Width of output window on destination buffer to be produced.
767 : * @param nDstYSize Height of output window on destination buffer to be
768 : produced.
769 : * @param pbInitialized Filled with boolean indicating if the buffer was
770 : initialized.
771 : *
772 : * @return Buffer capable for use as a warp operation output destination
773 : */
774 3386 : void *GDALWarpOperation::CreateDestinationBuffer(int nDstXSize, int nDstYSize,
775 : int *pbInitialized)
776 : {
777 :
778 : /* -------------------------------------------------------------------- */
779 : /* Allocate block of memory large enough to hold all the bands */
780 : /* for this block. */
781 : /* -------------------------------------------------------------------- */
782 3386 : const int nWordSize = GDALGetDataTypeSizeBytes(psOptions->eWorkingDataType);
783 :
784 3386 : void *pDstBuffer = VSI_MALLOC3_VERBOSE(
785 : cpl::fits_on<int>(nWordSize * psOptions->nBandCount), nDstXSize,
786 : nDstYSize);
787 3386 : if (pDstBuffer)
788 : {
789 3386 : auto eErr = InitializeDestinationBuffer(pDstBuffer, nDstXSize,
790 : nDstYSize, pbInitialized);
791 3386 : if (eErr != CE_None)
792 : {
793 2 : CPLFree(pDstBuffer);
794 2 : return nullptr;
795 : }
796 : }
797 3384 : return pDstBuffer;
798 : }
799 :
800 : /**
801 : * This method initializes a destination buffer for use with WarpRegionToBuffer.
802 : *
803 : * It is initialized based on the INIT_DEST settings.
804 : *
805 : * This method is called by CreateDestinationBuffer().
806 : * It is meant at being used by callers that have already allocated the
807 : * destination buffer without using CreateDestinationBuffer().
808 : *
809 : * @param pDstBuffer Buffer of size
810 : * GDALGetDataTypeSizeBytes(psOptions->eWorkingDataType) *
811 : * nDstXSize * nDstYSize * psOptions->nBandCount bytes.
812 : * @param nDstXSize Width of output window on destination buffer to be produced.
813 : * @param nDstYSize Height of output window on destination buffer to be
814 : * produced.
815 : * @param pbInitialized Filled with boolean indicating if the buffer was
816 : * initialized.
817 : * @since 3.10
818 : */
819 3485 : CPLErr GDALWarpOperation::InitializeDestinationBuffer(void *pDstBuffer,
820 : int nDstXSize,
821 : int nDstYSize,
822 : int *pbInitialized) const
823 : {
824 3485 : const int nWordSize = GDALGetDataTypeSizeBytes(psOptions->eWorkingDataType);
825 :
826 3485 : const GPtrDiff_t nBandSize =
827 3485 : static_cast<GPtrDiff_t>(nWordSize) * nDstXSize * nDstYSize;
828 :
829 : /* -------------------------------------------------------------------- */
830 : /* Initialize if requested in the options */
831 : /* -------------------------------------------------------------------- */
832 : const char *pszInitDest =
833 3485 : CSLFetchNameValue(psOptions->papszWarpOptions, "INIT_DEST");
834 :
835 3485 : if (pszInitDest == nullptr || EQUAL(pszInitDest, ""))
836 : {
837 384 : if (pbInitialized != nullptr)
838 : {
839 384 : *pbInitialized = FALSE;
840 : }
841 384 : return CE_None;
842 : }
843 :
844 3101 : if (pbInitialized != nullptr)
845 : {
846 1938 : *pbInitialized = TRUE;
847 : }
848 :
849 : CPLStringList aosInitValues(
850 6202 : CSLTokenizeStringComplex(pszInitDest, ",", FALSE, FALSE));
851 3101 : const int nInitCount = aosInitValues.Count();
852 :
853 8306 : for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
854 : {
855 5207 : double adfInitRealImag[2] = {0.0, 0.0};
856 : const char *pszBandInit =
857 5207 : aosInitValues[std::min(iBand, nInitCount - 1)];
858 :
859 5207 : if (EQUAL(pszBandInit, "NO_DATA"))
860 : {
861 709 : if (psOptions->padfDstNoDataReal == nullptr)
862 : {
863 1 : CPLError(CE_Failure, CPLE_AppDefined,
864 : "INIT_DEST was set to NO_DATA, but a NoData value was "
865 : "not defined.");
866 : }
867 : else
868 : {
869 708 : adfInitRealImag[0] = psOptions->padfDstNoDataReal[iBand];
870 708 : if (psOptions->padfDstNoDataImag != nullptr)
871 : {
872 601 : adfInitRealImag[1] = psOptions->padfDstNoDataImag[iBand];
873 : }
874 : }
875 : }
876 : else
877 : {
878 4498 : if (CPLStringToComplex(pszBandInit, &adfInitRealImag[0],
879 4498 : &adfInitRealImag[1]) != CE_None)
880 : {
881 2 : CPLError(CE_Failure, CPLE_AppDefined,
882 : "Error parsing INIT_DEST");
883 2 : return CE_Failure;
884 : }
885 : }
886 :
887 5205 : GByte *pBandData = static_cast<GByte *>(pDstBuffer) + iBand * nBandSize;
888 :
889 5205 : if (psOptions->eWorkingDataType == GDT_UInt8)
890 : {
891 9112 : memset(pBandData,
892 : std::max(
893 4556 : 0, std::min(255, static_cast<int>(adfInitRealImag[0]))),
894 : nBandSize);
895 : }
896 1293 : else if (!std::isnan(adfInitRealImag[0]) && adfInitRealImag[0] == 0.0 &&
897 1293 : !std::isnan(adfInitRealImag[1]) && adfInitRealImag[1] == 0.0)
898 : {
899 559 : memset(pBandData, 0, nBandSize);
900 : }
901 90 : else if (!std::isnan(adfInitRealImag[1]) && adfInitRealImag[1] == 0.0)
902 : {
903 90 : GDALCopyWords64(&adfInitRealImag, GDT_Float64, 0, pBandData,
904 90 : psOptions->eWorkingDataType, nWordSize,
905 90 : static_cast<GPtrDiff_t>(nDstXSize) * nDstYSize);
906 : }
907 : else
908 : {
909 0 : GDALCopyWords64(&adfInitRealImag, GDT_CFloat64, 0, pBandData,
910 0 : psOptions->eWorkingDataType, nWordSize,
911 0 : static_cast<GPtrDiff_t>(nDstXSize) * nDstYSize);
912 : }
913 : }
914 :
915 3099 : return CE_None;
916 : }
917 :
918 : /**
919 : * \fn void GDALWarpOperation::DestroyDestinationBuffer( void *pDstBuffer )
920 : *
921 : * This method destroys a buffer previously retrieved from
922 : * CreateDestinationBuffer
923 : *
924 : * @param pDstBuffer destination buffer to be destroyed
925 : *
926 : */
927 3384 : void GDALWarpOperation::DestroyDestinationBuffer(void *pDstBuffer)
928 : {
929 3384 : VSIFree(pDstBuffer);
930 3384 : }
931 :
932 : /************************************************************************/
933 : /* GDALCreateWarpOperation() */
934 : /************************************************************************/
935 :
936 : /**
937 : * @see GDALWarpOperation::Initialize()
938 : */
939 :
940 149 : GDALWarpOperationH GDALCreateWarpOperation(const GDALWarpOptions *psNewOptions)
941 : {
942 149 : GDALWarpOperation *poOperation = new GDALWarpOperation;
943 149 : if (poOperation->Initialize(psNewOptions) != CE_None)
944 : {
945 0 : delete poOperation;
946 0 : return nullptr;
947 : }
948 :
949 149 : return reinterpret_cast<GDALWarpOperationH>(poOperation);
950 : }
951 :
952 : /************************************************************************/
953 : /* GDALDestroyWarpOperation() */
954 : /************************************************************************/
955 :
956 : /**
957 : * @see GDALWarpOperation::~GDALWarpOperation()
958 : */
959 :
960 149 : void GDALDestroyWarpOperation(GDALWarpOperationH hOperation)
961 : {
962 149 : if (hOperation)
963 149 : delete static_cast<GDALWarpOperation *>(hOperation);
964 149 : }
965 :
966 : /************************************************************************/
967 : /* CollectChunkList() */
968 : /************************************************************************/
969 :
970 1226 : void GDALWarpOperation::CollectChunkList(int nDstXOff, int nDstYOff,
971 : int nDstXSize, int nDstYSize)
972 :
973 : {
974 : /* -------------------------------------------------------------------- */
975 : /* Collect the list of chunks to operate on. */
976 : /* -------------------------------------------------------------------- */
977 1226 : WipeChunkList();
978 1226 : CollectChunkListInternal(nDstXOff, nDstYOff, nDstXSize, nDstYSize);
979 :
980 : // Sort chunks from top to bottom, and for equal y, from left to right.
981 1226 : if (nChunkListCount > 1)
982 : {
983 56 : std::sort(pasChunkList, pasChunkList + nChunkListCount,
984 7201 : [](const GDALWarpChunk &a, const GDALWarpChunk &b)
985 : {
986 7201 : if (a.dy < b.dy)
987 3220 : return true;
988 3981 : if (a.dy > b.dy)
989 1415 : return false;
990 2566 : return a.dx < b.dx;
991 : });
992 : }
993 :
994 : /* -------------------------------------------------------------------- */
995 : /* Find the global source window. */
996 : /* -------------------------------------------------------------------- */
997 :
998 1226 : const int knIntMax = std::numeric_limits<int>::max();
999 1226 : const int knIntMin = std::numeric_limits<int>::min();
1000 1226 : int nSrcXOff = knIntMax;
1001 1226 : int nSrcYOff = knIntMax;
1002 1226 : int nSrcX2Off = knIntMin;
1003 1226 : int nSrcY2Off = knIntMin;
1004 1226 : double dfApproxAccArea = 0;
1005 3548 : for (int iChunk = 0; pasChunkList != nullptr && iChunk < nChunkListCount;
1006 : iChunk++)
1007 : {
1008 2322 : GDALWarpChunk *pasThisChunk = pasChunkList + iChunk;
1009 2322 : nSrcXOff = std::min(nSrcXOff, pasThisChunk->sx);
1010 2322 : nSrcYOff = std::min(nSrcYOff, pasThisChunk->sy);
1011 2322 : nSrcX2Off = std::max(nSrcX2Off, pasThisChunk->sx + pasThisChunk->ssx);
1012 2322 : nSrcY2Off = std::max(nSrcY2Off, pasThisChunk->sy + pasThisChunk->ssy);
1013 2322 : dfApproxAccArea +=
1014 2322 : static_cast<double>(pasThisChunk->ssx) * pasThisChunk->ssy;
1015 : }
1016 1226 : if (nSrcXOff < nSrcX2Off)
1017 : {
1018 1218 : const double dfTotalArea =
1019 1218 : static_cast<double>(nSrcX2Off - nSrcXOff) * (nSrcY2Off - nSrcYOff);
1020 : // This is really a gross heuristics, but should work in most cases
1021 1218 : if (dfApproxAccArea >= dfTotalArea * 0.80)
1022 : {
1023 1218 : GDALDataset::FromHandle(psOptions->hSrcDS)
1024 1218 : ->AdviseRead(nSrcXOff, nSrcYOff, nSrcX2Off - nSrcXOff,
1025 : nSrcY2Off - nSrcYOff, nDstXSize, nDstYSize,
1026 1218 : psOptions->eWorkingDataType, psOptions->nBandCount,
1027 1218 : psOptions->panSrcBands, nullptr);
1028 : }
1029 : }
1030 1226 : }
1031 :
1032 : /************************************************************************/
1033 : /* ChunkAndWarpImage() */
1034 : /************************************************************************/
1035 :
1036 : /**
1037 : * \fn CPLErr GDALWarpOperation::ChunkAndWarpImage(
1038 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize );
1039 : *
1040 : * This method does a complete warp of the source image to the destination
1041 : * image for the indicated region with the current warp options in effect.
1042 : * Progress is reported to the installed progress monitor, if any.
1043 : *
1044 : * This function will subdivide the region and recursively call itself
1045 : * until the total memory required to process a region chunk will all fit
1046 : * in the memory pool defined by GDALWarpOptions::dfWarpMemoryLimit.
1047 : *
1048 : * Once an appropriate region is selected GDALWarpOperation::WarpRegion()
1049 : * is invoked to do the actual work.
1050 : *
1051 : * @param nDstXOff X offset to window of destination data to be produced.
1052 : * @param nDstYOff Y offset to window of destination data to be produced.
1053 : * @param nDstXSize Width of output window on destination file to be produced.
1054 : * @param nDstYSize Height of output window on destination file to be produced.
1055 : *
1056 : * @return CE_None on success or CE_Failure if an error occurs.
1057 : */
1058 :
1059 1220 : CPLErr GDALWarpOperation::ChunkAndWarpImage(int nDstXOff, int nDstYOff,
1060 : int nDstXSize, int nDstYSize)
1061 :
1062 : {
1063 : /* -------------------------------------------------------------------- */
1064 : /* Collect the list of chunks to operate on. */
1065 : /* -------------------------------------------------------------------- */
1066 1220 : CollectChunkList(nDstXOff, nDstYOff, nDstXSize, nDstYSize);
1067 :
1068 : /* -------------------------------------------------------------------- */
1069 : /* Total up output pixels to process. */
1070 : /* -------------------------------------------------------------------- */
1071 1220 : double dfTotalPixels = 0.0;
1072 :
1073 3531 : for (int iChunk = 0; pasChunkList != nullptr && iChunk < nChunkListCount;
1074 : iChunk++)
1075 : {
1076 2311 : GDALWarpChunk *pasThisChunk = pasChunkList + iChunk;
1077 2311 : const double dfChunkPixels =
1078 2311 : pasThisChunk->dsx * static_cast<double>(pasThisChunk->dsy);
1079 :
1080 2311 : dfTotalPixels += dfChunkPixels;
1081 : }
1082 :
1083 : /* -------------------------------------------------------------------- */
1084 : /* Process them one at a time, updating the progress */
1085 : /* information for each region. */
1086 : /* -------------------------------------------------------------------- */
1087 1220 : double dfPixelsProcessed = 0.0;
1088 :
1089 3525 : for (int iChunk = 0; pasChunkList != nullptr && iChunk < nChunkListCount;
1090 : iChunk++)
1091 : {
1092 2311 : GDALWarpChunk *pasThisChunk = pasChunkList + iChunk;
1093 2311 : const double dfChunkPixels =
1094 2311 : pasThisChunk->dsx * static_cast<double>(pasThisChunk->dsy);
1095 :
1096 2311 : const double dfProgressBase = dfPixelsProcessed / dfTotalPixels;
1097 2311 : const double dfProgressScale = dfChunkPixels / dfTotalPixels;
1098 :
1099 2311 : CPLErr eErr = WarpRegion(
1100 : pasThisChunk->dx, pasThisChunk->dy, pasThisChunk->dsx,
1101 : pasThisChunk->dsy, pasThisChunk->sx, pasThisChunk->sy,
1102 : pasThisChunk->ssx, pasThisChunk->ssy, pasThisChunk->sExtraSx,
1103 : pasThisChunk->sExtraSy, dfProgressBase, dfProgressScale);
1104 :
1105 2311 : if (eErr != CE_None)
1106 6 : return eErr;
1107 :
1108 2305 : dfPixelsProcessed += dfChunkPixels;
1109 : }
1110 :
1111 1214 : WipeChunkList();
1112 :
1113 1214 : psOptions->pfnProgress(1.0, "", psOptions->pProgressArg);
1114 :
1115 1214 : return CE_None;
1116 : }
1117 :
1118 : /************************************************************************/
1119 : /* GDALChunkAndWarpImage() */
1120 : /************************************************************************/
1121 :
1122 : /**
1123 : * @see GDALWarpOperation::ChunkAndWarpImage()
1124 : */
1125 :
1126 149 : CPLErr GDALChunkAndWarpImage(GDALWarpOperationH hOperation, int nDstXOff,
1127 : int nDstYOff, int nDstXSize, int nDstYSize)
1128 : {
1129 149 : VALIDATE_POINTER1(hOperation, "GDALChunkAndWarpImage", CE_Failure);
1130 :
1131 : return reinterpret_cast<GDALWarpOperation *>(hOperation)
1132 149 : ->ChunkAndWarpImage(nDstXOff, nDstYOff, nDstXSize, nDstYSize);
1133 : }
1134 :
1135 : /************************************************************************/
1136 : /* ChunkThreadMain() */
1137 : /************************************************************************/
1138 :
1139 : struct ChunkThreadData
1140 : {
1141 : GDALWarpOperation *poOperation = nullptr;
1142 : GDALWarpChunk *pasChunkInfo = nullptr;
1143 : CPLJoinableThread *hThreadHandle = nullptr;
1144 : CPLErr eErr = CE_None;
1145 : double dfProgressBase = 0;
1146 : double dfProgressScale = 0;
1147 : CPLMutex *hIOMutex = nullptr;
1148 :
1149 : CPLMutex *hCondMutex = nullptr;
1150 : volatile int bIOMutexTaken = 0;
1151 : CPLCond *hCond = nullptr;
1152 :
1153 : CPLErrorAccumulator *poErrorAccumulator = nullptr;
1154 : };
1155 :
1156 11 : static void ChunkThreadMain(void *pThreadData)
1157 :
1158 : {
1159 11 : volatile ChunkThreadData *psData =
1160 : static_cast<volatile ChunkThreadData *>(pThreadData);
1161 :
1162 11 : GDALWarpChunk *pasChunkInfo = psData->pasChunkInfo;
1163 :
1164 : /* -------------------------------------------------------------------- */
1165 : /* Acquire IO mutex. */
1166 : /* -------------------------------------------------------------------- */
1167 11 : if (!CPLAcquireMutex(psData->hIOMutex, 600.0))
1168 : {
1169 0 : CPLError(CE_Failure, CPLE_AppDefined,
1170 : "Failed to acquire IOMutex in WarpRegion().");
1171 0 : psData->eErr = CE_Failure;
1172 : }
1173 : else
1174 : {
1175 11 : if (psData->hCond != nullptr)
1176 : {
1177 6 : CPLAcquireMutex(psData->hCondMutex, 1.0);
1178 6 : psData->bIOMutexTaken = TRUE;
1179 6 : CPLCondSignal(psData->hCond);
1180 6 : CPLReleaseMutex(psData->hCondMutex);
1181 : }
1182 :
1183 : auto oAccumulator =
1184 22 : psData->poErrorAccumulator->InstallForCurrentScope();
1185 11 : CPL_IGNORE_RET_VAL(oAccumulator);
1186 :
1187 22 : psData->eErr = psData->poOperation->WarpRegion(
1188 : pasChunkInfo->dx, pasChunkInfo->dy, pasChunkInfo->dsx,
1189 : pasChunkInfo->dsy, pasChunkInfo->sx, pasChunkInfo->sy,
1190 : pasChunkInfo->ssx, pasChunkInfo->ssy, pasChunkInfo->sExtraSx,
1191 11 : pasChunkInfo->sExtraSy, psData->dfProgressBase,
1192 11 : psData->dfProgressScale);
1193 :
1194 : /* --------------------------------------------------------------------
1195 : */
1196 : /* Release the IO mutex. */
1197 : /* --------------------------------------------------------------------
1198 : */
1199 11 : CPLReleaseMutex(psData->hIOMutex);
1200 : }
1201 11 : }
1202 :
1203 : /************************************************************************/
1204 : /* ChunkAndWarpMulti() */
1205 : /************************************************************************/
1206 :
1207 : /**
1208 : * \fn CPLErr GDALWarpOperation::ChunkAndWarpMulti(
1209 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize );
1210 : *
1211 : * This method does a complete warp of the source image to the destination
1212 : * image for the indicated region with the current warp options in effect.
1213 : * Progress is reported to the installed progress monitor, if any.
1214 : *
1215 : * Externally this method operates the same as ChunkAndWarpImage(), but
1216 : * internally this method uses multiple threads to interleave input/output
1217 : * for one region while the processing is being done for another.
1218 : *
1219 : * @param nDstXOff X offset to window of destination data to be produced.
1220 : * @param nDstYOff Y offset to window of destination data to be produced.
1221 : * @param nDstXSize Width of output window on destination file to be produced.
1222 : * @param nDstYSize Height of output window on destination file to be produced.
1223 : *
1224 : * @return CE_None on success or CE_Failure if an error occurs.
1225 : */
1226 :
1227 6 : CPLErr GDALWarpOperation::ChunkAndWarpMulti(int nDstXOff, int nDstYOff,
1228 : int nDstXSize, int nDstYSize)
1229 :
1230 : {
1231 6 : hIOMutex = CPLCreateMutex();
1232 6 : hWarpMutex = CPLCreateMutex();
1233 :
1234 6 : CPLReleaseMutex(hIOMutex);
1235 6 : CPLReleaseMutex(hWarpMutex);
1236 :
1237 6 : CPLCond *hCond = CPLCreateCond();
1238 6 : CPLMutex *hCondMutex = CPLCreateMutex();
1239 6 : CPLReleaseMutex(hCondMutex);
1240 :
1241 : /* -------------------------------------------------------------------- */
1242 : /* Collect the list of chunks to operate on. */
1243 : /* -------------------------------------------------------------------- */
1244 6 : CollectChunkList(nDstXOff, nDstYOff, nDstXSize, nDstYSize);
1245 :
1246 : /* -------------------------------------------------------------------- */
1247 : /* Process them one at a time, updating the progress */
1248 : /* information for each region. */
1249 : /* -------------------------------------------------------------------- */
1250 6 : ChunkThreadData volatile asThreadData[2] = {};
1251 6 : CPLErrorAccumulator oErrorAccumulator;
1252 18 : for (int i = 0; i < 2; ++i)
1253 : {
1254 12 : asThreadData[i].poOperation = this;
1255 12 : asThreadData[i].hIOMutex = hIOMutex;
1256 12 : asThreadData[i].poErrorAccumulator = &oErrorAccumulator;
1257 : }
1258 :
1259 6 : double dfPixelsProcessed = 0.0;
1260 6 : double dfTotalPixels = static_cast<double>(nDstXSize) * nDstYSize;
1261 :
1262 6 : CPLErr eErr = CE_None;
1263 22 : for (int iChunk = 0; iChunk < nChunkListCount + 1; iChunk++)
1264 : {
1265 17 : int iThread = iChunk % 2;
1266 :
1267 : /* --------------------------------------------------------------------
1268 : */
1269 : /* Launch thread for this chunk. */
1270 : /* --------------------------------------------------------------------
1271 : */
1272 17 : if (pasChunkList != nullptr && iChunk < nChunkListCount)
1273 : {
1274 11 : GDALWarpChunk *pasThisChunk = pasChunkList + iChunk;
1275 11 : const double dfChunkPixels =
1276 11 : pasThisChunk->dsx * static_cast<double>(pasThisChunk->dsy);
1277 :
1278 11 : asThreadData[iThread].dfProgressBase =
1279 11 : dfPixelsProcessed / dfTotalPixels;
1280 11 : asThreadData[iThread].dfProgressScale =
1281 11 : dfChunkPixels / dfTotalPixels;
1282 :
1283 11 : dfPixelsProcessed += dfChunkPixels;
1284 :
1285 11 : asThreadData[iThread].pasChunkInfo = pasThisChunk;
1286 :
1287 11 : if (iChunk == 0)
1288 : {
1289 6 : asThreadData[iThread].hCond = hCond;
1290 6 : asThreadData[iThread].hCondMutex = hCondMutex;
1291 : }
1292 : else
1293 : {
1294 5 : asThreadData[iThread].hCond = nullptr;
1295 5 : asThreadData[iThread].hCondMutex = nullptr;
1296 : }
1297 11 : asThreadData[iThread].bIOMutexTaken = FALSE;
1298 :
1299 11 : CPLDebug("GDAL", "Start chunk %d / %d.", iChunk, nChunkListCount);
1300 22 : asThreadData[iThread].hThreadHandle = CPLCreateJoinableThread(
1301 : ChunkThreadMain,
1302 11 : const_cast<ChunkThreadData *>(&asThreadData[iThread]));
1303 11 : if (asThreadData[iThread].hThreadHandle == nullptr)
1304 : {
1305 0 : CPLError(
1306 : CE_Failure, CPLE_AppDefined,
1307 : "CPLCreateJoinableThread() failed in ChunkAndWarpMulti()");
1308 0 : eErr = CE_Failure;
1309 0 : break;
1310 : }
1311 :
1312 : // Wait that the first thread has acquired the IO mutex before
1313 : // proceeding. This will ensure that the first thread will run
1314 : // before the second one.
1315 11 : if (iChunk == 0)
1316 : {
1317 6 : CPLAcquireMutex(hCondMutex, 1.0);
1318 12 : while (asThreadData[iThread].bIOMutexTaken == FALSE)
1319 6 : CPLCondWait(hCond, hCondMutex);
1320 6 : CPLReleaseMutex(hCondMutex);
1321 : }
1322 : }
1323 :
1324 : /* --------------------------------------------------------------------
1325 : */
1326 : /* Wait for previous chunks thread to complete. */
1327 : /* --------------------------------------------------------------------
1328 : */
1329 17 : if (iChunk > 0)
1330 : {
1331 11 : iThread = (iChunk - 1) % 2;
1332 :
1333 : // Wait for thread to finish.
1334 11 : CPLJoinThread(asThreadData[iThread].hThreadHandle);
1335 11 : asThreadData[iThread].hThreadHandle = nullptr;
1336 :
1337 11 : CPLDebug("GDAL", "Finished chunk %d / %d.", iChunk - 1,
1338 : nChunkListCount);
1339 :
1340 11 : eErr = asThreadData[iThread].eErr;
1341 :
1342 11 : if (eErr != CE_None)
1343 1 : break;
1344 : }
1345 : }
1346 :
1347 : /* -------------------------------------------------------------------- */
1348 : /* Wait for all threads to complete. */
1349 : /* -------------------------------------------------------------------- */
1350 18 : for (int iThread = 0; iThread < 2; iThread++)
1351 : {
1352 12 : if (asThreadData[iThread].hThreadHandle)
1353 0 : CPLJoinThread(asThreadData[iThread].hThreadHandle);
1354 : }
1355 :
1356 6 : CPLDestroyCond(hCond);
1357 6 : CPLDestroyMutex(hCondMutex);
1358 :
1359 6 : WipeChunkList();
1360 :
1361 6 : oErrorAccumulator.ReplayErrors();
1362 :
1363 6 : psOptions->pfnProgress(1.0, "", psOptions->pProgressArg);
1364 :
1365 12 : return eErr;
1366 : }
1367 :
1368 : /************************************************************************/
1369 : /* GDALChunkAndWarpMulti() */
1370 : /************************************************************************/
1371 :
1372 : /**
1373 : * @see GDALWarpOperation::ChunkAndWarpMulti()
1374 : */
1375 :
1376 0 : CPLErr GDALChunkAndWarpMulti(GDALWarpOperationH hOperation, int nDstXOff,
1377 : int nDstYOff, int nDstXSize, int nDstYSize)
1378 : {
1379 0 : VALIDATE_POINTER1(hOperation, "GDALChunkAndWarpMulti", CE_Failure);
1380 :
1381 : return reinterpret_cast<GDALWarpOperation *>(hOperation)
1382 0 : ->ChunkAndWarpMulti(nDstXOff, nDstYOff, nDstXSize, nDstYSize);
1383 : }
1384 :
1385 : /************************************************************************/
1386 : /* WipeChunkList() */
1387 : /************************************************************************/
1388 :
1389 4277 : void GDALWarpOperation::WipeChunkList()
1390 :
1391 : {
1392 4277 : CPLFree(pasChunkList);
1393 4277 : pasChunkList = nullptr;
1394 4277 : nChunkListCount = 0;
1395 4277 : nChunkListMax = 0;
1396 4277 : }
1397 :
1398 : /************************************************************************/
1399 : /* GetWorkingMemoryForWindow() */
1400 : /************************************************************************/
1401 :
1402 : /** Returns the amount of working memory, in bytes, required to process
1403 : * a warped window of source dimensions nSrcXSize x nSrcYSize and target
1404 : * dimensions nDstXSize x nDstYSize.
1405 : */
1406 4113 : double GDALWarpOperation::GetWorkingMemoryForWindow(int nSrcXSize,
1407 : int nSrcYSize,
1408 : int nDstXSize,
1409 : int nDstYSize) const
1410 : {
1411 : /* -------------------------------------------------------------------- */
1412 : /* Based on the types of masks in use, how many bits will each */
1413 : /* source pixel cost us? */
1414 : /* -------------------------------------------------------------------- */
1415 : int nSrcPixelCostInBits =
1416 4113 : GDALGetDataTypeSizeBits(psOptions->eWorkingDataType) *
1417 4113 : psOptions->nBandCount;
1418 :
1419 4113 : if (psOptions->pfnSrcDensityMaskFunc != nullptr)
1420 0 : nSrcPixelCostInBits += 32; // Float mask?
1421 :
1422 4113 : GDALRasterBandH hSrcBand = nullptr;
1423 4113 : if (psOptions->nBandCount > 0)
1424 : hSrcBand =
1425 4113 : GDALGetRasterBand(psOptions->hSrcDS, psOptions->panSrcBands[0]);
1426 :
1427 4113 : if (psOptions->nSrcAlphaBand > 0 || psOptions->hCutline != nullptr)
1428 118 : nSrcPixelCostInBits += 32; // UnifiedSrcDensity float mask.
1429 7990 : else if (hSrcBand != nullptr &&
1430 3995 : (GDALGetMaskFlags(hSrcBand) & GMF_PER_DATASET))
1431 6 : nSrcPixelCostInBits += 1; // UnifiedSrcValid bit mask.
1432 :
1433 4113 : if (psOptions->papfnSrcPerBandValidityMaskFunc != nullptr ||
1434 4113 : psOptions->padfSrcNoDataReal != nullptr)
1435 211 : nSrcPixelCostInBits += psOptions->nBandCount; // Bit/band mask.
1436 :
1437 4113 : if (psOptions->pfnSrcValidityMaskFunc != nullptr)
1438 0 : nSrcPixelCostInBits += 1; // Bit mask.
1439 :
1440 : /* -------------------------------------------------------------------- */
1441 : /* What about the cost for the destination. */
1442 : /* -------------------------------------------------------------------- */
1443 : int nDstPixelCostInBits =
1444 4113 : GDALGetDataTypeSizeBits(psOptions->eWorkingDataType) *
1445 4113 : psOptions->nBandCount;
1446 :
1447 4113 : if (psOptions->pfnDstDensityMaskFunc != nullptr)
1448 0 : nDstPixelCostInBits += 32;
1449 :
1450 4113 : if (psOptions->padfDstNoDataReal != nullptr ||
1451 2435 : psOptions->pfnDstValidityMaskFunc != nullptr)
1452 1678 : nDstPixelCostInBits += psOptions->nBandCount;
1453 :
1454 4113 : if (psOptions->nDstAlphaBand > 0)
1455 256 : nDstPixelCostInBits += 32; // DstDensity float mask.
1456 :
1457 4113 : const double dfTotalMemoryUse =
1458 4113 : (static_cast<double>(nSrcPixelCostInBits) * nSrcXSize * nSrcYSize +
1459 4113 : static_cast<double>(nDstPixelCostInBits) * nDstXSize * nDstYSize) /
1460 : 8.0;
1461 4113 : return dfTotalMemoryUse;
1462 : }
1463 :
1464 : /************************************************************************/
1465 : /* CollectChunkListInternal() */
1466 : /************************************************************************/
1467 :
1468 4402 : CPLErr GDALWarpOperation::CollectChunkListInternal(int nDstXOff, int nDstYOff,
1469 : int nDstXSize, int nDstYSize)
1470 :
1471 : {
1472 : /* -------------------------------------------------------------------- */
1473 : /* Compute the bounds of the input area corresponding to the */
1474 : /* output area. */
1475 : /* -------------------------------------------------------------------- */
1476 4402 : int nSrcXOff = 0;
1477 4402 : int nSrcYOff = 0;
1478 4402 : int nSrcXSize = 0;
1479 4402 : int nSrcYSize = 0;
1480 4402 : double dfSrcXExtraSize = 0.0;
1481 4402 : double dfSrcYExtraSize = 0.0;
1482 4402 : double dfSrcFillRatio = 0.0;
1483 : CPLErr eErr;
1484 : {
1485 4402 : CPLTurnFailureIntoWarningBackuper oBackuper;
1486 4402 : eErr = ComputeSourceWindow(nDstXOff, nDstYOff, nDstXSize, nDstYSize,
1487 : &nSrcXOff, &nSrcYOff, &nSrcXSize, &nSrcYSize,
1488 : &dfSrcXExtraSize, &dfSrcYExtraSize,
1489 : &dfSrcFillRatio);
1490 : }
1491 :
1492 4402 : if (eErr != CE_None)
1493 : {
1494 : const bool bErrorOutIfEmptySourceWindow =
1495 3 : CPLFetchBool(psOptions->papszWarpOptions,
1496 : "ERROR_OUT_IF_EMPTY_SOURCE_WINDOW", true);
1497 3 : if (bErrorOutIfEmptySourceWindow)
1498 : {
1499 3 : CPLError(CE_Warning, CPLE_AppDefined,
1500 : "Unable to compute source region for "
1501 : "output window %d,%d,%d,%d, skipping.",
1502 : nDstXOff, nDstYOff, nDstXSize, nDstYSize);
1503 : }
1504 : else
1505 : {
1506 0 : CPLDebug("WARP",
1507 : "Unable to compute source region for "
1508 : "output window %d,%d,%d,%d, skipping.",
1509 : nDstXOff, nDstYOff, nDstXSize, nDstYSize);
1510 : }
1511 : }
1512 :
1513 : /* -------------------------------------------------------------------- */
1514 : /* If we are allowed to drop no-source regions, do so now if */
1515 : /* appropriate. */
1516 : /* -------------------------------------------------------------------- */
1517 5527 : if ((nSrcXSize == 0 || nSrcYSize == 0) &&
1518 1125 : CPLFetchBool(psOptions->papszWarpOptions, "SKIP_NOSOURCE", false))
1519 492 : return CE_None;
1520 :
1521 : /* -------------------------------------------------------------------- */
1522 : /* Does the cost of the current rectangle exceed our memory */
1523 : /* limit? If so, split the destination along the longest */
1524 : /* dimension and recurse. */
1525 : /* -------------------------------------------------------------------- */
1526 : const double dfTotalMemoryUse =
1527 3910 : GetWorkingMemoryForWindow(nSrcXSize, nSrcYSize, nDstXSize, nDstYSize);
1528 :
1529 : // If size of working buffers need exceed the allow limit, then divide
1530 : // the target area
1531 : // Do it also if the "fill ratio" of the source is too low (#3120), but
1532 : // only if there's at least some source pixel intersecting. The
1533 : // SRC_FILL_RATIO_HEURISTICS warping option is undocumented and only here
1534 : // in case the heuristics would cause issues.
1535 : #if DEBUG_VERBOSE
1536 : CPLDebug("WARP",
1537 : "dst=(%d,%d,%d,%d) src=(%d,%d,%d,%d) srcfillratio=%.17g, "
1538 : "dfTotalMemoryUse=%.1f MB",
1539 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, nSrcXOff, nSrcYOff,
1540 : nSrcXSize, nSrcYSize, dfSrcFillRatio,
1541 : dfTotalMemoryUse / (1024 * 1024));
1542 : #endif
1543 870 : if ((dfTotalMemoryUse > psOptions->dfWarpMemoryLimit &&
1544 7820 : (nDstXSize > 2 || nDstYSize > 2)) ||
1545 3040 : (dfSrcFillRatio > 0 && dfSrcFillRatio < 0.5 &&
1546 350 : (nDstXSize > 100 || nDstYSize > 100) &&
1547 720 : CPLFetchBool(psOptions->papszWarpOptions, "SRC_FILL_RATIO_HEURISTICS",
1548 : true)))
1549 : {
1550 1589 : int nBlockXSize = 1;
1551 1589 : int nBlockYSize = 1;
1552 1589 : if (psOptions->hDstDS)
1553 : {
1554 1589 : GDALGetBlockSize(GDALGetRasterBand(psOptions->hDstDS, 1),
1555 : &nBlockXSize, &nBlockYSize);
1556 : }
1557 :
1558 1589 : int bStreamableOutput = CPLFetchBool(psOptions->papszWarpOptions,
1559 1589 : "STREAMABLE_OUTPUT", false);
1560 : const char *pszOptimizeSize =
1561 1589 : CSLFetchNameValue(psOptions->papszWarpOptions, "OPTIMIZE_SIZE");
1562 1589 : const bool bOptimizeSizeAuto =
1563 1589 : !pszOptimizeSize || EQUAL(pszOptimizeSize, "AUTO");
1564 : const bool bOptimizeSize =
1565 4436 : !bStreamableOutput &&
1566 97 : ((pszOptimizeSize && !bOptimizeSizeAuto &&
1567 1588 : CPLTestBool(pszOptimizeSize)) ||
1568 : // Auto-enable optimize-size mode if output region is at least
1569 : // 2x2 blocks large and the shapes of the source and target regions
1570 : // are not excessively different. All those thresholds are a bit
1571 : // arbitrary
1572 1491 : (bOptimizeSizeAuto && nSrcXSize > 0 && nDstYSize > 0 &&
1573 1259 : (nDstXSize > nDstYSize ? fabs(double(nDstXSize) / nDstYSize -
1574 506 : double(nSrcXSize) / nSrcYSize) <
1575 506 : 5 * double(nDstXSize) / nDstYSize
1576 753 : : fabs(double(nDstYSize) / nDstXSize -
1577 753 : double(nSrcYSize) / nSrcXSize) <
1578 753 : 5 * double(nDstYSize) / nDstXSize) &&
1579 1256 : nDstXSize / 2 >= nBlockXSize && nDstYSize / 2 >= nBlockYSize));
1580 :
1581 : // If the region width is greater than the region height,
1582 : // cut in half in the width. When we want to optimize the size
1583 : // of a compressed output dataset, do this only if each half part
1584 : // is at least as wide as the block width.
1585 1589 : bool bHasDivided = false;
1586 1589 : CPLErr eErr2 = CE_None;
1587 1589 : if (nDstXSize > nDstYSize &&
1588 658 : ((!bOptimizeSize && !bStreamableOutput) ||
1589 88 : (bOptimizeSize &&
1590 89 : (nDstXSize / 2 >= nBlockXSize || nDstYSize == 1)) ||
1591 1 : (bStreamableOutput && nDstXSize / 2 >= nBlockXSize &&
1592 1 : nDstYSize == nBlockYSize)))
1593 : {
1594 611 : bHasDivided = true;
1595 611 : int nChunk1 = nDstXSize / 2;
1596 :
1597 : // In the optimize size case, try to stick on target block
1598 : // boundaries.
1599 611 : if ((bOptimizeSize || bStreamableOutput) && nChunk1 > nBlockXSize)
1600 42 : nChunk1 = (nChunk1 / nBlockXSize) * nBlockXSize;
1601 :
1602 611 : int nChunk2 = nDstXSize - nChunk1;
1603 :
1604 611 : eErr = CollectChunkListInternal(nDstXOff, nDstYOff, nChunk1,
1605 : nDstYSize);
1606 :
1607 611 : eErr2 = CollectChunkListInternal(nDstXOff + nChunk1, nDstYOff,
1608 611 : nChunk2, nDstYSize);
1609 : }
1610 978 : else if (!(bStreamableOutput && nDstYSize / 2 < nBlockYSize))
1611 : {
1612 977 : bHasDivided = true;
1613 977 : int nChunk1 = nDstYSize / 2;
1614 :
1615 : // In the optimize size case, try to stick on target block
1616 : // boundaries.
1617 977 : if ((bOptimizeSize || bStreamableOutput) && nChunk1 > nBlockYSize)
1618 77 : nChunk1 = (nChunk1 / nBlockYSize) * nBlockYSize;
1619 :
1620 977 : const int nChunk2 = nDstYSize - nChunk1;
1621 :
1622 977 : eErr = CollectChunkListInternal(nDstXOff, nDstYOff, nDstXSize,
1623 : nChunk1);
1624 :
1625 977 : eErr2 = CollectChunkListInternal(nDstXOff, nDstYOff + nChunk1,
1626 : nDstXSize, nChunk2);
1627 : }
1628 :
1629 1589 : if (bHasDivided)
1630 : {
1631 1588 : if (eErr == CE_None)
1632 1588 : return eErr2;
1633 : else
1634 0 : return eErr;
1635 : }
1636 : }
1637 :
1638 : /* -------------------------------------------------------------------- */
1639 : /* OK, everything fits, so add to the chunk list. */
1640 : /* -------------------------------------------------------------------- */
1641 2322 : if (nChunkListCount == nChunkListMax)
1642 : {
1643 1387 : nChunkListMax = nChunkListMax * 2 + 1;
1644 1387 : pasChunkList = static_cast<GDALWarpChunk *>(
1645 1387 : CPLRealloc(pasChunkList, sizeof(GDALWarpChunk) * nChunkListMax));
1646 : }
1647 :
1648 2322 : pasChunkList[nChunkListCount].dx = nDstXOff;
1649 2322 : pasChunkList[nChunkListCount].dy = nDstYOff;
1650 2322 : pasChunkList[nChunkListCount].dsx = nDstXSize;
1651 2322 : pasChunkList[nChunkListCount].dsy = nDstYSize;
1652 2322 : pasChunkList[nChunkListCount].sx = nSrcXOff;
1653 2322 : pasChunkList[nChunkListCount].sy = nSrcYOff;
1654 2322 : pasChunkList[nChunkListCount].ssx = nSrcXSize;
1655 2322 : pasChunkList[nChunkListCount].ssy = nSrcYSize;
1656 2322 : pasChunkList[nChunkListCount].sExtraSx = dfSrcXExtraSize;
1657 2322 : pasChunkList[nChunkListCount].sExtraSy = dfSrcYExtraSize;
1658 :
1659 2322 : nChunkListCount++;
1660 :
1661 2322 : return CE_None;
1662 : }
1663 :
1664 : /************************************************************************/
1665 : /* WarpRegion() */
1666 : /************************************************************************/
1667 :
1668 : /**
1669 : * This method requests the indicated region of the output file be generated.
1670 : *
1671 : * Note that WarpRegion() will produce the requested area in one low level warp
1672 : * operation without verifying that this does not exceed the stated memory
1673 : * limits for the warp operation. Applications should take care not to call
1674 : * WarpRegion() on too large a region! This function
1675 : * is normally called by ChunkAndWarpImage(), the normal entry point for
1676 : * applications. Use it instead if staying within memory constraints is
1677 : * desired.
1678 : *
1679 : * Progress is reported from dfProgressBase to dfProgressBase + dfProgressScale
1680 : * for the indicated region.
1681 : *
1682 : * @param nDstXOff X offset to window of destination data to be produced.
1683 : * @param nDstYOff Y offset to window of destination data to be produced.
1684 : * @param nDstXSize Width of output window on destination file to be produced.
1685 : * @param nDstYSize Height of output window on destination file to be produced.
1686 : * @param nSrcXOff source window X offset (computed if window all zero)
1687 : * @param nSrcYOff source window Y offset (computed if window all zero)
1688 : * @param nSrcXSize source window X size (computed if window all zero)
1689 : * @param nSrcYSize source window Y size (computed if window all zero)
1690 : * @param dfProgressBase minimum progress value reported
1691 : * @param dfProgressScale value such as dfProgressBase + dfProgressScale is the
1692 : * maximum progress value reported
1693 : *
1694 : * @return CE_None on success or CE_Failure if an error occurs.
1695 : */
1696 :
1697 0 : CPLErr GDALWarpOperation::WarpRegion(int nDstXOff, int nDstYOff, int nDstXSize,
1698 : int nDstYSize, int nSrcXOff, int nSrcYOff,
1699 : int nSrcXSize, int nSrcYSize,
1700 : double dfProgressBase,
1701 : double dfProgressScale)
1702 : {
1703 0 : return WarpRegion(nDstXOff, nDstYOff, nDstXSize, nDstYSize, nSrcXOff,
1704 : nSrcYOff, nSrcXSize, nSrcYSize, 0, 0, dfProgressBase,
1705 0 : dfProgressScale);
1706 : }
1707 :
1708 : /**
1709 : * This method requests the indicated region of the output file be generated.
1710 : *
1711 : * Note that WarpRegion() will produce the requested area in one low level warp
1712 : * operation without verifying that this does not exceed the stated memory
1713 : * limits for the warp operation. Applications should take care not to call
1714 : * WarpRegion() on too large a region! This function
1715 : * is normally called by ChunkAndWarpImage(), the normal entry point for
1716 : * applications. Use it instead if staying within memory constraints is
1717 : * desired.
1718 : *
1719 : * Progress is reported from dfProgressBase to dfProgressBase + dfProgressScale
1720 : * for the indicated region.
1721 : *
1722 : * @param nDstXOff X offset to window of destination data to be produced.
1723 : * @param nDstYOff Y offset to window of destination data to be produced.
1724 : * @param nDstXSize Width of output window on destination file to be produced.
1725 : * @param nDstYSize Height of output window on destination file to be produced.
1726 : * @param nSrcXOff source window X offset (computed if window all zero)
1727 : * @param nSrcYOff source window Y offset (computed if window all zero)
1728 : * @param nSrcXSize source window X size (computed if window all zero)
1729 : * @param nSrcYSize source window Y size (computed if window all zero)
1730 : * @param dfSrcXExtraSize Extra pixels (included in nSrcXSize) reserved
1731 : * for filter window. Should be ignored in scale computation
1732 : * @param dfSrcYExtraSize Extra pixels (included in nSrcYSize) reserved
1733 : * for filter window. Should be ignored in scale computation
1734 : * @param dfProgressBase minimum progress value reported
1735 : * @param dfProgressScale value such as dfProgressBase + dfProgressScale is the
1736 : * maximum progress value reported
1737 : *
1738 : * @return CE_None on success or CE_Failure if an error occurs.
1739 : */
1740 :
1741 2322 : CPLErr GDALWarpOperation::WarpRegion(
1742 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize, int nSrcXOff,
1743 : int nSrcYOff, int nSrcXSize, int nSrcYSize, double dfSrcXExtraSize,
1744 : double dfSrcYExtraSize, double dfProgressBase, double dfProgressScale)
1745 :
1746 : {
1747 2322 : ReportTiming(nullptr);
1748 :
1749 : /* -------------------------------------------------------------------- */
1750 : /* Allocate the output buffer. */
1751 : /* -------------------------------------------------------------------- */
1752 2322 : int bDstBufferInitialized = FALSE;
1753 : void *pDstBuffer =
1754 2322 : CreateDestinationBuffer(nDstXSize, nDstYSize, &bDstBufferInitialized);
1755 2322 : if (pDstBuffer == nullptr)
1756 : {
1757 2 : return CE_Failure;
1758 : }
1759 :
1760 : /* -------------------------------------------------------------------- */
1761 : /* If we aren't doing fixed initialization of the output buffer */
1762 : /* then read it from disk so we can overlay on existing imagery. */
1763 : /* -------------------------------------------------------------------- */
1764 2320 : GDALDataset *poDstDS = GDALDataset::FromHandle(psOptions->hDstDS);
1765 2320 : if (!bDstBufferInitialized)
1766 : {
1767 384 : CPLErr eErr = CE_None;
1768 384 : if (psOptions->nBandCount == 1)
1769 : {
1770 : // Particular case to simplify the stack a bit.
1771 : // TODO(rouault): Need an explanation of what and why r34502 helps.
1772 360 : eErr = poDstDS->GetRasterBand(psOptions->panDstBands[0])
1773 720 : ->RasterIO(GF_Read, nDstXOff, nDstYOff, nDstXSize,
1774 : nDstYSize, pDstBuffer, nDstXSize, nDstYSize,
1775 360 : psOptions->eWorkingDataType, 0, 0, nullptr);
1776 : }
1777 : else
1778 : {
1779 24 : eErr = poDstDS->RasterIO(GF_Read, nDstXOff, nDstYOff, nDstXSize,
1780 : nDstYSize, pDstBuffer, nDstXSize,
1781 24 : nDstYSize, psOptions->eWorkingDataType,
1782 24 : psOptions->nBandCount,
1783 24 : psOptions->panDstBands, 0, 0, 0, nullptr);
1784 : }
1785 :
1786 384 : if (eErr != CE_None)
1787 : {
1788 0 : DestroyDestinationBuffer(pDstBuffer);
1789 0 : return eErr;
1790 : }
1791 :
1792 384 : ReportTiming("Output buffer read");
1793 : }
1794 :
1795 : /* -------------------------------------------------------------------- */
1796 : /* Perform the warp. */
1797 : /* -------------------------------------------------------------------- */
1798 : CPLErr eErr = nSrcXSize == 0
1799 2320 : ? CE_None
1800 1919 : : WarpRegionToBuffer(
1801 : nDstXOff, nDstYOff, nDstXSize, nDstYSize,
1802 1919 : pDstBuffer, psOptions->eWorkingDataType, nSrcXOff,
1803 : nSrcYOff, nSrcXSize, nSrcYSize, dfSrcXExtraSize,
1804 2320 : dfSrcYExtraSize, dfProgressBase, dfProgressScale);
1805 :
1806 : /* -------------------------------------------------------------------- */
1807 : /* Write the output data back to disk if all went well. */
1808 : /* -------------------------------------------------------------------- */
1809 2320 : if (eErr == CE_None)
1810 : {
1811 2315 : if (psOptions->nBandCount == 1)
1812 : {
1813 : // Particular case to simplify the stack a bit.
1814 2101 : eErr = poDstDS->GetRasterBand(psOptions->panDstBands[0])
1815 4202 : ->RasterIO(GF_Write, nDstXOff, nDstYOff, nDstXSize,
1816 : nDstYSize, pDstBuffer, nDstXSize, nDstYSize,
1817 2101 : psOptions->eWorkingDataType, 0, 0, nullptr);
1818 : }
1819 : else
1820 : {
1821 214 : eErr = poDstDS->RasterIO(GF_Write, nDstXOff, nDstYOff, nDstXSize,
1822 : nDstYSize, pDstBuffer, nDstXSize,
1823 214 : nDstYSize, psOptions->eWorkingDataType,
1824 214 : psOptions->nBandCount,
1825 214 : psOptions->panDstBands, 0, 0, 0, nullptr);
1826 : }
1827 :
1828 4630 : if (eErr == CE_None &&
1829 2315 : CPLFetchBool(psOptions->papszWarpOptions, "WRITE_FLUSH", false))
1830 : {
1831 0 : const CPLErr eOldErr = CPLGetLastErrorType();
1832 0 : const CPLString osLastErrMsg = CPLGetLastErrorMsg();
1833 0 : GDALFlushCache(psOptions->hDstDS);
1834 0 : const CPLErr eNewErr = CPLGetLastErrorType();
1835 0 : if (eNewErr != eOldErr ||
1836 0 : osLastErrMsg.compare(CPLGetLastErrorMsg()) != 0)
1837 0 : eErr = CE_Failure;
1838 : }
1839 2315 : ReportTiming("Output buffer write");
1840 : }
1841 :
1842 : /* -------------------------------------------------------------------- */
1843 : /* Cleanup and return. */
1844 : /* -------------------------------------------------------------------- */
1845 2320 : DestroyDestinationBuffer(pDstBuffer);
1846 :
1847 2320 : return eErr;
1848 : }
1849 :
1850 : /************************************************************************/
1851 : /* GDALWarpRegion() */
1852 : /************************************************************************/
1853 :
1854 : /**
1855 : * @see GDALWarpOperation::WarpRegion()
1856 : */
1857 :
1858 0 : CPLErr GDALWarpRegion(GDALWarpOperationH hOperation, int nDstXOff, int nDstYOff,
1859 : int nDstXSize, int nDstYSize, int nSrcXOff, int nSrcYOff,
1860 : int nSrcXSize, int nSrcYSize)
1861 :
1862 : {
1863 0 : VALIDATE_POINTER1(hOperation, "GDALWarpRegion", CE_Failure);
1864 :
1865 : return reinterpret_cast<GDALWarpOperation *>(hOperation)
1866 0 : ->WarpRegion(nDstXOff, nDstYOff, nDstXSize, nDstYSize, nSrcXOff,
1867 0 : nSrcYOff, nSrcXSize, nSrcYSize);
1868 : }
1869 :
1870 : /************************************************************************/
1871 : /* WarpRegionToBuffer() */
1872 : /************************************************************************/
1873 :
1874 : /**
1875 : * This method requests that a particular window of the output dataset
1876 : * be warped and the result put into the provided data buffer. The output
1877 : * dataset doesn't even really have to exist to use this method as long as
1878 : * the transformation function in the GDALWarpOptions is setup to map to
1879 : * a virtual pixel/line space.
1880 : *
1881 : * This method will do the whole region in one chunk, so be wary of the
1882 : * amount of memory that might be used.
1883 : *
1884 : * @param nDstXOff X offset to window of destination data to be produced.
1885 : * @param nDstYOff Y offset to window of destination data to be produced.
1886 : * @param nDstXSize Width of output window on destination file to be produced.
1887 : * @param nDstYSize Height of output window on destination file to be produced.
1888 : * @param pDataBuf the data buffer to place result in, of type eBufDataType.
1889 : * @param eBufDataType the type of the output data buffer. For now this
1890 : * must match GDALWarpOptions::eWorkingDataType.
1891 : * @param nSrcXOff source window X offset (computed if window all zero)
1892 : * @param nSrcYOff source window Y offset (computed if window all zero)
1893 : * @param nSrcXSize source window X size (computed if window all zero)
1894 : * @param nSrcYSize source window Y size (computed if window all zero)
1895 : * @param dfProgressBase minimum progress value reported
1896 : * @param dfProgressScale value such as dfProgressBase + dfProgressScale is the
1897 : * maximum progress value reported
1898 : *
1899 : * @return CE_None on success or CE_Failure if an error occurs.
1900 : */
1901 :
1902 1977 : CPLErr GDALWarpOperation::WarpRegionToBuffer(
1903 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize, void *pDataBuf,
1904 : GDALDataType eBufDataType, int nSrcXOff, int nSrcYOff, int nSrcXSize,
1905 : int nSrcYSize, double dfProgressBase, double dfProgressScale)
1906 : {
1907 1977 : return WarpRegionToBuffer(nDstXOff, nDstYOff, nDstXSize, nDstYSize,
1908 : pDataBuf, eBufDataType, nSrcXOff, nSrcYOff,
1909 : nSrcXSize, nSrcYSize, 0, 0, dfProgressBase,
1910 1977 : dfProgressScale);
1911 : }
1912 :
1913 : /**
1914 : * This method requests that a particular window of the output dataset
1915 : * be warped and the result put into the provided data buffer. The output
1916 : * dataset doesn't even really have to exist to use this method as long as
1917 : * the transformation function in the GDALWarpOptions is setup to map to
1918 : * a virtual pixel/line space.
1919 : *
1920 : * This method will do the whole region in one chunk, so be wary of the
1921 : * amount of memory that might be used.
1922 : *
1923 : * @param nDstXOff X offset to window of destination data to be produced.
1924 : * @param nDstYOff Y offset to window of destination data to be produced.
1925 : * @param nDstXSize Width of output window on destination file to be produced.
1926 : * @param nDstYSize Height of output window on destination file to be produced.
1927 : * @param pDataBuf the data buffer to place result in, of type eBufDataType.
1928 : * @param eBufDataType the type of the output data buffer. For now this
1929 : * must match GDALWarpOptions::eWorkingDataType.
1930 : * @param nSrcXOff source window X offset (computed if window all zero)
1931 : * @param nSrcYOff source window Y offset (computed if window all zero)
1932 : * @param nSrcXSize source window X size (computed if window all zero)
1933 : * @param nSrcYSize source window Y size (computed if window all zero)
1934 : * @param dfSrcXExtraSize Extra pixels (included in nSrcXSize) reserved
1935 : * for filter window. Should be ignored in scale computation
1936 : * @param dfSrcYExtraSize Extra pixels (included in nSrcYSize) reserved
1937 : * for filter window. Should be ignored in scale computation
1938 : * @param dfProgressBase minimum progress value reported
1939 : * @param dfProgressScale value such as dfProgressBase + dfProgressScale is the
1940 : * maximum progress value reported
1941 : *
1942 : * @return CE_None on success or CE_Failure if an error occurs.
1943 : */
1944 :
1945 3896 : CPLErr GDALWarpOperation::WarpRegionToBuffer(
1946 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize, void *pDataBuf,
1947 : // Only in a CPLAssert.
1948 : CPL_UNUSED GDALDataType eBufDataType, int nSrcXOff, int nSrcYOff,
1949 : int nSrcXSize, int nSrcYSize, double dfSrcXExtraSize,
1950 : double dfSrcYExtraSize, double dfProgressBase, double dfProgressScale)
1951 :
1952 : {
1953 3896 : const int nWordSize = GDALGetDataTypeSizeBytes(psOptions->eWorkingDataType);
1954 :
1955 3896 : CPLAssert(eBufDataType == psOptions->eWorkingDataType);
1956 :
1957 : /* -------------------------------------------------------------------- */
1958 : /* If not given a corresponding source window compute one now. */
1959 : /* -------------------------------------------------------------------- */
1960 3896 : if (nSrcXSize == 0 && nSrcYSize == 0)
1961 : {
1962 : // TODO: This taking of the warp mutex is suboptimal. We could get rid
1963 : // of it, but that would require making sure ComputeSourceWindow()
1964 : // uses a different pTransformerArg than the warp kernel.
1965 1778 : if (hWarpMutex != nullptr && !CPLAcquireMutex(hWarpMutex, 600.0))
1966 : {
1967 0 : CPLError(CE_Failure, CPLE_AppDefined,
1968 : "Failed to acquire WarpMutex in WarpRegion().");
1969 0 : return CE_Failure;
1970 : }
1971 : const CPLErr eErr =
1972 1778 : ComputeSourceWindow(nDstXOff, nDstYOff, nDstXSize, nDstYSize,
1973 : &nSrcXOff, &nSrcYOff, &nSrcXSize, &nSrcYSize,
1974 : &dfSrcXExtraSize, &dfSrcYExtraSize, nullptr);
1975 1778 : if (hWarpMutex != nullptr)
1976 0 : CPLReleaseMutex(hWarpMutex);
1977 1778 : if (eErr != CE_None)
1978 : {
1979 : const bool bErrorOutIfEmptySourceWindow =
1980 36 : CPLFetchBool(psOptions->papszWarpOptions,
1981 : "ERROR_OUT_IF_EMPTY_SOURCE_WINDOW", true);
1982 36 : if (!bErrorOutIfEmptySourceWindow)
1983 36 : return CE_None;
1984 0 : return eErr;
1985 : }
1986 : }
1987 :
1988 : /* -------------------------------------------------------------------- */
1989 : /* Prepare a WarpKernel object to match this operation. */
1990 : /* -------------------------------------------------------------------- */
1991 7720 : GDALWarpKernel oWK;
1992 :
1993 3860 : oWK.eResample = m_bIsTranslationOnPixelBoundaries ? GRA_NearestNeighbour
1994 3545 : : psOptions->eResampleAlg;
1995 3860 : oWK.eTieStrategy = psOptions->eTieStrategy;
1996 3860 : oWK.nBands = psOptions->nBandCount;
1997 3860 : oWK.eWorkingDataType = psOptions->eWorkingDataType;
1998 :
1999 3860 : oWK.pfnTransformer = psOptions->pfnTransformer;
2000 3860 : oWK.pTransformerArg = psOptions->pTransformerArg;
2001 :
2002 3860 : oWK.pfnProgress = psOptions->pfnProgress;
2003 3860 : oWK.pProgress = psOptions->pProgressArg;
2004 3860 : oWK.dfProgressBase = dfProgressBase;
2005 3860 : oWK.dfProgressScale = dfProgressScale;
2006 :
2007 3860 : oWK.papszWarpOptions = psOptions->papszWarpOptions;
2008 3860 : oWK.psThreadData = psThreadData;
2009 :
2010 3860 : oWK.padfDstNoDataReal = psOptions->padfDstNoDataReal;
2011 :
2012 : /* -------------------------------------------------------------------- */
2013 : /* Setup the source buffer. */
2014 : /* */
2015 : /* Eventually we may need to take advantage of pixel */
2016 : /* interleaved reading here. */
2017 : /* -------------------------------------------------------------------- */
2018 3860 : oWK.nSrcXOff = nSrcXOff;
2019 3860 : oWK.nSrcYOff = nSrcYOff;
2020 3860 : oWK.nSrcXSize = nSrcXSize;
2021 3860 : oWK.nSrcYSize = nSrcYSize;
2022 3860 : oWK.dfSrcXExtraSize = dfSrcXExtraSize;
2023 3860 : oWK.dfSrcYExtraSize = dfSrcYExtraSize;
2024 :
2025 : // Check for overflows in computation of nAlloc
2026 7117 : if (nSrcYSize > 0 &&
2027 3257 : ((static_cast<size_t>(nSrcXSize) >
2028 3257 : (std::numeric_limits<size_t>::max() - WARP_EXTRA_ELTS) / nSrcYSize) ||
2029 3257 : (static_cast<size_t>(nSrcXSize) * nSrcYSize + WARP_EXTRA_ELTS >
2030 3257 : std::numeric_limits<size_t>::max() /
2031 3257 : (nWordSize * psOptions->nBandCount))))
2032 : {
2033 0 : CPLError(CE_Failure, CPLE_AppDefined,
2034 : "WarpRegionToBuffer(): Integer overflow : nWordSize(=%d) * "
2035 : "(nSrcXSize(=%d) * nSrcYSize(=%d) + WARP_EXTRA_ELTS(=%d)) * "
2036 : "nBandCount(=%d)",
2037 : nWordSize, nSrcXSize, nSrcYSize, WARP_EXTRA_ELTS,
2038 0 : psOptions->nBandCount);
2039 0 : return CE_Failure;
2040 : }
2041 :
2042 3860 : const size_t nAlloc =
2043 3860 : nWordSize *
2044 3860 : (static_cast<size_t>(nSrcXSize) * nSrcYSize + WARP_EXTRA_ELTS) *
2045 3860 : psOptions->nBandCount;
2046 :
2047 3860 : oWK.papabySrcImage = static_cast<GByte **>(
2048 3860 : CPLCalloc(sizeof(GByte *), psOptions->nBandCount));
2049 3860 : oWK.papabySrcImage[0] = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nAlloc));
2050 :
2051 3860 : CPLErr eErr =
2052 3258 : nSrcXSize != 0 && nSrcYSize != 0 && oWK.papabySrcImage[0] == nullptr
2053 7118 : ? CE_Failure
2054 : : CE_None;
2055 :
2056 11313 : for (int i = 0; i < psOptions->nBandCount && eErr == CE_None; i++)
2057 7453 : oWK.papabySrcImage[i] =
2058 7453 : reinterpret_cast<GByte *>(oWK.papabySrcImage[0]) +
2059 7453 : nWordSize *
2060 7453 : (static_cast<GPtrDiff_t>(nSrcXSize) * nSrcYSize +
2061 7453 : WARP_EXTRA_ELTS) *
2062 7453 : i;
2063 :
2064 3860 : if (eErr == CE_None && nSrcXSize > 0 && nSrcYSize > 0)
2065 : {
2066 3247 : GDALDataset *poSrcDS = GDALDataset::FromHandle(psOptions->hSrcDS);
2067 3247 : if (psOptions->nBandCount == 1)
2068 : {
2069 : // Particular case to simplify the stack a bit.
2070 2023 : eErr = poSrcDS->GetRasterBand(psOptions->panSrcBands[0])
2071 4046 : ->RasterIO(GF_Read, nSrcXOff, nSrcYOff, nSrcXSize,
2072 2023 : nSrcYSize, oWK.papabySrcImage[0], nSrcXSize,
2073 2023 : nSrcYSize, psOptions->eWorkingDataType, 0, 0,
2074 : nullptr);
2075 : }
2076 : else
2077 : {
2078 1224 : eErr = poSrcDS->RasterIO(
2079 : GF_Read, nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
2080 1224 : oWK.papabySrcImage[0], nSrcXSize, nSrcYSize,
2081 1224 : psOptions->eWorkingDataType, psOptions->nBandCount,
2082 1224 : psOptions->panSrcBands, 0, 0,
2083 1224 : nWordSize * (static_cast<GPtrDiff_t>(nSrcXSize) * nSrcYSize +
2084 : WARP_EXTRA_ELTS),
2085 : nullptr);
2086 : }
2087 : }
2088 :
2089 3860 : ReportTiming("Input buffer read");
2090 :
2091 : /* -------------------------------------------------------------------- */
2092 : /* Initialize destination buffer. */
2093 : /* -------------------------------------------------------------------- */
2094 3860 : oWK.nDstXOff = nDstXOff;
2095 3860 : oWK.nDstYOff = nDstYOff;
2096 3860 : oWK.nDstXSize = nDstXSize;
2097 3860 : oWK.nDstYSize = nDstYSize;
2098 :
2099 3860 : oWK.papabyDstImage = reinterpret_cast<GByte **>(
2100 3860 : CPLCalloc(sizeof(GByte *), psOptions->nBandCount));
2101 :
2102 11305 : for (int i = 0; i < psOptions->nBandCount && eErr == CE_None; i++)
2103 : {
2104 7445 : oWK.papabyDstImage[i] =
2105 7445 : static_cast<GByte *>(pDataBuf) +
2106 7445 : i * static_cast<GPtrDiff_t>(nDstXSize) * nDstYSize * nWordSize;
2107 : }
2108 :
2109 : /* -------------------------------------------------------------------- */
2110 : /* Eventually we need handling for a whole bunch of the */
2111 : /* validity and density masks here. */
2112 : /* -------------------------------------------------------------------- */
2113 :
2114 : // TODO
2115 :
2116 : /* -------------------------------------------------------------------- */
2117 : /* Generate a source density mask if we have a source alpha band */
2118 : /* -------------------------------------------------------------------- */
2119 3860 : if (eErr == CE_None && psOptions->nSrcAlphaBand > 0 && nSrcXSize > 0 &&
2120 142 : nSrcYSize > 0)
2121 : {
2122 142 : CPLAssert(oWK.pafUnifiedSrcDensity == nullptr);
2123 :
2124 142 : eErr = CreateKernelMask(&oWK, 0 /* not used */, "UnifiedSrcDensity");
2125 :
2126 142 : if (eErr == CE_None)
2127 : {
2128 142 : int bOutAllOpaque = FALSE;
2129 284 : eErr = GDALWarpSrcAlphaMasker(
2130 142 : psOptions, psOptions->nBandCount, psOptions->eWorkingDataType,
2131 : oWK.nSrcXOff, oWK.nSrcYOff, oWK.nSrcXSize, oWK.nSrcYSize,
2132 142 : oWK.papabySrcImage, TRUE, oWK.pafUnifiedSrcDensity,
2133 : &bOutAllOpaque);
2134 142 : if (bOutAllOpaque)
2135 : {
2136 : #if DEBUG_VERBOSE
2137 : CPLDebug("WARP",
2138 : "No need for a source density mask as all values "
2139 : "are opaque");
2140 : #endif
2141 35 : CPLFree(oWK.pafUnifiedSrcDensity);
2142 35 : oWK.pafUnifiedSrcDensity = nullptr;
2143 : }
2144 : }
2145 : }
2146 :
2147 : /* -------------------------------------------------------------------- */
2148 : /* Generate a source density mask if we have a source cutline. */
2149 : /* -------------------------------------------------------------------- */
2150 3860 : if (eErr == CE_None && psOptions->hCutline != nullptr && nSrcXSize > 0 &&
2151 44 : nSrcYSize > 0)
2152 : {
2153 44 : const bool bUnifiedSrcDensityJustCreated =
2154 44 : (oWK.pafUnifiedSrcDensity == nullptr);
2155 44 : if (bUnifiedSrcDensityJustCreated)
2156 : {
2157 : eErr =
2158 44 : CreateKernelMask(&oWK, 0 /* not used */, "UnifiedSrcDensity");
2159 :
2160 44 : if (eErr == CE_None)
2161 : {
2162 44 : for (GPtrDiff_t j = 0;
2163 2929360 : j < static_cast<GPtrDiff_t>(oWK.nSrcXSize) * oWK.nSrcYSize;
2164 : j++)
2165 2929320 : oWK.pafUnifiedSrcDensity[j] = 1.0;
2166 : }
2167 : }
2168 :
2169 44 : int nValidityFlag = 0;
2170 44 : if (eErr == CE_None)
2171 44 : eErr = GDALWarpCutlineMaskerEx(
2172 44 : psOptions, psOptions->nBandCount, psOptions->eWorkingDataType,
2173 : oWK.nSrcXOff, oWK.nSrcYOff, oWK.nSrcXSize, oWK.nSrcYSize,
2174 44 : oWK.papabySrcImage, TRUE, oWK.pafUnifiedSrcDensity,
2175 : &nValidityFlag);
2176 44 : if (nValidityFlag == GCMVF_CHUNK_FULLY_WITHIN_CUTLINE &&
2177 : bUnifiedSrcDensityJustCreated)
2178 : {
2179 8 : VSIFree(oWK.pafUnifiedSrcDensity);
2180 8 : oWK.pafUnifiedSrcDensity = nullptr;
2181 : }
2182 : }
2183 :
2184 : /* -------------------------------------------------------------------- */
2185 : /* Generate a destination density mask if we have a destination */
2186 : /* alpha band. */
2187 : /* -------------------------------------------------------------------- */
2188 3860 : if (eErr == CE_None && psOptions->nDstAlphaBand > 0)
2189 : {
2190 1749 : CPLAssert(oWK.pafDstDensity == nullptr);
2191 :
2192 1749 : eErr = CreateKernelMask(&oWK, 0 /* not used */, "DstDensity");
2193 :
2194 1749 : if (eErr == CE_None)
2195 1749 : eErr = GDALWarpDstAlphaMasker(
2196 1749 : psOptions, psOptions->nBandCount, psOptions->eWorkingDataType,
2197 : oWK.nDstXOff, oWK.nDstYOff, oWK.nDstXSize, oWK.nDstYSize,
2198 1749 : oWK.papabyDstImage, TRUE, oWK.pafDstDensity);
2199 : }
2200 :
2201 : /* -------------------------------------------------------------------- */
2202 : /* If we have source nodata values create the validity mask. */
2203 : /* -------------------------------------------------------------------- */
2204 3860 : if (eErr == CE_None && psOptions->padfSrcNoDataReal != nullptr &&
2205 232 : nSrcXSize > 0 && nSrcYSize > 0)
2206 : {
2207 214 : CPLAssert(oWK.papanBandSrcValid == nullptr);
2208 :
2209 214 : bool bAllBandsAllValid = true;
2210 539 : for (int i = 0; i < psOptions->nBandCount && eErr == CE_None; i++)
2211 : {
2212 325 : eErr = CreateKernelMask(&oWK, i, "BandSrcValid");
2213 325 : if (eErr == CE_None)
2214 : {
2215 325 : double adfNoData[2] = {psOptions->padfSrcNoDataReal[i],
2216 325 : psOptions->padfSrcNoDataImag != nullptr
2217 325 : ? psOptions->padfSrcNoDataImag[i]
2218 325 : : 0.0};
2219 :
2220 325 : int bAllValid = FALSE;
2221 650 : eErr = GDALWarpNoDataMasker(
2222 325 : adfNoData, 1, psOptions->eWorkingDataType, oWK.nSrcXOff,
2223 : oWK.nSrcYOff, oWK.nSrcXSize, oWK.nSrcYSize,
2224 325 : &(oWK.papabySrcImage[i]), FALSE, oWK.papanBandSrcValid[i],
2225 : &bAllValid);
2226 325 : if (!bAllValid)
2227 188 : bAllBandsAllValid = false;
2228 : }
2229 : }
2230 :
2231 : // Optimization: if all pixels in all bands are valid,
2232 : // we don't need a mask.
2233 214 : if (bAllBandsAllValid)
2234 : {
2235 : #if DEBUG_VERBOSE
2236 : CPLDebug(
2237 : "WARP",
2238 : "No need for a source nodata mask as all values are valid");
2239 : #endif
2240 231 : for (int k = 0; k < psOptions->nBandCount; k++)
2241 136 : CPLFree(oWK.papanBandSrcValid[k]);
2242 95 : CPLFree(oWK.papanBandSrcValid);
2243 95 : oWK.papanBandSrcValid = nullptr;
2244 : }
2245 :
2246 : /* --------------------------------------------------------------------
2247 : */
2248 : /* If there's just a single band, then transfer */
2249 : /* papanBandSrcValid[0] as panUnifiedSrcValid. */
2250 : /* --------------------------------------------------------------------
2251 : */
2252 214 : if (oWK.papanBandSrcValid != nullptr && psOptions->nBandCount == 1)
2253 : {
2254 73 : oWK.panUnifiedSrcValid = oWK.papanBandSrcValid[0];
2255 73 : CPLFree(oWK.papanBandSrcValid);
2256 73 : oWK.papanBandSrcValid = nullptr;
2257 : }
2258 :
2259 : /* --------------------------------------------------------------------
2260 : */
2261 : /* Compute a unified input pixel mask if and only if all bands */
2262 : /* nodata is true. That is, we only treat a pixel as nodata if */
2263 : /* all bands match their respective nodata values. */
2264 : /* --------------------------------------------------------------------
2265 : */
2266 141 : else if (oWK.papanBandSrcValid != nullptr && eErr == CE_None)
2267 : {
2268 46 : bool bAtLeastOneBandAllValid = false;
2269 162 : for (int k = 0; k < psOptions->nBandCount; k++)
2270 : {
2271 116 : if (oWK.papanBandSrcValid[k] == nullptr)
2272 : {
2273 0 : bAtLeastOneBandAllValid = true;
2274 0 : break;
2275 : }
2276 : }
2277 :
2278 92 : const char *pszUnifiedSrcNoData = CSLFetchNameValue(
2279 46 : psOptions->papszWarpOptions, "UNIFIED_SRC_NODATA");
2280 84 : if (!bAtLeastOneBandAllValid && (pszUnifiedSrcNoData == nullptr ||
2281 38 : CPLTestBool(pszUnifiedSrcNoData)))
2282 : {
2283 45 : auto nMaskBits =
2284 45 : static_cast<GPtrDiff_t>(oWK.nSrcXSize) * oWK.nSrcYSize;
2285 :
2286 : eErr =
2287 45 : CreateKernelMask(&oWK, 0 /* not used */, "UnifiedSrcValid");
2288 :
2289 45 : if (eErr == CE_None)
2290 : {
2291 45 : CPLMaskClearAll(oWK.panUnifiedSrcValid, nMaskBits);
2292 :
2293 158 : for (int k = 0; k < psOptions->nBandCount; k++)
2294 : {
2295 113 : CPLMaskMerge(oWK.panUnifiedSrcValid,
2296 113 : oWK.papanBandSrcValid[k], nMaskBits);
2297 : }
2298 :
2299 : // If UNIFIED_SRC_NODATA is set, then we will ignore the
2300 : // individual nodata status of each band. If it is not set,
2301 : // both mechanism apply:
2302 : // - if panUnifiedSrcValid[] indicates a pixel is invalid
2303 : // (that is all its bands are at nodata), then the output
2304 : // pixel will be invalid
2305 : // - otherwise, the status band per band will be check with
2306 : // papanBandSrcValid[iBand][], and the output pixel will
2307 : // be valid
2308 45 : if (pszUnifiedSrcNoData != nullptr &&
2309 37 : !EQUAL(pszUnifiedSrcNoData, "PARTIAL"))
2310 : {
2311 123 : for (int k = 0; k < psOptions->nBandCount; k++)
2312 87 : CPLFree(oWK.papanBandSrcValid[k]);
2313 36 : CPLFree(oWK.papanBandSrcValid);
2314 36 : oWK.papanBandSrcValid = nullptr;
2315 : }
2316 : }
2317 : }
2318 : }
2319 : }
2320 :
2321 : /* -------------------------------------------------------------------- */
2322 : /* Generate a source validity mask if we have a source mask for */
2323 : /* the whole input dataset (and didn't already treat it as */
2324 : /* alpha band). */
2325 : /* -------------------------------------------------------------------- */
2326 : GDALRasterBandH hSrcBand =
2327 3860 : psOptions->nBandCount < 1
2328 3860 : ? nullptr
2329 3860 : : GDALGetRasterBand(psOptions->hSrcDS, psOptions->panSrcBands[0]);
2330 :
2331 3856 : if (eErr == CE_None && oWK.pafUnifiedSrcDensity == nullptr &&
2332 3713 : oWK.panUnifiedSrcValid == nullptr && psOptions->nSrcAlphaBand <= 0 &&
2333 3532 : (GDALGetMaskFlags(hSrcBand) & GMF_PER_DATASET)
2334 : // Need to double check for -nosrcalpha case.
2335 5 : && !(GDALGetMaskFlags(hSrcBand) & GMF_ALPHA) &&
2336 7718 : psOptions->padfSrcNoDataReal == nullptr && nSrcXSize > 0 &&
2337 2 : nSrcYSize > 0)
2338 :
2339 : {
2340 2 : eErr = CreateKernelMask(&oWK, 0 /* not used */, "UnifiedSrcValid");
2341 :
2342 2 : if (eErr == CE_None)
2343 2 : eErr = GDALWarpSrcMaskMasker(
2344 2 : psOptions, psOptions->nBandCount, psOptions->eWorkingDataType,
2345 : oWK.nSrcXOff, oWK.nSrcYOff, oWK.nSrcXSize, oWK.nSrcYSize,
2346 2 : oWK.papabySrcImage, FALSE, oWK.panUnifiedSrcValid);
2347 : }
2348 :
2349 : /* -------------------------------------------------------------------- */
2350 : /* If we have destination nodata values create the */
2351 : /* validity mask. We set the DstValid for any pixel that we */
2352 : /* do no have valid data in *any* of the source bands. */
2353 : /* */
2354 : /* Note that we don't support any concept of unified nodata on */
2355 : /* the destination image. At some point that should be added */
2356 : /* and then this logic will be significantly different. */
2357 : /* -------------------------------------------------------------------- */
2358 3860 : if (eErr == CE_None && psOptions->padfDstNoDataReal != nullptr)
2359 : {
2360 587 : CPLAssert(oWK.panDstValid == nullptr);
2361 :
2362 587 : const GPtrDiff_t nMaskBits =
2363 587 : static_cast<GPtrDiff_t>(oWK.nDstXSize) * oWK.nDstYSize;
2364 :
2365 587 : eErr = CreateKernelMask(&oWK, 0 /* not used */, "DstValid");
2366 : GUInt32 *panBandMask =
2367 587 : eErr == CE_None ? CPLMaskCreate(nMaskBits, true) : nullptr;
2368 :
2369 587 : if (eErr == CE_None && panBandMask != nullptr)
2370 : {
2371 1241 : for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
2372 : {
2373 712 : CPLMaskSetAll(panBandMask, nMaskBits);
2374 :
2375 712 : double adfNoData[2] = {psOptions->padfDstNoDataReal[iBand],
2376 712 : psOptions->padfDstNoDataImag != nullptr
2377 712 : ? psOptions->padfDstNoDataImag[iBand]
2378 712 : : 0.0};
2379 :
2380 712 : int bAllValid = FALSE;
2381 1424 : eErr = GDALWarpNoDataMasker(
2382 712 : adfNoData, 1, psOptions->eWorkingDataType, oWK.nDstXOff,
2383 : oWK.nDstYOff, oWK.nDstXSize, oWK.nDstYSize,
2384 712 : oWK.papabyDstImage + iBand, FALSE, panBandMask, &bAllValid);
2385 :
2386 : // Optimization: if there's a single band and all pixels are
2387 : // valid then we don't need a mask.
2388 712 : if (bAllValid && psOptions->nBandCount == 1)
2389 : {
2390 : #if DEBUG_VERBOSE
2391 : CPLDebug("WARP", "No need for a destination nodata mask as "
2392 : "all values are valid");
2393 : #endif
2394 58 : CPLFree(oWK.panDstValid);
2395 58 : oWK.panDstValid = nullptr;
2396 58 : break;
2397 : }
2398 :
2399 654 : CPLMaskMerge(oWK.panDstValid, panBandMask, nMaskBits);
2400 : }
2401 587 : CPLFree(panBandMask);
2402 : }
2403 : }
2404 :
2405 : /* -------------------------------------------------------------------- */
2406 : /* Release IO Mutex, and acquire warper mutex. */
2407 : /* -------------------------------------------------------------------- */
2408 3860 : if (hIOMutex != nullptr)
2409 : {
2410 11 : CPLReleaseMutex(hIOMutex);
2411 11 : if (!CPLAcquireMutex(hWarpMutex, 600.0))
2412 : {
2413 0 : CPLError(CE_Failure, CPLE_AppDefined,
2414 : "Failed to acquire WarpMutex in WarpRegion().");
2415 0 : return CE_Failure;
2416 : }
2417 : }
2418 :
2419 : /* -------------------------------------------------------------------- */
2420 : /* Optional application provided prewarp chunk processor. */
2421 : /* -------------------------------------------------------------------- */
2422 3860 : if (eErr == CE_None && psOptions->pfnPreWarpChunkProcessor != nullptr)
2423 0 : eErr = psOptions->pfnPreWarpChunkProcessor(
2424 0 : &oWK, psOptions->pPreWarpProcessorArg);
2425 :
2426 : /* -------------------------------------------------------------------- */
2427 : /* Perform the warp. */
2428 : /* -------------------------------------------------------------------- */
2429 3860 : if (eErr == CE_None)
2430 : {
2431 3856 : eErr = oWK.PerformWarp();
2432 3856 : ReportTiming("In memory warp operation");
2433 : }
2434 :
2435 : /* -------------------------------------------------------------------- */
2436 : /* Optional application provided postwarp chunk processor. */
2437 : /* -------------------------------------------------------------------- */
2438 3860 : if (eErr == CE_None && psOptions->pfnPostWarpChunkProcessor != nullptr)
2439 0 : eErr = psOptions->pfnPostWarpChunkProcessor(
2440 0 : &oWK, psOptions->pPostWarpProcessorArg);
2441 :
2442 : /* -------------------------------------------------------------------- */
2443 : /* Release Warp Mutex, and acquire io mutex. */
2444 : /* -------------------------------------------------------------------- */
2445 3860 : if (hIOMutex != nullptr)
2446 : {
2447 11 : CPLReleaseMutex(hWarpMutex);
2448 11 : if (!CPLAcquireMutex(hIOMutex, 600.0))
2449 : {
2450 0 : CPLError(CE_Failure, CPLE_AppDefined,
2451 : "Failed to acquire IOMutex in WarpRegion().");
2452 0 : return CE_Failure;
2453 : }
2454 : }
2455 :
2456 : /* -------------------------------------------------------------------- */
2457 : /* Write destination alpha if available. */
2458 : /* -------------------------------------------------------------------- */
2459 3860 : if (eErr == CE_None && psOptions->nDstAlphaBand > 0)
2460 : {
2461 1749 : eErr = GDALWarpDstAlphaMasker(
2462 1749 : psOptions, -psOptions->nBandCount, psOptions->eWorkingDataType,
2463 : oWK.nDstXOff, oWK.nDstYOff, oWK.nDstXSize, oWK.nDstYSize,
2464 1749 : oWK.papabyDstImage, TRUE, oWK.pafDstDensity);
2465 : }
2466 :
2467 : /* -------------------------------------------------------------------- */
2468 : /* Cleanup. */
2469 : /* -------------------------------------------------------------------- */
2470 3860 : CPLFree(oWK.papabySrcImage[0]);
2471 3860 : CPLFree(oWK.papabySrcImage);
2472 3860 : CPLFree(oWK.papabyDstImage);
2473 :
2474 3860 : if (oWK.papanBandSrcValid != nullptr)
2475 : {
2476 39 : for (int i = 0; i < oWK.nBands; i++)
2477 29 : CPLFree(oWK.papanBandSrcValid[i]);
2478 10 : CPLFree(oWK.papanBandSrcValid);
2479 : }
2480 3860 : CPLFree(oWK.panUnifiedSrcValid);
2481 3860 : CPLFree(oWK.pafUnifiedSrcDensity);
2482 3860 : CPLFree(oWK.panDstValid);
2483 3860 : CPLFree(oWK.pafDstDensity);
2484 :
2485 3860 : return eErr;
2486 : }
2487 :
2488 : /************************************************************************/
2489 : /* GDALWarpRegionToBuffer() */
2490 : /************************************************************************/
2491 :
2492 : /**
2493 : * @see GDALWarpOperation::WarpRegionToBuffer()
2494 : */
2495 :
2496 0 : CPLErr GDALWarpRegionToBuffer(GDALWarpOperationH hOperation, int nDstXOff,
2497 : int nDstYOff, int nDstXSize, int nDstYSize,
2498 : void *pDataBuf, GDALDataType eBufDataType,
2499 : int nSrcXOff, int nSrcYOff, int nSrcXSize,
2500 : int nSrcYSize)
2501 :
2502 : {
2503 0 : VALIDATE_POINTER1(hOperation, "GDALWarpRegionToBuffer", CE_Failure);
2504 :
2505 : return reinterpret_cast<GDALWarpOperation *>(hOperation)
2506 0 : ->WarpRegionToBuffer(nDstXOff, nDstYOff, nDstXSize, nDstYSize, pDataBuf,
2507 : eBufDataType, nSrcXOff, nSrcYOff, nSrcXSize,
2508 0 : nSrcYSize);
2509 : }
2510 :
2511 : /************************************************************************/
2512 : /* CreateKernelMask() */
2513 : /* */
2514 : /* If mask does not yet exist, create it. Supported types are */
2515 : /* the name of the variable in question. That is */
2516 : /* "BandSrcValid", "UnifiedSrcValid", "UnifiedSrcDensity", */
2517 : /* "DstValid", and "DstDensity". */
2518 : /************************************************************************/
2519 :
2520 2894 : CPLErr GDALWarpOperation::CreateKernelMask(GDALWarpKernel *poKernel, int iBand,
2521 : const char *pszType)
2522 :
2523 : {
2524 2894 : void **ppMask = nullptr;
2525 2894 : int nXSize = 0;
2526 2894 : int nYSize = 0;
2527 2894 : int nBitsPerPixel = 0;
2528 2894 : int nDefault = 0;
2529 2894 : int nExtraElts = 0;
2530 2894 : bool bDoMemset = true;
2531 :
2532 : /* -------------------------------------------------------------------- */
2533 : /* Get particulars of mask to be updated. */
2534 : /* -------------------------------------------------------------------- */
2535 2894 : if (EQUAL(pszType, "BandSrcValid"))
2536 : {
2537 325 : if (poKernel->papanBandSrcValid == nullptr)
2538 214 : poKernel->papanBandSrcValid = static_cast<GUInt32 **>(
2539 214 : CPLCalloc(sizeof(void *), poKernel->nBands));
2540 :
2541 325 : ppMask =
2542 325 : reinterpret_cast<void **>(&(poKernel->papanBandSrcValid[iBand]));
2543 325 : nExtraElts = WARP_EXTRA_ELTS;
2544 325 : nXSize = poKernel->nSrcXSize;
2545 325 : nYSize = poKernel->nSrcYSize;
2546 325 : nBitsPerPixel = 1;
2547 325 : nDefault = 0xff;
2548 : }
2549 2569 : else if (EQUAL(pszType, "UnifiedSrcValid"))
2550 : {
2551 47 : ppMask = reinterpret_cast<void **>(&(poKernel->panUnifiedSrcValid));
2552 47 : nExtraElts = WARP_EXTRA_ELTS;
2553 47 : nXSize = poKernel->nSrcXSize;
2554 47 : nYSize = poKernel->nSrcYSize;
2555 47 : nBitsPerPixel = 1;
2556 47 : nDefault = 0xff;
2557 : }
2558 2522 : else if (EQUAL(pszType, "UnifiedSrcDensity"))
2559 : {
2560 186 : ppMask = reinterpret_cast<void **>(&(poKernel->pafUnifiedSrcDensity));
2561 186 : nExtraElts = WARP_EXTRA_ELTS;
2562 186 : nXSize = poKernel->nSrcXSize;
2563 186 : nYSize = poKernel->nSrcYSize;
2564 186 : nBitsPerPixel = 32;
2565 186 : nDefault = 0;
2566 186 : bDoMemset = false;
2567 : }
2568 2336 : else if (EQUAL(pszType, "DstValid"))
2569 : {
2570 587 : ppMask = reinterpret_cast<void **>(&(poKernel->panDstValid));
2571 587 : nXSize = poKernel->nDstXSize;
2572 587 : nYSize = poKernel->nDstYSize;
2573 587 : nBitsPerPixel = 1;
2574 587 : nDefault = 0;
2575 : }
2576 1749 : else if (EQUAL(pszType, "DstDensity"))
2577 : {
2578 1749 : ppMask = reinterpret_cast<void **>(&(poKernel->pafDstDensity));
2579 1749 : nXSize = poKernel->nDstXSize;
2580 1749 : nYSize = poKernel->nDstYSize;
2581 1749 : nBitsPerPixel = 32;
2582 1749 : nDefault = 0;
2583 1749 : bDoMemset = false;
2584 : }
2585 : else
2586 : {
2587 0 : CPLError(CE_Failure, CPLE_AppDefined,
2588 : "Internal error in CreateKernelMask(%s).", pszType);
2589 0 : return CE_Failure;
2590 : }
2591 :
2592 : /* -------------------------------------------------------------------- */
2593 : /* Allocate if needed. */
2594 : /* -------------------------------------------------------------------- */
2595 2894 : if (*ppMask == nullptr)
2596 : {
2597 2894 : const GIntBig nBytes =
2598 : nBitsPerPixel == 32
2599 2894 : ? (static_cast<GIntBig>(nXSize) * nYSize + nExtraElts) * 4
2600 959 : : (static_cast<GIntBig>(nXSize) * nYSize + nExtraElts + 31) / 8;
2601 :
2602 2894 : const size_t nByteSize_t = static_cast<size_t>(nBytes);
2603 : #if SIZEOF_VOIDP == 4
2604 : if (static_cast<GIntBig>(nByteSize_t) != nBytes)
2605 : {
2606 : CPLError(CE_Failure, CPLE_OutOfMemory,
2607 : "Cannot allocate " CPL_FRMT_GIB " bytes", nBytes);
2608 : return CE_Failure;
2609 : }
2610 : #endif
2611 :
2612 2894 : *ppMask = VSI_MALLOC_VERBOSE(nByteSize_t);
2613 :
2614 2894 : if (*ppMask == nullptr)
2615 : {
2616 0 : return CE_Failure;
2617 : }
2618 :
2619 2894 : if (bDoMemset)
2620 959 : memset(*ppMask, nDefault, nByteSize_t);
2621 : }
2622 :
2623 2894 : return CE_None;
2624 : }
2625 :
2626 : /************************************************************************/
2627 : /* ComputeSourceWindowStartingFromSource() */
2628 : /************************************************************************/
2629 :
2630 : constexpr int DEFAULT_STEP_COUNT = 21;
2631 :
2632 1261 : void GDALWarpOperation::ComputeSourceWindowStartingFromSource(
2633 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize,
2634 : double *padfSrcMinX, double *padfSrcMinY, double *padfSrcMaxX,
2635 : double *padfSrcMaxY)
2636 : {
2637 1261 : const int nSrcRasterXSize = GDALGetRasterXSize(psOptions->hSrcDS);
2638 1261 : const int nSrcRasterYSize = GDALGetRasterYSize(psOptions->hSrcDS);
2639 1261 : if (nSrcRasterXSize == 0 || nSrcRasterYSize == 0)
2640 0 : return;
2641 :
2642 1261 : GDALWarpPrivateData *privateData = GetWarpPrivateData(this);
2643 1261 : if (privateData->nStepCount == 0)
2644 : {
2645 312 : int nStepCount = DEFAULT_STEP_COUNT;
2646 312 : std::vector<double> adfDstZ{};
2647 :
2648 : const char *pszSampleSteps =
2649 312 : CSLFetchNameValue(psOptions->papszWarpOptions, "SAMPLE_STEPS");
2650 312 : constexpr int knIntMax = std::numeric_limits<int>::max();
2651 312 : if (pszSampleSteps && !EQUAL(pszSampleSteps, "ALL"))
2652 : {
2653 0 : nStepCount = atoi(
2654 0 : CSLFetchNameValue(psOptions->papszWarpOptions, "SAMPLE_STEPS"));
2655 0 : nStepCount = std::max(2, nStepCount);
2656 : }
2657 :
2658 312 : const double dfStepSize = 1.0 / (nStepCount - 1);
2659 312 : if (nStepCount > knIntMax - 2 ||
2660 312 : (nStepCount + 2) > knIntMax / (nStepCount + 2))
2661 : {
2662 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps : %d",
2663 : nStepCount);
2664 0 : return;
2665 : }
2666 312 : const int nSampleMax = (nStepCount + 2) * (nStepCount + 2);
2667 :
2668 : try
2669 : {
2670 312 : privateData->abSuccess.resize(nSampleMax);
2671 312 : privateData->adfDstX.resize(nSampleMax);
2672 312 : privateData->adfDstY.resize(nSampleMax);
2673 312 : adfDstZ.resize(nSampleMax);
2674 : }
2675 0 : catch (const std::exception &)
2676 : {
2677 0 : return;
2678 : }
2679 :
2680 : /* --------------------------------------------------------------------
2681 : */
2682 : /* Setup sample points on a grid pattern throughout the source */
2683 : /* raster. */
2684 : /* --------------------------------------------------------------------
2685 : */
2686 312 : int iPoint = 0;
2687 7488 : for (int iY = 0; iY < nStepCount + 2; iY++)
2688 : {
2689 14040 : const double dfRatioY = (iY == 0) ? 0.5 / nSrcRasterYSize
2690 6864 : : (iY <= nStepCount)
2691 6864 : ? (iY - 1) * dfStepSize
2692 312 : : 1 - 0.5 / nSrcRasterYSize;
2693 172224 : for (int iX = 0; iX < nStepCount + 2; iX++)
2694 : {
2695 322920 : const double dfRatioX = (iX == 0) ? 0.5 / nSrcRasterXSize
2696 157872 : : (iX <= nStepCount)
2697 157872 : ? (iX - 1) * dfStepSize
2698 7176 : : 1 - 0.5 / nSrcRasterXSize;
2699 165048 : privateData->adfDstX[iPoint] = dfRatioX * nSrcRasterXSize;
2700 165048 : privateData->adfDstY[iPoint] = dfRatioY * nSrcRasterYSize;
2701 165048 : iPoint++;
2702 : }
2703 : }
2704 312 : CPLAssert(iPoint == nSampleMax);
2705 :
2706 : /* --------------------------------------------------------------------
2707 : */
2708 : /* Transform them to the output pixel coordinate space */
2709 : /* --------------------------------------------------------------------
2710 : */
2711 312 : psOptions->pfnTransformer(psOptions->pTransformerArg, FALSE, nSampleMax,
2712 : privateData->adfDstX.data(),
2713 : privateData->adfDstY.data(), adfDstZ.data(),
2714 : privateData->abSuccess.data());
2715 312 : privateData->nStepCount = nStepCount;
2716 : }
2717 :
2718 : /* -------------------------------------------------------------------- */
2719 : /* Collect the bounds, ignoring any failed points. */
2720 : /* -------------------------------------------------------------------- */
2721 1261 : const int nStepCount = privateData->nStepCount;
2722 1261 : const double dfStepSize = 1.0 / (nStepCount - 1);
2723 1261 : int iPoint = 0;
2724 : #ifdef DEBUG
2725 1261 : const size_t nSampleMax =
2726 1261 : static_cast<size_t>(nStepCount + 2) * (nStepCount + 2);
2727 1261 : CPL_IGNORE_RET_VAL(nSampleMax);
2728 1261 : CPLAssert(privateData->adfDstX.size() == nSampleMax);
2729 1261 : CPLAssert(privateData->adfDstY.size() == nSampleMax);
2730 1261 : CPLAssert(privateData->abSuccess.size() == nSampleMax);
2731 : #endif
2732 30264 : for (int iY = 0; iY < nStepCount + 2; iY++)
2733 : {
2734 56745 : const double dfRatioY = (iY == 0) ? 0.5 / nSrcRasterYSize
2735 : : (iY <= nStepCount)
2736 27742 : ? (iY - 1) * dfStepSize
2737 1261 : : 1 - 0.5 / nSrcRasterYSize;
2738 696072 : for (int iX = 0; iX < nStepCount + 2; iX++)
2739 : {
2740 667069 : if (privateData->abSuccess[iPoint] &&
2741 631551 : privateData->adfDstX[iPoint] >= nDstXOff &&
2742 459287 : privateData->adfDstX[iPoint] <= nDstXOff + nDstXSize &&
2743 1529120 : privateData->adfDstY[iPoint] >= nDstYOff &&
2744 230495 : privateData->adfDstY[iPoint] <= nDstYOff + nDstYSize)
2745 : {
2746 343551 : const double dfRatioX = (iX == 0) ? 0.5 / nSrcRasterXSize
2747 : : (iX <= nStepCount)
2748 168003 : ? (iX - 1) * dfStepSize
2749 7543 : : 1 - 0.5 / nSrcRasterXSize;
2750 175548 : double dfSrcX = dfRatioX * nSrcRasterXSize;
2751 175548 : double dfSrcY = dfRatioY * nSrcRasterYSize;
2752 175548 : *padfSrcMinX = std::min(*padfSrcMinX, dfSrcX);
2753 175548 : *padfSrcMinY = std::min(*padfSrcMinY, dfSrcY);
2754 175548 : *padfSrcMaxX = std::max(*padfSrcMaxX, dfSrcX);
2755 175548 : *padfSrcMaxY = std::max(*padfSrcMaxY, dfSrcY);
2756 : }
2757 667069 : iPoint++;
2758 : }
2759 : }
2760 : }
2761 :
2762 : /************************************************************************/
2763 : /* ComputeSourceWindowTransformPoints() */
2764 : /************************************************************************/
2765 :
2766 6704 : bool GDALWarpOperation::ComputeSourceWindowTransformPoints(
2767 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize, bool bUseGrid,
2768 : bool bAll, int nStepCount, bool bTryWithCheckWithInvertProj,
2769 : double &dfMinXOut, double &dfMinYOut, double &dfMaxXOut, double &dfMaxYOut,
2770 : int &nSamplePoints, int &nFailedCount)
2771 : {
2772 6704 : nSamplePoints = 0;
2773 6704 : nFailedCount = 0;
2774 :
2775 6704 : const double dfStepSize = bAll ? 0 : 1.0 / (nStepCount - 1);
2776 6704 : constexpr int knIntMax = std::numeric_limits<int>::max();
2777 6704 : int nSampleMax = 0;
2778 6704 : if (bUseGrid)
2779 : {
2780 1323 : if (bAll)
2781 : {
2782 0 : if (nDstXSize > knIntMax - 1 ||
2783 0 : nDstYSize > knIntMax / (nDstXSize + 1) - 1)
2784 : {
2785 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps");
2786 0 : return false;
2787 : }
2788 0 : nSampleMax = (nDstXSize + 1) * (nDstYSize + 1);
2789 : }
2790 : else
2791 : {
2792 1323 : if (nStepCount > knIntMax - 2 ||
2793 1323 : (nStepCount + 2) > knIntMax / (nStepCount + 2))
2794 : {
2795 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps : %d",
2796 : nStepCount);
2797 0 : return false;
2798 : }
2799 1323 : nSampleMax = (nStepCount + 2) * (nStepCount + 2);
2800 : }
2801 : }
2802 : else
2803 : {
2804 5381 : if (bAll)
2805 : {
2806 160 : if (nDstXSize > knIntMax / 2 - nDstYSize)
2807 : {
2808 : // Extremely unlikely !
2809 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps");
2810 0 : return false;
2811 : }
2812 160 : nSampleMax = 2 * (nDstXSize + nDstYSize);
2813 : }
2814 : else
2815 : {
2816 5221 : if (nStepCount > knIntMax / 4)
2817 : {
2818 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps : %d * 4",
2819 : nStepCount);
2820 0 : return false;
2821 : }
2822 5221 : nSampleMax = nStepCount * 4;
2823 : }
2824 : }
2825 :
2826 : int *pabSuccess =
2827 6704 : static_cast<int *>(VSI_MALLOC2_VERBOSE(sizeof(int), nSampleMax));
2828 : double *padfX = static_cast<double *>(
2829 6704 : VSI_MALLOC2_VERBOSE(sizeof(double) * 3, nSampleMax));
2830 6704 : if (pabSuccess == nullptr || padfX == nullptr)
2831 : {
2832 0 : CPLFree(padfX);
2833 0 : CPLFree(pabSuccess);
2834 0 : return false;
2835 : }
2836 6704 : double *padfY = padfX + nSampleMax;
2837 6704 : double *padfZ = padfX + static_cast<size_t>(nSampleMax) * 2;
2838 :
2839 : /* -------------------------------------------------------------------- */
2840 : /* Setup sample points on a grid pattern throughout the area. */
2841 : /* -------------------------------------------------------------------- */
2842 6704 : if (bUseGrid)
2843 : {
2844 1323 : if (bAll)
2845 : {
2846 0 : for (int iY = 0; iY <= nDstYSize; ++iY)
2847 : {
2848 0 : for (int iX = 0; iX <= nDstXSize; ++iX)
2849 : {
2850 0 : padfX[nSamplePoints] = nDstXOff + iX;
2851 0 : padfY[nSamplePoints] = nDstYOff + iY;
2852 0 : padfZ[nSamplePoints++] = 0.0;
2853 : }
2854 : }
2855 : }
2856 : else
2857 : {
2858 31752 : for (int iY = 0; iY < nStepCount + 2; iY++)
2859 : {
2860 59535 : const double dfRatioY = (iY == 0) ? 0.5 / nDstXSize
2861 : : (iY <= nStepCount)
2862 29106 : ? (iY - 1) * dfStepSize
2863 1323 : : 1 - 0.5 / nDstXSize;
2864 730296 : for (int iX = 0; iX < nStepCount + 2; iX++)
2865 : {
2866 1369300 : const double dfRatioX = (iX == 0) ? 0.5 / nDstXSize
2867 : : (iX <= nStepCount)
2868 669438 : ? (iX - 1) * dfStepSize
2869 30429 : : 1 - 0.5 / nDstXSize;
2870 699867 : padfX[nSamplePoints] = dfRatioX * nDstXSize + nDstXOff;
2871 699867 : padfY[nSamplePoints] = dfRatioY * nDstYSize + nDstYOff;
2872 699867 : padfZ[nSamplePoints++] = 0.0;
2873 : }
2874 : }
2875 : }
2876 : }
2877 : /* -------------------------------------------------------------------- */
2878 : /* Setup sample points all around the edge of the output raster. */
2879 : /* -------------------------------------------------------------------- */
2880 : else
2881 : {
2882 5381 : if (bAll)
2883 : {
2884 68927 : for (int iX = 0; iX <= nDstXSize; ++iX)
2885 : {
2886 : // Along top
2887 68767 : padfX[nSamplePoints] = nDstXOff + iX;
2888 68767 : padfY[nSamplePoints] = nDstYOff;
2889 68767 : padfZ[nSamplePoints++] = 0.0;
2890 :
2891 : // Along bottom
2892 68767 : padfX[nSamplePoints] = nDstXOff + iX;
2893 68767 : padfY[nSamplePoints] = nDstYOff + nDstYSize;
2894 68767 : padfZ[nSamplePoints++] = 0.0;
2895 : }
2896 :
2897 44154 : for (int iY = 1; iY < nDstYSize; ++iY)
2898 : {
2899 : // Along left
2900 43994 : padfX[nSamplePoints] = nDstXOff;
2901 43994 : padfY[nSamplePoints] = nDstYOff + iY;
2902 43994 : padfZ[nSamplePoints++] = 0.0;
2903 :
2904 : // Along right
2905 43994 : padfX[nSamplePoints] = nDstXOff + nDstXSize;
2906 43994 : padfY[nSamplePoints] = nDstYOff + iY;
2907 43994 : padfZ[nSamplePoints++] = 0.0;
2908 : }
2909 : }
2910 : else
2911 : {
2912 114862 : for (double dfRatio = 0.0; dfRatio <= 1.0 + dfStepSize * 0.5;
2913 109641 : dfRatio += dfStepSize)
2914 : {
2915 : // Along top
2916 109641 : padfX[nSamplePoints] = dfRatio * nDstXSize + nDstXOff;
2917 109641 : padfY[nSamplePoints] = nDstYOff;
2918 109641 : padfZ[nSamplePoints++] = 0.0;
2919 :
2920 : // Along bottom
2921 109641 : padfX[nSamplePoints] = dfRatio * nDstXSize + nDstXOff;
2922 109641 : padfY[nSamplePoints] = nDstYOff + nDstYSize;
2923 109641 : padfZ[nSamplePoints++] = 0.0;
2924 :
2925 : // Along left
2926 109641 : padfX[nSamplePoints] = nDstXOff;
2927 109641 : padfY[nSamplePoints] = dfRatio * nDstYSize + nDstYOff;
2928 109641 : padfZ[nSamplePoints++] = 0.0;
2929 :
2930 : // Along right
2931 109641 : padfX[nSamplePoints] = nDstXSize + nDstXOff;
2932 109641 : padfY[nSamplePoints] = dfRatio * nDstYSize + nDstYOff;
2933 109641 : padfZ[nSamplePoints++] = 0.0;
2934 : }
2935 : }
2936 : }
2937 :
2938 6704 : CPLAssert(nSamplePoints == nSampleMax);
2939 :
2940 : /* -------------------------------------------------------------------- */
2941 : /* Transform them to the input pixel coordinate space */
2942 : /* -------------------------------------------------------------------- */
2943 :
2944 138 : const auto RefreshTransformer = [this]()
2945 : {
2946 46 : if (GDALIsTransformer(psOptions->pTransformerArg,
2947 : GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
2948 : {
2949 0 : GDALRefreshGenImgProjTransformer(psOptions->pTransformerArg);
2950 : }
2951 46 : else if (GDALIsTransformer(psOptions->pTransformerArg,
2952 : GDAL_APPROX_TRANSFORMER_CLASS_NAME))
2953 : {
2954 46 : GDALRefreshApproxTransformer(psOptions->pTransformerArg);
2955 : }
2956 6750 : };
2957 :
2958 6704 : if (bTryWithCheckWithInvertProj)
2959 : {
2960 23 : CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ", "YES");
2961 23 : RefreshTransformer();
2962 : }
2963 6704 : psOptions->pfnTransformer(psOptions->pTransformerArg, TRUE, nSamplePoints,
2964 : padfX, padfY, padfZ, pabSuccess);
2965 6704 : if (bTryWithCheckWithInvertProj)
2966 : {
2967 23 : CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ", nullptr);
2968 23 : RefreshTransformer();
2969 : }
2970 :
2971 : /* -------------------------------------------------------------------- */
2972 : /* Collect the bounds, ignoring any failed points. */
2973 : /* -------------------------------------------------------------------- */
2974 1370660 : for (int i = 0; i < nSamplePoints; i++)
2975 : {
2976 1363950 : if (!pabSuccess[i])
2977 : {
2978 112179 : nFailedCount++;
2979 112179 : continue;
2980 : }
2981 :
2982 : // If this happens this is likely the symptom of a bug somewhere.
2983 1251770 : if (std::isnan(padfX[i]) || std::isnan(padfY[i]))
2984 : {
2985 : static bool bNanCoordFound = false;
2986 0 : if (!bNanCoordFound)
2987 : {
2988 0 : CPLDebug("WARP",
2989 : "ComputeSourceWindow(): "
2990 : "NaN coordinate found on point %d.",
2991 : i);
2992 0 : bNanCoordFound = true;
2993 : }
2994 0 : nFailedCount++;
2995 0 : continue;
2996 : }
2997 :
2998 1251770 : dfMinXOut = std::min(dfMinXOut, padfX[i]);
2999 1251770 : dfMinYOut = std::min(dfMinYOut, padfY[i]);
3000 1251770 : dfMaxXOut = std::max(dfMaxXOut, padfX[i]);
3001 1251770 : dfMaxYOut = std::max(dfMaxYOut, padfY[i]);
3002 : }
3003 :
3004 6704 : CPLFree(padfX);
3005 6704 : CPLFree(pabSuccess);
3006 6704 : return true;
3007 : }
3008 :
3009 : /************************************************************************/
3010 : /* ComputeSourceWindow() */
3011 : /************************************************************************/
3012 :
3013 : /** Given a target window starting at pixel (nDstOff, nDstYOff) and of
3014 : * dimension (nDstXSize, nDstYSize), compute the corresponding window in
3015 : * the source raster, and return the source position in (*pnSrcXOff, *pnSrcYOff),
3016 : * the source dimension in (*pnSrcXSize, *pnSrcYSize).
3017 : * If pdfSrcXExtraSize is not null, its pointed value will be filled with the
3018 : * number of extra source pixels in X dimension to acquire to take into account
3019 : * the size of the resampling kernel. Similarly for pdfSrcYExtraSize for the
3020 : * Y dimension.
3021 : * If pdfSrcFillRatio is not null, its pointed value will be filled with the
3022 : * the ratio of the clamped source raster window size over the unclamped source
3023 : * raster window size. When this ratio is too low, this might be an indication
3024 : * that it might be beneficial to split the target window to avoid requesting
3025 : * too many source pixels.
3026 : */
3027 6383 : CPLErr GDALWarpOperation::ComputeSourceWindow(
3028 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize, int *pnSrcXOff,
3029 : int *pnSrcYOff, int *pnSrcXSize, int *pnSrcYSize, double *pdfSrcXExtraSize,
3030 : double *pdfSrcYExtraSize, double *pdfSrcFillRatio)
3031 :
3032 : {
3033 : /* -------------------------------------------------------------------- */
3034 : /* Figure out whether we just want to do the usual "along the */
3035 : /* edge" sampling, or using a grid. The grid usage is */
3036 : /* important in some weird "inside out" cases like WGS84 to */
3037 : /* polar stereographic around the pole. Also figure out the */
3038 : /* sampling rate. */
3039 : /* -------------------------------------------------------------------- */
3040 6383 : int nStepCount = DEFAULT_STEP_COUNT;
3041 6383 : bool bAll = false;
3042 :
3043 : bool bUseGrid =
3044 6383 : CPLFetchBool(psOptions->papszWarpOptions, "SAMPLE_GRID", false);
3045 :
3046 : const char *pszSampleSteps =
3047 6383 : CSLFetchNameValue(psOptions->papszWarpOptions, "SAMPLE_STEPS");
3048 6383 : if (pszSampleSteps)
3049 : {
3050 94 : if (EQUAL(pszSampleSteps, "ALL"))
3051 : {
3052 94 : bAll = true;
3053 : }
3054 : else
3055 : {
3056 0 : nStepCount = atoi(pszSampleSteps);
3057 0 : nStepCount = std::max(2, nStepCount);
3058 : }
3059 : }
3060 6289 : else if (!bUseGrid)
3061 : {
3062 : // Detect if at least one of the 4 corner in destination raster fails
3063 : // to project back to source.
3064 : // Helps for long-lat to orthographic on areas that are partly in
3065 : // space / partly on Earth. Cf https://github.com/OSGeo/gdal/issues/9056
3066 : double adfCornerX[4];
3067 : double adfCornerY[4];
3068 5287 : double adfCornerZ[4] = {0, 0, 0, 0};
3069 5287 : int anCornerSuccess[4] = {FALSE, FALSE, FALSE, FALSE};
3070 5287 : adfCornerX[0] = nDstXOff;
3071 5287 : adfCornerY[0] = nDstYOff;
3072 5287 : adfCornerX[1] = nDstXOff + nDstXSize;
3073 5287 : adfCornerY[1] = nDstYOff;
3074 5287 : adfCornerX[2] = nDstXOff;
3075 5287 : adfCornerY[2] = nDstYOff + nDstYSize;
3076 5287 : adfCornerX[3] = nDstXOff + nDstXSize;
3077 5287 : adfCornerY[3] = nDstYOff + nDstYSize;
3078 5287 : if (!psOptions->pfnTransformer(psOptions->pTransformerArg, TRUE, 4,
3079 : adfCornerX, adfCornerY, adfCornerZ,
3080 5221 : anCornerSuccess) ||
3081 10508 : !anCornerSuccess[0] || !anCornerSuccess[1] || !anCornerSuccess[2] ||
3082 5221 : !anCornerSuccess[3])
3083 : {
3084 66 : bAll = true;
3085 : }
3086 : }
3087 :
3088 6383 : bool bTryWithCheckWithInvertProj = false;
3089 6383 : double dfMinXOut = std::numeric_limits<double>::infinity();
3090 6383 : double dfMinYOut = std::numeric_limits<double>::infinity();
3091 6383 : double dfMaxXOut = -std::numeric_limits<double>::infinity();
3092 6383 : double dfMaxYOut = -std::numeric_limits<double>::infinity();
3093 :
3094 6383 : int nSamplePoints = 0;
3095 6383 : int nFailedCount = 0;
3096 6383 : if (!ComputeSourceWindowTransformPoints(
3097 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, bUseGrid, bAll,
3098 : nStepCount, bTryWithCheckWithInvertProj, dfMinXOut, dfMinYOut,
3099 : dfMaxXOut, dfMaxYOut, nSamplePoints, nFailedCount))
3100 : {
3101 0 : return CE_Failure;
3102 : }
3103 :
3104 : // Use grid sampling as soon as a special point falls into the extent of
3105 : // the target raster.
3106 6383 : if (!bUseGrid && psOptions->hDstDS)
3107 : {
3108 13029 : for (const auto &xy : aDstXYSpecialPoints)
3109 : {
3110 17525 : if (0 <= xy.first &&
3111 1737 : GDALGetRasterXSize(psOptions->hDstDS) >= xy.first &&
3112 10427 : 0 <= xy.second &&
3113 796 : GDALGetRasterYSize(psOptions->hDstDS) >= xy.second)
3114 : {
3115 233 : bUseGrid = true;
3116 233 : bAll = false;
3117 233 : if (!ComputeSourceWindowTransformPoints(
3118 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, bUseGrid,
3119 : bAll, nStepCount, bTryWithCheckWithInvertProj,
3120 : dfMinXOut, dfMinYOut, dfMaxXOut, dfMaxYOut,
3121 : nSamplePoints, nFailedCount))
3122 : {
3123 0 : return CE_Failure;
3124 : }
3125 233 : break;
3126 : }
3127 : }
3128 : }
3129 :
3130 6383 : const int nRasterXSize = GDALGetRasterXSize(psOptions->hSrcDS);
3131 6383 : const int nRasterYSize = GDALGetRasterYSize(psOptions->hSrcDS);
3132 :
3133 : // Try to detect crazy values coming from reprojection that would not
3134 : // have resulted in a PROJ error. Could happen for example with PROJ
3135 : // <= 4.9.2 with inverse UTM/tmerc (Snyder approximation without sanity
3136 : // check) when being far away from the central meridian. But might be worth
3137 : // keeping that even for later versions in case some exotic projection isn't
3138 : // properly sanitized.
3139 6309 : if (nFailedCount == 0 && !bTryWithCheckWithInvertProj &&
3140 6309 : (dfMinXOut < -1e6 || dfMinYOut < -1e6 ||
3141 12696 : dfMaxXOut > nRasterXSize + 1e6 || dfMaxYOut > nRasterYSize + 1e6) &&
3142 23 : !CPLTestBool(CPLGetConfigOption("CHECK_WITH_INVERT_PROJ", "NO")))
3143 : {
3144 23 : CPLDebug("WARP",
3145 : "ComputeSourceWindow(): bogus source dataset window "
3146 : "returned. Trying again with CHECK_WITH_INVERT_PROJ=YES");
3147 23 : bTryWithCheckWithInvertProj = true;
3148 :
3149 : // We should probably perform the coordinate transformation in the
3150 : // warp kernel under CHECK_WITH_INVERT_PROJ too...
3151 23 : if (!ComputeSourceWindowTransformPoints(
3152 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, bUseGrid, bAll,
3153 : nStepCount, bTryWithCheckWithInvertProj, dfMinXOut, dfMinYOut,
3154 : dfMaxXOut, dfMaxYOut, nSamplePoints, nFailedCount))
3155 : {
3156 0 : return CE_Failure;
3157 : }
3158 : }
3159 :
3160 : /* -------------------------------------------------------------------- */
3161 : /* If we got any failures when not using a grid, we should */
3162 : /* really go back and try again with the grid. Sorry for the */
3163 : /* goto. */
3164 : /* -------------------------------------------------------------------- */
3165 6383 : if (!bUseGrid && nFailedCount > 0)
3166 : {
3167 65 : bUseGrid = true;
3168 65 : bAll = false;
3169 65 : if (!ComputeSourceWindowTransformPoints(
3170 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, bUseGrid, bAll,
3171 : nStepCount, bTryWithCheckWithInvertProj, dfMinXOut, dfMinYOut,
3172 : dfMaxXOut, dfMaxYOut, nSamplePoints, nFailedCount))
3173 : {
3174 0 : return CE_Failure;
3175 : }
3176 : }
3177 :
3178 : /* -------------------------------------------------------------------- */
3179 : /* If we get hardly any points (or none) transforming, we give */
3180 : /* up. */
3181 : /* -------------------------------------------------------------------- */
3182 6383 : if (nFailedCount > nSamplePoints - 5)
3183 : {
3184 : const bool bErrorOutIfEmptySourceWindow =
3185 39 : CPLFetchBool(psOptions->papszWarpOptions,
3186 : "ERROR_OUT_IF_EMPTY_SOURCE_WINDOW", true);
3187 39 : if (bErrorOutIfEmptySourceWindow)
3188 : {
3189 3 : CPLError(CE_Failure, CPLE_AppDefined,
3190 : "Too many points (%d out of %d) failed to transform, "
3191 : "unable to compute output bounds.",
3192 : nFailedCount, nSamplePoints);
3193 : }
3194 : else
3195 : {
3196 36 : CPLDebug("WARP", "Cannot determine source window for %d,%d,%d,%d",
3197 : nDstXOff, nDstYOff, nDstXSize, nDstYSize);
3198 : }
3199 39 : return CE_Failure;
3200 : }
3201 :
3202 6344 : if (nFailedCount > 0)
3203 39 : CPLDebug("GDAL",
3204 : "GDALWarpOperation::ComputeSourceWindow() %d out of %d "
3205 : "points failed to transform.",
3206 : nFailedCount, nSamplePoints);
3207 :
3208 : /* -------------------------------------------------------------------- */
3209 : /* In some cases (see https://github.com/OSGeo/gdal/issues/862) */
3210 : /* the reverse transform does not work at some points, so try by */
3211 : /* transforming from source raster space to target raster space and */
3212 : /* see which source coordinates end up being in the AOI in the target */
3213 : /* raster space. */
3214 : /* -------------------------------------------------------------------- */
3215 6344 : if (bUseGrid)
3216 : {
3217 1261 : ComputeSourceWindowStartingFromSource(nDstXOff, nDstYOff, nDstXSize,
3218 : nDstYSize, &dfMinXOut, &dfMinYOut,
3219 : &dfMaxXOut, &dfMaxYOut);
3220 : }
3221 :
3222 : /* -------------------------------------------------------------------- */
3223 : /* Early exit to avoid crazy values to cause a huge nResWinSize that */
3224 : /* would result in a result window wrongly covering the whole raster. */
3225 : /* -------------------------------------------------------------------- */
3226 6344 : if (dfMinXOut > nRasterXSize || dfMaxXOut < 0 || dfMinYOut > nRasterYSize ||
3227 5274 : dfMaxYOut < 0)
3228 : {
3229 1531 : *pnSrcXOff = 0;
3230 1531 : *pnSrcYOff = 0;
3231 1531 : *pnSrcXSize = 0;
3232 1531 : *pnSrcYSize = 0;
3233 1531 : if (pdfSrcXExtraSize)
3234 1531 : *pdfSrcXExtraSize = 0.0;
3235 1531 : if (pdfSrcYExtraSize)
3236 1531 : *pdfSrcYExtraSize = 0.0;
3237 1531 : if (pdfSrcFillRatio)
3238 949 : *pdfSrcFillRatio = 0.0;
3239 1531 : return CE_None;
3240 : }
3241 :
3242 : // For scenarios where warping is used as a "decoration", try to clamp
3243 : // source pixel coordinates to integer when very close.
3244 19252 : const auto roundIfCloseEnough = [](double dfVal)
3245 : {
3246 19252 : const double dfRounded = std::round(dfVal);
3247 19252 : if (std::fabs(dfRounded - dfVal) < 1e-6)
3248 15213 : return dfRounded;
3249 4039 : return dfVal;
3250 : };
3251 :
3252 4813 : dfMinXOut = roundIfCloseEnough(dfMinXOut);
3253 4813 : dfMinYOut = roundIfCloseEnough(dfMinYOut);
3254 4813 : dfMaxXOut = roundIfCloseEnough(dfMaxXOut);
3255 4813 : dfMaxYOut = roundIfCloseEnough(dfMaxYOut);
3256 :
3257 4813 : if (m_bIsTranslationOnPixelBoundaries)
3258 : {
3259 376 : CPLAssert(dfMinXOut == std::round(dfMinXOut));
3260 376 : CPLAssert(dfMinYOut == std::round(dfMinYOut));
3261 376 : CPLAssert(dfMaxXOut == std::round(dfMaxXOut));
3262 376 : CPLAssert(dfMaxYOut == std::round(dfMaxYOut));
3263 376 : CPLAssert(std::round(dfMaxXOut - dfMinXOut) == nDstXSize);
3264 376 : CPLAssert(std::round(dfMaxYOut - dfMinYOut) == nDstYSize);
3265 : }
3266 :
3267 : /* -------------------------------------------------------------------- */
3268 : /* How much of a window around our source pixel might we need */
3269 : /* to collect data from based on the resampling kernel? Even */
3270 : /* if the requested central pixel falls off the source image, */
3271 : /* we may need to collect data if some portion of the */
3272 : /* resampling kernel could be on-image. */
3273 : /* -------------------------------------------------------------------- */
3274 4813 : const int nResWinSize = m_bIsTranslationOnPixelBoundaries
3275 4813 : ? 0
3276 4437 : : GWKGetFilterRadius(psOptions->eResampleAlg);
3277 :
3278 : // Take scaling into account.
3279 : // Avoid ridiculous small scaling factors to avoid potential further integer
3280 : // overflows
3281 9626 : const double dfXScale = std::max(1e-3, static_cast<double>(nDstXSize) /
3282 4813 : (dfMaxXOut - dfMinXOut));
3283 9626 : const double dfYScale = std::max(1e-3, static_cast<double>(nDstYSize) /
3284 4813 : (dfMaxYOut - dfMinYOut));
3285 4813 : int nXRadius = dfXScale < 0.95
3286 4813 : ? static_cast<int>(ceil(nResWinSize / dfXScale))
3287 : : nResWinSize;
3288 4813 : int nYRadius = dfYScale < 0.95
3289 4813 : ? static_cast<int>(ceil(nResWinSize / dfYScale))
3290 : : nResWinSize;
3291 :
3292 : /* -------------------------------------------------------------------- */
3293 : /* Allow addition of extra sample pixels to source window to */
3294 : /* avoid missing pixels due to sampling error. In fact, */
3295 : /* fallback to adding a bit to the window if any points failed */
3296 : /* to transform. */
3297 : /* -------------------------------------------------------------------- */
3298 4813 : if (const char *pszSourceExtra =
3299 4813 : CSLFetchNameValue(psOptions->papszWarpOptions, "SOURCE_EXTRA"))
3300 : {
3301 81 : int nSrcExtra = cpl::strict_parse<int>(pszSourceExtra).value_or(-1);
3302 :
3303 81 : if (nSrcExtra < 0)
3304 : {
3305 : // no point raising CE_Failure because it will get converted into
3306 : // a warning at an outer scope
3307 1 : CPLError(CE_Warning, CPLE_IllegalArg,
3308 : "SOURCE_EXTRA must be a positive integer or zero.");
3309 1 : nSrcExtra = 0;
3310 : }
3311 :
3312 81 : nXRadius += nSrcExtra;
3313 81 : nYRadius += nSrcExtra;
3314 : }
3315 4732 : else if (nFailedCount > 0)
3316 : {
3317 33 : nXRadius += 10;
3318 33 : nYRadius += 10;
3319 : }
3320 :
3321 : /* -------------------------------------------------------------------- */
3322 : /* return bounds. */
3323 : /* -------------------------------------------------------------------- */
3324 : #if DEBUG_VERBOSE
3325 : CPLDebug("WARP",
3326 : "dst=(%d,%d,%d,%d) raw "
3327 : "src=(minx=%.17g,miny=%.17g,maxx=%.17g,maxy=%.17g)",
3328 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, dfMinXOut, dfMinYOut,
3329 : dfMaxXOut, dfMaxYOut);
3330 : #endif
3331 4813 : const int nMinXOutClamped = static_cast<int>(std::max(0.0, dfMinXOut));
3332 4813 : const int nMinYOutClamped = static_cast<int>(std::max(0.0, dfMinYOut));
3333 : const int nMaxXOutClamped = static_cast<int>(
3334 4813 : std::min(ceil(dfMaxXOut), static_cast<double>(nRasterXSize)));
3335 : const int nMaxYOutClamped = static_cast<int>(
3336 4813 : std::min(ceil(dfMaxYOut), static_cast<double>(nRasterYSize)));
3337 :
3338 : const double dfSrcXSizeRaw = std::max(
3339 14439 : 0.0, std::min(static_cast<double>(nRasterXSize - nMinXOutClamped),
3340 4813 : dfMaxXOut - dfMinXOut));
3341 : const double dfSrcYSizeRaw = std::max(
3342 14439 : 0.0, std::min(static_cast<double>(nRasterYSize - nMinYOutClamped),
3343 4813 : dfMaxYOut - dfMinYOut));
3344 :
3345 : // If we cover more than 90% of the width, then use it fully (helps for
3346 : // anti-meridian discontinuities)
3347 4813 : if (nMaxXOutClamped - nMinXOutClamped > 0.9 * nRasterXSize)
3348 : {
3349 1595 : *pnSrcXOff = 0;
3350 1595 : *pnSrcXSize = nRasterXSize;
3351 : }
3352 : else
3353 : {
3354 3218 : *pnSrcXOff =
3355 3218 : std::max(0, std::min(nMinXOutClamped - nXRadius, nRasterXSize));
3356 3218 : *pnSrcXSize =
3357 9654 : std::max(0, std::min(nRasterXSize - *pnSrcXOff,
3358 3218 : nMaxXOutClamped - *pnSrcXOff + nXRadius));
3359 : }
3360 :
3361 4813 : if (nMaxYOutClamped - nMinYOutClamped > 0.9 * nRasterYSize)
3362 : {
3363 1493 : *pnSrcYOff = 0;
3364 1493 : *pnSrcYSize = nRasterYSize;
3365 : }
3366 : else
3367 : {
3368 3320 : *pnSrcYOff =
3369 3320 : std::max(0, std::min(nMinYOutClamped - nYRadius, nRasterYSize));
3370 3320 : *pnSrcYSize =
3371 9960 : std::max(0, std::min(nRasterYSize - *pnSrcYOff,
3372 3320 : nMaxYOutClamped - *pnSrcYOff + nYRadius));
3373 : }
3374 :
3375 4813 : if (pdfSrcXExtraSize)
3376 4813 : *pdfSrcXExtraSize = *pnSrcXSize - dfSrcXSizeRaw;
3377 4813 : if (pdfSrcYExtraSize)
3378 4813 : *pdfSrcYExtraSize = *pnSrcYSize - dfSrcYSizeRaw;
3379 :
3380 : // Computed the ratio of the clamped source raster window size over
3381 : // the unclamped source raster window size.
3382 4813 : if (pdfSrcFillRatio)
3383 3653 : *pdfSrcFillRatio =
3384 7306 : static_cast<double>(*pnSrcXSize) * (*pnSrcYSize) /
3385 3653 : std::max(1.0, (dfMaxXOut - dfMinXOut + 2 * nXRadius) *
3386 3653 : (dfMaxYOut - dfMinYOut + 2 * nYRadius));
3387 :
3388 4813 : return CE_None;
3389 : }
3390 :
3391 : /************************************************************************/
3392 : /* ReportTiming() */
3393 : /************************************************************************/
3394 :
3395 12737 : void GDALWarpOperation::ReportTiming(const char *pszMessage)
3396 :
3397 : {
3398 12737 : if (!bReportTimings)
3399 12737 : return;
3400 :
3401 0 : const unsigned long nNewTime = VSITime(nullptr);
3402 :
3403 0 : if (pszMessage != nullptr)
3404 : {
3405 0 : CPLDebug("WARP_TIMING", "%s: %lds", pszMessage,
3406 0 : static_cast<long>(nNewTime - nLastTimeReported));
3407 : }
3408 :
3409 0 : nLastTimeReported = nNewTime;
3410 : }
|
