Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: High Performance Image Reprojector
4 : * Purpose: Thin Plate Spline transformer (GDAL wrapper portion)
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2004, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2011-2013, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_port.h"
15 : #include "thinplatespline.h"
16 :
17 : #include <stdlib.h>
18 : #include <string.h>
19 :
20 : #include <algorithm>
21 : #include <map>
22 : #include <utility>
23 :
24 : #include "cpl_atomic_ops.h"
25 : #include "cpl_conv.h"
26 : #include "cpl_error.h"
27 : #include "cpl_minixml.h"
28 : #include "cpl_multiproc.h"
29 : #include "cpl_string.h"
30 : #include "cpl_vsi.h" // CPLGetUsablePhysicalRAM()
31 : #include "gdal.h"
32 : #include "gdal_alg.h"
33 : #include "gdal_alg_priv.h"
34 : #include "gdal_priv.h"
35 : #include "gdalgenericinverse.h"
36 :
37 : CPL_C_START
38 : CPLXMLNode *GDALSerializeTPSTransformer(void *pTransformArg);
39 : void *GDALDeserializeTPSTransformer(CPLXMLNode *psTree);
40 : CPL_C_END
41 :
42 : struct TPSTransformInfo
43 : {
44 : GDALTransformerInfo sTI{};
45 :
46 : VizGeorefSpline2D *poForward{};
47 : VizGeorefSpline2D *poReverse{};
48 : bool bForwardSolved{};
49 : bool bReverseSolved{};
50 : double dfSrcApproxErrorReverse{};
51 :
52 : bool bReversed{};
53 : bool bForceBuiltinMethod = false;
54 :
55 : std::vector<gdal::GCP> asGCPs{};
56 :
57 : volatile int nRefCount{};
58 : };
59 :
60 : /************************************************************************/
61 : /* GDALCreateSimilarTPSTransformer() */
62 : /************************************************************************/
63 :
64 3 : static void *GDALCreateSimilarTPSTransformer(void *hTransformArg,
65 : double dfRatioX, double dfRatioY)
66 : {
67 3 : VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarTPSTransformer",
68 : nullptr);
69 :
70 3 : TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(hTransformArg);
71 :
72 3 : if (dfRatioX == 1.0 && dfRatioY == 1.0)
73 : {
74 : // We can just use a ref count, since using the source transformation
75 : // is thread-safe.
76 0 : CPLAtomicInc(&(psInfo->nRefCount));
77 : }
78 : else
79 : {
80 3 : auto newGCPs = psInfo->asGCPs;
81 21 : for (auto &gcp : newGCPs)
82 : {
83 18 : gcp.Pixel() /= dfRatioX;
84 18 : gcp.Line() /= dfRatioY;
85 : }
86 3 : psInfo = static_cast<TPSTransformInfo *>(GDALCreateTPSTransformer(
87 3 : static_cast<int>(newGCPs.size()), gdal::GCP::c_ptr(newGCPs),
88 3 : psInfo->bReversed));
89 : }
90 :
91 3 : return psInfo;
92 : }
93 :
94 : /************************************************************************/
95 : /* GDALCreateTPSTransformer() */
96 : /************************************************************************/
97 :
98 : /**
99 : * Create Thin Plate Spline transformer from GCPs.
100 : *
101 : * The thin plate spline transformer produces exact transformation
102 : * at all control points and smoothly varying transformations between
103 : * control points with greatest influence from local control points.
104 : * It is suitable for for many applications not well modeled by polynomial
105 : * transformations.
106 : *
107 : * Creating the TPS transformer involves solving systems of linear equations
108 : * related to the number of control points involved. This solution is
109 : * computed within this function call. It can be quite an expensive operation
110 : * for large numbers of GCPs. For instance, for reference, it takes on the
111 : * order of 10s for 400 GCPs on a 2GHz Athlon processor.
112 : *
113 : * TPS Transformers are serializable.
114 : *
115 : * The GDAL Thin Plate Spline transformer is based on code provided by
116 : * Gilad Ronnen on behalf of VIZRT Inc (http://www.visrt.com). Incorporation
117 : * of the algorithm into GDAL was supported by the Centro di Ecologia Alpina
118 : * (http://www.cealp.it).
119 : *
120 : * @param nGCPCount the number of GCPs in pasGCPList.
121 : * @param pasGCPList an array of GCPs to be used as input.
122 : * @param bReversed set it to TRUE to compute the reversed transformation.
123 : *
124 : * @return the transform argument or NULL if creation fails.
125 : */
126 :
127 5 : void *GDALCreateTPSTransformer(int nGCPCount, const GDAL_GCP *pasGCPList,
128 : int bReversed)
129 : {
130 5 : return GDALCreateTPSTransformerInt(nGCPCount, pasGCPList, bReversed,
131 5 : nullptr);
132 : }
133 :
134 0 : static void GDALTPSComputeForwardInThread(void *pData)
135 : {
136 0 : TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(pData);
137 0 : psInfo->bForwardSolved =
138 0 : psInfo->poForward->solve(psInfo->bForceBuiltinMethod) != 0;
139 0 : }
140 :
141 19 : void *GDALCreateTPSTransformerInt(int nGCPCount, const GDAL_GCP *pasGCPList,
142 : int bReversed, CSLConstList papszOptions)
143 :
144 : {
145 : /* -------------------------------------------------------------------- */
146 : /* Allocate transform info. */
147 : /* -------------------------------------------------------------------- */
148 19 : TPSTransformInfo *psInfo = new TPSTransformInfo();
149 :
150 19 : psInfo->asGCPs = gdal::GCP::fromC(pasGCPList, nGCPCount);
151 :
152 19 : psInfo->bReversed = CPL_TO_BOOL(bReversed);
153 19 : psInfo->poForward = new VizGeorefSpline2D(2);
154 19 : psInfo->poReverse = new VizGeorefSpline2D(2);
155 :
156 19 : memcpy(psInfo->sTI.abySignature, GDAL_GTI2_SIGNATURE,
157 : strlen(GDAL_GTI2_SIGNATURE));
158 19 : psInfo->sTI.pszClassName = "GDALTPSTransformer";
159 19 : psInfo->sTI.pfnTransform = GDALTPSTransform;
160 19 : psInfo->sTI.pfnCleanup = GDALDestroyTPSTransformer;
161 19 : psInfo->sTI.pfnSerialize = GDALSerializeTPSTransformer;
162 19 : psInfo->sTI.pfnCreateSimilar = GDALCreateSimilarTPSTransformer;
163 :
164 : /* -------------------------------------------------------------------- */
165 : /* Attach (non-redundant) points to the transformation. */
166 : /* -------------------------------------------------------------------- */
167 38 : std::map<std::pair<double, double>, int> oMapPixelLineToIdx;
168 38 : std::map<std::pair<double, double>, int> oMapXYToIdx;
169 2426 : for (int iGCP = 0; iGCP < nGCPCount; iGCP++)
170 : {
171 2407 : const double afPL[2] = {pasGCPList[iGCP].dfGCPPixel,
172 2407 : pasGCPList[iGCP].dfGCPLine};
173 2407 : const double afXY[2] = {pasGCPList[iGCP].dfGCPX,
174 2407 : pasGCPList[iGCP].dfGCPY};
175 :
176 : std::map<std::pair<double, double>, int>::iterator oIter(
177 : oMapPixelLineToIdx.find(
178 2407 : std::pair<double, double>(afPL[0], afPL[1])));
179 2407 : if (oIter != oMapPixelLineToIdx.end())
180 : {
181 3 : if (afXY[0] == pasGCPList[oIter->second].dfGCPX &&
182 1 : afXY[1] == pasGCPList[oIter->second].dfGCPY)
183 : {
184 1 : continue;
185 : }
186 : else
187 : {
188 1 : CPLError(CE_Warning, CPLE_AppDefined,
189 : "GCP %d and %d have same (pixel,line)=(%f,%f), "
190 : "but different (X,Y): (%f,%f) vs (%f,%f)",
191 1 : iGCP + 1, oIter->second, afPL[0], afPL[1], afXY[0],
192 1 : afXY[1], pasGCPList[oIter->second].dfGCPX,
193 1 : pasGCPList[oIter->second].dfGCPY);
194 : }
195 : }
196 : else
197 : {
198 2405 : oMapPixelLineToIdx[std::pair<double, double>(afPL[0], afPL[1])] =
199 : iGCP;
200 : }
201 :
202 2406 : oIter = oMapXYToIdx.find(std::pair<double, double>(afXY[0], afXY[1]));
203 2406 : if (oIter != oMapXYToIdx.end())
204 : {
205 1 : CPLError(CE_Warning, CPLE_AppDefined,
206 : "GCP %d and %d have same (x,y)=(%f,%f), "
207 : "but different (pixel,line): (%f,%f) vs (%f,%f)",
208 1 : iGCP + 1, oIter->second, afXY[0], afXY[1], afPL[0],
209 1 : afPL[1], pasGCPList[oIter->second].dfGCPPixel,
210 1 : pasGCPList[oIter->second].dfGCPLine);
211 : }
212 : else
213 : {
214 2405 : oMapXYToIdx[std::pair<double, double>(afXY[0], afXY[1])] = iGCP;
215 : }
216 :
217 2406 : bool bOK = true;
218 2406 : if (bReversed)
219 : {
220 0 : bOK &= psInfo->poReverse->add_point(afPL[0], afPL[1], afXY);
221 0 : bOK &= psInfo->poForward->add_point(afXY[0], afXY[1], afPL);
222 : }
223 : else
224 : {
225 2406 : bOK &= psInfo->poForward->add_point(afPL[0], afPL[1], afXY);
226 2406 : bOK &= psInfo->poReverse->add_point(afXY[0], afXY[1], afPL);
227 : }
228 2406 : if (!bOK)
229 : {
230 0 : GDALDestroyTPSTransformer(psInfo);
231 0 : return nullptr;
232 : }
233 : }
234 :
235 19 : psInfo->nRefCount = 1;
236 :
237 19 : psInfo->dfSrcApproxErrorReverse = CPLAtof(
238 : CSLFetchNameValueDef(papszOptions, "SRC_APPROX_ERROR_IN_PIXEL", "0"));
239 :
240 19 : int nThreads = 1;
241 19 : bool bForceBuiltinMethod = false;
242 : // Arbitrary threshold beyond which multithreading might be interesting,
243 : // and checking memory usage too...
244 19 : if (nGCPCount > 100)
245 : {
246 : const char *pszWarpThreads =
247 2 : CSLFetchNameValue(papszOptions, "NUM_THREADS");
248 2 : if (pszWarpThreads == nullptr)
249 2 : pszWarpThreads = CPLGetConfigOption("GDAL_NUM_THREADS", "1");
250 2 : if (EQUAL(pszWarpThreads, "ALL_CPUS"))
251 0 : nThreads = CPLGetNumCPUs();
252 : else
253 2 : nThreads = atoi(pszWarpThreads);
254 2 : nThreads = std::clamp(nThreads, 1, 2);
255 :
256 : // Do sanity checks w.r.t. available RAM
257 :
258 2 : const auto nRAM = CPLGetUsablePhysicalRAM();
259 2 : if (nRAM > 0)
260 : {
261 2 : const int nMatrixSize = nGCPCount + 3;
262 : #ifdef HAVE_ARMADILLO
263 : // Armadillo requires up to 3 matrices of size nMatrixSize x nMatrixSize
264 : // for each transformation direction
265 2 : constexpr int NUM_TEMP_MATRICES = 3;
266 2 : if (nMatrixSize >
267 2 : nRAM / (nThreads * NUM_TEMP_MATRICES * nMatrixSize *
268 : static_cast<int>(sizeof(double))))
269 : {
270 0 : if (nMatrixSize > nRAM / (NUM_TEMP_MATRICES * nMatrixSize *
271 : static_cast<int>(sizeof(double))))
272 : {
273 0 : CPLDebug("GDAL", "Not enough memory to use Armadillo "
274 : "solver for thinplatespline. Falling back "
275 : "to LU decomposition method");
276 0 : bForceBuiltinMethod = true;
277 : }
278 : else
279 : {
280 0 : nThreads = 1;
281 : }
282 : }
283 2 : if (bForceBuiltinMethod)
284 : #endif
285 : {
286 0 : if (nMatrixSize > nRAM / (nThreads * nMatrixSize *
287 : static_cast<int>(sizeof(double))))
288 : {
289 0 : nThreads = 1;
290 0 : if (nMatrixSize >
291 0 : nRAM / (nMatrixSize * static_cast<int>(sizeof(double))))
292 : {
293 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
294 : "thinplatespline: not enough memory. At least "
295 : "%u MB are required",
296 : static_cast<unsigned>(
297 0 : static_cast<uint64_t>(nMatrixSize) *
298 0 : nMatrixSize * sizeof(double) /
299 : (1024 * 1024)));
300 0 : GDALDestroyTPSTransformer(psInfo);
301 0 : return nullptr;
302 : }
303 : }
304 : }
305 : }
306 : }
307 :
308 19 : psInfo->bForceBuiltinMethod = bForceBuiltinMethod;
309 :
310 19 : if (nThreads == 2)
311 : {
312 : // Compute direct and reverse transforms in parallel.
313 : CPLJoinableThread *hThread =
314 0 : CPLCreateJoinableThread(GDALTPSComputeForwardInThread, psInfo);
315 0 : psInfo->bReverseSolved =
316 0 : psInfo->poReverse->solve(bForceBuiltinMethod) != 0;
317 0 : if (hThread != nullptr)
318 0 : CPLJoinThread(hThread);
319 : else
320 0 : psInfo->bForwardSolved =
321 0 : psInfo->poForward->solve(bForceBuiltinMethod) != 0;
322 : }
323 : else
324 : {
325 19 : psInfo->bForwardSolved =
326 19 : psInfo->poForward->solve(bForceBuiltinMethod) != 0;
327 19 : if (psInfo->bForwardSolved)
328 18 : psInfo->bReverseSolved =
329 18 : psInfo->poReverse->solve(bForceBuiltinMethod) != 0;
330 : }
331 :
332 19 : if (!psInfo->bForwardSolved || !psInfo->bReverseSolved)
333 : {
334 2 : GDALDestroyTPSTransformer(psInfo);
335 2 : return nullptr;
336 : }
337 :
338 17 : return psInfo;
339 : }
340 :
341 : /************************************************************************/
342 : /* GDALDestroyTPSTransformer() */
343 : /************************************************************************/
344 :
345 : /**
346 : * Destroy TPS transformer.
347 : *
348 : * This function is used to destroy information about a GCP based
349 : * polynomial transformation created with GDALCreateTPSTransformer().
350 : *
351 : * @param pTransformArg the transform arg previously returned by
352 : * GDALCreateTPSTransformer().
353 : */
354 :
355 19 : void GDALDestroyTPSTransformer(void *pTransformArg)
356 :
357 : {
358 19 : if (pTransformArg == nullptr)
359 0 : return;
360 :
361 19 : TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(pTransformArg);
362 :
363 19 : if (CPLAtomicDec(&(psInfo->nRefCount)) == 0)
364 : {
365 19 : delete psInfo->poForward;
366 19 : delete psInfo->poReverse;
367 :
368 19 : delete psInfo;
369 : }
370 : }
371 :
372 : /************************************************************************/
373 : /* GDALTPSTransform() */
374 : /************************************************************************/
375 :
376 : /**
377 : * Transforms point based on GCP derived polynomial model.
378 : *
379 : * This function matches the GDALTransformerFunc signature, and can be
380 : * used to transform one or more points from pixel/line coordinates to
381 : * georeferenced coordinates (SrcToDst) or vice versa (DstToSrc).
382 : *
383 : * @param pTransformArg return value from GDALCreateTPSTransformer().
384 : * @param bDstToSrc TRUE if transformation is from the destination
385 : * (georeferenced) coordinates to pixel/line or FALSE when transforming
386 : * from pixel/line to georeferenced coordinates.
387 : * @param nPointCount the number of values in the x, y and z arrays.
388 : * @param x array containing the X values to be transformed.
389 : * @param y array containing the Y values to be transformed.
390 : * @param z array containing the Z values to be transformed.
391 : * @param panSuccess array in which a flag indicating success (TRUE) or
392 : * failure (FALSE) of the transformation are placed.
393 : *
394 : * @return TRUE if all points have been successfully transformed.
395 : */
396 :
397 34806 : int GDALTPSTransform(void *pTransformArg, int bDstToSrc, int nPointCount,
398 : double *x, double *y, CPL_UNUSED double *z,
399 : int *panSuccess)
400 : {
401 34806 : VALIDATE_POINTER1(pTransformArg, "GDALTPSTransform", 0);
402 :
403 34806 : TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(pTransformArg);
404 :
405 118204 : for (int i = 0; i < nPointCount; i++)
406 : {
407 83398 : double xy_out[2] = {0.0, 0.0};
408 :
409 83398 : if (bDstToSrc)
410 : {
411 : // Compute initial guess
412 76166 : psInfo->poReverse->get_point(x[i], y[i], xy_out);
413 :
414 484456 : const auto ForwardTransformer = [](double xIn, double yIn,
415 : double &xOut, double &yOut,
416 : void *pUserData)
417 : {
418 484456 : double xyOut[2] = {0.0, 0.0};
419 484456 : TPSTransformInfo *l_psInfo =
420 : static_cast<TPSTransformInfo *>(pUserData);
421 484456 : l_psInfo->poForward->get_point(xIn, yIn, xyOut);
422 484456 : xOut = xyOut[0];
423 484456 : yOut = xyOut[1];
424 484456 : return true;
425 : };
426 :
427 : // Refine the initial guess
428 76166 : GDALGenericInverse2D(
429 76166 : x[i], y[i], xy_out[0], xy_out[1], ForwardTransformer, psInfo,
430 : xy_out[0], xy_out[1],
431 : /* computeJacobianMatrixOnlyAtFirstIter = */ true,
432 : /* toleranceOnOutputCoordinates = */ 0,
433 : psInfo->dfSrcApproxErrorReverse);
434 76166 : x[i] = xy_out[0];
435 76166 : y[i] = xy_out[1];
436 : }
437 : else
438 : {
439 7232 : psInfo->poForward->get_point(x[i], y[i], xy_out);
440 7232 : x[i] = xy_out[0];
441 7232 : y[i] = xy_out[1];
442 : }
443 83398 : panSuccess[i] = TRUE;
444 : }
445 :
446 34806 : return TRUE;
447 : }
448 :
449 : /************************************************************************/
450 : /* GDALSerializeTPSTransformer() */
451 : /************************************************************************/
452 :
453 2 : CPLXMLNode *GDALSerializeTPSTransformer(void *pTransformArg)
454 :
455 : {
456 2 : VALIDATE_POINTER1(pTransformArg, "GDALSerializeTPSTransformer", nullptr);
457 :
458 2 : TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(pTransformArg);
459 :
460 : CPLXMLNode *psTree =
461 2 : CPLCreateXMLNode(nullptr, CXT_Element, "TPSTransformer");
462 :
463 : /* -------------------------------------------------------------------- */
464 : /* Serialize bReversed. */
465 : /* -------------------------------------------------------------------- */
466 2 : CPLCreateXMLElementAndValue(
467 : psTree, "Reversed",
468 4 : CPLString().Printf("%d", static_cast<int>(psInfo->bReversed)));
469 :
470 : /* -------------------------------------------------------------------- */
471 : /* Attach GCP List. */
472 : /* -------------------------------------------------------------------- */
473 2 : if (!psInfo->asGCPs.empty())
474 : {
475 2 : GDALSerializeGCPListToXML(psTree, psInfo->asGCPs, nullptr);
476 : }
477 :
478 2 : if (psInfo->dfSrcApproxErrorReverse > 0)
479 : {
480 2 : CPLCreateXMLElementAndValue(
481 : psTree, "SrcApproxErrorInPixel",
482 4 : CPLString().Printf("%g", psInfo->dfSrcApproxErrorReverse));
483 : }
484 :
485 2 : return psTree;
486 : }
487 :
488 : /************************************************************************/
489 : /* GDALDeserializeTPSTransformer() */
490 : /************************************************************************/
491 :
492 3 : void *GDALDeserializeTPSTransformer(CPLXMLNode *psTree)
493 :
494 : {
495 : /* -------------------------------------------------------------------- */
496 : /* Check for GCPs. */
497 : /* -------------------------------------------------------------------- */
498 3 : CPLXMLNode *psGCPList = CPLGetXMLNode(psTree, "GCPList");
499 :
500 6 : std::vector<gdal::GCP> asGCPs;
501 3 : if (psGCPList != nullptr)
502 : {
503 3 : GDALDeserializeGCPListFromXML(psGCPList, asGCPs, nullptr);
504 : }
505 :
506 : /* -------------------------------------------------------------------- */
507 : /* Get other flags. */
508 : /* -------------------------------------------------------------------- */
509 3 : const int bReversed = atoi(CPLGetXMLValue(psTree, "Reversed", "0"));
510 :
511 3 : CPLStringList aosOptions;
512 : aosOptions.SetNameValue(
513 : "SRC_APPROX_ERROR_IN_PIXEL",
514 3 : CPLGetXMLValue(psTree, "SrcApproxErrorInPixel", nullptr));
515 :
516 : /* -------------------------------------------------------------------- */
517 : /* Generate transformation. */
518 : /* -------------------------------------------------------------------- */
519 3 : void *pResult = GDALCreateTPSTransformerInt(static_cast<int>(asGCPs.size()),
520 : gdal::GCP::c_ptr(asGCPs),
521 3 : bReversed, aosOptions.List());
522 :
523 6 : return pResult;
524 : }
|
