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