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
49 CPL_CVSID("$Id: gdal_tps.cpp 4e4327d7c0600596ce2a09bcce1784a55fa7c4e4 2020-09-07 21:50:04 +0200 Even Rouault $")
50
51 CPL_C_START
52 CPLXMLNode *GDALSerializeTPSTransformer( void *pTransformArg );
53 void *GDALDeserializeTPSTransformer( CPLXMLNode *psTree );
54 CPL_C_END
55
56 typedef struct
57 {
58 GDALTransformerInfo sTI;
59
60 VizGeorefSpline2D *poForward;
61 VizGeorefSpline2D *poReverse;
62 bool bForwardSolved;
63 bool bReverseSolved;
64
65 bool bReversed;
66
67 int nGCPCount;
68 GDAL_GCP *pasGCPList;
69
70 volatile int nRefCount;
71
72 } TPSTransformInfo;
73
74 /************************************************************************/
75 /* GDALCreateSimilarTPSTransformer() */
76 /************************************************************************/
77
78 static
GDALCreateSimilarTPSTransformer(void * hTransformArg,double dfRatioX,double dfRatioY)79 void* GDALCreateSimilarTPSTransformer( void *hTransformArg,
80 double dfRatioX, double dfRatioY )
81 {
82 VALIDATE_POINTER1( hTransformArg, "GDALCreateSimilarTPSTransformer", nullptr );
83
84 TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(hTransformArg);
85
86 if( dfRatioX == 1.0 && dfRatioY == 1.0 )
87 {
88 // We can just use a ref count, since using the source transformation
89 // is thread-safe.
90 CPLAtomicInc(&(psInfo->nRefCount));
91 }
92 else
93 {
94 GDAL_GCP *pasGCPList = GDALDuplicateGCPs( psInfo->nGCPCount,
95 psInfo->pasGCPList );
96 for( int i = 0; i < psInfo->nGCPCount; i++ )
97 {
98 pasGCPList[i].dfGCPPixel /= dfRatioX;
99 pasGCPList[i].dfGCPLine /= dfRatioY;
100 }
101 psInfo = static_cast<TPSTransformInfo *>(
102 GDALCreateTPSTransformer( psInfo->nGCPCount, pasGCPList,
103 psInfo->bReversed ));
104 GDALDeinitGCPs( psInfo->nGCPCount, pasGCPList );
105 CPLFree( pasGCPList );
106 }
107
108 return psInfo;
109 }
110
111 /************************************************************************/
112 /* GDALCreateTPSTransformer() */
113 /************************************************************************/
114
115 /**
116 * Create Thin Plate Spline transformer from GCPs.
117 *
118 * The thin plate spline transformer produces exact transformation
119 * at all control points and smoothly varying transformations between
120 * control points with greatest influence from local control points.
121 * It is suitable for for many applications not well modeled by polynomial
122 * transformations.
123 *
124 * Creating the TPS transformer involves solving systems of linear equations
125 * related to the number of control points involved. This solution is
126 * computed within this function call. It can be quite an expensive operation
127 * for large numbers of GCPs. For instance, for reference, it takes on the
128 * order of 10s for 400 GCPs on a 2GHz Athlon processor.
129 *
130 * TPS Transformers are serializable.
131 *
132 * The GDAL Thin Plate Spline transformer is based on code provided by
133 * Gilad Ronnen on behalf of VIZRT Inc (http://www.visrt.com). Incorporation
134 * of the algorithm into GDAL was supported by the Centro di Ecologia Alpina
135 * (http://www.cealp.it).
136 *
137 * @param nGCPCount the number of GCPs in pasGCPList.
138 * @param pasGCPList an array of GCPs to be used as input.
139 * @param bReversed set it to TRUE to compute the reversed transformation.
140 *
141 * @return the transform argument or NULL if creation fails.
142 */
143
GDALCreateTPSTransformer(int nGCPCount,const GDAL_GCP * pasGCPList,int bReversed)144 void *GDALCreateTPSTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
145 int bReversed )
146 {
147 return GDALCreateTPSTransformerInt(nGCPCount, pasGCPList, bReversed, nullptr);
148 }
149
GDALTPSComputeForwardInThread(void * pData)150 static void GDALTPSComputeForwardInThread( void *pData )
151 {
152 TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(pData);
153 psInfo->bForwardSolved = psInfo->poForward->solve() != 0;
154 }
155
GDALCreateTPSTransformerInt(int nGCPCount,const GDAL_GCP * pasGCPList,int bReversed,char ** papszOptions)156 void *GDALCreateTPSTransformerInt( int nGCPCount, const GDAL_GCP *pasGCPList,
157 int bReversed, char** papszOptions )
158
159 {
160 /* -------------------------------------------------------------------- */
161 /* Allocate transform info. */
162 /* -------------------------------------------------------------------- */
163 TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(
164 CPLCalloc(sizeof(TPSTransformInfo), 1));
165
166 psInfo->pasGCPList = GDALDuplicateGCPs( nGCPCount, pasGCPList );
167 psInfo->nGCPCount = nGCPCount;
168
169 psInfo->bReversed = CPL_TO_BOOL(bReversed);
170 psInfo->poForward = new VizGeorefSpline2D( 2 );
171 psInfo->poReverse = new VizGeorefSpline2D( 2 );
172
173 memcpy( psInfo->sTI.abySignature,
174 GDAL_GTI2_SIGNATURE,
175 strlen(GDAL_GTI2_SIGNATURE) );
176 psInfo->sTI.pszClassName = "GDALTPSTransformer";
177 psInfo->sTI.pfnTransform = GDALTPSTransform;
178 psInfo->sTI.pfnCleanup = GDALDestroyTPSTransformer;
179 psInfo->sTI.pfnSerialize = GDALSerializeTPSTransformer;
180 psInfo->sTI.pfnCreateSimilar = GDALCreateSimilarTPSTransformer;
181
182 /* -------------------------------------------------------------------- */
183 /* Attach (non-redundant) points to the transformation. */
184 /* -------------------------------------------------------------------- */
185 std::map< std::pair<double, double>, int > oMapPixelLineToIdx;
186 std::map< std::pair<double, double>, int > oMapXYToIdx;
187 for( int iGCP = 0; iGCP < nGCPCount; iGCP++ )
188 {
189 const double afPL[2] = {
190 pasGCPList[iGCP].dfGCPPixel,
191 pasGCPList[iGCP].dfGCPLine };
192 const double afXY[2] =
193 { pasGCPList[iGCP].dfGCPX, pasGCPList[iGCP].dfGCPY };
194
195 std::map< std::pair<double, double>, int >::iterator oIter(
196 oMapPixelLineToIdx.find(std::pair<double, double>(afPL[0],
197 afPL[1])));
198 if( oIter != oMapPixelLineToIdx.end() )
199 {
200 if( afXY[0] == pasGCPList[oIter->second].dfGCPX &&
201 afXY[1] == pasGCPList[oIter->second].dfGCPY )
202 {
203 continue;
204 }
205 else
206 {
207 CPLError(CE_Warning, CPLE_AppDefined,
208 "GCP %d and %d have same (pixel,line)=(%f,%f), "
209 "but different (X,Y): (%f,%f) vs (%f,%f)",
210 iGCP + 1, oIter->second,
211 afPL[0], afPL[1],
212 afXY[0], afXY[1],
213 pasGCPList[oIter->second].dfGCPX,
214 pasGCPList[oIter->second].dfGCPY);
215 }
216 }
217 else
218 {
219 oMapPixelLineToIdx[std::pair<double, double>(afPL[0], afPL[1])] =
220 iGCP;
221 }
222
223 oIter = oMapXYToIdx.find(std::pair<double, double>(afXY[0], afXY[1]));
224 if( oIter != oMapXYToIdx.end() )
225 {
226 CPLError(CE_Warning, CPLE_AppDefined,
227 "GCP %d and %d have same (x,y)=(%f,%f), "
228 "but different (pixel,line): (%f,%f) vs (%f,%f)",
229 iGCP + 1, oIter->second,
230 afXY[0], afXY[1],
231 afPL[0], afPL[1],
232 pasGCPList[oIter->second].dfGCPPixel,
233 pasGCPList[oIter->second].dfGCPLine);
234 }
235 else
236 {
237 oMapXYToIdx[std::pair<double, double>(afXY[0], afXY[1])] = iGCP;
238 }
239
240 bool bOK = true;
241 if( bReversed )
242 {
243 bOK &= psInfo->poReverse->add_point( afPL[0], afPL[1], afXY );
244 bOK &= psInfo->poForward->add_point( afXY[0], afXY[1], afPL );
245 }
246 else
247 {
248 bOK &= psInfo->poForward->add_point( afPL[0], afPL[1], afXY );
249 bOK &= psInfo->poReverse->add_point( afXY[0], afXY[1], afPL );
250 }
251 if( !bOK )
252 {
253 GDALDestroyTPSTransformer(psInfo);
254 return nullptr;
255 }
256 }
257
258 psInfo->nRefCount = 1;
259
260 int nThreads = 1;
261 if( nGCPCount > 100 )
262 {
263 const char* pszWarpThreads =
264 CSLFetchNameValue(papszOptions, "NUM_THREADS");
265 if( pszWarpThreads == nullptr )
266 pszWarpThreads = CPLGetConfigOption("GDAL_NUM_THREADS", "1");
267 if( EQUAL(pszWarpThreads, "ALL_CPUS") )
268 nThreads = CPLGetNumCPUs();
269 else
270 nThreads = atoi(pszWarpThreads);
271 }
272
273 if( nThreads > 1 )
274 {
275 // Compute direct and reverse transforms in parallel.
276 CPLJoinableThread* hThread =
277 CPLCreateJoinableThread(GDALTPSComputeForwardInThread, psInfo);
278 psInfo->bReverseSolved = psInfo->poReverse->solve() != 0;
279 if( hThread != nullptr )
280 CPLJoinThread(hThread);
281 else
282 psInfo->bForwardSolved = psInfo->poForward->solve() != 0;
283 }
284 else
285 {
286 psInfo->bForwardSolved = psInfo->poForward->solve() != 0;
287 psInfo->bReverseSolved = psInfo->poReverse->solve() != 0;
288 }
289
290 if( !psInfo->bForwardSolved || !psInfo->bReverseSolved )
291 {
292 GDALDestroyTPSTransformer(psInfo);
293 return nullptr;
294 }
295
296 return psInfo;
297 }
298
299 /************************************************************************/
300 /* GDALDestroyTPSTransformer() */
301 /************************************************************************/
302
303 /**
304 * Destroy TPS transformer.
305 *
306 * This function is used to destroy information about a GCP based
307 * polynomial transformation created with GDALCreateTPSTransformer().
308 *
309 * @param pTransformArg the transform arg previously returned by
310 * GDALCreateTPSTransformer().
311 */
312
GDALDestroyTPSTransformer(void * pTransformArg)313 void GDALDestroyTPSTransformer( void *pTransformArg )
314
315 {
316 if( pTransformArg == nullptr )
317 return;
318
319 TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(pTransformArg);
320
321 if( CPLAtomicDec(&(psInfo->nRefCount)) == 0 )
322 {
323 delete psInfo->poForward;
324 delete psInfo->poReverse;
325
326 GDALDeinitGCPs( psInfo->nGCPCount, psInfo->pasGCPList );
327 CPLFree( psInfo->pasGCPList );
328
329 CPLFree( pTransformArg );
330 }
331 }
332
333 /************************************************************************/
334 /* GDALTPSTransform() */
335 /************************************************************************/
336
337 /**
338 * Transforms point based on GCP derived polynomial model.
339 *
340 * This function matches the GDALTransformerFunc signature, and can be
341 * used to transform one or more points from pixel/line coordinates to
342 * georeferenced coordinates (SrcToDst) or vice versa (DstToSrc).
343 *
344 * @param pTransformArg return value from GDALCreateTPSTransformer().
345 * @param bDstToSrc TRUE if transformation is from the destination
346 * (georeferenced) coordinates to pixel/line or FALSE when transforming
347 * from pixel/line to georeferenced coordinates.
348 * @param nPointCount the number of values in the x, y and z arrays.
349 * @param x array containing the X values to be transformed.
350 * @param y array containing the Y values to be transformed.
351 * @param z array containing the Z values to be transformed.
352 * @param panSuccess array in which a flag indicating success (TRUE) or
353 * failure (FALSE) of the transformation are placed.
354 *
355 * @return TRUE.
356 */
357
GDALTPSTransform(void * pTransformArg,int bDstToSrc,int nPointCount,double * x,double * y,CPL_UNUSED double * z,int * panSuccess)358 int GDALTPSTransform( void *pTransformArg, int bDstToSrc,
359 int nPointCount,
360 double *x, double *y,
361 CPL_UNUSED double *z,
362 int *panSuccess )
363 {
364 VALIDATE_POINTER1( pTransformArg, "GDALTPSTransform", 0 );
365
366 TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(pTransformArg);
367
368 for( int i = 0; i < nPointCount; i++ )
369 {
370 double xy_out[2] = { 0.0, 0.0 };
371
372 if( bDstToSrc )
373 {
374 psInfo->poReverse->get_point( x[i], y[i], xy_out );
375 x[i] = xy_out[0];
376 y[i] = xy_out[1];
377 }
378 else
379 {
380 psInfo->poForward->get_point( x[i], y[i], xy_out );
381 x[i] = xy_out[0];
382 y[i] = xy_out[1];
383 }
384 panSuccess[i] = TRUE;
385 }
386
387 return TRUE;
388 }
389
390 /************************************************************************/
391 /* GDALSerializeTPSTransformer() */
392 /************************************************************************/
393
GDALSerializeTPSTransformer(void * pTransformArg)394 CPLXMLNode *GDALSerializeTPSTransformer( void *pTransformArg )
395
396 {
397 VALIDATE_POINTER1( pTransformArg, "GDALSerializeTPSTransformer", nullptr );
398
399 TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(pTransformArg);
400
401 CPLXMLNode *psTree = CPLCreateXMLNode(nullptr, CXT_Element, "TPSTransformer");
402
403 /* -------------------------------------------------------------------- */
404 /* Serialize bReversed. */
405 /* -------------------------------------------------------------------- */
406 CPLCreateXMLElementAndValue(
407 psTree, "Reversed",
408 CPLString().Printf( "%d", static_cast<int>(psInfo->bReversed) ) );
409
410 /* -------------------------------------------------------------------- */
411 /* Attach GCP List. */
412 /* -------------------------------------------------------------------- */
413 if( psInfo->nGCPCount > 0 )
414 {
415 GDALSerializeGCPListToXML( psTree,
416 psInfo->pasGCPList,
417 psInfo->nGCPCount,
418 nullptr );
419 }
420
421 return psTree;
422 }
423
424 /************************************************************************/
425 /* GDALDeserializeTPSTransformer() */
426 /************************************************************************/
427
GDALDeserializeTPSTransformer(CPLXMLNode * psTree)428 void *GDALDeserializeTPSTransformer( CPLXMLNode *psTree )
429
430 {
431 /* -------------------------------------------------------------------- */
432 /* Check for GCPs. */
433 /* -------------------------------------------------------------------- */
434 CPLXMLNode *psGCPList = CPLGetXMLNode( psTree, "GCPList" );
435 GDAL_GCP *pasGCPList = nullptr;
436 int nGCPCount = 0;
437
438 if( psGCPList != nullptr )
439 {
440 GDALDeserializeGCPListFromXML( psGCPList,
441 &pasGCPList,
442 &nGCPCount,
443 nullptr );
444 }
445
446 /* -------------------------------------------------------------------- */
447 /* Get other flags. */
448 /* -------------------------------------------------------------------- */
449 const int bReversed = atoi(CPLGetXMLValue(psTree, "Reversed", "0"));
450
451 /* -------------------------------------------------------------------- */
452 /* Generate transformation. */
453 /* -------------------------------------------------------------------- */
454 void *pResult =
455 GDALCreateTPSTransformer( nGCPCount, pasGCPList, bReversed );
456
457 /* -------------------------------------------------------------------- */
458 /* Cleanup GCP copy. */
459 /* -------------------------------------------------------------------- */
460 GDALDeinitGCPs( nGCPCount, pasGCPList );
461 CPLFree( pasGCPList );
462
463 return pResult;
464 }
465