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