1 : /******************************************************************************
2 : * $Id: overview.cpp 21356 2010-12-31 16:35:02Z rouault $
3 : *
4 : * Project: GDAL Core
5 : * Purpose: Helper code to implement overview support in different drivers.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2000, Frank Warmerdam
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "gdal_priv.h"
31 :
32 : CPL_CVSID("$Id: overview.cpp 21356 2010-12-31 16:35:02Z rouault $");
33 :
34 : typedef CPLErr (*GDALDownsampleFunction)
35 : ( int nSrcWidth, int nSrcHeight,
36 : GDALDataType eWrkDataType,
37 : void * pChunk,
38 : GByte * pabyChunkNodataMask,
39 : int nChunkXOff, int nChunkXSize,
40 : int nChunkYOff, int nChunkYSize,
41 : GDALRasterBand * poOverview,
42 : const char * pszResampling,
43 : int bHasNoData, float fNoDataValue,
44 : GDALColorTable* poColorTable,
45 : GDALDataType eSrcDataType);
46 :
47 : /************************************************************************/
48 : /* GDALDownsampleChunk32R_Near() */
49 : /************************************************************************/
50 :
51 : template <class T>
52 : static CPLErr
53 1413 : GDALDownsampleChunk32R_NearT( int nSrcWidth, int nSrcHeight,
54 : GDALDataType eWrkDataType,
55 : T * pChunk,
56 : GByte * pabyChunkNodataMask_unused,
57 : int nChunkXOff, int nChunkXSize,
58 : int nChunkYOff, int nChunkYSize,
59 : GDALRasterBand * poOverview,
60 : const char * pszResampling_unused,
61 : int bHasNoData_unused, float fNoDataValue_unused,
62 : GDALColorTable* poColorTable_unused,
63 : GDALDataType eSrcDataType)
64 :
65 : {
66 1413 : CPLErr eErr = CE_None;
67 :
68 : int nDstXOff, nDstXOff2, nDstYOff, nDstYOff2, nOXSize, nOYSize;
69 :
70 1413 : nOXSize = poOverview->GetXSize();
71 1413 : nOYSize = poOverview->GetYSize();
72 :
73 : /* -------------------------------------------------------------------- */
74 : /* Figure out the column to start writing to, and the first column */
75 : /* to not write to. */
76 : /* -------------------------------------------------------------------- */
77 1413 : nDstXOff = (int) (0.5 + (nChunkXOff/(double)nSrcWidth) * nOXSize);
78 1413 : nDstXOff2 = (int)
79 : (0.5 + ((nChunkXOff+nChunkXSize)/(double)nSrcWidth) * nOXSize);
80 :
81 1413 : if( nChunkXOff + nChunkXSize == nSrcWidth )
82 951 : nDstXOff2 = nOXSize;
83 :
84 1413 : int nDstXWidth = nDstXOff2 - nDstXOff;
85 :
86 : /* -------------------------------------------------------------------- */
87 : /* Allocate scanline buffer. */
88 : /* -------------------------------------------------------------------- */
89 :
90 1413 : T* pDstScanline = (T *) VSIMalloc(nDstXWidth * (GDALGetDataTypeSize(eWrkDataType) / 8));
91 1413 : int* panSrcXOff = (int*)VSIMalloc(nDstXWidth * sizeof(int));
92 :
93 1413 : if( pDstScanline == NULL || panSrcXOff == NULL )
94 : {
95 0 : CPLError( CE_Failure, CPLE_OutOfMemory,
96 : "GDALDownsampleChunk32R: Out of memory for line buffer." );
97 0 : VSIFree(pDstScanline);
98 0 : VSIFree(panSrcXOff);
99 0 : return CE_Failure;
100 : }
101 :
102 : /* -------------------------------------------------------------------- */
103 : /* Figure out the line to start writing to, and the first line */
104 : /* to not write to. In theory this approach should ensure that */
105 : /* every output line will be written if all input chunks are */
106 : /* processed. */
107 : /* -------------------------------------------------------------------- */
108 1413 : nDstYOff = (int) (0.5 + (nChunkYOff/(double)nSrcHeight) * nOYSize);
109 1413 : nDstYOff2 = (int)
110 : (0.5 + ((nChunkYOff+nChunkYSize)/(double)nSrcHeight) * nOYSize);
111 :
112 1413 : if( nChunkYOff + nChunkYSize == nSrcHeight )
113 413 : nDstYOff2 = nOYSize;
114 :
115 : /* ==================================================================== */
116 : /* Precompute inner loop constants. */
117 : /* ==================================================================== */
118 : int iDstPixel;
119 343301 : for( iDstPixel = nDstXOff; iDstPixel < nDstXOff2; iDstPixel++ )
120 : {
121 : int nSrcXOff;
122 :
123 341888 : nSrcXOff =
124 : (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
125 341888 : if ( nSrcXOff < nChunkXOff )
126 0 : nSrcXOff = nChunkXOff;
127 :
128 341888 : panSrcXOff[iDstPixel - nDstXOff] = nSrcXOff;
129 : }
130 :
131 : /* ==================================================================== */
132 : /* Loop over destination scanlines. */
133 : /* ==================================================================== */
134 103060 : for( int iDstLine = nDstYOff; iDstLine < nDstYOff2 && eErr == CE_None; iDstLine++ )
135 : {
136 : T *pSrcScanline;
137 : int nSrcYOff;
138 :
139 101647 : nSrcYOff = (int) (0.5 + (iDstLine/(double)nOYSize) * nSrcHeight);
140 101647 : if ( nSrcYOff < nChunkYOff )
141 140 : nSrcYOff = nChunkYOff;
142 :
143 101647 : pSrcScanline = pChunk + ((nSrcYOff-nChunkYOff) * nChunkXSize) - nChunkXOff;
144 :
145 : /* -------------------------------------------------------------------- */
146 : /* Loop over destination pixels */
147 : /* -------------------------------------------------------------------- */
148 23372312 : for( iDstPixel = 0; iDstPixel < nDstXWidth; iDstPixel++ )
149 : {
150 23270665 : pDstScanline[iDstPixel] = pSrcScanline[panSrcXOff[iDstPixel]];
151 : }
152 :
153 101647 : eErr = poOverview->RasterIO( GF_Write, nDstXOff, iDstLine, nDstXWidth, 1,
154 : pDstScanline, nDstXWidth, 1, eWrkDataType,
155 : 0, 0 );
156 : }
157 :
158 1413 : CPLFree( pDstScanline );
159 1413 : CPLFree( panSrcXOff );
160 :
161 1413 : return eErr;
162 : }
163 :
164 : static CPLErr
165 1413 : GDALDownsampleChunk32R_Near( int nSrcWidth, int nSrcHeight,
166 : GDALDataType eWrkDataType,
167 : void * pChunk,
168 : GByte * pabyChunkNodataMask_unused,
169 : int nChunkXOff, int nChunkXSize,
170 : int nChunkYOff, int nChunkYSize,
171 : GDALRasterBand * poOverview,
172 : const char * pszResampling_unused,
173 : int bHasNoData_unused, float fNoDataValue_unused,
174 : GDALColorTable* poColorTable_unused,
175 : GDALDataType eSrcDataType)
176 : {
177 1413 : if (eWrkDataType == GDT_Byte)
178 : return GDALDownsampleChunk32R_NearT(nSrcWidth, nSrcHeight,
179 : eWrkDataType,
180 : (GByte *) pChunk,
181 : pabyChunkNodataMask_unused,
182 : nChunkXOff, nChunkXSize,
183 : nChunkYOff, nChunkYSize,
184 : poOverview,
185 : pszResampling_unused,
186 : bHasNoData_unused, fNoDataValue_unused,
187 : poColorTable_unused,
188 1361 : eSrcDataType);
189 52 : else if (eWrkDataType == GDT_Float32)
190 : return GDALDownsampleChunk32R_NearT(nSrcWidth, nSrcHeight,
191 : eWrkDataType,
192 : (float *) pChunk,
193 : pabyChunkNodataMask_unused,
194 : nChunkXOff, nChunkXSize,
195 : nChunkYOff, nChunkYSize,
196 : poOverview,
197 : pszResampling_unused,
198 : bHasNoData_unused, fNoDataValue_unused,
199 : poColorTable_unused,
200 52 : eSrcDataType);
201 :
202 0 : CPLAssert(0);
203 0 : return CE_Failure;
204 : }
205 :
206 : /************************************************************************/
207 : /* GDALDownsampleChunk32R_Average() */
208 : /************************************************************************/
209 :
210 : template <class T, class Tsum>
211 : static CPLErr
212 149 : GDALDownsampleChunk32R_AverageT( int nSrcWidth, int nSrcHeight,
213 : GDALDataType eWrkDataType,
214 : T* pChunk,
215 : GByte * pabyChunkNodataMask,
216 : int nChunkXOff, int nChunkXSize,
217 : int nChunkYOff, int nChunkYSize,
218 : GDALRasterBand * poOverview,
219 : const char * pszResampling,
220 : int bHasNoData, float fNoDataValue,
221 : GDALColorTable* poColorTable,
222 : GDALDataType eSrcDataType)
223 :
224 : {
225 149 : CPLErr eErr = CE_None;
226 :
227 149 : int bBit2Grayscale = EQUALN(pszResampling,"AVERAGE_BIT2GRAYSCALE",13);
228 149 : if (bBit2Grayscale)
229 13 : poColorTable = NULL;
230 :
231 : int nDstXOff, nDstXOff2, nDstYOff, nDstYOff2, nOXSize, nOYSize;
232 : T *pDstScanline;
233 :
234 149 : T tNoDataValue = (T)fNoDataValue;
235 149 : if (!bHasNoData)
236 132 : tNoDataValue = 0;
237 :
238 149 : nOXSize = poOverview->GetXSize();
239 149 : nOYSize = poOverview->GetYSize();
240 :
241 : /* -------------------------------------------------------------------- */
242 : /* Figure out the column to start writing to, and the first column */
243 : /* to not write to. */
244 : /* -------------------------------------------------------------------- */
245 149 : nDstXOff = (int) (0.5 + (nChunkXOff/(double)nSrcWidth) * nOXSize);
246 149 : nDstXOff2 = (int)
247 : (0.5 + ((nChunkXOff+nChunkXSize)/(double)nSrcWidth) * nOXSize);
248 :
249 149 : if( nChunkXOff + nChunkXSize == nSrcWidth )
250 149 : nDstXOff2 = nOXSize;
251 :
252 149 : int nChunkRightXOff = MIN(nSrcWidth, nChunkXOff + nChunkXSize);
253 149 : int nDstXWidth = nDstXOff2 - nDstXOff;
254 :
255 : /* -------------------------------------------------------------------- */
256 : /* Allocate scanline buffer. */
257 : /* -------------------------------------------------------------------- */
258 :
259 149 : pDstScanline = (T *) VSIMalloc(nDstXWidth * (GDALGetDataTypeSize(eWrkDataType) / 8));
260 149 : int* panSrcXOffShifted = (int*)VSIMalloc(2 * nDstXWidth * sizeof(int));
261 :
262 149 : if( pDstScanline == NULL || panSrcXOffShifted == NULL )
263 : {
264 0 : CPLError( CE_Failure, CPLE_OutOfMemory,
265 : "GDALDownsampleChunk32R: Out of memory for line buffer." );
266 0 : VSIFree(pDstScanline);
267 0 : VSIFree(panSrcXOffShifted);
268 0 : return CE_Failure;
269 : }
270 :
271 : /* -------------------------------------------------------------------- */
272 : /* Figure out the line to start writing to, and the first line */
273 : /* to not write to. In theory this approach should ensure that */
274 : /* every output line will be written if all input chunks are */
275 : /* processed. */
276 : /* -------------------------------------------------------------------- */
277 149 : nDstYOff = (int) (0.5 + (nChunkYOff/(double)nSrcHeight) * nOYSize);
278 149 : nDstYOff2 = (int)
279 : (0.5 + ((nChunkYOff+nChunkYSize)/(double)nSrcHeight) * nOYSize);
280 :
281 149 : if( nChunkYOff + nChunkYSize == nSrcHeight )
282 89 : nDstYOff2 = nOYSize;
283 :
284 :
285 149 : int nEntryCount = 0;
286 149 : GDALColorEntry* aEntries = NULL;
287 149 : if (poColorTable)
288 : {
289 : int i;
290 2 : nEntryCount = poColorTable->GetColorEntryCount();
291 2 : aEntries = (GDALColorEntry* )CPLMalloc(sizeof(GDALColorEntry) * nEntryCount);
292 514 : for(i=0;i<nEntryCount;i++)
293 : {
294 512 : poColorTable->GetColorEntryAsRGB(i, &aEntries[i]);
295 : }
296 : }
297 :
298 : /* ==================================================================== */
299 : /* Precompute inner loop constants. */
300 : /* ==================================================================== */
301 : int iDstPixel;
302 149 : int bSrcXSpacingIsTwo = TRUE;
303 5159 : for( iDstPixel = nDstXOff; iDstPixel < nDstXOff2; iDstPixel++ )
304 : {
305 : int nSrcXOff, nSrcXOff2;
306 :
307 5010 : nSrcXOff =
308 : (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
309 5010 : if ( nSrcXOff < nChunkXOff )
310 0 : nSrcXOff = nChunkXOff;
311 5010 : nSrcXOff2 = (int)
312 : (0.5 + ((iDstPixel+1)/(double)nOXSize) * nSrcWidth);
313 :
314 5010 : if( nSrcXOff2 > nChunkRightXOff || iDstPixel == nOXSize-1 )
315 149 : nSrcXOff2 = nChunkRightXOff;
316 :
317 5010 : panSrcXOffShifted[2 * (iDstPixel - nDstXOff)] = nSrcXOff - nChunkXOff;
318 5010 : panSrcXOffShifted[2 * (iDstPixel - nDstXOff) + 1] = nSrcXOff2 - nChunkXOff;
319 5010 : if (nSrcXOff2 - nSrcXOff != 2)
320 783 : bSrcXSpacingIsTwo = FALSE;
321 : }
322 :
323 : /* ==================================================================== */
324 : /* Loop over destination scanlines. */
325 : /* ==================================================================== */
326 2359 : for( int iDstLine = nDstYOff; iDstLine < nDstYOff2 && eErr == CE_None; iDstLine++ )
327 : {
328 2210 : int nSrcYOff, nSrcYOff2 = 0;
329 :
330 2210 : nSrcYOff = (int) (0.5 + (iDstLine/(double)nOYSize) * nSrcHeight);
331 2210 : if ( nSrcYOff < nChunkYOff )
332 0 : nSrcYOff = nChunkYOff;
333 :
334 2210 : nSrcYOff2 =
335 : (int) (0.5 + ((iDstLine+1)/(double)nOYSize) * nSrcHeight);
336 :
337 2210 : if( nSrcYOff2 > nSrcHeight || iDstLine == nOYSize-1 )
338 89 : nSrcYOff2 = nSrcHeight;
339 2210 : if( nSrcYOff2 > nChunkYOff + nChunkYSize)
340 0 : nSrcYOff2 = nChunkYOff + nChunkYSize;
341 :
342 : /* -------------------------------------------------------------------- */
343 : /* Loop over destination pixels */
344 : /* -------------------------------------------------------------------- */
345 2210 : if (poColorTable == NULL)
346 : {
347 3105 : if (bSrcXSpacingIsTwo && nSrcYOff2 == nSrcYOff + 2 &&
348 : pabyChunkNodataMask == NULL && eWrkDataType == GDT_Byte)
349 : {
350 : /* Optimized case : no nodata, overview by a factor of 2 and regular x and y src spacing */
351 915 : T* pSrcScanlineShifted = pChunk + panSrcXOffShifted[0] + (nSrcYOff - nChunkYOff) * nChunkXSize;
352 43140 : for( iDstPixel = 0; iDstPixel < nDstXWidth; iDstPixel++ )
353 : {
354 : Tsum nTotal;
355 :
356 42225 : nTotal = pSrcScanlineShifted[0];
357 42225 : nTotal += pSrcScanlineShifted[1];
358 42225 : nTotal += pSrcScanlineShifted[nChunkXSize];
359 42225 : nTotal += pSrcScanlineShifted[1+nChunkXSize];
360 :
361 42225 : pDstScanline[iDstPixel] = (T) ((nTotal + 2) / 4);
362 42225 : pSrcScanlineShifted += 2;
363 : }
364 : }
365 : else
366 : {
367 1275 : nSrcYOff -= nChunkYOff;
368 1275 : nSrcYOff2 -= nChunkYOff;
369 :
370 48650 : for( iDstPixel = 0; iDstPixel < nDstXWidth; iDstPixel++ )
371 : {
372 47375 : int nSrcXOff = panSrcXOffShifted[2 * iDstPixel],
373 47375 : nSrcXOff2 = panSrcXOffShifted[2 * iDstPixel + 1];
374 :
375 : T val;
376 47375 : Tsum dfTotal = 0;
377 47375 : int nCount = 0, iX, iY;
378 :
379 149475 : for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
380 : {
381 336508 : for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
382 : {
383 234408 : val = pChunk[iX + iY *nChunkXSize];
384 234408 : if (pabyChunkNodataMask == NULL ||
385 : pabyChunkNodataMask[iX + iY *nChunkXSize])
386 : {
387 214308 : dfTotal += val;
388 214308 : nCount++;
389 : }
390 : }
391 : }
392 :
393 47375 : if( nCount == 0 )
394 4925 : pDstScanline[iDstPixel] = tNoDataValue;
395 42450 : else if (eWrkDataType == GDT_Byte)
396 42450 : pDstScanline[iDstPixel] = (T) ((dfTotal + nCount / 2) / nCount);
397 : else
398 0 : pDstScanline[iDstPixel] = (T) (dfTotal / nCount);
399 : }
400 : }
401 : }
402 : else
403 : {
404 20 : nSrcYOff -= nChunkYOff;
405 20 : nSrcYOff2 -= nChunkYOff;
406 :
407 220 : for( iDstPixel = 0; iDstPixel < nDstXWidth; iDstPixel++ )
408 : {
409 200 : int nSrcXOff = panSrcXOffShifted[2 * iDstPixel],
410 200 : nSrcXOff2 = panSrcXOffShifted[2 * iDstPixel + 1];
411 :
412 : T val;
413 200 : int nTotalR = 0, nTotalG = 0, nTotalB = 0;
414 200 : int nCount = 0, iX, iY;
415 :
416 600 : for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
417 : {
418 1200 : for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
419 : {
420 800 : val = pChunk[iX + iY *nChunkXSize];
421 800 : if (bHasNoData == FALSE || val != tNoDataValue)
422 : {
423 800 : int nVal = (int)val;
424 800 : if (nVal >= 0 && nVal < nEntryCount)
425 : {
426 800 : nTotalR += aEntries[nVal].c1;
427 800 : nTotalG += aEntries[nVal].c2;
428 800 : nTotalB += aEntries[nVal].c3;
429 800 : nCount++;
430 : }
431 : }
432 : }
433 : }
434 :
435 200 : if( nCount == 0 )
436 0 : pDstScanline[iDstPixel] = tNoDataValue;
437 : else
438 : {
439 200 : int nR = nTotalR / nCount, nG = nTotalG / nCount, nB = nTotalB / nCount;
440 : int i;
441 200 : Tsum dfMinDist = 0;
442 200 : int iBestEntry = 0;
443 51400 : for(i=0;i<nEntryCount;i++)
444 : {
445 : Tsum dfDist = (nR - aEntries[i].c1) * (nR - aEntries[i].c1) +
446 : (nG - aEntries[i].c2) * (nG - aEntries[i].c2) +
447 51200 : (nB - aEntries[i].c3) * (nB - aEntries[i].c3);
448 51200 : if (i == 0 || dfDist < dfMinDist)
449 : {
450 400 : dfMinDist = dfDist;
451 400 : iBestEntry = i;
452 : }
453 : }
454 200 : pDstScanline[iDstPixel] = (T)iBestEntry;
455 : }
456 : }
457 : }
458 :
459 2210 : eErr = poOverview->RasterIO( GF_Write, nDstXOff, iDstLine, nDstXWidth, 1,
460 : pDstScanline, nDstXWidth, 1, eWrkDataType,
461 : 0, 0 );
462 : }
463 :
464 149 : CPLFree( pDstScanline );
465 149 : CPLFree( aEntries );
466 149 : CPLFree( panSrcXOffShifted );
467 :
468 149 : return eErr;
469 : }
470 :
471 : static CPLErr
472 149 : GDALDownsampleChunk32R_Average( int nSrcWidth, int nSrcHeight,
473 : GDALDataType eWrkDataType,
474 : void * pChunk,
475 : GByte * pabyChunkNodataMask,
476 : int nChunkXOff, int nChunkXSize,
477 : int nChunkYOff, int nChunkYSize,
478 : GDALRasterBand * poOverview,
479 : const char * pszResampling,
480 : int bHasNoData, float fNoDataValue,
481 : GDALColorTable* poColorTable,
482 : GDALDataType eSrcDataType)
483 : {
484 149 : if (eWrkDataType == GDT_Byte)
485 : return GDALDownsampleChunk32R_AverageT<GByte, int>(nSrcWidth, nSrcHeight,
486 : eWrkDataType,
487 : (GByte *) pChunk,
488 : pabyChunkNodataMask,
489 : nChunkXOff, nChunkXSize,
490 : nChunkYOff, nChunkYSize,
491 : poOverview,
492 : pszResampling,
493 : bHasNoData, fNoDataValue,
494 : poColorTable,
495 149 : eSrcDataType);
496 0 : else if (eWrkDataType == GDT_Float32)
497 : return GDALDownsampleChunk32R_AverageT<float, double>(nSrcWidth, nSrcHeight,
498 : eWrkDataType,
499 : (float *) pChunk,
500 : pabyChunkNodataMask,
501 : nChunkXOff, nChunkXSize,
502 : nChunkYOff, nChunkYSize,
503 : poOverview,
504 : pszResampling,
505 : bHasNoData, fNoDataValue,
506 : poColorTable,
507 0 : eSrcDataType);
508 :
509 0 : CPLAssert(0);
510 0 : return CE_Failure;
511 : }
512 :
513 : /************************************************************************/
514 : /* GDALDownsampleChunk32R_Gauss() */
515 : /************************************************************************/
516 :
517 : static CPLErr
518 20 : GDALDownsampleChunk32R_Gauss( int nSrcWidth, int nSrcHeight,
519 : GDALDataType eWrkDataType,
520 : void * pChunk,
521 : GByte * pabyChunkNodataMask,
522 : int nChunkXOff, int nChunkXSize,
523 : int nChunkYOff, int nChunkYSize,
524 : GDALRasterBand * poOverview,
525 : const char * pszResampling,
526 : int bHasNoData, float fNoDataValue,
527 : GDALColorTable* poColorTable,
528 : GDALDataType eSrcDataType)
529 :
530 : {
531 20 : CPLErr eErr = CE_None;
532 :
533 20 : float * pafChunk = (float*) pChunk;
534 :
535 : /* -------------------------------------------------------------------- */
536 : /* Create the filter kernel and allocate scanline buffer. */
537 : /* -------------------------------------------------------------------- */
538 : int nDstXOff, nDstXOff2, nDstYOff, nDstYOff2, nOXSize, nOYSize;
539 : float *pafDstScanline;
540 20 : int nGaussMatrixDim = 3;
541 : const int *panGaussMatrix;
542 : static const int anGaussMatrix3x3[] ={
543 : 1,2,1,
544 : 2,4,2,
545 : 1,2,1
546 : };
547 : static const int anGaussMatrix5x5[] = {
548 : 1,4,6,4,1,
549 : 4,16,24,16,4,
550 : 6,24,36,24,6,
551 : 4,16,24,16,4,
552 : 1,4,6,4,1};
553 : static const int anGaussMatrix7x7[] = {
554 : 1,6,15,20,15,6,1,
555 : 6,36,90,120,90,36,6,
556 : 15,90,225,300,225,90,15,
557 : 20,120,300,400,300,120,20,
558 : 15,90,225,300,225,90,15,
559 : 6,36,90,120,90,36,6,
560 : 1,6,15,20,15,6,1};
561 :
562 20 : nOXSize = poOverview->GetXSize();
563 20 : nOYSize = poOverview->GetYSize();
564 20 : int nResYFactor = (int) (0.5 + (double)nSrcHeight/(double)nOYSize);
565 :
566 : // matrix for gauss filter
567 20 : if(nResYFactor <= 2 )
568 : {
569 20 : panGaussMatrix = anGaussMatrix3x3;
570 20 : nGaussMatrixDim=3;
571 : }
572 0 : else if (nResYFactor <= 4)
573 : {
574 0 : panGaussMatrix = anGaussMatrix5x5;
575 0 : nGaussMatrixDim=5;
576 : }
577 : else
578 : {
579 0 : panGaussMatrix = anGaussMatrix7x7;
580 0 : nGaussMatrixDim=7;
581 : }
582 :
583 : /* -------------------------------------------------------------------- */
584 : /* Figure out the column to start writing to, and the first column */
585 : /* to not write to. */
586 : /* -------------------------------------------------------------------- */
587 20 : nDstXOff = (int) (0.5 + (nChunkXOff/(double)nSrcWidth) * nOXSize);
588 : nDstXOff2 = (int)
589 20 : (0.5 + ((nChunkXOff+nChunkXSize)/(double)nSrcWidth) * nOXSize);
590 :
591 20 : if( nChunkXOff + nChunkXSize == nSrcWidth )
592 20 : nDstXOff2 = nOXSize;
593 :
594 :
595 20 : pafDstScanline = (float *) VSIMalloc((nDstXOff2 - nDstXOff) * sizeof(float));
596 20 : if( pafDstScanline == NULL )
597 : {
598 : CPLError( CE_Failure, CPLE_OutOfMemory,
599 0 : "GDALDownsampleChunk32R: Out of memory for line buffer." );
600 0 : return CE_Failure;
601 : }
602 :
603 : /* -------------------------------------------------------------------- */
604 : /* Figure out the line to start writing to, and the first line */
605 : /* to not write to. In theory this approach should ensure that */
606 : /* every output line will be written if all input chunks are */
607 : /* processed. */
608 : /* -------------------------------------------------------------------- */
609 20 : nDstYOff = (int) (0.5 + (nChunkYOff/(double)nSrcHeight) * nOYSize);
610 : nDstYOff2 = (int)
611 20 : (0.5 + ((nChunkYOff+nChunkYSize)/(double)nSrcHeight) * nOYSize);
612 :
613 20 : if( nChunkYOff + nChunkYSize == nSrcHeight )
614 20 : nDstYOff2 = nOYSize;
615 :
616 :
617 20 : int nEntryCount = 0;
618 20 : GDALColorEntry* aEntries = NULL;
619 20 : if (poColorTable)
620 : {
621 : int i;
622 2 : nEntryCount = poColorTable->GetColorEntryCount();
623 2 : aEntries = (GDALColorEntry* )CPLMalloc(sizeof(GDALColorEntry) * nEntryCount);
624 514 : for(i=0;i<nEntryCount;i++)
625 : {
626 512 : poColorTable->GetColorEntryAsRGB(i, &aEntries[i]);
627 : }
628 : }
629 :
630 20 : int nChunkRightXOff = MIN(nSrcWidth, nChunkXOff + nChunkXSize);
631 :
632 : /* ==================================================================== */
633 : /* Loop over destination scanlines. */
634 : /* ==================================================================== */
635 280 : for( int iDstLine = nDstYOff; iDstLine < nDstYOff2 && eErr == CE_None; iDstLine++ )
636 : {
637 : float *pafSrcScanline;
638 : GByte *pabySrcScanlineNodataMask;
639 260 : int nSrcYOff, nSrcYOff2 = 0, iDstPixel;
640 :
641 260 : nSrcYOff = (int) (0.5 + (iDstLine/(double)nOYSize) * nSrcHeight);
642 260 : nSrcYOff2 = (int) (0.5 + ((iDstLine+1)/(double)nOYSize) * nSrcHeight) + 1;
643 :
644 260 : if( nSrcYOff < nChunkYOff )
645 : {
646 0 : nSrcYOff = nChunkYOff;
647 0 : nSrcYOff2++;
648 : }
649 :
650 260 : int iSizeY = nSrcYOff2 - nSrcYOff;
651 260 : nSrcYOff = nSrcYOff + iSizeY/2 - nGaussMatrixDim/2;
652 260 : nSrcYOff2 = nSrcYOff + nGaussMatrixDim;
653 260 : if(nSrcYOff < 0)
654 0 : nSrcYOff = 0;
655 :
656 :
657 260 : if( nSrcYOff2 > nSrcHeight || iDstLine == nOYSize-1 )
658 20 : nSrcYOff2 = nSrcHeight;
659 260 : if( nSrcYOff2 > nChunkYOff + nChunkYSize)
660 0 : nSrcYOff2 = nChunkYOff + nChunkYSize;
661 :
662 260 : pafSrcScanline = pafChunk + ((nSrcYOff-nChunkYOff) * nChunkXSize);
663 260 : if (pabyChunkNodataMask != NULL)
664 150 : pabySrcScanlineNodataMask = pabyChunkNodataMask + ((nSrcYOff-nChunkYOff) * nChunkXSize);
665 : else
666 110 : pabySrcScanlineNodataMask = NULL;
667 :
668 : /* -------------------------------------------------------------------- */
669 : /* Loop over destination pixels */
670 : /* -------------------------------------------------------------------- */
671 4960 : for( iDstPixel = nDstXOff; iDstPixel < nDstXOff2; iDstPixel++ )
672 : {
673 : int nSrcXOff, nSrcXOff2;
674 :
675 4700 : nSrcXOff = (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
676 4700 : nSrcXOff2 = (int)(0.5 + ((iDstPixel+1)/(double)nOXSize) * nSrcWidth) + 1;
677 :
678 4700 : int iSizeX = nSrcXOff2 - nSrcXOff;
679 4700 : nSrcXOff = nSrcXOff + iSizeX/2 - nGaussMatrixDim/2;
680 4700 : nSrcXOff2 = nSrcXOff + nGaussMatrixDim;
681 4700 : if(nSrcXOff < 0)
682 0 : nSrcXOff = 0;
683 :
684 4700 : if( nSrcXOff2 > nChunkRightXOff || iDstPixel == nOXSize-1 )
685 260 : nSrcXOff2 = nChunkRightXOff;
686 :
687 4700 : if (poColorTable == NULL)
688 : {
689 4500 : double dfTotal = 0.0, val;
690 4500 : int nCount = 0, iX, iY;
691 4500 : int i = 0,j = 0;
692 4500 : const int *panLineWeight = panGaussMatrix;
693 :
694 17760 : for( j=0, iY = nSrcYOff; iY < nSrcYOff2;
695 : iY++, j++, panLineWeight += nGaussMatrixDim )
696 : {
697 52338 : for( i=0, iX = nSrcXOff; iX < nSrcXOff2; iX++,++i )
698 : {
699 39078 : val = pafSrcScanline[iX-nChunkXOff+(iY-nSrcYOff)*nChunkXSize];
700 71934 : if (pabySrcScanlineNodataMask == NULL ||
701 32856 : pabySrcScanlineNodataMask[iX-nChunkXOff+(iY-nSrcYOff)*nChunkXSize])
702 : {
703 18372 : int nWeight = panLineWeight[i];
704 18372 : dfTotal += val * nWeight;
705 18372 : nCount += nWeight;
706 : }
707 : }
708 : }
709 :
710 6714 : if (bHasNoData && nCount == 0)
711 : {
712 2214 : pafDstScanline[iDstPixel - nDstXOff] = fNoDataValue;
713 : }
714 : else
715 : {
716 2286 : if( nCount == 0 )
717 0 : pafDstScanline[iDstPixel - nDstXOff] = 0.0;
718 : else
719 2286 : pafDstScanline[iDstPixel - nDstXOff] = (float) (dfTotal / nCount);
720 : }
721 : }
722 : else
723 : {
724 : double val;
725 200 : int nTotalR = 0, nTotalG = 0, nTotalB = 0;
726 200 : int nTotalWeight = 0, iX, iY;
727 200 : int i = 0,j = 0;
728 200 : const int *panLineWeight = panGaussMatrix;
729 :
730 780 : for( j=0, iY = nSrcYOff; iY < nSrcYOff2;
731 : iY++, j++, panLineWeight += nGaussMatrixDim )
732 : {
733 2262 : for( i=0, iX = nSrcXOff; iX < nSrcXOff2; iX++,++i )
734 : {
735 1682 : val = pafSrcScanline[iX-nChunkXOff+(iY-nSrcYOff)*nChunkXSize];
736 1682 : if (bHasNoData == FALSE || val != fNoDataValue)
737 : {
738 1682 : int nVal = (int)val;
739 1682 : if (nVal >= 0 && nVal < nEntryCount)
740 : {
741 1682 : int nWeight = panLineWeight[i];
742 1682 : nTotalR += aEntries[nVal].c1 * nWeight;
743 1682 : nTotalG += aEntries[nVal].c2 * nWeight;
744 1682 : nTotalB += aEntries[nVal].c3 * nWeight;
745 1682 : nTotalWeight += nWeight;
746 : }
747 : }
748 : }
749 : }
750 :
751 200 : if (bHasNoData && nTotalWeight == 0)
752 : {
753 0 : pafDstScanline[iDstPixel - nDstXOff] = fNoDataValue;
754 : }
755 : else
756 : {
757 200 : if( nTotalWeight == 0 )
758 0 : pafDstScanline[iDstPixel - nDstXOff] = 0.0;
759 : else
760 : {
761 200 : int nR = nTotalR / nTotalWeight, nG = nTotalG / nTotalWeight, nB = nTotalB / nTotalWeight;
762 : int i;
763 200 : double dfMinDist = 0;
764 200 : int iBestEntry = 0;
765 51400 : for(i=0;i<nEntryCount;i++)
766 : {
767 102400 : double dfDist = (nR - aEntries[i].c1) * (nR - aEntries[i].c1) +
768 102400 : (nG - aEntries[i].c2) * (nG - aEntries[i].c2) +
769 204800 : (nB - aEntries[i].c3) * (nB - aEntries[i].c3);
770 51200 : if (i == 0 || dfDist < dfMinDist)
771 : {
772 400 : dfMinDist = dfDist;
773 400 : iBestEntry = i;
774 : }
775 : }
776 200 : pafDstScanline[iDstPixel - nDstXOff] =
777 200 : (float) iBestEntry;
778 : }
779 : }
780 : }
781 :
782 : }
783 :
784 : eErr = poOverview->RasterIO( GF_Write, nDstXOff, iDstLine, nDstXOff2 - nDstXOff, 1,
785 : pafDstScanline, nDstXOff2 - nDstXOff, 1, GDT_Float32,
786 260 : 0, 0 );
787 : }
788 :
789 20 : CPLFree( pafDstScanline );
790 20 : CPLFree( aEntries );
791 :
792 20 : return eErr;
793 : }
794 :
795 : /************************************************************************/
796 : /* GDALDownsampleChunk32R_Mode() */
797 : /************************************************************************/
798 :
799 : static CPLErr
800 18 : GDALDownsampleChunk32R_Mode( int nSrcWidth, int nSrcHeight,
801 : GDALDataType eWrkDataType,
802 : void * pChunk,
803 : GByte * pabyChunkNodataMask,
804 : int nChunkXOff, int nChunkXSize,
805 : int nChunkYOff, int nChunkYSize,
806 : GDALRasterBand * poOverview,
807 : const char * pszResampling,
808 : int bHasNoData, float fNoDataValue,
809 : GDALColorTable* poColorTable,
810 : GDALDataType eSrcDataType)
811 :
812 : {
813 18 : CPLErr eErr = CE_None;
814 :
815 18 : float * pafChunk = (float*) pChunk;
816 :
817 : /* -------------------------------------------------------------------- */
818 : /* Create the filter kernel and allocate scanline buffer. */
819 : /* -------------------------------------------------------------------- */
820 : int nDstXOff, nDstXOff2, nDstYOff, nDstYOff2, nOXSize, nOYSize;
821 : float *pafDstScanline;
822 :
823 18 : nOXSize = poOverview->GetXSize();
824 18 : nOYSize = poOverview->GetYSize();
825 :
826 : /* -------------------------------------------------------------------- */
827 : /* Figure out the column to start writing to, and the first column */
828 : /* to not write to. */
829 : /* -------------------------------------------------------------------- */
830 18 : nDstXOff = (int) (0.5 + (nChunkXOff/(double)nSrcWidth) * nOXSize);
831 : nDstXOff2 = (int)
832 18 : (0.5 + ((nChunkXOff+nChunkXSize)/(double)nSrcWidth) * nOXSize);
833 :
834 18 : if( nChunkXOff + nChunkXSize == nSrcWidth )
835 18 : nDstXOff2 = nOXSize;
836 :
837 :
838 18 : pafDstScanline = (float *) VSIMalloc((nDstXOff2 - nDstXOff) * sizeof(float));
839 18 : if( pafDstScanline == NULL )
840 : {
841 : CPLError( CE_Failure, CPLE_OutOfMemory,
842 0 : "GDALDownsampleChunk32R: Out of memory for line buffer." );
843 0 : return CE_Failure;
844 : }
845 :
846 : /* -------------------------------------------------------------------- */
847 : /* Figure out the line to start writing to, and the first line */
848 : /* to not write to. In theory this approach should ensure that */
849 : /* every output line will be written if all input chunks are */
850 : /* processed. */
851 : /* -------------------------------------------------------------------- */
852 18 : nDstYOff = (int) (0.5 + (nChunkYOff/(double)nSrcHeight) * nOYSize);
853 : nDstYOff2 = (int)
854 18 : (0.5 + ((nChunkYOff+nChunkYSize)/(double)nSrcHeight) * nOYSize);
855 :
856 18 : if( nChunkYOff + nChunkYSize == nSrcHeight )
857 18 : nDstYOff2 = nOYSize;
858 :
859 :
860 18 : int nEntryCount = 0;
861 18 : GDALColorEntry* aEntries = NULL;
862 18 : if (poColorTable)
863 : {
864 : int i;
865 2 : nEntryCount = poColorTable->GetColorEntryCount();
866 2 : aEntries = (GDALColorEntry* )CPLMalloc(sizeof(GDALColorEntry) * nEntryCount);
867 514 : for(i=0;i<nEntryCount;i++)
868 : {
869 512 : poColorTable->GetColorEntryAsRGB(i, &aEntries[i]);
870 : }
871 : }
872 :
873 18 : int nMaxNumPx = 0;
874 18 : float* pafVals = NULL;
875 18 : int* panSums = NULL;
876 :
877 18 : int nChunkRightXOff = MIN(nSrcWidth, nChunkXOff + nChunkXSize);
878 :
879 : /* ==================================================================== */
880 : /* Loop over destination scanlines. */
881 : /* ==================================================================== */
882 158 : for( int iDstLine = nDstYOff; iDstLine < nDstYOff2 && eErr == CE_None; iDstLine++ )
883 : {
884 : float *pafSrcScanline;
885 : GByte *pabySrcScanlineNodataMask;
886 140 : int nSrcYOff, nSrcYOff2 = 0, iDstPixel;
887 :
888 140 : nSrcYOff = (int) (0.5 + (iDstLine/(double)nOYSize) * nSrcHeight);
889 140 : if ( nSrcYOff < nChunkYOff )
890 0 : nSrcYOff = nChunkYOff;
891 :
892 : nSrcYOff2 =
893 140 : (int) (0.5 + ((iDstLine+1)/(double)nOYSize) * nSrcHeight);
894 :
895 140 : if( nSrcYOff2 > nSrcHeight || iDstLine == nOYSize-1 )
896 18 : nSrcYOff2 = nSrcHeight;
897 140 : if( nSrcYOff2 > nChunkYOff + nChunkYSize)
898 0 : nSrcYOff2 = nChunkYOff + nChunkYSize;
899 :
900 140 : pafSrcScanline = pafChunk + ((nSrcYOff-nChunkYOff) * nChunkXSize);
901 140 : if (pabyChunkNodataMask != NULL)
902 0 : pabySrcScanlineNodataMask = pabyChunkNodataMask + ((nSrcYOff-nChunkYOff) * nChunkXSize);
903 : else
904 140 : pabySrcScanlineNodataMask = NULL;
905 :
906 : /* -------------------------------------------------------------------- */
907 : /* Loop over destination pixels */
908 : /* -------------------------------------------------------------------- */
909 1340 : for( iDstPixel = nDstXOff; iDstPixel < nDstXOff2; iDstPixel++ )
910 : {
911 : int nSrcXOff, nSrcXOff2;
912 :
913 : nSrcXOff =
914 1200 : (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
915 1200 : if ( nSrcXOff < nChunkXOff )
916 0 : nSrcXOff = nChunkXOff;
917 : nSrcXOff2 = (int)
918 1200 : (0.5 + ((iDstPixel+1)/(double)nOXSize) * nSrcWidth);
919 :
920 1200 : if( nSrcXOff2 > nChunkRightXOff || iDstPixel == nOXSize-1 )
921 140 : nSrcXOff2 = nChunkRightXOff;
922 :
923 1950 : if (eSrcDataType != GDT_Byte || nEntryCount > 256)
924 : {
925 : /* I'm not sure how much sense it makes to run a majority
926 : filter on floating point data, but here it is for the sake
927 : of compatability. It won't look right on RGB images by the
928 : nature of the filter. */
929 750 : int nNumPx = (nSrcYOff2-nSrcYOff)*(nSrcXOff2-nSrcXOff);
930 750 : int iMaxInd = 0, iMaxVal = -1, iY, iX;
931 :
932 750 : if (nNumPx > nMaxNumPx)
933 : {
934 12 : pafVals = (float*) CPLRealloc(pafVals, nNumPx * sizeof(float));
935 12 : panSums = (int*) CPLRealloc(panSums, nNumPx * sizeof(int));
936 12 : nMaxNumPx = nNumPx;
937 : }
938 :
939 2550 : for( iY = nSrcYOff; iY < nSrcYOff2; ++iY )
940 : {
941 1800 : int iTotYOff = (iY-nSrcYOff)*nChunkXSize-nChunkXOff;
942 6600 : for( iX = nSrcXOff; iX < nSrcXOff2; ++iX )
943 : {
944 4800 : if (pabySrcScanlineNodataMask == NULL ||
945 0 : pabySrcScanlineNodataMask[iX+iTotYOff])
946 : {
947 4800 : float fVal = pafSrcScanline[iX+iTotYOff];
948 : int i;
949 :
950 : //Check array for existing entry
951 23148 : for( i = 0; i < iMaxInd; ++i )
952 23436 : if( pafVals[i] == fVal
953 4668 : && ++panSums[i] > panSums[iMaxVal] )
954 : {
955 420 : iMaxVal = i;
956 420 : break;
957 : }
958 :
959 : //Add to arr if entry not already there
960 4800 : if( i == iMaxInd )
961 : {
962 4380 : pafVals[iMaxInd] = fVal;
963 4380 : panSums[iMaxInd] = 1;
964 :
965 4380 : if( iMaxVal < 0 )
966 750 : iMaxVal = iMaxInd;
967 :
968 4380 : ++iMaxInd;
969 : }
970 : }
971 : }
972 : }
973 :
974 750 : if( iMaxVal == -1 )
975 0 : pafDstScanline[iDstPixel - nDstXOff] = fNoDataValue;
976 : else
977 750 : pafDstScanline[iDstPixel - nDstXOff] = pafVals[iMaxVal];
978 : }
979 : else /* if (eSrcDataType == GDT_Byte && nEntryCount < 256) */
980 : {
981 : /* So we go here for a paletted or non-paletted byte band */
982 : /* The input values are then between 0 and 255 */
983 450 : int anVals[256], nMaxVal = 0, iMaxInd = -1, iY, iX;
984 :
985 450 : memset(anVals, 0, 256*sizeof(int));
986 :
987 1450 : for( iY = nSrcYOff; iY < nSrcYOff2; ++iY )
988 : {
989 1000 : int iTotYOff = (iY-nSrcYOff)*nChunkXSize-nChunkXOff;
990 3400 : for( iX = nSrcXOff; iX < nSrcXOff2; ++iX )
991 : {
992 2400 : float val = pafSrcScanline[iX+iTotYOff];
993 2400 : if (bHasNoData == FALSE || val != fNoDataValue)
994 : {
995 2400 : int nVal = (int) val;
996 2400 : if ( ++anVals[nVal] > nMaxVal)
997 : {
998 : //Sum the density
999 : //Is it the most common value so far?
1000 948 : iMaxInd = nVal;
1001 948 : nMaxVal = anVals[nVal];
1002 : }
1003 : }
1004 : }
1005 : }
1006 :
1007 450 : if( iMaxInd == -1 )
1008 0 : pafDstScanline[iDstPixel - nDstXOff] = fNoDataValue;
1009 : else
1010 450 : pafDstScanline[iDstPixel - nDstXOff] = (float)iMaxInd;
1011 : }
1012 : }
1013 :
1014 : eErr = poOverview->RasterIO( GF_Write, nDstXOff, iDstLine, nDstXOff2 - nDstXOff, 1,
1015 : pafDstScanline, nDstXOff2 - nDstXOff, 1, GDT_Float32,
1016 140 : 0, 0 );
1017 : }
1018 :
1019 18 : CPLFree( pafDstScanline );
1020 18 : CPLFree( aEntries );
1021 18 : CPLFree( pafVals );
1022 18 : CPLFree( panSums );
1023 :
1024 18 : return eErr;
1025 : }
1026 :
1027 : /************************************************************************/
1028 : /* GDALDownsampleChunk32R_Cubic() */
1029 : /************************************************************************/
1030 :
1031 : static CPLErr
1032 48 : GDALDownsampleChunk32R_Cubic( int nSrcWidth, int nSrcHeight,
1033 : GDALDataType eWrkDataType,
1034 : void * pChunk,
1035 : GByte * pabyChunkNodataMask,
1036 : int nChunkXOff, int nChunkXSize,
1037 : int nChunkYOff, int nChunkYSize,
1038 : GDALRasterBand * poOverview,
1039 : const char * pszResampling,
1040 : int bHasNoData, float fNoDataValue,
1041 : GDALColorTable* poColorTable,
1042 : GDALDataType eSrcDataType)
1043 :
1044 : {
1045 :
1046 48 : CPLErr eErr = CE_None;
1047 :
1048 48 : float * pafChunk = (float*) pChunk;
1049 :
1050 : /* -------------------------------------------------------------------- */
1051 : /* Create the filter kernel and allocate scanline buffer. */
1052 : /* -------------------------------------------------------------------- */
1053 : int nDstXOff, nDstXOff2, nDstYOff, nDstYOff2, nOXSize, nOYSize;
1054 : float *pafDstScanline;
1055 :
1056 48 : nOXSize = poOverview->GetXSize();
1057 48 : nOYSize = poOverview->GetYSize();
1058 :
1059 : /* -------------------------------------------------------------------- */
1060 : /* Figure out the column to start writing to, and the first column */
1061 : /* to not write to. */
1062 : /* -------------------------------------------------------------------- */
1063 48 : nDstXOff = (int) (0.5 + (nChunkXOff/(double)nSrcWidth) * nOXSize);
1064 : nDstXOff2 = (int)
1065 48 : (0.5 + ((nChunkXOff+nChunkXSize)/(double)nSrcWidth) * nOXSize);
1066 :
1067 48 : if( nChunkXOff + nChunkXSize == nSrcWidth )
1068 48 : nDstXOff2 = nOXSize;
1069 :
1070 :
1071 48 : pafDstScanline = (float *) VSIMalloc((nDstXOff2 - nDstXOff) * sizeof(float));
1072 48 : if( pafDstScanline == NULL )
1073 : {
1074 : CPLError( CE_Failure, CPLE_OutOfMemory,
1075 0 : "GDALDownsampleChunk32R: Out of memory for line buffer." );
1076 0 : return CE_Failure;
1077 : }
1078 :
1079 : /* -------------------------------------------------------------------- */
1080 : /* Figure out the line to start writing to, and the first line */
1081 : /* to not write to. In theory this approach should ensure that */
1082 : /* every output line will be written if all input chunks are */
1083 : /* processed. */
1084 : /* -------------------------------------------------------------------- */
1085 48 : nDstYOff = (int) (0.5 + (nChunkYOff/(double)nSrcHeight) * nOYSize);
1086 : nDstYOff2 = (int)
1087 48 : (0.5 + ((nChunkYOff+nChunkYSize)/(double)nSrcHeight) * nOYSize);
1088 :
1089 48 : if( nChunkYOff + nChunkYSize == nSrcHeight )
1090 16 : nDstYOff2 = nOYSize;
1091 :
1092 :
1093 48 : int nEntryCount = 0;
1094 48 : GDALColorEntry* aEntries = NULL;
1095 48 : if (poColorTable)
1096 : {
1097 : int i;
1098 0 : nEntryCount = poColorTable->GetColorEntryCount();
1099 0 : aEntries = (GDALColorEntry* )CPLMalloc(sizeof(GDALColorEntry) * nEntryCount);
1100 0 : for(i=0;i<nEntryCount;i++)
1101 : {
1102 0 : poColorTable->GetColorEntryAsRGB(i, &aEntries[i]);
1103 : }
1104 : }
1105 :
1106 48 : int nChunkRightXOff = MIN(nSrcWidth, nChunkXOff + nChunkXSize);
1107 :
1108 : /* ==================================================================== */
1109 : /* Loop over destination scanlines. */
1110 : /* ==================================================================== */
1111 888 : for( int iDstLine = nDstYOff; iDstLine < nDstYOff2 && eErr == CE_None; iDstLine++ )
1112 : {
1113 : float *pafSrcScanline;
1114 : GByte *pabySrcScanlineNodataMask;
1115 840 : int nSrcYOff, nSrcYOff2 = 0, iDstPixel;
1116 :
1117 840 : nSrcYOff = (int) floor(((iDstLine+0.5)/(double)nOYSize) * nSrcHeight - 0.5)-1;
1118 840 : nSrcYOff2 = nSrcYOff + 4;
1119 840 : if(nSrcYOff < 0)
1120 8 : nSrcYOff = 0;
1121 840 : if(nSrcYOff < nChunkYOff)
1122 16 : nSrcYOff = nChunkYOff;
1123 :
1124 840 : if( nSrcYOff2 > nSrcHeight || iDstLine == nOYSize-1 )
1125 16 : nSrcYOff2 = nSrcHeight;
1126 840 : if( nSrcYOff2 > nChunkYOff + nChunkYSize)
1127 32 : nSrcYOff2 = nChunkYOff + nChunkYSize;
1128 :
1129 840 : pafSrcScanline = pafChunk + ((nSrcYOff-nChunkYOff) * nChunkXSize);
1130 840 : if (pabyChunkNodataMask != NULL)
1131 0 : pabySrcScanlineNodataMask = pabyChunkNodataMask + ((nSrcYOff-nChunkYOff) * nChunkXSize);
1132 : else
1133 840 : pabySrcScanlineNodataMask = NULL;
1134 :
1135 : /* -------------------------------------------------------------------- */
1136 : /* Loop over destination pixels */
1137 : /* -------------------------------------------------------------------- */
1138 57360 : for( iDstPixel = nDstXOff; iDstPixel < nDstXOff2; iDstPixel++ )
1139 : {
1140 : int nSrcXOff, nSrcXOff2;
1141 :
1142 56520 : nSrcXOff = (int) floor(((iDstPixel+0.5)/(double)nOXSize) * nSrcWidth - 0.5)-1;
1143 56520 : nSrcXOff2 = nSrcXOff + 4;
1144 :
1145 56520 : if(nSrcXOff < 0)
1146 600 : nSrcXOff = 0;
1147 :
1148 56520 : if( nSrcXOff2 > nChunkRightXOff || iDstPixel == nOXSize-1 )
1149 840 : nSrcXOff2 = nChunkRightXOff;
1150 :
1151 : // If we do not seem to have our full 4x4 window just
1152 : // do nearest resampling.
1153 62040 : if( nSrcXOff2 - nSrcXOff != 4 || nSrcYOff2 - nSrcYOff != 4 )
1154 : {
1155 5520 : int nLSrcYOff = (int) (0.5+(iDstLine/(double)nOYSize) * nSrcHeight);
1156 5520 : int nLSrcXOff = (int) (0.5+(iDstPixel/(double)nOXSize) * nSrcWidth);
1157 :
1158 5520 : if( nLSrcYOff < nChunkYOff )
1159 0 : nLSrcYOff = nChunkYOff;
1160 5520 : if( nLSrcYOff > nChunkYOff + nChunkYSize - 1 )
1161 0 : nLSrcYOff = nChunkYOff + nChunkYSize - 1;
1162 :
1163 5520 : pafDstScanline[iDstPixel - nDstXOff] =
1164 : pafChunk[(nLSrcYOff-nChunkYOff) * nChunkXSize
1165 5520 : + (nLSrcXOff - nChunkXOff)];
1166 : }
1167 : else
1168 : {
1169 : #define CubicConvolution(distance1,distance2,distance3,f0,f1,f2,f3) \
1170 : ( ( -f0 + f1 - f2 + f3) * distance3 \
1171 : + (2.0*(f0 - f1) + f2 - f3) * distance2 \
1172 : + ( -f0 + f2 ) * distance1 \
1173 : + f1 )
1174 :
1175 : int ic;
1176 : double adfRowResults[4];
1177 51000 : double dfSrcX = (((iDstPixel+0.5)/(double)nOXSize) * nSrcWidth);
1178 51000 : double dfDeltaX = dfSrcX - 0.5 - (nSrcXOff+1);
1179 51000 : double dfDeltaX2 = dfDeltaX * dfDeltaX;
1180 51000 : double dfDeltaX3 = dfDeltaX2 * dfDeltaX;
1181 :
1182 255000 : for ( ic = 0; ic < 4; ic++ )
1183 : {
1184 : float *pafSrcRow = pafSrcScanline +
1185 204000 : nSrcXOff-nChunkXOff+(nSrcYOff+ic-nSrcYOff)*nChunkXSize;
1186 :
1187 204000 : adfRowResults[ic] =
1188 1632000 : CubicConvolution(dfDeltaX, dfDeltaX2, dfDeltaX3,
1189 : pafSrcRow[0],
1190 : pafSrcRow[1],
1191 : pafSrcRow[2],
1192 1632000 : pafSrcRow[3] );
1193 : }
1194 :
1195 51000 : double dfSrcY = (((iDstLine+0.5)/(double)nOYSize) * nSrcHeight);
1196 51000 : double dfDeltaY = dfSrcY - 0.5 - (nSrcYOff+1);
1197 51000 : double dfDeltaY2 = dfDeltaY * dfDeltaY;
1198 51000 : double dfDeltaY3 = dfDeltaY2 * dfDeltaY;
1199 :
1200 51000 : pafDstScanline[iDstPixel - nDstXOff] = (float)
1201 408000 : CubicConvolution(dfDeltaY, dfDeltaY2, dfDeltaY3,
1202 : adfRowResults[0],
1203 : adfRowResults[1],
1204 : adfRowResults[2],
1205 408000 : adfRowResults[3] );
1206 : }
1207 : }
1208 :
1209 : eErr = poOverview->RasterIO( GF_Write, nDstXOff, iDstLine, nDstXOff2 - nDstXOff, 1,
1210 : pafDstScanline, nDstXOff2 - nDstXOff, 1, GDT_Float32,
1211 840 : 0, 0 );
1212 : }
1213 :
1214 48 : CPLFree( pafDstScanline );
1215 48 : CPLFree( aEntries );
1216 :
1217 48 : return eErr;
1218 : }
1219 :
1220 : /************************************************************************/
1221 : /* GDALDownsampleChunkC32R() */
1222 : /************************************************************************/
1223 :
1224 : static CPLErr
1225 8 : GDALDownsampleChunkC32R( int nSrcWidth, int nSrcHeight,
1226 : float * pafChunk, int nChunkYOff, int nChunkYSize,
1227 : GDALRasterBand * poOverview,
1228 : const char * pszResampling )
1229 :
1230 : {
1231 : int nDstYOff, nDstYOff2, nOXSize, nOYSize;
1232 : float *pafDstScanline;
1233 8 : CPLErr eErr = CE_None;
1234 :
1235 8 : nOXSize = poOverview->GetXSize();
1236 8 : nOYSize = poOverview->GetYSize();
1237 :
1238 8 : pafDstScanline = (float *) VSIMalloc(nOXSize * sizeof(float) * 2);
1239 8 : if( pafDstScanline == NULL )
1240 : {
1241 : CPLError( CE_Failure, CPLE_OutOfMemory,
1242 0 : "GDALDownsampleChunkC32R: Out of memory for line buffer." );
1243 0 : return CE_Failure;
1244 : }
1245 :
1246 : /* -------------------------------------------------------------------- */
1247 : /* Figure out the line to start writing to, and the first line */
1248 : /* to not write to. In theory this approach should ensure that */
1249 : /* every output line will be written if all input chunks are */
1250 : /* processed. */
1251 : /* -------------------------------------------------------------------- */
1252 8 : nDstYOff = (int) (0.5 + (nChunkYOff/(double)nSrcHeight) * nOYSize);
1253 : nDstYOff2 = (int)
1254 8 : (0.5 + ((nChunkYOff+nChunkYSize)/(double)nSrcHeight) * nOYSize);
1255 :
1256 8 : if( nChunkYOff + nChunkYSize == nSrcHeight )
1257 8 : nDstYOff2 = nOYSize;
1258 :
1259 : /* ==================================================================== */
1260 : /* Loop over destination scanlines. */
1261 : /* ==================================================================== */
1262 88 : for( int iDstLine = nDstYOff; iDstLine < nDstYOff2 && eErr == CE_None; iDstLine++ )
1263 : {
1264 : float *pafSrcScanline;
1265 : int nSrcYOff, nSrcYOff2, iDstPixel;
1266 :
1267 80 : nSrcYOff = (int) (0.5 + (iDstLine/(double)nOYSize) * nSrcHeight);
1268 80 : if( nSrcYOff < nChunkYOff )
1269 0 : nSrcYOff = nChunkYOff;
1270 :
1271 80 : nSrcYOff2 = (int) (0.5 + ((iDstLine+1)/(double)nOYSize) * nSrcHeight);
1272 80 : if( nSrcYOff2 > nSrcHeight || iDstLine == nOYSize-1 )
1273 8 : nSrcYOff2 = nSrcHeight;
1274 80 : if( nSrcYOff2 > nChunkYOff + nChunkYSize )
1275 0 : nSrcYOff2 = nChunkYOff + nChunkYSize;
1276 :
1277 80 : pafSrcScanline = pafChunk + ((nSrcYOff-nChunkYOff) * nSrcWidth) * 2;
1278 :
1279 : /* -------------------------------------------------------------------- */
1280 : /* Loop over destination pixels */
1281 : /* -------------------------------------------------------------------- */
1282 880 : for( iDstPixel = 0; iDstPixel < nOXSize; iDstPixel++ )
1283 : {
1284 : int nSrcXOff, nSrcXOff2;
1285 :
1286 800 : nSrcXOff = (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
1287 : nSrcXOff2 = (int)
1288 800 : (0.5 + ((iDstPixel+1)/(double)nOXSize) * nSrcWidth);
1289 800 : if( nSrcXOff2 > nSrcWidth )
1290 0 : nSrcXOff2 = nSrcWidth;
1291 :
1292 800 : if( EQUALN(pszResampling,"NEAR",4) )
1293 : {
1294 800 : pafDstScanline[iDstPixel*2] = pafSrcScanline[nSrcXOff*2];
1295 800 : pafDstScanline[iDstPixel*2+1] = pafSrcScanline[nSrcXOff*2+1];
1296 : }
1297 0 : else if( EQUAL(pszResampling,"AVERAGE_MAGPHASE") )
1298 : {
1299 0 : double dfTotalR = 0.0, dfTotalI = 0.0, dfTotalM = 0.0;
1300 0 : int nCount = 0, iX, iY;
1301 :
1302 0 : for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
1303 : {
1304 0 : for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
1305 : {
1306 : double dfR, dfI;
1307 :
1308 0 : dfR = pafSrcScanline[iX*2+(iY-nSrcYOff)*nSrcWidth*2];
1309 0 : dfI = pafSrcScanline[iX*2+(iY-nSrcYOff)*nSrcWidth*2+1];
1310 0 : dfTotalR += dfR;
1311 0 : dfTotalI += dfI;
1312 0 : dfTotalM += sqrt( dfR*dfR + dfI*dfI );
1313 0 : nCount++;
1314 : }
1315 : }
1316 :
1317 0 : CPLAssert( nCount > 0 );
1318 0 : if( nCount == 0 )
1319 : {
1320 0 : pafDstScanline[iDstPixel*2] = 0.0;
1321 0 : pafDstScanline[iDstPixel*2+1] = 0.0;
1322 : }
1323 : else
1324 : {
1325 0 : double dfM, dfDesiredM, dfRatio=1.0;
1326 :
1327 0 : pafDstScanline[iDstPixel*2 ] = (float) (dfTotalR/nCount);
1328 0 : pafDstScanline[iDstPixel*2+1] = (float) (dfTotalI/nCount);
1329 :
1330 0 : dfM = sqrt(pafDstScanline[iDstPixel*2 ]*pafDstScanline[iDstPixel*2 ]
1331 0 : + pafDstScanline[iDstPixel*2+1]*pafDstScanline[iDstPixel*2+1]);
1332 0 : dfDesiredM = dfTotalM / nCount;
1333 0 : if( dfM != 0.0 )
1334 0 : dfRatio = dfDesiredM / dfM;
1335 :
1336 0 : pafDstScanline[iDstPixel*2 ] *= (float) dfRatio;
1337 0 : pafDstScanline[iDstPixel*2+1] *= (float) dfRatio;
1338 : }
1339 : }
1340 0 : else if( EQUALN(pszResampling,"AVER",4) )
1341 : {
1342 0 : double dfTotalR = 0.0, dfTotalI = 0.0;
1343 0 : int nCount = 0, iX, iY;
1344 :
1345 0 : for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
1346 : {
1347 0 : for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
1348 : {
1349 0 : dfTotalR += pafSrcScanline[iX*2+(iY-nSrcYOff)*nSrcWidth*2];
1350 0 : dfTotalI += pafSrcScanline[iX*2+(iY-nSrcYOff)*nSrcWidth*2+1];
1351 0 : nCount++;
1352 : }
1353 : }
1354 :
1355 0 : CPLAssert( nCount > 0 );
1356 0 : if( nCount == 0 )
1357 : {
1358 0 : pafDstScanline[iDstPixel*2] = 0.0;
1359 0 : pafDstScanline[iDstPixel*2+1] = 0.0;
1360 : }
1361 : else
1362 : {
1363 0 : pafDstScanline[iDstPixel*2 ] = (float) (dfTotalR/nCount);
1364 0 : pafDstScanline[iDstPixel*2+1] = (float) (dfTotalI/nCount);
1365 : }
1366 : }
1367 : }
1368 :
1369 : eErr = poOverview->RasterIO( GF_Write, 0, iDstLine, nOXSize, 1,
1370 : pafDstScanline, nOXSize, 1, GDT_CFloat32,
1371 80 : 0, 0 );
1372 : }
1373 :
1374 8 : CPLFree( pafDstScanline );
1375 :
1376 8 : return eErr;
1377 : }
1378 :
1379 : /************************************************************************/
1380 : /* GDALRegenerateCascadingOverviews() */
1381 : /* */
1382 : /* Generate a list of overviews in order from largest to */
1383 : /* smallest, computing each from the next larger. */
1384 : /************************************************************************/
1385 :
1386 : static CPLErr
1387 22 : GDALRegenerateCascadingOverviews(
1388 : GDALRasterBand *poSrcBand, int nOverviews, GDALRasterBand **papoOvrBands,
1389 : const char * pszResampling,
1390 : GDALProgressFunc pfnProgress, void * pProgressData )
1391 :
1392 : {
1393 : /* -------------------------------------------------------------------- */
1394 : /* First, we must put the overviews in order from largest to */
1395 : /* smallest. */
1396 : /* -------------------------------------------------------------------- */
1397 : int i, j;
1398 :
1399 44 : for( i = 0; i < nOverviews-1; i++ )
1400 : {
1401 44 : for( j = 0; j < nOverviews - i - 1; j++ )
1402 : {
1403 :
1404 88 : if( papoOvrBands[j]->GetXSize()
1405 22 : * (float) papoOvrBands[j]->GetYSize() <
1406 22 : papoOvrBands[j+1]->GetXSize()
1407 22 : * (float) papoOvrBands[j+1]->GetYSize() )
1408 : {
1409 : GDALRasterBand * poTempBand;
1410 :
1411 0 : poTempBand = papoOvrBands[j];
1412 0 : papoOvrBands[j] = papoOvrBands[j+1];
1413 0 : papoOvrBands[j+1] = poTempBand;
1414 : }
1415 : }
1416 : }
1417 :
1418 : /* -------------------------------------------------------------------- */
1419 : /* Count total pixels so we can prepare appropriate scaled */
1420 : /* progress functions. */
1421 : /* -------------------------------------------------------------------- */
1422 22 : double dfTotalPixels = 0.0;
1423 :
1424 66 : for( i = 0; i < nOverviews; i++ )
1425 : {
1426 44 : dfTotalPixels += papoOvrBands[i]->GetXSize()
1427 88 : * (double) papoOvrBands[i]->GetYSize();
1428 : }
1429 :
1430 : /* -------------------------------------------------------------------- */
1431 : /* Generate all the bands. */
1432 : /* -------------------------------------------------------------------- */
1433 22 : double dfPixelsProcessed = 0.0;
1434 :
1435 66 : for( i = 0; i < nOverviews; i++ )
1436 : {
1437 : void *pScaledProgressData;
1438 : double dfPixels;
1439 : GDALRasterBand *poBaseBand;
1440 : CPLErr eErr;
1441 :
1442 44 : if( i == 0 )
1443 22 : poBaseBand = poSrcBand;
1444 : else
1445 22 : poBaseBand = papoOvrBands[i-1];
1446 :
1447 44 : dfPixels = papoOvrBands[i]->GetXSize()
1448 88 : * (double) papoOvrBands[i]->GetYSize();
1449 :
1450 : pScaledProgressData = GDALCreateScaledProgress(
1451 : dfPixelsProcessed / dfTotalPixels,
1452 : (dfPixelsProcessed + dfPixels) / dfTotalPixels,
1453 44 : pfnProgress, pProgressData );
1454 :
1455 : eErr = GDALRegenerateOverviews( (GDALRasterBandH) poBaseBand,
1456 : 1, (GDALRasterBandH *) papoOvrBands+i,
1457 : pszResampling,
1458 : GDALScaledProgress,
1459 44 : pScaledProgressData );
1460 44 : GDALDestroyScaledProgress( pScaledProgressData );
1461 :
1462 44 : if( eErr != CE_None )
1463 0 : return eErr;
1464 :
1465 44 : dfPixelsProcessed += dfPixels;
1466 :
1467 : /* we only do the bit2grayscale promotion on the base band */
1468 44 : if( EQUALN(pszResampling,"AVERAGE_BIT2GRAYSCALE",13) )
1469 8 : pszResampling = "AVERAGE";
1470 : }
1471 :
1472 22 : return CE_None;
1473 : }
1474 :
1475 : /************************************************************************/
1476 : /* GDALGetDownsampleFunction() */
1477 : /************************************************************************/
1478 :
1479 : static
1480 283 : GDALDownsampleFunction GDALGetDownsampleFunction(const char* pszResampling)
1481 : {
1482 283 : if( EQUALN(pszResampling,"NEAR",4) )
1483 160 : return GDALDownsampleChunk32R_Near;
1484 123 : else if( EQUALN(pszResampling,"AVER",4) )
1485 83 : return GDALDownsampleChunk32R_Average;
1486 40 : else if( EQUALN(pszResampling,"GAUSS",5) )
1487 22 : return GDALDownsampleChunk32R_Gauss;
1488 18 : else if( EQUALN(pszResampling,"MODE",4) )
1489 10 : return GDALDownsampleChunk32R_Mode;
1490 8 : else if( EQUALN(pszResampling,"CUBIC",5) )
1491 8 : return GDALDownsampleChunk32R_Cubic;
1492 : else
1493 : {
1494 : CPLError( CE_Failure, CPLE_AppDefined,
1495 : "GDALGetDownsampleFunction: Unsupported resampling method \"%s\".",
1496 0 : pszResampling );
1497 0 : return NULL;
1498 : }
1499 : }
1500 :
1501 : /************************************************************************/
1502 : /* GDALGetOvrWorkDataType() */
1503 : /************************************************************************/
1504 :
1505 253 : static GDALDataType GDALGetOvrWorkDataType(const char* pszResampling,
1506 : GDALDataType eSrcDataType)
1507 : {
1508 253 : if( (EQUALN(pszResampling,"NEAR",4) || EQUALN(pszResampling,"AVER",4)) &&
1509 : eSrcDataType == GDT_Byte)
1510 184 : return GDT_Byte;
1511 : else
1512 69 : return GDT_Float32;
1513 : }
1514 :
1515 : /************************************************************************/
1516 : /* GDALRegenerateOverviews() */
1517 : /************************************************************************/
1518 :
1519 : /**
1520 : * \brief Generate downsampled overviews.
1521 : *
1522 : * This function will generate one or more overview images from a base
1523 : * image using the requested downsampling algorithm. It's primary use
1524 : * is for generating overviews via GDALDataset::BuildOverviews(), but it
1525 : * can also be used to generate downsampled images in one file from another
1526 : * outside the overview architecture.
1527 : *
1528 : * The output bands need to exist in advance.
1529 : *
1530 : * The full set of resampling algorithms is documented in
1531 : * GDALDataset::BuildOverviews().
1532 : *
1533 : * This function will honour properly NODATA_VALUES tuples (special dataset metadata) so
1534 : * that only a given RGB triplet (in case of a RGB image) will be considered as the
1535 : * nodata value and not each value of the triplet independantly per band.
1536 : *
1537 : * @param hSrcBand the source (base level) band.
1538 : * @param nOverviewCount the number of downsampled bands being generated.
1539 : * @param pahOvrBands the list of downsampled bands to be generated.
1540 : * @param pszResampling Resampling algorithm (eg. "AVERAGE").
1541 : * @param pfnProgress progress report function.
1542 : * @param pProgressData progress function callback data.
1543 : * @return CE_None on success or CE_Failure on failure.
1544 : */
1545 : CPLErr
1546 268 : GDALRegenerateOverviews( GDALRasterBandH hSrcBand,
1547 : int nOverviewCount, GDALRasterBandH *pahOvrBands,
1548 : const char * pszResampling,
1549 : GDALProgressFunc pfnProgress, void * pProgressData )
1550 :
1551 : {
1552 268 : GDALRasterBand *poSrcBand = (GDALRasterBand *) hSrcBand;
1553 268 : GDALRasterBand **papoOvrBands = (GDALRasterBand **) pahOvrBands;
1554 : int nFullResYChunk, nWidth;
1555 : int nFRXBlockSize, nFRYBlockSize;
1556 : GDALDataType eType;
1557 : int bHasNoData;
1558 : float fNoDataValue;
1559 268 : GDALColorTable* poColorTable = NULL;
1560 :
1561 268 : if( pfnProgress == NULL )
1562 6 : pfnProgress = GDALDummyProgress;
1563 :
1564 268 : if( EQUAL(pszResampling,"NONE") )
1565 10 : return CE_None;
1566 :
1567 258 : GDALDownsampleFunction pfnDownsampleFn = GDALGetDownsampleFunction(pszResampling);
1568 258 : if (pfnDownsampleFn == NULL)
1569 0 : return CE_Failure;
1570 :
1571 : /* -------------------------------------------------------------------- */
1572 : /* Check color tables... */
1573 : /* -------------------------------------------------------------------- */
1574 365 : if ((EQUALN(pszResampling,"AVER",4)
1575 : || EQUALN(pszResampling,"MODE",4)
1576 : || EQUALN(pszResampling,"GAUSS",5)) &&
1577 107 : poSrcBand->GetColorInterpretation() == GCI_PaletteIndex)
1578 : {
1579 6 : poColorTable = poSrcBand->GetColorTable();
1580 6 : if (poColorTable != NULL)
1581 : {
1582 6 : if (poColorTable->GetPaletteInterpretation() != GPI_RGB)
1583 : {
1584 : CPLError(CE_Warning, CPLE_AppDefined,
1585 : "Computing overviews on palette index raster bands "
1586 : "with a palette whose color interpreation is not RGB "
1587 0 : "will probably lead to unexpected results.");
1588 0 : poColorTable = NULL;
1589 : }
1590 : }
1591 : else
1592 : {
1593 : CPLError(CE_Warning, CPLE_AppDefined,
1594 : "Computing overviews on palette index raster bands "
1595 0 : "without a palette will probably lead to unexpected results.");
1596 : }
1597 : }
1598 :
1599 :
1600 : /* If we have a nodata mask and we are doing something more complicated */
1601 : /* than nearest neighbouring, we have to fetch to nodata mask */
1602 : int bUseNoDataMask = (!EQUALN(pszResampling,"NEAR",4) &&
1603 258 : (poSrcBand->GetMaskFlags() & GMF_ALL_VALID) == 0);
1604 :
1605 : /* -------------------------------------------------------------------- */
1606 : /* If we are operating on multiple overviews, and using */
1607 : /* averaging, lets do them in cascading order to reduce the */
1608 : /* amount of computation. */
1609 : /* -------------------------------------------------------------------- */
1610 :
1611 : /* In case the mask made be computed from another band of the dataset, */
1612 : /* we can't use cascaded generation, as the computation of the overviews */
1613 : /* of the band used for the mask band may not have yet occured (#3033) */
1614 269 : if( (EQUALN(pszResampling,"AVER",4) || EQUALN(pszResampling,"GAUSS",5)) && nOverviewCount > 1
1615 11 : && !(bUseNoDataMask && poSrcBand->GetMaskFlags() != GMF_NODATA))
1616 : return GDALRegenerateCascadingOverviews( poSrcBand,
1617 : nOverviewCount, papoOvrBands,
1618 : pszResampling,
1619 : pfnProgress,
1620 22 : pProgressData );
1621 :
1622 : /* -------------------------------------------------------------------- */
1623 : /* Setup one horizontal swath to read from the raw buffer. */
1624 : /* -------------------------------------------------------------------- */
1625 : void *pChunk;
1626 236 : GByte *pabyChunkNodataMask = NULL;
1627 :
1628 236 : poSrcBand->GetBlockSize( &nFRXBlockSize, &nFRYBlockSize );
1629 :
1630 283 : if( nFRYBlockSize < 16 || nFRYBlockSize > 256 )
1631 47 : nFullResYChunk = 64;
1632 : else
1633 189 : nFullResYChunk = nFRYBlockSize;
1634 :
1635 236 : if( GDALDataTypeIsComplex( poSrcBand->GetRasterDataType() ) )
1636 8 : eType = GDT_CFloat32;
1637 : else
1638 228 : eType = GDALGetOvrWorkDataType(pszResampling, poSrcBand->GetRasterDataType());
1639 :
1640 236 : nWidth = poSrcBand->GetXSize();
1641 : pChunk =
1642 236 : VSIMalloc3((GDALGetDataTypeSize(eType)/8), nFullResYChunk, nWidth );
1643 236 : if (bUseNoDataMask)
1644 : {
1645 : pabyChunkNodataMask = (GByte *)
1646 22 : (GByte*) VSIMalloc2( nFullResYChunk, nWidth );
1647 : }
1648 :
1649 236 : if( pChunk == NULL || (bUseNoDataMask && pabyChunkNodataMask == NULL))
1650 : {
1651 0 : CPLFree(pChunk);
1652 0 : CPLFree(pabyChunkNodataMask);
1653 : CPLError( CE_Failure, CPLE_OutOfMemory,
1654 0 : "Out of memory in GDALRegenerateOverviews()." );
1655 :
1656 0 : return CE_Failure;
1657 : }
1658 :
1659 236 : fNoDataValue = (float) poSrcBand->GetNoDataValue(&bHasNoData);
1660 :
1661 : /* -------------------------------------------------------------------- */
1662 : /* Loop over image operating on chunks. */
1663 : /* -------------------------------------------------------------------- */
1664 236 : int nChunkYOff = 0;
1665 236 : CPLErr eErr = CE_None;
1666 :
1667 794 : for( nChunkYOff = 0;
1668 : nChunkYOff < poSrcBand->GetYSize() && eErr == CE_None;
1669 : nChunkYOff += nFullResYChunk )
1670 : {
1671 558 : if( !pfnProgress( nChunkYOff / (double) poSrcBand->GetYSize(),
1672 : NULL, pProgressData ) )
1673 : {
1674 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1675 0 : eErr = CE_Failure;
1676 : }
1677 :
1678 558 : if( nFullResYChunk + nChunkYOff > poSrcBand->GetYSize() )
1679 91 : nFullResYChunk = poSrcBand->GetYSize() - nChunkYOff;
1680 :
1681 : /* read chunk */
1682 558 : if (eErr == CE_None)
1683 : eErr = poSrcBand->RasterIO( GF_Read, 0, nChunkYOff, nWidth, nFullResYChunk,
1684 : pChunk, nWidth, nFullResYChunk, eType,
1685 558 : 0, 0 );
1686 558 : if (eErr == CE_None && bUseNoDataMask)
1687 46 : eErr = poSrcBand->GetMaskBand()->RasterIO( GF_Read, 0, nChunkYOff, nWidth, nFullResYChunk,
1688 : pabyChunkNodataMask, nWidth, nFullResYChunk, GDT_Byte,
1689 46 : 0, 0 );
1690 :
1691 : /* special case to promote 1bit data to 8bit 0/255 values */
1692 558 : if( EQUAL(pszResampling,"AVERAGE_BIT2GRAYSCALE") )
1693 : {
1694 : int i;
1695 :
1696 13 : if (eType == GDT_Float32)
1697 : {
1698 0 : float* pafChunk = (float*)pChunk;
1699 0 : for( i = nFullResYChunk*nWidth - 1; i >= 0; i-- )
1700 : {
1701 0 : if( pafChunk[i] == 1.0 )
1702 0 : pafChunk[i] = 255.0;
1703 : }
1704 : }
1705 13 : else if (eType == GDT_Byte)
1706 : {
1707 13 : GByte* pabyChunk = (GByte*)pChunk;
1708 168421 : for( i = nFullResYChunk*nWidth - 1; i >= 0; i-- )
1709 : {
1710 168408 : if( pabyChunk[i] == 1 )
1711 127437 : pabyChunk[i] = 255;
1712 : }
1713 : }
1714 : else
1715 0 : CPLAssert(0);
1716 : }
1717 545 : else if( EQUAL(pszResampling,"AVERAGE_BIT2GRAYSCALE_MINISWHITE") )
1718 : {
1719 : int i;
1720 :
1721 0 : if (eType == GDT_Float32)
1722 : {
1723 0 : float* pafChunk = (float*)pChunk;
1724 0 : for( i = nFullResYChunk*nWidth - 1; i >= 0; i-- )
1725 : {
1726 0 : if( pafChunk[i] == 1.0 )
1727 0 : pafChunk[i] = 0.0;
1728 0 : else if( pafChunk[i] == 0.0 )
1729 0 : pafChunk[i] = 255.0;
1730 : }
1731 : }
1732 0 : else if (eType == GDT_Byte)
1733 : {
1734 0 : GByte* pabyChunk = (GByte*)pChunk;
1735 0 : for( i = nFullResYChunk*nWidth - 1; i >= 0; i-- )
1736 : {
1737 0 : if( pabyChunk[i] == 1 )
1738 0 : pabyChunk[i] = 0;
1739 0 : else if( pabyChunk[i] == 0 )
1740 0 : pabyChunk[i] = 255;
1741 : }
1742 : }
1743 : else
1744 0 : CPLAssert(0);
1745 : }
1746 :
1747 1522 : for( int iOverview = 0; iOverview < nOverviewCount && eErr == CE_None; iOverview++ )
1748 : {
1749 1920 : if( eType == GDT_Byte || eType == GDT_Float32 )
1750 : eErr = pfnDownsampleFn(nWidth, poSrcBand->GetYSize(),
1751 : eType,
1752 : pChunk,
1753 : pabyChunkNodataMask,
1754 : 0, nWidth,
1755 : nChunkYOff, nFullResYChunk,
1756 : papoOvrBands[iOverview], pszResampling,
1757 : bHasNoData, fNoDataValue, poColorTable,
1758 956 : poSrcBand->GetRasterDataType());
1759 : else
1760 : eErr = GDALDownsampleChunkC32R(nWidth, poSrcBand->GetYSize(),
1761 : (float*)pChunk, nChunkYOff, nFullResYChunk,
1762 8 : papoOvrBands[iOverview], pszResampling);
1763 : }
1764 : }
1765 :
1766 236 : VSIFree( pChunk );
1767 236 : VSIFree( pabyChunkNodataMask );
1768 :
1769 : /* -------------------------------------------------------------------- */
1770 : /* Renormalized overview mean / stddev if needed. */
1771 : /* -------------------------------------------------------------------- */
1772 236 : if( eErr == CE_None && EQUAL(pszResampling,"AVERAGE_MP") )
1773 : {
1774 : GDALOverviewMagnitudeCorrection( (GDALRasterBandH) poSrcBand,
1775 : nOverviewCount,
1776 : (GDALRasterBandH *) papoOvrBands,
1777 0 : GDALDummyProgress, NULL );
1778 : }
1779 :
1780 : /* -------------------------------------------------------------------- */
1781 : /* It can be important to flush out data to overviews. */
1782 : /* -------------------------------------------------------------------- */
1783 570 : for( int iOverview = 0;
1784 : eErr == CE_None && iOverview < nOverviewCount;
1785 : iOverview++ )
1786 : {
1787 334 : eErr = papoOvrBands[iOverview]->FlushCache();
1788 : }
1789 :
1790 236 : if (eErr == CE_None)
1791 236 : pfnProgress( 1.0, NULL, pProgressData );
1792 :
1793 236 : return eErr;
1794 : }
1795 :
1796 :
1797 :
1798 : /************************************************************************/
1799 : /* GDALRegenerateOverviewsMultiBand() */
1800 : /************************************************************************/
1801 :
1802 : /**
1803 : * \brief Variant of GDALRegenerateOverviews, specialy dedicated for generating
1804 : * compressed pixel-interleaved overviews (JPEG-IN-TIFF for example)
1805 : *
1806 : * This function will generate one or more overview images from a base
1807 : * image using the requested downsampling algorithm. It's primary use
1808 : * is for generating overviews via GDALDataset::BuildOverviews(), but it
1809 : * can also be used to generate downsampled images in one file from another
1810 : * outside the overview architecture.
1811 : *
1812 : * The output bands need to exist in advance and share the same characteristics
1813 : * (type, dimensions)
1814 : *
1815 : * The resampling algorithms supported for the moment are "NEAREST", "AVERAGE"
1816 : * and "GAUSS"
1817 : *
1818 : * The pseudo-algorithm used by the function is :
1819 : * for each overview
1820 : * iterate on lines of the source by a step of deltay
1821 : * iterate on columns of the source by a step of deltax
1822 : * read the source data of size deltax * deltay for all the bands
1823 : * generate the corresponding overview block for all the bands
1824 : *
1825 : * This function will honour properly NODATA_VALUES tuples (special dataset metadata) so
1826 : * that only a given RGB triplet (in case of a RGB image) will be considered as the
1827 : * nodata value and not each value of the triplet independantly per band.
1828 : *
1829 : * @param nBands the number of bands, size of papoSrcBands and size of
1830 : * first dimension of papapoOverviewBands
1831 : * @param papoSrcBands the list of source bands to downsample
1832 : * @param nOverviews the number of downsampled overview levels being generated.
1833 : * @param papapoOverviewBands bidimension array of bands. First dimension is indexed
1834 : * by nBands. Second dimension is indexed by nOverviews.
1835 : * @param pszResampling Resampling algorithm ("NEAREST", "AVERAGE" or "GAUSS").
1836 : * @param pfnProgress progress report function.
1837 : * @param pProgressData progress function callback data.
1838 : * @return CE_None on success or CE_Failure on failure.
1839 : */
1840 :
1841 : CPLErr
1842 25 : GDALRegenerateOverviewsMultiBand(int nBands, GDALRasterBand** papoSrcBands,
1843 : int nOverviews,
1844 : GDALRasterBand*** papapoOverviewBands,
1845 : const char * pszResampling,
1846 : GDALProgressFunc pfnProgress, void * pProgressData )
1847 : {
1848 25 : CPLErr eErr = CE_None;
1849 : int iOverview, iBand;
1850 :
1851 25 : if( pfnProgress == NULL )
1852 0 : pfnProgress = GDALDummyProgress;
1853 :
1854 25 : if( EQUAL(pszResampling,"NONE") )
1855 0 : return CE_None;
1856 :
1857 : /* Sanity checks */
1858 25 : if (!EQUALN(pszResampling, "NEAR", 4) && !EQUAL(pszResampling, "AVERAGE") && !EQUAL(pszResampling, "GAUSS"))
1859 : {
1860 : CPLError(CE_Failure, CPLE_NotSupported,
1861 0 : "GDALRegenerateOverviewsMultiBand: pszResampling='%s' not supported", pszResampling);
1862 0 : return CE_Failure;
1863 : }
1864 :
1865 25 : GDALDownsampleFunction pfnDownsampleFn = GDALGetDownsampleFunction(pszResampling);
1866 25 : if (pfnDownsampleFn == NULL)
1867 0 : return CE_Failure;
1868 :
1869 25 : int nSrcWidth = papoSrcBands[0]->GetXSize();
1870 25 : int nSrcHeight = papoSrcBands[0]->GetYSize();
1871 25 : GDALDataType eDataType = papoSrcBands[0]->GetRasterDataType();
1872 63 : for(iBand=1;iBand<nBands;iBand++)
1873 : {
1874 76 : if (papoSrcBands[iBand]->GetXSize() != nSrcWidth ||
1875 38 : papoSrcBands[iBand]->GetYSize() != nSrcHeight)
1876 : {
1877 : CPLError(CE_Failure, CPLE_NotSupported,
1878 0 : "GDALRegenerateOverviewsMultiBand: all the source bands must have the same dimensions");
1879 0 : return CE_Failure;
1880 : }
1881 38 : if (papoSrcBands[iBand]->GetRasterDataType() != eDataType)
1882 : {
1883 : CPLError(CE_Failure, CPLE_NotSupported,
1884 0 : "GDALRegenerateOverviewsMultiBand: all the source bands must have the same data type");
1885 0 : return CE_Failure;
1886 : }
1887 : }
1888 :
1889 63 : for(iOverview=0;iOverview<nOverviews;iOverview++)
1890 : {
1891 38 : int nDstWidth = papapoOverviewBands[0][iOverview]->GetXSize();
1892 38 : int nDstHeight = papapoOverviewBands[0][iOverview]->GetYSize();
1893 98 : for(iBand=1;iBand<nBands;iBand++)
1894 : {
1895 180 : if (papapoOverviewBands[iBand][iOverview]->GetXSize() != nDstWidth ||
1896 120 : papapoOverviewBands[iBand][iOverview]->GetYSize() != nDstHeight)
1897 : {
1898 : CPLError(CE_Failure, CPLE_NotSupported,
1899 0 : "GDALRegenerateOverviewsMultiBand: all the overviews bands of the same level must have the same dimensions");
1900 0 : return CE_Failure;
1901 : }
1902 60 : if (papapoOverviewBands[iBand][iOverview]->GetRasterDataType() != eDataType)
1903 : {
1904 : CPLError(CE_Failure, CPLE_NotSupported,
1905 0 : "GDALRegenerateOverviewsMultiBand: all the overviews bands must have the same data type as the source bands");
1906 0 : return CE_Failure;
1907 : }
1908 : }
1909 : }
1910 :
1911 : /* First pass to compute the total number of pixels to read */
1912 25 : double dfTotalPixelCount = 0;
1913 63 : for(iOverview=0;iOverview<nOverviews;iOverview++)
1914 : {
1915 38 : nSrcWidth = papoSrcBands[0]->GetXSize();
1916 38 : nSrcHeight = papoSrcBands[0]->GetYSize();
1917 :
1918 38 : int nDstWidth = papapoOverviewBands[0][iOverview]->GetXSize();
1919 : /* Try to use previous level of overview as the source to compute */
1920 : /* the next level */
1921 38 : if (iOverview > 0 && papapoOverviewBands[0][iOverview - 1]->GetXSize() > nDstWidth)
1922 : {
1923 13 : nSrcWidth = papapoOverviewBands[0][iOverview - 1]->GetXSize();
1924 13 : nSrcHeight = papapoOverviewBands[0][iOverview - 1]->GetYSize();
1925 : }
1926 :
1927 38 : dfTotalPixelCount += (double)nSrcWidth * nSrcHeight;
1928 : }
1929 :
1930 25 : nSrcWidth = papoSrcBands[0]->GetXSize();
1931 25 : nSrcHeight = papoSrcBands[0]->GetYSize();
1932 :
1933 25 : GDALDataType eWrkDataType = GDALGetOvrWorkDataType(pszResampling, eDataType);
1934 :
1935 : /* If we have a nodata mask and we are doing something more complicated */
1936 : /* than nearest neighbouring, we have to fetch to nodata mask */
1937 : int bUseNoDataMask = (!EQUALN(pszResampling,"NEAR",4) &&
1938 25 : (papoSrcBands[0]->GetMaskFlags() & GMF_ALL_VALID) == 0);
1939 :
1940 25 : int* pabHasNoData = (int*)CPLMalloc(nBands * sizeof(int));
1941 25 : float* pafNoDataValue = (float*)CPLMalloc(nBands * sizeof(float));
1942 :
1943 88 : for(iBand=0;iBand<nBands;iBand++)
1944 : {
1945 63 : pabHasNoData[iBand] = FALSE;
1946 63 : pafNoDataValue[iBand] = (float) papoSrcBands[iBand]->GetNoDataValue(&pabHasNoData[iBand]);
1947 : }
1948 :
1949 : /* Second pass to do the real job ! */
1950 25 : double dfCurPixelCount = 0;
1951 63 : for(iOverview=0;iOverview<nOverviews && eErr == CE_None;iOverview++)
1952 : {
1953 38 : int iSrcOverview = -1; /* -1 means the source bands */
1954 :
1955 : int nDstBlockXSize, nDstBlockYSize;
1956 : int nDstWidth, nDstHeight;
1957 38 : papapoOverviewBands[0][iOverview]->GetBlockSize(&nDstBlockXSize, &nDstBlockYSize);
1958 38 : nDstWidth = papapoOverviewBands[0][iOverview]->GetXSize();
1959 38 : nDstHeight = papapoOverviewBands[0][iOverview]->GetYSize();
1960 :
1961 : /* Try to use previous level of overview as the source to compute */
1962 : /* the next level */
1963 38 : if (iOverview > 0 && papapoOverviewBands[0][iOverview - 1]->GetXSize() > nDstWidth)
1964 : {
1965 13 : nSrcWidth = papapoOverviewBands[0][iOverview - 1]->GetXSize();
1966 13 : nSrcHeight = papapoOverviewBands[0][iOverview - 1]->GetYSize();
1967 13 : iSrcOverview = iOverview - 1;
1968 : }
1969 :
1970 : /* Compute the chunck size of the source such as it will match the size of */
1971 : /* a block of the overview */
1972 38 : int nFullResXChunk = (nDstBlockXSize * nSrcWidth) / nDstWidth;
1973 38 : int nFullResYChunk = (nDstBlockYSize * nSrcHeight) / nDstHeight;
1974 :
1975 38 : void** papaChunk = (void**) CPLMalloc(nBands * sizeof(void*));
1976 38 : GByte* pabyChunkNoDataMask = NULL;
1977 136 : for(iBand=0;iBand<nBands;iBand++)
1978 : {
1979 98 : papaChunk[iBand] = VSIMalloc3(nFullResXChunk, nFullResYChunk, GDALGetDataTypeSize(eWrkDataType) / 8);
1980 98 : if( papaChunk[iBand] == NULL )
1981 : {
1982 0 : while ( --iBand >= 0)
1983 0 : CPLFree(papaChunk[iBand]);
1984 0 : CPLFree(papaChunk);
1985 0 : CPLFree(pabHasNoData);
1986 0 : CPLFree(pafNoDataValue);
1987 :
1988 : CPLError( CE_Failure, CPLE_OutOfMemory,
1989 0 : "GDALRegenerateOverviewsMultiBand: Out of memory." );
1990 0 : return CE_Failure;
1991 : }
1992 : }
1993 38 : if (bUseNoDataMask)
1994 : {
1995 4 : pabyChunkNoDataMask = (GByte*) VSIMalloc2(nFullResXChunk, nFullResYChunk);
1996 4 : if( pabyChunkNoDataMask == NULL )
1997 : {
1998 0 : for(iBand=0;iBand<nBands;iBand++)
1999 : {
2000 0 : CPLFree(papaChunk[iBand]);
2001 : }
2002 0 : CPLFree(papaChunk);
2003 0 : CPLFree(pabHasNoData);
2004 0 : CPLFree(pafNoDataValue);
2005 :
2006 : CPLError( CE_Failure, CPLE_OutOfMemory,
2007 0 : "GDALRegenerateOverviewsMultiBand: Out of memory." );
2008 0 : return CE_Failure;
2009 : }
2010 : }
2011 :
2012 : int nChunkYOff;
2013 : /* Iterate on destination overview, block by block */
2014 120 : for( nChunkYOff = 0; nChunkYOff < nSrcHeight && eErr == CE_None; nChunkYOff += nFullResYChunk )
2015 : {
2016 : int nYCount;
2017 82 : if (nChunkYOff + nFullResYChunk <= nSrcHeight)
2018 66 : nYCount = nFullResYChunk;
2019 : else
2020 16 : nYCount = nSrcHeight - nChunkYOff;
2021 :
2022 82 : if( !pfnProgress( dfCurPixelCount / dfTotalPixelCount,
2023 : NULL, pProgressData ) )
2024 : {
2025 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
2026 0 : eErr = CE_Failure;
2027 : }
2028 :
2029 : int nChunkXOff;
2030 318 : for( nChunkXOff = 0; nChunkXOff < nSrcWidth && eErr == CE_None; nChunkXOff += nFullResXChunk )
2031 : {
2032 : int nXCount;
2033 236 : if (nChunkXOff + nFullResXChunk <= nSrcWidth)
2034 220 : nXCount = nFullResXChunk;
2035 : else
2036 16 : nXCount = nSrcWidth - nChunkXOff;
2037 :
2038 : /* Read the source buffers for all the bands */
2039 928 : for(iBand=0;iBand<nBands && eErr == CE_None;iBand++)
2040 : {
2041 : GDALRasterBand* poSrcBand;
2042 692 : if (iSrcOverview == -1)
2043 558 : poSrcBand = papoSrcBands[iBand];
2044 : else
2045 134 : poSrcBand = papapoOverviewBands[iBand][iSrcOverview];
2046 : eErr = poSrcBand->RasterIO( GF_Read,
2047 : nChunkXOff, nChunkYOff,
2048 : nXCount, nYCount,
2049 : papaChunk[iBand],
2050 : nXCount, nYCount,
2051 692 : eWrkDataType, 0, 0 );
2052 : }
2053 :
2054 236 : if (bUseNoDataMask && eErr == CE_None)
2055 : {
2056 : GDALRasterBand* poSrcBand;
2057 4 : if (iSrcOverview == -1)
2058 4 : poSrcBand = papoSrcBands[0];
2059 : else
2060 0 : poSrcBand = papapoOverviewBands[0][iSrcOverview];
2061 4 : eErr = poSrcBand->GetMaskBand()->RasterIO( GF_Read,
2062 : nChunkXOff, nChunkYOff,
2063 : nXCount, nYCount,
2064 : pabyChunkNoDataMask,
2065 : nXCount, nYCount,
2066 4 : GDT_Byte, 0, 0 );
2067 : }
2068 :
2069 : /* Compute the resulting overview block */
2070 928 : for(iBand=0;iBand<nBands && eErr == CE_None;iBand++)
2071 : {
2072 : eErr = pfnDownsampleFn(nSrcWidth, nSrcHeight,
2073 : eWrkDataType,
2074 : papaChunk[iBand],
2075 : pabyChunkNoDataMask,
2076 : nChunkXOff, nXCount,
2077 : nChunkYOff, nYCount,
2078 692 : papapoOverviewBands[iBand][iOverview],
2079 : pszResampling,
2080 : pabHasNoData[iBand],
2081 : pafNoDataValue[iBand],
2082 : /*poColorTable*/ NULL,
2083 1384 : eDataType);
2084 : }
2085 : }
2086 :
2087 82 : dfCurPixelCount += (double)nYCount * nSrcWidth;
2088 : }
2089 :
2090 : /* Flush the data to overviews */
2091 136 : for(iBand=0;iBand<nBands;iBand++)
2092 : {
2093 98 : CPLFree(papaChunk[iBand]);
2094 98 : papapoOverviewBands[iBand][iOverview]->FlushCache();
2095 : }
2096 38 : CPLFree(papaChunk);
2097 38 : CPLFree(pabyChunkNoDataMask);
2098 :
2099 : }
2100 :
2101 25 : CPLFree(pabHasNoData);
2102 25 : CPLFree(pafNoDataValue);
2103 :
2104 25 : if (eErr == CE_None)
2105 25 : pfnProgress( 1.0, NULL, pProgressData );
2106 :
2107 25 : return eErr;
2108 : }
2109 :
2110 :
2111 : /************************************************************************/
2112 : /* GDALComputeBandStats() */
2113 : /************************************************************************/
2114 :
2115 : CPLErr CPL_STDCALL
2116 12 : GDALComputeBandStats( GDALRasterBandH hSrcBand,
2117 : int nSampleStep,
2118 : double *pdfMean, double *pdfStdDev,
2119 : GDALProgressFunc pfnProgress,
2120 : void *pProgressData )
2121 :
2122 : {
2123 12 : VALIDATE_POINTER1( hSrcBand, "GDALComputeBandStats", CE_Failure );
2124 :
2125 12 : GDALRasterBand *poSrcBand = (GDALRasterBand *) hSrcBand;
2126 : int iLine, nWidth, nHeight;
2127 12 : GDALDataType eType = poSrcBand->GetRasterDataType();
2128 : GDALDataType eWrkType;
2129 : int bComplex;
2130 : float *pafData;
2131 12 : double dfSum=0.0, dfSum2=0.0;
2132 12 : int nSamples = 0;
2133 :
2134 12 : if( pfnProgress == NULL )
2135 12 : pfnProgress = GDALDummyProgress;
2136 :
2137 12 : nWidth = poSrcBand->GetXSize();
2138 12 : nHeight = poSrcBand->GetYSize();
2139 :
2140 12 : if( nSampleStep >= nHeight || nSampleStep < 1 )
2141 0 : nSampleStep = 1;
2142 :
2143 12 : bComplex = GDALDataTypeIsComplex(eType);
2144 12 : if( bComplex )
2145 : {
2146 0 : pafData = (float *) VSIMalloc(nWidth * 2 * sizeof(float));
2147 0 : eWrkType = GDT_CFloat32;
2148 : }
2149 : else
2150 : {
2151 12 : pafData = (float *) VSIMalloc(nWidth * sizeof(float));
2152 12 : eWrkType = GDT_Float32;
2153 : }
2154 :
2155 12 : if( pafData == NULL )
2156 : {
2157 : CPLError( CE_Failure, CPLE_OutOfMemory,
2158 0 : "GDALComputeBandStats: Out of memory for buffer." );
2159 0 : return CE_Failure;
2160 : }
2161 :
2162 : /* -------------------------------------------------------------------- */
2163 : /* Loop over all sample lines. */
2164 : /* -------------------------------------------------------------------- */
2165 4498 : for( iLine = 0; iLine < nHeight; iLine += nSampleStep )
2166 : {
2167 : int iPixel;
2168 :
2169 4487 : if( !pfnProgress( iLine / (double) nHeight,
2170 : NULL, pProgressData ) )
2171 : {
2172 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
2173 0 : CPLFree( pafData );
2174 0 : return CE_Failure;
2175 : }
2176 :
2177 : CPLErr eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nWidth, 1,
2178 : pafData, nWidth, 1, eWrkType,
2179 4487 : 0, 0 );
2180 4487 : if ( eErr != CE_None )
2181 : {
2182 1 : CPLFree( pafData );
2183 1 : return eErr;
2184 : }
2185 :
2186 6214857 : for( iPixel = 0; iPixel < nWidth; iPixel++ )
2187 : {
2188 : float fValue;
2189 :
2190 6210371 : if( bComplex )
2191 : {
2192 : // Compute the magnitude of the complex value.
2193 :
2194 : fValue = (float)
2195 0 : sqrt(pafData[iPixel*2 ] * pafData[iPixel*2 ]
2196 0 : + pafData[iPixel*2+1] * pafData[iPixel*2+1]);
2197 : }
2198 : else
2199 : {
2200 6210371 : fValue = pafData[iPixel];
2201 : }
2202 :
2203 6210371 : dfSum += fValue;
2204 6210371 : dfSum2 += fValue * fValue;
2205 : }
2206 :
2207 4486 : nSamples += nWidth;
2208 : }
2209 :
2210 11 : if( !pfnProgress( 1.0, NULL, pProgressData ) )
2211 : {
2212 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
2213 0 : CPLFree( pafData );
2214 0 : return CE_Failure;
2215 : }
2216 :
2217 : /* -------------------------------------------------------------------- */
2218 : /* Produce the result values. */
2219 : /* -------------------------------------------------------------------- */
2220 11 : if( pdfMean != NULL )
2221 11 : *pdfMean = dfSum / nSamples;
2222 :
2223 11 : if( pdfStdDev != NULL )
2224 : {
2225 11 : double dfMean = dfSum / nSamples;
2226 :
2227 11 : *pdfStdDev = sqrt((dfSum2 / nSamples) - (dfMean * dfMean));
2228 : }
2229 :
2230 11 : CPLFree( pafData );
2231 :
2232 11 : return CE_None;
2233 : }
2234 :
2235 : /************************************************************************/
2236 : /* GDALOverviewMagnitudeCorrection() */
2237 : /* */
2238 : /* Correct the mean and standard deviation of the overviews of */
2239 : /* the given band to match the base layer approximately. */
2240 : /************************************************************************/
2241 :
2242 : CPLErr
2243 0 : GDALOverviewMagnitudeCorrection( GDALRasterBandH hBaseBand,
2244 : int nOverviewCount,
2245 : GDALRasterBandH *pahOverviews,
2246 : GDALProgressFunc pfnProgress,
2247 : void *pProgressData )
2248 :
2249 : {
2250 0 : VALIDATE_POINTER1( hBaseBand, "GDALOverviewMagnitudeCorrection", CE_Failure );
2251 :
2252 : CPLErr eErr;
2253 : double dfOrigMean, dfOrigStdDev;
2254 :
2255 : /* -------------------------------------------------------------------- */
2256 : /* Compute mean/stddev for source raster. */
2257 : /* -------------------------------------------------------------------- */
2258 : eErr = GDALComputeBandStats( hBaseBand, 2, &dfOrigMean, &dfOrigStdDev,
2259 0 : pfnProgress, pProgressData );
2260 :
2261 0 : if( eErr != CE_None )
2262 0 : return eErr;
2263 :
2264 : /* -------------------------------------------------------------------- */
2265 : /* Loop on overview bands. */
2266 : /* -------------------------------------------------------------------- */
2267 : int iOverview;
2268 :
2269 0 : for( iOverview = 0; iOverview < nOverviewCount; iOverview++ )
2270 : {
2271 0 : GDALRasterBand *poOverview = (GDALRasterBand *)pahOverviews[iOverview];
2272 : double dfOverviewMean, dfOverviewStdDev;
2273 : double dfGain;
2274 :
2275 : eErr = GDALComputeBandStats( pahOverviews[iOverview], 1,
2276 : &dfOverviewMean, &dfOverviewStdDev,
2277 0 : pfnProgress, pProgressData );
2278 :
2279 0 : if( eErr != CE_None )
2280 0 : return eErr;
2281 :
2282 0 : if( dfOrigStdDev < 0.0001 )
2283 0 : dfGain = 1.0;
2284 : else
2285 0 : dfGain = dfOrigStdDev / dfOverviewStdDev;
2286 :
2287 : /* -------------------------------------------------------------------- */
2288 : /* Apply gain and offset. */
2289 : /* -------------------------------------------------------------------- */
2290 0 : GDALDataType eWrkType, eType = poOverview->GetRasterDataType();
2291 : int iLine, nWidth, nHeight, bComplex;
2292 : float *pafData;
2293 :
2294 0 : nWidth = poOverview->GetXSize();
2295 0 : nHeight = poOverview->GetYSize();
2296 :
2297 0 : bComplex = GDALDataTypeIsComplex(eType);
2298 0 : if( bComplex )
2299 : {
2300 0 : pafData = (float *) VSIMalloc2(nWidth, 2 * sizeof(float));
2301 0 : eWrkType = GDT_CFloat32;
2302 : }
2303 : else
2304 : {
2305 0 : pafData = (float *) VSIMalloc2(nWidth, sizeof(float));
2306 0 : eWrkType = GDT_Float32;
2307 : }
2308 :
2309 0 : if( pafData == NULL )
2310 : {
2311 : CPLError( CE_Failure, CPLE_OutOfMemory,
2312 0 : "GDALOverviewMagnitudeCorrection: Out of memory for buffer." );
2313 0 : return CE_Failure;
2314 : }
2315 :
2316 0 : for( iLine = 0; iLine < nHeight; iLine++ )
2317 : {
2318 : int iPixel;
2319 :
2320 0 : if( !pfnProgress( iLine / (double) nHeight,
2321 : NULL, pProgressData ) )
2322 : {
2323 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
2324 0 : CPLFree( pafData );
2325 0 : return CE_Failure;
2326 : }
2327 :
2328 : poOverview->RasterIO( GF_Read, 0, iLine, nWidth, 1,
2329 : pafData, nWidth, 1, eWrkType,
2330 0 : 0, 0 );
2331 :
2332 0 : for( iPixel = 0; iPixel < nWidth; iPixel++ )
2333 : {
2334 0 : if( bComplex )
2335 : {
2336 0 : pafData[iPixel*2] *= (float) dfGain;
2337 0 : pafData[iPixel*2+1] *= (float) dfGain;
2338 : }
2339 : else
2340 : {
2341 0 : pafData[iPixel] = (float)
2342 0 : ((pafData[iPixel]-dfOverviewMean)*dfGain + dfOrigMean);
2343 :
2344 : }
2345 : }
2346 :
2347 : poOverview->RasterIO( GF_Write, 0, iLine, nWidth, 1,
2348 : pafData, nWidth, 1, eWrkType,
2349 0 : 0, 0 );
2350 : }
2351 :
2352 0 : if( !pfnProgress( 1.0, NULL, pProgressData ) )
2353 : {
2354 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
2355 0 : CPLFree( pafData );
2356 0 : return CE_Failure;
2357 : }
2358 :
2359 0 : CPLFree( pafData );
2360 : }
2361 :
2362 0 : return CE_None;
2363 : }
|
