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