diff --git a/alg/gdal_tps.cpp b/alg/gdal_tps.cpp index dbf53037df7b..e09e3fa9ebc7 100644 --- a/alg/gdal_tps.cpp +++ b/alg/gdal_tps.cpp @@ -62,6 +62,7 @@ typedef struct VizGeorefSpline2D *poReverse; bool bForwardSolved; bool bReverseSolved; + double dfSrcApproxErrorReverse; bool bReversed; @@ -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) { @@ -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]; } @@ -429,6 +436,13 @@ CPLXMLNode *GDALSerializeTPSTransformer(void *pTransformArg) nullptr); } + if (psInfo->dfSrcApproxErrorReverse > 0) + { + CPLCreateXMLElementAndValue( + psTree, "SrcApproxErrorInPixel", + CPLString().Printf("%g", psInfo->dfSrcApproxErrorReverse)); + } + return psTree; } @@ -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. */ diff --git a/alg/gdalgenericinverse.cpp b/alg/gdalgenericinverse.cpp index 6a3b7232acbf..cccc90ba22a9 100644 --- a/alg/gdalgenericinverse.cpp +++ b/alg/gdalgenericinverse.cpp @@ -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; @@ -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; } diff --git a/alg/gdalgenericinverse.h b/alg/gdalgenericinverse.h index 49ee999fd718..9798b7718add 100644 --- a/alg/gdalgenericinverse.h +++ b/alg/gdalgenericinverse.h @@ -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 diff --git a/apps/gdalwarp_lib.cpp b/apps/gdalwarp_lib.cpp index 6fa13e45feac..2b121191de07 100644 --- a/apps/gdalwarp_lib.cpp +++ b/apps/gdalwarp_lib.cpp @@ -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, @@ -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")))) { diff --git a/autotest/gcore/transformer.py b/autotest/gcore/transformer.py index 4cf92610cbb5..3fa5d65e33c5 100644 --- a/autotest/gcore/transformer.py +++ b/autotest/gcore/transformer.py @@ -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)