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 }

Generated at Sat Dec 21 14:02:00 2002 for GDAL by doxygen1.2.3-20001105 written by Dimitri van Heesch, © 1997-2000