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