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