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