1 : /******************************************************************************
2 : * $Id: gdal_tps.cpp 25304 2012-12-13 21:03:58Z rouault $
3 : *
4 : * Project: High Performance Image Reprojector
5 : * Purpose: Thin Plate Spline transformer (GDAL wrapper portion)
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2004, Frank Warmerdam <warmerdam@pobox.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 "thinplatespline.h"
31 : #include "gdal_alg.h"
32 : #include "cpl_conv.h"
33 : #include "cpl_string.h"
34 : #include "cpl_atomic_ops.h"
35 :
36 : CPL_CVSID("$Id: gdal_tps.cpp 25304 2012-12-13 21:03:58Z rouault $");
37 :
38 : CPL_C_START
39 : CPLXMLNode *GDALSerializeTPSTransformer( void *pTransformArg );
40 : void *GDALDeserializeTPSTransformer( CPLXMLNode *psTree );
41 : CPL_C_END
42 :
43 : typedef struct
44 : {
45 : GDALTransformerInfo sTI;
46 :
47 : VizGeorefSpline2D *poForward;
48 : VizGeorefSpline2D *poReverse;
49 :
50 : int bReversed;
51 :
52 : int nGCPCount;
53 : GDAL_GCP *pasGCPList;
54 :
55 : volatile int nRefCount;
56 :
57 : } TPSTransformInfo;
58 :
59 : /************************************************************************/
60 : /* GDALCloneTPSTransformer() */
61 : /************************************************************************/
62 :
63 0 : void* GDALCloneTPSTransformer( void *hTransformArg )
64 : {
65 0 : VALIDATE_POINTER1( hTransformArg, "GDALCloneTPSTransformer", NULL );
66 :
67 : TPSTransformInfo *psInfo =
68 0 : (TPSTransformInfo *) hTransformArg;
69 :
70 : /* We can just use a ref count, since using the source transformation */
71 : /* is thread-safe */
72 0 : CPLAtomicInc(&(psInfo->nRefCount));
73 :
74 0 : return psInfo;
75 : }
76 :
77 : /************************************************************************/
78 : /* GDALCreateTPSTransformer() */
79 : /************************************************************************/
80 :
81 : /**
82 : * Create Thin Plate Spline transformer from GCPs.
83 : *
84 : * The thin plate spline transformer produces exact transformation
85 : * at all control points and smoothly varying transformations between
86 : * control points with greatest influence from local control points.
87 : * It is suitable for for many applications not well modelled by polynomial
88 : * transformations.
89 : *
90 : * Creating the TPS transformer involves solving systems of linear equations
91 : * related to the number of control points involved. This solution is
92 : * computed within this function call. It can be quite an expensive operation
93 : * for large numbers of GCPs. For instance, for reference, it takes on the
94 : * order of 10s for 400 GCPs on a 2GHz Athlon processor.
95 : *
96 : * TPS Transformers are serializable.
97 : *
98 : * The GDAL Thin Plate Spline transformer is based on code provided by
99 : * Gilad Ronnen on behalf of VIZRT Inc (http://www.visrt.com). Incorporation
100 : * of the algorithm into GDAL was supported by the Centro di Ecologia Alpina
101 : * (http://www.cealp.it).
102 : *
103 : * @param nGCPCount the number of GCPs in pasGCPList.
104 : * @param pasGCPList an array of GCPs to be used as input.
105 : * @param bReversed set it to TRUE to compute the reversed transformation.
106 : *
107 : * @return the transform argument or NULL if creation fails.
108 : */
109 :
110 6 : void *GDALCreateTPSTransformer( int nGCPCount, const GDAL_GCP *pasGCPList,
111 : int bReversed )
112 :
113 : {
114 : TPSTransformInfo *psInfo;
115 : int iGCP;
116 :
117 : /* -------------------------------------------------------------------- */
118 : /* Allocate transform info. */
119 : /* -------------------------------------------------------------------- */
120 6 : psInfo = (TPSTransformInfo *) CPLCalloc(sizeof(TPSTransformInfo),1);
121 :
122 6 : psInfo->pasGCPList = GDALDuplicateGCPs( nGCPCount, pasGCPList );
123 6 : psInfo->nGCPCount = nGCPCount;
124 :
125 6 : psInfo->bReversed = bReversed;
126 6 : psInfo->poForward = new VizGeorefSpline2D( 2 );
127 12 : psInfo->poReverse = new VizGeorefSpline2D( 2 );
128 :
129 6 : strcpy( psInfo->sTI.szSignature, "GTI" );
130 6 : psInfo->sTI.pszClassName = "GDALTPSTransformer";
131 6 : psInfo->sTI.pfnTransform = GDALTPSTransform;
132 6 : psInfo->sTI.pfnCleanup = GDALDestroyTPSTransformer;
133 6 : psInfo->sTI.pfnSerialize = GDALSerializeTPSTransformer;
134 :
135 : /* -------------------------------------------------------------------- */
136 : /* Attach all the points to the transformation. */
137 : /* -------------------------------------------------------------------- */
138 29 : for( iGCP = 0; iGCP < nGCPCount; iGCP++ )
139 : {
140 : double afPL[2], afXY[2];
141 :
142 23 : afPL[0] = pasGCPList[iGCP].dfGCPPixel;
143 23 : afPL[1] = pasGCPList[iGCP].dfGCPLine;
144 23 : afXY[0] = pasGCPList[iGCP].dfGCPX;
145 23 : afXY[1] = pasGCPList[iGCP].dfGCPY;
146 :
147 23 : if( bReversed )
148 : {
149 0 : psInfo->poReverse->add_point( afPL[0], afPL[1], afXY );
150 0 : psInfo->poForward->add_point( afXY[0], afXY[1], afPL );
151 : }
152 : else
153 : {
154 23 : psInfo->poForward->add_point( afPL[0], afPL[1], afXY );
155 23 : psInfo->poReverse->add_point( afXY[0], afXY[1], afPL );
156 : }
157 : }
158 :
159 6 : psInfo->nRefCount = 1;
160 :
161 6 : psInfo->poForward->solve();
162 6 : psInfo->poReverse->solve();
163 :
164 6 : return psInfo;
165 : }
166 :
167 : /************************************************************************/
168 : /* GDALDestroyTPSTransformer() */
169 : /************************************************************************/
170 :
171 : /**
172 : * Destroy TPS transformer.
173 : *
174 : * This function is used to destroy information about a GCP based
175 : * polynomial transformation created with GDALCreateTPSTransformer().
176 : *
177 : * @param pTransformArg the transform arg previously returned by
178 : * GDALCreateTPSTransformer().
179 : */
180 :
181 6 : void GDALDestroyTPSTransformer( void *pTransformArg )
182 :
183 : {
184 6 : VALIDATE_POINTER0( pTransformArg, "GDALDestroyTPSTransformer" );
185 :
186 6 : TPSTransformInfo *psInfo = (TPSTransformInfo *) pTransformArg;
187 :
188 6 : if( CPLAtomicDec(&(psInfo->nRefCount)) == 0 )
189 : {
190 6 : delete psInfo->poForward;
191 6 : delete psInfo->poReverse;
192 :
193 6 : GDALDeinitGCPs( psInfo->nGCPCount, psInfo->pasGCPList );
194 6 : CPLFree( psInfo->pasGCPList );
195 :
196 6 : CPLFree( pTransformArg );
197 : }
198 : }
199 :
200 : /************************************************************************/
201 : /* GDALTPSTransform() */
202 : /************************************************************************/
203 :
204 : /**
205 : * Transforms point based on GCP derived polynomial model.
206 : *
207 : * This function matches the GDALTransformerFunc signature, and can be
208 : * used to transform one or more points from pixel/line coordinates to
209 : * georeferenced coordinates (SrcToDst) or vice versa (DstToSrc).
210 : *
211 : * @param pTransformArg return value from GDALCreateTPSTransformer().
212 : * @param bDstToSrc TRUE if transformation is from the destination
213 : * (georeferenced) coordinates to pixel/line or FALSE when transforming
214 : * from pixel/line to georeferenced coordinates.
215 : * @param nPointCount the number of values in the x, y and z arrays.
216 : * @param x array containing the X values to be transformed.
217 : * @param y array containing the Y values to be transformed.
218 : * @param z array containing the Z values to be transformed.
219 : * @param panSuccess array in which a flag indicating success (TRUE) or
220 : * failure (FALSE) of the transformation are placed.
221 : *
222 : * @return TRUE.
223 : */
224 :
225 1826 : int GDALTPSTransform( void *pTransformArg, int bDstToSrc,
226 : int nPointCount,
227 : double *x, double *y, double *z,
228 : int *panSuccess )
229 :
230 : {
231 1826 : VALIDATE_POINTER1( pTransformArg, "GDALTPSTransform", 0 );
232 :
233 : int i;
234 1826 : TPSTransformInfo *psInfo = (TPSTransformInfo *) pTransformArg;
235 :
236 4390 : for( i = 0; i < nPointCount; i++ )
237 : {
238 : double xy_out[2];
239 :
240 2564 : if( bDstToSrc )
241 : {
242 1423 : psInfo->poReverse->get_point( x[i], y[i], xy_out );
243 1423 : x[i] = xy_out[0];
244 1423 : y[i] = xy_out[1];
245 : }
246 : else
247 : {
248 1141 : psInfo->poForward->get_point( x[i], y[i], xy_out );
249 1141 : x[i] = xy_out[0];
250 1141 : y[i] = xy_out[1];
251 : }
252 2564 : panSuccess[i] = TRUE;
253 : }
254 :
255 1826 : return TRUE;
256 : }
257 :
258 : /************************************************************************/
259 : /* GDALSerializeTPSTransformer() */
260 : /************************************************************************/
261 :
262 1 : CPLXMLNode *GDALSerializeTPSTransformer( void *pTransformArg )
263 :
264 : {
265 1 : VALIDATE_POINTER1( pTransformArg, "GDALSerializeTPSTransformer", NULL );
266 :
267 : CPLXMLNode *psTree;
268 1 : TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(pTransformArg);
269 :
270 1 : psTree = CPLCreateXMLNode( NULL, CXT_Element, "TPSTransformer" );
271 :
272 : /* -------------------------------------------------------------------- */
273 : /* Serialize bReversed. */
274 : /* -------------------------------------------------------------------- */
275 : CPLCreateXMLElementAndValue(
276 : psTree, "Reversed",
277 1 : CPLString().Printf( "%d", psInfo->bReversed ) );
278 :
279 : /* -------------------------------------------------------------------- */
280 : /* Attach GCP List. */
281 : /* -------------------------------------------------------------------- */
282 1 : if( psInfo->nGCPCount > 0 )
283 : {
284 : int iGCP;
285 : CPLXMLNode *psGCPList = CPLCreateXMLNode( psTree, CXT_Element,
286 1 : "GCPList" );
287 :
288 5 : for( iGCP = 0; iGCP < psInfo->nGCPCount; iGCP++ )
289 : {
290 : CPLXMLNode *psXMLGCP;
291 4 : GDAL_GCP *psGCP = psInfo->pasGCPList + iGCP;
292 :
293 4 : psXMLGCP = CPLCreateXMLNode( psGCPList, CXT_Element, "GCP" );
294 :
295 4 : CPLSetXMLValue( psXMLGCP, "#Id", psGCP->pszId );
296 :
297 4 : if( psGCP->pszInfo != NULL && strlen(psGCP->pszInfo) > 0 )
298 0 : CPLSetXMLValue( psXMLGCP, "Info", psGCP->pszInfo );
299 :
300 : CPLSetXMLValue( psXMLGCP, "#Pixel",
301 4 : CPLString().Printf( "%.4f", psGCP->dfGCPPixel ) );
302 :
303 : CPLSetXMLValue( psXMLGCP, "#Line",
304 8 : CPLString().Printf( "%.4f", psGCP->dfGCPLine ) );
305 :
306 : CPLSetXMLValue( psXMLGCP, "#X",
307 8 : CPLString().Printf( "%.12E", psGCP->dfGCPX ) );
308 :
309 : CPLSetXMLValue( psXMLGCP, "#Y",
310 8 : CPLString().Printf( "%.12E", psGCP->dfGCPY ) );
311 :
312 4 : if( psGCP->dfGCPZ != 0.0 )
313 : CPLSetXMLValue( psXMLGCP, "#GCPZ",
314 0 : CPLString().Printf( "%.12E", psGCP->dfGCPZ ) );
315 : }
316 : }
317 :
318 1 : return psTree;
319 : }
320 :
321 : /************************************************************************/
322 : /* GDALDeserializeTPSTransformer() */
323 : /************************************************************************/
324 :
325 1 : void *GDALDeserializeTPSTransformer( CPLXMLNode *psTree )
326 :
327 : {
328 1 : GDAL_GCP *pasGCPList = 0;
329 1 : int nGCPCount = 0;
330 : void *pResult;
331 : int bReversed;
332 :
333 : /* -------------------------------------------------------------------- */
334 : /* Check for GCPs. */
335 : /* -------------------------------------------------------------------- */
336 1 : CPLXMLNode *psGCPList = CPLGetXMLNode( psTree, "GCPList" );
337 :
338 1 : if( psGCPList != NULL )
339 : {
340 1 : int nGCPMax = 0;
341 : CPLXMLNode *psXMLGCP;
342 :
343 : // Count GCPs.
344 5 : for( psXMLGCP = psGCPList->psChild; psXMLGCP != NULL;
345 : psXMLGCP = psXMLGCP->psNext )
346 4 : nGCPMax++;
347 :
348 1 : pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),nGCPMax);
349 :
350 5 : for( psXMLGCP = psGCPList->psChild; psXMLGCP != NULL;
351 : psXMLGCP = psXMLGCP->psNext )
352 : {
353 4 : GDAL_GCP *psGCP = pasGCPList + nGCPCount;
354 :
355 4 : if( !EQUAL(psXMLGCP->pszValue,"GCP") ||
356 : psXMLGCP->eType != CXT_Element )
357 0 : continue;
358 :
359 4 : GDALInitGCPs( 1, psGCP );
360 :
361 4 : CPLFree( psGCP->pszId );
362 4 : psGCP->pszId = CPLStrdup(CPLGetXMLValue(psXMLGCP,"Id",""));
363 :
364 4 : CPLFree( psGCP->pszInfo );
365 4 : psGCP->pszInfo = CPLStrdup(CPLGetXMLValue(psXMLGCP,"Info",""));
366 :
367 4 : psGCP->dfGCPPixel = atof(CPLGetXMLValue(psXMLGCP,"Pixel","0.0"));
368 4 : psGCP->dfGCPLine = atof(CPLGetXMLValue(psXMLGCP,"Line","0.0"));
369 :
370 4 : psGCP->dfGCPX = atof(CPLGetXMLValue(psXMLGCP,"X","0.0"));
371 4 : psGCP->dfGCPY = atof(CPLGetXMLValue(psXMLGCP,"Y","0.0"));
372 4 : psGCP->dfGCPZ = atof(CPLGetXMLValue(psXMLGCP,"Z","0.0"));
373 4 : nGCPCount++;
374 : }
375 : }
376 :
377 : /* -------------------------------------------------------------------- */
378 : /* Get other flags. */
379 : /* -------------------------------------------------------------------- */
380 1 : bReversed = atoi(CPLGetXMLValue(psTree,"Reversed","0"));
381 :
382 : /* -------------------------------------------------------------------- */
383 : /* Generate transformation. */
384 : /* -------------------------------------------------------------------- */
385 1 : pResult = GDALCreateTPSTransformer( nGCPCount, pasGCPList, bReversed );
386 :
387 : /* -------------------------------------------------------------------- */
388 : /* Cleanup GCP copy. */
389 : /* -------------------------------------------------------------------- */
390 1 : GDALDeinitGCPs( nGCPCount, pasGCPList );
391 1 : CPLFree( pasGCPList );
392 :
393 1 : return pResult;
394 : }
|
