Index: apps/gdalwarp.cpp =================================================================== --- apps/gdalwarp.cpp (revision 25854) +++ apps/gdalwarp.cpp (working copy) @@ -281,7 +281,7 @@ " srcfile* dstfile\n" "\n" "Available resampling methods:\n" - " near (default), bilinear, cubic, cubicspline, lanczos.\n" ); + " near (default), bilinear, cubic, cubicspline, lanczos, average, mode.\n" ); if( pszErrorMsg != NULL ) fprintf(stderr, "\nFAILURE: %s\n", pszErrorMsg); @@ -617,6 +617,10 @@ eResampleAlg = GRA_CubicSpline; else if ( EQUAL(argv[i], "lanczos") ) eResampleAlg = GRA_Lanczos; + else if ( EQUAL(argv[i], "average") ) + eResampleAlg = GRA_Average; + else if ( EQUAL(argv[i], "mode") ) + eResampleAlg = GRA_Mode; else { Usage(CPLSPrintf( "Unknown resampling method: \"%s\".", argv[i] )); Index: alg/gdalwarpkernel.cpp =================================================================== --- alg/gdalwarpkernel.cpp (revision 25854) +++ alg/gdalwarpkernel.cpp (working copy) @@ -44,7 +44,9 @@ 1, // Bilinear 2, // Cubic Convolution 2, // Cubic B-Spline - 3 // Lanczos windowed sinc + 3, // Lanczos windowed sinc + 0, // Average + 0, // Mode }; /* Used in gdalwarpoperation.cpp */ @@ -70,6 +72,7 @@ static CPLErr GWKNearestShort( GDALWarpKernel *poWK ); static CPLErr GWKNearestNoMasksFloat( GDALWarpKernel *poWK ); static CPLErr GWKNearestFloat( GDALWarpKernel *poWK ); +static CPLErr GWKAverageOrMode( GDALWarpKernel * ); /************************************************************************/ /* GWKJobStruct */ @@ -936,6 +939,12 @@ && eResample == GRA_NearestNeighbour ) return GWKNearestFloat( this ); + if( eResample == GRA_Average ) + return GWKAverageOrMode( this ); + + if( eResample == GRA_Mode ) + return GWKAverageOrMode( this ); + return GWKGeneralCase( this ); } @@ -4948,3 +4957,363 @@ CPLFree( pabSuccess ); } +/************************************************************************/ +/* GWKAverageOrMode() */ +/* */ +/************************************************************************/ + +static void GWKAverageOrModeThread(void* pData); + +static CPLErr GWKAverageOrMode( GDALWarpKernel *poWK ) +{ + return GWKRun( poWK, "GWKAverageOrMode", GWKAverageOrModeThread ); +} + +// overall logic based on GWKGeneralCaseThread() +static void GWKAverageOrModeThread( void* pData) +{ + GWKJobStruct* psJob = (GWKJobStruct*) pData; + GDALWarpKernel *poWK = psJob->poWK; + int iYMin = psJob->iYMin; + int iYMax = psJob->iYMax; + + int iDstY, iDstX, iSrcX, iSrcY, iDstOffset; + int nDstXSize = poWK->nDstXSize; + int nSrcXSize = poWK->nSrcXSize, nSrcYSize = poWK->nSrcYSize; + +/* -------------------------------------------------------------------- */ +/* Find out which algorithm to use (small optim.) */ +/* -------------------------------------------------------------------- */ + int nAlgo = 0; + if ( poWK->eResample == GRA_Average ) + { + nAlgo = 1; + } + else if( poWK->eResample == GRA_Mode ) + { + if ( poWK->eWorkingDataType != GDT_Byte ) // || nEntryCount > 256) + nAlgo = 2; + else + nAlgo = 3; + } + else + { + // other resample algorithms not permitted here + CPLDebug( "GDAL", "GDALWarpKernel():GWKAverageOrModeThread() ERROR, illegal resample" ); + return; + } + CPLDebug( "GDAL", "GDALWarpKernel():GWKAverageOrModeThread() using algo %d", nAlgo ); + +/* -------------------------------------------------------------------- */ +/* Allocate x,y,z coordinate arrays for transformation ... two */ +/* scanlines worth of positions. */ +/* -------------------------------------------------------------------- */ + double *padfX, *padfY, *padfZ; + double *padfX2, *padfY2, *padfZ2; + int *pabSuccess; + + padfX = (double *) CPLMalloc(sizeof(double) * nDstXSize); + padfY = (double *) CPLMalloc(sizeof(double) * nDstXSize); + padfZ = (double *) CPLMalloc(sizeof(double) * nDstXSize); + padfX2 = (double *) CPLMalloc(sizeof(double) * nDstXSize); + padfY2 = (double *) CPLMalloc(sizeof(double) * nDstXSize); + padfZ2 = (double *) CPLMalloc(sizeof(double) * nDstXSize); + pabSuccess = (int *) CPLMalloc(sizeof(int) * nDstXSize); + +/* ==================================================================== */ +/* Loop over output lines. */ +/* ==================================================================== */ + for( iDstY = iYMin; iDstY < iYMax; iDstY++ ) + { + +/* -------------------------------------------------------------------- */ +/* Setup points to transform to source image space. */ +/* -------------------------------------------------------------------- */ + for( iDstX = 0; iDstX < nDstXSize; iDstX++ ) + { + padfX[iDstX] = iDstX + poWK->nDstXOff; + padfY[iDstX] = iDstY + poWK->nDstYOff; + padfZ[iDstX] = 0.0; + padfX2[iDstX] = iDstX + 1.0 + poWK->nDstXOff; + padfY2[iDstX] = iDstY + 1.0 + poWK->nDstYOff; + padfZ2[iDstX] = 0.0; + } + +/* -------------------------------------------------------------------- */ +/* Transform the points from destination pixel/line coordinates */ +/* to source pixel/line coordinates. */ +/* -------------------------------------------------------------------- */ + poWK->pfnTransformer( psJob->pTransformerArg, TRUE, nDstXSize, + padfX, padfY, padfZ, pabSuccess ); + poWK->pfnTransformer( psJob->pTransformerArg, TRUE, nDstXSize, + padfX2, padfY2, padfZ2, pabSuccess ); + +/* ==================================================================== */ +/* Loop over pixels in output scanline. */ +/* ==================================================================== */ + for( iDstX = 0; iDstX < nDstXSize; iDstX++ ) + { + int iSrcOffset = 0; + +/* -------------------------------------------------------------------- */ +/* Do not try to apply transparent/invalid source pixels to the */ +/* destination. This currently ignores the multi-pixel input */ +/* of bilinear and cubic resamples. */ +/* -------------------------------------------------------------------- */ + double dfDensity = 1.0; + + // TODO - is this needed? + + // if( poWK->pafUnifiedSrcDensity != NULL ) + // { + // dfDensity = poWK->pafUnifiedSrcDensity[iSrcOffset]; + // if( dfDensity < 0.00001 ) + // continue; + // } + + // if( poWK->panUnifiedSrcValid != NULL + // && !(poWK->panUnifiedSrcValid[iSrcOffset>>5] + // & (0x01 << (iSrcOffset & 0x1f))) ) + // continue; + +/* ==================================================================== */ +/* Loop processing each band. */ +/* ==================================================================== */ + int iBand; + int bHasFoundDensity = FALSE; + + iDstOffset = iDstX + iDstY * nDstXSize; + for( iBand = 0; iBand < poWK->nBands; iBand++ ) + { + double dfBandDensity = 0.0; + double dfValueReal = 0.0; + double dfValueImag = 0.0; + double dfValueRealTmp = 0.0; + double dfValueImagTmp = 0.0; + +/* -------------------------------------------------------------------- */ +/* Collect the source value. */ +/* -------------------------------------------------------------------- */ + + double dfTotal = 0; + int nCount = 0; // count of pixels used to compute average/mode + int nCount2 = 0; // count of all pixels sampled, including nodata + int iSrcXMin, iSrcXMax,iSrcYMin,iSrcYMax; + + // compute corners in source crs + iSrcXMin = MAX( ((int) floor((padfX[iDstX] + 1e-10))) - poWK->nSrcXOff, 0 ); + iSrcXMax = MIN( ((int) ceil((padfX2[iDstX] + 1e-10))) - poWK->nSrcXOff, nSrcXSize ); + iSrcYMin = MAX( ((int) floor((padfY[iDstX] + 1e-10))) - poWK->nSrcYOff, 0 ); + iSrcYMax = MIN( ((int) ceil((padfY2[iDstX] + 1e-10))) - poWK->nSrcYOff, nSrcYSize ); + + // loop over source lines and pixels - 3 possible algorithms + + if ( nAlgo == 1 ) // poWK->eResample == GRA_Average + { + // this code adapted from GDALDownsampleChunk32R_AverageT() in gcore/overview.cpp + for( iSrcY = iSrcYMin; iSrcY < iSrcYMax; iSrcY++ ) + { + for( iSrcX = iSrcXMin; iSrcX < iSrcXMax; iSrcX++ ) + { + int iSrcOffset = iSrcX + iSrcY * nSrcXSize; + + if( poWK->panUnifiedSrcValid != NULL + && !(poWK->panUnifiedSrcValid[iSrcOffset>>5] + & (0x01 << (iSrcOffset & 0x1f))) ) + { + continue; + } + + nCount2++; + if ( GWKGetPixelValue( poWK, iBand, iSrcOffset, + &dfBandDensity, &dfValueRealTmp, &dfValueImagTmp ) && dfBandDensity > 0.0000000001 ) + { + // TODO - is this needed? + // if (pabyChunkNodataMask == NULL || + // pabyChunkNodataMask[iX + iY *nChunkXSize]) + dfTotal += dfValueRealTmp; + nCount++; + } + } + } + + if ( nCount > 0 ) + { + dfValueReal = dfTotal / nCount; + dfBandDensity = 1; + bHasFoundDensity = TRUE; + } + + } // GRA_Average + + else if ( nAlgo == 2 || nAlgo == 3 ) // poWK->eResample == GRA_Mode + { + // this code adapted from GDALDownsampleChunk32R_Mode() in gcore/overview.cpp + + // TODO add algorithm for Int16 / Int32 data type that uses less memory and cpu... + + if ( nAlgo == 2 ) // poWK->eWorkingDataType != GDT_Byte ) // || nEntryCount > 256) + { + /* I'm not sure how much sense it makes to run a majority + filter on floating point data, but here it is for the sake + of compatability. It won't look right on RGB images by the + nature of the filter. */ + int nNumPx = nSrcXSize * nSrcYSize; + int iMaxInd = 0, iMaxVal = -1; + + int nMaxNumPx = 0; + float* pafVals = NULL; + int* panSums = NULL; + + if (nNumPx > nMaxNumPx) + { + pafVals = (float*) CPLRealloc(pafVals, nNumPx * sizeof(float)); + panSums = (int*) CPLRealloc(panSums, nNumPx * sizeof(int)); + nMaxNumPx = nNumPx; + } + + for( iSrcY = iSrcYMin; iSrcY < iSrcYMax; iSrcY++ ) + { + for( iSrcX = iSrcXMin; iSrcX < iSrcXMax; iSrcX++ ) + { + iSrcOffset = iSrcX + iSrcY * nSrcXSize; + + if( poWK->panUnifiedSrcValid != NULL + && !(poWK->panUnifiedSrcValid[iSrcOffset>>5] + & (0x01 << (iSrcOffset & 0x1f))) ) + continue; + + nCount2++; + if ( GWKGetPixelValue( poWK, iBand, iSrcOffset, + &dfBandDensity, &dfValueRealTmp, &dfValueImagTmp ) && dfBandDensity > 0.0000000001 ) + { + float fVal = dfValueRealTmp; + int i; + + //Check array for existing entry + for( i = 0; i < iMaxInd; ++i ) + if( pafVals[i] == fVal + && ++panSums[i] > panSums[iMaxVal] ) + { + iMaxVal = i; + break; + } + + //Add to arr if entry not already there + if( i == iMaxInd ) + { + pafVals[iMaxInd] = fVal; + panSums[iMaxInd] = 1; + + if( iMaxVal < 0 ) + iMaxVal = iMaxInd; + + ++iMaxInd; + } + } + } + } + + if( iMaxVal != -1 ) + { + dfValueReal = pafVals[iMaxVal]; + dfBandDensity = 1; + bHasFoundDensity = TRUE; + } + } + + else // nAlgo == 3 / poWK->eWorkingDataType == GDT_Byte + { + /* So we go here for a paletted or non-paletted byte band */ + /* The input values are then between 0 and 255 */ + int anVals[256], nMaxVal = 0, iMaxInd = -1; + int nCount = 0; + memset(anVals, 0, 256*sizeof(int)); + + for( iSrcY = iSrcYMin; iSrcY < iSrcYMax; iSrcY++ ) + { + for( iSrcX = iSrcXMin; iSrcX < iSrcXMax; iSrcX++ ) + { + iSrcOffset = iSrcX + iSrcY * nSrcXSize; + + if( poWK->panUnifiedSrcValid != NULL + && !(poWK->panUnifiedSrcValid[iSrcOffset>>5] + & (0x01 << (iSrcOffset & 0x1f))) ) + continue; + + nCount2++; + if ( GWKGetPixelValue( poWK, iBand, iSrcOffset, + &dfBandDensity, &dfValueRealTmp, &dfValueImagTmp ) && dfBandDensity > 0.0000000001 ) + { + nCount++; + int nVal = (int) dfValueRealTmp; + if ( ++anVals[nVal] > nMaxVal) + { + //Sum the density + //Is it the most common value so far? + iMaxInd = nVal; + nMaxVal = anVals[nVal]; + } + } + } + } + + if( iMaxInd != -1 ) + { + dfValueReal = (float)iMaxInd; + dfBandDensity = 1; + bHasFoundDensity = TRUE; + } + } + + } // GRA_Mode + +/* -------------------------------------------------------------------- */ +/* We have a computed value from the source. Now apply it to */ +/* the destination pixel. */ +/* -------------------------------------------------------------------- */ + if ( bHasFoundDensity ) + { + // TODO should we compute dfBandDensity in fct of nCount/nCount2 , + // or use as a threshold to set the dest value? + // or fix gdalwarp crop_to_cutline to crop partially overlapping pixels + // dfBandDensity = (float) nCount / nCount2; + // if ( (float) nCount / nCount2 > 0.1 ) + GWKSetPixelValue( poWK, iBand, iDstOffset, + dfBandDensity, + dfValueReal, dfValueImag ); + } + } + + if (!bHasFoundDensity) + continue; + +/* -------------------------------------------------------------------- */ +/* Update destination density/validity masks. */ +/* -------------------------------------------------------------------- */ + GWKOverlayDensity( poWK, iDstOffset, dfDensity ); + + if( poWK->panDstValid != NULL ) + { + poWK->panDstValid[iDstOffset>>5] |= + 0x01 << (iDstOffset & 0x1f); + } + + } /* Next iDstX */ + +/* -------------------------------------------------------------------- */ +/* Report progress to the user, and optionally cancel out. */ +/* -------------------------------------------------------------------- */ + if (psJob->pfnProgress(psJob)) + break; + } + +/* -------------------------------------------------------------------- */ +/* Cleanup and return. */ +/* -------------------------------------------------------------------- */ + CPLFree( padfX ); + CPLFree( padfY ); + CPLFree( padfZ ); + CPLFree( pabSuccess ); +} + Index: alg/gdalwarper.cpp =================================================================== --- alg/gdalwarper.cpp (revision 25854) +++ alg/gdalwarper.cpp (working copy) @@ -968,6 +968,10 @@ pszAlgName = "CubicSpline"; else if( psWO->eResampleAlg == GRA_Lanczos ) pszAlgName = "Lanczos"; + else if( psWO->eResampleAlg == GRA_Average ) + pszAlgName = "Average"; + else if( psWO->eResampleAlg == GRA_Mode ) + pszAlgName = "Mode"; else pszAlgName = "Unknown"; @@ -1187,6 +1191,10 @@ psWO->eResampleAlg = GRA_CubicSpline; else if( EQUAL(pszValue,"Lanczos") ) psWO->eResampleAlg = GRA_Lanczos; + else if( EQUAL(pszValue,"Average") ) + psWO->eResampleAlg = GRA_Average; + else if( EQUAL(pszValue,"Mode") ) + psWO->eResampleAlg = GRA_Mode; else if( EQUAL(pszValue,"Default") ) /* leave as is */; else Index: alg/gdalwarpoperation.cpp =================================================================== --- alg/gdalwarpoperation.cpp (revision 25854) +++ alg/gdalwarpoperation.cpp (working copy) @@ -202,7 +202,9 @@ && psOptions->eResampleAlg != GRA_Bilinear && psOptions->eResampleAlg != GRA_Cubic && psOptions->eResampleAlg != GRA_CubicSpline - && psOptions->eResampleAlg != GRA_Lanczos ) + && psOptions->eResampleAlg != GRA_Lanczos + && psOptions->eResampleAlg != GRA_Average + && psOptions->eResampleAlg != GRA_Mode ) { CPLError( CE_Failure, CPLE_IllegalArg, "GDALWarpOptions.Validate()\n" Index: alg/gdalwarper.h =================================================================== --- alg/gdalwarper.h (revision 25854) +++ alg/gdalwarper.h (working copy) @@ -49,7 +49,9 @@ /*! Bilinear (2x2 kernel) */ GRA_Bilinear=1, /*! Cubic Convolution Approximation (4x4 kernel) */ GRA_Cubic=2, /*! Cubic B-Spline Approximation (4x4 kernel) */ GRA_CubicSpline=3, - /*! Lanczos windowed sinc interpolation (6x6 kernel) */ GRA_Lanczos=4 + /*! Lanczos windowed sinc interpolation (6x6 kernel) */ GRA_Lanczos=4, + /*! Average (variable kernel) */ GRA_Average=5, + /*! Mode (variable kernel) */ GRA_Mode=6 } GDALResampleAlg; typedef int