1 : /******************************************************************************
2 : * $Id: overview.cpp 18349 2009-12-19 14:45:41Z 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 18349 2009-12-19 14:45:41Z rouault $");
33 :
34 : typedef enum
35 : {
36 : GRM_Near = 0,
37 : GRM_Average = 1,
38 : GRM_Gauss = 2,
39 : GRM_Mode = 3,
40 : GRM_Cubic = 4,
41 : } GDALResamplingMethod;
42 :
43 : /************************************************************************/
44 : /* GDALDownsampleChunk32R() */
45 : /************************************************************************/
46 :
47 : static CPLErr
48 538 : GDALDownsampleChunk32R( int nSrcWidth, int nSrcHeight,
49 : float * pafChunk,
50 : GByte * pabyChunkNodataMask,
51 : int nChunkXOff, int nChunkXSize,
52 : int nChunkYOff, int nChunkYSize,
53 : GDALRasterBand * poOverview,
54 : const char * pszResampling,
55 : int bHasNoData, float fNoDataValue,
56 : GDALColorTable* poColorTable,
57 : GDALDataType eSrcDataType)
58 :
59 : {
60 : /* -------------------------------------------------------------------- */
61 : /* Check the input parameters. */
62 : /* TODO: I think that instead of passing the pszResampling string */
63 : /* and endless string comparisons this function should take */
64 : /* GDALResamplingMethod flag. The string to flag conversion */
65 : /* should be done on upper level. */
66 : /* -------------------------------------------------------------------- */
67 538 : CPLErr eErr = CE_None;
68 : GDALResamplingMethod eResampling;
69 :
70 538 : if ( EQUALN(pszResampling, "NEAR", 4) )
71 341 : eResampling = GRM_Near;
72 197 : else if ( EQUALN(pszResampling, "AVER", 4) )
73 111 : eResampling = GRM_Average;
74 86 : else if ( EQUALN(pszResampling, "GAUSS", 5) )
75 20 : eResampling = GRM_Gauss;
76 66 : else if ( EQUALN(pszResampling, "MODE", 4) )
77 18 : eResampling = GRM_Mode;
78 48 : else if ( EQUALN(pszResampling, "CUBIC", 6) )
79 48 : eResampling = GRM_Cubic;
80 : else
81 : {
82 : CPLError( CE_Failure, CPLE_AppDefined,
83 : "GDALDownsampleChunk32R: Unsupported resampling method \"%s\".",
84 0 : pszResampling );
85 0 : return CE_Failure;
86 : }
87 :
88 : /* -------------------------------------------------------------------- */
89 : /* Create the filter kernel and allocate scanline buffer. */
90 : /* -------------------------------------------------------------------- */
91 : int nDstXOff, nDstXOff2, nDstYOff, nDstYOff2, nOXSize, nOYSize;
92 : float *pafDstScanline;
93 538 : int nGaussMatrixDim = 3;
94 : const int *panGaussMatrix;
95 : static const int anGaussMatrix3x3[] ={
96 : 1,2,1,
97 : 2,4,2,
98 : 1,2,1
99 : };
100 : static const int anGaussMatrix5x5[] = {
101 : 1,4,6,4,1,
102 : 4,16,24,16,4,
103 : 6,24,36,24,6,
104 : 4,16,24,16,4,
105 : 1,4,6,4,1};
106 : static const int anGaussMatrix7x7[] = {
107 : 1,6,15,20,15,6,1,
108 : 6,36,90,120,90,36,6,
109 : 15,90,225,300,225,90,15,
110 : 20,120,300,400,300,120,20,
111 : 15,90,225,300,225,90,15,
112 : 6,36,90,120,90,36,6,
113 : 1,6,15,20,15,6,1};
114 :
115 538 : nOXSize = poOverview->GetXSize();
116 538 : nOYSize = poOverview->GetYSize();
117 538 : int nResYFactor = (int) (0.5 + (double)nSrcHeight/(double)nOYSize);
118 :
119 : // matrix for gauss filter
120 538 : if(nResYFactor <= 2 )
121 : {
122 326 : panGaussMatrix = anGaussMatrix3x3;
123 326 : nGaussMatrixDim=3;
124 : }
125 212 : else if (nResYFactor <= 4)
126 : {
127 186 : panGaussMatrix = anGaussMatrix5x5;
128 186 : nGaussMatrixDim=5;
129 : }
130 : else
131 : {
132 26 : panGaussMatrix = anGaussMatrix7x7;
133 26 : nGaussMatrixDim=7;
134 : }
135 :
136 : /* -------------------------------------------------------------------- */
137 : /* Figure out the column to start writing to, and the first column */
138 : /* to not write to. */
139 : /* -------------------------------------------------------------------- */
140 538 : nDstXOff = (int) (0.5 + (nChunkXOff/(double)nSrcWidth) * nOXSize);
141 : nDstXOff2 = (int)
142 538 : (0.5 + ((nChunkXOff+nChunkXSize)/(double)nSrcWidth) * nOXSize);
143 :
144 538 : if( nChunkXOff + nChunkXSize == nSrcWidth )
145 538 : nDstXOff2 = nOXSize;
146 :
147 :
148 538 : pafDstScanline = (float *) VSIMalloc((nDstXOff2 - nDstXOff) * sizeof(float));
149 538 : if( pafDstScanline == NULL )
150 : {
151 : CPLError( CE_Failure, CPLE_OutOfMemory,
152 0 : "GDALDownsampleChunk32R: Out of memory for line buffer." );
153 0 : return CE_Failure;
154 : }
155 :
156 : /* -------------------------------------------------------------------- */
157 : /* Figure out the line to start writing to, and the first line */
158 : /* to not write to. In theory this approach should ensure that */
159 : /* every output line will be written if all input chunks are */
160 : /* processed. */
161 : /* -------------------------------------------------------------------- */
162 538 : nDstYOff = (int) (0.5 + (nChunkYOff/(double)nSrcHeight) * nOYSize);
163 : nDstYOff2 = (int)
164 538 : (0.5 + ((nChunkYOff+nChunkYSize)/(double)nSrcHeight) * nOYSize);
165 :
166 538 : if( nChunkYOff + nChunkYSize == nSrcHeight )
167 232 : nDstYOff2 = nOYSize;
168 :
169 :
170 538 : int nEntryCount = 0;
171 538 : GDALColorEntry* aEntries = NULL;
172 538 : if (poColorTable)
173 : {
174 : int i;
175 6 : nEntryCount = poColorTable->GetColorEntryCount();
176 6 : aEntries = (GDALColorEntry* )CPLMalloc(sizeof(GDALColorEntry) * nEntryCount);
177 1542 : for(i=0;i<nEntryCount;i++)
178 : {
179 1536 : poColorTable->GetColorEntryAsRGB(i, &aEntries[i]);
180 : }
181 : }
182 :
183 538 : int nMaxNumPx = 0;
184 538 : float* pafVals = NULL;
185 538 : int* panSums = NULL;
186 :
187 : /* ==================================================================== */
188 : /* Loop over destination scanlines. */
189 : /* ==================================================================== */
190 11178 : for( int iDstLine = nDstYOff; iDstLine < nDstYOff2 && eErr == CE_None; iDstLine++ )
191 : {
192 : float *pafSrcScanline;
193 : GByte *pabySrcScanlineNodataMask;
194 10640 : int nSrcYOff, nSrcYOff2 = 0, iDstPixel;
195 :
196 10640 : if (eResampling == GRM_Gauss)
197 : {
198 260 : nSrcYOff = (int) (0.5 + (iDstLine/(double)nOYSize) * nSrcHeight);
199 260 : nSrcYOff2 = (int) (0.5 + ((iDstLine+1)/(double)nOYSize) * nSrcHeight) + 1;
200 :
201 260 : if( nSrcYOff < nChunkYOff )
202 : {
203 0 : nSrcYOff = nChunkYOff;
204 0 : nSrcYOff2++;
205 : }
206 :
207 260 : int iSizeY = nSrcYOff2 - nSrcYOff;
208 260 : nSrcYOff = nSrcYOff + iSizeY/2 - nGaussMatrixDim/2;
209 260 : nSrcYOff2 = nSrcYOff + nGaussMatrixDim;
210 260 : if(nSrcYOff < 0)
211 0 : nSrcYOff = 0;
212 : }
213 10380 : else if( eResampling == GRM_Cubic )
214 : {
215 840 : nSrcYOff = (int) floor(((iDstLine+0.5)/(double)nOYSize) * nSrcHeight - 0.5)-1;
216 840 : nSrcYOff2 = nSrcYOff + 4;
217 840 : if(nSrcYOff < 0)
218 8 : nSrcYOff = 0;
219 840 : if(nSrcYOff < nChunkYOff)
220 16 : nSrcYOff = nChunkYOff;
221 : }
222 : else
223 : {
224 9540 : nSrcYOff = (int) (0.5 + (iDstLine/(double)nOYSize) * nSrcHeight);
225 9540 : if ( nSrcYOff < nChunkYOff )
226 10 : nSrcYOff = nChunkYOff;
227 :
228 : nSrcYOff2 =
229 9540 : (int) (0.5 + ((iDstLine+1)/(double)nOYSize) * nSrcHeight);
230 : }
231 :
232 10640 : if( nSrcYOff2 > nSrcHeight || iDstLine == nOYSize-1 )
233 232 : nSrcYOff2 = nSrcHeight;
234 10640 : if( nSrcYOff2 > nChunkYOff + nChunkYSize)
235 48 : nSrcYOff2 = nChunkYOff + nChunkYSize;
236 :
237 10640 : pafSrcScanline = pafChunk + ((nSrcYOff-nChunkYOff) * nChunkXSize);
238 10640 : if (pabyChunkNodataMask != NULL)
239 950 : pabySrcScanlineNodataMask = pabyChunkNodataMask + ((nSrcYOff-nChunkYOff) * nChunkXSize);
240 : else
241 9690 : pabySrcScanlineNodataMask = NULL;
242 :
243 : /* -------------------------------------------------------------------- */
244 : /* Loop over destination pixels */
245 : /* -------------------------------------------------------------------- */
246 3458484 : for( iDstPixel = nDstXOff; iDstPixel < nDstXOff2; iDstPixel++ )
247 : {
248 : int nSrcXOff, nSrcXOff2;
249 :
250 3447844 : if (eResampling == GRM_Gauss)
251 : {
252 4700 : nSrcXOff = (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
253 4700 : nSrcXOff2 = (int)(0.5 + ((iDstPixel+1)/(double)nOXSize) * nSrcWidth) + 1;
254 :
255 4700 : int iSizeX = nSrcXOff2 - nSrcXOff;
256 4700 : nSrcXOff = nSrcXOff + iSizeX/2 - nGaussMatrixDim/2;
257 4700 : nSrcXOff2 = nSrcXOff + nGaussMatrixDim;
258 4700 : if(nSrcXOff < 0)
259 0 : nSrcXOff = 0;
260 : }
261 3443144 : else if( eResampling == GRM_Cubic )
262 : {
263 56520 : nSrcXOff = (int) floor(((iDstPixel+0.5)/(double)nOXSize) * nSrcWidth - 0.5)-1;
264 56520 : nSrcXOff2 = nSrcXOff + 4;
265 :
266 56520 : if(nSrcXOff < 0)
267 600 : nSrcXOff = 0;
268 : }
269 : else
270 : {
271 : nSrcXOff =
272 3386624 : (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
273 3386624 : if ( nSrcXOff < nChunkXOff )
274 0 : nSrcXOff = nChunkXOff;
275 : nSrcXOff2 = (int)
276 3386624 : (0.5 + ((iDstPixel+1)/(double)nOXSize) * nSrcWidth);
277 : }
278 :
279 3447844 : if( nSrcXOff2 > nSrcWidth || iDstPixel == nOXSize-1 )
280 10640 : nSrcXOff2 = nSrcWidth;
281 3447844 : if( nSrcXOff2 > nChunkXOff + nChunkXSize )
282 0 : nSrcXOff2 = nChunkXOff + nChunkXSize;
283 :
284 3447844 : if ( eResampling == GRM_Near )
285 : {
286 3332224 : pafDstScanline[iDstPixel - nDstXOff] = pafSrcScanline[nSrcXOff - nChunkXOff];
287 : }
288 115620 : else if( eResampling == GRM_Cubic )
289 : {
290 : // If we do not seem to have our full 4x4 window just
291 : // do nearest resampling.
292 62040 : if( nSrcXOff2 - nSrcXOff != 4 || nSrcYOff2 - nSrcYOff != 4 )
293 : {
294 5520 : int nLSrcYOff = (int) (0.5+(iDstLine/(double)nOYSize) * nSrcHeight);
295 5520 : int nLSrcXOff = (int) (0.5+(iDstPixel/(double)nOXSize) * nSrcWidth);
296 :
297 5520 : if( nLSrcYOff < nChunkYOff )
298 0 : nLSrcYOff = nChunkYOff;
299 5520 : if( nLSrcYOff > nChunkYOff + nChunkYSize - 1 )
300 0 : nLSrcYOff = nChunkYOff + nChunkYSize - 1;
301 :
302 5520 : pafDstScanline[iDstPixel - nDstXOff] =
303 : pafChunk[(nLSrcYOff-nChunkYOff) * nChunkXSize
304 5520 : + (nLSrcXOff - nChunkXOff)];
305 : }
306 : else
307 : {
308 : #define CubicConvolution(distance1,distance2,distance3,f0,f1,f2,f3) \
309 : ( ( -f0 + f1 - f2 + f3) * distance3 \
310 : + (2.0*(f0 - f1) + f2 - f3) * distance2 \
311 : + ( -f0 + f2 ) * distance1 \
312 : + f1 )
313 :
314 : int ic;
315 : double adfRowResults[4];
316 51000 : double dfSrcX = (((iDstPixel+0.5)/(double)nOXSize) * nSrcWidth);
317 51000 : double dfDeltaX = dfSrcX - 0.5 - (nSrcXOff+1);
318 51000 : double dfDeltaX2 = dfDeltaX * dfDeltaX;
319 51000 : double dfDeltaX3 = dfDeltaX2 * dfDeltaX;
320 :
321 255000 : for ( ic = 0; ic < 4; ic++ )
322 : {
323 : float *pafSrcRow = pafSrcScanline +
324 204000 : nSrcXOff-nChunkXOff+(nSrcYOff+ic-nSrcYOff)*nChunkXSize;
325 :
326 204000 : adfRowResults[ic] =
327 1632000 : CubicConvolution(dfDeltaX, dfDeltaX2, dfDeltaX3,
328 : pafSrcRow[0],
329 : pafSrcRow[1],
330 : pafSrcRow[2],
331 1632000 : pafSrcRow[3] );
332 : }
333 :
334 51000 : double dfSrcY = (((iDstLine+0.5)/(double)nOYSize) * nSrcHeight);
335 51000 : double dfDeltaY = dfSrcY - 0.5 - (nSrcYOff+1);
336 51000 : double dfDeltaY2 = dfDeltaY * dfDeltaY;
337 51000 : double dfDeltaY3 = dfDeltaY2 * dfDeltaY;
338 :
339 51000 : pafDstScanline[iDstPixel - nDstXOff] = (float)
340 408000 : CubicConvolution(dfDeltaY, dfDeltaY2, dfDeltaY3,
341 : adfRowResults[0],
342 : adfRowResults[1],
343 : adfRowResults[2],
344 408000 : adfRowResults[3] );
345 : }
346 : }
347 59100 : else if ( eResampling == GRM_Average )
348 : {
349 106200 : if (poColorTable == NULL
350 : || EQUALN(pszResampling,"AVERAGE_BIT2GRAYSCALE",13))
351 : {
352 : double val;
353 53000 : double dfTotal = 0.0;
354 53000 : int nCount = 0, iX, iY;
355 :
356 166300 : for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
357 : {
358 369504 : for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
359 : {
360 256204 : val = pafSrcScanline[iX-nChunkXOff+(iY-nSrcYOff)*nChunkXSize];
361 408204 : if (pabySrcScanlineNodataMask == NULL ||
362 152000 : pabySrcScanlineNodataMask[iX-nChunkXOff+(iY-nSrcYOff)*nChunkXSize])
363 : {
364 236904 : dfTotal += val;
365 236904 : nCount++;
366 : }
367 : }
368 : }
369 :
370 57800 : if (bHasNoData && nCount == 0)
371 : {
372 4800 : pafDstScanline[iDstPixel - nDstXOff] = fNoDataValue;
373 : }
374 : else
375 : {
376 48200 : if( nCount == 0 )
377 0 : pafDstScanline[iDstPixel - nDstXOff] = 0.0;
378 : else
379 48200 : pafDstScanline[iDstPixel - nDstXOff] = (float) (dfTotal / nCount);
380 : }
381 : }
382 : else
383 : {
384 : double val;
385 200 : int nTotalR = 0, nTotalG = 0, nTotalB = 0;
386 200 : int nCount = 0, iX, iY;
387 :
388 600 : for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
389 : {
390 1200 : for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
391 : {
392 800 : val = pafSrcScanline[iX-nChunkXOff+(iY-nSrcYOff)*nChunkXSize];
393 800 : if (bHasNoData == FALSE || val != fNoDataValue)
394 : {
395 800 : int nVal = (int)val;
396 800 : if (nVal >= 0 && nVal < nEntryCount)
397 : {
398 800 : nTotalR += aEntries[nVal].c1;
399 800 : nTotalG += aEntries[nVal].c2;
400 800 : nTotalB += aEntries[nVal].c3;
401 800 : nCount++;
402 : }
403 : }
404 : }
405 : }
406 :
407 200 : if (bHasNoData && nCount == 0)
408 : {
409 0 : pafDstScanline[iDstPixel - nDstXOff] = fNoDataValue;
410 : }
411 : else
412 : {
413 : CPLAssert( nCount > 0 );
414 200 : if( nCount == 0 )
415 0 : pafDstScanline[iDstPixel - nDstXOff] = 0.0;
416 : else
417 : {
418 200 : int nR = nTotalR / nCount, nG = nTotalG / nCount, nB = nTotalB / nCount;
419 : int i;
420 200 : double dfMinDist = 0;
421 200 : int iBestEntry = 0;
422 51400 : for(i=0;i<nEntryCount;i++)
423 : {
424 102400 : double dfDist = (nR - aEntries[i].c1) * (nR - aEntries[i].c1) +
425 102400 : (nG - aEntries[i].c2) * (nG - aEntries[i].c2) +
426 204800 : (nB - aEntries[i].c3) * (nB - aEntries[i].c3);
427 51200 : if (i == 0 || dfDist < dfMinDist)
428 : {
429 400 : dfMinDist = dfDist;
430 400 : iBestEntry = i;
431 : }
432 : }
433 200 : pafDstScanline[iDstPixel - nDstXOff] = iBestEntry;
434 : }
435 : }
436 : }
437 : }
438 5900 : else if ( eResampling == GRM_Mode )
439 : {
440 1950 : if (eSrcDataType != GDT_Byte || nEntryCount > 256)
441 : {
442 : /* I'm not sure how much sense it makes to run a majority
443 : filter on floating point data, but here it is for the sake
444 : of compatability. It won't look right on RGB images by the
445 : nature of the filter. */
446 750 : int nNumPx = (nSrcYOff2-nSrcYOff)*(nSrcXOff2-nSrcXOff);
447 750 : int iMaxInd = 0, iMaxVal = -1, iY, iX;
448 :
449 750 : if (nNumPx > nMaxNumPx)
450 : {
451 12 : pafVals = (float*) CPLRealloc(pafVals, nNumPx * sizeof(float));
452 12 : panSums = (int*) CPLRealloc(panSums, nNumPx * sizeof(int));
453 12 : nMaxNumPx = nNumPx;
454 : }
455 :
456 2550 : for( iY = nSrcYOff; iY < nSrcYOff2; ++iY )
457 : {
458 1800 : int iTotYOff = (iY-nSrcYOff)*nChunkXSize-nChunkXOff;
459 6600 : for( iX = nSrcXOff; iX < nSrcXOff2; ++iX )
460 : {
461 4800 : if (pabySrcScanlineNodataMask == NULL ||
462 0 : pabySrcScanlineNodataMask[iX+iTotYOff])
463 : {
464 4800 : float fVal = pafSrcScanline[iX+iTotYOff];
465 : int i;
466 :
467 : //Check array for existing entry
468 23148 : for( i = 0; i < iMaxInd; ++i )
469 23436 : if( pafVals[i] == fVal
470 4668 : && ++panSums[i] > panSums[iMaxVal] )
471 : {
472 420 : iMaxVal = i;
473 420 : break;
474 : }
475 :
476 : //Add to arr if entry not already there
477 4800 : if( i == iMaxInd )
478 : {
479 4380 : pafVals[iMaxInd] = fVal;
480 4380 : panSums[iMaxInd] = 1;
481 :
482 4380 : if( iMaxVal < 0 )
483 750 : iMaxVal = iMaxInd;
484 :
485 4380 : ++iMaxInd;
486 : }
487 : }
488 : }
489 : }
490 :
491 750 : if( iMaxVal == -1 )
492 0 : pafDstScanline[iDstPixel - nDstXOff] = fNoDataValue;
493 : else
494 750 : pafDstScanline[iDstPixel - nDstXOff] = pafVals[iMaxVal];
495 : }
496 : else /* if (eSrcDataType == GDT_Byte && nEntryCount < 256) */
497 : {
498 : /* So we go here for a paletted or non-paletted byte band */
499 : /* The input values are then between 0 and 255 */
500 450 : int anVals[256], nMaxVal = 0, iMaxInd = -1, iY, iX;
501 :
502 450 : memset(anVals, 0, 256*sizeof(int));
503 :
504 1450 : for( iY = nSrcYOff; iY < nSrcYOff2; ++iY )
505 : {
506 1000 : int iTotYOff = (iY-nSrcYOff)*nChunkXSize-nChunkXOff;
507 3400 : for( iX = nSrcXOff; iX < nSrcXOff2; ++iX )
508 : {
509 2400 : float val = pafSrcScanline[iX+iTotYOff];
510 2400 : if (bHasNoData == FALSE || val != fNoDataValue)
511 : {
512 2400 : int nVal = (int) val;
513 2400 : if ( ++anVals[nVal] > nMaxVal)
514 : {
515 : //Sum the density
516 : //Is it the most common value so far?
517 948 : iMaxInd = nVal;
518 948 : nMaxVal = anVals[nVal];
519 : }
520 : }
521 : }
522 : }
523 :
524 450 : if( iMaxInd == -1 )
525 0 : pafDstScanline[iDstPixel - nDstXOff] = fNoDataValue;
526 : else
527 450 : pafDstScanline[iDstPixel - nDstXOff] = iMaxInd;
528 : }
529 : }
530 4700 : else if ( eResampling == GRM_Gauss )
531 : {
532 4700 : if (poColorTable == NULL)
533 : {
534 4500 : double dfTotal = 0.0, val;
535 4500 : int nCount = 0, iX, iY;
536 4500 : int i = 0,j = 0;
537 4500 : const int *panLineWeight = panGaussMatrix;
538 :
539 17760 : for( j=0, iY = nSrcYOff; iY < nSrcYOff2;
540 : iY++, j++, panLineWeight += nGaussMatrixDim )
541 : {
542 52338 : for( i=0, iX = nSrcXOff; iX < nSrcXOff2; iX++,++i )
543 : {
544 39078 : val = pafSrcScanline[iX-nChunkXOff+(iY-nSrcYOff)*nChunkXSize];
545 71934 : if (pabySrcScanlineNodataMask == NULL ||
546 32856 : pabySrcScanlineNodataMask[iX-nChunkXOff+(iY-nSrcYOff)*nChunkXSize])
547 : {
548 18372 : int nWeight = panLineWeight[i];
549 18372 : dfTotal += val * nWeight;
550 18372 : nCount += nWeight;
551 : }
552 : }
553 : }
554 :
555 6714 : if (bHasNoData && nCount == 0)
556 : {
557 2214 : pafDstScanline[iDstPixel - nDstXOff] = fNoDataValue;
558 : }
559 : else
560 : {
561 2286 : if( nCount == 0 )
562 0 : pafDstScanline[iDstPixel - nDstXOff] = 0.0;
563 : else
564 2286 : pafDstScanline[iDstPixel - nDstXOff] = (float) (dfTotal / nCount);
565 : }
566 : }
567 : else
568 : {
569 : double val;
570 200 : int nTotalR = 0, nTotalG = 0, nTotalB = 0;
571 200 : int nTotalWeight = 0, iX, iY;
572 200 : int i = 0,j = 0;
573 200 : const int *panLineWeight = panGaussMatrix;
574 :
575 780 : for( j=0, iY = nSrcYOff; iY < nSrcYOff2;
576 : iY++, j++, panLineWeight += nGaussMatrixDim )
577 : {
578 2262 : for( i=0, iX = nSrcXOff; iX < nSrcXOff2; iX++,++i )
579 : {
580 1682 : val = pafSrcScanline[iX-nChunkXOff+(iY-nSrcYOff)*nChunkXSize];
581 1682 : if (bHasNoData == FALSE || val != fNoDataValue)
582 : {
583 1682 : int nVal = (int)val;
584 1682 : if (nVal >= 0 && nVal < nEntryCount)
585 : {
586 1682 : int nWeight = panLineWeight[i];
587 1682 : nTotalR += aEntries[nVal].c1 * nWeight;
588 1682 : nTotalG += aEntries[nVal].c2 * nWeight;
589 1682 : nTotalB += aEntries[nVal].c3 * nWeight;
590 1682 : nTotalWeight += nWeight;
591 : }
592 : }
593 : }
594 : }
595 :
596 200 : if (bHasNoData && nTotalWeight == 0)
597 : {
598 0 : pafDstScanline[iDstPixel - nDstXOff] = fNoDataValue;
599 : }
600 : else
601 : {
602 : CPLAssert( nTotalWeight > 0 );
603 200 : if( nTotalWeight == 0 )
604 0 : pafDstScanline[iDstPixel - nDstXOff] = 0.0;
605 : else
606 : {
607 200 : int nR = nTotalR / nTotalWeight, nG = nTotalG / nTotalWeight, nB = nTotalB / nTotalWeight;
608 : int i;
609 200 : double dfMinDist = 0;
610 200 : int iBestEntry = 0;
611 51400 : for(i=0;i<nEntryCount;i++)
612 : {
613 102400 : double dfDist = (nR - aEntries[i].c1) * (nR - aEntries[i].c1) +
614 102400 : (nG - aEntries[i].c2) * (nG - aEntries[i].c2) +
615 204800 : (nB - aEntries[i].c3) * (nB - aEntries[i].c3);
616 51200 : if (i == 0 || dfDist < dfMinDist)
617 : {
618 400 : dfMinDist = dfDist;
619 400 : iBestEntry = i;
620 : }
621 : }
622 200 : pafDstScanline[iDstPixel - nDstXOff] = iBestEntry;
623 : }
624 : }
625 : }
626 : } // end of gauss
627 : }
628 :
629 : eErr = poOverview->RasterIO( GF_Write, nDstXOff, iDstLine, nDstXOff2 - nDstXOff, 1,
630 : pafDstScanline, nDstXOff2 - nDstXOff, 1, GDT_Float32,
631 10640 : 0, 0 );
632 : }
633 :
634 538 : CPLFree( pafDstScanline );
635 538 : CPLFree( aEntries );
636 538 : CPLFree( pafVals );
637 538 : CPLFree( panSums );
638 :
639 538 : return eErr;
640 : }
641 :
642 : /************************************************************************/
643 : /* GDALDownsampleChunkC32R() */
644 : /************************************************************************/
645 :
646 : static CPLErr
647 0 : GDALDownsampleChunkC32R( int nSrcWidth, int nSrcHeight,
648 : float * pafChunk, int nChunkYOff, int nChunkYSize,
649 : GDALRasterBand * poOverview,
650 : const char * pszResampling )
651 :
652 : {
653 : int nDstYOff, nDstYOff2, nOXSize, nOYSize;
654 : float *pafDstScanline;
655 0 : CPLErr eErr = CE_None;
656 :
657 0 : nOXSize = poOverview->GetXSize();
658 0 : nOYSize = poOverview->GetYSize();
659 :
660 0 : pafDstScanline = (float *) VSIMalloc(nOXSize * sizeof(float) * 2);
661 0 : if( pafDstScanline == NULL )
662 : {
663 : CPLError( CE_Failure, CPLE_OutOfMemory,
664 0 : "GDALDownsampleChunkC32R: Out of memory for line buffer." );
665 0 : return CE_Failure;
666 : }
667 :
668 : /* -------------------------------------------------------------------- */
669 : /* Figure out the line to start writing to, and the first line */
670 : /* to not write to. In theory this approach should ensure that */
671 : /* every output line will be written if all input chunks are */
672 : /* processed. */
673 : /* -------------------------------------------------------------------- */
674 0 : nDstYOff = (int) (0.5 + (nChunkYOff/(double)nSrcHeight) * nOYSize);
675 : nDstYOff2 = (int)
676 0 : (0.5 + ((nChunkYOff+nChunkYSize)/(double)nSrcHeight) * nOYSize);
677 :
678 0 : if( nChunkYOff + nChunkYSize == nSrcHeight )
679 0 : nDstYOff2 = nOYSize;
680 :
681 : /* ==================================================================== */
682 : /* Loop over destination scanlines. */
683 : /* ==================================================================== */
684 0 : for( int iDstLine = nDstYOff; iDstLine < nDstYOff2 && eErr == CE_None; iDstLine++ )
685 : {
686 : float *pafSrcScanline;
687 : int nSrcYOff, nSrcYOff2, iDstPixel;
688 :
689 0 : nSrcYOff = (int) (0.5 + (iDstLine/(double)nOYSize) * nSrcHeight);
690 0 : if( nSrcYOff < nChunkYOff )
691 0 : nSrcYOff = nChunkYOff;
692 :
693 0 : nSrcYOff2 = (int) (0.5 + ((iDstLine+1)/(double)nOYSize) * nSrcHeight);
694 0 : if( nSrcYOff2 > nSrcHeight || iDstLine == nOYSize-1 )
695 0 : nSrcYOff2 = nSrcHeight;
696 0 : if( nSrcYOff2 > nChunkYOff + nChunkYSize )
697 0 : nSrcYOff2 = nChunkYOff + nChunkYSize;
698 :
699 0 : pafSrcScanline = pafChunk + ((nSrcYOff-nChunkYOff) * nSrcWidth) * 2;
700 :
701 : /* -------------------------------------------------------------------- */
702 : /* Loop over destination pixels */
703 : /* -------------------------------------------------------------------- */
704 0 : for( iDstPixel = 0; iDstPixel < nOXSize; iDstPixel++ )
705 : {
706 : int nSrcXOff, nSrcXOff2;
707 :
708 0 : nSrcXOff = (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
709 : nSrcXOff2 = (int)
710 0 : (0.5 + ((iDstPixel+1)/(double)nOXSize) * nSrcWidth);
711 0 : if( nSrcXOff2 > nSrcWidth )
712 0 : nSrcXOff2 = nSrcWidth;
713 :
714 0 : if( EQUALN(pszResampling,"NEAR",4) )
715 : {
716 0 : pafDstScanline[iDstPixel*2] = pafSrcScanline[nSrcXOff*2];
717 0 : pafDstScanline[iDstPixel*2+1] = pafSrcScanline[nSrcXOff*2+1];
718 : }
719 0 : else if( EQUAL(pszResampling,"AVERAGE_MAGPHASE") )
720 : {
721 0 : double dfTotalR = 0.0, dfTotalI = 0.0, dfTotalM = 0.0;
722 0 : int nCount = 0, iX, iY;
723 :
724 0 : for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
725 : {
726 0 : for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
727 : {
728 : double dfR, dfI;
729 :
730 0 : dfR = pafSrcScanline[iX*2+(iY-nSrcYOff)*nSrcWidth*2];
731 0 : dfI = pafSrcScanline[iX*2+(iY-nSrcYOff)*nSrcWidth*2+1];
732 0 : dfTotalR += dfR;
733 0 : dfTotalI += dfI;
734 0 : dfTotalM += sqrt( dfR*dfR + dfI*dfI );
735 0 : nCount++;
736 : }
737 : }
738 :
739 : CPLAssert( nCount > 0 );
740 0 : if( nCount == 0 )
741 : {
742 0 : pafDstScanline[iDstPixel*2] = 0.0;
743 0 : pafDstScanline[iDstPixel*2+1] = 0.0;
744 : }
745 : else
746 : {
747 0 : double dfM, dfDesiredM, dfRatio=1.0;
748 :
749 0 : pafDstScanline[iDstPixel*2 ] = (float) (dfTotalR/nCount);
750 0 : pafDstScanline[iDstPixel*2+1] = (float) (dfTotalI/nCount);
751 :
752 0 : dfM = sqrt(pafDstScanline[iDstPixel*2 ]*pafDstScanline[iDstPixel*2 ]
753 0 : + pafDstScanline[iDstPixel*2+1]*pafDstScanline[iDstPixel*2+1]);
754 0 : dfDesiredM = dfTotalM / nCount;
755 0 : if( dfM != 0.0 )
756 0 : dfRatio = dfDesiredM / dfM;
757 :
758 0 : pafDstScanline[iDstPixel*2 ] *= (float) dfRatio;
759 0 : pafDstScanline[iDstPixel*2+1] *= (float) dfRatio;
760 : }
761 : }
762 0 : else if( EQUALN(pszResampling,"AVER",4) )
763 : {
764 0 : double dfTotalR = 0.0, dfTotalI = 0.0;
765 0 : int nCount = 0, iX, iY;
766 :
767 0 : for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
768 : {
769 0 : for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
770 : {
771 0 : dfTotalR += pafSrcScanline[iX*2+(iY-nSrcYOff)*nSrcWidth*2];
772 0 : dfTotalI += pafSrcScanline[iX*2+(iY-nSrcYOff)*nSrcWidth*2+1];
773 0 : nCount++;
774 : }
775 : }
776 :
777 : CPLAssert( nCount > 0 );
778 0 : if( nCount == 0 )
779 : {
780 0 : pafDstScanline[iDstPixel*2] = 0.0;
781 0 : pafDstScanline[iDstPixel*2+1] = 0.0;
782 : }
783 : else
784 : {
785 0 : pafDstScanline[iDstPixel*2 ] = (float) (dfTotalR/nCount);
786 0 : pafDstScanline[iDstPixel*2+1] = (float) (dfTotalI/nCount);
787 : }
788 : }
789 : }
790 :
791 : eErr = poOverview->RasterIO( GF_Write, 0, iDstLine, nOXSize, 1,
792 : pafDstScanline, nOXSize, 1, GDT_CFloat32,
793 0 : 0, 0 );
794 : }
795 :
796 0 : CPLFree( pafDstScanline );
797 :
798 0 : return eErr;
799 : }
800 :
801 : /************************************************************************/
802 : /* GDALRegenerateCascadingOverviews() */
803 : /* */
804 : /* Generate a list of overviews in order from largest to */
805 : /* smallest, computing each from the next larger. */
806 : /************************************************************************/
807 :
808 : static CPLErr
809 12 : GDALRegenerateCascadingOverviews(
810 : GDALRasterBand *poSrcBand, int nOverviews, GDALRasterBand **papoOvrBands,
811 : const char * pszResampling,
812 : GDALProgressFunc pfnProgress, void * pProgressData )
813 :
814 : {
815 : /* -------------------------------------------------------------------- */
816 : /* First, we must put the overviews in order from largest to */
817 : /* smallest. */
818 : /* -------------------------------------------------------------------- */
819 : int i, j;
820 :
821 24 : for( i = 0; i < nOverviews-1; i++ )
822 : {
823 24 : for( j = 0; j < nOverviews - i - 1; j++ )
824 : {
825 :
826 48 : if( papoOvrBands[j]->GetXSize()
827 12 : * (float) papoOvrBands[j]->GetYSize() <
828 12 : papoOvrBands[j+1]->GetXSize()
829 12 : * (float) papoOvrBands[j+1]->GetYSize() )
830 : {
831 : GDALRasterBand * poTempBand;
832 :
833 0 : poTempBand = papoOvrBands[j];
834 0 : papoOvrBands[j] = papoOvrBands[j+1];
835 0 : papoOvrBands[j+1] = poTempBand;
836 : }
837 : }
838 : }
839 :
840 : /* -------------------------------------------------------------------- */
841 : /* Count total pixels so we can prepare appropriate scaled */
842 : /* progress functions. */
843 : /* -------------------------------------------------------------------- */
844 12 : double dfTotalPixels = 0.0;
845 :
846 36 : for( i = 0; i < nOverviews; i++ )
847 : {
848 24 : dfTotalPixels += papoOvrBands[i]->GetXSize()
849 48 : * (double) papoOvrBands[i]->GetYSize();
850 : }
851 :
852 : /* -------------------------------------------------------------------- */
853 : /* Generate all the bands. */
854 : /* -------------------------------------------------------------------- */
855 12 : double dfPixelsProcessed = 0.0;
856 :
857 36 : for( i = 0; i < nOverviews; i++ )
858 : {
859 : void *pScaledProgressData;
860 : double dfPixels;
861 : GDALRasterBand *poBaseBand;
862 : CPLErr eErr;
863 :
864 24 : if( i == 0 )
865 12 : poBaseBand = poSrcBand;
866 : else
867 12 : poBaseBand = papoOvrBands[i-1];
868 :
869 24 : dfPixels = papoOvrBands[i]->GetXSize()
870 48 : * (double) papoOvrBands[i]->GetYSize();
871 :
872 : pScaledProgressData = GDALCreateScaledProgress(
873 : dfPixelsProcessed / dfTotalPixels,
874 : (dfPixelsProcessed + dfPixels) / dfTotalPixels,
875 24 : pfnProgress, pProgressData );
876 :
877 : eErr = GDALRegenerateOverviews( (GDALRasterBandH) poBaseBand,
878 : 1, (GDALRasterBandH *) papoOvrBands+i,
879 : pszResampling,
880 : GDALScaledProgress,
881 24 : pScaledProgressData );
882 24 : GDALDestroyScaledProgress( pScaledProgressData );
883 :
884 24 : if( eErr != CE_None )
885 0 : return eErr;
886 :
887 24 : dfPixelsProcessed += dfPixels;
888 :
889 : /* we only do the bit2grayscale promotion on the base band */
890 24 : if( EQUALN(pszResampling,"AVERAGE_BIT2GRAYSCALE",13) )
891 4 : pszResampling = "AVERAGE";
892 : }
893 :
894 12 : return CE_None;
895 : }
896 :
897 : /************************************************************************/
898 : /* GDALRegenerateOverviews() */
899 : /************************************************************************/
900 :
901 : /**
902 : * \brief Generate downsampled overviews.
903 : *
904 : * This function will generate one or more overview images from a base
905 : * image using the requested downsampling algorithm. It's primary use
906 : * is for generating overviews via GDALDataset::BuildOverviews(), but it
907 : * can also be used to generate downsampled images in one file from another
908 : * outside the overview architecture.
909 : *
910 : * The output bands need to exist in advance.
911 : *
912 : * The full set of resampling algorithms is documented in
913 : * GDALDataset::BuildOverviews().
914 : *
915 : * This function will honour properly NODATA_VALUES tuples (special dataset metadata) so
916 : * that only a given RGB triplet (in case of a RGB image) will be considered as the
917 : * nodata value and not each value of the triplet independantly per band.
918 : *
919 : * @param hSrcBand the source (base level) band.
920 : * @param nOverviewCount the number of downsampled bands being generated.
921 : * @param pahOvrBands the list of downsampled bands to be generated.
922 : * @param pszResampling Resampling algorithm (eg. "AVERAGE").
923 : * @param pfnProgress progress report function.
924 : * @param pProgressData progress function callback data.
925 : * @return CE_None on success or CE_Failure on failure.
926 : */
927 : CPLErr
928 160 : GDALRegenerateOverviews( GDALRasterBandH hSrcBand,
929 : int nOverviewCount, GDALRasterBandH *pahOvrBands,
930 : const char * pszResampling,
931 : GDALProgressFunc pfnProgress, void * pProgressData )
932 :
933 : {
934 160 : GDALRasterBand *poSrcBand = (GDALRasterBand *) hSrcBand;
935 160 : GDALRasterBand **papoOvrBands = (GDALRasterBand **) pahOvrBands;
936 : int nFullResYChunk, nWidth;
937 : int nFRXBlockSize, nFRYBlockSize;
938 : GDALDataType eType;
939 : int bHasNoData;
940 : float fNoDataValue;
941 160 : GDALColorTable* poColorTable = NULL;
942 :
943 160 : if( pfnProgress == NULL )
944 4 : pfnProgress = GDALDummyProgress;
945 :
946 160 : if( EQUAL(pszResampling,"NONE") )
947 8 : return CE_None;
948 :
949 : /* -------------------------------------------------------------------- */
950 : /* Check color tables... */
951 : /* -------------------------------------------------------------------- */
952 219 : if ((EQUALN(pszResampling,"AVER",4)
953 : || EQUALN(pszResampling,"MODE",4)
954 : || EQUALN(pszResampling,"GAUSS",5)) &&
955 67 : poSrcBand->GetColorInterpretation() == GCI_PaletteIndex)
956 : {
957 6 : poColorTable = poSrcBand->GetColorTable();
958 6 : if (poColorTable != NULL)
959 : {
960 6 : if (poColorTable->GetPaletteInterpretation() != GPI_RGB)
961 : {
962 : CPLError(CE_Warning, CPLE_AppDefined,
963 : "Computing overviews on palette index raster bands "
964 : "with a palette whose color interpreation is not RGB "
965 0 : "will probably lead to unexpected results.");
966 0 : poColorTable = NULL;
967 : }
968 : }
969 : else
970 : {
971 : CPLError(CE_Warning, CPLE_AppDefined,
972 : "Computing overviews on palette index raster bands "
973 0 : "without a palette will probably lead to unexpected results.");
974 : }
975 : }
976 :
977 :
978 : /* If we have a nodata mask and we are doing something more complicated */
979 : /* than nearest neighbouring, we have to fetch to nodata mask */
980 : int bUseNoDataMask = (!EQUALN(pszResampling,"NEAR",4) &&
981 152 : (poSrcBand->GetMaskFlags() & GMF_ALL_VALID) == 0);
982 :
983 : /* -------------------------------------------------------------------- */
984 : /* If we are operating on multiple overviews, and using */
985 : /* averaging, lets do them in cascading order to reduce the */
986 : /* amount of computation. */
987 : /* -------------------------------------------------------------------- */
988 :
989 : /* In case the mask made be computed from another band of the dataset, */
990 : /* we can't use cascaded generation, as the computation of the overviews */
991 : /* of the band used for the mask band may not have yet occured (#3033) */
992 158 : if( (EQUALN(pszResampling,"AVER",4) || EQUALN(pszResampling,"GAUSS",5)) && nOverviewCount > 1
993 6 : && !(bUseNoDataMask && poSrcBand->GetMaskFlags() != GMF_NODATA))
994 : return GDALRegenerateCascadingOverviews( poSrcBand,
995 : nOverviewCount, papoOvrBands,
996 : pszResampling,
997 : pfnProgress,
998 12 : pProgressData );
999 :
1000 : /* -------------------------------------------------------------------- */
1001 : /* Setup one horizontal swath to read from the raw buffer. */
1002 : /* -------------------------------------------------------------------- */
1003 : float *pafChunk;
1004 140 : GByte *pabyChunkNodataMask = NULL;
1005 :
1006 140 : poSrcBand->GetBlockSize( &nFRXBlockSize, &nFRYBlockSize );
1007 :
1008 160 : if( nFRYBlockSize < 16 || nFRYBlockSize > 256 )
1009 20 : nFullResYChunk = 64;
1010 : else
1011 120 : nFullResYChunk = nFRYBlockSize;
1012 :
1013 140 : if( GDALDataTypeIsComplex( poSrcBand->GetRasterDataType() ) )
1014 0 : eType = GDT_CFloat32;
1015 : else
1016 140 : eType = GDT_Float32;
1017 :
1018 140 : nWidth = poSrcBand->GetXSize();
1019 : pafChunk = (float *)
1020 140 : VSIMalloc3((GDALGetDataTypeSize(eType)/8), nFullResYChunk, nWidth );
1021 140 : if (bUseNoDataMask)
1022 : {
1023 : pabyChunkNodataMask = (GByte *)
1024 17 : (GByte*) VSIMalloc2( nFullResYChunk, nWidth );
1025 : }
1026 :
1027 140 : if( pafChunk == NULL || (bUseNoDataMask && pabyChunkNodataMask == NULL))
1028 : {
1029 0 : CPLFree(pafChunk);
1030 0 : CPLFree(pabyChunkNodataMask);
1031 : CPLError( CE_Failure, CPLE_OutOfMemory,
1032 0 : "Out of memory in GDALRegenerateOverviews()." );
1033 :
1034 0 : return CE_Failure;
1035 : }
1036 :
1037 140 : fNoDataValue = (float) poSrcBand->GetNoDataValue(&bHasNoData);
1038 :
1039 : /* -------------------------------------------------------------------- */
1040 : /* Loop over image operating on chunks. */
1041 : /* -------------------------------------------------------------------- */
1042 140 : int nChunkYOff = 0;
1043 140 : CPLErr eErr = CE_None;
1044 :
1045 441 : for( nChunkYOff = 0;
1046 : nChunkYOff < poSrcBand->GetYSize() && eErr == CE_None;
1047 : nChunkYOff += nFullResYChunk )
1048 : {
1049 301 : if( !pfnProgress( nChunkYOff / (double) poSrcBand->GetYSize(),
1050 : NULL, pProgressData ) )
1051 : {
1052 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1053 0 : eErr = CE_Failure;
1054 : }
1055 :
1056 301 : if( nFullResYChunk + nChunkYOff > poSrcBand->GetYSize() )
1057 46 : nFullResYChunk = poSrcBand->GetYSize() - nChunkYOff;
1058 :
1059 : /* read chunk */
1060 301 : if (eErr == CE_None)
1061 : eErr = poSrcBand->RasterIO( GF_Read, 0, nChunkYOff, nWidth, nFullResYChunk,
1062 : pafChunk, nWidth, nFullResYChunk, eType,
1063 301 : 0, 0 );
1064 301 : if (eErr == CE_None && bUseNoDataMask)
1065 41 : eErr = poSrcBand->GetMaskBand()->RasterIO( GF_Read, 0, nChunkYOff, nWidth, nFullResYChunk,
1066 : pabyChunkNodataMask, nWidth, nFullResYChunk, GDT_Byte,
1067 41 : 0, 0 );
1068 :
1069 : /* special case to promote 1bit data to 8bit 0/255 values */
1070 301 : if( EQUAL(pszResampling,"AVERAGE_BIT2GRAYSCALE") )
1071 : {
1072 : int i;
1073 :
1074 39208 : for( i = nFullResYChunk*nWidth - 1; i >= 0; i-- )
1075 : {
1076 39204 : if( pafChunk[i] == 1.0 )
1077 23672 : pafChunk[i] = 255.0;
1078 : }
1079 : }
1080 297 : else if( EQUAL(pszResampling,"AVERAGE_BIT2GRAYSCALE_MINISWHITE") )
1081 : {
1082 : int i;
1083 :
1084 0 : for( i = nFullResYChunk*nWidth - 1; i >= 0; i-- )
1085 : {
1086 0 : if( pafChunk[i] == 1.0 )
1087 0 : pafChunk[i] = 0.0;
1088 0 : else if( pafChunk[i] == 0.0 )
1089 0 : pafChunk[i] = 255.0;
1090 : }
1091 : }
1092 :
1093 815 : for( int iOverview = 0; iOverview < nOverviewCount && eErr == CE_None; iOverview++ )
1094 : {
1095 514 : if( eType == GDT_Float32 )
1096 : eErr = GDALDownsampleChunk32R(nWidth, poSrcBand->GetYSize(),
1097 : pafChunk,
1098 : pabyChunkNodataMask,
1099 : 0, nWidth,
1100 : nChunkYOff, nFullResYChunk,
1101 : papoOvrBands[iOverview], pszResampling,
1102 : bHasNoData, fNoDataValue, poColorTable,
1103 514 : poSrcBand->GetRasterDataType());
1104 : else
1105 : eErr = GDALDownsampleChunkC32R(nWidth, poSrcBand->GetYSize(),
1106 : pafChunk, nChunkYOff, nFullResYChunk,
1107 0 : papoOvrBands[iOverview], pszResampling);
1108 : }
1109 : }
1110 :
1111 140 : VSIFree( pafChunk );
1112 140 : VSIFree( pabyChunkNodataMask );
1113 :
1114 : /* -------------------------------------------------------------------- */
1115 : /* Renormalized overview mean / stddev if needed. */
1116 : /* -------------------------------------------------------------------- */
1117 140 : if( eErr == CE_None && EQUAL(pszResampling,"AVERAGE_MP") )
1118 : {
1119 : GDALOverviewMagnitudeCorrection( (GDALRasterBandH) poSrcBand,
1120 : nOverviewCount,
1121 : (GDALRasterBandH *) papoOvrBands,
1122 0 : GDALDummyProgress, NULL );
1123 : }
1124 :
1125 : /* -------------------------------------------------------------------- */
1126 : /* It can be important to flush out data to overviews. */
1127 : /* -------------------------------------------------------------------- */
1128 348 : for( int iOverview = 0;
1129 : eErr == CE_None && iOverview < nOverviewCount;
1130 : iOverview++ )
1131 : {
1132 208 : eErr = papoOvrBands[iOverview]->FlushCache();
1133 : }
1134 :
1135 140 : if (eErr == CE_None)
1136 140 : pfnProgress( 1.0, NULL, pProgressData );
1137 :
1138 140 : return eErr;
1139 : }
1140 :
1141 :
1142 :
1143 : /************************************************************************/
1144 : /* GDALRegenerateOverviewsMultiBand() */
1145 : /************************************************************************/
1146 :
1147 : /**
1148 : * \brief Variant of GDALRegenerateOverviews, specialy dedicated for generating
1149 : * compressed pixel-interleaved overviews (JPEG-IN-TIFF for example)
1150 : *
1151 : * This function will generate one or more overview images from a base
1152 : * image using the requested downsampling algorithm. It's primary use
1153 : * is for generating overviews via GDALDataset::BuildOverviews(), but it
1154 : * can also be used to generate downsampled images in one file from another
1155 : * outside the overview architecture.
1156 : *
1157 : * The output bands need to exist in advance and share the same characteristics
1158 : * (type, dimensions)
1159 : *
1160 : * The resampling algorithms supported for the moment are "NEAREST", "AVERAGE"
1161 : * and "GAUSS"
1162 : *
1163 : * The pseudo-algorithm used by the function is :
1164 : * for each overview
1165 : * iterate on lines of the source by a step of deltay
1166 : * iterate on columns of the source by a step of deltax
1167 : * read the source data of size deltax * deltay for all the bands
1168 : * generate the corresponding overview block for all the bands
1169 : *
1170 : * This function will honour properly NODATA_VALUES tuples (special dataset metadata) so
1171 : * that only a given RGB triplet (in case of a RGB image) will be considered as the
1172 : * nodata value and not each value of the triplet independantly per band.
1173 : *
1174 : * @param nBands the number of bands, size of papoSrcBands and size of
1175 : * first dimension of papapoOverviewBands
1176 : * @param papoSrcBands the list of source bands to downsample
1177 : * @param nOverviews the number of downsampled overview levels being generated.
1178 : * @param papapoOverviewBands bidimension array of bands. First dimension is indexed
1179 : * by nBands. Second dimension is indexed by nOverviews.
1180 : * @param pszResampling Resampling algorithm ("NEAREST", "AVERAGE" or "GAUSS").
1181 : * @param pfnProgress progress report function.
1182 : * @param pProgressData progress function callback data.
1183 : * @return CE_None on success or CE_Failure on failure.
1184 : */
1185 :
1186 : CPLErr
1187 8 : GDALRegenerateOverviewsMultiBand(int nBands, GDALRasterBand** papoSrcBands,
1188 : int nOverviews,
1189 : GDALRasterBand*** papapoOverviewBands,
1190 : const char * pszResampling,
1191 : GDALProgressFunc pfnProgress, void * pProgressData )
1192 : {
1193 8 : CPLErr eErr = CE_None;
1194 : int iOverview, iBand;
1195 :
1196 8 : if( pfnProgress == NULL )
1197 0 : pfnProgress = GDALDummyProgress;
1198 :
1199 8 : if( EQUAL(pszResampling,"NONE") )
1200 0 : return CE_None;
1201 :
1202 : /* Sanity checks */
1203 8 : if (!EQUALN(pszResampling, "NEAR", 4) && !EQUAL(pszResampling, "AVERAGE") && !EQUAL(pszResampling, "GAUSS"))
1204 : {
1205 : CPLError(CE_Failure, CPLE_NotSupported,
1206 0 : "GDALRegenerateOverviewsMultiBand: pszResampling='%s' not supported", pszResampling);
1207 0 : return CE_Failure;
1208 : }
1209 :
1210 8 : int nSrcWidth = papoSrcBands[0]->GetXSize();
1211 8 : int nSrcHeight = papoSrcBands[0]->GetYSize();
1212 8 : GDALDataType eDataType = papoSrcBands[0]->GetRasterDataType();
1213 24 : for(iBand=1;iBand<nBands;iBand++)
1214 : {
1215 32 : if (papoSrcBands[iBand]->GetXSize() != nSrcWidth ||
1216 16 : papoSrcBands[iBand]->GetYSize() != nSrcHeight)
1217 : {
1218 : CPLError(CE_Failure, CPLE_NotSupported,
1219 0 : "GDALRegenerateOverviewsMultiBand: all the source bands must have the same dimensions");
1220 0 : return CE_Failure;
1221 : }
1222 16 : if (papoSrcBands[iBand]->GetRasterDataType() != eDataType)
1223 : {
1224 : CPLError(CE_Failure, CPLE_NotSupported,
1225 0 : "GDALRegenerateOverviewsMultiBand: all the source bands must have the same data type");
1226 0 : return CE_Failure;
1227 : }
1228 : }
1229 :
1230 16 : for(iOverview=0;iOverview<nOverviews;iOverview++)
1231 : {
1232 8 : int nDstWidth = papapoOverviewBands[0][iOverview]->GetXSize();
1233 8 : int nDstHeight = papapoOverviewBands[0][iOverview]->GetYSize();
1234 24 : for(iBand=1;iBand<nBands;iBand++)
1235 : {
1236 48 : if (papapoOverviewBands[iBand][iOverview]->GetXSize() != nDstWidth ||
1237 32 : papapoOverviewBands[iBand][iOverview]->GetYSize() != nDstHeight)
1238 : {
1239 : CPLError(CE_Failure, CPLE_NotSupported,
1240 0 : "GDALRegenerateOverviewsMultiBand: all the overviews bands of the same level must have the same dimensions");
1241 0 : return CE_Failure;
1242 : }
1243 16 : if (papapoOverviewBands[iBand][iOverview]->GetRasterDataType() != eDataType)
1244 : {
1245 : CPLError(CE_Failure, CPLE_NotSupported,
1246 0 : "GDALRegenerateOverviewsMultiBand: all the overviews bands must have the same data type as the source bands");
1247 0 : return CE_Failure;
1248 : }
1249 : }
1250 : }
1251 :
1252 : /* First pass to compute the total number of pixels to read */
1253 8 : double dfTotalPixelCount = 0;
1254 16 : for(iOverview=0;iOverview<nOverviews;iOverview++)
1255 : {
1256 8 : nSrcWidth = papoSrcBands[0]->GetXSize();
1257 8 : nSrcHeight = papoSrcBands[0]->GetYSize();
1258 :
1259 8 : int nDstWidth = papapoOverviewBands[0][iOverview]->GetXSize();
1260 : /* Try to use previous level of overview as the source to compute */
1261 : /* the next level */
1262 8 : if (iOverview > 0 && papapoOverviewBands[0][iOverview - 1]->GetXSize() > nDstWidth)
1263 : {
1264 0 : nSrcWidth = papapoOverviewBands[0][iOverview - 1]->GetXSize();
1265 0 : nSrcHeight = papapoOverviewBands[0][iOverview - 1]->GetYSize();
1266 : }
1267 :
1268 8 : dfTotalPixelCount += (double)nSrcWidth * nSrcHeight;
1269 : }
1270 :
1271 8 : nSrcWidth = papoSrcBands[0]->GetXSize();
1272 8 : nSrcHeight = papoSrcBands[0]->GetYSize();
1273 :
1274 : /* If we have a nodata mask and we are doing something more complicated */
1275 : /* than nearest neighbouring, we have to fetch to nodata mask */
1276 : int bUseNoDataMask = (!EQUALN(pszResampling,"NEAR",4) &&
1277 8 : (papoSrcBands[0]->GetMaskFlags() & GMF_ALL_VALID) == 0);
1278 :
1279 8 : int* pabHasNoData = (int*)CPLMalloc(nBands * sizeof(int));
1280 8 : float* pafNoDataValue = (float*)CPLMalloc(nBands * sizeof(float));
1281 :
1282 32 : for(iBand=0;iBand<nBands;iBand++)
1283 : {
1284 24 : pabHasNoData[iBand] = FALSE;
1285 24 : pafNoDataValue[iBand] = (float) papoSrcBands[iBand]->GetNoDataValue(&pabHasNoData[iBand]);
1286 : }
1287 :
1288 : /* Second pass to do the real job ! */
1289 8 : double dfCurPixelCount = 0;
1290 16 : for(iOverview=0;iOverview<nOverviews && eErr == CE_None;iOverview++)
1291 : {
1292 8 : int iSrcOverview = -1; /* -1 means the source bands */
1293 :
1294 : int nDstBlockXSize, nDstBlockYSize;
1295 : int nDstWidth, nDstHeight;
1296 8 : papapoOverviewBands[0][iOverview]->GetBlockSize(&nDstBlockXSize, &nDstBlockYSize);
1297 8 : nDstWidth = papapoOverviewBands[0][iOverview]->GetXSize();
1298 8 : nDstHeight = papapoOverviewBands[0][iOverview]->GetYSize();
1299 :
1300 : /* Try to use previous level of overview as the source to compute */
1301 : /* the next level */
1302 8 : if (iOverview > 0 && papapoOverviewBands[0][iOverview - 1]->GetXSize() > nDstWidth)
1303 : {
1304 0 : nSrcWidth = papapoOverviewBands[0][iOverview - 1]->GetXSize();
1305 0 : nSrcHeight = papapoOverviewBands[0][iOverview - 1]->GetYSize();
1306 0 : iSrcOverview = iOverview - 1;
1307 : }
1308 :
1309 : /* Compute the chunck size of the source such as it will match the size of */
1310 : /* a block of the overview */
1311 8 : int nFullResXChunk = (nDstBlockXSize * nSrcWidth) / nDstWidth;
1312 8 : int nFullResYChunk = (nDstBlockYSize * nSrcHeight) / nDstHeight;
1313 :
1314 8 : float** papafChunk = (float**) CPLMalloc(nBands * sizeof(void*));
1315 8 : GByte* pabyChunkNoDataMask = NULL;
1316 32 : for(iBand=0;iBand<nBands;iBand++)
1317 : {
1318 24 : papafChunk[iBand] = (float*) VSIMalloc3(nFullResXChunk, nFullResYChunk, sizeof(float));
1319 24 : if( papafChunk[iBand] == NULL )
1320 : {
1321 0 : while ( --iBand >= 0)
1322 0 : CPLFree(papafChunk[iBand]);
1323 0 : CPLFree(papafChunk);
1324 0 : CPLFree(pabHasNoData);
1325 0 : CPLFree(pafNoDataValue);
1326 :
1327 : CPLError( CE_Failure, CPLE_OutOfMemory,
1328 0 : "GDALRegenerateOverviewsMultiBand: Out of memory." );
1329 0 : return CE_Failure;
1330 : }
1331 : }
1332 8 : if (bUseNoDataMask)
1333 : {
1334 4 : pabyChunkNoDataMask = (GByte*) VSIMalloc2(nFullResXChunk, nFullResYChunk);
1335 4 : if( pabyChunkNoDataMask == NULL )
1336 : {
1337 0 : for(iBand=0;iBand<nBands;iBand++)
1338 : {
1339 0 : CPLFree(papafChunk[iBand]);
1340 : }
1341 0 : CPLFree(papafChunk);
1342 0 : CPLFree(pabHasNoData);
1343 0 : CPLFree(pafNoDataValue);
1344 :
1345 : CPLError( CE_Failure, CPLE_OutOfMemory,
1346 0 : "GDALRegenerateOverviewsMultiBand: Out of memory." );
1347 0 : return CE_Failure;
1348 : }
1349 : }
1350 :
1351 : int nChunkYOff;
1352 : /* Iterate on destination overview, block by block */
1353 16 : for( nChunkYOff = 0; nChunkYOff < nSrcHeight && eErr == CE_None; nChunkYOff += nFullResYChunk )
1354 : {
1355 : int nYCount;
1356 8 : if (nChunkYOff + nFullResYChunk <= nSrcHeight)
1357 0 : nYCount = nFullResYChunk;
1358 : else
1359 8 : nYCount = nSrcHeight - nChunkYOff;
1360 :
1361 8 : if( !pfnProgress( dfCurPixelCount / dfTotalPixelCount,
1362 : NULL, pProgressData ) )
1363 : {
1364 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1365 0 : eErr = CE_Failure;
1366 : }
1367 :
1368 : int nChunkXOff;
1369 16 : for( nChunkXOff = 0; nChunkXOff < nSrcWidth && eErr == CE_None; nChunkXOff += nFullResXChunk )
1370 : {
1371 : int nXCount;
1372 8 : if (nChunkXOff + nFullResXChunk <= nSrcWidth)
1373 0 : nXCount = nFullResXChunk;
1374 : else
1375 8 : nXCount = nSrcWidth - nChunkXOff;
1376 :
1377 : /* Read the source buffers for all the bands */
1378 32 : for(iBand=0;iBand<nBands && eErr == CE_None;iBand++)
1379 : {
1380 : GDALRasterBand* poSrcBand;
1381 24 : if (iSrcOverview == -1)
1382 24 : poSrcBand = papoSrcBands[iBand];
1383 : else
1384 0 : poSrcBand = papapoOverviewBands[iBand][iSrcOverview];
1385 : eErr = poSrcBand->RasterIO( GF_Read,
1386 : nChunkXOff, nChunkYOff,
1387 : nXCount, nYCount,
1388 24 : papafChunk[iBand],
1389 : nXCount, nYCount,
1390 48 : GDT_Float32, 0, 0 );
1391 : }
1392 :
1393 8 : if (bUseNoDataMask && eErr == CE_None)
1394 : {
1395 : GDALRasterBand* poSrcBand;
1396 4 : if (iSrcOverview == -1)
1397 4 : poSrcBand = papoSrcBands[0];
1398 : else
1399 0 : poSrcBand = papapoOverviewBands[0][iSrcOverview];
1400 4 : eErr = poSrcBand->GetMaskBand()->RasterIO( GF_Read,
1401 : nChunkXOff, nChunkYOff,
1402 : nXCount, nYCount,
1403 : pabyChunkNoDataMask,
1404 : nXCount, nYCount,
1405 4 : GDT_Byte, 0, 0 );
1406 : }
1407 :
1408 : /* Compute the resulting overview block */
1409 32 : for(iBand=0;iBand<nBands && eErr == CE_None;iBand++)
1410 : {
1411 : eErr = GDALDownsampleChunk32R(nSrcWidth, nSrcHeight,
1412 : papafChunk[iBand],
1413 : pabyChunkNoDataMask,
1414 : nChunkXOff, nXCount,
1415 : nChunkYOff, nYCount,
1416 24 : papapoOverviewBands[iBand][iOverview],
1417 : pszResampling,
1418 : pabHasNoData[iBand],
1419 : pafNoDataValue[iBand],
1420 : /*poColorTable*/ NULL,
1421 48 : eDataType);
1422 : }
1423 : }
1424 :
1425 8 : dfCurPixelCount += (double)nYCount * nSrcWidth;
1426 : }
1427 :
1428 : /* Flush the data to overviews */
1429 32 : for(iBand=0;iBand<nBands;iBand++)
1430 : {
1431 24 : CPLFree(papafChunk[iBand]);
1432 24 : papapoOverviewBands[iBand][iOverview]->FlushCache();
1433 : }
1434 8 : CPLFree(papafChunk);
1435 8 : CPLFree(pabyChunkNoDataMask);
1436 :
1437 : }
1438 :
1439 8 : CPLFree(pabHasNoData);
1440 8 : CPLFree(pafNoDataValue);
1441 :
1442 8 : if (eErr == CE_None)
1443 8 : pfnProgress( 1.0, NULL, pProgressData );
1444 :
1445 8 : return eErr;
1446 : }
1447 :
1448 :
1449 : /************************************************************************/
1450 : /* GDALComputeBandStats() */
1451 : /************************************************************************/
1452 :
1453 : CPLErr CPL_STDCALL
1454 11 : GDALComputeBandStats( GDALRasterBandH hSrcBand,
1455 : int nSampleStep,
1456 : double *pdfMean, double *pdfStdDev,
1457 : GDALProgressFunc pfnProgress,
1458 : void *pProgressData )
1459 :
1460 : {
1461 11 : VALIDATE_POINTER1( hSrcBand, "GDALComputeBandStats", CE_Failure );
1462 :
1463 11 : GDALRasterBand *poSrcBand = (GDALRasterBand *) hSrcBand;
1464 : int iLine, nWidth, nHeight;
1465 11 : GDALDataType eType = poSrcBand->GetRasterDataType();
1466 : GDALDataType eWrkType;
1467 : int bComplex;
1468 : float *pafData;
1469 11 : double dfSum=0.0, dfSum2=0.0;
1470 11 : int nSamples = 0;
1471 :
1472 11 : if( pfnProgress == NULL )
1473 11 : pfnProgress = GDALDummyProgress;
1474 :
1475 11 : nWidth = poSrcBand->GetXSize();
1476 11 : nHeight = poSrcBand->GetYSize();
1477 :
1478 11 : if( nSampleStep >= nHeight || nSampleStep < 1 )
1479 0 : nSampleStep = 1;
1480 :
1481 11 : bComplex = GDALDataTypeIsComplex(eType);
1482 11 : if( bComplex )
1483 : {
1484 0 : pafData = (float *) VSIMalloc(nWidth * 2 * sizeof(float));
1485 0 : eWrkType = GDT_CFloat32;
1486 : }
1487 : else
1488 : {
1489 11 : pafData = (float *) VSIMalloc(nWidth * sizeof(float));
1490 11 : eWrkType = GDT_Float32;
1491 : }
1492 :
1493 11 : if( pafData == NULL )
1494 : {
1495 : CPLError( CE_Failure, CPLE_OutOfMemory,
1496 0 : "GDALComputeBandStats: Out of memory for buffer." );
1497 0 : return CE_Failure;
1498 : }
1499 :
1500 : /* -------------------------------------------------------------------- */
1501 : /* Loop over all sample lines. */
1502 : /* -------------------------------------------------------------------- */
1503 1201 : for( iLine = 0; iLine < nHeight; iLine += nSampleStep )
1504 : {
1505 : int iPixel;
1506 :
1507 1191 : if( !pfnProgress( iLine / (double) nHeight,
1508 : NULL, pProgressData ) )
1509 : {
1510 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1511 0 : CPLFree( pafData );
1512 0 : return CE_Failure;
1513 : }
1514 :
1515 : CPLErr eErr = poSrcBand->RasterIO( GF_Read, 0, iLine, nWidth, 1,
1516 : pafData, nWidth, 1, eWrkType,
1517 1191 : 0, 0 );
1518 1191 : if ( eErr != CE_None )
1519 : {
1520 1 : CPLFree( pafData );
1521 1 : return eErr;
1522 : }
1523 :
1524 353030 : for( iPixel = 0; iPixel < nWidth; iPixel++ )
1525 : {
1526 : float fValue;
1527 :
1528 351840 : if( bComplex )
1529 : {
1530 : // Compute the magnitude of the complex value.
1531 :
1532 : fValue = (float)
1533 0 : sqrt(pafData[iPixel*2 ] * pafData[iPixel*2 ]
1534 0 : + pafData[iPixel*2+1] * pafData[iPixel*2+1]);
1535 : }
1536 : else
1537 : {
1538 351840 : fValue = pafData[iPixel];
1539 : }
1540 :
1541 351840 : dfSum += fValue;
1542 351840 : dfSum2 += fValue * fValue;
1543 : }
1544 :
1545 1190 : nSamples += nWidth;
1546 : }
1547 :
1548 10 : if( !pfnProgress( 1.0, NULL, pProgressData ) )
1549 : {
1550 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1551 0 : CPLFree( pafData );
1552 0 : return CE_Failure;
1553 : }
1554 :
1555 : /* -------------------------------------------------------------------- */
1556 : /* Produce the result values. */
1557 : /* -------------------------------------------------------------------- */
1558 10 : if( pdfMean != NULL )
1559 10 : *pdfMean = dfSum / nSamples;
1560 :
1561 10 : if( pdfStdDev != NULL )
1562 : {
1563 10 : double dfMean = dfSum / nSamples;
1564 :
1565 10 : *pdfStdDev = sqrt((dfSum2 / nSamples) - (dfMean * dfMean));
1566 : }
1567 :
1568 10 : CPLFree( pafData );
1569 :
1570 10 : return CE_None;
1571 : }
1572 :
1573 : /************************************************************************/
1574 : /* GDALOverviewMagnitudeCorrection() */
1575 : /* */
1576 : /* Correct the mean and standard deviation of the overviews of */
1577 : /* the given band to match the base layer approximately. */
1578 : /************************************************************************/
1579 :
1580 : CPLErr
1581 0 : GDALOverviewMagnitudeCorrection( GDALRasterBandH hBaseBand,
1582 : int nOverviewCount,
1583 : GDALRasterBandH *pahOverviews,
1584 : GDALProgressFunc pfnProgress,
1585 : void *pProgressData )
1586 :
1587 : {
1588 0 : VALIDATE_POINTER1( hBaseBand, "GDALOverviewMagnitudeCorrection", CE_Failure );
1589 :
1590 : CPLErr eErr;
1591 : double dfOrigMean, dfOrigStdDev;
1592 :
1593 : /* -------------------------------------------------------------------- */
1594 : /* Compute mean/stddev for source raster. */
1595 : /* -------------------------------------------------------------------- */
1596 : eErr = GDALComputeBandStats( hBaseBand, 2, &dfOrigMean, &dfOrigStdDev,
1597 0 : pfnProgress, pProgressData );
1598 :
1599 0 : if( eErr != CE_None )
1600 0 : return eErr;
1601 :
1602 : /* -------------------------------------------------------------------- */
1603 : /* Loop on overview bands. */
1604 : /* -------------------------------------------------------------------- */
1605 : int iOverview;
1606 :
1607 0 : for( iOverview = 0; iOverview < nOverviewCount; iOverview++ )
1608 : {
1609 0 : GDALRasterBand *poOverview = (GDALRasterBand *)pahOverviews[iOverview];
1610 : double dfOverviewMean, dfOverviewStdDev;
1611 : double dfGain;
1612 :
1613 : eErr = GDALComputeBandStats( pahOverviews[iOverview], 1,
1614 : &dfOverviewMean, &dfOverviewStdDev,
1615 0 : pfnProgress, pProgressData );
1616 :
1617 0 : if( eErr != CE_None )
1618 0 : return eErr;
1619 :
1620 0 : if( dfOrigStdDev < 0.0001 )
1621 0 : dfGain = 1.0;
1622 : else
1623 0 : dfGain = dfOrigStdDev / dfOverviewStdDev;
1624 :
1625 : /* -------------------------------------------------------------------- */
1626 : /* Apply gain and offset. */
1627 : /* -------------------------------------------------------------------- */
1628 0 : GDALDataType eWrkType, eType = poOverview->GetRasterDataType();
1629 : int iLine, nWidth, nHeight, bComplex;
1630 : float *pafData;
1631 :
1632 0 : nWidth = poOverview->GetXSize();
1633 0 : nHeight = poOverview->GetYSize();
1634 :
1635 0 : bComplex = GDALDataTypeIsComplex(eType);
1636 0 : if( bComplex )
1637 : {
1638 0 : pafData = (float *) VSIMalloc2(nWidth, 2 * sizeof(float));
1639 0 : eWrkType = GDT_CFloat32;
1640 : }
1641 : else
1642 : {
1643 0 : pafData = (float *) VSIMalloc2(nWidth, sizeof(float));
1644 0 : eWrkType = GDT_Float32;
1645 : }
1646 :
1647 0 : if( pafData == NULL )
1648 : {
1649 : CPLError( CE_Failure, CPLE_OutOfMemory,
1650 0 : "GDALOverviewMagnitudeCorrection: Out of memory for buffer." );
1651 0 : return CE_Failure;
1652 : }
1653 :
1654 0 : for( iLine = 0; iLine < nHeight; iLine++ )
1655 : {
1656 : int iPixel;
1657 :
1658 0 : if( !pfnProgress( iLine / (double) nHeight,
1659 : NULL, pProgressData ) )
1660 : {
1661 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1662 0 : CPLFree( pafData );
1663 0 : return CE_Failure;
1664 : }
1665 :
1666 : poOverview->RasterIO( GF_Read, 0, iLine, nWidth, 1,
1667 : pafData, nWidth, 1, eWrkType,
1668 0 : 0, 0 );
1669 :
1670 0 : for( iPixel = 0; iPixel < nWidth; iPixel++ )
1671 : {
1672 0 : if( bComplex )
1673 : {
1674 0 : pafData[iPixel*2] *= (float) dfGain;
1675 0 : pafData[iPixel*2+1] *= (float) dfGain;
1676 : }
1677 : else
1678 : {
1679 0 : pafData[iPixel] = (float)
1680 0 : ((pafData[iPixel]-dfOverviewMean)*dfGain + dfOrigMean);
1681 :
1682 : }
1683 : }
1684 :
1685 : poOverview->RasterIO( GF_Write, 0, iLine, nWidth, 1,
1686 : pafData, nWidth, 1, eWrkType,
1687 0 : 0, 0 );
1688 : }
1689 :
1690 0 : if( !pfnProgress( 1.0, NULL, pProgressData ) )
1691 : {
1692 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1693 0 : CPLFree( pafData );
1694 0 : return CE_Failure;
1695 : }
1696 :
1697 0 : CPLFree( pafData );
1698 : }
1699 :
1700 0 : return CE_None;
1701 : }
|
