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