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