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