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