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