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