00001 00002 00003 00004 00005 00006 00007 00008 00009 00010 00011 00012 00013 00014 00015 00016 00017 00018 00019 00020 00021 00022 00023 00024 00025 00026 00027 00028 00029 00030 00030 00030 00030 00031 00032 00033 00034 00035 00036 00037 00038 00039 00040 00041 00042 00043
00044
00045 #include "gdal_priv.h"
00046
00047
00048
00049
00050
00051 static CPLErr
00052 GDALDownsampleChunk32R( int nSrcWidth, int nSrcHeight,
00053 float * pafChunk, int nChunkYOff, int nChunkYSize,
00054 GDALRasterBand * poOverview,
00055 const char * pszResampling )
00056
00057 {
00058 int nDstYOff, nDstYOff2, nOXSize, nOYSize;
00059 float *pafDstScanline;
00060
00061 nOXSize = poOverview->GetXSize();
00062 nOYSize = poOverview->GetYSize();
00063
00064 pafDstScanline = (float *) CPLMalloc(nOXSize * sizeof(float));
00065
00066
00067
00068
00069
00070
00071
00072 nDstYOff = (int) (0.5 + (nChunkYOff/(double)nSrcHeight) * nOYSize);
00073 nDstYOff2 = (int)
00074 (0.5 + ((nChunkYOff+nChunkYSize)/(double)nSrcHeight) * nOYSize);
00075
00076 if( nChunkYOff + nChunkYSize == nSrcHeight )
00077 nDstYOff2 = nOYSize;
00078
00079
00080
00081
00082 for( int iDstLine = nDstYOff; iDstLine < nDstYOff2; iDstLine++ )
00083 {
00084 float *pafSrcScanline;
00085 int nSrcYOff, nSrcYOff2, iDstPixel;
00086
00087 nSrcYOff = (int) (0.5 + (iDstLine/(double)nOYSize) * nSrcHeight);
00088 if( nSrcYOff < nChunkYOff )
00089 nSrcYOff = nChunkYOff;
00090
00091 nSrcYOff2 = (int) (0.5 + ((iDstLine+1)/(double)nOYSize) * nSrcHeight);
00092 if( nSrcYOff2 > nSrcHeight || iDstLine == nOYSize-1 )
00093 nSrcYOff2 = nSrcHeight;
00094 if( nSrcYOff2 > nChunkYOff + nChunkYSize )
00095 nSrcYOff2 = nChunkYOff + nChunkYSize;
00096
00097 pafSrcScanline = pafChunk + ((nSrcYOff-nChunkYOff) * nSrcWidth);
00098
00099
00100
00101
00102 for( iDstPixel = 0; iDstPixel < nOXSize; iDstPixel++ )
00103 {
00104 int nSrcXOff, nSrcXOff2;
00105
00106 nSrcXOff = (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
00107 nSrcXOff2 = (int)
00108 (0.5 + ((iDstPixel+1)/(double)nOXSize) * nSrcWidth);
00109 if( nSrcXOff2 > nSrcWidth )
00110 nSrcXOff2 = nSrcWidth;
00111
00112 if( EQUALN(pszResampling,"NEAR",4) )
00113 {
00114 pafDstScanline[iDstPixel] = pafSrcScanline[nSrcXOff];
00115 }
00116 else if( EQUALN(pszResampling,"AVER",4) )
00117 {
00118 double dfTotal = 0.0;
00119 int nCount = 0, iX, iY;
00120
00121 for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
00122 {
00123 for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
00124 {
00125 dfTotal += pafSrcScanline[iX+(iY-nSrcYOff)*nSrcWidth];
00126 nCount++;
00127 }
00128 }
00129
00130 CPLAssert( nCount > 0 );
00131 if( nCount == 0 )
00132 {
00133 pafDstScanline[iDstPixel] = 0.0;
00134 }
00135 else
00136 pafDstScanline[iDstPixel] = dfTotal / nCount;
00137 }
00138 }
00139
00140 poOverview->RasterIO( GF_Write, 0, iDstLine, nOXSize, 1,
00141 pafDstScanline, nOXSize, 1, GDT_Float32,
00142 0, 0 );
00143 }
00144
00145 CPLFree( pafDstScanline );
00146
00147 return CE_None;
00148 }
00149
00150
00151
00152
00153
00154 static CPLErr
00155 GDALDownsampleChunkC32R( int nSrcWidth, int nSrcHeight,
00156 float * pafChunk, int nChunkYOff, int nChunkYSize,
00157 GDALRasterBand * poOverview,
00158 const char * pszResampling )
00159
00160 {
00161 int nDstYOff, nDstYOff2, nOXSize, nOYSize;
00162 float *pafDstScanline;
00163
00164 nOXSize = poOverview->GetXSize();
00165 nOYSize = poOverview->GetYSize();
00166
00167 pafDstScanline = (float *) CPLMalloc(nOXSize * sizeof(float) * 2);
00168
00169
00170
00171
00172
00173
00174
00175 nDstYOff = (int) (0.5 + (nChunkYOff/(double)nSrcHeight) * nOYSize);
00176 nDstYOff2 = (int)
00177 (0.5 + ((nChunkYOff+nChunkYSize)/(double)nSrcHeight) * nOYSize);
00178
00179 if( nChunkYOff + nChunkYSize == nSrcHeight )
00180 nDstYOff2 = nOYSize;
00181
00182
00183
00184
00185 for( int iDstLine = nDstYOff; iDstLine < nDstYOff2; iDstLine++ )
00186 {
00187 float *pafSrcScanline;
00188 int nSrcYOff, nSrcYOff2, iDstPixel;
00189
00190 nSrcYOff = (int) (0.5 + (iDstLine/(double)nOYSize) * nSrcHeight);
00191 if( nSrcYOff < nChunkYOff )
00192 nSrcYOff = nChunkYOff;
00193
00194 nSrcYOff2 = (int) (0.5 + ((iDstLine+1)/(double)nOYSize) * nSrcHeight);
00195 if( nSrcYOff2 > nSrcHeight || iDstLine == nOYSize-1 )
00196 nSrcYOff2 = nSrcHeight;
00197 if( nSrcYOff2 > nChunkYOff + nChunkYSize )
00198 nSrcYOff2 = nChunkYOff + nChunkYSize;
00199
00200 pafSrcScanline = pafChunk + ((nSrcYOff-nChunkYOff) * nSrcWidth);
00201
00202
00203
00204
00205 for( iDstPixel = 0; iDstPixel < nOXSize; iDstPixel++ )
00206 {
00207 int nSrcXOff, nSrcXOff2;
00208
00209 nSrcXOff = (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
00210 nSrcXOff2 = (int)
00211 (0.5 + ((iDstPixel+1)/(double)nOXSize) * nSrcWidth);
00212 if( nSrcXOff2 > nSrcWidth )
00213 nSrcXOff2 = nSrcWidth;
00214
00215 if( EQUALN(pszResampling,"NEAR",4) )
00216 {
00217 pafDstScanline[iDstPixel*2] = pafSrcScanline[nSrcXOff*2];
00218 pafDstScanline[iDstPixel*2+1] = pafSrcScanline[nSrcXOff*2+1];
00219 }
00220 else if( EQUALN(pszResampling,"AVER",4) )
00221 {
00222 double dfTotalR = 0.0, dfTotalI = 0.0;
00223 int nCount = 0, iX, iY;
00224
00225 for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
00226 {
00227 for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
00228 {
00229 dfTotalR += pafSrcScanline[iX*2+(iY-nSrcYOff)*nSrcWidth*2];
00230 dfTotalI += pafSrcScanline[iX*2+(iY-nSrcYOff)*nSrcWidth*2+1];
00231 nCount++;
00232 }
00233 }
00234
00235 CPLAssert( nCount > 0 );
00236 if( nCount == 0 )
00237 {
00238 pafDstScanline[iDstPixel*2] = 0.0;
00239 pafDstScanline[iDstPixel*2+1] = 0.0;
00240 }
00241 else
00242 {
00243 pafDstScanline[iDstPixel*2 ] = dfTotalR / nCount;
00244 pafDstScanline[iDstPixel*2+1] = dfTotalI / nCount;
00245 }
00246 }
00247 }
00248
00249 poOverview->RasterIO( GF_Write, 0, iDstLine, nOXSize, 1,
00250 pafDstScanline, nOXSize, 1, GDT_CFloat32,
00251 0, 0 );
00252 }
00253
00254 CPLFree( pafDstScanline );
00255
00256 return CE_None;
00257 }
00258
00259
00260
00261
00262
00263
00264
00265
00266 static CPLErr
00267 GDALRegenerateCascadingOverviews(
00268 GDALRasterBand *poSrcBand, int nOverviews, GDALRasterBand **papoOvrBands,
00269 const char * pszResampling,
00270 GDALProgressFunc pfnProgress, void * pProgressData )
00271
00272 {
00273
00274
00275
00276
00277 int i, j;
00278
00279 for( i = 0; i < nOverviews-1; i++ )
00280 {
00281 for( j = 0; j < nOverviews - i - 1; j++ )
00282 {
00283
00284 if( papoOvrBands[j]->GetXSize()
00285 * (float) papoOvrBands[j]->GetYSize() <
00286 papoOvrBands[j+1]->GetXSize()
00287 * (float) papoOvrBands[j+1]->GetYSize() )
00288 {
00289 GDALRasterBand * poTempBand;
00290
00291 poTempBand = papoOvrBands[j];
00292 papoOvrBands[j] = papoOvrBands[j+1];
00293 papoOvrBands[j+1] = poTempBand;
00294 }
00295 }
00296 }
00297
00298
00299
00300
00301
00302 double dfTotalPixels = 0.0;
00303
00304 for( i = 0; i < nOverviews; i++ )
00305 {
00306 dfTotalPixels += papoOvrBands[i]->GetXSize()
00307 * (double) papoOvrBands[i]->GetYSize();
00308 }
00309
00310
00311
00312
00313 double dfPixelsProcessed = 0.0;
00314
00315 for( i = 0; i < nOverviews; i++ )
00316 {
00317 void *pScaledProgressData;
00318 double dfPixels;
00319 GDALRasterBand *poBaseBand;
00320 CPLErr eErr;
00321
00322 if( i == 0 )
00323 poBaseBand = poSrcBand;
00324 else
00325 poBaseBand = papoOvrBands[i-1];
00326
00327 dfPixels = papoOvrBands[i]->GetXSize()
00328 * (double) papoOvrBands[i]->GetYSize();
00329
00330 pScaledProgressData = GDALCreateScaledProgress(
00331 dfPixelsProcessed / dfTotalPixels,
00332 (dfPixelsProcessed + dfPixels) / dfTotalPixels,
00333 pfnProgress, pProgressData );
00334
00335 eErr = GDALRegenerateOverviews( poBaseBand, 1, papoOvrBands + i,
00336 pszResampling,
00337 GDALScaledProgress,
00338 pScaledProgressData );
00339 GDALDestroyScaledProgress( pScaledProgressData );
00340
00341 if( eErr != CE_None )
00342 return eErr;
00343
00344 dfPixelsProcessed += dfPixels;
00345 }
00346
00347 return CE_None;
00348 }
00349
00350
00351
00352
00353 CPLErr
00354 GDALRegenerateOverviews( GDALRasterBand *poSrcBand,
00355 int nOverviews, GDALRasterBand **papoOvrBands,
00356 const char * pszResampling,
00357 GDALProgressFunc pfnProgress, void * pProgressData )
00358
00359 {
00360 int nFullResYChunk, nWidth;
00361 int nFRXBlockSize, nFRYBlockSize;
00362 GDALDataType eType;
00363
00364
00365
00366
00367
00368
00369 if( EQUALN(pszResampling,"AVER",4) && nOverviews > 1 )
00370 return GDALRegenerateCascadingOverviews( poSrcBand,
00371 nOverviews, papoOvrBands,
00372 pszResampling,
00373 pfnProgress,
00374 pProgressData );
00375
00376
00377
00378
00379 float *pafChunk;
00380
00381 poSrcBand->GetBlockSize( &nFRXBlockSize, &nFRYBlockSize );
00382
00383 if( nFRYBlockSize < 4 || nFRYBlockSize > 256 )
00384 nFullResYChunk = 32;
00385 else
00386 nFullResYChunk = nFRYBlockSize;
00387
00388 if( GDALDataTypeIsComplex( poSrcBand->GetRasterDataType() ) )
00389 eType = GDT_CFloat32;
00390 else
00391 eType = GDT_Float32;
00392
00393 nWidth = poSrcBand->GetXSize();
00394 pafChunk = (float *)
00395 VSIMalloc((GDALGetDataTypeSize(eType)/8) * nFullResYChunk * nWidth );
00396
00397 if( pafChunk == NULL )
00398 {
00399 CPLError( CE_Failure, CPLE_OutOfMemory,
00400 "Out of memory in GDALRegenerateOverviews()." );
00401
00402 return CE_Failure;
00403 }
00404
00405
00406
00407
00408 int nChunkYOff = 0;
00409
00410 for( nChunkYOff = 0;
00411 nChunkYOff < poSrcBand->GetYSize();
00412 nChunkYOff += nFullResYChunk )
00413 {
00414 if( !pfnProgress( nChunkYOff / (double) poSrcBand->GetYSize(),
00415 NULL, pProgressData ) )
00416 {
00417 CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
00418 return CE_Failure;
00419 }
00420
00421 if( nFullResYChunk + nChunkYOff > poSrcBand->GetYSize() )
00422 nFullResYChunk = poSrcBand->GetYSize() - nChunkYOff;
00423
00424
00425 poSrcBand->RasterIO( GF_Read, 0, nChunkYOff, nWidth, nFullResYChunk,
00426 pafChunk, nWidth, nFullResYChunk, eType,
00427 0, 0 );
00428
00429 for( int iOverview = 0; iOverview < nOverviews; iOverview++ )
00430 {
00431 if( eType == GDT_Float32 )
00432 GDALDownsampleChunk32R(nWidth, poSrcBand->GetYSize(),
00433 pafChunk, nChunkYOff, nFullResYChunk,
00434 papoOvrBands[iOverview], pszResampling);
00435 else
00436 GDALDownsampleChunkC32R(nWidth, poSrcBand->GetYSize(),
00437 pafChunk, nChunkYOff, nFullResYChunk,
00438 papoOvrBands[iOverview], pszResampling);
00439 }
00440 }
00441
00442 VSIFree( pafChunk );
00443
00444
00445
00446
00447 for( int iOverview = 0; iOverview < nOverviews; iOverview++ )
00448 papoOvrBands[iOverview]->FlushCache();
00449
00450 pfnProgress( 1.0, NULL, pProgressData );
00451
00452 return CE_None;
00453 }