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 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "cpl_port.h"
31 : #include "thinplatespline.h"
32 :
33 : #include <stdlib.h>
34 : #include <string.h>
35 : #include <map>
36 : #include <utility>
37 :
38 : #include "cpl_atomic_ops.h"
39 : #include "cpl_conv.h"
40 : #include "cpl_error.h"
41 : #include "cpl_minixml.h"
42 : #include "cpl_multiproc.h"
43 : #include "cpl_string.h"
44 : #include "gdal.h"
45 : #include "gdal_alg.h"
46 : #include "gdal_alg_priv.h"
47 : #include "gdal_priv.h"
48 : #include "gdalgenericinverse.h"
49 :
50 : CPL_C_START
51 : CPLXMLNode *GDALSerializeTPSTransformer(void *pTransformArg);
52 : void *GDALDeserializeTPSTransformer(CPLXMLNode *psTree);
53 : CPL_C_END
54 :
55 : struct TPSTransformInfo
56 : {
57 : GDALTransformerInfo sTI{};
58 :
59 : VizGeorefSpline2D *poForward{};
60 : VizGeorefSpline2D *poReverse{};
61 : bool bForwardSolved{};
62 : bool bReverseSolved{};
63 : double dfSrcApproxErrorReverse{};
64 :
65 : bool bReversed{};
66 :
67 : std::vector<gdal::GCP> asGCPs{};
68 :
69 : volatile int nRefCount{};
70 : };
71 :
72 : /************************************************************************/
73 : /* GDALCreateSimilarTPSTransformer() */
74 : /************************************************************************/
75 :
76 3 : static void *GDALCreateSimilarTPSTransformer(void *hTransformArg,
77 : double dfRatioX, double dfRatioY)
78 : {
79 3 : VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarTPSTransformer",
80 : nullptr);
81 :
82 3 : TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(hTransformArg);
83 :
84 3 : if (dfRatioX == 1.0 && dfRatioY == 1.0)
85 : {
86 : // We can just use a ref count, since using the source transformation
87 : // is thread-safe.
88 0 : CPLAtomicInc(&(psInfo->nRefCount));
89 : }
90 : else
91 : {
92 3 : auto newGCPs = psInfo->asGCPs;
93 21 : for (auto &gcp : newGCPs)
94 : {
95 18 : gcp.Pixel() /= dfRatioX;
96 18 : gcp.Line() /= dfRatioY;
97 : }
98 3 : psInfo = static_cast<TPSTransformInfo *>(GDALCreateTPSTransformer(
99 3 : static_cast<int>(newGCPs.size()), gdal::GCP::c_ptr(newGCPs),
100 3 : psInfo->bReversed));
101 : }
102 :
103 3 : return psInfo;
104 : }
105 :
106 : /************************************************************************/
107 : /* GDALCreateTPSTransformer() */
108 : /************************************************************************/
109 :
110 : /**
111 : * Create Thin Plate Spline transformer from GCPs.
112 : *
113 : * The thin plate spline transformer produces exact transformation
114 : * at all control points and smoothly varying transformations between
115 : * control points with greatest influence from local control points.
116 : * It is suitable for for many applications not well modeled by polynomial
117 : * transformations.
118 : *
119 : * Creating the TPS transformer involves solving systems of linear equations
120 : * related to the number of control points involved. This solution is
121 : * computed within this function call. It can be quite an expensive operation
122 : * for large numbers of GCPs. For instance, for reference, it takes on the
123 : * order of 10s for 400 GCPs on a 2GHz Athlon processor.
124 : *
125 : * TPS Transformers are serializable.
126 : *
127 : * The GDAL Thin Plate Spline transformer is based on code provided by
128 : * Gilad Ronnen on behalf of VIZRT Inc (http://www.visrt.com). Incorporation
129 : * of the algorithm into GDAL was supported by the Centro di Ecologia Alpina
130 : * (http://www.cealp.it).
131 : *
132 : * @param nGCPCount the number of GCPs in pasGCPList.
133 : * @param pasGCPList an array of GCPs to be used as input.
134 : * @param bReversed set it to TRUE to compute the reversed transformation.
135 : *
136 : * @return the transform argument or NULL if creation fails.
137 : */
138 :
139 5 : void *GDALCreateTPSTransformer(int nGCPCount, const GDAL_GCP *pasGCPList,
140 : int bReversed)
141 : {
142 5 : return GDALCreateTPSTransformerInt(nGCPCount, pasGCPList, bReversed,
143 5 : nullptr);
144 : }
145 :
146 0 : static void GDALTPSComputeForwardInThread(void *pData)
147 : {
148 0 : TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(pData);
149 0 : psInfo->bForwardSolved = psInfo->poForward->solve() != 0;
150 0 : }
151 :
152 19 : void *GDALCreateTPSTransformerInt(int nGCPCount, const GDAL_GCP *pasGCPList,
153 : int bReversed, char **papszOptions)
154 :
155 : {
156 : /* -------------------------------------------------------------------- */
157 : /* Allocate transform info. */
158 : /* -------------------------------------------------------------------- */
159 19 : TPSTransformInfo *psInfo = new TPSTransformInfo();
160 :
161 19 : psInfo->asGCPs = gdal::GCP::fromC(pasGCPList, nGCPCount);
162 :
163 19 : psInfo->bReversed = CPL_TO_BOOL(bReversed);
164 19 : psInfo->poForward = new VizGeorefSpline2D(2);
165 19 : psInfo->poReverse = new VizGeorefSpline2D(2);
166 :
167 19 : memcpy(psInfo->sTI.abySignature, GDAL_GTI2_SIGNATURE,
168 : strlen(GDAL_GTI2_SIGNATURE));
169 19 : psInfo->sTI.pszClassName = "GDALTPSTransformer";
170 19 : psInfo->sTI.pfnTransform = GDALTPSTransform;
171 19 : psInfo->sTI.pfnCleanup = GDALDestroyTPSTransformer;
172 19 : psInfo->sTI.pfnSerialize = GDALSerializeTPSTransformer;
173 19 : psInfo->sTI.pfnCreateSimilar = GDALCreateSimilarTPSTransformer;
174 :
175 : /* -------------------------------------------------------------------- */
176 : /* Attach (non-redundant) points to the transformation. */
177 : /* -------------------------------------------------------------------- */
178 38 : std::map<std::pair<double, double>, int> oMapPixelLineToIdx;
179 38 : std::map<std::pair<double, double>, int> oMapXYToIdx;
180 2426 : for (int iGCP = 0; iGCP < nGCPCount; iGCP++)
181 : {
182 2407 : const double afPL[2] = {pasGCPList[iGCP].dfGCPPixel,
183 2407 : pasGCPList[iGCP].dfGCPLine};
184 2407 : const double afXY[2] = {pasGCPList[iGCP].dfGCPX,
185 2407 : pasGCPList[iGCP].dfGCPY};
186 :
187 : std::map<std::pair<double, double>, int>::iterator oIter(
188 : oMapPixelLineToIdx.find(
189 2407 : std::pair<double, double>(afPL[0], afPL[1])));
190 2407 : if (oIter != oMapPixelLineToIdx.end())
191 : {
192 3 : if (afXY[0] == pasGCPList[oIter->second].dfGCPX &&
193 1 : afXY[1] == pasGCPList[oIter->second].dfGCPY)
194 : {
195 1 : continue;
196 : }
197 : else
198 : {
199 1 : CPLError(CE_Warning, CPLE_AppDefined,
200 : "GCP %d and %d have same (pixel,line)=(%f,%f), "
201 : "but different (X,Y): (%f,%f) vs (%f,%f)",
202 1 : iGCP + 1, oIter->second, afPL[0], afPL[1], afXY[0],
203 1 : afXY[1], pasGCPList[oIter->second].dfGCPX,
204 1 : pasGCPList[oIter->second].dfGCPY);
205 : }
206 : }
207 : else
208 : {
209 2405 : oMapPixelLineToIdx[std::pair<double, double>(afPL[0], afPL[1])] =
210 : iGCP;
211 : }
212 :
213 2406 : oIter = oMapXYToIdx.find(std::pair<double, double>(afXY[0], afXY[1]));
214 2406 : if (oIter != oMapXYToIdx.end())
215 : {
216 1 : CPLError(CE_Warning, CPLE_AppDefined,
217 : "GCP %d and %d have same (x,y)=(%f,%f), "
218 : "but different (pixel,line): (%f,%f) vs (%f,%f)",
219 1 : iGCP + 1, oIter->second, afXY[0], afXY[1], afPL[0],
220 1 : afPL[1], pasGCPList[oIter->second].dfGCPPixel,
221 1 : pasGCPList[oIter->second].dfGCPLine);
222 : }
223 : else
224 : {
225 2405 : oMapXYToIdx[std::pair<double, double>(afXY[0], afXY[1])] = iGCP;
226 : }
227 :
228 2406 : bool bOK = true;
229 2406 : if (bReversed)
230 : {
231 0 : bOK &= psInfo->poReverse->add_point(afPL[0], afPL[1], afXY);
232 0 : bOK &= psInfo->poForward->add_point(afXY[0], afXY[1], afPL);
233 : }
234 : else
235 : {
236 2406 : bOK &= psInfo->poForward->add_point(afPL[0], afPL[1], afXY);
237 2406 : bOK &= psInfo->poReverse->add_point(afXY[0], afXY[1], afPL);
238 : }
239 2406 : if (!bOK)
240 : {
241 0 : GDALDestroyTPSTransformer(psInfo);
242 0 : return nullptr;
243 : }
244 : }
245 :
246 19 : psInfo->nRefCount = 1;
247 :
248 19 : psInfo->dfSrcApproxErrorReverse = CPLAtof(
249 : CSLFetchNameValueDef(papszOptions, "SRC_APPROX_ERROR_IN_PIXEL", "0"));
250 :
251 19 : int nThreads = 1;
252 19 : if (nGCPCount > 100)
253 : {
254 : const char *pszWarpThreads =
255 2 : CSLFetchNameValue(papszOptions, "NUM_THREADS");
256 2 : if (pszWarpThreads == nullptr)
257 2 : pszWarpThreads = CPLGetConfigOption("GDAL_NUM_THREADS", "1");
258 2 : if (EQUAL(pszWarpThreads, "ALL_CPUS"))
259 0 : nThreads = CPLGetNumCPUs();
260 : else
261 2 : nThreads = atoi(pszWarpThreads);
262 : }
263 :
264 19 : if (nThreads > 1)
265 : {
266 : // Compute direct and reverse transforms in parallel.
267 : CPLJoinableThread *hThread =
268 0 : CPLCreateJoinableThread(GDALTPSComputeForwardInThread, psInfo);
269 0 : psInfo->bReverseSolved = psInfo->poReverse->solve() != 0;
270 0 : if (hThread != nullptr)
271 0 : CPLJoinThread(hThread);
272 : else
273 0 : psInfo->bForwardSolved = psInfo->poForward->solve() != 0;
274 : }
275 : else
276 : {
277 19 : psInfo->bForwardSolved = psInfo->poForward->solve() != 0;
278 19 : psInfo->bReverseSolved = psInfo->poReverse->solve() != 0;
279 : }
280 :
281 19 : if (!psInfo->bForwardSolved || !psInfo->bReverseSolved)
282 : {
283 2 : GDALDestroyTPSTransformer(psInfo);
284 2 : return nullptr;
285 : }
286 :
287 17 : return psInfo;
288 : }
289 :
290 : /************************************************************************/
291 : /* GDALDestroyTPSTransformer() */
292 : /************************************************************************/
293 :
294 : /**
295 : * Destroy TPS transformer.
296 : *
297 : * This function is used to destroy information about a GCP based
298 : * polynomial transformation created with GDALCreateTPSTransformer().
299 : *
300 : * @param pTransformArg the transform arg previously returned by
301 : * GDALCreateTPSTransformer().
302 : */
303 :
304 19 : void GDALDestroyTPSTransformer(void *pTransformArg)
305 :
306 : {
307 19 : if (pTransformArg == nullptr)
308 0 : return;
309 :
310 19 : TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(pTransformArg);
311 :
312 19 : if (CPLAtomicDec(&(psInfo->nRefCount)) == 0)
313 : {
314 19 : delete psInfo->poForward;
315 19 : delete psInfo->poReverse;
316 :
317 19 : delete psInfo;
318 : }
319 : }
320 :
321 : /************************************************************************/
322 : /* GDALTPSTransform() */
323 : /************************************************************************/
324 :
325 : /**
326 : * Transforms point based on GCP derived polynomial model.
327 : *
328 : * This function matches the GDALTransformerFunc signature, and can be
329 : * used to transform one or more points from pixel/line coordinates to
330 : * georeferenced coordinates (SrcToDst) or vice versa (DstToSrc).
331 : *
332 : * @param pTransformArg return value from GDALCreateTPSTransformer().
333 : * @param bDstToSrc TRUE if transformation is from the destination
334 : * (georeferenced) coordinates to pixel/line or FALSE when transforming
335 : * from pixel/line to georeferenced coordinates.
336 : * @param nPointCount the number of values in the x, y and z arrays.
337 : * @param x array containing the X values to be transformed.
338 : * @param y array containing the Y values to be transformed.
339 : * @param z array containing the Z values to be transformed.
340 : * @param panSuccess array in which a flag indicating success (TRUE) or
341 : * failure (FALSE) of the transformation are placed.
342 : *
343 : * @return TRUE.
344 : */
345 :
346 34790 : int GDALTPSTransform(void *pTransformArg, int bDstToSrc, int nPointCount,
347 : double *x, double *y, CPL_UNUSED double *z,
348 : int *panSuccess)
349 : {
350 34790 : VALIDATE_POINTER1(pTransformArg, "GDALTPSTransform", 0);
351 :
352 34790 : TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(pTransformArg);
353 :
354 112088 : for (int i = 0; i < nPointCount; i++)
355 : {
356 77298 : double xy_out[2] = {0.0, 0.0};
357 :
358 77298 : if (bDstToSrc)
359 : {
360 : // Compute initial guess
361 70066 : psInfo->poReverse->get_point(x[i], y[i], xy_out);
362 :
363 461315 : const auto ForwardTransformer = [](double xIn, double yIn,
364 : double &xOut, double &yOut,
365 : void *pUserData)
366 : {
367 461315 : double xyOut[2] = {0.0, 0.0};
368 461315 : TPSTransformInfo *l_psInfo =
369 : static_cast<TPSTransformInfo *>(pUserData);
370 461315 : l_psInfo->poForward->get_point(xIn, yIn, xyOut);
371 461315 : xOut = xyOut[0];
372 461315 : yOut = xyOut[1];
373 461315 : return true;
374 : };
375 :
376 : // Refine the initial guess
377 70066 : GDALGenericInverse2D(
378 70066 : x[i], y[i], xy_out[0], xy_out[1], ForwardTransformer, psInfo,
379 : xy_out[0], xy_out[1],
380 : /* computeJacobianMatrixOnlyAtFirstIter = */ true,
381 : /* toleranceOnOutputCoordinates = */ 0,
382 : psInfo->dfSrcApproxErrorReverse);
383 70066 : x[i] = xy_out[0];
384 70066 : y[i] = xy_out[1];
385 : }
386 : else
387 : {
388 7232 : psInfo->poForward->get_point(x[i], y[i], xy_out);
389 7232 : x[i] = xy_out[0];
390 7232 : y[i] = xy_out[1];
391 : }
392 77298 : panSuccess[i] = TRUE;
393 : }
394 :
395 34790 : return TRUE;
396 : }
397 :
398 : /************************************************************************/
399 : /* GDALSerializeTPSTransformer() */
400 : /************************************************************************/
401 :
402 2 : CPLXMLNode *GDALSerializeTPSTransformer(void *pTransformArg)
403 :
404 : {
405 2 : VALIDATE_POINTER1(pTransformArg, "GDALSerializeTPSTransformer", nullptr);
406 :
407 2 : TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(pTransformArg);
408 :
409 : CPLXMLNode *psTree =
410 2 : CPLCreateXMLNode(nullptr, CXT_Element, "TPSTransformer");
411 :
412 : /* -------------------------------------------------------------------- */
413 : /* Serialize bReversed. */
414 : /* -------------------------------------------------------------------- */
415 2 : CPLCreateXMLElementAndValue(
416 : psTree, "Reversed",
417 4 : CPLString().Printf("%d", static_cast<int>(psInfo->bReversed)));
418 :
419 : /* -------------------------------------------------------------------- */
420 : /* Attach GCP List. */
421 : /* -------------------------------------------------------------------- */
422 2 : if (!psInfo->asGCPs.empty())
423 : {
424 2 : GDALSerializeGCPListToXML(psTree, psInfo->asGCPs, nullptr);
425 : }
426 :
427 2 : if (psInfo->dfSrcApproxErrorReverse > 0)
428 : {
429 2 : CPLCreateXMLElementAndValue(
430 : psTree, "SrcApproxErrorInPixel",
431 4 : CPLString().Printf("%g", psInfo->dfSrcApproxErrorReverse));
432 : }
433 :
434 2 : return psTree;
435 : }
436 :
437 : /************************************************************************/
438 : /* GDALDeserializeTPSTransformer() */
439 : /************************************************************************/
440 :
441 3 : void *GDALDeserializeTPSTransformer(CPLXMLNode *psTree)
442 :
443 : {
444 : /* -------------------------------------------------------------------- */
445 : /* Check for GCPs. */
446 : /* -------------------------------------------------------------------- */
447 3 : CPLXMLNode *psGCPList = CPLGetXMLNode(psTree, "GCPList");
448 :
449 6 : std::vector<gdal::GCP> asGCPs;
450 3 : if (psGCPList != nullptr)
451 : {
452 3 : GDALDeserializeGCPListFromXML(psGCPList, asGCPs, nullptr);
453 : }
454 :
455 : /* -------------------------------------------------------------------- */
456 : /* Get other flags. */
457 : /* -------------------------------------------------------------------- */
458 3 : const int bReversed = atoi(CPLGetXMLValue(psTree, "Reversed", "0"));
459 :
460 3 : CPLStringList aosOptions;
461 : aosOptions.SetNameValue(
462 : "SRC_APPROX_ERROR_IN_PIXEL",
463 3 : CPLGetXMLValue(psTree, "SrcApproxErrorInPixel", nullptr));
464 :
465 : /* -------------------------------------------------------------------- */
466 : /* Generate transformation. */
467 : /* -------------------------------------------------------------------- */
468 3 : void *pResult = GDALCreateTPSTransformerInt(static_cast<int>(asGCPs.size()),
469 : gdal::GCP::c_ptr(asGCPs),
470 : bReversed, aosOptions.List());
471 :
472 6 : return pResult;
473 : }
|
