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