Skip to content

Commit

Permalink
Warper: fix shifted/aliased result with Lanczos resampling when XSCAL…
Browse files Browse the repository at this point in the history 
…E < 1 or YSCALE < 1

Fixes OSGeo#11042
  • Loading branch information
@rouault
rouault committed Oct 29, 2024
1 parent d75dee9 commit 406e5bb
Show file tree
Hide file tree
Showing 2 changed files with 193 additions and 49 deletions.
242 changes: 193 additions & 49 deletions alg/gdalwarpkernel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3501,11 +3501,19 @@ struct _GWKResampleWrkStruct
double *padfWeightsX;
bool *pabCalcX;

double *padfWeightsY; // Only used by GWKResampleOptimizedLanczos.
int iLastSrcX; // Only used by GWKResampleOptimizedLanczos.
int iLastSrcY; // Only used by GWKResampleOptimizedLanczos.
double dfLastDeltaX; // Only used by GWKResampleOptimizedLanczos.
double dfLastDeltaY; // Only used by GWKResampleOptimizedLanczos.
double *padfWeightsY; // Only used by GWKResampleOptimizedLanczos.
int iLastSrcX; // Only used by GWKResampleOptimizedLanczos.
int iLastSrcY; // Only used by GWKResampleOptimizedLanczos.
double dfLastDeltaX; // Only used by GWKResampleOptimizedLanczos.
double dfLastDeltaY; // Only used by GWKResampleOptimizedLanczos.
double dfCosPiXScale; // Only used by GWKResampleOptimizedLanczos.
double dfSinPiXScale; // Only used by GWKResampleOptimizedLanczos.
double dfCosPiXScaleOver3; // Only used by GWKResampleOptimizedLanczos.
double dfSinPiXScaleOver3; // Only used by GWKResampleOptimizedLanczos.
double dfCosPiYScale; // Only used by GWKResampleOptimizedLanczos.
double dfSinPiYScale; // Only used by GWKResampleOptimizedLanczos.
double dfCosPiYScaleOver3; // Only used by GWKResampleOptimizedLanczos.
double dfSinPiYScaleOver3; // Only used by GWKResampleOptimizedLanczos.

// Space for saving a row of pixels.
double *padfRowDensity;
Expand Down Expand Up @@ -3533,7 +3541,7 @@ static GWKResampleWrkStruct *GWKResampleCreateWrkStruct(GDALWarpKernel *poWK)
const int nYDist = (poWK->nYRadius + 1) * 2;

GWKResampleWrkStruct *psWrkStruct = static_cast<GWKResampleWrkStruct *>(
CPLMalloc(sizeof(GWKResampleWrkStruct)));
CPLCalloc(1, sizeof(GWKResampleWrkStruct)));

