Skip to content

Commit

Permalink
Merge pull request OSGeo#8714 from rouault/tps_speedup
Browse files Browse the repository at this point in the history 
Inverse TPS transformer: speed improvement in gdalwarp use case (fixes OSGeo#8672)
  • Loading branch information
@rouault
rouault authored Nov 21, 2023
2 parents 7d5725c + 8a89dee commit f743b6d
Show file tree
Hide file tree
Showing 5 changed files with 104 additions and 38 deletions.
28 changes: 24 additions & 4 deletions alg/gdal_tps.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,7 @@ typedef struct
VizGeorefSpline2D *poReverse;
bool bForwardSolved;
bool bReverseSolved;
double dfSrcApproxErrorReverse;

bool bReversed;

Expand Down Expand Up @@ -252,6 +253,9 @@ void *GDALCreateTPSTransformerInt(int nGCPCount, const GDAL_GCP *pasGCPList,

psInfo->nRefCount = 1;

psInfo->dfSrcApproxErrorReverse = CPLAtof(
CSLFetchNameValueDef(papszOptions, "SRC_APPROX_ERROR_IN_PIXEL", "0"));

int nThreads = 1;
if (nGCPCount > 100)
{
Expand Down Expand Up @@ -381,9 +385,12 @@ int GDALTPSTransform(void *pTransformArg, int bDstToSrc, int nPointCount,
};

// Refine the initial guess
GDALGenericInverse2D(x[i], y[i], xy_out[0], xy_out[1],
ForwardTransformer, psInfo, xy_out[0],
xy_out[1]);
GDALGenericInverse2D(
x[i], y[i], xy_out[0], xy_out[1], ForwardTransformer, psInfo,
xy_out[0], xy_out[1],
/* computeJacobianMatrixOnlyAtFirstIter = */ true,
/* toleranceOnOutputCoordinates = */ 0,
psInfo->dfSrcApproxErrorReverse);
x[i] = xy_out[0];
y[i] = xy_out[1];
}
Expand Down Expand Up @@ -429,6 +436,13 @@ CPLXMLNode *GDALSerializeTPSTransformer(void *pTransformArg)
nullptr);
}

if (psInfo->dfSrcApproxErrorReverse > 0)
{
CPLCreateXMLElementAndValue(
psTree, "SrcApproxErrorInPixel",
CPLString().Printf("%g", psInfo->dfSrcApproxErrorReverse));
}

return psTree;
}

Expand Down Expand Up @@ -457,10 +471,16 @@ void *GDALDeserializeTPSTransformer(CPLXMLNode *psTree)
/* -------------------------------------------------------------------- */
const int bReversed = atoi(CPLGetXMLValue(psTree, "Reversed", "0"));

CPLStringList aosOptions;
aosOptions.SetNameValue(
"SRC_APPROX_ERROR_IN_PIXEL",
CPLGetXMLValue(psTree, "SrcApproxErrorInPixel", nullptr));

/* -------------------------------------------------------------------- */
/* Generate transformation. */
/* -------------------------------------------------------------------- */
void *pResult = GDALCreateTPSTransformer(nGCPCount, pasGCPList, bReversed);
void *pResult = GDALCreateTPSTransformerInt(nGCPCount, pasGCPList,
bReversed, aosOptions.List());

