1 : /******************************************************************************
2 : * $Id: gdaltransformer.cpp 23156 2011-10-01 15:34:16Z rouault $
3 : *
4 : * Project: Mapinfo Image Warper
5 : * Purpose: Implementation of one or more GDALTrasformerFunc types, including
6 : * the GenImgProj (general image reprojector) transformer.
7 : * Author: Frank Warmerdam, warmerdam@pobox.com
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 2002, i3 - information integration and imaging
11 : * Fort Collin, CO
12 : *
13 : * Permission is hereby granted, free of charge, to any person obtaining a
14 : * copy of this software and associated documentation files (the "Software"),
15 : * to deal in the Software without restriction, including without limitation
16 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 : * and/or sell copies of the Software, and to permit persons to whom the
18 : * Software is furnished to do so, subject to the following conditions:
19 : *
20 : * The above copyright notice and this permission notice shall be included
21 : * in all copies or substantial portions of the Software.
22 : *
23 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 : * DEALINGS IN THE SOFTWARE.
30 : ****************************************************************************/
31 :
32 : #include "gdal_priv.h"
33 : #include "gdal_alg.h"
34 : #include "ogr_spatialref.h"
35 : #include "cpl_string.h"
36 : #include "gdal_alg_priv.h"
37 : #include "cpl_list.h"
38 : #include "cpl_multiproc.h"
39 :
40 : CPL_CVSID("$Id: gdaltransformer.cpp 23156 2011-10-01 15:34:16Z rouault $");
41 : CPL_C_START
42 : void *GDALDeserializeGCPTransformer( CPLXMLNode *psTree );
43 : void *GDALDeserializeTPSTransformer( CPLXMLNode *psTree );
44 : void *GDALDeserializeGeoLocTransformer( CPLXMLNode *psTree );
45 : void *GDALDeserializeRPCTransformer( CPLXMLNode *psTree );
46 : CPL_C_END
47 :
48 : static CPLXMLNode *GDALSerializeReprojectionTransformer( void *pTransformArg );
49 : static void *GDALDeserializeReprojectionTransformer( CPLXMLNode *psTree );
50 :
51 : static CPLXMLNode *GDALSerializeGenImgProjTransformer( void *pTransformArg );
52 : static void *GDALDeserializeGenImgProjTransformer( CPLXMLNode *psTree );
53 :
54 : static void GDALRefreshGenImgProjTransformer(void* hTransformArg);
55 :
56 : /************************************************************************/
57 : /* GDALTransformFunc */
58 : /* */
59 : /* Documentation for GDALTransformFunc typedef. */
60 : /************************************************************************/
61 :
62 : /*!
63 :
64 : \typedef int GDALTransformerFunc
65 :
66 : Generic signature for spatial point transformers.
67 :
68 : This function signature is used for a variety of functions that accept
69 : passed in functions used to transform point locations between two coordinate
70 : spaces.
71 :
72 : The GDALCreateGenImgProjTransformer(), GDALCreateReprojectionTransformer(),
73 : GDALCreateGCPTransformer() and GDALCreateApproxTransformer() functions can
74 : be used to prepare argument data for some built-in transformers. As well,
75 : applications can implement their own transformers to the following signature.
76 :
77 : \code
78 : typedef int
79 : (*GDALTransformerFunc)( void *pTransformerArg,
80 : int bDstToSrc, int nPointCount,
81 : double *x, double *y, double *z, int *panSuccess );
82 : \endcode
83 :
84 : @param pTransformerArg application supplied callback data used by the
85 : transformer.
86 :
87 : @param bDstToSrc if TRUE the transformation will be from the destination
88 : coordinate space to the source coordinate system, otherwise the transformation
89 : will be from the source coordinate system to the destination coordinate system.
90 :
91 : @param nPointCount number of points in the x, y and z arrays.
92 :
93 : @param x input X coordinates. Results returned in same array.
94 :
95 : @param y input Y coordinates. Results returned in same array.
96 :
97 : @param z input Z coordinates. Results returned in same array.
98 :
99 : @param panSuccess array of ints in which success (TRUE) or failure (FALSE)
100 : flags are returned for the translation of each point.
101 :
102 : @return TRUE if the overall transformation succeeds (though some individual
103 : points may have failed) or FALSE if the overall transformation fails.
104 :
105 : */
106 :
107 : /************************************************************************/
108 : /* GDALSuggestedWarpOutput() */
109 : /************************************************************************/
110 :
111 : /**
112 : * Suggest output file size.
113 : *
114 : * This function is used to suggest the size, and georeferenced extents
115 : * appropriate given the indicated transformation and input file. It walks
116 : * the edges of the input file (approximately 20 sample points along each
117 : * edge) transforming into output coordinates in order to get an extents box.
118 : *
119 : * Then a resolution is computed with the intent that the length of the
120 : * distance from the top left corner of the output imagery to the bottom right
121 : * corner would represent the same number of pixels as in the source image.
122 : * Note that if the image is somewhat rotated the diagonal taken isnt of the
123 : * whole output bounding rectangle, but instead of the locations where the
124 : * top/left and bottom/right corners transform. The output pixel size is
125 : * always square. This is intended to approximately preserve the resolution
126 : * of the input data in the output file.
127 : *
128 : * The values returned in padfGeoTransformOut, pnPixels and pnLines are
129 : * the suggested number of pixels and lines for the output file, and the
130 : * geotransform relating those pixels to the output georeferenced coordinates.
131 : *
132 : * The trickiest part of using the function is ensuring that the
133 : * transformer created is from source file pixel/line coordinates to
134 : * output file georeferenced coordinates. This can be accomplished with
135 : * GDALCreateGenImgProjTransformer() by passing a NULL for the hDstDS.
136 : *
137 : * @param hSrcDS the input image (it is assumed the whole input images is
138 : * being transformed).
139 : * @param pfnTransformer the transformer function.
140 : * @param pTransformArg the callback data for the transformer function.
141 : * @param padfGeoTransformOut the array of six doubles in which the suggested
142 : * geotransform is returned.
143 : * @param pnPixels int in which the suggest pixel width of output is returned.
144 : * @param pnLines int in which the suggest pixel height of output is returned.
145 : *
146 : * @return CE_None if successful or CE_Failure otherwise.
147 : */
148 :
149 :
150 : CPLErr CPL_STDCALL
151 26 : GDALSuggestedWarpOutput( GDALDatasetH hSrcDS,
152 : GDALTransformerFunc pfnTransformer,
153 : void *pTransformArg,
154 : double *padfGeoTransformOut,
155 : int *pnPixels, int *pnLines )
156 :
157 : {
158 26 : VALIDATE_POINTER1( hSrcDS, "GDALSuggestedWarpOutput", CE_Failure );
159 :
160 26 : double adfExtent[4] = { 0 };
161 :
162 : return GDALSuggestedWarpOutput2( hSrcDS, pfnTransformer, pTransformArg,
163 : padfGeoTransformOut, pnPixels, pnLines,
164 26 : adfExtent, 0 );
165 : }
166 :
167 :
168 634 : static int GDALSuggestedWarpOutput2_MustAdjustForRightBorder(
169 : GDALTransformerFunc pfnTransformer, void *pTransformArg,
170 : double* padfExtent, int nPixels, int nLines,
171 : double dfPixelSizeX, double dfPixelSizeY)
172 : {
173 : int nSamplePoints;
174 : double dfRatio;
175 : int bErr;
176 : int nBadCount;
177 634 : int abSuccess[21] = { 0 };
178 634 : double adfX[21] = { 0 };
179 634 : double adfY[21] = { 0 };
180 634 : double adfZ[21] = { 0 };
181 :
182 : //double dfMinXOut = padfExtent[0];
183 : //double dfMinYOut = padfExtent[1];
184 634 : double dfMaxXOut = padfExtent[2];
185 634 : double dfMaxYOut = padfExtent[3];
186 :
187 : // Take 20 steps
188 634 : nSamplePoints = 0;
189 13948 : for( dfRatio = 0.0; dfRatio <= 1.01; dfRatio += 0.05 )
190 : {
191 : // Ensure we end exactly at the end.
192 13314 : if( dfRatio > 0.99 )
193 634 : dfRatio = 1.0;
194 :
195 : // Along right
196 13314 : adfX[nSamplePoints] = dfMaxXOut;
197 13314 : adfY[nSamplePoints] = dfMaxYOut - dfPixelSizeY * dfRatio * nLines;
198 13314 : adfZ[nSamplePoints++] = 0.0;
199 : }
200 :
201 634 : bErr = FALSE;
202 634 : if( !pfnTransformer( pTransformArg, TRUE, nSamplePoints,
203 : adfX, adfY, adfZ, abSuccess ) )
204 : {
205 0 : bErr = TRUE;
206 : }
207 :
208 634 : if( !bErr && !pfnTransformer( pTransformArg, FALSE, nSamplePoints,
209 : adfX, adfY, adfZ, abSuccess ) )
210 : {
211 0 : bErr = TRUE;
212 : }
213 :
214 634 : nSamplePoints = 0;
215 634 : nBadCount = 0;
216 13948 : for( dfRatio = 0.0; !bErr && dfRatio <= 1.01; dfRatio += 0.05 )
217 : {
218 13314 : double expected_x = dfMaxXOut;
219 13314 : double expected_y = dfMaxYOut - dfPixelSizeY * dfRatio * nLines;
220 24964 : if (fabs(adfX[nSamplePoints] - expected_x) > dfPixelSizeX ||
221 11650 : fabs(adfY[nSamplePoints] - expected_y) > dfPixelSizeY)
222 1664 : nBadCount ++;
223 13314 : nSamplePoints ++;
224 : }
225 :
226 634 : return (nBadCount == nSamplePoints);
227 : }
228 :
229 :
230 622 : static int GDALSuggestedWarpOutput2_MustAdjustForBottomBorder(
231 : GDALTransformerFunc pfnTransformer, void *pTransformArg,
232 : double* padfExtent, int nPixels, int nLines,
233 : double dfPixelSizeX, double dfPixelSizeY)
234 : {
235 : int nSamplePoints;
236 : double dfRatio;
237 : int bErr;
238 : int nBadCount;
239 622 : int abSuccess[21] = { 0 };
240 622 : double adfX[21] = { 0 };
241 622 : double adfY[21] = { 0 };
242 622 : double adfZ[21] = { 0 };
243 :
244 622 : double dfMinXOut = padfExtent[0];
245 622 : double dfMinYOut = padfExtent[1];
246 : //double dfMaxXOut = padfExtent[2];
247 : //double dfMaxYOut = padfExtent[3];
248 :
249 : // Take 20 steps
250 622 : nSamplePoints = 0;
251 13684 : for( dfRatio = 0.0; dfRatio <= 1.01; dfRatio += 0.05 )
252 : {
253 : // Ensure we end exactly at the end.
254 13062 : if( dfRatio > 0.99 )
255 622 : dfRatio = 1.0;
256 :
257 : // Along right
258 13062 : adfX[nSamplePoints] = dfMinXOut + dfPixelSizeX * dfRatio * nPixels;
259 13062 : adfY[nSamplePoints] = dfMinYOut;
260 13062 : adfZ[nSamplePoints++] = 0.0;
261 : }
262 :
263 622 : bErr = FALSE;
264 622 : if( !pfnTransformer( pTransformArg, TRUE, nSamplePoints,
265 : adfX, adfY, adfZ, abSuccess ) )
266 : {
267 0 : bErr = TRUE;
268 : }
269 :
270 622 : if( !bErr && !pfnTransformer( pTransformArg, FALSE, nSamplePoints,
271 : adfX, adfY, adfZ, abSuccess ) )
272 : {
273 0 : bErr = TRUE;
274 : }
275 :
276 622 : nSamplePoints = 0;
277 622 : nBadCount = 0;
278 13684 : for( dfRatio = 0.0; !bErr && dfRatio <= 1.01; dfRatio += 0.05 )
279 : {
280 13062 : double expected_x = dfMinXOut + dfPixelSizeX * dfRatio * nPixels;
281 13062 : double expected_y = dfMinYOut;
282 24760 : if (fabs(adfX[nSamplePoints] - expected_x) > dfPixelSizeX ||
283 11698 : fabs(adfY[nSamplePoints] - expected_y) > dfPixelSizeY)
284 1364 : nBadCount ++;
285 13062 : nSamplePoints ++;
286 : }
287 :
288 622 : return (nBadCount == nSamplePoints);
289 : }
290 :
291 : /************************************************************************/
292 : /* GDALSuggestedWarpOutput2() */
293 : /************************************************************************/
294 :
295 : /**
296 : * Suggest output file size.
297 : *
298 : * This function is used to suggest the size, and georeferenced extents
299 : * appropriate given the indicated transformation and input file. It walks
300 : * the edges of the input file (approximately 20 sample points along each
301 : * edge) transforming into output coordinates in order to get an extents box.
302 : *
303 : * Then a resolution is computed with the intent that the length of the
304 : * distance from the top left corner of the output imagery to the bottom right
305 : * corner would represent the same number of pixels as in the source image.
306 : * Note that if the image is somewhat rotated the diagonal taken isnt of the
307 : * whole output bounding rectangle, but instead of the locations where the
308 : * top/left and bottom/right corners transform. The output pixel size is
309 : * always square. This is intended to approximately preserve the resolution
310 : * of the input data in the output file.
311 : *
312 : * The values returned in padfGeoTransformOut, pnPixels and pnLines are
313 : * the suggested number of pixels and lines for the output file, and the
314 : * geotransform relating those pixels to the output georeferenced coordinates.
315 : *
316 : * The trickiest part of using the function is ensuring that the
317 : * transformer created is from source file pixel/line coordinates to
318 : * output file georeferenced coordinates. This can be accomplished with
319 : * GDALCreateGenImgProjTransformer() by passing a NULL for the hDstDS.
320 : *
321 : * @param hSrcDS the input image (it is assumed the whole input images is
322 : * being transformed).
323 : * @param pfnTransformer the transformer function.
324 : * @param pTransformArg the callback data for the transformer function.
325 : * @param padfGeoTransformOut the array of six doubles in which the suggested
326 : * geotransform is returned.
327 : * @param pnPixels int in which the suggest pixel width of output is returned.
328 : * @param pnLines int in which the suggest pixel height of output is returned.
329 : * @param padfExtent Four entry array to return extents as (xmin, ymin, xmax, ymax).
330 : * @param nOptions Options, currently always zero.
331 : *
332 : * @return CE_None if successful or CE_Failure otherwise.
333 : */
334 :
335 : CPLErr CPL_STDCALL
336 574 : GDALSuggestedWarpOutput2( GDALDatasetH hSrcDS,
337 : GDALTransformerFunc pfnTransformer,
338 : void *pTransformArg,
339 : double *padfGeoTransformOut,
340 : int *pnPixels, int *pnLines,
341 : double *padfExtent, int nOptions )
342 :
343 : {
344 574 : VALIDATE_POINTER1( hSrcDS, "GDALSuggestedWarpOutput2", CE_Failure );
345 :
346 : /* -------------------------------------------------------------------- */
347 : /* Setup sample points all around the edge of the input raster. */
348 : /* -------------------------------------------------------------------- */
349 574 : int nSamplePoints = 0;
350 574 : int nInXSize = GDALGetRasterXSize( hSrcDS );
351 574 : int nInYSize = GDALGetRasterYSize( hSrcDS );
352 :
353 574 : if (pfnTransformer == GDALGenImgProjTransform)
354 : {
355 : /* In case CHECK_WITH_INVERT_PROJ has been modified */
356 574 : GDALRefreshGenImgProjTransformer(pTransformArg);
357 : }
358 :
359 : #define N_PIXELSTEP 50
360 574 : int nSteps = (int) (double(MIN(nInYSize, nInXSize)) / N_PIXELSTEP + .5);
361 574 : if (nSteps < 20)
362 562 : nSteps = 20;
363 :
364 : retry:
365 574 : int nSampleMax = (nSteps + 1)*(nSteps + 1);
366 574 : int *pabSuccess = NULL;
367 : double *padfX, *padfY, *padfZ;
368 : double *padfXRevert, *padfYRevert, *padfZRevert;
369 :
370 574 : double dfRatio = 0.0;
371 574 : double dfStep = 1. / nSteps;
372 :
373 574 : pabSuccess = (int *) VSIMalloc3(sizeof(int), nSteps + 1, nSteps + 1);
374 574 : padfX = (double *) VSIMalloc3(sizeof(double) * 3, nSteps + 1, nSteps + 1);
375 574 : padfXRevert = (double *) VSIMalloc3(sizeof(double) * 3, nSteps + 1, nSteps + 1);
376 574 : if (pabSuccess == NULL || padfX == NULL || padfXRevert == NULL)
377 : {
378 0 : CPLFree( padfX );
379 0 : CPLFree( padfXRevert );
380 0 : CPLFree( pabSuccess );
381 0 : if (nSteps > 20)
382 : {
383 0 : nSteps = 20;
384 0 : goto retry;
385 : }
386 0 : return CE_Failure;
387 : }
388 574 : padfY = padfX + nSampleMax;
389 574 : padfZ = padfX + nSampleMax * 2;
390 574 : padfYRevert = padfXRevert + nSampleMax;
391 574 : padfZRevert = padfXRevert + nSampleMax * 2;
392 :
393 :
394 : // Take N_STEPS steps
395 : int iStep;
396 12628 : for( iStep = 0; iStep <= nSteps; iStep ++ )
397 : {
398 12054 : dfRatio = (iStep == nSteps) ? 1.0 : iStep * dfStep;
399 :
400 : // Along top
401 12054 : padfX[nSamplePoints] = dfRatio * nInXSize;
402 12054 : padfY[nSamplePoints] = 0.0;
403 12054 : padfZ[nSamplePoints++] = 0.0;
404 :
405 : // Along bottom
406 12054 : padfX[nSamplePoints] = dfRatio * nInXSize;
407 12054 : padfY[nSamplePoints] = nInYSize;
408 12054 : padfZ[nSamplePoints++] = 0.0;
409 :
410 : // Along left
411 12054 : padfX[nSamplePoints] = 0.0;
412 12054 : padfY[nSamplePoints] = dfRatio * nInYSize;
413 12054 : padfZ[nSamplePoints++] = 0.0;
414 :
415 : // Along right
416 12054 : padfX[nSamplePoints] = nInXSize;
417 12054 : padfY[nSamplePoints] = dfRatio * nInYSize;
418 12054 : padfZ[nSamplePoints++] = 0.0;
419 : }
420 :
421 574 : CPLAssert( nSamplePoints == 4 * (nSteps + 1) );
422 :
423 574 : memset( pabSuccess, 1, sizeof(int) * nSampleMax );
424 :
425 : /* -------------------------------------------------------------------- */
426 : /* Transform them to the output coordinate system. */
427 : /* -------------------------------------------------------------------- */
428 574 : int nFailedCount = 0, i;
429 :
430 574 : if( !pfnTransformer( pTransformArg, FALSE, nSamplePoints,
431 : padfX, padfY, padfZ, pabSuccess ) )
432 : {
433 : CPLError( CE_Failure, CPLE_AppDefined,
434 : "GDALSuggestedWarpOutput() failed because the passed\n"
435 0 : "transformer failed." );
436 0 : CPLFree( padfX );
437 0 : CPLFree( padfXRevert );
438 0 : CPLFree( pabSuccess );
439 0 : return CE_Failure;
440 : }
441 :
442 48790 : for( i = 0; i < nSamplePoints; i++ )
443 : {
444 48216 : if( !pabSuccess[i] )
445 598 : nFailedCount++;
446 : }
447 :
448 : /* -------------------------------------------------------------------- */
449 : /* Check if the computed target coordinates are revertable. */
450 : /* If not, try the detailed grid sampling. */
451 : /* -------------------------------------------------------------------- */
452 574 : if (nFailedCount == 0 )
453 : {
454 562 : memcpy(padfXRevert, padfX, nSamplePoints * sizeof(double));
455 562 : memcpy(padfYRevert, padfY, nSamplePoints * sizeof(double));
456 562 : memcpy(padfZRevert, padfZ, nSamplePoints * sizeof(double));
457 562 : if( !pfnTransformer( pTransformArg, TRUE, nSamplePoints,
458 : padfXRevert, padfYRevert, padfZRevert, pabSuccess ) )
459 : {
460 0 : nFailedCount = 1;
461 : }
462 : else
463 : {
464 11758 : for( i = 0; nFailedCount == 0 && i < nSamplePoints; i++ )
465 : {
466 11196 : if( !pabSuccess[i] )
467 0 : nFailedCount++;
468 :
469 11196 : dfRatio = 0.0 + (i/4) * dfStep;
470 11196 : if (dfRatio>0.99)
471 424 : dfRatio = 1.0;
472 :
473 : double dfExpectedX, dfExpectedY;
474 11196 : if ((i % 4) == 0)
475 : {
476 3140 : dfExpectedX = dfRatio * nInXSize;
477 3140 : dfExpectedY = 0.0;
478 : }
479 8056 : else if ((i % 4) == 1)
480 : {
481 2686 : dfExpectedX = dfRatio * nInXSize;
482 2686 : dfExpectedY = nInYSize;
483 : }
484 5370 : else if ((i % 4) == 2)
485 : {
486 2686 : dfExpectedX = 0.0;
487 2686 : dfExpectedY = dfRatio * nInYSize;
488 : }
489 : else
490 : {
491 2684 : dfExpectedX = nInXSize;
492 2684 : dfExpectedY = dfRatio * nInYSize;
493 : }
494 :
495 21938 : if (fabs(padfXRevert[i] - dfExpectedX) > nInXSize / nSteps ||
496 10742 : fabs(padfYRevert[i] - dfExpectedY) > nInYSize / nSteps)
497 456 : nFailedCount ++;
498 : }
499 : }
500 : }
501 :
502 : /* -------------------------------------------------------------------- */
503 : /* If any of the edge points failed to transform, we need to */
504 : /* build a fairly detailed internal grid of points instead to */
505 : /* help identify the area that is transformable. */
506 : /* -------------------------------------------------------------------- */
507 574 : if( nFailedCount > 0 )
508 : {
509 : int iStep2;
510 : double dfRatio2;
511 468 : nSamplePoints = 0;
512 :
513 : // Take N_STEPS steps
514 10296 : for( iStep = 0; iStep <= nSteps; iStep ++ )
515 : {
516 9828 : dfRatio = (iStep == nSteps) ? 1.0 : iStep * dfStep;
517 :
518 216216 : for( iStep2 = 0; iStep2 <= nSteps; iStep2 ++ )
519 : {
520 206388 : dfRatio2 = (iStep2 == nSteps) ? 1.0 : iStep2 * dfStep;
521 :
522 : // Along top
523 206388 : padfX[nSamplePoints] = dfRatio2 * nInXSize;
524 206388 : padfY[nSamplePoints] = dfRatio * nInYSize;
525 206388 : padfZ[nSamplePoints++] = 0.0;
526 : }
527 : }
528 :
529 468 : CPLAssert( nSamplePoints == nSampleMax );
530 :
531 468 : if( !pfnTransformer( pTransformArg, FALSE, nSamplePoints,
532 : padfX, padfY, padfZ, pabSuccess ) )
533 : {
534 : CPLError( CE_Failure, CPLE_AppDefined,
535 : "GDALSuggestedWarpOutput() failed because the passed\n"
536 0 : "transformer failed." );
537 :
538 0 : CPLFree( padfX );
539 0 : CPLFree( padfXRevert );
540 0 : CPLFree( pabSuccess );
541 :
542 0 : return CE_Failure;
543 : }
544 : }
545 :
546 : /* -------------------------------------------------------------------- */
547 : /* Collect the bounds, ignoring any failed points. */
548 : /* -------------------------------------------------------------------- */
549 574 : double dfMinXOut=0, dfMinYOut=0, dfMaxXOut=0, dfMaxYOut=0;
550 574 : int bGotInitialPoint = FALSE;
551 :
552 574 : nFailedCount = 0;
553 215866 : for( i = 0; i < nSamplePoints; i++ )
554 : {
555 :
556 215292 : int x_i = i % (nSteps + 1);
557 215292 : int y_i = i / (nSteps + 1);
558 :
559 215292 : if (x_i > 0 && (pabSuccess[i-1] || pabSuccess[i]))
560 : {
561 204116 : double x_out_before = padfX[i-1];
562 204116 : double x_out_after = padfX[i];
563 204116 : int nIter = 0;
564 204116 : double x_in_before = (x_i - 1) * nInXSize * 1.0 / nSteps;
565 204116 : double x_in_after = x_i * nInXSize * 1.0 / nSteps;
566 204116 : int valid_before = pabSuccess[i-1];
567 204116 : int valid_after = pabSuccess[i];
568 :
569 : /* Detect discontinuity in target coordinates when the target x coordinates */
570 : /* change sign. This may be a false positive when the targe tx is around 0 */
571 : /* Dichotomic search to reduce the interval to near the discontinuity and */
572 : /* get a better out extent */
573 414440 : while ( (!valid_before || !valid_after ||
574 : x_out_before * x_out_after < 0) && nIter < 16 )
575 : {
576 6208 : double x = (x_in_before + x_in_after) / 2;
577 6208 : double y = y_i * nInYSize * 1.0 / nSteps;
578 6208 : double z= 0;
579 : //fprintf(stderr, "[%d] (%f, %f) -> ", nIter, x, y);
580 6208 : int bSuccess = TRUE;
581 6208 : if( !pfnTransformer( pTransformArg, FALSE, 1,
582 : &x, &y, &z, &bSuccess ) || !bSuccess )
583 : {
584 : //fprintf(stderr, "invalid\n");
585 486 : if (!valid_before)
586 : {
587 342 : x_in_before = (x_in_before + x_in_after) / 2;
588 : }
589 144 : else if (!valid_after)
590 : {
591 144 : x_in_after = (x_in_before + x_in_after) / 2;
592 : }
593 : else
594 0 : break;
595 : }
596 : else
597 : {
598 : //fprintf(stderr, "(%f, %f)\n", x, y);
599 :
600 5722 : if( !bGotInitialPoint )
601 : {
602 6 : bGotInitialPoint = TRUE;
603 6 : dfMinXOut = dfMaxXOut = x;
604 6 : dfMinYOut = dfMaxYOut = y;
605 : }
606 : else
607 : {
608 5716 : dfMinXOut = MIN(dfMinXOut,x);
609 5716 : dfMinYOut = MIN(dfMinYOut,y);
610 5716 : dfMaxXOut = MAX(dfMaxXOut,x);
611 5716 : dfMaxYOut = MAX(dfMaxYOut,y);
612 : }
613 :
614 8410 : if (!valid_before || x_out_before * x < 0)
615 : {
616 2688 : valid_after = TRUE;
617 2688 : x_in_after = (x_in_before + x_in_after) / 2;
618 2688 : x_out_after = x;
619 : }
620 : else
621 : {
622 3034 : valid_before = TRUE;
623 3034 : x_out_before = x;
624 3034 : x_in_before = (x_in_before + x_in_after) / 2;
625 : }
626 : }
627 6208 : nIter ++;
628 : }
629 : }
630 :
631 215292 : if( !pabSuccess[i] )
632 : {
633 1006 : nFailedCount++;
634 1006 : continue;
635 : }
636 :
637 214286 : if( !bGotInitialPoint )
638 : {
639 568 : bGotInitialPoint = TRUE;
640 568 : dfMinXOut = dfMaxXOut = padfX[i];
641 568 : dfMinYOut = dfMaxYOut = padfY[i];
642 : }
643 : else
644 : {
645 213718 : dfMinXOut = MIN(dfMinXOut, padfX[i]);
646 213718 : dfMinYOut = MIN(dfMinYOut, padfY[i]);
647 213718 : dfMaxXOut = MAX(dfMaxXOut, padfX[i]);
648 213718 : dfMaxYOut = MAX(dfMaxYOut, padfY[i]);
649 : }
650 : }
651 :
652 574 : if( nFailedCount > nSamplePoints - 10 )
653 : {
654 : CPLError( CE_Failure, CPLE_AppDefined,
655 : "Too many points (%d out of %d) failed to transform,\n"
656 : "unable to compute output bounds.",
657 0 : nFailedCount, nSamplePoints );
658 :
659 0 : CPLFree( padfX );
660 0 : CPLFree( padfXRevert );
661 0 : CPLFree( pabSuccess );
662 :
663 0 : return CE_Failure;
664 : }
665 :
666 574 : if( nFailedCount > 0 )
667 : CPLDebug( "GDAL",
668 : "GDALSuggestedWarpOutput(): %d out of %d points failed to transform.",
669 12 : nFailedCount, nSamplePoints );
670 :
671 : /* -------------------------------------------------------------------- */
672 : /* Compute the distance in "georeferenced" units from the top */
673 : /* corner of the transformed input image to the bottom left */
674 : /* corner of the transformed input. Use this distance to */
675 : /* compute an approximate pixel size in the output */
676 : /* georeferenced coordinates. */
677 : /* -------------------------------------------------------------------- */
678 : double dfDiagonalDist, dfDeltaX, dfDeltaY;
679 :
680 1136 : if( pabSuccess[0] && pabSuccess[nSamplePoints - 1] )
681 : {
682 562 : dfDeltaX = padfX[nSamplePoints-1] - padfX[0];
683 562 : dfDeltaY = padfY[nSamplePoints-1] - padfY[0];
684 : }
685 : else
686 : {
687 12 : dfDeltaX = dfMaxXOut - dfMinXOut;
688 12 : dfDeltaY = dfMaxYOut - dfMinYOut;
689 : }
690 :
691 574 : dfDiagonalDist = sqrt( dfDeltaX * dfDeltaX + dfDeltaY * dfDeltaY );
692 :
693 : /* -------------------------------------------------------------------- */
694 : /* Compute a pixel size from this. */
695 : /* -------------------------------------------------------------------- */
696 : double dfPixelSize;
697 :
698 : dfPixelSize = dfDiagonalDist
699 574 : / sqrt(((double)nInXSize)*nInXSize + ((double)nInYSize)*nInYSize);
700 :
701 574 : *pnPixels = (int) ((dfMaxXOut - dfMinXOut) / dfPixelSize + 0.5);
702 574 : *pnLines = (int) ((dfMaxYOut - dfMinYOut) / dfPixelSize + 0.5);
703 :
704 574 : double dfPixelSizeX = dfPixelSize;
705 574 : double dfPixelSizeY = dfPixelSize;
706 :
707 : double adfExtent[4];
708 574 : const double adfRatioArray[] = { 0, 0.001, 0.01, 0.1, 1 };
709 : size_t nRetry;
710 :
711 : #define N_ELEMENTS(x) (sizeof(x) / sizeof(x[0]))
712 :
713 : /* -------------------------------------------------------------------- */
714 : /* Check that the right border is not completely out of source */
715 : /* image. If so, adjust the x pixel size a bit in the hope it will */
716 : /* fit. */
717 : /* -------------------------------------------------------------------- */
718 638 : for( nRetry = 0; nRetry < N_ELEMENTS(adfRatioArray); nRetry ++ )
719 : {
720 : double dfTryPixelSizeX =
721 634 : dfPixelSizeX - dfPixelSizeX * adfRatioArray[nRetry] / *pnPixels;
722 634 : adfExtent[0] = dfMinXOut;
723 634 : adfExtent[1] = dfMaxYOut - (*pnLines) * dfPixelSizeY;
724 634 : adfExtent[2] = dfMinXOut + (*pnPixels) * dfTryPixelSizeX;
725 634 : adfExtent[3] = dfMaxYOut;
726 634 : if (!GDALSuggestedWarpOutput2_MustAdjustForRightBorder(
727 : pfnTransformer, pTransformArg,
728 : adfExtent, *pnPixels, *pnLines,
729 : dfTryPixelSizeX, dfPixelSizeY))
730 : {
731 570 : dfPixelSizeX = dfTryPixelSizeX;
732 570 : break;
733 : }
734 : }
735 :
736 : /* -------------------------------------------------------------------- */
737 : /* Check that the bottom border is not completely out of source */
738 : /* image. If so, adjust the y pixel size a bit in the hope it will */
739 : /* fit. */
740 : /* -------------------------------------------------------------------- */
741 626 : for( nRetry = 0; nRetry < N_ELEMENTS(adfRatioArray); nRetry ++ )
742 : {
743 : double dfTryPixelSizeY =
744 622 : dfPixelSizeY - dfPixelSizeY * adfRatioArray[nRetry] / *pnLines;
745 622 : adfExtent[0] = dfMinXOut;
746 622 : adfExtent[1] = dfMaxYOut - (*pnLines) * dfTryPixelSizeY;
747 622 : adfExtent[2] = dfMinXOut + (*pnPixels) * dfPixelSizeX;
748 622 : adfExtent[3] = dfMaxYOut;
749 622 : if (!GDALSuggestedWarpOutput2_MustAdjustForBottomBorder(
750 : pfnTransformer, pTransformArg,
751 : adfExtent, *pnPixels, *pnLines,
752 : dfPixelSizeX, dfTryPixelSizeY))
753 : {
754 570 : dfPixelSizeY = dfTryPixelSizeY;
755 570 : break;
756 : }
757 : }
758 :
759 :
760 : /* -------------------------------------------------------------------- */
761 : /* Recompute some bounds so that all return values are consistant */
762 : /* -------------------------------------------------------------------- */
763 574 : dfMaxXOut = dfMinXOut + (*pnPixels) * dfPixelSizeX;
764 574 : dfMinYOut = dfMaxYOut - (*pnLines) * dfPixelSizeY;
765 :
766 : /* -------------------------------------------------------------------- */
767 : /* Return raw extents. */
768 : /* -------------------------------------------------------------------- */
769 574 : padfExtent[0] = dfMinXOut;
770 574 : padfExtent[1] = dfMinYOut;
771 574 : padfExtent[2] = dfMaxXOut;
772 574 : padfExtent[3] = dfMaxYOut;
773 :
774 : /* -------------------------------------------------------------------- */
775 : /* Set the output geotransform. */
776 : /* -------------------------------------------------------------------- */
777 574 : padfGeoTransformOut[0] = dfMinXOut;
778 574 : padfGeoTransformOut[1] = dfPixelSizeX;
779 574 : padfGeoTransformOut[2] = 0.0;
780 574 : padfGeoTransformOut[3] = dfMaxYOut;
781 574 : padfGeoTransformOut[4] = 0.0;
782 574 : padfGeoTransformOut[5] = - dfPixelSizeY;
783 :
784 574 : CPLFree( padfX );
785 574 : CPLFree( padfXRevert );
786 574 : CPLFree( pabSuccess );
787 :
788 574 : return CE_None;
789 : }
790 :
791 : /************************************************************************/
792 : /* ==================================================================== */
793 : /* GDALGenImgProjTransformer */
794 : /* ==================================================================== */
795 : /************************************************************************/
796 :
797 : typedef struct {
798 :
799 : GDALTransformerInfo sTI;
800 :
801 : double adfSrcGeoTransform[6];
802 : double adfSrcInvGeoTransform[6];
803 :
804 : void *pSrcGCPTransformArg;
805 : void *pSrcRPCTransformArg;
806 : void *pSrcTPSTransformArg;
807 : void *pSrcGeoLocTransformArg;
808 :
809 : void *pReprojectArg;
810 :
811 : double adfDstGeoTransform[6];
812 : double adfDstInvGeoTransform[6];
813 :
814 : void *pDstGCPTransformArg;
815 :
816 : } GDALGenImgProjTransformInfo;
817 :
818 : /************************************************************************/
819 : /* GDALCreateGenImgProjTransformer() */
820 : /************************************************************************/
821 :
822 : /**
823 : * Create image to image transformer.
824 : *
825 : * This function creates a transformation object that maps from pixel/line
826 : * coordinates on one image to pixel/line coordinates on another image. The
827 : * images may potentially be georeferenced in different coordinate systems,
828 : * and may used GCPs to map between their pixel/line coordinates and
829 : * georeferenced coordinates (as opposed to the default assumption that their
830 : * geotransform should be used).
831 : *
832 : * This transformer potentially performs three concatenated transformations.
833 : *
834 : * The first stage is from source image pixel/line coordinates to source
835 : * image georeferenced coordinates, and may be done using the geotransform,
836 : * or if not defined using a polynomial model derived from GCPs. If GCPs
837 : * are used this stage is accomplished using GDALGCPTransform().
838 : *
839 : * The second stage is to change projections from the source coordinate system
840 : * to the destination coordinate system, assuming they differ. This is
841 : * accomplished internally using GDALReprojectionTransform().
842 : *
843 : * The third stage is converting from destination image georeferenced
844 : * coordinates to destination image coordinates. This is done using the
845 : * destination image geotransform, or if not available, using a polynomial
846 : * model derived from GCPs. If GCPs are used this stage is accomplished using
847 : * GDALGCPTransform(). This stage is skipped if hDstDS is NULL when the
848 : * transformation is created.
849 : *
850 : * @param hSrcDS source dataset, or NULL.
851 : * @param pszSrcWKT the coordinate system for the source dataset. If NULL,
852 : * it will be read from the dataset itself.
853 : * @param hDstDS destination dataset (or NULL).
854 : * @param pszDstWKT the coordinate system for the destination dataset. If
855 : * NULL, and hDstDS not NULL, it will be read from the destination dataset.
856 : * @param bGCPUseOK TRUE if GCPs should be used if the geotransform is not
857 : * available on the source dataset (not destination).
858 : * @param dfGCPErrorThreshold ignored/deprecated.
859 : * @param nOrder the maximum order to use for GCP derived polynomials if
860 : * possible. Use 0 to autoselect, or -1 for thin plate splines.
861 : *
862 : * @return handle suitable for use GDALGenImgProjTransform(), and to be
863 : * deallocated with GDALDestroyGenImgProjTransformer().
864 : */
865 :
866 : void *
867 94 : GDALCreateGenImgProjTransformer( GDALDatasetH hSrcDS, const char *pszSrcWKT,
868 : GDALDatasetH hDstDS, const char *pszDstWKT,
869 : int bGCPUseOK, double dfGCPErrorThreshold,
870 : int nOrder )
871 :
872 : {
873 94 : char **papszOptions = NULL;
874 : void *pRet;
875 :
876 94 : if( pszSrcWKT != NULL )
877 66 : papszOptions = CSLSetNameValue( papszOptions, "SRC_SRS", pszSrcWKT );
878 94 : if( pszDstWKT != NULL )
879 62 : papszOptions = CSLSetNameValue( papszOptions, "DST_SRS", pszDstWKT );
880 94 : if( !bGCPUseOK )
881 18 : papszOptions = CSLSetNameValue( papszOptions, "GCPS_OK", "FALSE" );
882 94 : if( nOrder != 0 )
883 : papszOptions = CSLSetNameValue( papszOptions, "MAX_GCP_ORDER",
884 0 : CPLString().Printf("%d",nOrder) );
885 :
886 94 : pRet = GDALCreateGenImgProjTransformer2( hSrcDS, hDstDS, papszOptions );
887 94 : CSLDestroy( papszOptions );
888 :
889 94 : return pRet;
890 : }
891 :
892 :
893 :
894 : /************************************************************************/
895 : /* InsertCenterLong() */
896 : /* */
897 : /* Insert a CENTER_LONG Extension entry on a GEOGCS to indicate */
898 : /* the center longitude of the dataset for wrapping purposes. */
899 : /************************************************************************/
900 :
901 90 : static CPLString InsertCenterLong( GDALDatasetH hDS, CPLString osWKT )
902 :
903 : {
904 90 : if( !EQUALN(osWKT.c_str(), "GEOGCS[", 7) )
905 42 : return osWKT;
906 :
907 48 : if( strstr(osWKT,"EXTENSION[\"CENTER_LONG") != NULL )
908 0 : return osWKT;
909 :
910 : /* -------------------------------------------------------------------- */
911 : /* For now we only do this if we have a geotransform since */
912 : /* other forms require a bunch of extra work. */
913 : /* -------------------------------------------------------------------- */
914 : double adfGeoTransform[6];
915 :
916 48 : if( GDALGetGeoTransform( hDS, adfGeoTransform ) != CE_None )
917 0 : return osWKT;
918 :
919 : /* -------------------------------------------------------------------- */
920 : /* Compute min/max longitude based on testing the four corners. */
921 : /* -------------------------------------------------------------------- */
922 : double dfMinLong, dfMaxLong;
923 48 : int nXSize = GDALGetRasterXSize( hDS );
924 48 : int nYSize = GDALGetRasterYSize( hDS );
925 :
926 : dfMinLong =
927 1008 : MIN(MIN(adfGeoTransform[0] + 0 * adfGeoTransform[1]
928 : + 0 * adfGeoTransform[2],
929 : adfGeoTransform[0] + nXSize * adfGeoTransform[1]
930 : + 0 * adfGeoTransform[2]),
931 : MIN(adfGeoTransform[0] + 0 * adfGeoTransform[1]
932 : + nYSize * adfGeoTransform[2],
933 : adfGeoTransform[0] + nXSize * adfGeoTransform[1]
934 1008 : + nYSize * adfGeoTransform[2]));
935 : dfMaxLong =
936 1152 : MAX(MAX(adfGeoTransform[0] + 0 * adfGeoTransform[1]
937 : + 0 * adfGeoTransform[2],
938 : adfGeoTransform[0] + nXSize * adfGeoTransform[1]
939 : + 0 * adfGeoTransform[2]),
940 : MAX(adfGeoTransform[0] + 0 * adfGeoTransform[1]
941 : + nYSize * adfGeoTransform[2],
942 : adfGeoTransform[0] + nXSize * adfGeoTransform[1]
943 1152 : + nYSize * adfGeoTransform[2]));
944 :
945 48 : if( dfMaxLong - dfMinLong > 360.0 )
946 0 : return osWKT;
947 :
948 : /* -------------------------------------------------------------------- */
949 : /* Insert center long. */
950 : /* -------------------------------------------------------------------- */
951 48 : OGRSpatialReference oSRS( osWKT );
952 48 : double dfCenterLong = (dfMaxLong + dfMinLong) / 2.0;
953 : OGR_SRSNode *poExt;
954 :
955 48 : poExt = new OGR_SRSNode( "EXTENSION" );
956 96 : poExt->AddChild( new OGR_SRSNode( "CENTER_LONG" ) );
957 48 : poExt->AddChild( new OGR_SRSNode( CPLString().Printf("%g",dfCenterLong) ));
958 :
959 48 : oSRS.GetRoot()->AddChild( poExt->Clone() );
960 48 : delete poExt;
961 :
962 : /* -------------------------------------------------------------------- */
963 : /* Convert back to wkt. */
964 : /* -------------------------------------------------------------------- */
965 48 : char *pszWKT = NULL;
966 48 : oSRS.exportToWkt( &pszWKT );
967 :
968 48 : osWKT = pszWKT;
969 48 : CPLFree( pszWKT );
970 :
971 48 : return osWKT;
972 : }
973 :
974 : /************************************************************************/
975 : /* GDALCreateGenImgProjTransformer2() */
976 : /************************************************************************/
977 :
978 : /**
979 : * Create image to image transformer.
980 : *
981 : * This function creates a transformation object that maps from pixel/line
982 : * coordinates on one image to pixel/line coordinates on another image. The
983 : * images may potentially be georeferenced in different coordinate systems,
984 : * and may used GCPs to map between their pixel/line coordinates and
985 : * georeferenced coordinates (as opposed to the default assumption that their
986 : * geotransform should be used).
987 : *
988 : * This transformer potentially performs three concatenated transformations.
989 : *
990 : * The first stage is from source image pixel/line coordinates to source
991 : * image georeferenced coordinates, and may be done using the geotransform,
992 : * or if not defined using a polynomial model derived from GCPs. If GCPs
993 : * are used this stage is accomplished using GDALGCPTransform().
994 : *
995 : * The second stage is to change projections from the source coordinate system
996 : * to the destination coordinate system, assuming they differ. This is
997 : * accomplished internally using GDALReprojectionTransform().
998 : *
999 : * The third stage is converting from destination image georeferenced
1000 : * coordinates to destination image coordinates. This is done using the
1001 : * destination image geotransform, or if not available, using a polynomial
1002 : * model derived from GCPs. If GCPs are used this stage is accomplished using
1003 : * GDALGCPTransform(). This stage is skipped if hDstDS is NULL when the
1004 : * transformation is created.
1005 : *
1006 : * Supported Options:
1007 : * <ul>
1008 : * <li> SRC_SRS: WKT SRS to be used as an override for hSrcDS.
1009 : * <li> DST_SRS: WKT SRS to be used as an override for hDstDS.
1010 : * <li> GCPS_OK: If false, GCPs will not be used, default is TRUE.
1011 : * <li> REFINE_MINIMUM_GCPS: The minimum amount of GCPs that should be available after the refinement.
1012 : * <li> REFINE_TOLERANCE: The tolernace that specifies when a GCP will be eliminated.
1013 : * <li> MAX_GCP_ORDER: the maximum order to use for GCP derived polynomials if
1014 : * possible. The default is to autoselect based on the number of GCPs.
1015 : * A value of -1 triggers use of Thin Plate Spline instead of polynomials.
1016 : * <li> METHOD: may have a value which is one of GEOTRANSFORM, GCP_POLYNOMIAL,
1017 : * GCP_TPS, GEOLOC_ARRAY, RPC to force only one geolocation method to be
1018 : * considered on the source dataset.
1019 : * <li> RPC_HEIGHT: A fixed height to be used with RPC calculations.
1020 : * <li> RPC_DEM: The name of a DEM file to be used with RPC calculations.
1021 : * <li> INSERT_CENTER_LONG: May be set to FALSE to disable setting up a
1022 : * CENTER_LONG value on the coordinate system to rewrap things around the
1023 : * center of the image.
1024 : * </ul>
1025 : *
1026 : * @param hSrcDS source dataset, or NULL.
1027 : * @param hDstDS destination dataset (or NULL).
1028 : * @param papszOptions NULL-terminated list of string options (or NULL).
1029 : *
1030 : * @return handle suitable for use GDALGenImgProjTransform(), and to be
1031 : * deallocated with GDALDestroyGenImgProjTransformer() or NULL on failure.
1032 : */
1033 :
1034 : void *
1035 670 : GDALCreateGenImgProjTransformer2( GDALDatasetH hSrcDS, GDALDatasetH hDstDS,
1036 : char **papszOptions )
1037 :
1038 : {
1039 : GDALGenImgProjTransformInfo *psInfo;
1040 : char **papszMD;
1041 : GDALRPCInfo sRPCInfo;
1042 670 : const char *pszMethod = CSLFetchNameValue( papszOptions, "METHOD" );
1043 : const char *pszValue;
1044 670 : int nOrder = 0, bGCPUseOK = TRUE, nMinimumGcps = -1, bRefine = FALSE;
1045 670 : double dfTolerance = 0.0;
1046 670 : const char *pszSrcWKT = CSLFetchNameValue( papszOptions, "SRC_SRS" );
1047 670 : const char *pszDstWKT = CSLFetchNameValue( papszOptions, "DST_SRS" );
1048 :
1049 670 : pszValue = CSLFetchNameValue( papszOptions, "MAX_GCP_ORDER" );
1050 670 : if( pszValue )
1051 0 : nOrder = atoi(pszValue);
1052 :
1053 670 : pszValue = CSLFetchNameValue( papszOptions, "GCPS_OK" );
1054 670 : if( pszValue )
1055 18 : bGCPUseOK = CSLTestBoolean(pszValue);
1056 :
1057 670 : pszValue = CSLFetchNameValue( papszOptions, "REFINE_MINIMUM_GCPS" );
1058 670 : if( pszValue )
1059 : {
1060 0 : if( atoi(pszValue) != -1)
1061 0 : nMinimumGcps = atoi(pszValue);
1062 : }
1063 :
1064 670 : pszValue = CSLFetchNameValue( papszOptions, "REFINE_TOLERANCE" );
1065 670 : if( pszValue )
1066 : {
1067 0 : dfTolerance = atof(pszValue);
1068 0 : bRefine = TRUE;
1069 : }
1070 :
1071 : /* -------------------------------------------------------------------- */
1072 : /* Initialize the transform info. */
1073 : /* -------------------------------------------------------------------- */
1074 : psInfo = (GDALGenImgProjTransformInfo *)
1075 670 : CPLCalloc(sizeof(GDALGenImgProjTransformInfo),1);
1076 :
1077 670 : strcpy( psInfo->sTI.szSignature, "GTI" );
1078 670 : psInfo->sTI.pszClassName = "GDALGenImgProjTransformer";
1079 670 : psInfo->sTI.pfnTransform = GDALGenImgProjTransform;
1080 670 : psInfo->sTI.pfnCleanup = GDALDestroyGenImgProjTransformer;
1081 670 : psInfo->sTI.pfnSerialize = GDALSerializeGenImgProjTransformer;
1082 :
1083 : /* -------------------------------------------------------------------- */
1084 : /* Get forward and inverse geotransform for the source image. */
1085 : /* -------------------------------------------------------------------- */
1086 670 : if( hSrcDS == NULL )
1087 : {
1088 22 : psInfo->adfSrcGeoTransform[0] = 0.0;
1089 22 : psInfo->adfSrcGeoTransform[1] = 1.0;
1090 22 : psInfo->adfSrcGeoTransform[2] = 0.0;
1091 22 : psInfo->adfSrcGeoTransform[3] = 0.0;
1092 22 : psInfo->adfSrcGeoTransform[4] = 0.0;
1093 22 : psInfo->adfSrcGeoTransform[5] = 1.0;
1094 : memcpy( psInfo->adfSrcInvGeoTransform, psInfo->adfSrcGeoTransform,
1095 22 : sizeof(double) * 6 );
1096 : }
1097 :
1098 1250 : else if( (pszMethod == NULL || EQUAL(pszMethod,"GEOTRANSFORM"))
1099 : && GDALGetGeoTransform( hSrcDS, psInfo->adfSrcGeoTransform )
1100 : == CE_None
1101 602 : && (psInfo->adfSrcGeoTransform[0] != 0.0
1102 0 : || psInfo->adfSrcGeoTransform[1] != 1.0
1103 0 : || psInfo->adfSrcGeoTransform[2] != 0.0
1104 0 : || psInfo->adfSrcGeoTransform[3] != 0.0
1105 0 : || psInfo->adfSrcGeoTransform[4] != 0.0
1106 0 : || ABS(psInfo->adfSrcGeoTransform[5]) != 1.0) )
1107 : {
1108 : GDALInvGeoTransform( psInfo->adfSrcGeoTransform,
1109 602 : psInfo->adfSrcInvGeoTransform );
1110 602 : if( pszSrcWKT == NULL )
1111 542 : pszSrcWKT = GDALGetProjectionRef( hSrcDS );
1112 : }
1113 :
1114 46 : else if( bGCPUseOK
1115 : && (pszMethod == NULL || EQUAL(pszMethod,"GCP_POLYNOMIAL") )
1116 : && GDALGetGCPCount( hSrcDS ) > 0 && nOrder >= 0 )
1117 : {
1118 32 : if(bRefine)
1119 : {
1120 : psInfo->pSrcGCPTransformArg =
1121 : GDALCreateGCPRefineTransformer( GDALGetGCPCount( hSrcDS ),
1122 : GDALGetGCPs( hSrcDS ), nOrder,
1123 0 : FALSE, dfTolerance, nMinimumGcps );
1124 : }
1125 : else
1126 : {
1127 : psInfo->pSrcGCPTransformArg =
1128 : GDALCreateGCPTransformer( GDALGetGCPCount( hSrcDS ),
1129 : GDALGetGCPs( hSrcDS ), nOrder,
1130 32 : FALSE );
1131 : }
1132 :
1133 32 : if( psInfo->pSrcGCPTransformArg == NULL )
1134 : {
1135 0 : GDALDestroyGenImgProjTransformer( psInfo );
1136 0 : return NULL;
1137 : }
1138 :
1139 32 : if( pszSrcWKT == NULL )
1140 32 : pszSrcWKT = GDALGetGCPProjection( hSrcDS );
1141 : }
1142 :
1143 14 : else if( bGCPUseOK
1144 : && GDALGetGCPCount( hSrcDS ) > 0
1145 : && nOrder <= 0
1146 : && (pszMethod == NULL || EQUAL(pszMethod,"GCP_TPS")) )
1147 : {
1148 : psInfo->pSrcTPSTransformArg =
1149 : GDALCreateTPSTransformer( GDALGetGCPCount( hSrcDS ),
1150 6 : GDALGetGCPs( hSrcDS ), FALSE );
1151 6 : if( psInfo->pSrcTPSTransformArg == NULL )
1152 : {
1153 0 : GDALDestroyGenImgProjTransformer( psInfo );
1154 0 : return NULL;
1155 : }
1156 :
1157 6 : if( pszSrcWKT == NULL )
1158 6 : pszSrcWKT = GDALGetGCPProjection( hSrcDS );
1159 : }
1160 :
1161 8 : else if( (pszMethod == NULL || EQUAL(pszMethod,"RPC"))
1162 : && (papszMD = GDALGetMetadata( hSrcDS, "RPC" )) != NULL
1163 : && GDALExtractRPCInfo( papszMD, &sRPCInfo ) )
1164 : {
1165 : psInfo->pSrcRPCTransformArg =
1166 6 : GDALCreateRPCTransformer( &sRPCInfo, FALSE, 0.1, papszOptions );
1167 6 : if( psInfo->pSrcRPCTransformArg == NULL )
1168 : {
1169 0 : GDALDestroyGenImgProjTransformer( psInfo );
1170 0 : return NULL;
1171 : }
1172 6 : if( pszSrcWKT == NULL )
1173 6 : pszSrcWKT = SRS_WKT_WGS84;
1174 : }
1175 :
1176 2 : else if( (pszMethod == NULL || EQUAL(pszMethod,"GEOLOC_ARRAY"))
1177 : && (papszMD = GDALGetMetadata( hSrcDS, "GEOLOCATION" )) != NULL )
1178 : {
1179 : psInfo->pSrcGeoLocTransformArg =
1180 2 : GDALCreateGeoLocTransformer( hSrcDS, papszMD, FALSE );
1181 :
1182 2 : if( psInfo->pSrcGeoLocTransformArg == NULL )
1183 : {
1184 0 : GDALDestroyGenImgProjTransformer( psInfo );
1185 0 : return NULL;
1186 : }
1187 2 : if( pszSrcWKT == NULL )
1188 2 : pszSrcWKT = CSLFetchNameValue( papszMD, "SRS" );
1189 : }
1190 :
1191 0 : else if( pszMethod != NULL )
1192 : {
1193 : CPLError( CE_Failure, CPLE_AppDefined,
1194 : "Unable to compute a %s based transformation between pixel/line\n"
1195 : "and georeferenced coordinates for %s.\n",
1196 0 : pszMethod, GDALGetDescription( hSrcDS ) );
1197 :
1198 0 : GDALDestroyGenImgProjTransformer( psInfo );
1199 0 : return NULL;
1200 : }
1201 :
1202 : else
1203 : {
1204 : CPLError( CE_Failure, CPLE_AppDefined,
1205 : "Unable to compute a transformation between pixel/line\n"
1206 : "and georeferenced coordinates for %s.\n"
1207 : "There is no affine transformation and no GCPs.",
1208 0 : GDALGetDescription( hSrcDS ) );
1209 :
1210 0 : GDALDestroyGenImgProjTransformer( psInfo );
1211 0 : return NULL;
1212 : }
1213 :
1214 : /* -------------------------------------------------------------------- */
1215 : /* Setup reprojection. */
1216 : /* -------------------------------------------------------------------- */
1217 670 : if( pszDstWKT == NULL && hDstDS != NULL )
1218 36 : pszDstWKT = GDALGetProjectionRef( hDstDS );
1219 :
1220 670 : if( pszSrcWKT != NULL && strlen(pszSrcWKT) > 0
1221 : && pszDstWKT != NULL && strlen(pszDstWKT) > 0
1222 : && !EQUAL(pszSrcWKT,pszDstWKT) )
1223 : {
1224 90 : CPLString osSrcWKT = pszSrcWKT;
1225 90 : if (hSrcDS
1226 : && CSLFetchBoolean( papszOptions, "INSERT_CENTER_LONG", TRUE ) )
1227 90 : osSrcWKT = InsertCenterLong( hSrcDS, osSrcWKT );
1228 :
1229 : psInfo->pReprojectArg =
1230 90 : GDALCreateReprojectionTransformer( osSrcWKT.c_str(), pszDstWKT );
1231 : }
1232 :
1233 : /* -------------------------------------------------------------------- */
1234 : /* Get forward and inverse geotransform for destination image. */
1235 : /* If we have no destination use a unit transform. */
1236 : /* -------------------------------------------------------------------- */
1237 670 : if( hDstDS )
1238 : {
1239 76 : GDALGetGeoTransform( hDstDS, psInfo->adfDstGeoTransform );
1240 : GDALInvGeoTransform( psInfo->adfDstGeoTransform,
1241 76 : psInfo->adfDstInvGeoTransform );
1242 : }
1243 : else
1244 : {
1245 594 : psInfo->adfDstGeoTransform[0] = 0.0;
1246 594 : psInfo->adfDstGeoTransform[1] = 1.0;
1247 594 : psInfo->adfDstGeoTransform[2] = 0.0;
1248 594 : psInfo->adfDstGeoTransform[3] = 0.0;
1249 594 : psInfo->adfDstGeoTransform[4] = 0.0;
1250 594 : psInfo->adfDstGeoTransform[5] = 1.0;
1251 : memcpy( psInfo->adfDstInvGeoTransform, psInfo->adfDstGeoTransform,
1252 594 : sizeof(double) * 6 );
1253 : }
1254 :
1255 670 : return psInfo;
1256 : }
1257 :
1258 : /************************************************************************/
1259 : /* GDALRefreshGenImgProjTransformer() */
1260 : /************************************************************************/
1261 :
1262 574 : void GDALRefreshGenImgProjTransformer(void* hTransformArg)
1263 : {
1264 : GDALGenImgProjTransformInfo *psInfo =
1265 574 : static_cast<GDALGenImgProjTransformInfo *>( hTransformArg );
1266 :
1267 574 : if (psInfo->pReprojectArg)
1268 : {
1269 56 : CPLXMLNode* psXML = GDALSerializeReprojectionTransformer(psInfo->pReprojectArg);
1270 56 : GDALDestroyReprojectionTransformer(psInfo->pReprojectArg);
1271 56 : psInfo->pReprojectArg = GDALDeserializeReprojectionTransformer(psXML);
1272 56 : CPLDestroyXMLNode(psXML);
1273 : }
1274 574 : }
1275 :
1276 : /************************************************************************/
1277 : /* GDALCreateGenImgProjTransformer3() */
1278 : /************************************************************************/
1279 :
1280 : /**
1281 : * Create image to image transformer.
1282 : *
1283 : * This function creates a transformation object that maps from pixel/line
1284 : * coordinates on one image to pixel/line coordinates on another image. The
1285 : * images may potentially be georeferenced in different coordinate systems,
1286 : * and may used GCPs to map between their pixel/line coordinates and
1287 : * georeferenced coordinates (as opposed to the default assumption that their
1288 : * geotransform should be used).
1289 : *
1290 : * This transformer potentially performs three concatenated transformations.
1291 : *
1292 : * The first stage is from source image pixel/line coordinates to source
1293 : * image georeferenced coordinates, and may be done using the geotransform,
1294 : * or if not defined using a polynomial model derived from GCPs. If GCPs
1295 : * are used this stage is accomplished using GDALGCPTransform().
1296 : *
1297 : * The second stage is to change projections from the source coordinate system
1298 : * to the destination coordinate system, assuming they differ. This is
1299 : * accomplished internally using GDALReprojectionTransform().
1300 : *
1301 : * The third stage is converting from destination image georeferenced
1302 : * coordinates to destination image coordinates. This is done using the
1303 : * destination image geotransform, or if not available, using a polynomial
1304 : * model derived from GCPs. If GCPs are used this stage is accomplished using
1305 : * GDALGCPTransform(). This stage is skipped if hDstDS is NULL when the
1306 : * transformation is created.
1307 : *
1308 : * @param pszSrcWKT source WKT (or NULL).
1309 : * @param padfSrcGeoTransform source geotransform (or NULL).
1310 : * @param pszDstWKT destination WKT (or NULL).
1311 : * @param padfDstGeoTransform destination geotransform (or NULL).
1312 : *
1313 : * @return handle suitable for use GDALGenImgProjTransform(), and to be
1314 : * deallocated with GDALDestroyGenImgProjTransformer() or NULL on failure.
1315 : */
1316 :
1317 : void *
1318 0 : GDALCreateGenImgProjTransformer3( const char *pszSrcWKT,
1319 : const double *padfSrcGeoTransform,
1320 : const char *pszDstWKT,
1321 : const double *padfDstGeoTransform )
1322 :
1323 : {
1324 : GDALGenImgProjTransformInfo *psInfo;
1325 :
1326 : /* -------------------------------------------------------------------- */
1327 : /* Initialize the transform info. */
1328 : /* -------------------------------------------------------------------- */
1329 : psInfo = (GDALGenImgProjTransformInfo *)
1330 0 : CPLCalloc(sizeof(GDALGenImgProjTransformInfo),1);
1331 :
1332 0 : strcpy( psInfo->sTI.szSignature, "GTI" );
1333 0 : psInfo->sTI.pszClassName = "GDALGenImgProjTransformer";
1334 0 : psInfo->sTI.pfnTransform = GDALGenImgProjTransform;
1335 0 : psInfo->sTI.pfnCleanup = GDALDestroyGenImgProjTransformer;
1336 0 : psInfo->sTI.pfnSerialize = GDALSerializeGenImgProjTransformer;
1337 :
1338 : /* -------------------------------------------------------------------- */
1339 : /* Get forward and inverse geotransform for the source image. */
1340 : /* -------------------------------------------------------------------- */
1341 0 : if( padfSrcGeoTransform )
1342 : {
1343 : memcpy( psInfo->adfSrcGeoTransform, padfSrcGeoTransform,
1344 0 : sizeof(psInfo->adfSrcGeoTransform) );
1345 : GDALInvGeoTransform( psInfo->adfSrcGeoTransform,
1346 0 : psInfo->adfSrcInvGeoTransform );
1347 : }
1348 : else
1349 : {
1350 0 : psInfo->adfSrcGeoTransform[0] = 0.0;
1351 0 : psInfo->adfSrcGeoTransform[1] = 1.0;
1352 0 : psInfo->adfSrcGeoTransform[2] = 0.0;
1353 0 : psInfo->adfSrcGeoTransform[3] = 0.0;
1354 0 : psInfo->adfSrcGeoTransform[4] = 0.0;
1355 0 : psInfo->adfSrcGeoTransform[5] = 1.0;
1356 : memcpy( psInfo->adfSrcInvGeoTransform, psInfo->adfSrcGeoTransform,
1357 0 : sizeof(double) * 6 );
1358 : }
1359 :
1360 : /* -------------------------------------------------------------------- */
1361 : /* Setup reprojection. */
1362 : /* -------------------------------------------------------------------- */
1363 0 : if( pszSrcWKT != NULL && strlen(pszSrcWKT) > 0
1364 : && pszDstWKT != NULL && strlen(pszDstWKT) > 0
1365 : && !EQUAL(pszSrcWKT, pszDstWKT) )
1366 : {
1367 : psInfo->pReprojectArg =
1368 0 : GDALCreateReprojectionTransformer( pszSrcWKT, pszDstWKT );
1369 : }
1370 :
1371 : /* -------------------------------------------------------------------- */
1372 : /* Get forward and inverse geotransform for destination image. */
1373 : /* If we have no destination matrix use a unit transform. */
1374 : /* -------------------------------------------------------------------- */
1375 0 : if( padfDstGeoTransform )
1376 : {
1377 : memcpy( psInfo->adfDstGeoTransform, padfDstGeoTransform,
1378 0 : sizeof(psInfo->adfDstGeoTransform) );
1379 : GDALInvGeoTransform( psInfo->adfDstGeoTransform,
1380 0 : psInfo->adfDstInvGeoTransform );
1381 : }
1382 : else
1383 : {
1384 0 : psInfo->adfDstGeoTransform[0] = 0.0;
1385 0 : psInfo->adfDstGeoTransform[1] = 1.0;
1386 0 : psInfo->adfDstGeoTransform[2] = 0.0;
1387 0 : psInfo->adfDstGeoTransform[3] = 0.0;
1388 0 : psInfo->adfDstGeoTransform[4] = 0.0;
1389 0 : psInfo->adfDstGeoTransform[5] = 1.0;
1390 : memcpy( psInfo->adfDstInvGeoTransform, psInfo->adfDstGeoTransform,
1391 0 : sizeof(double) * 6 );
1392 : }
1393 :
1394 0 : return psInfo;
1395 : }
1396 :
1397 : /************************************************************************/
1398 : /* GDALSetGenImgProjTransformerDstGeoTransform() */
1399 : /************************************************************************/
1400 :
1401 : /**
1402 : * Set GenImgProj output geotransform.
1403 : *
1404 : * Normally the "destination geotransform", or transformation between
1405 : * georeferenced output coordinates and pixel/line coordinates on the
1406 : * destination file is extracted from the destination file by
1407 : * GDALCreateGenImgProjTransformer() and stored in the GenImgProj private
1408 : * info. However, sometimes it is inconvenient to have an output file
1409 : * handle with appropriate geotransform information when creating the
1410 : * transformation. For these cases, this function can be used to apply
1411 : * the destination geotransform.
1412 : *
1413 : * @param hTransformArg the handle to update.
1414 : * @param padfGeoTransform the destination geotransform to apply (six doubles).
1415 : */
1416 :
1417 560 : void GDALSetGenImgProjTransformerDstGeoTransform(
1418 : void *hTransformArg, const double *padfGeoTransform )
1419 :
1420 : {
1421 560 : VALIDATE_POINTER0( hTransformArg, "GDALSetGenImgProjTransformerDstGeoTransform" );
1422 :
1423 : GDALGenImgProjTransformInfo *psInfo =
1424 560 : static_cast<GDALGenImgProjTransformInfo *>( hTransformArg );
1425 :
1426 560 : memcpy( psInfo->adfDstGeoTransform, padfGeoTransform, sizeof(double) * 6 );
1427 : GDALInvGeoTransform( psInfo->adfDstGeoTransform,
1428 560 : psInfo->adfDstInvGeoTransform );
1429 : }
1430 :
1431 : /************************************************************************/
1432 : /* GDALDestroyGenImgProjTransformer() */
1433 : /************************************************************************/
1434 :
1435 : /**
1436 : * GenImgProjTransformer deallocator.
1437 : *
1438 : * This function is used to deallocate the handle created with
1439 : * GDALCreateGenImgProjTransformer().
1440 : *
1441 : * @param hTransformArg the handle to deallocate.
1442 : */
1443 :
1444 754 : void GDALDestroyGenImgProjTransformer( void *hTransformArg )
1445 :
1446 : {
1447 754 : VALIDATE_POINTER0( hTransformArg, "GDALDestroyGenImgProjTransformer" );
1448 :
1449 : GDALGenImgProjTransformInfo *psInfo =
1450 754 : (GDALGenImgProjTransformInfo *) hTransformArg;
1451 :
1452 754 : if( psInfo->pSrcGCPTransformArg != NULL )
1453 42 : GDALDestroyGCPTransformer( psInfo->pSrcGCPTransformArg );
1454 :
1455 754 : if( psInfo->pSrcTPSTransformArg != NULL )
1456 8 : GDALDestroyTPSTransformer( psInfo->pSrcTPSTransformArg );
1457 :
1458 754 : if( psInfo->pSrcRPCTransformArg != NULL )
1459 6 : GDALDestroyRPCTransformer( psInfo->pSrcRPCTransformArg );
1460 :
1461 754 : if( psInfo->pSrcGeoLocTransformArg != NULL )
1462 4 : GDALDestroyGeoLocTransformer( psInfo->pSrcGeoLocTransformArg );
1463 :
1464 754 : if( psInfo->pDstGCPTransformArg != NULL )
1465 0 : GDALDestroyGCPTransformer( psInfo->pDstGCPTransformArg );
1466 :
1467 754 : if( psInfo->pReprojectArg != NULL )
1468 108 : GDALDestroyReprojectionTransformer( psInfo->pReprojectArg );
1469 :
1470 754 : CPLFree( psInfo );
1471 : }
1472 :
1473 : /************************************************************************/
1474 : /* GDALGenImgProjTransform() */
1475 : /************************************************************************/
1476 :
1477 : /**
1478 : * Perform general image reprojection transformation.
1479 : *
1480 : * Actually performs the transformation setup in
1481 : * GDALCreateGenImgProjTransformer(). This function matches the signature
1482 : * required by the GDALTransformerFunc(), and more details on the arguments
1483 : * can be found in that topic.
1484 : */
1485 :
1486 601348 : int GDALGenImgProjTransform( void *pTransformArg, int bDstToSrc,
1487 : int nPointCount,
1488 : double *padfX, double *padfY, double *padfZ,
1489 : int *panSuccess )
1490 : {
1491 : GDALGenImgProjTransformInfo *psInfo =
1492 601348 : (GDALGenImgProjTransformInfo *) pTransformArg;
1493 : int i;
1494 : double *padfGeoTransform;
1495 : void *pGCPTransformArg;
1496 : void *pRPCTransformArg;
1497 : void *pTPSTransformArg;
1498 : void *pGeoLocTransformArg;
1499 :
1500 : /* -------------------------------------------------------------------- */
1501 : /* Convert from src (dst) pixel/line to src (dst) */
1502 : /* georeferenced coordinates. */
1503 : /* -------------------------------------------------------------------- */
1504 601348 : if( bDstToSrc )
1505 : {
1506 357434 : padfGeoTransform = psInfo->adfDstGeoTransform;
1507 357434 : pGCPTransformArg = psInfo->pDstGCPTransformArg;
1508 357434 : pRPCTransformArg = NULL;
1509 357434 : pTPSTransformArg = NULL;
1510 357434 : pGeoLocTransformArg = NULL;
1511 : }
1512 : else
1513 : {
1514 243914 : padfGeoTransform = psInfo->adfSrcGeoTransform;
1515 243914 : pGCPTransformArg = psInfo->pSrcGCPTransformArg;
1516 243914 : pRPCTransformArg = psInfo->pSrcRPCTransformArg;
1517 243914 : pTPSTransformArg = psInfo->pSrcTPSTransformArg;
1518 243914 : pGeoLocTransformArg = psInfo->pSrcGeoLocTransformArg;
1519 : }
1520 :
1521 601348 : if( pGCPTransformArg != NULL )
1522 : {
1523 11602 : if( !GDALGCPTransform( pGCPTransformArg, FALSE,
1524 : nPointCount, padfX, padfY, padfZ,
1525 : panSuccess ) )
1526 0 : return FALSE;
1527 : }
1528 589746 : else if( pTPSTransformArg != NULL )
1529 : {
1530 1778 : if( !GDALTPSTransform( pTPSTransformArg, FALSE,
1531 : nPointCount, padfX, padfY, padfZ,
1532 : panSuccess ) )
1533 0 : return FALSE;
1534 : }
1535 587968 : else if( pRPCTransformArg != NULL )
1536 : {
1537 8 : if( !GDALRPCTransform( pRPCTransformArg, FALSE,
1538 : nPointCount, padfX, padfY, padfZ,
1539 : panSuccess ) )
1540 0 : return FALSE;
1541 : }
1542 587960 : else if( pGeoLocTransformArg != NULL )
1543 : {
1544 2 : if( !GDALGeoLocTransform( pGeoLocTransformArg, FALSE,
1545 : nPointCount, padfX, padfY, padfZ,
1546 : panSuccess ) )
1547 0 : return FALSE;
1548 : }
1549 : else
1550 : {
1551 17564906 : for( i = 0; i < nPointCount; i++ )
1552 : {
1553 : double dfNewX, dfNewY;
1554 :
1555 16976948 : if( padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL )
1556 : {
1557 1314 : panSuccess[i] = FALSE;
1558 1314 : continue;
1559 : }
1560 :
1561 16975634 : dfNewX = padfGeoTransform[0]
1562 16975634 : + padfX[i] * padfGeoTransform[1]
1563 16975634 : + padfY[i] * padfGeoTransform[2];
1564 16975634 : dfNewY = padfGeoTransform[3]
1565 16975634 : + padfX[i] * padfGeoTransform[4]
1566 16975634 : + padfY[i] * padfGeoTransform[5];
1567 :
1568 16975634 : padfX[i] = dfNewX;
1569 16975634 : padfY[i] = dfNewY;
1570 : }
1571 : }
1572 :
1573 : /* -------------------------------------------------------------------- */
1574 : /* Reproject if needed. */
1575 : /* -------------------------------------------------------------------- */
1576 601348 : if( psInfo->pReprojectArg )
1577 : {
1578 73290 : if( !GDALReprojectionTransform( psInfo->pReprojectArg, bDstToSrc,
1579 : nPointCount, padfX, padfY, padfZ,
1580 : panSuccess ) )
1581 0 : return FALSE;
1582 : }
1583 : else
1584 : {
1585 7986644 : for( i = 0; i < nPointCount; i++ )
1586 7458586 : panSuccess[i] = 1;
1587 : }
1588 :
1589 : /* -------------------------------------------------------------------- */
1590 : /* Convert dst (src) georef coordinates back to pixel/line. */
1591 : /* -------------------------------------------------------------------- */
1592 601348 : if( bDstToSrc )
1593 : {
1594 357434 : padfGeoTransform = psInfo->adfSrcInvGeoTransform;
1595 357434 : pGCPTransformArg = psInfo->pSrcGCPTransformArg;
1596 357434 : pRPCTransformArg = psInfo->pSrcRPCTransformArg;
1597 357434 : pTPSTransformArg = psInfo->pSrcTPSTransformArg;
1598 357434 : pGeoLocTransformArg = psInfo->pSrcGeoLocTransformArg;
1599 : }
1600 : else
1601 : {
1602 243914 : padfGeoTransform = psInfo->adfDstInvGeoTransform;
1603 243914 : pGCPTransformArg = psInfo->pDstGCPTransformArg;
1604 243914 : pRPCTransformArg = NULL;
1605 243914 : pTPSTransformArg = NULL;
1606 243914 : pGeoLocTransformArg = NULL;
1607 : }
1608 :
1609 601348 : if( pGCPTransformArg != NULL )
1610 : {
1611 12812 : if( !GDALGCPTransform( pGCPTransformArg, TRUE,
1612 : nPointCount, padfX, padfY, padfZ,
1613 : panSuccess ) )
1614 0 : return FALSE;
1615 : }
1616 588536 : else if( pTPSTransformArg != NULL )
1617 : {
1618 1862 : if( !GDALTPSTransform( pTPSTransformArg, TRUE,
1619 : nPointCount, padfX, padfY, padfZ,
1620 : panSuccess ) )
1621 0 : return FALSE;
1622 : }
1623 586674 : else if( pRPCTransformArg != NULL )
1624 : {
1625 8 : if( !GDALRPCTransform( pRPCTransformArg, TRUE,
1626 : nPointCount, padfX, padfY, padfZ,
1627 : panSuccess ) )
1628 0 : return FALSE;
1629 : }
1630 586666 : else if( pGeoLocTransformArg != NULL )
1631 : {
1632 518 : if( !GDALGeoLocTransform( pGeoLocTransformArg, TRUE,
1633 : nPointCount, padfX, padfY, padfZ,
1634 : panSuccess ) )
1635 0 : return FALSE;
1636 : }
1637 : else
1638 : {
1639 17283772 : for( i = 0; i < nPointCount; i++ )
1640 : {
1641 : double dfNewX, dfNewY;
1642 :
1643 16697624 : if( !panSuccess[i] )
1644 319966 : continue;
1645 :
1646 16377658 : dfNewX = padfGeoTransform[0]
1647 16377658 : + padfX[i] * padfGeoTransform[1]
1648 16377658 : + padfY[i] * padfGeoTransform[2];
1649 16377658 : dfNewY = padfGeoTransform[3]
1650 16377658 : + padfX[i] * padfGeoTransform[4]
1651 16377658 : + padfY[i] * padfGeoTransform[5];
1652 :
1653 16377658 : padfX[i] = dfNewX;
1654 16377658 : padfY[i] = dfNewY;
1655 : }
1656 : }
1657 :
1658 601348 : return TRUE;
1659 : }
1660 :
1661 : /************************************************************************/
1662 : /* GDALSerializeGenImgProjTransformer() */
1663 : /************************************************************************/
1664 :
1665 : static CPLXMLNode *
1666 12 : GDALSerializeGenImgProjTransformer( void *pTransformArg )
1667 :
1668 : {
1669 : char szWork[200];
1670 : CPLXMLNode *psTree;
1671 : GDALGenImgProjTransformInfo *psInfo =
1672 12 : (GDALGenImgProjTransformInfo *) pTransformArg;
1673 :
1674 12 : psTree = CPLCreateXMLNode( NULL, CXT_Element, "GenImgProjTransformer" );
1675 :
1676 : /* -------------------------------------------------------------------- */
1677 : /* Handle GCP transformation. */
1678 : /* -------------------------------------------------------------------- */
1679 12 : if( psInfo->pSrcGCPTransformArg != NULL )
1680 : {
1681 : CPLXMLNode *psTransformerContainer;
1682 : CPLXMLNode *psTransformer;
1683 :
1684 : psTransformerContainer =
1685 6 : CPLCreateXMLNode( psTree, CXT_Element, "SrcGCPTransformer" );
1686 :
1687 : psTransformer = GDALSerializeTransformer( GDALGCPTransform,
1688 6 : psInfo->pSrcGCPTransformArg);
1689 6 : if( psTransformer != NULL )
1690 6 : CPLAddXMLChild( psTransformerContainer, psTransformer );
1691 : }
1692 :
1693 : /* -------------------------------------------------------------------- */
1694 : /* Handle TPS transformation. */
1695 : /* -------------------------------------------------------------------- */
1696 6 : else if( psInfo->pSrcTPSTransformArg != NULL )
1697 : {
1698 : CPLXMLNode *psTransformerContainer;
1699 : CPLXMLNode *psTransformer;
1700 :
1701 : psTransformerContainer =
1702 2 : CPLCreateXMLNode( psTree, CXT_Element, "SrcTPSTransformer" );
1703 :
1704 : psTransformer =
1705 2 : GDALSerializeTransformer( NULL, psInfo->pSrcTPSTransformArg);
1706 2 : if( psTransformer != NULL )
1707 2 : CPLAddXMLChild( psTransformerContainer, psTransformer );
1708 : }
1709 :
1710 : /* -------------------------------------------------------------------- */
1711 : /* Handle GeoLoc transformation. */
1712 : /* -------------------------------------------------------------------- */
1713 4 : else if( psInfo->pSrcGeoLocTransformArg != NULL )
1714 : {
1715 : CPLXMLNode *psTransformerContainer;
1716 : CPLXMLNode *psTransformer;
1717 :
1718 : psTransformerContainer =
1719 0 : CPLCreateXMLNode( psTree, CXT_Element, "SrcGeoLocTransformer" );
1720 :
1721 : psTransformer =
1722 0 : GDALSerializeTransformer( NULL, psInfo->pSrcGeoLocTransformArg);
1723 0 : if( psTransformer != NULL )
1724 0 : CPLAddXMLChild( psTransformerContainer, psTransformer );
1725 : }
1726 :
1727 : /* -------------------------------------------------------------------- */
1728 : /* Handle RPC transformation. */
1729 : /* -------------------------------------------------------------------- */
1730 4 : else if( psInfo->pSrcRPCTransformArg != NULL )
1731 : {
1732 : CPLXMLNode *psTransformerContainer;
1733 : CPLXMLNode *psTransformer;
1734 :
1735 : psTransformerContainer =
1736 0 : CPLCreateXMLNode( psTree, CXT_Element, "SrcRPCTransformer" );
1737 :
1738 : psTransformer =
1739 0 : GDALSerializeTransformer( NULL, psInfo->pSrcRPCTransformArg);
1740 0 : if( psTransformer != NULL )
1741 0 : CPLAddXMLChild( psTransformerContainer, psTransformer );
1742 : }
1743 :
1744 : /* -------------------------------------------------------------------- */
1745 : /* Handle source geotransforms. */
1746 : /* -------------------------------------------------------------------- */
1747 : else
1748 : {
1749 : sprintf( szWork, "%.16g,%.16g,%.16g,%.16g,%.16g,%.16g",
1750 : psInfo->adfSrcGeoTransform[0],
1751 : psInfo->adfSrcGeoTransform[1],
1752 : psInfo->adfSrcGeoTransform[2],
1753 : psInfo->adfSrcGeoTransform[3],
1754 : psInfo->adfSrcGeoTransform[4],
1755 4 : psInfo->adfSrcGeoTransform[5] );
1756 4 : CPLCreateXMLElementAndValue( psTree, "SrcGeoTransform", szWork );
1757 :
1758 : sprintf( szWork, "%.16g,%.16g,%.16g,%.16g,%.16g,%.16g",
1759 : psInfo->adfSrcInvGeoTransform[0],
1760 : psInfo->adfSrcInvGeoTransform[1],
1761 : psInfo->adfSrcInvGeoTransform[2],
1762 : psInfo->adfSrcInvGeoTransform[3],
1763 : psInfo->adfSrcInvGeoTransform[4],
1764 4 : psInfo->adfSrcInvGeoTransform[5] );
1765 4 : CPLCreateXMLElementAndValue( psTree, "SrcInvGeoTransform", szWork );
1766 : }
1767 :
1768 : /* -------------------------------------------------------------------- */
1769 : /* Handle destination geotransforms. */
1770 : /* -------------------------------------------------------------------- */
1771 : sprintf( szWork, "%.16g,%.16g,%.16g,%.16g,%.16g,%.16g",
1772 : psInfo->adfDstGeoTransform[0],
1773 : psInfo->adfDstGeoTransform[1],
1774 : psInfo->adfDstGeoTransform[2],
1775 : psInfo->adfDstGeoTransform[3],
1776 : psInfo->adfDstGeoTransform[4],
1777 12 : psInfo->adfDstGeoTransform[5] );
1778 12 : CPLCreateXMLElementAndValue( psTree, "DstGeoTransform", szWork );
1779 :
1780 : sprintf( szWork, "%.16g,%.16g,%.16g,%.16g,%.16g,%.16g",
1781 : psInfo->adfDstInvGeoTransform[0],
1782 : psInfo->adfDstInvGeoTransform[1],
1783 : psInfo->adfDstInvGeoTransform[2],
1784 : psInfo->adfDstInvGeoTransform[3],
1785 : psInfo->adfDstInvGeoTransform[4],
1786 12 : psInfo->adfDstInvGeoTransform[5] );
1787 12 : CPLCreateXMLElementAndValue( psTree, "DstInvGeoTransform", szWork );
1788 :
1789 : /* -------------------------------------------------------------------- */
1790 : /* Do we have a reprojection transformer? */
1791 : /* -------------------------------------------------------------------- */
1792 12 : if( psInfo->pReprojectArg != NULL )
1793 : {
1794 : CPLXMLNode *psTransformerContainer;
1795 : CPLXMLNode *psTransformer;
1796 :
1797 : psTransformerContainer =
1798 0 : CPLCreateXMLNode( psTree, CXT_Element, "ReprojectTransformer" );
1799 :
1800 : psTransformer = GDALSerializeTransformer( GDALReprojectionTransform,
1801 0 : psInfo->pReprojectArg );
1802 0 : if( psTransformer != NULL )
1803 0 : CPLAddXMLChild( psTransformerContainer, psTransformer );
1804 : }
1805 :
1806 12 : return psTree;
1807 : }
1808 :
1809 : /************************************************************************/
1810 : /* GDALDeserializeGenImgProjTransformer() */
1811 : /************************************************************************/
1812 :
1813 84 : void *GDALDeserializeGenImgProjTransformer( CPLXMLNode *psTree )
1814 :
1815 : {
1816 : GDALGenImgProjTransformInfo *psInfo;
1817 : CPLXMLNode *psSubtree;
1818 :
1819 : /* -------------------------------------------------------------------- */
1820 : /* Initialize the transform info. */
1821 : /* -------------------------------------------------------------------- */
1822 : psInfo = (GDALGenImgProjTransformInfo *)
1823 84 : CPLCalloc(sizeof(GDALGenImgProjTransformInfo),1);
1824 :
1825 84 : strcpy( psInfo->sTI.szSignature, "GTI" );
1826 84 : psInfo->sTI.pszClassName = "GDALGenImgProjTransformer";
1827 84 : psInfo->sTI.pfnTransform = GDALGenImgProjTransform;
1828 84 : psInfo->sTI.pfnCleanup = GDALDestroyGenImgProjTransformer;
1829 84 : psInfo->sTI.pfnSerialize = GDALSerializeGenImgProjTransformer;
1830 :
1831 : /* -------------------------------------------------------------------- */
1832 : /* SrcGeotransform */
1833 : /* -------------------------------------------------------------------- */
1834 84 : if( CPLGetXMLNode( psTree, "SrcGeoTransform" ) != NULL )
1835 : {
1836 : sscanf( CPLGetXMLValue( psTree, "SrcGeoTransform", "" ),
1837 : "%lg,%lg,%lg,%lg,%lg,%lg",
1838 : psInfo->adfSrcGeoTransform + 0,
1839 : psInfo->adfSrcGeoTransform + 1,
1840 : psInfo->adfSrcGeoTransform + 2,
1841 : psInfo->adfSrcGeoTransform + 3,
1842 : psInfo->adfSrcGeoTransform + 4,
1843 70 : psInfo->adfSrcGeoTransform + 5 );
1844 :
1845 70 : if( CPLGetXMLNode( psTree, "SrcInvGeoTransform" ) != NULL )
1846 : {
1847 : sscanf( CPLGetXMLValue( psTree, "SrcInvGeoTransform", "" ),
1848 : "%lg,%lg,%lg,%lg,%lg,%lg",
1849 : psInfo->adfSrcInvGeoTransform + 0,
1850 : psInfo->adfSrcInvGeoTransform + 1,
1851 : psInfo->adfSrcInvGeoTransform + 2,
1852 : psInfo->adfSrcInvGeoTransform + 3,
1853 : psInfo->adfSrcInvGeoTransform + 4,
1854 70 : psInfo->adfSrcInvGeoTransform + 5 );
1855 :
1856 : }
1857 : else
1858 : GDALInvGeoTransform( psInfo->adfSrcGeoTransform,
1859 0 : psInfo->adfSrcInvGeoTransform );
1860 : }
1861 :
1862 : /* -------------------------------------------------------------------- */
1863 : /* Src GCP Transform */
1864 : /* -------------------------------------------------------------------- */
1865 84 : psSubtree = CPLGetXMLNode( psTree, "SrcGCPTransformer" );
1866 84 : if( psSubtree != NULL && psSubtree->psChild != NULL )
1867 : {
1868 : psInfo->pSrcGCPTransformArg =
1869 10 : GDALDeserializeGCPTransformer( psSubtree->psChild );
1870 : }
1871 :
1872 : /* -------------------------------------------------------------------- */
1873 : /* Src TPS Transform */
1874 : /* -------------------------------------------------------------------- */
1875 84 : psSubtree = CPLGetXMLNode( psTree, "SrcTPSTransformer" );
1876 84 : if( psSubtree != NULL && psSubtree->psChild != NULL )
1877 : {
1878 : psInfo->pSrcTPSTransformArg =
1879 2 : GDALDeserializeTPSTransformer( psSubtree->psChild );
1880 : }
1881 :
1882 : /* -------------------------------------------------------------------- */
1883 : /* Src GeoLoc Transform */
1884 : /* -------------------------------------------------------------------- */
1885 84 : psSubtree = CPLGetXMLNode( psTree, "SrcGeoLocTransformer" );
1886 84 : if( psSubtree != NULL && psSubtree->psChild != NULL )
1887 : {
1888 : psInfo->pSrcGeoLocTransformArg =
1889 2 : GDALDeserializeGeoLocTransformer( psSubtree->psChild );
1890 : }
1891 :
1892 : /* -------------------------------------------------------------------- */
1893 : /* Src RPC Transform */
1894 : /* -------------------------------------------------------------------- */
1895 84 : psSubtree = CPLGetXMLNode( psTree, "SrcRPCTransformer" );
1896 84 : if( psSubtree != NULL && psSubtree->psChild != NULL )
1897 : {
1898 : psInfo->pSrcRPCTransformArg =
1899 0 : GDALDeserializeRPCTransformer( psSubtree->psChild );
1900 : }
1901 :
1902 : /* -------------------------------------------------------------------- */
1903 : /* DstGeotransform */
1904 : /* -------------------------------------------------------------------- */
1905 84 : if( CPLGetXMLNode( psTree, "DstGeoTransform" ) != NULL )
1906 : {
1907 : sscanf( CPLGetXMLValue( psTree, "DstGeoTransform", "" ),
1908 : "%lg,%lg,%lg,%lg,%lg,%lg",
1909 : psInfo->adfDstGeoTransform + 0,
1910 : psInfo->adfDstGeoTransform + 1,
1911 : psInfo->adfDstGeoTransform + 2,
1912 : psInfo->adfDstGeoTransform + 3,
1913 : psInfo->adfDstGeoTransform + 4,
1914 84 : psInfo->adfDstGeoTransform + 5 );
1915 :
1916 84 : if( CPLGetXMLNode( psTree, "DstInvGeoTransform" ) != NULL )
1917 : {
1918 : sscanf( CPLGetXMLValue( psTree, "DstInvGeoTransform", "" ),
1919 : "%lg,%lg,%lg,%lg,%lg,%lg",
1920 : psInfo->adfDstInvGeoTransform + 0,
1921 : psInfo->adfDstInvGeoTransform + 1,
1922 : psInfo->adfDstInvGeoTransform + 2,
1923 : psInfo->adfDstInvGeoTransform + 3,
1924 : psInfo->adfDstInvGeoTransform + 4,
1925 52 : psInfo->adfDstInvGeoTransform + 5 );
1926 :
1927 : }
1928 : else
1929 : GDALInvGeoTransform( psInfo->adfDstGeoTransform,
1930 32 : psInfo->adfDstInvGeoTransform );
1931 : }
1932 :
1933 : /* -------------------------------------------------------------------- */
1934 : /* Reproject transformer */
1935 : /* -------------------------------------------------------------------- */
1936 84 : psSubtree = CPLGetXMLNode( psTree, "ReprojectTransformer" );
1937 84 : if( psSubtree != NULL && psSubtree->psChild != NULL )
1938 : {
1939 : psInfo->pReprojectArg =
1940 18 : GDALDeserializeReprojectionTransformer( psSubtree->psChild );
1941 : }
1942 :
1943 84 : return psInfo;
1944 : }
1945 :
1946 : /************************************************************************/
1947 : /* ==================================================================== */
1948 : /* GDALReprojectionTransformer */
1949 : /* ==================================================================== */
1950 : /************************************************************************/
1951 :
1952 : typedef struct {
1953 : GDALTransformerInfo sTI;
1954 :
1955 : OGRCoordinateTransformation *poForwardTransform;
1956 : OGRCoordinateTransformation *poReverseTransform;
1957 : } GDALReprojectionTransformInfo;
1958 :
1959 : /************************************************************************/
1960 : /* GDALCreateReprojectionTransformer() */
1961 : /************************************************************************/
1962 :
1963 : /**
1964 : * Create reprojection transformer.
1965 : *
1966 : * Creates a callback data structure suitable for use with
1967 : * GDALReprojectionTransformation() to represent a transformation from
1968 : * one geographic or projected coordinate system to another. On input
1969 : * the coordinate systems are described in OpenGIS WKT format.
1970 : *
1971 : * Internally the OGRCoordinateTransformation object is used to implement
1972 : * the reprojection.
1973 : *
1974 : * @param pszSrcWKT the coordinate system for the source coordinate system.
1975 : * @param pszDstWKT the coordinate system for the destination coordinate
1976 : * system.
1977 : *
1978 : * @return Handle for use with GDALReprojectionTransform(), or NULL if the
1979 : * system fails to initialize the reprojection.
1980 : **/
1981 :
1982 164 : void *GDALCreateReprojectionTransformer( const char *pszSrcWKT,
1983 : const char *pszDstWKT )
1984 :
1985 : {
1986 164 : OGRSpatialReference oSrcSRS, oDstSRS;
1987 : OGRCoordinateTransformation *poForwardTransform;
1988 :
1989 : /* -------------------------------------------------------------------- */
1990 : /* Ingest the SRS definitions. */
1991 : /* -------------------------------------------------------------------- */
1992 164 : if( oSrcSRS.importFromWkt( (char **) &pszSrcWKT ) != OGRERR_NONE )
1993 : {
1994 : CPLError( CE_Failure, CPLE_AppDefined,
1995 : "Failed to import coordinate system `%s'.",
1996 0 : pszSrcWKT );
1997 0 : return NULL;
1998 : }
1999 164 : if( oDstSRS.importFromWkt( (char **) &pszDstWKT ) != OGRERR_NONE )
2000 : {
2001 : CPLError( CE_Failure, CPLE_AppDefined,
2002 : "Failed to import coordinate system `%s'.",
2003 0 : pszSrcWKT );
2004 0 : return NULL;
2005 : }
2006 :
2007 : /* -------------------------------------------------------------------- */
2008 : /* Build the forward coordinate transformation. */
2009 : /* -------------------------------------------------------------------- */
2010 164 : poForwardTransform = OGRCreateCoordinateTransformation(&oSrcSRS,&oDstSRS);
2011 :
2012 164 : if( poForwardTransform == NULL )
2013 : // OGRCreateCoordinateTransformation() will report errors on its own.
2014 0 : return NULL;
2015 :
2016 : /* -------------------------------------------------------------------- */
2017 : /* Create a structure to hold the transform info, and also */
2018 : /* build reverse transform. We assume that if the forward */
2019 : /* transform can be created, then so can the reverse one. */
2020 : /* -------------------------------------------------------------------- */
2021 : GDALReprojectionTransformInfo *psInfo;
2022 :
2023 : psInfo = (GDALReprojectionTransformInfo *)
2024 164 : CPLCalloc(sizeof(GDALReprojectionTransformInfo),1);
2025 :
2026 164 : psInfo->poForwardTransform = poForwardTransform;
2027 : psInfo->poReverseTransform =
2028 164 : OGRCreateCoordinateTransformation(&oDstSRS,&oSrcSRS);
2029 :
2030 164 : strcpy( psInfo->sTI.szSignature, "GTI" );
2031 164 : psInfo->sTI.pszClassName = "GDALReprojectionTransformer";
2032 164 : psInfo->sTI.pfnTransform = GDALReprojectionTransform;
2033 164 : psInfo->sTI.pfnCleanup = GDALDestroyReprojectionTransformer;
2034 164 : psInfo->sTI.pfnSerialize = GDALSerializeReprojectionTransformer;
2035 :
2036 164 : return psInfo;
2037 : }
2038 :
2039 : /************************************************************************/
2040 : /* GDALDestroyReprojectionTransformer() */
2041 : /************************************************************************/
2042 :
2043 : /**
2044 : * Destroy reprojection transformation.
2045 : *
2046 : * @param pTransformArg the transformation handle returned by
2047 : * GDALCreateReprojectionTransformer().
2048 : */
2049 :
2050 164 : void GDALDestroyReprojectionTransformer( void *pTransformArg )
2051 :
2052 : {
2053 164 : VALIDATE_POINTER0( pTransformArg, "GDALDestroyReprojectionTransformer" );
2054 :
2055 : GDALReprojectionTransformInfo *psInfo =
2056 164 : (GDALReprojectionTransformInfo *) pTransformArg;
2057 :
2058 164 : if( psInfo->poForwardTransform )
2059 164 : delete psInfo->poForwardTransform;
2060 :
2061 164 : if( psInfo->poReverseTransform )
2062 164 : delete psInfo->poReverseTransform;
2063 :
2064 164 : CPLFree( psInfo );
2065 : }
2066 :
2067 : /************************************************************************/
2068 : /* GDALReprojectionTransform() */
2069 : /************************************************************************/
2070 :
2071 : /**
2072 : * Perform reprojection transformation.
2073 : *
2074 : * Actually performs the reprojection transformation described in
2075 : * GDALCreateReprojectionTransformer(). This function matches the
2076 : * GDALTransformerFunc() signature. Details of the arguments are described
2077 : * there.
2078 : */
2079 :
2080 73290 : int GDALReprojectionTransform( void *pTransformArg, int bDstToSrc,
2081 : int nPointCount,
2082 : double *padfX, double *padfY, double *padfZ,
2083 : int *panSuccess )
2084 :
2085 : {
2086 : GDALReprojectionTransformInfo *psInfo =
2087 73290 : (GDALReprojectionTransformInfo *) pTransformArg;
2088 : int bSuccess;
2089 :
2090 73290 : if( bDstToSrc )
2091 : bSuccess = psInfo->poReverseTransform->TransformEx(
2092 61034 : nPointCount, padfX, padfY, padfZ, panSuccess );
2093 : else
2094 : bSuccess = psInfo->poForwardTransform->TransformEx(
2095 12256 : nPointCount, padfX, padfY, padfZ, panSuccess );
2096 :
2097 73290 : return bSuccess;
2098 : }
2099 :
2100 : /************************************************************************/
2101 : /* GDALSerializeReprojectionTransformer() */
2102 : /************************************************************************/
2103 :
2104 : static CPLXMLNode *
2105 56 : GDALSerializeReprojectionTransformer( void *pTransformArg )
2106 :
2107 : {
2108 : CPLXMLNode *psTree;
2109 : GDALReprojectionTransformInfo *psInfo =
2110 56 : (GDALReprojectionTransformInfo *) pTransformArg;
2111 :
2112 56 : psTree = CPLCreateXMLNode( NULL, CXT_Element, "ReprojectionTransformer" );
2113 :
2114 : /* -------------------------------------------------------------------- */
2115 : /* Handle SourceCS. */
2116 : /* -------------------------------------------------------------------- */
2117 : OGRSpatialReference *poSRS;
2118 56 : char *pszWKT = NULL;
2119 :
2120 56 : poSRS = psInfo->poForwardTransform->GetSourceCS();
2121 56 : poSRS->exportToWkt( &pszWKT );
2122 56 : CPLCreateXMLElementAndValue( psTree, "SourceSRS", pszWKT );
2123 56 : CPLFree( pszWKT );
2124 :
2125 : /* -------------------------------------------------------------------- */
2126 : /* Handle DestinationCS. */
2127 : /* -------------------------------------------------------------------- */
2128 56 : poSRS = psInfo->poForwardTransform->GetTargetCS();
2129 56 : poSRS->exportToWkt( &pszWKT );
2130 56 : CPLCreateXMLElementAndValue( psTree, "TargetSRS", pszWKT );
2131 56 : CPLFree( pszWKT );
2132 :
2133 56 : return psTree;
2134 : }
2135 :
2136 : /************************************************************************/
2137 : /* GDALDeserializeReprojectionTransformer() */
2138 : /************************************************************************/
2139 :
2140 : static void *
2141 74 : GDALDeserializeReprojectionTransformer( CPLXMLNode *psTree )
2142 :
2143 : {
2144 74 : const char *pszSourceSRS = CPLGetXMLValue( psTree, "SourceSRS", NULL );
2145 74 : const char *pszTargetSRS= CPLGetXMLValue( psTree, "TargetSRS", NULL );
2146 74 : char *pszSourceWKT = NULL, *pszTargetWKT = NULL;
2147 74 : void *pResult = NULL;
2148 :
2149 74 : if( pszSourceSRS != NULL )
2150 : {
2151 74 : OGRSpatialReference oSRS;
2152 :
2153 74 : if( oSRS.SetFromUserInput( pszSourceSRS ) == OGRERR_NONE )
2154 74 : oSRS.exportToWkt( &pszSourceWKT );
2155 : }
2156 :
2157 74 : if( pszTargetSRS != NULL )
2158 : {
2159 74 : OGRSpatialReference oSRS;
2160 :
2161 74 : if( oSRS.SetFromUserInput( pszTargetSRS ) == OGRERR_NONE )
2162 74 : oSRS.exportToWkt( &pszTargetWKT );
2163 : }
2164 :
2165 148 : if( pszSourceWKT != NULL && pszTargetWKT != NULL )
2166 : {
2167 : pResult = GDALCreateReprojectionTransformer( pszSourceWKT,
2168 74 : pszTargetWKT );
2169 : }
2170 : else
2171 : {
2172 : CPLError( CE_Failure, CPLE_AppDefined,
2173 : "ReprojectionTransformer definition missing either\n"
2174 0 : "SourceSRS or TargetSRS definition." );
2175 : }
2176 :
2177 74 : CPLFree( pszSourceWKT );
2178 74 : CPLFree( pszTargetWKT );
2179 :
2180 74 : return pResult;
2181 : }
2182 :
2183 : /************************************************************************/
2184 : /* ==================================================================== */
2185 : /* Approximate transformer. */
2186 : /* ==================================================================== */
2187 : /************************************************************************/
2188 :
2189 : typedef struct
2190 : {
2191 : GDALTransformerInfo sTI;
2192 :
2193 : GDALTransformerFunc pfnBaseTransformer;
2194 : void *pBaseCBData;
2195 : double dfMaxError;
2196 :
2197 : int bOwnSubtransformer;
2198 : } ApproxTransformInfo;
2199 :
2200 : /************************************************************************/
2201 : /* GDALSerializeApproxTransformer() */
2202 : /************************************************************************/
2203 :
2204 : static CPLXMLNode *
2205 6 : GDALSerializeApproxTransformer( void *pTransformArg )
2206 :
2207 : {
2208 : CPLXMLNode *psTree;
2209 6 : ApproxTransformInfo *psInfo = (ApproxTransformInfo *) pTransformArg;
2210 :
2211 6 : psTree = CPLCreateXMLNode( NULL, CXT_Element, "ApproxTransformer" );
2212 :
2213 : /* -------------------------------------------------------------------- */
2214 : /* Attach max error. */
2215 : /* -------------------------------------------------------------------- */
2216 : CPLCreateXMLElementAndValue( psTree, "MaxError",
2217 6 : CPLString().Printf("%g",psInfo->dfMaxError) );
2218 :
2219 : /* -------------------------------------------------------------------- */
2220 : /* Capture underlying transformer. */
2221 : /* -------------------------------------------------------------------- */
2222 : CPLXMLNode *psTransformerContainer;
2223 : CPLXMLNode *psTransformer;
2224 :
2225 : psTransformerContainer =
2226 6 : CPLCreateXMLNode( psTree, CXT_Element, "BaseTransformer" );
2227 :
2228 : psTransformer = GDALSerializeTransformer( psInfo->pfnBaseTransformer,
2229 6 : psInfo->pBaseCBData );
2230 6 : if( psTransformer != NULL )
2231 6 : CPLAddXMLChild( psTransformerContainer, psTransformer );
2232 :
2233 6 : return psTree;
2234 : }
2235 :
2236 : /************************************************************************/
2237 : /* GDALCreateApproxTransformer() */
2238 : /************************************************************************/
2239 :
2240 : /**
2241 : * Create an approximating transformer.
2242 : *
2243 : * This function creates a context for an approximated transformer. Basically
2244 : * a high precision transformer is supplied as input and internally linear
2245 : * approximations are computed to generate results to within a defined
2246 : * precision.
2247 : *
2248 : * The approximation is actually done at the point where GDALApproxTransform()
2249 : * calls are made, and depend on the assumption that the roughly linear. The
2250 : * first and last point passed in must be the extreme values and the
2251 : * intermediate values should describe a curve between the end points. The
2252 : * approximator transforms and center using the approximate transformer, and
2253 : * then compares the true middle transformed value to a linear approximation
2254 : * based on the end points. If the error is within the supplied threshold
2255 : * then the end points are used to linearly approximate all the values
2256 : * otherwise the inputs points are split into two smaller sets, and the
2257 : * function recursively called till a sufficiently small set of points if found
2258 : * that the linear approximation is OK, or that all the points are exactly
2259 : * computed.
2260 : *
2261 : * This function is very suitable for approximating transformation results
2262 : * from output pixel/line space to input coordinates for warpers that operate
2263 : * on one input scanline at a time. Care should be taken using it in other
2264 : * circumstances as little internal validation is done, in order to keep things
2265 : * fast.
2266 : *
2267 : * @param pfnBaseTransformer the high precision transformer which should be
2268 : * approximated.
2269 : * @param pBaseTransformArg the callback argument for the high precision
2270 : * transformer.
2271 : * @param dfMaxError the maximum cartesian error in the "output" space that
2272 : * is to be accepted in the linear approximation.
2273 : *
2274 : * @return callback pointer suitable for use with GDALApproxTransform(). It
2275 : * should be deallocated with GDALDestroyApproxTransformer().
2276 : */
2277 :
2278 582 : void *GDALCreateApproxTransformer( GDALTransformerFunc pfnBaseTransformer,
2279 : void *pBaseTransformArg, double dfMaxError)
2280 :
2281 : {
2282 : ApproxTransformInfo *psATInfo;
2283 :
2284 582 : psATInfo = (ApproxTransformInfo*) CPLMalloc(sizeof(ApproxTransformInfo));
2285 582 : psATInfo->pfnBaseTransformer = pfnBaseTransformer;
2286 582 : psATInfo->pBaseCBData = pBaseTransformArg;
2287 582 : psATInfo->dfMaxError = dfMaxError;
2288 582 : psATInfo->bOwnSubtransformer = FALSE;
2289 :
2290 582 : strcpy( psATInfo->sTI.szSignature, "GTI" );
2291 582 : psATInfo->sTI.pszClassName = "GDALApproxTransformer";
2292 582 : psATInfo->sTI.pfnTransform = GDALApproxTransform;
2293 582 : psATInfo->sTI.pfnCleanup = GDALDestroyApproxTransformer;
2294 582 : psATInfo->sTI.pfnSerialize = GDALSerializeApproxTransformer;
2295 :
2296 582 : return psATInfo;
2297 : }
2298 :
2299 : /************************************************************************/
2300 : /* GDALApproxTransformerOwnsSubtransformer() */
2301 : /************************************************************************/
2302 :
2303 44 : void GDALApproxTransformerOwnsSubtransformer( void *pCBData, int bOwnFlag )
2304 :
2305 : {
2306 44 : ApproxTransformInfo *psATInfo = (ApproxTransformInfo *) pCBData;
2307 :
2308 44 : psATInfo->bOwnSubtransformer = bOwnFlag;
2309 44 : }
2310 :
2311 : /************************************************************************/
2312 : /* GDALDestroyApproxTransformer() */
2313 : /************************************************************************/
2314 :
2315 : /**
2316 : * Cleanup approximate transformer.
2317 : *
2318 : * Deallocates the resources allocated by GDALCreateApproxTransformer().
2319 : *
2320 : * @param pCBData callback data originally returned by
2321 : * GDALCreateApproxTransformer().
2322 : */
2323 :
2324 582 : void GDALDestroyApproxTransformer( void * pCBData )
2325 :
2326 : {
2327 582 : VALIDATE_POINTER0( pCBData, "GDALDestroyApproxTransformer" );
2328 :
2329 582 : ApproxTransformInfo *psATInfo = (ApproxTransformInfo *) pCBData;
2330 :
2331 582 : if( psATInfo->bOwnSubtransformer )
2332 44 : GDALDestroyTransformer( psATInfo->pBaseCBData );
2333 :
2334 582 : CPLFree( pCBData );
2335 : }
2336 :
2337 : /************************************************************************/
2338 : /* GDALApproxTransform() */
2339 : /************************************************************************/
2340 :
2341 : /**
2342 : * Perform approximate transformation.
2343 : *
2344 : * Actually performs the approximate transformation described in
2345 : * GDALCreateApproxTransformer(). This function matches the
2346 : * GDALTransformerFunc() signature. Details of the arguments are described
2347 : * there.
2348 : */
2349 :
2350 95044 : int GDALApproxTransform( void *pCBData, int bDstToSrc, int nPoints,
2351 : double *x, double *y, double *z, int *panSuccess )
2352 :
2353 : {
2354 95044 : ApproxTransformInfo *psATInfo = (ApproxTransformInfo *) pCBData;
2355 : double x2[3], y2[3], z2[3], dfDeltaX, dfDeltaY, dfError, dfDist, dfDeltaZ;
2356 : int nMiddle, anSuccess2[3], i, bSuccess;
2357 :
2358 95044 : nMiddle = (nPoints-1)/2;
2359 :
2360 : /* -------------------------------------------------------------------- */
2361 : /* Bail if our preconditions are not met, or if error is not */
2362 : /* acceptable. */
2363 : /* -------------------------------------------------------------------- */
2364 282314 : if( y[0] != y[nPoints-1] || y[0] != y[nMiddle]
2365 187270 : || x[0] == x[nPoints-1] || x[0] == x[nMiddle]
2366 : || psATInfo->dfMaxError == 0.0 || nPoints <= 5 )
2367 : {
2368 : return psATInfo->pfnBaseTransformer( psATInfo->pBaseCBData, bDstToSrc,
2369 2014 : nPoints, x, y, z, panSuccess );
2370 : }
2371 :
2372 : /* -------------------------------------------------------------------- */
2373 : /* Transform first, last and middle point. */
2374 : /* -------------------------------------------------------------------- */
2375 93030 : x2[0] = x[0];
2376 93030 : y2[0] = y[0];
2377 93030 : z2[0] = z[0];
2378 93030 : x2[1] = x[nMiddle];
2379 93030 : y2[1] = y[nMiddle];
2380 93030 : z2[1] = z[nMiddle];
2381 93030 : x2[2] = x[nPoints-1];
2382 93030 : y2[2] = y[nPoints-1];
2383 93030 : z2[2] = z[nPoints-1];
2384 :
2385 : bSuccess =
2386 : psATInfo->pfnBaseTransformer( psATInfo->pBaseCBData, bDstToSrc, 3,
2387 93030 : x2, y2, z2, anSuccess2 );
2388 93030 : if( !bSuccess || !anSuccess2[0] || !anSuccess2[1] || !anSuccess2[2] )
2389 : return psATInfo->pfnBaseTransformer( psATInfo->pBaseCBData, bDstToSrc,
2390 2892 : nPoints, x, y, z, panSuccess );
2391 :
2392 : /* -------------------------------------------------------------------- */
2393 : /* Is the error at the middle acceptable relative to an */
2394 : /* interpolation of the middle position? */
2395 : /* -------------------------------------------------------------------- */
2396 90138 : dfDeltaX = (x2[2] - x2[0]) / (x[nPoints-1] - x[0]);
2397 90138 : dfDeltaY = (y2[2] - y2[0]) / (x[nPoints-1] - x[0]);
2398 90138 : dfDeltaZ = (z2[2] - z2[0]) / (x[nPoints-1] - x[0]);
2399 :
2400 180276 : dfError = fabs((x2[0] + dfDeltaX * (x[nMiddle] - x[0])) - x2[1])
2401 180276 : + fabs((y2[0] + dfDeltaY * (x[nMiddle] - x[0])) - y2[1]);
2402 :
2403 90138 : if( dfError > psATInfo->dfMaxError )
2404 : {
2405 : #ifdef notdef
2406 : CPLDebug( "GDAL", "ApproxTransformer - "
2407 : "error %g over threshold %g, subdivide %d points.",
2408 : dfError, psATInfo->dfMaxError, nPoints );
2409 : #endif
2410 :
2411 : bSuccess =
2412 : GDALApproxTransform( psATInfo, bDstToSrc, nMiddle,
2413 14590 : x, y, z, panSuccess );
2414 :
2415 14590 : if( !bSuccess )
2416 0 : return FALSE;
2417 :
2418 : bSuccess =
2419 : GDALApproxTransform( psATInfo, bDstToSrc, nPoints - nMiddle,
2420 : x+nMiddle, y+nMiddle, z+nMiddle,
2421 14590 : panSuccess+nMiddle );
2422 :
2423 14590 : if( !bSuccess )
2424 0 : return FALSE;
2425 :
2426 14590 : return TRUE;
2427 : }
2428 :
2429 : /* -------------------------------------------------------------------- */
2430 : /* Error is OK since this is just used to compute output bounds */
2431 : /* of newly created file for gdalwarper. So just use affine */
2432 : /* approximation of the reverse transform. Eventually we */
2433 : /* should implement iterative searching to find a result within */
2434 : /* our error threshold. */
2435 : /* -------------------------------------------------------------------- */
2436 17119840 : for( i = nPoints-1; i >= 0; i-- )
2437 : {
2438 17044292 : dfDist = (x[i] - x[0]);
2439 17044292 : y[i] = y2[0] + dfDeltaY * dfDist;
2440 17044292 : x[i] = x2[0] + dfDeltaX * dfDist;
2441 17044292 : z[i] = z2[0] + dfDeltaZ * dfDist;
2442 17044292 : panSuccess[i] = TRUE;
2443 : }
2444 :
2445 75548 : return TRUE;
2446 : }
2447 :
2448 : /************************************************************************/
2449 : /* GDALDeserializeApproxTransformer() */
2450 : /************************************************************************/
2451 :
2452 : static void *
2453 42 : GDALDeserializeApproxTransformer( CPLXMLNode *psTree )
2454 :
2455 : {
2456 42 : double dfMaxError = atof(CPLGetXMLValue( psTree, "MaxError", "0.25" ));
2457 : CPLXMLNode *psContainer;
2458 42 : GDALTransformerFunc pfnBaseTransform = NULL;
2459 42 : void *pBaseCBData = NULL;
2460 :
2461 42 : psContainer = CPLGetXMLNode( psTree, "BaseTransformer" );
2462 :
2463 42 : if( psContainer != NULL && psContainer->psChild != NULL )
2464 : {
2465 : GDALDeserializeTransformer( psContainer->psChild,
2466 : &pfnBaseTransform,
2467 42 : &pBaseCBData );
2468 :
2469 : }
2470 :
2471 42 : if( pfnBaseTransform == NULL )
2472 : {
2473 : CPLError( CE_Failure, CPLE_AppDefined,
2474 0 : "Cannot get base transform for approx transformer." );
2475 0 : return NULL;
2476 : }
2477 : else
2478 : {
2479 : void *pApproxCBData = GDALCreateApproxTransformer( pfnBaseTransform,
2480 : pBaseCBData,
2481 42 : dfMaxError );
2482 42 : GDALApproxTransformerOwnsSubtransformer( pApproxCBData, TRUE );
2483 :
2484 42 : return pApproxCBData;
2485 : }
2486 : }
2487 :
2488 : /************************************************************************/
2489 : /* GDALApplyGeoTransform() */
2490 : /************************************************************************/
2491 :
2492 : /**
2493 : * Apply GeoTransform to x/y coordinate.
2494 : *
2495 : * Applies the following computation, converting a (pixel,line) coordinate
2496 : * into a georeferenced (geo_x,geo_y) location.
2497 : *
2498 : * *pdfGeoX = padfGeoTransform[0] + dfPixel * padfGeoTransform[1]
2499 : * + dfLine * padfGeoTransform[2];
2500 : * *pdfGeoY = padfGeoTransform[3] + dfPixel * padfGeoTransform[4]
2501 : * + dfLine * padfGeoTransform[5];
2502 : *
2503 : * @param padfGeoTransform Six coefficient GeoTransform to apply.
2504 : * @param dfPixel Input pixel position.
2505 : * @param dfLine Input line position.
2506 : * @param *pdfGeoX output location where geo_x (easting/longitude) location is placed.
2507 : * @param *pdfGeoY output location where geo_y (northing/latitude) location is placed.
2508 : */
2509 :
2510 6 : void CPL_STDCALL GDALApplyGeoTransform( double *padfGeoTransform,
2511 : double dfPixel, double dfLine,
2512 : double *pdfGeoX, double *pdfGeoY )
2513 : {
2514 6 : *pdfGeoX = padfGeoTransform[0] + dfPixel * padfGeoTransform[1]
2515 6 : + dfLine * padfGeoTransform[2];
2516 6 : *pdfGeoY = padfGeoTransform[3] + dfPixel * padfGeoTransform[4]
2517 6 : + dfLine * padfGeoTransform[5];
2518 6 : }
2519 :
2520 : /************************************************************************/
2521 : /* GDALInvGeoTransform() */
2522 : /************************************************************************/
2523 :
2524 : /**
2525 : * Invert Geotransform.
2526 : *
2527 : * This function will invert a standard 3x2 set of GeoTransform coefficients.
2528 : * This converts the equation from being pixel to geo to being geo to pixel.
2529 : *
2530 : * @param gt_in Input geotransform (six doubles - unaltered).
2531 : * @param gt_out Output geotransform (six doubles - updated).
2532 : *
2533 : * @return TRUE on success or FALSE if the equation is uninvertable.
2534 : */
2535 :
2536 1298 : int CPL_STDCALL GDALInvGeoTransform( double *gt_in, double *gt_out )
2537 :
2538 : {
2539 : double det, inv_det;
2540 :
2541 : /* we assume a 3rd row that is [1 0 0] */
2542 :
2543 : /* Compute determinate */
2544 :
2545 1298 : det = gt_in[1] * gt_in[5] - gt_in[2] * gt_in[4];
2546 :
2547 1298 : if( fabs(det) < 0.000000000000001 )
2548 0 : return 0;
2549 :
2550 1298 : inv_det = 1.0 / det;
2551 :
2552 : /* compute adjoint, and devide by determinate */
2553 :
2554 1298 : gt_out[1] = gt_in[5] * inv_det;
2555 1298 : gt_out[4] = -gt_in[4] * inv_det;
2556 :
2557 1298 : gt_out[2] = -gt_in[2] * inv_det;
2558 1298 : gt_out[5] = gt_in[1] * inv_det;
2559 :
2560 1298 : gt_out[0] = ( gt_in[2] * gt_in[3] - gt_in[0] * gt_in[5]) * inv_det;
2561 1298 : gt_out[3] = (-gt_in[1] * gt_in[3] + gt_in[0] * gt_in[4]) * inv_det;
2562 :
2563 1298 : return 1;
2564 : }
2565 :
2566 : /************************************************************************/
2567 : /* GDALSerializeTransformer() */
2568 : /************************************************************************/
2569 :
2570 26 : CPLXMLNode *GDALSerializeTransformer( GDALTransformerFunc pfnFunc,
2571 : void *pTransformArg )
2572 :
2573 : {
2574 26 : VALIDATE_POINTER1( pTransformArg, "GDALSerializeTransformer", NULL );
2575 :
2576 26 : GDALTransformerInfo *psInfo = static_cast<GDALTransformerInfo *>(pTransformArg);
2577 :
2578 26 : if( psInfo == NULL || !EQUAL(psInfo->szSignature,"GTI") )
2579 : {
2580 : CPLError( CE_Failure, CPLE_AppDefined,
2581 0 : "Attempt to serialize non-GTI transformer." );
2582 0 : return NULL;
2583 : }
2584 26 : else if ( psInfo->pfnSerialize == NULL )
2585 : {
2586 : CPLError( CE_Failure, CPLE_AppDefined,
2587 0 : "No serialization function available for this transformer." );
2588 0 : return NULL;
2589 : }
2590 : else
2591 : {
2592 26 : return psInfo->pfnSerialize( pTransformArg );
2593 : }
2594 : }
2595 :
2596 : /************************************************************************/
2597 : /* GDALRegisterTransformDeserializer() */
2598 : /************************************************************************/
2599 :
2600 : static CPLList* psListDeserializer = NULL;
2601 : static void* hDeserializerMutex = NULL;
2602 :
2603 : typedef struct
2604 : {
2605 : char* pszTransformName;
2606 : GDALTransformerFunc pfnTransformerFunc;
2607 : GDALTransformDeserializeFunc pfnDeserializeFunc;
2608 : } TransformDeserializerInfo;
2609 :
2610 1093 : void* GDALRegisterTransformDeserializer(const char* pszTransformName,
2611 : GDALTransformerFunc pfnTransformerFunc,
2612 : GDALTransformDeserializeFunc pfnDeserializeFunc)
2613 : {
2614 : TransformDeserializerInfo* psInfo =
2615 1093 : (TransformDeserializerInfo*)CPLMalloc(sizeof(TransformDeserializerInfo));
2616 1093 : psInfo->pszTransformName = CPLStrdup(pszTransformName);
2617 1093 : psInfo->pfnTransformerFunc = pfnTransformerFunc;
2618 1093 : psInfo->pfnDeserializeFunc = pfnDeserializeFunc;
2619 :
2620 1093 : CPLMutexHolderD(&hDeserializerMutex);
2621 1093 : psListDeserializer = CPLListInsert(psListDeserializer, psInfo, 0);
2622 :
2623 1093 : return psInfo;
2624 : }
2625 :
2626 : /************************************************************************/
2627 : /* GDALUnregisterTransformDeserializer() */
2628 : /************************************************************************/
2629 :
2630 1057 : void GDALUnregisterTransformDeserializer(void* pData)
2631 : {
2632 1057 : CPLMutexHolderD(&hDeserializerMutex);
2633 1057 : CPLList* psList = psListDeserializer;
2634 1057 : CPLList* psLast = NULL;
2635 2114 : while(psList)
2636 : {
2637 1057 : if (psList->pData == pData)
2638 : {
2639 : TransformDeserializerInfo* psInfo =
2640 1057 : (TransformDeserializerInfo*)pData;
2641 1057 : CPLFree(psInfo->pszTransformName);
2642 1057 : CPLFree(pData);
2643 1057 : if (psLast)
2644 0 : psLast->psNext = psList->psNext;
2645 : else
2646 1057 : psListDeserializer = NULL;
2647 1057 : CPLFree(psList);
2648 1057 : break;
2649 : }
2650 0 : psLast = psList;
2651 0 : psList = psList->psNext;
2652 1057 : }
2653 1057 : }
2654 :
2655 :
2656 : /************************************************************************/
2657 : /* GDALDeserializeTransformer() */
2658 : /************************************************************************/
2659 :
2660 126 : CPLErr GDALDeserializeTransformer( CPLXMLNode *psTree,
2661 : GDALTransformerFunc *ppfnFunc,
2662 : void **ppTransformArg )
2663 :
2664 : {
2665 126 : *ppfnFunc = NULL;
2666 126 : *ppTransformArg = NULL;
2667 :
2668 126 : CPLErrorReset();
2669 :
2670 126 : if( psTree == NULL || psTree->eType != CXT_Element )
2671 : CPLError( CE_Failure, CPLE_AppDefined,
2672 0 : "Malformed element in GDALDeserializeTransformer" );
2673 126 : else if( EQUAL(psTree->pszValue,"GenImgProjTransformer") )
2674 : {
2675 84 : *ppfnFunc = GDALGenImgProjTransform;
2676 84 : *ppTransformArg = GDALDeserializeGenImgProjTransformer( psTree );
2677 : }
2678 42 : else if( EQUAL(psTree->pszValue,"ReprojectionTransformer") )
2679 : {
2680 0 : *ppfnFunc = GDALReprojectionTransform;
2681 0 : *ppTransformArg = GDALDeserializeReprojectionTransformer( psTree );
2682 : }
2683 42 : else if( EQUAL(psTree->pszValue,"GCPTransformer") )
2684 : {
2685 0 : *ppfnFunc = GDALGCPTransform;
2686 0 : *ppTransformArg = GDALDeserializeGCPTransformer( psTree );
2687 : }
2688 42 : else if( EQUAL(psTree->pszValue,"TPSTransformer") )
2689 : {
2690 0 : *ppfnFunc = GDALTPSTransform;
2691 0 : *ppTransformArg = GDALDeserializeTPSTransformer( psTree );
2692 : }
2693 42 : else if( EQUAL(psTree->pszValue,"GeoLocTransformer") )
2694 : {
2695 0 : *ppfnFunc = GDALGeoLocTransform;
2696 0 : *ppTransformArg = GDALDeserializeGeoLocTransformer( psTree );
2697 : }
2698 42 : else if( EQUAL(psTree->pszValue,"RPCTransformer") )
2699 : {
2700 0 : *ppfnFunc = GDALRPCTransform;
2701 0 : *ppTransformArg = GDALDeserializeRPCTransformer( psTree );
2702 : }
2703 42 : else if( EQUAL(psTree->pszValue,"ApproxTransformer") )
2704 : {
2705 42 : *ppfnFunc = GDALApproxTransform;
2706 42 : *ppTransformArg = GDALDeserializeApproxTransformer( psTree );
2707 : }
2708 : else
2709 : {
2710 0 : GDALTransformDeserializeFunc pfnDeserializeFunc = NULL;
2711 : {
2712 0 : CPLMutexHolderD(&hDeserializerMutex);
2713 0 : CPLList* psList = psListDeserializer;
2714 0 : while(psList)
2715 : {
2716 : TransformDeserializerInfo* psInfo =
2717 0 : (TransformDeserializerInfo*)psList->pData;
2718 0 : if (strcmp(psInfo->pszTransformName, psTree->pszValue) == 0)
2719 : {
2720 0 : *ppfnFunc = psInfo->pfnTransformerFunc;
2721 0 : pfnDeserializeFunc = psInfo->pfnDeserializeFunc;
2722 0 : break;
2723 : }
2724 0 : psList = psList->psNext;
2725 0 : }
2726 : }
2727 :
2728 0 : if (pfnDeserializeFunc != NULL)
2729 : {
2730 0 : *ppTransformArg = pfnDeserializeFunc( psTree );
2731 : }
2732 : else
2733 : {
2734 : CPLError( CE_Failure, CPLE_AppDefined,
2735 : "Unrecognised element '%s' GDALDeserializeTransformer",
2736 0 : psTree->pszValue );
2737 : }
2738 : }
2739 :
2740 126 : return CPLGetLastErrorType();
2741 : }
2742 :
2743 : /************************************************************************/
2744 : /* GDALDestroyTransformer() */
2745 : /************************************************************************/
2746 :
2747 200 : void GDALDestroyTransformer( void *pTransformArg )
2748 :
2749 : {
2750 200 : GDALTransformerInfo *psInfo = (GDALTransformerInfo *) pTransformArg;
2751 :
2752 200 : if( psInfo == NULL || !EQUAL(psInfo->szSignature,"GTI") )
2753 : {
2754 : CPLError( CE_Failure, CPLE_AppDefined,
2755 0 : "Attempt to destroy non-GTI transformer." );
2756 : }
2757 : else
2758 200 : psInfo->pfnCleanup( pTransformArg );
2759 200 : }
2760 :
2761 : /************************************************************************/
2762 : /* GDALUseTransformer() */
2763 : /************************************************************************/
2764 :
2765 34 : int GDALUseTransformer( void *pTransformArg,
2766 : int bDstToSrc, int nPointCount,
2767 : double *x, double *y, double *z,
2768 : int *panSuccess )
2769 : {
2770 34 : GDALTransformerInfo *psInfo = (GDALTransformerInfo *) pTransformArg;
2771 :
2772 34 : if( psInfo == NULL || !EQUAL(psInfo->szSignature,"GTI") )
2773 : {
2774 : CPLError( CE_Failure, CPLE_AppDefined,
2775 0 : "Attempt to use non-GTI transformer." );
2776 0 : return FALSE;
2777 : }
2778 : else
2779 : return psInfo->pfnTransform( pTransformArg, bDstToSrc, nPointCount,
2780 34 : x, y, z, panSuccess );
2781 : }
2782 :
|