// Alloc space for saved X weights.
psWrkStruct->padfWeightsX =
Expand Down Expand Up @@ -3569,38 +3577,40 @@ static GWKResampleWrkStruct *GWKResampleCreateWrkStruct(GDALWarpKernel *poWK)
{
psWrkStruct->pfnGWKResample = GWKResampleOptimizedLanczos;

const double dfXScale = poWK->dfXScale;
if (dfXScale < 1.0)
if (poWK->dfXScale < 1)
{
int iMin = poWK->nFiltInitX;
int iMax = poWK->nXRadius;
while (iMin * dfXScale < -3.0)
iMin++;
while (iMax * dfXScale > 3.0)
iMax--;

for (int i = iMin; i <= iMax; ++i)
{
psWrkStruct->padfWeightsX[i - poWK->nFiltInitX] =
GWKLanczosSinc(i * dfXScale);
}
psWrkStruct->dfCosPiXScaleOver3 = cos(M_PI / 3 * poWK->dfXScale);
psWrkStruct->dfSinPiXScaleOver3 =
sqrt(1 - psWrkStruct->dfCosPiXScaleOver3 *
psWrkStruct->dfCosPiXScaleOver3);
// "Naive":
// const double dfCosPiXScale = cos( M_PI * dfXScale );
// const double dfSinPiXScale = sin( M_PI * dfXScale );
// but given that cos(3x) = 4 cos^3(x) - 3 cos(x) and x between 0 and M_PI
psWrkStruct->dfCosPiXScale = (4 * psWrkStruct->dfCosPiXScaleOver3 *
psWrkStruct->dfCosPiXScaleOver3 -
3) *
psWrkStruct->dfCosPiXScaleOver3;
psWrkStruct->dfSinPiXScale = sqrt(
1 - psWrkStruct->dfCosPiXScale * psWrkStruct->dfCosPiXScale);
}

const double dfYScale = poWK->dfYScale;
if (dfYScale < 1.0)
if (poWK->dfYScale < 1)
{
int jMin = poWK->nFiltInitY;
int jMax = poWK->nYRadius;
while (jMin * dfYScale < -3.0)
jMin++;
while (jMax * dfYScale > 3.0)
jMax--;

for (int j = jMin; j <= jMax; ++j)
{
psWrkStruct->padfWeightsY[j - poWK->nFiltInitY] =
GWKLanczosSinc(j * dfYScale);
}
psWrkStruct->dfCosPiYScaleOver3 = cos(M_PI / 3 * poWK->dfYScale);
psWrkStruct->dfSinPiYScaleOver3 =
sqrt(1 - psWrkStruct->dfCosPiYScaleOver3 *
psWrkStruct->dfCosPiYScaleOver3);
// "Naive":
// const double dfCosPiYScale = cos( M_PI * dfYScale );
// const double dfSinPiYScale = sin( M_PI * dfYScale );
// but given that cos(3x) = 4 cos^3(x) - 3 cos(x) and x between 0 and M_PI
psWrkStruct->dfCosPiYScale = (4 * psWrkStruct->dfCosPiYScaleOver3 *
psWrkStruct->dfCosPiYScaleOver3 -
3) *
psWrkStruct->dfCosPiYScaleOver3;
psWrkStruct->dfSinPiYScale = sqrt(
1 - psWrkStruct->dfCosPiYScale * psWrkStruct->dfCosPiYScale);
}
}
else
Expand Down Expand Up @@ -3840,11 +3850,86 @@ static bool GWKResampleOptimizedLanczos(const GDALWarpKernel *poWK, int iBand,

if (dfXScale < 1.0)
{
while (iMin * dfXScale < -3.0)
while ((iMin - dfDeltaX) * dfXScale < -3.0)
iMin++;
while (iMax * dfXScale > 3.0)
while ((iMax - dfDeltaX) * dfXScale > 3.0)
iMax--;
// padfWeightsX computed in GWKResampleCreateWrkStruct.

// clang-format off
/*
Naive version:
for (int i = iMin; i <= iMax; ++i)
{
psWrkStruct->padfWeightsX[i - poWK->nFiltInitX] =
GWKLanczosSinc((i - dfDeltaX) * dfXScale);
}
but given that:
GWKLanczosSinc(x):
if (dfX == 0.0)
return 1.0;
const double dfPIX = M_PI * dfX;
const double dfPIXoverR = dfPIX / 3;
const double dfPIX2overR = dfPIX * dfPIXoverR;
return sin(dfPIX) * sin(dfPIXoverR) / dfPIX2overR;
and
sin (a + b) = sin a cos b + cos a sin b.
cos (a + b) = cos a cos b - sin a sin b.
we can skip any sin() computation within the loop
*/
// clang-format on

if (iSrcX != psWrkStruct->iLastSrcX ||
dfDeltaX != psWrkStruct->dfLastDeltaX)
{
double dfX = (iMin - dfDeltaX) * dfXScale;

double dfPIXover3 = M_PI / 3 * dfX;
double dfCosOver3 = cos(dfPIXover3);
double dfSinOver3 = sin(dfPIXover3);

// "Naive":
// double dfSin = sin( M_PI * dfX );
// double dfCos = cos( M_PI * dfX );
// but given that cos(3x) = 4 cos^3(x) - 3 cos(x) and sin(3x) = 3 sin(x) - 4 sin^3 (x).
double dfSin = (3 - 4 * dfSinOver3 * dfSinOver3) * dfSinOver3;
double dfCos = (4 * dfCosOver3 * dfCosOver3 - 3) * dfCosOver3;

const double dfCosPiXScaleOver3 = psWrkStruct->dfCosPiXScaleOver3;
const double dfSinPiXScaleOver3 = psWrkStruct->dfSinPiXScaleOver3;
const double dfCosPiXScale = psWrkStruct->dfCosPiXScale;
const double dfSinPiXScale = psWrkStruct->dfSinPiXScale;
constexpr double THREE_PI_PI = 3 * M_PI * M_PI;
padfWeightsX[iMin - poWK->nFiltInitX] =
dfX == 0 ? 1.0 : THREE_PI_PI * dfSin * dfSinOver3 / (dfX * dfX);
for (int i = iMin + 1; i <= iMax; ++i)
{
dfX += dfXScale;
const double dfNewSin =
dfSin * dfCosPiXScale + dfCos * dfSinPiXScale;
const double dfNewSinOver3 = dfSinOver3 * dfCosPiXScaleOver3 +
dfCosOver3 * dfSinPiXScaleOver3;
padfWeightsX[i - poWK->nFiltInitX] =
dfX == 0
? 1.0
: THREE_PI_PI * dfNewSin * dfNewSinOver3 / (dfX * dfX);
const double dfNewCos =
dfCos * dfCosPiXScale - dfSin * dfSinPiXScale;
const double dfNewCosOver3 = dfCosOver3 * dfCosPiXScaleOver3 -
dfSinOver3 * dfSinPiXScaleOver3;
dfSin = dfNewSin;
dfCos = dfNewCos;
dfSinOver3 = dfNewSinOver3;
dfCosOver3 = dfNewCosOver3;
}

psWrkStruct->iLastSrcX = iSrcX;
psWrkStruct->dfLastDeltaX = dfDeltaX;
}
}
else
{
Expand All @@ -3860,17 +3945,18 @@ static bool GWKResampleOptimizedLanczos(const GDALWarpKernel *poWK, int iBand,
// following trigonometric formulas.

// TODO(schwehr): Move this somewhere where it can be rendered at
// LaTeX. sin(M_PI * (dfBase + k)) = sin(M_PI * dfBase) * cos(M_PI *
// k) + cos(M_PI * dfBase) * sin(M_PI * k) sin(M_PI * (dfBase + k))
// = dfSinPIBase * cos(M_PI * k) + dfCosPIBase * sin(M_PI * k)
// LaTeX.
// clang-format off
// sin(M_PI * (dfBase + k)) = sin(M_PI * dfBase) * cos(M_PI * k) +
// cos(M_PI * dfBase) * sin(M_PI * k)
// sin(M_PI * (dfBase + k)) = dfSinPIBase * cos(M_PI * k) + dfCosPIBase * sin(M_PI * k)
// sin(M_PI * (dfBase + k)) = dfSinPIBase * cos(M_PI * k)
// sin(M_PI * (dfBase + k)) = dfSinPIBase * (((k % 2) == 0) ? 1 :
// -1)
// sin(M_PI * (dfBase + k)) = dfSinPIBase * (((k % 2) == 0) ? 1 : -1)

// sin(M_PI / dfR * (dfBase + k)) = sin(M_PI / dfR * dfBase) *
// cos(M_PI / dfR * k) + cos(M_PI / dfR * dfBase) * sin(M_PI / dfR *
// k) sin(M_PI / dfR * (dfBase + k)) = dfSinPIBaseOverR * cos(M_PI /
// dfR * k) + dfCosPIBaseOverR * sin(M_PI / dfR * k)
// sin(M_PI / dfR * (dfBase + k)) = sin(M_PI / dfR * dfBase) * cos(M_PI / dfR * k) +
// cos(M_PI / dfR * dfBase) * sin(M_PI / dfR * k)
// sin(M_PI / dfR * (dfBase + k)) = dfSinPIBaseOverR * cos(M_PI / dfR * k) + dfCosPIBaseOverR * sin(M_PI / dfR * k)
// clang-format on

const double dfSinPIDeltaXOver3 = sin((-M_PI / 3.0) * dfDeltaX);
const double dfSin2PIDeltaXOver3 =
Expand Down Expand Up @@ -3916,11 +4002,69 @@ static bool GWKResampleOptimizedLanczos(const GDALWarpKernel *poWK, int iBand,

if (dfYScale < 1.0)
{
while (jMin * dfYScale < -3.0)
while ((jMin - dfDeltaY) * dfYScale < -3.0)
jMin++;
while (jMax * dfYScale > 3.0)
while ((jMax - dfDeltaY) * dfYScale > 3.0)
jMax--;
// padfWeightsY computed in GWKResampleCreateWrkStruct.

// clang-format off
/*
Naive version:
for (int j = jMin; j <= jMax; ++j)
{
psWrkStruct->padfWeightsY[j - poWK->nFiltInitY] =
GWKLanczosSinc((j - dfDeltaY) * dfYScale);
}
*/
// clang-format on

if (iSrcY != psWrkStruct->iLastSrcY ||
dfDeltaY != psWrkStruct->dfLastDeltaY)
{
double dfY = (jMin - dfDeltaY) * dfYScale;

double dfPIYover3 = M_PI / 3 * dfY;
double dfCosOver3 = cos(dfPIYover3);
double dfSinOver3 = sin(dfPIYover3);

// "Naive":
// double dfSin = sin( M_PI * dfY );
// double dfCos = cos( M_PI * dfY );
// but given that cos(3x) = 4 cos^3(x) - 3 cos(x) and sin(3x) = 3 sin(x) - 4 sin^3 (x).
double dfSin = (3 - 4 * dfSinOver3 * dfSinOver3) * dfSinOver3;
double dfCos = (4 * dfCosOver3 * dfCosOver3 - 3) * dfCosOver3;

const double dfCosPiYScaleOver3 = psWrkStruct->dfCosPiYScaleOver3;
const double dfSinPiYScaleOver3 = psWrkStruct->dfSinPiYScaleOver3;
const double dfCosPiYScale = psWrkStruct->dfCosPiYScale;
const double dfSinPiYScale = psWrkStruct->dfSinPiYScale;
constexpr double THREE_PI_PI = 3 * M_PI * M_PI;
padfWeightsY[jMin - poWK->nFiltInitY] =
dfY == 0 ? 1.0 : THREE_PI_PI * dfSin * dfSinOver3 / (dfY * dfY);
for (int j = jMin + 1; j <= iMax; ++j)
{
dfY += dfYScale;
const double dfNewSin =
dfSin * dfCosPiYScale + dfCos * dfSinPiYScale;
const double dfNewSinOver3 = dfSinOver3 * dfCosPiYScaleOver3 +
dfCosOver3 * dfSinPiYScaleOver3;
padfWeightsY[j - poWK->nFiltInitY] =
dfY == 0
? 1.0
: THREE_PI_PI * dfNewSin * dfNewSinOver3 / (dfY * dfY);
const double dfNewCos =
dfCos * dfCosPiYScale - dfSin * dfSinPiYScale;
const double dfNewCosOver3 = dfCosOver3 * dfCosPiYScaleOver3 -
dfSinOver3 * dfSinPiYScaleOver3;
dfSin = dfNewSin;
dfCos = dfNewCos;
dfSinOver3 = dfNewSinOver3;
dfCosOver3 = dfNewCosOver3;
}

psWrkStruct->iLastSrcY = iSrcY;
psWrkStruct->dfLastDeltaY = dfDeltaY;
}
}
else
{
Expand Down
Binary file modified autotest/alg/data/utmsmall_lanczos_2.tif
View file
Binary file not shown.

0 comments on commit 406e5bb

Please sign in to comment.