/src/gdal/alg/gdalwarper.cpp
Line | Count | Source (jump to first uncovered line) |
1 | | /****************************************************************************** |
2 | | * |
3 | | * Project: High Performance Image Reprojector |
4 | | * Purpose: Implementation of high level convenience APIs for warper. |
5 | | * Author: Frank Warmerdam, warmerdam@pobox.com |
6 | | * |
7 | | ****************************************************************************** |
8 | | * Copyright (c) 2003, Frank Warmerdam <warmerdam@pobox.com> |
9 | | * Copyright (c) 2008-2012, Even Rouault <even dot rouault at spatialys.com> |
10 | | * |
11 | | * SPDX-License-Identifier: MIT |
12 | | ****************************************************************************/ |
13 | | |
14 | | #include "cpl_port.h" |
15 | | #include "gdalwarper.h" |
16 | | |
17 | | #include <stdlib.h> |
18 | | #include <string.h> |
19 | | |
20 | | #include <algorithm> |
21 | | #include <cmath> |
22 | | #include <limits> |
23 | | |
24 | | #include "cpl_conv.h" |
25 | | #include "cpl_error.h" |
26 | | #include "cpl_float.h" |
27 | | #include "cpl_mask.h" |
28 | | #include "cpl_minixml.h" |
29 | | #include "cpl_progress.h" |
30 | | #include "cpl_string.h" |
31 | | #include "cpl_vsi.h" |
32 | | #include "gdal.h" |
33 | | #include "gdal_priv.h" |
34 | | #include "ogr_api.h" |
35 | | #include "ogr_core.h" |
36 | | #include "vrtdataset.h" // for VRTSerializeNoData |
37 | | |
38 | | #if (defined(__x86_64) || defined(_M_X64)) |
39 | | #include <emmintrin.h> |
40 | | #endif |
41 | | |
42 | | /************************************************************************/ |
43 | | /* GDALReprojectImage() */ |
44 | | /************************************************************************/ |
45 | | |
46 | | /** |
47 | | * Reproject image. |
48 | | * |
49 | | * This is a convenience function utilizing the GDALWarpOperation class to |
50 | | * reproject an image from a source to a destination. In particular, this |
51 | | * function takes care of establishing the transformation function to |
52 | | * implement the reprojection, and will default a variety of other |
53 | | * warp options. |
54 | | * |
55 | | * Nodata values set on destination dataset are taken into account. |
56 | | * |
57 | | * No metadata, projection info, or color tables are transferred |
58 | | * to the output file. Source overviews are not considered. |
59 | | * |
60 | | * For more advanced warping capabilities, consider using GDALWarp(). |
61 | | * |
62 | | * @param hSrcDS the source image file. |
63 | | * @param pszSrcWKT the source projection. If NULL the source projection |
64 | | * is read from from hSrcDS. |
65 | | * @param hDstDS the destination image file. |
66 | | * @param pszDstWKT the destination projection. If NULL the destination |
67 | | * projection will be read from hDstDS. |
68 | | * @param eResampleAlg the type of resampling to use. |
69 | | * @param dfWarpMemoryLimit the amount of memory (in bytes) that the warp |
70 | | * API is allowed to use for caching. This is in addition to the memory |
71 | | * already allocated to the GDAL caching (as per GDALSetCacheMax()). May be |
72 | | * 0.0 to use default memory settings. |
73 | | * @param dfMaxError maximum error measured in input pixels that is allowed |
74 | | * in approximating the transformation (0.0 for exact calculations). |
75 | | * @param pfnProgress a GDALProgressFunc() compatible callback function for |
76 | | * reporting progress or NULL. |
77 | | * @param pProgressArg argument to be passed to pfnProgress. May be NULL. |
78 | | * @param psOptions warp options, normally NULL. |
79 | | * |
80 | | * @return CE_None on success or CE_Failure if something goes wrong. |
81 | | * @see GDALWarp() |
82 | | */ |
83 | | |
84 | | CPLErr CPL_STDCALL GDALReprojectImage( |
85 | | GDALDatasetH hSrcDS, const char *pszSrcWKT, GDALDatasetH hDstDS, |
86 | | const char *pszDstWKT, GDALResampleAlg eResampleAlg, |
87 | | CPL_UNUSED double dfWarpMemoryLimit, double dfMaxError, |
88 | | GDALProgressFunc pfnProgress, void *pProgressArg, |
89 | | GDALWarpOptions *psOptions) |
90 | | |
91 | 0 | { |
92 | | /* -------------------------------------------------------------------- */ |
93 | | /* Setup a reprojection based transformer. */ |
94 | | /* -------------------------------------------------------------------- */ |
95 | 0 | void *hTransformArg = GDALCreateGenImgProjTransformer( |
96 | 0 | hSrcDS, pszSrcWKT, hDstDS, pszDstWKT, TRUE, 1000.0, 0); |
97 | |
|
98 | 0 | if (hTransformArg == nullptr) |
99 | 0 | return CE_Failure; |
100 | | |
101 | | /* -------------------------------------------------------------------- */ |
102 | | /* Create a copy of the user provided options, or a defaulted */ |
103 | | /* options structure. */ |
104 | | /* -------------------------------------------------------------------- */ |
105 | 0 | GDALWarpOptions *psWOptions = psOptions == nullptr |
106 | 0 | ? GDALCreateWarpOptions() |
107 | 0 | : GDALCloneWarpOptions(psOptions); |
108 | |
|
109 | 0 | psWOptions->eResampleAlg = eResampleAlg; |
110 | | |
111 | | /* -------------------------------------------------------------------- */ |
112 | | /* Set transform. */ |
113 | | /* -------------------------------------------------------------------- */ |
114 | 0 | if (dfMaxError > 0.0) |
115 | 0 | { |
116 | 0 | psWOptions->pTransformerArg = GDALCreateApproxTransformer( |
117 | 0 | GDALGenImgProjTransform, hTransformArg, dfMaxError); |
118 | |
|
119 | 0 | psWOptions->pfnTransformer = GDALApproxTransform; |
120 | 0 | } |
121 | 0 | else |
122 | 0 | { |
123 | 0 | psWOptions->pfnTransformer = GDALGenImgProjTransform; |
124 | 0 | psWOptions->pTransformerArg = hTransformArg; |
125 | 0 | } |
126 | | |
127 | | /* -------------------------------------------------------------------- */ |
128 | | /* Set file and band mapping. */ |
129 | | /* -------------------------------------------------------------------- */ |
130 | 0 | psWOptions->hSrcDS = hSrcDS; |
131 | 0 | psWOptions->hDstDS = hDstDS; |
132 | |
|
133 | 0 | int nSrcBands = GDALGetRasterCount(hSrcDS); |
134 | 0 | { |
135 | 0 | GDALRasterBandH hBand = GDALGetRasterBand(hSrcDS, nSrcBands); |
136 | 0 | if (hBand && GDALGetRasterColorInterpretation(hBand) == GCI_AlphaBand) |
137 | 0 | { |
138 | 0 | psWOptions->nSrcAlphaBand = nSrcBands; |
139 | 0 | nSrcBands--; |
140 | 0 | } |
141 | 0 | } |
142 | |
|
143 | 0 | int nDstBands = GDALGetRasterCount(hDstDS); |
144 | 0 | { |
145 | 0 | GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, nDstBands); |
146 | 0 | if (hBand && GDALGetRasterColorInterpretation(hBand) == GCI_AlphaBand) |
147 | 0 | { |
148 | 0 | psWOptions->nDstAlphaBand = nDstBands; |
149 | 0 | nDstBands--; |
150 | 0 | } |
151 | 0 | } |
152 | |
|
153 | 0 | GDALWarpInitDefaultBandMapping(psWOptions, std::min(nSrcBands, nDstBands)); |
154 | | |
155 | | /* -------------------------------------------------------------------- */ |
156 | | /* Set source nodata values if the source dataset seems to have */ |
157 | | /* any. Same for target nodata values */ |
158 | | /* -------------------------------------------------------------------- */ |
159 | 0 | for (int iBand = 0; iBand < psWOptions->nBandCount; iBand++) |
160 | 0 | { |
161 | 0 | GDALRasterBandH hBand = GDALGetRasterBand(hSrcDS, iBand + 1); |
162 | |
|
163 | 0 | int bGotNoData = FALSE; |
164 | 0 | double dfNoDataValue = GDALGetRasterNoDataValue(hBand, &bGotNoData); |
165 | 0 | if (bGotNoData) |
166 | 0 | { |
167 | 0 | GDALWarpInitSrcNoDataReal(psWOptions, -1.1e20); |
168 | 0 | psWOptions->padfSrcNoDataReal[iBand] = dfNoDataValue; |
169 | 0 | } |
170 | | |
171 | | // Deal with target band. |
172 | 0 | hBand = GDALGetRasterBand(hDstDS, iBand + 1); |
173 | |
|
174 | 0 | dfNoDataValue = GDALGetRasterNoDataValue(hBand, &bGotNoData); |
175 | 0 | if (bGotNoData) |
176 | 0 | { |
177 | 0 | GDALWarpInitDstNoDataReal(psWOptions, -1.1e20); |
178 | 0 | psWOptions->padfDstNoDataReal[iBand] = dfNoDataValue; |
179 | 0 | } |
180 | 0 | } |
181 | | |
182 | | /* -------------------------------------------------------------------- */ |
183 | | /* Set the progress function. */ |
184 | | /* -------------------------------------------------------------------- */ |
185 | 0 | if (pfnProgress != nullptr) |
186 | 0 | { |
187 | 0 | psWOptions->pfnProgress = pfnProgress; |
188 | 0 | psWOptions->pProgressArg = pProgressArg; |
189 | 0 | } |
190 | | |
191 | | /* -------------------------------------------------------------------- */ |
192 | | /* Create a warp options based on the options. */ |
193 | | /* -------------------------------------------------------------------- */ |
194 | 0 | GDALWarpOperation oWarper; |
195 | 0 | CPLErr eErr = oWarper.Initialize(psWOptions); |
196 | |
|
197 | 0 | if (eErr == CE_None) |
198 | 0 | eErr = oWarper.ChunkAndWarpImage(0, 0, GDALGetRasterXSize(hDstDS), |
199 | 0 | GDALGetRasterYSize(hDstDS)); |
200 | | |
201 | | /* -------------------------------------------------------------------- */ |
202 | | /* Cleanup. */ |
203 | | /* -------------------------------------------------------------------- */ |
204 | 0 | GDALDestroyGenImgProjTransformer(hTransformArg); |
205 | |
|
206 | 0 | if (dfMaxError > 0.0) |
207 | 0 | GDALDestroyApproxTransformer(psWOptions->pTransformerArg); |
208 | |
|
209 | 0 | GDALDestroyWarpOptions(psWOptions); |
210 | |
|
211 | 0 | return eErr; |
212 | 0 | } |
213 | | |
214 | | /************************************************************************/ |
215 | | /* GDALCreateAndReprojectImage() */ |
216 | | /* */ |
217 | | /* This is a "quicky" reprojection API. */ |
218 | | /************************************************************************/ |
219 | | |
220 | | /** Reproject an image and create the target reprojected image */ |
221 | | CPLErr CPL_STDCALL GDALCreateAndReprojectImage( |
222 | | GDALDatasetH hSrcDS, const char *pszSrcWKT, const char *pszDstFilename, |
223 | | const char *pszDstWKT, GDALDriverH hDstDriver, char **papszCreateOptions, |
224 | | GDALResampleAlg eResampleAlg, double dfWarpMemoryLimit, double dfMaxError, |
225 | | GDALProgressFunc pfnProgress, void *pProgressArg, |
226 | | GDALWarpOptions *psOptions) |
227 | | |
228 | 0 | { |
229 | 0 | VALIDATE_POINTER1(hSrcDS, "GDALCreateAndReprojectImage", CE_Failure); |
230 | | |
231 | | /* -------------------------------------------------------------------- */ |
232 | | /* Default a few parameters. */ |
233 | | /* -------------------------------------------------------------------- */ |
234 | 0 | if (hDstDriver == nullptr) |
235 | 0 | { |
236 | 0 | hDstDriver = GDALGetDriverByName("GTiff"); |
237 | 0 | if (hDstDriver == nullptr) |
238 | 0 | { |
239 | 0 | CPLError(CE_Failure, CPLE_AppDefined, |
240 | 0 | "GDALCreateAndReprojectImage needs GTiff driver"); |
241 | 0 | return CE_Failure; |
242 | 0 | } |
243 | 0 | } |
244 | | |
245 | 0 | if (pszSrcWKT == nullptr) |
246 | 0 | pszSrcWKT = GDALGetProjectionRef(hSrcDS); |
247 | |
|
248 | 0 | if (pszDstWKT == nullptr) |
249 | 0 | pszDstWKT = pszSrcWKT; |
250 | | |
251 | | /* -------------------------------------------------------------------- */ |
252 | | /* Create a transformation object from the source to */ |
253 | | /* destination coordinate system. */ |
254 | | /* -------------------------------------------------------------------- */ |
255 | 0 | void *hTransformArg = GDALCreateGenImgProjTransformer( |
256 | 0 | hSrcDS, pszSrcWKT, nullptr, pszDstWKT, TRUE, 1000.0, 0); |
257 | |
|
258 | 0 | if (hTransformArg == nullptr) |
259 | 0 | return CE_Failure; |
260 | | |
261 | | /* -------------------------------------------------------------------- */ |
262 | | /* Get approximate output definition. */ |
263 | | /* -------------------------------------------------------------------- */ |
264 | 0 | double adfDstGeoTransform[6] = {}; |
265 | 0 | int nPixels = 0; |
266 | 0 | int nLines = 0; |
267 | |
|
268 | 0 | if (GDALSuggestedWarpOutput(hSrcDS, GDALGenImgProjTransform, hTransformArg, |
269 | 0 | adfDstGeoTransform, &nPixels, |
270 | 0 | &nLines) != CE_None) |
271 | 0 | return CE_Failure; |
272 | | |
273 | 0 | GDALDestroyGenImgProjTransformer(hTransformArg); |
274 | | |
275 | | /* -------------------------------------------------------------------- */ |
276 | | /* Create the output file. */ |
277 | | /* -------------------------------------------------------------------- */ |
278 | 0 | GDALDatasetH hDstDS = GDALCreate( |
279 | 0 | hDstDriver, pszDstFilename, nPixels, nLines, GDALGetRasterCount(hSrcDS), |
280 | 0 | GDALGetRasterDataType(GDALGetRasterBand(hSrcDS, 1)), |
281 | 0 | papszCreateOptions); |
282 | |
|
283 | 0 | if (hDstDS == nullptr) |
284 | 0 | return CE_Failure; |
285 | | |
286 | | /* -------------------------------------------------------------------- */ |
287 | | /* Write out the projection definition. */ |
288 | | /* -------------------------------------------------------------------- */ |
289 | 0 | GDALSetProjection(hDstDS, pszDstWKT); |
290 | 0 | GDALSetGeoTransform(hDstDS, adfDstGeoTransform); |
291 | | |
292 | | /* -------------------------------------------------------------------- */ |
293 | | /* Perform the reprojection. */ |
294 | | /* -------------------------------------------------------------------- */ |
295 | 0 | CPLErr eErr = GDALReprojectImage( |
296 | 0 | hSrcDS, pszSrcWKT, hDstDS, pszDstWKT, eResampleAlg, dfWarpMemoryLimit, |
297 | 0 | dfMaxError, pfnProgress, pProgressArg, psOptions); |
298 | |
|
299 | 0 | GDALClose(hDstDS); |
300 | |
|
301 | 0 | return eErr; |
302 | 0 | } |
303 | | |
304 | | /************************************************************************/ |
305 | | /* GDALWarpNoDataMaskerT() */ |
306 | | /************************************************************************/ |
307 | | |
308 | | template <class T> |
309 | | static CPLErr GDALWarpNoDataMaskerT(const double *padfNoData, size_t nPixels, |
310 | | const T *pData, GUInt32 *panValidityMask, |
311 | | int *pbOutAllValid) |
312 | 0 | { |
313 | | // Nothing to do if value is out of range. |
314 | 0 | if (padfNoData[0] < cpl::NumericLimits<T>::min() || |
315 | 0 | padfNoData[0] > cpl::NumericLimits<T>::max() + 0.000001 || |
316 | 0 | padfNoData[1] != 0.0) |
317 | 0 | { |
318 | 0 | *pbOutAllValid = TRUE; |
319 | 0 | return CE_None; |
320 | 0 | } |
321 | | |
322 | 0 | const int nNoData = static_cast<int>(floor(padfNoData[0] + 0.000001)); |
323 | 0 | int bAllValid = TRUE; |
324 | 0 | for (size_t iOffset = 0; iOffset < nPixels; ++iOffset) |
325 | 0 | { |
326 | 0 | if (pData[iOffset] == nNoData) |
327 | 0 | { |
328 | 0 | bAllValid = FALSE; |
329 | 0 | CPLMaskClear(panValidityMask, iOffset); |
330 | 0 | } |
331 | 0 | } |
332 | 0 | *pbOutAllValid = bAllValid; |
333 | |
|
334 | 0 | return CE_None; |
335 | 0 | } Unexecuted instantiation: gdalwarper.cpp:CPLErr GDALWarpNoDataMaskerT<unsigned char>(double const*, unsigned long, unsigned char const*, unsigned int*, int*) Unexecuted instantiation: gdalwarper.cpp:CPLErr GDALWarpNoDataMaskerT<short>(double const*, unsigned long, short const*, unsigned int*, int*) Unexecuted instantiation: gdalwarper.cpp:CPLErr GDALWarpNoDataMaskerT<unsigned short>(double const*, unsigned long, unsigned short const*, unsigned int*, int*) |
336 | | |
337 | | /************************************************************************/ |
338 | | /* GDALWarpNoDataMasker() */ |
339 | | /* */ |
340 | | /* GDALMaskFunc for establishing a validity mask for a source */ |
341 | | /* band based on a provided NODATA value. */ |
342 | | /************************************************************************/ |
343 | | |
344 | | CPLErr GDALWarpNoDataMasker(void *pMaskFuncArg, int nBandCount, |
345 | | GDALDataType eType, int /* nXOff */, |
346 | | int /* nYOff */, int nXSize, int nYSize, |
347 | | GByte **ppImageData, int bMaskIsFloat, |
348 | | void *pValidityMask, int *pbOutAllValid) |
349 | | |
350 | 0 | { |
351 | 0 | const double *padfNoData = static_cast<double *>(pMaskFuncArg); |
352 | 0 | GUInt32 *panValidityMask = static_cast<GUInt32 *>(pValidityMask); |
353 | 0 | const size_t nPixels = static_cast<size_t>(nXSize) * nYSize; |
354 | |
|
355 | 0 | *pbOutAllValid = FALSE; |
356 | |
|
357 | 0 | if (nBandCount != 1 || bMaskIsFloat) |
358 | 0 | { |
359 | 0 | CPLError( |
360 | 0 | CE_Failure, CPLE_AppDefined, |
361 | 0 | "Invalid nBandCount or bMaskIsFloat argument in SourceNoDataMask"); |
362 | 0 | return CE_Failure; |
363 | 0 | } |
364 | | |
365 | 0 | CPLErr eErr = CE_None; |
366 | |
|
367 | 0 | switch (eType) |
368 | 0 | { |
369 | 0 | case GDT_Byte: |
370 | 0 | return GDALWarpNoDataMaskerT(padfNoData, nPixels, |
371 | 0 | *ppImageData, // Already a GByte *. |
372 | 0 | panValidityMask, pbOutAllValid); |
373 | | |
374 | 0 | case GDT_Int16: |
375 | 0 | return GDALWarpNoDataMaskerT( |
376 | 0 | padfNoData, nPixels, reinterpret_cast<GInt16 *>(*ppImageData), |
377 | 0 | panValidityMask, pbOutAllValid); |
378 | | |
379 | 0 | case GDT_UInt16: |
380 | 0 | return GDALWarpNoDataMaskerT( |
381 | 0 | padfNoData, nPixels, reinterpret_cast<GUInt16 *>(*ppImageData), |
382 | 0 | panValidityMask, pbOutAllValid); |
383 | | |
384 | 0 | case GDT_Float32: |
385 | 0 | { |
386 | 0 | const float fNoData = static_cast<float>(padfNoData[0]); |
387 | 0 | const float *pafData = reinterpret_cast<float *>(*ppImageData); |
388 | 0 | const bool bIsNoDataNan = CPL_TO_BOOL(std::isnan(fNoData)); |
389 | | |
390 | | // Nothing to do if value is out of range. |
391 | 0 | if (padfNoData[1] != 0.0) |
392 | 0 | { |
393 | 0 | *pbOutAllValid = TRUE; |
394 | 0 | return CE_None; |
395 | 0 | } |
396 | | |
397 | 0 | int bAllValid = TRUE; |
398 | 0 | for (size_t iOffset = 0; iOffset < nPixels; ++iOffset) |
399 | 0 | { |
400 | 0 | float fVal = pafData[iOffset]; |
401 | 0 | if ((bIsNoDataNan && std::isnan(fVal)) || |
402 | 0 | (!bIsNoDataNan && ARE_REAL_EQUAL(fVal, fNoData))) |
403 | 0 | { |
404 | 0 | bAllValid = FALSE; |
405 | 0 | CPLMaskClear(panValidityMask, iOffset); |
406 | 0 | } |
407 | 0 | } |
408 | 0 | *pbOutAllValid = bAllValid; |
409 | 0 | } |
410 | 0 | break; |
411 | | |
412 | 0 | case GDT_Float64: |
413 | 0 | { |
414 | 0 | const double dfNoData = padfNoData[0]; |
415 | 0 | const double *padfData = reinterpret_cast<double *>(*ppImageData); |
416 | 0 | const bool bIsNoDataNan = CPL_TO_BOOL(std::isnan(dfNoData)); |
417 | | |
418 | | // Nothing to do if value is out of range. |
419 | 0 | if (padfNoData[1] != 0.0) |
420 | 0 | { |
421 | 0 | *pbOutAllValid = TRUE; |
422 | 0 | return CE_None; |
423 | 0 | } |
424 | | |
425 | 0 | int bAllValid = TRUE; |
426 | 0 | for (size_t iOffset = 0; iOffset < nPixels; ++iOffset) |
427 | 0 | { |
428 | 0 | double dfVal = padfData[iOffset]; |
429 | 0 | if ((bIsNoDataNan && std::isnan(dfVal)) || |
430 | 0 | (!bIsNoDataNan && ARE_REAL_EQUAL(dfVal, dfNoData))) |
431 | 0 | { |
432 | 0 | bAllValid = FALSE; |
433 | 0 | CPLMaskClear(panValidityMask, iOffset); |
434 | 0 | } |
435 | 0 | } |
436 | 0 | *pbOutAllValid = bAllValid; |
437 | 0 | } |
438 | 0 | break; |
439 | | |
440 | 0 | default: |
441 | 0 | { |
442 | 0 | const int nWordSize = GDALGetDataTypeSizeBytes(eType); |
443 | |
|
444 | 0 | const bool bIsNoDataRealNan = |
445 | 0 | CPL_TO_BOOL(std::isnan(padfNoData[0])); |
446 | |
|
447 | 0 | eErr = CE_Failure; |
448 | 0 | double *padfWrk = static_cast<double *>( |
449 | 0 | VSI_MALLOC2_VERBOSE(nXSize, sizeof(double) * 2)); |
450 | 0 | if (padfWrk) |
451 | 0 | { |
452 | 0 | eErr = CE_None; |
453 | 0 | bool bAllValid = true; |
454 | 0 | for (int iLine = 0; iLine < nYSize; iLine++) |
455 | 0 | { |
456 | 0 | GDALCopyWords((*ppImageData) + nWordSize * iLine * nXSize, |
457 | 0 | eType, nWordSize, padfWrk, GDT_CFloat64, 16, |
458 | 0 | nXSize); |
459 | |
|
460 | 0 | for (int iPixel = 0; iPixel < nXSize; ++iPixel) |
461 | 0 | { |
462 | 0 | if (((bIsNoDataRealNan && |
463 | 0 | std::isnan(padfWrk[iPixel * 2])) || |
464 | 0 | (!bIsNoDataRealNan && |
465 | 0 | ARE_REAL_EQUAL(padfWrk[iPixel * 2], |
466 | 0 | padfNoData[0])))) |
467 | 0 | { |
468 | 0 | size_t iOffset = |
469 | 0 | iPixel + static_cast<size_t>(iLine) * nXSize; |
470 | |
|
471 | 0 | bAllValid = false; |
472 | 0 | CPLMaskClear(panValidityMask, iOffset); |
473 | 0 | } |
474 | 0 | } |
475 | 0 | } |
476 | 0 | *pbOutAllValid = bAllValid; |
477 | |
|
478 | 0 | VSIFree(padfWrk); |
479 | 0 | } |
480 | 0 | } |
481 | 0 | break; |
482 | 0 | } |
483 | | |
484 | 0 | return eErr; |
485 | 0 | } |
486 | | |
487 | | /************************************************************************/ |
488 | | /* GDALWarpSrcAlphaMasker() */ |
489 | | /* */ |
490 | | /* GDALMaskFunc for reading source simple 8bit alpha mask */ |
491 | | /* information and building a floating point density mask from */ |
492 | | /* it. */ |
493 | | /************************************************************************/ |
494 | | |
495 | | CPLErr GDALWarpSrcAlphaMasker(void *pMaskFuncArg, int /* nBandCount */, |
496 | | GDALDataType /* eType */, int nXOff, int nYOff, |
497 | | int nXSize, int nYSize, GByte ** /*ppImageData */, |
498 | | int bMaskIsFloat, void *pValidityMask, |
499 | | int *pbOutAllOpaque) |
500 | | |
501 | 0 | { |
502 | 0 | GDALWarpOptions *psWO = static_cast<GDALWarpOptions *>(pMaskFuncArg); |
503 | 0 | float *pafMask = static_cast<float *>(pValidityMask); |
504 | 0 | *pbOutAllOpaque = FALSE; |
505 | 0 | const size_t nPixels = static_cast<size_t>(nXSize) * nYSize; |
506 | | |
507 | | /* -------------------------------------------------------------------- */ |
508 | | /* Do some minimal checking. */ |
509 | | /* -------------------------------------------------------------------- */ |
510 | 0 | if (!bMaskIsFloat) |
511 | 0 | { |
512 | 0 | CPLAssert(false); |
513 | 0 | return CE_Failure; |
514 | 0 | } |
515 | | |
516 | 0 | if (psWO == nullptr || psWO->nSrcAlphaBand < 1) |
517 | 0 | { |
518 | 0 | CPLAssert(false); |
519 | 0 | return CE_Failure; |
520 | 0 | } |
521 | | |
522 | | /* -------------------------------------------------------------------- */ |
523 | | /* Read the alpha band. */ |
524 | | /* -------------------------------------------------------------------- */ |
525 | 0 | GDALRasterBandH hAlphaBand = |
526 | 0 | GDALGetRasterBand(psWO->hSrcDS, psWO->nSrcAlphaBand); |
527 | 0 | if (hAlphaBand == nullptr) |
528 | 0 | return CE_Failure; |
529 | | |
530 | | // Rescale. |
531 | 0 | const float inv_alpha_max = static_cast<float>( |
532 | 0 | 1.0 / CPLAtof(CSLFetchNameValueDef(psWO->papszWarpOptions, |
533 | 0 | "SRC_ALPHA_MAX", "255"))); |
534 | 0 | bool bOutAllOpaque = true; |
535 | |
|
536 | 0 | size_t iPixel = 0; |
537 | 0 | CPLErr eErr; |
538 | |
|
539 | 0 | #if (defined(__x86_64) || defined(_M_X64)) |
540 | 0 | GDALDataType eDT = GDALGetRasterDataType(hAlphaBand); |
541 | | // Make sure that pafMask is at least 8-byte aligned, which should |
542 | | // normally be always the case if being a ptr returned by malloc(). |
543 | 0 | if ((eDT == GDT_Byte || eDT == GDT_UInt16) && CPL_IS_ALIGNED(pafMask, 8)) |
544 | 0 | { |
545 | | // Read data. |
546 | 0 | eErr = GDALRasterIOEx( |
547 | 0 | hAlphaBand, GF_Read, nXOff, nYOff, nXSize, nYSize, pafMask, nXSize, |
548 | 0 | nYSize, eDT, static_cast<GSpacing>(sizeof(int)), |
549 | 0 | static_cast<GSpacing>(sizeof(int)) * nXSize, nullptr); |
550 | |
|
551 | 0 | if (eErr != CE_None) |
552 | 0 | return eErr; |
553 | | |
554 | | // Make sure we have the correct alignment before doing SSE |
555 | | // On Linux x86_64, the alignment should be always correct due |
556 | | // the alignment of malloc() being 16 byte. |
557 | 0 | const GUInt32 mask = (eDT == GDT_Byte) ? 0xff : 0xffff; |
558 | 0 | if (!CPL_IS_ALIGNED(pafMask, 16)) |
559 | 0 | { |
560 | 0 | pafMask[iPixel] = |
561 | 0 | (reinterpret_cast<GUInt32 *>(pafMask)[iPixel] & mask) * |
562 | 0 | inv_alpha_max; |
563 | 0 | if (pafMask[iPixel] >= 1.0f) |
564 | 0 | pafMask[iPixel] = 1.0f; |
565 | 0 | else |
566 | 0 | bOutAllOpaque = false; |
567 | 0 | iPixel++; |
568 | 0 | } |
569 | 0 | CPLAssert(CPL_IS_ALIGNED(pafMask + iPixel, 16)); |
570 | 0 | const __m128 xmm_inverse_alpha_max = _mm_load1_ps(&inv_alpha_max); |
571 | 0 | const float one_single = 1.0f; |
572 | 0 | const __m128 xmm_one = _mm_load1_ps(&one_single); |
573 | 0 | const __m128i xmm_i_mask = _mm_set1_epi32(mask); |
574 | 0 | __m128 xmmMaskNonOpaque0 = _mm_setzero_ps(); |
575 | 0 | __m128 xmmMaskNonOpaque1 = _mm_setzero_ps(); |
576 | 0 | __m128 xmmMaskNonOpaque2 = _mm_setzero_ps(); |
577 | 0 | for (; iPixel + 6 * 4 - 1 < nPixels; iPixel += 6 * 4) |
578 | 0 | { |
579 | 0 | __m128 xmm_mask0 = _mm_cvtepi32_ps(_mm_and_si128( |
580 | 0 | xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>( |
581 | 0 | pafMask + iPixel + 4 * 0)))); |
582 | 0 | __m128 xmm_mask1 = _mm_cvtepi32_ps(_mm_and_si128( |
583 | 0 | xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>( |
584 | 0 | pafMask + iPixel + 4 * 1)))); |
585 | 0 | __m128 xmm_mask2 = _mm_cvtepi32_ps(_mm_and_si128( |
586 | 0 | xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>( |
587 | 0 | pafMask + iPixel + 4 * 2)))); |
588 | 0 | __m128 xmm_mask3 = _mm_cvtepi32_ps(_mm_and_si128( |
589 | 0 | xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>( |
590 | 0 | pafMask + iPixel + 4 * 3)))); |
591 | 0 | __m128 xmm_mask4 = _mm_cvtepi32_ps(_mm_and_si128( |
592 | 0 | xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>( |
593 | 0 | pafMask + iPixel + 4 * 4)))); |
594 | 0 | __m128 xmm_mask5 = _mm_cvtepi32_ps(_mm_and_si128( |
595 | 0 | xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>( |
596 | 0 | pafMask + iPixel + 4 * 5)))); |
597 | 0 | xmm_mask0 = _mm_mul_ps(xmm_mask0, xmm_inverse_alpha_max); |
598 | 0 | xmm_mask1 = _mm_mul_ps(xmm_mask1, xmm_inverse_alpha_max); |
599 | 0 | xmm_mask2 = _mm_mul_ps(xmm_mask2, xmm_inverse_alpha_max); |
600 | 0 | xmm_mask3 = _mm_mul_ps(xmm_mask3, xmm_inverse_alpha_max); |
601 | 0 | xmm_mask4 = _mm_mul_ps(xmm_mask4, xmm_inverse_alpha_max); |
602 | 0 | xmm_mask5 = _mm_mul_ps(xmm_mask5, xmm_inverse_alpha_max); |
603 | 0 | xmmMaskNonOpaque0 = |
604 | 0 | _mm_or_ps(xmmMaskNonOpaque0, _mm_cmplt_ps(xmm_mask0, xmm_one)); |
605 | 0 | xmmMaskNonOpaque1 = |
606 | 0 | _mm_or_ps(xmmMaskNonOpaque1, _mm_cmplt_ps(xmm_mask1, xmm_one)); |
607 | 0 | xmmMaskNonOpaque2 = |
608 | 0 | _mm_or_ps(xmmMaskNonOpaque2, _mm_cmplt_ps(xmm_mask2, xmm_one)); |
609 | 0 | xmmMaskNonOpaque0 = |
610 | 0 | _mm_or_ps(xmmMaskNonOpaque0, _mm_cmplt_ps(xmm_mask3, xmm_one)); |
611 | 0 | xmmMaskNonOpaque1 = |
612 | 0 | _mm_or_ps(xmmMaskNonOpaque1, _mm_cmplt_ps(xmm_mask4, xmm_one)); |
613 | 0 | xmmMaskNonOpaque2 = |
614 | 0 | _mm_or_ps(xmmMaskNonOpaque2, _mm_cmplt_ps(xmm_mask5, xmm_one)); |
615 | 0 | xmm_mask0 = _mm_min_ps(xmm_mask0, xmm_one); |
616 | 0 | xmm_mask1 = _mm_min_ps(xmm_mask1, xmm_one); |
617 | 0 | xmm_mask2 = _mm_min_ps(xmm_mask2, xmm_one); |
618 | 0 | xmm_mask3 = _mm_min_ps(xmm_mask3, xmm_one); |
619 | 0 | xmm_mask4 = _mm_min_ps(xmm_mask4, xmm_one); |
620 | 0 | xmm_mask5 = _mm_min_ps(xmm_mask5, xmm_one); |
621 | 0 | _mm_store_ps(pafMask + iPixel + 4 * 0, xmm_mask0); |
622 | 0 | _mm_store_ps(pafMask + iPixel + 4 * 1, xmm_mask1); |
623 | 0 | _mm_store_ps(pafMask + iPixel + 4 * 2, xmm_mask2); |
624 | 0 | _mm_store_ps(pafMask + iPixel + 4 * 3, xmm_mask3); |
625 | 0 | _mm_store_ps(pafMask + iPixel + 4 * 4, xmm_mask4); |
626 | 0 | _mm_store_ps(pafMask + iPixel + 4 * 5, xmm_mask5); |
627 | 0 | } |
628 | 0 | if (_mm_movemask_ps( |
629 | 0 | _mm_or_ps(_mm_or_ps(xmmMaskNonOpaque0, xmmMaskNonOpaque1), |
630 | 0 | xmmMaskNonOpaque2))) |
631 | 0 | { |
632 | 0 | bOutAllOpaque = false; |
633 | 0 | } |
634 | 0 | for (; iPixel < nPixels; iPixel++) |
635 | 0 | { |
636 | 0 | pafMask[iPixel] = |
637 | 0 | (reinterpret_cast<GUInt32 *>(pafMask)[iPixel] & mask) * |
638 | 0 | inv_alpha_max; |
639 | 0 | if (pafMask[iPixel] >= 1.0f) |
640 | 0 | pafMask[iPixel] = 1.0f; |
641 | 0 | else |
642 | 0 | bOutAllOpaque = false; |
643 | 0 | } |
644 | 0 | } |
645 | 0 | else |
646 | 0 | #endif |
647 | 0 | { |
648 | | // Read data. |
649 | 0 | eErr = GDALRasterIO(hAlphaBand, GF_Read, nXOff, nYOff, nXSize, nYSize, |
650 | 0 | pafMask, nXSize, nYSize, GDT_Float32, 0, 0); |
651 | |
|
652 | 0 | if (eErr != CE_None) |
653 | 0 | return eErr; |
654 | | |
655 | | // TODO(rouault): Is loop unrolling by hand (r34564) actually helpful? |
656 | 0 | for (; iPixel + 3 < nPixels; iPixel += 4) |
657 | 0 | { |
658 | 0 | pafMask[iPixel] = pafMask[iPixel] * inv_alpha_max; |
659 | 0 | if (pafMask[iPixel] >= 1.0f) |
660 | 0 | pafMask[iPixel] = 1.0f; |
661 | 0 | else |
662 | 0 | bOutAllOpaque = false; |
663 | 0 | pafMask[iPixel + 1] = pafMask[iPixel + 1] * inv_alpha_max; |
664 | 0 | if (pafMask[iPixel + 1] >= 1.0f) |
665 | 0 | pafMask[iPixel + 1] = 1.0f; |
666 | 0 | else |
667 | 0 | bOutAllOpaque = false; |
668 | 0 | pafMask[iPixel + 2] = pafMask[iPixel + 2] * inv_alpha_max; |
669 | 0 | if (pafMask[iPixel + 2] >= 1.0f) |
670 | 0 | pafMask[iPixel + 2] = 1.0f; |
671 | 0 | else |
672 | 0 | bOutAllOpaque = false; |
673 | 0 | pafMask[iPixel + 3] = pafMask[iPixel + 3] * inv_alpha_max; |
674 | 0 | if (pafMask[iPixel + 3] >= 1.0f) |
675 | 0 | pafMask[iPixel + 3] = 1.0f; |
676 | 0 | else |
677 | 0 | bOutAllOpaque = false; |
678 | 0 | } |
679 | |
|
680 | 0 | for (; iPixel < nPixels; iPixel++) |
681 | 0 | { |
682 | 0 | pafMask[iPixel] = pafMask[iPixel] * inv_alpha_max; |
683 | 0 | if (pafMask[iPixel] >= 1.0f) |
684 | 0 | pafMask[iPixel] = 1.0f; |
685 | 0 | else |
686 | 0 | bOutAllOpaque = false; |
687 | 0 | } |
688 | 0 | } |
689 | | |
690 | 0 | *pbOutAllOpaque = bOutAllOpaque; |
691 | |
|
692 | 0 | return CE_None; |
693 | 0 | } |
694 | | |
695 | | /************************************************************************/ |
696 | | /* GDALWarpSrcMaskMasker() */ |
697 | | /* */ |
698 | | /* GDALMaskFunc for reading source simple 8bit validity mask */ |
699 | | /* information and building a one bit validity mask. */ |
700 | | /************************************************************************/ |
701 | | |
702 | | CPLErr GDALWarpSrcMaskMasker(void *pMaskFuncArg, int /* nBandCount */, |
703 | | GDALDataType /* eType */, int nXOff, int nYOff, |
704 | | int nXSize, int nYSize, GByte ** /*ppImageData */, |
705 | | int bMaskIsFloat, void *pValidityMask) |
706 | | |
707 | 0 | { |
708 | 0 | GDALWarpOptions *psWO = static_cast<GDALWarpOptions *>(pMaskFuncArg); |
709 | 0 | GUInt32 *panMask = static_cast<GUInt32 *>(pValidityMask); |
710 | | |
711 | | /* -------------------------------------------------------------------- */ |
712 | | /* Do some minimal checking. */ |
713 | | /* -------------------------------------------------------------------- */ |
714 | 0 | if (bMaskIsFloat) |
715 | 0 | { |
716 | 0 | CPLAssert(false); |
717 | 0 | return CE_Failure; |
718 | 0 | } |
719 | | |
720 | 0 | if (psWO == nullptr) |
721 | 0 | { |
722 | 0 | CPLAssert(false); |
723 | 0 | return CE_Failure; |
724 | 0 | } |
725 | | |
726 | | /* -------------------------------------------------------------------- */ |
727 | | /* Allocate a temporary buffer to read mask byte data into. */ |
728 | | /* -------------------------------------------------------------------- */ |
729 | 0 | GByte *pabySrcMask = |
730 | 0 | static_cast<GByte *>(VSI_MALLOC2_VERBOSE(nXSize, nYSize)); |
731 | 0 | if (pabySrcMask == nullptr) |
732 | 0 | { |
733 | 0 | return CE_Failure; |
734 | 0 | } |
735 | | |
736 | | /* -------------------------------------------------------------------- */ |
737 | | /* Fetch our mask band. */ |
738 | | /* -------------------------------------------------------------------- */ |
739 | 0 | GDALRasterBandH hMaskBand = nullptr; |
740 | 0 | GDALRasterBandH hSrcBand = |
741 | 0 | GDALGetRasterBand(psWO->hSrcDS, psWO->panSrcBands[0]); |
742 | 0 | if (hSrcBand != nullptr) |
743 | 0 | hMaskBand = GDALGetMaskBand(hSrcBand); |
744 | |
|
745 | 0 | if (hMaskBand == nullptr) |
746 | 0 | { |
747 | 0 | CPLAssert(false); |
748 | 0 | return CE_Failure; |
749 | 0 | } |
750 | | |
751 | | /* -------------------------------------------------------------------- */ |
752 | | /* Read the mask band. */ |
753 | | /* -------------------------------------------------------------------- */ |
754 | 0 | CPLErr eErr = GDALRasterIO(hMaskBand, GF_Read, nXOff, nYOff, nXSize, nYSize, |
755 | 0 | pabySrcMask, nXSize, nYSize, GDT_Byte, 0, 0); |
756 | |
|
757 | 0 | if (eErr != CE_None) |
758 | 0 | { |
759 | 0 | CPLFree(pabySrcMask); |
760 | 0 | return eErr; |
761 | 0 | } |
762 | | |
763 | | /* -------------------------------------------------------------------- */ |
764 | | /* Pack into 1 bit per pixel for validity. */ |
765 | | /* -------------------------------------------------------------------- */ |
766 | 0 | const size_t nPixels = static_cast<size_t>(nXSize) * nYSize; |
767 | 0 | for (size_t iPixel = 0; iPixel < nPixels; iPixel++) |
768 | 0 | { |
769 | 0 | if (pabySrcMask[iPixel] == 0) |
770 | 0 | CPLMaskClear(panMask, iPixel); |
771 | 0 | } |
772 | |
|
773 | 0 | CPLFree(pabySrcMask); |
774 | |
|
775 | 0 | return CE_None; |
776 | 0 | } |
777 | | |
778 | | /************************************************************************/ |
779 | | /* GDALWarpDstAlphaMasker() */ |
780 | | /* */ |
781 | | /* GDALMaskFunc for reading or writing the destination simple */ |
782 | | /* 8bit alpha mask information and building a floating point */ |
783 | | /* density mask from it. Note, writing is distinguished */ |
784 | | /* negative bandcount. */ |
785 | | /************************************************************************/ |
786 | | |
787 | | CPLErr GDALWarpDstAlphaMasker(void *pMaskFuncArg, int nBandCount, |
788 | | CPL_UNUSED GDALDataType /* eType */, int nXOff, |
789 | | int nYOff, int nXSize, int nYSize, |
790 | | GByte ** /*ppImageData */, int bMaskIsFloat, |
791 | | void *pValidityMask) |
792 | 0 | { |
793 | | /* -------------------------------------------------------------------- */ |
794 | | /* Do some minimal checking. */ |
795 | | /* -------------------------------------------------------------------- */ |
796 | 0 | if (!bMaskIsFloat) |
797 | 0 | { |
798 | 0 | CPLAssert(false); |
799 | 0 | return CE_Failure; |
800 | 0 | } |
801 | | |
802 | 0 | GDALWarpOptions *psWO = static_cast<GDALWarpOptions *>(pMaskFuncArg); |
803 | 0 | if (psWO == nullptr || psWO->nDstAlphaBand < 1) |
804 | 0 | { |
805 | 0 | CPLAssert(false); |
806 | 0 | return CE_Failure; |
807 | 0 | } |
808 | | |
809 | 0 | float *pafMask = static_cast<float *>(pValidityMask); |
810 | 0 | const size_t nPixels = static_cast<size_t>(nXSize) * nYSize; |
811 | |
|
812 | 0 | GDALRasterBandH hAlphaBand = |
813 | 0 | GDALGetRasterBand(psWO->hDstDS, psWO->nDstAlphaBand); |
814 | 0 | if (hAlphaBand == nullptr) |
815 | 0 | return CE_Failure; |
816 | | |
817 | 0 | size_t iPixel = 0; |
818 | | |
819 | | /* -------------------------------------------------------------------- */ |
820 | | /* Read alpha case. */ |
821 | | /* -------------------------------------------------------------------- */ |
822 | 0 | if (nBandCount >= 0) |
823 | 0 | { |
824 | 0 | const char *pszInitDest = |
825 | 0 | CSLFetchNameValue(psWO->papszWarpOptions, "INIT_DEST"); |
826 | | |
827 | | // Special logic for destinations being initialized on-the-fly. |
828 | 0 | if (pszInitDest != nullptr) |
829 | 0 | { |
830 | 0 | memset(pafMask, 0, nPixels * sizeof(float)); |
831 | 0 | return CE_None; |
832 | 0 | } |
833 | | |
834 | | // Rescale. |
835 | 0 | const float inv_alpha_max = static_cast<float>( |
836 | 0 | 1.0 / CPLAtof(CSLFetchNameValueDef(psWO->papszWarpOptions, |
837 | 0 | "DST_ALPHA_MAX", "255"))); |
838 | |
|
839 | 0 | #if (defined(__x86_64) || defined(_M_X64)) |
840 | 0 | const GDALDataType eDT = GDALGetRasterDataType(hAlphaBand); |
841 | | // Make sure that pafMask is at least 8-byte aligned, which should |
842 | | // normally be always the case if being a ptr returned by malloc(). |
843 | 0 | if ((eDT == GDT_Byte || eDT == GDT_UInt16) && |
844 | 0 | CPL_IS_ALIGNED(pafMask, 8)) |
845 | 0 | { |
846 | | // Read data. |
847 | 0 | const CPLErr eErr = GDALRasterIOEx( |
848 | 0 | hAlphaBand, GF_Read, nXOff, nYOff, nXSize, nYSize, pafMask, |
849 | 0 | nXSize, nYSize, eDT, static_cast<GSpacing>(sizeof(int)), |
850 | 0 | static_cast<GSpacing>(sizeof(int)) * nXSize, nullptr); |
851 | |
|
852 | 0 | if (eErr != CE_None) |
853 | 0 | return eErr; |
854 | | |
855 | | // Make sure we have the correct alignment before doing SSE |
856 | | // On Linux x86_64, the alignment should be always correct due |
857 | | // the alignment of malloc() being 16 byte. |
858 | 0 | const GUInt32 mask = (eDT == GDT_Byte) ? 0xff : 0xffff; |
859 | 0 | if (!CPL_IS_ALIGNED(pafMask, 16)) |
860 | 0 | { |
861 | 0 | pafMask[iPixel] = |
862 | 0 | (reinterpret_cast<GUInt32 *>(pafMask)[iPixel] & mask) * |
863 | 0 | inv_alpha_max; |
864 | 0 | pafMask[iPixel] = std::min(1.0f, pafMask[iPixel]); |
865 | 0 | iPixel++; |
866 | 0 | } |
867 | 0 | CPLAssert(CPL_IS_ALIGNED(pafMask + iPixel, 16)); |
868 | 0 | const __m128 xmm_inverse_alpha_max = _mm_load1_ps(&inv_alpha_max); |
869 | 0 | const float one_single = 1.0f; |
870 | 0 | const __m128 xmm_one = _mm_load1_ps(&one_single); |
871 | 0 | const __m128i xmm_i_mask = _mm_set1_epi32(mask); |
872 | 0 | for (; iPixel + 31 < nPixels; iPixel += 32) |
873 | 0 | { |
874 | 0 | __m128 xmm_mask0 = _mm_cvtepi32_ps(_mm_and_si128( |
875 | 0 | xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>( |
876 | 0 | pafMask + iPixel + 4 * 0)))); |
877 | 0 | __m128 xmm_mask1 = _mm_cvtepi32_ps(_mm_and_si128( |
878 | 0 | xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>( |
879 | 0 | pafMask + iPixel + 4 * 1)))); |
880 | 0 | __m128 xmm_mask2 = _mm_cvtepi32_ps(_mm_and_si128( |
881 | 0 | xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>( |
882 | 0 | pafMask + iPixel + 4 * 2)))); |
883 | 0 | __m128 xmm_mask3 = _mm_cvtepi32_ps(_mm_and_si128( |
884 | 0 | xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>( |
885 | 0 | pafMask + iPixel + 4 * 3)))); |
886 | 0 | __m128 xmm_mask4 = _mm_cvtepi32_ps(_mm_and_si128( |
887 | 0 | xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>( |
888 | 0 | pafMask + iPixel + 4 * 4)))); |
889 | 0 | __m128 xmm_mask5 = _mm_cvtepi32_ps(_mm_and_si128( |
890 | 0 | xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>( |
891 | 0 | pafMask + iPixel + 4 * 5)))); |
892 | 0 | __m128 xmm_mask6 = _mm_cvtepi32_ps(_mm_and_si128( |
893 | 0 | xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>( |
894 | 0 | pafMask + iPixel + 4 * 6)))); |
895 | 0 | __m128 xmm_mask7 = _mm_cvtepi32_ps(_mm_and_si128( |
896 | 0 | xmm_i_mask, _mm_load_si128(reinterpret_cast<__m128i *>( |
897 | 0 | pafMask + iPixel + 4 * 7)))); |
898 | 0 | xmm_mask0 = _mm_mul_ps(xmm_mask0, xmm_inverse_alpha_max); |
899 | 0 | xmm_mask1 = _mm_mul_ps(xmm_mask1, xmm_inverse_alpha_max); |
900 | 0 | xmm_mask2 = _mm_mul_ps(xmm_mask2, xmm_inverse_alpha_max); |
901 | 0 | xmm_mask3 = _mm_mul_ps(xmm_mask3, xmm_inverse_alpha_max); |
902 | 0 | xmm_mask4 = _mm_mul_ps(xmm_mask4, xmm_inverse_alpha_max); |
903 | 0 | xmm_mask5 = _mm_mul_ps(xmm_mask5, xmm_inverse_alpha_max); |
904 | 0 | xmm_mask6 = _mm_mul_ps(xmm_mask6, xmm_inverse_alpha_max); |
905 | 0 | xmm_mask7 = _mm_mul_ps(xmm_mask7, xmm_inverse_alpha_max); |
906 | 0 | xmm_mask0 = _mm_min_ps(xmm_mask0, xmm_one); |
907 | 0 | xmm_mask1 = _mm_min_ps(xmm_mask1, xmm_one); |
908 | 0 | xmm_mask2 = _mm_min_ps(xmm_mask2, xmm_one); |
909 | 0 | xmm_mask3 = _mm_min_ps(xmm_mask3, xmm_one); |
910 | 0 | xmm_mask4 = _mm_min_ps(xmm_mask4, xmm_one); |
911 | 0 | xmm_mask5 = _mm_min_ps(xmm_mask5, xmm_one); |
912 | 0 | xmm_mask6 = _mm_min_ps(xmm_mask6, xmm_one); |
913 | 0 | xmm_mask7 = _mm_min_ps(xmm_mask7, xmm_one); |
914 | 0 | _mm_store_ps(pafMask + iPixel + 4 * 0, xmm_mask0); |
915 | 0 | _mm_store_ps(pafMask + iPixel + 4 * 1, xmm_mask1); |
916 | 0 | _mm_store_ps(pafMask + iPixel + 4 * 2, xmm_mask2); |
917 | 0 | _mm_store_ps(pafMask + iPixel + 4 * 3, xmm_mask3); |
918 | 0 | _mm_store_ps(pafMask + iPixel + 4 * 4, xmm_mask4); |
919 | 0 | _mm_store_ps(pafMask + iPixel + 4 * 5, xmm_mask5); |
920 | 0 | _mm_store_ps(pafMask + iPixel + 4 * 6, xmm_mask6); |
921 | 0 | _mm_store_ps(pafMask + iPixel + 4 * 7, xmm_mask7); |
922 | 0 | } |
923 | 0 | for (; iPixel < nPixels; iPixel++) |
924 | 0 | { |
925 | 0 | pafMask[iPixel] = |
926 | 0 | (reinterpret_cast<GUInt32 *>(pafMask)[iPixel] & mask) * |
927 | 0 | inv_alpha_max; |
928 | 0 | pafMask[iPixel] = std::min(1.0f, pafMask[iPixel]); |
929 | 0 | } |
930 | 0 | } |
931 | 0 | else |
932 | 0 | #endif |
933 | 0 | { |
934 | | // Read data. |
935 | 0 | const CPLErr eErr = |
936 | 0 | GDALRasterIO(hAlphaBand, GF_Read, nXOff, nYOff, nXSize, nYSize, |
937 | 0 | pafMask, nXSize, nYSize, GDT_Float32, 0, 0); |
938 | |
|
939 | 0 | if (eErr != CE_None) |
940 | 0 | return eErr; |
941 | | |
942 | 0 | for (; iPixel < nPixels; iPixel++) |
943 | 0 | { |
944 | 0 | pafMask[iPixel] = pafMask[iPixel] * inv_alpha_max; |
945 | 0 | pafMask[iPixel] = std::min(1.0f, pafMask[iPixel]); |
946 | 0 | } |
947 | 0 | } |
948 | | |
949 | 0 | return CE_None; |
950 | 0 | } |
951 | | |
952 | | /* -------------------------------------------------------------------- */ |
953 | | /* Write alpha case. */ |
954 | | /* -------------------------------------------------------------------- */ |
955 | 0 | else |
956 | 0 | { |
957 | 0 | GDALDataType eDT = GDALGetRasterDataType(hAlphaBand); |
958 | 0 | const float cst_alpha_max = |
959 | 0 | static_cast<float>(CPLAtof(CSLFetchNameValueDef( |
960 | 0 | psWO->papszWarpOptions, "DST_ALPHA_MAX", "255"))) + |
961 | 0 | ((eDT == GDT_Byte || eDT == GDT_Int16 || eDT == GDT_UInt16 || |
962 | 0 | eDT == GDT_Int32 || eDT == GDT_UInt32) |
963 | 0 | ? 0.1f |
964 | 0 | : 0.0f); |
965 | |
|
966 | 0 | CPLErr eErr = CE_None; |
967 | |
|
968 | 0 | #if (defined(__x86_64) || defined(_M_X64)) |
969 | | // Make sure that pafMask is at least 8-byte aligned, which should |
970 | | // normally be always the case if being a ptr returned by malloc() |
971 | 0 | if ((eDT == GDT_Byte || eDT == GDT_Int16 || eDT == GDT_UInt16) && |
972 | 0 | CPL_IS_ALIGNED(pafMask, 8)) |
973 | 0 | { |
974 | | // Make sure we have the correct alignment before doing SSE |
975 | | // On Linux x86_64, the alignment should be always correct due |
976 | | // the alignment of malloc() being 16 byte |
977 | 0 | if (!CPL_IS_ALIGNED(pafMask, 16)) |
978 | 0 | { |
979 | 0 | reinterpret_cast<int *>(pafMask)[iPixel] = |
980 | 0 | static_cast<int>(pafMask[iPixel] * cst_alpha_max); |
981 | 0 | iPixel++; |
982 | 0 | } |
983 | 0 | CPLAssert(CPL_IS_ALIGNED(pafMask + iPixel, 16)); |
984 | 0 | const __m128 xmm_alpha_max = _mm_load1_ps(&cst_alpha_max); |
985 | 0 | for (; iPixel + 31 < nPixels; iPixel += 32) |
986 | 0 | { |
987 | 0 | __m128 xmm_mask0 = _mm_load_ps(pafMask + iPixel + 4 * 0); |
988 | 0 | __m128 xmm_mask1 = _mm_load_ps(pafMask + iPixel + 4 * 1); |
989 | 0 | __m128 xmm_mask2 = _mm_load_ps(pafMask + iPixel + 4 * 2); |
990 | 0 | __m128 xmm_mask3 = _mm_load_ps(pafMask + iPixel + 4 * 3); |
991 | 0 | __m128 xmm_mask4 = _mm_load_ps(pafMask + iPixel + 4 * 4); |
992 | 0 | __m128 xmm_mask5 = _mm_load_ps(pafMask + iPixel + 4 * 5); |
993 | 0 | __m128 xmm_mask6 = _mm_load_ps(pafMask + iPixel + 4 * 6); |
994 | 0 | __m128 xmm_mask7 = _mm_load_ps(pafMask + iPixel + 4 * 7); |
995 | 0 | xmm_mask0 = _mm_mul_ps(xmm_mask0, xmm_alpha_max); |
996 | 0 | xmm_mask1 = _mm_mul_ps(xmm_mask1, xmm_alpha_max); |
997 | 0 | xmm_mask2 = _mm_mul_ps(xmm_mask2, xmm_alpha_max); |
998 | 0 | xmm_mask3 = _mm_mul_ps(xmm_mask3, xmm_alpha_max); |
999 | 0 | xmm_mask4 = _mm_mul_ps(xmm_mask4, xmm_alpha_max); |
1000 | 0 | xmm_mask5 = _mm_mul_ps(xmm_mask5, xmm_alpha_max); |
1001 | 0 | xmm_mask6 = _mm_mul_ps(xmm_mask6, xmm_alpha_max); |
1002 | 0 | xmm_mask7 = _mm_mul_ps(xmm_mask7, xmm_alpha_max); |
1003 | | // Truncate to int. |
1004 | 0 | _mm_store_si128( |
1005 | 0 | reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 0), |
1006 | 0 | _mm_cvttps_epi32(xmm_mask0)); |
1007 | 0 | _mm_store_si128( |
1008 | 0 | reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 1), |
1009 | 0 | _mm_cvttps_epi32(xmm_mask1)); |
1010 | 0 | _mm_store_si128( |
1011 | 0 | reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 2), |
1012 | 0 | _mm_cvttps_epi32(xmm_mask2)); |
1013 | 0 | _mm_store_si128( |
1014 | 0 | reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 3), |
1015 | 0 | _mm_cvttps_epi32(xmm_mask3)); |
1016 | 0 | _mm_store_si128( |
1017 | 0 | reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 4), |
1018 | 0 | _mm_cvttps_epi32(xmm_mask4)); |
1019 | 0 | _mm_store_si128( |
1020 | 0 | reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 5), |
1021 | 0 | _mm_cvttps_epi32(xmm_mask5)); |
1022 | 0 | _mm_store_si128( |
1023 | 0 | reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 6), |
1024 | 0 | _mm_cvttps_epi32(xmm_mask6)); |
1025 | 0 | _mm_store_si128( |
1026 | 0 | reinterpret_cast<__m128i *>(pafMask + iPixel + 4 * 7), |
1027 | 0 | _mm_cvttps_epi32(xmm_mask7)); |
1028 | 0 | } |
1029 | 0 | for (; iPixel < nPixels; iPixel++) |
1030 | 0 | reinterpret_cast<int *>(pafMask)[iPixel] = |
1031 | 0 | static_cast<int>(pafMask[iPixel] * cst_alpha_max); |
1032 | | |
1033 | | // Write data. |
1034 | | // Assumes little endianness here. |
1035 | 0 | eErr = GDALRasterIOEx( |
1036 | 0 | hAlphaBand, GF_Write, nXOff, nYOff, nXSize, nYSize, pafMask, |
1037 | 0 | nXSize, nYSize, eDT, static_cast<GSpacing>(sizeof(int)), |
1038 | 0 | static_cast<GSpacing>(sizeof(int)) * nXSize, nullptr); |
1039 | 0 | } |
1040 | 0 | else |
1041 | 0 | #endif |
1042 | 0 | { |
1043 | 0 | for (; iPixel + 3 < nPixels; iPixel += 4) |
1044 | 0 | { |
1045 | 0 | pafMask[iPixel + 0] = static_cast<float>( |
1046 | 0 | static_cast<int>(pafMask[iPixel + 0] * cst_alpha_max)); |
1047 | 0 | pafMask[iPixel + 1] = static_cast<float>( |
1048 | 0 | static_cast<int>(pafMask[iPixel + 1] * cst_alpha_max)); |
1049 | 0 | pafMask[iPixel + 2] = static_cast<float>( |
1050 | 0 | static_cast<int>(pafMask[iPixel + 2] * cst_alpha_max)); |
1051 | 0 | pafMask[iPixel + 3] = static_cast<float>( |
1052 | 0 | static_cast<int>(pafMask[iPixel + 3] * cst_alpha_max)); |
1053 | 0 | } |
1054 | 0 | for (; iPixel < nPixels; iPixel++) |
1055 | 0 | pafMask[iPixel] = static_cast<float>( |
1056 | 0 | static_cast<int>(pafMask[iPixel] * cst_alpha_max)); |
1057 | | |
1058 | | // Write data. |
1059 | |
|
1060 | 0 | eErr = |
1061 | 0 | GDALRasterIO(hAlphaBand, GF_Write, nXOff, nYOff, nXSize, nYSize, |
1062 | 0 | pafMask, nXSize, nYSize, GDT_Float32, 0, 0); |
1063 | 0 | } |
1064 | 0 | return eErr; |
1065 | 0 | } |
1066 | 0 | } |
1067 | | |
1068 | | /************************************************************************/ |
1069 | | /* GDALWarpGetOptionList() */ |
1070 | | /************************************************************************/ |
1071 | | |
1072 | | /** Return a XML string describing options accepted by |
1073 | | * GDALWarpOptions::papszWarpOptions. |
1074 | | * |
1075 | | * @since 3.11 |
1076 | | */ |
1077 | | const char *GDALWarpGetOptionList(void) |
1078 | 0 | { |
1079 | 0 | return "<OptionList>" |
1080 | 0 | "<Option name='INIT_DEST' type='string' description='" |
1081 | 0 | "Numeric value or NO_DATA. This option forces the destination image " |
1082 | 0 | "to be initialized to the indicated value (for all bands) " |
1083 | 0 | "or indicates that it should be initialized to the NO_DATA value in " |
1084 | 0 | "padfDstNoDataReal/padfDstNoDataImag. If this value is not set the " |
1085 | 0 | "destination image will be read and overlaid.'/>" |
1086 | 0 | "<Option name='WRITE_FLUSH' type='boolean' description='" |
1087 | 0 | "This option forces a flush to disk of data after " |
1088 | 0 | "each chunk is processed. In some cases this helps ensure a serial " |
1089 | 0 | " writing of the output data otherwise a block of data may be " |
1090 | 0 | "written to disk each time a block of data is read for the input " |
1091 | 0 | "buffer resulting in a lot of extra seeking around the disk, and " |
1092 | 0 | "reduced IO throughput.' default='NO'/>" |
1093 | 0 | "<Option name='SKIP_NOSOURCE' type='boolean' description='" |
1094 | 0 | "Skip all processing for chunks for which there is no corresponding " |
1095 | 0 | "input data. This will disable initializing the destination " |
1096 | 0 | "(INIT_DEST) and all other processing, and so should be used " |
1097 | 0 | "carefully. Mostly useful to short circuit a lot of extra work " |
1098 | 0 | "in mosaicing situations. gdalwarp will automatically enable this " |
1099 | 0 | "option when it is assumed to be safe to do so.' default='NO'/>" |
1100 | | #ifdef undocumented |
1101 | | "<Option name='ERROR_OUT_IF_EMPTY_SOURCE_WINDOW' type='boolean' " |
1102 | | "description='By default, if the source window corresponding to the " |
1103 | | "current target window fails to be determined due to reprojection " |
1104 | | "errors, the warping fails. Setting this option to NO prevent such " |
1105 | | "failure from happening. The warped VRT mechanism automatically " |
1106 | | "sets it to NO.'/>" |
1107 | | #endif |
1108 | 0 | "<Option name='UNIFIED_SRC_NODATA' type='string-select' " |
1109 | 0 | "description='" |
1110 | 0 | "This setting determines how to take into account nodata values " |
1111 | 0 | "when there are several input bands. Consult " |
1112 | 0 | "GDALWarpOptions::papszWarpOptions documentation for more details.'>" |
1113 | 0 | " <Value>AUTO</Value>" |
1114 | 0 | " <Value>PARTIAL</Value>" |
1115 | 0 | " <Value>YES</Value>" |
1116 | 0 | " <Value>NO</Value>" |
1117 | 0 | "</Option>" |
1118 | 0 | "<Option name='CUTLINE' type='string' description='" |
1119 | 0 | "This may contain the WKT geometry for a cutline. It will be " |
1120 | 0 | "converted into a geometry by GDALWarpOperation::Initialize() and " |
1121 | 0 | "assigned to the GDALWarpOptions hCutline field. The coordinates " |
1122 | 0 | "must be expressed in source pixel/line coordinates. Note: this is " |
1123 | 0 | "different from the assumptions made for the -cutline option " |
1124 | 0 | "of the gdalwarp utility !'/>" |
1125 | 0 | "<Option name='CUTLINE_BLEND_DIST' type='float' description='" |
1126 | 0 | "This may be set with a distance in pixels which will be assigned " |
1127 | 0 | "to the dfCutlineBlendDist field in the GDALWarpOptions.'/>" |
1128 | 0 | "<Option name='CUTLINE_ALL_TOUCHED' type='boolean' description='" |
1129 | 0 | "This may be set to TRUE to enable ALL_TOUCHED mode when " |
1130 | 0 | "rasterizing cutline polygons. This is useful to ensure that that " |
1131 | 0 | "all pixels overlapping the cutline polygon will be selected, not " |
1132 | 0 | "just those whose center point falls within the polygon.' " |
1133 | 0 | "default='NO'/>" |
1134 | 0 | "<Option name='XSCALE' type='float' description='" |
1135 | 0 | "Ratio expressing the resampling factor (number of destination " |
1136 | 0 | "pixels per source pixel) along the target horizontal axis. The " |
1137 | 0 | "scale is used to determine the number of source pixels along the " |
1138 | 0 | "x-axis that are considered by the resampling algorithm. " |
1139 | 0 | "Equals to one for no resampling, below one for downsampling " |
1140 | 0 | "and above one for upsampling. This is automatically computed, " |
1141 | 0 | "for each processing chunk, and may thus vary among them, depending " |
1142 | 0 | "on the shape of output regions vs input regions. Such variations " |
1143 | 0 | "can be undesired in some situations. If the resampling factor " |
1144 | 0 | "can be considered as constant over the warped area, setting a " |
1145 | 0 | "constant value can lead to more reproducible pixel output.'/>" |
1146 | 0 | "<Option name='YSCALE' type='float' description='" |
1147 | 0 | "Same as XSCALE, but along the horizontal axis.'/>" |
1148 | 0 | "<Option name='OPTIMIZE_SIZE' type='boolean' description='" |
1149 | 0 | "This defaults to FALSE, but may be set to TRUE typically when " |
1150 | 0 | "writing to a compressed dataset (GeoTIFF with COMPRESS creation " |
1151 | 0 | "option set for example) for achieving a smaller file size. This " |
1152 | 0 | "is achieved by writing at once data aligned on full blocks of the " |
1153 | 0 | "target dataset, which avoids partial writes of compressed blocks " |
1154 | 0 | "and lost space when they are rewritten at the end of the file. " |
1155 | 0 | "However sticking to target block size may cause major processing " |
1156 | 0 | "slowdown for some particular reprojections. OPTIMIZE_SIZE mode " |
1157 | 0 | "is automatically enabled when it is safe to do so. " |
1158 | 0 | "As this parameter influences the shape of warping chunk, and by " |
1159 | 0 | "default the XSCALE and YSCALE parameters are computed per warping " |
1160 | 0 | "chunk, this parameter may influence the pixel output.' " |
1161 | 0 | "default='NO'/>" |
1162 | 0 | "<Option name='NUM_THREADS' type='string' description='" |
1163 | 0 | "Can be set to a numeric value or ALL_CPUS to set the number of " |
1164 | 0 | "threads to use to parallelize the computation part of the warping. " |
1165 | 0 | "If not set, computation will be done in a single thread..'/>" |
1166 | 0 | "<Option name='STREAMABLE_OUTPUT' type='boolean' description='" |
1167 | 0 | "This defaults to FALSE, but may be set to TRUE typically when " |
1168 | 0 | "writing to a streamed file. The gdalwarp utility automatically " |
1169 | 0 | "sets this option when writing to /vsistdout/ or a named pipe " |
1170 | 0 | "(on Unix). This option has performance impacts for some " |
1171 | 0 | "reprojections. Note: band interleaved output is " |
1172 | 0 | "not currently supported by the warping algorithm in a streamable " |
1173 | 0 | "compatible way.' default='NO'/>" |
1174 | 0 | "<Option name='SRC_COORD_PRECISION' type='float' description='" |
1175 | 0 | "Advanced setting. This defaults to 0, to indicate that no rounding " |
1176 | 0 | "of computing source image coordinates corresponding to the target " |
1177 | 0 | "image must be done. If greater than 0 (and typically below 1), " |
1178 | 0 | "this value, expressed in pixel, will be used to round computed " |
1179 | 0 | "source image coordinates. The purpose of this option is to make " |
1180 | 0 | "the results of warping with the approximated transformer more " |
1181 | 0 | "reproducible and not sensitive to changes in warping memory size. " |
1182 | 0 | "To achieve that, SRC_COORD_PRECISION must be at least 10 times " |
1183 | 0 | "greater than the error threshold. The higher the " |
1184 | 0 | "SRC_COORD_PRECISION/error_threshold ratio, the higher the " |
1185 | 0 | "performance will be, since exact reprojections must statistically " |
1186 | 0 | "be done with a frequency of " |
1187 | 0 | "4*error_threshold/SRC_COORD_PRECISION.' default='0'/>" |
1188 | 0 | "<Option name='SRC_ALPHA_MAX' type='float' description='" |
1189 | 0 | "Maximum value for the alpha band of the source dataset. If the " |
1190 | 0 | "value is not set and the alpha band has a NBITS metadata item, " |
1191 | 0 | "it is used to set SRC_ALPHA_MAX = 2^NBITS-1. Otherwise, if the " |
1192 | 0 | "value is not set and the alpha band is of type UInt16 " |
1193 | 0 | "(resp Int16), 65535 (resp 32767) is used. " |
1194 | 0 | "Otherwise, 255 is used.'/>" |
1195 | 0 | "<Option name='DST_ALPHA_MAX' type='float' description='" |
1196 | 0 | "Maximum value for the alpha band of the destination dataset. " |
1197 | 0 | "If the value is not set and the alpha band has a NBITS metadata " |
1198 | 0 | "item, it is used to set SRC_ALPHA_MAX = 2^NBITS-1. Otherwise, if " |
1199 | 0 | "the value is not set and the alpha band is of type UInt16 " |
1200 | 0 | "(resp Int16), 65535 (resp 32767) is used. " |
1201 | 0 | "Otherwise, 255 is used.'/>" |
1202 | 0 | "<Option name='SAMPLE_GRID' type='boolean' description='" |
1203 | 0 | "Setting this option to YES will force the sampling to " |
1204 | 0 | "include internal points as well as edge points which can be " |
1205 | 0 | "important if the transformation is esoteric inside out, or if " |
1206 | 0 | "large sections of the destination image are not transformable into " |
1207 | 0 | "the source coordinate system.' default='NO'/>" |
1208 | 0 | "<Option name='SAMPLE_STEPS' type='string' description='" |
1209 | 0 | "Modifies the density of the sampling grid. Increasing this can " |
1210 | 0 | "increase the computational cost, but improves the accuracy with " |
1211 | 0 | "which the source region is computed. This can be set to ALL to " |
1212 | 0 | "mean to sample along all edge points of the destination region " |
1213 | 0 | "(if SAMPLE_GRID=NO or not specified), or all points of the " |
1214 | 0 | "destination region if SAMPLE_GRID=YES.' default='21'/>" |
1215 | 0 | "<Option name='SOURCE_EXTRA' type='int' description='" |
1216 | 0 | "This is a number of extra pixels added around the source " |
1217 | 0 | "window for a given request, and by default it is 1 to take care " |
1218 | 0 | "of rounding error. Setting this larger will increase the amount of " |
1219 | 0 | "data that needs to be read, but can avoid missing source data.' " |
1220 | 0 | "default='1'/>" |
1221 | 0 | "<Option name='APPLY_VERTICAL_SHIFT' type='boolean' description='" |
1222 | 0 | "Force the use of vertical shift. This option is generally not " |
1223 | 0 | "necessary, except when using an explicit coordinate transformation " |
1224 | 0 | "(COORDINATE_OPERATION), and not specifying an explicit source and " |
1225 | 0 | "target SRS.'/>" |
1226 | 0 | "<Option name='MULT_FACTOR_VERTICAL_SHIFT' type='float' " |
1227 | 0 | "description='" |
1228 | 0 | "Multiplication factor for the vertical shift' default='1.0'/>" |
1229 | 0 | "<Option name='EXCLUDED_VALUES' type='string' " |
1230 | 0 | "description='" |
1231 | 0 | "Comma-separated tuple of values (thus typically \"R,G,B\"), that " |
1232 | 0 | "are ignored as contributing source pixels during resampling. " |
1233 | 0 | "The number of values in the tuple must be the same as the number " |
1234 | 0 | "of bands, excluding the alpha band. Several tuples of excluded " |
1235 | 0 | "values may be specified using the \"(R1,G1,B2),(R2,G2,B2)\" syntax." |
1236 | 0 | " Only taken into account by Average currently. This concept is a " |
1237 | 0 | "bit similar to nodata/alpha, but the main difference is that " |
1238 | 0 | "pixels matching one of the excluded value tuples are still " |
1239 | 0 | "considered as valid, when determining the target pixel " |
1240 | 0 | "validity/density.'/>" |
1241 | 0 | "<Option name='EXCLUDED_VALUES_PCT_THRESHOLD' type='float' " |
1242 | 0 | "min='0' max='100' description='" |
1243 | 0 | "Minimum percentage of source pixels that must be set at one of " |
1244 | 0 | "the EXCLUDED_VALUES to cause the excluded value, that is in " |
1245 | 0 | "majority among source pixels, to be used as the target pixel " |
1246 | 0 | "value. Only taken into account by Average currently.' " |
1247 | 0 | "default='50'/>" |
1248 | 0 | "<Option name='NODATA_VALUES_PCT_THRESHOLD' type='float' " |
1249 | 0 | "min='0' max='100' description='" |
1250 | 0 | "Minimum percentage of source pixels that must be at nodata (or " |
1251 | 0 | "alpha=0 or any other way to express transparent pixel) to cause " |
1252 | 0 | "the target pixel value to not be set. Default value is 100 (%), " |
1253 | 0 | "which means that a target pixel is not set only if all " |
1254 | 0 | "contributing source pixels are not set. Note that " |
1255 | 0 | "NODATA_VALUES_PCT_THRESHOLD is taken into account before " |
1256 | 0 | "EXCLUDED_VALUES_PCT_THRESHOLD. Only taken into account by Average " |
1257 | 0 | "currently.' default='100'/>" |
1258 | 0 | "<Option name='MODE_TIES' type='string-select' " |
1259 | 0 | "description='" |
1260 | 0 | "Strategy to use when breaking ties with MODE resampling. " |
1261 | 0 | "By default, the first value encountered will be used. " |
1262 | 0 | "Alternatively, the minimum or maximum value can be selected.' " |
1263 | 0 | "default='FIRST'>" |
1264 | 0 | " <Value>FIRST</Value>" |
1265 | 0 | " <Value>MIN</Value>" |
1266 | 0 | " <Value>MAX</Value>" |
1267 | 0 | "</Option>" |
1268 | 0 | "</OptionList>"; |
1269 | 0 | } |
1270 | | |
1271 | | /************************************************************************/ |
1272 | | /* ==================================================================== */ |
1273 | | /* GDALWarpOptions */ |
1274 | | /* ==================================================================== */ |
1275 | | /************************************************************************/ |
1276 | | |
1277 | | /** |
1278 | | * \var char **GDALWarpOptions::papszWarpOptions; |
1279 | | * |
1280 | | * A string list of additional options controlling the warp operation in |
1281 | | * name=value format. A suitable string list can be prepared with |
1282 | | * CSLSetNameValue(). |
1283 | | * |
1284 | | * The available options can also be retrieved programmatically with |
1285 | | * GDALWarpGetOptionList(). |
1286 | | * |
1287 | | * The following values are currently supported: |
1288 | | * <ul> |
1289 | | * <li>INIT_DEST=[value] or INIT_DEST=NO_DATA: This option forces the |
1290 | | * destination image to be initialized to the indicated value (for all bands) |
1291 | | * or indicates that it should be initialized to the NO_DATA value in |
1292 | | * padfDstNoDataReal/padfDstNoDataImag. If this value isn't set the |
1293 | | * destination image will be read and overlaid.</li> |
1294 | | * |
1295 | | * <li>WRITE_FLUSH=YES/NO: This option forces a flush to disk of data after |
1296 | | * each chunk is processed. In some cases this helps ensure a serial |
1297 | | * writing of the output data otherwise a block of data may be written to disk |
1298 | | * each time a block of data is read for the input buffer resulting in a lot |
1299 | | * of extra seeking around the disk, and reduced IO throughput. The default |
1300 | | * is NO.</li> |
1301 | | * |
1302 | | * <li>SKIP_NOSOURCE=YES/NO: Skip all processing for chunks for which there |
1303 | | * is no corresponding input data. This will disable initializing the |
1304 | | * destination (INIT_DEST) and all other processing, and so should be used |
1305 | | * carefully. Mostly useful to short circuit a lot of extra work in mosaicing |
1306 | | * situations. Starting with GDAL 2.4, gdalwarp will automatically enable this |
1307 | | * option when it is assumed to be safe to do so.</li> |
1308 | | * |
1309 | | * <li>UNIFIED_SRC_NODATA=YES/NO/PARTIAL: This setting determines |
1310 | | * how to take into account nodata values when there are several input bands. |
1311 | | * <ul> |
1312 | | * <li>When YES, all bands are considered as nodata if and only if, all bands |
1313 | | * match the corresponding nodata values. |
1314 | | * Note: UNIFIED_SRC_NODATA=YES is set by default, when called from gdalwarp |
1315 | | * / GDALWarp() with an explicit -srcnodata setting. |
1316 | | * |
1317 | | * Example with nodata values at (1, 2, 3) and target alpha band requested. |
1318 | | * <ul> |
1319 | | * <li>input pixel = (1, 2, 3) ==> output pixel = (0, 0, 0, 0)</li> |
1320 | | * <li>input pixel = (1, 2, 127) ==> output pixel = (1, 2, 127, 255)</li> |
1321 | | * </ul> |
1322 | | * </li> |
1323 | | * <li>When NO, nodata masking values is considered independently for each band. |
1324 | | * A potential target alpha band will always be valid if there are multiple |
1325 | | * bands. |
1326 | | * |
1327 | | * Example with nodata values at (1, 2, 3) and target alpha band requested. |
1328 | | * <ul> |
1329 | | * <li>input pixel = (1, 2, 3) ==> output pixel = (0, 0, 0, 255)</li> |
1330 | | * <li>input pixel = (1, 2, 127) ==> output pixel = (0, 0, 127, 255)</li> |
1331 | | * </ul> |
1332 | | * |
1333 | | * Note: NO was the default behavior before GDAL 3.3.2 |
1334 | | * </li> |
1335 | | * <li>When PARTIAL, or not specified at all (default behavior), |
1336 | | * nodata masking values is considered independently for each band. |
1337 | | * But, and this is the difference with NO, if for a given pixel, it |
1338 | | * evaluates to the nodata value of each band, the target pixel is |
1339 | | * considered as globally invalid, which impacts the value of a potential |
1340 | | * target alpha band. |
1341 | | * |
1342 | | * Note: PARTIAL is new to GDAL 3.3.2 and should not be used with |
1343 | | * earlier versions. The default behavior of GDAL < 3.3.2 was NO. |
1344 | | * |
1345 | | * Example with nodata values at (1, 2, 3) and target alpha band requested. |
1346 | | * <ul> |
1347 | | * <li>input pixel = (1, 2, 3) ==> output pixel = (0, 0, 0, 0)</li> |
1348 | | * <li>input pixel = (1, 2, 127) ==> output pixel = (0, 0, 127, 255)</li> |
1349 | | * </ul> |
1350 | | * </li> |
1351 | | * </ul> |
1352 | | * </li> |
1353 | | * |
1354 | | * <li>CUTLINE: This may contain the WKT geometry for a cutline. It will |
1355 | | * be converted into a geometry by GDALWarpOperation::Initialize() and assigned |
1356 | | * to the GDALWarpOptions hCutline field. The coordinates must be expressed |
1357 | | * in source pixel/line coordinates. Note: this is different from the |
1358 | | * assumptions made for the -cutline option of the gdalwarp utility !</li> |
1359 | | * |
1360 | | * <li>CUTLINE_BLEND_DIST: This may be set with a distance in pixels which |
1361 | | * will be assigned to the dfCutlineBlendDist field in the GDALWarpOptions.</li> |
1362 | | * |
1363 | | * <li>CUTLINE_ALL_TOUCHED: This defaults to FALSE, but may be set to TRUE |
1364 | | * to enable ALL_TOUCHED mode when rasterizing cutline polygons. This is |
1365 | | * useful to ensure that that all pixels overlapping the cutline polygon |
1366 | | * will be selected, not just those whose center point falls within the |
1367 | | * polygon.</li> |
1368 | | * |
1369 | | * <li>XSCALE: Ratio expressing the resampling factor (number of destination |
1370 | | * pixels per source pixel) along the target horizontal axis. |
1371 | | * The scale is used to determine the number of source pixels along the x-axis |
1372 | | * that are considered by the resampling algorithm. |
1373 | | * Equals to one for no resampling, below one for downsampling |
1374 | | * and above one for upsampling. This is automatically computed, for each |
1375 | | * processing chunk, and may thus vary among them, depending on the |
1376 | | * shape of output regions vs input regions. Such variations can be undesired |
1377 | | * in some situations. If the resampling factor can be considered as constant |
1378 | | * over the warped area, setting a constant value can lead to more reproducible |
1379 | | * pixel output.</li> |
1380 | | * |
1381 | | * <li>YSCALE: Same as XSCALE, but along the horizontal axis.</li> |
1382 | | * |
1383 | | * <li>OPTIMIZE_SIZE: This defaults to FALSE, but may be set to TRUE |
1384 | | * typically when writing to a compressed dataset (GeoTIFF with |
1385 | | * COMPRESS creation option set for example) for achieving a smaller |
1386 | | * file size. This is achieved by writing at once data aligned on full |
1387 | | * blocks of the target dataset, which avoids partial writes of |
1388 | | * compressed blocks and lost space when they are rewritten at the end |
1389 | | * of the file. However sticking to target block size may cause major |
1390 | | * processing slowdown for some particular reprojections. Starting |
1391 | | * with GDAL 3.8, OPTIMIZE_SIZE mode is automatically enabled when it is safe |
1392 | | * to do so. |
1393 | | * As this parameter influences the shape of warping chunk, and by default the |
1394 | | * XSCALE and YSCALE parameters are computed per warping chunk, this parameter may |
1395 | | * influence the pixel output. |
1396 | | * </li> |
1397 | | * |
1398 | | * <li>NUM_THREADS: (GDAL >= 1.10) Can be set to a numeric value or ALL_CPUS to |
1399 | | * set the number of threads to use to parallelize the computation part of the |
1400 | | * warping. If not set, computation will be done in a single thread.</li> |
1401 | | * |
1402 | | * <li>STREAMABLE_OUTPUT: (GDAL >= 2.0) This defaults to FALSE, but may |
1403 | | * be set to TRUE typically when writing to a streamed file. The |
1404 | | * gdalwarp utility automatically sets this option when writing to |
1405 | | * /vsistdout/ or a named pipe (on Unix). This option has performance |
1406 | | * impacts for some reprojections. Note: band interleaved output is |
1407 | | * not currently supported by the warping algorithm in a streamable |
1408 | | * compatible way.</li> |
1409 | | * |
1410 | | * <li>SRC_COORD_PRECISION: (GDAL >= 2.0). Advanced setting. This |
1411 | | * defaults to 0, to indicate that no rounding of computing source |
1412 | | * image coordinates corresponding to the target image must be |
1413 | | * done. If greater than 0 (and typically below 1), this value, |
1414 | | * expressed in pixel, will be used to round computed source image |
1415 | | * coordinates. The purpose of this option is to make the results of |
1416 | | * warping with the approximated transformer more reproducible and not |
1417 | | * sensitive to changes in warping memory size. To achieve that, |
1418 | | * SRC_COORD_PRECISION must be at least 10 times greater than the |
1419 | | * error threshold. The higher the SRC_COORD_PRECISION/error_threshold |
1420 | | * ratio, the higher the performance will be, since exact |
1421 | | * reprojections must statistically be done with a frequency of |
1422 | | * 4*error_threshold/SRC_COORD_PRECISION.</li> |
1423 | | * |
1424 | | * <li>SRC_ALPHA_MAX: (GDAL >= 2.2). Maximum value for the alpha band of the |
1425 | | * source dataset. If the value is not set and the alpha band has a NBITS |
1426 | | * metadata item, it is used to set SRC_ALPHA_MAX = 2^NBITS-1. Otherwise, if the |
1427 | | * value is not set and the alpha band is of type UInt16 (resp Int16), 65535 |
1428 | | * (resp 32767) is used. Otherwise, 255 is used.</li> |
1429 | | * |
1430 | | * <li>DST_ALPHA_MAX: (GDAL >= 2.2). Maximum value for the alpha band of the |
1431 | | * destination dataset. If the value is not set and the alpha band has a NBITS |
1432 | | * metadata item, it is used to set DST_ALPHA_MAX = 2^NBITS-1. Otherwise, if the |
1433 | | * value is not set and the alpha band is of type UInt16 (resp Int16), 65535 |
1434 | | * (resp 32767) is used. Otherwise, 255 is used.</li> |
1435 | | * </ul> |
1436 | | * |
1437 | | * Normally when computing the source raster data to |
1438 | | * load to generate a particular output area, the warper samples transforms |
1439 | | * 21 points along each edge of the destination region back onto the source |
1440 | | * file, and uses this to compute a bounding window on the source image that |
1441 | | * is sufficient. Depending on the transformation in effect, the source |
1442 | | * window may be a bit too small, or even missing large areas. Problem |
1443 | | * situations are those where the transformation is very non-linear or |
1444 | | * "inside out". Examples are transforming from WGS84 to Polar Stereographic |
1445 | | * for areas around the pole, or transformations where some of the image is |
1446 | | * untransformable. The following options provide some additional control |
1447 | | * to deal with errors in computing the source window: |
1448 | | * <ul> |
1449 | | * |
1450 | | * <li>SAMPLE_GRID=YES/NO: Setting this option to YES will force the sampling to |
1451 | | * include internal points as well as edge points which can be important if |
1452 | | * the transformation is esoteric inside out, or if large sections of the |
1453 | | * destination image are not transformable into the source coordinate |
1454 | | * system.</li> |
1455 | | * |
1456 | | * <li>SAMPLE_STEPS: Modifies the density of the sampling grid. The default |
1457 | | * number of steps is 21. Increasing this can increase the computational |
1458 | | * cost, but improves the accuracy with which the source region is |
1459 | | * computed. |
1460 | | * Starting with GDAL 3.7, this can be set to ALL to mean to sample |
1461 | | * along all edge points of the destination region (if SAMPLE_GRID=NO or not |
1462 | | * specified), or all points of the destination region if SAMPLE_GRID=YES.</li> |
1463 | | * |
1464 | | * <li>SOURCE_EXTRA: This is a number of extra pixels added around the source |
1465 | | * window for a given request, and by default it is 1 to take care of rounding |
1466 | | * error. Setting this larger will increase the amount of data that needs to |
1467 | | * be read, but can avoid missing source data.</li> |
1468 | | * <li>APPLY_VERTICAL_SHIFT=YES/NO: Force the use of vertical shift. |
1469 | | * This option is generally not necessary, except when using an explicit |
1470 | | * coordinate transformation (COORDINATE_OPERATION), and not specifying |
1471 | | * an explicit source and target SRS.</li> |
1472 | | * <li>MULT_FACTOR_VERTICAL_SHIFT: Multiplication factor for the vertical |
1473 | | * shift. Default 1.0</li> |
1474 | | * |
1475 | | * <li>EXCLUDED_VALUES: (GDAL >= 3.9) Comma-separated tuple of values |
1476 | | * (thus typically "R,G,B"), that are ignored as contributing source |
1477 | | * pixels during resampling. The number of values in the tuple must be the same |
1478 | | * as the number of bands, excluding the alpha band. |
1479 | | * Several tuples of excluded values may be specified using the |
1480 | | * "(R1,G1,B2),(R2,G2,B2)" syntax. |
1481 | | * Only taken into account by Average currently. |
1482 | | * This concept is a bit similar to nodata/alpha, but the main difference is |
1483 | | * that pixels matching one of the excluded value tuples are still considered |
1484 | | * as valid, when determining the target pixel validity/density. |
1485 | | * </li> |
1486 | | * |
1487 | | * <li>EXCLUDED_VALUES_PCT_THRESHOLD=[0-100]: (GDAL >= 3.9) Minimum percentage |
1488 | | * of source pixels that must be set at one of the EXCLUDED_VALUES to cause |
1489 | | * the excluded value, that is in majority among source pixels, to be used as the |
1490 | | * target pixel value. Default value is 50 (%). |
1491 | | * Only taken into account by Average currently.</li> |
1492 | | * |
1493 | | * <li>NODATA_VALUES_PCT_THRESHOLD=[0-100]: (GDAL >= 3.9) Minimum percentage |
1494 | | * of source pixels that must be at nodata (or alpha=0 or any other way to express |
1495 | | * transparent pixel) to cause the target pixel value to not be set. Default |
1496 | | * value is 100 (%), which means that a target pixel is not set only if all |
1497 | | * contributing source pixels are not set. |
1498 | | * Note that NODATA_VALUES_PCT_THRESHOLD is taken into account before |
1499 | | * EXCLUDED_VALUES_PCT_THRESHOLD. |
1500 | | * Only taken into account by Average currently.</li> |
1501 | | * |
1502 | | * <li>MODE_TIES=FIRST/MIN/MAX: (GDAL >= 3.11) Strategy to use when breaking |
1503 | | * ties with MODE resampling. By default, the first value encountered will be used. |
1504 | | * Alternatively, the minimum or maximum value can be selected.</li> |
1505 | | * |
1506 | | * </ul> |
1507 | | */ |
1508 | | |
1509 | | /************************************************************************/ |
1510 | | /* GDALCreateWarpOptions() */ |
1511 | | /************************************************************************/ |
1512 | | |
1513 | | /** Create a warp options structure. |
1514 | | * |
1515 | | * Must be deallocated with GDALDestroyWarpOptions() |
1516 | | */ |
1517 | | GDALWarpOptions *CPL_STDCALL GDALCreateWarpOptions() |
1518 | | |
1519 | 0 | { |
1520 | 0 | GDALWarpOptions *psOptions = |
1521 | 0 | static_cast<GDALWarpOptions *>(CPLCalloc(sizeof(GDALWarpOptions), 1)); |
1522 | |
|
1523 | 0 | psOptions->nBandCount = 0; |
1524 | 0 | psOptions->eResampleAlg = GRA_NearestNeighbour; |
1525 | 0 | psOptions->pfnProgress = GDALDummyProgress; |
1526 | 0 | psOptions->eWorkingDataType = GDT_Unknown; |
1527 | 0 | psOptions->eTieStrategy = GWKTS_First; |
1528 | |
|
1529 | 0 | return psOptions; |
1530 | 0 | } |
1531 | | |
1532 | | /************************************************************************/ |
1533 | | /* GDALDestroyWarpOptions() */ |
1534 | | /************************************************************************/ |
1535 | | |
1536 | | /** Destroy a warp options structure. */ |
1537 | | void CPL_STDCALL GDALDestroyWarpOptions(GDALWarpOptions *psOptions) |
1538 | | |
1539 | 0 | { |
1540 | 0 | if (psOptions == nullptr) |
1541 | 0 | return; |
1542 | | |
1543 | 0 | CSLDestroy(psOptions->papszWarpOptions); |
1544 | 0 | CPLFree(psOptions->panSrcBands); |
1545 | 0 | CPLFree(psOptions->panDstBands); |
1546 | 0 | CPLFree(psOptions->padfSrcNoDataReal); |
1547 | 0 | CPLFree(psOptions->padfSrcNoDataImag); |
1548 | 0 | CPLFree(psOptions->padfDstNoDataReal); |
1549 | 0 | CPLFree(psOptions->padfDstNoDataImag); |
1550 | 0 | CPLFree(psOptions->papfnSrcPerBandValidityMaskFunc); |
1551 | 0 | CPLFree(psOptions->papSrcPerBandValidityMaskFuncArg); |
1552 | |
|
1553 | 0 | if (psOptions->hCutline != nullptr) |
1554 | 0 | delete static_cast<OGRGeometry *>(psOptions->hCutline); |
1555 | |
|
1556 | 0 | CPLFree(psOptions); |
1557 | 0 | } |
1558 | | |
1559 | | #define COPY_MEM(target, type, count) \ |
1560 | 0 | do \ |
1561 | 0 | { \ |
1562 | 0 | if ((psSrcOptions->target) != nullptr && (count) != 0) \ |
1563 | 0 | { \ |
1564 | 0 | (psDstOptions->target) = \ |
1565 | 0 | static_cast<type *>(CPLMalloc(sizeof(type) * (count))); \ |
1566 | 0 | memcpy((psDstOptions->target), (psSrcOptions->target), \ |
1567 | 0 | sizeof(type) * (count)); \ |
1568 | 0 | } \ |
1569 | 0 | else \ |
1570 | 0 | (psDstOptions->target) = nullptr; \ |
1571 | 0 | } while (false) |
1572 | | |
1573 | | /************************************************************************/ |
1574 | | /* GDALCloneWarpOptions() */ |
1575 | | /************************************************************************/ |
1576 | | |
1577 | | /** Clone a warp options structure. |
1578 | | * |
1579 | | * Must be deallocated with GDALDestroyWarpOptions() |
1580 | | */ |
1581 | | GDALWarpOptions *CPL_STDCALL |
1582 | | GDALCloneWarpOptions(const GDALWarpOptions *psSrcOptions) |
1583 | | |
1584 | 0 | { |
1585 | 0 | GDALWarpOptions *psDstOptions = GDALCreateWarpOptions(); |
1586 | |
|
1587 | 0 | memcpy(psDstOptions, psSrcOptions, sizeof(GDALWarpOptions)); |
1588 | |
|
1589 | 0 | if (psSrcOptions->papszWarpOptions != nullptr) |
1590 | 0 | psDstOptions->papszWarpOptions = |
1591 | 0 | CSLDuplicate(psSrcOptions->papszWarpOptions); |
1592 | |
|
1593 | 0 | COPY_MEM(panSrcBands, int, psSrcOptions->nBandCount); |
1594 | 0 | COPY_MEM(panDstBands, int, psSrcOptions->nBandCount); |
1595 | 0 | COPY_MEM(padfSrcNoDataReal, double, psSrcOptions->nBandCount); |
1596 | 0 | COPY_MEM(padfSrcNoDataImag, double, psSrcOptions->nBandCount); |
1597 | 0 | COPY_MEM(padfDstNoDataReal, double, psSrcOptions->nBandCount); |
1598 | 0 | COPY_MEM(padfDstNoDataImag, double, psSrcOptions->nBandCount); |
1599 | | // cppcheck-suppress pointerSize |
1600 | 0 | COPY_MEM(papfnSrcPerBandValidityMaskFunc, GDALMaskFunc, |
1601 | 0 | psSrcOptions->nBandCount); |
1602 | 0 | psDstOptions->papSrcPerBandValidityMaskFuncArg = nullptr; |
1603 | |
|
1604 | 0 | if (psSrcOptions->hCutline != nullptr) |
1605 | 0 | psDstOptions->hCutline = |
1606 | 0 | OGR_G_Clone(static_cast<OGRGeometryH>(psSrcOptions->hCutline)); |
1607 | 0 | psDstOptions->dfCutlineBlendDist = psSrcOptions->dfCutlineBlendDist; |
1608 | |
|
1609 | 0 | return psDstOptions; |
1610 | 0 | } |
1611 | | |
1612 | | namespace |
1613 | | { |
1614 | | void InitNoData(int nBandCount, double **ppdNoDataReal, double dDataReal) |
1615 | 0 | { |
1616 | 0 | if (nBandCount <= 0) |
1617 | 0 | { |
1618 | 0 | return; |
1619 | 0 | } |
1620 | 0 | if (*ppdNoDataReal != nullptr) |
1621 | 0 | { |
1622 | 0 | return; |
1623 | 0 | } |
1624 | | |
1625 | 0 | *ppdNoDataReal = |
1626 | 0 | static_cast<double *>(CPLMalloc(sizeof(double) * nBandCount)); |
1627 | |
|
1628 | 0 | for (int i = 0; i < nBandCount; ++i) |
1629 | 0 | { |
1630 | 0 | (*ppdNoDataReal)[i] = dDataReal; |
1631 | 0 | } |
1632 | 0 | } |
1633 | | } // namespace |
1634 | | |
1635 | | /************************************************************************/ |
1636 | | /* GDALWarpInitDstNoDataReal() */ |
1637 | | /************************************************************************/ |
1638 | | |
1639 | | /** |
1640 | | * \brief Initialize padfDstNoDataReal with specified value. |
1641 | | * |
1642 | | * @param psOptionsIn options to initialize. |
1643 | | * @param dNoDataReal value to initialize to. |
1644 | | * |
1645 | | */ |
1646 | | void CPL_STDCALL GDALWarpInitDstNoDataReal(GDALWarpOptions *psOptionsIn, |
1647 | | double dNoDataReal) |
1648 | 0 | { |
1649 | 0 | VALIDATE_POINTER0(psOptionsIn, "GDALWarpInitDstNoDataReal"); |
1650 | 0 | InitNoData(psOptionsIn->nBandCount, &psOptionsIn->padfDstNoDataReal, |
1651 | 0 | dNoDataReal); |
1652 | 0 | } |
1653 | | |
1654 | | /************************************************************************/ |
1655 | | /* GDALWarpInitSrcNoDataReal() */ |
1656 | | /************************************************************************/ |
1657 | | |
1658 | | /** |
1659 | | * \brief Initialize padfSrcNoDataReal with specified value. |
1660 | | * |
1661 | | * @param psOptionsIn options to initialize. |
1662 | | * @param dNoDataReal value to initialize to. |
1663 | | * |
1664 | | */ |
1665 | | void CPL_STDCALL GDALWarpInitSrcNoDataReal(GDALWarpOptions *psOptionsIn, |
1666 | | double dNoDataReal) |
1667 | 0 | { |
1668 | 0 | VALIDATE_POINTER0(psOptionsIn, "GDALWarpInitSrcNoDataReal"); |
1669 | 0 | InitNoData(psOptionsIn->nBandCount, &psOptionsIn->padfSrcNoDataReal, |
1670 | 0 | dNoDataReal); |
1671 | 0 | } |
1672 | | |
1673 | | /************************************************************************/ |
1674 | | /* GDALWarpInitNoDataReal() */ |
1675 | | /************************************************************************/ |
1676 | | |
1677 | | /** |
1678 | | * \brief Initialize padfSrcNoDataReal and padfDstNoDataReal with specified |
1679 | | * value. |
1680 | | * |
1681 | | * @param psOptionsIn options to initialize. |
1682 | | * @param dNoDataReal value to initialize to. |
1683 | | * |
1684 | | */ |
1685 | | void CPL_STDCALL GDALWarpInitNoDataReal(GDALWarpOptions *psOptionsIn, |
1686 | | double dNoDataReal) |
1687 | 0 | { |
1688 | 0 | GDALWarpInitDstNoDataReal(psOptionsIn, dNoDataReal); |
1689 | 0 | GDALWarpInitSrcNoDataReal(psOptionsIn, dNoDataReal); |
1690 | 0 | } |
1691 | | |
1692 | | /************************************************************************/ |
1693 | | /* GDALWarpInitDstNoDataImag() */ |
1694 | | /************************************************************************/ |
1695 | | |
1696 | | /** |
1697 | | * \brief Initialize padfDstNoDataImag with specified value. |
1698 | | * |
1699 | | * @param psOptionsIn options to initialize. |
1700 | | * @param dNoDataImag value to initialize to. |
1701 | | * |
1702 | | */ |
1703 | | void CPL_STDCALL GDALWarpInitDstNoDataImag(GDALWarpOptions *psOptionsIn, |
1704 | | double dNoDataImag) |
1705 | 0 | { |
1706 | 0 | VALIDATE_POINTER0(psOptionsIn, "GDALWarpInitDstNoDataImag"); |
1707 | 0 | InitNoData(psOptionsIn->nBandCount, &psOptionsIn->padfDstNoDataImag, |
1708 | 0 | dNoDataImag); |
1709 | 0 | } |
1710 | | |
1711 | | /************************************************************************/ |
1712 | | /* GDALWarpInitSrcNoDataImag() */ |
1713 | | /************************************************************************/ |
1714 | | |
1715 | | /** |
1716 | | * \brief Initialize padfSrcNoDataImag with specified value. |
1717 | | * |
1718 | | * @param psOptionsIn options to initialize. |
1719 | | * @param dNoDataImag value to initialize to. |
1720 | | * |
1721 | | */ |
1722 | | void CPL_STDCALL GDALWarpInitSrcNoDataImag(GDALWarpOptions *psOptionsIn, |
1723 | | double dNoDataImag) |
1724 | 0 | { |
1725 | 0 | VALIDATE_POINTER0(psOptionsIn, "GDALWarpInitSrcNoDataImag"); |
1726 | 0 | InitNoData(psOptionsIn->nBandCount, &psOptionsIn->padfSrcNoDataImag, |
1727 | 0 | dNoDataImag); |
1728 | 0 | } |
1729 | | |
1730 | | /************************************************************************/ |
1731 | | /* GDALWarpResolveWorkingDataType() */ |
1732 | | /************************************************************************/ |
1733 | | |
1734 | | /** |
1735 | | * \brief If the working data type is unknown, this method will determine |
1736 | | * a valid working data type to support the data in the src and dest |
1737 | | * data sets and any noData values. |
1738 | | * |
1739 | | * @param psOptions options to initialize. |
1740 | | * |
1741 | | */ |
1742 | | void CPL_STDCALL GDALWarpResolveWorkingDataType(GDALWarpOptions *psOptions) |
1743 | 0 | { |
1744 | 0 | if (psOptions == nullptr) |
1745 | 0 | { |
1746 | 0 | return; |
1747 | 0 | } |
1748 | | /* -------------------------------------------------------------------- */ |
1749 | | /* If no working data type was provided, set one now. */ |
1750 | | /* */ |
1751 | | /* Ensure that the working data type can encapsulate any value */ |
1752 | | /* in the target, source, and the no data for either. */ |
1753 | | /* -------------------------------------------------------------------- */ |
1754 | 0 | if (psOptions->eWorkingDataType != GDT_Unknown) |
1755 | 0 | { |
1756 | 0 | return; |
1757 | 0 | } |
1758 | | |
1759 | 0 | psOptions->eWorkingDataType = GDT_Byte; |
1760 | | |
1761 | | // If none of the provided input nodata values can be represented in the |
1762 | | // data type of the corresponding source band, ignore them. |
1763 | 0 | if (psOptions->hSrcDS && psOptions->padfSrcNoDataReal) |
1764 | 0 | { |
1765 | 0 | int nCountInvalidSrcNoDataReal = 0; |
1766 | 0 | for (int iBand = 0; iBand < psOptions->nBandCount; iBand++) |
1767 | 0 | { |
1768 | 0 | GDALRasterBandH hSrcBand = GDALGetRasterBand( |
1769 | 0 | psOptions->hSrcDS, psOptions->panSrcBands[iBand]); |
1770 | |
|
1771 | 0 | if (hSrcBand && |
1772 | 0 | !GDALIsValueExactAs(psOptions->padfSrcNoDataReal[iBand], |
1773 | 0 | GDALGetRasterDataType(hSrcBand))) |
1774 | 0 | { |
1775 | 0 | nCountInvalidSrcNoDataReal++; |
1776 | 0 | } |
1777 | 0 | } |
1778 | 0 | if (nCountInvalidSrcNoDataReal == psOptions->nBandCount) |
1779 | 0 | { |
1780 | 0 | CPLFree(psOptions->padfSrcNoDataReal); |
1781 | 0 | psOptions->padfSrcNoDataReal = nullptr; |
1782 | 0 | CPLFree(psOptions->padfSrcNoDataImag); |
1783 | 0 | psOptions->padfSrcNoDataImag = nullptr; |
1784 | 0 | } |
1785 | 0 | } |
1786 | |
|
1787 | 0 | for (int iBand = 0; iBand < psOptions->nBandCount; iBand++) |
1788 | 0 | { |
1789 | 0 | if (psOptions->hDstDS != nullptr) |
1790 | 0 | { |
1791 | 0 | GDALRasterBandH hDstBand = GDALGetRasterBand( |
1792 | 0 | psOptions->hDstDS, psOptions->panDstBands[iBand]); |
1793 | |
|
1794 | 0 | if (hDstBand != nullptr) |
1795 | 0 | { |
1796 | 0 | psOptions->eWorkingDataType = |
1797 | 0 | GDALDataTypeUnion(psOptions->eWorkingDataType, |
1798 | 0 | GDALGetRasterDataType(hDstBand)); |
1799 | 0 | } |
1800 | 0 | } |
1801 | |
|
1802 | 0 | if (psOptions->hSrcDS != nullptr) |
1803 | 0 | { |
1804 | 0 | GDALRasterBandH hSrcBand = GDALGetRasterBand( |
1805 | 0 | psOptions->hSrcDS, psOptions->panSrcBands[iBand]); |
1806 | |
|
1807 | 0 | if (hSrcBand != nullptr) |
1808 | 0 | { |
1809 | 0 | psOptions->eWorkingDataType = |
1810 | 0 | GDALDataTypeUnion(psOptions->eWorkingDataType, |
1811 | 0 | GDALGetRasterDataType(hSrcBand)); |
1812 | 0 | } |
1813 | 0 | } |
1814 | |
|
1815 | 0 | if (psOptions->padfSrcNoDataReal != nullptr) |
1816 | 0 | { |
1817 | 0 | psOptions->eWorkingDataType = GDALDataTypeUnionWithValue( |
1818 | 0 | psOptions->eWorkingDataType, |
1819 | 0 | psOptions->padfSrcNoDataReal[iBand], false); |
1820 | 0 | } |
1821 | |
|
1822 | 0 | if (psOptions->padfSrcNoDataImag != nullptr && |
1823 | 0 | psOptions->padfSrcNoDataImag[iBand] != 0.0) |
1824 | 0 | { |
1825 | 0 | psOptions->eWorkingDataType = GDALDataTypeUnionWithValue( |
1826 | 0 | psOptions->eWorkingDataType, |
1827 | 0 | psOptions->padfSrcNoDataImag[iBand], true); |
1828 | 0 | } |
1829 | |
|
1830 | 0 | if (psOptions->padfDstNoDataReal != nullptr) |
1831 | 0 | { |
1832 | 0 | psOptions->eWorkingDataType = GDALDataTypeUnionWithValue( |
1833 | 0 | psOptions->eWorkingDataType, |
1834 | 0 | psOptions->padfDstNoDataReal[iBand], false); |
1835 | 0 | } |
1836 | |
|
1837 | 0 | if (psOptions->padfDstNoDataImag != nullptr && |
1838 | 0 | psOptions->padfDstNoDataImag[iBand] != 0.0) |
1839 | 0 | { |
1840 | 0 | psOptions->eWorkingDataType = GDALDataTypeUnionWithValue( |
1841 | 0 | psOptions->eWorkingDataType, |
1842 | 0 | psOptions->padfDstNoDataImag[iBand], true); |
1843 | 0 | } |
1844 | 0 | } |
1845 | |
|
1846 | 0 | const bool bApplyVerticalShift = CPLFetchBool( |
1847 | 0 | psOptions->papszWarpOptions, "APPLY_VERTICAL_SHIFT", false); |
1848 | 0 | if (bApplyVerticalShift && |
1849 | 0 | GDALDataTypeIsInteger(psOptions->eWorkingDataType)) |
1850 | 0 | { |
1851 | 0 | const double dfMultFactorVerticalShift = CPLAtof(CSLFetchNameValueDef( |
1852 | 0 | psOptions->papszWarpOptions, "MULT_FACTOR_VERTICAL_SHIFT", "1.0")); |
1853 | 0 | if (dfMultFactorVerticalShift != 1) |
1854 | 0 | { |
1855 | 0 | psOptions->eWorkingDataType = |
1856 | 0 | GDALDataTypeUnion(psOptions->eWorkingDataType, GDT_Float32); |
1857 | 0 | } |
1858 | 0 | } |
1859 | 0 | } |
1860 | | |
1861 | | /************************************************************************/ |
1862 | | /* GDALWarpInitDefaultBandMapping() */ |
1863 | | /************************************************************************/ |
1864 | | |
1865 | | /** |
1866 | | * \brief Init src and dst band mappings such that Bands[i] = i+1 |
1867 | | * for nBandCount |
1868 | | * Does nothing if psOptionsIn->nBandCount is non-zero. |
1869 | | * |
1870 | | * @param psOptionsIn options to initialize. |
1871 | | * @param nBandCount bands to initialize for. |
1872 | | * |
1873 | | */ |
1874 | | void CPL_STDCALL GDALWarpInitDefaultBandMapping(GDALWarpOptions *psOptionsIn, |
1875 | | int nBandCount) |
1876 | 0 | { |
1877 | 0 | if (psOptionsIn->nBandCount != 0) |
1878 | 0 | { |
1879 | 0 | return; |
1880 | 0 | } |
1881 | | |
1882 | 0 | psOptionsIn->nBandCount = nBandCount; |
1883 | |
|
1884 | 0 | psOptionsIn->panSrcBands = |
1885 | 0 | static_cast<int *>(CPLMalloc(sizeof(int) * psOptionsIn->nBandCount)); |
1886 | 0 | psOptionsIn->panDstBands = |
1887 | 0 | static_cast<int *>(CPLMalloc(sizeof(int) * psOptionsIn->nBandCount)); |
1888 | |
|
1889 | 0 | for (int i = 0; i < psOptionsIn->nBandCount; i++) |
1890 | 0 | { |
1891 | 0 | psOptionsIn->panSrcBands[i] = i + 1; |
1892 | 0 | psOptionsIn->panDstBands[i] = i + 1; |
1893 | 0 | } |
1894 | 0 | } |
1895 | | |
1896 | | /************************************************************************/ |
1897 | | /* GDALSerializeWarpOptions() */ |
1898 | | /************************************************************************/ |
1899 | | |
1900 | | CPLXMLNode *CPL_STDCALL GDALSerializeWarpOptions(const GDALWarpOptions *psWO) |
1901 | | |
1902 | 0 | { |
1903 | | /* -------------------------------------------------------------------- */ |
1904 | | /* Create root. */ |
1905 | | /* -------------------------------------------------------------------- */ |
1906 | 0 | CPLXMLNode *psTree = |
1907 | 0 | CPLCreateXMLNode(nullptr, CXT_Element, "GDALWarpOptions"); |
1908 | | |
1909 | | /* -------------------------------------------------------------------- */ |
1910 | | /* WarpMemoryLimit */ |
1911 | | /* -------------------------------------------------------------------- */ |
1912 | 0 | CPLCreateXMLElementAndValue( |
1913 | 0 | psTree, "WarpMemoryLimit", |
1914 | 0 | CPLString().Printf("%g", psWO->dfWarpMemoryLimit)); |
1915 | | |
1916 | | /* -------------------------------------------------------------------- */ |
1917 | | /* ResampleAlg */ |
1918 | | /* -------------------------------------------------------------------- */ |
1919 | 0 | const char *pszAlgName = nullptr; |
1920 | |
|
1921 | 0 | if (psWO->eResampleAlg == GRA_NearestNeighbour) |
1922 | 0 | pszAlgName = "NearestNeighbour"; |
1923 | 0 | else if (psWO->eResampleAlg == GRA_Bilinear) |
1924 | 0 | pszAlgName = "Bilinear"; |
1925 | 0 | else if (psWO->eResampleAlg == GRA_Cubic) |
1926 | 0 | pszAlgName = "Cubic"; |
1927 | 0 | else if (psWO->eResampleAlg == GRA_CubicSpline) |
1928 | 0 | pszAlgName = "CubicSpline"; |
1929 | 0 | else if (psWO->eResampleAlg == GRA_Lanczos) |
1930 | 0 | pszAlgName = "Lanczos"; |
1931 | 0 | else if (psWO->eResampleAlg == GRA_Average) |
1932 | 0 | pszAlgName = "Average"; |
1933 | 0 | else if (psWO->eResampleAlg == GRA_RMS) |
1934 | 0 | pszAlgName = "RootMeanSquare"; |
1935 | 0 | else if (psWO->eResampleAlg == GRA_Mode) |
1936 | 0 | pszAlgName = "Mode"; |
1937 | 0 | else if (psWO->eResampleAlg == GRA_Max) |
1938 | 0 | pszAlgName = "Maximum"; |
1939 | 0 | else if (psWO->eResampleAlg == GRA_Min) |
1940 | 0 | pszAlgName = "Minimum"; |
1941 | 0 | else if (psWO->eResampleAlg == GRA_Med) |
1942 | 0 | pszAlgName = "Median"; |
1943 | 0 | else if (psWO->eResampleAlg == GRA_Q1) |
1944 | 0 | pszAlgName = "Quartile1"; |
1945 | 0 | else if (psWO->eResampleAlg == GRA_Q3) |
1946 | 0 | pszAlgName = "Quartile3"; |
1947 | 0 | else if (psWO->eResampleAlg == GRA_Sum) |
1948 | 0 | pszAlgName = "Sum"; |
1949 | 0 | else |
1950 | 0 | pszAlgName = "Unknown"; |
1951 | |
|
1952 | 0 | CPLCreateXMLElementAndValue(psTree, "ResampleAlg", pszAlgName); |
1953 | | |
1954 | | /* -------------------------------------------------------------------- */ |
1955 | | /* Working Data Type */ |
1956 | | /* -------------------------------------------------------------------- */ |
1957 | 0 | CPLCreateXMLElementAndValue(psTree, "WorkingDataType", |
1958 | 0 | GDALGetDataTypeName(psWO->eWorkingDataType)); |
1959 | | |
1960 | | /* -------------------------------------------------------------------- */ |
1961 | | /* Name/value warp options. */ |
1962 | | /* -------------------------------------------------------------------- */ |
1963 | 0 | for (int iWO = 0; psWO->papszWarpOptions != nullptr && |
1964 | 0 | psWO->papszWarpOptions[iWO] != nullptr; |
1965 | 0 | iWO++) |
1966 | 0 | { |
1967 | 0 | char *pszName = nullptr; |
1968 | 0 | const char *pszValue = |
1969 | 0 | CPLParseNameValue(psWO->papszWarpOptions[iWO], &pszName); |
1970 | | |
1971 | | // EXTRA_ELTS is an internal detail that we will recover |
1972 | | // no need to serialize it. |
1973 | | // And CUTLINE is also serialized in a special way |
1974 | 0 | if (pszName != nullptr && !EQUAL(pszName, "EXTRA_ELTS") && |
1975 | 0 | !EQUAL(pszName, "CUTLINE")) |
1976 | 0 | { |
1977 | 0 | CPLXMLNode *psOption = |
1978 | 0 | CPLCreateXMLElementAndValue(psTree, "Option", pszValue); |
1979 | |
|
1980 | 0 | CPLCreateXMLNode(CPLCreateXMLNode(psOption, CXT_Attribute, "name"), |
1981 | 0 | CXT_Text, pszName); |
1982 | 0 | } |
1983 | |
|
1984 | 0 | CPLFree(pszName); |
1985 | 0 | } |
1986 | | |
1987 | | /* -------------------------------------------------------------------- */ |
1988 | | /* Source and Destination Data Source */ |
1989 | | /* -------------------------------------------------------------------- */ |
1990 | 0 | if (psWO->hSrcDS != nullptr) |
1991 | 0 | { |
1992 | 0 | CPLCreateXMLElementAndValue(psTree, "SourceDataset", |
1993 | 0 | GDALGetDescription(psWO->hSrcDS)); |
1994 | |
|
1995 | 0 | CSLConstList papszOpenOptions = |
1996 | 0 | GDALDataset::FromHandle(psWO->hSrcDS)->GetOpenOptions(); |
1997 | 0 | GDALSerializeOpenOptionsToXML(psTree, papszOpenOptions); |
1998 | 0 | } |
1999 | |
|
2000 | 0 | if (psWO->hDstDS != nullptr && |
2001 | 0 | strlen(GDALGetDescription(psWO->hDstDS)) != 0) |
2002 | 0 | { |
2003 | 0 | CPLCreateXMLElementAndValue(psTree, "DestinationDataset", |
2004 | 0 | GDALGetDescription(psWO->hDstDS)); |
2005 | 0 | } |
2006 | | |
2007 | | /* -------------------------------------------------------------------- */ |
2008 | | /* Serialize transformer. */ |
2009 | | /* -------------------------------------------------------------------- */ |
2010 | 0 | if (psWO->pfnTransformer != nullptr) |
2011 | 0 | { |
2012 | 0 | CPLXMLNode *psTransformerContainer = |
2013 | 0 | CPLCreateXMLNode(psTree, CXT_Element, "Transformer"); |
2014 | |
|
2015 | 0 | CPLXMLNode *psTransformerTree = GDALSerializeTransformer( |
2016 | 0 | psWO->pfnTransformer, psWO->pTransformerArg); |
2017 | |
|
2018 | 0 | if (psTransformerTree != nullptr) |
2019 | 0 | CPLAddXMLChild(psTransformerContainer, psTransformerTree); |
2020 | 0 | } |
2021 | | |
2022 | | /* -------------------------------------------------------------------- */ |
2023 | | /* Band count and lists. */ |
2024 | | /* -------------------------------------------------------------------- */ |
2025 | 0 | CPLXMLNode *psBandList = nullptr; |
2026 | |
|
2027 | 0 | if (psWO->nBandCount != 0) |
2028 | 0 | psBandList = CPLCreateXMLNode(psTree, CXT_Element, "BandList"); |
2029 | |
|
2030 | 0 | for (int i = 0; i < psWO->nBandCount; i++) |
2031 | 0 | { |
2032 | 0 | CPLXMLNode *psBand; |
2033 | |
|
2034 | 0 | psBand = CPLCreateXMLNode(psBandList, CXT_Element, "BandMapping"); |
2035 | 0 | if (psWO->panSrcBands != nullptr) |
2036 | 0 | CPLCreateXMLNode(CPLCreateXMLNode(psBand, CXT_Attribute, "src"), |
2037 | 0 | CXT_Text, |
2038 | 0 | CPLString().Printf("%d", psWO->panSrcBands[i])); |
2039 | 0 | if (psWO->panDstBands != nullptr) |
2040 | 0 | CPLCreateXMLNode(CPLCreateXMLNode(psBand, CXT_Attribute, "dst"), |
2041 | 0 | CXT_Text, |
2042 | 0 | CPLString().Printf("%d", psWO->panDstBands[i])); |
2043 | |
|
2044 | 0 | if (psWO->padfSrcNoDataReal != nullptr) |
2045 | 0 | { |
2046 | 0 | CPLCreateXMLElementAndValue( |
2047 | 0 | psBand, "SrcNoDataReal", |
2048 | 0 | VRTSerializeNoData(psWO->padfSrcNoDataReal[i], |
2049 | 0 | psWO->eWorkingDataType, 16) |
2050 | 0 | .c_str()); |
2051 | 0 | } |
2052 | |
|
2053 | 0 | if (psWO->padfSrcNoDataImag != nullptr) |
2054 | 0 | { |
2055 | 0 | if (std::isnan(psWO->padfSrcNoDataImag[i])) |
2056 | 0 | CPLCreateXMLElementAndValue(psBand, "SrcNoDataImag", "nan"); |
2057 | 0 | else |
2058 | 0 | CPLCreateXMLElementAndValue( |
2059 | 0 | psBand, "SrcNoDataImag", |
2060 | 0 | CPLString().Printf("%.16g", psWO->padfSrcNoDataImag[i])); |
2061 | 0 | } |
2062 | | // Compatibility with GDAL <= 2.2: if we serialize a SrcNoDataReal, |
2063 | | // it needs a SrcNoDataImag as well |
2064 | 0 | else if (psWO->padfSrcNoDataReal != nullptr) |
2065 | 0 | { |
2066 | 0 | CPLCreateXMLElementAndValue(psBand, "SrcNoDataImag", "0"); |
2067 | 0 | } |
2068 | |
|
2069 | 0 | if (psWO->padfDstNoDataReal != nullptr) |
2070 | 0 | { |
2071 | 0 | CPLCreateXMLElementAndValue( |
2072 | 0 | psBand, "DstNoDataReal", |
2073 | 0 | VRTSerializeNoData(psWO->padfDstNoDataReal[i], |
2074 | 0 | psWO->eWorkingDataType, 16) |
2075 | 0 | .c_str()); |
2076 | 0 | } |
2077 | |
|
2078 | 0 | if (psWO->padfDstNoDataImag != nullptr) |
2079 | 0 | { |
2080 | 0 | if (std::isnan(psWO->padfDstNoDataImag[i])) |
2081 | 0 | CPLCreateXMLElementAndValue(psBand, "DstNoDataImag", "nan"); |
2082 | 0 | else |
2083 | 0 | CPLCreateXMLElementAndValue( |
2084 | 0 | psBand, "DstNoDataImag", |
2085 | 0 | CPLString().Printf("%.16g", psWO->padfDstNoDataImag[i])); |
2086 | 0 | } |
2087 | | // Compatibility with GDAL <= 2.2: if we serialize a DstNoDataReal, |
2088 | | // it needs a SrcNoDataImag as well |
2089 | 0 | else if (psWO->padfDstNoDataReal != nullptr) |
2090 | 0 | { |
2091 | 0 | CPLCreateXMLElementAndValue(psBand, "DstNoDataImag", "0"); |
2092 | 0 | } |
2093 | 0 | } |
2094 | | |
2095 | | /* -------------------------------------------------------------------- */ |
2096 | | /* Alpha bands. */ |
2097 | | /* -------------------------------------------------------------------- */ |
2098 | 0 | if (psWO->nSrcAlphaBand > 0) |
2099 | 0 | CPLCreateXMLElementAndValue( |
2100 | 0 | psTree, "SrcAlphaBand", |
2101 | 0 | CPLString().Printf("%d", psWO->nSrcAlphaBand)); |
2102 | |
|
2103 | 0 | if (psWO->nDstAlphaBand > 0) |
2104 | 0 | CPLCreateXMLElementAndValue( |
2105 | 0 | psTree, "DstAlphaBand", |
2106 | 0 | CPLString().Printf("%d", psWO->nDstAlphaBand)); |
2107 | | |
2108 | | /* -------------------------------------------------------------------- */ |
2109 | | /* Cutline. */ |
2110 | | /* -------------------------------------------------------------------- */ |
2111 | 0 | if (psWO->hCutline != nullptr) |
2112 | 0 | { |
2113 | 0 | char *pszWKT = nullptr; |
2114 | 0 | if (OGR_G_ExportToWkt(static_cast<OGRGeometryH>(psWO->hCutline), |
2115 | 0 | &pszWKT) == OGRERR_NONE) |
2116 | 0 | { |
2117 | 0 | CPLCreateXMLElementAndValue(psTree, "Cutline", pszWKT); |
2118 | 0 | } |
2119 | 0 | CPLFree(pszWKT); |
2120 | 0 | } |
2121 | |
|
2122 | 0 | if (psWO->dfCutlineBlendDist != 0.0) |
2123 | 0 | CPLCreateXMLElementAndValue( |
2124 | 0 | psTree, "CutlineBlendDist", |
2125 | 0 | CPLString().Printf("%.5g", psWO->dfCutlineBlendDist)); |
2126 | |
|
2127 | 0 | return psTree; |
2128 | 0 | } |
2129 | | |
2130 | | /************************************************************************/ |
2131 | | /* GDALDeserializeWarpOptions() */ |
2132 | | /************************************************************************/ |
2133 | | |
2134 | | GDALWarpOptions *CPL_STDCALL GDALDeserializeWarpOptions(CPLXMLNode *psTree) |
2135 | | |
2136 | 0 | { |
2137 | 0 | CPLErrorReset(); |
2138 | | |
2139 | | /* -------------------------------------------------------------------- */ |
2140 | | /* Verify this is the right kind of object. */ |
2141 | | /* -------------------------------------------------------------------- */ |
2142 | 0 | if (psTree == nullptr || psTree->eType != CXT_Element || |
2143 | 0 | !EQUAL(psTree->pszValue, "GDALWarpOptions")) |
2144 | 0 | { |
2145 | 0 | CPLError(CE_Failure, CPLE_AppDefined, |
2146 | 0 | "Wrong node, unable to deserialize GDALWarpOptions."); |
2147 | 0 | return nullptr; |
2148 | 0 | } |
2149 | | |
2150 | | /* -------------------------------------------------------------------- */ |
2151 | | /* Create pre-initialized warp options. */ |
2152 | | /* -------------------------------------------------------------------- */ |
2153 | 0 | GDALWarpOptions *psWO = GDALCreateWarpOptions(); |
2154 | | |
2155 | | /* -------------------------------------------------------------------- */ |
2156 | | /* Warp memory limit. */ |
2157 | | /* -------------------------------------------------------------------- */ |
2158 | 0 | psWO->dfWarpMemoryLimit = |
2159 | 0 | CPLAtof(CPLGetXMLValue(psTree, "WarpMemoryLimit", "0.0")); |
2160 | | |
2161 | | /* -------------------------------------------------------------------- */ |
2162 | | /* resample algorithm */ |
2163 | | /* -------------------------------------------------------------------- */ |
2164 | 0 | const char *pszValue = CPLGetXMLValue(psTree, "ResampleAlg", "Default"); |
2165 | |
|
2166 | 0 | if (EQUAL(pszValue, "NearestNeighbour")) |
2167 | 0 | psWO->eResampleAlg = GRA_NearestNeighbour; |
2168 | 0 | else if (EQUAL(pszValue, "Bilinear")) |
2169 | 0 | psWO->eResampleAlg = GRA_Bilinear; |
2170 | 0 | else if (EQUAL(pszValue, "Cubic")) |
2171 | 0 | psWO->eResampleAlg = GRA_Cubic; |
2172 | 0 | else if (EQUAL(pszValue, "CubicSpline")) |
2173 | 0 | psWO->eResampleAlg = GRA_CubicSpline; |
2174 | 0 | else if (EQUAL(pszValue, "Lanczos")) |
2175 | 0 | psWO->eResampleAlg = GRA_Lanczos; |
2176 | 0 | else if (EQUAL(pszValue, "Average")) |
2177 | 0 | psWO->eResampleAlg = GRA_Average; |
2178 | 0 | else if (EQUAL(pszValue, "RootMeanSquare")) |
2179 | 0 | psWO->eResampleAlg = GRA_RMS; |
2180 | 0 | else if (EQUAL(pszValue, "Mode")) |
2181 | 0 | psWO->eResampleAlg = GRA_Mode; |
2182 | 0 | else if (EQUAL(pszValue, "Maximum")) |
2183 | 0 | psWO->eResampleAlg = GRA_Max; |
2184 | 0 | else if (EQUAL(pszValue, "Minimum")) |
2185 | 0 | psWO->eResampleAlg = GRA_Min; |
2186 | 0 | else if (EQUAL(pszValue, "Median")) |
2187 | 0 | psWO->eResampleAlg = GRA_Med; |
2188 | 0 | else if (EQUAL(pszValue, "Quartile1")) |
2189 | 0 | psWO->eResampleAlg = GRA_Q1; |
2190 | 0 | else if (EQUAL(pszValue, "Quartile3")) |
2191 | 0 | psWO->eResampleAlg = GRA_Q3; |
2192 | 0 | else if (EQUAL(pszValue, "Sum")) |
2193 | 0 | psWO->eResampleAlg = GRA_Sum; |
2194 | 0 | else if (EQUAL(pszValue, "Default")) |
2195 | 0 | /* leave as is */; |
2196 | 0 | else |
2197 | 0 | { |
2198 | 0 | CPLError(CE_Failure, CPLE_AppDefined, |
2199 | 0 | "Unrecognised ResampleAlg value '%s'.", pszValue); |
2200 | 0 | } |
2201 | | |
2202 | | /* -------------------------------------------------------------------- */ |
2203 | | /* Working data type. */ |
2204 | | /* -------------------------------------------------------------------- */ |
2205 | 0 | psWO->eWorkingDataType = GDALGetDataTypeByName( |
2206 | 0 | CPLGetXMLValue(psTree, "WorkingDataType", "Unknown")); |
2207 | | |
2208 | | /* -------------------------------------------------------------------- */ |
2209 | | /* Name/value warp options. */ |
2210 | | /* -------------------------------------------------------------------- */ |
2211 | 0 | for (CPLXMLNode *psItem = psTree->psChild; psItem != nullptr; |
2212 | 0 | psItem = psItem->psNext) |
2213 | 0 | { |
2214 | 0 | if (psItem->eType == CXT_Element && EQUAL(psItem->pszValue, "Option")) |
2215 | 0 | { |
2216 | 0 | const char *pszName = CPLGetXMLValue(psItem, "Name", nullptr); |
2217 | 0 | pszValue = CPLGetXMLValue(psItem, "", nullptr); |
2218 | |
|
2219 | 0 | if (pszName != nullptr && pszValue != nullptr) |
2220 | 0 | { |
2221 | 0 | psWO->papszWarpOptions = |
2222 | 0 | CSLSetNameValue(psWO->papszWarpOptions, pszName, pszValue); |
2223 | 0 | } |
2224 | 0 | } |
2225 | 0 | } |
2226 | | |
2227 | | /* -------------------------------------------------------------------- */ |
2228 | | /* Source Dataset. */ |
2229 | | /* -------------------------------------------------------------------- */ |
2230 | 0 | pszValue = CPLGetXMLValue(psTree, "SourceDataset", nullptr); |
2231 | |
|
2232 | 0 | if (pszValue != nullptr) |
2233 | 0 | { |
2234 | 0 | CPLXMLNode *psGeoLocNode = |
2235 | 0 | CPLSearchXMLNode(psTree, "GeoLocTransformer"); |
2236 | 0 | if (psGeoLocNode) |
2237 | 0 | { |
2238 | 0 | CPLCreateXMLElementAndValue(psGeoLocNode, "SourceDataset", |
2239 | 0 | pszValue); |
2240 | 0 | } |
2241 | |
|
2242 | 0 | CPLConfigOptionSetter oSetter("CPL_ALLOW_VSISTDIN", "NO", true); |
2243 | |
|
2244 | 0 | char **papszOpenOptions = GDALDeserializeOpenOptionsFromXML(psTree); |
2245 | 0 | psWO->hSrcDS = |
2246 | 0 | GDALOpenEx(pszValue, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, |
2247 | 0 | nullptr, papszOpenOptions, nullptr); |
2248 | 0 | CSLDestroy(papszOpenOptions); |
2249 | 0 | } |
2250 | | |
2251 | | /* -------------------------------------------------------------------- */ |
2252 | | /* Destination Dataset. */ |
2253 | | /* -------------------------------------------------------------------- */ |
2254 | 0 | pszValue = CPLGetXMLValue(psTree, "DestinationDataset", nullptr); |
2255 | |
|
2256 | 0 | if (pszValue != nullptr) |
2257 | 0 | { |
2258 | 0 | psWO->hDstDS = GDALOpenShared(pszValue, GA_Update); |
2259 | 0 | } |
2260 | | |
2261 | | /* -------------------------------------------------------------------- */ |
2262 | | /* First, count band mappings so we can establish the bandcount. */ |
2263 | | /* -------------------------------------------------------------------- */ |
2264 | 0 | CPLXMLNode *psBandTree = CPLGetXMLNode(psTree, "BandList"); |
2265 | |
|
2266 | 0 | int nBandCount = 0; |
2267 | 0 | CPLXMLNode *psBand = psBandTree ? psBandTree->psChild : nullptr; |
2268 | 0 | for (; psBand != nullptr; psBand = psBand->psNext) |
2269 | 0 | { |
2270 | 0 | if (psBand->eType != CXT_Element || |
2271 | 0 | !EQUAL(psBand->pszValue, "BandMapping")) |
2272 | 0 | continue; |
2273 | | |
2274 | 0 | nBandCount++; |
2275 | 0 | } |
2276 | |
|
2277 | 0 | GDALWarpInitDefaultBandMapping(psWO, nBandCount); |
2278 | | |
2279 | | /* ==================================================================== */ |
2280 | | /* Now actually process each bandmapping. */ |
2281 | | /* ==================================================================== */ |
2282 | 0 | int iBand = 0; |
2283 | |
|
2284 | 0 | psBand = psBandTree ? psBandTree->psChild : nullptr; |
2285 | |
|
2286 | 0 | for (; psBand != nullptr; psBand = psBand->psNext) |
2287 | 0 | { |
2288 | 0 | if (psBand->eType != CXT_Element || |
2289 | 0 | !EQUAL(psBand->pszValue, "BandMapping")) |
2290 | 0 | continue; |
2291 | | |
2292 | | /* -------------------------------------------------------------------- |
2293 | | */ |
2294 | | /* Source band */ |
2295 | | /* -------------------------------------------------------------------- |
2296 | | */ |
2297 | 0 | pszValue = CPLGetXMLValue(psBand, "src", nullptr); |
2298 | 0 | if (pszValue != nullptr) |
2299 | 0 | psWO->panSrcBands[iBand] = atoi(pszValue); |
2300 | | |
2301 | | /* -------------------------------------------------------------------- |
2302 | | */ |
2303 | | /* Destination band. */ |
2304 | | /* -------------------------------------------------------------------- |
2305 | | */ |
2306 | 0 | pszValue = CPLGetXMLValue(psBand, "dst", nullptr); |
2307 | 0 | if (pszValue != nullptr) |
2308 | 0 | psWO->panDstBands[iBand] = atoi(pszValue); |
2309 | |
|
2310 | 0 | const auto NormalizeValue = [](const char *pszValueIn, |
2311 | 0 | GDALDataType eDataType) -> double |
2312 | 0 | { |
2313 | 0 | if (eDataType == GDT_Float32 && |
2314 | 0 | CPLString().Printf( |
2315 | 0 | "%.16g", -std::numeric_limits<float>::max()) == pszValueIn) |
2316 | 0 | { |
2317 | 0 | return std::numeric_limits<float>::lowest(); |
2318 | 0 | } |
2319 | 0 | else if (eDataType == GDT_Float32 && |
2320 | 0 | CPLString().Printf("%.16g", |
2321 | 0 | std::numeric_limits<float>::max()) == |
2322 | 0 | pszValueIn) |
2323 | 0 | { |
2324 | 0 | return std::numeric_limits<float>::max(); |
2325 | 0 | } |
2326 | 0 | else |
2327 | 0 | { |
2328 | 0 | return CPLAtof(pszValueIn); |
2329 | 0 | } |
2330 | 0 | }; |
2331 | | |
2332 | | /* -------------------------------------------------------------------- |
2333 | | */ |
2334 | | /* Source nodata. */ |
2335 | | /* -------------------------------------------------------------------- |
2336 | | */ |
2337 | 0 | pszValue = CPLGetXMLValue(psBand, "SrcNoDataReal", nullptr); |
2338 | 0 | if (pszValue != nullptr) |
2339 | 0 | { |
2340 | 0 | GDALWarpInitSrcNoDataReal(psWO, -1.1e20); |
2341 | 0 | psWO->padfSrcNoDataReal[iBand] = |
2342 | 0 | NormalizeValue(pszValue, psWO->eWorkingDataType); |
2343 | 0 | } |
2344 | |
|
2345 | 0 | pszValue = CPLGetXMLValue(psBand, "SrcNoDataImag", nullptr); |
2346 | 0 | if (pszValue != nullptr) |
2347 | 0 | { |
2348 | 0 | GDALWarpInitSrcNoDataImag(psWO, 0); |
2349 | 0 | psWO->padfSrcNoDataImag[iBand] = CPLAtof(pszValue); |
2350 | 0 | } |
2351 | | |
2352 | | /* -------------------------------------------------------------------- |
2353 | | */ |
2354 | | /* Destination nodata. */ |
2355 | | /* -------------------------------------------------------------------- |
2356 | | */ |
2357 | 0 | pszValue = CPLGetXMLValue(psBand, "DstNoDataReal", nullptr); |
2358 | 0 | if (pszValue != nullptr) |
2359 | 0 | { |
2360 | 0 | GDALWarpInitDstNoDataReal(psWO, -1.1e20); |
2361 | 0 | psWO->padfDstNoDataReal[iBand] = |
2362 | 0 | NormalizeValue(pszValue, psWO->eWorkingDataType); |
2363 | 0 | } |
2364 | |
|
2365 | 0 | pszValue = CPLGetXMLValue(psBand, "DstNoDataImag", nullptr); |
2366 | 0 | if (pszValue != nullptr) |
2367 | 0 | { |
2368 | 0 | GDALWarpInitDstNoDataImag(psWO, 0); |
2369 | 0 | psWO->padfDstNoDataImag[iBand] = CPLAtof(pszValue); |
2370 | 0 | } |
2371 | |
|
2372 | 0 | iBand++; |
2373 | 0 | } |
2374 | | |
2375 | | /* -------------------------------------------------------------------- */ |
2376 | | /* Alpha bands. */ |
2377 | | /* -------------------------------------------------------------------- */ |
2378 | 0 | psWO->nSrcAlphaBand = atoi(CPLGetXMLValue(psTree, "SrcAlphaBand", "0")); |
2379 | 0 | psWO->nDstAlphaBand = atoi(CPLGetXMLValue(psTree, "DstAlphaBand", "0")); |
2380 | | |
2381 | | /* -------------------------------------------------------------------- */ |
2382 | | /* Cutline. */ |
2383 | | /* -------------------------------------------------------------------- */ |
2384 | 0 | const char *pszWKT = CPLGetXMLValue(psTree, "Cutline", nullptr); |
2385 | 0 | if (pszWKT) |
2386 | 0 | { |
2387 | 0 | char *pszWKTTemp = const_cast<char *>(pszWKT); |
2388 | 0 | OGRGeometryH hCutline = nullptr; |
2389 | 0 | OGR_G_CreateFromWkt(&pszWKTTemp, nullptr, &hCutline); |
2390 | 0 | psWO->hCutline = hCutline; |
2391 | 0 | } |
2392 | |
|
2393 | 0 | psWO->dfCutlineBlendDist = |
2394 | 0 | CPLAtof(CPLGetXMLValue(psTree, "CutlineBlendDist", "0")); |
2395 | | |
2396 | | /* -------------------------------------------------------------------- */ |
2397 | | /* Transformation. */ |
2398 | | /* -------------------------------------------------------------------- */ |
2399 | 0 | CPLXMLNode *psTransformer = CPLGetXMLNode(psTree, "Transformer"); |
2400 | |
|
2401 | 0 | if (psTransformer != nullptr && psTransformer->psChild != nullptr) |
2402 | 0 | { |
2403 | 0 | GDALDeserializeTransformer(psTransformer->psChild, |
2404 | 0 | &(psWO->pfnTransformer), |
2405 | 0 | &(psWO->pTransformerArg)); |
2406 | 0 | } |
2407 | | |
2408 | | /* -------------------------------------------------------------------- */ |
2409 | | /* If any error has occurred, cleanup else return success. */ |
2410 | | /* -------------------------------------------------------------------- */ |
2411 | 0 | if (CPLGetLastErrorType() != CE_None) |
2412 | 0 | { |
2413 | 0 | if (psWO->pTransformerArg) |
2414 | 0 | { |
2415 | 0 | GDALDestroyTransformer(psWO->pTransformerArg); |
2416 | 0 | psWO->pTransformerArg = nullptr; |
2417 | 0 | } |
2418 | 0 | if (psWO->hSrcDS != nullptr) |
2419 | 0 | { |
2420 | 0 | GDALClose(psWO->hSrcDS); |
2421 | 0 | psWO->hSrcDS = nullptr; |
2422 | 0 | } |
2423 | 0 | if (psWO->hDstDS != nullptr) |
2424 | 0 | { |
2425 | 0 | GDALClose(psWO->hDstDS); |
2426 | 0 | psWO->hDstDS = nullptr; |
2427 | 0 | } |
2428 | 0 | GDALDestroyWarpOptions(psWO); |
2429 | 0 | return nullptr; |
2430 | 0 | } |
2431 | | |
2432 | 0 | return psWO; |
2433 | 0 | } |
2434 | | |
2435 | | /************************************************************************/ |
2436 | | /* GDALGetWarpResampleAlg() */ |
2437 | | /************************************************************************/ |
2438 | | |
2439 | | /** Return a GDALResampleAlg from a string */ |
2440 | | bool GDALGetWarpResampleAlg(const char *pszResampling, |
2441 | | GDALResampleAlg &eResampleAlg, bool bThrow) |
2442 | 0 | { |
2443 | 0 | if (STARTS_WITH_CI(pszResampling, "near")) |
2444 | 0 | eResampleAlg = GRA_NearestNeighbour; |
2445 | 0 | else if (EQUAL(pszResampling, "bilinear")) |
2446 | 0 | eResampleAlg = GRA_Bilinear; |
2447 | 0 | else if (EQUAL(pszResampling, "cubic")) |
2448 | 0 | eResampleAlg = GRA_Cubic; |
2449 | 0 | else if (EQUAL(pszResampling, "cubicspline")) |
2450 | 0 | eResampleAlg = GRA_CubicSpline; |
2451 | 0 | else if (EQUAL(pszResampling, "lanczos")) |
2452 | 0 | eResampleAlg = GRA_Lanczos; |
2453 | 0 | else if (EQUAL(pszResampling, "average")) |
2454 | 0 | eResampleAlg = GRA_Average; |
2455 | 0 | else if (EQUAL(pszResampling, "rms")) |
2456 | 0 | eResampleAlg = GRA_RMS; |
2457 | 0 | else if (EQUAL(pszResampling, "mode")) |
2458 | 0 | eResampleAlg = GRA_Mode; |
2459 | 0 | else if (EQUAL(pszResampling, "max")) |
2460 | 0 | eResampleAlg = GRA_Max; |
2461 | 0 | else if (EQUAL(pszResampling, "min")) |
2462 | 0 | eResampleAlg = GRA_Min; |
2463 | 0 | else if (EQUAL(pszResampling, "med")) |
2464 | 0 | eResampleAlg = GRA_Med; |
2465 | 0 | else if (EQUAL(pszResampling, "q1")) |
2466 | 0 | eResampleAlg = GRA_Q1; |
2467 | 0 | else if (EQUAL(pszResampling, "q3")) |
2468 | 0 | eResampleAlg = GRA_Q3; |
2469 | 0 | else if (EQUAL(pszResampling, "sum")) |
2470 | 0 | eResampleAlg = GRA_Sum; |
2471 | 0 | else |
2472 | 0 | { |
2473 | 0 | if (bThrow) |
2474 | 0 | { |
2475 | 0 | throw std::invalid_argument("Unknown resampling method"); |
2476 | 0 | } |
2477 | 0 | else |
2478 | 0 | { |
2479 | 0 | CPLError(CE_Failure, CPLE_IllegalArg, |
2480 | 0 | "Unknown resampling method: %s.", pszResampling); |
2481 | 0 | return false; |
2482 | 0 | } |
2483 | 0 | } |
2484 | 0 | return true; |
2485 | 0 | } |