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.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 }