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