1 /******************************************************************************
2 * $Id: rasterio.cpp 29161 2015-05-06 10:18:19Z rouault $
3 *
4 * Project: GDAL Core
5 * Purpose: Contains default implementation of GDALRasterBand::IRasterIO()
6 * and supporting functions of broader utility.
7 * Author: Frank Warmerdam, warmerdam@pobox.com
8 *
9 ******************************************************************************
10 * Copyright (c) 1998, Frank Warmerdam
11 * Copyright (c) 2007-2014, Even Rouault <even dot rouault at mines-paris dot org>
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 "gdal_priv.h"
33 #include "gdal_vrt.h"
34 #include "vrtdataset.h"
35 #include "memdataset.h"
36 #include "gdalwarper.h"
37
38 // Define a list of "C++" compilers that have broken template support or
39 // broken scoping so we can fall back on the legacy implementation of
40 // GDALCopyWords
41 #define NOT_BROKEN_COMPILER \
42 (!(defined(_MSC_VER) && _MSC_VER <= 1200) && !defined(__BORLANDC__) && \
43 !defined(__SUNPRO_CC))
44
45 #if NOT_BROKEN_COMPILER
46 #include <stdexcept>
47 #include <limits>
48
49 // For now, work around MSVC++ 6.0's broken template support. If this value
50 // is not defined, the old GDALCopyWords implementation is used.
51 #define USE_NEW_COPYWORDS 1
52 #endif
53
54
55 CPL_CVSID("$Id: rasterio.cpp 29161 2015-05-06 10:18:19Z rouault $");
56
57 /************************************************************************/
58 /* IRasterIO() */
59 /* */
60 /* Default internal implementation of RasterIO() ... utilizes */
61 /* the Block access methods to satisfy the request. This would */
62 /* normally only be overridden by formats with overviews. */
63 /************************************************************************/
64
IRasterIO(GDALRWFlag eRWFlag,int nXOff,int nYOff,int nXSize,int nYSize,void * pData,int nBufXSize,int nBufYSize,GDALDataType eBufType,GSpacing nPixelSpace,GSpacing nLineSpace,GDALRasterIOExtraArg * psExtraArg)65 CPLErr GDALRasterBand::IRasterIO( GDALRWFlag eRWFlag,
66 int nXOff, int nYOff, int nXSize, int nYSize,
67 void * pData, int nBufXSize, int nBufYSize,
68 GDALDataType eBufType,
69 GSpacing nPixelSpace, GSpacing nLineSpace,
70 GDALRasterIOExtraArg* psExtraArg )
71
72 {
73 int nBandDataSize = GDALGetDataTypeSize( eDataType ) / 8;
74 int nBufDataSize = GDALGetDataTypeSize( eBufType ) / 8;
75 GByte *pabySrcBlock = NULL;
76 GDALRasterBlock *poBlock = NULL;
77 int nLBlockX=-1, nLBlockY=-1, iBufYOff, iBufXOff, iSrcY;
78
79 if( eRWFlag == GF_Write && eFlushBlockErr != CE_None )
80 {
81 CPLError(eFlushBlockErr, CPLE_AppDefined,
82 "An error occured while writing a dirty block");
83 CPLErr eErr = eFlushBlockErr;
84 eFlushBlockErr = CE_None;
85 return eErr;
86 }
87
88 /* ==================================================================== */
89 /* A common case is the data requested with the destination */
90 /* is packed, and the block width is the raster width. */
91 /* ==================================================================== */
92 if( nPixelSpace == nBufDataSize
93 && nLineSpace == nPixelSpace * nXSize
94 && nBlockXSize == GetXSize()
95 && nBufXSize == nXSize
96 && nBufYSize == nYSize )
97 {
98 // printf( "IRasterIO(%d,%d,%d,%d) rw=%d case 1\n",
99 // nXOff, nYOff, nXSize, nYSize,
100 // (int) eRWFlag );
101 CPLErr eErr = CE_None;
102
103 for( iBufYOff = 0; iBufYOff < nBufYSize; iBufYOff++ )
104 {
105 int nSrcByteOffset;
106
107 iSrcY = iBufYOff + nYOff;
108
109 if( iSrcY < nLBlockY * nBlockYSize
110 || iSrcY >= (nLBlockY+1) * nBlockYSize )
111 {
112 nLBlockY = iSrcY / nBlockYSize;
113 int bJustInitialize =
114 eRWFlag == GF_Write
115 && nXOff == 0 && nXSize == nBlockXSize
116 && nYOff <= nLBlockY * nBlockYSize
117 && nYOff + nYSize >= (nLBlockY+1) * nBlockYSize;
118
119 /* Is this a partial tile at right and/or bottom edges of */
120 /* the raster, and that is going to be completely written ? */
121 /* If so, don't load it from storage, but zeroized it so that */
122 /* the content outsize of the validity area is initialized */
123 int bMemZeroBuffer = FALSE;
124 if( eRWFlag == GF_Write && !bJustInitialize &&
125 nXOff == 0 && nXSize == nBlockXSize &&
126 nYOff <= nLBlockY * nBlockYSize &&
127 nYOff + nYSize == GetYSize() &&
128 (nLBlockY+1) * nBlockYSize > GetYSize() )
129 {
130 bJustInitialize = TRUE;
131 bMemZeroBuffer = TRUE;
132 }
133
134 if( poBlock )
135 poBlock->DropLock();
136
137 poBlock = GetLockedBlockRef( 0, nLBlockY, bJustInitialize );
138 if( poBlock == NULL )
139 {
140 CPLError( CE_Failure, CPLE_AppDefined,
141 "GetBlockRef failed at X block offset %d, "
142 "Y block offset %d", 0, nLBlockY );
143 eErr = CE_Failure;
144 break;
145 }
146
147 if( eRWFlag == GF_Write )
148 poBlock->MarkDirty();
149
150 pabySrcBlock = (GByte *) poBlock->GetDataRef();
151 if( pabySrcBlock == NULL )
152 {
153 poBlock->DropLock();
154 eErr = CE_Failure;
155 break;
156 }
157 if( bMemZeroBuffer )
158 {
159 memset(pabySrcBlock, 0,
160 nBandDataSize * nBlockXSize * nBlockYSize);
161 }
162 }
163
164 nSrcByteOffset = ((iSrcY-nLBlockY*nBlockYSize)*nBlockXSize + nXOff)
165 * nBandDataSize;
166
167 if( eDataType == eBufType )
168 {
169 if( eRWFlag == GF_Read )
170 memcpy( ((GByte *) pData) + (GPtrDiff_t)iBufYOff * nLineSpace,
171 pabySrcBlock + nSrcByteOffset,
172 (size_t)nLineSpace );
173 else
174 memcpy( pabySrcBlock + nSrcByteOffset,
175 ((GByte *) pData) + (GPtrDiff_t)iBufYOff * nLineSpace,
176 (size_t)nLineSpace );
177 }
178 else
179 {
180 /* type to type conversion */
181
182 if( eRWFlag == GF_Read )
183 GDALCopyWords( pabySrcBlock + nSrcByteOffset,
184 eDataType, nBandDataSize,
185 ((GByte *) pData) + (GPtrDiff_t)iBufYOff * nLineSpace,
186 eBufType, (int)nPixelSpace, nBufXSize );
187 else
188 GDALCopyWords( ((GByte *) pData) + (GPtrDiff_t)iBufYOff * nLineSpace,
189 eBufType, (int)nPixelSpace,
190 pabySrcBlock + nSrcByteOffset,
191 eDataType, nBandDataSize, nBufXSize );
192 }
193
194 if( psExtraArg->pfnProgress != NULL &&
195 !psExtraArg->pfnProgress(1.0 * (iBufYOff + 1) / nBufYSize, "",
196 psExtraArg->pProgressData) )
197 {
198 eErr = CE_Failure;
199 break;
200 }
201 }
202
203 if( poBlock )
204 poBlock->DropLock();
205
206 return eErr;
207 }
208
209 /* ==================================================================== */
210 /* Do we have overviews that would be appropriate to satisfy */
211 /* this request? */
212 /* ==================================================================== */
213 if( (nBufXSize < nXSize || nBufYSize < nYSize)
214 && GetOverviewCount() > 0 && eRWFlag == GF_Read )
215 {
216 int nOverview;
217 GDALRasterIOExtraArg sExtraArg;
218
219 GDALCopyRasterIOExtraArg(&sExtraArg, psExtraArg);
220
221 nOverview =
222 GDALBandGetBestOverviewLevel2(this, nXOff, nYOff, nXSize, nYSize,
223 nBufXSize, nBufYSize, &sExtraArg);
224 if (nOverview >= 0)
225 {
226 GDALRasterBand* poOverviewBand = GetOverview(nOverview);
227 if (poOverviewBand == NULL)
228 return CE_Failure;
229
230 return poOverviewBand->RasterIO( eRWFlag, nXOff, nYOff, nXSize, nYSize,
231 pData, nBufXSize, nBufYSize, eBufType,
232 nPixelSpace, nLineSpace, &sExtraArg );
233 }
234 }
235
236
237 if( eRWFlag == GF_Read &&
238 nBufXSize < nXSize / 100 && nBufYSize < nYSize / 100 &&
239 nPixelSpace == nBufDataSize &&
240 nLineSpace == nPixelSpace * nBufXSize &&
241 CSLTestBoolean(CPLGetConfigOption("GDAL_NO_COSTLY_OVERVIEW", "NO")) )
242 {
243 memset(pData, 0, (size_t)(nLineSpace * nBufYSize));
244 return CE_None;
245 }
246
247
248 /* ==================================================================== */
249 /* The second case when we don't need subsample data but likely */
250 /* need data type conversion. */
251 /* ==================================================================== */
252 int iSrcX;
253
254 if ( /* nPixelSpace == nBufDataSize
255 && */ nXSize == nBufXSize
256 && nYSize == nBufYSize )
257 {
258 // printf( "IRasterIO(%d,%d,%d,%d) rw=%d case 2\n",
259 // nXOff, nYOff, nXSize, nYSize,
260 // (int) eRWFlag );
261
262 /* -------------------------------------------------------------------- */
263 /* Loop over buffer computing source locations. */
264 /* -------------------------------------------------------------------- */
265 int nLBlockXStart, nXSpanEnd;
266
267 // Calculate starting values out of loop
268 nLBlockXStart = nXOff / nBlockXSize;
269 nXSpanEnd = nBufXSize + nXOff;
270
271 int nYInc = 0;
272 for( iBufYOff = 0, iSrcY = nYOff; iBufYOff < nBufYSize; iBufYOff+=nYInc, iSrcY+=nYInc )
273 {
274 GPtrDiff_t iBufOffset;
275 GPtrDiff_t iSrcOffset;
276 int nXSpan;
277
278 iBufOffset = (GPtrDiff_t)iBufYOff * (GPtrDiff_t)nLineSpace;
279 nLBlockY = iSrcY / nBlockYSize;
280 nLBlockX = nLBlockXStart;
281 iSrcX = nXOff;
282 while( iSrcX < nXSpanEnd )
283 {
284 int nXSpanSize;
285
286 nXSpan = (nLBlockX + 1) * nBlockXSize;
287 nXSpan = ( ( nXSpan < nXSpanEnd )?nXSpan:nXSpanEnd ) - iSrcX;
288 nXSpanSize = nXSpan * (int)nPixelSpace;
289
290 int bJustInitialize =
291 eRWFlag == GF_Write
292 && nYOff <= nLBlockY * nBlockYSize
293 && nYOff + nYSize >= (nLBlockY+1) * nBlockYSize
294 && nXOff <= nLBlockX * nBlockXSize
295 && nXOff + nXSize >= (nLBlockX+1) * nBlockXSize;
296
297 /* Is this a partial tile at right and/or bottom edges of */
298 /* the raster, and that is going to be completely written ? */
299 /* If so, don't load it from storage, but zeroized it so that */
300 /* the content outsize of the validity area is initialized */
301 int bMemZeroBuffer = FALSE;
302 if( eRWFlag == GF_Write && !bJustInitialize &&
303 nXOff <= nLBlockX * nBlockXSize &&
304 nYOff <= nLBlockY * nBlockYSize &&
305 (nXOff + nXSize >= (nLBlockX+1) * nBlockXSize ||
306 (nXOff + nXSize == GetXSize() && (nLBlockX+1) * nBlockXSize > GetXSize())) &&
307 (nYOff + nYSize >= (nLBlockY+1) * nBlockYSize ||
308 (nYOff + nYSize == GetYSize() && (nLBlockY+1) * nBlockYSize > GetYSize())) )
309 {
310 bJustInitialize = TRUE;
311 bMemZeroBuffer = TRUE;
312 }
313 // printf( "bJustInitialize = %d (%d,%d,%d,%d)\n",
314 // bJustInitialize,
315 // nYOff, nYSize,
316 // nLBlockY, nBlockYSize );
317 // bJustInitialize = FALSE;
318
319
320 /* -------------------------------------------------------------------- */
321 /* Ensure we have the appropriate block loaded. */
322 /* -------------------------------------------------------------------- */
323 poBlock = GetLockedBlockRef( nLBlockX, nLBlockY, bJustInitialize );
324 if( !poBlock )
325 {
326 CPLError( CE_Failure, CPLE_AppDefined,
327 "GetBlockRef failed at X block offset %d, "
328 "Y block offset %d", nLBlockX, nLBlockY );
329 return( CE_Failure );
330 }
331
332 if( eRWFlag == GF_Write )
333 poBlock->MarkDirty();
334
335 pabySrcBlock = (GByte *) poBlock->GetDataRef();
336 if( pabySrcBlock == NULL )
337 {
338 poBlock->DropLock();
339 return CE_Failure;
340 }
341 if( bMemZeroBuffer )
342 {
343 memset(pabySrcBlock, 0,
344 nBandDataSize * nBlockXSize * nBlockYSize);
345 }
346 /* -------------------------------------------------------------------- */
347 /* Copy over this chunk of data. */
348 /* -------------------------------------------------------------------- */
349 iSrcOffset = ((GPtrDiff_t)iSrcX - (GPtrDiff_t)nLBlockX*nBlockXSize
350 + ((GPtrDiff_t)(iSrcY) - (GPtrDiff_t)nLBlockY*nBlockYSize) * nBlockXSize)*nBandDataSize;
351 /* Fill up as many rows as possible for the loaded block */
352 int kmax = MIN(nBlockYSize - (iSrcY % nBlockYSize), nBufYSize - iBufYOff);
353 for(int k=0; k<kmax;k++)
354 {
355 if( eDataType == eBufType
356 && nPixelSpace == nBufDataSize )
357 {
358 if( eRWFlag == GF_Read )
359 memcpy( ((GByte *) pData) + iBufOffset + (GPtrDiff_t)k * nLineSpace,
360 pabySrcBlock + iSrcOffset, nXSpanSize );
361 else
362 memcpy( pabySrcBlock + iSrcOffset,
363 ((GByte *) pData) + iBufOffset + (GPtrDiff_t)k * nLineSpace, nXSpanSize );
364 }
365 else
366 {
367 /* type to type conversion */
368
369 if( eRWFlag == GF_Read )
370 GDALCopyWords( pabySrcBlock + iSrcOffset,
371 eDataType, nBandDataSize,
372 ((GByte *) pData) + iBufOffset + (GPtrDiff_t)k * nLineSpace,
373 eBufType, (int)nPixelSpace, nXSpan );
374 else
375 GDALCopyWords( ((GByte *) pData) + iBufOffset + (GPtrDiff_t)k * nLineSpace,
376 eBufType, (int)nPixelSpace,
377 pabySrcBlock + iSrcOffset,
378 eDataType, nBandDataSize, nXSpan );
379 }
380
381 iSrcOffset += nBlockXSize * nBandDataSize;
382 }
383
384 iBufOffset += nXSpanSize;
385 nLBlockX++;
386 iSrcX+=nXSpan;
387
388 poBlock->DropLock();
389 poBlock = NULL;
390 }
391
392 /* Compute the increment to go on a block boundary */
393 nYInc = nBlockYSize - (iSrcY % nBlockYSize);
394
395 if( psExtraArg->pfnProgress != NULL &&
396 !psExtraArg->pfnProgress(1.0 * MIN(nBufYSize, iBufYOff + nYInc) / nBufYSize, "",
397 psExtraArg->pProgressData) )
398 {
399 return CE_Failure;
400 }
401 }
402
403 return CE_None;
404 }
405
406 /* ==================================================================== */
407 /* Loop reading required source blocks to satisfy output */
408 /* request. This is the most general implementation. */
409 /* ==================================================================== */
410
411 double dfXOff, dfYOff, dfXSize, dfYSize;
412 if( psExtraArg->bFloatingPointWindowValidity )
413 {
414 dfXOff = psExtraArg->dfXOff;
415 dfYOff = psExtraArg->dfYOff;
416 dfXSize = psExtraArg->dfXSize;
417 dfYSize = psExtraArg->dfYSize;
418 CPLAssert(dfXOff - nXOff < 1.0);
419 CPLAssert(dfYOff - nYOff < 1.0);
420 CPLAssert(nXSize - dfXSize < 1.0);
421 CPLAssert(nYSize - dfYSize < 1.0);
422 }
423 else
424 {
425 dfXOff = nXOff;
426 dfYOff = nYOff;
427 dfXSize = nXSize;
428 dfYSize = nYSize;
429 }
430 /* -------------------------------------------------------------------- */
431 /* Compute stepping increment. */
432 /* -------------------------------------------------------------------- */
433 double dfSrcXInc, dfSrcYInc;
434 dfSrcXInc = dfXSize / (double) nBufXSize;
435 dfSrcYInc = dfYSize / (double) nBufYSize;
436 CPLErr eErr = CE_None;
437
438 // printf( "IRasterIO(%d,%d,%d,%d) rw=%d case 3\n",
439 // nXOff, nYOff, nXSize, nYSize,
440 // (int) eRWFlag );
441 if (eRWFlag == GF_Write)
442 {
443 /* -------------------------------------------------------------------- */
444 /* Write case */
445 /* Loop over raster window computing source locations in the buffer. */
446 /* -------------------------------------------------------------------- */
447 int iDstX, iDstY;
448 GByte* pabyDstBlock = NULL;
449
450 for( iDstY = nYOff; iDstY < nYOff + nYSize; iDstY ++)
451 {
452 GPtrDiff_t iBufOffset, iDstOffset;
453 iBufYOff = (int)((iDstY - nYOff) / dfSrcYInc);
454
455 for( iDstX = nXOff; iDstX < nXOff + nXSize; iDstX ++)
456 {
457 iBufXOff = (int)((iDstX - nXOff) / dfSrcXInc);
458 iBufOffset = (GPtrDiff_t)iBufYOff * (GPtrDiff_t)nLineSpace + iBufXOff * (GPtrDiff_t)nPixelSpace;
459
460 // FIXME: this code likely doesn't work if the dirty block gets flushed
461 // to disk before being completely written.
462 // In the meantime, bJustInitalize should probably be set to FALSE
463 // even if it is not ideal performance wise, and for lossy compression
464
465 /* -------------------------------------------------------------------- */
466 /* Ensure we have the appropriate block loaded. */
467 /* -------------------------------------------------------------------- */
468 if( iDstX < nLBlockX * nBlockXSize
469 || iDstX >= (nLBlockX+1) * nBlockXSize
470 || iDstY < nLBlockY * nBlockYSize
471 || iDstY >= (nLBlockY+1) * nBlockYSize )
472 {
473 nLBlockX = iDstX / nBlockXSize;
474 nLBlockY = iDstY / nBlockYSize;
475
476 int bJustInitialize =
477 nYOff <= nLBlockY * nBlockYSize
478 && nYOff + nYSize >= (nLBlockY+1) * nBlockYSize
479 && nXOff <= nLBlockX * nBlockXSize
480 && nXOff + nXSize >= (nLBlockX+1) * nBlockXSize;
481 /*int bMemZeroBuffer = FALSE;
482 if( !bJustInitialize &&
483 nXOff <= nLBlockX * nBlockXSize &&
484 nYOff <= nLBlockY * nBlockYSize &&
485 (nXOff + nXSize >= (nLBlockX+1) * nBlockXSize ||
486 (nXOff + nXSize == GetXSize() && (nLBlockX+1) * nBlockXSize > GetXSize())) &&
487 (nYOff + nYSize >= (nLBlockY+1) * nBlockYSize ||
488 (nYOff + nYSize == GetYSize() && (nLBlockY+1) * nBlockYSize > GetYSize())) )
489 {
490 bJustInitialize = TRUE;
491 bMemZeroBuffer = TRUE;
492 }*/
493 if( poBlock != NULL )
494 poBlock->DropLock();
495
496 poBlock = GetLockedBlockRef( nLBlockX, nLBlockY,
497 bJustInitialize );
498 if( poBlock == NULL )
499 {
500 return( CE_Failure );
501 }
502
503 poBlock->MarkDirty();
504
505 pabyDstBlock = (GByte *) poBlock->GetDataRef();
506 if( pabyDstBlock == NULL )
507 {
508 poBlock->DropLock();
509 return CE_Failure;
510 }
511 /*if( bMemZeroBuffer )
512 {
513 memset(pabyDstBlock, 0,
514 nBandDataSize * nBlockXSize * nBlockYSize);
515 }*/
516 }
517
518 /* -------------------------------------------------------------------- */
519 /* Copy over this pixel of data. */
520 /* -------------------------------------------------------------------- */
521 iDstOffset = ((GPtrDiff_t)iDstX - (GPtrDiff_t)nLBlockX*nBlockXSize
522 + ((GPtrDiff_t)iDstY - (GPtrDiff_t)nLBlockY*nBlockYSize) * nBlockXSize)*nBandDataSize;
523
524 if( eDataType == eBufType )
525 {
526 memcpy( pabyDstBlock + iDstOffset,
527 ((GByte *) pData) + iBufOffset, nBandDataSize );
528 }
529 else
530 {
531 /* type to type conversion ... ouch, this is expensive way
532 of handling single words */
533
534 GDALCopyWords( ((GByte *) pData) + iBufOffset, eBufType, 0,
535 pabyDstBlock + iDstOffset, eDataType, 0,
536 1 );
537 }
538 }
539
540 if( psExtraArg->pfnProgress != NULL &&
541 !psExtraArg->pfnProgress(1.0 * (iDstY - nYOff + 1) / nYSize, "",
542 psExtraArg->pProgressData) )
543 {
544 eErr = CE_Failure;
545 break;
546 }
547 }
548 }
549 else
550 {
551 if( psExtraArg->eResampleAlg != GRIORA_NearestNeighbour )
552 {
553 if( (psExtraArg->eResampleAlg == GRIORA_Cubic ||
554 psExtraArg->eResampleAlg == GRIORA_CubicSpline ||
555 psExtraArg->eResampleAlg == GRIORA_Bilinear ||
556 psExtraArg->eResampleAlg == GRIORA_Lanczos) &&
557 GetColorTable() != NULL )
558 {
559 CPLError(CE_Warning, CPLE_NotSupported,
560 "Resampling method not supported on paletted band. "
561 "Falling back to nearest neighbour");
562 }
563 else if( psExtraArg->eResampleAlg == GRIORA_Gauss &&
564 GDALDataTypeIsComplex( eDataType ) )
565 {
566 CPLError(CE_Warning, CPLE_NotSupported,
567 "Resampling method not supported on complex data type band. "
568 "Falling back to nearest neighbour");
569 }
570 else
571 {
572 return RasterIOResampled( eRWFlag,
573 nXOff, nYOff, nXSize, nYSize,
574 pData, nBufXSize, nBufYSize,
575 eBufType,
576 nPixelSpace, nLineSpace,
577 psExtraArg );
578 }
579 }
580
581 double dfSrcX, dfSrcY;
582 int nLimitBlockY = 0;
583 int bByteCopy = ( eDataType == eBufType && nBandDataSize == 1);
584 int nStartBlockX = -nBlockXSize;
585
586 /* -------------------------------------------------------------------- */
587 /* Read case */
588 /* Loop over buffer computing source locations. */
589 /* -------------------------------------------------------------------- */
590 for( iBufYOff = 0; iBufYOff < nBufYSize; iBufYOff++ )
591 {
592 GPtrDiff_t iBufOffset, iSrcOffset;
593
594 dfSrcY = (iBufYOff+0.5) * dfSrcYInc + dfYOff;
595 dfSrcX = 0.5 * dfSrcXInc + dfXOff;
596 iSrcY = (int) dfSrcY;
597
598 iBufOffset = (GPtrDiff_t)iBufYOff * (GPtrDiff_t)nLineSpace;
599
600 if( iSrcY >= nLimitBlockY )
601 {
602 nLBlockY = iSrcY / nBlockYSize;
603 nLimitBlockY = (nLBlockY + 1) * nBlockYSize;
604 nStartBlockX = -nBlockXSize; /* make sure a new block is loaded */
605 }
606 else if( (int)dfSrcX < nStartBlockX )
607 nStartBlockX = -nBlockXSize; /* make sure a new block is loaded */
608
609 GPtrDiff_t iSrcOffsetCst = (iSrcY - nLBlockY*nBlockYSize) * (GPtrDiff_t)nBlockXSize;
610
611 for( iBufXOff = 0; iBufXOff < nBufXSize; iBufXOff++, dfSrcX += dfSrcXInc )
612 {
613 iSrcX = (int) dfSrcX;
614 int nDiffX = iSrcX - nStartBlockX;
615
616 /* -------------------------------------------------------------------- */
617 /* Ensure we have the appropriate block loaded. */
618 /* -------------------------------------------------------------------- */
619 if( nDiffX >= nBlockXSize )
620 {
621 nLBlockX = iSrcX / nBlockXSize;
622 nStartBlockX = nLBlockX * nBlockXSize;
623 nDiffX = iSrcX - nStartBlockX;
624
625 if( poBlock != NULL )
626 poBlock->DropLock();
627
628 poBlock = GetLockedBlockRef( nLBlockX, nLBlockY,
629 FALSE );
630 if( poBlock == NULL )
631 {
632 eErr = CE_Failure;
633 break;
634 }
635
636 pabySrcBlock = (GByte *) poBlock->GetDataRef();
637 if( pabySrcBlock == NULL )
638 {
639 eErr = CE_Failure;
640 break;
641 }
642 }
643
644 /* -------------------------------------------------------------------- */
645 /* Copy over this pixel of data. */
646 /* -------------------------------------------------------------------- */
647 iSrcOffset = ((GPtrDiff_t)nDiffX + iSrcOffsetCst)*nBandDataSize;
648
649 if( bByteCopy )
650 {
651 ((GByte *) pData)[iBufOffset] = pabySrcBlock[iSrcOffset];
652 }
653 else if( eDataType == eBufType )
654 {
655 memcpy( ((GByte *) pData) + iBufOffset,
656 pabySrcBlock + iSrcOffset, nBandDataSize );
657 }
658 else
659 {
660 /* type to type conversion ... ouch, this is expensive way
661 of handling single words */
662 GDALCopyWords( pabySrcBlock + iSrcOffset, eDataType, 0,
663 ((GByte *) pData) + iBufOffset, eBufType, 0,
664 1 );
665 }
666
667 iBufOffset += (int)nPixelSpace;
668 }
669 if( eErr == CE_Failure )
670 break;
671
672 if( psExtraArg->pfnProgress != NULL &&
673 !psExtraArg->pfnProgress(1.0 * (iBufYOff + 1) / nBufYSize, "",
674 psExtraArg->pProgressData) )
675 {
676 eErr = CE_Failure;
677 break;
678 }
679 }
680 }
681
682 if( poBlock != NULL )
683 poBlock->DropLock();
684
685 return( eErr );
686 }
687
688 /************************************************************************/
689 /* GDALRasterIOTransformer() */
690 /************************************************************************/
691
692 typedef struct
693 {
694 double dfXOff, dfYOff;
695 double dfXRatioDstToSrc, dfYRatioDstToSrc;
696 } GDALRasterIOTransformerStruct;
697
GDALRasterIOTransformer(void * pTransformerArg,CPL_UNUSED int bDstToSrc,int nPointCount,double * x,double * y,CPL_UNUSED double * z,int * panSuccess)698 static int GDALRasterIOTransformer( void *pTransformerArg,
699 CPL_UNUSED int bDstToSrc, int nPointCount,
700 double *x, double *y, CPL_UNUSED double *z,
701 int *panSuccess )
702 {
703 CPLAssert(bDstToSrc);
704 GDALRasterIOTransformerStruct* psParams = (GDALRasterIOTransformerStruct*) pTransformerArg;
705 for(int i = 0; i < nPointCount; i++)
706 {
707 x[i] = x[i] * psParams->dfXRatioDstToSrc + psParams->dfXOff;
708 y[i] = y[i] * psParams->dfYRatioDstToSrc + psParams->dfYOff;
709 panSuccess[i] = TRUE;
710 }
711 return TRUE;
712 }
713
714 /************************************************************************/
715 /* RasterIOResampled() */
716 /************************************************************************/
717
RasterIOResampled(CPL_UNUSED GDALRWFlag eRWFlag,int nXOff,int nYOff,int nXSize,int nYSize,void * pData,int nBufXSize,int nBufYSize,GDALDataType eBufType,GSpacing nPixelSpace,GSpacing nLineSpace,GDALRasterIOExtraArg * psExtraArg)718 CPLErr GDALRasterBand::RasterIOResampled( CPL_UNUSED GDALRWFlag eRWFlag,
719 int nXOff, int nYOff, int nXSize, int nYSize,
720 void * pData, int nBufXSize, int nBufYSize,
721 GDALDataType eBufType,
722 GSpacing nPixelSpace, GSpacing nLineSpace,
723 GDALRasterIOExtraArg* psExtraArg )
724 {
725 CPLErr eErr = CE_None;
726
727 // Determine if we use warping resampling or overview resampling
728 int bUseWarp;
729 if( !GDALDataTypeIsComplex( eDataType ) )
730 bUseWarp = FALSE;
731 else
732 bUseWarp = TRUE;
733
734 double dfXOff, dfYOff, dfXSize, dfYSize;
735 if( psExtraArg->bFloatingPointWindowValidity )
736 {
737 dfXOff = psExtraArg->dfXOff;
738 dfYOff = psExtraArg->dfYOff;
739 dfXSize = psExtraArg->dfXSize;
740 dfYSize = psExtraArg->dfYSize;
741 CPLAssert(dfXOff - nXOff < 1.0);
742 CPLAssert(dfYOff - nYOff < 1.0);
743 CPLAssert(nXSize - dfXSize < 1.0);
744 CPLAssert(nYSize - dfYSize < 1.0);
745 }
746 else
747 {
748 dfXOff = nXOff;
749 dfYOff = nYOff;
750 dfXSize = nXSize;
751 dfYSize = nYSize;
752 }
753
754 double dfXRatioDstToSrc = dfXSize / nBufXSize;
755 double dfYRatioDstToSrc = dfYSize / nBufYSize;
756
757 /* Determine the coordinates in the "virtual" output raster to see */
758 /* if there are not integers, in which case we will use them as a shift */
759 /* so that subwindow extracts give the exact same results as entire raster */
760 /* scaling */
761 double dfDestXOff = dfXOff / dfXRatioDstToSrc;
762 int bHasXOffVirtual = FALSE;
763 int nDestXOffVirtual = 0;
764 if( fabs(dfDestXOff - (int)(dfDestXOff + 0.5)) < 1e-8 )
765 {
766 bHasXOffVirtual = TRUE;
767 dfXOff = nXOff;
768 nDestXOffVirtual = (int)(dfDestXOff + 0.5);
769 }
770
771 double dfDestYOff = dfYOff / dfYRatioDstToSrc;
772 int bHasYOffVirtual = FALSE;
773 int nDestYOffVirtual = 0;
774 if( fabs(dfDestYOff - (int)(dfDestYOff + 0.5)) < 1e-8 )
775 {
776 bHasYOffVirtual = TRUE;
777 dfYOff = nYOff;
778 nDestYOffVirtual = (int)(dfDestYOff + 0.5);
779 }
780
781 /* Create a MEM dataset that wraps the output buffer */
782 GDALDataset* poMEMDS = MEMDataset::Create("", nDestXOffVirtual + nBufXSize,
783 nDestYOffVirtual + nBufYSize, 0,
784 eBufType, NULL);
785 char szBuffer[64];
786 int nRet;
787
788 nRet = CPLPrintPointer(szBuffer, (GByte*)pData - nPixelSpace * nDestXOffVirtual
789 - nLineSpace * nDestYOffVirtual, sizeof(szBuffer));
790 szBuffer[nRet] = 0;
791 char** papszOptions = CSLSetNameValue(NULL, "DATAPOINTER", szBuffer);
792
793 papszOptions = CSLSetNameValue(papszOptions, "PIXELOFFSET",
794 CPLSPrintf(CPL_FRMT_GIB, (GIntBig)nPixelSpace));
795
796 papszOptions = CSLSetNameValue(papszOptions, "LINEOFFSET",
797 CPLSPrintf(CPL_FRMT_GIB, (GIntBig)nLineSpace));
798
799 poMEMDS->AddBand(eBufType, papszOptions);
800 CSLDestroy(papszOptions);
801 GDALRasterBandH hMEMBand = (GDALRasterBandH)poMEMDS->GetRasterBand(1);
802
803 /* Do the resampling */
804 if( bUseWarp )
805 {
806 VRTDatasetH hVRTDS = NULL;
807 GDALRasterBandH hVRTBand = NULL;
808 if( GetDataset() == NULL )
809 {
810 /* Create VRT dataset that wraps the whole dataset */
811 hVRTDS = VRTCreate(nRasterXSize, nRasterYSize);
812 VRTAddBand( hVRTDS, eDataType, NULL );
813 hVRTBand = GDALGetRasterBand(hVRTDS, 1);
814 VRTAddSimpleSource( (VRTSourcedRasterBandH)hVRTBand,
815 (GDALRasterBandH)this,
816 0, 0,
817 nRasterXSize, nRasterYSize,
818 0, 0,
819 nRasterXSize, nRasterYSize,
820 NULL, VRT_NODATA_UNSET );
821
822 /* Add a mask band if needed */
823 if( GetMaskFlags() != GMF_ALL_VALID )
824 {
825 ((GDALDataset*)hVRTDS)->CreateMaskBand(0);
826 VRTSourcedRasterBand* poVRTMaskBand =
827 (VRTSourcedRasterBand*)(((GDALRasterBand*)hVRTBand)->GetMaskBand());
828 poVRTMaskBand->
829 AddMaskBandSource( this,
830 0, 0,
831 nRasterXSize, nRasterYSize,
832 0, 0,
833 nRasterXSize, nRasterYSize);
834 }
835 }
836
837 GDALWarpOptions* psWarpOptions = GDALCreateWarpOptions();
838 psWarpOptions->eResampleAlg = (GDALResampleAlg)psExtraArg->eResampleAlg;
839 psWarpOptions->hSrcDS = (GDALDatasetH) (hVRTDS ? hVRTDS : GetDataset());
840 psWarpOptions->hDstDS = (GDALDatasetH) poMEMDS;
841 psWarpOptions->nBandCount = 1;
842 int nSrcBandNumber = (hVRTDS ? 1 : nBand);
843 int nDstBandNumber = 1;
844 psWarpOptions->panSrcBands = &nSrcBandNumber;
845 psWarpOptions->panDstBands = &nDstBandNumber;
846 psWarpOptions->pfnProgress = psExtraArg->pfnProgress ?
847 psExtraArg->pfnProgress : GDALDummyProgress;
848 psWarpOptions->pProgressArg = psExtraArg->pProgressData;
849 psWarpOptions->pfnTransformer = GDALRasterIOTransformer;
850 GDALRasterIOTransformerStruct sTransformer;
851 sTransformer.dfXOff = (bHasXOffVirtual) ? 0 : dfXOff;
852 sTransformer.dfYOff = (bHasYOffVirtual) ? 0 : dfYOff;
853 sTransformer.dfXRatioDstToSrc = dfXRatioDstToSrc;
854 sTransformer.dfYRatioDstToSrc = dfYRatioDstToSrc;
855 psWarpOptions->pTransformerArg = &sTransformer;
856
857 GDALWarpOperationH hWarpOperation = GDALCreateWarpOperation(psWarpOptions);
858 eErr = GDALChunkAndWarpImage( hWarpOperation,
859 nDestXOffVirtual, nDestYOffVirtual,
860 nBufXSize, nBufYSize );
861 GDALDestroyWarpOperation( hWarpOperation );
862
863 psWarpOptions->panSrcBands = NULL;
864 psWarpOptions->panDstBands = NULL;
865 GDALDestroyWarpOptions( psWarpOptions );
866
867 if( hVRTDS )
868 GDALClose(hVRTDS);
869 }
870 else
871 {
872 const char* pszResampling =
873 (psExtraArg->eResampleAlg == GRIORA_Bilinear) ? "BILINEAR" :
874 (psExtraArg->eResampleAlg == GRIORA_Cubic) ? "CUBIC" :
875 (psExtraArg->eResampleAlg == GRIORA_CubicSpline) ? "CUBICSPLINE" :
876 (psExtraArg->eResampleAlg == GRIORA_Lanczos) ? "LANCZOS" :
877 (psExtraArg->eResampleAlg == GRIORA_Average) ? "AVERAGE" :
878 (psExtraArg->eResampleAlg == GRIORA_Mode) ? "MODE" :
879 (psExtraArg->eResampleAlg == GRIORA_Gauss) ? "GAUSS" : "UNKNOWN";
880
881 int nKernelRadius;
882 GDALResampleFunction pfnResampleFunc =
883 GDALGetResampleFunction(pszResampling, &nKernelRadius);
884 CPLAssert(pfnResampleFunc);
885 GDALDataType eWrkDataType =
886 GDALGetOvrWorkDataType(pszResampling, eDataType);
887 int bHasNoData = FALSE;
888 float fNoDataValue = (float) GetNoDataValue(&bHasNoData);
889 if( !bHasNoData )
890 fNoDataValue = 0.0f;
891
892 int nDstBlockXSize = nBufXSize;
893 int nDstBlockYSize = nBufYSize;
894 int nFullResXChunk, nFullResYChunk;
895 while(TRUE)
896 {
897 nFullResXChunk = 3 + (int)(nDstBlockXSize * dfXRatioDstToSrc);
898 nFullResYChunk = 3 + (int)(nDstBlockYSize * dfYRatioDstToSrc);
899 if( (nDstBlockXSize == 1 && nDstBlockYSize == 1) ||
900 ((GIntBig)nFullResXChunk * nFullResYChunk <= 1024 * 1024) )
901 break;
902 /* When operating on the full width of a raster whose block width is */
903 /* the raster width, prefer doing chunks in height */
904 if( nFullResXChunk >= nXSize && nXSize == nBlockXSize && nDstBlockYSize > 1 )
905 nDstBlockYSize /= 2;
906 /* Otherwise cut the maximal dimension */
907 else if( nDstBlockXSize > 1 && nFullResXChunk > nFullResYChunk )
908 nDstBlockXSize /= 2;
909 else
910 nDstBlockYSize /= 2;
911 }
912
913 int nOvrFactor = MAX( (int)(0.5 + dfXRatioDstToSrc),
914 (int)(0.5 + dfYRatioDstToSrc) );
915 if( nOvrFactor == 0 ) nOvrFactor = 1;
916 int nFullResXSizeQueried = nFullResXChunk + 2 * nKernelRadius * nOvrFactor;
917 int nFullResYSizeQueried = nFullResYChunk + 2 * nKernelRadius * nOvrFactor;
918
919 void * pChunk =
920 VSIMalloc3((GDALGetDataTypeSize(eWrkDataType)/8),
921 nFullResXSizeQueried, nFullResYSizeQueried );
922 GByte * pabyChunkNoDataMask = NULL;
923
924 GDALRasterBand* poMaskBand = NULL;
925 int nMaskFlags = 0;
926 int bUseNoDataMask = FALSE;
927
928 poMaskBand = GetMaskBand();
929 nMaskFlags = GetMaskFlags();
930
931 bUseNoDataMask = ((nMaskFlags & GMF_ALL_VALID) == 0);
932 if (bUseNoDataMask)
933 {
934 pabyChunkNoDataMask = (GByte *)
935 (GByte*) VSIMalloc2( nFullResXSizeQueried, nFullResYSizeQueried );
936 }
937 if( pChunk == NULL || (bUseNoDataMask && pabyChunkNoDataMask == NULL) )
938 {
939 GDALClose(poMEMDS);
940 CPLFree(pChunk);
941 CPLFree(pabyChunkNoDataMask);
942 CPLError( CE_Failure, CPLE_OutOfMemory,
943 "Out of memory in RasterIO()." );
944 return CE_Failure;
945 }
946
947 int nTotalBlocks = ((nBufXSize + nDstBlockXSize - 1) / nDstBlockXSize) *
948 ((nBufYSize + nDstBlockYSize - 1) / nDstBlockYSize);
949 int nBlocksDone = 0;
950
951 int nDstYOff;
952 for( nDstYOff = 0; nDstYOff < nBufYSize && eErr == CE_None;
953 nDstYOff += nDstBlockYSize )
954 {
955 int nDstYCount;
956 if (nDstYOff + nDstBlockYSize <= nBufYSize)
957 nDstYCount = nDstBlockYSize;
958 else
959 nDstYCount = nBufYSize - nDstYOff;
960
961 int nChunkYOff = nYOff + (int) (nDstYOff * dfYRatioDstToSrc);
962 int nChunkYOff2 = nYOff + 1 + (int) ceil((nDstYOff + nDstYCount) * dfYRatioDstToSrc);
963 if( nChunkYOff2 > nRasterYSize )
964 nChunkYOff2 = nRasterYSize;
965 int nYCount = nChunkYOff2 - nChunkYOff;
966 CPLAssert(nYCount <= nFullResYChunk);
967
968 int nChunkYOffQueried = nChunkYOff - nKernelRadius * nOvrFactor;
969 int nChunkYSizeQueried = nYCount + 2 * nKernelRadius * nOvrFactor;
970 if( nChunkYOffQueried < 0 )
971 {
972 nChunkYSizeQueried += nChunkYOffQueried;
973 nChunkYOffQueried = 0;
974 }
975 if( nChunkYSizeQueried + nChunkYOffQueried > nRasterYSize )
976 nChunkYSizeQueried = nRasterYSize - nChunkYOffQueried;
977 CPLAssert(nChunkYSizeQueried <= nFullResYSizeQueried);
978
979 int nDstXOff;
980 for( nDstXOff = 0; nDstXOff < nBufXSize && eErr == CE_None;
981 nDstXOff += nDstBlockXSize )
982 {
983 int nDstXCount;
984 if (nDstXOff + nDstBlockXSize <= nBufXSize)
985 nDstXCount = nDstBlockXSize;
986 else
987 nDstXCount = nBufXSize - nDstXOff;
988
989 int nChunkXOff = nXOff + (int) (nDstXOff * dfXRatioDstToSrc);
990 int nChunkXOff2 = nXOff + 1 + (int) ceil((nDstXOff + nDstXCount) * dfXRatioDstToSrc);
991 if( nChunkXOff2 > nRasterXSize )
992 nChunkXOff2 = nRasterXSize;
993 int nXCount = nChunkXOff2 - nChunkXOff;
994 CPLAssert(nXCount <= nFullResXChunk);
995
996 int nChunkXOffQueried = nChunkXOff - nKernelRadius * nOvrFactor;
997 int nChunkXSizeQueried = nXCount + 2 * nKernelRadius * nOvrFactor;
998 if( nChunkXOffQueried < 0 )
999 {
1000 nChunkXSizeQueried += nChunkXOffQueried;
1001 nChunkXOffQueried = 0;
1002 }
1003 if( nChunkXSizeQueried + nChunkXOffQueried > nRasterXSize )
1004 nChunkXSizeQueried = nRasterXSize - nChunkXOffQueried;
1005 CPLAssert(nChunkXSizeQueried <= nFullResXSizeQueried);
1006
1007 /* Read the source buffers */
1008 eErr = RasterIO( GF_Read,
1009 nChunkXOffQueried, nChunkYOffQueried,
1010 nChunkXSizeQueried, nChunkYSizeQueried,
1011 pChunk,
1012 nChunkXSizeQueried, nChunkYSizeQueried,
1013 eWrkDataType, 0, 0, NULL );
1014
1015 int bSkipResample = FALSE;
1016 int bNoDataMaskFullyOpaque = FALSE;
1017 if (eErr == CE_None && bUseNoDataMask)
1018 {
1019 eErr = poMaskBand->RasterIO( GF_Read,
1020 nChunkXOffQueried,
1021 nChunkYOffQueried,
1022 nChunkXSizeQueried,
1023 nChunkYSizeQueried,
1024 pabyChunkNoDataMask,
1025 nChunkXSizeQueried,
1026 nChunkYSizeQueried,
1027 GDT_Byte, 0, 0, NULL );
1028
1029 /* Optimizations if mask if fully opaque or transparent */
1030 int nPixels = nChunkXSizeQueried * nChunkYSizeQueried;
1031 GByte bVal = pabyChunkNoDataMask[0];
1032 int i = 1;
1033 for( ; i < nPixels; i++ )
1034 {
1035 if( pabyChunkNoDataMask[i] != bVal )
1036 break;
1037 }
1038 if( i == nPixels )
1039 {
1040 if( bVal == 0 )
1041 {
1042 for(int j=0;j<nDstYCount;j++)
1043 {
1044 GDALCopyWords(&fNoDataValue, GDT_Float32, 0,
1045 (GByte*)pData + nLineSpace * (j + nDstYOff) +
1046 nDstXOff * nPixelSpace,
1047 eBufType, (int)nPixelSpace,
1048 nDstXCount);
1049 }
1050 bSkipResample = TRUE;
1051 }
1052 else
1053 {
1054 bNoDataMaskFullyOpaque = TRUE;
1055 }
1056 }
1057 }
1058
1059 if( !bSkipResample && eErr == CE_None )
1060 {
1061 eErr = pfnResampleFunc( dfXRatioDstToSrc,
1062 dfYRatioDstToSrc,
1063 dfXOff - nXOff, /* == 0 if bHasXOffVirtual */
1064 dfYOff - nYOff, /* == 0 if bHasYOffVirtual */
1065 eWrkDataType,
1066 pChunk,
1067 (bNoDataMaskFullyOpaque) ? NULL : pabyChunkNoDataMask,
1068 nChunkXOffQueried - ((bHasXOffVirtual) ? 0 : nXOff),
1069 nChunkXSizeQueried,
1070 nChunkYOffQueried - ((bHasYOffVirtual) ? 0 : nYOff),
1071 nChunkYSizeQueried,
1072 nDstXOff + nDestXOffVirtual,
1073 nDstXOff + nDestXOffVirtual + nDstXCount,
1074 nDstYOff + nDestYOffVirtual,
1075 nDstYOff + nDestYOffVirtual + nDstYCount,
1076 (GDALRasterBand *) hMEMBand,
1077 pszResampling,
1078 bHasNoData, fNoDataValue,
1079 GetColorTable(),
1080 eDataType );
1081 }
1082
1083 nBlocksDone ++;
1084 if( eErr == CE_None && psExtraArg->pfnProgress != NULL &&
1085 !psExtraArg->pfnProgress(1.0 * nBlocksDone / nTotalBlocks, "",
1086 psExtraArg->pProgressData) )
1087 {
1088 eErr = CE_Failure;
1089 }
1090 }
1091 }
1092
1093 CPLFree(pChunk);
1094 CPLFree(pabyChunkNoDataMask);
1095 }
1096
1097 GDALClose(poMEMDS);
1098
1099 return eErr;
1100 }
1101
1102 /************************************************************************/
1103 /* GDALSwapWords() */
1104 /************************************************************************/
1105
1106 /**
1107 * Byte swap words in-place.
1108 *
1109 * This function will byte swap a set of 2, 4 or 8 byte words "in place" in
1110 * a memory array. No assumption is made that the words being swapped are
1111 * word aligned in memory. Use the CPL_LSB and CPL_MSB macros from cpl_port.h
1112 * to determine if the current platform is big endian or little endian. Use
1113 * The macros like CPL_SWAP32() to byte swap single values without the overhead
1114 * of a function call.
1115 *
1116 * @param pData pointer to start of data buffer.
1117 * @param nWordSize size of words being swapped in bytes. Normally 2, 4 or 8.
1118 * @param nWordCount the number of words to be swapped in this call.
1119 * @param nWordSkip the byte offset from the start of one word to the start of
1120 * the next. For packed buffers this is the same as nWordSize.
1121 */
1122
GDALSwapWords(void * pData,int nWordSize,int nWordCount,int nWordSkip)1123 void CPL_STDCALL GDALSwapWords( void *pData, int nWordSize, int nWordCount,
1124 int nWordSkip )
1125
1126 {
1127 if (nWordCount > 0)
1128 VALIDATE_POINTER0( pData , "GDALSwapWords" );
1129
1130 int i;
1131 GByte *pabyData = (GByte *) pData;
1132
1133 switch( nWordSize )
1134 {
1135 case 1:
1136 break;
1137
1138 case 2:
1139 CPLAssert( nWordSkip >= 2 || nWordCount == 1 );
1140 for( i = 0; i < nWordCount; i++ )
1141 {
1142 GByte byTemp;
1143
1144 byTemp = pabyData[0];
1145 pabyData[0] = pabyData[1];
1146 pabyData[1] = byTemp;
1147
1148 pabyData += nWordSkip;
1149 }
1150 break;
1151
1152 case 4:
1153 CPLAssert( nWordSkip >= 4 || nWordCount == 1 );
1154 for( i = 0; i < nWordCount; i++ )
1155 {
1156 GByte byTemp;
1157
1158 byTemp = pabyData[0];
1159 pabyData[0] = pabyData[3];
1160 pabyData[3] = byTemp;
1161
1162 byTemp = pabyData[1];
1163 pabyData[1] = pabyData[2];
1164 pabyData[2] = byTemp;
1165
1166 pabyData += nWordSkip;
1167 }
1168 break;
1169
1170 case 8:
1171 CPLAssert( nWordSkip >= 8 || nWordCount == 1 );
1172 for( i = 0; i < nWordCount; i++ )
1173 {
1174 GByte byTemp;
1175
1176 byTemp = pabyData[0];
1177 pabyData[0] = pabyData[7];
1178 pabyData[7] = byTemp;
1179
1180 byTemp = pabyData[1];
1181 pabyData[1] = pabyData[6];
1182 pabyData[6] = byTemp;
1183
1184 byTemp = pabyData[2];
1185 pabyData[2] = pabyData[5];
1186 pabyData[5] = byTemp;
1187
1188 byTemp = pabyData[3];
1189 pabyData[3] = pabyData[4];
1190 pabyData[4] = byTemp;
1191
1192 pabyData += nWordSkip;
1193 }
1194 break;
1195
1196 default:
1197 CPLAssert( FALSE );
1198 }
1199 }
1200
1201 #ifdef USE_NEW_COPYWORDS
1202 // Place the new GDALCopyWords helpers in an anonymous namespace
1203 namespace {
1204 /************************************************************************/
1205 /* GetDataLimits() */
1206 /************************************************************************/
1207 /**
1208 * Compute the limits of values that can be placed in Tout in terms of
1209 * Tin. Usually used for output clamping, when the output data type's
1210 * limits are stable relative to the input type (i.e. no roundoff error).
1211 *
1212 * @param tMaxValue the returned maximum value
1213 * @param tMinValue the returned minimum value
1214 */
1215
1216 template <class Tin, class Tout>
GetDataLimits(Tin & tMaxValue,Tin & tMinValue)1217 inline void GetDataLimits(Tin &tMaxValue, Tin &tMinValue)
1218 {
1219 tMaxValue = std::numeric_limits<Tin>::max();
1220 tMinValue = std::numeric_limits<Tin>::min();
1221
1222 // Compute the actual minimum value of Tout in terms of Tin.
1223 if (std::numeric_limits<Tout>::is_signed && std::numeric_limits<Tout>::is_integer)
1224 {
1225 // the minimum value is less than zero
1226 if (std::numeric_limits<Tout>::digits < std::numeric_limits<Tin>::digits ||
1227 !std::numeric_limits<Tin>::is_integer)
1228 {
1229 // Tout is smaller than Tin, so we need to clamp values in input
1230 // to the range of Tout's min/max values
1231 if (std::numeric_limits<Tin>::is_signed)
1232 {
1233 tMinValue = static_cast<Tin>(std::numeric_limits<Tout>::min());
1234 }
1235 tMaxValue = static_cast<Tin>(std::numeric_limits<Tout>::max());
1236 }
1237 }
1238 else if (std::numeric_limits<Tout>::is_integer)
1239 {
1240 // the output is unsigned, so we just need to determine the max
1241 if (std::numeric_limits<Tout>::digits <= std::numeric_limits<Tin>::digits)
1242 {
1243 // Tout is smaller than Tin, so we need to clamp the input values
1244 // to the range of Tout's max
1245 tMaxValue = static_cast<Tin>(std::numeric_limits<Tout>::max());
1246 }
1247 tMinValue = 0;
1248 }
1249
1250 }
1251
1252 /************************************************************************/
1253 /* ClampValue() */
1254 /************************************************************************/
1255 /**
1256 * Clamp values of type T to a specified range
1257 *
1258 * @param tValue the value
1259 * @param tMax the max value
1260 * @param tMin the min value
1261 */
1262 template <class T>
ClampValue(const T tValue,const T tMax,const T tMin)1263 inline T ClampValue(const T tValue, const T tMax, const T tMin)
1264 {
1265 return tValue > tMax ? tMax :
1266 tValue < tMin ? tMin : tValue;
1267 }
1268
1269 /************************************************************************/
1270 /* CopyWord() */
1271 /************************************************************************/
1272 /**
1273 * Copy a single word, optionally rounding if appropriate (i.e. going
1274 * from the float to the integer case). Note that this is the function
1275 * you should specialize if you're adding a new data type.
1276 *
1277 * @param tValueIn value of type Tin; the input value to be converted
1278 * @param tValueOut value of type Tout; the output value
1279 */
1280
1281 template <class Tin, class Tout>
CopyWord(const Tin tValueIn,Tout & tValueOut)1282 inline void CopyWord(const Tin tValueIn, Tout &tValueOut)
1283 {
1284 Tin tMaxVal, tMinVal;
1285 GetDataLimits<Tin, Tout>(tMaxVal, tMinVal);
1286 tValueOut = static_cast<Tout>(ClampValue(tValueIn, tMaxVal, tMinVal));
1287 }
1288
1289 template <class Tin>
CopyWord(const Tin tValueIn,float & fValueOut)1290 inline void CopyWord(const Tin tValueIn, float &fValueOut)
1291 {
1292 fValueOut = (float) tValueIn;
1293 }
1294
1295 template <class Tin>
CopyWord(const Tin tValueIn,double & dfValueOut)1296 inline void CopyWord(const Tin tValueIn, double &dfValueOut)
1297 {
1298 dfValueOut = tValueIn;
1299 }
1300
CopyWord(const double dfValueIn,double & dfValueOut)1301 inline void CopyWord(const double dfValueIn, double &dfValueOut)
1302 {
1303 dfValueOut = dfValueIn;
1304 }
1305
CopyWord(const float fValueIn,float & fValueOut)1306 inline void CopyWord(const float fValueIn, float &fValueOut)
1307 {
1308 fValueOut = fValueIn;
1309 }
1310
CopyWord(const float fValueIn,double & dfValueOut)1311 inline void CopyWord(const float fValueIn, double &dfValueOut)
1312 {
1313 dfValueOut = fValueIn;
1314 }
1315
CopyWord(const double dfValueIn,float & fValueOut)1316 inline void CopyWord(const double dfValueIn, float &fValueOut)
1317 {
1318 fValueOut = static_cast<float>(dfValueIn);
1319 }
1320
1321 template <class Tout>
CopyWord(const float fValueIn,Tout & tValueOut)1322 inline void CopyWord(const float fValueIn, Tout &tValueOut)
1323 {
1324 float fMaxVal, fMinVal;
1325 GetDataLimits<float, Tout>(fMaxVal, fMinVal);
1326 tValueOut = static_cast<Tout>(
1327 ClampValue(fValueIn + 0.5f, fMaxVal, fMinVal));
1328 }
1329
1330 template <class Tout>
CopyWord(const double dfValueIn,Tout & tValueOut)1331 inline void CopyWord(const double dfValueIn, Tout &tValueOut)
1332 {
1333 double dfMaxVal, dfMinVal;
1334 GetDataLimits<double, Tout>(dfMaxVal, dfMinVal);
1335 tValueOut = static_cast<Tout>(
1336 ClampValue(dfValueIn + 0.5, dfMaxVal, dfMinVal));
1337 }
1338
CopyWord(const double dfValueIn,int & nValueOut)1339 inline void CopyWord(const double dfValueIn, int &nValueOut)
1340 {
1341 double dfMaxVal, dfMinVal;
1342 GetDataLimits<double, int>(dfMaxVal, dfMinVal);
1343 double dfValue = dfValueIn >= 0.0 ? dfValueIn + 0.5 :
1344 dfValueIn - 0.5;
1345 nValueOut = static_cast<int>(
1346 ClampValue(dfValue, dfMaxVal, dfMinVal));
1347 }
1348
CopyWord(const float fValueIn,short & nValueOut)1349 inline void CopyWord(const float fValueIn, short &nValueOut)
1350 {
1351 float fMaxVal, fMinVal;
1352 GetDataLimits<float, short>(fMaxVal, fMinVal);
1353 float fValue = fValueIn >= 0.0f ? fValueIn + 0.5f :
1354 fValueIn - 0.5f;
1355 nValueOut = static_cast<short>(
1356 ClampValue(fValue, fMaxVal, fMinVal));
1357 }
1358
CopyWord(const double dfValueIn,short & nValueOut)1359 inline void CopyWord(const double dfValueIn, short &nValueOut)
1360 {
1361 double dfMaxVal, dfMinVal;
1362 GetDataLimits<double, short>(dfMaxVal, dfMinVal);
1363 double dfValue = dfValueIn > 0.0 ? dfValueIn + 0.5 :
1364 dfValueIn - 0.5;
1365 nValueOut = static_cast<short>(
1366 ClampValue(dfValue, dfMaxVal, dfMinVal));
1367 }
1368
1369 // Roundoff occurs for Float32 -> int32 for max/min. Overload CopyWord
1370 // specifically for this case.
CopyWord(const float fValueIn,int & nValueOut)1371 inline void CopyWord(const float fValueIn, int &nValueOut)
1372 {
1373 if (fValueIn >= static_cast<float>(std::numeric_limits<int>::max()))
1374 {
1375 nValueOut = std::numeric_limits<int>::max();
1376 }
1377 else if (fValueIn <= static_cast<float>(std::numeric_limits<int>::min()))
1378 {
1379 nValueOut = std::numeric_limits<int>::min();
1380 }
1381 else
1382 {
1383 nValueOut = static_cast<int>(fValueIn > 0.0f ?
1384 fValueIn + 0.5f : fValueIn - 0.5f);
1385 }
1386 }
1387
1388 // Roundoff occurs for Float32 -> uint32 for max. Overload CopyWord
1389 // specifically for this case.
CopyWord(const float fValueIn,unsigned int & nValueOut)1390 inline void CopyWord(const float fValueIn, unsigned int &nValueOut)
1391 {
1392 if (fValueIn >= static_cast<float>(std::numeric_limits<unsigned int>::max()))
1393 {
1394 nValueOut = std::numeric_limits<unsigned int>::max();
1395 }
1396 else if (fValueIn <= static_cast<float>(std::numeric_limits<unsigned int>::min()))
1397 {
1398 nValueOut = std::numeric_limits<unsigned int>::min();
1399 }
1400 else
1401 {
1402 nValueOut = static_cast<unsigned int>(fValueIn + 0.5f);
1403 }
1404 }
1405
1406 /************************************************************************/
1407 /* GDALCopyWordsT() */
1408 /************************************************************************/
1409 /**
1410 * Template function, used to copy data from pSrcData into buffer
1411 * pDstData, with stride nSrcPixelStride in the source data and
1412 * stride nDstPixelStride in the destination data. This template can
1413 * deal with the case where the input data type is real or complex and
1414 * the output is real.
1415 *
1416 * @param pSrcData the source data buffer
1417 * @param nSrcPixelStride the stride, in the buffer pSrcData for pixels
1418 * of interest.
1419 * @param pDstData the destination buffer.
1420 * @param nDstPixelStride the stride in the buffer pDstData for pixels of
1421 * interest.
1422 * @param nWordCount the total number of pixel words to copy
1423 *
1424 * @code
1425 * // Assume an input buffer of type GUInt16 named pBufferIn
1426 * GByte *pBufferOut = new GByte[numBytesOut];
1427 * GDALCopyWordsT<GUInt16, GByte>(pSrcData, 2, pDstData, 1, numBytesOut);
1428 * @code
1429 * @note
1430 * This is a private function, and should not be exposed outside of rasterio.cpp.
1431 * External users should call the GDALCopyWords driver function.
1432 * @note
1433 */
1434
1435 template <class Tin, class Tout>
GDALCopyWordsT(const Tin * const pSrcData,int nSrcPixelStride,Tout * const pDstData,int nDstPixelStride,int nWordCount)1436 static void GDALCopyWordsT(const Tin* const pSrcData, int nSrcPixelStride,
1437 Tout* const pDstData, int nDstPixelStride,
1438 int nWordCount)
1439 {
1440 std::ptrdiff_t nDstOffset = 0;
1441
1442 const char* const pSrcDataPtr = reinterpret_cast<const char*>(pSrcData);
1443 char* const pDstDataPtr = reinterpret_cast<char*>(pDstData);
1444 for (std::ptrdiff_t n = 0; n < nWordCount; n++)
1445 {
1446 const Tin tValue = *reinterpret_cast<const Tin*>(pSrcDataPtr + (n * nSrcPixelStride));
1447 Tout* const pOutPixel = reinterpret_cast<Tout*>(pDstDataPtr + nDstOffset);
1448
1449 CopyWord(tValue, *pOutPixel);
1450
1451 nDstOffset += nDstPixelStride;
1452 }
1453 }
1454
1455 /************************************************************************/
1456 /* GDALCopyWordsComplexT() */
1457 /************************************************************************/
1458 /**
1459 * Template function, used to copy data from pSrcData into buffer
1460 * pDstData, with stride nSrcPixelStride in the source data and
1461 * stride nDstPixelStride in the destination data. Deals with the
1462 * complex case, where input is complex and output is complex.
1463 *
1464 * @param pSrcData the source data buffer
1465 * @param nSrcPixelStride the stride, in the buffer pSrcData for pixels
1466 * of interest.
1467 * @param pDstData the destination buffer.
1468 * @param nDstPixelStride the stride in the buffer pDstData for pixels of
1469 * interest.
1470 * @param nWordCount the total number of pixel words to copy
1471 *
1472 */
1473 template <class Tin, class Tout>
GDALCopyWordsComplexT(const Tin * const pSrcData,int nSrcPixelStride,Tout * const pDstData,int nDstPixelStride,int nWordCount)1474 inline void GDALCopyWordsComplexT(const Tin* const pSrcData, int nSrcPixelStride,
1475 Tout* const pDstData, int nDstPixelStride,
1476 int nWordCount)
1477 {
1478 std::ptrdiff_t nDstOffset = 0;
1479 const char* const pSrcDataPtr = reinterpret_cast<const char*>(pSrcData);
1480 char* const pDstDataPtr = reinterpret_cast<char*>(pDstData);
1481
1482 // Determine the minimum and maximum value we can have based
1483 // on the constraints of Tin and Tout.
1484 Tin tMaxValue, tMinValue;
1485 GetDataLimits<Tin, Tout>(tMaxValue, tMinValue);
1486
1487 for (std::ptrdiff_t n = 0; n < nWordCount; n++)
1488 {
1489 const Tin* const pPixelIn = reinterpret_cast<const Tin*>(pSrcDataPtr + n * nSrcPixelStride);
1490 Tout* const pPixelOut = reinterpret_cast<Tout*>(pDstDataPtr + nDstOffset);
1491
1492 CopyWord(pPixelIn[0], pPixelOut[0]);
1493 CopyWord(pPixelIn[1], pPixelOut[1]);
1494
1495 nDstOffset += nDstPixelStride;
1496 }
1497 }
1498
1499 /************************************************************************/
1500 /* GDALCopyWordsComplexOutT() */
1501 /************************************************************************/
1502 /**
1503 * Template function, used to copy data from pSrcData into buffer
1504 * pDstData, with stride nSrcPixelStride in the source data and
1505 * stride nDstPixelStride in the destination data. Deals with the
1506 * case where the value is real coming in, but complex going out.
1507 *
1508 * @param pSrcData the source data buffer
1509 * @param nSrcPixelStride the stride, in the buffer pSrcData for pixels
1510 * of interest, in bytes.
1511 * @param pDstData the destination buffer.
1512 * @param nDstPixelStride the stride in the buffer pDstData for pixels of
1513 * interest, in bytes.
1514 * @param nWordCount the total number of pixel words to copy
1515 *
1516 */
1517 template <class Tin, class Tout>
GDALCopyWordsComplexOutT(const Tin * const pSrcData,int nSrcPixelStride,Tout * const pDstData,int nDstPixelStride,int nWordCount)1518 inline void GDALCopyWordsComplexOutT(const Tin* const pSrcData, int nSrcPixelStride,
1519 Tout* const pDstData, int nDstPixelStride,
1520 int nWordCount)
1521 {
1522 std::ptrdiff_t nDstOffset = 0;
1523
1524 const Tout tOutZero = static_cast<Tout>(0);
1525
1526 const char* const pSrcDataPtr = reinterpret_cast<const char*>(pSrcData);
1527 char* const pDstDataPtr = reinterpret_cast<char*>(pDstData);
1528
1529 for (std::ptrdiff_t n = 0; n < nWordCount; n++)
1530 {
1531 const Tin tValue = *reinterpret_cast<const Tin* const>(pSrcDataPtr + n * nSrcPixelStride);
1532 Tout* const pPixelOut = reinterpret_cast<Tout*>(pDstDataPtr + nDstOffset);
1533 CopyWord(tValue, *pPixelOut);
1534
1535 pPixelOut[1] = tOutZero;
1536
1537 nDstOffset += nDstPixelStride;
1538 }
1539 }
1540
1541 /************************************************************************/
1542 /* GDALCopyWordsFromT() */
1543 /************************************************************************/
1544 /**
1545 * Template driver function. Given the input type T, call the appropriate
1546 * GDALCopyWordsT function template for the desired output type. You should
1547 * never call this function directly (call GDALCopyWords instead).
1548 *
1549 * @param pSrcData source data buffer
1550 * @param nSrcPixelStride pixel stride in input buffer, in pixel words
1551 * @param bInComplex input is complex
1552 * @param pDstData destination data buffer
1553 * @param eDstType destination data type
1554 * @param nDstPixelStride pixel stride in output buffer, in pixel words
1555 * @param nWordCount number of pixel words to be copied
1556 */
1557 template <class T>
GDALCopyWordsFromT(const T * const pSrcData,int nSrcPixelStride,bool bInComplex,void * pDstData,GDALDataType eDstType,int nDstPixelStride,int nWordCount)1558 inline void GDALCopyWordsFromT(const T* const pSrcData, int nSrcPixelStride, bool bInComplex,
1559 void *pDstData, GDALDataType eDstType, int nDstPixelStride,
1560 int nWordCount)
1561 {
1562 switch (eDstType)
1563 {
1564 case GDT_Byte:
1565 GDALCopyWordsT(pSrcData, nSrcPixelStride,
1566 static_cast<unsigned char*>(pDstData), nDstPixelStride,
1567 nWordCount);
1568 break;
1569 case GDT_UInt16:
1570 GDALCopyWordsT(pSrcData, nSrcPixelStride,
1571 static_cast<unsigned short*>(pDstData), nDstPixelStride,
1572 nWordCount);
1573 break;
1574 case GDT_Int16:
1575 GDALCopyWordsT(pSrcData, nSrcPixelStride,
1576 static_cast<short*>(pDstData), nDstPixelStride,
1577 nWordCount);
1578 break;
1579 case GDT_UInt32:
1580 GDALCopyWordsT(pSrcData, nSrcPixelStride,
1581 static_cast<unsigned int*>(pDstData), nDstPixelStride,
1582 nWordCount);
1583 break;
1584 case GDT_Int32:
1585 GDALCopyWordsT(pSrcData, nSrcPixelStride,
1586 static_cast<int*>(pDstData), nDstPixelStride,
1587 nWordCount);
1588 break;
1589 case GDT_Float32:
1590 GDALCopyWordsT(pSrcData, nSrcPixelStride,
1591 static_cast<float*>(pDstData), nDstPixelStride,
1592 nWordCount);
1593 break;
1594 case GDT_Float64:
1595 GDALCopyWordsT(pSrcData, nSrcPixelStride,
1596 static_cast<double*>(pDstData), nDstPixelStride,
1597 nWordCount);
1598 break;
1599 case GDT_CInt16:
1600 if (bInComplex)
1601 {
1602 GDALCopyWordsComplexT(pSrcData, nSrcPixelStride,
1603 static_cast<short *>(pDstData), nDstPixelStride,
1604 nWordCount);
1605 }
1606 else // input is not complex, so we need to promote to a complex buffer
1607 {
1608 GDALCopyWordsComplexOutT(pSrcData, nSrcPixelStride,
1609 static_cast<short *>(pDstData), nDstPixelStride,
1610 nWordCount);
1611 }
1612 break;
1613 case GDT_CInt32:
1614 if (bInComplex)
1615 {
1616 GDALCopyWordsComplexT(pSrcData, nSrcPixelStride,
1617 static_cast<int *>(pDstData), nDstPixelStride,
1618 nWordCount);
1619 }
1620 else // input is not complex, so we need to promote to a complex buffer
1621 {
1622 GDALCopyWordsComplexOutT(pSrcData, nSrcPixelStride,
1623 static_cast<int *>(pDstData), nDstPixelStride,
1624 nWordCount);
1625 }
1626 break;
1627 case GDT_CFloat32:
1628 if (bInComplex)
1629 {
1630 GDALCopyWordsComplexT(pSrcData, nSrcPixelStride,
1631 static_cast<float *>(pDstData), nDstPixelStride,
1632 nWordCount);
1633 }
1634 else // input is not complex, so we need to promote to a complex buffer
1635 {
1636 GDALCopyWordsComplexOutT(pSrcData, nSrcPixelStride,
1637 static_cast<float *>(pDstData), nDstPixelStride,
1638 nWordCount);
1639 }
1640 break;
1641 case GDT_CFloat64:
1642 if (bInComplex)
1643 {
1644 GDALCopyWordsComplexT(pSrcData, nSrcPixelStride,
1645 static_cast<double *>(pDstData), nDstPixelStride,
1646 nWordCount);
1647 }
1648 else // input is not complex, so we need to promote to a complex buffer
1649 {
1650 GDALCopyWordsComplexOutT(pSrcData, nSrcPixelStride,
1651 static_cast<double *>(pDstData), nDstPixelStride,
1652 nWordCount);
1653 }
1654 break;
1655 case GDT_Unknown:
1656 default:
1657 CPLAssert(FALSE);
1658 }
1659 }
1660
1661 } // end anonymous namespace
1662 #endif
1663
1664 /************************************************************************/
1665 /* GDALReplicateWord() */
1666 /************************************************************************/
1667
GDALReplicateWord(void * pSrcData,GDALDataType eSrcType,void * pDstData,GDALDataType eDstType,int nDstPixelStride,int nWordCount)1668 void GDALReplicateWord(void *pSrcData, GDALDataType eSrcType,
1669 void *pDstData, GDALDataType eDstType, int nDstPixelStride,
1670 int nWordCount)
1671 {
1672 /* ----------------------------------------------------------------------- */
1673 /* Special case when the source data is always the same value */
1674 /* (for VRTSourcedRasterBand::IRasterIO and VRTDerivedRasterBand::IRasterIO*/
1675 /* for example) */
1676 /* ----------------------------------------------------------------------- */
1677 /* Let the general translation case do the necessary conversions */
1678 /* on the first destination element */
1679 GDALCopyWords(pSrcData, eSrcType, 0,
1680 pDstData, eDstType, nDstPixelStride,
1681 1 );
1682
1683 /* Now copy the first element to the nWordCount - 1 following destination */
1684 /* elements */
1685 nWordCount--;
1686 GByte *pabyDstWord = ((GByte *)pDstData) + nDstPixelStride;
1687
1688 switch (eDstType)
1689 {
1690 case GDT_Byte:
1691 {
1692 if (nDstPixelStride == 1)
1693 {
1694 if (nWordCount > 0)
1695 memset(pabyDstWord, *(GByte*)pDstData, nWordCount);
1696 }
1697 else
1698 {
1699 GByte valSet = *(GByte*)pDstData;
1700 while(nWordCount--)
1701 {
1702 *pabyDstWord = valSet;
1703 pabyDstWord += nDstPixelStride;
1704 }
1705 }
1706 break;
1707 }
1708
1709 #define CASE_DUPLICATE_SIMPLE(enum_type, c_type) \
1710 case enum_type:\
1711 { \
1712 c_type valSet = *(c_type*)pDstData; \
1713 while(nWordCount--) \
1714 { \
1715 *(c_type*)pabyDstWord = valSet; \
1716 pabyDstWord += nDstPixelStride; \
1717 } \
1718 break; \
1719 }
1720
1721 CASE_DUPLICATE_SIMPLE(GDT_UInt16, GUInt16)
1722 CASE_DUPLICATE_SIMPLE(GDT_Int16, GInt16)
1723 CASE_DUPLICATE_SIMPLE(GDT_UInt32, GUInt32)
1724 CASE_DUPLICATE_SIMPLE(GDT_Int32, GInt32)
1725 CASE_DUPLICATE_SIMPLE(GDT_Float32,float)
1726 CASE_DUPLICATE_SIMPLE(GDT_Float64,double)
1727
1728 #define CASE_DUPLICATE_COMPLEX(enum_type, c_type) \
1729 case enum_type:\
1730 { \
1731 c_type valSet1 = ((c_type*)pDstData)[0]; \
1732 c_type valSet2 = ((c_type*)pDstData)[1]; \
1733 while(nWordCount--) \
1734 { \
1735 ((c_type*)pabyDstWord)[0] = valSet1; \
1736 ((c_type*)pabyDstWord)[1] = valSet2; \
1737 pabyDstWord += nDstPixelStride; \
1738 } \
1739 break; \
1740 }
1741
1742 CASE_DUPLICATE_COMPLEX(GDT_CInt16, GInt16)
1743 CASE_DUPLICATE_COMPLEX(GDT_CInt32, GInt32)
1744 CASE_DUPLICATE_COMPLEX(GDT_CFloat32, float)
1745 CASE_DUPLICATE_COMPLEX(GDT_CFloat64, double)
1746
1747 default:
1748 CPLAssert( FALSE );
1749 }
1750 }
1751
1752 /************************************************************************/
1753 /* GDALCopyWords() */
1754 /************************************************************************/
1755
1756 /**
1757 * Copy pixel words from buffer to buffer.
1758 *
1759 * This function is used to copy pixel word values from one memory buffer
1760 * to another, with support for conversion between data types, and differing
1761 * step factors. The data type conversion is done using the normal GDAL
1762 * rules. Values assigned to a lower range integer type are clipped. For
1763 * instance assigning GDT_Int16 values to a GDT_Byte buffer will cause values
1764 * less the 0 to be set to 0, and values larger than 255 to be set to 255.
1765 * Assignment from floating point to integer uses default C type casting
1766 * semantics. Assignment from non-complex to complex will result in the
1767 * imaginary part being set to zero on output. Assigment from complex to
1768 * non-complex will result in the complex portion being lost and the real
1769 * component being preserved (<i>not magnitidue!</i>).
1770 *
1771 * No assumptions are made about the source or destination words occuring
1772 * on word boundaries. It is assumed that all values are in native machine
1773 * byte order.
1774 *
1775 * @param pSrcData Pointer to source data to be converted.
1776 * @param eSrcType the source data type (see GDALDataType enum)
1777 * @param nSrcPixelStride Source pixel stride (i.e. distance between 2 words), in bytes
1778 * @param pDstData Pointer to buffer where destination data should go
1779 * @param eDstType the destination data type (see GDALDataType enum)
1780 * @param nDstPixelStride Destination pixel stride (i.e. distance between 2 words), in bytes
1781 * @param nWordCount number of words to be copied
1782 *
1783 *
1784 * @note
1785 * When adding a new data type to GDAL, you must do the following to
1786 * support it properly within the GDALCopyWords function:
1787 * 1. Add the data type to the switch on eSrcType in GDALCopyWords.
1788 * This should invoke the appropriate GDALCopyWordsFromT wrapper.
1789 * 2. Add the data type to the switch on eDstType in GDALCopyWordsFromT.
1790 * This should call the appropriate GDALCopyWordsT template.
1791 * 3. If appropriate, overload the appropriate CopyWord template in the
1792 * above namespace. This will ensure that any conversion issues are
1793 * handled (cases like the float -> int32 case, where the min/max)
1794 * values are subject to roundoff error.
1795 */
1796
1797 void CPL_STDCALL
GDALCopyWords(void * pSrcData,GDALDataType eSrcType,int nSrcPixelStride,void * pDstData,GDALDataType eDstType,int nDstPixelStride,int nWordCount)1798 GDALCopyWords( void * pSrcData, GDALDataType eSrcType, int nSrcPixelStride,
1799 void * pDstData, GDALDataType eDstType, int nDstPixelStride,
1800 int nWordCount )
1801
1802 {
1803 // Deal with the case where we're replicating a single word into the
1804 // provided buffer
1805 if (nSrcPixelStride == 0 && nWordCount > 1)
1806 {
1807 GDALReplicateWord(pSrcData, eSrcType, pDstData, eDstType, nDstPixelStride, nWordCount);
1808 return;
1809 }
1810
1811 #ifdef USE_NEW_COPYWORDS
1812
1813 int nSrcDataTypeSize = GDALGetDataTypeSize(eSrcType) / 8;
1814 // Let memcpy() handle the case where we're copying a packed buffer
1815 // of pixels.
1816 if (eSrcType == eDstType && nSrcPixelStride == nDstPixelStride &&
1817 nSrcPixelStride == nSrcDataTypeSize)
1818 {
1819 memcpy(pDstData, pSrcData, nWordCount * nSrcDataTypeSize);
1820 return;
1821 }
1822
1823 // Handle the more general case -- deals with conversion of data types
1824 // directly.
1825 switch (eSrcType)
1826 {
1827 case GDT_Byte:
1828 GDALCopyWordsFromT<unsigned char>(static_cast<unsigned char *>(pSrcData), nSrcPixelStride, false,
1829 pDstData, eDstType, nDstPixelStride,
1830 nWordCount);
1831 break;
1832 case GDT_UInt16:
1833 GDALCopyWordsFromT<unsigned short>(static_cast<unsigned short *>(pSrcData), nSrcPixelStride, false,
1834 pDstData, eDstType, nDstPixelStride,
1835 nWordCount);
1836 break;
1837 case GDT_Int16:
1838 GDALCopyWordsFromT<short>(static_cast<short *>(pSrcData), nSrcPixelStride, false,
1839 pDstData, eDstType, nDstPixelStride,
1840 nWordCount);
1841 break;
1842 case GDT_UInt32:
1843 GDALCopyWordsFromT<unsigned int>(static_cast<unsigned int *>(pSrcData), nSrcPixelStride, false,
1844 pDstData, eDstType, nDstPixelStride,
1845 nWordCount);
1846 break;
1847 case GDT_Int32:
1848 GDALCopyWordsFromT<int>(static_cast<int *>(pSrcData), nSrcPixelStride, false,
1849 pDstData, eDstType, nDstPixelStride,
1850 nWordCount);
1851 break;
1852 case GDT_Float32:
1853 GDALCopyWordsFromT<float>(static_cast<float *>(pSrcData), nSrcPixelStride, false,
1854 pDstData, eDstType, nDstPixelStride,
1855 nWordCount);
1856 break;
1857 case GDT_Float64:
1858 GDALCopyWordsFromT<double>(static_cast<double *>(pSrcData), nSrcPixelStride, false,
1859 pDstData, eDstType, nDstPixelStride,
1860 nWordCount);
1861 break;
1862 case GDT_CInt16:
1863 GDALCopyWordsFromT<short>(static_cast<short *>(pSrcData), nSrcPixelStride, true,
1864 pDstData, eDstType, nDstPixelStride,
1865 nWordCount);
1866 break;
1867 case GDT_CInt32:
1868 GDALCopyWordsFromT<int>(static_cast<int *>(pSrcData), nSrcPixelStride, true,
1869 pDstData, eDstType, nDstPixelStride,
1870 nWordCount);
1871 break;
1872 case GDT_CFloat32:
1873 GDALCopyWordsFromT<float>(static_cast<float *>(pSrcData), nSrcPixelStride, true,
1874 pDstData, eDstType, nDstPixelStride,
1875 nWordCount);
1876 break;
1877 case GDT_CFloat64:
1878 GDALCopyWordsFromT<double>(static_cast<double *>(pSrcData), nSrcPixelStride, true,
1879 pDstData, eDstType, nDstPixelStride,
1880 nWordCount);
1881 break;
1882 case GDT_Unknown:
1883 default:
1884 CPLAssert(FALSE);
1885 }
1886
1887 #else // undefined USE_NEW_COPYWORDS
1888 /* -------------------------------------------------------------------- */
1889 /* Special case when no data type translation is required. */
1890 /* -------------------------------------------------------------------- */
1891 if( eSrcType == eDstType )
1892 {
1893 int nWordSize = GDALGetDataTypeSize(eSrcType)/8;
1894 int i;
1895
1896 // contiguous blocks.
1897 if( nWordSize == nSrcPixelStride && nWordSize == nDstPixelStride )
1898 {
1899 memcpy( pDstData, pSrcData, nSrcPixelStride * nWordCount );
1900 return;
1901 }
1902
1903 GByte *pabySrc = (GByte *) pSrcData;
1904 GByte *pabyDst = (GByte *) pDstData;
1905
1906 // Moving single bytes.
1907 if( nWordSize == 1 )
1908 {
1909 if (nWordCount > 100)
1910 {
1911 /* ==================================================================== */
1912 /* Optimization for high number of words to transfer and some */
1913 /* typical source and destination pixel spacing : we unroll the */
1914 /* loop. */
1915 /* ==================================================================== */
1916 #define ASSIGN(_nSrcPixelStride, _nDstPixelStride, _k) \
1917 pabyDst[_nDstPixelStride * _k] = pabySrc[_nSrcPixelStride * _k]
1918 #define ASSIGN_LOOP(_nSrcPixelStride, _nDstPixelStride) \
1919 for( i = nWordCount / 16 ; i != 0; i-- ) \
1920 { \
1921 ASSIGN(_nSrcPixelStride, _nDstPixelStride, 0); \
1922 ASSIGN(_nSrcPixelStride, _nDstPixelStride, 1); \
1923 ASSIGN(_nSrcPixelStride, _nDstPixelStride, 2); \
1924 ASSIGN(_nSrcPixelStride, _nDstPixelStride, 3); \
1925 ASSIGN(_nSrcPixelStride, _nDstPixelStride, 4); \
1926 ASSIGN(_nSrcPixelStride, _nDstPixelStride, 5); \
1927 ASSIGN(_nSrcPixelStride, _nDstPixelStride, 6); \
1928 ASSIGN(_nSrcPixelStride, _nDstPixelStride, 7); \
1929 ASSIGN(_nSrcPixelStride, _nDstPixelStride, 8); \
1930 ASSIGN(_nSrcPixelStride, _nDstPixelStride, 9); \
1931 ASSIGN(_nSrcPixelStride, _nDstPixelStride, 10); \
1932 ASSIGN(_nSrcPixelStride, _nDstPixelStride, 11); \
1933 ASSIGN(_nSrcPixelStride, _nDstPixelStride, 12); \
1934 ASSIGN(_nSrcPixelStride, _nDstPixelStride, 13); \
1935 ASSIGN(_nSrcPixelStride, _nDstPixelStride, 14); \
1936 ASSIGN(_nSrcPixelStride, _nDstPixelStride, 15); \
1937 pabyDst += _nDstPixelStride * 16; \
1938 pabySrc += _nSrcPixelStride * 16; \
1939 } \
1940 nWordCount = nWordCount % 16;
1941
1942 if (nSrcPixelStride == 3 && nDstPixelStride == 1)
1943 {
1944 ASSIGN_LOOP(3, 1)
1945 }
1946 else if (nSrcPixelStride == 1 && nDstPixelStride == 3)
1947 {
1948 ASSIGN_LOOP(1, 3)
1949 }
1950 else if (nSrcPixelStride == 4 && nDstPixelStride == 1)
1951 {
1952 ASSIGN_LOOP(4, 1)
1953 }
1954 else if (nSrcPixelStride == 1 && nDstPixelStride == 4)
1955 {
1956 ASSIGN_LOOP(1, 4)
1957 }
1958 else if (nSrcPixelStride == 3 && nDstPixelStride == 4)
1959 {
1960 ASSIGN_LOOP(3, 4)
1961 }
1962 else if (nSrcPixelStride == 4 && nDstPixelStride == 3)
1963 {
1964 ASSIGN_LOOP(4, 3)
1965 }
1966 }
1967
1968 for( i = nWordCount; i != 0; i-- )
1969 {
1970 *pabyDst = *pabySrc;
1971 pabyDst += nDstPixelStride;
1972 pabySrc += nSrcPixelStride;
1973 }
1974 }
1975 else if (nWordSize == 2)
1976 {
1977 for( i = nWordCount; i != 0; i-- )
1978 {
1979 *(short*)pabyDst = *(short*)pabySrc;
1980 pabyDst += nDstPixelStride;
1981 pabySrc += nSrcPixelStride;
1982 }
1983 }
1984 else if (nWordSize == 4)
1985 {
1986 for( i = nWordCount; i != 0; i-- )
1987 {
1988 *(int*)pabyDst = *(int*)pabySrc;
1989 pabyDst += nDstPixelStride;
1990 pabySrc += nSrcPixelStride;
1991 }
1992 }
1993 else if (nWordSize == 8)
1994 {
1995 for( i = nWordCount; i != 0; i-- )
1996 {
1997 ((int*)pabyDst)[0] = ((int*)pabySrc)[0];
1998 ((int*)pabyDst)[1] = ((int*)pabySrc)[1];
1999 pabyDst += nDstPixelStride;
2000 pabySrc += nSrcPixelStride;
2001 }
2002 }
2003 else if (nWordSize == 16)
2004 {
2005 for( i = nWordCount; i != 0; i-- )
2006 {
2007 ((int*)pabyDst)[0] = ((int*)pabySrc)[0];
2008 ((int*)pabyDst)[1] = ((int*)pabySrc)[1];
2009 ((int*)pabyDst)[2] = ((int*)pabySrc)[2];
2010 ((int*)pabyDst)[3] = ((int*)pabySrc)[3];
2011 pabyDst += nDstPixelStride;
2012 pabySrc += nSrcPixelStride;
2013 }
2014 }
2015 else
2016 {
2017 CPLAssert(FALSE);
2018 }
2019
2020 return;
2021 }
2022
2023 /* ==================================================================== */
2024 /* General translation case */
2025 /* ==================================================================== */
2026 for( int iWord = 0; iWord < nWordCount; iWord++ )
2027 {
2028 void *pSrcWord, *pDstWord;
2029 double dfPixelValue=0.0, dfPixelValueI=0.0;
2030
2031 pSrcWord = static_cast<GByte *>(pSrcData) + iWord * nSrcPixelStride;
2032 pDstWord = static_cast<GByte *>(pDstData) + iWord * nDstPixelStride;
2033
2034 /* -------------------------------------------------------------------- */
2035 /* Fetch source value based on data type. */
2036 /* -------------------------------------------------------------------- */
2037 switch( eSrcType )
2038 {
2039 case GDT_Byte:
2040 {
2041 GByte byVal = *static_cast<GByte *>(pSrcWord);
2042 switch( eDstType )
2043 {
2044 case GDT_UInt16:
2045 *static_cast<GUInt16 *>(pDstWord) = byVal;
2046 continue;
2047 case GDT_Int16:
2048 *static_cast<GInt16 *>(pDstWord) = byVal;
2049 continue;
2050 case GDT_UInt32:
2051 *static_cast<GUInt32 *>(pDstWord) = byVal;
2052 continue;
2053 case GDT_Int32:
2054 *static_cast<GInt32 *>(pDstWord) = byVal;
2055 continue;
2056 case GDT_CInt16:
2057 {
2058 GInt16 *panDstWord = static_cast<GInt16 *>(pDstWord);
2059 panDstWord[0] = byVal;
2060 panDstWord[1] = 0;
2061 continue;
2062 }
2063 case GDT_CInt32:
2064 {
2065 GInt32 *panDstWord = static_cast<GInt32 *>(pDstWord);
2066 panDstWord[0] = byVal;
2067 panDstWord[1] = 0;
2068 continue;
2069 }
2070 default:
2071 break;
2072 }
2073 dfPixelValue = byVal;
2074 }
2075 break;
2076
2077 case GDT_UInt16:
2078 {
2079 GUInt16 nVal = *static_cast<GUInt16 *>(pSrcWord);
2080 switch( eDstType )
2081 {
2082 case GDT_Byte:
2083 {
2084 GByte byVal;
2085 if( nVal > 255 )
2086 byVal = 255;
2087 else
2088 byVal = static_cast<GByte>(nVal);
2089 *static_cast<GByte *>(pDstWord) = byVal;
2090 continue;
2091 }
2092 case GDT_Int16:
2093 if( nVal > 32767 )
2094 nVal = 32767;
2095 *static_cast<GInt16 *>(pDstWord) = nVal;
2096 continue;
2097 case GDT_UInt32:
2098 *static_cast<GUInt32 *>(pDstWord) = nVal;
2099 continue;
2100 case GDT_Int32:
2101 *static_cast<GInt32 *>(pDstWord) = nVal;
2102 continue;
2103 case GDT_CInt16:
2104 {
2105 GInt16 *panDstWord = static_cast<GInt16 *>(pDstWord);
2106 if( nVal > 32767 )
2107 nVal = 32767;
2108 panDstWord[0] = nVal;
2109 panDstWord[1] = 0;
2110 continue;
2111 }
2112 case GDT_CInt32:
2113 {
2114 GInt32 *panDstWord = static_cast<GInt32 *>(pDstWord);
2115 panDstWord[0] = nVal;
2116 panDstWord[1] = 0;
2117 continue;
2118 }
2119 default:
2120 break;
2121 }
2122 dfPixelValue = nVal;
2123 }
2124 break;
2125
2126 case GDT_Int16:
2127 {
2128 GInt16 nVal = *static_cast<GInt16 *>(pSrcWord);
2129 switch( eDstType )
2130 {
2131 case GDT_Byte:
2132 {
2133 GByte byVal;
2134 if( nVal > 255 )
2135 byVal = 255;
2136 else if (nVal < 0)
2137 byVal = 0;
2138 else
2139 byVal = static_cast<GByte>(nVal);
2140 *static_cast<GByte *>(pDstWord) = byVal;
2141 continue;
2142 }
2143 case GDT_UInt16:
2144 if( nVal < 0 )
2145 nVal = 0;
2146 *static_cast<GUInt16 *>(pDstWord) = nVal;
2147 continue;
2148 case GDT_UInt32:
2149 if( nVal < 0 )
2150 nVal = 0;
2151 *static_cast<GUInt32 *>(pDstWord) = nVal;
2152 continue;
2153 case GDT_Int32:
2154 *static_cast<GInt32 *>(pDstWord) = nVal;
2155 continue;
2156 case GDT_CInt16:
2157 {
2158 GInt16 *panDstWord = static_cast<GInt16 *>(pDstWord);
2159 panDstWord[0] = nVal;
2160 panDstWord[1] = 0;
2161 continue;
2162 }
2163 case GDT_CInt32:
2164 {
2165 GInt32 *panDstWord = static_cast<GInt32 *>(pDstWord);
2166 panDstWord[0] = nVal;
2167 panDstWord[1] = 0;
2168 continue;
2169 }
2170 default:
2171 break;
2172 }
2173 dfPixelValue = nVal;
2174 }
2175 break;
2176
2177 case GDT_Int32:
2178 {
2179 GInt32 nVal = *static_cast<GInt32 *>(pSrcWord);
2180 switch( eDstType )
2181 {
2182 case GDT_Byte:
2183 {
2184 GByte byVal;
2185 if( nVal > 255 )
2186 byVal = 255;
2187 else if (nVal < 0)
2188 byVal = 0;
2189 else
2190 byVal = nVal;
2191 *static_cast<GByte *>(pDstWord) = byVal;
2192 continue;
2193 }
2194 case GDT_UInt16:
2195 if( nVal > 65535 )
2196 nVal = 65535;
2197 else if( nVal < 0 )
2198 nVal = 0;
2199 *static_cast<GUInt16 *>(pDstWord) = nVal;
2200 continue;
2201 case GDT_Int16:
2202 if( nVal > 32767 )
2203 nVal = 32767;
2204 else if( nVal < -32768)
2205 nVal = -32768;
2206 *static_cast<GInt16 *>(pDstWord) = nVal;
2207 continue;
2208 case GDT_UInt32:
2209 if( nVal < 0 )
2210 nVal = 0;
2211 *static_cast<GUInt32 *>(pDstWord) = nVal;
2212 continue;
2213 case GDT_CInt16:
2214 {
2215 GInt16 *panDstWord = static_cast<GInt16 *>(pDstWord);
2216 if( nVal > 32767 )
2217 nVal = 32767;
2218 else if( nVal < -32768)
2219 nVal = -32768;
2220 panDstWord[0] = nVal;
2221 panDstWord[1] = 0;
2222 continue;
2223 }
2224 case GDT_CInt32:
2225 {
2226 GInt32 *panDstWord = static_cast<GInt32 *>(pDstWord);
2227 panDstWord[0] = nVal;
2228 panDstWord[1] = 0;
2229 continue;
2230 }
2231 default:
2232 break;
2233 }
2234 dfPixelValue = nVal;
2235 }
2236 break;
2237
2238 case GDT_UInt32:
2239 {
2240 GUInt32 nVal = *static_cast<GUInt32 *>(pSrcWord);
2241 switch( eDstType )
2242 {
2243 case GDT_Byte:
2244 {
2245 GByte byVal;
2246 if( nVal > 255 )
2247 byVal = 255;
2248 else
2249 byVal = nVal;
2250 *static_cast<GByte *>(pDstWord) = byVal;
2251 continue;
2252 }
2253 case GDT_UInt16:
2254 if( nVal > 65535 )
2255 nVal = 65535;
2256 *static_cast<GUInt16 *>(pDstWord) = nVal;
2257 continue;
2258 case GDT_Int16:
2259 if( nVal > 32767 )
2260 nVal = 32767;
2261 *static_cast<GInt16 *>(pDstWord) = nVal;
2262 continue;
2263 case GDT_Int32:
2264 if( nVal > 2147483647UL )
2265 nVal = 2147483647UL;
2266 *static_cast<GInt32 *>(pDstWord) = nVal;
2267 continue;
2268 case GDT_CInt16:
2269 {
2270 GInt16 *panDstWord = static_cast<GInt16 *>(pDstWord);
2271 if( nVal > 32767 )
2272 nVal = 32767;
2273 panDstWord[0] = nVal;
2274 panDstWord[1] = 0;
2275 continue;
2276 }
2277 case GDT_CInt32:
2278 {
2279 GInt32 *panDstWord = static_cast<GInt32 *>(pDstWord);
2280 if( nVal > 2147483647UL )
2281 nVal = 2147483647UL;
2282 panDstWord[0] = nVal;
2283 panDstWord[1] = 0;
2284 continue;
2285 }
2286 default:
2287 break;
2288 }
2289 dfPixelValue = nVal;
2290 }
2291 break;
2292
2293 case GDT_CInt16:
2294 {
2295 GInt16 *panSrcWord = static_cast<GInt16 *>(pSrcWord);
2296 GInt16 nVal = panSrcWord[0];
2297 switch( eDstType )
2298 {
2299 case GDT_Byte:
2300 {
2301 GByte byVal;
2302 if( nVal > 255 )
2303 byVal = 255;
2304 else if (nVal < 0)
2305 byVal = 0;
2306 else
2307 byVal = static_cast<GByte>(nVal);
2308 *static_cast<GByte *>(pDstWord) = byVal;
2309 continue;
2310 }
2311 case GDT_Int16:
2312 *static_cast<GInt16 *>(pDstWord) = nVal;
2313 continue;
2314 case GDT_UInt16:
2315 if( nVal < 0 )
2316 nVal = 0;
2317 *static_cast<GUInt16 *>(pDstWord) = nVal;
2318 continue;
2319 case GDT_UInt32:
2320 if( nVal < 0 )
2321 nVal = 0;
2322 *static_cast<GUInt32 *>(pDstWord) = nVal;
2323 continue;
2324 case GDT_Int32:
2325 *static_cast<GInt32 *>(pDstWord) = nVal;
2326 continue;
2327 case GDT_CInt32:
2328 {
2329 GInt32 *panDstWord = static_cast<GInt32 *>(pDstWord);
2330 panDstWord[0] = panSrcWord[0];
2331 panDstWord[1] = panSrcWord[1];
2332 continue;
2333 }
2334 default:
2335 break;
2336 }
2337 dfPixelValue = panSrcWord[0];
2338 dfPixelValueI = panSrcWord[1];
2339 }
2340 break;
2341
2342 case GDT_CInt32:
2343 {
2344 GInt32 *panSrcWord = static_cast<GInt32 *>(pSrcWord);
2345 GInt32 nVal = panSrcWord[0];
2346 switch( eDstType )
2347 {
2348 case GDT_Byte:
2349 {
2350 GByte byVal;
2351 if( nVal > 255 )
2352 byVal = 255;
2353 else if (nVal < 0)
2354 byVal = 0;
2355 else
2356 byVal = nVal;
2357 *static_cast<GByte *>(pDstWord) = byVal;
2358 continue;
2359 }
2360 case GDT_Int16:
2361 if( nVal > 32767 )
2362 nVal = 32767;
2363 else if( nVal < -32768)
2364 nVal = -32768;
2365 *static_cast<GInt16 *>(pDstWord) = nVal;
2366 continue;
2367 case GDT_UInt16:
2368 if( nVal > 65535 )
2369 nVal = 65535;
2370 else if( nVal < 0 )
2371 nVal = 0;
2372 *static_cast<GUInt16 *>(pDstWord) = nVal;
2373 continue;
2374 case GDT_UInt32:
2375 if( nVal < 0 )
2376 nVal = 0;
2377 *static_cast<GUInt32 *>(pDstWord) = nVal;
2378 continue;
2379 case GDT_Int32:
2380 *static_cast<GInt32 *>(pDstWord) = nVal;
2381 continue;
2382 case GDT_CInt16:
2383 {
2384 GInt16 *panDstWord = static_cast<GInt16 *>(pDstWord);
2385 if( nVal > 32767 )
2386 nVal = 32767;
2387 else if( nVal < -32768)
2388 nVal = -32768;
2389 panDstWord[0] = nVal;
2390 nVal = panSrcWord[1];
2391 if( nVal > 32767 )
2392 nVal = 32767;
2393 else if( nVal < -32768)
2394 nVal = -32768;
2395 panDstWord[1] = nVal;
2396 continue;
2397 }
2398 default:
2399 break;
2400 }
2401 dfPixelValue = panSrcWord[0];
2402 dfPixelValueI = panSrcWord[1];
2403 }
2404 break;
2405
2406 case GDT_Float32:
2407 {
2408 float fVal = *static_cast<float *>(pSrcWord);
2409 dfPixelValue = fVal;
2410 }
2411 break;
2412
2413 case GDT_Float64:
2414 {
2415 dfPixelValue = *static_cast<double *>(pSrcWord);
2416 }
2417 break;
2418
2419 case GDT_CFloat32:
2420 {
2421 float *pafSrcWord = static_cast<float *>(pSrcWord);
2422 dfPixelValue = pafSrcWord[0];
2423 dfPixelValueI = pafSrcWord[1];
2424 }
2425 break;
2426
2427 case GDT_CFloat64:
2428 {
2429 double *padfSrcWord = static_cast<double *>(pSrcWord);
2430 dfPixelValue = padfSrcWord[0];
2431 dfPixelValueI = padfSrcWord[1];
2432 }
2433 break;
2434
2435 default:
2436 CPLAssert( FALSE );
2437 }
2438
2439 /* -------------------------------------------------------------------- */
2440 /* Set the destination pixel, doing range clipping as needed. */
2441 /* -------------------------------------------------------------------- */
2442 switch( eDstType )
2443 {
2444 case GDT_Byte:
2445 {
2446 GByte *pabyDstWord = static_cast<GByte *>(pDstWord);
2447
2448 dfPixelValue += (float) 0.5;
2449
2450 if( dfPixelValue < 0.0 )
2451 *pabyDstWord = 0;
2452 else if( dfPixelValue > 255.0 )
2453 *pabyDstWord = 255;
2454 else
2455 *pabyDstWord = (GByte) dfPixelValue;
2456 }
2457 break;
2458
2459 case GDT_UInt16:
2460 {
2461 GUInt16 nVal;
2462
2463 dfPixelValue += 0.5;
2464
2465 if( dfPixelValue < 0.0 )
2466 nVal = 0;
2467 else if( dfPixelValue > 65535.0 )
2468 nVal = 65535;
2469 else
2470 nVal = (GUInt16) dfPixelValue;
2471
2472 *static_cast<GUInt16 *>(pDstWord) = nVal;
2473 }
2474 break;
2475
2476 case GDT_Int16:
2477 {
2478 GInt16 nVal;
2479
2480 dfPixelValue += 0.5;
2481
2482 if( dfPixelValue < -32768 )
2483 nVal = -32768;
2484 else if( dfPixelValue > 32767 )
2485 nVal = 32767;
2486 else
2487 nVal = (GInt16) floor(dfPixelValue);
2488
2489 *static_cast<GInt16 *>(pDstWord) = nVal;
2490 }
2491 break;
2492
2493 case GDT_UInt32:
2494 {
2495 GUInt32 nVal;
2496
2497 dfPixelValue += 0.5;
2498
2499 if( dfPixelValue < 0 )
2500 nVal = 0;
2501 else if( dfPixelValue > 4294967295U )
2502 nVal = 4294967295U;
2503 else
2504 nVal = (GInt32) dfPixelValue;
2505
2506 *static_cast<GUInt32 *>(pDstWord) = nVal;
2507 }
2508 break;
2509
2510 case GDT_Int32:
2511 {
2512 GInt32 nVal;
2513
2514 dfPixelValue += 0.5;
2515
2516 if( dfPixelValue < -2147483648.0 )
2517 nVal = INT_MIN;
2518 else if( dfPixelValue > 2147483647 )
2519 nVal = 2147483647;
2520 else
2521 nVal = (GInt32) floor(dfPixelValue);
2522
2523 *static_cast<GInt32 *>(pDstWord) = nVal;
2524 }
2525 break;
2526
2527 case GDT_Float32:
2528 {
2529 *static_cast<float *>(pDstWord) = static_cast<float>(dfPixelValue);
2530 }
2531 break;
2532
2533 case GDT_Float64:
2534 *static_cast<double *>(pDstWord) = dfPixelValue;
2535 break;
2536
2537 case GDT_CInt16:
2538 {
2539 GInt16 nVal;
2540 GInt16 *panDstWord = static_cast<GInt16 *>(pDstWord);
2541
2542 dfPixelValue += 0.5;
2543 dfPixelValueI += 0.5;
2544
2545 if( dfPixelValue < -32768 )
2546 nVal = -32768;
2547 else if( dfPixelValue > 32767 )
2548 nVal = 32767;
2549 else
2550 nVal = (GInt16) floor(dfPixelValue);
2551 panDstWord[0] = nVal;
2552
2553 if( dfPixelValueI < -32768 )
2554 nVal = -32768;
2555 else if( dfPixelValueI > 32767 )
2556 nVal = 32767;
2557 else
2558 nVal = (GInt16) floor(dfPixelValueI);
2559 panDstWord[1] = nVal;
2560 }
2561 break;
2562
2563 case GDT_CInt32:
2564 {
2565 GInt32 nVal;
2566 GInt32 *panDstWord = static_cast<GInt32 *>(pDstWord);
2567
2568 dfPixelValue += 0.5;
2569 dfPixelValueI += 0.5;
2570
2571 if( dfPixelValue < -2147483648.0 )
2572 nVal = INT_MIN;
2573 else if( dfPixelValue > 2147483647 )
2574 nVal = 2147483647;
2575 else
2576 nVal = (GInt32) floor(dfPixelValue);
2577
2578 panDstWord[0] = nVal;
2579
2580 if( dfPixelValueI < -2147483648.0 )
2581 nVal = INT_MIN;
2582 else if( dfPixelValueI > 2147483647 )
2583 nVal = 2147483647;
2584 else
2585 nVal = (GInt32) floor(dfPixelValueI);
2586
2587 panDstWord[1] = nVal;
2588 }
2589 break;
2590
2591 case GDT_CFloat32:
2592 {
2593 float *pafDstWord = static_cast<float *>(pDstWord);
2594 pafDstWord[0] = static_cast<float>(dfPixelValue);
2595 pafDstWord[1] = static_cast<float>(dfPixelValueI);
2596 }
2597 break;
2598
2599 case GDT_CFloat64:
2600 {
2601 double *padfDstWord = static_cast<double *>(pDstWord);
2602 padfDstWord[0] = dfPixelValue;
2603 padfDstWord[1] = dfPixelValueI;
2604 }
2605 break;
2606
2607 default:
2608 CPLAssert( FALSE );
2609 }
2610 } /* next iWord */
2611 #endif // defined USE_NEW_COPYWORDS
2612 }
2613
2614 /************************************************************************/
2615 /* GDALCopyBits() */
2616 /************************************************************************/
2617
2618 /**
2619 * Bitwise word copying.
2620 *
2621 * A function for moving sets of partial bytes around. Loosely
2622 * speaking this is a bitswise analog to GDALCopyWords().
2623 *
2624 * It copies nStepCount "words" where each word is nBitCount bits long.
2625 * The nSrcStep and nDstStep are the number of bits from the start of one
2626 * word to the next (same as nBitCount if they are packed). The nSrcOffset
2627 * and nDstOffset are the offset into the source and destination buffers
2628 * to start at, also measured in bits.
2629 *
2630 * All bit offsets are assumed to start from the high order bit in a byte
2631 * (ie. most significant bit first). Currently this function is not very
2632 * optimized, but it may be improved for some common cases in the future
2633 * as needed.
2634 *
2635 * @param pabySrcData the source data buffer.
2636 * @param nSrcOffset the offset (in bits) in pabySrcData to the start of the
2637 * first word to copy.
2638 * @param nSrcStep the offset in bits from the start one source word to the
2639 * start of the next.
2640 * @param pabyDstData the destination data buffer.
2641 * @param nDstOffset the offset (in bits) in pabyDstData to the start of the
2642 * first word to copy over.
2643 * @param nDstStep the offset in bits from the start one word to the
2644 * start of the next.
2645 * @param nBitCount the number of bits in a word to be copied.
2646 * @param nStepCount the number of words to copy.
2647 */
2648
GDALCopyBits(const GByte * pabySrcData,int nSrcOffset,int nSrcStep,GByte * pabyDstData,int nDstOffset,int nDstStep,int nBitCount,int nStepCount)2649 void GDALCopyBits( const GByte *pabySrcData, int nSrcOffset, int nSrcStep,
2650 GByte *pabyDstData, int nDstOffset, int nDstStep,
2651 int nBitCount, int nStepCount )
2652
2653 {
2654 VALIDATE_POINTER0( pabySrcData, "GDALCopyBits" );
2655
2656 int iStep;
2657 int iBit;
2658
2659 for( iStep = 0; iStep < nStepCount; iStep++ )
2660 {
2661 for( iBit = 0; iBit < nBitCount; iBit++ )
2662 {
2663 if( pabySrcData[nSrcOffset>>3]
2664 & (0x80 >>(nSrcOffset & 7)) )
2665 pabyDstData[nDstOffset>>3] |= (0x80 >> (nDstOffset & 7));
2666 else
2667 pabyDstData[nDstOffset>>3] &= ~(0x80 >> (nDstOffset & 7));
2668
2669
2670 nSrcOffset++;
2671 nDstOffset++;
2672 }
2673
2674 nSrcOffset += (nSrcStep - nBitCount);
2675 nDstOffset += (nDstStep - nBitCount);
2676 }
2677 }
2678
2679 /************************************************************************/
2680 /* GDALGetBestOverviewLevel() */
2681 /* */
2682 /* Returns the best overview level to satisfy the query or -1 if none */
2683 /* Also updates nXOff, nYOff, nXSize, nYSize and psExtraArg when */
2684 /* returning a valid overview level */
2685 /************************************************************************/
2686
GDALBandGetBestOverviewLevel(GDALRasterBand * poBand,int & nXOff,int & nYOff,int & nXSize,int & nYSize,int nBufXSize,int nBufYSize)2687 int GDALBandGetBestOverviewLevel(GDALRasterBand* poBand,
2688 int &nXOff, int &nYOff,
2689 int &nXSize, int &nYSize,
2690 int nBufXSize, int nBufYSize )
2691 {
2692 return GDALBandGetBestOverviewLevel2(poBand, nXOff, nYOff, nXSize, nYSize,
2693 nBufXSize, nBufYSize, NULL);
2694 }
2695
GDALBandGetBestOverviewLevel2(GDALRasterBand * poBand,int & nXOff,int & nYOff,int & nXSize,int & nYSize,int nBufXSize,int nBufYSize,GDALRasterIOExtraArg * psExtraArg)2696 int GDALBandGetBestOverviewLevel2(GDALRasterBand* poBand,
2697 int &nXOff, int &nYOff,
2698 int &nXSize, int &nYSize,
2699 int nBufXSize, int nBufYSize,
2700 GDALRasterIOExtraArg* psExtraArg )
2701 {
2702 double dfDesiredResolution;
2703 /* -------------------------------------------------------------------- */
2704 /* Compute the desired resolution. The resolution is */
2705 /* based on the least reduced axis, and represents the number */
2706 /* of source pixels to one destination pixel. */
2707 /* -------------------------------------------------------------------- */
2708 if( (nXSize / (double) nBufXSize) < (nYSize / (double) nBufYSize )
2709 || nBufYSize == 1 )
2710 dfDesiredResolution = nXSize / (double) nBufXSize;
2711 else
2712 dfDesiredResolution = nYSize / (double) nBufYSize;
2713
2714 /* -------------------------------------------------------------------- */
2715 /* Find the overview level that largest resolution value (most */
2716 /* downsampled) that is still less than (or only a little more) */
2717 /* downsampled than the request. */
2718 /* -------------------------------------------------------------------- */
2719 int nOverviewCount = poBand->GetOverviewCount();
2720 GDALRasterBand* poBestOverview = NULL;
2721 double dfBestResolution = 0;
2722 int nBestOverviewLevel = -1;
2723
2724 for( int iOverview = 0; iOverview < nOverviewCount; iOverview++ )
2725 {
2726 GDALRasterBand *poOverview = poBand->GetOverview( iOverview );
2727 if (poOverview == NULL)
2728 continue;
2729
2730 double dfResolution;
2731
2732 // What resolution is this?
2733 if( (poBand->GetXSize() / (double) poOverview->GetXSize())
2734 < (poBand->GetYSize() / (double) poOverview->GetYSize()) )
2735 dfResolution =
2736 poBand->GetXSize() / (double) poOverview->GetXSize();
2737 else
2738 dfResolution =
2739 poBand->GetYSize() / (double) poOverview->GetYSize();
2740
2741 // Is it nearly the requested resolution and better (lower) than
2742 // the current best resolution?
2743 if( dfResolution >= dfDesiredResolution * 1.2
2744 || dfResolution <= dfBestResolution )
2745 continue;
2746
2747 // Ignore AVERAGE_BIT2GRAYSCALE overviews for RasterIO purposes.
2748 const char *pszResampling =
2749 poOverview->GetMetadataItem( "RESAMPLING" );
2750
2751 if( pszResampling != NULL && EQUALN(pszResampling,"AVERAGE_BIT2",12))
2752 continue;
2753
2754 // OK, this is our new best overview.
2755 poBestOverview = poOverview;
2756 nBestOverviewLevel = iOverview;
2757 dfBestResolution = dfResolution;
2758 }
2759
2760 /* -------------------------------------------------------------------- */
2761 /* If we didn't find an overview that helps us, just return */
2762 /* indicating failure and the full resolution image will be used. */
2763 /* -------------------------------------------------------------------- */
2764 if( nBestOverviewLevel < 0 )
2765 return -1;
2766
2767 /* -------------------------------------------------------------------- */
2768 /* Recompute the source window in terms of the selected */
2769 /* overview. */
2770 /* -------------------------------------------------------------------- */
2771 int nOXOff, nOYOff, nOXSize, nOYSize;
2772 double dfXRes, dfYRes;
2773
2774 dfXRes = poBand->GetXSize() / (double) poBestOverview->GetXSize();
2775 dfYRes = poBand->GetYSize() / (double) poBestOverview->GetYSize();
2776
2777 nOXOff = MIN(poBestOverview->GetXSize()-1,(int) (nXOff/dfXRes+0.5));
2778 nOYOff = MIN(poBestOverview->GetYSize()-1,(int) (nYOff/dfYRes+0.5));
2779 nOXSize = MAX(1,(int) (nXSize/dfXRes + 0.5));
2780 nOYSize = MAX(1,(int) (nYSize/dfYRes + 0.5));
2781 if( nOXOff + nOXSize > poBestOverview->GetXSize() )
2782 nOXSize = poBestOverview->GetXSize() - nOXOff;
2783 if( nOYOff + nOYSize > poBestOverview->GetYSize() )
2784 nOYSize = poBestOverview->GetYSize() - nOYOff;
2785
2786 nXOff = nOXOff;
2787 nYOff = nOYOff;
2788 nXSize = nOXSize;
2789 nYSize = nOYSize;
2790
2791 if( psExtraArg && psExtraArg->bFloatingPointWindowValidity )
2792 {
2793 psExtraArg->dfXOff /= dfXRes;
2794 psExtraArg->dfXSize /= dfXRes;
2795 psExtraArg->dfYOff /= dfYRes;
2796 psExtraArg->dfYSize /= dfYRes;
2797 }
2798
2799 return nBestOverviewLevel;
2800 }
2801
2802
2803 /************************************************************************/
2804 /* OverviewRasterIO() */
2805 /* */
2806 /* Special work function to utilize available overviews to */
2807 /* more efficiently satisfy downsampled requests. It will */
2808 /* return CE_Failure if there are no appropriate overviews */
2809 /* available but it doesn't emit any error messages. */
2810 /************************************************************************/
2811
OverviewRasterIO(GDALRWFlag eRWFlag,int nXOff,int nYOff,int nXSize,int nYSize,void * pData,int nBufXSize,int nBufYSize,GDALDataType eBufType,GSpacing nPixelSpace,GSpacing nLineSpace,GDALRasterIOExtraArg * psExtraArg)2812 CPLErr GDALRasterBand::OverviewRasterIO( GDALRWFlag eRWFlag,
2813 int nXOff, int nYOff, int nXSize, int nYSize,
2814 void * pData, int nBufXSize, int nBufYSize,
2815 GDALDataType eBufType,
2816 GSpacing nPixelSpace, GSpacing nLineSpace,
2817 GDALRasterIOExtraArg* psExtraArg )
2818
2819
2820 {
2821 int nOverview;
2822 GDALRasterIOExtraArg sExtraArg;
2823
2824 GDALCopyRasterIOExtraArg(&sExtraArg, psExtraArg);
2825
2826 nOverview =
2827 GDALBandGetBestOverviewLevel2(this, nXOff, nYOff, nXSize, nYSize,
2828 nBufXSize, nBufYSize, &sExtraArg);
2829 if (nOverview < 0)
2830 return CE_Failure;
2831
2832 /* -------------------------------------------------------------------- */
2833 /* Recast the call in terms of the new raster layer. */
2834 /* -------------------------------------------------------------------- */
2835 GDALRasterBand* poOverviewBand = GetOverview(nOverview);
2836 if (poOverviewBand == NULL)
2837 return CE_Failure;
2838
2839 return poOverviewBand->RasterIO( eRWFlag, nXOff, nYOff, nXSize, nYSize,
2840 pData, nBufXSize, nBufYSize, eBufType,
2841 nPixelSpace, nLineSpace, &sExtraArg );
2842 }
2843
2844 /************************************************************************/
2845 /* GetBestOverviewLevel() */
2846 /* */
2847 /* Returns the best overview level to satisfy the query or -1 if none */
2848 /* Also updates nXOff, nYOff, nXSize, nYSize when returning a valid */
2849 /* overview level */
2850 /************************************************************************/
2851
2852 static
GDALDatasetGetBestOverviewLevel(GDALDataset * poDS,int & nXOff,int & nYOff,int & nXSize,int & nYSize,int nBufXSize,int nBufYSize,int nBandCount,int * panBandMap)2853 int GDALDatasetGetBestOverviewLevel(GDALDataset* poDS,
2854 int &nXOff, int &nYOff,
2855 int &nXSize, int &nYSize,
2856 int nBufXSize, int nBufYSize,
2857 int nBandCount, int *panBandMap)
2858 {
2859 int iBand, iOverview;
2860 int nOverviewCount = 0;
2861 GDALRasterBand *poFirstBand = NULL;
2862
2863 if (nBandCount == 0)
2864 return -1;
2865
2866 /* -------------------------------------------------------------------- */
2867 /* Check that all bands have the same number of overviews and */
2868 /* that they have all the same size and block dimensions */
2869 /* -------------------------------------------------------------------- */
2870 for( iBand = 0; iBand < nBandCount; iBand++ )
2871 {
2872 GDALRasterBand *poBand = poDS->GetRasterBand( panBandMap[iBand] );
2873 if (iBand == 0)
2874 {
2875 poFirstBand = poBand;
2876 nOverviewCount = poBand->GetOverviewCount();
2877 }
2878 else if (nOverviewCount != poBand->GetOverviewCount())
2879 {
2880 CPLDebug( "GDAL",
2881 "GDALDataset::GetBestOverviewLevel() ... "
2882 "mismatched overview count, use std method." );
2883 return -1;
2884 }
2885 else
2886 {
2887 for(iOverview = 0; iOverview < nOverviewCount; iOverview++)
2888 {
2889 GDALRasterBand* poOvrBand =
2890 poBand->GetOverview(iOverview);
2891 GDALRasterBand* poOvrFirstBand =
2892 poFirstBand->GetOverview(iOverview);
2893 if ( poOvrBand == NULL || poOvrFirstBand == NULL)
2894 continue;
2895
2896 if ( poOvrFirstBand->GetXSize() != poOvrBand->GetXSize() ||
2897 poOvrFirstBand->GetYSize() != poOvrBand->GetYSize() )
2898 {
2899 CPLDebug( "GDAL",
2900 "GDALDataset::GetBestOverviewLevel() ... "
2901 "mismatched overview sizes, use std method." );
2902 return -1;
2903 }
2904 int nBlockXSizeFirst=0, nBlockYSizeFirst=0;
2905 int nBlockXSizeCurrent=0, nBlockYSizeCurrent=0;
2906 poOvrFirstBand->GetBlockSize(&nBlockXSizeFirst, &nBlockYSizeFirst);
2907 poOvrBand->GetBlockSize(&nBlockXSizeCurrent, &nBlockYSizeCurrent);
2908 if (nBlockXSizeFirst != nBlockXSizeCurrent ||
2909 nBlockYSizeFirst != nBlockYSizeCurrent)
2910 {
2911 CPLDebug( "GDAL",
2912 "GDALDataset::GetBestOverviewLevel() ... "
2913 "mismatched block sizes, use std method." );
2914 return -1;
2915 }
2916 }
2917 }
2918 }
2919
2920 return GDALBandGetBestOverviewLevel2(poFirstBand,
2921 nXOff, nYOff, nXSize, nYSize,
2922 nBufXSize, nBufYSize, NULL);
2923 }
2924
2925 /************************************************************************/
2926 /* BlockBasedRasterIO() */
2927 /* */
2928 /* This convenience function implements a dataset level */
2929 /* RasterIO() interface based on calling down to fetch blocks, */
2930 /* much like the GDALRasterBand::IRasterIO(), but it handles */
2931 /* all bands at once, so that a format driver that handles a */
2932 /* request for different bands of the same block efficiently */
2933 /* (ie. without re-reading interleaved data) will efficiently. */
2934 /* */
2935 /* This method is intended to be called by an overridden */
2936 /* IRasterIO() method in the driver specific GDALDataset */
2937 /* derived class. */
2938 /* */
2939 /* Default internal implementation of RasterIO() ... utilizes */
2940 /* the Block access methods to satisfy the request. This would */
2941 /* normally only be overridden by formats with overviews. */
2942 /* */
2943 /* To keep things relatively simple, this method does not */
2944 /* currently take advantage of some special cases addressed in */
2945 /* GDALRasterBand::IRasterIO(), so it is likely best to only */
2946 /* call it when you know it will help. That is in cases where */
2947 /* data is at 1:1 to the buffer, and you know the driver is */
2948 /* implementing interleaved IO efficiently on a block by block */
2949 /* basis. Overviews will be used when possible. */
2950 /************************************************************************/
2951
2952 CPLErr
BlockBasedRasterIO(GDALRWFlag eRWFlag,int nXOff,int nYOff,int nXSize,int nYSize,void * pData,int nBufXSize,int nBufYSize,GDALDataType eBufType,int nBandCount,int * panBandMap,GSpacing nPixelSpace,GSpacing nLineSpace,GSpacing nBandSpace,GDALRasterIOExtraArg * psExtraArg)2953 GDALDataset::BlockBasedRasterIO( GDALRWFlag eRWFlag,
2954 int nXOff, int nYOff, int nXSize, int nYSize,
2955 void * pData, int nBufXSize, int nBufYSize,
2956 GDALDataType eBufType,
2957 int nBandCount, int *panBandMap,
2958 GSpacing nPixelSpace, GSpacing nLineSpace,
2959 GSpacing nBandSpace,
2960 GDALRasterIOExtraArg* psExtraArg )
2961
2962 {
2963 GByte **papabySrcBlock = NULL;
2964 GDALRasterBlock *poBlock = NULL;
2965 GDALRasterBlock **papoBlocks = NULL;
2966 int nLBlockX=-1, nLBlockY=-1, iBufYOff, iBufXOff, iSrcY, iBand;
2967 int nBlockXSize=1, nBlockYSize=1;
2968 CPLErr eErr = CE_None;
2969 GDALDataType eDataType = GDT_Byte;
2970
2971 CPLAssert( NULL != pData );
2972
2973 /* -------------------------------------------------------------------- */
2974 /* Ensure that all bands share a common block size and data type. */
2975 /* -------------------------------------------------------------------- */
2976
2977 for( iBand = 0; iBand < nBandCount; iBand++ )
2978 {
2979 GDALRasterBand *poBand = GetRasterBand( panBandMap[iBand] );
2980 int nThisBlockXSize, nThisBlockYSize;
2981
2982 if( iBand == 0 )
2983 {
2984 poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
2985 eDataType = poBand->GetRasterDataType();
2986 }
2987 else
2988 {
2989 poBand->GetBlockSize( &nThisBlockXSize, &nThisBlockYSize );
2990 if( nThisBlockXSize != nBlockXSize
2991 || nThisBlockYSize != nBlockYSize )
2992 {
2993 CPLDebug( "GDAL",
2994 "GDALDataset::BlockBasedRasterIO() ... "
2995 "mismatched block sizes, use std method." );
2996 return GDALDataset::IRasterIO( eRWFlag,
2997 nXOff, nYOff, nXSize, nYSize,
2998 pData, nBufXSize, nBufYSize,
2999 eBufType,
3000 nBandCount, panBandMap,
3001 nPixelSpace, nLineSpace,
3002 nBandSpace, psExtraArg );
3003 }
3004
3005 if( eDataType != poBand->GetRasterDataType()
3006 && (nXSize != nBufXSize || nYSize != nBufYSize) )
3007 {
3008 CPLDebug( "GDAL",
3009 "GDALDataset::BlockBasedRasterIO() ... "
3010 "mismatched band data types, use std method." );
3011 return GDALDataset::IRasterIO( eRWFlag,
3012 nXOff, nYOff, nXSize, nYSize,
3013 pData, nBufXSize, nBufYSize,
3014 eBufType,
3015 nBandCount, panBandMap,
3016 nPixelSpace, nLineSpace,
3017 nBandSpace, psExtraArg );
3018 }
3019 }
3020 }
3021
3022 /* ==================================================================== */
3023 /* In this special case at full resolution we step through in */
3024 /* blocks, turning the request over to the per-band */
3025 /* IRasterIO(), but ensuring that all bands of one block are */
3026 /* called before proceeding to the next. */
3027 /* ==================================================================== */
3028
3029 if( nXSize == nBufXSize && nYSize == nBufYSize )
3030 {
3031 GDALRasterIOExtraArg sDummyExtraArg;
3032 INIT_RASTERIO_EXTRA_ARG(sDummyExtraArg);
3033
3034 int nChunkYSize, nChunkXSize, nChunkXOff, nChunkYOff;
3035
3036 for( iBufYOff = 0; iBufYOff < nBufYSize; iBufYOff += nChunkYSize )
3037 {
3038 nChunkYSize = nBlockYSize;
3039 nChunkYOff = iBufYOff + nYOff;
3040 nChunkYSize = nBlockYSize - (nChunkYOff % nBlockYSize);
3041 if( nChunkYSize == 0 )
3042 nChunkYSize = nBlockYSize;
3043 if( nChunkYOff + nChunkYSize > nYOff + nYSize )
3044 nChunkYSize = (nYOff + nYSize) - nChunkYOff;
3045
3046 for( iBufXOff = 0; iBufXOff < nBufXSize; iBufXOff += nChunkXSize )
3047 {
3048 nChunkXSize = nBlockXSize;
3049 nChunkXOff = iBufXOff + nXOff;
3050 nChunkXSize = nBlockXSize - (nChunkXOff % nBlockXSize);
3051 if( nChunkXSize == 0 )
3052 nChunkXSize = nBlockXSize;
3053 if( nChunkXOff + nChunkXSize > nXOff + nXSize )
3054 nChunkXSize = (nXOff + nXSize) - nChunkXOff;
3055
3056 GByte *pabyChunkData;
3057
3058 pabyChunkData = ((GByte *) pData)
3059 + iBufXOff * nPixelSpace
3060 + (GPtrDiff_t)iBufYOff * nLineSpace;
3061
3062 for( iBand = 0; iBand < nBandCount; iBand++ )
3063 {
3064 GDALRasterBand *poBand = GetRasterBand(panBandMap[iBand]);
3065
3066 eErr =
3067 poBand->GDALRasterBand::IRasterIO(
3068 eRWFlag, nChunkXOff, nChunkYOff,
3069 nChunkXSize, nChunkYSize,
3070 pabyChunkData + (GPtrDiff_t)iBand * nBandSpace,
3071 nChunkXSize, nChunkYSize, eBufType,
3072 nPixelSpace, nLineSpace, &sDummyExtraArg );
3073 if( eErr != CE_None )
3074 return eErr;
3075 }
3076 }
3077
3078 if( psExtraArg->pfnProgress != NULL &&
3079 !psExtraArg->pfnProgress(1.0 * MAX(nBufYSize,iBufYOff + nChunkYSize) / nBufYSize, "", psExtraArg->pProgressData) )
3080 {
3081 return CE_Failure;
3082 }
3083 }
3084
3085 return CE_None;
3086 }
3087
3088 /* Below code is not compatible with that case. It would need a complete */
3089 /* separate code like done in GDALRasterBand::IRasterIO. */
3090 if (eRWFlag == GF_Write && (nBufXSize < nXSize || nBufYSize < nYSize))
3091 {
3092 return GDALDataset::IRasterIO( eRWFlag,
3093 nXOff, nYOff, nXSize, nYSize,
3094 pData, nBufXSize, nBufYSize,
3095 eBufType,
3096 nBandCount, panBandMap,
3097 nPixelSpace, nLineSpace,
3098 nBandSpace, psExtraArg );
3099 }
3100
3101 /* We could have a smarter implementation, but that will do for now */
3102 if( psExtraArg->eResampleAlg != GRIORA_NearestNeighbour &&
3103 (nBufXSize != nXSize || nBufYSize != nYSize) )
3104 {
3105 return GDALDataset::IRasterIO( eRWFlag,
3106 nXOff, nYOff, nXSize, nYSize,
3107 pData, nBufXSize, nBufYSize,
3108 eBufType,
3109 nBandCount, panBandMap,
3110 nPixelSpace, nLineSpace,
3111 nBandSpace, psExtraArg );
3112 }
3113
3114 /* ==================================================================== */
3115 /* Loop reading required source blocks to satisfy output */
3116 /* request. This is the most general implementation. */
3117 /* ==================================================================== */
3118
3119 int nBandDataSize = GDALGetDataTypeSize( eDataType ) / 8;
3120
3121 papabySrcBlock = (GByte **) CPLCalloc(sizeof(GByte*),nBandCount);
3122 papoBlocks = (GDALRasterBlock **) CPLCalloc(sizeof(void*),nBandCount);
3123
3124 /* -------------------------------------------------------------------- */
3125 /* Select an overview level if appropriate. */
3126 /* -------------------------------------------------------------------- */
3127 int nOverviewLevel = GDALDatasetGetBestOverviewLevel (this,
3128 nXOff, nYOff, nXSize, nYSize,
3129 nBufXSize, nBufYSize,
3130 nBandCount, panBandMap);
3131 if (nOverviewLevel >= 0)
3132 {
3133 GetRasterBand(panBandMap[0])->GetOverview(nOverviewLevel)->
3134 GetBlockSize( &nBlockXSize, &nBlockYSize );
3135 }
3136
3137 /* -------------------------------------------------------------------- */
3138 /* Compute stepping increment. */
3139 /* -------------------------------------------------------------------- */
3140 double dfSrcX, dfSrcY, dfSrcXInc, dfSrcYInc;
3141
3142 dfSrcXInc = nXSize / (double) nBufXSize;
3143 dfSrcYInc = nYSize / (double) nBufYSize;
3144
3145 /* -------------------------------------------------------------------- */
3146 /* Loop over buffer computing source locations. */
3147 /* -------------------------------------------------------------------- */
3148 for( iBufYOff = 0; iBufYOff < nBufYSize; iBufYOff++ )
3149 {
3150 GPtrDiff_t iBufOffset;
3151 int iSrcOffset;
3152
3153 dfSrcY = (iBufYOff+0.5) * dfSrcYInc + nYOff;
3154 iSrcY = (int) dfSrcY;
3155
3156 iBufOffset = (GPtrDiff_t)iBufYOff * (GPtrDiff_t)nLineSpace;
3157
3158 for( iBufXOff = 0; iBufXOff < nBufXSize; iBufXOff++ )
3159 {
3160 int iSrcX;
3161
3162 dfSrcX = (iBufXOff+0.5) * dfSrcXInc + nXOff;
3163
3164 iSrcX = (int) dfSrcX;
3165
3166 // FIXME: this code likely doesn't work if the dirty block gets flushed
3167 // to disk before being completely written.
3168 // In the meantime, bJustInitalize should probably be set to FALSE
3169 // even if it is not ideal performance wise, and for lossy compression
3170
3171 /* -------------------------------------------------------------------- */
3172 /* Ensure we have the appropriate block loaded. */
3173 /* -------------------------------------------------------------------- */
3174 if( iSrcX < nLBlockX * nBlockXSize
3175 || iSrcX >= (nLBlockX+1) * nBlockXSize
3176 || iSrcY < nLBlockY * nBlockYSize
3177 || iSrcY >= (nLBlockY+1) * nBlockYSize )
3178 {
3179 nLBlockX = iSrcX / nBlockXSize;
3180 nLBlockY = iSrcY / nBlockYSize;
3181
3182 int bJustInitialize =
3183 eRWFlag == GF_Write
3184 && nYOff <= nLBlockY * nBlockYSize
3185 && nYOff + nYSize >= (nLBlockY+1) * nBlockYSize
3186 && nXOff <= nLBlockX * nBlockXSize
3187 && nXOff + nXSize >= (nLBlockX+1) * nBlockXSize;
3188 /*int bMemZeroBuffer = FALSE;
3189 if( eRWFlag == GF_Write && !bJustInitialize &&
3190 nXOff <= nLBlockX * nBlockXSize &&
3191 nYOff <= nLBlockY * nBlockYSize &&
3192 (nXOff + nXSize >= (nLBlockX+1) * nBlockXSize ||
3193 (nXOff + nXSize == GetRasterXSize() && (nLBlockX+1) * nBlockXSize > GetRasterXSize())) &&
3194 (nYOff + nYSize >= (nLBlockY+1) * nBlockYSize ||
3195 (nYOff + nYSize == GetRasterYSize() && (nLBlockY+1) * nBlockYSize > GetRasterYSize())) )
3196 {
3197 bJustInitialize = TRUE;
3198 bMemZeroBuffer = TRUE;
3199 }*/
3200 for( iBand = 0; iBand < nBandCount; iBand++ )
3201 {
3202 GDALRasterBand *poBand = GetRasterBand( panBandMap[iBand]);
3203 if (nOverviewLevel >= 0)
3204 poBand = poBand->GetOverview(nOverviewLevel);
3205 poBlock = poBand->GetLockedBlockRef( nLBlockX, nLBlockY,
3206 bJustInitialize );
3207 if( poBlock == NULL )
3208 {
3209 eErr = CE_Failure;
3210 goto CleanupAndReturn;
3211 }
3212
3213 if( eRWFlag == GF_Write )
3214 poBlock->MarkDirty();
3215
3216 if( papoBlocks[iBand] != NULL )
3217 papoBlocks[iBand]->DropLock();
3218
3219 papoBlocks[iBand] = poBlock;
3220
3221 papabySrcBlock[iBand] = (GByte *) poBlock->GetDataRef();
3222 if( papabySrcBlock[iBand] == NULL )
3223 {
3224 eErr = CE_Failure;
3225 goto CleanupAndReturn;
3226 }
3227 /*if( bMemZeroBuffer )
3228 {
3229 memset(papabySrcBlock[iBand], 0,
3230 nBandDataSize * nBlockXSize * nBlockYSize);
3231 }*/
3232 }
3233 }
3234
3235 /* -------------------------------------------------------------------- */
3236 /* Copy over this pixel of data. */
3237 /* -------------------------------------------------------------------- */
3238 iSrcOffset = ((GPtrDiff_t)iSrcX - (GPtrDiff_t)nLBlockX*nBlockXSize
3239 + ((GPtrDiff_t)iSrcY - (GPtrDiff_t)nLBlockY*nBlockYSize) * nBlockXSize)*nBandDataSize;
3240
3241 for( iBand = 0; iBand < nBandCount; iBand++ )
3242 {
3243 GByte *pabySrcBlock = papabySrcBlock[iBand];
3244 GPtrDiff_t iBandBufOffset = iBufOffset + (GPtrDiff_t)iBand * (GPtrDiff_t)nBandSpace;
3245
3246 if( eDataType == eBufType )
3247 {
3248 if( eRWFlag == GF_Read )
3249 memcpy( ((GByte *) pData) + iBandBufOffset,
3250 pabySrcBlock + iSrcOffset, nBandDataSize );
3251 else
3252 memcpy( pabySrcBlock + iSrcOffset,
3253 ((GByte *)pData) + iBandBufOffset, nBandDataSize );
3254 }
3255 else
3256 {
3257 /* type to type conversion ... ouch, this is expensive way
3258 of handling single words */
3259
3260 if( eRWFlag == GF_Read )
3261 GDALCopyWords( pabySrcBlock + iSrcOffset, eDataType, 0,
3262 ((GByte *) pData) + iBandBufOffset,
3263 eBufType, 0, 1 );
3264 else
3265 GDALCopyWords( ((GByte *) pData) + iBandBufOffset,
3266 eBufType, 0,
3267 pabySrcBlock + iSrcOffset, eDataType, 0,
3268 1 );
3269 }
3270 }
3271
3272 iBufOffset += (int)nPixelSpace;
3273 }
3274 }
3275
3276 /* -------------------------------------------------------------------- */
3277 /* CleanupAndReturn. */
3278 /* -------------------------------------------------------------------- */
3279 CleanupAndReturn:
3280 CPLFree( papabySrcBlock );
3281 if( papoBlocks != NULL )
3282 {
3283 for( iBand = 0; iBand < nBandCount; iBand++ )
3284 {
3285 if( papoBlocks[iBand] != NULL )
3286 papoBlocks[iBand]->DropLock();
3287 }
3288 CPLFree( papoBlocks );
3289 }
3290
3291 return( CE_None );
3292 }
3293
3294 /************************************************************************/
3295 /* GDALCopyWholeRasterGetSwathSize() */
3296 /************************************************************************/
3297
GDALCopyWholeRasterGetSwathSize(GDALRasterBand * poSrcPrototypeBand,GDALRasterBand * poDstPrototypeBand,int nBandCount,int bDstIsCompressed,int bInterleave,int * pnSwathCols,int * pnSwathLines)3298 static void GDALCopyWholeRasterGetSwathSize(GDALRasterBand *poSrcPrototypeBand,
3299 GDALRasterBand *poDstPrototypeBand,
3300 int nBandCount,
3301 int bDstIsCompressed, int bInterleave,
3302 int* pnSwathCols, int *pnSwathLines)
3303 {
3304 GDALDataType eDT = poDstPrototypeBand->GetRasterDataType();
3305 int nSrcBlockXSize, nSrcBlockYSize;
3306 int nBlockXSize, nBlockYSize;
3307
3308 int nXSize = poSrcPrototypeBand->GetXSize();
3309 int nYSize = poSrcPrototypeBand->GetYSize();
3310
3311 poSrcPrototypeBand->GetBlockSize( &nSrcBlockXSize, &nSrcBlockYSize );
3312 poDstPrototypeBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
3313
3314 int nMaxBlockXSize = MAX(nBlockXSize, nSrcBlockXSize);
3315 int nMaxBlockYSize = MAX(nBlockYSize, nSrcBlockYSize);
3316
3317 int nPixelSize = (GDALGetDataTypeSize(eDT) / 8);
3318 if( bInterleave)
3319 nPixelSize *= nBandCount;
3320
3321 // aim for one row of blocks. Do not settle for less.
3322 int nSwathCols = nXSize;
3323 int nSwathLines = nBlockYSize;
3324
3325 const char* pszSrcCompression =
3326 poSrcPrototypeBand->GetMetadataItem("COMPRESSION", "IMAGE_STRUCTURE");
3327
3328 /* -------------------------------------------------------------------- */
3329 /* What will our swath size be? */
3330 /* -------------------------------------------------------------------- */
3331 /* When writing interleaved data in a compressed format, we want to be sure */
3332 /* that each block will only be written once, so the swath size must not be */
3333 /* greater than the block cache. */
3334 const char* pszSwathSize = CPLGetConfigOption("GDAL_SWATH_SIZE", NULL);
3335 int nTargetSwathSize;
3336 if( pszSwathSize != NULL )
3337 nTargetSwathSize = atoi(pszSwathSize);
3338 else
3339 {
3340 /* As a default, take one 1/4 of the cache size */
3341 nTargetSwathSize = (int)MIN(INT_MAX, GDALGetCacheMax64() / 4);
3342
3343 /* but if the minimum idal swath buf size is less, then go for it to */
3344 /* avoid unnecessarily abusing RAM usage */
3345 GIntBig nIdealSwathBufSize = (GIntBig)nSwathCols * nSwathLines * nPixelSize;
3346 if( pszSrcCompression != NULL && EQUAL(pszSrcCompression, "JPEG2000") &&
3347 (!bDstIsCompressed || ((nSrcBlockXSize % nBlockXSize) == 0 && (nSrcBlockYSize % nBlockYSize) == 0)) )
3348 {
3349 nIdealSwathBufSize = MAX(nIdealSwathBufSize,
3350 (GIntBig)nSwathCols * nSrcBlockYSize * nPixelSize);
3351 }
3352 if( nTargetSwathSize > nIdealSwathBufSize )
3353 nTargetSwathSize = (int)nIdealSwathBufSize;
3354 }
3355
3356 if (nTargetSwathSize < 1000000)
3357 nTargetSwathSize = 1000000;
3358
3359 /* But let's check that */
3360 if (bDstIsCompressed && bInterleave && nTargetSwathSize > GDALGetCacheMax64())
3361 {
3362 CPLError(CE_Warning, CPLE_AppDefined,
3363 "When translating into a compressed interleave format, the block cache size (" CPL_FRMT_GIB ") "
3364 "should be at least the size of the swath (%d) (GDAL_SWATH_SIZE config. option)", GDALGetCacheMax64(), nTargetSwathSize);
3365 }
3366
3367 #define IS_DIVIDER_OF(x,y) ((y)%(x) == 0)
3368 #define ROUND_TO(x,y) (((x)/(y))*(y))
3369
3370 /* if both input and output datasets are tiled, that the tile dimensions */
3371 /* are "compatible", try to stick to a swath dimension that is a multiple */
3372 /* of input and output block dimensions */
3373 if (nBlockXSize != nXSize && nSrcBlockXSize != nXSize &&
3374 IS_DIVIDER_OF(nBlockXSize, nMaxBlockXSize) &&
3375 IS_DIVIDER_OF(nSrcBlockXSize, nMaxBlockXSize) &&
3376 IS_DIVIDER_OF(nBlockYSize, nMaxBlockYSize) &&
3377 IS_DIVIDER_OF(nSrcBlockYSize, nMaxBlockYSize))
3378 {
3379 if (((GIntBig)nMaxBlockXSize) * nMaxBlockYSize * nPixelSize <=
3380 (GIntBig)nTargetSwathSize)
3381 {
3382 nSwathCols = nTargetSwathSize / (nMaxBlockYSize * nPixelSize);
3383 nSwathCols = ROUND_TO(nSwathCols, nMaxBlockXSize);
3384 if (nSwathCols == 0)
3385 nSwathCols = nMaxBlockXSize;
3386 if (nSwathCols > nXSize)
3387 nSwathCols = nXSize;
3388 nSwathLines = nMaxBlockYSize;
3389
3390 if (((GIntBig)nSwathCols) * nSwathLines * nPixelSize >
3391 (GIntBig)nTargetSwathSize)
3392 {
3393 nSwathCols = nXSize;
3394 nSwathLines = nBlockYSize;
3395 }
3396 }
3397 }
3398
3399 int nMemoryPerCol = nSwathCols * nPixelSize;
3400
3401 /* Do the computation on a big int since for example when translating */
3402 /* the JPL WMS layer, we overflow 32 bits*/
3403 GIntBig nSwathBufSize = (GIntBig)nMemoryPerCol * nSwathLines;
3404 if (nSwathBufSize > (GIntBig)nTargetSwathSize)
3405 {
3406 nSwathLines = nTargetSwathSize / nMemoryPerCol;
3407 if (nSwathLines == 0)
3408 nSwathLines = 1;
3409
3410 CPLDebug( "GDAL",
3411 "GDALCopyWholeRasterGetSwathSize(): adjusting to %d line swath "
3412 "since requirement (" CPL_FRMT_GIB " bytes) exceed target swath size (%d bytes) (GDAL_SWATH_SIZE config. option)",
3413 nSwathLines, (GIntBig)nBlockYSize * nMemoryPerCol, nTargetSwathSize);
3414 }
3415 // If we are processing single scans, try to handle several at once.
3416 // If we are handling swaths already, only grow the swath if a row
3417 // of blocks is substantially less than our target buffer size.
3418 else if( nSwathLines == 1
3419 || nMemoryPerCol * nSwathLines < nTargetSwathSize / 10 )
3420 {
3421 nSwathLines = MIN(nYSize,MAX(1,nTargetSwathSize/nMemoryPerCol));
3422
3423 /* If possible try to align to source and target block height */
3424 if ((nSwathLines % nMaxBlockYSize) != 0 && nSwathLines > nMaxBlockYSize &&
3425 IS_DIVIDER_OF(nBlockYSize, nMaxBlockYSize) &&
3426 IS_DIVIDER_OF(nSrcBlockYSize, nMaxBlockYSize))
3427 nSwathLines = ROUND_TO(nSwathLines, nMaxBlockYSize);
3428 }
3429
3430 if( pszSrcCompression != NULL && EQUAL(pszSrcCompression, "JPEG2000") &&
3431 (!bDstIsCompressed ||
3432 (IS_DIVIDER_OF(nBlockXSize, nSrcBlockXSize) &&
3433 IS_DIVIDER_OF(nBlockYSize, nSrcBlockYSize))) )
3434 {
3435 /* Typical use case: converting from Pleaiades that is 2048x2048 tiled */
3436 if (nSwathLines < nSrcBlockYSize)
3437 {
3438 nSwathLines = nSrcBlockYSize;
3439
3440 /* Number of pixels that can be read/write simultaneously */
3441 nSwathCols = nTargetSwathSize / (nSrcBlockXSize * nPixelSize);
3442 nSwathCols = ROUND_TO(nSwathCols, nSrcBlockXSize);
3443 if (nSwathCols == 0)
3444 nSwathCols = nSrcBlockXSize;
3445 if (nSwathCols > nXSize)
3446 nSwathCols = nXSize;
3447
3448 CPLDebug( "GDAL",
3449 "GDALCopyWholeRasterGetSwathSize(): because of compression and too high block,\n"
3450 "use partial width at one time");
3451 }
3452 else if ((nSwathLines % nSrcBlockYSize) != 0)
3453 {
3454 /* Round on a multiple of nSrcBlockYSize */
3455 nSwathLines = ROUND_TO(nSwathLines, nSrcBlockYSize);
3456 CPLDebug( "GDAL",
3457 "GDALCopyWholeRasterGetSwathSize(): because of compression, \n"
3458 "round nSwathLines to block height : %d", nSwathLines);
3459 }
3460 }
3461 else if (bDstIsCompressed)
3462 {
3463 if (nSwathLines < nBlockYSize)
3464 {
3465 nSwathLines = nBlockYSize;
3466
3467 /* Number of pixels that can be read/write simultaneously */
3468 nSwathCols = nTargetSwathSize / (nSwathLines * nPixelSize);
3469 nSwathCols = ROUND_TO(nSwathCols, nBlockXSize);
3470 if (nSwathCols == 0)
3471 nSwathCols = nBlockXSize;
3472 if (nSwathCols > nXSize)
3473 nSwathCols = nXSize;
3474
3475 CPLDebug( "GDAL",
3476 "GDALCopyWholeRasterGetSwathSize(): because of compression and too high block,\n"
3477 "use partial width at one time");
3478 }
3479 else if ((nSwathLines % nBlockYSize) != 0)
3480 {
3481 /* Round on a multiple of nBlockYSize */
3482 nSwathLines = ROUND_TO(nSwathLines, nBlockYSize);
3483 CPLDebug( "GDAL",
3484 "GDALCopyWholeRasterGetSwathSize(): because of compression, \n"
3485 "round nSwathLines to block height : %d", nSwathLines);
3486 }
3487 }
3488
3489 *pnSwathCols = nSwathCols;
3490 *pnSwathLines = nSwathLines;
3491 }
3492
3493 /************************************************************************/
3494 /* GDALDatasetCopyWholeRaster() */
3495 /************************************************************************/
3496
3497 /**
3498 * \brief Copy all dataset raster data.
3499 *
3500 * This function copies the complete raster contents of one dataset to
3501 * another similarly configured dataset. The source and destination
3502 * dataset must have the same number of bands, and the same width
3503 * and height. The bands do not have to have the same data type.
3504 *
3505 * This function is primarily intended to support implementation of
3506 * driver specific CreateCopy() functions. It implements efficient copying,
3507 * in particular "chunking" the copy in substantial blocks and, if appropriate,
3508 * performing the transfer in a pixel interleaved fashion.
3509 *
3510 * Currently the only papszOptions value supported are : "INTERLEAVE=PIXEL"
3511 * to force pixel interleaved operation and "COMPRESSED=YES" to force alignment
3512 * on target dataset block sizes to achieve best compression. More options may be supported in
3513 * the future.
3514 *
3515 * @param hSrcDS the source dataset
3516 * @param hDstDS the destination dataset
3517 * @param papszOptions transfer hints in "StringList" Name=Value format.
3518 * @param pfnProgress progress reporting function.
3519 * @param pProgressData callback data for progress function.
3520 *
3521 * @return CE_None on success, or CE_Failure on failure.
3522 */
3523
GDALDatasetCopyWholeRaster(GDALDatasetH hSrcDS,GDALDatasetH hDstDS,char ** papszOptions,GDALProgressFunc pfnProgress,void * pProgressData)3524 CPLErr CPL_STDCALL GDALDatasetCopyWholeRaster(
3525 GDALDatasetH hSrcDS, GDALDatasetH hDstDS, char **papszOptions,
3526 GDALProgressFunc pfnProgress, void *pProgressData )
3527
3528 {
3529 VALIDATE_POINTER1( hSrcDS, "GDALDatasetCopyWholeRaster", CE_Failure );
3530 VALIDATE_POINTER1( hDstDS, "GDALDatasetCopyWholeRaster", CE_Failure );
3531
3532 GDALDataset *poSrcDS = (GDALDataset *) hSrcDS;
3533 GDALDataset *poDstDS = (GDALDataset *) hDstDS;
3534 CPLErr eErr = CE_None;
3535
3536 if( pfnProgress == NULL )
3537 pfnProgress = GDALDummyProgress;
3538
3539 /* -------------------------------------------------------------------- */
3540 /* Confirm the datasets match in size and band counts. */
3541 /* -------------------------------------------------------------------- */
3542 int nXSize = poDstDS->GetRasterXSize(),
3543 nYSize = poDstDS->GetRasterYSize(),
3544 nBandCount = poDstDS->GetRasterCount();
3545
3546 if( poSrcDS->GetRasterXSize() != nXSize
3547 || poSrcDS->GetRasterYSize() != nYSize
3548 || poSrcDS->GetRasterCount() != nBandCount )
3549 {
3550 CPLError( CE_Failure, CPLE_AppDefined,
3551 "Input and output dataset sizes or band counts do not\n"
3552 "match in GDALDatasetCopyWholeRaster()" );
3553 return CE_Failure;
3554 }
3555
3556 /* -------------------------------------------------------------------- */
3557 /* Report preliminary (0) progress. */
3558 /* -------------------------------------------------------------------- */
3559 if( !pfnProgress( 0.0, NULL, pProgressData ) )
3560 {
3561 CPLError( CE_Failure, CPLE_UserInterrupt,
3562 "User terminated CreateCopy()" );
3563 return CE_Failure;
3564 }
3565
3566 /* -------------------------------------------------------------------- */
3567 /* Get our prototype band, and assume the others are similarly */
3568 /* configured. */
3569 /* -------------------------------------------------------------------- */
3570 if( nBandCount == 0 )
3571 return CE_None;
3572
3573 GDALRasterBand *poSrcPrototypeBand = poSrcDS->GetRasterBand(1);
3574 GDALRasterBand *poDstPrototypeBand = poDstDS->GetRasterBand(1);
3575 GDALDataType eDT = poDstPrototypeBand->GetRasterDataType();
3576
3577 /* -------------------------------------------------------------------- */
3578 /* Do we want to try and do the operation in a pixel */
3579 /* interleaved fashion? */
3580 /* -------------------------------------------------------------------- */
3581 int bInterleave = FALSE;
3582 const char *pszInterleave = NULL;
3583
3584 pszInterleave = poSrcDS->GetMetadataItem( "INTERLEAVE", "IMAGE_STRUCTURE");
3585 if( pszInterleave != NULL
3586 && (EQUAL(pszInterleave,"PIXEL") || EQUAL(pszInterleave,"LINE")) )
3587 bInterleave = TRUE;
3588
3589 pszInterleave = poDstDS->GetMetadataItem( "INTERLEAVE", "IMAGE_STRUCTURE");
3590 if( pszInterleave != NULL
3591 && (EQUAL(pszInterleave,"PIXEL") || EQUAL(pszInterleave,"LINE")) )
3592 bInterleave = TRUE;
3593
3594 pszInterleave = CSLFetchNameValue( papszOptions, "INTERLEAVE" );
3595 if( pszInterleave != NULL
3596 && (EQUAL(pszInterleave,"PIXEL") || EQUAL(pszInterleave,"LINE")) )
3597 bInterleave = TRUE;
3598 else if( pszInterleave != NULL && EQUAL(pszInterleave,"BAND") )
3599 bInterleave = FALSE;
3600
3601 /* If the destination is compressed, we must try to write blocks just once, to save */
3602 /* disk space (GTiff case for example), and to avoid data loss (JPEG compression for example) */
3603 int bDstIsCompressed = FALSE;
3604 const char* pszDstCompressed= CSLFetchNameValue( papszOptions, "COMPRESSED" );
3605 if (pszDstCompressed != NULL && CSLTestBoolean(pszDstCompressed))
3606 bDstIsCompressed = TRUE;
3607
3608 /* -------------------------------------------------------------------- */
3609 /* What will our swath size be? */
3610 /* -------------------------------------------------------------------- */
3611
3612 int nSwathCols, nSwathLines;
3613 GDALCopyWholeRasterGetSwathSize(poSrcPrototypeBand,
3614 poDstPrototypeBand,
3615 nBandCount,
3616 bDstIsCompressed, bInterleave,
3617 &nSwathCols, &nSwathLines);
3618
3619 int nPixelSize = (GDALGetDataTypeSize(eDT) / 8);
3620 if( bInterleave)
3621 nPixelSize *= nBandCount;
3622
3623 void *pSwathBuf = VSIMalloc3(nSwathCols, nSwathLines, nPixelSize );
3624 if( pSwathBuf == NULL )
3625 {
3626 CPLError( CE_Failure, CPLE_OutOfMemory,
3627 "Failed to allocate %d*%d*%d byte swath buffer in\n"
3628 "GDALDatasetCopyWholeRaster()",
3629 nSwathCols, nSwathLines, nPixelSize );
3630 return CE_Failure;
3631 }
3632
3633 CPLDebug( "GDAL",
3634 "GDALDatasetCopyWholeRaster(): %d*%d swaths, bInterleave=%d",
3635 nSwathCols, nSwathLines, bInterleave );
3636
3637 if( nSwathCols == nXSize && poSrcDS->GetDriver() != NULL &&
3638 EQUAL(poSrcDS->GetDriver()->GetDescription(), "ECW") )
3639 {
3640 poSrcDS->AdviseRead(0, 0, nXSize, nYSize, nXSize, nYSize, eDT, nBandCount, NULL, NULL);
3641 }
3642
3643 /* ==================================================================== */
3644 /* Band oriented (uninterleaved) case. */
3645 /* ==================================================================== */
3646 if( !bInterleave )
3647 {
3648 int iBand, iX, iY;
3649
3650 GDALRasterIOExtraArg sExtraArg;
3651 INIT_RASTERIO_EXTRA_ARG(sExtraArg);
3652
3653 int nTotalBlocks = nBandCount *
3654 ((nYSize + nSwathLines - 1) / nSwathLines) *
3655 ((nXSize + nSwathCols - 1) / nSwathCols);
3656 int nBlocksDone = 0;
3657
3658 for( iBand = 0; iBand < nBandCount && eErr == CE_None; iBand++ )
3659 {
3660 int nBand = iBand+1;
3661
3662 for( iY = 0; iY < nYSize && eErr == CE_None; iY += nSwathLines )
3663 {
3664 int nThisLines = nSwathLines;
3665
3666 if( iY + nThisLines > nYSize )
3667 nThisLines = nYSize - iY;
3668
3669 for( iX = 0; iX < nXSize && eErr == CE_None; iX += nSwathCols )
3670 {
3671 int nThisCols = nSwathCols;
3672
3673 if( iX + nThisCols > nXSize )
3674 nThisCols = nXSize - iX;
3675
3676 sExtraArg.pfnProgress = GDALScaledProgress;
3677 sExtraArg.pProgressData =
3678 GDALCreateScaledProgress( nBlocksDone / (double)nTotalBlocks,
3679 (nBlocksDone + 0.5) / (double)nTotalBlocks,
3680 pfnProgress,
3681 pProgressData );
3682 if( sExtraArg.pProgressData == NULL )
3683 sExtraArg.pfnProgress = NULL;
3684
3685 eErr = poSrcDS->RasterIO( GF_Read,
3686 iX, iY, nThisCols, nThisLines,
3687 pSwathBuf, nThisCols, nThisLines,
3688 eDT, 1, &nBand,
3689 0, 0, 0, &sExtraArg );
3690
3691 GDALDestroyScaledProgress( sExtraArg.pProgressData );
3692
3693 if( eErr == CE_None )
3694 eErr = poDstDS->RasterIO( GF_Write,
3695 iX, iY, nThisCols, nThisLines,
3696 pSwathBuf, nThisCols, nThisLines,
3697 eDT, 1, &nBand,
3698 0, 0, 0, NULL );
3699 nBlocksDone ++;
3700 if( eErr == CE_None
3701 && !pfnProgress( nBlocksDone / (double)nTotalBlocks,
3702 NULL, pProgressData ) )
3703 {
3704 eErr = CE_Failure;
3705 CPLError( CE_Failure, CPLE_UserInterrupt,
3706 "User terminated CreateCopy()" );
3707 }
3708 }
3709 }
3710 }
3711 }
3712
3713 /* ==================================================================== */
3714 /* Pixel interleaved case. */
3715 /* ==================================================================== */
3716 else /* if( bInterleave ) */
3717 {
3718 int iY, iX;
3719
3720 GDALRasterIOExtraArg sExtraArg;
3721 INIT_RASTERIO_EXTRA_ARG(sExtraArg);
3722
3723 int nTotalBlocks = ((nYSize + nSwathLines - 1) / nSwathLines) *
3724 ((nXSize + nSwathCols - 1) / nSwathCols);
3725 int nBlocksDone = 0;
3726
3727 for( iY = 0; iY < nYSize && eErr == CE_None; iY += nSwathLines )
3728 {
3729 int nThisLines = nSwathLines;
3730
3731 if( iY + nThisLines > nYSize )
3732 nThisLines = nYSize - iY;
3733
3734 for( iX = 0; iX < nXSize && eErr == CE_None; iX += nSwathCols )
3735 {
3736 int nThisCols = nSwathCols;
3737
3738 if( iX + nThisCols > nXSize )
3739 nThisCols = nXSize - iX;
3740
3741 sExtraArg.pfnProgress = GDALScaledProgress;
3742 sExtraArg.pProgressData =
3743 GDALCreateScaledProgress( nBlocksDone / (double)nTotalBlocks,
3744 (nBlocksDone + 0.5) / (double)nTotalBlocks,
3745 pfnProgress,
3746 pProgressData );
3747 if( sExtraArg.pProgressData == NULL )
3748 sExtraArg.pfnProgress = NULL;
3749
3750 eErr = poSrcDS->RasterIO( GF_Read,
3751 iX, iY, nThisCols, nThisLines,
3752 pSwathBuf, nThisCols, nThisLines,
3753 eDT, nBandCount, NULL,
3754 0, 0, 0, &sExtraArg );
3755
3756 GDALDestroyScaledProgress( sExtraArg.pProgressData );
3757
3758 if( eErr == CE_None )
3759 eErr = poDstDS->RasterIO( GF_Write,
3760 iX, iY, nThisCols, nThisLines,
3761 pSwathBuf, nThisCols, nThisLines,
3762 eDT, nBandCount, NULL,
3763 0, 0, 0, NULL );
3764
3765 nBlocksDone ++;
3766 if( eErr == CE_None
3767 && !pfnProgress( nBlocksDone / (double)nTotalBlocks,
3768 NULL, pProgressData ) )
3769 {
3770 eErr = CE_Failure;
3771 CPLError( CE_Failure, CPLE_UserInterrupt,
3772 "User terminated CreateCopy()" );
3773 }
3774 }
3775 }
3776 }
3777
3778 /* -------------------------------------------------------------------- */
3779 /* Cleanup */
3780 /* -------------------------------------------------------------------- */
3781 CPLFree( pSwathBuf );
3782
3783 return eErr;
3784 }
3785
3786
3787 /************************************************************************/
3788 /* GDALRasterBandCopyWholeRaster() */
3789 /************************************************************************/
3790
3791 /**
3792 * \brief Copy all raster band raster data.
3793 *
3794 * This function copies the complete raster contents of one band to
3795 * another similarly configured band. The source and destination
3796 * bands must have the same width and height. The bands do not have
3797 * to have the same data type.
3798 *
3799 * It implements efficient copying, in particular "chunking" the copy in
3800 * substantial blocks.
3801 *
3802 * Currently the only papszOptions value supported is : "COMPRESSED=YES" to
3803 * force alignment on target dataset block sizes to achieve best compression.
3804 * More options may be supported in the future.
3805 *
3806 * @param hSrcBand the source band
3807 * @param hDstBand the destination band
3808 * @param papszOptions transfer hints in "StringList" Name=Value format.
3809 * @param pfnProgress progress reporting function.
3810 * @param pProgressData callback data for progress function.
3811 *
3812 * @return CE_None on success, or CE_Failure on failure.
3813 */
3814
GDALRasterBandCopyWholeRaster(GDALRasterBandH hSrcBand,GDALRasterBandH hDstBand,char ** papszOptions,GDALProgressFunc pfnProgress,void * pProgressData)3815 CPLErr CPL_STDCALL GDALRasterBandCopyWholeRaster(
3816 GDALRasterBandH hSrcBand, GDALRasterBandH hDstBand, char **papszOptions,
3817 GDALProgressFunc pfnProgress, void *pProgressData )
3818
3819 {
3820 VALIDATE_POINTER1( hSrcBand, "GDALRasterBandCopyWholeRaster", CE_Failure );
3821 VALIDATE_POINTER1( hDstBand, "GDALRasterBandCopyWholeRaster", CE_Failure );
3822
3823 GDALRasterBand *poSrcBand = (GDALRasterBand *) hSrcBand;
3824 GDALRasterBand *poDstBand = (GDALRasterBand *) hDstBand;
3825 CPLErr eErr = CE_None;
3826
3827 if( pfnProgress == NULL )
3828 pfnProgress = GDALDummyProgress;
3829
3830 /* -------------------------------------------------------------------- */
3831 /* Confirm the datasets match in size and band counts. */
3832 /* -------------------------------------------------------------------- */
3833 int nXSize = poSrcBand->GetXSize(),
3834 nYSize = poSrcBand->GetYSize();
3835
3836 if( poDstBand->GetXSize() != nXSize
3837 || poDstBand->GetYSize() != nYSize )
3838 {
3839 CPLError( CE_Failure, CPLE_AppDefined,
3840 "Input and output band sizes do not\n"
3841 "match in GDALRasterBandCopyWholeRaster()" );
3842 return CE_Failure;
3843 }
3844
3845 /* -------------------------------------------------------------------- */
3846 /* Report preliminary (0) progress. */
3847 /* -------------------------------------------------------------------- */
3848 if( !pfnProgress( 0.0, NULL, pProgressData ) )
3849 {
3850 CPLError( CE_Failure, CPLE_UserInterrupt,
3851 "User terminated CreateCopy()" );
3852 return CE_Failure;
3853 }
3854
3855 GDALDataType eDT = poDstBand->GetRasterDataType();
3856
3857 /* If the destination is compressed, we must try to write blocks just once, to save */
3858 /* disk space (GTiff case for example), and to avoid data loss (JPEG compression for example) */
3859 int bDstIsCompressed = FALSE;
3860 const char* pszDstCompressed= CSLFetchNameValue( papszOptions, "COMPRESSED" );
3861 if (pszDstCompressed != NULL && CSLTestBoolean(pszDstCompressed))
3862 bDstIsCompressed = TRUE;
3863
3864 /* -------------------------------------------------------------------- */
3865 /* What will our swath size be? */
3866 /* -------------------------------------------------------------------- */
3867
3868 int nSwathCols, nSwathLines;
3869 GDALCopyWholeRasterGetSwathSize(poSrcBand,
3870 poDstBand,
3871 1,
3872 bDstIsCompressed, FALSE,
3873 &nSwathCols, &nSwathLines);
3874
3875 int nPixelSize = (GDALGetDataTypeSize(eDT) / 8);
3876
3877 void *pSwathBuf = VSIMalloc3(nSwathCols, nSwathLines, nPixelSize );
3878 if( pSwathBuf == NULL )
3879 {
3880 CPLError( CE_Failure, CPLE_OutOfMemory,
3881 "Failed to allocate %d*%d*%d byte swath buffer in\n"
3882 "GDALRasterBandCopyWholeRaster()",
3883 nSwathCols, nSwathLines, nPixelSize );
3884 return CE_Failure;
3885 }
3886
3887 CPLDebug( "GDAL",
3888 "GDALRasterBandCopyWholeRaster(): %d*%d swaths",
3889 nSwathCols, nSwathLines );
3890
3891 /* ==================================================================== */
3892 /* Band oriented (uninterleaved) case. */
3893 /* ==================================================================== */
3894
3895 int iX, iY;
3896
3897 for( iY = 0; iY < nYSize && eErr == CE_None; iY += nSwathLines )
3898 {
3899 int nThisLines = nSwathLines;
3900
3901 if( iY + nThisLines > nYSize )
3902 nThisLines = nYSize - iY;
3903
3904 for( iX = 0; iX < nXSize && eErr == CE_None; iX += nSwathCols )
3905 {
3906 int nThisCols = nSwathCols;
3907
3908 if( iX + nThisCols > nXSize )
3909 nThisCols = nXSize - iX;
3910
3911 eErr = poSrcBand->RasterIO( GF_Read,
3912 iX, iY, nThisCols, nThisLines,
3913 pSwathBuf, nThisCols, nThisLines,
3914 eDT, 0, 0, NULL );
3915
3916 if( eErr == CE_None )
3917 eErr = poDstBand->RasterIO( GF_Write,
3918 iX, iY, nThisCols, nThisLines,
3919 pSwathBuf, nThisCols, nThisLines,
3920 eDT, 0, 0, NULL );
3921
3922 if( eErr == CE_None
3923 && !pfnProgress(
3924 (iY+nThisLines) / (float) (nYSize),
3925 NULL, pProgressData ) )
3926 {
3927 eErr = CE_Failure;
3928 CPLError( CE_Failure, CPLE_UserInterrupt,
3929 "User terminated CreateCopy()" );
3930 }
3931 }
3932 }
3933
3934 /* -------------------------------------------------------------------- */
3935 /* Cleanup */
3936 /* -------------------------------------------------------------------- */
3937 CPLFree( pSwathBuf );
3938
3939 return eErr;
3940 }
3941
3942
3943 /************************************************************************/
3944 /* GDALCopyRasterIOExtraArg () */
3945 /************************************************************************/
3946
GDALCopyRasterIOExtraArg(GDALRasterIOExtraArg * psDestArg,GDALRasterIOExtraArg * psSrcArg)3947 void GDALCopyRasterIOExtraArg(GDALRasterIOExtraArg* psDestArg,
3948 GDALRasterIOExtraArg* psSrcArg)
3949 {
3950 INIT_RASTERIO_EXTRA_ARG(*psDestArg);
3951 if( psSrcArg )
3952 {
3953 psDestArg->eResampleAlg = psSrcArg->eResampleAlg;
3954 psDestArg->pfnProgress = psSrcArg->pfnProgress;
3955 psDestArg->pProgressData = psSrcArg->pProgressData;
3956 psDestArg->bFloatingPointWindowValidity = psSrcArg->bFloatingPointWindowValidity;
3957 if( psSrcArg->bFloatingPointWindowValidity )
3958 {
3959 psDestArg->dfXOff = psSrcArg->dfXOff;
3960 psDestArg->dfYOff = psSrcArg->dfYOff;
3961 psDestArg->dfXSize = psSrcArg->dfXSize;
3962 psDestArg->dfYSize = psSrcArg->dfYSize;
3963 }
3964 }
3965 }
3966