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