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 1748 : GDALWarpOperation::~GDALWarpOperation()
164 :
165 : {
166 : {
167 3496 : std::lock_guard<std::mutex> oLock(gMutex);
168 1748 : auto oItem = gMapPrivate.find(this);
169 1748 : if (oItem != gMapPrivate.end())
170 : {
171 260 : gMapPrivate.erase(oItem);
172 : }
173 : }
174 :
175 1748 : WipeOptions();
176 :
177 1748 : if (hIOMutex != nullptr)
178 : {
179 6 : CPLDestroyMutex(hIOMutex);
180 6 : CPLDestroyMutex(hWarpMutex);
181 : }
182 :
183 1748 : WipeChunkList();
184 1748 : if (psThreadData)
185 1744 : GWKThreadsEnd(psThreadData);
186 1748 : }
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 1752 : void GDALWarpOperation::WipeOptions()
204 :
205 : {
206 1752 : if (psOptions != nullptr)
207 : {
208 1748 : GDALDestroyWarpOptions(psOptions);
209 1748 : psOptions = nullptr;
210 : }
211 1752 : }
212 :
213 : /************************************************************************/
214 : /* ValidateOptions() */
215 : /************************************************************************/
216 :
217 1748 : int GDALWarpOperation::ValidateOptions()
218 :
219 : {
220 1748 : 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 1748 : 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 1748 : if (psOptions->eResampleAlg != GRA_NearestNeighbour &&
238 982 : psOptions->eResampleAlg != GRA_Bilinear &&
239 582 : 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 1748 : if (static_cast<int>(psOptions->eWorkingDataType) < 1 ||
259 1748 : 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 2001 : 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 1748 : 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 1748 : 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 1748 : 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 1748 : 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 4313 : for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
315 : {
316 5130 : if (psOptions->panSrcBands[iBand] < 1 ||
317 2565 : psOptions->panSrcBands[iBand] >
318 2565 : 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 5117 : if (psOptions->hDstDS != nullptr &&
326 2552 : (psOptions->panDstBands[iBand] < 1 ||
327 2552 : psOptions->panDstBands[iBand] >
328 2552 : 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 5117 : if (psOptions->hDstDS != nullptr &&
337 2552 : GDALGetRasterAccess(GDALGetRasterBand(
338 2552 : 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 1748 : 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 1748 : 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 1748 : 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 3496 : CPLStringList aosWO(CSLDuplicate(psOptions->papszWarpOptions));
374 : // A few internal/undocumented options
375 1748 : aosWO.SetNameValue("EXTRA_ELTS", nullptr);
376 1748 : aosWO.SetNameValue("USE_GENERAL_CASE", nullptr);
377 1748 : aosWO.SetNameValue("ERROR_THRESHOLD", nullptr);
378 1748 : aosWO.SetNameValue("ERROR_OUT_IF_EMPTY_SOURCE_WINDOW", nullptr);
379 1748 : aosWO.SetNameValue("MULT_FACTOR_VERTICAL_SHIFT_PIPELINE", nullptr);
380 1748 : aosWO.SetNameValue("SRC_FILL_RATIO_HEURISTICS", nullptr);
381 1748 : GDALValidateOptions(GDALWarpGetOptionList(), aosWO.List(), "option",
382 : "warp options");
383 : }
384 :
385 : const char *pszSampleSteps =
386 1748 : CSLFetchNameValue(psOptions->papszWarpOptions, "SAMPLE_STEPS");
387 1748 : 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 1748 : 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 1748 : 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 1748 : 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 1748 : 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 1748 : GDALGetRasterBand(psOptions->hSrcDS, psOptions->panSrcBands[0]);
442 1750 : 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 3496 : const bool bErrorOutIfEmptySourceWindow = CPLFetchBool(
454 1748 : psOptions->papszWarpOptions, "ERROR_OUT_IF_EMPTY_SOURCE_WINDOW", true);
455 2118 : 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 1748 : 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 1748 : static void SetTieStrategy(GDALWarpOptions *psOptions, CPLErr *peErr)
503 : {
504 1748 : if (const char *pszTieStrategy =
505 1748 : 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 1748 : }
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 1748 : GDALWarpOperation::Initialize(const GDALWarpOptions *psNewOptions,
559 : GDALTransformerFunc pfnTransformer,
560 : GDALTransformerArgUniquePtr psOwnedTransformerArg)
561 :
562 : {
563 : /* -------------------------------------------------------------------- */
564 : /* Copy the passed in options. */
565 : /* -------------------------------------------------------------------- */
566 1748 : if (psOptions != nullptr)
567 0 : WipeOptions();
568 :
569 1748 : CPLErr eErr = CE_None;
570 :
571 1748 : psOptions = GDALCloneWarpOptions(psNewOptions);
572 :
573 1748 : if (psOptions->pfnTransformer)
574 : {
575 785 : CPLAssert(pfnTransformer == nullptr);
576 785 : CPLAssert(psOwnedTransformerArg.get() == nullptr);
577 : }
578 : else
579 : {
580 963 : m_psOwnedTransformerArg = std::move(psOwnedTransformerArg);
581 963 : psOptions->pfnTransformer = pfnTransformer;
582 963 : psOptions->pTransformerArg = m_psOwnedTransformerArg.get();
583 : }
584 :
585 3496 : psOptions->papszWarpOptions =
586 1748 : 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 1748 : 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 1748 : GDALWarpResolveWorkingDataType(psOptions);
602 1748 : 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 1748 : if (psOptions->dfWarpMemoryLimit == 0.0)
612 : {
613 1488 : psOptions->dfWarpMemoryLimit = 64.0 * 1024 * 1024;
614 : }
615 :
616 : /* -------------------------------------------------------------------- */
617 : /* Are we doing timings? */
618 : /* -------------------------------------------------------------------- */
619 1748 : bReportTimings =
620 1748 : CPLFetchBool(psOptions->papszWarpOptions, "REPORT_TIMINGS", false);
621 :
622 : /* -------------------------------------------------------------------- */
623 : /* Support creating cutline from text warpoption. */
624 : /* -------------------------------------------------------------------- */
625 : const char *pszCutlineWKT =
626 1748 : CSLFetchNameValue(psOptions->papszWarpOptions, "CUTLINE");
627 :
628 1748 : 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 1748 : CSLFetchNameValue(psOptions->papszWarpOptions, "CUTLINE_BLEND_DIST");
642 1748 : if (pszBD)
643 0 : psOptions->dfCutlineBlendDist = CPLAtof(pszBD);
644 :
645 : /* -------------------------------------------------------------------- */
646 : /* Set SRC_ALPHA_MAX if not provided. */
647 : /* -------------------------------------------------------------------- */
648 1748 : if (psOptions->hSrcDS != nullptr && psOptions->nSrcAlphaBand > 0 &&
649 3597 : 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 1735 : if (psOptions->hDstDS != nullptr && psOptions->nDstAlphaBand > 0 &&
662 3870 : 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 1748 : if (!ValidateOptions())
675 0 : eErr = CE_Failure;
676 :
677 1748 : if (eErr != CE_None)
678 : {
679 4 : WipeOptions();
680 : }
681 : else
682 : {
683 3488 : psThreadData = GWKThreadsCreate(psOptions->papszWarpOptions,
684 1744 : psOptions->pfnTransformer,
685 1744 : psOptions->pTransformerArg);
686 1744 : 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 5232 : for (double dfY : {-89.9999, 89.9999})
698 : {
699 3488 : double dfX = 0;
700 3488 : if ((GDALIsTransformer(psOptions->pTransformerArg,
701 2828 : GDAL_APPROX_TRANSFORMER_CLASS_NAME) &&
702 2828 : GDALTransformLonLatToDestApproxTransformer(
703 6976 : psOptions->pTransformerArg, &dfX, &dfY)) ||
704 1984 : (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 1572 : std::pair<double, double>(dfX, dfY));
711 : }
712 : }
713 :
714 1744 : m_bIsTranslationOnPixelBoundaries =
715 3488 : GDALTransformIsTranslationOnPixelBoundaries(
716 1975 : psOptions->pfnTransformer, psOptions->pTransformerArg) &&
717 231 : CPLTestBool(
718 : CPLGetConfigOption("GDAL_WARP_USE_TRANSLATION_OPTIM", "YES"));
719 1744 : if (m_bIsTranslationOnPixelBoundaries)
720 : {
721 225 : CPLDebug("WARP",
722 : "Using translation-on-pixel-boundaries optimization");
723 : }
724 : }
725 :
726 3479 : if (eErr == CE_None && psOptions->hDstDS &&
727 1731 : 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 1748 : 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 2655 : 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 2655 : const int nWordSize = GDALGetDataTypeSizeBytes(psOptions->eWorkingDataType);
768 :
769 2655 : void *pDstBuffer = VSI_MALLOC3_VERBOSE(
770 : cpl::fits_on<int>(nWordSize * psOptions->nBandCount), nDstXSize,
771 : nDstYSize);
772 2655 : if (pDstBuffer)
773 : {
774 2655 : auto eErr = InitializeDestinationBuffer(pDstBuffer, nDstXSize,
775 : nDstYSize, pbInitialized);
776 2655 : if (eErr != CE_None)
777 : {
778 2 : CPLFree(pDstBuffer);
779 2 : return nullptr;
780 : }
781 : }
782 2653 : 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 2734 : CPLErr GDALWarpOperation::InitializeDestinationBuffer(void *pDstBuffer,
805 : int nDstXSize,
806 : int nDstYSize,
807 : int *pbInitialized) const
808 : {
809 2734 : const int nWordSize = GDALGetDataTypeSizeBytes(psOptions->eWorkingDataType);
810 :
811 2734 : const GPtrDiff_t nBandSize =
812 2734 : static_cast<GPtrDiff_t>(nWordSize) * nDstXSize * nDstYSize;
813 :
814 : /* -------------------------------------------------------------------- */
815 : /* Initialize if requested in the options */
816 : /* -------------------------------------------------------------------- */
817 : const char *pszInitDest =
818 2734 : CSLFetchNameValue(psOptions->papszWarpOptions, "INIT_DEST");
819 :
820 2734 : 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 2344 : if (pbInitialized != nullptr)
830 : {
831 1918 : *pbInitialized = TRUE;
832 : }
833 :
834 : CPLStringList aosInitValues(
835 4688 : CSLTokenizeStringComplex(pszInitDest, ",", FALSE, FALSE));
836 2344 : const int nInitCount = aosInitValues.Count();
837 :
838 5386 : for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
839 : {
840 3044 : double adfInitRealImag[2] = {0.0, 0.0};
841 : const char *pszBandInit =
842 3044 : aosInitValues[std::min(iBand, nInitCount - 1)];
843 :
844 3044 : 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 2354 : if (CPLStringToComplex(pszBandInit, &adfInitRealImag[0],
867 2354 : &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 3042 : GByte *pBandData = static_cast<GByte *>(pDstBuffer) + iBand * nBandSize;
876 :
877 3042 : if (psOptions->eWorkingDataType == GDT_UInt8)
878 : {
879 4806 : memset(pBandData,
880 : std::max(
881 2403 : 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 2342 : 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 2653 : void GDALWarpOperation::DestroyDestinationBuffer(void *pDstBuffer)
916 : {
917 2653 : VSIFree(pDstBuffer);
918 2653 : }
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 1213 : 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 1213 : WipeChunkList();
966 1213 : CollectChunkListInternal(nDstXOff, nDstYOff, nDstXSize, nDstYSize);
967 :
968 : // Sort chunks from top to bottom, and for equal y, from left to right.
969 1213 : 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 1213 : const int knIntMax = std::numeric_limits<int>::max();
987 1213 : const int knIntMin = std::numeric_limits<int>::min();
988 1213 : int nSrcXOff = knIntMax;
989 1213 : int nSrcYOff = knIntMax;
990 1213 : int nSrcX2Off = knIntMin;
991 1213 : int nSrcY2Off = knIntMin;
992 1213 : double dfApproxAccArea = 0;
993 3521 : for (int iChunk = 0; pasChunkList != nullptr && iChunk < nChunkListCount;
994 : iChunk++)
995 : {
996 2308 : GDALWarpChunk *pasThisChunk = pasChunkList + iChunk;
997 2308 : nSrcXOff = std::min(nSrcXOff, pasThisChunk->sx);
998 2308 : nSrcYOff = std::min(nSrcYOff, pasThisChunk->sy);
999 2308 : nSrcX2Off = std::max(nSrcX2Off, pasThisChunk->sx + pasThisChunk->ssx);
1000 2308 : nSrcY2Off = std::max(nSrcY2Off, pasThisChunk->sy + pasThisChunk->ssy);
1001 2308 : dfApproxAccArea +=
1002 2308 : static_cast<double>(pasThisChunk->ssx) * pasThisChunk->ssy;
1003 : }
1004 1213 : if (nSrcXOff < nSrcX2Off)
1005 : {
1006 1205 : const double dfTotalArea =
1007 1205 : static_cast<double>(nSrcX2Off - nSrcXOff) * (nSrcY2Off - nSrcYOff);
1008 : // This is really a gross heuristics, but should work in most cases
1009 1205 : if (dfApproxAccArea >= dfTotalArea * 0.80)
1010 : {
1011 1205 : GDALDataset::FromHandle(psOptions->hSrcDS)
1012 1205 : ->AdviseRead(nSrcXOff, nSrcYOff, nSrcX2Off - nSrcXOff,
1013 : nSrcY2Off - nSrcYOff, nDstXSize, nDstYSize,
1014 1205 : psOptions->eWorkingDataType, psOptions->nBandCount,
1015 1205 : psOptions->panSrcBands, nullptr);
1016 : }
1017 : }
1018 1213 : }
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 1207 : 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 1207 : CollectChunkList(nDstXOff, nDstYOff, nDstXSize, nDstYSize);
1055 :
1056 : /* -------------------------------------------------------------------- */
1057 : /* Total up output pixels to process. */
1058 : /* -------------------------------------------------------------------- */
1059 1207 : double dfTotalPixels = 0.0;
1060 :
1061 3504 : for (int iChunk = 0; pasChunkList != nullptr && iChunk < nChunkListCount;
1062 : iChunk++)
1063 : {
1064 2297 : GDALWarpChunk *pasThisChunk = pasChunkList + iChunk;
1065 2297 : const double dfChunkPixels =
1066 2297 : pasThisChunk->dsx * static_cast<double>(pasThisChunk->dsy);
1067 :
1068 2297 : dfTotalPixels += dfChunkPixels;
1069 : }
1070 :
1071 : /* -------------------------------------------------------------------- */
1072 : /* Process them one at a time, updating the progress */
1073 : /* information for each region. */
1074 : /* -------------------------------------------------------------------- */
1075 1207 : double dfPixelsProcessed = 0.0;
1076 :
1077 3498 : for (int iChunk = 0; pasChunkList != nullptr && iChunk < nChunkListCount;
1078 : iChunk++)
1079 : {
1080 2297 : GDALWarpChunk *pasThisChunk = pasChunkList + iChunk;
1081 2297 : const double dfChunkPixels =
1082 2297 : pasThisChunk->dsx * static_cast<double>(pasThisChunk->dsy);
1083 :
1084 2297 : const double dfProgressBase = dfPixelsProcessed / dfTotalPixels;
1085 2297 : const double dfProgressScale = dfChunkPixels / dfTotalPixels;
1086 :
1087 2297 : 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 2297 : if (eErr != CE_None)
1094 6 : return eErr;
1095 :
1096 2291 : dfPixelsProcessed += dfChunkPixels;
1097 : }
1098 :
1099 1201 : WipeChunkList();
1100 :
1101 1201 : psOptions->pfnProgress(1.0, "", psOptions->pProgressArg);
1102 :
1103 1201 : 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 4168 : void GDALWarpOperation::WipeChunkList()
1378 :
1379 : {
1380 4168 : CPLFree(pasChunkList);
1381 4168 : pasChunkList = nullptr;
1382 4168 : nChunkListCount = 0;
1383 4168 : nChunkListMax = 0;
1384 4168 : }
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 4075 : 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 4075 : GDALGetDataTypeSizeBits(psOptions->eWorkingDataType) *
1405 4075 : psOptions->nBandCount;
1406 :
1407 4075 : if (psOptions->pfnSrcDensityMaskFunc != nullptr)
1408 0 : nSrcPixelCostInBits += 32; // Float mask?
1409 :
1410 4075 : GDALRasterBandH hSrcBand = nullptr;
1411 4075 : if (psOptions->nBandCount > 0)
1412 : hSrcBand =
1413 4075 : GDALGetRasterBand(psOptions->hSrcDS, psOptions->panSrcBands[0]);
1414 :
1415 4075 : if (psOptions->nSrcAlphaBand > 0 || psOptions->hCutline != nullptr)
1416 129 : nSrcPixelCostInBits += 32; // UnifiedSrcDensity float mask.
1417 7892 : else if (hSrcBand != nullptr &&
1418 3946 : (GDALGetMaskFlags(hSrcBand) & GMF_PER_DATASET))
1419 6 : nSrcPixelCostInBits += 1; // UnifiedSrcValid bit mask.
1420 :
1421 4075 : if (psOptions->papfnSrcPerBandValidityMaskFunc != nullptr ||
1422 4075 : psOptions->padfSrcNoDataReal != nullptr)
1423 202 : nSrcPixelCostInBits += psOptions->nBandCount; // Bit/band mask.
1424 :
1425 4075 : 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 4075 : GDALGetDataTypeSizeBits(psOptions->eWorkingDataType) *
1433 4075 : psOptions->nBandCount;
1434 :
1435 4075 : if (psOptions->pfnDstDensityMaskFunc != nullptr)
1436 0 : nDstPixelCostInBits += 32;
1437 :
1438 4075 : if (psOptions->padfDstNoDataReal != nullptr ||
1439 2413 : psOptions->pfnDstValidityMaskFunc != nullptr)
1440 1662 : nDstPixelCostInBits += psOptions->nBandCount;
1441 :
1442 4075 : if (psOptions->nDstAlphaBand > 0)
1443 271 : nDstPixelCostInBits += 32; // DstDensity float mask.
1444 :
1445 4075 : const double dfTotalMemoryUse =
1446 4075 : (static_cast<double>(nSrcPixelCostInBits) * nSrcXSize * nSrcYSize +
1447 4075 : static_cast<double>(nDstPixelCostInBits) * nDstXSize * nDstYSize) /
1448 : 8.0;
1449 4075 : return dfTotalMemoryUse;
1450 : }
1451 :
1452 : /************************************************************************/
1453 : /* CollectChunkListInternal() */
1454 : /************************************************************************/
1455 :
1456 4387 : 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 4387 : int nSrcXOff = 0;
1465 4387 : int nSrcYOff = 0;
1466 4387 : int nSrcXSize = 0;
1467 4387 : int nSrcYSize = 0;
1468 4387 : double dfSrcXExtraSize = 0.0;
1469 4387 : double dfSrcYExtraSize = 0.0;
1470 4387 : double dfSrcFillRatio = 0.0;
1471 : CPLErr eErr;
1472 : {
1473 4387 : CPLTurnFailureIntoWarningBackuper oBackuper;
1474 4387 : eErr = ComputeSourceWindow(nDstXOff, nDstYOff, nDstXSize, nDstYSize,
1475 : &nSrcXOff, &nSrcYOff, &nSrcXSize, &nSrcYSize,
1476 : &dfSrcXExtraSize, &dfSrcYExtraSize,
1477 : &dfSrcFillRatio);
1478 : }
1479 :
1480 4387 : 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 5512 : 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 3895 : 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 7790 : (nDstXSize > 2 || nDstYSize > 2)) ||
1533 3025 : (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 2308 : if (nChunkListCount == nChunkListMax)
1630 : {
1631 1373 : nChunkListMax = nChunkListMax * 2 + 1;
1632 1373 : pasChunkList = static_cast<GDALWarpChunk *>(
1633 1373 : CPLRealloc(pasChunkList, sizeof(GDALWarpChunk) * nChunkListMax));
1634 : }
1635 :
1636 2308 : pasChunkList[nChunkListCount].dx = nDstXOff;
1637 2308 : pasChunkList[nChunkListCount].dy = nDstYOff;
1638 2308 : pasChunkList[nChunkListCount].dsx = nDstXSize;
1639 2308 : pasChunkList[nChunkListCount].dsy = nDstYSize;
1640 2308 : pasChunkList[nChunkListCount].sx = nSrcXOff;
1641 2308 : pasChunkList[nChunkListCount].sy = nSrcYOff;
1642 2308 : pasChunkList[nChunkListCount].ssx = nSrcXSize;
1643 2308 : pasChunkList[nChunkListCount].ssy = nSrcYSize;
1644 2308 : pasChunkList[nChunkListCount].sExtraSx = dfSrcXExtraSize;
1645 2308 : pasChunkList[nChunkListCount].sExtraSy = dfSrcYExtraSize;
1646 :
1647 2308 : nChunkListCount++;
1648 :
1649 2308 : 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 2308 : 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 2308 : ReportTiming(nullptr);
1736 :
1737 : /* -------------------------------------------------------------------- */
1738 : /* Allocate the output buffer. */
1739 : /* -------------------------------------------------------------------- */
1740 2308 : int bDstBufferInitialized = FALSE;
1741 : void *pDstBuffer =
1742 2308 : CreateDestinationBuffer(nDstXSize, nDstYSize, &bDstBufferInitialized);
1743 2308 : 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 2306 : GDALDataset *poDstDS = GDALDataset::FromHandle(psOptions->hDstDS);
1753 2306 : 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 2306 : ? CE_None
1788 1905 : : WarpRegionToBuffer(
1789 : nDstXOff, nDstYOff, nDstXSize, nDstYSize,
1790 1905 : pDstBuffer, psOptions->eWorkingDataType, nSrcXOff,
1791 : nSrcYOff, nSrcXSize, nSrcYSize, dfSrcXExtraSize,
1792 2306 : dfSrcYExtraSize, dfProgressBase, dfProgressScale);
1793 :
1794 : /* -------------------------------------------------------------------- */
1795 : /* Write the output data back to disk if all went well. */
1796 : /* -------------------------------------------------------------------- */
1797 2306 : if (eErr == CE_None)
1798 : {
1799 2301 : if (psOptions->nBandCount == 1)
1800 : {
1801 : // Particular case to simplify the stack a bit.
1802 2082 : eErr = poDstDS->GetRasterBand(psOptions->panDstBands[0])
1803 4164 : ->RasterIO(GF_Write, nDstXOff, nDstYOff, nDstXSize,
1804 : nDstYSize, pDstBuffer, nDstXSize, nDstYSize,
1805 2082 : 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 4602 : if (eErr == CE_None &&
1817 2301 : 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 2301 : ReportTiming("Output buffer write");
1828 : }
1829 :
1830 : /* -------------------------------------------------------------------- */
1831 : /* Cleanup and return. */
1832 : /* -------------------------------------------------------------------- */
1833 2306 : DestroyDestinationBuffer(pDstBuffer);
1834 :
1835 2306 : 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 2997 : 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 2997 : const int nWordSize = GDALGetDataTypeSizeBytes(psOptions->eWorkingDataType);
1942 :
1943 2997 : CPLAssert(eBufDataType == psOptions->eWorkingDataType);
1944 :
1945 : /* -------------------------------------------------------------------- */
1946 : /* If not given a corresponding source window compute one now. */
1947 : /* -------------------------------------------------------------------- */
1948 2997 : 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 5922 : GDALWarpKernel oWK;
1980 :
1981 2961 : oWK.eResample = m_bIsTranslationOnPixelBoundaries ? GRA_NearestNeighbour
1982 2711 : : psOptions->eResampleAlg;
1983 2961 : oWK.eTieStrategy = psOptions->eTieStrategy;
1984 2961 : oWK.nBands = psOptions->nBandCount;
1985 2961 : oWK.eWorkingDataType = psOptions->eWorkingDataType;
1986 :
1987 2961 : oWK.pfnTransformer = psOptions->pfnTransformer;
1988 2961 : oWK.pTransformerArg = psOptions->pTransformerArg;
1989 :
1990 2961 : oWK.pfnProgress = psOptions->pfnProgress;
1991 2961 : oWK.pProgress = psOptions->pProgressArg;
1992 2961 : oWK.dfProgressBase = dfProgressBase;
1993 2961 : oWK.dfProgressScale = dfProgressScale;
1994 :
1995 2961 : oWK.papszWarpOptions = psOptions->papszWarpOptions;
1996 2961 : oWK.psThreadData = psThreadData;
1997 :
1998 2961 : 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 2961 : oWK.nSrcXOff = nSrcXOff;
2007 2961 : oWK.nSrcYOff = nSrcYOff;
2008 2961 : oWK.nSrcXSize = nSrcXSize;
2009 2961 : oWK.nSrcYSize = nSrcYSize;
2010 2961 : oWK.dfSrcXExtraSize = dfSrcXExtraSize;
2011 2961 : oWK.dfSrcYExtraSize = dfSrcYExtraSize;
2012 :
2013 2961 : GInt64 nAlloc64 =
2014 2961 : nWordSize *
2015 2961 : (static_cast<GInt64>(nSrcXSize) * nSrcYSize + WARP_EXTRA_ELTS) *
2016 2961 : psOptions->nBandCount;
2017 : #if SIZEOF_VOIDP == 4
2018 : if (nAlloc64 != static_cast<GInt64>(static_cast<size_t>(nAlloc64)))
2019 : {
2020 : CPLError(CE_Failure, CPLE_AppDefined,
2021 : "Integer overflow : nSrcXSize=%d, nSrcYSize=%d", nSrcXSize,
2022 : nSrcYSize);
2023 : return CE_Failure;
2024 : }
2025 : #endif
2026 :
2027 2961 : oWK.papabySrcImage = static_cast<GByte **>(
2028 2961 : CPLCalloc(sizeof(GByte *), psOptions->nBandCount));
2029 5922 : oWK.papabySrcImage[0] =
2030 2961 : static_cast<GByte *>(VSI_MALLOC_VERBOSE(static_cast<size_t>(nAlloc64)));
2031 :
2032 2961 : CPLErr eErr =
2033 2942 : nSrcXSize != 0 && nSrcYSize != 0 && oWK.papabySrcImage[0] == nullptr
2034 5903 : ? CE_Failure
2035 : : CE_None;
2036 :
2037 7869 : for (int i = 0; i < psOptions->nBandCount && eErr == CE_None; i++)
2038 4908 : oWK.papabySrcImage[i] =
2039 4908 : reinterpret_cast<GByte *>(oWK.papabySrcImage[0]) +
2040 4908 : nWordSize *
2041 4908 : (static_cast<GPtrDiff_t>(nSrcXSize) * nSrcYSize +
2042 4908 : WARP_EXTRA_ELTS) *
2043 4908 : i;
2044 :
2045 2961 : if (eErr == CE_None && nSrcXSize > 0 && nSrcYSize > 0)
2046 : {
2047 2936 : GDALDataset *poSrcDS = GDALDataset::FromHandle(psOptions->hSrcDS);
2048 2936 : if (psOptions->nBandCount == 1)
2049 : {
2050 : // Particular case to simplify the stack a bit.
2051 1951 : eErr = poSrcDS->GetRasterBand(psOptions->panSrcBands[0])
2052 3902 : ->RasterIO(GF_Read, nSrcXOff, nSrcYOff, nSrcXSize,
2053 1951 : nSrcYSize, oWK.papabySrcImage[0], nSrcXSize,
2054 1951 : nSrcYSize, psOptions->eWorkingDataType, 0, 0,
2055 : nullptr);
2056 : }
2057 : else
2058 : {
2059 985 : eErr = poSrcDS->RasterIO(
2060 : GF_Read, nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
2061 985 : oWK.papabySrcImage[0], nSrcXSize, nSrcYSize,
2062 985 : psOptions->eWorkingDataType, psOptions->nBandCount,
2063 985 : psOptions->panSrcBands, 0, 0,
2064 985 : nWordSize * (static_cast<GPtrDiff_t>(nSrcXSize) * nSrcYSize +
2065 : WARP_EXTRA_ELTS),
2066 : nullptr);
2067 : }
2068 : }
2069 :
2070 2961 : ReportTiming("Input buffer read");
2071 :
2072 : /* -------------------------------------------------------------------- */
2073 : /* Initialize destination buffer. */
2074 : /* -------------------------------------------------------------------- */
2075 2961 : oWK.nDstXOff = nDstXOff;
2076 2961 : oWK.nDstYOff = nDstYOff;
2077 2961 : oWK.nDstXSize = nDstXSize;
2078 2961 : oWK.nDstYSize = nDstYSize;
2079 :
2080 2961 : oWK.papabyDstImage = reinterpret_cast<GByte **>(
2081 2961 : CPLCalloc(sizeof(GByte *), psOptions->nBandCount));
2082 :
2083 7861 : for (int i = 0; i < psOptions->nBandCount && eErr == CE_None; i++)
2084 : {
2085 4900 : oWK.papabyDstImage[i] =
2086 4900 : static_cast<GByte *>(pDataBuf) +
2087 4900 : i * static_cast<GPtrDiff_t>(nDstXSize) * nDstYSize * nWordSize;
2088 : }
2089 :
2090 : /* -------------------------------------------------------------------- */
2091 : /* Eventually we need handling for a whole bunch of the */
2092 : /* validity and density masks here. */
2093 : /* -------------------------------------------------------------------- */
2094 :
2095 : // TODO
2096 :
2097 : /* -------------------------------------------------------------------- */
2098 : /* Generate a source density mask if we have a source alpha band */
2099 : /* -------------------------------------------------------------------- */
2100 2961 : if (eErr == CE_None && psOptions->nSrcAlphaBand > 0 && nSrcXSize > 0 &&
2101 109 : nSrcYSize > 0)
2102 : {
2103 109 : CPLAssert(oWK.pafUnifiedSrcDensity == nullptr);
2104 :
2105 109 : eErr = CreateKernelMask(&oWK, 0 /* not used */, "UnifiedSrcDensity");
2106 :
2107 109 : if (eErr == CE_None)
2108 : {
2109 109 : int bOutAllOpaque = FALSE;
2110 218 : eErr = GDALWarpSrcAlphaMasker(
2111 109 : psOptions, psOptions->nBandCount, psOptions->eWorkingDataType,
2112 : oWK.nSrcXOff, oWK.nSrcYOff, oWK.nSrcXSize, oWK.nSrcYSize,
2113 109 : oWK.papabySrcImage, TRUE, oWK.pafUnifiedSrcDensity,
2114 : &bOutAllOpaque);
2115 109 : if (bOutAllOpaque)
2116 : {
2117 : #if DEBUG_VERBOSE
2118 : CPLDebug("WARP",
2119 : "No need for a source density mask as all values "
2120 : "are opaque");
2121 : #endif
2122 36 : CPLFree(oWK.pafUnifiedSrcDensity);
2123 36 : oWK.pafUnifiedSrcDensity = nullptr;
2124 : }
2125 : }
2126 : }
2127 :
2128 : /* -------------------------------------------------------------------- */
2129 : /* Generate a source density mask if we have a source cutline. */
2130 : /* -------------------------------------------------------------------- */
2131 2961 : if (eErr == CE_None && psOptions->hCutline != nullptr && nSrcXSize > 0 &&
2132 39 : nSrcYSize > 0)
2133 : {
2134 39 : const bool bUnifiedSrcDensityJustCreated =
2135 39 : (oWK.pafUnifiedSrcDensity == nullptr);
2136 39 : if (bUnifiedSrcDensityJustCreated)
2137 : {
2138 : eErr =
2139 39 : CreateKernelMask(&oWK, 0 /* not used */, "UnifiedSrcDensity");
2140 :
2141 39 : if (eErr == CE_None)
2142 : {
2143 39 : for (GPtrDiff_t j = 0;
2144 2879360 : j < static_cast<GPtrDiff_t>(oWK.nSrcXSize) * oWK.nSrcYSize;
2145 : j++)
2146 2879320 : oWK.pafUnifiedSrcDensity[j] = 1.0;
2147 : }
2148 : }
2149 :
2150 39 : int nValidityFlag = 0;
2151 39 : if (eErr == CE_None)
2152 39 : eErr = GDALWarpCutlineMaskerEx(
2153 39 : psOptions, psOptions->nBandCount, psOptions->eWorkingDataType,
2154 : oWK.nSrcXOff, oWK.nSrcYOff, oWK.nSrcXSize, oWK.nSrcYSize,
2155 39 : oWK.papabySrcImage, TRUE, oWK.pafUnifiedSrcDensity,
2156 : &nValidityFlag);
2157 39 : if (nValidityFlag == GCMVF_CHUNK_FULLY_WITHIN_CUTLINE &&
2158 : bUnifiedSrcDensityJustCreated)
2159 : {
2160 8 : VSIFree(oWK.pafUnifiedSrcDensity);
2161 8 : oWK.pafUnifiedSrcDensity = nullptr;
2162 : }
2163 : }
2164 :
2165 : /* -------------------------------------------------------------------- */
2166 : /* Generate a destination density mask if we have a destination */
2167 : /* alpha band. */
2168 : /* -------------------------------------------------------------------- */
2169 2961 : if (eErr == CE_None && psOptions->nDstAlphaBand > 0)
2170 : {
2171 906 : CPLAssert(oWK.pafDstDensity == nullptr);
2172 :
2173 906 : eErr = CreateKernelMask(&oWK, 0 /* not used */, "DstDensity");
2174 :
2175 906 : if (eErr == CE_None)
2176 906 : eErr = GDALWarpDstAlphaMasker(
2177 906 : psOptions, psOptions->nBandCount, psOptions->eWorkingDataType,
2178 : oWK.nDstXOff, oWK.nDstYOff, oWK.nDstXSize, oWK.nDstYSize,
2179 906 : oWK.papabyDstImage, TRUE, oWK.pafDstDensity);
2180 : }
2181 :
2182 : /* -------------------------------------------------------------------- */
2183 : /* If we have source nodata values create the validity mask. */
2184 : /* -------------------------------------------------------------------- */
2185 2961 : if (eErr == CE_None && psOptions->padfSrcNoDataReal != nullptr &&
2186 193 : nSrcXSize > 0 && nSrcYSize > 0)
2187 : {
2188 184 : CPLAssert(oWK.papanBandSrcValid == nullptr);
2189 :
2190 184 : bool bAllBandsAllValid = true;
2191 443 : for (int i = 0; i < psOptions->nBandCount && eErr == CE_None; i++)
2192 : {
2193 259 : eErr = CreateKernelMask(&oWK, i, "BandSrcValid");
2194 259 : if (eErr == CE_None)
2195 : {
2196 259 : double adfNoData[2] = {psOptions->padfSrcNoDataReal[i],
2197 259 : psOptions->padfSrcNoDataImag != nullptr
2198 259 : ? psOptions->padfSrcNoDataImag[i]
2199 259 : : 0.0};
2200 :
2201 259 : int bAllValid = FALSE;
2202 518 : eErr = GDALWarpNoDataMasker(
2203 259 : adfNoData, 1, psOptions->eWorkingDataType, oWK.nSrcXOff,
2204 : oWK.nSrcYOff, oWK.nSrcXSize, oWK.nSrcYSize,
2205 259 : &(oWK.papabySrcImage[i]), FALSE, oWK.papanBandSrcValid[i],
2206 : &bAllValid);
2207 259 : if (!bAllValid)
2208 168 : bAllBandsAllValid = false;
2209 : }
2210 : }
2211 :
2212 : // Optimization: if all pixels in all bands are valid,
2213 : // we don't need a mask.
2214 184 : if (bAllBandsAllValid)
2215 : {
2216 : #if DEBUG_VERBOSE
2217 : CPLDebug(
2218 : "WARP",
2219 : "No need for a source nodata mask as all values are valid");
2220 : #endif
2221 169 : for (int k = 0; k < psOptions->nBandCount; k++)
2222 90 : CPLFree(oWK.papanBandSrcValid[k]);
2223 79 : CPLFree(oWK.papanBandSrcValid);
2224 79 : oWK.papanBandSrcValid = nullptr;
2225 : }
2226 :
2227 : /* --------------------------------------------------------------------
2228 : */
2229 : /* If there's just a single band, then transfer */
2230 : /* papanBandSrcValid[0] as panUnifiedSrcValid. */
2231 : /* --------------------------------------------------------------------
2232 : */
2233 184 : if (oWK.papanBandSrcValid != nullptr && psOptions->nBandCount == 1)
2234 : {
2235 62 : oWK.panUnifiedSrcValid = oWK.papanBandSrcValid[0];
2236 62 : CPLFree(oWK.papanBandSrcValid);
2237 62 : oWK.papanBandSrcValid = nullptr;
2238 : }
2239 :
2240 : /* --------------------------------------------------------------------
2241 : */
2242 : /* Compute a unified input pixel mask if and only if all bands */
2243 : /* nodata is true. That is, we only treat a pixel as nodata if */
2244 : /* all bands match their respective nodata values. */
2245 : /* --------------------------------------------------------------------
2246 : */
2247 122 : else if (oWK.papanBandSrcValid != nullptr && eErr == CE_None)
2248 : {
2249 43 : bool bAtLeastOneBandAllValid = false;
2250 150 : for (int k = 0; k < psOptions->nBandCount; k++)
2251 : {
2252 107 : if (oWK.papanBandSrcValid[k] == nullptr)
2253 : {
2254 0 : bAtLeastOneBandAllValid = true;
2255 0 : break;
2256 : }
2257 : }
2258 :
2259 86 : const char *pszUnifiedSrcNoData = CSLFetchNameValue(
2260 43 : psOptions->papszWarpOptions, "UNIFIED_SRC_NODATA");
2261 81 : if (!bAtLeastOneBandAllValid && (pszUnifiedSrcNoData == nullptr ||
2262 38 : CPLTestBool(pszUnifiedSrcNoData)))
2263 : {
2264 42 : auto nMaskBits =
2265 42 : static_cast<GPtrDiff_t>(oWK.nSrcXSize) * oWK.nSrcYSize;
2266 :
2267 : eErr =
2268 42 : CreateKernelMask(&oWK, 0 /* not used */, "UnifiedSrcValid");
2269 :
2270 42 : if (eErr == CE_None)
2271 : {
2272 42 : CPLMaskClearAll(oWK.panUnifiedSrcValid, nMaskBits);
2273 :
2274 146 : for (int k = 0; k < psOptions->nBandCount; k++)
2275 : {
2276 104 : CPLMaskMerge(oWK.panUnifiedSrcValid,
2277 104 : oWK.papanBandSrcValid[k], nMaskBits);
2278 : }
2279 :
2280 : // If UNIFIED_SRC_NODATA is set, then we will ignore the
2281 : // individual nodata status of each band. If it is not set,
2282 : // both mechanism apply:
2283 : // - if panUnifiedSrcValid[] indicates a pixel is invalid
2284 : // (that is all its bands are at nodata), then the output
2285 : // pixel will be invalid
2286 : // - otherwise, the status band per band will be check with
2287 : // papanBandSrcValid[iBand][], and the output pixel will
2288 : // be valid
2289 42 : if (pszUnifiedSrcNoData != nullptr &&
2290 37 : !EQUAL(pszUnifiedSrcNoData, "PARTIAL"))
2291 : {
2292 123 : for (int k = 0; k < psOptions->nBandCount; k++)
2293 87 : CPLFree(oWK.papanBandSrcValid[k]);
2294 36 : CPLFree(oWK.papanBandSrcValid);
2295 36 : oWK.papanBandSrcValid = nullptr;
2296 : }
2297 : }
2298 : }
2299 : }
2300 : }
2301 :
2302 : /* -------------------------------------------------------------------- */
2303 : /* Generate a source validity mask if we have a source mask for */
2304 : /* the whole input dataset (and didn't already treat it as */
2305 : /* alpha band). */
2306 : /* -------------------------------------------------------------------- */
2307 : GDALRasterBandH hSrcBand =
2308 2961 : psOptions->nBandCount < 1
2309 2961 : ? nullptr
2310 2961 : : GDALGetRasterBand(psOptions->hSrcDS, psOptions->panSrcBands[0]);
2311 :
2312 2957 : if (eErr == CE_None && oWK.pafUnifiedSrcDensity == nullptr &&
2313 2853 : oWK.panUnifiedSrcValid == nullptr && psOptions->nSrcAlphaBand <= 0 &&
2314 2714 : (GDALGetMaskFlags(hSrcBand) & GMF_PER_DATASET)
2315 : // Need to double check for -nosrcalpha case.
2316 13 : && !(GDALGetMaskFlags(hSrcBand) & GMF_ALPHA) &&
2317 5920 : psOptions->padfSrcNoDataReal == nullptr && nSrcXSize > 0 &&
2318 2 : nSrcYSize > 0)
2319 :
2320 : {
2321 2 : eErr = CreateKernelMask(&oWK, 0 /* not used */, "UnifiedSrcValid");
2322 :
2323 2 : if (eErr == CE_None)
2324 2 : eErr = GDALWarpSrcMaskMasker(
2325 2 : psOptions, psOptions->nBandCount, psOptions->eWorkingDataType,
2326 : oWK.nSrcXOff, oWK.nSrcYOff, oWK.nSrcXSize, oWK.nSrcYSize,
2327 2 : oWK.papabySrcImage, FALSE, oWK.panUnifiedSrcValid);
2328 : }
2329 :
2330 : /* -------------------------------------------------------------------- */
2331 : /* If we have destination nodata values create the */
2332 : /* validity mask. We set the DstValid for any pixel that we */
2333 : /* do no have valid data in *any* of the source bands. */
2334 : /* */
2335 : /* Note that we don't support any concept of unified nodata on */
2336 : /* the destination image. At some point that should be added */
2337 : /* and then this logic will be significantly different. */
2338 : /* -------------------------------------------------------------------- */
2339 2961 : if (eErr == CE_None && psOptions->padfDstNoDataReal != nullptr)
2340 : {
2341 568 : CPLAssert(oWK.panDstValid == nullptr);
2342 :
2343 568 : const GPtrDiff_t nMaskBits =
2344 568 : static_cast<GPtrDiff_t>(oWK.nDstXSize) * oWK.nDstYSize;
2345 :
2346 568 : eErr = CreateKernelMask(&oWK, 0 /* not used */, "DstValid");
2347 : GUInt32 *panBandMask =
2348 568 : eErr == CE_None ? CPLMaskCreate(nMaskBits, true) : nullptr;
2349 :
2350 568 : if (eErr == CE_None && panBandMask != nullptr)
2351 : {
2352 1203 : for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
2353 : {
2354 693 : CPLMaskSetAll(panBandMask, nMaskBits);
2355 :
2356 693 : double adfNoData[2] = {psOptions->padfDstNoDataReal[iBand],
2357 693 : psOptions->padfDstNoDataImag != nullptr
2358 693 : ? psOptions->padfDstNoDataImag[iBand]
2359 693 : : 0.0};
2360 :
2361 693 : int bAllValid = FALSE;
2362 1386 : eErr = GDALWarpNoDataMasker(
2363 693 : adfNoData, 1, psOptions->eWorkingDataType, oWK.nDstXOff,
2364 : oWK.nDstYOff, oWK.nDstXSize, oWK.nDstYSize,
2365 693 : oWK.papabyDstImage + iBand, FALSE, panBandMask, &bAllValid);
2366 :
2367 : // Optimization: if there's a single band and all pixels are
2368 : // valid then we don't need a mask.
2369 693 : if (bAllValid && psOptions->nBandCount == 1)
2370 : {
2371 : #if DEBUG_VERBOSE
2372 : CPLDebug("WARP", "No need for a destination nodata mask as "
2373 : "all values are valid");
2374 : #endif
2375 58 : CPLFree(oWK.panDstValid);
2376 58 : oWK.panDstValid = nullptr;
2377 58 : break;
2378 : }
2379 :
2380 635 : CPLMaskMerge(oWK.panDstValid, panBandMask, nMaskBits);
2381 : }
2382 568 : CPLFree(panBandMask);
2383 : }
2384 : }
2385 :
2386 : /* -------------------------------------------------------------------- */
2387 : /* Release IO Mutex, and acquire warper mutex. */
2388 : /* -------------------------------------------------------------------- */
2389 2961 : if (hIOMutex != nullptr)
2390 : {
2391 11 : CPLReleaseMutex(hIOMutex);
2392 11 : if (!CPLAcquireMutex(hWarpMutex, 600.0))
2393 : {
2394 0 : CPLError(CE_Failure, CPLE_AppDefined,
2395 : "Failed to acquire WarpMutex in WarpRegion().");
2396 0 : return CE_Failure;
2397 : }
2398 : }
2399 :
2400 : /* -------------------------------------------------------------------- */
2401 : /* Optional application provided prewarp chunk processor. */
2402 : /* -------------------------------------------------------------------- */
2403 2961 : if (eErr == CE_None && psOptions->pfnPreWarpChunkProcessor != nullptr)
2404 0 : eErr = psOptions->pfnPreWarpChunkProcessor(
2405 0 : &oWK, psOptions->pPreWarpProcessorArg);
2406 :
2407 : /* -------------------------------------------------------------------- */
2408 : /* Perform the warp. */
2409 : /* -------------------------------------------------------------------- */
2410 2961 : if (eErr == CE_None)
2411 : {
2412 2957 : eErr = oWK.PerformWarp();
2413 2957 : ReportTiming("In memory warp operation");
2414 : }
2415 :
2416 : /* -------------------------------------------------------------------- */
2417 : /* Optional application provided postwarp chunk processor. */
2418 : /* -------------------------------------------------------------------- */
2419 2961 : if (eErr == CE_None && psOptions->pfnPostWarpChunkProcessor != nullptr)
2420 0 : eErr = psOptions->pfnPostWarpChunkProcessor(
2421 0 : &oWK, psOptions->pPostWarpProcessorArg);
2422 :
2423 : /* -------------------------------------------------------------------- */
2424 : /* Release Warp Mutex, and acquire io mutex. */
2425 : /* -------------------------------------------------------------------- */
2426 2961 : if (hIOMutex != nullptr)
2427 : {
2428 11 : CPLReleaseMutex(hWarpMutex);
2429 11 : if (!CPLAcquireMutex(hIOMutex, 600.0))
2430 : {
2431 0 : CPLError(CE_Failure, CPLE_AppDefined,
2432 : "Failed to acquire IOMutex in WarpRegion().");
2433 0 : return CE_Failure;
2434 : }
2435 : }
2436 :
2437 : /* -------------------------------------------------------------------- */
2438 : /* Write destination alpha if available. */
2439 : /* -------------------------------------------------------------------- */
2440 2961 : if (eErr == CE_None && psOptions->nDstAlphaBand > 0)
2441 : {
2442 906 : eErr = GDALWarpDstAlphaMasker(
2443 906 : psOptions, -psOptions->nBandCount, psOptions->eWorkingDataType,
2444 : oWK.nDstXOff, oWK.nDstYOff, oWK.nDstXSize, oWK.nDstYSize,
2445 906 : oWK.papabyDstImage, TRUE, oWK.pafDstDensity);
2446 : }
2447 :
2448 : /* -------------------------------------------------------------------- */
2449 : /* Cleanup. */
2450 : /* -------------------------------------------------------------------- */
2451 2961 : CPLFree(oWK.papabySrcImage[0]);
2452 2961 : CPLFree(oWK.papabySrcImage);
2453 2961 : CPLFree(oWK.papabyDstImage);
2454 :
2455 2961 : if (oWK.papanBandSrcValid != nullptr)
2456 : {
2457 27 : for (int i = 0; i < oWK.nBands; i++)
2458 20 : CPLFree(oWK.papanBandSrcValid[i]);
2459 7 : CPLFree(oWK.papanBandSrcValid);
2460 : }
2461 2961 : CPLFree(oWK.panUnifiedSrcValid);
2462 2961 : CPLFree(oWK.pafUnifiedSrcDensity);
2463 2961 : CPLFree(oWK.panDstValid);
2464 2961 : CPLFree(oWK.pafDstDensity);
2465 :
2466 2961 : return eErr;
2467 : }
2468 :
2469 : /************************************************************************/
2470 : /* GDALWarpRegionToBuffer() */
2471 : /************************************************************************/
2472 :
2473 : /**
2474 : * @see GDALWarpOperation::WarpRegionToBuffer()
2475 : */
2476 :
2477 0 : CPLErr GDALWarpRegionToBuffer(GDALWarpOperationH hOperation, int nDstXOff,
2478 : int nDstYOff, int nDstXSize, int nDstYSize,
2479 : void *pDataBuf, GDALDataType eBufDataType,
2480 : int nSrcXOff, int nSrcYOff, int nSrcXSize,
2481 : int nSrcYSize)
2482 :
2483 : {
2484 0 : VALIDATE_POINTER1(hOperation, "GDALWarpRegionToBuffer", CE_Failure);
2485 :
2486 : return reinterpret_cast<GDALWarpOperation *>(hOperation)
2487 0 : ->WarpRegionToBuffer(nDstXOff, nDstYOff, nDstXSize, nDstYSize, pDataBuf,
2488 : eBufDataType, nSrcXOff, nSrcYOff, nSrcXSize,
2489 0 : nSrcYSize);
2490 : }
2491 :
2492 : /************************************************************************/
2493 : /* CreateKernelMask() */
2494 : /* */
2495 : /* If mask does not yet exist, create it. Supported types are */
2496 : /* the name of the variable in question. That is */
2497 : /* "BandSrcValid", "UnifiedSrcValid", "UnifiedSrcDensity", */
2498 : /* "DstValid", and "DstDensity". */
2499 : /************************************************************************/
2500 :
2501 1925 : CPLErr GDALWarpOperation::CreateKernelMask(GDALWarpKernel *poKernel, int iBand,
2502 : const char *pszType)
2503 :
2504 : {
2505 1925 : void **ppMask = nullptr;
2506 1925 : int nXSize = 0;
2507 1925 : int nYSize = 0;
2508 1925 : int nBitsPerPixel = 0;
2509 1925 : int nDefault = 0;
2510 1925 : int nExtraElts = 0;
2511 1925 : bool bDoMemset = true;
2512 :
2513 : /* -------------------------------------------------------------------- */
2514 : /* Get particulars of mask to be updated. */
2515 : /* -------------------------------------------------------------------- */
2516 1925 : if (EQUAL(pszType, "BandSrcValid"))
2517 : {
2518 259 : if (poKernel->papanBandSrcValid == nullptr)
2519 184 : poKernel->papanBandSrcValid = static_cast<GUInt32 **>(
2520 184 : CPLCalloc(sizeof(void *), poKernel->nBands));
2521 :
2522 259 : ppMask =
2523 259 : reinterpret_cast<void **>(&(poKernel->papanBandSrcValid[iBand]));
2524 259 : nExtraElts = WARP_EXTRA_ELTS;
2525 259 : nXSize = poKernel->nSrcXSize;
2526 259 : nYSize = poKernel->nSrcYSize;
2527 259 : nBitsPerPixel = 1;
2528 259 : nDefault = 0xff;
2529 : }
2530 1666 : else if (EQUAL(pszType, "UnifiedSrcValid"))
2531 : {
2532 44 : ppMask = reinterpret_cast<void **>(&(poKernel->panUnifiedSrcValid));
2533 44 : nExtraElts = WARP_EXTRA_ELTS;
2534 44 : nXSize = poKernel->nSrcXSize;
2535 44 : nYSize = poKernel->nSrcYSize;
2536 44 : nBitsPerPixel = 1;
2537 44 : nDefault = 0xff;
2538 : }
2539 1622 : else if (EQUAL(pszType, "UnifiedSrcDensity"))
2540 : {
2541 148 : ppMask = reinterpret_cast<void **>(&(poKernel->pafUnifiedSrcDensity));
2542 148 : nExtraElts = WARP_EXTRA_ELTS;
2543 148 : nXSize = poKernel->nSrcXSize;
2544 148 : nYSize = poKernel->nSrcYSize;
2545 148 : nBitsPerPixel = 32;
2546 148 : nDefault = 0;
2547 148 : bDoMemset = false;
2548 : }
2549 1474 : else if (EQUAL(pszType, "DstValid"))
2550 : {
2551 568 : ppMask = reinterpret_cast<void **>(&(poKernel->panDstValid));
2552 568 : nXSize = poKernel->nDstXSize;
2553 568 : nYSize = poKernel->nDstYSize;
2554 568 : nBitsPerPixel = 1;
2555 568 : nDefault = 0;
2556 : }
2557 906 : else if (EQUAL(pszType, "DstDensity"))
2558 : {
2559 906 : ppMask = reinterpret_cast<void **>(&(poKernel->pafDstDensity));
2560 906 : nXSize = poKernel->nDstXSize;
2561 906 : nYSize = poKernel->nDstYSize;
2562 906 : nBitsPerPixel = 32;
2563 906 : nDefault = 0;
2564 906 : bDoMemset = false;
2565 : }
2566 : else
2567 : {
2568 0 : CPLError(CE_Failure, CPLE_AppDefined,
2569 : "Internal error in CreateKernelMask(%s).", pszType);
2570 0 : return CE_Failure;
2571 : }
2572 :
2573 : /* -------------------------------------------------------------------- */
2574 : /* Allocate if needed. */
2575 : /* -------------------------------------------------------------------- */
2576 1925 : if (*ppMask == nullptr)
2577 : {
2578 1925 : const GIntBig nBytes =
2579 : nBitsPerPixel == 32
2580 1925 : ? (static_cast<GIntBig>(nXSize) * nYSize + nExtraElts) * 4
2581 871 : : (static_cast<GIntBig>(nXSize) * nYSize + nExtraElts + 31) / 8;
2582 :
2583 1925 : const size_t nByteSize_t = static_cast<size_t>(nBytes);
2584 : #if SIZEOF_VOIDP == 4
2585 : if (static_cast<GIntBig>(nByteSize_t) != nBytes)
2586 : {
2587 : CPLError(CE_Failure, CPLE_OutOfMemory,
2588 : "Cannot allocate " CPL_FRMT_GIB " bytes", nBytes);
2589 : return CE_Failure;
2590 : }
2591 : #endif
2592 :
2593 1925 : *ppMask = VSI_MALLOC_VERBOSE(nByteSize_t);
2594 :
2595 1925 : if (*ppMask == nullptr)
2596 : {
2597 0 : return CE_Failure;
2598 : }
2599 :
2600 1925 : if (bDoMemset)
2601 871 : memset(*ppMask, nDefault, nByteSize_t);
2602 : }
2603 :
2604 1925 : return CE_None;
2605 : }
2606 :
2607 : /************************************************************************/
2608 : /* ComputeSourceWindowStartingFromSource() */
2609 : /************************************************************************/
2610 :
2611 : constexpr int DEFAULT_STEP_COUNT = 21;
2612 :
2613 1104 : void GDALWarpOperation::ComputeSourceWindowStartingFromSource(
2614 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize,
2615 : double *padfSrcMinX, double *padfSrcMinY, double *padfSrcMaxX,
2616 : double *padfSrcMaxY)
2617 : {
2618 1104 : const int nSrcRasterXSize = GDALGetRasterXSize(psOptions->hSrcDS);
2619 1104 : const int nSrcRasterYSize = GDALGetRasterYSize(psOptions->hSrcDS);
2620 1104 : if (nSrcRasterXSize == 0 || nSrcRasterYSize == 0)
2621 0 : return;
2622 :
2623 1104 : GDALWarpPrivateData *privateData = GetWarpPrivateData(this);
2624 1104 : if (privateData->nStepCount == 0)
2625 : {
2626 260 : int nStepCount = DEFAULT_STEP_COUNT;
2627 260 : std::vector<double> adfDstZ{};
2628 :
2629 : const char *pszSampleSteps =
2630 260 : CSLFetchNameValue(psOptions->papszWarpOptions, "SAMPLE_STEPS");
2631 260 : constexpr int knIntMax = std::numeric_limits<int>::max();
2632 260 : if (pszSampleSteps && !EQUAL(pszSampleSteps, "ALL"))
2633 : {
2634 0 : nStepCount = atoi(
2635 0 : CSLFetchNameValue(psOptions->papszWarpOptions, "SAMPLE_STEPS"));
2636 0 : nStepCount = std::max(2, nStepCount);
2637 : }
2638 :
2639 260 : const double dfStepSize = 1.0 / (nStepCount - 1);
2640 260 : if (nStepCount > knIntMax - 2 ||
2641 260 : (nStepCount + 2) > knIntMax / (nStepCount + 2))
2642 : {
2643 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps : %d",
2644 : nStepCount);
2645 0 : return;
2646 : }
2647 260 : const int nSampleMax = (nStepCount + 2) * (nStepCount + 2);
2648 :
2649 : try
2650 : {
2651 260 : privateData->abSuccess.resize(nSampleMax);
2652 260 : privateData->adfDstX.resize(nSampleMax);
2653 260 : privateData->adfDstY.resize(nSampleMax);
2654 260 : adfDstZ.resize(nSampleMax);
2655 : }
2656 0 : catch (const std::exception &)
2657 : {
2658 0 : return;
2659 : }
2660 :
2661 : /* --------------------------------------------------------------------
2662 : */
2663 : /* Setup sample points on a grid pattern throughout the source */
2664 : /* raster. */
2665 : /* --------------------------------------------------------------------
2666 : */
2667 260 : int iPoint = 0;
2668 6240 : for (int iY = 0; iY < nStepCount + 2; iY++)
2669 : {
2670 11700 : const double dfRatioY = (iY == 0) ? 0.5 / nSrcRasterYSize
2671 5720 : : (iY <= nStepCount)
2672 5720 : ? (iY - 1) * dfStepSize
2673 260 : : 1 - 0.5 / nSrcRasterYSize;
2674 143520 : for (int iX = 0; iX < nStepCount + 2; iX++)
2675 : {
2676 269100 : const double dfRatioX = (iX == 0) ? 0.5 / nSrcRasterXSize
2677 131560 : : (iX <= nStepCount)
2678 131560 : ? (iX - 1) * dfStepSize
2679 5980 : : 1 - 0.5 / nSrcRasterXSize;
2680 137540 : privateData->adfDstX[iPoint] = dfRatioX * nSrcRasterXSize;
2681 137540 : privateData->adfDstY[iPoint] = dfRatioY * nSrcRasterYSize;
2682 137540 : iPoint++;
2683 : }
2684 : }
2685 260 : CPLAssert(iPoint == nSampleMax);
2686 :
2687 : /* --------------------------------------------------------------------
2688 : */
2689 : /* Transform them to the output pixel coordinate space */
2690 : /* --------------------------------------------------------------------
2691 : */
2692 260 : psOptions->pfnTransformer(psOptions->pTransformerArg, FALSE, nSampleMax,
2693 : privateData->adfDstX.data(),
2694 : privateData->adfDstY.data(), adfDstZ.data(),
2695 : privateData->abSuccess.data());
2696 260 : privateData->nStepCount = nStepCount;
2697 : }
2698 :
2699 : /* -------------------------------------------------------------------- */
2700 : /* Collect the bounds, ignoring any failed points. */
2701 : /* -------------------------------------------------------------------- */
2702 1104 : const int nStepCount = privateData->nStepCount;
2703 1104 : const double dfStepSize = 1.0 / (nStepCount - 1);
2704 1104 : int iPoint = 0;
2705 : #ifdef DEBUG
2706 1104 : const size_t nSampleMax =
2707 1104 : static_cast<size_t>(nStepCount + 2) * (nStepCount + 2);
2708 1104 : CPL_IGNORE_RET_VAL(nSampleMax);
2709 1104 : CPLAssert(privateData->adfDstX.size() == nSampleMax);
2710 1104 : CPLAssert(privateData->adfDstY.size() == nSampleMax);
2711 1104 : CPLAssert(privateData->abSuccess.size() == nSampleMax);
2712 : #endif
2713 26496 : for (int iY = 0; iY < nStepCount + 2; iY++)
2714 : {
2715 49680 : const double dfRatioY = (iY == 0) ? 0.5 / nSrcRasterYSize
2716 : : (iY <= nStepCount)
2717 24288 : ? (iY - 1) * dfStepSize
2718 1104 : : 1 - 0.5 / nSrcRasterYSize;
2719 609408 : for (int iX = 0; iX < nStepCount + 2; iX++)
2720 : {
2721 584016 : if (privateData->abSuccess[iPoint] &&
2722 548544 : privateData->adfDstX[iPoint] >= nDstXOff &&
2723 404998 : privateData->adfDstX[iPoint] <= nDstXOff + nDstXSize &&
2724 1332800 : privateData->adfDstY[iPoint] >= nDstYOff &&
2725 200241 : privateData->adfDstY[iPoint] <= nDstYOff + nDstYSize)
2726 : {
2727 289168 : const double dfRatioX = (iX == 0) ? 0.5 / nSrcRasterXSize
2728 : : (iX <= nStepCount)
2729 141410 : ? (iX - 1) * dfStepSize
2730 6345 : : 1 - 0.5 / nSrcRasterXSize;
2731 147758 : double dfSrcX = dfRatioX * nSrcRasterXSize;
2732 147758 : double dfSrcY = dfRatioY * nSrcRasterYSize;
2733 147758 : *padfSrcMinX = std::min(*padfSrcMinX, dfSrcX);
2734 147758 : *padfSrcMinY = std::min(*padfSrcMinY, dfSrcY);
2735 147758 : *padfSrcMaxX = std::max(*padfSrcMaxX, dfSrcX);
2736 147758 : *padfSrcMaxY = std::max(*padfSrcMaxY, dfSrcY);
2737 : }
2738 584016 : iPoint++;
2739 : }
2740 : }
2741 : }
2742 :
2743 : /************************************************************************/
2744 : /* ComputeSourceWindowTransformPoints() */
2745 : /************************************************************************/
2746 :
2747 5794 : bool GDALWarpOperation::ComputeSourceWindowTransformPoints(
2748 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize, bool bUseGrid,
2749 : bool bAll, int nStepCount, bool bTryWithCheckWithInvertProj,
2750 : double &dfMinXOut, double &dfMinYOut, double &dfMaxXOut, double &dfMaxYOut,
2751 : int &nSamplePoints, int &nFailedCount)
2752 : {
2753 5794 : nSamplePoints = 0;
2754 5794 : nFailedCount = 0;
2755 :
2756 5794 : const double dfStepSize = bAll ? 0 : 1.0 / (nStepCount - 1);
2757 5794 : constexpr int knIntMax = std::numeric_limits<int>::max();
2758 5794 : int nSampleMax = 0;
2759 5794 : if (bUseGrid)
2760 : {
2761 1165 : if (bAll)
2762 : {
2763 0 : if (nDstYSize > knIntMax / (nDstXSize + 1) - 1)
2764 : {
2765 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps");
2766 0 : return false;
2767 : }
2768 0 : nSampleMax = (nDstXSize + 1) * (nDstYSize + 1);
2769 : }
2770 : else
2771 : {
2772 1165 : if (nStepCount > knIntMax - 2 ||
2773 1165 : (nStepCount + 2) > knIntMax / (nStepCount + 2))
2774 : {
2775 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps : %d",
2776 : nStepCount);
2777 0 : return false;
2778 : }
2779 1165 : nSampleMax = (nStepCount + 2) * (nStepCount + 2);
2780 : }
2781 : }
2782 : else
2783 : {
2784 4629 : if (bAll)
2785 : {
2786 154 : if (nDstXSize > (knIntMax - 2 * nDstYSize) / 2)
2787 : {
2788 : // Extremely unlikely !
2789 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps");
2790 0 : return false;
2791 : }
2792 154 : nSampleMax = 2 * (nDstXSize + nDstYSize);
2793 : }
2794 : else
2795 : {
2796 4475 : if (nStepCount > knIntMax / 4)
2797 : {
2798 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps : %d * 4",
2799 : nStepCount);
2800 0 : return false;
2801 : }
2802 4475 : nSampleMax = nStepCount * 4;
2803 : }
2804 : }
2805 :
2806 : int *pabSuccess =
2807 5794 : static_cast<int *>(VSI_MALLOC2_VERBOSE(sizeof(int), nSampleMax));
2808 : double *padfX = static_cast<double *>(
2809 5794 : VSI_MALLOC2_VERBOSE(sizeof(double) * 3, nSampleMax));
2810 5794 : if (pabSuccess == nullptr || padfX == nullptr)
2811 : {
2812 0 : CPLFree(padfX);
2813 0 : CPLFree(pabSuccess);
2814 0 : return false;
2815 : }
2816 5794 : double *padfY = padfX + nSampleMax;
2817 5794 : double *padfZ = padfX + nSampleMax * 2;
2818 :
2819 : /* -------------------------------------------------------------------- */
2820 : /* Setup sample points on a grid pattern throughout the area. */
2821 : /* -------------------------------------------------------------------- */
2822 5794 : if (bUseGrid)
2823 : {
2824 1165 : if (bAll)
2825 : {
2826 0 : for (int iY = 0; iY <= nDstYSize; ++iY)
2827 : {
2828 0 : for (int iX = 0; iX <= nDstXSize; ++iX)
2829 : {
2830 0 : padfX[nSamplePoints] = nDstXOff + iX;
2831 0 : padfY[nSamplePoints] = nDstYOff + iY;
2832 0 : padfZ[nSamplePoints++] = 0.0;
2833 : }
2834 : }
2835 : }
2836 : else
2837 : {
2838 27960 : for (int iY = 0; iY < nStepCount + 2; iY++)
2839 : {
2840 52425 : const double dfRatioY = (iY == 0) ? 0.5 / nDstXSize
2841 : : (iY <= nStepCount)
2842 25630 : ? (iY - 1) * dfStepSize
2843 1165 : : 1 - 0.5 / nDstXSize;
2844 643080 : for (int iX = 0; iX < nStepCount + 2; iX++)
2845 : {
2846 1205780 : const double dfRatioX = (iX == 0) ? 0.5 / nDstXSize
2847 : : (iX <= nStepCount)
2848 589490 : ? (iX - 1) * dfStepSize
2849 26795 : : 1 - 0.5 / nDstXSize;
2850 616285 : padfX[nSamplePoints] = dfRatioX * nDstXSize + nDstXOff;
2851 616285 : padfY[nSamplePoints] = dfRatioY * nDstYSize + nDstYOff;
2852 616285 : padfZ[nSamplePoints++] = 0.0;
2853 : }
2854 : }
2855 : }
2856 : }
2857 : /* -------------------------------------------------------------------- */
2858 : /* Setup sample points all around the edge of the output raster. */
2859 : /* -------------------------------------------------------------------- */
2860 : else
2861 : {
2862 4629 : if (bAll)
2863 : {
2864 68686 : for (int iX = 0; iX <= nDstXSize; ++iX)
2865 : {
2866 : // Along top
2867 68532 : padfX[nSamplePoints] = nDstXOff + iX;
2868 68532 : padfY[nSamplePoints] = nDstYOff;
2869 68532 : padfZ[nSamplePoints++] = 0.0;
2870 :
2871 : // Along bottom
2872 68532 : padfX[nSamplePoints] = nDstXOff + iX;
2873 68532 : padfY[nSamplePoints] = nDstYOff + nDstYSize;
2874 68532 : padfZ[nSamplePoints++] = 0.0;
2875 : }
2876 :
2877 43748 : for (int iY = 1; iY < nDstYSize; ++iY)
2878 : {
2879 : // Along left
2880 43594 : padfX[nSamplePoints] = nDstXOff;
2881 43594 : padfY[nSamplePoints] = nDstYOff + iY;
2882 43594 : padfZ[nSamplePoints++] = 0.0;
2883 :
2884 : // Along right
2885 43594 : padfX[nSamplePoints] = nDstXOff + nDstXSize;
2886 43594 : padfY[nSamplePoints] = nDstYOff + iY;
2887 43594 : padfZ[nSamplePoints++] = 0.0;
2888 : }
2889 : }
2890 : else
2891 : {
2892 98450 : for (double dfRatio = 0.0; dfRatio <= 1.0 + dfStepSize * 0.5;
2893 93975 : dfRatio += dfStepSize)
2894 : {
2895 : // Along top
2896 93975 : padfX[nSamplePoints] = dfRatio * nDstXSize + nDstXOff;
2897 93975 : padfY[nSamplePoints] = nDstYOff;
2898 93975 : padfZ[nSamplePoints++] = 0.0;
2899 :
2900 : // Along bottom
2901 93975 : padfX[nSamplePoints] = dfRatio * nDstXSize + nDstXOff;
2902 93975 : padfY[nSamplePoints] = nDstYOff + nDstYSize;
2903 93975 : padfZ[nSamplePoints++] = 0.0;
2904 :
2905 : // Along left
2906 93975 : padfX[nSamplePoints] = nDstXOff;
2907 93975 : padfY[nSamplePoints] = dfRatio * nDstYSize + nDstYOff;
2908 93975 : padfZ[nSamplePoints++] = 0.0;
2909 :
2910 : // Along right
2911 93975 : padfX[nSamplePoints] = nDstXSize + nDstXOff;
2912 93975 : padfY[nSamplePoints] = dfRatio * nDstYSize + nDstYOff;
2913 93975 : padfZ[nSamplePoints++] = 0.0;
2914 : }
2915 : }
2916 : }
2917 :
2918 5794 : CPLAssert(nSamplePoints == nSampleMax);
2919 :
2920 : /* -------------------------------------------------------------------- */
2921 : /* Transform them to the input pixel coordinate space */
2922 : /* -------------------------------------------------------------------- */
2923 :
2924 132 : const auto RefreshTransformer = [this]()
2925 : {
2926 44 : if (GDALIsTransformer(psOptions->pTransformerArg,
2927 : GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
2928 : {
2929 0 : GDALRefreshGenImgProjTransformer(psOptions->pTransformerArg);
2930 : }
2931 44 : else if (GDALIsTransformer(psOptions->pTransformerArg,
2932 : GDAL_APPROX_TRANSFORMER_CLASS_NAME))
2933 : {
2934 44 : GDALRefreshApproxTransformer(psOptions->pTransformerArg);
2935 : }
2936 5838 : };
2937 :
2938 5794 : if (bTryWithCheckWithInvertProj)
2939 : {
2940 22 : CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ", "YES");
2941 22 : RefreshTransformer();
2942 : }
2943 5794 : psOptions->pfnTransformer(psOptions->pTransformerArg, TRUE, nSamplePoints,
2944 : padfX, padfY, padfZ, pabSuccess);
2945 5794 : if (bTryWithCheckWithInvertProj)
2946 : {
2947 22 : CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ", nullptr);
2948 22 : RefreshTransformer();
2949 : }
2950 :
2951 : /* -------------------------------------------------------------------- */
2952 : /* Collect the bounds, ignoring any failed points. */
2953 : /* -------------------------------------------------------------------- */
2954 1222230 : for (int i = 0; i < nSamplePoints; i++)
2955 : {
2956 1216440 : if (!pabSuccess[i])
2957 : {
2958 110582 : nFailedCount++;
2959 110582 : continue;
2960 : }
2961 :
2962 : // If this happens this is likely the symptom of a bug somewhere.
2963 1105860 : if (std::isnan(padfX[i]) || std::isnan(padfY[i]))
2964 : {
2965 : static bool bNanCoordFound = false;
2966 0 : if (!bNanCoordFound)
2967 : {
2968 0 : CPLDebug("WARP",
2969 : "ComputeSourceWindow(): "
2970 : "NaN coordinate found on point %d.",
2971 : i);
2972 0 : bNanCoordFound = true;
2973 : }
2974 0 : nFailedCount++;
2975 0 : continue;
2976 : }
2977 :
2978 1105860 : dfMinXOut = std::min(dfMinXOut, padfX[i]);
2979 1105860 : dfMinYOut = std::min(dfMinYOut, padfY[i]);
2980 1105860 : dfMaxXOut = std::max(dfMaxXOut, padfX[i]);
2981 1105860 : dfMaxYOut = std::max(dfMaxYOut, padfY[i]);
2982 : }
2983 :
2984 5794 : CPLFree(padfX);
2985 5794 : CPLFree(pabSuccess);
2986 5794 : return true;
2987 : }
2988 :
2989 : /************************************************************************/
2990 : /* ComputeSourceWindow() */
2991 : /************************************************************************/
2992 :
2993 : /** Given a target window starting at pixel (nDstOff, nDstYOff) and of
2994 : * dimension (nDstXSize, nDstYSize), compute the corresponding window in
2995 : * the source raster, and return the source position in (*pnSrcXOff, *pnSrcYOff),
2996 : * the source dimension in (*pnSrcXSize, *pnSrcYSize).
2997 : * If pdfSrcXExtraSize is not null, its pointed value will be filled with the
2998 : * number of extra source pixels in X dimension to acquire to take into account
2999 : * the size of the resampling kernel. Similarly for pdfSrcYExtraSize for the
3000 : * Y dimension.
3001 : * If pdfSrcFillRatio is not null, its pointed value will be filled with the
3002 : * the ratio of the clamped source raster window size over the unclamped source
3003 : * raster window size. When this ratio is too low, this might be an indication
3004 : * that it might be beneficial to split the target window to avoid requesting
3005 : * too many source pixels.
3006 : */
3007 5483 : CPLErr GDALWarpOperation::ComputeSourceWindow(
3008 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize, int *pnSrcXOff,
3009 : int *pnSrcYOff, int *pnSrcXSize, int *pnSrcYSize, double *pdfSrcXExtraSize,
3010 : double *pdfSrcYExtraSize, double *pdfSrcFillRatio)
3011 :
3012 : {
3013 : /* -------------------------------------------------------------------- */
3014 : /* Figure out whether we just want to do the usual "along the */
3015 : /* edge" sampling, or using a grid. The grid usage is */
3016 : /* important in some weird "inside out" cases like WGS84 to */
3017 : /* polar stereographic around the pole. Also figure out the */
3018 : /* sampling rate. */
3019 : /* -------------------------------------------------------------------- */
3020 5483 : int nStepCount = DEFAULT_STEP_COUNT;
3021 5483 : bool bAll = false;
3022 :
3023 : bool bUseGrid =
3024 5483 : CPLFetchBool(psOptions->papszWarpOptions, "SAMPLE_GRID", false);
3025 :
3026 : const char *pszSampleSteps =
3027 5483 : CSLFetchNameValue(psOptions->papszWarpOptions, "SAMPLE_STEPS");
3028 5483 : if (pszSampleSteps)
3029 : {
3030 94 : if (EQUAL(pszSampleSteps, "ALL"))
3031 : {
3032 94 : bAll = true;
3033 : }
3034 : else
3035 : {
3036 0 : nStepCount = atoi(pszSampleSteps);
3037 0 : nStepCount = std::max(2, nStepCount);
3038 : }
3039 : }
3040 5389 : else if (!bUseGrid)
3041 : {
3042 : // Detect if at least one of the 4 corner in destination raster fails
3043 : // to project back to source.
3044 : // Helps for long-lat to orthographic on areas that are partly in
3045 : // space / partly on Earth. Cf https://github.com/OSGeo/gdal/issues/9056
3046 : double adfCornerX[4];
3047 : double adfCornerY[4];
3048 4535 : double adfCornerZ[4] = {0, 0, 0, 0};
3049 4535 : int anCornerSuccess[4] = {FALSE, FALSE, FALSE, FALSE};
3050 4535 : adfCornerX[0] = nDstXOff;
3051 4535 : adfCornerY[0] = nDstYOff;
3052 4535 : adfCornerX[1] = nDstXOff + nDstXSize;
3053 4535 : adfCornerY[1] = nDstYOff;
3054 4535 : adfCornerX[2] = nDstXOff;
3055 4535 : adfCornerY[2] = nDstYOff + nDstYSize;
3056 4535 : adfCornerX[3] = nDstXOff + nDstXSize;
3057 4535 : adfCornerY[3] = nDstYOff + nDstYSize;
3058 4535 : if (!psOptions->pfnTransformer(psOptions->pTransformerArg, TRUE, 4,
3059 : adfCornerX, adfCornerY, adfCornerZ,
3060 4475 : anCornerSuccess) ||
3061 9010 : !anCornerSuccess[0] || !anCornerSuccess[1] || !anCornerSuccess[2] ||
3062 4475 : !anCornerSuccess[3])
3063 : {
3064 60 : bAll = true;
3065 : }
3066 : }
3067 :
3068 5483 : bool bTryWithCheckWithInvertProj = false;
3069 5483 : double dfMinXOut = std::numeric_limits<double>::infinity();
3070 5483 : double dfMinYOut = std::numeric_limits<double>::infinity();
3071 5483 : double dfMaxXOut = -std::numeric_limits<double>::infinity();
3072 5483 : double dfMaxYOut = -std::numeric_limits<double>::infinity();
3073 :
3074 5483 : int nSamplePoints = 0;
3075 5483 : int nFailedCount = 0;
3076 5483 : if (!ComputeSourceWindowTransformPoints(
3077 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, bUseGrid, bAll,
3078 : nStepCount, bTryWithCheckWithInvertProj, dfMinXOut, dfMinYOut,
3079 : dfMaxXOut, dfMaxYOut, nSamplePoints, nFailedCount))
3080 : {
3081 0 : return CE_Failure;
3082 : }
3083 :
3084 : // Use grid sampling as soon as a special point falls into the extent of
3085 : // the target raster.
3086 5483 : if (!bUseGrid && psOptions->hDstDS)
3087 : {
3088 11322 : for (const auto &xy : aDstXYSpecialPoints)
3089 : {
3090 14653 : if (0 <= xy.first &&
3091 781 : GDALGetRasterXSize(psOptions->hDstDS) >= xy.first &&
3092 8034 : 0 <= xy.second &&
3093 317 : GDALGetRasterYSize(psOptions->hDstDS) >= xy.second)
3094 : {
3095 230 : bUseGrid = true;
3096 230 : bAll = false;
3097 230 : if (!ComputeSourceWindowTransformPoints(
3098 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, bUseGrid,
3099 : bAll, nStepCount, bTryWithCheckWithInvertProj,
3100 : dfMinXOut, dfMinYOut, dfMaxXOut, dfMaxYOut,
3101 : nSamplePoints, nFailedCount))
3102 : {
3103 0 : return CE_Failure;
3104 : }
3105 230 : break;
3106 : }
3107 : }
3108 : }
3109 :
3110 5483 : const int nRasterXSize = GDALGetRasterXSize(psOptions->hSrcDS);
3111 5483 : const int nRasterYSize = GDALGetRasterYSize(psOptions->hSrcDS);
3112 :
3113 : // Try to detect crazy values coming from reprojection that would not
3114 : // have resulted in a PROJ error. Could happen for example with PROJ
3115 : // <= 4.9.2 with inverse UTM/tmerc (Snyder approximation without sanity
3116 : // check) when being far away from the central meridian. But might be worth
3117 : // keeping that even for later versions in case some exotic projection isn't
3118 : // properly sanitized.
3119 5415 : if (nFailedCount == 0 && !bTryWithCheckWithInvertProj &&
3120 5415 : (dfMinXOut < -1e6 || dfMinYOut < -1e6 ||
3121 10902 : dfMaxXOut > nRasterXSize + 1e6 || dfMaxYOut > nRasterYSize + 1e6) &&
3122 22 : !CPLTestBool(CPLGetConfigOption("CHECK_WITH_INVERT_PROJ", "NO")))
3123 : {
3124 22 : CPLDebug("WARP",
3125 : "ComputeSourceWindow(): bogus source dataset window "
3126 : "returned. Trying again with CHECK_WITH_INVERT_PROJ=YES");
3127 22 : bTryWithCheckWithInvertProj = true;
3128 :
3129 : // We should probably perform the coordinate transformation in the
3130 : // warp kernel under CHECK_WITH_INVERT_PROJ too...
3131 22 : if (!ComputeSourceWindowTransformPoints(
3132 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, bUseGrid, bAll,
3133 : nStepCount, bTryWithCheckWithInvertProj, dfMinXOut, dfMinYOut,
3134 : dfMaxXOut, dfMaxYOut, nSamplePoints, nFailedCount))
3135 : {
3136 0 : return CE_Failure;
3137 : }
3138 : }
3139 :
3140 : /* -------------------------------------------------------------------- */
3141 : /* If we got any failures when not using a grid, we should */
3142 : /* really go back and try again with the grid. Sorry for the */
3143 : /* goto. */
3144 : /* -------------------------------------------------------------------- */
3145 5483 : if (!bUseGrid && nFailedCount > 0)
3146 : {
3147 59 : bUseGrid = true;
3148 59 : bAll = false;
3149 59 : if (!ComputeSourceWindowTransformPoints(
3150 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, bUseGrid, bAll,
3151 : nStepCount, bTryWithCheckWithInvertProj, dfMinXOut, dfMinYOut,
3152 : dfMaxXOut, dfMaxYOut, nSamplePoints, nFailedCount))
3153 : {
3154 0 : return CE_Failure;
3155 : }
3156 : }
3157 :
3158 : /* -------------------------------------------------------------------- */
3159 : /* If we get hardly any points (or none) transforming, we give */
3160 : /* up. */
3161 : /* -------------------------------------------------------------------- */
3162 5483 : if (nFailedCount > nSamplePoints - 5)
3163 : {
3164 : const bool bErrorOutIfEmptySourceWindow =
3165 39 : CPLFetchBool(psOptions->papszWarpOptions,
3166 : "ERROR_OUT_IF_EMPTY_SOURCE_WINDOW", true);
3167 39 : if (bErrorOutIfEmptySourceWindow)
3168 : {
3169 3 : CPLError(CE_Failure, CPLE_AppDefined,
3170 : "Too many points (%d out of %d) failed to transform, "
3171 : "unable to compute output bounds.",
3172 : nFailedCount, nSamplePoints);
3173 : }
3174 : else
3175 : {
3176 36 : CPLDebug("WARP", "Cannot determine source window for %d,%d,%d,%d",
3177 : nDstXOff, nDstYOff, nDstXSize, nDstYSize);
3178 : }
3179 39 : return CE_Failure;
3180 : }
3181 :
3182 5444 : if (nFailedCount > 0)
3183 33 : CPLDebug("GDAL",
3184 : "GDALWarpOperation::ComputeSourceWindow() %d out of %d "
3185 : "points failed to transform.",
3186 : nFailedCount, nSamplePoints);
3187 :
3188 : /* -------------------------------------------------------------------- */
3189 : /* In some cases (see https://github.com/OSGeo/gdal/issues/862) */
3190 : /* the reverse transform does not work at some points, so try by */
3191 : /* transforming from source raster space to target raster space and */
3192 : /* see which source coordinates end up being in the AOI in the target */
3193 : /* raster space. */
3194 : /* -------------------------------------------------------------------- */
3195 5444 : if (bUseGrid)
3196 : {
3197 1104 : ComputeSourceWindowStartingFromSource(nDstXOff, nDstYOff, nDstXSize,
3198 : nDstYSize, &dfMinXOut, &dfMinYOut,
3199 : &dfMaxXOut, &dfMaxYOut);
3200 : }
3201 :
3202 : /* -------------------------------------------------------------------- */
3203 : /* Early exit to avoid crazy values to cause a huge nResWinSize that */
3204 : /* would result in a result window wrongly covering the whole raster. */
3205 : /* -------------------------------------------------------------------- */
3206 5444 : if (dfMinXOut > nRasterXSize || dfMaxXOut < 0 || dfMinYOut > nRasterYSize ||
3207 4870 : dfMaxYOut < 0)
3208 : {
3209 956 : *pnSrcXOff = 0;
3210 956 : *pnSrcYOff = 0;
3211 956 : *pnSrcXSize = 0;
3212 956 : *pnSrcYSize = 0;
3213 956 : if (pdfSrcXExtraSize)
3214 956 : *pdfSrcXExtraSize = 0.0;
3215 956 : if (pdfSrcYExtraSize)
3216 956 : *pdfSrcYExtraSize = 0.0;
3217 956 : if (pdfSrcFillRatio)
3218 949 : *pdfSrcFillRatio = 0.0;
3219 956 : return CE_None;
3220 : }
3221 :
3222 : // For scenarios where warping is used as a "decoration", try to clamp
3223 : // source pixel coordinates to integer when very close.
3224 17952 : const auto roundIfCloseEnough = [](double dfVal)
3225 : {
3226 17952 : const double dfRounded = std::round(dfVal);
3227 17952 : if (std::fabs(dfRounded - dfVal) < 1e-6)
3228 14499 : return dfRounded;
3229 3453 : return dfVal;
3230 : };
3231 :
3232 4488 : dfMinXOut = roundIfCloseEnough(dfMinXOut);
3233 4488 : dfMinYOut = roundIfCloseEnough(dfMinYOut);
3234 4488 : dfMaxXOut = roundIfCloseEnough(dfMaxXOut);
3235 4488 : dfMaxYOut = roundIfCloseEnough(dfMaxYOut);
3236 :
3237 4488 : if (m_bIsTranslationOnPixelBoundaries)
3238 : {
3239 311 : CPLAssert(dfMinXOut == std::round(dfMinXOut));
3240 311 : CPLAssert(dfMinYOut == std::round(dfMinYOut));
3241 311 : CPLAssert(dfMaxXOut == std::round(dfMaxXOut));
3242 311 : CPLAssert(dfMaxYOut == std::round(dfMaxYOut));
3243 311 : CPLAssert(std::round(dfMaxXOut - dfMinXOut) == nDstXSize);
3244 311 : CPLAssert(std::round(dfMaxYOut - dfMinYOut) == nDstYSize);
3245 : }
3246 :
3247 : /* -------------------------------------------------------------------- */
3248 : /* How much of a window around our source pixel might we need */
3249 : /* to collect data from based on the resampling kernel? Even */
3250 : /* if the requested central pixel falls off the source image, */
3251 : /* we may need to collect data if some portion of the */
3252 : /* resampling kernel could be on-image. */
3253 : /* -------------------------------------------------------------------- */
3254 4488 : const int nResWinSize = m_bIsTranslationOnPixelBoundaries
3255 4488 : ? 0
3256 4177 : : GWKGetFilterRadius(psOptions->eResampleAlg);
3257 :
3258 : // Take scaling into account.
3259 : // Avoid ridiculous small scaling factors to avoid potential further integer
3260 : // overflows
3261 8976 : const double dfXScale = std::max(1e-3, static_cast<double>(nDstXSize) /
3262 4488 : (dfMaxXOut - dfMinXOut));
3263 8976 : const double dfYScale = std::max(1e-3, static_cast<double>(nDstYSize) /
3264 4488 : (dfMaxYOut - dfMinYOut));
3265 4488 : int nXRadius = dfXScale < 0.95
3266 4488 : ? static_cast<int>(ceil(nResWinSize / dfXScale))
3267 : : nResWinSize;
3268 4488 : int nYRadius = dfYScale < 0.95
3269 4488 : ? static_cast<int>(ceil(nResWinSize / dfYScale))
3270 : : nResWinSize;
3271 :
3272 : /* -------------------------------------------------------------------- */
3273 : /* Allow addition of extra sample pixels to source window to */
3274 : /* avoid missing pixels due to sampling error. In fact, */
3275 : /* fallback to adding a bit to the window if any points failed */
3276 : /* to transform. */
3277 : /* -------------------------------------------------------------------- */
3278 4488 : if (const char *pszSourceExtra =
3279 4488 : CSLFetchNameValue(psOptions->papszWarpOptions, "SOURCE_EXTRA"))
3280 : {
3281 80 : const int nSrcExtra = atoi(pszSourceExtra);
3282 80 : nXRadius += nSrcExtra;
3283 80 : nYRadius += nSrcExtra;
3284 : }
3285 4408 : else if (nFailedCount > 0)
3286 : {
3287 27 : nXRadius += 10;
3288 27 : nYRadius += 10;
3289 : }
3290 :
3291 : /* -------------------------------------------------------------------- */
3292 : /* return bounds. */
3293 : /* -------------------------------------------------------------------- */
3294 : #if DEBUG_VERBOSE
3295 : CPLDebug("WARP",
3296 : "dst=(%d,%d,%d,%d) raw "
3297 : "src=(minx=%.17g,miny=%.17g,maxx=%.17g,maxy=%.17g)",
3298 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, dfMinXOut, dfMinYOut,
3299 : dfMaxXOut, dfMaxYOut);
3300 : #endif
3301 4488 : const int nMinXOutClamped = static_cast<int>(std::max(0.0, dfMinXOut));
3302 4488 : const int nMinYOutClamped = static_cast<int>(std::max(0.0, dfMinYOut));
3303 : const int nMaxXOutClamped = static_cast<int>(
3304 4488 : std::min(ceil(dfMaxXOut), static_cast<double>(nRasterXSize)));
3305 : const int nMaxYOutClamped = static_cast<int>(
3306 4488 : std::min(ceil(dfMaxYOut), static_cast<double>(nRasterYSize)));
3307 :
3308 : const double dfSrcXSizeRaw = std::max(
3309 13464 : 0.0, std::min(static_cast<double>(nRasterXSize - nMinXOutClamped),
3310 4488 : dfMaxXOut - dfMinXOut));
3311 : const double dfSrcYSizeRaw = std::max(
3312 13464 : 0.0, std::min(static_cast<double>(nRasterYSize - nMinYOutClamped),
3313 4488 : dfMaxYOut - dfMinYOut));
3314 :
3315 : // If we cover more than 90% of the width, then use it fully (helps for
3316 : // anti-meridian discontinuities)
3317 4488 : if (nMaxXOutClamped - nMinXOutClamped > 0.9 * nRasterXSize)
3318 : {
3319 1528 : *pnSrcXOff = 0;
3320 1528 : *pnSrcXSize = nRasterXSize;
3321 : }
3322 : else
3323 : {
3324 2960 : *pnSrcXOff =
3325 2960 : std::max(0, std::min(nMinXOutClamped - nXRadius, nRasterXSize));
3326 2960 : *pnSrcXSize =
3327 8880 : std::max(0, std::min(nRasterXSize - *pnSrcXOff,
3328 2960 : nMaxXOutClamped - *pnSrcXOff + nXRadius));
3329 : }
3330 :
3331 4488 : if (nMaxYOutClamped - nMinYOutClamped > 0.9 * nRasterYSize)
3332 : {
3333 1412 : *pnSrcYOff = 0;
3334 1412 : *pnSrcYSize = nRasterYSize;
3335 : }
3336 : else
3337 : {
3338 3076 : *pnSrcYOff =
3339 3076 : std::max(0, std::min(nMinYOutClamped - nYRadius, nRasterYSize));
3340 3076 : *pnSrcYSize =
3341 9228 : std::max(0, std::min(nRasterYSize - *pnSrcYOff,
3342 3076 : nMaxYOutClamped - *pnSrcYOff + nYRadius));
3343 : }
3344 :
3345 4488 : if (pdfSrcXExtraSize)
3346 4488 : *pdfSrcXExtraSize = *pnSrcXSize - dfSrcXSizeRaw;
3347 4488 : if (pdfSrcYExtraSize)
3348 4488 : *pdfSrcYExtraSize = *pnSrcYSize - dfSrcYSizeRaw;
3349 :
3350 : // Computed the ratio of the clamped source raster window size over
3351 : // the unclamped source raster window size.
3352 4488 : if (pdfSrcFillRatio)
3353 3615 : *pdfSrcFillRatio =
3354 7230 : static_cast<double>(*pnSrcXSize) * (*pnSrcYSize) /
3355 3615 : std::max(1.0, (dfMaxXOut - dfMinXOut + 2 * nXRadius) *
3356 3615 : (dfMaxYOut - dfMinYOut + 2 * nYRadius));
3357 :
3358 4488 : return CE_None;
3359 : }
3360 :
3361 : /************************************************************************/
3362 : /* ReportTiming() */
3363 : /************************************************************************/
3364 :
3365 10917 : void GDALWarpOperation::ReportTiming(const char *pszMessage)
3366 :
3367 : {
3368 10917 : if (!bReportTimings)
3369 10917 : return;
3370 :
3371 0 : const unsigned long nNewTime = VSITime(nullptr);
3372 :
3373 0 : if (pszMessage != nullptr)
3374 : {
3375 0 : CPLDebug("WARP_TIMING", "%s: %lds", pszMessage,
3376 0 : static_cast<long>(nNewTime - nLastTimeReported));
3377 : }
3378 :
3379 0 : nLastTimeReported = nNewTime;
3380 : }
|
