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