1 : /******************************************************************************
2 : * $Id: gdalwarpkernel.cpp 18012 2009-11-13 06:20:52Z warmerdam $
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 18012 2009-11-13 06:20:52Z warmerdam $");
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 330 : int GWKGetFilterRadius(GDALResampleAlg eResampleAlg)
48 : {
49 330 : 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 330 : GDALWarpKernel::GDALWarpKernel()
471 :
472 : {
473 330 : eResample = GRA_NearestNeighbour;
474 330 : eWorkingDataType = GDT_Unknown;
475 330 : nBands = 0;
476 330 : nDstXOff = 0;
477 330 : nDstYOff = 0;
478 330 : nDstXSize = 0;
479 330 : nDstYSize = 0;
480 330 : nSrcXOff = 0;
481 330 : nSrcYOff = 0;
482 330 : nSrcXSize = 0;
483 330 : nSrcYSize = 0;
484 330 : dfXScale = 1.0;
485 330 : dfYScale = 1.0;
486 330 : dfXFilter = 0.0;
487 330 : dfYFilter = 0.0;
488 330 : nXRadius = 0;
489 330 : nYRadius = 0;
490 330 : nFiltInitX = 0;
491 330 : nFiltInitY = 0;
492 330 : pafDstDensity = NULL;
493 330 : pafUnifiedSrcDensity = NULL;
494 330 : panDstValid = NULL;
495 330 : panUnifiedSrcValid = NULL;
496 330 : papabyDstImage = NULL;
497 330 : papabySrcImage = NULL;
498 330 : papanBandSrcValid = NULL;
499 330 : pfnProgress = GDALDummyProgress;
500 330 : pProgress = NULL;
501 330 : dfProgressBase = 0.0;
502 330 : dfProgressScale = 1.0;
503 330 : pfnTransformer = NULL;
504 330 : pTransformerArg = NULL;
505 330 : papszWarpOptions = NULL;
506 330 : }
507 :
508 : /************************************************************************/
509 : /* ~GDALWarpKernel() */
510 : /************************************************************************/
511 :
512 330 : GDALWarpKernel::~GDALWarpKernel()
513 :
514 : {
515 330 : }
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 330 : CPLErr GDALWarpKernel::PerformWarp()
530 :
531 : {
532 : CPLErr eErr;
533 :
534 330 : if( (eErr = Validate()) != CE_None )
535 0 : return eErr;
536 :
537 : // See #2445 and #3079
538 330 : if (nSrcXSize <= 0 || nSrcYSize <= 0)
539 : {
540 : pfnProgress( dfProgressBase + dfProgressScale,
541 0 : "", pProgress );
542 0 : return CE_None;
543 : }
544 :
545 : /* -------------------------------------------------------------------- */
546 : /* Pre-calculate resampling scales and window sizes for filtering. */
547 : /* -------------------------------------------------------------------- */
548 :
549 330 : dfXScale = (double)nDstXSize / nSrcXSize;
550 330 : dfYScale = (double)nDstYSize / nSrcYSize;
551 :
552 330 : dfXFilter = anGWKFilterRadius[eResample];
553 330 : dfYFilter = anGWKFilterRadius[eResample];
554 :
555 : nXRadius = ( dfXScale < 1.0 ) ?
556 330 : (int)ceil( dfXFilter / dfXScale ) :(int)dfXFilter;
557 : nYRadius = ( dfYScale < 1.0 ) ?
558 330 : (int)ceil( dfYFilter / dfYScale ) : (int)dfYFilter;
559 :
560 : // Filter window offset depends on the parity of the kernel radius
561 330 : nFiltInitX = ((anGWKFilterRadius[eResample] + 1) % 2) - nXRadius;
562 330 : nFiltInitY = ((anGWKFilterRadius[eResample] + 1) % 2) - nYRadius;
563 :
564 : /* -------------------------------------------------------------------- */
565 : /* Set up resampling functions. */
566 : /* -------------------------------------------------------------------- */
567 330 : if( CSLFetchBoolean( papszWarpOptions, "USE_GENERAL_CASE", FALSE ) )
568 0 : return GWKGeneralCase( this );
569 :
570 330 : if( eWorkingDataType == GDT_Byte
571 : && eResample == GRA_NearestNeighbour
572 : && papanBandSrcValid == NULL
573 : && panUnifiedSrcValid == NULL
574 : && pafUnifiedSrcDensity == NULL
575 : && panDstValid == NULL
576 : && pafDstDensity == NULL )
577 35 : return GWKNearestNoMasksByte( this );
578 :
579 295 : if( eWorkingDataType == GDT_Byte
580 : && eResample == GRA_Bilinear
581 : && papanBandSrcValid == NULL
582 : && panUnifiedSrcValid == NULL
583 : && pafUnifiedSrcDensity == NULL
584 : && panDstValid == NULL
585 : && pafDstDensity == NULL )
586 10 : return GWKBilinearNoMasksByte( this );
587 :
588 285 : 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 275 : 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 265 : if( eWorkingDataType == GDT_Byte
607 : && eResample == GRA_NearestNeighbour )
608 10 : return GWKNearestByte( this );
609 :
610 255 : 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 241 : 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 233 : 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 225 : 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 206 : if( (eWorkingDataType == GDT_Int16 || eWorkingDataType == GDT_UInt16)
647 : && eResample == GRA_NearestNeighbour )
648 7 : return GWKNearestShort( this );
649 :
650 199 : 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 191 : if( eWorkingDataType == GDT_Float32
660 : && eResample == GRA_NearestNeighbour )
661 2 : return GWKNearestFloat( this );
662 :
663 189 : 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 330 : CPLErr GDALWarpKernel::Validate()
684 :
685 : {
686 330 : 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 330 : 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 1789850 : static void GWKOverlayDensity( GDALWarpKernel *poWK, int iDstOffset,
706 : double dfDensity )
707 : {
708 1789850 : if( dfDensity < 0.0001 || poWK->pafDstDensity == NULL )
709 1784048 : return;
710 :
711 5802 : poWK->pafDstDensity[iDstOffset] = (float)
712 5802 : ( 1.0 - (1.0 - dfDensity) * (1.0 - poWK->pafDstDensity[iDstOffset]) );
713 : }
714 :
715 : /************************************************************************/
716 : /* GWKSetPixelValue() */
717 : /************************************************************************/
718 :
719 282451 : static int GWKSetPixelValue( GDALWarpKernel *poWK, int iBand,
720 : int iDstOffset, double dfDensity,
721 : double dfReal, double dfImag )
722 :
723 : {
724 282451 : 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 282451 : 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 0 : && !((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 : 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 282451 : switch( poWK->eWorkingDataType )
847 : {
848 : case GDT_Byte:
849 279679 : CLAMP(GByte, 0.0, 255.0);
850 279679 : 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 282451 : return TRUE;
921 : }
922 :
923 : /************************************************************************/
924 : /* GWKGetPixelValue() */
925 : /************************************************************************/
926 :
927 479 : static int GWKGetPixelValue( GDALWarpKernel *poWK, int iBand,
928 : int iSrcOffset, double *pdfDensity,
929 : double *pdfReal, double *pdfImag )
930 :
931 : {
932 479 : GByte *pabySrc = poWK->papabySrcImage[iBand];
933 :
934 479 : if( poWK->panUnifiedSrcValid != NULL
935 0 : && !((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 0 : && poWK->papanBandSrcValid[iBand] != NULL
944 0 : && !((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 1925779 : static int GWKGetPixelRow( GDALWarpKernel *poWK, int iBand,
1026 : int iSrcOffset, int nHalfSrcLen,
1027 : double adfDensity[],
1028 : double adfReal[], double adfImag[] )
1029 : {
1030 : // We know that nSrcLen is even, so we can *always* unroll loops 2x
1031 1925779 : int nSrcLen = nHalfSrcLen * 2;
1032 1925779 : int bHasValid = FALSE;
1033 : int i;
1034 :
1035 : // Init the density
1036 9617385 : for ( i = 0; i < nSrcLen; i += 2 )
1037 : {
1038 7691606 : adfDensity[i] = 1.0;
1039 7691606 : adfDensity[i+1] = 1.0;
1040 : }
1041 :
1042 1925779 : 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 1925779 : if ( poWK->papanBandSrcValid != NULL
1067 0 : && poWK->papanBandSrcValid[iBand] != NULL)
1068 : {
1069 0 : for ( i = 0; i < nSrcLen; i += 2 )
1070 : {
1071 0 : if(poWK->papanBandSrcValid[iBand][(iSrcOffset+i)>>5]
1072 : & (0x01 << ((iSrcOffset+i) & 0x1f)))
1073 0 : bHasValid = TRUE;
1074 : else
1075 0 : adfDensity[i] = 0.0;
1076 :
1077 0 : if(poWK->papanBandSrcValid[iBand][(iSrcOffset+i+1)>>5]
1078 : & (0x01 << ((iSrcOffset+i+1) & 0x1f)))
1079 0 : bHasValid = TRUE;
1080 : else
1081 0 : adfDensity[i+1] = 0.0;
1082 : }
1083 :
1084 : // Reset or fail as needed
1085 0 : if ( bHasValid )
1086 0 : bHasValid = FALSE;
1087 : else
1088 0 : return FALSE;
1089 : }
1090 :
1091 : // Fetch data
1092 1925779 : switch( poWK->eWorkingDataType )
1093 : {
1094 : case GDT_Byte:
1095 : {
1096 1919280 : GByte* pSrc = (GByte*) poWK->papabySrcImage[iBand];
1097 1919280 : pSrc += iSrcOffset;
1098 9594840 : for ( i = 0; i < nSrcLen; i += 2 )
1099 : {
1100 7675560 : adfReal[i] = pSrc[i];
1101 7675560 : adfReal[i+1] = pSrc[i+1];
1102 : }
1103 1919280 : memset( adfImag, 0, nSrcLen * sizeof(double) );
1104 1919280 : break;
1105 : }
1106 :
1107 : case GDT_Int16:
1108 : {
1109 280 : GInt16* pSrc = (GInt16*) poWK->papabySrcImage[iBand];
1110 280 : pSrc += iSrcOffset;
1111 1314 : for ( i = 0; i < nSrcLen; i += 2 )
1112 : {
1113 1034 : adfReal[i] = pSrc[i];
1114 1034 : adfReal[i+1] = pSrc[i+1];
1115 : }
1116 280 : memset( adfImag, 0, nSrcLen * sizeof(double) );
1117 280 : break;
1118 : }
1119 :
1120 : case GDT_UInt16:
1121 : {
1122 691 : GUInt16* pSrc = (GUInt16*) poWK->papabySrcImage[iBand];
1123 691 : pSrc += iSrcOffset;
1124 2359 : for ( i = 0; i < nSrcLen; i += 2 )
1125 : {
1126 1668 : adfReal[i] = pSrc[i];
1127 1668 : adfReal[i+1] = pSrc[i+1];
1128 : }
1129 691 : memset( adfImag, 0, nSrcLen * sizeof(double) );
1130 691 : break;
1131 : }
1132 :
1133 : case GDT_Int32:
1134 : {
1135 691 : GInt32* pSrc = (GInt32*) poWK->papabySrcImage[iBand];
1136 691 : pSrc += iSrcOffset;
1137 2359 : for ( i = 0; i < nSrcLen; i += 2 )
1138 : {
1139 1668 : adfReal[i] = pSrc[i];
1140 1668 : adfReal[i+1] = pSrc[i+1];
1141 : }
1142 691 : memset( adfImag, 0, nSrcLen * sizeof(double) );
1143 691 : break;
1144 : }
1145 :
1146 : case GDT_UInt32:
1147 : {
1148 691 : GUInt32* pSrc = (GUInt32*) poWK->papabySrcImage[iBand];
1149 691 : pSrc += iSrcOffset;
1150 2359 : for ( i = 0; i < nSrcLen; i += 2 )
1151 : {
1152 1668 : adfReal[i] = pSrc[i];
1153 1668 : adfReal[i+1] = pSrc[i+1];
1154 : }
1155 691 : memset( adfImag, 0, nSrcLen * sizeof(double) );
1156 691 : break;
1157 : }
1158 :
1159 : case GDT_Float32:
1160 : {
1161 691 : float* pSrc = (float*) poWK->papabySrcImage[iBand];
1162 691 : pSrc += iSrcOffset;
1163 2359 : for ( i = 0; i < nSrcLen; i += 2 )
1164 : {
1165 1668 : adfReal[i] = pSrc[i];
1166 1668 : adfReal[i+1] = pSrc[i+1];
1167 : }
1168 691 : memset( adfImag, 0, nSrcLen * sizeof(double) );
1169 691 : break;
1170 : }
1171 :
1172 : case GDT_Float64:
1173 : {
1174 691 : double* pSrc = (double*) poWK->papabySrcImage[iBand];
1175 691 : pSrc += iSrcOffset;
1176 2359 : for ( i = 0; i < nSrcLen; i += 2 )
1177 : {
1178 1668 : adfReal[i] = pSrc[i];
1179 1668 : adfReal[i+1] = pSrc[i+1];
1180 : }
1181 691 : memset( adfImag, 0, nSrcLen * sizeof(double) );
1182 691 : break;
1183 : }
1184 :
1185 : case GDT_CInt16:
1186 : {
1187 691 : GInt16* pSrc = (GInt16*) poWK->papabySrcImage[iBand];
1188 691 : pSrc += 2 * iSrcOffset;
1189 2359 : for ( i = 0; i < nSrcLen; i += 2 )
1190 : {
1191 1668 : adfReal[i] = pSrc[2*i];
1192 1668 : adfImag[i] = pSrc[2*i+1];
1193 :
1194 1668 : adfReal[i+1] = pSrc[2*i+2];
1195 1668 : adfImag[i+1] = pSrc[2*i+3];
1196 : }
1197 691 : break;
1198 : }
1199 :
1200 : case GDT_CInt32:
1201 : {
1202 691 : GInt32* pSrc = (GInt32*) poWK->papabySrcImage[iBand];
1203 691 : pSrc += 2 * iSrcOffset;
1204 2359 : for ( i = 0; i < nSrcLen; i += 2 )
1205 : {
1206 1668 : adfReal[i] = pSrc[2*i];
1207 1668 : adfImag[i] = pSrc[2*i+1];
1208 :
1209 1668 : adfReal[i+1] = pSrc[2*i+2];
1210 1668 : adfImag[i+1] = pSrc[2*i+3];
1211 : }
1212 691 : break;
1213 : }
1214 :
1215 : case GDT_CFloat32:
1216 : {
1217 691 : float* pSrc = (float*) poWK->papabySrcImage[iBand];
1218 691 : pSrc += 2 * iSrcOffset;
1219 2359 : for ( i = 0; i < nSrcLen; i += 2 )
1220 : {
1221 1668 : adfReal[i] = pSrc[2*i];
1222 1668 : adfImag[i] = pSrc[2*i+1];
1223 :
1224 1668 : adfReal[i+1] = pSrc[2*i+2];
1225 1668 : adfImag[i+1] = pSrc[2*i+3];
1226 : }
1227 691 : break;
1228 : }
1229 :
1230 :
1231 : case GDT_CFloat64:
1232 : {
1233 691 : double* pSrc = (double*) poWK->papabySrcImage[iBand];
1234 691 : pSrc += 2 * iSrcOffset;
1235 2359 : for ( i = 0; i < nSrcLen; i += 2 )
1236 : {
1237 1668 : adfReal[i] = pSrc[2*i];
1238 1668 : adfImag[i] = pSrc[2*i+1];
1239 :
1240 1668 : adfReal[i+1] = pSrc[2*i+2];
1241 1668 : adfImag[i+1] = pSrc[2*i+3];
1242 : }
1243 691 : break;
1244 : }
1245 :
1246 : default:
1247 : CPLAssert(FALSE);
1248 0 : memset( adfDensity, 0, nSrcLen * sizeof(double) );
1249 0 : return FALSE;
1250 : }
1251 :
1252 1925779 : if( poWK->pafUnifiedSrcDensity == NULL )
1253 : {
1254 9617385 : for ( i = 0; i < nSrcLen; i += 2 )
1255 : {
1256 : // Take into account earlier calcs
1257 7691606 : if(adfDensity[i] > 0.000000001)
1258 : {
1259 7691606 : adfDensity[i] = 1.0;
1260 7691606 : bHasValid = TRUE;
1261 : }
1262 :
1263 7691606 : if(adfDensity[i+1] > 0.000000001)
1264 : {
1265 7691606 : adfDensity[i+1] = 1.0;
1266 7691606 : 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 1925779 : return bHasValid;
1287 : }
1288 :
1289 : /************************************************************************/
1290 : /* GWKGetPixelByte() */
1291 : /************************************************************************/
1292 :
1293 30805 : static int GWKGetPixelByte( GDALWarpKernel *poWK, int iBand,
1294 : int iSrcOffset, double *pdfDensity,
1295 : GByte *pbValue )
1296 :
1297 : {
1298 30805 : GByte *pabySrc = poWK->papabySrcImage[iBand];
1299 :
1300 31106 : if ( ( poWK->panUnifiedSrcValid != NULL
1301 301 : && !((poWK->panUnifiedSrcValid[iSrcOffset>>5]
1302 : & (0x01 << (iSrcOffset & 0x1f))) ) )
1303 : || ( poWK->papanBandSrcValid != NULL
1304 0 : && poWK->papanBandSrcValid[iBand] != NULL
1305 0 : && !((poWK->papanBandSrcValid[iBand][iSrcOffset>>5]
1306 : & (0x01 << (iSrcOffset & 0x1f)))) ) )
1307 : {
1308 0 : *pdfDensity = 0.0;
1309 0 : return FALSE;
1310 : }
1311 :
1312 30805 : *pbValue = pabySrc[iSrcOffset];
1313 :
1314 30805 : if ( poWK->pafUnifiedSrcDensity == NULL )
1315 15702 : *pdfDensity = 1.0;
1316 : else
1317 15103 : *pdfDensity = poWK->pafUnifiedSrcDensity[iSrcOffset];
1318 :
1319 30805 : return *pdfDensity != 0.0;
1320 : }
1321 :
1322 : /************************************************************************/
1323 : /* GWKGetPixelShort() */
1324 : /************************************************************************/
1325 :
1326 1488467 : static int GWKGetPixelShort( GDALWarpKernel *poWK, int iBand,
1327 : int iSrcOffset, double *pdfDensity,
1328 : GInt16 *piValue )
1329 :
1330 : {
1331 1488467 : GInt16 *pabySrc = (GInt16 *)poWK->papabySrcImage[iBand];
1332 :
1333 4462995 : if ( ( poWK->panUnifiedSrcValid != NULL
1334 0 : && !((poWK->panUnifiedSrcValid[iSrcOffset>>5]
1335 : & (0x01 << (iSrcOffset & 0x1f))) ) )
1336 : || ( poWK->papanBandSrcValid != NULL
1337 1487264 : && poWK->papanBandSrcValid[iBand] != NULL
1338 1487264 : && !((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 1203 : static int GWKGetPixelFloat( GDALWarpKernel *poWK, int iBand,
1360 : int iSrcOffset, double *pdfDensity,
1361 : float *pfValue )
1362 :
1363 : {
1364 1203 : float *pabySrc = (float *)poWK->papabySrcImage[iBand];
1365 :
1366 1203 : if ( ( poWK->panUnifiedSrcValid != NULL
1367 0 : && !((poWK->panUnifiedSrcValid[iSrcOffset>>5]
1368 : & (0x01 << (iSrcOffset & 0x1f))) ) )
1369 : || ( poWK->papanBandSrcValid != NULL
1370 0 : && poWK->papanBandSrcValid[iBand] != NULL
1371 0 : && !((poWK->papanBandSrcValid[iBand][iSrcOffset>>5]
1372 : & (0x01 << (iSrcOffset & 0x1f)))) ) )
1373 : {
1374 0 : *pdfDensity = 0.0;
1375 0 : return FALSE;
1376 : }
1377 :
1378 1203 : *pfValue = pabySrc[iSrcOffset];
1379 :
1380 1203 : if ( poWK->pafUnifiedSrcDensity == NULL )
1381 0 : *pdfDensity = 1.0;
1382 : else
1383 1203 : *pdfDensity = poWK->pafUnifiedSrcDensity[iSrcOffset];
1384 :
1385 1203 : return *pdfDensity != 0.0;
1386 : }
1387 :
1388 : /************************************************************************/
1389 : /* GWKBilinearResample() */
1390 : /* Set of bilinear interpolators */
1391 : /************************************************************************/
1392 :
1393 1008 : static int GWKBilinearResample( GDALWarpKernel *poWK, int iBand,
1394 : double dfSrcX, double dfSrcY,
1395 : double *pdfDensity,
1396 : double *pdfReal, double *pdfImag )
1397 :
1398 : {
1399 : // Save as local variables to avoid following pointers
1400 1008 : int nSrcXSize = poWK->nSrcXSize;
1401 1008 : int nSrcYSize = poWK->nSrcYSize;
1402 :
1403 1008 : int iSrcX = (int) floor(dfSrcX - 0.5);
1404 1008 : int iSrcY = (int) floor(dfSrcY - 0.5);
1405 1008 : int iSrcOffset = iSrcX + iSrcY * nSrcXSize;
1406 1008 : double dfRatioX = 1.5 - (dfSrcX - iSrcX);
1407 1008 : double dfRatioY = 1.5 - (dfSrcY - iSrcY);
1408 : double adfDensity[2], adfReal[2], adfImag[2];
1409 1008 : double dfAccumulatorReal = 0.0, dfAccumulatorImag = 0.0;
1410 1008 : double dfAccumulatorDensity = 0.0;
1411 1008 : double dfAccumulatorDivisor = 0.0;
1412 1008 : int bShifted = FALSE;
1413 :
1414 : // Shift so we don't overrun the array
1415 1008 : if( nSrcXSize * nSrcYSize == iSrcOffset + 1
1416 : || nSrcXSize * nSrcYSize == iSrcOffset + nSrcXSize + 1 )
1417 : {
1418 72 : bShifted = TRUE;
1419 72 : --iSrcOffset;
1420 : }
1421 :
1422 : // Get pixel row
1423 1008 : if ( iSrcY >= 0 && iSrcY < nSrcYSize
1424 : && iSrcOffset >= 0 && iSrcOffset < nSrcXSize * nSrcYSize
1425 : && GWKGetPixelRow( poWK, iBand, iSrcOffset, 1,
1426 : adfDensity, adfReal, adfImag ) )
1427 : {
1428 792 : double dfMult1 = dfRatioX * dfRatioY;
1429 792 : double dfMult2 = (1.0-dfRatioX) * dfRatioY;
1430 :
1431 : // Shifting corrected
1432 792 : if ( bShifted )
1433 : {
1434 72 : adfReal[0] = adfReal[1];
1435 72 : adfImag[0] = adfImag[1];
1436 72 : adfDensity[0] = adfDensity[1];
1437 : }
1438 :
1439 : // Upper Left Pixel
1440 1584 : if ( iSrcX >= 0 && iSrcX < nSrcXSize
1441 792 : && adfDensity[0] > 0.000000001 )
1442 : {
1443 792 : dfAccumulatorDivisor += dfMult1;
1444 :
1445 792 : dfAccumulatorReal += adfReal[0] * dfMult1;
1446 792 : dfAccumulatorImag += adfImag[0] * dfMult1;
1447 792 : dfAccumulatorDensity += adfDensity[0] * dfMult1;
1448 : }
1449 :
1450 : // Upper Right Pixel
1451 1422 : if ( iSrcX+1 >= 0 && iSrcX+1 < nSrcXSize
1452 630 : && adfDensity[1] > 0.000000001 )
1453 : {
1454 630 : dfAccumulatorDivisor += dfMult2;
1455 :
1456 630 : dfAccumulatorReal += adfReal[1] * dfMult2;
1457 630 : dfAccumulatorImag += adfImag[1] * dfMult2;
1458 630 : dfAccumulatorDensity += adfDensity[1] * dfMult2;
1459 : }
1460 : }
1461 :
1462 : // Get pixel row
1463 1008 : if ( iSrcY+1 >= 0 && iSrcY+1 < nSrcYSize
1464 : && iSrcOffset+nSrcXSize >= 0
1465 : && iSrcOffset+nSrcXSize < nSrcXSize * nSrcYSize
1466 : && GWKGetPixelRow( poWK, iBand, iSrcOffset+nSrcXSize, 1,
1467 : adfDensity, adfReal, adfImag ) )
1468 : {
1469 828 : double dfMult1 = dfRatioX * (1.0-dfRatioY);
1470 828 : double dfMult2 = (1.0-dfRatioX) * (1.0-dfRatioY);
1471 :
1472 : // Shifting corrected
1473 828 : if ( bShifted )
1474 : {
1475 36 : adfReal[0] = adfReal[1];
1476 36 : adfImag[0] = adfImag[1];
1477 36 : adfDensity[0] = adfDensity[1];
1478 : }
1479 :
1480 : // Lower Left Pixel
1481 1656 : if ( iSrcX >= 0 && iSrcX < nSrcXSize
1482 828 : && adfDensity[0] > 0.000000001 )
1483 : {
1484 828 : dfAccumulatorDivisor += dfMult1;
1485 :
1486 828 : dfAccumulatorReal += adfReal[0] * dfMult1;
1487 828 : dfAccumulatorImag += adfImag[0] * dfMult1;
1488 828 : dfAccumulatorDensity += adfDensity[0] * dfMult1;
1489 : }
1490 :
1491 : // Lower Right Pixel
1492 1476 : if ( iSrcX+1 >= 0 && iSrcX+1 < nSrcXSize
1493 648 : && adfDensity[1] > 0.000000001 )
1494 : {
1495 648 : dfAccumulatorDivisor += dfMult2;
1496 :
1497 648 : dfAccumulatorReal += adfReal[1] * dfMult2;
1498 648 : dfAccumulatorImag += adfImag[1] * dfMult2;
1499 648 : dfAccumulatorDensity += adfDensity[1] * dfMult2;
1500 : }
1501 : }
1502 :
1503 : /* -------------------------------------------------------------------- */
1504 : /* Return result. */
1505 : /* -------------------------------------------------------------------- */
1506 1008 : if ( dfAccumulatorDivisor == 1.0 )
1507 : {
1508 594 : *pdfReal = dfAccumulatorReal;
1509 594 : *pdfImag = dfAccumulatorImag;
1510 594 : *pdfDensity = dfAccumulatorDensity;
1511 594 : return TRUE;
1512 : }
1513 414 : else if ( dfAccumulatorDivisor < 0.00001 )
1514 : {
1515 0 : *pdfReal = 0.0;
1516 0 : *pdfImag = 0.0;
1517 0 : *pdfDensity = 0.0;
1518 0 : return FALSE;
1519 : }
1520 : else
1521 : {
1522 414 : *pdfReal = dfAccumulatorReal / dfAccumulatorDivisor;
1523 414 : *pdfImag = dfAccumulatorImag / dfAccumulatorDivisor;
1524 414 : *pdfDensity = dfAccumulatorDensity / dfAccumulatorDivisor;
1525 414 : return TRUE;
1526 : }
1527 : }
1528 :
1529 288782 : static int GWKBilinearResampleNoMasksByte( GDALWarpKernel *poWK, int iBand,
1530 : double dfSrcX, double dfSrcY,
1531 : GByte *pbValue )
1532 :
1533 : {
1534 288782 : double dfAccumulator = 0.0;
1535 288782 : double dfAccumulatorDivisor = 0.0;
1536 :
1537 288782 : int iSrcX = (int) floor(dfSrcX - 0.5);
1538 288782 : int iSrcY = (int) floor(dfSrcY - 0.5);
1539 288782 : int iSrcOffset = iSrcX + iSrcY * poWK->nSrcXSize;
1540 288782 : double dfRatioX = 1.5 - (dfSrcX - iSrcX);
1541 288782 : double dfRatioY = 1.5 - (dfSrcY - iSrcY);
1542 :
1543 : // Upper Left Pixel
1544 288782 : if( iSrcX >= 0 && iSrcX < poWK->nSrcXSize
1545 : && iSrcY >= 0 && iSrcY < poWK->nSrcYSize )
1546 : {
1547 281993 : double dfMult = dfRatioX * dfRatioY;
1548 :
1549 281993 : dfAccumulatorDivisor += dfMult;
1550 :
1551 : dfAccumulator +=
1552 281993 : (double)poWK->papabySrcImage[iBand][iSrcOffset] * dfMult;
1553 : }
1554 :
1555 : // Upper Right Pixel
1556 288782 : if( iSrcX+1 >= 0 && iSrcX+1 < poWK->nSrcXSize
1557 : && iSrcY >= 0 && iSrcY < poWK->nSrcYSize )
1558 : {
1559 285029 : double dfMult = (1.0-dfRatioX) * dfRatioY;
1560 :
1561 285029 : dfAccumulatorDivisor += dfMult;
1562 :
1563 : dfAccumulator +=
1564 285029 : (double)poWK->papabySrcImage[iBand][iSrcOffset+1] * dfMult;
1565 : }
1566 :
1567 : // Lower Right Pixel
1568 288782 : if( iSrcX+1 >= 0 && iSrcX+1 < poWK->nSrcXSize
1569 : && iSrcY+1 >= 0 && iSrcY+1 < poWK->nSrcYSize )
1570 : {
1571 288103 : double dfMult = (1.0-dfRatioX) * (1.0-dfRatioY);
1572 :
1573 288103 : dfAccumulatorDivisor += dfMult;
1574 :
1575 : dfAccumulator +=
1576 288103 : (double)poWK->papabySrcImage[iBand][iSrcOffset+1+poWK->nSrcXSize]
1577 288103 : * dfMult;
1578 : }
1579 :
1580 : // Lower Left Pixel
1581 288782 : if( iSrcX >= 0 && iSrcX < poWK->nSrcXSize
1582 : && iSrcY+1 >= 0 && iSrcY+1 < poWK->nSrcYSize )
1583 : {
1584 285054 : double dfMult = dfRatioX * (1.0-dfRatioY);
1585 :
1586 285054 : dfAccumulatorDivisor += dfMult;
1587 :
1588 : dfAccumulator +=
1589 285054 : (double)poWK->papabySrcImage[iBand][iSrcOffset+poWK->nSrcXSize]
1590 285054 : * dfMult;
1591 : }
1592 :
1593 : /* -------------------------------------------------------------------- */
1594 : /* Return result. */
1595 : /* -------------------------------------------------------------------- */
1596 : double dfValue;
1597 :
1598 288782 : if( dfAccumulatorDivisor < 0.00001 )
1599 : {
1600 0 : *pbValue = 0;
1601 0 : return FALSE;
1602 : }
1603 288782 : else if( dfAccumulatorDivisor == 1.0 )
1604 : {
1605 276809 : dfValue = dfAccumulator;
1606 : }
1607 : else
1608 : {
1609 11973 : dfValue = dfAccumulator / dfAccumulatorDivisor;
1610 : }
1611 :
1612 288782 : if ( dfValue < 0.0 )
1613 0 : *pbValue = 0;
1614 288782 : else if ( dfValue > 255.0 )
1615 925 : *pbValue = 255;
1616 : else
1617 287857 : *pbValue = (GByte)(0.5 + dfValue);
1618 :
1619 288782 : return TRUE;
1620 : }
1621 :
1622 272494 : static int GWKBilinearResampleNoMasksShort( GDALWarpKernel *poWK, int iBand,
1623 : double dfSrcX, double dfSrcY,
1624 : GInt16 *piValue )
1625 :
1626 : {
1627 272494 : double dfAccumulator = 0.0;
1628 272494 : double dfAccumulatorDivisor = 0.0;
1629 :
1630 272494 : int iSrcX = (int) floor(dfSrcX - 0.5);
1631 272494 : int iSrcY = (int) floor(dfSrcY - 0.5);
1632 272494 : int iSrcOffset = iSrcX + iSrcY * poWK->nSrcXSize;
1633 272494 : double dfRatioX = 1.5 - (dfSrcX - iSrcX);
1634 272494 : double dfRatioY = 1.5 - (dfSrcY - iSrcY);
1635 :
1636 : // Upper Left Pixel
1637 272494 : if( iSrcX >= 0 && iSrcX < poWK->nSrcXSize
1638 : && iSrcY >= 0 && iSrcY < poWK->nSrcYSize )
1639 : {
1640 266132 : double dfMult = dfRatioX * dfRatioY;
1641 :
1642 266132 : dfAccumulatorDivisor += dfMult;
1643 :
1644 : dfAccumulator +=
1645 266132 : (double)((GInt16 *)poWK->papabySrcImage[iBand])[iSrcOffset]
1646 266132 : * dfMult;
1647 : }
1648 :
1649 : // Upper Right Pixel
1650 272494 : if( iSrcX+1 >= 0 && iSrcX+1 < poWK->nSrcXSize
1651 : && iSrcY >= 0 && iSrcY < poWK->nSrcYSize )
1652 : {
1653 269168 : double dfMult = (1.0-dfRatioX) * dfRatioY;
1654 :
1655 269168 : dfAccumulatorDivisor += dfMult;
1656 :
1657 : dfAccumulator +=
1658 269168 : (double)((GInt16 *)poWK->papabySrcImage[iBand])[iSrcOffset+1] * dfMult;
1659 : }
1660 :
1661 : // Lower Right Pixel
1662 272494 : if( iSrcX+1 >= 0 && iSrcX+1 < poWK->nSrcXSize
1663 : && iSrcY+1 >= 0 && iSrcY+1 < poWK->nSrcYSize )
1664 : {
1665 272341 : double dfMult = (1.0-dfRatioX) * (1.0-dfRatioY);
1666 :
1667 272341 : dfAccumulatorDivisor += dfMult;
1668 :
1669 : dfAccumulator +=
1670 272341 : (double)((GInt16 *)poWK->papabySrcImage[iBand])[iSrcOffset+1+poWK->nSrcXSize]
1671 272341 : * dfMult;
1672 : }
1673 :
1674 : // Lower Left Pixel
1675 272494 : if( iSrcX >= 0 && iSrcX < poWK->nSrcXSize
1676 : && iSrcY+1 >= 0 && iSrcY+1 < poWK->nSrcYSize )
1677 : {
1678 269292 : double dfMult = dfRatioX * (1.0-dfRatioY);
1679 :
1680 269292 : dfAccumulatorDivisor += dfMult;
1681 :
1682 : dfAccumulator +=
1683 269292 : (double)((GInt16 *)poWK->papabySrcImage[iBand])[iSrcOffset+poWK->nSrcXSize]
1684 269292 : * dfMult;
1685 : }
1686 :
1687 : /* -------------------------------------------------------------------- */
1688 : /* Return result. */
1689 : /* -------------------------------------------------------------------- */
1690 272494 : if( dfAccumulatorDivisor == 1.0 )
1691 : {
1692 261669 : *piValue = (GInt16)(0.5 + dfAccumulator);
1693 261669 : return TRUE;
1694 : }
1695 10825 : else if( dfAccumulatorDivisor < 0.00001 )
1696 : {
1697 0 : *piValue = 0;
1698 0 : return FALSE;
1699 : }
1700 : else
1701 : {
1702 10825 : *piValue = (GInt16)(0.5 + dfAccumulator / dfAccumulatorDivisor);
1703 10825 : return TRUE;
1704 : }
1705 : }
1706 :
1707 : /************************************************************************/
1708 : /* GWKCubicResample() */
1709 : /* Set of bicubic interpolators using cubic convolution. */
1710 : /************************************************************************/
1711 :
1712 : #define CubicConvolution(distance1,distance2,distance3,f0,f1,f2,f3) \
1713 : ( ( -f0 + f1 - f2 + f3) * distance3 \
1714 : + (2.0*(f0 - f1) + f2 - f3) * distance2 \
1715 : + ( -f0 + f2 ) * distance1 \
1716 : + f1 )
1717 :
1718 558 : static int GWKCubicResample( GDALWarpKernel *poWK, int iBand,
1719 : double dfSrcX, double dfSrcY,
1720 : double *pdfDensity,
1721 : double *pdfReal, double *pdfImag )
1722 :
1723 : {
1724 558 : int iSrcX = (int) floor( dfSrcX - 0.5 );
1725 558 : int iSrcY = (int) floor( dfSrcY - 0.5 );
1726 558 : int iSrcOffset = iSrcX + iSrcY * poWK->nSrcXSize;
1727 558 : double dfDeltaX = dfSrcX - 0.5 - iSrcX;
1728 558 : double dfDeltaY = dfSrcY - 0.5 - iSrcY;
1729 558 : double dfDeltaX2 = dfDeltaX * dfDeltaX;
1730 558 : double dfDeltaY2 = dfDeltaY * dfDeltaY;
1731 558 : double dfDeltaX3 = dfDeltaX2 * dfDeltaX;
1732 558 : double dfDeltaY3 = dfDeltaY2 * dfDeltaY;
1733 : double adfValueDens[4], adfValueReal[4], adfValueImag[4];
1734 : double adfDensity[4], adfReal[4], adfImag[4];
1735 : int i;
1736 :
1737 : // Get the bilinear interpolation at the image borders
1738 558 : if ( iSrcX - 1 < 0 || iSrcX + 2 >= poWK->nSrcXSize
1739 : || iSrcY - 1 < 0 || iSrcY + 2 >= poWK->nSrcYSize )
1740 : return GWKBilinearResample( poWK, iBand, dfSrcX, dfSrcY,
1741 450 : pdfDensity, pdfReal, pdfImag );
1742 :
1743 540 : for ( i = -1; i < 3; i++ )
1744 : {
1745 2160 : if ( !GWKGetPixelRow(poWK, iBand, iSrcOffset + i * poWK->nSrcXSize - 1,
1746 : 2, adfDensity, adfReal, adfImag)
1747 432 : || adfDensity[0] < 0.000000001
1748 432 : || adfDensity[1] < 0.000000001
1749 432 : || adfDensity[2] < 0.000000001
1750 432 : || adfDensity[3] < 0.000000001 )
1751 : {
1752 : return GWKBilinearResample( poWK, iBand, dfSrcX, dfSrcY,
1753 0 : pdfDensity, pdfReal, pdfImag );
1754 : }
1755 :
1756 3456 : adfValueDens[i + 1] = CubicConvolution(dfDeltaX, dfDeltaX2, dfDeltaX3,
1757 3456 : adfDensity[0], adfDensity[1], adfDensity[2], adfDensity[3]);
1758 3456 : adfValueReal[i + 1] = CubicConvolution(dfDeltaX, dfDeltaX2, dfDeltaX3,
1759 3456 : adfReal[0], adfReal[1], adfReal[2], adfReal[3]);
1760 3456 : adfValueImag[i + 1] = CubicConvolution(dfDeltaX, dfDeltaX2, dfDeltaX3,
1761 3456 : adfImag[0], adfImag[1], adfImag[2], adfImag[3]);
1762 : }
1763 :
1764 :
1765 : /* -------------------------------------------------------------------- */
1766 : /* For now, if we have any pixels missing in the kernel area, */
1767 : /* we fallback on using bilinear interpolation. Ideally we */
1768 : /* should do "weight adjustment" of our results similarly to */
1769 : /* what is done for the cubic spline and lanc. interpolators. */
1770 : /* -------------------------------------------------------------------- */
1771 864 : *pdfDensity = CubicConvolution(dfDeltaY, dfDeltaY2, dfDeltaY3,
1772 : adfValueDens[0], adfValueDens[1],
1773 864 : adfValueDens[2], adfValueDens[3]);
1774 864 : *pdfReal = CubicConvolution(dfDeltaY, dfDeltaY2, dfDeltaY3,
1775 : adfValueReal[0], adfValueReal[1],
1776 864 : adfValueReal[2], adfValueReal[3]);
1777 864 : *pdfImag = CubicConvolution(dfDeltaY, dfDeltaY2, dfDeltaY3,
1778 : adfValueImag[0], adfValueImag[1],
1779 864 : adfValueImag[2], adfValueImag[3]);
1780 :
1781 108 : return TRUE;
1782 : }
1783 :
1784 278207 : static int GWKCubicResampleNoMasksByte( GDALWarpKernel *poWK, int iBand,
1785 : double dfSrcX, double dfSrcY,
1786 : GByte *pbValue )
1787 :
1788 : {
1789 278207 : int iSrcX = (int) floor( dfSrcX - 0.5 );
1790 278207 : int iSrcY = (int) floor( dfSrcY - 0.5 );
1791 278207 : int iSrcOffset = iSrcX + iSrcY * poWK->nSrcXSize;
1792 278207 : double dfDeltaX = dfSrcX - 0.5 - iSrcX;
1793 278207 : double dfDeltaY = dfSrcY - 0.5 - iSrcY;
1794 278207 : double dfDeltaX2 = dfDeltaX * dfDeltaX;
1795 278207 : double dfDeltaY2 = dfDeltaY * dfDeltaY;
1796 278207 : double dfDeltaX3 = dfDeltaX2 * dfDeltaX;
1797 278207 : double dfDeltaY3 = dfDeltaY2 * dfDeltaY;
1798 : double adfValue[4];
1799 : int i;
1800 :
1801 : // Get the bilinear interpolation at the image borders
1802 278207 : if ( iSrcX - 1 < 0 || iSrcX + 2 >= poWK->nSrcXSize
1803 : || iSrcY - 1 < 0 || iSrcY + 2 >= poWK->nSrcYSize )
1804 : return GWKBilinearResampleNoMasksByte( poWK, iBand, dfSrcX, dfSrcY,
1805 10574 : pbValue);
1806 :
1807 1338165 : for ( i = -1; i < 3; i++ )
1808 : {
1809 1070532 : int iOffset = iSrcOffset + i * poWK->nSrcXSize;
1810 :
1811 11775852 : adfValue[i + 1] = CubicConvolution(dfDeltaX, dfDeltaX2, dfDeltaX3,
1812 : (double)poWK->papabySrcImage[iBand][iOffset - 1],
1813 : (double)poWK->papabySrcImage[iBand][iOffset],
1814 : (double)poWK->papabySrcImage[iBand][iOffset + 1],
1815 11775852 : (double)poWK->papabySrcImage[iBand][iOffset + 2]);
1816 : }
1817 :
1818 267633 : double dfValue = CubicConvolution(dfDeltaY, dfDeltaY2, dfDeltaY3,
1819 : adfValue[0], adfValue[1], adfValue[2], adfValue[3]);
1820 :
1821 267633 : if ( dfValue < 0.0 )
1822 8 : *pbValue = 0;
1823 267625 : else if ( dfValue > 255.0 )
1824 11506 : *pbValue = 255;
1825 : else
1826 256119 : *pbValue = (GByte)(0.5 + dfValue);
1827 :
1828 267633 : return TRUE;
1829 : }
1830 :
1831 262207 : static int GWKCubicResampleNoMasksShort( GDALWarpKernel *poWK, int iBand,
1832 : double dfSrcX, double dfSrcY,
1833 : GInt16 *piValue )
1834 :
1835 : {
1836 262207 : int iSrcX = (int) floor( dfSrcX - 0.5 );
1837 262207 : int iSrcY = (int) floor( dfSrcY - 0.5 );
1838 262207 : int iSrcOffset = iSrcX + iSrcY * poWK->nSrcXSize;
1839 262207 : double dfDeltaX = dfSrcX - 0.5 - iSrcX;
1840 262207 : double dfDeltaY = dfSrcY - 0.5 - iSrcY;
1841 262207 : double dfDeltaX2 = dfDeltaX * dfDeltaX;
1842 262207 : double dfDeltaY2 = dfDeltaY * dfDeltaY;
1843 262207 : double dfDeltaX3 = dfDeltaX2 * dfDeltaX;
1844 262207 : double dfDeltaY3 = dfDeltaY2 * dfDeltaY;
1845 : double adfValue[4];
1846 : int i;
1847 :
1848 : // Get the bilinear interpolation at the image borders
1849 262207 : if ( iSrcX - 1 < 0 || iSrcX + 2 >= poWK->nSrcXSize
1850 : || iSrcY - 1 < 0 || iSrcY + 2 >= poWK->nSrcYSize )
1851 : return GWKBilinearResampleNoMasksShort( poWK, iBand, dfSrcX, dfSrcY,
1852 9186 : piValue);
1853 :
1854 1265105 : for ( i = -1; i < 3; i++ )
1855 : {
1856 1012084 : int iOffset = iSrcOffset + i * poWK->nSrcXSize;
1857 :
1858 11132924 : adfValue[i + 1] =CubicConvolution(dfDeltaX, dfDeltaX2, dfDeltaX3,
1859 : (double)((GInt16 *)poWK->papabySrcImage[iBand])[iOffset - 1],
1860 : (double)((GInt16 *)poWK->papabySrcImage[iBand])[iOffset],
1861 : (double)((GInt16 *)poWK->papabySrcImage[iBand])[iOffset + 1],
1862 11132924 : (double)((GInt16 *)poWK->papabySrcImage[iBand])[iOffset + 2]);
1863 : }
1864 :
1865 2024168 : *piValue = (GInt16)CubicConvolution(dfDeltaY, dfDeltaY2, dfDeltaY3,
1866 2024168 : adfValue[0], adfValue[1], adfValue[2], adfValue[3]);
1867 :
1868 253021 : return TRUE;
1869 : }
1870 :
1871 :
1872 : /************************************************************************/
1873 : /* GWKLanczosSinc() */
1874 : /************************************************************************/
1875 :
1876 : /*
1877 : * Lanczos windowed sinc interpolation kernel with radius r.
1878 : * /
1879 : * | sinc(x) * sinc(x/r), if |x| < r
1880 : * L(x) = | 1, if x = 0 ,
1881 : * | 0, otherwise
1882 : * \
1883 : *
1884 : * where sinc(x) = sin(PI * x) / (PI * x).
1885 : */
1886 :
1887 : #define GWK_PI 3.14159265358979323846
1888 3844908 : static double GWKLanczosSinc( double dfX, double dfR )
1889 : {
1890 3844908 : if ( fabs(dfX) > dfR )
1891 533142 : return 0.0;
1892 3311766 : if ( dfX == 0.0 )
1893 385 : return 1.0;
1894 :
1895 3311381 : double dfPIX = GWK_PI * dfX;
1896 3311381 : return ( sin(dfPIX) / dfPIX ) * ( sin(dfPIX / dfR) * dfR / dfPIX );
1897 : }
1898 : #undef GWK_PI
1899 :
1900 : /************************************************************************/
1901 : /* GWKBSpline() */
1902 : /************************************************************************/
1903 :
1904 4326797 : static double GWKBSpline( double x )
1905 : {
1906 4326797 : double xp2 = x + 2.0;
1907 4326797 : double xp1 = x + 1.0;
1908 4326797 : double xm1 = x - 1.0;
1909 :
1910 : // This will most likely be used, so we'll compute it ahead of time to
1911 : // avoid stalling the processor
1912 4326797 : double xp2c = xp2 * xp2 * xp2;
1913 :
1914 : // Note that the test is computed only if it is needed
1915 : return (((xp2 > 0.0)?((xp1 > 0.0)?((x > 0.0)?((xm1 > 0.0)?
1916 : -4.0 * xm1*xm1*xm1:0.0) +
1917 : 6.0 * x*x*x:0.0) +
1918 : -4.0 * xp1*xp1*xp1:0.0) +
1919 4326797 : xp2c:0.0) ) * 0.166666666666666666666;
1920 : }
1921 :
1922 :
1923 : typedef struct
1924 : {
1925 : // Space for saved X weights
1926 : double *padfWeightsX;
1927 : char *panCalcX;
1928 :
1929 : // Space for saving a row of pixels
1930 : double *padfRowDensity;
1931 : double *padfRowReal;
1932 : double *padfRowImag;
1933 : } GWKResampleWrkStruct;
1934 :
1935 87 : static GWKResampleWrkStruct* GWKResampleCreateWrkStruct(GDALWarpKernel *poWK)
1936 : {
1937 87 : int nXDist = ( poWK->nXRadius + 1 ) * 2;
1938 :
1939 : GWKResampleWrkStruct* psWrkStruct =
1940 87 : (GWKResampleWrkStruct*)CPLMalloc(sizeof(GWKResampleWrkStruct));
1941 :
1942 : // Alloc space for saved X weights
1943 87 : psWrkStruct->padfWeightsX = (double *)CPLCalloc( nXDist, sizeof(double) );
1944 87 : psWrkStruct->panCalcX = (char *)CPLMalloc( nXDist * sizeof(char) );
1945 :
1946 : // Alloc space for saving a row of pixels
1947 87 : psWrkStruct->padfRowDensity = (double *)CPLCalloc( nXDist, sizeof(double) );
1948 87 : psWrkStruct->padfRowReal = (double *)CPLCalloc( nXDist, sizeof(double) );
1949 87 : psWrkStruct->padfRowImag = (double *)CPLCalloc( nXDist, sizeof(double) );
1950 :
1951 87 : return psWrkStruct;
1952 : }
1953 :
1954 87 : static void GWKResampleDeleteWrkStruct(GWKResampleWrkStruct* psWrkStruct)
1955 : {
1956 87 : CPLFree( psWrkStruct->padfWeightsX );
1957 87 : CPLFree( psWrkStruct->panCalcX );
1958 87 : CPLFree( psWrkStruct->padfRowDensity );
1959 87 : CPLFree( psWrkStruct->padfRowReal );
1960 87 : CPLFree( psWrkStruct->padfRowImag );
1961 87 : CPLFree( psWrkStruct );
1962 87 : }
1963 :
1964 : /************************************************************************/
1965 : /* GWKResample() */
1966 : /************************************************************************/
1967 :
1968 279384 : static int GWKResample( GDALWarpKernel *poWK, int iBand,
1969 : double dfSrcX, double dfSrcY,
1970 : double *pdfDensity,
1971 : double *pdfReal, double *pdfImag,
1972 : const GWKResampleWrkStruct* psWrkStruct )
1973 :
1974 : {
1975 : // Save as local variables to avoid following pointers in loops
1976 279384 : int nSrcXSize = poWK->nSrcXSize;
1977 279384 : int nSrcYSize = poWK->nSrcYSize;
1978 :
1979 279384 : double dfAccumulatorReal = 0.0, dfAccumulatorImag = 0.0;
1980 279384 : double dfAccumulatorDensity = 0.0;
1981 279384 : double dfAccumulatorWeight = 0.0;
1982 279384 : int iSrcX = (int) floor( dfSrcX - 0.5 );
1983 279384 : int iSrcY = (int) floor( dfSrcY - 0.5 );
1984 279384 : int iSrcOffset = iSrcX + iSrcY * nSrcXSize;
1985 279384 : double dfDeltaX = dfSrcX - 0.5 - iSrcX;
1986 279384 : double dfDeltaY = dfSrcY - 0.5 - iSrcY;
1987 279384 : int eResample = poWK->eResample;
1988 :
1989 : double dfXScale, dfYScale;
1990 : double dfXFilter, dfYFilter;
1991 : int nXRadius, nFiltInitX;
1992 : int nYRadius, nFiltInitY;
1993 :
1994 279384 : dfXScale = poWK->dfXScale;
1995 279384 : dfYScale = poWK->dfYScale;
1996 279384 : nXRadius = poWK->nXRadius;
1997 279384 : nYRadius = poWK->nYRadius;
1998 279384 : nFiltInitX = poWK->nFiltInitX;
1999 279384 : nFiltInitY = poWK->nFiltInitY;
2000 279384 : dfXFilter = poWK->dfXFilter;
2001 279384 : dfYFilter = poWK->dfYFilter;
2002 :
2003 : int i, j;
2004 279384 : int nXDist = ( nXRadius + 1 ) * 2;
2005 :
2006 : // Space for saved X weights
2007 279384 : double *padfWeightsX = psWrkStruct->padfWeightsX;
2008 279384 : char *panCalcX = psWrkStruct->panCalcX;
2009 :
2010 : // Space for saving a row of pixels
2011 279384 : double *padfRowDensity = psWrkStruct->padfRowDensity;
2012 279384 : double *padfRowReal = psWrkStruct->padfRowReal;
2013 279384 : double *padfRowImag = psWrkStruct->padfRowImag;
2014 :
2015 : // Mark as needing calculation (don't calculate the weights yet,
2016 : // because a mask may render it unnecessary)
2017 279384 : memset( panCalcX, FALSE, nXDist * sizeof(char) );
2018 :
2019 : // Loop over pixel rows in the kernel
2020 2233398 : for ( j = nFiltInitY; j <= nYRadius; ++j )
2021 : {
2022 1954014 : int iRowOffset, nXMin = nFiltInitX, nXMax = nXRadius;
2023 : double dfWeight1;
2024 :
2025 : // Skip sampling over edge of image
2026 1954014 : if ( iSrcY + j < 0 || iSrcY + j >= nSrcYSize )
2027 30287 : continue;
2028 :
2029 : // Invariant; needs calculation only once per row
2030 1923727 : iRowOffset = iSrcOffset + j * nSrcXSize + nFiltInitX;
2031 :
2032 : // Make sure we don't read before or after the source array.
2033 1923727 : if ( iRowOffset < 0 )
2034 : {
2035 1316 : nXMin = nXMin - iRowOffset;
2036 1316 : iRowOffset = 0;
2037 : }
2038 1923727 : if ( iRowOffset + nXDist >= nSrcXSize*nSrcYSize )
2039 : {
2040 891 : nXMax = nSrcXSize*nSrcYSize - iRowOffset + nXMin - 1;
2041 891 : nXMax -= (nXMax-nXMin+1) % 2;
2042 : }
2043 :
2044 : // Get pixel values
2045 1923727 : if ( !GWKGetPixelRow( poWK, iBand, iRowOffset, (nXMax-nXMin+2)/2,
2046 : padfRowDensity, padfRowReal, padfRowImag ) )
2047 0 : continue;
2048 :
2049 : // Select the resampling algorithm
2050 1923727 : if ( eResample == GRA_CubicSpline )
2051 : // Calculate the Y weight
2052 : dfWeight1 = ( dfYScale < 1.0 ) ?
2053 : GWKBSpline(((double)j) * dfYScale) * dfYScale :
2054 1647 : GWKBSpline(((double)j) - dfDeltaY);
2055 1922080 : else if ( eResample == GRA_Lanczos )
2056 : dfWeight1 = ( dfYScale < 1.0 ) ?
2057 : GWKLanczosSinc(j * dfYScale, dfYFilter) * dfYScale :
2058 1922080 : GWKLanczosSinc(j - dfDeltaY, dfYFilter);
2059 : else
2060 0 : return FALSE;
2061 :
2062 : // Iterate over pixels in row
2063 15381135 : for (i = nXMin; i <= nXMax; ++i )
2064 : {
2065 : double dfWeight2;
2066 :
2067 : // Skip sampling at edge of image OR if pixel has zero density
2068 26721797 : if ( iSrcX + i < 0 || iSrcX + i >= nSrcXSize
2069 13264389 : || padfRowDensity[i-nXMin] < 0.000000001 )
2070 193019 : continue;
2071 :
2072 : // Make or use a cached set of weights for this row
2073 13264389 : if ( panCalcX[i-nFiltInitX] )
2074 : // Use saved weight value instead of recomputing it
2075 11339707 : dfWeight2 = dfWeight1 * padfWeightsX[i-nFiltInitX];
2076 : else
2077 : {
2078 : // Choose among possible algorithms
2079 1924682 : if ( eResample == GRA_CubicSpline )
2080 : // Calculate & save the X weight
2081 1854 : padfWeightsX[i-nFiltInitX] = dfWeight2 = (dfXScale < 1.0 ) ?
2082 : GWKBSpline((double)i * dfXScale) * dfXScale :
2083 1854 : GWKBSpline(dfDeltaX - (double)i);
2084 1922828 : else if ( eResample == GRA_Lanczos )
2085 : // Calculate & save the X weight
2086 1922828 : padfWeightsX[i-nFiltInitX] = dfWeight2 = (dfXScale < 1.0 ) ?
2087 : GWKLanczosSinc(i * dfXScale, dfXFilter) * dfXScale :
2088 1922828 : GWKLanczosSinc(i - dfDeltaX, dfXFilter);
2089 : else
2090 0 : return FALSE;
2091 :
2092 1924682 : dfWeight2 *= dfWeight1;
2093 1924682 : panCalcX[i-nFiltInitX] = TRUE;
2094 : }
2095 :
2096 : // Accumulate!
2097 13264389 : dfAccumulatorReal += padfRowReal[i-nXMin] * dfWeight2;
2098 13264389 : dfAccumulatorImag += padfRowImag[i-nXMin] * dfWeight2;
2099 13264389 : dfAccumulatorDensity += padfRowDensity[i-nXMin] * dfWeight2;
2100 13264389 : dfAccumulatorWeight += dfWeight2;
2101 : }
2102 : }
2103 :
2104 279384 : if ( dfAccumulatorWeight < 0.000001 || dfAccumulatorDensity < 0.000001 )
2105 : {
2106 0 : *pdfDensity = 0.0;
2107 0 : return FALSE;
2108 : }
2109 :
2110 : // Calculate the output taking into account weighting
2111 557870 : if ( dfAccumulatorWeight < 0.99999 || dfAccumulatorWeight > 1.00001 )
2112 : {
2113 278486 : *pdfReal = dfAccumulatorReal / dfAccumulatorWeight;
2114 278486 : *pdfImag = dfAccumulatorImag / dfAccumulatorWeight;
2115 278486 : *pdfDensity = dfAccumulatorDensity / dfAccumulatorWeight;
2116 : }
2117 : else
2118 : {
2119 898 : *pdfReal = dfAccumulatorReal;
2120 898 : *pdfImag = dfAccumulatorImag;
2121 898 : *pdfDensity = dfAccumulatorDensity;
2122 : }
2123 :
2124 279384 : return TRUE;
2125 : }
2126 :
2127 278207 : static int GWKCubicSplineResampleNoMasksByte( GDALWarpKernel *poWK, int iBand,
2128 : double dfSrcX, double dfSrcY,
2129 : GByte *pbValue, double *padfBSpline )
2130 :
2131 : {
2132 : // Commonly used; save locally
2133 278207 : int nSrcXSize = poWK->nSrcXSize;
2134 278207 : int nSrcYSize = poWK->nSrcYSize;
2135 :
2136 278207 : double dfAccumulator = 0.0;
2137 278207 : int iSrcX = (int) floor( dfSrcX - 0.5 );
2138 278207 : int iSrcY = (int) floor( dfSrcY - 0.5 );
2139 278207 : int iSrcOffset = iSrcX + iSrcY * nSrcXSize;
2140 278207 : double dfDeltaX = dfSrcX - 0.5 - iSrcX;
2141 278207 : double dfDeltaY = dfSrcY - 0.5 - iSrcY;
2142 :
2143 278207 : double dfXScale = poWK->dfXScale;
2144 278207 : double dfYScale = poWK->dfYScale;
2145 278207 : int nXRadius = poWK->nXRadius;
2146 278207 : int nYRadius = poWK->nYRadius;
2147 :
2148 278207 : GByte* pabySrcBand = poWK->papabySrcImage[iBand];
2149 :
2150 : // Politely refusing to process invalid coordinates or obscenely small image
2151 278207 : if ( iSrcX >= nSrcXSize || iSrcY >= nSrcYSize
2152 : || nXRadius > nSrcXSize || nYRadius > nSrcYSize )
2153 1 : return GWKBilinearResampleNoMasksByte( poWK, iBand, dfSrcX, dfSrcY, pbValue);
2154 :
2155 : // Loop over all rows in the kernel
2156 : int j, jC;
2157 1391030 : for ( jC = 0, j = 1 - nYRadius; j <= nYRadius; ++j, ++jC )
2158 : {
2159 : int iSampJ;
2160 : // Calculate the Y weight
2161 : double dfWeight1 = ( dfYScale < 1.0 ) ?
2162 : GWKBSpline((double)j * dfYScale) * dfYScale :
2163 1112824 : GWKBSpline((double)j - dfDeltaY);
2164 :
2165 : // Flip sampling over edge of image
2166 1112824 : if ( iSrcY + j < 0 )
2167 6700 : iSampJ = iSrcOffset - (iSrcY + j) * nSrcXSize;
2168 1106124 : else if ( iSrcY + j >= nSrcYSize )
2169 549 : iSampJ = iSrcOffset + (2*nSrcYSize - 2*iSrcY - j - 1) * nSrcXSize;
2170 : else
2171 1105575 : iSampJ = iSrcOffset + j * nSrcXSize;
2172 :
2173 : // Loop over all pixels in the row
2174 : int i, iC;
2175 5564120 : for ( iC = 0, i = 1 - nXRadius; i <= nXRadius; ++i, ++iC )
2176 : {
2177 : int iSampI;
2178 : double dfWeight2;
2179 :
2180 : // Flip sampling over edge of image
2181 4451296 : if ( iSrcX + i < 0 )
2182 26704 : iSampI = -iSrcX - i;
2183 4424592 : else if ( iSrcX + i >= nSrcXSize )
2184 2224 : iSampI = 2*nSrcXSize - 2*iSrcX - i - 1;
2185 : else
2186 4422368 : iSampI = i;
2187 :
2188 : // Make a cached set of GWKBSpline values
2189 4451296 : if( jC == 0 )
2190 : {
2191 : // Calculate & save the X weight
2192 1112824 : dfWeight2 = padfBSpline[iC] = ((dfXScale < 1.0 ) ?
2193 : GWKBSpline((double)i * dfXScale) * dfXScale :
2194 1112824 : GWKBSpline(dfDeltaX - (double)i));
2195 1112824 : dfWeight2 *= dfWeight1;
2196 : }
2197 : else
2198 3338472 : dfWeight2 = dfWeight1 * padfBSpline[iC];
2199 :
2200 : // Retrieve the pixel & accumulate
2201 4451296 : dfAccumulator += (double)pabySrcBand[iSampI+iSampJ] * dfWeight2;
2202 : }
2203 : }
2204 :
2205 278206 : if ( dfAccumulator < 0.0 )
2206 0 : *pbValue = 0;
2207 278206 : else if ( dfAccumulator > 255.0 )
2208 23 : *pbValue = 255;
2209 : else
2210 278183 : *pbValue = (GByte)(0.5 + dfAccumulator);
2211 :
2212 278206 : return TRUE;
2213 : }
2214 :
2215 262207 : static int GWKCubicSplineResampleNoMasksShort( GDALWarpKernel *poWK, int iBand,
2216 : double dfSrcX, double dfSrcY,
2217 : GInt16 *piValue, double *padfBSpline )
2218 :
2219 : {
2220 : //Save src size to local var
2221 262207 : int nSrcXSize = poWK->nSrcXSize;
2222 262207 : int nSrcYSize = poWK->nSrcYSize;
2223 :
2224 262207 : double dfAccumulator = 0.0;
2225 262207 : int iSrcX = (int) floor( dfSrcX - 0.5 );
2226 262207 : int iSrcY = (int) floor( dfSrcY - 0.5 );
2227 262207 : int iSrcOffset = iSrcX + iSrcY * nSrcXSize;
2228 262207 : double dfDeltaX = dfSrcX - 0.5 - iSrcX;
2229 262207 : double dfDeltaY = dfSrcY - 0.5 - iSrcY;
2230 :
2231 262207 : double dfXScale = poWK->dfXScale;
2232 262207 : double dfYScale = poWK->dfYScale;
2233 262207 : int nXRadius = poWK->nXRadius;
2234 262207 : int nYRadius = poWK->nYRadius;
2235 :
2236 : // Save band array pointer to local var; cast here instead of later
2237 262207 : GInt16* pabySrcBand = ((GInt16 *)poWK->papabySrcImage[iBand]);
2238 :
2239 : // Politely refusing to process invalid coordinates or obscenely small image
2240 262207 : if ( iSrcX >= nSrcXSize || iSrcY >= nSrcYSize
2241 : || nXRadius > nSrcXSize || nYRadius > nSrcYSize )
2242 1 : return GWKBilinearResampleNoMasksShort( poWK, iBand, dfSrcX, dfSrcY, piValue);
2243 :
2244 : // Loop over all pixels in the kernel
2245 : int j, jC;
2246 1311030 : for ( jC = 0, j = 1 - nYRadius; j <= nYRadius; ++j, ++jC )
2247 : {
2248 : int iSampJ;
2249 :
2250 : // Calculate the Y weight
2251 : double dfWeight1 = ( dfYScale < 1.0 ) ?
2252 : GWKBSpline((double)j * dfYScale) * dfYScale :
2253 1048824 : GWKBSpline((double)j - dfDeltaY);
2254 :
2255 : // Flip sampling over edge of image
2256 1048824 : if ( iSrcY + j < 0 )
2257 6180 : iSampJ = iSrcOffset - (iSrcY + j) * nSrcXSize;
2258 1042644 : else if ( iSrcY + j >= nSrcYSize )
2259 29 : iSampJ = iSrcOffset + (2*nSrcYSize - 2*iSrcY - j - 1) * nSrcXSize;
2260 : else
2261 1042615 : iSampJ = iSrcOffset + j * nSrcXSize;
2262 :
2263 : // Loop over all pixels in row
2264 : int i, iC;
2265 5244120 : for ( iC = 0, i = 1 - nXRadius; i <= nXRadius; ++i, ++iC )
2266 : {
2267 : int iSampI;
2268 : double dfWeight2;
2269 :
2270 : // Flip sampling over edge of image
2271 4195296 : if ( iSrcX + i < 0 )
2272 24624 : iSampI = -iSrcX - i;
2273 4170672 : else if(iSrcX + i >= nSrcXSize)
2274 144 : iSampI = 2*nSrcXSize - 2*iSrcX - i - 1;
2275 : else
2276 4170528 : iSampI = i;
2277 :
2278 : // Make a cached set of GWKBSpline values
2279 4195296 : if ( jC == 0 )
2280 : {
2281 : // Calculate & save the X weight
2282 1048824 : dfWeight2 = padfBSpline[iC] = ((dfXScale < 1.0 ) ?
2283 : GWKBSpline((double)i * dfXScale) * dfXScale :
2284 1048824 : GWKBSpline(dfDeltaX - (double)i));
2285 1048824 : dfWeight2 *= dfWeight1;
2286 : } else
2287 3146472 : dfWeight2 = dfWeight1 * padfBSpline[iC];
2288 :
2289 4195296 : dfAccumulator += (double)pabySrcBand[iSampI + iSampJ] * dfWeight2;
2290 : }
2291 : }
2292 :
2293 262206 : *piValue = (GInt16)(0.5 + dfAccumulator);
2294 :
2295 262206 : return TRUE;
2296 : }
2297 :
2298 : /************************************************************************/
2299 : /* GWKGeneralCase() */
2300 : /* */
2301 : /* This is the most general case. It attempts to handle all */
2302 : /* possible features with relatively little concern for */
2303 : /* efficiency. */
2304 : /************************************************************************/
2305 :
2306 189 : static CPLErr GWKGeneralCase( GDALWarpKernel *poWK )
2307 :
2308 : {
2309 : int iDstY;
2310 189 : int nDstXSize = poWK->nDstXSize, nDstYSize = poWK->nDstYSize;
2311 189 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
2312 189 : CPLErr eErr = CE_None;
2313 :
2314 : CPLDebug( "GDAL", "GDALWarpKernel()::GWKGeneralCase()\n"
2315 : "Src=%d,%d,%dx%d Dst=%d,%d,%dx%d",
2316 : poWK->nSrcXOff, poWK->nSrcYOff,
2317 : poWK->nSrcXSize, poWK->nSrcYSize,
2318 : poWK->nDstXOff, poWK->nDstYOff,
2319 189 : poWK->nDstXSize, poWK->nDstYSize );
2320 :
2321 189 : if( !poWK->pfnProgress( poWK->dfProgressBase, "", poWK->pProgress ) )
2322 : {
2323 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
2324 0 : return CE_Failure;
2325 : }
2326 :
2327 : /* -------------------------------------------------------------------- */
2328 : /* Allocate x,y,z coordinate arrays for transformation ... one */
2329 : /* scanlines worth of positions. */
2330 : /* -------------------------------------------------------------------- */
2331 : double *padfX, *padfY, *padfZ;
2332 : int *pabSuccess;
2333 :
2334 189 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
2335 189 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
2336 189 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
2337 189 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
2338 :
2339 189 : GWKResampleWrkStruct* psWrkStruct = NULL;
2340 189 : if (poWK->eResample == GRA_CubicSpline
2341 : || poWK->eResample == GRA_Lanczos )
2342 : {
2343 87 : psWrkStruct = GWKResampleCreateWrkStruct(poWK);
2344 : }
2345 :
2346 : /* ==================================================================== */
2347 : /* Loop over output lines. */
2348 : /* ==================================================================== */
2349 1838 : for( iDstY = 0; iDstY < nDstYSize && eErr == CE_None; iDstY++ )
2350 : {
2351 : int iDstX;
2352 :
2353 : /* -------------------------------------------------------------------- */
2354 : /* Setup points to transform to source image space. */
2355 : /* -------------------------------------------------------------------- */
2356 530372 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
2357 : {
2358 528723 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
2359 528723 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
2360 528723 : padfZ[iDstX] = 0.0;
2361 : }
2362 :
2363 : /* -------------------------------------------------------------------- */
2364 : /* Transform the points from destination pixel/line coordinates */
2365 : /* to source pixel/line coordinates. */
2366 : /* -------------------------------------------------------------------- */
2367 : poWK->pfnTransformer( poWK->pTransformerArg, TRUE, nDstXSize,
2368 1649 : padfX, padfY, padfZ, pabSuccess );
2369 :
2370 : /* ==================================================================== */
2371 : /* Loop over pixels in output scanline. */
2372 : /* ==================================================================== */
2373 530372 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
2374 : {
2375 : int iDstOffset;
2376 :
2377 528723 : if( !pabSuccess[iDstX] )
2378 0 : continue;
2379 :
2380 : /* -------------------------------------------------------------------- */
2381 : /* Figure out what pixel we want in our source raster, and skip */
2382 : /* further processing if it is well off the source image. */
2383 : /* -------------------------------------------------------------------- */
2384 : // We test against the value before casting to avoid the
2385 : // problem of asymmetric truncation effects around zero. That is
2386 : // -0.5 will be 0 when cast to an int.
2387 1057407 : if( padfX[iDstX] < poWK->nSrcXOff
2388 528684 : || padfY[iDstX] < poWK->nSrcYOff )
2389 29292 : continue;
2390 :
2391 : int iSrcX, iSrcY, iSrcOffset;
2392 :
2393 499431 : iSrcX = ((int) padfX[iDstX]) - poWK->nSrcXOff;
2394 499431 : iSrcY = ((int) padfY[iDstX]) - poWK->nSrcYOff;
2395 :
2396 : // If operating outside natural projection area, padfX/Y can be
2397 : // a very huge positive number, that becomes -2147483648 in the
2398 : // int trucation. So it is necessary to test now for non negativeness.
2399 499431 : if( iSrcX < 0 || iSrcX >= nSrcXSize || iSrcY < 0 || iSrcY >= nSrcYSize )
2400 217249 : continue;
2401 :
2402 282182 : iSrcOffset = iSrcX + iSrcY * nSrcXSize;
2403 :
2404 : /* -------------------------------------------------------------------- */
2405 : /* Do not try to apply transparent/invalid source pixels to the */
2406 : /* destination. This currently ignores the multi-pixel input */
2407 : /* of bilinear and cubic resamples. */
2408 : /* -------------------------------------------------------------------- */
2409 282182 : double dfDensity = 1.0;
2410 :
2411 282182 : if( poWK->pafUnifiedSrcDensity != NULL
2412 : && iSrcX >= 0 && iSrcY >= 0
2413 : && iSrcX < nSrcXSize && iSrcY < nSrcYSize )
2414 : {
2415 1203 : dfDensity = poWK->pafUnifiedSrcDensity[iSrcOffset];
2416 1203 : if( dfDensity < 0.00001 )
2417 1203 : continue;
2418 : }
2419 :
2420 280979 : if( poWK->panUnifiedSrcValid != NULL
2421 : && iSrcX >= 0 && iSrcY >= 0
2422 : && iSrcX < nSrcXSize && iSrcY < nSrcYSize
2423 0 : && !(poWK->panUnifiedSrcValid[iSrcOffset>>5]
2424 : & (0x01 << (iSrcOffset & 0x1f))) )
2425 0 : continue;
2426 :
2427 : /* ==================================================================== */
2428 : /* Loop processing each band. */
2429 : /* ==================================================================== */
2430 : int iBand;
2431 280979 : int bHasFoundDensity = FALSE;
2432 :
2433 280979 : iDstOffset = iDstX + iDstY * nDstXSize;
2434 561958 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
2435 : {
2436 280979 : double dfBandDensity = 0.0;
2437 280979 : double dfValueReal = 0.0;
2438 280979 : double dfValueImag = 0.0;
2439 :
2440 : /* -------------------------------------------------------------------- */
2441 : /* Collect the source value. */
2442 : /* -------------------------------------------------------------------- */
2443 281458 : if ( poWK->eResample == GRA_NearestNeighbour ||
2444 : nSrcXSize == 1 || nSrcYSize == 1)
2445 : {
2446 : GWKGetPixelValue( poWK, iBand, iSrcOffset,
2447 479 : &dfBandDensity, &dfValueReal, &dfValueImag );
2448 : }
2449 280500 : else if ( poWK->eResample == GRA_Bilinear )
2450 : {
2451 : GWKBilinearResample( poWK, iBand,
2452 558 : padfX[iDstX]-poWK->nSrcXOff,
2453 558 : padfY[iDstX]-poWK->nSrcYOff,
2454 : &dfBandDensity,
2455 1116 : &dfValueReal, &dfValueImag );
2456 : }
2457 279942 : else if ( poWK->eResample == GRA_Cubic )
2458 : {
2459 : GWKCubicResample( poWK, iBand,
2460 558 : padfX[iDstX]-poWK->nSrcXOff,
2461 558 : padfY[iDstX]-poWK->nSrcYOff,
2462 : &dfBandDensity,
2463 1116 : &dfValueReal, &dfValueImag );
2464 : }
2465 279384 : else if ( poWK->eResample == GRA_CubicSpline
2466 : || poWK->eResample == GRA_Lanczos )
2467 : {
2468 : GWKResample( poWK, iBand,
2469 279384 : padfX[iDstX]-poWK->nSrcXOff,
2470 279384 : padfY[iDstX]-poWK->nSrcYOff,
2471 : &dfBandDensity,
2472 558768 : &dfValueReal, &dfValueImag, psWrkStruct );
2473 : }
2474 :
2475 :
2476 : // If we didn't find any valid inputs skip to next band.
2477 280979 : if ( dfBandDensity < 0.0000000001 )
2478 0 : continue;
2479 :
2480 280979 : bHasFoundDensity = TRUE;
2481 :
2482 : /* -------------------------------------------------------------------- */
2483 : /* We have a computed value from the source. Now apply it to */
2484 : /* the destination pixel. */
2485 : /* -------------------------------------------------------------------- */
2486 : GWKSetPixelValue( poWK, iBand, iDstOffset,
2487 : dfBandDensity,
2488 280979 : dfValueReal, dfValueImag );
2489 :
2490 : }
2491 :
2492 280979 : if (!bHasFoundDensity)
2493 0 : continue;
2494 :
2495 : /* -------------------------------------------------------------------- */
2496 : /* Update destination density/validity masks. */
2497 : /* -------------------------------------------------------------------- */
2498 280979 : GWKOverlayDensity( poWK, iDstOffset, dfDensity );
2499 :
2500 280979 : if( poWK->panDstValid != NULL )
2501 : {
2502 : poWK->panDstValid[iDstOffset>>5] |=
2503 0 : 0x01 << (iDstOffset & 0x1f);
2504 : }
2505 :
2506 : } /* Next iDstX */
2507 :
2508 : /* -------------------------------------------------------------------- */
2509 : /* Report progress to the user, and optionally cancel out. */
2510 : /* -------------------------------------------------------------------- */
2511 1649 : if( !poWK->pfnProgress( poWK->dfProgressBase + poWK->dfProgressScale *
2512 : ((iDstY+1) / (double) nDstYSize),
2513 : "", poWK->pProgress ) )
2514 : {
2515 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
2516 0 : eErr = CE_Failure;
2517 : }
2518 : }
2519 :
2520 : /* -------------------------------------------------------------------- */
2521 : /* Cleanup and return. */
2522 : /* -------------------------------------------------------------------- */
2523 189 : CPLFree( padfX );
2524 189 : CPLFree( padfY );
2525 189 : CPLFree( padfZ );
2526 189 : CPLFree( pabSuccess );
2527 189 : if (psWrkStruct)
2528 87 : GWKResampleDeleteWrkStruct(psWrkStruct);
2529 :
2530 189 : return eErr;
2531 : }
2532 :
2533 : /************************************************************************/
2534 : /* GWKNearestNoMasksByte() */
2535 : /* */
2536 : /* Case for 8bit input data with nearest neighbour resampling */
2537 : /* without concerning about masking. Should be as fast as */
2538 : /* possible for this particular transformation type. */
2539 : /************************************************************************/
2540 :
2541 35 : static CPLErr GWKNearestNoMasksByte( GDALWarpKernel *poWK )
2542 :
2543 : {
2544 : int iDstY;
2545 35 : int nDstXSize = poWK->nDstXSize, nDstYSize = poWK->nDstYSize;
2546 35 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
2547 35 : CPLErr eErr = CE_None;
2548 :
2549 : CPLDebug( "GDAL", "GDALWarpKernel()::GWKNearestNoMasksByte()\n"
2550 : "Src=%d,%d,%dx%d Dst=%d,%d,%dx%d",
2551 : poWK->nSrcXOff, poWK->nSrcYOff,
2552 : poWK->nSrcXSize, poWK->nSrcYSize,
2553 : poWK->nDstXOff, poWK->nDstYOff,
2554 35 : poWK->nDstXSize, poWK->nDstYSize );
2555 :
2556 35 : if( !poWK->pfnProgress( poWK->dfProgressBase, "", poWK->pProgress ) )
2557 : {
2558 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
2559 0 : return CE_Failure;
2560 : }
2561 :
2562 : /* -------------------------------------------------------------------- */
2563 : /* Allocate x,y,z coordinate arrays for transformation ... one */
2564 : /* scanlines worth of positions. */
2565 : /* -------------------------------------------------------------------- */
2566 : double *padfX, *padfY, *padfZ;
2567 : int *pabSuccess;
2568 :
2569 35 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
2570 35 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
2571 35 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
2572 35 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
2573 :
2574 : /* ==================================================================== */
2575 : /* Loop over output lines. */
2576 : /* ==================================================================== */
2577 3729 : for( iDstY = 0; iDstY < nDstYSize && eErr == CE_None; iDstY++ )
2578 : {
2579 : int iDstX;
2580 :
2581 : /* -------------------------------------------------------------------- */
2582 : /* Setup points to transform to source image space. */
2583 : /* -------------------------------------------------------------------- */
2584 1746574 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
2585 : {
2586 1742880 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
2587 1742880 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
2588 1742880 : padfZ[iDstX] = 0.0;
2589 : }
2590 :
2591 : /* -------------------------------------------------------------------- */
2592 : /* Transform the points from destination pixel/line coordinates */
2593 : /* to source pixel/line coordinates. */
2594 : /* -------------------------------------------------------------------- */
2595 : poWK->pfnTransformer( poWK->pTransformerArg, TRUE, nDstXSize,
2596 3694 : padfX, padfY, padfZ, pabSuccess );
2597 :
2598 : /* ==================================================================== */
2599 : /* Loop over pixels in output scanline. */
2600 : /* ==================================================================== */
2601 1746574 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
2602 : {
2603 1742880 : if( !pabSuccess[iDstX] )
2604 153669 : continue;
2605 :
2606 : /* -------------------------------------------------------------------- */
2607 : /* Figure out what pixel we want in our source raster, and skip */
2608 : /* further processing if it is well off the source image. */
2609 : /* -------------------------------------------------------------------- */
2610 : // We test against the value before casting to avoid the
2611 : // problem of asymmetric truncation effects around zero. That is
2612 : // -0.5 will be 0 when cast to an int.
2613 3178222 : if( padfX[iDstX] < poWK->nSrcXOff
2614 1589011 : || padfY[iDstX] < poWK->nSrcYOff )
2615 120136 : continue;
2616 :
2617 : int iSrcX, iSrcY, iSrcOffset;
2618 :
2619 1469075 : iSrcX = ((int) padfX[iDstX]) - poWK->nSrcXOff;
2620 1469075 : iSrcY = ((int) padfY[iDstX]) - poWK->nSrcYOff;
2621 :
2622 : // If operating outside natural projection area, padfX/Y can be
2623 : // a very huge positive number, that becomes -2147483648 in the
2624 : // int trucation. So it is necessary to test now for non negativeness.
2625 1469075 : if( iSrcX < 0 || iSrcX >= nSrcXSize || iSrcY < 0 || iSrcY >= nSrcYSize )
2626 305331 : continue;
2627 :
2628 1163744 : iSrcOffset = iSrcX + iSrcY * nSrcXSize;
2629 :
2630 : /* ==================================================================== */
2631 : /* Loop processing each band. */
2632 : /* ==================================================================== */
2633 : int iBand;
2634 : int iDstOffset;
2635 :
2636 1163744 : iDstOffset = iDstX + iDstY * nDstXSize;
2637 :
2638 3916032 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
2639 : {
2640 2752288 : poWK->papabyDstImage[iBand][iDstOffset] =
2641 2752288 : poWK->papabySrcImage[iBand][iSrcOffset];
2642 : }
2643 : }
2644 :
2645 : /* -------------------------------------------------------------------- */
2646 : /* Report progress to the user, and optionally cancel out. */
2647 : /* -------------------------------------------------------------------- */
2648 3694 : if( !poWK->pfnProgress( poWK->dfProgressBase + poWK->dfProgressScale *
2649 : ((iDstY+1) / (double) nDstYSize),
2650 : "", poWK->pProgress ) )
2651 : {
2652 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
2653 0 : eErr = CE_Failure;
2654 : }
2655 : }
2656 :
2657 : /* -------------------------------------------------------------------- */
2658 : /* Cleanup and return. */
2659 : /* -------------------------------------------------------------------- */
2660 35 : CPLFree( padfX );
2661 35 : CPLFree( padfY );
2662 35 : CPLFree( padfZ );
2663 35 : CPLFree( pabSuccess );
2664 :
2665 35 : return eErr;
2666 : }
2667 :
2668 : /************************************************************************/
2669 : /* GWKBilinearNoMasksByte() */
2670 : /* */
2671 : /* Case for 8bit input data with bilinear resampling without */
2672 : /* concerning about masking. Should be as fast as possible */
2673 : /* for this particular transformation type. */
2674 : /************************************************************************/
2675 :
2676 10 : static CPLErr GWKBilinearNoMasksByte( GDALWarpKernel *poWK )
2677 :
2678 : {
2679 : int iDstY;
2680 10 : int nDstXSize = poWK->nDstXSize, nDstYSize = poWK->nDstYSize;
2681 10 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
2682 10 : CPLErr eErr = CE_None;
2683 :
2684 : CPLDebug( "GDAL", "GDALWarpKernel()::GWKBilinearNoMasksByte()\n"
2685 : "Src=%d,%d,%dx%d Dst=%d,%d,%dx%d",
2686 : poWK->nSrcXOff, poWK->nSrcYOff,
2687 : poWK->nSrcXSize, poWK->nSrcYSize,
2688 : poWK->nDstXOff, poWK->nDstYOff,
2689 10 : poWK->nDstXSize, poWK->nDstYSize );
2690 :
2691 10 : if( !poWK->pfnProgress( poWK->dfProgressBase, "", poWK->pProgress ) )
2692 : {
2693 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
2694 0 : return CE_Failure;
2695 : }
2696 :
2697 : /* -------------------------------------------------------------------- */
2698 : /* Allocate x,y,z coordinate arrays for transformation ... one */
2699 : /* scanlines worth of positions. */
2700 : /* -------------------------------------------------------------------- */
2701 : double *padfX, *padfY, *padfZ;
2702 : int *pabSuccess;
2703 :
2704 10 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
2705 10 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
2706 10 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
2707 10 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
2708 :
2709 : /* ==================================================================== */
2710 : /* Loop over output lines. */
2711 : /* ==================================================================== */
2712 703 : for( iDstY = 0; iDstY < nDstYSize && eErr == CE_None; iDstY++ )
2713 : {
2714 : int iDstX;
2715 :
2716 : /* -------------------------------------------------------------------- */
2717 : /* Setup points to transform to source image space. */
2718 : /* -------------------------------------------------------------------- */
2719 330036 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
2720 : {
2721 329343 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
2722 329343 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
2723 329343 : padfZ[iDstX] = 0.0;
2724 : }
2725 :
2726 : /* -------------------------------------------------------------------- */
2727 : /* Transform the points from destination pixel/line coordinates */
2728 : /* to source pixel/line coordinates. */
2729 : /* -------------------------------------------------------------------- */
2730 : poWK->pfnTransformer( poWK->pTransformerArg, TRUE, nDstXSize,
2731 693 : padfX, padfY, padfZ, pabSuccess );
2732 :
2733 : /* ==================================================================== */
2734 : /* Loop over pixels in output scanline. */
2735 : /* ==================================================================== */
2736 330036 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
2737 : {
2738 329343 : if( !pabSuccess[iDstX] )
2739 0 : continue;
2740 :
2741 : /* -------------------------------------------------------------------- */
2742 : /* Figure out what pixel we want in our source raster, and skip */
2743 : /* further processing if it is well off the source image. */
2744 : /* -------------------------------------------------------------------- */
2745 : // We test against the value before casting to avoid the
2746 : // problem of asymmetric truncation effects around zero. That is
2747 : // -0.5 will be 0 when cast to an int.
2748 658686 : if( padfX[iDstX] < poWK->nSrcXOff
2749 329343 : || padfY[iDstX] < poWK->nSrcYOff )
2750 0 : continue;
2751 :
2752 : int iSrcX, iSrcY, iSrcOffset;
2753 :
2754 329343 : iSrcX = ((int) padfX[iDstX]) - poWK->nSrcXOff;
2755 329343 : iSrcY = ((int) padfY[iDstX]) - poWK->nSrcYOff;
2756 :
2757 : // If operating outside natural projection area, padfX/Y can be
2758 : // a very huge positive number, that becomes -2147483648 in the
2759 : // int trucation. So it is necessary to test now for non negativeness.
2760 329343 : if( iSrcX < 0 || iSrcX >= nSrcXSize || iSrcY < 0 || iSrcY >= nSrcYSize )
2761 51136 : continue;
2762 :
2763 278207 : iSrcOffset = iSrcX + iSrcY * nSrcXSize;
2764 :
2765 : /* ==================================================================== */
2766 : /* Loop processing each band. */
2767 : /* ==================================================================== */
2768 : int iBand;
2769 : int iDstOffset;
2770 :
2771 278207 : iDstOffset = iDstX + iDstY * nDstXSize;
2772 :
2773 556414 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
2774 : {
2775 : GWKBilinearResampleNoMasksByte( poWK, iBand,
2776 278207 : padfX[iDstX]-poWK->nSrcXOff,
2777 278207 : padfY[iDstX]-poWK->nSrcYOff,
2778 834621 : &poWK->papabyDstImage[iBand][iDstOffset] );
2779 : }
2780 : }
2781 :
2782 : /* -------------------------------------------------------------------- */
2783 : /* Report progress to the user, and optionally cancel out. */
2784 : /* -------------------------------------------------------------------- */
2785 693 : if( !poWK->pfnProgress( poWK->dfProgressBase + poWK->dfProgressScale *
2786 : ((iDstY+1) / (double) nDstYSize),
2787 : "", poWK->pProgress ) )
2788 : {
2789 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
2790 0 : eErr = CE_Failure;
2791 : }
2792 : }
2793 :
2794 : /* -------------------------------------------------------------------- */
2795 : /* Cleanup and return. */
2796 : /* -------------------------------------------------------------------- */
2797 10 : CPLFree( padfX );
2798 10 : CPLFree( padfY );
2799 10 : CPLFree( padfZ );
2800 10 : CPLFree( pabSuccess );
2801 :
2802 10 : return eErr;
2803 : }
2804 :
2805 : /************************************************************************/
2806 : /* GWKCubicNoMasksByte() */
2807 : /* */
2808 : /* Case for 8bit input data with cubic resampling without */
2809 : /* concerning about masking. Should be as fast as possible */
2810 : /* for this particular transformation type. */
2811 : /************************************************************************/
2812 :
2813 10 : static CPLErr GWKCubicNoMasksByte( GDALWarpKernel *poWK )
2814 :
2815 : {
2816 : int iDstY;
2817 10 : int nDstXSize = poWK->nDstXSize, nDstYSize = poWK->nDstYSize;
2818 10 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
2819 10 : CPLErr eErr = CE_None;
2820 :
2821 : CPLDebug( "GDAL", "GDALWarpKernel()::GWKCubicNoMasksByte()\n"
2822 : "Src=%d,%d,%dx%d Dst=%d,%d,%dx%d",
2823 : poWK->nSrcXOff, poWK->nSrcYOff,
2824 : poWK->nSrcXSize, poWK->nSrcYSize,
2825 : poWK->nDstXOff, poWK->nDstYOff,
2826 10 : poWK->nDstXSize, poWK->nDstYSize );
2827 :
2828 10 : if( !poWK->pfnProgress( poWK->dfProgressBase, "", poWK->pProgress ) )
2829 : {
2830 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
2831 0 : return CE_Failure;
2832 : }
2833 :
2834 : /* -------------------------------------------------------------------- */
2835 : /* Allocate x,y,z coordinate arrays for transformation ... one */
2836 : /* scanlines worth of positions. */
2837 : /* -------------------------------------------------------------------- */
2838 : double *padfX, *padfY, *padfZ;
2839 : int *pabSuccess;
2840 :
2841 10 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
2842 10 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
2843 10 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
2844 10 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
2845 :
2846 : /* ==================================================================== */
2847 : /* Loop over output lines. */
2848 : /* ==================================================================== */
2849 703 : for( iDstY = 0; iDstY < nDstYSize && eErr == CE_None; iDstY++ )
2850 : {
2851 : int iDstX;
2852 :
2853 : /* -------------------------------------------------------------------- */
2854 : /* Setup points to transform to source image space. */
2855 : /* -------------------------------------------------------------------- */
2856 330036 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
2857 : {
2858 329343 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
2859 329343 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
2860 329343 : padfZ[iDstX] = 0.0;
2861 : }
2862 :
2863 : /* -------------------------------------------------------------------- */
2864 : /* Transform the points from destination pixel/line coordinates */
2865 : /* to source pixel/line coordinates. */
2866 : /* -------------------------------------------------------------------- */
2867 : poWK->pfnTransformer( poWK->pTransformerArg, TRUE, nDstXSize,
2868 693 : padfX, padfY, padfZ, pabSuccess );
2869 :
2870 : /* ==================================================================== */
2871 : /* Loop over pixels in output scanline. */
2872 : /* ==================================================================== */
2873 330036 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
2874 : {
2875 329343 : if( !pabSuccess[iDstX] )
2876 0 : continue;
2877 :
2878 : /* -------------------------------------------------------------------- */
2879 : /* Figure out what pixel we want in our source raster, and skip */
2880 : /* further processing if it is well off the source image. */
2881 : /* -------------------------------------------------------------------- */
2882 : // We test against the value before casting to avoid the
2883 : // problem of asymmetric truncation effects around zero. That is
2884 : // -0.5 will be 0 when cast to an int.
2885 658686 : if( padfX[iDstX] < poWK->nSrcXOff
2886 329343 : || padfY[iDstX] < poWK->nSrcYOff )
2887 0 : continue;
2888 :
2889 : int iSrcX, iSrcY, iSrcOffset;
2890 :
2891 329343 : iSrcX = ((int) padfX[iDstX]) - poWK->nSrcXOff;
2892 329343 : iSrcY = ((int) padfY[iDstX]) - poWK->nSrcYOff;
2893 :
2894 : // If operating outside natural projection area, padfX/Y can be
2895 : // a very huge positive number, that becomes -2147483648 in the
2896 : // int trucation. So it is necessary to test now for non negativeness.
2897 329343 : if( iSrcX < 0 || iSrcX >= nSrcXSize || iSrcY < 0 || iSrcY >= nSrcYSize )
2898 51136 : continue;
2899 :
2900 278207 : iSrcOffset = iSrcX + iSrcY * nSrcXSize;
2901 :
2902 : /* ==================================================================== */
2903 : /* Loop processing each band. */
2904 : /* ==================================================================== */
2905 : int iBand;
2906 : int iDstOffset;
2907 :
2908 278207 : iDstOffset = iDstX + iDstY * nDstXSize;
2909 :
2910 556414 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
2911 : {
2912 : GWKCubicResampleNoMasksByte( poWK, iBand,
2913 278207 : padfX[iDstX]-poWK->nSrcXOff,
2914 278207 : padfY[iDstX]-poWK->nSrcYOff,
2915 834621 : &poWK->papabyDstImage[iBand][iDstOffset] );
2916 : }
2917 : }
2918 :
2919 : /* -------------------------------------------------------------------- */
2920 : /* Report progress to the user, and optionally cancel out. */
2921 : /* -------------------------------------------------------------------- */
2922 693 : if( !poWK->pfnProgress( poWK->dfProgressBase + poWK->dfProgressScale *
2923 : ((iDstY+1) / (double) nDstYSize),
2924 : "", poWK->pProgress ) )
2925 : {
2926 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
2927 0 : eErr = CE_Failure;
2928 : }
2929 : }
2930 :
2931 : /* -------------------------------------------------------------------- */
2932 : /* Cleanup and return. */
2933 : /* -------------------------------------------------------------------- */
2934 10 : CPLFree( padfX );
2935 10 : CPLFree( padfY );
2936 10 : CPLFree( padfZ );
2937 10 : CPLFree( pabSuccess );
2938 :
2939 10 : return eErr;
2940 : }
2941 :
2942 : /************************************************************************/
2943 : /* GWKCubicSplineNoMasksByte() */
2944 : /* */
2945 : /* Case for 8bit input data with cubic spline resampling without */
2946 : /* concerning about masking. Should be as fast as possible */
2947 : /* for this particular transformation type. */
2948 : /************************************************************************/
2949 :
2950 10 : static CPLErr GWKCubicSplineNoMasksByte( GDALWarpKernel *poWK )
2951 :
2952 : {
2953 : int iDstY;
2954 10 : int nDstXSize = poWK->nDstXSize, nDstYSize = poWK->nDstYSize;
2955 10 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
2956 10 : CPLErr eErr = CE_None;
2957 :
2958 : CPLDebug( "GDAL", "GDALWarpKernel()::GWKCubicSplineNoMasksByte()\n"
2959 : "Src=%d,%d,%dx%d Dst=%d,%d,%dx%d",
2960 : poWK->nSrcXOff, poWK->nSrcYOff,
2961 : poWK->nSrcXSize, poWK->nSrcYSize,
2962 : poWK->nDstXOff, poWK->nDstYOff,
2963 10 : poWK->nDstXSize, poWK->nDstYSize );
2964 :
2965 10 : if( !poWK->pfnProgress( poWK->dfProgressBase, "", poWK->pProgress ) )
2966 : {
2967 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
2968 0 : return CE_Failure;
2969 : }
2970 :
2971 : /* -------------------------------------------------------------------- */
2972 : /* Allocate x,y,z coordinate arrays for transformation ... one */
2973 : /* scanlines worth of positions. */
2974 : /* -------------------------------------------------------------------- */
2975 : double *padfX, *padfY, *padfZ;
2976 : int *pabSuccess;
2977 :
2978 10 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
2979 10 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
2980 10 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
2981 10 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
2982 :
2983 10 : int nXRadius = poWK->nXRadius;
2984 10 : double *padfBSpline = (double *)CPLCalloc( nXRadius * 2, sizeof(double) );
2985 :
2986 : /* ==================================================================== */
2987 : /* Loop over output lines. */
2988 : /* ==================================================================== */
2989 703 : for( iDstY = 0; iDstY < nDstYSize && eErr == CE_None; iDstY++ )
2990 : {
2991 : int iDstX;
2992 :
2993 : /* -------------------------------------------------------------------- */
2994 : /* Setup points to transform to source image space. */
2995 : /* -------------------------------------------------------------------- */
2996 330036 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
2997 : {
2998 329343 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
2999 329343 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
3000 329343 : padfZ[iDstX] = 0.0;
3001 : }
3002 :
3003 : /* -------------------------------------------------------------------- */
3004 : /* Transform the points from destination pixel/line coordinates */
3005 : /* to source pixel/line coordinates. */
3006 : /* -------------------------------------------------------------------- */
3007 : poWK->pfnTransformer( poWK->pTransformerArg, TRUE, nDstXSize,
3008 693 : padfX, padfY, padfZ, pabSuccess );
3009 :
3010 : /* ==================================================================== */
3011 : /* Loop over pixels in output scanline. */
3012 : /* ==================================================================== */
3013 330036 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3014 : {
3015 329343 : if( !pabSuccess[iDstX] )
3016 0 : continue;
3017 :
3018 : /* -------------------------------------------------------------------- */
3019 : /* Figure out what pixel we want in our source raster, and skip */
3020 : /* further processing if it is well off the source image. */
3021 : /* -------------------------------------------------------------------- */
3022 : // We test against the value before casting to avoid the
3023 : // problem of asymmetric truncation effects around zero. That is
3024 : // -0.5 will be 0 when cast to an int.
3025 658686 : if( padfX[iDstX] < poWK->nSrcXOff
3026 329343 : || padfY[iDstX] < poWK->nSrcYOff )
3027 0 : continue;
3028 :
3029 : int iSrcX, iSrcY, iSrcOffset;
3030 :
3031 329343 : iSrcX = ((int) padfX[iDstX]) - poWK->nSrcXOff;
3032 329343 : iSrcY = ((int) padfY[iDstX]) - poWK->nSrcYOff;
3033 :
3034 : // If operating outside natural projection area, padfX/Y can be
3035 : // a very huge positive number, that becomes -2147483648 in the
3036 : // int trucation. So it is necessary to test now for non negativeness.
3037 329343 : if( iSrcX < 0 || iSrcX >= nSrcXSize || iSrcY < 0 || iSrcY >= nSrcYSize )
3038 51136 : continue;
3039 :
3040 278207 : iSrcOffset = iSrcX + iSrcY * nSrcXSize;
3041 :
3042 : /* ==================================================================== */
3043 : /* Loop processing each band. */
3044 : /* ==================================================================== */
3045 : int iBand;
3046 : int iDstOffset;
3047 :
3048 278207 : iDstOffset = iDstX + iDstY * nDstXSize;
3049 :
3050 556414 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
3051 : {
3052 : GWKCubicSplineResampleNoMasksByte( poWK, iBand,
3053 278207 : padfX[iDstX]-poWK->nSrcXOff,
3054 278207 : padfY[iDstX]-poWK->nSrcYOff,
3055 278207 : &poWK->papabyDstImage[iBand][iDstOffset],
3056 834621 : padfBSpline);
3057 : }
3058 : }
3059 :
3060 : /* -------------------------------------------------------------------- */
3061 : /* Report progress to the user, and optionally cancel out. */
3062 : /* -------------------------------------------------------------------- */
3063 693 : if( !poWK->pfnProgress( poWK->dfProgressBase + poWK->dfProgressScale *
3064 : ((iDstY+1) / (double) nDstYSize),
3065 : "", poWK->pProgress ) )
3066 : {
3067 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
3068 0 : eErr = CE_Failure;
3069 : }
3070 : }
3071 :
3072 : /* -------------------------------------------------------------------- */
3073 : /* Cleanup and return. */
3074 : /* -------------------------------------------------------------------- */
3075 10 : CPLFree( padfX );
3076 10 : CPLFree( padfY );
3077 10 : CPLFree( padfZ );
3078 10 : CPLFree( pabSuccess );
3079 10 : CPLFree( padfBSpline );
3080 :
3081 10 : return eErr;
3082 : }
3083 :
3084 : /************************************************************************/
3085 : /* GWKNearestByte() */
3086 : /* */
3087 : /* Case for 8bit input data with nearest neighbour resampling */
3088 : /* using valid flags. Should be as fast as possible for this */
3089 : /* particular transformation type. */
3090 : /************************************************************************/
3091 :
3092 10 : static CPLErr GWKNearestByte( GDALWarpKernel *poWK )
3093 :
3094 : {
3095 : int iDstY;
3096 10 : int nDstXSize = poWK->nDstXSize, nDstYSize = poWK->nDstYSize;
3097 10 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
3098 10 : CPLErr eErr = CE_None;
3099 :
3100 : CPLDebug( "GDAL", "GDALWarpKernel()::GWKNearestByte()\n"
3101 : "Src=%d,%d,%dx%d Dst=%d,%d,%dx%d",
3102 : poWK->nSrcXOff, poWK->nSrcYOff,
3103 : poWK->nSrcXSize, poWK->nSrcYSize,
3104 : poWK->nDstXOff, poWK->nDstYOff,
3105 10 : poWK->nDstXSize, poWK->nDstYSize );
3106 :
3107 10 : if( !poWK->pfnProgress( poWK->dfProgressBase, "", poWK->pProgress ) )
3108 : {
3109 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
3110 0 : return CE_Failure;
3111 : }
3112 :
3113 : /* -------------------------------------------------------------------- */
3114 : /* Allocate x,y,z coordinate arrays for transformation ... one */
3115 : /* scanlines worth of positions. */
3116 : /* -------------------------------------------------------------------- */
3117 : double *padfX, *padfY, *padfZ;
3118 : int *pabSuccess;
3119 :
3120 10 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3121 10 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3122 10 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3123 10 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
3124 :
3125 : /* ==================================================================== */
3126 : /* Loop over output lines. */
3127 : /* ==================================================================== */
3128 1121 : for( iDstY = 0; iDstY < nDstYSize && eErr == CE_None; iDstY++ )
3129 : {
3130 : int iDstX;
3131 :
3132 : /* -------------------------------------------------------------------- */
3133 : /* Setup points to transform to source image space. */
3134 : /* -------------------------------------------------------------------- */
3135 391732 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3136 : {
3137 390621 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
3138 390621 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
3139 390621 : padfZ[iDstX] = 0.0;
3140 : }
3141 :
3142 : /* -------------------------------------------------------------------- */
3143 : /* Transform the points from destination pixel/line coordinates */
3144 : /* to source pixel/line coordinates. */
3145 : /* -------------------------------------------------------------------- */
3146 : poWK->pfnTransformer( poWK->pTransformerArg, TRUE, nDstXSize,
3147 1111 : padfX, padfY, padfZ, pabSuccess );
3148 :
3149 : /* ==================================================================== */
3150 : /* Loop over pixels in output scanline. */
3151 : /* ==================================================================== */
3152 391732 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3153 : {
3154 : int iDstOffset;
3155 :
3156 390621 : if( !pabSuccess[iDstX] )
3157 0 : continue;
3158 :
3159 : /* -------------------------------------------------------------------- */
3160 : /* Figure out what pixel we want in our source raster, and skip */
3161 : /* further processing if it is well off the source image. */
3162 : /* -------------------------------------------------------------------- */
3163 : // We test against the value before casting to avoid the
3164 : // problem of asymmetric truncation effects around zero. That is
3165 : // -0.5 will be 0 when cast to an int.
3166 781218 : if( padfX[iDstX] < poWK->nSrcXOff
3167 390597 : || padfY[iDstX] < poWK->nSrcYOff )
3168 9740 : continue;
3169 :
3170 : int iSrcX, iSrcY, iSrcOffset;
3171 :
3172 380881 : iSrcX = ((int) padfX[iDstX]) - poWK->nSrcXOff;
3173 380881 : iSrcY = ((int) padfY[iDstX]) - poWK->nSrcYOff;
3174 :
3175 : // If operating outside natural projection area, padfX/Y can be
3176 : // a very huge positive number, that becomes -2147483648 in the
3177 : // int trucation. So it is necessary to test now for non negativeness.
3178 380881 : if( iSrcX < 0 || iSrcX >= nSrcXSize || iSrcY < 0 || iSrcY >= nSrcYSize )
3179 285079 : continue;
3180 :
3181 95802 : iSrcOffset = iSrcX + iSrcY * nSrcXSize;
3182 :
3183 : /* -------------------------------------------------------------------- */
3184 : /* Do not try to apply invalid source pixels to the dest. */
3185 : /* -------------------------------------------------------------------- */
3186 96203 : if( poWK->panUnifiedSrcValid != NULL
3187 401 : && !(poWK->panUnifiedSrcValid[iSrcOffset>>5]
3188 : & (0x01 << (iSrcOffset & 0x1f))) )
3189 100 : continue;
3190 :
3191 : /* -------------------------------------------------------------------- */
3192 : /* Do not try to apply transparent source pixels to the destination.*/
3193 : /* -------------------------------------------------------------------- */
3194 95702 : double dfDensity = 1.0;
3195 :
3196 95702 : if( poWK->pafUnifiedSrcDensity != NULL )
3197 : {
3198 90000 : dfDensity = poWK->pafUnifiedSrcDensity[iSrcOffset];
3199 90000 : if( dfDensity < 0.00001 )
3200 74897 : continue;
3201 : }
3202 :
3203 : /* ==================================================================== */
3204 : /* Loop processing each band. */
3205 : /* ==================================================================== */
3206 : int iBand;
3207 :
3208 20805 : iDstOffset = iDstX + iDstY * nDstXSize;
3209 :
3210 51610 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
3211 : {
3212 30805 : GByte bValue = 0;
3213 30805 : double dfBandDensity = 0.0;
3214 :
3215 : /* -------------------------------------------------------------------- */
3216 : /* Collect the source value. */
3217 : /* -------------------------------------------------------------------- */
3218 30805 : if ( GWKGetPixelByte( poWK, iBand, iSrcOffset, &dfBandDensity,
3219 : &bValue ) )
3220 : {
3221 30805 : if( dfBandDensity < 1.0 )
3222 : {
3223 1472 : if( dfBandDensity == 0.0 )
3224 : /* do nothing */;
3225 : else
3226 : {
3227 : /* let the general code take care of mixing */
3228 : GWKSetPixelValue( poWK, iBand, iDstOffset,
3229 : dfBandDensity, (double) bValue,
3230 1472 : 0.0 );
3231 : }
3232 : }
3233 : else
3234 : {
3235 29333 : poWK->papabyDstImage[iBand][iDstOffset] = bValue;
3236 : }
3237 : }
3238 : }
3239 :
3240 : /* -------------------------------------------------------------------- */
3241 : /* Mark this pixel valid/opaque in the output. */
3242 : /* -------------------------------------------------------------------- */
3243 20805 : GWKOverlayDensity( poWK, iDstOffset, dfDensity );
3244 :
3245 20805 : if( poWK->panDstValid != NULL )
3246 : {
3247 : poWK->panDstValid[iDstOffset>>5] |=
3248 702 : 0x01 << (iDstOffset & 0x1f);
3249 : }
3250 : } /* Next iDstX */
3251 :
3252 : /* -------------------------------------------------------------------- */
3253 : /* Report progress to the user, and optionally cancel out. */
3254 : /* -------------------------------------------------------------------- */
3255 1111 : if( !poWK->pfnProgress( poWK->dfProgressBase + poWK->dfProgressScale *
3256 : ((iDstY+1) / (double) nDstYSize),
3257 : "", poWK->pProgress ) )
3258 : {
3259 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
3260 0 : eErr = CE_Failure;
3261 : }
3262 : } /* Next iDstY */
3263 :
3264 : /* -------------------------------------------------------------------- */
3265 : /* Cleanup and return. */
3266 : /* -------------------------------------------------------------------- */
3267 10 : CPLFree( padfX );
3268 10 : CPLFree( padfY );
3269 10 : CPLFree( padfZ );
3270 10 : CPLFree( pabSuccess );
3271 :
3272 10 : return eErr;
3273 : }
3274 :
3275 : /************************************************************************/
3276 : /* GWKNearestNoMasksShort() */
3277 : /* */
3278 : /* Case for 16bit signed and unsigned integer input data with */
3279 : /* nearest neighbour resampling without concerning about masking. */
3280 : /* Should be as fast as possible for this particular */
3281 : /* transformation type. */
3282 : /************************************************************************/
3283 :
3284 14 : static CPLErr GWKNearestNoMasksShort( GDALWarpKernel *poWK )
3285 :
3286 : {
3287 : int iDstY;
3288 14 : int nDstXSize = poWK->nDstXSize, nDstYSize = poWK->nDstYSize;
3289 14 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
3290 14 : CPLErr eErr = CE_None;
3291 :
3292 : CPLDebug( "GDAL", "GDALWarpKernel()::GWKNearestNoMasksShort()\n"
3293 : "Src=%d,%d,%dx%d Dst=%d,%d,%dx%d",
3294 : poWK->nSrcXOff, poWK->nSrcYOff,
3295 : poWK->nSrcXSize, poWK->nSrcYSize,
3296 : poWK->nDstXOff, poWK->nDstYOff,
3297 14 : poWK->nDstXSize, poWK->nDstYSize );
3298 :
3299 14 : if( !poWK->pfnProgress( poWK->dfProgressBase, "", poWK->pProgress ) )
3300 : {
3301 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
3302 0 : return CE_Failure;
3303 : }
3304 :
3305 : /* -------------------------------------------------------------------- */
3306 : /* Allocate x,y,z coordinate arrays for transformation ... one */
3307 : /* scanlines worth of positions. */
3308 : /* -------------------------------------------------------------------- */
3309 : double *padfX, *padfY, *padfZ;
3310 : int *pabSuccess;
3311 :
3312 14 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3313 14 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3314 14 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3315 14 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
3316 :
3317 : /* ==================================================================== */
3318 : /* Loop over output lines. */
3319 : /* ==================================================================== */
3320 700 : for( iDstY = 0; iDstY < nDstYSize && eErr == CE_None; iDstY++ )
3321 : {
3322 : int iDstX;
3323 :
3324 : /* -------------------------------------------------------------------- */
3325 : /* Setup points to transform to source image space. */
3326 : /* -------------------------------------------------------------------- */
3327 328892 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3328 : {
3329 328206 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
3330 328206 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
3331 328206 : padfZ[iDstX] = 0.0;
3332 : }
3333 :
3334 : /* -------------------------------------------------------------------- */
3335 : /* Transform the points from destination pixel/line coordinates */
3336 : /* to source pixel/line coordinates. */
3337 : /* -------------------------------------------------------------------- */
3338 : poWK->pfnTransformer( poWK->pTransformerArg, TRUE, nDstXSize,
3339 686 : padfX, padfY, padfZ, pabSuccess );
3340 :
3341 : /* ==================================================================== */
3342 : /* Loop over pixels in output scanline. */
3343 : /* ==================================================================== */
3344 328892 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3345 : {
3346 : int iDstOffset;
3347 :
3348 328206 : if( !pabSuccess[iDstX] )
3349 63072 : continue;
3350 :
3351 : /* -------------------------------------------------------------------- */
3352 : /* Figure out what pixel we want in our source raster, and skip */
3353 : /* further processing if it is well off the source image. */
3354 : /* -------------------------------------------------------------------- */
3355 : // We test against the value before casting to avoid the
3356 : // problem of asymmetric truncation effects around zero. That is
3357 : // -0.5 will be 0 when cast to an int.
3358 530082 : if( padfX[iDstX] < poWK->nSrcXOff
3359 264948 : || padfY[iDstX] < poWK->nSrcYOff )
3360 186 : continue;
3361 :
3362 : int iSrcX, iSrcY, iSrcOffset;
3363 :
3364 264948 : iSrcX = ((int) padfX[iDstX]) - poWK->nSrcXOff;
3365 264948 : iSrcY = ((int) padfY[iDstX]) - poWK->nSrcYOff;
3366 :
3367 : // If operating outside natural projection area, padfX/Y can be
3368 : // a very huge positive number, that becomes -2147483648 in the
3369 : // int trucation. So it is necessary to test now for non negativeness.
3370 264948 : if( iSrcX < 0 || iSrcX >= nSrcXSize || iSrcY < 0 || iSrcY >= nSrcYSize )
3371 0 : continue;
3372 :
3373 264948 : iSrcOffset = iSrcX + iSrcY * nSrcXSize;
3374 :
3375 : /* ==================================================================== */
3376 : /* Loop processing each band. */
3377 : /* ==================================================================== */
3378 : int iBand;
3379 :
3380 264948 : iDstOffset = iDstX + iDstY * nDstXSize;
3381 :
3382 529896 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
3383 : {
3384 264948 : ((GInt16 *)poWK->papabyDstImage[iBand])[iDstOffset] =
3385 264948 : ((GInt16 *)poWK->papabySrcImage[iBand])[iSrcOffset];
3386 : }
3387 : }
3388 :
3389 : /* -------------------------------------------------------------------- */
3390 : /* Report progress to the user, and optionally cancel out. */
3391 : /* -------------------------------------------------------------------- */
3392 686 : if( !poWK->pfnProgress( poWK->dfProgressBase + poWK->dfProgressScale *
3393 : ((iDstY+1) / (double) nDstYSize),
3394 : "", poWK->pProgress ) )
3395 : {
3396 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
3397 0 : eErr = CE_Failure;
3398 : }
3399 : }
3400 :
3401 : /* -------------------------------------------------------------------- */
3402 : /* Cleanup and return. */
3403 : /* -------------------------------------------------------------------- */
3404 14 : CPLFree( padfX );
3405 14 : CPLFree( padfY );
3406 14 : CPLFree( padfZ );
3407 14 : CPLFree( pabSuccess );
3408 :
3409 14 : return eErr;
3410 : }
3411 :
3412 : /************************************************************************/
3413 : /* GWKBilinearNoMasksShort() */
3414 : /* */
3415 : /* Case for 16bit input data with cubic resampling without */
3416 : /* concerning about masking. Should be as fast as possible */
3417 : /* for this particular transformation type. */
3418 : /************************************************************************/
3419 :
3420 19 : static CPLErr GWKBilinearNoMasksShort( GDALWarpKernel *poWK )
3421 :
3422 : {
3423 : int iDstY;
3424 19 : int nDstXSize = poWK->nDstXSize, nDstYSize = poWK->nDstYSize;
3425 19 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
3426 19 : CPLErr eErr = CE_None;
3427 :
3428 : CPLDebug( "GDAL", "GDALWarpKernel()::GWKBilinearNoMasksShort()\n"
3429 : "Src=%d,%d,%dx%d Dst=%d,%d,%dx%d",
3430 : poWK->nSrcXOff, poWK->nSrcYOff,
3431 : poWK->nSrcXSize, poWK->nSrcYSize,
3432 : poWK->nDstXOff, poWK->nDstYOff,
3433 19 : poWK->nDstXSize, poWK->nDstYSize );
3434 :
3435 19 : if( !poWK->pfnProgress( poWK->dfProgressBase, "", poWK->pProgress ) )
3436 : {
3437 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
3438 0 : return CE_Failure;
3439 : }
3440 :
3441 : /* -------------------------------------------------------------------- */
3442 : /* Allocate x,y,z coordinate arrays for transformation ... one */
3443 : /* scanlines worth of positions. */
3444 : /* -------------------------------------------------------------------- */
3445 : double *padfX, *padfY, *padfZ;
3446 : int *pabSuccess;
3447 :
3448 19 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3449 19 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3450 19 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3451 19 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
3452 :
3453 : /* ==================================================================== */
3454 : /* Loop over output lines. */
3455 : /* ==================================================================== */
3456 654 : for( iDstY = 0; iDstY < nDstYSize && eErr == CE_None; iDstY++ )
3457 : {
3458 : int iDstX;
3459 :
3460 : /* -------------------------------------------------------------------- */
3461 : /* Setup points to transform to source image space. */
3462 : /* -------------------------------------------------------------------- */
3463 263942 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3464 : {
3465 263307 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
3466 263307 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
3467 263307 : padfZ[iDstX] = 0.0;
3468 : }
3469 :
3470 : /* -------------------------------------------------------------------- */
3471 : /* Transform the points from destination pixel/line coordinates */
3472 : /* to source pixel/line coordinates. */
3473 : /* -------------------------------------------------------------------- */
3474 : poWK->pfnTransformer( poWK->pTransformerArg, TRUE, nDstXSize,
3475 635 : padfX, padfY, padfZ, pabSuccess );
3476 :
3477 : /* ==================================================================== */
3478 : /* Loop over pixels in output scanline. */
3479 : /* ==================================================================== */
3480 263942 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3481 : {
3482 263307 : if( !pabSuccess[iDstX] )
3483 0 : continue;
3484 :
3485 : /* -------------------------------------------------------------------- */
3486 : /* Figure out what pixel we want in our source raster, and skip */
3487 : /* further processing if it is well off the source image. */
3488 : /* -------------------------------------------------------------------- */
3489 : // We test against the value before casting to avoid the
3490 : // problem of asymmetric truncation effects around zero. That is
3491 : // -0.5 will be 0 when cast to an int.
3492 526614 : if( padfX[iDstX] < poWK->nSrcXOff
3493 263307 : || padfY[iDstX] < poWK->nSrcYOff )
3494 0 : continue;
3495 :
3496 : int iSrcX, iSrcY, iSrcOffset;
3497 :
3498 263307 : iSrcX = ((int) padfX[iDstX]) - poWK->nSrcXOff;
3499 263307 : iSrcY = ((int) padfY[iDstX]) - poWK->nSrcYOff;
3500 :
3501 : // If operating outside natural projection area, padfX/Y can be
3502 : // a very huge positive number, that becomes -2147483648 in the
3503 : // int trucation. So it is necessary to test now for non negativeness.
3504 263307 : if( iSrcX < 0 || iSrcX >= nSrcXSize || iSrcY < 0 || iSrcY >= nSrcYSize )
3505 0 : continue;
3506 :
3507 263307 : iSrcOffset = iSrcX + iSrcY * nSrcXSize;
3508 :
3509 : /* ==================================================================== */
3510 : /* Loop processing each band. */
3511 : /* ==================================================================== */
3512 : int iBand;
3513 : int iDstOffset;
3514 :
3515 263307 : iDstOffset = iDstX + iDstY * nDstXSize;
3516 :
3517 526614 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
3518 : {
3519 263307 : GInt16 iValue = 0;
3520 : GWKBilinearResampleNoMasksShort( poWK, iBand,
3521 263307 : padfX[iDstX]-poWK->nSrcXOff,
3522 263307 : padfY[iDstX]-poWK->nSrcYOff,
3523 526614 : &iValue );
3524 263307 : ((GInt16 *)poWK->papabyDstImage[iBand])[iDstOffset] = iValue;
3525 : }
3526 : }
3527 :
3528 : /* -------------------------------------------------------------------- */
3529 : /* Report progress to the user, and optionally cancel out. */
3530 : /* -------------------------------------------------------------------- */
3531 635 : if( !poWK->pfnProgress( poWK->dfProgressBase + poWK->dfProgressScale *
3532 : ((iDstY+1) / (double) nDstYSize),
3533 : "", poWK->pProgress ) )
3534 : {
3535 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
3536 0 : eErr = CE_Failure;
3537 : }
3538 : }
3539 :
3540 : /* -------------------------------------------------------------------- */
3541 : /* Cleanup and return. */
3542 : /* -------------------------------------------------------------------- */
3543 19 : CPLFree( padfX );
3544 19 : CPLFree( padfY );
3545 19 : CPLFree( padfZ );
3546 19 : CPLFree( pabSuccess );
3547 :
3548 19 : return eErr;
3549 : }
3550 :
3551 : /************************************************************************/
3552 : /* GWKCubicNoMasksShort() */
3553 : /* */
3554 : /* Case for 16bit input data with cubic resampling without */
3555 : /* concerning about masking. Should be as fast as possible */
3556 : /* for this particular transformation type. */
3557 : /************************************************************************/
3558 :
3559 8 : static CPLErr GWKCubicNoMasksShort( GDALWarpKernel *poWK )
3560 :
3561 : {
3562 : int iDstY;
3563 8 : int nDstXSize = poWK->nDstXSize, nDstYSize = poWK->nDstYSize;
3564 8 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
3565 8 : CPLErr eErr = CE_None;
3566 :
3567 : CPLDebug( "GDAL", "GDALWarpKernel()::GWKCubicNoMasksShort()\n"
3568 : "Src=%d,%d,%dx%d Dst=%d,%d,%dx%d",
3569 : poWK->nSrcXOff, poWK->nSrcYOff,
3570 : poWK->nSrcXSize, poWK->nSrcYSize,
3571 : poWK->nDstXOff, poWK->nDstYOff,
3572 8 : poWK->nDstXSize, poWK->nDstYSize );
3573 :
3574 8 : if( !poWK->pfnProgress( poWK->dfProgressBase, "", poWK->pProgress ) )
3575 : {
3576 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
3577 0 : return CE_Failure;
3578 : }
3579 :
3580 : /* -------------------------------------------------------------------- */
3581 : /* Allocate x,y,z coordinate arrays for transformation ... one */
3582 : /* scanlines worth of positions. */
3583 : /* -------------------------------------------------------------------- */
3584 : double *padfX, *padfY, *padfZ;
3585 : int *pabSuccess;
3586 :
3587 8 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3588 8 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3589 8 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3590 8 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
3591 :
3592 : /* ==================================================================== */
3593 : /* Loop over output lines. */
3594 : /* ==================================================================== */
3595 533 : for( iDstY = 0; iDstY < nDstYSize && eErr == CE_None; iDstY++ )
3596 : {
3597 : int iDstX;
3598 :
3599 : /* -------------------------------------------------------------------- */
3600 : /* Setup points to transform to source image space. */
3601 : /* -------------------------------------------------------------------- */
3602 262732 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3603 : {
3604 262207 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
3605 262207 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
3606 262207 : padfZ[iDstX] = 0.0;
3607 : }
3608 :
3609 : /* -------------------------------------------------------------------- */
3610 : /* Transform the points from destination pixel/line coordinates */
3611 : /* to source pixel/line coordinates. */
3612 : /* -------------------------------------------------------------------- */
3613 : poWK->pfnTransformer( poWK->pTransformerArg, TRUE, nDstXSize,
3614 525 : padfX, padfY, padfZ, pabSuccess );
3615 :
3616 : /* ==================================================================== */
3617 : /* Loop over pixels in output scanline. */
3618 : /* ==================================================================== */
3619 262732 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3620 : {
3621 262207 : if( !pabSuccess[iDstX] )
3622 0 : continue;
3623 :
3624 : /* -------------------------------------------------------------------- */
3625 : /* Figure out what pixel we want in our source raster, and skip */
3626 : /* further processing if it is well off the source image. */
3627 : /* -------------------------------------------------------------------- */
3628 : // We test against the value before casting to avoid the
3629 : // problem of asymmetric truncation effects around zero. That is
3630 : // -0.5 will be 0 when cast to an int.
3631 524414 : if( padfX[iDstX] < poWK->nSrcXOff
3632 262207 : || padfY[iDstX] < poWK->nSrcYOff )
3633 0 : continue;
3634 :
3635 : int iSrcX, iSrcY, iSrcOffset;
3636 :
3637 262207 : iSrcX = ((int) padfX[iDstX]) - poWK->nSrcXOff;
3638 262207 : iSrcY = ((int) padfY[iDstX]) - poWK->nSrcYOff;
3639 :
3640 : // If operating outside natural projection area, padfX/Y can be
3641 : // a very huge positive number, that becomes -2147483648 in the
3642 : // int trucation. So it is necessary to test now for non negativeness.
3643 262207 : if( iSrcX < 0 || iSrcX >= nSrcXSize || iSrcY < 0 || iSrcY >= nSrcYSize )
3644 0 : continue;
3645 :
3646 262207 : iSrcOffset = iSrcX + iSrcY * nSrcXSize;
3647 :
3648 : /* ==================================================================== */
3649 : /* Loop processing each band. */
3650 : /* ==================================================================== */
3651 : int iBand;
3652 : int iDstOffset;
3653 :
3654 262207 : iDstOffset = iDstX + iDstY * nDstXSize;
3655 :
3656 524414 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
3657 : {
3658 262207 : GInt16 iValue = 0;
3659 : GWKCubicResampleNoMasksShort( poWK, iBand,
3660 262207 : padfX[iDstX]-poWK->nSrcXOff,
3661 262207 : padfY[iDstX]-poWK->nSrcYOff,
3662 524414 : &iValue );
3663 262207 : ((GInt16 *)poWK->papabyDstImage[iBand])[iDstOffset] = iValue;
3664 : }
3665 : }
3666 :
3667 : /* -------------------------------------------------------------------- */
3668 : /* Report progress to the user, and optionally cancel out. */
3669 : /* -------------------------------------------------------------------- */
3670 525 : if( !poWK->pfnProgress( poWK->dfProgressBase + poWK->dfProgressScale *
3671 : ((iDstY+1) / (double) nDstYSize),
3672 : "", poWK->pProgress ) )
3673 : {
3674 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
3675 0 : eErr = CE_Failure;
3676 : }
3677 : }
3678 :
3679 : /* -------------------------------------------------------------------- */
3680 : /* Cleanup and return. */
3681 : /* -------------------------------------------------------------------- */
3682 8 : CPLFree( padfX );
3683 8 : CPLFree( padfY );
3684 8 : CPLFree( padfZ );
3685 8 : CPLFree( pabSuccess );
3686 :
3687 8 : return eErr;
3688 : }
3689 :
3690 : /************************************************************************/
3691 : /* GWKCubicSplineNoMasksShort() */
3692 : /* */
3693 : /* Case for 16bit input data with cubic resampling without */
3694 : /* concerning about masking. Should be as fast as possible */
3695 : /* for this particular transformation type. */
3696 : /************************************************************************/
3697 :
3698 8 : static CPLErr GWKCubicSplineNoMasksShort( GDALWarpKernel *poWK )
3699 :
3700 : {
3701 : int iDstY;
3702 8 : int nDstXSize = poWK->nDstXSize, nDstYSize = poWK->nDstYSize;
3703 8 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
3704 8 : CPLErr eErr = CE_None;
3705 :
3706 : CPLDebug( "GDAL", "GDALWarpKernel()::GWKCubicSplineNoMasksShort()\n"
3707 : "Src=%d,%d,%dx%d Dst=%d,%d,%dx%d",
3708 : poWK->nSrcXOff, poWK->nSrcYOff,
3709 : poWK->nSrcXSize, poWK->nSrcYSize,
3710 : poWK->nDstXOff, poWK->nDstYOff,
3711 8 : poWK->nDstXSize, poWK->nDstYSize );
3712 :
3713 8 : if( !poWK->pfnProgress( poWK->dfProgressBase, "", poWK->pProgress ) )
3714 : {
3715 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
3716 0 : return CE_Failure;
3717 : }
3718 :
3719 : /* -------------------------------------------------------------------- */
3720 : /* Allocate x,y,z coordinate arrays for transformation ... one */
3721 : /* scanlines worth of positions. */
3722 : /* -------------------------------------------------------------------- */
3723 : double *padfX, *padfY, *padfZ;
3724 : int *pabSuccess;
3725 :
3726 8 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3727 8 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3728 8 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3729 8 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
3730 :
3731 8 : int nXRadius = poWK->nXRadius;
3732 : // Make space to save weights
3733 8 : double *padfBSpline = (double *)CPLCalloc( nXRadius * 2, sizeof(double) );
3734 :
3735 : /* ==================================================================== */
3736 : /* Loop over output lines. */
3737 : /* ==================================================================== */
3738 533 : for( iDstY = 0; iDstY < nDstYSize && eErr == CE_None; iDstY++ )
3739 : {
3740 : int iDstX;
3741 :
3742 : /* -------------------------------------------------------------------- */
3743 : /* Setup points to transform to source image space. */
3744 : /* -------------------------------------------------------------------- */
3745 262732 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3746 : {
3747 262207 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
3748 262207 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
3749 262207 : padfZ[iDstX] = 0.0;
3750 : }
3751 :
3752 : /* -------------------------------------------------------------------- */
3753 : /* Transform the points from destination pixel/line coordinates */
3754 : /* to source pixel/line coordinates. */
3755 : /* -------------------------------------------------------------------- */
3756 : poWK->pfnTransformer( poWK->pTransformerArg, TRUE, nDstXSize,
3757 525 : padfX, padfY, padfZ, pabSuccess );
3758 :
3759 : /* ==================================================================== */
3760 : /* Loop over pixels in output scanline. */
3761 : /* ==================================================================== */
3762 262732 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3763 : {
3764 262207 : if( !pabSuccess[iDstX] )
3765 0 : continue;
3766 :
3767 : /* -------------------------------------------------------------------- */
3768 : /* Figure out what pixel we want in our source raster, and skip */
3769 : /* further processing if it is well off the source image. */
3770 : /* -------------------------------------------------------------------- */
3771 : // We test against the value before casting to avoid the
3772 : // problem of asymmetric truncation effects around zero. That is
3773 : // -0.5 will be 0 when cast to an int.
3774 524414 : if( padfX[iDstX] < poWK->nSrcXOff
3775 262207 : || padfY[iDstX] < poWK->nSrcYOff )
3776 0 : continue;
3777 :
3778 : int iSrcX, iSrcY, iSrcOffset;
3779 :
3780 262207 : iSrcX = ((int) padfX[iDstX]) - poWK->nSrcXOff;
3781 262207 : iSrcY = ((int) padfY[iDstX]) - poWK->nSrcYOff;
3782 :
3783 : // If operating outside natural projection area, padfX/Y can be
3784 : // a very huge positive number, that becomes -2147483648 in the
3785 : // int trucation. So it is necessary to test now for non negativeness.
3786 262207 : if( iSrcX < 0 || iSrcX >= nSrcXSize || iSrcY < 0 || iSrcY >= nSrcYSize )
3787 0 : continue;
3788 :
3789 262207 : iSrcOffset = iSrcX + iSrcY * nSrcXSize;
3790 :
3791 : /* ==================================================================== */
3792 : /* Loop processing each band. */
3793 : /* ==================================================================== */
3794 : int iBand;
3795 : int iDstOffset;
3796 :
3797 262207 : iDstOffset = iDstX + iDstY * nDstXSize;
3798 :
3799 524414 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
3800 : {
3801 262207 : GInt16 iValue = 0;
3802 : GWKCubicSplineResampleNoMasksShort( poWK, iBand,
3803 262207 : padfX[iDstX]-poWK->nSrcXOff,
3804 262207 : padfY[iDstX]-poWK->nSrcYOff,
3805 : &iValue,
3806 524414 : padfBSpline);
3807 262207 : ((GInt16 *)poWK->papabyDstImage[iBand])[iDstOffset] = iValue;
3808 : }
3809 : }
3810 :
3811 : /* -------------------------------------------------------------------- */
3812 : /* Report progress to the user, and optionally cancel out. */
3813 : /* -------------------------------------------------------------------- */
3814 525 : if( !poWK->pfnProgress( poWK->dfProgressBase + poWK->dfProgressScale *
3815 : ((iDstY+1) / (double) nDstYSize),
3816 : "", poWK->pProgress ) )
3817 : {
3818 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
3819 0 : eErr = CE_Failure;
3820 : }
3821 : }
3822 :
3823 : /* -------------------------------------------------------------------- */
3824 : /* Cleanup and return. */
3825 : /* -------------------------------------------------------------------- */
3826 8 : CPLFree( padfX );
3827 8 : CPLFree( padfY );
3828 8 : CPLFree( padfZ );
3829 8 : CPLFree( pabSuccess );
3830 8 : CPLFree( padfBSpline );
3831 :
3832 8 : return eErr;
3833 : }
3834 :
3835 : /************************************************************************/
3836 : /* GWKNearestShort() */
3837 : /* */
3838 : /* Case for 32bit float input data with nearest neighbour */
3839 : /* resampling using valid flags. Should be as fast as possible */
3840 : /* for this particular transformation type. */
3841 : /************************************************************************/
3842 :
3843 7 : static CPLErr GWKNearestShort( GDALWarpKernel *poWK )
3844 :
3845 : {
3846 : int iDstY;
3847 7 : int nDstXSize = poWK->nDstXSize, nDstYSize = poWK->nDstYSize;
3848 7 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
3849 7 : CPLErr eErr = CE_None;
3850 :
3851 : CPLDebug( "GDAL", "GDALWarpKernel()::GWKNearestShort()\n"
3852 : "Src=%d,%d,%dx%d Dst=%d,%d,%dx%d",
3853 : poWK->nSrcXOff, poWK->nSrcYOff,
3854 : poWK->nSrcXSize, poWK->nSrcYSize,
3855 : poWK->nDstXOff, poWK->nDstYOff,
3856 7 : poWK->nDstXSize, poWK->nDstYSize );
3857 :
3858 7 : if( !poWK->pfnProgress( poWK->dfProgressBase, "", poWK->pProgress ) )
3859 : {
3860 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
3861 0 : return CE_Failure;
3862 : }
3863 :
3864 : /* -------------------------------------------------------------------- */
3865 : /* Allocate x,y,z coordinate arrays for transformation ... one */
3866 : /* scanlines worth of positions. */
3867 : /* -------------------------------------------------------------------- */
3868 : double *padfX, *padfY, *padfZ;
3869 : int *pabSuccess;
3870 :
3871 7 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3872 7 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3873 7 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
3874 7 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
3875 :
3876 : /* ==================================================================== */
3877 : /* Loop over output lines. */
3878 : /* ==================================================================== */
3879 2297 : for( iDstY = 0; iDstY < nDstYSize && eErr == CE_None; iDstY++ )
3880 : {
3881 : int iDstX;
3882 :
3883 : /* -------------------------------------------------------------------- */
3884 : /* Setup points to transform to source image space. */
3885 : /* -------------------------------------------------------------------- */
3886 1620626 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3887 : {
3888 1618336 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
3889 1618336 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
3890 1618336 : padfZ[iDstX] = 0.0;
3891 : }
3892 :
3893 : /* -------------------------------------------------------------------- */
3894 : /* Transform the points from destination pixel/line coordinates */
3895 : /* to source pixel/line coordinates. */
3896 : /* -------------------------------------------------------------------- */
3897 : poWK->pfnTransformer( poWK->pTransformerArg, TRUE, nDstXSize,
3898 2290 : padfX, padfY, padfZ, pabSuccess );
3899 :
3900 : /* ==================================================================== */
3901 : /* Loop over pixels in output scanline. */
3902 : /* ==================================================================== */
3903 1620626 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
3904 : {
3905 : int iDstOffset;
3906 :
3907 1618336 : if( !pabSuccess[iDstX] )
3908 0 : continue;
3909 :
3910 : /* -------------------------------------------------------------------- */
3911 : /* Figure out what pixel we want in our source raster, and skip */
3912 : /* further processing if it is well off the source image. */
3913 : /* -------------------------------------------------------------------- */
3914 : // We test against the value before casting to avoid the
3915 : // problem of asymmetric truncation effects around zero. That is
3916 : // -0.5 will be 0 when cast to an int.
3917 3236646 : if( padfX[iDstX] < poWK->nSrcXOff
3918 1618310 : || padfY[iDstX] < poWK->nSrcYOff )
3919 19528 : continue;
3920 :
3921 : int iSrcX, iSrcY, iSrcOffset;
3922 :
3923 1598808 : iSrcX = ((int) padfX[iDstX]) - poWK->nSrcXOff;
3924 1598808 : iSrcY = ((int) padfY[iDstX]) - poWK->nSrcYOff;
3925 :
3926 : // If operating outside natural projection area, padfX/Y can be
3927 : // a very huge positive number, that becomes -2147483648 in the
3928 : // int trucation. So it is necessary to test now for non negativeness.
3929 1598808 : if( iSrcX < 0 || iSrcX >= nSrcXSize || iSrcY < 0 || iSrcY >= nSrcYSize )
3930 110742 : continue;
3931 :
3932 1488066 : iSrcOffset = iSrcX + iSrcY * nSrcXSize;
3933 :
3934 : /* -------------------------------------------------------------------- */
3935 : /* Don't generate output pixels for which the source valid */
3936 : /* mask exists and is invalid. */
3937 : /* -------------------------------------------------------------------- */
3938 1488066 : if( poWK->panUnifiedSrcValid != NULL
3939 0 : && !(poWK->panUnifiedSrcValid[iSrcOffset>>5]
3940 : & (0x01 << (iSrcOffset & 0x1f))) )
3941 0 : continue;
3942 :
3943 : /* -------------------------------------------------------------------- */
3944 : /* Do not try to apply transparent source pixels to the destination.*/
3945 : /* -------------------------------------------------------------------- */
3946 1488066 : double dfDensity = 1.0;
3947 :
3948 1488066 : if( poWK->pafUnifiedSrcDensity != NULL )
3949 : {
3950 802 : dfDensity = poWK->pafUnifiedSrcDensity[iSrcOffset];
3951 802 : if( dfDensity < 0.00001 )
3952 401 : continue;
3953 : }
3954 :
3955 : /* ==================================================================== */
3956 : /* Loop processing each band. */
3957 : /* ==================================================================== */
3958 : int iBand;
3959 1487665 : iDstOffset = iDstX + iDstY * nDstXSize;
3960 :
3961 2976132 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
3962 : {
3963 1488467 : GInt16 iValue = 0;
3964 1488467 : double dfBandDensity = 0.0;
3965 :
3966 : /* -------------------------------------------------------------------- */
3967 : /* Collect the source value. */
3968 : /* -------------------------------------------------------------------- */
3969 1488467 : if ( GWKGetPixelShort( poWK, iBand, iSrcOffset, &dfBandDensity,
3970 : &iValue ) )
3971 : {
3972 1487752 : if( dfBandDensity < 1.0 )
3973 : {
3974 0 : if( dfBandDensity == 0.0 )
3975 : /* do nothing */;
3976 : else
3977 : {
3978 : /* let the general code take care of mixing */
3979 : GWKSetPixelValue( poWK, iBand, iDstOffset,
3980 : dfBandDensity, (double) iValue,
3981 0 : 0.0 );
3982 : }
3983 : }
3984 : else
3985 : {
3986 1487752 : ((GInt16 *)poWK->papabyDstImage[iBand])[iDstOffset] = iValue;
3987 : }
3988 : }
3989 : }
3990 :
3991 : /* -------------------------------------------------------------------- */
3992 : /* Mark this pixel valid/opaque in the output. */
3993 : /* -------------------------------------------------------------------- */
3994 1487665 : GWKOverlayDensity( poWK, iDstOffset, dfDensity );
3995 :
3996 1487665 : if( poWK->panDstValid != NULL )
3997 : {
3998 : poWK->panDstValid[iDstOffset>>5] |=
3999 0 : 0x01 << (iDstOffset & 0x1f);
4000 : }
4001 : } /* Next iDstX */
4002 :
4003 : /* -------------------------------------------------------------------- */
4004 : /* Report progress to the user, and optionally cancel out. */
4005 : /* -------------------------------------------------------------------- */
4006 2290 : if( !poWK->pfnProgress( poWK->dfProgressBase + poWK->dfProgressScale *
4007 : ((iDstY+1) / (double) nDstYSize),
4008 : "", poWK->pProgress ) )
4009 : {
4010 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
4011 0 : eErr = CE_Failure;
4012 : }
4013 : } /* Next iDstY */
4014 :
4015 : /* -------------------------------------------------------------------- */
4016 : /* Cleanup and return. */
4017 : /* -------------------------------------------------------------------- */
4018 7 : CPLFree( padfX );
4019 7 : CPLFree( padfY );
4020 7 : CPLFree( padfZ );
4021 7 : CPLFree( pabSuccess );
4022 :
4023 7 : return eErr;
4024 : }
4025 :
4026 : /************************************************************************/
4027 : /* GWKNearestNoMasksFloat() */
4028 : /* */
4029 : /* Case for 32bit float input data with nearest neighbour */
4030 : /* resampling without concerning about masking. Should be as fast */
4031 : /* as possible for this particular transformation type. */
4032 : /************************************************************************/
4033 :
4034 8 : static CPLErr GWKNearestNoMasksFloat( GDALWarpKernel *poWK )
4035 :
4036 : {
4037 : int iDstY;
4038 8 : int nDstXSize = poWK->nDstXSize, nDstYSize = poWK->nDstYSize;
4039 8 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
4040 8 : CPLErr eErr = CE_None;
4041 :
4042 : CPLDebug( "GDAL", "GDALWarpKernel()::GWKNearestNoMasksFloat()\n"
4043 : "Src=%d,%d,%dx%d Dst=%d,%d,%dx%d",
4044 : poWK->nSrcXOff, poWK->nSrcYOff,
4045 : poWK->nSrcXSize, poWK->nSrcYSize,
4046 : poWK->nDstXOff, poWK->nDstYOff,
4047 8 : poWK->nDstXSize, poWK->nDstYSize );
4048 :
4049 8 : if( !poWK->pfnProgress( poWK->dfProgressBase, "", poWK->pProgress ) )
4050 : {
4051 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
4052 0 : return CE_Failure;
4053 : }
4054 :
4055 : /* -------------------------------------------------------------------- */
4056 : /* Allocate x,y,z coordinate arrays for transformation ... one */
4057 : /* scanlines worth of positions. */
4058 : /* -------------------------------------------------------------------- */
4059 : double *padfX, *padfY, *padfZ;
4060 : int *pabSuccess;
4061 :
4062 8 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4063 8 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4064 8 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4065 8 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
4066 :
4067 : /* ==================================================================== */
4068 : /* Loop over output lines. */
4069 : /* ==================================================================== */
4070 533 : for( iDstY = 0; iDstY < nDstYSize && eErr == CE_None; iDstY++ )
4071 : {
4072 : int iDstX;
4073 :
4074 : /* -------------------------------------------------------------------- */
4075 : /* Setup points to transform to source image space. */
4076 : /* -------------------------------------------------------------------- */
4077 262732 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4078 : {
4079 262207 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
4080 262207 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
4081 262207 : padfZ[iDstX] = 0.0;
4082 : }
4083 :
4084 : /* -------------------------------------------------------------------- */
4085 : /* Transform the points from destination pixel/line coordinates */
4086 : /* to source pixel/line coordinates. */
4087 : /* -------------------------------------------------------------------- */
4088 : poWK->pfnTransformer( poWK->pTransformerArg, TRUE, nDstXSize,
4089 525 : padfX, padfY, padfZ, pabSuccess );
4090 :
4091 : /* ==================================================================== */
4092 : /* Loop over pixels in output scanline. */
4093 : /* ==================================================================== */
4094 262732 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4095 : {
4096 : int iDstOffset;
4097 :
4098 262207 : if( !pabSuccess[iDstX] )
4099 0 : continue;
4100 :
4101 : /* -------------------------------------------------------------------- */
4102 : /* Figure out what pixel we want in our source raster, and skip */
4103 : /* further processing if it is well off the source image. */
4104 : /* -------------------------------------------------------------------- */
4105 : // We test against the value before casting to avoid the
4106 : // problem of asymmetric truncation effects around zero. That is
4107 : // -0.5 will be 0 when cast to an int.
4108 524414 : if( padfX[iDstX] < poWK->nSrcXOff
4109 262207 : || padfY[iDstX] < poWK->nSrcYOff )
4110 0 : continue;
4111 :
4112 : int iSrcX, iSrcY, iSrcOffset;
4113 :
4114 262207 : iSrcX = ((int) padfX[iDstX]) - poWK->nSrcXOff;
4115 262207 : iSrcY = ((int) padfY[iDstX]) - poWK->nSrcYOff;
4116 :
4117 : // If operating outside natural projection area, padfX/Y can be
4118 : // a very huge positive number, that becomes -2147483648 in the
4119 : // int trucation. So it is necessary to test now for non negativeness.
4120 262207 : if( iSrcX < 0 || iSrcX >= nSrcXSize || iSrcY < 0 || iSrcY >= nSrcYSize )
4121 0 : continue;
4122 :
4123 262207 : iSrcOffset = iSrcX + iSrcY * nSrcXSize;
4124 :
4125 : /* ==================================================================== */
4126 : /* Loop processing each band. */
4127 : /* ==================================================================== */
4128 : int iBand;
4129 :
4130 262207 : iDstOffset = iDstX + iDstY * nDstXSize;
4131 :
4132 524414 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
4133 : {
4134 262207 : ((float *)poWK->papabyDstImage[iBand])[iDstOffset] =
4135 262207 : ((float *)poWK->papabySrcImage[iBand])[iSrcOffset];
4136 : }
4137 : }
4138 :
4139 : /* -------------------------------------------------------------------- */
4140 : /* Report progress to the user, and optionally cancel out. */
4141 : /* -------------------------------------------------------------------- */
4142 525 : if( !poWK->pfnProgress( poWK->dfProgressBase + poWK->dfProgressScale *
4143 : ((iDstY+1) / (double) nDstYSize),
4144 : "", poWK->pProgress ) )
4145 : {
4146 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
4147 0 : eErr = CE_Failure;
4148 : }
4149 : }
4150 :
4151 : /* -------------------------------------------------------------------- */
4152 : /* Cleanup and return. */
4153 : /* -------------------------------------------------------------------- */
4154 8 : CPLFree( padfX );
4155 8 : CPLFree( padfY );
4156 8 : CPLFree( padfZ );
4157 8 : CPLFree( pabSuccess );
4158 :
4159 8 : return eErr;
4160 : }
4161 :
4162 : /************************************************************************/
4163 : /* GWKNearestFloat() */
4164 : /* */
4165 : /* Case for 32bit float input data with nearest neighbour */
4166 : /* resampling using valid flags. Should be as fast as possible */
4167 : /* for this particular transformation type. */
4168 : /************************************************************************/
4169 :
4170 2 : static CPLErr GWKNearestFloat( GDALWarpKernel *poWK )
4171 :
4172 : {
4173 : int iDstY;
4174 2 : int nDstXSize = poWK->nDstXSize, nDstYSize = poWK->nDstYSize;
4175 2 : int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize;
4176 2 : CPLErr eErr = CE_None;
4177 :
4178 : CPLDebug( "GDAL", "GDALWarpKernel()::GWKNearestFloat()\n"
4179 : "Src=%d,%d,%dx%d Dst=%d,%d,%dx%d",
4180 : poWK->nSrcXOff, poWK->nSrcYOff,
4181 : poWK->nSrcXSize, poWK->nSrcYSize,
4182 : poWK->nDstXOff, poWK->nDstYOff,
4183 2 : poWK->nDstXSize, poWK->nDstYSize );
4184 :
4185 2 : if( !poWK->pfnProgress( poWK->dfProgressBase, "", poWK->pProgress ) )
4186 : {
4187 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
4188 0 : return CE_Failure;
4189 : }
4190 :
4191 : /* -------------------------------------------------------------------- */
4192 : /* Allocate x,y,z coordinate arrays for transformation ... one */
4193 : /* scanlines worth of positions. */
4194 : /* -------------------------------------------------------------------- */
4195 : double *padfX, *padfY, *padfZ;
4196 : int *pabSuccess;
4197 :
4198 2 : padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4199 2 : padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4200 2 : padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize);
4201 2 : pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize);
4202 :
4203 : /* ==================================================================== */
4204 : /* Loop over output lines. */
4205 : /* ==================================================================== */
4206 258 : for( iDstY = 0; iDstY < nDstYSize && eErr == CE_None; iDstY++ )
4207 : {
4208 : int iDstX;
4209 :
4210 : /* -------------------------------------------------------------------- */
4211 : /* Setup points to transform to source image space. */
4212 : /* -------------------------------------------------------------------- */
4213 131328 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4214 : {
4215 131072 : padfX[iDstX] = iDstX + 0.5 + poWK->nDstXOff;
4216 131072 : padfY[iDstX] = iDstY + 0.5 + poWK->nDstYOff;
4217 131072 : padfZ[iDstX] = 0.0;
4218 : }
4219 :
4220 : /* -------------------------------------------------------------------- */
4221 : /* Transform the points from destination pixel/line coordinates */
4222 : /* to source pixel/line coordinates. */
4223 : /* -------------------------------------------------------------------- */
4224 : poWK->pfnTransformer( poWK->pTransformerArg, TRUE, nDstXSize,
4225 256 : padfX, padfY, padfZ, pabSuccess );
4226 :
4227 : /* ==================================================================== */
4228 : /* Loop over pixels in output scanline. */
4229 : /* ==================================================================== */
4230 131328 : for( iDstX = 0; iDstX < nDstXSize; iDstX++ )
4231 : {
4232 : int iDstOffset;
4233 :
4234 131072 : if( !pabSuccess[iDstX] )
4235 0 : continue;
4236 :
4237 : /* -------------------------------------------------------------------- */
4238 : /* Figure out what pixel we want in our source raster, and skip */
4239 : /* further processing if it is well off the source image. */
4240 : /* -------------------------------------------------------------------- */
4241 : // We test against the value before casting to avoid the
4242 : // problem of asymmetric truncation effects around zero. That is
4243 : // -0.5 will be 0 when cast to an int.
4244 262118 : if( padfX[iDstX] < poWK->nSrcXOff
4245 131046 : || padfY[iDstX] < poWK->nSrcYOff )
4246 19528 : continue;
4247 :
4248 : int iSrcX, iSrcY, iSrcOffset;
4249 :
4250 111544 : iSrcX = ((int) padfX[iDstX]) - poWK->nSrcXOff;
4251 111544 : iSrcY = ((int) padfY[iDstX]) - poWK->nSrcYOff;
4252 :
4253 : // If operating outside natural projection area, padfX/Y can be
4254 : // a very huge positive number, that becomes -2147483648 in the
4255 : // int trucation. So it is necessary to test now for non negativeness.
4256 111544 : if( iSrcX < 0 || iSrcX >= nSrcXSize || iSrcY < 0 || iSrcY >= nSrcYSize )
4257 110742 : continue;
4258 :
4259 802 : iSrcOffset = iSrcX + iSrcY * nSrcXSize;
4260 :
4261 : /* -------------------------------------------------------------------- */
4262 : /* Do not try to apply invalid source pixels to the dest. */
4263 : /* -------------------------------------------------------------------- */
4264 802 : if( poWK->panUnifiedSrcValid != NULL
4265 0 : && !(poWK->panUnifiedSrcValid[iSrcOffset>>5]
4266 : & (0x01 << (iSrcOffset & 0x1f))) )
4267 0 : continue;
4268 :
4269 : /* -------------------------------------------------------------------- */
4270 : /* Do not try to apply transparent source pixels to the destination.*/
4271 : /* -------------------------------------------------------------------- */
4272 802 : double dfDensity = 1.0;
4273 :
4274 802 : if( poWK->pafUnifiedSrcDensity != NULL )
4275 : {
4276 802 : dfDensity = poWK->pafUnifiedSrcDensity[iSrcOffset];
4277 802 : if( dfDensity < 0.00001 )
4278 401 : continue;
4279 : }
4280 :
4281 : /* ==================================================================== */
4282 : /* Loop processing each band. */
4283 : /* ==================================================================== */
4284 : int iBand;
4285 :
4286 401 : iDstOffset = iDstX + iDstY * nDstXSize;
4287 :
4288 1604 : for( iBand = 0; iBand < poWK->nBands; iBand++ )
4289 : {
4290 1203 : float fValue = 0;
4291 1203 : double dfBandDensity = 0.0;
4292 :
4293 : /* -------------------------------------------------------------------- */
4294 : /* Collect the source value. */
4295 : /* -------------------------------------------------------------------- */
4296 1203 : if ( GWKGetPixelFloat( poWK, iBand, iSrcOffset, &dfBandDensity,
4297 : &fValue ) )
4298 : {
4299 1203 : if( dfBandDensity < 1.0 )
4300 : {
4301 0 : if( dfBandDensity == 0.0 )
4302 : /* do nothing */;
4303 : else
4304 : {
4305 : /* let the general code take care of mixing */
4306 : GWKSetPixelValue( poWK, iBand, iDstOffset,
4307 : dfBandDensity, (double) fValue,
4308 0 : 0.0 );
4309 : }
4310 : }
4311 : else
4312 : {
4313 1203 : ((float *)poWK->papabyDstImage[iBand])[iDstOffset]
4314 1203 : = fValue;
4315 : }
4316 : }
4317 : }
4318 :
4319 : /* -------------------------------------------------------------------- */
4320 : /* Mark this pixel valid/opaque in the output. */
4321 : /* -------------------------------------------------------------------- */
4322 401 : GWKOverlayDensity( poWK, iDstOffset, dfDensity );
4323 :
4324 401 : if( poWK->panDstValid != NULL )
4325 : {
4326 : poWK->panDstValid[iDstOffset>>5] |=
4327 0 : 0x01 << (iDstOffset & 0x1f);
4328 : }
4329 : }
4330 :
4331 : /* -------------------------------------------------------------------- */
4332 : /* Report progress to the user, and optionally cancel out. */
4333 : /* -------------------------------------------------------------------- */
4334 256 : if( !poWK->pfnProgress( poWK->dfProgressBase + poWK->dfProgressScale *
4335 : ((iDstY+1) / (double) nDstYSize),
4336 : "", poWK->pProgress ) )
4337 : {
4338 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
4339 0 : eErr = CE_Failure;
4340 : }
4341 : }
4342 :
4343 : /* -------------------------------------------------------------------- */
4344 : /* Cleanup and return. */
4345 : /* -------------------------------------------------------------------- */
4346 2 : CPLFree( padfX );
4347 2 : CPLFree( padfY );
4348 2 : CPLFree( padfZ );
4349 2 : CPLFree( pabSuccess );
4350 :
4351 2 : return eErr;
4352 : }
4353 :
|
