Main Page Class Hierarchy Compound List File List Compound Members File Members Related Pages

# overview.cpp

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