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