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