1 : /******************************************************************************
2 : * $Id: gdalwarpkernel.cpp 25304 2012-12-13 21:03:58Z rouault $
3 : *
4 : * Project: High Performance Image Reprojector
5 : * Purpose: Implementation of the GDALWarpKernel class. Implements the actual
6 : * image warping for a "chunk" of input and output imagery already
7 : * loaded into memory.
8 : * Author: Frank Warmerdam, warmerdam@pobox.com
9 : *
10 : ******************************************************************************
11 : * Copyright (c) 2003, Frank Warmerdam <warmerdam@pobox.com>
12 : *
13 : * Permission is hereby granted, free of charge, to any person obtaining a
14 : * copy of this software and associated documentation files (the "Software"),
15 : * to deal in the Software without restriction, including without limitation
16 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 : * and/or sell copies of the Software, and to permit persons to whom the
18 : * Software is furnished to do so, subject to the following conditions:
19 : *
20 : * The above copyright notice and this permission notice shall be included
21 : * in all copies or substantial portions of the Software.
22 : *
23 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 : * DEALINGS IN THE SOFTWARE.
30 : ****************************************************************************/
31 :
32 : #include "gdalwarper.h"
33 : #include "gdal_alg_priv.h"
34 : #include "cpl_string.h"
35 : #include "gdalwarpkernel_opencl.h"
36 : #include "cpl_atomic_ops.h"
37 : #include "cpl_multiproc.h"
38 :
39 : CPL_CVSID("$Id: gdalwarpkernel.cpp 25304 2012-12-13 21:03:58Z rouault $");
40 :
41 : static const int anGWKFilterRadius[] =
42 : {
43 : 0, // Nearest neighbour
44 : 1, // Bilinear
45 : 2, // Cubic Convolution
46 : 2, // Cubic B-Spline
47 : 3 // Lanczos windowed sinc
48 : };
49 :
50 : /* Used in gdalwarpoperation.cpp */
51 777 : int GWKGetFilterRadius(GDALResampleAlg eResampleAlg)
52 : {
53 777 : return anGWKFilterRadius[eResampleAlg];
54 : }
55 :
56 : #ifdef HAVE_OPENCL
57 : static CPLErr GWKOpenCLCase( GDALWarpKernel * );
58 : #endif
59 :
60 : static CPLErr GWKGeneralCase( GDALWarpKernel * );
61 : static CPLErr GWKNearestNoMasksByte( GDALWarpKernel *poWK );
62 : static CPLErr GWKBilinearNoMasksByte( GDALWarpKernel *poWK );
63 : static CPLErr GWKCubicNoMasksByte( GDALWarpKernel *poWK );
64 : static CPLErr GWKCubicSplineNoMasksByte( GDALWarpKernel *poWK );
65 : static CPLErr GWKNearestByte( GDALWarpKernel *poWK );
66 : static CPLErr GWKNearestNoMasksShort( GDALWarpKernel *poWK );
67 : static CPLErr GWKBilinearNoMasksShort( GDALWarpKernel *poWK );
68 : static CPLErr GWKCubicNoMasksShort( GDALWarpKernel *poWK );
69 : static CPLErr GWKCubicSplineNoMasksShort( GDALWarpKernel *poWK );
70 : static CPLErr GWKNearestShort( GDALWarpKernel *poWK );
71 : static CPLErr GWKNearestNoMasksFloat( GDALWarpKernel *poWK );
72 : static CPLErr GWKNearestFloat( GDALWarpKernel *poWK );
73 :
74 : /************************************************************************/
75 : /* GWKJobStruct */
76 : /************************************************************************/
77 :
78 : typedef struct _GWKJobStruct GWKJobStruct;
79 :
80 : struct _GWKJobStruct
81 : {
82 : void *hThread;
83 : GDALWarpKernel *poWK;
84 : int iYMin;
85 : int iYMax;
86 : volatile int *pnCounter;
87 : volatile int *pbStop;
88 : void *hCond;
89 : void *hCondMutex;
90 : int (*pfnProgress)(GWKJobStruct* psJob);
91 : void *pTransformerArg;
92 : } ;
93 :
94 : /************************************************************************/
95 : /* GWKProgressThread() */
96 : /************************************************************************/
97 :
98 : /* Return TRUE if the computation must be interrupted */
99 261 : static int GWKProgressThread(GWKJobStruct* psJob)
100 : {
101 261 : CPLAcquireMutex(psJob->hCondMutex, 1.0);
102 262 : (*(psJob->pnCounter)) ++;
103 262 : CPLCondSignal(psJob->hCond);
104 262 : int bStop = *(psJob->pbStop);
105 262 : CPLReleaseMutex(psJob->hCondMutex);
106 :
107 262 : return bStop;
108 : }
109 :
110 : /************************************************************************/
111 : /* GWKProgressMonoThread() */
112 : /************************************************************************/
113 :
114 : /* Return TRUE if the computation must be interrupted */
115 44726 : static int GWKProgressMonoThread(GWKJobStruct* psJob)
116 : {
117 44726 : GDALWarpKernel *poWK = psJob->poWK;
118 44726 : int nCounter = ++(*(psJob->pnCounter));
119 44726 : if( !poWK->pfnProgress( poWK->dfProgressBase + poWK->dfProgressScale *
120 : (nCounter / (double) psJob->iYMax),
121 : "", poWK->pProgress ) )
122 : {
123 1 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
124 1 : *(psJob->pbStop) = TRUE;
125 1 : return TRUE;
126 : }
127 44725 : return FALSE;
128 : }
129 :
130 : /************************************************************************/
131 : /* GWKGenericMonoThread() */
132 : /************************************************************************/
133 :
134 580 : static CPLErr GWKGenericMonoThread( GDALWarpKernel *poWK,
135 : void (*pfnFunc) (void *pUserData) )
136 : {
137 580 : volatile int bStop = FALSE;
138 580 : volatile int nCounter = 0;
139 :
140 : GWKJobStruct sThreadJob;
141 580 : sThreadJob.poWK = poWK;
142 580 : sThreadJob.pnCounter = &nCounter;
143 580 : sThreadJob.iYMin = 0;
144 580 : sThreadJob.iYMax = poWK->nDstYSize;
145 580 : sThreadJob.pbStop = &bStop;
146 580 : sThreadJob.hCond = NULL;
147 580 : sThreadJob.hCondMutex = NULL;
148 580 : sThreadJob.hThread = NULL;
149 580 : sThreadJob.pfnProgress = GWKProgressMonoThread;
150 580 : sThreadJob.pTransformerArg = poWK->pTransformerArg;
151 :
152 580 : pfnFunc(&sThreadJob);
153 :
154 580 : return !bStop ? CE_None : CE_Failure;
155 : }
156 :
157 : /************************************************************************/
158 : /* GWKRun() */
159 : /************************************************************************/
160 :
161 583 : static CPLErr GWKRun( GDALWarpKernel *poWK,
162 : const char* pszFuncName,
163 : void (*pfnFunc) (void *pUserData) )
164 :
165 : {
166 583 : int nDstYSize = poWK->nDstYSize;
167 :
168 : CPLDebug( "GDAL", "GDALWarpKernel()::%s()\n"
169 : "Src=%d,%d,%dx%d Dst=%d,%d,%dx%d",
170 : pszFuncName,
171 : poWK->nSrcXOff, poWK->nSrcYOff,
172 : poWK->nSrcXSize, poWK->nSrcYSize,
173 : poWK->nDstXOff, poWK->nDstYOff,
174 583 : poWK->nDstXSize, poWK->nDstYSize );
175 :
176 583 : if( !poWK->pfnProgress( poWK->dfProgressBase, "", poWK->pProgress ) )
177 : {
178 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
179 0 : return CE_Failure;
180 : }
181 :
182 583 : const char* pszWarpThreads = CSLFetchNameValue(poWK->papszWarpOptions, "NUM_THREADS");
183 : int nThreads;
184 583 : if (pszWarpThreads == NULL)
185 583 : pszWarpThreads = CPLGetConfigOption("GDAL_NUM_THREADS", "1");
186 583 : if (EQUAL(pszWarpThreads, "ALL_CPUS"))
187 1 : nThreads = CPLGetNumCPUs();
188 : else
189 582 : nThreads = atoi(pszWarpThreads);
190 583 : if (nThreads > 128)
191 0 : nThreads = 128;
192 583 : if (nThreads >= nDstYSize / 2)
193 165 : nThreads = nDstYSize / 2;
194 :
195 583 : if (nThreads <= 1)
196 : {
197 580 : return GWKGenericMonoThread(poWK, pfnFunc);
198 : }
199 : else
200 : {
201 : GWKJobStruct* pasThreadJob =
202 3 : (GWKJobStruct*)CPLCalloc(sizeof(GWKJobStruct), nThreads);
203 :
204 : /* -------------------------------------------------------------------- */
205 : /* Duplicate pTransformerArg per thread. */
206 : /* -------------------------------------------------------------------- */
207 : int i;
208 3 : int bTransformerCloningSuccess = TRUE;
209 :
210 11 : for(i=0;i<nThreads;i++)
211 : {
212 8 : pasThreadJob[i].pTransformerArg = GDALCloneTransformer(poWK->pTransformerArg);
213 8 : if( pasThreadJob[i].pTransformerArg == NULL )
214 : {
215 0 : CPLDebug("WARP", "Cannot deserialize transformer");
216 0 : bTransformerCloningSuccess = FALSE;
217 0 : break;
218 : }
219 : }
220 :
221 3 : if (!bTransformerCloningSuccess)
222 : {
223 0 : for(i=0;i<nThreads;i++)
224 : {
225 0 : if( pasThreadJob[i].pTransformerArg )
226 0 : GDALDestroyTransformer(pasThreadJob[i].pTransformerArg);
227 : }
228 0 : CPLFree(pasThreadJob);
229 :
230 : CPLDebug("WARP", "Cannot duplicate transformer function. "
231 0 : "Falling back to mono-thread computation");
232 0 : return GWKGenericMonoThread(poWK, pfnFunc);
233 : }
234 :
235 3 : void* hCond = CPLCreateCond();
236 3 : if (hCond == NULL)
237 : {
238 0 : for(i=0;i<nThreads;i++)
239 : {
240 0 : if( pasThreadJob[i].pTransformerArg )
241 0 : GDALDestroyTransformer(pasThreadJob[i].pTransformerArg);
242 : }
243 0 : CPLFree(pasThreadJob);
244 :
245 : CPLDebug("WARP", "Multithreading disabled. "
246 0 : "Falling back to mono-thread computation");
247 0 : return GWKGenericMonoThread(poWK, pfnFunc);
248 : }
249 :
250 3 : CPLDebug("WARP", "Using %d threads", nThreads);
251 :
252 3 : void* hCondMutex = CPLCreateMutex(); /* and take implicitely the mutex */
253 :
254 3 : volatile int bStop = FALSE;
255 3 : volatile int nCounter = 0;
256 :
257 : /* -------------------------------------------------------------------- */
258 : /* Lannch worker threads */
259 : /* -------------------------------------------------------------------- */
260 11 : for(i=0;i<nThreads;i++)
261 : {
262 8 : pasThreadJob[i].poWK = poWK;
263 8 : pasThreadJob[i].pnCounter = &nCounter;
264 8 : pasThreadJob[i].iYMin = (int)(((GIntBig)i) * nDstYSize / nThreads);
265 8 : pasThreadJob[i].iYMax = (int)(((GIntBig)(i + 1)) * nDstYSize / nThreads);
266 8 : pasThreadJob[i].pbStop = &bStop;
267 8 : pasThreadJob[i].hCond = hCond;
268 8 : pasThreadJob[i].hCondMutex = hCondMutex;
269 8 : pasThreadJob[i].pfnProgress = GWKProgressThread;
270 8 : pasThreadJob[i].hThread = CPLCreateJoinableThread( pfnFunc,
271 8 : (void*) &pasThreadJob[i] );
272 : }
273 :
274 : /* -------------------------------------------------------------------- */
275 : /* Report progress. */
276 : /* -------------------------------------------------------------------- */
277 253 : while(nCounter < nDstYSize)
278 : {
279 248 : CPLCondWait(hCond, hCondMutex);
280 :
281 248 : if( !poWK->pfnProgress( poWK->dfProgressBase + poWK->dfProgressScale *
282 : (nCounter / (double) nDstYSize),
283 : "", poWK->pProgress ) )
284 : {
285 1 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
286 1 : bStop = TRUE;
287 1 : break;
288 : }
289 : }
290 :
291 : /* Release mutex before joining threads, otherwise they will dead-lock */
292 : /* forever in GWKProgressThread() */
293 3 : CPLReleaseMutex(hCondMutex);
294 :
295 : /* -------------------------------------------------------------------- */
296 : /* Wait for all threads to complete and finish. */
297 : /* -------------------------------------------------------------------- */
298 11 : for(i=0;i<nThreads;i++)
299 : {
300 8 : CPLJoinThread(pasThreadJob[i].hThread);
301 8 : GDALDestroyTransformer(pasThreadJob[i].pTransformerArg);
302 : }
303 :
304 3 : CPLFree(pasThreadJob);
305 3 : CPLDestroyCond(hCond);
306 3 : CPLDestroyMutex(hCondMutex);
307 :
308 3 : return !bStop ? CE_None : CE_Failure;
309 : }
310 : }
311 :
312 : /************************************************************************/
313 : /* ==================================================================== */
314 : /* GDALWarpKernel */
315 : /* ==================================================================== */
316 : /************************************************************************/
317 :
318 : /**
319 : * \class GDALWarpKernel "gdalwarper.h"
320 : *
321 : * Low level image warping class.
322 : *
323 : * This class is responsible for low level image warping for one
324 : * "chunk" of imagery. The class is essentially a structure with all
325 : * data members public - primarily so that new special-case functions
326 : * can be added without changing the class declaration.
327 : *
328 : * Applications are normally intended to interactive with warping facilities
329 : * through the GDALWarpOperation class, though the GDALWarpKernel can in
330 : * theory be used directly if great care is taken in setting up the
331 : * control data.
332 : *
333 : * <h3>Design Issues</h3>
334 : *
335 : * My intention is that PerformWarp() would analyse the setup in terms
336 : * of the datatype, resampling type, and validity/density mask usage and
337 : * pick one of many specific implementations of the warping algorithm over
338 : * a continuim of optimization vs. generality. At one end there will be a
339 : * reference general purpose implementation of the algorithm that supports
340 : * any data type (working internally in double precision complex), all three
341 : * resampling types, and any or all of the validity/density masks. At the
342 : * other end would be highly optimized algorithms for common cases like
343 : * nearest neighbour resampling on GDT_Byte data with no masks.
344 : *
345 : * The full set of optimized versions have not been decided but we should
346 : * expect to have at least:
347 : * - One for each resampling algorithm for 8bit data with no masks.
348 : * - One for each resampling algorithm for float data with no masks.
349 : * - One for each resampling algorithm for float data with any/all masks
350 : * (essentially the generic case for just float data).
351 : * - One for each resampling algorithm for 8bit data with support for
352 : * input validity masks (per band or per pixel). This handles the common
353 : * case of nodata masking.
354 : * - One for each resampling algorithm for float data with support for
355 : * input validity masks (per band or per pixel). This handles the common
356 : * case of nodata masking.
357 : *
358 : * Some of the specializations would operate on all bands in one pass
359 : * (especially the ones without masking would do this), while others might
360 : * process each band individually to reduce code complexity.
361 : *
362 : * <h3>Masking Semantics</h3>
363 : *
364 : * A detailed explanation of the semantics of the validity and density masks,
365 : * and their effects on resampling kernels is needed here.
366 : */
367 :
368 : /************************************************************************/
369 : /* GDALWarpKernel Data Members */
370 : /************************************************************************/
371 :
372 : /**
373 : * \var GDALResampleAlg GDALWarpKernel::eResample;
374 : *
375 : * Resampling algorithm.
376 : *
377 : * The resampling algorithm to use. One of GRA_NearestNeighbour,
378 : * GRA_Bilinear, or GRA_Cubic.
379 : *
380 : * This field is required. GDT_NearestNeighbour may be used as a default
381 : * value.
382 : */
383 :
384 : /**
385 : * \var GDALDataType GDALWarpKernel::eWorkingDataType;
386 : *
387 : * Working pixel data type.
388 : *
389 : * The datatype of pixels in the source image (papabySrcimage) and
390 : * destination image (papabyDstImage) buffers. Note that operations on
391 : * some data types (such as GDT_Byte) may be much better optimized than other
392 : * less common cases.
393 : *
394 : * This field is required. It may not be GDT_Unknown.
395 : */
396 :
397 : /**
398 : * \var int GDALWarpKernel::nBands;
399 : *
400 : * Number of bands.
401 : *
402 : * The number of bands (layers) of imagery being warped. Determines the
403 : * number of entries in the papabySrcImage, papanBandSrcValid,
404 : * and papabyDstImage arrays.
405 : *
406 : * This field is required.
407 : */
408 :
409 : /**
410 : * \var int GDALWarpKernel::nSrcXSize;
411 : *
412 : * Source image width in pixels.
413 : *
414 : * This field is required.
415 : */
416 :
417 : /**
418 : * \var int GDALWarpKernel::nSrcYSize;
419 : *
420 : * Source image height in pixels.
421 : *
422 : * This field is required.
423 : */
424 :
425 : /**
426 : * \var int GDALWarpKernel::papabySrcImage;
427 : *
428 : * Array of source image band data.
429 : *
430 : * This is an array of pointers (of size GDALWarpKernel::nBands) pointers
431 : * to image data. Each individual band of image data is organized as a single
432 : * block of image data in left to right, then bottom to top order. The actual
433 : * type of the image data is determined by GDALWarpKernel::eWorkingDataType.
434 : *
435 : * To access the the pixel value for the (x=3,y=4) pixel (zero based) of
436 : * the second band with eWorkingDataType set to GDT_Float32 use code like
437 : * this:
438 : *
439 : * \code
440 : * float dfPixelValue;
441 : * int nBand = 1; // band indexes are zero based.
442 : * int nPixel = 3; // zero based
443 : * int nLine = 4; // zero based
444 : *
445 : * assert( nPixel >= 0 && nPixel < poKern->nSrcXSize );
446 : * assert( nLine >= 0 && nLine < poKern->nSrcYSize );
447 : * assert( nBand >= 0 && nBand < poKern->nBands );
448 : * dfPixelValue = ((float *) poKern->papabySrcImage[nBand-1])
449 : * [nPixel + nLine * poKern->nSrcXSize];
450 : * \endcode
451 : *
452 : * This field is required.
453 : */
454 :
455 : /**
456 : * \var GUInt32 **GDALWarpKernel::papanBandSrcValid;
457 : *
458 : * Per band validity mask for source pixels.
459 : *
460 : * Array of pixel validity mask layers for each source band. Each of
461 : * the mask layers is the same size (in pixels) as the source image with
462 : * one bit per pixel. Note that it is legal (and common) for this to be
463 : * NULL indicating that none of the pixels are invalidated, or for some
464 : * band validity masks to be NULL in which case all pixels of the band are
465 : * valid. The following code can be used to test the validity of a particular
466 : * pixel.
467 : *
468 : * \code
469 : * int bIsValid = TRUE;
470 : * int nBand = 1; // band indexes are zero based.
471 : * int nPixel = 3; // zero based
472 : * int nLine = 4; // zero based
473 : *
474 : * assert( nPixel >= 0 && nPixel < poKern->nSrcXSize );
475 : * assert( nLine >= 0 && nLine < poKern->nSrcYSize );
476 : * assert( nBand >= 0 && nBand < poKern->nBands );
477 : *
478 : * if( poKern->papanBandSrcValid != NULL
479 : * && poKern->papanBandSrcValid[nBand] != NULL )
480 : * {
481 : * GUInt32 *panBandMask = poKern->papanBandSrcValid[nBand];
482 : * int iPixelOffset = nPixel + nLine * poKern->nSrcXSize;
483 : *
484 : * bIsValid = panBandMask[iPixelOffset>>5]
485 : * & (0x01 << (iPixelOffset & 0x1f));
486 : * }
487 : * \endcode
488 : */
489 :
490 : /**
491 : * \var GUInt32 *GDALWarpKernel::panUnifiedSrcValid;
492 : *
493 : * Per pixel validity mask for source pixels.
494 : *
495 : * A single validity mask layer that applies to the pixels of all source
496 : * bands. It is accessed similarly to papanBandSrcValid, but without the
497 : * extra level of band indirection.
498 : *
499 : * This pointer may be NULL indicating that all pixels are valid.
500 : *
501 : * Note that if both panUnifiedSrcValid, and papanBandSrcValid are available,
502 : * the pixel isn't considered to be valid unless both arrays indicate it is
503 : * valid.
504 : */
505 :
506 : /**
507 : * \var float *GDALWarpKernel::pafUnifiedSrcDensity;
508 : *
509 : * Per pixel density mask for source pixels.
510 : *
511 : * A single density mask layer that applies to the pixels of all source
512 : * bands. It contains values between 0.0 and 1.0 indicating the degree to
513 : * which this pixel should be allowed to contribute to the output result.
514 : *
515 : * This pointer may be NULL indicating that all pixels have a density of 1.0.
516 : *
517 : * The density for a pixel may be accessed like this:
518 : *
519 : * \code
520 : * float fDensity = 1.0;
521 : * int nPixel = 3; // zero based
522 : * int nLine = 4; // zero based
523 : *
524 : * assert( nPixel >= 0 && nPixel < poKern->nSrcXSize );
525 : * assert( nLine >= 0 && nLine < poKern->nSrcYSize );
526 : * if( poKern->pafUnifiedSrcDensity != NULL )
527 : * fDensity = poKern->pafUnifiedSrcDensity
528 : * [nPixel + nLine * poKern->nSrcXSize];
529 : * \endcode
530 : */
531 :
532 : /**
533 : * \var int GDALWarpKernel::nDstXSize;
534 : *
535 : * Width of destination image in pixels.
536 : *
537 : * This field is required.
538 : */
539 :
540 : /**
541 : * \var int GDALWarpKernel::nDstYSize;
542 : *
543 : * Height of destination image in pixels.
544 : *
545 : * This field is required.
546 : */
547 :
548 : /**
549 : * \var GByte **GDALWarpKernel::papabyDstImage;
550 : *
551 : * Array of destination image band data.
552 : *
553 : * This is an array of pointers (of size GDALWarpKernel::nBands) pointers
554 : * to image data. Each individual band of image data is organized as a single
555 : * block of image data in left to right, then bottom to top order. The actual
556 : * type of the image data is determined by GDALWarpKernel::eWorkingDataType.
557 : *
558 : * To access the the pixel value for the (x=3,y=4) pixel (zero based) of
559 : * the second band with eWorkingDataType set to GDT_Float32 use code like
560 : * this:
561 : *
562 : * \code
563 : * float dfPixelValue;
564 : * int nBand = 1; // band indexes are zero based.
565 : * int nPixel = 3; // zero based
566 : * int nLine = 4; // zero based
567 : *
568 : * assert( nPixel >= 0 && nPixel < poKern->nDstXSize );
569 : * assert( nLine >= 0 && nLine < poKern->nDstYSize );
570 : * assert( nBand >= 0 && nBand < poKern->nBands );
571 : * dfPixelValue = ((float *) poKern->papabyDstImage[nBand-1])
572 : * [nPixel + nLine * poKern->nSrcYSize];
573 : * \endcode
574 : *
575 : * This field is required.
576 : */
577 :
578 : /**
579 : * \var GUInt32 *GDALWarpKernel::panDstValid;
580 : *
581 : * Per pixel validity mask for destination pixels.
582 : *
583 : * A single validity mask layer that applies to the pixels of all destination
584 : * bands. It is accessed similarly to papanUnitifiedSrcValid, but based
585 : * on the size of the destination image.
586 : *
587 : * This pointer may be NULL indicating that all pixels are valid.
588 : */
589 :
590 : /**
591 : * \var float *GDALWarpKernel::pafDstDensity;
592 : *
593 : * Per pixel density mask for destination pixels.
594 : *
595 : * A single density mask layer that applies to the pixels of all destination
596 : * bands. It contains values between 0.0 and 1.0.
597 : *
598 : * This pointer may be NULL indicating that all pixels have a density of 1.0.
599 : *
600 : * The density for a pixel may be accessed like this:
601 : *
602 : * \code
603 : * float fDensity = 1.0;
604 : * int nPixel = 3; // zero based
605 : * int nLine = 4; // zero based
606 : *
607 : * assert( nPixel >= 0 && nPixel < poKern->nDstXSize );
608 : * assert( nLine >= 0 && nLine < poKern->nDstYSize );
609 : * if( poKern->pafDstDensity != NULL )
610 : * fDensity = poKern->pafDstDensity[nPixel + nLine * poKern->nDstXSize];
611 : * \endcode
612 : */
613 :
614 : /**
615 : * \var int GDALWarpKernel::nSrcXOff;
616 : *
617 : * X offset to source pixel coordinates for transformation.
618 : *
619 : * See pfnTransformer.
620 : *
621 : * This field is required.
622 : */
623 :
624 : /**
625 : * \var int GDALWarpKernel::nSrcYOff;
626 : *
627 : * Y offset to source pixel coordinates for transformation.
628 : *
629 : * See pfnTransformer.
630 : *
631 : * This field is required.
632 : */
633 :
634 : /**
635 : * \var int GDALWarpKernel::nDstXOff;
636 : *
637 : * X offset to destination pixel coordinates for transformation.
638 : *
639 : * See pfnTransformer.
640 : *
641 : * This field is required.
642 : */
643 :
644 : /**
645 : * \var int GDALWarpKernel::nDstYOff;
646 : *
647 : * Y offset to destination pixel coordinates for transformation.
648 : *
649 : * See pfnTransformer.
650 : *
651 : * This field is required.
652 : */
653 :
654 : /**
655 : * \var GDALTransformerFunc GDALWarpKernel::pfnTransformer;
656 : *
657 : * Source/destination location transformer.
658 : *
659 : * The function to call to transform coordinates between source image
660 : * pixel/line coordinates and destination image pixel/line coordinates.
661 : * See GDALTransformerFunc() for details of the semantics of this function.
662 : *
663 : * The GDALWarpKern algorithm will only ever use this transformer in
664 : * "destination to source" mode (bDstToSrc=TRUE), and will always pass
665 : * partial or complete scanlines of points in the destination image as
666 : * input. This means, amoung other things, that it is safe to the the
667 : * approximating transform GDALApproxTransform() as the transformation
668 : * function.
669 : *
670 : * Source and destination images may be subsets of a larger overall image.
671 : * The transformation algorithms will expect and return pixel/line coordinates
672 : * in terms of this larger image, so coordinates need to be offset by
673 : * the offsets specified in nSrcXOff, nSrcYOff, nDstXOff, and nDstYOff before
674 : * passing to pfnTransformer, and after return from it.
675 : *
676 : * The GDALWarpKernel::pfnTransformerArg value will be passed as the callback
677 : * data to this function when it is called.
678 : *
679 : * This field is required.
680 : */
681 :
682 : /**
683 : * \var void *GDALWarpKernel::pTransformerArg;
684 : *
685 : * Callback data for pfnTransformer.
686 : *
687 : * This field may be NULL if not required for the pfnTransformer being used.
688 : */
689 :
690 : /**
691 : * \var GDALProgressFunc GDALWarpKernel::pfnProgress;
692 : *
693 : * The function to call to report progress of the algorithm, and to check
694 : * for a requested termination of the operation. It operates according to
695 : * GDALProgressFunc() semantics.
696 : *
697 : * Generally speaking the progress function will be invoked for each
698 : * scanline of the destination buffer that has been processed.
699 : *
700 : * This field may be NULL (internally set to GDALDummyProgress()).
701 : */
702 :
703 : /**
704 : * \var void *GDALWarpKernel::pProgress;
705 : *
706 : * Callback data for pfnProgress.
707 : *
708 : * This field may be NULL if not required for the pfnProgress being used.
709 : */
710 :
711 :
712 : /************************************************************************/
713 : /* GDALWarpKernel() */
714 : /************************************************************************/
715 :
716 585 : GDALWarpKernel::GDALWarpKernel()
717 :
718 : {
719 585 : eResample = GRA_NearestNeighbour;
720 585 : eWorkingDataType = GDT_Unknown;
721 585 : nBands = 0;
722 585 : nDstXOff = 0;
723 585 : nDstYOff = 0;
724 585 : nDstXSize = 0;
725 585 : nDstYSize = 0;
726 585 : nSrcXOff = 0;
727 585 : nSrcYOff = 0;
728 585 : nSrcXSize = 0;
729 585 : nSrcYSize = 0;
730 585 : dfXScale = 1.0;
731 585 : dfYScale = 1.0;
732 585 : dfXFilter = 0.0;
733 585 : dfYFilter = 0.0;
734 585 : nXRadius = 0;
735 585 : nYRadius = 0;
736 585 : nFiltInitX = 0;
737 585 : nFiltInitY = 0;
738 585 : pafDstDensity = NULL;
739 585 : pafUnifiedSrcDensity = NULL;
740 585 : panDstValid = NULL;
741 585 : panUnifiedSrcValid = NULL;
742 585 : papabyDstImage = NULL;
743 585 : papabySrcImage = NULL;
744 585 : papanBandSrcValid = NULL;
745 585 : pfnProgress = GDALDummyProgress;
746 585 : pProgress = NULL;
747 585 : dfProgressBase = 0.0;
748 585 : dfProgressScale = 1.0;
749 585 : pfnTransformer = NULL;
750 585 : pTransformerArg = NULL;
751 585 : papszWarpOptions = NULL;
752 585 : }
753 :
754 : /************************************************************************/
755 : /* ~GDALWarpKernel() */
756 : /************************************************************************/
757 :
758 585 : GDALWarpKernel::~GDALWarpKernel()
759 :
760 : {
761 585 : }
762 :
763 : /************************************************************************/
764 : /* PerformWarp() */
765 : /************************************************************************/
766 :
767 : /**
768 : * \fn CPLErr GDALWarpKernel::PerformWarp();
769 : *
770 : * This method performs the warp described in the GDALWarpKernel.
771 : *
772 : * @return CE_None on success or CE_Failure if an error occurs.
773 : */
774 :
775 584 : CPLErr GDALWarpKernel::PerformWarp()
776 :
777 : {
778 : CPLErr eErr;
779 :
780 584 : if( (eErr = Validate()) != CE_None )
781 0 : return eErr;
782 :
783 : // See #2445 and #3079
784 584 : if (nSrcXSize <= 0 || nSrcYSize <= 0)
785 : {
786 1 : if ( !pfnProgress( dfProgressBase + dfProgressScale,
787 : "", pProgress ) )
788 : {
789 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
790 0 : return CE_Failure;
791 : }
792 1 : return CE_None;
793 : }
794 :
795 : /* -------------------------------------------------------------------- */
796 : /* Pre-calculate resampling scales and window sizes for filtering. */
797 : /* -------------------------------------------------------------------- */
798 :
799 583 : dfXScale = (double)nDstXSize / nSrcXSize;
800 583 : dfYScale = (double)nDstYSize / nSrcYSize;
801 583 : if( nSrcXSize >= nDstXSize && nSrcXSize <= nDstXSize + 1 + 2 * anGWKFilterRadius[eResample] )
802 456 : dfXScale = 1;
803 583 : if( nSrcYSize >= nDstYSize && nSrcYSize <= nDstYSize + 1 + 2 * anGWKFilterRadius[eResample] )
804 462 : dfYScale = 1;
805 :
806 583 : dfXFilter = anGWKFilterRadius[eResample];
807 583 : dfYFilter = anGWKFilterRadius[eResample];
808 :
809 : nXRadius = ( dfXScale < 1.0 ) ?
810 583 : (int)ceil( dfXFilter / dfXScale ) :(int)dfXFilter;
811 : nYRadius = ( dfYScale < 1.0 ) ?
812 583 : (int)ceil( dfYFilter / dfYScale ) : (int)dfYFilter;
813 :
814 : // Filter window offset depends on the parity of the kernel radius
815 583 : nFiltInitX = ((anGWKFilterRadius[eResample] + 1) % 2) - nXRadius;
816 583 : nFiltInitY = ((anGWKFilterRadius[eResample] + 1) % 2) - nYRadius;
817 :
818 : /* -------------------------------------------------------------------- */
819 : /* Set up resampling functions. */
820 : /* -------------------------------------------------------------------- */
821 583 : if( CSLFetchBoolean( papszWarpOptions, "USE_GENERAL_CASE", FALSE ) )
822 0 : return GWKGeneralCase( this );
823 :
824 : #if defined(HAVE_OPENCL)
825 : if((eWorkingDataType == GDT_Byte
826 : || eWorkingDataType == GDT_CInt16
827 : || eWorkingDataType == GDT_UInt16
828 : || eWorkingDataType == GDT_Int16
829 : || eWorkingDataType == GDT_CFloat32
830 : || eWorkingDataType == GDT_Float32) &&
831 : (eResample == GRA_Bilinear
832 : || eResample == GRA_Cubic
833 : || eResample == GRA_CubicSpline
834 : || eResample == GRA_Lanczos) &&
835 : CSLFetchBoolean( papszWarpOptions, "USE_OPENCL", TRUE ))
836 : {
837 : CPLErr eResult = GWKOpenCLCase( this );
838 :
839 : // CE_Warning tells us a suitable OpenCL environment was not available
840 : // so we fall through to other CPU based methods.
841 : if( eResult != CE_Warning )
842 : return eResult;
843 : }
844 : #endif /* defined HAVE_OPENCL */
845 :
846 583 : if( eWorkingDataType == GDT_Byte
847 : && eResample == GRA_NearestNeighbour
848 : && papanBandSrcValid == NULL
849 : && panUnifiedSrcValid == NULL
850 : && pafUnifiedSrcDensity == NULL
851 : && panDstValid == NULL
852 : && pafDstDensity == NULL )
853 255 : return GWKNearestNoMasksByte( this );
854 :
855 328 : if( eWorkingDataType == GDT_Byte
856 : && eResample == GRA_Bilinear
857 : && papanBandSrcValid == NULL
858 : && panUnifiedSrcValid == NULL
859 : && pafUnifiedSrcDensity == NULL
860 : && panDstValid == NULL
861 : && pafDstDensity == NULL )
862 16 : return GWKBilinearNoMasksByte( this );
863 :
864 312 : if( eWorkingDataType == GDT_Byte
865 : && eResample == GRA_Cubic
866 : && papanBandSrcValid == NULL
867 : && panUnifiedSrcValid == NULL
868 : && pafUnifiedSrcDensity == NULL
869 : && panDstValid == NULL
870 : && pafDstDensity == NULL )
871 10 : return GWKCubicNoMasksByte( this );
872 :
873 302 : if( eWorkingDataType == GDT_Byte
874 : && eResample == GRA_CubicSpline
875 : && papanBandSrcValid == NULL
876 : && panUnifiedSrcValid == NULL
877 : && pafUnifiedSrcDensity == NULL
878 : && panDstValid == NULL
879 : && pafDstDensity == NULL )
880 10 : return GWKCubicSplineNoMasksByte( this );
881 :
882 292 : if( eWorkingDataType == GDT_Byte
883 : && eResample == GRA_NearestNeighbour )
884 28 : return GWKNearestByte( this );
885 :
886 264 : if( (eWorkingDataType == GDT_Int16 || eWorkingDataType == GDT_UInt16)
887 : && eResample == GRA_NearestNeighbour
888 : && papanBandSrcValid == NULL
889 : && panUnifiedSrcValid == NULL
890 : && pafUnifiedSrcDensity == NULL
891 : && panDstValid == NULL
892 : && pafDstDensity == NULL )
893 14 : return GWKNearestNoMasksShort( this );
894 :
895 250 : if( (eWorkingDataType == GDT_Int16 )
896 : && eResample == GRA_Cubic
897 : && papanBandSrcValid == NULL
898 : && panUnifiedSrcValid == NULL
899 : && pafUnifiedSrcDensity == NULL
900 : && panDstValid == NULL
901 : && pafDstDensity == NULL )
902 8 : return GWKCubicNoMasksShort( this );
903 :
904 242 : if( (eWorkingDataType == GDT_Int16 )
905 : && eResample == GRA_CubicSpline
906 : && papanBandSrcValid == NULL
907 : && panUnifiedSrcValid == NULL
908 : && pafUnifiedSrcDensity == NULL
909 : && panDstValid == NULL
910 : && pafDstDensity == NULL )
911 8 : return GWKCubicSplineNoMasksShort( this );
912 :
913 234 : if( (eWorkingDataType == GDT_Int16 )
914 : && eResample == GRA_Bilinear
915 : && papanBandSrcValid == NULL
916 : && panUnifiedSrcValid == NULL
917 : && pafUnifiedSrcDensity == NULL
918 : && panDstValid == NULL
919 : && pafDstDensity == NULL )
920 20 : return GWKBilinearNoMasksShort( this );
921 :
922 214 : if( (eWorkingDataType == GDT_Int16 || eWorkingDataType == GDT_UInt16)
923 : && eResample == GRA_NearestNeighbour )
924 7 : return GWKNearestShort( this );
925 :
926 207 : if( eWorkingDataType == GDT_Float32
927 : && eResample == GRA_NearestNeighbour
928 : && papanBandSrcValid == NULL
929 : && panUnifiedSrcValid == NULL
930 : && pafUnifiedSrcDensity == NULL
931 : && panDstValid == NULL
932 : && pafDstDensity == NULL )
933 8 : return GWKNearestNoMasksFloat( this );
934 :
935 199 : if( eWorkingDataType == GDT_Float32
936 : && eResample == GRA_NearestNeighbour )
937 6 : return GWKNearestFloat( this );
938 :
939 193 : return GWKGeneralCase( this );
940 : }
941 :
942 : /************************************************************************/
943 : /* Validate() */
944 : /************************************************************************/
945 :
946 : /**
947 : * \fn CPLErr GDALWarpKernel::Validate()
948 : *
949 : * Check the settings in the GDALWarpKernel, and issue a CPLError()
950 : * (and return CE_Failure) if the configuration is considered to be
951 : * invalid for some reason.
952 : *
953 : * This method will also do some standard defaulting such as setting
954 : * pfnProgress to GDALDummyProgress() if it is NULL.
955 : *
956 : * @return CE_None on success or CE_Failure if an error is detected.
957 : */
958 :
959 584 : CPLErr GDALWarpKernel::Validate()
960 :
961 : {
962 584 : if ( (size_t)eResample >=
963 : (sizeof(anGWKFilterRadius) / sizeof(anGWKFilterRadius[0])) )
964 : {
965 : CPLError( CE_Failure, CPLE_AppDefined,
966 0 : "Unsupported resampling method %d.", (int) eResample );
967 0 : return CE_Failure;
968 : }
969 :
970 : // Safety check for callers that would use GDALWarpKernel without using
971 : // GDALWarpOperation.
972 584 : if( (eResample == GRA_CubicSpline || eResample == GRA_Lanczos) &&
973 : atoi(CSLFetchNameValueDef(papszWarpOptions, "EXTRA_ELTS", "0") ) != WARP_EXTRA_ELTS )
974 : {
975 : CPLError( CE_Failure, CPLE_AppDefined,
976 : "Source arrays must have WARP_EXTRA_ELTS extra elements at their end. "
977 : "See GDALWarpKernel class definition. If this condition is fulfilled, "
978 0 : "define a EXTRA_ELTS=%d warp options", WARP_EXTRA_ELTS);
979 0 : return CE_Failure;
980 : }
981 :
982 584 : return CE_None;
983 : }
984 :
985 : /************************************************************************/
986 : /* GWKOverlayDensity() */
987 : /* */
988 : /* Compute the final density for the destination pixel. This */
989 : /* is a function of the overlay density (passed in) and the */
990 : /* original density. */
991 : /************************************************************************/
992 :
993 2238219 : static void GWKOverlayDensity( GDALWarpKernel *poWK, int iDstOffset,
994 : double dfDensity )
995 : {
996 2238219 : if( dfDensity < 0.0001 || poWK->pafDstDensity == NULL )
997 1789911 : return;
998 :
999 448308 : poWK->pafDstDensity[iDstOffset] = (float)
1000 448308 : ( 1.0 - (1.0 - dfDensity) * (1.0 - poWK->pafDstDensity[iDstOffset]) );
1001 : }
1002 :
1003 : /************************************************************************/
1004 : /* GWKSetPixelValue() */
1005 : /************************************************************************/
1006 :
1007 290867 : static int GWKSetPixelValue( GDALWarpKernel *poWK, int iBand,
1008 : int iDstOffset, double dfDensity,
1009 : double dfReal, double dfImag )
1010 :
1011 : {
1012 290867 : GByte *pabyDst = poWK->papabyDstImage[iBand];
1013 :
1014 : /* -------------------------------------------------------------------- */
1015 : /* If the source density is less than 100% we need to fetch the */
1016 : /* existing destination value, and mix it with the source to */
1017 : /* get the new "to apply" value. Also compute composite */
1018 : /* density. */
1019 : /* */
1020 : /* We avoid mixing if density is very near one or risk mixing */
1021 : /* in very extreme nodata values and causing odd results (#1610) */
1022 : /* -------------------------------------------------------------------- */
1023 290867 : if( dfDensity < 0.9999 )
1024 : {
1025 1472 : double dfDstReal, dfDstImag, dfDstDensity = 1.0;
1026 :
1027 1472 : if( dfDensity < 0.0001 )
1028 0 : return TRUE;
1029 :
1030 1472 : if( poWK->pafDstDensity != NULL )
1031 0 : dfDstDensity = poWK->pafDstDensity[iDstOffset];
1032 1472 : else if( poWK->panDstValid != NULL
1033 0 : && !((poWK->panDstValid[iDstOffset>>5]
1034 : & (0x01 << (iDstOffset & 0x1f))) ) )
1035 0 : dfDstDensity = 0.0;
1036 :
1037 : // It seems like we also ought to be testing panDstValid[] here!
1038 :
1039 1472 : switch( poWK->eWorkingDataType )
1040 : {
1041 : case GDT_Byte:
1042 1472 : dfDstReal = pabyDst[iDstOffset];
1043 1472 : dfDstImag = 0.0;
1044 1472 : break;
1045 :
1046 : case GDT_Int16:
1047 0 : dfDstReal = ((GInt16 *) pabyDst)[iDstOffset];
1048 0 : dfDstImag = 0.0;
1049 0 : break;
1050 :
1051 : case GDT_UInt16:
1052 0 : dfDstReal = ((GUInt16 *) pabyDst)[iDstOffset];
1053 0 : dfDstImag = 0.0;
1054 0 : break;
1055 :
1056 : case GDT_Int32:
1057 0 : dfDstReal = ((GInt32 *) pabyDst)[iDstOffset];
1058 0 : dfDstImag = 0.0;
1059 0 : break;
1060 :
1061 : case GDT_UInt32:
1062 0 : dfDstReal = ((GUInt32 *) pabyDst)[iDstOffset];
1063 0 : dfDstImag = 0.0;
1064 0 : break;
1065 :
1066 : case GDT_Float32:
1067 0 : dfDstReal = ((float *) pabyDst)[iDstOffset];
1068 0 : dfDstImag = 0.0;
1069 0 : break;
1070 :
1071 : case GDT_Float64:
1072 0 : dfDstReal = ((double *) pabyDst)[iDstOffset];
1073 0 : dfDstImag = 0.0;
1074 0 : break;
1075 :
1076 : case GDT_CInt16:
1077 0 : dfDstReal = ((GInt16 *) pabyDst)[iDstOffset*2];
1078 0 : dfDstImag = ((GInt16 *) pabyDst)[iDstOffset*2+1];
1079 0 : break;
1080 :
1081 : case GDT_CInt32:
1082 0 : dfDstReal = ((GInt32 *) pabyDst)[iDstOffset*2];
1083 0 : dfDstImag = ((GInt32 *) pabyDst)[iDstOffset*2+1];
1084 0 : break;
1085 :
1086 : case GDT_CFloat32:
1087 0 : dfDstReal = ((float *) pabyDst)[iDstOffset*2];
1088 0 : dfDstImag = ((float *) pabyDst)[iDstOffset*2+1];
1089 0 : break;
1090 :
1091 : case GDT_CFloat64:
1092 0 : dfDstReal = ((double *) pabyDst)[iDstOffset*2];
1093 0 : dfDstImag = ((double *) pabyDst)[iDstOffset*2+1];
1094 0 : break;
1095 :
1096 : default:
1097 0 : CPLAssert( FALSE );
1098 0 : dfDstDensity = 0.0;
1099 0 : return FALSE;
1100 : }
1101 :
1102 : // the destination density is really only relative to the portion
1103 : // not occluded by the overlay.
1104 1472 : double dfDstInfluence = (1.0 - dfDensity) * dfDstDensity;
1105 :
1106 : dfReal = (dfReal * dfDensity + dfDstReal * dfDstInfluence)
1107 1472 : / (dfDensity + dfDstInfluence);
1108 :
1109 : dfImag = (dfImag * dfDensity + dfDstImag * dfDstInfluence)
1110 1472 : / (dfDensity + dfDstInfluence);
1111 : }
1112 :
1113 : /* -------------------------------------------------------------------- */
1114 : /* Actually apply the destination value. */
1115 : /* */
1116 : /* Avoid using the destination nodata value for integer datatypes */
1117 : /* if by chance it is equal to the computed pixel value. */
1118 : /* -------------------------------------------------------------------- */
1119 :
1120 : #define CLAMP(type,minval,maxval) \
1121 : do { \
1122 : if (dfReal < minval) ((type *) pabyDst)[iDstOffset] = (type)minval; \
1123 : else if (dfReal > maxval) ((type *) pabyDst)[iDstOffset] = (type)maxval; \
1124 : else ((type *) pabyDst)[iDstOffset] = (minval < 0) ? (type)floor(dfReal + 0.5) : (type)(dfReal + 0.5); \
1125 : if (poWK->padfDstNoDataReal != NULL && \
1126 : poWK->padfDstNoDataReal[iBand] == (double)((type *) pabyDst)[iDstOffset]) \
1127 : { \
1128 : if (((type *) pabyDst)[iDstOffset] == minval) \
1129 : ((type *) pabyDst)[iDstOffset] = (type)(minval + 1); \
1130 : else \
1131 : ((type *) pabyDst)[iDstOffset] --; \
1132 : } } while(0)
1133 :
1134 290867 : switch( poWK->eWorkingDataType )
1135 : {
1136 : case GDT_Byte:
1137 287295 : CLAMP(GByte, 0.0, 255.0);
1138 287295 : break;
1139 :
1140 : case GDT_Int16:
1141 63 : CLAMP(GInt16, -32768.0, 32767.0);
1142 63 : break;
1143 :
1144 : case GDT_UInt16:
1145 252 : CLAMP(GUInt16, 0.0, 65535.0);
1146 252 : break;
1147 :
1148 : case GDT_UInt32:
1149 315 : CLAMP(GUInt32, 0.0, 4294967295.0);
1150 315 : break;
1151 :
1152 : case GDT_Int32:
1153 1115 : CLAMP(GInt32, -2147483648.0, 2147483647.0);
1154 1115 : break;
1155 :
1156 : case GDT_Float32:
1157 252 : ((float *) pabyDst)[iDstOffset] = (float) dfReal;
1158 252 : break;
1159 :
1160 : case GDT_Float64:
1161 315 : ((double *) pabyDst)[iDstOffset] = dfReal;
1162 315 : break;
1163 :
1164 : case GDT_CInt16:
1165 315 : if( dfReal < -32768 )
1166 0 : ((GInt16 *) pabyDst)[iDstOffset*2] = -32768;
1167 315 : else if( dfReal > 32767 )
1168 0 : ((GInt16 *) pabyDst)[iDstOffset*2] = 32767;
1169 : else
1170 315 : ((GInt16 *) pabyDst)[iDstOffset*2] = (GInt16) floor(dfReal+0.5);
1171 315 : if( dfImag < -32768 )
1172 0 : ((GInt16 *) pabyDst)[iDstOffset*2+1] = -32768;
1173 315 : else if( dfImag > 32767 )
1174 0 : ((GInt16 *) pabyDst)[iDstOffset*2+1] = 32767;
1175 : else
1176 315 : ((GInt16 *) pabyDst)[iDstOffset*2+1] = (GInt16) floor(dfImag+0.5);
1177 315 : break;
1178 :
1179 : case GDT_CInt32:
1180 315 : if( dfReal < -2147483648.0 )
1181 0 : ((GInt32 *) pabyDst)[iDstOffset*2] = (GInt32) -2147483648.0;
1182 315 : else if( dfReal > 2147483647.0 )
1183 0 : ((GInt32 *) pabyDst)[iDstOffset*2] = (GInt32) 2147483647.0;
1184 : else
1185 315 : ((GInt32 *) pabyDst)[iDstOffset*2] = (GInt32) floor(dfReal+0.5);
1186 315 : if( dfImag < -2147483648.0 )
1187 0 : ((GInt32 *) pabyDst)[iDstOffset*2+1] = (GInt32) -2147483648.0;
1188 315 : else if( dfImag > 2147483647.0 )
1189 0 : ((GInt32 *) pabyDst)[iDstOffset*2+1] = (GInt32) 2147483647.0;
1190 : else
1191 315 : ((GInt32 *) pabyDst)[iDstOffset*2+1] = (GInt32) floor(dfImag+0.5);
1192 315 : break;
1193 :
1194 : case GDT_CFloat32:
1195 315 : ((float *) pabyDst)[iDstOffset*2] = (float) dfReal;
1196 315 : ((float *) pabyDst)[iDstOffset*2+1] = (float) dfImag;
1197 315 : break;
1198 :
1199 : case GDT_CFloat64:
1200 315 : ((double *) pabyDst)[iDstOffset*2] = (double) dfReal;
1201 315 : ((double *) pabyDst)[iDstOffset*2+1] = (double) dfImag;
1202 315 : break;
1203 :
1204 : default:
1205 0 : return FALSE;
1206 : }
1207 :
1208 290867 : return TRUE;
1209 : }
1210 :
1211 : /************************************************************************/
1212 : /* GWKGetPixelValue() */
1213 : /************************************************************************/
1214 :
1215 1279 : static int GWKGetPixelValue( GDALWarpKernel *poWK, int iBand,
1216 : int iSrcOffset, double *pdfDensity,
1217 : double *pdfReal, double *pdfImag )
1218 :
1219 : {
1220 1279 : GByte *pabySrc = poWK->papabySrcImage[iBand];
1221 :
1222 1279 : if( poWK->panUnifiedSrcValid != NULL
1223 0 : && !((poWK->panUnifiedSrcValid[iSrcOffset>>5]
1224 : & (0x01 << (iSrcOffset & 0x1f))) ) )
1225 : {
1226 0 : *pdfDensity = 0.0;
1227 0 : return FALSE;
1228 : }
1229 :
1230 2879 : if( poWK->papanBandSrcValid != NULL
1231 800 : && poWK->papanBandSrcValid[iBand] != NULL
1232 800 : && !((poWK->papanBandSrcValid[iBand][iSrcOffset>>5]
1233 : & (0x01 << (iSrcOffset & 0x1f)))) )
1234 : {
1235 0 : *pdfDensity = 0.0;
1236 0 : return FALSE;
1237 : }
1238 :
1239 1279 : switch( poWK->eWorkingDataType )
1240 : {
1241 : case GDT_Byte:
1242 1 : *pdfReal = pabySrc[iSrcOffset];
1243 1 : *pdfImag = 0.0;
1244 1 : break;
1245 :
1246 : case GDT_Int16:
1247 1 : *pdfReal = ((GInt16 *) pabySrc)[iSrcOffset];
1248 1 : *pdfImag = 0.0;
1249 1 : break;
1250 :
1251 : case GDT_UInt16:
1252 4 : *pdfReal = ((GUInt16 *) pabySrc)[iSrcOffset];
1253 4 : *pdfImag = 0.0;
1254 4 : break;
1255 :
1256 : case GDT_Int32:
1257 867 : *pdfReal = ((GInt32 *) pabySrc)[iSrcOffset];
1258 867 : *pdfImag = 0.0;
1259 867 : break;
1260 :
1261 : case GDT_UInt32:
1262 67 : *pdfReal = ((GUInt32 *) pabySrc)[iSrcOffset];
1263 67 : *pdfImag = 0.0;
1264 67 : break;
1265 :
1266 : case GDT_Float32:
1267 4 : *pdfReal = ((float *) pabySrc)[iSrcOffset];
1268 4 : *pdfImag = 0.0;
1269 4 : break;
1270 :
1271 : case GDT_Float64:
1272 67 : *pdfReal = ((double *) pabySrc)[iSrcOffset];
1273 67 : *pdfImag = 0.0;
1274 67 : break;
1275 :
1276 : case GDT_CInt16:
1277 67 : *pdfReal = ((GInt16 *) pabySrc)[iSrcOffset*2];
1278 67 : *pdfImag = ((GInt16 *) pabySrc)[iSrcOffset*2+1];
1279 67 : break;
1280 :
1281 : case GDT_CInt32:
1282 67 : *pdfReal = ((GInt32 *) pabySrc)[iSrcOffset*2];
1283 67 : *pdfImag = ((GInt32 *) pabySrc)[iSrcOffset*2+1];
1284 67 : break;
1285 :
1286 : case GDT_CFloat32:
1287 67 : *pdfReal = ((float *) pabySrc)[iSrcOffset*2];
1288 67 : *pdfImag = ((float *) pabySrc)[iSrcOffset*2+1];
1289 67 : break;
1290 :
1291 : case GDT_CFloat64:
1292 67 : *pdfReal = ((double *) pabySrc)[iSrcOffset*2];
1293 67 : *pdfImag = ((double *) pabySrc)[iSrcOffset*2+1];
1294 67 : break;
1295 :
1296 : default:
1297 0 : *pdfDensity = 0.0;
1298 0 : return FALSE;
1299 : }
1300 :
1301 1279 : if( poWK->pafUnifiedSrcDensity != NULL )
1302 0 : *pdfDensity = poWK->pafUnifiedSrcDensity[iSrcOffset];
1303 : else
1304 1279 : *pdfDensity = 1.0;
1305 :
1306 1279 : return *pdfDensity != 0.0;
1307 : }
1308 :
1309 : /************************************************************************/
1310 : /* GWKGetPixelRow() */
1311 : /************************************************************************/
1312 :
1313 : /* It is assumed that adfImag[] is set to 0 by caller code for non-complex */
1314 : /* data-types. */
1315 :
1316 1675450 : static int GWKGetPixelRow( GDALWarpKernel *poWK, int iBand,
1317 : int iSrcOffset, int nHalfSrcLen,
1318 : double* padfDensity,
1319 : double adfReal[],
1320 : double* padfImag )
1321 : {
1322 : // We know that nSrcLen is even, so we can *always* unroll loops 2x
1323 1675450 : int nSrcLen = nHalfSrcLen * 2;
1324 1675450 : int bHasValid = FALSE;
1325 : int i;
1326 :
1327 1675450 : if( padfDensity != NULL )
1328 : {
1329 : // Init the density
1330 35648 : for ( i = 0; i < nSrcLen; i += 2 )
1331 : {
1332 18112 : padfDensity[i] = 1.0;
1333 18112 : padfDensity[i+1] = 1.0;
1334 : }
1335 :
1336 17536 : if ( poWK->panUnifiedSrcValid != NULL )
1337 : {
1338 0 : for ( i = 0; i < nSrcLen; i += 2 )
1339 : {
1340 0 : if(poWK->panUnifiedSrcValid[(iSrcOffset+i)>>5]
1341 : & (0x01 << ((iSrcOffset+i) & 0x1f)))
1342 0 : bHasValid = TRUE;
1343 : else
1344 0 : padfDensity[i] = 0.0;
1345 :
1346 0 : if(poWK->panUnifiedSrcValid[(iSrcOffset+i+1)>>5]
1347 : & (0x01 << ((iSrcOffset+i+1) & 0x1f)))
1348 0 : bHasValid = TRUE;
1349 : else
1350 0 : padfDensity[i+1] = 0.0;
1351 : }
1352 :
1353 : // Reset or fail as needed
1354 0 : if ( bHasValid )
1355 0 : bHasValid = FALSE;
1356 : else
1357 0 : return FALSE;
1358 : }
1359 :
1360 32536 : if ( poWK->papanBandSrcValid != NULL
1361 15000 : && poWK->papanBandSrcValid[iBand] != NULL)
1362 : {
1363 30000 : for ( i = 0; i < nSrcLen; i += 2 )
1364 : {
1365 15000 : if(poWK->papanBandSrcValid[iBand][(iSrcOffset+i)>>5]
1366 : & (0x01 << ((iSrcOffset+i) & 0x1f)))
1367 15000 : bHasValid = TRUE;
1368 : else
1369 0 : padfDensity[i] = 0.0;
1370 :
1371 15000 : if(poWK->papanBandSrcValid[iBand][(iSrcOffset+i+1)>>5]
1372 : & (0x01 << ((iSrcOffset+i+1) & 0x1f)))
1373 15000 : bHasValid = TRUE;
1374 : else
1375 0 : padfDensity[i+1] = 0.0;
1376 : }
1377 :
1378 : // Reset or fail as needed
1379 15000 : if ( bHasValid )
1380 15000 : bHasValid = FALSE;
1381 : else
1382 0 : return FALSE;
1383 : }
1384 : }
1385 :
1386 : // Fetch data
1387 1675450 : switch( poWK->eWorkingDataType )
1388 : {
1389 : case GDT_Byte:
1390 : {
1391 1668406 : GByte* pSrc = (GByte*) poWK->papabySrcImage[iBand];
1392 1668406 : pSrc += iSrcOffset;
1393 6611511 : for ( i = 0; i < nSrcLen; i += 2 )
1394 : {
1395 4943105 : adfReal[i] = pSrc[i];
1396 4943105 : adfReal[i+1] = pSrc[i+1];
1397 : }
1398 1668406 : break;
1399 : }
1400 :
1401 : case GDT_Int16:
1402 : {
1403 294 : GInt16* pSrc = (GInt16*) poWK->papabySrcImage[iBand];
1404 294 : pSrc += iSrcOffset;
1405 1096 : for ( i = 0; i < nSrcLen; i += 2 )
1406 : {
1407 802 : adfReal[i] = pSrc[i];
1408 802 : adfReal[i+1] = pSrc[i+1];
1409 : }
1410 294 : break;
1411 : }
1412 :
1413 : case GDT_UInt16:
1414 : {
1415 750 : GUInt16* pSrc = (GUInt16*) poWK->papabySrcImage[iBand];
1416 750 : pSrc += iSrcOffset;
1417 2232 : for ( i = 0; i < nSrcLen; i += 2 )
1418 : {
1419 1482 : adfReal[i] = pSrc[i];
1420 1482 : adfReal[i+1] = pSrc[i+1];
1421 : }
1422 750 : break;
1423 : }
1424 :
1425 : case GDT_Int32:
1426 : {
1427 750 : GInt32* pSrc = (GInt32*) poWK->papabySrcImage[iBand];
1428 750 : pSrc += iSrcOffset;
1429 2232 : for ( i = 0; i < nSrcLen; i += 2 )
1430 : {
1431 1482 : adfReal[i] = pSrc[i];
1432 1482 : adfReal[i+1] = pSrc[i+1];
1433 : }
1434 750 : break;
1435 : }
1436 :
1437 : case GDT_UInt32:
1438 : {
1439 750 : GUInt32* pSrc = (GUInt32*) poWK->papabySrcImage[iBand];
1440 750 : pSrc += iSrcOffset;
1441 2232 : for ( i = 0; i < nSrcLen; i += 2 )
1442 : {
1443 1482 : adfReal[i] = pSrc[i];
1444 1482 : adfReal[i+1] = pSrc[i+1];
1445 : }
1446 750 : break;
1447 : }
1448 :
1449 : case GDT_Float32:
1450 : {
1451 750 : float* pSrc = (float*) poWK->papabySrcImage[iBand];
1452 750 : pSrc += iSrcOffset;
1453 2232 : for ( i = 0; i < nSrcLen; i += 2 )
1454 : {
1455 1482 : adfReal[i] = pSrc[i];
1456 1482 : adfReal[i+1] = pSrc[i+1];
1457 : }
1458 750 : break;
1459 : }
1460 :
1461 : case GDT_Float64:
1462 : {
1463 750 : double* pSrc = (double*) poWK->papabySrcImage[iBand];
1464 750 : pSrc += iSrcOffset;
1465 2232 : for ( i = 0; i < nSrcLen; i += 2 )
1466 : {
1467 1482 : adfReal[i] = pSrc[i];
1468 1482 : adfReal[i+1] = pSrc[i+1];
1469 : }
1470 750 : break;
1471 : }
1472 :
1473 : case GDT_CInt16:
1474 : {
1475 750 : GInt16* pSrc = (GInt16*) poWK->papabySrcImage[iBand];
1476 750 : pSrc += 2 * iSrcOffset;
1477 2232 : for ( i = 0; i < nSrcLen; i += 2 )
1478 : {
1479 1482 : adfReal[i] = pSrc[2*i];
1480 1482 : padfImag[i] = pSrc[2*i+1];
1481 :
1482 1482 : adfReal[i+1] = pSrc[2*i+2];
1483 1482 : padfImag[i+1] = pSrc[2*i+3];
1484 : }
1485 750 : break;
1486 : }
1487 :
1488 : case GDT_CInt32:
1489 : {
1490 750 : GInt32* pSrc = (GInt32*) poWK->papabySrcImage[iBand];
1491 750 : pSrc += 2 * iSrcOffset;
1492 2232 : for ( i = 0; i < nSrcLen; i += 2 )
1493 : {
1494 1482 : adfReal[i] = pSrc[2*i];
1495 1482 : padfImag[i] = pSrc[2*i+1];
1496 :
1497 1482 : adfReal[i+1] = pSrc[2*i+2];
1498 1482 : padfImag[i+1] = pSrc[2*i+3];
1499 : }
1500 750 : break;
1501 : }
1502 :
1503 : case GDT_CFloat32:
1504 : {
1505 750 : float* pSrc = (float*) poWK->papabySrcImage[iBand];
1506 750 : pSrc += 2 * iSrcOffset;
1507 2232 : for ( i = 0; i < nSrcLen; i += 2 )
1508 : {
1509 1482 : adfReal[i] = pSrc[2*i];
1510 1482 : padfImag[i] = pSrc[2*i+1];
1511 :
1512 1482 : adfReal[i+1] = pSrc[2*i+2];
1513 1482 : padfImag[i+1] = pSrc[2*i+3];
1514 : }
1515 750 : break;
1516 : }
1517 :
1518 :
1519 : case GDT_CFloat64:
1520 : {
1521 750 : double* pSrc = (double*) poWK->papabySrcImage[iBand];
1522 750 : pSrc += 2 * iSrcOffset;
1523 2232 : for ( i = 0; i < nSrcLen; i += 2 )
1524 : {
1525 1482 : adfReal[i] = pSrc[2*i];
1526 1482 : padfImag[i] = pSrc[2*i+1];
1527 :
1528 1482 : adfReal[i+1] = pSrc[2*i+2];
1529 1482 : padfImag[i+1] = pSrc[2*i+3];
1530 : }
1531 750 : break;
1532 : }
1533 :
1534 : default:
1535 0 : CPLAssert(FALSE);
1536 0 : if( padfDensity )
1537 0 : memset( padfDensity, 0, nSrcLen * sizeof(double) );
1538 0 : return FALSE;
1539 : }
1540 :
1541 1675450 : if( padfDensity == NULL )
1542 1657914 : return TRUE;
1543 :
1544 17536 : if( poWK->pafUnifiedSrcDensity == NULL )
1545 : {
1546 35184 : for ( i = 0; i < nSrcLen; i += 2 )
1547 : {
1548 : // Take into account earlier calcs
1549 17880 : if(padfDensity[i] > 0.000000001)
1550 : {
1551 17880 : padfDensity[i] = 1.0;
1552 17880 : bHasValid = TRUE;
1553 : }
1554 :
1555 17880 : if(padfDensity[i+1] > 0.000000001)
1556 : {
1557 17880 : padfDensity[i+1] = 1.0;
1558 17880 : bHasValid = TRUE;
1559 : }
1560 : }
1561 : }
1562 : else
1563 : {
1564 464 : for ( i = 0; i < nSrcLen; i += 2 )
1565 : {
1566 232 : if(padfDensity[i] > 0.000000001)
1567 232 : padfDensity[i] = poWK->pafUnifiedSrcDensity[iSrcOffset+i];
1568 232 : if(padfDensity[i] > 0.000000001)
1569 176 : bHasValid = TRUE;
1570 :
1571 232 : if(padfDensity[i+1] > 0.000000001)
1572 232 : padfDensity[i+1] = poWK->pafUnifiedSrcDensity[iSrcOffset+i+1];
1573 232 : if(padfDensity[i+1] > 0.000000001)
1574 212 : bHasValid = TRUE;
1575 : }
1576 : }
1577 :
1578 17536 : return bHasValid;
1579 : }
1580 :
1581 : /************************************************************************/
1582 : /* GWKGetPixelByte() */
1583 : /************************************************************************/
1584 :
1585 725901 : static int GWKGetPixelByte( GDALWarpKernel *poWK, int iBand,
1586 : int iSrcOffset, double *pdfDensity,
1587 : GByte *pbValue )
1588 :
1589 : {
1590 725901 : GByte *pabySrc = poWK->papabySrcImage[iBand];
1591 :
1592 1101739 : if ( ( poWK->panUnifiedSrcValid != NULL
1593 375838 : && !((poWK->panUnifiedSrcValid[iSrcOffset>>5]
1594 : & (0x01 << (iSrcOffset & 0x1f))) ) )
1595 : || ( poWK->papanBandSrcValid != NULL
1596 0 : && poWK->papanBandSrcValid[iBand] != NULL
1597 0 : && !((poWK->papanBandSrcValid[iBand][iSrcOffset>>5]
1598 : & (0x01 << (iSrcOffset & 0x1f)))) ) )
1599 : {
1600 0 : *pdfDensity = 0.0;
1601 0 : return FALSE;
1602 : }
1603 :
1604 725901 : *pbValue = pabySrc[iSrcOffset];
1605 :
1606 725901 : if ( poWK->pafUnifiedSrcDensity == NULL )
1607 710798 : *pdfDensity = 1.0;
1608 : else
1609 15103 : *pdfDensity = poWK->pafUnifiedSrcDensity[iSrcOffset];
1610 :
1611 725901 : return *pdfDensity != 0.0;
1612 : }
1613 :
1614 : /************************************************************************/
1615 : /* GWKGetPixelShort() */
1616 : /************************************************************************/
1617 :
1618 1488467 : static int GWKGetPixelShort( GDALWarpKernel *poWK, int iBand,
1619 : int iSrcOffset, double *pdfDensity,
1620 : GInt16 *piValue )
1621 :
1622 : {
1623 1488467 : GInt16 *pabySrc = (GInt16 *)poWK->papabySrcImage[iBand];
1624 :
1625 4462995 : if ( ( poWK->panUnifiedSrcValid != NULL
1626 0 : && !((poWK->panUnifiedSrcValid[iSrcOffset>>5]
1627 : & (0x01 << (iSrcOffset & 0x1f))) ) )
1628 : || ( poWK->papanBandSrcValid != NULL
1629 1487264 : && poWK->papanBandSrcValid[iBand] != NULL
1630 1487264 : && !((poWK->papanBandSrcValid[iBand][iSrcOffset>>5]
1631 : & (0x01 << (iSrcOffset & 0x1f)))) ) )
1632 : {
1633 715 : *pdfDensity = 0.0;
1634 715 : return FALSE;
1635 : }
1636 :
1637 1487752 : *piValue = pabySrc[iSrcOffset];
1638 :
1639 1487752 : if ( poWK->pafUnifiedSrcDensity == NULL )
1640 1486549 : *pdfDensity = 1.0;
1641 : else
1642 1203 : *pdfDensity = poWK->pafUnifiedSrcDensity[iSrcOffset];
1643 :
1644 1487752 : return *pdfDensity != 0.0;
1645 : }
1646 :
1647 : /************************************************************************/
1648 : /* GWKGetPixelFloat() */
1649 : /************************************************************************/
1650 :
1651 1603 : static int GWKGetPixelFloat( GDALWarpKernel *poWK, int iBand,
1652 : int iSrcOffset, double *pdfDensity,
1653 : float *pfValue )
1654 :
1655 : {
1656 1603 : float *pabySrc = (float *)poWK->papabySrcImage[iBand];
1657 :
1658 2403 : if ( ( poWK->panUnifiedSrcValid != NULL
1659 0 : && !((poWK->panUnifiedSrcValid[iSrcOffset>>5]
1660 : & (0x01 << (iSrcOffset & 0x1f))) ) )
1661 : || ( poWK->papanBandSrcValid != NULL
1662 400 : && poWK->papanBandSrcValid[iBand] != NULL
1663 400 : && !((poWK->papanBandSrcValid[iBand][iSrcOffset>>5]
1664 : & (0x01 << (iSrcOffset & 0x1f)))) ) )
1665 : {
1666 40 : *pdfDensity = 0.0;
1667 40 : return FALSE;
1668 : }
1669 :
1670 1563 : *pfValue = pabySrc[iSrcOffset];
1671 :
1672 1563 : if ( poWK->pafUnifiedSrcDensity == NULL )
1673 360 : *pdfDensity = 1.0;
1674 : else
1675 1203 : *pdfDensity = poWK->pafUnifiedSrcDensity[iSrcOffset];
1676 :
1677 1563 : return *pdfDensity != 0.0;
1678 : }
1679 :
1680 : /************************************************************************/
1681 : /* GWKBilinearResample() */
1682 : /* Set of bilinear interpolators */
1683 : /************************************************************************/
1684 :
1685 8588 : static int GWKBilinearResample( GDALWarpKernel *poWK, int iBand,
1686 : double dfSrcX, double dfSrcY,
1687 : double *pdfDensity,
1688 : double *pdfReal, double *pdfImag )
1689 :
1690 : {
1691 : // Save as local variables to avoid following pointers
1692 8588 : int nSrcXSize = poWK->nSrcXSize;
1693 8588 : int nSrcYSize = poWK->nSrcYSize;
1694 :
1695 8588 : int iSrcX = (int) floor(dfSrcX - 0.5);
1696 8588 : int iSrcY = (int) floor(dfSrcY - 0.5);
1697 : int iSrcOffset;
1698 8588 : double dfRatioX = 1.5 - (dfSrcX - iSrcX);
1699 8588 : double dfRatioY = 1.5 - (dfSrcY - iSrcY);
1700 8588 : double adfDensity[2], adfReal[2], adfImag[2] = {0, 0};
1701 8588 : double dfAccumulatorReal = 0.0, dfAccumulatorImag = 0.0;
1702 8588 : double dfAccumulatorDensity = 0.0;
1703 8588 : double dfAccumulatorDivisor = 0.0;
1704 8588 : int bShifted = FALSE;
1705 :
1706 8588 : if (iSrcX == -1)
1707 : {
1708 0 : iSrcX = 0;
1709 0 : dfRatioX = 1;
1710 : }
1711 8588 : if (iSrcY == -1)
1712 : {
1713 150 : iSrcY = 0;
1714 150 : dfRatioY = 1;
1715 : }
1716 8588 : iSrcOffset = iSrcX + iSrcY * nSrcXSize;
1717 :
1718 : // Shift so we don't overrun the array
1719 8588 : if( nSrcXSize * nSrcYSize == iSrcOffset + 1
1720 : || nSrcXSize * nSrcYSize == iSrcOffset + nSrcXSize + 1 )
1721 : {
1722 111 : bShifted = TRUE;
1723 111 : --iSrcOffset;
1724 : }
1725 :
1726 : // Get pixel row
1727 8588 : if ( iSrcY >= 0 && iSrcY < nSrcYSize
1728 : && iSrcOffset >= 0 && iSrcOffset < nSrcXSize * nSrcYSize
1729 : && GWKGetPixelRow( poWK, iBand, iSrcOffset, 1,
1730 : adfDensity, adfReal, adfImag ) )
1731 : {
1732 8568 : double dfMult1 = dfRatioX * dfRatioY;
1733 8568 : double dfMult2 = (1.0-dfRatioX) * dfRatioY;
1734 :
1735 : // Shifting corrected
1736 8568 : if ( bShifted )
1737 : {
1738 111 : adfReal[0] = adfReal[1];
1739 111 : adfImag[0] = adfImag[1];
1740 111 : adfDensity[0] = adfDensity[1];
1741 : }
1742 :
1743 : // Upper Left Pixel
1744 17136 : if ( iSrcX >= 0 && iSrcX < nSrcXSize
1745 8568 : && adfDensity[0] > 0.000000001 )
1746 : {
1747 8552 : dfAccumulatorDivisor += dfMult1;
1748 :
1749 8552 : dfAccumulatorReal += adfReal[0] * dfMult1;
1750 8552 : dfAccumulatorImag += adfImag[0] * dfMult1;
1751 8552 : dfAccumulatorDensity += adfDensity[0] * dfMult1;
1752 : }
1753 :
1754 : // Upper Right Pixel
1755 16770 : if ( iSrcX+1 >= 0 && iSrcX+1 < nSrcXSize
1756 8202 : && adfDensity[1] > 0.000000001 )
1757 : {
1758 8202 : dfAccumulatorDivisor += dfMult2;
1759 :
1760 8202 : dfAccumulatorReal += adfReal[1] * dfMult2;
1761 8202 : dfAccumulatorImag += adfImag[1] * dfMult2;
1762 8202 : dfAccumulatorDensity += adfDensity[1] * dfMult2;
1763 : }
1764 : }
1765 :
1766 : // Get pixel row
1767 8588 : if ( iSrcY+1 >= 0 && iSrcY+1 < nSrcYSize
1768 : && iSrcOffset+nSrcXSize >= 0
1769 : && iSrcOffset+nSrcXSize < nSrcXSize * nSrcYSize
1770 : && GWKGetPixelRow( poWK, iBand, iSrcOffset+nSrcXSize, 1,
1771 : adfDensity, adfReal, adfImag ) )
1772 : {
1773 8372 : double dfMult1 = dfRatioX * (1.0-dfRatioY);
1774 8372 : double dfMult2 = (1.0-dfRatioX) * (1.0-dfRatioY);
1775 :
1776 : // Shifting corrected
1777 8372 : if ( bShifted )
1778 : {
1779 57 : adfReal[0] = adfReal[1];
1780 57 : adfImag[0] = adfImag[1];
1781 57 : adfDensity[0] = adfDensity[1];
1782 : }
1783 :
1784 : // Lower Left Pixel
1785 16744 : if ( iSrcX >= 0 && iSrcX < nSrcXSize
1786 8372 : && adfDensity[0] > 0.000000001 )
1787 : {
1788 8352 : dfAccumulatorDivisor += dfMult1;
1789 :
1790 8352 : dfAccumulatorReal += adfReal[0] * dfMult1;
1791 8352 : dfAccumulatorImag += adfImag[0] * dfMult1;
1792 8352 : dfAccumulatorDensity += adfDensity[0] * dfMult1;
1793 : }
1794 :
1795 : // Lower Right Pixel
1796 16432 : if ( iSrcX+1 >= 0 && iSrcX+1 < nSrcXSize
1797 8060 : && adfDensity[1] > 0.000000001 )
1798 : {
1799 8060 : dfAccumulatorDivisor += dfMult2;
1800 :
1801 8060 : dfAccumulatorReal += adfReal[1] * dfMult2;
1802 8060 : dfAccumulatorImag += adfImag[1] * dfMult2;
1803 8060 : dfAccumulatorDensity += adfDensity[1] * dfMult2;
1804 : }
1805 : }
1806 :
1807 : /* -------------------------------------------------------------------- */
1808 : /* Return result. */
1809 : /* -------------------------------------------------------------------- */
1810 8588 : if ( dfAccumulatorDivisor == 1.0 )
1811 : {
1812 8402 : *pdfReal = dfAccumulatorReal;
1813 8402 : *pdfImag = dfAccumulatorImag;
1814 8402 : *pdfDensity = dfAccumulatorDensity;
1815 8402 : return TRUE;
1816 : }
1817 186 : else if ( dfAccumulatorDivisor < 0.00001 )
1818 : {
1819 0 : *pdfReal = 0.0;
1820 0 : *pdfImag = 0.0;
1821 0 : *pdfDensity = 0.0;
1822 0 : return FALSE;
1823 : }
1824 : else
1825 : {
1826 186 : *pdfReal = dfAccumulatorReal / dfAccumulatorDivisor;
1827 186 : *pdfImag = dfAccumulatorImag / dfAccumulatorDivisor;
1828 186 : *pdfDensity = dfAccumulatorDensity / dfAccumulatorDivisor;
1829 186 : return TRUE;
1830 : }
1831 : }
1832 :
1833 289915 : static int GWKBilinearResampleNoMasksByte( GDALWarpKernel *poWK, int iBand,
1834 : double dfSrcX, double dfSrcY,
1835 : GByte *pbValue )
1836 :
1837 : {
1838 289915 : double dfAccumulator = 0.0;
1839 289915 : double dfAccumulatorDivisor = 0.0;
1840 :
1841 289915 : int iSrcX = (int) floor(dfSrcX - 0.5);
1842 289915 : int iSrcY = (int) floor(dfSrcY - 0.5);
1843 289915 : int iSrcOffset = iSrcX + iSrcY * poWK->nSrcXSize;
1844 289915 : double dfRatioX = 1.5 - (dfSrcX - iSrcX);
1845 289915 : double dfRatioY = 1.5 - (dfSrcY - iSrcY);
1846 :
1847 : // Upper Left Pixel
1848 289915 : if( iSrcX >= 0 && iSrcX < poWK->nSrcXSize
1849 : && iSrcY >= 0 && iSrcY < poWK->nSrcYSize )
1850 : {
1851 283063 : double dfMult = dfRatioX * dfRatioY;
1852 :
1853 283063 : dfAccumulatorDivisor += dfMult;
1854 :
1855 : dfAccumulator +=
1856 283063 : (double)poWK->papabySrcImage[iBand][iSrcOffset] * dfMult;
1857 : }
1858 :
1859 : // Upper Right Pixel
1860 289915 : if( iSrcX+1 >= 0 && iSrcX+1 < poWK->nSrcXSize
1861 : && iSrcY >= 0 && iSrcY < poWK->nSrcYSize )
1862 : {
1863 286086 : double dfMult = (1.0-dfRatioX) * dfRatioY;
1864 :
1865 286086 : dfAccumulatorDivisor += dfMult;
1866 :
1867 : dfAccumulator +=
1868 286086 : (double)poWK->papabySrcImage[iBand][iSrcOffset+1] * dfMult;
1869 : }
1870 :
1871 : // Lower Right Pixel
1872 289915 : if( iSrcX+1 >= 0 && iSrcX+1 < poWK->nSrcXSize
1873 : && iSrcY+1 >= 0 && iSrcY+1 < poWK->nSrcYSize )
1874 : {
1875 289188 : double dfMult = (1.0-dfRatioX) * (1.0-dfRatioY);
1876 :
1877 289188 : dfAccumulatorDivisor += dfMult;
1878 :
1879 : dfAccumulator +=
1880 289188 : (double)poWK->papabySrcImage[iBand][iSrcOffset+1+poWK->nSrcXSize]
1881 289188 : * dfMult;
1882 : }
1883 :
1884 : // Lower Left Pixel
1885 289915 : if( iSrcX >= 0 && iSrcX < poWK->nSrcXSize
1886 : && iSrcY+1 >= 0 && iSrcY+1 < poWK->nSrcYSize )
1887 : {
1888 286134 : double dfMult = dfRatioX * (1.0-dfRatioY);
1889 :
1890 286134 : dfAccumulatorDivisor += dfMult;
1891 :
1892 : dfAccumulator +=
1893 286134 : (double)poWK->papabySrcImage[iBand][iSrcOffset+poWK->nSrcXSize]
1894 286134 : * dfMult;
1895 : }
1896 :
1897 : /* -------------------------------------------------------------------- */
1898 : /* Return result. */
1899 : /* -------------------------------------------------------------------- */
1900 : double dfValue;
1901 :
1902 289915 : if( dfAccumulatorDivisor < 0.00001 )
1903 : {
1904 0 : *pbValue = 0;
1905 0 : return FALSE;
1906 : }
1907 289915 : else if( dfAccumulatorDivisor == 1.0 )
1908 : {
1909 274559 : dfValue = dfAccumulator;
1910 : }
1911 : else
1912 : {
1913 15356 : dfValue = dfAccumulator / dfAccumulatorDivisor;
1914 : }
1915 :
1916 289915 : if ( dfValue < 0.0 )
1917 0 : *pbValue = 0;
1918 289915 : else if ( dfValue > 255.0 )
1919 716 : *pbValue = 255;
1920 : else
1921 289199 : *pbValue = (GByte)(0.5 + dfValue);
1922 :
1923 289915 : return TRUE;
1924 : }
1925 :
1926 1714891 : static int GWKBilinearResampleNoMasksShort( GDALWarpKernel *poWK, int iBand,
1927 : double dfSrcX, double dfSrcY,
1928 : GInt16 *piValue )
1929 :
1930 : {
1931 1714891 : double dfAccumulator = 0.0;
1932 1714891 : double dfAccumulatorDivisor = 0.0;
1933 :
1934 1714891 : int iSrcX = (int) floor(dfSrcX - 0.5);
1935 1714891 : int iSrcY = (int) floor(dfSrcY - 0.5);
1936 1714891 : int iSrcOffset = iSrcX + iSrcY * poWK->nSrcXSize;
1937 1714891 : double dfRatioX = 1.5 - (dfSrcX - iSrcX);
1938 1714891 : double dfRatioY = 1.5 - (dfSrcY - iSrcY);
1939 :
1940 : // Upper Left Pixel
1941 1714891 : if( iSrcX >= 0 && iSrcX < poWK->nSrcXSize
1942 : && iSrcY >= 0 && iSrcY < poWK->nSrcYSize )
1943 : {
1944 1707355 : double dfMult = dfRatioX * dfRatioY;
1945 :
1946 1707355 : dfAccumulatorDivisor += dfMult;
1947 :
1948 : dfAccumulator +=
1949 1707355 : (double)((GInt16 *)poWK->papabySrcImage[iBand])[iSrcOffset]
1950 1707355 : * dfMult;
1951 : }
1952 :
1953 : // Upper Right Pixel
1954 1714891 : if( iSrcX+1 >= 0 && iSrcX+1 < poWK->nSrcXSize
1955 : && iSrcY >= 0 && iSrcY < poWK->nSrcYSize )
1956 : {
1957 1709182 : double dfMult = (1.0-dfRatioX) * dfRatioY;
1958 :
1959 1709182 : dfAccumulatorDivisor += dfMult;
1960 :
1961 : dfAccumulator +=
1962 1709182 : (double)((GInt16 *)poWK->papabySrcImage[iBand])[iSrcOffset+1] * dfMult;
1963 : }
1964 :
1965 : // Lower Right Pixel
1966 1714891 : if( iSrcX+1 >= 0 && iSrcX+1 < poWK->nSrcXSize
1967 : && iSrcY+1 >= 0 && iSrcY+1 < poWK->nSrcYSize )
1968 : {
1969 1713535 : double dfMult = (1.0-dfRatioX) * (1.0-dfRatioY);
1970 :
1971 1713535 : dfAccumulatorDivisor += dfMult;
1972 :
1973 : dfAccumulator +=
1974 1713535 : (double)((GInt16 *)poWK->papabySrcImage[iBand])[iSrcOffset+1+poWK->nSrcXSize]
1975 1713535 : * dfMult;
1976 : }
1977 :
1978 : // Lower Left Pixel
1979 1714891 : if( iSrcX >= 0 && iSrcX < poWK->nSrcXSize
1980 : && iSrcY+1 >= 0 && iSrcY+1 < poWK->nSrcYSize )
1981 : {
1982 1711682 : double dfMult = dfRatioX * (1.0-dfRatioY);
1983 :
1984 1711682 : dfAccumulatorDivisor += dfMult;
1985 :
1986 : dfAccumulator +=
1987 1711682 : (double)((GInt16 *)poWK->papabySrcImage[iBand])[iSrcOffset+poWK->nSrcXSize]
1988 1711682 : * dfMult;
1989 : }
1990 :
1991 : /* -------------------------------------------------------------------- */
1992 : /* Return result. */
1993 : /* -------------------------------------------------------------------- */
1994 1714891 : if( dfAccumulatorDivisor == 1.0 )
1995 : {
1996 1698761 : *piValue = (GInt16)(0.5 + dfAccumulator);
1997 1698761 : return TRUE;
1998 : }
1999 16130 : else if( dfAccumulatorDivisor < 0.00001 )
2000 : {
2001 0 : *piValue = 0;
2002 0 : return FALSE;
2003 : }
2004 : else
2005 : {
2006 16130 : *piValue = (GInt16)(0.5 + dfAccumulator / dfAccumulatorDivisor);
2007 16130 : return TRUE;
2008 : }
2009 : }
2010 :
2011 : /************************************************************************/
2012 : /* GWKCubicResample() */
2013 : /* Set of bicubic interpolators using cubic convolution. */
2014 : /************************************************************************/
2015 :
2016 : #define CubicConvolution(distance1,distance2,distance3,f0,f1,f2,f3) \
2017 : ( ( -f0 + f1 - f2 + f3) * distance3 \
2018 : + (2.0*(f0 - f1) + f2 - f3) * distance2 \
2019 : + ( -f0 + f2 ) * distance1 \
2020 : + f1 )
2021 :
2022 558 : static int GWKCubicResample( GDALWarpKernel *poWK, int iBand,
2023 : double dfSrcX, double dfSrcY,
2024 : double *pdfDensity,
2025 : double *pdfReal, double *pdfImag )
2026 :
2027 : {
2028 558 : int iSrcX = (int) floor( dfSrcX - 0.5 );
2029 558 : int iSrcY = (int) floor( dfSrcY - 0.5 );
2030 558 : int iSrcOffset = iSrcX + iSrcY * poWK->nSrcXSize;
2031 558 : double dfDeltaX = dfSrcX - 0.5 - iSrcX;
2032 558 : double dfDeltaY = dfSrcY - 0.5 - iSrcY;
2033 558 : double dfDeltaX2 = dfDeltaX * dfDeltaX;
2034 558 : double dfDeltaY2 = dfDeltaY * dfDeltaY;
2035 558 : double dfDeltaX3 = dfDeltaX2 * dfDeltaX;
2036 558 : double dfDeltaY3 = dfDeltaY2 * dfDeltaY;
2037 : double adfValueDens[4], adfValueReal[4], adfValueImag[4];
2038 558 : double adfDensity[4], adfReal[4], adfImag[4] = {0, 0, 0, 0};
2039 : int i;
2040 :
2041 : // Get the bilinear interpolation at the image borders
2042 558 : if ( iSrcX - 1 < 0 || iSrcX + 2 >= poWK->nSrcXSize
2043 : || iSrcY - 1 < 0 || iSrcY + 2 >= poWK->nSrcYSize )
2044 : return GWKBilinearResample( poWK, iBand, dfSrcX, dfSrcY,
2045 414 : pdfDensity, pdfReal, pdfImag );
2046 :
2047 720 : for ( i = -1; i < 3; i++ )
2048 : {
2049 2880 : if ( !GWKGetPixelRow(poWK, iBand, iSrcOffset + i * poWK->nSrcXSize - 1,
2050 : 2, adfDensity, adfReal, adfImag)
2051 576 : || adfDensity[0] < 0.000000001
2052 576 : || adfDensity[1] < 0.000000001
2053 576 : || adfDensity[2] < 0.000000001
2054 576 : || adfDensity[3] < 0.000000001 )
2055 : {
2056 : return GWKBilinearResample( poWK, iBand, dfSrcX, dfSrcY,
2057 0 : pdfDensity, pdfReal, pdfImag );
2058 : }
2059 :
2060 4608 : adfValueDens[i + 1] = CubicConvolution(dfDeltaX, dfDeltaX2, dfDeltaX3,
2061 4608 : adfDensity[0], adfDensity[1], adfDensity[2], adfDensity[3]);
2062 4608 : adfValueReal[i + 1] = CubicConvolution(dfDeltaX, dfDeltaX2, dfDeltaX3,
2063 4608 : adfReal[0], adfReal[1], adfReal[2], adfReal[3]);
2064 4608 : adfValueImag[i + 1] = CubicConvolution(dfDeltaX, dfDeltaX2, dfDeltaX3,
2065 4608 : adfImag[0], adfImag[1], adfImag[2], adfImag[3]);
2066 : }
2067 :
2068 :
2069 : /* -------------------------------------------------------------------- */
2070 : /* For now, if we have any pixels missing in the kernel area, */
2071 : /* we fallback on using bilinear interpolation. Ideally we */
2072 : /* should do "weight adjustment" of our results similarly to */
2073 : /* what is done for the cubic spline and lanc. interpolators. */
2074 : /* -------------------------------------------------------------------- */
2075 1152 : *pdfDensity = CubicConvolution(dfDeltaY, dfDeltaY2, dfDeltaY3,
2076 : adfValueDens[0], adfValueDens[1],
2077 1152 : adfValueDens[2], adfValueDens[3]);
2078 1152 : *pdfReal = CubicConvolution(dfDeltaY, dfDeltaY2, dfDeltaY3,
2079 : adfValueReal[0], adfValueReal[1],
2080 1152 : adfValueReal[2], adfValueReal[3]);
2081 1152 : *pdfImag = CubicConvolution(dfDeltaY, dfDeltaY2, dfDeltaY3,
2082 : adfValueImag[0], adfValueImag[1],
2083 1152 : adfValueImag[2], adfValueImag[3]);
2084 :
2085 144 : return TRUE;
2086 : }
2087 :
2088 278207 : static int GWKCubicResampleNoMasksByte( GDALWarpKernel *poWK, int iBand,
2089 : double dfSrcX, double dfSrcY,
2090 : GByte *pbValue )
2091 :
2092 : {
2093 278207 : int iSrcX = (int) floor( dfSrcX - 0.5 );
2094 278207 : int iSrcY = (int) floor( dfSrcY - 0.5 );
2095 278207 : int iSrcOffset = iSrcX + iSrcY * poWK->nSrcXSize;
2096 278207 : double dfDeltaX = dfSrcX - 0.5 - iSrcX;
2097 278207 : double dfDeltaY = dfSrcY - 0.5 - iSrcY;
2098 278207 : double dfDeltaX2 = dfDeltaX * dfDeltaX;
2099 278207 : double dfDeltaY2 = dfDeltaY * dfDeltaY;
2100 278207 : double dfDeltaX3 = dfDeltaX2 * dfDeltaX;
2101 278207 : double dfDeltaY3 = dfDeltaY2 * dfDeltaY;
2102 : double adfValue[4];
2103 : int i;
2104 :
2105 : // Get the bilinear interpolation at the image borders
2106 278207 : if ( iSrcX - 1 < 0 || iSrcX + 2 >= poWK->nSrcXSize
2107 : || iSrcY - 1 < 0 || iSrcY + 2 >= poWK->nSrcYSize )
2108 : return GWKBilinearResampleNoMasksByte( poWK, iBand, dfSrcX, dfSrcY,
2109 10570 : pbValue);
2110 :
2111 1338185 : for ( i = -1; i < 3; i++ )
2112 : {
2113 1070548 : int iOffset = iSrcOffset + i * poWK->nSrcXSize;
2114 :
2115 11776028 : adfValue[i + 1] = CubicConvolution(dfDeltaX, dfDeltaX2, dfDeltaX3,
2116 : (double)poWK->papabySrcImage[iBand][iOffset - 1],
2117 : (double)poWK->papabySrcImage[iBand][iOffset],
2118 : (double)poWK->papabySrcImage[iBand][iOffset + 1],
2119 11776028 : (double)poWK->papabySrcImage[iBand][iOffset + 2]);
2120 : }
2121 :
2122 267637 : double dfValue = CubicConvolution(dfDeltaY, dfDeltaY2, dfDeltaY3,
2123 : adfValue[0], adfValue[1], adfValue[2], adfValue[3]);
2124 :
2125 267637 : if ( dfValue < 0.0 )
2126 8 : *pbValue = 0;
2127 267629 : else if ( dfValue > 255.0 )
2128 11506 : *pbValue = 255;
2129 : else
2130 256123 : *pbValue = (GByte)(0.5 + dfValue);
2131 :
2132 267637 : return TRUE;
2133 : }
2134 :
2135 262207 : static int GWKCubicResampleNoMasksShort( GDALWarpKernel *poWK, int iBand,
2136 : double dfSrcX, double dfSrcY,
2137 : GInt16 *piValue )
2138 :
2139 : {
2140 262207 : int iSrcX = (int) floor( dfSrcX - 0.5 );
2141 262207 : int iSrcY = (int) floor( dfSrcY - 0.5 );
2142 262207 : int iSrcOffset = iSrcX + iSrcY * poWK->nSrcXSize;
2143 262207 : double dfDeltaX = dfSrcX - 0.5 - iSrcX;
2144 262207 : double dfDeltaY = dfSrcY - 0.5 - iSrcY;
2145 262207 : double dfDeltaX2 = dfDeltaX * dfDeltaX;
2146 262207 : double dfDeltaY2 = dfDeltaY * dfDeltaY;
2147 262207 : double dfDeltaX3 = dfDeltaX2 * dfDeltaX;
2148 262207 : double dfDeltaY3 = dfDeltaY2 * dfDeltaY;
2149 : double adfValue[4];
2150 : int i;
2151 :
2152 : // Get the bilinear interpolation at the image borders
2153 262207 : if ( iSrcX - 1 < 0 || iSrcX + 2 >= poWK->nSrcXSize
2154 : || iSrcY - 1 < 0 || iSrcY + 2 >= poWK->nSrcYSize )
2155 : return GWKBilinearResampleNoMasksShort( poWK, iBand, dfSrcX, dfSrcY,
2156 9182 : piValue);
2157 :
2158 1265125 : for ( i = -1; i < 3; i++ )
2159 : {
2160 1012100 : int iOffset = iSrcOffset + i * poWK->nSrcXSize;
2161 :
2162 11133100 : adfValue[i + 1] =CubicConvolution(dfDeltaX, dfDeltaX2, dfDeltaX3,
2163 : (double)((GInt16 *)poWK->papabySrcImage[iBand])[iOffset - 1],
2164 : (double)((GInt16 *)poWK->papabySrcImage[iBand])[iOffset],
2165 : (double)((GInt16 *)poWK->papabySrcImage[iBand])[iOffset + 1],
2166 11133100 : (double)((GInt16 *)poWK->papabySrcImage[iBand])[iOffset + 2]);
2167 : }
2168 :
2169 2024200 : *piValue = (GInt16)CubicConvolution(dfDeltaY, dfDeltaY2, dfDeltaY3,
2170 2024200 : adfValue[0], adfValue[1], adfValue[2], adfValue[3]);
2171 :
2172 253025 : return TRUE;
2173 : }
2174 :
2175 :
2176 : /************************************************************************/
2177 : /* GWKLanczosSinc() */
2178 : /************************************************************************/
2179 :
2180 : /*
2181 : * Lanczos windowed sinc interpolation kernel with radius r.
2182 : * /
2183 : * | sinc(x) * sinc(x/r), if |x| < r
2184 : * L(x) = | 1, if x = 0 ,
2185 : * | 0, otherwise
2186 : * \
2187 : *
2188 : * where sinc(x) = sin(PI * x) / (PI * x).
2189 : */
2190 :
2191 : #define GWK_PI 3.14159265358979323846
2192 :
2193 0 : static double GWKLanczosSinc( double dfX, double dfR )
2194 : {
2195 0 : if ( dfX == 0.0 )
2196 0 : return 1.0;
2197 :
2198 0 : const double dfPIX = GWK_PI * dfX;
2199 0 : const double dfPIXoverR = dfPIX / dfR;
2200 0 : const double dfPIX2overR = dfPIX * dfPIXoverR;
2201 0 : return sin(dfPIX) * sin(dfPIXoverR) / dfPIX2overR;
2202 : }
2203 :
2204 : //#undef GWK_PI
2205 :
2206 : /************************************************************************/
2207 : /* GWKBSpline() */
2208 : /************************************************************************/
2209 :
2210 4326896 : static double GWKBSpline( double x )
2211 : {
2212 4326896 : double xp2 = x + 2.0;
2213 4326896 : double xp1 = x + 1.0;
2214 4326896 : double xm1 = x - 1.0;
2215 :
2216 : // This will most likely be used, so we'll compute it ahead of time to
2217 : // avoid stalling the processor
2218 4326896 : double xp2c = xp2 * xp2 * xp2;
2219 :
2220 : // Note that the test is computed only if it is needed
2221 : return (((xp2 > 0.0)?((xp1 > 0.0)?((x > 0.0)?((xm1 > 0.0)?
2222 : -4.0 * xm1*xm1*xm1:0.0) +
2223 : 6.0 * x*x*x:0.0) +
2224 : -4.0 * xp1*xp1*xp1:0.0) +
2225 4326896 : xp2c:0.0) ) * 0.166666666666666666666;
2226 : }
2227 :
2228 :
2229 : /************************************************************************/
2230 : /* GWKResampleWrkStruct */
2231 : /************************************************************************/
2232 :
2233 : typedef struct _GWKResampleWrkStruct GWKResampleWrkStruct;
2234 :
2235 : typedef int (*pfnGWKResampleType) ( GDALWarpKernel *poWK, int iBand,
2236 : double dfSrcX, double dfSrcY,
2237 : double *pdfDensity,
2238 : double *pdfReal, double *pdfImag,
2239 : GWKResampleWrkStruct* psWrkStruct );
2240 :
2241 :
2242 : struct _GWKResampleWrkStruct
2243 : {
2244 : pfnGWKResampleType pfnGWKResample;
2245 :
2246 : // Space for saved X weights
2247 : double *padfWeightsX;
2248 : char *panCalcX;
2249 :
2250 : double *padfWeightsY; // only used by GWKResampleOptimizedLanczos
2251 : int iLastSrcX; // only used by GWKResampleOptimizedLanczos
2252 : int iLastSrcY; // only used by GWKResampleOptimizedLanczos
2253 : double dfLastDeltaX; // only used by GWKResampleOptimizedLanczos
2254 : double dfLastDeltaY; // only used by GWKResampleOptimizedLanczos
2255 :
2256 : // Space for saving a row of pixels
2257 : double *padfRowDensity;
2258 : double *padfRowReal;
2259 : double *padfRowImag;
2260 : };
2261 :
2262 : /************************************************************************/
2263 : /* GWKResampleCreateWrkStruct() */
2264 : /************************************************************************/
2265 :
2266 : static int GWKResample( GDALWarpKernel *poWK, int iBand,
2267 : double dfSrcX, double dfSrcY,
2268 : double *pdfDensity,
2269 : double *pdfReal, double *pdfImag,
2270 : GWKResampleWrkStruct* psWrkStruct );
2271 :
2272 : static int GWKResampleOptimizedLanczos( GDALWarpKernel *poWK, int iBand,
2273 : double dfSrcX, double dfSrcY,
2274 : double *pdfDensity,
2275 : double *pdfReal, double *pdfImag,
2276 : GWKResampleWrkStruct* psWrkStruct );
2277 :
2278 87 : static GWKResampleWrkStruct* GWKResampleCreateWrkStruct(GDALWarpKernel *poWK)
2279 : {
2280 87 : int nXDist = ( poWK->nXRadius + 1 ) * 2;
2281 :
2282 : GWKResampleWrkStruct* psWrkStruct =
2283 87 : (GWKResampleWrkStruct*)CPLMalloc(sizeof(GWKResampleWrkStruct));
2284 :
2285 : // Alloc space for saved X weights
2286 87 : psWrkStruct->padfWeightsX = (double *)CPLCalloc( nXDist, sizeof(double) );
2287 87 : psWrkStruct->panCalcX = (char *)CPLMalloc( nXDist * sizeof(char) );
2288 :
2289 87 : psWrkStruct->padfWeightsY = (double *)CPLCalloc( nXDist, sizeof(double) );
2290 87 : psWrkStruct->iLastSrcX = -10;
2291 87 : psWrkStruct->iLastSrcY = -10;
2292 87 : psWrkStruct->dfLastDeltaX = -10;
2293 87 : psWrkStruct->dfLastDeltaY = -10;
2294 :
2295 : // Alloc space for saving a row of pixels
2296 173 : if( poWK->pafUnifiedSrcDensity == NULL &&
2297 : poWK->panUnifiedSrcValid == NULL &&
2298 : poWK->papanBandSrcValid == NULL )
2299 : {
2300 86 : psWrkStruct->padfRowDensity = NULL;
2301 : }
2302 : else
2303 : {
2304 1 : psWrkStruct->padfRowDensity = (double *)CPLCalloc( nXDist, sizeof(double) );
2305 : }
2306 87 : psWrkStruct->padfRowReal = (double *)CPLCalloc( nXDist, sizeof(double) );
2307 87 : psWrkStruct->padfRowImag = (double *)CPLCalloc( nXDist, sizeof(double) );
2308 :
2309 137 : if( poWK->eResample == GRA_Lanczos &&
2310 : poWK->dfXFilter == 3.0 &&
2311 : poWK->dfYFilter == 3.0 )
2312 : {
2313 50 : psWrkStruct->pfnGWKResample = GWKResampleOptimizedLanczos;
2314 :
2315 50 : const double dfXScale = poWK->dfXScale;
2316 50 : if( dfXScale < 1.0 )
2317 : {
2318 0 : int iMin = poWK->nFiltInitX, iMax = poWK->nXRadius;
2319 0 : while( iMin * dfXScale < -3.0 )
2320 0 : iMin ++;
2321 0 : while( iMax * dfXScale > 3.0 )
2322 0 : iMax --;
2323 :
2324 0 : for(int i = iMin; i <= iMax; ++i)
2325 : {
2326 0 : psWrkStruct->padfWeightsX[i-poWK->nFiltInitX] =
2327 0 : GWKLanczosSinc(i * dfXScale, poWK->dfXFilter) * dfXScale;
2328 : }
2329 : }
2330 :
2331 50 : const double dfYScale = poWK->dfYScale;
2332 50 : if( dfYScale < 1.0 )
2333 : {
2334 0 : int jMin = poWK->nFiltInitY, jMax = poWK->nYRadius;
2335 0 : while( jMin * dfYScale < -3.0 )
2336 0 : jMin ++;
2337 0 : while( jMax * dfYScale > 3.0 )
2338 0 : jMax --;
2339 :
2340 0 : for(int j = jMin; j <= jMax; ++j)
2341 : {
2342 0 : psWrkStruct->padfWeightsY[j-poWK->nFiltInitY] =
2343 0 : GWKLanczosSinc(j * dfYScale, poWK->dfYFilter) * dfYScale;
2344 : }
2345 : }
2346 : }
2347 : else
2348 37 : psWrkStruct->pfnGWKResample = GWKResample;
2349 :
2350 87 : return psWrkStruct;
2351 : }
2352 :
2353 : /************************************************************************/
2354 : /* GWKResampleDeleteWrkStruct() */
2355 : /************************************************************************/
2356 :
2357 87 : static void GWKResampleDeleteWrkStruct(GWKResampleWrkStruct* psWrkStruct)
2358 : {
2359 87 : CPLFree( psWrkStruct->padfWeightsX );
2360 87 : CPLFree( psWrkStruct->padfWeightsY );
2361 87 : CPLFree( psWrkStruct->panCalcX );
2362 87 : CPLFree( psWrkStruct->padfRowDensity );
2363 87 : CPLFree( psWrkStruct->padfRowReal );
2364 87 : CPLFree( psWrkStruct->padfRowImag );
2365 87 : CPLFree( psWrkStruct );
2366 87 : }
2367 :
2368 : /************************************************************************/
2369 : /* GWKResample() */
2370 : /************************************************************************/
2371 :
2372 558 : static int GWKResample( GDALWarpKernel *poWK, int iBand,
2373 : double dfSrcX, double dfSrcY,
2374 : double *pdfDensity,
2375 : double *pdfReal, double *pdfImag,
2376 : GWKResampleWrkStruct* psWrkStruct )
2377 :
2378 : {
2379 : // Save as local variables to avoid following pointers in loops
2380 558 : const int nSrcXSize = poWK->nSrcXSize;
2381 558 : const int nSrcYSize = poWK->nSrcYSize;
2382 :
2383 558 : double dfAccumulatorReal = 0.0, dfAccumulatorImag = 0.0;
2384 558 : double dfAccumulatorDensity = 0.0;
2385 558 : double dfAccumulatorWeight = 0.0;
2386 558 : const int iSrcX = (int) floor( dfSrcX - 0.5 );
2387 558 : const int iSrcY = (int) floor( dfSrcY - 0.5 );
2388 558 : const int iSrcOffset = iSrcX + iSrcY * nSrcXSize;
2389 558 : const double dfDeltaX = dfSrcX - 0.5 - iSrcX;
2390 558 : const double dfDeltaY = dfSrcY - 0.5 - iSrcY;
2391 558 : const int eResample = poWK->eResample;
2392 :
2393 558 : const double dfXScale = poWK->dfXScale, dfYScale = poWK->dfYScale;
2394 558 : const double dfXFilter = poWK->dfXFilter, dfYFilter = poWK->dfYFilter;
2395 :
2396 : int i, j;
2397 558 : const int nXDist = ( poWK->nXRadius + 1 ) * 2;
2398 :
2399 : // Space for saved X weights
2400 558 : double *padfWeightsX = psWrkStruct->padfWeightsX;
2401 558 : char *panCalcX = psWrkStruct->panCalcX;
2402 :
2403 : // Space for saving a row of pixels
2404 558 : double *padfRowDensity = psWrkStruct->padfRowDensity;
2405 558 : double *padfRowReal = psWrkStruct->padfRowReal;
2406 558 : double *padfRowImag = psWrkStruct->padfRowImag;
2407 :
2408 : // Mark as needing calculation (don't calculate the weights yet,
2409 : // because a mask may render it unnecessary)
2410 558 : memset( panCalcX, FALSE, nXDist * sizeof(char) );
2411 :
2412 558 : CPLAssert( eResample == GRA_CubicSpline || eResample == GRA_Lanczos );
2413 :
2414 : // Skip sampling over edge of image
2415 558 : j = poWK->nFiltInitY;
2416 558 : int jMax= poWK->nYRadius;
2417 558 : if( iSrcY + j < 0 )
2418 108 : j = -iSrcY;
2419 558 : if( iSrcY + jMax >= nSrcYSize )
2420 216 : jMax = nSrcYSize - iSrcY - 1;
2421 :
2422 558 : int iMin = poWK->nFiltInitX, iMax = poWK->nXRadius;
2423 558 : if( iSrcX + iMin < 0 )
2424 108 : iMin = -iSrcX;
2425 558 : if( iSrcX + iMax >= nSrcXSize )
2426 216 : iMax = nSrcXSize - iSrcX - 1;
2427 :
2428 558 : const int bXScaleBelow1 = ( dfXScale < 1.0 );
2429 558 : const int bYScaleBelow1 = ( dfYScale < 1.0 );
2430 :
2431 558 : int iRowOffset = iSrcOffset + (j - 1) * nSrcXSize + iMin;
2432 :
2433 : // Loop over pixel rows in the kernel
2434 2358 : for ( ; j <= jMax; ++j )
2435 : {
2436 : double dfWeight1;
2437 :
2438 1800 : iRowOffset += nSrcXSize;
2439 :
2440 : // Get pixel values
2441 : // We can potentially read extra elements after the "normal" end of the source arrays,
2442 : // but the contract of papabySrcImage[iBand], papanBandSrcValid[iBand],
2443 : // panUnifiedSrcValid and pafUnifiedSrcDensity is to have WARP_EXTRA_ELTS
2444 : // reserved at their end.
2445 1800 : if ( !GWKGetPixelRow( poWK, iBand, iRowOffset, (iMax-iMin+2)/2,
2446 : padfRowDensity, padfRowReal, padfRowImag ) )
2447 0 : continue;
2448 :
2449 : // Select the resampling algorithm
2450 1800 : if ( eResample == GRA_CubicSpline )
2451 : // Calculate the Y weight
2452 : dfWeight1 = ( bYScaleBelow1 ) ?
2453 : GWKBSpline(((double)j) * dfYScale) * dfYScale :
2454 1800 : GWKBSpline(((double)j) - dfDeltaY);
2455 : else /*if ( eResample == GRA_Lanczos )*/
2456 : {
2457 0 : if( bYScaleBelow1 )
2458 0 : dfWeight1 = GWKLanczosSinc(j * dfYScale, dfYFilter) * dfYScale;
2459 : else
2460 0 : dfWeight1 = GWKLanczosSinc(j - dfDeltaY, dfYFilter);
2461 : }
2462 :
2463 : // Iterate over pixels in row
2464 7704 : for (i = iMin; i <= iMax; ++i )
2465 : {
2466 : double dfWeight2;
2467 :
2468 : // Skip sampling if pixel has zero density
2469 5904 : if ( padfRowDensity != NULL &&
2470 0 : padfRowDensity[i-iMin] < 0.000000001 )
2471 0 : continue;
2472 :
2473 : // Make or use a cached set of weights for this row
2474 5904 : if ( panCalcX[i-iMin] )
2475 : // Use saved weight value instead of recomputing it
2476 4104 : dfWeight2 = dfWeight1 * padfWeightsX[i-iMin];
2477 : else
2478 : {
2479 : // Choose among possible algorithms
2480 1800 : if ( eResample == GRA_CubicSpline )
2481 : // Calculate & save the X weight
2482 1800 : padfWeightsX[i-iMin] = dfWeight2 = ( bXScaleBelow1 ) ?
2483 : GWKBSpline((double)i * dfXScale) * dfXScale :
2484 1800 : GWKBSpline(dfDeltaX - (double)i);
2485 : else /*if ( eResample == GRA_Lanczos )*/
2486 : {
2487 : // Calculate & save the X weight
2488 0 : if( bXScaleBelow1 )
2489 0 : padfWeightsX[i-iMin] = dfWeight2 =
2490 0 : GWKLanczosSinc(i * dfXScale, dfXFilter) * dfXScale;
2491 : else
2492 0 : padfWeightsX[i-iMin] = dfWeight2 =
2493 0 : GWKLanczosSinc(i - dfDeltaX, dfXFilter);
2494 : }
2495 :
2496 1800 : dfWeight2 *= dfWeight1;
2497 1800 : panCalcX[i-iMin] = TRUE;
2498 : }
2499 :
2500 : // Accumulate!
2501 5904 : dfAccumulatorReal += padfRowReal[i-iMin] * dfWeight2;
2502 5904 : dfAccumulatorImag += padfRowImag[i-iMin] * dfWeight2;
2503 5904 : if( padfRowDensity != NULL )
2504 0 : dfAccumulatorDensity += padfRowDensity[i-iMin] * dfWeight2;
2505 5904 : dfAccumulatorWeight += dfWeight2;
2506 : }
2507 : }
2508 :
2509 558 : if ( dfAccumulatorWeight < 0.000001 ||
2510 : (padfRowDensity != NULL && dfAccumulatorDensity < 0.000001) )
2511 : {
2512 0 : *pdfDensity = 0.0;
2513 0 : return FALSE;
2514 : }
2515 :
2516 : // Calculate the output taking into account weighting
2517 882 : if ( dfAccumulatorWeight < 0.99999 || dfAccumulatorWeight > 1.00001 )
2518 : {
2519 324 : *pdfReal = dfAccumulatorReal / dfAccumulatorWeight;
2520 324 : *pdfImag = dfAccumulatorImag / dfAccumulatorWeight;
2521 324 : if( padfRowDensity != NULL )
2522 0 : *pdfDensity = dfAccumulatorDensity / dfAccumulatorWeight;
2523 : else
2524 324 : *pdfDensity = 1.0;
2525 : }
2526 : else
2527 : {
2528 234 : *pdfReal = dfAccumulatorReal;
2529 234 : *pdfImag = dfAccumulatorImag;
2530 234 : if( padfRowDensity != NULL )
2531 0 : *pdfDensity = dfAccumulatorDensity;
2532 : else
2533 234 : *pdfDensity = 1.0;
2534 : }
2535 :
2536 558 : return TRUE;
2537 : }
2538 :
2539 : /************************************************************************/
2540 : /* GWKResampleOptimizedLanczos() */
2541 : /************************************************************************/
2542 :
2543 278826 : static int GWKResampleOptimizedLanczos( GDALWarpKernel *poWK, int iBand,
2544 : double dfSrcX, double dfSrcY,
2545 : double *pdfDensity,
2546 : double *pdfReal, double *pdfImag,
2547 : GWKResampleWrkStruct* psWrkStruct )
2548 :
2549 : {
2550 : // Save as local variables to avoid following pointers in loops
2551 278826 : const int nSrcXSize = poWK->nSrcXSize;
2552 278826 : const int nSrcYSize = poWK->nSrcYSize;
2553 :
2554 278826 : double dfAccumulatorReal = 0.0, dfAccumulatorImag = 0.0;
2555 278826 : double dfAccumulatorDensity = 0.0;
2556 278826 : double dfAccumulatorWeight = 0.0;
2557 278826 : const int iSrcX = (int) floor( dfSrcX - 0.5 );
2558 278826 : const int iSrcY = (int) floor( dfSrcY - 0.5 );
2559 278826 : const int iSrcOffset = iSrcX + iSrcY * nSrcXSize;
2560 278826 : const double dfDeltaX = dfSrcX - 0.5 - iSrcX;
2561 278826 : const double dfDeltaY = dfSrcY - 0.5 - iSrcY;
2562 :
2563 278826 : const double dfXScale = poWK->dfXScale, dfYScale = poWK->dfYScale;
2564 :
2565 : // Space for saved X weights
2566 278826 : double *padfWeightsX = psWrkStruct->padfWeightsX;
2567 278826 : double *padfWeightsY = psWrkStruct->padfWeightsY;
2568 :
2569 : // Space for saving a row of pixels
2570 278826 : double *padfRowDensity = psWrkStruct->padfRowDensity;
2571 278826 : double *padfRowReal = psWrkStruct->padfRowReal;
2572 278826 : double *padfRowImag = psWrkStruct->padfRowImag;
2573 :
2574 : // Skip sampling over edge of image
2575 278826 : int jMin = poWK->nFiltInitY, jMax= poWK->nYRadius;
2576 278826 : if( iSrcY + jMin < 0 )
2577 12398 : jMin = -iSrcY;
2578 278826 : if( iSrcY + jMax >= nSrcYSize )
2579 934 : jMax = nSrcYSize - iSrcY - 1;
2580 :
2581 278826 : int iMin = poWK->nFiltInitX, iMax = poWK->nXRadius;
2582 278826 : if( iSrcX + iMin < 0 )
2583 11886 : iMin = -iSrcX;
2584 278826 : if( iSrcX + iMax >= nSrcXSize )
2585 934 : iMax = nSrcXSize - iSrcX - 1;
2586 :
2587 278826 : if( dfXScale < 1.0 )
2588 : {
2589 0 : while( iMin * dfXScale < -3.0 )
2590 0 : iMin ++;
2591 0 : while( iMax * dfXScale > 3.0 )
2592 0 : iMax --;
2593 : // padfWeightsX computed in GWKResampleCreateWrkStruct
2594 : }
2595 : else
2596 : {
2597 824284 : while( iMin - dfDeltaX < -3.0 )
2598 266632 : iMin ++;
2599 557652 : while( iMax - dfDeltaX > 3.0 )
2600 0 : iMax --;
2601 :
2602 278826 : if( iSrcX != psWrkStruct->iLastSrcX ||
2603 : dfDeltaX != psWrkStruct->dfLastDeltaX )
2604 : {
2605 : // Optimisation of GWKLanczosSinc(i - dfDeltaX) based on the following
2606 : // trigonometric formulas.
2607 :
2608 : //sin(GWK_PI * (dfBase + k)) = sin(GWK_PI * dfBase) * cos(GWK_PI * k) + cos(GWK_PI * dfBase) * sin(GWK_PI * k)
2609 : //sin(GWK_PI * (dfBase + k)) = dfSinPIBase * cos(GWK_PI * k) + dfCosPIBase * sin(GWK_PI * k)
2610 : //sin(GWK_PI * (dfBase + k)) = dfSinPIBase * cos(GWK_PI * k)
2611 : //sin(GWK_PI * (dfBase + k)) = dfSinPIBase * (((k % 2) == 0) ? 1 : -1)
2612 :
2613 : //sin(GWK_PI / dfR * (dfBase + k)) = sin(GWK_PI / dfR * dfBase) * cos(GWK_PI / dfR * k) + cos(GWK_PI / dfR * dfBase) * sin(GWK_PI / dfR * k)
2614 : //sin(GWK_PI / dfR * (dfBase + k)) = dfSinPIBaseOverR * cos(GWK_PI / dfR * k) + dfCosPIBaseOverR * sin(GWK_PI / dfR * k)
2615 :
2616 278826 : double dfSinPIDeltaXOver3 = sin((-GWK_PI / 3) * dfDeltaX);
2617 278826 : double dfSin2PIDeltaXOver3 = dfSinPIDeltaXOver3 * dfSinPIDeltaXOver3;
2618 : /* ok to use sqrt(1-sin^2) since GWK_PI / 3 * dfDeltaX < PI/2 */
2619 278826 : double dfCosPIDeltaXOver3 = sqrt(1 - dfSin2PIDeltaXOver3);
2620 278826 : double dfSinPIDeltaX = (3-4*dfSin2PIDeltaXOver3)*dfSinPIDeltaXOver3;
2621 278826 : const double dfInvPI2Over3 = 3.0 / (GWK_PI * GWK_PI);
2622 278826 : double dfInvPI2Over3xSinPIDeltaX = dfInvPI2Over3 * dfSinPIDeltaX;
2623 : double dfInvPI2Over3xSinPIDeltaXxm0d5SinPIDeltaXOver3 =
2624 278826 : -0.5 * dfInvPI2Over3xSinPIDeltaX * dfSinPIDeltaXOver3;
2625 278826 : const double dfSinPIOver3 = 0.8660254037844386;
2626 : double dfInvPI2Over3xSinPIDeltaXxSinPIOver3xCosPIDeltaXOver3 =
2627 278826 : dfSinPIOver3 * dfInvPI2Over3xSinPIDeltaX * dfCosPIDeltaXOver3;
2628 : double padfCst[] = {
2629 : dfInvPI2Over3xSinPIDeltaX * dfSinPIDeltaXOver3,
2630 : dfInvPI2Over3xSinPIDeltaXxm0d5SinPIDeltaXOver3 -
2631 : dfInvPI2Over3xSinPIDeltaXxSinPIOver3xCosPIDeltaXOver3,
2632 : dfInvPI2Over3xSinPIDeltaXxm0d5SinPIDeltaXOver3 +
2633 278826 : dfInvPI2Over3xSinPIDeltaXxSinPIOver3xCosPIDeltaXOver3 };
2634 :
2635 1934940 : for (int i = iMin; i <= iMax; ++i )
2636 : {
2637 1656114 : const double dfX = i - dfDeltaX;
2638 1656114 : if (dfX == 0.0)
2639 682 : padfWeightsX[i-poWK->nFiltInitX] = 1.0;
2640 : else
2641 1655432 : padfWeightsX[i-poWK->nFiltInitX] =
2642 1655432 : padfCst[(i + 3) % 3] / (dfX * dfX);
2643 : //CPLAssert(fabs(padfWeightsX[i-poWK->nFiltInitX] - GWKLanczosSinc(dfX, 3.0)) < 1e-10);
2644 : }
2645 :
2646 278826 : psWrkStruct->iLastSrcX = iSrcX;
2647 278826 : psWrkStruct->dfLastDeltaX = dfDeltaX;
2648 : }
2649 : }
2650 :
2651 278826 : if( dfYScale < 1.0 )
2652 : {
2653 0 : while( jMin * dfYScale < -3.0 )
2654 0 : jMin ++;
2655 0 : while( jMax * dfYScale > 3.0 )
2656 0 : jMax --;
2657 : // padfWeightsY computed in GWKResampleCreateWrkStruct
2658 : }
2659 : else
2660 : {
2661 823772 : while( jMin - dfDeltaY < -3.0 )
2662 266120 : jMin ++;
2663 557652 : while( jMax - dfDeltaY > 3.0 )
2664 0 : jMax --;
2665 :
2666 278826 : if( iSrcY != psWrkStruct->iLastSrcY ||
2667 : dfDeltaY != psWrkStruct->dfLastDeltaY )
2668 : {
2669 2364 : double dfSinPIDeltaYOver3 = sin((-GWK_PI / 3) * dfDeltaY);
2670 2364 : double dfSin2PIDeltaYOver3 = dfSinPIDeltaYOver3 * dfSinPIDeltaYOver3;
2671 : /* ok to use sqrt(1-sin^2) since GWK_PI / 3 * dfDeltaY < PI/2 */
2672 2364 : double dfCosPIDeltaYOver3 = sqrt(1 - dfSin2PIDeltaYOver3);
2673 2364 : double dfSinPIDeltaY = (3-4*dfSin2PIDeltaYOver3)*dfSinPIDeltaYOver3;
2674 2364 : const double dfInvPI2Over3 = 3.0 / (GWK_PI * GWK_PI);
2675 2364 : double dfInvPI2Over3xSinPIDeltaY = dfInvPI2Over3 * dfSinPIDeltaY;
2676 : double dfInvPI2Over3xSinPIDeltaYxm0d5SinPIDeltaYOver3 =
2677 2364 : -0.5 * dfInvPI2Over3xSinPIDeltaY * dfSinPIDeltaYOver3;
2678 2364 : const double dfSinPIOver3 = 0.8660254037844386;
2679 : double dfInvPI2Over3xSinPIDeltaYxSinPIOver3xCosPIDeltaYOver3 =
2680 2364 : dfSinPIOver3 * dfInvPI2Over3xSinPIDeltaY * dfCosPIDeltaYOver3;
2681 : double padfCst[] = {
2682 : dfInvPI2Over3xSinPIDeltaY * dfSinPIDeltaYOver3,
2683 : dfInvPI2Over3xSinPIDeltaYxm0d5SinPIDeltaYOver3 -
2684 : dfInvPI2Over3xSinPIDeltaYxSinPIOver3xCosPIDeltaYOver3,
2685 : dfInvPI2Over3xSinPIDeltaYxm0d5SinPIDeltaYOver3 +
2686 2364 : dfInvPI2Over3xSinPIDeltaYxSinPIOver3xCosPIDeltaYOver3 };
2687 :
2688 15547 : for ( int j = jMin; j <= jMax; ++j )
2689 : {
2690 13183 : const double dfY = j - dfDeltaY;
2691 13183 : if (dfY == 0.0)
2692 132 : padfWeightsY[j-poWK->nFiltInitY] = 1.0;
2693 : else
2694 13051 : padfWeightsY[j-poWK->nFiltInitY] =
2695 13051 : padfCst[(j + 3) % 3] / (dfY * dfY);
2696 : //CPLAssert(fabs(padfWeightsY[j-poWK->nFiltInitY] - GWKLanczosSinc(dfY, 3.0)) < 1e-10);
2697 : }
2698 :
2699 2364 : psWrkStruct->iLastSrcY = iSrcY;
2700 2364 : psWrkStruct->dfLastDeltaY = dfDeltaY;
2701 : }
2702 : }
2703 :
2704 278826 : int iRowOffset = iSrcOffset + (jMin - 1) * nSrcXSize + iMin;
2705 :
2706 : // If we have no density information, we can simply compute the
2707 : // accumulated weight.
2708 278826 : if( padfRowDensity == NULL )
2709 : {
2710 278826 : double dfRowAccWeight = 0.0;
2711 1934940 : for (int i = iMin; i <= iMax; ++i )
2712 : {
2713 1656114 : dfRowAccWeight += padfWeightsX[i-poWK->nFiltInitX];
2714 : }
2715 278826 : double dfColAccWeight = 0.0;
2716 1934940 : for ( int j = jMin; j <= jMax; ++j )
2717 : {
2718 1656114 : dfColAccWeight += padfWeightsY[j-poWK->nFiltInitY];
2719 : }
2720 278826 : dfAccumulatorWeight = dfRowAccWeight * dfColAccWeight;
2721 :
2722 278826 : if( !GDALDataTypeIsComplex(poWK->eWorkingDataType) )
2723 278578 : padfRowImag = NULL;
2724 : }
2725 :
2726 : // Loop over pixel rows in the kernel
2727 1934940 : for ( int j = jMin; j <= jMax; ++j )
2728 : {
2729 : double dfWeight1;
2730 :
2731 1656114 : iRowOffset += nSrcXSize;
2732 :
2733 : // Get pixel values
2734 : // We can potentially read extra elements after the "normal" end of the source arrays,
2735 : // but the contract of papabySrcImage[iBand], papanBandSrcValid[iBand],
2736 : // panUnifiedSrcValid and pafUnifiedSrcDensity is to have WARP_EXTRA_ELTS
2737 : // reserved at their end.
2738 1656114 : if ( !GWKGetPixelRow( poWK, iBand, iRowOffset, (iMax-iMin+2)/2,
2739 : padfRowDensity, padfRowReal, padfRowImag ) )
2740 0 : continue;
2741 :
2742 1656114 : dfWeight1 = padfWeightsY[j-poWK->nFiltInitY];
2743 :
2744 : // Iterate over pixels in row
2745 1656114 : if ( padfRowDensity != NULL )
2746 : {
2747 0 : for (int i = iMin; i <= iMax; ++i )
2748 : {
2749 : double dfWeight2;
2750 :
2751 : // Skip sampling if pixel has zero density
2752 0 : if ( padfRowDensity[i - iMin] < 0.000000001 )
2753 0 : continue;
2754 :
2755 : // Use a cached set of weights for this row
2756 0 : dfWeight2 = dfWeight1 * padfWeightsX[i- poWK->nFiltInitX];
2757 :
2758 : // Accumulate!
2759 0 : dfAccumulatorReal += padfRowReal[i - iMin] * dfWeight2;
2760 0 : dfAccumulatorImag += padfRowImag[i - iMin] * dfWeight2;
2761 0 : dfAccumulatorDensity += padfRowDensity[i - iMin] * dfWeight2;
2762 0 : dfAccumulatorWeight += dfWeight2;
2763 : }
2764 : }
2765 1656114 : else if( padfRowImag == NULL )
2766 : {
2767 1654938 : double dfRowAccReal = 0.0;
2768 11487773 : for (int i = iMin; i <= iMax; ++i )
2769 : {
2770 9832835 : double dfWeight2 = padfWeightsX[i- poWK->nFiltInitX];
2771 :
2772 : // Accumulate!
2773 9832835 : dfRowAccReal += padfRowReal[i - iMin] * dfWeight2;
2774 : }
2775 :
2776 1654938 : dfAccumulatorReal += dfRowAccReal * dfWeight1;
2777 : }
2778 : else
2779 : {
2780 1176 : double dfRowAccReal = 0.0;
2781 1176 : double dfRowAccImag = 0.0;
2782 7040 : for (int i = iMin; i <= iMax; ++i )
2783 : {
2784 5864 : double dfWeight2 = padfWeightsX[i- poWK->nFiltInitX];
2785 :
2786 : // Accumulate!
2787 5864 : dfRowAccReal += padfRowReal[i - iMin] * dfWeight2;
2788 5864 : dfRowAccImag += padfRowImag[i - iMin] * dfWeight2;
2789 : }
2790 :
2791 1176 : dfAccumulatorReal += dfRowAccReal * dfWeight1;
2792 1176 : dfAccumulatorImag += dfRowAccImag * dfWeight1;
2793 : }
2794 : }
2795 :
2796 278826 : if ( dfAccumulatorWeight < 0.000001 ||
2797 : (padfRowDensity != NULL && dfAccumulatorDensity < 0.000001) )
2798 : {
2799 0 : *pdfDensity = 0.0;
2800 0 : return FALSE;
2801 : }
2802 :
2803 : // Calculate the output taking into account weighting
2804 556970 : if ( dfAccumulatorWeight < 0.99999 || dfAccumulatorWeight > 1.00001 )
2805 : {
2806 278144 : const double dfInvAcc = 1.0 / dfAccumulatorWeight;
2807 278144 : *pdfReal = dfAccumulatorReal * dfInvAcc;
2808 278144 : *pdfImag = dfAccumulatorImag * dfInvAcc;
2809 278144 : if( padfRowDensity != NULL )
2810 0 : *pdfDensity = dfAccumulatorDensity * dfInvAcc;
2811 : else
2812 278144 : *pdfDensity = 1.0;
2813 : }
2814 : else
2815 : {
2816 682 : *pdfReal = dfAccumulatorReal;
2817 682 : *pdfImag = dfAccumulatorImag;
2818 682 : if( padfRowDensity != NULL )
2819 0 : *pdfDensity = dfAccumulatorDensity;
2820 : else
2821 682 : *pdfDensity = 1.0;
2822 : }
2823 :
2824 278826 : return TRUE;
2825 : }
2826 :
2827 278207 : static int GWKCubicSplineResampleNoMasksByte( GDALWarpKernel *poWK, int iBand,
2828 : double dfSrcX, double dfSrcY,
2829 : GByte *pbValue, double *padfBSpline )
2830 :
2831 : {
2832 : // Commonly used; save locally
2833 278207 : int nSrcXSize = poWK->nSrcXSize;
2834 278207 : int nSrcYSize = poWK->nSrcYSize;
2835 :
2836 278207 : double dfAccumulator = 0.0;
2837 278207 : int iSrcX = (int) floor( dfSrcX - 0.5 );
2838 278207 : int iSrcY = (int) floor( dfSrcY - 0.5 );
2839 278207 : int iSrcOffset = iSrcX + iSrcY * nSrcXSize;
2840 278207 : double dfDeltaX = dfSrcX - 0.5 - iSrcX;
2841 278207 : double dfDeltaY = dfSrcY - 0.5 - iSrcY;
2842 :
2843 278207 : double dfXScale = poWK->dfXScale;
2844 278207 : double dfYScale = poWK->dfYScale;
2845 278207 : int nXRadius = poWK->nXRadius;
2846 278207 : int nYRadius = poWK->nYRadius;
2847 :
2848 278207 : GByte* pabySrcBand = poWK->papabySrcImage[iBand];
2849 :
2850 : // Politely refusing to process invalid coordinates or obscenely small image
2851 278207 : if ( iSrcX >= nSrcXSize || iSrcY >= nSrcYSize
2852 : || nXRadius > nSrcXSize || nYRadius > nSrcYSize )
2853 1 : return GWKBilinearResampleNoMasksByte( poWK, iBand, dfSrcX, dfSrcY, pbValue);
2854 :
2855 : // Loop over all rows in the kernel
2856 : int j, jC;
2857 1391030 : for ( jC = 0, j = 1 - nYRadius; j <= nYRadius; ++j, ++jC )
2858 : {
2859 : int iSampJ;
2860 : // Calculate the Y weight
2861 : double dfWeight1 = ( dfYScale < 1.0 ) ?
2862 : GWKBSpline((double)j * dfYScale) * dfYScale :
2863 1112824 : GWKBSpline((double)j - dfDeltaY);
2864 :
2865 : // Flip sampling over edge of image
2866 1112824 : if ( iSrcY + j < 0 )
2867 6676 : iSampJ = iSrcOffset - (iSrcY + j) * nSrcXSize;
2868 1106148 : else if ( iSrcY + j >= nSrcYSize )
2869 556 : iSampJ = iSrcOffset + (2*nSrcYSize - 2*iSrcY - j - 1) * nSrcXSize;
2870 : else
2871 1105592 : iSampJ = iSrcOffset + j * nSrcXSize;
2872 :
2873 : // Loop over all pixels in the row
2874 : int i, iC;
2875 5564120 : for ( iC = 0, i = 1 - nXRadius; i <= nXRadius; ++i, ++iC )
2876 : {
2877 : int iSampI;
2878 : double dfWeight2;
2879 :
2880 : // Flip sampling over edge of image
2881 4451296 : if ( iSrcX + i < 0 )
2882 26704 : iSampI = -iSrcX - i;
2883 4424592 : else if ( iSrcX + i >= nSrcXSize )
2884 2224 : iSampI = 2*nSrcXSize - 2*iSrcX - i - 1;
2885 : else
2886 4422368 : iSampI = i;
2887 :
2888 : // Make a cached set of GWKBSpline values
2889 4451296 : if( jC == 0 )
2890 : {
2891 : // Calculate & save the X weight
2892 1112824 : dfWeight2 = padfBSpline[iC] = ((dfXScale < 1.0 ) ?
2893 : GWKBSpline((double)i * dfXScale) * dfXScale :
2894 1112824 : GWKBSpline(dfDeltaX - (double)i));
2895 1112824 : dfWeight2 *= dfWeight1;
2896 : }
2897 : else
2898 3338472 : dfWeight2 = dfWeight1 * padfBSpline[iC];
2899 :
2900 : // Retrieve the pixel & accumulate
2901 4451296 : dfAccumulator += (double)pabySrcBand[iSampI+iSampJ] * dfWeight2;
2902 : }
2903 : }
2904 :
2905 278206 : if ( dfAccumulator < 0.0 )
2906 0 : *pbValue = 0;
2907 278206 : else if ( dfAccumulator > 255.0 )
2908 19 : *pbValue = 255;
2909 : else
2910 278187 : *pbValue = (GByte)(0.5 + dfAccumulator);
2911 :
2912 278206 : return TRUE;
2913 : }
2914 :
2915 262207 : static int GWKCubicSplineResampleNoMasksShort( GDALWarpKernel *poWK, int iBand,
2916 : double dfSrcX, double dfSrcY,
2917 : GInt16 *piValue, double *padfBSpline )
2918 :
2919 : {
2920 : //Save src size to local var
2921 262207 : int nSrcXSize = poWK->nSrcXSize;
2922 262207 : int nSrcYSize = poWK->nSrcYSize;
2923 :
2924 262207 : double dfAccumulator = 0.0;
2925 262207 : int iSrcX = (int) floor( dfSrcX - 0.5 );
2926 262207 : int iSrcY = (int) floor( dfSrcY - 0.5 );
2927 262207 : int iSrcOffset = iSrcX + iSrcY * nSrcXSize;
2928 262207 : double dfDeltaX = dfSrcX - 0.5 - iSrcX;
2929 262207 : double dfDeltaY = dfSrcY - 0.5 - iSrcY;
2930 :
2931 262207 : double dfXScale = poWK->dfXScale;
2932 262207 : double dfYScale = poWK->dfYScale;
2933 262207 : int nXRadius = poWK->nXRadius;
2934 262207 : int nYRadius = poWK->nYRadius;
2935 :
2936 : // Save band array pointer to local var; cast here instead of later
2937 262207 : GInt16* pabySrcBand = ((GInt16 *)poWK->papabySrcImage[iBand]);
2938 :
2939 : // Politely refusing to process invalid coordinates or obscenely small image
2940 262207 : if ( iSrcX >= nSrcXSize || iSrcY >= nSrcYSize
2941 : || nXRadius > nSrcXSize || nYRadius > nSrcYSize )
2942 1 : return GWKBilinearResampleNoMasksShort( poWK, iBand, dfSrcX, dfSrcY, piValue);
2943 :
2944 : // Loop over all pixels in the kernel
2945 : int j, jC;
2946 1311030 : for ( jC = 0, j = 1 - nYRadius; j <= nYRadius; ++j, ++jC )
2947 : {
2948 : int iSampJ;
2949 :
2950 : // Calculate the Y weight
2951 : double dfWeight1 = ( dfYScale < 1.0 ) ?
2952 : GWKBSpline((double)j * dfYScale) * dfYScale :
2953 1048824 : GWKBSpline((double)j - dfDeltaY);
2954 :
2955 : // Flip sampling over edge of image
2956 1048824 : if ( iSrcY + j < 0 )
2957 6156 : iSampJ = iSrcOffset - (iSrcY + j) * nSrcXSize;
2958 1042668 : else if ( iSrcY + j >= nSrcYSize )
2959 36 : iSampJ = iSrcOffset + (2*nSrcYSize - 2*iSrcY - j - 1) * nSrcXSize;
2960 : else
2961 1042632 : iSampJ = iSrcOffset + j * nSrcXSize;
2962 :
2963 : // Loop over all pixels in row
2964 : int i, iC;
2965 5244120 : for ( iC = 0, i = 1 - nXRadius; i <= nXRadius; ++i, ++iC )
2966 : {
2967 : int iSampI;
2968 : double dfWeight2;
2969 :
2970 : // Flip sampling over edge of image
2971 4195296 : if ( iSrcX + i < 0 )
2972 24624 : iSampI = -iSrcX - i;
2973 4170672 : else if(iSrcX + i >= nSrcXSize)
2974 144 : iSampI = 2*nSrcXSize - 2*iSrcX - i - 1;
2975 : else
2976 4170528 : iSampI = i;
2977 :
2978 : // Make a cached set of GWKBSpline values
2979 4195296 : if ( jC == 0 )
2980 : {
2981 : // Calculate & save the X weight
2982 1048824 : dfWeight2 = padfBSpline[iC] = ((dfXScale < 1.0 ) ?
2983 : GWKBSpline((double)i * dfXScale) * dfXScale :
2984 1048824 : GWKBSpline(dfDeltaX - (double)i));
2985 1048824 : dfWeight2 *= dfWeight1;
2986 : } else
2987 3146472 : dfWeight2 = dfWeight1 * padfBSpline[iC];
2988 :
2989 4195296 : dfAccumulator += (double)pabySrcBand[iSampI + iSampJ] * dfWeight2;
2990 : }
2991 : }
2992 :
2993 262206 : *piValue = (GInt16)(0.5 + dfAccumulator);
2994 :
2995 262206 : return TRUE;
2996 : }
2997 :
2998 : /************************************************************************/
2999 : /* GWKOpenCLCase() */
3000 : /* */
3001 : /* This is identical to GWKGeneralCase(), but functions via */
3002 : /* OpenCL. This means we have vector optimization (SSE) and/or */
3003 : /* GPU optimization depending on our prefs. The code itsef is */
3004 : /* general and not optimized, but by defining constants we can */
3005 : /* make some pretty darn good code on the fly. */
3006 : /************************************************************************/
3007 :
3008 : #if defined(HAVE_OPENCL)
3009 : static CPLErr GWKOpenCLCase( GDALWarpKernel *poWK )
3010 : {
3011 : int iDstY, iBand;
3012 : int nDstXSize = poWK->nDstXSize, nDstYSize = poWK->nDstYSize;
3013 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
3014 : int nDstXOff = poWK->nDstXOff , nDstYOff = poWK->nDstYOff;
3015 : int nSrcXOff = poWK->nSrcXOff , nSrcYOff = poWK->nSrcYOff;
3016 : CPLErr eErr = CE_None;
3017 : struct oclWarper *warper;
3018 : cl_channel_type imageFormat;
3019 : int useImag = FALSE;
3020 : OCLResampAlg resampAlg;
3021 : cl_int err;
3022 :
3023 : switch ( poWK->eWorkingDataType )
3024 : {
3025 : case GDT_Byte:
3026 : imageFormat = CL_UNORM_INT8;
3027 : break;
3028 : case GDT_UInt16:
3029 : imageFormat = CL_UNORM_INT16;
3030 : break;
3031 : case GDT_CInt16:
3032 : useImag = TRUE;
3033 : case GDT_Int16:
3034 : imageFormat = CL_SNORM_INT16;
3035 : break;
3036 : case GDT_CFloat32:
3037 : useImag = TRUE;
3038 : case GDT_Float32:
3039 : imageFormat = CL_FLOAT;
3040 : break;
3041 : default:
3042 : // We don't support higher precision formats
3043 : CPLDebug( "OpenCL",
3044 : "Unsupported resampling OpenCL data type %d.",
3045 : (int) poWK->eWorkingDataType );
3046 : return CE_Warning;
3047 : }
3048 :
3049 : switch (poWK->eResample)
3050 : {
3051 : case GRA_Bilinear:
3052 : resampAlg = OCL_Bilinear;
3053 : break;
3054 : case GRA_Cubic:
3055 : resampAlg = OCL_Cubic;
3056 : break;
3057 : case GRA_CubicSpline:
3058 : resampAlg = OCL_CubicSpline;
3059 : break;
3060 : case GRA_Lanczos:
3061 : resampAlg = OCL_Lanczos;
3062 : break;
3063 : default:
3064 : // We don't support higher precision formats
3065 : CPLDebug( "OpenCL",
3066 : "Unsupported resampling OpenCL resampling alg %d.",
3067 : (int) poWK->eResample );
3068 : return CE_Warning;
3069 : }
3070 :
3071 : // Using a factor of 2 or 4 seems to have much less rounding error than 3 on the GPU.
3072 : // Then the rounding error can cause strange artifacting under the right conditions.
3073 : warper = GDALWarpKernelOpenCL_createEnv(nSrcXSize, nSrcYSize,
3074 : nDstXSize, nDstYSize,
3075 : imageFormat, poWK->nBands, 4,
3076 : useImag, poWK->papanBandSrcValid != NULL,
3077 : poWK->pafDstDensity,
3078 : poWK->padfDstNoDataReal,
3079 : resampAlg, &err );
3080 :
3081 : if(err != CL_SUCCESS || warper == NULL)
3082 : {
3083 : eErr = CE_Warning;
3084 : if (warper != NULL)
3085 : goto free_warper;
3086 : return eErr;
3087 : }
3088 :
3089 : CPLDebug( "GDAL", "GDALWarpKernel()::GWKOpenCLCase()\n"
3090 : "Src=%d,%d,%dx%d Dst=%d,%d,%dx%d",
3091 : nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
3092 : nDstXOff, nDstYOff, nDstXSize, nDstYSize );
3093 :
3094 : if( !poWK->pfnProgress( poWK->dfProgressBase, "", poWK->pProgress ) )
3095 : {
3096 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
3097 : eErr = CE_Failure;
3098 : goto free_warper;
3099 : }
3100 :
3101 : /* ==================================================================== */
3102 : /* Loop over bands. */
3103 : /* ==================================================================== */
3104 : for( iBand = 0; iBand < poWK->nBands; iBand++ ) {
3105 : if( poWK->papanBandSrcValid != NULL && poWK->papanBandSrcValid[iBand] != NULL) {
3106 : GDALWarpKernelOpenCL_setSrcValid(warper, (int *)poWK->papanBandSrcValid[iBand], iBand);
3107 : if(err != CL_SUCCESS)
3108 : {
3109 : CPLError( CE_Failure, CPLE_AppDefined,
3110 : "OpenCL routines reported failure (%d) on line %d.", (int) err, __LINE__ );
3111 : eErr = CE_Failure;
3112 : goto free_warper;
3113 : }
3114 : }
3115 :
3116 : err = GDALWarpKernelOpenCL_setSrcImg(warper, poWK->papabySrcImage[iBand], iBand);
3117 : if(err != CL_SUCCESS)
3118 : {
3119 : CPLError( CE_Failure, CPLE_AppDefined,
3120 : "OpenCL routines reported failure (%d) on line %d.", (int) err, __LINE__ );
3121 : eErr = CE_Failure;
3122 : goto free_warper;
3123 : }
3124 :
3125 : err = GDALWarpKernelOpenCL_setDstImg(warper, poWK->papabyDstImage[iBand], iBand);
3126 : if(err != CL_SUCCESS)
3127 : {
3128 : CPLError( CE_Failure, CPLE_AppDefined,
3129 : "OpenCL routines reported failure (%d) on line %d.", (int) err, __LINE__ );
3130 : eErr = CE_Failure;
3131 : goto free_warper;
3132 : }
3133 : }
3134 :
3135 : /* -------------------------------------------------------------------- */
3136 : /* Allocate x,y,z coordinate arrays for transformation ... one */
3137 : /* scanlines worth of positions. */
3138 : /* -------------------------------------------------------------------- */
3139 : double *padfX, *padfY, *padfZ;
3140 : int *pabSuccess;
3141 :
3142 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3143 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3144 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3145 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
3146 :
3147 : /* ==================================================================== */
3148 : /* Loop over output lines. */
3149 : /* ==================================================================== */
3150 : for( iDstY = 0; iDstY < nDstYSize && eErr == CE_None; ++iDstY )
3151 : {
3152 : int iDstX;
3153 :
3154 : /* ---------------------------------------------------------------- */
3155 : /* Setup points to transform to source image space. */
3156 : /* ---------------------------------------------------------------- */
3157 : for( iDstX = 0; iDstX < nDstXSize; ++iDstX )
3158 : {
3159 : padfX[iDstX] = iDstX + 0.5 + nDstXOff;
3160 : padfY[iDstX] = iDstY + 0.5 + nDstYOff;
3161 : padfZ[iDstX] = 0.0;
3162 : }
3163 :
3164 : /* ---------------------------------------------------------------- */
3165 : /* Transform the points from destination pixel/line coordinates*/
3166 : /* to source pixel/line coordinates. */
3167 : /* ---------------------------------------------------------------- */
3168 : poWK->pfnTransformer( poWK->pTransformerArg, TRUE, nDstXSize,
3169 : padfX, padfY, padfZ, pabSuccess );
3170 :
3171 : err = GDALWarpKernelOpenCL_setCoordRow(warper, padfX, padfY,
3172 : nSrcXOff, nSrcYOff,
3173 : pabSuccess, iDstY);
3174 : if(err != CL_SUCCESS)
3175 : {
3176 : CPLError( CE_Failure, CPLE_AppDefined,
3177 : "OpenCL routines reported failure (%d) on line %d.", (int) err, __LINE__ );
3178 : return CE_Failure;
3179 : }
3180 :
3181 : //Update the valid & density masks because we don't do so in the kernel
3182 : for( iDstX = 0; iDstX < nDstXSize && eErr == CE_None; iDstX++ )
3183 : {
3184 : double dfX = padfX[iDstX];
3185 : double dfY = padfY[iDstX];
3186 : int iDstOffset = iDstX + iDstY * nDstXSize;
3187 :
3188 : //See GWKGeneralCase() for appropriate commenting
3189 : if( !pabSuccess[iDstX] || dfX < nSrcXOff || dfY < nSrcYOff )
3190 : continue;
3191 :
3192 : int iSrcX = ((int) dfX) - nSrcXOff;
3193 : int iSrcY = ((int) dfY) - nSrcYOff;
3194 :
3195 : if( iSrcX < 0 || iSrcX >= nSrcXSize || iSrcY < 0 || iSrcY >= nSrcYSize )
3196 : continue;
3197 :
3198 : int iSrcOffset = iSrcX + iSrcY * nSrcXSize;
3199 : double dfDensity = 1.0;
3200 :
3201 : if( poWK->pafUnifiedSrcDensity != NULL
3202 : && iSrcX >= 0 && iSrcY >= 0
3203 : && iSrcX < nSrcXSize && iSrcY < nSrcYSize )
3204 : dfDensity = poWK->pafUnifiedSrcDensity[iSrcOffset];
3205 :
3206 : GWKOverlayDensity( poWK, iDstOffset, dfDensity );
3207 :
3208 : //Because this is on the bit-wise level, it can't be done well in OpenCL
3209 : if( poWK->panDstValid != NULL )
3210 : poWK->panDstValid[iDstOffset>>5] |= 0x01 << (iDstOffset & 0x1f);
3211 : }
3212 : }
3213 :
3214 : CPLFree( padfX );
3215 : CPLFree( padfY );
3216 : CPLFree( padfZ );
3217 : CPLFree( pabSuccess );
3218 :
3219 : err = GDALWarpKernelOpenCL_runResamp(warper,
3220 : poWK->pafUnifiedSrcDensity,
3221 : poWK->panUnifiedSrcValid,
3222 : poWK->pafDstDensity,
3223 : poWK->panDstValid,
3224 : poWK->dfXScale, poWK->dfYScale,
3225 : poWK->dfXFilter, poWK->dfYFilter,
3226 : poWK->nXRadius, poWK->nYRadius,
3227 : poWK->nFiltInitX, poWK->nFiltInitY);
3228 :
3229 : if(err != CL_SUCCESS)
3230 : {
3231 : CPLError( CE_Failure, CPLE_AppDefined,
3232 : "OpenCL routines reported failure (%d) on line %d.", (int) err, __LINE__ );
3233 : eErr = CE_Failure;
3234 : goto free_warper;
3235 : }
3236 :
3237 : /* ==================================================================== */
3238 : /* Loop over output lines. */
3239 : /* ==================================================================== */
3240 : for( iDstY = 0; iDstY < nDstYSize && eErr == CE_None; iDstY++ )
3241 : {
3242 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
3243 : {
3244 : int iDstX;
3245 : void *rowReal, *rowImag;
3246 : GByte *pabyDst = poWK->papabyDstImage[iBand];
3247 :
3248 : err = GDALWarpKernelOpenCL_getRow(warper, &rowReal, &rowImag, iDstY, iBand);
3249 : if(err != CL_SUCCESS)
3250 : {
3251 : CPLError( CE_Failure, CPLE_AppDefined,
3252 : "OpenCL routines reported failure (%d) on line %d.", (int) err, __LINE__ );
3253 : eErr = CE_Failure;
3254 : goto free_warper;
3255 : }
3256 :
3257 : //Copy the data from the warper to GDAL's memory
3258 : switch ( poWK->eWorkingDataType )
3259 : {
3260 : case GDT_Byte:
3261 : memcpy((void **)&(((GByte *)pabyDst)[iDstY*nDstXSize]),
3262 : rowReal, sizeof(GByte)*nDstXSize);
3263 : break;
3264 : case GDT_Int16:
3265 : memcpy((void **)&(((GInt16 *)pabyDst)[iDstY*nDstXSize]),
3266 : rowReal, sizeof(GInt16)*nDstXSize);
3267 : break;
3268 : case GDT_UInt16:
3269 : memcpy((void **)&(((GUInt16 *)pabyDst)[iDstY*nDstXSize]),
3270 : rowReal, sizeof(GUInt16)*nDstXSize);
3271 : break;
3272 : case GDT_Float32:
3273 : memcpy((void **)&(((float *)pabyDst)[iDstY*nDstXSize]),
3274 : rowReal, sizeof(float)*nDstXSize);
3275 : break;
3276 : case GDT_CInt16:
3277 : {
3278 : GInt16 *pabyDstI16 = &(((GInt16 *)pabyDst)[iDstY*nDstXSize]);
3279 : for (iDstX = 0; iDstX < nDstXSize; iDstX++) {
3280 : pabyDstI16[iDstX*2 ] = ((GInt16 *)rowReal)[iDstX];
3281 : pabyDstI16[iDstX*2+1] = ((GInt16 *)rowImag)[iDstX];
3282 : }
3283 : }
3284 : break;
3285 : case GDT_CFloat32:
3286 : {
3287 : float *pabyDstF32 = &(((float *)pabyDst)[iDstY*nDstXSize]);
3288 : for (iDstX = 0; iDstX < nDstXSize; iDstX++) {
3289 : pabyDstF32[iDstX*2 ] = ((float *)rowReal)[iDstX];
3290 : pabyDstF32[iDstX*2+1] = ((float *)rowImag)[iDstX];
3291 : }
3292 : }
3293 : break;
3294 : default:
3295 : // We don't support higher precision formats
3296 : CPLError( CE_Failure, CPLE_AppDefined,
3297 : "Unsupported resampling OpenCL data type %d.", (int) poWK->eWorkingDataType );
3298 : eErr = CE_Failure;
3299 : goto free_warper;
3300 : }
3301 : }
3302 : }
3303 : free_warper:
3304 : if((err = GDALWarpKernelOpenCL_deleteEnv(warper)) != CL_SUCCESS)
3305 : {
3306 : CPLError( CE_Failure, CPLE_AppDefined,
3307 : "OpenCL routines reported failure (%d) on line %d.", (int) err, __LINE__ );
3308 : return CE_Failure;
3309 : }
3310 :
3311 : return eErr;
3312 : }
3313 : #endif /* defined(HAVE_OPENCL) */
3314 :
3315 :
3316 : #define COMPUTE_iSrcOffset(_pabSuccess, _iDstX, _padfX, _padfY, _poWK, _nSrcXSize, _nSrcYSize) \
3317 : if( !_pabSuccess[_iDstX] ) \
3318 : continue; \
3319 : \
3320 : /* -------------------------------------------------------------------- */ \
3321 : /* Figure out what pixel we want in our source raster, and skip */ \
3322 : /* further processing if it is well off the source image. */ \
3323 : /* -------------------------------------------------------------------- */ \
3324 : /* We test against the value before casting to avoid the */ \
3325 : /* problem of asymmetric truncation effects around zero. That is */ \
3326 : /* -0.5 will be 0 when cast to an int. */ \
3327 : if( _padfX[_iDstX] < _poWK->nSrcXOff \
3328 : || _padfY[_iDstX] < _poWK->nSrcYOff ) \
3329 : continue; \
3330 : \
3331 : int iSrcX, iSrcY, iSrcOffset;\
3332 : \
3333 : iSrcX = ((int) (_padfX[_iDstX] + 1e-10)) - _poWK->nSrcXOff;\
3334 : iSrcY = ((int) (_padfY[_iDstX] + 1e-10)) - _poWK->nSrcYOff;\
3335 : \
3336 : /* If operating outside natural projection area, padfX/Y can be */ \
3337 : /* a very huge positive number, that becomes -2147483648 in the */ \
3338 : /* int trucation. So it is necessary to test now for non negativeness. */ \
3339 : if( iSrcX < 0 || iSrcX >= _nSrcXSize || iSrcY < 0 || iSrcY >= _nSrcYSize )\
3340 : continue;\
3341 : \
3342 : iSrcOffset = iSrcX + iSrcY * _nSrcXSize;
3343 :
3344 : /************************************************************************/
3345 : /* GWKGeneralCase() */
3346 : /* */
3347 : /* This is the most general case. It attempts to handle all */
3348 : /* possible features with relatively little concern for */
3349 : /* efficiency. */
3350 : /************************************************************************/
3351 :
3352 : static void GWKGeneralCaseThread(void* pData);
3353 :
3354 193 : static CPLErr GWKGeneralCase( GDALWarpKernel *poWK )
3355 : {
3356 193 : return GWKRun( poWK, "GWKGeneralCase", GWKGeneralCaseThread );
3357 : }
3358 :
3359 193 : static void GWKGeneralCaseThread( void* pData)
3360 :
3361 : {
3362 193 : GWKJobStruct* psJob = (GWKJobStruct*) pData;
3363 193 : GDALWarpKernel *poWK = psJob->poWK;
3364 193 : int iYMin = psJob->iYMin;
3365 193 : int iYMax = psJob->iYMax;
3366 :
3367 : int iDstY;
3368 193 : int nDstXSize = poWK->nDstXSize;
3369 193 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
3370 :
3371 : /* -------------------------------------------------------------------- */
3372 : /* Allocate x,y,z coordinate arrays for transformation ... one */
3373 : /* scanlines worth of positions. */
3374 : /* -------------------------------------------------------------------- */
3375 : double *padfX, *padfY, *padfZ;
3376 : int *pabSuccess;
3377 :
3378 193 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3379 193 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3380 193 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3381 193 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
3382 :
3383 193 : GWKResampleWrkStruct* psWrkStruct = NULL;
3384 193 : if (poWK->eResample == GRA_CubicSpline
3385 : || poWK->eResample == GRA_Lanczos )
3386 : {
3387 87 : psWrkStruct = GWKResampleCreateWrkStruct(poWK);
3388 : }
3389 :
3390 : /* ==================================================================== */
3391 : /* Loop over output lines. */
3392 : /* ==================================================================== */
3393 1917 : for( iDstY = iYMin; iDstY < iYMax; iDstY++ )
3394 : {
3395 : int iDstX;
3396 :
3397 : /* -------------------------------------------------------------------- */
3398 : /* Setup points to transform to source image space. */
3399 : /* -------------------------------------------------------------------- */
3400 533172 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3401 : {
3402 531448 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
3403 531448 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
3404 531448 : padfZ[iDstX] = 0.0;
3405 : }
3406 :
3407 : /* -------------------------------------------------------------------- */
3408 : /* Transform the points from destination pixel/line coordinates */
3409 : /* to source pixel/line coordinates. */
3410 : /* -------------------------------------------------------------------- */
3411 : poWK->pfnTransformer( psJob->pTransformerArg, TRUE, nDstXSize,
3412 1724 : padfX, padfY, padfZ, pabSuccess );
3413 :
3414 : /* ==================================================================== */
3415 : /* Loop over pixels in output scanline. */
3416 : /* ==================================================================== */
3417 533172 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3418 : {
3419 : int iDstOffset;
3420 :
3421 531448 : COMPUTE_iSrcOffset(pabSuccess, iDstX, padfX, padfY, poWK, nSrcXSize, nSrcYSize);
3422 :
3423 : /* -------------------------------------------------------------------- */
3424 : /* Do not try to apply transparent/invalid source pixels to the */
3425 : /* destination. This currently ignores the multi-pixel input */
3426 : /* of bilinear and cubic resamples. */
3427 : /* -------------------------------------------------------------------- */
3428 284907 : double dfDensity = 1.0;
3429 :
3430 284907 : if( poWK->pafUnifiedSrcDensity != NULL )
3431 : {
3432 1328 : dfDensity = poWK->pafUnifiedSrcDensity[iSrcOffset];
3433 1328 : if( dfDensity < 0.00001 )
3434 1299 : continue;
3435 : }
3436 :
3437 283608 : if( poWK->panUnifiedSrcValid != NULL
3438 0 : && !(poWK->panUnifiedSrcValid[iSrcOffset>>5]
3439 : & (0x01 << (iSrcOffset & 0x1f))) )
3440 0 : continue;
3441 :
3442 : /* ==================================================================== */
3443 : /* Loop processing each band. */
3444 : /* ==================================================================== */
3445 : int iBand;
3446 283608 : int bHasFoundDensity = FALSE;
3447 :
3448 283608 : iDstOffset = iDstX + iDstY * nDstXSize;
3449 573003 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
3450 : {
3451 289395 : double dfBandDensity = 0.0;
3452 289395 : double dfValueReal = 0.0;
3453 289395 : double dfValueImag = 0.0;
3454 :
3455 : /* -------------------------------------------------------------------- */
3456 : /* Collect the source value. */
3457 : /* -------------------------------------------------------------------- */
3458 290674 : if ( poWK->eResample == GRA_NearestNeighbour ||
3459 : nSrcXSize == 1 || nSrcYSize == 1)
3460 : {
3461 : GWKGetPixelValue( poWK, iBand, iSrcOffset,
3462 1279 : &dfBandDensity, &dfValueReal, &dfValueImag );
3463 : }
3464 288116 : else if ( poWK->eResample == GRA_Bilinear )
3465 : {
3466 : GWKBilinearResample( poWK, iBand,
3467 8174 : padfX[iDstX]-poWK->nSrcXOff,
3468 8174 : padfY[iDstX]-poWK->nSrcYOff,
3469 : &dfBandDensity,
3470 16348 : &dfValueReal, &dfValueImag );
3471 : }
3472 279942 : else if ( poWK->eResample == GRA_Cubic )
3473 : {
3474 : GWKCubicResample( poWK, iBand,
3475 558 : padfX[iDstX]-poWK->nSrcXOff,
3476 558 : padfY[iDstX]-poWK->nSrcYOff,
3477 : &dfBandDensity,
3478 1116 : &dfValueReal, &dfValueImag );
3479 : }
3480 279384 : else if ( poWK->eResample == GRA_CubicSpline
3481 : || poWK->eResample == GRA_Lanczos )
3482 : {
3483 : psWrkStruct->pfnGWKResample( poWK, iBand,
3484 279384 : padfX[iDstX]-poWK->nSrcXOff,
3485 279384 : padfY[iDstX]-poWK->nSrcYOff,
3486 : &dfBandDensity,
3487 838152 : &dfValueReal, &dfValueImag, psWrkStruct );
3488 : }
3489 :
3490 :
3491 : // If we didn't find any valid inputs skip to next band.
3492 289395 : if ( dfBandDensity < 0.0000000001 )
3493 0 : continue;
3494 :
3495 289395 : bHasFoundDensity = TRUE;
3496 :
3497 : /* -------------------------------------------------------------------- */
3498 : /* We have a computed value from the source. Now apply it to */
3499 : /* the destination pixel. */
3500 : /* -------------------------------------------------------------------- */
3501 : GWKSetPixelValue( poWK, iBand, iDstOffset,
3502 : dfBandDensity,
3503 289395 : dfValueReal, dfValueImag );
3504 :
3505 : }
3506 :
3507 283608 : if (!bHasFoundDensity)
3508 0 : continue;
3509 :
3510 : /* -------------------------------------------------------------------- */
3511 : /* Update destination density/validity masks. */
3512 : /* -------------------------------------------------------------------- */
3513 283608 : GWKOverlayDensity( poWK, iDstOffset, dfDensity );
3514 :
3515 283608 : if( poWK->panDstValid != NULL )
3516 : {
3517 : poWK->panDstValid[iDstOffset>>5] |=
3518 0 : 0x01 << (iDstOffset & 0x1f);
3519 : }
3520 :
3521 : } /* Next iDstX */
3522 :
3523 : /* -------------------------------------------------------------------- */
3524 : /* Report progress to the user, and optionally cancel out. */
3525 : /* -------------------------------------------------------------------- */
3526 1724 : if (psJob->pfnProgress(psJob))
3527 0 : break;
3528 : }
3529 :
3530 : /* -------------------------------------------------------------------- */
3531 : /* Cleanup and return. */
3532 : /* -------------------------------------------------------------------- */
3533 193 : CPLFree( padfX );
3534 193 : CPLFree( padfY );
3535 193 : CPLFree( padfZ );
3536 193 : CPLFree( pabSuccess );
3537 193 : if (psWrkStruct)
3538 87 : GWKResampleDeleteWrkStruct(psWrkStruct);
3539 193 : }
3540 :
3541 : /************************************************************************/
3542 : /* GWKNearestNoMasksByte() */
3543 : /* */
3544 : /* Case for 8bit input data with nearest neighbour resampling */
3545 : /* without concerning about masking. Should be as fast as */
3546 : /* possible for this particular transformation type. */
3547 : /************************************************************************/
3548 :
3549 : static void GWKNearestNoMasksByteThread(void* pData);
3550 :
3551 255 : static CPLErr GWKNearestNoMasksByte( GDALWarpKernel *poWK )
3552 : {
3553 255 : return GWKRun( poWK, "GWKNearestNoMasksByte", GWKNearestNoMasksByteThread );
3554 : }
3555 :
3556 255 : static void GWKNearestNoMasksByteThread(void* pData)
3557 : {
3558 255 : GWKJobStruct* psJob = (GWKJobStruct*) pData;
3559 255 : GDALWarpKernel *poWK = psJob->poWK;
3560 255 : int iYMin = psJob->iYMin;
3561 255 : int iYMax = psJob->iYMax;
3562 :
3563 : int iDstY;
3564 255 : int nDstXSize = poWK->nDstXSize;
3565 255 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
3566 :
3567 : /* -------------------------------------------------------------------- */
3568 : /* Allocate x,y,z coordinate arrays for transformation ... one */
3569 : /* scanlines worth of positions. */
3570 : /* -------------------------------------------------------------------- */
3571 : double *padfX, *padfY, *padfZ;
3572 : int *pabSuccess;
3573 :
3574 255 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3575 255 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3576 255 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3577 255 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
3578 :
3579 : /* ==================================================================== */
3580 : /* Loop over output lines. */
3581 : /* ==================================================================== */
3582 30684 : for( iDstY = iYMin; iDstY < iYMax; iDstY++ )
3583 : {
3584 : int iDstX;
3585 :
3586 : /* -------------------------------------------------------------------- */
3587 : /* Setup points to transform to source image space. */
3588 : /* -------------------------------------------------------------------- */
3589 8647071 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3590 : {
3591 8616642 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
3592 8616642 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
3593 8616642 : padfZ[iDstX] = 0.0;
3594 : }
3595 :
3596 : /* -------------------------------------------------------------------- */
3597 : /* Transform the points from destination pixel/line coordinates */
3598 : /* to source pixel/line coordinates. */
3599 : /* -------------------------------------------------------------------- */
3600 : poWK->pfnTransformer( psJob->pTransformerArg, TRUE, nDstXSize,
3601 30429 : padfX, padfY, padfZ, pabSuccess );
3602 :
3603 : /* ==================================================================== */
3604 : /* Loop over pixels in output scanline. */
3605 : /* ==================================================================== */
3606 8647071 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3607 : {
3608 8616642 : COMPUTE_iSrcOffset(pabSuccess, iDstX, padfX, padfY, poWK, nSrcXSize, nSrcYSize);
3609 :
3610 : /* ==================================================================== */
3611 : /* Loop processing each band. */
3612 : /* ==================================================================== */
3613 : int iBand;
3614 : int iDstOffset;
3615 :
3616 8331963 : iDstOffset = iDstX + iDstY * nDstXSize;
3617 :
3618 20252470 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
3619 : {
3620 11920507 : poWK->papabyDstImage[iBand][iDstOffset] =
3621 11920507 : poWK->papabySrcImage[iBand][iSrcOffset];
3622 : }
3623 : }
3624 :
3625 : /* -------------------------------------------------------------------- */
3626 : /* Report progress to the user, and optionally cancel out. */
3627 : /* -------------------------------------------------------------------- */
3628 30429 : if (psJob->pfnProgress(psJob))
3629 0 : break;
3630 : }
3631 :
3632 : /* -------------------------------------------------------------------- */
3633 : /* Cleanup and return. */
3634 : /* -------------------------------------------------------------------- */
3635 255 : CPLFree( padfX );
3636 255 : CPLFree( padfY );
3637 255 : CPLFree( padfZ );
3638 255 : CPLFree( pabSuccess );
3639 255 : }
3640 :
3641 : /************************************************************************/
3642 : /* GWKBilinearNoMasksByte() */
3643 : /* */
3644 : /* Case for 8bit input data with bilinear resampling without */
3645 : /* concerning about masking. Should be as fast as possible */
3646 : /* for this particular transformation type. */
3647 : /************************************************************************/
3648 :
3649 : static void GWKBilinearNoMasksByteThread(void* pData);
3650 :
3651 16 : static CPLErr GWKBilinearNoMasksByte( GDALWarpKernel *poWK )
3652 : {
3653 16 : return GWKRun( poWK, "GWKBilinearNoMasksByte", GWKBilinearNoMasksByteThread );
3654 : }
3655 :
3656 17 : static void GWKBilinearNoMasksByteThread(void* pData)
3657 : {
3658 17 : GWKJobStruct* psJob = (GWKJobStruct*) pData;
3659 17 : GDALWarpKernel *poWK = psJob->poWK;
3660 17 : int iYMin = psJob->iYMin;
3661 17 : int iYMax = psJob->iYMax;
3662 :
3663 : int iDstY;
3664 17 : int nDstXSize = poWK->nDstXSize;
3665 17 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
3666 :
3667 : /* -------------------------------------------------------------------- */
3668 : /* Allocate x,y,z coordinate arrays for transformation ... one */
3669 : /* scanlines worth of positions. */
3670 : /* -------------------------------------------------------------------- */
3671 : double *padfX, *padfY, *padfZ;
3672 : int *pabSuccess;
3673 :
3674 17 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3675 17 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3676 17 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3677 17 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
3678 :
3679 : /* ==================================================================== */
3680 : /* Loop over output lines. */
3681 : /* ==================================================================== */
3682 768 : for( iDstY = iYMin; iDstY < iYMax; iDstY++ )
3683 : {
3684 : int iDstX;
3685 :
3686 : /* -------------------------------------------------------------------- */
3687 : /* Setup points to transform to source image space. */
3688 : /* -------------------------------------------------------------------- */
3689 331234 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3690 : {
3691 330480 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
3692 330480 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
3693 330480 : padfZ[iDstX] = 0.0;
3694 : }
3695 :
3696 : /* -------------------------------------------------------------------- */
3697 : /* Transform the points from destination pixel/line coordinates */
3698 : /* to source pixel/line coordinates. */
3699 : /* -------------------------------------------------------------------- */
3700 : poWK->pfnTransformer( psJob->pTransformerArg, TRUE, nDstXSize,
3701 754 : padfX, padfY, padfZ, pabSuccess );
3702 :
3703 : /* ==================================================================== */
3704 : /* Loop over pixels in output scanline. */
3705 : /* ==================================================================== */
3706 331234 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3707 : {
3708 330480 : COMPUTE_iSrcOffset(pabSuccess, iDstX, padfX, padfY, poWK, nSrcXSize, nSrcYSize);
3709 :
3710 : /* ==================================================================== */
3711 : /* Loop processing each band. */
3712 : /* ==================================================================== */
3713 : int iBand;
3714 : int iDstOffset;
3715 :
3716 279344 : iDstOffset = iDstX + iDstY * nDstXSize;
3717 :
3718 558688 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
3719 : {
3720 : GWKBilinearResampleNoMasksByte( poWK, iBand,
3721 279344 : padfX[iDstX]-poWK->nSrcXOff,
3722 279344 : padfY[iDstX]-poWK->nSrcYOff,
3723 838032 : &poWK->papabyDstImage[iBand][iDstOffset] );
3724 : }
3725 : }
3726 :
3727 : /* -------------------------------------------------------------------- */
3728 : /* Report progress to the user, and optionally cancel out. */
3729 : /* -------------------------------------------------------------------- */
3730 754 : if (psJob->pfnProgress(psJob))
3731 3 : break;
3732 : }
3733 :
3734 : /* -------------------------------------------------------------------- */
3735 : /* Cleanup and return. */
3736 : /* -------------------------------------------------------------------- */
3737 17 : CPLFree( padfX );
3738 17 : CPLFree( padfY );
3739 17 : CPLFree( padfZ );
3740 17 : CPLFree( pabSuccess );
3741 17 : }
3742 :
3743 : /************************************************************************/
3744 : /* GWKCubicNoMasksByte() */
3745 : /* */
3746 : /* Case for 8bit input data with cubic resampling without */
3747 : /* concerning about masking. Should be as fast as possible */
3748 : /* for this particular transformation type. */
3749 : /************************************************************************/
3750 :
3751 : static void GWKCubicNoMasksByteThread(void* pData);
3752 :
3753 10 : static CPLErr GWKCubicNoMasksByte( GDALWarpKernel *poWK )
3754 : {
3755 10 : return GWKRun( poWK, "GWKCubicNoMasksByte", GWKCubicNoMasksByteThread );
3756 : }
3757 :
3758 10 : static void GWKCubicNoMasksByteThread( void* pData )
3759 :
3760 : {
3761 10 : GWKJobStruct* psJob = (GWKJobStruct*) pData;
3762 10 : GDALWarpKernel *poWK = psJob->poWK;
3763 10 : int iYMin = psJob->iYMin;
3764 10 : int iYMax = psJob->iYMax;
3765 :
3766 : int iDstY;
3767 10 : int nDstXSize = poWK->nDstXSize;
3768 10 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
3769 :
3770 : /* -------------------------------------------------------------------- */
3771 : /* Allocate x,y,z coordinate arrays for transformation ... one */
3772 : /* scanlines worth of positions. */
3773 : /* -------------------------------------------------------------------- */
3774 : double *padfX, *padfY, *padfZ;
3775 : int *pabSuccess;
3776 :
3777 10 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3778 10 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3779 10 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3780 10 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
3781 :
3782 : /* ==================================================================== */
3783 : /* Loop over output lines. */
3784 : /* ==================================================================== */
3785 703 : for( iDstY = iYMin; iDstY < iYMax; iDstY++ )
3786 : {
3787 : int iDstX;
3788 :
3789 : /* -------------------------------------------------------------------- */
3790 : /* Setup points to transform to source image space. */
3791 : /* -------------------------------------------------------------------- */
3792 330036 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3793 : {
3794 329343 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
3795 329343 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
3796 329343 : padfZ[iDstX] = 0.0;
3797 : }
3798 :
3799 : /* -------------------------------------------------------------------- */
3800 : /* Transform the points from destination pixel/line coordinates */
3801 : /* to source pixel/line coordinates. */
3802 : /* -------------------------------------------------------------------- */
3803 : poWK->pfnTransformer( psJob->pTransformerArg, TRUE, nDstXSize,
3804 693 : padfX, padfY, padfZ, pabSuccess );
3805 :
3806 : /* ==================================================================== */
3807 : /* Loop over pixels in output scanline. */
3808 : /* ==================================================================== */
3809 330036 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3810 : {
3811 329343 : COMPUTE_iSrcOffset(pabSuccess, iDstX, padfX, padfY, poWK, nSrcXSize, nSrcYSize);
3812 :
3813 : /* ==================================================================== */
3814 : /* Loop processing each band. */
3815 : /* ==================================================================== */
3816 : int iBand;
3817 : int iDstOffset;
3818 :
3819 278207 : iDstOffset = iDstX + iDstY * nDstXSize;
3820 :
3821 556414 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
3822 : {
3823 : GWKCubicResampleNoMasksByte( poWK, iBand,
3824 278207 : padfX[iDstX]-poWK->nSrcXOff,
3825 278207 : padfY[iDstX]-poWK->nSrcYOff,
3826 834621 : &poWK->papabyDstImage[iBand][iDstOffset] );
3827 : }
3828 : }
3829 :
3830 : /* -------------------------------------------------------------------- */
3831 : /* Report progress to the user, and optionally cancel out. */
3832 : /* -------------------------------------------------------------------- */
3833 693 : if (psJob->pfnProgress(psJob))
3834 0 : break;
3835 : }
3836 :
3837 : /* -------------------------------------------------------------------- */
3838 : /* Cleanup and return. */
3839 : /* -------------------------------------------------------------------- */
3840 10 : CPLFree( padfX );
3841 10 : CPLFree( padfY );
3842 10 : CPLFree( padfZ );
3843 10 : CPLFree( pabSuccess );
3844 10 : }
3845 :
3846 : /************************************************************************/
3847 : /* GWKCubicSplineNoMasksByte() */
3848 : /* */
3849 : /* Case for 8bit input data with cubic spline resampling without */
3850 : /* concerning about masking. Should be as fast as possible */
3851 : /* for this particular transformation type. */
3852 : /************************************************************************/
3853 :
3854 : static void GWKCubicSplineNoMasksByteThread(void* pData);
3855 :
3856 10 : static CPLErr GWKCubicSplineNoMasksByte( GDALWarpKernel *poWK )
3857 : {
3858 10 : return GWKRun( poWK, "GWKCubicSplineNoMasksByte", GWKCubicSplineNoMasksByteThread );
3859 : }
3860 :
3861 10 : static void GWKCubicSplineNoMasksByteThread( void* pData )
3862 :
3863 : {
3864 10 : GWKJobStruct* psJob = (GWKJobStruct*) pData;
3865 10 : GDALWarpKernel *poWK = psJob->poWK;
3866 10 : int iYMin = psJob->iYMin;
3867 10 : int iYMax = psJob->iYMax;
3868 :
3869 : int iDstY;
3870 10 : int nDstXSize = poWK->nDstXSize;
3871 10 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
3872 :
3873 : /* -------------------------------------------------------------------- */
3874 : /* Allocate x,y,z coordinate arrays for transformation ... one */
3875 : /* scanlines worth of positions. */
3876 : /* -------------------------------------------------------------------- */
3877 : double *padfX, *padfY, *padfZ;
3878 : int *pabSuccess;
3879 :
3880 10 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3881 10 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3882 10 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3883 10 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
3884 :
3885 10 : int nXRadius = poWK->nXRadius;
3886 10 : double *padfBSpline = (double *)CPLCalloc( nXRadius * 2, sizeof(double) );
3887 :
3888 : /* ==================================================================== */
3889 : /* Loop over output lines. */
3890 : /* ==================================================================== */
3891 703 : for( iDstY = iYMin; iDstY < iYMax; iDstY++ )
3892 : {
3893 : int iDstX;
3894 :
3895 : /* -------------------------------------------------------------------- */
3896 : /* Setup points to transform to source image space. */
3897 : /* -------------------------------------------------------------------- */
3898 330036 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3899 : {
3900 329343 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
3901 329343 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
3902 329343 : padfZ[iDstX] = 0.0;
3903 : }
3904 :
3905 : /* -------------------------------------------------------------------- */
3906 : /* Transform the points from destination pixel/line coordinates */
3907 : /* to source pixel/line coordinates. */
3908 : /* -------------------------------------------------------------------- */
3909 : poWK->pfnTransformer( psJob->pTransformerArg, TRUE, nDstXSize,
3910 693 : padfX, padfY, padfZ, pabSuccess );
3911 :
3912 : /* ==================================================================== */
3913 : /* Loop over pixels in output scanline. */
3914 : /* ==================================================================== */
3915 330036 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3916 : {
3917 329343 : COMPUTE_iSrcOffset(pabSuccess, iDstX, padfX, padfY, poWK, nSrcXSize, nSrcYSize);
3918 :
3919 : /* ==================================================================== */
3920 : /* Loop processing each band. */
3921 : /* ==================================================================== */
3922 : int iBand;
3923 : int iDstOffset;
3924 :
3925 278207 : iDstOffset = iDstX + iDstY * nDstXSize;
3926 :
3927 556414 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
3928 : {
3929 : GWKCubicSplineResampleNoMasksByte( poWK, iBand,
3930 278207 : padfX[iDstX]-poWK->nSrcXOff,
3931 278207 : padfY[iDstX]-poWK->nSrcYOff,
3932 278207 : &poWK->papabyDstImage[iBand][iDstOffset],
3933 834621 : padfBSpline);
3934 : }
3935 : }
3936 :
3937 : /* -------------------------------------------------------------------- */
3938 : /* Report progress to the user, and optionally cancel out. */
3939 : /* -------------------------------------------------------------------- */
3940 693 : if (psJob->pfnProgress(psJob))
3941 0 : break;
3942 : }
3943 :
3944 : /* -------------------------------------------------------------------- */
3945 : /* Cleanup and return. */
3946 : /* -------------------------------------------------------------------- */
3947 10 : CPLFree( padfX );
3948 10 : CPLFree( padfY );
3949 10 : CPLFree( padfZ );
3950 10 : CPLFree( pabSuccess );
3951 10 : CPLFree( padfBSpline );
3952 10 : }
3953 :
3954 : /************************************************************************/
3955 : /* GWKNearestByte() */
3956 : /* */
3957 : /* Case for 8bit input data with nearest neighbour resampling */
3958 : /* using valid flags. Should be as fast as possible for this */
3959 : /* particular transformation type. */
3960 : /************************************************************************/
3961 :
3962 : static void GWKNearestByteThread(void* pData);
3963 :
3964 28 : static CPLErr GWKNearestByte( GDALWarpKernel *poWK )
3965 : {
3966 28 : return GWKRun( poWK, "GWKNearestByte", GWKNearestByteThread );
3967 : }
3968 :
3969 32 : static void GWKNearestByteThread( void* pData )
3970 :
3971 : {
3972 32 : GWKJobStruct* psJob = (GWKJobStruct*) pData;
3973 32 : GDALWarpKernel *poWK = psJob->poWK;
3974 32 : int iYMin = psJob->iYMin;
3975 32 : int iYMax = psJob->iYMax;
3976 :
3977 : int iDstY;
3978 32 : int nDstXSize = poWK->nDstXSize;
3979 32 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
3980 :
3981 : /* -------------------------------------------------------------------- */
3982 : /* Allocate x,y,z coordinate arrays for transformation ... one */
3983 : /* scanlines worth of positions. */
3984 : /* -------------------------------------------------------------------- */
3985 : double *padfX, *padfY, *padfZ;
3986 : int *pabSuccess;
3987 :
3988 32 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3989 32 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3990 32 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3991 32 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
3992 :
3993 : /* ==================================================================== */
3994 : /* Loop over output lines. */
3995 : /* ==================================================================== */
3996 3572 : for( iDstY = iYMin; iDstY < iYMax; iDstY++ )
3997 : {
3998 : int iDstX;
3999 :
4000 : /* -------------------------------------------------------------------- */
4001 : /* Setup points to transform to source image space. */
4002 : /* -------------------------------------------------------------------- */
4003 1596029 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4004 : {
4005 1592489 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
4006 1592489 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
4007 1592489 : padfZ[iDstX] = 0.0;
4008 : }
4009 :
4010 : /* -------------------------------------------------------------------- */
4011 : /* Transform the points from destination pixel/line coordinates */
4012 : /* to source pixel/line coordinates. */
4013 : /* -------------------------------------------------------------------- */
4014 : poWK->pfnTransformer( psJob->pTransformerArg, TRUE, nDstXSize,
4015 3540 : padfX, padfY, padfZ, pabSuccess );
4016 :
4017 : /* ==================================================================== */
4018 : /* Loop over pixels in output scanline. */
4019 : /* ==================================================================== */
4020 1571262 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4021 : {
4022 : int iDstOffset;
4023 :
4024 1567723 : COMPUTE_iSrcOffset(pabSuccess, iDstX, padfX, padfY, poWK, nSrcXSize, nSrcYSize);
4025 :
4026 : /* -------------------------------------------------------------------- */
4027 : /* Do not try to apply invalid source pixels to the dest. */
4028 : /* -------------------------------------------------------------------- */
4029 732326 : if( poWK->panUnifiedSrcValid != NULL
4030 158681 : && !(poWK->panUnifiedSrcValid[iSrcOffset>>5]
4031 : & (0x01 << (iSrcOffset & 0x1f))) )
4032 32599 : continue;
4033 :
4034 : /* -------------------------------------------------------------------- */
4035 : /* Do not try to apply transparent source pixels to the destination.*/
4036 : /* -------------------------------------------------------------------- */
4037 541046 : double dfDensity = 1.0;
4038 :
4039 541046 : if( poWK->pafUnifiedSrcDensity != NULL )
4040 : {
4041 90000 : dfDensity = poWK->pafUnifiedSrcDensity[iSrcOffset];
4042 90000 : if( dfDensity < 0.00001 )
4043 74897 : continue;
4044 : }
4045 :
4046 : /* ==================================================================== */
4047 : /* Loop processing each band. */
4048 : /* ==================================================================== */
4049 : int iBand;
4050 :
4051 466149 : iDstOffset = iDstX + iDstY * nDstXSize;
4052 :
4053 1192046 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
4054 : {
4055 725901 : GByte bValue = 0;
4056 725901 : double dfBandDensity = 0.0;
4057 :
4058 : /* -------------------------------------------------------------------- */
4059 : /* Collect the source value. */
4060 : /* -------------------------------------------------------------------- */
4061 725901 : if ( GWKGetPixelByte( poWK, iBand, iSrcOffset, &dfBandDensity,
4062 : &bValue ) )
4063 : {
4064 725901 : if( dfBandDensity < 1.0 )
4065 : {
4066 1472 : if( dfBandDensity == 0.0 )
4067 : /* do nothing */;
4068 : else
4069 : {
4070 : /* let the general code take care of mixing */
4071 : GWKSetPixelValue( poWK, iBand, iDstOffset,
4072 : dfBandDensity, (double) bValue,
4073 1472 : 0.0 );
4074 : }
4075 : }
4076 : else
4077 : {
4078 724429 : poWK->papabyDstImage[iBand][iDstOffset] = bValue;
4079 : }
4080 : }
4081 : }
4082 :
4083 : /* -------------------------------------------------------------------- */
4084 : /* Mark this pixel valid/opaque in the output. */
4085 : /* -------------------------------------------------------------------- */
4086 466145 : GWKOverlayDensity( poWK, iDstOffset, dfDensity );
4087 :
4088 466145 : if( poWK->panDstValid != NULL )
4089 : {
4090 : poWK->panDstValid[iDstOffset>>5] |=
4091 1605 : 0x01 << (iDstOffset & 0x1f);
4092 : }
4093 : } /* Next iDstX */
4094 :
4095 : /* -------------------------------------------------------------------- */
4096 : /* Report progress to the user, and optionally cancel out. */
4097 : /* -------------------------------------------------------------------- */
4098 3539 : if (psJob->pfnProgress(psJob))
4099 0 : break;
4100 : } /* Next iDstY */
4101 :
4102 : /* -------------------------------------------------------------------- */
4103 : /* Cleanup and return. */
4104 : /* -------------------------------------------------------------------- */
4105 32 : CPLFree( padfX );
4106 32 : CPLFree( padfY );
4107 32 : CPLFree( padfZ );
4108 32 : CPLFree( pabSuccess );
4109 32 : }
4110 :
4111 : /************************************************************************/
4112 : /* GWKNearestNoMasksShort() */
4113 : /* */
4114 : /* Case for 16bit signed and unsigned integer input data with */
4115 : /* nearest neighbour resampling without concerning about masking. */
4116 : /* Should be as fast as possible for this particular */
4117 : /* transformation type. */
4118 : /************************************************************************/
4119 :
4120 : static void GWKNearestNoMasksShortThread(void* pData);
4121 :
4122 14 : static CPLErr GWKNearestNoMasksShort( GDALWarpKernel *poWK )
4123 : {
4124 14 : return GWKRun( poWK, "GWKNearestNoMasksShort", GWKNearestNoMasksShortThread );
4125 : }
4126 :
4127 14 : static void GWKNearestNoMasksShortThread( void* pData )
4128 :
4129 : {
4130 14 : GWKJobStruct* psJob = (GWKJobStruct*) pData;
4131 14 : GDALWarpKernel *poWK = psJob->poWK;
4132 14 : int iYMin = psJob->iYMin;
4133 14 : int iYMax = psJob->iYMax;
4134 :
4135 : int iDstY;
4136 14 : int nDstXSize = poWK->nDstXSize;
4137 14 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
4138 :
4139 : /* -------------------------------------------------------------------- */
4140 : /* Allocate x,y,z coordinate arrays for transformation ... one */
4141 : /* scanlines worth of positions. */
4142 : /* -------------------------------------------------------------------- */
4143 : double *padfX, *padfY, *padfZ;
4144 : int *pabSuccess;
4145 :
4146 14 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4147 14 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4148 14 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4149 14 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
4150 :
4151 : /* ==================================================================== */
4152 : /* Loop over output lines. */
4153 : /* ==================================================================== */
4154 700 : for( iDstY = iYMin; iDstY < iYMax; iDstY++ )
4155 : {
4156 : int iDstX;
4157 :
4158 : /* -------------------------------------------------------------------- */
4159 : /* Setup points to transform to source image space. */
4160 : /* -------------------------------------------------------------------- */
4161 328892 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4162 : {
4163 328206 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
4164 328206 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
4165 328206 : padfZ[iDstX] = 0.0;
4166 : }
4167 :
4168 : /* -------------------------------------------------------------------- */
4169 : /* Transform the points from destination pixel/line coordinates */
4170 : /* to source pixel/line coordinates. */
4171 : /* -------------------------------------------------------------------- */
4172 : poWK->pfnTransformer( psJob->pTransformerArg, TRUE, nDstXSize,
4173 686 : padfX, padfY, padfZ, pabSuccess );
4174 :
4175 : /* ==================================================================== */
4176 : /* Loop over pixels in output scanline. */
4177 : /* ==================================================================== */
4178 328892 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4179 : {
4180 : int iDstOffset;
4181 :
4182 328206 : COMPUTE_iSrcOffset(pabSuccess, iDstX, padfX, padfY, poWK, nSrcXSize, nSrcYSize);
4183 :
4184 : /* ==================================================================== */
4185 : /* Loop processing each band. */
4186 : /* ==================================================================== */
4187 : int iBand;
4188 :
4189 264948 : iDstOffset = iDstX + iDstY * nDstXSize;
4190 :
4191 529896 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
4192 : {
4193 264948 : ((GInt16 *)poWK->papabyDstImage[iBand])[iDstOffset] =
4194 264948 : ((GInt16 *)poWK->papabySrcImage[iBand])[iSrcOffset];
4195 : }
4196 : }
4197 :
4198 : /* -------------------------------------------------------------------- */
4199 : /* Report progress to the user, and optionally cancel out. */
4200 : /* -------------------------------------------------------------------- */
4201 686 : if (psJob->pfnProgress(psJob))
4202 0 : break;
4203 : }
4204 :
4205 : /* -------------------------------------------------------------------- */
4206 : /* Cleanup and return. */
4207 : /* -------------------------------------------------------------------- */
4208 14 : CPLFree( padfX );
4209 14 : CPLFree( padfY );
4210 14 : CPLFree( padfZ );
4211 14 : CPLFree( pabSuccess );
4212 14 : }
4213 :
4214 : /************************************************************************/
4215 : /* GWKBilinearNoMasksShort() */
4216 : /* */
4217 : /* Case for 16bit input data with cubic resampling without */
4218 : /* concerning about masking. Should be as fast as possible */
4219 : /* for this particular transformation type. */
4220 : /************************************************************************/
4221 :
4222 : static void GWKBilinearNoMasksShortThread(void* pData);
4223 :
4224 20 : static CPLErr GWKBilinearNoMasksShort( GDALWarpKernel *poWK )
4225 : {
4226 20 : return GWKRun( poWK, "GWKBilinearNoMasksShort", GWKBilinearNoMasksShortThread );
4227 : }
4228 :
4229 20 : static void GWKBilinearNoMasksShortThread( void* pData )
4230 :
4231 : {
4232 20 : GWKJobStruct* psJob = (GWKJobStruct*) pData;
4233 20 : GDALWarpKernel *poWK = psJob->poWK;
4234 20 : int iYMin = psJob->iYMin;
4235 20 : int iYMax = psJob->iYMax;
4236 :
4237 : int iDstY;
4238 20 : int nDstXSize = poWK->nDstXSize;
4239 20 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
4240 :
4241 : /* -------------------------------------------------------------------- */
4242 : /* Allocate x,y,z coordinate arrays for transformation ... one */
4243 : /* scanlines worth of positions. */
4244 : /* -------------------------------------------------------------------- */
4245 : double *padfX, *padfY, *padfZ;
4246 : int *pabSuccess;
4247 :
4248 20 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4249 20 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4250 20 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4251 20 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
4252 :
4253 : /* ==================================================================== */
4254 : /* Loop over output lines. */
4255 : /* ==================================================================== */
4256 1856 : for( iDstY = iYMin; iDstY < iYMax; iDstY++ )
4257 : {
4258 : int iDstX;
4259 :
4260 : /* -------------------------------------------------------------------- */
4261 : /* Setup points to transform to source image space. */
4262 : /* -------------------------------------------------------------------- */
4263 1707544 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4264 : {
4265 1705708 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
4266 1705708 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
4267 1705708 : padfZ[iDstX] = 0.0;
4268 : }
4269 :
4270 : /* -------------------------------------------------------------------- */
4271 : /* Transform the points from destination pixel/line coordinates */
4272 : /* to source pixel/line coordinates. */
4273 : /* -------------------------------------------------------------------- */
4274 : poWK->pfnTransformer( psJob->pTransformerArg, TRUE, nDstXSize,
4275 1836 : padfX, padfY, padfZ, pabSuccess );
4276 :
4277 : /* ==================================================================== */
4278 : /* Loop over pixels in output scanline. */
4279 : /* ==================================================================== */
4280 1707544 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4281 : {
4282 1705708 : COMPUTE_iSrcOffset(pabSuccess, iDstX, padfX, padfY, poWK, nSrcXSize, nSrcYSize);
4283 :
4284 : /* ==================================================================== */
4285 : /* Loop processing each band. */
4286 : /* ==================================================================== */
4287 : int iBand;
4288 : int iDstOffset;
4289 :
4290 1705708 : iDstOffset = iDstX + iDstY * nDstXSize;
4291 :
4292 3411416 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
4293 : {
4294 1705708 : GInt16 iValue = 0;
4295 : GWKBilinearResampleNoMasksShort( poWK, iBand,
4296 1705708 : padfX[iDstX]-poWK->nSrcXOff,
4297 1705708 : padfY[iDstX]-poWK->nSrcYOff,
4298 3411416 : &iValue );
4299 1705708 : ((GInt16 *)poWK->papabyDstImage[iBand])[iDstOffset] = iValue;
4300 : }
4301 : }
4302 :
4303 : /* -------------------------------------------------------------------- */
4304 : /* Report progress to the user, and optionally cancel out. */
4305 : /* -------------------------------------------------------------------- */
4306 1836 : if (psJob->pfnProgress(psJob))
4307 0 : break;
4308 : }
4309 :
4310 : /* -------------------------------------------------------------------- */
4311 : /* Cleanup and return. */
4312 : /* -------------------------------------------------------------------- */
4313 20 : CPLFree( padfX );
4314 20 : CPLFree( padfY );
4315 20 : CPLFree( padfZ );
4316 20 : CPLFree( pabSuccess );
4317 20 : }
4318 :
4319 : /************************************************************************/
4320 : /* GWKCubicNoMasksShort() */
4321 : /* */
4322 : /* Case for 16bit input data with cubic resampling without */
4323 : /* concerning about masking. Should be as fast as possible */
4324 : /* for this particular transformation type. */
4325 : /************************************************************************/
4326 :
4327 : static void GWKCubicNoMasksShortThread(void* pData);
4328 :
4329 8 : static CPLErr GWKCubicNoMasksShort( GDALWarpKernel *poWK )
4330 : {
4331 8 : return GWKRun( poWK, "GWKCubicNoMasksShort", GWKCubicNoMasksShortThread );
4332 : }
4333 :
4334 8 : static void GWKCubicNoMasksShortThread( void* pData )
4335 :
4336 : {
4337 8 : GWKJobStruct* psJob = (GWKJobStruct*) pData;
4338 8 : GDALWarpKernel *poWK = psJob->poWK;
4339 8 : int iYMin = psJob->iYMin;
4340 8 : int iYMax = psJob->iYMax;
4341 :
4342 : int iDstY;
4343 8 : int nDstXSize = poWK->nDstXSize;
4344 8 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
4345 :
4346 : /* -------------------------------------------------------------------- */
4347 : /* Allocate x,y,z coordinate arrays for transformation ... one */
4348 : /* scanlines worth of positions. */
4349 : /* -------------------------------------------------------------------- */
4350 : double *padfX, *padfY, *padfZ;
4351 : int *pabSuccess;
4352 :
4353 8 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4354 8 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4355 8 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4356 8 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
4357 :
4358 : /* ==================================================================== */
4359 : /* Loop over output lines. */
4360 : /* ==================================================================== */
4361 533 : for( iDstY = iYMin; iDstY < iYMax; iDstY++ )
4362 : {
4363 : int iDstX;
4364 :
4365 : /* -------------------------------------------------------------------- */
4366 : /* Setup points to transform to source image space. */
4367 : /* -------------------------------------------------------------------- */
4368 262732 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4369 : {
4370 262207 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
4371 262207 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
4372 262207 : padfZ[iDstX] = 0.0;
4373 : }
4374 :
4375 : /* -------------------------------------------------------------------- */
4376 : /* Transform the points from destination pixel/line coordinates */
4377 : /* to source pixel/line coordinates. */
4378 : /* -------------------------------------------------------------------- */
4379 : poWK->pfnTransformer( psJob->pTransformerArg, TRUE, nDstXSize,
4380 525 : padfX, padfY, padfZ, pabSuccess );
4381 :
4382 : /* ==================================================================== */
4383 : /* Loop over pixels in output scanline. */
4384 : /* ==================================================================== */
4385 262732 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4386 : {
4387 262207 : COMPUTE_iSrcOffset(pabSuccess, iDstX, padfX, padfY, poWK, nSrcXSize, nSrcYSize);
4388 :
4389 : /* ==================================================================== */
4390 : /* Loop processing each band. */
4391 : /* ==================================================================== */
4392 : int iBand;
4393 : int iDstOffset;
4394 :
4395 262207 : iDstOffset = iDstX + iDstY * nDstXSize;
4396 :
4397 524414 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
4398 : {
4399 262207 : GInt16 iValue = 0;
4400 : GWKCubicResampleNoMasksShort( poWK, iBand,
4401 262207 : padfX[iDstX]-poWK->nSrcXOff,
4402 262207 : padfY[iDstX]-poWK->nSrcYOff,
4403 524414 : &iValue );
4404 262207 : ((GInt16 *)poWK->papabyDstImage[iBand])[iDstOffset] = iValue;
4405 : }
4406 : }
4407 :
4408 : /* -------------------------------------------------------------------- */
4409 : /* Report progress to the user, and optionally cancel out. */
4410 : /* -------------------------------------------------------------------- */
4411 525 : if (psJob->pfnProgress(psJob))
4412 0 : break;
4413 : }
4414 :
4415 : /* -------------------------------------------------------------------- */
4416 : /* Cleanup and return. */
4417 : /* -------------------------------------------------------------------- */
4418 8 : CPLFree( padfX );
4419 8 : CPLFree( padfY );
4420 8 : CPLFree( padfZ );
4421 8 : CPLFree( pabSuccess );
4422 8 : }
4423 :
4424 : /************************************************************************/
4425 : /* GWKCubicSplineNoMasksShort() */
4426 : /* */
4427 : /* Case for 16bit input data with cubic resampling without */
4428 : /* concerning about masking. Should be as fast as possible */
4429 : /* for this particular transformation type. */
4430 : /************************************************************************/
4431 :
4432 : static void GWKCubicSplineNoMasksShortThread(void* pData);
4433 :
4434 8 : static CPLErr GWKCubicSplineNoMasksShort( GDALWarpKernel *poWK )
4435 : {
4436 8 : return GWKRun( poWK, "GWKCubicSplineNoMasksShort", GWKCubicSplineNoMasksShortThread );
4437 : }
4438 :
4439 8 : static void GWKCubicSplineNoMasksShortThread( void* pData )
4440 :
4441 : {
4442 8 : GWKJobStruct* psJob = (GWKJobStruct*) pData;
4443 8 : GDALWarpKernel *poWK = psJob->poWK;
4444 8 : int iYMin = psJob->iYMin;
4445 8 : int iYMax = psJob->iYMax;
4446 :
4447 : int iDstY;
4448 8 : int nDstXSize = poWK->nDstXSize;
4449 8 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
4450 :
4451 : /* -------------------------------------------------------------------- */
4452 : /* Allocate x,y,z coordinate arrays for transformation ... one */
4453 : /* scanlines worth of positions. */
4454 : /* -------------------------------------------------------------------- */
4455 : double *padfX, *padfY, *padfZ;
4456 : int *pabSuccess;
4457 :
4458 8 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4459 8 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4460 8 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4461 8 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
4462 :
4463 8 : int nXRadius = poWK->nXRadius;
4464 : // Make space to save weights
4465 8 : double *padfBSpline = (double *)CPLCalloc( nXRadius * 2, sizeof(double) );
4466 :
4467 : /* ==================================================================== */
4468 : /* Loop over output lines. */
4469 : /* ==================================================================== */
4470 533 : for( iDstY = iYMin; iDstY < iYMax; iDstY++ )
4471 : {
4472 : int iDstX;
4473 :
4474 : /* -------------------------------------------------------------------- */
4475 : /* Setup points to transform to source image space. */
4476 : /* -------------------------------------------------------------------- */
4477 262732 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4478 : {
4479 262207 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
4480 262207 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
4481 262207 : padfZ[iDstX] = 0.0;
4482 : }
4483 :
4484 : /* -------------------------------------------------------------------- */
4485 : /* Transform the points from destination pixel/line coordinates */
4486 : /* to source pixel/line coordinates. */
4487 : /* -------------------------------------------------------------------- */
4488 : poWK->pfnTransformer( psJob->pTransformerArg, TRUE, nDstXSize,
4489 525 : padfX, padfY, padfZ, pabSuccess );
4490 :
4491 : /* ==================================================================== */
4492 : /* Loop over pixels in output scanline. */
4493 : /* ==================================================================== */
4494 262732 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4495 : {
4496 262207 : COMPUTE_iSrcOffset(pabSuccess, iDstX, padfX, padfY, poWK, nSrcXSize, nSrcYSize);
4497 :
4498 : /* ==================================================================== */
4499 : /* Loop processing each band. */
4500 : /* ==================================================================== */
4501 : int iBand;
4502 : int iDstOffset;
4503 :
4504 262207 : iDstOffset = iDstX + iDstY * nDstXSize;
4505 :
4506 524414 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
4507 : {
4508 262207 : GInt16 iValue = 0;
4509 : GWKCubicSplineResampleNoMasksShort( poWK, iBand,
4510 262207 : padfX[iDstX]-poWK->nSrcXOff,
4511 262207 : padfY[iDstX]-poWK->nSrcYOff,
4512 : &iValue,
4513 524414 : padfBSpline);
4514 262207 : ((GInt16 *)poWK->papabyDstImage[iBand])[iDstOffset] = iValue;
4515 : }
4516 : }
4517 :
4518 : /* -------------------------------------------------------------------- */
4519 : /* Report progress to the user, and optionally cancel out. */
4520 : /* -------------------------------------------------------------------- */
4521 525 : if (psJob->pfnProgress(psJob))
4522 0 : break;
4523 : }
4524 :
4525 : /* -------------------------------------------------------------------- */
4526 : /* Cleanup and return. */
4527 : /* -------------------------------------------------------------------- */
4528 8 : CPLFree( padfX );
4529 8 : CPLFree( padfY );
4530 8 : CPLFree( padfZ );
4531 8 : CPLFree( pabSuccess );
4532 8 : CPLFree( padfBSpline );
4533 8 : }
4534 :
4535 : /************************************************************************/
4536 : /* GWKNearestShort() */
4537 : /* */
4538 : /* Case for 32bit float input data with nearest neighbour */
4539 : /* resampling using valid flags. Should be as fast as possible */
4540 : /* for this particular transformation type. */
4541 : /************************************************************************/
4542 :
4543 : static void GWKNearestShortThread(void* pData);
4544 :
4545 7 : static CPLErr GWKNearestShort( GDALWarpKernel *poWK )
4546 : {
4547 7 : return GWKRun( poWK, "GWKNearestShort", GWKNearestShortThread );
4548 : }
4549 :
4550 7 : static void GWKNearestShortThread(void* pData)
4551 : {
4552 7 : GWKJobStruct* psJob = (GWKJobStruct*) pData;
4553 7 : GDALWarpKernel *poWK = psJob->poWK;
4554 7 : int iYMin = psJob->iYMin;
4555 7 : int iYMax = psJob->iYMax;
4556 :
4557 : int iDstY;
4558 7 : int nDstXSize = poWK->nDstXSize;
4559 7 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
4560 :
4561 : /* -------------------------------------------------------------------- */
4562 : /* Allocate x,y,z coordinate arrays for transformation ... one */
4563 : /* scanlines worth of positions. */
4564 : /* -------------------------------------------------------------------- */
4565 : double *padfX, *padfY, *padfZ;
4566 : int *pabSuccess;
4567 :
4568 7 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4569 7 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4570 7 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4571 7 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
4572 :
4573 : /* ==================================================================== */
4574 : /* Loop over output lines. */
4575 : /* ==================================================================== */
4576 2297 : for( iDstY = iYMin; iDstY < iYMax; iDstY++ )
4577 : {
4578 : int iDstX;
4579 :
4580 : /* -------------------------------------------------------------------- */
4581 : /* Setup points to transform to source image space. */
4582 : /* -------------------------------------------------------------------- */
4583 1620626 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4584 : {
4585 1618336 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
4586 1618336 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
4587 1618336 : padfZ[iDstX] = 0.0;
4588 : }
4589 :
4590 : /* -------------------------------------------------------------------- */
4591 : /* Transform the points from destination pixel/line coordinates */
4592 : /* to source pixel/line coordinates. */
4593 : /* -------------------------------------------------------------------- */
4594 : poWK->pfnTransformer( psJob->pTransformerArg, TRUE, nDstXSize,
4595 2290 : padfX, padfY, padfZ, pabSuccess );
4596 :
4597 : /* ==================================================================== */
4598 : /* Loop over pixels in output scanline. */
4599 : /* ==================================================================== */
4600 1620626 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4601 : {
4602 : int iDstOffset;
4603 :
4604 1618336 : COMPUTE_iSrcOffset(pabSuccess, iDstX, padfX, padfY, poWK, nSrcXSize, nSrcYSize);
4605 :
4606 : /* -------------------------------------------------------------------- */
4607 : /* Don't generate output pixels for which the source valid */
4608 : /* mask exists and is invalid. */
4609 : /* -------------------------------------------------------------------- */
4610 1488066 : if( poWK->panUnifiedSrcValid != NULL
4611 0 : && !(poWK->panUnifiedSrcValid[iSrcOffset>>5]
4612 : & (0x01 << (iSrcOffset & 0x1f))) )
4613 0 : continue;
4614 :
4615 : /* -------------------------------------------------------------------- */
4616 : /* Do not try to apply transparent source pixels to the destination.*/
4617 : /* -------------------------------------------------------------------- */
4618 1488066 : double dfDensity = 1.0;
4619 :
4620 1488066 : if( poWK->pafUnifiedSrcDensity != NULL )
4621 : {
4622 802 : dfDensity = poWK->pafUnifiedSrcDensity[iSrcOffset];
4623 802 : if( dfDensity < 0.00001 )
4624 401 : continue;
4625 : }
4626 :
4627 : /* ==================================================================== */
4628 : /* Loop processing each band. */
4629 : /* ==================================================================== */
4630 : int iBand;
4631 1487665 : iDstOffset = iDstX + iDstY * nDstXSize;
4632 :
4633 2976132 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
4634 : {
4635 1488467 : GInt16 iValue = 0;
4636 1488467 : double dfBandDensity = 0.0;
4637 :
4638 : /* -------------------------------------------------------------------- */
4639 : /* Collect the source value. */
4640 : /* -------------------------------------------------------------------- */
4641 1488467 : if ( GWKGetPixelShort( poWK, iBand, iSrcOffset, &dfBandDensity,
4642 : &iValue ) )
4643 : {
4644 1487752 : if( dfBandDensity < 1.0 )
4645 : {
4646 0 : if( dfBandDensity == 0.0 )
4647 : /* do nothing */;
4648 : else
4649 : {
4650 : /* let the general code take care of mixing */
4651 : GWKSetPixelValue( poWK, iBand, iDstOffset,
4652 : dfBandDensity, (double) iValue,
4653 0 : 0.0 );
4654 : }
4655 : }
4656 : else
4657 : {
4658 1487752 : ((GInt16 *)poWK->papabyDstImage[iBand])[iDstOffset] = iValue;
4659 : }
4660 : }
4661 : }
4662 :
4663 : /* -------------------------------------------------------------------- */
4664 : /* Mark this pixel valid/opaque in the output. */
4665 : /* -------------------------------------------------------------------- */
4666 1487665 : GWKOverlayDensity( poWK, iDstOffset, dfDensity );
4667 :
4668 1487665 : if( poWK->panDstValid != NULL )
4669 : {
4670 : poWK->panDstValid[iDstOffset>>5] |=
4671 0 : 0x01 << (iDstOffset & 0x1f);
4672 : }
4673 : } /* Next iDstX */
4674 :
4675 : /* -------------------------------------------------------------------- */
4676 : /* Report progress to the user, and optionally cancel out. */
4677 : /* -------------------------------------------------------------------- */
4678 2290 : if (psJob->pfnProgress(psJob))
4679 0 : break;
4680 : } /* Next iDstY */
4681 :
4682 : /* -------------------------------------------------------------------- */
4683 : /* Cleanup and return. */
4684 : /* -------------------------------------------------------------------- */
4685 7 : CPLFree( padfX );
4686 7 : CPLFree( padfY );
4687 7 : CPLFree( padfZ );
4688 7 : CPLFree( pabSuccess );
4689 7 : }
4690 :
4691 : /************************************************************************/
4692 : /* GWKNearestNoMasksFloat() */
4693 : /* */
4694 : /* Case for 32bit float input data with nearest neighbour */
4695 : /* resampling without concerning about masking. Should be as fast */
4696 : /* as possible for this particular transformation type. */
4697 : /************************************************************************/
4698 :
4699 : static void GWKNearestNoMasksFloatThread(void* pData);
4700 :
4701 8 : static CPLErr GWKNearestNoMasksFloat( GDALWarpKernel *poWK )
4702 : {
4703 8 : return GWKRun( poWK, "GWKNearestNoMasksFloat", GWKNearestNoMasksFloatThread );
4704 : }
4705 :
4706 8 : static void GWKNearestNoMasksFloatThread( void* pData )
4707 :
4708 : {
4709 8 : GWKJobStruct* psJob = (GWKJobStruct*) pData;
4710 8 : GDALWarpKernel *poWK = psJob->poWK;
4711 8 : int iYMin = psJob->iYMin;
4712 8 : int iYMax = psJob->iYMax;
4713 :
4714 : int iDstY;
4715 8 : int nDstXSize = poWK->nDstXSize;
4716 8 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
4717 :
4718 : /* -------------------------------------------------------------------- */
4719 : /* Allocate x,y,z coordinate arrays for transformation ... one */
4720 : /* scanlines worth of positions. */
4721 : /* -------------------------------------------------------------------- */
4722 : double *padfX, *padfY, *padfZ;
4723 : int *pabSuccess;
4724 :
4725 8 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4726 8 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4727 8 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4728 8 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
4729 :
4730 : /* ==================================================================== */
4731 : /* Loop over output lines. */
4732 : /* ==================================================================== */
4733 533 : for( iDstY = iYMin; iDstY < iYMax; iDstY++ )
4734 : {
4735 : int iDstX;
4736 :
4737 : /* -------------------------------------------------------------------- */
4738 : /* Setup points to transform to source image space. */
4739 : /* -------------------------------------------------------------------- */
4740 262732 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4741 : {
4742 262207 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
4743 262207 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
4744 262207 : padfZ[iDstX] = 0.0;
4745 : }
4746 :
4747 : /* -------------------------------------------------------------------- */
4748 : /* Transform the points from destination pixel/line coordinates */
4749 : /* to source pixel/line coordinates. */
4750 : /* -------------------------------------------------------------------- */
4751 : poWK->pfnTransformer( psJob->pTransformerArg, TRUE, nDstXSize,
4752 525 : padfX, padfY, padfZ, pabSuccess );
4753 :
4754 : /* ==================================================================== */
4755 : /* Loop over pixels in output scanline. */
4756 : /* ==================================================================== */
4757 262732 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4758 : {
4759 : int iDstOffset;
4760 :
4761 262207 : COMPUTE_iSrcOffset(pabSuccess, iDstX, padfX, padfY, poWK, nSrcXSize, nSrcYSize);
4762 :
4763 : /* ==================================================================== */
4764 : /* Loop processing each band. */
4765 : /* ==================================================================== */
4766 : int iBand;
4767 :
4768 262207 : iDstOffset = iDstX + iDstY * nDstXSize;
4769 :
4770 524414 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
4771 : {
4772 262207 : ((float *)poWK->papabyDstImage[iBand])[iDstOffset] =
4773 262207 : ((float *)poWK->papabySrcImage[iBand])[iSrcOffset];
4774 : }
4775 : }
4776 :
4777 : /* -------------------------------------------------------------------- */
4778 : /* Report progress to the user, and optionally cancel out. */
4779 : /* -------------------------------------------------------------------- */
4780 525 : if (psJob->pfnProgress(psJob))
4781 0 : break;
4782 : }
4783 :
4784 : /* -------------------------------------------------------------------- */
4785 : /* Cleanup and return. */
4786 : /* -------------------------------------------------------------------- */
4787 8 : CPLFree( padfX );
4788 8 : CPLFree( padfY );
4789 8 : CPLFree( padfZ );
4790 8 : CPLFree( pabSuccess );
4791 8 : }
4792 :
4793 : /************************************************************************/
4794 : /* GWKNearestFloat() */
4795 : /* */
4796 : /* Case for 32bit float input data with nearest neighbour */
4797 : /* resampling using valid flags. Should be as fast as possible */
4798 : /* for this particular transformation type. */
4799 : /************************************************************************/
4800 :
4801 : static void GWKNearestFloatThread(void* pData);
4802 :
4803 6 : static CPLErr GWKNearestFloat( GDALWarpKernel *poWK )
4804 : {
4805 6 : return GWKRun( poWK, "GWKNearestFloat", GWKNearestFloatThread );
4806 : }
4807 :
4808 6 : static void GWKNearestFloatThread( void* pData )
4809 :
4810 : {
4811 6 : GWKJobStruct* psJob = (GWKJobStruct*) pData;
4812 6 : GDALWarpKernel *poWK = psJob->poWK;
4813 6 : int iYMin = psJob->iYMin;
4814 6 : int iYMax = psJob->iYMax;
4815 :
4816 : int iDstY;
4817 6 : int nDstXSize = poWK->nDstXSize;
4818 6 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
4819 :
4820 : /* -------------------------------------------------------------------- */
4821 : /* Allocate x,y,z coordinate arrays for transformation ... one */
4822 : /* scanlines worth of positions. */
4823 : /* -------------------------------------------------------------------- */
4824 : double *padfX, *padfY, *padfZ;
4825 : int *pabSuccess;
4826 :
4827 6 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4828 6 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4829 6 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4830 6 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
4831 :
4832 : /* ==================================================================== */
4833 : /* Loop over output lines. */
4834 : /* ==================================================================== */
4835 774 : for( iDstY = iYMin; iDstY < iYMax; iDstY++ )
4836 : {
4837 : int iDstX;
4838 :
4839 : /* -------------------------------------------------------------------- */
4840 : /* Setup points to transform to source image space. */
4841 : /* -------------------------------------------------------------------- */
4842 393984 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4843 : {
4844 393216 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
4845 393216 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
4846 393216 : padfZ[iDstX] = 0.0;
4847 : }
4848 :
4849 : /* -------------------------------------------------------------------- */
4850 : /* Transform the points from destination pixel/line coordinates */
4851 : /* to source pixel/line coordinates. */
4852 : /* -------------------------------------------------------------------- */
4853 : poWK->pfnTransformer( psJob->pTransformerArg, TRUE, nDstXSize,
4854 768 : padfX, padfY, padfZ, pabSuccess );
4855 :
4856 : /* ==================================================================== */
4857 : /* Loop over pixels in output scanline. */
4858 : /* ==================================================================== */
4859 393984 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4860 : {
4861 : int iDstOffset;
4862 :
4863 393216 : COMPUTE_iSrcOffset(pabSuccess, iDstX, padfX, padfY, poWK, nSrcXSize, nSrcYSize);
4864 :
4865 : /* -------------------------------------------------------------------- */
4866 : /* Do not try to apply invalid source pixels to the dest. */
4867 : /* -------------------------------------------------------------------- */
4868 1202 : if( poWK->panUnifiedSrcValid != NULL
4869 0 : && !(poWK->panUnifiedSrcValid[iSrcOffset>>5]
4870 : & (0x01 << (iSrcOffset & 0x1f))) )
4871 0 : continue;
4872 :
4873 : /* -------------------------------------------------------------------- */
4874 : /* Do not try to apply transparent source pixels to the destination.*/
4875 : /* -------------------------------------------------------------------- */
4876 1202 : double dfDensity = 1.0;
4877 :
4878 1202 : if( poWK->pafUnifiedSrcDensity != NULL )
4879 : {
4880 802 : dfDensity = poWK->pafUnifiedSrcDensity[iSrcOffset];
4881 802 : if( dfDensity < 0.00001 )
4882 401 : continue;
4883 : }
4884 :
4885 : /* ==================================================================== */
4886 : /* Loop processing each band. */
4887 : /* ==================================================================== */
4888 : int iBand;
4889 :
4890 801 : iDstOffset = iDstX + iDstY * nDstXSize;
4891 :
4892 2404 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
4893 : {
4894 1603 : float fValue = 0;
4895 1603 : double dfBandDensity = 0.0;
4896 :
4897 : /* -------------------------------------------------------------------- */
4898 : /* Collect the source value. */
4899 : /* -------------------------------------------------------------------- */
4900 1603 : if ( GWKGetPixelFloat( poWK, iBand, iSrcOffset, &dfBandDensity,
4901 : &fValue ) )
4902 : {
4903 1563 : if( dfBandDensity < 1.0 )
4904 : {
4905 0 : if( dfBandDensity == 0.0 )
4906 : /* do nothing */;
4907 : else
4908 : {
4909 : /* let the general code take care of mixing */
4910 : GWKSetPixelValue( poWK, iBand, iDstOffset,
4911 : dfBandDensity, (double) fValue,
4912 0 : 0.0 );
4913 : }
4914 : }
4915 : else
4916 : {
4917 1563 : ((float *)poWK->papabyDstImage[iBand])[iDstOffset]
4918 1563 : = fValue;
4919 : }
4920 : }
4921 : }
4922 :
4923 : /* -------------------------------------------------------------------- */
4924 : /* Mark this pixel valid/opaque in the output. */
4925 : /* -------------------------------------------------------------------- */
4926 801 : GWKOverlayDensity( poWK, iDstOffset, dfDensity );
4927 :
4928 801 : if( poWK->panDstValid != NULL )
4929 : {
4930 : poWK->panDstValid[iDstOffset>>5] |=
4931 400 : 0x01 << (iDstOffset & 0x1f);
4932 : }
4933 : }
4934 :
4935 : /* -------------------------------------------------------------------- */
4936 : /* Report progress to the user, and optionally cancel out. */
4937 : /* -------------------------------------------------------------------- */
4938 768 : if (psJob->pfnProgress(psJob))
4939 0 : break;
4940 : }
4941 :
4942 : /* -------------------------------------------------------------------- */
4943 : /* Cleanup and return. */
4944 : /* -------------------------------------------------------------------- */
4945 6 : CPLFree( padfX );
4946 6 : CPLFree( padfY );
4947 6 : CPLFree( padfZ );
4948 6 : CPLFree( pabSuccess );
4949 6 : }
4950 :
|
