00001 /******************************************************************************
00002 * $Id: overview_cpp-source.html,v 1.13 2002/12/21 19:13:13 warmerda Exp $
00003 *
00004 * Project: GDAL Core
00005 * Purpose: Helper code to implement overview support in different drivers.
00006 * Author: Frank Warmerdam, warmerda@home.com
00007 *
00008 ******************************************************************************
00009 * Copyright (c) 2000, Frank Warmerdam
00010 *
00011 * Permission is hereby granted, free of charge, to any person obtaining a
00012 * copy of this software and associated documentation files (the "Software"),
00013 * to deal in the Software without restriction, including without limitation
00014 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
00015 * and/or sell copies of the Software, and to permit persons to whom the
00016 * Software is furnished to do so, subject to the following conditions:
00017 *
00018 * The above copyright notice and this permission notice shall be included
00019 * in all copies or substantial portions of the Software.
00020 *
00021 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
00022 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
00023 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
00024 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
00025 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
00026 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
00027 * DEALINGS IN THE SOFTWARE.
00028 ******************************************************************************
00029 *
00030 * $Log: overview_cpp-source.html,v $
00030 * Revision 1.13 2002/12/21 19:13:13 warmerda
00030 * updated
00030 *
00031 * Revision 1.10 2002/07/09 20:33:12 warmerda
00032 * expand tabs
00033 *
00034 * Revision 1.9 2002/04/16 17:49:29 warmerda
00035 * Ensure dfRatio is assigned a value.
00036 *
00037 * Revision 1.8 2001/07/18 04:04:30 warmerda
00038 * added CPL_CVSID
00039 *
00040 * Revision 1.7 2001/07/16 15:21:46 warmerda
00041 * added AVERAGE_MAGPHASE option for complex images
00042 *
00043 * Revision 1.6 2001/01/30 22:32:42 warmerda
00044 * added AVERAGE_MP (magnitude preserving averaging) overview resampling type
00045 *
00046 * Revision 1.5 2000/11/22 18:41:45 warmerda
00047 * fixed bug in complex overview generation
00048 *
00049 * Revision 1.4 2000/08/18 15:25:06 warmerda
00050 * added cascading overview regeneration to speed up averaged overviews
00051 *
00052 * Revision 1.3 2000/07/17 17:08:45 warmerda
00053 * added support for complex data
00054 *
00055 * Revision 1.2 2000/06/26 22:17:58 warmerda
00056 * added scaled progress support
00057 *
00058 * Revision 1.1 2000/04/21 21:54:05 warmerda
00059 * New
00060 *
00061 */
00062
00063 #include "gdal_priv.h"
00064
00065 CPL_CVSID("$Id: overview_cpp-source.html,v 1.13 2002/12/21 19:13:13 warmerda Exp $");
00066
00067 /************************************************************************/
00068 /* GDALDownsampleChunk32R() */
00069 /************************************************************************/
00070
00071 static CPLErr
00072 GDALDownsampleChunk32R( int nSrcWidth, int nSrcHeight,
00073 float * pafChunk, int nChunkYOff, int nChunkYSize,
00074 GDALRasterBand * poOverview,
00075 const char * pszResampling )
00076
00077 {
00078 int nDstYOff, nDstYOff2, nOXSize, nOYSize;
00079 float *pafDstScanline;
00080
00081 nOXSize = poOverview->GetXSize();
00082 nOYSize = poOverview->GetYSize();
00083
00084 pafDstScanline = (float *) CPLMalloc(nOXSize * sizeof(float));
00085
00086 /* -------------------------------------------------------------------- */
00087 /* Figure out the line to start writing to, and the first line */
00088 /* to not write to. In theory this approach should ensure that */
00089 /* every output line will be written if all input chunks are */
00090 /* processed. */
00091 /* -------------------------------------------------------------------- */
00092 nDstYOff = (int) (0.5 + (nChunkYOff/(double)nSrcHeight) * nOYSize);
00093 nDstYOff2 = (int)
00094 (0.5 + ((nChunkYOff+nChunkYSize)/(double)nSrcHeight) * nOYSize);
00095
00096 if( nChunkYOff + nChunkYSize == nSrcHeight )
00097 nDstYOff2 = nOYSize;
00098
00099 /* ==================================================================== */
00100 /* Loop over destination scanlines. */
00101 /* ==================================================================== */
00102 for( int iDstLine = nDstYOff; iDstLine < nDstYOff2; iDstLine++ )
00103 {
00104 float *pafSrcScanline;
00105 int nSrcYOff, nSrcYOff2, iDstPixel;
00106
00107 nSrcYOff = (int) (0.5 + (iDstLine/(double)nOYSize) * nSrcHeight);
00108 if( nSrcYOff < nChunkYOff )
00109 nSrcYOff = nChunkYOff;
00110
00111 nSrcYOff2 = (int) (0.5 + ((iDstLine+1)/(double)nOYSize) * nSrcHeight);
00112 if( nSrcYOff2 > nSrcHeight || iDstLine == nOYSize-1 )
00113 nSrcYOff2 = nSrcHeight;
00114 if( nSrcYOff2 > nChunkYOff + nChunkYSize )
00115 nSrcYOff2 = nChunkYOff + nChunkYSize;
00116
00117 pafSrcScanline = pafChunk + ((nSrcYOff-nChunkYOff) * nSrcWidth);
00118
00119 /* -------------------------------------------------------------------- */
00120 /* Loop over destination pixels */
00121 /* -------------------------------------------------------------------- */
00122 for( iDstPixel = 0; iDstPixel < nOXSize; iDstPixel++ )
00123 {
00124 int nSrcXOff, nSrcXOff2;
00125
00126 nSrcXOff = (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
00127 nSrcXOff2 = (int)
00128 (0.5 + ((iDstPixel+1)/(double)nOXSize) * nSrcWidth);
00129 if( nSrcXOff2 > nSrcWidth )
00130 nSrcXOff2 = nSrcWidth;
00131
00132 if( EQUALN(pszResampling,"NEAR",4) )
00133 {
00134 pafDstScanline[iDstPixel] = pafSrcScanline[nSrcXOff];
00135 }
00136 else if( EQUALN(pszResampling,"AVER",4) )
00137 {
00138 double dfTotal = 0.0;
00139 int nCount = 0, iX, iY;
00140
00141 for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
00142 {
00143 for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
00144 {
00145 dfTotal += pafSrcScanline[iX+(iY-nSrcYOff)*nSrcWidth];
00146 nCount++;
00147 }
00148 }
00149
00150 CPLAssert( nCount > 0 );
00151 if( nCount == 0 )
00152 {
00153 pafDstScanline[iDstPixel] = 0.0;
00154 }
00155 else
00156 pafDstScanline[iDstPixel] = dfTotal / nCount;
00157 }
00158 }
00159
00160 poOverview->RasterIO( GF_Write, 0, iDstLine, nOXSize, 1,
00161 pafDstScanline, nOXSize, 1, GDT_Float32,
00162 0, 0 );
00163 }
00164
00165 CPLFree( pafDstScanline );
00166
00167 return CE_None;
00168 }
00169
00170 /************************************************************************/
00171 /* GDALDownsampleChunkC32R() */
00172 /************************************************************************/
00173
00174 static CPLErr
00175 GDALDownsampleChunkC32R( int nSrcWidth, int nSrcHeight,
00176 float * pafChunk, int nChunkYOff, int nChunkYSize,
00177 GDALRasterBand * poOverview,
00178 const char * pszResampling )
00179
00180 {
00181 int nDstYOff, nDstYOff2, nOXSize, nOYSize;
00182 float *pafDstScanline;
00183
00184 nOXSize = poOverview->GetXSize();
00185 nOYSize = poOverview->GetYSize();
00186
00187 pafDstScanline = (float *) CPLMalloc(nOXSize * sizeof(float) * 2);
00188
00189 /* -------------------------------------------------------------------- */
00190 /* Figure out the line to start writing to, and the first line */
00191 /* to not write to. In theory this approach should ensure that */
00192 /* every output line will be written if all input chunks are */
00193 /* processed. */
00194 /* -------------------------------------------------------------------- */
00195 nDstYOff = (int) (0.5 + (nChunkYOff/(double)nSrcHeight) * nOYSize);
00196 nDstYOff2 = (int)
00197 (0.5 + ((nChunkYOff+nChunkYSize)/(double)nSrcHeight) * nOYSize);
00198
00199 if( nChunkYOff + nChunkYSize == nSrcHeight )
00200 nDstYOff2 = nOYSize;
00201
00202 /* ==================================================================== */
00203 /* Loop over destination scanlines. */
00204 /* ==================================================================== */
00205 for( int iDstLine = nDstYOff; iDstLine < nDstYOff2; iDstLine++ )
00206 {
00207 float *pafSrcScanline;
00208 int nSrcYOff, nSrcYOff2, iDstPixel;
00209
00210 nSrcYOff = (int) (0.5 + (iDstLine/(double)nOYSize) * nSrcHeight);
00211 if( nSrcYOff < nChunkYOff )
00212 nSrcYOff = nChunkYOff;
00213
00214 nSrcYOff2 = (int) (0.5 + ((iDstLine+1)/(double)nOYSize) * nSrcHeight);
00215 if( nSrcYOff2 > nSrcHeight || iDstLine == nOYSize-1 )
00216 nSrcYOff2 = nSrcHeight;
00217 if( nSrcYOff2 > nChunkYOff + nChunkYSize )
00218 nSrcYOff2 = nChunkYOff + nChunkYSize;
00219
00220 pafSrcScanline = pafChunk + ((nSrcYOff-nChunkYOff) * nSrcWidth) * 2;
00221
00222 /* -------------------------------------------------------------------- */
00223 /* Loop over destination pixels */
00224 /* -------------------------------------------------------------------- */
00225 for( iDstPixel = 0; iDstPixel < nOXSize; iDstPixel++ )
00226 {
00227 int nSrcXOff, nSrcXOff2;
00228
00229 nSrcXOff = (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
00230 nSrcXOff2 = (int)
00231 (0.5 + ((iDstPixel+1)/(double)nOXSize) * nSrcWidth);
00232 if( nSrcXOff2 > nSrcWidth )
00233 nSrcXOff2 = nSrcWidth;
00234
00235 if( EQUALN(pszResampling,"NEAR",4) )
00236 {
00237 pafDstScanline[iDstPixel*2] = pafSrcScanline[nSrcXOff*2];
00238 pafDstScanline[iDstPixel*2+1] = pafSrcScanline[nSrcXOff*2+1];
00239 }
00240 else if( EQUAL(pszResampling,"AVERAGE_MAGPHASE") )
00241 {
00242 double dfTotalR = 0.0, dfTotalI = 0.0, dfTotalM = 0.0;
00243 int nCount = 0, iX, iY;
00244
00245 for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
00246 {
00247 for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
00248 {
00249 double dfR, dfI;
00250
00251 dfR = pafSrcScanline[iX*2+(iY-nSrcYOff)*nSrcWidth*2];
00252 dfI = pafSrcScanline[iX*2+(iY-nSrcYOff)*nSrcWidth*2+1];
00253 dfTotalR += dfR;
00254 dfTotalI += dfI;
00255 dfTotalM += sqrt( dfR*dfR + dfI*dfI );
00256 nCount++;
00257 }
00258 }
00259
00260 CPLAssert( nCount > 0 );
00261 if( nCount == 0 )
00262 {
00263 pafDstScanline[iDstPixel*2] = 0.0;
00264 pafDstScanline[iDstPixel*2+1] = 0.0;
00265 }
00266 else
00267 {
00268 double dfM, dfDesiredM, dfRatio=1.0;
00269
00270 pafDstScanline[iDstPixel*2 ] = dfTotalR / nCount;
00271 pafDstScanline[iDstPixel*2+1] = dfTotalI / nCount;
00272
00273 dfM = sqrt(pafDstScanline[iDstPixel*2 ]*pafDstScanline[iDstPixel*2 ]
00274 + pafDstScanline[iDstPixel*2+1]*pafDstScanline[iDstPixel*2+1]);
00275 dfDesiredM = dfTotalM / nCount;
00276 if( dfM != 0.0 )
00277 dfRatio = dfDesiredM / dfM;
00278
00279 pafDstScanline[iDstPixel*2 ] *= dfRatio;
00280 pafDstScanline[iDstPixel*2+1] *= dfRatio;
00281 }
00282 }
00283 else if( EQUALN(pszResampling,"AVER",4) )
00284 {
00285 double dfTotalR = 0.0, dfTotalI = 0.0;
00286 int nCount = 0, iX, iY;
00287
00288 for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
00289 {
00290 for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
00291 {
00292 dfTotalR += pafSrcScanline[iX*2+(iY-nSrcYOff)*nSrcWidth*2];
00293 dfTotalI += pafSrcScanline[iX*2+(iY-nSrcYOff)*nSrcWidth*2+1];
00294 nCount++;
00295 }
00296 }
00297
00298 CPLAssert( nCount > 0 );
00299 if( nCount == 0 )
00300 {
00301 pafDstScanline[iDstPixel*2] = 0.0;
00302 pafDstScanline[iDstPixel*2+1] = 0.0;
00303 }
00304 else
00305 {
00306 pafDstScanline[iDstPixel*2 ] = dfTotalR / nCount;
00307 pafDstScanline[iDstPixel*2+1] = dfTotalI / nCount;
00308 }
00309 }
00310 }
00311
00312 poOverview->RasterIO( GF_Write, 0, iDstLine, nOXSize, 1,
00313 pafDstScanline, nOXSize, 1, GDT_CFloat32,
00314 0, 0 );
00315 }
00316
00317 CPLFree( pafDstScanline );
00318
00319 return CE_None;
00320 }
00321
00322 /************************************************************************/
00323 /* GDALRegenerateCascadingOverviews() */
00324 /* */
00325 /* Generate a list of overviews in order from largest to */
00326 /* smallest, computing each from the next larger. */
00327 /************************************************************************/
00328
00329 static CPLErr
00330 GDALRegenerateCascadingOverviews(
00331 GDALRasterBand *poSrcBand, int nOverviews, GDALRasterBand **papoOvrBands,
00332 const char * pszResampling,
00333 GDALProgressFunc pfnProgress, void * pProgressData )
00334
00335 {
00336 /* -------------------------------------------------------------------- */
00337 /* First, we must put the overviews in order from largest to */
00338 /* smallest. */
00339 /* -------------------------------------------------------------------- */
00340 int i, j;
00341
00342 for( i = 0; i < nOverviews-1; i++ )
00343 {
00344 for( j = 0; j < nOverviews - i - 1; j++ )
00345 {
00346
00347 if( papoOvrBands[j]->GetXSize()
00348 * (float) papoOvrBands[j]->GetYSize() <
00349 papoOvrBands[j+1]->GetXSize()
00350 * (float) papoOvrBands[j+1]->GetYSize() )
00351 {
00352 GDALRasterBand * poTempBand;
00353
00354 poTempBand = papoOvrBands[j];
00355 papoOvrBands[j] = papoOvrBands[j+1];
00356 papoOvrBands[j+1] = poTempBand;
00357 }
00358 }
00359 }
00360
00361 /* -------------------------------------------------------------------- */
00362 /* Count total pixels so we can prepare appropriate scaled */
00363 /* progress functions. */
00364 /* -------------------------------------------------------------------- */
00365 double dfTotalPixels = 0.0;
00366
00367 for( i = 0; i < nOverviews; i++ )
00368 {
00369 dfTotalPixels += papoOvrBands[i]->GetXSize()
00370 * (double) papoOvrBands[i]->GetYSize();
00371 }
00372
00373 /* -------------------------------------------------------------------- */
00374 /* Generate all the bands. */
00375 /* -------------------------------------------------------------------- */
00376 double dfPixelsProcessed = 0.0;
00377
00378 for( i = 0; i < nOverviews; i++ )
00379 {
00380 void *pScaledProgressData;
00381 double dfPixels;
00382 GDALRasterBand *poBaseBand;
00383 CPLErr eErr;
00384
00385 if( i == 0 )
00386 poBaseBand = poSrcBand;
00387 else
00388 poBaseBand = papoOvrBands[i-1];
00389
00390 dfPixels = papoOvrBands[i]->GetXSize()
00391 * (double) papoOvrBands[i]->GetYSize();
00392
00393 pScaledProgressData = GDALCreateScaledProgress(
00394 dfPixelsProcessed / dfTotalPixels,
00395 (dfPixelsProcessed + dfPixels) / dfTotalPixels,
00396 pfnProgress, pProgressData );
00397
00398 eErr = GDALRegenerateOverviews( poBaseBand, 1, papoOvrBands + i,
00399 pszResampling,
00400 GDALScaledProgress,
00401 pScaledProgressData );
00402 GDALDestroyScaledProgress( pScaledProgressData );
00403
00404 if( eErr != CE_None )
00405 return eErr;
00406
00407 dfPixelsProcessed += dfPixels;
00408 }
00409
00410 return CE_None;
00411 }
00412
00413 /************************************************************************/
00414 /* GDALRegenerateOverviews() */
00415 /************************************************************************/
00416 CPLErr
00417 GDALRegenerateOverviews( GDALRasterBand *poSrcBand,
00418 int nOverviews, GDALRasterBand **papoOvrBands,
00419 const char * pszResampling,
00420 GDALProgressFunc pfnProgress, void * pProgressData )
00421
00422 {
00423 int nFullResYChunk, nWidth;
00424 int nFRXBlockSize, nFRYBlockSize;
00425 GDALDataType eType;
00426
00427 /* -------------------------------------------------------------------- */
00428 /* If we are operating on multiple overviews, and using */
00429 /* averaging, lets do them in cascading order to reduce the */
00430 /* amount of computation. */
00431 /* -------------------------------------------------------------------- */
00432 if( EQUALN(pszResampling,"AVER",4) && nOverviews > 1 )
00433 return GDALRegenerateCascadingOverviews( poSrcBand,
00434 nOverviews, papoOvrBands,
00435 pszResampling,
00436 pfnProgress,
00437 pProgressData );
00438
00439 /* -------------------------------------------------------------------- */
00440 /* Setup one horizontal swath to read from the raw buffer. */
00441 /* -------------------------------------------------------------------- */
00442 float *pafChunk;
00443
00444 poSrcBand->GetBlockSize( &nFRXBlockSize, &nFRYBlockSize );
00445
00446 if( nFRYBlockSize < 4 || nFRYBlockSize > 256 )
00447 nFullResYChunk = 32;
00448 else
00449 nFullResYChunk = nFRYBlockSize;
00450
00451 if( GDALDataTypeIsComplex( poSrcBand->GetRasterDataType() ) )
00452 eType = GDT_CFloat32;
00453 else
00454 eType = GDT_Float32;
00455
00456 nWidth = poSrcBand->GetXSize();
00457 pafChunk = (float *)
00458 VSIMalloc((GDALGetDataTypeSize(eType)/8) * nFullResYChunk * nWidth );
00459
00460 if( pafChunk == NULL )
00461 {
00462 CPLError( CE_Failure, CPLE_OutOfMemory,
00463 "Out of memory in GDALRegenerateOverviews()." );
00464
00465 return CE_Failure;
00466 }
00467
00468 /* -------------------------------------------------------------------- */
00469 /* Loop over image operating on chunks. */
00470 /* -------------------------------------------------------------------- */
00471 int nChunkYOff = 0;
00472
00473 for( nChunkYOff = 0;
00474 nChunkYOff < poSrcBand->GetYSize();
00475 nChunkYOff += nFullResYChunk )
00476 {
00477 if( !pfnProgress( nChunkYOff / (double) poSrcBand->GetYSize(),
00478 NULL, pProgressData ) )
00479 {
00480 CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
00481 return CE_Failure;
00482 }
00483
00484 if( nFullResYChunk + nChunkYOff > poSrcBand->GetYSize() )
00485 nFullResYChunk = poSrcBand->GetYSize() - nChunkYOff;
00486
00487 /* read chunk */
00488 poSrcBand->RasterIO( GF_Read, 0, nChunkYOff, nWidth, nFullResYChunk,
00489 pafChunk, nWidth, nFullResYChunk, eType,
00490 0, 0 );
00491
00492 for( int iOverview = 0; iOverview < nOverviews; iOverview++ )
00493 {
00494 if( eType == GDT_Float32 )
00495 GDALDownsampleChunk32R(nWidth, poSrcBand->GetYSize(),
00496 pafChunk, nChunkYOff, nFullResYChunk,
00497 papoOvrBands[iOverview], pszResampling);
00498 else
00499 GDALDownsampleChunkC32R(nWidth, poSrcBand->GetYSize(),
00500 pafChunk, nChunkYOff, nFullResYChunk,
00501 papoOvrBands[iOverview], pszResampling);
00502 }
00503 }
00504
00505 VSIFree( pafChunk );
00506
00507 /* -------------------------------------------------------------------- */
00508 /* Renormalized overview mean / stddev if needed. */
00509 /* -------------------------------------------------------------------- */
00510 if( EQUAL(pszResampling,"AVERAGE_MP") )
00511 {
00512 GDALOverviewMagnitudeCorrection( (GDALRasterBandH) poSrcBand,
00513 nOverviews,
00514 (GDALRasterBandH *) papoOvrBands,
00515 GDALDummyProgress, NULL );
00516 }
00517
00518 /* -------------------------------------------------------------------- */
00519 /* It can be important to flush out data to overviews. */
00520 /* -------------------------------------------------------------------- */
00521 for( int iOverview = 0; iOverview < nOverviews; iOverview++ )
00522 papoOvrBands[iOverview]->FlushCache();
00523
00524 pfnProgress( 1.0, NULL, pProgressData );
00525
00526 return CE_None;
00527 }
00528
00529 /************************************************************************/
00530 /* GDALComputeBandStats() */
00531 /************************************************************************/
00532
00533 CPLErr
00534 GDALComputeBandStats( GDALRasterBandH hSrcBand,
00535 int nSampleStep,
00536 double *pdfMean, double *pdfStdDev,
00537 GDALProgressFunc pfnProgress,
00538 void *pProgressData )
00539
00540 {
00541 GDALRasterBand *poSrcBand = (GDALRasterBand *) hSrcBand;
00542 int iLine, nWidth, nHeight;
00543 GDALDataType eType = poSrcBand->GetRasterDataType();
00544 GDALDataType eWrkType;
00545 int bComplex;
00546 float *pafData;
00547 double dfSum=0.0, dfSum2=0.0;
00548 int nSamples = 0;
00549
00550 nWidth = poSrcBand->GetXSize();
00551 nHeight = poSrcBand->GetYSize();
00552
00553 if( nSampleStep >= nHeight )
00554 nSampleStep = 1;
00555
00556 bComplex = GDALDataTypeIsComplex(eType);
00557 if( bComplex )
00558 {
00559 pafData = (float *) CPLMalloc(nWidth * 2 * sizeof(float));
00560 eWrkType = GDT_CFloat32;
00561 }
00562 else
00563 {
00564 pafData = (float *) CPLMalloc(nWidth * sizeof(float));
00565 eWrkType = GDT_Float32;
00566 }
00567
00568 /* -------------------------------------------------------------------- */
00569 /* Loop over all sample lines. */
00570 /* -------------------------------------------------------------------- */
00571 for( iLine = 0; iLine < nHeight; iLine += nSampleStep )
00572 {
00573 int iPixel;
00574
00575 if( !pfnProgress( iLine / (double) nHeight,
00576 NULL, pProgressData ) )
00577 {
00578 CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
00579 CPLFree( pafData );
00580 return CE_Failure;
00581 }
00582
00583 poSrcBand->RasterIO( GF_Read, 0, iLine, nWidth, 1,
00584 pafData, nWidth, 1, eWrkType,
00585 0, 0 );
00586
00587 for( iPixel = 0; iPixel < nWidth; iPixel++ )
00588 {
00589 float fValue;
00590
00591 if( bComplex )
00592 {
00593 // Compute the magnitude of the complex value.
00594
00595 fValue = sqrt(pafData[iPixel*2 ] * pafData[iPixel*2 ]
00596 + pafData[iPixel*2+1] * pafData[iPixel*2+1]);
00597 }
00598 else
00599 {
00600 fValue = pafData[iPixel];
00601 }
00602
00603 dfSum += fValue;
00604 dfSum2 += fValue * fValue;
00605 }
00606
00607 nSamples += nWidth;
00608 }
00609
00610 if( !pfnProgress( 1.0, NULL, pProgressData ) )
00611 {
00612 CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
00613 CPLFree( pafData );
00614 return CE_Failure;
00615 }
00616
00617 /* -------------------------------------------------------------------- */
00618 /* Produce the result values. */
00619 /* -------------------------------------------------------------------- */
00620 if( pdfMean != NULL )
00621 *pdfMean = dfSum / nSamples;
00622
00623 if( pdfStdDev != NULL )
00624 {
00625 double dfMean = dfSum / nSamples;
00626
00627 *pdfStdDev = sqrt((dfSum2 / nSamples) - (dfMean * dfMean));
00628 }
00629
00630 CPLFree( pafData );
00631
00632 return CE_None;
00633 }
00634
00635 /************************************************************************/
00636 /* GDALOverviewMagnitudeCorrection() */
00637 /* */
00638 /* Correct the mean and standard deviation of the overviews of */
00639 /* the given band to match the base layer approximately. */
00640 /************************************************************************/
00641
00642 CPLErr
00643 GDALOverviewMagnitudeCorrection( GDALRasterBandH hBaseBand,
00644 int nOverviewCount,
00645 GDALRasterBandH *pahOverviews,
00646 GDALProgressFunc pfnProgress,
00647 void *pProgressData )
00648
00649 {
00650 CPLErr eErr;
00651 double dfOrigMean, dfOrigStdDev;
00652
00653 /* -------------------------------------------------------------------- */
00654 /* Compute mean/stddev for source raster. */
00655 /* -------------------------------------------------------------------- */
00656 eErr = GDALComputeBandStats( hBaseBand, 2, &dfOrigMean, &dfOrigStdDev,
00657 pfnProgress, pProgressData );
00658
00659 if( eErr != CE_None )
00660 return eErr;
00661
00662 /* -------------------------------------------------------------------- */
00663 /* Loop on overview bands. */
00664 /* -------------------------------------------------------------------- */
00665 int iOverview;
00666
00667 for( iOverview = 0; iOverview < nOverviewCount; iOverview++ )
00668 {
00669 GDALRasterBand *poOverview = (GDALRasterBand *)pahOverviews[iOverview];
00670 double dfOverviewMean, dfOverviewStdDev;
00671 double dfGain;
00672
00673 eErr = GDALComputeBandStats( pahOverviews[iOverview], 1,
00674 &dfOverviewMean, &dfOverviewStdDev,
00675 pfnProgress, pProgressData );
00676
00677 if( eErr != CE_None )
00678 return eErr;
00679
00680 if( dfOrigStdDev < 0.0001 )
00681 dfGain = 1.0;
00682 else
00683 dfGain = dfOrigStdDev / dfOverviewStdDev;
00684
00685 /* -------------------------------------------------------------------- */
00686 /* Apply gain and offset. */
00687 /* -------------------------------------------------------------------- */
00688 GDALDataType eWrkType, eType = poOverview->GetRasterDataType();
00689 int iLine, nWidth, nHeight, bComplex;
00690 float *pafData;
00691
00692 nWidth = poOverview->GetXSize();
00693 nHeight = poOverview->GetYSize();
00694
00695 bComplex = GDALDataTypeIsComplex(eType);
00696 if( bComplex )
00697 {
00698 pafData = (float *) CPLMalloc(nWidth * 2 * sizeof(float));
00699 eWrkType = GDT_CFloat32;
00700 }
00701 else
00702 {
00703 pafData = (float *) CPLMalloc(nWidth * sizeof(float));
00704 eWrkType = GDT_Float32;
00705 }
00706
00707 for( iLine = 0; iLine < nHeight; iLine++ )
00708 {
00709 int iPixel;
00710
00711 if( !pfnProgress( iLine / (double) nHeight,
00712 NULL, pProgressData ) )
00713 {
00714 CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
00715 CPLFree( pafData );
00716 return CE_Failure;
00717 }
00718
00719 poOverview->RasterIO( GF_Read, 0, iLine, nWidth, 1,
00720 pafData, nWidth, 1, eWrkType,
00721 0, 0 );
00722
00723 for( iPixel = 0; iPixel < nWidth; iPixel++ )
00724 {
00725 if( bComplex )
00726 {
00727 pafData[iPixel*2] *= dfGain;
00728 pafData[iPixel*2+1] *= dfGain;
00729 }
00730 else
00731 {
00732 pafData[iPixel] = (pafData[iPixel]-dfOverviewMean)*dfGain
00733 + dfOrigMean;
00734
00735 }
00736 }
00737
00738 poOverview->RasterIO( GF_Write, 0, iLine, nWidth, 1,
00739 pafData, nWidth, 1, eWrkType,
00740 0, 0 );
00741 }
00742
00743 if( !pfnProgress( 1.0, NULL, pProgressData ) )
00744 {
00745 CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
00746 CPLFree( pafData );
00747 return CE_Failure;
00748 }
00749
00750 CPLFree( pafData );
00751 }
00752
00753 return CE_None;
00754 }
1.2.3-20001105 written by Dimitri van Heesch,
© 1997-2000