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 }