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