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.3 2000/10/20 20:21:34 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.3  2000/10/20 20:21:34  warmerda
    00030  * *** empty log message ***
    00030  *
    00031  * Revision 1.4  2000/08/18 15:25:06  warmerda
    00032  * added cascading overview regeneration to speed up averaged overviews
    00033  *
    00034  * Revision 1.3  2000/07/17 17:08:45  warmerda
    00035  * added support for complex data
    00036  *
    00037  * Revision 1.2  2000/06/26 22:17:58  warmerda
    00038  * added scaled progress support
    00039  *
    00040  * Revision 1.1  2000/04/21 21:54:05  warmerda
    00041  * New
    00042  *
    00043  */
    00044
    00045 #include "gdal_priv.h"
    00046
    00047 /************************************************************************/
    00048 /*                       GDALDownsampleChunk32R()                       */
    00049 /************************************************************************/
    00050
    00051 static CPLErr
    00052 GDALDownsampleChunk32R( int nSrcWidth, int nSrcHeight,
    00053                         float * pafChunk, int nChunkYOff, int nChunkYSize,
    00054                         GDALRasterBand * poOverview,
    00055                         const char * pszResampling )
    00056
    00057 {
    00058     int      nDstYOff, nDstYOff2, nOXSize, nOYSize;
    00059     float    *pafDstScanline;
    00060
    00061     nOXSize = poOverview->GetXSize();
    00062     nOYSize = poOverview->GetYSize();
    00063
    00064     pafDstScanline = (float *) CPLMalloc(nOXSize * sizeof(float));
    00065
    00066 /* -------------------------------------------------------------------- */
    00067 /*      Figure out the line to start writing to, and the first line     */
    00068 /*      to not write to.  In theory this approach should ensure that    */
    00069 /*      every output line will be written if all input chunks are       */
    00070 /*      processed.                                                      */
    00071 /* -------------------------------------------------------------------- */
    00072     nDstYOff = (int) (0.5 + (nChunkYOff/(double)nSrcHeight) * nOYSize);
    00073     nDstYOff2 = (int)
    00074         (0.5 + ((nChunkYOff+nChunkYSize)/(double)nSrcHeight) * nOYSize);
    00075
    00076     if( nChunkYOff + nChunkYSize == nSrcHeight )
    00077         nDstYOff2 = nOYSize;
    00078
    00079 /* ==================================================================== */
    00080 /*      Loop over destination scanlines.                                */
    00081 /* ==================================================================== */
    00082     for( int iDstLine = nDstYOff; iDstLine < nDstYOff2; iDstLine++ )
    00083     {
    00084         float *pafSrcScanline;
    00085         int   nSrcYOff, nSrcYOff2, iDstPixel;
    00086
    00087         nSrcYOff = (int) (0.5 + (iDstLine/(double)nOYSize) * nSrcHeight);
    00088         if( nSrcYOff < nChunkYOff )
    00089             nSrcYOff = nChunkYOff;
    00090
    00091         nSrcYOff2 = (int) (0.5 + ((iDstLine+1)/(double)nOYSize) * nSrcHeight);
    00092         if( nSrcYOff2 > nSrcHeight || iDstLine == nOYSize-1 )
    00093             nSrcYOff2 = nSrcHeight;
    00094         if( nSrcYOff2 > nChunkYOff + nChunkYSize )
    00095             nSrcYOff2 = nChunkYOff + nChunkYSize;
    00096
    00097         pafSrcScanline = pafChunk + ((nSrcYOff-nChunkYOff) * nSrcWidth);
    00098
    00099 /* -------------------------------------------------------------------- */
    00100 /*      Loop over destination pixels                                    */
    00101 /* -------------------------------------------------------------------- */
    00102         for( iDstPixel = 0; iDstPixel < nOXSize; iDstPixel++ )
    00103         {
    00104             int   nSrcXOff, nSrcXOff2;
    00105
    00106             nSrcXOff = (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
    00107             nSrcXOff2 = (int)
    00108                 (0.5 + ((iDstPixel+1)/(double)nOXSize) * nSrcWidth);
    00109             if( nSrcXOff2 > nSrcWidth )
    00110                 nSrcXOff2 = nSrcWidth;
    00111
    00112             if( EQUALN(pszResampling,"NEAR",4) )
    00113             {
    00114                 pafDstScanline[iDstPixel] = pafSrcScanline[nSrcXOff];
    00115             }
    00116             else if( EQUALN(pszResampling,"AVER",4) )
    00117             {
    00118                 double dfTotal = 0.0;
    00119                 int    nCount = 0, iX, iY;
    00120
    00121                 for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
    00122                 {
    00123                     for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
    00124                     {
    00125                         dfTotal += pafSrcScanline[iX+(iY-nSrcYOff)*nSrcWidth];
    00126                         nCount++;
    00127                     }
    00128                 }
    00129
    00130                 CPLAssert( nCount > 0 );
    00131                 if( nCount == 0 )
    00132                 {
    00133                     pafDstScanline[iDstPixel] = 0.0;
    00134                 }
    00135                 else
    00136                     pafDstScanline[iDstPixel] = dfTotal / nCount;
    00137             }
    00138         }
    00139
    00140         poOverview->RasterIO( GF_Write, 0, iDstLine, nOXSize, 1,
    00141                               pafDstScanline, nOXSize, 1, GDT_Float32,
    00142                               0, 0 );
    00143     }
    00144
    00145     CPLFree( pafDstScanline );
    00146
    00147     return CE_None;
    00148 }
    00149
    00150 /************************************************************************/
    00151 /*                       GDALDownsampleChunkC32R()                      */
    00152 /************************************************************************/
    00153
    00154 static CPLErr
    00155 GDALDownsampleChunkC32R( int nSrcWidth, int nSrcHeight,
    00156                          float * pafChunk, int nChunkYOff, int nChunkYSize,
    00157                          GDALRasterBand * poOverview,
    00158                          const char * pszResampling )
    00159
    00160 {
    00161     int      nDstYOff, nDstYOff2, nOXSize, nOYSize;
    00162     float    *pafDstScanline;
    00163
    00164     nOXSize = poOverview->GetXSize();
    00165     nOYSize = poOverview->GetYSize();
    00166
    00167     pafDstScanline = (float *) CPLMalloc(nOXSize * sizeof(float) * 2);
    00168
    00169 /* -------------------------------------------------------------------- */
    00170 /*      Figure out the line to start writing to, and the first line     */
    00171 /*      to not write to.  In theory this approach should ensure that    */
    00172 /*      every output line will be written if all input chunks are       */
    00173 /*      processed.                                                      */
    00174 /* -------------------------------------------------------------------- */
    00175     nDstYOff = (int) (0.5 + (nChunkYOff/(double)nSrcHeight) * nOYSize);
    00176     nDstYOff2 = (int)
    00177         (0.5 + ((nChunkYOff+nChunkYSize)/(double)nSrcHeight) * nOYSize);
    00178
    00179     if( nChunkYOff + nChunkYSize == nSrcHeight )
    00180         nDstYOff2 = nOYSize;
    00181
    00182 /* ==================================================================== */
    00183 /*      Loop over destination scanlines.                                */
    00184 /* ==================================================================== */
    00185     for( int iDstLine = nDstYOff; iDstLine < nDstYOff2; iDstLine++ )
    00186     {
    00187         float *pafSrcScanline;
    00188         int   nSrcYOff, nSrcYOff2, iDstPixel;
    00189
    00190         nSrcYOff = (int) (0.5 + (iDstLine/(double)nOYSize) * nSrcHeight);
    00191         if( nSrcYOff < nChunkYOff )
    00192             nSrcYOff = nChunkYOff;
    00193
    00194         nSrcYOff2 = (int) (0.5 + ((iDstLine+1)/(double)nOYSize) * nSrcHeight);
    00195         if( nSrcYOff2 > nSrcHeight || iDstLine == nOYSize-1 )
    00196             nSrcYOff2 = nSrcHeight;
    00197         if( nSrcYOff2 > nChunkYOff + nChunkYSize )
    00198             nSrcYOff2 = nChunkYOff + nChunkYSize;
    00199
    00200         pafSrcScanline = pafChunk + ((nSrcYOff-nChunkYOff) * nSrcWidth);
    00201
    00202 /* -------------------------------------------------------------------- */
    00203 /*      Loop over destination pixels                                    */
    00204 /* -------------------------------------------------------------------- */
    00205         for( iDstPixel = 0; iDstPixel < nOXSize; iDstPixel++ )
    00206         {
    00207             int   nSrcXOff, nSrcXOff2;
    00208
    00209             nSrcXOff = (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
    00210             nSrcXOff2 = (int)
    00211                 (0.5 + ((iDstPixel+1)/(double)nOXSize) * nSrcWidth);
    00212             if( nSrcXOff2 > nSrcWidth )
    00213                 nSrcXOff2 = nSrcWidth;
    00214
    00215             if( EQUALN(pszResampling,"NEAR",4) )
    00216             {
    00217                 pafDstScanline[iDstPixel*2] = pafSrcScanline[nSrcXOff*2];
    00218                 pafDstScanline[iDstPixel*2+1] = pafSrcScanline[nSrcXOff*2+1];
    00219             }
    00220             else if( EQUALN(pszResampling,"AVER",4) )
    00221             {
    00222                 double dfTotalR = 0.0, dfTotalI = 0.0;
    00223                 int    nCount = 0, iX, iY;
    00224
    00225                 for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
    00226                 {
    00227                     for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
    00228                     {
    00229                         dfTotalR += pafSrcScanline[iX*2+(iY-nSrcYOff)*nSrcWidth*2];
    00230                         dfTotalI += pafSrcScanline[iX*2+(iY-nSrcYOff)*nSrcWidth*2+1];
    00231                         nCount++;
    00232                     }
    00233                 }
    00234
    00235                 CPLAssert( nCount > 0 );
    00236                 if( nCount == 0 )
    00237                 {
    00238                     pafDstScanline[iDstPixel*2] = 0.0;
    00239                     pafDstScanline[iDstPixel*2+1] = 0.0;
    00240                 }
    00241                 else
    00242                 {
    00243                     pafDstScanline[iDstPixel*2  ] = dfTotalR / nCount;
    00244                     pafDstScanline[iDstPixel*2+1] = dfTotalI / nCount;
    00245                 }
    00246             }
    00247         }
    00248
    00249         poOverview->RasterIO( GF_Write, 0, iDstLine, nOXSize, 1,
    00250                               pafDstScanline, nOXSize, 1, GDT_CFloat32,
    00251                               0, 0 );
    00252     }
    00253
    00254     CPLFree( pafDstScanline );
    00255
    00256     return CE_None;
    00257 }
    00258
    00259 /************************************************************************/
    00260 /*                  GDALRegenerateCascadingOverviews()                  */
    00261 /*                                                                      */
    00262 /*      Generate a list of overviews in order from largest to           */
    00263 /*      smallest, computing each from the next larger.                  */
    00264 /************************************************************************/
    00265
    00266 static CPLErr
    00267 GDALRegenerateCascadingOverviews(
    00268     GDALRasterBand *poSrcBand, int nOverviews, GDALRasterBand **papoOvrBands,
    00269     const char * pszResampling,
    00270     GDALProgressFunc pfnProgress, void * pProgressData )
    00271
    00272 {
    00273 /* -------------------------------------------------------------------- */
    00274 /*      First, we must put the overviews in order from largest to       */
    00275 /*      smallest.                                                       */
    00276 /* -------------------------------------------------------------------- */
    00277     int   i, j;
    00278
    00279     for( i = 0; i < nOverviews-1; i++ )
    00280     {
    00281         for( j = 0; j < nOverviews - i - 1; j++ )
    00282         {
    00283
    00284             if( papoOvrBands[j]->GetXSize()
    00285                 * (float) papoOvrBands[j]->GetYSize() <
    00286                 papoOvrBands[j+1]->GetXSize()
    00287                 * (float) papoOvrBands[j+1]->GetYSize() )
    00288             {
    00289                 GDALRasterBand * poTempBand;
    00290
    00291                 poTempBand = papoOvrBands[j];
    00292                 papoOvrBands[j] = papoOvrBands[j+1];
    00293                 papoOvrBands[j+1] = poTempBand;
    00294             }
    00295         }
    00296     }
    00297
    00298 /* -------------------------------------------------------------------- */
    00299 /*      Count total pixels so we can prepare appropriate scaled         */
    00300 /*      progress functions.                                             */
    00301 /* -------------------------------------------------------------------- */
    00302     double       dfTotalPixels = 0.0;
    00303
    00304     for( i = 0; i < nOverviews; i++ )
    00305     {
    00306         dfTotalPixels += papoOvrBands[i]->GetXSize()
    00307             * (double) papoOvrBands[i]->GetYSize();
    00308     }
    00309
    00310 /* -------------------------------------------------------------------- */
    00311 /*      Generate all the bands.                                         */
    00312 /* -------------------------------------------------------------------- */
    00313     double      dfPixelsProcessed = 0.0;
    00314
    00315     for( i = 0; i < nOverviews; i++ )
    00316     {
    00317         void    *pScaledProgressData;
    00318         double  dfPixels;
    00319         GDALRasterBand *poBaseBand;
    00320         CPLErr  eErr;
    00321
    00322         if( i == 0 )
    00323             poBaseBand = poSrcBand;
    00324         else
    00325             poBaseBand = papoOvrBands[i-1];
    00326
    00327         dfPixels = papoOvrBands[i]->GetXSize()
    00328             * (double) papoOvrBands[i]->GetYSize();
    00329
    00330         pScaledProgressData = GDALCreateScaledProgress(
    00331             dfPixelsProcessed / dfTotalPixels,
    00332             (dfPixelsProcessed + dfPixels) / dfTotalPixels,
    00333             pfnProgress, pProgressData );
    00334
    00335         eErr = GDALRegenerateOverviews( poBaseBand, 1, papoOvrBands + i,
    00336                                         pszResampling,
    00337                                         GDALScaledProgress,
    00338                                         pScaledProgressData );
    00339         GDALDestroyScaledProgress( pScaledProgressData );
    00340
    00341         if( eErr != CE_None )
    00342             return eErr;
    00343
    00344         dfPixelsProcessed += dfPixels;
    00345     }
    00346
    00347     return CE_None;
    00348 }
    00349
    00350 /************************************************************************/
    00351 /*                      GDALRegenerateOverviews()                       */
    00352 /************************************************************************/
    00353 CPLErr
    00354 GDALRegenerateOverviews( GDALRasterBand *poSrcBand,
    00355                          int nOverviews, GDALRasterBand **papoOvrBands,
    00356                          const char * pszResampling,
    00357                          GDALProgressFunc pfnProgress, void * pProgressData )
    00358
    00359 {
    00360     int    nFullResYChunk, nWidth;
    00361     int    nFRXBlockSize, nFRYBlockSize;
    00362     GDALDataType eType;
    00363
    00364 /* -------------------------------------------------------------------- */
    00365 /*      If we are operating on multiple overviews, and using            */
    00366 /*      averaging, lets do them in cascading order to reduce the        */
    00367 /*      amount of computation.                                          */
    00368 /* -------------------------------------------------------------------- */
    00369     if( EQUALN(pszResampling,"AVER",4) && nOverviews > 1 )
    00370         return GDALRegenerateCascadingOverviews( poSrcBand,
    00371                                                  nOverviews, papoOvrBands,
    00372                                                  pszResampling,
    00373                                                  pfnProgress,
    00374                                                  pProgressData );
    00375
    00376 /* -------------------------------------------------------------------- */
    00377 /*      Setup one horizontal swath to read from the raw buffer.         */
    00378 /* -------------------------------------------------------------------- */
    00379     float *pafChunk;
    00380
    00381     poSrcBand->GetBlockSize( &nFRXBlockSize, &nFRYBlockSize );
    00382
    00383     if( nFRYBlockSize < 4 || nFRYBlockSize > 256 )
    00384         nFullResYChunk = 32;
    00385     else
    00386         nFullResYChunk = nFRYBlockSize;
    00387
    00388     if( GDALDataTypeIsComplex( poSrcBand->GetRasterDataType() ) )
    00389         eType = GDT_CFloat32;
    00390     else
    00391         eType = GDT_Float32;
    00392
    00393     nWidth = poSrcBand->GetXSize();
    00394     pafChunk = (float *)
    00395         VSIMalloc((GDALGetDataTypeSize(eType)/8) * nFullResYChunk * nWidth );
    00396
    00397     if( pafChunk == NULL )
    00398     {
    00399         CPLError( CE_Failure, CPLE_OutOfMemory,
    00400                   "Out of memory in GDALRegenerateOverviews()." );
    00401
    00402         return CE_Failure;
    00403     }
    00404
    00405 /* -------------------------------------------------------------------- */
    00406 /*      Loop over image operating on chunks.                            */
    00407 /* -------------------------------------------------------------------- */
    00408     int  nChunkYOff = 0;
    00409
    00410     for( nChunkYOff = 0;
    00411          nChunkYOff < poSrcBand->GetYSize();
    00412          nChunkYOff += nFullResYChunk )
    00413     {
    00414         if( !pfnProgress( nChunkYOff / (double) poSrcBand->GetYSize(),
    00415                           NULL, pProgressData ) )
    00416         {
    00417             CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
    00418             return CE_Failure;
    00419         }
    00420
    00421         if( nFullResYChunk + nChunkYOff > poSrcBand->GetYSize() )
    00422             nFullResYChunk = poSrcBand->GetYSize() - nChunkYOff;
    00423
    00424         /* read chunk */
    00425         poSrcBand->RasterIO( GF_Read, 0, nChunkYOff, nWidth, nFullResYChunk,
    00426                              pafChunk, nWidth, nFullResYChunk, eType,
    00427                              0, 0 );
    00428
    00429         for( int iOverview = 0; iOverview < nOverviews; iOverview++ )
    00430         {
    00431             if( eType == GDT_Float32 )
    00432                 GDALDownsampleChunk32R(nWidth, poSrcBand->GetYSize(),
    00433                                        pafChunk, nChunkYOff, nFullResYChunk,
    00434                                        papoOvrBands[iOverview], pszResampling);
    00435             else
    00436                 GDALDownsampleChunkC32R(nWidth, poSrcBand->GetYSize(),
    00437                                        pafChunk, nChunkYOff, nFullResYChunk,
    00438                                        papoOvrBands[iOverview], pszResampling);
    00439         }
    00440     }
    00441
    00442     VSIFree( pafChunk );
    00443
    00444 /* -------------------------------------------------------------------- */
    00445 /*      It can be important to flush out data to overviews.             */
    00446 /* -------------------------------------------------------------------- */
    00447     for( int iOverview = 0; iOverview < nOverviews; iOverview++ )
    00448         papoOvrBands[iOverview]->FlushCache();
    00449
    00450     pfnProgress( 1.0, NULL, pProgressData );
    00451
    00452     return CE_None;
    00453 }