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 }