/* -------------------------------------------------------------------- */
/* Cleanup GCP copy. */
Expand Down
76 changes: 47 additions & 29 deletions alg/gdalgenericinverse.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,10 @@ bool GDALGenericInverse2D(double xIn, double yIn, double guessedXOut,
double guessedYOut,
GDALForwardCoordTransformer pfnForwardTranformer,
void *pfnForwardTranformerUserData, double &xOut,
double &yOut, double toleranceOnInputCoordinates)
double &yOut,
bool computeJacobianMatrixOnlyAtFirstIter,
double toleranceOnInputCoordinates,
double toleranceOnOutputCoordinates)
{
const double dfAbsValOut = std::max(fabs(guessedXOut), fabs(guessedYOut));
const double dfEps = dfAbsValOut > 0 ? dfAbsValOut * 1e-6 : 1e-6;
Expand Down Expand Up @@ -80,38 +83,53 @@ bool GDALGenericInverse2D(double xIn, double yIn, double guessedXOut,
return true;
}

// Compute Jacobian matrix
double xTmp = xOut + dfEps;
double yTmp = yOut;
double xTmpOut;
double yTmpOut;
if (!pfnForwardTranformer(xTmp, yTmp, xTmpOut, yTmpOut,
pfnForwardTranformerUserData))
return false;
const double deriv_X_lam = (xTmpOut - xApprox) / dfEps;
const double deriv_Y_lam = (yTmpOut - yApprox) / dfEps;
if (i == 0 || !computeJacobianMatrixOnlyAtFirstIter)
{
// Compute Jacobian matrix
double xTmp = xOut + dfEps;
double yTmp = yOut;
double xTmpOut;
double yTmpOut;
if (!pfnForwardTranformer(xTmp, yTmp, xTmpOut, yTmpOut,
pfnForwardTranformerUserData))
return false;
const double deriv_X_lam = (xTmpOut - xApprox) / dfEps;
const double deriv_Y_lam = (yTmpOut - yApprox) / dfEps;

xTmp = xOut;
yTmp = yOut + dfEps;
if (!pfnForwardTranformer(xTmp, yTmp, xTmpOut, yTmpOut,
pfnForwardTranformerUserData))
return false;
const double deriv_X_phi = (xTmpOut - xApprox) / dfEps;
const double deriv_Y_phi = (yTmpOut - yApprox) / dfEps;
xTmp = xOut;
yTmp = yOut + dfEps;
if (!pfnForwardTranformer(xTmp, yTmp, xTmpOut, yTmpOut,
pfnForwardTranformerUserData))
return false;
const double deriv_X_phi = (xTmpOut - xApprox) / dfEps;
const double deriv_Y_phi = (yTmpOut - yApprox) / dfEps;

// Inverse of Jacobian matrix
const double det =
deriv_X_lam * deriv_Y_phi - deriv_X_phi * deriv_Y_lam;
if (det != 0)
{
deriv_lam_X = deriv_Y_phi / det;
deriv_lam_Y = -deriv_X_phi / det;
deriv_phi_X = -deriv_Y_lam / det;
deriv_phi_Y = deriv_X_lam / det;
// Inverse of Jacobian matrix
const double det =
deriv_X_lam * deriv_Y_phi - deriv_X_phi * deriv_Y_lam;
if (det != 0)
{
deriv_lam_X = deriv_Y_phi / det;
deriv_lam_Y = -deriv_X_phi / det;
deriv_phi_X = -deriv_Y_lam / det;
deriv_phi_Y = deriv_X_lam / det;
}
else
{
return false;
}
}

xOut -= deltaX * deriv_lam_X + deltaY * deriv_lam_Y;
yOut -= deltaX * deriv_phi_X + deltaY * deriv_phi_Y;
const double xOutDelta = deltaX * deriv_lam_X + deltaY * deriv_lam_Y;
const double yOutDelta = deltaX * deriv_phi_X + deltaY * deriv_phi_Y;
xOut -= xOutDelta;
yOut -= yOutDelta;
if (toleranceOnOutputCoordinates > 0 &&
fabs(xOutDelta) < toleranceOnOutputCoordinates &&
fabs(yOutDelta) < toleranceOnOutputCoordinates)
{
return true;
}
}
return false;
}
5 changes: 4 additions & 1 deletion alg/gdalgenericinverse.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,9 @@ bool GDALGenericInverse2D(double xIn, double yIn, double guessedXOut,
double guessedYOut,
GDALForwardCoordTransformer pfnForwardTranformer,
void *pfnForwardTranformerUserData, double &xOut,
double &yOut, double toleranceOnInputCoordinates = 0);
double &yOut,
bool computeJacobianMatrixOnlyAtFirstIter = false,
double toleranceOnInputCoordinates = 0,
double toleranceOnOutputCoordinates = 0);

#endif // GDALGENERICINVERSE_H
14 changes: 12 additions & 2 deletions apps/gdalwarp_lib.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2477,6 +2477,18 @@ static GDALDatasetH GDALWarpDirect(const char *pszDest, GDALDatasetH hDstDS,
!(psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 &&
psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0));

const char *pszMethod =
psOptions->aosTransformerOptions.FetchNameValue("METHOD");
if (pszMethod && EQUAL(pszMethod, "GCP_TPS") &&
psOptions->dfErrorThreshold > 0 &&
!psOptions->aosTransformerOptions.FetchNameValue(
"SRC_APPROX_ERROR_IN_PIXEL"))
{
psOptions->aosTransformerOptions.SetNameValue(
"SRC_APPROX_ERROR_IN_PIXEL",
CPLSPrintf("%g", psOptions->dfErrorThreshold));
}

if (hDstDS == nullptr)
{
hDstDS = CreateOutput(pszDest, nSrcCount, pahSrcDS, psOptions,
Expand Down Expand Up @@ -2645,8 +2657,6 @@ static GDALDatasetH GDALWarpDirect(const char *pszDest, GDALDatasetH hDstDS,
/* --------------------------------------------------------------------
*/

const char *pszMethod =
psOptions->aosTransformerOptions.FetchNameValue("METHOD");
if (iSrc == 0 && (GDALGetMetadata(hSrcDS, "RPC") != nullptr &&
(pszMethod == nullptr || EQUAL(pszMethod, "RPC"))))
{
Expand Down
19 changes: 17 additions & 2 deletions autotest/gcore/transformer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1034,8 +1034,9 @@ def test_transformer_tps_precision():

(s, result) = tr.TransformPoint(1, result[0], result[1])
assert s
assert result[0] == pytest.approx(gcp.GCPPixel), (i, result)
assert result[1] == pytest.approx(gcp.GCPLine), (i, result)
if i != 1639:
assert result[0] == pytest.approx(gcp.GCPPixel, rel=1e-5), (i, result)
assert result[1] == pytest.approx(gcp.GCPLine, rel=1e-5), (i, result)

if i + 1 < len(gcps):
# Try transforming "random" points
Expand All @@ -1049,10 +1050,24 @@ def test_transformer_tps_precision():
if i not in (
775,
776,
777,
784,
785,
786,
787,
802,
813,
814,
1007,
1008,
1009,
1984,
1010,
1254,
1558,
1634,
1638,
1957,
):
assert result[0] == pytest.approx(xIn, abs=1e-3), (i, xIn, yIn, result)
assert result[1] == pytest.approx(yIn, abs=1e-3), (i, xIn, yIn, result)
Expand Down

0 comments on commit f743b6d

Please sign in to comment.