Skip to content

Commit

Permalink
Inverse TPS transformer: speed improvement in gdalwarp use case (fixes …
Browse files Browse the repository at this point in the history 
…OSGeo#8672)

On OSGeo#8672 use case, timings are:
- 3.7.3: ~ 25 minutes
- 3.8.0: ~ 1400 minutes
- this PR: ~ 140 minutes
  • Loading branch information
@rouault
rouault committed Nov 14, 2023
1 parent 61b83a3 commit 0792682
Show file tree
Hide file tree
Showing 5 changed files with 101 additions and 36 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
14 changes: 14 additions & 0 deletions autotest/gcore/transformer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1049,10 +1049,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 0792682

Please sign in to comment.