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 1425 : GDALWarpOperation::GDALWarpOperation()
156 : : psOptions(nullptr), hIOMutex(nullptr), hWarpMutex(nullptr),
157 : nChunkListCount(0), nChunkListMax(0), pasChunkList(nullptr),
158 1425 : bReportTimings(FALSE), nLastTimeReported(0), psThreadData(nullptr)
159 : {
160 1425 : }
161 :
162 : /************************************************************************/
163 : /* ~GDALWarpOperation() */
164 : /************************************************************************/
165 :
166 1922 : GDALWarpOperation::~GDALWarpOperation()
167 :
168 : {
169 : {
170 2850 : std::lock_guard<std::mutex> oLock(gMutex);
171 1425 : auto oItem = gMapPrivate.find(this);
172 1425 : if (oItem != gMapPrivate.end())
173 : {
174 142 : gMapPrivate.erase(oItem);
175 : }
176 : }
177 :
178 1425 : WipeOptions();
179 :
180 1425 : if (hIOMutex != nullptr)
181 : {
182 5 : CPLDestroyMutex(hIOMutex);
183 5 : CPLDestroyMutex(hWarpMutex);
184 : }
185 :
186 1425 : WipeChunkList();
187 1425 : if (psThreadData)
188 1421 : GWKThreadsEnd(psThreadData);
189 1922 : }
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 1429 : void GDALWarpOperation::WipeOptions()
207 :
208 : {
209 1429 : if (psOptions != nullptr)
210 : {
211 1425 : GDALDestroyWarpOptions(psOptions);
212 1425 : psOptions = nullptr;
213 : }
214 1429 : }
215 :
216 : /************************************************************************/
217 : /* ValidateOptions() */
218 : /************************************************************************/
219 :
220 1425 : int GDALWarpOperation::ValidateOptions()
221 :
222 : {
223 1425 : 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 1425 : 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 1425 : if (psOptions->eResampleAlg != GRA_NearestNeighbour &&
241 774 : 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 1425 : if (static_cast<int>(psOptions->eWorkingDataType) < 1 ||
262 1425 : 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 1674 : 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 1425 : 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 1425 : 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 1425 : 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 1425 : 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 3316 : for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
319 : {
320 3782 : if (psOptions->panSrcBands[iBand] < 1 ||
321 1891 : psOptions->panSrcBands[iBand] >
322 1891 : 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 3769 : if (psOptions->hDstDS != nullptr &&
330 1878 : (psOptions->panDstBands[iBand] < 1 ||
331 1878 : psOptions->panDstBands[iBand] >
332 1878 : 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 3769 : if (psOptions->hDstDS != nullptr &&
341 1878 : GDALGetRasterAccess(GDALGetRasterBand(
342 1878 : 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 1425 : 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 1425 : 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 1425 : 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 1425 : CSLFetchNameValue(psOptions->papszWarpOptions, "SAMPLE_STEPS");
378 1425 : 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 1425 : 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 1425 : 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 1425 : 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 1425 : 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 2850 : const bool bErrorOutIfEmptySourceWindow = CPLFetchBool(
432 1425 : psOptions->papszWarpOptions, "ERROR_OUT_IF_EMPTY_SOURCE_WINDOW", true);
433 1773 : 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 1425 : 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 1425 : static void SetTieStrategy(GDALWarpOptions *psOptions, CPLErr *peErr)
481 : {
482 1425 : if (const char *pszTieStrategy =
483 1425 : 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 1425 : }
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 1425 : CPLErr GDALWarpOperation::Initialize(const GDALWarpOptions *psNewOptions)
530 :
531 : {
532 : /* -------------------------------------------------------------------- */
533 : /* Copy the passed in options. */
534 : /* -------------------------------------------------------------------- */
535 1425 : if (psOptions != nullptr)
536 0 : WipeOptions();
537 :
538 1425 : CPLErr eErr = CE_None;
539 :
540 1425 : psOptions = GDALCloneWarpOptions(psNewOptions);
541 2850 : psOptions->papszWarpOptions =
542 1425 : 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 1425 : 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 1425 : GDALWarpResolveWorkingDataType(psOptions);
558 1425 : 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 1425 : if (psOptions->dfWarpMemoryLimit == 0.0)
568 : {
569 1170 : psOptions->dfWarpMemoryLimit = 64.0 * 1024 * 1024;
570 : }
571 :
572 : /* -------------------------------------------------------------------- */
573 : /* Are we doing timings? */
574 : /* -------------------------------------------------------------------- */
575 1425 : bReportTimings =
576 1425 : CPLFetchBool(psOptions->papszWarpOptions, "REPORT_TIMINGS", false);
577 :
578 : /* -------------------------------------------------------------------- */
579 : /* Support creating cutline from text warpoption. */
580 : /* -------------------------------------------------------------------- */
581 : const char *pszCutlineWKT =
582 1425 : CSLFetchNameValue(psOptions->papszWarpOptions, "CUTLINE");
583 :
584 1425 : 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 1425 : CSLFetchNameValue(psOptions->papszWarpOptions, "CUTLINE_BLEND_DIST");
598 1425 : if (pszBD)
599 0 : psOptions->dfCutlineBlendDist = CPLAtof(pszBD);
600 :
601 : /* -------------------------------------------------------------------- */
602 : /* Set SRC_ALPHA_MAX if not provided. */
603 : /* -------------------------------------------------------------------- */
604 1425 : if (psOptions->hSrcDS != nullptr && psOptions->nSrcAlphaBand > 0 &&
605 2936 : 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 1412 : if (psOptions->hDstDS != nullptr && psOptions->nDstAlphaBand > 0 &&
618 3068 : 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 1425 : if (!ValidateOptions())
631 0 : eErr = CE_Failure;
632 :
633 1425 : if (eErr != CE_None)
634 : {
635 4 : WipeOptions();
636 : }
637 : else
638 : {
639 2842 : psThreadData = GWKThreadsCreate(psOptions->papszWarpOptions,
640 1421 : psOptions->pfnTransformer,
641 1421 : psOptions->pTransformerArg);
642 1421 : 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 4263 : for (double dfY : {-89.9999, 89.9999})
654 : {
655 2842 : double dfX = 0;
656 2842 : if ((GDALIsTransformer(psOptions->pTransformerArg,
657 2130 : GDAL_APPROX_TRANSFORMER_CLASS_NAME) &&
658 2130 : GDALTransformLonLatToDestApproxTransformer(
659 5684 : psOptions->pTransformerArg, &dfX, &dfY)) ||
660 1696 : (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 1421 : m_bIsTranslationOnPixelBoundaries =
671 2842 : GDALTransformIsTranslationOnPixelBoundaries(
672 1600 : psOptions->pfnTransformer, psOptions->pTransformerArg) &&
673 179 : CPLTestBool(
674 : CPLGetConfigOption("GDAL_WARP_USE_TRANSLATION_OPTIM", "YES"));
675 1421 : if (m_bIsTranslationOnPixelBoundaries)
676 : {
677 173 : CPLDebug("WARP",
678 : "Using translation-on-pixel-boundaries optimization");
679 : }
680 : }
681 :
682 1425 : 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 2323 : 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 2323 : const int nWordSize = GDALGetDataTypeSizeBytes(psOptions->eWorkingDataType);
709 :
710 2323 : void *pDstBuffer = VSI_MALLOC3_VERBOSE(
711 : cpl::fits_on<int>(nWordSize * psOptions->nBandCount), nDstXSize,
712 : nDstYSize);
713 2323 : if (pDstBuffer)
714 : {
715 2323 : InitializeDestinationBuffer(pDstBuffer, nDstXSize, nDstYSize,
716 : pbInitialized);
717 : }
718 2323 : 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 2394 : void GDALWarpOperation::InitializeDestinationBuffer(void *pDstBuffer,
741 : int nDstXSize,
742 : int nDstYSize,
743 : int *pbInitialized)
744 : {
745 2394 : const int nWordSize = GDALGetDataTypeSizeBytes(psOptions->eWorkingDataType);
746 :
747 2394 : const GPtrDiff_t nBandSize =
748 2394 : static_cast<GPtrDiff_t>(nWordSize) * nDstXSize * nDstYSize;
749 :
750 : /* -------------------------------------------------------------------- */
751 : /* Initialize if requested in the options */
752 : /* -------------------------------------------------------------------- */
753 : const char *pszInitDest =
754 2394 : CSLFetchNameValue(psOptions->papszWarpOptions, "INIT_DEST");
755 :
756 2394 : if (pszInitDest == nullptr || EQUAL(pszInitDest, ""))
757 : {
758 400 : if (pbInitialized != nullptr)
759 : {
760 400 : *pbInitialized = FALSE;
761 : }
762 400 : return;
763 : }
764 :
765 1994 : if (pbInitialized != nullptr)
766 : {
767 1591 : *pbInitialized = TRUE;
768 : }
769 :
770 : char **papszInitValues =
771 1994 : CSLTokenizeStringComplex(pszInitDest, ",", FALSE, FALSE);
772 1994 : const int nInitCount = CSLCount(papszInitValues);
773 :
774 4585 : for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
775 : {
776 2591 : double adfInitRealImag[2] = {0.0, 0.0};
777 : const char *pszBandInit =
778 2591 : papszInitValues[std::min(iBand, nInitCount - 1)];
779 :
780 2591 : if (EQUAL(pszBandInit, "NO_DATA") &&
781 467 : psOptions->padfDstNoDataReal != nullptr)
782 : {
783 443 : adfInitRealImag[0] = psOptions->padfDstNoDataReal[iBand];
784 443 : if (psOptions->padfDstNoDataImag != nullptr)
785 : {
786 369 : adfInitRealImag[1] = psOptions->padfDstNoDataImag[iBand];
787 : }
788 : }
789 : else
790 : {
791 2148 : CPLStringToComplex(pszBandInit, adfInitRealImag + 0,
792 : adfInitRealImag + 1);
793 : }
794 :
795 2591 : GByte *pBandData = static_cast<GByte *>(pDstBuffer) + iBand * nBandSize;
796 :
797 2591 : if (psOptions->eWorkingDataType == GDT_Byte)
798 : {
799 4290 : memset(pBandData,
800 : std::max(
801 2145 : 0, std::min(255, static_cast<int>(adfInitRealImag[0]))),
802 : nBandSize);
803 : }
804 889 : else if (!std::isnan(adfInitRealImag[0]) && adfInitRealImag[0] == 0.0 &&
805 889 : !std::isnan(adfInitRealImag[1]) && adfInitRealImag[1] == 0.0)
806 : {
807 380 : memset(pBandData, 0, nBandSize);
808 : }
809 66 : else if (!std::isnan(adfInitRealImag[1]) && adfInitRealImag[1] == 0.0)
810 : {
811 66 : GDALCopyWords64(&adfInitRealImag, GDT_Float64, 0, pBandData,
812 66 : psOptions->eWorkingDataType, nWordSize,
813 66 : 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 1994 : 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 2323 : void GDALWarpOperation::DestroyDestinationBuffer(void *pDstBuffer)
836 : {
837 2323 : VSIFree(pDstBuffer);
838 2323 : }
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 1074 : 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 1074 : WipeChunkList();
886 1074 : CollectChunkListInternal(nDstXOff, nDstYOff, nDstXSize, nDstYSize);
887 :
888 : // Sort chunks from top to bottom, and for equal y, from left to right.
889 1074 : 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 1074 : const int knIntMax = std::numeric_limits<int>::max();
907 1074 : const int knIntMin = std::numeric_limits<int>::min();
908 1074 : int nSrcXOff = knIntMax;
909 1074 : int nSrcYOff = knIntMax;
910 1074 : int nSrcX2Off = knIntMin;
911 1074 : int nSrcY2Off = knIntMin;
912 1074 : double dfApproxAccArea = 0;
913 3065 : for (int iChunk = 0; pasChunkList != nullptr && iChunk < nChunkListCount;
914 : iChunk++)
915 : {
916 1991 : GDALWarpChunk *pasThisChunk = pasChunkList + iChunk;
917 1991 : nSrcXOff = std::min(nSrcXOff, pasThisChunk->sx);
918 1991 : nSrcYOff = std::min(nSrcYOff, pasThisChunk->sy);
919 1991 : nSrcX2Off = std::max(nSrcX2Off, pasThisChunk->sx + pasThisChunk->ssx);
920 1991 : nSrcY2Off = std::max(nSrcY2Off, pasThisChunk->sy + pasThisChunk->ssy);
921 1991 : dfApproxAccArea +=
922 1991 : static_cast<double>(pasThisChunk->ssx) * pasThisChunk->ssy;
923 : }
924 1074 : if (nSrcXOff < nSrcX2Off)
925 : {
926 1068 : const double dfTotalArea =
927 1068 : static_cast<double>(nSrcX2Off - nSrcXOff) * (nSrcY2Off - nSrcYOff);
928 : // This is really a gross heuristics, but should work in most cases
929 1068 : if (dfApproxAccArea >= dfTotalArea * 0.80)
930 : {
931 1068 : GDALDataset::FromHandle(psOptions->hSrcDS)
932 1068 : ->AdviseRead(nSrcXOff, nSrcYOff, nSrcX2Off - nSrcXOff,
933 : nSrcY2Off - nSrcYOff, nDstXSize, nDstYSize,
934 1068 : psOptions->eWorkingDataType, psOptions->nBandCount,
935 1068 : psOptions->panSrcBands, nullptr);
936 : }
937 : }
938 1074 : }
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 1069 : 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 1069 : CollectChunkList(nDstXOff, nDstYOff, nDstXSize, nDstYSize);
975 :
976 : /* -------------------------------------------------------------------- */
977 : /* Total up output pixels to process. */
978 : /* -------------------------------------------------------------------- */
979 1069 : double dfTotalPixels = 0.0;
980 :
981 3050 : for (int iChunk = 0; pasChunkList != nullptr && iChunk < nChunkListCount;
982 : iChunk++)
983 : {
984 1981 : GDALWarpChunk *pasThisChunk = pasChunkList + iChunk;
985 1981 : const double dfChunkPixels =
986 1981 : pasThisChunk->dsx * static_cast<double>(pasThisChunk->dsy);
987 :
988 1981 : dfTotalPixels += dfChunkPixels;
989 : }
990 :
991 : /* -------------------------------------------------------------------- */
992 : /* Process them one at a time, updating the progress */
993 : /* information for each region. */
994 : /* -------------------------------------------------------------------- */
995 1069 : double dfPixelsProcessed = 0.0;
996 :
997 3045 : for (int iChunk = 0; pasChunkList != nullptr && iChunk < nChunkListCount;
998 : iChunk++)
999 : {
1000 1981 : GDALWarpChunk *pasThisChunk = pasChunkList + iChunk;
1001 1981 : const double dfChunkPixels =
1002 1981 : pasThisChunk->dsx * static_cast<double>(pasThisChunk->dsy);
1003 :
1004 1981 : const double dfProgressBase = dfPixelsProcessed / dfTotalPixels;
1005 1981 : const double dfProgressScale = dfChunkPixels / dfTotalPixels;
1006 :
1007 1981 : 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 1981 : if (eErr != CE_None)
1014 5 : return eErr;
1015 :
1016 1976 : dfPixelsProcessed += dfChunkPixels;
1017 : }
1018 :
1019 1064 : WipeChunkList();
1020 :
1021 1064 : psOptions->pfnProgress(1.0, "", psOptions->pProgressArg);
1022 :
1023 1064 : 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 3568 : void GDALWarpOperation::WipeChunkList()
1290 :
1291 : {
1292 3568 : CPLFree(pasChunkList);
1293 3568 : pasChunkList = nullptr;
1294 3568 : nChunkListCount = 0;
1295 3568 : nChunkListMax = 0;
1296 3568 : }
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 3491 : 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 3491 : int nSrcPixelCostInBits = GDALGetDataTypeSize(psOptions->eWorkingDataType) *
1316 3491 : psOptions->nBandCount;
1317 :
1318 3491 : if (psOptions->pfnSrcDensityMaskFunc != nullptr)
1319 0 : nSrcPixelCostInBits += 32; // Float mask?
1320 :
1321 3491 : GDALRasterBandH hSrcBand = nullptr;
1322 3491 : if (psOptions->nBandCount > 0)
1323 : hSrcBand =
1324 3491 : GDALGetRasterBand(psOptions->hSrcDS, psOptions->panSrcBands[0]);
1325 :
1326 3491 : if (psOptions->nSrcAlphaBand > 0 || psOptions->hCutline != nullptr)
1327 125 : nSrcPixelCostInBits += 32; // UnifiedSrcDensity float mask.
1328 6732 : else if (hSrcBand != nullptr &&
1329 3366 : (GDALGetMaskFlags(hSrcBand) & GMF_PER_DATASET))
1330 5 : nSrcPixelCostInBits += 1; // UnifiedSrcValid bit mask.
1331 :
1332 3491 : if (psOptions->papfnSrcPerBandValidityMaskFunc != nullptr ||
1333 3491 : psOptions->padfSrcNoDataReal != nullptr)
1334 156 : nSrcPixelCostInBits += psOptions->nBandCount; // Bit/band mask.
1335 :
1336 3491 : if (psOptions->pfnSrcValidityMaskFunc != nullptr)
1337 0 : nSrcPixelCostInBits += 1; // Bit mask.
1338 :
1339 : /* -------------------------------------------------------------------- */
1340 : /* What about the cost for the destination. */
1341 : /* -------------------------------------------------------------------- */
1342 3491 : int nDstPixelCostInBits = GDALGetDataTypeSize(psOptions->eWorkingDataType) *
1343 3491 : psOptions->nBandCount;
1344 :
1345 3491 : if (psOptions->pfnDstDensityMaskFunc != nullptr)
1346 0 : nDstPixelCostInBits += 32;
1347 :
1348 3491 : if (psOptions->padfDstNoDataReal != nullptr ||
1349 2260 : psOptions->pfnDstValidityMaskFunc != nullptr)
1350 1231 : nDstPixelCostInBits += psOptions->nBandCount;
1351 :
1352 3491 : if (psOptions->nDstAlphaBand > 0)
1353 269 : nDstPixelCostInBits += 32; // DstDensity float mask.
1354 :
1355 3491 : const double dfTotalMemoryUse =
1356 3491 : (static_cast<double>(nSrcPixelCostInBits) * nSrcXSize * nSrcYSize +
1357 3491 : static_cast<double>(nDstPixelCostInBits) * nDstXSize * nDstYSize) /
1358 : 8.0;
1359 3491 : return dfTotalMemoryUse;
1360 : }
1361 :
1362 : /************************************************************************/
1363 : /* CollectChunkListInternal() */
1364 : /************************************************************************/
1365 :
1366 3732 : 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 3732 : int nSrcXOff = 0;
1375 3732 : int nSrcYOff = 0;
1376 3732 : int nSrcXSize = 0;
1377 3732 : int nSrcYSize = 0;
1378 3732 : double dfSrcXExtraSize = 0.0;
1379 3732 : double dfSrcYExtraSize = 0.0;
1380 3732 : double dfSrcFillRatio = 0.0;
1381 : CPLErr eErr =
1382 3732 : ComputeSourceWindow(nDstXOff, nDstYOff, nDstXSize, nDstYSize, &nSrcXOff,
1383 : &nSrcYOff, &nSrcXSize, &nSrcYSize, &dfSrcXExtraSize,
1384 : &dfSrcYExtraSize, &dfSrcFillRatio);
1385 :
1386 3732 : 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 4612 : 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 3320 : 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 6640 : (nDstXSize > 2 || nDstYSize > 2)) ||
1427 2575 : (dfSrcFillRatio > 0 && dfSrcFillRatio < 0.5 &&
1428 284 : (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 1991 : if (nChunkListCount == nChunkListMax)
1524 : {
1525 1215 : nChunkListMax = nChunkListMax * 2 + 1;
1526 1215 : pasChunkList = static_cast<GDALWarpChunk *>(
1527 1215 : CPLRealloc(pasChunkList, sizeof(GDALWarpChunk) * nChunkListMax));
1528 : }
1529 :
1530 1991 : pasChunkList[nChunkListCount].dx = nDstXOff;
1531 1991 : pasChunkList[nChunkListCount].dy = nDstYOff;
1532 1991 : pasChunkList[nChunkListCount].dsx = nDstXSize;
1533 1991 : pasChunkList[nChunkListCount].dsy = nDstYSize;
1534 1991 : pasChunkList[nChunkListCount].sx = nSrcXOff;
1535 1991 : pasChunkList[nChunkListCount].sy = nSrcYOff;
1536 1991 : pasChunkList[nChunkListCount].ssx = nSrcXSize;
1537 1991 : pasChunkList[nChunkListCount].ssy = nSrcYSize;
1538 1991 : pasChunkList[nChunkListCount].sExtraSx = dfSrcXExtraSize;
1539 1991 : pasChunkList[nChunkListCount].sExtraSy = dfSrcYExtraSize;
1540 :
1541 1991 : nChunkListCount++;
1542 :
1543 1991 : 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 1991 : 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 1991 : ReportTiming(nullptr);
1630 :
1631 : /* -------------------------------------------------------------------- */
1632 : /* Allocate the output buffer. */
1633 : /* -------------------------------------------------------------------- */
1634 1991 : int bDstBufferInitialized = FALSE;
1635 : void *pDstBuffer =
1636 1991 : CreateDestinationBuffer(nDstXSize, nDstYSize, &bDstBufferInitialized);
1637 1991 : 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 1991 : GDALDataset *poDstDS = GDALDataset::FromHandle(psOptions->hDstDS);
1647 1991 : 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 3982 : CPLErr eErr = WarpRegionToBuffer(
1681 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, pDstBuffer,
1682 1991 : 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 1991 : if (eErr == CE_None)
1689 : {
1690 1986 : if (psOptions->nBandCount == 1)
1691 : {
1692 : // Particular case to simplify the stack a bit.
1693 1853 : eErr = poDstDS->GetRasterBand(psOptions->panDstBands[0])
1694 3706 : ->RasterIO(GF_Write, nDstXOff, nDstYOff, nDstXSize,
1695 : nDstYSize, pDstBuffer, nDstXSize, nDstYSize,
1696 1853 : 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 3972 : if (eErr == CE_None &&
1708 1986 : 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 1986 : ReportTiming("Output buffer write");
1719 : }
1720 :
1721 : /* -------------------------------------------------------------------- */
1722 : /* Cleanup and return. */
1723 : /* -------------------------------------------------------------------- */
1724 1991 : DestroyDestinationBuffer(pDstBuffer);
1725 :
1726 1991 : 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 2394 : 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 2394 : const int nWordSize = GDALGetDataTypeSizeBytes(psOptions->eWorkingDataType);
1833 :
1834 2394 : CPLAssert(eBufDataType == psOptions->eWorkingDataType);
1835 :
1836 : /* -------------------------------------------------------------------- */
1837 : /* If not given a corresponding source window compute one now. */
1838 : /* -------------------------------------------------------------------- */
1839 2394 : 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 4716 : GDALWarpKernel oWK;
1871 :
1872 2358 : oWK.eResample = m_bIsTranslationOnPixelBoundaries ? GRA_NearestNeighbour
1873 2181 : : psOptions->eResampleAlg;
1874 2358 : oWK.eTieStrategy = psOptions->eTieStrategy;
1875 2358 : oWK.nBands = psOptions->nBandCount;
1876 2358 : oWK.eWorkingDataType = psOptions->eWorkingDataType;
1877 :
1878 2358 : oWK.pfnTransformer = psOptions->pfnTransformer;
1879 2358 : oWK.pTransformerArg = psOptions->pTransformerArg;
1880 :
1881 2358 : oWK.pfnProgress = psOptions->pfnProgress;
1882 2358 : oWK.pProgress = psOptions->pProgressArg;
1883 2358 : oWK.dfProgressBase = dfProgressBase;
1884 2358 : oWK.dfProgressScale = dfProgressScale;
1885 :
1886 2358 : oWK.papszWarpOptions = psOptions->papszWarpOptions;
1887 2358 : oWK.psThreadData = psThreadData;
1888 :
1889 2358 : 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 2358 : oWK.nSrcXOff = nSrcXOff;
1898 2358 : oWK.nSrcYOff = nSrcYOff;
1899 2358 : oWK.nSrcXSize = nSrcXSize;
1900 2358 : oWK.nSrcYSize = nSrcYSize;
1901 2358 : oWK.dfSrcXExtraSize = dfSrcXExtraSize;
1902 2358 : oWK.dfSrcYExtraSize = dfSrcYExtraSize;
1903 :
1904 2358 : GInt64 nAlloc64 =
1905 2358 : nWordSize *
1906 2358 : (static_cast<GInt64>(nSrcXSize) * nSrcYSize + WARP_EXTRA_ELTS) *
1907 2358 : 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 2358 : oWK.papabySrcImage = static_cast<GByte **>(
1919 2358 : CPLCalloc(sizeof(GByte *), psOptions->nBandCount));
1920 4716 : oWK.papabySrcImage[0] =
1921 2358 : static_cast<GByte *>(VSI_MALLOC_VERBOSE(static_cast<size_t>(nAlloc64)));
1922 :
1923 2358 : CPLErr eErr =
1924 2055 : nSrcXSize != 0 && nSrcYSize != 0 && oWK.papabySrcImage[0] == nullptr
1925 4413 : ? CE_Failure
1926 : : CE_None;
1927 :
1928 5310 : for (int i = 0; i < psOptions->nBandCount && eErr == CE_None; i++)
1929 2952 : oWK.papabySrcImage[i] =
1930 2952 : reinterpret_cast<GByte *>(oWK.papabySrcImage[0]) +
1931 2952 : nWordSize *
1932 2952 : (static_cast<GPtrDiff_t>(nSrcXSize) * nSrcYSize +
1933 2952 : WARP_EXTRA_ELTS) *
1934 2952 : i;
1935 :
1936 2358 : if (eErr == CE_None && nSrcXSize > 0 && nSrcYSize > 0)
1937 : {
1938 2055 : GDALDataset *poSrcDS = GDALDataset::FromHandle(psOptions->hSrcDS);
1939 2055 : if (psOptions->nBandCount == 1)
1940 : {
1941 : // Particular case to simplify the stack a bit.
1942 1769 : eErr = poSrcDS->GetRasterBand(psOptions->panSrcBands[0])
1943 3538 : ->RasterIO(GF_Read, nSrcXOff, nSrcYOff, nSrcXSize,
1944 1769 : nSrcYSize, oWK.papabySrcImage[0], nSrcXSize,
1945 1769 : 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 2358 : ReportTiming("Input buffer read");
1962 :
1963 : /* -------------------------------------------------------------------- */
1964 : /* Initialize destination buffer. */
1965 : /* -------------------------------------------------------------------- */
1966 2358 : oWK.nDstXOff = nDstXOff;
1967 2358 : oWK.nDstYOff = nDstYOff;
1968 2358 : oWK.nDstXSize = nDstXSize;
1969 2358 : oWK.nDstYSize = nDstYSize;
1970 :
1971 2358 : oWK.papabyDstImage = reinterpret_cast<GByte **>(
1972 2358 : CPLCalloc(sizeof(GByte *), psOptions->nBandCount));
1973 :
1974 5308 : for (int i = 0; i < psOptions->nBandCount && eErr == CE_None; i++)
1975 : {
1976 2950 : oWK.papabyDstImage[i] =
1977 2950 : static_cast<GByte *>(pDataBuf) +
1978 2950 : 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 2358 : 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 2358 : 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 2358 : 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 2358 : if (eErr == CE_None && psOptions->padfSrcNoDataReal != nullptr &&
2077 122 : nSrcXSize > 0 && nSrcYSize > 0)
2078 : {
2079 122 : CPLAssert(oWK.papanBandSrcValid == nullptr);
2080 :
2081 122 : bool bAllBandsAllValid = true;
2082 290 : for (int i = 0; i < psOptions->nBandCount && eErr == CE_None; i++)
2083 : {
2084 168 : eErr = CreateKernelMask(&oWK, i, "BandSrcValid");
2085 168 : if (eErr == CE_None)
2086 : {
2087 168 : double adfNoData[2] = {psOptions->padfSrcNoDataReal[i],
2088 168 : psOptions->padfSrcNoDataImag != nullptr
2089 168 : ? psOptions->padfSrcNoDataImag[i]
2090 168 : : 0.0};
2091 :
2092 168 : int bAllValid = FALSE;
2093 336 : eErr = GDALWarpNoDataMasker(
2094 168 : adfNoData, 1, psOptions->eWorkingDataType, oWK.nSrcXOff,
2095 : oWK.nSrcYOff, oWK.nSrcXSize, oWK.nSrcYSize,
2096 168 : &(oWK.papabySrcImage[i]), FALSE, oWK.papanBandSrcValid[i],
2097 : &bAllValid);
2098 168 : 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 122 : 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 109 : for (int k = 0; k < psOptions->nBandCount; k++)
2113 59 : CPLFree(oWK.papanBandSrcValid[k]);
2114 50 : CPLFree(oWK.papanBandSrcValid);
2115 50 : 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 122 : 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 69 : 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 2358 : psOptions->nBandCount < 1
2200 2358 : ? nullptr
2201 2358 : : GDALGetRasterBand(psOptions->hSrcDS, psOptions->panSrcBands[0]);
2202 :
2203 2356 : if (eErr == CE_None && oWK.pafUnifiedSrcDensity == nullptr &&
2204 2276 : oWK.panUnifiedSrcValid == nullptr && psOptions->nSrcAlphaBand <= 0 &&
2205 2181 : (GDALGetMaskFlags(hSrcBand) & GMF_PER_DATASET)
2206 : // Need to double check for -nosrcalpha case.
2207 4716 : && !(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 2358 : if (eErr == CE_None && psOptions->padfDstNoDataReal != nullptr)
2230 : {
2231 627 : CPLAssert(oWK.panDstValid == nullptr);
2232 :
2233 627 : const GPtrDiff_t nMaskBits =
2234 627 : static_cast<GPtrDiff_t>(oWK.nDstXSize) * oWK.nDstYSize;
2235 :
2236 627 : eErr = CreateKernelMask(&oWK, 0 /* not used */, "DstValid");
2237 : GUInt32 *panBandMask =
2238 627 : eErr == CE_None ? CPLMaskCreate(nMaskBits, true) : nullptr;
2239 :
2240 627 : if (eErr == CE_None && panBandMask != nullptr)
2241 : {
2242 1200 : for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
2243 : {
2244 671 : CPLMaskSetAll(panBandMask, nMaskBits);
2245 :
2246 671 : double adfNoData[2] = {psOptions->padfDstNoDataReal[iBand],
2247 671 : psOptions->padfDstNoDataImag != nullptr
2248 671 : ? psOptions->padfDstNoDataImag[iBand]
2249 671 : : 0.0};
2250 :
2251 671 : int bAllValid = FALSE;
2252 1342 : eErr = GDALWarpNoDataMasker(
2253 671 : adfNoData, 1, psOptions->eWorkingDataType, oWK.nDstXOff,
2254 : oWK.nDstYOff, oWK.nDstXSize, oWK.nDstYSize,
2255 671 : 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 671 : 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 573 : CPLMaskMerge(oWK.panDstValid, panBandMask, nMaskBits);
2271 : }
2272 627 : CPLFree(panBandMask);
2273 : }
2274 : }
2275 :
2276 : /* -------------------------------------------------------------------- */
2277 : /* Release IO Mutex, and acquire warper mutex. */
2278 : /* -------------------------------------------------------------------- */
2279 2358 : 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 2358 : 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 2358 : if (eErr == CE_None)
2301 : {
2302 2356 : eErr = oWK.PerformWarp();
2303 2356 : ReportTiming("In memory warp operation");
2304 : }
2305 :
2306 : /* -------------------------------------------------------------------- */
2307 : /* Optional application provided postwarp chunk processor. */
2308 : /* -------------------------------------------------------------------- */
2309 2358 : 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 2358 : 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 2358 : 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 2358 : CPLFree(oWK.papabySrcImage[0]);
2342 2358 : CPLFree(oWK.papabySrcImage);
2343 2358 : CPLFree(oWK.papabyDstImage);
2344 :
2345 2358 : 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 2358 : CPLFree(oWK.panUnifiedSrcValid);
2352 2358 : CPLFree(oWK.pafUnifiedSrcDensity);
2353 2358 : CPLFree(oWK.panDstValid);
2354 2358 : CPLFree(oWK.pafDstDensity);
2355 :
2356 2358 : 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 1172 : CPLErr GDALWarpOperation::CreateKernelMask(GDALWarpKernel *poKernel, int iBand,
2392 : const char *pszType)
2393 :
2394 : {
2395 1172 : void **ppMask = nullptr;
2396 1172 : int nXSize = 0;
2397 1172 : int nYSize = 0;
2398 1172 : int nBitsPerPixel = 0;
2399 1172 : int nDefault = 0;
2400 1172 : int nExtraElts = 0;
2401 1172 : bool bDoMemset = true;
2402 :
2403 : /* -------------------------------------------------------------------- */
2404 : /* Get particulars of mask to be updated. */
2405 : /* -------------------------------------------------------------------- */
2406 1172 : if (EQUAL(pszType, "BandSrcValid"))
2407 : {
2408 168 : if (poKernel->papanBandSrcValid == nullptr)
2409 122 : poKernel->papanBandSrcValid = static_cast<GUInt32 **>(
2410 122 : CPLCalloc(sizeof(void *), poKernel->nBands));
2411 :
2412 168 : ppMask =
2413 168 : reinterpret_cast<void **>(&(poKernel->papanBandSrcValid[iBand]));
2414 168 : nExtraElts = WARP_EXTRA_ELTS;
2415 168 : nXSize = poKernel->nSrcXSize;
2416 168 : nYSize = poKernel->nSrcYSize;
2417 168 : nBitsPerPixel = 1;
2418 168 : nDefault = 0xff;
2419 : }
2420 1004 : 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 984 : 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 872 : else if (EQUAL(pszType, "DstValid"))
2440 : {
2441 627 : ppMask = reinterpret_cast<void **>(&(poKernel->panDstValid));
2442 627 : nXSize = poKernel->nDstXSize;
2443 627 : nYSize = poKernel->nDstYSize;
2444 627 : nBitsPerPixel = 1;
2445 627 : 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 1172 : if (*ppMask == nullptr)
2467 : {
2468 1172 : const GIntBig nBytes =
2469 : nBitsPerPixel == 32
2470 1172 : ? (static_cast<GIntBig>(nXSize) * nYSize + nExtraElts) * 4
2471 815 : : (static_cast<GIntBig>(nXSize) * nYSize + nExtraElts + 31) / 8;
2472 :
2473 1172 : 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 1172 : *ppMask = VSI_MALLOC_VERBOSE(nByteSize_t);
2484 :
2485 1172 : if (*ppMask == nullptr)
2486 : {
2487 0 : return CE_Failure;
2488 : }
2489 :
2490 1172 : if (bDoMemset)
2491 815 : memset(*ppMask, nDefault, nByteSize_t);
2492 : }
2493 :
2494 1172 : 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 142 : psOptions->pfnTransformer(psOptions->pTransformerArg, FALSE, nSampleMax,
2583 : privateData->adfDstX.data(),
2584 : privateData->adfDstY.data(), adfDstZ.data(),
2585 : privateData->abSuccess.data());
2586 142 : privateData->nStepCount = nStepCount;
2587 : }
2588 :
2589 : /* -------------------------------------------------------------------- */
2590 : /* Collect the bounds, ignoring any failed points. */
2591 : /* -------------------------------------------------------------------- */
2592 434 : const int nStepCount = privateData->nStepCount;
2593 434 : const double dfStepSize = 1.0 / (nStepCount - 1);
2594 434 : int iPoint = 0;
2595 : #ifdef DEBUG
2596 434 : const size_t nSampleMax =
2597 434 : static_cast<size_t>(nStepCount + 2) * (nStepCount + 2);
2598 434 : CPL_IGNORE_RET_VAL(nSampleMax);
2599 434 : CPLAssert(privateData->adfDstX.size() == nSampleMax);
2600 434 : CPLAssert(privateData->adfDstY.size() == nSampleMax);
2601 434 : CPLAssert(privateData->abSuccess.size() == nSampleMax);
2602 : #endif
2603 10416 : for (int iY = 0; iY < nStepCount + 2; iY++)
2604 : {
2605 19530 : const double dfRatioY = (iY == 0) ? 0.5 / nSrcRasterYSize
2606 : : (iY <= nStepCount)
2607 9548 : ? (iY - 1) * dfStepSize
2608 434 : : 1 - 0.5 / nSrcRasterYSize;
2609 239568 : for (int iX = 0; iX < nStepCount + 2; iX++)
2610 : {
2611 229586 : if (privateData->abSuccess[iPoint] &&
2612 220032 : privateData->adfDstX[iPoint] >= nDstXOff &&
2613 198208 : privateData->adfDstX[iPoint] <= nDstXOff + nDstXSize &&
2614 588116 : privateData->adfDstY[iPoint] >= nDstYOff &&
2615 138498 : privateData->adfDstY[iPoint] <= nDstYOff + nDstYSize)
2616 : {
2617 208886 : const double dfRatioX = (iX == 0) ? 0.5 / nSrcRasterXSize
2618 : : (iX <= nStepCount)
2619 102142 : ? (iX - 1) * dfStepSize
2620 4618 : : 1 - 0.5 / nSrcRasterXSize;
2621 106744 : double dfSrcX = dfRatioX * nSrcRasterXSize;
2622 106744 : double dfSrcY = dfRatioY * nSrcRasterYSize;
2623 106744 : *padfSrcMinX = std::min(*padfSrcMinX, dfSrcX);
2624 106744 : *padfSrcMinY = std::min(*padfSrcMinY, dfSrcY);
2625 106744 : *padfSrcMaxX = std::max(*padfSrcMaxX, dfSrcX);
2626 106744 : *padfSrcMaxY = std::max(*padfSrcMaxY, dfSrcY);
2627 : }
2628 229586 : iPoint++;
2629 : }
2630 : }
2631 : }
2632 :
2633 : /************************************************************************/
2634 : /* ComputeSourceWindowTransformPoints() */
2635 : /************************************************************************/
2636 :
2637 4705 : bool GDALWarpOperation::ComputeSourceWindowTransformPoints(
2638 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize, bool bUseGrid,
2639 : bool bAll, int nStepCount, bool bTryWithCheckWithInvertProj,
2640 : double &dfMinXOut, double &dfMinYOut, double &dfMaxXOut, double &dfMaxYOut,
2641 : int &nSamplePoints, int &nFailedCount)
2642 : {
2643 4705 : nSamplePoints = 0;
2644 4705 : nFailedCount = 0;
2645 :
2646 4705 : const double dfStepSize = bAll ? 0 : 1.0 / (nStepCount - 1);
2647 4705 : constexpr int knIntMax = std::numeric_limits<int>::max();
2648 4705 : int nSampleMax = 0;
2649 4705 : if (bUseGrid)
2650 : {
2651 494 : if (bAll)
2652 : {
2653 0 : if (nDstYSize > knIntMax / (nDstXSize + 1) - 1)
2654 : {
2655 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps");
2656 0 : return false;
2657 : }
2658 0 : nSampleMax = (nDstXSize + 1) * (nDstYSize + 1);
2659 : }
2660 : else
2661 : {
2662 494 : if (nStepCount > knIntMax - 2 ||
2663 494 : (nStepCount + 2) > knIntMax / (nStepCount + 2))
2664 : {
2665 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps : %d",
2666 : nStepCount);
2667 0 : return false;
2668 : }
2669 494 : nSampleMax = (nStepCount + 2) * (nStepCount + 2);
2670 : }
2671 : }
2672 : else
2673 : {
2674 4211 : if (bAll)
2675 : {
2676 151 : if (nDstXSize > (knIntMax - 2 * nDstYSize) / 2)
2677 : {
2678 : // Extremely unlikely !
2679 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps");
2680 0 : return false;
2681 : }
2682 151 : nSampleMax = 2 * (nDstXSize + nDstYSize);
2683 : }
2684 : else
2685 : {
2686 4060 : if (nStepCount > knIntMax / 4)
2687 : {
2688 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many steps : %d * 4",
2689 : nStepCount);
2690 0 : return false;
2691 : }
2692 4060 : nSampleMax = nStepCount * 4;
2693 : }
2694 : }
2695 :
2696 : int *pabSuccess =
2697 4705 : static_cast<int *>(VSI_MALLOC2_VERBOSE(sizeof(int), nSampleMax));
2698 : double *padfX = static_cast<double *>(
2699 4705 : VSI_MALLOC2_VERBOSE(sizeof(double) * 3, nSampleMax));
2700 4705 : if (pabSuccess == nullptr || padfX == nullptr)
2701 : {
2702 0 : CPLFree(padfX);
2703 0 : CPLFree(pabSuccess);
2704 0 : return false;
2705 : }
2706 4705 : double *padfY = padfX + nSampleMax;
2707 4705 : double *padfZ = padfX + nSampleMax * 2;
2708 :
2709 : /* -------------------------------------------------------------------- */
2710 : /* Setup sample points on a grid pattern throughout the area. */
2711 : /* -------------------------------------------------------------------- */
2712 4705 : if (bUseGrid)
2713 : {
2714 494 : if (bAll)
2715 : {
2716 0 : for (int iY = 0; iY <= nDstYSize; ++iY)
2717 : {
2718 0 : for (int iX = 0; iX <= nDstXSize; ++iX)
2719 : {
2720 0 : padfX[nSamplePoints] = nDstXOff + iX;
2721 0 : padfY[nSamplePoints] = nDstYOff + iY;
2722 0 : padfZ[nSamplePoints++] = 0.0;
2723 : }
2724 : }
2725 : }
2726 : else
2727 : {
2728 11856 : for (int iY = 0; iY < nStepCount + 2; iY++)
2729 : {
2730 22230 : const double dfRatioY = (iY == 0) ? 0.5 / nDstXSize
2731 : : (iY <= nStepCount)
2732 10868 : ? (iY - 1) * dfStepSize
2733 494 : : 1 - 0.5 / nDstXSize;
2734 272688 : for (int iX = 0; iX < nStepCount + 2; iX++)
2735 : {
2736 511290 : const double dfRatioX = (iX == 0) ? 0.5 / nDstXSize
2737 : : (iX <= nStepCount)
2738 249964 : ? (iX - 1) * dfStepSize
2739 11362 : : 1 - 0.5 / nDstXSize;
2740 261326 : padfX[nSamplePoints] = dfRatioX * nDstXSize + nDstXOff;
2741 261326 : padfY[nSamplePoints] = dfRatioY * nDstYSize + nDstYOff;
2742 261326 : padfZ[nSamplePoints++] = 0.0;
2743 : }
2744 : }
2745 : }
2746 : }
2747 : /* -------------------------------------------------------------------- */
2748 : /* Setup sample points all around the edge of the output raster. */
2749 : /* -------------------------------------------------------------------- */
2750 : else
2751 : {
2752 4211 : if (bAll)
2753 : {
2754 68258 : for (int iX = 0; iX <= nDstXSize; ++iX)
2755 : {
2756 : // Along top
2757 68107 : padfX[nSamplePoints] = nDstXOff + iX;
2758 68107 : padfY[nSamplePoints] = nDstYOff;
2759 68107 : padfZ[nSamplePoints++] = 0.0;
2760 :
2761 : // Along bottom
2762 68107 : padfX[nSamplePoints] = nDstXOff + iX;
2763 68107 : padfY[nSamplePoints] = nDstYOff + nDstYSize;
2764 68107 : padfZ[nSamplePoints++] = 0.0;
2765 : }
2766 :
2767 43329 : for (int iY = 1; iY < nDstYSize; ++iY)
2768 : {
2769 : // Along left
2770 43178 : padfX[nSamplePoints] = nDstXOff;
2771 43178 : padfY[nSamplePoints] = nDstYOff + iY;
2772 43178 : padfZ[nSamplePoints++] = 0.0;
2773 :
2774 : // Along right
2775 43178 : padfX[nSamplePoints] = nDstXOff + nDstXSize;
2776 43178 : padfY[nSamplePoints] = nDstYOff + iY;
2777 43178 : padfZ[nSamplePoints++] = 0.0;
2778 : }
2779 : }
2780 : else
2781 : {
2782 89320 : for (double dfRatio = 0.0; dfRatio <= 1.0 + dfStepSize * 0.5;
2783 85260 : dfRatio += dfStepSize)
2784 : {
2785 : // Along top
2786 85260 : padfX[nSamplePoints] = dfRatio * nDstXSize + nDstXOff;
2787 85260 : padfY[nSamplePoints] = nDstYOff;
2788 85260 : padfZ[nSamplePoints++] = 0.0;
2789 :
2790 : // Along bottom
2791 85260 : padfX[nSamplePoints] = dfRatio * nDstXSize + nDstXOff;
2792 85260 : padfY[nSamplePoints] = nDstYOff + nDstYSize;
2793 85260 : padfZ[nSamplePoints++] = 0.0;
2794 :
2795 : // Along left
2796 85260 : padfX[nSamplePoints] = nDstXOff;
2797 85260 : padfY[nSamplePoints] = dfRatio * nDstYSize + nDstYOff;
2798 85260 : padfZ[nSamplePoints++] = 0.0;
2799 :
2800 : // Along right
2801 85260 : padfX[nSamplePoints] = nDstXSize + nDstXOff;
2802 85260 : padfY[nSamplePoints] = dfRatio * nDstYSize + nDstYOff;
2803 85260 : padfZ[nSamplePoints++] = 0.0;
2804 : }
2805 : }
2806 : }
2807 :
2808 4705 : CPLAssert(nSamplePoints == nSampleMax);
2809 :
2810 : /* -------------------------------------------------------------------- */
2811 : /* Transform them to the input pixel coordinate space */
2812 : /* -------------------------------------------------------------------- */
2813 :
2814 132 : const auto RefreshTransformer = [this]()
2815 : {
2816 44 : if (GDALIsTransformer(psOptions->pTransformerArg,
2817 : GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
2818 : {
2819 0 : GDALRefreshGenImgProjTransformer(psOptions->pTransformerArg);
2820 : }
2821 44 : else if (GDALIsTransformer(psOptions->pTransformerArg,
2822 : GDAL_APPROX_TRANSFORMER_CLASS_NAME))
2823 : {
2824 44 : GDALRefreshApproxTransformer(psOptions->pTransformerArg);
2825 : }
2826 4749 : };
2827 :
2828 4705 : if (bTryWithCheckWithInvertProj)
2829 : {
2830 22 : CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ", "YES");
2831 22 : RefreshTransformer();
2832 : }
2833 4705 : psOptions->pfnTransformer(psOptions->pTransformerArg, TRUE, nSamplePoints,
2834 : padfX, padfY, padfZ, pabSuccess);
2835 4705 : if (bTryWithCheckWithInvertProj)
2836 : {
2837 22 : CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ", nullptr);
2838 22 : RefreshTransformer();
2839 : }
2840 :
2841 : /* -------------------------------------------------------------------- */
2842 : /* Collect the bounds, ignoring any failed points. */
2843 : /* -------------------------------------------------------------------- */
2844 829641 : for (int i = 0; i < nSamplePoints; i++)
2845 : {
2846 824936 : if (!pabSuccess[i])
2847 : {
2848 107923 : nFailedCount++;
2849 107923 : continue;
2850 : }
2851 :
2852 : // If this happens this is likely the symptom of a bug somewhere.
2853 717013 : if (std::isnan(padfX[i]) || std::isnan(padfY[i]))
2854 : {
2855 : static bool bNanCoordFound = false;
2856 0 : if (!bNanCoordFound)
2857 : {
2858 0 : CPLDebug("WARP",
2859 : "ComputeSourceWindow(): "
2860 : "NaN coordinate found on point %d.",
2861 : i);
2862 0 : bNanCoordFound = true;
2863 : }
2864 0 : nFailedCount++;
2865 0 : continue;
2866 : }
2867 :
2868 717013 : dfMinXOut = std::min(dfMinXOut, padfX[i]);
2869 717013 : dfMinYOut = std::min(dfMinYOut, padfY[i]);
2870 717013 : dfMaxXOut = std::max(dfMaxXOut, padfX[i]);
2871 717013 : dfMaxYOut = std::max(dfMaxYOut, padfY[i]);
2872 : }
2873 :
2874 4705 : CPLFree(padfX);
2875 4705 : CPLFree(pabSuccess);
2876 4705 : return true;
2877 : }
2878 :
2879 : /************************************************************************/
2880 : /* ComputeSourceWindow() */
2881 : /************************************************************************/
2882 :
2883 : /** Given a target window starting at pixel (nDstOff, nDstYOff) and of
2884 : * dimension (nDstXSize, nDstYSize), compute the corresponding window in
2885 : * the source raster, and return the source position in (*pnSrcXOff, *pnSrcYOff),
2886 : * the source dimension in (*pnSrcXSize, *pnSrcYSize).
2887 : * If pdfSrcXExtraSize is not null, its pointed value will be filled with the
2888 : * number of extra source pixels in X dimension to acquire to take into account
2889 : * the size of the resampling kernel. Similarly for pdfSrcYExtraSize for the
2890 : * Y dimension.
2891 : * If pdfSrcFillRatio is not null, its pointed value will be filled with the
2892 : * the ratio of the clamped source raster window size over the unclamped source
2893 : * raster window size. When this ratio is too low, this might be an indication
2894 : * that it might be beneficial to split the target window to avoid requesting
2895 : * too many source pixels.
2896 : */
2897 4399 : CPLErr GDALWarpOperation::ComputeSourceWindow(
2898 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize, int *pnSrcXOff,
2899 : int *pnSrcYOff, int *pnSrcXSize, int *pnSrcYSize, double *pdfSrcXExtraSize,
2900 : double *pdfSrcYExtraSize, double *pdfSrcFillRatio)
2901 :
2902 : {
2903 : /* -------------------------------------------------------------------- */
2904 : /* Figure out whether we just want to do the usual "along the */
2905 : /* edge" sampling, or using a grid. The grid usage is */
2906 : /* important in some weird "inside out" cases like WGS84 to */
2907 : /* polar stereographic around the pole. Also figure out the */
2908 : /* sampling rate. */
2909 : /* -------------------------------------------------------------------- */
2910 4399 : int nStepCount = DEFAULT_STEP_COUNT;
2911 4399 : bool bAll = false;
2912 :
2913 : bool bUseGrid =
2914 4399 : CPLFetchBool(psOptions->papszWarpOptions, "SAMPLE_GRID", false);
2915 :
2916 : const char *pszSampleSteps =
2917 4399 : CSLFetchNameValue(psOptions->papszWarpOptions, "SAMPLE_STEPS");
2918 4399 : if (pszSampleSteps)
2919 : {
2920 94 : if (EQUAL(pszSampleSteps, "ALL"))
2921 : {
2922 94 : bAll = true;
2923 : }
2924 : else
2925 : {
2926 0 : nStepCount = atoi(pszSampleSteps);
2927 0 : nStepCount = std::max(2, nStepCount);
2928 : }
2929 : }
2930 4305 : else if (!bUseGrid)
2931 : {
2932 : // Detect if at least one of the 4 corner in destination raster fails
2933 : // to project back to source.
2934 : // Helps for long-lat to orthographic on areas that are partly in
2935 : // space / partly on Earth. Cf https://github.com/OSGeo/gdal/issues/9056
2936 : double adfCornerX[4];
2937 : double adfCornerY[4];
2938 4117 : double adfCornerZ[4] = {0, 0, 0, 0};
2939 4117 : int anCornerSuccess[4] = {FALSE, FALSE, FALSE, FALSE};
2940 4117 : adfCornerX[0] = nDstXOff;
2941 4117 : adfCornerY[0] = nDstYOff;
2942 4117 : adfCornerX[1] = nDstXOff + nDstXSize;
2943 4117 : adfCornerY[1] = nDstYOff;
2944 4117 : adfCornerX[2] = nDstXOff;
2945 4117 : adfCornerY[2] = nDstYOff + nDstYSize;
2946 4117 : adfCornerX[3] = nDstXOff + nDstXSize;
2947 4117 : adfCornerY[3] = nDstYOff + nDstYSize;
2948 4117 : if (!psOptions->pfnTransformer(psOptions->pTransformerArg, TRUE, 4,
2949 : adfCornerX, adfCornerY, adfCornerZ,
2950 4060 : anCornerSuccess) ||
2951 8177 : !anCornerSuccess[0] || !anCornerSuccess[1] || !anCornerSuccess[2] ||
2952 4060 : !anCornerSuccess[3])
2953 : {
2954 57 : bAll = true;
2955 : }
2956 : }
2957 :
2958 4399 : bool bTryWithCheckWithInvertProj = false;
2959 4399 : double dfMinXOut = std::numeric_limits<double>::infinity();
2960 4399 : double dfMinYOut = std::numeric_limits<double>::infinity();
2961 4399 : double dfMaxXOut = -std::numeric_limits<double>::infinity();
2962 4399 : double dfMaxYOut = -std::numeric_limits<double>::infinity();
2963 :
2964 4399 : int nSamplePoints = 0;
2965 4399 : int nFailedCount = 0;
2966 4399 : if (!ComputeSourceWindowTransformPoints(
2967 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, bUseGrid, bAll,
2968 : nStepCount, bTryWithCheckWithInvertProj, dfMinXOut, dfMinYOut,
2969 : dfMaxXOut, dfMaxYOut, nSamplePoints, nFailedCount))
2970 : {
2971 0 : return CE_Failure;
2972 : }
2973 :
2974 : // Use grid sampling as soon as a special point falls into the extent of
2975 : // the target raster.
2976 4399 : if (!bUseGrid && psOptions->hDstDS)
2977 : {
2978 10604 : for (const auto &xy : aDstXYSpecialPoints)
2979 : {
2980 13984 : if (0 <= xy.first &&
2981 718 : GDALGetRasterXSize(psOptions->hDstDS) >= xy.first &&
2982 7651 : 0 <= xy.second &&
2983 300 : GDALGetRasterYSize(psOptions->hDstDS) >= xy.second)
2984 : {
2985 227 : bUseGrid = true;
2986 227 : bAll = false;
2987 227 : if (!ComputeSourceWindowTransformPoints(
2988 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, bUseGrid,
2989 : bAll, nStepCount, bTryWithCheckWithInvertProj,
2990 : dfMinXOut, dfMinYOut, dfMaxXOut, dfMaxYOut,
2991 : nSamplePoints, nFailedCount))
2992 : {
2993 0 : return CE_Failure;
2994 : }
2995 227 : break;
2996 : }
2997 : }
2998 : }
2999 :
3000 4399 : const int nRasterXSize = GDALGetRasterXSize(psOptions->hSrcDS);
3001 4399 : const int nRasterYSize = GDALGetRasterYSize(psOptions->hSrcDS);
3002 :
3003 : // Try to detect crazy values coming from reprojection that would not
3004 : // have resulted in a PROJ error. Could happen for example with PROJ
3005 : // <= 4.9.2 with inverse UTM/tmerc (Snyder approximation without sanity
3006 : // check) when being far away from the central meridian. But might be worth
3007 : // keeping that even for later versions in case some exotic projection isn't
3008 : // properly sanitized.
3009 4334 : if (nFailedCount == 0 && !bTryWithCheckWithInvertProj &&
3010 4334 : (dfMinXOut < -1e6 || dfMinYOut < -1e6 ||
3011 8737 : dfMaxXOut > nRasterXSize + 1e6 || dfMaxYOut > nRasterYSize + 1e6) &&
3012 22 : !CPLTestBool(CPLGetConfigOption("CHECK_WITH_INVERT_PROJ", "NO")))
3013 : {
3014 22 : CPLDebug("WARP",
3015 : "ComputeSourceWindow(): bogus source dataset window "
3016 : "returned. Trying again with CHECK_WITH_INVERT_PROJ=YES");
3017 22 : bTryWithCheckWithInvertProj = true;
3018 :
3019 : // We should probably perform the coordinate transformation in the
3020 : // warp kernel under CHECK_WITH_INVERT_PROJ too...
3021 22 : if (!ComputeSourceWindowTransformPoints(
3022 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, bUseGrid, bAll,
3023 : nStepCount, bTryWithCheckWithInvertProj, dfMinXOut, dfMinYOut,
3024 : dfMaxXOut, dfMaxYOut, nSamplePoints, nFailedCount))
3025 : {
3026 0 : return CE_Failure;
3027 : }
3028 : }
3029 :
3030 : /* -------------------------------------------------------------------- */
3031 : /* If we got any failures when not using a grid, we should */
3032 : /* really go back and try again with the grid. Sorry for the */
3033 : /* goto. */
3034 : /* -------------------------------------------------------------------- */
3035 4399 : if (!bUseGrid && nFailedCount > 0)
3036 : {
3037 57 : bUseGrid = true;
3038 57 : bAll = false;
3039 57 : if (!ComputeSourceWindowTransformPoints(
3040 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, bUseGrid, bAll,
3041 : nStepCount, bTryWithCheckWithInvertProj, dfMinXOut, dfMinYOut,
3042 : dfMaxXOut, dfMaxYOut, nSamplePoints, nFailedCount))
3043 : {
3044 0 : return CE_Failure;
3045 : }
3046 : }
3047 :
3048 : /* -------------------------------------------------------------------- */
3049 : /* If we get hardly any points (or none) transforming, we give */
3050 : /* up. */
3051 : /* -------------------------------------------------------------------- */
3052 4399 : if (nFailedCount > nSamplePoints - 5)
3053 : {
3054 : const bool bErrorOutIfEmptySourceWindow =
3055 38 : CPLFetchBool(psOptions->papszWarpOptions,
3056 : "ERROR_OUT_IF_EMPTY_SOURCE_WINDOW", true);
3057 38 : if (bErrorOutIfEmptySourceWindow)
3058 : {
3059 2 : CPLError(CE_Failure, CPLE_AppDefined,
3060 : "Too many points (%d out of %d) failed to transform, "
3061 : "unable to compute output bounds.",
3062 : nFailedCount, nSamplePoints);
3063 : }
3064 : else
3065 : {
3066 36 : CPLDebug("WARP", "Cannot determine source window for %d,%d,%d,%d",
3067 : nDstXOff, nDstYOff, nDstXSize, nDstYSize);
3068 : }
3069 38 : return CE_Failure;
3070 : }
3071 :
3072 4361 : if (nFailedCount > 0)
3073 31 : CPLDebug("GDAL",
3074 : "GDALWarpOperation::ComputeSourceWindow() %d out of %d "
3075 : "points failed to transform.",
3076 : nFailedCount, nSamplePoints);
3077 :
3078 : /* -------------------------------------------------------------------- */
3079 : /* In some cases (see https://github.com/OSGeo/gdal/issues/862) */
3080 : /* the reverse transform does not work at some points, so try by */
3081 : /* transforming from source raster space to target raster space and */
3082 : /* see which source coordinates end up being in the AOI in the target */
3083 : /* raster space. */
3084 : /* -------------------------------------------------------------------- */
3085 4361 : if (bUseGrid)
3086 : {
3087 434 : ComputeSourceWindowStartingFromSource(nDstXOff, nDstYOff, nDstXSize,
3088 : nDstYSize, &dfMinXOut, &dfMinYOut,
3089 : &dfMaxXOut, &dfMaxYOut);
3090 : }
3091 :
3092 : /* -------------------------------------------------------------------- */
3093 : /* Early exit to avoid crazy values to cause a huge nResWinSize that */
3094 : /* would result in a result window wrongly covering the whole raster. */
3095 : /* -------------------------------------------------------------------- */
3096 4361 : if (dfMinXOut > nRasterXSize || dfMaxXOut < 0 || dfMinYOut > nRasterYSize ||
3097 3752 : dfMaxYOut < 0)
3098 : {
3099 1012 : *pnSrcXOff = 0;
3100 1012 : *pnSrcYOff = 0;
3101 1012 : *pnSrcXSize = 0;
3102 1012 : *pnSrcYSize = 0;
3103 1012 : if (pdfSrcXExtraSize)
3104 1012 : *pdfSrcXExtraSize = 0.0;
3105 1012 : if (pdfSrcYExtraSize)
3106 1012 : *pdfSrcYExtraSize = 0.0;
3107 1012 : if (pdfSrcFillRatio)
3108 745 : *pdfSrcFillRatio = 0.0;
3109 1012 : return CE_None;
3110 : }
3111 :
3112 : // For scenarios where warping is used as a "decoration", try to clamp
3113 : // source pixel coordinates to integer when very close.
3114 13396 : const auto roundIfCloseEnough = [](double dfVal)
3115 : {
3116 13396 : const double dfRounded = std::round(dfVal);
3117 13396 : if (std::fabs(dfRounded - dfVal) < 1e-6)
3118 11141 : return dfRounded;
3119 2255 : return dfVal;
3120 : };
3121 :
3122 3349 : dfMinXOut = roundIfCloseEnough(dfMinXOut);
3123 3349 : dfMinYOut = roundIfCloseEnough(dfMinYOut);
3124 3349 : dfMaxXOut = roundIfCloseEnough(dfMaxXOut);
3125 3349 : dfMaxYOut = roundIfCloseEnough(dfMaxYOut);
3126 :
3127 3349 : if (m_bIsTranslationOnPixelBoundaries)
3128 : {
3129 182 : CPLAssert(dfMinXOut == std::round(dfMinXOut));
3130 182 : CPLAssert(dfMinYOut == std::round(dfMinYOut));
3131 182 : CPLAssert(dfMaxXOut == std::round(dfMaxXOut));
3132 182 : CPLAssert(dfMaxYOut == std::round(dfMaxYOut));
3133 182 : CPLAssert(std::round(dfMaxXOut - dfMinXOut) == nDstXSize);
3134 182 : CPLAssert(std::round(dfMaxYOut - dfMinYOut) == nDstYSize);
3135 : }
3136 :
3137 : /* -------------------------------------------------------------------- */
3138 : /* How much of a window around our source pixel might we need */
3139 : /* to collect data from based on the resampling kernel? Even */
3140 : /* if the requested central pixel falls off the source image, */
3141 : /* we may need to collect data if some portion of the */
3142 : /* resampling kernel could be on-image. */
3143 : /* -------------------------------------------------------------------- */
3144 3349 : const int nResWinSize = m_bIsTranslationOnPixelBoundaries
3145 3349 : ? 0
3146 3167 : : GWKGetFilterRadius(psOptions->eResampleAlg);
3147 :
3148 : // Take scaling into account.
3149 : // Avoid ridiculous small scaling factors to avoid potential further integer
3150 : // overflows
3151 6698 : const double dfXScale = std::max(1e-3, static_cast<double>(nDstXSize) /
3152 3349 : (dfMaxXOut - dfMinXOut));
3153 6698 : const double dfYScale = std::max(1e-3, static_cast<double>(nDstYSize) /
3154 3349 : (dfMaxYOut - dfMinYOut));
3155 3349 : int nXRadius = dfXScale < 0.95
3156 3349 : ? static_cast<int>(ceil(nResWinSize / dfXScale))
3157 : : nResWinSize;
3158 3349 : int nYRadius = dfYScale < 0.95
3159 3349 : ? static_cast<int>(ceil(nResWinSize / dfYScale))
3160 : : nResWinSize;
3161 :
3162 : /* -------------------------------------------------------------------- */
3163 : /* Allow addition of extra sample pixels to source window to */
3164 : /* avoid missing pixels due to sampling error. In fact, */
3165 : /* fallback to adding a bit to the window if any points failed */
3166 : /* to transform. */
3167 : /* -------------------------------------------------------------------- */
3168 3349 : if (CSLFetchNameValue(psOptions->papszWarpOptions, "SOURCE_EXTRA") !=
3169 : nullptr)
3170 : {
3171 160 : const int nSrcExtra = atoi(
3172 80 : CSLFetchNameValue(psOptions->papszWarpOptions, "SOURCE_EXTRA"));
3173 80 : nXRadius += nSrcExtra;
3174 80 : nYRadius += nSrcExtra;
3175 : }
3176 3269 : else if (nFailedCount > 0)
3177 : {
3178 25 : nXRadius += 10;
3179 25 : nYRadius += 10;
3180 : }
3181 :
3182 : /* -------------------------------------------------------------------- */
3183 : /* return bounds. */
3184 : /* -------------------------------------------------------------------- */
3185 : #if DEBUG_VERBOSE
3186 : CPLDebug("WARP",
3187 : "dst=(%d,%d,%d,%d) raw "
3188 : "src=(minx=%.17g,miny=%.17g,maxx=%.17g,maxy=%.17g)",
3189 : nDstXOff, nDstYOff, nDstXSize, nDstYSize, dfMinXOut, dfMinYOut,
3190 : dfMaxXOut, dfMaxYOut);
3191 : #endif
3192 3349 : const int nMinXOutClamped = static_cast<int>(std::max(0.0, dfMinXOut));
3193 3349 : const int nMinYOutClamped = static_cast<int>(std::max(0.0, dfMinYOut));
3194 : const int nMaxXOutClamped = static_cast<int>(
3195 3349 : std::min(ceil(dfMaxXOut), static_cast<double>(nRasterXSize)));
3196 : const int nMaxYOutClamped = static_cast<int>(
3197 3349 : std::min(ceil(dfMaxYOut), static_cast<double>(nRasterYSize)));
3198 :
3199 : const double dfSrcXSizeRaw = std::max(
3200 10047 : 0.0, std::min(static_cast<double>(nRasterXSize - nMinXOutClamped),
3201 3349 : dfMaxXOut - dfMinXOut));
3202 : const double dfSrcYSizeRaw = std::max(
3203 10047 : 0.0, std::min(static_cast<double>(nRasterYSize - nMinYOutClamped),
3204 3349 : dfMaxYOut - dfMinYOut));
3205 :
3206 : // If we cover more than 90% of the width, then use it fully (helps for
3207 : // anti-meridian discontinuities)
3208 3349 : if (nMaxXOutClamped - nMinXOutClamped > 0.9 * nRasterXSize)
3209 : {
3210 1308 : *pnSrcXOff = 0;
3211 1308 : *pnSrcXSize = nRasterXSize;
3212 : }
3213 : else
3214 : {
3215 2041 : *pnSrcXOff =
3216 2041 : std::max(0, std::min(nMinXOutClamped - nXRadius, nRasterXSize));
3217 2041 : *pnSrcXSize =
3218 6123 : std::max(0, std::min(nRasterXSize - *pnSrcXOff,
3219 2041 : nMaxXOutClamped - *pnSrcXOff + nXRadius));
3220 : }
3221 :
3222 3349 : if (nMaxYOutClamped - nMinYOutClamped > 0.9 * nRasterYSize)
3223 : {
3224 1190 : *pnSrcYOff = 0;
3225 1190 : *pnSrcYSize = nRasterYSize;
3226 : }
3227 : else
3228 : {
3229 2159 : *pnSrcYOff =
3230 2159 : std::max(0, std::min(nMinYOutClamped - nYRadius, nRasterYSize));
3231 2159 : *pnSrcYSize =
3232 6477 : std::max(0, std::min(nRasterYSize - *pnSrcYOff,
3233 2159 : nMaxYOutClamped - *pnSrcYOff + nYRadius));
3234 : }
3235 :
3236 3349 : if (pdfSrcXExtraSize)
3237 3349 : *pdfSrcXExtraSize = *pnSrcXSize - dfSrcXSizeRaw;
3238 3349 : if (pdfSrcYExtraSize)
3239 3349 : *pdfSrcYExtraSize = *pnSrcYSize - dfSrcYSizeRaw;
3240 :
3241 : // Computed the ratio of the clamped source raster window size over
3242 : // the unclamped source raster window size.
3243 3349 : if (pdfSrcFillRatio)
3244 3156 : *pdfSrcFillRatio =
3245 6312 : static_cast<double>(*pnSrcXSize) * (*pnSrcYSize) /
3246 3156 : std::max(1.0, (dfMaxXOut - dfMinXOut + 2 * nXRadius) *
3247 3156 : (dfMaxYOut - dfMinYOut + 2 * nYRadius));
3248 :
3249 3349 : return CE_None;
3250 : }
3251 :
3252 : /************************************************************************/
3253 : /* ReportTiming() */
3254 : /************************************************************************/
3255 :
3256 9091 : void GDALWarpOperation::ReportTiming(const char *pszMessage)
3257 :
3258 : {
3259 9091 : if (!bReportTimings)
3260 9091 : return;
3261 :
3262 0 : const unsigned long nNewTime = VSITime(nullptr);
3263 :
3264 0 : if (pszMessage != nullptr)
3265 : {
3266 0 : CPLDebug("WARP_TIMING", "%s: %lds", pszMessage,
3267 0 : static_cast<long>(nNewTime - nLastTimeReported));
3268 : }
3269 :
3270 0 : nLastTimeReported = nNewTime;
3271 : }
|
